引言

海洋核动力平台是海洋中的可移动式小型核电站，由船舶工程和核反应堆有机结合而成，可为海水淡化、能源开采和岛屿开发等海上工程作业提供高效而且清洁的能源。随着科技的进步，对海洋开发利用的关注程度不断提高，海洋核动力平台成为世界各国在核电领域和船舶领域的研发热点。有关资料显示，2016年以来，国内相关示范工程项目陆续设立，小型堆示范项目的初步设计工作也正在开展^[1]。

液体晃动是指贮液容器内的液体介质在外界载荷作用下产生运动的现象，环境载荷是由自然环境直接或间接产生进而作用在结构物上的载荷。海洋核动力平台长期在海洋上工作，工作海域的环境载荷对平台产生影响，在风载荷、波流载荷及波浪载荷单独或共同作用下，平台的晃动带动稳压器的晃动从而引起内部液体的晃动，对稳压器的整体结构受力及支承结构的疲劳寿命产生影响，并对反应堆的压力调节系统造成冲击。

稳压器是核电反应堆压力安全系统的主要设备，每个反应堆需设置一台稳压器，用来保持反应堆运行过程中的压力稳定，其对反应堆的安全运行至关重要。目前，常用的电加热式稳压器上部为汽空间、下部为液空间，在内部设置了喷雾器和电加热器，以完成对压力的调节和控制；同时设有泄压箱、安全阀等附件，用于一回路的压力调节和保护。当一回路压力有升高趋势时，通过喷雾器进行喷淋使水汽凝结以降压；当一回路压力有下降趋势时，通过加热器进行加热使水蒸发以升压。

研究液体晃动问题主要有理论分析^[2]、试验研究^[3-4]、数值模拟^[5]3种方法，对于海洋平台主要采用数值模拟的方法，在条件允许情况下辅以试验验证。目前已有文献研究得到了现有平台在渤海海域的水动力响应情况^[2]，但仍缺乏对平台上稳压器内液体晃动情况的研究。本文针对稳压器存在晃动的问题，首先通过查阅相关资料，确定海洋核动力平台在渤海海域的水动力响应；再运用Fluent软件建立流体仿真模型并定义稳压器的运动后，利用流体体积(VOF)模型进行数值模拟，通过研究不同工况下液体的晃动规律得到内部壁面平均静压强、内部压力等参数的变化情况及波动范围；最后对一种新型防晃液装置进行了仿真和分析，以验证其防晃效果。

1 稳压器的运动及有限元模型 1.1 平台工作环境及数据来源

海洋核动力平台使用单点系泊装置(single point mooring equipment)使其在海上固定。平台的主要参数见表 1。

海洋环境复杂多变，实际海浪、海流及风的波动具有不规则波的特点。平台经受不同风浪流组合作用，产生纵荡、横荡、垂荡、纵摇、横摇、艏摇6个自由度的动力响应。

文献[3]采用5种风浪流组合工况(表 2)，得到平台的运动响应(表 3)。采用渤海海域百年一遇的波浪工况，即有义波高4 m，有义周期11.04 s；海洋环境风速50 m/s；水流取渤海湾重现期为万年的流速2 m/s；运动形式为正弦。

1.2 欧拉刚体运动方程

一般小型核电反应堆稳压器相对于海洋核动力平台的尺寸很小，又都与平台刚性联接，因此，可把稳压器看成刚体，而且先不考虑液体的阻尼。

刚体的运动包括平动、定轴转动、平面运动、定点运动和一般运动等形式。在刚体上选定基点O后，刚体一般运动可以分解为刚体随基点的平动和相对基点平动系的定点运动。其中刚体随基点的平动用基点O相对于参考系原点O₀的矢径R₀来描述，刚体相对于基点平动系的定点运动用固连坐标系Oxyz关于平动坐标系OXYZ的方向余弦矩阵A描述，因此刚体的一般运动可以由矢径R₀和矩阵A完全确定^[6]。方向余弦矩阵A中包含3个相互独立的角度φ、ψ、θ，称为欧拉角。选取3-2-1欧拉角旋转方式，即先绕Oz轴转动，再绕Oy轴转动，最后绕Ox轴转动。刚体运动参考系如图 1所示。

