境界条件の設定

2021年6月20日

はじめに

境界条件の設定について。

使用バージョン

OpenFOAM v2012、v8

境界タイプの設定

OpenFOAM のメッシュファイル群 (constant/polyMesh) には boundary ファイルという境界情報を含むファイルがある。これは以下のような内容になっている。

5
(
    inlet
    {
        type            patch;
        nFaces          30;
        startFace       24170;
    }
    outlet
    {
        type            patch;
        nFaces          57;
        startFace       24200;
    }
    upperWall
    {
        type            wall;
        nFaces          223;
        startFace       24257;
    }
    lowerWall
    {
        type            wall;
        nFaces          250;
        startFace       24480;
    }
    frontAndBack
    {
        type            empty;
        inGroups        1(empty);
        nFaces          24450;
        startFace       24730;
    }
)

この情報はふつうメッシュ作成ソフトあるいはメッシュ変換ユーティリティが書き出すが、境界条件に合わせて境界のタイプを設定する必要がある。境界のタイプにはつぎのものがある。

  • patch (パッチ)
  • wall (壁)
  • symmetryPlane (対称面)
  • symmetry (対称面)
  • cyclic (周期境界)
  • cyclicAMI (不整合周期境界)
  • wedge (2 次元軸対称)
  • empty (2 次元)

パッチ (patch) というのは、入口や出口で使用されるタイプである。wall はそのまま壁条件であり、壁関数境界条件を設定する場合はこれを選択する必要がある。

symmetryPlane, symmetry は対称条件を与える面のタイプで、前者が平面、後者がそれ以外に適用されるものである (平面以外の対称面というのがよくわからんが)。基本的には前者を用いればよい。

cyclic, cyclicAMI については、対応する 2 面で設定する必要があり、それぞれ "neighbourPatch" で対応する境界を指定する。

2 次元問題の場合は、計算しない方向の両面に empty を指定する。2 次元軸対称については、実際の 3 次元モデルの 5°以下の楔形のモデルを作り、計算しない両面を wedge に設定すればよい。

v1706 以降では、オーバーセットメッシュ用の overset タイプがある。

境界条件の設定

境界条件の設定は 0 ディレクトリ内の各フィールドファイルで行う。必要なフィールドファイルはソルバーごとに異なる。乱流モデルに k-εモデルを用いる場合、simpleFoam、buoyantBoussinesqSimpleFoam、buoyantSimpleFoam で必要になるフィールドファイルはそれぞれつぎのようになる。

ソルバー必要なフィールドファイル
simpleFoam
U, p, k, epsilon, nut
buoyantBoussinesqSimpleFoam
U, p, p_rgh, T, k, epsilon, nut, alphat
buoyantSimpleFoam
U, p, p_rgh, T, k, epsilon, nut, alphat

v2.4.0 以前

ソルバー必要なフィールドファイル
simpleFoam
U, p, k, epsilon,
buoyantBoussinesqSimpleFoam
U, p, p_rgh, T, k, epsilon
buoyantSimpleFoam
U, p, p_rgh, T, k, epsilon

※v2.4.0 以前では、nut, alphat はあってもよい (なくても自動で作られる)。ただし、buoyantSimpleFoam では nut ではなく mut。

ただし、圧力の単位に注意する。simpleFoam と buoyantBoussinesqSimpleFoam の圧力は密度で割ったものであるが、buoyantSimpleFoam の圧力は本来の圧力である。

基本設定

U についてフィールドファイルの基本設定を示す。

U

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (1 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    wall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

ヘッダー

ヘッダーは「おまじない」と思えばよい。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}

フィールドデータの場合、class には以下のようなフィールドのクラスを書くようになっている。

  • volScalarField : スカラー場 (圧力や温度など)
  • volVectorField : ベクトル場 (速度など)
  • volSymmTensorField : 対称テンソル場 (レイノルズ応力など)

dimensions

dimensions で単位を設定する。次のような形式で各単位の指数を指定する。

[kg m s K mol A Cd]

"[0 1 -1 0 0 0 0]" は m/s を表している。つぎのような表現も可能である。

dimensions      [m/s];

internalField

フィールド内部の値を設定する。0 ディレクトリの場合、初期値を設定する。

boundaryField

境界条件を設定する。境界条件はつぎのような形式で設定する。

    inlet
    {
        type            fixedValue;
        value           uniform (10 0 0);
    }

メッシュの boundary に記述された境界名に対して、境界条件のタイプと設定を設定する。基本的な境界条件タイプを以下に挙げる。

