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

First clean-up of RedlichKwong (with some side effects)

Browse source
This commit is contained in:
Henrik Rusche 2011-02-12 19:41:42 +01:00
parent fa742cadc1
commit 4fd335f12a
13 changed files with 413 additions and 372 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -38,8 +38,6 @@ Germany
template<class MixtureType> template<class MixtureType>
void Foam::realGasHThermo<MixtureType>::calculate() void Foam::realGasHThermo<MixtureType>::calculate()
{ {
const scalarField& hCells = h_.internalField(); const scalarField& hCells = h_.internalField();
const scalarField& pCells = this->p_.internalField(); const scalarField& pCells = this->p_.internalField();
@ -56,10 +54,10 @@ void Foam::realGasHThermo<MixtureType>::calculate()
const typename MixtureType::thermoType& mixture_ = const typename MixtureType::thermoType& mixture_ =
this->cellMixture(celli); this->cellMixture(celli);
mixture_.TH(hCells[celli], TCells[celli],pCells[celli],rhoCells[celli]); mixture_.TH(hCells[celli], TCells[celli], pCells[celli], rhoCells[celli]);
psiCells[celli]=mixture_.psi(rhoCells[celli],TCells[celli]); psiCells[celli]=mixture_.psi(rhoCells[celli], TCells[celli]);
muCells[celli] = mixture_.mu(TCells[celli]); muCells[celli] = mixture_.mu(TCells[celli]);
alphaCells[celli] = mixture_.alpha(rhoCells[celli], TCells[celli]); alphaCells[celli] = mixture_.alpha(rhoCells[celli], TCells[celli]);
} }
@ -138,7 +136,7 @@ Foam::realGasHThermo<MixtureType>::realGasHThermo(const fvMesh& mesh)
( (
IOobject IOobject
( (
"rho1", "rho1",
mesh.time().timeName(), mesh.time().timeName(),
mesh, mesh,
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,

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

@ -105,7 +105,7 @@ public:
inline scalar azentricFactor() const; inline scalar azentricFactor() const;
inline scalar rhostd()const; inline scalar rhostd()const;
inline scalar rhocrit() const; inline scalar rhocrit() const;
inline scalar pReturn(const scalar rho, const scalar T) const; inline scalar p(const scalar rho, const scalar T) const;
//-Redlich Kwong factors //-Redlich Kwong factors
inline scalar a() const; inline scalar a() const;

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

@ -74,7 +74,7 @@ return azentricFactor_;
} }
//returns the pressure for a given density and temperature //returns the pressure for a given density and temperature
inline scalar aungierRedlichKwong::pReturn(const scalar rho,const scalar T) const inline scalar aungierRedlichKwong::p(const scalar rho,const scalar T) const
{ {
return this->RR*T/((this->W()/rho)-this->b()+this->c()) return this->RR*T/((this->W()/rho)-this->b()+this->c())
@ -338,7 +338,7 @@ molarVolume=this->W()/rho0;
i=0; i=0;
do do
{ {
molarVolume=molarVolumePrevIteration-((this->pReturn((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/ molarVolume=molarVolumePrevIteration-((this->p((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/
molarVolumePrevIteration),T)))/(pow(2,i)); molarVolumePrevIteration),T)))/(pow(2,i));
i++; i++;
if(i>8) if(i>8)
@ -347,10 +347,10 @@ molarVolume=this->W()/rho0;
//using bisection methode as backup, solution must be between rho=0.001 to rho=1500; //using bisection methode as backup, solution must be between rho=0.001 to rho=1500;
for(i=0;i<200;i++) for(i=0;i<200;i++)
{ {
f1= (this->pReturn(rho1,T)-p); f1= (this->p(rho1,T)-p);
f2= (this->pReturn(rho2,T)-p); f2= (this->p(rho2,T)-p);
rho3=(rho1+rho2)/2; rho3=(rho1+rho2)/2;
f3=(this->pReturn(rho3,T)-p); f3=(this->p(rho3,T)-p);
if ((f2<0 && f3>0)||(f2>0 &&f3<0)) if ((f2<0 && f3>0)||(f2>0 &&f3<0))
{ {
@ -379,7 +379,7 @@ molarVolume=this->W()/rho0;
} }
} }
} }
while(mag(this->pReturn((this->W()/molarVolume),T)-p)>mag(this->pReturn((this->W()/molarVolumePrevIteration),T)-p)); while(mag(this->p((this->W()/molarVolume),T)-p)>mag(this->p((this->W()/molarVolumePrevIteration),T)-p));
if (iter++ > maxIter_) if (iter++ > maxIter_)
{ {

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

@ -104,7 +104,7 @@ public:
inline scalar Tcrit() const; inline scalar Tcrit() const;
inline scalar azentricFactor() const; inline scalar azentricFactor() const;
inline scalar rhostd()const; inline scalar rhostd()const;
inline scalar pReturn(const scalar rho, const scalar T) const; inline scalar p(const scalar rho, const scalar T) const;
//-Redlich Kwong factors //-Redlich Kwong factors
inline scalar a() const; inline scalar a() const;

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

@ -70,7 +70,7 @@ return azentricFactor_;
} }
//returns the pressure for a given density and temperature //returns the pressure for a given density and temperature
inline scalar pengRobinson::pReturn(const scalar rho,const scalar T) const inline scalar pengRobinson::p(const scalar rho,const scalar T) const
{ {
return this->RR*T/((this->W()/rho)-this->b())-(this->a()*pow((1+this->n()*(1-pow((T/this->Tcrit()),0.5))),2))/(pow((this->W()/rho),2)+2*this->b()*(this->W()/rho)-pow(this->b(),2)); return this->RR*T/((this->W()/rho)-this->b())-(this->a()*pow((1+this->n()*(1-pow((T/this->Tcrit()),0.5))),2))/(pow((this->W()/rho),2)+2*this->b()*(this->W()/rho)-pow(this->b(),2));
} }
@ -349,7 +349,7 @@ molarVolume=this->W()/rho0;
do do
{ {
molarVolume=molarVolumePrevIteration-((this->pReturn((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/ molarVolume=molarVolumePrevIteration-((this->p((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/
molarVolumePrevIteration),T)))/(pow(2,i)); molarVolumePrevIteration),T)))/(pow(2,i));
@ -360,10 +360,10 @@ molarVolume=this->W()/rho0;
for(i=0;i<200;i++) for(i=0;i<200;i++)
{ {
f1= (this->pReturn(rho1,T)-p); f1= (this->p(rho1,T)-p);
f2= (this->pReturn(rho2,T)-p); f2= (this->p(rho2,T)-p);
rho3=(rho1+rho2)/2; rho3=(rho1+rho2)/2;
f3=(this->pReturn(rho3,T)-p); f3=(this->p(rho3,T)-p);
if ((f2<0 && f3>0)||(f2>0 &&f3<0)) if ((f2<0 && f3>0)||(f2>0 &&f3<0))
{ {
@ -394,7 +394,7 @@ molarVolume=this->W()/rho0;
} }
while(mag(this->pReturn((this->W()/molarVolume),T)-p)>mag(this->pReturn((this->W()/molarVolumePrevIteration),T)-p)); while(mag(this->p((this->W()/molarVolume),T)-p)>mag(this->p((this->W()/molarVolumePrevIteration),T)-p));
if (iter++ > maxIter_) if (iter++ > maxIter_)
{ {

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

@ -42,6 +42,18 @@ Germany
namespace Foam namespace Foam
{ {
/* * * * * * * * * * * * * * * Private static data * * * * * * * * * * * * * */
const scalar redlichKwong::rhoMin_
(
debug::tolerances("redlichKwongRhoMin", 1e-3)
);
const scalar redlichKwong::rhoMax_
(
debug::tolerances("redlichKwongRhoMax", 1500)
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
redlichKwong::redlichKwong(Istream& is) redlichKwong::redlichKwong(Istream& is)
@ -49,12 +61,12 @@ redlichKwong::redlichKwong(Istream& is)
specie(is), specie(is),
pcrit_(readScalar(is)), pcrit_(readScalar(is)),
Tcrit_(readScalar(is)), Tcrit_(readScalar(is)),
azentricFactor_(readScalar(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())))
{ {
is.check("redlichKwong::redlichKwong(Istream& is)"); is.check("redlichKwong::redlichKwong(Istream& is)");
rhostd_=this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R()));
rhoMax_=1500;
rhoMin_=0.001;
} }
@ -63,7 +75,7 @@ redlichKwong::redlichKwong(Istream& is)
Ostream& operator<<(Ostream& os, const redlichKwong& pg) Ostream& operator<<(Ostream& os, const redlichKwong& pg)
{ {
os << static_cast<const specie&>(pg)<< tab os << static_cast<const specie&>(pg)<< tab
<< pg.pcrit_ << tab<< pg.Tcrit_<< tab << pg.azentricFactor_; << pg.pcrit_ << tab<< pg.Tcrit_;
os.check("Ostream& operator<<(Ostream& os, const redlichKwong& st)"); os.check("Ostream& operator<<(Ostream& os, const redlichKwong& st)");
return os; return os;

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

@ -46,8 +46,6 @@ Germany
#include "specie.H" #include "specie.H"
#include "autoPtr.H" #include "autoPtr.H"
#include "word.H"
#include "scalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -61,8 +59,30 @@ namespace Foam
class redlichKwong class redlichKwong
: :
public specie 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)
//HR: Don't know, yet. Let's make these static for starters
static const scalar rhoMax_;
static const scalar rhoMin_;
public: public:
@ -71,18 +91,18 @@ public:
// Constructors // Constructors
//- Construct from components //- Construct from components
inline redlichKwong( inline redlichKwong
const specie& sp, (
scalar pcrit, const specie& sp,
scalar Tcrit, scalar pcrit,
scalar azentricFactor scalar Tcrit
); );
//- Construct from Istream //- Construct from Istream
redlichKwong(Istream&); redlichKwong(Istream&);
//- Construct as named copy //- Construct as named copy
inline redlichKwong(const word& name, redlichKwong&); inline redlichKwong(const word& name, const redlichKwong&);
//- Construct and return a clone //- Construct and return a clone
inline autoPtr<redlichKwong> clone() const; inline autoPtr<redlichKwong> clone() const;
@ -90,49 +110,73 @@ public:
// Selector from Istream // Selector from Istream
inline static autoPtr<redlichKwong> New(Istream& is); inline static autoPtr<redlichKwong> New(Istream& is);
//Member Variabels
scalar pcrit_;
scalar Tcrit_;
scalar rhostd_;
scalar azentricFactor_;
scalar rhoMax_; //should be read from the fvSolution file where rhoMax and rhoMin values must be define ( for rhoSimpleFoam)
scalar rhoMin_;
// Member functions // Member functions
inline scalar pcrit() const;
inline scalar Tcrit() const;
inline scalar rhostd()const;
inline scalar azentricFactor() const;
inline scalar pReturn(const scalar rho, const scalar T) const;
//-Redlich Kwong factors inline scalar rhostd() const;
inline scalar a() const;
inline scalar b() const; inline scalar p(const scalar rho, const scalar T) const;
// Derivatives
inline scalar dpdv(const scalar rho, const scalar T) const;
//derivatives
inline scalar dpdv(const scalar rho,const scalar T) const;
inline scalar dpdT(const scalar rho, const scalar T) const; inline scalar dpdT(const scalar rho, const scalar T) const;
inline scalar dvdT(const scalar rho,const scalar T) const;
inline scalar dvdp(const scalar rho, const scalar T) const;
inline scalar isobarExpCoef(const scalar rho,const scalar T) const;
inline scalar isothermalCompressiblity(const scalar rho,const scalar T) const;
inline scalar integral_d2pdT2_dv(const scalar rho,const scalar T) const ; // Used for cv
inline scalar d2pdv2(const scalar rho,const scalar T) const; // not Used At The Moment
inline scalar d2pdT2(const scalar rho,const scalar T) const; // not Used At The Moment
inline scalar d2pdvdT(const scalar rho,const scalar T) const; // not Used At The Moment
inline scalar d2vdT2(const scalar rho,const scalar T) const; // not Used At The Moment
inline scalar integral_p_dv(const scalar rho,const scalar T) const; //Used for internal Energy
inline scalar integral_dpdT_dv(const scalar rho,const scalar T) const; //Used for Entropy
//- Return density [kg/m^3] // rho0 is the starting point of the newton solver used to calculate rho inline scalar dvdT(const scalar rho, const scalar T) const;
inline scalar rho(const scalar p,const scalar T,const scalar rho0) const;
inline scalar rho(const scalar p,const scalar T) const; inline scalar dvdp(const scalar rho, const scalar T) const;
inline scalar isobarExpCoef(const scalar rho, const scalar T) const;
inline scalar isothermalCompressiblity
(
const scalar rho,
const scalar T
) const;
// Used for cv
inline scalar integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const;
//Used for internal Energy
inline scalar integral_p_dv(const scalar rho, const scalar T) const;
// Used for Entropy
inline scalar integral_dpdT_dv(const scalar rho, const scalar T) const;
// 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;
inline scalar d2pdvdT(const scalar rho, const scalar T) const;
inline scalar d2vdT2(const scalar rho, const scalar T) const;
//- Return density [kg/m^3]
// rho0 is the starting point of the newton solver used to calculate rho
inline scalar rho
(
const scalar p,
const scalar T,
const scalar rho0
) const;
inline scalar rho(const scalar p, const scalar T) const;
//- Return compressibility drho/dp [s^2/m^2] //- Return compressibility drho/dp [s^2/m^2]
inline scalar psi(const scalar rho, const scalar T) const; inline scalar psi(const scalar rho, const scalar T) const;
//- Return compression factor [] //- Return compression factor []
inline scalar Z(const scalar p,const scalar T,const scalar rho0) const; inline scalar Z
(
const scalar p,
const scalar T,
const scalar rho0
) const;
// Member operators // Member operators

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

@ -38,245 +38,36 @@ Germany
namespace Foam namespace Foam
{ {
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline scalar redlichKwong::pcrit()const
{
return pcrit_;
}
inline scalar redlichKwong::Tcrit()const
{
return Tcrit_;
}
// Returns the Azentric Factor (Acentric Factor)
inline scalar redlichKwong::azentricFactor() const
{
return azentricFactor_;
}
inline scalar redlichKwong::rhostd()const
{
return rhostd_;
}
//returns the pressure for a given density and temperature
inline scalar redlichKwong::pReturn(const scalar rho,const scalar T) const
{
return (this->RR *T/((this->W()/rho)-this->b()) - (this->a()/(sqrt(T)*(this->W()/rho)*((this->W()/rho)+this->b()))));
}
// Faktor a of the redlich Kwong equation of state
//(molar values)
inline scalar redlichKwong::a() const
{
return 0.42748*pow(this->RR,2)*pow(this->Tcrit(),2.5)/(this->pcrit());
}
// Faktor b of the redlich Kwong equation of state
//(molar values)
inline scalar redlichKwong::b() const
{
return 0.08664*this->RR*this->Tcrit()/this->pcrit();
}
//* * * * * * * * * * * * * Derivatives * * * * * * * * * * * //
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
{
return (this->a()*(pow(this->b(),3)-3*this->b()*pow((this->W()/rho),2)+2*pow((this->W()/rho),3))-this->RR*pow(T,1.5)*pow((this->W()/rho),2) *(pow(this->b(),2)+2*this->b()*(this->W()/rho)+pow((this->W()/rho),2))) / ( sqrt(T) * pow((this->W()/rho),2) * pow((this->b()+(this->W()/rho)),2) *pow( (this->b()-(this->W()/rho)),2) );
}
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar redlichKwong::dpdT(const scalar rho, const scalar T) const
{
return 0.5*this->a()/(pow(T,1.5)*(this->W()/rho)*(this->b()+(this->W()/rho)))-this->RR/(this->b()-(this->W()/rho));
}
//Real deviative dv/dT at constant pressure
//using implicit differentiation
//(molar values)
inline scalar redlichKwong::dvdT(const scalar rho,const scalar T) const
{
return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
}
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar redlichKwong::dvdp(const scalar rho,const scalar T) const
{
return 1/this->dpdv(rho,T);
}
//needed to calculate the internal energy
//(molar values)
inline scalar redlichKwong::integral_p_dv(const scalar rho,const scalar T) const
{
return this->RR*T*log((this->W()/rho)-this->b())+(this->a()*log(this->b()+(this->W()/rho)))/(this->b()*sqrt(T))-(this->a()*log(this->W()/rho))/(this->b()*sqrt(T));
}
//needed to calculate the entropy
//(molar values)
inline scalar redlichKwong::integral_dpdT_dv(const scalar rho,const scalar T) const
{
return this->RR*log((this->W()/rho)-this->b())
-(this->a()*log(this->b()+(this->W()/rho)))/(2*this->b()* pow(T,1.5) )
+(this->a()*log(this->W()/rho))/(2*this->b()* pow(T,1.5));
}
//* * * * * * * * * * * * * second order Derivative based functions * * * * * * * * * * * //
//(molar values)
inline scalar redlichKwong::d2pdT2(const scalar rho,const scalar T) const
{
return -0.75*this->a()/( pow(T,2.5)*(this->W()/rho)*(this->b()+(this->W()/rho)));
}
//(molar values)
inline scalar redlichKwong::d2pdv2(const scalar rho,const scalar T) const
{
return
(
2*(this->a()*(pow(this->b(),5)-3*pow(this->b(),3)*pow((this->W()/rho),2)-pow(this->b(),2)*pow((this->W()/rho),3)+6*this->b()*pow((this->W()/rho),4)-3*pow((this->W()/rho),5))+this->RR*pow(T,1.5)*pow((this->W()/rho),3)
*(pow(this->b(),3)+3*pow(this->b(),2)*(this->W()/rho)+3*this->b()*pow((this->W()/rho),2)+pow((this->W()/rho),3)))/( sqrt(T)* pow((this->W()/rho),3) *pow( (this->b()+(this->W()/rho)) ,3)* pow((this->W()/rho)-this->b(),3)));
}
//(molar values)
//using second Order implicit differentiation
inline scalar redlichKwong::d2vdT2(const scalar rho, const scalar T) const
{
return -(pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)+ pow(this->dpdv(rho,T),2) *this->d2pdT2(rho,T)- 2*this->dpdv(rho,T)*this->dpdT(rho,T)*this->d2pdvdT(rho,T))/( pow(this->dpdv(rho,T),3));
}
//(molar values)
inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const
{
return -(0.5*(this->a()*( pow(this->b(),3)-3*this->b()* pow((this->W()/rho),2) +2* pow((this->W()/rho),3) )+2*this->RR*pow(T,1.5)* pow((this->W()/rho),2)*( pow(this->b(),2) + 2*this->b()*(this->W()/rho)+ pow((this->W()/rho),2) )))/( pow(T,1.5)* pow((this->W()/rho),2)*pow(this->b()+(this->W()/rho),2)*pow(this->b()-(this->W()/rho),2));
}
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar redlichKwong::integral_d2pdT2_dv(const scalar rho,const scalar T) const
{
return 0.75*this->a()*log(this->b() + (this->W()/rho))/(pow(T,2.5)*this->b()) - 0.75*this->a()*log((this->W()/rho))/(pow(T,2.5)*this->b());
}
//* * * * * * * * * * * * * thermodynamic properties * * * * * * * * * * * //
//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
//(molar values)
inline scalar redlichKwong::isobarExpCoef(const scalar rho,const scalar T) const
{
return this->dvdT(rho, T)*rho/this->W();
}
//isothemal compressiblity kappa
//(molar values)
inline scalar redlichKwong::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();
}
// Construct from components // Construct from components
inline redlichKwong::redlichKwong inline redlichKwong::redlichKwong
( (
const specie& sp, const specie& sp,
scalar pcrit, scalar pcrit,
scalar Tcrit, scalar Tcrit
scalar azentricFactor
) )
: :
specie(sp), specie(sp),
pcrit_(pcrit), pcrit_(pcrit),
Tcrit_(Tcrit), Tcrit_(Tcrit),
azentricFactor_(azentricFactor) a_(0.42748*pow(this->RR,2)*pow(Tcrit_,2.5)/pcrit_),
{ b_(0.08664*this->RR*Tcrit_/pcrit_),
rhostd_=this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())); // Starting GUESS for the density by ideal gas law // Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R())))
} {}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct as named copy // Construct as named copy
// Starting GUESS for the density by ideal gas law inline redlichKwong::redlichKwong(const word& name, const redlichKwong& pg)
inline redlichKwong::redlichKwong(const word& name, redlichKwong& pg)
: :
specie(name, pg), specie(name, pg),
pcrit_(pg.pcrit_), pcrit_(pg.pcrit_),
Tcrit_(pg.Tcrit_), Tcrit_(pg.Tcrit_),
azentricFactor_(pg.azentricFactor_) a_(pg.a_),
{ b_(pg.b_),
pg.rhostd_=this->rho(Pstd,Tstd, (Pstd*this->W()/(Tstd*this->R()))); // Starting GUESS for the density by ideal gas law rhostd_(pg.rhostd_)
} {}
// Construct and return a clone // Construct and return a clone
@ -293,124 +84,320 @@ inline autoPtr<redlichKwong> redlichKwong::New(Istream& is)
} }
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline scalar redlichKwong::rhostd()const
{
return rhostd_;
}
//returns the pressure for a given density and temperature
inline scalar redlichKwong::p(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR*T/(Vm - b_) - a_/(sqrt(T)*Vm*(Vm + b_));
}
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return (a_*(pow(b_,3) - 3*b_*pow(Vm,2) + 2*pow(Vm,3))
- this->RR*pow(T,1.5)*pow(Vm,2)*(pow(b_,2) + 2*b_*Vm + pow(Vm,2)))
/(sqrt(T)*pow(Vm,2)*pow((b_ + Vm),2)*pow( (b_ - Vm),2));
}
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar redlichKwong::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return 0.5*a_/(pow(T,1.5)*Vm*(b_ + Vm))-this->RR/(b_ - Vm);
}
//Real deviative dv/dT at constant pressure
//using implicit differentiation
//(molar values)
inline scalar redlichKwong::dvdT(const scalar rho, const scalar T) const
{
return -this->dpdT(rho,T)/this->dpdv(rho,T);
}
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar redlichKwong::dvdp(const scalar rho, const scalar T) const
{
return 1/this->dpdv(rho,T);
}
//needed to calculate the internal energy
//(molar values)
inline scalar redlichKwong::integral_p_dv
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
return this->RR*T*log(Vm - b_)
+ (a_*log(b_ + Vm))/(b_*sqrt(T))
- (a_*log(Vm))/(b_*sqrt(T));
}
//needed to calculate the entropy
//(molar values)
inline scalar redlichKwong::integral_dpdT_dv
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
return this->RR*log(Vm - b_)
-(a_*log(b_ + Vm))/(2*b_*pow(T,1.5))
+(a_*log(Vm))/(2*b_*pow(T,1.5));
}
//(molar values)
inline scalar redlichKwong::d2pdT2(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return -0.75*a_/(pow(T,2.5)*Vm*(b_ + Vm));
}
//(molar values)
inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return
(
2*(
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)
)
+ this->RR*pow(T,1.5)*pow(Vm,3)*(pow(b_,3)
+ 3*pow(b_,2)*Vm
+ 3*b_*pow(Vm,2)+pow(Vm,3))
)
/(sqrt(T)*pow(Vm,3)*pow((b_ + Vm),3)*pow(Vm-b_,3))
);
}
//(molar values)
//using second Order implicit differentiation
inline scalar redlichKwong::d2vdT2(const scalar rho, const scalar T) const
{
return
-(
pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+ pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
- 2*this->dpdv(rho,T)*this->dpdT(rho,T)*this->d2pdvdT(rho,T)
)
/(pow(this->dpdv(rho,T),3));
}
//(molar values)
inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return
-(0.5*(
a_*(pow(b_,3) - 3*b_*pow(Vm,2) + 2*pow(Vm,3))
+ 2*this->RR*pow(T,1.5)*pow(Vm,2)*(pow(b_,2) + 2*b_*Vm + pow(Vm,2))
))
/(pow(T,1.5)*pow(Vm,2)*pow(b_ + Vm,2)*pow(b_ - Vm,2));
}
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar redlichKwong::integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
return 0.75*a_*log(b_ + Vm)/(pow(T,2.5)*b_)
- 0.75*a_*log(Vm)/(pow(T,2.5)*b_);
}
//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
//(molar values)
inline scalar redlichKwong::isobarExpCoef
(
const scalar rho,
const scalar T
) const
{
return this->dvdT(rho, T)*rho/this->W();
}
//isothemal compressiblity kappa
//(molar values)
inline scalar redlichKwong::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();
}
//- Return density [kg/m^3]on //- Return density [kg/m^3]on
inline scalar redlichKwong::rho(const scalar p,const scalar T,const scalar rho0) const inline scalar redlichKwong::rho
(
const scalar p,
const scalar T,
const scalar rho0
) const
{ {
scalar molarVolumePrevIteration; scalar molarVolumePrevIteration;
scalar molarVolume; scalar molarVolume;
int iter=0; int iter=0;
int maxIter_=400; int maxIter_=400;
scalar tol_=1e-8; scalar tol_=1e-8;
int i; scalar rho1=rhoMax_;
scalar rho1=rhoMax_, rho2=rhoMin_,rho3, f1,f2,f3; scalar rho2=rhoMin_;
molarVolume=this->W()/rho0; molarVolume=this->W()/rho0;
do do
{ {
molarVolumePrevIteration= molarVolume; molarVolumePrevIteration= molarVolume;
i=0; label i=0;
do do
{ {
molarVolume=molarVolumePrevIteration-((this->pReturn((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/ molarVolume=molarVolumePrevIteration
molarVolumePrevIteration),T)))/(pow(2,i)); -(
(this->p((this->W()/molarVolumePrevIteration),T) - p)
i++; /(this->dpdv((this->W()/molarVolumePrevIteration),T))
if(i>8) )/pow(2,i);
{
//using bisection methode as backup, solution must be between rho=0.001 to rho=1500; i++;
for(i=0;i<200;i++) if (i>8)
{ {
f1= (this->pReturn(rho1,T)-p); //using bisection methode as backup,
f2= (this->pReturn(rho2,T)-p); //solution must be between rho=0.001 to rho=1500;
rho3=(rho1+rho2)/2; for(i=0; i<200; i++)
f3=(this->pReturn(rho3,T)-p); {
scalar f1 = this->p(rho1,T) - p;
scalar f2 = this->p(rho2,T) - p;
scalar rho3 = (rho1 + rho2)/2;
scalar f3 = this->p(rho3,T) - p;
if ((f2<0 && f3>0)||(f2>0 &&f3<0)) if ((f2 < 0 && f3 > 0) || (f2 > 0 && f3 < 0))
{ {
rho1=rho3; rho1=rho3;
} }
else if ((f1<0 && f3>0)||(f1>0 &&f3<0)) else if ((f1 < 0 && f3 > 0)||(f1 > 0 && f3 < 0))
{ {
rho2=rho3; rho2=rho3;
} }
else
else {
{ rho2=(rho2 + rho3)/2;
rho2=(rho2+rho3)/2; }
}
if(mag(f3)<p*tol_) if(mag(f3) < p*tol_)
{ {
molarVolume=this->W()/rho3; molarVolume=this->W()/rho3;
molarVolumePrevIteration=this->W()/rho3; molarVolumePrevIteration=this->W()/rho3;
break; break;
} }
else
else {
{ molarVolumePrevIteration=this->W()/rho3;
molarVolumePrevIteration=this->W()/rho3; }
} }
} }
} }
} while
while(mag(this->pReturn((this->W()/molarVolume),T)-p)>mag(this->pReturn((this->W()/molarVolumePrevIteration),T)-p)); (
mag(this->p((this->W()/molarVolume),T) - p)
> mag(this->p((this->W()/molarVolumePrevIteration),T) - p)
if (iter++ > maxIter_) );
{
FatalErrorIn
(
"inline scalar redlichKwong::rho(const scalar p,const scalar T,const scalar rho0) const "
) << "Maximum number of iterations exceeded"
<< abort(FatalError);
}
if (iter++ > maxIter_)
{
FatalErrorIn
(
"inline scalar redlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
) << "Maximum number of iterations exceeded"
<< abort(FatalError);
}
} }
while(mag(molarVolumePrevIteration-molarVolume)>tol_*(this->W()/rho0)); while(mag(molarVolumePrevIteration-molarVolume) > tol_*(this->W()/rho0));
return this->W()/molarVolume;
return this->W()/molarVolume;
} }
//- Return density [kg/m^3]on //- Return density [kg/m^3]on
inline scalar redlichKwong::rho(const scalar p,const scalar T) const inline scalar redlichKwong::rho(const scalar p, const scalar T) const
{ {
// using perfect gas equation as starting point
scalar rho0=p/(this->R()*T); //using perfect gas equation as starting point return rho(p,T,p/(this->R()*T));
return rho(p,T,rho0);
} }
//- Return compressibility drho/dp [s^2/m^2] //- Return compressibility drho/dp [s^2/m^2]
inline scalar redlichKwong::psi(const scalar rho, const scalar T) const inline scalar redlichKwong::psi(const scalar rho, const scalar T) const
{ {
return -this->dvdp(rho,T)*pow(rho,2)/this->W(); return -this->dvdp(rho,T)*pow(rho,2)/this->W();
} }
//- Return compression factor [] //- Return compression factor []
inline scalar redlichKwong::Z( const scalar p, const scalar T,const scalar rho0) const inline scalar redlichKwong::Z
(
const scalar p,
const scalar T,
const scalar rho0
) const
{ {
return (p*this->rho(p,T,rho0))/(this->R()*T); return p*this->rho(p,T,rho0)/(this->R()*T);
} }
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void redlichKwong::operator+=(const redlichKwong& pg) inline void redlichKwong::operator+=(const redlichKwong& pg)
{ {
specie::operator+=(pg); specie::operator+=(pg);
} }
inline void redlichKwong::operator-=(const redlichKwong& pg)
inline void redlichKwong::operator-=(const redlichKwong& pg)
{ {
specie::operator-=(pg); specie::operator-=(pg);
} }
inline void redlichKwong::operator*=(const scalar s) inline void redlichKwong::operator*=(const scalar s)
{ {
specie::operator*=(s); specie::operator*=(s);

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

@ -103,7 +103,7 @@ public:
inline scalar Tcrit() const; inline scalar Tcrit() const;
inline scalar azentricFactor() const; inline scalar azentricFactor() const;
inline scalar rhostd()const; inline scalar rhostd()const;
inline scalar pReturn(const scalar rho, const scalar T) const; inline scalar p(const scalar rho, const scalar T) const;
//-Redlich Kwong factors //-Redlich Kwong factors
inline scalar a() const; inline scalar a() const;

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

@ -70,7 +70,7 @@ return azentricFactor_;
} }
//returns the pressure for a given density and temperature //returns the pressure for a given density and temperature
inline scalar soaveRedlichKwong::pReturn(const scalar rho,const scalar T) const inline scalar soaveRedlichKwong::p(const scalar rho,const scalar T) const
{ {
return (this->RR*T/((this->W()/rho)-this->b())-this->a()*pow((1+this->n()*(1-pow((T/this->Tcrit()),0.5))),2)/((this->W()/rho)*((this->W()/rho)+this->b()))); return (this->RR*T/((this->W()/rho)-this->b())-this->a()*pow((1+this->n()*(1-pow((T/this->Tcrit()),0.5))),2)/((this->W()/rho)*((this->W()/rho)+this->b())));
} }
@ -329,7 +329,7 @@ molarVolume=this->W()/rho0;
i=0; i=0;
do do
{ {
molarVolume=molarVolumePrevIteration-((this->pReturn((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/ molarVolume=molarVolumePrevIteration-((this->p((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/
molarVolumePrevIteration),T)))/(pow(2,i)); molarVolumePrevIteration),T)))/(pow(2,i));
i++; i++;
@ -338,10 +338,10 @@ molarVolume=this->W()/rho0;
//using bisection methode as backup, solution must be between rho=0.001 to rho=1500; //using bisection methode as backup, solution must be between rho=0.001 to rho=1500;
for(i=0;i<200;i++) for(i=0;i<200;i++)
{ {
f1= (this->pReturn(rho1,T)-p); f1= (this->p(rho1,T)-p);
f2= (this->pReturn(rho2,T)-p); f2= (this->p(rho2,T)-p);
rho3=(rho1+rho2)/2; rho3=(rho1+rho2)/2;
f3=(this->pReturn(rho3,T)-p); f3=(this->p(rho3,T)-p);
if ((f2<0 && f3>0)||(f2>0 &&f3<0)) if ((f2<0 && f3>0)||(f2>0 &&f3<0))
{ {
@ -371,7 +371,7 @@ molarVolume=this->W()/rho0;
} }
} }
} }
while(mag(this->pReturn((this->W()/molarVolume),T)-p)>mag(this->pReturn((this->W()/molarVolumePrevIteration),T)-p)); while(mag(this->p((this->W()/molarVolume),T)-p)>mag(this->p((this->W()/molarVolumePrevIteration),T)-p));
if (iter++ > maxIter_) if (iter++ > maxIter_)

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

