儲液容器內(nèi)液體晃動問題的動力學建模與模擬

賈善坡，趙友清，許成祥

(長江大學 城市建設學院，湖北 荊州434023)

0 引言

目前，充液系統(tǒng)液體晃動問題研究進一步熱化，主要表現(xiàn)在:微幅晃動問題研究，非線性大幅晃動問題研究，液固耦合晃動及考慮彈性壁的液固耦合晃動問題研究，液體晃動理論及等效力學模型研究等．液體晃動問題的解法主要有解析法和數(shù)值分析法． 解析法主要是采用Bessel 方程、Fourier 方程及Laplace 變換求解特殊形狀、邊界的線性流固耦合問題;數(shù)值計算方法同時將流體域和固體域離散化求解，如有限單元法、有限差分法、邊界元法和邊界積分法等．而對于復雜的三維晃動問題求解多使用有限單元法．具體代表如下:Housner［1］假設儲液罐壁是剛性的且置于剛性基礎上，使用Laplace 方程求解矩形和圓柱形儲液罐的基本模態(tài)． Wang Wei 等［2］使用FEM 分析了儲罐在不同邊界條件下液體晃動模態(tài)．賈善坡、鄭建華等［3-4］利用有限元軟件ANSYS 分析了固液耦合系統(tǒng)的模態(tài)及動態(tài)響應． 趙曉磊等［5］利用有限元軟件ADINA 對大型儲液罐進行了模態(tài)分析．包光偉等［6］采用VOF 方法對儲箱內(nèi)液體晃動進行了數(shù)值模擬分析． 尚春雨等［7］利用FLUENT 方法分析了液體晃動問題．劉富等［8］采用SPH 方法成功模擬了三維液體晃動產(chǎn)生的波浪翻卷． 李青等［9］推導出任意三維容器內(nèi)液體晃動的等效力學模型，并應用FLOW-3D 分析軟件驗證了該模型的正確性．

筆者使用聲學單元描述流體，流體自由表面單元“削皮”成膜單元，兩者通過彈簧單元連接，實現(xiàn)了位移與應力的兼容，并將計算得到的晃動模態(tài)及頻率與解析解對比，結(jié)果表明該方法是可行的，最后使用該方法分析了矩形儲液罐在不同液體高寬比時，其晃動頻率與液體深度、寬度的關系，并從中獲得矩形儲液罐液體晃動的一般規(guī)律．

1 液體晃動控制方程

1．1 聲學方程

流體視為聲學介質(zhì)，將其作為彈性介質(zhì)分析其晃動動力特性，對于可壓縮、無黏性和小擾動流體，考慮其材料阻尼的微幅運動平衡方程為［10］

式中:p 是流體動壓;x 是流體質(zhì)點的空間坐標;u·f是流體質(zhì)點的速度;u¨f是流體質(zhì)點加速度;ρf是流體的密度;γ 是體積曳力．

假設流體介質(zhì)是可壓縮、無黏性和線性的，計算流體動壓時考慮體積模量影響，其本構(gòu)方程為:

式中:Kf是流體的體積模量．

對可壓縮、無黏性流體，式(2)可寫為只含體積模量和體積應變線性方程:

p= -KfεV． (3)

式中:εV是流體體積應變．

1．2 流體動力學方程

對無黏性、可壓縮和小擾動的流體，由流體域Vf內(nèi)場方程聯(lián)立消去密度變量ρ 和流動速度vi得到以壓力擾動p 為場變量的波動方程:

式中:c0為流體中的聲速．

假設儲液罐壁為剛性且儲液罐放置在剛性基礎上，則可得到剛性邊界條件方程:

式中:n 為流體邊界外法線向量．

流體自由液面邊界條件方程:

2 有限元模型與數(shù)值方法

對流體域采用流場壓力p 作為基本變量，構(gòu)造插值函數(shù)Nk(x，y，z)，則流體域壓力分布為

式中:n 為流體單元節(jié)點總數(shù)．

對波動方程(4)和邊界條件方程(5)、(6)利用加權(quán)余量的Galerkin 法解方程:

式中:Vf為流體域;sb為剛性邊界面;sf自由液面邊界．

將式(7)代入式(8)，并考慮δp 的任意性，得到流體自由晃動有限元方程:+Kfp=0． (9)

式中:Mf為流體質(zhì)量矩陣;Kf流體剛度矩陣．

假定方程(9)解的形式p=φsinω(t-t0)并反代入，得到液體自由晃動的特征方程:

式中:ω1，ω2，…，ωn為n 個固有頻率;φ1，φ2，…，φn為其對應的n 個固有振型．

考慮到液體晃動系統(tǒng)自由度多且需要獲得大量固有模態(tài)，為提高其計算效率，采用蘭索斯法提取特征值．Lanczos 程序運行一系列迭代步，對每一次運行經(jīng)過譜變換:

［M］(K-σ［M］)-1［M］{Φ}=θ［M］{Φ}． (11)

式中:σ 為偏移值;θ 為特征值;{Φ}為特征向量．

方程(11)將快速收斂到所需要的特征值．方程(10)與方程(11)中，特征向量φ 與{Φ}是一致的，但固有頻率值ω 與θ 具有如下關系:

程序每運行一次Lanczos 命令生成一個Krylov 子空間序列，經(jīng)過一系列運算步計算出最近似的子空間特征向量． 理論上，Lanczos 的基本流程可以確定簡單的特征值，然而對大型矩陣特征值問題顯得太昂貴． 因此，筆者采用分塊Lanczos法［12］，通過創(chuàng)建一個正交向量塊，并利用每次Lanczos 步中的塊的大小增加Krylov 子空間的維數(shù)，這樣不僅可以自動計算大型矩陣特征值，還大大提高了計算效率．

3 數(shù)值算例

分析模型為一正方形截面的剛性儲液罐，其底部固定，邊長B =2 m，液面深度H =3 m，流體密度ρ=1 000 kg/m3，聲速c0=1 435 m/s，液體的體積模量2．06 GPa，重力加速度g=9．8 m/s2．

建立三維儲液容器液體晃動模型． 假設液體是無黏性、不可壓縮、無旋的理想液體，罐體錨固于剛性地基上，液體晃動為小波動．對液體域劃分為下部流體及其頂部自由液面，其中流體部分采用聲學單元AC3D8，頂部自由液面采用膜單元M3D4．在液體晃動問題中，重力作為恢復力是必不可少的，然而在聲學單元中不允許直接定義均布重力荷載，為了解決這個問題，引入重力彈簧附加至流體自由液面各個節(jié)點． 聲學單元僅有應力自由度，而彈簧單元僅有位移自由度，為了克服這種不兼容性，流體自由液面“削皮”成膜單元，膜單元沒有彎曲剛度，通過“TIE”命令實現(xiàn)與聲學單元的耦合．當液體晃動時，自由液面的每個微面均受到一個重力恢復力ρgA，A 為自由液面微面面積;其中取ρgA 的大小作為彈簧剛度(本例中流體橫截面被225 個節(jié)點分割，)． 另外，考慮到網(wǎng)格劃分對結(jié)果的影響，據(jù)文獻［11］知:空間步距Δx 需控制在Δx ＜λ/6 時，才能滿足波動數(shù)值模擬的精度要求，即一個波長內(nèi)一般不少于6 ～8個節(jié)點．

3．1 自由晃動結(jié)果分析

使用上述方法計算流體晃動頻率結(jié)果如圖1和圖2 所示．從圖中分析得到，流體自由晃動頻率包括低頻和高頻兩部分． 低頻部分對應于流體自由晃動時分布在自由液面的動壓力引起的自由液面波動;高頻部分對應于流體內(nèi)部動壓力波動．由于本例劃分網(wǎng)格時，在液體自由表面共布置了225 個節(jié)點，故分析中模態(tài)階數(shù)取其前300 階，其中前225 階為晃動頻率低頻部分，最大值為第225階ω=17．727 rad/s;晃動頻率高頻部分從第226 階ω=785．20 rad/s 到第300 階ω=9 331．70 rad/s．

圖1 流體自由晃動的低頻部分(1 ～70 階)Fig．1 Low-frequencies part of liquid sloshing(1 ～70)

為了驗證上述方法的準確性，首先計算出前9 階液體晃動頻率理論值，其次采用本文方法獲得液體晃動前9 階模態(tài)及頻率，最后與文獻［3］結(jié)果三者對比，如表1 所示．

圖2 流體自由晃動的高頻部分(226 ～300 階)Fig．2 High-frequencies part of liquid sloshing(226 ～300)

表1 流體自由晃動頻率計算結(jié)果比較Tab．1 Comparison of results for natural frequencies of liquid sloshing rad/s

由計算結(jié)果和理論解比較發(fā)現(xiàn)，誤差較小，控制在5%以內(nèi)，且所使用的算法比文獻［3］所得結(jié)果更接近理論值，由此可見，本文聲學模型算法的有效性和可靠性得到驗證．

圖3 流體自由表面波動模態(tài)圖Fig．3 Shapes of liquid system corresponding to the sloshing mode

3．2 晃動頻率與容器尺寸關系

