姚建喜 吳靜萍

(高性能船舶技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室1) 武漢 430063) (內(nèi)河航運(yùn)技術(shù)湖北省重點(diǎn)實(shí)驗(yàn)室2) 武漢 430063) (武漢理工大學(xué)交通學(xué)院3) 武漢 430063)

孤立波數(shù)值造波研究*

姚建喜1,2,3)吳靜萍1,2,3)

(高性能船舶技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室1)武漢 430063) (內(nèi)河航運(yùn)技術(shù)湖北省重點(diǎn)實(shí)驗(yàn)室2)武漢 430063) (武漢理工大學(xué)交通學(xué)院3)武漢 430063)

基于開源平臺(tái)OpenFOAM研究開發(fā)一種可模擬推板運(yùn)動(dòng)的動(dòng)邊界條件,實(shí)現(xiàn)推板造波數(shù)值模擬.以孤立波為研究對(duì)象,基于Goring孤立波造波理論以及Fenton九階孤立波解,確定兩種推板運(yùn)動(dòng)形式,并據(jù)此對(duì)孤立波造波和傳播過程進(jìn)行模擬.數(shù)值模擬結(jié)果表明，基于Goring方法所得到的孤立波后方水面凹陷現(xiàn)象較基于Fenton方法所得到的孤立波后方水面凹陷現(xiàn)象更為顯著,且波高衰減更快.將數(shù)值模擬結(jié)果與水槽試驗(yàn)結(jié)果進(jìn)行比較,對(duì)兩者之間的差別進(jìn)行分析和討論.

孤立波;OpenFOAM;推板造波;數(shù)值造波;水槽試驗(yàn)

0 引 言

近年來,海嘯頻發(fā)且近岸建筑物增多,研究海浪在近岸水域中的傳播機(jī)理和特性是預(yù)防災(zāi)難、減少損失的一個(gè)重要途徑.孤立波是一種淺水域有限振幅波,與海浪的傳播特性相似,因此通常以孤立波的研究結(jié)果來分析近岸的海浪.自1844年羅素在實(shí)驗(yàn)室中發(fā)現(xiàn)孤立波以來,人們對(duì)孤立波理論進(jìn)行了大量研究,其中如何在實(shí)驗(yàn)室中通過人工方式產(chǎn)生孤立波是研究其傳播特性的一個(gè)重要方面.

文獻(xiàn)[1-3]對(duì)實(shí)驗(yàn)室造孤立波的方法多有報(bào)道,其中Goring方法應(yīng)用較為廣泛.該方法是基于一階孤立波解的推板造波理論,適合波高較小的孤立波,當(dāng)波高較大時(shí),產(chǎn)生的孤立波后方會(huì)出現(xiàn)水面凹陷現(xiàn)象,且在傳播過程中波高發(fā)生衰減.Wu等[4]基于高階孤立波解對(duì)Goring方法進(jìn)行了改進(jìn),獲得的波形得到明顯改善.

通過近30年的發(fā)展,數(shù)值計(jì)算方法逐漸成熟且已經(jīng)成功地應(yīng)用于造波及波浪傳播過程的數(shù)值模擬研究,出現(xiàn)了所謂的“數(shù)值波浪水池”概念,數(shù)值計(jì)算方法的應(yīng)用更有利于觀察水波現(xiàn)象,是當(dāng)前研究波浪傳播特性的重要手段之一.本文基于開源計(jì)算流體動(dòng)力學(xué)平臺(tái)OpenFOAM研究開發(fā)一種動(dòng)邊界條件,使之能模擬推板線運(yùn)動(dòng),從而實(shí)現(xiàn)數(shù)值造波.根據(jù)Goring孤立波造波理論和Fenton九階孤立波解,確定兩種推板運(yùn)動(dòng)形式,對(duì)孤立波的造波及其傳播過程進(jìn)行數(shù)值模擬[5-10].將計(jì)算結(jié)果與水槽試驗(yàn)結(jié)果進(jìn)行比對(duì)以驗(yàn)證數(shù)值模擬的有效性.

1 孤立波推板造波理論

圖1為推板造波示意圖.對(duì)于給定水深h,當(dāng)推板以速度u沿水平方向運(yùn)動(dòng)時(shí),水面即興起波浪.推板運(yùn)動(dòng)速度是時(shí)間t或推板位置ξ的函數(shù).根據(jù)相關(guān)造波理論,當(dāng)推板的運(yùn)動(dòng)速度函數(shù)或位置函數(shù)確定后,便可驅(qū)動(dòng)推板作相應(yīng)運(yùn)動(dòng),獲得目標(biāo)波形.

