The analysis in this section follows the methodologies proposed in the GERG 2004 monograph published in 2007
System of Equations
Our residual vector $\mathbf{F}$ is equal to
$$F_i = \ln\phi(T,p,\mathbf{y})-\ln \phi(T,p,\mathbf{x})+\ln K_i=0, i=1,2,3... N$$$$F_{N+1} = \sum_{i=1}^{N}(y_i-x_i)=0$$where
$$x_i = \frac{z_i}{1-\beta+\beta K_i}$$and
$$y_i = \frac{K_iz_i}{1-\beta+\beta K_i}$$DO NOT NORMALIZE $x$ and $y$ !!!!
$$F_{N+2} = X_s - S = 0$$and the system to be solved is equal to
$$\mathbf{J}\mathbf{\Delta X}= -\mathbf{F}$$This is the trickiest part of this method. There are a lot of derivatives to implement, and we want to implement all of them analytically.
$$\frac{\partial F_i}{\partial \ln T} = T\left[ \left(\frac{\partial \ln \phi_i}{\partial T}\right)''_{p,\mathbf{n}} -\left(\frac{\partial \ln \phi_i}{\partial T}\right)'_{p,\mathbf{n}}\right]$$$$\frac{\partial F_i}{\partial \ln p} = p\left[ \left(\frac{\partial \ln \phi_i}{\partial p}\right)''_{T,\mathbf{n}} -\left(\frac{\partial \ln \phi_i}{\partial p}\right)'_{T,\mathbf{n}}\right]$$$$ \frac{\partial F_i}{\partial \ln K_j} = \frac{K_jz_j}{(1-\beta+\beta K_j)^2}[(1-\beta)\phi_{ij}''+\beta\phi_{ij}']+\zeta $$$\zeta$ is the Kronecker delta or $\zeta = 0$ for $i\neq j$ , and $\zeta = 0$ for $i=j$. Also
$$\phi_{ij} = n\left( \frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p}$$For the $F_{N+1}$ term,
$$\frac{\partial F_{N+1}}{\partial \ln K_j}=\frac{K_jz_j}{(1-\beta+\beta K_j)^2}$$and all other partials of $F_{N+1}$ in the Jacobian are zero. For the specified term
$$\frac{\partial F_{N+2}}{X_s}=1$$and all other partials of $F_{N+2}$ in the Jacobian are zero.
From GERG 2004 Monograph, Eqn 7.27:
$$\ln \phi_i = \left( \frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}-\ln Z$$and (Kunz, 2012, Table B4)
$$\left( \frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j} = \alpha^r + n\left( \frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}$$so
$$\ln \phi_i = \alpha^r + n\left( \frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}-\ln Z$$Another two alternatives have been proposed in A DENSITY MARCHING METHOD FOR CALCULATING PHASE ENVELOPES Gadhiraju Venkatarathnam, I&ECR, 2014
In this paper, density marching methods are proposed rather than methods that march in temperature, pressure, or K-factor.
The system of equations to be solved is similar to that of the GERG 2004 formulation, where the unknowns are $\ln(T)$, $\ln(p)$, and $\ln(K_i)$
(A1) - OK $$F_i = \ln K_i+\ln\phi(T,p,\mathbf{y})-\ln \phi(T,p,\mathbf{x})=0, i=1,2,3... N$$
(A2) - OK $$F_{N+1} = \sum_{i=1}^{N}\frac{z_i(K_i-1)}{1-\beta+\beta K_i}=0$$
(A3) - TYPO, should be $\ln(\rho)$ rather than $\rho$, and should be all on left-hand-side $$F_{N+2} = \ln \rho ''-\ln\rho''_{S} = 0$$
(A6) - TYPO, missing an n to multiply the terms $\left( \frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p}$ $$ \frac{\partial F_i}{\partial \ln K_j} = \frac{K_jz_j}{(1-\beta+\beta K_j)^2}[(1-\beta)\phi_{ij}''+\beta\phi_{ij}']+\zeta $$
$\zeta$ is the Kronecker delta or $\zeta = 0$ for $i\neq j$ , and $\zeta = 0$ for $i=j$.
(A7) - OK $$\frac{\partial F_i}{\partial \ln T} = T\left[ \left(\frac{\partial \ln \phi_i}{\partial T}\right)''_{p,\mathbf{n}} -\left(\frac{\partial \ln \phi_i}{\partial T}\right)'_{p,\mathbf{n}}\right]$$
(A8) - OK $$\frac{\partial F_i}{\partial \ln p} = p\left[ \left(\frac{\partial \ln \phi_i}{\partial p}\right)''_{T,\mathbf{n}} -\left(\frac{\partial \ln \phi_i}{\partial p}\right)'_{T,\mathbf{n}}\right]$$
(A9) - OK $$\frac{\partial F_{N+1}}{\partial \ln K_j}=\frac{K_jz_j}{(1-\beta+\beta K_j)^2}$$
(A11) - OK $$ \frac{\partial F_{N+2}}{\partial \ln K_j} = \frac{K_jz_j(1-\beta)\beta}{(1-\beta+\beta K_j)^2}\left(n\left(\frac{\partial \rho}{\partial n_j}\right)''_{T,p}\right)\left(\frac{1}{\rho''_{S}}\right)$$
(A12) - OK $$\frac{\partial F_{N+2}}{\partial \ln T}=\left(\frac{\partial \rho}{\partial T}\right)''_{p,n}\frac{T}{\rho''_{S}}$$
(A13) - OK $$\frac{\partial F_{N+2}}{\partial \ln p}=\left(\frac{\partial \rho}{\partial p}\right)''_{T,n}\frac{p}{\rho''_{S}}$$
Another two alternatives have been proposed in A DENSITY MARCHING METHOD FOR CALCULATING PHASE ENVELOPES Gadhiraju Venkatarathnam, I&ECR, 2014
In this paper, density marching methods are proposed rather than methods that march in temperature, pressure, or K-factor.
The system of equations to be solved is similar to that of the GERG 2004 formulation, where the unknowns are $\ln(T)$, $\ln(p)$, and $\ln(K_i)$
(A14) - OK $$F_i = \ln K_i+\ln\phi(T,\rho'',\mathbf{y})-\ln \phi(T,\rho',\mathbf{x})=0, i=1,2,3... N$$
(A15) - OK $$F_{N+1} = \sum_{i=1}^{N}\frac{z_i(K_i-1)}{1-\beta+\beta K_i}=0$$
(A16) - OK $$F_{N+2} = p(T,\rho',\mathbf{y})-p(T,\rho'',\mathbf{x}) = 0$$
(A17) - TYPO, missing an n to multiply the terms $\left( \frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p}$ $$ \frac{\partial F_i}{\partial \ln K_j} = \frac{K_jz_j}{(1-\beta+\beta K_j)^2}[(1-\beta)\phi_{ij}''+\beta\phi_{ij}']+\zeta $$
$\zeta$ is the Kronecker delta or $\zeta = 0$ for $i\neq j$ , and $\zeta = 0$ for $i=j$.
(A18) - OK $$\frac{\partial F_i}{\partial \ln T} = T\left[ \left(\frac{\partial \ln \phi_i}{\partial T}\right)''_{p,\mathbf{n}} -\left(\frac{\partial \ln \phi_i}{\partial T}\right)'_{p,\mathbf{n}}\right]$$
(A19) - OK $$\frac{\partial F_i}{\partial \ln \rho'} = -\rho'\left(\frac{\partial \ln \phi_i}{\partial \rho}\right)'_{T,n}$$
(A20) - OK $$\frac{\partial F_{N+1}}{\partial \ln K_j}=\frac{K_jz_j}{(1-\beta+\beta K_j)^2}$$
(A22) - TYPO Second derivative of ln(phi) with respect to rho' needs constraints, first needs to have the constraints in the right place $$ \frac{\partial F_{N+2}}{\partial \ln K_j} = \frac{RTK_jz_j}{(1-\beta+\beta K_j)^2}\left[(1-\beta)\left(\frac{\partial \ln \phi_i}{\partial \rho}\right)''_{T,n}+\beta\left(\frac{\partial \ln \phi_i}{\partial \rho}\right)'_{T,n}\right]$$
(A23) - TYPO Should be A23 $$\frac{\partial F_{N+2}}{\partial \ln T}=T\left[\left(\frac{\partial p}{\partial T}\right)'_{\rho',n}-\left(\frac{\partial p}{\partial T}\right)''_{\rho'',n} \right]$$
(A24) - TYPO Should be A24 $$\frac{\partial F_{N+2}}{\partial \ln p'}=\rho'\left(\frac{\partial p}{\partial \rho}\right)'_{T,n}$$
Three analyic derivatives are not provided in GERG and need to be rederived:
$$\left(\frac{\partial \ln \phi_i}{\partial \rho}\right)_{T,n}$$$$\left(\frac{\partial \ln \phi_i}{\partial T}\right)_{\rho,n}$$$$n\left(\frac{\partial \rho}{\partial n_j}\right)_{T,p}$$The last is for T,p marching, the first two are for density marching.
GERG 2007 Monograph Equation 7.32 gives
$$\left( \frac{\partial V}{\partial n_i}\right)_{T,p,n_j} = \dfrac{-\left(\dfrac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\dfrac{\partial p}{\partial V}\right)_{T,n}}$$expand the left hand side with
$V = vn = \dfrac{n}{\rho}$
n held constant in derivative, so get
$$\left( \frac{\partial V}{\partial n_i}\right)_{T,p,n_j} = n\left( \frac{\partial (1/\rho)}{\partial n_i}\right)_{T,p,n_j} = -\frac{n}{\rho^2}\left( \frac{\partial \rho}{\partial n_i}\right)_{T,p,n_j}$$so
$$ n\left( \frac{\partial \rho}{\partial n_i}\right)_{T,p,n_j} = -\rho^2\left( \frac{\partial V}{\partial n_i}\right)_{T,p,n_j} $$GERG 2007 Monograph 7.27 $$\ln \phi_i = \left(\frac{\partial n \alpha^r}{\partial n_i} \right)_{T,V,n_j} - \ln Z $$
GERG 2007 Monograph 7.34 $$\ln \left(\frac{f_i}{n_i}\right) = \ln\left(\frac{RT}{V}\right)+\left(\frac{\partial n\alpha^r}{\partial n_i} \right)_{T,V,n_j}$$
and $Z = (pV)/(nRT)$, thus
$$\ln \phi_i = \ln \left(\frac{f_i}{n_i}\right) - \ln\left(\frac{RT}{V}\right) - \ln \left(\frac{pV}{nRT}\right) ?= \ln \left(\frac{f_i}{n_i}\right) -\ln\left(\frac{p}{n}\right)$$$$\frac{\partial \left[\ln\left(\frac{p}{n}\right)\right]}{\partial T} = \frac{1}{pn}\left(\frac{\partial p}{\partial T}\right)_{\rho,n}$$Equivalent fugacities for the $i$-th component
$$ F_k^A = \ln f_i(T,p,\mathbf{x})-\ln f_i(T,p,\mathbf{y}) = 0\mbox{ for } k = i = 1...N $$Material balance
$$ F_k^B = \frac{z_i-x_i}{y_i-x_i}-\frac{z_{N-1}-x_{N-1}}{y_{N-1}-x_{N-1}}\mbox{ for }i=1..N-2; k = i+N; k = N+1..2N-2$$The independent variables to be obtained are
$$
To convert partial with $T$, $V$, $x_k$ held constant to one with $\tau$, $\delta$, $x_k$ held constant, use Gernert 3.118, or
$$ \frac{\partial}{\partial x_j} [Y]_{T,V,x_k} = \frac{\partial}{\partial x_j} [Y]_{\tau,\delta,x_k}+\left(\frac{\partial\delta}{\partial x_j}\right)_{T,V,x_k}\left.\frac{\partial Y}{\partial\delta}\right|_{\tau,\bar x}+\left(\frac{\partial\tau}{\partial x_j}\right)_{T,V,x_k}\left.\frac{\partial Y}{\partial\tau}\right|_{\delta,\bar x} $$To convert pressure,
$$ p=\rho R T(1+\delta \alpha_\delta) $$convert $\rho$ and $T$ to reduced variables
$$ p=\rho_r(\bar x)\delta R \frac{T_r(\bar x)}{\tau}(1+\delta \alpha_\delta) = \rho_r(\bar x)R \frac{T_r(\bar x)}{\tau}\delta (1+\delta \alpha_\delta)$$All the derivatives
$$ \frac{dp}{d\tau}\times\frac{d\tau}{dx_j} = -\rho_r(\bar x)\delta R \frac{T_r(\bar x)}{\tau^2}(1+\delta \alpha_\delta) \times \frac{1}{T}\frac{\partial T_r}{\partial x_j}|_{T,V,x_k}$$$$ \frac{dp}{d\delta}\times\frac{d\delta}{dx_j} = \rho_r(\bar x) R \frac{T_r(\bar x)}{\tau}[ (1+\delta \alpha_\delta)+ \delta(\alpha_\delta+\delta \alpha_{\delta\delta})] \times \frac{-\delta}{\rho_r}\frac{\partial \rho_r}{\partial x_j}|_{T,V,x_k}$$$$ \frac{\partial p}{\partial x_j}|_{\tau,\delta,x_k} = \frac{\delta R}{\tau}\left[\rho_r(\bar x)T_r(\bar x)(\delta \frac{\partial}{\partial x_j}[\alpha_\delta]_{\tau, \delta, x_k})+(1+\delta \alpha_\delta)\left(\rho_r \frac{\partial T_r}{\partial x_j}+T_r \frac{\partial \rho_r}{\partial x_j}\right)\right]$$with $\delta$ factored out
Turn back into normal variables $$ \frac{dp}{d\tau}\times\frac{d\tau}{dx_j} = -\frac{\rho R}{\tau}(1+\delta \alpha_\delta)\frac{\partial T_r}{\partial x_j}|_{T,V,x_k}$$
$$ \frac{dp}{d\delta}\times\frac{d\delta}{dx_j} = - R T[ (1+\delta \alpha_\delta)+ \delta(\alpha_\delta+\delta \alpha_{\delta\delta})]\delta\frac{\partial \rho_r}{\partial x_j}|_{T,V,x_k}$$$$ \frac{\partial p}{\partial x_j}|_{\tau,\delta,x_k} = R\left[\rho T (\delta \frac{\partial}{\partial x_j}[\alpha_\delta]_{\tau, \delta, x_k})+(1+\delta \alpha_\delta)\left(\frac{\rho}{\tau} \frac{\partial T_r}{\partial x_j}+\delta T \frac{\partial \rho_r}{\partial x_j}\right)\right]$$First term cancels with a term in the third one, yielding
$$ \frac{dp}{d\delta}\times\frac{d\delta}{dx_j} = -\delta R T[ (1+\delta \alpha_\delta)+ \delta(\alpha_\delta+\delta \alpha_{\delta\delta})]\frac{\partial \rho_r}{\partial x_j}|_{T,V,x_k}$$$$ \frac{\partial p}{\partial x_j}|_{\tau,\delta,x_k} = \rho R T (\delta \frac{\partial}{\partial x_j}[\alpha_\delta]_{\tau, \delta, x_k})+(1+\delta \alpha_\delta)\delta R T \frac{\partial \rho_r}{\partial x_j}$$First term in first line cancels with term at end of second line, yielding
$$ \frac{dp}{d\delta}\times\frac{d\delta}{dx_j} = -\delta R T[ \delta(\alpha_\delta+\delta \alpha_{\delta\delta})]\frac{\partial \rho_r}{\partial x_j}|_{T,V,x_k}$$$$ \frac{\partial p}{\partial x_j}|_{\tau,\delta,x_k} = \rho R T (\delta \frac{\partial}{\partial x_j}[\alpha_\delta]_{\tau, \delta, x_k})$$Total equation is
$$ \frac{\partial p}{\partial x_j}|_{T,V,x_k} = \rho R T \left(\delta \frac{\partial}{\partial x_j}[\alpha_\delta]_{\tau, \delta, x_k}-\frac{\delta}{\rho_r}(\alpha_\delta+\delta \alpha_{\delta\delta})\frac{\partial \rho_r}{\partial x_j}|_{T,V,x_k}\right)$$
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import CoolProp.CoolProp as CP
%matplotlib inline
In [2]:
p = np.linspace(1000, 12000, 100)
mix = 'REFPROP-MIX:Water[0.7]&Ethanol[0.3]'
rhoL = CP.Props('D','P',p,'Q',0,mix)
rhoV = CP.Props('D','P',p,'Q',1,mix)
In [3]:
plt.plot(rhoL,p,rhoV,p)
Out[3]:
In [3]:
K_i=\frac{y_i}{x_i}