童朝鋒，魏芷陽(yáng)，孟艷秋

(河海大學(xué) 港口海岸與近海工程學(xué)院，江蘇 南京 210098)

基于FLUENT流體力學(xué)模型建立數(shù)值波浪水槽是波浪與建筑物之間相互作用研究的經(jīng)典方法，周勤俊等[1]基于RANS方程和兩方程湍流模型，采用有限體積法，首次提出了適用于VOF方法的源造波、消波技術(shù);其后孫大鵬等[2]、趙艷等[3]采用計(jì)算流體力學(xué)模型成功利用源項(xiàng)、推板、速度邊界造波法建立基于FLUENT的數(shù)值波浪水槽，F(xiàn)LUENT流體力學(xué)模型開(kāi)始廣泛應(yīng)用。近年來(lái)國(guó)內(nèi)外新的數(shù)值造波方法層出不窮，Jacobsen等[4]將基于BOUSSINESQ水波方程的造波方法成功應(yīng)用到了開(kāi)源代碼OpenFOAM中取得了較好的造波效果;Wang等[5]基于CIP模型成功模擬了孤立波與水下建筑物的相互作用；Losada 等[6]、Higuera等[7]分別建立了基于COBRAS-UC模型、OpenFOAM模型的數(shù)值波浪水槽模型模擬了波浪與建筑物之間相互作用過(guò)程，取得了較好的效果。

相比之下，F(xiàn)LUENT作為最常用的計(jì)算流體力學(xué)模型，操作簡(jiǎn)便、界面友好，具有先進(jìn)的數(shù)值仿真和強(qiáng)大的前、后處理功能，在計(jì)算流體力學(xué)領(lǐng)域和實(shí)際工程中具有極強(qiáng)的應(yīng)用價(jià)值?；贔LUENT的數(shù)值造波方法根據(jù)造波原理主要有源項(xiàng)造波法、推板造波法和速度邊界造波法3種。為達(dá)到波浪的精準(zhǔn)數(shù)值仿真，須充分了解各種造波法造波效果以及達(dá)到最佳造波效果應(yīng)采用的網(wǎng)格尺寸，但目前尚無(wú)相關(guān)文獻(xiàn)對(duì)此進(jìn)行詳細(xì)對(duì)比分析。

為此，本文應(yīng)用FLUENT流體力學(xué)計(jì)算模型，采用上述3種造波法，分別建立不同網(wǎng)格尺寸和波陡條件下的數(shù)值波浪水槽，結(jié)合線性波和二階、三階Stokes波理論，對(duì)比不同方法的仿真結(jié)果及與理論波浪形態(tài)之差異，分析3種造波法的優(yōu)缺點(diǎn)，討論采用不同波浪方程、不同網(wǎng)格尺寸對(duì)仿真結(jié)果的影響，為利用數(shù)值波浪水槽開(kāi)展波浪與建筑物相互作用數(shù)值研究提供根據(jù)。

1 造波、消波原理及邊界條件設(shè)置

1.1 數(shù)值波浪仿真模型

波浪數(shù)值仿真采用FLUENT流體力學(xué)軟件，其軟件選取顯式或隱式差分格式，非耦合隱式算法、耦合顯式算法或耦合隱式算法等數(shù)值計(jì)算方法求解Navier-Stokes方程，可精確模擬層流、湍流、多相流等各類流體力學(xué)問(wèn)題，軟件提供用戶自定義函數(shù)便于二次開(kāi)發(fā)。本文數(shù)值波浪水槽模擬采用二階迎風(fēng)格式、PISO算法、k-ε模型，結(jié)合VOF模型和幾何重構(gòu)法進(jìn)行液面捕捉，獲得較為精準(zhǔn)的模擬結(jié)果。

1.2 造波原理及邊界條件

1.2.1源項(xiàng)造波法

