随机激励下液舱晃荡数值模拟分析

姜胜超，徐博，王子豪

1 大连理工大学 船舶工程学院，辽宁 大连 116024

2 中国船舶及海洋工程设计研究院，上海 200011

0 引 言

晃荡是一种非常常见的流体运动现象，通常发生在装有部分液体的液舱中，如航行中的液货船。当外界激励的频率接近液舱内液体振荡的固有频率，或是激励的幅值非常大时，很容易引起液舱内液体的剧烈晃荡，进而对舱壁或是舱顶产生强烈的冲击压力，从而造成结构破坏。晃动也具有重要的小尺度效应，如冲击压力、由气液混合物及因流体黏度和旋涡流动引起的机械能耗散。另外，液体晃荡与船舶运动之间还可发生强耦合。因此，了解液舱晃荡的基本物理原理对海上作业和结构安全具有重要意义。

在理论研究方面，Faltinsen[1]提出了一种基于势流理论的三阶晃荡模型的非线性分析方法，考虑了液舱大幅度运动的非线性效应。Faltinsen 和Timokha[2]完善了求解液体在有限液深矩形容器中非线性晃荡的半解析方法。Zhang 等[3]基于势流理论和摄动展开研究了三维液舱在纵荡和横荡耦合激励下的二阶晃荡问题，推导了理论解析解，发现当任意2 个激励分量的和频或是差频等于其中一个固有频率时，就会发生二阶共振。高速计算机的发展和计算流体力学（CFD）技术的成熟为分析晃荡运动提供了一种全新且高效的方法-数值模拟。Kim[4]建立了一个用于预测各种二、三维液舱冲击发生情况的有限差分模型，并在液舱顶部附近采用一个特殊的缓冲区进行冲击模拟，同时使用一个单值函数对晃动界面进行跟踪。随后，Kim 等[5]对该模型进行改进，将缓冲区扩展到罐顶附近的倾斜边界，减轻了以往对网格分辨率和时间增量的敏感性。Biswal 等[6]采用有限元方法对带有刚性隔板的二维矩形液舱中液体的非线性晃动响应进行了研究。Lu 等[7]针对有/无隔板矩形液舱内液体的晃动问题，建立了基于非惯性参考系的黏性流体模型，数值试验结果表明，耗散效应对带隔板液舱和无隔板液舱的晃动响应均有显著影响。由于物理耗散的作用，与固有频率相关的晃荡响应分量最终可以被阻尼耗散掉，但它们却完全保留在了势流解析解中，这就解释了势流理论对晃动振幅的估计普遍过高的问题。Liu 和Lin[8]基于两相流三维数值模型，研究了液舱晃荡过程中的黏性阻尼作用。肖凯隆等[9]基于多液舱晃荡和船舶运动耦合黏性流数值模型，对船舶横摇激励下的液舱晃荡特性进行了研究，并提出了多液舱的二维模型简化方法。

以上研究大多是针对规则激励，鲜有针对随机激励下液舱晃荡问题的研究，但相对于规则激励，随机激励更接近于海上波浪等对液舱作用的不规则性。本文拟在介绍使用黏性数值模型和随机激励原理的基础上，分析随机激励下液舱晃荡的瞬态效应和时间敏感性，以及随机激励下谱峰频率和幅值对液舱晃荡的影响。

1 基本原理与数值方法

本文将主要基于STAR-CCM+软件，开展强迫激励下液舱晃荡的数值模拟。考虑到不可压缩黏性两相流问题，因涉及网格运动，所以使用任意拉格朗日-欧拉（ALE）观点的Navier-Stokes （N-S）方程为：

式中：ρ 为流体密度；xi，xj分别为i，j方向的位移分量；ui，uj分别为i，j方向的速度分量；uim为网格运动速度，指在i方向的速度分量；t为时间；µ为流体动力学黏性系数；fi为单位体积流体所受到的体积力，本文中仅为重力，取为9.81 m/s2。

为使模型能够满足波浪破碎时的模拟要求，采用流体体积（VOF）法捕捉自由水面运动。定义流体相函数φ为

其满足ALE 观点下的边界面方程：

进而可确定两相流的密度及动力黏性系数分布。

