李巧生，唐軍，呂義港

( 1. 蘇州市航道管理處，江蘇 蘇州 215008；2. 大連理工大學(xué) 海岸和近海工程國家重點實驗室，遼寧 大連 116023；3. 嘉興市交通工程質(zhì)量安全管理服務(wù)中心，浙江 嘉興 314001)

1 引言

波浪是近岸環(huán)境中最重要的水動力因素之一，影響到水工建筑物的規(guī)劃和設(shè)計，也是近岸泥沙運動、岸線演變的主要動力。波浪從外海向近岸傳播過程中，受到地形、岸界等影響，發(fā)生折射、繞射、反射等變形。近岸波浪數(shù)值模擬是分析波浪傳播變形特性的主要手段，也是海岸工程水動力分析的理論技術(shù)基礎(chǔ)。緩坡方程波浪模型[1]綜合考慮了波浪傳播過程中的反射、折射和繞射等效應(yīng)，可有效地模擬近岸緩坡區(qū)域波浪的傳播變形，被廣泛應(yīng)用于近岸波浪場的模擬分析[2-6]。

近岸波浪數(shù)值模擬時，通過在規(guī)則矩形網(wǎng)格下對波浪模型離散求解。在天然近岸水域，受海岸建筑物（碼頭、防波堤等）、島嶼、不規(guī)則岸線等影響，波浪運動的水域通常比較復(fù)雜、不規(guī)則，采用矩形網(wǎng)格擬合必然會造成計算域邊界與實際邊界吻合較差，影響到數(shù)值模擬精度。為適應(yīng)工程實際中常有的不規(guī)則曲折邊界，研究者們在數(shù)值模型中引入了曲線網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格，但其求解的方程組不是對角占優(yōu)，網(wǎng)格偏離正交時會產(chǎn)生復(fù)雜的交叉項，所需內(nèi)存較大，計算時間長。因此，引入一種對復(fù)雜計算邊界適應(yīng)性好，并且計算過程簡單、穩(wěn)定收斂的數(shù)值方法是很有必要的。

近年來，非結(jié)構(gòu)化計算網(wǎng)格和四叉樹計算網(wǎng)格因?qū)?fù)雜計算邊界適應(yīng)性好，逐漸被用于水動力數(shù)值模擬[3-4,6]。Guo 等[7]基于非結(jié)構(gòu)化網(wǎng)格下的有限體積法，建立了河口潮流數(shù)值模型，模擬了錢塘江口的潮流及波浪變化，取得了較好的效果；林偉波和王義剛[8]采用非結(jié)構(gòu)化網(wǎng)格的海洋模型，建立了甌江口三維潮流數(shù)值模型，較好地模擬了甌江口潮流及波浪的時空分布特征；Zhang 等[9]應(yīng)用Quadtree 非結(jié)構(gòu)網(wǎng)格，采用有限體積方法離散方程，利用預(yù)條件的不完全Lu 分解方法加速其收斂性、提高計算效率；Wu 和Li[10]采用自適應(yīng)網(wǎng)格技術(shù)，考慮動水壓強的影響，建立了垂向二維開挖明渠的數(shù)值模型；華祖林等[11]基于四叉樹網(wǎng)格布設(shè)的思想，分層次對不同研究區(qū)域按不同尺度網(wǎng)格對近海水域進(jìn)行嵌套網(wǎng)格聯(lián)合布置，實現(xiàn)了在不同尺度的網(wǎng)格聯(lián)立求解，提高了數(shù)值計算效率，又保證了計算精度。本文采用有限體積法分別基于自適應(yīng)四叉樹計算網(wǎng)格和非結(jié)構(gòu)化三角形計算網(wǎng)格建立了近岸波浪橢圓型緩坡方程的數(shù)值模型，結(jié)合典型物理模型實驗結(jié)構(gòu)對所建模型分別進(jìn)行了分析驗證，并結(jié)合算例分析比較了兩種模型的計算精度和效率。計算結(jié)果表明，兩種數(shù)值模型均可有效地模擬近岸波浪的傳播變形；相對非結(jié)構(gòu)化三角形網(wǎng)格下的模型，基于自適應(yīng)四叉樹網(wǎng)格建立的數(shù)值模型在數(shù)值離散和求解過程中無需引入形函數(shù)、不產(chǎn)生復(fù)雜的交叉項，離散簡單，易于程序?qū)崿F(xiàn)，且節(jié)約計算存儲空間，計算效率高。

