Enhancement: Joule Thomson coefficient Addition about pyromat HOT 6 OPEN

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])

from pyromat.

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!

from pyromat.

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.

from pyromat.

Extremely Thank You

from pyromat.

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

from pyromat.

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

from pyromat.