式中：ρw，ρa分别为水和空气中的密度；µw，µa分别为水和空气中的动力黏性系数。在进行数据处理时，取φ= 0.5 的等值线作为液体的自由水面。

对两相流控制方程（式（1）～ 式（2））及边界面方程（式（4））采用有限体积法（FVM）进行离散，其中，时间离散采用欧拉格式，散度和梯度计算分别使用Gauss Vanleer 和Gauss linear 格式，扩散项使用Gauss linear corrected 格式。N-S 方程的求解 采 用PISO （ pressure implicit with splitting of operators）方法[10]，其中，速度方程可直接采用代数方法求解，压力方程则采用BI-CGSTAB[11]方法迭代求解。计算时，取容器顶盖中的点为参考压力点，设为0 Pa。

数值计算中，时间步长的选取需满足CFL（Courant-Friedrichs-Lewy）条件（式（6）），可通过克朗数Cr进行自动选取。

式中，umax，vmax，wmax分别为x，y，z方向的最大速度分量。计算中，克朗数Cr原则上取为1.0 即可满足CFL 条件要求，本文取Cr=0.20，用以保证自由水面做大振幅非线性运动时的计算精度。

2 数值模型的验证

采用上述数值模型对水平激励下矩形液舱内的流体晃荡问题进行研究。数值模型设置25 000 个网格，其中靠近自由液面和舱壁的区域设置为加密区，数值模型网格最小尺寸为5 mm，最大尺寸为10 mm。为了验证数值模型的正确性，对宽度W= 1 m、装载水深d= 0.5 m 的矩形液舱晃荡问题进行数值模拟，此时，液舱晃荡自振频率ω0= 5.316 rad/s，如图1 所示，液舱受到正弦激励作用：

图1 矩形液舱几何模型示意图Fig. 1 Schematic diagram of rectangular tank geometry model

式中：A为激励幅值；ω 为激励频率。

Faitinsen[12]的解析解以及Liu 和Lin[8]的数值结果均是在t= 0 时突然加速，并未使用缓冲函数。由于本节的目的是对本文数值模型进行验证，因此在数值模拟中也采用了同样的突然加速方法。选取工况1（A= 0.01 m，ω = ω0）和工况2（A= 0.000 4 m，ω = 0.95 ω0）这2 种激励工况进行计算，并与Faitinsen 的解析解及Liu 和Lin 数值模拟结果进行了对比，结果如图2 所示（图中，η 为波面高度）。从中可以看出，3 种方法均吻合较好。

图2 工况1 和工况2 下右舱壁波高历程线比较图Fig. 2 Comparison of the right bulkhead wave height at Case 1 and Case 2

进一步对液舱宽度W= 0.570 m、高度H=0.30 m、装载水深d= 0.150 m 的情况进行研究，此时，液舱晃荡自振频率ω0= 6.058 rad/s。仍考虑液舱做正弦运动的情况，选取工况3（A= 0.005 m，ω = 0.583ω0）和工况4（A= 0.005，ω = ω0）这2 种情况进行研究，数值模拟结果仍与Faitinsen[12]的解析解及Liu 和Lin[8]的实验与数值结果进行了对比，结果如图3 所示。从图3（a）中可以看出，对于非共振激励情况，解析方法、实验方法及数值结果与本文方法均吻合较好。对于共振激励情况，如图3（b）所示，可以看出在初始阶段采用3 种方法所得结果几乎重合，但随着时间的推移，当t＞2 s 时，共振激励条件下的实验结果出现了明显的波峰变高、波谷变坦的非线性自由水面特征，采用Liu 和Lin[8]的方法及本文数值方法可以对其进行准确描述，但解析解却不行，说明采用非线性液舱晃荡模型是有必要的。由上述对比可以发现，本文的数值结果与实验结果吻合较好，说明了本数值方法的正确性，能够进行下一步的数值模拟。

图3 工况3 和工况4 下右舱壁波高历程线比较图Fig. 3 Comparison of the right bulkhead wave height at Case 3 and Case 4

3 随机激励的位移谱

仿照JONSWAP 波浪谱给出水平激励的位移谱，即

其中：

式中：γ 为谱峰升高因子； σa， σb分别为谱峰左、右两侧的峰形参量；α 为Philips 参数；U为水面10 m高度处风速；k为常数，与谱峰升高因子对应；Hs为设计（有义）波高；X为风距，km；g为重力加速度，取9.81 m/s2；ωp为谱峰频率。

