李遇春，張 皓

（同濟(jì)大學(xué) 水利工程系，上海 200092）

晃動是一種常見的發(fā)生在貯液結(jié)構(gòu)中的物理現(xiàn)象。在土木、水利、化工、造船及航空航天工程中，常常遇到各種形狀的貯液結(jié)構(gòu)，如凈水池、水塔、水庫、渡槽、儲油罐、液化氣運輸船以及燃料箱等，為保證貯液結(jié)構(gòu)的可靠性，在結(jié)構(gòu)設(shè)計時，需要首先了解容器內(nèi)液體的晃動模態(tài)。

Dodge等［1－3］在他們的專著中對流體晃動的一般理論及工程應(yīng)用進(jìn)行了全面而系統(tǒng)的闡述。流體晃動模態(tài)分析方法有：變分法［4－5］，積分方程法［6］，保角映射法［7－8］、雙坐標(biāo)法［9］、其它解析方法［10］以及諸多的數(shù)值方法［11］，文獻(xiàn)［12－13］對三維晃動模態(tài)分析方法進(jìn)行了分析總結(jié)。

在土木、水利工程領(lǐng)域，許多三維晃動問題常常可按二維晃動問題進(jìn)行處理，例如：不同形狀的TLD（調(diào)頻液體阻尼器）、大型渡槽的內(nèi)流體的橫向晃動等都可以按二維晃動進(jìn)行分析。二維矩形截面內(nèi)液體晃動的模態(tài)問題已得到了充分的研究，然而對于非矩形截面內(nèi)的二維晃動模態(tài)問題研究不夠充分，所采用的計算方法不具一般性，僅用于解決特定的問題，例如：Hasheminejad等［8］采用保角映射的方法研究了半橢圓截面內(nèi)液體晃動模態(tài)的求解方法，但其計算公式只能適合半橢圓容器，最近Li等［14］基于Ritz方法提出了一個針對不同截面內(nèi)的二維晃動頻率的簡單近似解析方法，但計算精度有限，與其它解析與數(shù)值解的誤差在5%的范圍內(nèi)，若希望得到更高精度的頻率解時，這個方法不再適用，這個方法難以獲得可靠的晃動振型。

本文亦基于Ritz方法，通過引入Green公式，借助矩形截面晃動模態(tài)的精確解，針對任意二維截面形狀，采用統(tǒng)一的基函數(shù)，避免針對不同截面（由于邊界條件的不同）需要引入不同基函數(shù)的復(fù)雜性，使得Ritz方法在理論上適用于任意截面的二維晃動模態(tài)分析，可同時得到可靠的晃動頻率與振型，計算精度可以通過選取基函數(shù)的多寡進(jìn)行調(diào)整，可完全滿足工程計算要求。

1 二維自由晃動基本方程

如圖1所示為任意的二維（對稱）形狀截面，建立直角坐標(biāo)系OXZ，符號Ω，Sf和Sw分別表示靜止液體區(qū)域，液體自由表面和濕邊界。x軸與液體自由表面重合，z軸垂直向上并通過自由表面的中點。自由表面的寬度為2a，容器假定為剛體。

圖1 任意形狀截面Fig．1 An arbitrary-section

基于線性勢理論，不考慮表面張力作用，假定流體為無粘，無旋且不可壓縮的液體，假設(shè)液體自由表面的晃動幅度較小，那么其晃動問題可由下列的線性化方程來描述：

式中：（x，z，t）為速度勢函數(shù)，hf（x，z，t）為自由表面上的波高函數(shù)，t為時間，g為重力加速度，n為濕邊界的外法線，／n為濕邊界外法線方向的導(dǎo)數(shù)。

為了求解方程組（1）－（4），可以假定方程具有以下形式的解：

式中是晃動（圓）頻率，函數(shù)（x，z）和f（x，z）分別表示函數(shù) Φ（x，z，t）和 hf（x，z，t）的幅值，將表達(dá)式（5），（6）代入方程組（1）－（4），同時合并方程（3）（4），可以得到如下的特征值問題：

式中：特征值λ與晃動（圓）頻率ω之間的關(guān)系為：

根據(jù)圖1的液體區(qū)域Ω，引入平面Green定理：

