OpenFOAMを使い異方性多孔質体を解析する方法

By M.Sato

March 9, 2026

概要

前回のブログ記事において多孔質体の抵抗特性が異方性の場合を想定し、Fluentに実装する方法を説明した。今回はOpenFOAMを使い同様の計算を行う方法を説明する。

本記事で説明する内容は以下の通りである。

OpenFOAMに標準搭載されている多孔質体モデルの種類と使い方
異方性抵抗係数テンソルの座標変換による導出
DarcyForchheimerモデルおよびCodedSourceを使った9ケースの解析例
fvOptionsでの体積力指定における符号の注意点（バグ対応）

OpenFOAMにおける多孔質体モデル

流体中に置かれた多孔質体を計算に考慮する場合、流体には多孔質体による体積力が作用しているものとしてモデル化されることになる。

OpenFOAMには標準で以下の多孔質体モデルが組み込まれている。これらは基本的に直交異方性（等方性も含む）のモデルになっており主軸方向の抵抗係数を指定する。

DarcyForchheimer：流速に比例する抵抗項（粘性項）と流速の二乗に比例する項（慣性項）で定義 $$ \begin{align} S = -\left(\mu D + \frac{1}{2} \rho |v| F \right) v \notag \end{align} $$
fixedCoeff：流速の二次関数として抵抗を定義 $$ \begin{align} S = -\rho_{ref} \left(\alpha + \beta |v| \right) v \notag \end{align} $$
powerLaw：流速のべき関数として抵抗を定義 $$ \begin{align} S = -\rho C_0 |v|^{(C_1-1)/2} v \notag \end{align} $$
solidification：凝固流体を対象として、抵抗係数を温度関数として定義 $$ \begin{align} S = -\alpha \rho D(T) v \notag \end{align} $$

ここで、$S$は体積力、$v$は流速を表す。

これらの多孔質モデルを使用する場合、OpenFOAMのケースディレクトリ下のconstant/fvOptionsファイルにモデルに関するパラメータを指定することになる。

もし標準モデルでは対応できない場合、ユーザが定義するC++プログラムを使って体積力を指定することになる。

💡 標準モデルで対応可能なケース：抵抗が直交異方性であり、かつ座標系がデカルト座標系または円筒座標系であればDarcyForchheimerモデルを使って簡単に指定できる。それ以外の任意異方性・任意座標系の場合はCodedSource（C++コード記述）が必要となる。

異方性多孔質体を含む解析

ここではOpenFOAMの解析において異方性多孔質体の抵抗特性を指定する方法を説明する。 DarcyForchheimerモデルでは抵抗表現式は以下のように表される。

$$ \begin{align} S_i = -\left( \mu \sum_{j=1}^3 D_{ij} v_j + \frac{1}{2} \rho |v| \sum_{j=1}^3 F_{ij} v_j \right) \notag \end{align} $$

ここで$S_i$は$i$番目の運動量のソース項、$v_i$は$i$番目の流速、$|v|$は速度の大きさである。 $D_{ij},F_{ij}$は経験的に決められる抵抗係数テンソルである。上式をベクトル形式で表すと以下のようになる。

$$ \begin{align} \mathbf{S} &= -\left( \mu \mathbf{D} \mathbf{v} + \frac{1}{2} \rho |\mathbf{v}| \mathbf{F} \mathbf{v} \right) \notag \\ &= -\left( \mu \mathbf{D} + \frac{1}{2} \rho |\mathbf{v}| \mathbf{F} \right) \mathbf{v} \notag \\ & \equiv -\mathbf{A}\mathbf{v} \notag \end{align} $$

ここでは二次元モデルを対象にして極座標系に基づいて異方性の抵抗係数を指定する方法を説明する。まず座標系の回転によってベクトルの成分がどのように変換されるのかを示す。

図1 — 全体座標系と局所座標系

全体座標系と局所座標系でのベクトル成分の変換は次式で表される。

$$ \begin{align} \left( \begin{matrix} X \\ Y \end{matrix} \right) &= \left[ \begin{matrix} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{matrix} \right] \left( \begin{matrix} x \\ y \end{matrix} \right) = \mathbf{P} \left( \begin{matrix} x \\ y \end{matrix} \right) \notag \end{align} $$

