適用于曲邊界的淺水非線性波浪數值模型

張俊生，滕 斌

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

近岸及港口工程以及遠海淺灘上的人工島修建都必須考慮淺水波浪的作用和影響，而淺水非線性波浪數值模型則是其重要分析和計算工具。要做到計算準確、結果有效，必須滿足三個基本方面要求:1)能真實有效地模擬淺水非線性波浪;2)能在水深變化的大范圍水域內計算;3)能對曲邊界準確求解。這是淺水波浪特性、相關工程海域海底與邊界的現實需求。

Boussinesq 方程(以下簡稱為“B 方程”)是典型的淺水波浪數學模型，能夠準確描述淺水波浪的特性。Peregrine［1］推導出了適用于變水深三維問題的二維簡化B 方程。其二維簡化和變水深的特點使得B 方程成為分析淺水波浪的高效工具，能夠實現大面積水域內的波浪模擬。其后諸多學者不斷在色散性、非線性方面改進。Nwogu［2］和Beji 等［3］分別推導出了具有較高色散關系的二階方程。前者能保證方程線性化后計算的波浪相速度在O(μ4)的精度時仍與Airy 波相速度相匹配，后者的線性化色散關系能達到線性波浪理論色散關系的二階Padé 逼近，且方程形式簡單，利于有限元法的數值計算。Gobbi 等［4］和Madsen 等［5］諸多學者則提出了具有更高色散性、更強非線性的高階方程。

在數值方法上，有限差分法因求解簡單、程序實現方便被大量使用，Gobbi 和Madsen 便使用了有限差分法進行求解。但是，有限差分法對曲邊界的處理是非常困難的，而實際工程中曲邊界卻大量存在，因此對邊界有更好適應性的有限單元法得到了探討和應用。Li 等［6］、Walkley［7］和Zhao 等［8］采用了有限元方法求解方程。在采用等參單元的情況下，有限元法無法求解四階及以上導數，三階導數也只能利用中間變量降階為二階，然后利用格林公式轉化為空間導數只有一階的邊界積分求解，以上學者均如此處理。而若不采用等參單元，求解將非常困難。因此高階方程對有限元方法是不能適用的。但對于淺水問題而言，含有三階空間導數的Beji 改進型方程足以適用，且其形式簡單，便于有限元方法應用。

有限元法的優勢是對曲邊界的處理，但遺憾的是上述學者的計算卻并不理想。Li 等［6］、Zhao 等［8］和孫忠濱等［9］均基于Beji 改進型方程的有限元模型給出了圓柱面上淺水波浪爬高的計算結果，但與Isaacson［10］的試驗和解析計算結果相比，差別較大。有限元方法的優勢沒有得到體現。

文中建立的有限元模型利用類似于Engelman 等［11］提出的方法，在邊界上對方程進行了從x、y 方向到法線、切線方向的矢量分量變換，準確地利用了邊界條件求解。同時，利用Fenton［12］提出的非線性規則波浪解建立了準確的造波方法，從而構建了一個能夠準確求解的、適用于曲邊界的淺水非線性波浪模型。

1 數值模型的建立

1.1 控制方程

Beji 提出的改進型B 方程的形式:

其中，u = (u，v)為水深平均的水平速度矢量，η 為波面高度，d(x，y)為水深，g 為重力加速度，下標t 為時間導數，▽= (?/?x，?/?y)為二維梯度算子。Li 等［6］已驗證β = 1/5 時，方程的色散性最佳。

為求解式(2)中的三階空間導數，引入變量a = ?η/?x 和b = ?η/?y，得:

a 和b 的求解采用Li 等［6］的面積加權平均值方法，即

其中，Nei是圍繞i 節點的單元總數，Aij是圍繞i 節點所有單元中的第j 個單元的面積，(*)ij表示的是在第j個單元中求出的i 節點相關值，Ai是圍繞i 節點所有單元的面積和，(*)i為面積加權平均后求出的i 節點的值。1.3.2 節中將給出a 和b 在邊界上的求解方法。

1.2 空間離散

利用伽遼金法，對式(1)、(3)離散。在處理式(3)時，利用格林公式將含二階導數的面積積分轉化為只含一階導數的邊界積分。這樣，可用等參單元求解，在單元內任意點上有:

其中，f 為任意變量，fj為該變量在單元各節點上的值，Nj為形函數，K 為單元節點數。得離散方程:

各系數矩陣具體表達式:

為方便起見，若以x1、x2代替x、y，則Pij、Qij、Rij、Sij可統一表示為

其中，Pij:δ = 1，l = 1;Qij:δ = 1，l = 2;Rij:δ = 2，l = 2;Sij:δ = 2，l = 1。

1.3 邊界條件

1.3.1 入射邊界

文中采用Fenton［12］提出的非線性規則波浪解作為入射邊界條件，其具體表達式:

其中，Ψ = jk(xcosφ + ysinφ － ct)，k 為波數，c 為波速，φ 為波浪入射方向與x 方向夾角，U 為計算所得流速，Bj、Yj為兩組傅里葉展開系數，N 為展開多項式項數，文中取為20。k、c、U、Bj、Yj均可根據已知波浪參數由一套方程組求得，Fenton 給出了該方程組求解的方法。

1.3.2 全反射邊界

全反射邊界條件:

和

其中，n(nx，ny)為邊界外法向單位向量。該條件在計算中如何在曲邊界上得到滿足是計算準確性的關鍵。Li 等［6］在曲邊界上對所有的速度均取為零;Zhao 等［8］則在邊界上以式(11)直接替換動量方程;孫忠濱等［9］則根據邊界上法向矢量分量nx、ny大小的不同，以式(11)代替不同方向的動量方程，此法僅為避免nx、ny中趨于0 的量在求解方程組時做分母，本質上還是在曲邊界處直接求解x、y 方向的速度分量u、v，曲邊界處求解的準確性并沒有提高。

這里采用速度矢量分量變換的方法處理邊界上的計算。設直墻邊界處單位切向量為τ(τx，τy)，且滿足τ = k × n，對所求速度進行分量轉換，有:

其逆變換可寫為:

將式(14)代入式(7)、(8)，在方程組中直接求解直墻邊界上的un、uτ。對于這些邊界節點處的動量方程需做轉化，以式(11)直接代替x 方向動量方程，以切線方向動量方程代替y 方向動量方程。切線方向動量方程為:

其中，Mx、My為x、y 方向動量方程。

求解a = ?η/?x 和b = ?η/?y 時，在直墻邊界處同樣以求解?η/?n 和?η/?τ 替代之。如此，全反射邊界條件式(12)可直接應用，同時利用面積加權平均值求解?η/?τ:

當地優勢資源容易得到充分挖掘。農村當地優勢資源主要是基礎性生產要素，如獨具特色的地質地貌、歷史人物、宗教圣地、傳統工藝及民俗風情等自然資源與文人資源，農民工返鄉創業集群通過有形和無形的合約推動技術進步、金融互助、生產銷售和服務聯合等，實現優勢資源的有效整合，打造集群的競爭優勢，從而以點帶面帶動整個區域產業集群全面發展。

最后，在全反射邊界上的i 節點處可得:

1.3.3 吸收邊界

為避免波浪在開邊界處反射，在計算域的開邊界前布置一阻尼層吸收波浪。文中給出阻尼層為:

1.3.4 松弛區的使用

為了避免反射波浪對入射邊界的影響，采用Wang ＆Liu［13］類似的“松弛區方法”，如圖1 所示。此種情況下，入射邊界擴展為造波區，1.3.1 節中的入射條件解析式在整個造波區內給出。每時刻，松弛區內造波解析解與數值計算結果按照一定權重相加，以耗散掉反射到造波區的波浪，具體表達式:

其中，V 表示所計算的變量，下角標0 表示計算的有效結果，下角標I 表示造波解析解，下角標C 表示數值計算結果，ζ 為加權系數函數，本文采用式(18)。該方法可由后續算例證明行之有效。

圖1 松弛區布置示意Fig.1 Position of the relaxation zone

1.4 時間積分

采用四階Adams－Bashforth－Moulton 預報－校正方法進行時間積分計算。但為保證準確，前4 個時間步采用單步的四階Runge－Kutta 法進行計算。

若把時間積分方程表達為

其中，上角標(n)表示t = nΔt 時刻各變量所取或求得的值。預報出預估值后，進一步用四階的Adams－Moulton方法進行校正，具體表達式為

