RELAP5铅铋快堆模型拓展及验证

张家心，王成龙,*，赵寒冰，张大林，苏光辉，秋穗正，田文喜

(1.西安交通大学 核科学与技术学院，陕西 西安 710049；2.中国核动力研究设计院，四川 成都 610041)

铅铋快堆是6种第4代反应堆中应用前景广阔的一种，具有固有安全性高、自然循环能力强、可小型化的独特优势，可用于海洋核动力、移动式小型核电源、边远地区用电等方面[1]。同时快堆具备增殖、嬗变的功能，有助于核能可持续发展，改善能源结构[2]。因此，铅铋快堆的设计与安全分析成为目前核能研究的重要方向。热工水力分析程序在评估反应堆安全性时起到关键作用，目前铅铋快堆通用的热工水力分析程序几乎都是对其他反应堆系统分析程序进行修改，如RELAP5、RELAP5/MOD4、ATHLET、TRAC/AAA、SAS4/SASSYS和SIMMER等[3]，这些商用系统分析程序已经过校核与验证，对其进行二次开发只需完成对新添加的流体进行实验验证，保证铅铋流体模型的正确性和程序计算的准确性。

Ma等[4]指出不同程序计算结果的对比可保证程序计算的一致性，但不能取代实验验证，故使用TALL实验台架对TRAC/AAA进行瞬态实验验证。RELAP5系列程序是目前修改版本最多的程序。RELAP5-3D[5]是美国爱德华国家实验室推出的反应堆系统分析程序，该程序在RELAP5 ATHENA程序的基础上增加了氨、二氧化碳、联苯-联苯酚、甘油、铅铋合金(LBE)、锂铅合金、熔盐、钠和钠钾合金等工质流体和磁流体模型，其中液态金属非棒束通道采用Subbotin关系式，棒束通道间换热关系式采用Kazimi-Carelli关系式。RELAP5/MOD3系列则在不同单位均有修改，如比萨大学修改的RELAP5/MOD3.3(modified)[6]包含液态金属铅、铅铋合金、铅锂合金和钠的物性和壁面对流换热关系式，并通过CIRCE-ICE台架进行验证；印度理工学院对RELAP5/MOD4.0[7]进行铅铋流体改造，采用OECD/NEA推荐的铅铋物性，通过NACIE台架进行验证。

国内的通用系统分析程序通常在RELAP5的基础上进行二次开发，如南华大学开发RELAP5_LEAD程序并与FLUENT进行多尺度耦合，完成Code-to-Code的初步验证[8]。中国原子能科学研究院曾对RELAP5/MOD3进行铅铋快堆拓展与改造，但是研究成果尚未公开[3]。

本文对RELAP5开展铅铋流体适用性拓展和实验验证，添加的铅铋流体模型包括物性模型和换热模型，利用回路式NACIE-UP台架和池式CIRCE-ICE台架实验数据对程序进行验证。

1 液态铅铋合金模型

1.1 液相铅铋物性模型

液相铅铋物性模型主要采用OECD/NEA于2007年公布的《铅与铅铋合金物性手册》[9]，包括密度、比定压热容、黏度、热导率、等压膨胀系数、表面张力、声速。

密度ρ为：

ρ=11 096-1.323 6T

(1)

式中，T为流体温度，K。该公式在400～1 100 K温度范围内相对误差不超过0.8%。

比定压热容cp为：

cp=159-2.72×10-2T+7.12×10-6T2

(2)

该公式在400～1 100 K温度范围内相对误差不超过5%。

黏度η为：

(3)

该公式在400～1 100 K温度范围内相对误差不超过5%。

热导率λ为：

λ=3.61+1.517×10-2T-1.741×10-6T2

(4)

等压热膨胀率β为：

(5)

表面张力σ为：

σ=0.437-6.6×10-5T

(6)

声速C为：

C=1 773+0.104 9T-2.873×10-4T2

(7)

饱和压力ps与温度关系式采用Morita关系式[10]：

lnps=35.773+2.800 6×10-4T-

24 053/T-1.640 2lnT

(8)

