基于二階規(guī)則波模擬的數(shù)值波浪水槽優(yōu)化

曾 誠,尹雨然,周 婕,邱 斐,白一墨,王玲玲

(1.河海大學水利水電學院,江蘇 南京 210098; 2.河海大學力學與材料學院,江蘇 南京 211100)

近年來,隨著海洋資源的不斷開發(fā),海岸和近海工程受到社會各界的廣泛關(guān)注。波浪作為海洋環(huán)境中的重要因素,對海洋工程結(jié)構(gòu)的強度和穩(wěn)定性具有重要影響。數(shù)值模擬作為一種研究手段,在工程應用中得到了越來越多的重視[1]。目前,波浪研究中較為常用的數(shù)值模型是基于Navier-Stokes方程的黏性流數(shù)值波浪水槽[2]和基于光滑粒子流體動力學的數(shù)值波浪水槽[3]。與物理波浪水槽[4]相比,數(shù)值波浪水槽具有成本低、信息完整和模擬能力強等優(yōu)點。

數(shù)值造波和消波技術(shù)是數(shù)值波浪水槽構(gòu)建的關(guān)鍵技術(shù)。常用的造波方法包括速度邊界法[5]、推板造波法[6]和質(zhì)量源造波法[7]等,消波方法主要包括輻射邊界法、阻尼消波法和松弛區(qū)域法等。數(shù)值造波技術(shù)方面,徐剛等[8]采用速度邊界法建立了數(shù)值波浪水槽,研究了波浪衰減特性,并對數(shù)值波浪水槽的模擬精度和波浪衰減問題進行了分析。Higuera等[9]利用OpenFOAM建立了三維多推板數(shù)值波浪水槽,成功模擬生成了二維和三維聚焦波。Prasad等[10]采用推板造波法建立了數(shù)值波浪水槽,并對波能轉(zhuǎn)換器的性能進行了研究。Machado等[11]將推板和搖板2種造波方法進行了對比,認為推板造波法更準確。Lin等[12]將源造波法運用到基于Navier-Stokes方程的數(shù)值波浪水槽中,通過使用不同的源函數(shù),產(chǎn)生各種波序列。為消除波浪水槽邊界產(chǎn)生的反射波,保證水槽內(nèi)波浪穩(wěn)定,前人提出了多種消波方法。Milgram[13]在1970年最早提出了主動吸收造波理論。William等[14]在水槽末端設(shè)置斜坡來達到消波目的。Kim等[15]采用網(wǎng)格逐漸稀疏的消波方法建立了三維數(shù)值波浪水槽,產(chǎn)生了較好的消波效果。另外,還可采用出流邊界使波浪完全透射過出流邊界,從而避免發(fā)生波浪反射,其中Sommerfeld輻射邊界[16]具有較好的消波效果。Clément[17]對Sommerfeld輻射邊界條件進行了研究,指出由于該條件在時間和空間上是局部的,因此僅限于頻率已知的規(guī)則入射波和長波情況,而不適用于完全非定常入射波情況。鄒志利等[18-19]也對Sommerfeld輻射邊界條件進行了研究,發(fā)現(xiàn)單獨采用該條件的消波效果并不理想。Israeli等[20]提出了阻尼消波法,分析了消波系數(shù)對消波效果的影響,并嘗試通過輻射邊界條件和阻尼區(qū)消波相結(jié)合的方法提高消波效果。Baker等[21]在其二維邊界元模式中引入阻尼消波法,其主要思想是在消波區(qū)域?qū)ψ杂擅孢吔鐥l件或動量方程添加一阻尼耗散項,從而達到消波的目的。隨著阻尼消波法的不斷發(fā)展和完善,該方法成為了數(shù)值波浪水槽研究中廣泛運用的消波方法。Clément[17]指出,通過阻尼消波能夠有效消除高頻短波,但是對于低頻長波消波效果不太理想,需要延長阻尼區(qū)域,增加計算成本。Méndez等[22]針對阻尼消波法進行了研究,結(jié)果表明阻尼消波區(qū)設(shè)置不合理會導致平均水位的緩慢抬升。韓朋[23]基于VOF方法針對不規(guī)則波的數(shù)值模擬進行了阻尼消波研究,給出了不同波要素組合情況下阻尼系數(shù)的經(jīng)驗取值范圍。王國玉等[24]基于OpenFOAM求解器interDyMFoam,開發(fā)了數(shù)值波浪水槽,并將計算結(jié)果與試驗測量數(shù)據(jù)進行對比,驗證了數(shù)值造波和消波方法的可靠性。另外,Tian等[25-26]給出了阻尼項的不同表達形式。近年來,Jacobsen等[27]提出了松弛區(qū)域法,在波浪水槽前后端設(shè)置松弛區(qū)域,通過合理分配計算值和理論值之間的權(quán)重實現(xiàn)數(shù)值解與解析解之間的過渡,從而達到消波效果。田康等[28]基于OpenFOAM建立數(shù)值波浪水槽,并通過二次開發(fā)引入松弛區(qū)域法,取得了較好的消波效果。

