直立島式結(jié)構(gòu)物周圍波浪傳播的數(shù)值模擬

，　

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

1　問題的提出

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

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

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

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

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

2.1　控制方程

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

(1)

(2)

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

(3)

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

σ2=gktanh(kh)

(4)

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

(5)

(6)

(7)

F(η，Φ)=-gη

(8)

2.2　數(shù)值離散格式

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

Euler預測步：

(9)

(10)

式中：

(11)

(12)

Euler校正步：

(13)

(14)

式中：

(15)

(16)

3　初邊值條件

一般來說，邊界條件有3類，分別是：

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

(18)

(19)

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

(20)

(21)

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

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

(22)

(3)輻射邊界條件

(23)

(24)

4　模型驗證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5　結(jié)　語

針對直立島式結(jié)構(gòu)物附近波浪的傳播變形問題，本文基于時間關(guān)聯(lián)型緩坡方程和相應(yīng)的邊界條件，建立了具有二階精度的適合于復連通區(qū)域的波浪傳播的數(shù)值模擬模型。該模型不僅適用于變水深問題，克服了傳統(tǒng)解析解僅適用于等水深的局限性，而且適用于模擬線性波浪的時間和空間演化過程。對直立島式防波堤以及直立方柱周圍波浪傳播變形的數(shù)值模擬表明，模型可以用來模擬直立式結(jié)構(gòu)物影響下波浪傳播過程中發(fā)生的繞射和反射現(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]張洪生,丁平興,趙海虹.一般曲線坐標系下波浪傳播的數(shù)值模擬[J].海洋學報,2003，25(1): 110-119.

[11]趙明,滕斌.橢圓型緩坡方程的一個有效的有限元解[J].海洋學報,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]趙紅軍,張洪生,丁平興.一種求解時間關(guān)聯(lián)型緩坡方程的數(shù)值方法[J].上海交通大學學報,2006， 40(6): 1050-1054.