Code of the Week #15: RASModel kOmega

OpenFOAM
1

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())
            );
        }

参考文献

  1. Wilcox, D C. Turbulence Modeling for CFD.
2
3