ChEn-3170: Computational Methods in Chemical Engineering Fall 2018 UMass Lowell; Prof. V. F. de Almeida 03Dec2018
$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}} \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}} \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}} \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}} \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}} \newcommand{\Jmtrx}{\boldsymbol{\mathsf{J}}} \newcommand{\Smtrx}{\boldsymbol{\mathsf{S}}} \newcommand{\Xmtrx}{\boldsymbol{\mathsf{X}}} \newcommand{\Kmtrx}{\boldsymbol{\mathsf{K}}} \newcommand{\xvec}{\boldsymbol{x}} \newcommand{\avec}{\boldsymbol{\mathsf{a}}} \newcommand{\bvec}{\boldsymbol{\mathsf{b}}} \newcommand{\cvec}{\boldsymbol{\mathsf{c}}} \newcommand{\rvec}{\boldsymbol{\mathsf{r}}} \newcommand{\mvec}{\boldsymbol{\mathsf{m}}} \newcommand{\gvec}{\boldsymbol{\mathsf{g}}} \newcommand{\zerovec}{\boldsymbol{\mathsf{0}}} \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert} \newcommand{\abs}[1]{\left\lvert{#1}\right\rvert} \newcommand{\transpose}[1]{{#1}^\top} \DeclareMathOperator{\rank}{rank} \DeclareMathOperator{\gradx}{\nabla\!_{\xvec}} \newcommand{\Kcal}{\mathcal{K}} \newcommand{\Kcalvec}{\boldsymbol{\mathcal{K}}} \newcommand{\epsvec}{\boldsymbol{\varepsilon}} $
It is often asked how the chemical equilibrium species molar fractions involved in a reaction mechanism vary with temperature. The answer requires the evaluation of the thermodynamic equilibrium constant, $K_{\text{a},i}$, for a given reaction using the van't Hoff equation
\begin{equation*} \text{d}_T \ln K_{\text{a},i}= \frac{\Delta_\text{rxn} H_i(T)}{R\,T^2} \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m, \end{equation*}as a function of temperature using the definition
\begin{equation*} K_{\text{a},i}(T) := \exp{\biggl(\frac{-\Delta_\text{rxn} G^\circ_i(T)}{R\,T}\biggr)} \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m, \end{equation*}where $\Delta_\text{rxn} G_i^\circ$ is the standard state Gibbs energy change on the reaction. The standard state is the state of aggregation of the species involved in the reaction at the desired temperature and pressure of 1 bar. Here $K_{\text{a},i}$, by design, is a function of $T$ and the standard state (or state of aggregation) but not a function of pressure.
Integration of the van't Hoff equation with respect to temperature from the standard state $T=298.15$ K to any temperature $T$ gives
\begin{align*} \ln K_{\text{a},i}(T) - \ln K_{\text{a},i}(298.15) &= \int\limits_{298.15}^{T} \frac{\Delta_\text{rxn} H_i(T')}{R\,T'^2} \, \text{d}T' \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m, \\ \ln K_{\text{a},i}(T) &= \ln K_{\text{a},i}(298.15) + \int\limits_{298.15}^{T} \frac{\Delta_\text{rxn} H_i(T')}{R\,T'^2} \, \text{d}T' \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m, \\ K_{\text{a},i}(T) &= K_{\text{a},i}(298.15) \, \exp\biggl( \int\limits_{298.15}^{T} \frac{\Delta_\text{rxn} H_i(T')}{R\,T'^2} \, \text{d}T' \biggr) \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m, \\ K_{\text{a},i}(T) &= \exp{\biggl(\frac{-\Delta_\text{rxn} G^\circ_i(298.15)}{R\,298.15}\biggr)} \, \exp\biggl( \int\limits_{298.15}^{T} \frac{\Delta_\text{rxn} H_i(T')}{R\,T'^2} \, \text{d}T' \biggr) \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m . \end{align*}'''Equilibrium constant K_ai(T)'''
def k_a_cte( temp, delta_rxn_g_std_state_vec, delta_rxn_h_std_state_func_vec ):
assert temp > 0.0
n_reactions = delta_rxn_g_std_state_vec.size
assert delta_rxn_h_std_state_func_vec.size == n_reactions
from scipy.constants import gas_constant as r_gas # J/mol/K
from scipy.constants import calorie as calorie
r_gas /= calorie # cal/mol/K
import math
import numpy as np
k_a_298_vec = np.exp(-delta_rxn_g_std_state_vec/r_gas/298.15)
k_a_cte_vec = np.zeros(n_reactions)
import scipy.integrate as integrate
for i in range(n_reactions):
# quad returns a pair (integral,accuracy)
integral = integrate.quad(lambda T: delta_rxn_h_std_state_func_vec[i](T)/r_gas/T/T, 298.15, temp)
integral_error = integral[1]
assert abs(integral_error) <= 1e-5
k_a_cte_vec[i] = k_a_298_vec[i] * math.exp( integral[0] )
return k_a_cte_vec
The statement of chemical equilibrium obtained from the energy and entropy balance for a closed system is
\begin{equation*} K_{\text{a},i}(T) = \prod\limits_{j=1}^{n}\, a_j^{S_{i,j}} \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m , \end{equation*}where $a_j$ is the standard-state activity of the $j$th species which is a function of temperature, pressure, and standard state. For low and moderate pressure gas mixtures the activity is approximated by
\begin{equation*} a_j = \frac{x_j\,P}{P^\circ} , \end{equation*}where $x_j$ is the species molar fraction and $P^\circ = 1 \text{bar}$ is the standard state pressure. In this limit of pressure, the chemical equilibrium gives
\begin{equation*} K_{\text{a},i}(T) = \Bigl({\frac{P}{P^\circ}}\Bigr)^{\sum_j S_{i,j}} \, \prod\limits_{j=1}^{n}\, x_j^{S_{i,j}} \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m . \end{equation*}Again the left side of the equilibrium statement is only a function of temperature. Pressure is accounted for in the equilibrium through the activities of the species. If the $i$th reaction conserves moles then $\sum_j S_{i,j}=0$ and theres is no depedency on pressure at low to moderate pressures.
Equivalently, the chemical equilibrium statement can be written as
\begin{equation*} \Bigl({\frac{P}{P^\circ}}\Bigr)^{-\sum_j S_{i,j}} \, K_{\text{a},i}(T) = \prod\limits_{j=1}^{n}\, x_j^{S_{i,j}} \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m , \end{equation*}where using the molar fraction equilibrium function $K_{x,i}(T,P)$ produces the form
\begin{equation*} K_{x,i}(T,P) = \prod\limits_{j=1}^{n}\, x_j^{S_{i,j}} \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m , \end{equation*}with the definition
\begin{equation*} K_{x,i}(T,P) := \Bigl({\frac{P}{P^\circ}}\Bigr)^{-\sum_j S_{i,j}} \, K_{\text{a},i}(T) . \end{equation*}'''Molar equilibrium constant K_xi(T,P)'''
def k_x_cte( pressure, stoic_mtrx, k_a_cte_vec ):
assert temp > 0.0
assert pressure > 0.0
n_reactions = stoic_mtrx.shape[0]
assert k_a_cte_vec.size == n_reactions
import numpy as np
k_x_cte_vec = np.zeros(n_reactions)
for i in range(n_reactions):
k_x_cte_vec[i] = pressure**(-stoic_mtrx[i,:].sum()) * k_a_cte_vec[i]
return k_x_cte_vec
Determine the equilibrium compositions at pressures of 1 bar, 2 bar, and 3 bar, and varying temperature from 800 K to 2000 K for a pure stream of ethane. Discuss the results for simultaneous production of ethylene and acetylene from ethane.
Ethylene and Acetylene from Ethane | $\Delta_\text{rxn} H^\circ(T)$ [cal/mol] (T in K) | $\Delta_\text{rxn} G^\circ(298.15 K)$ [cal/mol] |
---|---|---|
C2H6 <=> C2H4 + H2 | $31\,094 + 6.101\,T - 1.46\times 10^{-3}\,T^2 -2.222\times 10^{-6}\,T^3 + 9.843\times 10^{-10}\,T^4$ | 24 142 |
C2H6 <=> C2H2 + 2 H2 | $69\,936 + 18.082\,T - 10^{-2}\,T^2 -6.617\times 10^{-7}\,T^3 + 1.236\times 10^{-9}\,T^4$ | 57 860 |
C2H4 <=> C2H2 + H2 | $38\,842 + 11.981\,T - 8.546\times 10^{-3}\,T^2 +1.561\times 10^{-6}\,T^3 + 2.517\times 10^{-10}\,T^4$ | 33 718 |
Through Newton's method the roots of multiple non-linear equations, $\Kcal\bigl(\xvec(\widehat{\epsvec})\bigr)=0$ can be computed for varying pressures and temperatures. In this case we write the equilibrium functions as a vector with components
\begin{equation*} \Kcal_i\bigl(\xvec\bigr)= K_{x,i}(T,P) - \prod\limits_{j=1}^n\, x_j^{S_{i,j}} = 0, \ \ \ \ \ \qquad \forall \ \ \ \ \ \qquad i=1,\ldots,m \end{equation*}where $i$ is the index of each reaction.
'''Equilibrium function vector'''
def keq_function( x_vec, k_x_cte_vec, stoic_mtrx ):
n_reactions = k_x_cte_vec.size
n_species = x_vec.size
assert n_reactions == stoic_mtrx.shape[0]
assert n_species == stoic_mtrx.shape[1]
prod_vec = np.ones(n_reactions,dtype=np.float64)
for i in range(n_reactions):
for j in range(n_species):
prod_vec[i] *= x_vec[j]**stoic_mtrx[i,j]
keq_vec = k_x_cte_vec - prod_vec
return keq_vec
The molar fractions vector is a function of the normalized extent of reactions vector
\begin{equation*} \xvec(\widehat{\epsvec}) = \frac{\xvec^{(0)} + \Smtrx^\top\widehat{\epsvec}\,}{1+\sum\limits_j (\Smtrx^\top\widehat{\epsvec})_j} , \end{equation*}where $\xvec^{(0)}$ is the vector of initial molar fractions (or reference).
'''Molar fractions function'''
def molar_fractions( ext_hat_vec, x_vec_0, stoic_mtrx ):
import numpy as np
assert abs(x_vec_0.sum()-1.0) <= 1e-10
assert np.all(x_vec_0>=0.0)
assert np.all(np.abs(ext_hat_vec)>=0.0),'ext_hat_vec=%r'%ext_hat_vec
denom = 1.0 + (stoic_mtrx.transpose()@ext_hat_vec).sum()
assert abs( denom ) >= 1e-9,'denom=%r; ext_hat_vec=%r; x_vec_0=%r'%(denom,ext_hat_vec,x_vec_0)
x_vec = ( x_vec_0 + stoic_mtrx.transpose()@ext_hat_vec ) / denom
return x_vec
The usage of Newton's method to compute an equilibrium molar fraction solution vector requires the value of $\Kcalvec\bigl(\xvec(\widehat{\epsvec})\bigr)$ at different values of $\widehat{\epsvec}$ and the total derivative
\begin{equation*} d_{\widehat{\epsvec}} \Kcalvec\bigl(\xvec(\widehat{\epsvec})\bigr) = \gradx\Kcalvec \, d_{\widehat{\epsvec}}\xvec , \end{equation*}which is matrix product. The molar fraction gradient is computed as a product of three matrices
\begin{equation*} \gradx\Kcal = - \Kmtrx\,\Smtrx\,\Xmtrx^{-1} , \end{equation*}that is,
\begin{equation*} \Kmtrx = \begin{pmatrix} K_{x,1} & 0 & \dots & 0 \\ 0 & K_{x,2} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & K_{x,m} \end{pmatrix}, \ \ \ \Smtrx = \begin{pmatrix} S_{1,1} & S_{1,2} & \dots & S_{1,n} \\ S_{2,1} & S_{2,2} & \dots & S_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ S_{m,1} & S_{m,2} & \dots & S_{m,n} \end{pmatrix}, \ \ \ \Xmtrx = \begin{pmatrix} x_1 & 0 & \dots & 0 \\ 0 & x_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & x_n \end{pmatrix} . % \end{equation*}'''Gradient wrt to molar fractions of the equilibrium function'''
def grad_x_k_function( x_vec, keq_cte_vec, stoic_mtrx ):
import numpy as np
# build the K matrix
k_mtrx = np.diag(keq_cte_vec)
# correct for division by a very small number
import numpy as np
x_vec_local = np.copy(x_vec)
max_x = x_vec_local.max()
for i in range(x_vec_local.size):
if x_vec_local[i] < 1e-8:
x_vec_local[i] = max_x # some reasonable mole fraction
x_vec_inv = 1.0/x_vec_local
x_mtrx_inv = np.diag(x_vec_inv)
grad_mtrx = - k_mtrx @ stoic_mtrx @ x_mtrx_inv
return grad_mtrx
The molar fraction derivative is
\begin{equation*} d_{\widehat{\epsvec}}\xvec = \frac{\Smtrx^\top}{1+\sum_j (\Smtrx^\top\widehat{\epsvec})_j} - \frac{\bigl(\xvec^{(0)} + \Smtrx^\top\widehat{\epsvec}\bigr)\otimes\bigl(\sum_j\Smtrx^\top_{j,\bullet}\bigr)}{\bigl(1+\sum_j (\Smtrx^\top\widehat{\epsvec})_j\bigr)^2} \end{equation*}'''Derivative of the molar fractions function wrt normalized extent of reaction'''
def d_ext_molar_fractions( ext_hat_vec, x_vec_0, stoic_mtrx ):
denom = 1.0 + (stoic_mtrx.transpose()@ext_hat_vec).sum()
denom2 = denom**2
dext_x_vec = stoic_mtrx.transpose() / denom \
- np.outer( x_vec_0 + stoic_mtrx.transpose()@ext_hat_vec, (stoic_mtrx.transpose()).sum(0) ) / \
denom2
return dext_x_vec
Given $\Kcalvec\bigl(\xvec(\widehat{\epsvec})\bigr)$, find the root
\begin{equation*} \Kcalvec\bigl(\xvec(\widehat{\epsvec})\bigr) = 0 , \end{equation*}using an iterative method based on the initial guess $\widehat{\epsvec}_0$. Compute the updates
\begin{equation*} d_{\widehat{\epsvec}} \Kcalvec\bigl(\xvec(\widehat{\epsvec})\bigr) \, \delta \widehat{\epsvec}_k = - \Kcalvec\bigl(\xvec(\widehat{\epsvec}_{k-1})\bigr) \ \qquad \ \forall \ \qquad \ k = 1,\ldots,k_\text{max} , \end{equation*}then compute the approximation to the root
\begin{equation*} \widehat{\epsvec}_k = \widehat{\epsvec}_{k-1} + \delta \widehat{\epsvec}_k \ \qquad \ \forall \ \qquad\ \ k = 1,\ldots,k_\text{max} , \end{equation*}until convergence, say, $\norm{\delta\widehat{\epsvec}_k} \le 10^{-8}$ and $\norm{\Kcal\bigl(\xvec(\widehat{\epsvec}_k)\bigr)} \le 10^{-8}$, or no convergence achieved , say $k>k_\text{max}$.
"""Newton's method"""
def newton_solve( x_vec_0, k_x_cte_vec, stoic_mtrx,
ext_hat_vec_0=None, k_max=30, tolerance=1.0e-10, verbose=True ):
#from chen_3170.toolkit import solve
import numpy as np
if ext_hat_vec_0 is None:
ext_hat_vec_0 = np.zeros(k_x_cte_vec.size,dtype=np.float64)
# Other initialization
delta_vec_k = 1e+10 * np.ones(ext_hat_vec_0.size,dtype=np.float64)
keq_vec_k = 1e+10 * np.ones(ext_hat_vec_0.size,dtype=np.float64) # equilibrium function initial value
ext_hat_vec = np.copy(ext_hat_vec_0)
if verbose is True:
print('\n')
print('******************************************************')
print(" Newton's Method Iterations ")
print('******************************************************')
print("k | K(e_k) | K'(e_k) | |del e_k| | e_k |convg|")
print('------------------------------------------------------')
import math
k = 0
while (np.linalg.norm(delta_vec_k) > tolerance or np.linalg.norm(keq_vec_k) > tolerance) and k <= k_max:
# compute the molar fractions
x_vec = molar_fractions( ext_hat_vec, x_vec_0, stoic_mtrx )
# compute the equilibrium function
keq_vec_k = keq_function( x_vec, k_x_cte_vec, stoic_mtrx )
# copute the molar fraction gradient of the equilibrium function
grad_x_k = grad_x_k_function( x_vec, k_x_cte_vec, stoic_mtrx )
# compute the extent of reaction derivative of the molar fraction
d_ext_x = d_ext_molar_fractions( ext_hat_vec, x_vec_0, stoic_mtrx )
# form the total derivative of the equilibrium function wrt the extent of reaction
d_ext_keq_k = grad_x_k @ d_ext_x # Jacobian matrix
# compute the Newton update
delta_vec_k_old = delta_vec_k
a_mtrx = d_ext_keq_k
b_vec = - keq_vec_k
#delta_vec_k = solve(a_mtrx, b_vec, zero_tol=1e-8)
delta_vec_k = np.linalg.solve(a_mtrx, b_vec) # using scipy linear solver
# compute the update to the root candidate
ext_hat_vec += delta_vec_k
if k > 0:
if np.linalg.norm(delta_vec_k) != 0.0 and np.linalg.norm(delta_vec_k_old) != 0.0:
convergence_factor = math.log(np.linalg.norm(delta_vec_k),10) / math.log(np.linalg.norm(delta_vec_k_old),10)
else:
convergence_factor = 0.0
else:
convergence_factor = 0.0
k = k + 1
if verbose is True:
print('%2i %+5.3e %+5.3e %+5.3e %+5.3e %5.2f'%\
(k,np.linalg.norm(keq_vec_k),np.linalg.norm(d_ext_keq_k),np.linalg.norm(delta_vec_k),np.linalg.norm(ext_hat_vec),convergence_factor))
if verbose is True:
print('******************************************************')
print('Root = ',ext_hat_vec)
return ext_hat_vec
Using the ethylene-acetylene from ethane reaction mechanism from the input file: data/ethylene-acetylene-rxn2.txt
!cat data/ethylene-acetylene-rxn2.txt
# # Production of ethylene and acetylene from ethane # ..str:float:str C2H6 <=> C2H4 + H2 : DeltaGstd = 24142 : DeltaHstd(T) = 31094+6.101*T-1.46e-3*T**2-2.222e-6*T**3+9.843e-10*T**4 C2H6 <=> C2H2 + 2 H2 : DeltaGstd = 57860 : DeltaHstd(T) = 69936+18.082*T-1.0e-2*T**2-6.617e-7*T**3+1.236e-9*T**4 C2H4 <=> C2H2 + H2 : DeltaGstd = 33718 : DeltaHstd(T) = 38842+11.981*T-8.546e-3*T**2+1.561e-6*T**3+2.517e-10*T**4
'''Import the ethylene-acetylene from ethane reaction mechanism'''
try:
from chen_3170.toolkit import reaction_mechanism
except ModuleNotFoundError:
assert False, 'You need to provide your own reaction_mechanism function here. Bailing out.'
# read species, reactions, equilibrium constants and build the stoichiometric matrix
(species, reactions, stoic_mtrx, delta_rxn_g_std_state, delta_rxn_h_std_state) = reaction_mechanism('data/ethylene-acetylene-rxn2.txt',shuffle=False)
'''Info on the data'''
print('species=',species)
from chen_3170.help import print_reactions
print('')
print_reactions(reactions)
print('')
from chen_3170.help import plot_matrix
plot_matrix(stoic_mtrx, title='Stoichiometric Matrix')
print('')
import numpy as np
np.set_printoptions(precision=3,threshold=100,edgeitems=5)
print('stoic_mtrx=\n',stoic_mtrx)
species= ['H2', 'C2H2', 'C2H4', 'C2H6'] r0 : C2H6 <=> C2H4 + H2 r1 : C2H6 <=> C2H2 + 2 H2 r2 : C2H4 <=> C2H2 + H2 n_reactions = 3 matrix shape = (3, 4)
stoic_mtrx= [[ 1. 0. 1. -1.] [ 2. 1. 0. -1.] [ 1. 1. -1. 0.]]
'''Create reactions energy deltas'''
import numpy as np
# Gibb's energy of reaction
delta_rxn_g_std_state_vec = np.array(delta_rxn_g_std_state)
# Enthalpy of reaction functions
delta_rxn_h_std_state_func_lst = list() # hold list of functions
# use eval() to evaluate a string as an expression and returns the result of the evaluation
# here the result returned is a lambda function for the argument T
for string in delta_rxn_h_std_state:
delta_rxn_h_std_state_func_lst.append( eval('lambda T:'+string) )
# create a numpy vector from the list of lambda functions
delta_rxn_h_std_state_func_vec = np.array(delta_rxn_h_std_state_func_lst)
# usage:
# delta...[i](T)
# will return the evaluation of the ith function for the argument T
The following reference molar fraction $\xvec^{(0)}$ is given:
Reference Molar Fraction | Parameter | Value |
---|---|---|
C2H6 | $x_\text{C2H6}$ | 1.0 |
C2H4 | $x_\text{C2H4}$ | 0.0 |
C2H2 | $x_\text{C2H2}$ | 0.0 |
H2 | $x_\text{H2}$ | 0.0 |
'''Set reference molar fraction'''
import numpy as np
x_dict_0 = {'C2H6':1.0,'C2H4':0.0,'C2H2':0.0,'H2':0.0}
'''Check the rank of the stoichiometric matrix'''
try:
from chen_3170.toolkit import matrix_rank
except ModuleNotFoundError:
assert False, 'You need to provide your own matrix_rank function here. Bailing out.'
rank = matrix_rank( stoic_mtrx )
print('stoic_mtrx m x n =',stoic_mtrx.shape)
print('stoic_mtrx rank =',rank)
stoic_mtrx m x n = (3, 4) stoic_mtrx rank = 2
'''Build the full-rank sub-mechanism reactions list'''
try:
from chen_3170.toolkit import sub_mechanisms
except ModuleNotFoundError:
assert False, 'You need to provide your own sub_mechanisms function here. Bailing out.'
sub_mechanisms = sub_mechanisms( species, reactions, stoic_mtrx )
# reactions = 3 # species = 4 rank of S = 2 # of all possible sub_mechanisms = 3 # of full-rank sub_mechanisms = 3
'''Plot full-rank sub-mechanism reactions list'''
try:
from chen_3170.toolkit import plot_reaction_mechanisms
except ModuleNotFoundError:
assert False, 'You need to provide your own plot_reaction_mechanisms function here. Bailing out.'
plot_reaction_mechanisms( sub_mechanisms )
Fixed $T$, $P$ equilibrium using the first sub-mechanism.
pressure = 1.0 # bar
temp = 800 # K
Select a reaction sub-mechanism; any sub-mechanism will do it but the initial guess for the extent of reaction in Newton's method will be critical depending on the choice of sub-mechanism. Just by inspecting the reaction sub-mechanisms and given the pure initial stream of ethane, one can guess the difficulty of finding a good initial guess.
'''Select a sub-mechanism'''
sub_mechanism_id = 2
sub_mechanism = sub_mechanisms[sub_mechanism_id]
print('species=',species)
for i in sub_mechanism:
print(i)
species= ['H2', 'C2H2', 'C2H4', 'C2H6'] (1, 2) ['C2H6 <=> C2H2 + 2 H2', 'C2H4 <=> C2H2 + H2'] [[ 2. 1. 0. -1.] [ 1. 1. -1. 0.]] 10.0
'''Setup a sub-mechanism'''
reaction_ids = sub_mechanism[0]
reactions = sub_mechanism[1]
stoic_mtrx = sub_mechanism[2]
print('species=',species)
print('')
print('reaction ids=',reaction_ids)
print('')
print_reactions(reactions)
print('')
import numpy as np
np.set_printoptions(precision=3,threshold=100,edgeitems=5)
print('stoic_mtrx=\n',stoic_mtrx)
delta_rxn_g_std_state_vec = delta_rxn_g_std_state_vec[np.array(reaction_ids)]
delta_rxn_h_std_state_func_vec = delta_rxn_h_std_state_func_vec[np.array(reaction_ids)]
print('')
print('delta_rxn_g_std_state_vec=',delta_rxn_g_std_state_vec)
print('')
print('delta_rxn_h_std_state=',delta_rxn_h_std_state)
species= ['H2', 'C2H2', 'C2H4', 'C2H6'] reaction ids= (1, 2) r0 : C2H6 <=> C2H2 + 2 H2 r1 : C2H4 <=> C2H2 + H2 n_reactions = 2 stoic_mtrx= [[ 2. 1. 0. -1.] [ 1. 1. -1. 0.]] delta_rxn_g_std_state_vec= [57860. 33718.] delta_rxn_h_std_state= ['31094+6.101*T-1.46e-3*T**2-2.222e-6*T**3+9.843e-10*T**4', '69936+18.082*T-1.0e-2*T**2-6.617e-7*T**3+1.236e-9*T**4', '38842+11.981*T-8.546e-3*T**2+1.561e-6*T**3+2.517e-10*T**4']
'''Setup a sub-mechanism'''
x_vec_0 = np.zeros(len(species),dtype=np.float64)
for i in range(len(species)):
x_vec_0[i] = x_dict_0[species[i]]
assert abs(np.sum(x_vec_0) - 1.0) <= 1e-12
assert np.all(x_vec_0 >=0.0)
k_a_cte_vec = k_a_cte( temp, delta_rxn_g_std_state_vec, delta_rxn_h_std_state_func_vec )
print('k_a_cte_vec =',k_a_cte_vec)
k_x_cte_vec = k_x_cte( pressure, stoic_mtrx, k_a_cte_vec )
print('k_x_cte_vec =',k_x_cte_vec)
k_a_cte_vec = [3.581e-08 7.852e-06] k_x_cte_vec = [3.581e-08 7.852e-06]
Plot equilibrium vector function for the first sub-mechanism.
'''Function: plot equilibrium function'''
def plot_keq_function( ext_hat_vec_min, ext_hat_vec_max, num_pts,
x_vec_0, k_x_cte_vec, stoic_mtrx,
ext_hat_vec_root=None ):
import matplotlib.pyplot as plt
plt.figure(1, figsize=(25, 5))
# number of functions to be plotted
num_keq_functions = k_x_cte_vec.size
import numpy as np
# allocate a matrix for all ext_hat_vec plot points (min to max)
ext_hat_plot_pts = np.zeros((num_pts,num_keq_functions))
# create all plot points for ext_hat_vec
for i in range(num_keq_functions):
ext_hat_plot_pts[:,i] = np.linspace(ext_hat_vec_min[i], ext_hat_vec_max[i], num_pts)
# allocate a matrix for all keq_function_vec plot points
keq_function_plot_values = np.zeros((num_pts,num_keq_functions))
# create all plot points for keq_function
for i in range(num_keq_functions):
for k in range(num_pts):
x_vec = molar_fractions( ext_hat_plot_pts[k,:], x_vec_0, stoic_mtrx )
keq_function_plot_values[k,:] = keq_function( x_vec, k_x_cte_vec, stoic_mtrx )
# Plot keq_function_vec in subplots
n_rows = 1
n_columns = num_keq_functions
# create colors for each equilibrium function
from chen_3170.help import color_map
colors = color_map(num_keq_functions)
for iplot in range(num_keq_functions):
plt.subplot(n_rows,n_columns, iplot+1)
color = colors[iplot]
plt.plot(ext_hat_plot_pts[:,iplot], keq_function_plot_values[:,iplot],
color=color,label='$K_{x,%i}=$%8.2e'%(reaction_ids[iplot],k_x_cte_vec[iplot]))
plt.xlabel(r'$\hat{\varepsilon}_{%i}$'%reaction_ids[iplot],fontsize=18)
plt.ylabel(r'$K_{%i}(\hat{\varepsilon})$'%reaction_ids[iplot],fontsize=18)
plt.title(reactions[iplot],fontsize=20)
plt.legend(loc='best',fontsize=12)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
(x_min,x_max) = plt.xlim()
dx = abs(x_max-x_min)
x_text = x_max+dx*1/100
(y_min,y_max) = plt.ylim()
dy = abs(y_max-y_min)
y_text = y_max - dy*0.1
for x,spc in zip(x_vec_0,species):
plt.text(x_text, y_text, r'$x_{%s}^{(0)}=$%8.2e'%(spc,x),fontsize=16)
y_text -= dy*0.15
if ext_hat_vec_root is not None:
e_root = ext_hat_vec_root[iplot]
plt.plot(e_root, 0.0,'r*',label='root',markersize=14)
(x_min,x_max) = plt.xlim()
dx = abs(x_max-x_min)
x_text = e_root + dx*0.01
(y_min,y_max) = plt.ylim()
dy = abs(y_max-y_min)
y_text = 0.0 + dy*0.01
plt.text(x_text, y_text, r'$\hat{\varepsilon}^*=$%8.2e'%e_root,fontsize=16)
plt.grid(True)
plt.show()
print('')
return
'''Plot equilibrium vector function'''
import numpy as np
num_keq_functions = k_x_cte_vec.size
ext_hat_vec_min = 1e-3*np.ones(num_keq_functions,dtype=np.float64)
ext_hat_vec_max = 1*np.ones(num_keq_functions,dtype=np.float64)
#ext_hat_vec_max[0] = -.49
#ext_hat_vec_max[1] = -0.55
n_pts = 100
plot_keq_function( ext_hat_vec_min, ext_hat_vec_max, n_pts,
x_vec_0, k_x_cte_vec, stoic_mtrx )
/Users/valmor_dealmeida/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:15: RuntimeWarning: divide by zero encountered in double_scalars from ipykernel import kernelapp as app
Solve for equilibrium for the first sub-mechanism.
'''Find root and equilibrium molar fractions'''
# sub-mechanism 0
if sub_mechanism_id == 0:
ext_hat_vec_0 = 1e-4 * np.ones(len(reactions),dtype=np.float64)
# sub-mechanism 1
if sub_mechanism_id == 1:
ext_hat_vec_0 = 1e-6 * np.ones(len(reactions),dtype=np.float64)
ext_hat_vec_0[1] = 0.0
# sub-mechanism 2
if sub_mechanism_id == 2:
ext_hat_vec_0 = 1e-5 * np.ones(len(reactions),dtype=np.float64)
ext_hat_vec_0[1] = -1e-5
k_max = 50
tolerance = 1.0e-8
ext_hat_vec = newton_solve( x_vec_0, k_x_cte_vec, stoic_mtrx,
ext_hat_vec_0, k_max, tolerance )
x_vec = molar_fractions( ext_hat_vec, x_vec_0, stoic_mtrx )
print('')
print('(T[K],P[bar])=(%4.1f,%4.1f)\n'%(temp,pressure))
print('Equilibrium mole fractions:\n')
for (x,spc) in zip(x_vec,species):
print('x_%s = %5.3e\n'%(spc,x))
assert np.all(x_vec >= 0.0)
assert abs(x_vec.sum() - 1.0) <= 1e-10
****************************************************** Newton's Method Iterations ****************************************************** k | K(e_k) | K'(e_k) | |del e_k| | e_k |convg| ------------------------------------------------------ 1 +7.852e-06 +2.221e+00 +5.000e-06 +1.118e-05 0.00 2 +7.148e-06 +4.503e+00 +1.013e-05 +2.105e-05 0.94 3 +1.577e-06 +3.383e+00 +1.319e-05 +3.416e-05 0.98 4 +7.330e-07 +2.756e+00 +1.959e-05 +5.369e-05 0.96 5 +4.720e-07 +2.320e+00 +2.974e-05 +8.338e-05 0.96 6 +3.441e-07 +2.020e+00 +4.516e-05 +1.285e-04 0.96 7 +2.545e-07 +1.817e+00 +6.827e-05 +1.967e-04 0.96 8 +1.849e-07 +1.682e+00 +1.028e-04 +2.996e-04 0.96 9 +1.323e-07 +1.592e+00 +1.546e-04 +4.541e-04 0.96 10 +9.491e-08 +1.532e+00 +2.321e-04 +6.863e-04 0.95 11 +6.979e-08 +1.493e+00 +3.483e-04 +1.035e-03 0.95 12 +5.407e-08 +1.466e+00 +5.225e-04 +1.557e-03 0.95 13 +4.500e-08 +1.448e+00 +7.837e-04 +2.341e-03 0.95 14 +4.019e-08 +1.435e+00 +1.175e-03 +3.516e-03 0.94 15 +3.780e-08 +1.426e+00 +1.761e-03 +5.277e-03 0.94 16 +3.662e-08 +1.419e+00 +2.636e-03 +7.913e-03 0.94 17 +3.598e-08 +1.413e+00 +3.935e-03 +1.185e-02 0.93 18 +3.544e-08 +1.407e+00 +5.838e-03 +1.769e-02 0.93 19 +3.467e-08 +1.400e+00 +8.543e-03 +2.623e-02 0.93 20 +3.315e-08 +1.391e+00 +1.213e-02 +3.835e-02 0.93 21 +3.005e-08 +1.378e+00 +1.607e-02 +5.443e-02 0.94 22 +2.417e-08 +1.363e+00 +1.834e-02 +7.276e-02 0.97 23 +1.497e-08 +1.346e+00 +1.517e-02 +8.793e-02 1.05 24 +5.336e-09 +1.332e+00 +6.526e-03 +9.446e-02 1.20 25 +6.237e-10 +1.326e+00 +8.188e-04 +9.528e-02 1.41 26 +8.155e-12 +1.326e+00 +1.080e-05 +9.529e-02 1.61 27 +1.386e-15 +1.326e+00 +1.836e-09 +9.529e-02 1.76 ****************************************************** Root = [ 0.067 -0.067] (T[K],P[bar])=(800.0, 1.0) Equilibrium mole fractions: x_H2 = 6.314e-02 x_C2H2 = 7.850e-06 x_C2H4 = 6.312e-02 x_C2H6 = 8.737e-01
'''Plot equilibrium vector function'''
import numpy as np
n_keq_functions = k_x_cte_vec.size
ext_hat_vec_min = -1e-4*np.ones(n_keq_functions,dtype=np.float64)
ext_hat_vec_max = ext_hat_vec
n_pts = 100
plot_keq_function( ext_hat_vec_min, ext_hat_vec_max, n_pts,
x_vec_0, k_x_cte_vec, stoic_mtrx, ext_hat_vec )
Varying temperature from $800$ K to $1800$ K for three values of pressure: $1$ bar, $2$ bar, and $3$ bar.
'''Loop over pressure and temperature'''
pressure_pts = [1,2,3]
# storage space
npts = 50
temp_pts = np.linspace(800,1800,npts)
x_vec_block = np.zeros((3,npts,len(species)))
k_x_cte_block = np.zeros((3,npts,delta_rxn_g_std_state_vec.size))
delta_rxn_h_std_state_mtrx = np.zeros((npts,delta_rxn_g_std_state_vec.size))
k_max = 50
tolerance = 1.0e-8
for p_pt in range(len(pressure_pts)):
pressure = pressure_pts[p_pt]
# initial guesses
# sub-mechanism 0
if sub_mechanism_id == 0:
ext_hat_vec_0 = 1e-4 * np.ones(len(reactions),dtype=np.float64)
# sub-mechanism 1
if sub_mechanism_id == 1:
ext_hat_vec_0 = 1e-6 * np.ones(len(reactions),dtype=np.float64)
ext_hat_vec_0[1] = 0.0
# sub-mechanism 2
if sub_mechanism_id == 2:
ext_hat_vec_0 = 1e-5 * np.ones(len(reactions),dtype=np.float64)
ext_hat_vec_0[1] = -1e-5
for t_pt in range(npts):
temp = temp_pts[t_pt]
# for plotting
delta_rxn_h_std_state_mtrx[t_pt,0] = delta_rxn_h_std_state_func_vec[0](T=temp)
delta_rxn_h_std_state_mtrx[t_pt,1] = delta_rxn_h_std_state_func_vec[1](T=temp)
k_a_cte_vec = k_a_cte( temp, delta_rxn_g_std_state_vec, delta_rxn_h_std_state_func_vec )
k_x_cte_vec = k_x_cte( pressure, stoic_mtrx, k_a_cte_vec )
k_x_cte_block[p_pt,t_pt,:] = k_x_cte_vec # for plotting
# solve
ext_hat_vec = newton_solve( x_vec_0, k_x_cte_vec, stoic_mtrx,
ext_hat_vec_0, k_max, tolerance, verbose=False )
x_vec = molar_fractions( ext_hat_vec, x_vec_0, stoic_mtrx )
print('')
print('(T[K],P[bar])=(%4.1f,%4.1f)\n'%(temp,pressure))
print('Equilibrium mole fractions:\n')
for (x,spc) in zip(x_vec,species):
print('x_%s = %5.3e\n'%(spc,x))
assert np.all(x_vec >= 0.0),'x_vec=%r'%x_vec
assert abs(x_vec.sum() - 1.0) <= 1e-10
x_vec_block[p_pt,t_pt,:] = x_vec
ext_hat_vec_0 = ext_hat_vec
(T[K],P[bar])=(800.0, 1.0) Equilibrium mole fractions: x_H2 = 6.314e-02 x_C2H2 = 7.850e-06 x_C2H4 = 6.312e-02 x_C2H6 = 8.737e-01 (T[K],P[bar])=(820.4, 1.0) Equilibrium mole fractions: x_H2 = 8.089e-02 x_C2H2 = 1.560e-05 x_C2H4 = 8.086e-02 x_C2H6 = 8.382e-01 (T[K],P[bar])=(840.8, 1.0) Equilibrium mole fractions: x_H2 = 1.018e-01 x_C2H2 = 2.998e-05 x_C2H4 = 1.018e-01 x_C2H6 = 7.964e-01 (T[K],P[bar])=(861.2, 1.0) Equilibrium mole fractions: x_H2 = 1.260e-01 x_C2H2 = 5.589e-05 x_C2H4 = 1.259e-01 x_C2H6 = 7.481e-01 (T[K],P[bar])=(881.6, 1.0) Equilibrium mole fractions: x_H2 = 1.531e-01 x_C2H2 = 1.013e-04 x_C2H4 = 1.529e-01 x_C2H6 = 6.938e-01 (T[K],P[bar])=(902.0, 1.0) Equilibrium mole fractions: x_H2 = 1.829e-01 x_C2H2 = 1.786e-04 x_C2H4 = 1.826e-01 x_C2H6 = 6.343e-01 (T[K],P[bar])=(922.4, 1.0) Equilibrium mole fractions: x_H2 = 2.147e-01 x_C2H2 = 3.072e-04 x_C2H4 = 2.141e-01 x_C2H6 = 5.709e-01 (T[K],P[bar])=(942.9, 1.0) Equilibrium mole fractions: x_H2 = 2.477e-01 x_C2H2 = 5.161e-04 x_C2H4 = 2.467e-01 x_C2H6 = 5.051e-01 (T[K],P[bar])=(963.3, 1.0) Equilibrium mole fractions: x_H2 = 2.809e-01 x_C2H2 = 8.478e-04 x_C2H4 = 2.792e-01 x_C2H6 = 4.390e-01 (T[K],P[bar])=(983.7, 1.0) Equilibrium mole fractions: x_H2 = 3.133e-01 x_C2H2 = 1.363e-03 x_C2H4 = 3.106e-01 x_C2H6 = 3.748e-01 (T[K],P[bar])=(1004.1, 1.0) Equilibrium mole fractions: x_H2 = 3.439e-01 x_C2H2 = 2.148e-03 x_C2H4 = 3.396e-01 x_C2H6 = 3.143e-01 (T[K],P[bar])=(1024.5, 1.0) Equilibrium mole fractions: x_H2 = 3.721e-01 x_C2H2 = 3.318e-03 x_C2H4 = 3.654e-01 x_C2H6 = 2.592e-01 (T[K],P[bar])=(1044.9, 1.0) Equilibrium mole fractions: x_H2 = 3.972e-01 x_C2H2 = 5.027e-03 x_C2H4 = 3.872e-01 x_C2H6 = 2.106e-01 (T[K],P[bar])=(1065.3, 1.0) Equilibrium mole fractions: x_H2 = 4.193e-01 x_C2H2 = 7.471e-03 x_C2H4 = 4.044e-01 x_C2H6 = 1.688e-01 (T[K],P[bar])=(1085.7, 1.0) Equilibrium mole fractions: x_H2 = 4.385e-01 x_C2H2 = 1.089e-02 x_C2H4 = 4.167e-01 x_C2H6 = 1.339e-01 (T[K],P[bar])=(1106.1, 1.0) Equilibrium mole fractions: x_H2 = 4.552e-01 x_C2H2 = 1.557e-02 x_C2H4 = 4.240e-01 x_C2H6 = 1.053e-01 (T[K],P[bar])=(1126.5, 1.0) Equilibrium mole fractions: x_H2 = 4.698e-01 x_C2H2 = 2.182e-02 x_C2H4 = 4.262e-01 x_C2H6 = 8.217e-02 (T[K],P[bar])=(1146.9, 1.0) Equilibrium mole fractions: x_H2 = 4.831e-01 x_C2H2 = 2.992e-02 x_C2H4 = 4.232e-01 x_C2H6 = 6.377e-02 (T[K],P[bar])=(1167.3, 1.0) Equilibrium mole fractions: x_H2 = 4.955e-01 x_C2H2 = 4.015e-02 x_C2H4 = 4.152e-01 x_C2H6 = 4.923e-02 (T[K],P[bar])=(1187.8, 1.0) Equilibrium mole fractions: x_H2 = 5.074e-01 x_C2H2 = 5.265e-02 x_C2H4 = 4.021e-01 x_C2H6 = 3.781e-02 (T[K],P[bar])=(1208.2, 1.0) Equilibrium mole fractions: x_H2 = 5.193e-01 x_C2H2 = 6.745e-02 x_C2H4 = 3.844e-01 x_C2H6 = 2.889e-02 (T[K],P[bar])=(1228.6, 1.0) Equilibrium mole fractions: x_H2 = 5.312e-01 x_C2H2 = 8.438e-02 x_C2H4 = 3.625e-01 x_C2H6 = 2.195e-02 (T[K],P[bar])=(1249.0, 1.0) Equilibrium mole fractions: x_H2 = 5.433e-01 x_C2H2 = 1.031e-01 x_C2H4 = 3.371e-01 x_C2H6 = 1.657e-02 (T[K],P[bar])=(1269.4, 1.0) Equilibrium mole fractions: x_H2 = 5.553e-01 x_C2H2 = 1.231e-01 x_C2H4 = 3.091e-01 x_C2H6 = 1.242e-02 (T[K],P[bar])=(1289.8, 1.0) Equilibrium mole fractions: x_H2 = 5.673e-01 x_C2H2 = 1.438e-01 x_C2H4 = 2.797e-01 x_C2H6 = 9.247e-03 (T[K],P[bar])=(1310.2, 1.0) Equilibrium mole fractions: x_H2 = 5.789e-01 x_C2H2 = 1.645e-01 x_C2H4 = 2.498e-01 x_C2H6 = 6.833e-03 (T[K],P[bar])=(1330.6, 1.0) Equilibrium mole fractions: x_H2 = 5.899e-01 x_C2H2 = 1.847e-01 x_C2H4 = 2.204e-01 x_C2H6 = 5.014e-03 (T[K],P[bar])=(1351.0, 1.0) Equilibrium mole fractions: x_H2 = 6.001e-01 x_C2H2 = 2.039e-01 x_C2H4 = 1.923e-01 x_C2H6 = 3.656e-03 (T[K],P[bar])=(1371.4, 1.0) Equilibrium mole fractions: x_H2 = 6.095e-01 x_C2H2 = 2.216e-01 x_C2H4 = 1.662e-01 x_C2H6 = 2.651e-03 (T[K],P[bar])=(1391.8, 1.0) Equilibrium mole fractions: x_H2 = 6.179e-01 x_C2H2 = 2.377e-01 x_C2H4 = 1.425e-01 x_C2H6 = 1.914e-03 (T[K],P[bar])=(1412.2, 1.0) Equilibrium mole fractions: x_H2 = 6.253e-01 x_C2H2 = 2.519e-01 x_C2H4 = 1.214e-01 x_C2H6 = 1.378e-03 (T[K],P[bar])=(1432.7, 1.0) Equilibrium mole fractions: x_H2 = 6.317e-01 x_C2H2 = 2.644e-01 x_C2H4 = 1.029e-01 x_C2H6 = 9.908e-04 (T[K],P[bar])=(1453.1, 1.0) Equilibrium mole fractions: x_H2 = 6.372e-01 x_C2H2 = 2.752e-01 x_C2H4 = 8.691e-02 x_C2H6 = 7.120e-04 (T[K],P[bar])=(1473.5, 1.0) Equilibrium mole fractions: x_H2 = 6.419e-01 x_C2H2 = 2.843e-01 x_C2H4 = 7.322e-02 x_C2H6 = 5.122e-04 (T[K],P[bar])=(1493.9, 1.0) Equilibrium mole fractions: x_H2 = 6.459e-01 x_C2H2 = 2.921e-01 x_C2H4 = 6.161e-02 x_C2H6 = 3.693e-04 (T[K],P[bar])=(1514.3, 1.0) Equilibrium mole fractions: x_H2 = 6.492e-01 x_C2H2 = 2.987e-01 x_C2H4 = 5.182e-02 x_C2H6 = 2.670e-04 (T[K],P[bar])=(1534.7, 1.0) Equilibrium mole fractions: x_H2 = 6.520e-01 x_C2H2 = 3.042e-01 x_C2H4 = 4.361e-02 x_C2H6 = 1.938e-04 (T[K],P[bar])=(1555.1, 1.0) Equilibrium mole fractions: x_H2 = 6.543e-01 x_C2H2 = 3.088e-01 x_C2H4 = 3.674e-02 x_C2H6 = 1.413e-04 (T[K],P[bar])=(1575.5, 1.0) Equilibrium mole fractions: x_H2 = 6.563e-01 x_C2H2 = 3.126e-01 x_C2H4 = 3.100e-02 x_C2H6 = 1.035e-04 (T[K],P[bar])=(1595.9, 1.0) Equilibrium mole fractions: x_H2 = 6.579e-01 x_C2H2 = 3.158e-01 x_C2H4 = 2.622e-02 x_C2H6 = 7.619e-05 (T[K],P[bar])=(1616.3, 1.0) Equilibrium mole fractions: x_H2 = 6.592e-01 x_C2H2 = 3.185e-01 x_C2H4 = 2.222e-02 x_C2H6 = 5.640e-05 (T[K],P[bar])=(1636.7, 1.0) Equilibrium mole fractions: x_H2 = 6.603e-01 x_C2H2 = 3.207e-01 x_C2H4 = 1.888e-02 x_C2H6 = 4.199e-05 (T[K],P[bar])=(1657.1, 1.0) Equilibrium mole fractions: x_H2 = 6.613e-01 x_C2H2 = 3.226e-01 x_C2H4 = 1.608e-02 x_C2H6 = 3.144e-05 (T[K],P[bar])=(1677.6, 1.0) Equilibrium mole fractions: x_H2 = 6.621e-01 x_C2H2 = 3.242e-01 x_C2H4 = 1.374e-02 x_C2H6 = 2.367e-05 (T[K],P[bar])=(1698.0, 1.0) Equilibrium mole fractions: x_H2 = 6.627e-01 x_C2H2 = 3.255e-01 x_C2H4 = 1.177e-02 x_C2H6 = 1.792e-05 (T[K],P[bar])=(1718.4, 1.0) Equilibrium mole fractions: x_H2 = 6.633e-01 x_C2H2 = 3.266e-01 x_C2H4 = 1.012e-02 x_C2H6 = 1.365e-05 (T[K],P[bar])=(1738.8, 1.0) Equilibrium mole fractions: x_H2 = 6.638e-01 x_C2H2 = 3.275e-01 x_C2H4 = 8.718e-03 x_C2H6 = 1.045e-05 (T[K],P[bar])=(1759.2, 1.0) Equilibrium mole fractions: x_H2 = 6.641e-01 x_C2H2 = 3.283e-01 x_C2H4 = 7.534e-03 x_C2H6 = 8.044e-06 (T[K],P[bar])=(1779.6, 1.0) Equilibrium mole fractions: x_H2 = 6.645e-01 x_C2H2 = 3.290e-01 x_C2H4 = 6.530e-03 x_C2H6 = 6.225e-06 (T[K],P[bar])=(1800.0, 1.0) Equilibrium mole fractions: x_H2 = 6.648e-01 x_C2H2 = 3.295e-01 x_C2H4 = 5.675e-03 x_C2H6 = 4.842e-06 (T[K],P[bar])=(800.0, 2.0) Equilibrium mole fractions: x_H2 = 4.553e-02 x_C2H2 = 3.925e-06 x_C2H4 = 4.552e-02 x_C2H6 = 9.089e-01 (T[K],P[bar])=(820.4, 2.0) Equilibrium mole fractions: x_H2 = 5.869e-02 x_C2H2 = 7.799e-06 x_C2H4 = 5.868e-02 x_C2H6 = 8.826e-01 (T[K],P[bar])=(840.8, 2.0) Equilibrium mole fractions: x_H2 = 7.444e-02 x_C2H2 = 1.499e-05 x_C2H4 = 7.441e-02 x_C2H6 = 8.511e-01 (T[K],P[bar])=(861.2, 2.0) Equilibrium mole fractions: x_H2 = 9.291e-02 x_C2H2 = 2.795e-05 x_C2H4 = 9.286e-02 x_C2H6 = 8.142e-01 (T[K],P[bar])=(881.6, 2.0) Equilibrium mole fractions: x_H2 = 1.142e-01 x_C2H2 = 5.065e-05 x_C2H4 = 1.141e-01 x_C2H6 = 7.717e-01 (T[K],P[bar])=(902.0, 2.0) Equilibrium mole fractions: x_H2 = 1.381e-01 x_C2H2 = 8.937e-05 x_C2H4 = 1.379e-01 x_C2H6 = 7.239e-01 (T[K],P[bar])=(922.4, 2.0) Equilibrium mole fractions: x_H2 = 1.645e-01 x_C2H2 = 1.538e-04 x_C2H4 = 1.642e-01 x_C2H6 = 6.711e-01 (T[K],P[bar])=(942.9, 2.0) Equilibrium mole fractions: x_H2 = 1.930e-01 x_C2H2 = 2.584e-04 x_C2H4 = 1.925e-01 x_C2H6 = 6.142e-01 (T[K],P[bar])=(963.3, 2.0) Equilibrium mole fractions: x_H2 = 2.230e-01 x_C2H2 = 4.248e-04 x_C2H4 = 2.221e-01 x_C2H6 = 5.545e-01 (T[K],P[bar])=(983.7, 2.0) Equilibrium mole fractions: x_H2 = 2.537e-01 x_C2H2 = 6.839e-04 x_C2H4 = 2.523e-01 x_C2H6 = 4.933e-01 (T[K],P[bar])=(1004.1, 2.0) Equilibrium mole fractions: x_H2 = 2.845e-01 x_C2H2 = 1.079e-03 x_C2H4 = 2.823e-01 x_C2H6 = 4.322e-01 (T[K],P[bar])=(1024.5, 2.0) Equilibrium mole fractions: x_H2 = 3.144e-01 x_C2H2 = 1.671e-03 x_C2H4 = 3.110e-01 x_C2H6 = 3.729e-01 (T[K],P[bar])=(1044.9, 2.0) Equilibrium mole fractions: x_H2 = 3.428e-01 x_C2H2 = 2.540e-03 x_C2H4 = 3.377e-01 x_C2H6 = 3.170e-01 (T[K],P[bar])=(1065.3, 2.0) Equilibrium mole fractions: x_H2 = 3.691e-01 x_C2H2 = 3.794e-03 x_C2H4 = 3.615e-01 x_C2H6 = 2.657e-01 (T[K],P[bar])=(1085.7, 2.0) Equilibrium mole fractions: x_H2 = 3.929e-01 x_C2H2 = 5.569e-03 x_C2H4 = 3.817e-01 x_C2H6 = 2.198e-01 (T[K],P[bar])=(1106.1, 2.0) Equilibrium mole fractions: x_H2 = 4.141e-01 x_C2H2 = 8.035e-03 x_C2H4 = 3.980e-01 x_C2H6 = 1.798e-01 (T[K],P[bar])=(1126.5, 2.0) Equilibrium mole fractions: x_H2 = 4.329e-01 x_C2H2 = 1.139e-02 x_C2H4 = 4.101e-01 x_C2H6 = 1.457e-01 (T[K],P[bar])=(1146.9, 2.0) Equilibrium mole fractions: x_H2 = 4.494e-01 x_C2H2 = 1.587e-02 x_C2H4 = 4.177e-01 x_C2H6 = 1.171e-01 (T[K],P[bar])=(1167.3, 2.0) Equilibrium mole fractions: x_H2 = 4.641e-01 x_C2H2 = 2.172e-02 x_C2H4 = 4.207e-01 x_C2H6 = 9.346e-02 (T[K],P[bar])=(1187.8, 2.0) Equilibrium mole fractions: x_H2 = 4.775e-01 x_C2H2 = 2.916e-02 x_C2H4 = 4.192e-01 x_C2H6 = 7.418e-02 (T[K],P[bar])=(1208.2, 2.0) Equilibrium mole fractions: x_H2 = 4.899e-01 x_C2H2 = 3.842e-02 x_C2H4 = 4.131e-01 x_C2H6 = 5.859e-02 (T[K],P[bar])=(1228.6, 2.0) Equilibrium mole fractions: x_H2 = 5.018e-01 x_C2H2 = 4.960e-02 x_C2H4 = 4.026e-01 x_C2H6 = 4.605e-02 (T[K],P[bar])=(1249.0, 2.0) Equilibrium mole fractions: x_H2 = 5.134e-01 x_C2H2 = 6.276e-02 x_C2H4 = 3.878e-01 x_C2H6 = 3.603e-02 (T[K],P[bar])=(1269.4, 2.0) Equilibrium mole fractions: x_H2 = 5.249e-01 x_C2H2 = 7.779e-02 x_C2H4 = 3.693e-01 x_C2H6 = 2.805e-02 (T[K],P[bar])=(1289.8, 2.0) Equilibrium mole fractions: x_H2 = 5.364e-01 x_C2H2 = 9.446e-02 x_C2H4 = 3.475e-01 x_C2H6 = 2.172e-02 (T[K],P[bar])=(1310.2, 2.0) Equilibrium mole fractions: x_H2 = 5.478e-01 x_C2H2 = 1.124e-01 x_C2H4 = 3.230e-01 x_C2H6 = 1.673e-02 (T[K],P[bar])=(1330.6, 2.0) Equilibrium mole fractions: x_H2 = 5.592e-01 x_C2H2 = 1.312e-01 x_C2H4 = 2.968e-01 x_C2H6 = 1.280e-02 (T[K],P[bar])=(1351.0, 2.0) Equilibrium mole fractions: x_H2 = 5.703e-01 x_C2H2 = 1.504e-01 x_C2H4 = 2.696e-01 x_C2H6 = 9.739e-03 (T[K],P[bar])=(1371.4, 2.0) Equilibrium mole fractions: x_H2 = 5.810e-01 x_C2H2 = 1.694e-01 x_C2H4 = 2.422e-01 x_C2H6 = 7.365e-03 (T[K],P[bar])=(1391.8, 2.0) Equilibrium mole fractions: x_H2 = 5.911e-01 x_C2H2 = 1.878e-01 x_C2H4 = 2.155e-01 x_C2H6 = 5.538e-03 (T[K],P[bar])=(1412.2, 2.0) Equilibrium mole fractions: x_H2 = 6.006e-01 x_C2H2 = 2.053e-01 x_C2H4 = 1.900e-01 x_C2H6 = 4.144e-03 (T[K],P[bar])=(1432.7, 2.0) Equilibrium mole fractions: x_H2 = 6.092e-01 x_C2H2 = 2.215e-01 x_C2H4 = 1.663e-01 x_C2H6 = 3.087e-03 (T[K],P[bar])=(1453.1, 2.0) Equilibrium mole fractions: x_H2 = 6.170e-01 x_C2H2 = 2.362e-01 x_C2H4 = 1.445e-01 x_C2H6 = 2.292e-03 (T[K],P[bar])=(1473.5, 2.0) Equilibrium mole fractions: x_H2 = 6.239e-01 x_C2H2 = 2.495e-01 x_C2H4 = 1.249e-01 x_C2H6 = 1.698e-03 (T[K],P[bar])=(1493.9, 2.0) Equilibrium mole fractions: x_H2 = 6.300e-01 x_C2H2 = 2.613e-01 x_C2H4 = 1.075e-01 x_C2H6 = 1.257e-03 (T[K],P[bar])=(1514.3, 2.0) Equilibrium mole fractions: x_H2 = 6.353e-01 x_C2H2 = 2.716e-01 x_C2H4 = 9.221e-02 x_C2H6 = 9.298e-04 (T[K],P[bar])=(1534.7, 2.0) Equilibrium mole fractions: x_H2 = 6.399e-01 x_C2H2 = 2.805e-01 x_C2H4 = 7.893e-02 x_C2H6 = 6.884e-04 (T[K],P[bar])=(1555.1, 2.0) Equilibrium mole fractions: x_H2 = 6.438e-01 x_C2H2 = 2.882e-01 x_C2H4 = 6.748e-02 x_C2H6 = 5.105e-04 (T[K],P[bar])=(1575.5, 2.0) Equilibrium mole fractions: x_H2 = 6.472e-01 x_C2H2 = 2.948e-01 x_C2H4 = 5.766e-02 x_C2H6 = 3.795e-04 (T[K],P[bar])=(1595.9, 2.0) Equilibrium mole fractions: x_H2 = 6.501e-01 x_C2H2 = 3.004e-01 x_C2H4 = 4.928e-02 x_C2H6 = 2.830e-04 (T[K],P[bar])=(1616.3, 2.0) Equilibrium mole fractions: x_H2 = 6.525e-01 x_C2H2 = 3.052e-01 x_C2H4 = 4.214e-02 x_C2H6 = 2.118e-04 (T[K],P[bar])=(1636.7, 2.0) Equilibrium mole fractions: x_H2 = 6.545e-01 x_C2H2 = 3.092e-01 x_C2H4 = 3.608e-02 x_C2H6 = 1.591e-04 (T[K],P[bar])=(1657.1, 2.0) Equilibrium mole fractions: x_H2 = 6.563e-01 x_C2H2 = 3.127e-01 x_C2H4 = 3.094e-02 x_C2H6 = 1.200e-04 (T[K],P[bar])=(1677.6, 2.0) Equilibrium mole fractions: x_H2 = 6.577e-01 x_C2H2 = 3.156e-01 x_C2H4 = 2.658e-02 x_C2H6 = 9.097e-05 (T[K],P[bar])=(1698.0, 2.0) Equilibrium mole fractions: x_H2 = 6.590e-01 x_C2H2 = 3.181e-01 x_C2H4 = 2.288e-02 x_C2H6 = 6.927e-05 (T[K],P[bar])=(1718.4, 2.0) Equilibrium mole fractions: x_H2 = 6.601e-01 x_C2H2 = 3.202e-01 x_C2H4 = 1.974e-02 x_C2H6 = 5.299e-05 (T[K],P[bar])=(1738.8, 2.0) Equilibrium mole fractions: x_H2 = 6.610e-01 x_C2H2 = 3.219e-01 x_C2H4 = 1.707e-02 x_C2H6 = 4.074e-05 (T[K],P[bar])=(1759.2, 2.0) Equilibrium mole fractions: x_H2 = 6.617e-01 x_C2H2 = 3.235e-01 x_C2H4 = 1.479e-02 x_C2H6 = 3.147e-05 (T[K],P[bar])=(1779.6, 2.0) Equilibrium mole fractions: x_H2 = 6.624e-01 x_C2H2 = 3.248e-01 x_C2H4 = 1.285e-02 x_C2H6 = 2.442e-05 (T[K],P[bar])=(1800.0, 2.0) Equilibrium mole fractions: x_H2 = 6.629e-01 x_C2H2 = 3.259e-01 x_C2H4 = 1.119e-02 x_C2H6 = 1.905e-05 (T[K],P[bar])=(800.0, 3.0) Equilibrium mole fractions: x_H2 = 3.750e-02 x_C2H2 = 2.617e-06 x_C2H4 = 3.750e-02 x_C2H6 = 9.250e-01 (T[K],P[bar])=(820.4, 3.0) Equilibrium mole fractions: x_H2 = 4.847e-02 x_C2H2 = 5.199e-06 x_C2H4 = 4.846e-02 x_C2H6 = 9.031e-01 (T[K],P[bar])=(840.8, 3.0) Equilibrium mole fractions: x_H2 = 6.168e-02 x_C2H2 = 9.996e-06 x_C2H4 = 6.166e-02 x_C2H6 = 8.767e-01 (T[K],P[bar])=(861.2, 3.0) Equilibrium mole fractions: x_H2 = 7.730e-02 x_C2H2 = 1.864e-05 x_C2H4 = 7.726e-02 x_C2H6 = 8.454e-01 (T[K],P[bar])=(881.6, 3.0) Equilibrium mole fractions: x_H2 = 9.544e-02 x_C2H2 = 3.377e-05 x_C2H4 = 9.537e-02 x_C2H6 = 8.092e-01 (T[K],P[bar])=(902.0, 3.0) Equilibrium mole fractions: x_H2 = 1.161e-01 x_C2H2 = 5.959e-05 x_C2H4 = 1.160e-01 x_C2H6 = 7.678e-01 (T[K],P[bar])=(922.4, 3.0) Equilibrium mole fractions: x_H2 = 1.393e-01 x_C2H2 = 1.026e-04 x_C2H4 = 1.391e-01 x_C2H6 = 7.216e-01 (T[K],P[bar])=(942.9, 3.0) Equilibrium mole fractions: x_H2 = 1.646e-01 x_C2H2 = 1.724e-04 x_C2H4 = 1.643e-01 x_C2H6 = 6.709e-01 (T[K],P[bar])=(963.3, 3.0) Equilibrium mole fractions: x_H2 = 1.919e-01 x_C2H2 = 2.835e-04 x_C2H4 = 1.913e-01 x_C2H6 = 6.165e-01 (T[K],P[bar])=(983.7, 3.0) Equilibrium mole fractions: x_H2 = 2.205e-01 x_C2H2 = 4.565e-04 x_C2H4 = 2.196e-01 x_C2H6 = 5.595e-01 (T[K],P[bar])=(1004.1, 3.0) Equilibrium mole fractions: x_H2 = 2.498e-01 x_C2H2 = 7.208e-04 x_C2H4 = 2.484e-01 x_C2H6 = 5.010e-01 (T[K],P[bar])=(1024.5, 3.0) Equilibrium mole fractions: x_H2 = 2.793e-01 x_C2H2 = 1.117e-03 x_C2H4 = 2.770e-01 x_C2H6 = 4.426e-01 (T[K],P[bar])=(1044.9, 3.0) Equilibrium mole fractions: x_H2 = 3.081e-01 x_C2H2 = 1.700e-03 x_C2H4 = 3.047e-01 x_C2H6 = 3.856e-01 (T[K],P[bar])=(1065.3, 3.0) Equilibrium mole fractions: x_H2 = 3.356e-01 x_C2H2 = 2.543e-03 x_C2H4 = 3.305e-01 x_C2H6 = 3.313e-01 (T[K],P[bar])=(1085.7, 3.0) Equilibrium mole fractions: x_H2 = 3.613e-01 x_C2H2 = 3.742e-03 x_C2H4 = 3.538e-01 x_C2H6 = 2.811e-01 (T[K],P[bar])=(1106.1, 3.0) Equilibrium mole fractions: x_H2 = 3.849e-01 x_C2H2 = 5.416e-03 x_C2H4 = 3.741e-01 x_C2H6 = 2.356e-01 (T[K],P[bar])=(1126.5, 3.0) Equilibrium mole fractions: x_H2 = 4.062e-01 x_C2H2 = 7.713e-03 x_C2H4 = 3.907e-01 x_C2H6 = 1.954e-01 (T[K],P[bar])=(1146.9, 3.0) Equilibrium mole fractions: x_H2 = 4.251e-01 x_C2H2 = 1.081e-02 x_C2H4 = 4.035e-01 x_C2H6 = 1.605e-01 (T[K],P[bar])=(1167.3, 3.0) Equilibrium mole fractions: x_H2 = 4.420e-01 x_C2H2 = 1.490e-02 x_C2H4 = 4.122e-01 x_C2H6 = 1.308e-01 (T[K],P[bar])=(1187.8, 3.0) Equilibrium mole fractions: x_H2 = 4.571e-01 x_C2H2 = 2.019e-02 x_C2H4 = 4.168e-01 x_C2H6 = 1.059e-01 (T[K],P[bar])=(1208.2, 3.0) Equilibrium mole fractions: x_H2 = 4.708e-01 x_C2H2 = 2.690e-02 x_C2H4 = 4.170e-01 x_C2H6 = 8.526e-02 (T[K],P[bar])=(1228.6, 3.0) Equilibrium mole fractions: x_H2 = 4.835e-01 x_C2H2 = 3.521e-02 x_C2H4 = 4.130e-01 x_C2H6 = 6.829e-02 (T[K],P[bar])=(1249.0, 3.0) Equilibrium mole fractions: x_H2 = 4.954e-01 x_C2H2 = 4.526e-02 x_C2H4 = 4.049e-01 x_C2H6 = 5.445e-02 (T[K],P[bar])=(1269.4, 3.0) Equilibrium mole fractions: x_H2 = 5.069e-01 x_C2H2 = 5.710e-02 x_C2H4 = 3.927e-01 x_C2H6 = 4.322e-02 (T[K],P[bar])=(1289.8, 3.0) Equilibrium mole fractions: x_H2 = 5.183e-01 x_C2H2 = 7.069e-02 x_C2H4 = 3.769e-01 x_C2H6 = 3.415e-02 (T[K],P[bar])=(1310.2, 3.0) Equilibrium mole fractions: x_H2 = 5.295e-01 x_C2H2 = 8.588e-02 x_C2H4 = 3.578e-01 x_C2H6 = 2.686e-02 (T[K],P[bar])=(1330.6, 3.0) Equilibrium mole fractions: x_H2 = 5.407e-01 x_C2H2 = 1.024e-01 x_C2H4 = 3.359e-01 x_C2H6 = 2.101e-02 (T[K],P[bar])=(1351.0, 3.0) Equilibrium mole fractions: x_H2 = 5.518e-01 x_C2H2 = 1.199e-01 x_C2H4 = 3.120e-01 x_C2H6 = 1.636e-02 (T[K],P[bar])=(1371.4, 3.0) Equilibrium mole fractions: x_H2 = 5.627e-01 x_C2H2 = 1.380e-01 x_C2H4 = 2.867e-01 x_C2H6 = 1.266e-02 (T[K],P[bar])=(1391.8, 3.0) Equilibrium mole fractions: x_H2 = 5.732e-01 x_C2H2 = 1.562e-01 x_C2H4 = 2.608e-01 x_C2H6 = 9.747e-03 (T[K],P[bar])=(1412.2, 3.0) Equilibrium mole fractions: x_H2 = 5.834e-01 x_C2H2 = 1.742e-01 x_C2H4 = 2.350e-01 x_C2H6 = 7.465e-03 (T[K],P[bar])=(1432.7, 3.0) Equilibrium mole fractions: x_H2 = 5.929e-01 x_C2H2 = 1.915e-01 x_C2H4 = 2.099e-01 x_C2H6 = 5.690e-03 (T[K],P[bar])=(1453.1, 3.0) Equilibrium mole fractions: x_H2 = 6.018e-01 x_C2H2 = 2.079e-01 x_C2H4 = 1.860e-01 x_C2H6 = 4.318e-03 (T[K],P[bar])=(1473.5, 3.0) Equilibrium mole fractions: x_H2 = 6.099e-01 x_C2H2 = 2.231e-01 x_C2H4 = 1.637e-01 x_C2H6 = 3.265e-03 (T[K],P[bar])=(1493.9, 3.0) Equilibrium mole fractions: x_H2 = 6.173e-01 x_C2H2 = 2.370e-01 x_C2H4 = 1.433e-01 x_C2H6 = 2.462e-03 (T[K],P[bar])=(1514.3, 3.0) Equilibrium mole fractions: x_H2 = 6.238e-01 x_C2H2 = 2.495e-01 x_C2H4 = 1.248e-01 x_C2H6 = 1.853e-03 (T[K],P[bar])=(1534.7, 3.0) Equilibrium mole fractions: x_H2 = 6.296e-01 x_C2H2 = 2.607e-01 x_C2H4 = 1.083e-01 x_C2H6 = 1.394e-03 (T[K],P[bar])=(1555.1, 3.0) Equilibrium mole fractions: x_H2 = 6.347e-01 x_C2H2 = 2.705e-01 x_C2H4 = 9.368e-02 x_C2H6 = 1.048e-03 (T[K],P[bar])=(1575.5, 3.0) Equilibrium mole fractions: x_H2 = 6.392e-01 x_C2H2 = 2.791e-01 x_C2H4 = 8.089e-02 x_C2H6 = 7.888e-04 (T[K],P[bar])=(1595.9, 3.0) Equilibrium mole fractions: x_H2 = 6.430e-01 x_C2H2 = 2.866e-01 x_C2H4 = 6.976e-02 x_C2H6 = 5.945e-04 (T[K],P[bar])=(1616.3, 3.0) Equilibrium mole fractions: x_H2 = 6.463e-01 x_C2H2 = 2.931e-01 x_C2H4 = 6.014e-02 x_C2H6 = 4.490e-04 (T[K],P[bar])=(1636.7, 3.0) Equilibrium mole fractions: x_H2 = 6.492e-01 x_C2H2 = 2.987e-01 x_C2H4 = 5.185e-02 x_C2H6 = 3.401e-04 (T[K],P[bar])=(1657.1, 3.0) Equilibrium mole fractions: x_H2 = 6.516e-01 x_C2H2 = 3.034e-01 x_C2H4 = 4.472e-02 x_C2H6 = 2.584e-04 (T[K],P[bar])=(1677.6, 3.0) Equilibrium mole fractions: x_H2 = 6.537e-01 x_C2H2 = 3.075e-01 x_C2H4 = 3.861e-02 x_C2H6 = 1.970e-04 (T[K],P[bar])=(1698.0, 3.0) Equilibrium mole fractions: x_H2 = 6.554e-01 x_C2H2 = 3.110e-01 x_C2H4 = 3.338e-02 x_C2H6 = 1.508e-04 (T[K],P[bar])=(1718.4, 3.0) Equilibrium mole fractions: x_H2 = 6.570e-01 x_C2H2 = 3.140e-01 x_C2H4 = 2.890e-02 x_C2H6 = 1.159e-04 (T[K],P[bar])=(1738.8, 3.0) Equilibrium mole fractions: x_H2 = 6.582e-01 x_C2H2 = 3.166e-01 x_C2H4 = 2.507e-02 x_C2H6 = 8.940e-05 (T[K],P[bar])=(1759.2, 3.0) Equilibrium mole fractions: x_H2 = 6.594e-01 x_C2H2 = 3.188e-01 x_C2H4 = 2.179e-02 x_C2H6 = 6.928e-05 (T[K],P[bar])=(1779.6, 3.0) Equilibrium mole fractions: x_H2 = 6.603e-01 x_C2H2 = 3.207e-01 x_C2H4 = 1.897e-02 x_C2H6 = 5.392e-05 (T[K],P[bar])=(1800.0, 3.0) Equilibrium mole fractions: x_H2 = 6.611e-01 x_C2H2 = 3.223e-01 x_C2H4 = 1.656e-02 x_C2H6 = 4.215e-05
import matplotlib.pyplot as plt
plt.figure(1, figsize=(8, 6))
from chen_3170.help import color_map
colors = color_map(len(pressure_pts))
line_style = ['-','--','-.',':'] # for species
for p_pt in range(len(pressure_pts)):
pressure = pressure_pts[p_pt]
color = colors[p_pt]
for i in range(len(species)):
plt.plot(temp_pts, x_vec_block[p_pt,:,i],line_style[i],color=color,label=r' P=%1.1f; $x_{%s}$'%(pressure,species[i]))
plt.xlabel(r'$T$ [K]',fontsize=16)
plt.ylabel(r'$x_j$',fontsize=16)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.legend(loc='center right',bbox_to_anchor=(1.4,0.5),fontsize=16)
plt.title('Temperature-Pressure Dependency',fontsize=20)
plt.grid(True)
plt.show()
print('')
Observations: the sweet spot for simultaneous production of ethylene and acetylene is at the cross over point. The cross over point at the lowest pressure occur at lower temperature, therefore lower pressure is more economical.
import matplotlib.pyplot as plt
(fig, ax1) = plt.subplots(1, figsize=(7, 7))
ax2 = ax1.twinx()
line_style = ['-','--','-.'] # pressure
for p_pt in range(len(pressure_pts)):
pressure = pressure_pts[p_pt]
ax1.plot(temp_pts,k_x_cte_block[p_pt,:,0],line_style[p_pt],color='blue',label=r'$P=%1.1f$'%pressure )
ax2.plot(temp_pts,k_x_cte_block[p_pt,:,1],line_style[p_pt],color='red',label=r'$P=%1.1f$'%pressure )
ax1.legend(loc='upper left',fontsize=12)
ax1.set_xlabel(r'$T$ [K]',fontsize=16)
ax1.set_ylabel(r'$K_{x,%i}(T,P)$'%reaction_ids[0],fontsize=16,color='blue')
ax1.tick_params(axis='y', labelcolor='blue', labelsize=14)
ax1.tick_params(axis='x', labelsize=14)
ax1.grid(True)
ax2.set_ylabel(r'$K_{x,%i}(T,P)$'%reaction_ids[1],fontsize=16,color='red')
ax2.tick_params(axis='y', labelcolor='red', labelsize=14)
ax2.tick_params(axis='x', labelsize=14)
ax2.legend(loc='lower left',fontsize=12)
#ax2.grid(True)
fig.tight_layout() # otherwise the right y-label is slightly clipped
#plt.grid(True)
plt.show()
print('')
import matplotlib.pyplot as plt
(fig, ax1) = plt.subplots(1, figsize=(7, 7))
ax2 = ax1.twinx()
ax1.plot(temp_pts,delta_rxn_h_std_state_mtrx[:,0],'b-',label=reactions[0] )
ax2.plot(temp_pts,delta_rxn_h_std_state_mtrx[:,1],'r-',label=reactions[1] )
ax1.legend(loc='upper left',fontsize=12)
ax1.set_xlabel(r'$T$ [K]',fontsize=16)
ax1.set_ylabel(r'$\Delta_{rxn,%i}H^\circ(T)$ [cal/mol]'%reaction_ids[0],fontsize=16,color='blue')
ax1.tick_params(axis='y', labelcolor='blue', labelsize=14)
ax1.tick_params(axis='x', labelsize=14)
ax1.grid(True)
ax2.set_ylabel(r'$\Delta_{rxn,%i}H^\circ(T)$ [cal/mol]'%reaction_ids[1],fontsize=16,color='red')
ax2.tick_params(axis='y', labelcolor='red', labelsize=14)
ax2.tick_params(axis='x', labelsize=14)
ax2.legend(loc='lower left',fontsize=12)
#ax2.grid(True)
fig.tight_layout() # otherwise the right y-label is slightly clipped
#plt.grid(True)
plt.show()
print('')