fixedValue
値の固定。
zeroGradient
勾配を 0 に固定。値の指定は必要ない。

特別な場合を除いて、基本的には type と value が 2 つとも必要である。value は境界条件の設定に使われる場合も使われない場合もあり、初期値として使われる場合もある。

境界名では、二重引用符で囲んで "..." のような書き方をすると、正規表現を使うことができる (v1.6 以降)。"inlet.*" と書くと "inlet-1" や "inlet-2" などに当てはまり、"(inlet-1|inlet-2)" と書くと "inlet-1" か "inlet-2" に当てはまる。".*" と書くとすべての境界に当てはまる。正規表現と直接指定が同時にある場合は、順番によらず直接指定が優先されるため、".*" で基本的な設定をしておき、それ以外は個別に設定するような書き方ができる。たとえば、次のような内容のファイルをテンプレートして用意しておけばよいかもしれない。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    ".*"
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

#includeEtc 文 (v2.4.0 以降) では、symmetryPlane や empty などの拘束型 (constraint) の境界条件の設定を読み込み、これらの設定を記述しなくてもよいようにしている。これは次のように書いたものと同じである。

#include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"

代表的なフィールドの基本設定

代表的なフィールドの基本設定を以下に示す。"基本設定" の中身は、type と value を表しており、追加のパラメタがあればそれを示している。

非圧縮性共通

フィールドクラス単位初期値基本設定 (壁の設定)
UvolVectorField[0 1 -1 0 0 0 0] (m/s)(0 0 0)fixedValue, (0 0 0)
pvolScalarField[0 2 -2 0 0 0 0] (m2/s2)0zeroGradient
TvolScalarField[0 0 0 1 0 0 0] (K)300zeroGradient
kvolScalarField[0 2 -2 0 0 0 0] (m2/s2)1kqRWallFunction, $internalField
epsilonvolScalarField[0 2 -3 0 0 0 0] (m2/s3)1epsilonWallFunction, $internalField
omegavolScalarField[0 0 -1 0 0 0 0] (1/s)1omegaWallFunction, $internalField
RvolSymmTensorField[0 2 -2 0 0 0 0] (m2/s2)(1 0 0 1 0 1)kqRWallFunction, $internalField
nutvolScalarField[0 2 -1 0 0 0 0] (m2/s)0nutkWallFunction, $internalField

※v4 以降は U の静止壁条件 (non-slip 条件) に noSlip が使える。

v2.4.0 以前

フィールドクラス単位初期値基本設定 (壁の設定)
UvolVectorField[0 1 -1 0 0 0 0] (m/s)(0 0 0)fixedValue, (0 0 0)
pvolScalarField[0 2 -2 0 0 0 0] (m2/s2)0zeroGradient
TvolScalarField[0 0 0 1 0 0 0] (K)300zeroGradient
kvolScalarField[0 2 -2 0 0 0 0] (m2/s2)1kqRWallFunction, $internalField
epsilonvolScalarField[0 2 -3 0 0 0 0] (m2/s3)1epsilonWallFunction, $internalField
omegavolScalarField[0 0 -1 0 0 0 0] (1/s)1omegaWallFunction, $internalField
RvolSymmTensorField[0 2 -2 0 0 0 0] (m2/s2)(1 0 0 1 0 1)kqRWallFunction, $internalField
nuSgsvolScalarField[0 2 -1 0 0 0 0] (m2/s)0zeroGradient

圧縮性共通

フィールドクラス単位初期値基本設定 (壁の設定)
UvolVectorField[0 1 -1 0 0 0 0] (m/s)(0 0 0)fixedValue, (0 0 0)
pvolScalarField[1 -1 -2 0 0 0 0] (kg/m-s2)101325zeroGradient
TvolScalarField[0 0 0 1 0 0 0] (K)300zeroGradient
kvolScalarField[0 2 -2 0 0 0 0] (m2/s2)1kqRWallFunction, $internalField
epsilonvolScalarField[0 2 -3 0 0 0 0] (m2/s3)1epsilonWallFunction, $internalField
omegavolScalarField[0 0 -1 0 0 0 0] (1/s)1omegaWallFunction, $internalField
RvolSymmTensorField[0 2 -2 0 0 0 0] (m2/s2)(1 0 0 1 0 1)kqRWallFunction, $internalField
nutvolScalarField[0 2 -1 0 0 0 0] (m2/s)0nutkWallFunction, $internalField
alphatvolScalarField[1 -1 -1 0 0 0 0] (kg/m-s)0zeroGradient

v2.4.0 以前

