三維大幅晃蕩的Youngs-VOF法模擬

劉永濤，馬寧，2，顧解忡，2

(1.上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240;2.上海交通大學(xué) 海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

伴隨著人類對清潔能源需求的增加以及造船科技的進(jìn)步，超大型液化天然氣船(LNG)自問世以來日益贏得人們的青睞.這類船舶具有大艙容的液體艙室，在復(fù)雜的船舶搖蕩運(yùn)動(dòng)激勵(lì)下，會(huì)產(chǎn)生劇烈的大幅晃蕩運(yùn)動(dòng)，導(dǎo)致液艙結(jié)構(gòu)損毀或運(yùn)動(dòng)性能惡化.因此晃蕩現(xiàn)象成為該類船舶的研究熱點(diǎn)問題.

Kim等［1-2］基于SOLA-VOF程序研究了二維和三維矩形液艙的晃蕩現(xiàn)象，其中自由表面為單值函數(shù)，忽略液面翻卷等劇烈局部流動(dòng)現(xiàn)象.Liu等［3］應(yīng)用VOF方法兩相流模型模擬了帶破碎自由液面的非線性三維晃蕩問題;Delorme［4］采用SPH方法和試驗(yàn)方法研究了二維強(qiáng)制橫搖激勵(lì)下矩形箱內(nèi)淺水晃蕩的砰擊壓強(qiáng)，發(fā)現(xiàn)極大值出現(xiàn)的頻率低于一階固有頻率;Marchandise等［5］應(yīng)用 level set方法，采用快速等值線法精確構(gòu)造界面，研究了低裝載矩形箱中不可壓縮兩相流動(dòng)問題;Kishev等［6］采用CIP方法對橫搖激勵(lì)下不同裝載深度二維矩形液箱的晃蕩及其壓強(qiáng)進(jìn)行了分析;Eswaran等［7］基于VOF法研究了帶制蕩結(jié)構(gòu)的液箱中的晃蕩問題;Godderidge等［8］對近共振頻率橫蕩激勵(lì)作用下二維矩形液箱的晃蕩壓強(qiáng)計(jì)算問題，考慮及箱中兩相的不同的壓縮特性，通過數(shù)值以及試驗(yàn)結(jié)果表明在產(chǎn)生行進(jìn)波的低裝載條件下，由于產(chǎn)生氣體卷入，壓縮模型的選取對抨擊壓強(qiáng)幅值有顯著影響.在國內(nèi)，Zhu等［9-10］采用VOF、Level-set方法預(yù)報(bào)了大幅晃蕩液箱的局部抨擊壓強(qiáng)，Pan等［11］基于MPS方法模擬了晃蕩液箱的自由液面以及局部壓強(qiáng)變化.

本文針對液箱在共振激勵(lì)頻率下的三維大幅晃蕩自由液面大變形現(xiàn)象，以及晃蕩抨擊壓強(qiáng)進(jìn)行數(shù)值模擬分析.鑒于大幅晃蕩過程中產(chǎn)生氣體卷入現(xiàn)象，而表面張力對卷入氣體的形態(tài)與運(yùn)動(dòng)的存在影響，進(jìn)而影響抨擊壓強(qiáng).因此本文也將就該因素的影響進(jìn)行探索分析.

1 液體晃蕩的數(shù)學(xué)模型

為便于描述船舶及其液艙內(nèi)液體的運(yùn)動(dòng)情況，引入2個(gè)坐標(biāo)系:1)固定坐標(biāo)系O0X0Y0Z0，即大地坐標(biāo)系;2)聯(lián)動(dòng)坐標(biāo)系OXYZ，它固結(jié)于船上，隨船舶運(yùn)動(dòng).如圖1所示.

圖1 船舶運(yùn)動(dòng)坐標(biāo)系Fig.1 Shipcoordinate systems

假定液體為不可壓縮的粘性流體，液體在聯(lián)動(dòng)坐標(biāo)系OXYZ中表達(dá)的控制方程為:連續(xù)方程:

動(dòng)量輸運(yùn)方程:

式中:μ為流體動(dòng)力粘性系數(shù)，F(xiàn)b為重力，F(xiàn)st為表面張力，而fb為流體控制域相對于固定坐標(biāo)系運(yùn)動(dòng)而引起的慣性力，可表示為

其中，α為聯(lián)動(dòng)坐標(biāo)系原點(diǎn)O相對慣性坐標(biāo)系的平移加速度，ω為箱體旋轉(zhuǎn)角速度，rp為聯(lián)動(dòng)坐標(biāo)系中流體質(zhì)點(diǎn)位置矢量.

表面張力模型選取為Lafaurie等［12］提出的CSS模型(continuous surface stress).該模型把表面張力轉(zhuǎn)化為質(zhì)量力，具體為

其中，C為流體體積份數(shù)，σ為表面張力.

2 數(shù)值計(jì)算模型

計(jì)算區(qū)域采用交錯(cuò)網(wǎng)格法來具體進(jìn)行差分離散，而連續(xù)方程(1)以及動(dòng)量輸運(yùn)方程(2)采用有限差分法對進(jìn)行離散.

2.1 計(jì)算網(wǎng)格

圖2為流場計(jì)算網(wǎng)格示意圖，在計(jì)算區(qū)域的每個(gè)邊界上，均設(shè)置一排虛擬網(wǎng)格(圖中用虛線表示).具體差分離散時(shí)，采用交錯(cuò)網(wǎng)格法，圖示就是一典型的交錯(cuò)差分網(wǎng)格.這里用(i，j，k)標(biāo)記網(wǎng)格，當(dāng)i=1，i=lmax時(shí)代表非物理網(wǎng)格，即虛擬網(wǎng)格.i=2，i-Im1=Imax-1為x向第一個(gè)和最后一個(gè)物理網(wǎng)格.以 δxi、δyj、δzk表示計(jì)算網(wǎng)格單元的長、寬、高.壓力 Pi，j，k和流體體積函數(shù) Fi，j，k位于單元體中心，速度分量 ui，j，k，vi，j，k與 wi，j，k分別位于單元體右側(cè)面、后面和頂面中心，如圖3所示.

圖2 計(jì)算網(wǎng)格Fig.2 Computation grid

圖3 各變量在計(jì)算網(wǎng)格中的位置Fig.3 Position of variables in cells

2.2 控制方程的數(shù)值解法

連續(xù)方程基于交錯(cuò)網(wǎng)格進(jìn)行中心差分離散.而動(dòng)量方程中對流項(xiàng)采用偏心差分格式(迎風(fēng)格式和中心差分格式的線性組合)離散，擴(kuò)散項(xiàng)則采用中心差分離散.

2.3 速度和壓力的計(jì)算

連續(xù)方程和動(dòng)量方程中，速度和壓力的求解采用超松弛迭代(SOR)法.該方法的基本過程是:先根據(jù)動(dòng)量離散方程得到新時(shí)刻的流場速度估算值，為保證該速度估算值滿足連續(xù)方程，需對速度和壓力同時(shí)進(jìn)行調(diào)整，從而形成逐次迭代過程，直至所有流場單元上滿足精度要求為止.最終迭代得到的u、v、w、p當(dāng)前值即為新時(shí)刻的 u、v、w、p.

2.4 Youngs-VOF重構(gòu)自由表面

該方法基于幾何學(xué)原理，在單個(gè)網(wǎng)格內(nèi)采用斜線段來近似界面.當(dāng)計(jì)算單元為自由面網(wǎng)格時(shí)(即0＜C＜1)，若已知其本身及其周圍網(wǎng)格的C值，則可計(jì)算出該網(wǎng)格內(nèi)界面的法向n().計(jì)算得到內(nèi)界面的法向n以后，則需確定該界面的具體位置.其計(jì)算過程可以邊長為1，且內(nèi)界面法向0＜的單位立方體為例說明.首先預(yù)先確定法向?yàn)閚的內(nèi)界面與立方體各角點(diǎn)相交時(shí)的空間位置及對應(yīng)的流體體積份數(shù)C值，圖4所示為6種典型情況.然后判斷本網(wǎng)格C值所處的范圍，并確定其具體位置.

