賈善坡， 趙友清， 許成祥， 姚華彥

（1.長江大學(xué) 城市建設(shè)學(xué)院，湖北 荊州 434023；2.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，安徽 合肥 230009）

儲罐液體晃動問題廣泛發(fā)生在石油、船舶、化工等領(lǐng)域，常因內(nèi)部儲存流體晃動過大造成容器破壞、重大經(jīng)濟(jì)損失及人員傷亡，因此，如何控制晃動問題引起研究者廣泛關(guān)注［1－2］。目前，常在容器內(nèi)布置隔板作為液體晃動阻尼器來控制液體晃動行為，如條形隔板、環(huán)形隔板及半月形隔板等［3］；文獻(xiàn)［4］運用特征函數(shù)法研究了豎直隔板對矩形儲液罐液體晃動特征頻率的影響；文獻(xiàn)［5］建立了底部安裝豎直隔板的矩形儲液罐數(shù)學(xué)模型，并進(jìn)行了數(shù)值解與實驗結(jié)果的對比研究；文獻(xiàn)［6］采用能量法研究了帶豎直隔板的矩形儲液罐在水平激勵作用下液體晃動阻尼；文獻(xiàn)［7］在線性勢流理論假設(shè)基礎(chǔ)上，使用有限元模擬分析了帶環(huán)形隔板的圓柱容器內(nèi)液體晃動阻尼；文獻(xiàn)［8］使用VOF技術(shù)研究了帶隔板立方儲液罐的防晃特性，并與實驗結(jié)果作了對比分析；文獻(xiàn)［9］采用分離變量法研究了水平簡諧激勵作用下帶環(huán)形隔板圓柱形儲液罐中液體晃動響應(yīng)；文獻(xiàn)［10］通過實驗和數(shù)值模擬研究了帶水平隔板及豎直隔板的矩形儲液罐動態(tài)阻尼特性。

相比于常規(guī)儲罐液體晃動問題，帶隔板的儲液罐液體晃動特性受隔板數(shù)量、半徑及布置位置等因素影響較大，因此，確定合適的隔板參數(shù)及選擇有效的分析方法是儲罐設(shè)計的關(guān)鍵。

本文以有限元軟件ABAQUS為平臺，建立帶隔板圓柱儲液罐三維晃動數(shù)值模型，獲得了帶環(huán)板和無環(huán)板時圓柱儲液罐液體晃動特征頻率、模態(tài)，分析了帶環(huán)板圓柱儲液罐液體晃動動力特性及液體晃動頻率與環(huán)板相對半徑、相對位置等參數(shù)的關(guān)系，并得出液體晃動的一般規(guī)律。

1 液體晃動控制方程

1.1 聲學(xué)方程

聲學(xué)介質(zhì)適合描述無黏、可壓縮的理想流體，其中無剪切應(yīng)力行為，考慮材料阻尼的微幅運動平衡方程［11］為：

其中，p為流體動壓力；x為流體質(zhì)點空間坐標(biāo)；為流體質(zhì)點的速度；為流體質(zhì)點加速度；ρf為流體的密度；γ為體積曳力?？紤]體積模量的聲學(xué)介質(zhì)動態(tài)壓力計算時，本構(gòu)方程為：

其中，Kf為流體體積模量；uf為流體質(zhì)點的位移。

假設(shè)液體為無黏、可壓縮、線性的理想流體，（2）式可表示為：

其中，εV為流體體積應(yīng)變，εV＝ε11＋ε22＋ε33。

1.2 流體動力學(xué)方程

設(shè)流體速度勢函數(shù)為φ（x，y，z，t），則由流體連續(xù)性方程得到不可壓縮性流體域V，滿足：

假設(shè)液體為無黏性理想流體，流體與罐壁具有不可滲透性條件，流體沿罐壁切面自由運動，流體沿罐壁法向速度分量等于0，則在罐壁濕表面?Vw上的邊界條件方程為：

流體自由表面?Vf上滿足運動學(xué)邊界條件，根據(jù)流體動力學(xué)方程及自由表面的伯努利方程，可得：

其中，h為流體自由面波高；g為重力加速度。

令φ（x，y，z，t）＝iωΦ（x，y，z）eiωt，h（x，y，z，t）＝H（x，y，z）eiωt代入（6）式，并聯(lián)立（4）式、（5）式，獲得自由晃動特征問題為：

其中，ω為特征頻率；Φ為對應(yīng)的特征模態(tài)。

引入泛函為：

則（7）式的自由晃動特征問題轉(zhuǎn)化為（8）式的極值解，δΨ＝0。對于三維儲液罐自由晃動問題，將解得的特征向量Φ（x，y，z）代入（8）式，則有：

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

液體自由晃動有限元方程為：

其中，Mf為液體質(zhì)量矩陣；Kf為液體剛度矩陣。

設(shè)（10）式解的表達(dá)式為：