式中：S＝Sw＋Sf，P與Q為封閉區(qū)域Ω內(nèi)連續(xù)可微的函數(shù)，若定義：

將式（12）代入方程（11），并利用方程（7）得到：

式中：s表示曲線坐標(biāo)，在邊界S上沿逆時針方向為正。τx＝cos（τ，x）及 τz＝cos（τ，z）為切向量 τ（見圖 1）的方向余弦，切向量與外法向量的方向余弦有如下的關(guān)系：

式中：nx＝cos（n，x）與 nz＝cos（n，z）為外法向量 n（見圖1）的方向余弦。將式（14）代入方程（13），將S上的積分分解為邊界Sf與Sw上的和，利用方程（8）和（9），最后可以得到：

式（15）為晃動頻率ω在二維情形下的Rayleigh商表達(dá)式。

2 Ritz求解方法

式（15）可以改寫為：

利用Ritz方法求解方程（16），可將函數(shù)表達(dá)為一組完備基函數(shù)的線性組合：

其中：αT＝（α1，…，αm）為待定系數(shù)向量

＝（φ1，…，φm）為基函數(shù)向量

將表達(dá)式（17）代入式（16）可以得到：

如果記：

則式（18）可以寫為：

方程（20）為廣義特征值問題，矩陣A，B為實對稱矩陣，若基函數(shù)向量 φT＝（φ1，…，φm）已知，解此特征值問題，可求得m個特征值與特征向量為λj，αTj（j＝1，2，…，m），根據(jù)方程（10）（17）可以得到流體系統(tǒng)的前m階晃動頻率與振型函數(shù)。

3 Ritz基函數(shù)的構(gòu)造

3.1 矩形截面基函數(shù)及其模態(tài)解

如圖2所示矩形截面的靜止流體寬度為2a0，深度為H。

圖2 矩形截面Fig．2 A rectangular section

引入無量綱變量 x＝X／a0，z＝Z／a0，無量綱液深h＝H／a0，構(gòu)造基函數(shù) φj（x，z），可將其分離為兩個獨立函數(shù) ξj（x）和 ηj（z）之積：

顯然：函數(shù)ξj（x）表示自由液面的第j階模態(tài)振型，基函數(shù) φj（x，z）應(yīng)滿足方程（7）～（9），將 φj（x，z）代入方程（7）～（9），求解之，可以得到基函數(shù) φj（x，z）及模態(tài)頻率ω2j如下：

式中：

可以發(fā)現(xiàn)式（22），（23）滿足 Laplace方程（7）及邊界條件式（8），（9），為特征值問題的精確解析解。其第j階自由表面的模態(tài)函數(shù)分別為（反對稱模態(tài)）及 cosκjx（正對稱模態(tài)）。

3.2 任意截面基函數(shù)及其模態(tài)解

如圖3，記 Ω0為矩形計算區(qū)域，Sf0和Sw0分別為其自由表面和濕邊界，其相對應(yīng)的基函數(shù)記為φi，特征值為需要求解的任意計算域Ω位于矩形區(qū)域Ω0中，并且兩者的自由表面Sf和Sf0相重合，為矩形區(qū)域Ω0與待求解域Ω之差，ΔSf為兩者自由表面Sf0以及Sf之差。

對于任意截面Ω內(nèi)的液體，若采用的Ritz方法求解，根據(jù)方程（19）（20），矩陣A以及B中元素aik和bik的表達(dá)式可以寫為以下形式：

圖3 任意截面與矩形截面Fig．3 Arbitrary and rectangular sections

引入第一格林（Green）公式的二維形式：

式中：S＝Sw＋Sf，u與 v為封閉區(qū)域 Ω內(nèi)連續(xù)可微的函數(shù)。將方程（26）應(yīng)用于矩形計算域Ω0，并設(shè)u＝φi，v＝φk代入式（26），并利用邊界條件（8）和（9）得：

考慮到指標(biāo)i和k具有可交換性，由式（27）可知，函數(shù)φi在自由表面Sf0上具有如下的正交性：

再次利用第一格林公式（26）并考慮邊界條件（8）（9）可以得到：

利用式（27）～（29），式（25）可以重新寫為以下形式：

其中

