Merge Request #36: corrections in the real gas library. Author: Christian Lucas. Merge: Henrik Rusche

Conflicts:
	src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.C
	src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacityI.H
	src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.C
	src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.H
	src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransportI.H
	tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/thermophysicalProperties
	tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/omega
	tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/omega
	tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties
This commit is contained in:
Henrik Rusche 2016-06-23 16:41:24 +02:00
commit ee238cd77e
27 changed files with 445 additions and 701 deletions

View file

@ -143,6 +143,7 @@ void Foam::calculateProperties_h
region=freesteam_region(S);
//CL:Liquid phase
if (region==1)
{
p=S.R1.p;
@ -168,6 +169,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;
@ -193,6 +195,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;
@ -234,9 +237,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;
@ -271,7 +275,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);
@ -282,7 +286,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))
@ -385,6 +389,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)
@ -398,6 +403,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)
@ -411,6 +417,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)
{
@ -433,9 +440,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;
@ -466,7 +474,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);
@ -476,7 +484,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))

View file

@ -126,7 +126,7 @@ void Foam::aungierRedlichKwong::write(Ostream& os) const
Foam::Ostream& Foam::operator<<(Ostream& os, const aungierRedlichKwong& ark)
{
os << static_cast<const specie&>(ark)<< token::SPACE
<< ark.pcrit_ << tab<< ark.Tcrit_<< tab<<ark.azentricFactor_<< tab<<ark.rhocrit_;
<< ark.pcrit_ << tab<< ark.Tcrit_<< tab<<ark.rhocrit_ << tab<<ark.azentricFactor_;
os.check("Ostream& operator<<(Ostream& os, const aungierRedlichKwong& st)");
return os;

View file

@ -96,7 +96,7 @@ class aungierRedlichKwong
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
@ -156,10 +156,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

View file

@ -136,12 +136,14 @@ inline scalar aungierRedlichKwong::azentricFactor() const
inline void aungierRedlichKwong::updateModelCoefficients(const scalar T)const
{
aSave=a0_*pow(T/Tcrit_, -n_);
daSave=-a0_*n_*pow(T/Tcrit_, -n_)/T;
d2aSave=a0_*(n_*n_ + n_)/(T*T)*pow(T/Tcrit_, -n_);
scalar TTc = pow(T/Tcrit_, -n_);
aSave = a0_*TTc;
daSave = -a0_*n_*TTc/T;
d2aSave = a0_*(n_*n_ + n_)/(T*T)*TTc;
//CL: saving the temperature at which the coefficients are valid
TSave=T;
TSave = T;
}
@ -224,6 +226,7 @@ inline scalar aungierRedlichKwong::n() const
inline scalar aungierRedlichKwong::p(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR()*T/(Vm - b_ + c_) - a(T)/(Vm*(Vm + b_));
}
@ -249,6 +252,7 @@ inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
inline scalar aungierRedlichKwong::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR()/(Vm - b_ + c_)-dadT(T)/(Vm*(Vm + b_));
}
@ -279,6 +283,7 @@ inline scalar aungierRedlichKwong::integral_p_dv
) const
{
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_;
}
@ -293,6 +298,7 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv
) const
{
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_;
}
@ -302,6 +308,7 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv
inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return -d2adT2(T)/(Vm*(Vm + b_));
}
@ -312,19 +319,27 @@ inline scalar aungierRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar cpVm = c_ + Vm;
scalar bpVm3 = pow(b_ + Vm, 3);
return -2*
(
a(T)*
(
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)
b5_ - 3*b4_*c_ + 3*b3_*pow(c_ - Vm, 2)
+ cpVm*
(
- b2_*(c2_ - 7*c_*Vm + Vm2)
+ cpVm*
(
3*b_*Vm*(2*Vm - c_)
- 3*Vm2*cpVm
)
)
)
+ this->RR()*T*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
+ this->RR()*T*Vm3*bpVm3
)
/(Vm3*pow(b_ - c_ - Vm, 3)*pow(b_ + Vm, 3));
/(Vm3*pow(b_ - c_ - Vm, 3)*bpVm3);
}
@ -336,15 +351,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;
}
@ -359,13 +376,17 @@ inline scalar aungierRedlichKwong::d2pdvdT
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar cpVm = c_ + Vm;
scalar cpVm2 = cpVm*cpVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return
(
dadT(T)*(b3_ - 2*b2_*c_ + b_*(c_ + Vm)*(c_ - 3*Vm)+2*Vm*pow(c_ + Vm, 2))
- this->RR()*Vm2*(b2_ + 2*b_*Vm + Vm2)
dadT(T)*(b3_ - 2*b2_*c_ + b_*cpVm*(c_ - 3*Vm) + 2*Vm*cpVm2)
- this->RR()*Vm2*bpVm2
)
/(Vm2*pow(b_ - c_ - Vm, 2)*pow(b_ + Vm, 2));
/(Vm2*pow(b_ - c_ - Vm, 2)*bpVm2);
}
@ -378,6 +399,7 @@ inline scalar aungierRedlichKwong::integral_d2pdT2_dv
) const
{
scalar Vm = this->W()/rho;
return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_;
}

