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

（1.武漢大學(xué) 水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢430072；2.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連116024）

容器固有頻率對(duì)液體晃蕩的影響＊

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

（1.武漢大學(xué) 水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢430072；2.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連116024）

針對(duì)矩形容器內(nèi)液體晃蕩問題，采用了時(shí)域高階邊界元方法建立自由水面滿足完全非線性邊界條件的數(shù)學(xué)模型。求解中采用混合歐拉－拉格朗日方法追蹤流體瞬時(shí)水面，運(yùn)用四階龍格庫塔方法更新下一時(shí)間步的波面和速度勢(shì)。通過將計(jì)算得到的波面結(jié)果與實(shí)驗(yàn)數(shù)據(jù)、解析解和已發(fā)表結(jié)果對(duì)比，吻合良好，驗(yàn)證了本方法的準(zhǔn)確性。進(jìn)而采用譜分析方法分析了波面時(shí)間歷程，得到容器各階固有頻率對(duì)液體晃蕩的影響。研究發(fā)現(xiàn)，基頻對(duì)液體晃蕩的影響最大，且非線性越強(qiáng)，更高階容器固有頻率的影響越大。

固有頻率；完全非線性；混合歐拉－拉格朗日方法；時(shí)域高階邊界元方法；譜分析

在航空航天、大型遠(yuǎn)洋船舶運(yùn)輸和超大型原油儲(chǔ)存罐等工程領(lǐng)域中，各種容器內(nèi)液體晃蕩問題普遍存在，且日益得到工業(yè)界和科學(xué)界的廣泛重視。尤其近些年來，由于特種船舶LNG、LPG貨輪和超大油輪的出現(xiàn)，艙內(nèi)液體晃蕩問題更是引起了人們極大的關(guān)注。當(dāng)容器裝有部分液體時(shí)，液體在外部激勵(lì)下會(huì)晃蕩起來，當(dāng)晃蕩頻率接近結(jié)構(gòu)固有頻率時(shí)，便會(huì)引發(fā)共振現(xiàn)象，對(duì)結(jié)構(gòu)安全造成極大的危害。因此，固有頻率對(duì)液體晃蕩的影響研究有著非常重要的工程應(yīng)用價(jià)值。

在液體晃蕩的物理特性、自振頻率下液體晃蕩的描述和減晃設(shè)施抑制液體晃蕩等方面，國內(nèi)外許多學(xué)者已開展了大量的研究。例如，Sriram等［1］和Frandsen等［2］研究了在外界激勵(lì)作用下，矩形容器中液體的二維晃蕩問題；陳科等［3］模擬了矩形貯箱內(nèi)液體在多種外界激勵(lì)作用下的系統(tǒng)響應(yīng)；包光偉等［4］計(jì)算了平放圓柱腔內(nèi)三維液體晃動(dòng)的固有頻率，并將矩形容器、球腔、帶“十”字隔板球形容器內(nèi)的液體三維晃蕩模擬結(jié)果與解析解、實(shí)驗(yàn)結(jié)果進(jìn)行了比較；岳寶增等［5－6］、王建軍等［7］利用有限元方法分別研究了在具有固有頻率的激勵(lì)作用下圓筒形貯腔和矩形容器中液體的大幅晃動(dòng)問題；Brswal等［8］采用有限元方法研究了矩形容器內(nèi)水平隔板的減晃作用；Liu等［9］采用VOF方法研究了矩形容器內(nèi)垂直隔板的減晃作用。

然而在以往有關(guān)容器自振頻率的晃蕩研究中，更多的是考慮外部激勵(lì)頻率與容器一階固有頻率相近或相等時(shí)的液體晃蕩問題，而一般情況下容器固有頻率對(duì)液體晃蕩的影響及高階固有頻率對(duì)液體晃蕩貢獻(xiàn)的分析還不是很多。本研究采用時(shí)域高階邊界元方法建立自由水面滿足完全非線性邊界條件的數(shù)學(xué)模型，在驗(yàn)證模型正確的基礎(chǔ)上，模擬研究容器各階固有頻率對(duì)液體晃蕩的影響，揭示一般情況下各階固有頻率對(duì)液體晃蕩結(jié)果貢獻(xiàn)的非線性現(xiàn)象。

