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

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

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

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

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

近期，Sun等[22]綜合分析了最高導數(shù)為2的雙層Boussinesq水波方程(Liu等[1])，發(fā)現(xiàn)方程中參數(shù)取值在0.13～0.25任意選擇時，水深在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方程的坐標示意Fig. 1 Schematic diagram of the two-layer Boussinesq model

在自由面上，滿足連續(xù)性條件及動力學邊界條件，對應的方程為：

(1)

(2)

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

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

(3)

(4)

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

結合兩層連接處的速度相等條件、水底運動學邊界條件、每層水體速度的泰勒展開式，并引入偽速度來表達，得到：

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模型的控制方程，具體方程推導過程見文獻[1]。1%色散誤差下，方程適應無因次水深到kh=19.7，見圖2。當α1=0.5時水體變?yōu)橐粚?，方?7)和(8)取消，上述雙層模型可退化為單層Boussinesq模型。

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

1.2 數(shù)值求解

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

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

(10)

(11)

式中：上標n為第n個時間步的標記，即對應n×Δt時刻的參數(shù)值；n+1對應(n+1)×Δt時刻的參數(shù)值；其他依此類推。

針對u0,uα*,uβ*,wβ,wα則分別由式(3)、(5)、(7) 、(9)、(8)求解，將u0,uα*,uβ*,wβ,wα作為未知變量移至方程的左側，其余項移至方程右側作為已知量。在空間離散上，對于方程左端未知量的一次、二次導數(shù)采用式(12)、(14)的二階精度格式；對于方程右端已知量的一次導數(shù)采用四階精度式(13)求解，二次導數(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)

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

1.3 入射邊界與出口邊界

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

內部造波法可避免在初始邊界采用固定造波板造波時，地形反射給造波板引起的二次反射問題，模型采用了Hsiao等[23]的內部造波法，其源項表達式為：

(18)

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

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

2 模型驗證與分析

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

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

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

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

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

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

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

圖6 x=1 000 m處波面速度u和w數(shù)值模擬與解析結果對比(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ù)值模擬與解析結果對比(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 橢圓形淺灘上的波浪傳播變形

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

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

(19)

橢圓形淺灘中心坐標為(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ù)值模擬空間步長為0.1 m×0.1 m，時間步長為0.01 s，為避免過小水深引起數(shù)值異常，設置最小水深為0.07 m。分別采用了雙層和單層Boussinesq模型計算，得到各個斷面的計算結果與試驗值對比見圖9。

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

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

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

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

(23)

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

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

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

圖10 工況1～3計算波幅與實測值的比較Fig. 10 Comparisons of the computed and measurement results of wave amplitudes under case 1～3

圖11 工況4～6計算波幅與實測值的比較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可見，對于波浪周期為1 s的工況，不同波高下的一次諧波和二次諧波的模擬結果與試驗值吻合均較好。對于波浪周期為2 s工況，在非線性較弱的工況2中，一次諧波的模擬值試驗值基本一致，而二次諧波和三次諧波的模擬值均略小于試驗值；當非線性增強時(工況5)，一次諧波和二次諧波的模擬值均略大于試驗值，三次諧波的模擬值與試驗值吻合較好。對于波浪周期為3 s工況，模型結果和試驗結果差異較大，目前所見的模型結果多存在這一現(xiàn)象(Chen和Liu[28])?？傮w而言，模型非線性會對數(shù)值結果產生一定影響，本文建立的雙層Boussinesq模型對強非線性波浪的演化具有較好的模擬精度。

3 結 語

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

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