基于四階完全非線性Boussinesq水波方程二維波浪傳播數值模型

房克照,鄒志利

(大連理工大學海岸和近海工程國家重點實驗室,遼寧 大連 116024)

基于四階完全非線性Boussinesq水波方程二維波浪傳播數值模型

房克照,鄒志利

(大連理工大學海岸和近海工程國家重點實驗室,遼寧 大連 116024)

建立基于四階完全非線性Boussinesq水波方程的二維波浪傳播數值模型。采用 Kennedy等提出的渦粘方法模擬波浪破碎。在矩形網格上對控制方程進行離散,采用高精度的數值格式對離散方程進行數值求解。對規則波在具有三維特征地形上的傳播過程進行了數值模擬,通過數值模擬結果與實驗結果的對比,對所建立的波浪傳播模型進行了驗證。同時,為了考察非線性對波浪傳播的影響,給出和上述模型具有同階色散性、變淺作用性能但僅具有二階完全非線性特征的波浪模型的數值結果。通過對比兩個模型的數值結果以及實驗數據,討論非線性在波浪傳播過程中的作用。研究結果表明,所建立的Boussinesq水波方程在深水范圍內不但具有較精確的色散性和變淺作用性能,而且具有四階完全非線性特征,適合模擬波浪在近岸水域的非線性運動。

Boussinesq水波方程;波浪傳播模型;波浪破碎;渦粘

Abstract:The present study develops a 2Dwave breakingmodel based on the fourth-order full nonlinear Boussinesq equations.This setof equations possesses higher accuracy in dispersion,shoaling and full nonlinearity characteristics,which enable it applicable in describingwater waves in nearshore zone.The discretizationis conducted on a rectangle grid system,higherorder finite difference formulasare used to approximate the spatial derivatives and high accuracy scheme is used to integrate equations in time.Kennedy et al′s eddy viscosity breakingmechanism is incorporated to mimicwave breaking.Themathematicalmodel is run to simulate waves transformation over typical three dimensional bathymetries.The numerical results are verified against available experimental data.The satisfactory agreement between numerical resultsand experimental data demonstrate the validation and efficiency of the presentwave breakingmodel.The effectof full nonlinearityonwave propagation is discussed by comparing the numerical resultswith those from Boussinesqmodelwith only the second-order full nonlinearity characteristics.

Key words:Boussinesqwave equations;wave propagationmodel;wave breaking;eddy viscosity

近幾十年來,基于Boussinesq水波方程的波浪模型得到了長足的發展并被廣泛用于模擬近岸區域各種波浪運動現象,見綜述性文獻[1]和[2]。建立這類波浪模型時,Boussiensq方程的性能是首要考慮的因素,方程的性能包括線性性能(色散性和變淺作用性能)和非線性性能(各階非線性傳遞函數以及波幅離散性能)。Boussinesq水波方程的發展過程很大程度上是提高上述各種性能的過程,起初諸多研究著重于提高方程的線性性能,代表性的工作有Madsen和Sorensen[3]、Nwogu[4]等,而后來的研究表明非線性在Boussinesq模型中有非常重要的作用。目前,具有強非線性特征的波浪模型不斷涌現,其中Madsen等[5]給出的方程以及Lynett和Liu[6]提出的多層模型頗具代表性。前者最大應用水深達到kh=40,但方程的數值求解過程涉及到大型稀疏矩陣的求解[7-8],即便對于一維問題控制方程多達6個,形成的稀疏矩陣帶寬至少為25[9];后者通過適當增加模型層數,也可以得到性能較好的方程,但隨著層數增加,控制方程增多,加大了數值求解工作的難度。相對而言,具有較低色散性精度但具有強非線性特征的方程數值求解簡便,仍不失為進行波浪數值模擬的有力工具。Wei等[10]給出了具有二階完全非線性特征的Boussinesq方程,基于此方程發展的二維模型FUNWAVE[11]得到了廣泛的應用。Madsen和Schaffer[12]給出了具有四階非線性特征的Boussinesq方程,其非線性達O(εμ4)階(這里ε為表征方程非線性的參數,定義為ε=A/h,A和h分別為特征波幅、特征水深;μ=kh為色散性參數,k為特征波數);Gobbi和Kirby[13]給出了具有四階完全非線性特征的Boussinesq方程FN4;Zou和Fang[14]也推導了具有四階完全非線性的Boussinesq方程BouN4D4。上述模型的相關研究成果均表明,具有強非線性特征的方程較僅具有弱非線性特征的模型有優越性,更適合模擬波浪的非線性運動。因此,建立具有強非線性特征的波浪模型非常必要,然而上述模型中,FUNWAVE雖然具有二階完全非線性特征,但方程色散性僅為精確色散關系的Pade[2,2]階近似,限制了模型在深水范圍的應用;Madsen和Schaffer[12]給出的方程僅具有弱非線性特征同時并未建立相應的數值模型;上述兩個具有完全非線性特征的方程(FN4和BouN4D4)雖然具有強非線性特征,但所建立的模型僅限于一維問題,也并未考慮波浪的破碎問題,而具有四階完全非線性特征的二維波浪模型更加鮮見。建立基于BouN4D4的二維波浪模型,該方程除具有四階完全非線性特征外,其色散性為精確色散關系的Pade[4,4]階近似,適用于深水波浪問題,同時具有良好的變淺作用性能。在矩形網格上對控制方程進行離散,建立高精度數值格式,考慮波浪的破碎問題。針對幾個典型三維地形上波浪進行數值模擬,通過對比分析數值結果和實驗數據,考察非線性特征在波浪傳播過程中的作用,同時驗證具有四階完全非線性特征波浪傳播模型的有效性。