基于上述文獻綜述,前人對非線性波模擬的數(shù)值波浪水槽優(yōu)化設(shè)置研究較少,本文基于CFD開源軟件包OpenFOAM建立了二維黏性流數(shù)值波浪水槽,綜合考慮不同波陡工況下(i=0.01～0.06)Stokes二階波的模擬精度和計算成本,從造、消波方式選擇、網(wǎng)格劃分、紊流模型選擇等方面對數(shù)值波浪水槽進行全面優(yōu)化,以期提高非線性波模擬的精度和效率。

1 數(shù)學模型

1.1 基本兩相流控制方程

基本兩相流控制方程主要包括不可壓縮流體的連續(xù)性方程(式(1))、動量守恒方程(式(2))以及基于體積函數(shù)(volume of fluid,VOF)法的相方程(式(3))。

?·U=0

(1)

(2)

(3)

式中:U=(u,w)為笛卡爾坐標系下的流速場;ρ為整個流場的密度;μ為動力黏滯系數(shù);g=(0,-9.81)為重力加速度;p為壓力場;α為水的體積分數(shù)。

當α=1,表示網(wǎng)格單元內(nèi)全部為水,當α=0,表示網(wǎng)格內(nèi)全為空氣,當0<α<1,代表該網(wǎng)格單元內(nèi)存在自由液面。通過對α的迭代求解,可根據(jù)式(4)進一步確定網(wǎng)格單元內(nèi)其他的特性參數(shù)。

φ=αφW+(1-α)φA

(4)

式中:φ為網(wǎng)格內(nèi)其他參數(shù),如密度、分子黏性系數(shù)等;下標W表示水,A表示空氣。

本文采用k-ωSST紊流模型封閉控制方程,采用PIMPLE算法進行速度-壓力解耦。其中,時間項采用一階隱式格式離散,對流項和黏性項采用二階迎風格式離散,壓力項、紊動能和耗散率項采用二階中心差分格式進行離散。

1.2 造波方法

1.2.1 推板造波法

在數(shù)值波浪水槽中采用推板造波法時,推板沖程S與目標波高H的關(guān)系為

(5)

其中k=2π/L

式中:k為波數(shù);L為波長;d為初始水深。

本文模擬的波浪中包括深水(水深大于波長一半)和有限水深(水深和波長關(guān)系為0.5>h/L>0.05)條件下的有限振幅波,因此采用Stokes非線性波理論。根據(jù)Madsen[29]造波理論,采用修正后的推板位移公式產(chǎn)生穩(wěn)定的Stokes二階波:

(6)

式中:xp為推板沖程;ω為圓頻率;T為波浪周期;t為時間。

本文理論解采用Stokes二階波波面方程:

(7)

式中x為沿波浪傳播方向的坐標。

1.2.2 速度邊界法

采用速度邊界法造波時,通過控制邊界處的水平和垂向速度產(chǎn)生目標波形,水平速度u和垂向速度w分別為

(8)

(9)

