趣味で計算流砂水理

趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

数値計算とか河川工学とかプログラミングのことを書いています

MENU

河川技術者向け基礎講座 準二次元不等流計算2/4:一般断面の不等流計算

教科書の「不等流計算」は理解できるが、実務で汎用される「準二次元不等流計算」がなかなか理解し難い方に向けて、記事を書きました。 多分相当わかりやすいと思いますので、河川行政に関わる方やコンサルの方に読んで頂きたいです。

全4回の第2回目です。


本記事はGitHub、nbviewer、Colab、Marp(スライド形式)でも公開しています。

※Colabの解説記事はこちら

※Marpの解説記事はこちら

a nbviewer Open In Colab a



はじめに

  • 下図に示す実河川ような横断面形状は一般断面や自然河道断面などと呼ばれる。前回は矩形断面を対象として不等流計算を行ったが今回は一般断面を対象とする。
  • このような実務者向けの計算テクニックは水理公式集、河川砂防技術基準などに一部記述はあるのものの参考書籍が少ない。
import matplotlib.pyplot as plt
L=[0, 5,93,100,200,206,294,300]
Z=[6, 3.5, 3.5,  0,  0,  3,  3,  6]
plt.plot(L,Z, c='k')
plt.show()                # 描画

png


基礎式

一般断面の不等流計算の基礎式は次式を用いる。

\begin{align}
& \dfrac{d}{dx}\left( \frac{Q^2}{2gA^2} + H \right) = -\dfrac{\tau}{\rho g A} 
\end{align}

ここに、Q:流量、H:水位、A:河積、\tau:コントロールボリュームに作用する力とする。

前項で示した矩形断面の式と比較して以下の変更を行っている。

  • エネルギー保存則をベースとした式形としたが、一般断面の計算では摩擦損失以外の抵抗を考慮するため、運動量保存則をベースとした式形とした。それにより、右辺の生成項の式形が変更されている。
  • h+z_b (圧力項+重力項)をHに書き換えている。

後者について補足すると、 例えば、下図のように、WL1、WL2と2つの水位が与えらた場合にそれぞれh,z_bは定義できるだろうか。 それぞれの定義は難しい上にz_bは水位によって変わることが理解できる。これが、一般断面の大きな特徴であり、水位によって理論的な(物理的に意味を持つ)河床高が変わる。 h,z_bの2変数を同時に解くことはできないため、h+z_b \equiv HとしてHのみを変数として計算を行なう。

import matplotlib.pyplot as plt
x=[0, 5,93,100,200,206,294,300]
y=[6, 3.5, 3.5,  0,  0,  3,  3,  6]
plt.plot(x,y,c='k')
plt.axhline(y=5,c='b',label='WL1')
plt.axhline(y=3.3,c='b',ls='dashed',label='WL2')
plt.legend()
plt.show()                # 描画

png

マニング則についても以下のとおり若干の変更を加える。

\begin{align}
 v &= \dfrac{1}{n}i_e^{1/2}R^{2/3} \\
 Q &= \dfrac{1}{n}i_e^{1/2}\dfrac{A^{5/3}}{S^{2/3}} \\
 Q &= K i_e^{1/2} \\
 K &\equiv \dfrac{A^{5/3}}{nS^{2/3}} 
\end{align}

ここに、Q:流量、H:水位、A:河積、v:流速、R=A/S:径深、S:潤辺、K:通水能、n:マニングの粗度係数とする。

水位A、潤辺S、粗度係数nは、水位Hの関数となる。式の簡略化のため、これらを合わせた通水能Kを用いることがある。


運動方程式の生成項が摩擦損失のみの場合、マニング則を用いると以下のとおりになる。

\begin{align}
 \dfrac{\tau}{\rho g A} = \dfrac{Q^2}{K^2}
\end{align}

限界水深、等流水深の定義

一般断面の基礎式を用いて、限界水深、等流水深を矩形断面と同様に定義で設定することは難しいため、便宜的に以下のとおりに設定する。

等流水深

全ての損失による水頭の勾配が河床勾配i_bと釣り合う状態を等流と定義してその水深を等流水深とする。 なお、一般断面では水深を用いないため、正確には等流時の水位(以降、等流水位と定義)となる。

\begin{align}
 \dfrac{\tau}{\rho g A} = i_b
\end{align}

生成項が摩擦損失のみの場合は、

\begin{align}
 \dfrac{Q^2}{K^2} = i_b
\end{align}

となり、これを満足する水位Hを反復法などにより求めれば良い。


