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 #include "StiffenedGas.hxx"
  20
  21 class SinglePhase : public ProblemFluid{
  22 public :
  23         /** \fn SinglePhase
  24          * \brief Constructor for the Navier-Stokes system
  25          * \param [in] phaseType : \ref Liquid or \ref Gas
  26          * \param [in] pressureEstimate : \ref around1bar or \ref around155bars
  27          * \param [in] int : mesh dimension
  28          * \param [in] bool : There are two possible equations of state for the fluid
  29          *  */
  30         SinglePhase(phaseType fluid, pressureEstimate pEstimate,int dim,bool useDellacherieEOS=false);
  31
  32         /** \fn setViscosity
  33          * \brief sets the viscosity
  34          * @param viscosite : value of the dynamic viscosity
  35          *       * */
  36         void setViscosityConstant( double viscosite ){
  37                 _fluides[0]->setViscosity(viscosite);
  38         };
  39
  40         //! system initialisation
  41         void initialize();
  42
  43         //fonctions d'echange de flux
  44         //      void getOutputField(const Vec &Flux, const string Champ, const int numBord)=0;//, PetscInt *indices_Flux, PetscInt *indices_Bord, const long range)=0;
  45         //      double trace(const int &numBord, Vec &out)=0;
  46         void testConservation();
  47
  48         /** \fn saveAllFields
  49          * \brief saves every interesting field in a separate file
  50          * @param boolean saveAllFields
  51          * */
  52         void saveAllFields(bool saveAllFields=true){
  53                 _saveAllFields=saveAllFields;
  54         }
  55
  56         void save();
  57
  58         /** \fn setIntletBoundaryCondition
  59          * \brief adds a new boundary condition of type Inlet
  60          * \details
  61          * \param [in] string : the name of the boundary
  62          * \param [in] double : the value of the temperature at the boundary
  63          * \param [in] double : the value of the x component of the velocity at the boundary
  64          * \param [in] double : the value of the y component of the velocity at the boundary
  65          * \param [in] double : the value of the z component of the velocity at the boundary
  66          * \param [out] void
  67          *  */
  68         void setInletBoundaryCondition(string groupName,double Temperature,double v_x=0, double v_y=0, double v_z=0){
  69                 _limitField[groupName]=LimitField(Inlet,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,-1);
  70         };
  71         /** \fn setIntletRotationBoundaryCondition
  72          * \brief adds a new boundary condition of type InletRotationVelocity
  73          * \details
  74          * \param [in] string : the name of the boundary
  75          * \param [in] double : the value of the temperature at the boundary
  76          * \param [in] double : the components of the rotational of the inlet velocity
  77          * \param [out] void
  78          *  */
  79         void setInletRotationBoundaryCondition(string groupName,double Temperature,double omega_x=0, double omega_y=0, double omega_z=0){
  80                 _limitField[groupName]=LimitField(InletRotationVelocity,0,vector<double>(1,omega_x),vector<double>(1,omega_y),vector<double>(1,omega_z),Temperature,-1,-1,-1);
  81         };
  82         /** \fn setIntletPressureBoundaryCondition
  83          * \brief adds a new boundary condition of type InletPressure
  84          * \details
  85          * \param [in] string : the name of the boundary
  86          * \param [in] double : the value of the pressure at the boundary
  87          * \param [in] double : the value of the temperature at the boundary
  88          * \param [in] double : the components of the rotational of the inlet velocity
  89          * \param [out] void
  90          *  */
  91         void setInletPressureBoundaryCondition(string groupName, double pressure,double Temperature,double omega_x=0, double omega_y=0, double omega_z=0){
  92                 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(1,omega_x),vector<double>(1,omega_y),vector<double>(1,omega_z),Temperature,-1,-1,-1);
  93         };
  94         /** \fn setIntletPressureBoundaryCondition
  95          * \brief adds a new boundary condition of type InletPressure taking into account the hydrostatic pressure variations
  96          * \details The pressure is not constant on the boundary but varies linearly with a slope given by the gravity vector
  97          * \param [in] string : the name of the boundary
  98          * \param [in] double : the value of the pressure at the boundary
  99          * \param [in] double : the value of the temperature at the boundary
 100          * \param [in] vector<double> : reference_point position on the boundary where the value Pressure will be imposed
 101          * \param [out] void
 102          *  */
 103         void setInletPressureBoundaryCondition(string groupName, double pressure,double Temperature, vector<double> reference_point){
 104                 /* On the boundary we have P-Pref=rho g\cdot(x-xref) hence P=Pref-g\cdot xref + g\cdot x */
 105                 pressure-=reference_point[0]*_GravityField3d[0];
 106                 if(_Ndim>1){
 107                         pressure-=reference_point[1]*_GravityField3d[1];
 108                         if(_Ndim>2)
 109                                 pressure-=reference_point[2]*_GravityField3d[2];
 110                 }
 111
 112                 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(0,0),vector<double>(0,0),vector<double>(0,0),Temperature,-1,-1,-1);
 113         };
 114         /** \fn setWallBoundaryCondition
 115          * \brief adds a new boundary condition of type Wall
 116          * \details
 117          * \param [in] string : the name of the boundary
 118          * \param [in] double : the value of the temperature at the boundary
 119          * \param [in] double : the value of the x component of the velocity at the boundary
 120          * \param [in] double : the value of the y component of the velocity at the boundary
 121          * \param [in] double : the value of the z component of the velocity at the boundary
 122          * \param [out] void
 123          *  */
 124         void setWallBoundaryCondition(string groupName,double Temperature=-1, double v_x=0, double v_y=0, double v_z=0){
 125                 if(Temperature!=-1)
 126                         _limitField[groupName]=LimitField(Wall,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,-1);
 127                 else
 128                         _limitField[groupName]=LimitField(Wall,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),getReferenceTemperature(),-1,-1,-1);
 129         };
 130
 131         /** \fn computeNewtonVariation
 132          * \brief Builds and solves the linear system to obtain the variation Vkp1-Vk in a Newton scheme using primitive variables
 133          * @param
 134          * */
 135         void computeNewtonVariation();
 136
 137         /** \fn iterateTimeStep
 138          * \brief calls computeNewtonVariation to perform one Newton iteration and tests the convergence
 139          * @param
 140          * @return boolean ok is true is the newton iteration gave a physically acceptable result
 141          * */
 142         bool iterateTimeStep(bool &ok);
 143
 144         //EOS functions
 145         StiffenedGas getFluidEOS(){ return *dynamic_cast<StiffenedGas*>(_fluides[0]); }
 146         double getReferencePressure()    { return _Pref; };
 147         double getReferenceTemperature() { return _Tref; };
 148
 149         /* get output fields for postprocessing or coupling */
 150         vector<string> getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement
 151         Field&         getOutputField(const string& nameField );//Renvoie un champs pour le postraitement
 152         Field& getPressureField();
 153         Field& getVelocityField();
 154         Field& getVelocityXField();
 155         Field& getVelocityYField();
 156         Field& getVelocityZField();
 157         Field& getTemperatureField();
 158         Field& getDensityField();
 159         Field& getMomentumField();
 160         Field& getTotalEnergyField();
 161         Field& getEnthalpyField();
 162         Field& getMachNumberField();
 163
 164 protected :
 165         double  _drho_sur_dp,   _drho_sur_dT;//derivatives of the density rho wrt cv, p, T
 166         double  _drhoE_sur_dp,  _drhoE_sur_dT;//derivatives of the total energy rho E wrt cv, p, T
 167         bool _useDellacherieEOS;
 168         double _Tref; //EOS reference temperature
 169         double _Pref; //EOS reference pressure
 170
 171         //!calcule l'etat de Roe de deux etats
 172         void convectionState( const long &i, const long &j, const bool &IsBord);
 173         //!calcule la matrice de convection de l'etat interfacial entre deux cellules voisinnes
 174         void convectionMatrices();
 175         //!Calcule le flux pour un état et une porosité et une normale donnés
 176         Vector convectionFlux(Vector U,Vector V, Vector normale, double porosity);
 177         //!calcule la matrice de diffusion de l'etat interface pour la diffusion
 178         void diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord);
 179         //!Computes the source vector associated to the cell i
 180         void sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i);
 181         //!Computes the pressure loss associated to the face ij
 182         void pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj);
 183         //!Computes the contribution of the porosity gradient associated to the face ij to the source term
 184         void porosityGradientSourceVector();
 185         //!Calcule la jacobienne de la CL convection
 186         void jacobian(const int &j, string nameOfGroup,double * normale);
 187         //!Calcule la jacobienne de la CL de diffusion
 188         void jacobianDiff(const int &j, string nameOfGroup);
 189         //!Calcule l'etat fictif a la frontiere
 190         void setBoundaryState(string nameOfGroup, const int &j,double *normale);// delete &nf Kieu
 191         //!Adds the contribution of diffusion to the RHS
 192         void addDiffusionToSecondMember(const int &i,const int &j,bool isBord);
 193         //!Computation of the Roe matrix
 194         void RoeMatrixConservativeVariables(double u_n, double total_enthalpy,Vector velocity, double k, double K);
 195         void convectionMatrixPrimitiveVariables( double rho, double u_n, double H,Vector velocity);
 196         //!Computation of the staggered Roe upwinding matrix in conservative variables
 197         void staggeredRoeUpwindingMatrixConservativeVariables(  double u_n, double total_enthalpy, Vector velocity, double k, double K);
 198         //!Computation of the staggered Roe upwinding matrix in primitive variables
 199         void staggeredRoeUpwindingMatrixPrimitiveVariables(double density, double u_n,double total_enthalpy, Vector velocity);
 200         /**** staggered VFFC scheme has some specific features (no need for Roe matrix), hence some specific functions ******/
 201         //!Computes the interfacial flux for the VFFC formulation of the staggered upwinding
 202         Vector staggeredVFFCFlux();
 203         //!Computes the matrices A^+ and A^- for the VFFC formulation of the staggered upwinding
 204         void staggeredVFFCMatricesConservativeVariables(double u_n);
 205         //!Computes the matrices A^+ and A^- for the VFFC formulation of the staggered upwinding using Primitive Variables
 206         void staggeredVFFCMatricesPrimitiveVariables(double u_n);
 207         //!Compute the corrected interfacial state for lowMach, pressureCorrection and staggered versions of the VFRoe formulation
 208         void applyVFRoeLowMachCorrections(bool isBord, string groupname="");
 209         //!Special preconditioner based on a matrix scaling strategy
 210         void computeScaling(double offset);
 211         //!Calcule les saut de valeurs propres pour la correction entropique
 212         void entropicShift(double* n);
 213         // Fonctions utilisant la loi d'etat
 214         void consToPrim(const double *Ucons, double* Vprim,double porosity=1);
 215         void primToCons(const double *V, const int &i, double *U, const int &j);
 216         void primToConsJacobianMatrix(double *V);
 217         /** \fn getDensityDerivatives
 218          * \brief Computes the partial derivatives of rho, and rho E with regard to the primitive variables  p and  T
 219          * @param pressure
 220          * @param temperature
 221          * @param square of the velocity vector
 222         */
 223         void getDensityDerivatives( double pressure, double temperature, double v2);
 224
 225         bool _saveAllFields;
 226         Field _Enthalpy, _Pressure, _Density, _Temperature, _Momentum, _TotalEnergy, _Vitesse, _VitesseX, _VitesseY, _VitesseZ, _MachNumber;
 227
 228         };
 229 #endif /* SINGLEPHASE_HXX_*/