儲液容器內(nèi)液體自由晃動的有限元分析

賈善坡 ， 許成祥

（1長江大學(xué) 城市建設(shè)學(xué)院，湖北 荊州 434023；2山東大學(xué) 巖土與結(jié)構(gòu)工程研究中心，濟(jì)南 250061）

1 引 言

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

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

2 基本方程

設(shè)液體為無旋無粘性的，計及流體的可壓縮性，則液體在域V內(nèi)滿足連續(xù)性方程、運動方程和狀態(tài)方程，即：

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

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

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

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

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

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

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

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

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

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

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

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

方程（8）的解假設(shè)為如下形式：

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

其中，固有頻率 ωj（j=1,2,…,n）對應(yīng)于固有模態(tài)列陣zj。

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

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

假設(shè)副自由度質(zhì)量表示的慣性力很小前提下，有：

就有：

由此得：

根據(jù)系統(tǒng)動能和勢能在縮減前后不變，有：

由此可導(dǎo)出縮減后質(zhì)量、剛度矩陣為：

因此，系統(tǒng)的振動方程變?yōu)椋?/p>

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

4 計算實例

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

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

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

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

4.1 自由晃動分析

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

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

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

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

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

4.2 晃動頻率與容器半經(jīng)的關(guān)系

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

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

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

4.3 晃動頻率與充液深度的關(guān)系

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

5 結(jié) 論

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

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

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

[3]廖樂康,張圣坤.承船廂二維非線性晃動的分析[J].船舶力學(xué),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]王照林,劉延柱.充液系統(tǒng)動力學(xué)[M].北京:科學(xué)出版社,2002.

[5]杜 穎,劉習(xí)軍,賈啟芬.液固耦合動力學(xué)問題的研究[J].機(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]包光偉,王振強(qiáng),張?zhí)煜?等.火箭推進(jìn)劑液體晃動關(guān)機(jī)響應(yīng)的數(shù)值仿真[J].宇航學(xué)報,2002,23(2):84-88.

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

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

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

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

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

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