限界水深

フルード数が1となる水位を限界流時の水位(以降、限界水位と定義)とする。

フルード数は、

\begin{align}
    Fr =  \dfrac{Q}{A\sqrt{gh}}
\end{align}

となるが、平方根の中に水深hが含まれるため、そのままでは計算できない。

そのため、水深hの代替として径深R=A/Sを使用することが多い。その他にはA/Bを用いることもある。

参考:FORUM8ソフトウェア:等流・不等流の計算・3DCAD Ver.9 Q&A


離散化

離散化は次式となる。 なお、i:上流側、i-1:下流側とする。

\begin{align}
  \left(\frac{Q^2}{2gA^2_i} + H_i  \right) 
-\left( \frac{Q^2}{2gA_{i-1}^2} + H_{i-1}  \right) 
= \dfrac{1}{2}\left(\dfrac{Q^2}{K_i^2} + \dfrac{Q^2}{K_{i-1}^{2}}\right)\Delta x
\end{align}

常流の場合、下流から逐次計算を行なうため未知数はH_iのみとなる。

数値計算方法

一般断面の場合、離散式の未知数H_iによる微分が難しいため、ニュートン法が使いづらい。そのため、二分法を使用することをおすすめする(もちろんニュートン法を使っても良い)。ただし、通常の二分法では安定的に計算することが難しいため、少し工夫が必要である。この点については次回詳述する。


数値計算方法2:断面特性(河積、潤辺)の計算方法

河道断面の座標を次図のように定義する。

A


河積A,潤辺S,径深Rは次式で計算する。

\begin{align}
S_i &= \sqrt{(L_{i+1}-L_i)^2 + (Z_{i+1}-Z_i)^2} \\
A_i &= \left[H-0.5(Z_{i+1}+Z_i)\right](L_{i+1}-L_i) \\
S &= \sum_{i=1}^{imax-1} \!\! S_i \\
A &= \sum_{i=1}^{imax-1} \!\! A_i \\
R &= \dfrac{ A }{ S }\\
\end{align}

ここでは、水際位置と座標定義点が同一であるが、実計算ではその限りではないため、水面と河道断面の交点を求めて計算する必要がある。詳細は後述するプログラムを参照いただきたい。


数値計算方法3:(断面内の流速が一定と仮定した)合成粗度係数の計算方法

ここまでは河道断面内で粗度係数が一定であることを前提としたが、上図のとおり、断面内で粗度係数が異なる場合について考える。

各潤辺の粗度の影響が及ぶ範囲を下図のように考えて、各領域の流速がu_iが全て等しいと仮定する(Einsteinの方法)。その場合、断面平均流速V=u_iとなる。断面全体および各潤辺領域でマニング則が成立すると仮定すると、次式が導出される。

A

出典:椿東一郎 水理学1 pp.147


\begin{align}
n' &= \left( \dfrac{  \displaystyle \sum_{i=1}^{imax-1} S_i n_i^{3/2}}{ \displaystyle \sum_{i=1}^{imax-1} S_i } \right)^{2/3} 
\end{align}

この方法は平均流速公式1a、n'は合成粗度係数と呼ばれる。

当然ではあるが水位によって合成粗度係数は変わる。


導出方法

断面全体の流速Vが各領域の流速u_iと一致すると仮定する。(Einsteinの方法)

マニング則を用いると次式が導かれる。

