劉紅兵 陳國(guó)明 呂 濤 黃 翱

(1. 中國(guó)石油大學(xué)(華東)海洋油氣裝備與安全技術(shù)研究中心 山東青島 266580； 2. 哈爾濱工程大學(xué)船舶工程學(xué)院 黑龍江哈爾濱 150001)

在復(fù)雜的海洋環(huán)境中，風(fēng)、浪、流通常是共同存在且相互交替運(yùn)動(dòng)，從而產(chǎn)生類似于“Doppler”效應(yīng)的風(fēng)-浪-流耦合作用[1]，導(dǎo)致流場(chǎng)結(jié)構(gòu)和波浪幅值、波高、周期等波浪要素發(fā)生改變，使得波浪發(fā)生折射乃至破碎[2]。因此，準(zhǔn)確模擬海洋風(fēng)-浪-流時(shí)空分布特征，探討風(fēng)-浪-流耦合作用機(jī)理以及耦合作用下風(fēng)-浪-流場(chǎng)分布特征，對(duì)于保證海洋工程結(jié)構(gòu)物安全運(yùn)行具有重要意義[3]。

目前，國(guó)內(nèi)外關(guān)于海洋環(huán)境水動(dòng)力性能的研究主要采用理論研究、實(shí)尺試驗(yàn)、模型試驗(yàn)和數(shù)值水池模擬等方法?；跓o(wú)黏無(wú)旋假設(shè)的勢(shì)流理論可以快速計(jì)算波浪傳播及其對(duì)海洋工程結(jié)構(gòu)物的作用，卻無(wú)法計(jì)算由于黏性效應(yīng)引起的一些物理現(xiàn)象，如風(fēng)-浪-流耦合遠(yuǎn)程傳播和漩渦演變等，具有較大的局限性[4]。實(shí)尺試驗(yàn)法受海洋環(huán)境影響因素較大，且試驗(yàn)成本高并難以控制，在海洋工程應(yīng)用中難以實(shí)現(xiàn)[5]。模型試驗(yàn)法主要采用一定的縮尺比例轉(zhuǎn)化關(guān)系，通過(guò)縮尺模型水動(dòng)力學(xué)性能反映整體結(jié)構(gòu)水動(dòng)力學(xué)性能，在海洋工程中應(yīng)用較為廣泛，但存在試驗(yàn)成本高、試驗(yàn)周期長(zhǎng)、試驗(yàn)縮尺模型和原型之間準(zhǔn)確轉(zhuǎn)換困難等缺陷，使其應(yīng)用仍受到一定限制[6]。隨著近年計(jì)算機(jī)技術(shù)高速發(fā)展，數(shù)值波浪水池以其低成本、短周期和高精度等特點(diǎn)，在海洋工程波浪誘導(dǎo)運(yùn)動(dòng)和波浪載荷預(yù)測(cè)等方面得到了廣泛應(yīng)用[7]。例如，Wang[8]和蘇曉杰 等[9]基于歐拉-拉格朗日算法追蹤流體瞬時(shí)自由液面的變化過(guò)程，建立了完全非線性的三維數(shù)值波浪水池。Alamian等[10]利用活塞運(yùn)動(dòng)造波法，在數(shù)值波浪水池中模擬了線性波和非線性波；Finnegan等[11]基于CFX軟件，建立了搖板式數(shù)值波浪水池；梁修峰 等[12]基于FLUENT軟件建立了一種面向海洋工程應(yīng)用的二維數(shù)值波浪水池，并對(duì)非線性Jonswap波浪譜進(jìn)行了模擬。但這些數(shù)值水池研究主要集中于單一波浪運(yùn)動(dòng)規(guī)律的研究，而對(duì)于風(fēng)-浪-流三者之間相互耦合作用研究較少。本文基于FLUENT軟件，利用UDF二次開(kāi)發(fā)，構(gòu)建了一種適于三維海洋工程的數(shù)值風(fēng)-浪-流水池，成功模擬了風(fēng)-浪-流之間的相互耦合作用，并探討了風(fēng)和海流對(duì)波浪的影響。

