0 引　言

当外部激励作用于储液容器时造成容器内部流体运动的现象称为晃荡，该现象具有很强的非线性和随机性，其产生的冲击力是液舱结构破坏和船舶失稳的主要原因之一。随着各种液化船的广泛使用，液舱晃荡现象越发成为人们关注的焦点。目前，对于晃荡的研究主要有理论分析方法、数值方法和实验方法3种方法。其中数值方法由于其适应性强、成本低、效率高的特点受到广泛欢迎。因此本文采用基于OpenFOAM的CFD技术进行数值模拟。

晃荡问题最早出现在航天航空领域，1966年Abramson^[1]对球型和矩形容器内的液体晃荡问题开展了相关研究，揭示了晃荡对于燃料箱稳性的影响。Abramosn是研究液舱晃荡问题的先驱，但由于当时条件的制约，该研究无论是结果精度还是条件考虑都有所欠缺。随着计算机和计算流体力学（CFD）的发展，晃荡研究的缺陷逐步完善。Faltinsen等^[2]对于二维情况下不同形状液舱的晃荡进行了分析，揭示了二维情况下方形液舱晃荡载荷特性。黄硕等^[3]对于二维液舱在非线性水动力影响下的运动特性进行了研究,得到了能量耗散系数的选取规律。孙士艳等^[4]对二维楔形储罐内液体进入平静水面的水动力问题进行了研究，揭示了二维楔形槽内液体进入平静水面的运动特性。张书谊等^[5]基于CFD软件Fluent探索了二维情况下矩形液舱内液体水动力特性，但是没有考虑液舱长度的影响。Kim^[6]和朱仁庆等^[7]对船舱内液体的三维晃荡特性进行了分析，揭示了三维情况下液舱晃荡特性。王庆丰等^[8]使用HydroSTAR对带液舱的FSRU船舶进行运动响应分析，揭示了液舱载液率和液舱布置对FSRU运动的影响。甄长文等^[9]对于三维共振频率下油船液舱舱内液体晃荡情况进行了模拟研究，发现当晃荡稳定后舱内液体以驻波与行波组合的方式运动。李裕龙等^[10]基于OpenFOAM对船舶与晃荡在波浪中的时域耦合运动进行研究，揭示了船舶在液舱晃荡时的横摇特性。Clearly^[11]和Bass等^[12]综合使用Matlab与CFD软件针对不同形式液舱液体晃荡进行了数值计算比较，并针对矩形舱室提出了抑制晃荡的方法。但上述研究皆默认水作为介质的情况，SoutoIglesias^[13]首先对不同宽度液舱在不同介质中的晃荡载荷进行研究，分析了液舱宽度以及介质变化引起的三维效应对冲击载荷的影响。在此之后，Zou^[14]系统地研究了介质粘性变化对于冲击载荷的影响，发现高粘度液体的冲击压力在黏性阻尼的作用下需要更长的时间才能达到稳定状态。Lee等^[15]模拟了流体粘度、密度比以及可压缩性对晃动载荷的影响等。

国内外已经取得了关于液舱晃荡的大量研究成果，但是很少有学者开展当液舱内流体状态发生变化时船舶的运动特性研究。本文将研究液舱内水冻成冰时，船舱的运动特性，基于对结构体的固有频率和共振特性分析，深入研究流体状态的发生改变时内部水或冰、船舱和外部波浪之间的耦合作用机理。

1 数学模型

本文采用OpenFOAM开源软件，建立液舱在波浪中运动的CFD数值计算模型，研究舱内流体力学特性对液舱运动和受力的影响。为进行对比分析，假设液舱内水被冻成冰块，忽略热力学热胀冷缩规律，在建模时让水体与冰块密度相同，质量相同，并且让冰块形状与液舱未发生运动时水体的初始形状相同。建立一个内部结构为矩形的二维液舱计算域示意图，如图1所示。在此基础上，定义一个笛卡尔空间固定坐标系 $ O - xyz $ ，坐标原点 $ O $ 固定在载水液舱质心处， $ x $ 轴沿水平方向， $ z $ 轴竖直向上， $ y $ 轴垂直于 $ O - xz $ 平面，如图2所示。液舱壁厚度a=0.25 m，舱宽b=2 m，舱高c=1 m，设置计算域长度d=42 m，水深e=10 m，液舱内水深h=0.35 m，液舱内底距坐标原点O距离为g=0.09 m，虚线代表防护栏用于防止外部液体进入液舱。假设液舱绕质心O旋转，且液舱只有绕 $ y $ 轴的旋转运动，其角速度用符号 $ \omega $ 表示。为了保证液舱在平衡位置时吃水为总体高度c的一半，液舱的质量设置为1250 kg。

