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

BUGFIX: Switching real gas tutorial back to Peng Robinson / Some clean-up

Browse source
This commit is contained in:
Henrik Rusche 2016-07-06 23:32:35 +02:00
parent bc6dca05f1
commit 244e6aa028
14 changed files with 125 additions and 174 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C
View file

@ -56,12 +56,12 @@ Foam::aungierRedlichKwong::aungierRedlichKwong(Istream& is)
//CL: Only uses the default values //CL: Only uses the default values
rhoMin_(1e-3), rhoMin_(1e-3),
rhoMax_(1500), rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0), aSave(0.0),
daSave(0.0), daSave(0.0),
d2aSave(0.0), d2aSave(0.0),
TSave(0.0) TSave(0.0),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{ {
is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)"); is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
} }
@ -91,12 +91,12 @@ Foam::aungierRedlichKwong::aungierRedlichKwong(const dictionary& dict)
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)), rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)), rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0), aSave(0.0),
daSave(0.0), daSave(0.0),
d2aSave(0.0), d2aSave(0.0),
TSave(0.0) TSave(0.0)
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
{ {
is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)"); is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
} }

6
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H
View file

@ -87,9 +87,6 @@ class aungierRedlichKwong
scalar rhoMin_; scalar rhoMin_;
scalar rhoMax_; scalar rhoMax_;
//- Density @STD, initialise after a0, b, c, n!
scalar rhostd_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture //CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures //CL: Variables must corrected for changing temperatures
mutable scalar aSave; mutable scalar aSave;
@ -99,6 +96,9 @@ class aungierRedlichKwong
//CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct
mutable scalar TSave; mutable scalar TSave;
//- Density @STD, initialise after a0, b, c, n!
scalar rhostd_;
//Protected functions //Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave) //CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const; inline void updateModelCoefficients(const scalar T) const;

23
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H
View file

@ -68,11 +68,11 @@ inline aungierRedlichKwong::aungierRedlichKwong(const word& name, const aungierR
c2_(pg.c2_), c2_(pg.c2_),
rhoMin_(pg.rhoMin_), rhoMin_(pg.rhoMin_),
rhoMax_(pg.rhoMax_), rhoMax_(pg.rhoMax_),
rhostd_(pg.rhostd_),
aSave(0.0), aSave(0.0),
daSave(0.0), daSave(0.0),
d2aSave(0.0), d2aSave(0.0),
TSave(0.0) TSave(0.0),
rhostd_(pg.rhostd_)
{} {}
@ -237,13 +237,16 @@ inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm; scalar Vm2 = Vm*Vm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
scalar cpVm = c_ + Vm;
return return
( (
a(T)*(b3_ - 2*b2_*c_ + b_*(c_ + Vm)*(c_ - 3*Vm) + 2*Vm*pow(c_ + Vm, 2)) a(T)*(b3_ - 2*b2_*c_ + b_*cpVm*(c_ - 3*Vm) + 2*Vm*cpVm*cpVm)
- this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2) - this->RR()*T*Vm2*bpVm2
) )
/(Vm2*pow(b_ - c_ - Vm, 2)*pow(b_ + Vm, 2)); /(Vm2*pow(b_ - c_ - Vm, 2)*bpVm2);
} }
@ -284,8 +287,7 @@ inline scalar aungierRedlichKwong::integral_p_dv
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return this->RR()*T*log(-b_ + c_ + Vm) + a(T)/b_*(log(b_ + Vm) - log(Vm)); return this->RR()*T*log(Vm - b_ + c_) + a(T)/b_*log((b_ + Vm)/Vm);
//return this->RR()*T*log(-b_ + c_ + Vm) + a(T)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_;
} }
@ -299,8 +301,7 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return this->RR()*log(-b_ + c_ + Vm) + dadT(T)/b_*(log(b_ + Vm) - log(Vm)); return this->RR()*log(Vm - b_ + c_) + dadT(T)/b_*log((b_ + Vm)/Vm);
//return this->RR()*log(-b_ + c_ + Vm) + dadT(T)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_;
} }
@ -309,7 +310,7 @@ inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return -d2adT2(T)/(Vm*(Vm + b_)); return -d2adT2(T)/(Vm*(b_ + Vm));
} }
@ -400,7 +401,7 @@ inline scalar aungierRedlichKwong::integral_d2pdT2_dv
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_; return d2adT2(T)/b_*log((b_ + Vm)/Vm);
} }

