Conduction, convection and radiation in OpenFOAM (Under construction)

This page is under construction. I will update by adding a description of each solver and model.



Conduction + Convection (Conjugate Heat Transfer)
  • chtMultiRegionSimpleFoam (Steady)
  • chtMultiRegionFoam (Transient)
+ Radiation

All the above solvers but laplacianFoam are able to deal with the radiative heat transfer. There are the following five (virtually four) models available in OpenFOAM. Their source code is located in src/thermophysicalModels/radiation/radiationModels and we can see the brief descriptions in the header file of each radiation class.

  • P1

Keywords: optical thickness [2]

  • fvDOM (Finite Volume Discrete Ordinates Method)

  • viewFactor

  • opaqueSolid

  • none

The settings of the radiation models are described in constant/radiationProperties file.


Solvers for heat transfer problems in OpenFOAM – buoyantBoussinesqPimpleFoam

There are several solvers for solving the heat transfer problems in OpenFOAM. In this post, I will try to give a description of how to derive the governing equations of the following two solvers:

  • buoyantBoussinesqSimpleFoam
  • buoyantBoussinesqPimpleFoam

In the presence of gravity body force, the mass and momentum conservation equations are:

\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \boldsymbol{u}\right) = 0, \tag{1} \label{eq:massEqn}

\frac{\partial \left(\rho \boldsymbol{u} \right)}{\partial t} + \nabla \cdot \left(\rho\boldsymbol{u}\boldsymbol{u}\right) = &-\nabla p + \rho \boldsymbol{g} + \nabla \cdot \left(2 \mu_{eff} D \left( \boldsymbol{u} \right) \right) \\
&-\nabla \left( \frac{2}{3}\mu_{eff}\left(\nabla \cdot \boldsymbol{u} \right) \right), \tag{2} \label{eq:momEqn}

where \(\boldsymbol{u}\) is the velocity field, \(p\) is the pressure field, \(\rho\) is the density field, and \(\boldsymbol{g}\) is the gravitational acceleration. The effective viscosity \(\mu_{eff}\) is the sum of the molecular and turbulent viscosity and the rate of strain tensor \(D\left(\boldsymbol{u}\right)\) is defined as \(D\left(\boldsymbol{u}\right) = \frac{1}{2}\left( \nabla \boldsymbol{u} + \left(\nabla \boldsymbol{u}\right)^{T} \right)\).

If both density and gravitational acceleration are constant, the gravitational force can be expressed in terms of gradient:

\rho \boldsymbol{g} = \nabla \left( \rho \boldsymbol{g} \cdot \boldsymbol{r} \right), \tag{3} \label{eq:Fgravity}

where \(\boldsymbol{r}\) is the position vector, so that the pressure gradient and gravity force can be lumped together as shown in the following equation:

\nabla p \;-\; \rho \boldsymbol{g} = \nabla \left( p \;-\; \rho \boldsymbol{g} \cdot \boldsymbol{r} \right). \tag{4} \label{eq:conversion}

Boussinesq approximation

Now we consider the case when the density is not constant. The Boussinesq approximation is valid when the variation of the density induced by the temperature change is small, as Ferziger and Peric (2001, pp.14-15) states that:

In flows accompanied by heat transfer, the fluid properties are normally functions of temperature. The variations may be small and yet be the cause of the fluid motion. If the density variation is not large, one may treat the density as constant in the unsteady and convection terms, and treat it as variable only in the gravitational term. This is called the Boussinesq approximation.

Hereafter, we denote the reference density by \(\rho_0\) at the reference temperature \(T_0\). If we replace \(\rho\) by \(\rho_0\) in the Eq.\eqref{eq:massEqn}\eqref{eq:momEqn} except for the gravitational term, then we get:

\nabla \cdot \boldsymbol{u} = 0, \tag{5} \label{eq:massEqn2}

