基于雙層Boussinesq方程的三維波浪數(shù)值模型及其驗(yàn)證

鄒文峰，王 平, ，劉忠波，房克照，孫家文，張寧川

(1. 大連理工大學(xué) 海岸及近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024; 2. 國(guó)家海洋環(huán)境監(jiān)測(cè)中心 國(guó)家環(huán)境保護(hù)海洋生態(tài)環(huán)境整治修復(fù)重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023; 3. 大連海事大學(xué) 交通運(yùn)輸工程學(xué)院，遼寧 大連 116026)

近幾十年來(lái)，Boussinesq水波理論得到了長(zhǎng)足發(fā)展，Madsen和S?rensen[2]、Nwogu[3]、Wei等[4]、鄒志利[5]、Madsen等[6]、Lynett和Liu[7]、Liu和Fang[8]都對(duì)其進(jìn)行了研究?；谏鲜龇匠痰哪Ｐ鸵脖怀晒τ糜谀M近岸區(qū)域的波浪運(yùn)動(dòng)，最新的Boussinesq水波方程的研究進(jìn)展可詳見(jiàn)Kirby[9]、張堯等[10]、孫家文等[11]的綜述。

傳統(tǒng)型式的Boussinesq水波方程將三維波浪問(wèn)題簡(jiǎn)化成二維問(wèn)題，較大程度上提高了數(shù)值模型的計(jì)算效率。當(dāng)前，基于Wei等[4]方程的FUNWAVE模型、基于Madsen和S?rensen[2]方程的MIKE21-BW模塊以及基于Lynett和Liu[7]方程的COULWAVE已得到較大程度的應(yīng)用。與這些國(guó)外流行的模型相比，國(guó)內(nèi)學(xué)者在一些基礎(chǔ)數(shù)值模型方面也取得了一些進(jìn)展，如柳淑學(xué)等[12]基于Beji和Nadaoka[13]的Boussinesq方程建立了有限元模型；劉忠波[14]、房克照[15]針對(duì)鄒志利[5]的Boussinesq水波方程建立了一系列有限差分?jǐn)?shù)值計(jì)算模型。不同于傳統(tǒng)型式的Boussinesq水波方程，F(xiàn)uhrman[16]、王本龍[17]基于Madsen等[6]的三維Boussinesq方程開(kāi)發(fā)了近岸波浪計(jì)算模型。

Liu等[1]給出了最高空間導(dǎo)數(shù)為2、3和5的多層Boussinesq水波方程，此多層方程具有最大的穩(wěn)定空間(Madsen和Fuhrman[18])。最近Liu和Fang[19]、Liu等[20]利用有限差分法開(kāi)發(fā)了最高空間導(dǎo)數(shù)為2和3的垂直二維模型，并將數(shù)值模型應(yīng)用到波浪與潛堤相互作用、雙色波群演化、非線性聚焦波演化以及海床運(yùn)動(dòng)興波問(wèn)題(Fang等[21])。

近期，Sun等[22]綜合分析了最高導(dǎo)數(shù)為2的雙層Boussinesq水波方程(Liu等[1])，發(fā)現(xiàn)方程中參數(shù)取值在0.13～0.25任意選擇時(shí)，水深在0

1 數(shù)值模型

1.1 雙層Boussinesq方程

Liu等[1]將水體分為兩層，分別考慮自由面處的變量η,uη,wη(z=η)、靜止水面處的變量u0,w0(z=0)、第一層中間位置處的變量uα,wα(z=zα)、第二層中間位置處的變量uβ,wβ(z=zβ)，水體分層及各變量如圖1所示。

圖1 雙層Boussinesq方程的坐標(biāo)示意Fig. 1 Schematic diagram of the two-layer Boussinesq model

在自由面上，滿足連續(xù)性條件及動(dòng)力學(xué)邊界條件，對(duì)應(yīng)的方程為：

(1)

(2)

