適于模擬不規(guī)則水域波浪的緩坡方程兩種數(shù)值模型比較

李巧生，唐軍，呂義港

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

1 引言

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

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

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

2 數(shù)值模型

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

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

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

（1）外海入射邊界

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

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

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

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

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

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

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

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

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

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

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

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

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

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

I

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

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

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

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

整理得：

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

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

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

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

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

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

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

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

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

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

對(duì)于式（10）中的第二項(xiàng)，對(duì)其控制體積分可用中心處的計(jì)算節(jié)點(diǎn)的變量來表示，

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

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

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

其中，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

橢圓形邊界定義為[17]

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

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

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

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

圖 9 橢圓形淺灘數(shù)值模擬相對(duì)波高等值線圖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í)驗(yàn)對(duì)兩種數(shù)值模型分別進(jìn)行了驗(yàn)證，并結(jié)合計(jì)算結(jié)果對(duì)比分析了兩種模型的計(jì)算精度和效率。計(jì)算結(jié)果表明，兩種數(shù)值模型均可有效地模擬近岸波浪的傳播變形；相對(duì)非結(jié)構(gòu)化三角形網(wǎng)格下的模型，基于自適應(yīng)四叉樹網(wǎng)格建立的數(shù)值模型在數(shù)值離散和求解過程中無需引入形函數(shù)、不產(chǎn)生復(fù)雜的交叉項(xiàng)，離散簡(jiǎn)單，特別是特別三級(jí)以內(nèi)的自適應(yīng)四叉樹網(wǎng)格易于程序?qū)崿F(xiàn)，且節(jié)約計(jì)算存儲(chǔ)空間，計(jì)算效率高。

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

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