海上核动力平台是小型核反应堆与船舶工程的有机结合, 可应用于诸多领域。其热工水力特性随海洋条件而变动, 海中的风、浪等复杂的水文环境造成船体倾斜、横摇、纵摇以及上下起伏, 使船体形成6个自由度的运动。运动改变了反应堆瞬态参数, 平台反应堆在动态事故和自然循环时, 热工水力分析更为复杂。因此, 研究不同海洋条件对反应堆瞬态参数的影响对我国海上反应堆的优化设计及安全分析具有重要的现实意义。

国外很早就开展了海洋条件下热工水力的特性研究。MURATA H等人^[1]先后对摇摆条件下的自然循环能力和自然循环换热进行了实验研究和理论分析; NARUKO Y等人^[2]以日本第一艘核舰船"陆奥号"的试验为基础, 对绕轴摆动等条件对舰船反应堆热工水力特性的影响进行了实验和数值研究, 开发了适于海洋条件下热工水力分析的RETRAN02/GRAV程序。在国内, 高璞珍等人^[3]在小型实验台架上进行了运动条件对热工水力特性影响的实验, 这方面的研究目前国内都以他提出相关研究的数学模型为基础。谭思超等人^[4]通过实验研究了摇摆情况下单相自然循环流动的传热能力。我国相关部门和许多科研人员正加强该方面的研究, 但我国目前还无自主知识产权的热工水力分析软件。

现有的热工水力系统分析程序仅适用于陆基条件, 摇摆条件下核反应堆研究虽取得了一定成果, 但很少有开展海洋条件下核动力平台事故时瞬态分析的研究, 对海洋条件下海上核动力平台的全面分析也有待加强。本文通过在动量方程中增加附加力模型, 并基于摇摆对传统热工水力系统各控制体的影响, 修改了程序模型坐标, 实现了程序的二次开发; 针对海洋小堆的特点建立了简化的系统模型, 并对正常运行时海洋条件的影响进行了计算分析; 对两种尺寸的破口事故开展了海洋条件下的模拟和分析。

摇摆条件下, 船相对于地球(惯性系)加速运动, 非惯性系下质量方程及需引入作用在流体上的附加力的动量方程分别为

$ \frac{\partial }{{\partial t}}({\alpha _k}{\rho _k}) + \frac{1}{A}\frac{\partial }{{\partial x}}({\alpha _k}{\rho _k}A{\kern 1pt} {\kern 1pt} {V_k}) = {F_k} $

(1)

$ \begin{array}{l} {\alpha _k}{\rho _k}A\frac{{\partial {v_k}}}{{\partial t}} + \frac{1}{2}{\alpha _k}{\rho _k}A\frac{{\partial v_k^2}}{{\partial x}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\alpha _k}A\frac{{\partial P}}{{\partial x}} + {\alpha _k}{\rho _k}({B_x} + {F_{{\rm{roll}},k}})A + {S_k} \end{array} $

(2)

式中:t——时间;

α_k, ρ_k, V_k——气相或液相占比、密度、速度, 气液两相下标k分别为g和f;

A——管道横截面积;

x——流动方向位置;

B_x——惯性坐标系下流动方向上质量力;

F_{roll, k}——附加力;

S_k——对流动时壁面及界面摩擦力、流体界面间动量互换、假想质量力求和。

根据附加力模型^[5]

$ {F_{{\rm{ roll }},k}} = {\alpha _k}{\rho _k}A{\kern 1pt} {a_{{\rm{ roll }},k}} $

(3)

$ \begin{array}{*{20}{l}} {{a_{{\rm{roll}}}} = {a_{\rm{t}}} + {a_{\rm{c}}} + {a_{\rm{u}}} = {\beta _{{\rm{roll}}}}r + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\omega _{{\rm{roll}}}}({\omega _{{\rm{roll}}}} \times r) + 2{\omega _{{\rm{roll}}}}u} \end{array} $

(4)

式中:a_t, a_c, a_u——摇摆产生的切向加速度、向心加速度和科氏加速度;

θ_roll, ω_roll, β_roll——摇摆的角位移、角速度和角加速度;

r——冷却剂至摇摆中的点的距离;

u——冷却剂相对惯性系的流速。

