随着我国“海洋强国”战略的提出，海洋开发对能源的需求日益迫切。浮动核电站作为一种安全、经济的核电新形式，能够有效地解决海上能源短缺的问题，开展相关的研发工作势在必行。利用系统分析程序评估海洋条件对反应堆安全的影响并制定相应的设计规范是我国浮动核电站研发和安审工作中亟待解决的问题。目前主流的反应堆系统分析程序如RELAP5、TRACE、RETRAN等均是基于陆基核电厂开发而来^[1]，无法针对浮动核电站开展系统安全分析，因此有必要开发适用于海洋条件的反应堆热工水力系统分析程序。

目前海洋条件下反应堆系统分析程序的开发主要可分为两类，即小型自编程序的开发和商用程序改进。由于小型自编程序^[2-5]往往模型算法简单且适用对象单一，难以满足浮动核电站的研发和安审工作的需要。对于商用程序的改进，日本Ishida等^[6]和韩国Jae-Hak Kim等^[7]分别基于RETRAN程序开发了海洋条件瞬态分析程序RETRAN-02/GRAV和RETRAN-03/INT，并采用实验验证了程序计算的准确性；我国鄢炳火等针对RELAP5程序进行了海洋条件的二次开发，但所开发程序针对摇摆运动，没有研究复杂海洋条件对反应堆安全的影响^[8]。谭长禄等基于RELAP5开发了海洋条件系统分析程序RELAP5/MC，但程序未开展实验验证，对复合海洋条件下反应堆系统安全分析的适用性不明^[9-10]。

针对当前相关研究的不足，本文建立了描述复合海洋条件下惯性加速度的数学模型，通过在RELAP5中添加海洋条件计算模块，开发了适用于复杂海洋条件下反应堆热工水力分析的系统分析程序。采用实验验证和程序间对比验证的方法，对程序在运动条件下的计算准确性进行验证。

1 海洋条件惯性加速度模型建立

海洋条件是指核动力船舶受海面风浪、船体机动等影响而产生的不同运动形式。通常采用倾斜、起伏、摇摆这三种基本运动形式来进行海洋条件的建模分析。为方便描述海洋条件下的流体运动，引入两个右手坐标系：固定于地球的世界坐标系O₀X₀Y₀Z₀和固定于运动控制体上的局部坐标系O₁X₁Y₁Z₁。

由于海洋条件下的局部坐标系O₁X₁Y₁Z₁是一种非惯性系，牛顿运动定律不成立，需要建立非惯性系下的动量守恒方程。式(1) 为非惯性系下积分形式的动量守恒方程表达式^[11]：

$\begin{align} & \iiint\limits_{{{\tau }'}}{\frac{{D}'{v}'}{Dt}}\rho \rm{d}{\tau }'=\iiint\limits_{{{\tau }'}}{\left[ \mathit{\boldsymbol{f}}-{{\mathit{\boldsymbol{a}}}_{0}}-\mathit{\boldsymbol{ }}\!\!\omega\!\!\rm{ }\times \left( \mathit{\boldsymbol{ }}\!\!\omega\!\!\rm{ }\times \mathit{\boldsymbol{{r}'}} \right)-\mathit{\boldsymbol{ }}\!\!\beta\!\!\rm{ }\times \mathit{\boldsymbol{{r}'}}- \right.} \\ & \quad \quad \quad \quad \quad \quad \left. 2\mathit{\boldsymbol{ }}\!\!\omega\!\!\rm{ }\times \mathit{\boldsymbol{{v}'}} \right]\rho \rm{d}{\tau }'+\iint\limits_{{{A}'}}{{{\mathit{\boldsymbol{p}}}_{n}}}\rm{d}{A}' \\ \end{align}$

(1)

同惯性系下的动量守恒方程相比，式(1) 增加了中括号中的后四项，分别为平移加速度、向心加速度、切向加速度和科氏加速度。因此，为应用海洋条件下流体的动量守恒方程，需要对船体运动所引入的附加惯性加速度进行数学建模。

1.1 倾斜运动

倾斜运动会引起反应堆几何姿态的改变，这不仅会造成容器中水位的变化，还会引起反应堆冷热源间有效高度差的改变，进而影响自然循环的热驱动力。因此，倾斜运动数学描述的重点为倾斜后系统水力学控制体竖直提升高度的计算。

