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

BUGFIX: pvdv function in Redlich Kwong model

Browse source
This commit is contained in:
Henrik Rusche 2016-08-19 10:43:17 +02:00
parent 65e6031bcb
commit 6ca6d7f302
3 changed files with 7 additions and 3 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -44,7 +44,8 @@ Foam::redlichKwong::redlichKwong(Istream& is)
Tcrit_(readScalar(is)),
a_(0.42748*pow(this->RR(), 2)*pow(Tcrit_, 2.5)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
b2_(pow(b_,2)),
b2_(b_*b_),
b3_(b2_*b_),
//CL: Only uses the default values
rhoMin_(1e-3),
rhoMax_(1500),
@ -67,7 +68,6 @@ Foam::redlichKwong::redlichKwong(const dictionary& dict)
b_(0.08664*this->RR()*Tcrit_/pcrit_),
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b5_(pow(b_,5)),
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
//CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin

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

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

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

@ -57,6 +57,7 @@ inline redlichKwong::redlichKwong(const word& name, const redlichKwong& rk)
a_(rk.a_),
b_(rk.b_),
b2_(rk.b2_),
b3_(rk.b3_),
rhoMin_(rk.rhoMin_),
rhoMax_(rk.rhoMax_),
rhostd_(rk.rhostd_)
@ -130,13 +131,15 @@ inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar T05 = pow(T, 0.5);
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return (a_*bmVm2*bpVm - this->RR()*T*T05*Vm2*bpVm2)/(T05*Vm2*bpVm2*bmVm2);
return (a_*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*T05*Vm2*bpVm2)
/(T05*Vm2*bpVm2*bmVm2);
}