式中z為沿水深方向坐標,向上為正。

1.3 消波方法

1.3.1 阻尼消波法

阻尼消波的原理為在動量方程中加入阻尼項,在數(shù)值水槽消波段實現(xiàn)動量衰減,從而消除邊界上的反射波,在式(2)中加入阻尼項ρθU后,動量方程如下:

(10)

式中:x0為消波段起點;x1為消波段終點;θ為阻尼函數(shù),常用的阻尼函數(shù)形式包括線性形式(r=1)和平方形式(r=2);θ1為阻尼系數(shù);r為阻尼項指數(shù)。

1.3.2 輻射邊界法

在水槽末端施加適當?shù)妮椛錀l件,能使波浪在出口邊界處無干擾或干擾很小的傳播出去。Sommerfeld輻射邊界條件的定義如下:

(11)

式中Cu、Cw分別為u和w在出口邊界處的值。

1.3.3 松弛區(qū)域法

在數(shù)值波浪水槽首端設(shè)置松弛區(qū)域可以有效消除計算區(qū)域內(nèi)部的波浪反射對入口的影響,在水槽末端設(shè)置松弛區(qū)域能夠防止出口邊界產(chǎn)生反射波,對研究區(qū)域產(chǎn)生影響。松弛區(qū)域法消波原理是在松弛區(qū)域內(nèi)引入一個消波因子αR:

(12)

式中XR為消波因子的局部坐標。

XR軸正方向始終指向計算域外,這樣使得消波因子在靠近流體內(nèi)部一側(cè)的變化非常緩慢,從而達到減少反射的目的。在松弛區(qū)域內(nèi),流速和體積分數(shù)等物理量(φ)由計算值(φc)和解析解(φa)加權(quán)得到:

φ=αRφc+(1-αR)φa

(13)

2 模型建立

2.1 模型區(qū)域設(shè)置

本文建立的二維數(shù)值波浪水槽計算區(qū)域、計算邊界及松弛區(qū)域內(nèi)松弛因子分布如圖1所示。水槽模型長度為6倍波長以上,其中造波區(qū)長度為5m,研究區(qū)長度為4～5倍波長,消波區(qū)為1～2倍波長,采用松弛區(qū)域法消波時,進口松弛區(qū)域長度為5m,根據(jù)初步模擬結(jié)果,出口松弛區(qū)域為1倍波長較為合理。模型頂部設(shè)置為大氣邊界,左側(cè)邊界為造波邊界,水槽底部設(shè)置為無滑移壁面,出口邊界根據(jù)消波方法的不同分別設(shè)置為無滑移壁面或輻射邊界。

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

2.2 網(wǎng)格無關(guān)性分析

為保證模擬精度,本文選取三套網(wǎng)格進行網(wǎng)格無關(guān)性分析,具體網(wǎng)格參數(shù)見表1。目標波高設(shè)置為H=0.1m,波長L=5.03m,采用推板造波,基于平方形式的阻尼函數(shù)消波。三套網(wǎng)格的計算結(jié)果如圖2所示。由圖2可知,mesh1的模擬值與理論值相差較大,模擬精度較低,而mesh2和mesh3的模擬結(jié)果相近,并且十分接近理論值,都能達到較高的模擬精度。為保證模擬精度的同時還節(jié)省計算資源,后續(xù)計算采用mesh2的網(wǎng)格設(shè)置方式,并且在波面附近進行網(wǎng)格加密,以達到更好的模擬效果。三套網(wǎng)格及后續(xù)模型的計算時間均為20T,計算時間步長取T/1000。

表1 網(wǎng)格參數(shù)

圖2 不同網(wǎng)格下的沿程波面

2.3 紊流模型對比

為比較紊流模型對黏性流數(shù)值波浪水槽計算結(jié)果的影響,本文對比分析了目前比較常用的k-ε紊流模型和k-ωSST紊流模型,模擬的波浪參數(shù)及采用的造波、消波方法與第2.2節(jié)相同。由圖3可知,k-ωSST紊流模型的模擬精度更接近理論值,而k-ε紊流模型的計算結(jié)果與理論值的偏差相對較大,波浪衰減較為明顯。因此,本文后續(xù)計算中均選擇k-ωSST模型。

