一種幾何VOF方法在液艙晃蕩流動模擬中的應用

李金龍， 尤云祥， 陳 科

(上海交通大學 海洋工程國家重點實驗室； 高新船舶與深海開發裝備協同創新中心， 上海 200240)

浮式天然氣生產儲卸船通過低溫的方式將在海上開采到的天然氣儲存在船體內部液艙.這些液艙沒有裝載率的限制，部分裝載的液艙在海洋波浪的激勵下產生晃蕩，影響船舶的運動、穩性及其結構安全，是工程上關注的熱點問題[1].近年來，隨著計算機技術的發展與進步，計算流體力學(CFD)變得越來越流行.CFD 的兩相流求解器在模擬液艙晃蕩流動時遇到的主要問題之一是如何處理自由液面的不連續性以及如何表達物理量在自由液面處的階躍[2].確定自由液面位置主要有界面追蹤和界面捕捉兩類方法[3].界面追蹤方法追蹤相界面的運動，例如 MAC (Marker-and-Cell)[4]、FT (Front-Tracking)[5]、ALE (Arbitrary Lagrangian-Eulerian)[6]以及σ-transformation 等方法[7].界面捕捉方法基于固定的Euler網格，借助一個標量場分布捕捉相界面的位置，允許較大的相界面變形、破碎以及合并等復雜拓撲變化.典型的界面捕捉方法包括 VOF和LS(Level Set)方法[3].VOF 方法采用體積分數率的分布捕捉兩相界面的位置，且能保證質量守恒，尤其適用于液艙晃蕩類封閉空間內的液體流動模擬計算.該方法又可以分為代數VOF和幾何VOF方法[2].代數VOF方法適用于任意三維網格，利用高精度格式和壓縮差分格式求解體積分數率輸運方程，以保證體積分數率的有界性，例如 OpenFOAM 中的 MULES (Multi-dimensional Universal Limi-ter with Explicit Solution).幾何 VOF 方法需進行界面重構，并計算網格單元和相界面的交線，再用重構出的幾何相界面更新體積分數率的分布.例如:SOLA (Solution Algorithm)-VOF 應用貢獻者-接受者通量近似的方法重構相界面形狀[8]；SLIC (Simple Line Interface Calculation)[9]和 PLIC (Piecewise Linear Interface Calculation)[10]通過求取當地相界面法向量重構幾何相界面.采用這些界面重構方法都能獲得明銳(相界面不連續性明顯)的兩相界面，但是這些方法通常僅適用于結構化網格，當推廣到三維網格時會引起復雜的幾何運算，而推廣到三維非結構化網格時則遇到了困難.最近被提出的isoAdvector方法構造了isoface重構幾何界面，并根據 isoface 的運動更新體積分數率的分布.相對于代數 VOF 方法而言，isoAdvector 方法能捕捉更為明銳和準確的相界面[2]，可適用于三維任意多面體網格.

本文采用動網格技術處理液艙運動，且避免引入液艙平移和旋轉慣性力求解非慣性坐標系下的動量方程，避免處理慣性力涉及的導數[11-12].采用isoAdvector 方法捕捉自由液面并考慮動網格的影響，引入了運動通量修正，提出了面-界面交線運動修正.采用一種近似公式，根據網格單元面上的運動通量重構出網格單元體心的剛體運動速度.基于不同的 VOF 方法對受迫非共振晃蕩、受迫共振晃蕩以及單次沖擊波面進行了模擬，并將模擬結果與試驗結果以及解析解進行了比較.此外，提出了新的估計兩相界面厚度的方法.

1 數學模型

基于開源平臺 OpenFOAM，采用有限體積法離散 RANS (Reynolds-averaged Navier-Stokes) 方程.

1.1 控制方程

兩相流被視為一種單一的連續介質，可采用一套控制方程來求解連續介質的流動.連續方程和動量方程采用張量表達：

(1)

(2)

式(2)的右端項涉及變換

