VOF與Level Set耦合的界面計算方法及其對波浪破碎過程的模擬

賈 偉，柳淑學，李金宣，范玉平

(大連理工大學 海岸和近海工程國家重點實驗室，大連 116024)

近年來，在海洋工程流體流動和波浪傳播過程的數值模擬中，經常采用兩相流模型(Two-Phase Fluid Model)，即將空氣和水視為兩種流動參數不同，但滿足相同控制方程的流體。該模型采用同一組控制方程同時求解兩相流體，并采用界面計算方法來捕捉兩相界面，因此也被稱為擬單流體法。界面計算方法在兩相流模擬中至關重要，VOF方法和Level Set方法是目前兩相流模型中較為常用的兩種界面處理方法。根據計算方法不同，VOF方法又分為代數方法和幾何方法。

最為常見的VOF代數方法是界面壓縮法，例如開源代碼OpenFOAM中采用的MULES(Multidimensional Universal Limiter with Explicit Solution)方法[1]。該方法通過限制器來混合一階與高階格式，并且在對流方程中添加壓縮項來保持界面的尖銳性。代數方法精確度相對較低，但是由于沒有幾何計算以及界面重構步驟，因此計算效率高。幾何方法計算較為精確，但是計算量大，并且往往對網格類型有所限制，其中最常見的是由Youngs[2]提出的PLIC(Piecewise Linear Interface Construction)方法。

VOF方法的優點是可以保持質量守恒，但由于采用了體積分數來表示流體在體積單元的占比，因此體積分數場不連續，無法精確計算界面處的曲率和法向量，也無法精確計算界面處的表面張力。Level Set方法采用Level Set函數表示某點與界面的距離，曲率和法向量的計算更為準確。但該方法質量不守恒，當網格較稀疏時尤為明顯。

因此，一些學者將VOF方法與Level Set方法進行耦合，從而利用兩種方法各自的優點并規避缺點。由于側重點和計算方法各不相同，各耦合方法也有較大差異。Sussman等[3]采用二階算子分裂算法，已知某時刻的體積分數、Level Set函數和速度場，根據該時刻各個單元體積面的通量求解下一時刻的體積分數場和Level Set函數。相比于Level Set方法，對于界面有尖角時，耦合方法計算效果好。Pijl等[4]通過體積分數與Level Set函數之間的代數關系，提出了二者之間簡化的耦合方法，從而避免了在每個時間步內同時計算體積分數和Level Set函數的對流方程。Sun等[5]提出了耦合的VOSET方法，該方法采用PLIC方法計算體積分數和重構界面，并提出了一種迭代的幾何計算方法，用于計算界面附近的Level Set函數，從而計算界面的曲率。Albadawi等[6]采用Level Set函數計算表面張力和界面曲率，通過耦合方法模擬了液體中的氣泡，并與實驗結果進行了對比。該方法適用于計算常接觸角問題，主要對表面張力進行了修正。Cao等[7]采用多維對流算法求解體積分數，采用類似Sun等的方法幾何迭代求解Level Set函數并進行耦合，該方法可用于非結構化網格。

在現有的數值波浪計算模型中，通常采用單一的VOF方法或Level Set方法計算界面，即存在無法精確計算界面法向量或無法保持質量守恒的缺點。此外，現有的大部分VOF與Level Set耦合方法不適用于非結構化網格，并且計算量較大。因此，本文建立了基于VOF與Level Set方法的耦合方法，在保證計算精度的同時提高計算效率。該方法可用于結構化網格及非結構化網格，并易于擴展到三維。隨后，采用該耦合方法計算了潰壩問題，驗證了方法的有效性，并進一步應用于波浪傳播過程中破碎問題，與實驗結果進行了對比，表明該方法可以較好地應用于波浪破碎過程的模擬。

1 數值方法

1.1 兩相流模型與界面處理方法

本文的數值模型基于OpenFOAM中的不可壓縮兩相流模型求解器interFoam進行開發，其控制方程為連續性方程以及N-S方程

·U=0

(1)

(2)

在interFoam求解器中采用VOF方法計算界面，即采用體積分數α作為控制參數，兩相交界處的流體性質參數可以表示為

ρ=αρw+(1-α)ρa

(3)

μ=αμw+(1-α)μa

(4)

式中：ρw、ρa分別為液體、氣體的密度，μw、μa分別為液體、氣體的動力粘滯系數。體積分數α滿足相分數方程

(5)

界面處的法向量可以表示為

(6)

在本文模型中，與Level Set方法進行了耦合。Level Set方法采用Level Set函數φ(x,t)表示t時刻x位置與界面的距離，在界面處函數值為0。用Γ(t)表示界面位置，則0時刻φ(x,0)的表達式為

(7)