式中：η為波面升高，uη和wη為自由水平面上的水平速度和垂向速度,uη=(uη，vη)，Ut為U在時(shí)間上的偏導(dǎo)，g為重力加速度，?為水平梯度算子。

自由面上的速度可利用靜水位處的速度來(lái)表達(dá)(Madsen等[6]；Liu和Fang[8])，以獲取更準(zhǔn)確的非線性特征，二者的關(guān)系式為：

(3)

(4)

式中：u0和w0表示定義在靜水位處的速度，可以通過(guò)將z=0代入第一層的速度表達(dá)式來(lái)確定。

結(jié)合兩層連接處的速度相等條件、水底運(yùn)動(dòng)學(xué)邊界條件、每層水體速度的泰勒展開(kāi)式，并引入偽速度來(lái)表達(dá)，得到：

u0=uα*-σ1?(?·uα*)+σ2?wα*-σ3?h(?·uα*)-σ4?h?2wα*

(5)

w0=wα*-σ1?2wα*-σ2?·uα*-σ3?h·?wα*+σ4?h·?(?·uα*)

(6)

(7)

(8)

(9)

方程(1)～(9)組成了雙層Boussinesq模型的控制方程，具體方程推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)[1]。1%色散誤差下，方程適應(yīng)無(wú)因次水深到kh=19.7，見(jiàn)圖2。當(dāng)α1=0.5時(shí)水體變?yōu)橐粚?，方?7)和(8)取消，上述雙層模型可退化為單層Boussinesq模型。

圖2 方程無(wú)因次相速度變化Fig. 2 Phase velocity of the two-layer Boussinesq model

1.2 數(shù)值求解

上述控制方程由η,uη,wη,u0,w0,uα*,wα,uβ*,wβ這9個(gè)變量和9個(gè)方程組成，方程的空間離散采用規(guī)則網(wǎng)格下的有限差分方法，在時(shí)間步內(nèi)采用預(yù)報(bào)—校正—迭代的形式進(jìn)行求解，直到計(jì)算誤差落入設(shè)定的誤差范圍內(nèi)，一次迭代過(guò)程分6步求解，數(shù)值計(jì)算流程詳見(jiàn)圖3。參數(shù)η,uη由方程(1)和(2)求解，在預(yù)報(bào)時(shí)采用三階Adams-Bashforth的時(shí)間格式，校正時(shí)采用四階Adams-Moulton格式，具體構(gòu)造為：

圖3 三維波浪模型數(shù)值計(jì)算流程Fig. 3 Numerical simulation process of the two-layer Boussinesq model

(10)

(11)

式中：上標(biāo)n為第n個(gè)時(shí)間步的標(biāo)記，即對(duì)應(yīng)n×Δt時(shí)刻的參數(shù)值；n+1對(duì)應(yīng)(n+1)×Δt時(shí)刻的參數(shù)值；其他依此類推。

針對(duì)u0,uα*,uβ*,wβ,wα則分別由式(3)、(5)、(7) 、(9)、(8)求解，將u0,uα*,uβ*,wβ,wα作為未知變量移至方程的左側(cè)，其余項(xiàng)移至方程右側(cè)作為已知量。在空間離散上，對(duì)于方程左端未知量的一次、二次導(dǎo)數(shù)采用式(12)、(14)的二階精度格式；對(duì)于方程右端已知量的一次導(dǎo)數(shù)采用四階精度式(13)求解，二次導(dǎo)數(shù)采用二階精度式(14)求解，具體格式為：

fx=(fi+1-fi-1)/(2Δx)+O(Δx)2

(12)

fx=(-fi+2+8fi+1-8fi-1+fi-2)/(12Δx)+O(Δx)4

(13)

fxx=(fi+1-2fi+fi-1)/(Δx)2+O(Δx)2

(14)

(15)

(16)

(17)

對(duì)于w0,wη由方程(6)和(4)直接求解，完成1～6步求解后即得到更新后的變量值，分析校正值和預(yù)報(bào)值誤差，若不滿足要求，則對(duì)預(yù)報(bào)值進(jìn)行更新，并重復(fù)1～6步計(jì)算調(diào)整后的校正值，直至其滿足誤差要求。