1 數值模型

1.1 Boussinesq方程

Zou和Fang[14]推導了具有四階完全非線性特征的Boussinesq水波方程,該組方程以波面η和水深平均速度ˉu表達,經擴展后(考慮波浪破碎、水低摩擦等效應)方程的二維無因次形式為

式中:ˉu為水深平均速度,η為波面升高,h為靜水深,d=h+εη為當地水深,g為重力加速度,▽=(?x,?y)為二維偏微分算子 ,α1=1/9,β1=1/945,α2=0.146 488,β2=0.001 996 為常數 ,其中α1、β1確定了方程的色散性,α2、β2決定了方程的變淺作用性能。上式中R定義為R=Rf+Rb+Rs+Rm,其中Rf代表底摩擦阻力,Rb、Rs、Rm分別為波浪破碎引起的耗散項、海綿吸收層和側混項,其表達式同FUNWAVE中一致,fs為內域波浪生成時引入的源函數,同Gobbi和 Kirby[13]給出的一致。

后文中將上述模型稱為BouN4D4,N4代表4階完全非線性,D4代表色散關系為精確色散關系的Pade[4,4]階近似。該方程保留了精確到O(μ4)階的所有非線性項,最高階達O(ε5μ4),使得方程具有四階完全非線性特征,更適合描述具有非線性特征的波浪運動現象,方程的推導以及性能分析過程見文獻[14]。忽略上述方程中所有O(μ4)階非線性,則得到僅具有二階完全非線性特征的方程,稱之為BouN2D4(N2代表2階完全非線性,D4代表色散關系為精確色散關系的Pade[4,4]階近似),其將用于和BouN4D4進行對比,以考察非線性在波浪傳播過程中的作用。

1.2 數值求解

有限差分方法由于具有直觀性和容易程序實現的優點,被廣泛用來數值求解Boussinesq方程,發展得比較完善,如FUNWAVE[11]以及 Gobbi和 Kirby[13]。因此,這里仍然采用有限差分方法對上述方程進行數值求解,但值得注意的是若采用傳統的手工編程方法對方程(1)～(7))進行程序實現非常繁瑣,而通過充分利用MAPLE強大的數學公式推導功能及其向Fortran語言的轉化功能,可實現快速、準確地撰寫程序代碼的目的。文中采用的數值求解格式同上述文獻基本類似,這里不再一一詳述,詳細過程可參考上述相關文獻,下面僅給出主要步驟。

1.2.1 方程整理

將方程中含時間導數的線性項移在方程左側,其它項移至方程右側:

注意,上式中ξ,U,V中包含了所有含ηt,,的線性項。E,Et,F,Ft以及G,Gt分別為質量方程,動量方程中剩余的所有空間項和時間項,r項為式(2)中R在x,y方向的分量。

1.2.2 方程離散及差分公式

在圖1所示的矩形網格上對控制方程進行空間離散,x,y方向網格標記分別為i(i=1,mi)和j(j=1,mj)。對于方程右端空間導數采用七點差分格式,內點采用中心格式,邊界點采用偏心格式,對于單方向的差分統一用如下公式表示

式中:δfl表示函數在f網格點l處的空間導數(i-3≤l≤i+3)。當l=i時為中心差分,當l≠i時為偏心差分,γ,αm為差分格式的系數,其具體取值以及各階差分格式精度參見附錄A.1。對于水平二維問題,還涉及到交叉導數項,其具體實施辦法在附錄A.2中列出。可以看出,對于任意一個內點(i,j),其差分近似涉及到周圍25個點,見圖1。對于方程左端線性時間項的系數采用五點中心差分格式,則按上述網格離散后,形成一個系數矩陣帶寬為五的線性方程組。對于方程右端對時間的一階導數項,采用 Gobbi和 Kirby[13]文中公式進行計算,這里不再給出。

圖1 空間離散示意以及差分模板Fig.1 The sketch of spatial discretization and difference stencil

1.2.3 時間積分、邊界處理等

按照上述網格以及差分公式對方程離散后,采用五階Adams-Bashforth預報,六階Adams-Moulton校正格式進行時間步進求解。對于波浪入射邊界、固壁邊界的處理同用 Gobbi和 Kirby[13]文中完全一致(注意文獻中為一維),不再給出。波浪傳播過程中,由于非線性或者地形突變容易引起高頻數值震蕩,計算中采用SG光滑技術[7]對數值結果進行光滑處理。

2 模型驗證和對比

將對上面建立的數值模型進行驗證,針對幾個典型三維地形進行波浪傳播的數值模擬。通過對比數值模擬結果和實驗結果對模型進行驗證;對比模型BouN4D4和BouN2D4的數值結果,討論非線性作用在波浪傳播過程中的作用。

2.1 圓型淺灘上波浪的傳播

Chawla和Kirby[15]進行了波浪跨通過圓型淺灘傳播的實驗(見圖2)。該實驗在水池中進行(長18 m,寬18.2m),水池中央放置圓型淺灘,其半徑為2.57m,淺灘中心位于(x=8.0 m,y=8.98 m)處,在七個不同斷面上進行波高測量。圖2為計算區域以及測量斷面示意圖,取計算域長23m,寬18.2 m,在計算域兩端分別設置2m、3m的海綿層,兩側面為固體邊界,空間網格步長為0.10 m,時間步長0.02 s,每0.1個周期對波面和速度應用一次SG光滑器。這里僅選取其中兩組規則波實驗進行模擬。其中,第一組波浪沒有發生破碎,其對應的波浪要素為水深0.45 m(圓型淺灘頂部水深0.08m),波高0.011 8m,周期1 s;第二組存在波浪破碎現象,對應的波浪要素為水深0.395m,波高0.02m,周期1 s。該實驗中除淺灘外均為常水深,這將導致波浪在淺灘后產生較強的波浪幅聚現象。

對波浪不破碎情況,模型運行40 s,取最后6 s的數據進行波幅計算(30 s以后波浪場基本達到穩定狀態)并將結果在圖3中給出。可以看出,BouN4D4和BouN2D4的數值結果基本重合,同時和實驗數據吻合良好,這表明這一算例對高階非線性的要求較低。沿斷面A-A,模型較準確地模擬了由于波浪淺化引起的波高增大現象以及淺灘后波高減小現象。還可以看出,淺灘處波高幾乎高達入射波高的2.7倍,這進一步表明該處有強烈的波浪匯聚現象。對于其它的測量斷面,模型也給出了較好的模擬結果,這表明其較準確地模擬了波浪的折射、繞射過程。從圖中還可以看出,由于淺灘的位置并未處于計算域中心位置,計算結果和實驗數據關于淺灘中心軸線并不完全對稱。

