OpenFOAMで実装されている低Re数型k-εモデルであるLamBremhorstKEについて説明する
keywords: OpenFOAM，乱流モデル，LamBremhorstKE
低レイノルズ数型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\) として計算される
LamBremhorstKE
OpenFOAMで実装されている LamBremhorstKEは次のようなモデルである
\begin{equation}
f_{\mu}=\left\{1-\exp{\left(-0.0165R_{y}\right)}\right\}^{2}\left(1+\frac{20.5}{R_{t}}\right)
\end{equation}
\begin{equation}
f_{1}=1+\left(\frac{0.05}{f_{\mu}}\right)^{3}, \qquad
f_{2}=1-\exp{\left(-R_{t}^{2}\right)}
\end{equation}
\begin{equation}
D=0, \qquad
E=0
\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における実装
LamBremhorstKEの実装を確認する
参考
≫LamBremhorstKE Class Reference
≫LamBremhorstKE.C
≫LamBremhorstKE.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}
void LamBremhorstKE::correctNut(const volScalarField& fMu)
{
nut_ = Cmu_*fMu*sqr(k_)/epsilon_;
nut_.correctBoundaryConditions();
}輸送方程式
\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(GName(), nut_*(twoSymm(tgradU()) && 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}}f_{1}P_{k}-C_{\varepsilon_{2}}f_{2}\tilde{\varepsilon}\right)\frac{\tilde{\varepsilon}}{k}+\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}-\tilde{\varepsilon}
+\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();
const volScalarField Rt(this->Rt());
const volScalarField fMu(this->fMu(Rt));
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
Ceps1_*f1(fMu)*G*epsilon_/k_
- fvm::Sp(Ceps2_*f2(Rt)*epsilon_/k_, epsilon_)
);
epsEqn.ref().relax();
epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
solve(epsEqn);
bound(epsilon_, epsilonMin_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::laplacian(DkEff(), k_)
==
G - fvm::Sp(epsilon_/k_, k_)
);
kEqn.ref().relax();
solve(kEqn);
bound(k_, kMin_); //- Return the effective diffusivity for k
tmp<volScalarField> DkEff() const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", nut_ + nu())
);
}
//- Return the effective diffusivity for epsilon
tmp<volScalarField> DepsilonEff() const
{
return tmp<volScalarField>
(
new volScalarField("DepsilonEff", nut_/sigmaEps_ + nu())
);
}補正項
\begin{equation}
R_{t}=\frac{k^2}{\nu\varepsilon}
\end{equation}
tmp<volScalarField> LamBremhorstKE::Rt() const
{
return sqr(k_)/(nu()*epsilon_);
}\begin{equation}
f_{\mu}=\left\{1-\exp{\left(-0.0165R_{y}\right)}\right\}^{2}\left(1+\frac{20.5}{R_{t}}\right)
\end{equation}
tmp<volScalarField> LamBremhorstKE::fMu(const volScalarField& Rt) const
{
tmp<volScalarField> Ry(sqrt(k_)*y_/nu());
return sqr(scalar(1) - exp(-0.0165*Ry))*(scalar(1) + 20.5/(Rt + SMALL));
}\begin{equation}
f_{1}=1+\left(\frac{0.05}{f_{\mu}}\right)^{3}
\end{equation}
tmp<volScalarField> LamBremhorstKE::f1(const volScalarField& fMu) const
{
return scalar(1) + pow3(0.05/(fMu + SMALL));
}\begin{equation}
f_{2}=1-\exp{\left(-R_{t}^{2}\right)}
\end{equation}
tmp<volScalarField> LamBremhorstKE::f2(const volScalarField& Rt) const
{
return scalar(1) - exp(-sqr(Rt));
}定数群
\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",
coeffDict_,
0.09
)
),
Ceps1_
(
dimensioned<scalar>::getOrAddToDict
(
"Ceps1",
coeffDict_,
1.44
)
),
Ceps2_
(
dimensioned<scalar>::getOrAddToDict
(
"Ceps2",
coeffDict_,
1.92
)
),
sigmaEps_
(
dimensioned<scalar>::getOrAddToDict
(
"alphaEps",
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}=\frac{2\nu k}{y^{2}}
\end{equation}
（渦粘性係数）
\begin{equation}
\nu_{t_{w}}=0
\end{equation}
まとめ
OpenFOAMで実装されている低Re数型k-εモデルであるLamBremhorstKEについて説明した
そのほかの乱流モデルについてはこちら
低レイノルズ数型k-εモデルについてもう少し詳しく知りたい人はこちら
コメント