1 數(shù)學(xué)模型

1.1 控制方程、初邊界條件

如圖1所示的矩形容器，長(zhǎng)為L(zhǎng)，水深為H。建立兩個(gè)坐標(biāo)系：O0X0Z0為固定的空間坐標(biāo)系，OXZ為隨容器運(yùn)動(dòng)的坐標(biāo)系。在初始時(shí)刻，兩坐標(biāo)系重合。Z0＝0位于靜水面上，且Z0軸向上為正，X0軸指向容器長(zhǎng)度方向。容器做平面運(yùn)動(dòng)，位移定義為：Xb＝［xb（t），zb（t）］。

假定整個(gè)流域內(nèi)流體為理想流體，則速度勢(shì)φ應(yīng)滿足Laplace方程：

在容器側(cè)壁，φ應(yīng)滿足：

圖1 定義圖Fig.1 Definition sketch

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

然后，利用下面的式（5）和（6），將式（3）和（4）式轉(zhuǎn)化到OXZ坐標(biāo)系中，寫成式（7）和（8）的形式。

式中，η＝η0－zb，是OXZ坐標(biāo)系中x確定時(shí)的自由表面。自由表面上φ隨時(shí)間的變化由下式控制：

那么，動(dòng)力學(xué)邊界條件寫成如下：

然后，將速度勢(shì)φ分離為：

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

初始條件φ＝0進(jìn)而可寫成如下形式：

初始波面η（x，0）根據(jù)實(shí)際模擬的波面情況來確定。

1.2 數(shù)值求解

在整個(gè)流域內(nèi)對(duì)速度勢(shì)應(yīng)用格林第二定理，可以得到如下邊界積分方程［10］：

式中，q為源點(diǎn)；P為場(chǎng)點(diǎn)；C為固角系數(shù)；Γ為流域的邊界，包括自由水面邊界和固體邊界。本研究用三節(jié)點(diǎn)高階邊界元離散計(jì)算域成一些曲線單元，單元內(nèi)任一點(diǎn)的幾何坐標(biāo)和速度勢(shì)等物理量可以用形狀函數(shù)插值得到，這與曲面單元內(nèi)的做法是相同的［11］。

積分方程（17）經(jīng)高階邊界元離散后，可建立如下線性方程組：

其中，A和B為系數(shù)矩陣。由于積分邊界是不斷地隨著時(shí)間變化的，在每一計(jì)算時(shí)刻都要重新建立系數(shù)矩陣，并且在每一計(jì)算時(shí)刻都要對(duì)方程進(jìn)行求解。

計(jì)算中認(rèn)為當(dāng)前時(shí)刻物面Sn上的速度勢(shì)法向?qū)?shù)和自由水面Sf上的速度勢(shì)是已知的，根據(jù)積分方程（18）計(jì)算當(dāng)前時(shí)刻物面Sn上的速度勢(shì)和自由水面Sf上的速度勢(shì)法向?qū)?shù)，然后應(yīng)用四階Runga－Kutta法，根據(jù)自由水面條件式（14）、（15）計(jì)算下一時(shí)刻的水質(zhì)點(diǎn)位置和自由水面Sf上的速度勢(shì)，再對(duì)自由水面重新劃分網(wǎng)格，重新應(yīng)用積分方程（18）計(jì)算下一時(shí)刻物面Sn上的速度勢(shì)和自由水面Sf上的速度勢(shì)法向?qū)?shù)。這樣計(jì)算周而復(fù)始，直到計(jì)算結(jié)束。

2 數(shù)值計(jì)算與分析

對(duì)于如圖1所示的矩形容器，其固有頻率如下式：

式中，n表示固有頻率階數(shù)；對(duì)于長(zhǎng)度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 模型驗(yàn)證