3 算例

3.1 共振橫蕩激勵(lì)下大幅晃蕩自由表面模擬

液箱尺寸(Hu［6］的試驗(yàn)):長 L=0.6 m，高度H=0.3m，寬度B=0.1m，裝液深度h=0.12m.液箱做η=ηA·sin(wt)的簡諧橫蕩運(yùn)動(dòng)，橫蕩幅值ηA為0.05 m，周期等于0.12 m裝液深度下的流體一階固有周期1.3 s.

圖4 自由面網(wǎng)格中的界面位置Fig.4 Free surface in a cellwith a given C

有效模擬大幅晃蕩現(xiàn)象是一個(gè)重要課題.本算例中，液箱在共振橫蕩激勵(lì)作用下，流體做劇烈晃蕩運(yùn)動(dòng)，產(chǎn)生明顯的翻卷、破碎、合并、沖頂?shù)却笞冃巫杂梢好孢\(yùn)動(dòng).圖5顯示了在不同時(shí)刻，自由液面的數(shù)值模擬結(jié)果(右圖)與試驗(yàn)結(jié)果(左圖，Hu［6］的試驗(yàn))的對比情況.從4個(gè)典型時(shí)刻0.1T～0.4T的對比圖可以看出，自由液面的數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果較為一致.

圖5 自由表面的試驗(yàn)結(jié)果與數(shù)值模擬結(jié)果對比Fig.5 Free surface configurations between experimental results and numerical results

3.2 共振橫蕩激勵(lì)下大幅晃蕩壓強(qiáng)模擬

液箱尺寸:長L=0.6 m，高度 H=0.3 m，寬度B=0.1 m，裝液深度 h=0.06 m.液箱做的簡諧橫蕩運(yùn)動(dòng)，橫蕩幅值ηA為0.06 m，周期等于0.06 m裝液深度下的流體一階固有周期1.7 s.

圖6顯示了在液箱右側(cè)壁面距底0.05 m處，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果(Hu［6］的試驗(yàn))的對比情況.圖6(a)為試驗(yàn)結(jié)果曲線，(b)為壓強(qiáng)時(shí)歷數(shù)值模擬結(jié)果，兩線分別表示不計(jì)表面張力結(jié)果，以及計(jì)入表面張力結(jié)果.圖6中3條曲線對比表明，壓強(qiáng)時(shí)歷數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果近似，而對比兩壓強(qiáng)時(shí)歷數(shù)值模擬結(jié)果可以發(fā)現(xiàn)，計(jì)入表面張力將對壓強(qiáng)時(shí)歷曲線的抨擊壓強(qiáng)峰值產(chǎn)生影響，使得抨擊壓強(qiáng)峰值降低，而除壓強(qiáng)峰值部分外，兩曲線保持一致;總體而言，考慮表面張力的壓強(qiáng)時(shí)歷曲線更相近于試驗(yàn)曲線.究其原因，分析抨擊發(fā)生過程發(fā)現(xiàn)，本算例抨擊壓強(qiáng)峰值發(fā)生于自由液面快速接觸壁面的時(shí)刻，此時(shí)測壓點(diǎn)處于液體表面位置，表面張力的液面穩(wěn)定作用，使得數(shù)值模擬抨擊壓強(qiáng)峰值降低.而在測壓點(diǎn)位于液面以上或以下時(shí)，表面張力基本無影響，壓強(qiáng)時(shí)歷曲線較為一致.

圖6 壓強(qiáng)時(shí)歷試驗(yàn)結(jié)果與數(shù)值計(jì)算結(jié)果對比Fig.6 Comparison of pressure time history between experimental results and numerical results

