境界適合格子による平面二次元河川流計算の式形についての考察（途中段階） - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

前回打合せ時に話した件です。ひたすら数式です。結論出ずです。

はじめに
パターン1：通常型
パターン2：反変成分型
パターン3：Delft型
パターンの比較検討(途中段階)
参考サイトなど

はじめに

河川流の平面二次元又は三次元計算では、境界適合格子（BFC：Boundary Fitting Coordinates）を用いることが一般的です。これは、境界条件の設定の容易さと解像度に対する計算格子数の少なさの点で他手法と比べて優位なためです。特に、後者は、河川流計算は他の流体計算と比べてイテレーションの回数が格段に多くなるため、計算資源の観点からも重要になります。

改めて、BFCのモデルを書いてみようと考えていますが、式形が複数パターン想定されるので、どれを採用するかの比較してみます。なお、一般座標系の河川流計算での軸は、流下方向を $\xi$ 、横断方向を $\eta$ と定義することが一般的です。

元農工研の浪平さん(...)の論文2012, 一般座標系における平面 2 次元流れと河床変動の数値シミュレーション手法に関する比較研究も参考になりますが、もう少しシンプルに整理して、以下の3パターンに分類しました。

通常型
反変成分型
Delft型

それぞれについては以降に詳述しています。

パターン1：通常型

デカルト座標系の連続式、運動方程式を一般座標系へ座標変換しただけの式形です。

デカルト座標系と同じ項数でシンプルな式系になります。

なお、下式は、「2001, 水理公式集例題プログラム集」または「1999, 水工学における計算機利用の講習会　講義集」の故長田さんのテキストから引用しています。

リンク

連続式

$\begin{align} \frac{\partial}{\partial t}\left(\frac{h}{J}\right) +\frac{\partial}{\partial \xi}\left(\frac{U h}{J}\right) +\frac{\partial}{\partial \eta}\left(\frac{V h}{J}\right)=0 \end{align}$

運動方程式

$\begin{align} & \frac{\partial}{\partial t}\left(\frac{M}{J}\right)+\frac{\partial}{\partial \xi}\left(\frac{U M}{J}\right)+\frac{\partial}{\partial \eta}\left(\frac{V M}{J}\right)=-g h\left(\frac{\xi_x}{J} \frac{\partial z_s}{\partial \xi}+\frac{\eta_x}{J} \frac{\partial z_s}{\partial \eta}\right)-\frac{\tau_{b x}}{\rho J} \\ & \quad+\frac{\xi_x}{J} \frac{\partial}{\partial \xi}\left(-\overline{u^{\prime 2}} h\right)+\frac{\xi_y}{J} \frac{\partial}{\partial \xi}\left(-\overline{u^{\prime} v^{\prime}} h\right)+\frac{\eta_x}{J} \frac{\partial}{\partial \eta}\left(-\overline{u^{\prime 2}} h\right)+\frac{\eta_y}{J} \frac{\partial}{\partial \eta}\left(-\overline{u^{\prime} v^{\prime}} h\right) \end{align}$

$\begin{align} & \frac{\partial}{\partial t}\left(\frac{N}{J}\right)+\frac{\partial}{\partial \xi}\left(\frac{U N}{J}\right)+\frac{\partial}{\partial \eta}\left(\frac{V M}{J}\right)=-g h\left(\frac{\xi_y}{J} \frac{\partial z_s}{\partial \xi}+\frac{\eta_y}{J} \frac{\partial z_s}{\partial \eta}\right)-\frac{\tau_{b y}}{\rho J} \\ & \quad+\frac{\xi_x}{J} \frac{\partial}{\partial \xi}\left(-\overline{u^{\prime} v^{\prime}} h\right)+\frac{\xi_y}{J} \frac{\partial}{\partial \xi}\left(-\overline{v^{\prime 2}} h\right)+\frac{\eta_x}{J} \frac{\partial}{\partial \eta}\left(-\overline{u^{\prime} v^{\prime}} h\right)+\frac{\eta_y}{J} \frac{\partial}{\partial \eta}\left(-\overline{v^{\prime 2}} h\right) \end{align}$