摇摆条件下, 船体偏离竖直状态的角度随时间变化, 且流体受附加力作用下, 舰船摇摆运动也通常为简谐运动。摇摆的角位移θ_roll、角速度ω_roll和摇摆的角加速度β_roll呈正弦变动。

惯性坐标系下, 运动的空间坐标在绕x轴, y轴, z轴摇动时改变, 物体垂直高度也会发生变化。摇摆条件是指船体绕某一轴线做往复运动, 会改变控制体在惯性坐标系下的空间坐标。若要描述摇摆后的控制体位置, 则需进行控制体坐标的旋转变换。起点为P和长为L的物体的三维坐标示意如图 1所示。

其方位角和倾斜角分别为φ和θ, 则P点坐标为

$ \left\{ {\begin{array}{*{20}{l}} {x = L\cos \theta \cos \varphi }\\ {y = L\cos \theta \sin \varphi }\\ {z = L\sin \theta } \end{array}} \right. $

(5)

此处只推导绕x轴的坐标变换, 如图 2所示。

点P在yz平面的投影是P_yz, Φ_xo为OP_yz和y轴的夹角, P_yz以向右旋转为正向绕x轴旋转得P_yz'。进行三角变换即可获得旋转后的新坐标(x_x, y_x, z_x), 写成矩阵形式为

$ {\left[ {\begin{array}{*{20}{c}} {{x_x}}\\ {{y_x}}\\ {{z_x}} \end{array}} \right]_t} = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&{\cos {\phi _x}}&{ - \sin {\phi _x}}\\ 0&{\sin {\phi _x}}&{\cos {\phi _x}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y\\ z \end{array}} \right] $

(6)

同理可证明(x, y, z)依次绕3个坐标轴分别旋转Φ_x, Φ_y, z角度后的新坐标(x', y', z')为

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{x^\prime }}\\ {{y^\prime }}\\ {{z^\prime }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos {\phi _z}}&{ - \sin {\phi _z}}&0\\ {\sin {\phi _z}}&{\cos {\phi _z}}&0\\ 0&0&1 \end{array}} \right] \times \left[ {\begin{array}{*{20}{c}} {\cos {\phi _y}}&0&{\sin {\phi _y}}\\ 0&1&0\\ { - \sin {\phi _y}}&0&{\cos {\phi _y}} \end{array}} \right] \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&{\cos {\phi _x}}&{ - \sin {\phi _x}}\\ 0&{\sin {\phi _x}}&{\cos {\phi _x}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} x\\ y\\ z \end{array}} \right] \end{array} $

(7)

热工水力系统程序中没有计算几何结构的绝对坐标, 只有用户输入的每个控制体在x, y, z 3个方向的长度和方位角, 但不影响摇摆条件的坐标变换。控制体在3个方向的长度可表示为$\left[ \begin{array}{l} {x_2} - {x_1}\\ {y_2} - {y_1}\\ {z_2} - {z_1} \end{array} \right]$, 而经坐标转换后$\left[ {\begin{array}{*{20}{l}} {{x_2}^\prime }\\ {{y_2}^\prime }\\ {{z_2}^\prime } \end{array}} \right] = \mathit{\boldsymbol{A}}\left[ {\begin{array}{*{20}{l}} {{x_2}}\\ {{y_2}}\\ {{z_2}} \end{array}} \right], \left[ {\begin{array}{*{20}{l}} {{x_1}^\prime }\\ {{y_1}^\prime }\\ {{z_1}^\prime } \end{array}} \right] = \mathit{\boldsymbol{A}}\left[ {\begin{array}{*{20}{l}} {{x_1}}\\ {{y_1}}\\ {{z_1}} \end{array}} \right]$, 则$\left[ \begin{array}{l} {x_2} - {x_1}\\ {y_2} - {y_1}\\ {z_2} - {z_1} \end{array} \right]$可通过相同的旋转矩阵A进行高度差转换。假设某点P=(x₀, y₀, z₀), 绕任一单位向量轴e=(e_x, e_y, e_z)旋转θ₀(平衡位置时的倾角)后, 新坐标P'=(x', y', z')可以由式(8)计算得到。

