OpenFOAM

【OpenFOAM】移流方程式で色々な離散化スキームをPythonで自動実行して試してみた

こんにちは(@t_kun_kamakiri

今回はOpenFOAMを使って移流方程式を解きたいと思います。その際に離散化スキームをどのように選択すると解の振る舞いが変わってくるのかを見ていきます。

前回の記事では移流方程式をの移流項を中心差分と風上差分で計算した結果を見ました。
中心差分は解が振動してぐちゃぐちゃになり、風上差分は解が拡散しながら右に進行する解の振る舞いをしました。

色々なスキームを試しますが理論的な解説はここではしません。「スキームを変えると解の振る舞いがどう変わるのか」を知ってもらえればと思います。
勉強会でも「スキームがわからない」という質問もたまに聞くため、ナビエストークス方程式のような難しい方程式からせず、まずは今回のような簡単な移流方程式から体験すれば良いのではないかと思います。
スキームを手で変えていくのは面倒なのでPythonを使って自動実行しますので、その辺も参考になればと思います。

本記事の内容
  • OpenFOAMの移流方程式を使って色々なスキームでの解の振る舞いを見る。
  • 計算の実行はPythonで自動化する

以前にも移流方程式を差分法を使ってFortranで解いたことがありますが、移流方程式は式が簡単で解析解があるにも関わらず数値的に解こうとすると、不安定だったり拡散してしまったり難しい問題があります。

1次元の移流方程式
\begin{align*}\frac{\partial T(x,t)}{\partial t}+u\frac{\partial T(x,t)}{\partial x}=0\tag{1}\end{align*}
※\(u\):一定速度
1次元の移流方程式は、下記の絵のように何かの物理量を形を変えずに一定速度\(u\)で運ぶ方程式を意味しています。
横軸は空間\(x\)ですが、縦軸は「密度」や「温度」だったり色々考えることができます。本記事では、温度\(T\)として扱うことにします。
Fortranを使った差分法の結果は以下のようになりました。
この結果を見ると中心差分は振動してしまい、前進差分と後退差分は流れの向きによっては安定だったり不安定だったりします。上の例のように右側に進行する($u>0$)場合は、安定に進みますが解が拡散してしまっています。

このように、流体の第二項の対流項$u\frac{\partial T}{\partial x}$の微分の離散化方法によっては、数値解が異なることがあります。異なると言っても、通常は流体のナビエストークスには流体粘性ので不連続解を分散させて解が振動しないようにしますし、不連続な解が出ることも稀なのでどのような離散化手法をとってもあまり解が変わりません。

ただ、衝撃波が出ると速度や圧力の不連続解が生じ、さらに粘性項より対流項が大きくなるので上記に近い方程式になります。

離散化スキームを以下のようにざっくり分けて特徴をまとめておきます。

離散化スキームざっくりまとめ

upwind_difference.pdf (xrea.com)

  • 1次精度:風上差分(upwind)→こちらは安定だが解が拡散してしまう
  • 2次精度:中心差分、QUICK→こちらは不連続解をシャープにとらえようとするが、数値振動が起こりやすい
  • TVDスキーム:minmod、vanLeer、vanAlbada、superBeeが有名。
    色々な種類がありますが、不連続部分で振動を抑えるようにしながら、解はシャープにとらえようとするいい感じのスキーム

離散化スキームの種類はむちゃくちゃありますので、まずはこれくらいざっくりわけて特徴を理解しておくとよいでしょう。

こんな方が対象
  • OpenFOAMの移流方程式を使って離散化スキームの基礎的な内容を理解したい方
  • OpenFOAMのチュートリアルから独自に検証するケースファイルを作りたい方

離散化スキームについて勉強したい方の参考になれば幸いです。

OpenFOAM-v2012(WSL2)
Python3.8.10(jupyter lab)

モデルの入手

移流方程式のモデルは前回の記事で作成しましたのでそちらをご参考ください。
サンプルモデルも用意したのでご利用ください。

フォルダ構成

フォルダ構成は以下のようになっています。

  • orgCase:元のケースファイル
  • schemeTest:スキームの一覧
  • Allrun:計算実行スクリプト
  • main.py:計算の自動実行スクリプト

ひとつずつ解説していきます。

orgCase:元のケースファイル

こちらはスキームを変えて実行するための元になるケースファイルです。
これをコピーしてスキームを変えていくというのを行います。

scalarTransportFoamは以下の偏微分方程式を解く非定常のソルバーです。

\begin{align*}
\frac{\partial T}{\partial t}+\nabla\cdot(\boldsymbol{u}T)-\nabla\cdot (D_{T}\nabla T)=\boldsymbol{S}_{T}\tag{1}
\end{align*}

本記事では(1)を少し簡単にするために、

  • 生成項を0:$\boldsymbol{S}_{T}=0$
  • 拡散項を0:$D_{T}=0$
  • 流速は一定:$\boldsymbol{u}=(0.1, 0, 0)$

と設定するようにします。

初期状態は以下のようになっています。

schemeTest:スキームの一覧

schemeTestというファイルにスキームの一覧を記載します。

schemeTest

これをケースファイルの「system/fvSchemes」の「divSchemes」を変えていきます。

system/fvSchemes

Allrun:計算実行スクリプト

計算を実行するためのスクリプトを作成します。

Allrun

OpenFOAMが用意しているスクリプトの「RunFunctions」には便利な関数が用意されているので利用すると良いでしょう。
例えば、「runApplication $(getApplication)」と書くと計算させたソルバを実行しながらlogファイルも作成してくれます。

main.py:計算の自動実行スクリプト

計算実行のパラメータスタディを行う前にまずはやりたいことを整理した方が良いでしょう。

構成を考える

  1. スキームの一覧ファイルを作成
  2. baseをコピーする
  3. スキームを入れ替える
  4. Allrunスクリプトを実行する(postProcess)
  5. グラフにまとめる
  6. 全部のグラフを2列ずつにしてまとめる

jupyter labでデバッグしながら数でまとめられるものはまとめてしまい、mainのプログラの記述を少なくするようにします。

結果はスキームごとに分けて「resultDir/savefig」に画像を作成し保存、全てのスキームでの結果を「graph_AllShemes.png」に画像を作成し保存しています。

pythonを実行します。

以下のように計算が進みます。

※WSLでplt.figure()などグラフ表示でエラーが出る場合は「export DISPLAY:=1.0」とするとエラーがなくなります。

結果を確認

結果は「resultDir/savefig」に画像として保存するようにしています。

例えばスキーム「linear(中心差分)」の結果は以下のようになります。
ぐちゃぐちゃですね。

例えばスキーム「upwind(風上差分)」の結果は以下のようになります。
安定していますが、拡散しています。

全てのスキームをまとめたものが「graph_AllShemes.png」です。

スキームは有名で実装が容易なものから聞いたことないものまで色々あります。
例えば以下のツイートのようなスキームは勉強不足で初めて聞きました💦

まとめ

スキームは試して結果を眺めてみないとわからないので、本記事のような簡単な移流方程式で色々と試してみました。
計算時間まではまとることをしていませんが、スキームによって計算時間に差が出るみたいなのでそちらも比較しても良いかもしれません。

Pythonの自動スクリプトなども参考になれば幸いです。

参考書

有限体積法の勉強にはこの本で間違いないです。
対流項目が中心差分で安定しないことを理論的に評価しておりとても参考になります。

数値流体力学 [第2版]

数値流体力学 [第2版]

H.K.Versteeg, W.Malalasekera
10,450円(02/29 01:32時点)
Amazonの情報を掲載しています

有限体積法の詳しい解説として英語ですが以下の書籍を挙げておきます。
OpenFOAMに沿った解説もされているのでOpenFOAMを勉強する人にとっては有用です。

The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM® an...

The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM® an…

Moukalled, F., Mangani, L., Darwish, M.
10,295円(02/28 22:53時点)
発売日: 2015/08/13
Amazonの情報を掲載しています

OpenFOAMの基本的な設定を学ぶことができる書籍として以下の参考書がお勧めです。

改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ(NextPublishing))

