実河川スケールの水深平均開水路乱流Depth Averaged Turbulenceの計算：まとめ - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

前回の記事（実河川スケールの水深平均開水路乱流Depth Averaged Turbulenceの計算：結果だけ - 趣味で計算流砂水理）の続きです。

モチベーション
計算モデル
結果と考察
式の導出
ニュートン流体の粘性応力
- 参考資料
ブシネスクの渦粘性近似
- 補足：レイノルズ応力の付加項
- 参考資料
河川流の乱れエネルギー
- 参考資料
Github
図化等の参考サイト

モチベーション

最近新しい河川流の平面二次元のプログラムを書いていて乱流モデルをどうしようかなと悩んでます。

いつもどおり、実河川のモデルなので、計測可能な地形データの解像度（ALB測量を使って水深平均二次元流計算をすると凄い結果がでた - 趣味で計算流砂水理参照）を考えるとメッシュサイズは1m程度が最小と考えています。

このサイズだと、乱流モデルではラフなモデルになるので0方程式モデルで十分と考えています。

河川流の平面二次元でよく使う0方程式モデルは次の２式です。（x方向の運動方程式）

$\begin{align} &\mathrm{case1} : \frac{\partial q_x}{\partial t}+\frac{\partial u q_x}{\partial x}+\frac{\partial v q_x}{\partial y}+gh\frac{\partial H}{\partial x}+\frac{\tau_{0x}}{\rho} - \frac{\partial }{\partial x} \left(h \nu_t \dfrac{\partial u}{\partial x}\right) - \frac{\partial }{\partial y} \left(h \nu_t \dfrac{\partial u}{\partial y}\right) = 0 \\ &\mathrm{case2} : \frac{\partial q_x}{\partial t}+\frac{\partial u q_x}{\partial x}+\frac{\partial v q_x}{\partial y}+gh\frac{\partial H}{\partial x}+\frac{\tau_{0x}}{\rho} - \frac{\partial }{\partial x} \left(2h \nu_t \dfrac{\partial u}{\partial x}- \dfrac{2}{3}hk \right ) - \frac{\partial }{\partial y}\left( h \nu_t \left(\dfrac{\partial u}{\partial y} + \dfrac{\partial v}{\partial x}\right) \right) = 0 \end{align}$

導出等は後述します。

そのため、２式の比較をしてみます。

計算モデル

水深平均平面二次元河川流モデル + ０方程式乱流モデル

離散化とかいつもどおりなのでソースを見てください。

計算水路の条件

解析範囲：流下方向1800m、横断方向90m
河床地形：流下方向1/1000、横断方向勾配なし
粗度係数：横断方向左岸から0～30m：0.10、30～90m：0.025

図中左から右が流下方向で粗度係数の分布はこのようになります。

f:id:SedimentHydraulics:20211001111620p:plain

計算条件

境界条件：流下方向：周期境界、横断方向：壁境界-スリップ条件（かなりいい加減）
初期条件：水深平均6m、流速：流下方向：等流計算で設定し0~0.1％の乱数を乗算、横断方向：0
計算条件：dx, dy ： 1m、dt：0.04（クーラン数0.2くらい）

プログラム

かなり汚いですが、とりあえず貼っておきます。

https://nbviewer.jupyter.org/github/computational-sediment-hyd/DepthAveragedTurbulence/blob/main/numericalModel.ipynb

結果と考察

計算結果として、計算開始5時間後の各caseの流速、渦度の空間分布を示します。渦度は $\dfrac{\partial v}{\partial x} - \dfrac{\partial u}{\partial y}$ のため、プラスが反時計回り、マイナスが時計回りになります。

拡大図はこちら

両ケースともに、粗度係数の境界付近に安定した渦が形成、維持された状態となり、その大きさ、間隔、等はほぼ同程度になります。なお、確認のため、粗度を一定とした場合の計算も行いましたが、10~20分くらいで初期の擾乱は消えて渦がない状態になりました。

次に計算開始から5時間後までの変化を10分間隔で図化したものをgifで示します（クリックで大きくなります）。

f:id:SedimentHydraulics:20210930232504g:plain

両ケースともにほぼ同じような結果ですが、若干case2のほうが安定しているように見えます。

今回の結論としては、0方程式ならどれでもO.K.という感じでしょうか。

熱が冷めないうちに、他の乱流モデルも書いてみようかなと思ってます。

式の導出

乱流は勉強不足で苦手なので間違ってるかも。それでも、真剣に書くとどんどん長くなるのでメモ程度です。

レイノルズ方程式は次のとおりである。

