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

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

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

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

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

1 數學模型

1.1 基本兩相流控制方程

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

?·U=0

(1)

(2)

(3)

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

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

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

(4)

式中:φ為網格內其他參數,如密度、分子黏性系數等;下標W表示水,A表示空氣。

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

1.2 造波方法

1.2.1 推板造波法

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

(5)

其中k=2π/L

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

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

(6)

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

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

(7)

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

1.2.2 速度邊界法

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

(8)

(9)

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

1.3 消波方法

1.3.1 阻尼消波法

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

(10)

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

1.3.2 輻射邊界法

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

(11)

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

1.3.3 松弛區(qū)域法

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

(12)

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

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

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

(13)

2 模型建立

2.1 模型區(qū)域設置

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

圖1 數值波浪水槽示意圖

2.2 網格無關性分析

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

表1 網格參數

圖2 不同網格下的沿程波面

2.3 紊流模型對比

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

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

3 計算結果分析

3.1 工況設置

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

表2 波浪參數

3.2 消波效果

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

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

表3 反射系數計算結果

3.3 造波效果

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

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

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

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

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

圖4 不同波陡條件下水槽內波高

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

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

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

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

4 結 論

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

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

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