${CDMATH_INCLUDES} #
${PETSC_INCLUDES} #
${CoreFlows_MODELS}/inc #
+ ${CoreFlows_MODELS}/inc/EOS #
) #
add_subdirectory (${CoreFlows_MODELS})
DiffusionEquation.hxx
StationaryDiffusionEquation.hxx
LinearElasticityModel.hxx
- Fluide.h
ProblemFluid.hxx
utilitaire_algebre.h
+ EOS/Fluide.h
+ EOS/StiffenedGas.hxx
+# EOS/IAPWS97.h
)
INSTALL(FILES ${include_models_HXX} DESTINATION include)
--- /dev/null
+#ifndef FLUIDE_H
+#define FLUIDE_H
+
+#include <string>
+#include "math.h"
+#include "CdmathException.hxx"
+
+/*! \class Fluide Fluide.hxx "Fluide.hxx"
+ * \brief Fluid thermodynamics properties
+ * \details Provides pressure, density, temperature, internal energy, enthalpy, viscosity and conductivity
+ */
+
+using namespace std;
+
+class Fluide{
+ protected:
+ double _mu, _lambda,_Cv, _Cp,_Tref,_gamma, _p0, _q,_dragCoeff;
+ public:
+ Fluide(){_mu=0; _lambda=0;_Cv=0;_Cp=0;_Tref=0;_dragCoeff=0;_gamma=0; _p0=0; _q=0;}
+ virtual ~Fluide(){};
+
+ //Stiffened gas equation of state
+ double getPressure(double rhoe,const double rho) {
+ return (_gamma - 1) * (rhoe - rho*_q) - _gamma*_p0;
+ };
+ double getPressureFromEnthalpy(double h,const double rho) {
+ return (_gamma - 1)/_gamma * rho * (h - _q) - _p0;
+ };
+ /*For the newton scheme in the IsothermalTwoFluid model */
+ double getPressureDerivativeRhoE() { return _gamma - 1; }
+ double getDensityFromEnthalpy(double p, double h)
+ {
+ return _gamma*(p+_p0)/((_gamma-1)*(h-_q));
+ }
+ double vitesseSonEnthalpie(double h) { return sqrt((_gamma-1)*h); };
+ double vitesseSonTemperature(const double T, const double rho)
+ {
+ double h= getEnthalpy(T,rho);
+ return vitesseSonEnthalpie(h);
+ }
+
+ double getViscosity(double T) {return _mu;};
+ double getConductivity(double T) {return _lambda;};
+ double getDragCoeffs(double T) { return _dragCoeff;};
+ void setViscosity(double mu) { _mu=mu;};
+ void setConductivity(double lambda) { _lambda= lambda;};
+ void setDragCoeffs(double dragCoeff) {_dragCoeff=dragCoeff;};
+ //return constants gamma, cp, cv, p0, q
+ double constante(string name)
+ {
+ if(name== "gamma")
+ return _gamma;
+ else if (name == "cv"||name == "Cv")
+ return _Cv;
+ else if (name == "cp"||name == "Cp")
+ return _Cp;
+ else if(name=="p0")
+ return _p0;
+ else if(name=="q")
+ return _q;
+ else
+ throw CdmathException("Unknown constant: "+name);
+ }
+ virtual double getDensity(double p, double T)=0;
+ virtual double getTemperatureFromPressure(const double p, const double rho)=0;
+ virtual double getTemperatureFromEnthalpy(const double h, const double rho)=0;
+ virtual double getInternalEnergy(double T, double rho)=0;
+ virtual double getEnthalpy(double T, double rho)=0;
+
+};
+
+#endif
--- /dev/null
+#ifndef IAPWS97_H
+#define IAPWS97_H
+
+#include <string>
+#include "CdmathException.hxx"
+#include "Fluide.h"
+
+/*! \class IAPWS97 IAPWS97.hxx "IAPWS97.hxx"
+ * \brief IAPWS97 release of water and steam thermodynamics properties (freesteam library)
+ * \details Provides pressure, density, temperature, internal energy, enthalpy, viscosity and conductivity
+ */
+
+using namespace std;
+
+//The IAPWS-IF97 law implemented by the freesteam team
+class FluideIAPWS97:public Fluide{
+ protected:
+
+ public:
+ FluideIAPWS97(){}
+ virtual ~FluideIAPWS97(){};
+
+ //Stiffened gas equation of state
+ double getPressure(double rhoe,const double rho) {
+ return (_gamma - 1) * (rhoe - rho*_q) - _gamma*_p0;
+ };
+ double getPressureFromEnthalpy(double h,const double rho) {
+ return (_gamma - 1)/_gamma * rho * (h - _q) - _p0;
+ };
+ //For the newton scheme in the IsothermalTwoFluid model
+ double getPressureDerivativeRhoE() { return _gamma - 1; }
+ double getDensityFromEnthalpy(double p, double h) { return rhomass_phmass(p,h); }
+ double vitesseSonEnthalpie(double h, double p) { return speed_sound(T_phmass(p,h), p); };
+ double vitesseSonTemperature(const double T, const double rho) { return speed_sound(T,rho); }
+
+ double getViscosity(double T, rho) {return visc(T, rho);};
+ double getConductivity(double T, double p) {return tcond(T,p,getDensity(T,p);};
+ //return constants gamma, cp, cv, p0, q
+ double constante(string name, double T, double p)
+ {
+ if(name== "gamma")
+ return cpmass(T,p)/cvmass(T,p);
+ else if (name == "cv"||name == "Cv")
+ return cvmass(T,p);
+ else if (name == "cp"||name == "Cp")
+ return cpmass(T,p);
+ else
+ throw CdmathException("Unknown constant: "+name);
+ }
+ double getDensity(double p, double T) { return rhomass(T,p); }
+ double getTemperatureFromPressure(const double p, const double rho) { return T_phmass(p,h); }
+ double getTemperatureFromEnthalpy(const double h, const double rho) { return T_phmass(p,h); }
+ double getInternalEnergy(double T, double rho) { return umass(T,rho); }
+ double getEnthalpy(double T, double rho) { return hmass(T,rho); }
+
+};
+
+#endif
--- /dev/null
+#ifndef STIFFENEDGAS_H
+#define STIFFENEDGAS_H
+
+#include <string>
+#include "math.h"
+#include "CdmathException.hxx"
+#include "Fluide.h"
+
+/*! \class StiffenedGas StiffenedGas.hxx "StiffenedGas.hxx"
+ * \brief Stiffened Gas law approximating water and steam : P = (gamma - 1) * (rho * e - rho * q) - gamma * p0
+ * \details Provides pressure, density, temperature, internal energy, enthalpy, viscosity and conductivity
+ */
+
+using namespace std;
+
+// A standard stiffened gas class
+
+/*! \class StiffenedGas Fluide.hxx "Fluide.hxx"
+ * \brief Class implementing a standard stiffened gas law between pressure, density and internal energy
+ * \details
+ */
+class StiffenedGas:public Fluide{
+ private:
+ double _e_ref;//Stiffened gas law : P=(gamma - 1) * rho e(T) - _gamma*_p0
+ public:
+ StiffenedGas():Fluide(){_e_ref=0;};
+ //Perfect gas EOS
+ StiffenedGas( double gamma, double cv, double T_ref, double e_ref);
+ //Stiffened gas law fitting reference pressure, density and sound speed
+ StiffenedGas(double rho_ref, double p_ref, double T_ref, double e_ref, double soundSpeed_ref, double heatCap_ref);
+
+ double getInternalEnergy(double T, double rho=0);
+ double getEnthalpy(double T, double rho);
+ double getTemperatureFromPressure(double p, double rho);
+ double getTemperatureFromEnthalpy(const double h, const double rho);
+ double getDensity(double p, double T);
+
+ // Functions used to compute the Roe matrix for the five equation model (Kieu)
+ /* get differential of the density rho = rho(P,e)
+ * wrt the pressure (const e) wrt the internal energy (const P) */
+ double getJumpDensPress(const double e_l, const double e_r);
+ double getJumpDensInternalEnergy(const double p_l,const double p_r,const double e_l,const double e_r);
+ double getJumpInternalEnergyTemperature();
+ double getDiffDensPress(const double e);
+ double getDiffDensInternalEnergy(const double p,const double e);
+ double getDiffInternalEnergyTemperature();
+ /* get differential of the density rho = rho(P,h)
+ * wrt the pressure (const h) wrt the enthalpy (const P) */
+ double getDiffDensEnthalpyPressconstant(const double p, const double h);
+ double getDiffDensPressEnthalpyconstant(const double h);
+};
+
+// S. Dellacherie stiffened gas class
+
+/*! \class StiffenedGasDellacherie Fluide.hxx "Fluide.hxx"
+ * \brief Class implementing a particular stiffened gas law including saturation properties
+ * \details
+ */
+class StiffenedGasDellacherie:public Fluide{
+ private:
+ double _h_ref;//Stiffened gas law according to S. Dellacherie : P=(gamma - 1) * rho (e(T)-q) - _gamma*_p0
+ public:
+ StiffenedGasDellacherie():Fluide(){_h_ref=0;};
+ /* Loi des gaz raidis avec coefficients imposés suivant S. Dellacherie*/
+ StiffenedGasDellacherie( double gamma, double p0, double q, double cv);
+
+ double getInternalEnergy(double T, double rho);
+ double getEnthalpy(double T, double rho=0);
+ double getTemperatureFromPressure(double p, double rho);
+ double getTemperatureFromEnthalpy(const double h, const double rho=0);
+ double getDensity(double p, double T);
+
+ };
+
+#endif
+++ /dev/null
-#ifndef FLUIDE_H
-#define FLUIDE_H
-
-#include <string>
-#include "math.h"
-#include "CdmathException.hxx"
-
-/*! \class Fluide Fluide.hxx "Fluide.hxx"
- * \brief Fluid thermodynamics properties
- * \details Provides pressure, density, temperature, internal energy, enthalpy, viscosity and conductivity
- */
-
-using namespace std;
-
-class Fluide{
- protected:
- double _mu, _lambda,_Cv, _Cp,_Tref,_gamma, _p0, _q,_dragCoeff;
- public:
- Fluide(){_mu=0; _lambda=0;_Cv=0;_Cp=0;_Tref=0;_dragCoeff=0;_gamma=0; _p0=0; _q=0;}
- virtual ~Fluide(){};
-
- //Stiffened gas equation of state
- double getPressure(double rhoe,const double rho) {
- return (_gamma - 1) * (rhoe - rho*_q) - _gamma*_p0;
- };
- double getPressureFromEnthalpy(double h,const double rho) {
- return (_gamma - 1)/_gamma * rho * (h - _q) - _p0;
- };
- /*For the newton scheme in the IsothermalTwoFluid model */
- double getPressureDerivativeRhoE() { return _gamma - 1; }
- double getDensityFromEnthalpy(double p, double h)
- {
- return _gamma*(p+_p0)/((_gamma-1)*(h-_q));
- }
- double vitesseSonEnthalpie(double h) { return sqrt((_gamma-1)*h); };
- double vitesseSonTemperature(const double T, const double rho)
- {
- double h= getEnthalpy(T,rho);
- return vitesseSonEnthalpie(h);
- }
-
- double getViscosity(double T) {return _mu;};
- double getConductivity(double T) {return _lambda;};
- double getDragCoeffs(double T) { return _dragCoeff;};
- void setViscosity(double mu) { _mu=mu;};
- void setConductivity(double lambda) { _lambda= lambda;};
- void setDragCoeffs(double dragCoeff) {_dragCoeff=dragCoeff;};
- //return constants gamma, cp, cv, p0, q
- double constante(string name)
- {
- if(name== "gamma")
- return _gamma;
- else if (name == "cv"||name == "Cv")
- return _Cv;
- else if (name == "cp"||name == "Cp")
- return _Cp;
- else if(name=="p0")
- return _p0;
- else if(name=="q")
- return _q;
- else
- throw CdmathException("Unknown constant: "+name);
- }
- virtual double getDensity(double p, double T)=0;
- virtual double getTemperatureFromPressure(const double p, const double rho)=0;
- virtual double getTemperatureFromEnthalpy(const double h, const double rho)=0;
- virtual double getInternalEnergy(double T, double rho)=0;
- virtual double getEnthalpy(double T, double rho)=0;
-
-};
-
-// A standard stiffened gas class
-
-/*! \class StiffenedGas Fluide.hxx "Fluide.hxx"
- * \brief Class implementing a standard stiffened gas law between pressure, density and internal energy
- * \details
- */
-class StiffenedGas:public Fluide{
- private:
- double _e_ref;//Stiffened gas law : P=(gamma - 1) * rho e(T) - _gamma*_p0
- public:
- StiffenedGas():Fluide(){_e_ref=0;};
- //Perfect gas EOS
- StiffenedGas( double gamma, double cv, double T_ref, double e_ref);
- //Stiffened gas law fitting reference pressure, density and sound speed
- StiffenedGas(double rho_ref, double p_ref, double T_ref, double e_ref, double soundSpeed_ref, double heatCap_ref);
-
- double getInternalEnergy(double T, double rho=0);
- double getEnthalpy(double T, double rho);
- double getTemperatureFromPressure(double p, double rho);
- double getTemperatureFromEnthalpy(const double h, const double rho);
- double getDensity(double p, double T);
-
- // Functions used to compute the Roe matrix for the five equation model (Kieu)
- /* get differential of the density rho = rho(P,e)
- * wrt the pressure (const e) wrt the internal energy (const P) */
- double getJumpDensPress(const double e_l, const double e_r);
- double getJumpDensInternalEnergy(const double p_l,const double p_r,const double e_l,const double e_r);
- double getJumpInternalEnergyTemperature();
- double getDiffDensPress(const double e);
- double getDiffDensInternalEnergy(const double p,const double e);
- double getDiffInternalEnergyTemperature();
- /* get differential of the density rho = rho(P,h)
- * wrt the pressure (const h) wrt the enthalpy (const P) */
- double getDiffDensEnthalpyPressconstant(const double p, const double h);
- double getDiffDensPressEnthalpyconstant(const double h);
-};
-
-// S. Dellacherie stiffened gas class
-
-/*! \class StiffenedGasDellacherie Fluide.hxx "Fluide.hxx"
- * \brief Class implementing a particular stiffened gas law including saturation properties
- * \details
- */
-class StiffenedGasDellacherie:public Fluide{
- private:
- double _h_ref;//Stiffened gas law according to S. Dellacherie : P=(gamma - 1) * rho (e(T)-q) - _gamma*_p0
- public:
- StiffenedGasDellacherie():Fluide(){_h_ref=0;};
- /* Loi des gaz raidis avec coefficients imposés suivant S. Dellacherie*/
- StiffenedGasDellacherie( double gamma, double p0, double q, double cv);
-
- double getInternalEnergy(double T, double rho);
- double getEnthalpy(double T, double rho=0);
- double getTemperatureFromPressure(double p, double rho);
- double getTemperatureFromEnthalpy(const double h, const double rho=0);
- double getDensity(double p, double T);
-
- };
-#endif
#define SINGLEPHASE_HXX_
#include "ProblemFluid.hxx"
+#include "StiffenedGas.hxx"
class SinglePhase : public ProblemFluid{
public :
* */
bool iterateTimeStep(bool &ok);
+ //EOS functions
+ StiffenedGas getFluidEOS(){ return *dynamic_cast<StiffenedGas*>(_fluides[0]); }
double getReferencePressure() { return _Pref; };
double getReferenceTemperature() { return _Tref; };
DiffusionEquation.cxx
StationaryDiffusionEquation.cxx
# LinearElasticityModel.cxx
- Fluide.cxx
ProblemFluid.cxx
utilitaire_algebre.cxx
+ EOS/StiffenedGas.cxx
+# EOS/IAPWS97.cxx
)
ADD_LIBRARY(CoreFlowsLibs SHARED ${src_models_CXX})
* Author: Michael Ndjinga
*/
#include "DriftModel.hxx"
+#include "StiffenedGas.hxx"
using namespace std;
--- /dev/null
+#include "StiffenedGas.hxx"
+#include <iostream>
+
+//Perfect gas EOS Loi d'etat gaz parfait
+StiffenedGas::StiffenedGas( double gamma, double cv, double T_ref, double e_ref): Fluide()
+{
+ if(gamma -1<=0)
+ throw CdmathException("StiffenedGas::StiffenedGas: gamma<1");
+ _gamma=gamma;
+ _Cv=cv;
+ _Cp=_gamma*_Cv;
+ _p0=0;
+ _Tref=T_ref;
+ _e_ref=e_ref;
+ _q=0;
+ cout<<"Perfect gas EOS P=(gamma - 1) * rho e with parameter"<< " gamma= " << _gamma<<endl;
+ 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;
+}
+//Stiffened gas fitted using sound speed
+StiffenedGas::StiffenedGas(double rho_ref, double p_ref, double T_ref, double e_ref, double c_ref, double cv_ref): Fluide()
+{
+ //Old formula
+ //_gamma=(1+sqrt(1+4*c_ref*c_ref/(cv_ref*T_ref)))/2;
+ //New formula
+ _e_ref=e_ref;
+ _gamma=1+c_ref*c_ref/(_e_ref+p_ref/rho_ref);
+ if(_gamma -1<=0)
+ throw CdmathException("StiffenedGas::setEOS: gamma<1");
+ _Tref=T_ref;
+ _Cv=cv_ref;
+ _Cp=_gamma*_Cv;
+ _p0=((_gamma-1)*rho_ref*_e_ref-p_ref)/_gamma;
+ _q=0;
+ cout<<"Stiffened gas EOS P=(gamma - 1) * rho (e(T)-q) - gamma*p0 with parameters"<< " p0= " << _p0 << " gamma= " << _gamma<<endl;
+ cout<<"Calibrated around pressure= "<<p_ref<< " Pa, density= "<< rho_ref<< " internal energy= " << e_ref << " J/Kg, sound speed= " << c_ref<< endl;
+ 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;
+}
+// Loi d'etat stiffened gas S. Dellacherie
+StiffenedGasDellacherie::StiffenedGasDellacherie( double gamma, double p0, double q, double cv_ref)
+{
+ if(gamma -1<=0)
+ throw CdmathException("StiffenedGas::StiffenedGas: gamma<1");
+ _gamma=gamma;
+ _Cv=cv_ref;
+ _Cp=_gamma*_Cv;
+ _p0=p0;
+ _q=q;
+ _Tref=0;
+ _h_ref=q;
+
+ cout<<"S. Dellacherie Stiffened gas EOS P=(gamma - 1)/gamma * rho (h(T)-q) - p0 with parameters"<< " p0=" << _p0 << " gamma= " << _gamma<< " q= " << _q<< endl;
+ 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;
+
+}
+double StiffenedGas::getInternalEnergy(double T, double rho)
+{
+ return _Cv*(T-_Tref)+_e_ref;
+}
+double StiffenedGasDellacherie::getInternalEnergy(double T, double rho)
+{
+ double h= getEnthalpy(T);//h=e+p/rho=e+(gamma-1)(e-q)-gamma p0/rho=gamma(e- p0/rho)-(gamma-1)q
+ return (h+(_gamma-1)*_q)/_gamma+_p0/rho;
+}
+double StiffenedGas::getEnthalpy(double T, double rho)
+{
+ double e=getInternalEnergy( T, rho);
+ return _gamma*(e-_p0/rho)-(_gamma-1)*_q;
+}
+double StiffenedGasDellacherie::getEnthalpy(double T, double rho)
+{
+ return _Cp*(T-_Tref)+_h_ref;
+}
+double StiffenedGas::getTemperatureFromPressure(const double p, const double rho)
+{
+ //p=(gamma-1)rho(e-q)-gamma p0
+ double e=_q +(p+_gamma*_p0)/((_gamma-1)*rho);
+ return (e-_e_ref)/_Cv + _Tref;
+}
+double StiffenedGasDellacherie::getTemperatureFromPressure(const double p, const double rho)
+{
+ //P=(gamma - 1)/gamma * rho (h(T)-q) - p0
+ double h=_q +_gamma*(p+_p0)/((_gamma-1)*rho);
+ return (h-_h_ref)/_Cp + _Tref;
+}
+
+double StiffenedGas::getTemperatureFromEnthalpy(const double h, const double rho)
+{
+ //h=e+p/rho et p=(gamma-1)rho(e-q)-gamma p0
+ double e=(h+(_gamma-1)*_q)/_gamma + _p0/rho;
+ return (e-_e_ref)/_Cv + _Tref;
+}
+double StiffenedGasDellacherie::getTemperatureFromEnthalpy(const double h, const double rho)
+{
+ return (h-_h_ref)/_Cp + _Tref;
+}
+
+double StiffenedGas::getDensity(double p, double T)
+{
+ return (p+_gamma*_p0)/((_gamma-1)*(getInternalEnergy(T)-_q));
+}
+double StiffenedGasDellacherie::getDensity(double p, double T)
+{
+ return _gamma*(p+_p0)/((_gamma-1)*(getEnthalpy(T)-_q));
+}
+
+double StiffenedGas::getJumpDensPress(const double e_l, const double e_r){
+ double inv_a_2;
+ inv_a_2 = (e_l+e_r)/((_gamma-1)*2*e_l*e_r);
+ return inv_a_2;
+}
+double StiffenedGas::getJumpDensInternalEnergy(const double p_l,const double p_r,const double e_l,const double e_r){
+ double b, p_inf;
+ p_inf = - _gamma*_p0;
+ b = (p_inf-0.5*(p_l+p_r))/((_gamma-1)*e_l*e_r);
+ return b;
+}
+double StiffenedGas::getJumpInternalEnergyTemperature(){
+ return _Cv;
+}
+// function to compute partial rho/ partial p (e constant), partial rho/partial e (p constant)
+// use to compute the Jacobian matrix
+double StiffenedGas::getDiffDensPress(const double e){
+ double inv_a_2;
+ inv_a_2 = 1/((_gamma-1)*e);
+ return inv_a_2;
+}
+double StiffenedGas::getDiffDensInternalEnergy(const double p,const double e){
+ double b, p_inf;
+ p_inf = - _gamma*_p0;
+ b = (p_inf-p)/((_gamma-1)*e*e);
+ return b;
+}
+double StiffenedGas::getDiffInternalEnergyTemperature(){
+ return _Cv;
+}
+// function to compute partial rho / partial h (p constant), partial rho/ partial p (h constant)
+// use to compute p_x stationary
+double StiffenedGas::getDiffDensEnthalpyPressconstant(const double p, const double h){
+ double p_inf = - _gamma*_p0;
+ return (p_inf - _gamma*p)/((_gamma-1)*h*h);
+}
+double StiffenedGas::getDiffDensPressEnthalpyconstant(const double h){
+ return _gamma/((_gamma-1)*h);
+}
*/
#include "FiveEqsTwoFluid.hxx"
#include <cstdlib>
+#include "StiffenedGas.hxx"
using namespace std;
+++ /dev/null
-#include "Fluide.h"
-#include <iostream>
-
-//Perfect gas EOS Loi d'etat gaz parfait
-StiffenedGas::StiffenedGas( double gamma, double cv, double T_ref, double e_ref): Fluide()
-{
- if(gamma -1<=0)
- throw CdmathException("StiffenedGas::StiffenedGas: gamma<1");
- _gamma=gamma;
- _Cv=cv;
- _Cp=_gamma*_Cv;
- _p0=0;
- _Tref=T_ref;
- _e_ref=e_ref;
- _q=0;
- cout<<"Perfect gas EOS P=(gamma - 1) * rho e with parameter"<< " gamma= " << _gamma<<endl;
- 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;
-}
-//Stiffened gas fitted using sound speed
-StiffenedGas::StiffenedGas(double rho_ref, double p_ref, double T_ref, double e_ref, double c_ref, double cv_ref): Fluide()
-{
- //Old formula
- //_gamma=(1+sqrt(1+4*c_ref*c_ref/(cv_ref*T_ref)))/2;
- //New formula
- _e_ref=e_ref;
- _gamma=1+c_ref*c_ref/(_e_ref+p_ref/rho_ref);
- if(_gamma -1<=0)
- throw CdmathException("StiffenedGas::setEOS: gamma<1");
- _Tref=T_ref;
- _Cv=cv_ref;
- _Cp=_gamma*_Cv;
- _p0=((_gamma-1)*rho_ref*_e_ref-p_ref)/_gamma;
- _q=0;
- cout<<"Stiffened gas EOS P=(gamma - 1) * rho (e(T)-q) - gamma*p0 with parameters"<< " p0= " << _p0 << " gamma= " << _gamma<<endl;
- cout<<"Calibrated around pressure= "<<p_ref<< " Pa, density= "<< rho_ref<< " internal energy= " << e_ref << " J/Kg, sound speed= " << c_ref<< endl;
- 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;
-}
-// Loi d'etat stiffened gas S. Dellacherie
-StiffenedGasDellacherie::StiffenedGasDellacherie( double gamma, double p0, double q, double cv_ref)
-{
- if(gamma -1<=0)
- throw CdmathException("StiffenedGas::StiffenedGas: gamma<1");
- _gamma=gamma;
- _Cv=cv_ref;
- _Cp=_gamma*_Cv;
- _p0=p0;
- _q=q;
- _Tref=0;
- _h_ref=q;
-
- cout<<"S. Dellacherie Stiffened gas EOS P=(gamma - 1)/gamma * rho (h(T)-q) - p0 with parameters"<< " p0=" << _p0 << " gamma= " << _gamma<< " q= " << _q<< endl;
- 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;
-
-}
-double StiffenedGas::getInternalEnergy(double T, double rho)
-{
- return _Cv*(T-_Tref)+_e_ref;
-}
-double StiffenedGasDellacherie::getInternalEnergy(double T, double rho)
-{
- double h= getEnthalpy(T);//h=e+p/rho=e+(gamma-1)(e-q)-gamma p0/rho=gamma(e- p0/rho)-(gamma-1)q
- return (h+(_gamma-1)*_q)/_gamma+_p0/rho;
-}
-double StiffenedGas::getEnthalpy(double T, double rho)
-{
- double e=getInternalEnergy( T, rho);
- return _gamma*(e-_p0/rho)-(_gamma-1)*_q;
-}
-double StiffenedGasDellacherie::getEnthalpy(double T, double rho)
-{
- return _Cp*(T-_Tref)+_h_ref;
-}
-double StiffenedGas::getTemperatureFromPressure(const double p, const double rho)
-{
- //p=(gamma-1)rho(e-q)-gamma p0
- double e=_q +(p+_gamma*_p0)/((_gamma-1)*rho);
- return (e-_e_ref)/_Cv + _Tref;
-}
-double StiffenedGasDellacherie::getTemperatureFromPressure(const double p, const double rho)
-{
- //P=(gamma - 1)/gamma * rho (h(T)-q) - p0
- double h=_q +_gamma*(p+_p0)/((_gamma-1)*rho);
- return (h-_h_ref)/_Cp + _Tref;
-}
-
-double StiffenedGas::getTemperatureFromEnthalpy(const double h, const double rho)
-{
- //h=e+p/rho et p=(gamma-1)rho(e-q)-gamma p0
- double e=(h+(_gamma-1)*_q)/_gamma + _p0/rho;
- return (e-_e_ref)/_Cv + _Tref;
-}
-double StiffenedGasDellacherie::getTemperatureFromEnthalpy(const double h, const double rho)
-{
- return (h-_h_ref)/_Cp + _Tref;
-}
-
-double StiffenedGas::getDensity(double p, double T)
-{
- return (p+_gamma*_p0)/((_gamma-1)*(getInternalEnergy(T)-_q));
-}
-double StiffenedGasDellacherie::getDensity(double p, double T)
-{
- return _gamma*(p+_p0)/((_gamma-1)*(getEnthalpy(T)-_q));
-}
-
-double StiffenedGas::getJumpDensPress(const double e_l, const double e_r){
- double inv_a_2;
- inv_a_2 = (e_l+e_r)/((_gamma-1)*2*e_l*e_r);
- return inv_a_2;
-}
-double StiffenedGas::getJumpDensInternalEnergy(const double p_l,const double p_r,const double e_l,const double e_r){
- double b, p_inf;
- p_inf = - _gamma*_p0;
- b = (p_inf-0.5*(p_l+p_r))/((_gamma-1)*e_l*e_r);
- return b;
-}
-double StiffenedGas::getJumpInternalEnergyTemperature(){
- return _Cv;
-}
-// function to compute partial rho/ partial p (e constant), partial rho/partial e (p constant)
-// use to compute the Jacobian matrix
-double StiffenedGas::getDiffDensPress(const double e){
- double inv_a_2;
- inv_a_2 = 1/((_gamma-1)*e);
- return inv_a_2;
-}
-double StiffenedGas::getDiffDensInternalEnergy(const double p,const double e){
- double b, p_inf;
- p_inf = - _gamma*_p0;
- b = (p_inf-p)/((_gamma-1)*e*e);
- return b;
-}
-double StiffenedGas::getDiffInternalEnergyTemperature(){
- return _Cv;
-}
-// function to compute partial rho / partial h (p constant), partial rho/ partial p (h constant)
-// use to compute p_x stationary
-double StiffenedGas::getDiffDensEnthalpyPressconstant(const double p, const double h){
- double p_inf = - _gamma*_p0;
- return (p_inf - _gamma*p)/((_gamma-1)*h*h);
-}
-double StiffenedGas::getDiffDensPressEnthalpyconstant(const double h){
- return _gamma/((_gamma-1)*h);
-}
*/
#include "IsothermalTwoFluid.hxx"
+#include "StiffenedGas.hxx"
using namespace std;
*/
#include "SinglePhaseStaggered.hxx"
+#include "StiffenedGas.hxx"
using namespace std;
+++ /dev/null
-//////////////////////////////////////////////////////////////////////////////\r
-//\r
-// File: fluid.h\r
-// Directory: \r
-// Version: \r
-//\r
-//////////////////////////////////////////////////////////////////////////////\r
-\r
-#ifndef fluid_included\r
-#define fluid_included\r
-\r
-#include <math.h>\r
-#include <Noms.h>\r
-#include <SFichier.h>\r
-#include <Process.h>\r
-\r
-class fluid\r
-{\r
-public :\r
- fluid() { fluide = 0; };\r
- void set_fluid(Nom &nom)\r
- {\r
- Noms liste(2); liste(0) = "Na"; liste(1) = "LBe"; \r
- fluide = liste.rang(nom);\r
- if (fluide == -1)\r
- {\r
- if (!Process::me()) Cerr << "fluide " << nom << " non reconnu parmi : " << liste << finl;\r
- Process::exit();\r
- }\r
- }\r
-\r
- //proprietes physiques\r
- /* Lois physiques du sodium issues de fits sur DEN/DANS/DM2S/STMF/LMES/RT/12-018/A */\r
- int fluide;\r
- #define Tk ((T) + 273.15)\r
- #define Ksil (fluide ? 3.022e-11 : 1e-9)\r
-\r
- /* inverse de la densite du liquide (hors pression) */\r
- inline double IRhoL(double T)\r
- {\r
- return fluide ? (9.03e-5 + Tk * (1.003e-8 + Tk * 2.01e-12))\r
- : (9.829728e-4 + Tk * (2.641186e-7 + Tk * (-3.340743e-11 + Tk * 6.680973e-14)));\r
- }\r
-\r
- /* sa derivee */\r
- inline double DTIRhoL(double T)\r
- {\r
- return fluide ? (1.003e-8 + 2.01e-12 * 2)\r
- : (2.641186e-7 + Tk * (-3.340743e-11 * 2 + Tk * 6.680973e-14 * 3));\r
- }\r
-\r
- /* densite du liquide */\r
- inline double RhoL(double T, double P)\r
- {\r
- return exp(Ksil * (P - 1e5)) / IRhoL(T);\r
- }\r
-\r
- /* enthalpie du liquide */\r
- inline double HL(double T, double P)\r
- {\r
- return (fluide ? 97.98e3 + Tk * (159 + Tk * (-2.72e-2 / 2 + Tk * 7.12e-6 / 3))\r
- : 2992600 / Tk - 365487.2 + Tk * (1658.2 + Tk * (-.42395 + Tk * 1.4847e-4)))\r
- + (IRhoL(T) - Tk * DTIRhoL(T)) * (1 - exp(Ksil*(1e5 - P))) / Ksil;\r
- }\r
- \r
- /* derivees */\r
- inline double DTHL(double T, double P)\r
- {\r
- return (fluide ? 159 + Tk * (-2.72e-2 + Tk * 7.12e-6)\r
- : -2992600 / (Tk*Tk) + 1658.2 + Tk * (-.42395 * 2 + Tk * 1.4847e-4 * 3))\r
- - Tk * (fluide ? 0 : -3.340743e-11 * 2 + Tk * 6.680973e-14 * 6) * (1 - exp(Ksil*(1e5 - P))) / Ksil;\r
- }\r
- \r
- /* inverses par methode de Newton */\r
- inline double TLh(double h, double T0, double P)\r
- {\r
- double T = T0, sec; int i;\r
- for (i = 0; i < 100 && dabs( sec = HL(T, P) - h ) > 1e-8; i++)\r
- T -= sec / DTHL(T, P);\r
- return i < 100 ? T : -1e10;\r
- }\r
- /* derivees */\r
- inline double DTRhoL(double T, double P)\r
- {\r
- return - exp(Ksil * (P - 1e5)) * DTIRhoL(T) / (IRhoL(T) * IRhoL(T));\r
- }\r
- inline double DPRhoL(double T, double P)\r
- {\r
- return Ksil * RhoL(T, P);\r
- }\r
- \r
- /* derivees */\r
- inline double DPHL(double T, double P)\r
- {\r
- return (IRhoL(T) - Tk * DTIRhoL(T)) * exp(Ksil * (1e5 - P));\r
- }\r
- \r
- /* energie volumique du liquide */\r
- inline double EL(double T, double P)\r
- {\r
- return HL(T, P) * RhoL(T, P);\r
- }\r
- /* derivees */\r
- inline double DTEL(double T, double P)\r
- {\r
- return DTHL(T, P) * RhoL(T, P) + HL(T, P) * DTRhoL(T, P);\r
- }\r
- inline double DPEL(double T, double P)\r
- {\r
- return DPHL(T, P) * RhoL(T, P) + HL(T, P) * DPRhoL(T, P);\r
- }\r
- \r
- inline double TLe(double e, double T0, double P)\r
- {\r
- double T = T0, sec; int i;\r
- for (i = 0; i < 100 && dabs( sec = EL(T, P) - e ) > 1e-3; i++)\r
- T -= sec / DTEL(T, P);\r
- return i < 100 ? T : -1e10;\r
- }\r
- \r
- /* viscosite du liquide*/\r
- inline double MuL( double T )\r
- {\r
- return fluide ? 4.94e-4 * exp(754.1 / Tk)\r
- : exp( -6.4406 - 0.3958 * log(Tk) + 556.835 / Tk );\r
- }\r
- /* conductivite du liquide */\r
- inline double LambdaL( double T )\r
- {\r
- return fluide ? 3.61 + Tk * (1.517e-2 - Tk * 1.741e-6)\r
- : max(124.67 + Tk * (-.11381 + Tk * (5.5226e-5 - Tk * 1.1848e-8)), 0.045);\r
- }\r
- /* tension superficielle */\r
- inline double SigmaL(double T)\r
- {\r
- return 0.2405 * pow(1 - Tk / 2503.7, 1.126);\r
- }\r
- /* Nusselt du liquide -> Skupinski */\r
- inline double NuL(double T, double P, double w, double Dh)\r
- {\r
- double Pe = RhoL(T, P) * dabs(w) * Dh * DTHL(T, P) / LambdaL(T);\r
- return 4.82 + 0.0185 * pow(Pe, 0.827);\r
- }\r
- /* temperature de saturation */\r
- inline double Tsat( double P )\r
- {\r
- double A = 7.8130, B = 11209, C = 5.249e5;\r
- return 2 * C / ( -B + sqrt(B*B + 4 * A * C - 4 * C * log(P / 1e6))) - 273.15;\r
- }\r
- /* sa derivee */\r
- inline double DTsat( double P )\r
- {\r
- double A = 7.8130, B = 11209, C = 5.249e5, Tsk = Tsat(P) + 273.15;\r
- return Tsk * Tsk / ( P * sqrt(B*B + 4 * A * C - 4 * C * log(P / 1e6)));\r
- }\r
- /* pression de vapeur saturante */\r
- inline double Psat( double T )\r
- {\r
- double A = 7.8130, B = 11209, C = 5.249e5;\r
- return 1e6 * exp(A - B / Tk - C / (Tk * Tk));\r
- }\r
- /* sa derivee */\r
- inline double DPsat( double T )\r
- {\r
- double B = 11209, C = 5.249e5;\r
- return (B / (Tk * Tk) + 2 * C / (Tk * Tk * Tk)) * Psat(T);\r
- }\r
- /* enthalpie massique de saturation */\r
- inline double Hsat( double P )\r
- {\r
- return HL(Tsat(P), P);\r
- }\r
- /* sa derivee */\r
- inline double DHsat( double P )\r
- {\r
- return DTsat(P) * DTHL(Tsat(P), P) + DPHL(Tsat(P), P);\r
- }\r
- /* chaleur latente, a prendre a Tsat */\r
- inline double Lvap( double P )\r
- {\r
- double Tc = 2503.7, Tsk = Tsat(P) + 273.15;\r
- return 3.9337e5 * ( 1 - Tsk / Tc) + 4.3986e6 * pow( 1 - Tsk / Tc, .29302);\r
- }\r
- \r
- /* sa derivee */\r
- inline double DLvap( double P )\r
- {\r
- double Tc = 2503.7, Tsk = Tsat(P) + 273.15;\r
- return DTsat(P) * (-3.9337e5 / Tc - 4.3986e6 * .29302 * pow( 1 - Tsk / Tc, .29302 - 1) / Tc);\r
- }\r
- \r
- /* densite de la vapeur : (gaz parfait) * f1(Tsat) * f2(DTsat)*/\r
- #define f1(Ts) (2.49121 + Ts * (-5.53796e-3 + Ts * (7.5465e-6 + Ts * (-4.20217e-9 + Ts * 8.59212e-13))))\r
- #define Df1(Ts) (-5.53796e-3 + Ts * (2 * 7.5465e-6 + Ts * (-3 * 4.20217e-9 + Ts * 4 * 8.59212e-13)))\r
- #define f2(dT) (1 + dT * (-5e-4 + dT * (6.25e-7 - dT * 4.1359e-25)))\r
- #define Df2(dT) (-5e-4 + dT * (2 * 6.25e-7 - dT * 3 * 4.1359e-25))\r
- inline double RhoV( double T, double P )\r
- {\r
- double Tsk = Tsat(P) + 273.15, dT = Tk - Tsk;\r
- return P * 2.7650313e-3 * f1(Tsk) * f2(dT) / Tk;\r
- }\r
- /* ses derivees */\r
- inline double DTRhoV( double T, double P )\r
- {\r
- double Tsk = Tsat(P) + 273.15, dT = Tk - Tsk;\r
- return P * 2.7650313e-3 * f1(Tsk) * (Df2(dT) - f2(dT) / Tk) / Tk;\r
- }\r
- inline double DPRhoV( double T, double P )\r
- {\r
- double Tsk = Tsat(P) + 273.15, dT = Tk - Tsk;\r
- return 2.7650313e-3 * (f1(Tsk) * f2(dT) + P * DTsat(P) * (Df1(Tsk) * f2(dT) - f1(Tsk) * Df2(dT))) / Tk;\r
- }\r
- #undef f1\r
- #undef Df1\r
- #undef f2\r
- #undef Df2\r
- /* inverse de la densite de la vapeur */\r
- inline double IRhoV( double T, double P )\r
- {\r
- return 1 / RhoV(T, P);\r
- }\r
- /* ses derivees */\r
- inline double DTIRhoV( double T, double P )\r
- {\r
- return -DTRhoV(T, P) / (RhoV(T, P) * RhoV(T, P));\r
- }\r
- inline double DPIRhoV( double T, double P )\r
- {\r
- return -DPRhoV(T, P) / (RhoV(T, P) * RhoV(T, P));\r
- }\r
- /* enthalpie de la vapeur */\r
- #define CpVs(Ts) (-1223.89 + Ts * (14.1191 + Ts * (-1.62025e-2 + Ts * 5.76923e-6)))\r
- #define DCpVs(Ts) (14.1191 + Ts * (- 2 * 1.62025e-2 + Ts * 3 * 5.76923e-6))\r
- inline double HV( double T, double P)\r
- {\r
- double Ts = Tsat(P), dT = T - Ts, Cp0 = 910, k = -4.6e-3;\r
- return HL(Ts, P) + Lvap(P) + Cp0 * dT + (Cp0 - CpVs(Ts)) * (1 - exp(k * dT)) / k;\r
- }\r
- /* ses derivees */\r
- inline double DTHV( double T, double P)\r
- {\r
- double Ts = Tsat(P), dT = T - Ts, Cp0 = 910, k = -4.6e-3;\r
- return Cp0 - (Cp0 - CpVs(Ts)) * exp(k * dT);\r
- }\r
- inline double DPHV( double T, double P)\r
- {\r
- double Ts = Tsat(P), dT = T - Ts, Cp0 = 910, k = -4.6e-3;\r
- return DPHL(Ts, P) + DLvap(P) + DTsat(P) * (DTHL(Ts, P) - Cp0 - DCpVs(Ts) * (1 - exp(k * dT)) / k + (Cp0 - CpVs(Ts)) * exp(k * dT));\r
- }\r
- #undef CpVs\r
- #undef DCpVs\r
- /* energie volumique de la vapeur */\r
- inline double EV( double T, double P)\r
- {\r
- return RhoV(T, P) * HV(T, P);\r
- }\r
- /* ses derivees */\r
- inline double DTEV(double T, double P)\r
- {\r
- return DTRhoV(T, P) * HV(T, P) + RhoV(T, P) * DTHV(T, P);\r
- }\r
- inline double DPEV(double T, double P)\r
- {\r
- return DPRhoV(T, P) * HV(T, P) + RhoV(T, P) * DPHV(T, P);\r
- }\r
- \r
- /* conductivite de la vapeur */\r
- inline double LambdaV( double T )\r
- {\r
- return 0.045;\r
- }\r
- \r
- /* viscosite de la vapeur */\r
- inline double MuV( double T )\r
- {\r
- return 2.5e-6 + 1.5e-8 * Tk;\r
- }\r
- \r
- /* Nusselt de la vapeur -> Dittus/ Boetler */\r
- inline double NuV(double T, double P, double w, double Dh)\r
- {\r
- double Re = RhoV(T, P) * dabs(w) * Dh / MuV(T), Pr = MuV(T) * DTHV(T, P) / LambdaV(T);\r
- return 0.023 * pow(Re, 0.8) * pow(Pr, 0.4);\r
- }\r
- \r
- inline double Rho( double T, double x, double P )\r
- {\r
- return 1 / ( (1-x)/RhoL(T, P) + x / RhoV(T, P) );\r
- }\r
- \r
- #undef Tk\r
- #undef Ksil\r
-};\r
-#endif\r
-#include "Fluide.h"
+#include "StiffenedGas.hxx"
#include <cstdlib>
#include <iostream>
#include "StationaryDiffusionEquation.hxx"
#include "SinglePhase.hxx"
#include "Fluide.h"
+#include "StiffenedGas.hxx"
%}
%include "StationaryDiffusionEquation.hxx"
%include "SinglePhase.hxx"
%include "Fluide.h"
+%include "StiffenedGas.hxx"