1 數(shù)值風(fēng)-浪-流水池?cái)?shù)學(xué)模型

1.1 控制方程

假定數(shù)值風(fēng)-浪-流水池所涉及的流體(空氣和海水)均為不可壓縮黏性流體，則數(shù)值風(fēng)-浪-流水池連續(xù)性方程和動(dòng)量方程分別為

(1)

(2)

波浪傳播過(guò)程中，其氣-液交界面時(shí)刻發(fā)生改變，因而需要準(zhǔn)確追蹤自由液面。本文采用VOF方法捕捉自由液面，即通過(guò)引入單元內(nèi)各流體體積分?jǐn)?shù)aq來(lái)捕捉氣液兩相自由液面。當(dāng)aq=0時(shí)，表示該網(wǎng)格單元內(nèi)沒(méi)有q相流體;當(dāng)aq=1時(shí)，表示該網(wǎng)格單元內(nèi)充滿q相流體;當(dāng)0

(3)

式(3)中：aq為q相流體在網(wǎng)格內(nèi)所占體積分?jǐn)?shù)，當(dāng)q=1、2時(shí)，aq分別表示空氣、水相體積分?jǐn)?shù)。

1.2 造波方法

數(shù)值造波方法對(duì)于數(shù)值水池準(zhǔn)確模擬波浪傳播特性具有重要作用，現(xiàn)階段主要的數(shù)值造波方法有仿物理造波法(推板造波法和搖板造波法)[13]、邊界造波法[14]和源項(xiàng)造波法[15]等,其中仿物理造波法原理簡(jiǎn)單，造波功能完善，但對(duì)計(jì)算區(qū)域內(nèi)網(wǎng)格質(zhì)量要求較高，計(jì)算時(shí)間長(zhǎng)；源項(xiàng)造波法產(chǎn)生的波浪效果好，波形穩(wěn)定，但計(jì)算量大，且不易收斂；邊界造波法通過(guò)在波浪入射邊界給定波浪水質(zhì)點(diǎn)速度和波面高度，從而生成數(shù)值波浪，具有較好的適用性并在海洋工程波浪領(lǐng)域應(yīng)用較為廣泛。因此，采用邊界造波法將風(fēng)-浪-流水池入口設(shè)置為風(fēng)-浪-流相互作用的起始位置，各速度場(chǎng)方程可由Airy線性波浪和梯度風(fēng)場(chǎng)理論描述，其中波面方程為

(4)

式(4)中：η為波面高度，m；H為波高，m；k為在2π相位角內(nèi)波的個(gè)數(shù)，rad/m；w為2π相位角內(nèi)波的振動(dòng)次數(shù)，rad/s；x為X方向位移,m。

風(fēng)-浪-流水池入口速度包括自由液面以下水質(zhì)點(diǎn)運(yùn)動(dòng)速度和自由液面以上空氣運(yùn)動(dòng)速度，其中入口水質(zhì)點(diǎn)速度由波浪水質(zhì)點(diǎn)速度和水流流動(dòng)速度兩部分組成，可根據(jù)波流共同作用下流速場(chǎng)理論進(jìn)行疊加[16]，即

ux=uwx+uc=

(5)

(6)

式(5)、(6)中：ux、uz分別為水質(zhì)點(diǎn)X和Z方向速度，m/s；uwx為波浪水質(zhì)點(diǎn)速度，m/s；z表示Z方向上的位移，m；uc為水流速度，m/s;T為波浪周期，s；d為水深，m。

自由液面之上的空氣流動(dòng)速度可采用指數(shù)梯度風(fēng)進(jìn)行定義[9]，即

(7)

1.3 消波方法

