库朗数是用来描述流体“数值流动”状态的一个无量纲参数,它关系到数值求解的稳定性。要了解库朗数的定义需要先知道什么是 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} \end{equation} \label{eq:courant_number} \]

这个参数被称为库朗数或 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|} \end{equation} \label{eq:courant_num} \]

其中,$\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
34
        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
25
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
16
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
41
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;
    }
}

参考资料

  1. http://web.mit.edu/16.90/BackUp/www/pdfs/Chapter14.pdf
  2. https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition
  3. https://www.openfoam.com/documentation/cpp-guide/html/guide-fos-field-courant-no.html
  4. https://www.cfd-online.com/Forums/main/2989-cfl-condition.html
  5. https://www.cfd-online.com/Forums/main/9374-timestep-unstructured-grid.html
  6. https://www.cfd-online.com/Forums/openfoam-pre-processing/121982-computing-courant-number.html