2020-11-03

Re:川幅

以下の前記事への意見です。

computational-sediment-hyd.hatenablog.jp

素晴らしい記事ありがとうございます。勉強になります。

$\frac{\partial B}{\partial x}$ の影響がかなり効きそうですね。非保存形（一般的な非保存形と混合しそう）の方が流量が保存されるというとカオス感が。。。

意見と今後の展開を書きます。

従来の静水の問題

従来の静水の問題、つまり、 $\frac{\partial H}{\partial x}=0$ で流れが静止している状況が表現できるかどうかの条件が保存形で書くとさらに複雑になります。

前記事の式を再記する。（この表現海外の論文でよく見るとけど、 $I$ の次元がバラバラなので個人的にはあまり好きじゃない。。。）

$\begin{align} 非保存形&: \dfrac{\partial Q}{\partial t} + \dfrac{\partial }{\partial x}\left(\dfrac{Q^2}{A}\right) = -gA\left(\dfrac{\partial H}{\partial x} + I_e\right) \\ 保存形&: \dfrac{\partial Q}{\partial t} + \dfrac{\partial }{\partial x} \left(\dfrac{Q^2}{A}+gI_1 \right) = -gA\left(\dfrac{\partial z}{\partial x} + I_e\right)+gI_2 \end{align}$

流速0とすると、

$\begin{align} 非保存形&: \dfrac{\partial Q}{\partial t} = -gA\dfrac{\partial H}{\partial x} \\ 保存形&: \dfrac{\partial Q}{\partial t} + g\dfrac{\partial I_1}{\partial x} = -gA\dfrac{\partial z}{\partial x} + gI_2 \end{align}$

となり、保存形において、 $\frac{\partial H}{\partial x}=0$ で流れが生じないためには、次式を満足する必要があります。

$\begin{align} g \dfrac{\partial I_1}{\partial x} = -gA\dfrac{\partial z}{\partial x} + gI_2 \end{align}$

この問題は非常に重要で、特に浅水流近似の平面二次元流れの場合（当然、 $I_2$ の項はない）は支配的になります。

そう考えていくと、結論としては浅水流近似では保存形の式は使えないことになります。

今後の展開

今後の展開としては、非保存形のままシステム形にする方法の可能性を検討するのはいかがでしょう。

と、そんなことを考えながらググっていたら、up-wind flux schemeのYingさんが非保存形のリーマンソルバーをやってました。

https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwjH7M2e6OPsAhWHdXAKHcfJC70QFjAQegQIChAC&url=http%3A%2F%2Fyadda.icm.edu.pl%2Fyadda%2Felement%2Fbwmeta1.element.baztech-article-BAT8-0017-0010%2Fc%2FYing.pdf&usg=AOvVaw0IOeWUNgAG_Eq8jz5E1_1R

時間を見つけてコーディングしてみます。

本業が忙しい。
大阪都構想否決が超落ち込む。

2020-11-01

川幅

久々に。。。

モチベーション

ひょんなことからTVD-MacCormack法を実装しようと思い、山地河川における河床変動の数値計算法（砂防学会の「急流河川の１次元河床変動（その２）」）を見て実装しようと思ったのですが、川幅に変化がある場合の計算が上手くいかなことを思い出したので少しまとめておきます。

基礎方程式

プリミティブ変数を通水断面積 $A$ と流量 $Q$ にした場合、連続式は次のようになるかと思います。

$\dfrac{\partial A}{\partial t} + \dfrac{\partial Q}{\partial x} = 0$

一方で運動方程式は、水面勾配を用いる場合と、河床勾配を用いる場合があります。ただし数式の表現は違いますが、同じものになります。前者と後者はそれぞれ次のようになります。

$\dfrac{\partial Q}{\partial t} + \dfrac{\partial }{\partial x}\left(\dfrac{Q^2}{A}\right) = -gA\left(\dfrac{\partial H}{\partial x} + I_e\right)$ ]