1.3 入射邊界與出口邊界

邊界造波法和內(nèi)部造波法是各類Boussinesq數(shù)值模型采用的兩種主要造波技術(shù)，邊界造波通?？梢赃x擇Stokes線性波理論、二階波理論或其他波浪理論等，給出邊界上的波面位移及速度值。

內(nèi)部造波法可避免在初始邊界采用固定造波板造波時(shí)，地形反射給造波板引起的二次反射問(wèn)題，模型采用了Hsiao等[23]的內(nèi)部造波法，其源項(xiàng)表達(dá)式為：

(18)

式中：Δt為數(shù)值模型中的時(shí)間步長(zhǎng)，T為入射波浪周期，Cg為波群速度，H為波高，C為波速，β=20/L2，L為波長(zhǎng)，x為網(wǎng)格點(diǎn)坐標(biāo)，x0為造波源點(diǎn)坐標(biāo)，ω為波浪角頻率，k為波數(shù)，t為數(shù)值模型當(dāng)前運(yùn)行時(shí)間。

對(duì)于出口邊界在本文所有的數(shù)值計(jì)算中，為提高數(shù)值模型穩(wěn)定性，數(shù)值水槽兩端均采用海綿層吸收條件，而兩側(cè)邊界采用反對(duì)稱的可滑移邊界條件；為使程序運(yùn)行穩(wěn)定，模型采用了Kirby等[24]給出的光滑濾波技術(shù)。

2 模型驗(yàn)證與分析

為了驗(yàn)證本文建立的模型，開(kāi)展了在深水規(guī)則波傳播演化、波浪在橢圓形淺灘地形的演化以及波浪在凹形淺灘上的演化三種情況的模擬研究，并與相關(guān)解析解和試驗(yàn)結(jié)果進(jìn)行了比較，驗(yàn)證模型的適應(yīng)性。

2.1 深水中規(guī)則波傳播過(guò)程

深水波浪傳播試驗(yàn)中波浪周期為 10 s，波高為 1.0 m，水深分別為 156 m和 468 m，其對(duì)應(yīng)的無(wú)因次水深kh分別為2π和6π。在數(shù)值模擬中采用邊界造波方式，空間網(wǎng)格大小為5 m，時(shí)間步長(zhǎng)為0.2 s。

圖4給出了在空間上模型和波面解析解的對(duì)比，圖5給出了x=1 000 m處波面歷時(shí)過(guò)程線對(duì)比，對(duì)比結(jié)果顯示，無(wú)論是kh=2π還是kh=6π，模型給出的波面演變結(jié)果均能與解析解較好地吻合，這反映出模型能夠精確再現(xiàn)規(guī)則波演化過(guò)程。

圖4 沿波浪傳播方向上的數(shù)值模擬與解析結(jié)果對(duì)比Fig. 4 Comparisons of the surface elevations between the numerical model and the analytical solution

圖5 x=1 000 m處波面過(guò)程數(shù)值模擬與解析結(jié)果對(duì)比Fig. 5 Comparisons of the surface elevations between the numerical model and the analytical solution at x=1 000 m

為考察三維模型對(duì)水平速度u及垂向速度w的模擬精度，針對(duì)kh=2π工況，圖6給出了x=1 000 m處波面位置的u和w，速度的模擬值與解析解吻合較好；圖7給出u和w的模擬值與解析解的對(duì)比，除第一層和第二層連接處附近略有差異外，其余深度的速度均能較好再現(xiàn)。針對(duì)kh=6π工況，盡管可以較好再現(xiàn)波面及波面處速度的演化過(guò)程，但由于其水深已超過(guò)本文雙層Boussinesq模型在速度分布特性上的適用范圍，其垂向速度累積誤差已超過(guò)6%，見(jiàn)圖8，這里不再給出對(duì)比結(jié)果。