圖1 推板造波示意圖

根據(jù)Goring孤立波造波理論,推板位置和速度的關(guān)系為

(1)

式中：c為波速；η為波面抬高且

η=Hs2

(2)

其中:H為目標(biāo)波高；s=sech [k(ξ-ct)].k和c的表達(dá)式為

(3)

(4)

當(dāng)水深和目標(biāo)波高確定后,將式(2)代入式(1)并采用數(shù)值方法進(jìn)行求解可獲得推板位置的時(shí)歷曲線.上述方法中的推板位置是根據(jù)一階孤立波解而得到的.

本文的研究同時(shí)采用Fenton推導(dǎo)的九階孤立波解來確定推板位置時(shí)歷曲線.η,k和c的九階解分別為

(5)

(6)

(7)

表1為式(5)～(7)中的系數(shù)ηi,ki和ci或其計(jì)算表達(dá)式. 將式(5)代入式(1),便可求解得到基于九階孤立波解的推板位置時(shí)歷曲線,求解算法可參考文獻(xiàn)[4],這里不再贅述.

表1 九階孤立波解的系數(shù)

2 數(shù)值造波方法

本研究通過求解納維-斯托克斯(Navier-Stokes,N-S)方程來模擬推板運(yùn)動(dòng),實(shí)現(xiàn)數(shù)值造波.這里僅考慮二維情況,在圖1的坐標(biāo)系0-xz下對(duì)流體動(dòng)力學(xué)控制方程進(jìn)行描述.基于不可壓縮假設(shè),流體連續(xù)性方程為

(8)

N-S方程為

(9)

(10)

式中：u和w為流體速度分量；ν為粘性系數(shù)；p為動(dòng)壓；ρ為密度；g為重力加速度.采用流體體積分?jǐn)?shù)法(volume of fluids,VoF)來捕捉自由面,體積分?jǐn)?shù)控制方程為

(1)

式中：F為流體體積分?jǐn)?shù),若F=1為計(jì)算單元充滿水,若F=0為計(jì)算單元充滿空氣,若0

式(9)和(10)中的粘性系數(shù)和密度為

ρ=Fρw+(1-F)ρe

(12)

ν=Fνw+(1-F)νe

(13)

式中：ρw和νw為水的密度和粘性系數(shù)；ρe和νe為空氣的密度和粘性系數(shù).

本文推板數(shù)值造波是在開源CFD平臺(tái)OpenFOAM上進(jìn)行的. OpenFOAM提供了基于有限體積方法的求解器以及多種方程離散格式.盡管OpenFOAM自帶基于網(wǎng)格變形方法的動(dòng)網(wǎng)格功能模塊,但該求解器尚不能用于模擬給定的推板運(yùn)動(dòng).因此,本研究在OpenFOAM的基礎(chǔ)上植入新的動(dòng)邊界條件,實(shí)現(xiàn)推板數(shù)值造波.新開發(fā)的邊界條件,可直接從文件中讀取推板位置時(shí)間歷程,數(shù)值模擬時(shí)推板按照文件中給定的位置時(shí)間歷程作線運(yùn)動(dòng).對(duì)于給定的水深和目標(biāo)波高,首先根據(jù)前述的Goring或Fenton方法確定推板位置時(shí)間歷程,并保存在一個(gè)文件中,然后將它作為求解器的輸入進(jìn)行數(shù)值造波模擬.

3 計(jì)算條件,網(wǎng)格和邊界條件

數(shù)值模擬時(shí),建立的二維虛擬波浪水池尺寸為:0 m

因?yàn)橛?jì)算區(qū)域是一個(gè)長(zhǎng)方形,所以數(shù)值模擬更適合采用結(jié)構(gòu)化網(wǎng)格,圖2為計(jì)算網(wǎng)格的局部視圖,圖中灰色代表水,白色代表空氣,圖中水面處于初始時(shí)的靜止?fàn)顟B(tài).整個(gè)計(jì)算區(qū)域由四個(gè)邊界組成,邊界條件的設(shè)置情況為：①在水池兩端(x=0 m和x=18 m)和水底(z=-0.3 m)邊界上設(shè)置壁面邊界條件；②在水池頂端(z=0.3 m)邊界上設(shè)置空氣入口邊界條件.

圖2 計(jì)算網(wǎng)格局部視圖

數(shù)值模擬時(shí),將水池一端壁面邊界看成是推板,初始時(shí)處于x=0 m的位置,并使之按給定的方式運(yùn)動(dòng)進(jìn)行數(shù)值造波.

4 網(wǎng)格和時(shí)間步長(zhǎng)依賴性分析

