非结构网格如何计算库朗数

库朗数是用来描述流体“数值流动”状态的一个无量纲参数，它关系到数值求解的稳定性。要了解库朗数的定义需要先知道什么是 CFL 条件。

CFL 条件

CFL 条件以 Courant、Friedrichs 和 Lewy 三人姓的首字母命名，最初是一个纯数学概念，用来证明某些偏微分方程解的存在性。CFL 条件指出：数值依赖域必须包含物理依赖域。

我们举一个简单的例子来说明什么是依赖域（domain of dependence）。

考虑一维对流方程

$$ \begin{equation} \frac{\partial \phi}{\partial t} + c \frac{\partial \phi}{\partial x} = 0 \end{equation} $$

其中 $\phi$ 是需要传输的物理量，$c$ 为传播速度，$t$ 和 $x$ 代表时间和空间。

采用有限差分法对上述方程进行离散，时间项采用显式方式处理，即对流项采用已知值。时间项和空间项均采用前向差分格式，得到

$$ \begin{equation} \frac{\phi_j^{n+1}-\phi_j^{n}}{\Delta t} + c \frac{\phi_{j+1}^{n}-\phi_j^{n}}{\Delta x} = 0 \end{equation} $$

其中上标 $n$ 和 $n+1$ 表示离散后的时间节点，下标 $j$ 和 $j+1$ 表示离散后的空间节点。

整理上式后得到

$$ \begin{equation} \phi_j^{n+1} = \left(1-c \frac{\Delta t}{\Delta x}\right)\phi_j^{n}- c \frac{\Delta t}{\Delta x}\phi_{j+1}^{n} \end{equation} $$

上式表示 $n+1$ 时刻 $j$ 节点的值通过 $n$ 时刻 $j$ 节点和 $j+1$节点的值求得。以此类推，最终得到一个如下图所示的黑点标注的区域，这个区域即为数值依赖域。

数值依赖域

如果 $c > 0$，那么从物理上看， $\phi$ 朝 $x$ 正方向输运。$\phi_j^{n+1}$ 的值应该受左侧（上游）节点的影响，这些上游节点组成的域被称为物理依赖域。实际上这种情况即为我们通常说的顺风（downwind）格式。这种情况下的数值依赖域不包含物理依赖域，因此不满足 CFL 条件，是不稳定的。

按照以上定义，当方程的时间项采用隐式方式处理时，整个离散域上的节点都会对 $\phi_j^{n+1}$ 产生影响。这种情况下的数值依赖域为整个离散域，包含所有物理依赖域，因此是无条件稳定的。

库朗数

在求解实际问题时，对流项通常采用迎风（upwind）格式离散。迎风格式的数值依赖域会随着速度的方向变化而改变，呈左右对称分布。定义以下参数

$$ \begin{equation} Co = |c|\frac{\Delta t}{\Delta x} \label{eq:courant_number} \end{equation} $$

这个参数被称为库朗数或 CFL 条件数。为了保证物理依赖域在数值依赖域以内，库朗数需满足 $Co < 1$。

以上是一维情况下的库朗数定义，二维情况下的库朗数定义如下

$$ \begin{equation} Co = |u| \frac{\Delta t}{\Delta x} + |v| \frac{\Delta t}{\Delta y} \end{equation} $$

上面讨论的都是有限差分法中的库朗数定义。对于基于有限体积法的三维非结构网格，网格单元的库朗数的定义如下

$$ \begin{equation} Co = \frac{\Delta t}{\Delta V} \frac{1}{2}\sum_f{|\phi_f|} \label{eq:courant_num} \end{equation} $$

其中，$\frac{1}{2}\sum_f{|\phi_f|}$ 代表 $\Delta t$ 时间内流入或流出网格单元的流量。

OpenFOAM 实现

OpenFOAM 3.0.x

