Ideal VLE and Flash Initialization (chemicals.flash_basic)

This module contains the ideal flash solver; two flash initialization routines; a vapor-liquid equilibrium constant correlation; a liquid-water equilibrium constant correlation, and a definition function to show the commonly used calculation frameworks.

For reporting bugs, adding feature requests, or submitting pull requests, please use the GitHub issue tracker.

Ideal Flash Function

chemicals.flash_basic.flash_ideal(zs, funcs, Tcs=None, T=None, P=None, VF=None)

PVT flash model using ideal, composition-independent equation. Solves the various cases of composition-independent models.

Capable of solving with two of T, P, and VF for the other one; that results in three solve modes, but for VF=1 and VF=0, there are additional solvers; for a total of seven solvers implemented.

The function takes a list of callables that take T in Kelvin as an argument, and return vapor pressure. The callables can include the effect of non-ideal pure component fugacity coefficients. For the (T, P) and (P, VF) cases, the Poynting correction factor can be easily included as well but not the (T, VF) case as the callable only takes T as an argument. Normally the Poynting correction factor is used with activity coefficient models with composition dependence.

Both flash_wilson and flash_Tb_Tc_Pc are specialized cases of this function and have the same functionality but with the model built right in.

Even when using more complicated models, this is useful for obtaining initial

This model uses flash_inner_loop to solve the Rachford-Rice problem.

Parameters

zslist[float]

Mole fractions of the phase being flashed, [-]

funcslist[Callable]

Functions to calculate ideal or real vapor pressures, take temperature in Kelvin and return pressure in Pa, [-]

Tcslist[float], optional

Critical temperatures of all species; uses as upper bounds and only for the case that T is not specified; if they are needed and not given, it is assumed a method solve_prop exists in each of funcs which will accept P in Pa and return temperature in K, [K]

Tfloat, optional

Temperature, [K]

Pfloat, optional

Pressure, [Pa]

VFfloat, optional

Molar vapor fraction, [-]

Returns

Tfloat

Temperature, [K]

Pfloat

Pressure, [Pa]

VFfloat

Molar vapor fraction, [-]

xslist[float]

Mole fractions of liquid phase, [-]

yslist[float]

Mole fractions of vapor phase, [-]

Notes

For the cases where VF is 1 or 0 and T is known, an explicit solution is used. For the same cases where P and VF are known, there is no explicit solution available.

There is an internal Tmax parameter, set to 50000 K; which, in the event of convergence of the Secant method, is used as a bounded for a bounded solver. It is used in the PVF solvers.

Examples

Basic case with four compounds, usingthe Antoine equation as a model and solving for vapor pressure:

>>> from chemicals import Antoine, Ambrose_Walton
>>> Tcs = [369.83, 425.12, 469.7, 507.6]
>>> Antoine_As = [8.92828, 8.93266, 8.97786, 9.00139]
>>> Antoine_Bs = [803.997, 935.773, 1064.84, 1170.88]
>>> Antoine_Cs = [-26.11, -34.361, -41.136, -48.833]
>>> Psat_funcs = []
>>> for i in range(4):
...     def Psat_func(T, A=Antoine_As[i], B=Antoine_Bs[i], C=Antoine_Cs[i]):
...         return Antoine(T, A, B, C)
...     Psat_funcs.append(Psat_func)
>>> zs = [.4, .3, .2, .1]
>>> T, P, VF, xs, ys = flash_ideal(T=330.55, P=1e6, zs=zs, funcs=Psat_funcs, Tcs=Tcs)
>>> round(VF, 10)
1.00817e-05

Similar case, using the Ambrose-Walton corresponding states method to estimate vapor pressures:

>>> Tcs = [369.83, 425.12, 469.7, 507.6]
>>> Pcs = [4248000.0, 3796000.0, 3370000.0, 3025000.0]
>>> omegas = [0.152, 0.193, 0.251, 0.2975]
>>> Psat_funcs = []
>>> for i in range(4):
...     def Psat_func(T, Tc=Tcs[i], Pc=Pcs[i], omega=omegas[i]):
...         return Ambrose_Walton(T, Tc, Pc, omega)
...     Psat_funcs.append(Psat_func)
>>> _, P, VF, xs, ys = flash_ideal(T=329.151, VF=0, zs=zs, funcs=Psat_funcs, Tcs=Tcs)
>>> round(P, 3)
1000013.343

