，　

(1.寧波市水利水電規(guī)劃設(shè)計(jì)研究院，浙江　寧波　315192； 2.中交第二航務(wù)工程勘察設(shè)計(jì)院有限公司，湖北　武漢　430071)

1　問(wèn)題的提出

波浪在近岸傳播時(shí)，由于受到地形、海底摩阻或固定式結(jié)構(gòu)物的影響,將會(huì)發(fā)生折射、繞射、反射、能量衰減或破碎等現(xiàn)象。關(guān)于固定式結(jié)構(gòu)物周?chē)牟ɡ藗鞑プ冃蔚难芯?早期的工作主要訴諸于解析法，其研究的水域水深定常，研究的結(jié)構(gòu)物的形狀簡(jiǎn)單、規(guī)則，如等水深水域內(nèi)無(wú)窮小厚度的島式防波堤、半無(wú)限防波堤、雙突堤等。

隨著波浪傳播數(shù)學(xué)模型的逐步成熟以及計(jì)算技術(shù)的迅速發(fā)展，數(shù)值模擬方法已成為研究波浪傳播問(wèn)題的主要方法之一。目前，較為活躍的數(shù)學(xué)模型有Boussinessq方程[1]和緩坡方程[2]。與Boussinessq方程相比，緩坡方程具有形式簡(jiǎn)單、數(shù)值求解方法多樣、數(shù)值計(jì)算量小的優(yōu)點(diǎn)，因此，緩坡方程被國(guó)內(nèi)外工程界廣泛應(yīng)用于研究近岸波浪的傳播變形問(wèn)題。

緩坡方程是由Berkholff基于線性勢(shì)流理論，利用小參數(shù)攝動(dòng)法推導(dǎo)得到的能夠考慮波浪聯(lián)合折射繞射影響的方程[2]，該方程具有完全頻散性，并且可以方便地?cái)U(kuò)展到含水流、弱非線性、波浪破碎和底摩阻引起的能量衰減和快速的地形變化的情況。自問(wèn)世至今，該方程在方程形式上和計(jì)算方法上，已得到了很大的發(fā)展。在方程形式上，有橢圓型緩坡方程[2]、時(shí)間關(guān)聯(lián)型緩坡方程[3]、拋物形緩坡方程[4]、一階雙曲型緩坡方程組[5-6]、等價(jià)的雙曲型控制方程組[7]、水波演化方程[8]。緩坡方程的各種形式及其相應(yīng)的數(shù)值求解方法各有其優(yōu)點(diǎn)和不足，對(duì)此，文獻(xiàn)[9-10]進(jìn)行了詳細(xì)的概括。各種緩坡方程的近似形式或等價(jià)形式，都是為了克服直接求解原型緩坡方程計(jì)算量大的困難，但隨著計(jì)算機(jī)軟硬件的發(fā)展，計(jì)算機(jī)的計(jì)算速度已不再像過(guò)去那樣嚴(yán)重地制約著數(shù)值計(jì)算，部分學(xué)者又對(duì)原型緩坡方程直接數(shù)值求解[11-12]。本文工作亦基于原型緩坡方程，未做任何近似和假定。

針對(duì)計(jì)算域中存在直立島式結(jié)構(gòu)物的復(fù)連通區(qū)域，基于時(shí)間關(guān)聯(lián)型緩坡方程和相應(yīng)的邊界條件，建立了在計(jì)算域中存在直立島式結(jié)構(gòu)物時(shí)波浪傳播的數(shù)值模擬模型。模型在方程形式上沒(méi)有任何近似，并且由于是對(duì)時(shí)間關(guān)聯(lián)型方程數(shù)值求解，所以模型適用于模擬波浪的時(shí)間和空間演化過(guò)程。在模型驗(yàn)證上，對(duì)直立島式防波堤及直立方柱周?chē)ɡ藗鞑プ冃芜M(jìn)行了數(shù)值模擬。結(jié)果表明，本文所建立的模型能夠有效地模擬計(jì)算域內(nèi)存在直立島式結(jié)構(gòu)物的波浪繞射和反射問(wèn)題。

2　控制方程和數(shù)值離散格式

2.1　控制方程

基于變分原理，Kirby[3](1993年)推導(dǎo)了時(shí)間關(guān)聯(lián)型緩坡方程，該方程為雙曲型偏微分方程組，其形式如下：

(1)

(2)

式中：g是重力加速度(m/s2)；η是復(fù)波面函數(shù)；Φ是復(fù)速度勢(shì)函數(shù)；h是局部水深(m)，κ是波數(shù)；σ是波浪角頻率(rad/s)；c、cg分別是波速和波群速(m/s)，其定義如下：

(3)

波數(shù)k和頻率σ滿(mǎn)足線性彌散關(guān)系，其表達(dá)式如下：

σ2=gktanh(kh)

(4)

方程(1)和(2)可以改寫(xiě)成如下形式：

(5)

(6)

(7)

F(η，Φ)=-gη

(8)

2.2　數(shù)值離散格式

方程(5)和(6)為時(shí)間關(guān)聯(lián)型方程，對(duì)該方程進(jìn)行數(shù)值求解，可在時(shí)間上采用Euler預(yù)測(cè)—較正—迭代格式，空間上采用三點(diǎn)差分格式[12]。該方法在時(shí)間和空間上都具有二階精度，其時(shí)間上的離散格式如下：

Euler預(yù)測(cè)步：

(9)

(10)

式中：

(11)

(12)

Euler校正步：

(13)

(14)

式中：

(15)

(16)

3　初邊值條件

一般來(lái)說(shuō)，邊界條件有3類(lèi)，分別是：

(1)入射邊界條件。在入射邊界處，允許反射波自由地流出入射邊界，其邊界條件可采用如下:

(18)

(19)

式中：n是邊界的外法向；η0和Φ0分別是入射邊界處的入射波復(fù)波面和復(fù)速度勢(shì)，其表達(dá)式如下：

(20)

(21)

式中：H0是入射波高(m)；α0是入射波向角(°)。

(2)對(duì)于直立式結(jié)構(gòu)物的邊界，可處理為固壁邊界條件

(22)

(3)輻射邊界條件

(23)

(24)

4　模型驗(yàn)證

4.1　單突堤附近的波浪傳播變形

單突堤附近的波高分布可按Sommerfeld’s繞射理論給出解析解(見(jiàn)圖1)。為了檢驗(yàn)本文數(shù)值模擬模型關(guān)于波浪繞射現(xiàn)象的計(jì)算效果，首先對(duì)單突堤附近波浪的傳播變形進(jìn)行了數(shù)值計(jì)算。

考慮線性單頻波正向入射，波高H為0.5 m，波周期T為6.0 s，數(shù)值計(jì)算的空間范圍為13L×13L，長(zhǎng)度為8L的單突堤位于y=0；空間步長(zhǎng)和時(shí)間步長(zhǎng)分別為：Δx=Δy=L/16，Δt=T/32，其中L為入射波波長(zhǎng)。模型下游邊界為輻射邊界條件，側(cè)邊界和單突堤按固壁全反射邊界條件處理。模型計(jì)算的數(shù)值波浪的時(shí)間為33T，波高由t=30T時(shí)刻到t=33T時(shí)刻的波面值按照均方根波高進(jìn)行計(jì)算[13]。圖2是數(shù)值計(jì)算的單突堤附近的相對(duì)波高分布圖。比較圖1和圖2可知，波浪在突堤前發(fā)生的反射現(xiàn)象和在突堤后發(fā)生的繞射現(xiàn)象，均可較為合理地由本文所建立的數(shù)值模擬模型給出。

圖1Sommerfeld’s繞射理論給出的單突堤附近相對(duì)波高的分布圖[12]單位：m

圖2　數(shù)值計(jì)算的單突堤附近相對(duì)波高的分布圖　單位：m

4.2　無(wú)限小厚度直立島式防波堤附近波浪的傳播

為驗(yàn)證數(shù)值模擬模型的有效性，本文首先對(duì)無(wú)窮小厚度直立島式防波堤附近波浪的傳播變形進(jìn)行數(shù)值模擬。模擬的計(jì)算域?yàn)?20 m×720 m，水深h為10 m。長(zhǎng)度為240 m的無(wú)限小厚度直立式防波堤位于斷面x=240 m的中部。一列線性單頻波沿x正向入射，入射波波高H為0.5 m，波周期為6.0 s。數(shù)值模擬模型的空間步長(zhǎng)和時(shí)間步長(zhǎng)分別為：△x=△y=L/16，△t=T/32，式中L為入射波波長(zhǎng)，m，T為波周期，S。模型計(jì)算的數(shù)值波浪的時(shí)間為33T，波高由t=30T時(shí)刻到t=33T時(shí)刻的波面值按照均方根波高計(jì)算[12]。

上帝有兩個(gè)住所，一個(gè)在天堂，另一個(gè)在感恩的心中。感激他人是表達(dá)你愛(ài)的一種方式。在我們的日常生活中，我們經(jīng)常得到父母、朋友、同事和陌生人的幫助。也許這是一件小事，幫你撿起掉落的筆，給你舉一個(gè)沉重的盒子，或者在公共汽車(chē)上讓給你一個(gè)座位。我們應(yīng)該感謝他們所做的一切。你給予的愛(ài)越多，你得到的愛(ài)就越多。

圖3是t=33T時(shí)刻時(shí)計(jì)算域內(nèi)的瞬時(shí)波面分布圖，圖4是計(jì)算域內(nèi)的相對(duì)波高分布圖。由圖3～4可知，由于直立式防波堤的作用，波浪在防波堤前的區(qū)域發(fā)生了明顯的反射，入射波和反射波相互作用穩(wěn)定后的相對(duì)波高最大可達(dá)2.0 m左右，在防波堤后，行成了一波浪掩護(hù)區(qū)域，在該區(qū)域內(nèi)，波浪的傳播發(fā)生了明顯的繞射現(xiàn)象。由于計(jì)算區(qū)域沿y=360 m中線對(duì)稱(chēng)，所以，圖3～4所示的數(shù)值結(jié)果亦顯示出了良好的對(duì)稱(chēng)性。數(shù)值結(jié)果顯示，本文所建立的數(shù)值模型能夠有效地模擬無(wú)窮小厚度島式防波堤周?chē)ɡ说膫鞑プ冃巍?/p>

圖3直立島式防波堤附近t=33T時(shí)刻的波面平面分布圖單位：m

圖4　直立島式防波堤附近相對(duì)波高平面分布圖　單位：m

4.3　直立方柱附近波浪的傳播

本節(jié)對(duì)直立方柱附近波浪的傳播變形進(jìn)行數(shù)值模擬。模擬的計(jì)算域?yàn)?8.339 m×18.339 m，水深h為0.4 m。尺寸為6.113 m ×6.133 m的方柱位于計(jì)算域的中部。波浪沿x負(fù)方向入射，入射波波高H為0.1 m，波周期為0.9 s。數(shù)值模型的空間步長(zhǎng)和時(shí)間步長(zhǎng)分別為：△x=△y=0.061 13 m，△t=0.015 s。模型計(jì)算的數(shù)值波浪的時(shí)間為50T，波高是通過(guò)計(jì)算t=40T時(shí)刻到t=50T時(shí)刻波面值的均方根得到的。

圖5是t=50T時(shí)刻時(shí)計(jì)算域內(nèi)瞬時(shí)波面的分布圖，圖6是計(jì)算域內(nèi)的相對(duì)波高分布圖。由圖5～6可知，在直立方柱的影響下，波浪的傳播發(fā)生了反射和繞射變形，在方柱的迎浪側(cè)和背浪側(cè)，分別形成了明顯的反射區(qū)域和繞射區(qū)域，在方柱的兩側(cè)，由于方柱側(cè)面與來(lái)波方向相同，所以，兩側(cè)波浪以前進(jìn)波的形式向前傳播。模型的數(shù)值結(jié)果表明，本文模型能夠?qū)Ψ街車(chē)ɡ说膫鞑プ冃螁?wèn)題進(jìn)行研究。

圖5直立方柱附近t=50T時(shí)刻的波面平面分布圖單位：m

圖6　直立方柱附近相對(duì)波高平面分布圖　單位：m

5　結(jié)　語(yǔ)

針對(duì)直立島式結(jié)構(gòu)物附近波浪的傳播變形問(wèn)題，本文基于時(shí)間關(guān)聯(lián)型緩坡方程和相應(yīng)的邊界條件，建立了具有二階精度的適合于復(fù)連通區(qū)域的波浪傳播的數(shù)值模擬模型。該模型不僅適用于變水深問(wèn)題，克服了傳統(tǒng)解析解僅適用于等水深的局限性，而且適用于模擬線性波浪的時(shí)間和空間演化過(guò)程。對(duì)直立島式防波堤以及直立方柱周?chē)ɡ藗鞑プ冃蔚臄?shù)值模擬表明，模型可以用來(lái)模擬直立式結(jié)構(gòu)物影響下波浪傳播過(guò)程中發(fā)生的繞射和反射現(xiàn)象。

參考文獻(xiàn)：

[1]Peregrine,D.H.Long Waves on a Beach[J].J.Fluid Mech.,1967(27): 815-827.

[2]Berkhoff J.C.W.Computation of combined refraction-diffraction[A].Proc 13th Conf on Coastal Eng[C].//Vancouver,Canada: ASCE,1972.

[3]Kirby,J.T.,Changhoon L.,and Chris R.Time-dependent solutions of the mild-slope wave equation[C].//Proceedings of the 23rd International Conference on Coastal Engineering,Venice,Italy.1993.

[5]Copeland GJM.A practical alternative to the mild-slope wave equation[J].Coastal Eng.,1985(9):125-149.

[6]Madsen,P.A.and Larsen,J.An Efficient Finite-Difference Approach to the Mild-Slope Equation[J].Coastal Eng.,1987(11):329-351.

[7]Ebersole,B.A.Refraction-diffraction model for linear water waves[J].J Wtrwy Port Coast and Ocean Eng.,1985，111(6):939-953.

[8]Li B..An evolution equation for water waves[J].Coastal Eng,1994(23):227-242.

[9]Panchang,V.G.，Pearce,B.R..Solution of the mild-slope wave problem by iteration[J].Appl Ocean Res,1991，13(4): 187-199.

[10]張洪生,丁平興,趙海虹.一般曲線坐標(biāo)系下波浪傳播的數(shù)值模擬[J].海洋學(xué)報(bào),2003，25(1): 110-119.

[11]趙明,滕斌.橢圓型緩坡方程的一個(gè)有效的有限元解[J].海洋學(xué)報(bào),2002 ，24(1): 117-123.

[12]Pan J N,Zuo Q H,Wang H C.Efficient numerical solution of the modified mild-slope equation[J].China Ocean Eng,2000,14(2): 161-174.

[13]趙紅軍,張洪生,丁平興.一種求解時(shí)間關(guān)聯(lián)型緩坡方程的數(shù)值方法[J].上海交通大學(xué)學(xué)報(bào),2006， 40(6): 1050-1054.