2 數(shù)值模型

2.1 橢圓型緩坡方程波浪模型

Berkhoff[1]提出的橢圓型緩坡方程考慮了波浪傳播過程中的折射、反射和繞射效應(yīng)，方程如下：

對于外海入射的簡諧波，如果入射波速度勢已知，則波浪邊界條件可由如下關(guān)系式給出：

（1）外海入射邊界

（2）下游和側(cè)邊界條件

圖 1 四叉樹網(wǎng)格不同分級單元間的通量計算Fig. 1 Flux across different graded cell interfaces

2.2 基于四叉樹網(wǎng)格的數(shù)值模型

2.2.1 四叉樹計算網(wǎng)格

四叉樹網(wǎng)格的生成方法可概括為以下3 個步驟[12-13]：

（1）初始化網(wǎng)格，剖分計算域。初步創(chuàng)建一組很粗但滿足要求的四叉樹網(wǎng)格。

（2）網(wǎng)格自適應(yīng)加密。根據(jù)網(wǎng)格單元屬性，遍歷整個四叉樹，依據(jù)計算精度不斷地遞歸細(xì)分，生成新的四叉樹網(wǎng)格，直到所有四叉樹網(wǎng)格的尺寸滿足必要的精度要求。本文中，以地形水深及波長為剖分標(biāo)準(zhǔn)。

（3）輸出結(jié)果。遍歷整個四叉樹并輸出網(wǎng)格屬性、單元編號、計算節(jié)點坐標(biāo)等。

自適應(yīng)四叉樹網(wǎng)格的生成可能使得相鄰單元間的級別不一（為確保過渡處變量的連續(xù)性，設(shè)定相鄰網(wǎng)格的級差最多為1），本文采用線性插值來估算不同分級單元間通量。四叉樹網(wǎng)格下相鄰計算節(jié)點的布置可歸納為圖1 所示的4 種情形。其中為虛擬單元為當(dāng)前計算單元，不同級別間單元的通量計算可以通過得到，虛擬單元上守恒物理量的q值按圖1 對應(yīng)的4 種情形計算[7-8]。

2.2.2 四叉樹網(wǎng)格下的數(shù)值模型離散

（1）同級別四叉樹網(wǎng)格下的方程離散

在自適應(yīng)四叉樹網(wǎng)格下結(jié)合有限體積法對式（1）在控制容積上積分得：

利用格林函數(shù)法，將式（4）中的第一項的面積分轉(zhuǎn)化為線積分并寫成求和的形式：

式中，A為控制體的面積為控制體到相鄰單元體邊的法向矢量；為控制體界面總個數(shù)，分別為控制體的西、南、東、北4 個方向上界面數(shù)。

I

對于式（4）中的第二項，可用控制單元 中心處的變量來表示對其控制容積積分

由式（5）、式（6）最后可得控制方程（1）的離散形式為

（2）不同級別四叉樹網(wǎng)格下的方程離散

以圖1a 為例推導(dǎo)式（4）的離散方程

整理得：

2.3 基于非結(jié)構(gòu)化三角形網(wǎng)格的數(shù)值模型

2.3.1 非結(jié)構(gòu)化三角形計算網(wǎng)格

本文采用Delaunay 三角形法[14-15]是生成非結(jié)構(gòu)化網(wǎng)格。網(wǎng)格生成可概括為以下5 個步驟：

（1）初始化網(wǎng)格。給定求解計算域邊界上的節(jié)點信息，將其作為初始的三角形頂點，形成一組很粗的但是滿足Delaunay 三角形要求的網(wǎng)格。

