帶填充結構圓柱形容器內液體橫向晃蕩特征研究

王肖飛，劉 俊，葉文斌

(1.大連理工大學 海岸與近海工程國家重點試驗室， 遼寧 大連 116024；2.大連理工大學 工程抗震研究所 建設工程學部水利工程學院， 遼寧 大連 116024)

液體晃蕩，廣泛存在于運輸、航天、土木、核電領域等諸多實際工程中，也是儲液罐設計主要考慮因素之一。對液體晃蕩的研究由來已久，基于小波振幅的假設，1954年Senda等[1]率先求得了圓柱形容器液體晃蕩的解析解。近幾年，圓柱形容器液體晃蕩在實驗及數(shù)值研究方面取得了很多的成果，并且FLURNT等流體力學仿真軟件也得到了廣泛的應用[2-3]。忽略流體黏性并假設流體不可壓縮，Tetsuya等[4]研究了地震激勵下帶浮頂圓柱形儲桶的晃蕩響應。Sygulski[5]在不可壓縮勢流假設的基礎上，采用邊界元方法研究了三維液體晃蕩問題。采用變分原理，Takahara等[6]分析了在俯仰激勵下圓柱形容器非線性晃蕩特性，并通過實驗驗證了非線性分析結果的有效性。Lee等[7]通過假設流動為小振幅理想晃動，分析了底部有間隙的環(huán)形圓槽內液體晃蕩的固有頻率。Zhen等[8]應用變量分離和拉普拉斯變換方法研究了圓柱形儲罐內分層液體的三維晃蕩運動。Chen等[9]采用正則化邊界積分法研究了二維和三維容器(包括圓柱形容器)中的非線性液體晃蕩問題。

比例邊界有限元方法(Scaled Boundary Finite Element Method，SBFEM)，由Wolf等[10]在20世紀90年代提出，該方法結合了有限元方法和邊界元方法的優(yōu)點，僅需離散計算域邊界，使計算維度減一。相對于邊界元方法，該方法無需基本解；相對于有限元方法，它具有在比例坐標下徑向解析的特性，從而提高了數(shù)值計算精度。許多學者已經成功將該方法應用于液體晃蕩領域。例如，Lin等[11]將比例邊界有限元方法應用在水平外部激勵下圓柱形或者軸對稱水箱液體晃蕩的分析。Liu等[12]基于線性勢流理論，研究了帶T形擋板的臥式圓柱罐內液體晃蕩問題。應用比例邊界有限元方法，Wang[13]研究了二維橢圓形容器內T形擋板對液體晃蕩響應的影響。本文基于線性勢流理論，應用變分原理推導出帶填充結構圓柱形容器內SBFEM控制方程并引入相應邊界條件進行求解。同時開展了在橫向外部激勵作用下相對波數(shù)、填充介質系數(shù)和填充半徑等參數(shù)對液體晃蕩特征影響的研究。

1 比例邊界有限元方程推導

如圖1所示，假設液體無黏且不可壓縮，圓柱形容器半徑為b，同軸圓柱形填充結構半徑為a，容器內液體深度為H。假設容器底部固定，容器壁為剛性。在x軸方向上作用外部水平位移激勵x=Ae-iωt，其中A為位移的幅值，ω為角頻率，t為時間，同時假設A足夠小，液體自由表面滿足線性勢流理論。

圖1帶填充結構圓柱形容器示意圖

將整個計算域劃分為填充子域Ω1和非填充子域Ω2。其中，填充結構內可采用“固體泡沫球”[14]等填充物，使得該填充結構滿足Sollitt等[15]在1972年所建立的流體與大孔隙多孔介質作用數(shù)學模型。

基于線性勢流理論，分離出時間因子，將速度勢表示為

Φ(x,y,z,t)=φ(x,y,z)e-iωt

(1)

復速度勢函數(shù)φ(x,y,z)在整個計算域內滿足拉普拉斯方程

(2)

及線性自由液面條件

(3)

容器邊界條件

(4)

其中g為重力加速度。

引入邊界條件，將變量z分離出來，速度勢函數(shù)φ(x,y,z)可表示為

(5)

上式右邊第一項代表傳播模態(tài)，第二項代表非傳播模態(tài)。k0和km(m=1,2,…,∞)分別為傳播模態(tài)和非傳播模態(tài)的波數(shù)，且滿足下述色散方程

ω2=gk0tanhk0H

(6)

ω2=-gkmtankmH，m=1,2,…,∞

