OpenFOAM

【OpenFOAMカスタマイズ】icoFoamに温度場輸送方程式を追加。

こんにちは(@t_kun_kamakiri)。

前回の記事でも書いた通り、OpenFOAMのカスタマイズできるように勉強をはじめています。

前回の記事はこちら

今はまだネット上にあるOpenFOAMのカスタマイズ記事を自分で理解しながら手を動かして進めているだけなので、記事の内容的に目新しいものはありませんのでご了承ください。

本記事の内容は、こちらの記事を参考にして自分でもまとめなおしています。
自分で理解して記事にしたりとアウトプットをすると理解が深まります。

それに、まとめなおしていってもOpenFOAMのバージョンや仕様が違うとソースコードの中身が違うので注意が必要ですね。
参考記事の方は、「OpenFOAM6(Foundation版)」ですが、僕の使用環境はESI版の「OpenFOAM-v1912」で行っていますし、何よりWSLでやっています。

【使用環境】
Windows10
WSL (Windows Subsystem for Linux)
OpenFOAM-v1912

今回は、上記の環境で行いましたが参考記事通りにすると問題なく動きました。

また、OpenFOAMのバージョンは古いですがDEXCSのOpenFOAMをカスタマイズしたときの資料もあります。

 

やりたいこと:ベースで解く方程式+追加する方程式

まずは、今回の目的(やりたいこと)を確認しておきましょう。

【やりたいこと】

icoFoamのソルバに温度場の輸送方程式を追加する。

では、偏微分方程式を確認します。

ナビエストークス方程式

\begin{align*}
\frac{\partial \boldsymbol{u}}{\partial t}+\big(\boldsymbol{u}\cdot\nabla\big)\boldsymbol{u}= -\frac{1}{\rho}\nabla p+\nu \nabla^2\boldsymbol{u} \tag{1}
\end{align*}

👆これがもともとicoFoamで解く偏微分方程式です
ここに、以下の温度場輸送方程式を導入します。

温度場輸送方程式

\begin{align*}
\frac{\partial \boldsymbol{T}}{\partial t}+\big(\boldsymbol{u}\cdot\nabla\big)\boldsymbol{T}= D_t \nabla^2\boldsymbol{T} \tag{2}
\end{align*}

\(D_t =\frac{k}{\rho C_p}\)と置いています。

  • \(k\):熱伝導率
  • \(C_p\):定圧比熱
  • \(\rho\):密度

 

icoFoamってもともとどんな用途のためのソルバなのか?って思った場合は以下の表を参考にすれば良いでしょう。

OpenFOAMの「非圧縮性流体ソルバ」を以下に示しておきます。

参考サイト

非圧縮性流体
adjointShapeOptimizationFoam定常/非ニュートン流体の乱流/随伴方程式を使用することで障害物による圧損を考慮してダクト形状を最適化
boundaryFoam定常/1次元乱流/通常は流入口での境界層条件の生成に使用
icoFoam非定常/ニュートン流体の層流
nonNewtonianIcoFoam非定常/非ニュートン流体の層流
pimpleFoam非定常/時間刻み幅大/(PISOとSIMPLEを組み合わせた)PIMPLEアルゴリズム
pimpleDyMFoam非定常/移動メッシュ/(PISOとSIMPLEを組み合わせた)PIMPLEアルゴリズム
SRFPimpleFoam非定常/時間刻み幅大/単一の回転領域
pisoFoam非定常/PISOアルゴリズム
shallowWaterFoam非定常/回転を伴う非粘性の浅水方程式
simpleFoam定常
porousSimpleFoam定常/多孔質体(陰解法・陽解法)有りの乱流/MRF(Multiple Reference Frame)機能を使用可能
SRFSimpleFoam定常/単一の回転領域

この中で今回は「icoFoam」という、非定常・非圧縮性の流れを扱う際に使うソルバをベースにします。

