PR

【OpenFOAMv2012】Smagorinskyの解説

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

keywords: OpenFOAM,乱流モデル,Smagorinsky

スポンサーリンク

Smagorinskyモデル

基礎方程式

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}

一般的なSmagorinkyモデル

参考書などに載っているSmagorinskyモデルでは,SGS渦粘性係数は次の式で定義される

\begin{equation}
\nu_{_{SGS}}=\left(C_{s}\Delta\right)^{2}\sqrt{2\overline{S}_{ij}\overline{S}_{ij}}
\end{equation}

\(\Delta\)はフィルター幅であり,\(C_{s}\)はこのモデルで与える必要のある唯一のモデル定数でSmagorinsky定数とよばれる

Smagorinsky定数の理論値は \(C_{s}=0.173\) であり,一様等方乱流では実験結果とよく一致することが確認されている

↓参考

OpenFOAMで実装されているSmagorinskyモデル

OpenFOAMで実装されているSmagorinskyモデルでは,SGS応力およびSGS渦粘性係数は次の式で定義される

\begin{align}
&\tau_{ij}=\frac{2}{3}k_{_{SGS}}\delta_{ij}-2\nu_{_{SGS}}\mathrm{dev}(\overline{S}_{ij}), \qquad
\mathrm{dev}(\overline{S}_{ij})=\overline{S}_{ij}-\frac{1}{3}S_{kk}\delta_{ij} \\ \\
&\nu_{_{SGS}}=C_{k}\Delta\sqrt{k_{_{SGS}}}
\end{align}

\(k_{_{SGS}}\)はSGSの運動エネルギー \(k_{_{SGS}}=\frac{1}{2}\left(\overline{u_{k}u_{k}}-\overline{u}_{k}\overline{u}_{k}\right)=\frac{1}{2}\tau_{kk}\) であり,以下の式で求められる

\begin{equation}
k_{_{SGS}}=\left(\frac{-b+\sqrt{b^{2}+4ac}}{2a}\right)^{2}
\end{equation}

それぞれの係数は次のように与えられる

\begin{align}
&a=\frac{C_{\varepsilon}}{\Delta}, \quad
b=\frac{2}{3}\mathrm{tr}(\overline{S}_{ij}), \quad
c=2C_{k}\Delta\left\{\mathrm{dev}(\overline{S}_{ij})\overline{S}_{ij}\right\} \\ \\
&C_{k}=0.094, \quad
C_{\varepsilon}=1.048
\end{align}

一見するとSmagorinskyモデルとは全く関係ない k Equationモデルのように見えるが,非圧縮性を仮定すると一般的なSmagorinkyモデルと同じものになる

\(k_{_{SGS}}\)の式

流れ場にフィルターをかけると,流れ全体の運動エネルギー \(\overline{k}=\frac{\overline{u_{i}u_{i}}}{2}\)(乱流エネルギーではない)はGS成分 \(k_{_{GS}}=\frac{\overline{u}_{i}\overline{u}_{i}}{2}\)とSGS成分 \(k_{_{SGS}}=\frac{\overline{u_{i}u_{i}}-\overline{u}_{i}\overline{u}_{i}}{2}\) に分けられる

\(k_{_{SGS}}\)の輸送方程式は次のように表される

\begin{align}
&\frac{\partial k_{_{SGS}}}{\partial t}
+\overline{u}_{j}\frac{\partial k_{_{SGS}}}{\partial x_{j}}
=
-\tau_{ij}\overline{S}_{ij}-\varepsilon_{_{SGS}} \\
&\qquad +\frac{\partial}{\partial x_{j}}
\left\{\overline{u}_{i}\tau_{ij}-\frac{\overline{u_{i}u_{i}u_{j}}+\overline{u}_{j}\overline{u_{i}u_{i}}}{2}-\frac{\overline{pu_{j}}-\overline{p}\overline{u}_{j}}{\rho}+\nu\frac{\partial k_{_{SGS}}}{\partial x_{j}}\right\}
\end{align}

ここで,右辺第1項 \(-\tau_{ij}S_{ij}\) はSGS運動エネルギーの生成項,右辺第2項\(\varepsilon_{_{SGS}}\)はSGS運動エネルギーの散逸項である

SGSにおいて運動エネルギーの生成と散逸がつりあっているという局所平衡(GSからSGSに入ってきたエネルギーがすべてSGSにおいて熱として散逸していること)を仮定すると,SGSにおける運動エネルギーの散逸率\(\varepsilon_{_{SGS}}\)は次のようになる

\begin{equation}
\varepsilon_{_{SGS}}=-\tau_{ij} \overline{S}_{ij}
\tag{i}
\end{equation}

