compressibleFoam ＆求解器编写思路探讨

nuaawqf

各位：
　　最近在研究这个求解器：
　　compressibleFoam，发现里面的求解思路和我之前采用fortran语言是一样的。思路就是按边来循环，分别对internalFace和boundaryFace进行处理，求得单元交界面上的通量。
　　控制方程如下
　　
其中一些量如下
　　 0_1481448837713_2016-12-11 17:30:24屏幕截图.png
　　 1_1481448837713_2016-12-11 17:30:38屏幕截图.png
　　 2_1481448837713_2016-12-11 17:31:04屏幕截图.png
　　下面这个是求内部面无粘通量的函数

/// Loop over internal faces and get face flux from L and R cell center
label leftCell , rightCell;
forAll( mesh.owner() , iface ) {
/// Get the left and right cell index
leftCell = mesh.owner()[iface];
rightCell = mesh.neighbour()[iface]; 
/// Approximate Riemann solver at interface 
scalar lambda = (*fluxSolver)( &rho[leftCell]  , &U[leftCell]  , &p[leftCell] , /// Left
     &rho[rightCell] , &U[rightCell] , &p[rightCell] ,          /// Right
     &massFlux[iface] , &momFlux[iface] , &energyFlux[iface] ,
     &nf[iface] );
/// Multiply with face area to get face flux
massFlux[iface] *= mesh.magSf()[iface];
momFlux[iface] *= mesh.magSf()[iface];
energyFlux[iface] *= mesh.magSf()[iface];
localDt[ leftCell ]  += lambda * mesh.magSf()[iface];
localDt[ rightCell ] += lambda * mesh.magSf()[iface];
}

下面这个是求物面边界条件的边界面上的无粘通量的函数

\*---------------------------------------------------------------------------*/
   /////////////////////////////////////////////////////////////////////////////////
   /// Remember BC's for Gas Dynamics solvers
   /// are based on characteristics and cannot be applied
   /// to individual variables. Rather it should be applied
   /// to the state vector as a whole. OpenFOAM framework
   /// does not allow us to do this without breaking the
   /// convention.
   /////////////////////////////////////////////////////////////////////////////////  
   /// Loop over boundary patches and determine the
   /// type of boundary condition of the patch
   /////////////////////////////////////////////////////////////////////////////////
   forAll ( mesh.boundaryMesh() , ipatch ) {
     scalar bflux[5];
     word BCTypePhysical = mesh.boundaryMesh().physicalTypes()[ipatch];
     word BCType         = mesh.boundaryMesh().types()[ipatch];
     word BCName         = mesh.boundaryMesh().names()[ipatch];
     const UList<label> &bfaceCells = mesh.boundaryMesh()[ipatch].faceCells();
   /////////////////////////////////////////////////////////////////////////////////
                                /// Slip Wall BC ///
   /////////////////////////////////////////////////////////////////////////////////
     if( BCTypePhysical == "slip" || BCTypePhysical == "symmetry" ) {
       forAll( bfaceCells  , iface ) {
         /// Extrapolate wall pressure
         scalar p_e = p[ bfaceCells[iface] ];
         vector normal = nf.boundaryField()[ipatch][iface];
         scalar face_area = mesh.magSf().boundaryField()[ipatch][iface];
         scalar lambda = std::fabs( U[ bfaceCells[iface] ] & normal ) +
                         std::sqrt ( gama * p_e / rho[ bfaceCells[iface] ] );
         momResidue[ bfaceCells[iface] ] += p_e * normal * face_area;
         localDt[ bfaceCells[iface] ]    += lambda * face_area;
       }
     }

之后我们将边界上的通量从owner单元减去，加到nighbour单元去

   /// Clear all residues
   massResidue = dimensionedScalar( "", massResidue.dimensions() , 0.0 );
   momResidue = dimensionedVector( "" , momResidue.dimensions() , vector( 0.0 , 0.0 , 0.0 ) );
   energyResidue = dimensionedScalar( "" , energyResidue.dimensions() ,  0.0 );
   /// Loop over each face and add the flux to left and subtract from the right
   forAll( mesh.owner() , iface ) {
     /// Store left and right cell reference
     leftCell = mesh.owner()[iface];
     rightCell = mesh.neighbour()[iface];
     /// Note that the normal vector to a face will point
     /// from the owner cell to the neighbour cell (L to R)
     /// Add to left cell
     massResidue[leftCell]    += massFlux[iface];
     momResidue[leftCell]     += momFlux[iface];
     energyResidue[leftCell]  += energyFlux[iface];
     /// Subtract from right cell
     massResidue[rightCell]   -= massFlux[iface];
     momResidue[rightCell]    -= momFlux[iface];
     energyResidue[rightCell] -= energyFlux[iface];
   }

原来编程就是一直按照这种思路来的，可是转到OpenFOAM 下一时适应不过来。看到官方的求解器，根本摸不清楚边界条件和有限体积是如何交互的。有人能解释一下这和我上述思路的差别和联系吗？
大家编写自己的求解器的时候都是怎么实现的？希望看透这一切联系的牛人能给予一点提示和点拨。跪谢！

CFD中文网

看到官方的求解器，根本摸不清楚边界条件和有限体积是如何交互的。有人能解释一下这和我上述思路的差别和联系吗？

粗略看了下，这个求解器完全没用openfoam的思路来，采用的他自己的思路，就像你说的，他在开发这个求解器的时候，可能是有其他的代码（比如同组的Fortran）转过来的。并且完全是面向过程的思想。

如果你要用OpenFOAM求解的思路，还是看OpenFOAM的求解器比较好他这个就是把求解的过程自己实现了，然而没有使用OpenFOAM自带的方法。

他的main函数已经写的很清楚了：

int main(int argc, char *argv[])
{
   #include "setRootCase.H"
   #include "createTime.H"
   #include "createMesh.H"
   #include "setInputValues.H"

   #include "createFields.H"
   #include "readFluxScheme.H"
      
   /// Time step loop
   /// Posts the non-blocking send/recv of fields
   long int iter = 0;
   while( runTime.loop() ) {

     /// 构造通量
     #include "constructFaceFlux.H"
     
     /// 有限体积离散
     #include "sumFlux.H"
     
     /// 边界通量修正
     #include "boundaryFlux.H"

     Info << "Iteration = " << ++iter << "  ";

     /// 矩阵计算
     #include "stateUpdateLTS.H"

      Info << " Max residue = " << rhoResidMax << endl;
     /// 输出结果
     runTime.write();
   }
 
   return 0;
}

CFD中文网

compressibleFoam ＆求解器编写思路探讨