$ \left[ {\begin{array}{*{20}{c}} {{x^\prime }}\\ {{y^\prime }}\\ {{z^\prime }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {e_x^2(1 - \cos {\theta _0}) + \cos {\theta _0}}&{{e_x}{e_y}(1 - \cos {\theta _0}) - {e_z}\sin {\theta _0}}\\ {{e_x}{e_y}(1 - \cos {\theta _0}) + {e_z}\sin {\theta _0}}&{e_y^2(1 - \cos {\theta _0}) + \cos {\theta _0}}\\ {{e_x}{n_z}(1 - \cos {\theta _0}) - {e_y}\sin {\theta _0}}&{{e_y}{e_z}(1 - \cos {\theta _0}) + {e_x}\sin {\theta _0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_0}}\\ {{y_0}}\\ {{z_0}} \end{array}} \right] $

(8)

程序中表示控制体x, y, z 3个方向的标高变量为h_x, h_y, h_z, 可对这3个变量进行修改。另外程序采用弧度制, 在计算中需注意弧度制与角度制的转换。程序需提供外部接口, 由用户输入旋转轴单位向量和倾斜角。

摇摆时, 对于几何坐标的变化只需在原坐标基础上左乘式(8)所示的矩阵, 同时摇摆运动使旋转角θ_roll随时间周期不断变化。

$ {\theta _{{\rm{ roll }}}} = {\theta _0} + {\theta _{\rm{a}}}\sin \frac{{2{\rm{ \mathsf{ π} }}t}}{{{T_{\rm{r}}}}} $

(9)

式中:θ_a, T_r——摇摆幅值和摇摆周期。

摇摆产生的附加作用力^[6]如式(3)和(4)所示。向量方向为流体流动方向。在对程序进行改造时, 可将F_{roll, k}加入动量方程等式右侧。式(4)中, 由于科氏加速度在流动方向上的分量为零, 因此一维程序中不予考虑。摇摆的角速度ω_roll和摇摆的角加速度β_roll的计算式分别为

$ {{\omega _{{\rm{ roll }}}} = \frac{{2{\rm{ \mathsf{ π} }}{\theta _{\rm{a}}}}}{{{T_{\rm{r}}}}}\cos \frac{{2{\rm{ \mathsf{ π} }}t}}{{{T_{\rm{r}}}}}} $

(10)

$ {{\beta _{{\rm{ roll }}}} = - \frac{{4{\kern 1pt} {\kern 1pt} {{\rm{ \mathsf{ π} }}^2}{\kern 1pt} {\theta _{\rm{a}}}}}{{T_{\rm{r}}^2}}\sin \frac{{2{\rm{ \mathsf{ π} }}t}}{{{T_{\rm{r}}}}}} $

(11)

若反应堆系统的位置矢量r=(x, y, z), 摇摆轴的单位矢量为e=(e_x, e_y, e_z), 则切向作用f_t和向心作用f_c在笛卡尔坐标轴上的分量形式^[6]为

$ {\left[ {\begin{array}{*{20}{l}} {{f_x}}\\ {{f_y}}\\ {{f_z}} \end{array}} \right]_{\rm{t}}} = {\beta _{{\rm{ roll }}}}\left[ {\begin{array}{*{20}{c}} 0&{{e_z}}&{ - {e_y}}\\ { - {e_z}}&0&{{e_x}}\\ {{e_y}}&{ - {e_x}}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{l}} x\\ y\\ z \end{array}} \right] $

(12)

$ {\left[ {\begin{array}{*{20}{c}} {{f_x}}\\ {{f_y}}\\ {{f_z}} \end{array}} \right]_{\rm{c}}} = \omega _{{\rm{roll}}}^2\left[ {\begin{array}{*{20}{c}} {e_y^2 + e_z^2}&{ - {e_x}{\kern 1pt} {e_y}}&{ - {e_x}{\kern 1pt} {e_z}}\\ { - {e_x}{\kern 1pt} {e_y}}&{e_x^2 + e_z^2}&{ - {e_y}{\kern 1pt} {e_z}}\\ { - {e_x}{\kern 1pt} {e_z}}&{ - {e_y}{\kern 1pt} {e_z}}&{e_x^2 + e_y^2} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y\\ z \end{array}} \right] $

(13)

