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モデルと等価である
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ついて説明した
そのほかの乱流モデルについてはこちら
コメント