通過聲學模型有限元法計算儲液罐液體晃動頻率，研究儲液罐液體在不同容器尺寸時(以容器截面為正方形為例)其晃動頻率的變化規(guī)律．在討論時，將液體深度H 作為一不變量，即液深H=3 m 恒定不變． 容器尺寸(正方形截面邊長a)在a=1，2，3，4，5，6 m 時，分別獲得前幾階自然晃動頻率，取前4 階自然頻率繪出f-a 變化曲線如圖4 所示．從圖中看出:液體前4 階晃動頻率受儲液罐體尺寸影響顯著;隨著儲液罐尺寸a 的增大，液體的自然晃動頻率均漸漸減小，可以說液體晃動頻率與儲液罐邊長基本上保持反比例關系;在儲液罐邊長由1 m 變化至2 m 時，前4 階晃動頻率的變化率均較大，隨后隨著儲液罐邊長的增大晃動頻率變化率較平緩．

圖4 液體晃動頻率與容器邊長的關系Fig．4 Liquid sloshing frequencies versus the length of container

3．3 晃動頻率與液體深度關系

用同樣的方法研究儲液罐液體晃動頻率與儲夜深度之間的變化規(guī)律．在討論時，約定儲液罐截面邊長a 為一恒定量，即a =2 m． 儲液深度H 為1，1．5，2，2．5，3，3．5，4，5，6 m 時，分別獲得前幾階自然晃動頻率，同樣取前4 階自然頻率為例繪出f-H 變化曲線如圖5．結(jié)果表明，第1 階和第2階頻率在充液深度由H =1 m 至H =1．5 m 時變化較大;在液深H=2 m 時，晃動頻率分別達到該階段最大值3．838 2 rad/s 和4．617 9 rad/s，隨后隨著充液深度增大，頻率基本無顯著變化．對于第3 階和第4 階頻率，在充夜深度變化期間其自然晃動頻率基本無明顯變化．總體來說，充液深度對液體晃動頻率影響沒有儲液罐體尺寸對其影響顯著．

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

4 結(jié)論

(1)由計算值與理論解及文獻［3］結(jié)果對比可知，以ABAQUS 軟件平臺編制的分析程序?qū)匦蝺σ汗奕S晃動動力特性分析研究是可行的．分析過程中流體采用ABAQUS 特有的聲單元，流體自由液面單元“削皮”成膜單元，通過彈簧單元連接兩者，實現(xiàn)了位移與應力的兼容，使得計算結(jié)果精度更高，更接近理論值．

(2)流體自由晃動頻率包括低頻和高頻兩部分，低頻部分對應于流體自由晃動時分布在自由液面的動壓力引起的自由液面波動，高頻部分對應于流體內(nèi)部動壓力波動;若要得到更多低階頻率需對橫截面網(wǎng)格細分;流體晃動頻率與儲液罐尺寸基本上保持反比例關系;充液深度對流體體晃動頻率影響沒有儲液罐體尺寸對其影響顯著．

(3)利用該方法可計算任意復雜形狀的三維晃動問題，但對流固耦合問題特別是流體-儲液罐-地基三者耦合問題有待做進一步研究．

［1］ HOUSNER G W．The dynamic behavior of water tanks［J］． Bulletin of Seismological Society of America ，1963，53(2):381 -389．

［2］ WANG Wei，LI Jun-feng，WANG Tian-shu． Modal analysis of liquid sloshing with different contact line boundary conditions using FEM［J］． Science Direct，2008(317):739 -759．

［3］ 賈善坡，許成祥，譚繼可．矩形容器內(nèi)液體三維晃動特性研究［J］． 水電能源科學，2012，30 (1):142-144．

［4］ 鄭建華，李金光，唐輝永． 立式圓柱形儲液罐的三維液固耦合模態(tài)分析［J］． 化工設計，2012，22(1):25 -27．

［5］ 趙曉磊，張榮花，張亮． 15 ×104m3浮放儲罐的模態(tài)分析［J］． 大連民族學院學報，2009，11(1):71-73．

［6］ 包光偉，王振強，張?zhí)煜瑁龋鸺七M劑液體晃動關機響應的數(shù)值仿真［J］． 宇航學報，2002，23(2):84 -88．

［7］ 尚春雨，趙金城． 用FLUENT 分析剛性容器內(nèi)液面晃動問題［J］．上海交通大學學報，2008，42(6):953-956．

［8］ 劉富，童明波，陳建平．基于SPH 方法的三維液體晃動數(shù)值模擬［J］． 南京航空航天大學學報，2010，42(1):122 -126．

［9］ 李青，馬興瑞，王天舒．非軸對稱貯箱液體晃動的等效力學模型［J］．宇航學報，2011，32(2):242 -249．

［10］ HKS． ABAQUS/Standard version 6．5 user’s manual［M］． Rhode Island，USA:Hibbit，Karlsson and Sorensen Inc． 2002．

［11］廖振鵬，工程波動理論導論［M］． 北京:科學出版社，2002．

［12］GRIMES R G，LEWIS J G，SIMON H D． A Shifted block lanczos algorithm for solving sparse symmetric generalized eigenproblems［J］．SIAM Journal on Matrix Analysis and Applications，1994(15):228 -272．