物性値の設定

2020年8月13日

はじめに

物性値の設定について。

使用バージョン

OpenFOAM v2006、v8

非圧縮性流体ソルバーの物性値の設定

simpleFoam などの非圧縮性流体ソルバーの場合、物性値は constant/transportProperties で設定する。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}

transportModel  Newtonian;

nu              nu [ 0 2 -1 0 0 0 0 ] 1e-05;

transportModel は物性の種類を設定する。以下のようなものが指定できる。

  • Newtonian
  • BirdCarreau
  • CrossPowerLaw
  • powerLaw
  • HerschelBulkley
  • Casson
  • strainRateFunction

ニュートン流体 (Newtonian) と非ニュートン流体 (Newtonian 以外) を選ぶことができる。それぞれ設定方法が異なる。

Newtonian

Newtonian の場合、nu で動粘性係数 (粘性係数/密度) を指定する。

nu              nu [ 0 2 -1 0 0 0 0 ] 1e-05; // [m2/s]

単純に次のように書いてもよい (古いバージョンだと不可)。

nu              [ 0 2 -1 0 0 0 0 ] 1e-05;

あるいは、さらに単純に

nu              1e-05;

BirdCarreau

Bird-Carreau モデルでは、動粘性係数は次式で表される。

$$\nu = \nu_\infty + (\nu_0 - \nu_\infty)\left\{ 1 + (k \dot{\gamma})^a \right\}^{(n - 1)/a}$$

\(a\) はデフォルトでは 2 である。設定は以下のように行う。

BirdCarreauCoeffs
{
    nu0         [ 0 2 -1 0 0 0 0 ] 0.1;
    nuInf       [ 0 2 -1 0 0 0 0 ] 1e-3;
    k           [ 0 0  1 0 0 0 0 ] 1;
    n           [ 0 0  0 0 0 0 0 ] 0.5;
}

CrossPowerLaw

Cross べき乗則モデルでは、動粘性係数は次式で表される。

$$\nu = \nu_\infty + \frac{\nu_0 - \nu_\infty}{1 + (m \dot{\gamma})^n}$$

設定は以下のように行う。

CrossPowerLawCoeffs
{
    nu0         [ 0 2 -1 0 0 0 0 ] 0.1;
    nuInf       [ 0 2 -1 0 0 0 0 ] 1e-3;
    m           [ 0 0  1 0 0 0 0 ] 1;
    n           [ 0 0  0 0 0 0 0 ] 0.5;
}

powerLaw

べき乗則モデルでは、動粘性係数は次式で表される。

$$\nu = k \dot{\gamma}^{n - 1} \hspace{8pt} \nu_\mathrm{min} \ge \nu \ge \nu_\mathrm{max}$$

設定は以下のように行う。

powerLawCoeffs
{
    k               [ 0 2 -1 0 0 0 0 ] 0.1;
    n               [ 0 0  0 0 0 0 0 ] 1;
    nuMin           [ 0 2 -1 0 0 0 0 ] 0.1;
    nuMax           [ 0 2 -1 0 0 0 0 ] 1;
}

HerschelBulkley

Herschel-Bulkley モデルはビンガム塑性とべき乗則が組み合わされたもので、動粘性係数は次式で表される。

$$\nu = \min(\nu_0, \tau_0/\dot{\gamma} + k \dot{\gamma}^{n - 1})$$

設定は以下のように行う。

HerschelBulkleyCoeffs
{
    k               [ 0 2 -1 0 0 0 0 ] 1e-3;
    n               [ 0 0  0 0 0 0 0 ] 1;
    tau0            [ 0 2 -2 0 0 0 0 ] 1;
    nu0             [ 0 2 -1 0 0 0 0 ] 0.1;
}

Casson

Casson モデルは血液などを表現するためのモデルで、動粘性係数は次式で表される。

$$\nu = (\sqrt{\tau_0/\dot{\gamma}} + \sqrt{m})^2 \hspace{8pt} \nu_\mathrm{min} \ge \nu \ge \nu_\mathrm{max}$$