忽略流场可压缩性，采用适用于不可压缩流的连续性方程和N-S方程来建立流场的基本控制方程。

$ \nabla \cdot {{\boldsymbol{U}}} = 0 ，$

(1)

$ \frac{{{\rm{D}}{\boldsymbol{U}}}}{{{\rm{D}}t}} = {{{F}}_b} - \frac{1}{\rho }\nabla p + \nu {\nabla ^2}{U}。$

(2)

式中： ${\boldsymbol{U}}$ 为速度矢量； $\nu $ 为流体运动学粘性系数； ${{F}_b}$ 为质量力； $\ \rho $ 是流体密度； $p$ 为压力。入口处设置5阶斯托克斯波入射条件，物面与水底满足无滑移边界条件。使用OpenFOAMv1812中的interFoam求解器，采用有限体积法（FVM）离散控制方程，流体体积法（VOF）捕捉自由面。

当物体运动时，液固交界面必然发生改变,如何处理物面附近网格则是一个难题。为此，本文在物体外部选定一个区域作为动网格区域，动网格求解器设置为dynamicMotionSolverFvMesh，该求解器是根据边界条件与扩散方程计算动网格区域网格的运动，不改变拓扑结构，较其他方法该方法适用于任何类型的小幅运动，并且能够保证较高的网格质量。动网格区域的内部边界设置为0.3，外部边界设置为1。图3为网格示意图，在布置网格时需要注意，在物体和自由面附近，采用精细网格，远离物体和自由面，网格逐步过渡到粗网格。

2 数值算例 2.1 收敛性分析

选择一个载水液舱对网格收敛性进行验证。液舱壁质量 $m = 1\ 250$ kg（不包含舱内水质量）、转动惯量 $ {I_{yy}} = {\text{756}}$ kg·m² 、旋转中心坐标位于 $ \left( {0,0,0} \right) $ 处，质心与旋转中心重合。外部波浪参数取波长 $\lambda = 14$ m，波高 $H = 0.2$ m。在数值计算过程中，计算结果的精度以及计算时间与网格密度息息相关。虽然增加网格数量可以提高计算结果的准确性，但会给计算机造成巨大的负担，甚至有时还会导致计算机崩溃。因此，在保证计算准确度的前提下，网格数量不宜过多。为此，一共选取3种网格布置方式，分别为粗糙网格44112、中等网格51664、精细网格83916，采用自适应时间步长，最大库朗数设置为0.5，计算结果如图5所示。由图可见，精细网格和中等网格的时历结果基本重合，而粗网格结果略有差异，说明本文的数值结果是单元格收敛的，中等网格的数值结果已经达到很好的数值计算精度，为此，选择中等网格作为本文的网格布置方式。

2.2 结果分析

在进行数值结果分析时，仍沿用收敛性分析中使用的船舱参数，即舱壁质量 $m = 1\ 250$ kg、转动惯量 $ {I_{yy}} = {\text{756}}$ kg·m²、旋转中心 $ \left( {0,0,0} \right) $ 。此时需注意这些参数不包括液舱内水的质量和转动惯量，也就是说液舱内水将以水动力的形式出现在运动方程之中。当水被冻成冰时，船体发生几点变化：1）将等价于水的冰质量加入到液舱壁质量之中，加入后冰舱总质量为2 300 kg；2）将等价于水的冰转动惯量加入到液舱壁转动惯量之中，加入后冰舱总转动惯量为 $ {I_{yy}} = {\text{970}}{\text{.58}} $ kg·m²；3）保证冰舱旋转中心与液舱旋转中心重合，但需要考虑水变成冰以后，冰舱质心将发生变化，冰舱质心为 $ \left(\text{0}，\text{0}，\text{0}\text{.03}\right) $ 。首先对液舱的固有频率进行分析，在静水中给液舱一个初始角位移，然后让其自由晃荡，直至平稳，然后对时历数据进行傅里叶频率分析，可以得到载水液舱1阶固有频率为0.27 s⁻¹（周期为3.7 s）、2阶固有频率为0.54 s⁻¹(周期为1.85 s)。载冰舱1阶固有频率为0.46 s⁻¹（周期为2.17 s）、2阶固有频率为0.92 s⁻¹（周期为1.08 s）。可以发现，水舱固有周期大于冰舱固有周期，这是由于冰块相当于是刚体，每个质点的加速度都相同，而对于水，每个质点的加速度都不相同，使得水舱即难加速度也难减速，最终导致水舱的固有周期大于冰舱。