其中d是與界面的距離。此外，在Level Set方法中，通常采用Heaviside函數光滑界面附近的密度、動力粘滯系數等不連續物理量，即

ρ(φ)=H(φ)ρw+(1-H(φ))ρa

(8)

μ(φ)=H(φ)μw+(1-H(φ))μa

(9)

其中H(φ)為Heaviside函數

(10)

式中：φ為Level Set函數，b為界面厚度，通常取b=1.5△x，其中△x為網格尺寸。

與體積分數相似，Level Set函數需要滿足對流方程

(11)

需要注意的是，通常在采用Level Set方法計算界面時，由式(11)計算所得的函數值并不能保證與該時刻點x和界面的距離相等，因此在每個時間步前需要對Level Set函數重新初始化[3]。本文采用了幾何迭代方法來計算Level Set函數，避免了求解對流方程及初始化過程。

由于Level Set函數是一個連續函數，因此可以較為精確地計算各點的法向量和曲率，即

(12)

(13)

在本文模型中，曲率采用式(13)計算，則控制方程式(2)可以表示為

(14)

其中ρ和μ等不連續物理量采用式(8)和式(9)進行計算。

1.2 耦合算法

本文采用幾何方法求解體積分數從而保證計算精度，只在界面附近通過幾何關系初始化Level Set函數，并迭代求解，從而避免求解Level Set函數的對流方程以及重新初始化方程，提高計算效率。該算法分為3個部分：求解體積分數的對流方程、界面重構和構造Level Set函數。

1.2.1 求解體積分數的對流方程

本文在計算體積分數的對流方程時，采用了isoAdvector方法，詳細內容可以參考Roenby等人的論文[8]，其開源代碼已植入到OpenFOAM-v1612+及之后的版本中。該方法的簡要步驟如下：

(1)尋找界面單元，即體積分數介于0與1之間的單元。

(2)對于界面單元的頂點，用共享該頂點的所有網格單元內的體積分數進行加權平均，權重系數為該節點到單元中心距離的倒數。該步驟可以將體積分數插值到頂點，假設某界面單元有N個頂點，其坐標分別表示為x1，x2，…，xN，相應的體積分數插值結果為f1，f2，…，fN。

(3)根據各頂點的f值，可以在單元內部通過線性插值構建關于f的等值線(二維問題)或等值面(三維問題)。如果該單元的某個邊(xm，xn)滿足fm

(15)

找到每條邊的交點并連接，即可得到等值面。

(4)找到等值面，使得該等值面按照該單元的體積分數αi將單元分成2個部分，即求f*，使α(f*)=αi。

首先根據幾何關系計算f1，f2，…，fN對應的體積分數值α(f1)，α(f2)，…，α(fN)。根據該單元的體積分數αi，在α(f1)，α(f2)……α(fN)中尋找最接近αi的兩個值α(fm)和α(fn)，使得α(fm)≤αi≤α(fn)。在此區間內，α與f之間的關系為三次多項式，可以在此區間內通過線性插值構造兩個等值面，由幾何關系求得相應的α，即可求解4×4矩陣得出多項式的系數，進而確定f*的值，并得到目標等值面。

(5)根據等值面在當前時間步的運動，即可通過幾何方法求出該單元下一時間步的體積分數。

1.2.2 界面重構

圖1 某網格單元內的等值線示意圖

需要注意的是，上述isoAdvector方法構建的等值面僅用于計算體積分數，并非真實的兩相界面。由于在構造Level Set函數時需要用到界面的幾何信息，而isoAdvector并未包含界面重構的步驟，因此在本模型中建立了簡便的重構算法：

(1)按照式(6)，根據體積分數場計算精度較低的法向量n′。

(2)對于頂點坐標為x1，x2，…，xN的界面單元，令d′ =n′·x，即可在單元內構建關于d′的等值面，且等值面均與法向量n′垂直，對于二維網格如圖2所示。

(3)此時，可以通過1.2.1(4)中的方法計算關于d′的等值面，使得α(d′)=αi。由于該界面與法向量垂直，且將該單元劃分為體積分數為αi與1-αi兩個部分，因此該等值面即該單元內的界面。

圖3 算法流程圖

1.2.3 構造Level Set函數

參考Sun 等[5]，Cao等[7]，Level Set與VOF的耦合過程采用了幾何迭代求解的方法。其簡要步驟如下：

(1)標記界面單元附近3△x區域內的體積單元，其中△x為網格尺寸。

(2)計算被標記區域內的Level Set函數。通過幾何關系，計算網格點與1.2.2節所得重構界面的距離d，即可按照式(7)構造Level Set函數。

(3)按照式(12)計算較為精確的法向量，從而按照1.2.2節的步驟重新構造更為精確的界面。構造Level Set函數和重構界面的過程可以迭代兩次，從而得到最終的界面。