由方程（30）可以看出，對于具有任意形狀的液體區(qū)域Ω，采用Ritz方法求解晃動模態(tài)時，可統(tǒng)一采用矩形區(qū)域（截面）的基函數(shù)（式（22））進(jìn)行求解，其廣義特征值矩陣的系數(shù)可按式（30）計算，從而避免了針對不同截面（由于邊界條件的不同）需分別構(gòu)造基函數(shù)的復(fù)雜性。

4 數(shù)值算例

4．1 半頂角為45°的三角形截面

文獻(xiàn)［2］給出了半頂角為45°的三角形渠道內(nèi)流體晃動頻率的近似解析解。截面尺寸如圖4所示，其中水深為H，自由液面的半寬值為a，容器側(cè)壁間半頂角θ為 45°。

圖4 半頂角為45°的三角形截面渠道Fig 4．Canal with 45°-walls

選取式（22）為Ritz方法的基函數(shù)（取150項），矩陣A與B的系數(shù)按式（30）計算，應(yīng)用Ritz方法得到歸一化的一階反對稱以及正對稱頻率分別為1．019 0、1．574 2。文獻(xiàn)［2］給出的近似解分別為1．000 0、1．524 4。本文 Ritz解答與文獻(xiàn)［2］的近似解最大相對誤差約3%，兩者吻合良好。

4．2 半橢圓截面

選取文獻(xiàn)［8］半橢圓截面液體作為計算實例，半橢圓容器無蓋板及有蓋板的的幾何尺寸如圖5所示。液體所占體積為容器容量的一半，符號a′，b′，L分別表示橢圓的半長軸、半短軸以及蓋板長度。符號a則是自由液面靜止時的半寬值。對于無蓋板的橢圓容器，蓋板長度L為零。

同樣選取式（22）為Ritz方法的基函數(shù)，矩陣A與B的系數(shù)按式（30）計算。

圖5 有蓋板以及無蓋板的半橢圓截面Fig．5 Baffled and un-baffled half elliptical sections

本文Ritz方法（基函數(shù)所取項數(shù)m為23）計算的前三階反對稱及正對稱正規(guī)化的晃動頻率與文獻(xiàn)［8］保角映射方法的結(jié)果都列于表1，通過比較發(fā)現(xiàn)，本文和文獻(xiàn)［8］結(jié)果的最大相對誤差約為1%，由于正規(guī)化頻率為頻率平方的倍數(shù)，且本文結(jié)果與文獻(xiàn)［8］的結(jié)果很接近，所以實際頻率的最大相對誤差大約為0．5%，兩者結(jié)果吻合很好。

為了考察此時Ritz方法的收斂性，取兩個無蓋板的半橢圓截面為例，截面尺寸分別為a′＝1．0 m，b′＝0．1 m及 a′＝1．0 m，b′＝0．8 m，短長軸的比值分別為b′／a′＝0．1，0．8。圖 6、7分別顯示了二個半橢圓容器內(nèi)流體前3階反對稱與對稱晃動頻率隨基函數(shù)項數(shù)m變化的圖線，由圖可見，對于本算例而言，Ritz方法收斂很快，只需取5～6項即可獲得較高的計算精度，可滿足工程使用要求。

通常情況下，采用本文Ritz方法計算，取較多的基函數(shù)項數(shù)可獲得較高的計算精度，實際計算中可根據(jù)所要求的計算精度來確定取多少項。

需要說明的是，文獻(xiàn)［8］的方法僅適用于半橢圓形狀，而本文方法不僅適用于半橢圓截面的頻率計算，同樣適用具有任意液深的橢圓截面容器。

表1 前三階正規(guī)化的反對稱以及正對稱晃動頻率Tab．1 The first three normalized anti-symmetric and symmetric sloshing frequencies Ωj＝g（j＝1，2，3）

表1 前三階正規(guī)化的反對稱以及正對稱晃動頻率Tab．1 The first three normalized anti-symmetric and symmetric sloshing frequencies Ωj＝g（j＝1，2，3）