10
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C
View file

@ -47,19 +47,15 @@ Foam::pengRobinson::pengRobinson(Istream& is)
b_(0.077796*this->RR()*Tcrit_/pcrit_), b_(0.077796*this->RR()*Tcrit_/pcrit_),
n_(0.37464 + 1.54226*azentricFactor_ - 0.26992*pow(azentricFactor_, 2)), n_(0.37464 + 1.54226*azentricFactor_ - 0.26992*pow(azentricFactor_, 2)),
b2_(b_*b_), b2_(b_*b_),
b3_(b2_*b_),
b4_(b3_*b_),
b5_(b4_*b_),
b6_(b5_*b_),
//CL: Only uses the default values //CL: Only uses the default values
rhoMin_(1e-3), rhoMin_(1e-3),
rhoMax_(1500), rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0), aSave(0.0),
daSave(0.0), daSave(0.0),
d2aSave(0.0), d2aSave(0.0),
TSave(0.0) TSave(0.0),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{ {
is.check("pengRobinson::pengRobinson(Istream& is)"); is.check("pengRobinson::pengRobinson(Istream& is)");
} }

10
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H
View file

@ -79,18 +79,11 @@ class pengRobinson
//CL: pow of constants b_ used in the code e.g. b2_=b*b; //CL: pow of constants b_ used in the code e.g. b2_=b*b;
scalar b2_; scalar b2_;
scalar b3_;
scalar b4_;
scalar b5_;
scalar b6_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function) //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMin_; scalar rhoMin_;
scalar rhoMax_; scalar rhoMax_;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture //CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures //CL: Variables must corrected for changing temperatures
mutable scalar aSave; mutable scalar aSave;
@ -100,6 +93,9 @@ class pengRobinson
//CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct
mutable scalar TSave; mutable scalar TSave;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//Protected functions //Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave) //CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const; inline void updateModelCoefficients(const scalar T) const;

84
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H
View file

