1 //============================================================================
4 * \author Michael NDJINGA
7 * \brief Four equation two phase flow drift model
9 //============================================================================
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
15 #ifndef DRIFTMODEL_HXX_
16 #define DRIFTMODEL_HXX_
18 #include "ProblemFluid.hxx"
20 class DriftModel : public ProblemFluid{
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
28 DriftModel( pressureEstimate pEstimate, int dim, bool useDellacherieEOS=true);
29 //! system initialisation
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;
36 void testConservation();
40 * \brief saves every interesting field in a separate file
41 * @param boolean saveAllFields
43 void saveAllFields(bool saveAllFields=true){
44 _saveAllFields=saveAllFields;
47 // Boundary conditions
48 /** \fn setIntletBoundaryCondition
49 * \brief adds a new boundary condition of type Inlet
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
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);
63 // Boundary conditions
64 /** \fn setIntletBoundaryCondition
65 * \brief adds a new boundary condition of type Inlet
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
75 void setInletEnthalpyBoundaryCondition(string groupName,double enthalpie, double v_x=0, double v_y=0, double v_z=0){
76 if(!_useDellacherieEOS)
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");
83 double Temperature =getMixtureTemperature( enthalpie);
84 double concentration=getMixtureConcentration( enthalpie);
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);
90 double getMixtureConcentration(double enthalpie)
91 {//On extrait la concentration à partir du titre thermodynamique
92 double titreThermo=(enthalpie-_hsatl)/(_hsatv-_hsatl);
95 else if(titreThermo>=1)
100 double getMixtureTemperature(double enthalpie)
102 if(!_useDellacherieEOS)
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");
109 //On extrait la température à partir du titre thermodynamique
110 double titreThermo=(enthalpie-_hsatl)/(_hsatv-_hsatl);
112 if(titreThermo<=_precision)
113 return _fluides[1]->getTemperatureFromEnthalpy(enthalpie, 0);
114 else if (titreThermo>1-_precision)
115 return _fluides[0]->getTemperatureFromEnthalpy(enthalpie, 0);
121 /** \fn setIntletPressureBoundaryCondition
122 * \brief adds a new boundary condition of type InletPressure
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
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);
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
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];
148 pressure-=reference_point[1]*_GravityField3d[1];
150 pressure-=reference_point[2]*_GravityField3d[2];
153 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(0,0),vector<double>(0,0),vector<double>(0,0),Temperature,-1,-1,concentration);
156 /** \fn setWallBoundaryCondition
157 * \brief adds a new boundary condition of type Wall
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
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);
170 /** \fn setInnerWallBoundaryCondition
171 * \brief adds a new boundary condition of type InnerWall (special treatment of the inner wall)
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
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);
184 /** \fn computeNewtonVariation
185 * \brief Builds and solves the linear system to obtain the variation Vkp1-Vk in a Newton scheme using primitive variables
188 void computeNewtonVariation();
190 /** \fn iterateTimeStep
191 * \brief calls computeNewtonVariation to perform one Newton iteration and tests the convergence
193 * @return boolean ok is true is the newton iteration gave a physically acceptable result
195 bool iterateTimeStep(bool &ok);
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
203 Field _Vitesse, _VitesseX, _VitesseY, _VitesseZ, _VoidFraction, _Concentration, _Pressure, _mixtureDensity, _Enthalpy, _Temperature, _DensiteLiquide, _DensiteVapeur, _EnthalpieLiquide, _EnthalpieVapeur;
205 bool _useDellacherieEOS;
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);
214 /** \fn convectionMatrices
215 * \brief calcule la matrice de convection de l'etat interfacial entre deux cellules voisinnes */
216 void convectionMatrices();
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);
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);
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);
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
255 void pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj);
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
272 void porosityGradientSourceVector();
275 * \brief matrice de gravite
278 void gravityMatrix(); */
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);
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);
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
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
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);
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"
332 void computeScaling(double offset);
334 // Fonctions utilisant la loi d'etat
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
342 void consToPrim(const double *U, double* V,double porosity=1);
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
351 void primToCons(const double *V, const int &i, double *U, const int &j);
353 /** \fn primToConsJacobianMatrix
354 * \brief computes the jacobian matrix of the cons->prim function
355 * @pram V : primitive vector state
357 void primToConsJacobianMatrix(double *V);
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);
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 */
372 Vector relative_velocity(double c_v, Vector mean_velocity, double rho_m)
374 Vector result(mean_velocity.getNumberOfRows());
375 for(int i=0;i<mean_velocity.getNumberOfRows();i++)
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);
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);
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
401 void getMixturePressureAndTemperature(double c_v, double rho_m, double rhom_em, double& P, double& T);
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
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
416 * @param square of the velocity vector
418 void getDensityDerivatives(double concentration, double pressure, double temperature, double v2);
420 #endif /* DRIFTMODEL_HXX_*/