摩擦項の半陰解法

前回の続きです。

1次元なのに発散することがあり、試行錯誤した結果、下記で安定しました。1次元だと思って舐めていました。

空間離散化

• 連続式: Upwind法（流れ方向に応じた風上差分）
• 運動量式: 中心差分

時間積分
• 連続式: 陽解法
• 運動量式: 半陰解法（Semi-implicit法）

効いたのは摩擦項への半陰解放の適用。摩擦項Sfの非線形性により陽解法では不安定になりやすいため、以下の線形化を実施しました。

Manning式のSfを現在の流速u_iの周りでTaylor展開：
Sf(u) ≈ Sf(u_i) + (u - u_i)·∂Sf/∂u|_{u_i}
     = Sf_const + Sf_coef·(u - u_i)

ここで：
Sf_const = n²·u_i·|u_i| / Rh^(4/3)          （定数項）
Sf_coef  = 2·n²·|u_i| / Rh^(4/3)            （線形係数）

運動量式の離散化
元の式：
∂u/∂t + u·∂u/∂x = -g·∂h/∂x + g(S₀ - Sf)
離散化：
(u_new - u_i)/dt + u_i·du_dx = -g·dh_dx + g·S₀ - g·Sf

線形化した摩擦項を代入して整理すると：
u_new·(1 + dt·g·Sf_coef) = u_i - dt·(u_i·du_dx + g·dh_dx - g·S₀ + g·Sf_const - g·sf_coef·u_i)
u_new = RHS / (1 + dt·g·Sf_coef)
RHS = u_i - dt·(u_i·du_dx + g·dh_dx - g·S₀ + g·Sf_const - g·sf_coef·u_i)

あとでプロと話していたら、1次元で不安定になりやすいのはやはり摩擦項で、半陰解法を用いて安定化させるのは実務でよくやる対応らしいです。1次元だから皆さん同じ方向に向かうのかもしれません。