.. sectionauthor:: 清水康行 ======================= 1次元河床変動計算 ======================= 掃流砂のみの場合 --------------------- 1次元流砂連続式(掃流砂) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 長方形断面水路における1次元流砂連続式を示します． .. math:: :label: 1d_conb \cfrac{\partial z_b}{\partial t}+\cfrac{1}{(1-\lambda)B} \cfrac{\partial(B q_b)}{\partial x}=0 ただし, :math:`B` は水路幅, :math:`\lambda` は河床材料の空隙率です． :eq:`1d_conb`　式の概念を :numref:`scon1d` に示します．図のように :math:`\Delta x` 離れた2断面で流砂の 収支を考え，この差が :math:`\Delta t` 時間継続した結果, 河床高が :math:`\Delta z_b` 変化するとしたものです． 空隙率の考え方は, 流砂量の収支によって生ずる河床変動量の体積は，空隙 :math:`\lambda` も含んで堆積(浸食)する ことを考慮したものです． .. _scon1d: .. figure:: images/07/scon1d.png :width: 80% : 1次元河床変動概念図(掃流砂) .. _stkd_qb: .. figure:: images/07/stkd_qb.png :width: 95% : スタカード格子(1次元河床変動計算) 掃流砂量に関する計算格子を :numref:`stkd_qb` のように配置した場合の :eq:`1d_conb` 式の差分表示は次式のように なります．なお，ここでの :math:`q_b` の計算点は :numref:`1d_stk_b` の :math:`u` の計算点, :math:`z_b` の計算点は :numref:`1d_stk_b` の :math:`h` の 計算点と同じです. .. math:: :label: 1d_qb_fd z_b(i)^{n+1}=z_b(i)^n-\cfrac{\Delta t}{(1-\lambda)B(i)\Delta x} \left\{q_b(i)B_u(i)-q_b(i-1)B_u(i-1)\right\} なお, :math:`B_u` は :math:`q_b` の計算点における :math:`B` の値です. また, q_b(i)は各計算点において, :eq:`MPM` 式や :eq:`Ashida` 式で求め, これら流砂量式で使用される :math:`\tau_\ast` は :eq:`tausta` 式で求めます． ここで, :eq:`tausta` 式で使用される :math:`I_e` はマニング式を用いる場合は, .. math:: :label: Ie_eq I_e= \cfrac{{n_m}^2 u^2}{h^{4/3}} で求められます． 計算手順 ^^^^^^^^^^^^^^^^^^^ 計算手順は, :numref:`flowchart5` の最後に :eq:`1d_qb_fd` 式の河床変動計算を 追加した以下のフローとなります. .. _flow_chart6: .. figure:: images/07/flow_chart6.png :width: 80% : 1次元河床変動計算フローチャート Pythonコード ^^^^^^^^^^^^^^^^^^^^ 上記計算手順に沿ったPythonコードはこちらです． https://github.com/YasuShimizu/1D_Flow_and_Bed_Deformation 計算の本体部分はこのレポジトリ内のmain.pyで，ここからこのレポジトリ内に記載された各々の 関数.py が呼び出されるようになっています． 計算条件は同じレポジトリ内のymlファイルで指定するようになっており下記４つの計算条件ファイルが含まれています． 1. 縦断方向中央部のマウンドがある河川の河床変動 : config_mound.yml 2. 縦断方向中央部に狭窄部がある河川の河床変動 : config_constriction.yml 3. 途中で勾配変化点を含む河川の河床変動: config_slope_change.yml 4. 上流から流砂の供給が無い河川（ダム下流）の河床変動: downstream_dam.yml 計算の実行は，例えば上記 1. の計算を実施したい場合はMinicondaのプロンプトで， > python main.py config_mound と，引数でファイル名の config_xxxxx.yml でxxxxの部分を指定してやると， その条件が記述されたymlファイルが呼び出されて実行されるようになっております． 計算結果 ^^^^^^^^^^^^^^^^^^^ 上記の4通りの条件による1次元河床変動計算結果を以下に示します． .. _exam_1: .. figure:: images/07/mound.gif :width: 95% : 縦断方向中央部のマウンドがある河川の河床変動計算結果 .. _exam_2: .. figure:: images/07/constriction.gif :width: 95% : 縦断方向中央部に狭窄部がある河川の河床変動計算結果 .. _exam_3: .. figure:: images/07/slope_change.gif :width: 95% : 途中で勾配変化点を含む河川の河床変動計算結果 .. _exam_4: .. figure:: images/07/downstream_dam.gif :width: 95% : 上流から流砂の供給が無い河川（ダム下流）の河床変動計算結果 浮遊砂による河床変動 ------------------------------- 1次元浮遊砂濃度の連続式 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ :numref:`sus_cont` に示すように，長方形断面流路中に間隔が :math:`\Delta x` のI, II断面を考え， 浮遊砂の連続条件を考えます. .. _sus_cont: .. figure:: images/07/sus_cont.png :width: 80% : 1次元浮遊砂の連続条件概念図 上記 :math:`\Delta x` 区間の見ずの体積を :math:`V`, 平均濃度を :math:``, 平均水深を :math:`h`, 区間内の平均水路幅を :math:`B` ,通過する浮遊砂量を :math:`Q_s` とすると， .. math:: :label: 1d_con_sus \cfrac{\partial(V)}{\partial t} \Delta t + \cfrac{\partial{Q_s}}{\partial x} \Delta x \Delta t =(q_{su}-w_f c_b) \Delta x \Delta t B ただし, :math:`q_su` は単位面積単位時間当たりの河床からの浮遊砂浮上量, :math:`w_f` は沈降速度, :math:`c_b` は河床における浮遊砂濃度です． また, .. math:: :label: vqs V=A\Delta x = BH \Delta x, \; \; \; Q_s=Q=Bhu なので， .. math:: :label: vqs_div \cfrac{\partial(V)}{\partial t}=B \Delta x \cfrac{\partial}{\partial t} (h) \\ \cfrac{\partial Q_S}{\partial x}=\cfrac{\partial}{\partial x}(Bhu) よって, :eq:`1d_con_sus` 式は .. math:: :label: 1d_cons_sus1 \cfrac{\partial(h)}{\partial t} B \cancel{\Delta x \Delta t} + \cfrac{\partial(Bhu)}{\partial x} \cancel{\Delta x \Delta t} = (q_{su}-w_f c_b) B \cancel{\Delta x \Delta t} さらに， .. math:: :label: 1d_cons_sus2 \cfrac{\partial h}{\partial t}+ h\cfrac{\partial}{\partial t}+ \cfrac{1}{B} \left\{ \cfrac{\partial(Bhu)}{\partial x}+ Bhu\cfrac{\partial}{\partial x} \right\} =q_{su}-w_f c_b .. math:: :label: 1d_cons_sus3 \cancelto{\mbox{連続式より}=0}{\left\{ \cfrac{\partial h}{\partial t} + \cfrac{1}{B} \cfrac{\partial(Buh)}{\partial x} \right\}} + h\cfrac{\partial}{\partial t} + hu \cfrac{\partial}{\partial x} =q_{su}-w_f c_b したがって, .. math:: :label: 1d_cons_sus4 \cfrac{\partial}{\partial t} + u \cfrac{\partial}{\partial x} =\cfrac{1}{h}(q_{su}-w_f c_b) なお, :eq:`1d_cons_sus4` 式において :math:`` と :math:`c_b` に :eq:`ave_c` 式の関係を用いると, .. math:: :label: 1d_cons_sus5 \cfrac{\partial}{\partial t} + u \cfrac{\partial}{\partial x} =\cfrac{1}{h}(q_{su}-\alpha w_f ) となります．ただし, .. math:: :label: 1d_cons_sus6 \alpha=\cfrac{\beta}{1-\exp(-\beta)}, \; \; \beta=\cfrac{w_f h}{\epsilon} です． さらに，浮遊砂濃度の拡散を考慮する場合, :eq:`1d_cons_sus5` の右辺に 拡散項を付加し，次式となります． .. math:: :label: 1d_cons_sus61 \cfrac{\partial}{\partial t} + u \cfrac{\partial}{\partial x} = D\cfrac{\partial ^2}{\partial x^2}+ \cfrac{1}{h}(q_{su}-\alpha w_f ) この結果， :eq:`1d_cons_sus61` 式が水深平均濃度の連続式で，これを差分法で計算することになります． :eq:`1d_cons_sus61` 式は, :eq:`sr1` 式で示した, Source項を含む移流方程式と同じ形で, :eq:`sr1` 式における :math:`f` を :math:`` に, :math:`G` を :math:`D\cfrac{\partial ^2}{\partial x^2}+\cfrac{1}{h}(q_{su}-\alpha w_f )` に置き換えたものとなり, 本テキストで用いてきたCIP法を用いることにより精度の高い計算が可能です． なお, :math:`` の計算では計算点を下図のように :math:`u` の計算点の 間に取りますので，新たに :math:`` の計算点における :math:`u` として :math:`u_x` :eq:`u_x` 式のように定義し，これを :eq:`1d_cons_sus5` 式の移流速度 :math:`u` として用いて計算に用います． .. _stkd_uc: .. figure:: images/07/stkd_uc.png :width: 100% : 濃度計算点の配置 .. math:: :label: u_x u_x(i)=\cfrac{u(i)+u(i-1)}{2} 1次元浮遊砂による河床の連続式 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 水深平均濃度 :math:`` が求まると，河床における濃度 :math:`c_b` も 同時に求まっていることになりますので，河床変動の計算は簡単です． .. _scon1d_s: .. figure:: images/07/scon1d_s.png :width: 75% : 浮遊砂による河床変動模式図 .. math:: :label: sus_zb \cfrac{\partial z_b}{\partial t}= -\cfrac{1}{1-\lambda}(q_{su}-\alpha w_f ) Pythonコード ^^^^^^^^^^^^^^^^^^^^ 上記計算手順に沿ったPythonコードはこちらです． https://github.com/YasuShimizu/Bed-deformation-with-suspended-load 計算結果 ^^^^^^^^^^^ 上記計算コードを用いて行った貯水地への浮遊砂の堆砂の計算例を :numref:`res` に示します． 計算条件は, 河床勾配 :math:`i=1/700` , 延長 :math:`L` =10km, 川幅 :math:`B` =100(m), Manningの粗度係数 :math:`n` =0.02, 河床材料の粒径 :math:`d` =0.1mm, の長方形断面に流量 :math:`Q` = 500m :math:`^3` /s が 流下し, 下流端の水位が :math:`H` =12mのダム貯水池を想定しています． なお，この計算例は `現場のための水理学(4) -浮遊砂と河床変動- `_ で例題として掲載されている計算例と同じ条件で，計算結果も同じものが得られてます． .. _res: .. figure:: images/07/res.gif :width: 95% : 貯水地への浮遊砂の堆砂の計算