Advanced Fluid Properties

Use of Extended Corresponding States for Transport Properties

For a limited selection of fluids, correlations are provided for the viscosity and the thermal conductivity. But for many fluids, no correlations are available, and therefore other methods must be employed. The extended corresponding states is a method of estimating the transport properties of a fluid by analogy with the transport properties of a fluid that are well defined.

Implementing the ECS method is quite a challenge, but CoolProp is one of the only fluid property databases that properly implements it. And the onlyopen-source package that does. A multi-step method is required, which is hopefully clearly laid out here.

To begin with, the reference fluid must be selected that the fluid of interest will be compared with. Ideally the shape of the molecules should be similar, but in practice, most fluids use R134a as the reference fluid since its thermodynamic and transport properties are well quantified with reference-quality correlations.

Once the reference fluid has been selected, the conformal state of the reference fluid must be determined. The conformal state is the state at which the transport properties of the reference fluid and the fluid of interest are (in theory) the same. In practice, at low densities the shape factors are assumed to be unity, and the conformal temperature and molar density are obtained from

T0=TTc0Tcj
¯ρ0=¯ρ¯ρ0c¯ρjc

Exact solution for the conformal temperature

If you have Helmholtz EOS for both the fluid and the reference fluid, you need to find a conformal temperature for the reference fluid that will yield the same compressibility factor and the residual Helmholtz energy

Zj(Tj,ρj)=Z0(T0,ρ0)
αrj(Tj,ρj)=αr0(T0,ρ0)

where “j” is for the fluid of interest, and the subscript “0” is for the reference fluid. The left side of each equation is already known from the equation of state. Literature suggests that solving for T0 and ρ0 directly is quite challenging. See McLinden 2000 or Klein 1997.

Alternatively, if the shape factors θ and ϕ are known, either from correlation or otherwise, the conformal temperature and density can be calculated directly.

T0=Tf=TTc0Tcjθ(Tj,ρj)
ρ0=ρh=ρρc0ρcjϕ(Tj,ρj)

Conversion from ideal gas term to Helmholtz energy term

Much of the time the coefficients for the ideal-gas part of the Helmholtz energy are given directly, but sometimes only the gas-specific heat is provided. Therefore you need to be able to go from specific heat to ideal-gas Helmholtz Energy. The ideal-gas Helmholtz energy is given by Equation 23 from Lemmon, 2004, Equations of State for Mixtures of R-32, R-125, R-134a, R-143a, and R-152a, J. Phys. Chem. Ref. Data, Vol. 33, No. 2, 2004 or

a0=RT+RTlnρTρ0T0+h00Ts00+TT0c0p(T)dTTTT0c0p(T)TdT

non-dimensionalizing

α0=a0RT=1+lnρTρ0T0+h00RTs00R+1RTTT0c0p(T)dT1RTT0c0p(T)TdT

Now we might want to do a change of variable in the integrals. If so, do a u-substitution in the integrals.

T=Tcτ

where

dT=Tcτ2dτ
α0=1+lnρTρ0T0+h00RTs00R+1RTττ0c0p(T)(Tcτ2dτ)1Rττ0c0p(τ)T(Tcτ2dτ)

Simplifying and factoring the τ term yields

α0=1+lnρTρ0T0+h00RTs00RτRττ0c0p(τ)τ2dτ+1Rττ0c0p(τ)τdτ

which finally yields the solution as of Equation 3 from Lemmon, 2003 (and others)

The specific-heat contribution can then be taken as a sum of the contributions

for a term of the form

c0pR=(B/T)2exp(B/T)(exp(B/T)1)2

the contribution is found from

