數(shù)值模擬推板式波浪水槽的關(guān)鍵參數(shù)取值探討

關(guān) 超 吳靜萍* 詹成勝 劉博宇

(武漢理工大學(xué)交通學(xué)院1) 武漢 430063) (中交機(jī)電工程局有限公司船舶及貿(mào)易事業(yè)部2) 北京 100088)

0 引 言

隨著計(jì)算機(jī)技術(shù)和計(jì)算流體動力學(xué)(computational fluid dynamics，CFD)的發(fā)展，從20世紀(jì)60年代開始至今，計(jì)算流體動力學(xué)取得了豐碩的成果，出現(xiàn)了眾多的CFD商業(yè)軟件，如PHOENICS，CFX，F(xiàn)LUENT，STAR CCM+，F(xiàn)LOW-3D，F(xiàn)IDIP等，使得數(shù)值模擬方法以其成本低、周期短、精度高等優(yōu)勢在流體力學(xué)學(xué)科的科學(xué)研究和工程項(xiàng)目中已經(jīng)得到認(rèn)可和被推廣使用[1-2].數(shù)值模擬方法在波浪運(yùn)動、波浪與結(jié)構(gòu)物相互作用問題中的應(yīng)用，因自由面捕捉技術(shù)、造波技術(shù)、消波技術(shù)、動網(wǎng)格技術(shù)等關(guān)鍵技術(shù)起步相對較晚，與物理波浪水槽并稱為數(shù)值波浪水槽，在海岸工程、港口與航道工程、船舶與海洋工程專業(yè)的波浪問題研究中起著重要作用.

數(shù)值波浪水槽造波有3種方法：仿物理(推板式或搖板式)造波法[3-4]、邊界造波法[5]和源造波法[6-7].對于水波問題，屬于非定常流動.推板式造波適用于淺水波及中等水深波，而搖板式造波則適用于深水波.數(shù)值模擬非定常的水波問題，計(jì)算相當(dāng)耗時(shí)，除了必須采用高性能計(jì)算機(jī)和并行技術(shù)之外，其網(wǎng)格尺度和時(shí)間步長對減少計(jì)算時(shí)間和保證計(jì)算精度起著關(guān)鍵作用.

網(wǎng)格尺度、時(shí)間步長的取值問題一直受到研究者重視.郭曉宇[8]對源造波的數(shù)值波浪水槽進(jìn)行驗(yàn)證時(shí)采用的波參數(shù)為波長L=2 m，波高H=0.05 m，水平方向網(wǎng)格尺度取L/20和L/40兩種，垂向波高范圍內(nèi)網(wǎng)格尺度取H/8，結(jié)果顯示兩者與解析波形均吻合良好；此外又通過一系列具有不同非線性程度的波浪輸入條件進(jìn)行驗(yàn)證，波長分別為L=30，60和90 m，波高分別為H=1，3和4 m，水平方向網(wǎng)格尺度為L/120，L/180，L/240，發(fā)現(xiàn)隨著波浪非線性程度的增大，單位波長需要較多網(wǎng)格節(jié)點(diǎn)才可以滿足足夠的精度.查晶晶等[9]用OpenFOAM實(shí)現(xiàn)數(shù)值波浪水槽造波，采用推板式造波法模擬波高H=0.02 m，周期T=1 s的線性波，所選網(wǎng)格尺度在30 m水池長度范圍內(nèi)的單元數(shù)為600，800，1 000，1 500，在垂向波高范圍內(nèi)的單元數(shù)為10，20，30，時(shí)間步長選取0.008，0.005，0.002，0.001 s進(jìn)行探討.通過與解析解比較，最終確定縱向網(wǎng)格、垂向網(wǎng)格和時(shí)間步長為1 500，20，0.002和800，10，0.002為最優(yōu)參數(shù)，并得出結(jié)論，計(jì)算結(jié)果的精度主要受波面區(qū)域矩形單元的長寬比影響，網(wǎng)格單元越細(xì)長，引入的數(shù)值黏性越大，由此造出的波與目標(biāo)波波幅相差太大.Finnegan等[10]采用CFX構(gòu)建二維搖板式數(shù)值波浪水槽，模擬線性深水波與有限水深波，討論網(wǎng)格依賴性，推薦了最小網(wǎng)格尺度計(jì)算式；推薦的時(shí)間步長為T/50.管陳等[11]采用FLUENT軟件實(shí)現(xiàn)搖板式造波，模擬波長L=2 m，波高分別為H=0.04，0.06和0.1 m的規(guī)則波，選取RNG雷諾時(shí)均模型，時(shí)間步長約為0.001 s，自由液面附近網(wǎng)格最小尺度為波高的1/10～1/20，水平方向網(wǎng)格尺度沿池長方向逐漸加大，并將數(shù)值結(jié)果與解析解進(jìn)行比較，發(fā)現(xiàn)同一波長下，波高越小，其誤差越小，模擬精度越高.

