938
I n d . Eng. C h e m . Res. 1987,26, 938-944
Dynamic Estimation of Bubble Parameters in a Fluidized Bed Subjected to Load Disturbances Dale C. Gyure* and David E. Clough Department of Chemical Engineering, University of Colorado, Boulder, Colorado 80309
Bubble frequency and bubble velocity in a fluidized bed are estimated from the cross-correlation function of pressure measurements during dynamic changes in fluidizing conditions. Estimation is based on a theoretical form of the cross-correlation function, a mathematical function of bubble frequency and velocity. The two estimation algorithms used in this work are sequential weighted least squares and a variation of Kalman filtering. The usefulness of both algorithms is demonstrated by measuring bubble frequency and bubble velocity in real time during dynamic changes in bed inventory and particle size distribution. Both algorithms successfully follow sudden changes in bubble parameters, and calculated values of bubble frequency and velocity are generally within 10% of values predicted by other methods. The most salient feature of a bubbling fluidized bed is its random upward movement of bubbles through the dense phase. Bubble-phase properties have a significant impact on fluidized bed performance, and bubble parameters such as frequency, velocity, size, and population density reflect the overall condition of both the bubble and the dense phases. Bubble parameter estimation during fluidized bed operation is of possible industrial interest as a means of tracking changes in fluidized bed hydrodynamic behavior. Changes in performance can be caused by unmeasurable disturbances in particle size or phase density or by measurable changes in air flow or working pressure. Fortunately, measurement of bubble-phase behavior via pressure measurements requires simple instrumentation commonly found or easily installed on most fluid beds. On-line analysis of bubbles is also valuable in the laboratory because bubble parameter estimation is equivalent to estimation of fundamental hydrodynamic properties of the bed that are not directly observable. The desirability of on-line bubble parameter estimation is evident in recent fluidization research. Robertson (1979) reported a narrow range of acceptable gas flow rates in a small-diameter, high-density fluid bed. Maintaining any degree of fluidization in the bubbling regime was impossible without identification of fluidization quality via pressure fluctuation amplitude (related to bubble size). Fan et al. (1981) investigated the effects of gas velocity, bed height, particle size, and distributor design on pressure fluctuation frequency and amplitude. They suggested that statistical analysis of pressure fluctuations could be used to monitor changes in bed conditions or to diagnose abnormalities in bed operation. Ligon and Amundson (1981a,b) modeled the stochastic nature of bubbles in an exothermic reacting fluidized bed with a fluctuating mass-transfer coefficient. Computer simulations suggested that conversion could be maximized by adjusting gas input to maintain small bubbles or low-variance,high-frequency pressure fluctuations. Monitoring applications are implied by Wen and Chen (19821, who noted fluctuations in entrainment rate corresponding to surface eruptions of bubbles, and Kennedy et al. (1981),who observed that tube bank stresses were experimentally proportional to bubble size. Bubble size measurement could be used to control entrainment or lengthen tube bank life. Previous workers have measured bubble-phase properties by using a statistical analysis of pressure fluctuations. *Present address: Coors Biotech Products Co., Ft. Collins, CO 80525.
0888-5885/87/2626-0938$01.50/0
Swinehart (1966) calculated the ”correlation average propogation velocity” off-line from the cross-correlation function of two sets of pressure measurements taken at the wall of a fluidized bed. More recently, Sitnai and Dent (1981) and Sitnai (1982) have shown that the cross-correlation function of four pressure measurements can be used to compute bubble frequency, velocity, rise time, and spacing. Data were collected from an incipiently fluidized bed injected with air to create singly rising noninteracting bubbles. Analysis of the data was completely off-line. Weimer and Quarderer (1986) employed Sitnai’s methods to compute bubble velocity and infer bubble size in highpressure fluidized beds. Fan et al. (1983) computed bubble velocities in real time from the “delay time” in the crosscorrelation function by using a correlation and probability analyzer. In like fashion, Ho et al. (1983a) studied the onset of slugging in a gassolid fluidized bed. Slugging was judged by bubble rise velocity and computed from the cross-correlation function by using Fan’s time-delay method. Ho et al. (1983b) and Yutani et al. (1983) have even used the autocorrelation function to study the periodicity or nonperiodicity of capacitance fluctuations near the distributor plate. In no cases have previous workers employed any technique for computing velocity other than repeated measurement of the first maximum in the cross-correlation function. Furthermore, frequency and velocity have never been computed simultaneously in real time. This work goes substantially beyond both batch and on-line calculations of single bubble parameters by presenting a method for computing both bubble frequency and velocity in real time. The methods were derived from advanced estimation theory and required simple instrumentation and reasonable real time computing capability.
Development of Estimation Algorithms Estimation of bubble frequency and velocity requires a functional model relating bubble parameters to output variables that are readily measurable during fluidized bed operation. The markedly periodic behavior of a laboratory fluidized bed and previous work suggested use of the cross-correlation function for measuring bubble parameters. Because the cross-correlation function is a statistical measure of the temporal relationship between two stationary time series, it naturally measures bubble velocity from two sets of pressure measurements taken at two positions parallel to the vertical axis of the bed. When coalescence no longer dominates bubble-phase hydrodynamics, pressure fluctuations due to rising bubbles gen0 1987 American Chemical Society
Ind. Eng. Chem. Res., Vol. 26, No. 5, 1987 939 erate approximately periodic pressure series a t each measurement location. The two series are displaced in time by an amount corresponding to the average bubble rise velocity. The cross-correlation function achieves its maximum value at a lag directly related to the time displacement of the two series. In the equation for the cross-correlation function,
collectively referred to as the state. x1 = nonideal amplitude of real sinusoidal cross-correlation functions; qualitatively measures the variance of bubble frequency and velocity (5a) x 2 = 27rfT, (rad)
(5b)
4 (phase lag, rad) In state variable vector notation, eq 2 becomes rxy(k)= h k ( ~ = ) X I COS ( k ~ 2~ 3 ) x3 =
X,and Yt are stationary input and output series, and E is the expectation operator (yielding expected value). a is the square root of series variance, w is the series mean, and the independent variable, k, is lag measured in integral sample intervals for discrete series. Suppose that two sinusoidal functions have frequency, f, but are separated by a phase lag, 4. If the two functions are sampled every T, seconds to give an input series, X,, and an output series, Y,, the theoretical cross-correlation function includes frequency and phase lag as independent variables.
r&)
= A cos (27rfkT, - 4)
(2)
The amplitude, A , is equal to unity, and T, is the discrete sample time. Because these sinusoidal series represent real pressure series caused by rising bubbles, the frequency in eq 2 represents the dominant bubble frequency between the two pressure taps. Furthermore, the phase shift, 4, is easily related to the time delay, 7, of the output wave relative to the input wave by (3)
The velocity of the bubble phase between the two taps is easily calculated Ubub
= D/7
(4)
where D is the known distance between the two pressure taps. This calculation is entirely consistent with the methods of Swinehart (1966), Oki et al. (1977), and Fan et al. (1983). Equation 2 is sinusoidal with an amplitude of unity because of the exact sinusoidal functions describing the input and output pressure series. These exact functions are used to derive a functional model of the cross-correlation function containing the parameters frequency, phase lag, and correlation lag. In a real system, the cross-correlation amplitude is always less than unity and decays as lag increases because of the stochastic nature of bubble frequency and velocity and because each pressure series is only approximately periodic. Real cross-correlation functions can, however, be strongly sinusoidal and similar to the theoretical form. Furthermore, the cross-correlation function consistently retains its sinusoidal character even for nonsinusoidal periodic input and output series. Gyure (1984) and Gyure and Clough (1986) have successfully used the cross-correlation function to calculate the frequency and phase lag of nonsinusoidal periodic waveforms and pressure waveforms caused by real bubbles.
Estimation by Sequential Weighted Least Squares With eq 2, the problem of estimating bubble frequency and velocity from sequential measurements of the crosscorrelation function can be stated as an optimal estimation problem. In this case, a three-component state vector is defined such that all state variables are dimensionless and of the order of unity. The elements of the vector, x,are
Bubble frequency, fbub, and velocity, quantities fbub
=
x2
2.lrT,
Ubub,
(bubbles/ s)
(5c) (6)
are derived (74
where D is the distance (30.5 cm) between the two pressure taps. Note that x1 does not appear in the equations for frequency and velocity even though it is needed to model real cross-correlation functions. The three-component state vector, x,was estimated at steady state using a measurement vector, y , containing random errors, v , independent of the state. The crosscorrelation functions at lags -8, -4, 0, +4, and +8 were chosen as the elements of the measurement vector, y , because rxy(k)deviated on the average least from a perfect sine wave at those lags. Mean zero measurement errors were needed for unbiased estimates of the state. The nonlinear model relating the cross-correlation function to the state is written for the specific lags of interest as [ yr (x-y8(]- 4 )
y = rxy(0) rxy(+4) r x y ( t 8)
=
C; ~$1[; I X, COS cos
(-434 ( - 8 ~ 2-x3)
xi c o s ( O x , - ~ ~ ) t x 1 cos ( t 4 x , 2 x3) X , COS (+ 8 ~ -, x,)
=b(x)+ v(8)
u,
u , ~ u,,
In order to compute a weighted least-squares estimate of the state, x, the following common quadratic error function must be minimized with respect to x. 1 J = -[(x - %)'*M-'(x - a) + (y - h ( ~ ) ) ~ - R --' hb ( ~ ) ) ] 2 (9) % is a state estimate prior to measurement. The weighting matrices are the state error covariance matrix, M, and the measurement error covariance matrix, R. The quadratic error function penalizes new estimates of x when they are both far from the previous estimate and yield poor predictions, h (x),for the measurements, y . The matrices, M and R,used in this work are given in eq 10. Note that 2.43 -0.02 0.26 0.16 0.27
1
0.16 0.27 2.74
0.000 0.298 -0.034 0.577 -0.067 0.346 -0.067 0.542 -0.128 0.346 -0.128 0.572 0.076 0.195 0.273 0.046 0.563 0.000
off-diagonal elements in both matrices were generally nonzero. A methodology for determining these matrices for any fluid bed can be found in Gyure (1984) and Gyure and Clough (1986). For a fluidized bed operating within
940
Ind. Eng. Chem. Res., Vol. 26, No. 5, 1987
Table I. Comparison of Bubble Frequency and Bubble Velocity Estimated before and after Changes in Fluidizing Conditions with Bubble Freauency and Bubble Velocity Computed by Other Methods frequency, s-'," frequency, s-l,* frequency, S - I , ~ velocity, cm/s,d velocity, cm/s from density theoretical time delay (this work) fluctuations prediction method exptl condition 0.97 f 0.04, 140 f 7 0.86 (-11%) 1.03 (+6%) 380 f 55 (+170%) blend C at 4.7Umf;H = 76 cm, Figure 5 1.08 (-8%) 1.15 (-2Yo) 127 f 12 (+7%) 1.18 f 0.05, 119 f 6 blend C at 4.7Ud; H = 61 cm, Figure 5 0.86 (-12%) 1.03 ( + 5 % ) 380 f 55 (+170%) 0.98 f 0.01, 139 f 3 blend C at 4.7Ud; H = 76 cm, Figure 6 1.08 (-8%) 1.15 (-2%) 127 f 12 ( + 5 % ) 1.17 f 0.02, 121 f 3 blend C at 4.7Ud; H = 61 cm, Figure 6 1.18 f 0.04, 131 f 8 1.18 (0%) 135 f 20 (+3%) blend D (before mixing); H = 61 cm, Figure 7 1.41 f 0.05, 149 f 5 1.27 (-10%) 138 f 5 (-7%) final blend (after mixing); H = 53 cm, Figure 7 "Frequency and velocity are reported as the average estimated value plus or minus the standard error, u. bBubble frequency was estimated from the dominant density fluctuation frequency measured 38 cm (15 in.) above the distributor plate. Detailed methods and calculations are given by Gyure (1984). 'Bubble frequency was estimated from the predictive equation of Verloop and Heertjes (1974), evaluated for 6 = 0.47, blend C, and 0.45 in the mixing experiment. 1 g2-€'12
[ I
f=zK T
The maximum locus of the cross correlation function, k,,, functions were measured earlier under the same conditions.
is computed by second-order polynomial interpolation. Cross-correlation
narrow steady-state limits, the matrics, M and R, need to be determined only once. Weighted least squares is ideally suited for computing a single optimal estimate of the state, x, given a priori information. Furthermore, the algorithm is easily adapted to a nonlinear model such as eq 8, according to Bryson and Ho (1975), for example. In this work, however, the weighted least-squares algorithm was adapted for sequential use by defining the state estimate prior to measurement as the previous estimate. That is, x = xi-1. The final form of the estimator given below is ideally suited for on-line computation. Partial justification for this extension of weighted least squares is provided by comparing the state error covariance matrix, M, and the covariance of estimated state variables, computed sequentially from real cross-correlation functions and eq 11. In general, the covariance of x,estimated by eq 11, was no larger than the a priori state error covariance, M. Further justification for this extension of weighted least squares can be seen in Table I, which compares average estimated bubble frequency and velocity to frequency and velocity computed by other methods. xi+l = xi
+ PiHiTR-l(yi+l - h (Xi))
pi-l = M-1
+ HiTR-lH,
(1la) (1lb)
ki+lis a state vector prediction after transition to stage i + 1 but before measurement in stage i + 1. Optimal estimates of the state after measurement in stage i + 1 are given by the Kalman filter. Xi+l
=4+1
+ Pi+1Hi+lTR-'(Yi+1- h(ki+l)) (12b)
P 1.+ 1 -l = Mi+l-l + Hi+lTR-lHi+l
(124
M ~ =+ Q~ ~ P +~riQiriT @ ~ ~
(12d)
Qi is the covariance matrix of the forcing vector, wi, and M and R are again the state and measurement error covariance matrices. In order to derive a second estimation algorithm for a bubbling fluidized bed operating at steady state, the assumption of stationary states is applied to the Kalman filter. That is, since average bubble frequency, velocity, and cross-correlation function amplitude are constant while operating at steady state, Qi = I, the identity matrix. Furthermore, in this variation of Kalman filtering, we neglect the influence of w on the states by setting Q equal to zero. Note that the assumption of stationary states leads to an estimation algorithm containing no process model. Also, no Q matrix needs to be specified. Thus, the second estimation algorithm is given by
xi+,= xi + Pi+lHi+lTR-l(yi+l - h(xi)) ax, ah-, Hi =
ax, ah, ax,
ax, ah.,
ax,
P 1.+ 1 -l = P[l
ah_,
ax, ah, ax,
ax,
ax,
ax,
ak, ax, at xl
Variation of Kalman Filtering with Exponential Forgetting A second estimation algorithm for the state vector, x, is derivable from the general Kalman filter. For the general case of nonstationary states driven by a known stochastic forcing function, state variable dynamics are given by 2i+1 =Q~X + ~riwi (124
+ Hi+lTR-lHi+l
(13a) (13b)
M and R are given by eq 11,and Hi+l is given by eq llc, again evaluated a t xi. Note the difference in subscript notation between the weighted least squares and the Kalman filter. This difference arises because weighted least-squares matrices are referenced to the past state estimates, and Pi and Hi are functions of xi. Kalman matrices are referenced to the stages, and Pi+land Hi+l are functions of ki+l. For stationary states, ki+lequals xi, and Pi+l and Hi+lare functions of xi. In practice, eq 13 gives rise to a zero-converging P matrice representing growing certainty in the estimated states with each additional measurement. This behavior is a natural consequence of removing the noise input, wi, from the stationary states model. Unfortunately, state estimation using a zero-converging P matrix leads to quiet estimates but exceedingly poor response to real changes in the states caused by a change in fluidizing conditions. To avoid this, exponential weighting of past data was in-
Ind. Eng. Chem. Res., Vol. 26, No. 5, 1987 941
....---.--._._-
4
0
0
10
20
0
IO
o
IO
Sequential satmate
IO
20
(*)
Of
Figure 1. State variables x 2 and x 3 estimated by Kalman filtering during a simulated step change in bed height. Estimator response is shown as a function of the forgetting factor, c2.
corporated into the version af Kaknan filtering above. Paralleling the development of recursive least squares given by Eykhoff (1974), we write
+ Hi+lTC-2R-1Hi+l0 < c2 I1
55-
z3 v3
0 -c z
(14)
5f
c2 is called the forgetting factor, and the kth past measurement contributes ck to the estimate of the current state. Note that for c2 = 1, all past measurements are equally weighted. Conversely, for c2