（2）不斷向區(qū)域內(nèi)部加點并與已有的點組成新的三角形。引入了兩個幾何參數(shù)—長度標(biāo)尺和外接圓無量綱半徑。循環(huán)以三角形的外心為新的節(jié)點，并利用上面的Delaunay 三角化方法生成新的三角形，直到所有三角形達(dá)到形狀和尺寸的要求。

（3）網(wǎng)格加密。根據(jù)計算的要求，可通過人為加密邊界上網(wǎng)格點的方法在局部區(qū)域加密網(wǎng)格。

（4）網(wǎng)格光順化處理。這里采用Laplace 光順化的方法，即移動每個內(nèi)部網(wǎng)格點使它成為包圍它的多邊形的重心。

（5）輸出結(jié)果。一般包括邊界網(wǎng)格點的個數(shù)及編號；三角形的個數(shù)和編號；三角形三個頂點的編號；所有網(wǎng)格點的坐標(biāo)。

2.3.2 Delaunay 三角形網(wǎng)格下的方程離散

結(jié)合有限體積法，計算節(jié)點布置在單元中心，控制容積即為網(wǎng)格單元。對式（1）在的第一項，利用格林函數(shù)法將面積分轉(zhuǎn)化為線積分并寫成求和形式，

對于式（10）中的第二項，對其控制體積分可用中心處的計算節(jié)點的變量來表示，

由式（10）、式（11）最后可得控制方程（1）的離散方程為

將式（13）代入式（12），就得到以 為計算中心節(jié)點的控制體的離散方程，整理后得到一般形式為

圖 2 格林函數(shù)法計算示意圖Fig. 2 Diagram of Green function method

其中，

2.4 方程迭代的收斂準(zhǔn)則

方程組的離散求解采用GPBiSTAB 法求解[2]。計算時，采用前后兩次迭代相對偏差小于某一固定值的方法作為停止方程組內(nèi)迭代求解的準(zhǔn)則[2]，本文中取

3 數(shù)值模型驗證及結(jié)果分析

3.1 圓形淺灘上的波浪傳播變形

Ito 和Tanimoto[16]于1972 年進(jìn)行了圓形淺灘地形上波浪傳播變形的模型實驗，并給出了3 個斷面的波高實測值。實驗中入射波高H*=0.010 4 m，周期T=0.511 s，波浪沿x軸正向入射。計算域平面示意圖及測量斷面的布置見圖3。在水深為0.15 m 的平底上布置一圓形淺灘，淺灘中心水深為0.05 m，計算域外水深為恒定值，淺灘上水深[16]為

圖 3 圓形淺灘計算域地形示意圖Fig. 3 Topography of a circular shoal zone

圖 4 圓形淺灘計算域網(wǎng)格圖Fig. 4 Numerical grids of a circular shoal zone

計算域取為3.2 m×2.4 m，相鄰的計算節(jié)點間距控制在0.03 m 左右，從而保證一個波長內(nèi)至少有8 個計算節(jié)點。起始邊界采用入射邊界條件，其他邊界取為輻射邊界條件。圖4 給出了兩種網(wǎng)格模型下整個計算區(qū)域網(wǎng)格圖。圖5 給出了兩種網(wǎng)格模型下整個計算區(qū)域數(shù)值模擬相對波高等值線圖，由圖可以看出，由于圓形淺灘的存在使得波浪發(fā)生折射和繞射，在圓形淺灘后明顯存在波高的匯聚區(qū)。圖6 給出了兩種模型在該地形下3 個斷面的數(shù)值解和實測值[16]的比較，可以看出各斷面計算值和實測值均符合的比較良好。表1 是相同計算機下兩種模型的計算效率比較，可以看出四叉樹計算網(wǎng)格下的網(wǎng)格數(shù)量明顯要小于三角形網(wǎng)格，且模型的計算時間優(yōu)于三角形計算網(wǎng)格。

圖 5 圓形淺灘數(shù)值模擬相對波高等值線圖Fig. 5 Contours of simulated wave height for wave propagation on a circular shoal

