0 引　言

液舱晃荡是指液舱中非满载液体在运动中冲击舱壁的现象。在船舶中，晃荡现象会对船舶的平衡性和稳定性产生不利影响，严重时可能会使船舶倾覆，因此，研究液舱晃荡对于工程的安全性十分重要。

研究晃荡问题的主要方法可分为理论研究、数值模拟、实验方法。在远离共振频率和小幅度晃荡运动中，可以采用线性势流理论对晃荡压力进行预报，然而对于处于大幅度和共振频率的晃荡运动中，自由液面表现出很强的非线性特征，线性势流理论也不再适用。Faltinsen^[1]和Ibrahim^[2]通过理论研究和实验方法，对晃荡的非线性现象和冲击压力展开研究，其成果对于研究晃荡有着重要意义。随着计算机的发展，数值模拟成为了研究晃荡的主要方法，对于模拟大幅度晃荡运动，自由液面捕捉技术是关键所在，常见的自由液面捕捉技术有ALE法（Arbitrary Lagrange-Euler）、VOF法（Volume of Fluid）、Level-Set法。周骥等^[3]通过ALE法准确模拟了液体自由液面运动，并适用于流固耦合问题的求解。孙丰等^[4]研究表明ALE法能很好捕捉物体着水时液体的飞溅现象。王勇学等^[5]使用VOF法研究圆柱液舱晃荡问题，数值模拟结果和实验结果大致吻合，证实了VOF法研究液舱晃荡的可靠性。邹昶方等^[6]使用VOF法研究液舱壁面受到的晃荡冲击压力，并进行三维模拟，对比实验结果，表明VOF法适用与三维数值模拟计算。Level-Set法能很好处理液体破碎翻卷等非线性运动，陶建华等^[7]利用Level-Set法有效的模拟溃坝波的传播，并给出了自由液面和流场速度结果。高玉丽等^[8]利用Level-Set法模拟计算瞬时溃坝案例，其计算结果和实验结果大致相同，可以证实Level-Set法在模拟液面破碎和融合和结构物体有拐角突变状态下，其计算结果的精确度具有可靠性。李师等^[9-10]通过VOF法研究了复合激励下液舱晃荡的压力特性以及激励频率对晃荡压力的影响。王德禹等^[11]提及了多个影响晃荡的因素和非线性特点，以及在激励幅值较大时，晃荡压力不总随外部激励的加强而加强，呈现出其“饱和特性”。众多学者也对抑制晃荡做了大量研究，邹昶方等^[12]根据相似原理设计了弹性支撑式液舱晃荡实验模型，探究了弹性支撑对晃荡压力和晃荡共振频率的影响，对比了刚、弹性支撑的晃荡压力。杜祥璞^[13]通过在壁面设置凸起结构，分析了液舱在不同装载率下，壁面结构形式对晃荡压力的影响。王晨屹等^[14]通过布置自由液面的多孔竖直构件，分析了液舱在低装载率下，复合式构件对晃荡压力抑制的效果，表明这种构件能减少晃荡对舱壁的冲击压力。Jiang等^[15]的研究表明液舱中增加垂直挡板高度可以更好的抑制晃荡现象，减小晃荡对舱壁的冲击压力。

目前大幅激励的晃荡现象研究甚少，本文针对二维矩形液舱在大幅单自由度和多自由度激励下引起的液舱晃荡问题，采用VOF法自由液面捕捉技术和CFD数值模拟方法，研究了液舱在不同装载率下，大幅横荡激励、大幅横荡-垂荡耦合激励的晃荡压力特性和变化规律，为大型LNG船液舱结构详细设计提供参考。

1 数值方法 1.1 数学模型

连续性方程是基于质量守恒原理，RANS方程是通过雷诺平均N-S方程进行时间平均得到，用来描述湍流运动^[14]，将2个方程结合起来可以提供液体晃荡运动行为的定量描述和预测。

1）连续性方程

$ \frac{\partial \rho }{\partial t}+\frac{\partial }{\partial {x}_{i}}\left(\rho {u}_{i}\right)=0 。$

(1)

式中：$ \rho $为密度，$ {u}_{i} $为$ i $方向的速度。

2）RANS方程

$ \rho \frac{\partial {\overline{u}}_{i}}{\partial t}+\rho {\overline{u}}_{j}\frac{\partial {\overline{u}}_{i}}{\partial {x}_{j}}=\rho {\overline{F}}_{i}-\frac{\partial \overline{p}}{\partial {x}_{i}}+\frac{\partial }{\partial {x}_{j}}\left(\mu \frac{\partial {\overline{u}}_{i}}{\partial {x}_{j}}-\rho \overline{{u}_{i}^{\mathrm{{'}}}{u}_{j}^{\mathrm{{'}}}}\right) 。$

