こんにちは、ぴよ工房を運営しているぴよ(@piy0_gadget)です！
この記事では、
- 二次元バックステップ流れの概要
- CFDコードの構築
- 結果の比較
について話していきます。
本記事では過去の記事である『二次元チャネル流のCFDコード構築』で作成したプログラムコードを用いて、二次元バックステップ流れを計算していきます。
-
-
【即実装】二次元チャネル流のCFDコード構築【コード例付き】
こんにちは、ぴよ工房を運営しているぴよ(@piy0_gadget)です！ この記事では、 二次元チャネル流のCFDコード構築 の解説をしていきます。 本記事の一番の目的は『即実装』です！ ...
続きを見る
どんな流れかイメージしやすいように動画を作ったのでまずはこちらをご覧ください。
(今回から新しく動画作成用のルーチンを追加しました！！)
流れが剥離して、地面で再付着する様子が見れますね！
今回はこのような流れを取り扱います！
(ちなみに、このバックステップ流れの計算コード構築が出来ると、応用して角柱まわりの流れの計算を簡単に行うことができカルマン渦の可視化をすることが出来ます。)
なお、本記事のコードはFortranで作っています。
バックステップ流れとは
動画を見ると流れが何となく分かったと思いますが、もう少し詳しく流れの説明をしたいと思います。
画像のように左端から流れが流入し、段差を過ぎると流れが剥離して、そして地面に流れが付着(再付着と一般的には言われます)する流れをバックステップ流れと言います。
文献を見ると、ダクト内の流れがバックステップ流れとなっているみたいですね。
ただ、「ダクト」と画像検索してもバックステップ流れっぽくはなかったです。おそらくですが、流れをできるだけ綺麗に保つために剥離しない構造にしてあるのだと思います。(実際どうなのかは分かりませんが・・・)
CFDコードの構築
それではCFDコードの構築をしていきます。
過去の記事で構築したコードを用いて、計算領域や格子、初期条件と境界条件をバックステップ流れに合わせて変更するだけなので簡単です。
ただし、物体回りを流れる計算に近いためその取扱いが少しだけ工夫が必要です。
基本となるコードはGithubで公開しているので、もしコードが無い場合はここから入手して下さい。
計算領域と格子数
計算領域と格子数は下記のようにします。
|
1 2 3 4 |
Real*8 ,Parameter :: Lx = 12.0d0 Real*8 ,Parameter :: Ly = 2.0d0 Integer,Parameter :: Nx = 121 Integer,Parameter :: Ny = 21 |
領域は計算の検証のため、カマキリさんの記事の計算結果と合わせました。
格子数は割と適当ですが、多すぎると計算負荷が大きくなるので少なめにしてあります。
物体の生成
今回の計算は一応物体まわりの流れになります。
そのため、物体の生成のために物体パラメータを下記のように定義します。
|
1 2 3 4 5 6 |
Real*8 ,Parameter :: Lx_sta = 0.0d0 Real*8 ,Parameter :: Lx_end = 2.0d0 Real*8 ,Parameter :: Ly_sta = 0.0d0 Real*8 ,Parameter :: Ly_end = 1.0d0 Integer i_Lxsta, i_Lxend Integer j_Lysta, j_Lyend |
物体パラメータを定義する意味としては、計算領域においてどの場所に物体が位置するかを決定するためです。
そして、このパラメータに沿って物体まわりにおける境界条件を計算します。
以下にパラメータの簡単な説明をまとめておきます。
| Lx_sta | x方向における物体の始点座標 |
|---|---|
| Lx_end | x方向における物体の終点座標 |
| Ly_sta | y方向における物体の始点座標 |
| Ly_end | y方向における物体の終点座標 |
| i_Lxsta | Lx_staに対応する配列番号[ i ] |
| i_Lxend | Lx_endに対応する配列番号[ i ] |
| j_Lysta | Ly_staに対応する配列番号[ j ] |
| j_Lyend | Ly_endに対応する配列番号[ j ] |
物体パラメータの計算
物体パラメータ(i_Lxsta,i_Lxend,j_Lysta,j_Lyend)を以下のサブルーチンで計算します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
Subroutine Calc_range use Definition Implicit none Integer i,j,k !x-range Do i=0,Nx If ( Lx_sta <= x(i) ) Then i_Lxsta = i Exit End If End Do Do i=0,Nx If ( Lx_end <= x(i) ) Then i_Lxend = i Exit End If End Do !y-range Do j=0,Ny If ( Ly_sta <= y(j) ) Then j_Lysta = j Exit End If End Do Do j=0,Ny If ( Ly_end <= y(j) ) Then j_Lyend = j Exit End If End Do Return End Subroutine |
初期条件と境界条件
初期条件は、初めは流れが無い状態にしたく下記のように定義しました。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
Subroutine Ini_Condition use Definition Implicit none Integer i,j,k !--- Initial_Condition --- Do j=0,Ny Do i=0,Nx U(i,j) = 0.0d0 V(i,j) = 0.0d0 P(i,j) = 0.0d0 End Do ; End Do Return End Subroutine |
また境界条件に関しては、左端から流入、右端で流出、上下は滑りなし条件としています。物体内に関しては速度が0で、物体壁は圧力勾配が0としています。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 |
Subroutine Bou_Condition use Definition Implicit none Integer i,j,k !--- Left --- Do j=0,Ny U(0,j) = 1.0d0 V(0,j) = 0.0d0 P(0,j) = P(1,j) End Do !--- Right --- Do j=0,Ny U(Nx,j) = U(Nx-1,j) V(Nx,j) = V(Nx-1,j) P(Nx,j) = P(Nx-1,j) End Do !--- Up --- Do i=0,Nx U(i,Ny) = 0.0d0 V(i,Ny) = 0.0d0 P(i,Ny) = P(i,Ny-1) End Do !--- Down --- Do i=0,Nx U(i,0) = 0.0d0 V(i,0) = 0.0d0 P(i,0) = P(i,1) End Do !--- in object --- Do j=j_Lysta,j_Lyend ; Do i=i_Lxsta,i_Lxend U(i,j) = 0.0d0 V(i,j) = 0.0d0 End Do ; End Do Do j=j_Lysta,j_Lyend P(i_Lxsta,j) = P(i_Lxsta+1,j) End Do Do i=i_Lxsta,i_Lxend P(i,j_Lysta) = P(i,j_Lysta+1) End Do Return End Subroutine |
出力
出力に関してもカマキリさんの記事を参考にコードを書きました。
下記にコードを載せておきますが、詳細はこちらの記事を参考にしてください！
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 |
Subroutine output_vtk(step) use Definition Implicit none Integer i,j,k Integer step, points character(len=40) filename Write(filename,"(a,i5.5,a)") "output/output",int(step),".vtk" Open(10,file=filename) Write(10,"('# vtk DataFile Version 3.0')") Write(10,"('test')") Write(10,"('ASCII ')") Write(10,"('DATASET STRUCTURED_GRID')") Write(10,"('DIMENSIONS ',3(1x,i3))") Nx+1, Ny+1, 1 Write(10,"('POINTS ',i9,' float')") (Nx+1)*(Ny+1)*1 Do j=0,Ny ; Do i=0,Nx Write(10,"(3(f9.4,1x))") i*dx, j*dy, 0.0d0 End Do ; End Do Write(10,"('POINT_DATA ',i9)") (Nx+1)*(Ny+1)*1 !date input Write(10,"('SCALARS U float')") Write(10,"('LOOKUP_TABLE default')") Do j=0,Ny ; Do i=0,Nx Write(10,*) U(i,j) End Do ; End Do Close(10) Return End Subroutine |
その他のパラメータ
その他のパラメータ(刻み時間、レイノルズ数)は以下のように定義します。
|
1 2 |
Real*8 ,Parameter :: dt = 0.005d0 Real*8 ,Parameter :: Re = 150.0d0 |
それでは計算を回してみましょう！！
計算結果
定常状態になったときの流れ方向速度Uの可視化を載せます。
やはり可視化すると気持ちいいですね！！
再付着位置もカマキリさんの結果と同じそうなので、間違ってはないかなーと思います。
まとめ
本記事では『二次元バックステップ流れのCFDコード構築』を行いました。
キャビティ流れもそうですが、雛型のコードがあるとかなり実装は楽ですね。
-
-
【数値計算】二次元キャビティ流れのCFDコード構築
こんにちは、ぴよ工房を運営しているぴよ(@piy0_gadget)です！ この記事では、 二次元キャビティ流れの概要 CFDコードの構築 結果の比較 について話していきます ...
続きを見る
物体パラメータの実装も出来たので、このコードを応用すれば更に計算の幅が広がります！
ぜひコードをいじって遊んでみて下さい！
疑問があればリプやDMで気軽にどうぞ！
https://twitter.com/piy0_gadget/status/1263602395478495237?s=20