From 4fd335f12a1d698fd15e09a90de10a7ac6253567 Mon Sep 17 00:00:00 2001 From: Henrik Rusche Date: Sat, 12 Feb 2011 19:41:42 +0100 Subject: [PATCH] First clean-up of RedlichKwong (with some side effects) --- .../psiThermo/realGasHThermo/realGasHThermo.C | 10 +- .../aungierRedlichKwong/aungierRedlichKwong.H | 2 +- .../aungierRedlichKwongI.H | 12 +- .../pengRobinson/pengRobinson.H | 2 +- .../pengRobinson/pengRobinsonI.H | 12 +- .../redlichKwong/redlichKwong.C | 22 +- .../redlichKwong/redlichKwong.H | 130 ++-- .../redlichKwong/redlichKwongI.H | 575 +++++++++--------- .../soaveRedlichKwong/soaveRedlichKwong.H | 2 +- .../soaveRedlichKwong/soaveRedlichKwongI.H | 12 +- .../nasaHeatCapacityPolynomial.C | 2 +- .../nasaHeatCapacityPolynomialI.H | 2 +- .../realGasConstTransportI.H | 2 +- 13 files changed, 413 insertions(+), 372 deletions(-) diff --git a/src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.C b/src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.C index a5f427949..18c6daa87 100755 --- a/src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.C +++ b/src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.C @@ -38,8 +38,6 @@ Germany template void Foam::realGasHThermo::calculate() { - - const scalarField& hCells = h_.internalField(); const scalarField& pCells = this->p_.internalField(); @@ -56,10 +54,10 @@ void Foam::realGasHThermo::calculate() const typename MixtureType::thermoType& mixture_ = this->cellMixture(celli); - mixture_.TH(hCells[celli], TCells[celli],pCells[celli],rhoCells[celli]); - psiCells[celli]=mixture_.psi(rhoCells[celli],TCells[celli]); + mixture_.TH(hCells[celli], TCells[celli], pCells[celli], rhoCells[celli]); + psiCells[celli]=mixture_.psi(rhoCells[celli], 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::realGasHThermo(const fvMesh& mesh) ( IOobject ( - "rho1", + "rho1", mesh.time().timeName(), mesh, IOobject::READ_IF_PRESENT, diff --git a/src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwong.H b/src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwong.H index de6967de3..3df68f310 100755 --- a/src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwong.H +++ b/src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwong.H @@ -105,7 +105,7 @@ public: inline scalar azentricFactor() const; inline scalar rhostd()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 inline scalar a() const; diff --git a/src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwongI.H index b13adb1bd..f4294f85a 100755 --- a/src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwongI.H +++ b/src/thermophysicalModels/specie/equationOfState/aungierRedlichKwong/aungierRedlichKwongI.H @@ -74,7 +74,7 @@ return azentricFactor_; } //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()) @@ -338,7 +338,7 @@ molarVolume=this->W()/rho0; i=0; 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)); i++; 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; for(i=0;i<200;i++) { - f1= (this->pReturn(rho1,T)-p); - f2= (this->pReturn(rho2,T)-p); + f1= (this->p(rho1,T)-p); + f2= (this->p(rho2,T)-p); 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)) { @@ -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_) { diff --git a/src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinson.H b/src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinson.H index 09013aa3a..ddd7b03ed 100755 --- a/src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinson.H +++ b/src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinson.H @@ -104,7 +104,7 @@ public: inline scalar Tcrit() const; inline scalar azentricFactor() 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 inline scalar a() const; diff --git a/src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinsonI.H b/src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinsonI.H index 45a04d45f..f2a07ac97 100755 --- a/src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinsonI.H +++ b/src/thermophysicalModels/specie/equationOfState/pengRobinson/pengRobinsonI.H @@ -70,7 +70,7 @@ return azentricFactor_; } //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)); } @@ -349,7 +349,7 @@ molarVolume=this->W()/rho0; 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)); @@ -360,10 +360,10 @@ molarVolume=this->W()/rho0; for(i=0;i<200;i++) { - f1= (this->pReturn(rho1,T)-p); - f2= (this->pReturn(rho2,T)-p); + f1= (this->p(rho1,T)-p); + f2= (this->p(rho2,T)-p); 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)) { @@ -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_) { diff --git a/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.C b/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.C index fd94b80b2..24fadc6f0 100755 --- a/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.C +++ b/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.C @@ -42,6 +42,18 @@ Germany namespace Foam { +/* * * * * * * * * * * * * * * Private static data * * * * * * * * * * * * * */ + +const scalar redlichKwong::rhoMin_ +( + debug::tolerances("redlichKwongRhoMin", 1e-3) +); + +const scalar redlichKwong::rhoMax_ +( + debug::tolerances("redlichKwongRhoMax", 1500) +); + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // redlichKwong::redlichKwong(Istream& is) @@ -49,12 +61,12 @@ redlichKwong::redlichKwong(Istream& is) specie(is), pcrit_(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)"); - 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) { os << static_cast(pg)<< tab - << pg.pcrit_ << tab<< pg.Tcrit_<< tab << pg.azentricFactor_; + << pg.pcrit_ << tab<< pg.Tcrit_; os.check("Ostream& operator<<(Ostream& os, const redlichKwong& st)"); return os; diff --git a/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.H b/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.H index 2be713d1d..3e4b8de98 100755 --- a/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.H +++ b/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwong.H @@ -46,8 +46,6 @@ Germany #include "specie.H" #include "autoPtr.H" -#include "word.H" -#include "scalar.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -61,8 +59,30 @@ namespace Foam 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) + //HR: Don't know, yet. Let's make these static for starters + static const scalar rhoMax_; + + static const scalar rhoMin_; public: @@ -71,18 +91,18 @@ public: // Constructors //- Construct from components - inline redlichKwong( - const specie& sp, - scalar pcrit, - scalar Tcrit, - scalar azentricFactor + inline redlichKwong + ( + const specie& sp, + scalar pcrit, + scalar Tcrit ); //- Construct from Istream redlichKwong(Istream&); //- Construct as named copy - inline redlichKwong(const word& name, redlichKwong&); + inline redlichKwong(const word& name, const redlichKwong&); //- Construct and return a clone inline autoPtr clone() const; @@ -90,49 +110,73 @@ public: // Selector from Istream inline static autoPtr 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 - 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 a() const; - inline scalar b() const; - + inline scalar rhostd() 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 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 rho(const scalar p,const scalar T,const scalar rho0) const; - inline scalar rho(const scalar p,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; + + // 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] inline scalar psi(const scalar rho, const scalar T) const; //- 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 diff --git a/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwongI.H index a3062d3e6..cdac92a68 100755 --- a/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwongI.H +++ b/src/thermophysicalModels/specie/equationOfState/redlichKwong/redlichKwongI.H @@ -38,245 +38,36 @@ Germany namespace Foam { - - - - -// * * * * * * * * * * * * * 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(); -} - - - - - - - - - +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Construct from components - inline redlichKwong::redlichKwong ( const specie& sp, - scalar pcrit, - scalar Tcrit, - scalar azentricFactor + scalar pcrit, + scalar Tcrit ) : - specie(sp), + specie(sp), pcrit_(pcrit), Tcrit_(Tcrit), - azentricFactor_(azentricFactor) -{ - rhostd_=this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())); // Starting GUESS for the density by ideal gas law - -} + 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()))) +{} -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - // Construct as named copy -// Starting GUESS for the density by ideal gas law -inline redlichKwong::redlichKwong(const word& name, redlichKwong& pg) +inline redlichKwong::redlichKwong(const word& name, const redlichKwong& pg) : specie(name, pg), pcrit_(pg.pcrit_), Tcrit_(pg.Tcrit_), - azentricFactor_(pg.azentricFactor_) -{ - pg.rhostd_=this->rho(Pstd,Tstd, (Pstd*this->W()/(Tstd*this->R()))); // Starting GUESS for the density by ideal gas law -} + a_(pg.a_), + b_(pg.b_), + rhostd_(pg.rhostd_) +{} // Construct and return a clone @@ -293,124 +84,320 @@ inline autoPtr redlichKwong::New(Istream& is) } + // * * * * * * * * * * * * * * * 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 -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 molarVolume; - int iter=0; + int iter=0; int maxIter_=400; scalar tol_=1e-8; - int i; - scalar rho1=rhoMax_, rho2=rhoMin_,rho3, f1,f2,f3; + scalar rho1=rhoMax_; + scalar rho2=rhoMin_; -molarVolume=this->W()/rho0; + molarVolume=this->W()/rho0; do { - molarVolumePrevIteration= molarVolume; - + molarVolumePrevIteration= molarVolume; - i=0; - do - { - molarVolume=molarVolumePrevIteration-((this->pReturn((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/ - molarVolumePrevIteration),T)))/(pow(2,i)); - - i++; - if(i>8) - { - //using bisection methode as backup, solution must be between rho=0.001 to rho=1500; - for(i=0;i<200;i++) - { - f1= (this->pReturn(rho1,T)-p); - f2= (this->pReturn(rho2,T)-p); - rho3=(rho1+rho2)/2; - f3=(this->pReturn(rho3,T)-p); + label i=0; + do + { + molarVolume=molarVolumePrevIteration + -( + (this->p((this->W()/molarVolumePrevIteration),T) - p) + /(this->dpdv((this->W()/molarVolumePrevIteration),T)) + )/pow(2,i); + + i++; + if (i>8) + { + //using bisection methode as backup, + //solution must be between rho=0.001 to rho=1500; + for(i=0; i<200; i++) + { + 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)) - { - rho1=rho3; - } - else if ((f1<0 && f3>0)||(f1>0 &&f3<0)) - { - rho2=rho3; - } - - else - { - rho2=(rho2+rho3)/2; - } + if ((f2 < 0 && f3 > 0) || (f2 > 0 && f3 < 0)) + { + rho1=rho3; + } + else if ((f1 < 0 && f3 > 0)||(f1 > 0 && f3 < 0)) + { + rho2=rho3; + } + else + { + rho2=(rho2 + rho3)/2; + } - if(mag(f3)W()/rho3; - molarVolumePrevIteration=this->W()/rho3; - break; - } - - else - { - molarVolumePrevIteration=this->W()/rho3; - } - } - } - } - while(mag(this->pReturn((this->W()/molarVolume),T)-p)>mag(this->pReturn((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(mag(f3) < p*tol_) + { + molarVolume=this->W()/rho3; + molarVolumePrevIteration=this->W()/rho3; + break; + } + else + { + molarVolumePrevIteration=this->W()/rho3; + } + } + } + } + while + ( + 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); + } } - while(mag(molarVolumePrevIteration-molarVolume)>tol_*(this->W()/rho0)); - - return this->W()/molarVolume; + while(mag(molarVolumePrevIteration-molarVolume) > tol_*(this->W()/rho0)); + return this->W()/molarVolume; } + //- 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 { - - scalar rho0=p/(this->R()*T); //using perfect gas equation as starting point - return rho(p,T,rho0); + // using perfect gas equation as starting point + return rho(p,T,p/(this->R()*T)); } + //- Return compressibility drho/dp [s^2/m^2] inline scalar redlichKwong::psi(const scalar rho, const scalar T) const { return -this->dvdp(rho,T)*pow(rho,2)/this->W(); } + //- 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 * * * * * * * * * * * * * // + inline void redlichKwong::operator+=(const redlichKwong& pg) { specie::operator+=(pg); } - inline void redlichKwong::operator-=(const redlichKwong& pg) + +inline void redlichKwong::operator-=(const redlichKwong& pg) { specie::operator-=(pg); } + inline void redlichKwong::operator*=(const scalar s) { specie::operator*=(s); diff --git a/src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwong.H b/src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwong.H index a6e8e728c..0d5f071d2 100755 --- a/src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwong.H +++ b/src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwong.H @@ -103,7 +103,7 @@ public: inline scalar Tcrit() const; inline scalar azentricFactor() 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 inline scalar a() const; diff --git a/src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwongI.H index 0cd134f43..633a18a04 100755 --- a/src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwongI.H +++ b/src/thermophysicalModels/specie/equationOfState/soaveRedlichKwong/soaveRedlichKwongI.H @@ -70,7 +70,7 @@ return azentricFactor_; } //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()))); } @@ -329,7 +329,7 @@ molarVolume=this->W()/rho0; i=0; 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)); i++; @@ -338,10 +338,10 @@ molarVolume=this->W()/rho0; //using bisection methode as backup, solution must be between rho=0.001 to rho=1500; for(i=0;i<200;i++) { - f1= (this->pReturn(rho1,T)-p); - f2= (this->pReturn(rho2,T)-p); + f1= (this->p(rho1,T)-p); + f2= (this->p(rho2,T)-p); 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)) { @@ -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_) diff --git a/src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C b/src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C index 0d7ada3d1..ddbe824fa 100755 --- a/src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C +++ b/src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C @@ -49,7 +49,7 @@ Foam::nasaHeatCapacityPolynomial::nasaHeatCapacityPolynomial(Is a7_(readScalar(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 } diff --git a/src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H b/src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H index 76ddafe9b..3f782752f 100755 --- a/src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H +++ b/src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H @@ -224,7 +224,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial::h const scalar T ) 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(); } diff --git a/src/thermophysicalModels/specie/transport/realGasConstTransport/realGasConstTransportI.H b/src/thermophysicalModels/specie/transport/realGasConstTransport/realGasConstTransportI.H index b64efcacf..f82d1ba50 100644 --- a/src/thermophysicalModels/specie/transport/realGasConstTransport/realGasConstTransportI.H +++ b/src/thermophysicalModels/specie/transport/realGasConstTransport/realGasConstTransportI.H @@ -92,7 +92,7 @@ inline scalar realGasConstTransport::alpha(const scalar rho,const scalar scalar deltaT = T - specie::Tstd; 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; }