CoreFlows/Models/inc/IsothermalTwoFluid.hxx

   1 //============================================================================
   2 /**
   3  * \file DriftModel.hxx
   4  * \author Michael NDJINGA, Kieu Nguyen
   5  * \version 1.0
   6  * \date 15 janv. 2015
   7  * \brief The isothermal two-fluid model
   8  * */
   9 //============================================================================
  10
  11 /*! \class IsothermalTwoFluid IsothermalTwoFluid.hxx "IsothermalTwoFluid.hxx"
  12  *  \brief Isothermal two-fluid model
  13  *  \details The model consists in two phasic mass equations, two phasic momentum equations, see \ref IsothermalPage for more details
  14  */
  15 #ifndef IsothermalTwoFluid_HXX_
  16 #define IsothermalTwoFluid_HXX_
  17
  18 #include "ProblemFluid.hxx"
  19
  20 class IsothermalTwoFluid : public ProblemFluid{
  21 public :
  22         /** \fn IsothermalTwoFluid
  23                          * \brief Constructor for four equation isothermal two-fluid model
  24                          * \param [in] pressureEstimate : \ref around1bar or \ref around155bars
  25                          * \param [in] int : mesh dimension
  26                          *  */
  27         IsothermalTwoFluid(pressureEstimate pEstimate, int dim);
  28         //initialisation du systeme
  29         void initialize();
  30
  31         void testConservation();
  32
  33         void save();
  34
  35         // Boundary conditions
  36         /** \fn setIntletBoundaryCondition
  37                          * \brief adds a new boundary condition of type Inlet
  38                          * \details
  39                          * \param [in] string : the name of the boundary
  40                          * \param [in] double : the value of the vapour volume fraction at the boundary
  41                          * \param [in] vector<double> : the values of the x component of the 2 velocities at the boundary
  42                          * \param [in] vector<double> : the values of the y component of the 2 velocities at the boundary
  43                          * \param [in] vector<double> : the values of the z component of the 2 velocities at the boundary
  44                          * \param [out] void
  45                          *  */
  46         void setInletBoundaryCondition(string groupName,double alpha, vector<double> v_x=vector<double>(3,0), vector<double> v_y=vector<double>(3,0), vector<double> v_z=vector<double>(3,0)){
  47                 _limitField[groupName]=LimitField(Inlet,-1,v_x,v_y,v_z,-1,-1,alpha,-1);
  48         };
  49         /** \fn setIntletPressureBoundaryCondition
  50                          * \brief adds a new boundary condition of type InletPressure
  51                          * \details
  52                          * \param [in] string : the name of the boundary
  53                          * \param [in] double : the value of the vapour volume fraction at the boundary
  54                          * \param [in] double : the value of the vapour Pressure at the boundary
  55                          * \param [out] void
  56                          *  */
  57         void setInletPressureBoundaryCondition(string groupName,double alpha, double Pressure){
  58                 _limitField[groupName]=LimitField(InletPressure,Pressure,vector<double>(0,0),vector<double>(0,0),vector<double>(0,0),-1,-1,alpha,-1);
  59         };
  60         /** \fn setWallBoundaryCondition
  61                          * \brief adds a new boundary condition of type Wall
  62                          * \details
  63                          * \param [in] string : the name of the boundary
  64                          * \param [in] vector<double> : the values of the x component of the 2 velocities at the boundary
  65                          * \param [in] vector<double> : the values of the y component of the 2 velocities at the boundary
  66                          * \param [in] vector<double> : the values of the z component of the 2 velocities at the boundary
  67                          * \param [out] void
  68                          *  */
  69         void setWallBoundaryCondition(string groupName,vector<double> v_x=vector<double>(3,0), vector<double> v_y=vector<double>(3,0), vector<double> v_z=vector<double>(3,0)){
  70                 _limitField[groupName]=LimitField(Wall,-1,v_x,v_y,v_z,-1,-1,-1,-1);
  71         };
  72         /** \fn setIntPressCoeff
  73                          * \brief sets a value for the interfacial pressure default coefficient
  74                          * \details
  75                          * \param [in] double : the value for the interfacial pressure default coefficient
  76                          * \param [out] void
  77                          *  */
  78         void setIntPressCoeff(double delta){
  79                 _intPressCoeff=delta;
  80         }
  81
  82 protected :
  83         Field _Vitesse1,_Vitesse2;
  84         double _Temperature, _internalEnergy1, _internalEnergy2, _guessalpha;
  85         bool _afficheg2press, _afficheg2alpha;
  86         double _intPressCoeff;
  87         //!calcule l'etat de Roe de deux etats
  88         void convectionState( const long &i, const long &j, const bool &IsBord);
  89         //!calcule la matrice de convection de l'etat interfacial entre deux cellules voisinnes
  90         void convectionMatrices();
  91         //!calcule la matrice de diffusion de l'etat interface pour la diffusion
  92         //!Calcule le flux pour un état et une porosité et une normale donnés
  93         Vector convectionFlux(Vector U,Vector V, Vector normale, double porosity);
  94         void diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord);
  95         //!Ajoute au second membre la contribution de la gravite, chgt phase, chauffage et frottement
  96         void sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i);
  97         //!Computes the pressure loss associated to the face ij
  98         void pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj);
  99         //!Computes the contribution of the porosity gradient associated to the face ij to the source term
 100         void porosityGradientSourceVector();
 101         //matrice de gravite
 102         // void gravityMatrix();
 103         //!Calcule la jacobienne de la CL convection
 104         void jacobian(const int &j, string nameOfGroup,double * normale);
 105         //!Calcule la jacobienne de la CL de diffusion
 106         void jacobianDiff(const int &j, string nameOfGroup);
 107         //!Calcule l'etat fictif a� la frontiere
 108         void setBoundaryState(string nameOfGroup, const int &j,double *normale);
 109         //!Ajoute au second membre la contribution de la diffusion
 110         void addDiffusionToSecondMember(const int &i,const int &j,bool isBoundary);
 111         //!Computes the interfacial flux for the VFFC formulation of the staggered upwinding
 112         Vector staggeredVFFCFlux();
 113         //!Compute the corrected interfacial state for lowMach, pressureCorrection and staggered versions of the VFRoe formulation
 114         void applyVFRoeLowMachCorrections(bool isBord, string groupname="");
 115         //!remplit les vecteurs de scaling
 116         void computeScaling(double offset);
 117         //!Calcule les saut de valeurs propres pour la correction entropique
 118         void entropicShift(double* n);
 119
 120         // Functions of equations of states
 121         void consToPrim(const double *Ucons, double* Vprim,double porosity=1);
 122         void primToCons(const double *V, const int &i, double *U, const int &j);
 123         void primToConsJacobianMatrix(double *V);
 124         double ecartPression(double m1,double m2, double alpha, double e1, double e2);
 125         double ecartPressionDerivee(double m1,double m2, double alpha, double e1, double e2);
 126
 127         double intPressDef(double alpha, double ur_n, double rho1, double rho2);
 128 };
 129
 130 #endif /* IsothermalTwoFluid_HXX_ */