由于附加力计算需用到旋转力矩, 用来计算反应堆系统相对以旋转轴起点为原点的空间坐标, 热工水力系统程序内部中并不计算坐标, 因此需要增加这部分内容。摇摆时需要提供的用户接口是旋转轴单位向量、摇摆的振幅和周期。

由于缺乏实验数据, 受已公开文献中相关计算模型参数完整性的限制, 本实验中仅通过对比相似系统的计算结果变化趋势验证改造的程序。

目前在压水堆核动力系统中, 反应堆参数的不同及应用背景的不同等因素, 可导致其系统结构存在些许差别, 但总的构成基本相同, 即一回路通常由1个反应堆和2~4个环路组成。2个环路的各参数、运行方式、瞬态特性等均为独立模块, 可模拟双环路对称运行瞬态或不对称运行瞬态, 在事故瞬态分析时也可采用完好环路与事故环路的两环路分析方法, 以更好地模拟双环路实际的各种工况。对于双环路以上的多环路情况, 可将所有完好回路等效为一个环路, 另一环路为事故环路, 或将参数及运行、事故方式相同的环路等效为一个环路。这种可等效为2个环路进行仿真的实际工况覆盖了压水堆核动力运行及事故瞬态的很大范围, 极大地扩充了程序的功能和使用范围。为分析海洋条件对海上核动力小堆热工水力参数的影响, 本文采用改造相关程序建立了2个回路的简化小堆系统, 简化的海洋小堆建模节点如图 3所示。

图 3中, 100P和热构件模拟堆芯部分, 106和403P及对应的206和503P模拟小堆上常用的直流式蒸汽发生器, 109和209为2个回路的主泵, 水力学部件401-407和501-507模拟二回路。

简化后海洋小堆设计参数如表 1所示。

首先, 采用改造的热工水力系统程序对主泵驱动的满功率运行的正常工况进行计算, 分析摇摆运动对流场和温度场的影响。主要研究绕y轴摇摆; 取摇摆幅值30°, 周期12 s的横摇运动。有无摇摆时瞬态参数对比如图 4所示。

由图 4可知, 在无海洋条件下的各参数均在系统启动后逐渐趋于稳定直至无变化, 与陆地上的情况一致。在未发生事故或无外界因素影响时, 核反应堆系统稳定工作, 各热工水力学参数是恒定不变的。

绕y轴摇摆时, 系统流量呈周期性变动, 如图 4(b)所示, 横摇对堆芯总流量的影响小于对两个支路影响, 堆芯总流量波动幅值远小于环路内流量波动幅值; 堆芯总流量变动周期近似摇摆周期的1/2, 海洋小堆两侧支路内的流量波动周期是摇摆运动的周期, 流量相对于摇摆角度存在一定的滞后, 即相位差, 且2个支路内流量近似反相波动。实验中模拟多组工况发现, 绕轴摇摆时流量变化近似正弦波动, 随周期减小, 摇摆振幅增加, 流量波形会渐渐偏离正弦函数。其主要原因是摇摆引起的附加力导致在流动方向上的附加加速度会对流体产生附加压降, 造成系统参数出现周期性的波动。附加压降又取决于加速度(切向加速度a_t和向心加速度a_c), 其运动符合正弦函数规律。此外, 绕y轴摆动对堆芯和蒸发器温度的影响并不明显, 温度波形较复杂, 如图 4(c)所示。

本文对绕x轴摆动的工况也进行了计算。由于环路无垂直方向的高度变化, 因此对流量、温度等参数影响很小。

失水事故是具有代表性的核反应堆设计基准事故。虽然海洋条件引起的附加作用压头远小于泵的驱动压头对流场和温度场的影响, 可能导致海洋条件对由主泵作用下的强制循环的热工水力特性影响较小, 但对事故工况后的停堆或停泵情况, 以及如破口事故初期的临界流、事故后期的再淹没等特殊现象的影响未知。本文为测试海洋条件对反应堆事故流场、温度场以及安全性的影响, 选择了两种尺寸的破口事故进行计算分析。海洋条件下具体参数选取可以参考文献[7]中对舰船多自由度摇摆的规定。

