OpenFOAM类库介绍（三）隐式扩散项

简介

openFOAM中的扩散项(Laplacian)有显式和隐式两种，显式的Laplacian将返回一个场，而隐式的Laplacian将返回一个矩阵fvMatrix。在求解热传导问题时，我们采用隐式的Laplacian。
$$\nabla \cdot \left(\Gamma \nabla T\right)$$
本文将根据笔者的理解，列出隐式的Laplacian所需的主要函数，并给出详细注释。

fvmLaplacian函数

隐式的Laplacian的主要代码是fvmLaplacian()函数，位于gaussLaplacianScheme.C文件内，该函数将完成所有计算流程，代码如下

template<class Type, class GType>//泛型，Type是待求场vf的类型，GType是系数gamma的类型
 tmp<fvMatrix<Type>>//返回一个tmp临时对象，指向一个矩阵，如果求解这个矩阵就可以得到Laplacian方程的解
 gaussLaplacianScheme<Type, GType>::fvmLaplacian
 (
     const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,//gamma,是一个储存在单元面上的场
     const GeometricField<Type, fvPatchField, volMesh>& vf//待求场vf,是一个储存在单元中心的场
 )
 {
     const fvMesh& mesh = this->mesh();//获取网格
 
     const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());//计算单元面单位法向量场，这个场储存在单元面上
 
     const surfaceVectorField SfGamma(mesh.Sf() & gamma);//计算单元面矢量与gamma的内积，这个场储存在单元面上
     const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
     (
         SfGamma & Sn//计算SfGamma与Sn的内积，得到标量场。事实上就是面积乘以gamma
     );
     const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);//计算两个向量场的差值，但是根据实际计算，SfGammaCorr几乎处处为0，不知有何意义。如有知晓，请在评论区留言！
 
     tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected
     (//调用fvmLaplacianUncorrected函数完成矩阵构造，这个函数将在下文解释
         SfGammaSn,
         this->tsnGradScheme_().deltaCoeffs(vf),//deltaCoeffs是1/delta,作
         vf
     );
     fvMatrix<Type>& fvm = tfvm.ref();//拿取矩阵的引用
 
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tfaceFluxCorrection
         = gammaSnGradCorr(SfGammaCorr, vf);//调用gammaSnGradCorr函数来计算面修正量，由于SfGammaCorr是0，不知有何意义。如有知晓，请在评论区留言！
 
     if (this->tsnGradScheme_().corrected())//判断是否需要非正交修正
     {
         tfaceFluxCorrection.ref() +=
             SfGammaSn*this->tsnGradScheme_().correction(vf);//使用correction函数完成非正交修正，并将修正量叠加在tfaceFluxCorrection上
     }
 
     fvm.source() -= mesh.V()*fvc::div(tfaceFluxCorrection())().primitiveField();//使用面积分的方式，将tfaceFluxCorrection转换到每个单元的源项
 
     if (mesh.fluxRequired(vf.name()))
     {
         fvm.faceFluxCorrectionPtr() = tfaceFluxCorrection.ptr();
     }
 
     return tfvm;//返回矩阵
 }

openfoamwiki中给出了一些解释
代码中的deltaCoeffs仅与网格有关，作为矩阵系数的一部分
上述代码中调用了三个子函数：fvmLaplacianUncorrected，gammaSnGradCorr，corrected。接下来将依次介绍

fvmLaplacianUncorrected函数

fvmLaplacianUncorrected函数可以构造非正交修正前的矩阵，在阅读代码之前应该了解openFOAM矩阵的存储格式LDU，fvmLaplacianUncorrected的代码如下：

template<class Type, class GType>
 tmp<fvMatrix<Type>>
 gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
 (
     const surfaceScalarField& gammaMagSf,//gamma*S的标量场，数据储存在面心
     const surfaceScalarField& deltaCoeffs,//系数，数据储存在面心
     const GeometricField<Type, fvPatchField, volMesh>& vf//待求场，数据储存在体心
 )
 {
     tmp<fvMatrix<Type>> tfvm//新建一个矩阵
     (
         new fvMatrix<Type>
         (
             vf,
             deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
         )
     );
     fvMatrix<Type>& fvm = tfvm.ref();//获取矩阵的指针
     
     //内部面对矩阵的贡献
     fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();//上三角元素，每个面上deltaCoeffs与gammaMagSf的乘积。上下三角是对称的。
     fvm.negSumDiag();//计算对角线元素，同一行非对角线元素的相反数

	//边界面对矩阵的贡献
     forAll(vf.boundaryField(), patchi)//遍历每一个patch
     {
         const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];//待求场的边界条件
         const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];//在边界上gamma*S的标量场
         const fvsPatchScalarField& pDeltaCoeffs =
             deltaCoeffs.boundaryField()[patchi];//在边界上的deltaCoeffs场
 
         if (pvf.coupled())//如果是耦合边界条件（特殊情况）
         {
             fvm.internalCoeffs()[patchi] =
                 pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);//计算边界条件给对角线元素的贡献
             fvm.boundaryCoeffs()[patchi] =
                -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);//计算边界条件给源项（右端项）的贡献
         }
         else//如果不是耦合边界条件（更一般的情况）
         {
             fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();//计算边界条件给对角线元素的贡献
             fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();//计算边界条件给源项（右端项）的贡献
         }
     }
 
     return tfvm;//返回矩阵
 }

