PR

【OpenFOAMv2012】upwindの解説

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

スポンサーリンク

有限体積法による離散化

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

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}\)の具体的な値を決定するのは偏微分方程式の離散化に用いられた離散化スキームであり,離散化スキームの精度や性質によって計算の安定性や結果の精度は大きく変化する

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

風上差分

風上差分(upwind difference:UD)は1次精度であり,数値粘性が非常に大きく,計算の安定性は大きく向上するが計算結果の精度は非常に悪い

計算が難しい複雑な流れ場をRANSで計算するときに,速度の対流項に二次精度以上の風上差分を設定したうえで,\(k\)や\(\varepsilon\)の対流項に風上差分を設定するのはまだ許される

風上差分では,質量流束\(F=\mathbf{u}_{f}\cdot\mathbf{S}_{f}\) の符号で上流側を判断し,上流側の値をそのまま持ってくることによって\(\phi_{f}\)を次のように補完する

\begin{equation}
\phi_{f}=\left\{\begin{array}{ll}
\phi_{_{P}} & \left(F>0\right) \\
\phi_{_{N}} & \left(F<0\right) \\
\end{array}\right.
\end{equation}

↓詳しくはこちら

OpenFOAMにおける実装

OpenFOAMでは,風上差分は以下のように実装される

\begin{align}
&\phi_{f}=w\left(\phi_{_{P}}-\phi_{_{N}}\right)+\phi_{_{N}} \tag{1} \\ \\
&w=\left\{\begin{array}{ll}
1 & \left(F>0\right) \\
0 & \left(F<0\right) \\
\end{array}\right. \tag{2}
\end{align}

参考
User Guide: Upwind divergence scheme - OpenFOAM

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

surfaceInterpolationScheme.C:離散化スキームの親クラス
upwind.H:子クラス

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

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

重み係数

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

\begin{equation}
w=\left\{\begin{array}{ll}
1 & \left(F>0\right) \\
0 & \left(F<0\right) \\
\end{array}\right. \tag{2}
\end{equation}

//- Return the interpolation weighting factors
tmp<surfaceScalarField> weights() const
{
    return pos0(this->faceFlux_);
}

≫upwind.H

pos0() は Scalar.H で次のように定義されている

//- Return 1 if s is greater_equal zero, or otherwise 0
inline Scalar pos0(const Scalar s)
{
    return (s >= 0)? 1: 0;
}

≫Scalar.H

補間式

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

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

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

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

const Field<Type>& vfi = vf; // phi
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]]);
}

≫surfaceInterpolationScheme.C

最後に書かれている sfi[fi] = Sfi[fi] & (lambda[fi]*(vfi[P[fi]] - vfi[N[fi]]) + vfi[N[fi]]) が補間式(1)に該当する

surfaceInterpolationScheme の構造

surfaceInterpolationScheme では,まず引数として\(\phi\)と\(w\)をとり,境界面での補間を返す関数である interpolate() が呼ばれていると思われる

//- 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>&
);

≫surfaceInterpolationScheme.H

interpolate() の中では,さらに dotInterpolate(geometricOneField(), vf, tlambdas) を呼んでいる

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);
}

.Hファイルのコメントから,vf が\(\phi\),tlambdasが重み係数\(w\)であると考えられる

geometricOneField() は\(1\)であり,非圧縮流れにおける密度などに使われる

Detailed Description

A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for objects which are known to be one at compile-time.

Used for example as the density argument to a function written for compressible to be used for incompressible flow.

Definition at line 55 of file geometricOneField.H.

OpenFOAM API Guide geometricOneField Class Reference

dotInterpolate() は surfaceInterpolationScheme.H で次のように定義されている

//- 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
);

dotted with given field Sf の意味は分からないが,とりあえず引数として\(1\)と\(\phi\)と\(w\)をとり,境界面での補間を返す関数である

まとめ

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

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

↓おすすめ記事

コメント