圖3 不同紊流模型的沿程波面模擬結(jié)果

3 計算結(jié)果分析

3.1 工況設(shè)置

本文對6種不同波陡的Stokes二階波進行模擬計算,波陡范圍為0.01～0.06。各工況的波浪參數(shù)及模擬計算中采用阻尼消波法時的θ1見表2。由表2可知,θ1隨著i的增大而增大。

表2 波浪參數(shù)

3.2 消波效果

對表2中所示不同波陡工況,采用推板造波法產(chǎn)生Stokes二階波,分別采用輻射邊界法、阻尼消波法和松弛區(qū)域法進行模擬計算,對其消波效果進行對比分析。其中,在采用阻尼消波法時,對采用線性形式阻尼函數(shù)(r=1)和平方形式函數(shù)(r=2)的消波效果進行對比分析。本文采用Ursell等[30]提出反射系數(shù)定量評價消波效果,一般認為反射系數(shù)在5%以內(nèi)可以滿足消波要求。

由表3可知,阻尼消波的2種阻尼函數(shù)形式的反射系數(shù)均在5%以內(nèi),都具有較好的消波效果,但是隨著波陡的增加,平方形式的阻尼函數(shù)優(yōu)勢變得更加明顯。而采用輻射邊界法的消波效果僅僅在波陡較小的情況下(i=0.01和i=0.02),反射系數(shù)才滿足要求,而在其他波陡情況下(i=0.03～0.06),反射系數(shù)均大于5%,采用松弛區(qū)域法計算得到的反射系數(shù)均保持在較小的水平,但是在大波陡工況下,松弛區(qū)域法不如平方形式的阻尼函數(shù)優(yōu)勢明顯。在不同波陡情況下,線性形式和平方形式的阻尼函數(shù)能夠保證較好的消波效果,而平方形式的阻尼函數(shù)相較于線性形式,消波效果上更有優(yōu)勢,并且從阻尼在消波區(qū)域內(nèi)的分布來看,平方形式的阻尼函數(shù)更有利于保證研究區(qū)域內(nèi)的模擬精度,阻尼分布更加合理。后續(xù)進行造波效果分析時,消波方法均采用基于平方形式阻尼函數(shù)的阻尼消波法。

表3 反射系數(shù)計算結(jié)果

3.3 造波效果

為對模型造波效果進行對比分析,對表2中所示不同波陡工況分別采用推板造波法和速度邊界法進行Stokes二階波的模擬計算,消波方法均采用基于平方形式阻尼函數(shù)的阻尼消波法。針對不同工況計算結(jié)果,由于所計算的工況波高差距較大,為了得到更合理的誤差分析結(jié)果,首先將計算值與理論值通過min-max標準化方法進行歸一化處理,再分別計算不同工況下的均方根誤差,計算結(jié)果見表4。

表4 均方根誤差及計算耗時

由表4可知,在波陡較小的情況下,波浪的非線性較弱,利用速度邊界法具有較高的模擬精度,在波陡為0.01的情況下,速度邊界法的模擬精度高于推板造波法。但是隨著波陡的增加,波浪的非線性逐漸增強,速度邊界逐漸變得不適用,模擬精度逐漸低于推板造波法。在波陡為0.05和0.06時,2種造波方法產(chǎn)生的波浪與理論值的均方根誤差相差較大,在模擬精度上推板造波法的優(yōu)勢較為明顯。

本文的模擬計算在河海大學高性能計算平臺上進行,每個工況利用48核心并行計算,CPU類型為2顆Intel Xeon 6240R(24核2.4GHz)。雖然推板造波能夠產(chǎn)生較為精確的波形,但是在計算時間方面,采用推板造波法的計算耗時是速度邊界法的5～8倍,尤其是在波陡較大,所需網(wǎng)格數(shù)量較多的情況下,計算耗時的差別會更加明顯(表4)。