$\dfrac{\partial Q}{\partial t} + \dfrac{\partial }{\partial x} \left(\dfrac{Q^2}{A}+gI_1 \right) = -gA\left(\dfrac{\partial z}{\partial x} + I_e\right)+gI_2$

実河川への適用を考えると、前者の方が使いやすいですが、後者は所謂、保存系であり、システム系の方程式の解法（FDS、FVS、HLL系等）が使えます。この記事では前者を非保存系、後者を保存系と呼ぶことにします。何度も言いますが、数学的にこの二つの式は等価です。ただし離散化して解くと同じ答えにならない場合があります。

前者と後者を比較すると、後者では $I_1$ と $I_2$ が新たに現れています。 $gI_1$ は断面に働く静水圧、 $gI_2$ は川幅の変化とかサイドスロープの勾配の変化を表現する項になります。

多くの論文は単位幅だったり一様水路だったりするので、 $I_1$ は表現されているものの、 $I_2$ は矩形であれば0となるのでネグってることが多いです。次に数式でみてみます。

$I_1$ と $I_2$ の数式による表現

それぞれ表現すると次のようになります。

$gI_1 = g\int_0^{h}\left(h-\eta\right)bd\eta$

$gI_2 =g\int_0^{h}\left(h-\eta\right)\dfrac{\partial b}{\partial x}d\eta$

$b$ は流下距離 $x$ と河床から水面までの距離 $\eta$ の関数です。

$I_1$ は実河川でも求められますが、 $I_2$ はかなりめんどくさいです。ここでは矩形断面を想定し、河床から水面まで幅は一定( $=B$ )とします。

$gI_1 = gB\int_0^{h}\left(h-\eta\right)d\eta=gB\dfrac{h^2}{2}$

$gI_2 = g\dfrac{\partial B}{\partial x}\int_0^{h}\left(h-\eta\right)d\eta =g\dfrac{\partial B}{\partial x} \dfrac{h^2}{2}$

$I_2$ は一様水路の場合は0 $\left(\partial B/\partial x=0\right)$ になります。

ちなみに台形断面の場合は $I_2$ は比較的簡単に求められます。河床における川幅を $B_b$ 、サイドスロープの勾配を $S_L$ とすると $b(x,\eta)$ は次のように表現できます。

$b(x,\eta) = B_b + 2S_L \eta$

$gI_1 =g\int_0^{h}\left(h-\eta\right) \left(B_b + 2S_L \eta\right)d\eta = gB\dfrac{h^2}{2}+\dfrac{S_L}{3} h^3$

$gI_2 =g\int_0^{h}\left(h-\eta\right)\dfrac{\partial }{\partial x}\left(\dfrac{\partial B_b}{\partial x} + 2 \eta \dfrac{\partial S_L}{\partial x}\right)d\eta=g\dfrac{\partial B}{\partial x} \dfrac{h^2}{2} + \dfrac{h^3}{3} \dfrac{\partial S_L}{\partial x}$

「急流河川の１次元河床変動（その２）」の事例を記載の保存系の方程式で解く場合、 $I_2$ を考慮していないため再現できません。久々に解いたので忘れていたのもありますが、ここら辺の書き方がモヤッとしているため少し整理しようというのがモチベーションにつながっています。

保存系のシミュレーション

前述の「急流河川の１次現河床変動（その２）」の再現をしてみます。矩形の非一様断面の流れの解析です。運動方程式は次のようになります。

$\dfrac{\partial Q}{\partial t} + \dfrac{\partial }{\partial x} \left(\dfrac{Q^2}{A}+gB\dfrac{h^2}{2}\right) = -gA\left(\dfrac{\partial z}{\partial x} + I_e\right)+g\dfrac{\partial B}{\partial x} \dfrac{h^2}{2}$

連続式を再掲すると次の通りです。

$\dfrac{\partial A}{\partial t} + \dfrac{\partial Q}{\partial x} = 0$