以所有計算節點的相鄰兩次校正結果的相對誤差的加權平均值作為結果收斂的判別標準，表達式:

其中，上角標(n+1)表示當前一次校正結果，上角標(n+1)'表示上一次校正結果。文中取Δf = 10－5為時間積分收斂的判別標準。

求解速度變量的傳統做法［6，8－9］為在方程(7)、(8)中分別求解u、v，求解其中的一個變量時，另一個則作為已知處理。于是，兩套方程組求得的u、v 需彼此反復迭代，直至結果收斂。這種解法效率不高，且易發散。在解方程時把方程(7)、(8)視做一個方程組同時求解u、v，利用1.3.2 節中對邊界的處理方法，直接求解邊界處法向速度和切向速度。方程組數值解法上，采用“重開始廣義極小殘量法”(GMRES(k)，k 為計算中重開始迭代步數，文中取k=100)，該方法效率高，計算準，收斂快。Brussino ＆ Sonnad［14］給出了GMRES 及其它一些線性方程組數值解法的運算過程和效率比較。u、v 的收斂則分別以滿足式(23)來判別。

計算時間積分時，CFL 條件需要滿足:

其中，c 為波速，Δlmin為最小單元尺度。

淺水非線性波浪波峰窄而尖、波谷平而寬，必須有足夠的網格才能捕捉到有效的波形，否則計算失去準確性。Ursell［15］的研究已明確指出參數Ur= H/(k2d3)是描述淺水波浪波形的唯一參數，是衡量非線性的標準，被稱為“Ursell 數”。為方便使用，從Ur= 0.05 到16.20，每隔0.05 取一個波形，共324 個，通過對Fenton解析解和本文數值結果進行分析，擬合出淺水規則波一個波長所需劃分單元數的經驗公式:

后續算例驗證了該式的準確性。

2 數值模型的驗證

2.1 造波的驗證

為驗證本文模型模擬非線性波浪的效果，現計算兩組、八個波形，并與Fenton 方法［12］的解析解進行對比。計算域為12 m×2 m，水深d = 0.08 m，在計算域的右側鋪設一個入射波長寬度的阻尼層，吸收波浪。具體波浪參數如表1 所示。

表1 造波測試的波浪參數Tab.1 Parameters for the test of wave generation

兩組算例的波高不同，波長對應、依次增加，形成了非線性不同的8 個規則波浪，與解析解的比較結果如圖2 所示，該圖給出的本文結果均為時間達到20 個入射波周期時(t = 20T)的計算波面。

由圖2 可以明顯地看到，不論非線性的大小，本模型的結果與解析解都十分吻合，且計算穩定，20 個周期后，算出的波面仍然光滑、準確。其中，圖2(h)的Ursell 數最大，此時波浪非線性已非常強，由圖2(h)可以看出波谷非常寬且平，貼近靜水面，但計算仍能準確地模擬出該波浪。由此可見，本模型的造波方法行之有效，數值計算準確，求解結果穩定。同時，亦可看到采用的阻尼層有顯著的消波效果。

圖2 非線性規則波計算結果與解析解比較Fig.2 Comparison between the present result and the analytical result for nonlinear waves

圖3 橢圓形淺灘上波浪傳播計算域Fig.3 Computational domain over the elliptical shoal

2.2 波浪在變化地形上傳播的驗證

為驗證本文模型適用于計算大面積水深變化水域的波浪模擬，對Berkhoff 等［16］的橢圓型淺灘上波浪傳播的試驗進行數值計算。計算域為長24.5 m、寬20 m 的矩形區域，模擬的海底地形為斜坡上隆起橢圓半球，地形具體表達式為:

1)不考慮橢圓半球，斜坡上的水深為

其中，x' = xcos20° － ysin20°，y' = xsin20° + ycos20° 。

詳細計算域如圖3 所示，圖中標出了x－y 與x' －y'坐標系、地形等高線以及造波區、松弛區和阻尼層的布置位置。入射波浪為波高HI= 0.046 4 m、周期T = 1 s 的規則波。