源項(xiàng)造波法原理是通過(guò)UDF二次編譯功能在造波源區(qū)域給動(dòng)量方程添加動(dòng)量源項(xiàng)s(x,y,t)，使造波源區(qū)域中水質(zhì)點(diǎn)具有運(yùn)動(dòng)速度，并在重力作用下水面升降形成波浪。根據(jù)相關(guān)文獻(xiàn)[8-9]，造波源區(qū)域高度宜取水深的23，或使造波源區(qū)域上方水體高度約為1倍波高，寬度宜取波長(zhǎng)的150。實(shí)際計(jì)算時(shí)可選取大致合適的尺寸，再不斷試算調(diào)整選取最佳造波源尺寸。

數(shù)值水浪水槽如圖1所示。造波源區(qū)域設(shè)在入口處，波浪沿x方向傳播，右端為消波區(qū)，消波區(qū)長(zhǎng)度取1～2倍波長(zhǎng)。水槽入口、出口及下部邊界均取固壁邊界，上部為壓力出口。

圖1 數(shù)值波浪水槽

造波源函數(shù)根據(jù)以下公式確定：

(1)

(2)

式中：s(x,y,t)為造波源區(qū)域要添加的造波源項(xiàng)；u(x,y,t)為t時(shí)刻目標(biāo)波浪在坐標(biāo)(x,y)處的水質(zhì)點(diǎn)水平速度，根據(jù)線性波理論或Stokes波理論取值；u、v分別為x、y方向速度；Δx為水平方向網(wǎng)格長(zhǎng)度。

1.2.2推板造波法

推板造波法原理與物理模型中推板式造波機(jī)相同，其數(shù)值波浪水槽與圖1相同，入口處固壁邊界添加動(dòng)網(wǎng)格以模擬推波板運(yùn)動(dòng)。

推板運(yùn)動(dòng)速度根據(jù)Schaffer等[10-11]推得的推板位移變化與波面高度之間解析關(guān)系確定，其水平運(yùn)動(dòng)速度函數(shù)u(t)如下：

(3)

(4)

式中：η為波面升高，根據(jù)Stokes波理論或線性波理論確定其值；ω為角頻率；c0為造波板運(yùn)動(dòng)與給定波浪間的傳遞函數(shù)；k為波數(shù)；d為水深。

1.2.3速度邊界造波法

速度邊界造波法將入口邊界水質(zhì)點(diǎn)賦予已知相應(yīng)波水平和垂向周期性振蕩速度，促使水質(zhì)點(diǎn)運(yùn)動(dòng)產(chǎn)生波浪，其數(shù)值水槽與圖1相同。數(shù)值水槽的造波區(qū)入口處邊界設(shè)為速度入口邊界，水槽底部和右側(cè)為固壁邊界，上部為壓力出口邊界。

數(shù)值水槽入口邊界處水質(zhì)點(diǎn)水平向和垂向速度：

u=u(x,y,t)

(5)

v=v(x,y,t)

(6)

式中：u，v分別為入口邊界水平、豎直方向速度；u(x,y,t)、v(x,y,t)分別為在t時(shí)刻坐標(biāo)(x,y)處水平、豎向速度，其值根據(jù)線性波理論或Stokes波理論確定。

1.3 消波方法

采用動(dòng)量源消波。在消波區(qū)，水體動(dòng)量方程添加阻尼源項(xiàng)：

(7)

(8)

式中：ρ為水密度；p為水壓力；μ為消波系數(shù)，為使波浪平緩過(guò)渡，μ隨x由0遞增，即:

μ=a(x-x0)

(9)

式中：a為常數(shù)，根據(jù)經(jīng)驗(yàn)其取值范圍為1～50；x0為消波區(qū)起始位置水平坐標(biāo)。

2 模型建立

2.1 網(wǎng)格及模擬工況

波浪要素見(jiàn)表1，波況1～3保持水深和波高不變，只改變周期達(dá)到改變波陡的目的。波況1～3采用的水槽總長(zhǎng)280 m、高10 m、水深6 m，消波區(qū)長(zhǎng)30 m。為方便比較，考慮盡可能少的影響條件，波況4采用總長(zhǎng)140 m、高8 m、水深6 m的深水?dāng)?shù)值水槽，探究非線性項(xiàng)對(duì)造波效果的影響，水槽消波區(qū)長(zhǎng)20 m。

