Fortran 数値計算 流体力学

【即実装】二次元チャネル流のCFDコード構築【コード例付き】

2020年5月18日

 

こんにちは、ぴよ工房を運営しているぴよ(@piy0_gadget)です!

この記事では、

  • 二次元チャネル流のCFDコード構築

の解説をしていきます。

本記事の一番の目的は『即実装』です!

そのため、詳しい原理の説明はせずに進めていきます。原理・原則を理解することは大切ですが、まずは大まかな概要を把握するために、この記事を参考にしてコードの構築をしてください。

ただし、ただコードの真似して実装しても自分のスキルとして身につかないため、各ルーチンの計算では一体何をやっているのか・この計算の目的とその結果は何か、について丁寧に説明します。

つまり、大まかな理解で実装して、詳しい知識は後から肉付けをしてもらう狙いです。

このチャネル流のコードの構築ができれば、そのコードを基にして様々な拡張(流れを三次元にしたり、境界条件を変えれば他の流れ場もできます)が出来るため、まずは本記事でコード実装を目指しましょう!

本記事ではFortranでコードを書いています。Fortranのインストール方法は過去の記事(Fortranのインストール方法)で解説しているため、インストールしてない場合は参考にしてみて下さい!

また、本記事で作成したコードはGitHubにて公開もしているので、プログラム全体を見てみたいって人は参考にしてみて下さい!

 

それではやっていこー!!

計算の概要

いきなりチャネル流のコード実装に取り組む前に、まずは計算の概要について話していきたいと思います。

必要最低限の事しか話さないので、さらに知りたい人は各自調べてください。

チャネル流とは

まず、計算対象にしているチャネル流の説明をします。とは言っても図を見た方がイメージ出来ますよね。


流れ場のイメージは出来ましたか?

左から流れが流入し、上下が板(壁)の間を流れ、右へ流れが流出します。比較的単純な流れ場ですが、流入・流出条件や壁面境界条件が含まれており、参考になる箇所は多いかと思います。

ちなみに、チャネル流とは平行平板間流れとも呼ばれます。言葉通りですね笑

支配方程式( 質量保存とナビエ・ストークス方程式 )

計算の支配方程式は『質量保存の式』『(非圧縮性)ナビエ・ストークス方程式』です。流体の勉強をやったことがある人なら聞いたことがあるかもしれませんね。

支配方程式の分かりやすい説明は『このサイト』が参考になると思います。

下記に支配方程式を記述しますが、今回は二次元であるため式が少ないことに注意して下さい。
なお、流速のx,y方向の成分をそれぞれu,vとしています。

  • 質量保存の式

  • ナビエ・ストークス方程式

さらに、この式から無次元化を施します。無次元化をすることで、流れを支配するパラメータによって流れ場がどのように変化するかを調査できます。

代表長さと代表速度で無次元化をすると以下の式になります。この式変形についてはいずれやりたいと思います(多分)。式中の『*』は無次元であることを表しています。また、計算しやすくするため

として、圧力を再定義します。これは、今回の計算は非圧縮性を仮定しており密度ρが一定であるため、こういった操作が出来ます。

  • 質量保存の式

  • ナビエ・ストークス方程式

ここで『Re』はレイノルズ数と言います。イメージはRe数(レイノルズ数と読みます)が大きいほど勢いがある流れで、小さいほど粘っこい流れ(どろどろ・ねっちょりとかの方がイメージしやすいですかね)です。

本記事の計算では、無次元化された支配方程式を使っていきます!

計算のフローチャート

次は、計算のフローチャートを載せておきます。このフローチャートによって、計算の流れや、それぞれのルーチンで何をやっているのかといった理解が進むかと思います。ここで、フローチャート内の『NSE』は『Navier-Stokes Equations』の頭文字を取っています。

計算は『Fractional Step Method』という方法を採用してます。この手法の概要だけ説明すると、初めに支配方程式の圧力項を無視して計算し流速を求め、次に圧力方程式を解いて速度の修正をする、といった手法になります。

おそらく、「圧力項を無視?」「速度の修正?」といった疑問は出るかと思いますが、あえてここでは割愛します。そのため、今は「Fractional Step Methodという手法があり、それは圧力項を無視して仮の流速を求め、無視した圧力を考慮するために圧力方程式を解いて流速を修正している」といった認識をしておいて下さい。

ちなみに、この手法を採用したのは取り扱いが簡単だからです。

フローチャートのそれぞれについて軽く説明していきます。

 

1.格子生成

流れ場を計算するため、格子を作っていきます。イメージ的には配列に座標値を与える感じです。

