OpenFOAM类库介绍（六）对流项高分辨率格式理论

高分辨率格式

简介

一般地，有限体积法将物理量存储在体心。但是如果需要完成对流项的离散，那么必须将物理量插值到面心。如果使用线性插值（一般也称作中心差分）会导致在较大贝克莱数下发生非物理震荡。所以人们提出了一阶迎风格式，直接用上游的体心物理量作为面心物理量。一阶迎风格式不会发生非物理振荡，但是精度很低。为了提高精度，人们又提出了大量高阶格式（HO），遗憾的是，这些高阶格式并不能保证物理单调性。

在高阶格式（HO）的基础上，人们提出了高分辨率格式（HR）。其中以TVD格式为主要代表（也是OpenFOAM所采用的）。其原理是设置一个通量限制器（limiter），当物理量变化剧烈时退回一阶迎风格式，当物理量变化平缓时返回高阶格式，有效避免振荡。

理论推导

如下图所示，考虑一个单元面

在OpenFOAM中任何一个内部面都有一个owner单元（标记为O）和一个neighbour单元（标记为N），面的法向量由O指向N

面心的物理量可以由体心线性表示
$${\varphi }_{\mathrm{f}}=\lambda {\varphi }_O+\left(1-\lambda \right){\varphi }_N\text{\hspace{1em}(1)}$$式中$\lambda$称为权重，在OpenFOAM中也是用$\lambda$这个希腊字母的。

特殊地，对于均匀网格，如果采用中心差分（CD），那么$\lambda ={\lambda }_{\mathrm{CD}}=0.5$。

高分辨率格式使用通量限制器$\psi$控制权重，其计算方式为:
（1）如果流向为从O到N，那么
$$\lambda ={\lambda }_{\mathrm{CD}}\psi +\left(1-\psi \right),O\to N$$
（2）如果流向为从N到O，那么
$$\lambda ={\lambda }_{\mathrm{CD}}\psi ,N\to O$$式中${\lambda }_{\mathrm{CD}}$是使用中心差分（CD）得到的权重。需要注意，通量限制器$\psi$是一个场，任何一个单元面都应该有自己的$\psi$值。

接下来将权重表达式带入式(1)中，化简得到
（1）如果流向为从O到N，那么
$${\varphi }_{\mathrm{f}}={\varphi }_O+\frac{\psi }{2}\left({\varphi }_N-{\varphi }_O\right)$$
（1）如果流向为从N到O，那么
$${\varphi }_{\mathrm{f}}={\varphi }_N+\frac{\psi }{2}\left({\varphi }_O-{\varphi }_N\right)$$
其物理意义是非常明确的，第一项表示一阶迎风，即上游的体心物理量，由于一阶迎风有较大的数值粘性，第二项做出修正，该修正可以使数值粘性减小，最终达到二阶精度。

通量限制器$\psi$的要求是：当物理量变化剧烈时退回一阶迎风格式，当物理量变化平缓时返回高阶格式，那么需要引入一个无量纲数$r$来检测物理量变化是否剧烈，定量描述剧烈程度。不同的高阶格式即就是设计不同的$\psi (r)$函数。这里给出$r$的定义
$$r=\frac{{\varphi }_{cen}-{\varphi }_{up}}{{\varphi }_{down-}{\varphi }_{cen}}\text{\hspace{1em}(2)}$$
下角标$cen$表示所考察面的上游单元，$up$是$cen$的上游单元，$down$是$cen$的下游单元