4 結(jié)論

針對低裝載共振激勵(lì)作用下液箱大幅晃蕩問題，本文應(yīng)用Youngs-VOF法對自由液面大變形現(xiàn)象進(jìn)行了數(shù)值模擬，就液面大變形引起的抨擊壓強(qiáng)進(jìn)行了時(shí)歷預(yù)報(bào)，并與試驗(yàn)結(jié)果進(jìn)行了對驗(yàn)證.此外，就表面張力因素對抨擊壓強(qiáng)數(shù)值預(yù)報(bào)結(jié)果的影響進(jìn)行了對比分析.結(jié)果表明:

1)在共振激勵(lì)頻率條件下，液箱流體運(yùn)動(dòng)出現(xiàn)明顯的大幅晃蕩現(xiàn)象，如翻卷、破碎、沖頂?shù)茸杂梢好娲笞冃维F(xiàn)象，并因此對液箱局部產(chǎn)生抨擊作用;

2)對于三維大幅晃蕩，結(jié)合直角網(wǎng)格以及Youngs-VOF的自由液面重構(gòu)方法能夠有效模擬液面翻卷、破碎等自由液面大變形現(xiàn)象，并能有效預(yù)報(bào)液箱局部抨擊壓強(qiáng);

3)計(jì)入表面張力影響，有助于降低抨擊壓強(qiáng)峰值的數(shù)值模擬結(jié)果，使數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果更一致.

［1］KImY.Numerical simulation of sloshing flows with impact load［J］.Applied Ocean Research，2001，23(1):53-62.

［2］KImY，SHIN Y，LEE K.Numerical study on slosh-induced impactpressures on three dimensional prismatic tanks［J］.Applied Ocean Research，2004，26(5):213-226.

［3］LIU Dongming，LIN Pengzhi.A numerical study of threedimensional liquid sloshing in tanks［J］.Journal of Computational Physics，2008，227(8):3921-3939.

［4］DELORME L，COLAGROSSIA.A set of canonical problems ins loshing，Part I:pressure field in forced roll-comparison between experimental results and SPH［J］.Ocean Engineering，2009，36(2):168-178.

［5］MARCHANDISE E，REMACLE JF.A stabilized finite elementmethod using a discontinuous level set approach for solving two phase incompressible flows［J］.Journal of Computational Physics，2006 ，219(2):780-800.

［6］KISHEV Z R，HU Changhong，KASHIWAGIM.Numerical simulation of violent sloshing by a CIP-based method［J］.Journal of Marine Science Technology，2006，11(2):111-122.

［7］ESWARAN M，SAHA U K，MAITY D.Effect of baffles on a partially filled cubic tank:Numerical simulation and experimental validation［J］.Computers and Structures，2009，87(3):198-205.

［8］GODDERIDGE B，TURNOCK S，EARL C，et al.The effect of fluid compressibility on the simulation of sloshing impacts［J］.Ocean Engineering，2009，36(8):578-587.

［9］ZHU Renqing，WANG Zhidong，F(xiàn)ANG Zhiyong.Prediction of pressure induced by large amplitude sloshing in a tank［J］.Journal of ShipMechanics，2004，8(6):71-78.

［10］ZHU Renqing，F(xiàn)ANG Zhiyong，ZHANG Zhaogang，et al.Level-setmethod for predicting impact pressure induced by violent sloshing in a tank［J］.Journal of ShipMechanics，2008，12(6):344-351.

［11］PAN Xujie，ZHANGHuaixin，LU Yuntao.Numerical simulation of viscous liquid sloshing by moving-particle semiimplicitmethod［J］.Journal of Marine Science and Application，2008，7(3):184-189.

［12］LAFAURIE B，NARDONE C，SCARDOVELLIR，et al.Modelling merging and fragmentation in multiphase flows with SURFER［J］.Journal of Computational Physics，1994，113(1):134-147.