これをTVD-MacCormack法で解いてみます。河川の解析では、Causon型が好まれているようですが、今回はGarciaのTVDを使ってみました。なお、Causonの原著を読んでいないので詳しくはわからないのですが、どうも任意に設定できるパラメータがあるようです。日本の河川の論文ではそのパラメータを0.5として表現している模様です。

え、なぜMacCormack法か？これは他のスキームと比較して明確でない！理由があります。

集中格子かつ等間隔のため、右辺の１階微分項を中央差分で解いた場合、[tex:(z{i+1} - z{i-1})/2\Delta x]となり $z_i$ が使われず気持ち悪い。
2Stepの場合、前進と後進をそれぞれ実施するため、上記の課題が克服できてるのではないか。

真面目に考えていないですが、こんな感じです。

TVD-MacCormack法

TVD-MacCormack法は連続式と運動方程式をベクトル表示して表現すると次のとおりです。

$\dfrac{\partial U}{\partial t} + \dfrac{\partial E}{\partial x} = S$

$U = \left(A,Q\right)^T$

$E = \left(Q,Q^2/A + gBh^2/2\right)$

$S = \left(0,-g\left(\partial z / \partial x +I_e + gh^2/2(\partial B /\partial x)\right)\right)$

GarciaのTVD-MacCormack法はCausonとは前進、後進が逆のようです。表現すると次のとおりです。

$U_i^{pred} = U_i - \dfrac{\Delta t}{\Delta x}\left(E_{i+1}^n - E_{i}^{n}\right) + \Delta t S_{i}^n$

$U_i^{corr} = U_i - \dfrac{\Delta t}{\Delta x}\left(E_{i+1}^{pred} - E_{i}^{pred}\right) + \Delta t S_{i}^{pred}$

$U_i^{ave} = \dfrac{U_i^{pred} + U_i^{corr}}{2}$

$U_i^{n+1} =U_i^{ave} + \dfrac{\Delta t}{\Delta x}\left(D_{i+1/2}^n-D_{i-1/2}^n\right)$

4行目の右辺第２項がTVD項になります。注意点は生成項 $S$ です。 $x$ に関する1階微分の項は予測子、修正子に合わせて前進、後進を使い分けています。

$D$ の求め方については書くのがだるいので、次を参考にしてください。 Roeの平均を使ってます。

Shallow Water Hydraulics

作者:Castro-Orgaz, Oscar,Hager, Willi H.
発売日: 2020/12/04
メディア: ペーパーバック

教えていただいたnumbaを使いました。当初は普通にfor文を使わずに書いてましたが、遅すぎたのでFortranチックにわざわざ書き換えてnumbaで高速化しました。

川幅の変化を考慮した結果

結果はこんな感じです。 Q=0.002m3/s、ks=0.018、下流端0.058mのケースを解いています。

f:id:SedimentHydraulics:20201101234856j:plain

$I_2$ の項を考慮する、しないを比較してみました。当たり前といえば当たり前ですが・・・。

急拡、急縮の前後で水位を正しく表現できません。また幅を考慮していない場合、フルード数は1を超えています。元記事では水位のグラフのみしか示されていないのでなんとも言えませんが。

非保存系で解いた場合

非保存系でも解いてみました。保存系では集中格子で解いたので、千鳥格子にはせず、UpwindFluxSchemeを使いました。

TVD-MacCormackと似たような結果になります。

流量を比較すると保存系で解く場合はあまりよろしくありません。

結論

素直に非保存系を基礎方程式として解くのが良いのではないでしょうか？

論文等で１次元の洪水や河床変動解析においてFDSやHLL系等の保存系の基礎方程式を使って解いている場合、 $gI_2$ の川幅の変化をどう取り扱っているのかを確認した方が良いです。もちろん、現場のデータでは格子間隔が大きく（数十m〜200m程度）、川幅が大きく変わるようなデータにならない可能性はあり、川幅の影響を考慮しなくてもそれなりの結果となる場合もあるでしょうが。

