二維矩形彈性液艙內(nèi)液體晃蕩數(shù)值模擬研究

， ，

(大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,遼寧 大連 116024)

以往對(duì)液體晃蕩問(wèn)題的研究均視液艙為剛性液艙[1-3]，這與實(shí)際情況不相符。實(shí)際上液艙都具有一定的彈性。當(dāng)液艙結(jié)構(gòu)為彈性時(shí)，艙內(nèi)液體與結(jié)構(gòu)將會(huì)相互耦合作用，這將會(huì)急劇地改變晃蕩沖擊壓力與液體自由液面的非線性運(yùn)動(dòng)。尤其當(dāng)晃蕩沖擊壓力的作用時(shí)間與彈性結(jié)構(gòu)的濕周期相近或者處于同一量級(jí)時(shí)，較之剛性結(jié)構(gòu)，彈性結(jié)構(gòu)的動(dòng)態(tài)響應(yīng)與液體的壓力可能會(huì)急劇增大，這時(shí)就要研究其中的流固耦合和水彈性問(wèn)題了[4-6]。本文采用ADINA-FSI模塊來(lái)研究縱蕩激勵(lì)下二維矩形彈性液艙內(nèi)的液體晃蕩問(wèn)題。

1 數(shù)值模型

1.1 基本假定

假設(shè)流體為不可壓縮的粘性流體，結(jié)構(gòu)為線彈性材料。流體與結(jié)構(gòu)在流固耦合接觸面上法向不直接脫離，只沿切向滑動(dòng)，即流體與結(jié)構(gòu)在接觸面上任意時(shí)刻法向位移都是相等。

1.2 控制方程

1.2.1 結(jié)構(gòu)運(yùn)動(dòng)基本方程

線彈性結(jié)構(gòu)小變形時(shí)在直角坐標(biāo)系中的運(yùn)動(dòng)方程為

(1)

應(yīng)變率滿足

(2)

Fis——結(jié)構(gòu)體積力分量；

ρs——結(jié)構(gòu)質(zhì)量密度；

1.2.2 流體域內(nèi)N-S方程和連續(xù)方程

流體域內(nèi)ALE描述二維不可壓縮粘性流體的N-S方程和連續(xù)方程為

(3)

式中：u，v——x和y方向的速度；

gx，gy——x和y方向的體積加速度；

p——流體壓力；

ρ——流體密度；

υ——流體運(yùn)動(dòng)粘滯系數(shù)。

1.3 邊界條件

流固耦合模型包括三部分：流體域、結(jié)構(gòu)域和流固耦合接觸面，因此邊界條件由流固耦合接觸面和流體自由表面組成。

1.3.1 流固耦合接觸面上的邊界條件

1)運(yùn)動(dòng)學(xué)條件。在流固耦合接觸面上必須滿足位移協(xié)調(diào)性條件，用拉格朗日坐標(biāo)系可以表示為

(4)

下劃線表示變量?jī)H定義在流固耦合接觸面上。

由式(4)可推導(dǎo)出速度協(xié)調(diào)性條件。

固壁上的不可滑移和可滑移條件為

(5)

n——流固耦合接觸面的法向向量，方向向外。

2)動(dòng)力學(xué)條件。

(6)

1.3.2 流體自由表面條件

自由液面是一個(gè)移動(dòng)的邊界，流體需滿足運(yùn)動(dòng)學(xué)條件和動(dòng)力學(xué)條件。

1)運(yùn)動(dòng)學(xué)條件。它要求組成自由表面的流體質(zhì)點(diǎn)不能流出液面，或者說(shuō)自由液面的法向速度要和流體的法向速度相同，這一條件對(duì)于不同的表面追蹤方法，會(huì)有不同的數(shù)學(xué)表達(dá)式。

2)動(dòng)力學(xué)條件。必須滿足表面應(yīng)力條件，若不計(jì)表面張力，則表面應(yīng)力條件為

(7)

式中：un——垂直于自由表面的法向速度(對(duì)流體域來(lái)講，以外法向?yàn)檎?；

ut——切向速度；

p0——表面大氣壓強(qiáng)。

1.4 流固耦合系統(tǒng)方程

(8)

式中：Ff,Fs——流體與結(jié)構(gòu)相應(yīng)的有限元方程。

ADINA在求解流固耦合問(wèn)題時(shí)有直接耦合求解和迭代耦合求解兩種方法。本文在研究液體晃蕩問(wèn)題時(shí)采用直接耦合法，同時(shí)求解流體運(yùn)動(dòng)和結(jié)構(gòu)的響應(yīng)。直接耦合求解時(shí)采用Newton-Raphson求解耦合方程。

(9)

其中：λd，λτ——流體方程位移松弛因子和結(jié)構(gòu)方程應(yīng)力松弛因子。

求解過(guò)程概括如下：設(shè)初始解X0=Xt，對(duì)k=1,2,3…進(jìn)行迭代，求t+Δt時(shí)刻的解Xt+Δt。

1)整理流體、結(jié)構(gòu)和FSI界面方程組，得到耦合的系統(tǒng)方程組。

