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

姚建喜 吳靜萍

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

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

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

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

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

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

0 引 言

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

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

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

1 孤立波推板造波理論

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

圖1 推板造波示意圖

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

(1)

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

η=Hs2

(2)

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

(3)

(4)

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

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

(5)

(6)

(7)

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

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

2 數(shù)值造波方法

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

(8)

N-S方程為

(9)

(10)

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

(1)

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

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

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

(12)

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

(13)

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

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

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

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

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

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

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

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

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

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

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

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

圖4 時間步長依賴性分析

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

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

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

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

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

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

圖6 比較計算結(jié)果和試驗結(jié)果

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

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

6 結(jié) 論

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

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

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

[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é)報(交通科學(xué)與工程版)，2016，40(2)：238-244，250.

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

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

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

[10]王飛.內(nèi)孤立波作用下小尺度豎直圓柱體的水動力特性研究[D].青島：中國海洋大學(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

*國家自然科學(xué)基金項目(51609188,51609187,51609186)、內(nèi)河航運技術(shù)湖北省重點實驗室基金項目(NHHY201502)、中央高校基本科研業(yè)務(wù)費專項資金(2016IVB007,2017IVB006)資助

U675.91

10.3963/j.issn.2095-3844.2017.04.014

姚建喜(1985—)：男,博士,講師,主要研究領(lǐng)域為船舶與海洋結(jié)構(gòu)物水動力性能