2.2.1 波高为0.2 m时水舱和冰舱运动时历对比分析

图6给出了波高为0.2 m时不同波长条件下水舱和冰舱的运动时历曲线。波长分别取为8.4 m，11.2 m，14 m，16.8 m和22.4 m，基于5阶无限水深斯托克斯波波长与周期的色散关系，可以得到周期分别为2.294，2.662，2.982，3.291和3.782 s。在初始时刻，水舱和冰舱均处于静止状态，在波浪激励下开始运动，然后逐渐达到稳定。观察图6，可以发现一个规律，随着波长的增加，载冰舱运动幅值逐渐减小，而载水舱幅值逐渐增大。与波长相比，0.2 m波高对应的最大波陡只有0.024。因此，可以基于线性共振理论对问题进行分析。由图5可知，载水舱固有频率为3.7 s，在图6所计算的所有波长中，最大波长所对应的波浪周期最接近于载水舱固有周期，从共振角度出发，载水舱运动幅值将越来越大，图6的计算结果完全符合关于水舱的分析。而对于冰，结果则恰恰相反，即波长越大，幅值越低，这是由于冰的固有周期是2.17 s，与最小波长所对应的固有周期最接近。

2.2.2 波高为0.6 m时水舱和冰舱运动时历对比分析

图7给出了波高为0.6 m时不同波长条件下水舱和冰舱的运动时历曲线。波长分别取为8.4 m，11.2 m，14 m，16.8 m和22.4 m。随着波陡的增加，此时非线性波的周期分别为2.098，2.483，2.888，3.199和3.734 s。与图6类似，由于冰舱固有周期接近于小波长固有周期，液舱固有周期接近于大波长固有周期。因此对于小波长情况冰舱运动幅值大，对于大波长情况液舱运动幅值大。此外，计算结果也表现出一些非线性特性，与图6相比，图7的波高是图6的3倍，但在相同工况下，图7的运动幅值并未达到图6的3倍，约为2倍左右。其次由于水舱的固有周期较大，图中没有表现出明显的高频共振特征。但对比冰舱，由于其固有频率较高，周期较短，运动时历结果表现出高频运动特征，尤其是图7(e)，由于此时波浪周期约为固有周期的1.7倍，因此运动时历结果的波峰处出现双峰特征。由图6和图7可知，相同质量的水和冰，运动特性差异极大，水舱和冰舱最主要的变化特征是固有频率的变化，因此在实际工程计算中，应尽量避免液舱固有频率接近波浪频率，可极大降低船舶航行风险。

2.2.3 液舱内水的晃荡特性分析

图8给出以上2个算例的液舱内部自由面高度变化情况，取液舱内部监测点 $ x = 0.7 $ m处的自由液面变化情况。由图可见，当波高是0.2 m时，自由面高度出现稳定的周期性变化特征，因此可以判断，此时液舱内部晃荡比较平稳。当波高提高到0.6 m，非线性特性逐渐变得显著，自由面开始出现高频振动特征。在图8(c)和图8(d)中，可以看到明显的双峰特征；在图8(d)和图8(f)中，在一个周期内明显看到3个小波峰，这些都体现了液舱内部流体的高频率振动特征，也体现出在非线性效应明显的情况下，液舱内部流体的高频特征更容易体现出来。但也可以看到，无论波高为多少，内部自由液面的主频仍然与入射波浪频率一致。

3 结　语

本文基于OpenFOAM开源软件，从连续性方程和N-S方程出发，采用动网格技术，将载冰舱与载水液舱运动特性进行对比，研究液舱内水运动对于液舱运动特性的影响。

1）相同形状和相同密度水舱固有周期大于冰舱固有周期，这是由于冰块相当于是刚体，每个质点的加速度都是相同的，而对于水，每个质点的加速度都不相同，这使得水舱即难加速度也难减速，最终导致水舱的固有周期大于冰舱。

2）受船舱固有周期与波浪周期相对关系的影响，在本文算例中，当波长逐步增加时，载冰舱运动幅值逐渐减小，而载水舱幅值逐渐增大，其深层次的原因是，船舱固有频率和波浪频率接近时，船舱将发生共振，此时运动幅值最大。

3）当波高不大时，液舱内流体随波浪柔和振荡，随着波高的增加，非线性效应逐渐显现，液舱内流体的高频振动特征也随之出现，在一个波浪周期内，流体质点可能出现双峰或三峰振动特征。