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

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

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

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

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

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

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

1 數(shù)值方法

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

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

·U=0

(1)

(2)

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

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

(3)

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

(4)

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

(5)

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

(6)

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

(7)

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

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

(8)

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

(9)

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

(10)

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

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

(11)

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

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

(12)

(13)

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

(14)

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

1.2 耦合算法

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

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

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

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

(2)對(duì)于界面單元的頂點(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個(gè)頂點(diǎn)，其坐標(biāo)分別表示為x1，x2，…，xN，相應(yīng)的體積分?jǐn)?shù)插值結(jié)果為f1，f2，…，fN。

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

(15)

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

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

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

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

1.2.2 界面重構(gòu)

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

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

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

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

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

圖3 算法流程圖

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

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

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

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

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

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

2 數(shù)值模擬驗(yàn)證

圖4 潰壩算例的布置圖

Fig.4 Layout of dam break problem

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

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

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

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

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 不同時(shí)刻的水體自由表面圖

Fig.6 Snapshots of the water interfaces at different times

3 波浪傳播破碎問(wèn)題的模擬

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

3.1 波浪在斜坡上的破碎

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

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

在數(shù)值波浪水槽中，計(jì)算域總長(zhǎng)度為20 m，高度0.57 m，在水槽邊界采用速度入口法[11]進(jìn)行造波。整個(gè)計(jì)算域內(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的波面歷時(shí)曲線

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

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

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í)驗(yàn)觀察所得的破碎點(diǎn)橫坐標(biāo)。

圖9 模擬算例1的波浪破碎過(guò)程

Fig.9 Snapshots of breaking process in case 1

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

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

圖10 算例2的波面歷時(shí)曲線

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

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

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

圖12 算例3的波面歷時(shí)曲線

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

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

圖13 算例4的波面歷時(shí)曲線

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

4 結(jié)論

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

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