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

，

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

液體晃蕩是極為復(fù)雜的一種流體運(yùn)動(dòng)形式，它具有高度的非線性和隨機(jī)性。在以往的研究中，液體晃蕩問(wèn)題大部分的研究工作主要集中在剛性液艙數(shù)值模擬方面[1-3]，然而實(shí)際工程中所使用的液艙都有一定的彈性，假定液艙是剛性結(jié)構(gòu)，忽略了液艙結(jié)構(gòu)的變形，同時(shí)也忽略了液艙變形與內(nèi)部流體的相互作用。近年來(lái)有部分學(xué)者對(duì)彈性液艙的液體晃蕩問(wèn)題進(jìn)行了一些研究[4-6]，對(duì)于彈性液艙晃蕩，必須研究晃蕩過(guò)程中帶自由液面流體的非線性動(dòng)力學(xué)問(wèn)題及其結(jié)構(gòu)變形的相互耦合問(wèn)題。

1 數(shù)值方法

假設(shè)矩形液艙為線彈性結(jié)構(gòu)，液體為不可壓縮的粘性流體。流固耦合接觸面上，流體和結(jié)構(gòu)不直接脫離，即流體與結(jié)構(gòu)在接觸面上任意時(shí)刻法向位移都相等。

采用有限元法對(duì)流體域連續(xù)方程、運(yùn)動(dòng)方程、質(zhì)量守恒方程等進(jìn)行離散求解；采用隱式積分法對(duì)結(jié)構(gòu)域幾何方程、運(yùn)動(dòng)方程、物理方程等進(jìn)行積分求解[6]。

ADINA在求解流固耦合問(wèn)題時(shí)有直接耦合求解和迭代耦合求解。本文采用直接耦合法進(jìn)行計(jì)算，將流體域、結(jié)構(gòu)域和耦合作用構(gòu)造在同一控制方程——流固耦合系統(tǒng)方程中，在單一時(shí)間步內(nèi)，對(duì)流固耦合系統(tǒng)方程所有變量同時(shí)進(jìn)行求解。采用位移-壓力格式描述流固耦合系統(tǒng)。

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

(1)

式中:X——耦合系統(tǒng)的解向量，X=(Xf,Xs)。

其中：Xf,Xs——定義在流體和結(jié)構(gòu)節(jié)點(diǎn)上的解向量。

因此流固耦合面上結(jié)構(gòu)的位移和流體應(yīng)力分別表示為ds=ds(Xs)，τf=τf(Xf)。Ff,Fs分別是流體與結(jié)構(gòu)相應(yīng)的有限元方程。

對(duì)于單獨(dú)的流體或結(jié)構(gòu)模型來(lái)說(shuō)，其線性方程如下。

1)流體方程。

(2)

2)結(jié)構(gòu)方程。

(3)

將流體方程和結(jié)構(gòu)方程放在一個(gè)方程組中

(4)

流體應(yīng)力和結(jié)構(gòu)位移分別使用了應(yīng)力松弛因子λτ(0<λτ<1)和位移松弛因子λd(0<λd<1)。

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

分別整理流體方程和結(jié)構(gòu)方組，得到耦合的方程組。

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

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

2 模型建立與驗(yàn)證

2.1 模型建立

采用有限元軟件ADINA分別建立三維矩形剛性模型與對(duì)應(yīng)的彈性模型。

2.1.1 三維剛性液艙模型

矩形液艙計(jì)算模型示意見(jiàn)圖1。

圖1 ADINA液艙計(jì)算模型示意

模型尺寸為L(zhǎng)=2a=0.5 m，B=2b=0.1 m，H=0.5 m，載液率r為16%和50%。

液艙厚度δ=12 mm，密度ρ=1.185×103kg/m3，泊松比ν=0.35，彈性模量E=2.1×1011Pa。流體采用自來(lái)水，粘度ν=1.005×10-3Pa·s，密度ρ=1×103kg/m3。在ADINA-Structure中建立結(jié)構(gòu)模型，將結(jié)構(gòu)定義為Shell 單元，設(shè)定FSI流固耦合邊界條件；同時(shí)在ADINA-CFD中建立流體模型，劃分為3-D fluid 單元，設(shè)定與結(jié)構(gòu)對(duì)應(yīng)的FSI流固耦合邊界條件和自由表面邊界條件，然后在ADINA-FSI中進(jìn)行耦合計(jì)算，ADINA網(wǎng)格見(jiàn)圖2。