空间位置坐标为(x₀, y₀, z₀)的控制体倾斜后的新坐标(x₁, y₁, z₁)可以用三维旋转矩阵表示：

$\begin{array}{l} \left[ {\begin{array}{*{20}{l}} {{x_1}}\\ {{y_1}}\\ {{z_1}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {\cos r} & { - \sin r} & 0\\ {\sin r} & {\cos r} & 0\\ {\;0} & {\;0} & 1 \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {\cos \beta } & 0 & {\sin \beta }\\ {\;0} & 1 & {\;0}\\ { - \sin \beta } & 0 & {\cos \beta } \end{array}} \right].\\ \quad \quad \quad \left[ {\begin{array}{*{20}{l}} 1 & {\;0} & {\;0}\\ 0 & {\cos \alpha } & { - \sin \alpha }\\ 0 & {\sin \alpha } & {\cos \alpha } \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{x_0}}\\ {{y_0}}\\ {{z_0}} \end{array}} \right] \end{array}$

(2)

式中：α、β、γ分别为控制体绕世界坐标系X₀、Y₀、Z₀轴倾斜的角度。计算得到的控制体新坐标中的z₁即为倾斜后控制体在世界坐标系Z₀轴方向上的竖直提升高度h_z。

1.2 起伏运动

起伏是指船舶周期性的上下平动，它不会改变系统的几何姿态。通常情况下，起伏运动可被简化描述为一种简谐变化规律。在世界坐标系中，起伏条件下流体的惯性加速度a_heaving可以表示为

${\mathit{\boldsymbol{a}}_{{\rm{heaving}}}} = - \mathit{\boldsymbol{a}}\left( t \right) = \frac{{4{{\rm{ \mathsf{ π} }}^2}{x_m}}}{{{T^2}}}\sin \left( {\frac{{2{\rm{ \mathsf{ π} }}t}}{T}} \right)\mathit{\boldsymbol{k}}$

(3)

式中：x_m为起伏位移的最大值，T为起伏周期，k为世界坐标系中Z₀轴正向的单位向量。

起伏运动的本质为一种竖直方向上的变加速直线运动，对于一般形式的直线加速运动也可以采用同样的方式表示：

${\mathit{\boldsymbol{a}}_{{\rm{line}}}} = - {a_x}\left( t \right)\mathit{\boldsymbol{i}} - {a_y}\left( t \right)\mathit{\boldsymbol{j}} - {a_z}\left( t \right)\mathit{\boldsymbol{k}}$

(4)

式中：a_x、a_y、a_z分别为船体沿世界坐标系X₀、Y₀、Z₀轴方向的直线加速度。

1.3 摇摆运动

摇摆是指船舶绕某一轴线进行的往复运动，它不仅会引起系统几何姿态的变化，还会引入附加的切向加速度和向心加速度。在核动力船舶实际运行时，其摇摆运动可能发生在任意方向上，为模拟船舶受海洋条件影响的实际运动情况，需建立绕任意轴摇摆的惯性加速度模型。

如图 1所示，在向量V=(x₀, y₀, z₀)绕任意单位向量n=(n_x, n_y, n_z)旋转θ角度过程中，任意时刻新向量V′的空间坐标(x_t, y_t, z_t)可以由式(5) 计算的得到^[12]：

$\left[ {\begin{array}{*{20}{l}} {{x_t}}\\ {{y_t}}\\ {{z_t}} \end{array}} \right] = \\ \left[ {\begin{array}{*{20}{l}} {n_x^2\left( {1 - \cos {\theta _t}} \right) + \cos {\theta _t}} & {{n_x}{n_y}\left( {1 - \cos {\theta _t}} \right) - {n_z}\sin {\theta _t}} & {{n_x}{n_z}\left( {1 - \cos {\theta _t}} \right) - {n_z}\sin {\theta _t}}\\ {{n_x}{n_y}\left( {1 - \cos {\theta _t}} \right) - {n_z}\sin {\theta _t}} & {n_y^2\left( {1 - \cos {\theta _t}} \right) + \cos {\theta _t}} & {{n_y}{n_z}\left( {1 - \cos {\theta _t}} \right) - {n_x}\sin {\theta _t}}\\ {{n_x}{n_z}\left( {1 - \cos {\theta _t}} \right) - {n_y}\sin {\theta _t}} & {{n_y}{n_z}\left( {1 - \cos {\theta _t}} \right) - {n_x}\sin {\theta _t}} & {n_z^2\left( {1 - \cos {\theta _t}} \right) + \cos {\theta _t}} \end{array}} \right]\\\left[ {\begin{array}{*{20}{l}} {{x_0}}\\ {{y_0}}\\ {{z_0}} \end{array}} \right]$