為驗(yàn)證數(shù)值模型的準(zhǔn)確性，首先考慮一靜止容器具有余弦型初始波面η＝Acos（πx／L）作為算例。其中，A是初始波面的波幅，并且引入波陡ε＝Aω12／g來衡量非線性的強(qiáng)弱，容器長(zhǎng)和靜水深分別定義為L(zhǎng)＝1.0m和H＝0.5m。經(jīng)開展數(shù)值收斂性試驗(yàn)選定時(shí)間步長(zhǎng)△t＝0.015s，空間步長(zhǎng)△x＝△z＝0.025m。

圖2給出了不同波陡ε情況下本文數(shù)值結(jié)果與線性解析解、二階解析解以及其它數(shù)值方法結(jié)果的對(duì)比。當(dāng)ε＝0.001 4時(shí)，非線性的影響很小，波面時(shí)間歷程的波峰和波谷呈現(xiàn)良好的對(duì)稱性，數(shù)值結(jié)果和線性解、二階解析解均吻合良好。當(dāng)ε增大到0.14時(shí)，非線性的影響明顯加強(qiáng)。波面時(shí)間歷程的波峰和波谷不再對(duì)稱，波峰變高，波谷變平緩。此時(shí)，數(shù)值結(jié)果和二階解析解仍然吻合良好；線性解析解因未考慮高階項(xiàng)和數(shù)值結(jié)果、二階解析解間出現(xiàn)了很大的差距。當(dāng)ε繼續(xù)增大到0.288時(shí)，非線性的作用更加明顯。波峰變的更加高陡，波谷變的更加平緩。由于液體晃蕩的強(qiáng)非線性明顯，數(shù)值結(jié)果和二階解析解不能很好吻合，相位、幅值都出現(xiàn)了差距，說明對(duì)于強(qiáng)非線性問題二階解析解已不能很高的描述，這是其忽略掉的高階項(xiàng)導(dǎo)致的。而本文數(shù)值結(jié)果與Frandsen的有限差分法計(jì)算的完全非線性數(shù)值結(jié)果［12］卻仍然吻合良好，說明本文數(shù)值模型對(duì)具有初始波面情況下強(qiáng)非線性液體晃蕩可以準(zhǔn)確的模擬。

圖2 容器左側(cè)壁點(diǎn)自由水面的波面時(shí)間歷程圖及與其它方法結(jié)果對(duì)比Fig.2 Time series of wave elevation of free water surface at the left wall of the tank and comparisons with other solutions

接下來考慮容器做水平運(yùn)動(dòng)xb（t）＝asin（ωt）的情況。其中，a是運(yùn)動(dòng)幅值，ω是運(yùn)動(dòng)角頻率，容器長(zhǎng)和靜水深分別為L(zhǎng)＝2.0m，H＝1.0m，初始波面為靜止水面。根據(jù)數(shù)值收斂性試驗(yàn)選定時(shí)間步長(zhǎng)△t＝0.015s，空間步長(zhǎng)△x＝△z＝0.025m。針對(duì)a＝0.018 6m，ω＝0.999ω1的情況進(jìn)行了模擬。

圖3給出2個(gè)時(shí)刻水槽波面分布情況，同時(shí)也給出了本文數(shù)值結(jié)果與線性解析解、實(shí)驗(yàn)結(jié)果［13］以及有限單元法數(shù)值結(jié)果［14］對(duì)比情況。由圖可知，本研究數(shù)值結(jié)果和實(shí)驗(yàn)結(jié)果以及有限單元法數(shù)值結(jié)果吻合的很好，但對(duì)于自振頻率下的晃動(dòng)問題，線性解析解已經(jīng)不能很好的描述液體晃蕩的波面情況。

圖4給出了容器左側(cè)壁處自由水面的波面時(shí)間歷程圖及與有限單元法方法的對(duì)比。由圖可知，本文數(shù)值結(jié)果與有限單元法結(jié)果吻合的仍然很好。由于容器水平運(yùn)動(dòng)頻率ω非常接近容器一階固有頻率，液體晃蕩的幅值隨時(shí)間逐漸增大，而且其非線性特性越發(fā)明顯。