2.1.2 三維彈性液艙模型

為研究液艙彈性對(duì)晃蕩的影響，將模型左右側(cè)艙壁和頂蓋的壁厚改為2 mm，其余保持不變。

圖2 ADINA網(wǎng)格

2.2 模型驗(yàn)證

2.2.1 液體固有頻率的分析

縱蕩激勵(lì)下、液艙內(nèi)晃蕩液體的固有頻率為fn。對(duì)于給定幾何尺寸的液艙而言，晃蕩液體的固有頻率可以由Abramson線性理論公式進(jìn)行估算：

(5)

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

h——液艙內(nèi)液深；

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

運(yùn)用ADINA建立三維矩形液艙模型來(lái)計(jì)算該液艙內(nèi)晃蕩液體的固有頻率并將其與Abramson線性公式計(jì)算結(jié)果進(jìn)行比較，結(jié)果見(jiàn)表1。

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

注：fElastic為2.1.2彈性模型的計(jì)算頻率；fRigid為2.1.1剛性模型的計(jì)算頻率。

由表1可見(jiàn)，剛性液艙內(nèi)液體各階固有頻率與理論計(jì)算所得液體各階固有頻率非常接近；同時(shí)由于流固耦合作用的影響，當(dāng)液艙為低液深時(shí)，彈性液艙內(nèi)液體一階固有頻率大于剛性液艙內(nèi)液體的一階固有頻率；反之當(dāng)液艙為高液深時(shí)，彈性液艙內(nèi)液體一階固有頻率小于剛性液艙內(nèi)液體的一階固有頻率。

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

Faltinsen針對(duì)水平激勵(lì)下的二維矩形液艙內(nèi)的液體晃蕩，提出了一個(gè)基于勢(shì)流理論的線性解析解，廣泛地應(yīng)用于各種數(shù)值模型的驗(yàn)證。本文將該解析解由二維情況擴(kuò)展到三維情況[7]，可得到自由液面位移η。

(6)

(7)

采用2.1.1模型，根據(jù)2.2分析得知：載液率r為50%時(shí),液艙內(nèi)液體一階固有頻率為fRigid=1.197 Hz。給定外界激勵(lì)振幅A=0.004 m，采用兩種工況對(duì)A點(diǎn)處自由液面的波高歷時(shí)曲線進(jìn)行對(duì)比分析，見(jiàn)圖3。

圖3 A點(diǎn)處自由液面波高歷時(shí)曲

由圖3比較可知，采用ADINA模型計(jì)算所得到的數(shù)值解與Faltinsen解析解吻合較好，說(shuō)明本文的ADINA模型有效。

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

為了驗(yàn)證三維彈性模型的正確性及可行性，對(duì)模型進(jìn)行物理模型實(shí)驗(yàn)，并將結(jié)果與ADINA模型模擬結(jié)果進(jìn)行對(duì)比。

將壓力傳感器布設(shè)在液艙的右側(cè)艙壁上，標(biāo)簽“PN”代表壓力傳感器PN處的晃蕩壓力數(shù)據(jù)，N=1，2，3，…見(jiàn)圖4。其中P1，P2，P3號(hào)壓力傳感器分別距離底板80、125、250 mm。本文只針對(duì) P1，P2，P3號(hào)壓力點(diǎn)進(jìn)行了壓力分析。

圖4 壓力傳感器P布置示意

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

對(duì)上述模型分別按表2各工況進(jìn)行模擬，各工況激勵(lì)振幅A=0.004 m。

表2 計(jì)算工況

注：fA，fE為彈性液艙內(nèi)液體的一階固有頻率，且由于r=50%時(shí)，ADIAN計(jì)算模型的一階固有頻率fA=1.194 Hz與對(duì)應(yīng)的實(shí)驗(yàn)?zāi)Ｐ鸵浑A固有頻率fE=1.090 Hz具有一定的差距，為了達(dá)到共振效果，ADINA模型采用的外界激勵(lì)頻率fA=1.194 Hz。

分析A點(diǎn)處自由液面的波高歷時(shí)曲線，并且將計(jì)算結(jié)果與實(shí)驗(yàn)值進(jìn)行對(duì)比，見(jiàn)圖5。

圖5 A點(diǎn)處自由液面的波高歷時(shí)曲線

在數(shù)值模擬過(guò)程中觀察到波面的三維情況，圖6和圖7中分別給出了晃蕩液體模型不同時(shí)刻的三維波面圖及晃蕩液體模型在A和C點(diǎn)的波高歷時(shí)曲線對(duì)比圖。