(5)

将向量V=(x₀, y₀, z₀)分解成与单位向量n平行的向量V_ll=(V·n)n和垂直的向量V_⊥=V－V_ll，则向量V距旋转中心的距离R=V_⊥。

同理，旋转后的向量V′=(x_t, y_t, z_t)也可以分解成平行于n的V_ll和垂直于n的V′_⊥，则该时刻与V′_⊥同向的单位向量n_c可以计算得到；又因为切向加速度垂直于向量n和向量n_c，则与切向加速度同向的单位向量n_t可以表示为

${\mathit{\boldsymbol{n}}_t} = {\mathit{\boldsymbol{n}}_c} \times \mathit{\boldsymbol{n}}$

(6)

向心加速度矢量a_c和切向加速度矢量a_t可以表示为

$\begin{array}{l} {\mathit{\boldsymbol{a}}_c} = - {w^2}R{\mathit{\boldsymbol{n}}_c} = \left( {{a_{cx}},{a_{cy}},{a_{cz}}} \right)\\ {\mathit{\boldsymbol{a}}_c} = \beta R{\mathit{\boldsymbol{n}}_t} = \left( {{a_{tx}},{a_{ty}},{a_{tz}}} \right) \end{array}$

(7)

所以绕任意轴摇摆控制体的惯性加速度就可以表示为

${\mathit{\boldsymbol{a}}_{{\rm{rolling}}}} = - \left( {{a_{tx}} + {a_{cz}}} \right)\mathit{\boldsymbol{i}} - \left( {{a_{ty}} + {a_{cy}}} \right)\mathit{\boldsymbol{j}} - \left( {{a_{tz}} + {a_{cz}}} \right)\mathit{\boldsymbol{k}}$

(8)

1.4 复合海洋运动

基于以上典型运动形式的海洋条件描述，可以组合得到多运动形式耦合的复合海洋条件惯性加速度模型：

${\mathit{\boldsymbol{a}}_{{\rm{couple}}}} = {\mathit{\boldsymbol{a}}_{{\rm{line}}}} + {\mathit{\boldsymbol{a}}_{{\rm{heaving}}}} + {\mathit{\boldsymbol{a}}_{{\rm{rolling}}}}$

(9)

复合海洋条件所引起的系统控制体空间位置的改变，仍可采用式(5) 计算得到。

2 海洋条件反应堆系统分析程序开发

RELAP5程序是美国Idaho国家实验室开发的一款用于轻水堆设计安全评审的最佳估算瞬态分析程序，采用两流体、非均匀、非平衡的两相流模型，能够模拟包括破口事故、全场断电等几乎所有的核电厂运行瞬态及安全事故，在世界范围内得到了广泛的认可和应用^[13]。本文基于建立的海洋条件惯性加速度模型，通过修改RELAP5程序的动量方程并添加海洋条件计算模块，实现对RELAP5程序的二次开发。

RELAP5中采用的扩展形式的气相和液相动量方

程^[13]分别为

$\begin{array}{l} {\alpha _g}{\rho _g}A\frac{{\partial v_g^2}}{{\partial t}} + \frac{1}{2}{\alpha _g}{\rho _g}A\frac{{\partial v_g^2}}{{\partial x}} = - {\alpha _g}A\frac{{\partial P}}{{\partial x}} + {\alpha _g}{\rho _g}{B_x}A - \\ \quad \quad \quad \left( {{\alpha _g}{\rho _g}A} \right){\rm{FWG}} \cdot {v_g} + {\mathit{\Gamma }_g}A\left( {{v_{gI}} - {v_g}} \right) - \\ \quad \quad \quad \quad \quad \quad \left( {{\alpha _g}{\rho _g}A} \right){\rm{FIG}}\left( {{v_g} - {v_f}} \right) - \\ \quad \quad C{\alpha _g}{\alpha _f}{\rho _m}A\left[ {\frac{{\partial \left( {{v_g} - {v_f}} \right)}}{{\partial t}} + {v_f}\frac{{\partial {v_g}}}{{\partial x}} - {v_g}\frac{{\partial {v_f}}}{{\partial x}}} \right]\\ {\alpha _f}{\rho _f}A\frac{{\partial v_f^2}}{{\partial t}} + \frac{1}{2}{\alpha _f}{\rho _f}A\frac{{\partial v_f^2}}{{\partial x}} = - {\alpha _f}A\frac{{\partial P}}{{\partial x}} + {\alpha _f}{\rho _f}{B_x}A - \\ \quad \quad \quad \left( {{\alpha _f}{\rho _f}A} \right){\rm{FWF}}{v_f} + {\mathit{\Gamma }_g}A\left( {{v_{fI}} - {v_f}} \right) - \\ \quad \quad \quad \quad \quad \left( {{\alpha _f}{\rho _f}A} \right){\rm{FIF}}\left( {{v_f} - {v_g}} \right) - \\ C{\alpha _g}{\alpha _f}{\rho _m}A\left[ {\frac{{\partial \left( {{v_f} - {v_g}} \right)}}{{\partial t}} + {v_g}\frac{{\partial {v_f}}}{{\partial x}} - {v_f}\frac{{\partial {v_g}}}{{\partial x}}} \right] \end{array}$

