OSCILLATORY OPERATION OF JACKETED TUBULAR REACTORS K. S. C H A N G ’ A N D S. G. BANKOFF
Dejartment of Chemical Engineering, h’orthwestern University, Evanston, Ill.
60201
Optimal control is studied of the cooling jacket water temperature of an empty tubular reactor in which consecutive (Bilous-Amundson) reactions are carried out, subject to sinusoidal perturbations of the feed concentration. In this particular system, second-order effects on the desired outlet concentration were small, and could b e either positive or negative, depending upon the frequency and parametric effects. The improvements seemed to be principally due to the zero-frequency effect, based upon the curvature of the steady-state response surface to inlet concentration, using optimal cooling water temperatures. Dynamic effects decreased the performance. Parametric studies were made of the effects of heat transfer coefficient and heat of reaction.
others, Horn (Horn, 1967, 1968; Horn and Lin, Douglas and Gaitonde, 1967, 1968; Douglas and Rippin, 1966) have demonstrated that significant improvements in the time-average value of a suitable profit function can be obtained by oscillatory, rather than steady-state, operation of lumped-parameter systems. Three effects can be distinguished: (1) a zerofrequency effect, whose sign can be determined from the curvature of the steady-state response surface in the neighborhood of the steady-state optimum; (2) a finite-frequency effect, which depends upon the system dynamics and may exhibit sharp resonances if two or more inputs are varied simultaneously; and ( ( 3 ) an infinite frequency effect, which depends upon a theorem due to Warga (1962). This essentially states that points in the convex hull of the attainable nonconvex set of stat? points can be reached by a relaxed control, corresponding to infinitely fast switching between two or more points in’ the admissible set of controls. The system cannot, of course, respond to infinitely fast switching, and hence remains a t a steady state which is a linear combination of the steady states corresponding to the selected controls, and which may give a larger value of the profit function than any stationary admissible control. The difficulty, of course, is that infinitely fast switching is not possible, and a t finite frequencies the full dy:namic equations have to be taken into account. For distributed-parameter systems the corresponding armament is not yet availa’ble, and a direct variational attack is the approach used in i.his work. A strong maximum principle for distributed-parameter systems was first developed by Butkovskii (1961, 1963) and has since been the subject of a number of other stuldies (Chang, 1967; Chaudhuri, 1965; Denn, 1966; Denn et (21.: 1966; Jackson, 1964, 1966; Katz, 1964; Lure, 1963 ; Sirazetdinov, 1964; Volin and Ostrovskii, 1964; Wang, 1964; Wang and Tung, 1964). After a simple formidation of the necessary conditions for optimality, one of the inputs to the system, which is here an empty tubular reactor in which consecutive reactions A .--, B 4 C take place, is perturbed sinusoidally, and the optimal, time-dependent values of another input, designated as the control, are determined. For the Bilous-Amundson system (1956), small improvements in the time-average yield of the MONG
A. 1967) and Douglas (Douglas, 1967;
desired compound are found compared to the optimum steadystate yield a t the same average value of the perturbed input. Statement of Problem
Consider a tubular reactor in which the first-order consecutive reactions A + ! B +! C take place. Two cases are considered, corresponding to the presence or absence of heat sources due to chemical reaction. At time t = 0, a disturbance in the inlet feed concentration takes place, in response to which the best temperature control, u ( t ) , is to be chosen in order to maximize the time-average yield of component B over a finite time interval 0 _< t _< t,. I t is assumed that the jacket water cooling temperature can be varied as a function of time, and that its velocity is constant and sufficiently large so that spatial variations in the jacket temperature can be ignored. Also, the heat capacity of the jacket is considered to be small compared to that of the reactor, so that the jacket dynamics do not have to be considered. The governing system of partial differential equations for the process is av1
-=
bt
bv2 - _at
klVl
-
bv 1 bX
-v
-klvl
k2v2
- v
bv 2
bX
subject to the initial and boundary conditions: V((O?X)
=
a&),
V&O)
=
IC&),
2
= 1, 2, 3
(2)
where v l ( t , x ) and ~ . ‘ e ( t , x ) are, respectively, the concentrations of A and B; v s ( t , x ) and v ( t ) , the temperature and the velocity of the reactant fluid in the tube; b l = (-AH,/Cp p ) and bz = (-AH2/Cp p ) where ( - A H , ) is the heat of reaction per gram mole of v i , and C, and p are heat capacity and density, respectively; CY = r C, p/2U, a heat exchange parameter, where r is the tube radius and U is the over-all heat transfer ~ (--E,IRT) (i = 1,2), rate constants; coefficient; k 6 = k , exp and a&) and G i ( t ) are, in general, prescribed functions of x and of t, respectively. We wish to maximize the cumulative yield of u 2 in the interval 0 5 t _< t,, corresponding to minimizing the functional
1 Present address, Department of Chemical Engineering, University of Waterloo, Waterloo, Ontario, Canada.
VOL. 7
NO. 4
NOVEMBER 1 9 6 8
633
A computation scheme (Chang, 1967) based upon the method of characteristics (Sirazetdinov, 1964) is particularly advantageous here, since both the state and adjoint equations have the same principal part, V . This is equivalent to saying that, in a plug-flow reactor with no axial diffusion, a marked simplification can be obtained by using a material derivative, d / d s = d / b t Y a/bx, where s is the residence time for a fluid particle. T o solve this problem by a variational method, we first write the state equations (Equations 1) in the form
+
dv - = ds
+ v vz = f ( v , u )
Ut
together with the initial conditions (Equation 2). Alternatively, the periodic boundary conditions (Equations 12 and 13) may be employed. These conditions are the analog of those of Horn and Lin (1 967) for optimum periodic lumpedparameter processes. The interest in this work was in examining short-term, transient effects, as well as the ultimate optimum periodic control. Hence the boundary conditions (Equations 9 and 10) were used, with $ ( t ) a periodic function. I n particular, the adjoint equations become
(4)
where vT = ( V I , V Z , U S ) ; the letter subscripts denote partial differentiations; and f is a vector function of class C2. A set ) ~ (pl, p z , p 3 ) is now introduced, of adjoint functions, p ( ~ , t = together with a Hamiltonian functim
H = pTf
(5)
The state equations are now appended to the objective function (Equation 3) to form a Lagrange function, A = A(v,u,p) :
A =
L'' ' + 1"iL
[H - p T ( u ,
cTv dt
1 x=L
+
with the final-time and boundary conditions: Y
v , ) ] d x dt (6)
where cT is a constant row vector, which in the present case is (0,l:O). For the objective function to be stationary to small perturbations in control, it is necessary that the Lagrange function be stationary to arbitrary weak variations in its arguments : 612 = 0 = L " c T 6v
-
I
x=L
1'' v
pT 6 v
dt
I,
it'iL
+
[ ( H p- utT
dt
x=L
-
x=o
+
6v
dx
(7)
t=O
VPzT =
+l(t) = ulg
+ a sin w t
$&) = vzS
- a sin u t ;
0
(8)
subject to the final-time and boundary conditions :
v
= =
61 =
b1 = together with Equations 2 and 4 and the stationarity requirement:
lL
H,dx = 0
7,x);
vlS = vZs = a
=
a
= =
7
Alterpatively, it may be presumed that the process is strictly periodic in time, so that
+
+
u(t) = u(t
7);
p(t,x) = p(t
+
ups =
1 gram mole per liter
+
v3s
where v l S and v z s are the steady-state inlet concentrations; vSs is the corresponding steady-state optimal feed temperature, determined by a control space gradient method; a is the amplitude of the oscillation; and w is the angular velocity. The parameter values given by Bilous and Amundson (1956) were used: klo = 0.535 X 10" min.-'; kzo = 0.461 X lo1* min.-'; E1 = 18 kcal./gram-mole; E2 = 30 kcal./grammole; R = 2.0 cal./(gram-mole) ( O F ) . Other values which were chosen were :
L
v(t,x) = u ( t
vlS
(18)
$4)=
l L p T I'="
(17)
For oscillatory operation, a steady sinusoidal disturbance is introduced through the inlet conditions:
- Y vzT)
where integration by parts has been employed. This, in turn, requires that the Lagrange multipliers satisfy the adjoint equations:
H, +pzT
p ( t , , x ) = 0 ; vpT(t,L) = (0,1,0)
1 unit of length 0.1 unit of length/min. 0, bz = 0 (case without heat generation) 100, 6 2 = -50 (case with heat generation) 0.80 gram-mole/liter 0.20 gram-mole/liter 0.15 gram-mole/liter 0.5, 3.0, and 10.0 10, 15, 30, and 45 minutes
I t is readily shown that a control gradient algorithm is obtained by relaxing Equation 1 and choosing
7,x)
(1 2) for all t and x . I n particular, the inlet conditions are assumed to be a specified periodic function: v(t,O) = v ( t
+
7,O)
=
$(t) =
$0 + 7)
(1 3)
The integrated terms in Equation 7 will now vanish, provided that the final conditions (Equations 9 and 10) are satisfied, 634
l&EC FUNDAMENTALS
where 6 u ( t ) refers to the control increment, and E > 0 is a step-size parameter. The time-mesh length was chosen to be 0.5 minute, except for T = 10 minutes, where At was taken to be 0.2 minute. A fourth-order Runge-Kutta method was used for integrating along a characteristic line with step size As = 0.01 minute. Values of u ( t ) between time-mesh points were obtained by
0700~Jy 337.50
0710
0
0690 >
331 501
I
I
I
I
,
I
w
00
0 680
I
I
2
7
TIME, min. 0 670 00
20
40
6Q
00
Figure 3.
100
TIME, min.
Effect of frequency of optimal control 01
= 3.0, bl = 0,ba = 0
Figure 1. Optimal control and yield with oscillating feed concentration 7
= 10 minutes,
bl
= 0, bz = 0
- O P T I M A L CONTROL YIELD --- CONSTANT CONTROL YIELD
0 710-
~
30 min
I
0.700 0 .
n
+
v
3
0 330.0
z 0 W
067OL
I
I
I
I
00
I
U'IIWI
I
-r
I
I
I
I
w
J
27r Y
TIME,min.
Figure 4. Yields with optimal time-varying control and with optimal steady-state control CY
0.670L 0.0
60
120
180
240
30.0
TIME, min. Figure 2. Optimal control and yield with feed oscillation 7
= 3 0 minutes, bl = 0, bz = 0
interpolation. As usual, the state equations (Equation 4) were integrated forward, followed by the backward integration of the adjoint Equations 14 to 16. The best steady-state jacket temperature was taken as the initial guess for u ( t ) . Reactor without Heat (Generation
Optimal controls and yields for b l = bz = 0 are shown in Figures 1 and 2. T h e results for 7 = 45 minutes were substantially the same as for = 30 minutes. I t is clear that dynamic effects are absent a t such slow oscillations, so that the improvement can be ascribed to the curvature of the steadystate response curve to inlet concentration variations. As one might expect, poor heat transfer coefficients (large a ) result in large amplitude of the cooling-water temperature variations. There is some effect of a on the yield curves a t the lowest period (10 minutes), but this is negligible at larger periods.
= 3.0,bi = 0,bz = 0
O n e cycle of the optimal control is plotted for different periods in Figure 3. The jacket temperature exhibits almost perfect sinusoidal cycles with, however, frequency-dependent amplitudes, which is a nonlinear effect. If the state equations and the objective function are linear and initial condition is sinusoidal in time, it is readily verified that the solution is everywhere sinusoidal in time, with a space-dependent phase and amplitude. I n this case only the state Equation 1 for u s is linear. The amplitude increases as the period approaches the residence time of 10 minutes, suggesting the appearance of dynamic effects. For r / r l = 1, the control peaks a t t / ~ , 0.5, showing the averaging effect of the control. For 7/r7 > 1, the peak occurs a t larger values, reflecting the expected decrease in phase lag at lower frequencies. Yield curves u z ( t , l ) for r = 10, 15, 30, and 45 minutes are plotted for a = 3.0 in Figure 4. The dashed lines represent the yields obtained from the application to the oscillatory reactor of the best cooling water temperature (constant control) for the steady-state reactor with feed concentrations equal to the time-mean inlet values. These curves are sine waves, since the constant-temperature system is linear. The optimal-control yield curves are located slightly above the corresponding constant-control yield curves, the improvement apparently maximizing when r / r r = 2-3, where r1 is the residence time.
-
VOL. 7
NO. 4
N O V E M B E R1 9 6 8
635
Y 0 . 1 I
c
a -I
0
a
I2
0
V
n
0.690
-I
w>
0,670
0.650 0.0
T I M E , min.
Figure 5. Optimal control and yield with feed oscillation T
60
12.0 18.0 TIME, min.
24,O
30 0
Figure 7. Optimal control and yield with feed oscillation
= 10 minutes, bl = 100, bl = - 5 0
r
= 30 minutes, bl = 100, bt = - 5 0
335.701
1
0 Y
-I
0 U
IF
I-
z
0,690
n
0
u
-I
w
0670 0.650 00
3.0
6.0
9.0
12.0
15.0
TIME ,min. Figure 6. Optimal control and yield with feed oscillation T
= 15 minutes, bl = 100, bz = - 5 0
Reactor with Heat Generation
The optimal controls and yields are obtained as before, with now b l = 100 and b p = -50 (Figures 5 to 7 ) . The effect of heat generation here is to shift the phase of the peak control by about a half cycle. The optimal control is further from a sine wave shape. I n this case none of the state equations is linear. For larger C Y , the amplitude for a given CY is much larger than before, since more control action is now necessary. The optimal yield curves for different CY and the same period are now far apart, the difference being dependent on the values of a. Moreover, the poorer the heat transfer coefficient, the lower is the optimal yield. This phenomenon is clearly seen in Figure 5. However, as the period gradually increases, the optimal yield curves of a = 3.0 and of CY = 10 approach each other (Figure G), and finally become almost one curve, as shown in Figure 7 . The optimal yield curve for CY = 0.5 is still far above the other two. This implies 636
I&EC FUNDAMENTALS
TIME, min. Figure 8.
Effect of frequency of optimal control CY
= 0.5, bi = 100,bz = - 5 0
that, if the jacket heat transfer coefficient is sufficiently large, the adverse effect of heat generation can be greatly reduced. The optimal yield curves are now strongly distorted, especially for small 7 , as can be observed with a = 3.0 and 10 in Figures 5 and 6. For a given period the distortion of the yield curve becomes more noticeable as the heat transfer becomes poorer. Again the phase lag decreases with the frequency, but is higher than in the absence of heat sources. At the same time the optimal yield amplitudes are larger than those for the corresponding optimal yield curves without heat generation, reflecting the adverse influence of heat generation from the point ofview of system stability. The optimal jacket temperatures are plotted for CY = 0.5 and 3.0 in Figures 8 and 9, respectively. For a = 0.5, because of the good heat transfer, the curves are similar to those without heat generation. There is an averaging effect when T = 10 minutes, the residence time, causing the temperature
327'00K
/
Figure 9. control
Effect of frequency of optimal
a = 3.0,bi = 100,b2 = - 5 0
TIME, min.
Figure 10. (Y
t 0.650 0.0
Optimal control yield
= 3.0, bi = 100,bn = - 5 0
I I
I
I
7 -
I
I
I
W
I
27 W
TIME, min.
OPTIMAL CONTROL YIELD CONSTANT CONTROL YIELD
IO rnin.
0.665 Figure 1 1. Yields with optimal time-varying control and with optimal steady-state control OL
0.0
7r
W
= 3.0, bi = 100,bz = - 5 0
3 W
TIME, min. VOL. 7
NO. 4
NOVEMBER 1 9 6 8
637
0.710 0 6950 0
-
*
0.5
0
.-c
0.705
5
0.700
1 0 0
a=3Q
I
n
->
t
a > a
0,5950
05750
0.690
I
B 0 z 8
- CONSTANT CONTROL YIELD --- OPTIMAL CONTROL YIELD
1 -0
0.5550.
0.695
g
0.6150
W
mi
z
0.6350
W
$
A
I
10.0 I
I
I
I
I
0.685
k
E
0.680
1
0 I
0.65
L
0.75 0.80 0.85 INLET CONCENTRATION OF A
0.70
0.90
Figure 13. Inlet concentration vs. optimal yield of B for optimal space-dependent steady-state control Midpoint of upper chord represents cumulative yield for infinitely fast switching between peak concentrations; lower chord infinitely slow switching between same two values
to be nearly constant. With poorer heat transfer, the optimal controls tend to be more anticipatory, resulting in a phase shift of nearly 180'. The phase again increases and the amplitude decreases as T increases. The pattern of the optimal controls for CY = 10 is almost the same as that in Figure 9 for CY = 3.0, but with much larger amplitudes and lower mean values. Figure 10 shows the optimal yield curves with a = 3.0. Because of the smaller averaging effect, the amplitude increases with increasing period. The nonlinear effects are clearly visible a t the higher frequencies. For CY = 3.0, the comparison of the yields from the optimal controls with those from the constant control is shown in Figure 11 for 10- and 30-minute periods. The difference in the two yields is significant, showing a synergistic effect for two input time functions as compared to only one (inlet concentration). Similar results are obtained with other values of CY and T , as summarized in Figure 12. The difference is particularly marked with poor heat transfer.
0'6950
4 t
STEADY STATE YIELD
--- OPTIMAL OSCILLATORY YIELD
Comparison of Yield of Oscillatory and Steady-State Feed
Figure 13 shows the steady-state optimal yield response surface with the best space-dependent control, computed by a gradient procedure. This then represents an upper limit to the zero-frequency response surface, since the averaging effects due to the use of a control dependent upon time only cannot improve the yield. For a zero-frequency oscillation with amplitude of 0.125 around a mean inlet concentration of 0.775 gram mole per liter of A, the increase is seen to be roughly 0.1% based upon the crude estimate that the concentration lies, on the average, half the time at the upper midpoint (0.8375 gram mole per liter) and half at the lower midpoint (0.7125). If, on the other hand, the switching between the extreme values could be performed with infinite frequency, the result would be a steady state on the upper chord in Figure 13, corresponding to an improvement of about 0.3%. 638
l&EC FUNDAMENTALS
0,6700
10.0
20.0
x). 0
40.0
PERIOD, min. Figure 14. Optimal oscillating yield compared to optimal steady-state yield bl = 100, bz = -50.
In both cases control is space-independent
A comparison of the time-average optimal yields with the optimal steady-state yields for the same average inlet concentration is shown in Figure 14. I t is seen that the increment may be negative for T / T ~= 1, but for T / T ~2 3 it reaches 0.1 to 0.30j0, depending upon the heat transfer coefficient. The dynamic effects, when the oscillation period approximates
the residence time, here decrease the time-average performance, in view of the unwanted averaging effect of the spatially uniform control. For 7/rr > 1, the decrease is monotonic, indicating the probable absence of resonance effects for this particular system, excited by this choice of input perturbations. Nomenclature = = = =
a
bi, b i c
C2
c,
= = = AH2 =
El, E2
H AH1, J
k l , kp klo, kg0
L
pl,
R
p2, pa
r S
T t t,
U U
vi, X
u g , u3
= = = = = = = = = = = = = = =
amplitude of input concentration wave
( - A H i ) / c p p , (-AHz)/C,p constant vector class of twice continuously differentiable functions specific heat activation energies Hamiltonian function heats of reaction performance index rate constants pre-exponential factors length of reactor adjoint variables gas constant radius of tube parameter along characteristic curve (particle residence time) temperature, u3 time final time over-all heat transfer coefficient control variable (jacket water temperature) state variables (concentration of component A, concentration of B, reactant temperature) distance firom inlet of reactor
GREEKLETTERS = heat exchange parameter, rCpp/2U = variation = iteration step size = velocity of reacting fluid = density = period, 2 ~ / w = residence time = frequency
A @I, @2,
= @3
@I, @p, @ 3
Lagrange function
= initial conditions =
boundary conditions
SUBSCRIPT S
= steady-state
literature Cited
Bilous, O., Amundson, N. R., Chem. Eng. Sci. 5 , 81-92, 115-26 (1956). Butkovskii, A. G., Autom. Remote Control 22, 13, 1156, 1429 (1961); 24, 1106 (1963). Chang, K. S., Ph. D. thesis, Northwestern University, Evanston, Ill., 1967. Chaudhuri, A. K., Intern. J . Control 2, 365 (1965). Denn, M. M., Intern. J . Control 4, 167 (1966). Denn, M. M., Gray, R. D., Jr., Ferron, J. R., IND.ENG.CHEM. FUNDAMENTALS 5 , 59 (1966). Douglas, J. M., Ind. Eng. Chem. Process Design Develop. 6 , 43-8 (1967 ). Douglas; J. M., Gaitonde, N. Y., A. I. Ch. E. meeting, St. Louis, Mo.. 1968. Douglas, J. M., Gaitonde, N. Y . , IND.ENG.CHEM.FUNDAMENTALS 6,265-76 (1967). Douglas, J. M., Rippin, D. W., Chem. Eng. Sci. 21, 305-15 (1966). Horn. F. J. M.. A. I. Ch. E. meetinc. St. Louis. Mo.. 1968. Horn; F. J . G., Ind. Eng. Chem. T r a m s De& Dwelip. 6 , 30-5 (19671. Horn, F: 3. M., Lin, R. C., Ind. Eng. Chem. Process Design Deoelop. 6, 21-30 (1967). Jackson, R., I. Ch. E.-A. I. Ch. E. Joint Meeting, London, 1964. Jackson, R., Intern. J . Control 4, 127 (1966). Katz. S..J . Elect. Control 16. 189 (1964). Lure; K: A,, J . Appl. Math.'Mech.'27, 1284 (1963). Sirazetdinov, T. K., Autom. Remote Control 25, 431 (1964). Volin, Y . M., Ostrovskii, G. M., Autom. Remote Control 25, 1414-20 (1964). Wang, P. K. C., I E E E Trans. Autom. Control AC-9, 13 (1964). Wang, P. K. C., Tung, F., J . Basic Eng. Trans. A S M E D86, 67-79 (1964). Warga, J., J . Math. Anal. A j p l . 4, 111, 129 (1962). RECEIVED for review October 13, 1967 ACCEPTEDApril 19, 1968 Work supported by the National Science Foundation (N.S.F. Grant GK-1126), together with the I.P.A.C. System Laboratory, Northwestern University (Office of Naval Research Contract Number NO0014-67-A-0356-0003 Modification AB).
VOL. 7
NO. 4
NOVEMBER 1 9 6 8
639