圖6 Case2液體模型在t=0.0、6.04、12.4及21.2 s時(shí)的三維波面圖

圖7 A、B點(diǎn)處自由液面波高歷時(shí)曲線對(duì)比圖

綜合分析圖5～7，對(duì)比不同載液率液艙在共振頻率下液艙晃蕩自由液面波高歷時(shí)曲線及波面圖可以得出以下結(jié)論。

1)從圖5～7可見(jiàn)，基于有限元數(shù)值模擬軟件ADINA的數(shù)值模擬結(jié)果與實(shí)驗(yàn)值基本吻合，驗(yàn)證了彈性模型的有效性。

2)從波高歷時(shí)曲線圖5中可以看出，液艙內(nèi)液體晃蕩劇烈，液體在壁面形成爬高。圖5a)和圖5b)分別為低、高液深時(shí)A點(diǎn)處波高歷時(shí)曲線，從圖中可以看出在兩種載液水平下，波高歷時(shí)曲線的形式有所不同。低液深液體晃蕩自由表面運(yùn)動(dòng)比高液深時(shí)非線性現(xiàn)象更加明顯，高液深共振時(shí)，波面出現(xiàn)橢圓余弦波。

3)從圖6a)可以看出，當(dāng)t=0 s時(shí)，波面靜止；當(dāng)t=6.04s時(shí)，由圖6b)中波面右側(cè)壁看出，波面中間凹陷，隨著晃蕩的劇烈，如圖6c)和圖6d)所示，波面的三維現(xiàn)象更加明顯。圖7a)、7b)為Case1、Case2在A和B點(diǎn)處波高歷史曲線對(duì)比圖?？梢钥闯?，各曲線的峰值及歷時(shí)曲線形狀具有區(qū)別，得出本文采用的模型在晃蕩時(shí)具有一定三維現(xiàn)象，屬有效的三維模型。

3.2 彈性液艙內(nèi)晃蕩液體的動(dòng)水壓強(qiáng)

為了研究裝載率對(duì)彈性液艙內(nèi)液體晃蕩水動(dòng)力特性，仍然采用2.1.2的液艙模型，對(duì)液艙內(nèi)液體各觀測(cè)點(diǎn)進(jìn)行動(dòng)力分析，與實(shí)驗(yàn)值進(jìn)行對(duì)比，各工況下的動(dòng)水壓強(qiáng)歷時(shí)曲線見(jiàn)圖8～10。

圖8 P1點(diǎn)處動(dòng)水壓強(qiáng)歷時(shí)曲線

圖9 P2點(diǎn)處動(dòng)水壓強(qiáng)歷時(shí)曲線

圖10 P3點(diǎn)處動(dòng)水壓強(qiáng)歷時(shí)曲線

綜合分析圖8～10，對(duì)比不同載液率液艙在共振頻率下各觀測(cè)點(diǎn)動(dòng)水壓強(qiáng)歷時(shí)曲線可以看出：

1)ADINA模型動(dòng)水壓強(qiáng)模擬結(jié)果與實(shí)驗(yàn)值基本相符，可見(jiàn)模型驗(yàn)證具有有效性。

2)圖8a)為載液率為16%時(shí)P1觀測(cè)點(diǎn)的動(dòng)水壓強(qiáng)歷時(shí)曲線，圖8b)，圖9，圖10為載液率為50%時(shí)P1，P2，P3觀測(cè)點(diǎn)的動(dòng)水壓強(qiáng)歷時(shí)曲線，可以看出彈性液艙晃蕩動(dòng)水壓強(qiáng)歷時(shí)曲線出現(xiàn)雙峰現(xiàn)象。其中圖8a)和圖10分別為低液深和高液深時(shí)自由表面附近點(diǎn)所受的動(dòng)水壓強(qiáng)，從圖中可以明顯地看出沖擊現(xiàn)象；圖8b)和圖9為高液深水面以下點(diǎn)P1和P2點(diǎn)的動(dòng)水壓強(qiáng)歷時(shí)曲線，從圖中可以看出壓力雙峰中前后兩峰值相差不大，但后一峰值的持續(xù)時(shí)間大于前一峰值。且對(duì)比圖8b)、圖9和圖10可知，動(dòng)水壓強(qiáng)隨液面到底部逐漸減小，這點(diǎn)與理論相符。

