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

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

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

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

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

1 液體晃動(dòng)控制方程

1.1 聲學(xué)方程

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

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

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

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

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

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

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

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

流體自由表面?Vf上滿足運(yùn)動(dòng)學(xué)邊界條件，根據(jù)流體動(dòng)力學(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）式，獲得自由晃動(dòng)特征問題為：

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

引入泛函為：

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

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

液體自由晃動(dòng)有限元方程為：

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

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

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

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

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

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

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

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

3 算例分析

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

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

3.1 驗(yàn)證實(shí)例

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

對(duì)于液體自由晃動(dòng)的三維問題，假定液體是理想液體，無(wú)黏性、不可壓縮、無(wú)旋，則液體晃動(dòng)的固有頻率為：

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

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

圖2 無(wú)環(huán)板儲(chǔ)液罐液體自由表面晃動(dòng)模態(tài)圖

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

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

3.2 防晃板罐液體自由晃動(dòng)分析

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

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

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

圖4 流體自由晃動(dòng)頻率

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

［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］ 丁遂亮.充液儲(chǔ)箱液體三維晃動(dòng)及其動(dòng)力分析［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)形隔板的圓柱儲(chǔ)箱內(nèi)液體晃動(dòng)阻尼分析［J］.力學(xué)學(xué)報(bào)，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］ 王佳棟，周 叮，劉偉慶.水平激勵(lì)下帶環(huán)形剛性隔板圓柱形儲(chǔ)液罐中流體的晃動(dòng)響應(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òng)理論導(dǎo)論［M］.北京：科學(xué)出版社，2002：10－63.