View file

@ -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

View file

@ -116,9 +116,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;
@ -223,7 +225,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));
/pow(b_ - Vm, 6);
}
@ -233,6 +235,7 @@ inline scalar pengRobinson::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
return this->RR()/(Vm - b_) - dadT(T)/(Vm2 + 2*b_*Vm - b2_);
}
@ -283,6 +286,7 @@ inline scalar pengRobinson::integral_dpdT_dv
{
scalar Vm = this->W()/rho;
scalar root2=pow(2, 0.5);
return
- root2*dadT(T)*log(b_*(1 - root2) + Vm)/(4*b_)
+ this->RR()*log(Vm - b_) + root2*dadT(T)*log(b_*(root2 + 1) + Vm)/(4*b_);
@ -294,7 +298,8 @@ inline scalar pengRobinson::d2pdT2(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
return -d2adT2(T)/(Vm2 + 2*b_*Vm-b2_);
return -d2adT2(T)/(Vm2 + 2*b_*Vm - b2_);
}
@ -319,7 +324,7 @@ 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));
/pow(b_ - Vm, 9);
}
@ -331,15 +336,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;
}
@ -362,7 +369,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));
/pow(b_ - Vm, 6);
}

View file

@ -121,6 +121,7 @@ inline scalar redlichKwong::b() const
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_));
}
@ -132,12 +133,14 @@ 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);
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*bpVm2
)
/(sqrt(T)*Vm2*pow((b_ + Vm), 2)*pow( (b_ - Vm), 2));
/(T05*Vm2*bpVm2*pow(b_ - Vm, 2));
}
@ -146,7 +149,8 @@ inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
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);
return 0.5*a_/(pow(T, 1.5)*Vm*(b_ + Vm)) - this->RR()/(b_ - Vm);
}
@ -176,6 +180,7 @@ inline scalar redlichKwong::integral_p_dv
) const
{
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));
}
@ -190,6 +195,7 @@ inline scalar redlichKwong::integral_dpdT_dv
) const
{
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_);
}
@ -211,15 +217,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*pow((b_ + Vm), 3)*pow(Vm - b_, 3))
);
}
@ -228,15 +235,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;
}

View file

@ -83,11 +83,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),

View file

@ -92,7 +92,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
@ -144,10 +144,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

View file

@ -114,12 +114,14 @@ 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;
TSave = T;
}
@ -221,6 +223,7 @@ inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
inline scalar soaveRedlichKwong::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR()/(Vm - b_) - dadT(T)/(Vm*(Vm + b_));
}
@ -251,6 +254,7 @@ inline scalar soaveRedlichKwong::integral_p_dv
) const
{
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_;
}
@ -266,6 +270,7 @@ inline scalar soaveRedlichKwong::integral_dpdT_dv
) const
{
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_;
}
@ -275,6 +280,7 @@ inline scalar soaveRedlichKwong::integral_dpdT_dv
inline scalar soaveRedlichKwong::d2pdT2(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return -d2adT2(T)/(Vm*(Vm + b_));
}
@ -305,15 +311,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;
}
@ -347,6 +355,7 @@ inline scalar soaveRedlichKwong::integral_d2pdT2_dv
) const
{
scalar Vm = this->W()/rho;
return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_;
}

View file

@ -38,14 +38,15 @@ template<class equationOfState>
Foam::constantHeatCapacity<equationOfState>::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)");

View file

@ -108,7 +108,7 @@ class constantHeatCapacity
{
// Private data
//CL: spec. cp
//CL: specific values
scalar Cp0_;
//CL: molar values
scalar cp0_;

View file

@ -169,7 +169,6 @@ inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::cv0
const scalar T
) const
{
return this->cp0(T)-this->RR();
}
@ -183,10 +182,14 @@ inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::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
// 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
(
@ -278,7 +281,7 @@ inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::s
- integral_dpdT_dv_std
+ this->s0(T)
- s0_std
- this->RR*log(T/this->Tstd())
- this->RR()*log(T/this->Tstd())
+ this->integral_dpdT_dv(rho, T);
}

