OpenFOAMで実装されている低Re数型k-εモデルであるLaunderSharmaKEについて説明する
keywords: OpenFOAM，乱流モデル，LaunderSharmaKE
低レイノルズ数型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}\)をどのように求めるかが渦粘性モデルの唯一にして最も重要な点である
一般的な低レイノルズ数型モデルは次のように表される
（渦粘性係数）
\begin{equation}
\nu_{t}=C_{\mu}f_{\mu}\frac{k^{2}}{\varepsilon}
\end{equation}
（\(k\)の輸送方程式）
\begin{equation}
\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{equation}
（\(\varepsilon\)の輸送方程式）
\begin{equation}
\frac{\partial \tilde{\varepsilon}}{\partial t}+\overline{u}_{j}\frac{\partial \tilde{\varepsilon}}{\partial x_{j}}
=\left(C_{\varepsilon_{1}}f_{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\}
\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}
\(f_{\mu}\)は\(\nu_{t}\)に対する減衰関数，\(f_{1}\)および\(f_{2}\)は\(\varepsilon\)の生産項および消滅項における補正関数，\(D\)および\(E\)は壁面補正項である
粘性散逸率の真の値は壁面補正項\(D\)を用いて \(\varepsilon=\tilde{\varepsilon}+D\) として計算される
LaunderSharmaKE
OpenFOAMで実装されているLaunderSharmaKEは次のようなモデルである
\begin{equation}
f_{\mu}=\exp{\left[\frac{-3.4}{\left\{1+\left(\frac{k^2}{\nu\varepsilon}/50\right)\right\}^{2}}\right]}, \qquad
f_{1}=1, \qquad
f_{2}=1-0.3\exp\left\{-\left(\frac{k^2}{\nu\varepsilon}\right)^{2}\right\}
\end{equation}
\begin{equation}
D=2\nu\left(\frac{\partial \sqrt{k}}{\partial t}\right)^{2}, \qquad
E=2\nu\nu_{t}\left(\frac{\partial^{2} \overline{u}}{\partial y^{2}}\right)^{2}
\end{equation}
\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}
OpenFOAMにおける実装
LaunderSharmaKEの実装を確認する
参考
≫LaunderSharmaKE< BasicTurbulenceModel > Class Template Reference
≫LaunderSharmaKE.C
≫LaunderSharmaKE.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 LaunderSharmaKE<BasicTurbulenceModel>::correctNut()
{
this->nut_ = Cmu_*fMu()*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);
volScalarField G(this->GName(), nut*(tgradU() && dev(twoSymm(tgradU()))));
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}
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(alpha, rho, epsilon_)
+ fvm::div(alphaRhoPhi, epsilon_)
- fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
==
C1_*alpha*rho*G*epsilon_/k_
- fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha*rho*divU, epsilon_)
- fvm::Sp(C2_*f2()*alpha*rho*epsilon_/k_, epsilon_)
+ alpha*rho*E
+ 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_ + D)/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}
f_{\mu}=\exp{\left[\frac{-3.4}{\left\{1+\left(\frac{k^2}{\nu\varepsilon}/50\right)\right\}^{2}}\right]}
\end{equation}
template<class BasicTurbulenceModel>
tmp<volScalarField> LaunderSharmaKE<BasicTurbulenceModel>::fMu() const
{
return exp(-3.4/sqr(scalar(1) + sqr(k_)/(this->nu()*epsilon_)/50.0));
}\begin{equation}
f_{2}=1-0.3\exp\left\{-\left(\frac{k^2}{\nu\varepsilon}\right)^{2}\right\}
\end{equation}
template<class BasicTurbulenceModel>
tmp<volScalarField> LaunderSharmaKE<BasicTurbulenceModel>::f2() const
{
return
scalar(1)
- 0.3*exp(-min(sqr(sqr(k_)/(this->nu()*epsilon_)), scalar(50)));
}\begin{align}
E&=2\nu\nu_{t}\left(\frac{\partial^{2} \overline{u}}{\partial y^{2}}\right)^{2} \\
\\
D&=2\nu\left(\frac{\partial \sqrt{k}}{\partial t}\right)^{2}
\end{align}
volScalarField E(2.0*this->nu()*nut*fvc::magSqrGradGrad(U));
volScalarField D(2.0*this->nu()*magSqr(fvc::grad(sqrt(k_))));定数群
\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
)
),壁面境界条件
低Re数型モデルを使うときは，第1格子点を粘性低層（\(y^{+} \leq 2 \sim 4\)）に配置して， 滑りなしの壁面境界条件を計算する
（速度）
\begin{equation}
u_{i}=0
\end{equation}
（圧力）
\begin{equation}
\frac{\partial p}{\partial n}=0
\end{equation}
（乱流エネルギー）
\begin{equation}
k_{w}=0
\end{equation}
（乱流散逸率）
\begin{equation}
\tilde{\varepsilon}_{w}=0
\end{equation}
（渦粘性係数）
\begin{equation}
\nu_{t_{w}}=0
\end{equation}
Launder-Sharmaモデルは補正項として\(D\)項を用いているので， \(\varepsilon_{w}=\tilde{\varepsilon}_{w}+D=\frac{2\nu k}{y^{2}}\)より \(\tilde{\varepsilon}\)の壁面境界条件として \(\tilde{\varepsilon}_{w}=0\) を設定する必要がある
まとめ
OpenFOAMで実装されている低Re数型k-εモデルであるLaunderSharmaKEについて説明した
そのほかの乱流モデルについてはこちら
低レイノルズ数型k-εモデルについてもう少し詳しく知りたい人はこちら
コメント