改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ(NextPublishing))

川畑 真一
1,320円(02/28 07:45時点)
発売日: 2022/04/15
Amazonの情報を掲載しています
OpenFOAMによる熱移動と流れの数値解析(第2版)

OpenFOAMによる熱移動と流れの数値解析(第2版)

3,520円(02/29 00:46時点)
Amazonの情報を掲載しています

本書はESI版ではなくFoundation版でありOpenFOAMv2.xで書かれたやや古い内容です。その点に注意してESI版のお使いのversionとの違いを丁寧に調べることができる方であれば、十分参考になる書籍です。

OpenFOAMプログラミング

OpenFOAMプログラミング

Mari´c,Tomislav, H¨opken,Jens, Mooney,Kyle
8,250円(02/28 11:51時点)
Amazonの情報を掲載しています

OpenFOAM全般の設定などまとめられた書籍として以下のものがあります。
こちらもFoundation版なので、ESI版のお使いのversionとの違いを丁寧に調べることができる方であれば、十分参考になる書籍です。

OpenFOAMライブラリリファレンス

OpenFOAMライブラリリファレンス

株式会社テラバイト, 人見 大輔
16,500円(02/28 16:48時点)
Amazonの情報を掲載しています

こちらのオープンCAE学会から技術書典にOpenFOAMでメッシュ作成【PDF版】が出ています。難しい内容ではありますが、OpenFOAMのメッシュに関する内容が1000円で学べます。

関連記事もどうぞ

COMMENT

メールアドレスが公開されることはありません。 が付いている欄は必須項目です