あとは河床勾配の評価は難しいですね。

$I_2$ については、こうやって近似しよう、みたいな論文はみたことがありますが、実河川のデータでどの程度の近似精度があるのかはわかりません。

以上から、１次元の計算であれば、煩雑なスキームを使わずシンプルに解ける方が良いとの結論です。

2020-10-09

Rouse分布とLane-Kalinske分布を比較する

Rouse分布の近似式であるLane-Kalinske分布って結構荒っぽくないかなと昔から気になっていたので検証してみました。

ついでに式の導出も手書きしか無かったのでTeXにまとめておきました。

土砂濃度の鉛直分布の導出
2式の比較
- 濃度分布の比較
- 水深平均濃度の比較
  - 水深平均濃度式の導出
  - 比較
まとめ
参考文献

土砂濃度の鉛直分布の導出

基本的な考え方

浮遊砂の濃度 $C$ の移流拡散方程式は次のようになる。

$\begin{align} \dfrac{\partial C}{\partial t} + u\dfrac{\partial C}{\partial x} + v\dfrac{\partial C}{\partial y} + w\dfrac{\partial C}{\partial z} = \dfrac{\partial }{\partial x} \left( \epsilon_{sx} \dfrac{\partial C}{\partial x} \right) + \dfrac{\partial }{\partial y} \left( \epsilon_{sy} \dfrac{\partial C}{\partial y} \right) + \dfrac{\partial }{\partial z} \left( \epsilon_{sz} \dfrac{\partial C}{\partial z} \right) + w_0 \dfrac{\partial C}{\partial z} \end{align}$

定常、等流、平衡場を考えて、 $\partial C /\partial x = \partial C /\partial y = 0$ と $v=w=0$ を与えると上式は以下のようになる。

$\begin{align} \dfrac{\partial }{\partial z} \left( \epsilon_{sz} \dfrac{\partial C}{\partial z} \right) + w_0 \dfrac{\partial C}{\partial z} = 0 \end{align}$

上式を河床から水面まで積分する。両者の条件より積分定数は0となる（詳細は省略）。

$\begin{align} \epsilon_{sz} \dfrac{\partial C}{\partial z} + w_0 C = 0 \end{align}$

さらに上式を変数分離形にして基準点の高さ $a$ から水面まで積分する。

$\begin{align} \int \dfrac{d C}{C} dz = \int \dfrac{-w_0}{\epsilon_{sz}} dz \end{align}$

基準点の高さ $a$ の濃度を $C_a$ として整理すると次式が導かれる。

$\begin{align} \dfrac{C}{C_a} = \mathrm{exp} \left( - \int^z_a \dfrac{w_0}{\epsilon_{sz}} dz \right) \end{align}$

Rouse分布

せん断力の河床から水面までの分布は次式となる。

$\begin{align} \tau = \tau_0 \left(1-\dfrac{z}{h}\right) \end{align}$

乱流のせん断応力は次式で示される。( $\nu \gg \epsilon_{z}$ のため、 $\nu$ は無視する)

$\begin{align} \tau = \rho \epsilon_{z} \dfrac{d u}{d z} \end{align}$

2式の釣り合いより

$\begin{align} \epsilon_{z} = \dfrac{ u_*^2 \left(1-\dfrac{z}{h}\right)}{\dfrac{d u}{d z}} \end{align}$

対数則の微分（対数則の元）は、

$\begin{align} \dfrac{du}{dz} = \dfrac{u_*}{\kappa z} \end{align}$

のため、 $\epsilon_{z}$ は、

$\begin{align} \epsilon_{z} = u_* \kappa z \left(1-\dfrac{z}{h}\right) \end{align}$