あるいは

$$ \begin{align} \left( \begin{matrix} x \\ y \end{matrix} \right) &= \left[ \begin{matrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{matrix} \right] \left( \begin{matrix} X \\ Y \end{matrix} \right) = \mathbf{R} \left( \begin{matrix} X \\ Y \end{matrix} \right) \notag \end{align} $$

ここで座標値変換マトリックスを$\mathbf{R}$（あるいは$\mathbf{P}$）と表すことにする。 $\mathbf{R}$（あるいは$\mathbf{P}$）は直交行列であり、逆行列は転置行列になる（$\mathbf{R}^{-1}=\mathbf{R}^T$、$\mathbf{P}^{-1}=\mathbf{P}^T$）。

流速ベクトル$\mathbf{v}$、多孔質体による体積力ベクトル$\mathbf{S}$は以下の座標変換で表される。ここで添字$G$、$L$はそれぞれ全体座標系、局所座標系を表している。

$$ \mathbf{v_{L}}= \left( \begin{matrix} v_X \\ v_Y \end{matrix} \right) = \mathbf{P} \left( \begin{matrix} v_x \\ v_y \end{matrix} \right) =\mathbf{P} \mathbf{v_{G}} $$

$$ \mathbf{S_{L}}= \left( \begin{matrix} S_X \\ S_Y \end{matrix} \right) = \mathbf{P} \left( \begin{matrix} S_x \\ S_y \end{matrix} \right) =\mathbf{P} \mathbf{S_{G}} $$

次に抵抗係数であるが、局所座標系で表すと体積力と流速には次の関係式が成り立つとする（ここで抵抗係数$\mathbf{A}$は$\mathbf{D},\mathbf{F}$をまとめて表現したものと考えることができる）。局所座標系は多孔質体の主軸方向を基準にとるものと解釈することにする。$A_{12}=A_{21}=0$の場合が直交異方性に対応している。

$$ \begin{align} \mathbf{S_{L}}= \left( \begin{matrix} S_X \\ S_Y \end{matrix} \right) = \left[ \begin{matrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix} \right] \left( \begin{matrix} v_X \\ v_Y \end{matrix} \right) = \mathbf{A_{L}} \mathbf{v_{L}} \notag \end{align} $$

上式を使うと全体座標系での多孔質体の抵抗係数テンソル$\mathbf{A_G}$は以下のように変換されることがわかる。

$$ \begin{align} \mathbf{A}_{G} &= \mathbf{P}^T \mathbf{A}_{L} \mathbf{P} \notag \\ &= \mathbf{R} \mathbf{A}_{L} \mathbf{R}^T \notag \end{align} $$

多孔質体の抵抗特性が極座標系で記述される場合を考える。

図2 — 極座標系における多孔質体の抵抗特性

この場合、$\mathbf{A_G}$を具体的に記述すると以下のようになる。なおここでは簡単のため$C=\cos\theta,S=\sin\theta$とおいた。

$$ \mathbf{A}_{G} = $$ $$ \left[ \begin{matrix} C^2A_{11}+S^2A_{22}-CS(A_{12}+A_{21}) & CS(A_{11}-A_{22})+C^2A_{12}-S^2A_{21} \\ CS(A_{11}-A_{22})+C^2A_{21}-S^2A_{12} & S^2A_{11}+C^2A_{22}+CS(A_{12}+A_{21}) \end{matrix} \right] $$

解析例

OpenFOAMでは、もし多孔質体の抵抗特性が直交異方性であり、かつその空間分布がデカルト座標系あるいは円筒座標系である場合には、標準機能を使って簡単に多孔質特性を指定することが可能である（以下の例題参照）。

しかし汎用性をもたせて任意の座標方向に任意の異方性を指定したい場合にはfvOptionsファイル内に抵抗力をC++プログラム（CodedSourceタイプ）として記述する必要がある。

ここでは、以下のケースについて解析を行った。

case	model	resistance	type	numerics
1	1	isotropic（等方性）	DarcyForchheimer
2-1	1	isotropic（等方性）	CodedSource	impl
2-2	1	isotropic（等方性）	CodedSource	expl
3-1	1	orthotropic（直交異方性）	DarcyForchheimer
3-2	1	orthotropic（直交異方性）	CodedSource	expl
3-3	1	orthotropic（直交異方性）	CodedSource	impl+expl
4-1	1	orthotropic（直交異方性）	DarcyForchheimer
4-2	1	orthotropic（直交異方性）	CodedSource	impl+expl
5	2	anisotropic（完全異方性）	CodedSource	impl+expl

注：typeはfvOptionsファイルでの指定方法を表している。numericsは体積力をコード内で指定する場合の数値処理方法を表している。CodedSourceタイプの場合、自分で体積力をコーディングすることになるが、implは体積力を陰的に扱うことを意味し、explは体積力を陽的に扱うことを意味している。

❗ case-1〜4は直交異方性（主軸方向に対して独立な抵抗係数）を対象としており、case-5のみ完全異方性（非対角成分を持つ抵抗係数テンソル）を扱っている。完全異方性は標準機能では指定できないため、CodedSourceが必須となる。

以下に陰的、陽的扱いについて簡単に説明する。体積力 $S$ を、現在のステップの変数 $\phi$ を用いて以下のように線形近似して考える。

$$ \begin{align} S &= S_p \phi + S_u \notag \end{align} $$

$S_p$（陰的項）: 変数 $\phi$ に比例する係数。行列の対角成分に加えられる。
$S_u$（陽的項）: 変数に依存しない定数項。行列の右辺（ソースベクトル）に加えられる。

以下に解析に使用した2つのモデルの説明図を示す。モデル１では平行流れの中にリング状の多孔質体が置かれた場合を想定している。

モデル２は中心にある円形の境界から流入した流体がリング状の多孔質体を通過する場合を想定している。

図3 — モデル1（平行流れ中のリング状多孔質体）

図4 — モデル2（円形境界から流入する流れ）

これらのモデルは二次元場であるが、OpenFOAMでは厚み方向に1セル分を押し出して３次元メッシュとして解析を行う必要がある。以下にモデル１、モデル2のメッシュを示す。

図5 — モデル1のメッシュ

図6 — モデル2のメッシュ

case-1

このケースは等方性の多孔質体を計算した例である。これは基準となる計算であり、以下に示すcase-2と比較するために計算を行った。

計算で得られた流速ベクトル、圧力コンター図を以下に示す。流体は多孔質体を避けるように流動していることがわかる。

図7 — case-1, 2の流速ベクトル図

図8 — case-1の計算結果（流速・圧力コンター）

fvOptionsファイルは以下のように記述した。多孔質体の抵抗はDarcyForchheimerモデルを使って指定している。

porosity1
{
    type            explicitPorositySource;
    active          true;

    explicitPorositySourceCoeffs
    {
        type            DarcyForchheimer;
        selectionMode   cellZone;
        cellZone        porous;

        DarcyForchheimerCoeffs
        {
            d   d [0 -2 0 0 0 0 0] (100 100 100);
            f   f [0 -1 0 0 0 0 0] (0 0 0);

            coordinateSystem
            {
                origin  (0 0 0);
                e3      (0 0 1);
                e1      (1 0 0);
            }
        }
    }
}

case-2-1

このケースはcase-1と同じ条件であるがCodedSourceを使って多孔質体の抵抗特性を指定したものである。計算結果を見るとcase-1と全く同じ流れ場になっているため結果図は示していない（case-1参照）。

fvOptionsファイルを以下に示す。ここでは体積力をeqn.diag()項に組み込んで陰的に扱っている。

porousResistance
{
    type            vectorCodedSource;
    active          true;
    name            porousDarcySource;

    vectorCodedSourceCoeffs
    {
        selectionMode   cellZone;
        cellZone        porous;   // メッシュで定義したセルゾーン名

        fields          (U);

        codeInclude
        #{
            #include "fvMesh.H"
            #include "volFields.H"
        #};

        codeCorrect
        #{
            // 何もしない
        #};

        codeAddSup
        #{
            // 圧縮性流れの場合は codeAddSupRho を使うこと
            // ここでは非圧縮性を想定

            const Time& time = mesh_.time();
            const volVectorField& U = eqn.psi();

            // 多孔質パラメータ
            const scalar mu    = 1.0;  // 粘性係数 [Pa·s]
            const scalar rho   = 1.0;  // 密度 [kg/m³]
            const scalar alpha = 100;  // 粘性抵抗係数 [1/m²]
            const scalar beta  =   0;  // 慣性抵抗係数 [1/m]

            // セルゾーンのセルインデックスを取得
            const labelList& cells = this->cells();

            forAll(cells, i)
            {
                const label cellI = cells[i];
                const scalar& V = mesh_.V()[cellI];    // セル体積

                const scalar magU = mag(U[cellI]);
                const scalar Scoeff = mu * alpha + 0.5 * rho * magU * beta;

                // 線形化して implicit に追加するのが安定
                eqn.diag()[cellI] -= Scoeff * mesh_.V()[cellI];
            }
        #};

        codeConstrain
        #{
            // 何もしない
        #};

        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude \
            -I$(LIB_SRC)/meshTools/lnInclude
        #};

        codeLibs
        #{
            -lfiniteVolume \
            -lmeshTools
        #};
    }
}

case-2-2

このケースは基本的にcase-2-1と同じであるが体積力をeqn.source()項に組み込んで陽的に扱っている。計算結果を見るとcase-1と全く同じ流れ場になっているため結果図は示していない（case-1参照）。

fvOptionsファイルを以下に示す。

porousResistance
{
    type            vectorCodedSource;
    active          true;
    name            porousDarcySource;

    vectorCodedSourceCoeffs
    {
        selectionMode   cellZone;
        cellZone        porous;   // メッシュで定義したセルゾーン名

        fields          (U);

        codeInclude
        #{
            #include "fvMesh.H"
            #include "volFields.H"
        #};

        codeCorrect
        #{
            // 何もしない
        #};

        codeAddSup
        #{
            // 圧縮性流れの場合は codeAddSupRho を使うこと
            // ここでは非圧縮性を想定

            const Time& time = mesh_.time();
            const volVectorField& U = eqn.psi();

            // 多孔質パラメータ
            const scalar mu    = 1.0;  // 粘性係数 [Pa·s]
            const scalar rho   = 1.0;  // 密度 [kg/m³]
            const scalar alpha = 100;  // 粘性抵抗係数 [1/m²]
            const scalar beta  =   0;  // 慣性抵抗係数 [1/m]

            // セルゾーンのセルインデックスを取得
            const labelList& cells = this->cells();

            forAll(cells, i)
            {
                const label cellI = cells[i];
                const scalar& V = mesh_.V()[cellI];    // セル体積

                const scalar magU = mag(U[cellI]);
                const scalar Scoeff = mu * alpha + 0.5 * rho * magU * beta;

                // 陽的に追加
                eqn.source()[cellI] += Scoeff * mesh_.V()[cellI] * U[cellI];
            }
        #};

        codeConstrain
        #{
            // 何もしない
        #};

        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude \
            -I$(LIB_SRC)/meshTools/lnInclude
        #};

        codeLibs
        #{
            -lfiniteVolume \
            -lmeshTools
        #};
    }
}

case-3-1

💡 case-3以降は直交異方性の多孔質体を扱う。円筒座標系で異方性を指定する場合、coordinateSystemのtypeをcylindricalに設定することがポイントである。

このケースは主軸#2方向（円周方向）の抵抗を大きくしたケースであり、その結果主軸#1方向（半径方向）の流れが顕著になる。計算結果を見ると多孔質体リングの対称面付近で流速が大きくなっていることがわかる。

図9 — case-3の流速ベクトル図

図10 — case-3の計算結果（流速・圧力コンター）

fvOptionsファイルにはDarcyForchheimerモデルを指定した。このケースでは異方性は円筒座標系に基づいて指定されるため、coordinateSystemのtypeをcylindricalにしている。

fvOptionsファイルを以下に示す（粘性抵抗項dで指定する3成分は$D_{rr}, D_{\theta\theta}, D_{zz}$の順になっている）。

porosity1
{
    type            explicitPorositySource;
    active          true;

    explicitPorositySourceCoeffs
    {
        type            DarcyForchheimer;
        selectionMode   cellZone;
        cellZone        porous;

        DarcyForchheimerCoeffs
        {
            d   d [0 -2 0 0 0 0 0] (1 100 0);
            f   f [0 -1 0 0 0 0 0] (0 0 0);

            coordinateSystem
            {
                type    cylindrical;
                origin  (0 0 0);
                e3      (0 0 1);
                e1      (1 0 0);
            }
        }
    }
}

case-3-2

このケースはcase-3-1と同じ内容であるが、FluentにおけるUDF記述と同様に流体力成分を、抵抗力テンソルを使って、素朴な方法で、成分ごとに個別に記述したものである。

なお、計算結果はcase-3-1と全く同じ流れ場になっているため結果図は示していない（case-3-1参照）。

fvOptionsファイルを以下に示す。

porousResistance
{
    type            vectorCodedSource;
    active          true;
    name            porousDarcySource;

    vectorCodedSourceCoeffs
    {
        selectionMode   cellZone;
        cellZone        porous;   // メッシュで定義したセルゾーン名

        fields          (U);

        codeInclude
        #{
            #include "fvMesh.H"
            #include "volFields.H"
        #};

        codeCorrect
        #{
            // 何もしない
        #};

        codeAddSup
        #{
            // 圧縮性流れの場合は codeAddSupRho を使うこと
            // ここでは非圧縮性を想定

            const Time& time = mesh_.time();
            const volVectorField& U = eqn.psi();

            // 多孔質パラメータ
            const scalar mu    = 1.0;  // 粘性係数 [Pa·s]
            const scalar rho   = 1.0;  // 密度 [kg/m³]

            // セルゾーンのセルインデックスを取得
            const labelList& cells = this->cells();

            // ソース項を eqn に追加（負号 = 抵抗）
            forAll(cells, i)
            {
                const label cellI = cells[i];
                const scalar& V = mesh_.V()[cellI];    // セル体積
                const vector& C = mesh_.C()[cellI];    // セル重心座標

                const scalar x = C[vector::X];
                const scalar y = C[vector::Y];

                const scalar rr = sqrt(x*x+y*y);
                const scalar cs = x/rr;
                const scalar sn = y/rr;

                // 粘性抵抗係数テンソル [1/m²]
                const scalar A11 =   1.0;
                const scalar A12 =   0.0;
                const scalar A21 =   0.0;
                const scalar A22 = 100.0;

                const scalar d11 = cs*cs*A11 + sn*sn*A22 - cs*sn*(A12+A21);
                const scalar d12 = cs*sn*(A11-A22) + cs*cs*A12 - sn*sn*A21;
                const scalar d21 = cs*sn*(A11-A22) + cs*cs*A21 - sn*sn*A12;
                const scalar d22 = sn*sn*A11 + cs*cs*A22 + cs*sn*(A12+A21);

                // 各方向の速度成分
                const scalar Ux = U[cellI][vector::X];
                const scalar Uy = U[cellI][vector::Y];

                const scalar srcx = mu*( +d11*Ux + d12*Uy );
                const scalar srcy = mu*( +d21*Ux + d22*Uy );

                eqn.source()[cellI][vector::X] += V * srcx;
                eqn.source()[cellI][vector::Y] += V * srcy;
                eqn.source()[cellI][vector::Z] += 0;
            }
        #};

        codeConstrain
        #{
            // 何もしない
        #};

        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude \
            -I$(LIB_SRC)/meshTools/lnInclude
        #};

        codeLibs
        #{
            -lfiniteVolume \
            -lmeshTools
        #};
    }
}

case-3-3

このケースはcase-3-2と同じ内容であるが、個別のテンソル成分をコードに記述するのではなく、全体の抵抗係数テンソルをコード内で座標変換して計算するようにしたものである。なお体積力は、陰的成分（eqn.diag()）と陽的成分（eqn.source()）の両方をコード内で指定している。

💡 impl+explの組み合わせが最も汎用的な記述方法である。等方成分をeqn.diag()（陰的・安定化）で、異方性補正分をeqn.source()（陽的）で指定することで、収束性と精度を両立できる。

なお、計算結果はcase-3-1と全く同じ流れ場になっているため結果図は示していない（case-3-1参照）。

fvOptionsファイルを以下に示す。

porousResistance
{
    type            vectorCodedSource;
    active          true;
    name            porousDarcySource;

    vectorCodedSourceCoeffs
    {
        selectionMode   cellZone;
        cellZone        porous;   // メッシュで定義したセルゾーン名

        fields          (U);

        codeInclude
        #{
            #include "fvMesh.H"
            #include "volFields.H"
        #};

        codeCorrect
        #{
            // 何もしない
        #};

        codeAddSup
        #{
            // 圧縮性流れの場合は codeAddSupRho を使うこと
            // ここでは非圧縮性を想定

            const Time& time = mesh_.time();
            const volVectorField& U = eqn.psi();

            // 多孔質パラメータ
            const scalar mu    = 1.0;  // 粘性係数 [Pa·s]
            const scalar rho   = 1.0;  // 密度 [kg/m³]

            // 抵抗テンソル（対角異方性の例）
            const tensor alpha
            (
                1,     0,     0,
                0,     100,   0,
                0,     0,     1
            );
            const tensor beta
            (
                0,     0,     0,
                0,     0,     0,
                0,     0,     0
            );

            // セルゾーンのセルインデックスを取得
            const labelList& cells = this->cells();

            // ソース項を eqn に追加（負号 = 抵抗）
            forAll(cells, i)
            {
                const label cellI = cells[i];

                const scalar& V = mesh_.V()[cellI];    // セル体積
                const vector& C = mesh_.C()[cellI];    // セル重心座標

                const scalar x = C[vector::X];
                const scalar y = C[vector::Y];

                const scalar rr = sqrt(x*x+y*y);
                const scalar cs = x/rr;
                const scalar sn = y/rr;

                // z軸周り回転行列
                const tensor R(cs, -sn, 0,   sn, cs, 0,   0, 0, 1);

                // 回転済みテンソル
                const tensor alphaR = R & alpha & R.T();
                const tensor betaR  = R & beta  & R.T();

                const scalar magU = mag(U[cellI]);
                // fixedCoeff と同じ構造
                const tensor Cd    = rho*(alphaR + betaR*magU);
                const scalar isoCd = tr(Cd);

                // implicit部（等方・安定化）
                eqn.diag()[cellI]   -= V * isoCd;

                // explicit部（異方性補正）
                eqn.source()[cellI] += V * ((Cd - I*isoCd) & U[cellI]);
            }
        #};

        codeConstrain
        #{
            // 何もしない
        #};

        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude \
            -I$(LIB_SRC)/meshTools/lnInclude
        #};

        codeLibs
        #{
            -lfiniteVolume \
            -lmeshTools
        #};
    }
}

case-4-1

このケースは主軸#1方向（半径方向）の抵抗を大きくしたケースであり、その結果主軸#2方向（円周方向）の流れが顕著になる。計算結果を見ると多孔質体リングの対称面付近で流速が遅くなり、側面付近で流速が速くなっていることがわかる。

図11 — case-4の流速ベクトル図

図12 — case-4の計算結果（流速・圧力コンター）

fvOptionsファイルにはDarcyForchheimerモデルを指定した。このケースでは異方性は円筒座標系に基づいて指定されるため、coordinateSystemのtypeをcylindricalにしている。fvOptionsファイルを以下に示す。

porosity1
{
    type            explicitPorositySource;
    active          true;

    explicitPorositySourceCoeffs
    {
        type            DarcyForchheimer;
        selectionMode   cellZone;
        cellZone        porous;

        DarcyForchheimerCoeffs
        {
            d   d [0 -2 0 0 0 0 0] (100 1 0);
            f   f [0 -1 0 0 0 0 0] (0 0 0);

            coordinateSystem
            {
                type    cylindrical;
                origin  (0 0 0);
                e3      (0 0 1);
                e1      (1 0 0);
            }
        }
    }
}

case-4-2

このケースは基本的にcase-4-1と同じであるが、CodedSourceを使って多孔質体の抵抗特性を指定したものである。なお、計算結果はcase-4-1と全く同じ流れ場になっているため結果図は示していない（case-4-1参照）。

fvOptionsファイルを以下に示す。

porousResistance
{
    type            vectorCodedSource;
    active          true;
    name            porousDarcySource;

    vectorCodedSourceCoeffs
    {
        selectionMode   cellZone;
        cellZone        porous;   // メッシュで定義したセルゾーン名

        fields          (U);

        codeInclude
        #{
            #include "fvMesh.H"
            #include "volFields.H"
        #};

        codeCorrect
        #{
            // 何もしない
        #};

        codeAddSup
        #{
            // 圧縮性流れの場合は codeAddSupRho を使うこと
            // ここでは非圧縮性を想定

            const Time& time = mesh_.time();
            const volVectorField& U = eqn.psi();

            // 多孔質パラメータ
            const scalar mu    = 1.0;  // 粘性係数 [Pa·s]
            const scalar rho   = 1.0;  // 密度 [kg/m³]

            // 抵抗テンソル（対角異方性の例）
            const tensor alpha
            (
                100,   0,     0,
                0,     1,     0,
                0,     0,     1
            );
            const tensor beta
            (
                0,     0,     0,
                0,     0,     0,
                0,     0,     0
            );

            // セルゾーンのセルインデックスを取得
            const labelList& cells = this->cells();

            // ソース項を eqn に追加（負号 = 抵抗）
            forAll(cells, i)
            {
                const label cellI = cells[i];

                const scalar& V = mesh_.V()[cellI];    // セル体積
                const vector& C = mesh_.C()[cellI];    // セル重心座標

                const scalar x = C[vector::X];
                const scalar y = C[vector::Y];

                const scalar rr = sqrt(x*x+y*y);
                const scalar cs = x/rr;
                const scalar sn = y/rr;

                // z軸周り回転行列
                const tensor R(cs, -sn, 0,   sn, cs, 0,   0, 0, 1);

                // 回転済みテンソル
                const tensor alphaR = R & alpha & R.T();
                const tensor betaR  = R & beta  & R.T();

                const scalar magU = mag(U[cellI]);
                // fixedCoeff と同じ構造
                const tensor Cd    = rho*(alphaR + betaR*magU);
                const scalar isoCd = tr(Cd);

                // implicit部（等方・安定化）
                eqn.diag()[cellI]   -= V * isoCd;

                // explicit部（異方性補正）
                eqn.source()[cellI] += V * ((Cd - I*isoCd) & U[cellI]);
            }
        #};

        codeConstrain
        #{
            // 何もしない
        #};

        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude \
            -I$(LIB_SRC)/meshTools/lnInclude
        #};

        codeLibs
        #{
            -lfiniteVolume \
            -lmeshTools
        #};
    }
}

case-5

❗ このケースは完全異方性（非対称テンソル）を扱う唯一の例であり、DarcyForchheimer標準モデルでは再現不可能な流れ場を生成する。CodedSourceによる実装の真価が発揮されるケースである。

このケースは完全異方性の多孔質特性を指定した例であり、流速方向と抵抗方向が異なっている。ここでは流速方向に対して右側に体積力が作用するとした。計算結果を見ると多孔質体に沿って時計まわりに循環するような流れが生じていることがわかる。

図13 — case-5の流速ベクトル図

図14 — case-5の計算結果（流速・圧力コンター）

fvOptionsファイルを以下に示す。

porousResistance
{
    type            vectorCodedSource;
    active          true;
    name            porousDarcySource;

    vectorCodedSourceCoeffs
    {
        selectionMode   cellZone;
        cellZone        porous;   // メッシュで定義したセルゾーン名

        fields          (U);

        codeInclude
        #{
            #include "fvMesh.H"
            #include "volFields.H"
        #};

        codeCorrect
        #{
            // 何もしない
        #};

        codeAddSup
        #{
            // 圧縮性流れの場合は codeAddSupRho を使うこと
            // ここでは非圧縮性を想定

            const Time& time = mesh_.time();
            const volVectorField& U = eqn.psi();

            // 多孔質パラメータ
            const scalar mu    = 1.0;  // 粘性係数 [Pa·s]
            const scalar rho   = 1.0;  // 密度 [kg/m³]

            // 抵抗テンソル（完全異方性の例）
            const tensor alpha
            (
                1,     -100,  0,
                100,   1,     0,
                0,     0,     1
            );
            const tensor beta
            (
                0,     0,     0,
                0,     0,     0,
                0,     0,     0
            );

            // セルゾーンのセルインデックスを取得
            const labelList& cells = this->cells();

            // ソース項を eqn に追加（負号 = 抵抗）
            forAll(cells, i)
            {
                const label cellI = cells[i];
                const scalar& V = mesh_.V()[cellI];    // セル体積
                const vector& C = mesh_.C()[cellI];    // セル重心座標

                const scalar x = C[vector::X];
                const scalar y = C[vector::Y];

                const scalar rr = sqrt(x*x+y*y);
                const scalar cs = x/rr;
                const scalar sn = y/rr;

                // z軸周り回転行列
                const tensor R(cs, -sn, 0,   sn, cs, 0,   0, 0, 1);

                // 回転済みテンソル
                const tensor alphaR = R & alpha & R.T();
                const tensor betaR  = R & beta  & R.T();

                const scalar magU = mag(U[cellI]);
                // fixedCoeff と同じ構造
                const tensor Cd    = rho*(alphaR + betaR*magU);
                const scalar isoCd = tr(Cd);

                // implicit部（等方・安定化）
                eqn.diag()[cellI]   -= V * isoCd;

                // explicit部（異方性補正）
                eqn.source()[cellI] += V * ((Cd - I*isoCd) & U[cellI]);
            }
        #};

        codeConstrain
        #{
            // 何もしない
        #};

        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude \
            -I$(LIB_SRC)/meshTools/lnInclude
        #};

        codeLibs
        #{
            -lfiniteVolume \
            -lmeshTools
        #};
    }
}

OpenFOAMのバグ？

⚠️ CodedSourceで体積力を指定する際、符号を直感と逆にしないと収束しないという問題がある。現バージョンのOpenFOAMを使用する場合は、以下の説明をよく読んで符号の扱いに注意してほしい。この問題はOpenFOAM開発者に報告済みである。

CodedSourceを使って体積力を指定する例では、奇妙な記述になっていることを気づいた読者がいるかも知れない。

例えばcase-5では、連立方程式に以下のような記述で体積力を追加している。

                // implicit部（等方・安定化）
                eqn.diag()[cellI]   -= V * isoCd;
                // explicit部（異方性補正）
                eqn.source()[cellI] += V * ((Cd - I*isoCd) & U[cellI]);

本来ならば以下の様に記述するべき（eqn.diag()は連立方程式の係数行列の対角成分を対角優位にする）であるが、収束しないため符号を逆にしている。実際このようにすることで正解に収束することが確認できている。

                // implicit部（等方・安定化）
                eqn.diag()[cellI]   += V * isoCd;
                // explicit部（異方性補正）
                eqn.source()[cellI] -= V * ((Cd - I*isoCd) & U[cellI]);

これについては、OpenFOAMの不具合として開発者に連絡している。

現バージョンのOpenFOAMでfvOptionsファイルで体積力を指定する場合には、上記のように符号を逆転する必要があることに注意してほしい。

まとめ

✅ OpenFOAMを使って異方性多孔質体を通過する流れを計算する設定例を示した。

✅ デカルト座標系、円筒座標系の場合、直交異方性は標準機能を使って簡単に指定できることを示した。

✅ コード内で体積力を指定する方法を使って、任意の座標系に基づき任意の抵抗力成分（テンソル）を指定できることを示した。

弊社の解析事例

弊社の流体解析事例については、下記のリンクからご覧ください。

👉 解析事例のご紹介 OpenFOAM

👉 流体解析実績及び事例