容器固有频率对液体晃荡的影响＊

宁德志，宋伟华，滕 斌，卓晓玲

（1.武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉430072；2.大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连116024）

容器固有频率对液体晃荡的影响＊

宁德志1，2，宋伟华2，滕 斌2，卓晓玲2

（1.武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉430072；2.大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连116024）

针对矩形容器内液体晃荡问题，采用了时域高阶边界元方法建立自由水面满足完全非线性边界条件的数学模型。求解中采用混合欧拉－拉格朗日方法追踪流体瞬时水面，运用四阶龙格库塔方法更新下一时间步的波面和速度势。通过将计算得到的波面结果与实验数据、解析解和已发表结果对比，吻合良好，验证了本方法的准确性。进而采用谱分析方法分析了波面时间历程，得到容器各阶固有频率对液体晃荡的影响。研究发现，基频对液体晃荡的影响最大，且非线性越强，更高阶容器固有频率的影响越大。

固有频率；完全非线性；混合欧拉－拉格朗日方法；时域高阶边界元方法；谱分析

在航空航天、大型远洋船舶运输和超大型原油储存罐等工程领域中，各种容器内液体晃荡问题普遍存在，且日益得到工业界和科学界的广泛重视。尤其近些年来，由于特种船舶LNG、LPG货轮和超大油轮的出现，舱内液体晃荡问题更是引起了人们极大的关注。当容器装有部分液体时，液体在外部激励下会晃荡起来，当晃荡频率接近结构固有频率时，便会引发共振现象，对结构安全造成极大的危害。因此，固有频率对液体晃荡的影响研究有着非常重要的工程应用价值。

在液体晃荡的物理特性、自振频率下液体晃荡的描述和减晃设施抑制液体晃荡等方面，国内外许多学者已开展了大量的研究。例如，Sriram等［1］和Frandsen等［2］研究了在外界激励作用下，矩形容器中液体的二维晃荡问题；陈科等［3］模拟了矩形贮箱内液体在多种外界激励作用下的系统响应；包光伟等［4］计算了平放圆柱腔内三维液体晃动的固有频率，并将矩形容器、球腔、带“十”字隔板球形容器内的液体三维晃荡模拟结果与解析解、实验结果进行了比较；岳宝增等［5－6］、王建军等［7］利用有限元方法分别研究了在具有固有频率的激励作用下圆筒形贮腔和矩形容器中液体的大幅晃动问题；Brswal等［8］采用有限元方法研究了矩形容器内水平隔板的减晃作用；Liu等［9］采用VOF方法研究了矩形容器内垂直隔板的减晃作用。

然而在以往有关容器自振频率的晃荡研究中，更多的是考虑外部激励频率与容器一阶固有频率相近或相等时的液体晃荡问题，而一般情况下容器固有频率对液体晃荡的影响及高阶固有频率对液体晃荡贡献的分析还不是很多。本研究采用时域高阶边界元方法建立自由水面满足完全非线性边界条件的数学模型，在验证模型正确的基础上，模拟研究容器各阶固有频率对液体晃荡的影响，揭示一般情况下各阶固有频率对液体晃荡结果贡献的非线性现象。

1 数学模型

1.1 控制方程、初边界条件

如图1所示的矩形容器，长为L，水深为H。建立两个坐标系：O0X0Z0为固定的空间坐标系，OXZ为随容器运动的坐标系。在初始时刻，两坐标系重合。Z0＝0位于静水面上，且Z0轴向上为正，X0轴指向容器长度方向。容器做平面运动，位移定义为：Xb＝［xb（t），zb（t）］。

假定整个流域内流体为理想流体，则速度势φ应满足Laplace方程：

在容器侧壁，φ应满足：

图1 定义图Fig.1 Definition sketch

式中，η0表示自由表面；g表示重力加速度。

然后，利用下面的式（5）和（6），将式（3）和（4）式转化到OXZ坐标系中，写成式（7）和（8）的形式。

式中，η＝η0－zb，是OXZ坐标系中x确定时的自由表面。自由表面上φ随时间的变化由下式控制：

那么，动力学边界条件写成如下：

然后，将速度势φ分离为：

式中，u、w分别为U在x、z方向上的速度分量。将式（11）代入式（1）、（2）、（8）、（10）可得：

初始条件φ＝0进而可写成如下形式：

初始波面η（x，0）根据实际模拟的波面情况来确定。

1.2 数值求解

在整个流域内对速度势应用格林第二定理，可以得到如下边界积分方程［10］：

式中，q为源点；P为场点；C为固角系数；Γ为流域的边界，包括自由水面边界和固体边界。本研究用三节点高阶边界元离散计算域成一些曲线单元，单元内任一点的几何坐标和速度势等物理量可以用形状函数插值得到，这与曲面单元内的做法是相同的［11］。