采取Walkley［7］的做法:共計算波浪傳播50 個周期，取第44、45 周期的結果進行分析，得每一點的波高H，求出相對波高H/HI。共取出x=1、3、5、7、9 m 和y = －2、0、2 m 八個斷面的結果與Berkhoff 等的試驗結果及Walkley 的數值計算結果進行比較，如圖4 所示。本文結果與另二者吻合較好，說明本文模型能準確地模擬出大面積區域地形變化對波浪傳播的影響，計算穩定，顯現出了波浪的折射、繞射效果。圖5 給出了t=45 s(45 個周期)時，波浪傳播的透視圖，圖中同時顯示了地形，可以很明顯地觀察出地形對波面的影響。

圖4 橢圓形淺灘上波浪傳播計算結果比較Fig.4 Comparison of wave heights in various sections

2.3 曲邊界計算的驗證

為驗證本文模型對曲邊界計算的準確性，現對波浪在圓柱上的爬高進行計算。Isaacson［10］給出了四組淺水波浪爬高試驗，是由兩個不同大小的圓柱和兩個不同非線性(Ursell 數不同)的波浪進行的不同組合，Isaacson 同時也給出了其演算的解析解。對于數值方法而言，此問題的計算就是對曲邊界計算能力的考驗。正如引言所說，絕大多數B 方程的數值模型采用了有限差分法，因而喪失了對曲邊界準確計算的能力，故回避了該問題。在為數不多的有限元方法中，Li 等［6］、Zhao 等［8］和孫忠濱等［9］進行了圓柱上波浪爬高的計算，但只對照Isaacson 試驗中的第四組給出了結果，且計算結果并不理想。

文中對Isaacson 的四組試驗都進行了計算，表2 給出四組試驗的參數。

圖5 橢圓形淺灘上波浪傳播透視圖Fig.5 Perspective of the wave elevation over the elliptical shoal

表2 圓柱面波浪爬高試驗參數Tab.2 Parameters for the test of wave run-up on a cylinder

因測試波浪為平行于x 方向入射，故為了提高計算效率，利用對稱性，只取半個計算域簡化計算。如圖6 所示，半個圓柱放置于一側邊墻處，計算域長10.5 m，寬4.0 m，造波區、松弛區和阻尼層分別置于計算域的兩側，且寬度皆為入射波的波長，θ 作為描述圓周的參數。

在分析計算圓柱面波浪爬高時，必須排除圓柱繞射波浪在圖6 所示的上側邊墻反射后形成的反射波對波浪爬高結果的影響。考慮到兩個測試波浪的波長(1.570 8 m 和1.208 3 m)，該計算域4 m 的寬度可以保證至少4 個周期的爬高結果是不受干擾的，因此，文中取前4 個周期的計算結果進行分析。在網格劃分方面，圓柱周圍需要加密，半圓周分為80 個單元，如圖7 所示。

圖6 圓柱面波浪爬高計算域Fig.6 Computational domain in the test of wave run up on a cylinder

圖7 圓柱周圍網格劃分示意Fig.7 Mesh around the cyliner

圖8 圓柱面上波浪爬高計算結果比較Fig.8 Comparison of results of wave run－up on a cylinder

圖9 波浪在圓柱周圍繞射透視圖Fig.9 Perspective of the wave diffraction around a cylinder

圖8 給出了本文計算結果、Isaacson 試驗和解析計算結果以及線性理論結果，可以很明顯地看出在四組試驗中本文結果與Isaacson 的計算結果都吻合地很好，且在第Ⅰ、Ⅱ、Ⅳ組試驗中二者與Isaacson 的試驗結果也很接近，只第Ⅲ組相差較大。但值得注意的是，第Ⅲ組的測試波浪與第Ⅳ組相同，但圓柱尺寸比第Ⅳ組小，因此，兩者的爬高曲線必然不一樣，最大爬高自然亦不一樣，但給出的第Ⅲ組試驗結果卻與第Ⅳ組相仿，最大爬高也基本一致。在第Ⅳ組試驗結果中，同時給出了Li 等、Zhao 等和孫忠濱等的計算結果，可以看到三者的結果與Isaacson 的試驗和計算結果相比，偏差較大，相反，本文結果卻吻合地很好，特別是與Isaacson 的解析計算結果吻合非常理想。這充分體現出本文對曲邊界的處理是準確和合理的。圖9 給出了本文計算波浪傳播的透視圖，可以明顯地看出波浪在圓柱周圍的繞射和在圓柱上的爬高，波面光滑、穩定，說明了本文模型的計算有很好的收斂性和穩定性。

