.. sectionauthor:: 清水康行 .. _2dh_free: =============================== 2次元自由水面非定常流 =============================== 本章では, 自由水面を有する流れの2次元非定常流れの計算法について解説します. ここで扱う流れ場は水深方向の流速分布や鉛直方向の流速成分は考慮していないことから, 水深スケールが平面スケールに対して十分に小さい,すなわち,水深が相対的に浅い流れということで, 2次元浅水流とも呼ばれます. 2次元浅水流の基礎式は以下の連続式と, .. math:: :label: cnd \frac{\partial h}{\partial t}+\frac{\partial (uh)}{\partial x}+ \frac{\partial (vh)}{\partial y}=0 運動方程式から構成されます. .. math:: :label: mt_1 \frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+ v \frac{\partial u}{\partial y}= -g \frac{\partial H}{\partial x} -\frac{\tau_x}{\rho h}+ \frac{\partial}{\partial x} (\nu_t \frac{\partial u}{\partial x})+ \frac{\partial}{\partial y} (\nu_t \frac{\partial u}{\partial y}) .. math:: :label: mt_2 \frac{\partial v}{\partial t}+u \frac{\partial v}{\partial x}+ v \frac{\partial v}{\partial y}= -g \frac{\partial H}{\partial y} -\frac{\tau_y}{\rho h}+ \frac{\partial}{\partial x} (\nu_t \frac{\partial v}{\partial x})+ \frac{\partial}{\partial y} (\nu_t \frac{\partial v}{\partial y}) ただし, :math:`x, y` は互いに直交する平面座標軸, :math:`t` は時間, :math:`u, v` は :math:`x, y` 方向の水深平均流速, :math:`h` は水深, :math:`H` は水位, :math:`g` は重力加速度, :math:`\tau_x, \tau_y` は :math:`x, y` 方向の河床せん断力, :math:`\rho` は水の密度, :math:`\nu_t` は渦動粘性係数です. 抵抗則にマニング則を用いると, 河床せん断力は下記で表されます. .. math:: :label: frac \frac{\tau_x}{\rho h}= \frac{gn^2}{h^{4/3}} u\sqrt{u^2+v^2}, \ \ \ \frac{\tau_y}{\rho h}= \frac{gn^2}{h^{4/3}} v\sqrt{u^2+v^2} ただし, :math:`n` はマニングの粗度係数です. :eq:`frac` 式の導出についてはAppendixの :ref:`tau_exp` を参照ください. 渦動粘性係数は,ゼロ方程式モデルを採用し,下記で表されます. .. math:: :label: eddy \nu_t = \frac{\kappa}{6} u_\ast h ただし,:math:`u_\ast` は摩擦速度, .. math:: :label: fric_vel u_\ast=\sqrt{ghI_f} :math:`I_f` は摩擦勾配で次式で表されます. .. math:: :label: fric_slope I_f=\frac{n^2 (u^2+v^2)}{h^{4/3}} 基礎式と分離解法の適用 ========================= 自由水面を持つ2次元非定常流の計算では, :eq:`cnd` 式から水深 :math:`h` を, :math:`x` 方向の運動方程式 :eq:`mt_1` 式から :math:`x` 方向の流速 :math:`u` を, :math:`y` 方向の運動方程式 :eq:`mt_1` 式から :math:`y` 方向の流速 :math:`v` を算出します. なお, 水位 :math:`H` は河床高 :math:`\eta` が与えられた場合には :math:`H=\eta+h` で :math:`h` と一意的に対応しますので :math:`h` の代わりに :math:`H` 未知量と考えても同じことです. 1次元の時と同じように :eq:`mt_1` 式および :eq:`mt_2` 式に分離解法を適用します. ほとんど 同じ形なので :eq:`mt_1` 式の場合のみ説明します. :eq:`mt_1` 式に :math:`H=\eta+h` の関係を用いて 再記しておきます. .. math:: :label: mt_11 {{\partial u}\over{\partial t}}+u {{\partial u}\over{\partial x}}+v{{\partial u}\over{\partial y}} = -g \left( {{\partial h}\over{\partial x}}+{{\partial \eta}\over{\partial x}} \right) -{{g n_m^2 u\sqrt{u^2+v^2}}\over{h^{4/3}}} +\nu_t \left[ {{\partial^2 u}\over{\partial x^2}} + {{\partial^2 u }\over{\partial y^2}} \right] なお本式は拡散項も含むため, 移流項に関する部分, 圧力項と摩擦項に関する部分, および拡散項に関する部分 の3つに分離することとします. .. math:: :label: 2b_1 {{\partial u}\over{\partial t}} = -g \displaystyle{\left( {{\partial h}\over{\partial x}} + {{\partial \eta}\over{\partial x}} \right) }  -{{g n_m^2 u\sqrt{u^2+v^2}}\over{h^{4/3}}} .. math:: :label: 2b_2 {{\partial u}\over{\partial t}} = \nu_t \displaystyle{\left[ {{\partial^2 u }\over{\partial x^2}} +{{\partial^2 u}\over{\partial y^2}} \right] } .. math:: :label: 2b_3 {{\partial u}\over{\partial t}}+ u{{\partial u}\over{\partial x}} +v{{\partial u}\over{\partial y}} = 0 :math:`n` を時間ステップとして, :eq:`mt_11` 式で, :math:`u^n \rightarrow u^{n+1}` の計算をする 代わりに, 中間に :math:`\overline{u}` と :math:`\widetilde{u}` を設けて :eq:`2b_1` 式で :math:`u^n \rightarrow \overline{u}` を, :eq:`2b_2` 式で :math:`\overline{u} \rightarrow \widetilde{u}` を, :eq:`2b_3` 式で :math:`\widetilde{u} \rightarrow u^{n+1}` を計算することにします. なお, :eq:`2b_1` 式には :math:`h` が含まれているため :eq:`cnd` 式の 連続式も同時に満たす必要があります. これを :math:`h^n \rightarrow \overline{h}` と 表しておきます. これに対して :eq:`2b_2` 式, :eq:`2b_3` 式には :math:`h` が含まれていませんので, :math:`h^{n+1}` はこの :math:`\overline{h}` をそのまま用いることとします. 即ち, .. math:: :label: 2b_4 h^{n+1} = \overline{h} ということです. :math:`v` についてもまったく同様です. :math:`v` に関する運動方程式も 以下のように分離して計算します. .. math:: :label: 2b_5 {{\partial v}\over{\partial t}} = -g \displaystyle{\left( {{\partial h}\over{\partial y}} + {{\partial \eta}\over{\partial y}} \right) }  -{{g n_m^2 v\sqrt{u^2+v^2}}\over{h^{4/3}}} .. math:: :label: 2b_6 {{\partial v}\over{\partial t}} = \nu_t \displaystyle{\left[ {{\partial^2 v }\over{\partial x^2}} +{{\partial^2 v}\over{\partial y^2}} \right] } .. math:: :label: 2b_7 \cfrac{\partial v}{\partial t}+u\cfrac{\partial v}{\partial x} +v\cfrac{\partial v}{\partial y} = 0 以上を整理すると, .. list-table:: 2次元流れの分離解法 :header-rows: 1 * - Phase名 - 変数の更新方法 - - - 使用する式 * - Non-Advection Phase I - :math:`u^n \rightarrow \overline{u}` - :math:`v^n \rightarrow \overline{v}` - :math:`h^n \rightarrow \overline{h} \rightarrow h^{n+1}` - :eq:`2b_1` , :eq:`2b_5` , :eq:`cnd` , :eq:`2b_4` 式 * - Non-Advection Phase II - :math:`\overline{u} \rightarrow \widetilde{u}` - :math:`\overline{v} \rightarrow \widetilde{v}` - - :eq:`2b_2` , :eq:`2b_6` 式 * - Advection Phase - :math:`\widetilde{u} \rightarrow u^{n+1}` - :math:`\widetilde{v} \rightarrow v^{n+1}` - - :eq:`2b_3` , :eq:`2b_7` 式 ということになり, 以上でCIP法を適用する準備が整いました. Non Advection Phase I ========================= :eq:`2b_1` 式で :math:`u^n \rightarrow \overline{u}` を, :eq:`2b_5` 式で :math:`v^n \rightarrow \overline{n}` を計算するので :eq:`2b_1` 式および :eq:`2b_5` 式を以下のように表しておきます. .. math:: :label: 2d_bb1 {{\overline{u}-u^n}\over{\Delta t}} = -g \displaystyle{ \left( {{\partial \overline{h}}\over{\partial x}} + {{\partial \eta}\over{\partial x}} \right) } -{{g n_m^2 \overline{u} \sqrt{\overline{u}^2 +\overline{v}^2}} \over{\overline{h}^{4/3}}} .. math:: :label: 2d_bb2 {{\overline{v}-v^n}\over{\Delta t}} = -g \displaystyle{ \left( {{\partial \overline{h}}\over{\partial y}} + {{\partial \eta}\over{\partial y}} \right) } -{{g n_m^2 \overline{v} \sqrt{\overline{u}^2 +\overline{v}^2}} \over{\overline{h}^{4/3}}} これより, .. math:: :label: 2d_bb3 \overline{u}= u^n - g \Delta t {{\partial \overline{h}}\over{\partial x}} - g \Delta t {{\partial \eta}\over{\partial x}} - g \Delta t {{n_m^2 \overline{u}\sqrt{\overline{u}^2 +\overline{v}^2}} \over{\overline{h}^{4/3}}} .. math:: :label: 2d_bb4 \overline{v}= v^n - g \Delta t {{\partial \overline{h}}\over{\partial y}} - g \Delta t {{\partial \eta}\over{\partial y}} - g \Delta t {{n_m^2 \overline{v}\sqrt{\overline{u}^2 +\overline{v}^2}} \over{\overline{h}^{4/3}}} 一方, :eq:`cnd` 式で :math:`h^n \rightarrow \overline{h}` を求めるので, .. math:: :label: 2d_bb5 {{\overline{h}-h^n}\over{\Delta t}} + {{\partial}\over{\partial x}} (\overline{u} \overline{h}) + {{\partial}\over{\partial y}} (\overline{v} \overline{h})=0 を満たす必要があります. 従って, :eq:`2d_bb3` :math:`\sim` :eq:`2d_bb5` 式を連立して得られる :math:`\overline{h}` , :math:`\overline{u}` および :math:`\overline{v}` を求める必要があります. CIP法で流体の計算を行う場合は一般に :math:`u` , :math:`v` と :math:`h` の計算点を交互に配置するスタッカード格子が用いられます. ここでは, :numref:`plane_uvij` に示すような計算点の配置を採用することにします. .. _plane_uvij: .. figure:: images/05/2d_uv_location.png :width: 100% : 計算格子の配置法 :numref:`plane_uvij` の計算点の配置を考慮して :eq:`2d_bb3` 式を差分表示すると, .. math:: :label: 2d_cc1 \overline{u}(i,j) = & u^n(i,j)-g\Delta t {{\overline{h}(i+1,j) - \overline{h}(i,j)}\over{\Delta x}} -g\Delta t {{\eta(i+1,j) - \eta(i,j)}\over{\Delta x}} \\ & - g \Delta t {{n_m^2 \overline{u}(i,j) \sqrt{\overline{u}(i,j)^2+\overline{v_{up}}^2} }\over{\overline{h_{up}}^{4/3}}} ここで, 下付き添え字 :math:`_{up}` は :math:`u` の計算点における値であることを表し, :numref:`plane_uvij` の計算点の配置によれば, .. math:: :label: 2d_cc2 \overline{v_{up}}={{\overline{v}(i,j)+\overline{v}(i,j-1) +\overline{v}(i+1,j)+\overline{v}(i+1,j-1)}\over 4} .. math:: :label: 2d_cc3 \overline{h_{up}}={{\overline{h}(i,j)+\overline{h}(i+1,j)}\over 2} となります. 同様に :eq:`2d_bb4` 式を差分表示すると, .. math:: :label: 2d_cc4 \overline{v}(i,j) = & v^n(i,j)-g\Delta t {{\overline{h}(i,j+1) - \overline{h}(i,j)}\over{\Delta y}} -g\Delta t {{\eta(i,j+1) - \eta(i,j)}\over{\Delta y}} \\ & - g \Delta t {{n_m^2 \overline{v}(i,j) \sqrt{\overline{u_{vp}}^2+\overline{v}(i,j)^2} }\over{\overline{h_{vp}}^{4/3}}} となります. 下付き添え字 :math:`_{vp}` は :math:`u` の計算点における値であることを表し, .. math:: :label: 2d_cc5 \overline{u_{vp}}={{\overline{u}(i,j)+\overline{u}(i-1,j) +\overline{u}(i,j+1)+\overline{u}(i-1,j+1)}\over 4} .. math:: :label: 2d_cc6 \overline{h_{vp}}={{\overline{h}(i,j)+\overline{h}(i,j+1)}\over 2} で表されます. :eq:`2d_cc1` 式および :eq:`2d_cc4` 式を :eq:`2d_bb5` 式の 左辺第2項および第3項に代入するために, .. math:: :label: 2d_cc7 \overline{q_u}=\overline{u} \overline{h}, \; \; \; \overline{q_v}=\overline{v} \overline{h} を定義し( :math:`q_u`nおよび :math:`q_v` の計算点はそれぞれ :math:`u` および :math:`v` と同じとする) 計算格子点の配置を考慮して, .. math:: :label: 2d_cc8 \overline{q_x}(i,j)= \overline{u}(i,j) {{\overline{h}(i+1,j)+\overline{h}(i,j)}\over{2}} .. math:: :label: 2d_cc9 \overline{q_y}(i,j)= \overline{v}(i,j) {{\overline{h}(i,j+1)+\overline{h}(i,j)}\over{2}} および, :math:`q_x` , :math:`q_y` を用いた :eq:`2d_bb5` 式の差分式 .. math:: :label: 2d_cc10 \overline{h}(i) =h^n(i,j) + \left[ {{\overline{q_x}(i,j)-\overline{q_x}(i-1,j)}\over{\Delta x}} +{{\overline{q_y}(i,j)-\overline{q_y}(i,j-1)}\over{\Delta y}} \right]\Delta t より :math:`\overline{h}` を求めることが出来ます. ただし, :eq:`2d_cc10` 式 には右辺の :math:`\overline{q_x}` , :math:`\overline{q_y}` にも :math:`\overline{h}` が陰的な形で含まれるため, 繰り返し計算が 必要になります. この手順を :numref:`2d_flow_chart` に示します. .. _2d_flow_chart: .. figure:: images/04/2d_flow_chart.png :width: 70% : Non-Advection Phase I における :math:`\overline{h}(i,j)` , :math:`\overline{u}(i,j)` および :math:`\overline{v}(i,j)` の計算手順 この繰り返し計算を通して得られる :math:`\overline{h}(i)` , :math:`\overline{u}(i,j)` および :math:`\widetilde{v}(i,j)` が **Non-Advection Phase I** の計算結果です. **Non-Advection Phase I** における :math:`\overline{u}(i,j)` と :math:`\overline{v}(i,j)` が求まった段階で, **Non-Advection Phase I** における :math:`\displaystyle{{\partial u}\over{\partial x}}` , :math:`\displaystyle{{\partial u}\over{\partial y}}` , :math:`\displaystyle{{\partial v}\over{\partial x}}` および :math:`\displaystyle{{\partial v}\over{\partial y}}` すなわち, :math:`\overline{\displaystyle{{\partial u}\over{\partial x}}}` , :math:`\overline{\displaystyle{{\partial u}\over{\partial y}}}` , :math:`\overline{\displaystyle{{\partial v}\over{\partial x}}}` および :math:`\overline{\displaystyle{{\partial v}\over{\partial y}}}` を求めておく必要があります. これは1次元の時にやったのと全く同様に以下の式で 計算できます. .. math:: :label: tp1 \overline{{\partial u}\over{\partial x}}(i,j)= {{\partial u}\over{\partial x}}^n(i,j) +{{1}\over{2\Delta t}} \left[ \overline{u}(i+1,j)-u^n(i+1,j)-\overline{u}(i-1,j)+u^n(i-1,j) \right] .. math:: :label: tp2 \overline{{\partial u}\over{\partial y}}(i,j)= {{\partial u}\over{\partial y}}^n(i,j) +{{1}\over{2\Delta t}} \left[ \overline{u}(i,j+1)-u^n(i,j+1)-\overline{u}(i,j-1)+u^n(i,j-1) \right] .. math:: :label: tp3 \overline{{\partial v}\over{\partial x}}(i,j)= {{\partial v}\over{\partial x}}^n(i,j) +{{1}\over{2\Delta t}} \left[ \overline{v}(i+1,j)-v^n(i+1,j)-\overline{v}(i-1,j)+v^n(i-1,j) \right] .. math:: :label: tp4 \overline{{\partial v}\over{\partial y}}(i,j)= {{\partial v}\over{\partial y}}^n(i,j) +{{1}\over{2\Delta t}} \left[ \overline{v}(i,j+1)-v^n(i,j+1)-\overline{v}(i,j-1)+v^n(i,j-1) \right] Non Advection Phase II =============================== ここでは拡散項の計算を行います. 拡散項は :eq:`2b_2` 式および :eq:`2b_6` 式を用いて, :math:`\overline{u} \rightarrow \widetilde{u}` および :math:`\overline{v} \rightarrow \widetilde{v}` の計算を行います. 拡散項は単純に中央差分が適用できますので, :eq:`2b_2` 式および :eq:`2b_6` 式 の差分表示は以下のようになります. .. math:: \widetilde{u}(i,j) = \overline{u}(i,j) + \Delta t \nu_t \left[ {{\overline{u}(i+1,j)-2\overline{u}(i,j)+\overline{u}(i-1,j)}\over{\Delta x^2}} \right. .. math:: :label: df2_1 \left. +{{\overline{u}(i,j+1)-2\overline{u}(i,j)+\overline{u}(i,j-1)}\over{\Delta y^2}}\right] .. math:: \widetilde{v}(i,j) = \overline{v}(i,j) + \Delta t \nu_t \left[ {{\overline{v}(i+1,j)-2\overline{v}(i,j)+\overline{v}(i-1,j)}\over{\Delta x^2}} \right. .. math:: :label: df2_2 \left. +{{\overline{v}(i,j+1)-2\overline{v}(i,j)+\overline{v}(i,j-1)}\over{\Delta y^2}}\right] Advection Phase ========================== Advection Phaseで計算する式は :eq:`2b_3` 式と :eq:`2b_7` 式です. これは, 前章の2次元移流方程式( :ref:`2d_advection_eq` )そのものです. :eq:`2d_advection` 式で, :math:`f` を :math:`u` としたのが :eq:`2b_3` 式で, :math:`f` を :math:`v` としたのが :eq:`2b_7` 式ですので前章の内容をそのまま適用できます. ただし, 今回は :math:`u` と :math:`v` の計算点の配置方法が異なりますのでこの点だけ注意する 必要があります. たとえば, :eq:`2d_3` 式にCIP法を適用する場合, 左辺第3項に現れる :math:`v` は :math:`u` の計算点のものを使う必要があり, :eq:`2d_cc2` 式で行ったような周辺の値からの換算が必要になります. 同様に :eq:`2b_7` 式の左辺第2項に現れる :math:`u` は :math:`v` の計算点の :math:`u` を使用する必要があり, 同様に :eq:`2d_cc5` 式のような変換が必要になります. これらを考慮して, **Advection Phase** の計算は以下のように式で行われます. Advection Phaseの *u* -------------------------------- .. math:: u(i,j)^{n+1}= \left[ (a_1 X + c_1 Y + e_1) X + g_1 Y + \overline{{\partial u}\over{\partial x}}(i,j) \right] X + .. math:: :label: 2d_adv_1 \left[ (b_1 Y + d_1 X + f_1) Y + \overline{{\partial u}\over{\partial y}}(i,j) \right] Y + \widetilde{u}(i,j) ただし, .. math:: :label: 2d_adv_2 X=-\widetilde{u}(i,j)\Delta t, \; \; \; Y=-\widetilde{v_{up}}(i,j)\Delta t .. math:: :label: 2d_adv_3 \widetilde{v_{up}}={{\widetilde{v}(i,j)+\widetilde{v}(i,j-1) +\widetilde{v}(i+1,j)+\widetilde{v}(i+1,j-1)}\over 4} なお :eq:`2d_adv_1` 式中の :math:`a_1 \sim g_1` は前章の :eq:`coefs` 式で :math:`f` を :math:`\widetilde{u}` に,:math:`f_x` および :math:`f_y` を :math:`u_x` および :math:`u_y` に 置き換えたものを使用します.ただし, .. math:: :label: 2d_adv_4 u_x=\cfrac{\partial u}{\partial x}, \; \; \; u_y=\cfrac{\partial u}{\partial x} です. Advection Phaseの *v* ------------------------------ .. math:: v(i,j)^{n+1}= \left[ (a_1 X + c_1 Y + e_1) X + g_1 Y + \overline{{\partial v}\over{\partial x}}(i,j) \right] X + .. math:: :label: 2d_adv_5 \left[ (b_1 Y + d_1 X + f_1) Y + \overline{{\partial v}\over{\partial y}}(i,j) \right] Y + \widetilde{v}(i,j) ただし, .. math:: :label: ad_adv_6 X=-\widetilde{u_{vp}}(i,j)\Delta t, \; \; \; Y=-\widetilde{v}(i,j)\Delta tp3 .. math:: :label: 2d_2dv_7 \widetilde{u_{vp}}={{\widetilde{u}(i,j)+\widetilde{u}(i-1,j) +\widetilde{u}(i,j+1)+\widetilde{u}(i-1,j+1)}\over 4} なお :eq:`2d_adv_5` 式中の :math:`a_1 \sim g_1` は :eq:`coefs` 式で :math:`f` を :math:`\widetilde{v}` に, :math:`f_x` および :math:`f_y` を :math:`v_x` および :math:`v_y` に 置き換えたものを使用します. ただし, .. math:: :label: 2d_adv_8 v_x=\cfrac{\partial v}{\partial x}, \; \; \; v_y=\cfrac{\partial v}{\partial x} です. Advection Phaseにおける微分量の更新 ------------------------------------------ 前章の :eq:`bibun2_2` 式および :eq:`bibun2_3` 式を用いた微分量の更新が必要となります. 確認のためにこれらの式を再記しておきます. .. math:: :label: bibun2_21 {{\partial f_x}\over{\partial t}} + u{{\partial f_x}\over{\partial x}} + v{{\partial f_x}\over{\partial y}} = -\left( f_x {{\partial u}\over{\partial x}} + f_y {{\partial v}\over{\partial x}} \right) .. math:: :label: bibun2_31 {{\partial f_y}\over{\partial t}} + u{{\partial f_y}\over{\partial x}} + v{{\partial f_y}\over{\partial y}} = -\left( f_x {{\partial u}\over{\partial y}} + f_y {{\partial v}\over{\partial y}} \right) まず, :math:`bibun2_21` 式で :math:`f` が :math:`u` の場合を考えます. 左辺には分離解法によるCIP法を適用できますので, 2段階に分けて, CIPによる移流部分とそれ以外に分けて計算します. たとえば, :math:`u` の :math:`x` 微分に関する更新を, .. math:: :label: h_advec_1 \displaystyle{\overline{{\partial u}\over{\partial x}} \longrightarrow \widehat{{\partial u}\over{\partial x}} \longrightarrow {{\partial u}\over{\partial x}}^{n+1}} と表示することにしますと, 移流部分は, .. math:: :label: h_advec_2 \widehat{{\partial u}\over{\partial x}}(i,j)= \left[3 a_1 X + 2(c_1 Y +e_1) \right] X + (d_1 Y +g_1) Y +\overline{{\partial u}\over{\partial x}}(i,j) それ以外の部分は .. math:: {{\partial u}\over{\partial x}}^{n+1}(i,j)= \widehat{{\partial u}\over{\partial x}}(i,j) .. math:: :label: h_advec_3 -\left[ \widehat{{\partial u}\over{\partial x}}(i,j) {{\widetilde{u}(i+1,j)-\widetilde{u}(i-1,j)}\over{2 \Delta x}} +\widehat{{\partial u}\over{\partial y}}(i,j) {{\widetilde{v_{up}}(i+1,j)-\widetilde{v_{up}}(i-1,j)}\over{2 \Delta x}} \right] \Delta t によって計算することが出来ます( :math:`a_1, c_1, e_1, d_1, g_1` は前節を参照願います). 同様に, :math:`u` の :math:`y` に関する微分量の更新は, .. math:: :label: h_advec_4 \widehat{{\partial u}\over{\partial y}}(i,j)= \left[3 b_1 Y + 2(d_1 X +f_1) \right] Y + (c_1 X +g_1) X +\overline{{\partial u}\over{\partial y}}(i,j) .. math:: {{\partial u}\over{\partial y}}^{n+1}(i,j)= \widehat{{\partial u}\over{\partial y}}(i,j) .. math:: :label: h_advec_5 -\left[ \widehat{{\partial u}\over{\partial x}}(i,j) {{\widetilde{u}(i,j+1)-\widetilde{u}(i,j-1)}\over{2 \Delta y}} +\widehat{{\partial u}\over{\partial y}}(i,j) {{\widetilde{v_{up}}(i,j+1)-\widetilde{v_{up}}(i,j-1)}\over{2 \Delta y}} \right] \Delta t によって計算することが出来ます. なお, :eq:`h_advec_3` 式および :eq:`h_advec_5` 式右辺第2項と第3項の :math:`\displaystyle{\widehat{{\partial u}\over{\partial x}}}` および :math:`\displaystyle{\widehat{{\partial u}\over{\partial y}}}` は, それぞれ直前の :eq:`h_advec_2` 式や :eq:`h_advec_4` 式で計算されるものを用いるよりは, 別途, :math:`\widetilde{u}` を用いて計算し直したほうが計算の安定性が高まる ことが知られております. この場合, :eq:`h_advec_3` 式と :eq:`h_advec_5` 式の代わりに以下の式を用いるということになります. .. math:: {{\partial u}\over{\partial x}}^{n+1}(i,j)= \widehat{{\partial u}\over{\partial x}}(i,j) -\left\{ {{[\widetilde{u}(i+1,j)-\widetilde{u}(i-1,j)]}\over{2 \Delta x}} {{[\widetilde{u}(i+1,j)-\widetilde{u}(i-1,j)]}\over{2 \Delta x}} \right. .. math:: :label: h_advec_6 \left. + {{[\widetilde{u}(i,j+1)-\widetilde{u}(i,j-1)]}\over{2 \Delta y}} {{[\widetilde{v_{up}}(i+1,j)-\widetilde{v_{up}}(i-1,j)]}\over{2 \Delta x}} \right\} \Delta t .. math:: {{\partial u}\over{\partial y}}^{n+1}(i,j)= \widehat{{\partial u}\over{\partial y}}(i,j) -\left\{ {{[\widetilde{u}(i+1,j)-\widetilde{u}(i-1,j)]}\over{2 \Delta x}} {{[\widetilde{u}(i,j+1)-\widetilde{u}(i,j-1)]}\over{2 \Delta y}} \right. .. math:: :label: h_advec_7 \left. + {{[\widetilde{u}(i,j+1)-\widetilde{u}(i,j-1)]}\over{2 \Delta y}} {{[\widetilde{v_{up}}(i,j+1)-\widetilde{v_{up}}(i,j-1)]}\over{2 \Delta y}} \right\} \Delta t 演習問題(7) ^^^^^^^^^^^^^^ .. list-table:: * - Advection Phaseにおける :math:`\displaystyle{\cfrac{\partial v}{\partial x}}` , :math:`\displaystyle{\cfrac{\partial v}{\partial y}}` の計算法を説明してください Pythonコードと計算例 ================================ Pythonコード -------------------- 以上で説明した2次元非定常自由水面流れの計算用Pythonコードをこちらに示します. https://github.com/YasuShimizu/2DH_Python 計算条件に関するパラメーターはこのレポジトリ内のconfig.ymlに記述します. なお,このコードでは本章で説明した以外に水路側壁(河岸)の摩擦を考慮出来るようにして あります.水路側壁(河岸)の摩擦の計算方法に関する解説を,Appendixの :ref:`bank_friction` に 記載します. また,このコードでは水路中の障害物を考慮出来るようになっています.構造物の個数および位置は 同じレポジトリ中のobst.datというテキストファイルに記述するようになっており, 以下に示す計算例では,障害物お個数が1個で :math:`x` 方向に :math:`i=` 10~14番目のセル, :math:`y` 方向に :math:`j=` 6~14番目のセルが障害物であるという設定になってます. これによって,直線水路の上流側中央部に長方形の障害物があって,これによって流れが分断され 後方に間欠的なカルマン渦が発生する様子が計算されています. .. _2dh_anim: .. figure:: images/05/2dh_animation.gif :width: 100% : 計算結果のアニメーション(流速ベクトルと流速による渦度) NaysMiniについて ===================== `NaysMini `_ は,本章で解説した 直交座標2次元モデルをiRICのGUIと連動させて使用するためのソルバで, ソースコードもすべて公開 されております. `NaysMiniのマニュアル `_ も公開されておりますので, 併せてご覧ください.