但是遗憾的是，对于非结构网格，我们可以很容易知道$cen$和$down$单元，但是很难知道$up$单元是谁，所以需要用线性近似的方式得到 ${\varphi }_{up}$值：
$${\varphi }_{up}={\varphi }_{down}-2d_{cen-down}\cdot {\left(\nabla \varphi \right)}_{cen}$$
式中$d$是一个距离向量，从$cen$的体心指向$down$体心，将上式带入式(2)可得
$$r=\frac{2d_{cen-down}\cdot {\left(\nabla \varphi \right)}_{cen}}{{\varphi }_{down-}{\varphi }_{cen}}-1$$
考虑到OpenFoam的网格格式（owner单元和neigbour)：
（1）如果流向为从O到N，那么将$cen$替换为O，将$down$替换为N，可得
$$r=\frac{2d_{ON}\cdot {\left(\nabla \varphi \right)}_O}{{\varphi }_{N-}{\varphi }_O}-1$$
（2）如果流向为从N到O，那么将$cen$替换为N，将$down$替换为O，可得
$$r=\frac{2d_{NO}\cdot {\left(\nabla \varphi \right)}_N}{{\varphi }_{O-}{\varphi }_N}-1=\frac{2d_{ON}\cdot {\left(\nabla \varphi \right)}_N}{{\varphi }_{N-}{\varphi }_O}-1$$

代码简介

OpenFoam中对r的定义位于NVDTVD.H文件

scalar r
(
    const scalar faceFlux,
    const scalar phiP,//这里P就是owner
    const scalar phiN,
    const vector& gradcP,//这里P就是owner
    const vector& gradcN,
    const vector& d
) const
{
    scalar gradf = phiN - phiP;

    scalar gradcf;

    if (faceFlux > 0)
    {
        gradcf = d & gradcP;
    }
    else
    {
        gradcf = d & gradcN;
    }

    if (mag(gradcf) >= 1000*mag(gradf))
    {
        return 2*1000*sign(gradcf)*sign(gradf) - 1;
    }
    else
    {
        return 2*(gradcf/gradf) - 1;
    }
}

一阶迎风格式对$\psi$的定义位于upwind.H文件，是均匀的0场

 //- Return the interpolation limiter
 virtual tmp<surfaceScalarField> limiter
 (
     const GeometricField<Type, fvPatchField, volMesh>&
 ) const
 {
     return surfaceScalarField::New
     (
         "upwindLimiter",
         this->mesh(),
         dimensionedScalar(dimless, 0)
     );
 }

SuperBee 格式对$\psi$的定义位于SuperBee.H

scalar limiter
(
    const scalar cdWeight,
    const scalar faceFlux,
    const typename LimiterFunc::phiType& phiP,
    const typename LimiterFunc::phiType& phiN,
    const typename LimiterFunc::gradPhiType& gradcP,
    const typename LimiterFunc::gradPhiType& gradcN,
    const vector& d
) const
{
    scalar r = LimiterFunc::r
    (
        faceFlux, phiP, phiN, gradcP, gradcN, d
    );

    return max(max(min(2*r, 1), min(r, 2)), 0);
}

Minmod 格式对$\psi$的定义位于Minmod.H

scalar limiter
(
    const scalar cdWeight,
    const scalar faceFlux,
    const typename LimiterFunc::phiType& phiP,
    const typename LimiterFunc::phiType& phiN,
    const typename LimiterFunc::gradPhiType& gradcP,
    const typename LimiterFunc::gradPhiType& gradcN,
    const vector& d
) const
{
    scalar r = LimiterFunc::r
    (
        faceFlux, phiP, phiN, gradcP, gradcN, d
    );

    return max(min(r, 1), 0);
}

vanLeer格式对$\psi$的定义位于vanLeer.H

scalar limiter
(
    const scalar,
    const scalar faceFlux,
    const typename LimiterFunc::phiType& phiP,
    const typename LimiterFunc::phiType& phiN,
    const typename LimiterFunc::gradPhiType& gradcP,
    const typename LimiterFunc::gradPhiType& gradcN,
    const vector& d
) const
{
    scalar r = LimiterFunc::r
    (
        faceFlux, phiP, phiN, gradcP, gradcN, d
    );

    return (r + mag(r))/(1 + mag(r));
}

类库的具体结构和代码注释将在之后的文章中给出