表1 數(shù)值水槽波況參數(shù)

采用a、b、 c、d、e共5種結(jié)構(gòu)化網(wǎng)格，其中網(wǎng)格a水平分辨率Δx=0.4 m，垂向分辨率Δy=0.02 m；網(wǎng)格b的Δx=0.4 m，Δy從液面附近向兩邊遞增，最小為0.02 m，最大為0.1 m；網(wǎng)格c的Δx=0.4 m、Δy=0.05 m；網(wǎng)格d的Δx=0.4 m、Δy=0.12 m；網(wǎng)格e的Δx=0.1 m，Δy從液面附近向兩邊遞增，最小為0.004 m，最大為0.06 m。

為探究垂向網(wǎng)格分辨率對(duì)模擬結(jié)果的影響，波況1采用a、b、c、d4種網(wǎng)格分別建立數(shù)值波浪水槽。網(wǎng)格e用于波況4深水情況下微幅波的模擬計(jì)算。波況1～3分別采用3種造波法和線性波、二階Stokes波兩種波浪方程進(jìn)行模擬，以探究造波法對(duì)造波效果的影響。波況4采用3種造波法和線性波、二階、三階Stokes波方程進(jìn)行模擬，研究數(shù)值造波中非線性作用的影響。算例模擬時(shí)間步長(zhǎng)取0.01 s，每步最多迭代40次，計(jì)算總時(shí)間超過(guò)80 s。

2.2 造波和消波驗(yàn)證

使用波況1(波陡Hλ=0.043)、網(wǎng)格b、二階波浪方程，則3種造波方法得出的x=λ處水面變化歷時(shí)曲線如圖2所示。3種模擬結(jié)果都在時(shí)間t=10 s左右達(dá)到穩(wěn)定，之后水面隨時(shí)間規(guī)律性變化，且與理論值吻合良好。模擬得出的x=λ處一個(gè)周期內(nèi)速度和動(dòng)壓強(qiáng)與理論值對(duì)比見(jiàn)圖3。取波峰處時(shí)間為T(mén)4,波谷處時(shí)間為3T4，可以看出模擬值與理論值基本吻合，說(shuō)明基于FLUENT建立的數(shù)值波浪水槽可用于波浪與結(jié)構(gòu)物之間相互作用的研究。

圖2 3種造波法x=λ處液面歷時(shí)曲線

波況1使用推板造波法、二階波浪方程,t=60 s時(shí)水槽末端水質(zhì)點(diǎn)速度矢量分布見(jiàn)圖4。結(jié)果顯示在非消波區(qū)水質(zhì)點(diǎn)運(yùn)動(dòng)速度較大，消波區(qū)水質(zhì)點(diǎn)運(yùn)動(dòng)速度逐漸減小，最后變?yōu)?，說(shuō)明消波區(qū)消波效果良好。

圖3 推板造波法x=λ處水質(zhì)點(diǎn)速度、動(dòng)壓強(qiáng)與理論值對(duì)比

圖4 t=60 s時(shí)水質(zhì)點(diǎn)速度矢量分布

3 結(jié)果對(duì)比與分析

3.1 網(wǎng)格分辨率對(duì)數(shù)值造波效果影響

采用上述3種造波法、二階波浪方程和4種不同垂向分辨率的網(wǎng)格分別建立數(shù)值波浪水槽，模擬仿真波況1的波浪。為便于辨析某一位置處模擬波面與理論波面貼合程度，將自由水面和波高的理論值均取該點(diǎn)處模擬波的平均水面和波高，采用如下均方根相對(duì)誤差計(jì)算一個(gè)周期內(nèi)模擬值與理論值相對(duì)差異，評(píng)價(jià)模擬仿真效果：

(10)

