Code of the Week #14: RASModel kEpsilon

rheoworks · 2019 年1 月 2 日 15:29

#201901

新年好！

新年第一弹，我们开讲湍流模型 k–\varepsilon 。在 Code of the Week #12: turbulence on 中，我们提到了一般湍流模型的设置。

今天，我们来看看湍流模型的由来。

首先，大家对 N-S 方程都已经不再陌生了。那湍流是什么呢？根据 wikipedia: 湍流，湍流是一种流动状态， “当流速很小时，流体分层流动，互不混合，称为层流，或称为片流；逐渐增加流速，流体的流线开始出现波浪状的摆动，摆动的频率及振幅随流速的增加而增加，此种流况称为过渡流；当流速增加到很大时，流线不再清楚可辨，流场中有许多小漩涡，称为湍流，又称为乱流（日本及港澳台用字）、扰流或紊流”。

简单地说，湍流就是当流体的惯性力影响大于黏滞力时，流动由规则分层明显的层流变为不规则运动的一种流动的状态。N-S 方程目前来说，只要网格足够细，是可以描述湍流的。然而，其计算量之大，是一般人或公司不能承受的，而也不是每一个问题都需要获得如此细致的解析。要想获得有效的／有用的解，或是较宏观的解，就需要将湍流模型化。主要有两种方法：雷诺平均法、LES 法。

我们先讲雷诺平均法。其基本思想就是把速度分为：平均速度和脉动速度之和，即

u=\overline{u} + u'

平均速度定义为某一时段的平均

\overline{u} = \frac{1}{\Delta t}\int_t^{t+\Delta t} u(t) d t

对于时间间隔 \Delta t 相对于湍流随机脉动的周期来说足够大，而相对于流场的各个时均量的缓慢变化来说又足够小。

对于连续性方程 (以不可压缩为例)

\nabla\cdot U=0

式中， U=(u,\, v,\, w) ，对其做时均处理，有

