儲液容器內(nèi)液體晃蕩的非線性動力學(xué)分析＊

李文盛，趙友清，賈善坡，2，3，王 凱，譚繼可

（1.長江大學(xué)城市建設(shè)學(xué)院，湖北 荊州434023；2.山東大學(xué)巖土與結(jié)構(gòu)工程研究中心，山東 濟(jì)南250061；3.中國科學(xué)院武漢巖土力學(xué)研究所巖土力學(xué)與工程國家重點實驗室，湖北 武漢430071）

液體晃蕩問題廣泛存在石油、核工程、船舶、化工等領(lǐng)域，具有非常復(fù)雜的流體運動現(xiàn)象，呈現(xiàn)很強的非線性和隨機(jī)性。當(dāng)外界激勵頻率接近儲液系統(tǒng)的固有頻率時，內(nèi)部儲存液體晃動將變得更劇烈并導(dǎo)致容器破壞，不僅造成重大經(jīng)濟(jì)損失及人員傷亡，還造成更嚴(yán)重的次生災(zāi)害。因此，對液體晃蕩和控制晃蕩的研究引起了廣泛關(guān)注［1－2］。

目前，對液體晃蕩問題的數(shù)值研究方法主要有MAC法、VOF法、有限單元法（混合插值FEM、分布FEM、ALE有限元方法）及邊界元方法等，并取得了一些重要成果。A.A.Amsden等［3］采用MAC方法研究了帶有自由液面的大幅晃動問題，但這種方法計算存儲量大、穩(wěn)定性差，V.Armenio［4］提出了一種改進(jìn)的 MAC法；包光偉等［5］、J.H.Jung等［6］采用VOF方法對液體晃動問題進(jìn)行了數(shù)值計算；劉永濤等［7］基于VOF法得到液體晃蕩運動的數(shù)值計算模型；周宏等［8］利用任意拉格朗日歐拉法（ALE）實現(xiàn)了自由液面的跟蹤，并對晃動問題進(jìn)行了分析討論；R.Sygulski［9］假設(shè)流體為無黏、不可壓縮、微幅晃動的理想流體，采用邊界元方法研究了三維液體的晃動問題；C.Z.Wang等［10］基于流體速度勢理論對流體非線性晃蕩問題進(jìn)行了有限元分析。

1 非線性波動理論模型

使用流體動力學(xué)非線性波動理論求解儲液系統(tǒng)流體自由液面大幅波動問題，將三維模型簡化為二維幾何非線性平面應(yīng)變問題，本文中運用有限元軟件ABAQUS［11］，采用狀態(tài)方程進(jìn)行求解。在非線性波動理論模型有限元網(wǎng)格中，流體部分采用4節(jié)點平面應(yīng)變單元CPE4R，每個節(jié)點具有2個位移自由度；剛體部分采用R2D2單元；剛體與流體相互作用使用法向接觸，其中剛體表面為主面，流體表面為從面；儲液容器壁高足夠大以保證晃動時液體沒有溢出。

為了得到穩(wěn)定的數(shù)值解，假設(shè)流體材料為幾乎不可壓縮的，且黏性很小。流體體積響應(yīng)通過us－up線性狀態(tài)方程和牛頓黏性剪切模型建模。在儲液系統(tǒng)模型中，為確定材料體積強度和壓力提供了一個水動力模型作為液體密度和能量（質(zhì)量內(nèi)能）的函數(shù)。Mie－Grüneisen狀態(tài)方程中，利用激波速度us和粒子速度up的線性關(guān)系，有：

式中：ζ＝1－ρ0／ρ，ρ0是擾動前流體密度，ρ是擾動引起變化的流體密度。

模型中，非線性通過體積應(yīng)變計算，由于水的體積彈性模量很大（2.07GPa），幾乎不可壓縮，根據(jù)文獻(xiàn)［11］，取體積彈性模量比真實值小2～3個數(shù)量級，即2.07MPa，也滿足了介質(zhì)的不可壓縮性。水的剪切黏度ν應(yīng)足夠小，因為水無黏性，一個大的黏度將導(dǎo)致過硬反應(yīng)。由牛頓流體Navier－Poisson定律，有：

式中：s為偏應(yīng)力，ep為偏應(yīng)變率。