式中：S為均方根相對(duì)誤差；H為波高；ηi為i時(shí)刻理論液面高度；η′i為i時(shí)刻模擬液面高度；n為一個(gè)周期內(nèi)的數(shù)據(jù)數(shù)量。

不同造波法所造波浪平均水面飄移隨網(wǎng)格分辨率變化如圖5a)所示。源項(xiàng)造波法模擬仿真的水面漂移絕對(duì)值受網(wǎng)格尺寸變化影響較大，液面附近垂向網(wǎng)格分辨率Δy取0.02 m時(shí)波面漂移絕對(duì)值為0 m，隨Δy增加，迅速增加到0.07 m左右；推板造波法模擬結(jié)果在Δy≤0.05 m時(shí)水面漂移絕對(duì)值集中在0.015 m附近，當(dāng)Δy=0.12 m時(shí)偏移量急劇增加至0.06 m；速度邊界造波法模擬仿真結(jié)果隨網(wǎng)格步長(zhǎng)增加波面漂移量變化不明顯，基本保持在0.02 m左右。

液面高度相對(duì)差異隨網(wǎng)格分辨率變化如圖5b)所示。采用源項(xiàng)造波法得出的模擬波形與理論波形之間誤差最大，隨網(wǎng)格變化誤差增加最快；采用推板造波法得出的模擬波形與理論波形吻合最好，相對(duì)誤差均較小；使用速度邊界法模擬值在Δy≤0.05 m時(shí)誤差較小，Δy=0.12 m時(shí)誤差急劇增加。

注：波況1，Hλ=0.043,x=λ處；Δy為0.02、0.02～0.10、0.05、0.12 m分別對(duì)應(yīng)網(wǎng)格a、b、c、d。

圖5采用不同造波方法的水面漂移絕對(duì)值和仿真與理論液面高度相對(duì)差異隨不同網(wǎng)格分辨率變化

波高沿程變化如圖6所示?？梢钥闯觯捎貌煌直媛实木W(wǎng)格建立數(shù)值波浪水槽，隨波浪傳播距離增加，采用不同網(wǎng)格模擬得出的波高值產(chǎn)生差異。其中使用源項(xiàng)造波法造波得出的波高差異最大，使用推板造波法得出的波高差異最小，采用速度邊界造波法得出的沿程波高在網(wǎng)格垂向分辨率Δy≤0.05 m時(shí)基本相等，使用Δy=0.12 m的網(wǎng)格時(shí)得出的沿程波高與使用其他網(wǎng)格得出的模擬值差異較大。波高隨傳播距離增加持續(xù)衰減，衰減到一定程度后逐漸穩(wěn)定，這是由數(shù)值耗散和水體黏性影響產(chǎn)生的正?，F(xiàn)象，實(shí)際應(yīng)用中可以增加起始波高，在波高趨于穩(wěn)定的位置放置建筑物。

圖6 波況1使用不同造波法沿程波高對(duì)比

綜上，垂向網(wǎng)格步長(zhǎng)越小模擬結(jié)果越精確，源項(xiàng)造波法模擬仿真波浪效果對(duì)垂向網(wǎng)格分辨率變化比較敏感，液面附近垂向網(wǎng)格分辨率取Δy=0.02 m，即Δy=H50(H為波高)時(shí)即可取得較為精確的結(jié)果；推板造波法和速度邊界造波法模擬結(jié)果對(duì)網(wǎng)格分辨率變化敏感性相對(duì)較弱，垂向網(wǎng)格分辨率Δy=0.05 m，即Δy=H20時(shí)模擬效果較好。

3.2 波陡對(duì)數(shù)值造波效果影響

根據(jù)上述網(wǎng)格對(duì)3種造波法造波效果影響分析結(jié)果，采用網(wǎng)格b對(duì)波況1～3模擬計(jì)算，討論波陡對(duì)數(shù)值造波效果影響。計(jì)算分別采用了3種造波法和2種波浪方程。