對于波浪破碎情況,模型運行50 s,取最后10 s數據進行均方根波高計算,圖4給出了數值結果同實驗數據的對比(圖中H0為入射波高)。可以看出,BouN4D4和BouN2D4給出的數值結果相近,但較波浪不破碎情況差異加大,這可能是由于入射波浪非線性增強所致。在斷面E-E之前,BouN4D4在y=10m區域附近給出的幅值略高于BouN2D4并且更接近實驗數據;在淺灘頂部(斷面F-F),兩模型模擬的波高均高于實驗值但BouN4D4較BouN2D4和實驗數據吻合程度稍好。總體而言,兩個模型的數值表現良好,但較不破碎情況差。這可能是由于以下幾個原因:1)用于模擬波浪破碎的渦粘公式是經驗性的,不能很好地描述波浪破碎這一復雜物理現象,同時由于實驗中波浪破碎時具有三維特征,而采用的渦粘方法在數值實現時僅能沿某幾個固定方向對波浪破碎進行捕捉,存在誤差;2)波浪破碎時波面會變得非常陡峭,因此為了更好地刻畫波面,計算時所用的空間、時間步長要比較小,但出于計算時間和對比考慮,這里參數的取值同第一個算例一致。

圖5給出了針對這兩個算例t=35 s時的波面瞬時圖。可以看出,由于波浪的淺化、折射、繞射以及非線性的相互作用,最終的流場呈現復雜的三維形態。其中,波浪不破碎和破碎兩種情況存在以下共同特征:波峰在淺灘上發生了彎曲,這是由于波浪折射所致,同時由于波浪淺化作用波高增大;在淺灘末端出現了明顯的波浪匯聚現象;由于淺灘的存在其后部都存在波浪掩護區;觀察波峰線,都明顯存在波浪繞射現象。不同的是:波浪破碎情況下,淺灘上波面極其陡峭并且最大值出現在淺灘上,而波浪不破碎時淺灘上波面變化較緩且最大值出現在淺灘末端。

圖2 Chawla和kirby實驗布置以及測量斷面分布Fig.2 Experiment setup and measurement transects(Chawla and Kirby)

圖3 各測量斷面波高分布Fig.3 Wave height comparisons along different transects

圖4 各測量斷面波高分布Fig.4 Wave height comparisons along different transects

圖5 穩態后流場波高分布示意(BouN4D4)Fig.5 Elevation contoursof the steadywave fields(BouN4D4)

圖6給出了流場的時均流速圖。可以看到,在潛堤頂端,由于波浪破碎引起強烈的水流,幅值明顯大于其它區域波浪運動的水流速度。這同一般的波浪破碎理論相吻合,波浪破碎總是伴隨著流的產生,實驗過程中也觀測到了這一現象[15]。

2.2 凸型淺灘上波浪的傳播

Whalin[16]在波浪水槽中進行了一系列波浪在淺灘上的非線性折射和繞射物理實驗,實驗水槽尺度為25.603m×6.096m。實驗地形由以下公式給出式中:所有變量單位為米。圖7為計算域以及地形示意圖。在圖示地形上,波浪會發生變淺幅聚作用,由于非線性作用高次諧波變得尤為突出。因此,該實驗被廣泛地用來驗證各類波浪模型[7-8,18]。實驗可按照波浪周期不同分為三組,各組波浪要素及計算參數取值見表1,從表中可以看出,對應于某一周期,入射波浪的非線性依次增強,因此通過對這些工況的模擬,可以考察模型非線性精度對數值模擬結果的影響。下面,將模型BouN4D4和BouN2D4分別用于模擬Whalin的實驗并進行對比分析。

圖6 流場時均流速分布Fig.6 Phase-averaged velocity field

圖7 計算域及地形示意Fig.7 The computation domain and bathymetry sketch

表1 Whalin實驗波浪要素及數值模擬網格尺寸Tab.1 Wave conditions for Whalin experiment and grid size in simulations