自由液面波高历时线可以由大量的正弦波线性叠加得到：

式中：ωi为i个线性正弦波的圆频率；Nω为所有线性正弦波的数量；Ai，φi分别为每个线性波的波幅和相位（φi为随机变量，均匀分布在0～2π 区间内）。其中，Ai可通过下式得到：

式中，Δω 为频率间隔。由于高频对产生的波没有贡献，故在一定频率处截断。

与上述随机波浪的产生类似，将下式作为液舱的水平激励，可以表示为

式中：x(t)为液舱受到的随机水平激励；Nω为所有线性余弦波的数量，本文取Nω=200，这足以描述液舱的随机激励位移。本节将以式（8）中的变量Hs和ωp作为参数来控制所产生随机水平激励的激励幅值和谱峰频率。

4 随机激励下的瞬态效应和时间敏感性

对于简谐激励，由于在初始时刻突然加入了激励速度，会导致瞬态效应的产生，故在固有频率处会产生晃荡响应分量，从能量的角度分析，由于实际的物理耗散，固有频率处的能量最终会被黏性阻尼耗散掉。因此，为了避免瞬态效应对数值模拟结果的影响，可以对激励位移历程线做缓冲函数的处理，使V(t)=v(t)·Rm（v(t) 为激励速度，V(t)为处理后的激励速度，Rm为缓冲函数），从而使舱内液体开始缓慢地晃动。

式中，Tm为缓冲时间（作用时间100 s），一般为周期的整数倍。

下面，将讨论瞬态效应对随机激励的影响。选取工况1（ωp=0.1ω0，Hs=0.01d）和工况2（ωp=0.65ω0，Hs=0.01d）这2 组具有代表性的工况进行分析。图4 所示为在工况1 的激励条件下有/无缓冲函数处理的右舱壁波高历时线对比图。由图中可以看出，经过缓冲处理的右舱壁波高缓慢增大，缓冲结束后2 条曲线的差别较大，主要表现为缓冲处理后波的周期明显变大，说明此时液舱的晃荡频率发生了变化。进一步采用快速傅里叶变换（FFT）算法进行处理，得到右舱壁波高在100～500 s时间段的能量谱，如图5 所示。由图5（a）可以看到，对于未经缓冲处理的液舱，在激励位移谱峰频率所在的0～2 rad/s 范围内，观察到有一个宽峰，但该峰值的量级较小，不过在一阶固有频率处观察到有一个明显的峰值。对比图5（b）发现，经过缓冲处理后的液舱仅在位移谱峰频率所在的0～2 rad/s 范围内观察到有一个宽峰，这和图5（a）类似，但在其他各阶固有频率处却未出现峰值。这说明图5（a）中一阶固有频率处出现的峰值是由瞬态效应导致的。

图4 工况1 下右舱壁波高历时线Fig. 4 Right bulkhead wave height history at Case 1

图5 工况1 下右舱壁波高在100～500 s 时间段的FFT 谱分析Fig. 5 Fast Fourier transform power spectrum analysis of right bulkhead wave height during 100～500 s at Case 1

图6 所示为在工况2 的激励条件下有/无缓冲函数处理的右舱壁波高历时线对比图。进一步采用FFT 算法进行处理，得到右舱壁波高在100～500 s时间段的能量谱如图7 所示。从图7（a）中可以看到，对于未经缓冲处理的液舱，能量谱中出现了多个峰值，主要对应于液舱的各阶固有频率。值得注意的是，一阶固有频率处的峰值远大于其他各处，在能量谱中占据主导地位。经对比图7（b）发现，经过缓冲处理后的液舱其各阶固有频率处的峰值均大幅降低，说明图7（a）中过高的峰值是由瞬态效应所导致，致使固有频率处对应的晃荡响应分量偏大。

图6 工况2 下右舱壁波高历时线Fig. 6 Right bulkhead wave height history of Case 2

图7 工况2 右舱壁波高100～500 s 时间段FFT 谱分析Fig. 7 Power spectrum analysis of right bulkhead wave height during 100～500 s at Case 2