非圧縮性の流体は基本的には温度場は解きません。
なぜなら、知りたい物理量は「流速\(\boldsymbol{u}\)」と「圧力\(p\)」の2つであり、解くべき方程式「運動量保存則」と「質量保存則」の2つです。

なので、方程式自体が閉じいるので温度を解く必要がないですし、密度変化が小さい非圧縮性の仮定が成り立つような状況において温度変化による密度への影響というのはとても小さいので温度を考える必要がないというのもあります。

こちらの記事もどうぞ

温度変化することによる密度変化を考慮しないといけない場合は、有名なのはブシネスク近似と呼ばれるものですが、ここでは深い話はしません。

以下の2点を頭においてくださればと思います。

  1. 温度が流れに沿って移動する方程式(温度の輸送方程式)を追加する。
  2. 温度が変化したことによる流速への影響は考慮していない。

全体の流れを確認

OpenFOAMのファイル構造は詳しくは理解していないですが、わからないなりに理解をしながら進めたいと思っています。

まず、今回やることの全体像は以下です。

各ステップの説明は、その都度する予定です。

今は、全体像を把握しておくだけで良いでしょう。

ライブラリの修正

  • 自身でカスタマイズする部分のライブラリは、環境変数「$WM_PROJECT_USER_DIR」というフォルダに作成していきます。
  • OpenFOAMの本体のライブラリは、環境変数「 $WM_PROJECT_DIR」にあります。

環境変数は前回の記事☟でまとめてありますので、確認しながら進めていってください。

前回の記事はこちら

 

環境変数が九人出来たら、実際にどのフォルダを指しているのかをWindows上でフォルダ移動して中身を確認してみればよいです。
絵の右側のフォルダに自身のカスタマイズしたいライブラリやソルバをぽいぽい入れていきます。

OpenFOAM本体と同じようなフォルダ構造を意識しながら、今回は「 /home/kamakiri/OpenFOAM/kamakiri-v1912/applications/solvers/incompressible/」を自身のカスタマイズしたフォルダに作成します。

作業ディレクトリを作成する

前回の記事でも書きましたが、WSL上でOpenFOAMのフォルダ・ファイル作成・コピー・削除を行うにはLinux上から行うのが安全です。

前回の記事はこちら

Windows上でフォルダ作成してもファイルに権限が付与されていないなど、権限を付与したりする手間があるので以下ではLinuxのコメンドベースでフォルダを操作していきます。

【やること】

  1. フォルダ移動:「/home/kamakiri/OpenFOAM/kamakiri-v1912」へ移動
  2. ファイル作成:「/solvers/incompressible」を作成

 

ターミナル上で以下のコマンドを打ちます。

ファイルできているか確認してみましょう。

フォルダの権限も確認しておきます。

ちゃんと権限が付与されています。

ベースコードを作業ディレクトリに複製する

では、ベースとなるソースコードを自身の作業ディレクトリにコピーしてソースコードを編集する準備をします。

【やること】

  1. ベースのコードを自身のカスタマイズフォルダへコピーする
  2. フォルダ名「icoFoam」を「my_icoFoam」に変更
  3. 「icoFoam.C」を「my_icoFoam.C」に変更

 

ターミナル上で以下のコマンドを打ちます。

フォルダがコピーされているか確認してみましょう。

コピーできていますね。

この「icoFoam」を「my_icoFoam」という名前に変えます。

ターミナル上で以下のコマンドを打ちます。

これでフォルダ名が変わりました。

では、「icoFoam.C」を「my_icoFoam.C」に変更

ファイル名が変わっていることを確認しましょう。

温度場輸送方程式のヘッダーファイルを準備

次に、温度場輸送方程式をC++で追加します。

温度場輸送方程式

\begin{align*}
\frac{\partial \boldsymbol{T}}{\partial t}+\big(\boldsymbol{u}\cdot\nabla\big)\boldsymbol{T}= D_t \nabla^2\boldsymbol{T} \tag{2}
\end{align*}