\frac{\partial \left(\rho_0 \boldsymbol{u} \right)}{\partial t} + \nabla \cdot \left(\rho_0\boldsymbol{u}\boldsymbol{u}\right) = -\nabla p + \rho \boldsymbol{g} + \nabla \cdot \left(2 \mu_{eff} D \left( \boldsymbol{u} \right) \right). \tag{6} \label{eq:momEqn2}

Then divide both sides of Eq.\eqref{eq:momEqn2} by \(\rho_0\) and we get:

\frac{\partial \boldsymbol{u}}{\partial t} + \nabla \cdot \left(\boldsymbol{u}\boldsymbol{u}\right) = -\frac{1}{\rho_0} \left( \nabla p \;-\; \rho\boldsymbol{g} \right) + \nabla \cdot \left(2 \nu_{eff} D \left( \boldsymbol{u} \right) \right). \tag{7} \label{eq:momEqn3}

Here, the density \(\rho\) in the gravitational term is expressed by the linear function of the temperature \(T\):

\rho \approx \rho_0\left[ 1 \;-\; \beta\left( T \;-\; T_0 \right) \right]. \tag{8} \label{eq:rhoEqn}

This relation is derived from the definition of the volumetric thermal expansion coefficient \(\beta\):

\beta \equiv -\frac{1}{\rho}\frac{\partial \rho}{\partial T}
\approx -\frac{1}{\rho_0}\frac{\rho \;-\; \rho_0}{T \;-\; T_0}. \tag{9} \label{eq:beta}

We have to know when these approximations are valid. Ferziger and Peric (2001, p.15) states that:

This approximation introduces errors of the order of 1% if the temperature differences are below e.g. 2℃ for water and 15℃ for air. The error may be more substantial when temperature differences are larger; the solution may even be qualitatively wrong.

If temperature differences are larger, we should use buoyantSimpleFoam or buoyantPimpleFoam without Boussinesq approximation.

In terms of the implementation, the pressure gradient and gravity force terms are rearranged in the following form:

&-\nabla \left(\frac{p}{\rho_0}\right) + \left(\frac{\rho}{\rho_0}\right)\boldsymbol{g} \\
=&-\nabla \left(\frac{p \;- \rho \boldsymbol{g} \cdot \boldsymbol{r}}{\rho_0} + \frac{\rho \boldsymbol{g} \cdot \boldsymbol{r}}{\rho_0} \right) + \left(\frac{\rho}{\rho_0}\right)\boldsymbol{g} \tag{10a}\\
=&-\nabla p_{rgh} – \left( \boldsymbol{g} \cdot \boldsymbol{r} \right) \nabla \left(\frac{\rho}{\rho_0}\right), \tag{10b} \label{eq:final}

where we define a new variable \(p_{rgh} = \left(p-\rho\boldsymbol{g}\cdot\boldsymbol{r}\right)/\rho_0\) and finally we get:

\frac{\partial \boldsymbol{u}}{\partial t} + \nabla \cdot \left(\boldsymbol{u}\boldsymbol{u}\right) =\; &- \nabla p_{rgh} – \left( \boldsymbol{g} \cdot \boldsymbol{r} \right) \nabla \left(\frac{\rho}{\rho_0}\right) \\ &+ \nabla \cdot \left(2 \nu_{eff} D \left( \boldsymbol{u} \right) \right). \tag{11} \label{eq:momEqn4}

This momentum equation \eqref{eq:momEqn4} is implemented in UEqn.H as shown in the following code:

The Eq. \eqref{eq:final} is calculated using the fvc::reconstruct function from the face flux fields and the variable rhok is defined in createFields.H as \(\rho/\rho_0\):


[1] Joel H. Ferziger and Milovan Peric (2001): Computational Methods for Fluid Dynamics

Please wait till the next update!

reconstruct 演算子について #2

前回の投稿 の続きとして,reconstruct 演算子をもう少し詳しく見ていきます.前回の投稿で掲載した reconstruct 演算子の定義式は,数式では次式\eqref{eq:reconstruct}のように表現できます.

$$ reconstruct(ssf) = \left(\sum_{f}^{}{\left({\bf n}_f \otimes {\bf S}_f\right)} \right)^{-1}\left(\sum_{f}^{} {\bf n}_f ssf\right) \tag{1} \label{eq:reconstruct} $$