圖4為6種波陡工況下,采用2種造波方法的沿程波面計算結(jié)果與理論值的對比。波陡過小或過大都不容易產(chǎn)生與理論值吻合程度較高的結(jié)果。在小波陡工況(W1和W2)下,雖然計算值與理論值的均方根誤差表現(xiàn)較小,但是在波形方面,2種造波方法所產(chǎn)生的波浪與理論值都存在一定的偏差,而此時速度邊界法所產(chǎn)生的波形在精度方面較有優(yōu)勢。在大波陡工況(W5和W6)下,波浪非線性較強,速度邊界法的優(yōu)勢不明顯,推板造波法具有較高的精度。此外,在W3和W4工況下,2種造波方法所產(chǎn)生的波形與理論值吻合較好,采用速度邊界法和推板造波法均能達到較高的模擬精度。

圖4 不同波陡條件下水槽內(nèi)波高

為分析推板造波法和速度邊界法所產(chǎn)生的波浪隨時間變化的情況,在距離水槽造波邊界2倍波長處設(shè)置監(jiān)測點,監(jiān)測波高隨時間變化的情況。圖5為采用2種造波方法達到穩(wěn)定波高的時間序列(15s后)與理論值的比較情況。由圖5可知,采用2種造波方法時監(jiān)測點處波面歷時曲線與理論值基本吻合,表現(xiàn)為小波陡工況(W1、W2和W3)下吻合程度較高,大波陡工況(W4、W5和W6)下吻合較差,誤差主要體現(xiàn)在波峰處。另外,隨著計算時間的增長,采用速度邊界法和阻尼消波法的組合方式時出現(xiàn)了較為明顯的水面抬升現(xiàn)象,如圖5(b)～(f)所示。波陡越大,水面抬升越明顯,在W6工況(i=0.06)下,波面最大約抬升了波高的10%。出現(xiàn)水面抬升現(xiàn)象的原因是由于大波陡工況下,波浪的非線性較強,在一個波浪周期內(nèi)邊界上輸入與抽取的速度不等,會引起質(zhì)量守恒方面的問題,隨著計算時間的增加,就產(chǎn)生了波面抬升問題[31]。然而,采用推板造波法和阻尼消波法的組合方式時,在各個波陡工況下,均能保持波面穩(wěn)定,沒有出現(xiàn)水面抬升現(xiàn)象。

圖5 不同波陡條件下x=2L處波面歷時曲線

在進一步的模擬研究中發(fā)現(xiàn),模擬非線性較強波浪時,如采用速度邊界法造波,則必須采取相應的修正措施。例如,采用松弛區(qū)域法,在數(shù)值水槽的前后端添加松弛區(qū),既能保證質(zhì)量守恒又消除了二次反射波。圖6為大波陡工況(i=0.06)下采用速度邊界法造波時,分別采用阻尼消波法和松弛區(qū)域法消波的計算結(jié)果與理論值的波面歷時曲線對比,從圖6可看出,松弛區(qū)域法能夠有效避免波面抬升現(xiàn)象,模擬結(jié)果與理論值符合較好。

圖6 大波陡工況(i=0.06)時速度邊界法修正效果對比

4 結(jié) 論

a.消波效果方面,采用基于平方形式阻尼函數(shù)的阻尼消波法具有較好的消波效果。采用松弛區(qū)域法可有效消除大波陡工況下采用速度邊界法導致的水面抬升現(xiàn)象。

b.造波效果方面,速度邊界法在小波陡工況下模擬精度較高,計算耗時較少。推板造波法在大波陡工況下模擬精度較高,波形穩(wěn)定性較好。

c.根據(jù)本文二階規(guī)則波的數(shù)值模擬計算結(jié)果,對于小波陡(i=0.01～0.03)工況采用速度邊界法造波和松弛區(qū)域法消波的組合方式具有較高的模擬精度,對于大波陡(i=0.04～0.06)工況采用推板造法造波和基于平方形式阻尼函數(shù)的阻尼消波法消波的組合方式具有較高的模擬精度。