kOmega 是湍流模型中应用广泛的一个模型。在上一节 Code of the Week,做了 kEpsilon 模型的初步讲解。本次,对 kOmega 作初步的解析。
k 方程已经在 Code of the Week #14: RASModel kEpsilon 作了些解释。现在要做的是选择另外一个参数-长度,Kolmogorov 选择第二个参数为 \omega ,单位湍动能的耗散率,其方程与 k 极其相似,虽然没有什么物理意义,Kolmogorov 还把 \omega 称为单位时间和体积的能量的耗散率;Kolmogorov 的 \omega 方程为
\rho\frac{\partial\omega}{\partial t}+\rho U_j\frac{\partial\omega}{\partial x_j}
=
-\beta\rho\omega^2 + \frac{\partial}{\partial x_j}\left(\sigma\mu_t\frac{\partial\omega}{\partial x_j} \right)
后来,经过 Wilcox and Alber, Spalding 等不懈努力,形成了现今的 k-\omega 模型,Wilcox 的模型如下:
\rho\frac{\partial k}{\partial t} + \rho U_j\frac{\partial k}{\partial x_j}
=
\tau_{ij}\frac{\partial U_i}{\partial x_j} - C_\mu \rho k \omega
+\frac{\partial}{\partial x_j}\left[(\mu +\alpha_k \mu_t)\frac{\partial k}{\partial x_j} \right]
\omega 方程
\rho\frac{\partial\omega}{\partial t} + \rho U_j\frac{\partial \omega}{\partial x_j}
=
\alpha\frac{\omega}{k}\tau_{ij}\frac{\partial U_i}{\partial x_j} - \beta \rho \omega^2
+\frac{\partial}{\partial x_j}\left[(\mu + \alpha_\omega \mu_t)\frac{\partial\omega}{\partial x_j} \right]
其中, \mu_t=\rho k/\omega, \varepsilon=C_\mu \omega k, l=\dfrac{k^{1/2}}{\omega} ,其他参数为
\alpha = 5/9,\, C_\mu=0.009,\, \beta = 3/40,\,\alpha_k=0.5,\,\alpha_\omega=0.5
在 kOmega 模型中,照例,最关键的函数是 correct(),位于 kOmega.C (foam-extend-3.1 版本) 中
void kOmega::correct()
{
// Bound in case of topological change
// HJ, 22/Aug/2007
if (mesh_.changing())
{
bound(k_, k0_);
bound(omega_, omega0_);
}
RASModel::correct();
if (!turbulence_)
{
return;
}
volScalarField G("RASModel::G", nut_*2*magSqr(symm(fvc::grad(U_)))); // tau_{ij} 雷诺应力
// Update omega and G at the wall
omega_.boundaryField().updateCoeffs();
// Turbulence specific dissipation rate equation
tmp<fvScalarMatrix> omegaEqn // omega 方程
(
fvm::ddt(omega_) // 左边第 1 项
+ fvm::div(phi_, omega_) // 左边第 2 项
+ fvm::SuSp(-fvc::div(phi_), omega_)
- fvm::laplacian(DomegaEff(), omega_) // omega 方程右边第 3 项
==
alpha_*G*omega_/k_ // 右边第 1 项
- fvm::Sp(beta_*omega_, omega_) // 右边第 2 项
);
omegaEqn().relax(); // 对方程进行松弛操作
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// omegaEqn().completeAssembly();
solve(omegaEqn);
bound(omega_, omega0_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
+ fvm::SuSp(-fvc::div(phi_), k_)
- fvm::laplacian(DkEff(), k_)
==
G
- fvm::Sp(Cmu_*omega_, k_)
);
kEqn().relax();
solve(kEqn);
bound(k_, k0_);
// Re-calculate viscosity
nut_ = k_/(omega_ + omegaSmall_);
nut_.correctBoundaryConditions();
}
DkEff()
为湍动能耗散系数,为 k 方程中的 \mu +\alpha_k \mu_t
tmp<volScalarField> DkEff() const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", alphaK_*nut_ + nu())
);
}
DomegaEff()
为 \omega 耗散系数,为 \omega 方程中的 \mu +\alpha_\omega \mu_t
tmp<volScalarField> DomegaEff() const
{
return tmp<volScalarField>
(
new volScalarField("DomegaEff", alphaOmega_*nut_ + nu())
);
}
参考文献
- Wilcox, D C. Turbulence Modeling for CFD.