此外，雖然一些學(xué)者研究數(shù)值波浪水槽時(shí)未對相關(guān)數(shù)值造波技術(shù)參數(shù)進(jìn)行依賴性的討論，但在文獻(xiàn)中給出了其取值.嚴(yán)汝建等[12]使用FLUENT軟件基于搖板式造波，數(shù)值模擬波參數(shù)為L=8 m，H=0.2 m的規(guī)則波，水平方向和自由液面附近垂向網(wǎng)格尺度分別為L/20和H/10，并選取RNGk-ε湍流模型.李宏偉[13]使用FLUENT建立二維數(shù)值水槽，采用搖板式造波法模擬參數(shù)為L=8 m，H=0.4 m，T=2.267 8 s的規(guī)則波，水槽長度方向上網(wǎng)格尺度為L/100，自由液面附近進(jìn)行加密，時(shí)間步長為0.001 s，采用RNGk-ε湍流模型.黃華等[14]利用推板式造波法模擬二維二階Stokes波，波浪參數(shù)為L=5.783 2 m，H=0.129 m，T=2 s，給出的網(wǎng)格尺度水平和豎直方向分別約為L/111和H/21，時(shí)間步長為T/400.廉靜靜等[15]采用FLUENT軟件，通過搖板式造波模擬參數(shù)為L=8 m，T=2.284 8 s的規(guī)則波時(shí)，時(shí)間步長選為0.001 s，采用k-ε湍流模型.

本文以某流體力學(xué)實(shí)驗(yàn)室的小型推板式波浪水槽為原型，基于商業(yè)軟件FLUENT建立推板式造波的二維數(shù)值波浪水槽.選取波浪水槽有代表性的波浪參數(shù)造波.考慮水波理論多數(shù)基于勢流理論，而真實(shí)的流體具有黏性，所以假設(shè)本文模擬的波浪運(yùn)動為層流流動.首先討論了沿波長方向不同網(wǎng)格尺度、自由面附近沿波高方向上不同網(wǎng)格尺度以及不同時(shí)間步長對計(jì)算精度和計(jì)算耗時(shí)的影響，通過將數(shù)值結(jié)果與解析解的比較，綜合考慮波高誤差、單位距離波高相對衰減量和計(jì)算耗時(shí)，推薦綜合性能較優(yōu)的網(wǎng)格尺度、時(shí)間步長；然后將基于推薦的網(wǎng)格尺度、時(shí)間步長的數(shù)值結(jié)果與試驗(yàn)結(jié)果進(jìn)行比較，兩種吻合程度相當(dāng)令人滿意.

1 數(shù)值方法的基本控制方程

本文中二維波浪流動屬于黏性、不可壓縮和非定常類型，其波浪傳播方向?yàn)閤軸正向，波浪高度方向向上為z軸正向，則流體連續(xù)性方程為

(1)

動量方程為

(2)

(3)

式中：u為x方向速度；w為z方向速度；t為時(shí)間；ρ為流體密度；p為壓強(qiáng)；ν為運(yùn)動黏度；g為重力加速度.

采用VOF方法捕捉自由液面，控制方程為

(4)

式中：F為體積分?jǐn)?shù)的函數(shù)，F(xiàn)=1為單元體全為水，F(xiàn)=0為單元體全為空氣，0

2 數(shù)值計(jì)算方法

