Pythonで書かれたソルバー[NaysMini]のページへようこそ!

はじめに

NaysMiniはiRICのソルバーであるNays2DHの流れの計算部分のみを直交座標系 で表し, 一般公開用にPython言語で書き直したものである。 ソースコードがすべてオープンになっているので、開水路の流れの数値計算をPythonコードと 一緒に学習用教材として利用にも適している。

なお、Python言語によるNaysMiniのiRICへの搭載に際しては、(株)River Link の旭一岳氏、(株)みすほ情報総研 の井上敬介氏に多大なるご協力を頂いた。また、2次元流れ計算のPythonコードについては北海道大学大学院生の 高橋広大君と岡安努君より協力を頂いた。ここに記して謝意を表する。

2021年6月22日 清水康行

基礎式

NaysMiniでは以下の2次元浅水流式を用いる。

\[\frac{\partial h}{\partial t}+\frac{\partial (uh)}{\partial x}+ \frac{\partial (vh)}{\partial y}=0\]

\[\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})\]

\[\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})\]

ただし, \(x, y\) は互いに直交する平面座標軸, \(t\) は時間, \(u, v\) は \(x, y\) 方向の水深平均流速, \(h\) は水深, \(H\) は水位, \(g\) は重力加速度, \(\tau_x, \tau_y\) は \(x, y\) 方向の河床せん断力, \(\rho\) は水の密度, \(\nu_t\) は渦動粘性係数である。

抵抗則にマニング則を用いると, 河床せん断力は下記で表される。

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

ただし, \(n\) はマニングの粗度係数. 渦動粘性係数は、ゼロ方程式モデルを採用し、下記で表す。

\[\nu_t = \frac{\kappa}{6} u_\ast H\]

ただし、\(u_\ast\) は摩擦速度

\[u_\ast=\sqrt{ghI_f}, \ \ \ I_f=\frac{n^2 (u^2+v^2)}{h^{4/3}}\]

ただし、\(I_f\) は摩擦勾配である。

数値計算法

計算は直交格子のセルセンターに水深の計算点、格子の辺の流速の計算点と配置する スタカード格子による差分法で行われる。タイムステップ毎に、 運動方程式は移流部分を非移流部分に分ける分離解法を 用い、移流部分はCIP法を、非移流部分は運動方程式と連立して \(u, v, h\) が同時に満たされる ように繰り返し計算を行う。

計算手順

計算手順は下記の通りである。

1. 計算格子の作成

2. 障害物の挿入

3. 計算条件の設定

4. 計算の実施

5. 計算結果の表示

計算事例

ここでは直線水路における障害物を含む流れの計算例を示す。

iRICの起動とソルバーの選択

iRICのオープニング画面で、[新しいプロジェクト]を選択(Figure 1 )

: 新しいプロジェクト

Figure 2 ソルバーの選択Windowで[NaysMini(Pythonで書かれたNays2dHの簡易版)]を選択し、 [OK]OKを押す。

: ソルバーの選択

計算格子の作成

メインメニューから[格子]→[格子生成アルゴリズムの選択]を選ぶ

: 格子生成(1)

[格子生成アルゴリズムの選択]ウィンドウで[簡易直線・蛇行水路生成ツール]を選び[OK]を押す。

: 格子生成(2)

[格子生成]ウィンドウの[水路形状]タブでパラメータを Figure 5 のように設定し、[格子生成]を押す。

: 格子生成(3)

例によってしつこく「マッピングしますか？」と聞かれるので、[はい]を押す。

: 格子生成(4)

Figure 7 のような格子が作成される。

: 格子生成完了

障害物セルの指定

ここでは流路の中央上流側に長方形の障害物を設置する。 Figure 8 に示すように, オブジェクトブラウザーの[セルの属性]および[障害物セル]に☑マークを入れ、 Figure 9 に示すように障害物にしていしたいセルをマウスで囲い、 右クリックして現れる[障害物セルの編集]ウィンドウで、[障害物セル]に指定し[OK]を押す。

:障害物の設定(1)

:障害物の設定(2)

:障害物の設定(完了)

計算条件の設定

メインメニューから[計算条件]→[設定]を選択。(Figure 11 )

:計算条件の設定(1)

[計算条件]ウィンドウで、[水理条件および物理定数]を Figure 12 のように設定する。

:計算条件の設定(2)

[計算条件]ウィンドウで、[境界条件]を Figure 13 のように設定する。

:計算条件の設定(3)

[計算条件]ウィンドウで、[時間および繰り返し計算パラメーター]]を Figure 14 のように設定し、[保存して閉じる]をクリック。

:計算条件の設定(4)

