Skip to content

Commit

Permalink
Add routines for generalized phase equilibrium
Browse files Browse the repository at this point in the history
Generally working! Need to add additional tests after implementing the properties that require ideal-gas part
  • Loading branch information
ianhbell committed Aug 22, 2024
1 parent f74c747 commit f5fe546
Show file tree
Hide file tree
Showing 5 changed files with 681 additions and 0 deletions.
193 changes: 193 additions & 0 deletions doc/source/algorithms/GeneralizedPhaseEquilibria.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "b0a74d8c",
"metadata": {},
"source": [
"# Generalized Phase Equilibrium\n",
"\n",
"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.\n",
"\n",
"## Theory\n",
"\n",
"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:\n",
"\n",
"* Bubble-point at specified $T$ or $p$ (saturated liquid for two phases in equilibrium)\n",
"* Dew-point at specified $T$ or $p$ (saturated vapor for two phases in equilibrium)\n",
"* TP flash\n",
"* PH flash \n",
"* XY flash (the more general case of any two specified thermodynamic variables selected from $T$,$p$,$\\rho$,$h$,$s$,$u$)\n",
"\n",
"The independent variables in this formulation are:\n",
"\n",
"* The temperature $T$ (same in all phases)\n",
"* 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\n",
"* The molar phase fraction for each phase (amount of substance in the phase divided by the total amount of substance)\n",
"\n",
"Thus there are $(N+1)\\pi+1$ independent variables if $N$ is the number of components and $\\pi$ the number of phases\n",
"\n",
"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.\n",
"\n",
"The specification equations that must always be satisfied are:\n",
"\n",
"* Equality of fugacity of all components in all phases\n",
"* Equal pressures in all phases\n",
"* Material balances\n",
"* Summation of molar phase fractions should be 1.0\n",
"\n",
"This leaves two additional specification options, to be selected from:\n",
"\n",
"* $T$\n",
"* $p$\n",
"* $\\beta$ of a phase\n",
"* $v$ (or equivalently $\\rho$) of the overall system\n",
"* $h$, $s$, $u$ (WIP)\n",
"\n",
"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."
]
},
{
"cell_type": "markdown",
"id": "c7622c7c",
"metadata": {},
"source": [
"## Implementation"
]
},
{
"cell_type": "raw",
"id": "f0ca0b22",
"metadata": {
"raw_mimetype": "text/restructuredtext"
},
"source": [
"The overall class is called :py:meth:`GeneralizedPhaseEquilibrium <teqp.teqp.phaseequil.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."
]
},
{
"cell_type": "markdown",
"id": "f13f6ac3",
"metadata": {},
"source": [
"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:\n",
"$$\n",
"\\mathbf{J}\\Delta \\mathbf{x} = -\\mathbf{r}\n",
"$$\n",
"and for better stability you can take smaller steps with \n",
"$$\n",
"\\mathbf{J}\\Delta \\mathbf{x} = -\\omega\\mathbf{r}\n",
"$$\n",
"with $0<\\omega<1$ and update $\\mathbf{x}$ with \n",
"$$\n",
"\\mathbf{x}_{\\rm new} = \\mathbf{x}_{\\rm old} + \\mathbf{\\Delta x}\n",
"$$\n",
"After calling the call method, the residual and Jacobian can be obtained from the ``res.r`` and ``res.J`` attributes, respectively."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "5da5ee18",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2.3283064365386963e-09\n",
"1.3969838619232178e-08\n",
"1.3504177331924438e-08\n",
"1.3969838619232178e-08\n",
"1.7229467630386353e-08\n",
"5.820766091346741e-08\n",
"1.7229467630386353e-08\n",
"4.7031790018081665e-08\n",
"4.330649971961975e-08\n",
"4.6100467443466187e-08\n",
"7.35744833946228e-08\n",
"3.213062882423401e-08\n",
"1.862645149230957e-09\n",
"1.862645149230957e-09\n",
"4.656612873077393e-10\n",
"4.6100467443466187e-08\n",
"1.862645149230957e-09\n",
"1.3969838619232178e-08\n",
"1.3969838619232178e-09\n",
"1.5366822481155396e-08\n"
]
}
],
"source": [
"# Here is an example of generating phase equilibrium data\n",
"# from an isobaric trace with the Peng-Robinson EOS and\n",
"# then calculating the residuals with the new routine. \n",
"#\n",
"# In this case no iteration\n",
"# is required since the residuals should all be close to zero\n",
"# because polishing is enabled by default in the tracing routine\n",
"import numpy as np \n",
"import pandas\n",
"\n",
"import teqp\n",
"from teqp import phaseequil as pe\n",
"\n",
"# Get a good result from the tracing\n",
"Tc_K = [190.564, 154.581]\n",
"pc_Pa = [4599200, 5042800]\n",
"acentric = [0.011, 0.022]\n",
"model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)\n",
"T = 170.0 # [K] # Note: above Tc of the second component\n",
"rhoL0, rhoV0 = model.superanc_rhoLV(T, ifluid=0) # start off at pure of the first component\n",
"p0 = rhoL0*model.get_R(np.array([1.0,0]))*T*(1+model.get_Ar01(T, rhoL0, np.array([1.0,0])))\n",
"j = model.trace_VLE_isobar_binary(p0, T, np.array([rhoL0, 0]), np.array([rhoV0, 0]))\n",
"df = pandas.DataFrame(j) # Now as a data frame\n",
"\n",
"for ir, row in df.iterrows():\n",
" # Only do every fifth, and skip the first because it is infinite dilution\n",
" # which requires special treatment\n",
" if (ir+1)%5 != 0: continue \n",
" \n",
" # The initial values of the variables\n",
" T = row['T / K']\n",
" rhovecs = [row['rhoL / mol/m^3'], row['rhoV / mol/m^3']]\n",
" betas = [1.0, 0.0]\n",
" unpacked = pe.UnpackedVariables(T, rhovecs, betas)\n",
" \n",
" zbulk = rhovecs[0]/np.sum(rhovecs[0])\n",
" specs = [\n",
" pe.TSpecification(T),\n",
" pe.BetaSpecification(1.0, 0) # Bubble point calculation, first phase (liquid) with index 0 is the whole mixture\n",
" # or you could consider instead to specify pressure, or ...\n",
" ]\n",
" gpe = pe.GeneralizedPhaseEquilibrium(model, zbulk, unpacked, specs)\n",
" \n",
" Xinit = unpacked.pack()\n",
" gpe.call(Xinit)\n",
" print(np.max(np.abs(gpe.res.r)))"
]
}
],
"metadata": {
"celltoolbar": "Raw Cell Format",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
1 change: 1 addition & 0 deletions doc/source/algorithms/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ Algorithms
:maxdepth: 2
:caption: Contents:

GeneralizedPhaseEquilibria
VLE
VLLE
VLLE-p
Expand Down
Loading

0 comments on commit f5fe546

Please sign in to comment.