【OpenFOAMv2012】dynamicKEqnの解説

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

keywords: OpenFOAM，乱流モデル，dynamicKEqn

dynamicKEqn
OpenFOAMにおける実装
まとめ

dynamicKEqn

基礎方程式

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}

\(\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{D}_{ij}
\end{equation}

\(\overline{D}_{ij}\)はひずみ速度テンソルである

\begin{equation}
\overline{D}_{ij}=\frac{1}{2}\left(\frac{\partial \overline{u}_{i}}{\partial x_{j}}+\frac{\partial \overline{u}_{j}}{\partial x_{i}}\right)
\end{equation}

dynamicKEqn

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

\begin{equation}
\nu_{_{SGS}}=C_{k}\overline{\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}}}{\overline{\Delta}}
+\frac{\partial}{\partial x_{j}}\left\{\left(\nu+\nu_{_{SGS}}\right)\frac{\partial k_{_{SGS}}}{\partial x_{j}}\right\}
\end{equation}

\(C_{k}\)および\(C_{\varepsilon}\)は次のように与えられる

\begin{align}
&C_{k}=-\frac{1}{2}\frac{L_{ij}^{a}M_{ij}}{M_{kl}M_{kl}} \\ \\

&C_{\varepsilon}=\frac{\left(\nu+\nu_{_{SGS}}\right)\widetilde{\Delta}}{k_{_{test}}^{3/2}}
\left(\widetilde{\frac{\partial \overline{u}_{i}}{\partial x_{j}}\frac{\partial \overline{u}_{i}}{\partial x_{j}}}-\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\frac{\partial \widetilde{\overline{u}}_{i}}{\partial x_{j}}\right)
\end{align}

\begin{align}
&L_{ij}^{a}=\mathrm{dev}(\widetilde{\overline{u}_{i}\overline{u}_{j}}
-\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}) \\ \\
&M_{ij}=-\widetilde{\Delta}\sqrt{k_{_{test}}}\widetilde{\overline{D}}_{ij}, \qquad
k_{_{test}}=\frac{L_{ii}}{2}
\end{align}

ここで，\(\overline{G}\)はグリッドフィルター，\(\tilde{G}\)はテストフィルターである

参考
≫W-W Kim and S. Menon. A new dynamic one-equation subgrid-scale model for large eddy simulations. In 33rd Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 1995.

以下のサイトにあるPDFのpp.184-193が上記の論文
≫A New Approach to Validate Subgrid Models in Complex High Reynolds Number Flows.

OpenFOAMにおける実装

dynamicKEqnの実装を確認する

OpenFOAMのdelta()で参照しているのはグリッドフィルター幅\(\overline{\Delta}\)であり，dynamicKEqnの実装では\(\tilde{\Delta}=2\overline{\Delta}\)としている

参考
≫OpenFOAM User Guide Dynamic k eqn
≫dynamicKEqn.C
≫dynamicKEqn.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}\overline{\Delta}\sqrt{k_{_{SGS}}}
\end{equation}

dynamicKEqn.C ll.111-123

template<class BasicTurbulenceModel>
void dynamicKEqn<BasicTurbulenceModel>::correctNut
(
    const volSymmTensorField& D,
    const volScalarField& KK
)
{
    this->nut_ = Ck(D, KK)*sqrt(k_)*this->delta();
    this->nut_.correctBoundaryConditions();
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}

\(k_{_{SGS}}\)

\begin{equation}
\overline{D}_{ij}=\mathrm{dev}\left\{\frac{1}{2}\left(\frac{\partial \overline{u}_{i}}{\partial x_{j}}+\frac{\partial \overline{u}_{j}}{\partial x_{i}}\right)\right\}
\end{equation}

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

\begin{align}
k_{_{test}}=\frac{1}{2}L_{ii}=\frac{1}{2}\left(\widetilde{\overline{u}_{i}\overline{u}_{i}}
-\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{i}\right)
\end{align}

dynamicKEqn.C ll.279-287

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

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

    volScalarField KK(0.5*(filter_(magSqr(U)) - magSqr(filter_(U))));
    KK.max(dimensionedScalar("small", KK.dimensions(), SMALL));

\begin{equation}
\frac{\partial k_{_{SGS}}}{\partial t}
+\overline{u}_{j}\frac{\partial k_{_{SGS}}}{\partial x_{j}}
-\frac{\partial}{\partial x_{j}}\left\{\left(\nu+\nu_{_{SGS}}\right)\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}}}{\overline{\Delta}}
\end{equation}

dynamicKEqn.C ll.289-306

    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(Ce(D, KK)*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_);

dynamicKEqn.H ll.187-194

        //- 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{align}
&L_{ij}^{a}=\mathrm{dev}(\widetilde{\overline{u}_{i}\overline{u}_{j}}
-\widetilde{\overline{u}}_{i}\widetilde{\overline{u}}_{j}) \\ \\
&M_{ij}=-2\overline{\Delta}\sqrt{k_{_{test}}}\widetilde{\overline{D}}_{ij} \\ \\
&C_{k}=-\frac{1}{2}\frac{L_{ij}^{a}M_{ij}}{M_{kl}M_{kl}} \\ \\
\end{align}

dynamicKEqn.C ll.41-75

template<class BasicTurbulenceModel>
volScalarField dynamicKEqn<BasicTurbulenceModel>::Ck
(
    const volSymmTensorField& D,
    const volScalarField& KK
) const
{
    const volSymmTensorField LL
    (
        simpleFilter_(dev(filter_(sqr(this->U_)) - (sqr(filter_(this->U_)))))
    );

    const volSymmTensorField MM
    (
        simpleFilter_
        (
            -2.0*this->delta()*sqrt
            (
                max(KK, dimensionedScalar(KK.dimensions(), Zero))
            )*filter_(D)
        )
    );

    const volScalarField Ck
    (
        simpleFilter_(0.5*(LL && MM))
       /(
            simpleFilter_(magSqr(MM))
          + dimensionedScalar("small", sqr(MM.dimensions()), VSMALL)
        )
    );

    tmp<volScalarField> tfld = 0.5*(mag(Ck) + Ck);
    return tfld();
}

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

dynamicKEqn.C ll.78-93

template<class BasicTurbulenceModel>
volScalarField dynamicKEqn<BasicTurbulenceModel>::Ce
(
    const volSymmTensorField& D,
    const volScalarField& KK
) const
{
    const volScalarField Ce
    (
        simpleFilter_(this->nuEff()*(filter_(magSqr(D)) - magSqr(filter_(D))))
       /simpleFilter_(pow(KK, 1.5)/(2.0*this->delta()))
    );

    tmp<volScalarField> tfld = 0.5*(mag(Ce) + Ce);
    return tfld();
}

\(\frac{\partial \overline{u}_{i}}{\partial x_{j}}\frac{\partial \overline{u}_{i}}{\partial x_{j}}=D_{ij}D_{ij}\)としている原因は不明（等しくはないはず）

まとめ

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

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