Case with fugacities in the liquid phase, vapor phase, activity coefficients in the liquid phase, and Poynting correction factors.

>>> Tcs = [647.14, 514.0]
>>> Antoine_As = [10.1156, 10.3368]
>>> Antoine_Bs = [1687.54, 1648.22]
>>> Antoine_Cs = [-42.98, -42.232]
>>> gammas = [1.1, .75]
>>> fugacities_gas = [.995, 0.98]
>>> fugacities_liq = [.9999, .9998]
>>> Poyntings = [1.000001, .999999]
>>> zs = [.5, .5]
>>> funcs = []
>>> for i in range(2):
...     def K_over_P(T, A=Antoine_As[i], B=Antoine_Bs[i], C=Antoine_Cs[i], fl=fugacities_liq[i],
...                  fg=fugacities_gas[i], gamma=gammas[i], poy=Poyntings[i]):
...         return Antoine(T, A, B, C)*gamma*poy*fl/fg
...     funcs.append(K_over_P)
>>> _, _, VF, xs, ys = flash_ideal(zs, funcs, Tcs=Tcs, P=1e5, T=364.0)
>>> VF, xs, ys
(0.5108639717, [0.55734934039, 0.44265065960], [0.44508982795, 0.554910172040])

Note that while this works for PT composition independent flashes - an outer iterating loop is needed for composition dependence!

Flash Initialization

chemicals.flash_basic.flash_wilson(zs, Tcs, Pcs, omegas, T=None, P=None, VF=None)

PVT flash model using Wilson’s equation - useful for obtaining initial guesses for more rigorous models, or it can be used as its own model. Capable of solving with two of T, P, and VF for the other one; that results in three solve modes, but for VF=1 and VF=0, there are additional solvers; for a total of seven solvers implemented.

This model uses flash_inner_loop to solve the Rachford-Rice problem.

$$K_i = \frac{P_c}{P} \exp\left(5.37(1+\omega)\left[1 - \frac{T_c}{T} \right]\right)$$

Parameters

zslist[float]

Mole fractions of the phase being flashed, [-]

Tcslist[float]

Critical temperatures of all species, [K]

Pcslist[float]

Critical pressures of all species, [Pa]

omegaslist[float]

Acentric factors of all species, [-]

Tfloat, optional

Temperature, [K]

Pfloat, optional

Pressure, [Pa]

VFfloat, optional

Molar vapor fraction, [-]

Returns

Tfloat

Temperature, [K]

Pfloat

Pressure, [Pa]

VFfloat

Molar vapor fraction, [-]

xslist[float]

Mole fractions of liquid phase, [-]

yslist[float]

Mole fractions of vapor phase, [-]

Notes

For the cases where VF is 1 or 0 and T is known, an explicit solution is used. For the same cases where P and VF are known, there is no explicit solution available.

There is an internal Tmax parameter, set to 50000 K; which, in the event of convergence of the Secant method, is used as a bounded for a bounded solver. It is used in the PVF solvers. This typically allows pressures up to 2 GPa to be converged to. However, for narrow-boiling mixtures, the PVF failure may occur at much lower pressures.

Examples

>>> Tcs = [305.322, 540.13]
>>> Pcs = [4872200.0, 2736000.0]
>>> omegas = [0.099, 0.349]
>>> zs = [0.4, 0.6]
>>> flash_wilson(zs=zs, Tcs=Tcs, Pcs=Pcs, omegas=omegas, T=300, P=1e5)
(300, 100000.0, 0.422194532936, [0.02093881508003, 0.979061184919], [0.918774185622, 0.0812258143])

chemicals.flash_basic.flash_Tb_Tc_Pc(zs, Tbs, Tcs, Pcs, T=None, P=None, VF=None)

PVT flash model using a model published in [1], which provides a PT surface using only each compound’s boiling temperature and critical temperature and pressure. This is useful for obtaining initial guesses for more rigorous models, or it can be used as its own model. Capable of solving with two of T, P, and VF for the other one; that results in three solve modes, but for VF=1 and VF=0, there are additional solvers; for a total of seven solvers implemented.

This model uses flash_inner_loop to solve the Rachford-Rice problem.

$$K_i = \frac{P_{c,i}^{\left(\frac{1}{T} - \frac{1}{T_{b,i}} \right) / \left(\frac{1}{T_{c,i}} - \frac{1}{T_{b,i}} \right)}}{P}$$

Parameters

zslist[float]

Mole fractions of the phase being flashed, [-]

