一種能準確模擬波浪復雜演化過程的非靜力學數值模型

陳善群，吳 昊，湯潤超

（安徽工程大學 建筑工程學院，安徽 蕪湖241000）

0 引 言

眾所周知，波浪由深水區向沿岸淺水區傳播過程中會出現波浪折射、波浪衍射、波浪淺化及波浪破碎等復雜現象。如何準確模擬這些復雜現象一直是港口、海岸及近海工程等相關研究領域的熱點問題。

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

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

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

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

1 控制方程

本文選取三維笛卡爾坐標系下的不可壓縮Navier-Stokes方程組作為非靜力學數值模型的控制方程，表達式如下：

為準確貼合計算域不規則底部及自由面，本文采用Lin和Li[24]提出的σ坐標系對笛卡爾坐標系進行坐標變換，變換的表達式為：

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

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

2 控制方程的離散求解

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

2.1 時間步迭代

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

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

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

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

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

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

式中：C為Courant常數，為保證數值模型計算的精度和穩定性，C取值0.5。

2.2 空間離散

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

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

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

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

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

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

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

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

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

2.3 邊界條件

（1）自由面邊界條件

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

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

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

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

（3）其余邊界條件

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

3 算例驗證

3.1 淹沒式堤壩上波浪傳播

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

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

3.2 孤立波沿斜坡爬高

3.3 海底滑坡誘發海嘯

滑坡體輪廓曲線為截斷雙曲正割函數，滿足關系式：

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

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

據由Enet和Grilli的實驗數據可擬合得出滑坡體的運動方程，關系式為：

4 結 論

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