本文建立二維數(shù)值波浪水槽.水槽長為19 m(后4 m為消波段)，高為0.8 m，造波水深為0.5 m.生成目標(biāo)波形的參數(shù)為L=1.5 m，T=0.995 1 s，H=0.06 m(為波峰致波谷之間的高度差)，屬于二階Stokes波.

采用建模軟件Gambit進(jìn)行數(shù)值波浪水槽的模型建立及網(wǎng)格劃分，網(wǎng)格均為結(jié)構(gòu)化網(wǎng)格.水槽底部設(shè)置邊界層網(wǎng)格，自由液面附近進(jìn)行網(wǎng)格加密，加密總高度為波高的2倍.數(shù)值波浪水槽與物理波浪水槽一樣，需采要取一定的措施消除水槽末端池壁反射波的影響.數(shù)值消波技術(shù)主要采用人工阻尼區(qū)[16]、輻射邊界條件[17]和主動吸收[18]，或者在數(shù)值波浪水槽末端使用稀疏網(wǎng)格進(jìn)行數(shù)值耗散[19].本文在消波段采用漸變的稀疏網(wǎng)格耗散波能，從而消除水槽末端反射波的影響.

數(shù)值波浪水槽上方設(shè)為壓力出口邊界條件；兩端及底部設(shè)為壁面邊界條件，其中左邊模擬物理推板造波，賦予推板一個(gè)UDF函數(shù)，采用layer的動網(wǎng)格方法，使其按照給定的方式進(jìn)行運(yùn)動.

數(shù)值波浪水槽模型和邊界條件見圖1，網(wǎng)格劃分見圖2.

圖1 數(shù)值波浪水槽模型和邊界條件

圖2 數(shù)值波浪水槽網(wǎng)格劃分

3 數(shù)值波浪水槽模擬計(jì)算

3.1 時(shí)間步長變化

時(shí)間步長的選取會影響數(shù)值計(jì)算結(jié)果的準(zhǔn)確度和計(jì)算耗時(shí).本文僅考慮波浪傳播方向及自由液面附近網(wǎng)格尺度，首先取計(jì)算網(wǎng)格尺度Δx=L/70，Δz=H/10，考慮黏性，采用層流模型，取Δt=T/2 000，T/1 000，T/500三種時(shí)間步長進(jìn)行數(shù)值模擬.

圖3為t=19 s時(shí)沿池長方向的波形空間分布，z為波面抬高.表1為x=5 m與x=8 m處時(shí)間步長Δt變化時(shí)數(shù)值波高與解析波高比較的相對誤差，以及單位距離波高相對衰減量.

圖3 時(shí)間步長變化對數(shù)值模擬波浪傳播的影響

Δt網(wǎng)格單元長寬比H1/m(x=5 m)波高相對誤差|H1-H|H×100%H2/m(x=8 m)波高相對誤差|H2-H|H×100%單位距離波高相對衰減量/%T/2 0003.570.054 88.670.053 710.500.61T/1 0003.570.056 46.000.056 26.330.11T/5003.570.057 14.830.058 03.330.50

由表1可知，隨著時(shí)間步長的增大，數(shù)值結(jié)果與解析解的誤差反而更小.由圖3的波形及單位距離波高相對衰減量來看，當(dāng)Δt=T/1 000時(shí)，波形的穩(wěn)定性相對較好.綜合考慮計(jì)算精度和波形穩(wěn)定性，Δt=T/1 000是合適的時(shí)間步長選擇.

3.2 網(wǎng)格尺度變化

網(wǎng)格尺度會對計(jì)算結(jié)果的精度產(chǎn)生影響，因此需對計(jì)算網(wǎng)格進(jìn)行依賴性測試.首先討論波長方向上網(wǎng)格尺度對數(shù)值造波的影響.取自由面附近波浪高度方向上網(wǎng)格尺度為Δz=H/10不變，波浪傳播方向(水平方向)上網(wǎng)格尺度分別取Δx=L/50，L/60，L/70，L/80，L/90，L/100，L/110.

圖4為在t=19 s時(shí)波形的空間分布.表2為x=5 m與x=8 m處網(wǎng)格尺度Δx變化時(shí)數(shù)值波高與解析波高比較的相對誤差，以及單位距離波高相對衰減量.

