OpenFOAM

【OpenFOAMカスタマイズ】温度依存する温度拡散係数に変更

こんにちは。
本記事ではOpenFOAMのカスタマイズのためのメモを残しておきます。

今回は熱拡散方程式(熱伝導方程式とも呼ぶ)の以下の式を扱います。

\begin{align*}
\frac{\partial T}{\partial t} = \nabla \cdot (D\nabla T)\tag{1}
\end{align*}

OpenFOAMのソースコードでも以下のようになっており上記のイメージと結びつきやすいでしょう。

今回は熱拡散方程式の温度拡散係数に温度依存性をさせたものにします。

今回の内容

温度拡散係数の温度依存性を追加

\begin{align*}
DT=A_{1}T^2+A_{2} T +A_{3}\tag{2}
\end{align*}

本記事は以下の記事を参考にしています。
以下の記事の内容は、数年前に書かれているため年月が経つと今のOpenFOAMのバージョンでのコードと変わってきています。

今のバージョンOpenFOAM-v2112で解説をし直したのがこの記事です。

本記事がOpenFOAMのソースコードの理解の助けになれば幸いです_(._.)_

ラプラシアンソルバの編集

まずはlaplacianFoamソルバのベースとなるコードを「laplacianFoamDepT」という名前のフォルダにしてコピーします。

「$FOAM_APP」は環境変数で「/opt/OpenFOAM/OpenFOAM-v2012/applications」とパスの設定がされているので、各自の環境に合わせて読みかえてください。

※カスタマイズするものによっては「$WM_PROJECT_USER_DIR = /home/kamakiri/OpenFOAM/kamakiri-v2012」などの環境変数に応じたフォルダでカスタマイズソースを作成する方がOpenFOAMのファイル構成で扱いやすいかもしれません。今回はその辺は気にしなくても良いです。

以下のコマンドで「laplacianFoamDepT」フォルダに移動します。

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

以下のコマンドでファイル名を変えておきます。

温度拡散係数の温度依存の追加

温度拡散係数はデフォルトでは温度依存の定数ですが、温度依存する物性値というのはおおくあります。
今回は温度拡散係数の温度依存を以下のように2次の項までの展開した形に変更を行います。

\begin{align*}
DT=A_{1}T^2+A_{2} T +A_{3}\tag{2}
\end{align*}

「createFields.H」の「DT」を以下のように変更します。

createFields.H

まずはtransportPropertiesという名前でIOdictionary型を用意します。

「runTime.constant()」は文字列のconstantを返します。つまり、OpenFOAMのケースのconstantフォルダに当たります。constantの中の「transportProperties」ファイルを読み込むようにしており、「IOobject::MUST_READ」でファイルの読み込みを、「IOobject::NO_WRITE」ファイルの書き込みを行わないように設定しています。

つぎに、温度拡散係数の各項の係数を定義します。

A1, A2, A3はスカラーを持った量として定義する必要があります

  • $A_{1}T^2$
  • $A_{2} T$
  • $A_{3}$

の各項は温度拡散係数と同じ次元$m^2/s$を持つようにするため、

  • $A_{1}$:$m^2/K^2 s$
  • $A_{2}$:$m^2/Ks$
  • $A_{3}$:$m^2/s$

の単位を持つように設定する必要があります。
ここで単純に値を設定して「1.0T*T + 2.0*T +3.0」などとしても上手くいきそうな感じはしますが、実際は異なる物理量の足し算や引き算は間違いのモノですのでOpenFOAMでは単位がそろったもの同士の足し算引き算を行うようになっています。
この単位の設定は後ほどケースフォルダの「constant/transportProperties」で行います。

「volScalarField」という型で変数名DTとして定義します。

ここでは初期値を設定しているだけです。
createFields.HはlaplacianFoamDep.Cで時間発展のループの前にインクルードで読み込んでいるため、 A1*T*T+A2*T+A3と設定しても初期値の温度でDTの値を計算しているにすぎません。

ちなみに、GeometricFieldをtypedefでvolScalarFieldに定義しなおしたものなので詳細を知りたい場合はGeometricFieldを見ることになります。
引数的に以下のコンストラクタが呼び出されているのかと思います。

IOobjectの中で「IOobject::MUST_READ」となっているため、“DT”という名前で「runTime.timeName()←これは出力した時間のフォルダ名」から値を読み込みみます。
そして「 IOobject::AUTO_WRITE」となっているので、controlDict内で設定した出力時間に応じてDTの値を出力するようにしています。
では、温度依存の計算を行うコードをlaplacianFoamDep.Cに追加します。
laplacianFoamDep.C