考虑简化的破口事故情况, 假想1 000 s时在反应堆冷管段发生破口, 破口直径为5 cm和10 cm, 属小中破口事故。在图 3控制体110左侧加入时间控制体601和trip阀602表示发生破口; 在图 3控制体110和210下方分别加入时间控制体603和703及时间连接604和704表示安注系统的注入。两种尺寸破口事故时冷热段管道半径设定相同, 为0.549 9 m。为简化系统, 未考虑辅助系统如化容系统的压力补偿。假设当破口事故发生后, 压力降低到12.82 MPa下时依次发生停堆、停泵和安注, 时间上有3.5 s的延迟。试验中模拟了系统压力、堆芯和环路的质量流量、冷却剂温度、堆芯包壳温度、泵的转速、破口处质量流量、安注系统流量等参数, 下面将基于失水事故模型进行分析。

为了便于后续海洋条件的分析, 此处安注流量选择10 kg/s。图 5和图 6分别为无海洋条件和摇摆时的5 cm破口瞬态参数图。

破口事故时, 因流体质量和能量耗损, 压力骤降, 随后降低的速率减缓, 下降到一定程度后, 在稳定值内出现细微波动。无海洋条件时, 加热段出口、冷段1出入口、冷段2出入口冷却剂温度下降到2 200 s后有短暂的上升过程, 加热段出口的温升最大, 两冷段的入口温升幅度远高于出口, 但温升后都未超过事故发生前的温度, 见图 5(a)。无海洋条件时, 节点10的堆芯包壳温度下降到2 200 s后温度突增至峰值, 约800 K, 后骤降到与其他节点相同的温度值400 K左右; 节点9温度也出现相似现象, 但时间极短, 温升小, 见图 5(b)。这是因为事故中堆芯裸露, 且在安注系统当前10 kg/s流量下, 堆芯顶端会烧干, 导致包壳温度在冷却过程中有所上升, 但随着安注系统的投入, 最终可冷却至饱和温度。包壳峰值温度(Peak Cladding Temperature, PCT)是失水事故安全性的衡量标准, PCT主要取决于堆芯的裸露深度和裸露持续时间。摇摆时, 破口事故中高温冷却剂流出, 安注系统投入使冷却剂总温度下降, 当流量减少到一定值时堆芯裸露, 包壳温度迅速提高, 摇摆使堆芯裸露时间变长, 冷却剂温度飞升, 随着安注系统注水使包壳再次被淹没, 包壳温度下降, 冷却剂温度随之下降, 见图 6。其他破口事故瞬态参数与无海洋条件时变化趋势相同。

安注流量选择30 kg/s进行分析。图 7和图 8分别为无海洋条件和摇摆时的10 cm破口瞬态参数图。

对比破口直径为5 cm时的数据可知, 整体参数变化趋势基本未变, 但破口直径的增大使同一条件下的瞬态参数都有上升, 见图 7; 当破口直径增大, 破口处质量流量越大, 系统压力变化越快, 压力容器内装水量快速减少, 因泄压较快, 所以压力容器内冷却剂发生闪蒸较早。闪蒸所产生的蒸汽通过破口流出, 临界流降低, 破口处降压速度减缓, 堆芯空泡份额变大, 堆芯传热恶化, 使得堆芯包壳温度上升。之后, 由于安注系统的投入, 使得压力容器内蒸汽凝结, 压力减小, 堆芯水位上涨, 包壳温度降低, 见图 7(b)和图 8(b)。破口直径越大, 堆芯包壳峰值温度越高, 对反应堆各参数的影响也越大; 海洋条件下, 摇摆幅度越大, 反应堆各参数变化越明显。

本文开展了摇摆条件下海洋小堆的程序改造和建模分析, 并研究了其在事故工况下的特性。主要结论如下。

(1) 在靠密度差驱动的自然循环中, 海洋条件影响较为明显。摇摆会造成流量周期性的变化且2个支路流量波动存在近180°的相位差, 波峰波谷相互补偿, 使总流量变化减小。

(2) 在强制循环为主的正常运行工况下, 泵的驱动压头远大于海洋条件附加力的作用, 因此海洋条件对流场和温度场的影响较小。

(3) 对于破口事故, 摇摆作用带来的周期附加力会对原流量进行补偿或叠加, 所以影响更为复杂。海洋条件对温度的影响较小, 但周期性变化的摇摆工况会带来较大的流量波动, 造成剧烈的不稳定性。