CoreFlows/Models/inc/SinglePhase.hxx

   1 //============================================================================
   2 /**
   3  * \file DriftModel.hxx
   4  * \author Michael NDJINGA, Kieu Nguyen
   5  * \version 1.0
   6  * \date 01 janv. 2016
   7  * \brief The compressible Navier-Stokes equations
   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 SINGLEPHASE_HXX_
  16 #define SINGLEPHASE_HXX_
  17
  18 #include "ProblemFluid.hxx"
  19
  20 class SinglePhase : public ProblemFluid{
  21 public :
  22         /** \fn SinglePhase
  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         SinglePhase(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         void testConservation();
  37
  38         void save();
  39
  40         /** \fn setIntletBoundaryCondition
  41          * \brief adds a new boundary condition of type Inlet
  42          * \details
  43          * \param [in] string : the name of the boundary
  44          * \param [in] double : the value of the temperature at the boundary
  45          * \param [in] double : the value of the x component of the velocity at the boundary
  46          * \param [in] double : the value of the y component of the velocity at the boundary
  47          * \param [in] double : the value of the z component of the velocity at the boundary
  48          * \param [out] void
  49          *  */
  50         void setInletBoundaryCondition(string groupName,double Temperature,double v_x=0, double v_y=0, double v_z=0){
  51                 _limitField[groupName]=LimitField(Inlet,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,-1);
  52         };
  53         /** \fn setIntletRotationBoundaryCondition
  54          * \brief adds a new boundary condition of type InletRotationVelocity
  55          * \details
  56          * \param [in] string : the name of the boundary
  57          * \param [in] double : the value of the temperature at the boundary
  58          * \param [in] double : the components of the rotational of the inlet velocity
  59          * \param [out] void
  60          *  */
  61         void setInletRotationBoundaryCondition(string groupName,double Temperature,double omega_x=0, double omega_y=0, double omega_z=0){
  62                 _limitField[groupName]=LimitField(InletRotationVelocity,0,vector<double>(1,omega_x),vector<double>(1,omega_y),vector<double>(1,omega_z),Temperature,-1,-1,-1);
  63         };
  64         /** \fn setIntletPressureBoundaryCondition
  65          * \brief adds a new boundary condition of type InletPressure
  66          * \details
  67          * \param [in] string : the name of the boundary
  68          * \param [in] double : the value of the pressure at the boundary
  69          * \param [in] double : the value of the temperature at the boundary
  70          * \param [in] double : the components of the rotational of the inlet velocity
  71          * \param [out] void
  72          *  */
  73         void setInletPressureBoundaryCondition(string groupName, double pressure,double Temperature,double omega_x=0, double omega_y=0, double omega_z=0){
  74                 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(1,omega_x),vector<double>(1,omega_y),vector<double>(1,omega_z),Temperature,-1,-1,-1);
  75         };
  76         /** \fn setIntletPressureBoundaryCondition
  77          * \brief adds a new boundary condition of type InletPressure taking into account the hydrostatic pressure variations
  78          * \details The pressure is not constant on the boundary but varies linearly with a slope given by the gravity vector
  79          * \param [in] string : the name of the boundary
  80          * \param [in] double : the value of the pressure at the boundary
  81          * \param [in] double : the value of the temperature at the boundary
  82          * \param [in] vector<double> : reference_point position on the boundary where the value Pressure will be imposed
  83          * \param [out] void
  84          *  */
  85         void setInletPressureBoundaryCondition(string groupName, double pressure,double Temperature, vector<double> reference_point){
  86                 /* On the boundary we have P-Pref=rho g\cdot(x-xref) hence P=Pref-g\cdot xref + g\cdot x */
  87                 pressure-=reference_point[0]*_GravityField3d[0];
  88                 if(_Ndim>1){
  89                         pressure-=reference_point[1]*_GravityField3d[1];
  90                         if(_Ndim>2)
  91                                 pressure-=reference_point[2]*_GravityField3d[2];
  92                 }
  93
  94                 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(0,0),vector<double>(0,0),vector<double>(0,0),Temperature,-1,-1,-1);
  95         };
  96         /** \fn setWallBoundaryCondition
  97          * \brief adds a new boundary condition of type Wall
  98          * \details
  99          * \param [in] string : the name of the boundary
 100          * \param [in] double : the value of the temperature at the boundary
 101          * \param [in] double : the value of the x component of the velocity at the boundary
 102          * \param [in] double : the value of the y component of the velocity at the boundary
 103          * \param [in] double : the value of the z component of the velocity at the boundary
 104          * \param [out] void
 105          *  */
 106         void setWallBoundaryCondition(string groupName,double Temperature=-1, double v_x=0, double v_y=0, double v_z=0){
 107                 if(Temperature!=-1)
 108                         _limitField[groupName]=LimitField(Wall,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,-1);
 109                 else
 110                         _limitField[groupName]=LimitField(Wall,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),getReferenceTemperature(),-1,-1,-1);
 111         };
 112
 113         /** \fn computeNewtonVariation
 114          * \brief Builds and solves the linear system to obtain the variation Vkp1-Vk in a Newton scheme using primitive variables
 115          * @param
 116          * */
 117         void computeNewtonVariation();
 118
 119         /** \fn iterateTimeStep
 120          * \brief calls computeNewtonVariation to perform one Newton iteration and tests the convergence
 121          * @param
 122          * @return boolean ok is true is the newton iteration gave a physically acceptable result
 123          * */
 124         bool iterateTimeStep(bool &ok);
 125
 126         double getReferencePressure()    { return _Pref; };
 127         double getReferenceTemperature() { return _Tref; };
 128
 129 protected :
 130         Field _Vitesse;
 131         double  _drho_sur_dp,   _drho_sur_dT;//derivatives of the density rho wrt cv, p, T
 132         double  _drhoE_sur_dp,  _drhoE_sur_dT;//derivatives of the total energy rho E wrt cv, p, T
 133         bool _useDellacherieEOS;
 134         double _Tref; //EOS reference temperature
 135         double _Pref; //EOS reference pressure
 136
 137         //!calcule l'etat de Roe de deux etats
 138         void convectionState( const long &i, const long &j, const bool &IsBord);
 139         //!calcule la matrice de convection de l'etat interfacial entre deux cellules voisinnes
 140         void convectionMatrices();
 141         //!Calcule le flux pour un état et une porosité et une normale donnés
 142         Vector convectionFlux(Vector U,Vector V, Vector normale, double porosity);
 143         //!calcule la matrice de diffusion de l'etat interface pour la diffusion
 144         void diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord);
 145         //!Computes the source vector associated to the cell i
 146         void sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i);
 147         //!Computes the pressure loss associated to the face ij
 148         void pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj);
 149         //!Computes the contribution of the porosity gradient associated to the face ij to the source term
 150         void porosityGradientSourceVector();
 151         //!Calcule la jacobienne de la CL convection
 152         void jacobian(const int &j, string nameOfGroup,double * normale);
 153         //!Calcule la jacobienne de la CL de diffusion
 154         void jacobianDiff(const int &j, string nameOfGroup);
 155         //!Calcule l'etat fictif a la frontiere
 156         void setBoundaryState(string nameOfGroup, const int &j,double *normale);// delete &nf Kieu
 157         //!Adds the contribution of diffusion to the RHS
 158         void addDiffusionToSecondMember(const int &i,const int &j,bool isBord);
 159         //!Computation of the Roe matrix
 160         void RoeMatrixConservativeVariables(double u_n, double total_enthalpy,Vector velocity, double k, double K);
 161         void convectionMatrixPrimitiveVariables( double rho, double u_n, double H,Vector velocity);
 162         //!Computation of the staggered Roe upwinding matrix in conservative variables
 163         void staggeredRoeUpwindingMatrixConservativeVariables(  double u_n, double total_enthalpy, Vector velocity, double k, double K);
 164         //!Computation of the staggered Roe upwinding matrix in primitive variables
 165         void staggeredRoeUpwindingMatrixPrimitiveVariables(double density, double u_n,double total_enthalpy, Vector velocity);
 166         /**** staggered VFFC scheme has some specific features (no need for Roe matrix), hence some specific functions ******/
 167         //!Computes the interfacial flux for the VFFC formulation of the staggered upwinding
 168         Vector staggeredVFFCFlux();
 169         //!Computes the matrices A^+ and A^- for the VFFC formulation of the staggered upwinding
 170         void staggeredVFFCMatricesConservativeVariables(double u_n);
 171         //!Computes the matrices A^+ and A^- for the VFFC formulation of the staggered upwinding using Primitive Variables
 172         void staggeredVFFCMatricesPrimitiveVariables(double u_n);
 173         //!Compute the corrected interfacial state for lowMach, pressureCorrection and staggered versions of the VFRoe formulation
 174         void applyVFRoeLowMachCorrections(bool isBord, string groupname="");
 175         //!Special preconditioner based on a matrix scaling strategy
 176         void computeScaling(double offset);
 177         //!Calcule les saut de valeurs propres pour la correction entropique
 178         void entropicShift(double* n);
 179         // Fonctions utilisant la loi d'etat
 180         void consToPrim(const double *Ucons, double* Vprim,double porosity=1);
 181         void primToCons(const double *V, const int &i, double *U, const int &j);
 182         void primToConsJacobianMatrix(double *V);
 183         /** \fn getDensityDerivatives
 184          * \brief Computes the partial derivatives of rho, and rho E with regard to the primitive variables  p and  T
 185          * @param pressure
 186          * @param temperature
 187          * @param square of the velocity vector
 188         */
 189         void getDensityDerivatives( double pressure, double temperature, double v2);
 190
 191 };
 192 #endif /* SINGLEPHASE_HXX_*/