RELAP5的求解过程需要等温压缩率，在实验数据缺乏且没有参考文献的情况下，这里采用了理论推导的方式计算等温压缩率。由热力学性质可知：

(9)

式中：γ为绝热指数；v为比体积；p为压力；下标s和T分别表示等熵过程和等温过程。

对于铅铋合金流体γ≈1，因此：

(10)

式中：βs为等熵压缩率；βT为等温压缩率。

焓和内能采用定义式计算：

(11)

u=h-pv

(12)

式中：h为比焓；u为比内能。

1.2 气相铅铋物性模型

铅铋气相的物性尚未查到实验关系式，且铅铋物性很高，一般不会达到沸腾的程度，气相也不会参与计算，因此铅铋气相物性采用理想气体模型，这样既可满足程序两流体方程的需要，又可避免改变源代码结构，减少代码修改工作量。

《铅与铅铋合金物性手册》(2007版)[9]仅给出1.01×105Pa下的铅铋合金标准汽化潜热hfg0。RELAP5通过内能求解温度，为保证温度和汽化潜热模型计算准确，采用以下简化和假设：1) 气相比焓根据等压线下的饱和气比焓计算；2) 饱和气比焓在饱和线上采用理想气体模型。饱和气比焓可按式(13)～(15)计算。

(13)

(14)

(15)

将1×105Pa下气相焓值作为参考焓值，可得气相比焓与比内能公式。

(16)

ug=cVT+937 693.243 8

(17)

1.3 液态铅铋换热模型

目前在RELAP5中添加圆管换热关系式和棒束换热关系式。圆管换热关系式选择Xu Cheng关系式和Subbotin关系式[11]，这两个换热关系式均在目前的铅铋快堆系统分析程序中普遍应用。

Xu Cheng关系式为：

Nu=A+0.018Pe

(18)

(19)

Subbotin关系式为：

Nu=5.0+0.025Pe0.8

(20)

Subbotin关系式适用范围为200≤Pe≤1 150。

棒束换热模型采用Borishanski关系式：

(21)

式中：A为与贝克莱数Pe有关的数；NuL为层流努塞尔数；p为栅距；d为单棒外径。

1.4 流动阻力

RELAP5本身无法根据管道类型计算壁面摩擦系数，所以棒束通道和特殊部件的阻力系数通常在内置壁面摩擦系数的基础上增加局部阻力系数。NACIE-UP台架加热段采用绕丝棒束，Forgione等[12]采用CFD方法对Cheng-Todreas模型和RELAP5内置壁面摩擦模型进行对比，给出了NACIE-UP台架加热段燃料棒模拟装置(FPS)增加的局部阻力系数K。

(22)

CIRCE-ICE台架流量计采用文丘里流量计，Narcisi等[13]使用RELAP5-3D时给出了文丘里流量计的局部损失系数KVenturi。

(23)

2 验证实验

2.1 NACIE-UP台架

NACIE-UP台架是欧洲为了开展池式整体性实验，为CIRCE-ICE台架提供技术支持建立的回路式台架。图1a为NACIE-UP结构示意图，其主要部件包括FPS、上升段、膨胀箱、换热器(HX)、下降段和流量计。

图1 NACIE-UP台架结构示意图(a)和节点图(b)Fig.1 Structure scheme (a) and node graph (b) of NACIE-UP facility

1) ADP00实验工况

ADP00是NACIE-UP台架用于研究燃料棒非均匀加热局部影响的实验工况，其实验工况稳态参数列于表1[14]。该实验为无保护失流事故，第1阶段通过注气口向上升管下端注入氩气，在上升管中形成LBE-氩气两相流，依靠相间曳力和上升段与下降段密度差建立稳定的强迫循环。当系统达到稳定状态后，开始实验第2阶段，保持FPS总功率不变，停止注气，实验装置逐渐过渡至自然循环，达到新的稳态。实验过程中保持二次侧给水压力为1.6×106Pa，进口温度为170 ℃，总体积流量为10 m3/h。

表1 ADP00工况稳态参数Table 1 Steady-state parameter of ADP00 case