(11)

式中：等式右边各项依次为：压力梯度项、体积力项、壁面摩擦项、界面质量传输导致的动量传递项、界面摩擦阻力项、虚拟质量力项。

在RELAP5中，重力等体积力是通过B_x引入动量方程的。考虑到海洋条件带来的附加惯性力与重力具有相似的性质^{[6, 9]}，可以将附加惯性力沿世界坐标系分解后，参照与重力相同的处理方式引入RELAP5的动量方程，则处理后的等效加速度a可以表示为

$\mathit{\boldsymbol{a}}{\rm{ = }} - {a_x}\mathit{\boldsymbol{i}} - {a_y}\mathit{\boldsymbol{j}} + \left( {{g_0} - {a_z}} \right)\mathit{\boldsymbol{k}}$

(12)

式中：g₀为重力加速度，a_x、a_y、a_z分别为海洋条件附加惯性加速度在世界坐标系X₀、Y₀、Z₀轴正向上的分量。

图 2为修改后的RELAP5瞬态计算流程示意图。海洋条件参数及相关数据的输入在时间控制模块之后执行，海洋条件附加惯性力的引入及其相关计算在RELAP5瞬态计算中的水力学计算模块中执行。

由于目前针对海洋条件下流体的流动传热、阻力特性、CHF等的研究还处于机理性研究阶段，相关的理论模型尚不完善。因此本研究未进行海洋条件下本构关系式的修改，依然采用RELAP5原有的本构关系式。

3 程序验证

本文采用哈尔滨工程大学开展的摇摆条件下单相强迫循环实验^[14]对引入海洋条件的RELAP5程序在运动条件下强迫循环工况计算的准确性进行验证；同时利用同类型程序RETRAN-02/GRAV^[6]和RELAP5/MC^[9]，通过建立相似的自然循环回路模型，对比各程序的计算结果来验证本文程序在运动条件下自然循环工况计算的准确性。

3.1 摇摆条件下单相强迫循环实验验证

该实验主要研究了摇摆运动下竖直窄矩形通道内单相水低速流动时流量的波动特征。图 3为实验回路的系统布置简图。主实验回路被固定在摇摆实验台上，通过橡胶软管与主实验回路相连的水箱、离心水泵等设备不参与摇摆运动，具体实验过程及系统相关参数详见文献[14]。图 4为该实验回路的RELAP5节点划分示意图。控制体尺寸按实验装置实际尺寸进行设置，系统运行状态为常温常压，Y轴为摇摆轴。为尽可能模拟实验真实情况，水泵被布置在摇摆上，并通过合理设置稳压器的空间坐标来消除摇摆的影响。

计算过程中，通过调节RELAP5程序中的阀门开度及各部件的阻力参数，使系统流量和压降等参数的计算值与静态实验值相匹配。设置与实验工况相同的摇摆运动参数(周期、振幅)进行计算，并将计算结果与实验结果进行对比分析。本文选用表 1所示三种不同的实验工况进行分析验证。

图 5为三组实验工况的RELAP5计算结果与实验数据对比。计算得到的回路流量波动趋势与实验数据符合较好，流量波动周期及相位与实际情况基本相同，但在流量波动极值处存在偏差。综合考虑RELAP5对实验台架进行简化建模过程中引入的误差，可以认为本程序对于摇摆条件下强迫循环工况计算的适用性较好。

