OpenFOAM类库介绍（四）对流项

简介

对流项的处理计算流体力学的核心问题，OpenFOAM分别提供了显式对流项和隐式对流项的计算方式。
$$\nabla \cdot \left(\rho UT\right)$$下文将根据笔者的认识，给出函数调用关系以及详细注释

隐式对流项

gaussConvectionScheme.C中的fvmDiv函数包含了计算流程，首先新建一个矩阵，然后考虑内部单元面对矩阵的贡献，最后考虑边界面对矩阵的贡献

template<class Type>
 tmp<fvMatrix<Type>>//返回一个tmp临时对象，这个对象指向矩阵
 gaussConvectionScheme<Type>::fvmDiv
 (
     const surfaceScalarField& faceFlux,//面通量，在CFD求解器中一般为phi
     const GeometricField<Type, fvPatchField, volMesh>& vf//自变量场
 ) const
 {
     tmp<surfaceScalarField> tweights = tinterpScheme_().weights(vf);//根据插值格式，调用weights函数，给出权重场，这个权重场的数值储存在面心
     const surfaceScalarField& weights = tweights();//获取权重场的引用
 
     tmp<fvMatrix<Type>> tfvm
     (//新建一个矩阵，用于输出
         new fvMatrix<Type>
         (
             vf,
             faceFlux.dimensions()*vf.dimensions()
         )
     );
     fvMatrix<Type>& fvm = tfvm.ref();//获取矩阵的引用
 	//内部单元面对矩阵的贡献
     fvm.lower() = -weights.primitiveField()*faceFlux.primitiveField();//矩阵的下三角元素，由权重和通量的乘积构成，注意有一个负号
     fvm.upper() = fvm.lower() + faceFlux.primitiveField();//矩阵的上三角元素，由(1-权重)和通量的乘积构成
     fvm.negSumDiag();//对角线元素是同一行的非对角线元素之和的相反数
 //边界面对矩阵的贡献
     forAll(vf.boundaryField(), patchi)
     {	
     		//提取对应场的边界
         const fvPatchField<Type>& psf = vf.boundaryField()[patchi];
         const fvsPatchScalarField& patchFlux = faceFlux.boundaryField()[patchi];
         const fvsPatchScalarField& pw = weights.boundaryField()[patchi];
 
         fvm.internalCoeffs()[patchi] = patchFlux*psf.valueInternalCoeffs(pw);//边界条件给对角线元素的贡献
         fvm.boundaryCoeffs()[patchi] = -patchFlux*psf.valueBoundaryCoeffs(pw);//边界条件给源项（右端项）的贡献
     }
 
     if (tinterpScheme_().corrected())//如果corrected()是true，需要使用correction()函数修正。比如surfaceInterpolationScheme的派生类skewCorrected就有correction()函数
     {
         fvm += fvc::surfaceIntegrate(faceFlux*tinterpScheme_().correction(vf));
     }
 
     return tfvm;//返回矩阵
 }

上述代码中包含weights()函数，weights()的定义在surfaceInterpolationScheme.H中

 //- Return the interpolation weighting factors for the given field
 virtual tmp<surfaceScalarField> weights
 (
     const GeometricField<Type, fvPatchField, volMesh>&
 ) const = 0;

该函数是一个纯虚函数，必须由surfaceInterpolationScheme的派生类实现，不同的实现方式将决定了不同的插值格式。对于中心差分格式（线性），其定义在linear.H文件中，代码如下

 //- Return the interpolation weighting factors
 tmp<surfaceScalarField> weights
 (
     const GeometricField<Type, fvPatchField, volMesh>&
 ) const
 {
     return this->mesh().surfaceInterpolation::weights();
 }

有趣的是，返回值是来自网格的权重函数weights()，这个函数属于surfaceInterpolation类，具体代码见surfaceInterpolation.C. 而上文中的weights()函数在surfaceInterpolationScheme类中声明，在linear类中定义。它们是不同的。

同理，对于其他差分格式如下风格式(downwind)，其派生类在downwind.H中，代码如下

 //- Return the interpolation weighting factors
 virtual tmp<surfaceScalarField> weights
 (
     const GeometricField<Type, fvPatchField, volMesh>&
 ) const
 {
     return neg(faceFlux_);
 }

显式对流项

gaussConvectionScheme.C中的fvcDiv函数给出了计算显式对流项的流程

 template<class Type>
 tmp<GeometricField<Type, fvPatchField, volMesh>>//返回量是一个场，而不是矩阵，所以称为显式
 gaussConvectionScheme<Type>::fvcDiv
 (
     const surfaceScalarField& faceFlux,//面通量，数据储存在面心
     const GeometricField<Type, fvPatchField, volMesh>& vf//自变量场
 ) const
 {
     tmp<GeometricField<Type, fvPatchField, volMesh>> tConvection
     (
         fvc::surfaceIntegrate(flux(faceFlux, vf))//使用面积分得到源场的散度，也就是对流项的场
     );
 
     tConvection.ref().rename
     (
         "convection(" + faceFlux.name() + ',' + vf.name() + ')'
     );
 
     return tConvection;//返回得到的场
 }

上述代码中调用了位于fvcSurfaceIntegrate.C的surfaceIntegrate函数，这个函数用于求通量围绕控制体的面积积分，代码比较简单，这里不再列出。
上述代码还调用了flux函数，flux函数依然位于gaussConvectionScheme.C，其定义如下：

 template<class Type>
 tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
 gaussConvectionScheme<Type>::flux
 (
     const surfaceScalarField& faceFlux,
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
     return faceFlux*interpolate(faceFlux, vf);
 }