フィールドクラス単位初期値基本設定 (壁の設定)
UvolVectorField[0 1 -1 0 0 0 0] (m/s)(0 0 0)fixedValue, (0 0 0)
pvolScalarField[1 -1 -2 0 0 0 0] (kg/m-s2)101325zeroGradient
TvolScalarField[0 0 0 1 0 0 0] (K)300zeroGradient
kvolScalarField[0 2 -2 0 0 0 0] (m2/s2)1compressible::kqRWallFunction, $internalField
epsilonvolScalarField[0 2 -3 0 0 0 0] (m2/s3)1compressible::epsilonWallFunction, $internalField
omegavolScalarField[0 0 -1 0 0 0 0] (1/s)1compressible::omegaWallFunction, $internalField
RvolSymmTensorField[0 2 -2 0 0 0 0] (m2/s2)(1 0 0 1 0 1)compressible::kqRWallFunction, $internalField
muSgsvolScalarField[1 -1 -1 0 0 0 0] (kg/m-s)0zeroGradient
alphaSgsvolScalarField[1 -1 -1 0 0 0 0] (kg/m-s)0zeroGradient

buoyantBoussinesqSimpleFoam

フィールドクラス単位初期値基本設定 (壁の設定)
pvolScalarField[0 2 -2 0 0 0 0] (m2/s2)0calculated, $internalField
p_rghvolScalarField[0 2 -2 0 0 0 0] (m2/s2)0fixedFluxPressure, $internalField, rho rhok
alphatvolScalarField[0 2 -1 0 0 0 0] (m2/s)0alphatJayatillekeWallFunction, $internalField, Prt 0.85

buoyantSimpleFoam

フィールドクラス単位初期値基本設定 (壁の設定)
pvolScalarField[1 -1 -2 0 0 0 0] (kg/m-s2)101325calculated, $internalField
p_rghvolScalarField[1 -1 -2 0 0 0 0] (kg/m-s2)101325fixedFluxPressure, $internalField
alphatvolScalarField[1 -1 -1 0 0 0 0] (kg/m-s)1compressible::alphatJayatillekeWallFunction, $internalField, Prt 0.85

ここで、R (レイノルズ応力) の volSymmTensorField は対称テンソル場で、対称テンソルは 6 成分からなり、各成分を (xx xy xz yy yz zz) の順番で指定する。

0 ディレクトリのバックアップ

0 ディレクトリの中身はたまに標準ユーティリティなどにより書き換えられてしまうことがあるため、基本的な設定をしたあとに、"0.org" などといった名前でオリジナルをコピーしておくとよい。

$ cp -r 0 0.org

入口・出口条件

単純な条件

単純な入口・出口条件は、つぎのようになる。

入口

フィールド設定
UfixedValue, (ux, uy, uz)
p/p_rghzeroGradient
その他のスカラー値fixedValue, 値指定

出口

フィールド設定
UzeroGradient/advective
p/p_rghfixedValue, 値指定 (0 か 101325)
その他のスカラー値zeroGradient

逆流を考慮した出口条件

つぎのような設定により、逆流を考慮した出口条件を設定できる。ただし、逆流は計算上好ましくないので、極力逆流が起こらないように出口を設定したほうがよい。

フィールド設定
UpressureInletOutletVelocity, (0 0 0)
p/p_rghfixedValue, 値指定 (0 か 101325)、totalPressure、prghPressure
その他のスカラー値inletOutlet, inletValue 値指定

pressureInletOutletVelocity は圧力から入出力速度を決めるもので、つぎのように設定する。

    outlet
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }

totalPressure は全圧を指定するもので、つぎのように設定する。

    outlet
    {
        type            totalPressure;
        value           uniform 0;
        gamma           1.4;
        p0              uniform 0;
    }

p0 で圧力を指定、gamma は比熱比である。

浮力を考慮した圧力を指定したい場合は、prghPressure を用いる。

    outlet
    {
        type            prghPressure;
        p               uniform 101325;
        value           uniform 101325;
    }

浮力で流れが上部境界から抜けていく場合などは、この境界条件を用いる。

inletOutlet はつぎのように設定する。

    outlet
    {
        type            inletOutlet;
        value           uniform 300;
        inletValue      uniform 300;
    }

流速の条件

流入速度を指定するために fixedValue では不便な場合がある。以下に fixedValue 以外の流入境界条件を挙げる。

面の法線方向の流速の指定

    inlet
    {
        type            surfaceNormalFixedValue;
        refValue        uniform 10; // [m/s]
    }

refValue で面の法線方向の流速の大きさを指定。

