compressibleFoam &求解器编写思路探讨



  • 各位:
      最近在研究这个求解器:
      compressibleFoam,发现里面的求解思路和我之前采用fortran语言是一样的。思路就是按边来循环,分别对internalFace和boundaryFace进行处理,求得单元交界面上的通量。
      控制方程如下
      0_1481448815813_equation1.png
    其中一些量如下
      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 下一时适应不过来。看到官方的求解器,根本摸不清楚边界条件和有限体积是如何交互的。有人能解释一下这和我上述思路的差别和联系吗?
    大家编写自己的求解器的时候都是怎么实现的?希望看透这一切联系的牛人能给予一点提示和点拨。跪谢!



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

    粗略看了下,这个求解器完全没用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 中国 的连接断开,我们正在尝试重连,请耐心等待