容器固有頻率對液體晃蕩的影響＊

寧德志，宋偉華，滕 斌，卓曉玲

（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

（杜素蘭 編輯）