通过对比图8(a) ～图8 (d)中的数值模拟结果可以发现，经过缓冲处理后，液舱在缓冲结束后外部激励频率 ωe处的峰值随着时间的推移其能量级的大小保持相同，说明达到了稳定阶段，而未经缓冲处理的液舱则需要经过更长时间的物理耗散才能消除瞬态效应的影响，最终只保留了激励频率 ωe处的晃荡响应分量。

图8 不同时间段的能量谱对比Fig. 8 Energy spectrum comparison at different time periods

通过以上分析可以发现，瞬态效应对随机激励的液舱晃荡影响较大，尤其是对于激励频率远离固有频率的工况，瞬态效应会在各阶固有频率处激起较大的峰值，而由于一些峰值的量级与未经过缓冲处理的相比差异较大，会影响进一步的分析，因此对于随机激励的模拟，对激励位移历时线进行缓冲函数处理十分有必要。

由于随机激励的随机性，随着时间的推移，可能会对液舱晃荡的能量和运动造成影响，因此同简谐激励相比，需要进行更长时间的数值模拟才能得到稳定、可靠的分析结果。但过长的模拟时间也意味着计算资源的浪费，因此有必要分析随机激励下的时间敏感性以确定合理的模拟时间。不考虑缓冲处理的前100 s，将1 300 s 的数值模拟结果分为100～500 s，500～900 s，900～1 300 s 三段，分析时间对随机激励的影响。对3 组工况下不同时间段右侧舱壁波高历时曲线概率密度图进行对比，结果如图9～图11 所示。由图可以看出，3 个时间段的概率密度图基本符合高斯分布曲线，高斯分布位置函数为ζ(t) = 0，即静止波面处。而且通过对比，还发现概率密度图的区间和概率分布近似一致，这说明不同时间段的随机激励模拟结果具有相同的分布规律。

图9 右舱壁波高概率密度分布图（ω p=0.65ω0 ， Hs=0.01d ）Fig. 9 Probability density distribution of right bulkhead wave height (ω p=0.65ω0 ， Hs=0.01d )

图10 右舱壁波高概率密度分布图（ω p=ω0 ， Hs=0.01d ）Fig. 10 Probability density distribution of right bulkhead wave height (ω p=ω0 ， Hs=0.01d )

图11 右舱壁波高概率密度分布图（ω p=1.5ω0 ， Hs=0.01d ）Fig. 11 Probability density distribution of right bulkhead wave height (ω p=1.5ω0 ， Hs=0.01d )

又由表1 所示可以看出，在3 种有代表性的频率下，3 个时间段的平均幅值Aave和对应的周期T均十分接近，因此，模拟500 s 时的结果进行分析是合理、可靠的，后文的分析都将在该时间段内。

表1 振幅周期统计分析Table 1 Statistical analysis of amplitude and period

5 随机激励谱峰频率分析

图12 所示为不同谱峰频率下的右舱壁波高历时线，激励幅值均设置为Hs=0.01d。从图中可以看出，随机激励的随机性波高历时线呈现出明显的非线性，这与简谐激励有所不同。因此，将进一步对波高历时线进行快速傅里叶变换，以分析随机激励下不同谱峰频率的液舱晃荡。

图12 右舱壁波高历时线 (Hs = 0.01d)Fig. 12 Right bulkhead wave height history (Hs = 0.01d)

采用图1 所示液舱模型，首先考虑3 种工况，即ωp分别 为0.1ω0(0.53 rad/s)，0.35ω0(1.86 rad/s)和ω0(5.31 rad/s)，激励幅值统一设置为Hs= 0.01d。右舱壁波高的能量谱如图13 所示。从图13（a）中可以看出，在0～2 rad/s 频率范围内出现了一个宽峰，但该峰值的量级较小，结合激励位移能量谱可知，该峰值正好位于谱峰频率附近，处在主能量范围内。值得注意的是，在一阶固有频率处未出现峰值，这和以往关于随机激励的结论不同，可能的原因是以往关于随机激励的研究没有排除瞬态效应的影响，导致对固有频率处的晃荡响应估计过高。从图13（b）中可以看出，在一阶固有频率处有一个较大的峰值，在二阶和三阶固有频率ω1，ω2处也可以看到2 个较小的峰值。这是因为上述工况的激励谱峰频率ωp比ω1和ω2更接近于ω0，所以ω0处的峰值对能量贡献最大。通过对比发现，随着ωp的增大，激励频率处的峰值消失，在各阶固有频率处出现了较大的峰值，且一阶固有频率处的峰值占据主导地位；同时，在其他频率处也出现了明显的峰值，如在一阶固有频率的1.5 倍和2 倍处，当ωp逐渐接近一阶固有频率时，各处的峰值均明显增大，舱内液体晃荡得越来越剧烈。最终，当ωp= ω0时，会发生共振现象。

