張俊生，滕 斌

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

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

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

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

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

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

1 數(shù)值模型的建立

1.1 控制方程

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

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

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

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

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

1.2 空間離散

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

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

各系數(shù)矩陣具體表達式:

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

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

1.3 邊界條件

1.3.1 入射邊界

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

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

1.3.2 全反射邊界

全反射邊界條件:

和

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

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

其逆變換可寫為:

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

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

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

當?shù)貎?yōu)勢資源容易得到充分挖掘。農(nóng)村當?shù)貎?yōu)勢資源主要是基礎(chǔ)性生產(chǎn)要素，如獨具特色的地質(zhì)地貌、歷史人物、宗教圣地、傳統(tǒng)工藝及民俗風情等自然資源與文人資源，農(nóng)民工返鄉(xiāng)創(chuàng)業(yè)集群通過有形和無形的合約推動技術(shù)進步、金融互助、生產(chǎn)銷售和服務(wù)聯(lián)合等，實現(xiàn)優(yōu)勢資源的有效整合，打造集群的競爭優(yōu)勢，從而以點帶面帶動整個區(qū)域產(chǎn)業(yè)集群全面發(fā)展。

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

1.3.3 吸收邊界

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

1.3.4 松弛區(qū)的使用

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

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

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

1.4 時間積分

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

若把時間積分方程表達為

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

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

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

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

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

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

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

后續(xù)算例驗證了該式的準確性。

2 數(shù)值模型的驗證

2.1 造波的驗證

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

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

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

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

圖2 非線性規(guī)則波計算結(jié)果與解析解比較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］的橢圓型淺灘上波浪傳播的試驗進行數(shù)值計算。計算域為長24.5 m、寬20 m 的矩形區(qū)域，模擬的海底地形為斜坡上隆起橢圓半球，地形具體表達式為:

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

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

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

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

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

2.3 曲邊界計算的驗證

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

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

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

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

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

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

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

圖7 圓柱周圍網(wǎng)格劃分示意Fig.7 Mesh around the cyliner

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

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

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

3 結(jié) 語

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

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

［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］孫忠濱，柳淑學(xué)，李金宣.基于改進Boussinesq 方程三角形網(wǎng)格有限元模型［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.