沈 遠(yuǎn)，陳相宇，張希疆，葉 云，倪云林*

(1.舟山市自然資源測(cè)繪設(shè)計(jì)中心，浙江 舟山 316021；2.浙江海洋大學(xué)海洋工程裝備學(xué)院，浙江 舟山 316022)

0 引言

波浪在從外海向近岸傳播過(guò)程中，遇到出露海面的島嶼就會(huì)發(fā)生繞射，從而引起波高在島嶼岸線及周邊海域的分布變化。準(zhǔn)確模擬波高的分布情況對(duì)于海岸及近海工程的規(guī)劃設(shè)計(jì)、施工安全和運(yùn)營(yíng)管理具有十分重要的意義。緩坡類方程是一種用于計(jì)算近海海域波浪傳播變形的數(shù)學(xué)模型，可以綜合考慮波浪的折射與繞射，具有計(jì)算效率高、計(jì)算精度好等優(yōu)點(diǎn)。

原始的緩坡方程由荷蘭Delft水力學(xué)研究所的BERKHOFF[1]在1972年推導(dǎo)得到。他基于線性波浪理論和緩坡假定，通過(guò)沿水深積分，得到了一個(gè)二維橢圓型偏微分方程，并將其命名為緩坡方程(mild slope equation)。該方程形式簡(jiǎn)單，應(yīng)用范圍廣，因此引起了許多學(xué)者的研究興趣。BOOIJ[2]利用緩坡方程計(jì)算了斜坡海底對(duì)波浪的反射，闡明了對(duì)于坡度小于1∶3的緩坡而言，緩坡方程的計(jì)算結(jié)果具有足夠的精度。KIRBY[3]將沙波地形分解為一個(gè)緩變地形和一個(gè)小振幅快變地形的疊加，推導(dǎo)了擴(kuò)展型緩坡方程(extended mild slope equation)，從而改進(jìn)了緩坡條件的限制，拓寬了緩坡方程在斜坡小振幅快變地形上的應(yīng)用。CHAMBERLAIN et al[4]在推導(dǎo)過(guò)程中保留了坡度平方項(xiàng)(sloped-squared term)和地形曲率項(xiàng)(bottom curvature term)，得到了修正型緩坡方程(modified mild slope equation),不僅提高了緩坡方程在陡變和快變地形上的計(jì)算精度，而且進(jìn)一步拓寬了緩坡方程的應(yīng)用范圍。

許多學(xué)者應(yīng)用緩坡類方程研究了波浪在島嶼地形上繞射的問(wèn)題。HOMMA[5]給出了波浪經(jīng)過(guò)拋物形淺灘上圓形島(即Homma島)的級(jí)數(shù)解。LIU et al[6]采用Padé逼近方法，推導(dǎo)了有限水深波經(jīng)過(guò)Homma島情況下緩坡方程的逼近解析解。隨后，LIN et al[7]、LIU et al[8]、NIU et al[9]、KUO et al[10]采用解析方法開展了波浪在島嶼地形上繞射的研究。解析解的提出為數(shù)值解提供了驗(yàn)證方法，但通常只有在某些規(guī)則的島嶼地形下才能得到。在數(shù)值研究方面，劉紅玲[11]以洪廣文推導(dǎo)的非定常線性波緩坡方程為基礎(chǔ)，采用經(jīng)典的ADI法開展了島嶼水域波浪場(chǎng)的數(shù)值模擬。張軍 等[12]提出了島體邊界條件，利用緩坡方程建立了計(jì)算域中存在直立島式結(jié)構(gòu)物時(shí)波浪傳播的數(shù)值模擬模型。本文將采用混合元方法對(duì)修正型緩坡方程進(jìn)行數(shù)值求解，在對(duì)其進(jìn)行驗(yàn)證的基礎(chǔ)上，詳細(xì)探討工程尺度背景下，波浪入射周期、淺灘形狀參數(shù)和島嶼尺寸等因素對(duì)波浪繞射的影響。

1 控制方程和數(shù)值解法

1.1 控制方程

BERKHOFF推導(dǎo)得到的緩坡方程為

▽·(ccg▽?duì)?+k2ccgφ=0

(1)

(2)

6kh(2kh+2sinh2kh)(cosh22kh-2cosh2kh+3)]

(3)

分別表示地形曲率項(xiàng)系數(shù)和坡度平方項(xiàng)系數(shù)。由此得到修正型緩坡方程：

▽·(ccg▽?duì)?+{k2ccg+g[u1(h)▽2h+

(4)

1.2 數(shù)值解法

圖1 內(nèi)域和外域劃分示意圖

(5)

其中:

(6)

(7)

(8)

內(nèi)域?yàn)榈匦巫兓瘏^(qū)域，用六節(jié)點(diǎn)三角形單元進(jìn)行劃分，速度勢(shì)φ可以寫成：

(9)

其中：φi是三角形單元各節(jié)點(diǎn)速度勢(shì)，Ni是三角形單元各節(jié)點(diǎn)的形函數(shù)。

通過(guò)等參變換并采用Galerkin方法可建立單元Ω內(nèi)的矩陣方程：

(10)

式中：

(11)

(12)

(13)

(14)

式中：Nj、Nm是三角形單元各節(jié)點(diǎn)的形函數(shù)，hm是三角形單元各節(jié)點(diǎn)的水深，Li為內(nèi)、外域公共邊界s上的線單元形函數(shù)，n為內(nèi)、外域公共邊界s上的單位法向向量。

集合內(nèi)域的各單元矩陣可得：

(15)

式中：q為總節(jié)點(diǎn)數(shù)，l為公共邊界上的節(jié)點(diǎn)數(shù)。

通過(guò)內(nèi)域和外域公共邊界上的速度勢(shì)連續(xù)和法向速度連續(xù)的匹配條件，可以建立方程求解的線性方程組，并采用高斯消去法求解。詳細(xì)求解過(guò)程可參考文獻(xiàn)[13]。

2 模型驗(yàn)證與討論

2.1 模型驗(yàn)證

KUO et al[10]利用Bessel函數(shù)求解了波浪在三維圓形島地形上繞射的長(zhǎng)波方程解析解。如圖2所示，三維圓形島由半徑為r0的圓形淺灘和其上半徑為r1的圓柱形島嶼組成，距離淺灘中心r處的水深為

圖2 三維圓形島圖示

(16)

式中：h0是三維圓形島周圍的常水深，β是淺灘形狀參數(shù)。計(jì)算時(shí)，KUO et al[10]取r0=0.2π，r1=0.1π，h0=0.02，β先后取1和2，入射波浪波高為Hi，周期T=2，對(duì)應(yīng)的波長(zhǎng)L=0.883。

本文同樣計(jì)算了上述參數(shù)條件下的兩個(gè)算例(算例一中β取1，算例二中β取2)，對(duì)比計(jì)算得到的y=0斷面上的相對(duì)波高[即沿程波高H(x)與入射波高Hi的幅值比]和KUO et al的解析解(圖3)。由圖可見(jiàn)，在淺灘形狀參數(shù)β=1和2的兩種情況中，本文的數(shù)值解與KUO et al的解析解均吻合很好，證明了本文數(shù)學(xué)模型的正確性。

圖3 圓形島y=0斷面上相對(duì)波高本文數(shù)值解與KUO et al解析解的對(duì)比

2.2 討論

針對(duì)波浪在三維圓形島地形上繞射的問(wèn)題，本節(jié)主要討論波浪入射周期、淺灘形狀參數(shù)和島嶼尺寸等因素對(duì)波浪繞射的影響，并按工程尺度選取計(jì)算參數(shù)。

2.2.1 波浪入射周期的影響分析

首先，取r0=200π m，r1=100π m，h0=20 m，β=1；考慮長(zhǎng)周期重力波沿x軸入射，周期T先后取60 s、80 s、100 s，對(duì)應(yīng)的波長(zhǎng)L=836.860 m、1 117.659 m、1 398.135 m，h0/L=0.024、0.018、0.014。在不同周期的波浪入射情況下，y=0斷面上的相對(duì)波高分布如圖4所示，可以看出，隨著入射波浪周期的變大，島嶼迎浪側(cè)相對(duì)波高最大值有所減小，最小值有所增大；島嶼背浪側(cè)的相對(duì)波高隨之減小。但隨著波浪遠(yuǎn)離島嶼，不同周期波浪的相對(duì)波高差異減小。

圖4 不同周期波浪在圓形島y=0斷面上的相對(duì)波高分布

進(jìn)一步研究不同周期的波浪在島嶼岸線上的相對(duì)波高分布，結(jié)果繪于圖5。對(duì)于不同周期的波浪而言，相對(duì)波高均關(guān)于x軸對(duì)稱分布。當(dāng)T=60 s時(shí)，相對(duì)波高最大值為2.493，位于θ=180°處(即迎浪點(diǎn))，相對(duì)波高最小值為0.452，位于θ=25°和335°處(即背浪點(diǎn)兩側(cè)25°)；當(dāng)T=80 s時(shí)，相對(duì)波高最大值為2.314，位于θ=161°和199°處(即迎浪點(diǎn)兩側(cè)19°)，相對(duì)波高最小值為0.533，位于θ=30°和330°處(即背浪點(diǎn)兩側(cè)30°)；當(dāng)T=100 s時(shí)，相對(duì)波高最大值為2.273，位于θ=180°處(即迎浪點(diǎn))，相對(duì)波高最小值為0.596，位于θ=35°和325°處(即背浪點(diǎn)兩側(cè)35°)。這說(shuō)明隨著入射波浪周期的變大，相對(duì)波高最大值隨之減小，發(fā)生最大值的位置位于迎浪點(diǎn)及其附近；相對(duì)波高最小值則隨之增大，發(fā)生最小值的位置有遠(yuǎn)離背浪點(diǎn)的趨勢(shì)。

圖5 不同周期波浪在圓形島岸線上的相對(duì)波高分布

2.2.2 淺灘形狀參數(shù)的影響分析

接著，取r0=200π m，r1=100π m，h0=20 m，T=60 s，淺灘形狀參數(shù)β分別取0.5、0.8、1.0、1.5。入射波浪同樣沿x軸正方向傳播，不同淺灘形狀系數(shù)情況下，y=0斷面上的相對(duì)波高分布如圖6所示。可以看出，隨著淺灘形狀參數(shù)β的增大，島嶼迎浪側(cè)相對(duì)波高振蕩幅度加劇、最大值增加，島嶼背浪側(cè)的相對(duì)波高也隨之增加，這和淺灘形狀參數(shù)增大，島嶼周圍水深減小，波浪淺水變形增強(qiáng)密切相關(guān)。

圖6 不同淺灘形狀情況下圓形島y=0斷面上的相對(duì)波高分布

同樣，研究了不同淺灘形狀參數(shù)情況下，波浪在島嶼岸線上的相對(duì)波高分布，結(jié)果繪于圖7。由此可見(jiàn)，島嶼岸線上的相對(duì)波高隨著淺灘形狀參數(shù)β的增大而增加，且相對(duì)波高的大小分布差異也隨之增加；當(dāng)β=0.5、0.8和1.0時(shí)，相對(duì)波高最大值均位于迎浪點(diǎn)(θ=180°)，分別為2.133、2.350和2.493；相對(duì)波高最小值均位于背浪點(diǎn)兩側(cè)約25°處(θ=24.5°和335.5°)，分別為0.412、0.435和0.452。當(dāng)β=1.5 時(shí)，相對(duì)波高最大值為2.736，但不再位于迎浪點(diǎn)，而位于迎浪點(diǎn)兩側(cè)約25°處(θ=155.5°和 204.5°)；最小值為0.502，位于背浪點(diǎn)兩側(cè)約22°處(θ=21.8°和338.2°)。

圖7 不同淺灘形狀情況下波浪在圓形島岸線上的相對(duì)波高分布

2.2.3 島嶼尺寸的影響分析

最后，取r0=200π m，h0=20 m，T=60 s，β=1.0，圓形島半徑r1分別取50π m、100π m、150π m，r1/r0=0.25、0.50、0.75。入射波浪依舊沿x軸正方向傳播，不同島嶼尺寸情況下，y=0斷面上的相對(duì)波高分布繪于圖8。由圖可知，隨著島嶼尺寸的減小，島嶼周邊的水深減小，波浪的淺水變形加劇，導(dǎo)致島嶼迎浪點(diǎn)及后方背浪側(cè)的相對(duì)波高增加。從不同尺寸島嶼岸線上的相對(duì)波高分布情況(圖9)也可以看出，隨著島嶼尺寸的減小，島嶼岸線上的相對(duì)波高總體上隨之增加；當(dāng)r1=50π m、100π m、150π m時(shí)，即r1/r0=0.25、0.50、0.75，相對(duì)波高最大值分別為 3.147、2.493、2.076，均位于迎浪點(diǎn)(θ=180°)；最小值分別為0.803、0.515、0.269，并呈現(xiàn)出向背浪點(diǎn)靠近的變化趨勢(shì)，即最小值發(fā)生位置分別為θ=31.8°和328.2°、θ=27.3°和332.7°以及θ=19.8°和340.2°。

圖8 不同島嶼尺寸情況下圓形島y=0斷面上的相對(duì)波高分布

圖9 波浪在不同尺寸圓形島岸線上的相對(duì)波高分布

3 結(jié)論

本文通過(guò)混合元方法求解修正型緩坡方程，開展了波浪在三維圓形島地形上繞射的研究，探討了在工程尺度情況下，波浪入射周期、淺灘形狀參數(shù)和島嶼尺寸等因素對(duì)波浪繞射的影響，結(jié)果表明：

(1)在沿波浪傳播方向上，圓形島迎浪側(cè)的相對(duì)波高振蕩幅度隨入射波周期的減小、淺灘形狀參數(shù)的增大和島嶼尺寸的減小而增大。波浪經(jīng)過(guò)圓形島，其背浪側(cè)的相對(duì)波高同樣隨入射波周期的減小、淺灘形狀參數(shù)的增大和島嶼尺寸的減小而增大，并且隨著傳播距離的增加而趨于一致。

(2)圓形島嶼岸線上的相對(duì)波高總體上隨入射波周期的減小、淺灘形狀參數(shù)的增大和島嶼尺寸的減小而增大；相對(duì)波高最大值大多數(shù)發(fā)生在迎浪點(diǎn)，個(gè)別發(fā)生在迎浪點(diǎn)兩側(cè)20°～25°處；相對(duì)波高最小值發(fā)生在背浪點(diǎn)兩側(cè)30°附近。這是由于當(dāng)波浪的波峰靠近圓形島時(shí)，迎浪點(diǎn)及兩側(cè)的水位抬高，相對(duì)波高達(dá)到最大值；水流沿圓形島兩側(cè)流向背浪點(diǎn)并最終相遇，動(dòng)能轉(zhuǎn)為勢(shì)能，相對(duì)波高在背浪點(diǎn)及附近發(fā)生最小值。