メッシュによっては面の方向が領域の内側ではなく外側に向いている場合があり、注意が必要。

体積流量の指定

    inlet
    {
        type            flowRateInletVelocity;
        volumetricFlowRate  constant 2.54e-4; // [m3/s]
        value           uniform (0 0 0);
    }

質量流量の指定

    inlet
    {
        type            flowRateInletVelocity;
        massFlowRate    constant 2.54e-4; // [kg/s]
        rhoInlet        1; // [kg/m3]
        value           uniform (0 0 0);
    }

非圧縮性流体ソルバーの場合、rhoInlet で流入する流体の密度を指定する。圧縮性ソルバーの場合は rho の値が用いられる。

乱流流入条件

   inlet
   {
       type             turbulentInlet;
       referenceField   uniform (10 0 0); // 流速 [m/s]
       fluctuationScale (0.02 0.01 0.01); // 変動スケール
   }

変動を含んだ流入条件を設定する。

乱流諸量の条件

入口の $k$ や $\varepsilon$、$\omega$ などにどのような値を指定すべきかは、一般的にははっきりしない。実験値が得られない場合、適当に見積もる必要がある。$k$、$\varepsilon$、$\omega$ はそれぞれ次式で表される。

$$k = \frac{3}{2} (UI)^2$$ $$\varepsilon = \frac{C_{\mu}^{3/4} k^{3/2}}{L_m}$$ $$\omega = \frac{\varepsilon}{C_\mu k} = \frac{k^{1/2}}{C_{\mu}^{1/4} L_m}$$ $$C_{\mu} = 0.09$$

※$\omega$ の $C_{\mu}$ は、乱流粘性が $k$-$\varepsilon$ と $k$-$\omega$ で一致するとしたときに出てくる。

$I$ は乱流強度 (turbulent intensity)、$L_m$ は混合長 (mixing length) である。乱流強度は、十分に発達した流れでは数%の値をとる。混合長は、ダクト内の十分発達した流れでは次式で近似できる。

$$L_m = 0.07 L$$

$L$ は代表長さで、円形のダクトであれば直径、そうでなければ水力直径 (4×断面積/断面周長) とする。

レイノルズ応力 $R$ については、対角成分に $(3/2)k$ を入れておけばよい。

$k$ や$\varepsilon$ などを直接指定するのではなく、乱流強度や混合長で指定するための特別な境界条件がある。

$k$

    inlet
    {
        type            turbulentIntensityKineticEnergyInlet;
        intensity       0.05; // 5 %
        value           uniform 1;
    }

intensity に乱流強度を設定する。乱流強度が 5 % ならば 0.05 を指定する。

$\varepsilon$

    inlet
    {
        type            turbulentMixingLengthDissipationRateInlet;
        mixingLength    0.005; // [m]
        value           uniform 1;
    }

mixingLength に混合長を設定する。

$\omega$

    inlet
    {
        type            turbulentMixingLengthFrequencyInlet;
        mixingLength    0.005; // [m]
        k               k;
        value           uniform 1;
    }

mixingLength に混合長を設定する。

v2.4.x 以前では、buoyantSimpleFoam などの圧縮性流体ソルバーの場合、εとωについては、タイプ名の頭に "compressible::" を付ける (たとえば "compressible::turbulentMixingLengthDissipationRateInlet")。

壁面速度の設定

壁面速度の設定は fixedValue でそのまま指定すればよい。いくつか特殊な設定について以下で述べる。

スリップ

壁面をスリップ条件にしたい場合は、U で slip を指定する。

    wall
    {
        type            slip;
    }

それ以外の変数については zeroGradient にする。

回転速度の設定

円筒の壁が回転するような場合は、つぎのようにして回転速度を設定する。

    wall
    {
        type            rotatingWallVelocity;
        orgin           (0 0 0);     // 回転軸の位置
        axis            (0 0 1);     // 回転軸ベクトル
        omega           constant 10; // 角速度 [rad/s]
    }

速度 v を指定したい場合、円筒半径を r とすると、omega = v/r を指定すればよい。

熱の境界条件の設定

熱 (温度) の境界条件としては、断熱である zeroGradient、温度指定である fixedValue 以外に以下のようなものがある。

熱流束の設定

壁面に熱流束を設定するには externalWallHeatFluxTemperature を用いる。

    wall
    {
        type            externalWallHeatFluxTemperature;
        mode            flux;
        kappaMethod     fluidThermo;
        q               uniform 100; // 熱流束 [W/m2]
        value           uniform 300; // 初期温度 [K]
    }

