Incompressible Fluid FLow Solver For Equil-Spaced Rectilinear Mesh

October 7, 2022 · View on GitHub

等間隔直交格子のための非圧縮性流体ソルバ.

現在の状態

移流項と粘性項の評価を選択可能にした. 陰解法には対応出来ていない.
並列化に対応.(OpenMP)

概要

等間隔直交格子で流れ計算を行う.

離散化手法

有限体積法(コロケート格子)

移流項: 1次精度風上法
粘性項: 2次精度中心差分法
時間精度:
- solver_fs_t : 1次精度Euler陽解法
- solver_fs_imp_t : 粘性項のみクランク･ニコルソン法
- その他, ユーザが各自で用意する.

解法

フラクショナルステップ法

解析手順

ビルド方法

プログラムの構成

本ソルバはライブラリ群libとメインプログラムに分けられている. libにあるファイル群は数値シミュレーションに必要なモジュール群であり, これらから静的ライブラリを生成する.

メインプログラムは各々のケースに応じて用意し, コンパイル時に上記のライブラリとリンクすることで実行ファイルを生成する.

プログラムの取得

gitの場合 git clone https://github.com/hirotori/incomp3D.git

gitがない場合はDownload Zipから.

ライブラリのビルド

libのビルドにはcmakeを利用する.

まず初めに, ビルド用のディレクトリを作成しておく.
次にconfigureを行う. コンパイラがgfortranの場合は以下のようにする.

cmake .. -DCMAKE_Fortran_COMPILER=gfortran

コンパイラオプションを適用する場合は上に-DCMAKE_BUILD_TYPE=を指定する. デバッグモードはdebug, リリースの場合はreleaseとする. OpenMPによる並列化を行う場合は-Duse_OpenMP=ONとする. 指定しない場合は並列化は行われない.
cmakeによるconfigureが完了後, makeコマンドを実行することでコンパイルを実施する. 終了後アーカイブlibincomp3d.a, libincomp3d_util.aが生成される.

メインプログラムのコンパイル

例題として3つのサンプルケースを用意している.ここでは例としてsampleのコンパイル方法を説明する.

./srcにはメインプログラムが入っている. 中身は本ライブラリの使用例で, 解析ケース, アルゴリズム, シミュレーターの3つのユーザ派生型でメインプログラムを構成している.

メインプログラムのビルドはmakeで実施することになっている. サンプルのMakefileを用意しているのでそちらを参考に.

注意ライブラリとメインプログラムのコンパイラは同じである必要がある. libのコンパイルをgfortranで行った場合, メインプログラムのコンパイルも同じコンパイラで実施すること. そうしないとコンパイルエラーが出る.

プログラムの実行

解析にあたって, まず計算条件設定ファイルconfig.txtを編集する.

メッシュの生成

ファイルの2~3行目.

33 33 2  !grid points (imax,jmax,kmax)
3.141592653 3.141592653 0.01  !length (x,y,z)

grid points は計算領域の格子点の個数を指定する.

lengthは計算領域の寸法を指定する.

計算条件

注釈通りに設定する.

境界条件

次節を参照.

境界条件

境界条件は圧力, 速度それぞれに指定する. 現在指定可能な条件は境界値固定, 境界勾配固定, 周期境界, バッファ付き周期境界の4種類である.

境界条件はプログラム内部で以下のように番号付けされている.

bc_fix_val = 1              !固定値境界
bc_fix_grad = 2             !固定勾配
bc_fix_periodic = 3         !周期境界
bc_fix_periodic_buffer = 4  !バッファ付き周期境界

境界条件はコンフィグファイルconfig.txtで指定する. フォーマットは以下の通りとなる.

3 3 0.0 0.0 0.0 0.0  !B.C. for i = 1    face. [bc_u, bc_p, u, v, w, p]
3 3 0.0 0.0 0.0 0.0  !B.C. for i = imax face. [bc_u, bc_p, u, v, w, p]
3 3 0.0 0.0 0.0 0.0  !B.C. for j = 1    face. [bc_u, bc_p, u, v, w, p]
3 3 1.0 0.0 0.0 0.0  !B.C. for j = jmax face. [bc_u, bc_p, u, v, w, p]
2 2 0.0 0.0 0.0 0.0  !B.C. for k = 1    face. [bc_u, bc_p, u, v, w, p]
2 2 0.0 0.0 0.0 0.0  !B.C. for k = kmax face. [bc_u, bc_p, u, v, w, p]