整體而言，上述方法可以視為一種改進的VOSET方法。在計算體積分數時，原始的VOSET方法采用了PLIC方法，該方法適用于結構化網格，拓展到非結構化網格時較為復雜，因此本文采用了isoAdvector算法，該算法計算速度較快，并且適用于各類網格；之后參考isoAdvector中計算等值面的方法，建立了重構界面的簡便算法；采用VOSET方法中的迭代幾何方法計算Level Set函數，從而精確計算界面處的法向量及曲率，避免了求解Level Set函數的對流方程及初始化方程。該算法的計算流程如圖3所示。

2 數值模擬驗證

圖4 潰壩算例的布置圖

Fig.4 Layout of dam break problem

圖5 水體前端的無因次位置與無因次時間的關系

Fig.5 Relationship between dimensionless locations of water front and dimensionless time

潰壩問題是經典的兩相流算例，因此通過模擬該問題，與Martin 等[9]的實驗結果進行對比。圖4為潰壩算例的布置圖，整個計算域是邊長為4L的正方形，其中L=0.146 m。在計算域左下方有一塊寬為L，高為2L的靜止水體，當計算開始后，水體由于重力作用開始流動。水的密度為1.0×103kg/m3，表面張力系數為0.075 5 N/m，重力加速度為9.8 m/s2。

在數值模擬時，采用均勻的正方形網格，共100×100個。采用圖4所示的xOy坐標系，將運動水體的最前端與y軸的距離用x表示。將x和時間t無因次化，即令X*=x/L,t*=t·(2g/L)0.5，則計算所得X*與t*的關系如圖5所示。圖中同時給出了Martin 等[9]的實驗結果以及Sun 等[5]和Cao等[7]的VOSET模型計算結果。從圖中可以看出，本文模型的計算結果與實驗結果及文獻中VOSET模型的計算結果均較為接近，因此可以說明本文的界面計算方法是有效的。作為示例，圖6給出了潰壩過程中，典型時刻的瞬時水體自由表面圖。

6-at=0.0 s 6-bt=0.075 s 6-ct=0.15 s 6-dt=0.225 s 6-et=0.30 s 6-ft=0.375 s

圖6 不同時刻的水體自由表面圖

Fig.6 Snapshots of the water interfaces at different times

3 波浪傳播破碎問題的模擬

波浪破碎是波浪傳播過程中的重要現象，由于波浪破碎過程的復雜性和隨機性，其界面往往較為復雜。本文模擬了規則波在斜坡地形及島礁地形上的破碎問題，并與實驗結果進行了對比。

3.1 波浪在斜坡上的破碎

圖7 波浪在斜坡上破碎實驗布置圖(單位：m)

該實驗在大連理工大學海岸和近海工程國家重點實驗室的三維水池中進行[10]，水池長度40 m，寬度24 m，實驗水深0.4 m，地形為1:15的斜坡，坡腳與造波機的距離為10 m，坡頂高度0.4 m，從而保證波浪破碎發生在斜坡上，浪高儀及地形布置如圖7，波浪參數如表1所示。由于數值計算中僅采用了正向規則波，因此簡化為二維問題進行計算。

在數值波浪水槽中，計算域總長度為20 m，高度0.57 m，在水槽邊界采用速度入口法[11]進行造波。整個計算域內采用尺寸相同的矩形網格，網格尺寸為△x=0.01 m，△y=0.007 5 m，在坡面與網格交界處應用OpenFOAM中的snappyHexMesh方法，將網格切割從而使網格與坡面貼合。湍流采用k-ωSST模型進行模擬。表1為模擬波浪的參數。

表1 斜坡地形波浪參數

圖8 算例1的波面歷時曲線

針對算例1所模擬的波面歷時曲線與試驗測量結果的對比如圖8所示，圖中右上角標注數字是圖7中的浪高儀編號，圖中虛線為實驗結果，實線為數值計算結果。1號浪高儀位于平底區域，可以看出數值模型中的入射波浪與實驗入射波浪基本一致。由4號浪高儀到13號浪高儀，波浪受淺化作用影響，波浪的非線性及波面的不對稱性逐漸增強，13號浪高儀處為波浪的破碎點，之后進入破碎區，波高迅速減小。由實驗觀察可知，該波浪為卷破波，波面翻卷之后拍打坡面，造成液滴飛濺，因此在破碎區內，波浪的紊流特性較為明顯，波面隨機性極強，同時破碎波浪存在的浪花及水氣混合使得波浪測量結果存在一定的誤差。整體而言，該數值模型能夠較為準確地模擬波浪淺化過程，也能夠較為準確地預測波浪破碎位置及破碎波高。

圖9顯示了算例1在不同時刻的水面變化情況，圖中灰色區域為空氣，黑色區域為水體。根據實驗觀察，13號浪高儀為波浪的破碎點(即圖中的垂線所在位置，距離造波機14.25 m處)，即卷破波傳播到該點時，波前近乎為垂直面。由圖示模擬結果可見，卷破波的波前在t=29.9 s時幾乎垂直，而在t=29.92 s時已微微前卷，說明數值模型所預測的破碎點與實驗結果一致。此外，圖9展示了卷破波的波前翻卷，與水體接觸之后飛濺(splash-up)的過程，說明該界面計算方法能夠較為準確地計算波浪破碎過程。

9-at=29.82 s 9-bt=29.9 s 9-ct=29.92 s

9-dt=30 s 9-et=30.08 s 9-ft=30.16 s

注：圖中垂線的橫坐標為實驗觀察所得的破碎點橫坐標。

圖9 模擬算例1的波浪破碎過程

Fig.9 Snapshots of breaking process in case 1

算例2的計算結果如圖10所示。其入射波高與第1組基本相同，但相較于第1組波浪周期增大，波長變大且水深不變，因此波浪非線性明顯增強。在坡腳附近，由于波浪的非線性的影響，波浪已經存在較為明顯的變形。之后波浪在斜坡地形上發生淺化變形，直到13號浪高儀處發生破碎。與第1組波浪類似，本模型較好地計算了波浪的破碎過程，包括破碎位置、破碎波高等與實驗結果基本一致。

3.2 波浪在島礁地形上的破碎

圖10 算例2的波面歷時曲線

實驗在大連理工大學海岸和近海工程國家重點實驗室的波流水槽內進行[12]，實驗的水槽長度為69 m，寬0.8 m，深1.8 m，地形及浪高儀布置如圖11所示。島礁地形簡化為坡度為1:5的斜坡，坡頂高度為0.5 m，水深0.715 m。模擬波浪的周期為1.6 s，其他參數如表2所示。與3.1部分的算例不同，該地形為簡化的島礁地形，波浪在經過斜坡地形之后進入頂部的礁坪段，即波浪在淺化變形之后再在平面上傳播，地形對波浪的影響更為復雜。選取的2組算例均在波浪傳播到礁坪之后發生破碎，分別為崩破波和卷破波。

表2 島礁地形波浪實驗參數

圖11 波浪在簡化島礁地形上破碎實驗布置圖(單位：m)

圖12 算例3的波面歷時曲線

在數值波浪水槽中，計算域總長度為45 m，高度0.9 m，在水槽邊界采用速度入口法進行造波。整個計算域內的初始網格均為尺寸相同的正方形，△x=△y=0.012 m，在坡面與網格交界處采用SnappyHexMesh方法，將網格切割從而使網格與坡面貼合。

算例3的數值波面歷時曲線與實驗結果的對比如圖12所示。可以看出，數值模擬結果與實驗結果吻合較好。在圖中1號、3號浪高儀處，波浪在平底上傳播，波形較為穩定；在6號、8號、10號浪高儀處，波浪在島礁斜坡上傳播，由于水深快速變淺，波浪淺化作用明顯，波高逐漸增大，波浪的非線性、波面的不對稱性逐漸增強；12號浪高儀處于坡頂的礁坪段，波高繼續增大，并最終在13號浪高儀處發生波浪破碎；14號、15號、16號浪高儀位于波浪破碎點之后，該區域內湍流極強，波面的隨機性也較強，在該區域內數值結果與實驗結果存在一些誤差，但整體上基本一致。

圖13 算例4的波面歷時曲線

算例4的計算結果如圖13所示。可以看出，與之前組次類似，數值結果與實驗結果整體上吻合較好，可以較為準確地計算波浪的淺化作用(9號)、破碎點(12號)的位置以及破碎波高。總之，根據以上4個算例，采用本文的界面計算方法可以較為準確地模擬在復雜斜坡上的波浪傳播破碎問題。

4 結論

本文建立了一種VOF方法與Level Set方法耦合的界面計算方法，其中體積分數采用幾何方法——isoAdvector方法進行計算，可保證計算精度；建立了簡化的界面重構算法，從而便于采用幾何方法建立Level Set函數；通過迭代方法計算更為準確的Level Set函數及法向量。該算法是VOSET方法的一種改進，保持VOF方法及Level Set方法各自的優點，并盡可能地保證算法的簡潔以及非結構化網格的適用性。此外，雖然本文僅考慮了二維模型，但該方法的各個步驟均易于推廣到三維。

基于潰壩問題的求解，證明了該計算方法的有效性。并進一步應用于波浪在斜坡地形、島礁地形上的傳播破碎問題，數值模擬結果與實驗結果吻合較好，說明該方法可以很好地模擬復雜地形上的波浪傳播破碎問題。