【やること】

  1. 「TEeq.H」に温度場輸送方程式のコードを追加する。

 

ターミナル上で以下のコマンドを打って「TEeq.H」ファイルを作成します。

「TEeq.H」ファイルの中身は以下のようにしておきます。

「TEeq.H」

編集は、Windows上のテキストエディタを使って良いです。
僕はVisual Stdio Codeを使っています。

※ちゃんと編集がLinux上で反映されているか心配な方は、「vi」コマンドでファイルの中身を確認してみても良いです。

いちおうちゃんとファイルの中身は反映されています。

※後々説明しますが、これを「my_icoFoam.C」でIncludeします。

「Dt(拡散係数)」と「T(温度)」の定義を行う

【やること】

  1. 「createField.H」に「Dt(拡散係数)」と「T(温度)」の変数定義。

 

「createField.H」
※「//Add here…」~「//Add End」までを追加します。

残念なことにC++の知識がないため、説明ができないですが「Dt(拡散係数)」と「T(温度)」の変数定義を行ったんだな~っていうのはわかりました。

ソースコードの修正

「createField.H」も「TEeq.H」もソースコードと言えばそうなのですが、今から修正する「my_icoFoam.C」は偏微分方程式を解いているメインのコードで、「createField.H」と「TEeq.H」をヘッダーファイルとして読み込んでいるファイルなので、いちおう区別しておきました。

 

【やること】

  1. 「my_icoFoam.C」に「TEeq.H」を読み込む。

 

「my_icoFoam.C」の中に、温度場輸送方程式を解く部分「TEeq.H」を読み込むようにすれば良いです。

my_icoFoam.C

 

「TEeq.H」を速度Uの偏微分方程式を解いている真下に追加しました。

C++のコードは理解していないですが、どこが時間発展のループを回りしている部分かはわかります。

これで一応ソースコードができました。

次に、この作成したソースコードをコンパイルする必要があります。

ソースコードのコンパイル

ソースコードは人間が読み書きできるプログラミング言語で書かれているのですが、それを機械が理解できるいわゆる「機械語」ってやつに変換する必要があります。

それは「Make」フォルダの中のファイルを編集することで行うことができます。

【やること】

  1. 「Make」フォルダの中の「files」を編集する。
  2. ソースコードのコンパイルして実行ファイルを作成

 

files

わかりにくいですが、以下の内容を書いていると理解しました。

  • どのソースコードをコンパイルする?:my_icoFoam.C
  • どこにコンパイルしてどこに何という名前で実行ファイルを保存する?:$(FOAM_USER_APPBIN)/my_icoFoam

「$FOAM_USER_APPBIN」という環境変数がわからない場合は、「echo $FOAM_USER_APPBIN」で調べてみましょう。

環境変数
$FOAM_USER_APPBIN/home/kamakiri/OpenFOAM/kamakiri-v1912/platforms/linux64Gcc63DPInt32Opt/bin

どうやら「/home/kamakiri/OpenFOAM/kamakiri-v1912/platforms/linux64Gcc63DPInt32Opt/bin」に「my_icoFoam」という実行ファイルが作成されるようです。

ソースコードのコンパイルには、「my_icoFoam」上で「wmake」とします。

ターミナルで以下のコンパイルを打ちます。

すると以下のような表示がされます。

これで「my_icoFoam」という名前の実行ファイルが作成できたと思います。

実行ファイルが確認できたかを確認するには「which」コマンドで実行ファイルの場所の保存場所を確認してみればよいでしょう。

どうやら「/home/kamakiri/OpenFOAM/kamakiri-v1912/platforms/linux64Gcc63DPInt32Opt/bin/」の中に「my_icoFoam」という実行ファイルが作成されているようですね。

これは、filesで「EXE = $(FOAM_USER_APPBIN)/my_icoFoam」と設定したことによってこのようになったということです。


