鄒學(xué)鋒 朱良生 張善舉

(華南理工大學(xué)土木與交通學(xué)院，廣東廣州510640)

珊瑚礁廣泛分布于太平洋和印度洋的亞熱帶地區(qū)，其中也分布于我國南海地區(qū)。珊瑚礁的特點(diǎn)是礁坪水深很淺，常伴有瀉湖、凹地，并且礁前斜坡陡峭。近年來，為了開發(fā)南海的自然資源和提高南海的經(jīng)濟(jì)價(jià)值，我國在南海諸島上建立了多個(gè)重要港口，它們多建立在珊瑚礁地形上，按地形分有瀉湖港、礁緣港、島緣港、連島港和人工島港等，我國在南沙群島上施工的幾處較大的人工島都包含港口［1］。在島礁區(qū)周圍、島礁之間都已經(jīng)形成了一些習(xí)慣性的航道。波浪是港口工程里常見的水動(dòng)力現(xiàn)象，當(dāng)波浪傳播到礁面時(shí)，波浪在斜坡上發(fā)生淺水變形，波高增大，影響周圍水域內(nèi)船舶的航行。另外，波浪在礁緣或者礁坪上破碎，發(fā)生增水，對(duì)礁坪后方的港口建筑物產(chǎn)生威脅。因此，研究島礁的水動(dòng)力特征，對(duì)港口的建設(shè)和船舶的航行都有著十分重要的意義。

波浪破碎一直是港口工程領(lǐng)域內(nèi)的重要課題。一方面，波浪破碎會(huì)給海床、海岸建筑物施加額外的作用力，帶動(dòng)泥沙輸運(yùn)，改變著海岸線的面貌;另一方面，波浪破碎后引起的增水和波高變化對(duì)堤岸和港口的設(shè)計(jì)與建造有著主要的影響。波浪破碎的地方稱為破波區(qū)［2］。波浪在破碎前的形態(tài)發(fā)生劇烈變化，并可能在破波區(qū)發(fā)生多次破碎，最終大部分能量轉(zhuǎn)化為湍流動(dòng)能。

Galvin［3］將波浪破碎分為3種類型:崩破波、卷破波、激破波。Battjes［4］用 Iribarren數(shù) ξ0來判斷波浪破碎類型，ξ0＜0.5時(shí)為崩破波，0.5＜ξ0＜3.3時(shí)為卷破波，ξ0＞3.3時(shí)為激破波。在波浪破碎過程中，波浪開始形成渦旋，湍流的產(chǎn)生是其中的控制性因素之一。關(guān)于崩破波和卷破波的實(shí)驗(yàn)和數(shù)學(xué)模型一直以來都在不斷發(fā)展。任冰等［5］通過物理實(shí)驗(yàn)?zāi)M了規(guī)則波在陡峭珊瑚礁地形上的傳播和破碎變化過程;李紹武等［6］通過物理實(shí)驗(yàn)研究了不規(guī)則波況下岸礁礁坪上的增水;諸裕良等［7］通過物理實(shí)驗(yàn)?zāi)M了復(fù)合坡度下珊瑚礁地形上的波浪破碎。另一方面，數(shù)學(xué)模型能提供波浪破碎時(shí)的具體水動(dòng)力學(xué)變量值，有些在實(shí)驗(yàn)中不易獲取，如流速和湍流動(dòng)能分布等。在模擬波浪破碎的數(shù)學(xué)模型中，主要有Bousniessq類模型，基于納維埃-斯托克斯方程的非靜壓模型和VOF模型，Bousniessq方程是垂向積分方程，模擬波浪破碎時(shí)需要在動(dòng)量方程中加入一個(gè)耗散項(xiàng)來將破碎過程參數(shù)化，但是參數(shù)不易選?。?］。而非靜壓模型不能處理破波區(qū)的摻氣作用和水體分離，VOF模型相較于前兩者能較好地模擬考慮空氣作用時(shí)的波浪破碎，以及波浪破碎時(shí)水氣界面上的湍流分布。Zhang等［9］使用改進(jìn)的Boussinesq方程模擬了極陡珊瑚礁地形上的波浪傳播;房克照等［10］將非靜壓模型分解為靜壓步和非靜壓步兩個(gè)過程混合求解，處理了海岸區(qū)域常見的波浪破碎和水-陸動(dòng)邊界問題。在以往的波浪破碎模型中一些學(xué)者采用了單相流模型［2，8，10-11］，但是也有些學(xué)者認(rèn)為單相流模型無法考慮波浪破碎中卷入空氣的能量耗散，會(huì)造成湍流動(dòng)能整體偏大［12-14］。Xie［15］在兩相流中采用了 k-ε 模型模擬了崩破波和卷破波;Jacobsen等［16］在兩相流中使用k-ω模擬了波浪破碎。k-ω SST模擬存在反向壓力梯度的流動(dòng)場(chǎng)有較好的表現(xiàn)［17］。不過，在以往的波浪破碎模型中，使用標(biāo)準(zhǔn)的湍流模型會(huì)造成波浪在水面處的過度衰減，本文采用k-ω SST湍流模型，并按照Brown等［18］的方法建立了兩相湍流模型，模擬了波浪在珊瑚礁水域的波高、水位和湍流動(dòng)能的變化規(guī)律，數(shù)學(xué)模型使用Yao等［19］的物模結(jié)果進(jìn)行驗(yàn)證。