圖6 x=1 000 m處波面速度u和w數(shù)值模擬與解析結(jié)果對(duì)比(kh=2π)Fig. 6 Comparisons of the horizontal and vertical velocities at the wave surface between the numerical model and the analytical solution (kh=2π)

圖7 x=1 000 m處沿水深速度u和w數(shù)值模擬與解析結(jié)果對(duì)比(kh=2π)Fig. 7 Comparisons of the horizontal and vertical velocities along the water depth between the numerical model and the analytical solution (kh=2π)

圖8 不同水深條件下的速度累積誤差情況Fig. 8 Accuracy of the horizontal and vertical velocity profiles

2.2 橢圓形淺灘上的波浪傳播變形

為研究波浪傳播過(guò)程中的折射、繞射、淺化以及非線性作用，Berkhoff等[25]設(shè)計(jì)了波浪在斜坡上疊加橢圓形淺灘的物理模型試驗(yàn)，相關(guān)試驗(yàn)數(shù)據(jù)已廣泛用于驗(yàn)證各種波浪數(shù)值模型，如基于N-S方程、各類緩坡方程和不同形式的Boussinesq方程的模型。

在數(shù)值模擬中，計(jì)算區(qū)域?yàn)?9 m×20 m，斜坡坡度為1∶50，波浪周期為1.0 s，波高為0.046 4 m，波浪正向入射，斜坡梯度方向與波浪入射方向的夾角為20°，內(nèi)部造波源位于x=4 m處的位置。斜坡旋轉(zhuǎn)后的坐標(biāo)與計(jì)算坐標(biāo)的關(guān)系為：

(19)

橢圓形淺灘中心坐標(biāo)為(0, 0)，淺灘邊界定義為：

(x1/3.0)2+(y1/4.0)2=1.0

(20)

平底區(qū)域及斜坡上的水深hs(單位：mm)為:

(21)

淺灘上水深h(單位：mm)為:

h=hs-0.5[1-(x1/3.75)2-(y1/5.0)2]0.5+0.3

(22)

數(shù)值模擬空間步長(zhǎng)為0.1 m×0.1 m，時(shí)間步長(zhǎng)為0.01 s，為避免過(guò)小水深引起數(shù)值異常，設(shè)置最小水深為0.07 m。分別采用了雙層和單層Boussinesq模型計(jì)算，得到各個(gè)斷面的計(jì)算結(jié)果與試驗(yàn)值對(duì)比見(jiàn)圖9。

從圖9可以看出，單層和雙層Boussinesq(BT)數(shù)值模型在8個(gè)斷面的計(jì)算結(jié)果與試驗(yàn)結(jié)果都基本吻合，在D3～D5斷面上模擬結(jié)果略小于試驗(yàn)數(shù)據(jù)，在Wei等[4]的數(shù)值研究中也出現(xiàn)類似現(xiàn)象，可能是由于強(qiáng)非線性作用，降低了波浪匯聚區(qū)的波幅值(Kirby和Dalrymple[26])。多層Boussinesq模型的非線性適用范圍更廣，其在D4～D5斷面上計(jì)算的波高值略小于單層模型。

圖9 數(shù)值結(jié)果與實(shí)測(cè)數(shù)據(jù)的比較Fig. 9 Comparisons of the wave heights between the numerical model and the measurement results

2.3 凹形淺灘上的波浪傳播變形

Whalian[27]設(shè)計(jì)了波浪在凹形淺灘上的非線性折射和繞射物理試驗(yàn)，在該地形上波浪會(huì)發(fā)生變淺匯聚作用，且淺水非線性作用產(chǎn)生了高次諧波過(guò)程，因此該試驗(yàn)被廣泛地用來(lái)驗(yàn)證各類波浪模型。試驗(yàn)水槽長(zhǎng)25.6 m，寬為6.096 m，試驗(yàn)的地形由下式給出：

(23)

其中，G=[y(6.096-y)]0.5，0≤y≤6.096。