圖8給出了T=1 s時,模型BouN4D4和BouN2D4的計算結果同實驗數據的對比。可以看出,數值結果同實驗數據吻合良好。其中,當H=0.019 4 m時,兩個模型的數值結果幾乎沒有差別,都較好地模擬了各次諧波沿空間的變化,但隨著波浪非線性的增強,兩個模型的差異變大,即當H=0.039 m時,BouN2D4給出的數值結果在x大于15 m時高于BouN4D4的數值結果,而BouN4D4和實驗數據更為接近。由于兩模型的色散性相同因此這一差別反映了非線性對數值結果的影響。圖9給出了T=2 s時,模型BouN4D4和BouN2D4的計算結果同實驗數據的對比。可以看出,隨著波浪非線性增強模型表現出的差異增大,這與上一算例一致,不同的是,兩個模型的數值結果差別較大。具體的,BouN2D4模擬的一、二次諧波幅值均高于實驗數據和BouN4D4的結果而后者同實驗數據更為吻合,其中一次諧波差別最大。對于三次諧波,BouN2D4給出的結果同實驗數據吻合良好而BouN4D4模擬的幅值略低于實驗數據,但考慮到三次諧波的幅值非常小(毫米量級),這一不足可以忽略。由于兩個模型的色散性以及變淺性能相同,因此這一差別反映了非線性對數值結果的影響。圖9也給出了T=3 s時,模型BouN4D4和BouN2D4的計算結果同實驗數據的對比。可以看出,模型結果同實驗數據差別較大,目前所見各類模型普遍存在這一現象[7-8,18]。不同于上述兩個算例,兩模型的數值結果幾乎沒有差別,這是由于該算例值最小,兩模型的非線性性能曲線在此處差別較小所致。因此,對于Whalin的實驗,當周期較小,非線性較大時,模型的非線性性能對數值結果的影響較大,總體而言,BouN4D4模擬的數值結果較BouN2D4更接近實驗數據。

圖8T=1.0 s時諧波幅值對比Fig.8 Comparisonsof harmonic amplitude forwave periodT=1.0 s

圖9 不同時刻諧波幅值對比Fig.9 Comparisonsof harmonic amplitude forwave periodT=2.0 s(left)andT=3.0 s(right)

3 結 語

建立基于四階完全非線性Boussinesq水波方程的二維波浪破碎模型。選用的控制方程BouN4D4具有強非線性特征同時具有較好的色散性和變淺作用性能,更適合于描述波浪的非線性運動。在矩形網格上對上述方程進行空間離散,采用高階導數近似方程中的時、空間導數,用高精度ABM預報校正格式進行時間積分,建立了高精度的數值求解格式,同時通過采用Kennedy等的渦粘破碎方法處理波浪破碎問題。

將所建立的具有四階完全非線性特征的波浪傳播模型BouN4D4和僅具有二階完全非線性特征的模型BouN2D4用于模擬規則波在具有三維特征地形上的傳播過程。通過對比兩個模型的數值結果以及相關的實驗數據,驗證了所建立模型的有效性,討論分析了非線性在波浪傳播過程中的作用。結果表明:1)實驗數據同數值結果吻合良好,這表明本研究所建立的模型可以較準確的模擬波浪在三維地形上的傳播、破碎過程;2)入射波浪非線性較強時,模型的強非線性特征發揮較大作用,具有抑制波幅過度增長的作用。

[1] Kirby J T.Boussinesqmodels and application to nearshorewave propagation,surfzone processes and wave-induced currents[M]∥advances in Coastal Engineering.Lakhan V C(ed),Elsevier,2003:1-41.

[2] 李孟國,王正林,蔣德才.近岸波浪傳播變形數學模型的研究與進展[J].海洋工程,2002,20(4):43-57.

[3] Madsen PA,Sφrensen O R.A new form of the Boussinesq equationswith improved linear dispersion characteristics.Part 2:A slowlyvarying bathymetry[J].Coast.Engrg.,1992,18:183-204.

[4] Nwogu O.An alternative form of the Boussinesq equations for nearshorewave propagation[J].J.Wtrwy,Port,Coast.and Oc.Engrg.,1993,119:618-638.

[5] Madsen PA,Bingham H B,Schaffer H A.Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves:derivation and analysis[J].Proc.R.Soc.London.A,2003,459(2033):1075-1104.

[6] Lynett P,Liu PL F.A two layer approach to wavemodelling[J].Proc.R.Soc.London.A,2004,460:2637-2669.