となる。 $\epsilon_{sz}$ は $\epsilon_{z}$ に比例するとして、 $\epsilon_{sz}=\beta \epsilon_{z}$ と定義し、前述の濃度分布式に代入したものがRouse分布である。

$\begin{align} \dfrac{C}{C_a} &= \left( \dfrac{h-z}{z} \dfrac{a}{h-a} \right)^{\dfrac{w_0}{\beta \kappa u_*}} % Z &= \dfrac{w_0}{\beta \kappa u_*} \end{align}$

ここに、 $\kappa$ :カルマン定数, $a$ :基準面の河床から高さ( $=0.05h$ ), $\beta$ :土砂の拡散係数と水の渦動粘性係数の比（ $\epsilon_s/\epsilon$ ）である。

参考：式展開

$\begin{align} \dfrac{C}{C_a} = \mathrm{exp} \left( - \int^z_a \dfrac{w_0}{ \beta u_* \kappa z \left(1-\dfrac{z}{h}\right)} dz \right) \end{align}$

$t = \dfrac{h-z}{z}, \dfrac{dt}{dx} = - \dfrac{h}{z^2}$ として、置換積分を行う。

$\begin{align} \dfrac{C}{C_a} &= \mathrm{exp} \left( \dfrac{w_0}{ \beta u_* \kappa } \int^z_a \dfrac{1}{t} dt \right) \\ &= \mathrm{exp} \left( \dfrac{w_0}{ \beta u_* \kappa } \left[\log_e{ \left( \dfrac{h-z}{z} \right)} \right]^z_a \right) \\ &= \left( \dfrac{h-z}{z} \dfrac{a}{h-a} \right)^{\dfrac{w_0}{\beta \kappa u_*}} \end{align}$

Lane-Kalinske分布

Rouse分布で用いた次の $\epsilon_{z}$ の代用として水深平均した $\overline{\epsilon_{z}}$ を用いる。

$\begin{align} \epsilon_{z} = u_* \kappa z \left(1-\dfrac{z}{h}\right) \end{align}$

$\overline{\epsilon_{z}}$ は次のようになる。

$\begin{align} \overline{\epsilon_{z}} = \dfrac{1}{h} \int^h_0 u_* \kappa z \left(1-\dfrac{z}{h}\right) dz = \dfrac{1}{6} \kappa u_* h \end{align}$

この式を用いてRouse分布と同様に式展開を行うと次のLane-Kalinske分布が得られる。

$\begin{align} \dfrac{C}{C_a} &= \mathrm{exp} \left( -6\dfrac{w_0}{\beta \kappa u_*} \dfrac{z-a}{h} \right) % Z &= \dfrac{w_0}{\beta \kappa u_*} \end{align}$

本式はRouse分布の近似式といえる。

2式の比較

濃度分布の比較

以下の条件下で2式の比較を行う。

河道条件：水深5.0m、河床勾配1/1000、粒径は0.1m、0.5mm、2mmの3ケース
基準面濃度式:板倉・岸の式
沈降速度式：Rubeyの式

両式を図化すると以下のようになる。

Lane-Kalinske分布はあくまでも近似式なのでRouse分布が基本となる。 Lane-Kalinske分布は次項に示す水深平均濃度の算出に数値計算ではよく使用されるが、 Rouse分布との差は、特に鉛直方向に分布を持つ場合に結構大きいと思う。

水深平均濃度の比較

水深平均濃度式の導出

Rouse分布は積分できないため、Lane-Kalinske分布を使用して水深平均濃度の計算式を導出する。Lane-Kalinske分布による水深平均濃度は次のとおりとなる。

$\begin{align} \overline{C} &= \dfrac{1}{h} \int^h_{0} C \, dz \\ &= \dfrac{1}{h} \int^h_{0} C_a \mathrm{exp} \left( -6\dfrac{w_0}{\beta \kappa u_*} \dfrac{z-a}{h} \right) \, dz \\ &= \dfrac{C_a}{h} \dfrac{\beta \kappa u_* h}{6 w_0} \left[1 - \mathrm{exp}\left(-\dfrac{6 w_0}{\beta \kappa u_* h}h\right) \right] \\ \end{align}$

