OpenFOAM

OpenFOAM 中的边界条件（二）

本篇在上一篇的基础上来解读 OpenFOAM 中的基础边界条件。基础边界条件一般包括三类，一是Dirichlet 边界，二是 Neumann 边界，三是混合 Dirichlet 和 Neumann 的边界。

1. fixedValue

这个是 OpenFOAM 中的 Dirichlet 边界条件。

  • 构造函数
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    template<class Type>
    fixedValueFvPatchField<Type>::fixedValueFvPatchField
    (
        const fvPatch& p,
        const DimensionedField<Type, volMesh>& iF,
        const dictionary& dict
    )
    :
        fvPatchField<Type>(p, iF, dict, true)
    {}

熟悉 OpenFOAM 的人都知道， fixedValue 这个边界条件需要用 value 关键字来指定边界的值。value 这个关键字是通过 DimensionedField 类来处理的。 DimensionedField 这个类将读取 value 关键字对应的场的值用来初始化边界上的值。

  • coefficients
    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
    template<class Type>
    tmp<Field<Type> > fixedValueFvPatchField<Type>::valueInternalCoeffs
    (
        const tmp<scalarField>&
    ) const
    {
        return tmp<Field<Type> >
        (
            new Field<Type>(this->size(), pTraits<Type>::zero)
        );
    }

    template<class Type>
    tmp<Field<Type> > fixedValueFvPatchField<Type>::valueBoundaryCoeffs
    (
        const tmp<scalarField>&
    ) const
    {
        return *this;
    }
    这里 "*this" 表示类本身，即当前边界上的值。这个值在上面的构造函数中进行了初始化，所以，可以理解为 valueBoundaryCoeffs 函数返回的正是关键字 "value" 所对应的值。

    template<class Type>
    tmp<Field<Type> > fixedValueFvPatchField<Type>::gradientInternalCoeffs() const
    {
        return -pTraits<Type>::one*this->patch().deltaCoeffs();
    }

    template<class Type>
    tmp<Field<Type> > fixedValueFvPatchField<Type>::gradientBoundaryCoeffs() const
    {
        return this->patch().deltaCoeffs()*(*this);
    }

$$
\begin{align}
valueInternalCoeffs & = 0 \\
valueBoundaryCoeffs & = value \\
gradientInternalCoeffs & = -delta \\
gradientBoundaryCoeffs & = delta* value
\end{align}
$$

其中 $delta$ 为面心与面所属网格中心的距离的倒数。
从上述系数，可以知道，fixedValue 边界条件对边界的值和梯度值的计算为如下：
$$
\begin{align}
x_p & = value \\
\nabla x_p & = - delta \cdot x_C + delta \cdot value = (value - x_C) \cdot delta
\end{align}
$$
这与预期是一致的。

2. zeroGradient

这个是 OpenFOAM 中的一种特殊的 Neumann 边界条件，即边界的梯度为零。

  • evaluate 函数
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    template<class Type>
    void zeroGradientFvPatchField<Type>::evaluate(const Pstream::commsTypes)
    {
        if (!this->updated())
        {
            this->updateCoeffs();
        }

        fvPatchField<Type>::operator==(this->patchInternalField());
        fvPatchField<Type>::evaluate();
    }

注意，这里的 operator== 与 operator= 的作用是一样的，都是赋值运算，而不是比较。

  • coefficients
    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
    template<class Type>
    tmp<Field<Type> > zeroGradientFvPatchField<Type>::valueInternalCoeffs
    (
        const tmp<scalarField>&
    ) const
    {
        return tmp<Field<Type> >
        (
            new Field<Type>(this->size(), pTraits<Type>::one)
        );
    }

    template<class Type>
    tmp<Field<Type> > zeroGradientFvPatchField<Type>::valueBoundaryCoeffs
    (
        const tmp<scalarField>&
    ) const
    {
        return tmp<Field<Type> >
        (
            new Field<Type>(this->size(), pTraits<Type>::zero)
        );
    }

    template<class Type>
    tmp<Field<Type> > zeroGradientFvPatchField<Type>::gradientInternalCoeffs() const
    {
        return tmp<Field<Type> >
        (
            new Field<Type>(this->size(), pTraits<Type>::zero)
        );
    }

    template<class Type>
    tmp<Field<Type> > zeroGradientFvPatchField<Type>::gradientBoundaryCoeffs() const
    {
        return tmp<Field<Type> >
        (
            new Field<Type>(this->size(), pTraits<Type>::zero)
        );
    }

$$
\begin{align}
valueInternalCoeffs & = 1 \\
valueBoundaryCoeffs & = 0 \\
gradientInternalCoeffs & = 0 \\
gradientBoundaryCoeffs & = 0
\end{align}
$$

从上述系数，可以知道，fixedValue 边界条件对边界的值和梯度值的计算为如下：

$$
\begin{align}
x_p & = x_C \\
\nabla x_p & = 0
\end{align}
$$
这与预期是一致的。

3. fixedGradient

这个是 OpenFOAM 中的 Neumann 边界条件，可以指定边界上的梯度值。

  • 构造函数
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    template<class Type>
    fixedGradientFvPatchField<Type>::fixedGradientFvPatchField
    (
        const fvPatch& p,
        const DimensionedField<Type, volMesh>& iF,
        const dictionary& dict
    )
    :
        fvPatchField<Type>(p, iF, dict),
        gradient_("gradient", dict, p.size())
    {
        evaluate();
    }

需要读取关键字 “gradient” 对应的值来初始化变量 gradient_。

  • evaluate 函数
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    template<class Type>
    void fixedGradientFvPatchField<Type>::evaluate(const Pstream::commsTypes)
    {
        if (!this->updated())
        {
            this->updateCoeffs();
        }

        Field<Type>::operator=
        (
            this->patchInternalField() + gradient_/this->patch().deltaCoeffs()
        );

        fvPatchField<Type>::evaluate();
    }

$$
x_p = x_C + \frac{gradient} {delta}
$$
其中 $delta$ 为面心与面所属网格中心的距离的倒数。

  • coefficients
    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
    template<class Type>
    tmp<Field<Type> > fixedGradientFvPatchField<Type>::valueInternalCoeffs
    (
        const tmp<scalarField>&
    ) const
    {
        return tmp<Field<Type> >(new Field<Type>(this->size(), pTraits<Type>::one));
    }

    template<class Type>
    tmp<Field<Type> > fixedGradientFvPatchField<Type>::valueBoundaryCoeffs
    (
        const tmp<scalarField>&
    ) const
    {
        return gradient()/this->patch().deltaCoeffs();
    }

    template<class Type>
    tmp<Field<Type> > fixedGradientFvPatchField<Type>::
    gradientInternalCoeffs() const
    {
        return tmp<Field<Type> >
        (
            new Field<Type>(this->size(), pTraits<Type>::zero)
        );
    }

    template<class Type>
    tmp<Field<Type> > fixedGradientFvPatchField<Type>::
    gradientBoundaryCoeffs() const
    {
        return gradient();
    }

$$
\begin{align}
valueInternalCoeffs & = 1 \\
valueBoundaryCoeffs & = \tfrac{gradient}{delta} \\
gradientInternalCoeffs & = 0 \\
gradientBoundaryCoeffs & = gradient
\end{align}
$$

4. mixed

这是 OpenFOAM 中的混合边界条件。

  • 构造函数
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    template<class Type>
    mixedFvPatchField<Type>::mixedFvPatchField
    (
        const fvPatch& p,
        const DimensionedField<Type, volMesh>& iF,
        const dictionary& dict
    )
    :
        fvPatchField<Type>(p, iF, dict),
        refValue_("refValue", dict, p.size()),
        refGrad_("refGradient", dict, p.size()),
        valueFraction_("valueFraction", dict, p.size())
    {
        evaluate();
    }

需要读取三个参数。

  • evaluate
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    template<class Type>
    void mixedFvPatchField<Type>::evaluate(const Pstream::commsTypes)
    {
        if (!this->updated())
        {
            this->updateCoeffs();
        }

        Field<Type>::operator=
        (
            valueFraction_*refValue_
          +
            (1.0 - valueFraction_)*
            (
                this->patchInternalField()
              + refGrad_/this->patch().deltaCoeffs()
            )
        );
        fvPatchField<Type>::evaluate();
    }

$$
x_p = valueFraction \cdot refValue + (1-valueFraction) \cdot (x_C + \frac{refGrad}{delta})
$$

  • coefficients
    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
    template<class Type>
    tmp<Field<Type> > mixedFvPatchField<Type>::valueInternalCoeffs
    (
        const tmp<scalarField>&
    ) const
    {
        return Type(pTraits<Type>::one)*(1.0 - valueFraction_);
    }

    template<class Type>
    tmp<Field<Type> > mixedFvPatchField<Type>::valueBoundaryCoeffs
    (
        const tmp<scalarField>&
    ) const
    {
        return
             valueFraction_*refValue_
           + (1.0 - valueFraction_)*refGrad_/this->patch().deltaCoeffs();
    }

    template<class Type>
    tmp<Field<Type> > mixedFvPatchField<Type>::gradientInternalCoeffs() const
    {
        return -Type(pTraits<Type>::one)*valueFraction_*this->patch().deltaCoeffs();
    }

    template<class Type>
    tmp<Field<Type> > mixedFvPatchField<Type>::gradientBoundaryCoeffs() const
    {
        return
            valueFraction_*this->patch().deltaCoeffs()*refValue_
          + (1.0 - valueFraction_)*refGrad_;
    }

$$
\begin{align}
valueInternalCoeffs & = 1-valueFraction \\
valueBoundaryCoeffs & = valueFraction \cdot refValue + (1-valueFraction) \cdot \tfrac{refGrad}{delta} \\
gradientInternalCoeffs & = -valueFraction \cdot delta \\
gradientBoundaryCoeffs & = valueFraction \cdot refValue \cdot delta + (1-valueFraction) \cdot refGrad
\end{align}
$$

附注：本篇中所有的下标 $p$ 都表示当前边界（present boundary patch），下标 $C$ 表示当前边界所属的网格的中心。