积分方程（17）经高阶边界元离散后，可建立如下线性方程组：

其中，A和B为系数矩阵。由于积分边界是不断地随着时间变化的，在每一计算时刻都要重新建立系数矩阵，并且在每一计算时刻都要对方程进行求解。

计算中认为当前时刻物面Sn上的速度势法向导数和自由水面Sf上的速度势是已知的，根据积分方程（18）计算当前时刻物面Sn上的速度势和自由水面Sf上的速度势法向导数，然后应用四阶Runga－Kutta法，根据自由水面条件式（14）、（15）计算下一时刻的水质点位置和自由水面Sf上的速度势，再对自由水面重新划分网格，重新应用积分方程（18）计算下一时刻物面Sn上的速度势和自由水面Sf上的速度势法向导数。这样计算周而复始，直到计算结束。

2 数值计算与分析

对于如图1所示的矩形容器，其固有频率如下式：

式中，n表示固有频率阶数；对于长度L＝1.0m，水深H＝0.5m的矩形容器，可得到其各阶固有频率：ω1＝5.32s－1，ω2＝7.84s－1，ω3＝9.61s－1，ω2／ω1＝1.47，ω3／ω1＝1.81。

2.1 模型验证

为验证数值模型的准确性，首先考虑一静止容器具有余弦型初始波面η＝Acos（πx／L）作为算例。其中，A是初始波面的波幅，并且引入波陡ε＝Aω12／g来衡量非线性的强弱，容器长和静水深分别定义为L＝1.0m和H＝0.5m。经开展数值收敛性试验选定时间步长△t＝0.015s，空间步长△x＝△z＝0.025m。

图2给出了不同波陡ε情况下本文数值结果与线性解析解、二阶解析解以及其它数值方法结果的对比。当ε＝0.001 4时，非线性的影响很小，波面时间历程的波峰和波谷呈现良好的对称性，数值结果和线性解、二阶解析解均吻合良好。当ε增大到0.14时，非线性的影响明显加强。波面时间历程的波峰和波谷不再对称，波峰变高，波谷变平缓。此时，数值结果和二阶解析解仍然吻合良好；线性解析解因未考虑高阶项和数值结果、二阶解析解间出现了很大的差距。当ε继续增大到0.288时，非线性的作用更加明显。波峰变的更加高陡，波谷变的更加平缓。由于液体晃荡的强非线性明显，数值结果和二阶解析解不能很好吻合，相位、幅值都出现了差距，说明对于强非线性问题二阶解析解已不能很高的描述，这是其忽略掉的高阶项导致的。而本文数值结果与Frandsen的有限差分法计算的完全非线性数值结果［12］却仍然吻合良好，说明本文数值模型对具有初始波面情况下强非线性液体晃荡可以准确的模拟。

图2 容器左侧壁点自由水面的波面时间历程图及与其它方法结果对比Fig.2 Time series of wave elevation of free water surface at the left wall of the tank and comparisons with other solutions

接下来考虑容器做水平运动xb（t）＝asin（ωt）的情况。其中，a是运动幅值，ω是运动角频率，容器长和静水深分别为L＝2.0m，H＝1.0m，初始波面为静止水面。根据数值收敛性试验选定时间步长△t＝0.015s，空间步长△x＝△z＝0.025m。针对a＝0.018 6m，ω＝0.999ω1的情况进行了模拟。

图3给出2个时刻水槽波面分布情况，同时也给出了本文数值结果与线性解析解、实验结果［13］以及有限单元法数值结果［14］对比情况。由图可知，本研究数值结果和实验结果以及有限单元法数值结果吻合的很好，但对于自振频率下的晃动问题，线性解析解已经不能很好的描述液体晃荡的波面情况。

图4给出了容器左侧壁处自由水面的波面时间历程图及与有限单元法方法的对比。由图可知，本文数值结果与有限单元法结果吻合的仍然很好。由于容器水平运动频率ω非常接近容器一阶固有频率，液体晃荡的幅值随时间逐渐增大，而且其非线性特性越发明显。

综合对容器静止和运动两种情况的验证可知，本文数值模型是准确可靠的，能够很好的描述液体在容器内的晃荡过程。

2.2 容器固有频率影响分析

在以往的研究中，更多的是关于容器运动频率与容器一阶固有频率相等时液体晃荡的模拟和研究，而一般情况下容器固有频率对液体晃荡影响及高阶固有频率贡献的分析还不是很多。本研究针对固定容器内具有不同初始波面情况进行研究。选取容器长和水深比L／H＝2。