(7)

計算域中φ0(x,y)和φm(x,y)分別滿足二維Helmholtz方程及修正的Helmholtz方程

(8)

(9)

用φⅠ和φⅡ分別代表計算子域Ω1、Ω2的速度總勢，在r=a處，兩個子域之間的耦合邊界條件可表示為

εφⅠ,n=-φⅡ,n

(10)

(s+if)φⅠ=φⅡ

(11)

其中：ε為多孔介質孔隙度系數(shù)；f為線性阻力系數(shù)；s為慣性系數(shù)。

容器邊壁處邊界條件表示為

φⅡ,n=-iωAcosθ，r=b

(12)

其中：θ為域內任一點與 軸正向所成的角度。

采用變分原理推導比例邊界有限元方程。φ0(x,y)和φm(x,y)在邊界上的法向導數(shù)分別滿足

(13)

(14)

方程(8)、式(9)的等效泛函形式可以表示為

(15)

(16)

圖2比例邊界有限元計算方法示意圖

(17)

(18)

在比例坐標系統(tǒng)中拉普拉斯算子表示為：

(19)

根據(jù)等參變換概念，速度勢函數(shù)可采用插值函數(shù)N(s)進行離散：

φ(ξ,s)=N(s)φ(ξ)=φT(ξ)NT(s)

(20)

將式(19)、式(20)代入式(15)，可得傳播模態(tài)SBFEM控制方程及邊界條件：

(21)

(22)

(23)

式中系數(shù)矩陣

(24)

E1=0

(25)

(26)

同理得非傳播模態(tài)SBFEM控制方程及邊界條件

(27)

(28)

(29)

2 比例邊界有限元方程求解

域Ω1內速度勢可以重新表示為：

(30)

式(21)、式(27)為齊次貝塞爾微分方程，其解的形式可表達為：

(31)

(32)

其中：Tj為n階向量；c0j和cmj為系數(shù)；Jrj(ξ)為第一類貝塞爾函數(shù)；Irj(ξ)為第一類修正的貝塞爾函數(shù)?？紤]貝塞爾函數(shù)和第一類修正貝塞爾函數(shù)的性質并將式(31)代入式(21)，將式(32)代入式(27)，可得

(33)

(34)

同時考慮c0jJrj(ξ)、cmjIrj(ξ)的任意性可得特征方程

(35)

類似地，Ω2內速度勢可以表示為

(36)

其中φⅡ0(ξ)和φⅡm(ξ)分別滿足比例邊界有限元控制方程式(21)和式(27)，以及邊界條件式(22)、式(23)、式(28)和式(29)。方程式(21)和式(27)的解的形式可表達為：

=T(J(ξ)Α0+Y(ξ)B0)

(37)

=T(I(ξ)Αm+K(ξ)Bm)

(38)

其中：a0j、b0j、amj和bmj為系數(shù)；Yrj(ξ)為第二類貝塞爾函數(shù)；Krj(ξ)為第二類修正的貝塞爾函數(shù)。

(39)

(40)

當r=a，由式(9)得：

(41)

將式(31)、式(32)代入式(23)，將式(37)、式(38)代入式(28)，聯(lián)立后式(11)、式(41)表示為：

(s+if)J0aC0=J0aA0+Y0aB0

(42)

(s+if)ImaCm=ImaAm+KmaBm

(43)

(44)

(45)

聯(lián)立式(31)、式(32)、式(37)、式(38)、式(39)、式(40)、式(42)、式(43)、式(44)、式(45)，可求得所有系數(shù)矩陣，代入式(30)、式(36)可求得Ω1、Ω2的速度勢。

水平向晃蕩力Fx以通過對z進行積分得到

(46)

3 計算結果及分析

為了驗證SBFEM方法在填充結構圓柱形容器內液體晃蕩問題的精確性和高效性，將本方法計算結果與解析解[16]進行對比。選取計算參數(shù)如下：填充結構半徑a=0(相當于無填充結構)，容器半徑b=1.0 m，重力加速度g=9.81 m/s2，液體密度ρ=1.0 g/cm3，容器內液面深度H=2.0 m，A=0.1 m，波數(shù)k=0.6。采用4個、12個、20個二階3節(jié)點單元對容器邊壁有限元離散。圖3給出了SBFEM計算結果與解析解歸一化后容器邊壁處自由液面液面高程(η/A)的對比。從圖3中可以看出，在較少的離散單元下，即可得到與解析解相吻合的計算結果，從而驗證了本方法具有很高的精度及計算效率。

