【OpenFOAMv2012】kEqnの解説

OpenFOAMで実装されている1方程式LESモデルであるkEqnについて説明する

keywords: OpenFOAM，乱流モデル，kEqn

k Equationモデル
OpenFOAMにおける実装
まとめ

k Equationモデル

基礎方程式

LESの基礎方程式はフィルターを施した連続の式とNavier-Stokes方程式である

フィルターをかけた基礎方程式は次のように表される

（連続の式）

\begin{equation}
\frac{\partial \overline{u}_{i}}{\partial x_{i}}=0
\end{equation}

（Navier-Stokes方程式）

\begin{equation}
\frac{\partial \overline{u}_{i}}{\partial t}
+\frac{\partial \left(\overline{u}_{i}\overline{u}_{j}\right)}{\partial x_{j}}
=
-\frac{1}{\rho}\frac{\partial \overline{P}}{\partial x_{i}}
+\frac{\partial}{\partial x_{j}}\left\{\nu_{e}\left(\frac{\partial \overline{u}_{i}}{\partial x_{j}}+\frac{\partial \overline{u}_{j}}{\partial x_{i}}\right)\right\}
\end{equation}

\begin{equation}
\nu_{e}=\nu+\nu_{_{SGS}}, \qquad \overline{P}=\overline{p}+\frac{2}{3}\rho k_{_{SGS}}
\end{equation}

\(\nu_{_{SGS}}\)はSGS渦粘性係数であり，\(\overline{S}_{ij}\)はひずみ速度テンソル \(\overline{S}_{ij}=\frac{1}{2}\left(\frac{\partial \overline{u}_{i}}{\partial x_{j}}+\frac{\partial \overline{u}_{j}}{\partial x_{i}}\right)\) である

\(\overline{u}_{i}\)，\(\overline{p}\)以外の変数はSGS渦粘性係数\(\nu_{_{SGS}}\)ただ1つであり，この\(\nu_{_{SGS}}\)をどのように求めるかがSGS渦粘性モデルの唯一にして最も重要な点である

上式において，SGSの乱れの効果を表すSGS応力\(\tau_{ij}\)は渦粘性近似によって次のようにモデル化されている

\begin{equation}
\tau_{ij}=\overline{u_{i}u_{j}}-\overline{u}_{i}\overline{u}_{j}
=\frac{2}{3}k_{_{SGS}}\delta_{ij}-2\nu_{_{SGS}}\overline{S}_{ij}
\end{equation}

k Equationモデル

k Equationモデルでは，SGS渦粘性係数\(\nu_{_{SGS}}\)は次の式で定義される

\begin{equation}
\nu_{_{SGS}}=C_{k}\Delta\sqrt{k_{_{SGS}}}
\end{equation}

\(k_{_{SGS}}\)は以下の輸送方程式を解くことによって求められる

\begin{equation}
\frac{\partial k_{_{SGS}}}{\partial t}
+\overline{u}_{j}\frac{\partial k_{_{SGS}}}{\partial x_{j}}
=
-\tau_{ij}\frac{\partial \overline{u}_{i}}{\partial x_{j}}
-C_{\varepsilon}\frac{k_{_{SGS}}^{\frac{3}{2}}}{\Delta}
+\frac{\partial}{\partial x_{j}}\left\{\left(\nu+\nu_{_{SGS}}\right)\frac{\partial k_{_{SGS}}}{\partial x_{j}}\right\}
\end{equation}

モデル定数は以下の通り

\begin{equation}
C_{k}=0.094, \qquad C_{\varepsilon}=1.048
\end{equation}

参考
≫Yoshizawa, A., Horiuti, K. A Statistically-Derived Subgrid-Scale Kinetic Energy Model for the Large-Eddy Simulation of Turbulent Flows (1985) Journal of the Physical Society of Japan, 54 (8), pp. 2834-2839.

OpenFOAMにおける実装

kEqnの実装を確認する

参考
≫OpenFOAM User Guide k Equation
≫kEqn.C
≫kEqn.H

実際に計算が行われるときは，simpleFOAMやpisoFOAMなどのソルバーで turbulence->correct() が呼び出されているので，乱流モデルのプログラムは correct() → correctNut() → correctNonlinearStress() の順番で実行されている

simpleFoam.C ll.97-113

