CoreFlows/Models/src/Fluide.cxx

   1 #include "Fluide.h"
   2 #include <iostream>
   3
   4 //Perfect gas EOS Loi d'etat gaz parfait
   5 StiffenedGas::StiffenedGas( double gamma, double cv, double T_ref, double e_ref): Fluide()
   6 {
   7         if(gamma -1<=0)
   8                 throw CdmathException("StiffenedGas::StiffenedGas: gamma<1");
   9         _gamma=gamma;
  10         _Cv=cv;
  11         _Cp=_gamma*_Cv;
  12         _p0=0;
  13         _Tref=T_ref;
  14         _e_ref=e_ref;
  15         _q=0;
  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;
  18 }
  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()
  21 {
  22         //Old formula
  23         //_gamma=(1+sqrt(1+4*c_ref*c_ref/(cv_ref*T_ref)))/2;
  24         //New formula
  25         _e_ref=e_ref;
  26         _gamma=1+c_ref*c_ref/(_e_ref+p_ref/rho_ref);
  27         if(_gamma -1<=0)
  28                 throw CdmathException("StiffenedGas::setEOS: gamma<1");
  29         _Tref=T_ref;
  30         _Cv=cv_ref;
  31         _Cp=_gamma*_Cv;
  32         _p0=((_gamma-1)*rho_ref*_e_ref-p_ref)/_gamma;
  33         _q=0;
  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;
  37 }
  38 // Loi d'etat stiffened gas S. Dellacherie
  39 StiffenedGasDellacherie::StiffenedGasDellacherie( double gamma, double p0, double q, double cv_ref)
  40 {
  41         if(gamma -1<=0)
  42                 throw CdmathException("StiffenedGas::StiffenedGas: gamma<1");
  43         _gamma=gamma;
  44         _Cv=cv_ref;
  45         _Cp=_gamma*_Cv;
  46         _p0=p0;
  47         _q=q;
  48         _Tref=0;
  49         _h_ref=q;
  50
  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;
  53
  54 }
  55 double StiffenedGas::getInternalEnergy(double T, double rho)
  56 {
  57         return _Cv*(T-_Tref)+_e_ref;
  58 }
  59 double StiffenedGasDellacherie::getInternalEnergy(double T, double rho)
  60 {
  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;
  63 }
  64 double StiffenedGas::getEnthalpy(double T, double rho)
  65 {
  66         double e=getInternalEnergy( T, rho);
  67         return _gamma*(e-_p0/rho)-(_gamma-1)*_q;
  68 }
  69 double StiffenedGasDellacherie::getEnthalpy(double T, double rho)
  70 {
  71         return _Cp*(T-_Tref)+_h_ref;
  72 }
  73 double StiffenedGas::getTemperatureFromPressure(const double  p, const double rho)
  74 {
  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;
  78 }
  79 double StiffenedGasDellacherie::getTemperatureFromPressure(const double  p, const double rho)
  80 {
  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;
  84 }
  85
  86 double StiffenedGas::getTemperatureFromEnthalpy(const double  h, const double rho)
  87 {
  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;
  91 }
  92 double StiffenedGasDellacherie::getTemperatureFromEnthalpy(const double  h, const double rho)
  93 {
  94         return (h-_h_ref)/_Cp + _Tref;
  95 }
  96
  97 double StiffenedGas::getDensity(double p, double T)
  98 {
  99         return (p+_gamma*_p0)/((_gamma-1)*(getInternalEnergy(T)-_q));
 100 }
 101 double StiffenedGasDellacherie::getDensity(double p, double T)
 102 {
 103         return _gamma*(p+_p0)/((_gamma-1)*(getEnthalpy(T)-_q));
 104 }
 105
 106 double StiffenedGas::getJumpDensPress(const double e_l, const double e_r){
 107         double inv_a_2;
 108     inv_a_2 = (e_l+e_r)/((_gamma-1)*2*e_l*e_r);
 109         return inv_a_2;
 110 }
 111 double StiffenedGas::getJumpDensInternalEnergy(const double p_l,const double p_r,const double e_l,const double e_r){
 112         double b, p_inf;
 113         p_inf = - _gamma*_p0;
 114         b = (p_inf-0.5*(p_l+p_r))/((_gamma-1)*e_l*e_r);
 115         return b;
 116 }
 117 double StiffenedGas::getJumpInternalEnergyTemperature(){
 118         return _Cv;
 119 }
 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){
 123         double inv_a_2;
 124     inv_a_2 = 1/((_gamma-1)*e);
 125         return inv_a_2;
 126 }
 127 double StiffenedGas::getDiffDensInternalEnergy(const double p,const double e){
 128         double b, p_inf;
 129         p_inf = - _gamma*_p0;
 130         b = (p_inf-p)/((_gamma-1)*e*e);
 131         return b;
 132 }
 133 double StiffenedGas::getDiffInternalEnergyTemperature(){
 134         return _Cv;
 135 }
 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);
 141 }
 142 double StiffenedGas::getDiffDensPressEnthalpyconstant(const double h){
 143         return _gamma/((_gamma-1)*h);
 144 }