\begin{align}
V &= u_i \\
\dfrac{1}{n'}\dfrac{A^{2/3}}{S^{2/3}} i_e^{1/2}&=
\dfrac{1}{n_i}\dfrac{A_i^{2/3}}{S_i^{2/3}} i_e^{1/2}\\
\dfrac{1}{n'}\dfrac{A^{2/3}}{S^{2/3}} &=
\dfrac{1}{n_i}\dfrac{A_i^{2/3}}{S_i^{2/3}}  \\
\dfrac{1}{n'^{3/2}}\dfrac{A}{S} &=
\dfrac{1}{n_i^{3/2}}\dfrac{A_i}{S_i} \\
A_i & = \dfrac{n_i^{3/2} S_i}{n'^{3/2} S }A
\end{align}

A = \displaystyle \sum A_i, S = \displaystyle \sum S_iのため、

\begin{align}
A &= \displaystyle \sum A_i \\
  &= \displaystyle \sum \dfrac{n_i^{3/2} S_i}{n'^{3/2} S }A \\
  &= \dfrac{A}{n'^{3/2} S } \displaystyle \sum n_i^{3/2} S_i \\
n' &= \left( \dfrac{  \displaystyle \sum  n_i^{3/2}S_i}{ \displaystyle \sum S_i } \right)^{2/3} 
\end{align}

となり、合成粗度が導かれる。


断面特性を計算するプログラム例

ある河川断面の座標が次図のとおりとする。

import matplotlib.pyplot as plt
L=[0, 5,93,100,200,206,294,300]
Z=[6, 3.5, 3.5,  0,  0,  3,  3,  6]
plt.plot(L,Z, c='k')
plt.show()                # 描画

png

また、粗度係数は次のとおりとする。(個数は測点間で定義するためi-1個)

n=[0.041, 0.041,0.030,0.030,0.030,0.040,0.040]

断面特性を計算するプログラムは次のように記述できる。

def H2ABSKn(l, z, n, H):
    A, B, S, SN = float(0), float(0), float(0), float(0)
    
    for i in range(1, len(l)):
        dx = l[i] - l[i-1]    
        dy = z[i] - z[i-1]    
        hb, hf = H - z[i-1], H - z[i]
        
        if hb <= float(0) :
            if hf > float(0) :
                dx_dh = dx / (hf - hb)
                B += hf * dx_dh
                A += 0.5 * hf * hf * dx_dh
                Sp = hf * np.sqrt( dx_dh * dx_dh + 1.0)
                S +=  Sp
                SN += Sp * n[i-1]**1.5
        elif hf <= float(0) :
            if hb > float(0) :
                dx_dh = dx / (hf - hb)
                B -= hb * dx_dh
                A -= 0.5 * hb * hb * dx_dh
                Sp = hb * np.sqrt(dx_dh * dx_dh + 1.0)
                S += Sp
                SN += Sp * n[i-1]**1.5
        else :
            B += dx
            A += 0.5 * dx * (hf + hb)
            Sp = np.sqrt(dx**2 + dy**2)
            S += Sp
            SN += Sp * n[i-1]**1.5
            
    if S <= float(0):
        nd = float(0)
        K = float(0)
    else:
        nd = (SN/S)**(2.0/3.0)
        K = A**(5.0/3.0)/nd/S**(2.0/3.0)
        
    return A, B, S, K, nd

  • 水位が5.0の場合、各種断面諸元は次のとおりとなる。
import numpy as np
A, B, S, K, nd = H2ABSKn(np.array(L), np.array(Z), np.array(n), 5.0)
print('河積:{}\n水面幅:{}\n潤辺:{}\n通水能:{}\n合成粗度係数:{}'.format(A, B, S, K, nd))
河積:858.0
水面幅:296.0
潤辺:298.3606797749979
通水能:47342.84520415623
合成粗度係数:0.03664910724429057

不等流計算のプログラム例

Coming soon.

参考文献

関連記事

computational-sediment-hyd.hatenablog.jp

河川技術者向け基礎講座 準二次元不等流計算1/4:不等流計算の基礎

教科書の「不等流計算」は理解できるが、実務で汎用される「準二次元不等流計算」がなかなか理解し難い方に向けて、記事を書きました。 多分相当わかりやすいと思いますので、河川行政に関わる方やコンサルの方に読んで頂きたいです。

全4回の第1回目です。


本記事はGitHub、nbviewer、Colab、Marp(スライド形式)でも公開しています。

※Colabの解説記事はこちら

※Marpの解説記事はこちら

a nbviewer Open In Colab a



基礎式

矩形近似した断面の不等流計算の基礎式はベルヌーイの式とマニング則より次式となる。

\begin{align}
& \dfrac{d}{dx}\left( \frac{q^2}{2gh^2} + h + z_b \right) = -i_e \\
& q = \dfrac{1}{n}i_e^{1/2}h^{5/3} \\
\end{align}

ここに、q:単位幅流量、h:水深、z_b:河床高、i_e:エネルギー勾配、n:マニングの粗度係数、g:重力加速度とする。


等流水深h_0、限界水深h_cは次式となる。

\begin{align}
h_0&=\left(\dfrac{q^2n^2}{i_b}\right)^{3/10} \\
h_c&=\left(\dfrac{q^2}{g}\right)^{1/3}
\end{align}

ここに、i_b = - \dfrac{d z_b}{dx}:河床勾配とする。


離散化

ベルヌーイの式とマニング則を連立させて離散化すると次式となる。 なお、i:上流側、i-1:下流側とする。

\begin{align}
\left(\frac{q^2}{2gh^2_i} + h_i + {z_b}_i \right) 
-\left( \frac{q^2}{2gh_{i-1}^2} + h_{i-1} + {z_b}_{i-1} \right) 
= \dfrac{1}{2}\left(\dfrac{q^2n^2}{h_i^{10/3}} + \dfrac{q^2n^2}{h_{i-1}^{10/3}}\right)\Delta x
\end{align}

エネルギー勾配は上下流断面の平均値を用いている。

常流の場合、下流から逐次計算を行なうため未知数はh_iのみとなる。


水面形方程式

前出のベルヌーイの式を変形して次の水面形方程式(水深の微分に関する方程式)を得る。

\begin{align}
 \dfrac{d}{dx}\left( \dfrac{q^2}{2gh^2} + h + z_b \right) &= -i_e \\
- \dfrac{q^2}{gh^3}\dfrac{dh}{dx} + \dfrac{dh}{dx} + \dfrac{d z_b}{dx} &= -i_e \\
\dfrac{dh}{dx} &= \dfrac{i_b -i_e }{ 1 - \dfrac{q^2}{gh^3}}  \\
 &= \dfrac{i_b -i_e }{ 1 - {Fr}^2} 
\end{align}

ここに、Fr = \dfrac{v}{\sqrt{gh}} = \dfrac{q}{\sqrt{gh^3}}:フルード数とする。


水面形方程式の分子=0は\dfrac{dh}{dx}=0のため等流状態を示し、前出の等流水深が得られる(等流水深の定義)。 また、分母=0から限界水深が得られ(限界水深の定義)、\dfrac{dh}{dx}=\inftyのため水面形が得られない。この地点は支配断面と呼ばれる。

上式をh_0,h_cを用いて展開すると次式が得られる。(式展開は例えばここを参照)

\begin{align}
\dfrac{dh}{dx} &= i_b \dfrac{1-\left(\dfrac{h_0}{h}\right)^{10/3}}{1-\left(\dfrac{h_c}{h}\right)^3 }
\end{align}

不等流計算を行なう上で本式より以下を理解しておく必要がある。


水面形の追跡方向

  • Fr \lt 1またはh_c \lt hまたはv\lt\sqrt{gh}は、常流
  • Fr \gt 1またはh_c \gt hまたはv\gt\sqrt{gh}は、射流

と定義される。 \sqrt{gh}は(微小振幅波理論による水面波の)波速であり、射流はv\gt\sqrt{gh}のため水面変化が下流のみに影響するが、常流では上流にも影響をすることが示されている。 そのため、不等流計算の水面形の追跡方向は、常流は下流から上流、射流は上流から下流となる。


水面形の概形

h,h_0,h_cの関係より\dfrac{dh}{dx}の符号が決まるため水面形の概形が把握できる。教科書にみられる緩勾配水路(M水路、h_0 \gt h_c)や急勾配水路(S水路、h_0 \lt h_c)などの水面形はこの考え方による描くことができる。

A

出典:椿東一郎 水理学1 pp.151


数値計算方法

未知数h_iは直接計算できなため反復法によって近似解を求める。ここでは、ニュートン法を用いた。(参考:ニュートン法とは

ニュートン法アルゴリズムに沿うと、計算式は次のとおりとなる。

\begin{align}
f &= \left(\frac{q^2}{2gh^2_i} + h_i + {z_b}_i \right) 
-\left( \frac{q^2}{2gh_{i-1}^2} + h_{i-1} + {z_b}_{i-1} \right) 
 - \dfrac{1}{2}\left(\dfrac{q^2n^2}{h_i^{10/3}} + \dfrac{q^2n^2}{h_{i-1}^{10/3}}\right)\Delta x \\
\dfrac{df}{dh_i} &= -\frac{q^2}{gh^3_i} + 1 
+ \dfrac{5}{3}\dfrac{q^2n^2}{h_i^{13/3}} \Delta x
\end{align}
\begin{align}
h^{new}_i &= h_i - \dfrac{f}{\dfrac{df}{dh_i}}
% f &= \left(\frac{q^2}{2gh^2_i} + h_i + {z_b}_i \right) 
% -\left( \frac{q^2}{2gh_{i-1}^2} + h_{i-1} + {z_b}_{i-1} \right) 
%  + \dfrac{1}{2}\left(\dfrac{q^2n^2}{h_i^{10/3}} + \dfrac{q^2n^2}{h_{i-1}^{10/3}}\right)\Delta x \\
% \dfrac{df}{dh_i} &= -\frac{q^2}{gh^3_i} + 1 
% - \dfrac{5}{3}\dfrac{q^2n^2}{h_i^{13/3}} \Delta x
\end{align}

サンプルコード

単位幅流量0.5\mathrm{m^3/s/m}、マニングの粗度係数0.02、河床勾配1/500、水路延長50\mathrm{m}の矩形水路を対象に下流端水位を限界水深としたときの水面形を計算せよ。なお、\Delta xは0.1\mathrm{m}とする。


import numpy as np

q = 0.5
n = 0.02
ib = 1/500
g = 9.8
dx = 0.1

h0 = (q**2*n**2/ib)**0.3 #等流水深
hc = (q**2/g)**(1/3) # 限界水深

L = np.arange(0,50.01,dx) #追加距離の配列
zb = L*ib #河床高の配列
h = np.zeros_like(L) #水深の配列

h[0] = hc #下流端条件
for i in range(1,len(h)):
    h[i] = h[i-1] #収束計算の初期値:一つ下流側の断面の水深
    f = 1.0 #仮値
    dfdh = 1.0 #仮値
    while np.abs(f/dfdh) > 10**(-8): # 反復計算の収束条件
        f = q**2/2.0/g/h[i]**2 + h[i] + zb[i] \
          -(q**2/2.0/g/h[i-1]**2 + h[i-1] + zb[i-1]) \
          - 0.5*(q**2*n**2/h[i]**(10/3) + q**2*n**2/h[i-1]**(10/3))*dx
        dfdh = -q**2/g/h[i]**3 + 1 + 5/3*q**2*n**2/h[i]**(13/3)*dx
        h[i] -= f/dfdh

import matplotlib.pyplot as plt

plt.plot(L,zb, label='zb')
plt.plot(L,h+zb, label='water level')
plt.plot(L,zb+h0, label='h0')
plt.plot(L,zb+hc, label='hc')
plt.legend()              # 凡例の表示
plt.show()                # 描画


補足:常流、射流混在流れの数値計算方法

実河川では射流が出現する範囲は局所的であり、射流の水面形が求めれることはそう多くない。

常流、射流の混在流れの計算方法については例えば以下が参考となる。

かなり煩雑な計算が必要な上、わずかな条件の差で跳水位置が変わるなど実務的に取り扱いが難しい。 また、射流部では\Delta xを十分に小さくとることが必要となる。


実務では、常流区間の中で一部が射流であれば限界水深に置き換えて計算する場合が多い。また、射流部のフルード数が大きい場合は、次に示すようなリミッターを設定することにより水面形を滑らかにすることもある。

計算負荷が大きくても問題がない場合は不定流計算による定常解により水面形を与える方法も考えられる。

常流、射流混在流れの水面形は計算手法によって大きく異なってしまうため、計算の際は十分に検討することが必要である。

参考文献

関連記事

computational-sediment-hyd.hatenablog.jp

(未解決)測量の高解像度を意識して、浅水流一次元計算でdxをかなり小さくすると上手くいかない件について

今年中に何とかしたかったのですが出来そうに案件です。

浅水流の計算(1D,2D共通)で、メッシュサイズを小さくしすぎると、ある閾値以下になると急に水位が大きくなると問題について、簡易なテストケースにより分析してみました。原因は未解明のままです。

放置しておくと厄介な問題になりそうな気がしています。


TL;DR
  • dxをどんどん小さくして0.5mとした場合に不定流計算の結果が怪しいっぽい
  • その原因は重力項?河床高+水深⇒水位として扱う影響か?
  • この問題は河床変動計算を行う場合、影響が大きいのでは?


モチベーション

近年、直轄管理河川では、河床の面的測量を順次実施されており、0.5~1.0m間隔の河床地形が得られる状況が整いつつあります。 それに合わせて、計算格子をより小さくする必要があると考えられます。

そこで、平面二次元計算でメッシュサイズを小さくしていくとある閾値で急激に水位が上昇する結果が出たため、詳細について一次元で検証してみました。

河床縦断形状のモデル化

一次元計算のためのモデル地形を生成しました。 1/400の河床勾配に基本に、縦断方向に0.5m間隔の格子点を設けて、各点に-10~10cmの乱数を加算した以下の形状を想定しました。

  • 全体図
  • 拡大図
  • 河床地形の生成 aOpen In Colab

  • 図化 a Open In Colab

dxの差異による一次元計算結果の比較

計算方法は過去記事と同様に、不等流計算、不定流計算(コロケート格子)、不定流計算(スタッガード格子)の3手法でテストしました。

3手法ともにdx=0.5、10、100mの3ケースの比較を行っています。dx=10、100mの河床地形は前述の0.5m間隔の格子点より縦断距離が一致する地点を抽出することにより作成しました。

なお、不定流計算では定常となるまで(実時間で3時間程度)繰り返し計算を行い、縦断的に流量が一定であることを確認している。

computational-sediment-hyd.hatenablog.jp

不等流計算

dx=0.5、10、100mの3ケースの平均的な水位は同程度であり、dxが小さいほど河床地形に応答して細かい水位変化が生じていることが分かる。

  • 全体図
  • 拡大図

a Open In Colab

不定流計算(コロケート格子)

dx=10、100mの2ケースと比較して、dx=0.5mの平均水位が0.2m程度大きくなっている。

  • 全体図
  • 拡大図

a Open In Colab

不定流計算(スタッガード格子)

dx=10、100mの2ケースと比較して、dx=0.5mの平均水位が0.1m程度大きくなっている。

  • 全体図
  • 拡大図

a Open In Colab

各計算手法の差異

dx=0.5mにおける3計算手法の水位を比較した。

  • 水位は不等流<不定流:スタッガード<不定流:コロケートの順で大きくなっている。前項の結果を踏まえると、不定流の結果が真値からずれていると想定される。
  • 不等流は、河床の凹凸に逆位相となるような常流の特徴を示しているが、不定流ではその特徴が見られない。
  • 不定流:スタッガードは、河床の凹凸に対して訛ったような水面形を示している。
  • 不定流:コロケートは、不等流と比較して、水面形が若干下流側に位相がずれるような形状となっている。

  • 全体図

  • 拡大図
  • 計算 a Open In Colab
  • 図化 a Open In Colab

まとめ

  • 不定流計算では計算格子が小さくなると計算誤差が大きくなる傾向にある。一方、不等流計算ではその傾向は見られない。
  • その要因は不明であるが、重力項(河床勾配)の影響ではないかと思われる。
  • 一般的にスキームチェックを行う場合、河床勾配はフラットで行うことが多いが、本記事のような急変地形でのテストも必要と考える。

GitHub

github.com

参考書籍

関連記事

  • グラフの作成はHoloviewsを使用

computational-sediment-hyd.hatenablog.jp

  • はてブロへのhtmlグラフの埋め込み

computational-sediment-hyd.hatenablog.jp

computational-sediment-hyd.hatenablog.jp

祝:いつの間にかColabの標準ライブラリにholoviewsが追加されていました(2023-12-30)

いつの間にかGoogle Colaboratoryの標準ライブラリにholoviewsが追加されていました。

これでさらにholoviewsが使いやすいなります。

computational-sediment-hyd.hatenablog.jp

computational-sediment-hyd.hatenablog.jp

備忘録:Pythonでどうしようもないテキストファイルの読み込み方

データを扱っていると、どうしようもないテキストファイルってあると思いますが、それを読み込むときのテンプレートです。

# とりあえずエンコードを指定して開く
with open(filepath, mode='r', encoding='shiftJIS') as response:
 #   skip header n回
    for _ in range(n) : l = next(response) 
    
    #try文でファイルの終わりまで(errorがでるまで)ループ
    try:
        while True:
            #一行読んで改行を置換    
            l = next(response) 
            ll = l.replace('\n','')
            
            # 区切り文字で区切る場合
            xs = ll.split(',')
            # 固定長の場合:lstrip(),rstrip()もあり
            x = ll[2:10].strip()
           
            # 何らかの条件で読み込みを終了する場合
            if (条件) : break
                
    except StopIteration:
        pass

境界適合格子による平面二次元河川流計算の式形についての考察(途中段階)

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



はじめに

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

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

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

  1. 通常型
  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と比べて、項数が増えて複雑な式系になります。

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

\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}
\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}

ここに、h_1, h_2, h_3はメトリックス係数とする。

パターンの比較検討(途中段階)

正直メリット、デメリットは全然わからないです。。。もう少し考えます。

  • パターン1, パターン3は移流項の風上化などが難しいように思われる。
  • パターン2は乱流モデルを複雑にした場合などは取り扱いが大変。
  • 非保存形での記述の場合、特に陽解法では離散化方法に注意が必要。
  • ヤコビアン等の座標変換に関する変数の空間微分は中央差分?

参考サイトなど