賈 偉，柳淑學(xué)，李金宣，范玉平

(大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗室，大連 116024)

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

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

VOF方法的優(yōu)點(diǎn)是可以保持質(zhì)量守恒，但由于采用了體積分?jǐn)?shù)來表示流體在體積單元的占比，因此體積分?jǐn)?shù)場不連續(xù)，無法精確計算界面處的曲率和法向量，也無法精確計算界面處的表面張力。Level Set方法采用Level Set函數(shù)表示某點(diǎn)與界面的距離，曲率和法向量的計算更為準(zhǔn)確。但該方法質(zhì)量不守恒，當(dāng)網(wǎng)格較稀疏時尤為明顯。

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

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

1 數(shù)值方法

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

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

·U=0

(1)

(2)

在interFoam求解器中采用VOF方法計算界面，即采用體積分?jǐn)?shù)α作為控制參數(shù)，兩相交界處的流體性質(zhì)參數(shù)可以表示為

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

(3)

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

(4)

式中：ρw、ρa(bǔ)分別為液體、氣體的密度，μw、μa分別為液體、氣體的動力粘滯系數(shù)。體積分?jǐn)?shù)α滿足相分?jǐn)?shù)方程

(5)

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

(6)

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

(7)

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

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

(8)

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

(9)

其中H(φ)為Heaviside函數(shù)

(10)

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

與體積分?jǐn)?shù)相似，Level Set函數(shù)需要滿足對流方程

(11)

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

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

(12)

(13)

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

(14)

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

1.2 耦合算法

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

1.2.1 求解體積分?jǐn)?shù)的對流方程

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

(1)尋找界面單元，即體積分?jǐn)?shù)介于0與1之間的單元。

(2)對于界面單元的頂點(diǎn)，用共享該頂點(diǎn)的所有網(wǎng)格單元內(nèi)的體積分?jǐn)?shù)進(jìn)行加權(quán)平均，權(quán)重系數(shù)為該節(jié)點(diǎn)到單元中心距離的倒數(shù)。該步驟可以將體積分?jǐn)?shù)插值到頂點(diǎn)，假設(shè)某界面單元有N個頂點(diǎn)，其坐標(biāo)分別表示為x1，x2，…，xN，相應(yīng)的體積分?jǐn)?shù)插值結(jié)果為f1，f2，…，fN。

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

(15)

找到每條邊的交點(diǎn)并連接，即可得到等值面。

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

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

(5)根據(jù)等值面在當(dāng)前時間步的運(yùn)動，即可通過幾何方法求出該單元下一時間步的體積分?jǐn)?shù)。

1.2.2 界面重構(gòu)

圖1 某網(wǎng)格單元內(nèi)的等值線示意圖

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

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

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

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

圖3 算法流程圖

1.2.3 構(gòu)造Level Set函數(shù)

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

(1)標(biāo)記界面單元附近3△x區(qū)域內(nèi)的體積單元，其中△x為網(wǎng)格尺寸。

(2)計算被標(biāo)記區(qū)域內(nèi)的Level Set函數(shù)。通過幾何關(guān)系，計算網(wǎng)格點(diǎn)與1.2.2節(jié)所得重構(gòu)界面的距離d，即可按照式(7)構(gòu)造Level Set函數(shù)。

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

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

2 數(shù)值模擬驗證

圖4 潰壩算例的布置圖

Fig.4 Layout of dam break problem

圖5 水體前端的無因次位置與無因次時間的關(guān)系

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

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

在數(shù)值模擬時，采用均勻的正方形網(wǎng)格，共100×100個。采用圖4所示的xOy坐標(biāo)系，將運(yùn)動水體的最前端與y軸的距離用x表示。將x和時間t無因次化，即令X*=x/L,t*=t·(2g/L)0.5，則計算所得X*與t*的關(guān)系如圖5所示。圖中同時給出了Martin 等[9]的實(shí)驗結(jié)果以及Sun 等[5]和Cao等[7]的VOSET模型計算結(jié)果。從圖中可以看出，本文模型的計算結(jié)果與實(shí)驗結(jié)果及文獻(xiàn)中VOSET模型的計算結(jié)果均較為接近，因此可以說明本文的界面計算方法是有效的。作為示例，圖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 波浪傳播破碎問題的模擬

波浪破碎是波浪傳播過程中的重要現(xiàn)象，由于波浪破碎過程的復(fù)雜性和隨機(jī)性，其界面往往較為復(fù)雜。本文模擬了規(guī)則波在斜坡地形及島礁地形上的破碎問題，并與實(shí)驗結(jié)果進(jìn)行了對比。

3.1 波浪在斜坡上的破碎

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

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

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

表1 斜坡地形波浪參數(shù)

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

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

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

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

注：圖中垂線的橫坐標(biāo)為實(shí)驗觀察所得的破碎點(diǎn)橫坐標(biāo)。

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

Fig.9 Snapshots of breaking process in case 1

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

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

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

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

表2 島礁地形波浪實(shí)驗參數(shù)

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

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

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

算例3的數(shù)值波面歷時曲線與實(shí)驗結(jié)果的對比如圖12所示?？梢钥闯觯瑪?shù)值模擬結(jié)果與實(shí)驗結(jié)果吻合較好。在圖中1號、3號浪高儀處，波浪在平底上傳播，波形較為穩(wěn)定；在6號、8號、10號浪高儀處，波浪在島礁斜坡上傳播，由于水深快速變淺，波浪淺化作用明顯，波高逐漸增大，波浪的非線性、波面的不對稱性逐漸增強(qiáng)；12號浪高儀處于坡頂?shù)慕钙憾?，波高繼續(xù)增大，并最終在13號浪高儀處發(fā)生波浪破碎；14號、15號、16號浪高儀位于波浪破碎點(diǎn)之后，該區(qū)域內(nèi)湍流極強(qiáng)，波面的隨機(jī)性也較強(qiáng)，在該區(qū)域內(nèi)數(shù)值結(jié)果與實(shí)驗結(jié)果存在一些誤差，但整體上基本一致。

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

算例4的計算結(jié)果如圖13所示?？梢钥闯?，與之前組次類似，數(shù)值結(jié)果與實(shí)驗結(jié)果整體上吻合較好，可以較為準(zhǔn)確地計算波浪的淺化作用(9號)、破碎點(diǎn)(12號)的位置以及破碎波高?？傊?，根據(jù)以上4個算例，采用本文的界面計算方法可以較為準(zhǔn)確地模擬在復(fù)雜斜坡上的波浪傳播破碎問題。

4 結(jié)論

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

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