采用數(shù)值風(fēng)-浪-流水池進(jìn)行模擬時(shí)，當(dāng)波浪傳播至水池末端時(shí)，波浪會(huì)產(chǎn)生反射，對(duì)前進(jìn)波浪產(chǎn)生干擾。因此，為準(zhǔn)確獲得高質(zhì)量、高精度的波浪，需對(duì)反射波進(jìn)行消波處理。常用的消波方法主要有多孔介質(zhì)消波法和動(dòng)量阻尼消波法，其中多孔介質(zhì)消波法通過(guò)在消波段設(shè)置一個(gè)多孔介質(zhì)區(qū)域，實(shí)現(xiàn)消波目的，但關(guān)于多孔介質(zhì)模型的選取及參數(shù)的設(shè)置較為復(fù)雜；動(dòng)量阻尼消波法通過(guò)在消波段的流體動(dòng)量方程中添加一個(gè)黏性阻尼衰減項(xiàng)，從而實(shí)現(xiàn)消波目的。消波段動(dòng)量方程[10]為

(8)

(9)

式(8)、(9)中：θ為動(dòng)量阻尼消波系數(shù)；x0和x1分別為消波段前端和尾端X方向坐標(biāo)，m；θ1為動(dòng)量阻尼消波經(jīng)驗(yàn)系數(shù)，取值10。

2 數(shù)值風(fēng)-浪-流水池模擬及分析

2.1 網(wǎng)格模型及參數(shù)設(shè)置

采用GAMBIT建立三維數(shù)值風(fēng)-浪-流水池幾何模型并劃分網(wǎng)格，其整體尺寸為：長(zhǎng)50 m、寬20 m、高15 m、水深為10 m、右端消波段長(zhǎng)度為8 m，如圖1所示。整個(gè)數(shù)值風(fēng)-浪-流水池網(wǎng)格采用漸進(jìn)式六面體網(wǎng)格，越靠近水池底端、右端和上端邊界處網(wǎng)格越稀疏，同時(shí)對(duì)自由液面附近兩倍波高區(qū)域進(jìn)行加密處理，便于更好地捕捉波形變化。

圖1 數(shù)值風(fēng)-浪-流水池幾何模型網(wǎng)格劃分

數(shù)值風(fēng)-浪-流水池左邊界設(shè)置為速度入口，氣液兩相速度值和體積組分?jǐn)?shù)可根據(jù)氣液兩相運(yùn)動(dòng)速度場(chǎng)和波面運(yùn)動(dòng)方程并采用UDF進(jìn)行定義;右邊界設(shè)置為自由出流邊界，上邊界設(shè)置為對(duì)稱邊界，其余邊界為固壁邊界。造波過(guò)程中，采用三維非定常分離隱式求解器求解，其中湍流模型為RNGk-ε模型，壓力速度之間的耦合作用采用PISO算法，離散格式為中心差分格式，壓力方程采用加權(quán)體積力格式，動(dòng)量方程、湍動(dòng)能和湍動(dòng)能耗散率方程均采用二階迎風(fēng)格式。壓力參考值為一個(gè)標(biāo)準(zhǔn)大氣壓，重力加速度為-9.81 m/s2，時(shí)間步長(zhǎng)為0.01 s。

2.2 模型驗(yàn)證

為了驗(yàn)證數(shù)值風(fēng)-浪-流水池的準(zhǔn)確性和適用性，采用線性Airy波浪對(duì)其進(jìn)行驗(yàn)證。設(shè)定模擬的線性Airy波浪參數(shù)為：波長(zhǎng)L=8 m，波高H=0.40 m，波陡λ=0.05，周期T=2.263 6 s，波速c=3.534 2 m/s。圖2為該線性Airy波浪生成和傳播過(guò)程示意圖，上部藍(lán)色部分為空氣，下部紅色部分為水，中間綠色部分為波面，可以看出t=20 s時(shí)數(shù)值水池已形成穩(wěn)定的Airy波面，同時(shí)水池右端采用動(dòng)量阻尼消波效果較好，消波段波面基本保持不變。

