過去記事「（未解決）測量の高解像度を意識して、浅水流一次元計算でdxをかなり小さくすると上手くいかない件について」のスタッガード格子の問題の解決

反響の多かった過去記事：（未解決）測量の高解像度を意識して、浅水流一次元計算でdxをかなり小さくすると上手くいかない件について - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learningについて、スタッガード格子の離散化の問題点を修正したため投稿します。

課題：dxが小さい場合に一次元不定流計算の結果がおかしい

過去記事で以下の課題が見られました。

一次元浅水流計算でdxを小さくしていくと、dx=0.5mの場合、計算水位が相対的に大きくなる。（下図、再掲）

dxが小さい場合に不等流計算と不定流計算に差異が生じる。（下図、再掲）

この課題のうち、スタッガードスキームについて解決したため、以降にまとめます。

一般断面の一次元不定流計算の基礎式

一般断面の一次元不定流計算の基礎式は以下となります。(上：質量保存則（連続式）、下：運動量保存則)

$\begin{align} &\dfrac{\partial A}{\partial t}+\dfrac{\partial Q}{\partial x} = 0 \\ &\frac{\partial Q}{\partial t} + \frac{\partial }{\partial x}\left( \frac{\beta Q^2}{A}\right)+gA\frac{\partial H}{\partial x}+gAI_e = 0 \\ &I_e = \dfrac{T}{\rho g A}\\ \end{align}$

$I_e$ が河床の摩擦損失のみの場合、マニング則を用いると以下となります。

$\begin{align} I_e = \dfrac{n^2 V|V|}{ R^{4/3}} \end{align}$

一般断面の場合、 $A \neq B \cdot h$ 、 $H \neq h + z_b$ のため、上記のような式形となります。

基礎式のスタッガードスキームによる離散化

前回記事

前回の記事では、水理公式集例題プログラム集(13年版)【例題 2-9】常流・射流混在流れのコードを参考に以下の離散式を用いました。

リンク

$\begin{align} &\frac{A_{i}^{n+1}-A_{i}^n}{\Delta t}+\frac{Q_{i+1/2}^n-Q_{i-1/2}^n}{\Delta x}=0 \\ &\frac{Q_{i+1/2}^{n+1}-Q_{i+1/2}^n}{\Delta t}+\frac{u_{i+1}^n Q_{i+1/2+a}^n-u_{i}^n Q_{i-1/2+b}^n}{\Delta x}+g A_i^{n+1} \frac{H_{i+1}^{n+1}-H_{i}^{n+1}}{\Delta x}=-\left(g A_{i+1/2}^{n+1} \dfrac{n^2 V_{i+1/2}^n|V_{i+1/2}^n|}{ {R_{i+1/2}^{n+1}}^{4/3}}\right) \end{align}$

$\begin{align} a=\left\{\begin{array}{l}0 ; u_{i+1}^n \geq 0 \\ 1 ; u_{i+1}^n<0\end{array} , b=\left\{\begin{array}{l}0 ; u_{i}^n \geq 0 \\ 1 ; u_{i}^n<0\end{array} , u_{i}^n=\frac{u_{i+1/2}^n+u_{i-1/2}^n}{2} , u_{i+1/2}^n=\frac{Q_{i+1/2}^n}{A_{i+c}^{n+1}} , c=\left\{\begin{array}{l}0 ; Q_{i+1/2}^{n} \geq 0 \\ 1 ; Q_{i+1/2}^n<0\end{array} \right.\right.\right. \end{align}$

本記事の修正点

今回の修正した点は、最後の式の $u_{i}^n=\dfrac{u_{i+1/2}^n+u_{i-1/2}^n}{2},u_{i+1/2}^n=\dfrac{Q_{i+1/2}^n}{A_{i+c}^n}$ で示すハーフセルの流速の計算方法です。
風上化は行わずに中心差分に変更しました。修正した式は以下の通りとなります。

$\begin{align} a=\left\{\begin{array}{l}0 ; u_{i+1}^n \geq 0 \\ 1 ; u_{i+1}^n<0\end{array} , b=\left\{\begin{array}{l}0 ; u_{i}^n \geq 0 \\ 1 ; u_{i}^n<0\end{array} , u_i^n=\frac{Q_{i+1/2}^n + Q_{i-1/2}^n}{2A_{i}^{n+1}} \right.\right. \end{align}$