為保證模擬精度,首先研究分析波形對(duì)網(wǎng)格疏密程度和時(shí)間步長(zhǎng)的依賴性,這里僅考慮基于Fenton方法得到的波形.為此,生成三個(gè)疏密程度不同的結(jié)構(gòu)化網(wǎng)格,網(wǎng)格1、網(wǎng)格2和網(wǎng)格3.這三個(gè)網(wǎng)格分別由300×10,600×20和1 200×40個(gè)計(jì)算單元組成,且沿水池長(zhǎng)度和高度方向是等間距的. 對(duì)于網(wǎng)格依賴性分析,數(shù)值模擬的時(shí)間步長(zhǎng)為0.001 s.圖3為在孤立波傳播方向上不同位置時(shí)的波面抬高時(shí)歷曲線.

圖3 網(wǎng)格依賴性分析

由圖3可知，在同一位置,隨著網(wǎng)格加密波面抬高峰值和其出現(xiàn)時(shí)間的差異變小;采用三個(gè)網(wǎng)格計(jì)算得到的波面抬高峰值出現(xiàn)時(shí)間之間的差異隨離推板初始位置的距離而逐漸變大,例如，在x=2 m位置時(shí),三條曲線峰值出現(xiàn)時(shí)間較接近,而在x=11 m位置時(shí)差異較大;使用網(wǎng)格3計(jì)算得到的波面抬高峰值與目標(biāo)波高非常接近,表明網(wǎng)格3的密度能保證數(shù)值模擬精度.

對(duì)于時(shí)間步長(zhǎng)依賴性分析,數(shù)值模擬時(shí)采用網(wǎng)格3,且時(shí)間步長(zhǎng)分別為0.001,0.002和0.004 s值模擬結(jié)果見圖4.由圖4可知,時(shí)間步長(zhǎng)對(duì)波面抬高峰值有較大影響,時(shí)間步長(zhǎng)越大孤立波在傳播過程中衰減越快.

圖4 時(shí)間步長(zhǎng)依賴性分析

基于以上網(wǎng)格和時(shí)間步長(zhǎng)依賴性分析,發(fā)現(xiàn)使用網(wǎng)格3和時(shí)間步長(zhǎng)0.001 s時(shí)的數(shù)值模擬結(jié)果已經(jīng)達(dá)到足夠精度.

5 數(shù)值結(jié)果分析和驗(yàn)證

首先根據(jù)Goring方法確定前文計(jì)算工況時(shí)的推板運(yùn)動(dòng)時(shí)間歷程,并將時(shí)間歷程作為求解器的輸入進(jìn)行造波數(shù)值模擬.數(shù)值模擬時(shí)使用網(wǎng)格3,時(shí)間步長(zhǎng)取為0.001 s. 圖5將基于Goring和Fenton方法得到的在孤立波傳播方向上不同位置時(shí)的波面抬高時(shí)歷曲線進(jìn)行了比較.采用兩種方法得到的波面抬高的一個(gè)主要差別在于孤立波傳播過程中Goring方法得到的波形后方水面凹陷更為顯著.值得注意的是在x=11 m時(shí)基于Fenton方法得到的波形后方水面凹陷現(xiàn)象已經(jīng)不明顯.另一方面,由圖5可知，Goring方法得到的波高隨時(shí)間的推移衰減更快,導(dǎo)致波速變得更慢.

圖5 比較基于Goring和Fenton方法的計(jì)算結(jié)果

圖6將基于Goring和Fenton方法得到的在x=1.55 m和x=3.95 m位置時(shí)的波面抬高時(shí)歷曲線與試驗(yàn)曲線進(jìn)行了比較.水槽試驗(yàn)是在流體力學(xué)實(shí)驗(yàn)室完成的,水槽尺寸為18 m×0.6 m×0.8 m.

由圖6可知，計(jì)算結(jié)果與試驗(yàn)結(jié)果表現(xiàn)出良好的一致性,但與計(jì)算曲線相比,試驗(yàn)曲線的峰值都較小,且在孤立波后方出現(xiàn)了更劇烈的水面凹陷和振蕩現(xiàn)象.計(jì)算結(jié)果和試驗(yàn)結(jié)果之間的差異估計(jì)由以下原因造成:

圖6 比較計(jì)算結(jié)果和試驗(yàn)結(jié)果

1) 數(shù)值造波僅考慮二維情況,忽略了水槽側(cè)壁的影響,這可能是導(dǎo)致波高發(fā)生衰減的原因之一.

