THESIS J LIBRARY MlChlgah State fl , ‘J" uvv‘r '— ~— This is to certify that the thesis entitled YIELD OPTIMIZATION ANALYSIS OF COMPLEX CHEMICAL REACTIONS IN ISOTHERMAL CHEMICAL REACTORS UNDER FORCED PERIODIC OPERATION presented by Richard Dale Skeirik has been accepted towards fulfillment of the requirements for M. S . degree in Chemical En qi neering m flew / Major professor Date Feb. 10, 1981 0-7639 1/ OVERDUE FINES: 25¢ per day per ital RETQRNING LIQRARY MATERIALS: P'Iacc in book return to remove charge from circulation records YIELD OPTIMIZATION ANALYSIS OF COMPLEX CHEMICAL REACTIONS IN ISOTHERMAL CHEMICAL REACTORS UNDER FORCED PERIODIC OPERATION By Richard Daie Skeirik A THESIS Submitted to Michigan State University in partial fuifiiiment of the requirements for the degree of MASTER OF SCIENCE Department of Chemical Engineering 1981 ABSTRACT YIELD OPTIMIZATION ANALYSIS OF COMPLEX CHEMICAL REACTIONS IN ISOTHERMAL CHEMICAL REACTORS UNDER FORCED PERIODIC OPERATION By Richard Dale Skeirik Periodic forcing of feed concentration to isothermal CSTR's with Van de Vusse kinetics has been examined. Square wave and sinusoidal forcing to homogeneous reactors and square wave forcing to a CSTR with one dimensional diffusion inside the catalyst pellets have been considered. For the homogenous reactors, analytic or semi- analytic solutions have been developed for species concentrations as well as time averages. The heterogeneous reactor was treated numerically. All cases showed a decrease in yield of the inter- mediate and an increase in yield of the side product. The effects of periodic operation increase as the forcing frequency decreases, hence the maximum effect of forcing can be achieved by a smaller steady state reactor. TABLE OF CONTENTS Chapter I. PERIODIC SQUARE WAVE FEED CONCENTRATION TO A HOMOGENEOUS ISOTHERMAL CSTR WITH VAN DE VUSSE REACTION NETWORK .................. l The Choice of the System The Method of Comparison The Periodic Scheme First Order Kinetics Van de Vusse Kinetics Calculative Results Conclusions 11. PERIODIC SINUSOIDAL FEED CONCENTRATION TO A HOMOGENEOUS ISOTHERMAL CSTR WITH A VAN DE VUSSE REACTION NETWORK .................. 40 The Method of Comparison Consecutive Parallel Reactions Van de Vusse Kinetics Numerical Results Conclusion III. PERIODIC SQUARE WAVE FEED CONCENTRATION TO AN ISOTHERMAL CSTR WITH A VAN DE VUSSE REACTION OCCURRING IN POROUS SLAB CATALYST PARTICLES ........ 82 Numerical Results Conclusion LIST OF REFERENCES ....................... 104 ii LIST OF FIGURES Chapter 1 Page Figure l.(a-d) Behavior of the Van de Vusse system for a case of intermediate cycling frequency. K1=K2=K3=],eo=]. .............. 24 2.(a-d) Behavior of the Van de Vusse system for a case of slow cycling frequency. K1 = K2 = K3 = l, 90 = l0. .................... 29 3.(a-d) Behavior of the Van de Vusse system for a case of rapid cycling frequency. K; = K2 = K3=],eo=.3 .................... 34 Chapter 2 Figure l.(a-g) Behavior of Van de Vusse system for Tkl = k2/k1= Tkng = A =1, E = .1 ............ 60 2.(a-g) Behavior of Van de Vusse system for Tkl = k2/k1 = Tk3A0 = A = l, c = l ............. 63 3.(a-b) Behavior of Van de Vusse system for 1k; = kz/kl = Tk3Ao = A = l, e = .75 ............ 66 4.(a-b) Behavior of Van de Vusse system for Tk1 = k2/k1 = TkaAo = A = l, e = .25 ............ 66 5.(a-b) Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k2/k1 =(T Tk3Ao = 5 A = T e = .25. % error = C - . pert num Cnum lOO ....................... 69 6.(a-b) Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = kz/kl = l rkaAo = 5 A - l c = .75. % error = (C - c ~1oo .......... PeTt. . '3‘“? ....... 69 8.(a-b) 9.(a—b) lO.(a-b) ll.(a-b) 12.(a-b) l3.(a-b) l4.(a-b) Chapter 3 Figure 1. Deviation of perturbation solution from numerical simulation for Van de Vusse system Tk1 = k2/k1 = 5 TksAo = I A = l 25. % error = (C - C )/ . pert num Cnum lOO ....................... Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k2/k1 = 5 TkaAo = I A = l c = .75. % error = (C - C )/ cnum 100 .......... Peft. . PUT ....... Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = kz/kl = TkaAO 5 A = T E = .25. % error = (Cpert'o Cnum)/Cnum°lOO ........... Deviation of perturbation solution from numerical simulation for Van de Vusse system 1k; = kZ/kl = Tk3 A0 =5 A = l c = .75. % error = (C Cnum)/Cnum°100 ........... Deviation of perturbation solution from numerical simulation for Van de Vusse system rkl = k2/k1 = TkaA0C =l A = 5105 = .25. % error = (Cpert'o Cnu mnu)/C ........... pert Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k 2/k1 = IkgA0C =l =5 0c = .75. % error = (Cpert'o Cnu m)/Cnu ........... Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k2/k1 = TkaAo 1 A = .2 E = .25. A error = (Cpert'ocnum)/Cm100 ........... Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = kz/kl = Tkng l A = .2 c = .75. % error = (Cpert'ocnum)/Cm100 ........... Bulk phase reactant concentration vs. time for square wave cycling. ¢=l. v=l. B=l. K3=l. h=l. K2=l. A=l. w=l ............. iv 72 76 77 90 Bulk phase intermediate concentration vs. time for square wave cycling. ¢=l. v=l. B=1. K3=1. h=l. K2=l. A=I. w=l ....... Bulk phase plane behavior for square wave cycling. ¢=l. v=l. B=l. K3=l. h=l. 1. K2: . A= Reactant concentration profile within the pellet at five selected points in the limit cycle. ¢=l. v=l. B=l. K3=l. h=l. K2=T. A=1. w=l ............. Intermediate concentration profile within the pellet at five selected pointed in the limit cycle. ¢=l. v=l. B=l. K3=l. h=l. K2=l. A=I. w=1 ............. w=l ................ 92 96 98 LIST OF TABLES Chapter 1 l. Enhancement of time average periodic concentrations over steady state concen- trations. Approximation to consecutive reactions when K320. Approximation to second order reaction when K120, K2=O ...... Chapter 2 1. Percent deviation of ap roximate time average concentrations Eequations (24)] from numerically simulated time average concentrations; percent enhancement of approximate time average concentrations over steady state concentrations ........ Chapter 3 1. Percent increase in time average of periodic bulk phase concentrations over steady state concentrations for square wave perturbation ................ vi Page 20 ..... TOO NOTATIONS--CHAPTER 1 Concentration of reactant species Reference concentration Concentration of intermediate species B/Ao Time average b for periodic limit cycle Product species Constants, defined by eqns. (l4c),(20c-20h)(200)(l6e) Concentration of product species D/Ao Time average d for periodic limit cycle Locus of zero lepe in b-f or d-f phase plane A/Ao Initial condition at beginning of Region l Initial condition at beginning of Region 2 Time average f for periodic limit cycle Represents behavior of b(f) or d(f) in proof Dimensionless rate constant eqns (4), (12) Reaction rate constant Volumetric flow rate Time Period of l/2 of feed cycle CSTR volume Subscripts f Feed condition i Denotes interchangability of Regions l and 2 ss Steady state 1 Portion of feed cycle in which dimensionless feed concentration is 2.0 2 Portion of feed cycle in which dimensionless feed concentration is 0.0 w Approaching infinity Greek Symbols jeo 60, 360, 660,... n60 0, 280, 460,... d(neo) Defined by eqn (l4d) d(n60)m Defined by eqn (25f) d(jao) Defined by eqn (14e) 6(j60)m Defined by eqn (259) B Constant in eqns (18) O t/T eo to/T E Defined by eqns (2) I V/g ¢i Constants, defined by eqns (20i)-(20n) wi Substitution parameter viii NOTATIONS--CHAPTER 2 Concentration of reactant species Constant, defined by eqn (lSn) Constants, defined by eqn (l5u) Constants, defined by eqns (6) Concentration of intermediate species B/Ao Time average b Constants, defined by eqn (l5v) Concentration of product species C/Ao Time average c Constants, defined by eqns (lO), (l7d) Constants, defined by eqn (15w) Concentration of product species D/A0 Time average d A/A0 Time average f Integration constants Reaction rate constant Volumetric flow rate ix P. Periodic functions t Time V CSTR volume Greek Symbols a Defined by (lSn) “i Constants, defined by eqn (l5x) Bi Constants, defined by eqn (l5p) Yi Constants, defined by eqn (l5q) c Adjustable parameter controlling the amplitude of the input disturbance--perturbation parameter for Van de Vusse kinetics ci Defined by eqns (195)-(l9w) n1 Constants, defined by eqn (lSr) B t/I A1,i Defined by eqn (lSn) A wT “i Constant, defined by eqn (lBe) ii Constant, defined by eqn (le) v1 Constants, defined by eqn (15$) Oi Constants, defined by eqn (l5t) T V/q Xi Constants defined by eqns (235)-(23w) Oi Substitution parameter w Frequency of input disturbance u£,i Defined by eqn (18d) Subscripts f Feed condition 0 Reference concentration 5 Steady state m Approaching infinity xi NOTATIONS--CHAPTER 3 Concentration of reactant species Spatial distance within catalyst particle Effective diffusivity of species A Effective diffusivity of species 8 Time Effective rate constant based on concentration in bulk phase adjacent to surface Volumetric flow rate External surface of one catalyst pellet perpendicular to mass flux Number of catalyst particles in CSTR Concentration of intermediate species l/2 thickness of catalyst slab A/Ao, eqn (3a) kng/ki, eqn (3a) B/A0 kz/k1 D /0 EA 98 xii Greek Symbols E .(‘u A superscripts A Subscripts 0 f Fraction of CSTR volume not occupied by catalyst pellets Thiele modulus, eqn (3a) X/L t D /L2 8A vs 0 /qL2 eA De npA/qL A k L D 9 / 9A Dimensionless frequency of input disturbance bulk phase Reference concentration Feed xiii CHAPTER I PERIODIC SQUARE WAVE FEED CONCENTRATION TO A HOMOGENEOUS ISOTHERMAL CSTR WITH VAN DE VUSSE REACTION NETWORK 2 Chemical reactors under forced periodic operation often give a time average output which is different from the corresponding steady state. In general, this is true whenever some of the output variables depend in a non-linear way on the forcing input variable. For those selectivity reactions which exhibit an absolute maximum yield or selectivity under steady state operation, forced periodic operation is a promising method of improving on steady state performance. A large segment of the work done in this area has centered on non- isothermal systems, which allows the exploitation of the exponential non-linearity of the Arrhenius rate expression. However, the complexity of the simultaneous heat and mass balances combined with at least one exponential non-linearity virtually eliminates the possibility of obtaining an analytic solution. Hence the treatments have been either approximate, usually through the use of some kind of perturbation technique, or numerical. While Douglas and coworkers [l] [2] [3] [4], Renken [5] [6], and Farhadpour and Gibalaro [l4] have studied specific reacting systems under forced periodic Operation, no one has yet obtained a fully analytic prediction of system behavior. In contrast, this work treats complex kinetics occurring in an iso- thermal CSTR with a square wave feed concentration and provides fully analytic solutions for the reactant and product concentrations as well as their time averages. THE CHOICE OF THE SYSTEM In deciding on a complex kinetics system, the number of possible choices is infinite. Some important characteristics led to the choice 3 of the Van de Vusse system, A+B+C, 2A+D, for this work. As Van de Vusse proposed, the system cannot be easily classified by general reactor selection guidelines. The intermediate B exhibits an absolute maximum concentration under steady state operation [8], thus this is a system which explores the real potential of periodic operation. Probably most important from the perspective of enhancing reactor performance, under some reactor conditions,the intermediate B increases through a transient maximum, even though it is already above the steady state it is approaching. Hence periodic operation provides an inviting method to capitalize on this phenomenon. Any initial optimism about this aspect of the system must be tempered, however, by the fact that B also goes through a transient minimum even when below the eventual steady state [7]. THE METHOD OF COMPARISON In order to draw any meaningful conclusions from the study of periodic Operation of chemical reactors, comparison against a corresponding steady state is necessary. In addition, the comparison is subject to certain restrictions. (l) Residence Time: (a) For those products for which steady state operation can be made more effective by simply using a longer residence time (for example, the reactant A in the Van de Vusse kinetics scheme), the comparison should be based on the same residence time for periodic and steady state operation. (b) For those species which exhibit an absolute maximum concen- tration with respect to residence time in steady state mode, the comparison should be between the optimum steady state concentration and the optimum periodic concentration, regard- less of the periodic process residence time. 4 (2) Equal time average molar flow rates of reactant species into the reactor should be used in both periodic and steady state mode. (3) Initial transients are ignored. For steady state (i.e., non- periodic) mode, only the steady state concentrations should be considered. For the periodic mode, time averages should be taken only after the effluent concentrations achieve a stable repetitive cycle (limit cycle). (4) For the feed concentration disturbance used here, volumetric flow rate in the periodic process should not vary with time. THE PERIODIC SCHEME The input disturbance used here consists of a symmetric square wave reactant feed concentration. The feed oscillates between dimensionless concentrations of 0.0 and 2.0. Each phase of the cycle lasts for a period to, giving the inlet disturbance an overall period of Zto. The steady state against which the periodic process is compared is simply the same process with an invariant dimensionless feed concentration of l.0. To simplify the solution, the process is arbitrarily defined to begin operation at time t=O with a feed concentration of 2.0 at the beginning of a feed cycle. The square wave input disturbance has several important character- istics. First, the input concentration during either half of the feed cycle is constant. The transient mass balance for the reactant species during either portion of the feed cycle can therefore be represented by an autonomous ordinary differential equation. If some other disturbance, such as a sine wave, were introduced, the reactant mass balance would become nonautonomous, and an analytic solution would most likely be impossible. 5 Second, because of the different inlet molar flow rates, the reactant mass balances for the two portions of the feed cycle are different, and there does not appear to be any way to combine the two. The mathematical treatment of the problem thus becomes a series of initial value problems, in which the final value of the solution in one region becomes the initial value for the solution in the next region. In order to do a rigorous analysis of the periodic process, the governing mass balances should be solvable, which for this kinetics virtually requires the system to be autonomous. However, the treatment of the system as a series of initial value problems presents some serious difficulties. The solution, although its value is always continuous, is not represented by a single function. During one portion of the feed cycle, the solution is given by one expression; during the other portion, by a different expression. It is not possible to conclude from inspection of the solution that the output will exhibit a stable periodicity as time becomes large. Hence, even after the solutions to the system concentrations as a function of time are known, a great deal of analysis remains to deter- mine how the system will behave as time progresses. Once it is shown that the system goes to a stable, repetitive cycle, or limit cycle, the con- centration functions can be integrated over one cycle to obtain the time average concentrations. The original objective of this work was to examine Van de Vusse kinetics. However, because of the apparent difficulty of proving the limiting behavior of the system, we first examined simple first order kinetics, and then extended the techniques of that analysis to Van de Vusse kinetics. In each case the approach to the problem was the same. First, full transient solutions for the species in question were obtained. Then, by 6 a qualitative examination of the governing differential equations, the limiting behavior of the system was established. This qualitative information about the system allowed calculation of the asymptotic limits of the output concentrations at the end of each region of the feed cycle. The transient solution for the ouptut concentrations were then integrated over the time period of one feed cycle, using the limiting values as initial and final conditions, giving analytical expressions for the time averages of these quantities. Definitions The treatment of the periodic process necessitates the definition of some terminology. Since there will be separate solutions to the mass balances for each part of the feed cycle, the two are distinguished by subscripts. The subscript l is used to designate that half of the feed cycle in which the dimensionless feed concentration is 2.0; this is identified as Region l. The subscript 2 designates the portion of the cycle with a feed concentration of 0.0; call this Region 2. Moreover, the solution takes the form of a series of functions, each of which depends on the value of the solution at the end of the previous region. Since each portion of the feed cycle lasts for a duration 60 (=to/T), and since we have defined the process to start at t=0 at the beginning of Region 1, then Region T will always end and Region 2 begin on an odd multiple of 60. Similarly, Region 2 will end and Region 1 begin on even multiples of 60. The symbols n00 0,260,460,... (13) jeo = 60,360,560,... (1b) allow a general representation of the beginning of either region. 7 The "initial" condition for any concentration function in a region is expressed by the subscript representing the region and the nomencla— ture for the beginning of the region i.e., eqns (l). Thus the solution to f(6) in Region l, beginning at n60, has its initial condition represented by f1(neo). Similarly f2(6) would have the initial condition f2(jeo). Although the solutions depend on time, it is only the time elapsed since the beginning of the region which is important. The substitution - e-neo ; Region l (2a) 4* l - O-jeo ; Region 2 (2b) 01 l is defined. FIRST ORDER KINETICS Although a first order reaction k1 A + B taking place in an isothermal CSTR shows an increase in conversion under periodic Operation, the methodology used for this case illustrates the methods used for Van de Vusse kinetics. The transient mass balance on the reactant species is given by qu - qA - k1VA = -—- (3) The problem is treated as two separate differential equations for the two regions of the feed cycle. Substituting the dimensionless variables - A_. = X. = I. = f-AO 1 q a T K1 1k, (4) we obtain the following two differential equations 0. f1 8 = 2 - (l+K1)f1 (5a) 0. 31—2—1: - (1+K1)f1 (5b) The solutions to equations (5) are given by f1(€) = fSSI- [f551-f1(ne°)] eXP [’(ITK1)€] (6a) f2(350) eXP ['(ITK1)€] (5b) f2(5) The solution to equation (5a) with the time derivative set equal to zero is System Behavior The behavior of the system for large a can be established by a qualitative examination of the governing differential equations, (5). First, ggl-< O, and f2 approaches zero asymptotically. Also giL <0,f1>f551 d6 - O ’ f1 - f$51 > O , f] < fSSI and f1 approaches fS asymptotically. $1 df dfz If f(6=0)>fsw both—L dB and do are O and dT < 0 for e>ec. 551 Hence we have an oscillation for f. df 551’ dB oscillation is immediately established. Here Be becomes immaterial. dfz For the case f(e=O)< f >0 and 36—'< O for all 6 and hence 9 At any time 6* (at the start of some cycle or half a cycle) during the oscillation, one of three events could occur, viz. f (e*+260) §=f(e*). By the property of uniqueness of ordinary differential equations, it becomes evident that f(6*)§f(6*+260)§f(6*+4eo), etc. Therefore, the system produces a characteristic of either increasing, decreasing, or ). unchanging sequence of initial conditions. Furthermore since fc(0,fSSl the sequence of initial conditions must approach an asymptotic limit. Hence, a limit cycle is established. The uniqueness of the limit cycle is discussed next. Asymptotic Limits and Time Averages Denoting the asymptotic limits at n60 and jeoby the symbols f1(n60)m and f2(jOo)m, we make use of the fact that if f2(jeo)co is the initial condition for f2(E), then its value after a time 80 must equal f1(neo)w. Similarly f1(n60)oo used as initial condition for Region 1 over an interval 60 must return the value f2(jeo)m. This gives f1(050)m = f2(jeo)co eXP [’(1+K1)eo] (7a) f2(.leo)co = [SSI’LT f1(neo)w]9XP['(]+K1)eo] (7b) The solution of this linear system is unique and is easily found to be 551' 7 _ l-EX - T+K1 O f‘(ne°)w ' f$51 [2 sinh l+K1 60 (88) f (.6 ) = f 9X 1+K1 O ’1 .1 (8b) 2 3 ° w 551 2 sinh l+K1 60] Once these limits are known, we can integrate over one full cycle using these values as initial conditions to determine the time average f. We define the time average, f as (:°f.e‘c“€+f551k. 3573.2.[an.a(e.>1+c.w.[a.a(neo)1} <15a> \ 3:0 where j 0 - 'Czj€_ ”Cui'. - _ C5j+1 1. e J9 CSJ f 1 WI[€OJOQJ = . (15b) c..Or]Ee'C“g ; c5j = -l c ?c % w.(a.a.a) ; c. finite '. 2 -- j-O V2(€aa) = _ (15C) g;_e-Cu€ 1n [1-016 ng . C = m C 'ng ’ 5 2 1-oe r j ('+1)a c a Cl “C2 J _ " u . ° _ C6j+1 [9 9 J 9 C63 f 1} V3(€.J.G1 = < g (15d) [o3(c.-c2)Ee'c“€ ; cej = -l J mm = b2(je.)e'°r’=-K.c.aneowaaneon (ise) 13 ( - __113:11-F Z 0W5( .i a) ; K2? 1 Wu15sa) = < (15f) 'C15 C15 ,- Le 1n §___:££% ; K2 = 1 C1 1 + Kaa I-K3a|£ (2+1) )C1€_ e-Cug] ; C72 f -1 _ c7 t+ l WS19329a) - (159) L'Kaa]Q K1(K2'1)€e-C“€ ; C72 = -1 -C K3f:s -f m d1(g)=d1(n60)e 5' ‘—E:—i' (166) 2:2 QQC -C 2' -E WG(€OQOQ) = E:E%T [e 2 9_e J (16b) 1/C1oo : = - -, -€_ C1 'Q G M.) (image—(13,961 2- foaming 2 amen] (16c) _ _- 2+2-1/C1 _ _ W7(€,Q,Q) = E+%K:]1/C1 [8 (£+2)CI€'e a] = . = ___. = c2 , = C u 1 + K2. Cs Cu . Ce C2-Cu . C7 'ETIfjfg) (159) System Behavior (16d) The qualitative behavior of f is identical to the first order case, therefore the solution to f will oscillate regularly, and approach asymptotic limits either from above or below, except in the case where the asymptotic limit is established immediately. The proof of the limit cycles for 6(6) and d(e) rest on phase plane arguments. Since the form Of equations (l3c) and (l3d) are similar, the phase plane behavior of b(f) and d(f) are also similar. 14 The differential equations governing the behavior in the phase plane of b(f) and d(f) are Obtained by dividing equations (l3e) and (l3d) by either (l3a) or (l3b), depending on which region is being considered. Points in the phase plane where gg-Or gg-become zero are found by setting the numerators of these expressions to zero. This is equivalent to solving for the locus of steady state solutions to (l3c) and (l3d). The loci of zero slope in the phase plane (loci of steady state values) are given by KlfS 2 b5 = (__W'TKZ , d5 = K3fs (17) For the purpose of the proof that follows it is sufficient to note that both loci are monotone increasing and pass through the origin. Hence. the proofs for the two cases are similar. The four differential equations in the phase plane can be represented by dgj _ F(f11'9i d—f.‘ ‘ ___.”. H” 1’2 (‘8) 1 B__'_|_ dO where 9 represents either b or d, F(f) is given by equations (l7), and B is a positive constant. The numerators of equations (18) will be positive for g < F(f) and negative for g > F(f). Since we know the signs of 351 after oscillation of f has begun, the signs of (l8) are known Region l (l9a) 15 1? < 0 3 92 < F(fz) 0. f2 Region 2 (19b) 0. T\) d > O 5 92 > F(fZ) —h 2 Also, from equation (18) and equations (l3a,b), it follows that > O , i = l 91” (19c.d) /\ O o _a. II N since F(fi) is monotone increasing. The behavior of 9 will depend on the behavior of f, hence the limit cycle for 9 cannot precede that of f. Moreover, after f has reached its limit cycle, 9 need not have yet achieved its limit cycle. At f(ew), 9 could be in any Of three regions of the phase plane, viz., gF(T2(Jeo)m) 5 F2 and 9€[F1.F2]- If 9>F21 I. 2.22.20. 0 -C 6 - . c“ 'e A 0) ' K1C10(Jeo)wW17 (24C) W13 18 r 00 OCHOC ” 2015(j) ; C5 finite 1. u 2 j=0 = 1’“ c2 c2 d(neolwu _ C151HL1-0(n90)gJ‘ = m L 50 C15 c2 C6 1 f . . - d(n002i ltc13+]) _ l-e c160 ; csj f _1 1 CsJ+1 1+1 c2 c» W16(j) = > d(nQ )i (c -c;) ec“e° -c e . ”° “ -eoe “ ° ; C6.) = -l K Cu eCu / - ,5 . '1 f 0. n0 j A. 1 :23- ; J f 0 ‘l e-c‘060 \ ( 901w - . . _ C5341” . ’ C1, 9 C53 f 1 W15(1) = < b 60 , J - 0 _ > - ‘Cueo 01(neo)J [L - 60 e'c“e°] ; csj = -l K m c. 2 f \ _—1?2:11-F2 0W1e( 2) ; K2 2 l W17 =fi > %%i In 1+K3&(ng)w 1 1n l+K3&(jeo)mc;. ; K2=1 kl 1/C13+K3(1(j90)oo CTK301jeO)m 1+K3a(j80)co J K - . 2 (2+1) “C060 1 [Lfiflyfififllel_ 1'Clu _ 1'9 . _ cyt + l (2+l1c1 c. ’ c7£ # 1 W18(£) =< 'Cueo > [-K3&(j60) 1R K1(K2-Il [1- 9 e c161; C72 ' '1 w Cu Cu F J a = W13 + W29 200 K3f2 _ 9 = d110901m11‘e-60) ' "‘EEEL <§2[1‘90‘e e01+2(1+C3) (1+C312 % (1‘11W21( £=2 £1 Q=1 020W211i) + (24d) (24e) (24f) (24g) (24h) (25a) (25b) 2 1 i-e2C290 , -6 W21(£) = 0(neo)mC2 RCZ-l RC2 T e 0'] (25C) cl['a(j90)w]]/Cl KgC1-])/C1 wzo = d2(j50)m(l'e-60) - §(m+i)u,2(m) (25d) m 0 (m+2)c1+l ' - - " +2)C160 _ C1 ’0 90 .K3] C1 1-9 (m '60_ W22(m) - ——£——&%fi—%§7L1 _ 1 [ (m +23c1 + e l (25e) _ 2K3fl(n€Q)m + Cl’Cz 0(ne°)m ' 2K3f1(neolm + c; + c2 (25f) and amen, = 'f2(3'“°)°“ (259) K3f2(j60)m + C1 CALCULATIVE RESULTS All of the solutions presented are exact analytic solutions. Unfortunately, the complexity of the results has frustrated all efforts to draw a patent conclusions. This reduces the analysis of the results to the examination of the trends calculated from the analytic expressions. Numerical simulation of the process was performed as a check. The analytic predictions match the numerical simulation with deviations on the order of the accuracy of the simulation technique using the IMSL subroutine DVOGER. Table l shows the predicted enhancement of the time average outlet concentrations over the corresponding steady state outlet concentrations. As mentioned in the section on the method of comparison, the comparisons for f and d are made with the periodic residence time equal to that of the steady state (nonperiodic) process, while b is compared with the optimum steady state bS [8]. 20 TABLE l. Enhancement of time average periodic concentrations over steady state concentrations. Approximation to consecutive reactions when K320. Approximation to second order reaction when K120, K220. 21 m.w c.e ow. o.- o.o ——. o.po v.“ llllIll ll 2 o.cmn .o— o.-- .— oc. u p. ~.opn .o— —.¢ . .— 00. u —. m.m u OF N@. i .— Noo. u —. o .o— o .— o _. Koo. .op v—o. .p moo. —. Koo. .c. v—o. .— Noo. —. u 6: mmu\~mmu-muv .o— .o— .o— 08 08 Q ( O). O), o? .— .— .— .o_ .o. .o— 1.1.." -| it; - 'll] .o_ .o_ .op .o— .— .op .o_ .o— p. .o— .op .o— .op .op .o_ .— ._ .o— .o— .— _. .op .o— ._ .op .— .op .o_ .p .p .o— .op _. .F .o— .op .op .o_ .p .p .— .cp .F .p p. .o— .— .p .o— .— .op .p .— .— .op .p —. .— .o— .— .o_ .— .— .o— .p .— .p .o— p. .p .p .o_ .o_ .— .— .p .F .— .p .— —. .p .— .— 3. FM «x ~¥ .— w_nmh 22 Several trends can be distinguished from the data. Except for very small deviations in the case of consecutive reactions, conversion is enhanced for all of the constants examined. Frustrating the objective of the study is the result that periodic operation appears to be always detrimental to the yield of intermediate b. The only encouraging result is that the time average concentration of d is enhanced for all of the constants examined. The drop in f ranges from O-Zl.5%, while b drops 0- 84.3%. d, however, is enhanced by as much as 95% over the steady state. d is most enhanced when the kinetic constants favor the formation of b and c. This is logical, since the rate of formation of d is second order in f, while the formation of b is first order in f. By increasing the feed concentration, the amount of f in the reactor is increased. The rate of reaction of d grows as the square of that for b, hence production of d should rise. If, however, the kinetic constants_already favor the formation of d, then mostly d would be formed without cycling, making the possible enhancement small. The effect of cycling on enhancement increases as the period 60 increases. The influence of large 60 is shown in Figure 2. Notice the concentrations spend a large fraction of their time close to the respective steady state for the region. This configuration approaches a single CSTR with feed concentration equal to twice that of the steady state process being compared, operated for one half the total length of time. The approximation to consecutive kinetics where K3m0 has two dis- tinctive characteristics. First, there is virtually no enhancement in f (Table l). The apparent loss of production of b comes about because this comparison is made against the optimum b5. When compared to the steady state with the same residence time as the periodic process, the 23 periodic time averages closely match bs' Second, there is no change in enhancement (b or f) as 60 is varied. Both of these trends should be expected, since the consecutive system is completely linear. Figures l, 2, and 3 show the system behavior for 3 different combinations of kinetic constants and cycling frequencies. Figure l is a case of intermediate cycling frequency with 60 = l. Figure 2 shows a system which gives 5l% enhancement of d. This is slow cycling, with 60 = l0. As mentioned, the concentrations are near their respective steady states most of the time. It is also interesting that the trajectory in the b-f phase plane need not stay close to the locus of steady states. This behavior is expected for quasi-steady state operation where the change in the process input occurs so slowly that the system always stays near the steady state corresponding to the inlet condition. The rapid switching of feed concentration, however, makes this case different from quasi-steady state. Figure 3 shows a case of rapid cycling with 60 = .3. The concentrations are rapidly changing most of the time, and many cycles are required for the system to approach a limit cycle. The magnitudes of the oscillations are small compared to slower cycling. CONCLUSIONS The analytic expressions for time average outlet concentrations resulting from square wave cycling of feed concentration to an isothermal CSTR with Van de Vusse kinetics are apparently too complicated to yield qualitative a priori predictions. Analysis of trends resulting from various system constants shows that time average concentrations of f and b appear always less than the corresponding steady state concentrations. Time average of f compares with the steady state process having the same residence time as the periodic process while the time average of b is 24 FIGURE l.(a-d) Behavior of the Van de Vusse system for a case of intermediate cycling frequency. K1 = K2 = K3 = l, 90:]. 25 ov.w om.m ow.v oo.v P wmou>o :32 ow.mw ov.N P ow.~ ow.o oo. o va._ mmszu 08'0 09' OV'U UZ'O DO'U I 00' 26 mmou>u :32 omum Aav._ mm=wHu r 91'0 80'0 DU' VZ'D O 28' OV'U 27 mmgu>o :32 om.mw Auv.P mmzufiu ov.m _ ow. # om.o oo. oo-d:3 4T 91'0 80‘0 r VZ‘U O 28' ov- 28 om. fl ov. om. fi 00. H AuV.F ”gamed om. o ov. ON. 0 OD. [Jo-(f3 80'0 VZ'O SI 0 ZS' OV‘O 29 FIGURE 2.(a-d) Behavior of the Van de Vusse system for a case of slow cycling frequency. K, = K2 = K3 = l, 80 = l0. 30 DD. V om. m 00. m om.N F mmoo>o 2:2 oo.m om.fi co.“ om.o oo. P — AmV.N mmonu 00-0:3 0 r r r 08'0 09‘0 ov~o OZ‘ OO'I 3l 00. v om. m 00. m om.N r wwou»u :32 H:U.N om.~ L Do. H om. o co. Anv.m mmawHu 00' 80'0 QI'O tZ‘D 0 28' UV' 32 wm4u>o :32 oo.¢ om.m oo.m om.m oo.m om.~ oo.~ om.o oo.pu F _ P _ _l _ Pl b o c \ e o 0 foo mm Tmu Z 0 f0 p. 0 ID ..V 0 r0 ms nU AUV.N mmstu 33 00.~ 0¢.~ Ac0.~ mmszu 00.0 OD'D 80'0 F0 r ZS'O VZ'U 91' OV'D 34 FIGURE 3.(a-d) Behavior of the Van de Vusse system for a case of rapid cycling frequency. K1 = K2 = K3 = l, 60 = .3. 35 mmqu>0 :32 oo.mw Amv.m mmaaHa 00.v 00.N 00. oz-o‘3 OO‘I 08'0 09'0 UV'O OZ'I 36 0000xo :32 nzuum Anv.m mmsweu oo~ocD OZ'O SI'O OI’O SO'O SZ'U 37 OO.©H p 00.xfl OO.NH _ 00.3H 0000>0 :02 00mm Auv.m mmsuHa 00.0 P 00.v 00.N 00. oz-o SIJO otjo 90-0 00-0“3 SZ'O 38 00. H 00. 0 00. 0 05. Hc0.m mmszu 00. 0 0N. SO'E) [NJ-d3 01 V0 31' O 02' SZ'O 39 always compared with the optimum steady state. Time average concentra- tion of d appears always greater than the corresponding steady state with residence time equal to the periodic process. The subcase of consecutive kinetics shows no change under periodic operation. The subcase of second order kinetics has its product distribution shifted toward d and away from f. Hence square wave cycling of feed concentration appears detrimental whenever b is desired, and the optimum steady state process should be used. When d is desired, square wave cycling of feed concentration can give enhancement close to l00% over steady state operation. Periodic operation of porous catalytic systems is only beginning to be explored [9] [10] [ll] [l2]. The interaction of diffusion within a catalyst pellet with cycled feed concentration to the bathing medium adds a new dimension to the problems so far considered. The coupling of the mass balances for the medium and the catalyst pellet, combined with a complex kinetics makes this a challenging problem, and it is currently being studied by these authors. CHAPTER 2 PERIODIC SINUSOIDAL FEED CONCENTRATION TO A HOMOGENEOUS ISOTHERMAL CSTR WITH A VAN OE VUSSE REACTION NETWORK 4o 41 Chemical reactors operated under forced periodic operation sometimes give time average outlet concentrations which differ from the corresponding steady state concentrations. This is generally observed when the governing mass balances contain a non-linear term, either from non-linear kinetics, or, for non-isothermal systems, from the exponential temperature dependence of one or more kinetic rate constants. A good deal of theoretical research in this area has focused on kinetic schemes in which conversion is the only economically important criterion, probably because this approach avoids the complicated expressions which result from more complex kinetics. Con- version can usually be increased, however, by simply using a larger reactor at steady state, and this provides a much simpler remedy. A potential application of periodic reactor operation is in complex kinetics for which the yield of a desired intermediate product exhibits an absolute steady state maximum. A simple example is the consecutive reaction A + B + C. More difficult problems, which are considered here, are the kinetics scheme proposed by Van de Vusse A + B + C and 2A + D, and the consecutive-parallel reactions A + B + C and A + D, in both of which the desired product is the intermediate B. Douglas and coworkers have done a substantial amount of work on isothermal systems in this area. Douglas and Rippin [l] considered a sinusoidal feed concentration input to an isothermal CSTR with a second order, irreversible reaction. Using numerical simulation, they were able to show a very small increase in conversion for small amplitude variations in feed composition. Douglas [2] considered the same problem using a perturbation approach in which he introduced an arbitrary perturbation parameter. He was able to show good agreement of his analytical and numerical work, reinforcing the conclusions in [l]. Dorawala and Douglas [3] considered both parallel and consecutive reactions, using a sinusoidal flow 42 rate disturbance. They were able to show small (.3%, .02%) improvements in yield over the steady state system. Other investigators [4,5,6,7,lO] have looked at different aspects of forced periodic operation of reactors. Bailey [8] provides a review of work through 1973. We have presented [l3] analysis of the Van de Vusse system taking place in an isothermal CSTR under square wave cycling of the feed concentration. The results presented there were exact analytic predictions of the time average outlet concentrations which, unfortunately, frustrated all efforts to draw a priori conclusions as to the direction of change. In this work, using a very accurate approximation, we are able to draw a priori conclusions as to the direction of change of the time average concentrations from the steady state concentrations. This is most valuable, since it gives the relative effectiveness of cycling independent of the system constants, without resorting to numerical calculations. THE METHOD OF COMPARISON In order to make a meaningful comparison of steady state operation to periodic operation, several conditions must be observed. (1) Residence time: For those products for which steady state operation can be made as effective as desired by simply using a longer residence time (for example, the reactant A in both the consecutive- parallel and Van de Vusse schemes), the comparison is based on the same residence time for both periodic and steady state operation. For those products which exhibit an absolute maximum concentration with respect to residence time in steady state mode, the periodic residence time which produces the largest time average concentration is used. Any increase over the steady state maximum, regardless of residence time, is clearly beneficial. 43 (2) Equal time average molar flow rates of reactant species into the reactor must be used in periodic and steady state mode. (3) Initial transients are ignored. For the steady state mode, only the steady state concentrations are considered. For the periodic mode, time averages are taken only after the effluent concentrations achieve a stable repetitive cycle (limit cycle). (4) Volumetric flow rate does not vary with time in the periodic process. CONSECUTIVE PARALLEL REACTIONS Basic Equations For consecutive parallel reactions taking place in an isothermal CSTR, the transient mass balances for A and B are: - ea qu - qA - (k1+k3)AV - v dt (la) qB - qB - (k B-k A)V = v 5‘3 (lb) f 2 1 dt Here we take the feed concentration of B to be zero. The disturbance of the inlet concentration of reactant is given by Af = A0 [l + e sin(wt)] (2) where c is a parameter adjustable between zero and one and w is the frequency of the input disturbance. Introducing the dimensionless variables: f 2—,b ,AEwT (3) O , 6 III .o|< Ill rile-v- —B— T A0 and substituting equation (2) gives for equations (l) %%.+ [1 + 1(k1+k3)]f = l + a sin AO (3a) %%-+ (l + Tkz)b = 1k1f (3b) with initial conditions 6 = O: f = f(O); b = 0 (3C) 44 Analytic Solutions The solutions to (3) are easily found to be f(9) = _l__+ c{a2251n AB - A cos A6} + -60 “ 4 622 agz + A2 a31 eXP ( .25) ( ) _ l E _ 2 ~ _ b(O) ' Tkl [311312 + (a§2+A2)(a711+A7) {5611322 A ) Sin A6 3316 EXP ("3115); 611:31 A(a,,+a22) cos A6 + ( ) +a32exp(-a116) (5) a3leXP “3229 . , a a a11_a12 11¢ 1 where =1 + (k +k )' =1 + k - a =f(O) - 1 + E“ 322 - T 1 3 a all - T 2’ 31 T 322 532+AY (6a) ( ) 0 B 311 = 312 _ 1 EA 611+612 + C . 632 ' - Tk1 311312 + (a32+f~.2)(a¥1+/F) __311‘322 9 all #312 (6b) The solutions are periodic with the same frequency as the input disturbance. Evaluation of Enhancement The time average limit cycle concentrations are found by taking the limits of (4) and (5) as 6 becomes large, integrating over one cycle, and dividing by the cycle time A/Zn. The terms containing the constants in both equations vanish as 9 becomes large, and the contributions from the integrals of the sin and cos terms are zero, giving the time average reactant concentrations .._1_ f - azz (7a) _ Ik, 5 - 311322 (7b) The right hand sides of equations (7) are just the steady state values of f and b corresponding to e = 0. Hence we see that for consecutive-parallel reactions, sinusoidal forcing of the inlet concentration gives time average outlet concentrations equal to the corresponding steady state values. 45 VAN DE VUSSE KINETICS Basic Equations The reaction scheme k1 k2 k3 A + B + C , A + A + D was proposed by Van de Vusse [9] as a simple example that defied general reactor selection guidelines. The transient mass balances for this kinetics in an isothermal CSTR are v 31—? = qu - qA - k1VA - k3VA2 (8a) dB , V aE'- QBf - QB ' k2VB + k1VA (8b) v git:- = qcf - qC + szV (8c) v g? = qof - qD + kgAZV (8d) Imposing the condition in (2) and using the dimensionless variables in (3) and bf = cf = df = O we obtain gg+cfi+$zlf2=i+esinxxe (9a) $5- + ctb = mf (9b) {1—2 + c = tkzb (9c) g—EU d = Ik3A0f2 (9d) where C1 5 (1+Tk1): C2 5 2Tk3Ao, C3 5 [CI+2C2]]/2: Cu 5 (1+Tk2) (10) with initial conditions 6 = O: f = f(O), b = O, c = 0, d = 0 (ll) 46 Perturbation Solution An analytic solution to (9a) is apparently not possible. In order to develop an approximate one, the solution is assumed to be com- posed of a non-periodic term and a sum of periodic terms: f(6) = F1(9) + 6P1(6) + €2P2(9) (l2) This imposes the restriction on (9a) that e be small enough so that term on the order of c3 are negligible. Substituting (l2) into (9a) and collecting terms of like order in 5 gives the following three differential equations: dF 351 + {c1 + czrlip, = sin no (l3b) 3:2 + {C1 + C2F1}P2 = 'P12 (13C) Note that the zero order equation returns the transient mass balance for an unperturbed (e=O) feed. The solutions to equations (l3) are straightforward but are fairly tedious. Defining for convenience the steady state values corresponding to e = O: k f Tzk k f = (cg-c1) = T 1 s = 1 2 s = Tk3A0f2 fS c2 , bS Cu , cS -——7::———, dS 5 (l4) we obtain [ll] for the solutions to equations [l3] = 2c3qe'C36 F1(6) f5 + c2(]-ae_c3e) (15a) P1(9) (l 1-c36)2{W1(9)'20¢2(9)+02W3(9)+1le-ae} (15b) -ae P= ‘ {1 ‘ae ( + 2 (l-ae-c39)2 29 ’Wu 6) Ws(9) ‘ W6(9)} (15C) 47 a sin A6 - A cos AB WAG) = a2 + A2 (15d) W2(9) = e-C39 (a-%;)C:;g ABA; A cos Ag (l5e) W3(6) = e-2c36[(a-%§3%3§;n+Aiz- A cos A6] (15f) w,(e) = ?O(i+1)a‘{w7(i.e)+wa(i.e)-w9(i.e)} (159) 1: 5 2+] A 9,1 'C3(£+1-1)9 w7(1,e) - zg-l) [B.** ZA ][AZ +4A2 ][%i“Ae(A2 i sin AB - Q: , a 2 - 3A . ; A2,1 # 0 ) 2A cos A6) + 1’1 ..2. ‘39. - in We ,Az’i - 0 — (15h) . 5 2+] e-c3(t+i+l)e we(1.6) = Z](-l) YR (YR .+4Az cosAB(A; 1 cos A6 + 2A sinAB) 2= ,1 ” - 2 + 2,1 ’ (l5i) k 2A266-66,A£ i = 0 -c (2+i-l)e 5 ”29 3 . 2 . ¢,(i,e) = Z (_])R+l 2A Sln AB (lSJ) £=l w5(e) = 211% (i+l)oi{ E (-1)£+] ?;[a+ca§f:igl)]ilo(i.2.e)> (l5k) i=0 £=l 2,1‘3 w10(i,2,6) v [(AQ i-a)sinAe-AcosAB]-oi[(kfi i-a)cosAe+AsinAe] (l5l) It W6(e) = IIiEO(l+l)a '3F77TSET (l5m) and c f O +c -c . = ciféog+ci+ci 9 a=C39 AQ,1 = C3(2-Q-1) (lSn) 48 e, = 5%, 82 = 2515,, B, = 25,6, + 5%, e, = 2516,, 85 = 6% (15p) Y1 = 55. Y2 = 25252, Y3 = 25252 + 5%, yt = 25252, Y5 = 6% (lSq) O1 = 25152: D2 = 2(5152 + 5251), D3 = 2(5152 + 5152 + 5152) , U4 = 2(5152 + 5251) : D5 = 25152 (15?) v1 = 51, v2 = 51, v3 = El (155) 01 = 52, 02 = 52, 03 = 52 (l5t) 51 = 80;, 52 = A01 (ISU) 51 = O; 52 = ZoozA (l5v) E; = -o2o3c3, Ez = 0203A (15w) a; = a3 = (az+A2)-1 , o; = l/A2 (15x) Evaluation of the constants 11 and 12 requires the use of the initial condition as well as the differential equations. Since F1(O) = f(O) for all e, we have P1(O) + cP2(O) = 0. Setting a = O in equations (13b), (l3c) gives relations between PI(O), P§(O) and P1(O), P2(O). Eliminating Pf(O), P§(0), and P2(O) gives P1(O) = 0, from which P2(O) = 0. Since P1 only contains 11, and P2 only contains 12, their values are readily found to be 11 = 20Wz(0)'W1(0)'02W3(0) (ASY) 12 = w,(0) - ws(0) + w6(0) (152) 49 The Solutions to b, c, and d In solving the differential equations for b, c, and d, no assumptions about the form of the solutions were made. For b, the expression for f derived above was substituted directly into the differential equation and solved. Because of the complexity of the square of the full solution for f, the full transient solution for d was not obtained. The limiting d equation (6+w) was obtained by first taking the limit of f as 6+w, squaring the limiting form, and using that form in the differential equation for d. The b equation was used directly in the differential equation to find the solution to c. As a consistency check, the limiting c equation was evaluated by the limiting method used to find the d equation, and this result was compared with the limit as e+w of the full c solution. As expected, the two forms showed perfect agreement. The full transient b solution, which in comparison makes the solution of f look.easy, turns out to be of the form mm=1nwxm+eHW)+aAW)-€“eui (m) . . fs 2c wnere 61(6) = E:'+ “Ej‘ W11(e) (173) .% GJ+1 W12(9.j) 3 C3 f Cu 3'0 (l7b) a - O zz-e C3 ln[ec3e-o] ; c3 = c, W11(9) C5 W12(9:j (17C) L e(C5'C4)e ; C5 f 0 ) = 'Cue Be ; c5 = 0 C5 = Cu ' (5+1)C3 (17d) 50 and ( ) m . 3 C3 1+1“) 6 P3(9) = Z (i+])01‘z '])£+1 w: + A2 W13(isiae) +11W1u(1s9) (185) 1:0 :1 2:1 w13(2.i,e) = (u£w£,i-OQA)51nAO+(Ap£+fi£w 2 1)cos AB (l8b) e-(cai+a)6 . Cu'Cai-a ; c“ f C3] + a Wlu(is6) = (18C) ee'c“e ; c, = c3i + a w£,i = Cu-C3(£+i']) (18d) U1 = 51, U2 = 519 U3 = 61 (189) and E1 = 52: R2 = 529 R3 = 62 (18f) _ w k , , w e . i+k . Pt(6)-12£a (k+])W15(k96)+Z a(1+l)(k+l)a w1e(k.i,6) (l9a) k=O k=O i=0 where e-(a+c3k)e Cu’a'C3k , Cu f a+C3k W1s(k,9) = (19b) ee-c,6 ; ct = a + c3k 5 3 w16(k91se) )= Z (])£+]{—w17('k, 2,],9)'W18(k,2,l,9)}+211Z(-1)2+1 M= 2:1 WI9(k!£!ise) ' I?W20(k9iae) (19C) re[C1(k,2,i)-c,]e . . , CTIKfiai)z+4AZ W21(ka2a1,6); c10s2.1)#o A I r (19d) ‘Cue 6 sin 2A9 . e [-- -——————— ; k,2,i = 0 K 2 4A C1( ) 2 w21(k .£,i,e) =C-%§:§—;)+ sinAB[c1(k,£,i)sinAO-2AcosAe] (l9e) W17(k,£,i,6) wle(k,2,i,e) =< SI e CS(£91) “C49 (C5(R,l) W2 W2 reukan-coi k , 1,9 1 k,£,i +4A 7(£,1,6) ; A£,l=0 T W25(k,2,l,6)];g1( (k Hi 1 )f0 T (21(k92’31) O W22(k92al:5)=Cu(£,l)W21(k,£,l,9)+Y£A E MW23(k 1,1 9)+¢2u(k,£,1 6) (199) 2A2 W23(ka£,l,5)= :‘TETET‘ w2q(ka£9196)=hg2( $25(k ,9. ,i,6) = 2(12C3(£,1)J w29(k,2,i,e) ; A£,1 = 0 ‘ J , e-ae 1 . iT(k,£,i)-a [6' ;1(k,g,1)-a ; C1(k,1,l) f a 1 WZ8(ks£9i,e) = 2 . 9"9-§1(k’£’1)e ; C1(k,£,i) = a K 9 sin 2A8 sin 2A9 w26(2.i,e) = §“(£’i)[§" 4A 1 + Y2A£,1[§‘ 4A.‘“J 2 --£-L--2 cos 2A8 + 2“ ;3(g ,1)e 2,1 . 2 , _ 427(A,l,9) = éaéflfll) cos 2A8 - 2A a 2 1) e a6[e + g] $19(k9Q91196) = w29(k9£9126) W30(k:Qai99) )+cosAe[§1(k, ,t, i)cosAe+2A sin A6] £,i)[c1(k,£,i) sin ZAB-ZA cos 2A6] [ l ' A At i§1(k,£,i) 2,1 l f 0 e-[a+ca(k+£+i-])]e[@29(ks£9136)'w30(k,£9196)] (11 i-a)? + A2 E1(k,2,i)-a]7 + A2 [v 2.93m 'a) ' ORA] [(C1(k9£:i) -a) sinAe-AcosAe] = [V£A+o£(A£,1-a)] [(;,(k.2.i)-a)cosAe+AsinAe] (19h) (l9i) (191) (19k) (191) (19m) (19n) (l9p) (l9q) (l9f) 52 e-[a+c3(i+k)le ‘ (a+c.u[c.-a-c.(i+k)1 ; ct # a + c3(i+k) wzo(k.i.e) - > (l9r) 6e-C“e . 13—1-E377' ; C» = a + C3(l+k) J and c1(k,£,i) = cu-Ca(k+fi+i-l) (195) n}. . mm) = YE-Bi- ‘Lz—K’i (l9t) n A . 43H.” = Y9. +8. + 3;" (1%) n A . ;u(9.,l) = Aid. (6; + ifil) (19V) c5(£,i) = A: i + 4A2 (19w) The double summation appearing in equation (193) is a condensation of a quadruple summation which occurs when the term l/(l-oe’cBe) is written in series form and squared. The constant of integration 13 was found by setting a = 0 and b(6=0) = 0. The result is simply 13 = 31(0) + €P3(0) + €2Pu(0) (20) The simplicity of equation (20) belies the triple summations which it contains. The full transient solution for C(B) was evaluated by substitution of equations (16)-(20) into (9c) and integrating. The result is almost twice as lengthy as b(e). It was used here as a check to the limiting form (6+w) of c(e) found by using the limiting form of b(e) in the differential equation. Because the full transient c equation is almost twice as lengthy as b(e), and because we are largely interested in the limiting behavior of the components, it is not included here. 53 Limit Cycle Forms From the perspective of analyzing process performance, we are interested in looking at the behavior of the system over a long time period. It can be seen from equations (lS)-(l9) that both f(e) and b(e) consist of sin and cos terms, many of which are damped out as 6+w (this is easier to see when the equations are written on one line; however, this requires a sheet of paper almost l2 feet long). As e+w all of the exponentials become zero or one, and the only 6 dependence remaining arises from the trigonometric functions. The functions oscillate regularly, and both functions complete one cycle in the longest period of their component functions, Zw/A. Moreover, since C(8) and d(e) are obtained by integrating an exponential times either f2(e) or b(e), both C(8) and d(e) will retain the cyclic properties of f(e) and b(e). As 9+m, many of the terms in f(e) and b(e) drop out, and we obtain the simpler forms: mew) = f, (2la) P1(8w) = a Slna99+-AQ COS A9 (2lb) 81 T fliéigg. , . 2A2 p2(ew) _ - Af,o+ 4 AY SinAe(A1,o SinAe-ZAcosAe)+ A1.o _ 1 . 2A2 nlsinZAG XETETZK7' cosA6(A1,ocosAe+2ASinAe)+ A1.o+ 2A (21c) b5 31(em) = $13- (22a) p3(em) = (“lwii°'P’Ai§::Af'£9“L+plwl’°) C°SA9 (22b) 2A2 0 2 o Infar- 4’ c..51n AB-ASlHZAGJ 9 m - 2A(CE+4A2) - A1so+4A (22C) .0 5' A CD v ' W33(9) = £2( 54 W31(9) = X1 + W33(9) c? + 4A2 The limiting forms of the c and d equations are: where and where C(eco W37 W38 Wul ’ ) = C5 + T2k1k2EEP5(9m) + €2P6(9m)] _ p3. sinAe-was cosAe (Hie TA?) (HA7) V I V I ' Was + $37 siner + Was sin 2A6 = d5 + Ik3AO [6P7(€m) e2P3(6m)] 2fs[(a-A2)sinAe-A(l+a)cosAe] w) ' (a2+A2) (l+A7) . 2 . - Wulsln AG’qu Sln2A6 + Wug 8V I = (U1w1,o ‘fi1A) ’ A(AU1+fi1uHm) = A(ulwlfi-Q1A) + (AUI + filwlfl) 2A ’ X3 + (T:E377'[A(X~+Xz) + X5] = l 4A2 T35X7 [Xv-4AX5- (7§§%:ZKT)J l 4A3x 1+4KY’[(A§,,+ZAY7' AX“ ’ X5] 1g2_A2)_4aA2 - (a2+A?)2 + 2fs [Xus'4flwue] (l+4A3) (22d) l,0){A[(A10+ct)sin2Ae-2Ac052AGJ-ctklfisin2 A9} (22e) (23a) (23b) (23c) (23d) 55 az-A2 +a A LX7373K;771"+ 2fs[w“5A + w“6] (23m) w“2 = 1+4A2 A2‘};§i:§§§2*za) + 2fS[2A2w.s+2Aw..J was = (l+4A2) + 2f5w~u (23”) w““= (Ai:1+4A2) [ZAziifl’0) + Y‘A"°] (23p) Was = 2‘7: +W C2“ 0) (23(1) 1).. = fi—fifim‘r on 0) <23r) = 2"‘223(1’0)<3§AR2 + .11) 11-1—2— <23» x. = $21—59 (23t) X3 ‘ c.t’c‘214.-a ' #27:“! (23”) X“ "' Maggi? Ciliéifiil (23V) A(>\1:0+C~)X2 + “1 (23w) X5 = A1so+4A2 ETEE:ZKTT All four equations are limit cycles with longest period 2n/A. In each case the first order correction term (:1) is composed of first powers of sinAe and cosAe. The second order corrections (52), contain sinZAe, sinzxe and constants. Analysis of Time Average Concentrations The time average concentrations are found by integrating the concen- tration functions over one period and dividing by the length of the period. The longest period appearing in an equation determines its overall period, which here is 2n/A for all four equations. Notice that since the first 56 order terms are all composed of sines and cosines of A6, their integral over one period is zero. Hence the first order terms taken alone do not give a time average concentration which differs from the steady state concentration. Similarly, the sin2A6 terms which appear in the second order (52) terms have a zero integral over one cycle. However, the sin2A6 term yields a constant when integrated, and this combined with the constants appearing in the second order terms allow the time average concen- trations to take on a value different frOm the steady state. Considering the complexity of the full transient solutions, the values of the time average concentrations, with dedicated manipulation, can be put into remarkably simple forms. After integration and simplification, we obtain for the time averages: 2 _ E f ’ fs ' 2a(a2+A2) (24a) 5 = b TklAc2[(a2+A2)cEa2+2A“(a2+ c3+4A2)] (24b) 5 ' (AROMA?) (62+qu ' _ €2T2k1k2(Yt+81) c - cS - 2C“A1,0 (24c) 0.! ll 0.. + 2 _2_C'+ C -2 C3 S e on 251W}— (24d) lUl of the time average concentrations are less than the steady state except 6, which varies depending on the magnitude of the constants c1, c2, c3. Qualitatively, it appears that only when c1 and c2 are both small will a be less than ds' From these results, it is clear that if’B or C is the desired product, sinusoidal fluctuation of the feed concentration (when 53 is negligible) is detrimental to product formation. If D is desired, then for moderate 57 values of c1 and c2, the average concentration of D from the reactor may be greater under sinusoidal concentration perturbation than at steady state. The magnitude of any enhancement can be calculated from equation (24d). If B is desired, then our aim is to exceed the maximum steady state b. However, 5 is always less than bS at the same residence time. Since bS at any residence time is less than or equal to the maximum steady state bg we are guaranteed that the steady state process is more effective. The objective of enhancing the yield b in Van de Vusse kinetics by a small (53:0) sinusoidal perturbation of the feed can never be achieved. NUMERICAL RESULTS The heuristic nature of the perturbation solution to f(e) requires verification. Since equations (9) are not highly non-linear, they are easily integrated numerically with high accuracy. In this case the SOIU- tion to equations (9) was simulated using the IMSL subroutine DVOGER on MSU's Control Data Corporation Cyber l70 model 750 computer. The results are presented in terms of the dimensionless groups Tkl, kZ/kl, Tkng, A, and c. The initial conditions used for all calculations were f(0) = l and b(0)=c(0)=d(0)=0. In the formation of the perturbation solution, the assumption was made that terms on the order of 53 were negligible. However, for some kinetic constants the solution gives reasonable accuracy for 5 up to the physical maximum of c=l (since c>l requires the feed concentration to be negative). The equations used to evaluate the perturbation solu- tions were the full transient solutions for f and b, equations (l5), (l6)-(l9), and the limit cycle forms of c and d, equations (23). The perturbation solutions were found to show good agreement with the numerical simulation for the values Tkl = k2/k1 = TkaAo = A = 1. These values of 58 the constants will be used as a starting point to show how the deviation of the approximate solution from numerical simulation changes with the various system constants. Figures l and 2 show the results for the constants tk, = kz/k1= Ikng = A = l. For 6 = .l and e = l. Figures le and 2e show the b-f phase plane. As expected, the limit cycle is much larger for a larger perturbation. Virtually all of the calculations of concentration vs time showed behavior qualitatively similar to that in Figures 1 and 2. Figures 3 and 4 and parts e and f of Figures l and 2 show the percent deviation of the approximate solution from the numerical solution for 5 values of .l, .25, .75 and l. The error increases significantly with a, but for e=l remains within about 115% except for periodic spikes. The Spikes are caused by a relatively constant error when the value of the function itself gets small enough for the error to be significant. Considering the initial assumption that c3 is negligible, the agreement is quite surprising. Figures 5 and 6 show the sensitivity of the solution accuracy to increasing Tkng to 5. The agreement is still good for c - .25, but increases sharply for c = .75. Figures 7 and 8 show the error for Tk1 = k2/k1 = 5. Even for c = .75, the error remains small, showing a small sensitivity to these two constants. Figures 9 and lo show what happens when 1k; = k2/k1 = Tszo = 5. The error increases sharply even from 1k; = kz/kl = l; TkaAo = 5, which already showed a substantial error. Hence, although the solution accuracy always deteriorates as Tkng is increased, the effect of Tk] and k2/k1 may be large or small. Finally, the effect of increasing A is shown in Figures ll and l2. Decreasing A is shown in Figures l3 and l4. The error remains small 59 except for A = .2 and e = .75, where the d equation spikes up to almost twice the numerical simulation. Why this should occur for this slow cycling is not clear. Notice that in many cases, the initial response of the perturbation solution to f drops below the numerical simulation, giving a sharp down- ward spike in the error. The b solution initially tends to climb above the numerical, giving a sharp upward spike in the error. Exactly why this occurs is certainly not clear from the huge equations, but the perturbation solution apparently cannot follow the rapid changes in the initial region. It may be that higher order terms contain higher harmonics which are quickly damped, hence they could add (or subtract) an initial spurt to the solution before they die out. It may also be due in part to the initial inaccuracy in the numerical simulation, since a fixed step size is required to obtain the time average concentration by inte- gration of the numerical simulation. A typical plot over 3 cycles required about l.5 central processing seconds on the CDC l70 model 750. This involves about 90 evaluations of the b and f equations which totalled an average of about .l5 CP seconds, or l0% of the total computing time. A reasonable estimate is at least l.0 CP seconds were required for the numerical simulation, with the remainder used to evaluate c and d. Hence, the perturbation solution, while much more laborious to program, runs much more efficiently in spite of the triple summations that must be evaluated in each call. The evaluation of enhancement by numerical simulation requires that the differential equations be solved and then integrated over time, in this case an average of about I cp second for each set of constants. The predictions from the perturbation expansion, equations (24), are much simplen requiring a small fraction of the time and memory. Hence 60 FIGURE l.(a-g) Behavior of Van de Vusse system for Tkl = kz/kl = Tk3AO = A = 1, E = .1. FIGURE 1. pert 4.00 v 3‘50 CYCLES .00 NUM 3T .00 IE'O £330 Sl'G LC' of»; L) gfi (\‘I - St.) L) C) ( J x.) G" (.3 C‘ Q ( ) (_1 ‘1 ff -00 '1' C1. L L3 .50 rte 62 mu40»u :02 00.0. 09v oo.v om.n oo.m om.~ oo.m. P b b L b P 0 .0d 9..” 3 u .0 3 TDA 0 9N n u N 10A mus 0 d 3 HO 1 r0 3 ' A 07 v 00. wm00>u :32 m 00. 00.~ omw_ 00w_ 00.0 00 Aw mad ow 0 U .1: l: l (l( vmuA OI 0 N n Us 0 10A mxb 0 d 3 HO fi 1 0 r o H. 0 04.0 00. P... w 00. omuo 00w0 0c. c. . p C 0 to U 9 :2. a I... 8 ind m T T ' 0macwucou . F mmszm E 63 FIGURE 2.(a-g) Behavior of Van de Vusse system for Tkl = kz/kl = Tk3A0 = A = l, c = l. 64 mm00>0 :02 . . . 0000>u :02 oo.m on.. no». omwn sown omww oa.nu oo.n om.~ oouu owns on“. capo oo.pu m m I... v0 ..,. 0 8 no as! m r. :0 r. 30 w m 0040>0 :02 m000»0 :02 oopm cave cow. omwn oowm omww oo.m. 00. cm. cow“ emu. eon. oauo oo.pu E q 9 0 7.3 0 1. 1. 0 7 V 0 I1 .2 to; a I... -0 r 9 79 o “to. u TCO ADV T... AMV 79 0.- 1s 0 0 .N mmszm 65 00.0 «040,9 =32 . can. can. can" gown ooww oo.~. mu... mun a. 0 '3 r_40 0030 183d SA HUN A OO'OZ OO'OZ' 0000>0 :02 om._ o oown omun 00u~ .. 0w0 00.0 00.0. h z i mUd on“ . m .flnu o. n. “N n“ ”H 11A 0 we. 0 nYd a. M. fi I. z z w 0 m 0““.— 00 — 00.0 00 0 0w 0 0~.0 00.9U - P h b .P P o o O I! as... Lu Z 0 cmmnv an to r 0 E r. casewpcou - N mmaqu A 1.... V Amv FIGURE 3.(a-b) FIGURE 4.(a-b) 66 Behavior of Van de Vusse system for Tkl Tk3A0 = A = l, c = .75. Behavior of Van de Vusse system for 1k, Tkng = A = l, c = .25. kz/k1 kz/k1 FIGURE 4. FIGURE 3. .00 1150 'NUH CYCLES T 1.00 n U jo 00. .1 0931 a 183d SA HUN A (U V 3 2.00 2.50 1150 NUH CYCLES .00 r I 0.50 0.00 on it 007a 0039 00'!- 183d SA HUN 'A30 83d A (U V 67 0:”o o: d? ' DE 6- 05' 133d SA HUN A30 383d A J: V 1 ° > 0 5.00 .00 .50 3‘50 NUH CYCLES - 0 T 0 m 2.50 -00 0055 oo-f— oo-s- lHBd SA HUN A30 oc-s- 3333 FIGURE 5.(a-b) FIGURE 6.(a-b) 68 Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k2/k1 = l Ck3A°l005 A = l e = .25. % error = (Cpert - Cnum)/ num Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k2/k1 = l Ek3A01305 A = l c = .75. % error = (Cpert - cnum)/ num’ FIGURE 6. FIGURE 5. 1.50 NUH CYCLES 1.00 0.00 oo-dEVS 00-62 00-61 oo-ox— 183d SA HUN ‘A30 83d ('6 v 3.00 2.50 2.00 LSD NUH CYCLES I 1.00 Y 0-50 0.00 oo-r 9 003' 0030 00"- 1833 SA HUN 'ABU 33d A (U V 69 4.50 5.00 .00 4 1 3&0 NUH CYCLES .00 1 3 2.50 D O oo-éT—— 0039 0030 oo-e- 183d SA HUN A30 383d (M -00 4.50 4.00 I 300 Nun CYCLES .00 O 0 oc’I'S oe’o oe{o 02'0-N 183d SA HON A30 383d (b) FIGURE 7.(a-b) FIGURE 8.(a-b) 70 Deviation of perturbation solution from numerical simulation for Van de Vusse system rkl = kz/kl = 5 Ek3A01301 A = l e = .25. % error = (cpert - cnum)/ num° Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k2/k1 = 5 Ck3A0100] A = l e = .75. % error = (Cpert - Chum)/ num' FIGURE 8. FIGURE 7. 3.00 T 1.50 NU" CYCLES 1 1.00 0.50 .00 093210 09’: 1836 SA HGN A (U v 0 .0- V 09’ l o 09 A30 83d 2-50 3.00 .00 2 F L50 NUN CYCLES T l.00 71 I 00‘ A .D V s 00" oojo oo'v 183d SA HFN A90 333d V/, .00 4 1 3.50 NUM CYCLFS .00 Y c 0r? 9 153d SA kGN l’ , A 13 V )t‘ ‘LJ FIGURE 9.(a-b) FIGURE lO.(a-b) 72 Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k2/k1 = Ek3A01305 A = I e = .25. % error = (Cpert - Cnum)/ num Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k2/k1 = TkaA - 5 A = l e = .75. % error = (Cpert - Cnum o Cnum-lOO. )/ 73 0000»0 :02 0040»0 :02 oo.m om.. oo.. omun oo.n an.“ ca.mu oo.m om.. oo.. ovum on.” om.~ oo.mu pl b P b P 7: ._ 00.0 00.0 0040»0 :02 om._ co P b b b b P L 0000»0 :02 co." om.~ oa.~ om._ oo.. on.o oo.o P b P P F b 00.0 om.~ 00.N 00‘0 09’0 ‘A30 83d ‘A30 83d 02'0' oz-o 183d SA unN OO'O Y oo-e 183d SA unN 1 A (U V I r 00'91 A (U V 09'0 .0F umszm .0 mmame FIGURE ll.(a-b) FIGURE 12.(a-b) 74 Deviation of perturbation solution from numerical simulation for Van de Vusse system 1k; = kz/kl = Ck3A9I00] A = 5 c = .25. % error = (Cpert - Cnum)/ num ' Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = kz/kl = Ck3A0100] A = 5 c = .75. % error = (Cpert - Cnum)/ num mm..0>0 :02 mmfi§u :02 75 00. 0 om. h 00. 0 0m. 0 00. 0. 09.0 omwm 00.r0 00.1. 00.... 00% 00. 0. b r b r pound Mud mm” 00 3 3 . U i . nU _. in .. Bu 0 v n 0 H H . A 0A I| .l/ . U To 9 pl? 708 0 u 0 \ 3 a 3 3 8 8 I. I. rl r0 0 0 00.65 :02 mgowsw 202% 8.. o 8. o . ... 8.... 8 .Lo. 8. w... 8. on 8.... . OF 8% 8. o. D 0 o .4” my a DH 88 0 0 3 3 1.0A #TMUA o. . o 0 0 N 9N 0 0 H H 10A 10A rs .nus 0 0 d d 3 0 3 no 8 Vol. A3 -01 .b we 0 ' .NP 0000: . Z 0000: A .0 v FIGURE l3.(a-b) FIGURE l4.(a-b) 76 Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = k2/k1 = Ek3A9130l A = .2 e = .25. % error = (Cpert - Cnum num ' Deviation of perturbation solution from numerical simulation for Van de Vusse system Tkl = kz/kl = Ek3A01301 A = .2 c = .75. % error = (Cpert - Cnum num' 00. 77 00. m00u>u :02 ompn o m on. oo.. o.n om.~ oo.y b F F S Dd .3 0 .8 u 03 .U &u3. 0A 0 N m Q\ 19A 1C. o 08 1. .8 l .m 0 U 0 mm00»0 :02 % 00. Down 0mw. oow— omwo 00.0. m.” 3. 0 0 .1. .m.n NU. ON 0 H [alum \\ 7.9 ma a .3 .8 l .m PU EU .vF 020000 00. 00. 0000»0 :02 MW 00 v 00». owwn 00wm onwm 00.N z 0 0 / o 1. U 0 0 //I\\\\\lII/////l\\\\\llI/////l\\\\\ll///.2 A o 0 u l' q 0 0000>0 :02 m 00.0 00.~ 00.. 00.. 00.0 00.0 F b L p F P .P U» o 0 0 V 4a \ o a .t 0 O .H D 0 .mp mm0uHu ‘A30 838 1838 SA HGN 30 E 78 the perturbation solution gives an accurate and convenient prediction of the time average system behavior. Table 1 shows the predictions from the perturbation solutions for those values of the system constants which appear in the figures. The first 4 columns give the percent error in the average concentrations calculated from equations (24) when compared to the numerical simulation. The largest error is 5.43%, including all the cases where the error in the functions themselves were large. The second four columns give the predicted enhancement in concentration. Again, f, E, and 3 compare to the steady state with the same residence time as the periodic process while 5 compares to the optimum steady state b given by DeVera and Varma [l2]. Several trends can be discerned from the data in Table l. These are very similar to the behavior of the Van de Vusse system under square wave feed concentration [l3]. As with square wave cycling, the effect of periodic Operation on the time averages increases as the period of the feed disturbance increases, or as A decreases. This trend implies that the most effective cycling for this case is quasi-steady state cycling, where the process input changes very slowly. In that case, the process responds quickly enough to always be near the steady state corresponding to the instantaneous value of the input. If the effects of cycling are desirable, it may be more advantageous to run several steady state processes in parallel, since the performance of a quasi-steady state process can be approached in this manner. The form of the analytic expressions (24) show that conversion is always enhanced, yields of b and c are diminished, and yield of d is almost always enhanced. These a priori conclusions reinforce the same trends in the data for the square wave feed, although a priori predictions were not possible in that case. 79 TABLE 1. Percent deviation of approximate time average concentra- tions [equations (24)] from numerically simulated time average concentrations; percent enhancement of approximate time average concentrations over steady state concentra- tions. .mm_0 a macaw xcmmum E:E_uao ou moccasou n--umco_cma mm mama mocmcwmwc mEmm saw: macaw xvmmum ou mcwaeoo u .u .waa a to .u .9 .c - Pu. 8C) ¢.v_ wo.N- mm.m- ma.~- qu. ooF.- mv.m oo_.- ma. N. ._ ._ .. cc.F _mm.- w~a.- _mm.- coo. _oo.- Nam. _oo.- mm. N. ._ .F .P _m.m ~N~.- NNN.- NNN.- oso.- opo.- Noo.- _oo.- ma. .m ._ ._ .— oom. .mo.- .wo.- _wo.- NNc.- m_o.- Nooo.- _ooo.- mm. .m .p ._ ._ m._~ mc¢.- N.¢o- mv¢.- m~.~- mm.p- om.~- mm..- m“. ._ .m .m .m mm.m ovo.- m.eo- oqo.- Pmm.- wo..- o-.- wa_.- mm. ._ .m .m .m N.¢N woo.- m.mc- moo.- «No. mooo.- o...- wooo.- ms. ._ ._ .m .m e~.~ «No.- m.mc- «No.- moo. coco. _N_.- oooo.- mm. ._ ._ .m .m m.__ N¢N.- mm.~- No~.- om.m- mo.m- Nc.N- mo.m- wk. ._ .m ._ ._ _m._ mmo.- mem.- mmo.- Nao.- Nmm.- NN_.- 0mm.- mm. .F .m ._ .— o.- c~.e- mo.o- QN.Q- Na“. m_m.- N... m_N.- ._ ._ ._ ._ ._ m.N_ oo.~- ov.m- @o.~- 9cm. voo.- Ame. eco.- me. ._ ._ ._ ._ no._ N¢N.- mum.- Nm~.- eoo.- xcoo.- ch.- mcoo.- mm. ._ ._ .F ._ m-.- Neo.- _eo.- Neo.- Ncoo.- mcco.- m_o. o:co.- _. ._ ._ .P a u a c m w m w a < «.8 _x\.x _xp 00_ . m _0\Am .0-ucmm.wuv 00_ . Ea: _0\AE=: _0. atom muv - -zauemmmmamwnmmmwmmmi. .scaa._ a==._- tau 0 m> 0 .cowumw>m0 ucmocma .— m—nmh 81 One other trend is observed which also appeared in square wave cycling. The effects of periodic Operation were most pronounced when the kinetic constants favor the formation of 8 over the formation of D. This is in fact logical. If the kinetic constants favor the formation of D, then little 8 or C would be produced, making it impossible to sizeably shift the product distribution from B and C toward D. CONCLUSION Sinusoidal concentration fluctuations to an isothermal CSTR with either consecutive parallel or Van de Vusse kinetics is ineffective in increasing the yield of intermediate product 8, or the output of product c. For Van de Vusse kinetics, however, output of the side product D can be sizeably enhanced by up to 25% by such an input disturbance. The perturbation solution given above provides a convenient and calculationally simple method of approximating the time average output of the reactor. CHAPTER 3 PERIODIC SQUARE WAVE FEED CONCENTRATION TO AN ISOTHERMAL CSTR WITH A VAN DE VUSSE REACTION OCCURRING IN POROUS SLAB CATALYST PARTICLES 82 83 Periodic operation Of chemical reactors has its most significant impact in the area of heterogeneously catalyzed systems, which are by far most important from an industrial viewpoint. The complexity Of the mathematics evolving from the treatment Of such systems, however, is immense. For this reason, research in this area has either relied heavily on extreme simplifying assumptions, has used approximate or numerical approaches, or both. Bailey has explored the theoretical aspects of periodically operated catalytic processes. His early work emphasized a control oriented approach to the tOpic and neglected diffusional or mass transfer effects entirely [l,2,3,4,5,6]. His later work considers a less idealized case in which a porous catalyst pellet is examined under the influence of a periodic boundary condition [7,8]. He was able to show increases in selectivity using numerical and approximate analytic techniques. Others have studied similar problems from a theoretical viewpoint. Oruzheinkov et al. [9] used an approach similar to the early work of Bailey and concluded that selectivity in a non-porous catalyst could be improved by periodic Operation. Rice and coworkers [lO.ll] looked at the transient, non-periodic behavior of a porous catalyst both with and without external mass transfer limitations. Experimental work in the area has centered on problems in which diffusional effects are absent [12,l3,l4], with the exception of Unni, et al., who showed an increase in conversion from concentration cycling [15]. The Model Unlike the other models, e.g., Bailey, we consider here a case where the perturbation is caused by an external forcing function. A CSTR in which porous catalyst pellets are contained is fed with a stream whose 84 concentration varies in a square wave. The pellets are considered to be one-dimensional. A complex reaction network named after Van de Vusse A+A53O takes place on the active surface within the catalyst with no reaction volume change. The diffusivities Of the various species may be unequal but are considered constant. Mass transfer coefficients are assumed equal for all species. Adsorption-desorption effects are lumped into effective rate constants which give the rate Of surface reaction in terms Of the adjacent fluid phase concentrations. Diffusional resistance within the pellet is considered important and an external mass transfer resistance is included. Components A and B will be considered. The Basic Equations Two mass balances are involved for each component: the balance with the catalyst slab and the balance on the fluid medium in which the particles are bathed, which will be referred to as the bulk phase. For component A in the slab we have 32A 2_ g5 DeA §;7-- k1A - k3A - at (la) and the in the bulk phase “ .81: . . 9A (M q Af(t) qA DeAax (x L,t)Anp Ve dt with boundary conditions “ _ _ , _ _ , BA _ _ , 3A _ _ “ _ A(t-O)-O, A(x,t-0)-0, §;-(x—0)-0, D ——(x-L)-kq[A-A(x-L)] (lc) eAax 85 For component B in the slab, 82B DEB-377+k12=gA-I(B flIICD and in the bulk “ BB _ _ 51_ q8f- qB- DeB SE-(x—L,t)Anp - Ve dt and boundary conditions “_-. ==.8_B,=.§_f_3.==“_. B(t-O)-O, B(x,t O) O, ax(x 0) O, DeBax (x l) kng B(x l)] Defining the following dimensionless variables, A : lez ' :._X_- U:A_. K :k3AO V- De 3 _L, uO-As 3-kl A t0 0 E = 9A . b : §_, , z 5;. h . _EA . L ’ - A ’ “2 T k; ’ D 0 e B VeO D n A 9 ep : A : A :th i-‘a‘z— 5- qL V-oe x x A A A A B u:- b:— A0 A0 the mass balances can be made dimensionless as follows 2U -¢ lu+Ka u ]= QJOJ CI): ”a A “ Bu _ _ .gg uf - u - B-gg (g-l,e) - o dB with dimensionless boundary conditions (a= o>=o o;-§! (a=1) = vlfi- u=o;b(;.e=0)=o.—(a=0)=o. <5=1)=hv[B-b(a=m (Sc) Q) J“ Q) {'1 0' Equations (4) are coupled but independent of equations (5); equations (5) are unidirectionally coupled to equations (4). The periodic forcing function Of interest is the square wave 0; tc[O,n/A) uf = ; uf(t+2n/A)=uf(t); bf = O (6) 2; tE[T/A,2W/A) There is no intermediate b in the feed. The steady state anologs of equations (4) and (5) are obtained by setting the time derivatives equal to zero and letting uf equal one. Of the 4 steady state and 4 transient equations, only the steady state counterparts Of the reactant balances (4) were found amenable to solution. The solution to u(§) takes on one Of three simple forms depending on whether u(§=0) >< = l/2K3, or l/2K3 > l. Since the evaluation of u(E=O) requires a knowledge Of the solution form, it is convenient to determine the value of ¢0(=/$%)which corresponds to u(g=0) = l/2K3. By maximum principles and the uniqueness of solution, it can be shown that 935%:9)< 0, so that the value of o determines which form of the solution applies. The result is for O > $0, or if l/(2K3)>l 2 - u(E) = c + (uo-c) dc2{‘/Eg§%igi-E)i} ; ¢>¢o (7a) The value of uo is found for ¢>¢o by solving the following equations simultaneously: 87 _£§2—- /(ul-Uo)(ul-b)(ul-C) = 135v [l-ull (7b) u1 = c + [uo-c] «ZR/Mfg”.C7 } (7c) where b+ C_ '(Uo+3/2K3)i/(%0+3/2K3)2'4U0(“0+3/2K3) (7d) The notation b+,c_ in equation (7d) implies that b corresponds to the plus sign and c to the minus sign. FOE O = $0, the solution takes on a simple form um = uou+3 tanz (g—azn ; we (8) l-cn[/2K3/3g2 Oil l+cn[/2K3/3g2 ¢€J °A, ¢<¢O (96) Here, dc and cn are the Jacobian elliptic functions. Again uo is found for ¢<¢o by the simultaneous solution Of 2 equations, viz., + l-cnfi/2K37392 $] ,A (9b) l+cn[vZK2739 $1 and equation (7b). The quantities needed for evaluation of (9b) are U1 = U0 A2 = (bl'Uo)2 + a12 (9C) _ 2 9 =1/ms bl = mzi 9 af = ' b4C (9d) where b and c are again defined by (7d). When l/(2K3)>l, O0 need not be calculated; otherwise the value of $0 is found by solving _ , 1/2 _ ./—123'< 2 tan {MW—l ”2&5] [ufil/KJ [1.4/sz = OIL—:31 (lOa) 88 for u;. The result is substituted into the following equation and solved for o0: u. = i/2K3{i+3tan2(-¢‘—g)} (lOb) The bulk phase concentration is given by a simple algebraic relationship: “ _ l + BVU(§EJ) u - l+Bv (ll) The steady state counterparts to equations (5) can be solved by variation Of parameters in combination with equations (7-lO). However, the integrals arising from such a treatment contain expressions Of the at at e e 0 0 form J EOE—BE dt or I EEEYEf-dt which cannot be expressed in terms Of a finite series of elementary functions or, to our knowledge, any infinite but convergent series. The periodic equations (4,5,6) are nonlinear partial differential equations and are presumed to have no exact solution. As Of yet no applicable approximate solutions have been found, although such work is in progress. NUMERICAL RESULTS The steady state forms, which give second order ordinary differential equations (ODE's) with boundary values, were solved using the iterative McGuinness method [19]. The periodic forms were treated by the method Of lines [20], in which the spatial derivatives are discretized, yielding a coupled system of ODE's with respect to time. These are solved simultaneously with the ODE governing the bulk phase concentration. All of the ODE's encountered were solved with the IMSL routine DGEAR using Adam's method. The steady state solution for u(g) was checked against numerical simulation and gave arbitrarily close agreement within the accuracy of DGEAR for all three solution regions. Thereafter, in all steady state 89 calculations, including the numerical simulation Of b(g), the analytic expression for u(§) was used. Table 1 summarizes the results of the numerical solution Of the periodic process for square wave forcing. As can be seen for the results with h = 5., the accuracy of the results change dramatically with the number Of spatial increments used. Close convergence Of the results with increasing number Of increments is elusive, because computing time increases very rapidly. In some cases as much as l25 CP seconds on MSU's Cyber 750 model l7O were needed to perform a simulation with 25 spatial increments. Rows l and 2 Of Table l, marked Ref, give the results when all Of > the system constants are equal to one. Row 2 with 20 spatial increments is the more accurate and provides a reference against which the other results can be compared. In each case one constant was varied and all others were held fixed at one. The column on the left shows which constant was varied. Figures l through 5 show the system behavior for the case where all constants are equal to one. There is an abrupt break in the slope of O(B) (Fig. l), whereas 8(a) (Fig. 2) shows no abrupt changes in slope. Figure 3 shows the phase plane for the bulk phase concentrations. Figures 4 and 5 show the concentration profiles within the pellet at five different points in the cycle. The horizontal sections on the right of each profile show the bulk phase concentrations. The diagonal lines connecting the values at the outside of the catalyst pellet to those in the bulk represent a physical gradient which imparts an impression of the concentration difference across the external boundary layer. It is interesting that the bulk phase behavior, Figures 1, 2, 3, is quite similar to that Observed for the square wave forced homogeneous system [l6]. 90 FIGURE 1. Bulk phase reactant concentration vs. time for square wave cycling. ¢=l. v=l. B=l. K3=l. h=l. K2=l. A=l. w=l. 91 0000>0 m omuw 0 000202 comm onus 00.~ 00.0 00. ._ 0000Hm oo-dD O?‘O T OB'O 92 FIGURE 2. Bulk phase intermediate concentration vs time for square wave cycling. ¢=l. v=l. B=l. K3=l. h=l. K2= . A=l. w=l. 93 0000>0 “.0 m~umz0 Dm.N Do. . p .N 0000: 2 nzw.~ oo-rF VO'O 94 FIGURE 3. Bulk phase plane behavior for square wave cycling. ¢=l. v=l. B=l. K3=l. h=l. K2=l. A=l. w=l. 95 2000 0 Dmufi O¢u~ DN.~ DD.~ Dm.o D©.D O¢.D ON.D OO.DU h L . __I . . . 0 U fnu mu 7 To mo 8 8 8 ton hsfil 3“” r0 m ID no 0 .m mm0wH0 96 FIGURE 4. Reactant concentration profile within the pellet at five selected points in the limit cycle. ¢=l. v=l. B=l. K3=l. h=l. K2=l. A=l. h=l. 97 muchwHo mwmu-zHo 00.0 0v." 00.0 00.0 00.0 00.0 0¢.0 00.0 00.00 p p p p P r p p o 0 0 (ill) Till) 0 r 0 IO m. 0 vu.nu 3 0 .... 0 13 .o 0 .v 0000H0 98 FIGURE 5. Intermediate concentration profile within the pellet at five selected pointed in the limit cycle. ¢=l. v=l. B=l. K3=l. h=l. K2=l. A=l. h=l. 99 0020500 00.31200 . 00nd ocua owufi 00mg 00u0 00u0 0vu0 00u0 00 Du m. 8 / I0 .0 8 ID a ununu m 0 r. e 10 r. 8 .0 000000 100 TABLE 1. Percent increase in time average of periodic bulk phase concentrations over steady state concentrations for square wave perturbation. lOl TABLE 1. c 5 v K, h K; 7. : ~51+ 5.0 ’.—b 1. 1. 1. 1. 1. 1. 1. 1. 5 2.43 - 2.03 RBI 1. 1. 1. 1. 1. 1. 1. 1. 20 - .074 - 3.82 Dec .1 1. 1. 1. 1. 1. 1. 1. 5 .340 -19.1 a .1 1. 1. 1. 1. 1. 1. 1. 20 .026 - 4.97 .4 1. 1. 1. 1. 1. 1. 1. 5 1.77 -12 0 Inc 1. 1. 1. 1. 1. 1. 1. 5. 5 3.30 2.11 . 1. 1. 1. 1. 1. 1. 1. 5. 10 1.43 .777 1 1 1 1. 1 1 1 5 20 .637 - .112 Dec 1. 1. 1. 1. 1. . 1. 1. 8 1.09 - 6.62 K; 1. 1. 1. 1. 1. .2 1. 1. 16 .119 - 5.47 Dec 1. 1. 1. .2 1. 1. 1. 1. 10 1.39 .438 k, 1. 1. 1. .2 1. 1. 1. 1. 20 .557 - .421 1. 1. 1. 1. .1 1. 1. 1. 12 .440 - 2.68 1. 1. 1. 1. .2 1. 1. 1. 1o .71 - 2.4 vary- 1. 1. 1. 1. 5. 1. 1. 1. 5 2.49 - 9.68 ‘“9 1. 1. 1. 1. 5 1. 1. 1. 15 .17 - 6.09 h 1. 1. 1. 1. 5. 1. 1. 1. 25 - .o7 - 5.6 1. 1. 1. 1. 1o. 1. 1. 1. 15 .044 - 1.24 lnc v 1. 10. 1. 1. 1. 1. 1. 1. 15 - .377 - 5.86 +NSI is the number of spatial increments used in the method of lines. 102 Discussion Over a narrow range Of physical parameter variation, the preliminary results presented here are less than encouraging for the effort tO increase yield by periodic Operation. Those parameters related to mass transfer resistance, namely o, v, and h, do not seem to combine in any way to enhance time average yield. For instance, decreasing 0, which has the effect of rendering diffusion of A more significant in comparison to reaction of A, appears tO be detrimental to both time average conversion and time average yield. An increase in O can be interpreted in several ways, e.g., the most direct is an increased residence time. This results in a decrease in both conversion and yield. The same effect is produced by a decrease in either K2 or K3. Variation Of h, the ratio of diffusion coefficients for A and B, shows no systematic pattern in its effect on periodic Operation and is not beneficial to yield in any case. It should be noted that changing h should not affect the time average of A. This is closely born out for those cases with l5 or more Spatial increments. Lastly, an increase in v, reflecting lesser external mass transfer resistance relative to diffusion, markedly decreases average yield, but has a slight favorable effect on average conversion. It appears, for those cases considered, that steady state Operation is superior. While no definitive instances of yield enhancement are presented in this work, this effort should not be abandoned. The results presented here reflect only a few combinations of system constants, and the possible variations of eight system constants are enormous. Moreover, the probable error in this calculation via the method of lines is as high as lOO% for the lower values of enhancement. In fact, the results in Table l illustrate that in some cases the enhancement changes sign as the number Of Spatial increments is increased. lO3 It was not possible to achieve strict convergence because Of the excessive computation time involved. Clearly a more efficient numerical technique is necessary to allow a larger range of system constants to be examined. Orthogonal collocation, as used by Bailey [8], may be appropriate. CONCLUSION For the heterogeneously catalysed Van de Vusse reaction occurring in slab catalyst particles in an isothermal CSTR, periodic square wave forcing of the input concentration may not be effective in enhancing yield relative to steady state. Preliminary approximate results show mainly detrimental effects on both yield of intermediate 8 and conversion of reactant A. Additional work,specifically on mathematical analysis Of the integrO-differential equation resulting from the model.is necessary. LIST OF REFERENCES 104 10. ll. l2. l3. l4. REFERENCES--CHAPTER l Dorawala T. G. and Douglas J. M., A.I.Ch.E. J. l97l lZ_974. Ritter A. B. and Douglas J. M., Ind. Engng Chem. Fundls l97O 9 2l. Douglas J. M., Ind. Engng Chem. Proc. Des. Dev. l967 §_43. Douglas J. M. and Rippin D. w. T., Chem. Engng Sci. l966 21_305. Renken A., Chem. Engng Sci. l972 21_l925. Nandrey C. and Renken A., Chem. Engng Sci. l977 §2 448. DeVera A. L. and Varma A., Chem. Engng Sci. l979 §4_l377. DeVera A. L. and Varma A., Chem. Engng J. l979 l]_l63. Bailey J. E. and Horn F. J. M., Chem. Engng Sci. l972 2]_l09. Lee C. K. and Bailey J. E., Chem. Engng Sci. l974 gg_1157. Towler B. F. and Rice R. 6., Chem. Engng Sci. l974 22_1828. Imison 8. w. and Rice R. 6., Chem. Engng Sci. l975 §9_l42l. Bailey J. E., Chem. Engng Commun. l973 1_lll. Farhadpour F. A. and Gibilaro L. 6., Chem. Engng Sci. l975 39 997. 105 10. 11. 12. 13. REFERENCES--CHAPTER 2 Douglas J. M. and Rippin D. N. T., Chem. Engng Sci 1966 21_305. Douglas J. M., Ind Engng Chem. Proc. Des. Dev. 1967 §_43. Dorawala T. G. and Douglas J. M., A.I.Ch.E. J. 1971 11.974- Ritter A. B. and Douglas J. M., Ind. Engng Chem. Fundls 1970 9_21. Renken A., Chem. Engng Sci. 1972 21_1925. Nandrey C. and Renken A., Chem. Engng Sci 1977 §2_448. Farhadpour F. A. and Gibilaro L. G., Chem. Engng Sci 1975 §9_997. Bailey J. E., Chem. Engng Commun. 1973 l_lll. Van de Vusse J. G., Chem. Engng Sci. 1964 19_994. Beek J., A.I.Ch.E. J. 1972 l§_228. DeVera A. L. and Varma A., Chem. Engng Sci. 1979 §4_l377. DeVera A. L. and Varma A., Chem. Engng J. 1979 lZ_l63. Skeirik R. D. and DeVera A. L., Chem. Engng Sci. submitted (1980). 106 10. REFERENCES--CHAPTER 3 Horn, F. J. M., and Bailey, J. E., "An Application of the Theorem Of Relaxed Control to the Problem of Increasing Catalyst Selectivity,“ J. Opt. Thry Appl. 2, 441 (1968). Bailey, J. E., and Horn, F., "Catalyst Selectivity Under Steady-State and Dynamic Operation," Deutschen Bunsen Gesselschaft fur Physikalische Chemie Berichte (Ber. Bunseges. Phys. Chem.) 22 274 (1969). Bailey, J. E., and Horn, F. J. M., "Catalyst Selectivity Under Steady- State and Dynamic Operation: An Investigation of Several Kinetic Mechanisms,” Deutschen Bunsen Gesselschaft fur Physikalische Chemie Berichte (Ber. Bunseges. Phys. Chem.) 23, 611 (1970). Bailey, J. E., and Horn, F. J. M., "Oscillatory Operation Of Jacketed Tubular Reactors," Ind. Engng Chem. Fundls. 2, 299 (1970). Bailey, J. E., and Horn, F. J. M., “Improvement of the Performance Of a Fixed-Bed Catalytic Reactor by Relaxed Steady State Operation," AIChE J. 11, 550 (1971). Bailey, J. E., Horn, F. J. M. and Lin, R. C., "Cyclic Operation of Reaction Systems: Effects Of Heat and Mass Transfer Resistance,“ AIChE J. 11, 818 (1971). Bailey, J. E., and Horn, F. J. M., "Cyclic Operation Of Reaction Systems: The Influence of Diffusion on Catalyst Selectivity," Chem. Engng Sci. 22, 109 (1972). Lee, C. K., and Bailey, J. E., "Diffusion Waves and Selectivity Modifications in Cyclic Operation of a Porous Catalyst," Chem. Engng Sci. 22, 1157 (1974). Oruzheinkov, A. 1., et al., ""Selectivity Of Catalytic Processes Under Forced Cyclic Operation," React. Kinet. Catal. Lett. 19, 341 (1979). Towler, B. F., and Rice, R. G., “A Note on the Response Of a CSTR to a Spherical Catalyst Pellet," Chem. Engng Sci. 22, 1828 (1974). 107 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 108 Imison, B. W. and Rice, R. G., "The Transient Response of a CSTR to a Spherical Catalyst Pellet with Mass Transfer Resistance," Chem. Engng Sci. 29, 1421 (1975). Lehr, C. G., Vurchak, S., and Kabel, R. L., "Response of a Tubular Heterogeneous Catalytic Reactor to a Step Increase in Flow Rate," AIChE J. 13, 627 (1968). Denis, G. H., and Kabel, R. L., "The Effect on Conversion of Flow Rate Variations in a Heterogeneous Catalytic Reactor," AIChE J. 1g, 972 (1970). Al-Taie, A. S., and Kershenbaum, L. 5., "Effect Of Periodic Operation on the Selectivity of Catalytic Reactions," Chem. React. Engng- Houston (Weekman, V. W. and Luss, D., Ed.) ACS Sympos. Ser. (NO. 65) p. 513 (1978). Unni, M. P., Hudgins, R. R., and Silveston, P. L., "Influence Of Cycling on the Rate of Oxidation of $02 over a Vanadia Catalyst,“ Can. J. Chem. Engng, §l, 623 (1973). Skeirik, R. D., and DeVera, A. L., "Periodic Operation of Chemical Reactors 1: Yield Optimization Analysis on the Forced Square Wave Feed Concentration to an Isothermal CSTR with Complex Kinetics," submitted to Chem. Engng Sci. (1980). Abdul-Kareem, H. K., Hudgins, R. R., and Silveston, P. L., "Forced Cycling of the Catalytic Oxidation of CO over a V205 Catalyst-I, Concentration Cycling," Chem. Engng Sci. ;§ 2077 (1980). Abdul-Kareem, H. K., Hudgins, R. R. and Silveston, P. L., "Forced Cycling of the Catalytic Oxidation of CO over a V205 Catalyst-II, Temperature Cycling," Chem. Engng. Sci. 22, 2085 (1980). McGinnis, P. H., Jr., "Numerical Solution of Boundary Value Nonlinear ODE's,“ Chem. Engng Prog. Symp. Ser. §l_(NO. 55) 2 (1965). Ames, W. F., "Numerical Methods for Partial Differential Equations," p. 302, 2nd Ed. Academic Press, New York (1977).