(3)

gi=0

1.2 幾何 VOF 方法——isoAdvector

isoAdvector 方法分為兩個步驟：① 根據當前的體積分數率場，在每一個界面網格單元(該單元的體積分數率介于0和1之間)重構出一系列的幾何平面(isoface)，每一個 isoface將其所在的界面單元分割為兩個子單元，這兩個子單元分別被兩相流體填充，其體積比和該網格單元的體積分數率一致.isoface的集合可以視為兩相界面的一種表達，但該集合所表示的幾何面在各個網格單元交界處可能存在不連續性[2].② 利用速度場預測這些 isofaces 的運動，從而更新體積分數率場.此更新操作不涉及求解線性代數方程組的問題，避免了選取高精度和壓縮差分格式.由于本文利用動網格技術來處理液艙的運動，所以需要對 isoAdvector 方法進行下列兩個方面的修正.

首先通過對體積分數率更新公式的重新推導，引入第一個修正.定義相指示函數[2]

(4)

式中：x是空間位置坐標；ρ1和ρ2分別為液相和氣相的密度.如果x=(xi,xj,xk)處于液相，則H=1；如果x處于氣相，則H=0.

連續介質的質量守恒方程為

(5)

式中：u是連續介質的流速.設V是控制體體積，?V是控制體的邊界，n是邊界的單位外法向量.對式(5)關于控制體進行體積分得

(6)

分別采用雷諾輸運定理、散度定理對式(6)的第1和2項進行變換得

(7)

式中：ub是控制體邊界的移動速度；dS是控制體面上積分微元.將式(4)代入式(7)得

(8)

式中：Fj為單元面；j為單元面的索引；B為單元面總數.液相體積分數率定義為

(9)

將式(9)代入式(8)，并在等式兩端對時間積分，得到在該控制體處體積分數率的更新公式：

α1(t+Δt)=α1(t)-

(10)

(11)

圖1 網格的剛體運動和單元面掃過的體積Fig.1 The solid-body motion of the mesh and the volume swept by a cell face

圖2 面-界面交線運動修正Fig.2 The moving velocity corrections for face-interface intersection line

對于任意一個網格單元i，有如下等式近似成立：

(12)

下文中，經過運動通量修正和面-界面運動修正的 isoAdvector 方法稱為 isoAdvector-c2，僅有運動通量修正的 isoAdvector 方法稱為isoAdvector-c1.和幾何 VOF 方法不同，代數 VOF 方法不涉及任何關于幾何體的操作，而是求解關于體積分數率的輸運方程[13].

1.3 界面厚度估計

1.4 數值模擬的相關設置

由于計算域是封閉的空間，所有的邊界都是移動的壁面.采用浮力修正的k-ωSST 湍流模型[17]模擬湍流的影響.用 Crank-Nicolson 格式離散瞬態項，用高精度(High Resolution, HR)格式vanLeer 離散對流項，用Gauss linear 格式處理拉普拉斯算子的離散，用 PIMPLE 算法解耦求解 RANS 方程，用Gauss Seidel 方法求解線性方程組.求解壓力泊松方程時，GAMG(Generalized geometric-algebraic multi-grid) 和 PCG(Preconditioned conjugate gradient) 方法分別用于第1次和第2次壓力修正時方程組的求解.波高的提取步驟見文獻[11].

2 計算結果與討論

2.1 受迫非共振晃蕩

模型液艙長L=570 mm，寬S=310 mm，高D=300 mm[11].初始液位高度h0=150 mm.建立固結于大地的坐標系，初始時刻坐標系原點位于液艙底部中心，如圖3所示.波高監測點在xOy平面上的坐標是 (0, 0.265) m.

圖3 非共振晃蕩模型液艙和波高監測點布置(mm)Fig.3 The model tank and the wave probe of non-resonant sloshing (mm)

液艙受到簡諧強迫激勵：

y=-asin(ωt)

