4 //Perfect gas EOS Loi d'etat gaz parfait
5 StiffenedGas::StiffenedGas( double gamma, double cv, double T_ref, double e_ref): Fluide()
8 throw CdmathException("StiffenedGas::StiffenedGas: gamma<1");
16 cout<<"Perfect gas EOS P=(gamma - 1) * rho e with parameter"<< " gamma= " << _gamma<<endl;
17 cout<<"Linearised internal energy law e(T)= e_ref+ cv_ref (T-Tref), around temperature Tref= "<< _Tref<<" K, internal energy e_ref= "<<_e_ref<<" J/Kg, with specific heat cv_ref= "<< _Cv<<" J/Kg/K"<<endl;
19 //Stiffened gas fitted using sound speed
20 StiffenedGas::StiffenedGas(double rho_ref, double p_ref, double T_ref, double e_ref, double c_ref, double cv_ref): Fluide()
23 //_gamma=(1+sqrt(1+4*c_ref*c_ref/(cv_ref*T_ref)))/2;
26 _gamma=1+c_ref*c_ref/(_e_ref+p_ref/rho_ref);
28 throw CdmathException("StiffenedGas::setEOS: gamma<1");
32 _p0=((_gamma-1)*rho_ref*_e_ref-p_ref)/_gamma;
34 cout<<"Stiffened gas EOS P=(gamma - 1) * rho (e(T)-q) - gamma*p0 with parameters"<< " p0= " << _p0 << " gamma= " << _gamma<<endl;
35 cout<<"Calibrated around pressure= "<<p_ref<< " Pa, density= "<< rho_ref<< " internal energy= " << e_ref << " J/Kg, sound speed= " << c_ref<< endl;
36 cout<<"Linearised internal energy law e(T)= e_ref+ cv_ref (T-Tref), around temperature Tref= "<< _Tref<<" K, internal energy e_ref= "<<_e_ref<<" J/Kg, specific heat cv_ref= "<< _Cv<<" J/Kg/K"<<endl;
38 // Loi d'etat stiffened gas S. Dellacherie
39 StiffenedGasDellacherie::StiffenedGasDellacherie( double gamma, double p0, double q, double cv_ref)
42 throw CdmathException("StiffenedGas::StiffenedGas: gamma<1");
51 cout<<"S. Dellacherie Stiffened gas EOS P=(gamma - 1)/gamma * rho (h(T)-q) - p0 with parameters"<< " p0=" << _p0 << " gamma= " << _gamma<< " q= " << _q<< endl;
52 cout<<"Linearised internal energy law h(T)= h_ref+ cp_ref (T-Tref) around temperature "<< _Tref<<" K, enthalpy "<<_h_ref<<" J/Kg,, specific heat cp_ref= "<< _Cp<<" J/Kg/K"<<endl;
55 double StiffenedGas::getInternalEnergy(double T, double rho)
57 return _Cv*(T-_Tref)+_e_ref;
59 double StiffenedGasDellacherie::getInternalEnergy(double T, double rho)
61 double h= getEnthalpy(T);//h=e+p/rho=e+(gamma-1)(e-q)-gamma p0/rho=gamma(e- p0/rho)-(gamma-1)q
62 return (h+(_gamma-1)*_q)/_gamma+_p0/rho;
64 double StiffenedGas::getEnthalpy(double T, double rho)
66 double e=getInternalEnergy( T, rho);
67 return _gamma*(e-_p0/rho)-(_gamma-1)*_q;
69 double StiffenedGasDellacherie::getEnthalpy(double T, double rho)
71 return _Cp*(T-_Tref)+_h_ref;
73 double StiffenedGas::getTemperatureFromPressure(const double p, const double rho)
75 //p=(gamma-1)rho(e-q)-gamma p0
76 double e=_q +(p+_gamma*_p0)/((_gamma-1)*rho);
77 return (e-_e_ref)/_Cv + _Tref;
79 double StiffenedGasDellacherie::getTemperatureFromPressure(const double p, const double rho)
81 //P=(gamma - 1)/gamma * rho (h(T)-q) - p0
82 double h=_q +_gamma*(p+_p0)/((_gamma-1)*rho);
83 return (h-_h_ref)/_Cp + _Tref;
86 double StiffenedGas::getTemperatureFromEnthalpy(const double h, const double rho)
88 //h=e+p/rho et p=(gamma-1)rho(e-q)-gamma p0
89 double e=(h+(_gamma-1)*_q)/_gamma + _p0/rho;
90 return (e-_e_ref)/_Cv + _Tref;
92 double StiffenedGasDellacherie::getTemperatureFromEnthalpy(const double h, const double rho)
94 return (h-_h_ref)/_Cp + _Tref;
97 double StiffenedGas::getDensity(double p, double T)
99 return (p+_gamma*_p0)/((_gamma-1)*(getInternalEnergy(T)-_q));
101 double StiffenedGasDellacherie::getDensity(double p, double T)
103 return _gamma*(p+_p0)/((_gamma-1)*(getEnthalpy(T)-_q));
106 double StiffenedGas::getJumpDensPress(const double e_l, const double e_r){
108 inv_a_2 = (e_l+e_r)/((_gamma-1)*2*e_l*e_r);
111 double StiffenedGas::getJumpDensInternalEnergy(const double p_l,const double p_r,const double e_l,const double e_r){
113 p_inf = - _gamma*_p0;
114 b = (p_inf-0.5*(p_l+p_r))/((_gamma-1)*e_l*e_r);
117 double StiffenedGas::getJumpInternalEnergyTemperature(){
120 // function to compute partial rho/ partial p (e constant), partial rho/partial e (p constant)
121 // use to compute the Jacobian matrix
122 double StiffenedGas::getDiffDensPress(const double e){
124 inv_a_2 = 1/((_gamma-1)*e);
127 double StiffenedGas::getDiffDensInternalEnergy(const double p,const double e){
129 p_inf = - _gamma*_p0;
130 b = (p_inf-p)/((_gamma-1)*e*e);
133 double StiffenedGas::getDiffInternalEnergyTemperature(){
136 // function to compute partial rho / partial h (p constant), partial rho/ partial p (h constant)
137 // use to compute p_x stationary
138 double StiffenedGas::getDiffDensEnthalpyPressconstant(const double p, const double h){
139 double p_inf = - _gamma*_p0;
140 return (p_inf - _gamma*p)/((_gamma-1)*h*h);
142 double StiffenedGas::getDiffDensPressEnthalpyconstant(const double h){
143 return _gamma/((_gamma-1)*h);