圖4 水平網(wǎng)格尺度變化對數(shù)值模擬波浪傳播的影響

Δx網(wǎng)格單元長寬比H1/m(x=5 m)波高相對誤差|H1-H|H×100%H2/m(x=8 m)波高相對誤差|H2-H|H×100%單位距離波高相對衰減量/%L/505.000.054 49.330.053 910.170.28L/604.170.055 87.000.055 57.500.17L/703.570.056 46.000.056 26.330.11L/803.130.057 05.000.056 95.170.06L/902.780.057 24.670.057 14.830.05L/1002.500.057 44.330.057 34.500.06L/1102.270.057 54.170.057 44.330.05

由圖4與表2可知，當(dāng)波浪傳播方向上網(wǎng)格尺度Δx=L/80時(shí)，數(shù)值結(jié)果與解析解的誤差已比較小.

保持Δx=L/80不變，取波浪高度方向上、在自由面附近加密的網(wǎng)格尺度分別為Δz=H/10,H/20,H/30,H/40.圖5為在t=19 s時(shí)波形的空間分布.表3為x=5 m與x=8 m處網(wǎng)格尺度Δz變化時(shí)數(shù)值波高與解析波高比較的相對誤差，以及單位距離波高相對衰減量.

圖5 豎直網(wǎng)格尺度變化對數(shù)值模擬波浪傳播的影響

Δz網(wǎng)格單元長寬比H1/m(x=5 m)波高相對誤差|H1-H|H×100%H2/m(x=8 m)波高相對誤差|H2-H|H×100%單位距離波高相對衰減量/%H/103.130.05705.000.056 95.170.06H/206.250.056 85.330.056 46.000.22H/309.380.056 85.330.056 46.000.22H/4012.500.056 95.170.056 46.000.28

由圖5和表3可知，波浪高度方向上網(wǎng)格尺度為Δz=H/10時(shí)數(shù)值計(jì)算結(jié)果基本趨近解析解.一味減小Δz并不能增加計(jì)算精度，要考慮合適的網(wǎng)格單元長寬比.

綜合上述工作，本文波浪數(shù)值模擬的計(jì)算網(wǎng)格尺度為Δx=L/80，Δz=H/10已經(jīng)能夠達(dá)到較好的精度.

4 與物理試驗(yàn)結(jié)果比較

依據(jù)數(shù)值模擬推板式波浪水槽波浪傳播的關(guān)鍵參數(shù)取值分析，所推薦的網(wǎng)格尺度Δx=L/80，Δz=H/10和時(shí)間步長Δt=T/1 000，選取層流模型，數(shù)值模擬波長為1.5 m、波高為0.06 m的波，將此數(shù)值計(jì)算結(jié)果與物理試驗(yàn)比較，進(jìn)一步驗(yàn)證該數(shù)值算法的可靠性.物理試驗(yàn)水槽為武漢理工大學(xué)流體力學(xué)實(shí)驗(yàn)室小型推板式波浪水槽，長為18 m、寬為0.6 m、高為0.8 m、水深為0.5 m.試驗(yàn)采用伺服電機(jī)控制推板造波，在水槽x=5 m與x=10 m處各布置一個(gè)浪高儀，用以測量波高隨時(shí)間的變化.

在數(shù)值計(jì)算中，同樣在x=5 m與10 m處設(shè)置虛擬浪高儀，得出隨時(shí)間變化的波形圖.

圖6為數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果的比較.表4為x=5 m與10 m處數(shù)值平均波高與試驗(yàn)平均波高的比較及相對誤差(圖6方框內(nèi)平均波高).

圖6 數(shù)值結(jié)果與試驗(yàn)結(jié)果的比較

H1/m(x=5 m)H2/m(x=10 m)試驗(yàn)結(jié)果 0.057 20.055 9數(shù)值結(jié)果 0.056 50.055 3相對誤差/%1.221.07

由圖6與表4可知，數(shù)值模擬的波形與試驗(yàn)波形吻合程度很好，其變化趨勢表現(xiàn)出良好的一致性，驗(yàn)證了該數(shù)值算法的有效性.

5 結(jié) 論