無蓋板 （a′＝1．0 m，b′＝0．8 m）反對稱 正對稱本文 文獻(xiàn)［8］ 本文 文獻(xiàn)［8］Ω1 1．125 10 1．113 90 2．664 79 2．667 65 Ω2 4．132 64 4．120 74 5．562 31 5．578 62 Ω3 6．979 68 6．973 76 8．392 50 8．408 76有蓋板 （a′＝1．0 m，b′＝0．05 m，L＝0．05 m）反對稱 正對稱本文 文獻(xiàn)［8］ 本文 文獻(xiàn)［8］Ω1 0．028 2 0．028 2 0．101 2 0．101 6 Ω2 0．220 4 0．219 5 0．380 7 0．380 6 Ω3 0．588 9 0．582 5 0．828 7 0．822 4

圖6 半橢圓 a′＝1．0 m，b′＝0．1 m頻率收斂圖線Fig 6．Convergence curves of frequencies（half elliptical section a′＝1．0 m，b′＝0．1 m）

圖7 半橢圓 a′＝1．0 m，b′＝0．8 m頻率收斂圖線Fig．7 Convergence curves of frequencies（half elliptical section a′＝1．0 m，b′＝0．8）

5 結(jié) 論

基于Ritz方法，本文引入Green公式，將矩形截面晃動模態(tài)的精確解作為基函數(shù)，針對任意二維截面形狀，建立了特征值問題Ritz求解方法的統(tǒng)一計算格式，避免不同截面（由于邊界條件的不同）需要引入不同基函數(shù)的復(fù)雜性，使得Ritz方法在理論上適用于任意截面的二維晃動模態(tài)分析，可同時得到模態(tài)頻率與振型解答，采用本文方法對半頂角為45°的三角形、半橢圓形截面的晃動模態(tài)進(jìn)行求解，本文計算結(jié)果與其它解析方法的結(jié)果吻合良好，其計算精度完全滿足工程要求。

［1］Dodge F T．The New Dynamic Behavior of Liquids in Moving Containers［R］．Southwest Research Institute，San Antonio，TX，2000．

［2］Ibrahim R A． Liquid sloshing dynamics： theory and applications，［M］．Cambridge：Cambridge University Press，2005．

［3］Faltinsen O M，Timokha A N．Sloshing［M］．Cambridge：Cambridge University Press，2009．

［4］Lawrence H R，Wang C J，Reddy R B．Variational solution of fuel sloshing modes jet propulsion，1958，28（11）：728－736．

［5］Lukovsky L A．Variational methods of solving dynamic problems for fluid-containing bodies［J］． International Applied Mechanics，2004，40（10）：1092－1128．

［6］Budiansky B．Sloshing of liquids in circular canals and spherical tanks［J］．Journal of Aerospace Sciences，1960，27（3）：161－173．

［7］Fox D W，Kuttler J R．Sloshing frequencies，zeitschrift für angewandte mathematik und physik ZAMP，1983，34（5），668－696．

［8］Hasheminejad SM，Aghabeigi M．Liquid sloshing in half-full horizontal elliptical tanks［J］． Journal of Sound and Vibration，2009，324（1－2）：332－349．

［9］Mciver P．Sloshing frequencies for cylindrical and spherical containers filled to an arbitrary depth［J］．Journal of Fluid Mechanics，1989，201：243－257．

［10］周叮．截面任意形狀的柱形和環(huán)形容器內(nèi)液體晃動特性的精確解法［J］．工程力學(xué)，1995，12（2）：58－64．ZHOU Ding．An exact solution of liquid sloshing in cylindrical and annular tanks with arbitrary cross-sections［J］．Engineering Mechanics，1995，12（2）：58－64．

［11］Ibrahim R A，Pilipchuk V N，Ikeda T．Recent advances in liquid sloshing dynamics［J］．Applied Mechanics Reviews，2001，54（2）：133－199．

［12］Moiseev N N，Petrov A A．The calculation of free oscillations of a liquid in a motionless container［J］．Advances in Applied Mechanics，1966，9：91－154．

［13］Ohayon R，Morand J P．Fluid-structure interaction：applied numerical methods［M］． New York： John Wiley ＆Sons，1995．

［14］Li Y，Wang Z．An approximate analytical solution of sloshing frequencies for a liquid in various shape aqueducts［J］．Shock and Vibration，http：11dx．doi．org／10．1155／2014．672648