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

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

$$ \begin{equation} \int_V \nabla \cdot \mathbf F dV = \oint_S \mathbf n \cdot \mathbf F dS \end{equation} $$

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

标量场的梯度

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

$$ \begin{equation} \nabla \cdot (\mathbf c T) = \mathbf c \cdot \nabla T + \cancel{T \nabla \cdot \mathbf c} \end{equation} $$

对上式积分，有

$$ \begin{equation} \int_V \nabla \cdot (\mathbf c T) dV = \int_V \mathbf c \cdot \nabla T dV \label{eq:div_scalar} \end{equation} $$

对方程 $\eqref{eq:div_scalar}$ 左端项运用高斯定理，并将 $\mathbf c$ 移至积分外，可得

$$ \begin{equation} \int_V \nabla \cdot (\mathbf c T) dV = \oint_S \mathbf n \cdot (\mathbf c T) d S  = \mathbf c \cdot \oint_S \mathbf n T d S \end{equation} $$

将方程 $\eqref{eq:div_scalar}$ 右端项中的 $ \mathbf c $ 移至积分外，有

$$ \begin{equation} \int_V \mathbf c \cdot \nabla T dV = \mathbf c \cdot \int_V \nabla T dV \end{equation} $$

于是有

$$ \begin{equation} \mathbf c \cdot \int_V \nabla T dV = \mathbf c \cdot \oint_S \mathbf n T d S \end{equation} $$

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

$$ \begin{equation} \boxed{\int_V \nabla T dV = \oint_S \mathbf n T d S = \oint_S d \mathbf S T} \end{equation} $$

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

矢量场的梯度

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

$$ \begin{equation} \mathbf a \otimes \mathbf b \equiv \mathbf a \cdot \mathbf b^T \end{equation} $$

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

$$ \begin{equation} \nabla \otimes \mathbf U = \nabla \cdot \mathbf U^T = (\nabla U_1 ; \nabla U_2 ; \nabla U_3) \end{equation} $$

对其积分，可得

$$ \begin{equation} \int_V \nabla \otimes \mathbf U dV = \int_V (\nabla U_1 ; \nabla U_2 ; \nabla U_3) dV = \oint_S (\mathbf n U_1 ; \mathbf n U_2 ; \mathbf n U_3) d S = \oint_S \mathbf n \otimes \mathbf U d S \end{equation} $$

化简后得到：

$$ \begin{equation} \boxed{\int_V \nabla \otimes \mathbf U dV = \oint_S \mathbf n \otimes \mathbf U d S = \oint_S d \mathbf S \otimes \mathbf U} \end{equation} $$

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

有限体积近似

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

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

$$ \begin{equation} \int_V (\mathbf x - \mathbf x_P) dV = 0 \end{equation} $$

假定场量 $ \phi $ 在控制体单元空间内线性变化，控制体内任意一点 $ \mathbf x $ 处的 $ \phi $，记为 $ \phi_\mathbf x$，满足

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

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

对控制体积分，可得

$$ \begin{equation} \int_V \phi_\mathbf x dV = \int_V \left[\phi_P + (\nabla \phi)_P \cdot (\mathbf x - \mathbf x_P)\right] dV = \phi_P V \end{equation} $$

上式表示在控制体内线性变化的物理量 $ \phi $ 对控制体积分等于体心的值 $ \phi_P $ 乘以控制体体积 $ V $，即

$$ \begin{equation} \phi_P = \frac{1}{V}\int_V \phi_\mathbf x dV \end{equation} $$

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

$$ \begin{equation} (\nabla T)_ % workaround for hugo P = \frac{1}{V} \int_V (\nabla T) dV = \frac{1}{V} \oint_S d \mathbf S T = \frac{1}{V} \sum\limits_{face} \mathbf S_{face} T_{face} \end{equation} $$

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

$$ \begin{equation} (\nabla \otimes \mathbf U)_ % workaround for hugo P = \frac{1}{V} \sum\limits_{face} \mathbf S_{face} \otimes \mathbf U_{face} \end{equation} $$

OpenFOAM 的实现

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

操作符重载

利用 operator * 实现标量与矢量相乘，可用于计算 $\mathbf S_{face} T_{face}$ ：

 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 * 实现矢量与矢量相乘，可用于计算 $\mathbf S_{face} \otimes \mathbf U_{face}$：

 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.