カルマン渦を計算機上で再現したい！！
今回のモチベーションはこれにつきます。
物体まわりを通りすぎる流れは規則的な渦列を発生させ、この渦列をカルマン渦と言いますが、度々、話題に上がっている現象なので馴染みがある人も多いと思います。
気象衛星からもカルマン渦を確認することができ、地球規模で発生しているので面白いですよね。
下の画像は気象衛星ひまわりから得られた写真です。
2019年12月12日に済州(チェジュ)島付近で発生していて、こちらのサイトから見ることが出来ます。
このサイトではリアルタイムに雲の様子を見ることもできて非常に楽しいですよ。
それでは、このカルマン渦を数値計算で再現していきます。
計算には以前の記事で用いたFortranコードを用いて実装していきます。
実装の詳細は記事の中で説明するので、ぜひ一緒にコードを構築してカルマン渦を可視化してみましょう！
-
-
【数値計算】二次元バックステップ流れのCFDコード構築【動画付】
こんにちは、ぴよ工房を運営しているぴよ(@piy0_gadget)です！ この記事では、 二次元バックステップ流れの概要 CFDコードの構築 結果の比較 について話していきます。 本記事 ...
続きを見る
結果にはカルマン渦への遷移動画もあるよ！！
Contents
計算の概要
今回の計算の概要図です。
計算領域内に角柱を設置し、左から流れが入り角柱まわりを流れ、右へ流れが出ていきます。
角柱の下流でカルマン渦が発生するはずなので、早く計算して挙動を見たいですね！！
CFDコードの構築
それではCFDコードの構築をしていきます。
バックステップ流れで構築したコードを用いて計算するため、計算領域・格子数、物体生成、初期条件・境界条件を今回の計算に合わせて変更するだけなので簡単です。
基本となるコードはGithubにて公開しているので、コードが無い場合はここから入手しコードを書き換え下さい。
さらに詳しく
計算領域と格子数
計算領域と格子数は下記のようにします。
|
1 2 3 4 |
Real*8 ,Parameter :: Lx = 25.0d0 Real*8 ,Parameter :: Ly = 15.0d0 Integer,Parameter :: Nx = 251 Integer,Parameter :: Ny = 151 |
領域の大きさは割と適当ですが、小さすぎるとカルマン渦が計算領域内に収まらず流れを正しく再現できない可能性があります。
ただ、計算負荷の観点から少し小さめに設定しました。
格子数は以前の計算(キャビティ流れ、バックステップ流れ)と比べると多めにしてあります。
これは計算領域が広くなり、それに伴い格子幅も広くなったため格子数を増やしました。
さらに詳しく
物体生成
それでは計算領域内に角柱を設置します。
今回の計算では角柱の計算領域内における位置を与えて、角柱を設置するようにしてあります。
この与える値(角柱の位置)を本記事では物体パラメータと呼んでいます。
物体パラメータは下記のように与えます。
|
1 2 3 4 5 6 |
Real*8 ,Parameter :: Lx_sta = Lx/5.0d0 - 0.5d0 Real*8 ,Parameter :: Lx_end = Lx/5.0d0 + 0.5d0 Real*8 ,Parameter :: Ly_sta = Ly/2.0d0 - 0.5d0 Real*8 ,Parameter :: Ly_end = Ly/2.0d0 + 0.5d0 Integer i_Lxsta, i_Lxend Integer j_Lysta, j_Lyend |
Lx/5.0d0、Ly/2.0d0は角柱の中心位置を表しており、その中心から0.5離れた位置を角柱の外縁としています。
今回は角柱のx方向の長さは1で、y方向の長さも1であることを表しています。
そのため、アスペクト比(縦横比)は1:1となっており、この0.5を書き換えることで任意のアスペクト比の角柱を生成することが出来ます。
物体パラメータの簡単な説明をまとめておきますので、理解の助けとなれば幸いです。
| 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 ] |
初期条件と境界条件
初期条件は下記のように定義しました。
|
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) = 1.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) = 1.0d0 V(i,Ny) = 0.0d0 P(i,Ny) = P(i,Ny-1) End Do !--- Down --- Do i=0,Nx U(i,0) = 1.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 |
Real*8 ,Parameter :: dt = 0.01d0 Real*8 ,Parameter :: Re = 100.0d0 |
ただ、今回はレイノルズ数によって流れがどのように変わるか見てみたかったのでレイノルズ数150の条件でも計算しました。
それでは計算を回してみましょう！！
計算結果
カルマン渦が再現できているか見てみましょう！
なお、可視化は主流方向速度を表しています。
レイノルズ数100における計算結果
動画を見るとカルマン渦を再現できていそうですね！！！
遷移の前半は流れがじわーっと下流まで伸びていて、しばらく流れは安定していますが段々と流れがぐわんぐわんして最終的にはカルマン渦が形成されてるのが分かります。
カルマン渦が形成されてからは、流れがうねうねして面白いですね。
レイノルズ数150における計算結果
こちらもカルマン渦を再現できました。
動画を見る限り遷移の傾向はレイノルズ数が100と同じっぽいですね。
ただし、やはりレイノルズ数が大きい方がカルマン渦の形成にかかる時間は短いです。
これはレイノルズ数が小さいと粘性の影響で、流れが周囲へじわじわと拡散してしまうことが原因だと思います。
(もし間違っていたらTwitterで教えて下さい)
レイノルズ数100と150の結果の比較
結果の比較をするためRe=100,150の7000stepにおける可視化を載せておきます。
これを見ると、レイノルズ数が大きい方が流速は大きくなっていますね(色がより濃いことから)。
挙動の違いとしては、レイノルズ数が大きいとカルマン渦がより丸みを帯びていますね。これは、レイノルズ数が小さいと粘性の影響で渦の端が周囲へ拡散しているから形を保てずこういった形状になるのだと思います。
あと、レイノルズ数が大きい方が渦が細かくなっています。
渦の縦方向の大きさは同程度ですね。これは角柱のアスペクト比が影響してそうです。
あと気になる点としては、角柱の前方ではガタガタしている点です。原因は色々あると思いますが、一番は格子数が足りてないことですかね。
まぁ、何はともあれカルマン渦を計算機上で再現できました！！
まとめ
本記事では『角柱まわりの数値計算』を行いました。
カルマン渦という実際の現象を、数値的にも求めることが出来るのは本当に面白いですね。
アスペクト比を変更して流れの挙動をみたり、物体パラメータを増やして複数の角柱で計算してみても面白そうです。
ぜひコードをいじって遊んでみて下さい！
計算コードの詳しい説明は以前の記事で解説してあるので、興味ある方はこちらもご覧あれ！
さらに詳しく
疑問があればリプやDMで気軽にどうぞ！