1) 本文數(shù)值波浪水槽中推薦采用：波浪傳播方向網(wǎng)格尺度取Δx=L/80，自由液面附近波高方向網(wǎng)格加密尺度取Δz=H/10，時(shí)間步長取Δt=T/1 000；網(wǎng)格尺度確定還要滿足一定的單元長寬比.

2) 數(shù)值模擬的結(jié)果與物理試驗(yàn)結(jié)果吻合程度很好，驗(yàn)證了本文數(shù)值算法的有效性.

在展開新型浮式防波堤衰減長周期波性能研究時(shí)遭遇物理波浪水槽長度尺度限制的問題.本文工作為發(fā)揮數(shù)值波浪水槽的優(yōu)勢，利用數(shù)值模擬手段輔助物理試驗(yàn)進(jìn)行研究奠定了基礎(chǔ).也為讀者應(yīng)用數(shù)值波浪水槽進(jìn)行研究提供有益參考.

參考文獻(xiàn)

[1] 王福軍.計(jì)算流體動力學(xué)分析:CFD軟件原理與應(yīng)用[M].北京:清華大學(xué)出版社，2004.

[2] 韓占忠.FLUENT：流體工程仿真計(jì)算實(shí)例與應(yīng)用[M].北京:北京理工大學(xué)出版社，2010.

[4] PRASAD D D, AHMED M R, LEE Y H, et al. Validation of a piston type wave-maker using Numerical Wave Tank [J]. Ocean Engineering, 2017,131:57-67.

[5] 劉霞，譚國煥，王大國.基于邊界造波法的二階Stokes波的數(shù)值生成[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào)(自然科學(xué)版)，2010，29(1):107-111.

[6] LIN P, LIU P, L F. Internal wave-maker for navier-stokes equations models[J]. Journal of Waterway, Port, Coastal and Ocean Engineering, 1999,125(4):207-217.

[7] WEI G, KIRBY J T, SINHA A. Generation of waves in boussinesq models using a source function method[J]. Coastal Engineering, 1999,36(4):271-299.

[8] 郭曉宇.數(shù)值波浪水槽及其應(yīng)用研究[D].上海：上海交通大學(xué)，2011.

[9] 查晶晶，萬德成.用OpenFOAM實(shí)現(xiàn)數(shù)值水池造波和消波[J].海洋工程,2011,29(3):1-12.

[10] FINNEGAN W, GOGGINS J. Numerical simulation of linear water waves and wave-structure interaction [J]. Ocean Engineering, 2012,43(4):23-31.

[11] 管陳，董國祥，金允龍.三維數(shù)值波浪水池造波技術(shù)研究[J].上海船舶運(yùn)輸科學(xué)研究所學(xué)報(bào),2013,36(2):11-15.

[12] 嚴(yán)汝建，龐永杰，李宏偉，等.深水池造波系統(tǒng)數(shù)值造波仿真研究[J].哈爾濱工程大學(xué)學(xué)報(bào),2010,31(1):32-41.

[13] 李宏偉.數(shù)值水池造波方法研究[D].哈爾濱：哈爾濱工程大學(xué)，2011.

[14] 黃華，鄧冰，陳昱松，等.數(shù)值波浪水槽構(gòu)建與二階Stokes波仿真[J].系統(tǒng)仿真學(xué)報(bào),2012,24(1):227-231.

[15] 廉靜靜，尹勇，楊曉，等.基于黏性流船舶數(shù)值波浪水池造波和消波方法研究[J].船舶力學(xué),2013(1):56-62.

[16] 趙西增，孫昭晨，梁書秀.非線性波浪水槽內(nèi)波浪的產(chǎn)生和傳播[J].船舶力學(xué),2010,14(3):203-216.

[17] ORLANSKI I. A simple boundary condition for unbounded hyperbolic flows[J]. J Comput Phys, 1976,21(3):251-269.

[18] DEEPAK D P,MOHAMMED R A,YOUNG H L,et al. Validation of a piston type wave-maker using Numerical Wave Tank [J]. Ocean Engineering, 2017,131:57-67.

[19] PARK J C, UNO Y, SATO T, et al. Numerical reproduction of fully nonlinear multi-directional waves by a viscous 3D numerical wave tank[J]. Ocean Engineering, 2004,31(11):1549-1565.