海洋工程數(shù)值風(fēng)-浪-流水池模擬及分析*

劉紅兵 陳國明 呂 濤 黃 翱

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

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

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

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

1.1 控制方程

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

(1)

(2)

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

(3)

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

1.2 造波方法

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

(4)

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

風(fēng)-浪-流水池入口速度包括自由液面以下水質(zhì)點(diǎn)運(yùn)動速度和自由液面以上空氣運(yùn)動速度，其中入口水質(zhì)點(diǎn)速度由波浪水質(zhì)點(diǎn)速度和水流流動速度兩部分組成，可根據(jù)波流共同作用下流速場理論進(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。

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

(7)

1.3 消波方法

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

(8)

(9)

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

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

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

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

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

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

2.2 模型驗(yàn)證

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

圖2 Airy波浪生成和傳播過程

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

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

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

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

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

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

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

圖4 浪-流耦合作用下x=2 m處波面時程對比

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

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

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

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

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

3 結(jié)論

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

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

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