將（12）式代入（10）式，得到液體自由晃動特征方程為：

其中，ω1，ω2，ω3，…，ωn為n個特征頻率；φ1，φ2，φ3，…，φn為相應(yīng)的n個特征振型。

由（13）式提取特征值時，考慮液體晃動系統(tǒng)自由度多，同時需要得到大量特征模態(tài)，為了減少求解時間，使其快速收斂，使用 Lanczos法［12－13］提取特征值，即Lanczos譜變換為：

其中，σ為偏移值，由Gershgorin圓幾何平均值確定；θ為特征值；Φ為特征向量。特征值ω與θ關(guān)系為：

Lanczos程序每次運行創(chuàng)建一序列Krylov子空間，由一系列運算步通過執(zhí)行子空間維數(shù)增長，獲得最逼近的子空間特征向量。Lanczos這種檢測丟失和強(qiáng)制運行的基本模式，理論上只可以計算簡單特征值問題，而對液體晃動這類含大型矩陣特征值問題是費時的。分塊Lanczos法［14］能較好地解決這一問題。

以本文所述方法創(chuàng)建三維圓柱型儲液罐液體晃動模型，流體域采用ABAQUS特有的聲單元AC3D8，自由液面采用膜單元M3D4。重力是液體晃動必不可少的恢復(fù)力，而聲學(xué)單元不可直接附加均布重力載荷，因此，在聲單元和膜單元之間引入彈簧單元，同時也實現(xiàn)了位移與應(yīng)力的兼容。最終液體晃動時，劃分的自由液面微面受到一重力恢復(fù)力ρgdA（dA為微面面積），分析時取ρgdA數(shù)值大小作為引入的單個彈簧剛度。此外，考慮到網(wǎng)格密度對數(shù)值模擬精確度影響，在劃分有限元模型網(wǎng)格時，一個波長內(nèi)布置節(jié)點不少于6～8個，使空間步距Δx＜λ／6，可保證波動數(shù)值模擬準(zhǔn)確［15］。

3 算例分析

為了驗證本文方法的準(zhǔn)確性，首先建立一個無環(huán)板的圓柱儲液罐實例，求解出數(shù)值解，再與同假設(shè)條件下求解出的理論解對比驗證。驗證方法可行后，研究圓柱儲液罐在不同環(huán)板數(shù)量、半徑寬度及位置時的液體晃動特性，如圖1所示。

圖1 帶環(huán)形隔板的圓柱儲液罐

3.1 驗證實例

某圓柱儲液罐，半徑R為40m，液體深度H為30m，密度ρf為1 000kg／m3，聲速c0為1 435m／s，液體體積模量為2.06×109Pa。假設(shè)罐壁剛性，且固定在剛性地基上。

對于液體自由晃動的三維問題，假定液體是理想液體，無黏性、不可壓縮、無旋，則液體晃動的固有頻率為：

其中，ωi為液體晃動第i階頻率；λi為一階Bessel函數(shù)導(dǎo)數(shù)的第i個根。

使用本文所述方法獲得液體晃動前8階模態(tài)，如圖2所示。由（16）式求解圓柱儲液罐液體自由晃動前幾階頻率理論解，結(jié)果見表1所列。

圖2 無環(huán)板儲液罐液體自由表面晃動模態(tài)圖

表1中，m和n分別是橫截面內(nèi)2個相互垂 直方向上的半波數(shù)，本文結(jié)果和理論值誤差在2.5%以內(nèi)，說明了本文方法的有效性和可靠性。

表1 流體自由晃動頻率計算結(jié)果比較 rad／s

3.2 防晃板罐液體自由晃動分析

假設(shè)圓柱儲液罐半徑R為40m，液體深度H為30m；罐內(nèi)置剛性環(huán)板，厚度為0.008m，相對半 徑Ri／R為 0.7，即 環(huán) 板 寬 度 為 12m；h／H＝0.5，即環(huán)板離自由液面高度h為15m；有限元網(wǎng)格密度、材料參數(shù)及邊界條件均一致。計算得到帶單層環(huán)板儲液罐液體晃動前8階模態(tài)及前30階特征頻率，如圖3、圖4所示。

由圖4結(jié)果可以看出，流體自由晃動頻率包括低頻和高頻2部分。低頻部分對應(yīng)于流體自由晃動時分布在自由液面的動壓力引起的自由液面波動；高頻部分對應(yīng)于流體內(nèi)部動壓力波動。加環(huán)板的圓柱儲液罐與不加環(huán)板時相比，由于環(huán)板作用，每階晃動頻率不同程度減小，環(huán)板在儲液罐中起到很好的防晃作用。

圖3 單層環(huán)板儲液罐液體自由表面晃動模態(tài)圖

圖4 流體自由晃動頻率

3.3 晃動頻率與防晃板半徑的關(guān)系

