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