，

(大連理工大學(xué) 海岸和近海工程國家重點實驗室，遼寧 大連 116023)

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

1 數(shù)值方法

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

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

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

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

(1)

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

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

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

對于單獨的流體或結(jié)構(gòu)模型來說，其線性方程如下。

1)流體方程。

(2)

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

(3)

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

(4)

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

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

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

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

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

2 模型建立與驗證

2.1 模型建立

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

2.1.1 三維剛性液艙模型

矩形液艙計算模型示意見圖1。

圖1 ADINA液艙計算模型示意

模型尺寸為L=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。流體采用自來水，粘度ν=1.005×10-3Pa·s，密度ρ=1×103kg/m3。在ADINA-Structure中建立結(jié)構(gòu)模型，將結(jié)構(gòu)定義為Shell 單元，設(shè)定FSI流固耦合邊界條件；同時在ADINA-CFD中建立流體模型，劃分為3-D fluid 單元，設(shè)定與結(jié)構(gòu)對應(yīng)的FSI流固耦合邊界條件和自由表面邊界條件，然后在ADINA-FSI中進行耦合計算，ADINA網(wǎng)格見圖2。

2.1.2 三維彈性液艙模型

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

圖2 ADINA網(wǎng)格

2.2 模型驗證

2.2.1 液體固有頻率的分析

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

(5)

式中：L——液艙長度；

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

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

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

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

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

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

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

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

(6)

(7)

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

圖3 A點處自由液面波高歷時曲

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

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

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

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

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

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

對上述模型分別按表2各工況進行模擬，各工況激勵振幅A=0.004 m。

表2 計算工況

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

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

圖5 A點處自由液面的波高歷時曲線

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

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

圖7 A、B點處自由液面波高歷時曲線對比圖

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

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

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

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

3.2 彈性液艙內(nèi)晃蕩液體的動水壓強

為了研究裝載率對彈性液艙內(nèi)液體晃蕩水動力特性，仍然采用2.1.2的液艙模型，對液艙內(nèi)液體各觀測點進行動力分析，與實驗值進行對比，各工況下的動水壓強歷時曲線見圖8～10。

圖8 P1點處動水壓強歷時曲線

圖9 P2點處動水壓強歷時曲線

圖10 P3點處動水壓強歷時曲線

綜合分析圖8～10，對比不同載液率液艙在共振頻率下各觀測點動水壓強歷時曲線可以看出：

1)ADINA模型動水壓強模擬結(jié)果與實驗值基本相符，可見模型驗證具有有效性。

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

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

圖11 P1,PC1點處動水壓強歷時曲線

圖12 P3,PC3點處動水壓強歷時曲線

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

4 結(jié)論

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

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

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

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

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

[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é)報,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].無錫：中國船舶科學(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.