1 //============================================================================
4 * \author Michael NDJINGA, Kieu Nguyen
7 * \brief The compressible Navier-Stokes equations
9 //============================================================================
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
15 #ifndef SINGLEPHASE_HXX_
16 #define SINGLEPHASE_HXX_
18 #include "ProblemFluid.hxx"
19 #include "StiffenedGas.hxx"
21 class SinglePhase : public ProblemFluid{
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
30 SinglePhase(phaseType fluid, pressureEstimate pEstimate,int dim,bool useDellacherieEOS=false);
33 * \brief sets the viscosity
34 * @param viscosite : value of the dynamic viscosity
36 void setViscosityConstant( double viscosite ){
37 _fluides[0]->setViscosity(viscosite);
40 //! system initialisation
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();
49 * \brief saves every interesting field in a separate file
50 * @param boolean saveAllFields
52 void saveAllFields(bool saveAllFields=true){
53 _saveAllFields=saveAllFields;
58 /** \fn setIntletBoundaryCondition
59 * \brief adds a new boundary condition of type Inlet
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
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);
71 /** \fn setIntletRotationBoundaryCondition
72 * \brief adds a new boundary condition of type InletRotationVelocity
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
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);
82 /** \fn setIntletPressureBoundaryCondition
83 * \brief adds a new boundary condition of type InletPressure
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
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);
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
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];
107 pressure-=reference_point[1]*_GravityField3d[1];
109 pressure-=reference_point[2]*_GravityField3d[2];
112 _limitField[groupName]=LimitField(InletPressure,pressure,vector<double>(0,0),vector<double>(0,0),vector<double>(0,0),Temperature,-1,-1,-1);
114 /** \fn setWallBoundaryCondition
115 * \brief adds a new boundary condition of type Wall
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
124 void setWallBoundaryCondition(string groupName,double Temperature=-1, double v_x=0, double v_y=0, double v_z=0){
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);
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);
131 /** \fn computeNewtonVariation
132 * \brief Builds and solves the linear system to obtain the variation Vkp1-Vk in a Newton scheme using primitive variables
135 void computeNewtonVariation();
137 /** \fn iterateTimeStep
138 * \brief calls computeNewtonVariation to perform one Newton iteration and tests the convergence
140 * @return boolean ok is true is the newton iteration gave a physically acceptable result
142 bool iterateTimeStep(bool &ok);
145 StiffenedGas getFluidEOS(){ return *dynamic_cast<StiffenedGas*>(_fluides[0]); }
146 double getReferencePressure() { return _Pref; };
147 double getReferenceTemperature() { return _Tref; };
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();
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
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
221 * @param square of the velocity vector
223 void getDensityDerivatives( double pressure, double temperature, double v2);
226 Field _Enthalpy, _Pressure, _Density, _Temperature, _Momentum, _TotalEnergy, _Vitesse, _VitesseX, _VitesseY, _VitesseZ, _MachNumber;
229 #endif /* SINGLEPHASE_HXX_*/