OpenFOAM

【OpenFOAM】移流方程式で離散化スキームの勉強をする

こんにちは(@t_kun_kamakiri

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

本記事の内容

OpenFOAMの移流方程式を使って中心差分と風上差分の結果比較を行います。

以前にも移流方程式を差分法を使って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$)場合は、安定に進みますが解が拡散してしまっています。
こんな方が対象
  • OpenFOAMの移流方程式を使って離散化スキームの基礎的な内容を理解したい方
  • OpenFOAMのチュートリアルから独自に検証するケースファイルを作りたい方

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

チュートリアルをコピーする

スカラー輸送方程式のチュートリアルをコピーします。

「$FOAM_TUTORIALS」は環境変数で「/opt/OpenFOAM/OpenFOAM-v2012/tutorials」と割り当てられています(個人の環境に依る)。

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

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

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

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

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

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

では、フォルダを移動します。

メッシュ生成

メッシュ生成はblockMeshを使います。

「system/blockMeshDict」ファイルを以下のように編集します。

system/blockMeshDict

以下のコマンドを実行してメッシュを作ります。

このように1.0×0.1×0.1の計算領域を作りました。

境界条件の設定

今回はxの正の方向へ流れが一様に0.1m/sで流れるようにします。

UT
XMin
type fixedValue;
value uniform (0.1 0 0);
or
type zeroGradient;
type zeroGradient;
XMax
type zeroGradient;
type zeroGradient;
YMin
type slip;
type slip;
YMax
type slip;
type slip;
ZMin
type empty;
type empty;
ZMax
type empty;
type empty;

Uがデフォルトでは既に速度の分布があるファイルになっているので、以下のように書き換えます。

0/U

流速は一定:$\boldsymbol{u}=(0.1, 0, 0)$です。
これは時間が進行しても常に同じ値です。

次に温度境界条件を設定します。

0/T

初期状態の設定

「setFields」で温度分布の初期条件を作ります。

このように$x$方向0.1~0.2だけ$T=1$とし、その他を$T=0$となるように設定します。

設定ファイルはsetFieldsDictに書きますが、デフォルトではファイルが無いので適当なチュートリアルから持ってくる必要があります。
有名なのはinterFoamのdamBreakにsetFieldsDictがありますが、それを知らない場合は以下のコマンドでsetFieldsDictが使われているチュートリアルを探すことができます。

コマンドを実行すると以下のようにずら~っと候補が出てきます。

適当に以下からコピーをしてきます。

以下のように編集します。

system/setFieldsDict

デフォルトは「volScalarFieldValue T 0.0」とし、「box ( 0.1 -1 -1 ) ( 0.2 1 1 );」の領域だけを「volScalarFieldValue T 1」という値にします。

では、以下のコマンドで初期状態を設定します。

物性値の設定

今回は「拡散項を0:$D_{T}=0$」とします。物性値の設定は「constant/transportProperties」で行います。

constant/transportProperties

 

これにより、(2)のスカラー輸送方程式は移流方程式になります。

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\)で運ぶ方程式を意味しています。
数値計算の結果が正しく解けていれば、上のように形を変えずに右側へ温度分布が進行していくはずです。

離散化スキームの設定

次に、離散化スキームを設定します。離散化スキームは「system/fvSchemes」で行います。
今回は移流方程式の$u\frac{\partial T(x,t)}{\partial x}$を移流項といいます。
この移流項の離散化方法によって大きく結果が変わることを確認します。
まずは、linear(中心差分)から見ていきます。

system/fvSchemes

その他のスキームは別の記事で詳しく確認したいと思います。

計算制御の設定

解析時間や時間刻みは「system/controlDict」で設定を行います。

  • 計算時間:5(秒)まで
  • 時間刻み:0.001(秒)

system/controlDict

 

行列ソルバの設定

今回は特に設定を変えずにいきます。

system/fvSolution

解析実行(中心差分)

では計算を実行します。

計算が終わったらParaViewで確認しましょう。
中心差分だと温度分布の形を保ったまま右へ進行するどころか、振動してとても不安定な結果になってしまいました。

これはこちらの記事の内容と同じ結果になっていますね。

OpenFOAMは有限体積法でFortranでの結果は差分法ですが、移流項を中心差分で離散化するとどちらも振動する結果となっています。

ナビエストークス方程式には拡散項などがあるので振動も抑えられてしまいますが、発散項の中心差分での離散化は単独だと振動する解が出てきます。移流項を中心差分で解くと空間離散化の2次精度とは言え、振動する解となってしまうことは十分注意すべき点です。

解析実行(風上差分)

linear(中心差分で)は振動する解が得られたので今度は安定な差分法である風上差分を使います。以下「system/fvSchemes」でコメントアウトして風上差分を設定します。

system/fvSchemes

Twitterにアップしたような拡散しながら右へ進行する解となりました。
こちらも記事の内容と同じ結果になっていますね。

まとめ

今回はスカラー輸送方程式のチュートリアルを変更して移流方程式にして、離散化スキームによる結果を見ました。

微分を代数的に解くために離散化しているため、完ぺきな解法というものはなく一長一短ということですね。今回は中心差分と風上差分のみ紹介しましたが、その他に「線形風上差分」、「QUICK」、「TVD」など色々とあります。このようにスキームに対してどのような結果になるのかを知るためには離散化手法を知ることが大事ですが、何より計算で確かめてみることが良いでしょう。

有限体積法の「中心差分」による離散化はこちらにまとめています。1次元の熱拡散方程式ですが、参考になるでしょう。

次回はスキームを色々変えて結果の違いを確かめたいと思います。複数ケースの計算を行うためスキームを変えて計算するプログラムを作って計算します。

対流項の発散スキームの比較」にシェルスクリプトで自動で計算するプログラムがあります。シェルスクリプトはちょっと慣れないのでPythonで行いたいと思います。

参考書

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

数値流体力学 [第2版]

数値流体力学 [第2版]

H.K.Versteeg, W.Malalasekera
10,450円(03/30 22:20時点)
Amazonの情報を掲載しています

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

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

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

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

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

3,520円(03/30 21:05時点)
Amazonの情報を掲載しています

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

OpenFOAMプログラミング

OpenFOAMプログラミング

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

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

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

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

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

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

COMMENT

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