Generalized Phase Equilibrium

New in version 0.22 is a set of generalized routines for building the residual and Jacobian for carrying out phase equilibrium for mixtures of arbitrary number of phases and components.

Theory

There is a need for generalized phase equilibrium routines capable of handling mixtures with an arbitrary number of phases and an arbitrary number of components. Furthermore, it is desired to handle generically a wide range of phase equilibrium specification problems:

  • Bubble-point at specified \(T\) or \(p\) (saturated liquid for two phases in equilibrium)

  • Dew-point at specified \(T\) or \(p\) (saturated vapor for two phases in equilibrium)

  • TP flash

  • PH flash

  • XY flash (the more general case of any two specified thermodynamic variables selected from \(T\),\(p\),\(\rho\),\(h\),\(s\),\(u\))

The independent variables in this formulation are:

  • The temperature \(T\) (same in all phases)

  • The molar concentrations of all components in all phases, one \(\vec\rho\) per phase. Molar concentration \(\rho_i=x_i\rho\) where \(x_i\) is the mole fraction

  • The molar phase fraction for each phase (amount of substance in the phase divided by the total amount of substance)

Thus there are \((N+1)\pi+1\) independent variables if \(N\) is the number of components and \(\pi\) the number of phases

Note that pressure is NOT an independent variable. It is enforced to be equal between the phases as a specification, since the models in teqp do NOT have pressure as one of the independent variables and it is more natural to consider densities as independent variables.

The specification equations that must always be satisfied are:

  • Equality of fugacity of all components in all phases

  • Equal pressures in all phases

  • Material balances

  • Summation of molar phase fractions should be 1.0

This leaves two additional specification options, to be selected from:

  • \(T\)

  • \(p\)

  • \(\beta\) of a phase

  • \(v\) (or equivalently \(\rho\)) of the overall system

  • \(h\), \(s\), \(u\) (WIP)

Note that you must provide guess values for the values for all phases. In some cases the guess values will be trivial, for instance if you are specifying the temperature as a specification equation, you know exactly the right guess value for temperature. In other cases the guesses (especially for molar concentration) are very difficult to come by. As a user, you are responsible to find some way to determine what starting values to use as obtaining them is situation dependent.

Implementation

The overall class is called GeneralizedPhaseEquilibrium() and takes in the model, the bulk composition, the initial set of variables (used to initialize the number of phases and components), and the two additional specification equations.

Methods are available for building the residual vector and the Jacobian matrix. In the call method, no special treatment is done of the entries in the Jacobian that might be required (the use of logarithmic molar concentrations, etc.) or special handling of components that have zero mole fractions. To iterate towards the true phase equilibrium solution, you can use conventional Newton iterations:

\[\mathbf{J}\Delta \mathbf{x} = -\mathbf{r}\]

and for better stability you can take smaller steps with

\[\mathbf{J}\Delta \mathbf{x} = -\omega\mathbf{r}\]

with \(0<\omega<1\) and update \(\mathbf{x}\) with

\[\mathbf{x}_{\rm new} = \mathbf{x}_{\rm old} + \mathbf{\Delta x}\]

After calling the call method, the residual and Jacobian can be obtained from the res.r and res.J attributes, respectively.

[1]:

# Here is an example of generating phase equilibrium data
# from an isobaric trace with the Peng-Robinson EOS and
# then calculating the residuals with the new routine.
#
# In this case no iteration
# is required since the residuals should all be close to zero
# because polishing is enabled by default in the tracing routine
import numpy as np
import pandas

import teqp
from teqp import phaseequil as pe

# Get a good result from the tracing
Tc_K = [190.564, 154.581]
pc_Pa = [4599200, 5042800]
acentric = [0.011, 0.022]
model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)
T = 170.0 # [K] # Note: above Tc of the second component
rhoL0, rhoV0 = model.superanc_rhoLV(T, ifluid=0) # start off at pure of the first component
p0 = rhoL0*model.get_R(np.array([1.0,0]))*T*(1+model.get_Ar01(T, rhoL0, np.array([1.0,0])))
j = model.trace_VLE_isobar_binary(p0, T, np.array([rhoL0, 0]), np.array([rhoV0, 0]))
df = pandas.DataFrame(j) # Now as a data frame

for ir, row in df.iterrows():
    # Only do every fifth, and skip the first because it is infinite dilution
    # which requires special treatment
    if (ir+1)%5 != 0: continue

    # The initial values of the variables
    T = row['T / K']
    rhovecs = [row['rhoL / mol/m^3'], row['rhoV / mol/m^3']]
    betas = [1.0, 0.0]
    unpacked = pe.UnpackedVariables(T, rhovecs, betas)

    zbulk = rhovecs[0]/np.sum(rhovecs[0])
    specs = [
        pe.TSpecification(T),
        pe.BetaSpecification(1.0, 0) # Bubble point calculation, first phase (liquid) with index 0 is the whole mixture
        # or you could consider instead to specify pressure, or ...
    ]
    gpe = pe.GeneralizedPhaseEquilibrium(model, zbulk, unpacked, specs)

    Xinit = unpacked.pack()
    gpe.call(Xinit)
    print(np.max(np.abs(gpe.res.r)))

1.7763568394002505e-15
2.8405338525772095e-08
1.7695128917694092e-08
1.5366822481155396e-08
1.816079020500183e-08
4.330649971961975e-08
1.5832483768463135e-08
1.4901161193847656e-08
1.5366822481155396e-08
1.4901161193847656e-08
1.6298145055770874e-08
1.862645149230957e-09
2.7939677238464355e-08
9.033828973770142e-08
2.7939677238464355e-08
2.8405338525772095e-08
4.7497451305389404e-08
2.8870999813079834e-08
1.7229467630386353e-08
4.656612873077393e-10