よくわからない場合は、後述のコードを参考に『手を動かして、コードを実装して、可視化して、何をやっているか』を考えてください。可視化するとこんな感じです。

ちなみに、格子はレギュラー格子としています。これは格子の交点に変数を定義しているものです。

 

2.初期条件

流れ場全体にわたって主流方向速度に『1』という値を与えます。これは、流速が代表速度と一致していることを表しています。

 

3.境界条件

今回は差分法という手法を用いているため、境界の値が必要となります。境界というのは、今回の流れ場でいう流入・流速・壁面の事を指します。「なぜ差分法だと境界の値が必要となるの?」みたいな質問は置いといて、とりあえずコードの実装をしてみて下さい。

 

4.NSEの各項を計算

NSE(ナビエストークス方程式)の各項を計算します。とは言っても、やるのは微分の計算をするだけです。微分とはつまり差分で、そのために離散化して・・・みたいなごちゃごちゃしたことは無視です。

とりあえず、コードを実装してみて雰囲気つかんでください。(詳しいことは、後日まとめたいと思います。)

また、各項を計算と言いましたが圧力項はしなくて良いです。それがFractional Step Methodなので。

 

5.速度(仮)の計算

これは簡単です。加速度に時間を掛けたら速度が出ますよね(次元を見れば明らか)。そういう感じです。

ちなみに、ここでいう時間というのは時間刻みとか微小時間とか色々な名前で呼ばれるものです。要は、ちょっと時間が経ったときに流速はどうなるか、を計算して求めます。そして、そのちょっとを何千回と繰り返して流れ場がどうなるかを見ます。

この『ちょっと』の時間刻みをどのように決定すれば良いかは、制限があったりしますが割愛です。(クーラン数というのが関わってきます。)

 

6.圧力方程式を計算

ここのコード実装が一番の難所かもしれません。

行っている計算は圧力方程式・・・つまり、ポアソン方程式を解いています。ここで、ポアソン方程式というのは、以下の形をした方程式のことです。

  • ポアソン方程式

 

この計算で解くべきポアソン方程式も示しておきます。

  • 圧力方程式

 

この方程式を解いて圧力を求めます。

