Ind. Eng. Chem. Fundam. 1986, 25, 621-627
621
Algorithmic Synthesis of Chemical Reactor Networks Using Mathematical Programming Luke E. K. Achenle and Lorenz T. Blegler’ Department of Chemical Engineering, Carnegie-Mellon University, Pittsburgh, Pennsylvania 152 13
Given a kinetic mechanism and expressions for the reaction rate, one is frequently concerned with appropriate choices for mixing patterns and heat addition/removal in a reactor design. This study presents a novel nonlinear programming (NLP) formulation for optimally generating this reactor network. Building on earlier ideas for adjoint networks, we also include networks with ideal and nonideal reactors a s well as complex mixing patterns. In addiiion, by controlling heat additionhemoval in the network, this work also extends to nonisothermal reactors. Model and adjoint equations of this formulation form a two-point boundary value problem that interfaces with an efficient optimization strategy. Decisions representing network structure, reactor type, and the amount of heat addition are made through continuous parameters in the model. The method is therefore fairly general and can be applied to large kinetic mechanisms with almost any objective function. Literature examples are presented and solved in order to demonstrate the effectiveness of the approach.
1. Introduction
Over the past 20 years, process synthesis has developed into an active research area. An extensive review has been given by Nishida et al. (1981), where powerful methods for synthesizing heat-exchanger networks and separation sequences are discussed. However, relatively few studies deal with the synthesis of reactor networks. While this is a classical area of research for chemical reactor design, reactor network synthesis also represents a difficult problem because, unlike other synthesis problems, reactors are usually described by differential equation models; heuristics and intuitive methods can only deal with limited cases. Because of the nature of these problems, variational concepts have been successful in describing the sensitivity of reactor networks to decision variables. Horn and Tsai (1967) and Jackson (1968) developed adjoint networks for reactor systems. To account for local mixing, Ravimohan (1971) extended Jackson’s superstructure to include CSTR’s. However, in his formulation, addition of a continuous stirred-tank reactor (CSTR) to the network was a discrete decision and the CSTR volume could not be considered automatically by the optimization algorithm. As discussed later, we avoid this difficulty by representing CSTR’s and plug-flow reactors (PFR’s) as limiting cases of continuous decision variables. While the above studies dealt with general objective functions based on final product rates, Chitra and Govind (1981, 1985) dealt exclusively with reactor networks that maximize yield. After classifying different reaction mechanisms from previous literature studies and applying heuristic strategies for simpler cases, they apply a direct search procedure to optimize a serial network of recycle reactors. Finally, Conti and Paterson (1985) point out that to improve overall “process” yield, reactor selectivity, not reactor yield, should be the criterion for deriving a reactor network. Using heuristics derived from experience, they derive an overall process that is more profitable than that with reactors based on maximum yield. In order to avoid discrete decisions on reactor types, as observed in Ravimohan (1971), Paynter and Haskins (1970) used an axial dispersion reactor model, in which the dimensionless dispersion coefficient, D, was a continuous function of position in the reactor. They concluded that 0196-4313/86/1025-0621$01.50/0
at the optimal solution the control variable u (the derivative of D with respect to position along the reactor) was at a bound; however, this conclusion ignores constraints on D, a state variable. Also, it can be shown that, at least for a zero-order reaction, D need not be a t a bound to satisfy the optimality conditions. In the constant-dispersion model (CDM) discussed here, D is a piecewise constant control variable; the optimization strategy is similar to an approach used by Sargent and Sullivan (1977). It is well-known that in the limit as D approaches zero (infinity) the differential equations of CDM reduce to the PFR (algebraic CSTR equations). The CSTR limit can be proved formally by a perturbation argument; the PFR limit, on the other hand, can be derived through a singular perturbation approach (see Nayfeh, 1973; O’Malley, 1974). Thus, there is justification for using CDM. This paper presents a nonlinear programming formulation for generating complex reactor networks. Following the concepts proposed by Jackson, an optimization problem is developed with split fractions, source points, and sink points as continuous parameters that determine the mixing pattern (i.e., possible bypasses and recycles). With adiabatic and nonisothermal operation in mind, the inlet temperature and reactor jacket temperature (assumed constant) are also included as continuous parameters. In the constant-dispersion model, the reactor type is determined by its corresponding dispersion coefficient. A PFR (CSTR) corresponds to the lower (upper) bound on the dispersion coefficient. Any other value of this coefficient characterizes a reactor with an intermediate degree of mixing. From this framework a system of adjoint equations is developed by using optimal control theory to calculate gradients for a nonlinear programming algorithm. In the limit as the dispersion coefficient goes to zero (the PFR limit), the adjoint system of the isothermal version of the formulation coincides with Jackson’s formulation. Thus, the resulting NLP problem has the ability to handle arbitrary, nonlinear reaction mechanisms and can generate optimal reactor networks based on a wide range of quantifiable objective functions such as yield and selectivity. Optimal reactor networks are generated efficiently by coupling COLSYS (Ascher et al., 1979), a code for solving 0 1986 American Chemical Society
622 0 3.r
,
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986 ./rl
heat e x c h a n g e
--i---A “i
___
t%
il
I
U
t
N
heat exchange
+P-
*Npl,“l
e--+----i n o-i/
4
;\t
P P
Figure 1. Network of PFR’s.
Figure 2. ADR and heat-exchanger combination.
two-point boundary value problems with successive quadratic programming (SQP)to solve the optimization problem. The interface for the combined algorithm allows for synthesis using multiple reaction paths that interact with each other via side streams. Interaction of the reactor network with the rest of the flow sheet is also possible. Moreover, the nonlinear programming framework is flexible enough to accommodate additional side conditions and constraints. The current formulation is presented for nonisothermal homogeneous reactions, which will be specialized to treat isothermal and adiabatic systems. In the next section we define the reactor synthesis problem and describe the Jackson superstructure, as well as the changes we introduce into it. We next point out the assumptions inherent in the mass and energy balances for our basic reactor model. The model will then be coupled with the network superstructure. A summary of the ODE’S that determine the adjoint variables (Lagrange multipliers, corresponding to differential equations that change with position along the reactor) then follows. An optimization algorithm which uses these adjoint variables to select the optimal decision variables in the network is then presented. Finally, some test examples will be discussed to demonstrate the effectiveness of this approach. 2. Problem Definition and General Reactor Network The reactor network synthesis problem can be stated as follows: Given the reaction mechanism, the kinetic expressions, and an objective function (e.g., reactor yield, process yield, selectivity, profit, etc.), what is the optimal reactor network that maximizes this objective? Our objective functions include those that depend directly on concentrations of all species a t the exit of the reactor network. Remember that this form can be made sufficiently general by introducing new state variables and equations that account for heat duties, volumes, etc. In an earlier approach Jackson (1968) postulated a network of reactors and dealt specifically with plug-flow reactor models (PFR). An objective function was optimized with respect to the source or sink positions and the split fractions. In our approach we consider the general nonisothermal reaction with an axial dispersion reactor instead of a PFR, allowing for interstage cooling or heating. We now consider the Jackson superstructure and the changes that we have made to it in more detail. The reactor superstructure given in Figure 1 has a number of smaller structures embedded within it. Our goal is to select the best substructure that optimizes the objective function. In figure 1the segment bounded by (t,,,, ti+,,,) represents the ith reactor of the mth reaction path. The point t,,,, is a source point with a separation matrix uL,,which from physical considerations may be assumed diagonal. Thus, an element ukm of uL,mis the fraction of species k that is diverted from the exit of reactor (1 - 1, rn) to the entrance of reactor (j’, n). Also, a simple splitting of a stream corresponds to u = CUI, where cy is less than one and I is the identity matrix. Other values of u can be used to describe separation operations. Point tl,nis a sink point linked to the source point f bv a side stream. In addition, ~~
t:,m
-.,
//-
t;,r
t,.l,rn
.1+:,rr
Figure 3. Partial representation of the superstructure.
\
\
4\i\‘
f A
PFK
\
A
PFK
> e
Figure 4. (top) Initial and optimum structures for example la. (bottom) Initial and optimum structures for example Ib.
one can think of point tk+l,pas a “degenerate” source point with u = 0; in this respect to,, (an origin) and tN,,,, (a terminus) are examples of degenerate source points. Without loss of generality, we fix the origin a t zero and allow the terminus to vary. This network has one obvious shortcoming in that it does not influence or control mixing on the microscale. Moreover, Jackson’s approach is limited to reactors that are either isothermal or adiabatic or those for which the reaction temperature can be expressed solely as a function of concentrations. In the next section we introduce micromixing by replacing the PFR model by an axial dispersion reactor (ADR), a tubular reactor with mixing in the axial direction. Mixing in an ADR is characterized by the dimensionless dispersion coefficient, D, whose lower (upper) bound corresponds to a PFR (CSTR). Throughout the rest of this paper D (the inverse of the Peclet number) will simply be called the “dispersion coefficient”. With nonisothermal reactions in mind, interstage cooling or heating will be used. Heat exchange with a jacketed reactor or heat exchanger is modeled by specifying a constant reactor jacket temperature and a dimensionless heat-transfer coefficient.
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986 L
PFR
623
-
PFR
and energy terms given by
q R = hAs(Ts - T) Figure 5. Initial and optimum structures for example 2. ?FR
PFR
VIP
PIE
k is the number of stoichiometrically independent reactions. For convenience, we nondimensionalize eq 1-3 by introducing the following variables and relations: t = y/L, 0 = T/To, DT = K,/uL D = D/uL,
PiR
-
Figure 6. Initial and optimum structures for example 3.
3. Mass and Energy Balances The following quantities will be defined for the development in the following sections: u = average axial velocity in the reactor D = turbulent diffusion coefficient Ks = thermal diffusivity in the axial direction y = position along reactor T = temperature C = vector of concentrations of reacting species R = vector of reaction rate expressions p = average density of reaction mixture C, = heat capacity per mole of reaction mixture qE = heat exchanged with surroundings qR = heat produced by reaction h = overall heat-transfer coefficient based on bulk-wall temperature difference A, = effective area for heat transfer through wall of reactor T, = coolant or heat source temperature (assumed constant) To = reference temperature rp = rate of reference reactant of pth independent reaction AHp = heat of reaction (based on one reactant) of pth independent reaction Q = volumetric flow rate in reactor Qo = volumetric flow rate entering reactor network Co = concentration of one reactant (used as a reference) in initial feed (note that Co is scalar) L = reference length Consider the nonisothermal axial dispersion reactor model (ADR) in Figure 2 with constant temperature for the heating or cooling jacket. We now restrict the system to homogeneous reactions only. Then a t steady state, assuming radial concentration and temperature gradients are negligible and that the density of the reaction mixture takes on a constant average value, the species continuity equations and the energy balance are those given by Carberry (1976) dC d2C u- D= R(C,T) ytm Iy IY;+~,,, ( l a ) dy dy2
with boundary conditions (due to Danckwerts, 1953)
(3b)
r
X = QC/(QoCo),
W = T,/To
=L~A,/(~Pc~ G )( x, , ~ )= ( L / u T , ~ c , ) ~ ~
F(X$) = QLR(x,O)/ (uQoCo) (4) We assume also that DT a D, i.e., let DT = tD, where t is a positive constant. To determine the adjoint system, it is convenient to rewrite eq 1 as the following system of first-order differential equations dX/dt = Z
tlm It I
(54
dZ/dt = [ Z - F(X,0)]/ D
(5b)
d8/dt = V
(5c)
dV/dt = [V - G(X,B) - I'(W
-
@)]/(ED) (5d)
with boundary conditions x(t$J
= x(t,m)
+ Dz(tlm)
(64
z(G+l,rn) =0 0(tLm) = 0(t$ + tDV(t$J
(6b)
V(tT+1,m) = 0
(64
(6c)
4. The Optimization Problem Consider M reaction paths such that the mth path has N sources and/or sinks (excluding the origin). For the reactor (i, m) shown in Figure 3 the boundary conditions can be written as
x(t&J
= x(ti,m) + Di,mz(t&)
where if ti,mis a sink: X(ti,,) = X(t;,)
+ ci,mX(f;m)
if ti,mis a source: X(ti,,) = (I - ~ i , ~ ) X ( t ; ~ ) The tilde will be used to relate corresponding sink and source points, t,,n and ti,m,and will be used to differentiate between a sink and a source when the context does not make it clear. For example ti,m= and Zi,, = t,,? in Figure 1; ci,,, is the separation matrix of the corresponding source. Since the model must allow for preheating or precooling of the feed to reactor (i, m),we make the inlet temperature 8(tiim)a decision variable by setting i t equal to a decision variable, Y ~ , ~In. this study we follow Jackson and use an objective function that depends only on the sum of the flows a t the end ( t = tNm)of the reaction paths. More explicitly, we define the vector p by M
P =
m=l
X(tkm,m)
Given a superstructure composed of the segments shown
624
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986
in Figure 2, the following optimization problem can be posed: maximize (or minimize) J = 9 ( p ) subject to the state equations for tzm 5 t 5
dX/dt = 2 dZ/dt = [ Z - F(X.O)]/D
i = 0 , ..., ( N ,
(7a) - 1)
(7b)
d6/dt = V
(7c)
dV/dt = [V - G(X,O) - r ( W - O)]/(tD)
(7d)
and boundary conditions for each segment (t,,,, t,+l,m) if ti,mis the origin: X(t$,,) = X(t,,)
+ D,,,Z(t&,J
(7e)
if ti,mis a sink:
X(tzm) = x(t,m) + ci,mX(zi.m) + Di,J'(t{m) (7f) if ti,mis a source: (7g) x(t{m) = (1 - ui,m)X(tLm)+ Di,mz(t{mt,) z(t;+l,m) = 0 fl(t{m)
=
7i.m
+ tDi,mv(t&J
V(tY+1,m) = 0 to,m
=0
(7h)
(7k)
5. A S u m m a r y of the Adjoint System
To solve the above nonlinear program (NLP), a strategy is proposed in which the decision variables, which represent the degrees of freedom in the system, will be chosen iteratively using a gradient-based optimization algorithm. To apply this method, we need to develop relations for adjoint variables (Lagrange multipliers associated with the differential equations) which reflect the sensitivity of the objective function to changes in the decision variables tNm,,, ti,m,ui,,, Y~,,, Wi,,, ri,,,and Di,m. The derivation of this system of equations is carried out by using the standard approach of adjoining the objective function, J , with the state equations and boundary conditions via position-dependent multipliers, XA, AB, &, AD, and constant multipliers, v, A, p , respectively, to form a Hamiltonian function. Employing variational techniques, as outlined in Bryson and Ho (1975) for example, and defining these multipliers so that the Hamiltonian is directly influenced by the decision variables result in the following system of adjoint equations dXA_-_ dt 1 D
+ (C,G/t)X,]
for tzm I t 5
= (1 - "i+l,rn)TXA(tAl,m)+
o?+,,mXA(~~l,m)
(8i) AC(tL+l,m)
(7j)
where the dimensionless feed, X ( t , , ) , is fixed at a given value. Note that in Figure 3 ti,mis a source and is a sink point. Note also that the introduction of yi,, makes it unnecessary to consider temperature relations at source and sink points. Finally, from physical considerations the following bounds on the decision variables are also imposed: (a) 0 I rmin (adiabatic) I I?;,, I rmax (isothermal at temperature Wi,,) < m (bound on the heat-transfer coefficient). (b) 0 IWminIWi,mI W,,, I m (bounds on the temperature of cooling or heating jacket). (c) 0 IDmin(PFR)I Di,mI D,,(CSTR) C m (bounds on the dispersion coefficient). (d) 0 I7min I Y ~ I, Y,,~ ~ < m (bounds on inlet temperatures). (e) 0 Iui,m II, 0 I ti,m5 t,Lrm,m < m.
-[(~,F)XB
XA(t;+l,m)
(7i)
(8a)
=0
(83
The gradient of the objective function, J, with respect to the decision variables tNm,m, ti,,, ~ i , Wi,,, ~ , ri,m, and Di,,, is given in terms of the adjoint variables as M
dJ =
C H(tkm,,) dt,,,,
m=l
+
where H is the Hamiltonian function H = Xi2 + Xg[Z - F(X,6)]/D + AcV + X D [ V - G(X,fl) - r ( W - O)]/(tD)
c,H = -(r/(tD))XD V r H = X D ( ~- W)/(tD) YDH = -Xg[Z - F(X,O)]/ D 2 -
XD[V- G(X,O) + r(8- W)]/(cD2)
6. Proposed Algorithm Let Y be the vector of decision variables tNm,,, t,,,, u,,,, Y ~ , W,,,, ~ , f l , m , and Q m . Then in terms of Y the NLP takes the simpler form max (or min) 9 = 9(Y) (a) subject to (b) Ymm 5 Y 5 Ymax (C) S,,, 5 A Y I S ,,, where A is the coefficient matrix. The linear constraint set c may involve some or all the decision variables. For example, we may require the source and sink points not to cross or source and sink points to be located before the terminus. Note also that we can replace (c) by a more general nonlinear constraint
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986 625
S(Y)5 Sm,
(C*) at the expense of making the algorithm only slightly more complicated. Based on the above derivation, the following algorithm is proposed: (i) Choose a superstructure of reactor networks, as shown in Figure 1. (ii) Provide an initial guess for the decision variables ti,m, ui,m!.Yi,m, Wi,w I’i,rn, and Di,m* (111) For each reactor segment (ti,rn, ti+l,rn) solve the state equations (7), which are a set of nonlinear two-point boundary value ODE’S. Here we use COLSYS (Ascher et al., 1979), a TPBVP ODE solver that uses orthogonal collocation on finite elements. Store X,2, 8, V, and the objective function, J. (iv) With the state variables from (iii), solve the linear adjoint equations (8) with COLSYS. (v) Use eq 9 to calculate the derivatives of the objective function with respect to the above decision variables. (vi) Calculate a search direction for the decision variables by setting up and solving the following quadratic program at Yk: max (or min) Q = V @ ( Y k ) T d + k ‘/zd;fBkdk Smin
5
dk
subject to
Smin-AYk 5 Adk 5
smmAYk
Here Bk is a positive definite matrix constructed by a quasi-Newton updating formula; it serves as an approximation to the Hessian of @(Yk). With the search direction d k , update Y as Yk+l = Yk + akdk,where f f k is a step size selected so that a sufficient decrease in @ is found at Yk+l. This successive quadratic programming (SQP) algorithm is presented by Han (1977) and Powell (1977). A complete description of an updated algorithm is given in Eiegler and Cuthrell (1985). (vii) Stop if lldkll < 6, where 6 is a specified tolerance. The necessary optimality conditions are satisfied for the NLP. If lldkll > 6, go to step iii. Note that, depending on the general nature of the objective function, J, there could be multiple local optima for the NLP. To gain confidence that a global optimum has been found, restart the algorithm from different initial points. However, unless the objective function exhibits special features, such as convexity in the decision variables, there is no guarantee that a global optimum can be found.
7. Test Examples In order to illustrate the properties and performance of the above algorithm, we consider the following test problems. All are drawn from the literature and involve interesting and nontrivial kinetic mechanisms. Our aim is to show that the above algorithm finds optimal reactor networks efficiently and, in some cases, yields better solutions than those reported in the literature. Due to space limitations, presentation of the examples will be brief. The reaction models for these examples can be found in the corresponding references. Example 1. Van de Vusse Reaction. kl
A-B-C
It J k,
D
k2
The isothermal version of the constant dispersion model has been tested on the Van de Vusse reaction (1964) (see Chitra et al., 1981). We now consider two different cases for this mechanism in order to demonstrate the algorithm and properties of the superstructure. For this problem the objective function is the yield of B and the kinetic parameters are kl = 10 s-l (first order), k2 = l s-l (first order), k, = 1 L/(g-mol s) (second order), and feed flow rate (pure A) = 100 L/s. The reaction rate expressions are given in Chitra and Govind (1981). The first case below was also considered by Chitra and Govind while the second involves a more complex superstructure consisting of two reaction paths. ( a )Feed Concentration of A = 0.58 g-mol/L. Here, the initially guessed reactor configuration has two CSTR’s in series. The final optimum structure is a single PFR (see Figure 4, top) with the following characteristics (Chitra and Govind’s result is shown in parentheses): maximum B yield = 0.4368 g-mol/L (0.4362 g-mol/L) VpFR = 29.65 L DCsTR = 50
DpFR = 0.001
( b ) Feed Concentration of A = 0.58 g-mol/L. The initially guessed configuration has two reaction paths with two CSTRs in each. In addition, the two paths have one side stream connecting them. The final optimum structure has two PFR’s in series. All of B is removed from the reaction path at the exit of the first PFR (see Figure 4, bottom) with the following characteristics: maximum B yield = 0.4969 g-mol/L u = diag [0, 1, 01 VpF,, = 4.412 L Vp,,, = 6.503 L
Example 2. Selectivity Optimization. The following isothermal reaction is due to Trambouze et al. (1959) (see also Figure 5):
A
xB - c
The kinetic constants and other parameters are kl = 0.025 g-mol/(L min) (zero order), k2 = 0.2 min-’ (first order), k3 = 0.4 L/(g-mol min) (second order), feed flow rate (pure A) = 100 L/min, A,, = 1g-mol/L, DpFR = 0.001, and Dc~TR = 20. The objective function assumed is the selectivity of A to C defined by Xc/(l - X A ) ,where X (the dimensionless concentration) is defined in eq 4 of section 3. Paynter et al. (1970) solved this example using a dispersion model in which D varied continuously with the position along the reactor. They reported an optimal configuration of one CSTR in series with a PFR. For comparison, their results are shown in parentheses. Our initially guessed reactor configuration is two PFR’s in series. Two CSTRs in series are found to be optimum. The results are maximum selectivity = 0.4999 (0.4951, TCSTR,J = 7.721 min, and TCSTR,2 = 0.0975 min (TCSTR = 9.1 min, rpFR= 0.9 min) where 7 is the residence time. In addition, the accuracy of our results was checked by solving the CSTR equations with the above values for residence times. Example 3. Adiabatic Naphthalene Oxidation (Chitra et al., 1985). This is a highly exothermic reaction
626
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986
A-B
kl
k4
where A = naphthalene, B = naphthaquinone, C = phthalic anhydride, and D = COz + HzO. The objective function is the yield of C with rate constants given as k, = k2 = 2.0 X lo1, exp(-38000/RT) h-l, k, = 8.15 X l O I 7 exp(-50000/RT) h-l, k, = 2.1 X lo5 exp(-20000/RT) h-l, [PI, P2, &, p,] = [0.12, 0.43, 0.36, 0.741, where Pi = (aHi)A,/(PCPTfee,-J = dimensionless heat of reaction, R = 1.987 cal/(g-mol K), and T = temperature in K. All reactions are first-order with the following information: feed flow rate (pure A) = los L/h, feed concentration of A (A,) = 1.0 g-mol/L, and Tfeed= 1018 K. Tfeed was assumed to be fixed, and no attempt was made to optimize it. The initial reactor configuration has four PFR’s in series and allows for intercooling between reactors. The final optimum structure has one PFR with an exit temperature of 1481 K. The initial and final configurations are given in Figure 6 with maximum C yield = 0.9999 g-mol/L (0.995 g-mol/L), VpFR = 9.5 L, DCsTR = 50, and DpFR = 0.001. (Chitra and Govind’s results are given in parentheses.) In example 1 the results compare favorably with those of Chitra and Govind. In general, we found slight improvements in the objective function. Moreover, the fact that our results were obtained from starting points far from the optimum makes our approach even more viable. In example lb, use of two reaction paths and one side stream led to a smaller total reactor volume (10.92 compared to 29.65 L) and a larger B yield (0.4969 vs. 0.4362 g-mol/L). Note that the theoretical bound on the B yield is 0.58 g-mol/L. The results in (lb) are not surprising because one way to increase the B yield is to prevent it from reacting to form C by removing B from the mixture as soon as it is formed. This seems to suggest that we include more side streams in the superstructure to make removal of B possible. Example 2 was also reworked by using an initially guessed configuration of two reaction paths, each with two P F R s in series. The two paths were connected by one side stream (just as in example lb). However, the optimal configuration and the objective did not change. Thus, one reaction path is probably adequate for the chosen objective function, selectivity of A to C. This example also shows the validity of our approach. Finally, for the naphthalene oxidation reactions in example 3, lz, and k, are much larger than k,. As a result, the C to D reaction is negligible. Since all the other reactions lead to the production of C, it is to be expected that there will be virtually complete conversion of A to C. All of the examples require from four to ten iterations of the SQP optimization algorithm to converge. The most time-consuming part of these computations lies in solving the two-point boundary value problems with COLSYS (see step iii in the algorithm). COLSYS applies Newton’s method to collocation equations over finite elements, and for certain values of the decision variables it may be difficult for this method to converge. One reason for this difficulty is that for D (the PFR limit) the solution has a “boundary layer”, a small region near the reactor exit in which the derivative of concentration with respect to position rapidly goes to zero in order to satisfy the right-
end boundary condition. For COLSYS to converge properly, the mesh selection algorithm has to be able to recognize this boundary layer and tailor the mesh size in that segment. To deal with nonconvergence, we tried a combination of the following: (i) change the initial mesh sent to COLSYS, (ii) change the minimum step size that is used in the Newton iteration step in COLSYS, and (iii) rescale the problem (presently done by inspection) by a judicious choice of reference length, temperature, and so on. Since most of the computation time was spent in solving the nonlinear system differential equations, we are looking for better ways to scale the ODE’s and ways to accelerate convergence in COLSYS. Also, we intend to consider modifying the formulation of the algorithm and using other TPBVP methods. 8. Conclusions and Recommendations for Future Work Of the advantages of the constant-dispersion model, the ones that stand out are the following: 1. The adjoint variables that are needed to update the parameters in the optimization are calculated relatively quickly by using COLSYS, once the state equations are solved. This is because the adjoint system of ODE’s is linear. 2. A wide range of objective functions and kinetic mechanisms can be handled with this approach. 3. Allowance is made for considering very general nonisothermal systems. 4. It is very easy to modify this formulation to handle nonlinear side constraints and other process parameters such as the reactor volume. The main drawbacks with the constant dispersion model are the following: 1. We do not yet know how to interpret reactor networks that include segments with intermediate dispersion coefficients. 2. The choices of Dmin(Dmm), the lower (upper) bound on the dispersion coefficient, are somewhat arbitrary and can have a big impact on the performance of COLSYS. In the above examples, however, we checked our bounds on D with the actual CSTR and PFR models and found excellent agreement. As a result of these drawbacks, a second model has been developed which uses a recycle reactor, instead of a dispersion reactor, as the basic unit of the reactor network. The goal is to find some comparison between the two formulations. This will be reported in a later paper. The objective function is also being expanded to include the split fractions and the reactor volumes. Thus, it will be possible to incorporate the cost of separation at source points and the cost of reactor volumes into the objective. Also since most of the computation time (about 85%) is spent in solving the nonlinear ODE’s of the state variables, it is of utmost importance to find ways of accelerating convergence in COLSYS. The final question that has not been posed so far is, how does one find a superstructure? The present strategy is to use a trial-and-error approach based on heuristics. A more systematic procedure that incorporates guidelines based on characteristics of the specific reaction system will be studied in the future.
Acknowledgment Financial support under NSF Grant No. ENG-8451058 is gratefully acknowledged. A detailed derivation of the adjoint system and proofs for the PFR/CSTR limits can be found in Achenie and Biegler (1986), which is available on request from the authors.
627
Ind. Eng. Chem. Fundam. 1986, 25, 627-633
Literature Cited Achenie, L. E. K.; Biegler, L. T. EDRC Report, 1986; Carnegie-Mellon University, Pittsburgh, PA. Ascher, U.; Christiansen, J.; Russell, R. D. Mafh. Compuf. 1979, 33, 659-679. Bielger. L. T.; Cuthrell, J. E. Comput. Chem. Eng. 1985, 9 ,257-267. Bryson, A. E.: Ho, Y.-C. Applied Optimal Control; Hemisphere: New york, 1975. Carberry. J. J. Chemical and Catalytic Reaction Engineering; McGraw-Hill: New York, 1976. Chitra, S.P.; Govind, R . Chem. Eng. Sci. 1981, 36, 1219-1225. Chitra, S.P.: Govind, R . AIChE J. 1985, 3 1 . 177-193. Conti, G. A. P.; Paterson, W. R . Inst. Chem. Eng. symp. Ser. $985, ,v0,
Jackson, R . J. O D ~Theor. . A m / . 1968, 2. 240-259 Neyfeh, A. Pertuibation MethGs; Wiley: New York, 1973. Nishida, N.; Stephanopoulos, G.: Westerberg, A. W. AIChE J. 1981, 27, 321-350. O'Malley, .~R. E. Introduction to Singular Perturbations; Academic: New York, 1974. PaYnter, J. D.; Haskins, D. E. Chem. Eng. Sci. 1970, 25, 1415-1422, Powell, M. J. D. Lect. Notes Math. 1977, No. 630, 144-155. Ravimohan, A. L. J . Opt. Theor. 1971, 8 , 204-211. Sargent, R. W. H.; Sullivan, G. R. Proc. 8th I F I P Conf. O p t . Tech. 1977, Part 2, 158-168. Trambouze, P. J.; Piret, E. L. AIChE J. 1959, 5 , 384-390. Van de Vusse, J. G. Chem. Eng. Sci. 1964, 79,994-999.
92 . - , 391-402 -.
Danckwerts. P . V . Chem. Eng. Sci. 1953, 2, 1-18. Han, S.-P. J. Opt. Theor. Appl. 1977, 22, 297-309. Horn, F. J. M.; Tsai, M. J. J . Opt. Theo. Appl. 1967, 7 , 131-145.
Received for review June 18, 1986 Accepted July 18, 1986
Changes in Shear-Induced Hemolysis due to Blood Additives: Synthetic Biopolymers and Asthma Drugs Robert C. Lljana,+Joseph M. Monroe,$and Mlchael C. Williams" Chemical Engineering Department, University o f California, Berkeley, Berkeley, California 94720
Mechanical fragility of human red blood cells was evaluated by shearing the blood between rotating polyethylene disks and measuring the timedependent release of hemoglobin. Several blood additives were tested for their effect on this hemolysis: (a) two synthetic polynucleotides (polyadenylate and polycytidylate) as biopolymers with purine and pyrimidine moieties, respectively; (b) two low-molecular-weight purine derivatives (the drug theophylline and uric acid). It was found that polyadenylate always increased hemolysis, polycytidylate often reduced it, theophylline always reduced it, and uric acid was always ineffective. Drug localization data on theophylline showed a large uptake of the additive by cell membranes, the degree of hemolysis protection being proportional to the mass of the drug absorbed. Supplementary cell characterization by resistive pulse spectroscopy documented that the protective drugs caused cell volumes to increase, deformability to decrease, and osmotic fragility to decrease. Chemical and mechanical mechanisms for changes in mechanical fragility are proposed.
Introduction Flow of blood through artificial organs and implanted conduits leads to various forms of blood degradation. The release of erythrocyte hemoglobin is often used as a general measure of red cell damage. In previous work, this laboratory has examined how flow-induced hemolysis is influenced by hydrodynamic factors (Shapiro and Williams, 1970) and by the roughness (Monroe et al., 1981) and surface composition of conduit materials (Offeman and Williams 1979a; Monroe et al., 1980: Offeman and Williams, 1976). Reviews of the topic have also appeared (Bernstein, 1971; Blackshear, 1972; Hellums and Brown, 1977). Few investigations have been directed toward the chemical features of blood that dispose it to hemolyze in shear flow. Hemolysis tendencies have been correlated with preshear plasma chemistry (Offeman and Williams, 1976), statistical correlations being found with elevated concentrations of lactate dehydrogenase (LDH), very low density lipoproteins, and uric acid. A subsequent study (Monroe et al., 1980) of plasma chemistry changes during in vitro shear-induced hemolysis revealed a drop in glucose and rise in uric acid (both presumably due to shear-induced metabolic activity), while Ca2*and Na+ entered the cells and LDH was released. Procter and Gamble Co., Miami Valley Laboratories, Cincinnati, OH. *Union Oil Company, Los Angeles, CA.
The concept that chemical additives to the plasma might reduce hemolysis by some sort of protective action on erythrocytes is relatively new. Blood preservatives have been designed to preserve essential features of cell viability but not with attention to the tendency to hemolyze in subsequent shear flow. The introduction of adenosine into some blood preservatives did, however, reduce such hemolysis (Tadano et al., 1977). It seems likely that other additives could have similar or superior results, and this possibility was the principal motivation for our extended work on the effects of blood additives on shear-induced hemolysis. Studies of plasma proteins, used for precoating synthetic surfaces and added to red cell suspensions, indicated that these substances also afforded hemolytic protection (Nichols and Williams, 1976). Their role in the precoating tests was clear: they transformed the synthetics into surfaces resembling biological ones which induced less chemical trauma in red cells contacting them. The beneficial result of adding proteins to the fluid environment of erythrocytes was less clear. (The most successful was y-globulin.) In a medium deprived of proteins, enrichment by protein supplementation is a step toward restoring a normal (physiological) environment and, thus, returning the cell condition to one of normal hemolytic resistance. Still, this picture fails to address the specifics of chemical interaction between protein and cell that lead to hemolytic resistance. The present study, therefore, was also designed to investigate the role of simple homopolymers, synthesized from a single bio-"mer", having far fewer chemical com-
0196-4313/86/1025-0627$01.50/00 1986 American Chemical Society
