【OpenFOAMv2012】linearの解説

OpenFOAMで実装されている離散化スキームである linear について解説する

有限体積法による離散化

コンピュータで偏微分方程式を解くためには，計算格子上に偏微分方程式を離散化する必要がある

OpenFOAMで用いられている離散化手法である有限体積法を用いると，輸送方程式，発散，勾配は次のように離散化できる

（輸送方程式）

\begin{align}
&\frac{\partial}{\partial t}(\rho\phi)
+\nabla\cdot(\rho\mathbf{u}\phi)
=s
+\nabla\cdot\left(\Gamma\nabla \phi\right) \\
\\ \rightarrow \quad
&\frac{\partial}{\partial t}(\rho_{_{P}}\phi_{_{P}}V_{_{P}})
+\sum_{f}{F\phi_{f}}
=s_{_{P}}V_{_{P}}
+\sum_{f}{\Gamma_{f}\left|S_{f}\right|\frac{\phi_{_{N}}-\phi_{_{P}}}{|\mathbf{x}_{_{N}}-\mathbf{x}_{_{P}}|}}
\end{align}

（発散）

\begin{equation}
\nabla\cdot\phi \quad \rightarrow \quad \sum_{f}\mathbf{S}_{f}\cdot\phi_{f}
\end{equation}

（勾配）

\begin{equation}
\nabla\phi \quad \rightarrow \quad \sum_{f}\mathbf{S}_{f}\phi_{f}
\end{equation}

ここで，添え字\(_P\)，\(_N\)，\(_f\) はそれぞれ注目セル\(P\)，隣接セル\(N\)，注目セルと隣接セルの境界面\(f\)における代表値であることを表し，\(\phi\)は何かしらの物理量，\(V\)はセル体積，\(\Gamma\)は拡散係数である

また，\(|S_{f}|\)は注目セルと隣接セルの境界面の面積，\(\mathbf{S}_{f}\)は境界面に垂直で境界面の面積\(|S_{f}|\)を大きさに持つ法線ベクトルであり，\(\sum_{f}\)はセルがもつすべての面に対する総和を表す

上記3つの離散化を見てみると未知数は\(\phi_{f}\)のみであり，有限体積法の離散化は「境界面上の代表値\(\phi_{f}\)をどのように補間するか」という点に集約されていることがわかる

セルが不規則に配置されている非構造格子では，\(\phi_{f}\)の離散化に使える値は注目セル\(P\)の代表値\(\phi_{_{P}}\)と隣接セル\(N\)の代表値\(\phi_{_{N}}\)の2つのみである

\(\phi_{f}\)を\(\phi_{_P}\)と\(\phi_{_N}\)によって補間すると次の連立方程式が得られる

\begin{equation}
\mathbf{A}_{_{P}}\phi_{_{P}}+\sum{\mathbf{A}_{_N}\phi_{_N}}
=\mathbf{b}
\end{equation}

\(\mathbf{A}\)や\(\mathbf{b}\)の具体的な値を決定するのは偏微分方程式の離散化に用いられた離散化スキームであり，離散化スキームの精度や性質によって計算の安定性や結果の精度は大きく変化する

↓離散化の導出についてはこちら

有限体積法による偏微分方程式の離散化②

一般的な移流拡散方程式を有限体積法によって離散化する方法について説明する

mtkbirdman.com

2021.12.14

中心差分

中心差分（central difference：CD）は2次精度であり，対流項，発散，勾配，すべてに使用できる有限体積法で最もポピュラーな離散化スキームである

中心差分では，\(\phi_{f}\)を次のように補完する

\begin{equation}
\phi_{f}=w\phi_{_{P}}+\left(1-w\right)\phi_{_{N}}, \qquad
w=\frac{\left|\mathbf{x}_{f}-\mathbf{x}_{_{N}}\right|}{\left|\mathbf{x}_{_{P}}-\mathbf{x}_{_{N}}\right|}
\end{equation}

↓詳しくはこちら

有限体積法における離散化スキームの解説と導出

有限体積法における離散化スキームについての解説と導出を行う

mtkbirdman.com

2021.12.14

OpenFOAMにおける実装

OpenFOAMでは，中心差分は以下のように実装される

\begin{align}
&\phi_{f}=w\left(\phi_{_{P}}-\phi_{_{N}}\right)+\phi_{_{N}} \tag{1} \\ \\
&w=\frac{\left|\mathbf{x}_{f}-\mathbf{x}_{_{N}}\right|}{\left|\mathbf{x}_{_{P}}-\mathbf{x}_{_{N}}\right|} \tag{2}
\end{align}