1TTT0(B/T)2exp(B/T)(exp(B/T)1)2dTTT0(B/T)2exp(B/T)(exp(B/T)1)21TdT
1T[Bexp(B/T)1|TT0[BT(1exp(B/T)1+1)log[exp(B/T)1]|TT0dT

Factor out a B, First two terms cancel, leaving

[BTlog[exp(B/T)1]|TT0dT
[log[exp(B/T)1]BT|TT0dT
log[exp(B/T)1]BT(log[exp(B/T0)1]BT0)

or in terms of τ

log[exp(Bτ/Tc)1]BτTc(log[exp(Bτ0/Tc)1]Bτ0Tc)

for a term of the form

c0pR=c

the contribution is found from

1TTT0cdTTT0cTdT
cT(TT0)clog(T/T0)

or in terms of τ

ccT0τTc+clog(τ/τ0)

for a term of the form

c0pR=cTt,t0

the contribution is found from

1TTT0cTtdTTT0cTtTdT
cT(Tt+1t+1Tt+10t+1)c(TttTt0t)
cTt(1t+11t)cTt+10T(t+1)+cTt0t

or in terms of τ

cTtcτt(1t+11t)cTt+10τTc(t+1)+cTt0t

These terms can be summarized by the following table:

c0pRTermα0TermClassNameak(bk/T)2exp(bk/T)(exp(bk/T)1)2akln[1exp(bkτTc)]phi0_Planck_Einstein(a,b/Tc,[iStart,iEnd])ac(b/T)2exp(b/T)(cexp(b/T)+1)2aln[c+exp(bτTc)]phi0_Planck_Einstein2(a,b/Tc,c)yuckakτbkphi0_power(a,b,[iStart,iEnd])aaaττ0+aln(ττ0)phi0_cp0_constant(a,Tc,T0)a1+a2(a3/Tsinh(a3/T))2+a4(a5/Tcosh(a5/T))2yuckphi0_cp0_AlyLee(a,Tc,T0,R)n/alog(δ)+a1+a2τphi0_lead(a1,a2)n/aalogτphi0_logtau(a)

If the reference enthalpy is known, you can determine the constants from

h0RT=τ[(α0τ)δ+(αrτ)δ]+δ(αrδ)τ+1
(α0τ)δ=1τ(h0RTδ(αrδ)τ1)(αrτ)δ

For the specific heat The two integral terms are

τRττ0c0pτ2dτ+1Rττ0c0pτdτ

First derivative

ddτ[τRττ0c0pτ2dτ+1Rττ0c0pτdτ]=c0pτR1Rττ0c0pτ2dτ+c0pτR=1Rττ0c0pτ2dτ

Second Derivative

d2dτ2[τRττ0c0pτ2dτ+1Rττ0c0pτdτ]=ddτ[1Rττ0c0pτ2dτ]=c0pτ2R

Converting Bender and mBWR EOS

If the EOS is of the form

pρRT=Z(T,ρ)=1+iniTsiρri+iniTsiρriexp(γi(ρρc)2)

To convert to standard power form in CoolProp, use

δidiaiτtiδdi1=iniTsiρri=ini(Tcτ)si(ρcδ)ri=iniTsicρricτsiδri

The left-hand-side is the derivative of the residual Helmholtz energy with respect to delta times the reduced density since

pρRT=1+δ(αrδ)τ

where

δ:di1+1=ridi=ri
τ:ti=si
c:diai=niTsicρric
p=ρRT+iniTsiρri+iniTsiρriexp(γi(ρρc)2)(Eq3.28)
pρRT=1+iniRTsi1ρri1+iniRTsi1ρri1exp(γi(ρρc)2)
δidiaiτtiδdi1=iniR(Tcτ)si1(ρcδ)ri1=iniRTsi1cρri1cτ(si1)δri1
δ:1+di1=ri1
τ:ti=(si1)
c:diai=niRTsi1cρri1c

In the Bender EOS, for the exponential part you have terms that can be converted to reduced form

aiδdiτtiexp(γδ2)

which yields the terms in the following table (from Span, 2000)

\multicolumn4cFrom Bender\multicolumn4cPower termiditiγiniditiγi1423γn14/(2γ)+n17/(2γ2)0301524γn15/(2γ)+n17/(2γ2)0401625γn16/(2γ)+n17/(2γ2)0501743γn14/(2γ)n17/(2γ2)03γ1844γn15/(2γ)n18/(2γ2)04γ1945γn16/(2γ)n19/(2γ2)05γ20n17/(2γ)23γ21n18/(2γ)24γ22n19/(2γ)25γ

Warning

If the terms in the EOS are in terms of T and ρ rather than τ and δ, make sure to multiply appropriately by the critical densities in the exponential term. For instance in Polt paper, the first constant should be n14ρ2c/(2γ)+n17ρ4c/(2γ2)/T3c Be careful!