儲液容器內液體自由晃動的有限元分析

賈善坡 ， 許成祥

（1長江大學 城市建設學院，湖北 荊州 434023；2山東大學 巖土與結構工程研究中心，濟南 250061）

1 引 言

液體的晃動問題在船舶、航天、石化、高層建筑減振和城市供水等領域均有應用，具有廣泛的工程背景。在液固耦合系統動力學的研究中，液體晃動動力學是必須的基礎。容器液體晃動和控制問題的研究受國內外研究人員的廣泛重視，主要研究液體自由晃動的固有頻率，在激勵下液體的受迫晃動以及它對貯箱的作用力，貯箱與液體晃動之間的耦合動力學問題，液體晃動的等效力學模型，液體晃動研究的工程方法，液體大幅度晃動等課題[1-4]。到目前為止，僅求出了少數特殊形狀容器內液體晃動特性的解析解，對于一般形狀的容器，多采用數值方法求解。目前，出現一批數值求解液體晃動的成功方法，如有限差分方法、有限元法、邊界積分法和邊界元法等[5-7]。具有代表性的研究方法有：王照林等[4]針對液體非線性晃動的一些定性理論以及穩定性問題進行了較為系統的研究；李遇春等[8]利用邊界元方法對液體非線性晃動問題進行了數值分析；包光偉等[9]采用VOF方法對液體晃動進行了數值仿真；周宏等[10]采用任意拉格朗日—歐拉（簡稱ALE）方法對帶有自由液面的晃動問題進行了分析討論；尚春雨等[11]提出用FLUENT計算液體晃動問題的方法；劉富等[12]采用SPH方法對棱形液艙在外加激勵作用下，不同充液比工況所對應的艙內液體晃動進行了三維數值模擬，并將其與實驗進行對比，兩者吻合較好，同時成功地模擬出液體晃動產生的波浪翻卷和破碎；李青等[13]推導了任意三維貯箱內液體晃動的等效力學模型，采用有限元方法建立計算等效力學模型參數的數值算法，利用矩陣相似變換的性質分離液體零頻以消除液體剛度矩陣的奇異性，并采用商用分析軟件FLOW-3D驗證了等效力學模型的有效性。

容器內液體的晃動問題，可以歸結為兩種類型：一是特征值問題，求解容器內液體晃動的自然頻率，是進一步研究液體晃動問題的基礎。第二類型是研究在外界擾動下，容器內液體的反應。數值模擬液體的晃動實際上就是用數值方法求解帶有自由邊界的非定常流體的動力學問題，由于自由液面的位置是未知的，而且自由液面邊界條件是復雜的非線性方程，因此，具有自由液面的流體的流動是用數值方法求解的最困難的問題之一。當貯液結構（容器）具有大的空間剛度時，可忽略器壁的彈性變形，將容器近似為絕對剛體。本文建立了可壓縮液體自由晃動的特征方程，并提出了相應的有限元計算方法，計算過程中假定容器壁面是剛性的，液體在固—液界面上沿法線方向的分速度為零，在晃動自由面上忽略了液體的表面張力。

2 基本方程

設液體為無旋無粘性的，計及流體的可壓縮性，則液體在域V內滿足連續性方程、運動方程和狀態方程，即：

式中，上標“·”表示對時間的求導，vi是流體擾動的流體速度分量，ρ0是擾動前流體的密度，ρ和p分別是擾動引起的密度變化和流場壓力變化，c0是流體中的聲速，表示為：

式中k為流體的體積模量。

將容器壁簡化為剛性的，流體在容器濕表面?Vw（固壁邊界）上滿足不可滲透性條件，在自由面?Vf（自由面邊界）上滿足運動學邊界條件和動力學邊界條件，則三維容器內液體晃動特征問題滿足下列方程和邊界條件[1，4]：

式中，g是重力加速度，V、?Vf和?Vw分別是液體域、液體自由面邊界和腔壁邊界。

3 液體晃動問題的有限元模型

對流體采用壓力格式，則單元內的壓力分布可表示為：

同樣可以形式上將液體域V上的解函數寫成如下形式：

式中，（p1,p2,…,pn）是液體域V上的所有節點變量。

對（3）式利用Galerkin法形成求解方程，即：

對（6）式中的V∫δp（ p,ii）dV項進行分部積分，可得：

