OpenFOAM 中的有限体积离散

介绍

这篇文章中你将会学到：

基于任意多面体网格的有限体积离散方法
时间项的处理方式
对流项的处理方式
源项的处理方式
边界条件的处理方式
OpenFOAM 中对有限体积离散的代码实现

基本概念

OpenFOAM 采用基于同位网格（collocated grid）的有限体积法进行离散，所有的变量都存储在控制体网格内部的同一个点上。

(图片来自 Hrvoje Jasak¹)

考虑有如图所示的任意多面体网格单元 $V_P$ ，$P$ 为其体心，定义如下

$$ \begin{equation} \int_{V_P} (\mathbf x - \mathbf x _P) dV = 0 \end{equation} $$

基于以上定义，假设物理量在控制体内线性变化，位于单元内 $\mathbf x$ 处的任意物理量可表示为

$$ \begin{equation} \phi (\mathbf x) = \phi_P + (\mathbf x - \mathbf x_P) \cdot (\nabla \phi)_P \end{equation} $$

将上式对单元进行体积分，可得

$$ \begin{equation} \begin{aligned} \int_{V_P} \phi (\mathbf x) dV &= \int_{V_ % workaround for hugo P} [\phi_ % workaround for hugo P + (\mathbf x - \mathbf x_ % workaround for hugo P) \cdot (\nabla \phi)_ % workaround for hugo P ] dV \newline &= \phi_P V_ % workaround for hugo P + (\nabla \phi)_ % workaround for hugo P \cdot \int_{V_ % workaround for hugo P} (\mathbf x - \mathbf x_ % workaround for hugo P) dV \newline & = \phi_P V_P \end{aligned} \end{equation} $$

对流扩散方程

考虑一般形式的对流扩散方程，其形式可以表示如下

$$ \begin{equation} \underbrace{\frac{\partial \phi}{\partial t}}_ % workaround for hugo {\text{时间项}} + \underbrace{\nabla \cdot (\phi \mathbf u)}_ % workaround for hugo {\text{对流项}} - \underbrace{\nabla \cdot (\gamma \nabla \phi)}_ % workaround for hugo {\text{扩散项或粘性项}} = \underbrace{S_\phi(\phi)}_ % workaround for hugo {\text{源项}} \end{equation} $$

将其对单元做体积分，可得

$$ \begin{equation} \int_{V_P} \frac{\partial \phi}{\partial t} dV + \int_{V_P} \nabla\cdot (\phi \mathbf u) dV - \int_{V_P} \nabla \cdot (\gamma \nabla \phi) dV = \int_{V_P} S_\phi(\phi) dV \end{equation} $$

可通过高斯定理将对流项和扩散项的体积分转为面积分，上式变为

$$ \begin{equation} \int_{V_P} \frac{\partial \phi}{\partial t} dV + \sum_{f} [\mathbf S_f \cdot (\phi \mathbf u)_ % workaround for hugo f] - \sum_{f} [\mathbf S_f \cdot (\gamma \nabla \phi)_ % workaround for hugo f] = \int_{V_P} S_\phi(\phi) dV \label{eq:int-diff-conv} \end{equation} $$

时间项的处理

为了表述方便，我们定义 $\phi^o$ 和 $\phi^n$ 为旧时刻和新时刻的值（通常省略 $n$）。以Euler格式为例， $\eqref{eq:int-diff-conv}$ 中的时间项可以写作

$$ \begin{equation} \int_{V_P} \frac{\partial \phi}{\partial t} dV =\frac{\phi^n - \phi^o}{\Delta t} V_P \end{equation} $$

对流项的处理

式 $\eqref{eq:int-diff-conv}$ 中的下标 $f$ 代表面单元上的物理量，这些值通过与相邻单元插值得到。考虑下图所示的两个相邻控制体单元 $V_P$ 和 $V_N$ ，它们共享面 $f$， $\mathbf S_f$ 为面法向向量，其大小为 $f$ 的面积。

(图片来自 OpenFOAM Foundation²)

面心物理量可由体心物理量插值得到，以二阶精度的中心差分格式为例

$$ \begin{equation} \phi_f = \lambda \phi_P + (1 - \lambda) \phi_N \end{equation} $$

这里 $\lambda = \overline{fN} / \overline{PN}$ 。

于是 $\eqref{eq:int-diff-conv}$ 中的对流项可以表示为

$$ \begin{equation} \sum_{f} [\mathbf S_f \cdot (\phi \mathbf u)_ % workaround for hugo f] = \underbrace{\sum_{f} [\lambda (\mathbf S_f \cdot \mathbf u_P)] \phi_P}_ {\text{P 单元}} + \underbrace{\sum_{f} [(1 - \lambda) (\mathbf S_f \cdot \mathbf u_N) \phi_N]}_ % workaround for hugo {\text{与 P 相邻的所有单元}} \label{eq:conv_discre} \end{equation} $$

扩散项的处理

扩散项由于包含面法向 $\mathbf S_f$ 与梯度项 $(\nabla \phi)_f$ 的点乘，处理起来比对流项更为复杂。

先看简单的情况，考虑正交网格， $\mathbf S_f$ 与 $\overline{PN}$ 平行，扩散项可以写作

$$ \begin{equation} -\sum_{f} [\mathbf S_f \cdot (\gamma \nabla \phi)_ % workaround for hugo f] = -\sum_{f} \gamma_f | \mathbf S_f | \frac{\phi_N - \phi_P}{|\overline{PN}|} \end{equation} $$

对于非正交网格，则需对上式进行修正，考虑如下图所示的非正交网格，面法向量表示为 $\mathbf S_f = \Delta + \mathbf k$ ，其中 $\Delta$ 与 $\overline{PN}$ 平行，$\mathbf k$ 与面 $f$ 平行

(图片来自 Hrvoje Jasak)

于是有

$$ \begin{equation} \mathbf S_f \cdot (\nabla \phi)_ % workaround for hugo f = \underbrace{|\Delta| \frac{\phi_N - \phi_P}{|\overline{PN}|}}_ % workaround for hugo {\text{正交部分}} + \underbrace{\mathbf k \cdot (\nabla \phi)_ % workaround for hugo f}_{\text{非正交修正部分}} \end{equation} $$

非正交修正将产生两部分：正交部分和非正交修正部分。OpenFOAM 将正交部分进行隐式处理，将其放入系数矩阵中，而将非正交部分显式处理，放入源项中。这样一来便涉及到两个问题

显式处理的非正交部分采用的是旧迭代步（或时间步）的值，存在信息滞后
显式处理的非正交部分放入源项中，将产生一个负的源项，若太大将导致求解崩溃

对于第一个问题，通常采用迭代方式来减小旧值和新值的差距。如在 system/fvSchemes 里通过指定关键字 nNonOrthogonalCorrectors 来指定求解压力泊松方程时需要进行非正交修正的次数。

对于第二个问题，则通过对放入源项的非正交部分 $\mathbf k \cdot (\nabla \phi)_f$ 施加限制器 $\lambda$ 来保证数值稳定，定义如下

$$ \begin{equation} \mathbf k \cdot (\nabla \phi)_f < \lambda|\mathbf S_f| \frac{\phi_N - \phi_P}{|\overline{PN}|} \end{equation} $$

源项的处理

源项通常随时间和空间变化，较为复杂。为了确保数值求解的稳定性和有界性，通常对其做线性化处理

$$ \begin{equation} S_\phi(\phi) = Su (\phi) + Sp (\phi) \phi \end{equation} $$

$\eqref{eq:int-diff-conv}$ 中的源项可以写作

$$ \begin{equation} \int_{V_P}S_\phi (\phi) dV = Su V_P + Sp V_P \phi_P \label{eq:linearized_source_term} \end{equation} $$

注意到这里为止，对流项和扩散项离散后的物理量（$\phi$ 和 $\mathbf u$）均为新时刻的值，即省略了上角标 $n$ 。

将 $\phi_N$ 和 $\phi_P$ 的系数都相加，同时将时间项中 $\phi^o$ 移至方程右侧，记 $r_P = Su V_P + Sp V_P \phi_P + \frac{\phi^o V_P}{\Delta t}$ ，于是方程 $\eqref{eq:int-diff-conv}$ 变成了

$$ \begin{equation} a_P \phi_P + \sum_f a_N \phi_N = r_P \end{equation} $$

联立所有单元的离散方程得到一个大型稀疏线性方程组，写作

$$ \begin{equation} \mathbf A \pmb \phi = \mathbf R \label{eq:discre_coeffs_matrix} \end{equation} $$

$\mathbf A$ 被称为系数矩阵 ，通常被分为三部分，分别是下三角元素 $\mathbf L$ 、对角元素 $\mathbf D$ 和上三角元素 $\mathbf U$ 。

边界条件的处理

常见的边界条件有两种：一种是给定值的 Dirichlet 边界条件，另一种是给定梯度的 Neumann 边界条件。