1 數(shù)值模型

1.1 控制方程

本文采用OpenFOAM中interFOAM的兩相流求解器進(jìn)行模擬，并使用Higuera等［20］的速度邊界造波和主動(dòng)消波邊界來實(shí)現(xiàn)入口的造波和出口消波?？刂品匠倘缦?

式中:U為速度矢量;t為時(shí)間;p*為動(dòng)壓力;ρ為計(jì)算單元體積密度;g為重力加速度矢量;X=(x，y，z)，為位置矢量;σ為表面張力系數(shù);κ=，為界面的平均曲度;α為體積分?jǐn)?shù)函數(shù);μeff為有效動(dòng)力黏度，μeff=μ +ρνt，μ 為分子動(dòng)力黏度，νt為湍流運(yùn)動(dòng)黏度，由所選的湍流模型決定。流體體積方程為

其中，α為計(jì)算單元中水的體積分?jǐn)?shù)，是個(gè)有界值，即0≤α≤1。

在OpenFOAM里沒用壓縮差分格式計(jì)算自由表面，而是加入了人為壓縮項(xiàng)·［Urα (1-α)］，使得空氣與水的交界面保持尖銳，因此修改后的α方程變?yōu)?/p>

式中，Ur=min［cαU，max(U)］，與交界面垂直，cα為自定義壓縮系數(shù)，決定了空氣和水體交界面的壓縮程度，增大cα的數(shù)值可以加大交界面的壓縮性，使交界面更加尖銳，但同時(shí)也會(huì)產(chǎn)生更大的速度梯度，不利于數(shù)值穩(wěn)定［13］，這里取默認(rèn)值cα=1。體積分?jǐn)?shù)α方程的有界性通過MULES(多維通用顯式求解限制器)來完成，它通過離散散度項(xiàng)中的限制因子來確保α最后的值位于0到1之間。

1.2 湍流模型

在這里，使用標(biāo)準(zhǔn)k-ω SST模型模擬了波浪在珊瑚礁上的破碎，并根據(jù)Brown等［18］的建議將單元格密度顯式地加入到方程中，建立了兩相湍流模型。

標(biāo)準(zhǔn)k-ω SST模型的方程如下:

式中:k為湍流動(dòng)能;ν為分子運(yùn)動(dòng)黏度;ω為湍流內(nèi)能轉(zhuǎn)化率;Pk為湍流動(dòng)能源項(xiàng)，其定義見本文2.1節(jié);b1=，S為流體應(yīng)變率，F(xiàn)2為混合函數(shù)［17-18］;常系數(shù)β*=0.09;σω、σk、β、γ是混合變量，由混合表達(dá)式 (7)計(jì)算，默認(rèn)值見表1，

F1為混合函數(shù)，詳細(xì)定義參見文獻(xiàn) ［17-18］。

表1 混合函數(shù)的默認(rèn)值Table 1 Default values for blender functions

進(jìn)行顯式密度后的k-ω SST模型的方程如下:

1.3 邊界條件

左側(cè)邊界采用速度邊界入口，湍流動(dòng)能入口邊界采用Lin等［8］的方法，即k=(cI)2，其中，c和I分別為波形傳播速度和湍流強(qiáng)度。計(jì)算區(qū)域內(nèi)部的湍流動(dòng)能初始值與入口處相同，入口邊界湍流內(nèi)能轉(zhuǎn)化率，其中，這里采用Lin等［8］的取值0.0025，底部邊界的k和ω使用壁面函數(shù)［17］。底部邊界的速度設(shè)置為Dirichlet邊界條件，壓力和體積分?jǐn)?shù)為Neumann邊界條件。頂部的速度、壓力和體積分?jǐn)?shù)都為Dirichlet-Neumann混合邊界條件。邊界條件的具體設(shè)置詳見表2。

1.3.1 造波邊界

這里使用Higuera等［20-21］的速度邊界造波法，在入口處使用速度邊界造波，并能主動(dòng)吸收從島礁反射回來的波浪。在OpenFOAM中，為了方便造波，假定重力方向?yàn)閆軸負(fù)方向，同一個(gè)造波邊界面的最低點(diǎn)處在同一水平線上，但不同的造波邊界面可以在不同水平線上。通過不同的波浪理論，能夠計(jì)算出任意坐標(biāo)位置上自由水面的高度和質(zhì)點(diǎn)速度，稱為Dirichlet邊界條件，而壓力采用Open-FOAM里的NeuMann邊界條件來計(jì)算。1.3.2 消波邊界

表2 OpenFOAM中邊界條件的設(shè)置Table 2 Settings of boundary conditions in OpenFOAM

模型在入口邊界和出口邊界使用主動(dòng)消波。主動(dòng)消波是物模與數(shù)模中其中一個(gè)重要的消波方法。在寬闊的海域里，波浪反射的地方距離研究區(qū)域甚遠(yuǎn)。但在物模與數(shù)模中受到空間的限制，如有限長(zhǎng)度的波浪水槽和受計(jì)算機(jī)性能限制的有限網(wǎng)格數(shù)的計(jì)算單元體反射的波浪會(huì)影響研究區(qū)域的實(shí)驗(yàn)結(jié)果。以前普遍采用的是在波浪出口邊界設(shè)置海綿層、多孔介質(zhì)等被動(dòng)消波區(qū)域，或者是人為增加一個(gè)消波阻力。被動(dòng)消波存在缺點(diǎn)，如對(duì)于長(zhǎng)波無法在消波層上消除，因?yàn)樗黄扑椋虼巳匀粫?huì)發(fā)生反射。在數(shù)模中，學(xué)者多使用基于雷諾應(yīng)力平均的納維埃-斯托克斯方程的多孔介質(zhì)和松弛區(qū)域的消波方法［22-23］。另一種方法是主動(dòng)消波［24］，這種方法最初是在推波板造波的物模實(shí)驗(yàn)中實(shí)現(xiàn)，該造波機(jī)通過測(cè)量波高來調(diào)整移動(dòng)速度，以此產(chǎn)生需要的波浪，同時(shí)防止反射至造波處的波浪二次反射影響造波。這種主動(dòng)消波的方法也適用于只用來消波的出口邊界。在數(shù)模中，由于波浪邊界條件固定，所以這種方法通過加入校正速度邊界的方法來實(shí)現(xiàn)。在二維波浪數(shù)模中，為了消除反射波，邊界必須產(chǎn)生一個(gè)與入射波的速度相同，但方向相反、指向計(jì)算區(qū)域的波浪，這也是主動(dòng)消波法的實(shí)現(xiàn)原理。本次數(shù)模在入口、出口處設(shè)置主動(dòng)消波。

2 數(shù)模設(shè)置及數(shù)值結(jié)果分析

本文采用Yao等［19］的實(shí)驗(yàn)進(jìn)行數(shù)模的驗(yàn)證。實(shí)驗(yàn)在長(zhǎng)36m、寬0.55m、高0.6m的波浪水槽中進(jìn)行，其中水深0.45m，在水槽一端設(shè)置平移式推波板造波，在距離推波板16.35m處設(shè)置1∶6的斜坡，斜坡由PVC材料制作，斜坡長(zhǎng)2.1m，斜坡后接水平平臺(tái)，水平平臺(tái)長(zhǎng)9.9 m、距離水槽底部0.35m。在波浪水槽的兩端均設(shè)置多孔介質(zhì)材料以減少波浪反射的影響，具體設(shè)置詳見文獻(xiàn)［19］。

本文的數(shù)模計(jì)算區(qū)域示意圖如圖1所示。二維數(shù)值波浪水槽長(zhǎng)25m、高0.75m。選定與浪高儀相同位置的12個(gè)測(cè)點(diǎn)的水體體積分?jǐn)?shù)計(jì)算波高，計(jì)算區(qū)域的坐標(biāo)原點(diǎn)選在造波邊界入口與底部邊界的交點(diǎn)處。12個(gè)測(cè)點(diǎn)的坐標(biāo)值如表3所示。Open-FOAM有兩種網(wǎng)格劃分工具，分別為blockMesh和snappyHexMesh。blockMesh用于創(chuàng)建由8個(gè)頂點(diǎn)、6個(gè)邊界面組成的“塊”?！皦K”的生成需要定義8個(gè)頂點(diǎn)的坐標(biāo)位置和邊界面的邊界條件，這里使用blockMesh生成結(jié)構(gòu)矩形背景網(wǎng)格。snappyHexMesh使用三角形構(gòu)成的STL文件切割blockMesh生成的六面體。這里使用snappyHexMesh生成斜坡和礁坪面，并移除斜坡以下的非計(jì)算區(qū)域的網(wǎng)格。整個(gè)計(jì)算區(qū)域橫向網(wǎng)格單元長(zhǎng)度為0.02m，縱向網(wǎng)格長(zhǎng)度為0.01m，在水面處附近進(jìn)行橫向加密，加密后的橫向單元長(zhǎng)度為0.01m，縱向單元長(zhǎng)度為0.005m。斜坡位置上的網(wǎng)格示意圖如圖1所示。坐標(biāo)原點(diǎn)選取在測(cè)點(diǎn)的位置見圖1。

圖1 數(shù)值模型計(jì)算區(qū)域和浪高儀位置Fig.1 Computational domains of numerical model and the position of wave gauges

表3 浪高儀的坐標(biāo)位置Table 3 Coordinates of wave gauges

本文研究的規(guī)則波的入射波高為0.095 m，周期為1.25s，為斯托克斯二階波［25］，計(jì)算得Iribarren數(shù)為2.36，因此屬于卷破波。模型計(jì)算耗時(shí)60s。

2.1 湍流動(dòng)能與雷諾應(yīng)力

在雷諾應(yīng)力模型中，湍流動(dòng)能表示為

式中，U′為脈動(dòng)流速。在k-ω SST模型中，湍流動(dòng)力黏度表示為

湍流動(dòng)能的源項(xiàng)表示為

對(duì)于波陡較大的波浪 (波高與波長(zhǎng)的比值H/L＞0.07)，在水面處容易產(chǎn)生非常大的湍流動(dòng)能，這是因?yàn)榭諝馀c水有很大的密度比 (ρa(bǔ)/ρw=1∶1000)。

2.1.1 深水區(qū)湍流動(dòng)能

圖2分別是標(biāo)準(zhǔn)k-ω SST模型和兩相k-ω SST模型下4個(gè)波長(zhǎng)范圍內(nèi)的湍流動(dòng)能分布云圖，圖中的左端是入口邊界，黑色線表示是水面。標(biāo)準(zhǔn)k-ω SST模型在波浪破碎前段的水面處存在著湍流動(dòng)能過度偏大的現(xiàn)象，這會(huì)引起波浪在傳播過程中波能過度消耗，使得傳播到破碎帶前的波高預(yù)估不足;兩相k-ω SST模型在水面處限制了湍流動(dòng)能發(fā)生不合理的增長(zhǎng)，因而也降低了湍流動(dòng)力黏度，避免在使用湍流模型后波高發(fā)生衰減。

2.1.2 破波帶湍流動(dòng)能

圖3分別是標(biāo)準(zhǔn)k-ω SST模型和兩相k-ω SST模型下斜坡上的湍流分布云圖，其中t/T為無量綱化時(shí)間。以41.5s開始為例，兩種模型在破碎前的自由表面湍流動(dòng)能分布有很大的差異。在波浪破碎位置上，兩相k-ω SST模型的波浪破碎前的湍流動(dòng)能很小，破碎后的湍流動(dòng)能從破波點(diǎn)開始急劇增大，在波峰前端位置形成具有較大湍流動(dòng)能的區(qū)域。破碎后的湍流動(dòng)能繼續(xù)向礁岸傳播，并在礁坪上的傳播過程中衰減。此外，在波浪破碎前，兩種模型在右端礁坪區(qū)域的湍流動(dòng)能都比較大，這表明在下一個(gè)波浪發(fā)生破碎時(shí)，前一個(gè)波浪產(chǎn)生的湍流動(dòng)能還未完全耗散。

圖2 深水區(qū)湍流動(dòng)能云圖Fig.2 Turbulence energy contours in the deep water area

2.1.3 雷諾應(yīng)力

其中I為單位矩陣。因此，在雷諾應(yīng)力模型中，計(jì)算雷諾應(yīng)力改為計(jì)算湍流動(dòng)能和湍流運(yùn)動(dòng)黏度，并由式 (12)確定雷諾應(yīng)力。雷諾應(yīng)力的大小=。圖4為雷諾應(yīng)力的分布云圖，從圖中可以看出，波浪在發(fā)生卷破后水面附近的雷諾應(yīng)力增大，雷諾應(yīng)力的空間分布類似于湍流動(dòng)能，主要集中在破碎位置的波峰處。

2.2 波高與平均水位

圖5為波高和平均水位的數(shù)模模擬結(jié)果和實(shí)驗(yàn)值的對(duì)比。圖中的RMSE代表均方根誤差，可按下式計(jì)算:

式中:N為測(cè)量點(diǎn)總數(shù);Qic為第i個(gè)測(cè)量點(diǎn)的變量模擬值;Qi

圖3 標(biāo)準(zhǔn)k-ω SST模型和兩相k-ω SST模型下斜坡上的湍流動(dòng)能云圖 (僅水體)Fig.3 Turbulence energy contours on the slope with standard kω SST and multiphase k-ω SST models(Only for water phase)

m為第i個(gè)測(cè)量點(diǎn)的變量實(shí)驗(yàn)值。圖5(a)給出了波浪在標(biāo)準(zhǔn)k-ω SST模型和兩相k-ω SST模型下的數(shù)模結(jié)果與文獻(xiàn) ［19］的實(shí)驗(yàn)結(jié)果沿珊瑚礁方向傳播的波高分布。其中，RMSE文獻(xiàn)、RMSE標(biāo)準(zhǔn)、RMSE兩相分別為由文獻(xiàn) ［19］、標(biāo)準(zhǔn)k-ω SST模型、兩相k-ω模型獲得的均方根誤差。從圖中可以看出，兩個(gè)模型在破碎前段波高都與實(shí)驗(yàn)值相差不大，波高的沿程起伏變化是由于礁前陡坡對(duì)波浪有反射作用，實(shí)際波面是由入射波和反射波疊加而成，因而在深水區(qū)形成了不完全駐波。值得注意的是，相較于兩相k-ω SST模型，標(biāo)準(zhǔn)k-ω SST模型在礁前發(fā)生了明顯的衰減，這是考慮了湍流模型后，湍流動(dòng)能在水面處產(chǎn)生了不合理的源項(xiàng)引起的。當(dāng)波浪到達(dá)斜坡頂端附近時(shí)，波高達(dá)到最大，波浪發(fā)生破碎，波能急劇耗散，波高急劇變小。由于底部是光滑面，底摩阻引起的波能損失較小，而在礁坪上，波浪破碎后的波高由當(dāng)?shù)厮顩Q定。Yao等［19］的數(shù)模中模擬湍流采用的是混合長(zhǎng)度模型，從圖5(a)中可以看出，該模型在破碎前的波高預(yù)測(cè)值較實(shí)驗(yàn)值和k-ω SST模型都偏大，這是由于混合長(zhǎng)度模型低估了破碎前的湍流強(qiáng)度，使得湍流運(yùn)動(dòng)黏度νt偏小［26］。圖5(b)為波浪在標(biāo)準(zhǔn)k-ω SST模型和兩相k-ω SST模型下的數(shù)模結(jié)果與實(shí)驗(yàn)結(jié)果沿珊瑚礁方向傳播的平均水位分布。平均水位通過對(duì)自由表面的歷時(shí)線取時(shí)間平均求得。以靜水位為縱軸零點(diǎn)，正值表示增水，反之為減水。從圖中可以看出，兩種模型在深水區(qū)的模擬結(jié)果相差不大，在礁前都發(fā)生了減水，平均水位在礁坪邊緣處達(dá)到了最低值，兩相k-ω SST模型的礁前減水更明顯一些。數(shù)值模型的水位回升比實(shí)驗(yàn)值緩慢一些，這也是由波高衰減得比較緩慢引起的。而在礁坪區(qū)域，兩相k-ω SST模型的增水比標(biāo)準(zhǔn)k-ω SST模型要更大一些。

圖4 標(biāo)準(zhǔn)k-ω SST模型和兩相k-ω SST模型下斜坡上的雷諾應(yīng)力云圖 (僅水體)Fig.4 Reynolds stress contours on the slope with standard k-ω SST and multiphase k-ω SST models(Only for water phase)

圖5 數(shù)值模型與實(shí)驗(yàn)值的波高和平均水位對(duì)比Fig.5 Comparisons of wave height and mean water level between numerical and experimental results

2.3 波浪傳播變形