View file

@ -45,12 +45,13 @@ Foam::nasaHeatCapacityPolynomial<equationOfState>::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)");

View file

@ -157,15 +157,23 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::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_)
)
)
)
)
);
}
@ -188,15 +196,23 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s0
const scalar T
) const
{
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
);
scalar oneOverT = 1/T;
return this->RR()*
(
oneOverT*(- 0.5*this->a1_*oneOverT - this->a2_)
+ this->a3_*log(T)
+ T*
(
this->a4_ + T*
(
0.5*this->a5_ + T*
(
this->a6_/3 + T*(0.25*this->a7_)
)
)
)
);
}
@ -207,15 +223,13 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cp0
const scalar T
) const
{
scalar oneOverT = 1/T;
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)
oneOverT*(this->a1_*oneOverT + this->a2_)
+ this->a3_
+ T*(this->a4_ + T*(this->a5_ + T*(this->a6_ + T*this->a7_)))
);
}
@ -241,19 +255,21 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::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
return
min
// 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)
/min(this->dpdv(rho, T),-1)
this->cv(rho, T)
- T*pow(this->dpdT(rho, T), 2)/min(this->dpdv(rho, T), -1)
)
);
}
@ -267,7 +283,10 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::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));
return fabs
(
this->cv(rho,T)-T*pow((this->dpdT(rho, T)),2)/min(this->dpdv(rho, T), -1)
);
}
@ -292,7 +311,9 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::h
const scalar T
) const
{
return this->e(rho,T)+this->p(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();
}
@ -308,11 +329,12 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::e
{
return
(
-this->Tstd()*integral_dpdT_dv_std
+integral_p_dv_std
+this->e0(T)-e0_std
+T*this->integral_dpdT_dv(rho,T)
-this->integral_p_dv(rho,T)
- this->Tstd()*integral_dpdT_dv_std
+ integral_p_dv_std
+ this->e0(T)
- e0_std
+ T*this->integral_dpdT_dv(rho, T)
- this->integral_p_dv(rho, T)
);
}
@ -328,10 +350,14 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::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)
);
}

View file

@ -131,7 +131,7 @@ inline Foam::realGasSpecieThermo<thermo>::realGasSpecieThermo
template<class thermo>
inline Foam::scalar Foam::realGasSpecieThermo<thermo>::gamma(const scalar rho, const scalar T ) const
{
return -1/(rho*this->p(rho,T))*this->cp(rho,T)/this->cv(rho,T)*this->dpdv(rho,T);
return -1/(rho*this->p(rho, T))*this->cp(rho, T)/this->cv(rho, T)*this->dpdv(rho, T);
}
@ -166,7 +166,7 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::Cv( const scalar rho, con
template<class thermo>
inline Foam::scalar Foam::realGasSpecieThermo<thermo>::H(const scalar rho, const scalar T) const
{
return this->h(rho, T)/this->W();
return this->h(rho, T)/this->W();
}
@ -187,7 +187,7 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::E(const scalar rho, const
template<class thermo>
inline Foam::scalar Foam::realGasSpecieThermo<thermo>::G(const scalar rho, const scalar T) const
{
return this->g(rho, T)/this->W();
return this->g(rho, T)/this->W();
}
@ -207,12 +207,12 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiH
) const
{
scalar beta=this->isobarExpCoef(rho,T);
scalar beta=this->isobarExpCoef(rho, T);
return
-(
(T*beta*beta-beta)/this->Cp(rho,T)
-this->isothermalCompressiblity(rho,T)*rho
(T*beta*beta-beta)/this->Cp(rho, T)
- this->isothermalCompressiblity(rho, T)*rho
);
}
@ -227,20 +227,18 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiE
) const
{
scalar V = 1/rho;
scalar cp=this->Cp(rho,T);
scalar beta=this->isobarExpCoef(rho,T);
scalar cp = this->Cp(rho, T);
scalar beta = this->isobarExpCoef(rho, T);
return
-(
(
T*pow(beta,2)*V
-this->isothermalCompressiblity(rho,T)*cp
)
/
(
cp*V
-beta*this->p(rho,T)*pow(V,2)
)
T*beta*beta*V - this->isothermalCompressiblity(rho, T)*cp
)
/
(
cp*V - beta*this->p(rho, T)*V*V
)
);
}

View file