2)求解耦合系統(tǒng)的線性化系統(tǒng)方程組并更新得到的解。計(jì)算應(yīng)力和位移殘量并與給定的迭代容差作比較。如果得到的解不收斂，且沒(méi)有達(dá)到FSI迭代的最大次數(shù)，則回到第一步繼續(xù)下一個(gè)迭代；否則，程序停止并顯示不收斂的信息。

3)保存并輸出流體和結(jié)構(gòu)的結(jié)果。

2 ADINA仿真結(jié)果分析

2.1 二維剛性液艙模型

考慮一個(gè)足夠長(zhǎng)的液艙的橫向運(yùn)動(dòng)，將其作為一個(gè)二維問(wèn)題處理，忽略沿長(zhǎng)度方向的流動(dòng)及液艙結(jié)構(gòu)縱向變形的影響。矩形液艙縮比模型見(jiàn)圖1，尺寸為長(zhǎng)L=1.0 m，高H=10 m，靜止水深為h=0.5 m，液艙壁厚δ=10 mm，結(jié)構(gòu)的材料參數(shù)為ρ=7.8×103kg/m3，彈性模量E=210 GPa，泊松比v=0.3；流體采用室溫下的自來(lái)水，密度為ρ=1 000 kg/m3，粘度ν=1.005×10-3Pa·m。在ADINA-Structure中建立二維矩形結(jié)構(gòu)模型，假定為平面應(yīng)變問(wèn)題，劃分為4節(jié)點(diǎn)2-D Solid單元，設(shè)定FSI流固耦合邊界條件；同時(shí)在ADINA-CFD中建立流體模型，劃分為4節(jié)點(diǎn)2-D流體單元，設(shè)定與結(jié)構(gòu)對(duì)應(yīng)的FSI流固耦合邊界條件和自由表面邊界條件，然后在ADINA-FSI中進(jìn)行耦合計(jì)算。

a)模型尺寸 b)ADINA網(wǎng)格

圖1矩形液艙模型

2.2 液艙內(nèi)液體晃蕩固有頻率分析

縱蕩激勵(lì)下液艙內(nèi)晃蕩液體固有頻率為ωn。當(dāng)液艙所受激勵(lì)的頻率等于艙內(nèi)晃蕩液體的固有頻率時(shí)，液體將會(huì)產(chǎn)生大幅晃蕩現(xiàn)象。對(duì)于給定幾何尺寸的液艙而言，晃蕩液體的固有頻率由液深決定，一般由Abramson線性理論公式進(jìn)行估算。

(10)

式中：l——液艙長(zhǎng)度；

n——模態(tài)數(shù)(n=0,1,2,…)。

將本文模型的相關(guān)尺寸代入式(10)，求得液艙內(nèi)液體晃蕩的固有頻率fn(fn=2π/ωn)。然后，采用ADINA來(lái)計(jì)算該模型液體晃蕩的固有頻率fn，將其與線性公式計(jì)算結(jié)果進(jìn)行比較，基本一致，比較結(jié)果見(jiàn)表1。

表1 不同模態(tài)下液體晃蕩固有頻率

2.3 ADINA模型驗(yàn)證

2.3.1 Faltinsen解析解

Faltinsen針對(duì)水平激勵(lì)下二維矩形液艙內(nèi)的液體晃蕩，提出了一個(gè)基于勢(shì)流理論的線性解析解，廣泛地應(yīng)用于各種數(shù)值模型的驗(yàn)證。

二維矩形液艙內(nèi)靜止水深為h，液艙長(zhǎng)度為2a，遭受水平激勵(lì)，所遭受的水平周期激勵(lì)為

ue=-Acosωt

(11)

式中：ue——液艙激勵(lì)速度函數(shù)；

A——速度幅值,A=bω；

b——激勵(lì)位移幅值；

ω——激勵(lì)的圓頻率。

由速度勢(shì)函數(shù)φ，根據(jù)Faltinsen[6]線性解析解，可得到自由面位移η

(12)

(13)

Kn=

原點(diǎn)位于液艙中心、靜止水平面處。

2.3.2 ADINA數(shù)值解與解析解的比較

采用2.1模型，根據(jù)2.2分析得知：其一階固有頻率為f0=0.846 Hz。考慮兩種工況,分析x=a=0.5 m處自由液面的波高歷時(shí)曲線。

1)給定外界激勵(lì)頻率f=0.5f0，位移激勵(lì)振幅A=0.01 m，求解時(shí)間t=10 s。由圖2a)可知，ADINA計(jì)算所得數(shù)值解與解析解吻合較好。

2)給定共振頻率附近f=0.95f0，為避免大幅共振，將位移激勵(lì)振幅調(diào)小為A=0.000 4 m。由圖2b)可知，盡管此次位移激勵(lì)振幅僅為前一情況的4%，但是由于共振效應(yīng)，當(dāng)t=10 s時(shí)，其晃蕩波高幅值也達(dá)到了前者的幅值量級(jí)。

綜上所述，由1)、2)比較可知，采用ADINA模型計(jì)算所得數(shù)值解與Faltinsen解析解吻合較好，說(shuō)明本文ADINA模型有效。

