【OpenFOAMv2012】kEpsilonの解説

OpenFOAMで実装されている2方程式RANSモデルであるkEpsilonついて説明する

keywords: OpenFOAM,乱流モデル,kEpsilon

スポンサーリンク

標準k-εモデル

RANSにおける渦粘性モデルの基礎方程式は次の式で表される

\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_{t}, \qquad \overline{P}=\overline{p}+\frac{2}{3}\rho k
\end{equation}

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

k-εモデルの中でも最も基本となる「標準k-εモデル」とよばれるモデルは,乱流エネルギー\(k\)と乱流散逸率\(\varepsilon\)を用いて次のように渦粘性係数を計算している

\begin{equation}
\nu_{t}=C_{\mu}\frac{k^{2}}{\varepsilon}
\end{equation}

乱流エネルギー\(k\)と乱流散逸率\(\varepsilon\)はそれぞれの輸送方程式を解くことで求められる

\begin{equation}
\frac{\partial k}{\partial t}+\overline{u}_{j}\frac{\partial k}{\partial x_{j}}
=P_{k}-\varepsilon
+\frac{\partial}{\partial x_{j}}\left\{\left(\frac{\nu_{t}}{\sigma_{k}}+\nu\right)\frac{\partial k}{\partial x_{j}}\right\}
\end{equation}

\begin{equation}
\frac{\partial \varepsilon}{\partial t}+\overline{u}_{j}\frac{\partial \varepsilon}{\partial x_{j}}
=\left(C_{\varepsilon_{1}}P_{k}-C_{\varepsilon_{2}}\varepsilon\right)\frac{\varepsilon}{k}+\frac{\partial}{\partial x_{j}}\left\{\left(\frac{\nu_{t}}{\sigma_{\varepsilon}}+\nu\right)\frac{\partial \varepsilon}{\partial x_{j}}\right\}
\end{equation}

\(P_{k}\)は乱流エネルギー\(k\)の生成項で次の式で表される(\(S_{ij}=\frac{1}{2}\left(\frac{\partial \overline{u}_{i}}{\partial x_{j}}+\frac{\partial \overline{u}_{j}}{\partial x_{i}}\right)\))

\begin{equation}
P_{k}=2\nu_{t}S_{ij}\frac{\partial \overline{u}_{i}}{\partial x_{j}}
\end{equation}

ここで,モデル定数は次の通り

\begin{equation}
C_{\mu}=0.09, \quad \sigma_{k}=1.0, \quad \sigma_{\varepsilon}=1.3, \quad
C_{\varepsilon_{1}}=1.44, \quad C_{\varepsilon_{2}}=1.92
\end{equation}

参考
Launder, B.E., Spalding, D.B. The numerical computation of turbulent flows (1974) Computer Methods in Applied Mechanics and Engineering, 3 (2), pp. 269-289.

OpenFOAMにおける実装

kEpsilonの実装を確認する

参考
OpenFOAM User Guide k-epsilon
kEpsilon< BasicTurbulenceModel > Class Template Reference
kEpsilon.C
kEpsilon.H

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

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_{t}=C_{\mu}f_{\mu}\frac{k^{2}}{\varepsilon}
\end{equation}

template<class BasicTurbulenceModel>
void kEpsilon<BasicTurbulenceModel>::correctNut()
{
    this->nut_ = Cmu_*sqr(k_)/epsilon_;
    this->nut_.correctBoundaryConditions();
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}

輸送方程式

\begin{equation}
P_{k}=2\nu_{t}S_{ij}\frac{\partial \overline{u}_{i}}{\partial x_{j}}
\end{equation}

    tmp<volTensorField> tgradU = fvc::grad(U);
    const volScalarField::Internal GbyNu
    (
        this->type() + ":GbyNu",
        tgradU().v() && dev(twoSymm(tgradU().v()))
    );
    const volScalarField::Internal G(this->GName(), nut()*GbyNu);
    tgradU.clear();