圖 6 圓形淺灘地形上波高數(shù)值解與實測值[16]的比較Fig. 6 Comparison between numerical simulated and measured[16] wave heights on a circular shoal

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

Berkhoff 等[17]于1982 年進(jìn)行了橢圓形淺灘的模型實驗，實驗中入射波高H*=0.023 2 m，周期T=1.0 s，波浪沿x軸正向入射。實驗給出了相對波高等值線圖以及8 個斷面的實測資料。

計算域平面示意圖及測量斷面的布置見圖7。在均勻的斜坡上布置一橢圓形淺灘，斜坡梯度為1∶50，斜坡梯度方向與波浪入射邊界的法向夾角為20°，水

表 1 圓形淺灘地形上兩種計算網(wǎng)格下的模型效率比較Table 1 Comparison of the efficiency for the model in two kind of numerical grids on a circular shoal

圖 7 橢圓形淺灘計算域地形Fig. 7 Topography of a elliptic shoal zone

深由旋轉(zhuǎn)后的新坐標(biāo)確定[17]

橢圓形邊界定義為[17]

平底區(qū)及斜坡上的水深為[17]

而橢圓形淺灘上的水深為[17]

圖 8 橢圓形淺灘計算域網(wǎng)格圖Fig. 8 Numerical grids of an elliptic shoal zone

計算域取21.5 m×20.0 m，相鄰的計算節(jié)點間距控制在0.12 m 左右，從而保證一個波長內(nèi)至少有8 個計算節(jié)點。起始邊界采用入射邊界條件，其他邊界取為輻射邊界條件。圖8 是兩種模型的計算網(wǎng)格圖，圖9和圖10 給出了整個計算區(qū)域及8 個斷面相對波高等值線圖。由圖9 可以看出由于淺灘的存在使得波浪發(fā)生折射和繞射，在橢圓形淺灘后明顯存在波高的匯聚區(qū)。由圖10 可以看出兩種模型計算結(jié)果和實測值[17]都符合較好，而四叉樹網(wǎng)格下橢圓型緩坡方程的數(shù)值模型計算結(jié)果和實驗結(jié)果總體吻合更好。表2 是相同計算機下的兩種模型效率比較，可以看出四叉樹計算網(wǎng)格下的網(wǎng)格數(shù)量明顯要小于三角形網(wǎng)格，且模型的計算時間優(yōu)于三角形計算網(wǎng)格。

圖 9 橢圓形淺灘數(shù)值模擬相對波高等值線圖Fig. 9 Contours of simulated wave height for wave propagation on an elliptic shoal

圖 10

4 結(jié)論

本文分析比較了兩種適于不規(guī)則水域波浪模擬的橢圓型緩坡方程數(shù)值模型。兩種模型均采用有限體積法離散，分別基于四叉樹網(wǎng)格和非結(jié)構(gòu)化三角形網(wǎng)格建立。首先結(jié)合近岸波浪傳播的典型物理模型實驗對兩種數(shù)值模型分別進(jìn)行了驗證，并結(jié)合計算結(jié)果對比分析了兩種模型的計算精度和效率。計算結(jié)果表明，兩種數(shù)值模型均可有效地模擬近岸波浪的傳播變形；相對非結(jié)構(gòu)化三角形網(wǎng)格下的模型，基于自適應(yīng)四叉樹網(wǎng)格建立的數(shù)值模型在數(shù)值離散和求解過程中無需引入形函數(shù)、不產(chǎn)生復(fù)雜的交叉項，離散簡單，特別是特別三級以內(nèi)的自適應(yīng)四叉樹網(wǎng)格易于程序?qū)崿F(xiàn)，且節(jié)約計算存儲空間，計算效率高。

圖 10 橢圓形淺灘地形上波高數(shù)值解與實測值[17]的比較Fig. 10 Comparison between numerical simulated and measured [17] wave heights on an elliptic shoal

表 2 橢圓形淺灘地形上兩種計算網(wǎng)格下的模型效率比較Table 2 Comparison of the efficiency for the model in two kind of numerical grids on an elliptic shoal