変更したのは以下の部分です。

では、できたら実行コマンドを作ります。

Make/files

Make/options

以下のコマンドでコンパイルを行います。

「$FOAM_USER_APPBIN」にlaplacianFoamDepができていたら成功です。

チュートリアルの編集

ケースファイルは前回の記事同様です。
以下に設定ファイルを載せておきます。

ラプラシアンソルバを使った適当なチュートリアルをコピーしてきます。

このチュートリアルはblockMeshが無いので適当に持ってきます。

blockMeshDictは適当に設定します。

system/blockMeshDict

では、blockMeshを実行します。

ParaViewでモデルを読み込んで確認しましょう。

今回はわかりやすく1×1×1(m3)でx方向に5分割したモデルにしました。

初期条件として「0」フォルダにはTを用意します。

0/T

前回の記事でDT設定ファイルが必要でしたがcreateField.HでDTの定義は以下のようになっているため、0/DTは必要ありません。

注目は「IOobject::NO_READ」の部分です。
「runTime.timeName()」は毎時間出力したフォルダを指しますが、「IOobject::NO_READ」となっていることで、そこにDTが無くてもエラーとなりません。「A1*T*T+A2*T+A3」で次元を持った初期値を設定しているイメージです。

温度拡散係数の係数をconstant/transportPropertiesで行います。

constant/transportProperties

DTは値だけですが、A1,A2,A3は次元と値を設定しています。
少し気持ちの悪い感じがしますので、A1,A2,A3も値だけを設定すれば計算ができるように後ほど変更します。実際は、A1,A2,A3の次元を頭に入れて計算するなんてユーザー目線からすると面倒なのでプログラムの中で自動的に次元も考慮するようにしたいものですよね。
後ほど変更するとして、今はこのまま進めていきます。

またプログラムが正常に動作しているかを確認するために、A1,A2,A3の値はわかりやすい数値にしました。

以下設定をします。

system/controlDict

system/fvSchemes

system/fvSolution

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

ParaViewで結果を確認しましょう。
特に温度分布から温度拡散係数を手計算で求めた値と一致しているかを確認します。

問題はなさそうです。

コードを編集する

createFields.Hが以下のように書いていたせいでconstant/transportPropertiesにA1, A2, A3の次元を書く必要がありました。

createFields.H(1部抜粋)

constant/transportProperties

これではDTは値だけ書けばよく、次元は意識しなくても良いのにA1,A2,A3は次元と値も書く必要があり少々面倒に感じます。
できれば次元を意識することなくA1,A2,A3の値のみを設定したいです。
そのためには、以下のようにcreateFields.HでA1,A2,A3の次元も一緒に定義しておけば良いです。

createFields.H(1部抜粋)


※dimensionSet(kg, m, s, K, mol, A, Cd)のように関数の引数に次数を記述します。
※ちなみにdimensionedScalarは以下のようにdimensioned型にテンプレートでscalar(実数)を定義しなおしています。

以下のようにコンストラクタで引数を設定します。


変更したら再度wmakeでコンパイルしなおします。

そうするとconstant/transportProperties以下のようにシンプルな記述で済みます。

constant/transportProperties

以上のように設定できたら再度計算をして結果を確認してみましょう。
先ほどと同じ結果だったらOKです。

補足(コンパイル時のwarning)

※コンパイルすると以下のようなwarningが出ます。

どうやら

のような書き方は非推奨となっているようです。

以下のような書き方で良いようです。

以上のように書き換えて再度コンパイルするとwarningは消えます。
解析結果も先ほどの結果と同じです。

まとめ

温度拡散係数に温度依存性を持たせた熱拡散方程式を解くようにカスタマイズしました。

特に温度依存の係数は適当に設定したので意味のないカスタマイズかもしれませんが、OpenFOAMのカスタマイズの練習としては良い例題かなと思います。

参考書

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

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

3,520円(01/31 22:36時点)
Amazonの情報を掲載しています
OpenFOAMプログラミング

OpenFOAMプログラミング

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

Foundation版で書かれていますが、OpenFOAMの型定義や関数などを調べるのによくまとまった辞書です。

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

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

株式会社テラバイト, 人見 大輔
16,500円(01/31 16:27時点)
Amazonの情報を掲載しています
関連記事もどうぞ

COMMENT

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