exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

real gas EOS for mixtures

Browse source
This commit is contained in:
Christian Lucas 2012-05-28 20:55:24 +02:00 committed by Henrik Rusche
parent e8b0b4968a
commit 47c93e3447
18 changed files with 649 additions and 793 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
src/thermophysicalModels/basic/basicThermo/basicThermo.H
View file

@ -102,28 +102,6 @@ protected:
//- Construct as copy (not implemented)
basicThermo(const basicThermo&);
//- Return the enthalpy field boundary types by interrogating the
// temperature field boundary types
// for real gases
wordList hRealBoundaryTypes();
//- Correct the enthalpy field boundaries
// for real gases
void hRealBoundaryCorrection(volScalarField& h);
// Internal energy
//- Return the internal energy field boundary types by
// interrogating the temperature field boundary types
// for real gases
wordList eRealBoundaryTypes();
//- Correct the internal energy field boundaries
// for real gases
void eRealBoundaryCorrection(volScalarField& e);
public:
//- Runtime type information

2
src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.C
View file

@ -157,8 +157,6 @@ Foam::realGasHThermo<MixtureType>::realGasHThermo(const fvMesh& mesh)
mesh,
dimensionSet(1, -5, 2, 0, 0)
)
{
scalarField& hCells = h_.internalField();

3
src/thermophysicalModels/specie/Make/files
View file

@ -16,5 +16,8 @@ $(equationOfState)/aungierRedlichKwong/aungierRedlichKwong.C
$(equationOfState)/pengRobinson/pengRobinson.C
$(equationOfState)/redlichKwong/redlichKwong.C
$(equationOfState)/soaveRedlichKwong/soaveRedlichKwong.C
$(equationOfState)/mixtureRedlichKwong/mixtureRedlichKwong.C
$(equationOfState)/mixturePengRobinson/mixturePengRobinson.C
$(equationOfState)/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwong.C
LIB = $(FOAM_LIBBIN)/libspecie

13
src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwong.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,8 +31,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "aungierRedlichKwong.H"
@ -43,7 +41,6 @@ Germany
namespace Foam
{
/* * * * * * * * * * * * * * * Private static data * * * * * * * * * * * * * */
const scalar aungierRedlichKwong::rhoMin_
@ -70,7 +67,7 @@ aungierRedlichKwong::aungierRedlichKwong(Istream& is)
c_(this->RR*Tcrit_/(pcrit_+(a_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_),
n_(0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())))
rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
{
is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
}
@ -78,10 +75,10 @@ aungierRedlichKwong::aungierRedlichKwong(Istream& is)
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const aungierRedlichKwong& pg)
Ostream& operator<<(Ostream& os, const aungierRedlichKwong& ark)
{
os << static_cast<const specie&>(pg)<< tab
<< pg.pcrit_ << tab<< pg.Tcrit_<< tab<<pg.azentricFactor_<< tab<<pg.rhocrit_;
os << static_cast<const specie&>(ark)<< tab
<< ark.pcrit_ << tab<< ark.Tcrit_<< tab<<ark.azentricFactor_<< tab<<ark.rhocrit_;
os.check("Ostream& operator<<(Ostream& os, const aungierRedlichKwong& st)");
return os;

64
src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwong.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,6 +28,11 @@ Class
Description
Aungier Redlich Kwong equation of state.
Paper:
Title: A Fast, Accurate Real Gas Equation of State for Fluid Dynamic Analysis Applications
Authors: R. H. Aungier
Journal: Journal of Fluids Engineering, Vol.117, 1995
SourceFiles
aungierRedlichKwongI.H
aungierRedlichKwong.C
@ -67,17 +72,18 @@ class aungierRedlichKwong
scalar rhocrit_;
//Aungier Redlich Kwong factors
scalar a_;
scalar b_;
scalar c_;
scalar n_;
mutable scalar a_;
mutable scalar b_;
mutable scalar c_;
mutable scalar n_;
//Density @STD, initialise after a, b!
scalar rhostd_;
// Static Private data
//Static Private data
//should be read from the fvSolution file where rhoMax and rhoMin values must be define ( for rhoSimpleFoam)
//HR: Don't know, yet. Let's make these static for starters
static const scalar rhoMax_;
static const scalar rhoMin_;
@ -113,9 +119,10 @@ public:
inline scalar rhostd() const;
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
//first order derivatives
//CL: first order derivatives
inline scalar dpdv(const scalar rho, const scalar T) const;
inline scalar dpdT(const scalar rho, const scalar T) const;
@ -132,14 +139,14 @@ public:
const scalar T
) const;
// Used for cv
//CL: Used for cv
inline scalar integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const ;
// second order derivatives, not Used At The Moment
//CL: second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho, const scalar T) const;
inline scalar d2pdT2(const scalar rho, scalar T) const;
@ -148,10 +155,10 @@ public:
inline scalar d2vdT2(const scalar rho, const scalar T) const;
//Used for internal Energy
//CL: Used for internal Energy
inline scalar integral_p_dv(const scalar rho, const scalar T) const;
//Used for Entropy
//CL: Used for Entropy
inline scalar integral_dpdT_dv(const scalar rho, const scalar T) const;
//- Return density [kg/m^3]
@ -172,45 +179,14 @@ public:
inline scalar Z
(
const scalar p,
const scalar T,
const scalar T,
const scalar rho0
) const;
// Member operators
inline void operator+=(const aungierRedlichKwong&);
inline void operator-=(const aungierRedlichKwong&);
inline void operator*=(const scalar);
/* // Friend operators
inline friend aungierRedlichKwong operator+
(
const aungierRedlichKwong&,
const aungierRedlichKwong&
);
inline friend aungierRedlichKwong operator-
(
const aungierRedlichKwong&,
const aungierRedlichKwong&
);
inline friend aungierRedlichKwong operator*
(
const scalar s,
const aungierRedlichKwong&
);
inline friend aungierRedlichKwong operator==
(
const aungierRedlichKwong&,
const aungierRedlichKwong&
);
*/
// Friend operators
// Ostream Operator

93
src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwongI.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,7 +28,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "aungierRedlichKwong.H"
@ -38,8 +37,6 @@ Germany
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
@ -62,7 +59,7 @@ inline aungierRedlichKwong::aungierRedlichKwong
c_(this->RR*Tcrit_/(pcrit_+(a_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_),
n_(0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())))
rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
{}
@ -103,6 +100,7 @@ inline scalar aungierRedlichKwong::rhostd()const
return rhostd_;
}
//returns the pressure for a given density and temperature
inline scalar aungierRedlichKwong::p(const scalar rho,const scalar T) const
{
@ -163,7 +161,6 @@ inline scalar aungierRedlichKwong::integral_p_dv
}
//needed to calculate the entropy
//(molar values)
inline scalar aungierRedlichKwong::integral_dpdT_dv
@ -186,8 +183,6 @@ inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
}
//(molar values)
inline scalar aungierRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
{
@ -198,8 +193,6 @@ inline scalar aungierRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
}
//(molar values)
//using second order implicit differentiation
inline scalar aungierRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
@ -214,8 +207,6 @@ inline scalar aungierRedlichKwong::d2vdT2(const scalar rho, const scalar T) cons
}
//(molar values)
inline scalar aungierRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
{
@ -225,8 +216,6 @@ inline scalar aungierRedlichKwong::d2pdvdT(const scalar rho, const scalar T) con
}
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar aungierRedlichKwong::integral_d2pdT2_dv
@ -292,8 +281,9 @@ inline scalar aungierRedlichKwong::rho(
i++;
if (i>8)
{
//using bisection methode as backup,
//solution must be between rho=0.001 to rho=1500;
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
//CL: if not, change rhoMax_ and rhoMin_
for(i=0; i<200; i++)
{
scalar f1 = this->p(rho1,T) - p;
@ -347,6 +337,7 @@ inline scalar aungierRedlichKwong::rho(
return this->W()/molarVolume;
}
//- Return density [kg/m^3]
inline scalar aungierRedlichKwong::rho(const scalar p,const scalar T) const
{
@ -354,88 +345,24 @@ inline scalar aungierRedlichKwong::rho(const scalar p,const scalar T) const
return rho(p,T,p/(this->R()*T));
}
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar aungierRedlichKwong::psi(const scalar rho, const scalar T) const
{
return -this->dvdp(rho,T)*pow(rho,2)/this->W();
}
//- Return compression factor []
inline scalar aungierRedlichKwong::Z( const scalar p, const scalar T,const scalar rho0) const
{
return (p*this->rho(p,T,rho0))/(this->R()*T);
return p/(this->R()*T*this->rho(p,T,rho0));
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void aungierRedlichKwong::operator+=(const aungierRedlichKwong& pg)
{
specie::operator+=(pg);
}
inline void aungierRedlichKwong::operator-=(const aungierRedlichKwong& pg)
{
specie::operator-=(pg);
}
inline void aungierRedlichKwong::operator*=(const scalar s)
{
specie::operator*=(s);
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
//****************not working**************//
/*
inline aungierRedlichKwong operator+
(
const aungierRedlichKwong& pg1,
const aungierRedlichKwong& pg2
)
{
return aungierRedlichKwong
(
static_cast<const specie&>(pg1)
+ static_cast<const specie&>(pg2)
);
}
inline aungierRedlichKwong operator-
(
const aungierRedlichKwong& pg1,
const aungierRedlichKwong& pg2
)
{
return aungierRedlichKwong
(
static_cast<const specie&>(pg1)
- static_cast<const specie&>(pg2)
);
}
inline aungierRedlichKwong operator*
(
const scalar s,
const aungierRedlichKwong& pg
)
{
return aungierRedlichKwong(s*static_cast<const specie&>(pg));
}
inline aungierRedlichKwong operator==
(
const aungierRedlichKwong& pg1,
const aungierRedlichKwong& pg2
)
{
return pg2 - pg1;
}
*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

14
src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinson.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -61,22 +61,22 @@ pengRobinson::pengRobinson(Istream& is)
pcrit_(readScalar(is)),
Tcrit_(readScalar(is)),
azentricFactor_(readScalar(is)),
a_(0.457235*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
a0_(0.457235*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
b_(0.077796*this->RR*Tcrit_/pcrit_),
n_(0.37464+1.54226*azentricFactor_-0.26992*pow(azentricFactor_,2)),
rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())))
TSave(0.0),
rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
{
is.check("pengRobinson::pengRobinson(Istream& is)");
}
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const pengRobinson& pg)
Ostream& operator<<(Ostream& os, const pengRobinson& pr)
{
os << static_cast<const specie&>(pg)<< tab
<< pg.pcrit_ << tab<< pg.Tcrit_<< tab << pg.azentricFactor_;
os << static_cast<const specie&>(pr)<< tab
<< pr.pcrit_ << tab<< pr.Tcrit_<< tab << pr.azentricFactor_;
os.check("Ostream& operator<<(Ostream& os, const pengRobinson& st)");
return os;

105
src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinson.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,6 +28,12 @@ Class
Description
Peng Robinson equation of state.
Paper:
Title: A New Two-Constant Equation of State
Authors: Ding-Yu Peng and Donald B. Robinson
Journal: Ind. Eng. Chem., Fundam., Vol. 15, No. 1, 1976
SourceFiles
pengRobinsonI.H
pengRobinson.C
@ -39,7 +45,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef pengRobinson_H
@ -63,19 +68,6 @@ class pengRobinson
public specie
{
//Member Variabels
scalar pcrit_;
scalar Tcrit_;
scalar azentricFactor_;
//-Peng Robinson factors
scalar a_;
scalar b_;
scalar n_;
//- Density @STD, initialise after a, b!
scalar rhostd_;
// Static Private data
//should be read from the fvSolution file where rhoMax and rhoMin values must be defined (for rhoSimpleFoam)
@ -83,6 +75,34 @@ class pengRobinson
static const scalar rhoMax_;
static const scalar rhoMin_;
protected:
// Protected data
scalar pcrit_;
scalar Tcrit_;
scalar azentricFactor_;
//-Peng Robinson factors
scalar n_;
scalar a0_;
mutable scalar b_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures
mutable scalar aSave;
mutable scalar daSave;
mutable scalar d2aSave;
//CL: save the temperature for which the save coefficients (amix,dadTmix,d2adT2mix) are correct
mutable scalar TSave;
//- Density @STD, initialise after a0, b!
mutable scalar rhostd_;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;
public:
@ -90,11 +110,9 @@ public:
// Constructors
//- Construct from components
inline pengRobinson(
const specie& sp,
scalar pcrit,
scalar Tcrit,
scalar azentricFactor
inline pengRobinson
(
const specie& sp
);
//- Construct from Istream
@ -109,15 +127,28 @@ public:
// Selector from Istream
inline static autoPtr<pengRobinson> New(Istream& is);
// Member functions
inline scalar rhostd()const;
//CL: Model coefficient a(T)
inline scalar a(const scalar T)const;
//CL: temperature deriviative of model coefficient a(T)
inline scalar dadT(const scalar T)const;
//CL: second order temperature deriviative of model coefficient a(T)
inline scalar d2adT2(const scalar T)const;
inline scalar a0()const;
inline scalar b()const;
inline scalar n()const;
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
//first order derivatives
//CL: first order derivatives
inline scalar dpdv(const scalar rho,const scalar T) const;
inline scalar dpdT(const scalar rho, const scalar T) const;
@ -138,14 +169,14 @@ public:
const scalar T
) const;
// Used for cv
//CL: Used for cv
inline scalar integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const ;
// second order derivatives, not Used At The Moment
//CL: second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho,const scalar T) const;
inline scalar d2pdT2(const scalar rho,const scalar T) const;
@ -154,10 +185,10 @@ public:
inline scalar d2vdT2(const scalar rho,const scalar T) const;
//Used for internal Energy
//CL: Used for internal Energy
inline scalar integral_p_dv(const scalar rho,const scalar T) const;
//Used for Entropy
//CL: Used for Entropy
inline scalar integral_dpdT_dv(const scalar rho,const scalar T) const;
//- Return density [kg/m^3]
@ -185,13 +216,16 @@ public:
// Member operators
/*
inline void operator+=(const pengRobinson&);
inline void operator-=(const pengRobinson&);
inline void operator*=(const scalar);
*/
/* // Friend operators
// Friend operators
inline friend pengRobinson operator+
(
@ -199,25 +233,12 @@ public:
const pengRobinson&
);
inline friend pengRobinson operator-
(
const pengRobinson&,
const pengRobinson&
);
inline friend pengRobinson operator*
(
const scalar s,
const pengRobinson&
);
inline friend pengRobinson operator==
(
const pengRobinson&,
const pengRobinson&
);
*/
// Ostream Operator
friend Ostream& operator<<(Ostream&, const pengRobinson&);

332
src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinsonI.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,7 +28,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "pengRobinson.H"
@ -38,40 +37,31 @@ Germany
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
inline pengRobinson::pengRobinson
(
const specie& sp,
scalar pcrit,
scalar Tcrit,
scalar azentricFactor
const specie& sp
)
:
specie(sp),
pcrit_(pcrit),
Tcrit_(Tcrit),
azentricFactor_(azentricFactor),
a_(0.457235*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
b_(0.077796*this->RR*Tcrit_/pcrit_),
n_(0.37464+1.54226*azentricFactor_-0.26992*pow(azentricFactor_,2)),
rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())))
TSave(0)
{}
// Construct as named copy
inline pengRobinson::pengRobinson(const word& name, const pengRobinson& pg)
inline pengRobinson::pengRobinson(const word& name, const pengRobinson& pr)
:
specie(name, pg),
pcrit_(pg.pcrit_),
Tcrit_(pg.Tcrit_),
azentricFactor_(pg.azentricFactor_),
a_(pg.a_),
b_(pg.b_),
n_(pg.n_),
rhostd_(pg.rhostd_)
specie(name, pr),
pcrit_(pr.pcrit_),
Tcrit_(pr.Tcrit_),
azentricFactor_(pr.azentricFactor_),
a0_(pr.a0_),
b_(pr.b_),
n_(pr.n_),
rhostd_(pr.rhostd_),
TSave(0)
{}
@ -91,19 +81,98 @@ inline autoPtr<pengRobinson> pengRobinson::New(Istream& is)
// * * * * * * * * * * * * * Member Functions * * * * * * * * * * * //
inline void pengRobinson::updateModelCoefficients(const scalar T)const
{
aSave=a0_*pow(1+n_*(1-pow(T/Tcrit_,0.5)),2);
daSave=a0_*n_*(n_*sqrt(T/Tcrit_)-n_-1)*sqrt(T/Tcrit_)/T;
d2aSave=a0_*n_*(n_+1)*sqrt(T/Tcrit_)/(2*pow(T,2));
//CL: saving the temperature at which the coefficients are valid
TSave=T;
}
inline scalar pengRobinson::rhostd()const
{
return rhostd_;
}
//CL: Model coefficient a(T)
inline scalar pengRobinson::a(const scalar T)const
{
//CL: check if a has already been calculated for this temperature
if(TSave==T)
{
return aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return aSave;
}
}
//CL: temperature deriviative of model coefficient a(T)
inline scalar pengRobinson::dadT(const scalar T)const
{
// check if a has already been calculated for this temperature
if(TSave==T)
{
return daSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return daSave;
}
}
//CL: second order temperature deriviative of model coefficient a(T)
inline scalar pengRobinson::d2adT2(const scalar T)const
{
// check if a has already been calculated for this temperature
if(TSave==T)
{
return d2aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return d2aSave;
}
}
inline scalar pengRobinson::a0()const
{
return a0_;
}
inline scalar pengRobinson::b()const
{
return b_;
}
inline scalar pengRobinson::n()const
{
return n_;
}
//returns the pressure for a given density and temperature
inline scalar pengRobinson::p(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR*T/(Vm-b_)
-(a_*pow((1+n_*(1-pow((T/Tcrit_),0.5))),2))
/(pow(Vm,2)+2*b_*Vm-pow(b_,2));
return this->RR*T/(Vm-b_)-a(T)/(pow(Vm,2)+2*b_*Vm-pow(b_,2));
}
@ -112,20 +181,19 @@ inline scalar pengRobinson::p(const scalar rho,const scalar T) const
inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return
-(
4*a_*n_*Tcrit_*(b_-Vm)*(pow(b_,2)
-pow(Vm,2))*(n_+1)*pow((T/Tcrit_),0.5)
+Tcrit_*(-2*a_*pow((n_+1),2)*(pow(b_,3)-pow(b_,2)*Vm
-b_*pow(Vm,2)+pow(Vm,3))
+this->RR*T*(pow(b_,4)-4*pow(b_,3)*Vm
+2*pow(b_,2)*pow(Vm,2)
+4*b_*pow(Vm,3)+pow(Vm,4)))
-2*a_*pow(n_,2)*T*(pow(b_,3)-pow(b_,2)*Vm
-b_*pow(Vm,2)+pow(Vm,3))
)
/(Tcrit_*pow((b_-Vm),2)*pow((pow(b_,2)
-2*b_*Vm-pow(Vm,2)),2));
return(
2*a(T)*
(
pow(b_,3)-pow(b_,2)*Vm-b_*pow(Vm,2)+pow(Vm,3)
)
-this->RR*T*
(
pow(b_,4)-4*pow(b_,3)*Vm+2*pow(b_,2)*pow(Vm,2)
+4*b_*pow(Vm,3)+pow(Vm,4)
)
)
/(pow(b_-Vm,2)*pow(pow(b_,2)-2*b_*Vm-pow(Vm,2),2));
}
@ -134,15 +202,8 @@ inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
inline scalar pengRobinson::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return
(
-a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
/(T*(pow(b_,2)
-2*b_*Vm-pow(Vm,2)))
+a_*pow(n_,2)/(Tcrit_*(pow(b_,2)
-2*b_*Vm-pow(Vm,2)))
-this->RR/(b_-Vm)
);
return this->RR/(Vm-b_)-dadT(T)/(pow(Vm,2)+2*b_*Vm-pow(b_,2));
}
@ -167,14 +228,8 @@ inline scalar pengRobinson::dvdp(const scalar rho,const scalar T) const
inline scalar pengRobinson::integral_p_dv(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return pow(2,0.5)*a_*(2*n_*Tcrit_*(n_+1)*pow(T/Tcrit_,0.5)
-Tcrit_*(pow(n_,2)+2*n_+1)-pow(n_,2)*T)
*log(b_*(1-pow(2,0.5))+Vm)/(4*b_*Tcrit_)
+this->RR*T*log(Vm-b_)
-pow(2,0.5)*a_*(2*n_*Tcrit_*(n_+1)*pow(T/Tcrit_,0.5)
-Tcrit_*(pow(n_,2)+2*n_+1)
-pow(n_,2)*T)*log(b_*(pow(2,0.5)+1)+Vm)
/(4*b_*Tcrit_);
return -pow(2,0.5)*a(T)*log(b_*(1-pow(2,0.5))+Vm)/(4*b_)+this->RR*T*log(Vm-b_)+pow(2,0.5)*a(T)*log(b_*(pow(2,0.5)+1)+Vm)/(4*b_);
}
@ -183,13 +238,8 @@ inline scalar pengRobinson::integral_p_dv(const scalar rho,const scalar T) const
inline scalar pengRobinson::integral_dpdT_dv(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return (pow(2,0.5)*a_*n_*(n_+1)*pow(T/Tcrit_,0.5)/(4*b_*T)
-pow(2,0.5)*a_*pow(n_,2)/(4*b_*Tcrit_))
*log(b_*(1-pow(2,0.5))+Vm)
+this->RR*log(Vm-b_)
+(pow(2,0.5)*a_*pow(n_,2)/(4*b_*Tcrit_)
-pow(2,0.5)*a_*n_*(n_+1)*pow(T/Tcrit_,0.5)/(4*b_*T))
*log(b_*(pow(2,0.5)+1)+Vm);
return -pow(2,0.5)*dadT(T)*log(b_*(1-pow(2,0.5))+Vm)/(4*b_)+this->RR*log(Vm-b_)+pow(2,0.5)*dadT(T)*log(b_*(pow(2,0.5)+1)+Vm)/(4*b_);
}
@ -197,10 +247,8 @@ inline scalar pengRobinson::integral_dpdT_dv(const scalar rho,const scalar T) co
inline scalar pengRobinson::d2pdT2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return
a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
/(2*pow(T,2)*(pow(b_,2)
-2*b_*Vm-pow(Vm,2)));
return -d2adT2(T)/(pow(Vm,2)+2*b_*Vm-pow(b_,2));
}
@ -208,51 +256,20 @@ inline scalar pengRobinson::d2pdT2(const scalar rho,const scalar T) const
inline scalar pengRobinson::d2pdv2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return
2*(
2*a_*n_*Tcrit_*(b_-Vm)*(n_+1)*
(
5*pow(b_,4)-4*pow(b_,3)*Vm
-4*pow(b_,2)*pow(Vm,2)
+3*pow(Vm,4)
)
*pow(T/Tcrit_,0.5)
+Tcrit_*
(
this->RR*T*
(
pow(b_,6)
-6*pow(b_,5)*Vm
+9*pow(b_,4)*pow(Vm,2)
+4*pow(b_,3)*pow(Vm,3)
-9*pow(b_,2)*pow(Vm,4)
-6*b_*pow(Vm,5)
-pow(Vm,6)
)
-a_*pow((n_+1),2)*
(
5*pow(b_,5)
-9*pow(b_,4)*Vm
+4*pow(b_,2)*pow(Vm,3)
+3*b_*pow(Vm,4)-3*pow(Vm,5)
)
)
-a_*pow(n_,2)*T*
(
5*pow(b_,5)
-9*pow(b_,4)*Vm
+4*pow(b_,2)*pow(Vm,3)
+3*b_*pow(Vm,4)
-3*pow(Vm,5)
)
)
/
(
Tcrit_*
pow((pow(b_,2)
-2*b_*Vm-pow(Vm,2)),3)
*pow((Vm-b_),3)
);
return 2*
(
a(T)*
(
5*pow(b_,5)-9*pow(b_,4)*Vm+4*pow(b_,2)*pow(Vm,3)+3*b_*pow(Vm,4)-3*pow(Vm,5)
)
-this->RR*T*
(
pow(b_,6)-6*pow(b_,5)*Vm+9*pow(b_,4)*pow(Vm,2)+4*pow(b_,3)*pow(Vm,3)
-9*pow(b_,2)*pow(Vm,4)-6*b_*pow(Vm,5)-pow(Vm,6)
)
)
/(pow(b_-Vm,3)*pow(pow(b_,2)-2*b_*Vm-pow(Vm,2),3));
}
@ -274,11 +291,19 @@ inline scalar pengRobinson::d2vdT2(const scalar rho, const scalar T) const
inline scalar pengRobinson::d2pdvdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return -2*a_*n_*(b_+Vm)*(n_+1)*pow(T/Tcrit_,0.5)
/(T*pow((pow(b_,2)-2*b_*Vm-pow(Vm,2)),2))
+(2*a_*pow(n_,2)*(b_+Vm))
/(Tcrit_*pow(pow(b_,2)-2*b_*Vm-pow(Vm,2),2))
-this->RR/(pow(b_-Vm,2));
return(
2*dadT(T)*
(
pow(b_,3)-pow(b_,2)*Vm-b_*pow(Vm,2)+pow(Vm,3)
)
-this->RR*
(
pow(b_,4)-4*pow(b_,3)*Vm+2*pow(b_,2)*pow(Vm,2)
+4*b_*pow(Vm,3)+pow(Vm,4)
)
)
/(pow(b_-Vm,2)*pow(pow(b_,2)-2*b_*Vm-pow(Vm,2),2));
}
@ -287,10 +312,8 @@ inline scalar pengRobinson::d2pdvdT(const scalar rho, const scalar T) const
inline scalar pengRobinson::integral_d2pdT2_dv(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return pow(2,0.5)*a_*n_*(n_+1)*pow(T/Tcrit_,0.5)
*log(b_*(pow(2,0.5)+1)+Vm)/(8*b_*pow(T,2))
-pow(2,0.5)*a_*n_*(n_+1)*pow(T/Tcrit_,0.5)
*log(b_*(1-pow(2,0.5))+Vm)/(8*b_*pow(T,2));
return pow(2,0.5)*d2adT2(T)*log(b_*(pow(2,0.5)+1)+Vm)/(4*b_)-pow(2,0.5)*d2adT2(T)*log(b_*(1-pow(2,0.5))+Vm)/(4*b_);
}
@ -307,7 +330,8 @@ inline scalar pengRobinson::isobarExpCoef(const scalar rho,const scalar T) cons
inline scalar pengRobinson::isothermalCompressiblity(const scalar rho,const scalar T) const
{
return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
//also possible : return -this->dvdp(rho,T)*rho/this->W();
//CL: also possible
//CL: return -this->dvdp(rho,T)*rho/this->W();
}
@ -345,8 +369,9 @@ inline scalar pengRobinson::rho(
i++;
if (i>8)
{
//using bisection methode as backup,
//solution must be between rho=0.001 to rho=1500;
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
//CL: if not, change rhoMax_ and rhoMin_
for(i=0; i<200; i++)
{
scalar f1 = this->p(rho1,T) - p;
@ -419,20 +444,21 @@ inline scalar pengRobinson::psi(const scalar rho, const scalar T) const
//- Return compression factor []
inline scalar pengRobinson::Z( const scalar p, const scalar T,const scalar rho0) const
{
return (p*this->rho(p,T,rho0))/(this->R()*T);
return p/(this->R()*T*this->rho(p,T,rho0));
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void pengRobinson::operator+=(const pengRobinson& pg)
/*
inline void pengRobinson::operator+=(const pengRobinson& pr)
{
specie::operator+=(pg);
specie::operator+=(pr);
}
inline void pengRobinson::operator-=(const pengRobinson& pg)
inline void pengRobinson::operator-=(const pengRobinson& pr)
{
specie::operator-=(pg);
specie::operator-=(pr);
}
inline void pengRobinson::operator*=(const scalar s)
@ -440,59 +466,31 @@ inline void pengRobinson::operator*=(const scalar s)
specie::operator*=(s);
}
*/
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
//****************not working**************//
/*
inline pengRobinson operator+
(
const pengRobinson& pg1,
const pengRobinson& pg2
const pengRobinson& pr1,
const pengRobinson& pr2
)
{
return pengRobinson
(
static_cast<const specie&>(pg1)
+ static_cast<const specie&>(pg2)
static_cast<const specie&>(pr1)
+ static_cast<const specie&>(pr2)
);
}
inline pengRobinson operator-
(
const pengRobinson& pg1,
const pengRobinson& pg2
)
{
return pengRobinson
(
static_cast<const specie&>(pg1)
- static_cast<const specie&>(pg2)
);
}
inline pengRobinson operator*
(
const scalar s,
const pengRobinson& pg
const pengRobinson& pr
)
{
return pengRobinson(s*static_cast<const specie&>(pg));
return pengRobinson(s*static_cast<const specie&>(pr));
}
inline pengRobinson operator==
(
const pengRobinson& pg1,
const pengRobinson& pg2
)
{
return pg2 - pg1;
}
*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

11
src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,7 +31,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "redlichKwong.H"
@ -64,7 +63,7 @@ redlichKwong::redlichKwong(Istream& is)
a_(0.42748*pow(this->RR,2)*pow(Tcrit_,2.5)/pcrit_),
b_(0.08664*this->RR*Tcrit_/pcrit_),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R())))
rhostd_(this->rho(Pstd, Tstd, Pstd/(Tstd*this->R())))
{
is.check("redlichKwong::redlichKwong(Istream& is)");
}
@ -72,10 +71,10 @@ redlichKwong::redlichKwong(Istream& is)
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const redlichKwong& pg)
Ostream& operator<<(Ostream& os, const redlichKwong& rk)
{
os << static_cast<const specie&>(pg)<< tab
<< pg.pcrit_ << tab<< pg.Tcrit_;
os << static_cast<const specie&>(rk)<< tab
<< rk.pcrit_ << tab<< rk.Tcrit_;
os.check("Ostream& operator<<(Ostream& os, const redlichKwong& st)");
return os;

54
src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -38,7 +38,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef redlichKwong_H
@ -60,18 +59,6 @@ class redlichKwong
:
public specie
{
// Private data
scalar pcrit_;
scalar Tcrit_;
//-Redlich Kwong factors
scalar a_;
scalar b_;
//- Density @STD, initialise after a, b!
scalar rhostd_;
// Static Private data
//should be read from the fvSolution file where rhoMax and rhoMin values must be defined (for rhoSimpleFoam)
@ -79,18 +66,26 @@ class redlichKwong
static const scalar rhoMax_;
static const scalar rhoMin_;
protected:
// Protected data
scalar pcrit_;
scalar Tcrit_;
//-Redlich Kwong factors
mutable scalar a_;
mutable scalar b_;
//- Density @STD, initialise after a, b!
mutable scalar rhostd_;
public:
// Constructors
//- Construct from components
inline redlichKwong
(
const specie& sp,
scalar pcrit,
scalar Tcrit
const specie& sp
);
//- Construct from Istream
@ -109,6 +104,10 @@ public:
inline scalar rhostd() const;
inline scalar a() const;
inline scalar b() const;
inline scalar p(const scalar rho, const scalar T) const;
@ -173,16 +172,16 @@ public:
const scalar rho0
) const;
/*
// Member operators
inline void operator+=(const redlichKwong&);
inline void operator-=(const redlichKwong&);
inline void operator*=(const scalar);
*/
/* // Friend operators
// Friend operators
inline friend redlichKwong operator+
(
@ -190,25 +189,12 @@ public:
const redlichKwong&
);
inline friend redlichKwong operator-
(
const redlichKwong&,
const redlichKwong&
);
inline friend redlichKwong operator*
(
const scalar s,
const redlichKwong&
);
inline friend redlichKwong operator==
(
const redlichKwong&,
const redlichKwong&
);
*/
// Ostream Operator
friend Ostream& operator<<(Ostream&, const redlichKwong&);

105
src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwongI.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,7 +28,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "redlichKwong.H"
@ -43,30 +42,22 @@ namespace Foam
// Construct from components
inline redlichKwong::redlichKwong
(
const specie& sp,
scalar pcrit,
scalar Tcrit
const specie& sp
)
:
specie(sp),
pcrit_(pcrit),
Tcrit_(Tcrit),
a_(0.42748*pow(this->RR,2)*pow(Tcrit_,2.5)/pcrit_),
b_(0.08664*this->RR*Tcrit_/pcrit_),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R())))
specie(sp)
{}
// Construct as named copy
inline redlichKwong::redlichKwong(const word& name, const redlichKwong& pg)
inline redlichKwong::redlichKwong(const word& name, const redlichKwong& rk)
:
specie(name, pg),
pcrit_(pg.pcrit_),
Tcrit_(pg.Tcrit_),
a_(pg.a_),
b_(pg.b_),
rhostd_(pg.rhostd_)
specie(name, rk),
pcrit_(rk.pcrit_),
Tcrit_(rk.Tcrit_),
a_(rk.a_),
b_(rk.b_),
rhostd_(rk.rhostd_)
{}
@ -84,16 +75,26 @@ inline autoPtr<redlichKwong> redlichKwong::New(Istream& is)
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline scalar redlichKwong::rhostd()const
{
return rhostd_;
}
inline scalar redlichKwong::a()const
{
return a_;
}
inline scalar redlichKwong::b()const
{
return b_;
}
//returns the pressure for a given density and temperature
inline scalar redlichKwong::p(const scalar rho, const scalar T) const
{
@ -300,8 +301,9 @@ inline scalar redlichKwong::rho
i++;
if (i>8)
{
//using bisection methode as backup,
//solution must be between rho=0.001 to rho=1500;
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
//CL: if not, change rhoMax_ and rhoMin_
for(i=0; i<200; i++)
{
scalar f1 = this->p(rho1,T) - p;
@ -379,22 +381,22 @@ inline scalar redlichKwong::Z
const scalar rho0
) const
{
return p*this->rho(p,T,rho0)/(this->R()*T);
return p/(this->R()*T*this->rho(p,T,rho0));
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
/*
inline void redlichKwong::operator+=(const redlichKwong& pg)
inline void redlichKwong::operator+=(const redlichKwong& rk)
{
specie::operator+=(pg);
specie::operator+=(rk);
}
inline void redlichKwong::operator-=(const redlichKwong& pg)
inline void redlichKwong::operator-=(const redlichKwong& rk)
{
specie::operator-=(pg);
specie::operator-=(rk);
}
@ -403,59 +405,32 @@ inline void redlichKwong::operator*=(const scalar s)
specie::operator*=(s);
}
*/
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
//****************not working**************////
/*
inline redlichKwong operator+
(
const redlichKwong& pg1,
const redlichKwong& pg2
const redlichKwong& rk1,
const redlichKwong& rk2
)
{
return redlichKwong
(
static_cast<const specie&>(pg1)
+ static_cast<const specie&>(pg2)
static_cast<const specie&>(rk1)
+ static_cast<const specie&>(rk2)
);
}
inline redlichKwong operator-
(
const redlichKwong& pg1,
const redlichKwong& pg2
)
{
return redlichKwong
(
static_cast<const specie&>(pg1)
- static_cast<const specie&>(pg2)
);
}
inline redlichKwong operator*
(
const scalar s,
const redlichKwong& pg
const redlichKwong& rk
)
{
return redlichKwong(s*static_cast<const specie&>(pg));
return redlichKwong(s*static_cast<const specie&>(rk));
}
inline redlichKwong operator==
(
const redlichKwong& pg1,
const redlichKwong& pg2
)
{
return pg2 - pg1;
}
*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

12
src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwong.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -60,11 +60,11 @@ soaveRedlichKwong::soaveRedlichKwong(Istream& is)
pcrit_(readScalar(is)),
Tcrit_(readScalar(is)),
azentricFactor_(readScalar(is)),
a_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/(pcrit_)),
a0_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/(pcrit_)),
b_(0.08664*this->RR*Tcrit_/pcrit_),
n_(0.48508+1.55171*azentricFactor_-0.15613*pow(azentricFactor_,2)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())))
rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
{
is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)");
}
@ -72,10 +72,10 @@ soaveRedlichKwong::soaveRedlichKwong(Istream& is)
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const soaveRedlichKwong& pg)
Ostream& operator<<(Ostream& os, const soaveRedlichKwong& srk)
{
os << static_cast<const specie&>(pg)<< tab
<< pg.pcrit_ << tab<< pg.Tcrit_<<tab<<pg.azentricFactor_;
os << static_cast<const specie&>(srk)<< tab
<< srk.pcrit_ << tab<< srk.Tcrit_<<tab<<srk.azentricFactor_;
os.check("Ostream& operator<<(Ostream& os, const soaveRedlichKwong& st)");
return os;

105
src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwong.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,6 +28,11 @@ Class
Description
Soave Redlich Kwong equation of state.
Paper:
Title: Equilibrium Constants from a Modified Redlich-Kwong Equation of State
Authors: G. Soave
Journal: Chemical Engineering Science, vol. 27(6), 1972, pp. 1197-1203
SourceFiles
soaveRedlichKwongI.H
soaveRedlichKwong.C
@ -61,20 +66,6 @@ class soaveRedlichKwong
public specie
{
// Private data
scalar pcrit_;
scalar Tcrit_;
scalar azentricFactor_;
//-Soave Redlich Kwong factors
scalar a_;
scalar b_;
scalar n_;
//- Density @STD, initialise after a, b!
scalar rhostd_;
// Static Private data
//should be read from the fvSolution file where rhoMax and rhoMin values must be defined (for rhoSimpleFoam)
@ -82,18 +73,46 @@ class soaveRedlichKwong
static const scalar rhoMax_;
static const scalar rhoMin_;
protected:
// Protected data
scalar pcrit_;
scalar Tcrit_;
scalar azentricFactor_;
//-Soave Redlich Kwong
scalar n_;
scalar a0_;
mutable scalar b_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures
mutable scalar aSave;
mutable scalar daSave;
mutable scalar d2aSave;
//CL: save the temperature for which the save coefficients (amix,dadTmix,d2adT2mix) are correct
mutable scalar TSave;
//- Density @STD, initialise after a0, b!
mutable scalar rhostd_;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;
public:
// Constructors
//- Construct from components
inline soaveRedlichKwong(
const specie& sp,
scalar pcrit,
scalar Tcrit,
scalar azentricFactor
);
inline soaveRedlichKwong
(
const specie& s
);
//- Construct from Istream
soaveRedlichKwong(Istream&);
@ -110,10 +129,25 @@ public:
// Member functions
inline scalar rhostd()const;
//CL: Model coefficient a(T)
inline scalar a(const scalar T)const;
//CL: temperature deriviative of model coefficient a(T)
inline scalar dadT(const scalar T)const;
//CL: second order temperature deriviative of model coefficient a(T)
inline scalar d2adT2(const scalar T)const;
inline scalar a0()const;
inline scalar b()const;
inline scalar n()const;
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
//first order derivatives
//CL: first order derivatives
inline scalar dpdv(const scalar rho,const scalar T) const;
inline scalar dpdT(const scalar rho, const scalar T) const;
@ -134,10 +168,10 @@ public:
const scalar T
) const;
// Used for cv
//CL: Used for cv
inline scalar integral_d2pdT2_dv(const scalar rho,const scalar T) const ;
// second order derivatives, not Used At The Moment
//CL: second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho,const scalar T) const;
inline scalar d2pdT2(const scalar rho,const scalar T) const;
@ -146,7 +180,7 @@ public:
inline scalar d2vdT2(const scalar rho,const scalar T) const;
//Used for internal Energy
//CL: Used for internal Energy
inline scalar integral_p_dv(const scalar rho,const scalar T) const;
//Used for Entropy
@ -177,13 +211,15 @@ public:
// Member operators
/*
inline void operator+=(const soaveRedlichKwong&);
inline void operator-=(const soaveRedlichKwong&);
inline void operator*=(const scalar);
*/
/* // Friend operators
// Friend operators
inline friend soaveRedlichKwong operator+
(
@ -191,25 +227,12 @@ public:
const soaveRedlichKwong&
);
inline friend soaveRedlichKwong operator-
(
const soaveRedlichKwong&,
const soaveRedlichKwong&
);
inline friend soaveRedlichKwong operator*
(
const scalar s,
const soaveRedlichKwong&
);
inline friend soaveRedlichKwong operator==
(
const soaveRedlichKwong&,
const soaveRedlichKwong&
);
*/
// Ostream Operator
friend Ostream& operator<<(Ostream&, const soaveRedlichKwong&);

318
src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwongI.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -43,38 +43,26 @@ namespace Foam
// Construct from components
inline soaveRedlichKwong::soaveRedlichKwong
(
const specie& sp,
scalar pcrit,
scalar Tcrit,
scalar azentricFactor
const specie& sp
)
:
specie(sp),
pcrit_(pcrit),
Tcrit_(Tcrit),
azentricFactor_(azentricFactor),
a_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/(pcrit_)),
b_(0.08664*this->RR*Tcrit_/pcrit_),
n_(0.48508+1.55171*azentricFactor_-0.15613*pow(azentricFactor_,2)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())))
TSave(0)
{}
// Construct as named copy
inline soaveRedlichKwong::soaveRedlichKwong(const word& name,const soaveRedlichKwong& pg)
inline soaveRedlichKwong::soaveRedlichKwong(const word& name,const soaveRedlichKwong& srk)
:
specie(name, pg),
pcrit_(pg.pcrit_),
Tcrit_(pg.Tcrit_),
azentricFactor_(pg.azentricFactor_),
a_(pg.a_),
b_(pg.b_),
n_(pg.n_),
rhostd_(pg.rhostd_)
specie(name, srk),
pcrit_(srk.pcrit_),
Tcrit_(srk.Tcrit_),
azentricFactor_(srk.azentricFactor_),
a0_(srk.a0_),
b_(srk.b_),
n_(srk.n_),
rhostd_(srk.rhostd_),
TSave(0)
{}
@ -94,21 +82,101 @@ inline autoPtr<soaveRedlichKwong> soaveRedlichKwong::New(Istream& is)
// * * * * * * * * * * * * * Member Functions * * * * * * * * * * * //
inline void soaveRedlichKwong::updateModelCoefficients(const scalar T)const
{
aSave=a0_*pow(1+n_*(1-pow(T/Tcrit_,0.5)),2);
daSave=a0_*n_*(n_*sqrt(T/Tcrit_)-n_-1)*sqrt(T/Tcrit_)/T;
d2aSave=a0_*n_*(n_+1)*sqrt(T/Tcrit_)/(2*pow(T,2));
//CL: saving the temperature at which the coefficients are valid
TSave=T;
}
inline scalar soaveRedlichKwong::rhostd()const
{
return rhostd_;
}
//CL: Model coefficient a(T)
inline scalar soaveRedlichKwong::a(const scalar T)const
{
//CL: check if a has already been calculated for this temperature
if(TSave==T)
{
return aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return aSave;
}
}
//CL: temperature deriviative of model coefficient a(T)
inline scalar soaveRedlichKwong::dadT(const scalar T)const
{
// check if a has already been calculated for this temperature
if(TSave==T)
{
return daSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return daSave;
}
}
//CL: second order temperature deriviative of model coefficient a(T)
inline scalar soaveRedlichKwong::d2adT2(const scalar T)const
{
// check if a has already been calculated for this temperature
if(TSave==T)
{
return d2aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return d2aSave;
}
}
inline scalar soaveRedlichKwong::a0()const
{
return a0_;
}
inline scalar soaveRedlichKwong::b()const
{
return b_;
}
inline scalar soaveRedlichKwong::n()const
{
return n_;
}
//returns the pressure for a given density and temperature
inline scalar soaveRedlichKwong::p(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return
(
this->RR*T/(Vm-b_)
-a_*pow((1+n_*(1-pow((T/Tcrit_),0.5))),2)
/(Vm*(Vm+b_))
-a(T)/(Vm*(Vm+b_))
);
}
@ -118,27 +186,13 @@ inline scalar soaveRedlichKwong::p(const scalar rho,const scalar T) const
inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return
-(
2*a_*n_*Tcrit_*(b_-Vm)*(pow(b_,2)
+b_*Vm-2*pow(Vm,2))*(n_+1)*pow((T/Tcrit_),0.5)
+Tcrit_*
(
this->RR*T*pow(Vm,2)*(pow(b_,2)
+2*b_*Vm+pow(Vm,2))
-a_*
(pow(b_,3)-3*b_*pow(Vm,2)+
2*pow(Vm,3))
*pow((n_+1),2)
)
-a_*pow(n_,2)*T*(pow(b_,3)
-3*b_*pow(Vm,2)+2*pow(Vm,3))
)
/
(
pow(Vm,2)*Tcrit_*pow((b_+Vm),2)
*pow((b_-Vm),2)
);
(
a(T)*(pow(b_,3)-3*b_*pow(Vm,2)+2*pow(Vm,3))
-this->RR*T*pow(Vm,2)*(pow(b_,2)+2*b_*Vm+pow(Vm,2))
)
/(pow(Vm,2)*pow(b_+Vm,2)*pow(b_-Vm,2));
}
@ -147,11 +201,12 @@ inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
inline scalar soaveRedlichKwong::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
/(T*Vm*(b_+Vm))
-a_*pow(n_,2)
/(Vm*Tcrit_*(b_+Vm))
-this->RR/(b_-Vm);
return
(
this->RR/(Vm-b_)
-dadT(T)/(Vm*(Vm+b_))
);
}
@ -177,19 +232,8 @@ inline scalar soaveRedlichKwong::dvdp(const scalar rho,const scalar T) const
inline scalar soaveRedlichKwong::integral_p_dv(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR*T*log(Vm-b_)
-(
a_*(2*n_*Tcrit_*(n_+1)*pow(T/Tcrit_,0.5)
-Tcrit_*(pow(n_,2)+2*n_+1)
-pow(n_,2)*T)*log(b_+Vm)
)
/(b_*Tcrit_)
+a_*
(
2*n_*Tcrit_*(n_+1)*pow(T/Tcrit_,0.5)
-Tcrit_*(pow(n_,2)+2*n_+1)-pow(n_,2)*T
)
*log(Vm)/(b_*Tcrit_);
return this->RR*T*log(Vm-b_)+a(T)*log(b_+Vm)/b_-a(T)*log(Vm)/b_;
}
@ -198,15 +242,8 @@ inline scalar soaveRedlichKwong::integral_p_dv(const scalar rho,const scalar T)
inline scalar soaveRedlichKwong::integral_dpdT_dv(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR*log(Vm-b_)
+(
pow(n_,2)*a_/(b_*Tcrit_)
-a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
/(b_*T)
)*log(b_+Vm)
+(a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
/(b_*T)
-a_*pow(n_,2)/(b_*Tcrit_))*log(Vm);
return this->RR*log(Vm-b_)+dadT(T)*log(b_+Vm)/b_-dadT(T)*log(Vm)/b_;
}
@ -214,8 +251,8 @@ inline scalar soaveRedlichKwong::integral_dpdT_dv(const scalar rho,const scalar
inline scalar soaveRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return -a_*n_*(n_+1)*pow(T/Tcrit_,0.5)
/(2*pow(T,2)*Vm*(b_+Vm));
return -d2adT2(T)/(Vm*(Vm+b_));
}
@ -223,45 +260,20 @@ inline scalar soaveRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
inline scalar soaveRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return
2*(
2*a_*n_*Tcrit_*(b_-Vm)*
2*
(
a(T)*
(
pow(b_,4)+pow(b_,3)*Vm
-2*pow(b_,2)*pow(Vm,2)
-3*b_*pow(Vm,3)
+3*pow(Vm,4)
pow(b_,5)-3*pow(b_,3)*pow(Vm,2)-pow(b_,2)*pow(Vm,3)+6*b_*pow(Vm,4)-3*pow(Vm,5)
)
*(n_+1)*pow(T/Tcrit_,0.5)
-Tcrit_*
+this->RR*T*pow(Vm,3)*
(
a_*
(
pow(b_,5)
-3*pow(b_,3)*pow(Vm,2)
-pow(b_,2)*pow(Vm,3)
+6*b_*pow(Vm,4)
-3*pow(Vm,5)
)*pow((n_+1),2)
+this->RR*T*pow(Vm,3)*
(
pow(b_,3)
+3*pow(b_,2)*Vm
+3*b_*pow(Vm,2)
+pow(Vm,3)
)
)
-a_*pow(n_,2)*T*
(
pow(b_,5)
-3*pow(b_,3)*pow(Vm,2)
-pow(b_,2)*pow(Vm,3)
+6*b_*pow(Vm,4)
-3*pow(Vm,5)
pow(b_,3)+3*pow(b_,2)*Vm+3*b_*pow(Vm,2)+pow(Vm,3)
)
)
/(pow(Vm,3)*Tcrit_*pow((Vm+b_),3)
*pow((b_-Vm),3));
/(pow(Vm,3)*pow(b_+Vm,3)*pow(Vm-b_,3));
}
@ -283,19 +295,13 @@ inline scalar soaveRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
inline scalar soaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return -a_*n_*(b_+2*Vm)*(n_+1)*pow(T/Tcrit_,0.5)
/(T*pow(Vm,2)*pow((b_+Vm),2))
-(
this->RR*pow(Vm,2)*Tcrit_*(pow(b_,2)
+2*b_*Vm+pow(Vm,2))
-a_*pow(n_,2)*
(
pow(b_,3)
-3*b_*pow(Vm,2)
+2*pow(Vm,3)
)
return
(
dadT(T)*(pow(b_,3)-3*b_*pow(Vm,2)+2*pow(Vm,3))
-this->RR*pow(Vm,2)*(pow(b_,2)+2*b_*Vm+pow(Vm,2))
)
/(pow(Vm,2)*Tcrit_*pow((b_+Vm),2)*pow((b_-Vm),2));
/(pow(Vm,2)*pow(b_+Vm,2)*pow(b_-Vm,2));
}
@ -304,10 +310,8 @@ inline scalar soaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
inline scalar soaveRedlichKwong::integral_d2pdT2_dv(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return a_*n_*(n_+1)*pow(T/Tcrit_,0.5)*log(b_+Vm)
/(2*b_*pow(T,2))
-a_*n_*(n_+1)*pow(T/Tcrit_,0.5)*log(Vm)
/(2*b_*pow(T,2));
return d2adT2(T)*log(b_+Vm)/b_-d2adT2(T)*log(Vm)/b_;
}
@ -361,8 +365,9 @@ inline scalar soaveRedlichKwong::rho(
i++;
if (i>8)
{
//using bisection methode as backup,
//solution must be between rho=0.001 to rho=1500;
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
//CL: if not, change rhoMax_ and rhoMin_
for(i=0; i<200; i++)
{
scalar f1 = this->p(rho1,T) - p;
@ -419,8 +424,7 @@ inline scalar soaveRedlichKwong::rho(
//- Return density [kg/m^3]on
inline scalar soaveRedlichKwong::rho(const scalar p,const scalar T) const
{
//using perfect gas equation as starting point
//CL: using perfect gas equation as starting point
return rho(p,T,p/(this->R()*T));
}
@ -435,19 +439,19 @@ inline scalar soaveRedlichKwong::psi(const scalar rho, const scalar T) const
//- Return compression factor []
inline scalar soaveRedlichKwong::Z( const scalar p, const scalar T,const scalar rho0) const
{
return (p*this->rho(p,T,rho0))/(this->R()*T);
return p/(this->R()*T*this->rho(p,T,rho0));
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void soaveRedlichKwong::operator+=(const soaveRedlichKwong& pg)
/*
inline void soaveRedlichKwong::operator+=(const soaveRedlichKwong& srk)
{
specie::operator+=(pg);
specie::operator+=(srk);
}
inline void soaveRedlichKwong::operator-=(const soaveRedlichKwong& pg)
inline void soaveRedlichKwong::operator-=(const soaveRedlichKwong& srk)
{
specie::operator-=(pg);
specie::operator-=(srk);
}
inline void soaveRedlichKwong::operator*=(const scalar s)
@ -455,59 +459,31 @@ inline void soaveRedlichKwong::operator*=(const scalar s)
specie::operator*=(s);
}
*/
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
//**************** not working**************//
/*
inline soaveRedlichKwong operator+
(
const soaveRedlichKwong& pg1,
const soaveRedlichKwong& pg2
const soaveRedlichKwong& srk1,
const soaveRedlichKwong& srk2
)
{
return soaveRedlichKwong
(
static_cast<const specie&>(pg1)
+ static_cast<const specie&>(pg2)
static_cast<const specie&>(srk1)
+ static_cast<const specie&>(srk2)
);
}
inline soaveRedlichKwong operator-
(
const soaveRedlichKwong& pg1,
const soaveRedlichKwong& pg2
)
{
return soaveRedlichKwong
(
static_cast<const specie&>(pg1)
- static_cast<const specie&>(pg2)
);
}
inline soaveRedlichKwong operator*
(
const scalar s,
const soaveRedlichKwong& pg
const soaveRedlichKwong& srk
)
{
return soaveRedlichKwong(s*static_cast<const specie&>(pg));
return soaveRedlichKwong(s*static_cast<const specie&>(srk));
}
inline soaveRedlichKwong operator==
(
const soaveRedlichKwong& pg1,
const soaveRedlichKwong& pg2
)
{
return pg2 - pg1;
}
*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

9
src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C
View file

@ -28,7 +28,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "nasaHeatCapacityPolynomial.H"
@ -65,12 +64,12 @@ template<class equationOfState>
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const nasaHeatCapacityPolynomial<equationOfState>& ct
const nasaHeatCapacityPolynomial<equationOfState>& np
)
{
os << static_cast<const equationOfState&>(ct) << tab
<< ct.a1_ << tab<< ct.a2_ << tab << ct.a3_ << tab << ct.a4_ << tab << ct.a5_ << tab << ct.a6_ << tab << ct.a7_ ;
os.check("Ostream& operator<<(Ostream& os, const nasaHeatCapacityPolynomial& ct)");
os << static_cast<const equationOfState&>(np) << tab
<< np.a1_ << tab<< np.a2_ << tab << np.a3_ << tab << np.a4_ << tab << np.a5_ << tab << np.a6_ << tab << np.a7_ ;
os.check("Ostream& operator<<(Ostream& os, const nasaHeatCapacityPolynomial& np)");
return os;
}

164
src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H
View file

@ -70,22 +70,22 @@ template<class equationOfState>
inline Foam::nasaHeatCapacityPolynomial<equationOfState>::nasaHeatCapacityPolynomial
(
const word& name,
const nasaHeatCapacityPolynomial& ct
const nasaHeatCapacityPolynomial& np
)
:
equationOfState(name, ct),
a1_(ct.a1_),
a2_(ct.a2_),
a3_(ct.a3_),
a4_(ct.a4_),
a5_(ct.a5_),
a6_(ct.a6_),
a7_(ct.a7_),
e0_std(ct.e0_std),
s0_std(ct.s0_std),
integral_p_dv_std(ct.integral_p_dv_std),
integral_dpdT_dv_std(ct.integral_dpdT_dv_std),
cp_std(ct.cp_std)
equationOfState(name, np),
a1_(np.a1_),
a2_(np.a2_),
a3_(np.a3_),
a4_(np.a4_),
a5_(np.a5_),
a6_(np.a6_),
a7_(np.a7_),
e0_std(np.e0_std),
s0_std(np.s0_std),
integral_p_dv_std(np.integral_p_dv_std),
integral_dpdT_dv_std(np.integral_dpdT_dv_std),
cp_std(np.cp_std)
{}
@ -316,47 +316,47 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s
template<class equationOfState>
inline void Foam::nasaHeatCapacityPolynomial<equationOfState>::operator+=
(
const nasaHeatCapacityPolynomial<equationOfState>& ct
const nasaHeatCapacityPolynomial<equationOfState>& np
)
{
scalar molr1 = this->nMoles();
equationOfState::operator+=(ct);
equationOfState::operator+=(np);
molr1 /= this->nMoles();
scalar molr2 = ct.nMoles()/this->nMoles();
scalar molr2 = np.nMoles()/this->nMoles();
a1_ = molr1*a1_ + molr2*ct.a1_;
a2_ = molr1*a2_ + molr2*ct.a2_;
a3_ = molr1*a3_ + molr2*ct.a3_;
a4_ = molr1*a4_ + molr2*ct.a4_;
a5_ = molr1*a5_ + molr2*ct.a5_;
a6_ = molr1*a6_ + molr2*ct.a6_;
a7_ = molr1*a7_ + molr2*ct.a7_;
a1_ = molr1*a1_ + molr2*np.a1_;
a2_ = molr1*a2_ + molr2*np.a2_;
a3_ = molr1*a3_ + molr2*np.a3_;
a4_ = molr1*a4_ + molr2*np.a4_;
a5_ = molr1*a5_ + molr2*np.a5_;
a6_ = molr1*a6_ + molr2*np.a6_;
a7_ = molr1*a7_ + molr2*np.a7_;
}
template<class equationOfState>
inline void Foam::nasaHeatCapacityPolynomial<equationOfState>::operator-=
(
const nasaHeatCapacityPolynomial<equationOfState>& ct
const nasaHeatCapacityPolynomial<equationOfState>& np
)
{
scalar molr1 = this->nMoles();
nasaHeatCapacityPolynomial::operator-=(ct);
nasaHeatCapacityPolynomial::operator-=(np);
molr1 /= this->nMoles();
scalar molr2 = ct.nMoles()/this->nMoles();
scalar molr2 = np.nMoles()/this->nMoles();
a1_ = molr1*a1_ - molr2*ct.a1_;
a2_ = molr1*a2_ - molr2*ct.a2_;
a3_ = molr1*a3_ - molr2*ct.a3_;
a4_ = molr1*a4_ - molr2*ct.a4_;
a5_ = molr1*a5_ - molr2*ct.a5_;
a6_ = molr1*a6_ - molr2*ct.a6_;
a7_ = molr1*a7_ - molr2*ct.a7_;
a1_ = molr1*a1_ - molr2*np.a1_;
a2_ = molr1*a2_ - molr2*np.a2_;
a3_ = molr1*a3_ - molr2*np.a3_;
a4_ = molr1*a4_ - molr2*np.a4_;
a5_ = molr1*a5_ - molr2*np.a5_;
a6_ = molr1*a6_ - molr2*np.a6_;
a7_ = molr1*a7_ - molr2*np.a7_;
}
@ -365,33 +365,33 @@ inline void Foam::nasaHeatCapacityPolynomial<equationOfState>::operator-=
template<class equationOfState>
inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator+
(
const nasaHeatCapacityPolynomial<equationOfState>& ct1,
const nasaHeatCapacityPolynomial<equationOfState>& ct2
const nasaHeatCapacityPolynomial<equationOfState>& np1,
const nasaHeatCapacityPolynomial<equationOfState>& np2
)
{
equationOfState eofs
(
static_cast<const equationOfState&>(ct1)
+ static_cast<const equationOfState&>(ct2)
static_cast<const equationOfState&>(np1)
+ static_cast<const equationOfState&>(np2)
);
return nasaHeatCapacityPolynomial<equationOfState>
(
eofs,
ct1.nMoles()/eofs.nMoles()*ct1.a1_
+ ct2.nMoles()/eofs.nMoles()*ct2.a1_,
ct1.nMoles()/eofs.nMoles()*ct1.a2_
+ ct2.nMoles()/eofs.nMoles()*ct2.a2_,
ct1.nMoles()/eofs.nMoles()*ct1.a3_
+ ct2.nMoles()/eofs.nMoles()*ct2.a3_,
ct1.nMoles()/eofs.nMoles()*ct1.a4_
+ ct2.nMoles()/eofs.nMoles()*ct2.a4_,
ct1.nMoles()/eofs.nMoles()*ct1.a5_
+ ct2.nMoles()/eofs.nMoles()*ct2.a5_,
ct1.nMoles()/eofs.nMoles()*ct1.a6_
+ ct2.nMoles()/eofs.nMoles()*ct2.a6_,
ct1.nMoles()/eofs.nMoles()*ct1.a7_
+ ct2.nMoles()/eofs.nMoles()*ct2.a7_
np1.nMoles()/eofs.nMoles()*np1.a1_
+ np2.nMoles()/eofs.nMoles()*np2.a1_,
np1.nMoles()/eofs.nMoles()*np1.a2_
+ np2.nMoles()/eofs.nMoles()*np2.a2_,
np1.nMoles()/eofs.nMoles()*np1.a3_
+ np2.nMoles()/eofs.nMoles()*np2.a3_,
np1.nMoles()/eofs.nMoles()*np1.a4_
+ np2.nMoles()/eofs.nMoles()*np2.a4_,
np1.nMoles()/eofs.nMoles()*np1.a5_
+ np2.nMoles()/eofs.nMoles()*np2.a5_,
np1.nMoles()/eofs.nMoles()*np1.a6_
+ np2.nMoles()/eofs.nMoles()*np2.a6_,
np1.nMoles()/eofs.nMoles()*np1.a7_
+ np2.nMoles()/eofs.nMoles()*np2.a7_
);
}
@ -399,33 +399,33 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator+
template<class equationOfState>
inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator-
(
const nasaHeatCapacityPolynomial<equationOfState>& ct1,
const nasaHeatCapacityPolynomial<equationOfState>& ct2
const nasaHeatCapacityPolynomial<equationOfState>& np1,
const nasaHeatCapacityPolynomial<equationOfState>& np2
)
{
equationOfState eofs
(
static_cast<const equationOfState&>(ct1)
- static_cast<const equationOfState&>(ct2)
static_cast<const equationOfState&>(np1)
- static_cast<const equationOfState&>(np2)
);
return nasaHeatCapacityPolynomial<equationOfState>
(
eofs,
ct1.nMoles()/eofs.nMoles()*ct1.a1_
- ct2.nMoles()/eofs.nMoles()*ct2.a1_,
ct1.nMoles()/eofs.nMoles()*ct1.a2_
- ct2.nMoles()/eofs.nMoles()*ct2.a2_,
ct1.nMoles()/eofs.nMoles()*ct1.a3_
- ct2.nMoles()/eofs.nMoles()*ct2.a3_,
ct1.nMoles()/eofs.nMoles()*ct1.a4_
- ct2.nMoles()/eofs.nMoles()*ct2.a4_,
ct1.nMoles()/eofs.nMoles()*ct1.a5_
- ct2.nMoles()/eofs.nMoles()*ct2.a5_,
ct1.nMoles()/eofs.nMoles()*ct1.a6_
- ct2.nMoles()/eofs.nMoles()*ct2.a6_,
ct1.nMoles()/eofs.nMoles()*ct1.a7_
- ct2.nMoles()/eofs.nMoles()*ct2.a7_
np1.nMoles()/eofs.nMoles()*np1.a1_
- np2.nMoles()/eofs.nMoles()*np2.a1_,
np1.nMoles()/eofs.nMoles()*np1.a2_
- np2.nMoles()/eofs.nMoles()*np2.a2_,
np1.nMoles()/eofs.nMoles()*np1.a3_
- np2.nMoles()/eofs.nMoles()*np2.a3_,
np1.nMoles()/eofs.nMoles()*np1.a4_
- np2.nMoles()/eofs.nMoles()*np2.a4_,
np1.nMoles()/eofs.nMoles()*np1.a5_
- np2.nMoles()/eofs.nMoles()*np2.a5_,
np1.nMoles()/eofs.nMoles()*np1.a6_
- np2.nMoles()/eofs.nMoles()*np2.a6_,
np1.nMoles()/eofs.nMoles()*np1.a7_
- np2.nMoles()/eofs.nMoles()*np2.a7_
);
}
@ -434,19 +434,19 @@ template<class equationOfState>
inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator*
(
const scalar s,
const nasaHeatCapacityPolynomial<equationOfState>& ct
const nasaHeatCapacityPolynomial<equationOfState>& np
)
{
return nasaHeatCapacityPolynomial<equationOfState>
(
s*static_cast<const equationOfState&>(ct),
ct.a1_,
ct.a2_,
ct.a3_,
ct.a4_,
ct.a5_,
ct.a6_,
ct.a7_
s*static_cast<const equationOfState&>(np),
np.a1_,
np.a2_,
np.a3_,
np.a4_,
np.a5_,
np.a6_,
np.a7_
);
}
@ -454,11 +454,11 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator*
template<class equationOfState>
inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator==
(
const nasaHeatCapacityPolynomial<equationOfState>& ct1,
const nasaHeatCapacityPolynomial<equationOfState>& ct2
const nasaHeatCapacityPolynomial<equationOfState>& np1,
const nasaHeatCapacityPolynomial<equationOfState>& np2
)
{
return ct2 - ct1;
return np2 - np1;
}

6
src/thermophysicalModels/specie/thermo/realGasSpecieThermo/realGasSpecieThermoI.H
View file

@ -76,7 +76,7 @@ inline void Foam::realGasSpecieThermo<thermo>::T
do
{
//CL: using a stabilizing newton solver
//CL: if the solve is diverging, the time step is reduced until the solver converges
//CL: if the solve is diverging, the step is reduced until the solver converges
Tnew = Test - ((this->*F)(rho,Test) - f)/(this->*dFdT)(rho,Test)/(pow(2,i));
i++;
}while
@ -139,14 +139,14 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::gamma(const scalar rho,
template<class thermo>
inline Foam::scalar Foam::realGasSpecieThermo<thermo>::g(const scalar rho, const scalar T ) const
{
return this->h(rho, this->p(rho,T)) - T*this->s(rho, this->p(rho,T));
return this->h(rho, T) - T*this->s(rho, T);
}
template<class thermo>
inline Foam::scalar Foam::realGasSpecieThermo<thermo>::a(const scalar rho, const scalar T ) const
{
return this->e(rho,this->p(rho,T)) - T*this->s(rho,this->p(rho,T));
return this->e(rho,T ) - T*this->s(rho, T);
}