在液體晃動問題中，重力作為恢復(fù)力是必不可少的，因此在初始狀態(tài)將重力賦于所有流體單元?；蝿硬ǜ邼M足：

液體晃動動態(tài)平衡方程：

羅四強說著拉起阿里的手，帶著他朝花壇走去。剛走幾步，另一側(cè)的悼念廳又有哀樂響起。阿里又是一個寒噤，大聲道：“姆媽醒了!”說話間，甩開羅四強的手，又朝著側(cè)廳奔跑。羅四強急追好幾步，才抓到他。

式中：üM為節(jié)點相對加速度矢量；MNM為對角化的一致質(zhì)量矩陣，MNM＝∫VρNNNMdV，NN為單元形函數(shù)矢量；IN＝∫V（βN：σ）dV，σ為單元應(yīng)力；PN＝∫SNMfds＋∫VNNFdV，f為面力，F(xiàn) 為體力。

當(dāng)儲液系統(tǒng)受到外加激勵¨Xg時，上式可寫為：

利用體積黏性壓力抑制最高頻率振蕩單元，儲液系統(tǒng)體積黏性線性表達(dá)式為［12］：

式中：b1為阻尼系數(shù)；cd為疏密波的速度；Le為單元特征長度；在方程（5）時，體積黏力也應(yīng)包含在內(nèi)。

2 算 例

模型為剛性矩形容器，底部固定于剛性基礎(chǔ)。模型相關(guān)參數(shù)如下：底部截面寬度2a＝3.0m，剛壁高度H＝4.5m，容器內(nèi)液體深h＝4.0m，密度ρ＝1 000kg／m3，聲速c0＝1 438.75m／s（體積模量2.07GPa），黏滯系數(shù)0.001Pa·s，重力加速度g＝9.81m／s2，如圖1所示。

圖1 非線性波動理論模型Fig.1 The 2Dmodel considering nonlinear wave theory

二維儲液容器自由晃動問題中，濕面在x、y方向是相互獨立的，對Navier－Strokes方程分離變量即可獲得流體自由晃動的固有頻率。由于容器內(nèi)流體材料為水，黏性足夠小，可視為不可壓縮、無黏、無旋的理想流體，則流體晃動的固有頻率［13］：

自由液面波動形狀：

式（8）～（9）表明，液體自由液面波動形狀為正弦或余弦曲線。

2.1 數(shù)值驗證

為了證明本文中所述有限單元法的準(zhǔn)確性，將水平正弦波激勵作用下二維矩形容器內(nèi)液體晃動問題作為驗證實例。在底部平行x方向施加水平正弦波激勵¨ug＝Agsinωit，其中A為激勵幅值（實例中取為0.01）。對每次施加的加速度頻率ωi，獲得其前4個周期內(nèi)對應(yīng)的液體自由表面晃動波高的最大值ηmax，如圖2所示。由于所施加的外激勵頻率ωi與系統(tǒng)固有頻率相等時系統(tǒng)即發(fā)生共振現(xiàn)象，而流體的共振主要表現(xiàn)為自由液面晃動時產(chǎn)生最大波高，通過波高峰值對應(yīng)的頻率可確定系統(tǒng)的固有頻率。所以分別找出前三階固有頻率：第1階，固有頻率ω1＝3.20rad／s，η1，max＝0.201m；第2階，固有頻率ω2＝4.51rad／s，η2，max＝0.025m；第3階，固有頻率ω3＝5.60rad／s，η3，max＝0.029m。

為驗證上述所得結(jié)果，通過式（7）計算前3階液體自由晃動時固有頻率解析解（見表1）。由表1可知，本文中采用的有限單元法具有有效性和可靠性。

圖2 儲液系統(tǒng)晃動波高Fig.2 Surface wave displacement for the tank－liquid system

表1 液體自由晃動頻率Table 1 Frequencies of liquid sloshing

2.2 水平正弦激勵下的液體晃蕩

由圖2可知，在不同頻率的正弦激勵作用下引起的液體晃動最大波高，其中第1階固有頻率引起的液體晃動最大波高是第2、3階固有頻率的6～8倍?？梢?，儲液容器液體晃動響應(yīng)主要由第1階固有頻率決定。圖3為第1階固有頻率和第3階固有頻率激勵作用時自由液面點B的波高曲線，可見，激勵共振頻率決定了液體晃動波形形狀。