@ -60,17 +60,13 @@ inline pengRobinson::pengRobinson(const word& name, const pengRobinson& pr)
b_(pr.b_), b_(pr.b_),
n_(pr.n_), n_(pr.n_),
b2_(pr.b2_), b2_(pr.b2_),
b3_(pr.b3_),
b4_(pr.b4_),
b5_(pr.b5_),
b6_(pr.b6_),
rhoMin_(pr.rhoMin_), rhoMin_(pr.rhoMin_),
rhoMax_(pr.rhoMax_), rhoMax_(pr.rhoMax_),
rhostd_(pr.rhostd_),
aSave(0.0), aSave(0.0),
daSave(0.0), daSave(0.0),
d2aSave(0.0), d2aSave(0.0),
TSave(0.0) TSave(0.0),
rhostd_(pr.rhostd_)
{} {}
@ -123,7 +119,7 @@ inline void pengRobinson::updateModelCoefficients(const scalar T) const
d2aSave = a0_*n_*(n_ + 1)*TTc/(2*T*T); d2aSave = a0_*n_*(n_ + 1)*TTc/(2*T*T);
//CL: saving the temperature at which the coefficients are valid //CL: saving the temperature at which the coefficients are valid
TSave=T; TSave = T;
} }
@ -131,7 +127,7 @@ inline void pengRobinson::updateModelCoefficients(const scalar T) const
inline scalar pengRobinson::a(const scalar T) const inline scalar pengRobinson::a(const scalar T) const
{ {
//CL: check if a has already been calculated for this temperature //CL: check if a has already been calculated for this temperature
if(TSave==T) if(TSave == T)
{ {
return aSave; return aSave;
} }
@ -148,7 +144,7 @@ inline scalar pengRobinson::a(const scalar T) const
inline scalar pengRobinson::dadT(const scalar T) const inline scalar pengRobinson::dadT(const scalar T) const
{ {
// check if a has already been calculated for this temperature // check if a has already been calculated for this temperature
if(TSave==T) if(TSave == T)
{ {
return daSave; return daSave;
} }
@ -165,7 +161,7 @@ inline scalar pengRobinson::dadT(const scalar T) const
inline scalar pengRobinson::d2adT2(const scalar T) const inline scalar pengRobinson::d2adT2(const scalar T) const
{ {
// check if a has already been calculated for this temperature // check if a has already been calculated for this temperature
if(TSave==T) if(TSave == T)
{ {
return d2aSave; return d2aSave;
} }
@ -210,22 +206,15 @@ inline scalar pengRobinson::p(const scalar rho, const scalar T) const
inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
return scalar bmVm = b_ - Vm;
( scalar bmVm2 = bmVm*bmVm;
2*a(T)* scalar bpVm = b_ + Vm;
(
b3_ - b2_*Vm - b_*Vm2 + Vm3 scalar t1 = bmVm2 - 2*Vm*Vm;
) scalar t12 = t1*t1;
- this->RR()*T*
( return (2*a(T)*bmVm2*bpVm - this->RR()*T*t12)/(bmVm2*t12);
b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4
)
)
/pow(b_ - Vm, 8);
} }
@ -285,7 +274,7 @@ inline scalar pengRobinson::integral_dpdT_dv
) const ) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar root2=pow(2, 0.5); scalar root2 = pow(2, 0.5);
return return
- root2*dadT(T)*log(b_*(1 - root2) + Vm)/(4*b_) - root2*dadT(T)*log(b_*(1 - root2) + Vm)/(4*b_)
@ -308,23 +297,16 @@ inline scalar pengRobinson::d2pdv2(const scalar rho, const scalar T) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm; scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
scalar Vm5 = Vm4*Vm;
scalar Vm6 = Vm5*Vm;
return 2* scalar bmVm = b_ - Vm;
( scalar bmVm2 = bmVm*bmVm;
a(T)* scalar bmVm3 = bmVm2*bmVm;
(
5*b5_ - 9*b4_*Vm + 4*b2_*Vm3 + 3*b_*Vm4 - 3*Vm5 scalar t1 = bmVm2 - 2*Vm2;
) scalar t12 = t1*t1;
- this->RR()*T* scalar t13 = t12*t1;
(
b6_ - 6*b5_*Vm + 9*b4_*Vm2 + 4*b3_*Vm3 - 9*b2_*Vm4 - 6*b_*Vm5 - Vm6 return 2*(a(T)*bmVm3*(5*b_ + 6*b_*Vm + 3*Vm2) - this->RR()*T*t13)/(bmVm*t1);
)
)
/pow(b_ - Vm, 9);
} }
@ -360,16 +342,16 @@ inline scalar pengRobinson::d2pdvdT
) const ) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
return scalar bmVm = b_ - Vm;
( scalar bmVm2 = bmVm*bmVm;
2*dadT(T)*(b3_ - b2_*Vm - b_*Vm2 + Vm3) scalar bpVm = b_ + Vm;
- this->RR()*(b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4)
) scalar t1 = bmVm2 - 2*Vm*Vm;
/pow(b_ - Vm, 6); scalar t12 = t1*t1;
scalar t14 = t12*t12;
return (2*dadT(T)*bmVm2*bpVm - this->RR()*t14)/(bmVm2*t12);
} }

2
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C
View file

@ -45,8 +45,6 @@ Foam::redlichKwong::redlichKwong(Istream& is)
a_(0.42748*pow(this->RR(), 2)*pow(Tcrit_, 2.5)/pcrit_), a_(0.42748*pow(this->RR(), 2)*pow(Tcrit_, 2.5)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_), b_(0.08664*this->RR()*Tcrit_/pcrit_),
b2_(pow(b_,2)), b2_(pow(b_,2)),
b3_(pow(b_,3)),
b5_(pow(b_,5)),
//CL: Only uses the default values //CL: Only uses the default values
rhoMin_(1e-3), rhoMin_(1e-3),
rhoMax_(1500), rhoMax_(1500),