一方で,k Equationモデルでは\(\varepsilon\)は次元解析によって次のように計算される

\begin{equation}
\varepsilon_{_{SGS}}=C_{\varepsilon}\frac{k_{_{SGS}}^{\frac{3}{2}}}{\Delta}
\tag{ii}
\end{equation}

(i)と(ii)より

\begin{equation}
\tau_{ij} \overline{S}_{ij}+C_{\varepsilon}\frac{k_{_{SGS}}^{\frac{3}{2}}}{\Delta}=0
\end{equation}

SGS応力 \(\tau_{ij}=\frac{2}{3}k_{_{SGS}}\delta_{ij}-2\nu_{_{SGS}}\mathrm{dev}(\overline{S}_{ij})\) およびSGS渦粘性係数 \(\nu_{_{SGS}}=C_{k}\Delta\sqrt{k_{_{SGS}}}\) を代入して

\begin{align}
&\left\{\frac{2}{3}k_{_{SGS}}\delta_{ij}-2\nu_{_{SGS}}\mathrm{dev}(\overline{S}_{ij})\right\}\overline{S}_{ij}+C_{\varepsilon}\frac{k_{_{SGS}}^{\frac{3}{2}}}{\Delta}=0 \\
\rightarrow \quad
&\left\{\frac{2}{3}k_{_{SGS}}\delta_{ij}-2C_{k}\Delta\sqrt{k_{_{SGS}}}\mathrm{dev}(\overline{S}_{ij})\right\}\overline{S}_{ij}+C_{\varepsilon}\frac{k_{_{SGS}}^{\frac{3}{2}}}{\Delta}=0 \\
\rightarrow \quad
&\frac{2}{3}k_{_{SGS}}\overline{S}_{ij}\delta_{ij}-2C_{k}\Delta\sqrt{k_{_{SGS}}}\mathrm{dev}(\overline{S}_{ij})\overline{S}_{ij}+C_{\varepsilon}\frac{k_{_{SGS}}^{\frac{3}{2}}}{\Delta}=0 \\
\rightarrow \quad
&\sqrt{k_{_{SGS}}}\left\{\frac{2}{3}\sqrt{k_{_{SGS}}}\overline{S}_{kk}-2C_{k}\Delta\mathrm{dev}(\overline{S}_{ij})\overline{S}_{ij}+\frac{C_{\varepsilon}}{\Delta}k_{_{SGS}}\right\}=0 \\
\rightarrow \quad
&\frac{C_{\varepsilon}}{\Delta}k_{_{SGS}}+\frac{2}{3}\sqrt{k_{_{SGS}}}\overline{S}_{kk}-2C_{k}\Delta\mathrm{dev}(\overline{S}_{ij})\overline{S}_{ij}=0 \\
\rightarrow \quad
&\frac{C_{\varepsilon}}{\Delta}\left(\sqrt{k_{_{SGS}}}\right)^{2}+\frac{2}{3}\overline{S}_{kk}\sqrt{k_{_{SGS}}}-2C_{k}\Delta\mathrm{dev}(\overline{S}_{ij})\overline{S}_{ij}=0 \\
\rightarrow \quad
&a\left(\sqrt{k_{_{SGS}}}\right)^{2}+b\sqrt{k_{_{SGS}}}-c=0 \\
\end{align}

よって

\begin{equation}
k_{_{SGS}}=\left(\frac{-b+\sqrt{b^{2}+4ac}}{2a}\right)^{2}
\end{equation}

\begin{equation}
a=\frac{C_{\varepsilon}}{\Delta}, \quad
b=\frac{2}{3}\mathrm{tr}(\overline{S}_{ij}), \quad
c=2C_{k}\Delta\left\{\mathrm{dev}(\overline{S}_{ij})\overline{S}_{ij}\right\}
\end{equation}

Smagorinkyモデルの\(k_{_{SGS}}\)の計算式は,SGSにおける運動エネルギーの生成と散逸の局所平衡を仮定して導出されている

本当にSmagorinskyモデルと同じものなのか

非圧縮性を仮定すると,連続の式より

\begin{equation}
\frac{\partial \overline{u}_{i}}{\partial x_{i}}=0
\quad \rightarrow \quad
\overline{S}_{ii}=\frac{1}{2}\left(\frac{\partial \overline{u}_{i}}{\partial x_{i}}+\frac{\partial \overline{u}_{i}}{\partial x_{i}}\right)=0
\end{equation}

なので

\begin{equation}
\mathrm{dev}(\overline{S}_{ij})=S_{ij}-\frac{1}{3}S_{kk}\delta_{ij}=S_{ij}
\end{equation}

