賈善坡 ， 許成祥

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

1 引 言

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

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

2 基本方程

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

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

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

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

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

3 液體晃動(dòng)問題的有限元模型

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

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

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

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

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

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

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

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

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

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

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

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

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

就有：

由此得：

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

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

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

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

4 計(jì)算實(shí)例

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

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

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

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

4.1 自由晃動(dòng)分析

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

[2]許成祥,賈善坡.儲(chǔ)罐結(jié)構(gòu)有限元?jiǎng)恿δＰ托拚皡?shù)辨識(shí)研究[J].武漢理工大學(xué)學(xué)報(bào),2010,32(9):286-290.

[3]廖樂康,張圣坤.承船廂二維非線性晃動(dòng)的分析[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)動(dòng)力學(xué)[M].北京:科學(xué)出版社,2002.

[5]杜 穎,劉習(xí)軍,賈啟芬.液固耦合動(dòng)力學(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]李遇春,樓夢(mèng)麟.渡槽中流體非線性晃動(dòng)的邊界元模擬[J].地震工程與工程振動(dòng),2000,20(2):51-56.

[9]包光偉,王振強(qiáng),張?zhí)煜?等.火箭推進(jìn)劑液體晃動(dòng)關(guān)機(jī)響應(yīng)的數(shù)值仿真[J].宇航學(xué)報(bào),2002,23(2):84-88.

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

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

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

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

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

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