(13)

式中：a=0.005 m為激勵振幅，a/L=0.008 77；ω為受迫激勵頻率，ω/ω1=0.583，ω1為固有頻率[11].

在初始液面附近進行網格局部加密，最小網格尺寸約為 0.005 m.采用自適應步長進行計算，最大庫朗數是 0.2，初始時間步長為 0.000 1 s.分別采用3種 VOF 方法(isoAdvector-c2，isoAdvector-c1和MULES方法)捕捉自由液面的位置，并將結果與試驗結果[11]以及解析解[18]進行對比，如圖4所示，圖中η為波高.由圖4可見，對于遠離固有頻率的非共振晃蕩，波高振幅較小.3種VOF方法獲得的波高時歷均與解析解以及試驗結果較為吻合.

圖4 監測點處波高時歷對比Fig.4 The comparison of free-surface elevations at the probe

2.2 受迫共振晃蕩

模型試驗所用的液艙是邊長L=0.6 m 的立方體，建立固結于大地的坐標系，在初始時刻，坐標系的原點在液艙底部中心[19].初始液位高度h0=0.3 m.在 w3 處設置波高監測儀，如圖5所示.

圖5 共振晃蕩模型液艙和波高監測儀布置俯視圖(mm)Fig.5 Top view of the model tank and the wave probe of resonant sloshing (mm)

圖6 波高儀 w3 處波高的網格收斂性驗證 Fig.6 The verification of grid convergence for the free-surface elevations at wave probe w3

液艙受迫激勵，受迫運動為簡諧振蕩運動：

x=asin(ωt)

(14)

式中：a/L=0.008 71；ω/ω1,0=1.037，ω1,0是一階固有頻率[11].

表1 不同網格的基本信息Tab.1 The information about the different meshes

為進一步進行定量比較，引入數值計算結果和試驗數據間的方均根誤差E：

(15)

2.2.2自由液面位置對比 分別采用3種 VOF 方法(isoAdvector-c2、isoAdvector-c1和 MULES方法)捕捉自由液面的位置，并將結果與試驗結果[19]進行對比，如圖7所示.由圖可見：波峰的絕對值大于波谷的絕對值，表現為非線性特征；t<6 s時，各種 VOF 方法獲得的波高時歷和試驗值吻合較好；t>6 s時，isoAdvector-c1 方法得到的波高時歷和試驗值相差較遠，而isoAdvector-c2方法得到的波高時歷和試驗值較為吻合，這說明面-界面交線運動修正十分重要.3種方法得到的波高時歷與試驗結果的方均根誤差見表2.可見，由MULES得到的Eηw3比isoAdvector-c2方法得到的Eηw3大，說明isoAdvector-c2方法的計算結果和試驗值更為接近.MULES和isoAdvector-c2方法的界面厚度見表3，其中R為界面厚度與由MULES方法得到的界面厚度的比值.可以看出， isoAdvector-c2方法得到的界面厚度大約為MULES方法得到的界面厚度的1/2.

圖7 w3 波高儀處波高時歷Fig.7 The time histories of free-surface elevations at wave probe w3

表2 波高儀w3處波高和縱向力的方均根誤差

Tab.2 The RMSE of the free-surface elevations at wave probe w3 and the longitudinal force

方法Eηw3/cmEFx/NMULES2.07033.775isoAdvector-c21.51728.902isoAdvector-c15.01368.840

表3 界面厚度Tab.3 The interface thickness

圖8 晃蕩縱向力Fig.8 The longitudinal force