首先考虑初始波面为η＝Acos（πx／L），波陡ε分别为0.005，0.05，0.1和0.2时的4种情况下的液体晃荡问题。如上一节那样通过完全非线性时域模拟得到容器左侧壁处自由水面波面时间历程，然后采用傅里叶变换的方法对无量纲化的波面时间历程η／A进行谱分析，结果如图5所示。由图可知，当ε＝0.005时，非线性很弱，起主导作用的只有一阶固有频率ω1；随着波陡增大，当ε＝0.05和0.1时，除一阶固有频率ω1之外，二阶固有频率ω2和二倍频2ω1的作用突显出来。当波陡继续增大到0.2时，除上述3个频率外，三阶固有频率ω3也有贡献。同时，各固有频率间的相互作用加强，一阶、二阶固有频率间的和频ω1＋ω2、差频ω2－ω1也显现了出来。比较图5a～d可知，一阶固有频率ω1在各波陡情况下都起着主导作用，对液体晃荡的影响最大，其它固有频率和倍频、和频、差频的影响是次要的。然而，随着波陡的增大，非线性增强，次要频率的影响作用也逐渐增强，不容忽视。譬如图5d，除基频固有频率外，其它高阶固有频率、倍频、和频、差频的总贡献已经达到27%。

图5 不同波陡情况下容器左侧壁波面时间历程的谱分析结果Fig.5 Spectral analysis results for the time series of free water surface at the left wall of the tank with different wave slopes

图5中只分析了容器中一点波面变化受固有频率影响情况，为了更清楚的看出各固有频率对容器内水面上各点液体晃荡的影响，图6给出了各波陡情况下ω1、ω2和2ω1频率所对应的波面幅值分布情况。从图6a中可以看出，ω1频率所对应的波面幅值左右对称，在矩形容器两侧的幅值最大，在容器中间位置处的幅值最小，数值为0。当波陡增大时，ω1频率所对应的波面幅值并没有变化。这与图5中不同波陡情况下ω1对应的谱密度相同的结果是一致的。在图6b中，2ω1频率所对应的波面幅值以X＝L／2对称分布，且波陡越大，波面幅值也随之增大。这与图5中不同波陡情况下2ω2对应能量谱密度值的变化趋势是相同的。在图6c中，ω2频率所对应的波面幅值在波陡ε＝0.005，即非线性很弱时呈对称性分布，与图6a相同；然而随着波陡的增大，非线性增强，ω2频率所对应的波面幅值的分布也随着发生了变化，对称点向左偏移，同时对称性失衡。对比分析图6可知，容器两侧波面贡献主要取决于基频固有频率，而中间位置处基频固有频率贡献为0，主要取决于高阶固有频率和基频的多倍频。

为了使问题更具一般性，下面考虑固定容器初始波面为直线型η＝A×（2.0×x／L－1.0）的情况，容器尺寸同前例。

考虑波陡A／H＝0.005，0.025和0.037的3种情况，对无量纲化的波面时间历程η／A采用傅里叶变换方法进行谱分析，结果如图7所示。从图中可以得到与图5相似的规律，基频固有频率贡献仍然占主导作用，而随着非线性的增强，更多的高阶固有频率开始发生贡献。

3 结 论

本研究利用时域高阶边界元方法建立二维矩形容器中液体晃荡问题的完全非线性数学模型。通过对具有初始波面的静止容器和做正弦运动的容器内液体晃荡问题的模拟，及与解析解、试验数据和其它数值方法结果的对比，验证了本研究模型的正确性。进而通过傅里叶变换方法对容器左侧壁处波面时间历程进行谱分析，定量研究了固有频率对液体晃荡的影响。研究表明：除占主导地位的基频固有频率贡献外，随着波陡的增大，非线性增强，其它阶次固有频率的作用也相继显现；同时，波陡的增大使各固有频率间的相互作用加强，激发出了倍频、和频以及差频。在本研究中，除基频外的其它频率对波面的贡献最大可达到27%，这在工程应用中应引起足够的重视。

（References）：

［1］SRIRAM V，SANNASIRAJ S A，SUNDAR V.Numerical simulation of 2Dsloshing waves due to horizontal and vertical random excitation［J］.Applied Ocean Research，2006，28：19－32.

［2］FRANDSEN J B.Sloshing motions in excited tanks［J］.Journal of Computational Physics，2004，196：53－87.

［3］CHEN K，LI J F，WANG T S.Nonlinear dynamics modeling and analysis of liquid sloshing in rectangle tank［J］.Acta Mechanica Sinica，2005，37（3）：339－345.陈科，李俊峰，王天舒.矩形贮箱内液体非线性晃动动力学建模与分析 ［J］.力学学报，2005，37（3）：339－345.