ここで,\({\bf n}_f\) は,フェイスの単位法線ベクトルを,\({\bf S}_f\left(=|{\bf S}_f|{\bf n}_f\right)\) は,フェイスの面積を大きさとし,向きが \({\bf n}_f\) と等しい法線ベクトルをそれぞれ表します.また,surfaceSum 演算子は,各セルにおいてそのセルのフェイスの値の総和を計算します(fvcSurfaceIntegrate.C).


\eqref{eq:reconstruct}式中で使用されている2つのベクトル \(\bf a\) と \(\bf b\) の二項積\(\otimes\) (dyadic product) は次式\eqref{eq:dyad}のように計算されます.

$$ {\bf a} \otimes {\bf b} = \left(
a_1b_1 & a_1b_2 & a_1b_3 \\
a_2b_1 & a_2b_2 & a_2b_3 \\
a_3b_1 & a_3b_2 & a_3b_3
\right) \tag{2} \label{eq:dyad} $$

この二項積とベクトル \(\bf c\) との積を計算すると,次式\eqref{eq:D}のように変形できます.

\begin{align} \left({\bf a} \otimes {\bf b}\right) {\bf c} &= \left(
a_1b_1 & a_1b_2 & a_1b_3 \\
a_2b_1 & a_2b_2 & a_2b_3 \\
a_3b_1 & a_3b_2 & a_3b_3
c_1 \\
c_2 \\
\right) \\
& = {\bf a}\left({\bf b} \cdot {\bf c}\right) \tag{3} \label{eq:D}



では本題です.関係式\eqref{eq:D}を使用して,reconstruct 演算子を調べていきます.ここで,\eqref{eq:reconstruct}式中の以下の行列の部分を

$$ D \equiv \sum_{f}^{}{\left({\bf n}_f \otimes {\bf S}_f\right)} \tag{4} \label{eq:Dop} $$

として,密度 \(\rho\) の勾配 \(\nabla \rho\) にこれを作用させてみます.

D \nabla \rho &= \sum_{f}^{} {\left({\bf n}_f\left({\bf S}_f \cdot \nabla \rho \right)\right)} \tag{5a}\\
&= \sum_{f}^{} {\left({\bf n}_f\left(|{\bf S}_f|{\bf n}_f \cdot \nabla \rho \right)\right)} \tag{5b}\\
&\approx \sum_{f}^{} {\left({\bf n}_f\left(|{\bf S}_f|{\bf n}_f \cdot \left(\nabla \rho\right)_f \right)\right)} \tag{5c}\\
& = \sum_{f}^{} {\left({\bf n}_f\left(|{\bf S}_f|snGrad\left(\rho\right)_f \right)\right)} \tag{5d} \label{eq:conversion}

$$ \nabla \rho \approx D^{-1}\left(\sum_{f}^{} {\left({\bf n}_f|{\bf S}_f|snGrad\left(\rho\right)_f \right)}\right) \tag{6} \label{eq:final} $$


[1] CFD Online Forums: fvc::reconstruct( ) algorithm

reconstruct 演算子について #1

OpenFOAM のソースコードを眺めていると,fvc::reconstruct という演算子を目にします.例えば,buoyantSimpleFoam の UEqn.H ではこの演算子が下記のように使用されています.

これは,具体的にどのような演算を行っているのでしょうか?こんな時は,ソースコードを確認しましょう.ソースコードを参照して実装方法をすぐに確認できるのが,オープンソースの1つのメリットですね.この演算子のソースファイルは,fvcReconstruct.C です.このファイルの 76 行目に演算の定義が記述されています.

私は一見しただけでは何を行っているのか分かりませんでした.ヘッダーファイル fvcReconstruct.H も確認してみます.

