引言

成功實現一個穩定且精確的水動力學模型，關鍵在于妥善處理源項和邊界條件。這兩個環節是數值格式產生非物理振蕩和誤差的主要來源。本章將詳細介紹“守恒-平衡”（well-balanced）格式的核心技術，以及通過“虛擬單元”實現各類物理邊界條件的方法，并最終給出一個完整的算法流程。

1. 守恒-平衡格式：精確維持“靜湖”狀態

• 1.1 源項離散化的挑戰
  一個核心挑戰是數值格式能否精確維持“靜湖”（lake at rest）狀態。在這種狀態下，流速為零，水面是平的（$h+z=const$）。此時，動量方程中的壓力梯度項與河床坡度源項應精確抵消。然而，簡單的離散方式會打破這種平衡，導致靜止的水體產生虛假的流動。一個能夠精確保持“靜湖”狀態的格式被稱為滿足C特性（C-property）或守恒-平衡的格式。

• 1.2 靜水重構法 (Hydrostatic Reconstruction)
  靜水重構法是實現守恒-平衡格式的先進技術。其核心思想是修改輸入到黎曼求解器的左右狀態，使得在“靜湖”條件下，通量梯度與經過特殊離散的源項能夠精確平衡。

  • 實現步驟:
    1. 定義界面高程: 在界面 $x_{i+1/2}$ 處，定義界面河床高程為兩側單元河床高程的最大值：$z_{i+1/2}=\max (z_i,z_{i+1})$。
    2. 重構界面水深: 在界面左右兩側，分別重構水深值，確保水深非負：
       $h_{L,i+1/2}=\max (0,h_i+z_i?z_{i+1/2})$
       $h_{R,i+1/2}=\max (0,h_{i+1}+z_{i+1}?z_{i+1/2})$
    3. 構造重構狀態: 使用重構后的水深和原有的單元平均流速，構造用于黎曼求解器的新左右狀態 ${\mathbf{U}}_{L,i+1/2}^{?}$ 和 ${\mathbf{U}}_{R,i+1/2}^{?}$。
    4. 計算中間通量: 將重構后的狀態代入近似黎曼求解器（如HLLC），得到一個中間通量 ${\mathbf{F}}_{i+1/2}^{?}$。
    5. 離散源項: 采用一種與重構過程相匹配的特殊方式離散源項，以確保在“靜湖”條件下，通量梯度與源項精確平衡。

• 1.3 摩阻坡降的離散化
  摩阻項在流速為零時為零，不參與靜水重構。在動態模擬中，它通常在單元中心進行顯式評估。

2. 通過虛擬單元處理邊界條件

在FVM中，邊界條件通常通過在計算域兩端設置“虛擬單元”（Ghost Cells）來實現。這些單元內的物理量根據邊界條件人為設定，為邊界處的黎曼求解器提供“外部”狀態。靜水重構法同樣適用于邊界處的通量計算。

• 上游邊界 (給定流量 $Q_{in}(t)$):

  • 虛擬單元設置: 動量根據給定流量設定 $q_0=Q_{in}(t)/B_0$；水深通常從計算域內部進行零階外插 $h_0=h_1$；河床高程也外插 $z_0=z_1$。

• 下游邊界 (給定水位 ${\eta }_{out}(t)$):

  • 虛擬單元設置: 水深由給定水位計算 $h_{N+1}={\eta }_{out}(t)?z_{N+1}$；動量從內部外插 $q_{N+1}=q_N$。

• 固壁邊界 (反射邊界):

  • 物理條件: 法向流速為零 $u=0$。
  • 虛擬單元設置: 動量反號反射 $q_0=?q_1$；水深對稱反射 $h_0=h_1$。

3. 時間積分與穩定性

• 時間步進: 簡單的顯式歐拉前差格式即可實現時間積分，為了提高精度，也可以采用更高階的Runge-Kutta方法。
• CFL穩定性條件: 顯式格式的穩定性受到CFL條件的嚴格限制。其物理本質是，在一個時間步 $\Delta t$ 內，數值信息的傳播距離不能超過一個空間步長 $\Delta x$。
  $\Delta t\le C_r\frac{\Delta x}{{\max }_i(|u_i|+c_i)}$
  其中 $C_r$ 是Courant數，通常取0.5到0.9之間。
• 自適應時間步長: 在每個時間步，根據當前流場計算出全局最大波速，然后用CFL條件動態計算下一個時間步長。這是一種高效且安全的策略。

以下二階龍格-庫塔法（Heun法）的代碼實現：

def solve_one_step_rk2(U_n, dx, dt, boundary_conditions, geo_params):"""使用二階龍ge-Kutta法（Heun法）執行一個完整的時間步求解。Args:U_n (np.ndarray): n 時刻的守恒量數組dx (float): 空間步長dt (float): 時間步長boundary_conditions: 邊界條件信息geo_params: 河道幾何參數Returns:np.ndarray: n+1 時刻的守恒量數組"""# 這是一個高層函數，它調用了之前章節定義的各種子函數# calculate_rhs 封裝了設置虛擬單元、靜水重構、黎曼求解、計算源項和通量梯度的所有空間離散步驟# --- 階段一: 預測 ---# 根據 U_n 計算時間導數 K1K1 = calculate_rhs(U_n, dx, boundary_conditions, geo_params)# 計算中間狀態 U*U_star = U_n + dt * K1# --- 階段二: 校正 ---# 根據中間狀態 U* 計算時間導數 K2K2 = calculate_rhs(U_star, dx, boundary_conditions, geo_params)# --- 更新最終解 ---# 使用兩個導數的平均值來更新U_np1 = U_n + 0.5 * dt * (K1 + K2)return U_np1

4. 算法總成與實踐建議

綜合以上所有技術，一個完整的時間步進算法流程（偽代碼）如下：

1. 計算時間步長: 遍歷所有單元，找到全局最大特征波速 ${\lambda }_{max}$，并根據CFL條件計算 $\Delta t$。
2. 設置虛擬單元: 根據各類邊界條件，設定計算域兩端虛擬單元的狀態 $(h,q,z)$。
3. 計算界面通量: 遍歷所有內部和邊界界面：
   • 執行靜水重構，得到用于黎曼求解器的左右狀態 $U_L^{?},U_R^{?}$。
   • 調用HLLC黎曼求解器，計算中間通量 $F^{?}$。
4. 更新單元守恒量: 遍歷所有物理單元：
   • 計算守恒-平衡源項（重力坡度項）和摩阻源項。
   • 使用前向歐拉法或龍格-庫塔法，根據通量梯度和源項更新單元的守恒狀態 $U_i^{n+1}$。
5. 更新時間：$t=t+\Delta t$。

在實踐中，強烈推薦使用 HLLC 求解器，因為它在精度、魯棒性和計算成本之間取得了最佳平衡，是現代通用一維水動力模型的首選。

5. 總結與展望

通過四篇文章的旅程，我們系統地構建了一個基于顯式有限體積法，采用HLLC求解器和靜水重構技術的先進一維水動力模型。在投入實際應用前，模型必須經過嚴格的驗證與確認（Verification and Validation），包括通過“靜湖試驗”、“干河床潰壩波”和“駝峰上的跨臨界流”等基準算例來檢驗代碼的正確性。

希望這個系列能為您提供一個堅實的理論基礎和清晰的實踐指南。