［4］BAO G W，WANG Z W.Finite element method for eigen problem of liquid 3Dsloshing［J］.Chinese Quarterly of Mechanics，2003，24（2）：185－190.包光伟，王政伟.液体三维晃动特征问题的有限元数值计算方法 ［J］.力学季刊，2003，24（2）：185－190.

［5］YUE B Z.Three dimensional large amplitude liquid sloshing under pitching excitation［J］.Acta Mechanica Sinica，2005，37（2）：199－203.岳宝增.信仰激励下三维液体大幅晃动问题研究 ［J］.力学学报，2005，37（2）：199－203.

［6］YUE B Z，LIU Y Z，WANG Z L.ALE finite element method for three－dimensional large ampligude liquid sloshing using fractional step method［J］.Chinese Journal of Applied Mechanics，2001，18（1）：110－115.岳宝增，刘延柱，王照林.三维液体非线性晃动动力学特性的数值模拟 ［J］.应用力学学报，2001，18（1）：110－115.

［7］WANG J J，LU M W，ZHANG X，et al.A classical Galerkin finite element method for the large fluid sloshing［J］.Journal of hydrodynamics：Ser.A，2001，16（3）：390－395.王建军，陆明万，张雄，等.自由液面流体大晃动有限元方法 ［J］.水动力学研究与进展：A辑，2001，16（3）：390－395.

［8］BISWAL K C，BHATTACHARYYA S K，SINHA P K.Nonlinear sloshing in partially liquid filled containers with baffles［J］.Int.J.Numer.Meth.Engng，2006，68：317－337.

［9］LIU D M，LIN P Z.Three－dimensional liquid sloshing in a tank with baffles［J］.Ocean Engineering，2009，36：202－212.

［10］LI Y C，TENG B.Wave action on maritime structhures：the second edition［M］.Beijing：Ocean Press，2002.李玉成，滕斌.波浪对海上建筑物的作用：第2版 ［M］.北京：海洋出版社，2002.

［11］NING D Z.The application of fast multipole boundary element method for fully nonlinear water wave problems［D］.Dalian：Dalian University of Technology，2005.宁德志.快速多极子边界元方法在完全非线性水波问题中的应用 ［D］.大连：大连理工大学，2005.

［12］FRANDSEN J B.Sloshing motions in excited tanks［J］.Journal of Computational Physics，2004，196：53－87.

［13］OKAMOTO T，KAWAHARA M.Two－dimensional sloshing analysis by Lagrangian finite element method［J］.International Journal for Numerical Methods in Fluids，1990，11：453－477.

［14］WU G X，MA O W，TAYLOR R E.Numerical simulation of sloshing waves in a 3Dtank based on a finite element method［J］.Applied Ocean Research，1998，20：337－355.

Effect of Natural Frequencies of Container on Fluid Sloshing

NING De－zhi1，2，SONG Wei－hua2，TENG Bin2，ZHUO Xiao－ling2
（1.State Key Laboratory of Water Resources Hydropower Engineering Science，Wuhan University，Wuhan 430072，China；2.State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian 116024，China）

To solve the problem of liquid sloshing in a rectangular container，a time－domain higher－order boundary element method was adopted to establish a mathematical model with fully nonlinear boundary conditions satisfied by free water surface.In the solving process，the mixed Eulerian－Lagrangian technique was applied to track the transient liquid surface and the 4th－order Runge－Kutta method was used to refresh wave elevation and velocity potential on the free water surface at each time step.The calculated results were compared with experimental data，analytical solutions and published results respectively.Good agreements among them were obtained and the accuracy of present model was verified.After the Fourier Transformation method was adopted to study the time series of wave elevation，the contributions from various order natural frequencies of container for sloshing were analyzed.It shows that the fundamental frequency contributes strongest effect on the sloshing，and the higher－order natural frequencies have more effect with the nonlinearity increased.

natural frequency；fully nonlinear；Mixed Eulerian－Lagrangian method；Time－domain higherorder boundary element method；spectrum analysis

December 29，2010

P75

A

1671－6647（2012）01－0045－09

2010－12－29

国家自然科学创新群体基金——海洋环境灾害作用与结构安全防护（50921001）；国家自然科学基金面上项目——极值波浪与水流混合与海洋结构物作用的模拟研究（5179028）；水资源与水电工程科学国家重点实验室开放基金——复杂边界条件下溃坝问题的实时模拟研究（2009B057）

宁德志（1975－），男，黑龙江五常人，副教授，主要从事海洋工程中非线性波浪与结构物作用方面研究.E－mail：dzning＠dlut.edu.cn

（杜素兰 编辑）