Category: Japanese

前回の投稿 の続きとして，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}式中で使用されている２つのベクトル \(\bf a\) と \(\bf b\) の二項積\(\otimes\) (dyadic product) は次式\eqref{eq:dyad}のように計算されます．

$$ {\bf a} \otimes {\bf b} = \left(
\begin{array}{ccc}
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
\end{array}
\right) \tag{2} \label{eq:dyad} $$

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

\begin{align} \left({\bf a} \otimes {\bf b}\right) {\bf c} &= \left(
\begin{array}{ccc}
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
\end{array}
\right)\left(
\begin{array}{c}
c_1 \\
c_2 \\
c_3
\end{array}
\right) \\
& = {\bf a}\left({\bf b} \cdot {\bf c}\right) \tag{3} \label{eq:D}
\end{align}

準備は以上です．

では本題です．関係式\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\) にこれを作用させてみます．

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

ここで，１つ目の変形で関係式\eqref{eq:D}を使用しています．そして，３つ目の変形に近似が含まれています．これを\eqref{eq:reconstruct}式に合わせて書き直すと，次式が得られます．
$$ \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

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

これは，具体的にどのような演算を行っているのでしょうか？こんな時は，ソースコードを確認しましょう．ソースコードを参照して実装方法をすぐに確認できるのが，オープンソースの１つのメリットですね．この演算子のソースファイルは，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}\) を使用することの数値計算上の利点に関する説明につなげられたらと考えています．では！