導出は作成中

この式形は、以下のWuさんの式形（論文：2004, Depth-Averaged Two-Dimensional Numerical Modeling of Unsteady Flow and Nonuniform Sediment Transport in Open Channelsより引用。以下の書籍にも記載があります。）に近いと思います。

リンク

連続式・運動方程式

$\begin{align} \frac{\partial}{\partial t}(J \rho h \phi) +\frac{\partial}{\partial \xi_m}\left(J \rho h\left(\hat{u}_m \phi-\Gamma_\phi \alpha_j^m \alpha_j^n \frac{\partial \phi}{\partial \xi_n}\right)\right)=J \rho h S \phi \end{align}$

ここに、 $\phi=(1,U,V)$ とする。

パターン2：反変成分型

パターン1の一般座標系の運動方程式を連立させて、反変成分の単位幅流量を運動方程式のプリミティブ変数になるように式変形したものです。

パターン1と比べて、項数が増えて複雑な式系になります。

運動方程式

$\begin{align} & \frac{\partial}{\partial t}\left(\frac{Q^{\xi}}{J}\right)+\frac{\partial}{\partial \xi}\left(\frac{U Q^{\xi}}{J}\right)+\frac{\partial}{\partial \eta}\left(\frac{V Q^{\xi}}{J}\right)-\frac{M}{J}\left(U \frac{\partial \xi_x}{\partial \xi}+V \frac{\partial \xi_x}{\partial \eta}\right)-\frac{N}{J}\left(U \frac{\partial \xi_y}{\partial \xi}+V \frac{\partial \xi_y}{\partial \eta}\right) \\ & =-g h\left(\frac{\xi_x^2+\xi_y^2}{J} \frac{\partial z_s}{\partial \xi}+\frac{\xi_x \eta_x+\xi_y \eta_y}{J} \frac{\partial z_s}{\partial \eta}\right)-\frac{\tau_b^{\xi}}{\rho J}+\frac{\xi_x^2}{J} \frac{\partial}{\partial \xi}\left(-\overline{u^{\prime 2}} h\right)+\frac{\xi_x \eta_x}{J} \frac{\partial}{\partial \eta}\left(-\overline{u^{\prime 2}} h\right) \\ & +\frac{\xi_y^2}{J} \frac{\partial}{\partial \xi}\left(-\overline{v^{\prime 2}} h\right)+\frac{\xi_y \eta_y}{J} \frac{\partial}{\partial \eta}\left(-\overline{v^{\prime 2}} h\right)+\frac{\xi_x \eta_y+\xi_y \eta_x}{J} \frac{\partial}{\partial \eta}\left(-\overline{u^{\prime} v^{\prime}} h\right)+\frac{2 \xi_x \xi_y}{J} \frac{\partial}{\partial \xi}\left(-\overline{u^{\prime} v^{\prime}} h\right) \end{align}$

$\begin{align} & \frac{\partial}{\partial t}\left(\frac{Q^\eta}{J}\right)+\frac{\partial}{\partial \xi}\left(\frac{U Q^\eta}{J}\right)+\frac{\partial}{\partial \eta}\left(\frac{V Q^\eta}{J}\right)-\frac{M}{J}\left(U \frac{\partial \eta_x}{\partial \xi}+V \frac{\partial \eta_x}{\partial \eta}\right)-\frac{N}{J}\left(U \frac{\partial \eta_y}{\partial \xi}+V \frac{\partial \eta_y}{\partial \eta}\right) \\ & =-g h\left(\frac{\xi_x \eta_x+\xi_y \eta_y}{J} \frac{\partial z_s}{\partial \xi}+\frac{\eta_x^2+\eta_y^2}{J} \frac{\partial z_s}{\partial \eta}\right)-\frac{\tau_b^\eta}{\rho J}+\frac{\xi_x \eta_x}{J} \frac{\partial}{\partial \xi}\left(-\overline{u^{\prime 2}} h\right)+\frac{\eta_x^2}{J} \frac{\partial}{\partial \eta}\left(-\overline{u^{\prime 2}} h\right) \\ & +\frac{\xi_y \eta_y}{J} \frac{\partial}{\partial \xi}\left(-\overline{v^{\prime 2}} h\right)+\frac{\eta_y^2}{J} \frac{\partial}{\partial \eta}\left(-\overline{v^{\prime 2}} h\right)+\frac{\xi_x \eta_y+\xi_y \eta_x}{J} \frac{\partial}{\partial \xi}\left(-\overline{u^{\prime} v^{\prime}} h\right)+\frac{2 \eta_x \eta_y}{J} \frac{\partial}{\partial \xi}\left(-\overline{u^{\prime} v^{\prime}} h\right) \end{align}$

