基于Newton-Raphson算法的COSINE热工水力程序安注箱模块研发

陈 俊，沙会娥，吴照国 ，孔晓宁，尹铁男，林 萌,杨燕华,

(1.国家电投集团科学技术研究院有限公司 国家能源核电软件重点实验室， 北京 102209；2.上海交通大学 核能科学与工程学院，上海 200240)

核电设计软件用于核电厂设计、安全分析及运行管理，是支撑我国先进核电技术持续创新发展的关键重大共性技术。COSINE软件包是国家核电技术有限公司北京软件技术中心自主研发的堆芯设计和安全分析一体化软件包，其中热工水力安全分析程序包括系统程序、子通道程序及安全壳程序。

非能动安全注入系统在非LOCA事故情况下，可对反应堆冷却剂系统(RCS)进行补水和硼化，在LOCA事故下可对RCS进行安全注入[1]。安注箱属中压安注，是核电厂重要的专设安全设施，目前国内外针对他的实验分析和模拟分析尚不成熟[2-4]。

安注箱模型主要与系统程序联立计算。与其他程序对安注箱模型特殊设置不同[2，3,5-7]，COSINE系统程序安注箱模块有独立的守恒方程和求解算法，既可以与COSINE系统程序联立求解，也可以安注箱模块独立程序求解。

本文首先介绍COSINE系统程序安注箱模块的物理模型和数学算法；其次分析了COSINE安注箱模型与参考程序RELAP5结果的对比计算。

1 安注箱物理模型及算法设计

现有程序假设认为安注箱气体空间(如氮气)内没有水蒸气；实际上气空间与液空间处于动态平衡，气相空间里存在部分通过相变转换的蒸汽组分。本文假设气体空间中存在蒸汽，并以饱和状态均匀分散在气相空间中；如果蒸汽含量足够小，可以通过设置或者计算体积份额来判断或者设置是否考虑蒸汽相。水蒸气和不凝气体共同占据气相的体积，并具备相同的体积份额。由于液体比热容大，密度变化可以忽略不计，整个安注过程中状态基本不变；本文假设安注过程中液相温度维持不变。

1.1 守恒方程

1.1.1 质量守恒方程

质量守恒方程主要考虑不凝气体连续方程和液相连续方程。在一般电厂中，不凝气体都采用氮气，此处设计不凝气体为后续扩展用；用户可以选择氮气，也可以选其他不凝气体作为安注气体。相应的质量守恒方程表述如下。

不凝气体的质量守恒方程：

(1)

液相质量守恒方程：

(2)

式中：αn——不凝气体体积份额；

ρn——不凝气体密度；

ρf——液相密度；

V——安注箱体积；

AL——安注箱流通面积；

Γf——液相质量变化率；

u——安注速度。

1.1.2 能量守恒方程

由于液相在安注过程中状态基本保持不变，认为是恒温过程，只需建立气相能量守恒方程：

(3)

式中：hn——不凝气体焓；

hv——蒸汽焓；

pTK——安注箱压力；

H——气液相换热系数；

QW——壁面传热量。

1.1.3 动量守恒方程

安注过程可能安注气相，也可能安注液相，或者气液混合相。假设安注过程中气液两相具备相同流速，由于蒸汽和不凝气体密度接近，为简化计算，小量蒸汽部分用不凝气体替代，相应的动量守恒方程可以表示为：

(4)

式中：Lf,TK——安注箱内液相高度；

Lf,L——安注管内液柱长度；

Lg,L——安注管内气体长度；

λ——摩擦因子；

KL——局部损失系数。

1.2 数值算法

守恒方程采用半隐式离散方法。首先根据动量守恒方程得到速度与压力之间的关系式，并利用此关系式得到临时速度，用于质量守恒方程和能量守恒方程的更新计算。速度的压力表达式可以表述为：

(5)

式中：C——速度系数；

A——方程显式项；

B——压力系数。

其次联立质量守恒方程和能量守恒方程，建立质量守恒方程和能量守恒方程的线性方程组；利用Newton-Raphson算法计算下一时刻气相体积份额，不凝气体焓，以及不凝气体压力。该线性方程组包括不凝气体质量残差方程，液相质量残差方程以及气相能量残差方程。质能方程求解时选取气体体积份额，不凝气体焓以及不凝气体压力为主变量。相应的控制体线性方程组可以表示为：

(6)

式中：εc,n——不凝气体质量残差表达式；

εc,f——液相质量残差表达式；

εe,g——气相能量残差表达式；

δαn——气相体积份额变化率；

δhn——不凝气体焓变化率；

δpn——不凝气体压力变化率。

2 结果与讨论

根据安注箱物理模型和数学算法，将COSINE计算结果与RELAP5相对比。为排除系统等其他因素对安注箱计算程序影响，程序系统一回路处理为边界参数输入安注箱。安注箱主要有球形和柱形，本文选取柱形安注箱为对象。

选取定背压工况和变背压工况作为安注箱边界条件进行结果验证，其柱形安注箱几何形状如图1所示，几何参数如表1所示。

表1 柱形安注箱几何参数Table 1 Geometry of cylindrical accumulator

图1 柱形安注箱示意图Fig.1 Tank of accumulator

2.1 定背压工况

