陳善群，吳 昊，湯潤超

（安徽工程大學(xué) 建筑工程學(xué)院，安徽 蕪湖241000）

0 引 言

眾所周知，波浪由深水區(qū)向沿岸淺水區(qū)傳播過程中會出現(xiàn)波浪折射、波浪衍射、波浪淺化及波浪破碎等復(fù)雜現(xiàn)象。如何準(zhǔn)確模擬這些復(fù)雜現(xiàn)象一直是港口、海岸及近海工程等相關(guān)研究領(lǐng)域的熱點問題。

目前比較成熟的波浪模擬數(shù)值模型有兩大類別。一類是基于Boussinesq方程的波浪數(shù)值模型，該類模型具有非線性和色散特性，能準(zhǔn)確有效地模擬波浪的演化過程，尤其在淺水區(qū)的波浪演化過程的模擬上已體現(xiàn)出較大優(yōu)勢[1-2]。目前該類模型的研究拓展主要集中于色散精度的提高[3-4]、湍流模型的加入[5]、自由面溶質(zhì)的影響及輸運[6]等幾個方面。這些研究拓展的實現(xiàn)需在原數(shù)值模型內(nèi)加入大量的附加控制方程，使得控制方程組的求解難度明顯增加，一定程度上阻礙了該類波浪數(shù)值模型的持續(xù)性發(fā)展。另一類是直接求解Navier-Stokes方程組結(jié)合自由面追蹤技術(shù)來模擬波浪的演化過程。該類模型具有適用性強，精度較高等優(yōu)點，已在工程實際中得到大量使用[7-9]。但研究者們發(fā)現(xiàn)該類模型同樣存在明顯缺陷，主要體現(xiàn)在：（1）消耗大量計算機時，難以拓展到大尺度計算域；（2）自由面壓力邊界的準(zhǔn)確求解難度較大，直接影響自由面速度場的準(zhǔn)確性；（3）計算網(wǎng)格一般采用笛卡爾坐標(biāo)系，難以準(zhǔn)確貼合自由面和不規(guī)則區(qū)域，直接影響求解準(zhǔn)確性等幾個方面。

綜合上述兩類波浪數(shù)值模型的優(yōu)點與不足，研究者們另行發(fā)展了一類非靜力學(xué)波浪數(shù)值模型[10-15]。該類模型的主要特點有：（1）建立自由面豎直方向水位與橫向坐標(biāo)的函數(shù)關(guān)系，減少自由面追蹤消耗的大量計算機時；（2）控制方程采用非靜力學(xué)方程，方程中的壓力被分解成靜力學(xué)壓力和非靜力學(xué)壓力兩部分，可結(jié)合有限體積或有限差分方法對控制方程進(jìn)行離散求解。研究表明，該類模型非常適用于自由面壓力邊界的準(zhǔn)確求解，并在淺水區(qū)波浪演化過程、自由面湍流運動以及溶質(zhì)輸運等方面的數(shù)值計算中都取得了較好的效果。

另外，針對早期非靜力學(xué)數(shù)值模型中存在的豎直方向網(wǎng)格數(shù)較少使得波浪的色散特性難以被準(zhǔn)確評估的難點問題，研究者們普遍認(rèn)為豎直方向頂層計算單元壓力的準(zhǔn)確求解是解決問題的關(guān)鍵。為此，Stelling和Zijlema[16]提出了一種Keller-box網(wǎng)格劃分方法替代常用的交錯網(wǎng)格，使得壓力處于計算單元界面位置而不再處于計算單元中心位置。這樣做無需近似，自由面的壓力可精確設(shè)置為零。Yuan和Wu[17]提出了一種交錯網(wǎng)格框架下的積分方法，直接避免了豎直方向頂層計算單元的靜力學(xué)壓力設(shè)定問題。Young和Wu[18]對類Boussinesq型方程進(jìn)行解析近似求解，結(jié)合參考速度，得到了豎直方向頂層計算單元非靜力學(xué)壓力分布。以上研究成果已在一定程度上妥善解決了豎直方向網(wǎng)格數(shù)較少使得波浪的色散特性難以被準(zhǔn)確評估的難點問題。然而，早期非靜力學(xué)數(shù)值模型在模擬碎波帶波浪演化及破碎波浪沿斜坡爬高等波浪復(fù)雜演化過程時遇到的難以準(zhǔn)確模擬的問題仍沒有得到妥善解決。Zijlema和Stelling[19]則認(rèn)為需在非靜力學(xué)模型中結(jié)合適當(dāng)?shù)募げú蹲剿惴▉韺崿F(xiàn)對波浪復(fù)雜演化過程的準(zhǔn)確模擬。與此同時，Godunov型格式[20-23]作為一類成熟的激波捕捉算法已在含有間斷流動的計算流體力學(xué)問題中得到廣泛使用。大量算例證明，Godunov型格式無需添加任何判據(jù)，就具有直接判定并準(zhǔn)確捕捉波浪破碎位置的能力，被認(rèn)為是準(zhǔn)確模擬波浪復(fù)雜演化過程的理想數(shù)值算法。