綜合對(duì)容器靜止和運(yùn)動(dòng)兩種情況的驗(yàn)證可知，本文數(shù)值模型是準(zhǔn)確可靠的，能夠很好的描述液體在容器內(nèi)的晃蕩過程。

2.2 容器固有頻率影響分析

在以往的研究中，更多的是關(guān)于容器運(yùn)動(dòng)頻率與容器一階固有頻率相等時(shí)液體晃蕩的模擬和研究，而一般情況下容器固有頻率對(duì)液體晃蕩影響及高階固有頻率貢獻(xiàn)的分析還不是很多。本研究針對(duì)固定容器內(nèi)具有不同初始波面情況進(jìn)行研究。選取容器長(zhǎng)和水深比L／H＝2。

首先考慮初始波面為η＝Acos（πx／L），波陡ε分別為0.005，0.05，0.1和0.2時(shí)的4種情況下的液體晃蕩問題。如上一節(jié)那樣通過完全非線性時(shí)域模擬得到容器左側(cè)壁處自由水面波面時(shí)間歷程，然后采用傅里葉變換的方法對(duì)無量綱化的波面時(shí)間歷程η／A進(jìn)行譜分析，結(jié)果如圖5所示。由圖可知，當(dāng)ε＝0.005時(shí)，非線性很弱，起主導(dǎo)作用的只有一階固有頻率ω1；隨著波陡增大，當(dāng)ε＝0.05和0.1時(shí)，除一階固有頻率ω1之外，二階固有頻率ω2和二倍頻2ω1的作用突顯出來。當(dāng)波陡繼續(xù)增大到0.2時(shí)，除上述3個(gè)頻率外，三階固有頻率ω3也有貢獻(xiàn)。同時(shí)，各固有頻率間的相互作用加強(qiáng)，一階、二階固有頻率間的和頻ω1＋ω2、差頻ω2－ω1也顯現(xiàn)了出來。比較圖5a～d可知，一階固有頻率ω1在各波陡情況下都起著主導(dǎo)作用，對(duì)液體晃蕩的影響最大，其它固有頻率和倍頻、和頻、差頻的影響是次要的。然而，隨著波陡的增大，非線性增強(qiáng)，次要頻率的影響作用也逐漸增強(qiáng)，不容忽視。譬如圖5d，除基頻固有頻率外，其它高階固有頻率、倍頻、和頻、差頻的總貢獻(xiàn)已經(jīng)達(dá)到27%。

圖5 不同波陡情況下容器左側(cè)壁波面時(shí)間歷程的譜分析結(jié)果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中只分析了容器中一點(diǎn)波面變化受固有頻率影響情況，為了更清楚的看出各固有頻率對(duì)容器內(nèi)水面上各點(diǎn)液體晃蕩的影響，圖6給出了各波陡情況下ω1、ω2和2ω1頻率所對(duì)應(yīng)的波面幅值分布情況。從圖6a中可以看出，ω1頻率所對(duì)應(yīng)的波面幅值左右對(duì)稱，在矩形容器兩側(cè)的幅值最大，在容器中間位置處的幅值最小，數(shù)值為0。當(dāng)波陡增大時(shí)，ω1頻率所對(duì)應(yīng)的波面幅值并沒有變化。這與圖5中不同波陡情況下ω1對(duì)應(yīng)的譜密度相同的結(jié)果是一致的。在圖6b中，2ω1頻率所對(duì)應(yīng)的波面幅值以X＝L／2對(duì)稱分布，且波陡越大，波面幅值也隨之增大。這與圖5中不同波陡情況下2ω2對(duì)應(yīng)能量譜密度值的變化趨勢(shì)是相同的。在圖6c中，ω2頻率所對(duì)應(yīng)的波面幅值在波陡ε＝0.005，即非線性很弱時(shí)呈對(duì)稱性分布，與圖6a相同；然而隨著波陡的增大，非線性增強(qiáng)，ω2頻率所對(duì)應(yīng)的波面幅值的分布也隨著發(fā)生了變化，對(duì)稱點(diǎn)向左偏移，同時(shí)對(duì)稱性失衡。對(duì)比分析圖6可知，容器兩側(cè)波面貢獻(xiàn)主要取決于基頻固有頻率，而中間位置處基頻固有頻率貢獻(xiàn)為0，主要取決于高階固有頻率和基頻的多倍頻。