对定背压边界工况，设置安注箱初始压力为6 MPa，边界出口背压分别设置为1 MPa和3 MPa两种情况。安注箱的表面粗糙度设置为0.000 6，局部阻力损失系数设置为20.0。

图2、图3给出不同定背压条件下安注箱流量随时间变化趋势图。虚线为COSINE计算结果，实线为RELAP5计算结果(下同)。当背压为1 MPa时(见图2)，RELAP5计算的安注箱安注入量在110 s是趋近于零，而COSINE在120 s之后；当背压为3 MPa时(见图3)，RELAP5安注箱注入量在25 s左右趋于0，而COSINE在30 s左右趋于0。另对比图2、图3中的刚开始几秒，RELAP会发生剧烈震荡，COSINE计算结果稳定性要稍好于RELAP5。

图2 安注箱流量变化趋势图(背压为1 MPa)Fig.2 Flow rate of accumulator(Pback=1 MPa)

图3 安注箱流量变化趋势图(背压为3 MPa)Fig.3 Flow Rate of accumulator(Pback=3 MPa)

图4、图5给出不同定背压下安注箱内压力随时间变化趋势图。当背压为1 MPa时(见图4)，RELAP5和COSINE安注箱压力在120 s左右趋于零；当背压为3 MPa时(见图5)，RELAP5和COSINE安注箱压力在30 s左右趋于0。

图4 安注箱压力变化趋势图(背压为1 MPa)Fig.4 Tank pressure of accumulator(Pback=1 MPa)

图5 安注箱压力变化趋势图(背压为3 MPa)Fig.5 Tank pressure of accumulator(Pback=3 MPa)

造成COSINE与RELAP5计算结果之间差异的主要可能原因为RELAP5采用了三个守恒方程：质量守恒方程，动量守恒方程和能量守恒方程；而COSINE的安注箱模型采用四个守恒方程，两个质量守恒方程，动量守恒方程及能量守恒方程。具体分析如下：

(1)质量守恒方程

RELAP5质量守恒方程认为空气中气相不考虑任何相变，其质量为常数，不考虑流动，忽略一切其他因素影响：

Mc,n=const=ρc,nV

(7)

本文安注箱程序采用如公式(1)、公式(2)所示的两个质量守恒方程。首先，COSINE程序安注箱模块考虑液体的相变，当相变量很小时两者差别并不大，一旦发生相变则RELAP5会产生较大误差。

其次，COSINE考虑质量守恒方程的对流作用；由于RELAP5没有考虑对流项的影响，如果系统初始给的速度不合理，程序计算容易发散；COSINE则不会遇到此问题。因此考虑对流项可使程序不同状态之间变换更具平滑性和稳定性。

最后，COSINE还给出了安注箱液相质量守恒方程，是程序计算更具适用性。

(2)能量守恒方程

RELAP5的能量守恒方程没有考虑蒸汽部分的能量变化，而COSINE将不凝气体(如氮气等)和蒸汽同时考虑成气相，并合并成统一能量守恒方程。此改动与质量守恒方程考虑蒸汽的原因一致。

(3)求解算法

RELAP5将质量守恒方程、能量守恒方程以及动量守恒方程，结合相应的状态方程将三个方程化简，结合相应的假设折合成一等价公式，并利用此公式进行计算。COSINE利用选取的四个守恒方程，在基于Newton-Raphson算法的基础上结合四个守恒方程设计具备自主知识产权的求解方法；此方法利用现有的守恒方程，不进行任何假设，采用纯数学的方法设计，具有通用性和扩展性。

2.2 变背压工况

图6给出背压随时间的变化趋势，其背压从初始4.38 MPa变化到150 s附近的0.9 MPa。

图6 安注箱背压变化趋势图Fig.6 Back pressure of accumulator

变背压工况下安注箱注入量随时间的变化趋势如图7所示。从结果上看两者变化趋势一致，RELAP5计算的波动性要大于COSINE计算结果。具体原因如前面所述，守恒方程选取、求解方法不同以及假设等因素所致。

图7 安注箱注入量变化趋势Fig.7 Flow rate of accumulator(variable back pressure)

图8给出了安注箱压力随时间变化趋势。COSINE与RELAP5计算趋势一致，但过程中存在差异，相对误差并不大。

图8 安注箱压力变化趋势图Fig.8 Pressure of accumulator(variable back pressure)

3 结论

本文建立安注箱的物理模型和数学算法，并在此基础上编写成程序。通过对该模型的分析，以及不同条件下与RELAP5的计算结果对比分析，其结论如下：

(1)该安注箱程序既可作为安注箱设备程序独立计算，也可以嵌入COSINE程序作为安全箱模块与COSINE程序整体耦合计算；

(2)通过设计不同工况下算例，将此程序与现有的系统程序RELAP5计算结果进行对比，整体趋势吻合较好，在局部情况下较RELAP5更具稳定性和收敛性。造成此差异的可能主要原因有以下两点：

1)两者守恒方程的选取和模型简化；

2)两者算法的设计原理不同。

本文的计算算法是COSINE热工水力程序的外扩，为具备自主知识产权的COSINE热工水力程序的开发和应用奠定了一定的基础。

致谢

感谢项目国家科技重大专项“核动力厂安全分析用计算机软件评估基准题及共享平台开发”资助(2019ZX06005001)。