mode で "power" を指定すれば、熱流束 [W/m2] の代わりに熱量 [W] を指定できる。

    wall
    {
        type            externalWallHeatFluxTemperature;
        mode            power;
        kappaMethod     fluidThermo;
        Q               100; // 熱量 [W]
        value           uniform 300;
    }

externalWallHeatFluxTemperature は buoyantBoussinesqSimpleFoam などでは利用できない。Boussinesq 近似と併用したい場合は、buoyantSimpleFoam などで状態方程式に Boussinesq を用いる (v4 以上)。

turbulentHeatFluxTemperature

v4 以前であれば、turbulentHeatFluxTemperature という境界条件も使える。

buoyantSimpleFoam などではつぎのように設定する。

    wall
    {
        type            compressible::turbulentHeatFluxTemperature;
        heatSource      flux;        // flux (熱流束) または power (熱量)
        q               uniform 100; // 熱流束 [W/m2] または 熱量 [W]
        kappaMethod     fluidThermo;
        kappa           none;
        value           uniform 300;
    }

v3.0 以前だと、kappaMethod, kappa の部分は以下のように書く。

        kappa           fluidThermo;
        kappaName       none;

buoyantBoussinesqSimpleFoam などではつぎのように設定する。

    wall
    {
        type            turbulentHeatFluxTemperature;
        heatSource      flux;
        q               uniform 100;
        alphaEff        alphaEff;
        value           uniform 300;
    }

ただし、constant/transportProperties で Cp0 (密度×比熱) を設定する必要がある。

壁面熱伝達の設定

壁面熱伝達を設定するには、externalWallHeatFluxTemperature の mode で "coefficient" を指定する。

    wall
    {
        type            externalWallHeatFluxTemperature;
        mode            coefficient;
        kappaMethod     fluidThermo;
        Ta              constant 500; // 外部温度 [K]
        h               constant 10;  // 熱伝達係数 [W/m2-K]
        value           uniform 300;
    }

外部の輻射の影響を考慮したい場合は、emissivity で放射率を指定する。

        emissivity      0.8; // 放射率

また、熱伝達率の異なる板が複数枚積層している場合を表すこともできる。

        thicknessLayers (0.1 0.2 0.3 0.4); // 厚み [m]
        kappaLayers     (1 2 3 4);         // 熱伝導率 [W/m-K]

また、wallHeatTransfer という境界条件もある。

    wall
    {
        type            wallHeatTransfer;
        Tinf            uniform 500;  // 外部温度 [K]
        alphaWall       uniform 10;   // 熱伝達係数 [m/m2-K]
        value           uniform 300;
    }

これは buoyantBoussinesqSimpleFoam などでは利用できない。

表による境界条件の設定

fixedProfile により表形式で分布を持った境界条件を表現できる。

    inlet
    {
        type            fixedProfile;
        profile         csvFile;
        profileCoeffs
        {
            file            "k.csv";
            nHeaderLine     0;
            refColumn       0;
            componentColumns (1);
            separator       ",";
            mergeSeparators no;
        }
        direction       (0 0 1);
        origin          0;
    }

ここでは profile に "csvFile" を指定することで、CSV ファイルを指定するタイプを選択している。profileCoeffs の file で CSV ファイルを指定している。CSV ファイルの中身は次のようになっている。

k.csv

0.01,0.327
0.013,0.367
0.025,0.405
...

式による境界条件の設定

OpenCFD 版では、exprFixedValue により式による境界条件を表現できる。

    inlet
    {
        type            exprFixedValue;
        valueExpr       "-10*face()/area()";
        value           uniform (0 0 0);
    }

valueExpr で式を記述する。これは #eval と同じ形式を受け付けるが、ここでは境界特有の関数を用いることができる。上の例では、面の法線方向に大きさ 10 の速度を設定している。face() は境界面の面積ベクトル、area() はその面積である。マイナスが付いているのは、境界面の面積ベクトルは領域の外側を向いているためである。

    inlet
    {
        type            exprFixedValue;
        h               0.0254;                                                  
        valueExpr       "-10*sin(pi()*pos().y()/$h)*face()/area()"; 
        value           uniform (0 0 0);
    }

流速分布を与えることもできる。pos().y() は面の y 座標である。

非定常計算で速度の時間変動を与えることもできる。

    inlet
    {
        type            exprFixedValue;
        T               100;
        valueExpr       "-10*sin(2*pi()*time()/$T)*face()/area()"; 
        value           uniform (0 0 0);
    }

time() は時刻である。