冒頭の Description の項目に,Reconstruct volField from a face flux field という記述があります.最初に挙げた buoyantSimpleFoam での使用例を確認すると,fvc::reconstruct(引数) の引数は,フェイス(セル界面)で値をもっています(face flux field).具体的には,

  • snGrad(rho) は,フェイスの法線方向の密度の勾配を,
  • mesh.magSf() は,フェイスの面積を

それぞれ計算するので,両方とも値はフェイスで保持しています.この例では,OpenFOAM の用語を使用すると,surfaceScalarField を引数をとして受け取り,セル中心で値をもつ volVectorField を返しています.


次回の投稿では,簡単な例でこの演算子を書き下して,数式を使用した説明にチャレンジしたいと思います.最終的には,重力のような体積力を考慮した計算において,圧力の変数として \(p\) ではなく \(p_{rgh}\) を使用することの数値計算上の利点に関する説明につなげられたらと考えています.では!

fixedFluxPressure 境界条件

OpenFOAM で自然対流を伴う熱流体計算や混相流 (Multiphase) 計算のように重力を考慮する計算を行う場合には,流れ場を表す変数として静圧 \(p\) そのものではなく,次式\eqref{eq:prgh}で定義される変数 \(p_{rgh}\) を使用します.

p_{rgh} = p\;-\;\rho{\bf g} \cdot {\bf x} \tag{1} \label{eq:prgh}

ここで,\(\rho\) は流体の密度を,\(\boldsymbol{g}\) は重力加速度ベクトルを,\(\boldsymbol{x}\) は位置ベクトルをそれぞれ表します.OpenFOAM で p_rgh と表されるこの変数の境界条件として,以前は,buoyantPressure 条件が使用されていましたが,v2.3 以降では,fixedFluxPressure 条件が使用されています.今回の投稿では,この2つの境界条件について取り上げます.

fixedFluxPressure 境界条件

セル界面の流速 \(\boldsymbol{U}_f\) については,次式\eqref{eq:2}が成り立ちます.

{\bf U}_f \cdot {\bf S}_f = \left(\frac{{\bf H}}{A_p}\right)_f \cdot {\bf S}_f\;-\;\left(\frac{1}{A_p}\right)_f (\nabla p_{rgh})_f \cdot {\bf S}_f \tag{2} \label{eq:2}


(\nabla p_{rgh})_f \cdot {\bf n}_f = \frac{1}{|{\bf S}_f|\left(\frac{1}{A_p}\right)_f} \left(\left(\frac{{\bf H}}{A_p}\right)_f \cdot {\bf S}_f\;-\;{\bf U}_f \cdot {\bf S}_f \right) \tag{3} \label{eq:3}

fixedFluxPressure 境界条件は,その境界で課されている流速 \(\boldsymbol{U}\) に対する境界条件から決まる流束を満たすように,境界の法線方向の圧力勾配を式\eqref{eq:3}の右辺により計算します.

この境界条件の updateCoeffs メンバ関数の 137行目の勾配値 snGradp は,最新の開発版では,constrainPressure 関数 (constrainPressure.C) で計算されます.式\eqref{eq:3}の右辺の式が,以下のソースコードの 68-73 行目に対応しています.

それ以前のバージョンでは,各ソルバーの pEqn.H において,境界での勾配計算が行われています.例えば,v3.0.x の interFoam の場合には,以下の部分でこの計算を行っています.

buoyantPressure 境界条件

この境界条件は,2.3.x 以降のバージョンには実装されていませんが,ここで簡単に触れたいと思います.この境界条件は,変数 p_rgh の境界法線方向 \({\bf n}\) の勾配を次式\eqref{eq:4}により計算します.

$$ \frac{\partial (p_{rgh})}{\partial {\bf n}} =  -({\bf g} \cdot {\bf x}) \frac{\partial \rho}{\partial {\bf n}} \tag{4} \label{eq:4} $$

したがって,静圧 \(p\) の境界法線方向の勾配は次式\eqref{eq:5}で与えられます.

$$ \frac{\partial p}{\partial {\bf n}} = \frac{\partial (p_{rgh})}{\partial {\bf n}}\;+\;\frac{\partial (rgh)}{\partial {\bf n}} = \rho ({\bf g} \cdot {\bf n}) \tag{5} \label{eq:5} $$

よって,法線方向 \({\bf n}\) が重力ベクトル \({\bf g}\) に直交するような境界においては,静圧 \(p\) に関して勾配0条件を意味します.

$$ \frac{\partial p}{\partial {\bf n}} = 0 \tag{6} \label{eq:6} $$


$$ \frac{\partial p}{\partial {\bf n}} = \pm \rho |{\bf g}| \tag{7} \label{eq:7} $$

式\eqref{eq:4}は,buoyantPressureFvPatchScalarField クラスのメンバ関数 updateCoeffs に実装されています.


最新の開発版では fixedFluxExtrapolatedPressure という名前の境界条件が新たに実装されています.fixedFluxPressure 境界条件が,重力という体積力を伴う計算において,圧力勾配を適切に設定するための条件であるのに対し,この新しい境界条件は体積力の有無に関わらずより汎用的な圧力境界条件として設計されているようです.まだテストの段階であり,実用レベルに至るには更なる検証計算が必要とのことです.この境界条件についても調査して別の投稿で取り上げたいと思います.

OpenFOAM のベクトル演算

OpenFOAM には,ベクトルに対する様々な演算子が用意されています.次のプログラムを実行して,それぞれの演算子の意味を確認しましょう.

Test Program


Meaning of the Operations
    • cmptMax

$$cmptMax({\bf a}) = max (a_{1},\;a_{2},\;a_{3})$$

    • cmptMin

$$cmptMin({\bf a}) = min (a_{1},\;a_{2},\;a_{3})$$

    • cmptSum

$$cmptSum({\bf a}) = a_{1} + a_{2} + a_{3}$$

    • cmptAv

$$cmptAv({\bf a}) = \frac{a_{1} + a_{2} + a_{3}}{3}$$

    • cmptProduct

$$cmptProduct({\bf a}) = a_{1} \times a_{2} \times a_{3}$$

    • cmptMag

$$cmptMag({\bf a}) = (|a_{1}|\;\;|a_{2}|\;\;|a_{3}|)$$

    • cmptMultiply

$$cmptMultiply({\bf b},\;{\bf c}) = (b_{1}\times c_{1}\;\;b_{2}\times c_{2}\;\;b_{3}\times c_{3})$$

    • max

$$max({\bf b},\;{\bf c}) = (max(b_{1},\;c_{1})\;\;max(b_{2},\;c_{2})\;\;max(b_{3},\;c_{3}))$$

    • min

$$min({\bf b},\;{\bf c}) = (min(b_{1},\;c_{1})\;\;min(b_{2},\;c_{2})\;\;min(b_{3},\;c_{3}))$$

これらのベクトル演算子は,VectorSpaceI.H ファイルで実装されています.

境界に隣接するセルからなる cellZone の生成

ある境界に隣接するセルのみから構成される cellZone は,setSet ユーティリティを使用して,以下のような手順で生成することができます.

1行目のコマンドで,境界名が duct の境界フェイスから,ductwall という名前の faceSet を生成しています.2行目のコマンドでは,この faceSet に属するフェイスの owner セルから構成される cellSet を生成し,layer1Cells と名付けています.最後に3行目のコマンドで,この cellSet から,layer1Zone という名前の cellZone を生成しています.


ここで,2行目のコマンドで使用している faceToCell オプションは,指定した faceSet に属する faceneighbour セル,owner セル,またはその両方のセル (any) から構成される cellSet を生成することができます.また,all オプションを指定した場合には,全てのフェイスが指定した faceSet に属するセルから構成される cellSet を生成します.境界フェイスの場合,owner セルしか存在しないので,ownerany のどちらを指定しても同じ cellSet が生成されます.

実装は,faceToCell.C ファイルで確認できます.

OpenFOAM では,それぞれのセルのフェイスの両側に位置する2つのセルを,owner セルと neighbour セルと呼んでいます(境界フェイスの場合は,neighbour セルのみ存在します).ここで,セル番号が小さい方のセルが owner セルです.ここら辺の話題については,こちらのスライドで触れています.

