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

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

(1.上海交通大學 船舶海洋與建筑工程學院，上海 200240;2.上海交通大學 海洋工程國家重點實驗室，上海 200240)

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

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

本文針對液箱在共振激勵頻率下的三維大幅晃蕩自由液面大變形現象，以及晃蕩抨擊壓強進行數值模擬分析.鑒于大幅晃蕩過程中產生氣體卷入現象，而表面張力對卷入氣體的形態與運動的存在影響，進而影響抨擊壓強.因此本文也將就該因素的影響進行探索分析.

1 液體晃蕩的數學模型

為便于描述船舶及其液艙內液體的運動情況，引入2個坐標系:1)固定坐標系O0X0Y0Z0，即大地坐標系;2)聯動坐標系OXYZ，它固結于船上，隨船舶運動.如圖1所示.

圖1 船舶運動坐標系Fig.1 Shipcoordinate systems

假定液體為不可壓縮的粘性流體，液體在聯動坐標系OXYZ中表達的控制方程為:連續方程:

動量輸運方程:

式中:μ為流體動力粘性系數，Fb為重力，Fst為表面張力，而fb為流體控制域相對于固定坐標系運動而引起的慣性力，可表示為

其中，α為聯動坐標系原點O相對慣性坐標系的平移加速度，ω為箱體旋轉角速度，rp為聯動坐標系中流體質點位置矢量.

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

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

2 數值計算模型

計算區域采用交錯網格法來具體進行差分離散，而連續方程(1)以及動量輸運方程(2)采用有限差分法對進行離散.

2.1 計算網格

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

圖2 計算網格Fig.2 Computation grid

圖3 各變量在計算網格中的位置Fig.3 Position of variables in cells

2.2 控制方程的數值解法

連續方程基于交錯網格進行中心差分離散.而動量方程中對流項采用偏心差分格式(迎風格式和中心差分格式的線性組合)離散，擴散項則采用中心差分離散.

2.3 速度和壓力的計算

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

2.4 Youngs-VOF重構自由表面

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

3 算例

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

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

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

有效模擬大幅晃蕩現象是一個重要課題.本算例中，液箱在共振橫蕩激勵作用下，流體做劇烈晃蕩運動，產生明顯的翻卷、破碎、合并、沖頂等大變形自由液面運動.圖5顯示了在不同時刻，自由液面的數值模擬結果(右圖)與試驗結果(左圖，Hu［6］的試驗)的對比情況.從4個典型時刻0.1T～0.4T的對比圖可以看出，自由液面的數值模擬結果與試驗結果較為一致.

圖5 自由表面的試驗結果與數值模擬結果對比Fig.5 Free surface configurations between experimental results and numerical results

3.2 共振橫蕩激勵下大幅晃蕩壓強模擬

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

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

圖6 壓強時歷試驗結果與數值計算結果對比Fig.6 Comparison of pressure time history between experimental results and numerical results

4 結論

針對低裝載共振激勵作用下液箱大幅晃蕩問題，本文應用Youngs-VOF法對自由液面大變形現象進行了數值模擬，就液面大變形引起的抨擊壓強進行了時歷預報，并與試驗結果進行了對驗證.此外，就表面張力因素對抨擊壓強數值預報結果的影響進行了對比分析.結果表明:

1)在共振激勵頻率條件下，液箱流體運動出現明顯的大幅晃蕩現象，如翻卷、破碎、沖頂等自由液面大變形現象，并因此對液箱局部產生抨擊作用;

2)對于三維大幅晃蕩，結合直角網格以及Youngs-VOF的自由液面重構方法能夠有效模擬液面翻卷、破碎等自由液面大變形現象，并能有效預報液箱局部抨擊壓強;

3)計入表面張力影響，有助于降低抨擊壓強峰值的數值模擬結果，使數值模擬結果與試驗結果更一致.

［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，FANG Zhiyong.Prediction of pressure induced by large amplitude sloshing in a tank［J］.Journal of ShipMechanics，2004，8(6):71-78.

［10］ZHU Renqing，FANG 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.