OpenFOAMで実装されている1方程式渦粘性モデルであるSpalartAllmarasについて説明する
keywords: OpenFOAM，乱流モデル，SpalartAllmaras
Spalart-Allmarasモデル
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}\)をどのように求めるかが渦粘性モデルの唯一にして最も重要な点である
Spalart-Allmarasモデルにおいて，渦粘性係数は次の式で与えられる
\begin{equation}
\nu_{t}=\tilde{\nu}f_{v_{1}}, \qquad
f_{v_{1}}=\frac{\chi^{3}}{\chi^{3}+C_{V_{1}}^{3}}, \qquad
\chi=\frac{\tilde{\nu}}{\nu}
\end{equation}
\(\tilde{\nu}\)は次の輸送方程式によって求められる
\begin{equation}
\frac{\partial \tilde{\nu}}{\partial t}
+\overline{u}_{i}\frac{\partial \tilde{\nu}}{\partial x_{i}}
=
C_{b_{1}}\tilde{S}\tilde{\nu}
-C_{w_{1}}f_{w}\left(\frac{\tilde{\nu}}{y}\right)^{2}
+\frac{1}{\sigma_{\nu_{t}}}\left\{\frac{\partial}{\partial x_{j}}
\left(\tilde{\nu}+\nu\right)\frac{\partial \tilde{\nu}}{\partial x_{j}}
+C_{b_{2}}\frac{\partial \tilde{\nu}}{\partial x_{j}}\frac{\partial \tilde{\nu}}{\partial x_{j}}\right\}
\end{equation}
ここで，\(y\)は壁面からの距離である
\(\tilde{S}\)は次の式で与えられる
\begin{equation}
\tilde{S}=\mathrm{max}\left[
\sqrt{2}\Omega+f_{v_{2}}\frac{\nu}{\left(\kappa d\right)^{2}}, \,
\sqrt{2}C_{S}\Omega
\right]
\end{equation}
\begin{equation}
\Omega=\sqrt{\Omega_{ij}\Omega_{ij}}, \qquad
\Omega_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}-\frac{\partial u_{j}}{\partial x_{i}}\right)
\end{equation}
\begin{equation}
f_{v_{2}}=1-\frac{\chi}{1.0+\chi f_{v_{1}}}
\end{equation}
\(f_{w}\)は次の式で与えられる
\begin{equation}
f_{w}=g\left(\frac{1+C_{w_{3}}^{6}}{g^{6}+C_{w_{3}}^{6}}\right)^{\frac{1}{6}}
\end{equation}
\begin{equation}
g=r+C_{w_{2}}\left(r^{6}-r\right), \qquad
r=\mathrm{min}\left[\frac{\tilde{\nu}}{\tilde{S}\left(\kappa d\right)^{2}},\,10\right]
\end{equation}
モデル定数は次の通り
\begin{align}
&\sigma_{\nu_{t}}=2/3, \quad
\kappa=0.41, \quad
C_{b_{1}}=0.1355, \quad
C_{b_{2}}=0.622, \quad \\
&C_{w_{1}}=\frac{C_{b_{1}}}{\kappa^{2}}+\frac{1+C_{b_{2}}}{\sigma_{\nu_{t}}}, \quad
C_{w_{2}}=0.3, \quad
C_{w_{3}}=2, \quad
C_{V_{1}}=7.1, \quad
C_{S}=0.3
\end{align}
OpenFOAMにおける実装
SpalartAllmarasの実装を確認する
参考
≫SpalartAllmaras< BasicTurbulenceModel > Class Template Reference
≫SpalartAllmaras.C
≫SpalartAllmaras.H
ただし，OpenFOAMで実装されているモデルは論文中のモデルから以下の変更を加えている
- trip-termを実装していない
- \(\tilde{S}\)の非物理的な挙動を防ぐために制限を課している
The model is implemented without the trip-term and hence the ft2 term is not needed.
It is necessary to limit the Stilda generation term as the model generates unphysical results if this term becomes negative which occurs for complex flow. Several approaches have been proposed to limit Stilda but it is not clear which is the most appropriate. Here the limiter proposed by Spalart is implemented in which Stilda is clipped at Cs*Omega with the default value of Cs = 0.3.
SpalartAllmaras< BasicTurbulenceModel > Class Template Reference
実際に計算が行われるときは，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}
\chi=\frac{\tilde{\nu}}{\nu}
\end{equation}
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::chi() const
{
return nuTilda_/this->nu();
}\begin{equation}
f_{v_{1}}=\frac{\chi^{3}}{\chi^{3}+C_{V_{1}}^{3}}
\end{equation}
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fv1
(
const volScalarField& chi
) const
{
const volScalarField chi3(pow3(chi));
return chi3/(chi3 + pow3(Cv1_));
}\begin{equation}
\nu_{t}=\tilde{\nu}f_{v_{1}}
\end{equation}
template<class BasicTurbulenceModel>
void SpalartAllmaras<BasicTurbulenceModel>::correctNut
(
const volScalarField& fv1
)
{
this->nut_ = nuTilda_*fv1;
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);
BasicTurbulenceModel::correctNut();
}
template<class BasicTurbulenceModel>
void SpalartAllmaras<BasicTurbulenceModel>::correctNut()
{
correctNut(fv1(this->chi()));
}輸送方程式
\begin{equation}
\frac{\partial \tilde{\nu}}{\partial t}
+\overline{u}_{i}\frac{\partial \tilde{\nu}}{\partial x_{i}}
=
C_{b_{1}}\tilde{S}\tilde{\nu}
-C_{w_{1}}f_{w}\left(\frac{\tilde{\nu}}{y}\right)^{2}
+\frac{1}{\sigma_{\nu_{t}}}\left\{\frac{\partial}{\partial x_{j}}
\left(\tilde{\nu}+\nu\right)\frac{\partial \tilde{\nu}}{\partial x_{j}}
+C_{b_{2}}\frac{\partial \tilde{\nu}}{\partial x_{j}}\frac{\partial \tilde{\nu}}{\partial x_{j}}\right\}
\end{equation}
整理して
\begin{equation}
\frac{\partial \tilde{\nu}}{\partial t}
+\overline{u}_{i}\frac{\partial \tilde{\nu}}{\partial x_{i}}
-\frac{\partial}{\partial x_{j}}\left\{
\frac{\partial }{\partial x_{j}}\left(
\frac{\tilde{\nu}+\nu}{\sigma_{\nu_{t}}}\tilde{\nu}
\right)\right\}
-\frac{C_{b_{2}}}{\sigma_{\nu_{t}}}\left(\frac{\partial \tilde{\nu}}{\partial x_{j}}\right)^{2}
=
C_{b_{1}}\tilde{S}\tilde{\nu}
-C_{w_{1}}f_{w}\left(\frac{\tilde{\nu}}{y}\right)^{2}
\end{equation}
tmp<fvScalarMatrix> nuTildaEqn
(
fvm::ddt(alpha, rho, nuTilda_)
+ fvm::div(alphaRhoPhi, nuTilda_)
- fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
- Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
==
Cb1_*alpha*rho*Stilda*nuTilda_
- fvm::Sp(Cw1_*alpha*rho*fw(Stilda)*nuTilda_/sqr(y_), nuTilda_)
+ fvOptions(alpha, rho, nuTilda_)
);
nuTildaEqn.ref().relax();
fvOptions.constrain(nuTildaEqn.ref());
solve(nuTildaEqn);
fvOptions.correct(nuTilda_);
bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
nuTilda_.correctBoundaryConditions();template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::DnuTildaEff() const
{
return tmp<volScalarField>
(
new volScalarField("DnuTildaEff", (nuTilda_ + this->nu())/sigmaNut_)
);
}生成項（Production term）
\begin{equation}
f_{v_{2}}=1-\frac{\chi}{1.0+\chi f_{v_{1}}}
\end{equation}
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fv2
(
const volScalarField& chi,
const volScalarField& fv1
) const
{
return 1.0 - chi/(1.0 + chi*fv1);
}\begin{equation}
\Omega=\sqrt{2\Omega_{ij}\Omega_{ij}}, \qquad
\Omega_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}-\frac{\partial u_{j}}{\partial x_{i}}\right)
\end{equation}
\begin{equation}
\tilde{S}=\mathrm{max}\left[
\Omega+f_{v_{2}}\frac{\nu}{\left(\kappa d\right)^{2}}, \,
C_{S}\Omega
\right]
\end{equation}
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::Stilda
(
const volScalarField& chi,
const volScalarField& fv1
) const
{
volScalarField Omega(::sqrt(2.0)*mag(skew(fvc::grad(this->U_))));
return
(
max
(
Omega
+ fv2(chi, fv1)*nuTilda_/sqr(kappa_*y_),
Cs_*Omega
)
);
}消滅項（Destruction term）
\begin{equation}
r=\mathrm{min}\left[\frac{\tilde{\nu}}{\tilde{S}\left(\kappa d\right)^{2}},\,10\right]
\end{equation}
\begin{equation}
g=r+C_{w_{2}}\left(r^{6}-r\right)
\end{equation}
\begin{equation}
f_{w}=g\left(\frac{1+C_{w_{3}}^{6}}{g^{6}+C_{w_{3}}^{6}}\right)^{\frac{1}{6}}
\end{equation}
template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fw
(
const volScalarField& Stilda
) const
{
volScalarField r
(
min
(
nuTilda_
/(
max
(
Stilda,
dimensionedScalar("SMALL", Stilda.dimensions(), SMALL)
)
*sqr(kappa_*y_)
),
scalar(10)
)
);
r.boundaryFieldRef() == 0.0;
const volScalarField g(r + Cw2_*(pow6(r) - r));
return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
}定数群
\begin{align}
&\sigma_{\nu_{t}}=2/3, \quad
\kappa=0.41, \quad
C_{b_{1}}=0.1355, \quad
C_{b_{2}}=0.622, \quad \\
&C_{w_{1}}=\frac{C_{b_{1}}}{\kappa^{2}}+\frac{1+C_{b_{2}}}{\sigma_{\nu_{t}}}, \quad
C_{w_{2}}=0.3, \quad
C_{w_{3}}=2, \quad
C_{V_{1}}=7.1, \quad
C_{S}=0.3
\end{align}
sigmaNut_
(
dimensioned<scalar>::getOrAddToDict
(
"sigmaNut",
this->coeffDict_,
0.66666
)
),
kappa_
(
dimensioned<scalar>::getOrAddToDict
(
"kappa",
this->coeffDict_,
0.41
)
),
Cb1_
(
dimensioned<scalar>::getOrAddToDict
(
"Cb1",
this->coeffDict_,
0.1355
)
),
Cb2_
(
dimensioned<scalar>::getOrAddToDict
(
"Cb2",
this->coeffDict_,
0.622
)
),
Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
Cw2_
(
dimensioned<scalar>::getOrAddToDict
(
"Cw2",
this->coeffDict_,
0.3
)
),
Cw3_
(
dimensioned<scalar>::getOrAddToDict
(
"Cw3",
this->coeffDict_,
2.0
)
),
Cv1_
(
dimensioned<scalar>::getOrAddToDict
(
"Cv1",
this->coeffDict_,
7.1
)
),
Cs_
(
dimensioned<scalar>::getOrAddToDict
(
"Cs",
this->coeffDict_,
0.3
)
),
境界条件
速度と圧力の壁面境界条件はすべりなし条件 \(u_{i}=0\)，\(\frac{\partial p}{\partial n}=0\) を設定する
\(\tilde{\nu}\)の境界条件は壁面境界条件として \(\tilde{\nu}=0\) を設定する
自由流境界としては\(\tilde{\nu}=0\)が理想だが，以下の条件でもよい
\begin{equation}
\tilde{\nu}<\frac{\nu}{10}
\end{equation}
初期条件は自由流境界と同様の基準で決定する
The wall boundary condition is \(\tilde{\nu}=0\). In the freestream 0 is best, provided numerical errors do not \(\tilde{\nu}\) to negative values near the edge of the boundary layer (the exact solution cannot go negative). Value below \(\nu/10\) will be acceptable. The same applies to the initial condition.
A One-Equation Turbulence Model for Aerodynamic Flows," Recherche Aerospatiale pp.20-21
まとめ
OpenFOAMで実装されている1方程式渦粘性モデルであるSpalartAllmarasについて説明した
↓関連記事
コメント