圖2 Airy波浪生成和傳播過(guò)程

分別在數(shù)值風(fēng)-浪-流水池中不同位置處設(shè)置浪高儀，實(shí)時(shí)監(jiān)測(cè)不同位置處波面時(shí)程。x=2 m和x=18 m位置處波面時(shí)程如圖3所示,可以看出，兩個(gè)位置處波面時(shí)程均較為穩(wěn)定，且與理論波面時(shí)程變化規(guī)律基本一致,其中x=2 m處數(shù)值波面時(shí)程與理論波面時(shí)程吻合度較高,而x=18 m處數(shù)值波面時(shí)程峰值略小于理論波面峰值,這主要是由于波浪傳播過(guò)程中由于黏滯耗能等影響而導(dǎo)致波浪最大幅值呈現(xiàn)一定程度的衰減。

圖3 波浪數(shù)值波形和理論波形比較

進(jìn)一步統(tǒng)計(jì)該數(shù)值模擬波浪的波長(zhǎng)、波幅和周期等特征，并與理論值進(jìn)行比較分析(表1)，可以看出：數(shù)值波長(zhǎng)平均值為8.301 m，與理論波長(zhǎng)相對(duì)誤差(取絕對(duì)值，余同)為3.76%；x=2 m處，模擬波幅為0.196 m，與理論波幅相對(duì)誤差為2.00%；隨著x取值逐漸遠(yuǎn)離造波邊界，模擬波幅與理論波幅相對(duì)誤差絕對(duì)值逐漸增大，這主要是由于波浪傳播過(guò)程中會(huì)存在逐步的衰減作用，但波幅相對(duì)誤差均小于10%，且不同位置處波浪周期與理論周期相對(duì)誤差均小于1%，充分驗(yàn)證了該數(shù)值風(fēng)-浪-流水池模擬波浪的可靠性。

表1 波長(zhǎng)、波幅和周期模擬值與理論值相對(duì)誤差統(tǒng)計(jì)

2.3 浪-流耦合作用模擬及分析

根據(jù)式(5)～(10)，分別設(shè)定浪-流耦合作用下造波和消波方法，其中速度入口邊界通過(guò)UDF進(jìn)行定義(包括流速場(chǎng)和風(fēng)場(chǎng))，消波段通過(guò)增加UDF添加動(dòng)量阻尼源項(xiàng)進(jìn)行消波，然后利用數(shù)值風(fēng)-浪-流水池模擬浪-流耦合作用。為了便于比較，設(shè)定浪-流耦合作用幾何模型、邊界條件、物理屬性以及波浪參數(shù)等均與模型驗(yàn)證時(shí)保持一致；海流速度共設(shè)置4組，分別為uc=0 、-0.1 、 0.1 、0.3 m/s。

圖4為4種不同海流流速下x=2 m處波面時(shí)程曲線，可以看出，不同流速下波面時(shí)程曲線形狀基本保持一致，但流速越大波面時(shí)程形成穩(wěn)定波形所需時(shí)間越短，即uc=0.3 m/s時(shí)波形最先達(dá)到穩(wěn)定，uc=-0.1 m/s時(shí)波浪最后達(dá)到穩(wěn)定狀態(tài)。這主要是由于波流同向時(shí)，海流對(duì)波浪傳播具有推動(dòng)作用，導(dǎo)致波浪傳播速度加快；而波流逆向時(shí)，海流對(duì)波浪傳播具有阻礙作用。

圖4 浪-流耦合作用下x=2 m處波面時(shí)程對(duì)比

