田正林，孫昭晨，梁書秀

(大連理工大學 海岸及近海工程國家重點實驗室，大連 116024)

淺水中質(zhì)量源造波方法

田正林，孫昭晨，梁書秀

(大連理工大學 海岸及近海工程國家重點實驗室，大連 116024)

針對質(zhì)量源造波過程中源區(qū)上方水體高度難以確定的問題提出了改進方法，將源區(qū)高度值考慮到源強方程中，避免了造波過程中源區(qū)上方水體高度值多次試算帶來的不便。該方法彌補了質(zhì)量源造波方法在水深小于1.5 m的淺水中所造波浪與理論波浪之間差距較大的缺陷，尤其是水平網(wǎng)格步長與源區(qū)高度值相比小于10的情況，推廣了質(zhì)量源造波方法在淺水中的應(yīng)用。文中驗證了波高對所造波浪的影響以及時間步長對淺水中所造波浪的影響，得出結(jié)論：當波高與水深的比值小于0.2時所造波浪波高與理論值的擬合較好；當時間步長小于0.005 s時所造波浪無明顯的上下漂浮現(xiàn)象，波形穩(wěn)定，符合理論波形要求。

質(zhì)量源造波；淺水造波；源區(qū)高度；改進方法

數(shù)值波浪水槽是港口海岸工程數(shù)值模擬的基礎(chǔ)，造波系統(tǒng)與消波系統(tǒng)是數(shù)值波浪水槽的重要組成部分，造波系統(tǒng)造波能力及消波系統(tǒng)消波能力又是數(shù)值波浪水槽可用性的關(guān)鍵。大多數(shù)值模擬采用加長計算域或是波動沒有達到準穩(wěn)定狀態(tài)時結(jié)束模擬避免二次反射波的影響[1]。為了得到長時間穩(wěn)定的數(shù)值模擬，對數(shù)值水槽的要求是一能夠降低二次反射波，二所造目標波浪不受反射波的影響，流體域內(nèi)部造波方法能夠很好地解決該問題。Larsen和Dancy[2]在水域內(nèi)部沿一定長度的直線段上通過擾動水面生成目標波浪，避免了二次反射波的影響。

目前，流體內(nèi)域造波主要有動量源造波和質(zhì)量源造波。動量源造波是將目標波浪作為人工源項加到動量方程中生成所需波浪場，質(zhì)量方程保持不變，自由水面以下區(qū)域均屬造波源區(qū)，造波過程中在整個造波源區(qū)內(nèi)施加動量源方程；質(zhì)量源造波是將目標波浪作為人工源項加到質(zhì)量方程中生成所需波浪場，造波過程中并不是自由水面以下區(qū)域均屬造波源區(qū)，源區(qū)頂部與自由水面之間存在一定高度無質(zhì)量源的水體，正因為此部分水體的存在使得源區(qū)高度的確定變得繁瑣。Lin和Liu[3]通過試算得出質(zhì)量源造波源區(qū)中心到靜水面之間的距離為13～12水深時所造波浪比較理想。周玲玲，丁全林，蘭慶琳等[4]指出質(zhì)量源距離水面13水深時的模擬結(jié)果較為精確，并將此處定義為質(zhì)量源的標準位置。余勇飛[5]在深水數(shù)值水槽造波時質(zhì)量源區(qū)上方水體高度采用一倍波高時所造波浪精度較高，與理論波浪匹配較好。李宏偉[6]在造孤立波時質(zhì)量源區(qū)上方水體高度采用半倍波高時，所造孤立波波形與理論波形匹配較好。Chen和Hsiao[7]在FLOW-3D中開發(fā)了水體內(nèi)域造波質(zhì)量源區(qū)域的設(shè)計方法，詳細描述了質(zhì)量源位置的確定公式，為質(zhì)量源造波方法的應(yīng)用提供了便利。

已有文獻質(zhì)量源區(qū)上方水體高度大多采用試算法確定，只給出經(jīng)驗高度區(qū)間，并且試算法的成功快慢取決于個人經(jīng)驗，對于經(jīng)驗不足者尤顯耗時。筆者造波過程中發(fā)現(xiàn)質(zhì)量源區(qū)上方水體高度為一倍波高時在水深大于1.5 m的深水情況下所造波浪極好，但在水深小于1.5 m的淺水中所造波浪并不理想，尤其是在水平網(wǎng)格步長與源區(qū)高度值相比小于10時所造波浪波高與目標波浪波高相差較大，必須不斷地加大或減小源區(qū)上方水體高度以減小差距。而通常情況下需要造多種波高的目標波浪，這樣每次變化波高時都需要不斷調(diào)整源區(qū)上方水體高度，既費時也是沒有必要的。故本文針對質(zhì)量源造波源區(qū)上方水體高度難以確定的問題，提出改進方法，避免了多次試算帶來的麻煩，進一步完善了質(zhì)量源造波方法。

1 數(shù)值模型

基于FLUENT求解器，求解雷諾平均N-S方程，采用雷諾應(yīng)力粘性模型封閉方程組。采用有限體積法離散方程，PISO壓力速度耦合算法求解離散方程，對流項采用三階MUSCL格式，用改進的QUICK方法跟蹤自由表面。邊界條件為左右側(cè)為對稱邊界條件，底部為墻邊界條件，上部為壓力進口邊界條件。以下僅就文中的控制方程、造波方法及消波方法進行描述。

1.1控制方程

張量形式的連續(xù)性方程和動量方程如下

(1)

(2)

對于二維問題i,j=1,2；ui,uj為速度；p為壓強；ρ為密度；g為重力加速度。粘性應(yīng)力張量τij為

(3)

μ為分子運動粘性。

1.2造波方法

數(shù)值波浪水槽的質(zhì)量源設(shè)置于計算域內(nèi)，質(zhì)量源區(qū)距離左側(cè)消波區(qū)至少一倍波長的距離可以很好地消除波浪二次反射(Lin & Liu[1])，本文造波源區(qū)距離左側(cè)消波區(qū)2L(L為波長)。

造波源區(qū)內(nèi)質(zhì)量守恒方程為

(4)

s(xi,t)為質(zhì)量源強方程，其大小由設(shè)計的入射波波面方程決定，造波源區(qū)內(nèi)非零，造波源區(qū)外為零。理論上，當質(zhì)量源寬度遠小于波長時，可作點源處理，假定質(zhì)量源引起的質(zhì)量增減對應(yīng)于入射波波面的變化，并且源區(qū)內(nèi)質(zhì)點的水平速度等于設(shè)計波浪的水平速度，因為質(zhì)量源可以看做均勻分布的點源，所以在造波源域Ω內(nèi)每個網(wǎng)格中應(yīng)滿足下列關(guān)系[5-6，8-9]

s(xi,t)ΔxΔyΔt=2u(xi,t)ΔyΔt

(5)

化簡式(5)得下式

(6)

式中：xi為“源強”的水平位置；Δx為網(wǎng)格的水平尺度；Δy為網(wǎng)格的垂向尺度；右端系數(shù)2表示質(zhì)量源產(chǎn)生傳播方向相反的兩列波。

二維二階Stocks 波的水平速度方程如下

(7)

將式(7)代入(6)得源強方程如下

(8)

式中：H為波高；k為波數(shù)；ω為頻率；D為水深。

以下通過對造波源區(qū)域源強方程的改進，可以得到更適合淺水區(qū)域的造波方法。造波源區(qū)簡圖如圖1所示，造波源區(qū)底部為墻邊界并無水體向下噴出，考慮單位時間內(nèi)造波源區(qū)其他3個方向噴出的水體質(zhì)量守恒可推得

s(t)·Δx·Hs=2u·Hs+Δx·v

(9)

式中：u為造波源區(qū)流體質(zhì)點水平速度；v為造波源區(qū)流體質(zhì)點垂向速度；Hs為造波源區(qū)高度?；?9)可得

s(t)=2uΔx+vHs

(10)

圖1 源造波區(qū)簡圖Fig.1 Sketch of the source wave-maker zone

質(zhì)量源造波假定源區(qū)內(nèi)各點產(chǎn)生的水平速度等于波生成對應(yīng)點產(chǎn)生的水平速度，造波源區(qū)上方水體依靠水體自重形成的恢復(fù)力產(chǎn)生波浪，因此源區(qū)上方至自由水面高度的確定比較關(guān)鍵。仔細研究方程(10)可以發(fā)現(xiàn)：造波過程中u和v的量級相當，由于Δx一般小于波長的5%，非常小，而Hs一般大于23倍的水深，對于水深大于1.5 m的深水Δx相比Hs更為小量，至少小一個量級，式中的第二項對源強的貢獻較小可以忽略，但對于水深小于1.5 m的淺水，第二項會隨著水深的減小而增大，其對造波源的貢獻也會逐漸增大，如果第二項繼續(xù)忽略會影響所造波浪的波高。數(shù)值波浪水槽造波時如果水深小于1.5 m，Hs小于1.0 m，Δx又小于0.1 m，可以看出此時第一項與第二項相比在10以內(nèi)，不足一個量級，如果水深繼續(xù)減小，Hs繼續(xù)減小，式中第二項對源強的貢獻會更大，不應(yīng)忽略。根據(jù)Broren.M.,Larsen.J.[10]描述，源造波的原理類似于“水下爆炸”原理，那么源爆炸瞬間質(zhì)點各個方向的速度大小應(yīng)該相等。所以質(zhì)量源向外噴出的流體速度各方向應(yīng)該大小相等，因此u=v，最終可得源強方程為

(11)

周玲玲，丁全林，蘭慶琳等[4]指出為避免數(shù)值振蕩，質(zhì)量源初始時刻應(yīng)設(shè)置為0，二階Stocks波令初始相位為-π2，并在計算歷程開始時刻設(shè)置軟啟動緩慢施加質(zhì)量源驅(qū)動波面變化。本文采用此種做法，最終施加的源強方程如下

圖2 源區(qū)右側(cè)數(shù)值水槽示意圖Fig.2 Right side schematic diagram of numerical wave tank

(12)

式中：N為緩變控制參數(shù)；t0為緩變時間。

1.3消浪層

為了吸收二次反射波和不需要的源造波，數(shù)值波浪水槽左右端應(yīng)布置至少一倍波長的海綿吸波層耗散波能。本文采用的吸收方程是由Niels and David[11]提出的橢圓方程(13)，x0表示消波層起始點x坐標；xs表示消波層長度，方程從海綿層起點1逐漸減小到海綿層終點0。

圖3 實際波形與理論波形的對比Fig.3 Comparisons of present wave elevation and theory wave elevation

海綿層消波方程μ(X)

(13)

數(shù)值波浪水槽源區(qū)右半側(cè)示意圖，如圖2所示。

2 數(shù)值模擬結(jié)果

本文數(shù)值波浪水槽長36 m，高0.7 m，水深0.5 m，質(zhì)量源高度Hs為水深的23，△x小于波長的5%。以下對所造波浪結(jié)果、波形沿程變化、波高對目標波浪的影響以及時間步長對波形的影響進行描述。

圖4 波形沿程變化圖Fig.4 Wave elevation displacement of present model with distance from the source zone

2.1實際波形與理論波形對比

為了驗證所造波浪的有效性，給出水深D=0.5 m，波高H=0.06 m，周期T分別為0.8 s、1.0 s、1.2 s、1.4 s時的實造波浪波形與理論波浪波形的對比(圖3)。從圖3中可以看出：不同周期時的實際波形與理論波浪波形匹配較好，雖波谷處振幅值比理論值略小但仍都可以滿足精度要求。

2.2波浪沿程變化情況

為了說明波浪沿程變化情況，給出了水深D=0.5 m，波高H=0.06 m，周期T分別為0.8 s、1.0 s、1.2 s、1.4 s時造波源區(qū)(約5倍波長處)右側(cè)波浪10倍波長范圍內(nèi)的波浪沿程變化情況(圖4)。從圖4可以看出：所造波浪沿程衰減較小，滿足目標波浪精度要求。

圖5 所造波浪波高與理論波高的對比Fig.5 Comparisons of the targeted wave height and theoretical wave height

2.3波高對目標波浪的影響

為了檢驗改進的質(zhì)量源造波方法所造不同波高的波浪與理論波浪的差距，本部分計算了D=0.5 m，T=1.2 s，波高H分別為0.06 m、0.08 m、0.10 m、0.12 m造波源區(qū)(約5倍波長處)右側(cè)波浪10倍波長范圍內(nèi)的波浪隨時間的變化(圖5)。從圖5可以看出：隨著波高的增大，所造波浪波高與理論波高的差距在增大；波高H=0.10 m時所造波浪波高只是在波峰處略小于理論值，波谷處仍與理論值擬合較好，但是當波高H=0.12 m時所造波浪波峰和波谷處均與理論值有偏差，所造波浪效果變差。由此可知當波高與水深的比值超過0.2時所造波浪波高與理論值的擬合較差，所以此種造波方法適合于波高與水深的比值在0.2以下的情況。

2.4時間步長的影響

圖6 波形隨時間變化圖Fig.6 Time step effects on wave surface

為了檢驗時間步長對改進的質(zhì)量源造波方法所造波浪的影響，給出了水深D=0.5 m，周期T=1.2 s，波高H=0.06 m時所造波浪波形隨時間的變化。時間步長驗證之前需要對網(wǎng)格進行反復(fù)調(diào)整驗證，當網(wǎng)格連續(xù)3次調(diào)整后的波形無較大變化后停止，進行時間步長驗證。時間步長分別取0.01 s、0.005 s、0.002 s、0.000 2 s進行驗證。對比所得波面由圖6可以看出：當dt=0.01 s、dt=0.005 s時，隨著時間的增長波面出現(xiàn)了幅度較大的上下飄動現(xiàn)象，而當dt=0.002 s、dt=0.000 2 s時，波面已無明顯的上下飄動現(xiàn)象，波面較穩(wěn)定。

3 結(jié)論

將質(zhì)量源區(qū)高度值考慮到造波源強方程之中，量化源區(qū)上方水體高度，克服了水深小于1.5 m時造波過程中源區(qū)上方水體高度難以確定的困難，得到了適合于淺水的質(zhì)量源造波方法。該方法可以在水深小于1.5 m的淺水中造出比較理想的波浪，更適合于水平網(wǎng)格步長相比于源區(qū)高度值小10倍的情況。文中驗證了波高對所造波浪的影響、時間步長對淺水中所造波浪的影響及波浪沿程變化，得出結(jié)論：當波高與水深的比值小于0.2時所造波浪波高與理論波高值的擬合較好；當時間步長小于0.005 s時所造波浪無明顯的上下漂浮現(xiàn)象，波形穩(wěn)定，符合理論波形要求，波浪沿程衰減較小。

[1]Lin P, Liu L F. Discussion of "Vertical variation of the flow across the surf zone" [Coast. Eng. 45 (2002) 169-198][J]. Coastal Engineering, 2004, 50(3):161-164.

[2]Larsen J, Dancy H. Open boundaries in short wave simulations-A new approach[J]. Coastal Engineering, 1983, 7(3):285-297.

[3]Lin P, Liu L F. Internal Wave-maker for Navier-stokes Equations Models[J].Journal of Waterway, Port, Coastal and Ocean Engineering, 1999,125(4): 207-215.

[4]周玲玲,丁全林,蘭慶琳,等.基于LB方法的質(zhì)量源造波的數(shù)值波浪水槽[J].水電能源科學，2015,33(3):99-103. ZHOU L L,DING Q L,LAN Q L, et al.Numerical Wave Tank with Mass Source Generating Waves Based on LB Method[J]. China Civil Engineering Journal, 2015, 33(3): 99-103.

[5]余勇飛.波浪的數(shù)值模擬及其流場結(jié)構(gòu)分析[D].哈爾濱：哈爾濱工業(yè)大學，2013.

[6]李宏偉.造波理論與方法研究[D].哈爾濱：哈爾濱工程大學，2013.

[7]Chen Y L, Hsiao S C. Generation of 3D water waves using mass source wavemaker applied to Navier-Stokes model[J]. Coastal Engineering, 2016, 109:76-95.

[8]房卓,張寧川,臧志鵬.透空式梳式防波堤的數(shù)值模擬和波浪透射系數(shù)的研究[J].水道港口,2011,32(2)：86-93. FANG Z, ZHANG N C, ZANG Z P. Numerical simulations of open comb-type breakwater and research on its wave transmission coefficient[J].Journal of Waterway and Harbor,2011,32(2):86-93.

[9]FANG Z, ZHANG N C, ZANG Z P. Experimental and Numerical Study on Hydrodynamic Performance of Impermeable Comb-Type Breakwater[J]. Journal of Ship Mechanics, 2012, 16(6):632-645.

[10]Brorsen M, Larsen J. Source generation of nonlinear gravity waves with the boundary integral equation method[J]. Coastal Engineering, 1987, 11(2):93-113.

[11]Jacobsen N G, Fuhrman D R, Fredsφe J. A wave generation toolbox for the open-source CFD library: OpenFoam[J]. International Journal for Numerical Methods in Fluids, 2012, 70(9):1 073-1 088.

The method of mass source wavemaker in shallow water

TIANZheng-lin,SUNZhao-chen,LIANGShu-xiu

(StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian116024,China)

In view of the problem that the distance between the top of the source region and the still water level is difficult to determine in the process of mass source wavemaker, the improved method was proposed in this paper. In the mass source functions took the height of the source region into account, the inconvenience of many trials to get the height of the mass source region can be eliminated in the process of wavemaker. This method remedies the defect that the difference between targeted wave and theoretical wave in the shallow water whose depth less than 1.5 m, especially horizontal grid step is an order of magnitude 10 smaller than the source region height, the mass source wave-maker method is promoted in the shallow water.The effects of wave height and the time step on simulated waves were studied. The conclusions show that the fitting between made wave and theoretical wave is good when the ratio of wave height and water depth is less than 0.2;when the time step is less than 0.005 s, the floating phenomenon of the simulated waves is well controlled, waves surface is stable and conforms to the requirements of the wave theory.

mass source wavemaker; shallow-water wavemaker; height of source region; improved method

TV 139.2

：A

：1005-8443(2017)04-0325-05

2017-01-13；

：2017-04-06

國家自然科學基金資助項目(51279028)；國家海洋局公益性行業(yè)科研專項(201405025-1)

田正林(1985-),男，黑龍江省哈爾濱人，博士研究生，主要從事波浪與結(jié)構(gòu)物相互作用的研究。

Biography： TIAN Zheng-lin(1985-),male, doctor student.