参考
≫User Guide: Linear divergence scheme - OpenFOAM

実装は主に以下のファイルに記述されている

≫surfaceInterpolationScheme.C：離散化スキームの親クラス
≫linear.H：子クラス

surfaceInterpolationScheme.C で補間式(1)が定義され，linear.H で重み係数\(w\)が定義されている

それでは1つずつ確認していく

重み係数

重み係数は linear.H の weights で次のように定義されている

\begin{equation}
w=\frac{\left|\mathbf{x}_{f}-\mathbf{x}_{_{N}}\right|}{\left|\mathbf{x}_{_{P}}-\mathbf{x}_{_{N}}\right|} \tag{2}
\end{equation}

//- Return the interpolation weighting factors
tmp<surfaceScalarField> weights
(
    const GeometricField<Type, fvPatchField, volMesh>&
) const
{
    return this->mesh().surfaceInterpolation::weights();
}

//- Return reference to linear difference weighting factors
virtual const surfaceScalarField& weights() const;

const Foam::surfaceScalarField& Foam::surfaceInterpolation::weights() const
{
    if (!weights_.valid())
    {
        weights_.set(geometry().weights().ptr());
    }

    return weights_();
}

補間式

補間式は surfaceInterpolationScheme.C の中の dotInterpolate() で次のように定義されている

\begin{align}
&\phi_{f}=w\left(\phi_{_{P}}-\phi_{_{N}}\right)+\phi_{_{N}}
\end{align}

プログラムでは \(\mathrm{vf}=\phi\)，\(\mathrm{lambdas}=w\) という変数が使われている

typedef typename Foam::innerProduct<typename SFType::value_type, Type>::type
    RetType;

const surfaceScalarField& lambdas = tlambdas(); // 重み係数

const Field<Type>& vfi = vf; // 物理量
const scalarField& lambda = lambdas;

const fvMesh& mesh = vf.mesh();
const labelUList& P = mesh.owner(); // 注目セルPのインデックス
const labelUList& N = mesh.neighbour(); // 隣接セルNのインデックス

tmp<GeometricField<RetType, fvsPatchField, surfaceMesh>> tsf // 補完した値を入れるための変数
(

    new GeometricField<Type, fvsPatchField, surfaceMesh>
    (
        IOobject
        (
            "interpolate("+vf.name()+')',
            vf.instance(),
            vf.db()
        ),
        mesh,
        vf.dimensions()
    )
);
GeometricField<RetType, fvsPatchField, surfaceMesh>& sf = tsf.ref();

Field<RetType>& sfi = sf.primitiveFieldRef(); // 境界面上で補間された値

const typename SFType::Internal& Sfi = Sf(); // 1

for (label fi=0; fi<P.size(); fi++) // 境界面上の値を補間
{
    sfi[fi] = Sfi[fi] & (lambda[fi]*(vfi[P[fi]] - vfi[N[fi]]) + vfi[N[fi]]);
}

//- Return the face-interpolate of the given cell field
//  with the given weighting factors
static tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
interpolate
(
    const GeometricField<Type, fvPatchField, volMesh>&,
    const tmp<surfaceScalarField>&
);

template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
Foam::surfaceInterpolationScheme<Type>::interpolate
(
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const tmp<surfaceScalarField>& tlambdas
)
{
    return dotInterpolate(geometricOneField(), vf, tlambdas);
}

//- Return the face-interpolate of the given cell field
//  with the given weighting factors dotted with given field Sf
template<class SFType>
static tmp
<
    GeometricField
    <
        typename innerProduct<typename SFType::value_type, Type>::type,
        fvsPatchField,
        surfaceMesh
    >
>
dotInterpolate
(
    const SFType& Sf,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const tmp<surfaceScalarField>& tlambdas
);

まとめ

OpenFOAMで実装されている離散化スキームである linear について解説した

↓その他の離散化スキーム

【まとめ】OpenFOAMの非圧縮ソルバーまとめ

OpenFOAMで使う非圧縮ソルバーについて解説している記事のまとめ

mtkbirdman.com

2021.12.15

↓おすすめ記事

数値流体力学のおすすめ参考書

大学院の研究でCFDを使うことになったときに参考にした書籍を紹介する

mtkbirdman.com

2021.11.04

【まとめ】大学院生のための数値流体力学シリーズ

研究でCFDを使うことになった大学院生の役に立つような説明を行う記事のまとめ

mtkbirdman.com

2021.12.15