さらに、

$\begin{align} \epsilon_s &= \beta \epsilon \\ \epsilon &= \dfrac{\kappa u_* h}{6} \end{align}$

として式変形を行うと、次のとおりよく使用される式形になる。

$\begin{align} \overline{C} &= C_a \dfrac{\epsilon_s}{h w_0} \left[1 - \mathrm{exp}\left(-\dfrac{h w_0}{\epsilon_s} \right) \right] \end{align}$

河川流は浅水流近似で取り扱うことが多いため本式は多用される。

比較

以下の3つの方法による水深平均濃度の比較を行う。

上式Lane-Kalinske分布の積分式
Lane-Kalinske分布の数値積分
Rouse分布の数値積分

計算条件は前述の濃度分布の比較と同様とする。また、数値積分は相対水深の5%から100%まで95分割して計算する。

計算したは、水深平均土砂濃度(ppm)は次のとおりである。

粒径	Lane-Kalinske:水深平均	Lane-Kalinske:数値積分	Rouse:数値積分
0.1mm	25375.833528263323	24489.938287379184	23800.833750921876
0.5mm	670.4528478433683	682.3398848314391	513.0807965609013
2mm	11.852228352927673	12.45126260315744	6.631464513654046

やっぱり結構違う。

数値積分は水深の5%から100％の範囲で積分、水深平均式は河床からのため、Lane-Kalinskeの比較ではその差が生じている。
Rouse分布との差は非常に大きい。特に濃度勾配の大きい2mmの場合は決定的な差異が生じている。

まとめ

これからはRouseの数値積分を使おう。

参考文献

河村三郎：土砂水理学1 pp.271-

リンク

グラフの作成はHoloviewsを使用

computational-sediment-hyd.hatenablog.jp

はてブロへのhtmlグラフの埋め込みは以下を参照

computational-sediment-hyd.hatenablog.jp

2020-10-09

jupyterのTeXをはてブロのmarkdownに変換するpythonスクリプト:2024/1/3時点

jupyterでTeXを書くことが多いのですが、はてなブログ（markdown記法）に貼ると上手く行かないので変換用のスクリプトを書いておきました。

多分awkとかperlでもできるけど苦手なのでpythonで。

はてブロ用のTeXに対応するため、以下の処理を行う。 ※数式はalign環境のみ対応

コードブロック

コードブロックは処理から除く

align環境

数式はalign環境のみ対応。文中$$にも対応。
align環境は、「$$\nbegin{align} ～ end{align}$$\n」から「 '<'div align="center"'>' '['tex:\begin{align} ～\end{align}]' '<'/div'>' 」に変換する。
align環境の内の'[', ']'はエラーが出るため、エスケープする。

文中数式

文中$$は、$ ～ $は、'['tex: ～ ']'に変換する。
文中$$の内の'_','^', '*'はエラーが出るため、エスケープする。

まだまだ不具合、不足があると思うので逐次修正していきます。

gist21c28dae7f7f5e6cb2c3ede7041c4ce0

はてなブログの数式は全然うまくいかないけどとりあえずなんとかなりそうです。

参考サイト

atmarkit.itmedia.co.jp

messefor.hatenablog.com

リンク

2020-10-06

周期境界は配列の循環シフトを使えばソースコードがスッキリするかも

循環シフトとは

配列の最初と最後を循環してシフトするものです。 pythonのnumpy.rollで書くとこんな感じになります。

import numpy as np

x = np.arange(5)
x # array([0, 1, 2, 3, 4]) 

np.roll(x,2)
# array([4, 0, 1, 2, 3])

np.roll(x,-1)
# array([1, 2, 3, 4, 0])