圖6 比較了兩種數(shù)學(xué)模型下，在具有代表性的浪高儀位置 (#1、#3、#5、#7、#9、#11)處的自由表面的時(shí)間序列。從整體上看，兩種模型在連續(xù)6個(gè)周期內(nèi)的波形基本保持一致，這說明k-ω SST模型在模擬波浪破碎的過程時(shí)，不僅在破碎前有很好的穩(wěn)定性，在破碎后的區(qū)域也能夠保持波形穩(wěn)定。浪高儀#1位于破碎前段，表現(xiàn)為波峰和波谷對(duì)稱的規(guī)則波，兩種數(shù)學(xué)模型與實(shí)驗(yàn)值基本相符。當(dāng)波浪傳播到浪高儀#3位置時(shí)，底部為礁前斜坡，淺化作用使得單個(gè)周期內(nèi)的波峰和波谷開始出現(xiàn)不對(duì)稱特征，具體表現(xiàn)為波峰向前傾斜、波谷變長(zhǎng)。值得注意的是，標(biāo)準(zhǔn)k-ω SST模型的波高明顯偏低，而且出現(xiàn)了相位側(cè)移。浪高儀#5位于破碎區(qū)內(nèi)，波高急劇變小，波峰和波谷嚴(yán)重不對(duì)稱，實(shí)驗(yàn)結(jié)果出現(xiàn)明顯的鋸齒狀，周期間的波形差異較大，而兩個(gè)數(shù)值模型在破波區(qū)內(nèi)的波能耗散比較緩慢，與實(shí)驗(yàn)值相比都偏大，這表明數(shù)值模型在破碎帶的能量耗散得比實(shí)驗(yàn)值慢。此外，本文采用體積平均法計(jì)算波高，這在破碎帶外計(jì)算波高確實(shí)可行，但在破碎區(qū)，由于流體體積法 (VOF)計(jì)算出來的波面附近含有強(qiáng)烈的破碎和摻氣作用，水面處的一些單元格水體發(fā)生飛濺，空氣也卷入到了水體中，這一區(qū)域的波高計(jì)算方法仍然有待探討。在浪高儀#7、#9處，波浪破碎停止，在礁坪區(qū)域重新生成了波高較小的段波。由于礁坪區(qū)域水深變淺，受非線性作用影響而呈現(xiàn)顯著的非線性特征，具體表現(xiàn)為這兩處的波峰變短抬升且向前傾，波谷變長(zhǎng)且更加平坦，不對(duì)稱特征非常明顯。

圖6 浪高儀#1、#3、#5、#7、#9、#11處波面高度的時(shí)間序列Fig.6 Time-series of the surface elevations at six wave gauges#1，#3，#5，#7，#9，#11

2.4 流速分布

圖7 是兩種湍流模型在浪高儀#7處一個(gè)周期內(nèi)不同時(shí)刻的流速垂向分布圖。以波浪開始發(fā)生破碎的時(shí)刻為零點(diǎn)，時(shí)間間隔為0.1 s。可以看出，在浪高儀#7處底部存在著明顯的下層逆流，且兩相k-ω SST模型的逆向流速要更大一些。此外，水平流速梯度在破碎前后的方向相反，底部水體的流速很小，并沿著水面的方向流速逐漸增大。另外，k-ω SST模型下的最大流速略大于標(biāo)準(zhǔn)k-ω SST模型，位置也更接近礁坪面。

圖7 標(biāo)準(zhǔn)k-ω SST模型與兩相k-ω SST的水平速度垂向分布比較Fig.7 Comparisons of vertical distributions of horizontal velocity between standard k-ω SST model and multiphase k-ω SST model

3 結(jié)論

本文使用OpenFOAM中的標(biāo)準(zhǔn)k-ω SST模型和兩相k-ω SST模型模擬了波浪在珊瑚礁地形上的破碎過程，并對(duì)數(shù)值模型的波高、水位、流速進(jìn)行了分析?？傮w上看，數(shù)模結(jié)果與實(shí)驗(yàn)結(jié)果大體一致，主要結(jié)論如下:

(1)兩相k-ω SST湍流模型相比標(biāo)準(zhǔn)k-ω SST模型對(duì)波高、水位的模擬更加準(zhǔn)確。密度修正后的模型在波浪破碎前的層流區(qū)域能有效減輕湍流動(dòng)能的過度產(chǎn)生，使得波浪在傳播過程中不會(huì)發(fā)生不合理的明顯衰減，也使得對(duì)破碎前的波高預(yù)測(cè)更加準(zhǔn)確。此外，兩相k-ω SST湍流模型在礁坪上的增水更加明顯，更接近實(shí)驗(yàn)值。

(2)在破碎區(qū)域，在自由表面附近產(chǎn)生了明顯的湍流動(dòng)能且主要集中在波峰處，這說明在波浪破碎區(qū)域應(yīng)當(dāng)使用湍流模型，且選擇適當(dāng)?shù)耐牧髂Ｐ蛯?duì)模擬結(jié)果精度十分重要。

(3)使用兩相流的模擬結(jié)果表明，在波浪的破碎過程中，不管是上層空氣還是底部水體，都存在著明顯的逆向梯度流，使用k-ω SST模型計(jì)算得到的波面具有較好的穩(wěn)定性。