CoreFlows/Models/inc/FiveEqsTwoFluid.hxx

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