为进一步研究更高的谱峰频率ωp对响应能量谱的影响，选取1.35ω0(7.17 rad/s)，1.65ω0(8.76 rad/s)和ω1(7.83 rad/s)这3 个谱峰频率模拟了随机激励下的舱内液体晃动问题。右舱壁波高的能量谱如图14 所示。从图中可以看出，二阶固有频率处的峰值是逐渐增大的，当峰值频率增大至接近二阶固有频率时，该处的峰值在能量谱中占主导地位。与此相反的是，一阶频率处的峰值是随着谱峰频率ωp的增大而逐渐减小的。据此可以预测，ωp接近哪一阶的固有频率，该频率处的峰值就将占主导地位。在图14 (a)中，峰值频率处于ω0和ω1之间，但更接近于ω0，所以在2 个固有频率处有2 个大的峰值，其中ω0处的峰值更大一点。从图中还可以看出，峰值也会出现在一阶固有频率的1.5 倍或2 倍处。通过对比图13（c）和图14（c）可以发现，当ωp= ω0时响应能量谱峰值的量级更大，因此同ωp= ω1相比，ωp= ω0时舱内液体的晃荡将更加剧烈和危险。然而，当ωp= ω1或是更高谱峰频率激励时，可能会引起舱内液体的高频振动。

图13 右舱壁波高能量谱（ω p ≤ω0 , Hs=0.01d）Fig. 13 Power spectrum of right bulkhead wave height (ω p ≤ω0 , Hs=0.01d)

图14 右舱壁波高能量谱（ω p ＞ω0 , Hs=0.01d）Fig. 14 Power spectrum of right bulkhead wave height (ω p ＞ω0 , Hs=0.01d)

6 随机激励幅值分析

图15 所示为固定谱峰频率ωp= ω0，激励幅值分别为Hs=0.002d，0.006d，0.01d时得到的3 种激励幅值下的右舱壁波高历时线对比。从图中可以看到，随着激励幅值的逐渐增大，出现了波峰变尖、波谷变坦的非线性现象，这和简谐激励时增大激励幅值观察到的现象相同。值得注意的是，在谱峰频率ωp= ω0下对液舱持续激励，由于输入的位移是随机分布的，故输出的波高没有出现持续增大的共振现象。

图15 右舱壁波高历时线 (ωp = ω0)Fig. 15 Right bulkhead wave height history (ωp = ω0)

在随机晃动问题的数值模拟中，由于波和力随时间的变化非常复杂，因此要得到非常准确的概率密度和能量谱结果很困难。但若随机激励的模拟时间合理，仍可以观察到一些普遍的现象，并得到有意义的结果。模拟中存在的非线性效应，甚至是在共振频率上存在的随机相位，都表明随着有义激励幅值的增加，输出的结果在分布上很可能会表现出非平稳趋势。图16 所示为Hs= 0.01d时右舱壁波高的概率密度分布图。从图中可以看出，对应的波高历程线上有一个较小的波谷极小值η = -18Hs和较大的波峰极大值η = 30Hs，因此它不是关于静止水面上下对称的，从图中可以看到，它是非高斯分布的。结合上面的分析可知，输出响应在本质上并非与输入激励线性对应，这种偏离高斯行为的情况是合理的。由于与随机平稳过程存在一些大致相似的特征，因此可以对偏差进行量化处理。

图16 右舱壁波高概率密度分布图（H s=0.01d）Fig. 16 Probability density distribution of right bulkhead wave height ( H s=0.01d)

