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)_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)_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 可以是 scalarvector。通过 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.