\begin{align}
\frac{\partial \tilde{\varepsilon}}{\partial t}+\overline{u}_{j}\frac{\partial \tilde{\varepsilon}}{\partial x_{j}}
&=\left(C_{\varepsilon_{1}}P_{k}-C_{\varepsilon_{2}}f_{2}\tilde{\varepsilon}\right)\frac{\tilde{\varepsilon}}{k}+E+\frac{\partial}{\partial x_{j}}\left\{\left(\frac{\nu_{t}}{\sigma_{\varepsilon}}+\nu\right)\frac{\partial \tilde{\varepsilon}}{\partial x_{j}}\right\} \\
&\\
\frac{\partial k}{\partial t}+\overline{u}_{j}\frac{\partial k}{\partial x_{j}}
&=P_{k}-\left(\tilde{\varepsilon}+D\right)
+\frac{\partial}{\partial x_{j}}\left\{\left(\frac{\nu_{t}}{\sigma_{k}}+\nu\right)\frac{\partial k}{\partial x_{j}}\right\}
\end{align}

    // Update epsilon and G at the wall
    epsilon_.boundaryFieldRef().updateCoeffs();

    // Dissipation equation
    tmp<fvScalarMatrix> epsEqn
    (
        fvm::ddt(alpha, rho, epsilon_)
      + fvm::div(alphaRhoPhi, epsilon_)
      - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
     ==
        C1_*alpha()*rho()*GbyNu*Cmu_*k_()
      - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_)
      - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
      + epsilonSource()
      + fvOptions(alpha, rho, epsilon_)
    );

    epsEqn.ref().relax();
    fvOptions.constrain(epsEqn.ref());
    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
    solve(epsEqn);
    fvOptions.correct(epsilon_);
    bound(epsilon_, this->epsilonMin_);

    // Turbulent kinetic energy equation
    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(alpha()*rho()*epsilon_()/k_(), k_)
      + kSource()
      + fvOptions(alpha, rho, k_)
    );

    kEqn.ref().relax();
    fvOptions.constrain(kEqn.ref());
    solve(kEqn);
    fvOptions.correct(k_);
    bound(k_, this->kMin_);
        //- Return the effective diffusivity for k
        tmp<volScalarField> DkEff() const
        {
            return tmp<volScalarField>
            (
                new volScalarField
                (
                    "DkEff",
                    (this->nut_/sigmak_ + this->nu())
                )
            );
        }

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

定数群

\begin{equation}
C_{\mu}=0.09, \quad
C_{\varepsilon_{1}}=1.44, \quad
C_{\varepsilon_{2}}=1.92, \quad
\sigma_{k}=1.0, \quad
\sigma_{\varepsilon}=1.3
\end{equation}

    Cmu_
    (
        dimensioned<scalar>::getOrAddToDict
        (
            "Cmu",
            this->coeffDict_,
            0.09
        )
    ),
    C1_
    (
        dimensioned<scalar>::getOrAddToDict
        (
            "C1",
            this->coeffDict_,
            1.44
        )
    ),
    C2_
    (
        dimensioned<scalar>::getOrAddToDict
        (
            "C2",
            this->coeffDict_,
            1.92
        )
    ),
    C3_
    (
        dimensioned<scalar>::getOrAddToDict
        (
            "C3",
            this->coeffDict_,
            0
        )
    ),
    sigmak_
    (
        dimensioned<scalar>::getOrAddToDict
        (
            "sigmak",
            this->coeffDict_,
            1.0
        )
    ),
    sigmaEps_
    (
        dimensioned<scalar>::getOrAddToDict
        (
            "sigmaEps",
            this->coeffDict_,
            1.3
        )
    ),

まとめ

OpenFOAMで実装されている2方程式RANSモデルであるkEpsilonついて説明した

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

コメント

タイトルとURLをコピーしました