Enhancement: Joule Thomson coefficient Addition

[defencedog]

We know for ideal gases JT coeff is zero; but can we add this parameter or somehow calculate it?

Important for HVAC community.

It is possible to express the Joule-Thomson coefficient in terms of measurable properties

https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Physical_Chemistry_(Fleming)/04:_Putting_the_First_Law_to_Work/4.05:_The_Joule-Thomson_Effect

[chmarti1]

I'm glad you raise this topic - it's been on my list of things to add in future releases. In the meantime, the easiest way to calculate the Joule-Thomson coefficient is numerically.

Most texts recommend calculating the JT coefficient from the partial derivative of volume with respect to temperature at constant pressure (thermal expansion coefficient). The easiest way to calculate that without coding your own methods is to perturb specific volume at constant pressure.

>>> dT=.01; T=300.; dp=1e-5; p=1.01325
>>> alpha = (n2.v(T=T+dT,p=p) - n2.v(T=T,p=p)) / dT
>>> alpha
array([0.00293603])

I don't recommend using this approach in PYroMat, though. It's much easier to calculate the derivatives of enthalpy with respect to temperature and pressure directly. This example calculates the JT coefficient in Kelvin per bar.

>>> dhdp = (n2.h(p=p+dp,T=T) - n2.h(p=p,T=T)) / dp
>>> dhdt = (n2.h(p=p,T=T+dt) - n2.h(p=p,T=T)) / dt      # Or just use dhdt = n2.cp(p=p,T=T)
>>> mu = -dhdp / dhdt
>>> mu
array([0.21195838])

These methods are a little inelegant, but they work well. The back-end doesn't use these kinds of methods - it calculates these derivatives analytically from the EOS coefficients, so it's more efficient.

[chmarti1]

There's another comment worth making here. The Joule-Thomson coefficient is classically used to approximate the temperature after the isenthalpic throttle in a refrigeration cycle. PYroMat can do that job without the intermediate step of calculating "mu."

In this example, we get the temperature after a throttling operation from 300K, 20bar down to 1bar.

>>> r = pm.get('mp.C2H2F4')
>>> h1 = r.h(T=300., p=20.)
>>> T2 = r.T(h=h1, p=1.)
>>> T2
array([246.78990841])

[defencedog]

Extremely Thank You

[defencedog]

@chmarti1 Kindly guide on what basis are these values dT=.01, dp=1e-5 ...small but how small?

[chmarti1]

Good question. As a percentage of the value being perturbed, the "delta" should be larger than the machine error in the calculation (1e-6 is the WORST error PYroMat usually tolerates in numerical inversion problems, 1e-15 is more typical), but small enough that changes in the perturbation does not significantly change the value of the slope.

For example, here are three ways of calculating the "same" number:

>>> n2 = pm.get('mp.N2')
>>> n2.cp(T=T,p=p)    # First, what is the "True" value?
array([1.04135039])
>>> T = 300; dT = .01; p=1.01325; dp = 1e-5
>>> dhdt = (n2.h(T=T+dT,p=p) - n2.h(T=T,p=p)) / dT    # Now, get the same value numerically
>>> dhdt
array([1.04135041])    # Pretty good!
>>> dT = .02
>>> dhdt = (n2.h(T=T+dT,p=p) - n2.h(T=T,p=p)) / dT    # Now, repeat with a larger dT
>>> dhdt
array([1.04135043])    # Still pretty good!

[defencedog]

@chmarti1 perturbation does not significantly change the value of the slope. Done & Done, thanks