圖3 不同諧頻下自由液面點B的波面響應(yīng)Fig.3 Free surface elevation of poit B for harmonic frequencies

2.2.1 第1階激振頻率時動力響應(yīng)分析

算例中，A＝0.01，第1階的固有頻率ω1＝3.20rad／s，總時間為10s，分析步長Δt＝0.01s，水的黏滯系數(shù)足夠小，取ν＝0.001。圖4（a）為點A、O和B在第1階激振頻率液體晃蕩時自由液面波高曲線，A、B處液面波動隨著周期數(shù)增多而逐漸加強，意味著系統(tǒng)在該頻率作用下發(fā)生了明顯的低頻共振現(xiàn)象，施加激勵的頻率即為系統(tǒng)固有頻率值。液面中心點O處液體隨時間變化幾乎沒有明顯波動；左壁面點A液體晃動波高在t＝8.82s時出現(xiàn)最大值ηmax＝0.183m；右壁面點B液體晃動波高在t＝9.80s時出現(xiàn)最大值ηmax＝0.201m。觀察點A、B處的波高曲線，兩點的波動性態(tài)具有相反性，且液體晃動具有明顯的非線性，主要表現(xiàn)在液體晃動的幅值向上的要大于向下的。

為了研究激勵作用下不同液深處流體的晃動響應(yīng)，在流體右界面（流體與剛壁交界面）由自由液面向下每隔0.5m依次取節(jié)點。圖4（b）給出計算獲得的點B、C、D、E、F、G液體晃動波高曲線，由圖可以看出，波高幅值隨液體深度的增加逐漸遞減，這說明了水平正弦激勵下底部固定的儲液容器的振動是基本的梁型振動。

圖4 第1階頻率作用下液體晃動波高曲線Fig.4 Variation in time of the surface wave in the first sloshing mode

圖5給出了非線性分析時流體第1階晃動模態(tài)圖及對應(yīng)液體晃動位移矢量圖，由圖可以看出，自適應(yīng)網(wǎng)格技術(shù)使非線性流體動力響應(yīng)分析保持很好的結(jié)構(gòu)網(wǎng)格，整個計算過程中單元網(wǎng)格沒有出現(xiàn)過大的畸形，保證了使用該方法計算結(jié)果的穩(wěn)定與準(zhǔn)確。

圖5 第1階液體晃動模態(tài)圖和液體晃動位移矢量圖Fig.5 Liquid shapes corresponding to the first sloshing mode and displacement vector plot

2.2.2 不同幅值激勵時動力響應(yīng)分析

取水平正弦激勵¨ug＝0.01gsinω1t和¨ug＝0.005gsinω1t，圖6分別給出了兩種不同幅值激勵作用下自由液面點A、B處液體晃動響應(yīng)曲線。由圖可知，液體晃動波高曲線與水平正弦激勵頻率基本相同，說明在不同幅值激勵作用下液體晃動響應(yīng)頻率不隨幅值變化；同時液體晃動上下波動的幅值出現(xiàn)明顯不對稱性，向上的波動高于向下的波動，并且相差程度隨著幅值的增加不斷增大。這些現(xiàn)象恰恰證明了液體晃動具有強的非線性。

圖6 不同幅值激勵下液體晃動的波高曲線Fig.6 Free surface elevation of liquid for different amplitudes

2.3 EI地震波作用下的液體晃蕩

采用地震反應(yīng)分析最常用的EI Centro地震波，最大峰值加速度為35cm／s2，加速度時間間隔為0.01s，總時長20s，圖7（a）為該地震波加速度曲線，圖7（b）為EI Centro地震波作用下液體自由液面相對兩點A、B的晃動波高曲線。從圖可以看出，晃動時點A處的波高在t＝17.80s時出現(xiàn)最大值ηmax＝9.3cm；晃動時點B處的波高在t＝15.10s時出現(xiàn)最大值ηmax＝6.9cm；在開始晃動的前12s，液面晃動波高幅值較小，隨后液面晃動波高幅值逐漸變大并在17.80s后減小，可見，地震波對儲液容器液體晃動反應(yīng)需要一個激發(fā)過程，液體晃動波高幅值并不同隨地震加速度幅值的變大而變大、減小而減小，最大液面波高幅值往往發(fā)生在激勵加速度最大幅值以后，這是因為輸入的地震波激勵加速度雖減小了，但儲液系統(tǒng)中輸入的地震波能量還在增加。同時還可以看出，點A、B的波動性態(tài)具有相反性，且液體晃動具有明顯的非線性，主要表現(xiàn)在液體晃動的幅值向上的大于向下的。