境界面の値(u,v,w,p)は, 固定値境界の場合は速度/圧力として, 固定勾配の場合は速度勾配/圧力勾配として扱われる.

数値シミュレーションで扱う境界条件の組み合わせは限られている. 以下では主な組み合わせについて説明する.

1. 流入条件

境界面に固定された速度（一様流）を与える. 通常,圧力はNeumann境界となる.

例: 一様流入

速度 $\bm{u}=(1,0,0)$ で流入する条件の場合は以下のように指定する.

1 2 1.0 0.0 0.0 0.0 ![u_bc, p_bc, u,v,w,p]

圧力をゼロ勾配とするため, pに0.0を指定している.

2. 流出条件

対流流出条件を与える. ここでは流れ方向の速度の微分がゼロとする. 圧力は境界面で固定する.

2 1 0.0 0.0 0.0 1.0 ![u_bc, p_bc, u,v,w,p]

速度にゼロ勾配を与えるため, u,v,wを全てゼロにしている.

注意1 勾配規定条件は面に垂直な勾配として考える. そのため, 例えばi=1,imax の面に対しては境界速度成分はuのみが有効. ただv,wにゼロ以外を入れると狂うので全てゼロにしてください.

注意2 圧力を任意に指定することが出来るが, 普通は圧力の初期値と対応させるべき.

3. 滑り無し条件

滑り無し条件を課す. 境界面で速度をゼロに固定する. 圧力は流入条件と同様にNeumann境界となる.

1 2 0.0 0.0 0.0 0.0 ![u_bc, p_bc, u,v,w,p]

注意速度はゼロ以外でも可能. 物理的に適切な値を考えて入れてください.

4. 滑り条件

速度をゼロ勾配, 圧力もゼロ勾配で実現する.

5. 周期境界条件

2つの境界面に指定することで周期境界条件を達成する.

注意1 可能なペアは i=1とi=imx, j=1とi=jmx, k=1とk=kmxの三通りのみ. 注意2 圧力については普通は周期境界条件とする.

6. 非定常流出条件

まだ作ってない. 移流方程式に基づいて流出条件を設定する.

7. バッファ付き周期境界条件

周期境界として参照する速度を任意の場所に指定する.

注意例えばi=1をこれにしたとき, i=imxについては周期境界ではなく流出条件とするべき. 例:

4 2 0.0 0.0 0.0 0.0  !B.C. for i = 1    face. [bc_u, bc_p, u, v, w, p]
2 1 0.0 0.0 0.0 1.0  !B.C. for i = imax face. [bc_u, bc_p, u, v, w, p]

その他の注意

サンプルケースの管理

サンプルケースはgitで管理されているため, git cloneしたリポジトリのサンプルケースを色々いじってからgit pullするとマージコンフリクトが生じてpull出来ない可能性がある. これを避けるには

stashする. 変更を退避させて, pullしてから元に戻す. それでもローカルでコンフリクトはおきるが
(推奨) sampleフォルダをコピーして違うところで作業する. これが一番確実. その場合, Makefile内のインクルードパスなどは各自で調整する.(絶対パスに置き換えるなどして)

cmake

コンパイルオプションを変更したい場合は, 再度ライブラリのビルド手順を行う. その際に-DCMAKE_BUILD_TYPEフラグを変更する.

上記が上手くいってなさそうな場合は, 一旦buildディレクトリを削除してから, 再度mkdir build→ビルド.

ソルバの仕様

速度半陰解法(fs_solver_imp_t)の反復の閾値が今のところ1.0e-16となっている. どの問題もここまで必要かは不明. ただ, 現在の解法ではsample3のポアズイユ流の圧力分布において振動が除去される.

閾値など変更したい場合は型束縛手続き setting_parameterで変更する.