Material Balances XV
Contents
Material Balances XV#
This lecture continues with a discussion of yield, selectivity, and optimization.
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
from scipy.integrate import solve_ivp
from scipy.interpolate import interp1d
Example Problem 01#
This is a slightly modified version of Example 8.4 from Fogler
The following pair of liquid-phase reactions occur in series. Each reaction is first order in its reactant.
By “series” or “sequential” reactions we mean that a product of the first reaction is a reactant in a second reaction, and so on. These are very common systems in Chemical Reaction Engineering – especially in the petrochemical industry, where we often make reactive things from inert alkanes. Quite often, the reactive things we make then react away to form things we don’t want. Selectivity control in series reactions requires careful attention to operating conditions and reaction times (space times for flow reactors).
You carry out these reactions in a CSTR operating at steady state. You may assume the density of the fluid does not change with reaction and so the inlet volumetric flowrate is equal to the outlet volumetric flowrate. You have the following data available:
Find the space time required to maximize the yield of B (with respect to A) in this CSTR.
Solution to Example Problem 01#
To get yields, selectivities, and conversions, we need to write material balances on all species. This allows us to determine the exit molar flowrates of all species. This problem is a bit unique in that we can solve balances for A, B, and C sequentially. This is enabled by the fact that the balance on A is independent of species B; solving the balance on A enables us to solve the balance on B; the balance on B is independent of species C; and solving the balance on B enables us to solve the balance on C. See below.
Material Balance on A:#
We aren’t given a ton of information except for the feed concentration of A. Presumably, the feed concentrations of B and C are both zero. We are also told we want to work with a space time, which is defined for a CSTR as:
I’m going to convert my balance so that it is written in terms of concentration and space time instead of molar flowrate and volume. I can express molar flowrates in terms of concentrations and volumetric flowrates:
We know the intensive production rate of A:
Since the reactions are first order in the reactant:
Making these substitutions, we find:
We are told that we can assume constant density in this system, which means (at steady state)
Thus:
Dividing by \(Q_f\), we express the balance in terms of a space time and concentrations.
And we can solve this for \(C_A\), the exit concentration of A from the CSTR:
This is an important result. It gives us the exit concentration of A as a function of space time. Now we move on to a balance on B:
Balance on B#
We can replace molar flowrates with the product of concentration and volumetric flowrate, noting again that \(Q = Q_f\) for this problem:
We can define the intensive production rate of B, \(R_B\):
We know the reactions are both first order in their reactant, so:
And we can substitute this into the material balance:
Dividing through by \(Q_f\):
Noting that \(C_{B,f} = 0\):
And we can solve this for \(C_B\):
The above is expressed in terms of \(C_A\). We can substitute that expression in to get:
This is another important equation as it gives us \(C_B\) as a function of \(\tau\).
Balance on C#
Technically, we don’t need this to solve for yield/selectivity to B, but I’m curious to see how it behaves, so we’ll go ahead and write a balance on C:
Converting to a function of concentration:
Defining the production rate for C:
We get:
Dividing by \(Q_f\):
Noting that \(C_f = 0\):
Which we can solve for \(C_C\):
And we can substitute the expression for \(C_B\) here to get:
There we have our concentrations expressed as 3 functions of tau. These are developed as lambda functions below, which are plotted as a function of tau.
k1 = 0.5 #1/h
k2 = 0.2 #1/h
CAf = 20 #mol/L
CA = lambda tau: CAf/(1 + k1*tau)
CB = lambda tau: k1*CAf*tau/(1+k1*tau)/(1+k2*tau)
CC = lambda tau: k1*k2*CAf*tau**2/(1+k1*tau)/(1+k2*tau)
tauset = np.linspace(0, 30, 200)
plt.plot(tauset, CA(tauset), label = 'CA')
plt.plot(tauset, CB(tauset), label = 'CB')
plt.plot(tauset, CC(tauset), label = 'CC')
plt.title('Species concentration vs. CSTR space time', fontsize = 14)
plt.xlabel('tau (h)', fontsize = 14)
plt.ylabel('Concentration (mol/L)', fontsize = 14)
plt.xlim(0, max(tauset))
plt.ylim(0, 20)
plt.legend()
plt.show()
Defining Conversion, Selectivities, and Yields#
With the concentrations of each species defined, we can now consider Conversion, Yields and Selectivities.
Conversion#
We define the fractional conversion of A as:
Although we have no way to solve for molar flowrates (since Q is unknown), we can express them in terms of concentrations, noting symbolically that \(Q = Q_f\).
We can factor and cancel the volumetric flowrate, giving us the result below, which is always true for a constant density fluid:
Yields#
For a flow reactor, we would generally define the yield of B with respect to A as:
Which simplfies for this problem to:
We actually don’t know the molar flowrates since we don’t know the size (Q, V) of this system, but we can express these in terms of concentrations, noting that since density is constant \(Q = Q_f\):
This reduces to:
When we go to calculate a yield of C with respect to A, we have to note that by summing the two reactions, we can get the overall stoichiometry of converting A into C:
With that, we can use a similar approach to that above to define the yield to C in terms of concentrations:
Selectivities#
We can define the selectivity of species B with respect to A:
Substituting stoichiometric coefficients, expressing in terms of concentrations, and noting that \(Q = Q_f\):
Factoring and cancelling the volumetric flowrate (which we can do in cases where it is constant):
You can follow similar steps to define a selectivity to species C with respect to the reactant A:
We’ll write lambda functions for each of these; note that I’m re-using my concentration functions defined above.
YB = lambda tau: CB(tau)/CAf
YC = lambda tau: CC(tau)/CAf
SB = lambda tau: CB(tau)/(CAf - CA(tau))
SC = lambda tau: CC(tau)/(CAf - CA(tau))
XA = lambda tau: (CAf - CA(tau))/CAf
Now we can plot everything for a graphical analysis; note, this will throw a warning because we divde by zero at some point. This occurs at 0 tau, where \(C_A = C_{A,f}\), causing a zero in the denominator of our selectivity expressions. In this case, Python just reports the error and moves on, but just know that there is an infinity or a NAN value in our selectivity data. In some cases, this may cause us problems, so we should know where the error is coming from and how we can deal with it.
#This just prints out the first few selectivities to show the NAN ("not a number") error
#print(np.array([SB(tauset)[0:5], SC(tauset)[0:5]]).T)
plt.figure(1)
plt.plot(tauset, XA(tauset), label = 'Conversion of A')
plt.plot(tauset, SB(tauset), label = 'Selectivity to B')
plt.plot(tauset, SC(tauset), label = 'Selectivity to C')
plt.xlabel('Space time (h)', fontsize = 14)
plt.ylabel('Conversion/Selectivity', fontsize = 14)
plt.xlim(0,30)
plt.xticks(fontsize = 11)
plt.ylim(0,1)
plt.yticks(fontsize = 11)
plt.legend()
plt.show(1)
plt.figure(2)
plt.plot(tauset, XA(tauset), label = 'Conversion of A')
plt.plot(tauset, YB(tauset), label = 'Yield to B')
plt.plot(tauset, YC(tauset), label = 'Yield to C')
plt.xlabel('Space time (h)', fontsize = 14)
plt.ylabel('Conversion/Yield', fontsize = 14)
plt.xlim(0,30)
plt.xticks(fontsize = 11)
plt.ylim(0,1)
plt.yticks(fontsize = 11)
plt.legend()
plt.show(2)
C:\Users\jqbond\AppData\Local\Temp\ipykernel_15788\3616170257.py:3: RuntimeWarning: invalid value encountered in true_divide
SB = lambda tau: CB(tau)/(CAf - CA(tau))
C:\Users\jqbond\AppData\Local\Temp\ipykernel_15788\3616170257.py:4: RuntimeWarning: invalid value encountered in true_divide
SC = lambda tau: CC(tau)/(CAf - CA(tau))
Maximizing the Yield of B#
An important question is clearly what space time do we need to use to maximize the yield to species B (the intermediate in the series reaction). From the graph of Yield vs. tau, we can see this occurs at a space time of about 3 hours. To get more precise, we can either solve it analytically by taking a derivative of \(Y_{B/A}\) with respect to tau; after some calculus and algebra, we get:
To find the optimum value of tau, we solve the following expression for \(\tau\):
I did this by hand to find:
Which we can plug into a calculator to get:
Or, you know, we could just create an objective function (return negative yield of B at a specific volume) and minimize that with an optimization routine. See below. You get the same answer.
obj = lambda tau: -1*YB(tau) #this objective function returns -1*YB for a given value of tau
optsol = opt.minimize_scalar(obj) #this varies tau until we hit minimum value of -1*YB
print(optsol, '\n')
print(f'The maximum yield to species B is {-1*optsol.fun:3.2f} which occurs at a space time of {optsol.x:3.2f} hours')
fun: -0.3752470442573563
message: '\nOptimization terminated successfully;\nThe returned value satisfies the termination criteria\n(using xtol = 1.48e-08 )'
nfev: 17
nit: 11
success: True
x: 3.162277615990361
The maximum yield to species B is 0.38 which occurs at a space time of 3.16 hours
Example Problem 02#
This is a slightly modified version of Example 8.5 from Fogler
The gas-phase reactions given below are carried out in an isobaric (constant pressure) Packed Bed Reactor
Both reactions follow elementary rate laws (though neither reaction is truly an elementary reaction).
We have the following data available:
plot \(X_A\), \(Y_{C/A}\), \(Y_{D/A}\), \(S_{C/A}\), and \(S_{D/A}\) as a function of catalyst mass (up to 1000kg of catalyst).
what catalyst mass would you recommend in order to maximize the selectivity toward C in this reactor?
what catalyst mass would you recommend in order to maximize the yield of C in this reactor?
Solution to Example Problem 02#
Start by writing material balances on each species. Packed bed balances are almost identical to PFR balances, except we use mass of catalyst instead of volume as the independent variable (the one that defines the size of the reactor), and we write our intensive rate law in units of moles per mass of catalyst per time. Hence:
Where the prime (\(^\prime\)) superscript just means that intensive rates are specified per unit mass of catalyst. Otherwise, intensive production rates are defined as usual.
These reactions have elementary kinetics, so:
We define concentrations in terms of flowrates (dependent variables) as usual for a tubular reactor:
Since this is a gas-phase reaction at low pressure, we can define Q using the ideal gas law:
And since we know the reactor is isobaric, we can calculate the total pressure from the total feed concentration of species:
At this point, everything on the right-hand side of our ODEs is defined in terms of our dependent variables, so we’re good to solve with solve_ivp.
We do want to define conversion, selectivity, and yield in terms of molar flowrates since our problem statement asks us to plot these quantities.
Conversion#
We define conversion as usual:
Selectivities#
Based on the stoichiometry of the first reaction, we can see that 1 mole of C requires 1 mole of A. The overall selectivity to C would therefore be defined as:
The stoichiometric relationship between species A and species D is less clear from the two reactions given; however, if we sum them up:
With that, we can see that producing 1 mole of D will consume 5 moles of A, so:
Yields#
We use the same insights about stoichiometry discussed above in order to define yields to C and D:
def P02(W, var):
FA, FB, FC, FD = var
FT = FA + FB + FC + FD
k1 = 100 #L^3/mol^2/kg/min
k2 = 500 #L^5/mol^4/kg/min
CAf = 0.2 #mol/L
T = 573 #K
R = 0.08206 #L*atm/mol/K
P = CAf*R*T
Q = FT*R*T/P
CA = FA/Q
CB = FB/Q
CC = FC/Q
r1 = k1*CA*CB**2
r2 = k2*CA**2*CC**3
RA = -r1 -2*r2
RB = -2*r1
RC = r1 - 3*r2
RD = r2
return [RA, RB, RC, RD]
FAf = 10 #mol/min
FBf = 10 #mol/min
FCf = 0 #mol/min
FDf = 0 #mol/min
var0 = [FAf, FBf, FCf, FDf]
Wspan = (0, 1000) #kg catalyst
answer = solve_ivp(P02, Wspan, var0, atol = 1e-3, rtol = 1e-3)
W = answer.t
FA = answer.y[0]
FB = answer.y[1]
FC = answer.y[2]
FD = answer.y[3]
XA = (FAf - FA)/FAf
SCA = FC/(FAf - FA)
SDA = 5*FD/(FAf - FA)
YCA = FC/FAf
YDA = 5*FD/FAf
plt.figure(1)
plt.plot(W, XA, label = 'Conversion of A')
plt.xlim(0, 1000)
plt.ylim(0, 1)
plt.xlabel('Mass of Catalyst (kg)', fontsize = 14)
plt.xticks(fontsize = 11)
plt.ylabel('Conversion', fontsize = 14)
plt.yticks(fontsize = 11)
plt.legend()
plt.figure(2)
plt.plot(W, SCA, label = 'Selectivity to C')
plt.plot(W, SDA, label = 'Selectivity to D')
plt.xlim(0, 1000)
plt.ylim(0, 1)
plt.xlabel('Mass of Catalyst (kg)', fontsize = 14)
plt.xticks(fontsize = 11)
plt.ylabel('Selectivity', fontsize = 14)
plt.yticks(fontsize = 11)
plt.legend()
plt.figure(3)
plt.scatter(W, YCA, label = 'Yield to C')
plt.plot(W, YDA, label = 'Yield to D')
plt.xlim(0, 1000)
plt.ylim(0, 1)
plt.xlabel('Mass of Catalyst (kg)', fontsize = 14)
plt.xticks(fontsize = 11)
plt.ylabel('Yield', fontsize = 14)
plt.yticks(fontsize = 11)
plt.legend()
plt.show()
C:\Users\jqbond\AppData\Local\Temp\ipykernel_15788\3189117652.py:20: RuntimeWarning: invalid value encountered in true_divide
SCA = FC/(FAf - FA)
C:\Users\jqbond\AppData\Local\Temp\ipykernel_15788\3189117652.py:21: RuntimeWarning: invalid value encountered in true_divide
SDA = 5*FD/(FAf - FA)
Optimizing Yields and Selectivities#
Here, we have the classic problem in maximizing selectivity/yield to the intermediate species in a sequence of reactions. Selectivity is maximizing by operating as close to zero conversion as possible. This is because, in this limit, the concentration of species C is near zero, so the rate of the second reaction is near zero. Unfortunately, since yield is the product of conversion and selectivity, the yield of C also approaches zero at zero conversion. To produce a viable quantity of product C, you’ll have to tolerate some loss of selectivity – The point of maximum yield will be the mass of catalyst where conversion has increased enough to produce an appreciable quantity of C, but not so much that we are producing mostly D. As shown below, this occurs at 140 kg of catalyst, where we find a maximum yield of B of approximately 40%.
YIELD = interp1d(W, YCA) #This gives me yield as a continuous function of W
#OBJ = lambda W: YIELD(W)*-1 #I want to minimize negative yield.
def OBJ(W):
if 0 <= W <= 1000:
print([W, YIELD(W)])
else:
print(f'Error: {W}kg of catalyst is outside interpolation range')
return -1*YIELD(W)
answer2 = opt.minimize_scalar(OBJ)
W_opt = answer2.x
print(f'The maximum yield of C is {YIELD(W_opt):0.2f} at a catalyst loading of {W_opt:0.0f}kg')
[0.0, array(0.)]
[1.0, array(0.0098198)]
[2.6180339999999998, array(0.02451275)]
[17.892964975441355, array(0.14608862)]
[42.60832264135863, array(0.25468483)]
[54.8926563161009, array(0.30430925)]
[74.76912586917882, array(0.34880455)]
[84.8177098060513, array(0.36144965)]
[101.07666026776482, array(0.38190987)]
Error: 1873.302260594539kg of catalyst is outside interpolation range
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_15788\1450459284.py in <cell line: 10>()
8 return -1*YIELD(W)
9
---> 10 answer2 = opt.minimize_scalar(OBJ)
11 W_opt = answer2.x
12
~\anaconda3\lib\site-packages\scipy\optimize\_minimize.py in minimize_scalar(fun, bracket, bounds, args, method, tol, options)
884 return method(fun, args=args, bracket=bracket, bounds=bounds, **options)
885 elif meth == 'brent':
--> 886 return _minimize_scalar_brent(fun, bracket, args, **options)
887 elif meth == 'bounded':
888 if bounds is None:
~\anaconda3\lib\site-packages\scipy\optimize\_optimize.py in _minimize_scalar_brent(func, brack, args, xtol, maxiter, disp, **unknown_options)
2496 full_output=True, maxiter=maxiter, disp=disp)
2497 brent.set_bracket(brack)
-> 2498 brent.optimize()
2499 x, fval, nit, nfev = brent.get_result(full_output=True)
2500
~\anaconda3\lib\site-packages\scipy\optimize\_optimize.py in optimize(self)
2266 # set up for optimization
2267 func = self.func
-> 2268 xa, xb, xc, fa, fb, fc, funcalls = self.get_bracket_info()
2269 _mintol = self._mintol
2270 _cg = self._cg
~\anaconda3\lib\site-packages\scipy\optimize\_optimize.py in get_bracket_info(self)
2240 ### carefully DOCUMENT any CHANGES in core ##
2241 if brack is None:
-> 2242 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
2243 elif len(brack) == 2:
2244 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],
~\anaconda3\lib\site-packages\scipy\optimize\_optimize.py in bracket(func, xa, xb, args, grow_limit, maxiter)
2790 elif (w - wlim)*(wlim - xc) >= 0.0:
2791 w = wlim
-> 2792 fw = func(*((w,) + args))
2793 funcalls += 1
2794 elif (w - wlim)*(xc - w) > 0.0:
~\AppData\Local\Temp\ipykernel_15788\1450459284.py in OBJ(W)
6 else:
7 print(f'Error: {W}kg of catalyst is outside interpolation range')
----> 8 return -1*YIELD(W)
9
10 answer2 = opt.minimize_scalar(OBJ)
~\anaconda3\lib\site-packages\scipy\interpolate\_polyint.py in __call__(self, x)
76 """
77 x, x_shape = self._prepare_x(x)
---> 78 y = self._evaluate(x)
79 return self._finish_y(y, x_shape)
80
~\anaconda3\lib\site-packages\scipy\interpolate\_interpolate.py in _evaluate(self, x_new)
703 y_new = self._call(self, x_new)
704 if not self._extrapolate:
--> 705 below_bounds, above_bounds = self._check_bounds(x_new)
706 if len(y_new) > 0:
707 # Note fill_value must be broadcast up to the proper size
~\anaconda3\lib\site-packages\scipy\interpolate\_interpolate.py in _check_bounds(self, x_new)
735 "range.")
736 if self.bounds_error and above_bounds.any():
--> 737 raise ValueError("A value in x_new is above the interpolation "
738 "range.")
739
ValueError: A value in x_new is above the interpolation range.