ちなみにポアソン方程式を数値的に解く方法は色々ありますが、今回はガウス・ザイデル法を用いています。これも詳しい説明は(以下略

 

7.速度の修正

圧力方程式を解いたら、『5.速度(仮)の計算』で求めた速度を修正します。

 

以上が一連の計算の概要となります。かなり詳細を省いているため疑問点も多いかと思いますが、一度コードの実装をして挙動を確認した後に、個別に詳細を調べるのが良いと思います。

それではコードの構築をしていきましょう!

コード構築

まずは、モジュールでの宣言を載せておきます。

こんな感じです!

それぞれの文字について軽く説明していきますね。

Lx x方向長さ
Ly y方向長さ
Nx x方向の格子数
Ny y方向の格子数
ltime 計算繰り返し数
dt 刻み時間
Re レイノルズ数
r1_dt 刻み時間の逆数
r1_Re レイノルズ数の逆数
x x方向座標
y y方向座標
dx x方向の格子幅
dy y方向の格子幅
U x方向流速
V y方向流速
P 圧力
dU 流速Uの微分
dV 流速Vの微分

 

1.格子生成

 

まずは、格子生成していきます

格子生成のコードはこんな感じになります。解説をしていきますが、コードを見て何となく分かりますかね。ここで『x』と『y』は格子の座標を、『dx』と『dy』は格子間の幅の大きさを表します。

さて、コードを作成したらそれが妥当か確認しましょう。下が出力用のコードです。

これを格子生成のプログラムの後に加えて出力してみましょう。出力方法についてもここで説明します。

 

まずは作成したプログラムを実行してください。するとフォルダ内に『12_xyGrid.txt』が作られます。

次に『Paraview』というソフトをインストールして下さい。インストール方法は割愛ですが、リンク先がダウンロードページとなっているので、そこからzipファイルをダウンロードして進めて下さい。

インストールしたら『Paraview』を起動してください。起動したら、『Paraview』の画面まで『12_xyGrid.txt』をドラッグ・ドロップして下さい。すると下のような画面になると思います。

 

ここで、データをインプットするために『Apply』(緑色の箇所)を押しますがその前に3点やってほしいことがあります。

  1. 左下の『Have Headers』のチェックを外す
  2. 『Field Delimiter Characters』にて空白(スペース)を入力
  3. 『Merge Consecutive Delimiters』にチェックを入れる

です。この3点の操作をしたあとに『Apply』を押してください。

 

次にデータを可視化させます。まずは下の画像を見てください。

ツールバーにある『Filters』⇒『Alphabetical』⇒『Table To Points』を選択して下さい。

 

選択すると下のように、『Table To Points1』が追加されると思います。

ここで、また『Apply』を選択する前に2点やってほしいことがあります。

  1. 『Y Column』を『Field 1』に変更
  2. 『2D Points』にチェック
  3. 『Merge Consecutive Delimiters』にチェックを入れる

です。そうすると下の画像のようになるので、右側の『SpreadSheetView1』の『×』を押して閉じて下さい。

 

そして、『Table To Points1』の左側(目のマーク)をクリックして下さい。すると計算領域が可視化されると思います。

 

良い感じですね!これで構築した格子生成コードの確認が出来ました!

2.初期条件

 

『x方向流速のU』と『y方向流速のV』、『圧力P』に初期条件を与えます。とは言っても簡単なのでさらっといきましょう。

まずはコード例を示します。

コードはこれだけです。『x方向流速U』だけ流れている状況になります。

さて、これを可視化してみましょう。手順は先ほどの『格子生成』と同じなので、『14_Ini_Condition.txt』をさっきと同じように格子の可視化まで行ってください。

格子の可視化まで出来たら、画像のように『Solid Color』の『Field 2』を選択して下さい。この『Field 2』が『流速U』に対応してるので、『流速U』の可視化になります。ちなみに、『Field 3』と『Field 4』は『V』と『P』に対応しています。これはテキストへの書き出す際の順番に対応しています。

ただ、これだと可視化しているか分からない場合もあるので可視化の範囲を変えます。画像のように『Rescale to Custom Data Range』を選択して下さい。

 

選択すると範囲を指定できるので、今回は最小値を0として最大値を1とします。

 

入力したら『Rescale』を選択して下さい。すると下の画像のようになります。

 

真っ赤になっていて、流速Uが値を持っていることが分かりますね!

3.境界条件

 

次は境界条件です。今回はチャネル流の計算なので、流入・流出・壁面条件を与えます。

以下にコード例を示します。

境界条件は、基本的に『値を与える(ディリクレ境界条件)』『勾配一定(ノイマン境界条件)』を与えます。

今回の計算では以下のような境界条件となります。

  • 流入でディリクレ境界条件(流速1)
  • 流出でノイマン境界条件(勾配一定)
  • 壁面はディリクレ境界条件(流速0)
  • 圧力はノイマン境界条件(勾配一定)

さて、可視化してみましょう。手順は今までと同じです。

これは『流速U』の可視化になります。少し見にくいですが、壁面で流速が0となっているのが分かりますね!

4.NSEの各項を計算

 

NSE(Navier-Stokes Equations)の各項を計算します。とは言っても、やることは微分の計算です。正確には一階と二階の2次精度中心差分です。差分式の導出は割愛しますが、これは今後記事にしたいと思います(多分)。差分式だけ載せておきます。

  • 一階の2次精度中心差分

  • 二階の2次精度中心差分

 

さて、いきなりNSEの各項を計算しても良いですが、その前に簡単な微分の計算をしてコードの妥当性を確かめたいと思います。

確かめる計算は、sinxの一階微分と二階微分です。

  • sinxの一階微分⇒cosx
  • sinxの二階微分⇒-sinx

となるため、比較は非常に簡単だと思います。

それでは、コード例を示します。

ここで、コードの妥当性を確かめるため、14行目にて『U』にsinxを与えてあります。

また、19~20行目の『Ux』と『Uxx』はそれぞれ『U』のx方向の一階微分と二階微分を表しています。

それでは実行して、『16_Ux_Uxx.txt』を作成して下さい。作成したら適当なソフトでグラフを作ってみてください。ちなみに、おすすめのソフトは『Sma4』です。ドラッグ・ドロップで簡単にグラフを作れるので良いですよ。愛用しています。(リンク先からダウンロードできます)

さて、一階微分のグラフは以下のようです。

二階微分もこんな感じ。

ここで、白丸が計算したデータで、実線がSma4内で作成した関数を描画したものです。どちらの微分も正しく計算できていますね!

それでは、この作成したコードを基にしてNSEの各項を計算していきます。コードを以下に示します。

少しだけ説明します。

22~23行目にある『dU、dV』の右辺第一項は粘性項を表し、第二項は対流項を表します。

『dU、dV』の式がなぜこのような形となるか(なぜ圧力項を無視しているのか)は、『Fractional Step Method』について調べると分かると思います。ここについても、そのうち記事にしたいなーとは思ってます。

後日記事に…が多いですが、あくまでも『即実装』が目的なので悪しからず。

5.速度(仮)を計算

 

ここは速攻で終わります。先ほど計算した『dU、dV』に時間刻み『dt』を掛けて、前の時間の流速に足してあげるだけです。

コードはこんな感じです。簡単!!

6.圧力方程式を計算

 

色々考えましたが、ここは非常に参考になるサイト(世界一易しいPoisson方程式シミュレーション)があるため、このサイトを基にコードの構築をしていきます。

以下にコード例を載せておきます。

基本的には世界一易しいPoisson方程式シミュレーションに沿ってコードを構築しているため、詳しい解説はそちらをご覧ください。

ただし、圧力方程式を解くようにコードを変えてあるので何をしているのか概要だけ説明します。

まず、このルーチンで解きたい方程式は以下のようになります。

この式の右辺は分かりやすいように、コード内では『f』という文字に値を入れています。コードの行数で言うと、19~20行目が右辺の計算をしています。

あとは繰り返し『圧力』と『誤差』を計算し、『誤差』が『収束判定をクリア』したらこの繰り返し計算を終了し、『圧力』を決定します。

その『収束判定』に対応するのがコードの31行目です。

本当は、ここでも妥当性の確認をするべきですが僕自身の知識不足もあり割愛します。そのため、計算を回してみて挙動が怪しかったらバグ取りに取り組むという形を取りたいと思います。

非常に申し訳なく、また自分自身も悔しいため更に勉強に励みます。

ただ結果論になりますが、このコードで計算を動かしてみたところそれらしい結果になったため、このコードで計算を回して下さい。

7.速度の修正

先ほど求めた圧力を用いて速度の修正を行います。コードは以下のようになります。

この計算コードは、先ほどの圧力方程式のコードに付け加えて下さい。

これにて、一通りの計算は終了です。お疲れさまでした。

後は各々のサブルーチンを組み合わせて計算をするだけです!!

計算開始

計算を動かす前に、上記のようにサブルーチンを組み合わせて計算できるようにしておいて下さい。計算の流れとしては、初めに示したフローチャート通りなのでちゃちゃっと実装できると思います。

また、計算終了時に結果を出力するために19~27行のプログラムを付け加えて下さい。

  • 『100_output.txt』はParaviewで流れ場を可視化をするためのファイル
  • 『101_output.txt』は下流で流れ場の流速分布をグラフ化するためのファイル

さて、それでは実行してみましょう!!!

計算結果

計算は上手く回りましたかね。今回はステップ数(ltime)を1000としているため、若干計算時間がかかるかと思いますが、それでもすぐだと思います。

それでは、構築した流れ場の『可視下図』『下流の流速分布』を見てみましょう!!

(ちなみにRe数は200で、格子数は(Nx,Ny)=(150,50)です。)

Paraviewでの可視化

綺麗に可視化できていますね!

この可視化を見ると、壁面付近は流速が遅く、中央部では流速が早くなっていることが確認できます。チャネル流っぽくなっておりコードは正しそうです!

なお、この可視化をするにあたって基本的な操作は今までと同じですが、より綺麗に表示させるために、ツールバーの『Filters』⇒『Alphabetical』⇒『Delaunay 2D』を追加しています。『Table To Points』を追加後に、『Delaunay 2D』を追加することで綺麗な可視化が出来ます。

それでは、次は下流における鉛直方向の流速分布を見てみましょう。

下流における流速分布

これも非常に綺麗に分布しています!

如何にもチャネル流って感じで気持ちが良いですね!!

まとめ

本記事ではチャネル流の数値計算について解説しました。どうですかね、無事に計算結果を可視化することが出来ましたかね。

計算できたなら、Re数を変えて挙動の違いを見てみたり、境界条件や初期条件を変えて流れ場の様子を確認してみて下さい。他にも、出力を毎ステップ(や100ステップ毎)にしたりと色々と遊べると思いますので、色々いじってみて下さい!

とりあえず、一通りコードの構築が出来たと思うので、今後はそれぞれのルーチンで何をやっているか詳しく勉強していくと理解が深まり、より楽しくなると思います。

更に知っておくべき内容として以下の項目が挙げられますので、参考にして下さい!

  1. 格子の種類(スタガード格子など)
  2. CFL条件
  3. Fractional Step Method
  4. 差分法
  5. 時間積分
  6. ポアソン方程式
  7. 境界条件(対流流出条件など)

長くなりましたが、今回の記事が誰かの参考になったら幸いです。(もし間違っている箇所があればこっそり教えて下さい)

 

Twitter(piy0_gadget)もやっているので是非フォローしてね!
疑問があればリプやDMで気軽にどうぞ!

https://twitter.com/piy0_gadget/status/1251079790548807684?s=20

-Fortran, 数値計算, 流体力学

© 2020 ぴよ工房 Powered by AFFINGER5