本文擬在早期非靜力學(xué)波浪模擬數(shù)值模型的基礎(chǔ)之上，結(jié)合研究者們對數(shù)值模型的改進(jìn)方法，并在數(shù)值模型中附加Godunov型格式以實現(xiàn)激波捕捉，發(fā)展一種能準(zhǔn)確模擬波浪復(fù)雜演化過程的非靜力學(xué)數(shù)值模型。該數(shù)值模型控制方程采用笛卡爾坐標(biāo)系下的不可壓縮Navier-Stokes方程組。為準(zhǔn)確模擬計算域不規(guī)則底部與自由面之間水體的運動，使用σ坐標(biāo)系對笛卡爾坐標(biāo)系進(jìn)行坐標(biāo)變換，后續(xù)采用有限體積和有限差分混合方法結(jié)合Godunov型格式對控制方程進(jìn)行離散求解。為便于結(jié)合Godunov型格式，數(shù)值模型中速度變量定義在計算單元的中心位置，同時將非靜力學(xué)壓力定義在計算單元豎直方向界面位置。為進(jìn)一步提高數(shù)值模型的時空分辨率，計算單元界面的通量計算采用HLL黎曼求解器，時間步迭代采用二階非線性強穩(wěn)定性保持龍格庫塔（SSP Runge-Kutta）格式。后續(xù)擬將非靜力學(xué)數(shù)值模型用于經(jīng)典波浪復(fù)雜演化問題的數(shù)值模擬，并與解析解或?qū)嶒灲Y(jié)果進(jìn)行對比，驗證該數(shù)值模型的準(zhǔn)確性與適用性。

1 控制方程

本文選取三維笛卡爾坐標(biāo)系下的不可壓縮Navier-Stokes方程組作為非靜力學(xué)數(shù)值模型的控制方程，表達(dá)式如下：

為準(zhǔn)確貼合計算域不規(guī)則底部及自由面，本文采用Lin和Li[24]提出的σ坐標(biāo)系對笛卡爾坐標(biāo)系進(jìn)行坐標(biāo)變換，變換的表達(dá)式為：

將變換式（3）代入笛卡爾坐標(biāo)系下的控制方程（（1）式）后推導(dǎo)得到σ坐標(biāo)系（x，y，σ）下的控制方程，表達(dá)式為：

連續(xù)性方程（（5）式）通過積分結(jié)合豎直方向速度ω需滿足的邊界條件，可得σ坐標(biāo)系下自由面運動控制方程，表達(dá)式如下：

2 控制方程的離散求解

本文采用有限體積和有限差分混合方法結(jié)合Godunov型格式對控制方程（（6）式）進(jìn)行離散求解。為便于結(jié)合Godunov型格式，選用一種另類的交錯網(wǎng)格框架來劃分計算域，其中速度變量定義在計算單元的中心位置，非靜力學(xué)壓力定義在計算單元豎直方向界面位置，具體如圖2所示。以下章節(jié)為具體闡述控制方程的離散求解過程。

2.1 時間步迭代

為獲得二階時間精度，本文的時間步迭代選取Gottlieb等[25]提出的兩段式SSP Runge-Kutta格式，具體過程如下：

第一階段，采用一種常見的一階精度兩步投影法對（6）式進(jìn)行時間離散，獲取中間量U（）1：

式中：Un代表時間步n時U的取值；U*為兩步投影法計算格式中的中間量；U（）1為經(jīng)過第一階段迭代計算后得到的中間量。