設定は以下のように行う (血液のパラメータ)。

CassonCoeffs
{
    // blood
    m           [ 0 2 -1 0 0 0 0 ] 3.934986e-6;
    tau0        [ 0 2 -2 0 0 0 0 ] 2.9032e-6;
    nuMin       [ 0 2 -1 0 0 0 0 ] 3.9047e-6;
    nuMax       [ 0 2 -1 0 0 0 0 ] 13.333e-6;
}

strainRateFunction

strainRateFunction では一般のひずみ速度関数を設定できる。

strainRateFunctionCoeffs
{
    function polynomial ((0 0.1) (1 1.3));
}

Arrhenius 型

(v1712 以上) 温度依存を考慮したアレニウス型のモデルが用意されている。

  • ArrheniusNewtonian
  • ArrheniusBirdCarreau
  • ArrheniusCrossPowerLaw
  • ArrheniusHerschelBulkley
  • ArrheniusCasson

それぞれの粘性係数モデルに以下の係数が掛けられる。

$$\nu_a = \exp(-\alpha(T - T_\alpha))$$

パラメータとして alpha, Talpha が要求される。

非圧縮性流体ソルバー用の設定なので、基本的に温度計算はされないが、これは function object の energyTransport を用いた温度計算を併用することが前提されている。

Boussinesq 近似ソルバー

buoyantBoussinesqSimpleFoam のような浮力に Boussinesq 近似を用いたソルバーは、非圧縮性流体でありながら浮力を考慮できる。これらのソルバーには、浮力と熱伝導に関する設定が必要である。

// Thermal expansion coefficient
beta            beta [0 0 0 -1 0 0 0] 3e-03; // 体積膨張係数 [1/K]

// Reference temperature
TRef            TRef [0 0 0 1 0 0 0] 300; // 参照温度 [K]

// Laminar Prandtl number
Pr              Pr [0 0 0 0 0 0 0] 0.9; // 層流プラントル数

// Turbulent Prandtl number
Prt             Prt [0 0 0 0 0 0 0] 0.7; // 乱流プラントル数

熱拡散率 \(\alpha\) が動粘性係数 \(\nu\) とプラントル数 \(Pr\) (粘性係数×比熱/熱伝導率) から次式で計算される

$$\alpha = \frac{\nu}{Pr}$$

圧縮性流体ソルバーの物性値の設定

buoyantSimpleFoam などの圧縮性流体ソルバーの場合は、物性値を constant/thermophysicalProperties で設定する。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      thermophysicalProperties;
}

thermoType
{
    type                heRhoThermo;
    mixture             pureMixture;
    specie              specie;
    equationOfState     perfectGas;
    transport           const;
    thermo              hConst;
    energy              sensibleEnthalpy;
}

mixture
{
    specie
    {
        molWeight       28.9;
    }
    thermodynamics
    {
        Cp              1000;
        Hf              0;
    }
    transport
    {
        mu              1.8e-05;
        Pr              0.7;
    }
}

thermoType は物性の種類を設定する。いくつかの項目があるが、それぞれで設定するのではなく、可能な項目の組み合わせがソルバーごとに決まっている。また、設定できたとしても、ソルバーがそれに関する計算に対応しているとは限らない (あくまで物性値を扱う枠組みにすぎないため)。

type

熱物理モデル (thermophysical model) の設定。heRhoThermo や hePsiThermo などを指定。Rho と Psi の違いは、密度を直接扱うか、圧縮率 (compressibility) から計算するかの違いである。Psi は理想気体を想定しており、Rho はより一般的である。

mixture

化学種の扱い方の設定。pureMixture, multiComponentMixture, reactingMixture などを指定。それぞれ単一化学種、複数化学種、化学種の反応を扱うものだが、実際に計算でそれらが扱えるかどうかはソルバーによる。

pureMixture