2.2.3整體水動力載荷對比 在整個液艙各個表面對熱力學壓力和黏性作用力進行積分，可得到晃蕩引起的整體水動力載荷.整體載荷在x方向的分量稱為縱向力Fx.上述3種VOF方法計算得到的Fx如圖8所示.可以看出：isoAdvector-c1 方法得到的Fx與試驗值不符，原因在于采用該方法捕捉的自由液面位置和試驗值偏差很大；而isoAdvector-c2 方法得到的縱向力和試驗值吻合較好.同樣引入方均根誤差來衡量MULES和isoAdvector-c2方法計算Fx的精度.將式(15)中的波高替換為Fx即得到縱向力方均根誤差EFx，結果見表2.可以看出，采用 isoAdvector-c2 方法計算得到的EFx比采用 MULES 方法計算得到的EFx更小.

2.3 單次沖擊波面

單次沖擊波面(SIW)[20]試驗的模型液艙如圖9所示，圖中縮尺比為1∶20.初始液位裝載率是20%，其一階共振周期T1=2.469 3 s.沿y軸進行特別的直線激勵，從而產生SIW.受迫激勵位移如圖10所示，中間一段強迫位移是1/2個簡諧波，對應的周期T=2.469 2 s，略小于T1.1/2個簡諧波激勵幅值a=244 mm,a/L=0.125.

圖9 SIW 模型液艙(mm)Fig.9 The SIW model tank (mm)

圖10 SIW 外激勵受迫位移時歷Fig.10 The time history of the forced displacement under the SIW excitation

對初始液面附近和相界面翻卷的區域進行局部網格加密，最小網格間距為 0.005 m.采用自適應步長計算，最大庫朗數為 0.2，初始時間步長設為 0.000 1 s. isoAdvector-c2和MULES 方法捕捉到的翻卷波形見圖11.試驗波面翻卷見文獻[20]中的圖16(b)，文獻[20]中的圖17(d)標出了波面形狀，波面之外的飛濺是 Kelvin-Helmholtz 相界面不穩定性造成的[20].由圖11可見，isoAdvector-c2和MULES方法捕捉得到的α1=0.5 等值線和試驗的波面等值線非常相似.用isoAdvector-c2方法得到的α1=0.001 等值線和α1=0.999 等值線的間距比MULES方法得到的間距小，這和表3中的結果一致，說明isoAdvector-c2捕捉的相界面更加明銳.

圖11 SIW 沖擊前波浪翻卷對比Fig.11 The comparison of the overturning waves before the impact

圖12 SIW 沖擊的三維視圖(m)Fig.12 The 3D views for the SIW (m)

波面沖擊壁面后引起破碎，沖擊前后的三維視圖對比如圖12所示.由圖12可見：在波峰后側，MULES捕捉得到的波面出現了褶皺，這是代數VOF方法應用壓縮格式保證界面明銳所帶來的負面影響；而isoAdvector-c2方法捕捉到的波面更為平滑，波峰后側沒有出現褶皺.

3 結論

基于幾何VOF方法 isoAdvector引入運動通量修正，提出面-界面交線運動修正，并采用一種近似公式，根據各個網格單元面的運動通量重構出各個單元體心的剛體運動速度，得到如下結論：

(1) 在模擬受迫非共振晃蕩時，MULES、isoAdvector-c1和 isoAdvector-c2方法計算得到的波高均與試驗結果和解析解吻合較好.在模擬受迫共振晃蕩時，相比于 isoAdvector-c1，isoAdvector-c2方法因增加了面-界面交線運動修正，得到的波高和縱向力的結果與試驗值吻合得更好.

(2) 相比代數 VOF 方法 MULES，采用 isoAdvector-c2 方法得到的波高和縱向力相對于試驗值的方均根誤差更小.

(3) 采用修正的 isoAdvector 方法捕捉到的兩相界面的厚度大約是代數 VOF 方法 MULES所得界面厚度的1/2，兩相界面更為明銳.

(4) 在模擬單次沖擊波面時，isoAdvector-c2和 MULES 方法均能捕捉到和試驗波面相似的波面形狀，能很好地模擬波面的翻卷和破碎；基于 MULES 方法捕捉的波面在波峰后側出現褶皺，而 isoAdvector-c2 方法捕捉到的波面沒有褶皺，更加平滑.