ここに、

$\begin{align} Q^{\xi} &= \xi_x M+\xi_y N \\ Q^\eta &= \eta_x M+\eta_y N \\ \tau_b^{\xi} &= \xi_x \tau_{b x}+\xi_y \tau_{b y} \\ \tau_b^\eta &= \eta_x \tau_{b x}+\eta_y \tau_{b y} \end{align}$

導出は作成中

この式形は日本国内で人気があり（水理公式集例題プログラム集の影響？）、様々な場所で用いられています。

例えば、汎用ソフトウェアiRICの式形は以下のとおりです。（iRIC Software Nays2DH Solver Manual より引用。 $\xi$ 方向の運動方程式のみ。)

$\begin{align} \frac{\partial u^{\xi}}{\partial t}+u^{\xi} \frac{\partial u^{\xi}}{\partial \xi} & +u^\eta \frac{\partial u^{\xi}}{\partial \eta}+\alpha_1 u^{\xi} u^{\xi}+\alpha_2 u^{\xi} u^\eta+\alpha_3 u^\eta u^\eta= \\ & -g\left(\left(\xi_x^2+\xi_y^2\right) \frac{\partial H}{\partial \xi}+\left(\xi_x \eta_x+\xi_y \eta_y\right) \frac{\partial H}{\partial \eta}\right) \\ & -\left(C_f+\frac{1}{2} C_D a_s h_v\right) \frac{u^{\xi}}{h J} \sqrt{\left(\eta_y u^{\xi}-\xi_y u^\eta\right)^2+\left(-\eta_x u^{\xi}+\xi_x u^\eta\right)^2}+D^{\xi} \\ \end{align}$

パターン3：Delft型

Delft型と定義していますが、汎用ソフトウェアDelft3Dで使用されている式型です。

以下は、三次元（σ座標系）の式ですが、平面二次元でもほぼ同形となります。 (Delft3D-FLOW User Manualより引用)

基本的にはパターン1と同じですが、 $J$ を使用せずに展開しています。詳しく分析できていませんが、多分格子生成段階の関係でこのような形になっていると思われます。

連続式

$\begin{align} \frac{\partial \zeta}{\partial t}+\frac{1}{\sqrt{G_{\xi \xi}} \sqrt{G_{\eta \eta}}} \frac{\partial\left((d+\zeta) U \sqrt{G_{\eta \eta}}\right)}{\partial \xi}+\frac{1}{\sqrt{G_{\xi \xi}} \sqrt{G_{\eta \eta}}} \frac{\partial\left((d+\zeta) V \sqrt{G_{\xi \xi}}\right)}{\partial \eta}=(d+\zeta) Q \end{align}$

運動方程式