圖3容器邊壁液面高程與解析解的對比

用本文提出的方法，對四組不同填充結構系數(shù)情況下容器內晃蕩問題進行了研究。系數(shù)ε、f取值分別為：(1)ε=0.20、f=5.0；(2)ε=0.45、f=2.0；(3)ε=0.80、f=2.0；(4)ε=1.00、f=0.0。其中慣性系數(shù)取定值s=1.0。對于流體而言，即流域內沒有填充物時，孔隙率系數(shù)ε=1.0，阻力系數(shù)f=0，慣性力系數(shù)s=1.0。在此基礎上研究了相對波數(shù)kb、填充結構半徑a、孔隙度系數(shù)ε和線性阻力系數(shù)f等參數(shù)對容器內晃蕩力及自由液面高程的影響，以下為數(shù)值計算的結果。

圖4描述了填充半徑a=0.3、a=0.5時相對波數(shù)kb變化對系統(tǒng)晃蕩力Fx的影響。從圖中可以看出，隨著相對波數(shù)的增加，F(xiàn)x呈現(xiàn)出先增大后減小之后再緩慢增大的趨勢；四組參數(shù)中，隨著ε的增加及f的減小，F(xiàn)x峰值先減小后增大，峰值對應的相對波數(shù)隨之增大；且當ε=0.45、f=2.0時所對應的Fx的峰值最小。

圖4系統(tǒng)晃蕩力隨相對波數(shù)的變化

圖5描述了當填充系數(shù)為ε=0.45、f=2.0時，考慮四種填充半徑a(取值分別為0.1、0.3、0.5、0.7)條件下系統(tǒng)晃蕩力Fx及x正方向上容器邊壁處液面高程η/A隨相對波數(shù)的變化。由圖中可以看出，填充半徑a增大時，F(xiàn)x峰值及η/A峰值均減小，F(xiàn)x及η/A峰值對應相對波數(shù)也隨之減?。划斕畛浒霃絘由0.1增加0.5時，F(xiàn)x峰值及η/A峰值減小幅度較大，當填充半徑a由0.5增加到0.7時，F(xiàn)x峰值及η/A峰值減小趨勢變緩。以上情況表明，填充結構可以有效減緩液體晃蕩的程度，減小晃蕩產生的沖擊壓力對容器的破壞。同時，相對于改變填充結構參數(shù)，改變容器內填充區(qū)域的范圍能取得更好的減晃效果，在實際工程中也易于實現(xiàn)。

圖6描述了填充系數(shù)為ε=0.8、f=0.2和兩種波數(shù)(kb=1.0，2.0)條件下，填充半徑a變化對晃蕩力的影響。圖中可以看出，隨著填充半徑的增加，系統(tǒng)晃蕩力與容器壁處晃蕩力變化趨勢一致，數(shù)值差別很??；填充半徑a增加時，填充結構邊壁處晃蕩力變化較為平穩(wěn)，且始終遠小于系統(tǒng)晃蕩力與容器壁處晃蕩力。同時，由圖6(a)可知，在k=1.0時，隨著晃蕩半徑的增加，系統(tǒng)晃蕩力、容器壁處晃蕩力及填充結構邊壁處晃蕩力增加較為平穩(wěn)；由圖6(b)可知，在k=2.0時，隨著填充半徑a的增加，系統(tǒng)晃蕩

圖5 不同填充半徑下，相對波數(shù)對晃蕩的影響

圖6晃蕩力隨填充半徑a的變化

力與容器邊壁晃蕩力先減小后增大，并且在a=0.6左右時晃蕩力值取得最小。

4 結 論

(1) 與解析解結果對比表明，比例邊界有限元方法在很少的離散單元即可得到很高的計算精度，驗證了本文采用方法在液體晃蕩問題研究的精確性及高效性。

(2) 填充半徑a一定時，孔隙度系數(shù)ε和線性阻力系數(shù)f的變化，對無量綱晃蕩力峰值及峰值對應相對波數(shù)有顯著的影響；所研究四組填充參數(shù)情況中，采用ε=0.45、f=2.0時，晃蕩力峰值取得最小值。

(3) 隨著填充半徑a的增加，無量綱晃蕩力峰值及峰值對應相對波數(shù)都隨之減小；當填充半徑a大于0.5并繼續(xù)增大時，峰值減小趨勢變緩，即此時增大填充區(qū)域后減晃的效果不明顯。