while (simple.loop())
{
    Info<< "Time = " << runTime.timeName() << nl << endl;

    // --- Pressure-velocity SIMPLE corrector
    {
        #include "UEqn.H"
        #include "pEqn.H"
    }

    laminarTransport.correct();
    turbulence->correct();

    runTime.write();

    runTime.printExecutionTime(Info);
    }

≫applications/solvers/incompressible/simpleFoam/simpleFoam.C

渦粘性係数

\begin{equation}
\nu_{_{SGS}}=C_{k}\Delta\sqrt{k_{_{SGS}}}
\end{equation}

kEqn.C ll.41-49

template<class BasicTurbulenceModel>
void kEqn<BasicTurbulenceModel>::correctNut()
{
    this->nut_ = Ck_*sqrt(k_)*this->delta();
    this->nut_.correctBoundaryConditions();
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}

\(k_{_{SGS}}\)

\begin{equation}
G=2\nu_{_{SGS}}\overline{S}_{ij}\frac{\partial \overline{u}_{i}}{\partial x_{j}}
\end{equation}

\begin{equation}
\frac{\partial k_{_{SGS}}}{\partial t}
+\overline{u}_{j}\frac{\partial k_{_{SGS}}}{\partial x_{j}}
-\nu_{e}\frac{\partial}{\partial x_{j}}\left\{\frac{\partial k_{_{SGS}}}{\partial x_{j}}\right\}
=
G
-\frac{2}{3}\frac{\partial \overline{u}_{i}}{\partial x_{i}}k
-C_{\varepsilon}\frac{k_{_{SGS}}^{\frac{3}{2}}}{\Delta}
\end{equation}

kEqn.C ll.199-222

    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));

    tmp<volTensorField> tgradU(fvc::grad(U));
    volScalarField G(this->GName(), nut*(tgradU() && dev(twoSymm(tgradU()))));
    tgradU.clear();

    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(alpha, rho, k_)
      + fvm::div(alphaRhoPhi, k_)
      - fvm::laplacian(alpha*rho*DkEff(), k_)
     ==
        alpha*rho*G
      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
      - fvm::Sp(this->Ce_*alpha*rho*sqrt(k_)/this->delta(), k_)
      + kSource()
      + fvOptions(alpha, rho, k_)
    );

    kEqn.ref().relax();
    fvOptions.constrain(kEqn.ref());
    solve(kEqn);
    fvOptions.correct(k_);
    bound(k_, this->kMin_);

kEqn.H ll.157-164

    //- Return the effective diffusivity for k
    tmp<volScalarField> DkEff() const
    {
        return tmp<volScalarField>
        (
            new volScalarField("DkEff", this->nut_ + this->nu())
        );
    }

\(-\frac{2}{3}\frac{\partial \overline{u}_{i}}{\partial x_{i}}k\)は圧縮性の影響（\(\frac{\partial \overline{u}_{i}}{\partial x_{i}} \neq 0\)）で出てきた項である

\begin{equation}
-\tau_{ij}\frac{\partial \overline{u}_{i}}{\partial x_{j}}
=-\left(\frac{2}{3}k\delta_{ij}-2\nu_{e}\overline{S}_{ij}\right)
\frac{\partial \overline{u}_{i}}{\partial x_{j}}
=-\frac{2}{3}\frac{\partial \overline{u}_{k}}{\partial x_{k}}k
+2\nu_{e}\overline{S}_{ij}\frac{\partial \overline{u}_{i}}{\partial x_{j}}
\end{equation}

定数群

\begin{equation}
C_{k}=0.094
\end{equation}

kEqn.C ll.107-115

    Ck_
    (
        dimensioned<scalar>::getOrAddToDict
        (
            "Ck",
            this->coeffDict_,
            0.094
        )
    )

\begin{equation}
C_{\varepsilon}=1.048
\end{equation}

Smagoriskyの親クラスであるLESeddyViscosityで定義されている
≫LESeddyViscosity.C

kEqn.C ll.66-74

    Ce_
    (
        dimensioned<scalar>::getOrAddToDict
        (
            "Ce",
            this->coeffDict_,
            1.048
        )
    )

まとめ

OpenFOAMで実装されている1方程式LESモデルであるkEqnについて説明した

そのほかの乱流モデルについてはこちら