库朗数是用来描述流体“数值流动”状态的一个无量纲参数,它关系到数值求解的稳定性。要了解库朗数的定义需要先知道什么是 CFL 条件。
CFL 条件
CFL 条件以 Courant、Friedrichs 和 Lewy 三人姓的首字母命名,最初是一个纯数学概念,用来证明某些偏微分方程解的存在性。CFL 条件指出:数值依赖域必须包含物理依赖域。
我们举一个简单的例子来说明什么是依赖域(domain of dependence)。
考虑一维对流方程
其中 是需要传输的物理量, 为传播速度, 和 代表时间和空间。
采用有限差分法对上述方程进行离散,时间项采用显式方式处理,即对流项采用已知值。时间项和空间项均采用前向差分格式,得到
其中上标 和 表示离散后的时间节点,下标 和 表示离散后的空间节点。
整理上式后得到
上式表示 时刻 节点的值通过 时刻 节点和 节点的值求得。以此类推,最终得到一个如下图所示的黑点标注的区域,这个区域即为数值依赖域。
如果 ,那么从物理上看, 朝 正方向输运。 的值应该受左侧(上游)节点的影响,这些上游节点组成的域被称为物理依赖域。实际上这种情况即为我们通常说的顺风(downwind)格式。这种情况下的数值依赖域不包含物理依赖域,因此不满足 CFL 条件,是不稳定的。
按照以上定义,当方程的时间项采用隐式方式处理时,整个离散域上的节点都会对 产生影响。这种情况下的数值依赖域为整个离散域,包含所有物理依赖域,因此是无条件稳定的。
库朗数
在求解实际问题时,对流项通常采用迎风(upwind)格式离散。迎风格式的数值依赖域会随着速度的方向变化而改变,呈左右对称分布。定义以下参数
这个参数被称为库朗数或 CFL 条件数。为了保证物理依赖域在数值依赖域以内,库朗数需满足 。
以上是一维情况下的库朗数定义,二维情况下的库朗数定义如下
上面讨论的都是有限差分法中的库朗数定义。对于基于有限体积法的三维非结构网格,网格单元的库朗数的定义如下
其中, 代表 时间内流入或流出网格单元的流量。
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
的量纲:若为 ,则为质量流量,需要将计算得到的结果除以密度。若为 ,则为体积流量,可依据公式 直接计算。
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;
}
}
|
参考资料
- http://web.mit.edu/16.90/BackUp/www/pdfs/Chapter14.pdf
- https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition
- https://www.openfoam.com/documentation/cpp-guide/html/guide-fos-field-courant-no.html
- https://www.cfd-online.com/Forums/main/2989-cfl-condition.html
- https://www.cfd-online.com/Forums/main/9374-timestep-unstructured-grid.html
- https://www.cfd-online.com/Forums/openfoam-pre-processing/121982-computing-courant-number.html