Tbslist[float]

Boiling temperatures of all species, [K]

Tcslist[float]

Critical temperatures of all species, [K]

Pcslist[float]

Critical pressures of all species, [Pa]

Tfloat, optional

Temperature, [K]

Pfloat, optional

Pressure, [Pa]

VFfloat, optional

Molar vapor fraction, [-]

Returns

Tfloat

Temperature, [K]

Pfloat

Pressure, [Pa]

VFfloat

Molar vapor fraction, [-]

xslist[float]

Mole fractions of liquid phase, [-]

yslist[float]

Mole fractions of vapor phase, [-]

Notes

For the cases where VF is 1 or 0 and T is known, an explicit solution is used. For the same cases where P and VF are known, there is no explicit solution available.

There is an internal Tmax parameter, set to 50000 K; which, in the event of convergence of the Secant method, is used as a bounded for a bounded solver. It is used in the PVF solvers. This typically allows pressures up to 2 MPa to be converged to. Failures may still occur for other conditions.

This model is based on [1], which aims to estimate dew and bubble points using the same K value formulation as used here. While this implementation uses a numerical solver to provide an exact bubble/dew point estimate, [1] suggests a sequential substitution and flowchart based solver with loose tolerances. That model was also implemented, but found to be slower and less reliable than this implementation.

References

1(1,2,3)

