Code of the Week #32: type wedge

#201919

本次,我们来谈谈 wedge 这个特殊的边界。 在 OpenFOAM 中,如果使用二维轴对称,要求在第三个方向使用一个网格的厚度,形成楔形的网格,同时在边界上使用以下设置:

boundaryNAME
{
    type            wedge;
}

对于这样类似的网格的要求是[1]:

  • 一对 wedge 边界须对等地横跨一个坐标平面
  • 分离一对 wedge 边界的网格要求是
    • 仅有一层网格
    • 在中间平面的法向,两边界上的点必须处于同一坐标值
  • 颗粒之类的将归到两个wedge 之间的中心平面

‘wedgeFvPatchField’ 继承自 transformFvPatchField< Type >,在 wedgeFvPatchField.H 中见不到 valueInternalCoeff 之类的,因此这个已经在 src/finiteVolume/fields/fvPatchFields/basic/transform/transformFvPatchField.C 中定义,

 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::transformFvPatchField<Type>::valueInternalCoeffs
 (
     const tmp<scalarField>&
 ) const
 {
     return pTraits<Type>::one - snGradTransformDiag();
 }
 
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::transformFvPatchField<Type>::valueBoundaryCoeffs
 (
     const tmp<scalarField>&
 ) const
 {
     return
         *this
       - cmptMultiply
         (
             valueInternalCoeffs(this->patch().weights()),
             this->patchInternalField()
         );
 }
 
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::transformFvPatchField<Type>::gradientInternalCoeffs() const
 {
     return -this->patch().deltaCoeffs()*snGradTransformDiag();
 }
 
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::transformFvPatchField<Type>::gradientBoundaryCoeffs() const
 {
     return
         snGrad()
       - cmptMultiply(gradientInternalCoeffs(), this->patchInternalField());
 }

snGrad, snGradTransformDiag 函数于 transformFvPatchField.H 中定义了虚函数

             //- Return gradient at boundary
             virtual tmp<Field<Type>> snGrad() const = 0;
 
             //- Return face-gradient transform diagonal
             virtual tmp<Field<Type>> snGradTransformDiag() const = 0;

因此,这两个函数是 wedge 边界的关键。

通过对上面几个函数的考查,可以看到,对于 wedge 边界的几个系数设置如下:

Coefficient function Expression
A_1 valueInternalCoeffs 1-\nabla_f
B_1 valueBoundaryCoeffs -(1-\nabla_f)\phi_c
A_2 gradientInternalCoeffs -\text{deltaCoeffs()}\nabla_f
B_2 gradientBoundaryCoeffs \nabla_{nf}-(-\text{deltaCoeffs()}\nabla_f)\times \phi_c

其中 deltaCoeffs 为边界面到该面所属 cell 的中心之距离的倒数。 \nabla_{nf} 指的是 snGrad(), \nabla_f 指的是 snGradTransformDiag()。从 wedgeFvPatchField.C 中查到,

template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::wedgeFvPatchField<Type>::snGrad() const
{
    const Field<Type> pif(this->patchInternalField());

    return
    (
        transform(refCast<const wedgeFvPatch>(this->patch()).cellT(), pif) - pif
    )*(0.5*this->patch().deltaCoeffs());
}

template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::wedgeFvPatchField<Type>::snGradTransformDiag() const
{
    const diagTensor diagT =
        0.5*diag(I - refCast<const wedgeFvPatch>(this->patch()).cellT());

    const vector diagV(diagT.xx(), diagT.yy(), diagT.zz());

    return tmp<Field<Type> >
    (
        new Field<Type>
        (
            this->size(),
            transformMask<Type>
            (
                pow
                (
                    diagV,
                    pTraits<typename powProduct<vector, pTraits<Type>::rank>
                    ::type>::zero
                )
            )
        )
    );
}

faceT() 为中心平面的法向到 patch 的法向变换的旋转张量。cellT() 为 faceT() 的自身内积。
faceT (和 cellT) 为 wedgeFvPatch.H 中的函数,返回了 wedgePolyPatchfaceT_ 张量和 cellT_。看到 wedgePolyPatch.C 中有

         faceT_ = rotationTensor(centreNormal_, n_);
         cellT_ = faceT_ & faceT_;

rotationTensor 的作用是返回两个平面的旋转张量,见 transform.H

 //- Rotational transformation tensor from vector n1 to n2
 inline tensor rotationTensor
 (
     const vector& n1,
     const vector& n2
 )
 {
     const scalar s = n1 & n2;
     const vector n3 = n1 ^ n2;
     const scalar magSqrN3 = magSqr(n3);
 
     // n1 and n2 define a plane n3
     if (magSqrN3 > SMALL)
     {
         // Return rotational transformation tensor in the n3-plane
         return
             s*I
           + (1 - s)*sqr(n3)/magSqrN3
           + (n2*n1 - n1*n2);
     }
     // n1 and n2 are contradirectional
     else if (s < 0)
     {
         // Return mirror transformation tensor
         return I + 2*n1*n2;
     }
     // n1 and n2 are codirectional
     else
     {
         // Return null transformation tensor
         return I;
     }
 }

如果边界是标量,则返回相同的标量值。若边界是向量,则使中心向量旋转到边界向量。其源代码部分 transform.H

 //- No-op rotational transform of a scalar
 inline scalar transform(const tensor&, const scalar s)
 {
     return s;
 }
 
 
...
 
 
 //- Use rotational tensor to transform a vector.
 //  Same as (rot & v)
 template<class Cmpt>
 inline Vector<Cmpt> transform(const tensor& tt, const Vector<Cmpt>& v)
 {
     return (tt & v);
 }

所以,对于 snGrad() 相应的表达式为

\nabla_{nf}=(\phi_f-\phi_c)\times 0.5\times\text{deltaCoeffs()}

:joy:小作业:写出 snGradTransformDiag() 的表达式

参考文献

  1. https://www.openfoam.com/documentation/guides/latest/doc/guide-bcs-constraint-wedge.html
  2. https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1wedgeFvPatch.html
1 Like