CoreFlows/Models/inc/SinglePhaseStaggered.hxx

   1 //============================================================================
   2 /**
   3  * \file DriftModel.hxx
   4  * \author Michael NDJINGA
   5  * \version 1.0
   6  * \date 01 janv. 2018
   7  * \brief The compressible Navier-Stokes equations with an ICE scheme on staggered meshes
   8  * */
   9 //============================================================================
  10
  11 /*! \class SinglePhase SinglePhase.hxx "SinglePhase.hxx"
  12  *  \brief The compressible Navier-Stokes equations
  13  *  \details The model consists in one mass, one momentum and one energy equation, see \ref NSModelsPage for more details
  14  */
  15 #ifndef SINGLEPHASESTAGGERED_HXX_
  16 #define SINGLEPHASESTAGGERED_HXX_
  17
  18 #include "ProblemFluid.hxx"
  19
  20 class SinglePhaseStaggered : public ProblemFluid{
  21 public :
  22         /** \fn SinglePhaseStaggered
  23          * \brief Constructor for the Navier-Stokes system
  24          * \param [in] phaseType : \ref Liquid or \ref Gas
  25          * \param [in] pressureEstimate : \ref around1bar or \ref around155bars
  26          * \param [in] int : mesh dimension
  27          * \param [in] bool : There are two possible equations of state for the fluid
  28          *  */
  29         SinglePhaseStaggered(phaseType fluid, pressureEstimate pEstimate,int dim,bool useDellacherieEOS=false);
  30         //! system initialisation
  31         void initialize();
  32
  33         //fonctions d'echange de flux
  34         //      void getOutputField(const Vec &Flux, const string Champ, const int numBord)=0;//, PetscInt *indices_Flux, PetscInt *indices_Bord, const long range)=0;
  35         //      double trace(const int &numBord, Vec &out)=0;
  36
  37         void save();
  38
  39         /** \fn setIntletBoundaryCondition
  40          * \brief adds a new boundary condition of type Inlet
  41          * \details
  42          * \param [in] string : the name of the boundary
  43          * \param [in] double : the value of the temperature at the boundary
  44          * \param [in] double : the value of the x component of the velocity at the boundary
  45          * \param [in] double : the value of the y component of the velocity at the boundary
  46          * \param [in] double : the value of the z component of the velocity at the boundary
  47          * \param [out] void
  48          *  */
  49         void setInletBoundaryCondition(string groupName,double Temperature,double v_x=0, double v_y=0, double v_z=0){
  50                 _limitField[groupName]=LimitField(Inlet,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,-1);
  51         };
  52         /** \fn setIntletPressureBoundaryCondition
  53          * \brief adds a new boundary condition of type InletPressure
  54          * \details
  55          * \param [in] string : the name of the boundary
  56          * \param [in] double : the value of the pressure at the boundary
  57          * \param [in] double : the value of the temperature at the boundary
  58          * \param [out] void
  59          *  */
  60         void setInletPressureBoundaryCondition(string groupName, double pressure,double Temperature){
  61                 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(0,0),vector<double>(0,0),vector<double>(0,0),Temperature,-1,-1,-1);
  62         };
  63         /** \fn setIntletPressureBoundaryCondition
  64          * \brief adds a new boundary condition of type InletPressure taking into account the hydrostatic pressure variations
  65          * \details The pressure is not constant on the boundary but varies linearly with a slope given by the gravity vector
  66          * \param [in] string : the name of the boundary
  67          * \param [in] double : the value of the pressure at the boundary
  68          * \param [in] double : the value of the temperature at the boundary
  69          * \param [in] vector<double> : reference_point position on the boundary where the value Pressure will be imposed
  70          * \param [out] void
  71          *  */
  72         void setInletPressureBoundaryCondition(string groupName, double pressure,double Temperature, vector<double> reference_point){
  73                 /* On the boundary we have P-Pref=rho g\cdot(x-xref) hence P=Pref-g\cdot xref + g\cdot x */
  74                 pressure-=reference_point[0]*_GravityField3d[0];
  75                 if(_Ndim>1){
  76                         pressure-=reference_point[1]*_GravityField3d[1];
  77                         if(_Ndim>2)
  78                                 pressure-=reference_point[2]*_GravityField3d[2];
  79                 }
  80
  81                 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(0,0),vector<double>(0,0),vector<double>(0,0),Temperature,-1,-1,-1);
  82         };
  83         /** \fn setWallBoundaryCondition
  84          * \brief adds a new boundary condition of type Wall
  85          * \details
  86          * \param [in] string : the name of the boundary
  87          * \param [in] double : the value of the temperature at the boundary
  88          * \param [in] double : the value of the x component of the velocity at the boundary
  89          * \param [in] double : the value of the y component of the velocity at the boundary
  90          * \param [in] double : the value of the z component of the velocity at the boundary
  91          * \param [out] void
  92          *  */
  93         void setWallBoundaryCondition(string groupName,double Temperature,double v_x, double v_y=0, double v_z=0){
  94                 _limitField[groupName]=LimitField(Wall,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,-1);
  95         };
  96
  97         /** \fn computeNewtonVariation
  98          * \brief Builds and solves the linear system to obtain the variation Vkp1-Vk in a Newton scheme using primitive variables
  99          * @param
 100          * */
 101         void computeNewtonVariation();
 102
 103         /** \fn iterateTimeStep
 104          * \brief calls computeNewtonVariation to perform one Newton iteration and tests the convergence
 105          * @param
 106          * @return boolean ok is true is the newton iteration gave a physically acceptable result
 107          * */
 108         bool iterateTimeStep(bool &ok);
 109
 110         void computeVelocityMCells(const Field& velocity,
 111                       Field& velocityMCells)
 112 protected :
 113         Field _Vitesse;
 114         double  _drho_sur_dp,   _drho_sur_dT;//derivatives of the density rho wrt cv, p, T
 115         double  _drhoE_sur_dp,  _drhoE_sur_dT;//derivatives of the total energy rho E wrt cv, p, T
 116         bool _useDellacherieEOS;
 117
 118         //!calcule l'etat de Roe de deux etats
 119         void convectionState( const long &i, const long &j, const bool &IsBord);
 120         //!calcule la matrice de convection de l'etat interfacial entre deux cellules voisinnes
 121         void convectionMatrices();
 122         //!Calcule le flux pour un état et une porosité et une normale donnés
 123         Vector convectionFlux(Vector U,Vector V, Vector normale, double porosity);
 124         //!Computes the source vector associated to the cell i
 125         void sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i);
 126         //!Computes the pressure loss associated to the face ij
 127         void pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj);
 128         //!Computes the contribution of the porosity gradient associated to the face ij to the source term
 129         void porosityGradientSourceVector();
 130         //!Calcule la jacobienne de la CL convection
 131         void jacobian(const int &j, string nameOfGroup,double * normale);
 132         //!Calcule l'etat fictif a la frontiere
 133         void setBoundaryState(string nameOfGroup, const int &j,double *normale);// delete &nf Kieu
 134         //!Adds the contribution of diffusion to the RHS
 135         void convectionMatrixPrimitiveVariables( double rho, double u_n, double H,Vector velocity);
 136         /** \fn getDensityDerivatives
 137          * \brief Computes the partial derivatives of rho, and rho E with regard to the primitive variables  p and  T
 138          * @param pressure
 139          * @param temperature
 140          * @param square of the velocity vector
 141         */
 142         void getDensityDerivatives( double pressure, double temperature, double v2);
 143
 144 };
 145 #endif /* SINGLEPHASESTAGGERED_HXX_*/