OpenFOAM

OpenFOAM 中的单相流湍流模型之SpalartAllmaras

本篇简要分析不可压缩的 SpalartAllmaras 模型的代码。主要内容包括模型输运方程的代码说明，以及一些使用方面的细节。

湍流模型代码实例

这一部分分析几个 OpenFOAM 里自带的湍流模型， 并探讨修改或者添加新湍流模型的方法。

1 SpalartAllmaras 模型

1.1. 模型分析

SpalartAllmaras 模型是一方程模型，它只需要求解一个输运方程。OpenFOAM 中的 SpalartAllmaras 模型是在原始论文[1] 的基础上，引入了 Ashford [2] 对这个模型的修正。下面来看这个模型在 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
9
10
tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
{
    const volScalarField nuEff_(nuEff());

    return
    (
      - fvm::laplacian(nuEff_, U)
      - fvc::div(nuEff_*dev(T(fvc::grad(U))))
    );
}

这个函数，写成公式是如下形式
$$
-\nabla \cdot (\nu_{eff}\nabla U)-\nabla \cdot [\nu_{eff}dev(\nabla ^TU)]
$$
OpenFOAM 中的 dev运算符定义为 $dev(T)=T-\frac{1}{3}tr(T)I$ ， tr 为张量的迹，即主对角元素之和所以。对于张量 $\nabla ^TU$ ，$tr(\nabla ^TU) =\nabla \cdot U$，因此上式可以写为【注一】
$$
-\nabla \cdot (\nu_{eff}\nabla U)-\nabla \cdot [\nu_{eff}(\nabla^TU-\frac{1}{3}(\nabla \cdot U))]
$$
上述 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();

这句代码负责湍流粘度的更新。这里又涉及到一个湍流类的成员函数： corretc。这个函数是湍流模型类的核心，因为湍流模型输运方程的求解，湍流粘度的更新都在这个函数里面。下面来看 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 ，这个部分的每一行写成公式，分别如下：
$$
\frac{\partial \tilde{\nu}}{\partial t} \\
\nabla \cdot (U \tilde{\nu}) \\
-\nabla \cdot [D_{\tilde{\nu}_{eff}} \nabla \tilde{\nu}] \\
-\frac{C_{b2}}{\sigma_{\nu_t}} \cdot |\nabla \tilde{\nu}|^2 \\
C_{b1}\cdot \tilde{S} \cdot \tilde{\nu} \\
-\frac{C_{w1}\cdot f_w\cdot \tilde{\nu}}{d^2} \cdot \tilde{\nu}\\
$$

最后两项都包含了需要求解的量 $\tilde{\nu}$， 但是处理方式却不一样，最后一项，其实是将原本的 $\tilde{\nu} ^2$ 项进行了线性化，并用 fvm::Sp 操作符进行隐式处理。而倒数第二项，则是作的显式处理。差别在于，倒数第二项仅仅是用上一步得到的 $\tilde{\nu}$ 代入，然后这一项被加到 $Ax=b$ 中的 $b$ 部分。而隐式处理则会对系数矩阵 $A$ 有贡献。
输运方程中涉及到的变量和函数，全部都在 SpalartAllmaras 类中有定义，这里就不再单独分析了，仅将涉及到的公式列出如下
$$
D_{\tilde{\nu}_{eff}}=\frac{\nu+\tilde{\nu}}{\sigma_{\nu_t}}
$$
$$
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) \\
r=min\left[ \frac{\tilde{\nu}}{max(\tilde{S}, SMALL)\cdot \kappa^2 d^2},10 \right]
$$

$$
\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)
$$

$f_{v2}$ 和 $f_{v3}$ 有点特殊，因为涉及到了 Ashford 对原始模型的修正。根据作者的论文，修正的目的是防止 $\tilde{S}$ 出现负值。
如果不引入 Ashford 修正，则 $f_{v3}=1$；引入以后，则
$$
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}
$$
对于 $f_{v2}$，如果不引入 Ashford 修正，则
$$
f_{v2}=\frac{1}{(1+\frac{\chi}{c_{v2}})^3}
$$
若引入，则
$$
f_{v2}=1-\frac{\chi}{1+\chi f_{v1}}
$$
此外，
$$
f_{v1}=\frac{\chi ^3}{\chi ^3 + c_{v1}^3}, \quad \chi=\frac{\tilde{\nu}}{\nu}, \quad \nu_t = \tilde{\nu}
f_{v1}$$
$d$ 是与壁面的距离，有一个专门的类 wallDist 来计算它。
其他的没有提到的变量，则都是模型常数了。具体的值都在 SpalartAllmaras 类的构造函数中定义了，这里就不再赘述了。

