二維晃動模態的統一Ritz計算格式

李遇春，張 皓

（同濟大學 水利工程系，上海 200092）

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

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

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

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

1 二維自由晃動基本方程

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

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

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

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

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

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

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

根據圖1的液體區域Ω，引入平面Green定理：

式中：S＝Sw＋Sf，P與Q為封閉區域Ω內連續可微的函數，若定義：

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

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

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

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

2 Ritz求解方法

式（15）可以改寫為：

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

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

＝（φ1，…，φm）為基函數向量

將表達式（17）代入式（16）可以得到：

如果記：

則式（18）可以寫為：

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

3 Ritz基函數的構造

3.1 矩形截面基函數及其模態解

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

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

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

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

式中：

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

3.2 任意截面基函數及其模態解

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

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

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

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

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

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

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

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

其中

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

4 數值算例

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

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

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

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

4．2 半橢圓截面

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

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

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

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

為了考察此時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分別顯示了二個半橢圓容器內流體前3階反對稱與對稱晃動頻率隨基函數項數m變化的圖線，由圖可見，對于本算例而言，Ritz方法收斂很快，只需取5～6項即可獲得較高的計算精度，可滿足工程使用要求。

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

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

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

表1 前三階正規化的反對稱以及正對稱晃動頻率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）反對稱 正對稱本文 文獻［8］ 本文 文獻［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）反對稱 正對稱本文 文獻［8］ 本文 文獻［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 結 論

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

［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］周叮．截面任意形狀的柱形和環形容器內液體晃動特性的精確解法［J］．工程力學，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