CoreFlows/Models/inc/DriftModel.hxx

   1 //============================================================================
   2 /**
   3  * \file DriftModel.hxx
   4  * \author Michael NDJINGA
   5  * \version 1.0
   6  * \date 01 janv. 2016
   7  * \brief Four equation two phase flow drift model
   8  * */
   9 //============================================================================
  10
  11 /*! \class DriftModel DriftModel.hxx "DriftModel.hxx"
  12  *  \brief Four equation two phase flow drift model
  13  *  \details One total mass equation, one vapour mass equation, one total momentum equation, one total energy equation, see \ref DriftModelPage for more details
  14  */
  15 #ifndef DRIFTMODEL_HXX_
  16 #define DRIFTMODEL_HXX_
  17
  18 #include "ProblemFluid.hxx"
  19
  20 class DriftModel : public ProblemFluid{
  21 public :
  22         /** \fn DriftModel
  23          * \brief Constructor for the four equation drift model
  24          * \param [in] pressureEstimate : \ref around1bar or \ref around155bars
  25          * \param [in] int : mesh dimension
  26          * \param [in] bool : There are two possible equations of state for each phase
  27          *  */
  28         DriftModel( pressureEstimate pEstimate, int dim, bool useDellacherieEOS=true);
  29         //! system initialisation
  30         void initialize();
  31
  32         //fonctions d'echange de flux
  33         //      void getOutputField(const Vec &Flux, const string Champ, const int numBord)=0;//, PetscInt *indices_Flux, PetscInt *indices_Bord, const long range)=0;
  34         //      double trace(const int &numBord, Vec &out)=0;
  35
  36         void testConservation();
  37         void save();
  38
  39         /** \fn saveAllFields
  40          * \brief saves every interesting field in a separate file
  41          * @param boolean saveAllFields
  42          * */
  43         void saveAllFields(bool saveAllFields=true){
  44                 _saveAllFields=saveAllFields;
  45         }
  46
  47         // Boundary conditions
  48         /** \fn setIntletBoundaryCondition
  49          * \brief adds a new boundary condition of type Inlet
  50          * \details
  51          * \param [in] string : the name of the boundary
  52          * \param [in] double : the value of the temperature at the boundary
  53          * \param [in] double : the value of the concentration at the boundary
  54          * \param [in] double : the value of the x component of the velocity at the boundary
  55          * \param [in] double : the value of the y component of the velocity at the boundary
  56          * \param [in] double : the value of the z component of the velocity at the boundary
  57          * \param [out] void
  58          *  */
  59         void setInletBoundaryCondition(string groupName,double Temperature,double concentration, double v_x=0, double v_y=0, double v_z=0){
  60                 _limitField[groupName]=LimitField(Inlet,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,concentration);
  61         };
  62
  63         // Boundary conditions
  64         /** \fn setIntletBoundaryCondition
  65          * \brief adds a new boundary condition of type Inlet
  66          * \details
  67          * \param [in] string : the name of the boundary
  68          * \param [in] double : the value of the enthalpy at the boundary
  69          * \param [in] double : the value of the concentration at the boundary
  70          * \param [in] double : the value of the x component of the velocity at the boundary
  71          * \param [in] double : the value of the y component of the velocity at the boundary
  72          * \param [in] double : the value of the z component of the velocity at the boundary
  73          * \param [out] void
  74          *  */
  75         void setInletEnthalpyBoundaryCondition(string groupName,double enthalpie, double v_x=0, double v_y=0, double v_z=0){
  76                 if(!_useDellacherieEOS)
  77                 {
  78                         cout<<"!!!!!!!!!!!!!!!!!!!!!!!! DriftModel::setInletEnthalpyBoundaryCondition only available for Dellacherie EOS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
  79                         *_runLogFile<<"!!!!!!!!!!!!!!!!!!!!!!!! DriftModel::setInletEnthalpyBoundaryCondition only available for Dellacherie EOS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
  80                         throw CdmathException("DriftModel::setInletEnthalpyBoundaryCondition only available for Dellacherie EOS");
  81                 }
  82                 else{
  83                         double Temperature  =getMixtureTemperature( enthalpie);
  84                         double concentration=getMixtureConcentration( enthalpie);
  85
  86                         _limitField[groupName]=LimitField(Inlet,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,concentration);
  87                 }
  88         };
  89
  90         double getMixtureConcentration(double enthalpie)
  91         {//On extrait la concentration à partir du titre thermodynamique
  92                 double titreThermo=(enthalpie-_hsatl)/(_hsatv-_hsatl);
  93                 if(titreThermo<=0)
  94                         return 0.;
  95                 else if(titreThermo>=1)
  96                         return 1.;
  97                 else
  98                         return  titreThermo;
  99         }
 100         double getMixtureTemperature(double enthalpie)
 101         {
 102                 if(!_useDellacherieEOS)
 103                 {
 104                         cout<<"!!!!!!!!!!!!!!!!!!!!!!!! DriftModel::setInletEnthalpyBoundaryCondition only available for Dellacherie EOS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 105                         *_runLogFile<<"!!!!!!!!!!!!!!!!!!!!!!!! DriftModel::setInletEnthalpyBoundaryCondition only available for Dellacherie EOS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 106                         throw CdmathException("DriftModel::setInletEnthalpyBoundaryCondition only available for Dellacherie EOS");
 107                 }
 108                 else{
 109                         //On extrait la température à partir du titre thermodynamique
 110                 double titreThermo=(enthalpie-_hsatl)/(_hsatv-_hsatl);
 111
 112                 if(titreThermo<=_precision)
 113                         return _fluides[1]->getTemperatureFromEnthalpy(enthalpie, 0);
 114                 else if (titreThermo>1-_precision)
 115                         return _fluides[0]->getTemperatureFromEnthalpy(enthalpie, 0);
 116                 else
 117                         return _Tsat;
 118                 }
 119         }
 120
 121         /** \fn setIntletPressureBoundaryCondition
 122          * \brief adds a new boundary condition of type InletPressure
 123          * \details
 124          * \param [in] string : the name of the boundary
 125          * \param [in] double : the value of the pressure at the boundary
 126          * \param [in] double : the value of the temperature at the boundary
 127          * \param [in] double : the value of the concentration at the boundary
 128          * \param [out] void
 129          *  */
 130         void setInletPressureBoundaryCondition(string groupName, double pressure,double Temperature,double concentration){
 131                 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(0,0),vector<double>(0,0),vector<double>(0,0),Temperature,-1,-1,concentration);
 132         };
 133
 134         /** \fn setIntletPressureBoundaryCondition
 135          * \brief adds a new boundary condition of type InletPressure taking into account the hydrostatic pressure variations
 136          * \details The pressure is not constant on the boundary but varies linearly with a slope given by the gravity vector
 137          * \param [in] string : the name of the boundary
 138          * \param [in] double : the value of the pressure at the boundary
 139          * \param [in] double : the value of the temperature at the boundary
 140          * \param [in] double : the value of the concentration at the boundary
 141          * \param [in] vector<double> : reference_point position on the boundary where the value Pressure will be imposed
 142          * \param [out] void
 143          *  */
 144         void setInletPressureBoundaryCondition(string groupName, double pressure,double Temperature,double concentration, vector<double> reference_point){
 145                 /* On the boundary we have P-Pref=rho g\cdot(x-xref) hence P=Pref-g\cdot xref + g\cdot x */
 146                 pressure-=reference_point[0]*_GravityField3d[0];
 147                 if(_Ndim>1){
 148                         pressure-=reference_point[1]*_GravityField3d[1];
 149                         if(_Ndim>2)
 150                                 pressure-=reference_point[2]*_GravityField3d[2];
 151                 }
 152
 153                 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(0,0),vector<double>(0,0),vector<double>(0,0),Temperature,-1,-1,concentration);
 154         };
 155
 156         /** \fn setWallBoundaryCondition
 157          * \brief adds a new boundary condition of type Wall
 158          * \details
 159          * \param [in] string : the name of the boundary
 160          * \param [in] double : the value of the temperature at the boundary
 161          * \param [in] double : the value of the x component of the velocity at the boundary
 162          * \param [in] double : the value of the y component of the velocity at the boundary
 163          * \param [in] double : the value of the z component of the velocity at the boundary
 164          * \param [out] void
 165          *  */
 166         void setWallBoundaryCondition(string groupName,double Temperature,double v_x, double v_y=0, double v_z=0){
 167                 _limitField[groupName]=LimitField(Wall,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,-1);
 168         };
 169
 170         /** \fn setInnerWallBoundaryCondition
 171          * \brief adds a new boundary condition of type InnerWall (special treatment of the inner wall)
 172          * \details
 173          * \param [in] string : the name of the boundary
 174          * \param [in] double : the value of the temperature at the boundary
 175          * \param [in] double : the value of the x component of the velocity at the boundary
 176          * \param [in] double : the value of the y component of the velocity at the boundary
 177          * \param [in] double : the value of the z component of the velocity at the boundary
 178          * \param [out] void
 179          *  */
 180         void setInnerWallBoundaryCondition(string groupName,double Temperature,double v_x, double v_y=0, double v_z=0){
 181                 _limitField[groupName]=LimitField(InnerWall,-1,vector<double>(1,v_x),vector<double>(1,v_y),vector<double>(1,v_z),Temperature,-1,-1,-1);
 182         };
 183
 184         /** \fn computeNewtonVariation
 185          * \brief Builds and solves the linear system to obtain the variation Vkp1-Vk in a Newton scheme using primitive variables
 186          * @param
 187          * */
 188         void computeNewtonVariation();
 189
 190         /** \fn iterateTimeStep
 191          * \brief calls computeNewtonVariation to perform one Newton iteration and tests the convergence
 192          * @param
 193          * @return boolean ok is true is the newton iteration gave a physically acceptable result
 194          * */
 195         bool iterateTimeStep(bool &ok);
 196
 197 protected :
 198         double _khi, _ksi, _kappa;//mixture pressure derivatives with regard to rhom, cm and rhom_em
 199         double _drho_sur_dcv,   _drho_sur_dp,   _drho_sur_dT;//derivatives of the total density rho wrt cv, p, T
 200         double _drhocv_sur_dcv, _drhocv_sur_dp, _drhocv_sur_dT;//derivatives of the gas partial density rho cv wrt cv, p, T
 201         double _drhoE_sur_dcv,  _drhoE_sur_dp,  _drhoE_sur_dT;//derivatives of the total energy rho E wrt cv, p, T
 202
 203         Field _Vitesse, _VitesseX, _VitesseY, _VitesseZ, _VoidFraction, _Concentration, _Pressure, _mixtureDensity, _Enthalpy, _Temperature, _DensiteLiquide, _DensiteVapeur, _EnthalpieLiquide, _EnthalpieVapeur;
 204         bool _saveAllFields;
 205         bool _useDellacherieEOS;
 206
 207         /** \fn convectionState
 208          * \brief calcule l'etat de Roe de deux etats
 209          * @param i,j sont des entiers qui correspondent aux numeros des cellules à gauche et à droite de l'interface
 210          * @param IsBord est un booléen qui dit si la cellule i est sur le bord
 211          * @return  l'état de Roe de i et j*/
 212         void convectionState( const long &i, const long &j, const bool &IsBord);
 213
 214         /** \fn convectionMatrices
 215          * \brief calcule la matrice de convection de l'etat interfacial entre deux cellules voisinnes */
 216         void convectionMatrices();
 217
 218         /** \fn Flux
 219          * \brief Computes the convection flux F(U) projected on a vector n
 220          * @param U is the conservative variable vector (total density, vapour partial density, total momentum, total energy)
 221          * @param V is the extended primitive variable vector (vapour mass concentration, pressure, mixture velocity, temperature, vapour volume fraction, relative velocity, vapour density, liquid density mixture enthalpy)
 222          * @param normal is a unit vector giving the direction where the convection flux matrix F(U) is projected
 223          * @param porosity is the ration of the volume occupied by the fluid in the cell (default value is 1)
 224          * @return The convection flux projected in the direction given by the normal vector: F(U)*normal */
 225         Vector convectionFlux(Vector U,Vector V, Vector normale, double porosity);
 226
 227         /** \fn diffusionStateAndMatrices
 228          * \brief calcule la matrice de diffusion de l'etat interface pour la diffusion
 229          * @param i,j sont des entiers qui correspondent aux indices  des cellules de gauche et droite respectivement
 230          * @param IsBord: bollean telling if (i,j) is a boundary face
 231          * @return la matrice de diffusion de l'etat interface pour la diffusion calculé */
 232         void diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord);
 233
 234         /** \fn sourceVector
 235          * \brief Ajoute au second membre la contribution de la gravité
 236          * @param Si vecteur de double correspond au Snd membre
 237          * @param Ui vecteur de double contient les variables conservatives
 238          * @param Vi vecteur de double contient les variables primitives
 239          * @param i int representing the cell number. Used for reading power field
 240          * @return Snd membre mis à jour en ajoutant la contribution de la gravité */
 241         void sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i);
 242
 243         //Computes the pressure loss associated to the face ij
 244         /** \fn pressureLossVector
 245          * \brief Computes the contribution of pressure loss terms in the source term
 246          * \Details pure virtual function, overloaded by each model
 247          * @param pressureLoss output vector containing the pressure loss contributions
 248          * @param K, input pressure loss coefficient
 249          * @param Ui input primitive vectors
 250          * @param Vi input conservative vectors
 251          * @param Uj input primitive vectors
 252          * @param Vj input conservative vectors
 253          * @return
 254          * */
 255         void pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj);
 256
 257         //Computes the contribution of the porosity gradient associated to the face ij to the source term
 258         /** \fn porosityGradientSourceVector
 259          * \brief Computes the contribution of the porosity variation in the source term
 260          * \Details pure virtual function, overloaded by each model
 261          * @param porosityGradientVector output vector containing the porosity variation contributions to the source term
 262          * @param Ui input primitive vectors celli
 263          * @param Vi input conservative vectors celli
 264          * @param porosityi input porosity value celli
 265          * @param deltaxi input diameter celli
 266          * @param Uj input primitive vectors cellj
 267          * @param Vj input conservative vectors cellj
 268          * @param porosityj input porosity value cellj
 269          * @param deltaxj input diameter cellj
 270          * @return
 271          * */
 272         void porosityGradientSourceVector();
 273
 274         /* \fn gravityMatrix
 275          * \brief matrice de gravite
 276          * @param
 277          * @return
 278          void gravityMatrix(); */
 279
 280         /** \fn jacobian
 281          * \brief Calcule la jacobienne de la CL convection
 282          * @param j est un entier constant correpond  à l'indice de la cellule sur le bord
 283          * @param nameOfGroup est une chaine de caractère qui correspond au nom du bord
 284          *               * @return Calcule la jacobienne de la CL convection calculé*/
 285         void jacobian(const int &j, string nameOfGroup,double * normale);
 286
 287         /** \fn jacobianDiff
 288          * \brief Calcule la jacobienne de la CL de diffusion
 289          * @param j est un entier constant correpond  à l'indice de la cellule sur le bord
 290          * @param nameOfGroup est une chaine de caractère qui correspond au nom du bord
 291          * @return la jacobienne de la CL de diffusion calculé */
 292         void jacobianDiff(const int &j, string nameOfGroup);
 293
 294         /** \fn setBoundaryState
 295          * \brief Calcule l'etat fictif a la frontière
 296          * @param j est un entier constant correpond  à l'indice de la cellule sur le bord
 297          * @param nameOfGroup est une chaine de caractère qui correspond au nom du bord
 298          * @return l'etat fictif a la frontière calculé */
 299         void setBoundaryState(string nameOfGroup, const int &j,double * normale);// delete &nf Kieu
 300
 301         /** \fn addDiffusionToSecondMember
 302          * \brief Compute the contribution of the diffusion operator to the right hand side of the system
 303          * \Details this function is pure virtual, and overloaded in each physical model class
 304          * @param i left cell number
 305          * @param j right cell number
 306          * @param boolean isBoundary is true for a boundary face (i,j) and false otherwise
 307          * */
 308         void addDiffusionToSecondMember(const int &i,const int &j,bool isBoundary);
 309         //!Computation of the Roe matrix
 310         void RoeMatrixConservativeVariables(double concentration, double umn,double kinetic_energy, double total_mixture_enthalpy,Vector velocity);
 311         void convectionMatrixPrimitiveVariables(double concentration, double mixture_density, double umn, double total_mixture_enthalpy, Vector vitesse);
 312         //!Computation of the staggered Roe upwinding matrix in conservative variables
 313         void staggeredRoeUpwindingMatrixConservativeVariables(double concentration, double umn,double  kinetic_energy, double total_mixture_enthalpy, Vector velocity);
 314         //!Computation of the staggered Roe upwinding matrix in primitive variables
 315         void staggeredRoeUpwindingMatrixPrimitiveVariables(double concentration, double mixture_density, double umn,double total_mixture_enthalpy, Vector velocity);
 316         /**** staggered VFFC scheme has some specific features (no need for Roe matrix), hence some specific functions ******/
 317         //!Computes the interfacial flux for the VFFC formulation of the staggered upwinding
 318         Vector staggeredVFFCFlux();
 319         //!Computes the matrices A^+ and A^- for the VFFC formulation of the staggered upwinding
 320         void staggeredVFFCMatricesConservativeVariables(double u_mn);
 321         //!Computes the matrices A^+ and A^- for the VFFC formulation of the staggered upwinding using Primitive Variables
 322         void staggeredVFFCMatricesPrimitiveVariables(double u_mn);
 323         //!Compute the corrected interfacial state for lowMach, pressureCorrection and staggered versions of the VFRoe formulation
 324         void applyVFRoeLowMachCorrections(bool isBord, string groupname="");
 325         //!Calcule les saut de valeurs propres pour la correction entropique
 326         void entropicShift(double* n);
 327
 328         /** \fn computeScaling
 329          * \brief Special preconditioner based on a matrix scaling strategy
 330          * @param offset est un double, correspond au la plus grande valeure propre "Maxvp"
 331          */
 332         void computeScaling(double offset);
 333
 334         // Fonctions utilisant la loi d'etat
 335
 336         /** \fn consToPrim
 337          * \brief computes the primitive vector state from a conservative vector state
 338          * @param Ucons : conservative variable vector
 339          * @pram Vprim : primitive variable vector
 340          * @param porosity is the porosity coefficient in case of a porous modeling
 341          * */
 342         void consToPrim(const double *U, double* V,double porosity=1);
 343
 344         /** \fn primToCons
 345          * \brief computes the conservative vector state from a primitive vector state
 346          * @param U : conservative variable vector, may contain several states
 347          * @pram V : primitive variable vector, may contain several states
 348          * @param i : index of the conservative state in the vector U
 349          * @param j : index of the primitive state in the vector V
 350          *       */
 351         void primToCons(const double *V, const int &i, double *U, const int &j);
 352
 353         /** \fn primToConsJacobianMatrix
 354          * \brief computes the jacobian matrix of the cons->prim function
 355          * @pram V : primitive vector state
 356          *       */
 357         void primToConsJacobianMatrix(double *V);
 358
 359         /** \fn computeExtendedPrimState
 360          * \brief Computes extended primitive variable vector
 361          * @param V the primitive vector
 362          * @return The vector (vapour mass concentration, pressure, mixture velocity, temperature, vapour volume fraction, relative velocity, vapour density, liquid density mixture enthalpy) */
 363         Vector computeExtendedPrimState(double *V);
 364
 365         /** \fn relative_velocity
 366          * \brief Computes the relative velocity between the two phases
 367          * @param c_v is the vapour mass concentration: alpha_v rho_v/(alpha_v rho_v + alpha_l rho_l)
 368          * @param mean_velocity is the mixture velocity (alpha_v rho_v u_v+ alpha_l rho_l u_l)/(alpha_v rho_v + alpha_l rho_l)
 369          * @param rho_m is the mixture density : alpha_v rho_v + alpha_l rho_l
 370          * @return The relative velocity vector u_v-u_l */
 371
 372         Vector relative_velocity(double c_v, Vector mean_velocity, double rho_m)
 373         {
 374                 Vector result(mean_velocity.getNumberOfRows());
 375                 for(int i=0;i<mean_velocity.getNumberOfRows();i++)
 376                         result(i)=0;
 377                 return result;
 378         }
 379         /** \fn getMixtureTemperature
 380          * \brief Computes the temperature of the two phase mixture from its pressure, mixture density and Gas concentration
 381          * @param c_v is the vapour mass concentration: alpha_v rho_v/(alpha_v rho_v + alpha_l rho_l)
 382          * @param pressureure is the mixture pressure
 383          * @param rho_m is the mixture density : alpha_v rho_v + alpha_l rho_l
 384          * @return The mixture temperature*/
 385         double getMixtureTemperature(double c_v, double rho_m, double pressure);
 386
 387         /** \fn getMixturePressure
 388          * \brief Computes the pressure of the two phase mixture from its temperature, mixture density and Gas concentration
 389          * @param c_v is the vapour mass concentration: alpha_v rho_v/(alpha_v rho_v + alpha_l rho_l)
 390          * @param temperature is the mixture temperature
 391          * @param rho_m is the mixture density : alpha_v rho_v + alpha_l rho_l
 392          * @return The mixture pressure */
 393         double getMixturePressure(double c_v, double rho_m, double temperature);
 394
 395         /** \fn getMixturePressureAndTemperature
 396          * \brief Computes the pressure and temperature of the two phase mixture from its mixture density, Gas concentration and internal energy
 397          * @param c_v is the vapour mass concentration: alpha_v rho_v/(alpha_v rho_v + alpha_l rho_l)
 398          * @param rho_m is the mixture density : alpha_v rho_v + alpha_l rho_l
 399          * @param rhom_em is the mixture internal energy
 400          * @return P and T */
 401         void getMixturePressureAndTemperature(double c_v, double rho_m, double rhom_em, double& P, double& T);
 402
 403         /** \fn getMixturePressureDerivatives
 404          * \brief Computes the partial derivatives khi, ksi and kappa of the pressure of the two phase mixture with regard to the total mass rho, partial gas mass rho c_v and internal energy rho e
 405          * @param m_v is the vapour partial density: alpha_v rho_v)
 406          * @param m_l is the vapour partial density: alpha_l rho_l)
 407          * @param pressure is the mixture pressure
 408          * @param temperature is the mixture temperature
 409          */
 410         void getMixturePressureDerivatives(double m_v, double m_l, double pression, double Tm);
 411         /** \fn getDensityDerivatives
 412          * \brief Computes the partial derivatives of rho, rho c_v and rho E with regard to the primitive variables c_v, p, T
 413          * @param concentration
 414          * @param pressure
 415          * @param temperature
 416          * @param square of the velocity vector
 417          */
 418         void getDensityDerivatives(double concentration, double pressure, double temperature, double v2);
 419 };
 420 #endif /* DRIFTMODEL_HXX_*/
 421
 422
 423