為了使問題更具一般性，下面考慮固定容器初始波面為直線型η＝A×（2.0×x／L－1.0）的情況，容器尺寸同前例。

考慮波陡A／H＝0.005，0.025和0.037的3種情況，對(duì)無量綱化的波面時(shí)間歷程η／A采用傅里葉變換方法進(jìn)行譜分析，結(jié)果如圖7所示。從圖中可以得到與圖5相似的規(guī)律，基頻固有頻率貢獻(xiàn)仍然占主導(dǎo)作用，而隨著非線性的增強(qiáng)，更多的高階固有頻率開始發(fā)生貢獻(xiàn)。

3 結(jié) 論

本研究利用時(shí)域高階邊界元方法建立二維矩形容器中液體晃蕩問題的完全非線性數(shù)學(xué)模型。通過對(duì)具有初始波面的靜止容器和做正弦運(yùn)動(dòng)的容器內(nèi)液體晃蕩問題的模擬，及與解析解、試驗(yàn)數(shù)據(jù)和其它數(shù)值方法結(jié)果的對(duì)比，驗(yàn)證了本研究模型的正確性。進(jìn)而通過傅里葉變換方法對(duì)容器左側(cè)壁處波面時(shí)間歷程進(jìn)行譜分析，定量研究了固有頻率對(duì)液體晃蕩的影響。研究表明：除占主導(dǎo)地位的基頻固有頻率貢獻(xiàn)外，隨著波陡的增大，非線性增強(qiáng)，其它階次固有頻率的作用也相繼顯現(xiàn)；同時(shí)，波陡的增大使各固有頻率間的相互作用加強(qiáng)，激發(fā)出了倍頻、和頻以及差頻。在本研究中，除基頻外的其它頻率對(duì)波面的貢獻(xiàn)最大可達(dá)到27%，這在工程應(yīng)用中應(yīng)引起足夠的重視。

（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.陳科，李俊峰，王天舒.矩形貯箱內(nèi)液體非線性晃動(dòng)動(dòng)力學(xué)建模與分析 ［J］.力學(xué)學(xué)報(bào)，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.包光偉，王政偉.液體三維晃動(dòng)特征問題的有限元數(shù)值計(jì)算方法 ［J］.力學(xué)季刊，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.岳寶增.信仰激勵(lì)下三維液體大幅晃動(dòng)問題研究 ［J］.力學(xué)學(xué)報(bào)，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.岳寶增，劉延柱，王照林.三維液體非線性晃動(dòng)動(dòng)力學(xué)特性的數(shù)值模擬 ［J］.應(yīng)用力學(xué)學(xué)報(bào)，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.王建軍，陸明萬，張雄，等.自由液面流體大晃動(dòng)有限元方法 ［J］.水動(dòng)力學(xué)研究與進(jìn)展：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.李玉成，滕斌.波浪對(duì)海上建筑物的作用：第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.寧德志.快速多極子邊界元方法在完全非線性水波問題中的應(yīng)用 ［D］.大連：大連理工大學(xué)，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－h(huán)ua2，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

國家自然科學(xué)創(chuàng)新群體基金——海洋環(huán)境災(zāi)害作用與結(jié)構(gòu)安全防護(hù)（50921001）；國家自然科學(xué)基金面上項(xiàng)目——極值波浪與水流混合與海洋結(jié)構(gòu)物作用的模擬研究（5179028）；水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室開放基金——復(fù)雜邊界條件下潰壩問題的實(shí)時(shí)模擬研究（2009B057）

寧德志（1975－），男，黑龍江五常人，副教授，主要從事海洋工程中非線性波浪與結(jié)構(gòu)物作用方面研究.E－mail：dzning＠dlut.edu.cn

（杜素蘭 編輯）