2
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H
View file

@ -69,8 +69,6 @@ class redlichKwong
//CL: pow of constants b_ used in the code e.g. b2_=b*b; //CL: pow of constants b_ used in the code e.g. b2_=b*b;
scalar b2_; scalar b2_;
scalar b3_;
scalar b5_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function) //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMin_; scalar rhoMin_;

61
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H
View file

@ -57,8 +57,6 @@ inline redlichKwong::redlichKwong(const word& name, const redlichKwong& rk)
a_(rk.a_), a_(rk.a_),
b_(rk.b_), b_(rk.b_),
b2_(rk.b2_), b2_(rk.b2_),
b3_(rk.b3_),
b5_(rk.b5_),
rhoMin_(rk.rhoMin_), rhoMin_(rk.rhoMin_),
rhoMax_(rk.rhoMax_), rhoMax_(rk.rhoMax_),
rhostd_(rk.rhostd_) rhostd_(rk.rhostd_)
@ -132,15 +130,13 @@ inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm; scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar T05 = pow(T, 0.5); scalar T05 = pow(T, 0.5);
scalar bpVm2 = pow(b_ + Vm, 2); scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return return (a_*bmVm2*bpVm - this->RR()*T*T05*Vm2*bpVm2)/(T05*Vm2*bpVm2*bmVm2);
(
a_*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*T05*Vm2*bpVm2
)
/(T05*Vm2*bpVm2*pow(b_ - Vm, 2));
} }
@ -181,8 +177,7 @@ inline scalar redlichKwong::integral_p_dv
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return this->RR()*T*log(Vm - b_) + a_/(b_*sqrt(T))*(log(b_ + Vm) - log(Vm)); return this->RR()*T*log(Vm - b_) + a_/(b_*sqrt(T))*log((b_ + Vm)/Vm);
//return this->RR()*T*log(Vm - b_) + (a_*log(b_ + Vm))/(b_*sqrt(T)) - (a_*log(Vm))/(b_*sqrt(T));
} }
@ -196,8 +191,7 @@ inline scalar redlichKwong::integral_dpdT_dv
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return this->RR()*log(Vm - b_) - a_/(2*b_*pow(T, 1.5))*(log(b_ + Vm) - log(Vm)); return this->RR()*log(Vm - b_) - a_/(2*b_*pow(T, 1.5))*log((b_ + Vm)/Vm);
//return this->RR()*log(Vm - b_) - (a_*log(b_ + Vm))/(2*b_*T15_) + (a_*log(Vm))/(2*b_*T15_);
} }
@ -205,6 +199,7 @@ inline scalar redlichKwong::integral_dpdT_dv
inline scalar redlichKwong::d2pdT2(const scalar rho, const scalar T) const inline scalar redlichKwong::d2pdT2(const scalar rho, const scalar T) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return -0.75*a_/(pow(T, 2.5)*Vm*(b_ + Vm)); return -0.75*a_/(pow(T, 2.5)*Vm*(b_ + Vm));
} }
@ -215,19 +210,16 @@ inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm; scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm; scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
scalar Vm5 = Vm4*Vm;
scalar T05 = pow(T, 0.5); scalar T05 = pow(T, 0.5);
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bmVm3 = bmVm2*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
scalar bpVm3 = bpVm2*bpVm;
return return 2*(a_*bmVm3*(b2_ + 3*b_*Vm + 3*Vm2) + this->RR()*T*T05*Vm3*bpVm3)
( /(T05*Vm3*bpVm3*bmVm3);
2*
(
a_*(b5_ - 3*b3_*Vm2 - b2_*Vm3 + 6*b_*Vm4 - 3*Vm5)
+ this->RR()*T*T05*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
)
/(T05*Vm3*pow((b_ + Vm), 3)*pow(Vm - b_, 3))
);
} }
@ -256,18 +248,14 @@ inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm; scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar T15_ = pow(T, 1.5); scalar T15_ = pow(T, 1.5);
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return return -(0.5*a_*bmVm2*(b_ + 2*Vm) + this->RR()*T15_*Vm2*bpVm2)
-( /(T15_*Vm2*bpVm2*bmVm2);
0.5*
(
a_*(b3_ - 3*b_*Vm2 + 2*Vm3)
+ 2*this->RR()*T15_*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
)
/(T15_*Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
} }
@ -279,9 +267,10 @@ inline scalar redlichKwong::integral_d2pdT2_dv
const scalar T const scalar T
) const ) const
{ {
scalar T25_=pow(T,2.5);
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return 0.75*a_*log(b_ + Vm)/(T25_*b_) - 0.75*a_*log(Vm)/(T25_*b_); scalar T25 = pow(T, 2.5);
return 0.75*a_/(T25*b_)*log((b_ + Vm)/Vm);
} }

