Skip to content
Snippets Groups Projects
Commit 750c9adf authored by Ian Bell's avatar Ian Bell
Browse files

Switch to newer API

parent a6787c13
Branches
No related tags found
No related merge requests found
%% Cell type:markdown id:8218498b tags: %% Cell type:markdown id:8218498b tags:
# Phase equilibria # Phase equilibria
Two basic approaches are implemented in teqp: Two basic approaches are implemented in teqp:
* Iterative calculations given guess values * Iterative calculations given guess values
* Tracing along iso-curves (constant temperature, etc.) powered by the isochoric thermodynamics formalism * Tracing along iso-curves (constant temperature, etc.) powered by the isochoric thermodynamics formalism
%% Cell type:code id:c0b4e863 tags: %% Cell type:code id:c0b4e863 tags:
``` python ``` python
import teqp import teqp
import numpy as np import numpy as np
import pandas import pandas
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
teqp.__version__ teqp.__version__
``` ```
%% Cell type:markdown id:e57e532b tags: %% Cell type:markdown id:e57e532b tags:
## Iterative Phase Equilibria ## Iterative Phase Equilibria
%% Cell type:markdown id:1baf38e1 tags: %% Cell type:markdown id:1baf38e1 tags:
### Pure fluid ### Pure fluid
For a pure fluid, phase equilibrium between two phases is defined by equating the pressures and Gibbs energies in the two phases. This represents a 2D non-linear rootfinding problem. Newton's method can be used for the rootfinding, and in teqp, automatic differentiation is used to obtain the necessary Jacobian matrix so the implementation is quite efficient. For a pure fluid, phase equilibrium between two phases is defined by equating the pressures and Gibbs energies in the two phases. This represents a 2D non-linear rootfinding problem. Newton's method can be used for the rootfinding, and in teqp, automatic differentiation is used to obtain the necessary Jacobian matrix so the implementation is quite efficient.
The method requires guess values, which are the densities of the liquid and vapor densities. In some cases, ancillary or superancillary equations have been developed which provide curves of guess densities as a function of temperature. The method requires guess values, which are the densities of the liquid and vapor densities. In some cases, ancillary or superancillary equations have been developed which provide curves of guess densities as a function of temperature.
For a pure fluid, you can use the ``pure_VLE_T`` method to carry out the iteration. For a pure fluid, you can use the ``pure_VLE_T`` method to carry out the iteration.
%% Cell type:raw id:f0ca0b22 tags: %% Cell type:raw id:f0ca0b22 tags:
The Python method is here: :py:meth:`pure_VLE_T <teqp.teqp.pure_VLE_T>` The Python method is here: :py:meth:`pure_VLE_T <teqp.teqp.pure_VLE_T>`
%% Cell type:code id:2674227c tags: %% Cell type:code id:2674227c tags:
``` python ``` python
# Instantiate the model # Instantiate the model
model = teqp.canonical_PR([300], [4e6], [0.1]) model = teqp.canonical_PR([300], [4e6], [0.1])
T = 250 # [K], Temperature to be used T = 250 # [K], Temperature to be used
# Here we use the superancillary to get guess values (actually these are more # Here we use the superancillary to get guess values (actually these are more
# accurate than the results we will obtain from iteration!) # accurate than the results we will obtain from iteration!)
rhoL0, rhoV0 = model.superanc_rhoLV(T) rhoL0, rhoV0 = model.superanc_rhoLV(T)
display('guess:', [rhoL0, rhoV0]) display('guess:', [rhoL0, rhoV0])
# Carry out the iteration, return the liquid and vapor densities # Carry out the iteration, return the liquid and vapor densities
# The guess values are perturbed to make sure the iteration is actually # The guess values are perturbed to make sure the iteration is actually
# changing the values # changing the values
teqp.pure_VLE_T(model, T, rhoL0*0.98, rhoV0*1.02, 10) model.pure_VLE_T(T, rhoL0*0.98, rhoV0*1.02, 10)
``` ```
%% Cell type:markdown id:f8805ae1 tags: %% Cell type:markdown id:f8805ae1 tags:
### Binary Mixture ### Binary Mixture
%% Cell type:markdown id:76bccf19 tags: %% Cell type:markdown id:76bccf19 tags:
For a binary mixture, the approach is roughly similar to that of a pure fluid. The pressure is equated between phases, and the chemical potentials of each component in each phase are forced to be the same. For a binary mixture, the approach is roughly similar to that of a pure fluid. The pressure is equated between phases, and the chemical potentials of each component in each phase are forced to be the same.
Again, the user is required to provide guess values, in this case molar concentrations in each phase, and a Newton method is implemented to solve for the phase equilibrium. The analytical Jacobian is obtained from automatic differentiation. Again, the user is required to provide guess values, in this case molar concentrations in each phase, and a Newton method is implemented to solve for the phase equilibrium. The analytical Jacobian is obtained from automatic differentiation.
The ``mix_VLE_Tx`` function is the binary mixture analog to ``pure_VLE_T`` for pure fluids. The ``mix_VLE_Tx`` function is the binary mixture analog to ``pure_VLE_T`` for pure fluids.
%% Cell type:raw id:eef189fd tags: %% Cell type:raw id:eef189fd tags:
The Python method is here: :py:meth:`mix_VLE_Tx <teqp.teqp.mix_VLE_Tx>` The Python method is here: :py:meth:`mix_VLE_Tx <teqp.teqp.mix_VLE_Tx>`
%% Cell type:code id:b12bd318 tags: %% Cell type:code id:b12bd318 tags:
``` python ``` python
zA = np.array([0.01, 0.99]) zA = np.array([0.01, 0.99])
model = teqp.canonical_PR([300,310], [4e6,4.5e6], [0.1, 0.2]) model = teqp.canonical_PR([300,310], [4e6,4.5e6], [0.1, 0.2])
model1 = teqp.canonical_PR([300], [4e6], [0.1]) model1 = teqp.canonical_PR([300], [4e6], [0.1])
T = 273.0 # [K] T = 273.0 # [K]
# start off at pure of the first component # start off at pure of the first component
rhoL0, rhoV0 = model1.superanc_rhoLV(T) rhoL0, rhoV0 = model1.superanc_rhoLV(T)
# then we shift to the given composition in the first phase # then we shift to the given composition in the first phase
# to get guess values # to get guess values
rhovecA0 = rhoL0*zA rhovecA0 = rhoL0*zA
rhovecB0 = rhoV0*zA rhovecB0 = rhoV0*zA
# carry out the iteration # carry out the iteration
code, rhovecA, rhovecB = teqp.mix_VLE_Tx(model, T, rhovecA0, rhovecB0, zA, code, rhovecA, rhovecB = model.mix_VLE_Tx(T, rhovecA0, rhovecB0, zA,
1e-10, 1e-10, 1e-10, 1e-10, # stopping conditions 1e-10, 1e-10, 1e-10, 1e-10, # stopping conditions
10 # maximum number of iterations 10 # maximum number of iterations
) )
code, rhovecA, rhovecB code, rhovecA, rhovecB
``` ```
%% Cell type:markdown id:d72ef08e tags: %% Cell type:markdown id:d72ef08e tags:
You can (and should) check the value of the return code to make sure the iteration succeeded. Do not rely on the numerical value of the enumerated return codes! You can (and should) check the value of the return code to make sure the iteration succeeded. Do not rely on the numerical value of the enumerated return codes!
%% Cell type:markdown id:f4e3f914 tags: %% Cell type:markdown id:f4e3f914 tags:
# Tracing (isobars and isotherms) # Tracing (isobars and isotherms)
When it comes to mixture thermodynamics, as soon as you add another component to a pure component to form a binary mixture, the complexity of the thermodynamics entirely changes. For that reason, mixture iterative calculations for mixtures are orders of magnitude more difficult to carry out. Asymmetric mixtures can do all sorts of interesting things that are entirely unlike those of pure fluids, and the algorithms are therefore much, much more complicated. Formulating phase equilibrium problems is not much more complicated than for pure fluids, but the most challenging aspect is to obtain good guess values from which to start an iterative routine, and the difficulty of this problem increases with the complexity of the mixture thermodynamics. When it comes to mixture thermodynamics, as soon as you add another component to a pure component to form a binary mixture, the complexity of the thermodynamics entirely changes. For that reason, mixture iterative calculations for mixtures are orders of magnitude more difficult to carry out. Asymmetric mixtures can do all sorts of interesting things that are entirely unlike those of pure fluids, and the algorithms are therefore much, much more complicated. Formulating phase equilibrium problems is not much more complicated than for pure fluids, but the most challenging aspect is to obtain good guess values from which to start an iterative routine, and the difficulty of this problem increases with the complexity of the mixture thermodynamics.
Ulrich Deiters and Ian Bell have developed a number of algorithms for tracing phase equilibrium solutions as the solution of ordinary differential equations rather than carrying out iterative routines for a given state point. The advantage of the tracing calculations is that they can often be initiated at a state point that is entirely known, for instance the pure fluid endpoint for a subcritical isotherm or isobar. Ulrich Deiters and Ian Bell have developed a number of algorithms for tracing phase equilibrium solutions as the solution of ordinary differential equations rather than carrying out iterative routines for a given state point. The advantage of the tracing calculations is that they can often be initiated at a state point that is entirely known, for instance the pure fluid endpoint for a subcritical isotherm or isobar.
%% Cell type:raw id:e0097771 tags: %% Cell type:raw id:e0097771 tags:
The Python method is here: :py:meth:`trace_VLE_isotherm_binary <teqp.teqp.trace_VLE_isotherm_binary>` The Python method is here: :py:meth:`trace_VLE_isotherm_binary <teqp.teqp.trace_VLE_isotherm_binary>`
%% Cell type:markdown id:63902dba tags: %% Cell type:markdown id:63902dba tags:
The C++ implementation returns a string in JSON format, which can be conveniently operated upon, for instance after converting the returned data structure to a ``pandas.DataFrame``. A simple example of plotting a subcritical isotherm for a "boring" mixture is presented here: The C++ implementation returns a string in JSON format, which can be conveniently operated upon, for instance after converting the returned data structure to a ``pandas.DataFrame``. A simple example of plotting a subcritical isotherm for a "boring" mixture is presented here:
%% Cell type:code id:49dcba2b tags: %% Cell type:code id:49dcba2b tags:
``` python ``` python
model = teqp.canonical_PR([300,310], [4e6,4.5e6], [0.1, 0.2]) model = teqp.canonical_PR([300,310], [4e6,4.5e6], [0.1, 0.2])
model1 = teqp.canonical_PR([300], [4e6], [0.1]) model1 = teqp.canonical_PR([300], [4e6], [0.1])
T = 273.0 # [K] T = 273.0 # [K]
rhoL0, rhoV0 = model1.superanc_rhoLV(T) # start off at pure of the first component rhoL0, rhoV0 = model1.superanc_rhoLV(T) # start off at pure of the first component
j = teqp.trace_VLE_isotherm_binary(model, T, np.array([rhoL0, 0]), np.array([rhoV0, 0])) j = model.trace_VLE_isotherm_binary(T, np.array([rhoL0, 0]), np.array([rhoV0, 0]))
display(str(j)[0:100]+'...') # The first few bits of the data display(str(j)[0:100]+'...') # The first few bits of the data
df = pandas.DataFrame(j) # Now as a data frame df = pandas.DataFrame(j) # Now as a data frame
df.head(3) df.head(3)
``` ```
%% Cell type:code id:9aecca78 tags: %% Cell type:code id:9aecca78 tags:
``` python ``` python
plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6) plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6)
plt.plot(df['xV_0 / mole frac.'], df['pL / Pa']/1e6) plt.plot(df['xV_0 / mole frac.'], df['pL / Pa']/1e6)
plt.gca().set(xlabel='$x_1,y_1$ / mole frac.', ylabel='p / MPa') plt.gca().set(xlabel='$x_1,y_1$ / mole frac.', ylabel='p / MPa')
plt.show() plt.show()
``` ```
%% Cell type:markdown id:a9c9fe55 tags: %% Cell type:markdown id:a9c9fe55 tags:
Isn't that exciting! Isn't that exciting!
You can also provide an optional set of flags to the function to control other behaviors of the function, and switch between simple Euler and adaptive RK45 integration (the default) You can also provide an optional set of flags to the function to control other behaviors of the function, and switch between simple Euler and adaptive RK45 integration (the default)
%% Cell type:raw id:264c5123 tags: %% Cell type:raw id:264c5123 tags:
The options class is here: :py:meth:`TVLEOptions <teqp.teqp.TVLEOptions>` The options class is here: :py:meth:`TVLEOptions <teqp.teqp.TVLEOptions>`
%% Cell type:markdown id:d4193110 tags: %% Cell type:markdown id:d4193110 tags:
Supercritical isotherms work approximately in the same manner Supercritical isotherms work approximately in the same manner
%% Cell type:code id:c5b925ad tags: %% Cell type:code id:c5b925ad tags:
``` python ``` python
Tc_K = [190.564, 154.581] Tc_K = [190.564, 154.581]
pc_Pa = [4599200, 5042800] pc_Pa = [4599200, 5042800]
acentric = [0.011, 0.022] acentric = [0.011, 0.022]
model = teqp.canonical_PR(Tc_K, pc_Pa, acentric) model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)
model1 = teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]]) model1 = teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]])
T = 170.0 # [K] # Note: above Tc of the second component T = 170.0 # [K] # Note: above Tc of the second component
rhoL0, rhoV0 = model1.superanc_rhoLV(T) # start off at pure of the first component rhoL0, rhoV0 = model1.superanc_rhoLV(T) # start off at pure of the first component
j = teqp.trace_VLE_isotherm_binary(model, T, np.array([rhoL0, 0]), np.array([rhoV0, 0])) j = model.trace_VLE_isotherm_binary(T, np.array([rhoL0, 0]), np.array([rhoV0, 0]))
df = pandas.DataFrame(j) # Now as a data frame df = pandas.DataFrame(j) # Now as a data frame
plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6) plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6)
plt.plot(df['xV_0 / mole frac.'], df['pL / Pa']/1e6) plt.plot(df['xV_0 / mole frac.'], df['pL / Pa']/1e6)
plt.gca().set(xlabel='$x_1,y_1$ / mole frac.', ylabel='p / MPa') plt.gca().set(xlabel='$x_1,y_1$ / mole frac.', ylabel='p / MPa')
plt.show() plt.show()
``` ```
%% Cell type:markdown id:936c5f55 tags: %% Cell type:markdown id:936c5f55 tags:
As of version 0.10.0, isobar tracing has been added to ``teqp``. It operates in fundamentally the same fashion as the isotherm tracing and the same recommendations about starting at a pure fluid apply As of version 0.10.0, isobar tracing has been added to ``teqp``. It operates in fundamentally the same fashion as the isotherm tracing and the same recommendations about starting at a pure fluid apply
%% Cell type:raw id:81eacaf0 tags: %% Cell type:raw id:81eacaf0 tags:
The tracer function class is here: :py:meth:`trace_VLE_isobar_binary <teqp.teqp.trace_VLE_isobar_binary>` The tracer function class is here: :py:meth:`trace_VLE_isobar_binary <teqp.teqp.trace_VLE_isobar_binary>`
%% Cell type:code id:109bc65b tags: %% Cell type:code id:109bc65b tags:
``` python ``` python
Tc_K = [190.564, 154.581] Tc_K = [190.564, 154.581]
pc_Pa = [4599200, 5042800] pc_Pa = [4599200, 5042800]
acentric = [0.011, 0.022] acentric = [0.011, 0.022]
model = teqp.canonical_PR(Tc_K, pc_Pa, acentric) model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)
model1 = teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]]) model1 = teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]])
T = 170.0 # [K] # Note: above Tc of the second component T = 170.0 # [K] # Note: above Tc of the second component
rhoL0, rhoV0 = model1.superanc_rhoLV(T) # start off at pure of the first component rhoL0, rhoV0 = model1.superanc_rhoLV(T) # start off at pure of the first component
p0 = rhoL0*model1.get_R(np.array([1.0]))*T*(1+model1.get_Ar01(T, rhoL0, np.array([1.0]))) p0 = rhoL0*model1.get_R(np.array([1.0]))*T*(1+model1.get_Ar01(T, rhoL0, np.array([1.0])))
j = teqp.trace_VLE_isobar_binary(model, p0, T, np.array([rhoL0, 0]), np.array([rhoV0, 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 df = pandas.DataFrame(j) # Now as a data frame
plt.plot(df['xL_0 / mole frac.'], df['T / K']) plt.plot(df['xL_0 / mole frac.'], df['T / K'])
plt.plot(df['xV_0 / mole frac.'], df['T / K']) plt.plot(df['xV_0 / mole frac.'], df['T / K'])
plt.gca().set(xlabel='$x_1,y_1$ / mole frac.', ylabel='T / K') plt.gca().set(xlabel='$x_1,y_1$ / mole frac.', ylabel='T / K')
plt.show() plt.show()
``` ```
......
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment