OpenFOAM 中的梯度计算：高斯法

OpenFOAM 采用同位网格离散计算域，速度和压力等物理量均存储在有限控制体单元的体心。对于同位网格，高斯法是最常用的计算梯度的方法。

高斯定理可以将体积分转化为封闭曲面积分：

$\begin{matrix} (1) & \int_{V} \nabla \cdot F d V = \oint_{S} n \cdot F d S \end{matrix}$

其中， $F$ 为任意矢量， $V$ 为空间体积， $S$ 为空间 $V$ 的表面， $n$ 为 $d S$ 的单位法向量。

标量场的梯度

对于标量场 $T$ ，要计算其梯度 $\nabla T$ ，可以假设一任意常矢量 $c$ ，有以下公式：

$\begin{matrix} (2) & \nabla \cdot (c T) = c \cdot \nabla T + T \nabla \cdot c \end{matrix}$

对上式积分，有

$\begin{matrix} (3) & \int_{V} \nabla \cdot (c T) d V = \int_{V} c \cdot \nabla T d V \end{matrix}$

对方程 $(3)$ 左端项运用高斯定理，并将 $c$ 移至积分外，可得

$\begin{matrix} (4) & \int_{V} \nabla \cdot (c T) d V = \oint_{S} n \cdot (c T) d S = c \cdot \oint_{S} n T d S \end{matrix}$

将方程 $(3)$ 右端项中的 $c$ 移至积分外，有

$\begin{matrix} (5) & \int_{V} c \cdot \nabla T d V = c \cdot \int_{V} \nabla T d V \end{matrix}$

于是有

$\begin{matrix} (6) & c \cdot \int_{V} \nabla T d V = c \cdot \oint_{S} n T d S \end{matrix}$

由于 $c$ 是任意常矢量，所以我们可以得到：

$\begin{matrix} (7) & \int_{V} \nabla T d V = \oint_{S} n T d S = \oint_{S} d S T \end{matrix}$

可见对标量场梯度的体积分可以通过高斯定理转化为面积分。

矢量场的梯度

对于矢量场 $U$ ，其梯度可以表示为 $\nabla U$ 或 $\nabla \otimes U$ 。其中 $\otimes$ 为张量积符号，其定义如下：

$\begin{matrix} (8) & a \otimes b \equiv a \cdot b^{T} \end{matrix}$

于是 $\nabla \otimes U$ 可以表示为

$\begin{matrix} (9) & \nabla \otimes U = \nabla \cdot U^{T} = (\nabla U_{1}; \nabla U_{2}; \nabla U_{3}) \end{matrix}$

对其积分，可得

$\begin{matrix} (10) & \int_{V} \nabla \otimes U d V = \int_{V} (\nabla U_{1}; \nabla U_{2}; \nabla U_{3}) d V = \oint_{S} (n U_{1}; n U_{2}; n U_{3}) d S = \oint_{S} n \otimes U d S \end{matrix}$

化简后得到：

$\begin{matrix} (11) & \int_{V} \nabla \otimes U d V = \oint_{S} n \otimes U d S = \oint_{S} d S \otimes U \end{matrix}$

可见对矢量场梯度的体积分也可以通过高斯定理转化为面积分。

有限体积近似

上面我们得到了标量场和矢量场梯度的体积分，而有限体积法将连续的空间域划分成有限数量个控制体，如何进一步计算梯度在体心处的值？这里可以采用有限体积近似中的线性假设。

对于任意有限控制体 $V$ ，其体心（形心） $P$ 的定义满足：

$\begin{matrix} (12) & \int_{V} (x - x_{P}) d V = 0 \end{matrix}$

假定场量 $ϕ$ 在控制体单元空间内线性变化，控制体内任意一点 $x$ 处的 $ϕ$ ，记为 $ϕ_{x}$ ，满足

$\begin{matrix} (13) & ϕ_{x} \approx ϕ_{P} + (\nabla ϕ)_{P} \cdot (x - x_{P}) \end{matrix}$

上式可以看作略去二阶及更高阶项的泰勒展开形式，因此具有二阶精度。

对控制体积分，可得

$\begin{matrix} (14) & \int_{V} ϕ_{x} d V = \int_{V} [ϕ_{P} + (\nabla ϕ)_{P} \cdot (x - x_{P})] d V = ϕ_{P} V \end{matrix}$

上式表示在控制体内线性变化的物理量 $ϕ$ 对控制体积分等于体心的值 $ϕ_{P}$ 乘以控制体体积 $V$ ，即

$\begin{matrix} (15) & ϕ_{P} = \frac{1}{V} \int_{V} ϕ_{x} d V \end{matrix}$

因此，标量 $T$ 在体心处的梯度可以表示为：

$\begin{matrix} (16) & (\nabla T)_{P} = \frac{1}{V} \int_{V} (\nabla T) d V = \frac{1}{V} \oint_{S} d S T = \frac{1}{V} \sum_{f a c e} S_{f a c e} T_{f a c e} \end{matrix}$

同理，矢量 $U$ 在控制体体心处的梯度可以表示为：

$\begin{matrix} (17) & (\nabla \otimes U)_{P} = \frac{1}{V} \sum_{f a c e} S_{f a c e} \otimes U_{f a c e} \end{matrix}$

OpenFOAM 的实现

注意到 $T$ 和 $U$ 的梯度计算具有相同的形式，因此 OpenFOAM 采用操作符重载（operator overloading）和类模板（class template）实现泛型编程（generic programming），减少代码重复率。下面以 OpenFOAM 7 为例说明。

操作符重载

利用 operator * 实现标量与矢量相乘，可用于计算 $S_{f a c e} T_{f a c e}$ ：

 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


// src/OpenFOAM/primitives/VectorSpace/VectorSpaceI.H
template<class Form, class Cmpt, int nCmpt>
inline Form operator*
(
    scalar s,
    const VectorSpace<Form, Cmpt, nCmpt>& vs
)
{
    Form v;
    VectorSpaceOps<nCmpt,0>::opSV(v, s, vs, multiplyOp3<Cmpt, scalar, Cmpt>());
    return v;
}


template<class Form, class Cmpt, int nCmpt>
inline Form operator*
(
    const VectorSpace<Form, Cmpt, nCmpt>& vs,
    scalar s
)
{
    Form v;
    VectorSpaceOps<nCmpt,0>::opVS(v, vs, s, multiplyOp3<Cmpt, Cmpt, scalar>());
    return v;
}

利用 operator * 实现矢量与矢量相乘，可用于计算 $S_{f a c e} \otimes U_{f a c e}$ ：

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


// src/OpenFOAM/primitives/Tensor/TensorI.H
template<class Cmpt>
inline typename outerProduct<Vector<Cmpt>, Vector<Cmpt> >::type
operator*(const Vector<Cmpt>& v1, const Vector<Cmpt>& v2)
{
    return Tensor<Cmpt>
    (
        v1.x()*v2.x(), v1.x()*v2.y(), v1.x()*v2.z(),
        v1.y()*v2.x(), v1.y()*v2.y(), v1.y()*v2.z(),
        v1.z()*v2.x(), v1.z()*v2.y(), v1.z()*v2.z()
    );
}

类模板

通过 calcGrad 函数计算体心梯度：

 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


// src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
template<class Type>
Foam::tmp
<
    Foam::GeometricField
    <
        typename Foam::outerProduct<Foam::vector, Type>::type,
        Foam::fvPatchField,
        Foam::volMesh
    >
>
Foam::fv::gaussGrad<Type>::calcGrad
(
    const GeometricField<Type, fvPatchField, volMesh>& vsf,
    const word& name
) const
{
    typedef typename outerProduct<vector, Type>::type GradType;

    tmp<GeometricField<GradType, fvPatchField, volMesh>> tgGrad
    (
        gradf(tinterpScheme_().interpolate(vsf), name)  // 第一个参数为插值得到的面心场
    );
    GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad();

    correctBoundaryConditions(vsf, gGrad);

    return tgGrad;
}

这个函数的类模板参数 Type 可以是 scalar 或 vector。通过 Foam::outerProduct<Foam::vector, Type>::type 得到梯度的类型。

若是 Type = scalar，则其梯度类型为 vector；
若是 Type = vector，则其梯度类型为 tensor。

上面的代码调用了 gradf 函数，计算体心梯度的工作实际上由这个函数完成。在执行 gradf 函数前，已经通过 tinterpScheme_().interpolate(vsf) 将体心场插值计算面心场，将并面心场作为参数 ssf 传给 gradf 函数。然后 gradf 函数将面心值与面法向量相乘、求和，并除以控制体体积，得到体心梯度：

 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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81


// src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
template<class Type>
Foam::tmp
<
    Foam::GeometricField
    <
        typename Foam::outerProduct<Foam::vector, Type>::type,
        Foam::fvPatchField,
        Foam::volMesh
    >
>
Foam::fv::gaussGrad<Type>::gradf
(
    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf, // 插值得到的面心场
    const word& name
)
{
    typedef typename outerProduct<vector, Type>::type GradType;

    const fvMesh& mesh = ssf.mesh();

    tmp<GeometricField<GradType, fvPatchField, volMesh>> tgGrad
    (
        new GeometricField<GradType, fvPatchField, volMesh>
        (
            IOobject
            (
                name,
                ssf.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensioned<GradType>
            (
                "0",
                ssf.dimensions()/dimLength,
                Zero
            ),
            extrapolatedCalculatedFvPatchField<GradType>::typeName
        )
    );
    GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad.ref();

    const labelUList& owner = mesh.owner();
    const labelUList& neighbour = mesh.neighbour();
    const vectorField& Sf = mesh.Sf();

    Field<GradType>& igGrad = gGrad;
    const Field<Type>& issf = ssf;

    forAll(owner, facei)
    {
        GradType Sfssf = Sf[facei]*issf[facei]; // 面法向量乘以面心值

        igGrad[owner[facei]] += Sfssf; // 求和
        igGrad[neighbour[facei]] -= Sfssf; // 求和
    }

    forAll(mesh.boundary(), patchi)
    {
        const labelUList& pFaceCells =
            mesh.boundary()[patchi].faceCells();

        const vectorField& pSf = mesh.Sf().boundaryField()[patchi];

        const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];

        forAll(mesh.boundary()[patchi], facei)
        {
            igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
        }
    }

    igGrad /= mesh.V(); // 除以控制体体积

    gGrad.correctBoundaryConditions();

    return tgGrad;
}

参考资料

Jasak, H. (2009). Finite Volume Discretisation with Polyhedral Cell Support.

Contents