圖2 x=a處自由液面波高歷時(shí)曲線比較

3 彈性液艙內(nèi)的液體晃蕩

3.1 二維彈性液艙模型

為了研究艙壁彈性對(duì)液艙內(nèi)液體晃蕩的影響，將2.1二維剛性液艙的壁厚減小為3和5 mm(液艙模型其它參數(shù)不變)，并且將結(jié)果與10 mm進(jìn)行比較。

3.2 液艙內(nèi)液體晃蕩固有頻率分析

表2 不同模態(tài)下液體晃蕩固有頻率

3.3 彈性液艙內(nèi)晃蕩液體的自由液面運(yùn)動(dòng)

在位移激勵(lì)振幅A=0.002 m的情況下，給定5個(gè)不同的激勵(lì)頻率,分析彈性液艙內(nèi)不同結(jié)構(gòu)剛度、不同激勵(lì)頻率下晃蕩液體的自由液面運(yùn)動(dòng)情況。頻率分布為：f=0.80f0；f=0.90f0；f=0.95f0；f=1.00f0；f=1.10f0，所對(duì)應(yīng)的自由液面的波高歷時(shí)曲線見(jiàn)圖3。

圖3 x=a處自由液面波高歷時(shí)曲線

由圖3可以看出，結(jié)構(gòu)壁厚為δ=3 mm時(shí),液體自由液面波高曲線圖中毛刺較多。這是由于液艙艙壁結(jié)構(gòu)為彈性且彈性較大，在晃蕩過(guò)程中與液體有強(qiáng)烈的流固耦合作用，由于流體對(duì)艙壁的沖擊壓力，彈性艙壁產(chǎn)生了變形，影響了自由液面的波動(dòng)。隨著艙壁剛度的提高，自由液面波高曲線圖中的毛刺情況有所好轉(zhuǎn)。這是由于結(jié)構(gòu)剛度提高，艙壁變形較小，對(duì)自由液面影響較小。

另外發(fā)現(xiàn)，在共振頻率附近(f=0.95f0、f=1.00f0)，自由液面波高曲線的毛刺較小，結(jié)構(gòu)的彈性對(duì)晃蕩液體的波動(dòng)影響較小；而其它非共振頻率時(shí)自由液面波高曲線的毛刺較多，結(jié)構(gòu)的彈性對(duì)晃蕩液體的波動(dòng)影響較大。

圖3f)為位移激勵(lì)振幅為A=0.002 m下3個(gè)不同剛度、5個(gè)不同頻率下液艙內(nèi)晃蕩液體的自由液面波高曲線的最大值。可以看出，由于流固耦合作用，彈性液艙內(nèi)液體的一階固有頻率要小于剛性液艙內(nèi)液體的一階固有頻率，與表2采用ADINA計(jì)算所得的晃蕩液體的固有頻率的比較結(jié)果一致。因此，晃蕩液體的共振頻率已不能按照剛性液艙內(nèi)的頻率來(lái)計(jì)算，在結(jié)構(gòu)設(shè)計(jì)時(shí)應(yīng)給予足夠重視。

4 結(jié)論

1)通過(guò)與解析解的比較，基于有限元軟件ADINA建立的二維矩形液艙晃蕩模型的數(shù)值模擬是有效的。

2)彈性液艙內(nèi)晃蕩液體與艙壁有強(qiáng)烈的流固耦合作用，由于流體對(duì)艙壁的沖擊壓力，彈性艙壁產(chǎn)生了變形，影響了自由液面的波動(dòng)情況；隨著結(jié)構(gòu)剛度的增大，影響逐漸減小。

4)由于流固耦合作用的影響，彈性液艙內(nèi)液體的一階固有頻率要小于剛性液艙內(nèi)液體的一階固有頻率。

[1] LIU Dongming, LIN Pengzhi. A numerical study of three-dimensional liquid sloshing in tanks[J].Journal of Computational Physics,2008(227):3921-3939.

[2] AKYILDIZA H, UNAL E. Experimental investigation of pressure distribution on a rectangular tank due to the liquid sloshing[J].Ocean Engineering,2005(32):1503-1516.

[3] NASAR T, SANNASIRAJ S A. Experimental study of liquid sloshing dynamics in a barge carrying tank[J].Fluid Dynamics Research,2008(40):427-458.

[4] LEE DONG Y. Study on sloshing in cargo tanks including hydro elastic effects[J]. Journal of Marine Scienceand Technology,1999(4):27-34.

[5] CHO J R, PARK S W, KIM H, et al. Hydroelastic analysis of insulation containment of LNG carrier by global-local approach[J].International Journal For Numerical methods in engineering,2008(76):749-774.

[6] ESWARAN M, SAHA U K, MAITY D. Effect of baffles on a partially filled cubic tank: Numerical simulation and experimental validation[J].Computers & Structures, 2009(87):198-205.

[7] WU G X, MA Q W, TAYLOR R E. Numerical simulation of sloshing waves in a 3D tank based on a finite element method[J].Applied Ocean Research,1998(20): 337-355.