いまいち何に使うのかわかりませんが(ビット演算とかでしょうか)、これを使えば周期境界の離散式がスッキリかけると思います。

ちなみにfortranだとcshiftです。

python

numpy.roll

https://note.nkmk.me/python-numpy-roll/

fortran

cshift

https://www.sci.hokudai.ac.jp/~inaz/doc/B/Fortran90/node19.html

循環シフトを用いた周期境界のコーディング

以下の単純な移流方程式を対象とする。

$\begin{align} \dfrac{\partial u }{\partial t } &= - c\dfrac{\partial u }{\partial x } \\ c &= \rm{const} \end{align}$

計算条件として、 $c=1$ および $u$ の初期値を次式で与える。

$\begin{align} u_0 = \mathrm{e}^{-(x-ct)^2} \end{align}$

また、 $dt=0.01,dx=0.1$ とする。

x = np.linspace(-5,5,101, endpoint=True )
u = np.exp(-x**2)

c = 1.0
dt = 0.01
dx = 0.1

周期境界で $u$ の時間発展を中央差分で解くコードは次のようになる。

ut = np.zeros_like(u)
for i in range(len(u)):
    if i == len(u)-1:
        ut[i] = u[i] - dt * c * 0.5 * (u[0] - u[i-1])/dx
    else:
        ut[i] = u[i] - dt * c * 0.5 * (u[i+1] - u[i-1])/dx

pythonの場合、-1のインデックスは配列の末尾を示すため、配列の最後の計算だけをif文で分岐させる必要がある。

または、スライシングを使って次のようになる。

ut = np.zeros_like(u)
ut[1:len(u)-1] = u[1:len(u)-1] - dt * c * 0.5 * (u[2:len(u)] - u[0:len(u)-2])/dx

i = 0 
ut[i] = u[i] - dt * c * 0.5 * (u[i+1] - u[i-1])/dx

i = len(u)-1
ut[i] = u[i] - dt * c * 0.5 * (u[0] - u[i-1])/dx

境界部分は上手くスライシングできないため、別に書く必要がある。

次に、このコードを循環シフトを使って書くと次のようになる。

ut = u - dt * c * 0.5*( np.roll(u,-1) - np.roll(u,1) )/dx

if文が不要で、ループがいらないためインデックスも不要である。

roll関数が若干わかりにくいので好みはわかれるがワンライナー好きとしてはこちらの方が良い。

計算例

先程の移流方程式を循環シフトを使って中央差分と後退差分(風上差分)で計算してみた。

10秒後の計算結果を示す。なお、周期境界のため真値は初期値と同じになる。

import numpy as np
import holoviews as hv

hv.extension('bokeh')

x = np.linspace(-5,5,101, endpoint=True )
uini = np.exp(-x**2)
c = 1.0
dt = 0.01
dx = 0.1

#中央差分
u = np.copy(uini)
for i in range(1000):
    ut = u - dt * c * 0.5*( np.roll(u,-1) - np.roll(u,1) )/dx
    u = np.copy(ut)
uc = np.copy(ut)

#風上差分
u = np.copy(uini)
for i in range(1000):
    ut = u - dt * c * (u - np.roll(u,1))/dx
    u = np.copy(ut)
ub = np.copy(ut)

g = hv.Curve((x,uini),label='initial & true') * hv.Curve((x,uc), label='central') * hv.Curve((x,ub), label='upwind')
g.options(width=500, legend_position='right')

f:id:SedimentHydraulics:20201006193508p:plain

まとめ

循環シフトはディープコピーになると思うので、速度面のデメリットはあるかもしれないが、ソースコードがすっきりするので、特に周期境界のときは役立つのではないかと思います。

gist

Jupyter Notebook Viewer

2020-10-01

河川合流の美しい風景

年に1回くらいこの記事が見たくなる。

karapaia.com

計算してみたいな。平面二次元で行けるのかな。三次元じゃないと無理かな。

2020-09-29

安定河道の自動設計をやりたい