計算の実行

メインメニューから[計算]→[実行]を選択。

:計算の実行(1)

プロジェクトを保存するか聞かれるので通常は[はい]を選択して、プロジェクトを保存する。

:計算の実行(2)

Figure 17 のコンソールウィンドウが出て計算が実行される。

:計算の実行(3)

計算が終了すると「ソルバーの計算が終了しました」と表示されるので、[OK]を押す。

:計算の実行(4)

計算結果の表示

メインメニューから[計算結果]→[新しい可視化ウィンドウ(2D)を開く] を選択する。

:計算結果の表示(1)

[可視化ウィンドウ(2D)]が表示される。

:計算結果の表示(2)

[オブジェクトブラウザー]で[セル属性]と[障害物セル(障害物セル)]に☑マークを入れると、 障害物の部分が Figure 21 のように色が変わる。

:計算結果の表示(3)

[オブジェクトブラウザー]で[スカラー(格子点)]と[Vorticity]に☑マークを入れて、[Vorticity]を右クリックして、 [プロパティ]を選択する( Figure 22 )。なお。[Vorticity]は 渦度 のことで、次式で求められる水平渦の強度である。

\[\Omega = \frac{\partial u}{\partial y}-\frac{\partial v}{\partial x}\]

:計算結果の表示(4)

[スカラー設定]ウィンドウで、[物理量:], [値の範囲]の[自動]の前にある☑マークを外し、 [最大値]および[最小値]をそれぞれ[2]と[-2]に設定する( Figure 23 ).

:計算結果の表示(5)

同じ[スカラー設定]ウィンドウで、[カラーマップ]の[手動]ボタンを押して[設定]をクリックして 現れる[カスタムカラー]ウィンドウで[種類]を[3色]にして、[最大値を]赤に、最小値を[青] にして[OK]を押す。

:計算結果の表示(6)

[スカラー設定]ウィンドウに戻って、[半透明]の前の☑ボックスからチェックをはずして、[OK]を押す。

:計算結果の表示(7)

オブジェクトブラウザーの[スカラー]の[Vorticity]をクリックし、 カラーマップをドラッグして適当な位置に移動する。さらに、オブジェクトブラウザーで、 [時刻]を右クリックして[プロパティ]を表示し。[フォント]のサイズを適当に大きくする。

:計算結果の表示(9)

メインメニューから[アニメーション]→[実行]を選択すると、渦度とパーティクルの動画が始まる。

:計算結果の表示(10)

同様に、下記 Figure 28 の手順で、流速ベクトル、流速コンターのアニメーションを 表示できる。

:計算結果の表示(11)

アニメーションファイルの作成

メインメニューから[ファイル]→[連続スナップショット/動画/Google Earth出力]を選択

:アニメーションのファイルの作成(1)

[イントロダクション]が表示されるので[次へ]を押す。

:アニメーションのファイルの作成(2)

[ウィンドウの選択]が表示されるので、アニメーションにしたいウィンドウを選んで[次へ]を押す・

:アニメーションのファイルの作成(3)

[ファイル属性]ウィンドウが表示されるので、ファイルを保存するフォルダ、静止画像の形式、連番の桁数など指定して [次へ]を押す。

:アニメーションのファイルの作成(4)

[動画ファイルの設定]ウィンドウで、[動画ファイルを出力する]に☑マークを入れる。 動画は[動画の長さ]もしくは[1秒あたりのフレーム数]を指定できるが、 ここでは[動画の長さ]を[10秒]に指定し、[次へ]を押す。

:アニメーションのファイルの作成(5)

[タイムステップ設定]では開始時間、終了時間、間引きなどの指定が出来るので、適宜指定して[次へ]を 押す。

:アニメーションのファイルの作成(6)

iRICでは、アニメーションをGoogle Earthに出力することが出来る。必要であれば [Google Earthに出力する]に☑マークを入れる。不要であれば☐マークは入れないで、 [次へ]を押す。

:アニメーションのファイルの作成(7)

最後に生成されるフィルの一覧が表示されるので、確認して良ければ[完了]を押す。

:アニメーションのファイルの作成(8)

保存の進行状況がパーセント表示で示され、終了すると、 指定したフォルダに連番のイメージファイルと動画ファイル(mp4)形式が生成される。

:アニメーションのファイルの作成(9)

動画をmp4以外の形式に変換したい場合は ffmpeg など が便利である。例えば、mp4からgifアニメに変換する場合はコンソール画面で:

ffmpeg -i img.mp4 -vf scale=640:-1 -r 30 -loop 0 img.gif

打ち込めば下記のアニメーションgifが生成される。

:アニメーションファイル