儲液容器內液體晃蕩的非線性動力學分析＊

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

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

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

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

1 非線性波動理論模型

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

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

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

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

式中：s為偏應力，ep為偏應變率。

在液體晃動問題中，重力作為恢復力是必不可少的，因此在初始狀態將重力賦于所有流體單元。晃動波高滿足：

液體晃動動態平衡方程：

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

式中：üM為節點相對加速度矢量；MNM為對角化的一致質量矩陣，MNM＝∫VρNNNMdV，NN為單元形函數矢量；IN＝∫V（βN：σ）dV，σ為單元應力；PN＝∫SNMfds＋∫VNNFdV，f為面力，F 為體力。

當儲液系統受到外加激勵¨Xg時，上式可寫為：

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

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

2 算 例

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

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

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

自由液面波動形狀：

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

2.1 數值驗證

為了證明本文中所述有限單元法的準確性，將水平正弦波激勵作用下二維矩形容器內液體晃動問題作為驗證實例。在底部平行x方向施加水平正弦波激勵¨ug＝Agsinωit，其中A為激勵幅值（實例中取為0.01）。對每次施加的加速度頻率ωi，獲得其前4個周期內對應的液體自由表面晃動波高的最大值ηmax，如圖2所示。由于所施加的外激勵頻率ωi與系統固有頻率相等時系統即發生共振現象，而流體的共振主要表現為自由液面晃動時產生最大波高，通過波高峰值對應的頻率可確定系統的固有頻率。所以分別找出前三階固有頻率：第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。

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

圖2 儲液系統晃動波高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倍。可見，儲液容器液體晃動響應主要由第1階固有頻率決定。圖3為第1階固有頻率和第3階固有頻率激勵作用時自由液面點B的波高曲線，可見，激勵共振頻率決定了液體晃動波形形狀。

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

2.2.1 第1階激振頻率時動力響應分析

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

為了研究激勵作用下不同液深處流體的晃動響應，在流體右界面（流體與剛壁交界面）由自由液面向下每隔0.5m依次取節點。圖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階晃動模態圖及對應液體晃動位移矢量圖，由圖可以看出，自適應網格技術使非線性流體動力響應分析保持很好的結構網格，整個計算過程中單元網格沒有出現過大的畸形，保證了使用該方法計算結果的穩定與準確。

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

2.2.2 不同幅值激勵時動力響應分析

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

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

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

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

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

3 結 論

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

［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］包光偉，王振強，張天翔.火箭推進劑液體晃動關機響應的數值仿真［J］.宇航學報，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］劉永濤，馬寧，顧解忡.各種激勵作用下液艙內液體晃蕩的計算與分析［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有限元模擬的網格更新方法［J］.力學學報，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］梅強中.水波動力學［M］.北京：科學出版社，1984.