倪云林，滕斌，叢龍飛

(1.大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024； 2. 浙江海洋大學(xué) 港航與交通運(yùn)輸工程學(xué)院，浙江 舟山 316022)

?

修正型緩坡方程的有限元模型

倪云林1,2，滕斌1*，叢龍飛1

(1.大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024； 2. 浙江海洋大學(xué) 港航與交通運(yùn)輸工程學(xué)院，浙江 舟山 316022)

與緩坡方程相比，修正型緩坡方程增加了地形曲率項(xiàng)和坡度平方項(xiàng)，從而提高了數(shù)值求解的復(fù)雜性。本文將計(jì)算域劃分為內(nèi)域和外域，內(nèi)域?yàn)樗钭兓瘏^(qū)域，使用修正型緩坡方程，其中的地形曲率項(xiàng)和坡度平方項(xiàng)可用有限單元各節(jié)點(diǎn)的水深信息和單元插值函數(shù)表示，外域?yàn)樗詈愣▍^(qū)，速度勢(shì)滿(mǎn)足Helmholtz方程，通過(guò)內(nèi)外域的邊界匹配建立有限元方程，并用高斯消去法求解。進(jìn)而分別模擬了波浪傳過(guò)Homma島和圓形淺灘的變形，其結(jié)果與相關(guān)的解析解和實(shí)驗(yàn)數(shù)據(jù)吻合良好，證明了本文有限元模型的正確性。同時(shí)，通過(guò)與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比也明顯看出，在地形坡度較陡的情況下，修正型緩坡方程較緩坡方程具有更高的計(jì)算精度。

修正型緩坡方程；有限元；Homma島；圓形淺灘

1 引言

波浪在由深海向淺水傳播的過(guò)程中，受海底地形的影響，會(huì)發(fā)生淺化、折射、繞射等一系列的變形。為了研究這一過(guò)程中波浪要素的變化規(guī)律，許多學(xué)者建立和改進(jìn)了各種數(shù)學(xué)模型。緩坡方程(mild slope equation, MSE)就是其中的一種，通過(guò)沿水深積分的方法將三維問(wèn)題轉(zhuǎn)化為二維問(wèn)題，具有方程形式簡(jiǎn)單且應(yīng)用范圍較廣的優(yōu)點(diǎn)。

緩坡方程最初由Berkhoff(1972年)提出，又稱(chēng)為折射-繞射聯(lián)合模型，是在緩坡假設(shè)的基礎(chǔ)上，用線性波浪理論研究緩變、不可滲透地形上波浪的傳播變形問(wèn)題[1]。關(guān)于該方程的使用條件，Booij(1983年)證明緩坡方程在底坡小于1∶3的情況下具有足夠的精度[2]，但在底坡大于1∶3或地形波動(dòng)的情況下，緩坡方程則不能準(zhǔn)確描述波浪的傳播變形。為了改善“緩坡假設(shè)”的條件限制，Kirby(1986年)最先對(duì)緩坡方程進(jìn)行了改進(jìn)，他將劇變地形分解成一個(gè)緩變地形和一個(gè)小振幅波動(dòng)地形的疊加，導(dǎo)出了擴(kuò)展型緩坡方程(extended mild slope equation，EMSE)，提高了地形波動(dòng)情況下緩坡方程的計(jì)算精度[3]。之后，Chamberlain和Porter[4](1995年)，Chandrasekera和Cheung[5](1997年)通過(guò)引入地形曲率項(xiàng)和坡度平方項(xiàng)對(duì)緩坡方程進(jìn)行了修正，得到了修正型緩坡方程(modified mild slope equation, MMSE)，提高了陡變地形和劇變地形情況下緩坡方程的計(jì)算精度，從而大大拓寬了緩坡方程的應(yīng)用范圍。

目前，有不少文獻(xiàn)對(duì)修正型緩坡方程的求解方法進(jìn)行了論述。例如，Suh等將修正型緩坡方程轉(zhuǎn)化為一階雙曲型方程組，采用有限差分法進(jìn)行求解[6]。Silva等采用中心差分格式離散控制方程和邊界條件來(lái)求解修正型緩坡方程[7]。有限差分法具有數(shù)學(xué)概念直觀、表達(dá)簡(jiǎn)單、發(fā)展成熟等優(yōu)點(diǎn)，但差分網(wǎng)格通常為矩形，在邊界不規(guī)則或形狀復(fù)雜時(shí)精度降低。因此，本文將采用有限元法直接對(duì)修正型緩坡方程進(jìn)行求解，解決地形曲率項(xiàng)和坡度平方項(xiàng)的空間離散問(wèn)題。在求解過(guò)程中，將整個(gè)計(jì)算域劃分為內(nèi)域和外域，內(nèi)域是水深變化的有限區(qū)域，可用有限元離散，外域是水深恒定的無(wú)限區(qū)域，速度勢(shì)滿(mǎn)足Helmholtz方程，其解為級(jí)數(shù)形式，通過(guò)內(nèi)外域公共邊界上的匹配條件建立并求解有限元方程[8—11]。

2 控制方程

針對(duì)水深變化海床上的波浪運(yùn)動(dòng)問(wèn)題，Chamberlain和Porter[4]采用變分原理推導(dǎo)了修正型緩坡方程：

▽·(ccg▽?duì)?+[k2ccg+f1(kh)g▽2h

+f2(kh)gk(▽h)2]φ=0.

(1)

式中，▽為水平梯度算子；φ為速度勢(shì)的空間分量；c和cg分別為波速和群速度；k為波數(shù)；h為水深；g為重力加速度；f1(kh)和f2(kh)分別為地形曲率項(xiàng)系數(shù)和坡度平方項(xiàng)系數(shù)，其表達(dá)式如式(2)、(3)所示[5]。

f1(kh)=[-4khcosh(kh)+sinh(3kh)+sinh(kh)

+8(kh)2sinh(kh)]/{8cosh3(kh)[2kh+sinh(2kh)]}

(2)

16(kh)3sinh(2kh)-9sinh2(2kh)cosh(2kh)+

12(kh)[1+2sinh4(kh)][kh+sinh(2kh)]}.

(3)

上述兩個(gè)系數(shù)均為kh的函數(shù)。從圖1可以看出，f1(kh)和f2(kh)在有限水深區(qū)數(shù)值較大，而隨著水深的增大，兩者均趨于0。

圖1 f1和f2隨相對(duì)水深kh的變化Fig.1 Variation of coefficients f1 and f2 with kh

當(dāng)忽略地形曲率項(xiàng)和坡度平方項(xiàng)時(shí)，修正型緩坡方程就簡(jiǎn)化為Berkhoff的緩坡方程：

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

(4)

當(dāng)水深不變時(shí)，則轉(zhuǎn)換為Helmholtz方程：

▽2φ+k2φ=0.

(5)

3 有限元方程的建立

3.1 計(jì)算域的劃分

如圖2所示，將計(jì)算域劃分為內(nèi)域和外域。其中，內(nèi)域Ω1是地形變化區(qū)域，劃分六節(jié)點(diǎn)三角形單元網(wǎng)格，用單元插值函數(shù)作為權(quán)函數(shù)以建立有限元方程；外域Ω2是水深恒定區(qū)域，不劃分網(wǎng)格，用Hankel函數(shù)作為權(quán)函數(shù)，并應(yīng)用格林公式將外域的面積分轉(zhuǎn)化為邊界Ba上的線積分以建立方程式；Ba為內(nèi)外域的公共邊界，通常情況下取為圓周，通過(guò)Ba上的匹配條件可進(jìn)行有限元方程的求解；Bb為物面邊界；B∞為無(wú)窮遠(yuǎn)邊界。

圖2 計(jì)算域劃分Fig.2 Division of computational domain

3.2 邊界條件

設(shè)內(nèi)域Ω1的速度勢(shì)為φ，則在物面邊界Bb上滿(mǎn)足

(6)

(7)

(8)

式中，r為Ω2域內(nèi)x-y平面上某點(diǎn)到原點(diǎn)O的距離；φ為r與x軸正方向的夾角；Hn(kr)為n階第1類(lèi)Hankel函數(shù)；An(n=0,1,2，…)和Bn(n=1,2，…)為待定系數(shù)。在n=N處截?cái)?，并令G1=H0(kr)，G2i=Hi(kr)cosiφ，G2i+1=Hi(kr)siniφ，α1=A0，α2i=Ai，α2i+1=Bi(i=1,2,…，N)，式(8)可簡(jiǎn)寫(xiě)成：

(9)

式中，M=2N+1。則無(wú)窮遠(yuǎn)處輻射邊界條件可以表示為:

(10)

在公共邊界Ba上滿(mǎn)足：

(11)

3.3 內(nèi)域有限元方程

設(shè)φj為單元各節(jié)點(diǎn)速度勢(shì)，Nj為單元插值函數(shù)，則單元內(nèi)的速度勢(shì)：

(12)

利用單元各節(jié)點(diǎn)的水深信息和單元插值函數(shù)，坡度平方項(xiàng)可表示為:

(13)

(14)

(15)

類(lèi)似的，地形曲率項(xiàng)

(16)

取單元插值函數(shù)Ni為權(quán)函數(shù)，利用Galerkin方法在單元內(nèi)建立積分方程

[▽·(ccg▽?duì)?+(k2ccg+f1(kh)g▽2h

+f2(kh)gk(▽h)2)φ]NidΩ=0.

(17)

引入三節(jié)點(diǎn)線單元插值函數(shù)Li，對(duì)邊界Ba上的單元而言，Li為該單元外側(cè)邊界上的單元插值函數(shù)。將邊界Ba上的速度勢(shì)分解為入射勢(shì)φI和繞射勢(shì)φS，應(yīng)用格林公式，式(17)可以重新寫(xiě)成:

(18)

結(jié)合式(12)、(15)和(16)，可以寫(xiě)成單元矩陣形式

(19)

其中,

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

(20)

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

3.4 外域方程式[12—13]

(21)

應(yīng)用第二格林公式，式(21)可以寫(xiě)成:

(22)

(23)

(24)

3.5 總矩陣方程

根據(jù)式(20)和式(24)，可以形成總矩陣方程

[K](q+M)×(q+M){φ}(q+M)×1={B}(q+M)×1.

(25)

式中，

；

由于Ba為圓周，所以[K22]可用解析表達(dá)式表示

[K22]=2πkRccgdiag{H0′(kR)H0(kR),

H1′(kR)H1(kR),…,HM′(kR)HM(kR)},

(26)

式中,R為圓周Ba的半徑。

最后，本文采用列主元高斯消去法求解總矩陣方程式(25)。

4 結(jié)果驗(yàn)證和討論

為了驗(yàn)證本文修正型緩坡方程的有限元模型，首先計(jì)算了長(zhǎng)波經(jīng)過(guò)Homma島的繞射問(wèn)題，計(jì)算結(jié)果與Zhai等[14]的解析解進(jìn)行了比較；接著計(jì)算了圓形淺灘上波浪的繞射問(wèn)題，計(jì)算結(jié)果與Suh等[6]的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了比較，并進(jìn)一步探討了較陡地形對(duì)不同周期波浪的影響。

4.1 長(zhǎng)波過(guò)Homma島的繞射問(wèn)題

Homma島由旋轉(zhuǎn)拋物體及上部圓柱體組合而成。如圖3所示，上部圓柱體的半徑a=10 km，旋轉(zhuǎn)拋物體的底面半徑b=30 km。距島嶼中心r(單位：km)處的水深(單位：km)為:

(27)

圖3 Homma島地形Fig.3 Sketch of Homma Island

圖4 岸線圓周上修正型緩坡方程計(jì)算的相對(duì)波高和解析解的對(duì)比Fig.4 Comparison of calculated wave height by MMSE with analytical solutions along the shoreline

4.2 波浪過(guò)圓形淺灘的繞射問(wèn)題

為了進(jìn)一步驗(yàn)證較陡地形上修正型緩坡方程有限元模型的正確性和應(yīng)用范圍，本文將修正型緩坡方程有限元解、緩坡方程的有限元解與Sub等[6]的實(shí)驗(yàn)結(jié)果進(jìn)行比較。實(shí)驗(yàn)設(shè)計(jì)的圓形淺灘地形如圖5所示，距淺灘中心r(0≤r≤R，單位：m)處的水深(單位：m)為:

(28)

圖5 圓形淺灘地形Fig.5 Sketch of circular shoal

圖6是修正型緩坡方程和緩坡方程計(jì)算的y=0斷面上相對(duì)波高與實(shí)驗(yàn)結(jié)果的對(duì)比。從圖中可以看出，相比緩坡方程，修正型緩坡方程的計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)更加吻合，這說(shuō)明在地形較陡的情況下，地形曲率項(xiàng)和坡度平方項(xiàng)對(duì)波浪的傳播變形會(huì)起到十分重要的作用。通過(guò)比較圖6a～c，并結(jié)合表1可以發(fā)現(xiàn)，在相對(duì)水深較大時(shí)，波高最大值發(fā)生在灘頂后方，而隨著相對(duì)水深的減小，波高的最大值有所增加，發(fā)生最大值的位置向前方移動(dòng)，緩坡方程的計(jì)算結(jié)果在灘頂后方也表現(xiàn)出更大的誤差，這說(shuō)明長(zhǎng)波較短波更易受到地形曲率和坡度平方的影響。

5 結(jié)論

(1)本文建立了修正型緩坡方程的有限元模型，并應(yīng)用該模型計(jì)算了波浪過(guò)Homma島和圓形淺灘的繞射問(wèn)題，計(jì)算結(jié)果與翟西媛等的解析解和Sub等的實(shí)驗(yàn)結(jié)果吻合良好，證明了本文有限元模型的正確性。

圖6 y=0斷面修正型緩坡方程和緩坡方程計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比Fig.6 Comparison of MMSE and MSE results with experimental data along y=0 section

波浪參數(shù)最大相對(duì)波高Hmax/H0最大相對(duì)波高發(fā)生位置x/R實(shí)驗(yàn)結(jié)果MMSE計(jì)算結(jié)果MSE計(jì)算結(jié)果實(shí)驗(yàn)結(jié)果MMSE計(jì)算結(jié)果MSE計(jì)算結(jié)果T=1.259s,kh0=0.9651.2621.2561.239-0.220-0.178-0.022T=0.791s,kh0=2.0031.2521.2351.3130.6670.7110.800T=0.636s,kh0=3.0031.1971.1821.2311.3251.0000.889

(2)通過(guò)本文修正型緩坡方程計(jì)算結(jié)果、緩坡方程計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)的三者對(duì)比，可以看出，對(duì)于較陡地形上的波浪變形問(wèn)題，修正型緩坡方程較緩坡方程具有更高的計(jì)算精度，這從側(cè)面證明了本文的有限元模型具有更寬的適用范圍。

(3)針對(duì)長(zhǎng)波入射的情況，Homma島上部圓柱周?chē)淖畲蟛ǜ卟皇前l(fā)生在迎浪點(diǎn)，而是對(duì)稱(chēng)分布于迎浪點(diǎn)左右兩側(cè)，這與直立圓柱的情形不同。

(4)在地形坡度較陡的情況下，地形曲率項(xiàng)和坡度平方項(xiàng)會(huì)對(duì)波浪的傳播變形起到不可忽略的影響，并且隨著波浪周期的增大，波浪受地形曲率和坡度平方的影響也隨之增大。

[1] Berkhoff J C W. Computation of Combined Refraction-diffraction[M]. Delft Hydraulics Laboratory, 1974.

[2] Booij N. A note on the accuracy of the mild-slope equation[J]. Coastal Engineering, 1983, 7(3): 191-203.

[3] Kirby J T. A general wave equation for waves over rippled beds[J]. Journal of Fluid Mechanics, 1986, 162: 171-186.

[4] Chamberlain P G, Porter D. The modified mild-slope equation[J]. Journal of Fluid Mechanics, 1995, 291: 393-407.

[5] Chandrasekera C N, Cheung K F. Extended linear refraction-diffraction model[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 1997, 123(5): 280-286.

[6] Suh K D, Lee C, Park Y H, et al. Experimental verification of horizontal two-dimensional modified mild-slope equation model[J]. Coastal Engineering, 2001, 44(1): 1-12.

[7] Silva R, Borthwick A G L, Taylor R E. Numerical implementation of the harmonic modified mild-slope equation[J]. Coastal engineering, 2005, 52(5): 391-407.

[8] 李孟國(guó), 蔣德才. 關(guān)于波浪緩坡方程的研究[J]. 海洋通報(bào), 1999, 18(4): 70-92.

Li Mengguo, Jiang Decai. A review on the study of mild-slope equation[J]. Marine Science Bulletin, 1999,18(4): 70-92.

[9] Chen H S, Mei C C. Oscillations and wave forces in a man-made harbor in the open sea[C]//Symposium on Naval Hydrodynamics, 10th, Proceedings, Pap and Discuss, Cambridge, Mass, June 24-28, 1974. 1976.

[10] Tsay T K, Liu P L F. A finite element model for wave refraction and diffraction[J]. Applied Ocean Research, 1983, 5(1): 30-37.

[11] Houston J R. Combined refraction and diffraction of short waves using the finite element method[J]. Applied Ocean Research, 1981, 3(4): 163-170.

[12] 趙明. 波浪作用下建筑物周?chē)哪嗌硾_刷及海床演變[D]. 大連: 大連理工大學(xué), 2002.

Zhao Ming. The local scour and topographical change around offshore structures under wave action[D]. Dalian: Dalian University of Technology, 2002.

[13] 趙明, 滕斌. 緩坡方程的有限元解[J]. 大連理工大學(xué)學(xué)報(bào), 2000, 40(1): 117-119.

Zhao Ming, Teng Bin. Finite element solutions for mild slope equation[J]. Journal of Dalian University of Technology, 2000, 40(1): 117-119.

[14] Zhai X Y, Liu H W, Xie J J. Analytic study to wave scattering by a general Homma island using the explicit modified mild-slope equation[J]. Applied Ocean Research, 2013, 43: 175-183.

[15] Chau F P, Taylor R E. Second-order wave diffraction by a vertical cylinder[J]. Journal of Fluid Mechanics, 1992, 240(1): 571-599.

FEM model of the modified mild slope equation

Ni Yunlin1,2，Teng Bin1，Cong Longfei1

(1.StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian116024,China；2.SchoolofPortandTransportationEngineering,ZhejiangOceanUniversity,Zhoushan316022,China)

Compared with the mild slope equation, the bottom curvature and slope-squared terms are contained in the modified mild slope equation, which increases the complexity of solving the equation numerically. In this paper, the physical domain is divided into a finite inner region and an infinite outer region. The inner region of a variable depth is studied with the modified mild slope equation, and the outer region has a constant water depth, the velocity potential satisfying the Helmholtz equation. The bottom curvature and slope-squared terms in the equation are evaluated from the input bathymetry at each node of every element using the interpolation functions. With the boundary matching method and Gaussian elimination technique, the finite element equations are established and solved. Then, wave transformation over a Homma Island and a circular shoal is simulated respectively, and the results agree well with analytical solutions and experimental data, verifying the validity of the FEM model in this paper. Meanwhile, the comparison between the numerical and experimental results indicates the modified mild slope equation can provide more accurate predictions than the mild slope equation for wave propagation over relatively steep bathymetry.

the modified mild slope equation; finite element method; Homma Island; circular shoal

10.3969/j.issn.0253-4193.2017.01.011

2016-03-07；

2016-05-17。

國(guó)家自然科學(xué)基金(51379032，51490672)。

倪云林(1986—)，男，浙江省舟山市人，講師，博士生，主要從事波浪對(duì)海上建筑物作用的研究工作。E-mail:nylzjou@126.com

*通信作者：滕斌，教授，主要從事波浪對(duì)海上建筑物的作用研究。E-mail: bteng@dlut.edu.cn

P753

A

0253-4193(2017)01-0104-07

倪云林, 滕斌, 叢龍飛. 修正型緩坡方程的有限元模型[J]. 海洋學(xué)報(bào), 2017, 39(1):104-110,

Ni Yunlin，Teng Bin，Cong Longfei. FEM model of the modified mild slope equation[J]. Haiyang Xuebao, 2017, 39(1): 104-110, doi: 10.3969/j.issn.0253-4193.2017.01.011