1.2. 一些细节
  1. turbulence_ 这个变量
    correct 函数中，只有当 turbulence_ 为 true 时，下面求解输运方程的部分代码才会执行。这个变量又是从 RASModel 类继承过来的，其初始化的语句为 turbulence_(lookup("turbulence")) ，这意味着，这个变量的值是从 RASProperties 文件中读取的（详见文末的 RASProperties 示例）。
  2. 模型常数，除了使用默认值，还可以自己指定（注：除非很有把握，一般不要随便修改模型常数）。具体方法是，在 RASProperties 文件里，建立一个名为 SpalartAllmarasCoeffs 的 subDict，详见文末的示例。
  3. bound 函数
    这个函数的定义在头文件 bound.H 里，其功能是限制 $\tilde{\nu}$ 的最小值。不过这种限制其实没有多少物理意义，纯粹是数值技巧。其中用到了 fvc::average 函数（见我的另一篇博文）。
  4. read 成员函数
    这个成员函数，初看之下似乎跟模型常数的初始化有关，实则不然。模型常数的初始化在类的构造函数部分完成， coeffDict_ 这个变量在基类 RASModel 的构造函数中先用 type()+Coeffs (这里是 SpalartAllmarasCoeffs )这个 subDict 初始化了，然后，在 SpalartAllmaras 类的构造函数中，每一个模型常数都是先去 coeffDict_ 中查找是否用用户自定义的值，如果没有，则用模型的默认值，并且将这个默认值添加到 coeffDict_ 这个字典中。
  5. 其他成员函数
    SpalartAllmaras 显然不需要用到湍动能 k 以及湍动能耗散 epsilon ，但是在 SpalartAllmaras.C 中确有这两个函数的定义，而且返回值是0。这是因为，在基类 turbulenceModel 中， k()， epsilon() 函数被声明为纯虚函数，如果不在 SpalartAllmaras 中重新定义这两个函数，则 SpalartAllmaras 也将是包括纯虚函数的抽象类了。在 C++ 中，抽象类是无法建立对象的。
    devReff， R ， divDevRhoReff 这几个成员函数，是湍流模型提供给其他模型调用的函数，比如，在计算表面应力的 forces 类中就需要调用湍流模型的 devReff 函数。
附： RASProperties 文件示例
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
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      RASProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

RASModel        SpalartAllmaras ;
turbulence      on; //若为 off，则湍流模型的输运方程将不会求解！
printCoeffs     on; // 输出模型常数

SpalartAllmarasCoeffs // 设置模型常数(仅为示例，具体数值不要较真)
{
    Cw2       0.20;
    Cw3       0.10;
}

// ************************************************************************* //

【注一】：N-S 方程中，这一项应该是
$$
-\nabla \cdot (\nu_{eff}\nabla U)-\nabla \cdot [\nu_{eff}(\nabla^TU-\frac{2}{3}(\nabla \cdot U))]
$$
而且，在可压缩的湍流模型中， divDevReff 函数用的是 dev2 函数，所以其定义的形式就跟 N-S 方程一致了。但不知道为什么不可压缩的湍流模型为什么要用 dev 函数。虽然对于不可压缩的情形，由于 $\nabla \cdot U=0$，用 dev 或者 dev2 应该关系不大，但这种做法仍然很费解。在OpenFOAM-3.0版本中，可压和不可压用的都是 dev2 函数了。

参考：

[1] P. Spalart and S. Allmaras. “A one-equation turbulence model for aerodynamic flows”. Technical Report AIAA-92-0439. American Institute of Aeronautics and Astronautics. 1992.

[2] “An Unstructured Grid Generation and Adaptive Solution Technique for High Reynolds Number Compressible Flows”, G.A. Ashford, Ph.D. thesis, University of Michigan, 1996.