数値計算 流体力学

【数値計算】スタガード格子の差分法【サンプルコード付】

2020年6月4日

 

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

この記事では、

  • 格子のメリット・デメリット
  • スタガード格子の差分法
  • 計算コードの構築

について話していきます。

本記事の一番の目的はスタガード格子の差分の計算コード実装です。

 

ネットや本を見ても、スタガード格子とは何かといった説明はあっても、どうやって計算コードに組み込むかを説明してないことが多いと思います。

そのため、本記事ではスタガード格子をどのようにして実装するかに焦点を当てたいと思います。

実際にFortranでコードを書いて挙動を確認しているので、ぜひ参考にしてみて下さい!

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

 

格子の種類(レギュラー格子・スタガード格子・コロケート格子)

実は...格子の種類と言うと少し語弊があり、正確には格子の変数配置の種類です。今回は便宜上、格子の種類という題にしました。

 

それぞれの特徴を一言で言うと以下の通りです。

  1. レギュラー格子は実装が簡単だが計算が不安定
  2. スタガード格子は実装が難しいが計算が安定
  3. コロケート格子は一般座標系でよく用いられる

それでは、もう少し具体的にそれぞれの格子について説明していきます。

 

1.レギュラー格子

レギュラー格子はスカラー値(圧力・密度など)流速の位置を、格子が直交する点に定義します。

上の図を見ると直感的に分かると思います。

定義位置がすべて同じであるため、差分の扱いは非常に簡単になりますが、その一方で計算が不安定になりやすいといったデメリットがあります。

この計算の不安定さを解決するために、スタガード格子が生み出されました。

 

2.スタガード格子

スタガード格子はスカラー値(圧力・密度など)は格子の中心に、流速は速度成分が面に垂直に交わるように定義します。

これも上の図をみると直感的に分かると思いますが、れぞれの定義点が異なっているのでレギュラー格子と比べて複雑になっていますね。これが、コードの実装を難しくさせている要因となっています。

ただし、メリットも勿論あります。図のように隣り合う圧力の間に流速があるため、圧力差によって流れが駆動されるという、物理的に正確な定義の配置となっており、これがより自然な流れを表現し、計算が安定になります。

 

3.コロケート格子

コロケート格子はスカラー値(圧力・密度など)流速の位置を、格子の中心に定義します。

これも上の図を見ると、直感的に分かると思います。もっと正確には流束成分を面に垂直になるように定義しますが、これがコロケート格子が一般座標系においてよく用いられている理由です。流束成分は一般座標系の方程式にあるからです。

ただ、コロケート格子は保存性が成り立たないといった指摘があり、修正コロケート格子が開発されています。(詳しくはググってください)

スタガード格子の差分

スタガード格子の差分方法について解説していきます。

図のように、格子(添え字が[i,j])とその格子の添え字に対応する流速成分を定義します。

まず、[i,j+1/2]におけるx方向の差分を取ってみますが、注意してほしいのが、流速uは差分をとる点([i,j+1/2])と定義点が一致していますが、流速vは差分をとる点([i,j+1/2])と定義点が一致していない事です。

この注意点を意識して読み進めて下さい。

 

1.流速uの差分

ここで、u[i+1/2,j+1/2] ,u[i+1/2,j-1/2]はともに定義点にない値ですよね。そのため(x方向の)隣り合う点から値を補間する必要があります。具体的には以下の通りです。

それではこの値を差分式に代入してみましょう。すると以下のようになります。

 

これで、流速uのx方向差分が出来ました。

 

 

2.流速vの差分

ここで、v[i+1/2,j+1/2] ,v[i+1/2,j-1/2]はともに定義点にない値ですよね。そのため(y方向の)隣り合う点から値を補間する必要があります。具体的には以下の通りです。

この値を差分式に代入してみます。

 

これで、流速vのx方向差分が出来ました。

 

スタガード格子のコード実装

さて、スタガード格子の実装です。計算コードは、フローチャートに沿って構築していき、実際にsin関数を微分して挙動の確認をします。

 

モジュール宣言

まずはモジュールの宣言を載せておきます。モジュール???って人は、過去に記事を書いているので参考にしてみて下さい。

【Fortran応用編】サブルーチンとモジュール

こんにちは、ぴよ工房を運営しているぴよ(@piy0_gadget)です! この記事では、 サブルーチン(Subroutine)の使い方 モジュール(Module)の使い方 について話していきます。 サ ...

続きを見る

 

 

1.格子生成

差分を取るためには座標の値が必要になるので、まずは格子を作っていきます。

 

2.初期条件

流速u,vに初期条件を与えます。今回は差分の結果を見たいのでsin関数を与えます。

 

3.差分の計算

準備が整ったので差分を取っていきます。とは言っても難しくはなく、前述した通りにコードを実装するだけです。

今回用いる差分式をもう一度載せておきますね。

 

また、添え字が分数だと扱いにくいので以下のように再定義します。

 

これをコードに落とし込むと以下のようになります。

 

4.出力

計算した結果を出力して、正しく差分が出来ているか確認しましょう。今回、初期条件で与えた関数がsinなので、差分した値はcosの関数になっているはずですよね。

以下は出力プログラムです。

 

さて、これをグラフにして本当にcos関数になっているか比較してみます。

 

〇が計算した値で、実線が理論解となっています。

まとめ

本記事では、スタガード格子の差分式の導出、計算コードの構築』の解説をしました。

格子の説明をしているサイトはあっても実際にコード構築までしているサイトは珍しいと思うので、ぜひ参考にしてもらえると幸いです。

またスタガード格子の差分を用いて、前回のチャネル流の記事のコードを書き換えるのも面白いと思いますので、ぜひ取り組んでみて下さい。

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

  こんにちは、ぴよ工房を運営しているぴよ(@piy0_gadget)です! この記事では、 二次元チャネル流のCFDコード構築 の解説をしていきます。 本記事の一番の目的は『即実装』です! ...

続きを見る

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

 

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

-数値計算, 流体力学

© 2020 ぴよ工房 Powered by AFFINGER5