pureMixture の場合、上記のように mixture で物性値を指定する。物性値をどのように指定するかは、thermoType の設定による。

multiComponentMixture

multiComponentMixture の場合、species で化学種を指定、その化学種ごとに物性値を指定する。

species
(
    air
    H2O
);

air
{
    specie
    {
        molWeight       28.9;
    }
    thermodynamics
    {
        Cp              1000;
        Hf              0;
    }
    transport
    {
        mu              1.8e-05;
        Pr              0.7;
    }
}

H2O
{
    ...
}

化学種ごとにフィールドファイルが必要である。

reactingMixture

reactingMixture の場合、つぎのように指定する。

inertSpecie N2;

chemistryReader foamChemistryReader;

foamChemistryFile "$FOAM_CASE/constant/reactions";

foamChemistryThermoFile "$FOAM_CASE/constant/thermo.compressibleGas";

inertSpecie で不活性の化学種 (基本となる化学種) を指定する。chemistryReader は設定方法の設定で、foamChemistryReader と chemkinReader が指定できる。foamChemistryReader の場合は foamChemistryFile、foamChemistryThermoFile でつぎのような内容のファイルを指定する。

reactions

species
(
    O2
    H2O
    CH4
    CO2
    N2
);

reactions
{
    methaneReaction
    {
        type         irreversibleArrheniusReaction;
        reaction     "CH4 + 2O2 = CO2 + 2H2O";
        A            5.2e16;
        beta         0;
        Ta           14906;
    }
}

thermo.compressibleGas

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      thermo.compressibleGas;
}

O2
{
    specie
    {
        molWeight       31.9988;
    }
    thermodynamics
    {
        Tlow            200;
        Thigh           5000;
        Tcommon         1000;
        highCpCoeffs    ( 3.69758 0.00061352 -1.25884e-07 1.77528e-11 -1.13644e-15 -1233.93 3.18917 );
        lowCpCoeffs     ( 3.21294 0.00112749 -5.75615e-07 1.31388e-09 -8.76855e-13 -1005.25 6.03474 );
    }
    transport
    {
        As              1.67212e-06;
        Ts              170.672;
    }
}

H2O
{
    ...
}

...

multiComponentMixture との違いは、反応を設定するための reactions があるところである。

CHEMKIN のファイルを OpenFOAM の形式に変換する chemkinToFoam というユーティリティがある。

specie

つねに "specie" を指定する。入力としてはモル質量 (物質量) を指定する。

    specie
    {
        molWeight       28.9; // モル質量
    }

OpenFOAM 4.x までは、モル数の指定が必要である。

    specie
    {
        nMoles          1;    // モル数
        molWeight       28.9; // モル質量
    }

equationOfState

状態方程式の設定を行う。以下のようなものを指定できる。

rhoConst
密度を指定する。
    equationOfState
    {
        rho             1000; // 密度 [kg/m3]
    }
perfectGas
密度を完全気体 (理想気体) の状態方程式から計算する。 密度は次式で計算される。 $$\rho = \frac{p}{RT}$$ ここで R は分子量で割られた気体定数である。
incompressiblePerfectGas
密度を完全気体 (理想気体) の状態方程式から計算するが、圧力には参照圧力を用いる。
    equationOfState
    {
        pRef             101325; // 参照圧力 [Pa]
    }
密度は次式で計算される。 $$\rho = \frac{p_\mathrm{ref}}{RT}$$
Boussinesq
(OpenFOAM 4.x 以上) 密度の計算に Boussinesq 近似を用いる。
    equationOfState
    {
        rho0 1.17; // 参照密度 [kg/m3]
        T0   300; // 参照温度 [K]
        beta 3e-03; // 体積膨張係数 [1/K]
    }