由三維波面圖可知，液艙側(cè)壁中間部分與同一水平線對(duì)應(yīng)各角點(diǎn)的所受水體沖擊有所不同，為了進(jìn)一步分析了解各角點(diǎn)與對(duì)應(yīng)側(cè)壁面中間部分的動(dòng)水壓強(qiáng)的區(qū)別，增加了對(duì)側(cè)壁各角點(diǎn)的壓力觀測(cè)點(diǎn)PC1、PC3，其在ADINA中的壓強(qiáng)傳感器分布見(jiàn)圖1，其中PC1、PC3號(hào)壓力傳感器分別與相應(yīng)的P1和P3號(hào)壓力傳感器處于同一水平位置。對(duì)液艙內(nèi)液體各觀測(cè)點(diǎn)進(jìn)行動(dòng)水壓強(qiáng)分析，各工況下的側(cè)壁中點(diǎn)與角點(diǎn)動(dòng)水壓強(qiáng)歷時(shí)曲線對(duì)比見(jiàn)圖11～12。

圖11 P1,PC1點(diǎn)處動(dòng)水壓強(qiáng)歷時(shí)曲線

圖12 P3,PC3點(diǎn)處動(dòng)水壓強(qiáng)歷時(shí)曲線

從圖11～12可以看出，由于三維現(xiàn)象的存在，各角點(diǎn)的壓強(qiáng)值略大于各壁面中點(diǎn)的壓強(qiáng)。

4 結(jié)論

1)通過(guò)與解析解和實(shí)驗(yàn)值的比較，基于有限元軟件ADINA建立的三維矩形彈性液艙晃蕩數(shù)值模型是正確及有效的。

2)由于流固耦合作用的影響，當(dāng)液艙內(nèi)為低液深時(shí)，彈性液艙內(nèi)液體一階固有頻率大于剛性液艙內(nèi)液體的一階固有頻率；反之當(dāng)液艙內(nèi)為高液深時(shí)，彈性液艙內(nèi)液體一階固有頻率小于剛性液艙內(nèi)液體的一階固有頻率。

3)兩種載液水平下，波高歷時(shí)曲線的形式有所不同。低液深液體晃蕩自由表面運(yùn)動(dòng)比高液深時(shí)非線性現(xiàn)象更加明顯，且自由液面歷時(shí)曲線出現(xiàn)了雙峰。高液深共振時(shí)，波面出現(xiàn)橢圓余弦波。

4)三維彈性液艙液體晃蕩劇烈沖擊艙壁，使得自由表面出現(xiàn)三維效應(yīng)；另外，動(dòng)水壓強(qiáng)歷時(shí)曲線出現(xiàn)雙峰現(xiàn)象，前后兩峰值雖相差不大，但后一峰值的持續(xù)時(shí)間大于前一峰值，且在同一水平面上各角點(diǎn)的沖擊壓力略大于艙壁中點(diǎn)，為后續(xù)研究提供了參考依據(jù)。

本文對(duì)三維矩形彈性液艙內(nèi)液體晃蕩的模擬只做了小幅振蕩研究，而大幅振蕩特別當(dāng)液面出現(xiàn)翻卷、破碎等復(fù)雜情況時(shí)液體對(duì)艙壁抨擊力將會(huì)更大，結(jié)構(gòu)穩(wěn)定性將會(huì)更低，將更加具有危害性。在今后的研究中應(yīng)對(duì)大幅振蕩進(jìn)行探討。

[1] FALTINSEN O M.A numerical nonlinear method of sloshing in tanks with two-dimensional flow[J].Journal of Ship Research,1978(3):193-202.

[2] AKYILDIZ H,UNAL E.Sloshing in a three-dimensional rectangular tank-numerical simulation and experiment validation[J].Ocean Engineering,2006(23):2135-2149.

[3] 朱仁慶.盛液容器內(nèi)液體二維晃蕩的數(shù)值模擬[J].華東船舶工業(yè)學(xué)院學(xué)報(bào),1998,12(2):14-12.

[4] GRADINSCAK M.Liquid sloshing in containers with flexibility[D].Australia:Engineering and Science Victoria University Melbourne,2009.

[5] 朱仁慶.液體晃蕩及其結(jié)構(gòu)的相互作用[D].無(wú)錫：中國(guó)船舶科學(xué)研究中心,2001.

[6] 張秋艷，任 冰，蔣梅榮.二維矩形彈性液艙內(nèi)液體晃蕩的數(shù)值模擬研究[J].船海工程,2012,(4):11-15.

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