0 引　言

晃荡现象是液舱内液体在外界激励下产生的复杂流动行为^[1]，具有显著的非线性特征，难以通过简单线性模型准确预测。当外界激励频率接近液体的固有谐振频率时，液体进入不稳定状态，容易引发剧烈晃荡，对液舱舱壁产生局部冲击载荷，影响结构安全。此外，由于液舱晃荡过程中存在非线性共振行为，在不同波浪激励下，液体的动能将通过多种非线性机制传递，从而进一步增加了晃荡过程的控制和预测难度。近年来，由于晃荡问题与船舶设计及航行安全的直接关系，已成为流体力学与船舶海洋工程领域的研究热点之一，同时受到了工程界的广泛关注^{[2 − 3]}。

与贮液容器在水平激励下产生的横波晃荡不同，施加纵向激励时，系统会引发参数性晃动。当外界激励符合特定周期性函数时，液体自由表面会出现共振现象，此时所激发的液面波称为法拉第波^[4]。基于马修方程^[5]的稳定性求解问题，研究表明在固定激励系数下纵向激励的频率对应即次谐波、谐波和超谐波3种不稳定模式。这些模式下液体自由表面易失稳晃动，表现出复杂的动力学特征，特别是在液舱重心发生变化的情况下。因此进一步探讨水平激励与纵向法拉第波晃荡的叠加效应，对于深入理解复杂环境下液舱的非线性动力响应具有重要意义。

液舱晃荡时，自由液面的响应是一种受重力影响的非线性时变波动现象，其动力学行为与外界激励频率密切相关^[6]。为了有效识别晃荡的动态变化及其内部的非线性关系，使用合适的时频域分析工具至关重要。快速傅里叶变换（Fast Fourier Transform, FFT）、小波变换（Wavelet Transform, WT）、小波二阶相干谱（Wavelet Bicoherence Spectrum, WBS）等工具，能够从时频角度^{[7 − 9]}揭示晃荡现象的复杂动力学特征，捕捉晃荡的瞬时频率成分，为研究晃荡的共振行为及非线性动力学提供了全新的分析思路。

响应表现上，于曰旻^[10]从波高幅值的角度开展了高载液率下矩形液舱晃荡的一阶横波共振研究，结果表明，高载液液舱共振时表现出液面破碎等非线性特征。王立时等^[11]通过施加纵向次谐波成功激发了二维容器的一至四阶振型，并从中提取了液体的固有频率。谢梦斌等^[12]基于FSRU进行了横荡和纵荡的数值研究，评估了晃荡所引起流场扰动，这一过程也体现了液体晃荡的非线性特征。张会霞等^[13]开展了横摇-横荡复合激励的晃荡压力特性研究，分析了低载液率下晃荡载荷的时域与频域特征。响应机理上，Liu等^[8]研究了密度分层液体自由液面在法拉第波下的响应机理，一定程度上解释了界面波的晃荡行为。薛米安等^[14]则从时频域分析的角度探索了密度分层液体的非线性晃荡机理。Zhang等^[15]进行了实验和数值仿真的研究，通过对波面信号进行快速傅里叶变化，分析了不同周期及幅度激励下的响应机理。刘玉龙等^[16]开展了不同频域入射波与液舱晃荡非线性作用机理的研究，解释了入射波频域与液舱运动幅值的关系。

综上所述，本文主要聚焦于矩形液舱横波及法拉第波共振时刻的研究场景，结合数值仿真的方法，针对液舱共振时的波面幅值与载荷响应特点，开展了基于小波变换和小波二阶相干谱等的时频域分析，以更好地解释晃荡过程中的动态响应与非线性过程。

1 理论分析与数值模型 1.1 理论分析 1.1.1 二维矩形液舱固有频率计算

液舱的理论固有频率基于Faltinsen^[17]方法得来，其计算公式为：