\overline{\nabla\cdot (\overline{U} + U')}=\nabla\cdot \overline U + \nabla\cdot \overline{U'}=0

对于时均值来讲，有

\overline{U'}=0,\, \overline{\overline{U}}=0,\, \overline{\overline{U}+U'}=\overline{U},\,\\ \frac{\overline{\partial U}}{\partial t } = \frac{\partial \overline{U}}{\partial t },\,\, \frac{\overline{\partial U_i}}{\partial j}=\frac{\partial \overline{U_i}}{\partial j},\,\, \frac{\overline{\partial U'_i}}{\partial j}=\frac{\partial \overline{U'_i}}{\partial j}=0

因而，有

\nabla\cdot \overline U =0\\ \nabla\cdot \overline{U'}=0

因此，可以说，湍流的时均值也满足连续性方程。

对于动量方程来说，

\begin{equation} \frac{\partial U}{\partial t} + \nabla\cdot(UU) - \nabla\cdot(\nu\nabla U) = -\nabla \frac{p}{\rho} \end{equation} \tag{}

将其写成脉动形式，并做平均处理，有

\begin{equation} \frac{\overline{\partial (\overline U + U')}}{\partial t} + \nabla\cdot((\overline{\overline U + U'})(\overline{\overline U + U'})) - \nabla\cdot(\nu\nabla (\overline{\overline U + U'})) = -\nabla \frac{\overline{\overline p + p'}}{\rho} \end{equation} \tag{}

先取 x 方向的来说明，

\begin{equation} \frac{\overline{\partial (\overline u + u')}}{\partial t} + \nabla\cdot((\overline{\overline u + u'})(\overline{\overline U + U'})) - \nu(\nabla^2 (\overline{\overline U + U'})) = \frac{\partial \overline{\overline p + p'}}{\rho\partial x} \end{equation} \tag{}

化简后，有

\begin{equation} \frac{\partial \overline u }{\partial t} + \nabla\cdot(\overline u \overline U) + \nabla\cdot{\overline {u' U'}} - \nu(\nabla^2 \overline{U}) = \frac{\partial \overline p }{\rho \partial x} \end{equation} \tag{}

将脉动量和耗散项移动到等式右边，有

\begin{equation} \frac{\partial \overline u }{\partial t} + \nabla\cdot(\overline u \overline U) = \frac{\partial \overline p }{\rho \partial x}+ \nabla\cdot(v \nabla\overline U - \overline {u' U'}) \end{equation} \tag{}

写成全张量形式，有

\begin{equation} \frac{\partial \overline U }{\partial t} + \nabla\cdot(\overline U \overline U) = \frac{1}{\rho }\nabla \overline p + \nabla\cdot(v \nabla\overline U - \overline {U' U'}) \end{equation} \tag{}

若有其他变量 \phi ，其湍流方程为

\begin{equation} \frac{\partial \overline \phi }{\partial t} + \nabla\cdot(\overline \phi \overline U) = \frac{1}{\rho }\nabla \overline p+ \nabla\cdot(\varGamma \nabla\overline \phi - \overline {\phi' U'}) \end{equation} \tag{}

由上述方程可以看出，一次项保持不变，而二次项产生了包含脉动值乘积的附加项 \rho \overline{U' U'} ，该项代表了由湍流（脉动）引起的能量转移 (应力、热流密度等)，其中， \rho \overline{U' U'} 称为湍流应力或 Reynolds 应力， \overline {\phi' U'} 称为 turbulent scalar flux。

因此，整个方程组有 5 个平均量 (\overline U,\,\overline p,\,\overline\phi)，6 个脉动值乘积的附加项 (\overline{U' U'})，要使这些方程组封闭可求解，必须补充这 6 个附加量的关系，并且不能引入新的变量。而湍流模型则是把这些脉动值的附件项与时均量联系起来的一些特定关系式，依靠理论与经验的结合，引进一系列模型假设，建立一组描写湍流平均量的封闭方程组的理论计算方法。Reynolds 应力的主要贡献来自大尺度脉动，而大尺度脉动的性质及结果和流动的边界条件密切相关，因此，雷诺应力的封闭模型不是普适的，也就是说，没有适用于一切复杂流动的统一湍流模型。

为了使方程组封闭从而可求解，求得 Reynolds 应力，根据 Boussinesq 假设，即 Reynolds 应力正比于时均速度的应变，系数为涡粘系数 (Eddy viscosity, \mu_t)

\rho \overline{u_i' u_j'} = \mu_t\left( \frac{\partial \overline u_i}{\partial x_j} + \frac{\partial \overline u_j}{\partial x_i}\right) -\frac{2}{3}\delta_{ij}\left(\rho k +\mu_t \frac{\partial \overline u_m}{\partial \overline x_m} \right)

式中， k=\dfrac{1}{2}\sum\overline {u_i'u_i'} 为湍动能，涡粘系数不是流体属性，而是湍流的特征量，随流体流动的位置和状态而改变。尽管这一假说在细节上并不成立，但它很容易应用，并能获得很多流动合理的解。

使用最简单的模型，湍流可由两个参数表征：湍动能 k 或速度 q=\sqrt{2k} ，和湍流长度 L ，

\mu_t = C_\mu \rho q L

式中， C_\mu 是无量纲数。

基于这一假说的模型称为涡粘模型 (Eddy Viscosity Model, EVM)，最简单的涡粘模型就是混合长度模型 (mixed-length model)， k 可从平均速度场近似得来 q=L\dfrac{\partial u}{\partial y} ， L 是空间的一个函数。准确的 L 分布对于简单流动来说是可行的，但是不能用于复杂的三维流动。混合长度模型因此仅能用于较简单的流动。混合长度模型与称为零方程模型 (zero-equation models)。

模化湍流很困难，但是，可不可以从偏微分方程中去计算呢？比如，用湍动能来表征湍流的速度，

\frac{\partial \rho k}{\partial t} + \frac{\partial (\rho \overline u_j k)}{\partial x_j} = \frac{\partial}{\partial x_j}\left(\mu\frac{\partial k}{\partial x_j} \right) -\frac{\partial}{\partial x_j}\left(\frac{\rho}{2}\overline{u'_j u'_i u'_i} +\overline{p' u'_j} \right) -\rho\overline{u'_j u'_i }\frac{\partial \overline u_i}{\partial x_j} - \mu \overline{\frac{\partial u'_i}{\partial x_k} \frac{\partial u'_i}{\partial x_k}}

等式右边最后一项代表的是密度与耗散率 \varepsilon 的乘积，耗散率指的是湍流转化为内能的速率。

等式右边第二项代表湍动能的扩散，实质上是脉动本身输运脉动，可有如下假设

-\left(\frac{\rho}{2}\overline{u'_j u'_i u'_i} +\overline{p' u'_j} \right)=\frac{\mu_t}{\sigma_k}\frac{\partial k}{\partial x_j}

式中， \sigma_k 是湍流 Prandtl 数，近似为 1。在一些更为复杂的模型中， \mu_t 成为了张量。

右边第三项代表湍动能与平均流的乘积的变化率，动能从平均流到湍流的转化，如果用涡粘假设来估计 Reynolds 应力，有

P_r = -\rho\overline{u'_j u'_i }\frac{\partial \overline u_i}{\partial x_j} =\mu_t\left( \frac{\partial \overline u_i}{\partial x_j} + \frac{\partial \overline u_j}{\partial x_i}\right) \frac{\partial \overline u_i}{\partial x_j}

至此，湍动能的方程就完备了。

上面已经提到对于湍流的表征，除了速度，还需要长度。对于长度用什么参数来表征，可以有很多选择。最广泛的是耗散率 \varepsilon

\varepsilon=\frac{k^{3/2}}{L}

该方程是基于惯性能转化的估计。对于 \varepsilon

\frac{\partial (\rho \varepsilon)}{\partial t} +\frac{\partial (\rho u_j \varepsilon)}{\partial x_j} = C_{\varepsilon1} P_k \frac{\epsilon}{k} -\rho C_{\varepsilon2}\frac{\varepsilon^2}{k} +\frac{\partial}{\partial x_j}\left(\frac{\mu_t}{\sigma_\varepsilon} \frac{\partial \varepsilon}{\partial x_j} \right)

式 () 和 () 就称为 k–\varepsilon 模型。有五个参数，广泛使用的值为

C_\mu = 0.09;\, C_{\varepsilon1}=1.44;\, C_{\varepsilon2}=1.92;\, \sigma_k=1.0;\,\sigma_\varepsilon=1.3

使用这一湍流模型时，原动量方程中的粘度 \mu 被替换成 \mu_\text{eff} = \mu + \mu_t ，再外加两个 k 和 \varepsilon 的方程求解。

OpenFOAM 中 k–\varepsilon 实现

由于增加了方程及系数，就需要定义变量，在 OpenFOAM 中，湍流模型位于 src/turbulenceModels，先来看一看 turulenceModel 声明，文件是 src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H

/*---------------------------------------------------------------------------*\
Description
    Abstract base class for incompressible turbulence models
    (RAS, LES and laminar).
...
\*---------------------------------------------------------------------------*/
namespace Foam
{

// Forward declarations
class fvMesh;

namespace incompressible
{

/*---------------------------------------------------------------------------*\
                           Class turbulenceModel Declaration
\*---------------------------------------------------------------------------*/

class turbulenceModel 
{

protected:                 // 受保护数据，取的是地址

    // Protected data

        const Time& runTime_;
        const fvMesh& mesh_;

        const volVectorField& U_;
        const surfaceScalarField& phi_;

        transportModel& transportModel_;


private:                                    // 私有成员

    // Private Member Functions

        //- Disallow default bitwise copy construct
        turbulenceModel(const turbulenceModel&); // 构造函数

        //- Disallow default bitwise assignment
        void operator=(const turbulenceModel&);


public:                      // 公开成员

    //- Runtime type information
    TypeName("turbulenceModel");        // 运行时类型名字


    // Declare run-time New selection table // 声明运行时新建列表

#ifndef SWIG
        declareRunTimeNewSelectionTable
        (
            autoPtr,
            turbulenceModel,
            turbulenceModel,
            (
                const volVectorField& U,
                const surfaceScalarField& phi,
                transportModel& lamTransportModel
            ),
            (U, phi, lamTransportModel)
        );
#endif


    // Constructors

        //- Construct from components
        turbulenceModel
        (
            const volVectorField& U,
            const surfaceScalarField& phi,
            transportModel& lamTransportModel
        );


    // Selectors 选择器

        //- Return a reference to the selected turbulence model
        static autoPtr<incompressible::turbulenceModel> New
        (
            const volVectorField& U,
            const surfaceScalarField& phi,
            transportModel& lamTransportModel
        );


    // Destructor

        virtual ~turbulenceModel()
        {}


    // Member Functions

        //- Access function to velocity field
        inline const volVectorField& U() const
        {
            return U_;
        }

        //- Access function to flux field
        inline const surfaceScalarField& phi() const
        {
            return phi_;
        }

        //- Access function to incompressible transport model
        inline transportModel& transport() const
        {
            return transportModel_;
        }

        //- Return the laminar viscosity
        const volScalarField& nu() const  // 返回粘度
        {
            return transportModel_.nu();
        }

        //- Return the turbulence viscosity // 返回湍流粘度
        virtual tmp<volScalarField> nut() const = 0;

        //- Return the effective viscosity // 返回有效粘度
        virtual tmp<volScalarField> nuEff() const = 0;

        //- Return the turbulence kinetic energy // 湍动能，在每个湍流模型中重新定义
        virtual tmp<volScalarField> k() const = 0;

        //- Return the turbulence kinetic energy dissipation rate // 耗散率，在每个湍流模型中重新定义
        virtual tmp<volScalarField> epsilon() const = 0;

        //- Return the Reynolds stress tensor // 雷诺应力张量，在每个湍流模型中重新定义
        virtual tmp<volSymmTensorField> R() const = 0;

        //- Return the effective stress tensor including the laminar stress // 应力张量，含层流应力 (平均应力)
        virtual tmp<volSymmTensorField> devReff() const = 0;

        //- Return the source term for the momentum equation // 返回源项，这一项加入到动量方程中
        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const = 0;

        //- Solve the turbulence equations and correct the turbulence viscosity // 求解湍流方程，修正湍流粘度
        virtual void correct() = 0;

        //- Read turbulenceProperties dictionary // 读取 turbulenceProperties 配置
        virtual bool read() = 0;
};


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

} // End namespace incompressible
} // End namespace Foam

对于每一个单独的湍流模型，以上成员函数均须重新定义。比如 laminar 模型，其定义如下 (laminar.C)

tmp<volScalarField> laminar::nut() const
{
    return tmp<volScalarField>  // 临时体心标量场
    (
        new volScalarField
        (
            IOobject
            (
                "nut",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh_,
            dimensionedScalar("nut", nu().dimensions(), 0.0) // 返回 0
        )
    );
}


tmp<volScalarField> laminar::nuEff() const
{
    return tmp<volScalarField>(new volScalarField("nuEff", nu()));  // 返回 nu
}


tmp<volScalarField> laminar::k() const
{
    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                "k",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh_,
            dimensionedScalar("k", sqr(U_.dimensions()), 0.0)  // 返回 0
        )
    );
}


tmp<volScalarField> laminar::epsilon() const
{
    return tmp<volScalarField> 
    (
        new volScalarField
        (
            IOobject
            (
                "epsilon",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh_,
            dimensionedScalar
            (
                "epsilon", sqr(U_.dimensions())/dimTime, 0.0  // 返回 0
            )
        )
    );
}

若是 k – \varepsilon 模型，其增加的方程在 correct() 函数里面 (文件是 kEpsilon.C)

void kEpsilon::correct()
{
    // Bound in case of topological change
    // HJ, 22/Aug/2007
    if (mesh_.changing())
    {
        bound(k_, k0_);
        bound(epsilon_, epsilon0_);
    }

    RASModel::correct();

    if (!turbulence_)
    {
        return;
    }

    volScalarField G("RASModel::G", nut_*2*magSqr(symm(fvc::grad(U_))));

    // Update epsilon and G at the wall
    epsilon_.boundaryField().updateCoeffs();

    // Dissipation equation
    tmp<fvScalarMatrix> epsEqn         // epsilon 方程
    (
        fvm::ddt(epsilon_)
      + fvm::div(phi_, epsilon_)
      + fvm::SuSp(-fvc::div(phi_), epsilon_)
      - fvm::laplacian(DepsilonEff(), epsilon_)
     ==
        C1_*G*epsilon_/k_
      - fvm::Sp(C2_*epsilon_/k_, epsilon_)
    );

    epsEqn().relax();

    // No longer needed: matrix completes at the point of solution
    // HJ, 17/Apr/2012
//     epsEqn().completeAssembly();

    solve(epsEqn);
    bound(epsilon_, epsilon0_);


    // Turbulent kinetic energy equation
    tmp<fvScalarMatrix> kEqn             // k 方程
    (
        fvm::ddt(k_)
      + fvm::div(phi_, k_)
      + fvm::SuSp(-fvc::div(phi_), k_)
      - fvm::laplacian(DkEff(), k_)
     ==
        G
      - fvm::Sp(epsilon_/k_, k_)
    );

    kEqn().relax();
    solve(kEqn);
    bound(k_, k0_);


    // Re-calculate viscosity    重新计算涡粘
    nut_ = Cmu_*sqr(k_)/epsilon_;
    nut_.correctBoundaryConditions();
}

除了 k–\varepsilon 模型，还有 k–\omega 等模型。

雷诺应力模型

涡粘模型是基于湍流各项同性假设的，对于各向异性的复杂流动 (如大曲率流动、强旋流、冲击流等) 通常得不到好的结果。这时，雷诺应力模型 (Reynolds Stress Model, RSM) 应运而生，在更高阶的应力上封闭 RANS 方程。通常做法是求解六个未知应力分量的方程及其对应的耗散率方程，因而可以较好地模拟具有强烈非线性的复杂流动。然而它需要求解许多额外的方程，计算量相比涡粘模型大得多，因此在工程中应用不多。常见的雷诺应力模型有 LRR (Launder-Reece-Rodi) 和 SSG (Speziale-Sarkar-Gatski) 模型。

参考文献

Ferziger J H, Peric M. Computational Methods for Fluid Dynamics, 3rd, rev. ed., Berlin: Springer, 2002