Ind.
Eng. Chem. Res. 1987,26, 325-333
Tulsa, OK, April 15-18, 1984; paper S P E 12644. Heller, J. P.; Lien, C. L.; Kuntamukkula, M. S. Presented a t the Annual Technical Conference and Exhibition of the Society of Petroleum Engineers, New Orleans, LA, Sept 26-29,1982; paper SPE 11233. Hirasaki, G. J.; Lawson, J. B. Presented at the Annual Technical Conference and Exhibition of the Society of Petroleum Engineers, San Francisco, CA, Oct 5-8, 1983; paper SPE 12129. Holm, L. W. Trans. AIME 1959, 216, 225. Kraynik, A. M. Presented at the Annual Meeting of the Society of Rheology, Evanston, IL, Oct 24-28, 1982; paper B-10. Krug, J. A. Oil Gas J. 1975, 73(0ct 61, 74. Lincicome, J. D. World Oil 1984, 198, 29. Marsden, S.S.; Khan, S. A. SOC. Pet. Eng. J. 1966, 6, 17. Mast, R. F. Presented at the Annual Meeting of the Society of Petroleum Engineers, San Antonio, TX, Oct 8-11, 1972; paper SPE 3997. Merrill, E. W. In Modern Chemical Engineering; Acrivos, A., Ed.; Reinhold: New York, 1963; Vol. I, pp 141-196. Minssieux, L. J . Pet. Technol. 1974, 26, 100. Mitchell, B. J. Oil Gas J . 1971, 69(Sept 6), 96. Mooney, M. J . Rheol. (Easton, Pa.) 1931, 2, 210. Okpobiri, G. A.; Ikoku, C. U. J. Energy Resour. Technol. 1983,105, 542. Patton, J. T.; Belkin, H.; Holbrook, S.; Kuntamukkula, M. report 78 MC 03259, Oct 1979; US Department of Energy. Patton, J. T.; Holbrook, S. T.; Hsu, W. SOC.Pet. Eng. J . 1983, 23, 456.
325
Patton, J. T.; Patton, J. T., Jr.; Kuntamukkula, Tf, S.; Holbrook, S. T. Polym. Prep (Am. Chem. SOC.,Diu. Polym. Chem.) 1981,22, 2, 46. Princen, H. M. J. Colloid Interface Sci. 1983, 91, 160. Prud’homme, R. K. Presented at the Annual Meeting of the Society of Rheology, Louisville, KY, Oct 11-15, 1981; paper E-7. Rand, P. B. US Patent 4 442 018, April 10, 1984. Raza, S.H. SOC.Pet. Eng. J . 1970, 10, 328. Raza, S. H.; Marsden, S. S. SOC.Pet. Eng. J . 1967, 7, 359. Rosen, M. J. Surfactants and Interfacial Phenomena; Wiley: New York, 1978. Sanghani, V.; Ikoku, C. U. J. Energy Resour. Technol. 1983,105, 362. Sebba, F. Chem. Ind. 1984, IO(May 21), 367. Sibree, J. 0. Trans. Faraday SOC.1934, 30, 325. Wellington, S., Shell Development Co., Houston, TX, personal communication, 1982. Wendorff, C. L.; Ainley, R. B. Presented at the Annual Fall Meeting of the Society of Petroleum Engineers, San Antonio, TX, Oct 5-7, 1981; paper SPE 10257. Wenzel, H. G.; Brungraber, R. J.; Stelson, T. E. J . Mater. 1970, 5, 2, 396. Wenzel, H. G., Jr.; Stelson, T. E.; Brungraber, R. J. J . Eng. Mech. Am. SOC.Ciu. Eng. 1967, 93(EM6), 153.
Received for review February 22, 1985 Revised manuscript received March 31, 1986 Accepted August 27, 1986
Nonlinear and Adaptive Parameter Estimation Methods for Tubular Reactors Willi Rutzler Chemical Engineering Program, T h e University of Alabama in Huntsville, Huntsville, Alabama 35899
The problem of nonlinear pseudo-steady-state parameter estimation for tubular reactors is addressed. I t is shown t h a t constant-gain, nonlinear, Kalman filters are feasible tools for estimating catalyst activity and the heat-transfer coefficient. Inclusion of a concentration measurement in addition t o temperature measurements brought a large improvement in the correlation of t,he estimates of two parameters. An adaptation algorithm is presented t o estimate one parameter for cases where process and measurement noises are unknown. It finds the optimal filter gain from data of suboptimal filter performance. 1. Introduction The control of chemical processes is often hindered by the unavailability of good process models. If accurate models are available, the methods of classical and modern control theory can be applied to get smooth operation. In the chemical process industries, however, the modern methods of multivariable control have not been widely used, partly due to the fact that good dynamic models for a chemical process are difficult to obtain. One approach taken in the chemical process industries has been to develop controllers based on models derived from experimental plant testing, such as step or impulse responses (Richalet et al., 1978; Cutler and Ramaker, 1979). Reviews of these methods are available (Martin, 1981; Garcia and Morari, 1982). The advantage of using stepresponse and impulse-response methods is that they can be obtained experimentally and that no knowledge of the physicochemical phenomena is required. The major drawback in using such models is that no advantage is taken of available information about the process. Recent control work in research laboratories (Fisher and Seborg, 1976; Silva et al., 1979; Bonvin et al., 1980, 1983; Wong et al., 1983) indicates, however, that multivariable control could be a powerful tool provided a good model of the process is available. This implies that 0888-5885/87/2626-0325$01.50/0
a good working knowledge of the physicochemical phenomena must be available and that parameters are known with reasonable accuracy. Tubular reactors are particularly difficult processes to control, since they are distributed parameter systems having both spatial and temporal coordinates. Since they constitute an important component in many chemical plants, their control has received considerable attention in the chemical engineering literature. Jutan et al. (1977), Silva et al. (1979), Bonvin et al. (1980, 1983),and Wong et al. (1983) found that the application of modern control theory, i.e., multivariable control theory in the time domain, critically depends on the availability of accurate models for the process. I t is a rare case, however, when a model can be built from theory alone. Often only the structure of the process is known, while parameters are poorly known and change with time. In a packed-bed tubular reactor, which is used as an example in this work, the structure of the model is well understood, but parameters such as catalyst activity or heat-transfer coefficient are unknown and change with time. In previous work on tubular reactor control, parameters either were estimated before a control experiment was done (Jutan et al., 1977; Silva et al., 1979; Wong et al., 0 1987 American Chemical Society
326 Ind. Eng. Chem. Res., Vol. 26, No. 2, 1987 1983) or were assumed constant (Georgakis et al., 1977).
The availability of powerful estimators could provide the accurate models needed for the use of multivariable controllers and could therefore be the missing link for the practical application of the vast literature on multivariable control to industrial reactors. Nonlinear pseudo-steady-state estimators for parameter estimation on tubular reactors have been presented by Gavalas and Seinfeld (1969) and by Clough and Ramirez (1979). By use of a pseudo-steady-state approximation, the reactor is considered to be a t steady state for the purpose of parameter estimation. This is justified if the parameters change slowly, compared to the dynamics of the reactor. Goldman and Sargent (1971) estimate slowly changing parameters in a distillation column and in a tubular reactor with a pseudo-steady-state estimator. They model the parameters as constants and find that the filter does not converge. They overcome this problem by employing a “forgetting factor“ to discard old data. Gavalas and Seinfeld (1969) also observe poor convergence when they use a parameter modeled as a constant. In the present work it is shown that these convergence problems can be dealt with by modeling the parameters as random walk processes and by using a constant-gain, nonlinear filter to estimate the slowly changing parameters. Working with a constant gain becomes possible because the parameters are modeled as random-walk processes instead of as constants. Thus, it is shown that such a seemingly minor change in the problem formulation is of considerable significance. One drawback in the use of recursive filters is that their optimality depends on design parameters that are often unavailable or not known very accurately; i.e., the filter design equations require knowledge of the process-noise and measurement-noise covariance matrices. If these design parameters are inaccurate, the filter will be suboptimal. It is possible, however, to use the data from the performance of a suboptimal filter in an adaptation scheme to obtain better design parameters (Mehra, 1970). Carew and Belanger (1973) present an adaptation algorithm to directly compute the optimal filter gain from suboptimal filter data. In the present work this algorithm is applied to the adaptive estimation of a reaction rate parameter in a packed-bed reactor. 2. Statement of Problem The problem addressed is the estimation of parameters in a packed-bed tubular reactor for the production of phthalic anhydride from o-xylene. A static model of the reactor is developed, and two key parameters, the catalyst activity and the overall heat-transfer coefficient, are assumed to be unknown and changing with time. They need to be estimated in order to obtain an accurate reactor model. The objective of this work is to show that pseudostatic, nonlinear, and adaptive estimation schemes are feasible tools to be employed for reactor-modeling purposes. 2.1. Modeling of the Reactor. Several models for the phthalic anhydride reactor have been presented in the literature (Froment, 1967; Stewart and Sorenson, 1972). For the purpose of the estimation studies in this work, a pseudohomogeneous and radially lumped reactor model is developed. Carberry and White (1969) note that even for a highly exothermic reactor the catalyst particles remain nearly isothermal. Intraparticle concentration gradients are small if the particles are sufficiently small, while for larger catalyst particles internal concentration gradients must be accounted for. For the reactor considered here, the particles are assumed to be sufficiently small to allow
B h
n
M”
Figure 1. Kinetic scheme for the oxidation of o-xylene.
the pseudohomogeneous approximation to be valid. The reaction scheme used has been reported by Froment (1967) and consists of one parallel and two consecutive reactions (Figure 1). Oxygen is present in large excess, and therefore, the reaction kinetics follow a pseudo-firstorder law. The reaction rates for the three reactions are rl
*
= k 1re-E1/RpT*CA*
(1)
The pseudohomogeneous static reactor model consists of three ordinary differential equations representing two material balances and an energy balance. dCA* V -= p ( 1 - c)(-rl* - r3*) (4) dz dCB* D - = p(1 - c)(rl* - rz*) (5) dz dT* PFpp = 2u ~ (-1c ) ( - M l r l * - AH2?-**- M3r3*)- f ( T * - TJ (6)
With the initial conditions at z = 0 CA* = CA*in CB*
=0
T* = T*in We define the following dimensionless variables and parameters
where
pc, = (1 - 4PSCPS + CPfcpf Equations 4-6 can be written as
Ind. Eng. Chem. Res., Vol. 26, No. 2, 1987 327
.
with initial conditions a t [ = 0 of =
CA
CAin
CB = 0 T = Ti,
Ziji = Ziji-1 - Zi/i-lHT(HZi/i-lHT+ R)-'HZi/i-l zi+l/i= Zi/i
2.2. Formulation of the Estimation Problem. The parameters p and v are considered to change slowly with time, and therefore, updated values for their estimates must be computed at intervals. The parameter estimation problem will be formulated in such a way that a nonlinear recursive estimator can be used. The parameters are written as a vector
[:I
k =
The parameter vector, k, is modeled as a random-walk process, i.e., ki+l
=
ki
+ wi
(10)
where wi is a random, white sequence with zero mean and covariance Q, i.e.,
E [ w ~ w= ~Q~ ] E [ w J = 0 E[wiwnT]= 0 for i z n The subscript i denotes discrete time steps, and E[ ] denotes the expectation operator. Possible measurements include the temperature a t different locations along the length of the reactor. Therefore, the equation relating the parameters and the measurements is a nonlinear equation which is defined implicitly by eq 7-9, i.e., E[uiuiT] = R yi = @ki) + ui E[ui]= 0 E[uiunT]= 0 for if n
(11)
yi is a vector of measurements, and h is a nonlinear function relating the parameters and the measurements through eq 7-9. Equations 10 and 11constitute the model equations for the estimation problem. The measurement-noise covariance matrix, R, is taken as a diagonal matrix with the variances of the measurement errors as diagonal elements. Q is taken as a diagonal matrix. The diagonal elements are referred to as process-noise variances and reflect the change in magnitude that the parameters are likely to experience within one time interval. The performance of the estimator will depend on how well Q is chosen. 3. Pseudo-Steady-State Filter The estimation problem stated in eq 10 and 11 can be solved through a constant-gain, nonlinear estimator. The gain is computed by using a linearized version of eq 11and linear Kalman filter theory (Swerling, 1959; Kalman, 1960). Linearizing eq 1'1 yields (12) yi = Hki + ~i H is the Jacobian matrix obtained by linearizing eq 11and the elements of H represent the sensitivities of the measurements to changes in the parameters, i.e.,
- -
Hj, =
(134
k.l + l / l . = k l./ l.
:1 I
The equations for the nonlinear filter are obtained by using linear filter theory for the computation of the filter gain (eq 13c and 13d) and using the nonlinear model in the update equations (eq 13a and 13b)
+Q
(13~) (13d)
where Ki = Zi/i-lHTIHZi/i-lHT+ R1-l. kiji denotes _estimates a t time i, using the measurements to time i. kili-l denotes an estimate a t time i, using the measurements to time i - 1. and Zili-l are the corresponding error covariance matrices. For fixed H, R, and Q, the algorithm consisting of eq 13c and 13d is independent of the measurements. I t converges to a steady-state Zi l and ZiliVl, provided the system of eq 10 and 12 is observahle. A proof of this is available in Anderson and Moore (1979). The constant filter gain becomes
K
+
= Zi/i-lHTIHEi,i-lHT R]-'
(14)
The possibility of using a constant gain is a direct result of modeling the parameters as a random-walk process. If the parameters were modeled as constants, the constant gain would be zero and no measurements would be processed. If the algorithm (eq 13c and 13d) is used and the parameters are modeled as constants, Ki 0 as i and filter saturation occurs; i.e., new measurements are ignored.
- -
4. Problem of Poor Observability Consider a system
= Anxnki (154 Yi = Hmxnki (15b) The following theorem is one of the cornerstones of the linear system theory and is available, e.g., in Kwakernaak and Sivan (1972). Theorem: System 15 is observable if, and only if, the matrix ki+l
I
H A J
spans the n-dimensional space. For a system of the particular form used in this work, A = I and therefore
=
E1
The condition for observability can be expressed in terms of H, because for 0 to span the n-dimensional space it is necessary and sufficient that H spans the n-dimensional space (Le., rank H = n). Therefore, we can distinguish three types of unobservability: (a) One or more columns of H are zero, and therefore, rank H < n; zero columns of H imply that parameters corresponding to zero columns do not affect any of the measurements and therefore cannot be estimated. The only way to remedy the situation is to add one or more measurements until all columns of H contain at least one nonzero element. (b) Two columns of H are linearly dependent; this type of unobservability manifests itself in the impossibility of estimating the two parameters simultaneously. (c) The number of measurements is less than the number of parameters to be estimated. In this case the system
328 Ind. Eng. Chem. Res., Vol. 26, No. 2, 1987
S1
'1
I
j
l
l
I
P I
0. m
10
1
I
m
I
1c 2c
24 00
20 00
I
1
10 00
so TO
0.00
TIME
I
1 1O.W
20.M
I
I
30.00
4O.W
1
I 50.00
60.02
TIME
F i g u r e 2. Estimation of catalyst activity, six temperature measurements, q = 0.0001. (1)Catalyst activity p,; (2) estimate of catalyst activity, p,,,.
is always unobservable since rank H < n. The only remedy is to add additional measurements. In the estimation problem for the phthalic anhydride reactor considered above, the columns of H are almost linearly dependent if only temperature measurements are used. It is therefore expected that the estimates of the catalyst activity parameter, p, and the heat-transfer coefficient, u, will be highly correlated.
Figure 3. Estimation of catalyst activity, six temperature measurements, q = 0.0004. l and 2 are the same as in Figure 2.
8
1
5. Numerical Results for the Estimation of Reactor Parameters The sensitivity matrix, H, of eq 12 is computed by solving eq 7-9 numerically for the three cases (1) p = 1.0 u = 1.0 (2)
p =
(3)
p
1.01
= 1.0
u
= 1.0
u =
1.01
The partial derivatives, ( d T / d p ) , and (dT/du),, are computed numerically a t six temperature measurement locations ([ = 0.025,0.050,0.075,0.125,0.175,and 0.250). The partial derivatives of the product concentration, (dCB/dp), and (dCB/dv),, are computed numerically at the reactor exit. Matrix H contains the partial derivatives with respect to the parameters as its elements. For temperature measurements only, the first column contains ( d T / d p ) , and the second column contains ( d T / d v ) , a t six measurement locations.
H
=
[ I 0.025 0.051 0.078 0.109
-0.010
-0.030 -0.055 -0.091 0.074 -0.080 0.010 -0.036
In order to test the performance of the constant-gain filters, a sinusoidal change of the catalyst activity was used as a test function. Gaussian measurement-noise is simulated with a random-number generator. 5.1. Estimation of the Catalyst Activity. For the estimation of one parameter, the process-noise matrix, Q , becomes a scalar, q. Figures 2 and 3 show graphs of the catalyst activity and its estimate vs. time for different design values, q = 0.0001 and 0.0004. The measurement covariance matrix for all calculations is R = diag(0.000005), which corresponds to uncorrelated measurement errors
om
I
IO M
I 20 M
I
30 M
I
10 00
I 50 00
I 6000
TIME
Figure 4. Estimation of catalyst activity and heat-transfer coefficient, six temperature measurements, Q = 0.000011. (1) Catalyst ~ ;heat-transfer activity p,; (2) estimate of catalyst activity I ~ ,(3) coefficient u t ; (4) estimate of heat-transfer coefficient C,,,.
with standard deviations of 1.5 K. One unit of time corresponds to 20 min, a time long enough for the reactor to reach a pseudo steady state. For q = 0.0001, the estimates are smooth but the estimator tracks poorly. For q = 0.0004, excellent tracking is obtained without excessive noise. Figure 4 shows graphs for the simultaneous estimation of the catalyst activity, p , and the overall heat-transfer coefficient, u. p is assumed to follow a sinusoidal change, while u remains constant (this fact, of course, is not known t o the estimator). As expected the estimates deteriorate if two parameters are estimated simultaneously and a large positive correlation of the estimation errors is obtained; i.e., if the estimate of the catalyst activity is too low, very likely the estimate of the heat-transfer coefficient will also be too low. Since the two columns of the sensitivity matrix, H, are almost linearly dependent (Le., an increase in catalyst activity has a similar effect on the temperature profile as a decrease in the heat-transfer Coefficient), it is extremely difficult for the estimator to separate the effects of the two parameters on the temperature measurements due to poor observability.
Ind. Eng. Chem. Res., Vol. 26, No. 2, 1987 329
']
Table I. Estimation of Catalyst Activity 4 r Zili-1 0.000 005 0.195 X 0.0001 0.538 X 0.0004 0.000 005 0.000005 0.105 X 0.0009 r 2wr.S = (HW'H)-' weighted least-squares 0.000 005 0.186 x 10-3
8
Tables I and I1 list the values of the error covariance matrix ZL/L-l,for the different cases m well as for a weighted least-squares estimator. Table I1 also lists the correlation coefficient and the values of det (Zi/L.-l), a measure of the area of the confidence region of the parameter estimates and, therefore, a natural criterion for the quality of the estimates (Box and Lucas, 1959; Draper and Hunter, 1966). For comparison, the corresponding values are also calculated for a non-Bayesian weighted least-squares estimator for which ZwLs = (HTR-'H)-'
1
1 '
0
10.00
I 20.00
I 10.00
I
1 50.m
i0.W
1 60.00
TIME
Figure 5. Estimation of catalyst activity and heat-transfer coefficient, six temperature measurements and one concentration measurement, Q = 0.000011. 1-4 are the same as in Figure 4.
For lower values of q, the filter performs better, while for larger values of q the least-squares estimator seems to be better. Higher values of q, however, imply quickly changing parameters, and the weighted least-squares estimator could become infeasible because of the limited time available to perform the iterative calculations on-line. Since catalytic reactors have slow dynamics, a weighted least-squares estimator is feasible, but for reactors with faster dynamics, such as empty tubular reactors, the filtering approach becomes more attractive. The main advantage of the filter is that information is carried from one time step to the next and, therefore, the temporal redundancy of information is exploited. The correlation coefficient 2
w=-
ULIU UpQU
is drastically reduced for the filter compared to the weighted least-squares estimators ( p = 0.66 and 0.56 vs. p = 0.96). Adding one concentration measurement improves the performance of the estimators considerably (Figure 5), and it is demonstrated that having one temperature and one concentration measurement is superior to having six temperature measurements (Figure 6).
T
lo
I
10.00
,
I
30.M
20.02
1
9.00
4O.W
Figure 6. Estimation of catalyst activity and heat-transfer coefficient, one temperature and one concentration measurement, Q = 0.00011. (1)Catalyst activity p,; (2) estimate of catalyst activity b+; (3) heat-transfer coefficient u,; (4) estimate of heat-transfer coefficient
surement-noise are not known perfectly. In such a situation the estimator performance will be suboptimal. Heffes (1966) studied the effect of inaccurate knowledge of the process and noise covariance matrices, Q and R, on the filter performance of discrete systems, and Nishimura (1967) studied the problem for continuous systems. Several authors have proposed schemes to obtain estimates of the noise statistics from suboptimal filter data. Mehra (1970) solved the problem of estimating Q and R from correlation data of the innovation sequence. Carew and Belanger (1973) developed an alternate algorithm to com-
6. Adaptive Estimation The optimality of the pseudo-steady-state estimates depends on the accuracy of the information used in designing the estimators. The Kalman filter provides the optimal estimate if the model is linear and accurate and if all noise statistics (i.e., process-noise and measurement-noise) are known exactly. In a realistic situation, the information available is not perfect. The process model is not perfect and often nonlinear, and the statistics of the process-noise and mea-
Table 11. Estimation of Catalyst Activity and Heat-Transfer Coefficient 9 0.0001
r 0.000 005
0.0004
0.000 005
weighted least squares
P
~,/,-l
0.442 X 0.321 X 0.1021 x 10-2 0.6185 X r 0.000 005
,
60.00
TIME
0.321 X 0.538 X 0.6184 X 0.1206 X
lo-'
ZwLs = (HTR-'H)-' 0.254 x 0.282 x 10-2 0.282 X lo-' 0.338 X
1$/1-11
0.658
1.34 x 10-7
0.557
8.48 X P
0.962
IZWLSl
6.32 x
330 Ind. Eng. Chem. Res., Vol. 26, No. 2, 1987
pute the optimal filter gain from correlation data of the innovation sequence. In the current work, the adaptation algorithm of Carew and Belanger (1973) is employed to determine the optimal filter gain from data of suboptimal filters. 6.1. Adaptive Filter. In the following we will consider a one-dimensional system ki+l = ki + wi E[wi]= 0 E[wi2]= q (16a) yi = hki + ui
E[ui2]= r
E[ui]= 0
(16b)
where k is a parameter to be estimated, y is a measurement, and wi and ui are white sequences with zero mean and covariances, q and r, respectively. Following Carew and Belanger (1973), an adaptation algorithm is derived. The optimal estimate is given by the Kalman filter k i + l j i = hili-1 + K(yi - hki,,-,) (17) where K is the optimal steady-state gain which is assumed to be unknown. A suboptimal filter with gain K* is k*j+l/i = k*iji-1 + K*(yi - hk*i/i-l) (18) where * denotes suboptimal gain and estimate. Rewriting eq 18 and adding and subtracting K*hkiji-l yields k*.l + l ./ l = k*i/i-l K*si + K*h(kiIi-,- k*iii-l) (19)
+
where si = y i - hki is the innovation sequence for the optimal filter. Su6iracting eq 19 from eq 17 yields k.l + l / l . - k*i+lji = (1- K*h)(kijj-l - k*i,i-l) + (K - K * ) s ~ (20) Squaring and applying the expected value operator, one obtains E[(Ai+l/i- ~z*i+l/i)~l = (1 - K*h)2E[(ki/i-1 - k*i/i-J2] + ( K - K * ) 2 E [ ~ i+ 2] 2(1 - K*h)(K - K*)E[si(ki/i_, - k*i/i-l)] (21) The innovation sequence of ,an optimal filter is a white sequence, uncorrelated with (ki,i-l - k*i/i-l),and therefore, the last term in eq 21 is zero. Defining (224 P*i+l/i = E[(hi+l,i - k * i + l j i ~ ~ l and
w = E[s?]
(22b)
one obtains p*j+l/i=
+
(1- K*h)2p*i/i-1 ( K - K * ) W
(22c)
Since the filter is a t steady state, ~ * i + ~ =/ ip*i,i-1 = jj*+1 p*i/i-1 (1- K*h)2p*iji-1 + ( K - K*)W (22d) The innovation sequence for the suboptimal filter is defined as si* = y . - hk*.1/1-1 (22e) and the correlations of the innovation process are defined as c, = E { U i * V i , * ) (23) For the scalar system, Co and C1 will be needed for the algorithm to calculate the optimal gain. A derivation of Co and C, is given in Appendix B. We have Co = W + h2p*i/i-l (24) C1 = h2(1 - K*h)p*i/i-l
+ h(K
-
K*)W
Combining eq 24 and 25, we obtain C, + hK*Co = h[p*i,i-lh+ KW]
(25) (26)
I
1o.m
f
20.m
I
so.m TIME
I
4n.m
I
5D.m
I
6O.m
Figure 7. Estimation of catalyst activity, one temperature measurement a t [ = 0.125, filter gain K* = l. l and 2 are the same as in Figure 2.
Therefore, the following adaptation algorithm is obtained from eq 24, W = Co - h2jj*i/i-1
(274
+ hK*Co)/h -jj*iji-lh]/W
(27b)
from eq 26, K = [(C, and from eq 22d, p*i/i-1 = (1 - K*h)'P*i,i-l
+ ( K - K*)'W
(27~)
A good choice to start the algorithm is p*+-, = 0 since p*+, = 0 if K = K*. Convergence is guaranteed if the model is perfect, linear, and observable and if good estimates of Co and C1 are available (Carew and Belanger, 1973). 6.2. Adaptive Estimation of the Catalyst Activity in the Phthalic Anhydride Reactor. Equation 27 was applied to calculate the optimal gain for a pseudosteady-state filter employed to estimate the catalyst activity in the phthalic anhydride reactor. A temperature measurement is located at 5 = 0.125, and the actual value of the variance of the measurement-noise is u = 0.000 005. The variance of the measurement-noise is only used to simulate the measurement errors with a random-number generator, but it is assumed unknown to the filter designer. The equations for the pseudostatic scalar system are kI+l = ki + wi (28) yi = f(ki) + u; f is a nonlinear function relating the measurement and the parameter through eq 7-9. A nonlinear constant-gain filter was employed to estimate the parameter hi. Several values were used for the filter gain. Figure 7 shows the simulation results for the case of K = 1, which is lower than the optimal gain. The correlation covariances, Co and C1, were approximated from finite samples of the suboptimal filtering results as fin
Ind. Eng. Chem. Res., Vol. 26, No. 2, 1987 331 Table 111. Results of Filter Gain Adaptation C1 K* CO 0.4571 X 0.5428 X 0.5 0.3461 X 1 0.4390 X 0.5350 X 10”’ 4 0.1748 X 0.6195 X 10” 0.1542 X 6 0.1538 X lo4 -0.2956 X 6.36 -0.2656 X 0.1597 X 10”’ 8 -0.6411 X 0.1839 X 10 -0.6742 X 0.8770 X 14
poulos are gratefully acknowledged.
~
K
Nomenclature
6.56 6.54 6.36 6.36 6.38 6.44
A = linear system matrix CA = dimensionless concentration of A CB = dimensionless concentration of B CA* = concentration of A, kmol/m3 CB* = concentration of B, kmol/m3 CAI, = dimensionless inlet concentration of A CB,,= dimensionless inlet concentration of B CA*ln= inlet concentration of A, kmol/m3 CB*ln = inlet concentration of B, kmol/m3 C, = correlation of innovation sequence for time difference ,j
C, = estimate of correlation of innovation sequence for time difference j Cref = reference concentration, kmol/m3 cPr= fluid heat capacity, J/(kg K) c = catalyst heat capacity, J/(kg K) = activation energies, J/kmol; j = 1-3 f = nonlinear function defined by eq 28 H = sensitivity matrix h = nonlinear function defined by eq 11 h = scalar sensitivity i = discrete time step K = optimal steady-state filter gain for adaptive filter K = optimal steady-state filter gain K, = Kalman filter gain at time i K = optimal steady state filter gain for adaptive filter K* = suboptimal steady-state filter gain for adaptive filter k = parameter vector k,/,:l = parameter estimate at time i using measurements to time i - 1 A,/,, = parameter estimate at time i using measurements to time
6
81
n.m
1
1o.m
I
20
I 30.00
w
I 40.00
I 50 00
I
60 00
TJME
Figure 8. Estimation of catalyst activity, one temperature measurement a t [ = 0.125, filter gain K* = 6.36. l and 2 are the same as in Figure 2.
e,
e,
The values of and are substituted in eq 27, and the linearized version of eq 28 is used, i.e., h = df/dk The results are presented in Table 111. Suboptimal gains larger than 14 and smaller than 0.5 lead to unstable filters. For K* = 0.5 and 14, the adaptation algorithm did not converge. For suboptimal gains over the range from K* = 1to K* = 10, the algorithm found the optimal gain with good accuracy. The performance of the filter found through adaptation is presented in Figure 8. 7. Conclusion
It was shown that parameter estimation with a constant-gain, nonlinear filter is a feasible tool for packed-bed reactor modeling. The catalyst activity is estimated from several temperature measurements. A concentration measurement must be included in order to successfully estimate the heat-transfer coefficient and the catalyst activity. Furthermore, it was shown that if measurement- and process-noise variances are unknown, the optimal filter gain can be found from suboptimal filter data for a scalar system (one parameter, one measurement). More work is needed to establish the feasibility of adaptive estimation for simultaneous estimation of several reactor parameters. Acknowledgment
This paper is based on work performed at the University of Minnesota for the author’s Ph.D. thesis. Help and guidance from the thesis advisor Prof. George Stephano-
1
k*lll-l = suboptimal estimate at time i using measurements to time i - 1 k,,,T = suboptimal estimate at time i using measurements to time i k, = reaction rate frequency factor, l/s; j = 1-3 L = reactor length, m m = number of parameters n = number of measurements 0 = observability matrix p* = steady state value of P * , + ~ / , P*,+~,,= as defined in eq 22a Q = process-noise covariance matrix q = scalar process-noise R = measurement-noise covariance matrix R, = gas constant r,* = reaction rates, kmol/(m3 9); j = 1-3 r = scalar measurement-noise T = reactor radius, m s, = innovation sequence of optimal filter s,* = innovation sequence of suboptimal filter T = dimensionless temperature T , = dimensionless coolant temperature T,, = dimensionless inlet temperature T* = temperature, K T,* = coolant temperature, K T*,, = inlet temperature, K Tref= reference temperature, K t* = time, s U = overall heat-transfer coefficient, J/(m2 s K) U,,,, = nominal overall heat-transfer coefficient, J/(m2s K) D = fluid velocity, m/s u, = measurement-noise W = as defined in eq 22b w , = process-noise y , = measurement at time i z = reactor length coordinate Greek Symbols AH, = heat of reactions, j = 1-3 e
= void volume fraction
332 Ind. Eng. Chem. Res., Vol. 26, No. 2, 1987
dimensionless activation energies, j = 1-3 = dimensionless heat of reactions, j = 1-3 p = catalyst activity pi/, = estimate of at time i using measurements to time i yj = K~
Gi,,L-l = estimate of
at time i using measurements to time
i-1 v = normalized dimensionless heat-transfer coefficient ?i,i-l = estimate of v at time i using measurements to time i
-1 ?i,i = estimate of v at time i using measurements to time i $ j = dimensionless reaction rate constants, j = 1-3 pf = fluid density, kg/m3 p s = catalyst density, kg/m3 Z,,, = error covariance at time i using measurements to time 2
Zi,!-l = error covariance at time i using measurements to time 2-1
Ziji = sLeady-state error covariance matrix at time i using measurements to time i
2i,i-l= steady-state error covariance matrix at time i using measurements to time i - 1
hs = error covariance for weighted least squares estimators up = standard deviation of catalyst activity estimate uV = standard deviation of heat-transfer coefficient estimate
ugU*= T
5
covariance of parameter estimates
= heat capacity ratio = reactor length coordinate
Ro = nominal dimensionless heat-transfer coefficient R = dimensionless heat-transfer coefficient w
= correlation coefficient
Appendix A. Numerical Values of the Parameters of the Phthalic Anhydride Reactor L=4m t
= 0.35
k,’ = 2.418
X
lo9 s-l
k2’ = 6.706
X
lo9 s-’
k,’ = 1.013 X
lo9 s-]
E , = 1.129 X
los J/kmol
E z = 1.313 X
lo8 J/kmol
E , = 1.196 X
los J/kmol
AH, = -1.285
X
lo9 J/kmol
AHz = -3.276
X
lo9 J/kmol
AH, = -4.561
X
lo9 J/kmol
Cref= 0.000 1811 kmol/m3 Tref= 628 K ps
= 2000 kg/m3
pf =
0.582 kg/m3
c p , = 836 J/(kg K)
cpr= 1045 J/(kg K) F =
0.135 m
Literature Cited Anderson, B. D.; Moore, J. B. Optimal Filtering; Prentice Hall: Englewood Cliffs, NJ, 1979. Bonvin, D.; Rinker, R. G.; Mellichamp, D. A. Chem. Eng. Sci. 1980, 35, 603. Bonvin, D.; Rinker, R. G.; Mellichamp, D. A. Chem. Eng. Sci. 1983, 38, 233, 245. Box, C. E. P.; Lucas, H. L. Biometrika 1959, 46, 77. Carberry, J. J.; White, D. Ind. Eng. Chem. 1969, 61, 7. Carew, B.; Belanger, P. R. IEEE Trans. Autom. Control 1973, AC18, 582. Clough, D. E.; Ramirez, W. F. Presented at the 72nd AIChE Annual Meeting, San Francisco, 1979. Cutler, C. F.; Ramaker, B. L. Presented at the 86th AIChE Annual Meeting, April 1979; also in Joint Automatic Control Conference Proceedings, San Francisco, 1980. Draper, N. R.; Hunter, W. G. Biometrika 1966, 53, 525. Fisher, D. G.; Seborg, D. E. Multivariable Computer Control, A Case Study; North Holland-American Elsevier: New York, 1976. Froment, G. F. Ind. Eng. Chem. 1967,59, 2, 18. Garcia, C.; Morari, M. I n d . Eng. Chem. Process Des. Dev. 1982,21, 308.
Gavalas, G. R.; Seinfeld, T. H. Chem. Eng. Sci. 1969, 24, 625. Georgakis, C.; Aris, R.; Amundson, N. R. Chem. Eng. Sci. 1977, 32, 1359 (three parts). Goldmann, S. F.; Sargent, R. W. H. Chem. Eng. Sci. 1971,26,1532. Heffes, H. IEEE Trans. Autom. Control 1966, AC-11, 541. Jutan, A.; Tremblay, J. P.; MacGregor, J. F.; Wright, T. D. AIChE J. 1977, 23, 732. Kalman, R. E. J. Basic Eng. 1960, 82, 35. Kwakernaak, H.; Sivan, R. Linear Optimal Control; Wiley-Interscience: New York, 1972. Martin, G. D. AIChE J . 1981, 27, 748. Mehra, R. K. IEEE Trans. Autom. Control 1970, AC-15, 175. Nishimura, T. IEEE Trans. Autom. Control 1967, AC-12, 268. Richalet, J. A.; Rault, A.; Testud, J. D.; Papon, J. Automatica 1978, 14, 413.
Ind. Eng. Chem. Res. 1987,26, 333-336 Silva, J. M.; Wallman, P. H.; Foss, A. S. Ind. Eng. Chem. Fundam. 1979, 18, 383. Stewart, W. E.; Sorensen, J. P. Paper presented a t the Proceedings of the 5th European Symposium on Chemical Reaction Engineering, Amsterdam, 1972. Swerling, P. J. J. Astronaut. Sei. 1959, 6 , 46.
333
Wong, C.; Bonvin, D.; Mellichamp, D. A.; Rinker, R. G. Chem. Eng. Sei. 1983, 38, 619.
Received for review March 25, 1985 Revised manuscript received February 13, 1986 Accepted March 28, 1986
Cooling of a Stretching Sheet in a Viscous Flow Binay K. Dutta Chemical Engineering Department, Calcutta University, Calcutta 700 009, India
Anadi S. Gupta* Mathematics Department, Indian Institute of Technology, Kharagpur 721 302, India
A hot, flat sheet issues from a thin slit into a n incompressible viscous fluid and is subsequently stretched in its own plane a t a velocity proportional to the distance from the slit. T h e cooling of the sheet in this extension zone is studied by solving the coupled heat-transfer problem in the fluid and the sheet. Variation of the sheet temperature with distance from the slit is found for several values of the Prandtl number and stretching speeds. I t is shown t h a t for a fixed Prandtl number, the surface temperature decreases with a n increase in the stretching speed. 1. Introduction Boundary-layer flow on moving solid surfaces is an important type of flow arising in a number of technical processes. This problem was first studied by Sakiadis (1961). The momentum and heat and mass transfers from such surfaces are determined by the structure of this layer. Due to the entrainment of the ambient fluid, this boundary layer is different from that in the Blasius flow past a flat plate. Erickson et al. (1966) extended this problem to the case when the transverse velocity a t the moving surface is non-zero and is such that similar solutions exist. These investigations have a bearing on the problem of a polymer sheet extruded continuously from a die and is based on the tacit assumption that the sheet is inextensible. However, situations may also arise in the polymer industry where one has to deal with flow over a stretching sheet. For example, in a melt-spinning process, the extrudate from the die is generally drawn and simultaneously stretched into a filament or sheet, which is then solidified through rapid quenching or gradual cooling by direct contact with water or chilled metal rolls. In fact, stretching imparts a unidirectional orientation to the extrudate, thereby improving its mechanical properties (Fisher, 1976). Several metallurgical processes involve the cooling of continuous strips or filaments by drawing them through a quiescent fluid, and in the process of drawing, these strips are sometimes stretched. Mention may be made of drawing, annealing, and tinning of copper wires. In all these cases, the properties of the final product depend to a great extent on the rate of cooling. McCormack and Crane (1973) gave a similar solution in a closed analytical form for steady two-dimensional, boundary-layer flow over a sheet which issues from a slit and is subsequently stretched at a velocity proportional to the distance from the slit. The corresponding nonsimilar solution for the same problem in the presence of a uniform free stream velocity was obtained by Danberg and Fansler (1976). Gupta and Gupta (1977) analyzed heat and mass transfer in the boundary layer over a stretching sheet held at constant temperature, the sheet being subject to suction or blowing. Similar studies for the temperature distribution in the flow a viscous, incompressible fluid caused by the stretching of an impermeable, nonporous sheet were made by Crane (1970) for prescribed surface temperature variation and 0888-5885/81/2626-0333$01.50/0
by Dutta et al. (1985) for uniform surface heat flux. On the other hand, Vleggaar (1977) investigated the cooling of a sheet that issues from a slit and is stretched at a velocity proportional to the distance from the slit by solving the conjugate problem of heat transfer inside the sheet and the fluid. However, this study suffers from the serious defect that the temperature distribution in the boundary layer was found for a constant surface temperature of the sheet, which was subsequently used to determine the longitudinal variation of the same surface temperature. The aim of this paper is to remove this anomaly through a rigorous analysis of the above conjugate heat-transfer problem. The effects of the Prandtl number and the stretching rate on the surface temperature are also investigated. 2. Mathematical Formulation Consider a broad, flat sheet issuing from a long, thin slit at x = 0 and y = 0 and, subsequently, being stretched a t a velocity proportional to the distance from the slit as in Figure 1. Assuming the boundary layer approximations, the equations of continuity, momentum, and energy in the steady flow near the sheet are in usual notations
(3) where the effects of natural convection and viscous dissipation are neglected in (3) and the fluid properties are assumed to remain constant over the temperature range of interest. The sheet material (of density p and specific heat C,) gradually cools down along the direction of motion, and a heat balance equation in the sheet gives the temperature, T2,of the sheet as the solution of
0 1987 American Chemical Society