圖7 EI Centro地震波和液面點A、B液體晃動的波高曲線Fig.7 EI Centro seismic wave and free surface displacement curve at points Aand B

3 結(jié) 論

使用有限元法，基于非線性波動理論建立儲液容器液體晃動數(shù)學(xué)模型，通過對其施加水平簡諧激勵得到液體晃動的固有頻率和模態(tài)，研究了二維剛性儲液系統(tǒng)非線性晃動響應(yīng)特性。該方法亦可適用于任意復(fù)雜形狀三維彈性容器內(nèi)液體的晃蕩問題。由數(shù)值實例獲得的結(jié)果表明，該方法求得容器內(nèi)液體晃動的固有頻率具有較高的精確度，外加激勵頻率及幅值對液體晃蕩影響較大，表現(xiàn)了明顯的非線性。

［1］Hasheminejad S M，Aghabeigi M.Sloshing characteristics in half－full horizontal elliptical tanks with vertical baffles［J］.Applied Mathematical Modelling，2012，36（1）：57－71.

［2］Xue Mi－an，Lin Peng－zhi.Numerical study of ring baffle effects on reducing violent liquid sloshing［J］.Computers＆ Fluids，2011，52：116－129.

［3］Amsden A A，Harlow F H.A simplified MAC technique for incompressible fluid flow calculations［J］.Journal of Computational Physics，1970，6（2）：322－325.

［4］Armenio V.An improved MAC method（SIMAC）for unsteady high－reynolds free surface flows［J］.International Journal for Numerical Methods in Fluids，1997，24（2）：185－214.

［5］包光偉，王振強，張?zhí)煜?火箭推進(jìn)劑液體晃動關(guān)機(jī)響應(yīng)的數(shù)值仿真［J］.宇航學(xué)報，2002，23（2）：84－88.Bao Guang－wei，Wang Zhen－qiang，Zhang Tian－xiang.Numerical simulation of slosh of the propellant fuel during the period the shut－down of the rocket［J］.Jourual of Astronautics，2002，23（2）：84－88.

［6］Jung J H，Yoon H S，Lee C Y，et al.Effect of the vertical baffle height on the liquid sloshing in a three－dimensional rectangular tank［J］.Ocean Engineering，2012，44：79－89.

［7］劉永濤，馬寧，顧解忡.各種激勵作用下液艙內(nèi)液體晃蕩的計算與分析［J］.船海工程，2009，38（5）：7－12.Liu Yong－tao，Ma Ning，Gu Xie－chong.Calculation and analysis of liquid sloshing loads in tanks under different kinds of stimulations［J］.Ship and Ocean Engineering，2009，38（5）：7－12.

［8］周宏，李俊峰，王天舒.用 ALE有限元模擬的網(wǎng)格更新方法［J］.力學(xué)學(xué)報，2008，40（2）：267－272.Zhou Hong，Li Jun－feng，Wang Tian－shu.Mesh update algorithm in ale finite method within free surface flow［J］.Chinese Journal of Theoretical and Applied Mechanics，2008，40（2）：267－272.

［9］Sygulski R.Boundary element analysis of liquid sloshing in baffled tanks［J］.Engineering Analysis with Boundary Elements，2011，35（8）：978－983.

［10］Wang C Z，Khoo B C.Finite element analysis of two－dimensional nonlinear sloshing problems in random excitations［J］.Ocean Engineering，2005，32（2）：107－133.

［11］ABAQUS／explicit user’s manual：Version 6.4［M］.Rhode Island：Hibbit，Karlsson and Sorensen Inc，2002.

［12］ABAQUS／theory user’s manual：Version 6.4［M］.Rhode Island：Hibbit，Karlsson and Sorensen Inc，2002.

［13］梅強中.水波動力學(xué)［M］.北京：科學(xué)出版社，1984.