$\begin{align} \dfrac{\partial u_i}{\partial t} + u_j\dfrac{\partial u_i}{\partial x_j} = X_i - \dfrac{1}{\rho} \dfrac{\partial P_i}{\partial x_i} + \dfrac{1}{\rho} \dfrac{\partial \tau_{ij}}{\partial x_j} + \dfrac{1}{\rho} \dfrac{\partial \tau^{'}_{ij}}{\partial x_j} % + \nu \left( \dfrac{\partial v}{\partial x} + \dfrac{\partial u}{\partial y} \right) \end{align}$

粘性応力 $\tau_{ij}$ 、レイノルズ応力 $\tau^{'}_{ij}$ に後述のモデル式を導入し、若干の式変形を行なうと次式となる。

$\begin{align} \dfrac{\partial u_i}{\partial t} + u_j\dfrac{\partial u_i}{\partial x_j} = X_i - \dfrac{1}{\rho} \dfrac{\partial P_i}{\partial x_i} + \dfrac{\partial }{\partial x_j}\left( (\nu + \nu_t) \left( \dfrac{\partial u_j}{\partial x_i} + \dfrac{\partial u_i}{\partial x_j} \right) - \dfrac{2}{3} k \delta_{ij} \right) \end{align}$

$\nu \ll \nu_t$ のため、 $(\nu+\nu_t) \fallingdotseq \nu_t$ として水深平均化すると前出のcase2となる。

また、 $k$ の項を無視し、 $\nu_t$ を一定として式変形すると次式となる。（わかりやすくx方向のみ）

$\begin{align} \dfrac{\partial u_i}{\partial t} + u_j\dfrac{\partial u_i}{\partial x_j} = X_i - \dfrac{1}{\rho} \dfrac{\partial P_i}{\partial x_i} + \nu_t \left( \dfrac{\partial^2 u}{\partial x^2} + \dfrac{\partial^2 u}{\partial y^2} + \dfrac{\partial^2 u}{\partial z^2} \right) + \nu_t \dfrac{\partial }{\partial x} \left( \dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} + \dfrac{\partial u}{\partial z} \right) \end{align}$

連続式を考慮してこの式を水深平均化するとcase1となる。

ニュートン流体の粘性応力

ニュートン流体の粘性応力 $\tau_{ij}$ (せん断応力+垂直応力)は、速度勾配、粘性係数を用いて次のように示される。

$\begin{align} \tau_{ij} &= \mu \left( \dfrac{\partial u_j}{\partial x_i} + \dfrac{\partial u_i}{\partial x_j} \right) \\ \dfrac{\tau_{ij}}{\rho} &= \nu \left( \dfrac{\partial u_j}{\partial x_i} + \dfrac{\partial u_i}{\partial x_j} \right) \end{align}$

ここに、流体の粘性係数 $\mu[\rm{Pa·s}$ ] 、流体の密度を $\rho[\rm{kg/m^3}$ ]、動粘性係数 $\nu=\mu/\rho[\rm{m^2/s}$ ]とする。(わかりやすいように単位を記載)

参考資料

古典的な流体力学、水理学の教科書。なんでもO.K.。例えば以下。

リンク

webだとこのあたりhttp://lemons.k.u-tokyo.ac.jp/SATO/lecture/FMC/6.pdf

ブシネスクの渦粘性近似

ブシネスクの渦粘性近似では、レイノルズ応力 $\tau^{'}_{ij}$ は上記の粘性応力の類推から次式で表す。

$\begin{align} \tau^{'}_{ij} &= -\rho\overline{u^{'}_i u^{'}_j} = \mu_t \left( \dfrac{\partial u_j}{\partial x_i} + \dfrac{\partial u_i}{\partial x_j} \right) - \dfrac{2}{3}\rho k \delta_{ij} \\ \dfrac{\tau^{'}_{ij}}{\rho} &= -\overline{u^{'}_i u^{'}_j} = \nu_t \left( \dfrac{\partial u_j}{\partial x_i} + \dfrac{\partial u_i}{\partial x_j} \right) - \dfrac{2}{3} k \delta_{ij} \\ \end{align}$

ここに、渦粘性係数 $\mu_t[\rm{Pa·s}$ ] 、流体の密度を $\rho[\rm{kg/m^3}$ ]、渦動粘性係数 $\nu_t=\mu_t/\rho[\rm{m^2/s}$ ]とする。(わかりやすいように単位を記載)

上式の右辺の最後の項は次項に示す乱れエネルギー $k = \dfrac{\overline{u^{'}_i u^{'}_i}}{2}$ を満たすために付加している。(補足参照)

補足：レイノルズ応力の付加項

レイノルズ応力を

$\begin{align} -\overline{u^{'}_i u^{'}_j} = \nu_t \left( \dfrac{\partial u_j}{\partial x_i} + \dfrac{\partial u_i}{\partial x_j} \right) % - \dfrac{2}{3}\rho k \delta_{ij} \end{align}$

と仮定して、乱れエネルギーの式に代入すると次式となり $k$ が0となる。

$\begin{align} k=-\nu_{t}\left(\dfrac{\partial u}{\partial x}+\dfrac{\partial v}{\partial y}+\dfrac{\partial w}{\partial z}\right) \\ \dfrac{\partial u}{\partial x}+\dfrac{\partial v}{\partial y}+\dfrac{\partial w}{\partial z}=0 \end{align}$

これを満たすためにレイノルズ応力に $- \dfrac{2}{3}\rho k \delta_{ij}$ を付加する。

参考資料

書籍

リンク

web
- http://penguinitis.g1.xrea.com/study/note/turbulence_model.pdf#page=14
- http://www.opencae.or.jp/wp-content/uploads/2015/06/%E7%AC%AC5%E5%9B%9E_%E8%BC%AA%E8%AC%9B_%E6%95%B0%E5%80%A4%E6%B5%81%E4%BD%93%E5%8A%9B%E5%AD%A6.pdf

河川流の乱れエネルギー

乱れエネルギー $k$ は速度の乱れ成分 $\overline{u^{'}},\overline{v^{'}},\overline{v^{'}}$ を用いて次のように定義される。

$\begin{align} k = \dfrac{1}{2}( \overline{u^{'} u^{'}} + \overline{v^{'} v^{'}} + \overline{w^{'} w^{'}}) \end{align}$

$\overline{u^{'}},\overline{v^{'}},\overline{v^{'}}$ の水深方向分布は、禰津らの実験結果より次のように示される。

$\begin{align} \dfrac{ \sqrt{ \overline{u^{'}}^2 } }{U_*} &= {\rm A_u}e^{-z/h} \\ \dfrac{ \sqrt{ \overline{v^{'}}^2 } }{U_*} &= {\rm A_v}e^{-z/h} \\ \dfrac{ \sqrt{ \overline{w^{'}}^2 } }{U_*} &= {\rm A_w}e^{-z/h} \end{align}$