14
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C
View file

@ -47,17 +47,15 @@ Foam::soaveRedlichKwong::soaveRedlichKwong(Istream& is)
b_(0.08664*this->RR()*Tcrit_/pcrit_), b_(0.08664*this->RR()*Tcrit_/pcrit_),
n_(0.48508 + 1.55171*azentricFactor_ - 0.15613*pow(azentricFactor_, 2)), n_(0.48508 + 1.55171*azentricFactor_ - 0.15613*pow(azentricFactor_, 2)),
b2_(b_*b_), b2_(b_*b_),
b3_(b2_*b_),
b5_(b2_*b3_),
//CL: Only uses the default values //CL: Only uses the default values
rhoMin_(1e-3), rhoMin_(1e-3),
rhoMax_(1500), rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0), aSave(0.0),
daSave(0.0), daSave(0.0),
d2aSave(0.0), d2aSave(0.0),
TSave(0.0) TSave(0.0),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{ {
is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)"); is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)");
} }
@ -83,12 +81,12 @@ soaveRedlichKwong::soaveRedlichKwong(const dictionary& dict)
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)), rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)), rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0), aSave(0.0),
daSave(0.0), daSave(0.0),
d2aSave(0.0), d2aSave(0.0),
TSave(0.0) TSave(0.0),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{ {
is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)"); is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)");
} }

8
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H
View file

@ -76,16 +76,11 @@ class soaveRedlichKwong
//CL: pow of constants b_ used in the code e.g. b2_=b*b; //CL: pow of constants b_ used in the code e.g. b2_=b*b;
scalar b2_; scalar b2_;
scalar b3_;
scalar b5_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function) //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMin_; scalar rhoMin_;
scalar rhoMax_; scalar rhoMax_;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture //CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures //CL: Variables must corrected for changing temperatures
mutable scalar aSave; mutable scalar aSave;
@ -95,6 +90,9 @@ class soaveRedlichKwong
//CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct
mutable scalar TSave; mutable scalar TSave;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//Protected functions //Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave) //CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const; inline void updateModelCoefficients(const scalar T) const;

59
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H
View file