考慮到δp的任意性，將（5）式代入（7）式，則可得到液體自由晃動問題的有限元方程，

式中，p為流體節點壓力向量，Mf和Kf分別為流體質量矩陣和流體剛度矩陣，具體表示式為：

方程（8）的解假設為如下形式：

式中，z是n階向量，ω是向量z的振動頻率，t為時間變量，t0是由初始條件確定的時間常數。將（11）式代入方程（8），可以導出液體自由晃動問題的特征方程：

其中，固有頻率 ωj（j=1,2,…,n）對應于固有模態列陣zj。

對于三維液體晃動問題，系統自由度相當多，計算中耗時、耗費及求解困難。因此，有必要采用縮減技術來分析三維液體晃動問題，縮減技術可以減少系統的自由度，節省計算和分析的耗費[14-15]。將p寫為p=[pm, ps]，（8）式可寫為：

式中，下標m表示保留的主自由度，而s表示被縮減掉的副自由度。

假設副自由度質量表示的慣性力很小前提下，有：

就有：

由此得：

根據系統動能和勢能在縮減前后不變，有：

由此可導出縮減后質量、剛度矩陣為：

因此，系統的振動方程變為：

平衡方程組數目由原來的n降為所選擇的主自由度數m了。

4 計算實例

考慮一圓柱形的剛性含液容器，容器底部固定（見圖1）。主要幾何及力學參數為：油罐半徑R=40m，充液高度h=30m，流體密度ρ0=1 000kg/m3，重力加速度g=9.8m/s2，聲速 c0=1 435m/s（流體體積模量 2.06GPa）。

假定儲罐是剛性的，液體是可壓縮的，地面運動為水平平移運動，沒有旋轉分量。液體對流晃動模態是基于Laplace等式的第二項的線性解求得[1]：

圖1 圓柱形容器內液體晃動的有限元模型Fig.1 The finite element model of liquid sloshing in cylindrical container

式中， fi為液體晃動第 i階頻率（單位：rad/s）；λi為一階Bessel函數導數的第i個根；g為重力加速度；R為儲罐半徑；h為液體高度。

4.1 自由晃動分析

通過（22）式獲得的前幾階液體自然晃動頻率的理論解如表1所示。采用本文所介紹的有限元法計算液體的自然晃動頻率及液面晃動形式，兩種計算方法獲得的結果比較如表1，有限元所得的液面晃動形式如圖2所示。

在表1中，m和n分別代表橫截面內兩個相互垂直方向上的半波數，由表1的計算值和理論解比較，結果令人滿意，可見本文所提出的有限元法是十分有效和可靠的。

表1 剛性容器內流體自由晃動頻率計算結果比較（單位：rad/s）Tab.1 Results calculated by finite element method and theoretic solution

圖2 流體自由表面波動模態圖Fig.2 Modal characteristics of liquid sloshing in cylindrical container

從特征值計算得到的頻率值以及相應的模態分析可知，流體頻率實際上是由低頻和高頻兩部分組成的。低頻部分對應于流體自由表面波動，此時流體動壓力僅在自由表面附近不為零，在流體內部為零。高頻部分則對應于由于流體可壓縮性引起的內部動壓力波動，而在自由表面處動壓力為零。如果在計算中不考慮流體自由表面波動，則特征值計算的結果將只能得到高頻部分。如果要想獲得更多階的低頻，則需將橫截面內的網格加密。另外，還對不同壓縮性（即不同聲速）時的情況進行了計算，結果發現流體壓縮性對自由表面波動部分的頻率基本無影響，而對高頻部分（對應于內部壓力波動）的影響則很大，而且頻率值與聲速近似成正比關系，這意味著當考慮流體為不可壓縮（即聲速無窮大）時，這部分頻率也將為無窮大。

4.2 晃動頻率與容器半經的關系

為了揭示液體晃動頻率與容器半經的關系，保持充液深度h=30m不變，改變容器半經，得出了前幾階自然晃動頻率，以前4階固有頻率為例，畫出f-R關系曲線，如圖3所示，研究液體晃動頻率與容器半徑之間的關系，可以看出，隨著容器半徑的增大，液體晃動頻率逐漸減小，液體晃動頻率與容器半徑基本上成反比例關系。

