關 超 吳靜萍* 詹成勝 劉博宇

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

0 引 言

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

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

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

此外，雖然一些學者研究數(shù)值波浪水槽時未對相關數(shù)值造波技術參數(shù)進行依賴性的討論，但在文獻中給出了其取值.嚴汝建等[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，自由液面附近進行加密，時間步長為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，時間步長為T/400.廉靜靜等[15]采用FLUENT軟件，通過搖板式造波模擬參數(shù)為L=8 m，T=2.284 8 s的規(guī)則波時，時間步長選為0.001 s，采用k-ε湍流模型.

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

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

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

(1)

動量方程為

(2)

(3)

式中：u為x方向速度；w為z方向速度；t為時間；ρ為流體密度；p為壓強；ν為運動黏度；g為重力加速度.

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

(4)

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

2 數(shù)值計算方法

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

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

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

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

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

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

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

3.1 時間步長變化

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

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

圖3 時間步長變化對數(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ù)值結(jié)果與解析解的誤差反而更小.由圖3的波形及單位距離波高相對衰減量來看，當Δt=T/1 000時，波形的穩(wěn)定性相對較好.綜合考慮計算精度和波形穩(wěn)定性，Δt=T/1 000是合適的時間步長選擇.

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

網(wǎng)格尺度會對計算結(jié)果的精度產(chǎn)生影響，因此需對計算網(wǎng)格進行依賴性測試.首先討論波長方向上網(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時波形的空間分布.表2為x=5 m與x=8 m處網(wǎng)格尺度Δx變化時數(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可知，當波浪傳播方向上網(wǎng)格尺度Δx=L/80時，數(shù)值結(jié)果與解析解的誤差已比較小.

保持Δx=L/80不變，取波浪高度方向上、在自由面附近加密的網(wǎng)格尺度分別為Δz=H/10,H/20,H/30,H/40.圖5為在t=19 s時波形的空間分布.表3為x=5 m與x=8 m處網(wǎng)格尺度Δz變化時數(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ù)值計算結(jié)果基本趨近解析解.一味減小Δz并不能增加計算精度，要考慮合適的網(wǎng)格單元長寬比.

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

參考文獻

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

[2] 韓占忠.FLUENT：流體工程仿真計算實例與應用[M].北京:北京理工大學出版社，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].遼寧工程技術大學學報(自然科學版)，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ù)值波浪水槽及其應用研究[D].上海：上海交通大學，2011.

[9] 查晶晶，萬德成.用OpenFOAM實現(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ù)值波浪水池造波技術研究[J].上海船舶運輸科學研究所學報,2013,36(2):11-15.

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

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

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

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

[16] 趙西增，孫昭晨，梁書秀.非線性波浪水槽內(nèi)波浪的產(chǎn)生和傳播[J].船舶力學,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.