數(shù)值模型中考慮不同周期及波高條件下的波浪演化過(guò)程，工況1～3的波浪周期為1 s、2 s、3 s，對(duì)應(yīng)的波高為0.019 4 m、0.015 0 m、0.013 6 m；工況4～6的波浪周期為1 s、2 s、3 s，對(duì)應(yīng)的波高為0.039 0 m、0.029 8 m、0.029 2 m。隨著入射波高的增大，波浪的非線性逐漸增強(qiáng)。

在數(shù)值模擬中，數(shù)值水槽設(shè)置長(zhǎng)35～45 m，采用內(nèi)部造波時(shí)，周期為1 s和2 s、3 s工況的源位置分別位于x=5 m和x=10 m處，時(shí)間步長(zhǎng)為0.02 s，空間步長(zhǎng)為0.1 m×0.101 6 m。圖10給出了雙層Boussinesq模型下工況1～3計(jì)算結(jié)果和試驗(yàn)結(jié)果的比較，圖11給出了工況4～6計(jì)算結(jié)果和試驗(yàn)結(jié)果對(duì)比，圖12給出了t=40 s時(shí)工況1～3的波面分布。

圖10 工況1～3計(jì)算波幅與實(shí)測(cè)值的比較Fig. 10 Comparisons of the computed and measurement results of wave amplitudes under case 1～3

圖11 工況4～6計(jì)算波幅與實(shí)測(cè)值的比較Fig. 11 Comparison of the computed and the measurement results of wave amplitudes under case 4～6

圖12 不同工況下的波面分布Fig. 12 Wave surface distribution under different cases

從圖10和11可見(jiàn)，對(duì)于波浪周期為1 s的工況，不同波高下的一次諧波和二次諧波的模擬結(jié)果與試驗(yàn)值吻合均較好。對(duì)于波浪周期為2 s工況，在非線性較弱的工況2中，一次諧波的模擬值試驗(yàn)值基本一致，而二次諧波和三次諧波的模擬值均略小于試驗(yàn)值；當(dāng)非線性增強(qiáng)時(shí)(工況5)，一次諧波和二次諧波的模擬值均略大于試驗(yàn)值，三次諧波的模擬值與試驗(yàn)值吻合較好。對(duì)于波浪周期為3 s工況，模型結(jié)果和試驗(yàn)結(jié)果差異較大，目前所見(jiàn)的模型結(jié)果多存在這一現(xiàn)象(Chen和Liu[28])?？傮w而言，模型非線性會(huì)對(duì)數(shù)值結(jié)果產(chǎn)生一定影響，本文建立的雙層Boussinesq模型對(duì)強(qiáng)非線性波浪的演化具有較好的模擬精度。

3 結(jié) 語(yǔ)

對(duì)Liu等[1]給出的最高導(dǎo)數(shù)為2的雙層Boussinesq水波方程在矩形網(wǎng)格上進(jìn)行了空間離散，時(shí)間積分上采用混合4階Adams-Bashforth-Moulton的預(yù)報(bào)—校正格式，得到了雙層Boussinesq方程下的三維波浪數(shù)值模型，模型具有較好的色散性和非線性。

針對(duì)水深在kh<10深水波浪，模型可以很好地給出波面、波面處的速度以及速度的垂向分布，具有較好的計(jì)算精度。針對(duì)橢圓形淺灘的物理試驗(yàn)，除波浪匯聚區(qū)的波幅值略小于實(shí)測(cè)值外，模型很好地刻畫(huà)了復(fù)雜地形上的波浪折射、繞射、淺化以及非線性作用。針對(duì)凹形淺灘地形上波浪淺水非線性作用產(chǎn)生的高次諧波過(guò)程，在周期1 s和2 s的工況中，模型均能很好地得到各次諧波值，且隨著非線性的增大，高次諧波更為顯著。通過(guò)不同典型算例的驗(yàn)證，本文建立的雙層Boussinesq模型對(duì)強(qiáng)非線性波浪的演化具有較好的模擬精度。