圖3 液體晃動頻率與容器半徑的關系Fig.3 Liquid sloshing frequencies versus the radius of container

圖4 液體晃動頻率與液體深度的關系Fig.4 Liquid sloshing frequencies versus the depth of liquid

4.3 晃動頻率與充液深度的關系

為了揭示液體晃動頻率與充液深度的關系，保持容器半經R=40m不變，改變充液深度，得出了前幾階自然晃動的頻率。圖4表示出半徑不變的容器其內部液體自然晃動的頻率，隨著液面高度變化而發生變化的趨勢，橫坐標表示深度h，縱坐標表示各階頻率值。計算結果表明，容器內液體的充滿程度即液面的高度h對液體晃動頻率也有很大的影響，液體深度對晃動頻率的影響很明顯。總體上說，對于第1階頻率來說，當充液深度達到40m時，頻率達到最大值，隨后，隨著充液深度的增大，而頻率逐漸減小；對第2階頻率來說，當充液深度達到20m時，頻率達到最大值，隨后，隨著充液深度的增大，而頻率逐漸減小；對第3、4階頻率來說，隨著充液深度的增大而頻率逐漸減小。

5 結 論

本文采用有限元法建立了圓柱形容器內液體晃動問題的數值計算方法，假定容器壁面是剛性的，在晃動自由面上忽略了液體的表面張力，采用縮減法來分析容器內液體的晃動動力特性問題，該方法可適用于復雜形狀容器內液體的晃動問題。由計算實例的結果可知，應用有限元法計算能夠求得容器內液體自然晃動的各階頻率和振型，且具有較高的精確度，該方法可以應用于三維任意形狀容器內液體的晃動問題求解。在設計中如何考慮液體晃動的影響以及如何抑制液體晃動，仍是一個值得研究的課題，此外，將液體晃動問題與結構彈性體振動問題統一起來就可以應用有限元法求解液體-結構耦合問題。

[1]賈善坡.大型儲罐地基變形特性研究及儲液晃動分析[D].東營:中國石油大學(華東),2006.

[2]許成祥,賈善坡.儲罐結構有限元動力模型修正及參數辨識研究[J].武漢理工大學學報,2010,32(9):286-290.

[3]廖樂康,張圣坤.承船廂二維非線性晃動的分析[J].船舶力學,2006,10(2):47-55.Liao Lekang,Zhang Shengkun.Analysis on 2-D nonlinear slosh of ship-chamber[J].Journal of Ship Mechanics,2006,10(2):47-55.

[4]王照林,劉延柱.充液系統動力學[M].北京:科學出版社,2002.

[5]杜 穎,劉習軍,賈啟芬.液固耦合動力學問題的研究[J].機床與液壓,2004(11):9-12.

[6]Liu Dongming,Lin Pengzhi.A numerical study of three-dimensional liquid sloshing in tanks[J].Journal of Computational Physics,2008(227):3921-3939.

[7]Mitra S,Sinhamahapatra K P.2D simulation of fluid-structure interaction using finite element method[J].Finite Elements in Analysis and Design,2008(45):52-59.

[8]李遇春,樓夢麟.渡槽中流體非線性晃動的邊界元模擬[J].地震工程與工程振動,2000,20(2):51-56.

[9]包光偉,王振強,張天翔,等.火箭推進劑液體晃動關機響應的數值仿真[J].宇航學報,2002,23(2):84-88.

[10]周 宏,李俊峰,王天舒.用ALE有限元模擬的網格更新方法[J].力學學報,2008,40(2):267-272.

[11]尚春雨,趙金城.用FLUENT分析剛性容器內液面晃動問題[J].上海交通大學學報,2008,42(6):953-956.

[12]劉 富,童明波,陳建平.基于SPH方法的三維液體晃動數值模擬[J].南京航空航天大學學報,2010,42(1):122-126.

[13]李 青,馬興瑞,王天舒.非軸對稱貯箱液體晃動的等效力學模型[J].宇航學報,2011,32(2):242-249.

[14]趙玉成,張亞紅,白長青,許慶余.固體火箭發動機模態分析的縮聚方法[J].固體火箭技術,2004,27(2):98-100.

[15]朱安文,曲廣吉,高耀南.航天器結構動力模型修正中的縮聚方法[J].中國空間科學技術,2003(2):6-10.