前回の記事で椿先生の本¹を参照してときに、ペラペラのめくっていると安定河道の話があったのでそれを元に考えてみました。

河道の安定条件のモデル化

簡単のため、一次元の不等流場、断面形状は矩形断面を想定する。また、河床材料は単一粒径で掃流運動のみを考える。

河道が安定、つまり動的平衡状態であることは、各断面を通過する流砂量 $Q_b$ が縦断方向で一定となることを示している。

流砂量式の一般形を次のように定義する。

$\begin{align} q_{b*} &= \alpha \tau_*^{\beta} = \alpha \left( \dfrac{h i_e}{\rho_{s}d} \right)^{\beta} \end{align}$

エネルギー勾配はマニング則を適用すると次式となる。

$\begin{align} i_e &= \left( \dfrac{Qn}{Bh^{5/3}} \right)^2 \end{align}$

これらの式を用いると、断面積分した流砂量 $Q_b$ は次のように示される。

$\begin{align} Q_b &= B q_{b} = B \sqrt{\rho_s g d^3}q_{b*} \\ & = B^{1-2\beta} h^{-7\beta/3} \alpha \rho_s^{0.5-\beta} g^{0.5} d^{1.5-\beta} Q^{2\beta}n^{2\beta} \end{align}$

マニングの粗度係数が縦断方向に変化がないとすると、上式より各断面の流砂量が一定であることは $B^{1-2\beta} h^{-7\beta/3}$ が一定であることがわかる。

安定河道の河床高、川幅の計算方法

不等流計算の基礎式は次式である。

$\begin{align} \dfrac{d }{dx}\left( \dfrac{Q^2}{2gh^2B^2} + h + z_b \right) &= i_e \end{align}$

川幅が既知の場合は、安定河道を考慮した縦断的な河床高は、上式を変形した次式で示される。

$\begin{align} -\dfrac{d z_b}{dx} = \dfrac{d }{dx} \left( \dfrac{Q^2}{2gh^2B^2} + h \right) - i_e \end{align}$

水深は前項の安定条件より川幅のみの関数で示される。よって、本式を逐次計算することにより安定河道の河床高が計算できる。

次に、河床高が既知である場合、縦断的な川幅を示す式は次のように示される。

$\begin{align} \dfrac{d B}{dx} = \dfrac{\left(1 - \dfrac{Q^2}{gh^3B^2} \right)\dfrac{d h}{dx} + \dfrac{d z_b}{dx} - i_e }{\dfrac{Q^2}{gh^2B^3}} \end{align}$

河床高と同様に上式を解くことにより安定条件を満たす川幅を計算することが可能である。なお、水深が川幅の関数のため陰的に計算する必要がある。

式だけでまだコーディングはしてません。どうでしょうか。

最近、この椿先生の本や岸・黒木の論文など過去の文献を見ていると本当に気付きが多いです。

アーカイブを再発掘するだけで新しいアイデアが生まれそうな気がします。

椿:水理学演習下 pp.212-215
水理学演習下 POD版
- 作者:正夫, 荒木,東一郎, 椿
- 発売日: 2017/11/17
- メディア: 単行本
↩

従来の静水の問題

今後の展開

モチベーション

基礎方程式

との数式による表現

保存系のシミュレーション

TVD-MacCormack法

川幅の変化を考慮した結果

非保存系で解いた場合

結論

土砂濃度の鉛直分布の導出

基本的な考え方

Rouse分布

参考：式展開

Lane-Kalinske分布

2式の比較

濃度分布の比較

水深平均濃度の比較

水深平均濃度式の導出

比較

まとめ

参考文献

コードブロック

align環境

文中数式

参考サイト

循環シフトとは

循環シフトを用いた周期境界のコーディング

計算例

まとめ

gist

河道の安定条件のモデル化

安定河道の河床高、川幅の計算方法

$I_1$ と $I_2$ の数式による表現