对于 Dirichlet 边界，边界上的面单元的值等于给定值： $\phi_f = \phi_b$ 。

对流项

$$ \begin{equation} \mathbf S_f \cdot (\phi \mathbf u)_f = \mathbf S_b \cdot (\phi \mathbf u)_b \end{equation} $$

这项只影响 $\mathbf R$ 。

扩散项

$$ \begin{equation} \mathbf S_f \cdot (\gamma \nabla \phi)_f = \mathbf S_b \cdot (\gamma \nabla \phi)_b = \gamma_b |\mathbf S_b |\frac{\phi_b - \phi_P}{\mathbf d_b} \end{equation} $$

$\phi_P$ 影响 $\mathbf D$，而 $\phi_b$ 影响 $\mathbf R$ 。

对于 Neumann 边界，边界上的面单元的梯度等于给定值：$\frac{\mathbf S_b}{|\mathbf S_b|}(\nabla \phi)_b = g_b$ 。

对流项

边界面上的值 $\phi_b = \phi_P + \mathbf d_b \cdot (\nabla \phi)_b = \phi_P + | \mathbf d_b | g_b$ 。

$$ \begin{equation} \mathbf S_f \cdot (\phi \mathbf u)_f = \mathbf S_b \cdot \mathbf u_b (\phi_P + |\mathbf d_b| g_b) \end{equation} $$

对流项将影响 $\mathbf D$ 和 $\mathbf R$ 。

扩散项

$$ \begin{equation} \mathbf S_f \cdot (\gamma \nabla \phi)_f = \gamma_b |\mathbf S_f| g_b \end{equation} $$

扩散项只影响 $\mathbf R$ 。

代码实现

系数矩阵

描述系数矩阵的类为 lduMatrix 及其子类 fvMatrix，涉及到的主要函数有：

Foam::lduMatrix<Type>::diag()
Foam::lduMatrix<Type>::lower()
Foam::lduMatrix<Type>::upper()
Foam::fvMatrix<Type>::D()
Foam::fvMatrix<Type>::A()
Foam::fvMatrix<Type>::source()

lower() 返回下三角元素 $\mathbf L$
diag() 返回不是对角元素 $\mathbf D$，而是未考虑边界条件影响的对角元素
upper() 返回上三角元素 $\mathbf U$

D() 返回的是对角元素 $\mathbf D$，通过 addCmptAvBoundaryDiag 施加边界条件影响

1
2
3
4
5
6
7
8


// OpenFOAM-dev/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C#L703
template<class Type>
Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
{
    tmp<scalarField> tdiag(new scalarField(diag()));
    addCmptAvBoundaryDiag(tdiag.ref());
    return tdiag;
}

A() 返回的不是 $\eqref{eq:discre_coeffs_matrix}$ 中的 $\mathbf A$ ，是 $\mathbf D$ 除以网格体积之后的值。

 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