其中还调用了同文件中的interpolate函数，该函数的定义如下：

 template<class Type>
 tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
 gaussConvectionScheme<Type>::interpolate
 (
     const surfaceScalarField&,
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
     return tinterpScheme_().interpolate(vf);
 }

其中又调用了surfaceInterpolationScheme类中的interpolate函数

 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
 Foam::surfaceInterpolationScheme<Type>::interpolate
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
     if (surfaceInterpolation::debug)
     {
         InfoInFunction
             << "Interpolating "
             << vf.type() << " "
             << vf.name()
             << " from cells to faces"
             << endl;
     }
 
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsf
         = interpolate(vf, weights(vf));
 
     if (corrected())
     {
         tsf.ref() += correction(vf);
     }
 
     return tsf;
 }

其中再调用了同文件下interpolate的重载函数（有两个参数），代码如下：

 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
 Foam::surfaceInterpolationScheme<Type>::interpolate
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf,
     const tmp<surfaceScalarField>& tlambdas
 )
 {
     return dotInterpolate(geometricOneField(), vf, tlambdas);
 }

其中再调用了同文件下的dotInterpolate函数，代码如下：

template<class Type>
 template<class SFType>
 Foam::tmp
 <
     Foam::GeometricField//返回一个场的临时对象，指向场
     <
         typename Foam::innerProduct<typename SFType::value_type, Type>::type,//返回类型通过内积推导得到
         Foam::fvsPatchField,
         Foam::surfaceMesh
     >
 >
 Foam::surfaceInterpolationScheme<Type>::dotInterpolate
 (
     const SFType& Sf,//是单位1
     const GeometricField<Type, fvPatchField, volMesh>& vf,//自变量场
     const tmp<surfaceScalarField>& tlambdas//权重场
 )
 {
     if (surfaceInterpolation::debug)
     {
         InfoInFunction
             << "Interpolating "
             << vf.type() << " "
             << vf.name()
             << " from cells to faces "
                "without explicit correction"
             << endl;
     }
 
     typedef typename Foam::innerProduct<typename SFType::value_type, Type>::type
         RetType;
 
     const surfaceScalarField& lambdas = tlambdas();
 
     const Field<Type>& vfi = vf;
     const scalarField& lambda = lambdas;
 
     const fvMesh& mesh = vf.mesh();
     const labelUList& P = mesh.owner();//内部面的owner列表
     const labelUList& N = mesh.neighbour();//内部表面的neighbour列表
 
     tmp<GeometricField<RetType, fvsPatchField, surfaceMesh>> tsf
     (
         GeometricField<RetType, fvsPatchField, surfaceMesh>::New
         (//新建一个场
             "interpolate("+vf.name()+')',
             mesh,
             Sf.dimensions()*vf.dimensions()
         )
     );
     GeometricField<RetType, fvsPatchField, surfaceMesh>& sf = tsf.ref();//提取引用
 
     Field<RetType>& sfi = sf.primitiveFieldRef();//只保留这个场的内部，舍弃边界条件
 
     const typename SFType::Internal& Sfi = Sf();
 
     for (label fi=0; fi<P.size(); fi++)//遍历每个内部面
     {
         sfi[fi] = Sfi[fi] & (lambda[fi]*(vfi[P[fi]] - vfi[N[fi]]) + vfi[N[fi]]);//注意：Sfi[fi]是单位1,所以等效于(lambda[fi]*(vfi[P[fi]] - vfi[N[fi]]) + vfi[N[fi]])，作用是使用加权平均方式得到面上物理量
     }
 
     // Interpolate across coupled patches using given lambdas
 
     typename GeometricField<RetType, fvsPatchField, surfaceMesh>::
         Boundary& sfbf = sf.boundaryFieldRef();//提取边界
 
     forAll(lambdas.boundaryField(), pi)//遍历每个patch
     {
         const fvsPatchScalarField& pLambda = lambdas.boundaryField()[pi];
         const typename SFType::Patch& pSf = Sf.boundaryField()[pi];
         fvsPatchField<RetType>& psf = sfbf[pi];
 
         if (vf.boundaryField()[pi].coupled())//如果是coupled类型边界
         {
             psf =
                 pSf
               & (
                     pLambda*vf.boundaryField()[pi].patchInternalField()
                   + (1.0 - pLambda)*vf.boundaryField()[pi].patchNeighbourField()
                 );//使用边界面的owner和neigbour共同计算边界处的数值，注意：如果不是coupled类型，那么patchNeighbourField()函数不存在
         }
         else//如果不是coupled类型边界
         {
             psf = pSf & vf.boundaryField()[pi];//边界上的数值即为vf的边界
         }
     }
 
     tlambdas.clear();
 
     return tsf;//返回tsf场作为结果
 }

总结

• 隐式对流项的主函数为fvmDiv，需要调用weights函数，不同格式对应不同的weights函数。隐式对流项主要用于方程求解。
• 显式对流项的主函数为fvcDiv，需要调用surfaceIntegrate函数和flux函数，而flux函数调用了三层interpolate函数和一层dotInterpolate，也需要用到weights函数
• 如果需要修改对流项格式，可以从surfaceInterpolationScheme类派生新的类，并重载weights函数