OpenFOAM 3.0.x 计算库朗数的代码有两份：一份是 utility，编译后生成可执行文件，可以直接用 Co 命令调用；另一份是 function object，编译后生成动态链接库，通过 function objects 机制调用。二者的代码分别位于：

applications/utilities/postProcessing/velocityField/Co/
src/postProcessing/functionObjects/utilities/CourantNo/

utility 中的代码如下：

 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


        if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
        {
            Info<< "    Calculating compressible Co" << endl;

            Info<< "    Reading rho" << endl;
            volScalarField rho
            (
                IOobject
                (
                    "rho",
                    runTime.timeName(),
                    mesh,
                    IOobject::MUST_READ
                ),
                mesh
            );

            Co.dimensionedInternalField() =
                (0.5*runTime.deltaT())
               *fvc::surfaceSum(mag(phi))().dimensionedInternalField()
               /(rho*mesh.V());
            Co.correctBoundaryConditions();
        }
        else if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
        {
            Info<< "    Calculating incompressible Co" << endl;

            Co.dimensionedInternalField() =
                (0.5*runTime.deltaT())
               *fvc::surfaceSum(mag(phi))().dimensionedInternalField()
               /mesh.V();
            Co.correctBoundaryConditions();
        }

首先判断 phi 的量纲：若为 $kg \cdot s^{-1}$，则为质量流量，需要将计算得到的结果除以密度。若为 $m^3 \cdot s^{-1}$，则为体积流量，可依据公式 $\eqref{eq:courant_num}$ 直接计算。

function object 中的代码如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24


void Foam::CourantNo::execute()
{
    if (active_)
    {
        const fvMesh& mesh = refCast<const fvMesh>(obr_);

        const surfaceScalarField& phi =
            mesh.lookupObject<surfaceScalarField>(phiName_);

        volScalarField& Co =
            const_cast<volScalarField&>
            (
                mesh.lookupObject<volScalarField>(type())
            );

        Co.dimensionedInternalField() = byRho
        (
            (0.5*mesh.time().deltaT())
           *fvc::surfaceSum(mag(phi))().dimensionedInternalField()
           /mesh.V()
        );
        Co.correctBoundaryConditions();
    }
}

其中 byRho 函数中会判断库朗数的量纲，若量纲为密度的量纲dimDensity，则说明通量是质量通量，需要将得到的库朗数除以密度：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


Foam::tmp<Foam::volScalarField::DimensionedInternalField>
Foam::CourantNo::byRho
(
    const tmp<volScalarField::DimensionedInternalField>& Co
) const
{
    if (Co().dimensions() == dimDensity)
    {
        return Co/obr_.lookupObject<volScalarField>(rhoName_);
    }
    else
    {
        return Co;
    }
}

OpenFOAM 4.x 及后续版本

OpenFOAM 4.x 及后续版本中库朗数的计算代码与 OpenFOAM 3.0.x 中的 function object 差别不大，代码位于

src/functionObjects/field/CourantNo/

相应代码如下：

 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


bool Foam::functionObjects::CourantNo::calc()
{
    if (foundObject<surfaceScalarField>(fieldName_))
    {
        const surfaceScalarField& phi =
            lookupObject<surfaceScalarField>(fieldName_);

        tmp<volScalarField> tCo
        (
            new volScalarField
            (
                IOobject
                (
                    resultName_,
                    mesh_.time().timeName(),
                    mesh_
                ),
                mesh_,
                dimensionedScalar("0", dimless, 0.0),
                zeroGradientFvPatchScalarField::typeName
            )
        );

        tCo->ref() =
            byRho
            (
                (0.5*mesh_.time().deltaT())
               *fvc::surfaceSum(mag(phi))()()
               /mesh_.V()
            );

        tCo->correctBoundaryConditions();

        return store(resultName_, tCo);
    }
    else
    {
        return false;
    }
}

非结构网格如何计算库朗数

Contents

CFL 条件

库朗数

OpenFOAM 实现

OpenFOAM 3.0.x

OpenFOAM 4.x 及后续版本

参考资料