Kandula, Vamshi Krishna, John C. Telotte, and F. Carl Knopf. “It`s Not as Easy as It Looks: Revisiting Peng—Robinson Equation of State Convergence Issues for Dew Point, Bubble Point and Flash Calculations.” International Journal of Mechanical Engineering Education 41, no. 3 (July 1, 2013): 188-202. https://doi.org/10.7227/IJMEE.41.3.2.

Examples

>>> Tcs = [305.322, 540.13]
>>> Pcs = [4872200.0, 2736000.0]
>>> Tbs = [184.55, 371.53]
>>> zs = [0.4, 0.6]
>>> flash_Tb_Tc_Pc(zs=zs, Tcs=Tcs, Pcs=Pcs, Tbs=Tbs, T=300, P=1e5)
(300, 100000.0, 0.3807040748145, [0.0311578430365, 0.968842156963], [0.9999999998827, 1.1729141887e-10])

Equilibrium Constants

chemicals.flash_basic.K_value(P=None, Psat=None, phi_l=None, phi_g=None, gamma=None, Poynting=1.0)

Calculates the equilibrium K-value assuming Raoult’s law, or an equation of state model, or an activity coefficient model, or a combined equation of state-activity model.

The calculation procedure will use the most advanced approach with the provided inputs:

• If P, Psat, phi_l, phi_g, and gamma are provided, use the combined approach.

• If P, Psat, and gamma are provided, use the modified Raoult’s law.

• If phi_l and phi_g are provided, use the EOS only method.

• If P and Psat are provided, use Raoult’s law.

Definitions:

$$K_i=\frac{y_i}{x_i}$$

Raoult’s law:

$$K_i = \frac{P_{i}^{sat}}{P}$$

Activity coefficient, no EOS (modified Raoult’s law):

$$K_i = \frac{\gamma_i P_{i}^{sat}}{P}$$

Equation of state only:

$$K_i = \frac{\phi_i^l}{\phi_i^v} = \frac{f_i^l y_i}{f_i^v x_i}$$

Combined approach (liquid reference fugacity coefficient is normally calculated the saturation pressure for it as a pure species; vapor fugacity coefficient calculated normally):

$$K_i = \frac{\gamma_i P_i^{sat} \phi_i^{l,ref}}{\phi_i^v P}$$

Combined approach, with Poynting Correction Factor (liquid molar volume in the integral is for i as a pure species only):

$$K_i = \frac{\gamma_i P_i^{sat} \phi_i^{l, ref} \exp\left[\frac{ \int_{P_i^{sat}}^P V_i^l dP}{RT}\right]}{\phi_i^v P}$$

Parameters

Pfloat

System pressure, optional

Psatfloat

Vapor pressure of species i, [Pa]

phi_lfloat

Fugacity coefficient of species i in the liquid phase, either at the system conditions (EOS-only case) or at the saturation pressure of species i as a pure species (reference condition for the combined approach), optional [-]

phi_gfloat

Fugacity coefficient of species i in the vapor phase at the system conditions, optional [-]

gammafloat

Activity coefficient of species i in the liquid phase, optional [-]

Poyntingfloat

Poynting correction factor, optional [-]

Returns

Kfloat

Equilibrium K value of component i, calculated with an approach depending on the provided inputs [-]

Notes

The Poynting correction factor is normally simplified as follows, due to a liquid’s low pressure dependency:

$$K_i = \frac{\gamma_i P_i^{sat} \phi_i^{l, ref} \exp\left[\frac{V_l (P-P_i^{sat})}{RT}\right]}{\phi_i^v P}$$

References

1

Gmehling, Jurgen, Barbel Kolbe, Michael Kleiber, and Jurgen Rarey. Chemical Thermodynamics for Process Simulation. 1st edition. Weinheim: Wiley-VCH, 2012.

2

Skogestad, Sigurd. Chemical and Energy Process Engineering. 1st edition. Boca Raton, FL: CRC Press, 2008.

Examples

Raoult’s law:

>>> K_value(101325, 3000.)
0.029607698001480384

Modified Raoult’s law:

>>> K_value(P=101325, Psat=3000, gamma=0.9)
0.026646928201332347

EOS-only approach:

>>> K_value(phi_l=1.6356, phi_g=0.88427)
1.8496613025433408

Gamma-phi combined approach:

>>> K_value(P=1E6, Psat=1938800, phi_l=1.4356, phi_g=0.88427, gamma=0.92)
2.8958055544121137

Gamma-phi combined approach with a Poynting factor:

>>> K_value(P=1E6, Psat=1938800, phi_l=1.4356, phi_g=0.88427, gamma=0.92,
... Poynting=0.999)
2.8929097488577016

chemicals.flash_basic.Wilson_K_value(T, P, Tc, Pc, omega)

Calculates the equilibrium K-value for a component using Wilson’s heuristic mode. This is very useful for initialization of stability tests and flashes.

$$K_i = \frac{P_c}{P} \exp\left(5.37(1+\omega)\left[1 - \frac{T_c}{T} \right]\right)$$

Parameters

Tfloat

System temperature, [K]

Pfloat

System pressure, [Pa]

Tcfloat

Critical temperature of fluid [K]

Pcfloat

Critical pressure of fluid [Pa]

omegafloat

Acentric factor for fluid, [-]

Returns

Kfloat

Equilibrium K value of component, calculated via the Wilson heuristic [-]

Notes

There has been little literature exploration of other formlulas for the same purpose. This model may be useful even for activity coefficient models.

Note the K-values are independent of composition; the correlation is applicable up to 3.5 MPa.

A description for how this function was generated can be found in [2].

References

1

Wilson, Grant M. “A Modified Redlich-Kwong Equation of State, Application to General Physical Data Calculations.” In 65th National AIChE Meeting, Cleveland, OH, 1969.

2

Peng, Ding-Yu, and Donald B. Robinson. “Two and Three Phase Equilibrium Calculations for Systems Containing Water.” The Canadian Journal of Chemical Engineering, December 1, 1976. https://doi.org/10.1002/cjce.5450540620.

Examples

Ethane at 270 K and 76 bar:

>>> Wilson_K_value(270.0, 7600000.0, 305.4, 4880000.0, 0.098)
0.2963932297479371

The “vapor pressure” predicted by this equation can be calculated by multiplying by pressure:

>>> Wilson_K_value(270.0, 7600000.0, 305.4, 4880000.0, 0.098)*7600000.0
2252588.546084322

chemicals.flash_basic.PR_water_K_value(T, P, Tc, Pc)

Calculates the equilibrium K-value for a component against water according to the Peng and Robinson (1976) heuristic.

$$K_i = 10^6 \frac{P_{ri}}{T_{ri}}$$

Parameters

Tfloat

System temperature, [K]

Pfloat

System pressure, [Pa]

Tcfloat

Critical temperature of chemical [K]

Pcfloat

Critical pressure of chemical [Pa]

Returns

Kfloat

Equilibrium K value of component with water as the other phase ( not as the reference), calculated via this heuristic [-]

Notes

Note the K-values are independent of composition.

References

1

Peng, Ding-Yu, and Donald B. Robinson. “Two and Three Phase Equilibrium Calculations for Systems Containing Water.” The Canadian Journal of Chemical Engineering, December 1, 1976. https://doi.org/10.1002/cjce.5450540620.

Examples

Octane at 300 K and 1 bar:

>>> PR_water_K_value(300, 1e5, 568.7, 2490000.0)
76131.19143239626