平台的欧拉角形式运动方程为

$ \left\{ \begin{array}{l} \varphi = \frac{{5.54{\rm{ \mathsf{ π} }}}}{{{\rm{}}180}}{\rm{sin}}\left( {\frac{{2{\rm{ \mathsf{ π} }}}}{{11.04}}t} \right)\\ \psi = \frac{{0.55{\rm{ \mathsf{ π} }}}}{{180}}{\rm{sin}}\left( {\frac{{2{\rm{ \mathsf{ π} }}}}{{11.04}}t} \right)\\ \theta = \frac{{10.99{\rm{ \mathsf{ π} }}}}{{180}}{\rm{sin}}\left( {\frac{{2{\rm{ \mathsf{ π} }}}}{{11.04}}t} \right) \end{array} \right. $

(1)

方向余弦矩阵为

$ \mathit{\boldsymbol{A}} = ({\mathit{\boldsymbol{A}}_1}\;\;{\mathit{\boldsymbol{A}}_2}\;\;{\mathit{\boldsymbol{A}}_3}) $

(2)

式中

$ {A_1} = \left( {\begin{array}{*{20}{c}} {{\rm{cos}}\psi {\rm{cos}}\theta }\\ {-{\rm{sin}}\psi {\rm{cos}}\theta }\\ {-{\rm{sin}}\theta } \end{array}} \right) $

$ {A_2} = \left( {\begin{array}{*{20}{c}} {{\rm{sin}}\psi {\rm{cos}}\varphi + {\rm{cos}}\psi {\rm{sin}}\theta {\rm{sin}}\varphi }\\ {{\rm{cos}}\psi {\rm{cos}}\varphi-{\rm{sin}}\psi {\rm{sin}}\theta {\rm{sin}}\varphi }\\ {{\rm{cos}}\theta {\rm{sin}}\varphi } \end{array}} \right) $

$ {A_3} = \left( {\begin{array}{*{20}{c}} {{\rm{sin}}\psi {\rm{sin}}\varphi-{\rm{cos}}\psi {\rm{sin}}\theta {\rm{cos}}\varphi }\\ {{\rm{cos}}\psi {\rm{sin}}\varphi + {\rm{sin}}\psi {\rm{sin}}\theta {\rm{cos}}\varphi }\\ {{\rm{cos}}\theta {\rm{cos}}\varphi } \end{array}} \right) $

刚体定点转动角速度在固定坐标系中的坐标阵为

$ \mathit{\boldsymbol{\omega }} = \left[\begin{array}{l} {\omega _X}\\ {\omega _Y}\\ {\omega _Z} \end{array} \right] $

(3)

式中

$ \left\{ \begin{array}{l} {\omega _X} = \dot \varphi-\dot \psi {\rm{sin}}\theta \\ {\omega _Y} = \dot \psi {\rm{cos}}\varphi {\rm{cos}}\theta-\dot \theta {\rm{sin}}\varphi \\ {\omega _Z} = \dot \psi {\rm{sin}}\varphi {\rm{cos}}\theta-\dot \theta {\rm{sin}}\varphi {\rm{tan}}\varphi + \dot \theta {\rm{sec}}\varphi \end{array} \right. $

这样通过欧拉刚体运动方程，定义了稳压器的刚体运动规律。

1.3 有限元模型

Fluent系列软件中多相流模型包括VOF模型、Eulerian模型以及Mixture模型。本文使用多相流模型VOF方法模拟自由液面晃动，使用用户自定义函数(UDF)中的DEFINE_CG_MOTION宏来定义稳压器的刚体运动。在此宏中，通过用户提供每个时间步长的线速度和角速度来指定Fluent中某一特定区域的运动，Fluent利用这些速度值来更新模型的节点位置。

稳压器有限元分析模型参数见表 4。稳压器材料为SA50 8- 3钢，材料属性见表 5^[7]。容器内部工作压力15.4 MPa，工作温度345 ℃。内部介质为水和蒸汽，水占容积的60%，蒸汽占容积的40%，黏性不可压。在工作状态下，水的密度为592.25 kg/m³，蒸气的密度为103.59 kg/m³。

流体仿真模拟对计算机性能要求较高，需要对有限元模型进行适当简化。在模拟不同工况下液体的晃动规律时，省略掉稳压器内部细小的构件不会影响整体液体晃动规律的比较，但有利于加速模拟。定义初始模型在软件中OYZ平面的正法向侧为前壁面，另一侧为后壁面，采用四面体方法进行有限元网格划分，共划分网格106 789个，正交品质(orthogonal quality)为0.23，网格质量较好。网格划分结果如图 2所示。

进行数值模拟计算时，网格的数量和质量会对计算结果产生影响，因此需要在计算前验证网格无关性，尽量排除网格对结果的影响。

定义水和蒸汽的分布域后，求稳压器内初始绝对压力的分布，并与理论值比较。定义操作密度为103.59 kg/m³，即以蒸汽相的密度为基准，这样可以改善计算的精度，提高模拟速度。理论绝对压力最大值(p_max)为

$ \begin{array}{l} \;\;\;\;\;\;{p_{{\rm{max}}}} = {p_{{\rm{op}}}} + \rho gh = 101\;325 + \left( {103.59-592.25} \right)\\ \left( {-9.81} \right) \times 3.42{\rm{ }} = 1.18 \times {10^5}\;{\rm{Pa}} \end{array} $

式中，p_op为操作压力，ρ为相对密度，h为液柱高度。

数值模拟出容器内绝对压力由顶部到底部呈线性分布，底部绝对压力最大值为1.18×10⁵ Pa，与理论值吻合，故流体网格合适。

2 数值模拟结果及分析

参考相关理论和文献[8]，对Fluent的求解器进行如下设置：瞬态模型；压力基求解器；标准的k-ε模型；Emplicits显式算法；Implicit Body Force隐式体力；Absolute绝对速度；1s t- Order Implicit一阶隐式格式；Body Force Weighted Pressure；First Order Upwind Momentum。设时间步长0.011 04 s，步数5 000步，即模拟5个周期。

2.1 自由液面

根据稳压器的水位控制系统，当稳压器的水位值高出正常工况的10%时，发出高液位警报，此时备用电加热器启动，通过加热使部分水蒸发；当水位值低于正常工况的5%时，发出低水位警报。晃动时如果液体的晃动超出规定范围，则触发报警系统。

根据软件的坐标系统和原点设置，稳压器模型在初始状态下的水位Y值为420 mm，高水位Y值为762 mm，低水位Y值为249 mm，如图 3所示。

由于稳压器模型存在倾斜和摇摆运动，而Fluent软件中的坐标系是固定的，因此警戒水位应为相对于初始液面的水位高度，即h=762-420=342 mm。

液体晃动过程中，每10步采集一次稳压器内部自由液面晃动数据，观察不同工况下自由液面变化情况，结果如图 4所示。

由图 4看出，工况1和工况2下液面晃动没有超过警戒水位；工况3启动时相对液面最高达1 574.54 mm，后续相对最高液面392.29 mm，液面晃动均超过了警戒水位；而工况4、5启动时相对液面最高达870.00 mm和967.97 mm，液面晃动超过了警戒水位，后续相对最高液面分别为275.46 mm和274.91 mm，没有超过警戒水位。说明在渤海海域工作的稳压器会因晃动引起水位控制系统的频繁开启。

2.2 壁面静压强的变化过程

在液体的晃动过程中，液体与璧面产生相互作用力，作用力垂直于物体表面，由流体与璧面的压差产生，在Fluent中为静压强^[7]。通过数值模拟记录不同工况下的晃动全程中稳压器前、后壁面静压强平均值(p)如图 5所示。分析图中壁面静压强的变化曲线，并结合自由液面的晃动过程模拟结果(图 4)，可以看出在运动的初始阶段液体的晃动比较剧烈，随着流动时间的延长，液体的晃动逐渐呈规律变化。

稳压器壁面所受冲击力的大小为静压强在面积上的积分，5种工况下稳压器前后壁面所受冲击力平均值见表 6，可见液体晃动对前后壁面造成的冲击力大致相等。

2.3 稳压器内部压力的波动

根据稳压器的运行设定，当稳压器内部压力上升至15.6 MPa时，喷雾器开启运行；压力达到15.93 MPa时，喷淋阀全开；压力达到16.12 MPa时，泄压阀开启；压力高达16.37 MPa时，反应堆紧急停堆。稳压器内部压力低于额定工作压力时，电加热器开启运行；压力15.24 MPa时，备用组电加热器投入运行；压力14.7 MPa时，则发出低压预警信号；低至12.25 MPa时，反应堆紧急停堆。

定义Fluent中操作压力为工作压力，即15.4 MPa。以Fluent中操作压力为参考值，计算所得压力为相对压力值，这样可以提高软件计算精度和速度。

在运动的启动阶段，稳压器内部的瞬间压力较大，随着流动时间的增加，稳压器内部的压力逐渐呈正弦规律变化。在Fluent中静压强即为表压。提取5种工况下启动时(第一个记录点0.011 04 s)的瞬间压力值和规律时的压力峰值，如表 7所示。可以看出在晃动最剧烈的工况3下，启动时的瞬间压力超过了0.2 MPa，即晃动的影响可能会造成喷雾器的开启。

3 防晃液装置数值模拟

通过对不同工况下稳压器内液体晃动的数值模拟发现，若只考虑环境载荷，液体的晃动在可控范围内；但由于海洋工况的复杂多变，在冲击载荷作用下稳压器内液体晃动会较剧烈，所以可以考虑在稳压器内设置防晃液装置以应对极端情况。

稳压器喷雾器的喷淋过程以及电加热器的加热过程都需要保证其内部液体分布均匀，同时还要考虑到内部的汽液转换过程。郝天才等^[9]提出了一种新型的防晃液装置，形式如图 6。该防晃液装置与稳压器壁面采用焊接或者焊接加螺栓联接的形式；在水平挡板上均布有圆孔，保证了汽液蒸发凝结过程的顺利发生。

将防晃液装置安装在初始自由液面处后，通过数值模拟，得到工况3下内部液体不同时刻的晃动效果如图 7所示(黑色表示液相，灰色表示气相)。由图可以直观看出，液体晃动十分平稳。

数值模拟结果如表 8所示。可以看出前壁面最大静压强由6 512 Pa下降到4 490 Pa，后壁面最大静压强由8 204 Pa下降到4 966 Pa，内部静压强由29 813 Pa下降到20 730 Pa，表明稳压器结构整体受力峰值减小，壁面受力更加平均。启动时容器内部瞬间压力为0.11 MPa，证明防晃液装置有效抑制了液体的晃动。

4 结论及展望

(1) 通过模拟得出在平台与风载荷呈180°、与波浪载荷呈180°、与波流载荷呈90°时，稳压器内部液体的晃动最为剧烈。

(2) 同工况下内部介质的晃动表现为初始时较大，属于冲击响应；随后逐渐呈规律、平稳变化；波动形式主要为驻波和行进波，在受冲击时会有水跃现象。

(3) 极端工况下，稳压器内部压力的波动仍在可控范围内，但瞬态冲击情况下液体晃动较为剧烈。因此，在防晃液装置研究中，应以瞬态响应为研究重点。

(4) 初步模拟分析稳压器内部液体防晃装置，表明该型防晃液装置可以有效抑制液体晃动幅度。本文模拟分析结果将为海洋核动力平台反应堆系统的安全运行、稳压器结构强度设计、内部防晃液装置的优化及流固耦合分析等提供一定的理论基础。