@ -1,56 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
-------------------------------------------------------------------------------
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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "constRealGasTransport.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Thermo>
Foam::constRealGasTransport<Thermo>::constRealGasTransport(Istream& is)
:
Thermo(is),
mu_(readScalar(is)),
rPr_(1.0/readScalar(is))
{
is.check("constRealGasTransport::constRealGasTransport(Istream& is)");
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
template<class Thermo>
Foam::Ostream& Foam::operator<<(Ostream& os, const constRealGasTransport<Thermo>& ct)
{
operator<<(os, static_cast<const Thermo&>(ct));
os << tab << ct.mu_ << tab << 1.0/ct.rPr_;
os.check("Ostream& operator<<(Ostream&, const constRealGasTransport&)");
return os;
}
// ************************************************************************* //

View file

@ -1,202 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
-------------------------------------------------------------------------------
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 <http://www.gnu.org/licenses/>.
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 Thermo> class constRealGasTransport;
template<class Thermo>
inline constRealGasTransport<Thermo> operator+
(
const constRealGasTransport<Thermo>&,
const constRealGasTransport<Thermo>&
);
template<class Thermo>
inline constRealGasTransport<Thermo> operator-
(
const constRealGasTransport<Thermo>&,
const constRealGasTransport<Thermo>&
);
template<class Thermo>
inline constRealGasTransport<Thermo> operator*
(
const scalar,
const constRealGasTransport<Thermo>&
);
template<class Thermo>
inline constRealGasTransport<Thermo> operator==
(
const constRealGasTransport<Thermo>&,
const constRealGasTransport<Thermo>&
);
template<class Thermo>
Ostream& operator<<
(
Ostream&,
const constRealGasTransport<Thermo>&
);
/*---------------------------------------------------------------------------*\
Class constRealGasTransport Declaration
\*---------------------------------------------------------------------------*/
template<class Thermo>
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+ <Thermo>
(
const constRealGasTransport&,
const constRealGasTransport&
);
friend constRealGasTransport operator- <Thermo>
(
const constRealGasTransport&,
const constRealGasTransport&
);
friend constRealGasTransport operator* <Thermo>
(
const scalar,
const constRealGasTransport&
);
friend constRealGasTransport operator== <Thermo>
(
const constRealGasTransport&,
const constRealGasTransport&
);
// Ostream Operator
friend Ostream& operator<< <Thermo>
(
Ostream&,
const constRealGasTransport&
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "constRealGasTransportI.H"
#ifdef NoRepository
# include "constRealGasTransport.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -1,217 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
-------------------------------------------------------------------------------
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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Thermo>
inline Foam::constRealGasTransport<Thermo>::constRealGasTransport
(
const Thermo& t,
const scalar mu,
const scalar Pr
)
:
Thermo(t),
mu_(mu),
rPr_(1.0/Pr)
{}
template<class Thermo>
inline Foam::constRealGasTransport<Thermo>::constRealGasTransport
(
const word& name,
const constRealGasTransport& ct
)
:
Thermo(name, ct),
mu_(ct.mu_),
rPr_(ct.rPr_)
{}
// Construct and return a clone
template<class thermo>
inline autoPtr<constRealGasTransport<Thermo> > constTransport<thermo>::clone
() const
{
return autoPtr<constRealGasTransport<Thermo> >
(
new constTransport<thermo>(*this)
);
}
// Selector from Istream
template<class thermo>
inline autoPtr<constRealGasTransport<Thermo> > constTransport<thermo>::New
(
Istream& is
)
{
return autoPtr<constRealGasTransport<Thermo> >
(
new constRealGasTransport<thermo>(is)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Thermo>
inline Foam::scalar Foam::constRealGasTransport<Thermo>::mu(const scalar) const
{
return mu_;
}
// CL: for real gas thermo
// Thermal conductivity [W/mK]
template<class thermo>
inline Foam::scalar Foam::constRealGasTransport<thermo>::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<class thermo>
inline Foam::scalar Foam::constRealGasTransport<thermo>::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<class Thermo>
inline Foam::constRealGasTransport<Thermo>& Foam::constRealGasTransport<Thermo>::operator=
(
const constRealGasTransport<Thermo>& ct
)
{
Thermo::operator=(ct);
mu_ = ct.mu_;
rPr_ = ct.rPr_;
return *this;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template<class Thermo>
inline Foam::constRealGasTransport<Thermo> Foam::operator+
(
const constRealGasTransport<Thermo>& ct1,
const constRealGasTransport<Thermo>& ct2
)
{
Thermo t
(
static_cast<const Thermo&>(ct1) + static_cast<const Thermo&>(ct2)
);
scalar molr1 = ct1.nMoles()/t.nMoles();
scalar molr2 = ct2.nMoles()/t.nMoles();
return constRealGasTransport<Thermo>
(
t,
molr1*ct1.mu_ + molr2*ct2.mu_,
molr1*ct1.rPr_ + molr2*ct2.rPr_
);
}
template<class Thermo>
inline Foam::constRealGasTransport<Thermo> Foam::operator-
(
const constRealGasTransport<Thermo>& ct1,
const constRealGasTransport<Thermo>& ct2
)
{
Thermo t
(
static_cast<const Thermo&>(ct1) - static_cast<const Thermo&>(ct2)
);
scalar molr1 = ct1.nMoles()/t.nMoles();
scalar molr2 = ct2.nMoles()/t.nMoles();
return constRealGasTransport<Thermo>
(
t,
molr1*ct1.mu_ - molr2*ct2.mu_,
molr1*ct1.rPr_ - molr2*ct2.rPr_
);
}
template<class Thermo>
inline Foam::constRealGasTransport<Thermo> Foam::operator*
(
const scalar s,
const constRealGasTransport<Thermo>& ct
)
{
return constRealGasTransport<Thermo>
(
s*static_cast<const Thermo&>(ct),
ct.mu_,
ct.rPr_
);
}
template<class Thermo>
inline Foam::constRealGasTransport<Thermo> Foam::operator==
(
const constRealGasTransport<Thermo>& ct1,
const constRealGasTransport<Thermo>& ct2
)
{
return ct2 - ct1;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -18,11 +18,11 @@ FoamFile
//CL: List of possible real gas models
thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//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<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
mixture CO2 1 44.01 73.773e5 304.13 467.6 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
@ -30,24 +30,19 @@ mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<redlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 839 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<aungierRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 467.6 0.22394 839 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<constantHeatCapacity<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<constantHeatCapacity<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<constantHeatCapacity<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<constantHeatCapacity<soaveRedlichKwong>>>>>;
@ -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 sutherlandTransport 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 coefficents --> coefficients, see paper cited in class
// .... --> two coefficents for sutherlandTransport or for the constRealGasTransport model
// 839 --> constant cp of perfect gas (used for constantHeatCapacity)
// *********************************************************************************************************************** //

View file

@ -51,11 +51,4 @@ maxDeltaT 0.01;
runTimeModifiable yes;
libs
(
"libIAPWSThermo.so"
"libfreesteam.so"
);
// ************************************************************************* //

View file

@ -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;
}
}
// ************************************************************************* //

View file

@ -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);
}
}
// ************************************************************************* //

View file

@ -3,7 +3,7 @@
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{
@ -11,35 +11,46 @@ FoamFile
format ascii;
class volScalarField;
location "0";
object omega;
object epsilon;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 -1 0 0 0 0];
dimensions [0 2 -3 0 0 0 0];
internalField uniform 2.6;
internalField uniform 200;
boundaryField
{
movingWall
outlet
{
type compressible::omegaWallFunction;
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;
value uniform 2.6;
}
fixedWalls
{
type compressible::omegaWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 2.6;
}
frontAndBack
{
type empty;
}
}

View file

@ -3,7 +3,7 @@
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{
@ -11,35 +11,40 @@ FoamFile
format ascii;
class volScalarField;
location "0";
object omega;
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 -1 0 0 0 0];
dimensions [0 2 -2 0 0 0 0];
internalField uniform 2.6;
internalField uniform 1;
boundaryField
{
movingWall
outlet
{
type compressible::omegaWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 2.6;
type zeroGradient;
}
fixedWalls
inlet1
{
type compressible::omegaWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 2.6;
type turbulentIntensityKineticEnergyInlet;
intensity 0.05;
U U;
phi phi;
value uniform 1;
}
frontAndBack
inlet2
{
type empty;
type turbulentIntensityKineticEnergyInlet;
intensity 0.05;
U U;
phi phi;
value uniform 1;
}
wall
{
type compressible::kqRWallFunction;
value uniform 1;
}
}

View file

@ -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;
}
}
// ************************************************************************* //

View file

@ -18,36 +18,31 @@ FoamFile
//CL: List of possible real gas models
thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//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<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<redlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 839 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<aungierRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 467.6 0.22394 839 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<constantHeatCapacity<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<constantHeatCapacity<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<constantHeatCapacity<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<constantHeatCapacity<soaveRedlichKwong>>>>>;
@ -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 coefficents --> coefficients, see paper cited in class
// .... --> two coefficents for sutherlandRealGasTransport or for the constRealGasTransport model
// 839 --> constant cp of perfect gas (used for constantHeatCapacity)
// *********************************************************************************************************************** //