係数\(a\),\(b\),\(c\)について

\begin{align}
\left\{\begin{array}{l}
&a=\frac{C_{\varepsilon}}{\Delta} \\
&b=\frac{2}{3}\mathrm{tr}(\overline{S}_{ij})=0 \\
&c=2C_{k}\Delta\left\{\mathrm{dev}(\overline{S}_{ij})\overline{S}_{ij}\right\}
=C_{k}\Delta\left(2\overline{S}_{ij}\overline{S}_{ij}\right)
\end{array}\right.
\end{align}

これらの係数を\(k_{_{SGS}}\)に代入すると

\begin{equation}
k_{_{SGS}}=\left(\frac{-b+\sqrt{b^{2}+4ac}}{2a}\right)^{2}
=\left(\frac{\sqrt{4ac}}{2a}\right)^{2}=\frac{c}{a}
=\frac{C_{k}}{C_{\varepsilon}}\Delta^{2}\left(2\overline{S}_{ij}\overline{S}_{ij}\right)
\end{equation}

上記の\(k_{_{SGS}}\)をSGS渦粘性係数 \(\nu_{_{SGS}}=C_{k}\Delta\sqrt{k_{_{SGS}}}\)に代入すると

\begin{equation}
\nu_{_{SGS}}=C_{k}\Delta\sqrt{\frac{C_{k}}{C_{\varepsilon}}\Delta^{2}\left(2\overline{S}_{ij}\overline{S}_{ij}\right)}
=C_{k}\sqrt{\frac{C_{k}}{C_{\varepsilon}}}\Delta^{2}\sqrt{2\overline{S}_{ij}\overline{S}_{ij}}
\end{equation}

定数を \(C_{s}^{2}=C_{k}\sqrt{\frac{C_{k}}{C_{\varepsilon}}}\) とすれば,一般的なSmagorinskyモデルの表記

\begin{equation}
\nu_{_{SGS}}=\left(C_{s}\Delta\right)^{2}\sqrt{2\overline{S}_{ij}\overline{S}_{ij}}
\end{equation}

になる

ちなみに,Smagorinsky定数は \(C_{s}=C_{k}^{\frac{3}{4}}C_{\varepsilon}^{-\frac{1}{4}} \simeq 0.167\)である

k Equationモデルの表式で表されたSmagorinskyモデルは,非圧縮性を仮定すると\(C_{s}=0.167\)の一般的なSmagorisnkyモデルと等価である

参考
Smagorinsky SGS model in OpenFOAM _ CFD WITH A MISSION

参考
J. Smagorinsky. General Circulation Experiments with the Primitive Equations I. the Basic Experiment*. Monthly Weather Review, 91(3):99–164, 1963.

OpenFOAMにおける実装

Smagorinskyの実装を確認する

参考
OpenFOAM User Guide Smagorinsky
Smagorinsky.C
Smagorinsky.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_{_{SGS}}=C_{k}\Delta\sqrt{k_{_{SGS}}}
\end{equation}

template<class BasicTurbulenceModel>
void Smagorinsky<BasicTurbulenceModel>::correctNut()
{
    volScalarField k(this->k(fvc::grad(this->U_)));

    this->nut_ = Ck_*this->delta()*sqrt(k);
    this->nut_.correctBoundaryConditions();
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}

\(k_{_{SGS}}\)

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

\begin{equation}
a=\frac{C_{\varepsilon}}{\Delta}, \quad
b=\frac{2}{3}\mathrm{tr}(\overline{S}_{ij}), \quad
c=2C_{k}\Delta\left\{\mathrm{dev}(\overline{S}_{ij})\overline{S}_{ij}\right\}
\end{equation}

\begin{equation}
k_{_{SGS}}=\left(\frac{-b+\sqrt{b^{2}+4ac}}{2a}\right)^{2}
\end{equation}

\(\overline{S}_{ij}\)には\(D_{ij}\)という文字が使われている

template<class BasicTurbulenceModel>
tmp<volScalarField> Smagorinsky<BasicTurbulenceModel>::k
(
    const tmp<volTensorField>& gradU
) const
{
    volSymmTensorField D(symm(gradU));

    volScalarField a(this->Ce_/this->delta());
    volScalarField b((2.0/3.0)*tr(D));
    volScalarField c(2*Ck_*this->delta()*(dev(D) && D));

    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName("k", this->alphaRhoPhi_.group()),
                this->runTime_.timeName(),
                this->mesh_
            ),
            sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a))
        )
    );
}

定数群

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

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

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

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

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

まとめ

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

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

コメント