需要指出的是，边界条件正式通过gradientInternalCoeffs()和gradientBoundaryCoeffs()引入矩阵的。前者会在对角线元素的基础上添加边界贡献，后者会在右端项的基础上添加边界贡献。不同的边界条件计算公式不同，这里列出了Dirichlet边界条件的代码：

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

前者翻译成数学公式为（deltaCoeffs用$Co$表示）
$$-Co$$
后者翻译成数学公式为（边界自变量用$T_b$表示）
$$Co\cdot T_b$$

gammaSnGradCorr函数

调用gammaSnGradCorr函数来计算面修正量，由于SfGammaCorr是0，不知有何意义。如有知晓，请在评论区留言！

 template<class Type, class GType>
 tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>//返回一个场，数据储存在面心
 gaussLaplacianScheme<Type, GType>::gammaSnGradCorr
 (
     const surfaceVectorField& SfGammaCorr,//一个矢量场，数据储存在面心
     const GeometricField<Type, fvPatchField, volMesh>& vf//待求的场
 )
 {
     const fvMesh& mesh = this->mesh();//提取网格
 
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tgammaSnGradCorr
     (//新建一个场，用于返回数据
         GeometricField<Type, fvsPatchField, surfaceMesh>::New
         (
             "gammaSnGradCorr("+vf.name()+')',
             mesh,
             SfGammaCorr.dimensions()
            *vf.dimensions()*mesh.deltaCoeffs().dimensions()
         )
     );
 
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
         tgammaSnGradCorr.ref().replace
         (
             cmpt,
             fvc::dotInterpolate(SfGammaCorr, fvc::grad(vf.component(cmpt)))//将vf求梯度后插值到单元面上，再与SfGammaCorr点乘
         );
     }
 
     return tgammaSnGradCorr;//返回场
 }

corrected函数

corrected并不在gaussLaplacianScheme.C文件内，而在correctedSnGrad.C文件内，用于计算非正交导致的通量，代码如下

template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
 Foam::fv::correctedSnGrad<Type>::correction
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
     const fvMesh& mesh = this->mesh();//提取网格
 
     // 构造一个场，数据储存在面心
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tssf
     (
         GeometricField<Type, fvsPatchField, surfaceMesh>::New
         (
             "snGradCorr("+vf.name()+')',
             mesh,
             vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions()
         )
     );
     GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tssf.ref();
 
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
         ssf.replace
         (
             cmpt,
             correctedSnGrad<typename pTraits<Type>::cmptType>(mesh)
            .fullGradCorrection(vf.component(cmpt))//对vf的每个分量调用fullGradCorrection子函数
         );
     }
 
     return tssf;
 }

fullGradCorrection子函数的定义如下：

 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
 Foam::fv::correctedSnGrad<Type>::fullGradCorrection
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
     const fvMesh& mesh = this->mesh();
 
     // 构造一个场，数据储存在面心
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tssf =
         linear<typename outerProduct<vector, Type>::type>(mesh).dotInterpolate
         (//对于每个面：将vf的梯度线性插值到面心后与网格的非正交矢量作乘积
             mesh.nonOrthCorrectionVectors(),
             gradScheme<Type>::New
             (
                 mesh,
                 mesh.gradScheme("grad(" + vf.name() + ')')
             )().grad(vf, "grad(" + vf.name() + ')')
         );
     tssf.ref().rename("snGradCorr(" + vf.name() + ')');
 
     return tssf;//返回场
 }

最终，corrected得到了每个面上的通量$grad(T)\ast t$，t是非正交修正矢量，仅与网格有关，在构建网格时已经生成。

总结

隐式的Laplacian主要流程为fvmLaplacian函数，而其中又包含三个子函数fvmLaplacianUncorrected，gammaSnGradCorr，corrected。最主要的计算工作（矩阵合成）由fvmLaplacianUncorrected承担，非正交修正由corrected承担。