• 返回主页
  • OpenFOAM中的SpalartAllmaras模型

    作者:邱小平
    编辑:Y.M.Z
    本篇简要分析不可压缩的SpalartAllmaras模型的代码。主要内容包括模型输运方程的代码说明,以及一些使用方面的细节。

    模型分析

    SpalartAllmaras模型是一方程模型,它只需要求解一个输运方程。OpenFOAM中的SpalartAllmaras模型是在Spalart and Allmaras的论文的基础上,引入了Ashford对这个模型的修正。下面来看这个模型在OpenFOAM中的实现,代码位于src/turbulenceModels/incompressible/RAS/SpalartAllmaras下。

    首先是头文件SpalartAllmaras.H。头文件开始,声明了一个类:SpalartAllmaras

    1
    2
    3
    class SpalartAllmaras
    :
        public RASModel
    

    说明这个类是公有继承RASModel类的。接下来是一些数据成员以及成员函数的声明,这里无需赘述,直接来看函数的具体定义吧。函数的具体定义,在SpalartAllmaras.C里(注意到 OpenFOAM的代码几乎都是这样讲声明和定义分开,这样做的目的,是为了防止重复声明的问题,更详细的说明,见我的另一篇博文)。注意,湍流模型的定义,最终是为了在求解器中进行调用的,所以,先来回顾一下湍流模型贡献给求解器中动量方程中那一项:

    1
    turbulence->divDevReff(U)
    

    这一项,不同求解器可能会有不同,但大致的形式是这样的。很显然,这个divDevReff应该是湍流模型类的成员函数,实际也的确如此,SpalartAllmaras.C有这个函数的定义:

    1
    2
    3
    4
    5
    6
    7
    8
    tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
    {
        return
        (
          - fvm::laplacian(nuEff(), U)
          - fvc::div(nuEff()*dev(T(fvc::grad(U))))
        );
    }
    

    这个函数,写成公式是如下形式 \begin{equation} -\nabla \cdot (\nu_{eff}\nabla \mathbf{U})-\nabla \cdot [\nu_{eff}\dev(\nabla ^\mathrm{T}\bf{U})] \end{equation} OpenFOAM中dev运算符定义为$\dev(\bfU)=\bfU-\frac{1}{3}\tr(\bfU)\bfI$,$\tr$为张量的迹,即主对角元素之和所以。对于张量$\nabla ^\mathrm{T}\bfU$,有$\tr(\nabla ^\mathrm{T}\bfU) =\nabla \cdot \bfU$,因此上式可以写为 \begin{equation} -\nabla \cdot (\nu_{eff}\nabla \bfU)-\nabla \cdot [\nu_{eff}(\nabla^\rT\bfU-\frac{1}{3}(\nabla \cdot \bfU))] \end{equation} 上述divDevReff函数中,变量nuEff_d的值是函数nuEff()的返回值,但是,在SpalartAllmaras.C中却未见nuEff()函数的定义。这里就要记住了,C++ 的类继承的一个特性是派生类会从基类继承其中定义的成员函数。所以看到本类中没定义的函数时,don’t panic,往基类找就是了。翻看基类RASModel的代码,可以很容易从RASModel.C 中找出来nuEff()的定义:

    1
    2
    3
    4
    5
    6
    7
    virtual tmp<volScalarField> nuEff() const  
    {
        return tmp<volScalarField>
        (
            new volScalarField("nuEff", nut() + nu())
        );
    }
    

    函数nut()的定义可以从SpalartAllmaras.H中找到,可是nu()呢?虽然从名字可以猜到它返回的是分子粘度(可见变量命名的重要性啊),可是,它的定义在哪里?RASModel.C也没有啊… 答案是,继续向上找,RASModel类没有,那就去更上一层的基类,turbulenceModel。果然,在turbulenceModel.C中能找到其定义:

    1
    2
    3
    4
    inline tmp<volScalarField> nu() const
        {
            return transportModel_.nu();
        }
    

    transportModel_又是什么鬼?这个不是本篇博文的内容,该回到正题了,这里应该关注的是湍流粘度nut_ 。湍流模型不可能仅仅是给求解器贡献那一项,湍流粘度应该也随着求解器的运行而更新。是这样的,求解器中的

    1
    turbulence->correct();
    

    这句代码负责湍流粘度的更新。这里又涉及到一个湍流类的成员函数: correct。这个函数是湍流模型类的核心,因为湍流模型输运方程的求解,湍流粘度的更新都在这个函数里面。下面来看SpalartAllmaras模型的correct函数:

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    void SpalartAllmaras::correct()
    {
        RASModel::correct();
        if (!turbulence_)
        {
            // Re-calculate viscosity
            nut_ = nuTilda_*fv1(this->chi());
            nut_.correctBoundaryConditions();
    
            return;
        }
        if (mesh_.changing())
        {
            d_.correct();
        }
    
        const volScalarField chi(this->chi());
        const volScalarField fv1(this->fv1(chi));
    
        const volScalarField Stilda
        (
            fv3(chi, fv1)*::sqrt(2.0)*mag(skew(fvc::grad(U_)))
          + fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_)
        );
    
        tmp<fvScalarMatrix> nuTildaEqn
        (
            fvm::ddt(nuTilda_)
          + fvm::div(phi_, nuTilda_)
          - fvm::laplacian(DnuTildaEff(), nuTilda_)
          - Cb2_/sigmaNut_*magSqr(fvc::grad(nuTilda_))
         ==
            Cb1_*Stilda*nuTilda_
          - fvm::Sp(Cw1_*fw(Stilda)*nuTilda_/sqr(d_), nuTilda_)
        );
    
        nuTildaEqn().relax();
        solve(nuTildaEqn);
        bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
        nuTilda_.correctBoundaryConditions();
    
        // Re-calculate viscosity
        nut_.internalField() = fv1*nuTilda_.internalField();
        nut_.correctBoundaryConditions();
    }
    

    先看核心的输运方程,nuTildaEqn,这个部分的每一行写成公式,分别如下: \begin{equation} \frac{\partial \tilde{\nu}}{\partial t} \end{equation} \begin{equation} \nabla \cdot (\bfU \tilde{\nu}) \end{equation} \begin{equation} -\nabla \cdot [D_{\tilde{\nu}_{eff}} \nabla \tilde{\nu}] \end{equation} \begin{equation} -\frac{C_{b2}}{\sigma_{\nu_t}} \cdot |\nabla \tilde{\nu}|^2 \end{equation} \begin{equation} C_{b1}\cdot \tilde{S} \cdot \tilde{\nu} \end{equation} \begin{equation} -\frac{C_{w1}\cdot f_w\cdot \tilde{\nu}}{d^2} \cdot \tilde{\nu} \end{equation} 最后两项都包含了需要求解的量$\tilde{\nu}$, 但是处理方式却不一样,最后一项,其实是将原本的$\tilde{\nu}^2$项进行了线性化,并用fvm::Sp操作符进行隐式处理。而倒数第二项,则是作的显式处理。差别在于,倒数第二项仅仅是用上一步得到的$\tilde{\nu}^2$代入,然后这一项被加到$Ax=b$中的$b$ 部分。而隐式处理则会对系数矩阵$A$有贡献。

    输运方程中涉及到的变量和函数,全部都在SpalartAllmaras类中有定义,这里就不再单独分析了,仅将涉及到的公式列出如下 \begin{equation} D_{\tilde{\nu}_{eff}}=\frac{\nu+\tilde{\nu}}{\sigma_{\nu_t}} \end{equation} \begin{equation} f_w=g\cdot \left [ \frac{1+C_{w3}^6}{g^6+C_{w3}^6} \right]^{1/6},\quad g=r+C_{w2}\cdot(r^6-r) \end{equation} \begin{equation} r=\min\left[ \frac{\tilde{\nu}}{\max(\tilde{S}, \mathrm{SMALL})\cdot \kappa^2 d^2},10 \right] \end{equation} \begin{equation} \tilde{S}=f_{v3}\cdot\sqrt{2\Omega_{ij}:\Omega_{ij}}+\frac{f_{v2}\tilde{\nu}}{\kappa^2 d^2},\quad \Omega_{ij}=\frac{1}{2}(\nabla U -\nabla ^T U) \end{equation} $f_{v2}$和$f_{v3}$有点特殊,因为涉及到了Ashford对原始模型的修正。根据作者的论文,修正的目的是防止$\tilde{S}$ 出现负值。如果不引入Ashford修正,在$f_{v3}=1$引入以后,则 \begin{equation}f_{v3}=\frac{(1+\chi \cdot f_{v1})}{c_{v2}}\cdot \frac{3+\frac{\chi}{c_{v2}} + (\frac{\chi}{c_{v2}})^2 }{(1+\frac{\chi}{c_{v2}})^3} \end{equation} 对于$f_{v2}$,如果不引入Ashford修正,则 \begin{equation}f_{v2}=\frac{1}{(1+\frac{\chi}{c_{v2}})^3}\end{equation} 若引入,则 \begin{equation}f_{v2}=1-\frac{\chi}{1+\chi f_{v1}}\end{equation} 此外, \begin{equation} f_{v1}=\frac{\chi ^3}{\chi ^3 + c_{v1}^3}, \quad \chi=\frac{\tilde{\nu}}{\nu}, \quad \nu_t = \tilde{\nu} f_{v1} \end{equation} $d$是与壁面的距离。

    一些细节

  • correct函数中,只有当turbulence_true 时,下面求解输运方程的部分代码才会执行。这个变量又是从RASModel类继承过来的,其初始化的语句为turbulence_(lookup("turbulence")),这意味着,这个变量的值是从RASProperties文件中读取的(详见文末的 RASProperties 示例)。

  • 模型常数除了使用默认值,还可以自己指定(注:除非很有把握,一般不要随便修改模型常数)。具体方法是,在RASProperties文件里,建立一个名为SpalartAllmarasCoeffs 的subDict。

  • bound 函数的定义在头文件bound.H 里,其功能是限制$\tilde{\nu}$的最小值。不过这种限制其实没有多少物理意义,纯粹是数值技巧。其中用到了fvc::average 函数(见OpenFOAM中的div与snGrad操作符)。

  • read 成员函数初看之下似乎跟模型常数的初始化有关,实则不然。模型常数的初始化在类的构造函数部分完成,coeffDict_这个变量在基类RASModel的构造函数中先用type()+Coeffs (这里是 SpalartAllmarasCoeffs )这个subDict初始化了,然后,在SpalartAllmaras类的构造函数中,每一个模型常数都是先去coeffDict_中查找是否用用户自定义的值,如果没有,则用模型的默认值,并且将这个默认值添加到coeffDict_这个字典中。

  • 对于其他成员函数,SpalartAllmaras显然不需要用到湍动能k以及湍动能耗散epsilon ,但是在SpalartAllmaras.C中确有这两个函数的定义,而且返回值是0。这是因为,在基类 turbulenceModel 中,k()epsilon()函数被声明为纯虚函数,如果不在SpalartAllmaras中重新定义这两个函数,则SpalartAllmaras也将是包括纯虚函数的抽象类了。在 C++ 中,抽象类是无法建立对象的。devReffRdivDevRhoReff这几个成员函数,是湍流模型提供给其他模型调用的函数,比如,在计算表面应力的forces类中就需要调用湍流模型的devReff函数。

  • 东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国