二維晃動(dòng)模態(tài)的統(tǒng)一Ritz計(jì)算格式

李遇春，張 皓

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

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

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

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

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

1 二維自由晃動(dòng)基本方程

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

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

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

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

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

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

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

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

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

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

式中：s表示曲線坐標(biāo)，在邊界S上沿逆時(shí)針?lè)较驗(yàn)檎＆觴＝cos（τ，x）及 τz＝cos（τ，z）為切向量 τ（見(jiàn)圖 1）的方向余弦，切向量與外法向量的方向余弦有如下的關(guān)系：

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

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

2 Ritz求解方法

式（15）可以改寫(xiě)為：

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

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

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

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

如果記：

則式（18）可以寫(xiě)為：

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

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

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

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

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

引入無(wú)量綱變量 x＝X／a0，z＝Z／a0，無(wú)量綱液深h＝H／a0，構(gòu)造基函數(shù) φj（x，z），可將其分離為兩個(gè)獨(dú)立函數(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），為特征值問(wèn)題的精確解析解。其第j階自由表面的模態(tài)函數(shù)分別為（反對(duì)稱(chēng)模態(tài)）及 cosκjx（正對(duì)稱(chēng)模態(tài)）。

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

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

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

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

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

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

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

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

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

其中

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

4 數(shù)值算例

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

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

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

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

4．2 半橢圓截面

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

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

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

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

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

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

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

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

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

無(wú)蓋板 （a′＝1．0 m，b′＝0．8 m）反對(duì)稱(chēng) 正對(duì)稱(chēng)本文 文獻(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）反對(duì)稱(chēng) 正對(duì)稱(chēng)本文 文獻(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公式，將矩形截面晃動(dòng)模態(tài)的精確解作為基函數(shù)，針對(duì)任意二維截面形狀，建立了特征值問(wèn)題Ritz求解方法的統(tǒng)一計(jì)算格式，避免不同截面（由于邊界條件的不同）需要引入不同基函數(shù)的復(fù)雜性，使得Ritz方法在理論上適用于任意截面的二維晃動(dòng)模態(tài)分析，可同時(shí)得到模態(tài)頻率與振型解答，采用本文方法對(duì)半頂角為45°的三角形、半橢圓形截面的晃動(dòng)模態(tài)進(jìn)行求解，本文計(jì)算結(jié)果與其它解析方法的結(jié)果吻合良好，其計(jì)算精度完全滿足工程要求。

［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)液體晃動(dòng)特性的精確解法［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