模擬結(jié)果在x=25 m處水面漂移絕對(duì)值見(jiàn)圖7?？梢钥闯觯ǘ赶嗤瑫r(shí)采用速度邊界法得出的水面漂移絕對(duì)值最大，源項(xiàng)法最?。痪筒ɡ朔匠潭?，其他條件相同時(shí)，采用不同波浪方程得出的水面漂移絕對(duì)值基本相等，說(shuō)明采用的波浪方程不是影響水面漂移的主要因素。

圖7 不同波陡條件下采用不同波浪方程模擬x=25 m處水面漂移絕對(duì)值

模擬結(jié)果在x=25 m處一個(gè)周期內(nèi)液面高度與理論值的相對(duì)誤差S見(jiàn)圖8?？梢钥闯?，其他條件相同時(shí)，采用不同波浪方程仿真得出的結(jié)果與二階Stokes波理論值吻合均較好。說(shuō)明造波源不管采用哪種波浪方程，在波浪傳播過(guò)程中，模擬結(jié)果都會(huì)受邊界的非線性作用影響，使波浪呈現(xiàn)出非線性波特性。

圖8 不同波陡和波浪方程模擬液面與理論液面誤差

仿真得出的沿程平均波高對(duì)比見(jiàn)圖9?？梢钥闯觯捎盟俣冗吔缭觳ǚǖ贸龅难爻唐骄ǜ咦畲?，采用源項(xiàng)法最??；其他條件相同，波陡增大波高衰減速度變快；只改變波浪方程對(duì)比模擬結(jié)果，使用二階波浪方程得出的沿程平均波高普遍高于使用一階方程得出的結(jié)果，波陡越大兩者差異越不明顯。

圖9 采用不同造波法和波浪方程模擬沿程平均波高對(duì)比

3.3 數(shù)值造波中的非線性作用影響

根據(jù)上述波陡對(duì)數(shù)值造波效果影響結(jié)果，無(wú)論使用哪種波浪方程建立數(shù)值波浪水槽，在波浪傳播過(guò)程中都會(huì)使模擬波面變化接近二階Stokes波理論值。增加波況4(Hλ=0.021)，深水?dāng)?shù)值水槽算例，研究采用不同波浪方程對(duì)造波效果影響。計(jì)算分別采用一階、二階、三階波浪方程,其水面漂移絕對(duì)值見(jiàn)表2。可以看出，造波方法相同時(shí)，使用不同波浪方程得出的水面漂移絕對(duì)值完全相等。模擬波面與理論波面相對(duì)誤差如圖10所示，造波方法相同時(shí)采用的波浪方程階數(shù)越高模擬值與同階理論值之間相對(duì)誤差越小；低階波傳播過(guò)程中受非線性影響，模擬值與相應(yīng)理論值相對(duì)誤差反而大。高階波波形已經(jīng)體現(xiàn)了非線性項(xiàng)的變形作用，傳播過(guò)程中反而變形較小，其模擬值與理論值之間誤差也較小。

表2 波況4在深水情況下水面漂移絕對(duì)值

圖10 波況4在深水條件下模擬波面與理論波面相對(duì)誤差

4 結(jié)論

1)采用源項(xiàng)造波法造波，波高沿水槽衰減最快；推板造波法采用動(dòng)網(wǎng)格技術(shù)，計(jì)算效率相對(duì)低，數(shù)值造波波形與理論波形最接近；采用速度邊界造波法計(jì)算效率最高，數(shù)值波高沿水槽衰減速度最小。

2)網(wǎng)格越細(xì)模擬結(jié)果越精確，采用源項(xiàng)造波法模擬結(jié)果對(duì)網(wǎng)格尺寸變化最敏感，采用推板造波法和速度邊界造波法模擬結(jié)果對(duì)網(wǎng)格尺寸變化較不敏感。

3)數(shù)值水槽計(jì)算中，波陡越大波高衰減越快；數(shù)值波浪模擬顯示，在數(shù)值水槽的波浪傳播有非線性作用影響，導(dǎo)致即使造波是線性波，傳播到一定距離后，波面也有非線性變形。