$ \left\{ \begin{gathered} {\omega _m}^2 = g{k_m}\tanh ({k_m}h)，\\ {k_m} = m{\text{π /}}L。\\ \end{gathered} \right. $

(1)

式中：$ g $为重力加速度；$ h $为液体高度；$ L $为液舱长度；$ {\omega _m} $和$ {{\text{k}}_m} $和分别为圆频率和波数；$ m $为晃荡模态的阶数。当$ m $=1时，对应的是液舱的一阶固有频率。

1.1.2 基于Matlab的Morlet小波变换与小波二阶相干谱

小波变换是一种时频分析方法^[18]，能够同时表征信号的时间和频率信息。相较于傅里叶变换，Morlet小波具有良好的频率分辨率，适用于分析非平稳或局部振荡类信号，能揭示晃荡中主频、谐波及非线性作用的演化过程。本文中小波变换基于Matlab中的cwt函数实现，函数依赖于有限离散采样点进行，并进行了归一化处理，因此有其离散化表达：

$ W[\psi {\text{(}}a,b{\text{)}}] =\frac{1}{{\sqrt {\left| a \right|} }}\sum\limits_{n = 1}^N x [n]{\psi ^*}\left(\frac{{n - b}}{a}\right)\Delta t。$

(2)

式中：$ W[\psi {\text{(}}a,b{\text{)}}] $为母小波函数$ \psi $在尺度$ a $和位置$ b $下的小波变换系数的结果；$ {\psi ^*} $为母小波函数的共轭复数；$ x[n] $为输入的$ N $个数据中第$ n $个离散数据；$ \mathrm{\Delta } $$ t $为信号的时间间隔。

小波二阶相干谱是一种结合小波变换与双谱分析的非线性信号处理方法，其核心思想是利用小波变换在时间和频率域的局部化特性，结合双谱分析的非线性相干性检测能力，提取信号中的非线性相互作用信息，从而有效揭示系统动态行为中的非线性耦合和相干现象，适用于晃荡信号中的非线性特征分析。其在Matlab中的离散化表达为：

$ {B({a_1},{a_2}) = \sum\nolimits_b W [\psi ]({a_1},b) \cdot W[\psi ]({a_2},b) \cdot W[\psi *]({a_1} + {a_2},b)。} $

(3)

式中：$ B({a_1},{a_2}) $为在尺度$ {a_1} $和$ {a_2} $下的二阶相关谱。

1.2 数值模型

二维矩形液舱的仿真计算基于Fluent 2024 R2进行，网格划分基于Mesh软件的四边形结构化网格。矩形液舱的边长建模设置为20 cm。本文中液体的载液率设置为60%，理由如下：一是实际生活中液舱通常会避免在低载液率下运行，因此过低的载液率无法代表一般工况；二是由于本文的工况设置围绕固有频率进行，在更高的载液率下液体时常出现“冲顶”现象，干扰波面数据的收集，因此本文选取60%载液率作为代表工况。

计算基于显示体积函数法（Volume of Fluid Method, VOF）的空气-水两相流数值模型，同时结合标准k-ε湍流模型及标准壁面函数，求解器使用压力隐式算子分裂法（Pressure Implicit with Splitting of Operators, PISO），压力插值格式选用压力交错格式（Pressure Staggering Option, PRESTO），动量及湍流计算均使用二阶迎风格式。为提高计算精度，计算过程中采用自适应时间步长法，总计算时间为20 s，全局库朗数设定为1.5。在数值计算过程中，液舱壁面则假定为刚性结构。

液舱动网格的设置基于Fluent自定义函数功能中的“CG_MOTION”函数，沿X/Y方向做往复运动。液舱激励可以表示为：

$ \left\{ \begin{gathered} x/{{y}} = A\sin (\omega t)，\\ {v_{x/y}} = \omega \cdot A\cos (\omega t)。\\ \end{gathered} \right. $

(4)

式中：$ A $为振幅，大小固定为0.5 cm；$ \omega $为运动周期的圆频率，rad/s。

横波激励的一阶解析固有频率$ {\omega _0} $由式(1)计算得到，并以此为频率中心，参照文献[10]中的频率取值，以0.03 Hz间隔扩展至7个横波工况；纵向激励设置根据马修方程有3种不稳定模式^[5]，其中超谐波激励为0.667 $ {\omega _0} $、谐波激励为$ {\omega _0} $、次谐波激励为2 $ {\omega _0} $，具体见表1。

2 数值模型验证

首先针对本文所使用的数值模型，开展网格无关性的验证。使用1.2中的二维矩形液舱，采用相同条件算的幅值（0.5 cm）和频率激励*（$ \omega $=12.1305 rad/s）进行分析。在保持全局库朗数（Courant Friedrichs Lewy condition, CFL）为1.5的条件下，提出如下3种网格：粗网格D₁的尺寸为5 mm，中等网格D₂为1 mm，细网格D₃为0.1 mm。3种网格的计算结果如图1所示，可知，除最初开始的1 s内，3种网格的波面幅值稍有差异，其余时刻下基于D₂和D₃网格得到的结果非常接近，即使采用粗网格D₁，其计算精度与细网格D₃的相对误差最大也不超过5%，因此，为了节省计算时间和保证计算精度，本文选取D₂尺寸的网格开展晃荡研究。

其次，进行模型计算精度的分析。基于Liu等^[19]与张珍^[20]的物理实验建立数值模型，其设置条件：文献[19]为简谐振动，激励频率0.583$ {\omega _0} $（$ {\omega _0} $=6.057 rad/s），激励振幅0.5 cm；文献[20]设置为50%装载率，简谐振动，激励频率1 $ {\omega _0} $（$ {\omega _0} $=6.626 rad/s），激励振幅1 cm的情况下，将数值计算结果与实验测量值进行对比。由于法拉第波激励下涉及主参数振动，将数值模拟主参数晃荡过程中自由液面的衰减历程进行快速傅里叶变化，得到的自然频率结果其数值与文献[11]一致。

由图2可知，数值模型的结果与实验数据吻合较好，模拟精度高，从而证明了该模型可以用来继续开展液舱晃荡的相关研究。

3 结果与分析 3.1 施加水平激励下波面及压力响应

由于液舱处于一阶共振频率附近时晃荡的能量较大，工程上较为关心，为更好地理解液舱共振时的动态行为，首先围绕解析一阶固有频率附近，分析施加水平激励下的波面及压力的响应情况。

如图3和图4所示，在施加横向激励时，液舱的波面响应幅值与压力极值的变化趋势为从共振频率中心向两侧衰减。通过分析横向激励下的波面及压力响应关系，发现共振频率略低于理论解析值。在本文设定的仿真工况下，矩形液舱的横波共振中心频率为圆频率11.7261 rad/s(1.87 Hz)，与解析值12.1305 rad/s(1.93 Hz)相比存在3.6%的偏差，并且处于中心频率时，晃荡过程出现了波面破碎的情况，针对这一现象，分析认为：非线性效应限制了共振效应下波幅的无限增长，波面的破碎，降低了系统的响应频率，从而使得此时外界激励在低于解析频率下便引发了强烈共振，蒋梅荣等^[21]也有类似的实验观测结果。图5为频率$ f $=1.93/1.87 Hz的水平激励下，液舱波面及压力的时序图。与其他工况相比，强水平晃荡工况的发生需要更长时间来积累能量，并具有较长的持续周期。液体受惯性与重力的影响，在液舱左右运动过程中与推动液面的舱壁和阻碍重力下落的舱壁发生2次碰撞，从而导致液舱侧壁的晃荡压力出现了明显的双峰现象，这与Liu等^[22]的观测结果一致。此外，在强水平晃荡工况下，两峰的峰曲线之间出现了纠缠现象。

3.2 施加纵向激励下波面及压力响应情况

如图6所示，液舱在受到纵向超谐波和谐波激励作用后，仅在液舱运动初期的数秒内产生了响应，随后外界激励便不再足以引发系统的动力失稳，液面进入一个持续衰减的过程，最终达到微小振幅的稳定响应状态。相比之下，纵向次谐波激励下，液舱表现出明显的主参数激振现象，液面沿舱壁迅速爬升至舱顶，并伴随有一次较大的压力响应。

通过对比图7在2种激励方式下强晃荡时波面和压力的响应情况，可以发现，在横向激励下，波面与压力的变化趋势在时间上保持一致；而在纵向激励的动力学特征上，压力的极值响应的时间点上晚于波幅，表现出更为复杂的非线性行为。

3.3 基于小波变换及二阶谱的非线性相互作用分析

由于在自适应步长下，采集的数据在时间序列上并非均匀分布，因此在进行时频域分析之前，需对数据进行插值处理。为此，本文采用分段立方插值PCHIP（Piecewise Cubic Hermite Interpolating Polynomial）方法对数据进行时序间隔为0.005 s的插值处理，以保证进行时频域分析时采样率稳定为200 Hz。插值后的效果如图8所示。

由图9中系列快速傅里叶变换的频谱图可知，当横荡液舱激励从一阶共振频率的中心向高频和低频方向变化时，波幅频谱呈现出2个显著特征：首先，此时波幅频谱具有2个主要峰值，分别对应于激励的一阶和二阶频率；其次，是其能量谱密度主要集中在一阶频率处。表明在横波晃荡过程中，波面的振动趋势与外部激励频率趋向保持一致。横荡液舱的压力频谱也展现出2个特点：一是压力频谱中主要出现3个峰值，分别对应于激励的一阶、二阶和三阶频率；二是压力频谱的能量谱密度集中在二阶频率，这与实验中观察到的压力双峰现象相吻合。而在共振频率中心，压力频谱表现出不同的特征：其能量密度峰值较为显著地低于外界激励频率，且在强共振条件下，压力频谱能量向低阶频率转移，位于一阶波幅频谱处的能量占比增加。

图9为系列激励下的Morlet小波变换图，图例代表归一化的频率强度，在系列横波晃荡频率下，Morlet小波变换的时频图服从相似的逻辑，其图谱能量密度与快速傅里叶变化的到频谱能量分布一致。同时，时频图揭示了晃荡过程中二阶频与三阶频之间存在的连续转换关系。值得注意的是，在共振中心时，小波时频图显示液舱一阶频谱的能量密度显著增强，在中心时刻达到了与二、三阶频相当的程度，同时以一、二、三阶频为核心的能量转化表现出更明显的周期性。

图10为主要尺度范围下的小波双谱，其图例色阶颜色及数值代表了不同频率之间的相干性和非线性耦合强度，1表示完全相关，0表示无相互作用。其结果表明横波激励下压力频谱的一、二、三阶频率之间存在双相干耦合关系。以水平晃荡的共振中心，即外界激励f=1.87 Hz时为例，小波双谱中存在2个强中心及1个次中心，自耦合区（3.64 Hz，3.64 Hz）对应着晃荡过程中的主要非线性特征频率，这表明压力频谱图中能量谱最高的3.64 Hz成分自身产生了强烈的非线性相互作用，导致该频率成分的能量增强并由此产生了新的四阶频率特征7.28 Hz。与此同时，随着外界激励逐步接近共振中心，双谱图表现出愈发强烈的自耦合性，自耦合性增强说明此时系统涉及二次或三次谐波现象，一阶与二阶频率之间的成分经过自耦与互耦的非线性相互作用生成了更高阶的和差频成分。

图11为纵向激励下液舱波高及压力响应的快速傅里叶变换频谱及Morlet小波变换图。结果显示，在超谐波及谐波激励下，系统的响应较小，波幅频谱主要集中在外界激励频率附近，压力频谱则呈现出近似二阶频的特征。当次谐波激励发生时，液面幅值随时间指数增长，从微小幅度逐步扩大，并最终导致系统动力失稳，短时间内产生强烈的晃荡效应。在此工况下，波幅频谱的主要能量集中在横荡共振时的强晃荡频率附近，而压力频谱则以二阶频和四阶频成分为主。此外，从小波变换时频图可以观察到，纵向激励引发的晃荡现象，其频谱能量主要集中在失稳及动力失稳的初始阶段，这反映出激励作用在早期对系统的显著影响。

图12为纵向激励下法拉第波的小波双谱。在外界激励频率为3.86 Hz的次谐波激励状态下，双谱中对应最大压力能量频谱的自耦区域最为明亮，这表明压力的二阶频成分在频谱能量中占据主导地位。同时，所在的自耦合区引发了强烈的非线性作用，这可以解释为何在频谱图中，四阶频成分的占比高于三阶频。而与横荡时的情况不同，在外界激励频率为1.93 Hz的谐波激励状态下，尽管在频谱图和时频图中难以直接观察到1 Hz的成分，但双谱图表明该成分在自耦合区外也参与了非线性作用。而在激励频率为1.29 Hz的超谐波激励状态下，除压力二阶频成分表现强烈外，其他非线性作用由于能量占比较低或持续时间极短，难以清晰辨识。从双谱图亮点的整体趋势来看，交叉耦合区的相干现象逐渐向两端扩散，表明非线性作用随着激励频率的减少逐渐减弱，且相干性在更广泛的频率范围内趋向减小。

4 结　语

本文开展了矩形液舱晃荡过程的数值仿真模拟，利用快速傅里叶变换、小波变换、小波双干谱等时频域工具对矩形液舱在横波与法拉第纵波激励影响下的晃荡响应展开分析，得到的结论如下：

1）波高及压力响应表现：在横波激励下，矩形液舱的一阶共振中心频率低于理论解析值，且压力与波高的峰值响应时刻同步。法拉第波激励下，超谐波和次谐波激励难以维持系统的持续晃荡，强烈晃荡主要由次谐波引发的主参数共振导致。同时，压力峰值响应时刻滞后于波高峰值，且不稳定阶段的晃荡强度显著高于稳定阶段。横向与纵向晃荡的压力响应均表现出双峰特征，强共振时双峰曲线存在明显的纠缠现象。

2）波高及压力时序信号的时频域表现：在频域分布上，波幅频谱的能量谱主要集中在一阶处，且与外界激励同步；而压力频谱的能量谱则主要集中在二阶处，并且在接近共振时伴随着一阶成分占比增加。在小波时频分析中，晃荡初期的不稳定阶段非线性作用最强，能量分布于一至三阶频，且在二阶与三阶频率之间存在明显的能量转移。

3）压力时序信号的相干性表现：随着横波激励接近共振中心频率，一阶频与二阶频在自耦区和交叉耦合区表现出显著的相干性，并形成了三、四阶频谱的和差频成分。此外，自耦区的强相干频率与能量谱的峰值频率一致。在纵向法拉第波共振中，随着激励频率从次谐波向谐波变化，自耦合区外的相干性逐渐减弱，最终集中在压力频谱的能量中心，次谐波激励下的非线性耦合作用尤为显著，导致了强烈的晃荡与压力响应，而在谐波与超谐波激励下，系统趋于稳定。