@ -60,15 +60,13 @@ inline soaveRedlichKwong::soaveRedlichKwong(const word& name, const soaveRedlich
b_(srk.b_), b_(srk.b_),
n_(srk.n_), n_(srk.n_),
b2_(srk.b2_), b2_(srk.b2_),
b3_(srk.b3_),
b5_(srk.b5_),
rhoMin_(srk.rhoMin_), rhoMin_(srk.rhoMin_),
rhoMax_(srk.rhoMax_), rhoMax_(srk.rhoMax_),
rhostd_(srk.rhostd_),
aSave(0.0), aSave(0.0),
daSave(0.0), daSave(0.0),
d2aSave(0.0), d2aSave(0.0),
TSave(0.0) TSave(0.0),
rhostd_(srk.rhostd_)
{} {}
@ -129,7 +127,7 @@ inline void soaveRedlichKwong::updateModelCoefficients(const scalar T) const
inline scalar soaveRedlichKwong::a(const scalar T) const inline scalar soaveRedlichKwong::a(const scalar T) const
{ {
//CL: check if a has already been calculated for this temperature //CL: check if a has already been calculated for this temperature
if(TSave==T) if(TSave == T)
{ {
return aSave; return aSave;
} }
@ -146,7 +144,7 @@ inline scalar soaveRedlichKwong::a(const scalar T) const
inline scalar soaveRedlichKwong::dadT(const scalar T) const inline scalar soaveRedlichKwong::dadT(const scalar T) const
{ {
// check if a has already been calculated for this temperature // check if a has already been calculated for this temperature
if(TSave==T) if(TSave == T)
{ {
return daSave; return daSave;
} }
@ -163,7 +161,7 @@ inline scalar soaveRedlichKwong::dadT(const scalar T) const
inline scalar soaveRedlichKwong::d2adT2(const scalar T) const inline scalar soaveRedlichKwong::d2adT2(const scalar T) const
{ {
// check if a has already been calculated for this temperature // check if a has already been calculated for this temperature
if(TSave==T) if(TSave == T)
{ {
return d2aSave; return d2aSave;
} }
@ -208,13 +206,12 @@ inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm; scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm; scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return return(a(T)*bmVm2*(b_ + 2*Vm) - this->RR()*T*Vm2*bpVm2)/(Vm2*bpVm2*bmVm2);
(
a(T)*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
/(Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
} }
@ -255,8 +252,7 @@ inline scalar soaveRedlichKwong::integral_p_dv
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return this->RR()*T*log(Vm - b_) + a(T)/b_*(log(b_ + Vm) - log(Vm)); return this->RR()*T*log(Vm - b_) + a(T)/b_*log((b_ + Vm)/Vm);
//return this->RR()*T*log(Vm - b_) + a(T)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_;
} }
@ -271,8 +267,7 @@ inline scalar soaveRedlichKwong::integral_dpdT_dv
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return this->RR()*log(Vm - b_) + dadT(T)/b_*(log(b_ + Vm) - log(Vm)); return this->RR()*log(Vm - b_) + dadT(T)/b_*log((b_ + Vm)/Vm);
//return this->RR()*log(Vm - b_) + dadT(T)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_;
} }
@ -291,15 +286,15 @@ inline scalar soaveRedlichKwong::d2pdv2(const scalar rho, const scalar T) const
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm; scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm; scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm; scalar bmVm = b_ - Vm;
scalar Vm5 = Vm4*Vm; scalar bmVm2 = bmVm*bmVm;
scalar bmVm3 = bmVm2*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
scalar bpVm3 = bpVm2*bpVm;
return 2* return 2*(a(T)*bmVm3*(b2_ + 3*b_*Vm + 3*Vm2) + this->RR()*T*Vm3*bpVm3)
( /(Vm3*bpVm3*bmVm3);
a(T)*(b5_ - 3*b3_*Vm2 - b2_*Vm3 + 6*b_*Vm4 - 3*Vm5)
+ this->RR()*T*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
)
/(Vm3*pow(b_ + Vm, 3)*pow(Vm - b_, 3));
} }
@ -336,13 +331,13 @@ inline scalar soaveRedlichKwong::d2pdvdT
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm; scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm; scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return return (dadT(T)*bmVm2*(b_ + 2*Vm) - this->RR()*Vm2*bmVm2)
( /(Vm2*bpVm2*bmVm2);
dadT(T)*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
/(Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
} }
@ -356,7 +351,7 @@ inline scalar soaveRedlichKwong::integral_d2pdT2_dv
{ {
scalar Vm = this->W()/rho; scalar Vm = this->W()/rho;
return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_; return d2adT2(T)/b_*log((b_ + Vm)/Vm);
} }

8
tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties
View file

@ -18,14 +18,14 @@ FoamFile
//CL: List of possible real gas models //CL: List of possible real gas models
thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>; //thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116; //mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>; //thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 467.6 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1.4792e-06 116; //mixture CO2 1 44.01 73.773e5 304.13 467.6 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>; thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116; mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>; //thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116; //mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;

2
tutorials/compressible/realFluidPisoFoam/ras/t-junction/system/controlDict
View file

@ -45,7 +45,7 @@ timePrecision 10;
adjustTimeStep yes; adjustTimeStep yes;
maxCo 0.25; maxCo 0.1;
maxDeltaT 1e-2; maxDeltaT 1e-2;