无标题 （04/21/2021）

未分类

CFD中非结构网格的数据格式与代数方程组的装配

OpenFOAM，Fluent，Star都采取了类似的数据结构，FVM中很多体积分通过散度定理转变为了面积分，在装配代数方程组的时候，只需要考虑存储在面单元上的变量的符号即可。

以不可压流动中的连续方程为例：
$\int_V{\nabla \cdot \vec{U}}\,dV = \int_{\partial V}{\vec{U}_f \cdot d \vec{A}_f} = \sum_f{\vec{U}_f} \cdot \vec{A}_f = \sum_f{\phi}_f$
其中$\vec{A}_f$是控制体的面积矢量，其方向必须向外。（控制体一般就是几何体单元（格心格式），CFX例外，CFX使用节点存储变量，重新构造了以节点为中心的控制体（格点格式））

OpenFOAM中$\phi$存储于面单元，上述积分代码实现为fvc::div(phi).具体计算就是进行面循环，分别将$\phi$加进左右单元。

fvc::div算子源代码

 template<class Type>
 tmp<GeometricField<Type, fvPatchField, volMesh>>
 div
 (
     const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
 )
 {
     return tmp<GeometricField<Type, fvPatchField, volMesh>>
     (
         new GeometricField<Type, fvPatchField, volMesh>
         (
             "div("+ssf.name()+')',
             fvc::surfaceIntegrate(ssf)
         )
     );
 }

fvc::surfaceIntegrate算子源代码

 template<class Type>
 tmp<GeometricField<Type, fvPatchField, volMesh>>
 surfaceIntegrate
 (
     const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
 )
 {
     const fvMesh& mesh = ssf.mesh();
     tmp<GeometricField<Type, fvPatchField, volMesh>> tvf
     (
         new GeometricField<Type, fvPatchField, volMesh>
         (
             IOobject
             (
                 "surfaceIntegrate("+ssf.name()+')',
                 ssf.instance(),
                 mesh,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
             mesh,
             dimensioned<Type>(ssf.dimensions()/dimVol, Zero),
             extrapolatedCalculatedFvPatchField<Type>::typeName
         )
     );
     GeometricField<Type, fvPatchField, volMesh>& vf = tvf.ref();
     surfaceIntegrate(vf.primitiveFieldRef(), ssf);
     vf.correctBoundaryConditions();
     return tvf;
 }

具体实现

 template<class Type>
 void surfaceIntegrate
 (
     Field<Type>& ivf,
     const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
 )
 {
     const fvMesh& mesh = ssf.mesh();
     const labelUList& owner = mesh.owner();            // lduAddr().lowerAddr()
     const labelUList& neighbour = mesh.neighbour();
     const Field<Type>& issf = ssf;
     // 内场面循环，一次循环即可完成计算，如果基于体单元循环，每个面其实循环了两次
     forAll(owner, facei)
     {
         ivf[owner[facei]] += issf[facei];
         ivf[neighbour[facei]] -= issf[facei];
     }
     forAll(mesh.boundary(), patchi)
     {
         const labelUList& pFaceCells =
             mesh.boundary()[patchi].faceCells();
         const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];
         forAll(mesh.boundary()[patchi], facei)
         {
             ivf[pFaceCells[facei]] += pssf[facei];
         }
     }
     ivf /= mesh.Vsc();
 }