From e0f776a6544d230f3d955d5447701d5863d7c38f Mon Sep 17 00:00:00 2001 From: Christian Lucas Date: Wed, 22 Jun 2016 20:34:48 +0200 Subject: [PATCH] corrections in real gas library --- .../IAPWSThermo/IAPWS-IF97.C | 20 +- .../aungierRedlichKwong/aungierRedlichKwong.C | 2 +- .../aungierRedlichKwong/aungierRedlichKwong.H | 6 +- .../aungierRedlichKwongI.H | 33 +-- .../pengRobinson/pengRobinson.H | 6 +- .../pengRobinson/pengRobinsonI.H | 28 ++- .../redlichKwong/redlichKwongI.H | 22 +- .../soaveRedlichKwong/soaveRedlichKwong.C | 5 - .../soaveRedlichKwong/soaveRedlichKwong.H | 6 +- .../soaveRedlichKwong/soaveRedlichKwongI.H | 24 +- .../constantHeatCapacity.C | 5 +- .../constantHeatCapacity.H | 2 +- .../constantHeatCapacityI.H | 24 +- .../nasaHeatCapacityPolynomial.C | 5 +- .../nasaHeatCapacityPolynomialI.H | 76 +++--- .../realGasSpecieThermoI.H | 4 +- .../constRealGas/constRealGasTransport.C | 56 ----- .../constRealGas/constRealGasTransport.H | 202 ---------------- .../constRealGas/constRealGasTransportI.H | 217 ------------------ .../constant/thermophysicalProperties | 35 ++- .../ras/backStep_IAPWS97/system/controlDict | 7 - .../ras/cavity_IAPWS97/0.org/omega | 47 ---- .../ras/cavity_IAPWS97/0/omega | 47 ---- .../realFluidPisoFoam/ras/t-junction/0.org/T | 48 ++++ .../realFluidPisoFoam/ras/t-junction/0.org/U | 47 ++++ .../ras/t-junction/0.org/epsilon | 58 +++++ .../realFluidPisoFoam/ras/t-junction/0.org/k | 52 +++++ .../realFluidPisoFoam/ras/t-junction/0.org/p | 46 ++++ .../constant/thermophysicalProperties | 39 ++-- 29 files changed, 439 insertions(+), 730 deletions(-) delete mode 100644 src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.C delete mode 100644 src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.H delete mode 100644 src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransportI.H delete mode 100644 tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/omega delete mode 100644 tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/omega create mode 100755 tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/T create mode 100644 tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/U create mode 100644 tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/epsilon create mode 100644 tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/k create mode 100644 tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/p diff --git a/src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.C b/src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.C index 415a8e555..38c7c32dc 100644 --- a/src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.C +++ b/src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.C @@ -144,6 +144,7 @@ void Foam::calculateProperties_h region=freesteam_region(S); + //CL:Liquid phase if (region==1) { p=S.R1.p; @@ -169,6 +170,7 @@ void Foam::calculateProperties_h lambda=freesteam_k_rhoT(rho,T); alpha=lambda/cp; //Cl: Important info -->alpha= thermal diffusivity time density } + //CL:vapor phase else if (region==2) { p=S.R2.p; @@ -194,6 +196,7 @@ void Foam::calculateProperties_h lambda=freesteam_k_rhoT(rho,T); alpha=lambda/cp; //Cl: Important info -->alpha= thermal diffusivity time density } + //CL: supercritial fluid else if (region==3) { scalar gamma,cv; @@ -235,9 +238,10 @@ void Foam::calculateProperties_h alpha=lambda/cp; //Cl: Important info -->alpha= thermal diffusivity time density } + //inside the vapor dome else if (region==4) { - scalar rhov,rhol,betav,betal,kappav,kappal,vv,vl,cpl,hl,hv,cp; + scalar rhov,rhol,betav,betal,kappav,kappal,vv,vl,cpl,cpv,hl,hv,cp; scalar dvldp,dvvdp,dhldp,dhvdp; scalar dpdT,dvdh,dvdp,dxdp; @@ -272,7 +276,7 @@ void Foam::calculateProperties_h betav=freesteam_region2_alphav_pT(Sv.R2.p,Sv.R2.T); cpl=freesteam_region1_cp_pT(Sl.R1.p,Sl.R1.T); - //cpv=freesteam_region2_cp_pT(Sv.R2.p,Sv.R2.T); + cpv=freesteam_region2_cp_pT(Sv.R2.p,Sv.R2.T); hl=freesteam_region1_h_pT(Sl.R1.p,Sl.R1.T); hv=freesteam_region2_h_pT(Sv.R2.p,Sv.R2.T); @@ -283,7 +287,7 @@ void Foam::calculateProperties_h dvvdp=betav*vv/dpdT-kappav*vv; dhldp=vl*(1-betal*Sl.R1.T)+cpl/dpdT; - dhvdp=vv*(1-betav*Sv.R2.T)+cpl/dpdT; + dhvdp=vv*(1-betav*Sv.R2.T)+cpv/dpdT; dxdp=-dhldp/(hv-hl) +(h-hl)/((hv-hl)*(hv-hl)) @@ -386,6 +390,7 @@ Foam::scalar Foam::psiH(SteamState S) region=freesteam_region(S); + //CL:liquid phase if (region==1) { //Cl: note: in FreeStream, beta=1/V*(dV/dP)_P=const is called alphaV (in this region) @@ -399,6 +404,7 @@ Foam::scalar Foam::psiH(SteamState S) //CL: psiH=(drho/dp)_h=const psiH=-((S.R1.T*beta*beta-beta)/cp-kappa*rho); } + //CL:vapor phase else if (region==2) { //Cl: note: in FreeStream, beta=1/V*(dV/dP)_P=const is called alphaV (in this region) @@ -412,6 +418,7 @@ Foam::scalar Foam::psiH(SteamState S) //CL: psiH=(drho/dp)_h=const psiH=-((S.R2.T*beta*beta-beta)/cp-kappa*rho); } + //CL:supercritical fluid else if (region==3) { @@ -434,9 +441,10 @@ Foam::scalar Foam::psiH(SteamState S) psiH=-((S.R3.T*beta*beta-beta)/cp-kappa*rho); } + //inside the vapor dome else if (region==4) { - scalar rhov,rhol,betav,betal,kappav,kappal,vv,vl,cpl,hl,hv,h,p; + scalar rhov,rhol,betav,betal,kappav,kappal,vv,vl,cpl,cpv,hl,hv,h,p; scalar dvldp,dvvdp,dhldp,dhvdp; scalar dpdT,dvdp,dxdp; @@ -467,7 +475,7 @@ Foam::scalar Foam::psiH(SteamState S) betav=freesteam_region2_alphav_pT(Sv.R2.p,Sv.R2.T); cpl=freesteam_region1_cp_pT(Sl.R1.p,Sl.R1.T); - //cpv=freesteam_region2_cp_pT(Sv.R2.p,Sv.R2.T); + cpv=freesteam_region2_cp_pT(Sv.R2.p,Sv.R2.T); hl=freesteam_region1_h_pT(Sl.R1.p,Sl.R1.T); hv=freesteam_region2_h_pT(Sv.R2.p,Sv.R2.T); @@ -477,7 +485,7 @@ Foam::scalar Foam::psiH(SteamState S) dvvdp=betav*vv/dpdT-kappav*vv; dhldp=vl*(1-betal*Sl.R1.T)+cpl/dpdT; - dhvdp=vv*(1-betav*Sv.R2.T)+cpl/dpdT; + dhvdp=vv*(1-betav*Sv.R2.T)+cpv/dpdT; dxdp=-dhldp/(hv-hl) +(h-hl)/((hv-hl)*(hv-hl)) diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C index 2c4716ce6..1fd9fe926 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C @@ -127,7 +127,7 @@ void Foam::aungierRedlichKwong::write(Ostream& os) const Foam::Ostream& Foam::operator<<(Ostream& os, const aungierRedlichKwong& ark) { os << static_cast(ark)<< token::SPACE - << ark.pcrit_ << tab<< ark.Tcrit_<< tab<W()/rho; + return this->RR()*T/(Vm - b_ + c_) - a(T)/(Vm*(Vm + b_)); } @@ -238,10 +241,10 @@ inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const 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_*(c_ + Vm)*(c_ - 3*Vm) + 2*Vm*(c_ + Vm)*(c_ + Vm)) - this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2) ) - /(Vm2*pow(b_ - c_ - Vm, 2)*pow(b_ + Vm, 2)); + /(Vm2*(b_ - c_ - Vm)*(b_ - c_ - Vm)*(b_ + Vm)*(b_ + Vm)); } @@ -320,12 +323,12 @@ inline scalar aungierRedlichKwong::d2pdv2(const scalar rho,const scalar T) const ( b5_ - 3*b4_*c_ + 3*b3_*(c2_ - c_*Vm - Vm2) - b2_*(c_ + Vm)*(c2_ - 7*c_*Vm+Vm2) - + 3*b_*Vm*pow(c_ + Vm, 2)*(2*Vm - c_) - - 3*Vm2*pow(c_ + Vm, 3) + + 3*b_*Vm*(c_ + Vm)*(c_ + Vm)*(2*Vm - c_) + - 3*Vm2*(c_ + Vm)*(c_ + Vm)*(c_ + Vm) ) + this->RR()*T*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3) ) - /(Vm3*pow(b_ - c_ - Vm, 3)*pow(b_ + Vm, 3)); + /(Vm3*(b_ - c_ - Vm)*(b_ - c_ - Vm)*(b_ - c_ - Vm)*(b_ + Vm)*(b_ + Vm)*(b_ + Vm)); } @@ -337,15 +340,17 @@ inline scalar aungierRedlichKwong::d2vdT2 const scalar T ) const { - scalar dpdT2=this->dpdT(rho, T)*this->dpdT(rho, T); - scalar dpdv2=this->dpdv(rho, T)*this->dpdv(rho, T); - scalar dpdv3=dpdv2*this->dpdv(rho, T); + scalar dpdT_=this->dpdT(rho, T); + scalar dpdT2=dpdT_*dpdT_; + scalar dpdv_=this->dpdv(rho, T); + scalar dpdv2=dpdv_*dpdv_; + scalar dpdv3=dpdv2*dpdv_; return -( dpdT2*this->d2pdv2(rho, T) + dpdv2*this->d2pdT2(rho, T) - - 2*this->dpdv(rho, T)*this->dpdT(rho, T)*this->d2pdvdT(rho, T) + - 2*dpdv_*dpdT_*this->d2pdvdT(rho, T) ) /dpdv3; } @@ -363,10 +368,10 @@ inline scalar aungierRedlichKwong::d2pdvdT return ( - dadT(T)*(b3_ - 2*b2_*c_ + b_*(c_ + Vm)*(c_ - 3*Vm)+2*Vm*pow(c_ + Vm, 2)) + dadT(T)*(b3_ - 2*b2_*c_ + b_*(c_ + Vm)*(c_ - 3*Vm)+2*Vm*(c_ + Vm)*(c_ + Vm)) - this->RR()*Vm2*(b2_ + 2*b_*Vm + Vm2) ) - /(Vm2*pow(b_ - c_ - Vm, 2)*pow(b_ + Vm, 2)); + /(Vm2*(b_ - c_ - Vm)*(b_ - c_ - Vm)*(b_ + Vm)*(b_ + Vm)); } diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H index 90739c075..e23614451 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H @@ -97,7 +97,7 @@ class pengRobinson mutable scalar daSave; mutable scalar d2aSave; - //CL: save the temperature for which the save coefficients (amix,dadTmix,d2adT2mix) are correct + //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct mutable scalar TSave; //Protected functions @@ -149,10 +149,10 @@ public: //CL: Model coefficient a(T) inline scalar a(const scalar T) const; - //CL: temperature deriviative of model coefficient a(T) + //CL: temperature derivative of model coefficient a(T) inline scalar dadT(const scalar T) const; - //CL: second order temperature deriviative of model coefficient a(T) + //CL: second order temperature derivative of model coefficient a(T) inline scalar d2adT2(const scalar T) const; //CL: Equation of state diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H index 92ec2f30b..e0dcc9a12 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H @@ -117,9 +117,11 @@ inline scalar pengRobinson::Tcrit() const inline void pengRobinson::updateModelCoefficients(const scalar T) const { - aSave=a0_*pow(1 + n_*(1 - pow(T/Tcrit_, 0.5)), 2); - daSave=a0_*n_*(n_*sqrt(T/Tcrit_) - n_ - 1)*sqrt(T/Tcrit_)/T; - d2aSave=a0_*n_*(n_ + 1)*sqrt(T/Tcrit_)/(2*T*T); + scalar TTc_=sqrt(T/Tcrit_); + + aSave=a0_*(1 + n_*(1 - TTc_))*(1 + n_*(1 - TTc_)); + daSave=a0_*n_*(n_*TTc_ - n_ - 1)*TTc_/T; + d2aSave=a0_*n_*(n_ + 1)*TTc_/(2*T*T); //CL: saving the temperature at which the coefficients are valid TSave=T; @@ -224,7 +226,7 @@ inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4 ) ) - /(pow(b_ - Vm, 2)*pow(b2_ - 2*b_*Vm - Vm2, 2)); + /((b_ - Vm)*(b_ - Vm)*(b2_ - 2*b_*Vm - Vm2)*(b2_ - 2*b_*Vm - Vm2)); } @@ -320,7 +322,11 @@ inline scalar pengRobinson::d2pdv2(const scalar rho, const scalar T) const b6_ - 6*b5_*Vm + 9*b4_*Vm2 + 4*b3_*Vm3 - 9*b2_*Vm4 - 6*b_*Vm5 - Vm6 ) ) - /(pow(b_ - Vm, 3)*pow(b2_ - 2*b_*Vm - Vm2, 3)); + / + ( + (b_ - Vm)*(b_ - Vm)*(b_ - Vm)*(b2_ - 2*b_*Vm - Vm2) + *(b2_ - 2*b_*Vm - Vm2)*(b2_ - 2*b_*Vm - Vm2) + ); } @@ -332,15 +338,17 @@ inline scalar pengRobinson::d2vdT2 const scalar T ) const { - scalar dpdT2=this->dpdT(rho, T)*this->dpdT(rho, T); - scalar dpdv2=this->dpdv(rho, T)*this->dpdv(rho, T); - scalar dpdv3=dpdv2*this->dpdv(rho, T); + scalar dpdT_=this->dpdT(rho, T); + scalar dpdT2=dpdT_*dpdT_; + scalar dpdv_=this->dpdv(rho, T); + scalar dpdv2=dpdv_*dpdv_; + scalar dpdv3=dpdv2*dpdv_; return -( dpdT2*this->d2pdv2(rho, T) + dpdv2*this->d2pdT2(rho, T) - - 2*this->dpdv(rho, T)*this->dpdT(rho, T)*this->d2pdvdT(rho, T) + - 2*dpdv_*dpdT_*this->d2pdvdT(rho, T) ) /dpdv3; } @@ -363,7 +371,7 @@ inline scalar pengRobinson::d2pdvdT 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, 2)*pow(b2_ - 2*b_*Vm - Vm2, 2)); + /((b_ - Vm)*(b_ - Vm)*(b2_ - 2*b_*Vm - Vm2)*(b2_ - 2*b_*Vm - Vm2)); } diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H index 253b1fb00..9cf9c1721 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H @@ -133,12 +133,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); return ( - a_*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*pow(T, 1.5)*Vm2*(b2_ + 2*b_*Vm + Vm2) + a_*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*T05_*Vm2*(b2_ + 2*b_*Vm + Vm2) ) - /(sqrt(T)*Vm2*pow((b_ + Vm), 2)*pow( (b_ - Vm), 2)); + /(T05_*Vm2*(b_ + Vm)*(b_ + Vm)*(b_ - Vm)*(b_ - Vm)); } @@ -212,15 +213,16 @@ inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const scalar Vm3 = Vm2*Vm; scalar Vm4 = Vm3*Vm; scalar Vm5 = Vm4*Vm; + scalar T05_ = pow(T, 0.5); return ( 2* ( a_*(b5_ - 3*b3_*Vm2 - b2_*Vm3 + 6*b_*Vm4 - 3*Vm5) - + this->RR()*pow(T, 1.5)*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3) + + this->RR()*T*T05_*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3) ) - /(sqrt(T)*Vm3*pow((b_ + Vm), 3)*pow(Vm - b_, 3)) + /(T05_*Vm3*(b_ + Vm)*(b_ + Vm)*(b_ + Vm)*(Vm - b_)*(Vm - b_)*(Vm - b_)) ); } @@ -229,15 +231,17 @@ inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const //using second order implicit differentiation inline scalar redlichKwong::d2vdT2(const scalar rho, const scalar T) const { - scalar dpdT2=this->dpdT(rho, T)*this->dpdT(rho, T); - scalar dpdv2=this->dpdv(rho, T)*this->dpdv(rho, T); - scalar dpdv3=dpdv2*this->dpdv(rho, T); + scalar dpdT_=this->dpdT(rho, T); + scalar dpdT2=dpdT_*dpdT_; + scalar dpdv_=this->dpdv(rho, T); + scalar dpdv2=dpdv_*dpdv_; + scalar dpdv3=dpdv2*dpdv_; return -( dpdT2*this->d2pdv2(rho, T) + dpdv2*this->d2pdT2(rho, T) - - 2*this->dpdv(rho, T)*this->dpdT(rho, T)*this->d2pdvdT(rho, T) + - 2*dpdv_*dpdT_*this->d2pdvdT(rho, T) ) /dpdv3; } @@ -259,7 +263,7 @@ inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const + 2*this->RR()*T15_*Vm2*(b2_ + 2*b_*Vm + Vm2) ) ) - /(T15_*Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2)); + /(T15_*Vm2*(b_ + Vm)*(b_ + Vm)*(b_ - Vm)*(b_ - Vm)); } diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C index 1ad68fba6..3b86beb24 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C @@ -84,11 +84,6 @@ 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)), - //CL: rhoMin and rhoMax are only used as boundaries for the bisection method (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 - 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), diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H index 4114abfb0..d4b246a06 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H @@ -93,7 +93,7 @@ class soaveRedlichKwong mutable scalar daSave; mutable scalar d2aSave; - //CL: save the temperature for which the save coefficients (amix,dadTmix,d2adT2mix) are correct + //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct mutable scalar TSave; //Protected functions @@ -145,10 +145,10 @@ public: //CL: Model coefficient a(T) inline scalar a(const scalar T) const; - //CL: temperature deriviative of model coefficient a(T) + //CL: temperature derivative of model coefficient a(T) inline scalar dadT(const scalar T) const; - //CL: second order temperature deriviative of model coefficient a(T) + //CL: second order temperature derivative of model coefficient a(T) inline scalar d2adT2(const scalar T) const; //CL: Equation of state diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H index 58bd76efc..82c369d5a 100755 --- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H +++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H @@ -115,9 +115,11 @@ inline scalar soaveRedlichKwong::Tcrit() const inline void soaveRedlichKwong::updateModelCoefficients(const scalar T) const { - aSave=a0_*pow(1 + n_*(1 - pow(T/Tcrit_, 0.5)), 2); - daSave=a0_*n_*(n_*sqrt(T/Tcrit_) - n_ - 1)*sqrt(T/Tcrit_)/T; - d2aSave=a0_*n_*(n_ + 1)*sqrt(T/Tcrit_)/(2*T*T); + scalar TTc_=pow(T/Tcrit_, 0.5); + + aSave=a0_*(1 + n_*(1 - TTc_))*(1 + n_*(1 - TTc_)); + daSave=a0_*n_*(n_*TTc_ - n_ - 1)*TTc_/T; + d2aSave=a0_*n_*(n_ + 1)*TTc_/(2*T*T); //CL: saving the temperature at which the coefficients are valid TSave=T; @@ -213,7 +215,7 @@ inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const ( 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)); + /(Vm2*(b_ + Vm)*(b_ + Vm)*(b_ - Vm)*(b_ - Vm)); } @@ -294,7 +296,7 @@ inline scalar soaveRedlichKwong::d2pdv2(const scalar rho, const scalar T) const 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)); + /(Vm3*(b_ + Vm)*(b_ + Vm)*(b_ + Vm)*(Vm - b_)*(Vm - b_)*(Vm - b_)); } @@ -306,15 +308,17 @@ inline scalar soaveRedlichKwong::d2vdT2 const scalar T ) const { - scalar dpdT2=this->dpdT(rho, T)*this->dpdT(rho, T); - scalar dpdv2=this->dpdv(rho, T)*this->dpdv(rho, T); - scalar dpdv3=dpdv2*this->dpdv(rho, T); + scalar dpdT_=this->dpdT(rho, T); + scalar dpdT2=dpdT_*dpdT_; + scalar dpdv_=this->dpdv(rho, T); + scalar dpdv2=dpdv_*dpdv_; + scalar dpdv3=dpdv2*dpdv_; return -( dpdT2*this->d2pdv2(rho, T) + dpdv2*this->d2pdT2(rho, T) - - 2*this->dpdv(rho, T)*this->dpdT(rho, T)*this->d2pdvdT(rho, T) + - 2*dpdv_*dpdT_*this->d2pdvdT(rho, T) ) /dpdv3; } @@ -335,7 +339,7 @@ inline scalar soaveRedlichKwong::d2pdvdT ( 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)); + /(Vm2*(b_ + Vm)*(b_ + Vm)*(b_ - Vm)*(b_ - Vm)); } diff --git a/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.C b/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.C index a9d178362..b8daa981e 100755 --- a/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.C +++ b/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.C @@ -38,14 +38,15 @@ template Foam::constantHeatCapacity::constantHeatCapacity(Istream& is) : equationOfState(is), + //CL: reads specific perfect gas heat capacity Cp0_(readScalar(is)), cp0_(Cp0_*this->W()), - //values for some need terms at std + //CL: values for some need terms at std e0_std(e0(this->Tstd())), s0_std(s0(this->Tstd())), integral_p_dv_std(this->integral_p_dv(this->rhostd(),this->Tstd())), integral_dpdT_dv_std(this->integral_dpdT_dv(this->rhostd(),this->Tstd())), - // cp @ STD (needed to limit cp for stability + //CL: cp @ STD (needed to limit cp for stability cp_std(this->cp_nonLimited(this->rhostd(),this->Tstd())) { is.check("constantHeatCapacity::constantHeatCapacity(Istream& is)"); diff --git a/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.H b/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.H index 438b386f2..4748e486a 100755 --- a/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.H +++ b/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.H @@ -108,7 +108,7 @@ class constantHeatCapacity { // Private data - //CL: spec. cp + //CL: specific values scalar Cp0_; //CL: molar values scalar cp0_; diff --git a/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacityI.H b/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacityI.H index f851be3ef..46d8837ac 100755 --- a/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacityI.H +++ b/src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacityI.H @@ -183,10 +183,15 @@ inline Foam::scalar Foam::constantHeatCapacity::cp const scalar T ) const { - // Problem --> dpdv(rho,T) is =0 at some points within the vapour dome. To increase stability, (dp/dv) has to be limited - // cp can be negative within the vapor dome. To avoid this nonphysical result, the absolute value is used. - // within the vapourdome and at the critical point, cp increases to very high values --> infinity, - // this would decrease the stability, so cp will be limited to 20 times the cp @ STD + scalar dpdT_=this->dpdT(rho, T); + + // CL: within the vapor dome, some stabilty problems might occur due to unphysical results (EOS not valite in vapor dome) + // CL: In order to find the right solution, the newton solver might go into the vapor dome (in an intermediate step) + // CL: Therefore, to avoid a crash, the following stabilisation step should be used + // CL: dpdv(rho,T) is =0 at some points within the vapour dome. To increase stability, (dp/dv) has to be limited + // CL: cp can be negative within the vapor dome. To avoid this nonphysical result, the absolute value is used. + // CL: within the vapourdome and at the critical point, cp increases to very high values --> infinity, + // CL: this would decrease the stability, so cp will be limited to 20 times the cp @ STD return min ( @@ -194,7 +199,7 @@ inline Foam::scalar Foam::constantHeatCapacity::cp fabs ( this->cv(rho,T) - - T*pow((this->dpdT(rho, T)), 2)/min(this->dpdv(rho, T), -1) + - T*dpdT_*dpdT_/min(this->dpdv(rho, T), -1) ) ); } @@ -209,9 +214,11 @@ inline Foam::scalar Foam::constantHeatCapacity::cp_nonLimited const scalar T ) const { + scalar dpdT_=this->dpdT(rho, T); + return fabs(this->cv(rho, T) - - T*pow((this->dpdT(rho, T)), 2)/min(this->dpdv(rho, T), -1)); + - T*dpdT_*dpdT_/min(this->dpdv(rho, T), -1)); } @@ -278,9 +285,8 @@ inline Foam::scalar Foam::constantHeatCapacity::s { return - integral_dpdT_dv_std - + this->s0(T) - - s0_std - - this->RR*log(T/this->Tstd()) + + this->s0(T) - s0_std + - this->RR()*log(T/this->Tstd()) + this->integral_dpdT_dv(rho, T); } diff --git a/src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C b/src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C index 20af3b65d..2654ded65 100755 --- a/src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C +++ b/src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C @@ -46,12 +46,13 @@ Foam::nasaHeatCapacityPolynomial::nasaHeatCapacityPolynomial(Is a5_(readScalar(is)), a6_(readScalar(is)), a7_(readScalar(is)), - //values for some need terms at std + //CL: values at std needed to calculate enthalpy and entropy + //CL: enthalpy and entropy = 0 @ std e0_std(e0(this->Tstd())), s0_std(s0(this->Tstd())), integral_p_dv_std(this->integral_p_dv(this->rhostd(),this->Tstd())), integral_dpdT_dv_std(this->integral_dpdT_dv(this->rhostd(),this->Tstd())), - //cp @ STD (needed to limit cp for stability + //cp @ STD (needed to limit cp for stability) cp_std(this->cp_nonLimited(this->rhostd(),this->Tstd())) { is.check("nasaHeatCapacityPolynomial::nasaHeatCapacityPolynomial(Istream& is)"); diff --git a/src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H b/src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H index dbdd3016c..11f204e8b 100755 --- a/src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H +++ b/src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H @@ -158,15 +158,15 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial::h0 ) const { return - this->RR()*T* + this->RR()* ( - -this->a1_*pow(T,-2) - +this->a2_*log(T)/T - +this->a3_ - +0.5*this->a4_*T - +(this->a5_*pow(T,2))/3 - +(this->a6_*pow(T,3))/4 - +(this->a7_*pow(T,4))/5 + - this->a1_/T + + this->a2_*log(T) + + T*(this->a3_ + + T*(0.5*this->a4_ + + T*(this->a5_/3 + + T*(0.25*this->a6_ + + T*(0.2*this->a7_))))) ); } @@ -191,12 +191,14 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial::s0 { return this->RR()* ( - this->a1_*(-1)/(2*pow(T,2)) - -this->a2_/T+this->a3_*log(T) - +this->a4_*T - +(this->a5_*pow(T,2))/2 - +(this->a6_*pow(T,3))/3 - +(this->a7_*pow(T,4))/4 + 1/T*( + - 0.5*this->a1_/T + - this->a2_) + + this->a3_*log(T) + + T*(this->a4_ + + T*(0.5*this->a5_ + + T*(this->a6_/3 + + T*(0.25*this->a7_)))) ); } @@ -210,13 +212,14 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial::cp0 { return this->RR()* ( - this->a1_*1/pow(T,2) - +this->a2_*1/T - +this->a3_ - +this->a4_*T - +this->a5_*pow(T,2) - +this->a6_*pow(T,3) - +this->a7_*pow(T,4) + 1/T*( + this->a1_*1/T + + this->a2_) + + this->a3_ + + T*(this->a4_ + + T*(this->a5_ + + T*(this->a6_ + + T*(this->a7_)))) ); } @@ -242,18 +245,22 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial::cp const scalar T ) const { - // Problem --> dpdv(rho,T) is =0 at some points within the vapour dome. To increase stability, (dp/dv) has to be limited - // cp can be negative within the vapor dome. To avoid this nonphysical result, the absolute value is used. - // within the vapourdome and at the critical point, cp increases to very high values --> infinity, - // this would decrease the stability, so cp will be limited to 20 times the cp @ STD + scalar dpdT_=this->dpdT(rho, T); + + // CL: within the vapor dome, some stabilty problems might occur due to unphysical results (EOS not valite in vapor dome) + // CL: In order to find the right solution, the newton solver might go into the vapor dome (in an intermediate step) + // CL: Therefore, to avoid a crash, the following stabilisation step should be used + // CL: dpdv(rho,T) is =0 at some points within the vapour dome. To increase stability, (dp/dv) has to be limited + // CL: cp can be negative within the vapor dome. To avoid this nonphysical result, the absolute value is used. + // CL: within the vapourdome and at the critical point, cp increases to very high values --> infinity, + // CL: this would decrease the stability, so cp will be limited to 20 times the cp @ STD return min ( cp_std*20, fabs ( - this->cv(rho,T) - -T*pow((this->dpdT(rho, T)),2) + this->cv(rho,T) - T*dpdT_*dpdT_ /min(this->dpdv(rho, T),-1) ) ); @@ -268,7 +275,9 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial::cp_nonLim const scalar T ) const { - return fabs(this->cv(rho,T)-T*pow((this->dpdT(rho, T)),2)/min(this->dpdv(rho, T),-1)); + scalar dpdT_=this->dpdT(rho, T); + + return fabs(this->cv(rho,T)-T*dpdT_*dpdT_/min(this->dpdv(rho, T),-1)); } @@ -329,10 +338,13 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial::s const scalar T ) const { - return -integral_dpdT_dv_std - +(this->s0(T)-s0_std) - -this->RR()*log(T/this->Tstd()) - + this->integral_dpdT_dv(rho,T); + return + ( + -integral_dpdT_dv_std + +this->s0(T)-s0_std + -this->RR()*log(T/this->Tstd()) + + this->integral_dpdT_dv(rho,T) + ); } diff --git a/src/thermophysicalModels/specie/thermo/realGasThermo/realGasSpecieThermo/realGasSpecieThermoI.H b/src/thermophysicalModels/specie/thermo/realGasThermo/realGasSpecieThermo/realGasSpecieThermoI.H index c8834608e..37fe57a89 100755 --- a/src/thermophysicalModels/specie/thermo/realGasThermo/realGasSpecieThermo/realGasSpecieThermoI.H +++ b/src/thermophysicalModels/specie/thermo/realGasThermo/realGasSpecieThermo/realGasSpecieThermoI.H @@ -234,13 +234,13 @@ inline Foam::scalar Foam::realGasSpecieThermo::psiE return -( ( - T*pow(beta,2)*V + T*beta*beta*V -this->isothermalCompressiblity(rho,T)*cp ) / ( cp*V - -beta*this->p(rho,T)*pow(V,2) + -beta*this->p(rho,T)*V*V ) ); } diff --git a/src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.C b/src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.C deleted file mode 100644 index 1892116c0..000000000 --- a/src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.C +++ /dev/null @@ -1,56 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | - \\ / A nd | For copyright notice see file Copyright - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "constRealGasTransport.H" -#include "IOstreams.H" - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -Foam::constRealGasTransport::constRealGasTransport(Istream& is) -: - Thermo(is), - mu_(readScalar(is)), - rPr_(1.0/readScalar(is)) -{ - is.check("constRealGasTransport::constRealGasTransport(Istream& is)"); -} - - -// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // - -template -Foam::Ostream& Foam::operator<<(Ostream& os, const constRealGasTransport& ct) -{ - operator<<(os, static_cast(ct)); - os << tab << ct.mu_ << tab << 1.0/ct.rPr_; - - os.check("Ostream& operator<<(Ostream&, const constRealGasTransport&)"); - - return os; -} - - -// ************************************************************************* // diff --git a/src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.H b/src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.H deleted file mode 100644 index 7b73a4df1..000000000 --- a/src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.H +++ /dev/null @@ -1,202 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | - \\ / A nd | For copyright notice see file Copyright - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - Foam::constRealGasTransport - -Description - Constant properties Transport package. - Templated into a given thermodynamics package (needed for thermal - conductivity). - -SourceFiles - constRealGasTransportI.H - constRealGasTransport.C - -\*---------------------------------------------------------------------------*/ - -#ifndef constRealGasTransport_H -#define constRealGasTransport_H - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of friend functions and operators - -template class constRealGasTransport; - -template -inline constRealGasTransport operator+ -( - const constRealGasTransport&, - const constRealGasTransport& -); - -template -inline constRealGasTransport operator- -( - const constRealGasTransport&, - const constRealGasTransport& -); - -template -inline constRealGasTransport operator* -( - const scalar, - const constRealGasTransport& -); - -template -inline constRealGasTransport operator== -( - const constRealGasTransport&, - const constRealGasTransport& -); - -template -Ostream& operator<< -( - Ostream&, - const constRealGasTransport& -); - - -/*---------------------------------------------------------------------------*\ - Class constRealGasTransport Declaration -\*---------------------------------------------------------------------------*/ - -template -class constRealGasTransport -: - public Thermo -{ - // Private data - - //- Constant dynamic viscosity [Pa.s] - scalar mu_; - - //- Reciprocal Prandtl Number [] - scalar rPr_; - - - // Private Member Functions - - //- Construct from components - inline constRealGasTransport - ( - const Thermo& t, - const scalar mu, - const scalar Pr - ); - - -public: - - // Constructors - - //- Construct as named copy - inline constRealGasTransport(const word&, const constRealGasTransport&); - - //- Construct from Istream - constRealGasTransport(Istream&); - - - // Member functions - - //- Dynamic viscosity [kg/ms] - inline scalar mu(const scalar T) const; - - //- Thermal conductivity [W/mK] - inline scalar kappa(const scalar rho, const scalar T) const; - - //- Thermal diffusivity for enthalpy [kg/ms] - inline scalar alpha(const scalar rho, const scalar T) const; - - // Species diffusivity - //inline scalar D(const scalar T) const; - - - // Member operators - - inline constRealGasTransport& operator= - ( - const constRealGasTransport& - ); - - - // Friend operators - - friend constRealGasTransport operator+ - ( - const constRealGasTransport&, - const constRealGasTransport& - ); - - friend constRealGasTransport operator- - ( - const constRealGasTransport&, - const constRealGasTransport& - ); - - friend constRealGasTransport operator* - ( - const scalar, - const constRealGasTransport& - ); - - friend constRealGasTransport operator== - ( - const constRealGasTransport&, - const constRealGasTransport& - ); - - - // Ostream Operator - - friend Ostream& operator<< - ( - Ostream&, - const constRealGasTransport& - ); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#include "constRealGasTransportI.H" - -#ifdef NoRepository -# include "constRealGasTransport.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransportI.H b/src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransportI.H deleted file mode 100644 index 24211f178..000000000 --- a/src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransportI.H +++ /dev/null @@ -1,217 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | - \\ / A nd | For copyright notice see file Copyright - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -inline Foam::constRealGasTransport::constRealGasTransport -( - const Thermo& t, - const scalar mu, - const scalar Pr -) -: - Thermo(t), - mu_(mu), - rPr_(1.0/Pr) -{} - - -template -inline Foam::constRealGasTransport::constRealGasTransport -( - const word& name, - const constRealGasTransport& ct -) -: - Thermo(name, ct), - mu_(ct.mu_), - rPr_(ct.rPr_) -{} - - -// Construct and return a clone -template -inline autoPtr > constTransport::clone -() const -{ - return autoPtr > - ( - new constTransport(*this) - ); -} - - -// Selector from Istream -template -inline autoPtr > constTransport::New -( - Istream& is -) -{ - return autoPtr > - ( - new constRealGasTransport(is) - ); -} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -inline Foam::scalar Foam::constRealGasTransport::mu(const scalar) const -{ - return mu_; -} - - -// CL: for real gas thermo -// Thermal conductivity [W/mK] -template -inline Foam::scalar Foam::constRealGasTransport::kappa(const scalar rho, const scalar T) const -{ - return this->Cp(rho, T)*mu(T)*rPr_; -} - - -// CL: for real gas thermo -// Thermal diffusivity for enthalpy [kg/ms] -template -inline Foam::scalar Foam::constRealGasTransport::alpha(const scalar rho, const scalar T) const -{ - scalar Cp_ = this->Cp(rho, T); - - scalar deltaT = T - specie::Tstd(); - scalar CpBar = - (deltaT*(this->H(rho, T) - this->H(this->rhostd(), specie::Tstd())) + Cp_)/(sqr(deltaT) + 1); - - return Cp_*mu(T)*rPr_/CpBar; -} - - -// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // - -template -inline Foam::constRealGasTransport& Foam::constRealGasTransport::operator= -( - const constRealGasTransport& ct -) -{ - Thermo::operator=(ct); - - mu_ = ct.mu_; - rPr_ = ct.rPr_; - - return *this; -} - - -// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // - -template -inline Foam::constRealGasTransport Foam::operator+ -( - const constRealGasTransport& ct1, - const constRealGasTransport& ct2 -) -{ - Thermo t - ( - static_cast(ct1) + static_cast(ct2) - ); - - scalar molr1 = ct1.nMoles()/t.nMoles(); - scalar molr2 = ct2.nMoles()/t.nMoles(); - - return constRealGasTransport - ( - t, - molr1*ct1.mu_ + molr2*ct2.mu_, - molr1*ct1.rPr_ + molr2*ct2.rPr_ - ); -} - - -template -inline Foam::constRealGasTransport Foam::operator- -( - const constRealGasTransport& ct1, - const constRealGasTransport& ct2 -) -{ - Thermo t - ( - static_cast(ct1) - static_cast(ct2) - ); - - scalar molr1 = ct1.nMoles()/t.nMoles(); - scalar molr2 = ct2.nMoles()/t.nMoles(); - - return constRealGasTransport - ( - t, - molr1*ct1.mu_ - molr2*ct2.mu_, - molr1*ct1.rPr_ - molr2*ct2.rPr_ - ); -} - - -template -inline Foam::constRealGasTransport Foam::operator* -( - const scalar s, - const constRealGasTransport& ct -) -{ - return constRealGasTransport - ( - s*static_cast(ct), - ct.mu_, - ct.rPr_ - ); -} - - -template -inline Foam::constRealGasTransport Foam::operator== -( - const constRealGasTransport& ct1, - const constRealGasTransport& ct2 -) -{ - return ct2 - ct1; -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/thermophysicalProperties b/tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/thermophysicalProperties index cc2797e09..f635849cc 100755 --- a/tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/thermophysicalProperties +++ b/tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/thermophysicalProperties @@ -18,11 +18,11 @@ 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 0.22394 467.6 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; @@ -30,24 +30,19 @@ mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0. //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 839 1.4792e-06 116; - //thermoType realGasHThermo>>>>; - +//mixture CO2 1 44.01 73.773e5 304.13 467.6 0.22394 839 1.4792e-06 116; //thermoType realGasHThermo>>>>; //thermoType realGasHThermo>>>>; - -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; @@ -61,8 +56,8 @@ mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0. // 304.13 --> critical temperatur // 0.22394 --> acentric factor (not needed for redlich kwong) // 467.6 --> critical density (only for aungier redlich kwong) -// 49436.5054 --> 2.849677801e-13 --> 7 heat capacity polynomial coefficent's +// 49436.5054 --> 2.849677801e-13 --> 7 heat capacity polynomial coefficent's --> coefficients, see paper cited in class // .... --> two coefficent's for sutherlandTransport or for the constRealGasTransport model -// 839 --> perfect gas heat capacity Cp0 (J/kgK), needed for constantHeatCapacity +// 839 --> constant cp of perfect gas (used for constantHeatCapacity) // *********************************************************************************************************************** // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/system/controlDict b/tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/system/controlDict index 0c49e76d9..bd5435758 100755 --- a/tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/system/controlDict +++ b/tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/system/controlDict @@ -51,11 +51,4 @@ maxDeltaT 0.01; runTimeModifiable yes; -libs -( - "libIAPWSThermo.so" - "libfreesteam.so" -); - - // ************************************************************************* // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/omega b/tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/omega deleted file mode 100644 index 63e1b3d75..000000000 --- a/tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/omega +++ /dev/null @@ -1,47 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: 2.1.x | -| \\ / A nd | Web: www.OpenFOAM.org | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class volScalarField; - location "0"; - object omega; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -dimensions [0 0 -1 0 0 0 0]; - -internalField uniform 2.6; - -boundaryField -{ - movingWall - { - type compressible::omegaWallFunction; - Cmu 0.09; - kappa 0.41; - E 9.8; - value uniform 2.6; - } - fixedWalls - { - type compressible::omegaWallFunction; - Cmu 0.09; - kappa 0.41; - E 9.8; - value uniform 2.6; - } - frontAndBack - { - type empty; - } -} - - -// ************************************************************************* // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/omega b/tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/omega deleted file mode 100644 index 63e1b3d75..000000000 --- a/tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/omega +++ /dev/null @@ -1,47 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: 2.1.x | -| \\ / A nd | Web: www.OpenFOAM.org | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class volScalarField; - location "0"; - object omega; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -dimensions [0 0 -1 0 0 0 0]; - -internalField uniform 2.6; - -boundaryField -{ - movingWall - { - type compressible::omegaWallFunction; - Cmu 0.09; - kappa 0.41; - E 9.8; - value uniform 2.6; - } - fixedWalls - { - type compressible::omegaWallFunction; - Cmu 0.09; - kappa 0.41; - E 9.8; - value uniform 2.6; - } - frontAndBack - { - type empty; - } -} - - -// ************************************************************************* // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/T b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/T new file mode 100755 index 000000000..e9ec126c0 --- /dev/null +++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/T @@ -0,0 +1,48 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object T; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 333.15; + +boundaryField +{ + outlet + { + type zeroGradient; + } + + inlet1 + { + type fixedValue; + value uniform 273.15; + } + + inlet2 + { + type fixedValue; + value uniform 333.15; + } + + wall + { + type zeroGradient; + } +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/U b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/U new file mode 100644 index 000000000..d5527191c --- /dev/null +++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/U @@ -0,0 +1,47 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 1 0); + +boundaryField +{ + outlet + { + type zeroGradient; + } + + inlet1 + { + type fixedValue; + value uniform (0 1 0); + } + + inlet2 + { + type fixedValue; + value uniform (0 -1 0); + } + + wall + { + type fixedValue; + value uniform (0 0 0); + } +} + +// ************************************************************************* // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/epsilon b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/epsilon new file mode 100644 index 000000000..785a65be2 --- /dev/null +++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/epsilon @@ -0,0 +1,58 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 4.0 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | For copyright notice see file Copyright | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object epsilon; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -3 0 0 0 0]; + +internalField uniform 200; + +boundaryField +{ + outlet + { + type inletOutlet; + inletValue uniform 200; + value uniform 200; + } + inlet1 + { + type compressible::turbulentMixingLengthDissipationRateInlet; + mixingLength 0.005; + phi phi; + k k; + value uniform 200; + } + inlet2 + { + type compressible::turbulentMixingLengthDissipationRateInlet; + mixingLength 0.005; + phi phi; + k k; + value uniform 200; + } + wall + { + type compressible::epsilonWallFunction; + refValue uniform 0; + value uniform 200; + Cmu 0.09; + kappa 0.41; + E 9.8; + } +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/k b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/k new file mode 100644 index 000000000..4e70c50f9 --- /dev/null +++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/k @@ -0,0 +1,52 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 4.0 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | For copyright notice see file Copyright | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object k; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + outlet + { + type zeroGradient; + } + inlet1 + { + type turbulentIntensityKineticEnergyInlet; + intensity 0.05; + U U; + phi phi; + value uniform 1; + } + inlet2 + { + type turbulentIntensityKineticEnergyInlet; + intensity 0.05; + U U; + phi phi; + value uniform 1; + } + wall + { + type compressible::kqRWallFunction; + value uniform 1; + } +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/p b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/p new file mode 100644 index 000000000..6c0f9da9a --- /dev/null +++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/0.org/p @@ -0,0 +1,46 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -2 0 0 0 0]; + +internalField uniform 80e5; + +boundaryField +{ + outlet + { + type fixedMeanValue; + meanValue 8e6; + value uniform 8e+06; + } + + inlet1 + { + type zeroGradient; + } + + inlet2 + { + type zeroGradient; + } + + wall + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties index d284c4cc8..f778483dd 100755 --- a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties +++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties @@ -18,36 +18,31 @@ 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 0.22394 467.6 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>>>>; -//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; - -//thermoType realGasHThermo>>>>; -//mixture CO2 1 44.01 73.773e5 304.13 839 1.4792e-06 116; - //thermoType realGasHThermo>>>>; - +//mixture CO2 1 44.01 73.773e5 304.13 467.6 0.22394 839 1.4792e-06 116; //thermoType realGasHThermo>>>>; //thermoType realGasHThermo>>>>; - -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; -//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; +//thermoType realGasHThermo>>>>; @@ -61,8 +56,8 @@ mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0. // 304.13 --> critical temperatur // 0.22394 --> acentric factor (not needed for redlich kwong) // 467.6 --> critical density (only for aungier redlich kwong) -// 49436.5054 --> 2.849677801e-13 --> 7 heat capacity polynomial coefficent's -// .... --> two coefficent's for sutherlandRealGasTransport or for the constRealGasTransport model -// 839 --> perfect gas heat capacity Cp0 (J/kgK), needed for constantHeatCapacity +// 49436.5054 --> 2.849677801e-13 --> 7 heat capacity polynomial coefficent's --> coefficients, see paper cited in class +// .... --> two coefficent's for sutherlandTransport or for the constRealGasTransport model +// 839 --> constant cp of perfect gas (used for constantHeatCapacity) // *********************************************************************************************************************** //