當(dāng)罐內(nèi)分別布置不同半徑環(huán)板且置于不同位置時，研究液體晃動基頻的變化規(guī)律。假設(shè)圓柱儲液罐半徑R為1m，液體深度H為1m，罐內(nèi)置環(huán)板為剛性，厚度為0.004m，有限元網(wǎng)格密度、材料參數(shù)及邊界條件均一致。環(huán)板布置于7個不同相對位置，即h／H分別為 0.01、0.05、0.20、0.50、0.70、0.80、0.90；且在每個對應(yīng)位置時，考慮環(huán)板的6個不同相對半徑，即Ri／R分別為0.1、0.3、0.5、0.7、0.8、0.9，分別獲得前幾階液體晃動頻率。

取第1階晃動頻率分析，繪出基頻與環(huán)板相對半徑曲線圖及基頻與環(huán)板相對位置曲線圖，如圖5、圖6所示。其中，無環(huán)板時圓柱罐液體晃動基頻數(shù)值解ω為3.960 5rad／s，解析解ω為3.990 8rad／s，誤差為0.76%。

圖5 液體晃動基頻與環(huán)板相對半徑關(guān)系

圖6 液體晃動基頻與環(huán)板相對位置關(guān)系

3.4 晃動頻率與防晃板位置的關(guān)系

由圖6可以看出，環(huán)板所處位置對液體晃動基頻有較大影響。隨著環(huán)板位置趨向自由液面時，液體晃動基頻均不同程度減小，特別在h／H＜0.20范圍內(nèi)液體晃動基頻變化率最大；h／H＝0.70時，液體晃動基頻達(dá)到階段最大值，隨后隨著環(huán)板位置繼續(xù)趨向底面時，液體晃動基頻基本無明顯變化，尤其當(dāng)h／H＞0.70、Ri／R＞0.7時，液體晃動基頻幾乎相當(dāng)于不設(shè)置環(huán)板的圓柱罐液體晃動基頻，即ω為3.960 5rad／s。

4 結(jié) 論

（1）由計算值與理論值對比可知，采用本文方法對帶環(huán)板的圓柱儲液罐三維晃動動力特性分析研究是可行的。

（2）環(huán)板對圓柱罐液體具有很好防晃特性，隨著環(huán)板相對半徑的減小，液體晃動基頻均減小，且在Ri／R＜0.7時，對不同位置環(huán)板的液體基頻均影響顯著，在Ri／R＞0.9時，對不同位置環(huán)板的液體基頻均影響不明顯。

（3）隨著環(huán)板位置趨向自由液面時，液體晃動基頻均不同程度減小，h／H＜0.50時，環(huán)板的防晃特性最優(yōu)，在h／H＝0.70時，液體晃動基頻達(dá)到階段最大值，隨著環(huán)板位置繼續(xù)趨向底面時，液體晃動基頻基本無明顯變化。

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

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

［3］ 丁遂亮.充液儲箱液體三維晃動及其動力分析［D］.上海：上海交通大學(xué)，2003.

［4］ Evans D V，Mclver P.Resonant frequencies in a container with a vertical baffle［J］.Fluid Mech，1987，175：295－307.

［5］ Armenio A，Rocca M L.On the analysis of sloshing of water in rectangular containers：numerical study and experimental validation［J］.Ocean Eng，1996，23：705－739.

［6］ Isaacson M，Premasiri S.Hydrodynamic damping due to baffles in a rectangular tank［J］.Civil Eng，2001，28：608－616.

［7］ 楊 蔓，李俊峰，王天舒，等.帶環(huán)形隔板的圓柱儲箱內(nèi)液體晃動阻尼分析［J］.力學(xué)學(xué)報，2006，38（5）：660－667.

［8］ Eswaran M，Saha U K，Maity D.Effect of baffles on a partially filled cubic tank：numerical simulation and experimental validation［J］.Comput Struct，2009，87：198－205.

［9］ 王佳棟，周 叮，劉偉慶.水平激勵下帶環(huán)形剛性隔板圓柱形儲液罐中流體的晃動響應(yīng)［J］.力學(xué)季刊，2011，32（2）：166－172.

［10］ Goudarzi M A，Sabbagh－Yazdi S R.Analytical and experimental evaluation on the effectiveness of upper mounted baffles with respect to commonly used baffles［J］.Ocean Eng，2012，42：205－217.

［11］ ABAQUS／Standard version 6.5User's Manual［M］.Rhode Island：Hibbitt，Karlsson and Sorensen Inc，2002.

［12］ Parlett B N，Nour－Omid B.Toward a black box Lanczos program［J］.Computer Physics Communications，1989，53：169－179.

［13］ Ericsson T，Ruhe A.The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems［J］.Mathematics of Computation，1980，35：1251－1268.

［14］ 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.

［15］ 廖振鵬.工程波動理論導(dǎo)論［M］.北京：科學(xué)出版社，2002：10－63.