根据前面的分析，非线性随机晃荡响应曲线具有更陡的波峰和更平坦的波谷，这就决定了其分布不同于线性随机波，因其偏离了高斯（正态）分布。对于线性波，波高或水动力的偏差应该为0。为了说明晃动波或水动力的非线性是如何影响所述偏差大小的，本文计算得到了3 个不同激励幅值Hs下的波高和水动力的标准差，以此分析非线性对偏差的影响。右舱壁波高和水平力的标准差如表2 所示。从中可以看出，随着Hs的增加，波高和力的偏移趋势随之增大，说明非线性越强，偏移越大。值得注意的是，这种增长似乎与Hs呈线性相关的趋势。进一步的分析表明，输出响应结果表现出了线性缓慢变化的非平稳过程特征。

表2 标准差统计分析Table 2 Statistical analysis of standard deviation

对于非线性响应，响应结果分布是不关于平均水线对称的，可以采用偏离度λ3来描述这种不对称。

式中：n为总的计算时间内的样本数量； σ为标准差；ηi为第i个时间点的波高；为所有波高的平均值。就线性波而言，其偏离度为0。这里，将以上面的3 个算例为例，通过计算Hs= 0.002d，0.006d，0.01d得到其对应的偏离度分别为0.063，0.115，0.283，这表明随着晃荡响应的非线性增强，波峰和波谷会越来越不对称。

根据矩形液舱的固有频率计算公式，可以计算得到液舱的一、二阶固有频率分别为5.31 和7.83 rad/s。采用FFT 算法对上述3 个算例的右侧舱壁波高进行处理，得到对应的能量谱如图17 所示。从中可以明显地看到，在一阶固有频率附近出现了一个非常明显的窄峰，此外还观察到有一个小得多的峰值出现在二阶固有频率处，其幅度至少低了一个数量级，在三阶固有频率附近还出现了一个更小的峰值。严格来说，在第i（i= 4，5，···）阶频率上也应该有相应的峰值。这些后续的峰值在目前的图中无法观察到，因为它们不在激励谱的主要能量范围内，其峰值的量级非常小，比一阶固有频率小了几个数量级。值得注意的是，在一阶和二阶固有频率范围内同样也出现了一些峰值，如1.5ω0处。此外从图中还可以观察到，在一阶固有频率的2 倍处（10.62 rad/s）有一个明显的峰值，这可能是因为Hs的增加导致了非线性的增强。由上述结果可以推断，能量主要集中在一阶固有频率附近的窄带内。对比图17（a）～图17（c）可以看出，随着激励幅值Hs的增加，对应峰的峰值也会增大。

图17 右舱壁波高能量谱Fig. 17 Power spectrum of right bulkhead wave height

7 结 论

本文以无隔板液舱为研究对象，首先对所采用液舱的数值模型以及随机激励数值模拟的前处理、数值计算和后处理分析过程予以介绍，接着分析瞬态效应对随机激励的影响，并进一步进行时间敏感性分析，研究了时间对随机激励的影响，最后在此基础上研究了随机激励幅值和随机激励谱峰频率对随机响应的影响，主要得到如下结论：

1） 瞬态效应对随机激励的液舱晃荡影响较大，尤其是对于激励频率远离固有频率的工况，瞬态效应会在各阶固有频率处激起较大的峰值，其中一些峰值的量级和未经过缓冲处理的相比差异较大，这会对进一步的分析产生影响。因此，对于随机激励的模拟，尤其是对激励频率远离固有频率的工况，对激励历时线进行缓冲函数处理非常有必要。

2） 随机激励谱峰频率的研究表明，当激励谱峰频率从低频向高频移动时，在远离固有频率工况下，在激励谱峰频率范围内会出现一个峰值，而非出现在固有频率处；随着谱峰频率的增大，激励频率处的峰值会消失，而在各阶固有频率处则会出现较大的峰值，并且一阶固有频率处的峰值将占据主导地位；当谱峰频率远离一阶固有频率向更高频移动，在接近第i阶固有频率时，该频率处的峰值将占据主导地位。

3） 随机激励幅值的研究表明，随着激励幅值的逐渐增大，会出现波峰变尖、波谷变坦的非线性现象，较强的非线性波将导致概率密度相对于平均水线更加不对称，这表明随着波非线性的增强，概率密度将逐渐偏离正态分布；随着激励幅值的逐渐增大，还将导致能量谱中的峰变得更大，更多的峰出现在固有频率处，甚至是固有频率的2 倍处。因而非线性的强弱可以用标准差和偏离度来定量描述。