neighbour のスペルちょっと違和感ありますか?OpenFOAM Code Style Guide の最後の Orthography の項でも言及されているように,OpenFOAM はもともと,英国産のソフトウェアなので,今でも英国式のスペリングが好まれています.他には,centre なども頻出しますね.

cellSource を使用した cellZone におけるセル値の演算処理

OpenFOAM の Function Objects 機能を活用することで,結果処理の効率化が可能です.その1つである cellSource を使用すると,ユーザーが指定した cellZone (セルの集合) に属するセルの変数値に対して,セル体積を重みとした加重平均や体積積分などの処理を計算実行中に行って,値をモニタリングすることが可能です.例えば,porous という名前の cellZone において体積を重みとした温度の加重平均値を計算したい場合には,ケースの system/controlDict ファイルに下記のような設定を行います.

実行する演算処理を operation に指定します.選択可能な 11 個のオプションが,cellSource.C ファイルで確認できます.

それぞれの演算オプションの実装は,cellSourceTemplates.C ファイルで確認できます.

ここで,V はセル体積を表しているので,volIntegrate オプションは,体積積分を計算します.volAverage は,この体積積分値を cellZone の合計体積で割ることで,体積を重みとした加重平均値を計算します.これに対して,average は,単なる算術平均を計算します.また,名前に weighted が付いているオプションでは,加重平均の重みに使用する変数を weightField に指定します.sumMag オプションで使用されている cmptMag 演算子についてはこちらの記事で取り上げています.




Structure of src Directory in OpenFOAM part1

OpenFOAM の各種ライブラリのソースコードは,src ディレクトリ に設置されています.例えば,2015年のオープン CAE シンポジウムの講習会で取り上げた境界条件のライブラリ libfiniteVolume のソースコードは,src/finiteVolume/ 以下に設置されています.第1回目の今回は,以下の4つのディレクトリを取り上げます.


既存のソルバーのソースコードを直接編集することなく支配方程式にソース項(ポーラスメディアの抵抗ソースや熱交換器の体積熱ソースなど)を追加するなどのカスタマイズを行うことができる fvOptions 機能のライブラリ libfvOptions のソースコードが設置されています.各ソルバーが fvOptions 機能に対応しているかどうかは,ソルバーのソースコードで確認できます.下記は,simpleFoam の例です.

fvOptions で使用可能なタイプのリストを,こちらのページにまとめています.


セル番号の付け替え(Ordering)を行うための手法(Cuthill-McKeeSloan アルゴリズムなど)のソースコードが設置されています.renumberMesh ユーティリティは,このライブラリ librenumberMethods を使用して,セル番号の付け替えにより係数行列のバンド幅を小さくすることで,連立方程式の反復解法の計算時間の短縮化を図ります.デフォルトの手法である Cuthill-McKee 以外の手法を使用するためには,設定ファイルを用意し,実行時に -dict オプションでそのパスを指定します.

  • CuthillMcKee: デフォルトの手法
  • Sloan: boost C++ ライブラリの Sloan アルゴリズムを使用
  • manual
  • random
  • spring
  • structured

最初の2つの手法は,特別なパラメータ等のコントロールなしに使用することができます.それぞれ,バンド幅を小さくするのか,profile を小さくするのか等の特徴があります.


OpenFOAM では,メッシュの移動変形機能を Dynamic Mesh という名前で呼んでいます.並進・回転などの剛体運動,1次元方向の伸び縮み運動など,トポロジー変化を伴わないメッシュの移動変形のライブラリのソースコードが設置されています.次のスライドで詳しく説明しています.

[slideshare id=31811977&doc=openfoamdynamicmesh20140302-140302022414-phpapp01&w=500]



上記のリンクのように,OpenFOAM のリリースノートには,新機能や機能変更に関して比較的詳しい説明が記載されています.自分が使用しているバージョンに対応する情報を収集するためにも,良くチェックするようにしたいですね.