Code of the Week #32: type wedge

rheoworks · 2019 年5 月 6 日 16:00

#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 中的函数，返回了 wedgePolyPatch 的 faceT_ 张量和 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()}

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

Code of the Week #32: type wedge

参考文献