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

Mixture sublevel

parent e97aaad9
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
doc/source/algorithms/critical_curves.ipynb
+ 1
− 1
View file @ 1638c830
%% Cell type:code id:d88bffed tags: %% Cell type:code id:d88bffed tags:
``` python ``` python
import teqp import teqp
teqp.__version__ teqp.__version__
``` ```
%% Cell type:markdown id:8218498b tags: %% Cell type:markdown id:8218498b tags:
# Critical curves & points # Critical curves & points
## Pure Fluids ## Pure Fluids
Solving for the critical point involves finding the temperature and density that make Solving for the critical point involves finding the temperature and density that make
$$ $$
\left(\frac{\partial p}{\partial \rho}\right)_T = \left(\frac{\partial^2 p}{\partial \rho^2}\right)_T = 0 \left(\frac{\partial p}{\partial \rho}\right)_T = \left(\frac{\partial^2 p}{\partial \rho^2}\right)_T = 0
$$ $$
by 2D non-linear rootfinding. Newton steps are taken, and the analytic Jacobian is used (thanks to the ability to do derivatives with automatic differentiation). This is all handily wrapped up in the by 2D non-linear rootfinding. Newton steps are taken, and the analytic Jacobian is used (thanks to the ability to do derivatives with automatic differentiation). This is all handily wrapped up in the
``solve_pure_critical`` method which requires the user to provide guess values for temperature and density ``solve_pure_critical`` method which requires the user to provide guess values for temperature and density
%% Cell type:code id:46657a96 tags: %% Cell type:code id:46657a96 tags:
``` python ``` python
# Values taken from http://dx.doi.org/10.6028/jres.121.011 # Values taken from http://dx.doi.org/10.6028/jres.121.011
modelPR = teqp.canonical_PR([190.564], [4599200], [0.011]) modelPR = teqp.canonical_PR([190.564], [4599200], [0.011])
# Solve for the critical point from a point close to the critical point # Solve for the critical point from a point close to the critical point
T0 = 192.0 T0 = 192.0
# Critical compressibility factor of P-R is 0.307401308698.. (see https://doi.org/10.1021/acs.iecr.1c00847) # Critical compressibility factor of P-R is 0.307401308698.. (see https://doi.org/10.1021/acs.iecr.1c00847)
rhoc = (4599200/(8.31446261815324*190.564))/0.3074 rhoc = (4599200/(8.31446261815324*190.564))/0.3074
rho0 = rhoc*1.2345 # Perturb to make sure we are doing something in the solver rho0 = rhoc*1.2345 # Perturb to make sure we are doing something in the solver
teqp.solve_pure_critical(modelPR, T0, rho0) teqp.solve_pure_critical(modelPR, T0, rho0)
``` ```
%% Cell type:markdown id:e15eeded tags: %% Cell type:markdown id:e15eeded tags:
# Mixtures ## Mixtures
A pure fluid has a single vapor-liquid critical point, but mixtures are different: A pure fluid has a single vapor-liquid critical point, but mixtures are different:
* They may have multiple (or zero!) critical points for a given mixture composition * They may have multiple (or zero!) critical points for a given mixture composition
* The critical curves may not emanate from the pure fluid endpoints * The critical curves may not emanate from the pure fluid endpoints
When it comes to critical points, intuition from pure fluids is not helpful, or sometimes even counter-productive. When it comes to critical points, intuition from pure fluids is not helpful, or sometimes even counter-productive.
teqp has methods for working with the critical loci of binary mixtures (only binary mixtures, for now) and especially, methods for tracing the critical curves emanating from the pure fluid endpoints. teqp has methods for working with the critical loci of binary mixtures (only binary mixtures, for now) and especially, methods for tracing the critical curves emanating from the pure fluid endpoints.
The tracing method in teqp is based explicitly on the isochoric thermodynamics formalism introduced by Ulrich Deiters and Sergio Quinones-Cisneros. It uses the Helmholtz energy density as the fundamental potential and all other properties are derived from it. For critical curves it is based upon the integration of sets of ordinary differential equations; the differential equations are in the form of derivatives of the molar concentrations of each component in the mixture with respect to an integration variable. The set of ODE is then integrated. The tracing method in teqp is based explicitly on the isochoric thermodynamics formalism introduced by Ulrich Deiters and Sergio Quinones-Cisneros. It uses the Helmholtz energy density as the fundamental potential and all other properties are derived from it. For critical curves it is based upon the integration of sets of ordinary differential equations; the differential equations are in the form of derivatives of the molar concentrations of each component in the mixture with respect to an integration variable. The set of ODE is then integrated.
Here is an example of the construction of the critical curves emanating from the pure fluid endpoints for the mixture nitrogen + ethane. Here is an example of the construction of the critical curves emanating from the pure fluid endpoints for the mixture nitrogen + ethane.
%% Cell type:code id:81619ae2 tags: %% Cell type:code id:81619ae2 tags:
``` python ``` python
import timeit import timeit
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import pandas import pandas
import teqp import teqp
def get_critical_curve(ipure): def get_critical_curve(ipure):
""" Return curve as pandas DataFrame """ """ Return curve as pandas DataFrame """
names = ['Nitrogen', 'Ethane'] names = ['Nitrogen', 'Ethane']
model = teqp.build_multifluid_model(names, teqp.get_datapath()) model = teqp.build_multifluid_model(names, teqp.get_datapath())
T0 = model.get_Tcvec()[ipure] T0 = model.get_Tcvec()[ipure]
rho0 = np.array([1.0/model.get_vcvec()[ipure]]*2) rho0 = np.array([1.0/model.get_vcvec()[ipure]]*2)
rho0[1-ipure] = 0 rho0[1-ipure] = 0
o = teqp.TCABOptions() o = teqp.TCABOptions()
o.init_dt = 1.0 # step in the parameter o.init_dt = 1.0 # step in the parameter
o.rel_err = 1e-8 o.rel_err = 1e-8
o.abs_err = 1e-5 o.abs_err = 1e-5
o.integration_order = 5 o.integration_order = 5
o.calc_stability = True o.calc_stability = True
o.polish = True o.polish = True
curveJSON = teqp.trace_critical_arclength_binary(model, T0, rho0, '', o) curveJSON = teqp.trace_critical_arclength_binary(model, T0, rho0, '', o)
df = pandas.DataFrame(curveJSON) df = pandas.DataFrame(curveJSON)
rhotot = df['rho0 / mol/m^3']+df['rho1 / mol/m^3'] rhotot = df['rho0 / mol/m^3']+df['rho1 / mol/m^3']
df['z0 / mole frac.'] = df['rho0 / mol/m^3']/rhotot df['z0 / mole frac.'] = df['rho0 / mol/m^3']/rhotot
return df return df
fig, ax = plt.subplots(1,1,figsize=(7, 6)) fig, ax = plt.subplots(1,1,figsize=(7, 6))
tic = timeit.default_timer() tic = timeit.default_timer()
for ipure in [1,0]: for ipure in [1,0]:
df = get_critical_curve(ipure) df = get_critical_curve(ipure)
first_unstable = np.argmax(~df['locally stable']) first_unstable = np.argmax(~df['locally stable'])
df = df.iloc[0:(first_unstable if first_unstable else len(df))] df = df.iloc[0:(first_unstable if first_unstable else len(df))]
line, = plt.plot(df['T / K'], df['p / Pa']/1e6, '-') line, = plt.plot(df['T / K'], df['p / Pa']/1e6, '-')
plt.plot(df['T / K'].iloc[0], df['p / Pa'].iloc[0]/1e6, 'd', plt.plot(df['T / K'].iloc[0], df['p / Pa'].iloc[0]/1e6, 'd',
color=line.get_color()) color=line.get_color())
elap = timeit.default_timer()-tic elap = timeit.default_timer()-tic
plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / MPa', plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / MPa',
xlim=(100, 350), ylim=(1, 1e3)) xlim=(100, 350), ylim=(1, 1e3))
plt.yscale('log') plt.yscale('log')
plt.tight_layout(pad=0.2) plt.tight_layout(pad=0.2)
plt.gcf().text(0,0,f'time: {elap:0.1f} s', ha='left', va='bottom', fontsize=7) plt.gcf().text(0,0,f'time: {elap:0.1f} s', ha='left', va='bottom', fontsize=7)
plt.gcf().text(1,0,f'teqp: {teqp.__version__}', ha='right', va='bottom', fontsize=7); plt.gcf().text(1,0,f'teqp: {teqp.__version__}', ha='right', va='bottom', fontsize=7);
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment