波浪在不平坦海岸底床上傳播的平面二維數(shù)學(xué)模型

彭延建，劉應(yīng)中，時(shí)鐘

（上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 港口與海岸工程系，海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200030）

波浪在不平坦海岸底床上傳播的平面二維數(shù)學(xué)模型

彭延建，劉應(yīng)中，時(shí)鐘*

（上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 港口與海岸工程系，海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200030）

基 于線(xiàn)性波 浪 理論，利用攝動(dòng) 理 論 中 的 多 重 尺度 法 ， 以坡度 ε 為小 參 數(shù) ， Liu 和 Shi[1]推導(dǎo)出 波 浪 在不 平 坦底床傳播的波幅方程，以及一階勢(shì)函數(shù)的解析解。為驗(yàn)證該模型是否可用于平面二維波浪傳播問(wèn)題，利用有限差分法，分別對(duì)波幅方程和邊界條件進(jìn)行離散和求解，得到因折射、繞射引起的波幅和波向的變化。通過(guò)引入繞射因子，使該 模型的應(yīng)用范圍更為廣泛，并且計(jì)算精度得到提高。將計(jì)算得到 的波向和波 高與 Berkhoff[2]、Berkhoff 等[3]橢圓形淺灘試驗(yàn)值進(jìn)行比較，表明該模型可以應(yīng)用于平面二維問(wèn)題。

多重尺度法；波浪傳播；不平坦底床

0 引言

波浪在不平坦海岸底床上的傳播是波浪研究的重要課題之一。波浪在從深水區(qū)向近岸傳播過(guò)程中，因?yàn)楹０兜匦巫兓戎T多因素，會(huì)發(fā)生淺化、折射與繞射等現(xiàn)象，波況也隨之變化，比如：波長(zhǎng)變短、波速變慢、波高先減小后增大、波向也會(huì)趨于與岸線(xiàn)垂直。本文擬對(duì)線(xiàn)性波浪在不平坦海岸底床上的傳播進(jìn)行初步研究，對(duì)海岸工程學(xué)科具有科學(xué)意義。

對(duì)波浪在不平底床傳播問(wèn)題的研究[4-9]，大體有兩種研究方法：

第一種是基于緩坡方程和 Boussinesq 方程等已經(jīng)較為成熟的模型，尋求數(shù)學(xué)模型和數(shù)值求解方法的改進(jìn)，比如在緩坡方程中引入底摩擦、波浪破碎、非線(xiàn)性、波流相互作用等效應(yīng)，以及對(duì)緩坡方程的簡(jiǎn)化和近似，例如雙曲形緩坡方程、拋物形緩坡方程等，國(guó)內(nèi)外眾多學(xué)者在這方面的努力使基于緩坡方程和 Boussinesq 方程的模型在計(jì)算精度、計(jì)算速度、穩(wěn)定性等方面有了較大的提高，已經(jīng)能夠較為成功地應(yīng)用于某些工程案例。這類(lèi)模型也存在一些缺點(diǎn)，比如有的只適合大面積計(jì)算區(qū)域，有的只適合小面積計(jì)算區(qū)域，求解短周期波浪場(chǎng)需要較高的網(wǎng)格密度，需要較大的計(jì)算機(jī)內(nèi)存，求解過(guò)程也較為耗時(shí)等。

第二種是理論解析的方法，例如通過(guò)直接求解底部邊界條件、線(xiàn)性自由表面條件來(lái)推導(dǎo)波浪勢(shì)函數(shù)的解，或通過(guò)攝動(dòng)展開(kāi)求波浪勢(shì)函數(shù)的近似解。這類(lèi)模型存在的問(wèn)題是，給出的解析解形式一般是局限在某些特定條件下，例如二維的斜坡床底等特定水深地形。

根據(jù)以上論述，對(duì)波浪在不平底床傳播問(wèn)題的研究，多數(shù)研究還是基于緩坡方程和Boussinesq 方程，進(jìn)行模型的改進(jìn)和數(shù)值方法完善，而在進(jìn)一步發(fā)展波浪在不平岸灘傳播的數(shù)學(xué)模型方面的研究較少。所以，進(jìn)一步發(fā)展波浪在不平底床傳播的數(shù)學(xué)模型，對(duì)波浪在不平底床傳播的研究是有意義的。

本文的目的為通過(guò)引入繞射因子，進(jìn)一步擴(kuò)展 Liu 和 Shi[1]提出的模型，并對(duì)其進(jìn)行驗(yàn)證。

1 數(shù)學(xué)模型

本文采用的是 Liu 和 Shi[1]中的數(shù)學(xué)模型，為了全文的完整性，本章將數(shù)學(xué)模型的推導(dǎo)過(guò)程在此做簡(jiǎn)要介紹。本文數(shù)學(xué)模型的推導(dǎo)基于經(jīng)典的線(xiàn)性波理論。以笛卡爾坐標(biāo)系為參考坐標(biāo)，水平x軸位于靜水面，與岸線(xiàn)的交點(diǎn)為坐標(biāo)原點(diǎn)，指向外海。z軸垂直于靜水面，向上為正。靜水水深為（x，y）的函數(shù)，即 h=h（x，y）。假設(shè)外海入射波為單色波，無(wú)窮遠(yuǎn)處水深不變，比如 z=-h0，h0是常數(shù)。其速度勢(shì)函數(shù)如下：

式中：φ 為速度勢(shì)；g為重力加速度；a0為深海波幅；i為虛數(shù)單位，ω 為波浪頻率；k0為深海波數(shù)。

ζ為波面函數(shù)，代表自由水位高程，表達(dá)式為：

波浪的頻率 ω 與波數(shù) k0滿(mǎn)足色散關(guān)系式：

波浪在向近岸傳播過(guò)程中會(huì)因地形的變化發(fā)生變形效應(yīng)，引起折射和繞射，如果地形變化劇烈還可能發(fā)生反射。

假設(shè)海底地形為緩變海底，（|▽h|/（kh））為一小量ε，即在一個(gè)波長(zhǎng)范圍內(nèi)，h的變化很小。引入兩種尺度，長(zhǎng)尺度 （x?，y?）= ε （x，y） 和 短尺度（x，y），前者對(duì)應(yīng)于地形變化的尺度，后者對(duì)應(yīng)于波浪傳播的尺度。利用多重尺度法，可以得到：

z方向只有一種尺度。假設(shè)解具有以下形式：

其中：φ0為零階近似解；φ1為一階近似解。將這些關(guān)系代入通常的線(xiàn)性波浪問(wèn)題的控制方程和邊界條件，可以得出：φ0滿(mǎn)足拉普拉斯方程和自由表面邊界條件及底部邊界條件：

φ1滿(mǎn)足：

無(wú)窮遠(yuǎn)處：

式 （6d） 說(shuō)明，在無(wú)窮遠(yuǎn)處只存在零階近似解。此處的一階近似解是在線(xiàn)性波理論框架下考慮緩坡的高階近似。

從式 （6a） ~ （6d） 可解得一階近似解。對(duì) x求偏導(dǎo)：

對(duì) y，y?的推導(dǎo)和 x，x?類(lèi)似。線(xiàn)性色散關(guān)系兩邊分別對(duì) x?和y?求偏導(dǎo)，得到：

代入式（6a）~（6d），一階近似解的方程和條件變?yōu)椋?/p>

無(wú)窮遠(yuǎn)處，φ1=0。

其中：

以上推導(dǎo)中的參數(shù)均為關(guān)于水平尺度（x，y）的函數(shù)。假設(shè)一階解的形式為：

將式 （10） 代入式 （9a），比較等式兩邊的系數(shù)，得到：

需要指出的是，系數(shù) αi，（i=1，2，3，4，5）均為實(shí)數(shù)，零階近似解的波幅函數(shù)也在實(shí)數(shù)域內(nèi)，所以一階勢(shì)函數(shù) φ1= φa1meiS，幅值函數(shù) φa1m也是實(shí)數(shù)。

結(jié)合色散關(guān)系式，將式 （10） 代入式 （9b），自由表面邊界條件為：

將 αi（i=1，2，5）代入式 （12）：

將 A ，C1，C2的表達(dá)式代入式 （13），得到波幅方程：

這里的波幅函數(shù)是零階近似解的波幅。

2 計(jì)算方法

在每 個(gè)網(wǎng) 格點(diǎn) ，k（x?，y?，t）與當(dāng)?shù)?深度 h（x?，y?，t）的關(guān)系是如同 h 為常數(shù)一樣的色散關(guān)系。根據(jù)已知每個(gè)網(wǎng)格點(diǎn)的水深，及波浪在此傳播過(guò)程中周期不變，可求得每個(gè)網(wǎng)格點(diǎn)的波數(shù) k。波數(shù)無(wú)旋性條件還可寫(xiě)成如下形式：

式 （15） 為一階的非線(xiàn)性偏微分方程，通過(guò)有限差分可以求得其數(shù)值解。

差分格式采用以下形式：

可用迭代法求解。w為權(quán)重因子，α為松弛因子，用以加速數(shù)值迭代之速度。當(dāng) w=1.0 時(shí)，方程為全隱格式。因?yàn)槭?（16） 中，x 方向?yàn)橄蚯安罘?，y方向?yàn)橹行牟罘?，引進(jìn)權(quán)重因子 w，考慮第 i行 y方向中心差分的影響。松弛因子 α的作用是考慮 j-1 列和 j+1 列的影響。波浪傳播方向的初值用 Snell定律求解，然后用迭代的方法求解式 （16），當(dāng)兩次求解的波向角相差小于某個(gè)小參數(shù)δ時(shí)，迭代停止。

下面為構(gòu)造方程式 （14） 的 Crank-Nicolson[10]格式離散方程：

k1（i，j）代表 k1在 （i，j）處的值，k2（i，j）代表 k2在 （i，j）處的值。式 （17） 在波數(shù) k 的離散和波幅的 y 方向做的 Crank-Nicolson[10]離散，x 方向是從外海邊界到海岸方向做步進(jìn)求解 （a0為波幅邊界條件），計(jì)算格式實(shí)現(xiàn)了穩(wěn)定，編程的工作量也相對(duì)減輕。

y方向透射邊界條件：

3 計(jì)算實(shí)例

荷蘭 Berkhoff[2]、Berkhoff等[3]的橢圓形淺灘試驗(yàn)是緩坡方程被引用最多的試驗(yàn)，試驗(yàn)地形和測(cè)量斷面如圖1所示。為了驗(yàn)證本文模型，根據(jù)上節(jié)中的方法，將本文模型應(yīng)用于上述算例。采用等間距網(wǎng)格，Δx 和 Δy 均取為 0.25m。

圖1 Berkhoff[2]、Berkhoff 等[3]橢圓形淺灘試驗(yàn)水深等值線(xiàn) （單位：m）及 8個(gè)波幅測(cè)量斷面分布Fig.1 Contoursofwater depth（m）and distribution of eight sections for wave am p litudemeasurements of Berkhoff[2]、Berkhoff etal.[3]elliptic shoal test

根據(jù)上節(jié)中波浪傳播方向的求解方法，計(jì)算得到的越過(guò) Berkhoff[2]、Berkhoff等[3]試驗(yàn)橢圓形淺灘的波向角計(jì)算結(jié)果如圖 2 所示。因?yàn)閷?duì) Berkhoff[2]、Berkhoff等[3]淺灘試驗(yàn)而言，波浪傳播方向的變化主要發(fā)生在淺灘之后，淺灘之前和靠近側(cè)邊界的地方波向基本沒(méi)有變化，所以，在圖2中，只選取了部分區(qū)域，比如在圖 2 （a） 中，取 5m

圖2 本模型計(jì)算得到的波浪越過(guò) Berkhoff[2]、Berkhoff等[3]試驗(yàn)中橢圓形淺灘的波向量分布Fig.2 Distributionsofwave vectorscalculated by the presentm odelover the elliptic shoalof Berkhoff[2]or Berkhoff etal.[3]’s test

因 為 波 浪 傳 播 方 向 很 難 測(cè) 量 ， Berkhoff[2]、Berkhoff 等[3]試驗(yàn)中 沒(méi)有提供波向角的試驗(yàn)數(shù)據(jù)。本文模型計(jì)算得到的波向角結(jié)果，以及其他模型的波向角結(jié)果，均不能通過(guò)試驗(yàn)數(shù)據(jù)驗(yàn)證精確性。鑒于以上原因，本文通過(guò)與其他模型的波向角計(jì)算結(jié)果比較，分析計(jì)算結(jié)果的合理性。Maa 等[11]針對(duì)應(yīng)用較為廣泛的幾種模型，如 REF/DIF-1、RCPWAVE 等，做了歸類(lèi)和比較，給出了各種不同模型關(guān)于波向角的計(jì)算結(jié)果?？偨Y(jié)各個(gè)模型波向角計(jì)算結(jié)果的特點(diǎn)，波向的改變主要發(fā)生在淺灘后，并且向淺灘后的中心線(xiàn)相聚，這與本文的結(jié)果圖 2 中 （b） 區(qū)域是符合的。

位于淺灘后的8個(gè)斷面的計(jì)算值和試驗(yàn)值的對(duì)比如圖 3。對(duì)于 Berkhoff[2]、Berkhoff等[3]的特定地形，本文計(jì)算值給出了一個(gè)較為合理的結(jié)果，波高分布與試驗(yàn)值較為相似。斷面SⅠ的計(jì)算值與試驗(yàn)值非常吻合；斷面SⅡ的計(jì)算值與試驗(yàn)值在淺灘兩側(cè)吻合的較好，在淺灘上計(jì)算值要高于試驗(yàn)值；斷面 SⅢ、斷面 SⅣ和斷面 SⅤ位于淺灘后，由于折射作用，會(huì)發(fā)生波能的相聚，波高增大，某些區(qū)域達(dá)到了入射波高的兩倍多，本文計(jì)算值給出了一個(gè)在波高和變化趨勢(shì)上非常接近的結(jié)果，但是沒(méi)有反映出兩側(cè)波高“減小—增大—減小”波動(dòng)的過(guò)程。在向岸方向的3個(gè)斷面中，斷面SⅥ計(jì)算值與試驗(yàn)值吻合程度不夠理想，斷面 SⅦ和斷面 SⅧ計(jì)算值與試驗(yàn)值吻合較好。

圖3 本模型計(jì)算值與 Berkhoff[2]、Berkhoff 等[3]橢圓形淺灘試驗(yàn) 8 個(gè)斷面（SⅠ、SⅡ、SⅢ、SⅣ、SⅤ、SⅥ、SⅦ和 SⅧ） 的相對(duì)波高試驗(yàn)值比較Fig.3 Com parison between calculated relativewave heightsby the presentmodeland experimental data for eight sections（SⅠ，SⅡ，SⅢ，SⅣ，SⅤ，SⅥ，SⅦ，and SⅧ） of Berkhoff[2]、Berkhoff etal.[3]'selliptic shoal test

本文數(shù)值結(jié)果的不足及其原因：

1） 如圖 3 所示，在剛過(guò)淺灘的地方，斷面SⅡ與斷面 SⅢ之間，計(jì)算值略高于試驗(yàn)值，這可能是由于本文模型為線(xiàn)性模型，沒(méi)有考慮非線(xiàn)性作用。如果考慮非線(xiàn)性作用，能夠得到更好的結(jié)果。在本模型以后的研究中，可以引入非線(xiàn)性色散關(guān)系來(lái)考慮非線(xiàn)性效應(yīng)。

2）與試驗(yàn)數(shù)據(jù)相比，在淺灘的后面，即入射邊界 10m 之后，斷面 SⅣ和斷面 SⅤ缺少因繞射而引起的波高在淺灘兩側(cè)的波動(dòng)，顯得平滑。原因是因?yàn)橐肓斯獬毯瘮?shù)方程。

本文的目的是對(duì) Liu 和 Shi[1]提出的 數(shù)學(xué)模式進(jìn)行數(shù)值計(jì)算驗(yàn)證和擴(kuò)展，使其可以應(yīng)用于平面二維波浪變形問(wèn)題，本文重點(diǎn)在于驗(yàn)證 Liu 和Shi[1]模型應(yīng)用與平面二維問(wèn)題的有效性。模擬的波高結(jié)果在部分?jǐn)嗝嫔系挠?jì)算精度尚需提高，今后，可以將計(jì)算結(jié)果與其它模型進(jìn)行定量對(duì)比。

4 結(jié)語(yǔ)

本文基于 Liu 和 Shi[1]的數(shù)學(xué)模型，引入繞射因子，利用有限差分法，對(duì)數(shù)學(xué)模型的波幅方程和邊界條件進(jìn)行數(shù)值離散，空間步進(jìn)求解波浪由外海向近岸傳播時(shí)，因?yàn)檎凵洹⒗@射和淺化等波浪變形現(xiàn)象，引起的波高、波向等波浪要素的變化。將本文的數(shù)值模 型應(yīng)用于經(jīng)典的 Berkhoff[2]、Berkhoff等[3]橢圓形淺灘算例，并與其試驗(yàn)結(jié)果做比較驗(yàn)證，證實(shí)本文模型求解二維地形下的波浪折射、繞射是可行的，具有方程求解簡(jiǎn)單、可適用于大面積海域的優(yōu)點(diǎn)。并且，光程函數(shù)的引入，使 Liu 和 Shi[1]的模型應(yīng)用范圍更廣泛，在求解平面二維問(wèn)題時(shí)精度得到提高。

[1]LIU Y Z，SHIJZ.A theoretical formulation for wave propagations overuneven bottoms[J].Ocean Engineering，2008，35（3-4）：426-432.

[2]BERKHOFF J C W.Computation of combined refractiondiffraction[C]//Proceedingsof the 13th international conference on coastalengineering.Vancouver，ASCE，1972：471-490.

[3]BERKHOFF J C W，BOOIJ N，RADDER A C.Verification of numerical wave propagation models for simple harmonic linear waterwaves[J].Coastal Engineering，1982（6）：255-279.

[4]ECKARTC.The propagation of gravitywaves from deep to shallow water[M]//Gravitywaves.Washington DC：National Bureau of StandardsCircular 521，1952：165-173.

[5]BIESEL F.Study ofwave progression in waterofgradually varying depth[M]//Gravity waves.Washington D C：National Bureau of Standards Circular521，1952：243-253.

[6]BATTJESJA.Refraction ofwaterwaves[J].Journalof theWaterways，Harbors and Coastal Engineering Division ASCE，1968，WW4：437-451.

[7]MEIC C，LEMEHAUTE B.Note on the equations of long waves over an uneven bottom[J].Journal of Geophysical Research，1968，71：393-400.

[8]MEIC C，TLAPA G A，EAGLESON PS.An asymptotic theory for water waves on beaches of mild slope[J].JournalofGeophysical Research，1968，73：4 555-4 560.

[9]LIU P L F，DINGEMANS M W.Derivation of the third-order evolution equations for weakly nonlinear water waves propagating over uneven bottoms[J].WaveMotion，1989，11：41-64.

[10]CRANK J，NICOLSON P.A practical method for numerical evaluation ofsolutions ofpartialdifferentialequationsof the heatconduction type[J].Proceedings of the Cambridge Philosophical Society，1947，43：50-67.

[11]MAA JP Y，HSU TW，TSAIC H，JUANG W J.Comparison of wave refraction and diffractionmodels[J].Journal ofCoastal Research，2000，16：1 073-1 082.

A two-dimensional horizontal numericalmodel of wave propagation over an uneven coastal bottom

PENGYan-jian，LIUYing-zhong， SHIJohn Z.
（SchoolofNavalArchitecture，Ocean and CivilEngineering，Shanghai Jiao Tong University，Shanghai200030，China）

Based on the linearwave theory，by using themultiple-scalemethod and a small slope parameter ε，Liu and Shi[1]obtain the first-order approximate solutions of the wave potential function，aswell aswave amplitude equation， for the problem of wave propagation over an uneven coastal bottom.This study introduces the diffraction factor，which expands Liu and Shi[1]＇s application range and improves itsaccuracy.Both wave amplitude equation and boundary conditions are discretized and solved by the finite differencemethod.Variations in wave amp litude and direction ofwave propagation caused by the refraction and diffraction are obtained.The validation study hasbeenmade useing of theexperimental results of Berkhoff[2]or Berkhoff et al[3]＇selliptic shoal test.It is suggested that the presentmodelbe applied to two-dimensionalwave propagation problems.

multiple-scalemethod；wave propagation；uneven bottom

TV139.2

A

1003-3688（2014）01-0008-06

10.7640/zggw js201401002

2013-04-12

2013-09-02

海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室自主研究課題 （GKZD010056-9）

彭延建 （1982 — ），男，山東德州市人，碩士，工程師，港口

海岸及近海工程專(zhuān)業(yè)。

* 通訊作者：時(shí)鐘，E-mail：zshi@sjtu.edu.cn