OpenFOAM

【OpenFOAM】移流方程式で色々な離散化スキーム(解の振る舞いと計算時間)

こんにちは(@t_kun_kamakiri

前回の記事では移流方程式における移流項の離散化スキームを変えて結果がどう変わるかを検討しましたが、今回もその続きです。

前回の記事でも離散化スキームについて書きましたのでご参考ください。

今回は前回の記事以上に色々なスキームを試してみました。OpenFOAMに用意されているスキームは色々あるので自分で実装するよりも簡単にスキームの検討ができて便利ですよね。

勉強会などで「スキームなどを勉強したいけど、どうしたら良いのか」という質問もたまにあります。そんなときはナビエストークス方程式のような複雑な方程式ではなく、比較的簡単で理論解もある移流方程式で色々なスキームを試してみるのが良いでしょう。

本記事では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*}

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

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

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

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

upwind_difference.pdf (xrea.com)

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

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

本記事の内容
  • OpenFOAMでの移流項での離散化スキームについて結果比較
  • Pythonで計算を自動化

※離散化の理論的な解説は割愛します。

本記事は前回の記事とは以下の点で異なります。

  • スカラー輸送方程式(scalarTransportFoam)に計算時間を追加
  • 前回の記事よりも多く色々なスキームを試した
  • グラフに計算時間を載せるようにした

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

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

モデルの入手

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

モデルを入手して本記事の内容に沿って編集していくと再現できるでしょう。

フォルダ構成

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

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

スカラー輸送方程式(scalarTransportFoam)に計算時間を追加

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)$

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

  • $\Delta x = \frac{1}{100} = 0.01$[m]
  • $\Delta t = 0.01$[s]

なのでクーラン数$C=\frac{u\Delta t}{\Delta x }=0.1$としています。
クーラン数1だとSuperBeeでエラー出たので少し小さめにしておきました。

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

理論解は以下のように形を変えずに速度に乗って進行します。

しかし、微分方程式を離散化して代数的に解くためどうしても近似的な解になってしまいます。どのような離散化スキームを使うとどういう結果になるのかを知ることが本記事の目的であります。

もう一つの目的は、離散化スキームによって計算時間が異なってくるため計算時間を調査することです。しかし、困ったことにscalarTransportFoamはデフォルトで計算時間を出力する関数が使われていないため、計算時間がわかりませんでした。

そこで計算時間を知るために以下の方法を考えました。

  1. system/controlDictにcodeを埋め込む
    →上手くいかなかったので断念(自分がわかっていないだけなのでわかる方いたらコメント頂けると助かります_(._.)_)
    追記:「time -f ‘,%U’ scalarTransportFoam > log.scalarTransportFoam」でコマンドを実行するのに使用したユーザーCPU時間が計測できるようです。
    なので②のようにCPU計測のための関数を追加しなくても良かったです…
  2. カスタマイズする
    関数「runTime.printExecutionTime(Info); //Add」を追加するだけなのでこちらを採用します。

②の方法を採用しソルバにちょこっと関数を追加することで計算時間を出力するようにしました。

では、以下のコマンドでソルバをコピーします。

今回は使い捨てコードでローカルにコピーしても動作するコードなのでローカルにコピーして編集します。

フォルダ移動します。

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

ファイルの名前を変えます。

MyscalarTransportFoam.Cの時間のループの中に「runTime.printExecutionTime(Info); //Add」を追加して計算時間を出力するようにします。

MyscalarTransportFoam.C

そして実行ファイルを出力する先を指定してコマンドに名前を付けます。

Make/files

以下のコマンドコンパイルします。

これで「\$FOAM_USER_APPBIN」にMyscalarTransportFoamという実行ファイルが生成されています。ちなみに環境変数は「\$FOAM_USER_APPBIN=/home/ユーザ名/OpenFOAM/kamakiri-v2012/platforms/linux64Gcc63DPInt32Opt/bin」となっています。

では、計算時間を出力してくれるかを確認しましょう。
フォルダを移動して、

以下のコマンドで計算を実行します。

計算のログに「ExecutionTime = 0.05 s 」などが出ていればオッケーです。

計算時間の桁数はもう少し多くほしいですが・・・・


複数回計算させて平均を取るかもう少し大きなモデルで解析をするかなど考えないといけませんね。

ひとまず正常に動作しているようであれば、初期状態に戻しておきましょう。
以下のコマンドでメッシュ状態はそのままで計算した結果だけを削除することができます。

では、Pythonコードのあるフォルダに戻ります。

前回の記事よりも多く色々なスキームを試した

schemeTestというファイルにスキームの一覧を記載します。
前回の記事よりスキームは多くなっています。

schemeTest

グラフに計算時間を載せるようにした

今回は計算時間も測定してグラフに載せるようにしました。
計算時間はあくまで1回の計算をさせた場合での結果です。他のアプリの使用状態によってはCPUの使用状態も変わってくるので、本来は複数回計算させた際の平均値を取るのが良いでしょう。

Pythonを実行します。

Pythonで「system/fvSchemes」の「divSchemes」スキームを変えています。

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

pythonを実行します。

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

結果を確認

結果は「resultDir/savefig」に画像として保存するようにしていますのでじっくり眺めることができます。
全体のグラフの比較は「graph_Allschemes.png」で行うことができます。

graph_Allschemes.png

やはり複数回計算させた結果を計算時間とするべきでしょうね。なので、今回の計算時間は参考程度ということです。

 アニメーション作成

アニメーションにした方がわかりやすいかと思ったのでアニメーションにしてみました。

移流方程式というナビエストークス方程式よりは単純な形をした一般解のある偏微分方程式ですが、移流項$u\frac{\partial T}{\partial x}$の部分のスキーム違いで結果が大きく変わりますね。
特に初期状態の矩形波の変化が激しい部分の微分を離散化が難しく、離散化スキームをどれにするかによって結果が変わるということがわかりました。これがもしsin波のような滑らかな初期状態だったり、拡散項があるとあまり結果に違いがでないかもしれません….(試す必要がありますね)

アニメーション作成用のコードは即興で作ったのでお粗末ですがご参考に…
先ほどのPythonコードに以下の関数を追加して実行すればOK。

まとめ

スキームは試して結果を眺めてみないとわからないので、本記事のような簡単な移流方程式で色々と試してみました。
今回は結果と計算時間だけを一覧で眺めれるようにしました。スキームによっては解が安定しなかったり、安定していても拡散する傾向にあったりします。離散化スキームの理論的な話は割愛しましたが、勉強するとなぜそのようなことが起こるのかがわかってきて理解が深まります。

参考書

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

数値流体力学 [第2版]

数値流体力学 [第2版]

H.K.Versteeg, W.Malalasekera
10,450円(11/21 01:22時点)
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円(11/21 12:23時点)
発売日: 2015/08/13
Amazonの情報を掲載しています

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

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

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

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

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

3,520円(11/21 07:56時点)
Amazonの情報を掲載しています

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

OpenFOAMプログラミング

OpenFOAMプログラミング

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

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

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

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

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

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

関連記事もどうぞ

COMMENT

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