@ -49,7 +49,7 @@ Foam::nasaHeatCapacityPolynomial<equationOfState>::nasaHeatCapacityPolynomial(Is
a7_(readScalar(is)) a7_(readScalar(is))
{ {
is.check("nasaHeatCapacityPolynomial::nasaHeatCapacityPolynomial(Istream& is)"); is.check("nasaHeatCapacityPolynomial::nasaHeatCapacityPolynomial(Istream& is)");
cp_std=this->cp_nonLimited(this->rhostd_,this->Tstd); // cp @ STD (needed to limit cp for stability cp_std=this->cp_nonLimited(this->rhostd(),this->Tstd); // cp @ STD (needed to limit cp for stability
} }

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

@ -224,7 +224,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::h
const scalar T const scalar T
) const ) const
{ {
return this->e(rho,T)+this->pReturn(rho,T)/rho*this->W()-this->Pstd/this->rhostd()*this->W(); return this->e(rho,T)+this->p(rho,T)/rho*this->W()-this->Pstd/this->rhostd()*this->W();
} }

2
src/thermophysicalModels/specie/transport/realGasConstTransport/realGasConstTransportI.H
View file

@ -92,7 +92,7 @@ inline scalar realGasConstTransport<thermo>::alpha(const scalar rho,const scalar
scalar deltaT = T - specie::Tstd; scalar deltaT = T - specie::Tstd;
scalar CpBar = scalar CpBar =
(deltaT*(this->H(rho,T) - this->H(this->rhostd_,specie::Tstd)) + Cp_)/(sqr(deltaT) + 1); (deltaT*(this->H(rho,T) - this->H(this->rhostd(),specie::Tstd)) + Cp_)/(sqr(deltaT) + 1);
return Cp_*mu(T)*rPr/CpBar; return Cp_*mu(T)*rPr/CpBar;
} }