第二階段，利用上階段采用的投影法對（8）式再次進(jìn)行迭代計算（兩段式的兩步投影法構(gòu)成二階非線性強穩(wěn)定性保持龍格庫塔（SSP Runge-Kutta）格式），獲取時間步n+1時U的取值Un+1：

（9）式與（11）式中的Sp為非靜力學(xué)壓力源項，非靜力學(xué)壓力p需通過求解Poisson方程得到，Poisson方程的推導(dǎo)與求解將在后續(xù)章節(jié)中闡述。

（8）式與（10）式中的 Δt為自適應(yīng)時間步長，滿足 Courant-Friedrichs-Lewy （CFL）條件，Δt的表達(dá)式為：

式中：C為Courant常數(shù)，為保證數(shù)值模型計算的精度和穩(wěn)定性，C取值0.5。

2.2 空間離散

本文選取Liang和Marche[26]提出的二階Godunov型有限體積方法對（8）式與（10）式進(jìn)行空間離散求解，具體過程如下：

首先對通量F和G以及源項Sh進(jìn)行修正，將D=h+η代入源項Sh的表達(dá)式，對Sh進(jìn)行簡單分解：

將通量F和G分別與（14）式和（15）式等式右邊的第一項相減后，可得修正后的通量F和G以及源項Sh的表達(dá)式：

接下來對計算單元i中的守恒變量U進(jìn)行分段線性空間重構(gòu)，以x方向為例：

守恒變量U在計算單元界面i+1/2左側(cè)和右側(cè)的估值分別為：

從（9）式與（11）式中可以看出，非靜力學(xué)壓力p是求解非靜力學(xué)速度場（u，v， ）w 的前提，兩者之間存在明顯的依賴關(guān)系，表達(dá)式為：

式中：m=1，2代表兩段式SSP Runge-Kutta格式中的階段數(shù)。

求解非靜力學(xué)壓力p所需的Poisson方程由連續(xù)性方程（（5）式）演化得來，先將（3）式和（4）式代入（5）式，得出連續(xù)性方程（（5）式）的坐標(biāo)變換形式：

采用二階中心有限差分方法對（22）式進(jìn)行空間離散，得到對應(yīng)的線性方程：

2.3 邊界條件

（1）自由面邊界條件

不考慮附加風(fēng)效應(yīng)，自由面切向應(yīng)力應(yīng)等于零，同時自由面豎直方向速度w滿足運動學(xué)邊界條件，以及自由面的非靜力學(xué)壓力p在求解Poisson方程時應(yīng)設(shè)置為零，即：

（2）計算域底部邊界條件

計算域底部法向速度w滿足運動學(xué)邊界條件，同時水平方向速度滿足滑移邊界條件，以及計算域底部剪切應(yīng)力滿足粗糙固壁邊界條件，即：

計算域底部非靜力學(xué)壓力p滿足Neumann邊界條件，即：

（3）其余邊界條件

豎直方向邊界是固壁時，滿足滑移邊界條件，即法向速度和切向應(yīng)力設(shè)置為零，法向壓力梯度設(shè)置為零；豎直方向邊界為波浪入口邊界時，自由面高程和速度大小通過解析解來施加；側(cè)向邊界滿足滑移邊界條件。

3 算例驗證

3.1 淹沒式堤壩上波浪傳播

為驗證非靜力學(xué)數(shù)值模型對波浪淺化等復(fù)雜演化過程模擬的精確性與適用性，選取Beji和Battjes[28]關(guān)于淹沒式堤壩上周期性波浪傳播過程的水槽實驗作為研究對象。根據(jù)Beji和Battjes的實驗裝置，建立相應(yīng)的計算域，如圖3所示。其中，計算域分為波浪水槽和消波區(qū)兩部分。波浪水槽長25.0 m，高0.6 m，水槽中靜水高度初始為0.4 m，淹沒式堤壩外形為梯形，高度為0.3 m，堤壩向波坡面坡度1：20，背波坡面坡度1：10。波浪水槽中初始水位距堤壩頂部0.1 m。消波區(qū)位于波浪水槽右部，長10.0 m，采用Larsen和Dancy[29]提出的海綿層消波技術(shù)。數(shù)值水槽左邊壁為周期性波浪入口邊界，選取周期T=2.02 s，振幅A=1.0 cm和周期T=1.25 s，振幅A=1.2 cm兩種規(guī)則波形條件。數(shù)值水槽內(nèi)共設(shè)置a-f六個水位監(jiān)測點，距離水槽左邊壁的距離xa-xf分別為10.5 m、12.5 m、13.5 m、14.5 m、15.7 m和17.3 m。計算域水平方向網(wǎng)格數(shù)為1750，豎直方向網(wǎng)格數(shù)為3，計算單元總數(shù)為5250。