该实验工况的RELAP5节点模型如图1b所示。130控制体为FPS活性区，170控制体为换热器。换热器结构如图2所示。一、二次侧中间设有304不锈钢粉末夹层，目前对该结构不锈钢粉末热导率的估计值在0.5～2.0 W/(m·K)之间，本文采用1.11～1.49 W/(m·K)。

图2 NACIE-UP换热器示意图[15]Fig.2 Scheme of NACIE-UP HX[15]

2) TEST-3实验工况

TEST-3工况是功率和质量流量过渡实验，从高功率、高质量流量的强迫循环过渡到中低功率、低质量流量的自然循环。初始稳态下，FPS总功率为100 kW，气体流量为20 NL/min；过渡时FPS功率降低至10 kW，功率变化率为10 kW/s，同时关闭注气装置，气体流量约1 s后降至0；在过渡之后，建立新的稳定状态并维持一段时间。实验过程保持二次侧给水压力为1.6×106Pa，进口温度为170 ℃，总容积流量为10 m3/h[15]。TEST-3与ADP00工况采用相同的RELAP5节点模型，但TEST-3工况采用七管的管壳式换热器。

2.2 CIRCE-ICE台架

在CIRCE-ICE台架进行了两组混合对流与热分层实验(TEST Ⅰ和TEST Ⅱ)[16]，程序选择TEST Ⅰ实验数据进行验证。混合对流与热分层实验本质为保护性失流事故，即在获得稳定初始工况参数后，在指定时刻(实验的第7 h时)停止注入氩气，降低功率，强迫循环逐渐过渡至自然循环，同时停止给水，启动余热排出热交换器(DHX)。实验进程参数列于表2。

RELAP5的CIRCE-ICE节点图如图3所示。其中130控制体为上铅池；410控制体为DHX壳侧；310控制体为HX壳侧；180和190控制体为下铅池；250和260控制体分别为进料管和FPS；210控制体为上升管；520和620控制体分别为HX和DHX二次侧。

表2 实验进程Table 2 Experimental schedule

图3 CIRCE-ICE节点图Fig.3 Node graph of CIRCE-ICE

2.3 结果与分析

1) ADP00工况计算结果

ADP00工况主要对FPS进出口的冷却剂温度进行对比，选择FPS之外的测点TP101和TP102作为验证对象。ADP00的流量对比如图4a所示。由图4a可看出，流量曲线整体趋势与实验值相似，稳定自然循环的流量为1.32 kg/s，实验值为1.31 kg/s。

图4b为ADP00工况温度曲线。由图4b可看出，初态与末态的相对误差不超过3%，瞬态过程出口温度计算值与实验值趋势相同，相对误差不超过5%。但瞬态进口温度计算值与实验值相差较大，主要原因为：1) 附加结构热惯性的影响；2) 数据采集控制系统的调节过程未知。因为NACIE-UP管道直径仅有6.35 cm，台架附加的各种固定支撑结构和保温、绝热等材料能明显影响其热惯性，但这些参数无法具体估计，只能通过合理假设获得相近结果。另外控制系统在瞬态过程中对各阀门开度、阻力件、功率的调控过程未知，因此瞬态过程误差要比稳态情况的大，但是最大相对误差不超过5%。瞬态初期FPS进出口温度模拟较好，其中进口温度在瞬态初期有上升回落过程，这是因为下水平段的流量计采用热流量计，热流量计的加热功率使进口温度呈现上升回落趋势，这一过程的计算结果与实验值变化趋势相同。

图4 ADP00工况流量(a)和温度(b)曲线Fig.4 Flow rate (a) and temperature (b) curves of ADP00 case

2) TEST-3工况计算结果

TEST-3工况认为测点TP101和TP102的温度代表FPS进出口温度。图5a示出TEST-3工况的流量曲线。该实验只有初始稳态，自然循环过渡过程未完全建立稳态。由图5a可看出，计算流量与实验流量变化趋势相同，相对误差不超过2%。在瞬态初期，计算结果稍比实验结果波动剧烈，原因可能是RELAP5节点模型的阻力与实际阻力不完全一致。