[7] Fuhramn D.Numerical solutionsof Boussinesq equations for fully nonlinear and extremely dispersivewaterwaves[D].Denmark:Technical Univ.of Denmark,2006.

[8] 王本龍.基于高階Boussinesq方程的海岸破波帶數學模型研究[D].上海:上海交通大學,2005.

[9] 張洪生,馮文靜,王亞玲,等.非線性波傳播的新型數值模擬模型及其實驗驗證[J].海洋學報,2007,29(4):137-147.

[10] Wei G,Kirby J T,Grilli ST.A fully nonlinear Boussinesqmodel for surfacewaves.Part I.Highly nonlinear unsteadywaves[J].J.Fluid Mech.,1995,294:71-92.

[11] Kirby J T,Wei G,Chen Q,et al.FUNWAVE 1.0,Fully nonlinear Boussinesq wavemodel documentation and user′smanual[R].CACR search report,University of Delaware,1998:CACR-98-06.

[12] Madsen PA,Sch?ffer H A.Higher-order Boussinesq-type formulation for surface gravity waves:Derivation and analysis[J].Phil.Trans.R.Soc.Lond.A,1998,356:3123-3184.

[13] Gobbi F,Kirby J T.Wave evolution over submerged sills:testsof a high-order Boussinesqmodel[J].Coast.Engrg.,1999,37(1):57-96.

[14] Zou Z L,Fang K Z.Alternative forms of the higher-order Boussinesq equations:Derivations and validations[J].Coast.Engrg.,2008,55(6):506-521.

[15] Chawla A,Kirby J T.Wave transformationover a submerged shoal[R].CACR research report,University of Delaware,1996:CACR-96-03.

[16] Whalin RW.the limitof applicability of linearwave refraction theory on a convergence zone[R].Res.Rep.,Vicksburg,MS:U.S.Army Corpsof Engineers,Waterways Expt.Station.1971:Res.Rep.H-71.

[17] Kennedy A B,Chen Q,Kirby J T,et al.Boussinesqmodeling of wave transformation,breaking,and runup[J].I:1D.J.Wtrwy,Port,Coast.and Oc.Engrg.,2000,47:39-47.

[18] Chan Y,Liu PL F.Modified Boussinesq equations and associated parabolicmodels forwaterwave propagation[J].J.Fluid Mech.,1995,288:351-381.

附錄:

1 單方向差分格式

導數階數/m 所在點 γ αi?3 αi?2 αi?1 αi αi±1 αi±2 αi±3i90 1i?1i?2i?3±60-12-98-10-147-24-77 360-45-35 150-450 80-100 400 45-30 50-225-15 72 1-12-10 2 3 4ii?1i?2i?3i i?1i?2i?3i i?1i?2i?3 180±8 6 2-12 137 812 1-1-15-49-151 7 35-27 228-147-3 132-8-85 6 232 12-18-84-186 270-420-255 5 265 13 35-83-461-39 21 171 411-490 200 470-5 080 0-48 64 496 56-4-184-484 270 15-285 2 970-13 29-29-307-39-9 111 321-27-12 93-972 8-88 104 12 6-36-114 22-13 137-11-1-15-1-151 7i4-505 5i?1i?2i?3±2-1-3-5-7 16 28 40-35-65-95 40 80 120-25-55-85-482 0 32 1-1-3-5

2 混合導數差分格式

混合導數差分格式可以表示為

式中:m和n分別表示沿x,y方向導數的階數;Fm,n為系數矩陣,分別為

對于F1,2,F2,3,F1,3,F1,4則可通過對上面公式的簡單變換得到。

A 2D numericalwavemodel based on fourth order fully nonlinear Boussinesq equations

FANG Ke-zhao,ZOU Zhi-li
(The State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116024,China)

TV148

A

1005-9865(2011)01-0032-09

2010-04-07

中央高校基本科研業務費專項資金資助項目;大連理工大學海岸和近海工程國家重點實驗室、河海大學海岸災害及防護教育部重點實驗室開放基金資助項目;國家自然科學基金資助項目(51009018)

房克照(1980-),男,山東日照人,講師,博士,從事海岸工程研究。E-mail:kfang@dlut.edu.cn