2) 數(shù)值模擬時(shí),將虛擬水池的一端作為推板,而水槽試驗(yàn)中的推板是獨(dú)立于水槽壁面的,且推板與水槽側(cè)壁和水底存在間隙,盡管間隙很小,但試驗(yàn)時(shí)仍然觀察到水從間隙中濺出.數(shù)值模擬引入的簡(jiǎn)化估計(jì)是導(dǎo)致計(jì)算結(jié)果和試驗(yàn)結(jié)果之間差異的主要原因.

6 結(jié) 論

1) 數(shù)值模擬得到的波高與目標(biāo)波高非常接近,且試驗(yàn)結(jié)果和計(jì)算結(jié)果的比較表明,數(shù)值模擬是有效的.

2) 與Goring方法相比,采用Fenton方法得到的孤立波波高衰減及其后方水面凹陷現(xiàn)象得到改善.

3) 計(jì)算結(jié)果與試驗(yàn)結(jié)果之間的差異估計(jì)由數(shù)值模擬時(shí)引入的簡(jiǎn)化造成,后續(xù)研究工作可考慮真實(shí)工況以進(jìn)一步驗(yàn)證本文方法的精度.

[1]HSIAO S C, LIN T C. Tsunami-like solitary waves impinging and overtopping an impermeable seawall: experiment and RANS modeling [J].Coastal Engineering, 2010,57:759-774.

[2]CHEN Y Y, YANG J H. An experimental study of steep solitary wave reflection at a vertical wall [J].European Journal of Mechanics B/Fluids, 2015,49:20-28.

[3]CHEN Y S, YEH H. Laboratory experimental on counter-propagating collisions of solitary waves[J]. Journal of Fluid Mechanics, 2014(7):577-596.

[4]WU N J, HSIAO S C. The study on solitary waves generated by a piston-type wave maker[J]. Ocean Engineering, 2016(11):114-129.

[5]FENTON J. A ninth-order solution for the solitary wave [J].Journal of Fluid Mechanics,1972,53(2):257-271.

[6]李邦華，鄭向遠(yuǎn)，李偉，等.波浪水槽中S to kes五階波的數(shù)值生成[J].武漢理工大學(xué)學(xué)報(bào)(交通科學(xué)與工程版)，2016，40(2)：238-244，250.

[7]張莉，郭海燕，李效民.南海內(nèi)孤立波作用下頂張力立管極值響應(yīng)研究[J].振動(dòng)與沖擊，2013(10)：858-865.

[8]關(guān)暉，魏崗，杜輝.內(nèi)孤立波與潛艇相互作用的水動(dòng)力學(xué)特性[J].解放軍理工大學(xué)學(xué)報(bào)(自然科學(xué)版)，2012(5)：478-486.

[9]徐鑫哲.內(nèi)波生成機(jī)理及二維內(nèi)波數(shù)值水槽模型研究[D].哈爾濱：哈爾濱工程大學(xué)，2012.

[10]王飛.內(nèi)孤立波作用下小尺度豎直圓柱體的水動(dòng)力特性研究[D].青島：中國(guó)海洋大學(xué)，2012.

A Study on Numerical Generation of Solitary Wave

YAO Jianxi1,2,3)WU Jingping1,2,3)
(KeyLaboratoryofHighPerformanceShipTechnologyofMinistryofEducation,Wuhan430063,China)1)

(HubeiKeyLaboratoryofInlandShippingTechnology,Wuhan430063,China)2)(SchoolofTransportation,WuhanUniversityofTechnology,Wuhan430063,China)3)

To realize the simulation of wave generations by a piston-type wave maker, a dynamic boundary condition is developed to simulate the motion of a paddle based on the open source platform OpenFOAM. Taking the solitary wave as an example, two types of motion which are derived from the Goring’s theory for generation of solitary wave and the Fenton’s ninth order solution of solitary wave are determined. Based on that, the generations of solitary wave and its propagation are studied. The numerical results show that the depression of the free surface behind the Goring-based solitary wave is more significant than that behind the Fenton-based wave and the Goring-based wave height decays more quickly.

solitary wave; OpenFOAM; piston-type wave maker; numerical wave generation; flume experiment

2017-05-22

*國(guó)家自然科學(xué)基金項(xiàng)目(51609188,51609187,51609186)、內(nèi)河航運(yùn)技術(shù)湖北省重點(diǎn)實(shí)驗(yàn)室基金項(xiàng)目(NHHY201502)、中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(2016IVB007,2017IVB006)資助

U675.91

10.3963/j.issn.2095-3844.2017.04.014

姚建喜(1985—)：男,博士,講師,主要研究領(lǐng)域?yàn)榇芭c海洋結(jié)構(gòu)物水動(dòng)力性能