图5b为TEST-3工况温度曲线。由图5b可看出，该工况未达到自然循环稳态，但初始稳态的温度计算结果与实验相比，相对误差不到1%。瞬态过程误差比稳态略大，但整体趋势保持一致，最大相对误差不超过5%。其中瞬态初期进口温度依然出现上升回落现象，这与APD00工况的原因相同，因为热流量计具有一定热功率。

图5 TEST-3工况流量(a)和温度(b)曲线Fig.5 Flow rate (a) and temperature (b) curves of TEST-3 case

图6 TEST-3工况流量计算结果对比Fig.6 Calculation result comparison of flow rate of TEST-3 case

图6示出TEST-3工况流量计算结果对比，除本文开发的系统分析程序外，其余4个程序分别为ENEA的CATHARE、GRS的ATHLET、罗马大学的RELAP5-3D和比萨大学改进的RELAP5/MOD3.3(modified)[16]。从图6可看出，本文程序与CATHARE和ATHLET程序计算结果最接近，其他程序流量计算误差比本程序误差大。

图7示出TEST-3工况进出口温度计算结果对比。由图7a可见，进口温度变化趋势几乎相同，最大相对偏差不超过10%。由图7b可见，出口温度变化趋势相似，但ATHLET与RELAP5/MOD3.3(modified)在瞬态开始瞬间出现了明显尖峰，一般来说可能是由于节点模型划分和边界条件简化过程造成。从稳态和瞬态结果来看，本文程序计算结果要优于其他4个程序。

3) TEST Ⅰ工况计算结果

图8示出TEST Ⅰ工况的流量和温度曲线。由图8a可看出，初始稳态和瞬态后期流量符合较好，但是在瞬态初期有一定误差，这是因为TEST Ⅰ工况没有公开瞬态控制过程、DHX和HX的边界条件，因此RELAP5节点模型和边界条件设定可能与实验存在一定偏差。整体上看，当功率按照实验所给的参数，程序可计算得到相似的变化趋势。

由图8b可见：初始稳态的温度与实验值符合较好，相对误差不超过2%；瞬态过程与实验值误差较大，但相对误差在10%以内。与回路式台架相比，瞬态误差偏大的原因为：1) 实验条件和边界条件不完整可能造成误差，该实验并未给出HX给水温度和DHX空气进口温度，给水温度作为热阱对一次侧温度有着明显影响，本文根据其设计参数进行计算；2) 池式台架在铅池内可能存在明显的搅混等三维热工水力现象，搅混等三维热工水力现象可能会对池内换热过程有影响，导致一维系统程序对温度分布计算不准确，但趋势相同；3) 铅池部分包括铅池、HX壳侧、DHX壳侧3个并联通道，3个通道阻力特性对各部分流量分配有着显著影响，但各通道阻力特性未知，本程序未对各通道添加额外的阻力系数，因此HX和DHX换热量可能与实验不完全相同。DHX换热能力不足以排除余热，进而LBE温度持续升高。

图7 TEST-3工况进出口温度计算结果对比Fig.7 Calculation result comparison of inlet and outlet temperatures of TEST-3 case

图8 TEST Ⅰ工况流量(a)和温度(b)曲线Fig.8 Flow rate (a) and temperature (b) curves of TEST Ⅰ case

3 结论

本文对商业系统分析程序RELAP5进行液态铅铋合金流体适用性拓展，并对修改后的程序进行实验验证，主要结论如下。

1) 本文修改后的RELAP5程序可实现对液态铅铋合金物性和流动换热过程计算，可实现回路式系统和池式系统的模拟。

2) 对于NACIE-UP台架，本文程序计算的温度的稳态相对误差在2%以内，瞬态相对误差在5%以内，流量计算结果与实验结果趋势一致，精度符合系统分析程序的需求；与其他系统分析程序相比，计算结果趋势相同，相对偏差在10%以内，且受各程序节点模型和结构参数准确性的影响。

3) 对于CIRCE-ICE台架，本文程序计算的温度的稳态相对误差在2%以内，瞬态相对误差在10%以内。