fixedFluxPressure 境界条件

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

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

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

fixedFluxPressure 境界条件

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

\begin{equation}
{\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}
\end{equation}

この式\eqref{eq:2}を変形すると，次式\eqref{eq:3}が得られます．

\begin{equation}
(\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}
\end{equation}

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 境界条件が，重力という体積力を伴う計算において，圧力勾配を適切に設定するための条件であるのに対し，この新しい境界条件は体積力の有無に関わらずより汎用的な圧力境界条件として設計されているようです．まだテストの段階であり，実用レベルに至るには更なる検証計算が必要とのことです．この境界条件についても調査して別の投稿で取り上げたいと思います．