Nomenclature:

$d$ means a derivative ALONG the saturation line,
$\partial$ means a partial derivative AT the saturation line (or anywhere in the single phase region).

References:

Krafcik and Velasco, DOI 10.1119/1.4858403
Thorade and Saadat, DOI 10.1007/s12665-013-2394-z
https://en.wikipedia.org/wiki/Product_rule
https://en.wikipedia.org/wiki/Quotient_rule
https://en.wikipedia.org/wiki/Triple_product_rule

Clausius-Clapeyron

Clausius-Clapeyron p/T

\begin{equation} \frac{dp}{dT} = \frac{s''-s'}{v''-v'} \end{equation}

derivative of Clausius-Clapeyron \begin{equation} \begin{split} \frac{d^2 p}{d T^2} &= \frac{\left(\frac{ds''}{dT}-\frac{ds'}{dT}\right) (v''-v')}{(v''-v')^2}

\frac{(s''-s') \left(\frac{dv''}{dT}-\frac{dv'}{dT}\right)}{(v''-v')^2} \ &= \frac{\left(\frac{ds''}{dT}-\frac{ds'}{dT}\right)}{(v''-v')}
\frac{dp}{dT} \frac{\left(\frac{dv''}{dT}-\frac{dv'}{dT}\right)}{(v''-v')} \end{split} \end{equation}

Clausius-Clapeyron T/p \begin{equation} \frac{dT}{dp} = \frac{v''-v'}{s''-s'} \end{equation}

derivative of Clausius-Clapeyron \begin{equation} \begin{split} \frac{d^2 T}{d p^2} &= \frac{\left(\frac{dv''}{dp}-\frac{dv'}{dp}\right) (s''-s')}{(s''-s')^2}

\frac{(v''-v') \left(\frac{ds''}{dp}-\frac{ds'}{dp}\right)}{(s''-s')^2} \ &= \frac{\left(\frac{dv''}{dp}-\frac{dv'}{dp}\right)}{(s''-s')}
\frac{dT}{dp} \frac{\left(\frac{ds''}{dp}-\frac{ds'}{dp}\right)}{(s''-s')} \end{split} \end{equation}

where \begin{equation} \begin{split} \frac{ds}{dT} &= \left(\frac{\partial s}{\partial T}\right)_p + \left(\frac{\partial s}{\partial p}\right)_T \frac{dp}{dT}\\ \frac{dv}{dT} &= \left(\frac{\partial v}{\partial T}\right)_p + \left(\frac{\partial v}{\partial p}\right)_T \frac{dp}{dT}\\ \frac{dv}{dp} &= \left(\frac{\partial v}{\partial p}\right)_T + \left(\frac{\partial v}{\partial T}\right)_p \frac{dT}{dp}\\ \frac{ds}{dp} &= \left(\frac{\partial s}{\partial p}\right)_T + \left(\frac{\partial s}{\partial T}\right)_p \frac{dT}{dp} \end{split} \end{equation}

Temporary Names

Introduce temporary names for some of the partial derivatives wrt $p$ and $T$: \begin{equation} \begin{split} A &= \left(\frac{\partial \rho}{\partial T}\right)_p \\ B &= \left(\frac{\partial \rho}{\partial p}\right)_T \\ C &= \left(\frac{\partial s}{\partial T}\right)_p \\ E &= \left(\frac{\partial s}{\partial p}\right)_T \\ G &= \left(\frac{\partial h}{\partial T}\right)_p \\ K &= \left(\frac{\partial h}{\partial p}\right)_T \end{split} \end{equation}

Introduce temporary names for some of the partial derivatives wrt $\rho$ and $T$: \begin{equation} \begin{split} W &= \left(\frac{\partial p}{\partial T}\right)_{\rho} \\ X &= \left(\frac{\partial p}{\partial \rho}\right)_T \\ Y &= \left(\frac{\partial s}{\partial T}\right)_{\rho} \\ Z &= \left(\frac{\partial s}{\partial \rho}\right)_T \end{split} \end{equation}

Introduce temporary names for some of the derivatives along the saturation line: \begin{equation} \begin{split} M &= \frac{d \rho}{d h} = \frac{{d \rho}/{dT}}{{dh}/{dT}} \\ N &= \frac{d s}{d h} = \frac{{ds}/{dT}}{{dh}/{dT}} \end{split} \end{equation}

Now the rest is not too hard, but with intermediate steps it is quite long.

First example: $d^2 \rho / dT^2$

The corresponding first derivative can be written in two ways: \begin{equation} \begin{split} \frac{d \rho}{dT} &= \left(\frac{\partial \rho}{\partial T}\right)_p + \left(\frac{\partial \rho}{\partial p}\right)_T \frac{d p}{dT} = A + B \frac{dp}{dT}\\ &= {\frac{dp}{dT} \left(\frac{\partial p}{\partial T}\right)_{\rho}} / {\left(\frac{\partial p}{\partial \rho}\right)_T} = {\frac{dp}{dT} W}/X \end{split} \end{equation} Both ways can be used as a starting point for the second derivative: \begin{equation} \begin{split} \frac{d^2 \rho}{dT^2} &= \frac{dA}{dT} + \frac{dB}{dT}\frac{dp}{dT} + B\frac{d^2p}{dT^2}\\ &= \frac{\left(\frac{d^2p}{dT^2} - \frac{dW}{dT}\right)X - \left(\frac{dp}{dT}-W\right)\frac{dX}{dT}}{X^2} \end{split} \end{equation} The first variant seems to be nicer, because the first deriv did not have a fraction in it.
The derivatives along the saturation line that appeared can be written as (again, using two different approaches): \begin{equation} \begin{split} \frac{dA}{dT} &= \left(\frac{\partial A}{\partial T}\right)_p + \left(\frac{\partial A}{\partial p}\right)_T \frac{d p}{dT}\\ \frac{dB}{dT} &= \left(\frac{\partial B}{\partial T}\right)_p + \left(\frac{\partial B}{\partial p}\right)_T \frac{d p}{dT}\\ \frac{dW}{dT} &= \left(\frac{\partial W}{\partial T}\right)_{\rho} + \left(\frac{\partial W}{\partial \rho}\right)_T \frac{d\rho}{dT}\\ \frac{dX}{dT} &= \left(\frac{\partial X}{\partial T}\right)_{\rho} + \left(\frac{\partial X}{\partial \rho}\right)_T \frac{d\rho}{dT}\\ \end{split} \end{equation} and then the partial derivatives of A, B, W and X wrt p and T can be written as: \begin{equation} \begin{split} \left(\frac{\partial A}{\partial T}\right)_p &= \left(\frac{\partial^2 \rho}{\partial T^2}\right)_p \\ \left(\frac{\partial A}{\partial p}\right)_T &= \left(\frac{\partial^2 \rho}{\partial p \partial T}\right) \\ \left(\frac{\partial B}{\partial T}\right)_p &= \left(\frac{\partial^2 \rho}{\partial T \partial p}\right) = \left(\frac{\partial A}{\partial p}\right)_T \\ \left(\frac{\partial B}{\partial p}\right)_T &= \left(\frac{\partial^2 \rho}{\partial p^2}\right)_T \\ \left(\frac{\partial W}{\partial T}\right)_{\rho} &= \left(\frac{\partial^2 p}{\partial T^2}\right)_{\rho} \\ \left(\frac{\partial W}{\partial \rho}\right)_T &= \left(\frac{\partial^2 p}{\partial \rho \partial T}\right) \\ \left(\frac{\partial X}{\partial T}\right)_{\rho} &= \left(\frac{\partial^2 p}{\partial T \partial \rho}\right) = \left(\frac{\partial W}{\partial \rho}\right)_T \\ \left(\frac{\partial X}{\partial \rho}\right)_T &= \left(\frac{\partial^2 p}{\partial {\rho}^2}\right)_T \end{split} \end{equation} So this time, the second way looks maybe nicer, because it uses partial derivatives wrt to density and temperature. On the other hand, all second partial derivatives are already implemented in Coolprop and the pT derivatives are internally rewritten as dT derivatives.

Second example: $d^2s/dT^2$

The corresponding first derivative can be written in two ways: \begin{equation} \begin{split} \frac{ds}{dT} &= \left(\frac{\partial s}{\partial T}\right)_p + \left(\frac{\partial s}{\partial p}\right)_T \frac{dp}{dT} = C + E\frac{dp}{dT}\\ &= \left(\frac{\partial s}{\partial T}\right)_{\rho} + \left(\frac{\partial s}{\partial \rho}\right)_T \frac{d \rho}{dT} = Y + Z\frac{d \rho}{dT} \end{split} \end{equation} Both can be used as starting point for the second derivatives. \begin{split} \frac{d^2 s}{dT^2} &= \frac{dC}{dT} + \frac{dE}{dT}\frac{dp}{dT} + E\frac{d^2p}{dT^2}\\ &= \frac{dY}{dT} + \frac{dZ}{dT}\frac{dp}{dT} + Z\frac{d^2 \rho}{dT^2} \end{split} Now, which one is nicer to work with? Not sure. Just as shown above, the derivatives along the saturation line that appeared can be written as (again, using two different approaches): \begin{equation} \begin{split} \frac{dC}{dT} &= \left(\frac{\partial C}{\partial T}\right)_p + \left(\frac{\partial C}{\partial p}\right)_T \frac{d p}{dT}\\ \frac{E}{dT} &= \left(\frac{\partial E}{\partial T}\right)_p + \left(\frac{\partial E}{\partial p}\right)_T \frac{d p}{dT}\\ \frac{dY}{dT} &= \left(\frac{\partial Y}{\partial T}\right)_{\rho} + \left(\frac{\partial Y}{\partial \rho}\right)_T \frac{d\rho}{dT}\\ \frac{dZ}{dT} &= \left(\frac{\partial Z}{\partial T}\right)_{\rho} + \left(\frac{\partial Z}{\partial \rho}\right)_T \frac{d\rho}{dT}\\ \end{split} \end{equation} Rewriting the the partial derivatives works just like shown above.

Third example: $d^2{\rho}/dp^2$

The corresponding first derivative can be written in many ways, here is one: \begin{equation} \frac{d \rho}{dp} = \left(\frac{\partial \rho}{\partial p}\right)_T + \left(\frac{\partial \rho}{\partial T}\right)_p \frac{dT}{dp} = B + A\frac{dT}{dp} \end{equation}

The second derivative is then: \begin{equation} \frac{d^2 \rho}{dp^2} = \frac{dB}{dp} + \frac{dA}{dp}\frac{dT}{dp} + A\frac{d^2T}{dp^2} \end{equation} The derivatives wrt p along the saturation line can be written as: \begin{equation} \begin{split} \frac{dB}{dp} &= \left(\frac{\partial B}{\partial p}\right)_T + \left(\frac{\partial B}{\partial T}\right)_p \frac{d T}{dp}\\ \frac{dA}{dp} &= \left(\frac{\partial A}{\partial p}\right)_T + \left(\frac{\partial A}{\partial T}\right)_p \frac{d T}{dp} \end{split} \end{equation} Rewriting the partial derivatives works just like shown above.

Fourth example: $d^2s/dp^2$

The corresponding first derivative can be written in many ways, here is one: \begin{equation} \frac{ds}{dp} = \left(\frac{\partial s}{\partial p}\right)_T + \left(\frac{\partial s}{\partial T}\right)_p \frac{dT}{dp} = E + C\frac{dT}{dp} \end{equation}

The second derivative is then: \begin{equation} \frac{d^2 s}{dp^2} = \frac{dE}{dp} + \frac{dC}{dp}\frac{dT}{dp} + C\frac{d^2T}{dp^2} \end{equation} Rewriting the derivatives along the sat line and the partial derivatives works just like shown above.

Fifth example: $d^2h/dT^2$

First: \begin{equation} \frac{dh}{dT} = \left(\frac{\partial h}{\partial T}\right)_p + \left(\frac{\partial h}{\partial p}\right)_T \frac{dp}{dT} = G + K\frac{dp}{dT} \end{equation}

Second: \begin{equation} \frac{d^2 h}{dT^2} = \frac{dG}{dT} + \frac{dK}{dT}\frac{dp}{dT} + K\frac{d^2p}{dT^2} \end{equation}

Sixth example: $d^2{\rho}/dh^2$

This is a bit different, but also not difficult.
First: \begin{equation} M = \frac{d\rho}{dh} = \frac{{d\rho}/{dT}}{{dh}/{dT}} \end{equation}

Intermediate: \begin{equation} \frac{dM}{dT} = \frac{\frac{d^2 \rho}{dT^2}\frac{dh}{dT}-\frac{d \rho}{dT}\frac{d^2h}{dT^2}}{\left(\frac{dh}{dT}\right)^2} \end{equation}

Second: \begin{equation} \frac{d^2 \rho}{dh^2} = \frac{dM}{dh} = \frac{{dM}/{dT}}{{dh}/{dT}} = \frac{\frac{d^2 \rho}{dT^2}\frac{dh}{dT}-\frac{d \rho}{dT}\frac{d^2h}{dT^2}}{\left(\frac{dh}{dT}\right)^3} \end{equation}

Seventh example: $d^2s/dh^2$

First: \begin{equation} N = \frac{ds}{dh} = \frac{{ds}/{dT}}{{dh}/{dT}} \end{equation}

Intermediate: \begin{equation} \frac{dN}{dT} = \frac{\frac{d^2s}{dT^2}\frac{dh}{dT}-\frac{ds}{dT}\frac{d^2h}{dT^2}}{\left(\frac{dh}{dT}\right)^2} \end{equation}

Second: \begin{equation} \frac{d^2s}{dh^2} = \frac{dN}{dh} = \frac{{dN}/{dT}}{{dh}/{dT}} = \frac{\frac{d^2s}{dT^2}\frac{dh}{dT}-\frac{ds}{dT}\frac{d^2h}{dT^2}}{\left(\frac{dh}{dT}\right)^3} \end{equation}

Eigth example: $d^2s/(dh dp)$

First: \begin{equation} N = \frac{ds}{dh} = \frac{{ds}/{dT}}{{dh}/{dT}} \end{equation}

Intermediate: \begin{equation} \frac{dN}{dT} = \frac{\frac{d^2s}{dT^2}\frac{dh}{dT}-\frac{ds}{dT}\frac{d^2h}{dT^2}}{\left(\frac{dh}{dT}\right)^2} \end{equation}

Second: \begin{equation} \frac{d^2s}{dh dp} = \frac{dN}{dp} = \frac{{dN}/{dT}}{{dp}/{dT}} = \frac{\frac{d^2s}{dT^2}\frac{dh}{dT}-\frac{ds}{dT}\frac{d^2h}{dT^2}}{\left(\frac{dh}{dT}\right)^2 \left(\frac{dp}{dT}\right)} \end{equation}



In [6]:

    
# load some bits and pieces
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

import CoolProp as CP
from CoolProp.CoolProp import PropsSI

# Check: CoolProp version
print(CP.__version__)
print(CP.__gitrevision__)

# Constants
eps = 1e-9
kilo = 1e3
Mega = 1e6
golden = (1 + 5 ** 0.5) / 2
nPoints = 1000
width = 12.5

# helper function: get slope from two sorted arrays
def numSlopeAr(xAr, yAr):
 deltaX = np.ones(len(xAr))
 deltaY = np.ones(len(yAr))
 slopeAr = np.ones(len(yAr))
 for index in range(1,len(xAr)-1):
     deltaX[index] = xAr[index-1] - xAr[index+1]
     deltaY[index] = yAr[index-1] - yAr[index+1]
     slopeAr[index] = deltaY[index]/deltaX[index]
 # inaccurate, but who cares?
 slopeAr[0]=slopeAr[1]
 slopeAr[-1]=slopeAr[-2]  
 return slopeAr

# helper function: draw tangent from array
def drawTangent(xAr, yAr, slopeAr, index: int):
 xVal = xAr[index]
 yVal = yAr[index]
 slopeVal = slopeAr[index]
 # line length as percentage of length of x axis
 xRange = 0.15*abs(max(xAr)-min(xAr))
 xRangeAr = np.arange(xVal-xRange, xVal+xRange)
 tang = yVal + slopeVal*(xRangeAr-xVal)
 plt.plot(xVal,yVal,'om',xRangeAr,tang,'-k')









    



5.1.1
12f006445f234e572e64cc820146ab5d2c2a9d10



In [7]:

    
# All calculations happen in this cell
# All following cells are used for plotting only
# Run this cell twice to get rid of error message

# Set FluidName
FluidName = 'Propane'

# Triple and critical data
T_crt = PropsSI('TCRIT',FluidName)
T_trp = PropsSI('TTRIPLE',FluidName)
p_crt = PropsSI('PCRIT',FluidName)
p_trp = PropsSI('PTRIPLE',FluidName)
d_crt = PropsSI('RHOCRIT',FluidName)
d_trp_liq = PropsSI('D','T',T_trp,'Q',0,FluidName)
d_trp_vap  = PropsSI('D','T',T_trp,'Q',1,FluidName)
print("T_crt = " + str(T_crt))
print("T_trp = " + str(T_trp))

# Properties at saturation, liq and vap
# All the way to the crt point or keep some distance?
T_sat = np.linspace(T_trp, T_crt, num=nPoints)
p_sat = CP.CoolProp.PropsSI('P','T',T_sat,'Q',0,FluidName)
d_sat_liq = CP.CoolProp.PropsSI('D','T',T_sat,'Q',0,FluidName)
d_sat_vap = CP.CoolProp.PropsSI('D','T',T_sat,'Q',1,FluidName)
v_sat_liq = 1/d_sat_liq
v_sat_vap = 1/d_sat_vap
s_sat_liq = CP.CoolProp.PropsSI('S','T',T_sat,'Q',0,FluidName)
s_sat_vap = CP.CoolProp.PropsSI('S','T',T_sat,'Q',1,FluidName)
h_sat_liq = CP.CoolProp.PropsSI('H','T',T_sat,'Q',0,FluidName)
h_sat_vap = CP.CoolProp.PropsSI('H','T',T_sat,'Q',1,FluidName)
u_sat_liq = CP.CoolProp.PropsSI('U','T',T_sat,'Q',0,FluidName)
u_sat_vap = CP.CoolProp.PropsSI('U','T',T_sat,'Q',1,FluidName)

# Clausius-Clapeyron
# at crt point, vap=liq, so D becomes zero
Ds_vl = (s_sat_vap-s_sat_liq)
Dv_vl = (v_sat_vap-v_sat_liq)
dp_dT_q =  Ds_vl/Dv_vl 
dT_dp_q = Dv_vl/Ds_vl # = 1/dp_dT_q
dp_dT_qn = numSlopeAr(T_sat,p_sat)
dT_dp_qn = numSlopeAr(p_sat,T_sat)

# derivs wrt pT, single-phase AT the sat liq line
dd_dp_Tl = CP.CoolProp.PropsSI('d(Dmass)/d(P)|T','D',d_sat_liq,'T',T_sat,FluidName)
dd_dT_pl = CP.CoolProp.PropsSI('d(Dmass)/d(T)|P','D',d_sat_liq,'T',T_sat,FluidName)
ds_dp_Tl = CP.CoolProp.PropsSI('d(Smass)/d(P)|T','D',d_sat_liq,'T',T_sat,FluidName)
ds_dT_pl = CP.CoolProp.PropsSI('d(Smass)/d(T)|P','D',d_sat_liq,'T',T_sat,FluidName)
dh_dp_Tl = CP.CoolProp.PropsSI('d(Hmass)/d(P)|T','D',d_sat_liq,'T',T_sat,FluidName)
dh_dT_pl = CP.CoolProp.PropsSI('d(Hmass)/d(T)|P','D',d_sat_liq,'T',T_sat,FluidName)
dv_dp_Tl = -dd_dp_Tl/d_sat_liq**2
dv_dT_pl = -dd_dT_pl/d_sat_liq**2

# derivs wrt pT, single-phase AT the sat vap line
dd_dp_Tv = CP.CoolProp.PropsSI('d(Dmass)/d(P)|T','D',d_sat_vap,'T',T_sat,FluidName)
dd_dT_pv = CP.CoolProp.PropsSI('d(Dmass)/d(T)|P','D',d_sat_vap,'T',T_sat,FluidName)
ds_dp_Tv = CP.CoolProp.PropsSI('d(Smass)/d(P)|T','D',d_sat_vap,'T',T_sat,FluidName)
ds_dT_pv = CP.CoolProp.PropsSI('d(Smass)/d(T)|P','D',d_sat_vap,'T',T_sat,FluidName)
dh_dp_Tv = CP.CoolProp.PropsSI('d(Hmass)/d(P)|T','D',d_sat_vap,'T',T_sat,FluidName)
dh_dT_pv = CP.CoolProp.PropsSI('d(Hmass)/d(T)|P','D',d_sat_vap,'T',T_sat,FluidName)
dv_dp_Tv = -dd_dp_Tv/d_sat_vap**2
dv_dT_pv = -dd_dT_pv/d_sat_vap**2

# derivs wrt dT, single-phase AT the sat liq line
dp_dd_Tl = CP.CoolProp.PropsSI('d(P)/d(Dmass)|T','D',d_sat_liq,'T',T_sat,FluidName)
dp_dT_dl = CP.CoolProp.PropsSI('d(P)/d(T)|Dmass','D',d_sat_liq,'T',T_sat,FluidName)
ds_dd_Tl = CP.CoolProp.PropsSI('d(Smass)/d(D)|T','D',d_sat_liq,'T',T_sat,FluidName)
ds_dT_dl = CP.CoolProp.PropsSI('d(Smass)/d(T)|D','D',d_sat_liq,'T',T_sat,FluidName)
dh_dd_Tl = CP.CoolProp.PropsSI('d(Hmass)/d(D)|T','D',d_sat_liq,'T',T_sat,FluidName)
dh_dT_dl = CP.CoolProp.PropsSI('d(Hmass)/d(T)|D','D',d_sat_liq,'T',T_sat,FluidName)
dp_dv_Tl = -d_sat_liq**2 * dp_dd_Tl
dp_dT_vl = dp_dT_dl

# derivs wrt dT, single-phase AT the sat vap line
dp_dd_Tv = CP.CoolProp.PropsSI('d(P)/d(Dmass)|T','D',d_sat_vap,'T',T_sat,FluidName)
dp_dT_dv = CP.CoolProp.PropsSI('d(P)/d(T)|Dmass','D',d_sat_vap,'T',T_sat,FluidName)
ds_dd_Tv = CP.CoolProp.PropsSI('d(Smass)/d(D)|T','D',d_sat_vap,'T',T_sat,FluidName)
ds_dT_dv = CP.CoolProp.PropsSI('d(Smass)/d(T)|D','D',d_sat_vap,'T',T_sat,FluidName)
dh_dd_Tv = CP.CoolProp.PropsSI('d(Hmass)/d(D)|T','D',d_sat_vap,'T',T_sat,FluidName)
dh_dT_dv = CP.CoolProp.PropsSI('d(Hmass)/d(T)|D','D',d_sat_vap,'T',T_sat,FluidName)
dp_dv_Tv = -d_sat_vap**2 * dp_dd_Tv
dp_dT_vv = dp_dT_dv

# two ways for derivs along the sat line: using derivs wrt dT or pT

# derivs wrt T ALONG the sat liq line (that is, at q=const=0)
dd_dT_ql = (dp_dT_q - dp_dT_dl)/dp_dd_Tl
dv_dT_ql = dv_dT_pl + dv_dp_Tl*dp_dT_q
ds_dT_ql = ds_dT_dl + ds_dd_Tl*dd_dT_ql
#ds_dT_ql = ds_dT_pl + ds_dp_Tl*dp_dT_q
dh_dT_ql = dh_dT_pl + dh_dp_Tl*dp_dT_q

# derivs wrt T ALONG the sat vap line (that is, at q=const=1)
dd_dT_qv = (dp_dT_q - dp_dT_dv)/dp_dd_Tv
dv_dT_qv = (dp_dT_q - dp_dT_dv)/dp_dd_Tv
dv_dT_qv = dv_dT_pv + dv_dp_Tv*dp_dT_q
ds_dT_qv = ds_dT_dv + ds_dd_Tv*dd_dT_qv
dh_dT_qv = dh_dT_pv + dh_dp_Tv*dp_dT_q

# derivs wrt p ALONG the sat liq line (that is, at q=const=0)
dd_dp_ql = dd_dp_Tl + dd_dT_pl*dT_dp_q
ds_dp_ql = ds_dp_Tl + ds_dT_pl*dT_dp_q
dv_dp_ql = dv_dp_Tl + dv_dT_pl*dT_dp_q

# derivs wrt p ALONG the sat vap line (that is, at q=const=1)
dd_dp_qv = dd_dp_Tv + dd_dT_pv*dT_dp_q
ds_dp_qv = ds_dp_Tv + ds_dT_pv*dT_dp_q
dv_dp_qv = dv_dp_Tv + dv_dT_pv*dT_dp_q

# derivs wrt h ALONG the sat line
dd_dh_ql = dd_dT_ql/dh_dT_ql
dd_dh_qv = dd_dT_qv/dh_dT_qv
ds_dh_ql = ds_dT_ql/dh_dT_ql
ds_dh_qv = ds_dT_qv/dh_dT_qv

# derivs of Clausius-Clapeyron
# d2p_dT2_q = ((ds_dT_qv-ds_dT_ql)*Dv_vl - Ds_vl*(dv_dT_qv-dv_dT_ql)) / Dv_vl**2
d2p_dT2_q = (ds_dT_qv-ds_dT_ql)/Dv_vl - dp_dT_q*(dv_dT_qv-dv_dT_ql)/Dv_vl
d2T_dp2_q = (dv_dp_qv-dv_dp_ql)/Ds_vl - dT_dp_q*(ds_dp_qv-ds_dp_ql)/Ds_vl
d2p_dT2_qn = numSlopeAr(T_sat,dp_dT_q)
d2T_dp2_qn = numSlopeAr(p_sat,dT_dp_q)

# second derivs AT the sat line
# AT the liq line
# wrt pT 
d2d_dT2_pl = np.empty(nPoints)
d2d_dp2_Tl = np.empty(len(T_sat))
d2d_dpTl   = np.empty(len(T_sat))
d2s_dT2_pl = np.empty(len(T_sat))
d2s_dp2_Tl = np.empty(len(T_sat))
d2s_dpTl   = np.empty(len(T_sat))
d2h_dT2_pl = np.empty(len(T_sat))
d2h_dp2_Tl = np.empty(len(T_sat))
d2h_dpTl   = np.empty(len(T_sat))
# wrt dT
d2s_dT2_dl = np.empty(len(T_sat))
d2s_dd2_Tl = np.empty(len(T_sat))
d2s_ddTl   = np.empty(len(T_sat))
# AT the vap line
# wrt pT 
d2d_dT2_pv = np.empty(len(T_sat))
d2d_dp2_Tv = np.empty(len(T_sat))
d2d_dpTv   = np.empty(len(T_sat))
d2s_dT2_pv = np.empty(len(T_sat))
d2s_dp2_Tv = np.empty(len(T_sat))
d2s_dpTv   = np.empty(len(T_sat))
d2h_dT2_pv = np.empty(len(T_sat))
d2h_dp2_Tv = np.empty(len(T_sat))
d2h_dpTv   = np.empty(len(T_sat))
# wrt dT
d2s_dT2_dv = np.empty(len(T_sat))
d2s_dd2_Tv = np.empty(len(T_sat))
d2s_ddTv   = np.empty(len(T_sat))

HEOS = CP.AbstractState("HEOS", FluidName)
for idx in range(0,len(T_sat)):
    # AT the liq line
    HEOS.update(CP.QT_INPUTS, 0, T_sat[idx])    
    # wrt pT
    d2d_dT2_pl[idx] = HEOS.second_partial_deriv(CP.iDmass, CP.iT, CP.iP, CP.iT, CP.iP)
    d2d_dp2_Tl[idx] = HEOS.second_partial_deriv(CP.iDmass, CP.iP, CP.iT, CP.iP, CP.iT)
    d2d_dpTl[idx] = HEOS.second_partial_deriv(CP.iDmass, CP.iP, CP.iT, CP.iT, CP.iP)
    d2s_dT2_pl[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iT, CP.iP, CP.iT, CP.iP)
    d2s_dp2_Tl[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iP, CP.iT, CP.iP, CP.iT)
    d2s_dpTl[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iP, CP.iT, CP.iT, CP.iP)
    d2h_dT2_pl[idx] = HEOS.second_partial_deriv(CP.iHmass, CP.iT, CP.iP, CP.iT, CP.iP)
    d2h_dp2_Tl[idx] = HEOS.second_partial_deriv(CP.iHmass, CP.iP, CP.iT, CP.iP, CP.iT)
    d2h_dpTl[idx] = HEOS.second_partial_deriv(CP.iHmass, CP.iP, CP.iT, CP.iT, CP.iP)
    # wrt dT
    d2s_dT2_dl[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iT, CP.iDmass, CP.iT, CP.iDmass)
    d2s_dd2_Tl[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iDmass, CP.iT, CP.iDmass, CP.iT)
    d2s_ddTl[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iDmass, CP.iT, CP.iT, CP.iDmass)
    
    # AT the vap line
    HEOS.update(CP.QT_INPUTS, 1, T_sat[idx])    
    # wrt pT
    d2d_dT2_pv[idx] = HEOS.second_partial_deriv(CP.iDmass, CP.iT, CP.iP, CP.iT, CP.iP)
    d2d_dp2_Tv[idx] = HEOS.second_partial_deriv(CP.iDmass, CP.iP, CP.iT, CP.iP, CP.iT)
    d2d_dpTv[idx] = HEOS.second_partial_deriv(CP.iDmass, CP.iP, CP.iT, CP.iT, CP.iP)
    d2s_dT2_pv[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iT, CP.iP, CP.iT, CP.iP)
    d2s_dp2_Tv[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iP, CP.iT, CP.iP, CP.iT)
    d2s_dpTv[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iP, CP.iT, CP.iT, CP.iP)
    d2h_dT2_pv[idx] = HEOS.second_partial_deriv(CP.iHmass, CP.iT, CP.iP, CP.iT, CP.iP)
    d2h_dp2_Tv[idx] = HEOS.second_partial_deriv(CP.iHmass, CP.iP, CP.iT, CP.iP, CP.iT)
    d2h_dpTv[idx] = HEOS.second_partial_deriv(CP.iHmass, CP.iP, CP.iT, CP.iT, CP.iP)
    # wrt dT
    d2s_dT2_dv[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iT, CP.iDmass, CP.iT, CP.iDmass)
    d2s_dd2_Tv[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iDmass, CP.iT, CP.iDmass, CP.iT)
    d2s_ddTv[idx] = HEOS.second_partial_deriv(CP.iSmass, CP.iDmass, CP.iT, CP.iT, CP.iDmass)    
    

# calculate 2nd derivs ALONG the sat line (analytically and numerically)

# A=dd_dT_p and B=dd_dp_T
# liq side
dA_dT_ql = d2d_dT2_pl + d2d_dpTl*dp_dT_q
dB_dT_ql = d2d_dpTl + d2d_dp2_Tl*dp_dT_q
d2d_dT2_qla = dA_dT_ql + dB_dT_ql*dp_dT_q + dd_dp_Tl*d2p_dT2_q
d2d_dT2_qln = numSlopeAr(T_sat, dd_dT_ql)
# vap side
dA_dT_qv = d2d_dT2_pv + d2d_dpTv*dp_dT_q
dB_dT_qv = d2d_dpTv + d2d_dp2_Tv*dp_dT_q
d2d_dT2_qva = dA_dT_qv + dB_dT_qv*dp_dT_q + dd_dp_Tv*d2p_dT2_q
d2d_dT2_qvn = numSlopeAr(T_sat, dd_dT_qv)
#print(d2d_dT2_qla - d2d_dT2_qln)
#print(d2d_dT2_qva - d2d_dT2_qln)

# C=ds_dT_p and E=ds_dp_T
# liq side
dC_dT_ql = d2s_dT2_dl + d2s_ddTl*dd_dT_ql
dE_dT_ql = d2s_ddTl + d2s_dd2_Tl*dd_dT_ql
d2s_dT2_qla = dC_dT_ql + dE_dT_ql*dd_dT_ql + ds_dd_Tl*d2d_dT2_qla
d2s_dT2_qln = numSlopeAr(T_sat, ds_dT_ql)
# vap side
dC_dT_qv = d2s_dT2_dv + d2s_ddTv*dd_dT_qv
dE_dT_qv = d2s_ddTv + d2s_dd2_Tv*dd_dT_qv
d2s_dT2_qva = dC_dT_qv + dE_dT_qv*dd_dT_qv + ds_dd_Tv*d2d_dT2_qva
d2s_dT2_qvn = numSlopeAr(T_sat, ds_dT_qv)
#print(d2s_dT2_qla - d2s_dT2_qln)
#print(d2s_dT2_qva - d2s_dT2_qvn)

# B=dd_dp_T and A=dd_dT_p
# liq side
dB_dp_ql = d2d_dp2_Tl + d2d_dpTl*dT_dp_q
dA_dp_ql = d2d_dpTl + d2d_dT2_pl*dT_dp_q
d2d_dp2_qla = dB_dp_ql + dA_dp_ql*dT_dp_q + dd_dT_pl*d2T_dp2_q
d2d_dp2_qln = numSlopeAr(T_sat, ds_dp_ql)
# vap side
dB_dp_qv = d2d_dp2_Tv + d2d_dpTv*dT_dp_q
dA_dp_qv = d2d_dpTv + d2d_dT2_pv*dT_dp_q
d2d_dp2_qva = dB_dp_qv + dA_dp_qv*dT_dp_q + dd_dT_pv*d2T_dp2_q
d2d_dp2_qvn = numSlopeAr(T_sat, ds_dp_qv)
#print(d2d_dp2_qla - d2d_dp2_qln)
#print(d2d_dp2_qva - d2d_dp2_qvn)

# E=ds_dp_T and C=ds_dT_p
# liq side
dE_dp_ql = d2s_dp2_Tl + d2s_dpTl*dT_dp_q
dC_dp_ql = d2s_dpTl + d2s_dT2_pl*dT_dp_q
d2s_dp2_qla = dE_dp_ql + dC_dp_ql*dT_dp_q + ds_dT_pl*d2T_dp2_q
d2s_dp2_qln = numSlopeAr(p_sat, ds_dp_ql)
# vap side
dE_dp_qv = d2s_dp2_Tv + d2s_dpTv*dT_dp_q
dC_dp_qv = d2s_dpTv + d2s_dT2_pv*dT_dp_q
d2s_dp2_qva = dE_dp_qv + dC_dp_qv*dT_dp_q + ds_dT_pv*d2T_dp2_q
d2s_dp2_qvn = numSlopeAr(p_sat, ds_dp_qv)
#print(d2s_dp2_qla - d2s_dp2_qln)
#print(d2s_dp2_qva - d2s_dp2_qvn)

# G=dh_dT_p and K=dh_dp_T
# liq side
dG_dT_ql = d2h_dT2_pl + d2h_dpTl*dp_dT_q
dK_dT_ql = d2h_dpTl + d2h_dp2_Tl*dp_dT_q
d2h_dT2_qla = dG_dT_ql + dK_dT_ql*dp_dT_q + dh_dp_Tl*d2p_dT2_q
d2h_dT2_qln = numSlopeAr(T_sat, dh_dT_ql)
# vap side
dG_dT_qv = d2h_dT2_pv + d2h_dpTv*dp_dT_q
dK_dT_qv = d2h_dpTv + d2h_dp2_Tv*dp_dT_q
d2h_dT2_qva = dG_dT_qv + dK_dT_qv*dp_dT_q + dh_dp_Tv*d2p_dT2_q
d2h_dT2_qvn = numSlopeAr(T_sat, dh_dT_qv)
#print(d2h_dT2_qla - d2h_dT2_qln)
#print(d2h_dT2_qva - d2h_dT2_qvn)

# liq side
d2d_dh2_qla = (d2d_dT2_qla*dh_dT_ql - dd_dT_ql*d2h_dT2_qla) / dh_dT_ql**3
d2d_dh2_qln = numSlopeAr(h_sat_liq, dd_dh_ql)
# vap side
d2d_dh2_qva = (d2d_dT2_qva*dh_dT_qv - dd_dT_qv*d2h_dT2_qva) / dh_dT_qv**3
d2d_dh2_qvn = numSlopeAr(h_sat_vap, dd_dh_qv)
#print(d2d_dh2_qla - d2d_dh2_qln)
#print(d2d_dh2_qva - d2d_dh2_qvn)

# liq side
d2s_dh2_qla = (d2s_dT2_qla*dh_dT_ql - ds_dT_ql*d2h_dT2_qla) / dh_dT_ql**3
d2s_dh2_qln = numSlopeAr(h_sat_liq, ds_dh_ql)
# vap side
d2s_dh2_qva = (d2s_dT2_qva*dh_dT_qv - ds_dT_qv*d2h_dT2_qva) / dh_dT_qv**3
d2s_dh2_qvn = numSlopeAr(h_sat_vap, ds_dh_qv)
#print(d2s_dh2_qla - d2s_dh2_qln)
#print(d2s_dh2_qva - d2s_dh2_qvn)









    



T_crt = 369.89
T_trp = 85.525



In [8]:

    
# vapor pressure, Clausius-Clapeyron

# pick a point and set figure size
mySatPoint = int(nPoints-50)
print("T_sat = " + str(T_sat[mySatPoint]))
print("p_sat = " + str(p_sat[mySatPoint]))
plt.figure(figsize=(width,width*3/2/golden))

plt.subplot(3,2,1)
plt.plot(T_sat, p_sat, color='blue')
#plt.yscale('log')
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.grid(b=True, linestyle=':')
plt.minorticks_on()
plt.xlabel('Temperature in K')
plt.ylabel('Pressure in Pa')
drawTangent(T_sat, p_sat, dp_dT_q, mySatPoint)
#print(dp_dT_q - dp_dT_qn)

plt.subplot(3,2,2)
plt.plot(p_sat, T_sat, color='blue')
#plt.yscale('log')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.grid(b=True, linestyle=':')
plt.minorticks_on()
plt.xlabel('Pressure in Pa')
plt.ylabel('Temperature in K')
drawTangent(p_sat, T_sat, dT_dp_q, mySatPoint)
#print(dT_dp_q - dT_dp_qn)

plt.subplot(3,2,3)
plt.plot(T_sat, dp_dT_q, color='blue')
#plt.yscale('log')
plt.grid(b=True, linestyle=':')
plt.minorticks_on()
plt.xlabel('Temperature in K')
plt.ylabel('dp/dT in Pa/K')
drawTangent(T_sat, dp_dT_q, d2p_dT2_q, mySatPoint)
#print(d2p_dT2_q - d2p_dT2_qn)

plt.subplot(3,2,4)
plt.plot(p_sat, dT_dp_q, color='blue')
plt.yscale('log')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.grid(b=True, linestyle=':')
plt.minorticks_on()
plt.xlabel('Pressure in Pa')
plt.ylabel('dT/dp in K/Pa')
drawTangent(p_sat, dT_dp_q, d2T_dp2_q, mySatPoint)
#print(d2T_dp2_q - d2T_dp2_qn)


plt.subplot(3,2,5)
plt.plot(T_sat, d2p_dT2_q, color='blue')
plt.grid(b=True, linestyle=':')
plt.minorticks_on()
plt.xlabel('Temperature in K')
plt.ylabel('d2p/dT2 in Pa/K2')

plt.subplot(3,2,6)
plt.plot(p_sat, abs(d2T_dp2_q), color='blue')
plt.yscale('log')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.grid(b=True, linestyle=':')
plt.minorticks_on()
plt.xlabel('Pressure in Pa')
plt.ylabel('|d2T/dp2| in K/Pa2')









    



T_sat = 355.942167167
p_sat = 3298979.183






    Out[8]:





<matplotlib.text.Text at 0xda6d748>



In [9]:

    
# Pick a point and set figure size
mySatPoint = int(nPoints-50)
print("T_sat = " + str(T_sat[mySatPoint]))
print("p_sat = " + str(p_sat[mySatPoint]))
plt.figure(figsize=(width,width*3/2/golden))

# d versus T
plt.subplot(3,2,1)
plt.plot(T_sat, d_sat_liq, color='blue')
plt.plot(T_sat, d_sat_vap, color='red')
#plt.plot(T_sat, (d_sat_liq+d_sat_vap)/2, color='black')
plt.xlabel('Temperature in K')
plt.ylabel('Density in kg/m³')
drawTangent(T_sat, d_sat_liq, dd_dT_ql, mySatPoint)
drawTangent(T_sat, d_sat_vap, dd_dT_qv, mySatPoint)

# s versus T
plt.subplot(3,2,2)
plt.plot(T_sat, s_sat_liq, color='blue')
plt.plot(T_sat, s_sat_vap, color='red')
plt.xlabel('Temperature in K')
plt.ylabel('Specific entropy in J/(kg·K)')
drawTangent(T_sat, s_sat_liq, ds_dT_ql, mySatPoint)
drawTangent(T_sat, s_sat_vap, ds_dT_qv, mySatPoint)

# d versus p
plt.subplot(3,2,3)
plt.plot(p_sat, d_sat_liq, color='blue')
plt.plot(p_sat, d_sat_vap, color='red')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel('Pressure in Pa')
plt.ylabel('Density in kg/m³')
drawTangent(p_sat, d_sat_liq, dd_dp_ql, mySatPoint)
drawTangent(p_sat, d_sat_vap, dd_dp_qv, mySatPoint)
#print(dd_dp_ql - numSlopeAr(p_sat, d_sat_liq))
#print(dd_dp_qv - numSlopeAr(p_sat, d_sat_vap))

# s versus p
plt.subplot(3,2,4)
plt.plot(p_sat, s_sat_liq, color='blue')
plt.plot(p_sat, s_sat_vap, color='red')
#plt.yscale('log')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel('Pressure in Pa')
plt.ylabel('Specific entropy in J/(kg·K)')
drawTangent(p_sat, s_sat_liq, ds_dp_ql, mySatPoint)
drawTangent(p_sat, s_sat_vap, ds_dp_qv, mySatPoint)
#print(ds_dp_ql - numSlopeAr(p_sat, s_sat_liq))
#print(ds_dp_qv - numSlopeAr(p_sat, s_sat_vap))

# d versus h
plt.subplot(3,2,5)
plt.plot(h_sat_liq, d_sat_liq, color='blue')
plt.plot(h_sat_vap, d_sat_vap, color='red')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel('Specific enthalpy in J/kg')
plt.ylabel('Density in kg/m³')
drawTangent(h_sat_liq, d_sat_liq, dd_dh_ql, mySatPoint)
drawTangent(h_sat_vap, d_sat_vap, dd_dh_qv, mySatPoint)

# s versus h
plt.subplot(3,2,6)
plt.plot(h_sat_liq, s_sat_liq, color='blue')
plt.plot(h_sat_vap, s_sat_vap, color='red')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel('Specific enthalpy in J/kg')
plt.ylabel('Specific entropy in J/(kg·K)')
drawTangent(h_sat_liq, s_sat_liq, ds_dh_ql, mySatPoint)
drawTangent(h_sat_vap, s_sat_vap, ds_dh_qv, mySatPoint)









    



T_sat = 355.942167167
p_sat = 3298979.183



In [10]:

    
# Pick a point and set figure size
mySatPoint = int(nPoints-50)
print("T_sat = " + str(T_sat[mySatPoint]))
print("p_sat = " + str(p_sat[mySatPoint]))
plt.figure(figsize=(width,width*3/2/golden))

# dd_dT versus T
plt.subplot(3,2,1)
plt.plot(T_sat, dd_dT_ql, color='blue')
plt.plot(T_sat, dd_dT_qv, color='red')
plt.xlabel('Temperature in K')
plt.ylabel('dd/dT in kg/m³/K')
drawTangent(T_sat, dd_dT_ql, d2d_dT2_qla, mySatPoint)
drawTangent(T_sat, dd_dT_qv, d2d_dT2_qva, mySatPoint)

# ds_dT versus T
plt.subplot(3,2,2)
plt.plot(T_sat, ds_dT_ql, color='blue')
plt.plot(T_sat, ds_dT_qv, color='red')
plt.xlabel('Temperature in K')
plt.ylabel('ds/dT in J//(kg·K)/K')
drawTangent(T_sat, ds_dT_ql, d2s_dT2_qla, mySatPoint)
drawTangent(T_sat, ds_dT_qv, d2s_dT2_qva, mySatPoint)

# dd_dp versus p
plt.subplot(3,2,3)
#plt.plot(p_sat, dd_dp_ql, color='blue')
plt.plot(p_sat, dd_dp_qv, color='red')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel('Pressure in Pa')
plt.ylabel('dd/dp in kg/m³/Pa')
#drawTangent(p_sat, dd_dp_ql, d2d_dp2_qla, mySatPoint)
drawTangent(p_sat, dd_dp_qv, d2d_dp2_qva, mySatPoint)

# ds_dp versus p
plt.subplot(3,2,4)
plt.plot(p_sat, ds_dp_ql, color='blue')
#plt.plot(p_sat, ds_dp_qv, color='red')
#plt.yscale('log')
plt.ylim([0,10*min(ds_dp_ql)])
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel('Pressure in Pa')
plt.ylabel('ds/dp in J//(kg·K)/Pa')
drawTangent(p_sat, ds_dp_ql, d2s_dp2_qla, mySatPoint)
#drawTangent(p_sat, ds_dp_qv, d2s_dp2_qva, mySatPoint)

# dd_dh versus h
plt.subplot(3,2,5)
plt.plot(h_sat_liq, dd_dh_ql, color='blue')
#plt.plot(h_sat_vap, dd_dh_qv, color='red')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel('Specific enthalpy in J/kg')
plt.ylabel('dd/dh')
drawTangent(h_sat_liq, dd_dh_ql, d2d_dh2_qla, mySatPoint)
#drawTangent(h_sat_vap, dd_dh_qv, d2d_dh2_qva, mySatPoint)

# ds_dh versus h
plt.subplot(3,2,6)
plt.plot(h_sat_liq, ds_dh_ql, color='blue')
#plt.plot(h_sat_vap, ds_dh_qv, color='red')
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel('Specific enthalpy in J/kg')
plt.ylabel('ds/dh')
drawTangent(h_sat_liq, ds_dh_ql, d2s_dh2_qla, mySatPoint)
#drawTangent(h_sat_vap, ds_dh_qv, d2s_dh2_qva, mySatPoint)









    



T_sat = 355.942167167
p_sat = 3298979.183