※ちなみに「wmake」をすると以下のような、フォルダとファイルが作成されています。

「wmake」したときにエラーなどが生じてコンパイルをやり直したい場合は、こちらのファイルを削除する必要があります。
それには「wclean」と打てばフォルダが消えます。


以上で、OpenFOAMでのC++で書かれたソースコードの修正を行うことができました。

これで、温度場の輸送方程式も解くことができるようになりましたが・・・

OpenFOAMのケースファイルの方も「Dt(拡散係数)」「T(温度)」の設定を追加しなくてはいけません。

ケースファイルの修正

では、ケースファイルに「DT(拡散係数)」と「T(温度)」を追加します。

【やること】

  1. OpenFOAMチュートリアルのケースファイルを自身の作業ディレクトリにコピーする。
  2. T(温度)の境界条件の設定
  3. Dt(拡散係数)の定義
  4. 温度場の解き方の設定(fvSchemes、fvSolution)

 

まずは、「OpenFOAMチュートリアルのケースファイルを自身の作業ディレクトリにコピー」をします。

以下のように、ターミナルに打つことでOpenFOAMのチュートリアルを自身の作業ディレクトリにコピーすることができます。

その際には、以下の環境変数を利用すると便利です。

  • $FOAM_TUTORIALS
  • $FOAM_RUN

自身の作業ディレクトリは「20200515_my_icoFoam_cavity」としています。

コピーできているかを確認してみましょう。
まず「20200515_my_icoFoam_cavity」というフォルダができています。

そして、「20200515_my_icoFoam_cavity」の中身は、

こうなっています。

OpenFOAMのチュートリアルをコピーできています。

次に、この作業ディレクトリに移ります。

と打てば、「/home/kamakiri/OpenFOAM/kamakiri-v1912/run」に移動することができます。

あとは、「cd」コマンドでフォルダ移動を行ってください。

T(温度)の境界条件の設定

温度の設定は、「0」というフォルダの中の「p」ファイルを「T」という名前でコピーします。

ファイルの中身は以下のようにします。

T

 

Dt(拡散係数)の定義

拡散係数の定義は「constant」の中の「transportProperties」というファイルを編集します。

constant/transportProperties

「Dt [0 2 -1 0 0 0 0] 0.002;」という定義を追加しました。

次に、温度の輸送方程式をどのように解くかという設定を追加します。
そのためには、「system/fvSchemes」と「system/fvSolution」の2つのファイルを編集します。

system/fvSchemes

「div(phi,T) Gauss upwind;」としています。

これは\(\nabla\cdot(\boldsymbol{U}T)\)の移流項を解く部分の設定です。

system/fvSolution

温度場の解き方を追加しました。

これでケースファイルの編集は終わりです。

計算実行

では、計算を実行します。

【やること】

  1. 「blockMesh」でメッシュ作成
  2. 「my_icoFoam」で計算実行

 

まずは、blockMeshでメッシュ作成します。

問題が無ければ、以下のようになります。

 

次に、計算を実行します。
自身で作成した温度場輸送方程式も解く「my_icoFoam」という実行ファイルを実行します。

エラーが無ければこのように、「ダーーー」っと計算が流れるのがわかります。

計算がうまくいっていれば以下のように「0.1」「0.2」「0.3」「0.4」「0.5」というフォルダができています。

では、温度場が速度に応じて輸送しているのかをparaviewで確認してみましょう。

結果確認

 

結果確認にはparaviewを使用します。

OpenFOAMの結果を見るためには、以下のように「.foam」という拡張子のファイルを用意してして、それをparaviewで読み込めば結果を表示することができます。

なので、以下のように「a.foam」というからファイルを作成して、paraviewdで結果を表示してみます。

説明だとわかりにくいので動画を用意しました。

paraviewの表示にも温度「T」が追加されて、温度が速度に乗って解かれているのが確認できますね(^^)/