3 結 語

基于Beji 的改進型方程，利用有限單元法建立了一個波浪數值模型。本模型的核心是建立了準確的曲邊界計算方法，并配以速度分量同時在一個方程組求解的處理方法，利用高效的“重開始廣義極小殘量法”準確地求解曲邊界的波浪計算問題，真正意義上體現出了有限單元法對曲邊界求解的計算優勢。

通過三個算例的模型驗證，顯現出了本文模型:1)能準確地模擬淺水非線性波浪;2)能在大面積水底地形變化的水域上，準確模擬波浪的折射、繞射;3)能對曲邊界進行準確的計算。可見，本文模型是一個能有效而準確地適用于曲邊界的淺水非線性波浪模型，能夠應用于具有不規則邊界的、海底變化的大面積水域波浪模擬等實際工程問題，也為橢圓余弦波與物體相互作用的進一步分析和研究工作提供了可靠的數值計算工具。

［1］PEREGRINE D H.Long waves on a beach［J］.Journal of Fluid Mechanics，1967，27:815－827.

［2］NWOGU O.Alternative form of Boussinesq equations for nearshore wave propagation［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，1993，119(6):618－638.

［3］BEJI S，NADAOKA K.A formal derivation and numerical modeling of the improved Boussinesq equations for varying depth［J］.Ocean Engineering，1996，23(8):691－704.

［4］GOBBI M F，KIRBY J T.A fully nonlinear Boussinesq model for surface waves.Part 2:Extension to O(kh)4［J］.Journal of Fluid Mechanics，2000，405:181－210.

［5］MADSEN P A，BINGHAM H B，SCH?FFER H A.Boussinesq－type formulations for fully nonlinear and extremely dispersive water waves:derivation and analysis［J］.Proceedings of the Royal Society of London，Series A，2002，459:1075－1104.

［6］LI Y S，LIU S X，YU Y X，et al.Numerical modeling of Boussinesq equations by finite element method ［J］.Coastal Engineering，1999，37:97－122.

［7］WALKLEY M.A numerical method for extended Boussinesq shallow－water wave equations ［D］.Leeds:The University of Leeds，1999:110－132.

［8］ZHAO M，TENG B，LIU S X.Numerical simulation of improved Boussinesq equations by a finite element method［J］.Journal of Hydrodynamics，Ser.B，2003(4):31－40.

［9］孫忠濱，柳淑學，李金宣.基于改進Boussinesq 方程三角形網格有限元模型［J］.海洋工程，2010，28(3):50－58.(SUN Zhongbin，LIU Shuxue，LI Jinxuan.Triangular element FEM model based on modified Boussinesq equations［J］.The Ocean Engineering，2010，28(3):50－58.(in Chinese))

［10］ISAACSON M de St Q.Wave run－up around large circular cylinder［J］.Journal of the Waterway，Port，Coastal and Ocean Engineering，ASCE，1978，104(1):69－79.

［11］ENGELMAN M S，SANI R L.The implementation of normal and/or tangential boundary conditions in finite element codes for incompressible fluid flow［J］.International Journal for Numerical Methods in Fluids，1982，2:225－238.

［12］FENTON J D.The numerical solution of steady water wave problems［J］.Computers ＆ Geosciences，1988，14:357－368.

［13］WANG B L，LIU H.Higher order boussinesq－type equations for water waves on uneven bottom［J］.Applied Mathematics and Mechanics，2005，26(6):774－784.

［14］BRUSSINO G，SONNAD V.A comparison of direct and preconditioned iterative techniques for sparse，unsymmetric systems of linear equations［J］.International Journal for Numerical Methods in Engineering，1989，28:801－815.

［15］URSELL F.The long－wave paradox in the theory of gravity waves［J］.Proceedings of the Cambridge Philosophical Society，1953，49:685－694.

［16］BERKHOFF J C W，BOOY N，RADDER A C.Verification of numerical wave propagation models for simple harmonic linear water waves［J］.Coastal Engineering，1982，6:255－279.