田正林，孫昭晨，梁書(shū)秀

(大連理工大學(xué) 海岸及近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024)

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

田正林，孫昭晨，梁書(shū)秀

(大連理工大學(xué) 海岸及近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024)

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

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

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

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

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

1 數(shù)值模型

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

1.1控制方程

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

(1)

(2)

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

(3)

μ為分子運(yùn)動(dòng)粘性。

1.2造波方法

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

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

(4)

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

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

(5)

化簡(jiǎn)式(5)得下式

(6)

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

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

(7)

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

(8)

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

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

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

(9)

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

s(t)=2uΔx+vHs

(10)

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

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

(11)

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

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

(12)

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

1.3消浪層

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

圖3 實(shí)際波形與理論波形的對(duì)比Fig.3 Comparisons of present wave elevation and theory wave elevation

海綿層消波方程μ(X)

(13)

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

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

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

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

2.1實(shí)際波形與理論波形對(duì)比

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

2.2波浪沿程變化情況

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

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

2.3波高對(duì)目標(biāo)波浪的影響

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

2.4時(shí)間步長(zhǎng)的影響

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

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

3 結(jié)論

將質(zhì)量源區(qū)高度值考慮到造波源強(qiáng)方程之中，量化源區(qū)上方水體高度，克服了水深小于1.5 m時(shí)造波過(guò)程中源區(qū)上方水體高度難以確定的困難，得到了適合于淺水的質(zhì)量源造波方法。該方法可以在水深小于1.5 m的淺水中造出比較理想的波浪，更適合于水平網(wǎng)格步長(zhǎng)相比于源區(qū)高度值小10倍的情況。文中驗(yàn)證了波高對(duì)所造波浪的影響、時(shí)間步長(zhǎng)對(duì)淺水中所造波浪的影響及波浪沿程變化，得出結(jié)論：當(dāng)波高與水深的比值小于0.2時(shí)所造波浪波高與理論波高值的擬合較好；當(dāng)時(shí)間步長(zhǎng)小于0.005 s時(shí)所造波浪無(wú)明顯的上下漂浮現(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].水電能源科學(xué)，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ù)值模擬及其流場(chǎng)結(jié)構(gòu)分析[D].哈爾濱：哈爾濱工業(yè)大學(xué)，2013.

[6]李宏偉.造波理論與方法研究[D].哈爾濱：哈爾濱工程大學(xué)，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

國(guó)家自然科學(xué)基金資助項(xiàng)目(51279028)；國(guó)家海洋局公益性行業(yè)科研專(zhuān)項(xiàng)(201405025-1)

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

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