圖4和圖5分別給出了周期T=2.02 s，振幅A=1.0 cm和周期T=1.25 s，振幅A=1.2 cm兩種波形條件下a-f六個水位監(jiān)測點自由面高程數(shù)值計算與實驗結(jié)果的對比。從圖中可以看出，兩種波形條件下，非靜力學(xué)數(shù)值模型對六個監(jiān)測點的自由面高程預(yù)測結(jié)果與Beji和Battjes的實驗結(jié)果整體吻合較好。另外，從圖4（c）-（f）中可以看出，對于淹沒式堤壩上波浪傳播過程中出現(xiàn)非線性淺水波這一復(fù)雜演化現(xiàn)象，非靜力學(xué)數(shù)值模型同樣能較準(zhǔn)確地模擬。

3.2 孤立波沿斜坡爬高

3.3 海底滑坡誘發(fā)海嘯

滑坡體輪廓曲線為截斷雙曲正割函數(shù)，滿足關(guān)系式：

式中：ut=1.17 m/s，為滑坡體的自由沉降速度；a0=1.12 m/s2，為滑坡體初始加速度。

圖9給出了海底滑坡誘發(fā)海嘯過程中四個水位監(jiān)測點 （x0，y0）自由面高程數(shù)值計算與實驗結(jié)果的對比。圖中tC=0.900+7.07d，為文獻(xiàn)[31]給出的經(jīng)驗參考時間。從圖中可以看出，非靜力學(xué)數(shù)值模型對四個監(jiān)測點的自由面高程預(yù)測結(jié)果與Enet和Grilli的實驗結(jié)果整體吻合較好。對于海底滑坡誘發(fā)海嘯過程，非靜力學(xué)數(shù)值模型能較準(zhǔn)確地數(shù)值模擬。

據(jù)由Enet和Grilli的實驗數(shù)據(jù)可擬合得出滑坡體的運動方程，關(guān)系式為：

4 結(jié) 論

針對早期非靜力學(xué)波浪模擬數(shù)值模型的不足之處，結(jié)合研究者們對數(shù)值模型的改進(jìn)方法，并在數(shù)值模型中附加激波捕捉Godunov型格式，本文發(fā)展了一種能準(zhǔn)確模擬波浪復(fù)雜演化過程的非靜力學(xué)數(shù)值模型。該數(shù)值模型選取σ坐標(biāo)系下的不可壓縮Navier-Stokes方程組作為控制方程，利用σ坐標(biāo)系貼合不規(guī)則底部與自由面之間的計算域；為便于結(jié)合Godunov型格式，將速度變量定義在計算單元的中心位置，同時將非靜力學(xué)壓力定義在計算單元豎直方向界面位置；采用有限體積和有限差分混合方法結(jié)合Godunov型格式對控制方程進(jìn)行空間離散；為進(jìn)一步提高時空分辨率，利用HLL黎曼求解器求解計算單元之間的通量以實現(xiàn)激波捕捉，并采用二階非線性強穩(wěn)定性保持龍格庫塔（SSP Runge-Kutta格式）進(jìn)行時間步迭代。將非靜力學(xué)數(shù)值模型用于數(shù)值模擬淹沒式堤壩上波浪傳播、孤立波沿斜坡爬高和海底滑坡誘發(fā)海嘯三個波浪復(fù)雜演化過程，數(shù)值計算結(jié)果與相應(yīng)的實驗結(jié)果較為吻合，說明該非靜力學(xué)模型對波浪淺化、波浪破碎、滑坡誘發(fā)海嘯等復(fù)雜現(xiàn)象的數(shù)值模擬具有較高的精確性和適用性，可望為波浪復(fù)雜演化過程的數(shù)值模擬提供一種較為準(zhǔn)確有效的研究手段。