まとめ

  1. 簡単に温度場の輸送方程式をカスタマイズすることができました。
  2. 温度場が変化したことによる速度場への影響変化については解いていません。
  3. 今回はOpenFOAM v19-12でのカスタマイズを行ったが、バージョンが変わったことによる違いはわからないので今後の追跡調査が必要です。
  4. 思い通りのカスタマイズを行うには、OpenFOAMのファイル構成とC++の知識が必要であります。

OpenFOAMの古いバージョンですが、pisoFoamに温度場を追加したケースもやってみました。

C++の勉強について

自分自身がC++の知識が全くないので、いちから勉強する必要があります。

C++は、C言語にオブジェクト指向(クラスなどが使える)を追加して拡張された言語であるとのことなので、まずはC言語から勉強をはじめています。

「Cの絵本」はとてもわかりやすくて、すらすら読めます。
ページ数も多くないので、僕はこれを1週間くらいで読み終えました。

Cの絵本 第2版 C言語が好きになる新しい9つの扉

Cの絵本 第2版 C言語が好きになる新しい9つの扉

アンク
1,518円(10/30 02:04時点)
Amazonの情報を掲載しています

C言語の環境構築はWSLで使えるようにしています。エディタはVisual Stdio Codeです。
環境構築については、こちらのQiitaの記事が参考になります。

C言語もむちゃくちゃ詳しくなる必要はないと思いますが、一応コードを書いて計算実行したりして理解をしていきました。
「Cの絵本」だけでは詳しいところまでは学べないので、以下の参考書も手元に置いています。

新・明解C言語 入門編 (明解シリーズ)

新・明解C言語 入門編 (明解シリーズ)

柴田 望洋
3,860円(10/30 02:04時点)
Amazonの情報を掲載しています

ページ数も多く詳しく書かれていますし、説明もとてもわかりやすいです。

そして、C言語がそこそこ理解できたらいよいよC++の勉強を始めていきます(👈今、ここ)
同じく「C++の絵本」を買って読もうと思ったのですが、あまりサンプルコードが載っていなかったので読むのをやめました(「Cの絵本」はサンプルコードあったのに・・)

以下の、「猫でもわかるシリーズ」で勉強を進めています。

猫でもわかるC++プログラミング 第2版 (猫でもわかるプログラミング)

猫でもわかるC++プログラミング 第2版 (猫でもわかるプログラミング)

粂井 康孝
2,420円(10/30 02:04時点)
Amazonの情報を掲載しています

これが結構わかりやすいですし、サンプルコードも程よい量で載っています。

あとは、やっぱり詳しいところまで知りたい場合は以下の参考書に頼ります。

新・明解C++入門編 (明解シリーズ)

新・明解C++入門編 (明解シリーズ)

柴田 望洋
3,025円(10/30 16:06時点)
Amazonの情報を掲載しています

あとは、Twitterで教えてもらったこれとかを手元に置いています。

ロベールのC++入門講座

ロベールのC++入門講座

ロベール
4,054円(10/30 12:51時点)
Amazonの情報を掲載しています

では、OpenFOAMのC++はかなり特殊らしいのですが、C++の基礎知識を身につけないことにはどうしようもないのでまずは2週間程度でC++の基礎を身につけたいと思います(^^)/

OpenFOAMの参考書

OpenFOAMの参考書と言えば、こちらのサイトを書籍にまとめなおした☟これがおすすめです。

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

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

3,520円(10/30 04:54時点)
Amazonの情報を掲載しています

OpenFOAMの使用にあたっては「どうやってOpenFOAMをインストールするの?」というところでつまずきます。

そんな時は、以下の書籍をおすすめします。

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

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

川畑 真一
3,740円(10/30 16:06時点)
Amazonの情報を掲載しています

インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。
著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。

関連記事もどうぞ

POSTED COMMENT

  1. uppoi より:

    v2006でもできました。ありがとうございます!

COMMENT

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