$\begin{align} & \frac{\partial u}{\partial t}+\frac{u}{\sqrt{G_{\xi \xi}}} \frac{\partial u}{\partial \xi}+\frac{v}{\sqrt{G_{\eta \eta}}} \frac{\partial u}{\partial \eta}+\frac{\omega}{d+\zeta} \frac{\partial u}{\partial \sigma}-\frac{v^2}{\sqrt{G_{\xi \xi}} \sqrt{G_{\eta \eta}}} \frac{\partial \sqrt{G_{\eta \eta}}}{\partial \xi} \\ \nonumber &+\frac{1}{\sqrt{G_{\xi \xi}} \sqrt{G_{\eta \eta}}} \frac{\partial \sqrt{G_{\xi \xi}}}{\partial \eta}-f v=-\frac{1}{\rho_0 \sqrt{G_{\xi \xi}}} P_{\xi}+F_{\xi} \\ \nonumber &+\frac{1}{(d+\zeta)^2} \frac{\partial}{\partial \sigma}\left(\nu_V \frac{\partial u}{\partial \sigma}\right)+M_{\xi} \\ &\frac{\partial v}{\partial t}+\frac{u}{\sqrt{G_{\xi \xi}}} \frac{\partial v}{\partial \xi}+\frac{v}{\sqrt{G_{\eta \eta}}} \frac{\partial v}{\partial \eta}+\frac{\omega}{d+\zeta} \frac{\partial v}{\partial \sigma}+\frac{u v}{\sqrt{G_{\xi \xi}} \sqrt{G_{\eta \eta}}} \frac{\partial \sqrt{G_{\eta \eta}}}{\partial \xi} \\ \nonumber &-\frac{u^2}{\sqrt{G_{\xi \xi}} \sqrt{G_{\eta \eta}}} \frac{\partial \sqrt{G_{\xi \xi}}}{\partial \eta}+f u=-\frac{1}{\rho_0 \sqrt{G_{\eta \eta}}} P_\eta+F_\eta \\ \nonumber &+\frac{1}{(d+\zeta)^2} \frac{\partial}{\partial \sigma}\left(\nu_V \frac{\partial v}{\partial \sigma}\right)+M_\eta \end{align}$

導出は作成中

この式形は劉さんのもの近いと思います。（論文1991,Study on Sediment Transport and Bed Evolution in Compound Channelsより引用。三次元。 $x_1$ 方向の運動方程式のみ。）

当然同じ系譜のこの方も同じ式形を使っています。

連続式

$\begin{align} \frac{1}{h_1 h_2 h_3}\left(\frac{\partial}{\partial x_1}\left(h_2 h_3 u_1\right)+\frac{\partial}{\partial x_2}\left(h_1 h_3 u_2\right)+\frac{\partial}{\partial x_3}\left(h_1 h_2 u_3\right)\right)=0 \end{align}$

運動方程式( $x_1$ 方向のみ)

$\begin{align} &\frac{u_1}{h_1} \frac{\partial u_1}{\partial x_1}+\frac{u_2}{h_2} \frac{\partial u_1}{\partial x_2}+\frac{u_3}{h_3} \frac{\partial u_1}{\partial x_3}+\frac{u_1 u_2}{h_1 h_2} \frac{\partial h_1}{\partial x_2}+\frac{u_1 u_3}{h_1 h_3} \frac{\partial h_1}{\partial x_3}-\frac{u_2^2}{h_1 h_2} \frac{\partial h_2}{\partial x_1}-\frac{u_3^2}{h_1 h_3} \frac{\partial h_3}{\partial x_1} \\ = & -\frac{1}{h_1} \frac{\partial}{\partial x_1}\left(\frac{p}{\rho}+\Omega\right) -\frac{1}{h_1 h_2 h_3}\left( \frac{\partial}{\partial x_1}\left(h_2 h_3 \overline{u_1^{\prime} u_1^{\prime}}\right) +\frac{\partial}{\partial x_2}\left(h_1 h_3 \overline{u_1^{\prime} u_2^{\prime}}\right) +\frac{\partial}{\partial x_3}\left(h_1 h_2 \overline{u_1^{\prime} u_3^{\prime}}\right) \right) \\ & - \frac{\overline{u_1^{\prime} u_2^{\prime}}}{h_1 h_2} \frac{\partial h_1}{\partial x_2} - \frac{\overline{u_1^{\prime} u_3^{\prime}}}{h_1 h_3} \frac{\partial h_1}{\partial x_3} + \frac{\overline{u_2^{\prime} u_2^{\prime}}}{h_1 h_2} \frac{\partial h_2}{\partial x_1} + \frac{\overline{u_3^{\prime} u_3^{\prime}}}{h_1 h_3} \frac{\partial h_3}{\partial x_1} \end{align}$