分別提取4種不同流速下的波長(zhǎng)和波幅，并同波流耦合作用下波長(zhǎng)和波幅理論計(jì)算公式[16]進(jìn)行比較分析，結(jié)果見(jiàn)表2。由表2知，4種流速下數(shù)值波長(zhǎng)和波幅與理論值吻合較好，其中數(shù)值波長(zhǎng)和理論波長(zhǎng)最大相對(duì)誤差為3.76%，數(shù)值波幅和理論波幅最大相對(duì)誤差為5.35%。對(duì)比4種流速下波浪峰值大小可知，uc=-0.1 m/s時(shí)波浪峰值最大，uc=0.3 m/s時(shí)波浪峰值最小，這主要是由于波流同向時(shí)，波浪傳播速度增大，波形被拉長(zhǎng)，導(dǎo)致波浪峰值減?。欢髂嫦驎r(shí)，波浪傳播速度減小，波形被壓縮，導(dǎo)致波浪峰值增大所致。

表2 浪-流耦合作用下波長(zhǎng)和波幅模擬值與理論值

2.4 風(fēng)-浪-流耦合作用模擬及分析

實(shí)際海洋環(huán)境中通常風(fēng)、浪、流三者是共同存在的，進(jìn)一步基于該數(shù)值風(fēng)-浪-流水池模擬風(fēng)-浪-流耦合作用。根據(jù)式(7)在造波邊界自由液面上方空氣區(qū)域施加梯度風(fēng)風(fēng)場(chǎng)，同時(shí)海流速度uc取值0.1 m/s，波浪及其他物理參數(shù)與上述完全相同。圖5為考慮波浪、浪-流耦合以及風(fēng)-浪-流耦合等3種情況下，x=2 m位置處波面時(shí)程對(duì)比，可以看出，風(fēng)-浪-流耦合作用下波面達(dá)到穩(wěn)定時(shí)間和浪-流耦合作用下達(dá)到穩(wěn)定時(shí)間差不多，但比單純波浪達(dá)到穩(wěn)定時(shí)間快；風(fēng)-浪-流耦合作用下波面峰值附近呈現(xiàn)一定的傾角，這主要是由于海風(fēng)使得波面峰值附近波浪運(yùn)動(dòng)速度加快，導(dǎo)致波峰發(fā)生傾斜，加速波浪發(fā)生碎浪現(xiàn)象。進(jìn)一步統(tǒng)計(jì)風(fēng)-浪-流耦合作用下波浪波長(zhǎng)和波幅，其平均值分別為8.836 m和0.185 m，相比于浪-流耦合作用下的波長(zhǎng)和波幅分別增大2.45%和減小3.14%。

圖5 風(fēng)-浪-流耦合作用下x=2 m處波面時(shí)程對(duì)比

3 結(jié)論

1) 基于FLUENT軟件，通過(guò)UDF二次開(kāi)發(fā)，構(gòu)建了適于海洋工程的數(shù)值風(fēng)-浪-流水池，實(shí)現(xiàn)了數(shù)值風(fēng)-浪-流水池造波和消波功能，成功模擬了浪-流以及風(fēng)-浪-流相互耦合作用，模擬結(jié)果與理論結(jié)果吻合度較好，驗(yàn)證了該數(shù)值風(fēng)-浪-流水池具有較好的精度。

2) 浪-流耦合作用數(shù)值模擬結(jié)果表明，海流對(duì)于波浪運(yùn)動(dòng)參數(shù)具有較大影響。波流同向時(shí)，海流導(dǎo)致波浪傳播速度和波長(zhǎng)增大，幅值減??；波流逆向時(shí)，波浪傳播速度和波長(zhǎng)均減小，波幅增大，且海流速度越大，波浪參數(shù)增大或減小的幅度越明顯。

3) 風(fēng)-浪-流耦合作用數(shù)值模擬結(jié)果表明，風(fēng)會(huì)導(dǎo)致波浪傳播過(guò)程中波長(zhǎng)增大，波幅減小，同時(shí)會(huì)導(dǎo)致波峰呈現(xiàn)一定傾角，加速波浪發(fā)生破碎。