密度は次式で計算される。 $$\rho = \rho_0 \{1 - \beta (T - T_0)\}$$
PengRobinsonGas
(OpenFOAM 2.3.x 以上) Peng-Robinson の状態方程式を用いる。
    equationOfState
    {
        Tc 126.2;   // 臨界温度 [K]
        Vc 0.089;   // 臨界体積 [m3/kmol] (= 分子量 [kg/kmol] / 臨界密度 [kg/m3])
        Pc 3.39e6;  // 臨界圧力 [Pa]
        omega 0.04; // 偏心因子 (acentric factor) [-]
    }
perfectFluid
密度を参照密度と理想気体の状態方程式から求める。

    equationOfState
    {
        R               3000; // モル質量当たりの気体定数に相当する定数 [J/kg-K]
        rho0            1025; // 参照密度 [kg/m3]
    }
密度は次式で計算される。 $$\rho = \rho_0 + \frac{p}{RT}$$
icoPolynomial
密度を多項式で表す。
    equationOfState
    {
        rhoCoeffs<8>    ( 4.0097 -0.016954 3.3057e-05 -3.0042e-08 1.0286e-11 0 0 0 );
    }
係数はそれぞれ低次の項から指定する。

transport

輸送特性を設定する。つぎのようなものが指定できる。

const
粘性係数とプラントル数 (粘性係数×比熱/熱伝導率) を設定する。熱伝導率はプラントル数から計算される。
    transport
    {
        mu              1.8e-05; // 粘性係数 [Pa-s]
        Pr              0.7;     // プラントル数
    }
polynomial
粘性係数と熱伝導率を多項式で設定する。
    transport
    {
          muCoeffs<8>     ( 1.5061e-06 6.16e-08 -1.819e-11 0 0 0 0 0 );
          kappaCoeffs<8>  ( 0.0025219 8.506e-05 -1.312e-08 0 0 0 0 0 );
    }
係数はそれぞれ低次の項から指定する。
sutherland
Sutherland の式による設定。
    transport
    {
        As              1.67212e-06;
        Ts              170.672;
    }
Sutherland の式により、粘性係数は次式で計算される。 $$\mu = \frac{As \sqrt{T}}{1 + T_s/T}$$ 熱伝導率は次の modified Eucken correlation で計算される。 $$k = \mu C_v (1.32 + 1.77 R/C_v)$$ ここで Cv は単位質量あたりの定積比熱 [J/kg-K] である。

thermo

熱特性を設定する。つぎのようなものを指定できる。

hConst
比熱と標準生成エンタルピーを設定する。
    thermodynamics
    {
        Cp              1000; // 比熱 [J/kg-K]
        Hf              0;    // 標準生成エンタルピー [J/kg]
    }
hPolynomial
比熱を温度の多項式で設定する。
    thermodynamics
    {
        Hf              0; // 標準生成エンタルピー [J/kg]
        Sf              0; // 標準エントロピー [J/kg-K]
        CpCoeffs<8>     ( 948.76 0.39171 -0.00095999 1.393e-06 -6.2029e-10 0 0 0 );
    }
係数 CpCoeffs は低次の項から指定する。
janaf
JANAF 表形式で設定する。
    thermodynamics
    {
        Tlow            200;
        Thigh           5000;
        Tcommon         1000;
        highCpCoeffs    ( 3.69758 0.00061352 -1.25884e-07 1.77528e-11 -1.13644e-15 -1233.93 3.18917 );
        lowCpCoeffs     ( 3.21294 0.00112749 -5.75615e-07 1.31388e-09 -8.76855e-13 -1005.25 6.03474 );
    }
比熱を温度の多項式で表す。温度の範囲と多項式の係数を指定する。比熱は次式で計算される。 $$C_p = R (((a_4 T + a_2) T + a_1) T + a_0)$$

energy

エネルギーの種類を設定する。以下のようなものを指定する。

  • sensibleEnthalpy
  • sensibleInternalEnergy
  • absoluteEnthalpy
  • absoluteInternalEnergy

圧力の時間微分項の考慮

物性ではないが、sensibleEnthalpy の場合、エネルギー方程式で用いられる圧力の時間微分項の有無を指定できる (実際にはソルバー側の対応が必要)。

dpdt no;