3.2 程序对比验证

本节采用日本Ishida等开发的RETRAN-02/GRAV^[6]与国内谭长禄等开发的RELAP5/MC程序^[9-10]进行对比验证。受公开文献中相关计算模型参数完整性的限制，只能通过对比相似系统的热工水力参数受海洋条件影响的变化趋势来验证本文程序。综合考虑文献[6, 10]中所采用的自然循环回路形式，建立如图 6所示的并联双回路单相自然循环计算模型，具体参数详见表 2。

3.2.1 倾斜条件

图 7为RETRAN-02/GRAV、RELAP5/MC与本文程序对不同横斜角度(绕Y轴)下并联双回路自然循环系统的计算结果对比，其中图 7(b)中的实心点代表程序计算结果，空心点为Iyori等实验结果^[15]。可以看出，三种程序得到的系统流量和温度随横摇角度变化趋势基本相同：随着倾斜角度的增加，支路的自然循环流量随着其相对位置的上升或下降而增加或减小，流经加热器的总流量逐渐减少。这主要是由倾斜造成的回路冷热源高度差的改变而引起的自然循环热驱动压头的变化所导致的。由于加热功率恒定，各回路的出口温度与其流量的变化保持对应关系。需要说明的是，对比图 7(e)与图 7(f)，当倾斜角度大于-50°时，RELAP5/MC计算得到的环路1冷段出口温度出现了跃升。这主要是因为文献[9]中采用的自然循环回路尺寸及系统参数与本文模型并不完全相同，当倾斜角度接近-60°时，文献[9]中相对位置下降的环路1出现了零流量，导致系统换热能力下降而造成温度上升。

考虑到各程序采用的计算回路尺寸等参数不同所造成的偏差，可以认为本文程序对倾斜条件下自然循环工况计算具有较好的适用性。

3.2.2 起伏条件

图 8为本文程序计算得到的双回路自然循环系统在起伏条件下的系统流量波动情况。可以看出，在竖直方向变加速度的影响下，自然循环流量围绕其静态值呈正/余弦规律波动。同竖直方向的加速度变化相比，流量的波动存在约1.7 s的延迟，这主要是由于流道阻力造成的。同两个支路相比，加热段的流量波动幅度更大。相同的变化趋势在Ishida等的计算中也被得到^[6]。

起伏运动参数的变化会对自然循环系统的热工水力特性产生影响。图 9(a)和图 9(b)分别为本文程序和RELAP5/MC计算得到的起伏幅值对自然循环流量的影响情况。两种程序得到的影响规律基本相同：随着起伏运动幅值的增加，自然循环系统各回路的流量波动幅度逐步变大。

3.2.3 摇摆条件

图 10和图 11分别为本文程序和RELAP5/MC计算得到的横摇条件下(绕Y轴)双回路自然循环系统的流量波动情况。两种程序得到的流量变化趋势基本相同，即摇摆条件下加热段的流量几何保持恒定，两个并联环路的流量波动相位相反，且冷热段流量的变化相对于摇摆角度存在延迟。

摇摆运动参数的变化会对自然循环系统的热工水力特性产生影响。图 12和图 13分别为本文程序计算得到的横摇周期和幅值对双回路自然循环系统流量波动的影响情况。

从图中可以看出，支路流量的波动幅度随周期的增加而减小，并在周期约为40 s时振荡幅度达到最小值。随着摇摆幅值的增加，支路流量的平均值减小，但流量波动幅度增大。相似的变化趋势在谭长禄等的计算^[10]中也被得到。通过程序间对比验证的方式，对本程序在运动条件下自然循环工况计算的适用性进行了评价。结果表明，相同工况下程序间计算结果的变化趋势符合良好，证实了本文海洋条件模型建立和程序修改的正确性。

4 结论

1) 本文通过建立不同运动条件下流体的惯性加速度模型，基于RELAP5程序开发得到适用于复杂海洋条件的反应堆热工水力系统分析程序。

2) 本文程序计算结果与摇摆运动下竖直通道内单相水低速流动实验结果比较吻合，不同运动条件下程序的结果趋势与RETRAN-02/GRAV和RELAP5/MC程序符合较好。

未来在海洋条件本构关系式模型修改和实验验证方面还需开展进一步的工作。经过后续不断改进验证，该程序可成为我国浮动核电站设计研发及安全评审的有力工具。