(2)

式中：$ {\overline{u}}_{i} $为时均速度，$ {u}_{i}^{'} $为脉动速度，$ \overline{p} $为微元所受总压强，$ {\overline{F}}_{i} $为微元所受质量力，$ \mu $为动力黏度。

3）Standard k-ε方程

Standard k-ε方程是基于湍流的动力学原理和经验参数建立的，k代表湍动能，描述了液体流动的涡旋强度，ε代表湍动耗散率，表示了涡旋的衰减速率。Standard k-ε模型通过考虑方程组的湍流能量平衡和耗散率的传输方程，模拟了湍流运动和相互作用。

$\begin{split} \frac{\partial }{\partial t}\left(\rho k\right)+&\frac{\partial }{\partial {x}_{i}}\left(\rho k{u}_{i}\right)=\frac{\partial }{\partial {x}_{i}}\left[\left(\mu +\frac{{\mu }_{t}}{{\sigma }_{k}}\right)\frac{\partial k}{\partial {x}_{j}}\right]+\\ &{G}_{k}+{G}_{b}+\rho \epsilon -{Y}_{M}+{S}_{k} ，\end{split}$

(3)

$\begin{split} \frac{\partial }{\partial t}\left(\rho \varepsilon \right)+&\frac{\partial }{\partial {x}_{i}}\left(\rho \varepsilon {u}_{i}\right)=\frac{\partial }{\partial {x}_{j}}\left[\left(\mu +\frac{{\mu }_{t}}{{\sigma }_{\varepsilon }}\right)\frac{\partial \varepsilon }{\partial {x}_{j}}\right]+\\ &{C}_{1\varepsilon }\frac{\varepsilon }{k}\left({G}_{k}+{C}_{3\varepsilon }{G}_{b}\right)-{C}_{2\varepsilon }\rho \frac{{\varepsilon }^{2}}{k}+{S}_{\varepsilon }。\end{split} $

(4)

式中：$ {G}_{k} $、$ {G}_{b} $为平均速度梯度和浮力引起的湍流能k的派生项；$ {C}_{1\varepsilon } $、$ {C}_{2\varepsilon } $、$ {C}_{3\varepsilon } $为经验常数；$ {\mu }_{t} $为湍流黏度。

1.2 自由液面捕捉技术

VOF法基于体积分数概念，对液体和气体在计算域内的体积进行守恒，考虑了多相流的连续方程和Navier-Stokes方程，适用与液体和气体之间的复杂界面问题^[16] ，其控制方程为：

$ \begin{array}{c}\displaystyle\frac{\partial }{\partial t}\left({F}_{i}\right)+\displaystyle\frac{\partial }{\partial x}\left(u{F}_{i}\right)+\displaystyle\frac{\partial }{\partial y}\left(v{F}_{i}\right)=0，\\ {\displaystyle\sum }_{i=1}^{N}{F}_{i}=1。\end{array} $

(5)

式中：$ {F}_{i} $为体积分数函数，$ {F}_{i} $=1表明单元格全为液体，$ {F}_{i} $=0表明单元格全为空气，0＜$ {F}_{i} $＜1，表明单元格有空气和液体；N表示多相流相的数量。

2 数值模型验证

为验证数值模型的可靠性，本文选择文献[13]中的实验模型，该实验模型长1 m，宽0.1 m，高1 m，装载高度0.125 m。横荡激励为：

$ {x=A{\mathrm{sin}}(2{\text π} ft)}。$

(6)

式中：A为激励幅值，A=0.1 m；f为激励频率，f=0.4 Hz。

图1为该液舱晃荡过程自由液面演化的试验结果与数值模拟结果。由图1(a)可知，液舱内出现行波并向右舱壁面传播，波峰前缘波形也在不断变化，直到波峰前缘撞击到右侧舱壁，此时液体的速度急剧降低，一分部动能转化为压力能，另一部分动能使得液体沿着右侧壁面向上爬升，转变为势能，随后液体在重力作用下回落，势能又转化为动能，并形成行波向左舱壁面传播。由图1可知，晃荡过程自由液面演化的数值解与试验结果^[13]高度一致，这表明本文采用的数值模拟方法研究晃荡现象可信。

3 网格独立性验证

液舱晃荡涉及流体复杂运动和非线性响应，通过验证网格独立性来评估数值模拟结果对网格的依赖，以增加结果的可靠性。当模拟结果相对不受网格数量干扰，则可确定最优网格尺寸。为验证网格相关性，本模型设计8套网格，数量分别为1万、2.25万、3.24万、4万、4.84万、6.25万、7.84万和9万。外激励为横荡，激励规律如式(6)所示；液舱固有频率参考式(7)，得出$ {{f}}_{0} $=0.76 Hz；中等、大幅激励参考文献[17]，取A/L=0.1，因此激励幅值为0.1 m；装载率为30%；监测压力点为左舱壁自由液面处，晃荡压力取每个激励周期的压力峰值的算术平均值；时间步长取0.001 s。

矩形液舱固有频率公式：

$ {{f}}_{{n}}=\frac{1}{2}\sqrt{{g}{n}\mathrm{tan}{h}\left(\frac{{{\text{π}} }{h}}{L}\right)/{{\text{π}}}{L}}。$

(7)

式中：n为阶数；g为重力加速度；h为载液高度；L为自由液面宽度；$ {f}_{n} $为液体的固有频率。

由图2可知，网格数量为1万的模型的晃荡压力为4580 Pa，数量为4万的晃荡压力为4809 Pa，两者相差4.7%，网格为6.25万的晃荡压力为4839 Pa，与4万网格模型晃荡压力相差0.7%，网格数量为9万的模型的晃荡压力为4815 Pa，与4万网格模型的晃荡压力相差0.1%，表明网格数量在4万以上的模型对晃荡结果依赖程度很小，且时间步长0.001 s也能较好反应液体的运动情况。因此，本文后续的数值模拟计算将选定网格数量为4万和时间步长为0.001 s。

4 单自由度大幅激励下的液舱晃荡特性

矩形液舱晃荡最剧烈是在液体最低阶固有频率处，其理论计算公式见式(7)。研究表明^{[18 − 19]}液舱中液体的最低阶固有频率的理论值、数值模拟值、实验值会有差别，影响因素包括液体粘性、晃荡的非线性程度等。由于液体固有频率理论值是基于无旋、无粘及小幅晃荡条件得出，因此对图3模型进行晃荡数值模拟计算，通过扫频的方式确定舱内液体的最低阶固有频率，具体工况见表1，外激励公式见式(6)，压力监测点布置在液舱左舱壁自由液面处。

由图4(a)、图4(b)可知，晃荡压力随着激励频率由小逐渐增大而接近固有频率时，呈快速增大趋势，在达到最大值后，随着激励频率继续增大，晃荡压力会缓慢减小。由图4(c)可知，晃荡压力随着激励频率由小逐渐增大而接近固有频率时，呈缓慢增大趋势，并在达到最大值后，随着激励频率继续增大，晃荡压力会快速减小。由图4(d)、图4(e)可知，晃荡压力随着激励频率由小逐渐增大而接近固有频率时，呈快速增大趋势，并在达到最大值后，随着激励频率继续增大，晃荡压力会快速减小。

由图4(f)可知，装载率为15%、25%、50%、65%及75%时液体一阶共振频率分别为0.58、0.62、0.99、0.97及0.94 Hz。对于单自由度的大幅激励液舱晃荡，在低装载下（15%），液体的一阶固有频率与理论计算值一致，25%装载下，液体的一阶固有频率小于理论计算频率，而中、高装载，液体的一阶固有频率大于理论计算频率。对于基于线性势流理论得到的液体一阶共振频率随着装载率的增加而增加，但是对于大幅激励而言，液体一阶共振频率随着装载率的增加则表现为先增加后略有减小。

图5表明液体在共振频率下进行大幅激励，自由面运动剧烈，出现冲顶、飞溅、翻卷和破碎等非线性行为，且大量空气被卷入晃荡液体中，连续的大气泡在液体的剪切作用下被粉碎为小气泡，甚至可能对壁面产生射流冲击。

图6为25%装载率下晃荡压力时程图，结合图4(b)可知，当液舱受到的激励频率小于等于$\overline f_1 $(0.6 Hz)和大于等于$\overline f_2 $(0.74 Hz)时，每个激励周期内压力的峰值大致相等，呈现规律性的压力变化，当液舱受到的激励频率在$\overline f_1 $和$\overline f_2 $之间时，每个激励周期内的压力峰值呈非线性变化，尽管激励频率在[$\overline f_1 $, $\overline f_2 $]范围内或其范围外，晃荡压力都表现为明显的冲击特点，但是其峰值压力幅值规律相差较大，激励频率在($\overline f_1 $, $\overline f_2 $)范围内冲击压力幅值可能满足某种统计规律。以上表明，即使液舱受到大幅激励，但在受到不同的激励频率时也会表现出不同的晃荡压力特性，尤其在共振状态，晃荡压力的冲击幅值表现出极强的非线性特性。

大幅激励下的晃荡液体共振频率及晃荡压力不能通过基于线性势流理论的相关公式得到。小幅激励、中等幅度激励及大幅激励的界定由液舱尺寸决定^[18]。对于本文所研究的如图3所示的液舱模型，其大幅横荡激励的具体参数如表2所示，外激励频率取液舱内液体一阶共振频率数值模拟值，横荡位移公式见式(6)。

由图7可知，25%装载率的液舱在受到小于或等于0.08 m激励幅值时，每个激励周期的压力峰值大致相等。而当液舱受到大于或等于0.09 m激励幅值时，每个激励周期内的晃荡压力峰值变化很大，且晃荡压力出现快速增大特性，这表明0.09 m激励幅值是液舱中的液体进入大幅度晃荡的临界值，砰击载荷的非线性表现的非常强烈，此时，长期预报时的冲击特性比小、中等激励情况发生根本性的改变。

由图8可知，液舱自由液面处的晃荡冲击压力随着激励幅值的增大呈现3种特征，在小幅激励下，晃荡压力是线性增大，随着激励幅值增大到中等激励幅值后，晃荡压力呈快速增大特性；继续增大激励幅值，则晃荡压力在一定范围内呈饱和状态。液舱在15%装载下，晃荡压力最大值为37.8 kPa，在75%装载下，晃荡压力最大值为13.1 kPa，晃荡压力相差188.5%，表明液舱装载率与晃荡压力负相关，即随着装载率提高，最大晃荡压力会减小。

5 大幅横荡-垂荡联合激励下液舱晃荡特性

当激励幅值在0.05～0.08 m时，液体将出现自由液面大幅度晃荡现象，后续研究大幅度晃荡现象的激励幅值取大于0.08 m。为研究大幅度横荡-垂荡作用下的晃荡特性，设计表3工况进行数值模拟，其中液舱装载率为15%、25%、50%、65%、75%，激励频率分别取0.58、0.62、0.99、0.97、0.94 Hz。

图9中，晃荡压力取值是响应进入稳定阶段后的20～350 s内每个激励周期的压力峰值的算术平均值，压力测点位置是液舱自由液面处左舱壁。可知，在大幅激励下，液舱受横荡-垂荡激励和只受横荡激励的晃荡压力总体变化趋势相似，但是，25%和50%装载情况下对于0.08 m的横荡激励，如果考虑垂荡，则复合激励下的晃荡压力明显大于单自由度横荡激励情况。液舱在低装载工况下，受大幅度横荡-垂荡激励的液舱晃荡压力最大值和受大幅度横荡激励的液舱晃荡压力最大值大致相同，而随着装载率达到中、高装载率，受大幅度横荡-垂荡激励的液舱的晃荡压力最大值会明显高于受大幅度横荡激励的液舱。此外，液舱在低装载工况下，随着激励幅值变化，受横荡-垂荡作用的液舱和受横荡作用的液舱的晃荡压力值大小交替变化，而随着装载率提高到中、高装载率，受大幅度横荡-垂荡激励的液舱的晃荡压力总大于受大幅度横荡激励的液舱。

由图10可知，Case4和Case5压力周期吻合，表明在共振频率和大幅度横荡激励下引入大幅度垂荡激励并不会改变晃荡压力周期，液舱在低装载工况下，Case4的晃荡压力大于Case5的晃荡压力，而在中、高装载下Case5的晃荡压力大于Case4的晃荡压力，对于压力次峰，Case4和Case5的次峰大致吻合，表明引入的大幅度垂荡激励不会改变次峰值。

6 结　语

1）基于线性势流理论得到液体一阶共振频率随装载率的增加而增加，但对大幅激励而言，液体一阶共振频率随着装载率的增加则表现为先增加后略有减小。

2）随着激励幅值由小增大，晃荡压力-激励幅值曲线呈现出由线性增大到非线性快速增大，再下降到在一定范围内所呈现的“饱和特性”；大幅度横荡激励下，低装载液舱晃荡冲击压力大于高装载率液舱。

3）对比受大幅横荡-垂荡耦合激励的液舱和受大幅横荡激励的液舱，其压力特性在低装载下相似，中、高装载下，受大幅横荡-垂荡激励的液舱晃荡压力更高，而引入大幅度垂荡激励不会改变晃荡压力周期和压力次峰。