// OpenFOAM-dev/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C#L736
template<class Type>
Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
{
    tmp<volScalarField> tAphi
    (
        new volScalarField
        (
            IOobject
            (
                "A("+psi_.name()+')',
                psi_.instance(),
                psi_.mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi_.mesh(),
            dimensions_/psi_.dimensions()/dimVol,
            extrapolatedCalculatedFvPatchScalarField::typeName
        )
    );

    tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
    tAphi.ref().correctBoundaryConditions();

    return tAphi;
}

时间项

时间项产生的系数矩阵由 fvm::ddt 系列函数计算得到。以 icoFoam 中使用的 fvm::ddt(U) 为例，对应的实现如下

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


// OpenFOAM-dev/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
template<class Type>
tmp<fvMatrix<Type>>
ddt
(
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return fv::ddtScheme<Type>::New
    (
        vf.mesh(),
        vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
    ).ref().fvmDdt(vf);
}

这里的 New 函数利用 RTS 机制构造字典文件中 ddt(U) 指定的时间离散格式，获得系数矩阵。

对流项

对流项产生的系数矩阵由 fvm::div 系列函数计算得到。同样以 icoFoam 中使用的 fvm::div(phi, U) 为例，对应的实现如下

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11


// OpenFOAM-dev/src/finiteVolume/finiteVolume/fvm/fvmDiv.C#L75
template<class Type>
tmp<fvMatrix<Type>>
div
(
    const surfaceScalarField& flux,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return fvm::div(flux, vf, "div("+flux.name()+','+vf.name()+')');
}

实际上调用的是 fvm::div(phi, U, "div(phi, U)") ，对应的实现如下

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


// OpenFOAM-dev/src/finiteVolume/finiteVolume/fvm/fvmDiv.C#L43
template<class Type>
tmp<fvMatrix<Type>>
div
(
    const surfaceScalarField& flux,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const word& name
)
{
    return fv::convectionScheme<Type>::New
    (
        vf.mesh(),
        flux,
        vf.mesh().divScheme(name)
    )().fvmDiv(flux, vf);
}

和时间项一样，这里也利用 RTS机制，读取字典中 "div(phi, U)" 对应的关键字，构造相应的对流格式，获得系数矩阵。

扩散项

扩散项产生的系数矩阵有 fvm::laplacian 系列函数得到。以 icoFoam 中使用的 fvm::laplacian(nu, U) 为例，对应的实现如下

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


// OpenFOAM-dev/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C#L179
template<class Type, class GType>
tmp<fvMatrix<Type>>
laplacian
(
    const dimensioned<GType>& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const GeometricField<GType, fvsPatchField, surfaceMesh> Gamma
    (
        IOobject
        (
            gamma.name(),
            vf.instance(),
            vf.mesh(),
            IOobject::NO_READ
        ),
        vf.mesh(),
        gamma
    );

    return fvm::laplacian(Gamma, vf);
}

这里做了一次转换，用类型为 dimensionedScalar 的 nu 构造了一个 surfaceScalarField 类型的场量，最终调用了以下实现

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


// OpenFOAM-dev/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C#L206
template<class Type, class GType>
tmp<fvMatrix<Type>>
laplacian
(
    const GeometricField<GType, fvPatchField, volMesh>& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const word& name
)
{
    return fv::laplacianScheme<Type, GType>::New
    (
        vf.mesh(),
        vf.mesh().laplacianScheme(name)
    ).ref().fvmLaplacian(gamma, vf);
}

同样是利用 New 函数，读取字典中 "laplacian(nu, U)" 对应的关键字，构造相应的 laplacian 格式，获得系数矩阵。

源项

式 $\eqref{eq:linearized_source_term}$ 中线性化后的源项右侧有包括 $\phi_P$ 的导数项，需要做隐式处理。OpenFOAM 采用 fvm::su ， fvm::sp 和 fvm::susp 将源项的影响施加到系数矩阵上。其中 fvm::su 用来处理纯显式项 $Su V_P$，fvm::sp 用来处理纯隐式项 $Sp V_P \phi_P$，fvm::susp 则用来处理既包含显式项又包含隐式项的情况。

fvm::su 将 su 的值乘以网格体积后施加到 source() 上：

 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


// OpenFOAM/src/finiteVolume/finiteVolume/fvm/fvmSup.C#L32
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Su
(
    const DimensionedField<Type, volMesh>& su,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    tmp<fvMatrix<Type>> tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            dimVol*su.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm.ref();

    fvm.source() -= mesh.V()*su.field();

    return tfvm;
}

fvm::sp 将 sp 的值乘以网格体积后施加到对角元素 diag() 上：

 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


// OpenFOAM/src/finiteVolume/finiteVolume/fvm/fvmSup.C#L98
template<class Type>template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Sp
(
    const volScalarField::Internal& sp,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    tmp<fvMatrix<Type>> tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            dimVol*sp.dimensions()*vf.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm.ref();

    fvm.diag() += mesh.V()*sp.field();

    return tfvm;
}

fvm::susp 较复杂，它为了数值稳定性考虑，把源项里的正负成分（对应源和汇）分开处理，将大于 0 的部分乘以网格体积后施加到对角元素 diag() 上，小于 0 的部分乘以网格体积和 $\phi_P^o$ 后施加到 source() 上：

 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


// OpenFOAM/src/finiteVolume/finiteVolume/fvm/fvmSup.C#L190
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::SuSp
(
    const volScalarField::Internal& susp,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    tmp<fvMatrix<Type>> tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            dimVol*susp.dimensions()*vf.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm.ref();

    fvm.diag() += mesh.V()*max(susp.field(), scalar(0));

    fvm.source() -= mesh.V()*min(susp.field(), scalar(0))
        *vf.primitiveField();

    return tfvm;
}

参考资料

Jasak, H. (2009). Finite Volume Discretisation with Polyhedral Cell Support. ↩︎
OpenFOAM Foundation (2015). OpenFOAM Programmers' Guide Version 3.0.0. ↩︎

Contents

介绍