From 244e6aa02817e5141bb0dba9f91c8eadee6b17aa Mon Sep 17 00:00:00 2001 From: Henrik Rusche Date: Wed, 6 Jul 2016 23:32:35 +0200 Subject: [PATCH] BUGFIX: Switching real gas tutorial back to Peng Robinson / Some clean-up --- .../aungierRedlichKwong/aungierRedlichKwong.C | 10 +-- .../aungierRedlichKwong/aungierRedlichKwong.H | 6 +- .../aungierRedlichKwongI.H | 23 ++--- .../pengRobinson/pengRobinson.C | 10 +-- .../pengRobinson/pengRobinson.H | 10 +-- .../pengRobinson/pengRobinsonI.H | 84 ++++++++----------- .../redlichKwong/redlichKwong.C | 2 - .../redlichKwong/redlichKwong.H | 2 - .../redlichKwong/redlichKwongI.H | 61 ++++++-------- .../soaveRedlichKwong/soaveRedlichKwong.C | 14 ++-- .../soaveRedlichKwong/soaveRedlichKwong.H | 8 +- .../soaveRedlichKwong/soaveRedlichKwongI.H | 59 ++++++------- .../constant/thermophysicalProperties | 8 +- .../ras/t-junction/system/controlDict | 2 +- 14 files changed, 125 insertions(+), 174 deletions(-) diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C index d1164ec30..983244872 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C @@ -56,12 +56,12 @@ Foam::aungierRedlichKwong::aungierRedlichKwong(Istream& is) //CL: Only uses the default values rhoMin_(1e-3), 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), daSave(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)"); } @@ -91,12 +91,12 @@ Foam::aungierRedlichKwong::aungierRedlichKwong(const dictionary& dict) //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)), 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), daSave(0.0), d2aSave(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)"); } diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H index 472ccf69c..264dae224 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H @@ -87,9 +87,6 @@ class aungierRedlichKwong scalar rhoMin_; 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 must corrected for changing temperatures mutable scalar aSave; @@ -99,6 +96,9 @@ class aungierRedlichKwong //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct mutable scalar TSave; + //- Density @STD, initialise after a0, b, c, n! + scalar rhostd_; + //Protected functions //CL: function updates the coefficients (aSave, daSave, d2aSave) inline void updateModelCoefficients(const scalar T) const; diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H index 887717b8b..f4e37cfa6 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H @@ -68,11 +68,11 @@ inline aungierRedlichKwong::aungierRedlichKwong(const word& name, const aungierR c2_(pg.c2_), rhoMin_(pg.rhoMin_), rhoMax_(pg.rhoMax_), - rhostd_(pg.rhostd_), aSave(0.0), daSave(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 Vm2 = Vm*Vm; + scalar bpVm = b_ + Vm; + scalar bpVm2 = bpVm*bpVm; + scalar cpVm = c_ + Vm; return ( - a(T)*(b3_ - 2*b2_*c_ + b_*(c_ + Vm)*(c_ - 3*Vm) + 2*Vm*pow(c_ + Vm, 2)) - - this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2) + a(T)*(b3_ - 2*b2_*c_ + b_*cpVm*(c_ - 3*Vm) + 2*Vm*cpVm*cpVm) + - 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; - return this->RR()*T*log(-b_ + c_ + Vm) + a(T)/b_*(log(b_ + Vm) - log(Vm)); - //return this->RR()*T*log(-b_ + c_ + Vm) + a(T)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_; + return this->RR()*T*log(Vm - b_ + c_) + a(T)/b_*log((b_ + Vm)/Vm); } @@ -299,8 +301,7 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv { scalar Vm = this->W()/rho; - return this->RR()*log(-b_ + c_ + Vm) + dadT(T)/b_*(log(b_ + Vm) - log(Vm)); - //return this->RR()*log(-b_ + c_ + Vm) + dadT(T)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_; + return this->RR()*log(Vm - b_ + c_) + dadT(T)/b_*log((b_ + Vm)/Vm); } @@ -309,7 +310,7 @@ inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const { 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; - return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_; + return d2adT2(T)/b_*log((b_ + Vm)/Vm); } diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C index 1beb1cf01..719920f50 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C @@ -47,19 +47,15 @@ Foam::pengRobinson::pengRobinson(Istream& is) b_(0.077796*this->RR()*Tcrit_/pcrit_), n_(0.37464 + 1.54226*azentricFactor_ - 0.26992*pow(azentricFactor_, 2)), b2_(b_*b_), - b3_(b2_*b_), - b4_(b3_*b_), - b5_(b4_*b_), - b6_(b5_*b_), //CL: Only uses the default values rhoMin_(1e-3), 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), daSave(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)"); } diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H index 0829ae642..49c75b0d7 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H @@ -79,18 +79,11 @@ class pengRobinson //CL: pow of constants b_ used in the code e.g. b2_=b*b; 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) scalar rhoMin_; 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 must corrected for changing temperatures mutable scalar aSave; @@ -100,6 +93,9 @@ class pengRobinson //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct mutable scalar TSave; + //- Density @STD, initialise after a0, b! + scalar rhostd_; + //Protected functions //CL: function updates the coefficients (aSave, daSave, d2aSave) inline void updateModelCoefficients(const scalar T) const; diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H index 39d74d245..26147d68e 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H @@ -60,17 +60,13 @@ inline pengRobinson::pengRobinson(const word& name, const pengRobinson& pr) b_(pr.b_), n_(pr.n_), b2_(pr.b2_), - b3_(pr.b3_), - b4_(pr.b4_), - b5_(pr.b5_), - b6_(pr.b6_), rhoMin_(pr.rhoMin_), rhoMax_(pr.rhoMax_), - rhostd_(pr.rhostd_), aSave(0.0), daSave(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); //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 { //CL: check if a has already been calculated for this temperature - if(TSave==T) + if(TSave == T) { return aSave; } @@ -148,7 +144,7 @@ inline scalar pengRobinson::a(const scalar T) const inline scalar pengRobinson::dadT(const scalar T) const { // check if a has already been calculated for this temperature - if(TSave==T) + if(TSave == T) { return daSave; } @@ -165,7 +161,7 @@ inline scalar pengRobinson::dadT(const scalar T) const inline scalar pengRobinson::d2adT2(const scalar T) const { // check if a has already been calculated for this temperature - if(TSave==T) + if(TSave == T) { 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 { scalar Vm = this->W()/rho; - scalar Vm2 = Vm*Vm; - scalar Vm3 = Vm2*Vm; - scalar Vm4 = Vm3*Vm; - return - ( - 2*a(T)* - ( - b3_ - b2_*Vm - b_*Vm2 + Vm3 - ) - - this->RR()*T* - ( - b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4 - ) - ) - /pow(b_ - Vm, 8); + scalar bmVm = b_ - Vm; + scalar bmVm2 = bmVm*bmVm; + scalar bpVm = b_ + Vm; + + scalar t1 = bmVm2 - 2*Vm*Vm; + scalar t12 = t1*t1; + + return (2*a(T)*bmVm2*bpVm - this->RR()*T*t12)/(bmVm2*t12); } @@ -285,7 +274,7 @@ inline scalar pengRobinson::integral_dpdT_dv ) const { scalar Vm = this->W()/rho; - scalar root2=pow(2, 0.5); + scalar root2 = pow(2, 0.5); return - 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 Vm2 = Vm*Vm; - scalar Vm3 = Vm2*Vm; - scalar Vm4 = Vm3*Vm; - scalar Vm5 = Vm4*Vm; - scalar Vm6 = Vm5*Vm; - return 2* - ( - a(T)* - ( - 5*b5_ - 9*b4_*Vm + 4*b2_*Vm3 + 3*b_*Vm4 - 3*Vm5 - ) - - this->RR()*T* - ( - b6_ - 6*b5_*Vm + 9*b4_*Vm2 + 4*b3_*Vm3 - 9*b2_*Vm4 - 6*b_*Vm5 - Vm6 - ) - ) - /pow(b_ - Vm, 9); + scalar bmVm = b_ - Vm; + scalar bmVm2 = bmVm*bmVm; + scalar bmVm3 = bmVm2*bmVm; + + scalar t1 = bmVm2 - 2*Vm2; + scalar t12 = t1*t1; + scalar t13 = t12*t1; + + return 2*(a(T)*bmVm3*(5*b_ + 6*b_*Vm + 3*Vm2) - this->RR()*T*t13)/(bmVm*t1); } @@ -360,16 +342,16 @@ inline scalar pengRobinson::d2pdvdT ) const { scalar Vm = this->W()/rho; - scalar Vm2 = Vm*Vm; - scalar Vm3 = Vm2*Vm; - scalar Vm4 = Vm3*Vm; - return - ( - 2*dadT(T)*(b3_ - b2_*Vm - b_*Vm2 + Vm3) - - this->RR()*(b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4) - ) - /pow(b_ - Vm, 6); + scalar bmVm = b_ - Vm; + scalar bmVm2 = bmVm*bmVm; + scalar bpVm = b_ + Vm; + + scalar t1 = bmVm2 - 2*Vm*Vm; + scalar t12 = t1*t1; + scalar t14 = t12*t12; + + return (2*dadT(T)*bmVm2*bpVm - this->RR()*t14)/(bmVm2*t12); } diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C index 07f0ddc03..ea8af3415 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C @@ -45,8 +45,6 @@ Foam::redlichKwong::redlichKwong(Istream& is) a_(0.42748*pow(this->RR(), 2)*pow(Tcrit_, 2.5)/pcrit_), b_(0.08664*this->RR()*Tcrit_/pcrit_), b2_(pow(b_,2)), - b3_(pow(b_,3)), - b5_(pow(b_,5)), //CL: Only uses the default values rhoMin_(1e-3), rhoMax_(1500), diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H index 300afd1a4..73473af84 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H @@ -69,8 +69,6 @@ class redlichKwong //CL: pow of constants b_ used in the code e.g. b2_=b*b; scalar b2_; - scalar b3_; - scalar b5_; //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function) scalar rhoMin_; diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H index 9511a3733..d5d40c9a4 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H @@ -57,8 +57,6 @@ inline redlichKwong::redlichKwong(const word& name, const redlichKwong& rk) a_(rk.a_), b_(rk.b_), b2_(rk.b2_), - b3_(rk.b3_), - b5_(rk.b5_), rhoMin_(rk.rhoMin_), rhoMax_(rk.rhoMax_), rhostd_(rk.rhostd_) @@ -132,15 +130,13 @@ 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 bpVm2 = pow(b_ + Vm, 2); + scalar bmVm = b_ - Vm; + scalar bmVm2 = bmVm*bmVm; + scalar bpVm = b_ + Vm; + scalar bpVm2 = bpVm*bpVm; - return - ( - a_*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*T05*Vm2*bpVm2 - ) - /(T05*Vm2*bpVm2*pow(b_ - Vm, 2)); + return (a_*bmVm2*bpVm - this->RR()*T*T05*Vm2*bpVm2)/(T05*Vm2*bpVm2*bmVm2); } @@ -181,8 +177,7 @@ inline scalar redlichKwong::integral_p_dv { 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_*log(b_ + Vm))/(b_*sqrt(T)) - (a_*log(Vm))/(b_*sqrt(T)); + return this->RR()*T*log(Vm - b_) + a_/(b_*sqrt(T))*log((b_ + Vm)/Vm); } @@ -196,8 +191,7 @@ inline scalar redlichKwong::integral_dpdT_dv { 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_*log(b_ + Vm))/(2*b_*T15_) + (a_*log(Vm))/(2*b_*T15_); + return this->RR()*log(Vm - b_) - a_/(2*b_*pow(T, 1.5))*log((b_ + Vm)/Vm); } @@ -205,6 +199,7 @@ inline scalar redlichKwong::integral_dpdT_dv 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)); } @@ -215,19 +210,16 @@ inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const scalar Vm = this->W()/rho; scalar Vm2 = Vm*Vm; scalar Vm3 = Vm2*Vm; - scalar Vm4 = Vm3*Vm; - scalar Vm5 = Vm4*Vm; 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 - ( - 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)) - ); + return 2*(a_*bmVm3*(b2_ + 3*b_*Vm + 3*Vm2) + this->RR()*T*T05*Vm3*bpVm3) + /(T05*Vm3*bpVm3*bmVm3); } @@ -256,18 +248,14 @@ inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const { scalar Vm = this->W()/rho; scalar Vm2 = Vm*Vm; - scalar Vm3 = Vm2*Vm; scalar T15_ = pow(T, 1.5); + scalar bmVm = b_ - Vm; + scalar bmVm2 = bmVm*bmVm; + scalar bpVm = b_ + Vm; + scalar bpVm2 = bpVm*bpVm; - return - -( - 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)); + return -(0.5*a_*bmVm2*(b_ + 2*Vm) + this->RR()*T15_*Vm2*bpVm2) + /(T15_*Vm2*bpVm2*bmVm2); } @@ -279,9 +267,10 @@ inline scalar redlichKwong::integral_d2pdT2_dv const scalar T ) const { - scalar T25_=pow(T,2.5); 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); } diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C index 449744218..1b2a60056 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C @@ -47,17 +47,15 @@ Foam::soaveRedlichKwong::soaveRedlichKwong(Istream& is) b_(0.08664*this->RR()*Tcrit_/pcrit_), n_(0.48508 + 1.55171*azentricFactor_ - 0.15613*pow(azentricFactor_, 2)), b2_(b_*b_), - b3_(b2_*b_), - b5_(b2_*b3_), //CL: Only uses the default values rhoMin_(1e-3), 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), daSave(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)"); } @@ -83,12 +81,12 @@ soaveRedlichKwong::soaveRedlichKwong(const dictionary& dict) //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)), 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), daSave(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)"); } diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H index a33b8da83..adec662f1 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H @@ -76,16 +76,11 @@ class soaveRedlichKwong //CL: pow of constants b_ used in the code e.g. b2_=b*b; scalar b2_; - scalar b3_; - scalar b5_; //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function) scalar rhoMin_; 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 must corrected for changing temperatures mutable scalar aSave; @@ -95,6 +90,9 @@ class soaveRedlichKwong //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct mutable scalar TSave; + //- Density @STD, initialise after a0, b! + scalar rhostd_; + //Protected functions //CL: function updates the coefficients (aSave, daSave, d2aSave) inline void updateModelCoefficients(const scalar T) const; diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H index 9a478d085..106fd1f6f 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H @@ -60,15 +60,13 @@ inline soaveRedlichKwong::soaveRedlichKwong(const word& name, const soaveRedlich b_(srk.b_), n_(srk.n_), b2_(srk.b2_), - b3_(srk.b3_), - b5_(srk.b5_), rhoMin_(srk.rhoMin_), rhoMax_(srk.rhoMax_), - rhostd_(srk.rhostd_), aSave(0.0), daSave(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 { //CL: check if a has already been calculated for this temperature - if(TSave==T) + if(TSave == T) { return aSave; } @@ -146,7 +144,7 @@ inline scalar soaveRedlichKwong::a(const scalar T) const inline scalar soaveRedlichKwong::dadT(const scalar T) const { // check if a has already been calculated for this temperature - if(TSave==T) + if(TSave == T) { return daSave; } @@ -163,7 +161,7 @@ inline scalar soaveRedlichKwong::dadT(const scalar T) const inline scalar soaveRedlichKwong::d2adT2(const scalar T) const { // check if a has already been calculated for this temperature - if(TSave==T) + if(TSave == T) { return d2aSave; } @@ -208,13 +206,12 @@ inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const { scalar Vm = this->W()/rho; 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 - ( - 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)); + return(a(T)*bmVm2*(b_ + 2*Vm) - this->RR()*T*Vm2*bpVm2)/(Vm2*bpVm2*bmVm2); } @@ -255,8 +252,7 @@ inline scalar soaveRedlichKwong::integral_p_dv { 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)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_; + return this->RR()*T*log(Vm - b_) + a(T)/b_*log((b_ + Vm)/Vm); } @@ -271,8 +267,7 @@ inline scalar soaveRedlichKwong::integral_dpdT_dv { 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)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_; + return this->RR()*log(Vm - b_) + dadT(T)/b_*log((b_ + Vm)/Vm); } @@ -291,15 +286,15 @@ inline scalar soaveRedlichKwong::d2pdv2(const scalar rho, const scalar T) const scalar Vm = this->W()/rho; scalar Vm2 = Vm*Vm; scalar Vm3 = Vm2*Vm; - scalar Vm4 = Vm3*Vm; - scalar Vm5 = Vm4*Vm; + 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 2* - ( - 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)); + return 2*(a(T)*bmVm3*(b2_ + 3*b_*Vm + 3*Vm2) + this->RR()*T*Vm3*bpVm3) + /(Vm3*bpVm3*bmVm3); } @@ -336,13 +331,13 @@ inline scalar soaveRedlichKwong::d2pdvdT { scalar Vm = this->W()/rho; 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 - ( - 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)); + return (dadT(T)*bmVm2*(b_ + 2*Vm) - this->RR()*Vm2*bmVm2) + /(Vm2*bpVm2*bmVm2); } @@ -356,7 +351,7 @@ inline scalar soaveRedlichKwong::integral_d2pdT2_dv { 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); } diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties index 6df25c688..74073cbc4 100755 --- a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties +++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties @@ -18,14 +18,14 @@ FoamFile //CL: List of possible real gas models -thermoType realGasHThermo>>>>; -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>>>>; +//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>>>>; //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>>>>; -//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>>>>; +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>>>>; //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; diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/system/controlDict b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/system/controlDict index 9a3b44200..f3e29a18e 100644 --- a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/system/controlDict +++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/system/controlDict @@ -45,7 +45,7 @@ timePrecision 10; adjustTimeStep yes; -maxCo 0.25; +maxCo 0.1; maxDeltaT 1e-2;