Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Time series analysis of voidage signals in a bubbling fluidised bed Shen, Chengyu 1999

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1999-0442.pdf [ 4.89MB ]
JSON: 831-1.0058937.json
JSON-LD: 831-1.0058937-ld.json
RDF/XML (Pretty): 831-1.0058937-rdf.xml
RDF/JSON: 831-1.0058937-rdf.json
Turtle: 831-1.0058937-turtle.txt
N-Triples: 831-1.0058937-rdf-ntriples.txt
Original Record: 831-1.0058937-source.json
Full Text

Full Text

T I M E SERIES A N A L Y S I S OF V O I D A G E S I G N A L S IN A B U B B L I N G FLUTDISED B E D by C H E N G Y U S H E N B.Eng., Nanjing University of Chemical Technology, 1987 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF M A S T E R OF A P P L I E D S C I E N C E in T H E F A C U L T Y OF G R A D U A T E STUDIES D E P A R T M E N T OF C H E M I C A L E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF B R I T I S H C O L U M B I A J U N E 1999 © Chengyu Shen, 1999 jr u i_ — ': s — -3 -5« • T H u , 1 <+ : e 1 A C T — C O M ,• / p . o in presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, ) agree that the Library shall make it, freely available for reference and study. 1 further agree * that permission for. extensive . . .copying of this thesis for scholarly purposes may • be< granted by the head . of my ! department or . by . his or her representatives. It : ts . understood - that copying or •'.''publication of this thesis for financial gain shall not be allowed . without my written . permission. ' ' ' Department of QAM^^C £ ^ > Q & ^ ^ The University of British Columbia Vancouver, : Canada — — DE-6 (2/88) Abstract The dynamics of voidage signals in bubbling fluidised beds are investigated both experimentally and numerically. Experimental voidage signals were obtained by using an optical fibre probe at different positions in a freely bubbling fluidised bed of diameter of 150 mm. Numerical simulated voidage signals were also predicted using the Cliff and Grace bubble coalescence model. Analysis of the experimental signals demonstrates that deterministic chaos can arise in bubbling fluidised beds. Voidage signals collected at different geometrical positions indicate that bubble size and frequency influence the chaotic characteristics o f the voidage signals. In the bottom, the voidage signals were generally more chaotic than in the top section, while near the bed surface, with voidage signals becoming less chaotic because of the bubble coalescence which has already occurred. The numerical study used the three-dimensional version of the bubble coalescence model to simulate optical probe signals in a gas-solid fluidised bed. These predictions demonstrate that deterministic chaos can arise from non-linear bubble interactions in bubbling fluidised beds. It is found that both bubble frequency and initial bubble size influence the chaotic characteristics of the simulated voidage fraction signal, with bubble frequency having a stronger effect than the initial bubble size. II T A B L E O F C O N T E N T S Abstract JJ Table of Contents III List of Tables V I List of Figures V U Acknowledgement I X Chapter 1 Bubbling Fluidised Beds and Their Properties 1 1.1. Bubbles in Fluidised Beds 1 1.2. Simulation of Bubbling Fluidised Beds Hydrodynamics 4 1.2.1 Bubble Interactions 4 1.2.2 Distribution of Gas between Phases 6 1.2.3. Maximum Stable Bubble Size 7 1.3. Conclusion 8 Chapter 2. Introduction to Deterministic Chaos 9 2.1 Linear and Non-linear Systems 9 2.2. Non-Linearity and Chaos 10 2.2.1 A Model of Convecting Fluids: The Lorenz Model 12 2.2.2.Logistic M A P (Devaney, 1989). 16 2.3. Reasons for Quantifying Chaotic Behaviour 16 2.4. Time Series Analysis 18 2.4.1. Time History 19 2.4.2. Strange Attractor 19 2.4.3. Autocorrelation Function 21 III 2.4.4. Fractal Dimension 22 2.5. Previous Application of Deterministic Chaos Theory to Fluidised Beds 30 2.5.1. Fluidisation Regime Identification 31 2.5.2. Characterisation of Hydrodynamic Structure 32 2.5.3. Fluidisation Quality Control 34 2.5.4. Fluidised Bed Scale-Up 35 2.6 Scope of Thesis 37 Chapter 3. Deterministic Chaos Analysis of Experimental Voidage Signals in a Freely Bubbl ing Fluidised Bed 40 3.1 Introduction 40 3.2. Experimental Equipment and Probe Calibration 41 3.3 Experimental Voidage Distribution for a Freely Bubbl ing Fluidised Bed 48 3.4 Deterministic Analysis of Optical Fibre Signals 49 3.4.1. Analysis Conditions 49 3.5 Determinstic Chaos Analysis o f Re-Constructed Optical Fiber Signals 53 3.6 Summary 54 Chapter 4 Bubble Coalescence Simulation Based on Clift and Grace Model and Dynamic Behaviour Analysis 61 4.1 Introduction 61 4.2 Model l ing of Three-Dimensional Bubbl ing Fluidised Bed 61 4.2.1 Coalescence 67 4.2.2 Numerical Solution 67 4.3 Signal Simulation 69 IY 4.4. Chaos Analytical Method 70 4.4.1 Phase-Space Plots 70 4.4.2 Lyapunov Exponent 70 4.4.3. Autocorrelation Function 70 4.4.4 Correlation Dimension 71 4.5 Sample Results and Discussion 71 Chapter 5. Conclusions and Recommendations 5.1 Conclusion 5.2 Recommendations Nomenclature References Appendix 3-Dimension Fluidised Bed Simulation Fortran Computer Code 78 78 79 80 82 89 V L I S T O F T A B L E S Table 2.1 Relative Advantages and Disadvantage of Voidage and Pressure Measurements 39 Table 3.1. Optical Fibre Probe Sample Locations in the Free Bubbl ing Fluidised Bed and Chaotic Properties 51 Table 3.2. Lypunov Exponent o f Optical Fibre Signals Collected at Different Locations in a Fluidised Bed 53 Table 4.1 Interaction Coefficients for the Prediction of Bubble Mot ion in Three Dimensions (adapted from Johnsson's thesis, 1973) 66 Table 4.2 Summary o f results o f numerical simulations using three-dimensional bubble coalescence model 77 VI L I S T O F F I G U R E S 1 Figure 1.1 F low patterns in gas-solids fluidised beds Figure 2.1 Lorenz attractor 15 Figure 2.2 Bifurcation diagram of logistic difference equation 17 Figure 3.1 Optical fibre probe calibration equipment proposed by Issangya et. al. (1996) 43 Figure 3.2 Calibration curve for optical fibre probe 44 Figure 3.3 Optical fibre probe signal at the height of 110 mm 45 Figure 3.4 Optical fibre probe signal at the height of 220 mm 46 Figure 3.5 Optical fibre probe signal at the height of 55 mm 47 Figure 3.6 Radial voidage profile at U=0.2 m/s for different height 48 Figure 3.7 Autocorrelation function for case OM1 51 Figure 3.8 Phase space plot for case OM1 52 Figure 3.9 Correlation integrals for optical fiber probe signal sampled under freely bubbling conditions at a sampling frequency of 200 H z 54 Figure 3.10 Autocorrelation function for superimposed signals o f Lorenz attractor and Logmap 55 Figure 3.11 Correlation integrals for a superimposed signal consisting o f Lorenz attractor and logistic map, a simple one dimensional map (time delay =1) 56 Figure 3.12 Correlation integrals for a superimposed signal consisting of Lorenz attractor and logistic map, a simple one dimensional map (time delay =10) 57 VII Figure 3.13 Correlation integrals for a superimposed signal consisting of Lorenz attractor and logistic map, a simple one dimensional map (time delay =10) 58 Figure 3.14 Correlation integrals for reconstructed optical fibre signal sampled under freely bubbling conditions, assuming all dense phase voidages equal to the voidage at minimum fluidisation, sampling frequency = 200 H z (time delay =10) 59 Figure 3.15 Correlation integrals for reconstructed optical fibre signal sampled under freely bubbling conditions, assuming all bubble phase voidages equal to 1, sampling frequency = 200 H z (time delay =15) 60 Figure 4.1 Phase space plot of the simulated voidage fraction signal 73 Figure 4.2 Correlation integrals for simulated voidage signal . 74 Figure 4.3 Correlation integrals for simulated voidage signal 75 Figure 4.4 Correlation integrals for simulated voidage signal 76 VII I A C K N O W L E D G E M E N T S I wish to thank Dr. John R. Grace for his support and supervision of my work. I would l ike to thank Dr. Bruce D. Bowen for his kind suggestion and help. I also wish to thank Dr. Junping Zhang for his assistance on many occasions. I finally wish to thank my wife, May Zhou for her consistent encouragement and support. I X Chapter 1 Bubbling Fluidized Beds and Their Properties The introduction of gas from the bottom of a column containing solid particles via a gas distributor can cause the particles to be fluidized. As shown in Figure 1.1, with increasing gas velocity, several f low patterns/regimes have been identified (i.e., fixed bed, delayed bubbling or particulate fluidization, bubbling fluidization, slugging fluidization, turbulent fluidization, fast fluidization and dilute pneumatic conveying). PNEUMATIC CONVEYING pELAYED ' V " BUBBLING AGGREGATIVE FLUIDIZATION' INCREASING U, < Figure 1.1 Flow patterns in gas-solids fluidized beds (Grace, 1986) 1.1. Bubbles in Fluidized Beds In a bubbling gas fluidized bed, the superficial gas velocity is between the minimum bubbling velocity and the minimum slugging velocity. A dense bubbling fluidized bed 1 contains regions of low solids density, usually called voids or bubbles. The surrounding region of higher density is called the emulsion or dense phase. Voids form at or near the distributor, grow mostly by coalescence, and rise to the surface. The top surface of a bubbling fluidized bed is well defined, with bubbles breaking through periodically. There are irregular pressure fluctuations of appreciable amplitude, caused by bubble formation, passage and bursting, as well as by surging of the bed surface. In many ways, a bubbling fluidized bed behaves like a bubble column where gas bubbles rise through a liquid of low viscosity. The shapes of bubbles are generally spherical-cap as in liquids where surface tension forces are relatively small. For both systems, bubbles tend to coalesce to give large bubbles. The interaction of a chain of bubbles leads to a higher rise velocity in both cases. Wal l effects act to retard bubbles in fluidized beds as in liquids. Experiments in bubbling fluidized beds indicate that most gas in excess of that needed to just fluidize the bed passes through the bed as bubbles and that the emulsion phase remains close to that at minimum fluidizing conditions. Unlike gas-liquid systems, there is convective flow(often called throughflow) of gas through the bubbles in fluidized beds, this being made possible by the permeable nature of the dense phase and the pressure gradient across the bubbles. The rise velocity of single bubbles in fluidized beds has been measured experimentally ( e.g. Davidson et. al., 1959; Harrison and Leung, 1962; Rowe et. al., 1965; Toei et. al., 1965) and is often expressed as ubr = on\\(gdb)m ( i ) 2 where db is the diameter of sphere having the same volume as the spherical-cap bubble and where wall effects are negligible. The rise velocity of bubbles depends on the same factors as for spherical-cap bubbles in liquids. For bubbles in freely bubbling fluidized bed, the bubble velocity can be expressed approximately by Ub = Uo ~Umf + Ubr (2) where the (u 0 -u„ i f ) accounts approximately for bubble interactions. As a first approximation the emulsion is often treated as a liquid of low or negligible viscosity. The so-called "two -phase theory of fluidization" is often commonly assumed, which says that the excess gas above that required for minimum fluidization travels through the bed as bubbles. The voidage of a bed, not counting the bubbles, remains close to Smf. A t minimum fluidizing conditions, the solids are nearly quiescent. At higher gas velocities the rising bubbles cause churning and mixing of solids. Davidson and Harrison (1963) developed a model which accounted for the movement of both gas and solids and the pressure distribution about rising bubbles in two-and three-dimensional beds. Collins (1965) and Stewart (1968) assumed kidney-shaped bubbles with indented bases with properties otherwise similar to those of the bubbles studied by Davidson and Harrison (1963). Jackson (1971) retained the spherical bubble, but allowed the voidage of the emulsion phase to vary. Murray (1965) developed a somewhat more rigorous model. Stewart (1968) summarized and commented on these studies. 3 The particle motion generated by bubbles causes rapid solids mixing, temperature uniformity and high heat transfer. These are valuable properties of industrial gas fluidized beds. Thus the distribution of bubbles in the fluidized beds is of considerable practical interest. It is well known, particularly for Geldart (1973) Group B materials, that bubble coalescence causes bubbles to migrate inward from the walls of the bed, and to set up regions of high bubble activity associated with characteristic particle circulation patterns. The origins of the non-uniformity in bubble flow have been identified and explained by Grace and Harrison (1968). 1.2. Simulation of Bubbling Fluidized Beds Hydrodynamics In order to simulate bubble distribution in fluidized beds, it is essential to consider not only the particle motion around a single bubble and a small number of interacting bubbles, but also the behavior of a large number of bubbles over a period of time which is much longer than the inverse of the bubble frequency. As well, some features are identified as particularly important, such as bubble interaction and coalescence, the size and frequency of bubbles formed at the gas distributor, the maximum stable bubble size and bubble splitting. 1.2.1 Bubble interactions Bubble interactions are of key importance in determining such factors as bubble growth, interphase mass transfer, and lateral motion of particles and gas. To a first approximation, the motion of each bubble in an interacting pair can be obtained by superposition, taking the vector sum of the rise velocity in isolation and the velocity at the position of the nose caused by the other bubble (Clift & Grace 1970, 1971). For relatively short-range interactions, the velocity field caused by each bubble, outside its own boundary and its wake, can be represented as a first approximation by a doublet. The : wake is assumed to be stagnant relative to the bubble with which it is moving. These approximations lead to a reasonable description of the trajectories of coalescing bubbles and of the boundaries of the region within which a smaller following bubble overtakes a larger one ( Clift & Grace 1971; Clift & Grace 1974). This approach has been extended to bubble chains (Clift & Grace 1972) and to freely bubbling beds ( Johnsson et. al. 1974; Nguyen et. al. 1976; Farokhalee et. al. 1989) and can be used to simulate the development of bubble populations in fluidized beds, given the initial distribution at the grid. When two bubbles interact and coalesce, the volume of the resulting composite bubble formed is generally larger than the sum of the volumes of the original pair. This growth is associated with the fact that bubbles in fluidized beds have permeable boundaries. It also implies a net transfer of gas from the dense phase to the bubble phase. Recent work by Yates and Cheesman (1992) has provided a possible new insight into why this growth of coalescence occurs. Using X-ray observations to measure voidage variations around bubbles in a Geldart group A material in a three-dimensional bed, they showed that the voidage is somewhat higher in the wake and much higher in a "shell" surrounding the bubble, reaching 0.7 at the bubble surface. The gas volume in the "shell" can be larger than in the "visible" bubble. They showed that, i f the thickness and voidage profile in the "shell" are independent of bubble size, and if the total gas volume associated 5 with the bubbles is conserved on coalescence, the new composite bubble grows by an amount which is fully consistent with experimental observations. 1.2.2 Distribution of Gas between Phases The two-phase theory proposed by Toomey and Johnstone (1952), postulating that the gas in excess of that required for incipient fluidization enters the bubble phase, may be written as U ~ Umf =GBI A (2) where U= superficial gas velocity Umf= superficial gas velocity at minimum fluidization G B = "visible" gas flow, i.e. f low due to translation of voids A= bed cross-sectional area However, it is well documented that the visible bubble flow rate is, in practice, less than (U-Umf) (e.g., see Nguyen and Leung, 1972; Grace and Clift, 1974; Yates and Rowe, 1983; Valenzuela and Glicksman, 1985; Olowson and Almstedt, 1992). Three possible reasons for the deficiency in bubble flow rate have been identified ( Grace and Clift, 1974): 1. interstitial gas flow in the particulate phase greater than its value at minimum fluidization; 2. gas "throughflow" within the bubbles probably enhanced by bubble interactions 3. regions of high interstitial voidage convected with the bubbles. 6 1.2.3. Maximum Stable Bubble Size Rowe (1971) observed that the roof of a rising bubble sometimes develops a downward cusp, which then tends to grow rapidly, causing the bubble to split vertically. When this knifing or curtain action slices off a small daughter bubble, the latter is usually reabsorbed almost immediately by the larger bubble. When the two bubbles formed are of nearly equal size, the larger one often first grows at the expense of the smaller one, which is then absorbed by the larger bubble. In some cases, the faster-rising larger bubble is able to pull away from the smaller bubble, leaving two discrete bubbles. In beds of fine particles, recoalescence is less frequent than in large particle beds. Toei et al. (1967) measured the frequency of bubble splitting in two-dimensional beds of uniform particles ranging in size from 210 to 360 pm (Umf=3-10 cm/s) and concluded that the frequency of splitting is inversely proportional to Umf, ranging from 3 to 20 s"1, and is almost independent of bubble size. In this study, the bubble frequency was 3 to 20 s"1 and the bubble diameter ranged from 3 to 13 cm. One consequence of coalescence and splitting is that a maximum stable bubble size may be present in the bed. Since bubble splitting is more frequent in Geldart type A solids than in beds of larger particles, the maximum bubble size is larger in beds of coarse particles than in beds of fine particles. 7 1.3. Conclusion Key properties of bubbling fluidized beds are summarized based on the literature. In Chapter 4 these findings are used in developing a "discrete-bubble Langrangian model" to simulate the dynamic character of bubbling fluidized bed hydrodynamics with reasonable computer machine time. 8 Chapter 2 . Introduction to Deterministic Chaos 2 . 1 Linear and Non-Linear Systems In geometry, linearity refers to Euclidean objects: lines, planes, (flat) three-dimensional space, etc. These objects appear the same no matter how we examine them. A non-linear object, a sphere for example, looks different on different scales—when looked at closely enough it looks like a plane, while from a far enough distance it looks like a point. In algebra, linearity is defined in terms of functions that have the property f(x+y) = f(x)+f(y) and f(ax) = af(x), where a is constant. Non-linear is defined as the negation of linear. This means that the change in f(x+y) may be out of proportion to the change in input x or y. The result may be more than linear, as when a diode begins to pass current, or less than linear, as when finite resources limit Malthusian population growth. In such cases the fundamental simplifying tools of linear analysis are no longer available. Linearity is rather special, and few, i f any, real systems are truly linear. Some problems are profitably studied as linear approximations to the real situation. For example, Hooke's law, the linear law of elasticity (strain proportional to stress) is approximately valid for small amplitude motion of a pendulum implying that its period is independent of amplitude. However, as the amplitude increases the period also increases, a fundamental effect of non-linearity in the pendulum equations. 9 Non-linear systems exhibit surprising and complex effects that would never be anticipated by a scientist trained only in linear techniques. Prominent examples include turbulence, bifurcation and chaos. Non-linearity has its most profound effects on dynamic systems. Further, while scientists can enumerate linear objects, non-linear ones are non-denumerable, and, as yet, mostly unclassified. Scientists currently have no general techniques (and very few special ones) for telling whether a particular non-linear system will exhibit the complexity of chaos, or the relative simplicity of order. Since they cannot yet subdivide non-linear science into proper subfields, non-linear science exists as a whole. Non-linear science has applications to a wide variety of fields, from mathematics, physics, biology, and chemistry, to engineering, economics, and medicine. This is one of its most exciting aspects—that it brings researchers from many disciplines together with a common language. A non-linear system is one whose time evolution equations are non-linear: that is, the dynamic variables describing the properties of the system (for example, position, velocity, acceleration, pressure, etc.) appear in the equations in a non-linear form, e.g. as products or raised to powers other than one. 2.2. Non-Linearity and Chaos Sudden and dramatic changes in non- linear systems may give rise to the complex behaviour called chaos. The noun chaos is used to describe the time behaviour of a system 10 when that behaviour is aperiodic and is apparently random or "noisy". Many textbooks (e.g. Hilborn, 1990) have been written in recent years providing general descriptions of chaotic systems and their characterisation. In general, three ingredients are needed to determine the behaviour of a system: 1. time-evolution equations; 2. values of the parameters describing the system; 3. initial conditions. A system is said to be deterministic i f knowledge of the time-evolution equations, the parameters that describe the system and the initial conditions are sufficient in principle to determine completely the subsequent behaviour of the system. The obvious problem is how to reconcile this underlying determinism with the apparent randomness and unpredictability of chaotic systems. Traditional scientists and engineers see most real systems like fluidised beds as having complex and random like behaviour. They try to explain some types of behaviour either on the basis of "noise" or based on "complexity". The noise argument implies that complex behaviour is due to the influence of uncontrolled outside effects. According to the complexity argument, most real systems are made of billions of small units, like particles, atoms and molecules, and it is not surprising that this leads to fluctuations and randomness in the overall behaviour of the system. Sometimes scientists explain that these complex systems have many degrees offreedom. The meaning of this notion is explained below. 11 2 . 2 . 1 A Model of Convecting Fluids: The Lorenz Model In the early sixties, a meteorologist named Edward Lorenz experimented with computer simulations of weather on a relatively primitive computer, His program used twelve recursive equations to simulate rudimentary aspects of weather; he entered several variables into his program each time he ran it, and watched to see what types of weather patterns such initial conditions would generate. Lorenz tried to recreate an interesting weather pattern, one he had seen previously, by re-entering the values the computer had previously calculated and reported. However, when he ran the program again, his results were different from the initial run. After checking the two plots, however, he realized his "error"; on his previous computer printout, the one he had used to enter the initial conditions into the computer for the second trial run, the figures were printed with three significant digits. In the program, all values were calculated to six significant digits. Lorenz had assumed that the difference, only one part in a thousand, would be inconsequential; however, due to the recursive nature of the equations, little errors would first cause tiny variations, which would then affect the next calculation a bit more, etc., which would affect the output of the next run even more. The final result of a long string of recursive calculations would lead to a weather pattern totally different from the expected values. The term "sensitive dependence on initial conditions" was coined to describe the phenomenon that small changes in a recursive system can drastically change the results of running that system. A term Lorenz coined to describe sensitive dependence on 12 initial conditions is the "butterfly effect." This is another thought experiment which is hardly testable: imagine that there exist two earths, so that an incorporeal observer could compare events on one earth to another. N o w imagine that both earths are identical except for one fact; in one, a butterfly flaps its wings somewhere in South America, and in the other, this butterfly remains still. One might think that such a small discrepancy between the two earths would be inconsequential; after all, nobody was there, nobody could even notice the butterfly's wings flapping, and air currents would be affected only in a minor way by such a miniscule event. After a period of time which is impossible to calculate, however, the weather patterns of the two earths would be totally different. Why? Because of the difference caused by the flap of one butterfly's wings! The miniscule event affected air currents around that butterfly in a very miniscule fashion, true; but those tiny air currents affected in turn slightly larger air currents, which affected still larger air currents, and the small difference in air flow between the two earths exponentially increases to become a large difference. The wind patterns on the two earths, which started out otherwise identical and had every reason to remain identical in a nice deterministic manner, would now be different in every way. Eventually, as time multiplies the differences between the two earths, completely different weather patterns would emerge, all because of the fact that on one of the earths a butterfly in South America decided to flap its wings, while on another of the earths, the same butterfly did not. The differences might not cause any major catastrophic events immediately; the thought experiment does not suggest that murdering a butterfly could cause a 13 hurricane in a few minutes or a tornado in a few hours. However, the air currents and wind patterns would be different. It is thus completely impossible, even in theory, to perform long-term weather prediction in any accurate manner. Unless a computer could be constructed which could monitor each individual atom on earth, even the smallest undetected anomaly could affect the weather in profound ways. Edward Lorenz's weather model involved 12 differential equations. He decided to look for complex behavior in a simpler set of equations. To do so he looked toward the phenomenon of rolling fluid convection. The physical model is simple: place a gas in a solid rectangular box with a heat source on the bottom. He simplified governing equations from fluid dynamics and heat transfer and ended up with a set of three nonlinear equations: dX/dt = P * (Y - X ) (1) dY/dt = R * ( X - Y - X Z ) (2) dZ/dt = X Y - B Y (3) Where X , Y and Z are the coordinate function of an element of fluid. Here P is the Prandtl number representing the ratio of the fluid viscosity to its thermal conductivity, R represents the difference in temperature between the top and bottom of the system, and B is the ratio of the width to height o f the box used to hold the system. The values Lorenz used are P = 10, R = 28, B = 8/3. On the surface these three equations seem simple to solve. However, they represent an extremely complicated dynamical system. If one plots the results in three dimensions, Figure 2 .1, called the Lorenz attractor, is obtained. (X_plot, Y_plot, Z_plot) Figure 2 . 1 . Lorenz attractor 15 2.2.2.Logistic M A P (Devaney, 1989). The so-called Logistic equation is used in biology to model animal populations. The equation is Xn+l=cXn(l-Xn) (4) where X is the population (between 0 and 1) and c is a growth constant.The subscript denotes generation. Iteration of this equation yields the period doubling route to chaos. Note the period means the years the population repeats. For c between 1 and 3, the population settles to a fixed value. A t c=3, the period doubles to 2; one year the population is very high, causing a decline to a low population the next year, leading to a high population the following year. At c=3.45, the period doubles again to 4, meaning the population has a four year cycle. The period keeps doubling, faster and faster, at 3.54, 3.564, 3.569, and so forth. At c=3.57, chaos occurs; the population never settles to a fixed period. For most c values between 3.57 and 4, the population is chaotic, but there are also periodic regions. For any fixed period, there is some c value that wil l yield that period. The bifurcation diagram of the logistic difference equation is shown below in Figure 2.2. 2.3. Reasons for Quant i fy ing Chaot ic Behaviour Reasons for quantifying chaotic behaviour are as follows (Hilborn, 1990): 1. The quantifiers may help distinguish chaotic behaviour from noisy behaviour; 16 2. The quantifiers may help determine how many variables are needed to model the dynamics of the systems; 3. The quantifiers may help us sort systems into universality classes; 4. Changes in the quantifiers may be linked to important changes in system behaviour, such as changes in flow regime. Steady state values of pn Value of r BIFURCATION DIAGRAM Pn+l=rpn(1-pn) S ta r t i ng va lue p = .02 p„ converges to a limit (r = 2.50) Figure 2.2. Bifurcation diagram of logistic difference equation 17 There are two different, but related, types of quantification or description. One is a dynamic (time-dependent) description. The other is a geometric description. Both are discussed and used in this thesis to explore bubbling fluidisation systems. 2.4. Time Series Analysis Time series analysis is the application of dynamical systems techniques to a data series, usually obtained by "measuring" the value of a single observable variable as a function of time. The major tool is "delay coordinate embedding" which creates a phase space portrait from a single data series. It seems remarkable at first, but one can reconstruct a picture equivalent (topologically) to the full Lorenz attractor in three-dimensional space by measuring only one of its coordinates, say x(t), and plotting the delay coordinates (x(t), x(t+x), x(t+2x)) for a fixed value of x. The idea of using derivatives or delay coordinates in time series modelling is not new. It goes back at least to the work of Yule who used an autoregressive (AR) model to make a predictive model for the sunspot cycle. A R models are nothing more than delay coordinates used with a linear model. Delays, derivatives, principal components, and a variety of other methods of reconstruction have been widely used in time series analysis since the early 1950s, and are described in several hundred books. The new aspects raised by dynamic systems theory are : (i) a geometric view of temporal behaviour, and (ii) the existence of "geometric invariants", such as dimensions and Lyapunov exponents. The central question is not whether delay coordinates are useful for time series analysis, but rather whether reconstruction methods preserve the geometry and the geometric invariants of dynamical systems. Some common methods and numerical algorithms of time series analysis are as follows: 2.4.1. T ime History Time history is a common way to classify a system: Unlike a periodic or quasiperiodic signal that features regularity, a chaotic signal exhibits irregularity with respect to time. A chaotic signal with different time resolutions is such that successive magnification will uncover more and more details about the fluctuating signals. This means that chaotic signals for different time scales exhibit self-similarity. However, since the irregularity can also be due to randomness, distinction between chaotic and random signals can be realised by fractal dimension analysis, as discussed below. 2.4.2. Strange Att ractor Phase space is the collection of possible states of a dynamic system. A phase space can be finite (e.g. for the ideal coin toss we have two states -heads and tails), countably infinite (e.g. state variables as integers), or uncountably infinite (e.g. state variables as real numbers). Implicit in this notion is that a particular state in phase space specifies the 19 system completely; it is all we need to know about the system to have complete knowledge of the immediate future. Thus the phase space of the planar pendulum is two-dimensional, consisting of the position (angle) and angular velocity. According to Newton, specification of these two variables uniquely determines the subsequent motion of the pendulum. Note that i f we have a non-autonomous system, where the map or vector field depends explicitly on time (e.g. a model for plant growth depending on solar flux), then according to our definition of phase space, we must include time as a phase space coordinate—since one must specify a specific time (e.g. 3 P M on Tuesday) to know the subsequent motion. Thus dz/dt = F(z,t) is a dynamical system in the phase space consisting of z and t. The path in phase space traced out by a solution of an initial value problem is called an orbit or trajectory of the dynamical system. If the state variables take real values in a continuum, the orbit of a continuous-time system is a curve, while the orbit of a discrete-time system is a sequence of points. A strange attractor is the limit set o f a chaotic trajectory. A strange attractor is topologically distinct from a periodic orbit or a limit cycle. A strange attractor can be considered a fractal attractor. A n example of a strange attractor is the Henon attractor. Consider a volume in phase space defined by all the initial conditions a system may have. For a dissipative system, this volume shrinks as the system evolves in time (Liouville's Theorem). If the system is sensitive to initial conditions, the trajectories of the points defining initial conditions move apart in some directions, closer together in others, 20 but there wil l be a net shrinkage in volume. Ultimately, all points lie along a fine line of zero volume. This is the strange attractor. A l l initial points in phase space which ultimately land on the attractor form a Basin o f Attraction. A strange attractor results i f a system is sensitive to initial conditions and is not conservative. In practice experiments are commonly done with single-variable time series measurements made at evenly spaced time intervals. Reconstruction of an attractor from such a single-variable time series using the time delay method provides a way to investigate the real multi-dimensional system in nature. 2.4.3. Autocorrelation Function The autocorrelation function, sometimes called the correlation function, can be defined as 1 N-T r W = ^ (5) The autocorrelation function is obtained by multiplying each x(t) by x (t- x) and summing the result over all data points. The sum is then plotted as a function of time delay x. This gives a measure of how strongly data points are correlated with neighbouring points. The value of x at which the autocorrelation function remains essentially zero is defined as the correlation time. Highly random data have no correlation, and for these data the correlation function drops abruptly to zero, implying a small correlation time. Highly 21 correlated data on the other hand have an amplitude that decreases only slowly. Chaotic data tend to show little correlation except when the time delay is very small. 2.4.4. Fracta l Dimension The fractal dimension describes the complexity of a chaotic attractor, and, in principle, provides a measure of the freedom of the system or the least number of independent variables needed to describe the dynamics on the attractor (Barnsley, 1988; Parker and Chua, 1989; Moon, 1992; Mullin,1993). and Box Count ing Method Among many definitions of fractal dimension, there are three in common usage -capacity, information and correlation dimensions. These fractal dimensions can be calculated by the Box-Counting method. In a d-dimensional state space where the attractor lies, the box counting method works to cover the attractor with d-dimensional hypercubes of linear size r. The minimum number of hypercubes, N(r), needed to cover the attractor is then counted. Naturally, N(r) increases as r is decreased. From the limit as r tends to zero, we have the following definitions: c a p ~<> l n ( l / r ) v ' 22 2>-ln# Di = l i m - ^ - (7) r-*° Inr w Af(r) 2 D c = l im—^ (8) Here pi is the relative frequency (or probability) with which a typical trajectory enters the ith volume element, i.e. p;=N;/N. Clearly, the capacity dimension D c a p , is a purely geometric measure of the attractor, while the other two dimensions, D i and D c , account for the non-uniform distribution of the attractor in the phase space. Differences occur between these three dimensions, but we always have D c a p > D i > D c . Two methods have been developed (Grassberger and Procaccia, 1983a; Badii and Politi, 1985) to calculate fractal dimensions in the case where a 2-dimensional phase space is considered: Method of Grassberger and Procaccia(1983a) This method is the most widely used algorithm for calculating the correlation dimension. In this approach, a vector with M elements is first determined as follows: x(t,), x(t1+r), x(t i + 2 r), X 0-i+(d-l ) r ) x = x(t,) = • x(tj), x(tj+r), x(tj+2r), . . . X(tj+(d-l)r) (9) x(tM)_ x(tM), x(tM + r), x(tM + 2 r), X 0 - M + ( d - ! ) r ) 23 The M points can be determined either randomly or successively beginning with a starting point. Then the distances between pairs of points, say 5ij= calculated using the maximum norm x-x, are <?,=||X(0-X(r,)|| = maxfl*(0-x(fy)|, \x(ti+r)-x(tJ+t)\, |x(fl+(,_1)T) - x ' t M , _ ^ } (10) The correlation dimension is D c then defined as C(r) oc rD< 01) where the correlation function, C(r), is given by i M C(r) = l im [YHJr - X, - X. (12) and Hf is the Heaviside function defined by Hf(x) = 1 for x > 0 0 for x < 0 (13) The correlation integral C(r) may be considered as the probability that two points on the attractor lie within a cell of size r. Since 8ij=6ji, calculations for only M ( M - l ) / 2 pairs (i=l, M ; j= l , M-i) are needed. Hence 24 2 M M-i C(r) = l im -*«> M(M - 1 ) (14) i=l ;=1 Note that equation (9) is only valid for r reasonably small compared to the dimensions of the attractor, but large enough to reduce the contribution of instrumental uncertainty. Grassberger et al. (1991) provide a list of practical issues that have to be dealt with when devising an algorithm for the calculation of the correlation integral. Summation over all pairs (i, j) increases the statistical accuracy of C(r). Calculation (Bai et. al., 1997) has shown that a very reasonable average can be obtained if M is taken as 100 points or so, resulting in considerable economy of computation time. In our approach, all pairs i, j were considered and M=1000 was usually used. To increase the statistical accuracy, an average C(r) corresponding to a specific distance r is usually based on 5x10 6 to 5x10 7 calculated distances. Repeated selection of the vector on the attractor gave a standard deviation of the values of less than 1%. Finally, the correlation dimension can be determined from the slope of the linear portion of the log-log plot of C(r) versus r with a correlation coefficient higher than 0.98. Method of Badii and Politi (19851 A n alternative method, called the nearest-neighbour method was proposed by Badii and Politi (1985, 1987). In this approach, a central point X . is selected on the 25 reconstructed attractor. A subset of k points, denoted X}. (j=l, 2 , k and k<N), is also chosen at random from the original set of N points. W e then define 5i as the single minimum of the distances from Xj to X}, i.e. the distance to the nearest neighbour: 5, = min Xt - Xj (15) To obtain a statistically useful value of the minimum distance, this calculation is repeated for 2000 sets of randomly selected central points from which an average value of 5 is obtained. The dimension of the attractor is then obtained from For the same data, the correlation dimension obtained by the nearest neighbour method should be consistent with that obtained by the correlation integral method. The method of Badii and Politi (1985) gives more extensive linear portions than that of Grassberger and Proccaccia (1983 a), thus improving the reliability of obtaining slopes and hence D c . Moreover, the Badii and Politi (1985) method requires less computation time than that of Grassberger and Proccaccia (1983 a). Suppose we have measured a time series of N points (xi, x 2 , . . . ,x N), where N is typically of order 104 to 105. From this time series the attractor is reconstructed using delay time coordinates based on Takens' theorem (1981). To determine the Lyapunov log£ = -—-logk (16) c 26 exponents, one can conceptually imagine a reference trajectory in phase space. As the initial condition, for an arbitrary point on this trajectory, X i = X i ( ' o ) = { x r > x ; + r . x I + 2 r . •••> x I + ( d - i ) r } 0?) while for a point on a trajectory (preferably the nearest neighbour) Xj = Xj(r0) = {x j ; x j + r , x j + 2 r , x j + ( d . ] ) r j (18) The two points are initially close to one another in phase space. When the time advances t e t=kAt to ti=to+tet (where At is the sampling time interval and k is an integer number), X;(t\) = Xj(/0 + tet) = Xj+k = : | x i + k , x i + k + r , x i + k + 2 r , x i + k + ( d _ n r | (19) and Xj(^|) = Xj(?0 + tet) — Xj+k =|xj+k, X j + k + r , X j + k + 2 r , X j + k + ( d . 1 ) r | (20) The distance between these two points, 5jj(t,), is calculated based on equation (16). <w=|*,('i )-*,(',)| 1/2 (21) For a chaotic attractor, the distance grows exponentially with time, i.e., W = <W2A(',"') (22) 27 Hence X= l i m — ^ l o g 2 ^ ^ - ( 2 3 ) Base 2 is used instead of e so that the Lyapunov exponent has units of bits per unit time (i.e. bits/s). It is clear that one measurement is insufficient since the divergence of chaotic orbits can only be locally exponential. Instead the calculations must be averaged over different regions of phase space. This average can be represented by — m t ^ - t ^ & ( 2 4 ) If the time evolves along the trajectory sequentially, we have m(t k-t k.i)=tm-t0 and equation (24) becomes m lm l0 k=l °ij(lk-lJ In practice, we have a finite time series. The use of a finite number of experimental data points does not allow us to probe the desired infinitesimal length scales of an attractor. However, for the practical case of limited data one can estimate the Lyapunov exponent A,et for a specific selected evolution time tet. Practical methods for estimating Lyapunov exponents from a time series are available in the literature (Sano and Sawada, 1985; Wol f et al., 1985; Sato et al., 1987; 2 8 Zeng et al., 1991, 1992; Parlitz, 1992;). Wo l f et al. (1985) also provided a computer programs for calculating the Lyapunov exponent corresponding to a fixed time evolution tet=kAt. The calculation is initiated by carrying out an exhaustive search of the data file to locate the nearest neighbor to the first point (the fiducial point). The distance, 8y(to), between the current pair of points is next calculated. The current pair of points is then propagated k steps (each of step length te t) through the attractor, and its final separation, 5ij (ti), is computed. The logarithm of the ratio o f the final to the initial separation o f this pair updates a running average rate of orbit divergence. When 5 y (ti) becomes too large, i.e. the growth departs from exponential behaviour, a replacement step is attempted. In this step, the program looks for a new nearest data point, 5 y(ti), as the "new" initial condition. If more than one candidate point is found, the point defining the smallest angular change, G(ti), is used for replacement. If an adequate replacement point cannot be found, the points that were being used are retained. This procedure is repeated until the reference trajectory has traversed the entire data file, at which point we estimate 1 8 (I ) K = — L - Z l o § 2 (26) where m is the total number of replacement steps and the time step, te t=tn-tn-i, is held constant. In practice, maximizing the evolution time is highly desirable as it both-reduces the frequency with which orientation errors are made and greatly diminishes the cost of calculations. To find an adequate evolution time, a series of calculations is recommended 29 for various evolution times. A constant Xv i.e. the largest Lyapunov exponent of the attractor, is approached as tet is increased beyond 500 to 700. Hay et al. (1995) suggested the following empirical function to determine the Lyapunov exponent: K=J- + *>x (27) However, this method depends greatly on the distribution of data points calculated, as they give different weights to the parameter fitting. Therefore, it is generally reasonable to find Xi by averaging values of Xei obtained at large tet. 2.5. Previous Appl icat ion of Deterministic Chaos Theory to Flu id ized Beds The theoretical foundations of chaos are on a firm footing. Chaotic behavior has been studied by analytical, computational, and experimental means. Manifestations of chaos involve a wide range of systems: mechanical, electrical, and optical devices, hydrodynamic processes of various length scales, transport processes, biological and chemical reactions, etc. However, the findings have often been a posteriori, i.e., explaining observed complex behavior and demonstrating that the roots of the complexity can be traced back to a deterministic explanation. 30 Despite these advances, applications of the concepts of chaos to engineering and design of new processes and products are very few. Nevertheless, the theoretical concepts can lead to new applications of interest in chemical engineering. For example, Ottino (1992) improved the basic understanding of stirring and mixing by means of chaos theory. Some of the basic developments are ready for implementation in practical situations. The most obvious possibilities are the use of basic concepts for the design of new mixing devices. Chang and Sen (1992) applied chaos to the improvement of mixing and heat transfer. They used a mathematical model to find the optimal forcing frequency of mixing equipment for maximization of chaotic enhancement. The application of deterministic chaos theory to fluidized beds has been pioneered by several research groups (e.g. Daw et al., 1990; Daw and Halow, 1991, 1992, 1993; Skrzycke et al., 1993; Schouten and van den Bleek et al., 1991, 1992; Schouten et al., 1992; van der Stappen et al., 1993a, 1993b, 1994; van den Bleek and Schouten, 1993a, 1993b; Tarn and Devine, 1991; Hay et al., 1995). Although Tam and Devine (1990) were unable to identify a saturation dimension with increasing embedding dimension up to 12, other researchers have established that gas-solids fluidized beds are deterministic chaotic systems, at least for some operating conditions. This opens up the possibility of extracting important information, previously unknown about fluidized beds, by employing chaotic time series analysis. Such information might also be useful for monitoring, control, design and scale-up of fluidized bed reactors. In this section, possible applications of deterministic chaos to fluidization are discussed briefly. 31 2.5.1. Fluidization regime identification In addition to conventional measures used to distinguish between flow regimes such as the standard deviation and probability distribution of absolute or differential pressure fluctuations, flow characteristics of different fluidization regimes can potentially be recognized by chaos analysis. Experimental and modeling work have been performed to classify the fluidization regimes using deterministic chaos analysis (e.g. Schouten et al., 1992; Skrzycke et al.,1993). For steady state operation of a packed bed, it is expected that the correlation dimension based on voidage fluctuations should be zero because of the absence of solids motion and hence of voidage variations. Above the point of minimum fluidization, the correlation dimension increase with increasing superficial gas velocity, suggesting an increased degree of complexity and possibly chaos. Daw and Halow (1991) analyzed the bubbling and slugging regimes based on bubble-bubble interactions. The correlation dimension based on differential pressure fluctuation signals ranged from 2.5 to 2.8. Skrzycke et al. (1993) found that the local level of chaos in a fluidized bed reaches a minimum just beyond a regime transition. For example, minima in Lyapunov exponent (Bai et al., 1996) were found at the transition from bubbling to turbulent fluidization. The correlation dimension based on pressure fluctuation measurements has been found to range from 2 to 5 in bubbling and slugging fluidized beds (Schouten and van den Bleek et al., 1991, 1992; Schouten et al., 1992; van den Bleek and Schouten, 1993; Hay et al., 1995). In a two-dimensional phase space, Bai et al. (1996) found that the capacity dimension exhibits a minimum at the onset of the turbulent fluidization regime. 32 2.5.2. Characterization of hydrodynamic structure The hydrodynamics of fluidized beds have been extensively investigated using conventional statistical tools such as probability distribution, standard deviation and F F T spectrum. These tools, however, provide no time-dependent information. In contrast, chaos analysis focuses on the dynamic behavior of fluidized beds, thereby providing a powerful tool to explore the underlying structure . Previous studies have demonstrated that the dynamic behavior of fluidized beds, as expressed by the correlation dimension, is dependent on operating variables such as the gas velocity, particle properties, column diameter and bed depth (van der Stappen et al., 1993; Skrzycke et al., 1993). Phase-space trajectories have characterized aspects of bubble behaviour such as bubble formation, breakup and coalescence in bubbling and slugging beds (Daw and Halow, 1993). The bubble properties are also reflected by the correlation dimension (Daw and Halow, 1991; van den Stappen et al., 1993, 1995; Hay et al., 1995 ). Chaotic analysis may also permit detection of underlying flow structures from temporal data, this being impossible using conventional statistical means. Superposition of chaotic signals may occur, reflecting separate processes co-existing in the system (Franca et al., 1991; Izrar and Lusseyran, 1993; Bai et al., 1996). In general, the gas-solids flow in low- and high-velocity fluidized beds is extremely heterogeneous. Low-velocity fluidized beds (bubbling and slugging flow regimes) consist of two distinct phases, a dense emulsion phase and a dilute void phase (Grace, 1984), while particles in high velocity (fast fluidized) 33 beds can be considered to exist in two distinguishable forms, as dispersed particles and as clusters (Grace, 1990). The flow patterns o f the characteristic two phases should be distinguishable from experimental pressure fluctuation signals. Properties o f the two phases, including their relative fluctuating magnitudes, complexities and sensitivities to initial conditions, can then be resolved by chaotic analysis ( Ba i et al., 1996). Studies (Bai et al., 1996) have shown that the dense emulsion phase should generate small scale rapid oscillations, thereby yielding a higher correlation dimension and a quicker rate of information loss. The motion of the bubble phase results in large-scale low-frequency pressure fluctuations, producing a lower correlation. 2.5.3. Flu id izat ion Qual i ty Cont ro l Deterministic chaos analysis provides information on dynamic fluidized beds including the number of degrees of freedom and sensitivity to changes of initial conditions. This assists in monitoring the "quality" of fluidization by specifying the number of variables required to characterize the system, and allowing the operating conditions to be adjusted to improve operations. There are two possible strategies for controlling the chaotic behaviour of fluidized beds: global control and high-speed control. The basic approach in global control is to monitor time-averaged properties over tens of seconds or minutes of the phase-space trajectories reconstructed from one variable (e.g. pressure). Once desired limits of the time-averaged chaotic properties are known (e.g., dimension, principal component eigenvalues, or general appearance of the visualized trajectory), action can be taken to maintain the fluidization system within the desired range by 34 adjusting critical operating variables that define the dynamic state o f the fluidized bed. Such a global control system has been demonstrated in laboratory simulations of a batch fluidized bed reactor processing Geldart group D (i.e. relatively coarse) uranium particles (Daw and Halow, 1993), who found that the second principal component o f the pressure drop signal provided a quantitative indicator of fluidization quality that very quickly (e.g. in less than 5-10 seconds) detects changes in fluidization behaviour due to gas flow disturbances, particle attrition or changes in gas and particle properties. Using the resulting deviation of the eigenvalue from a set point as an error signal, it was shown that a PI controller could provide a continuously updated correction signal to the control valve. High-speed control focuses on short-time-scale behaviour. By constantly and promptly updating the information on the current state of the system, it should be possible to make relatively accurate short-term predictions of the fluidization state. The key is the ability to make extremely rapid (e.g. milliseconds) updates on the current dynamic state, predictions about its future course, and decisions regarding the appropriate control perturbations. Such capabilities should soon be within the range of computers and analytical algorithms. Tests should be carried out to test such control strategies and to compare them with more conventional techniques. 3 5 2.5.4. F lu id ized Bed Scale-Up The scaling laws of fluidized beds have been extensively investigated (e.g. Glicksman, 1984, 1988; Glicksman et al., 1994). These scaling laws have been developed based on fundamental or empirical hydrodynamic models. (For a detailed description, see Gogolek and Grace, 1995). Schouten and van den Bleek (1991) noted that these scaling laws do not include time-dependent information on fluidized bed hydrodynamics, van den Bleek and Schouten (1993) suggested that i f two beds are properly scaled, the rate of information change in both systems must be the same. A major advantage of this approach is that the same particles can be used in small and large columns. Thus difficulties of finding proper particles or fluids in previous scaling laws are avoided. Although a large volume of research has been done since the end of the last decade, a very interesting topic opened by Bai et al.(1996) has not yet been settled. From the discovery of double-sloped correlation integral curves of pressure signals obtained from fluidized beds first noted by Tam and Devine (1992), Bai et al. suggested that fluidized beds are high dimensional multi-fractal dynamic systems. Because the gas-solids flow in gas fluidized beds is heterogeneous, bubbling fluidized beds usually consist of two distinct phases, a dense emulsion phase and a dilute bubble phase. These two phases make separate and distinguishable contributions to the pressure fluctuation signals. The dense emulsion phase should generate small scale rapid oscillations, thereby yielding a higher correlation dimension and a quicker rate of loss of information. On the other hand, the motion of the bubble phase results in large scale low-frequency pressure fluctuations, 36 producing a lower correlation dimension. More direct proof is needed to support this suggestion. Further research would be simplified i f pressure fluctuations could be separated into two parts, one corresponding to the dense phase and the other to the dilute phase. Time series analysis of those separated signals could help to validate this suggestion. However, because of the complexity of fluidized systems, pressure fluctuation separation is very difficult. Some important features of pressure and voidage measurements are summarised in Table 2.1. 2.6 Scope of Remainder of Thesis In the present research project, experimental and predicted void data generated from fluidised beds are used to explore the structure of chaotic hydrodynamics of fluidised beds. In the work described in this dissertation, experimental voidage data were collected using optical fibre probes and investigated using deterministic chaos theory. Local voidage signals can be divided into two parts, one corresponding to the bubble phase, and the other to the dense phase. This bipartite nature of voidage signals is an advantage over pressure signals for research on fluidised bed chaotic hydrodynamics. The experimental work is described in Chapter 3. In addition to the experimental work, simulation of three-dimensional fluidised beds using the coalescence model of Clift and Grace (1970, 1971, 1972) has been used to 37 predict time series voidage data. The idea is to show whether chaotic behaviour can be predicted based on the behaviour of bubbles in fluidized beds. Details of these simulations are given in Chapter 4. Conclusions and recommendations for future work are summarised in Chapter 5. 38 Table 2.1 Relative advantages and disadvantages of voidage and pressure measurements P R E S S U R E • sensors available for a wide range of operating conditions (including high temperature and pressure) • sensors commercially available • produce global information about bubbles in fluidized bed • high response frequency • smooth, continuous time signals • physical interpretation difficult • choice between absolute pressure and differential pressure; latter gives information which is more localized V O I D A G E • optical probes have been developed for cold models only • probes are commercially available • local information only • very rapid response rate • signals are 'binary' in dense beds or spiky for small measuring volumes • information depends on measurement volume, which is usually very small 39 Chapter 3. Deterministic Chaos Analysis of Experimental Voidage Signals in a Freely Bubbling Fluidised Bed 3.1 Introduction A number of experimental parameters have been utilised to show chaotic phenomena in fluidised beds including pressure and voidage fluctuations. A comparison of the pros and cons of voidage signals and pressure signals was provided in Chapter 2. Fluidised beds are known to behave as chaotic systems for most operating conditions. Different signals show different, but related, chaotic phenomena. Most previous research on deterministic chaos in fluidised beds has been based on pressure signals. The work presented in this chapter focuses on voidage signals obtained from an optical fibre probe in a freely bubbling fluidised bed. According to previous work in our laboratory (Zhou, 1994) the particles are large enough to give a linear relationship between the optical fibre probe output and local voidage. The overall hydrodynamics of freely bubbling fluidised beds of Geldart Class B particles can be predicted with reasonable accuracy. Voidage signals from optical fibre probes are very localised signals and can be used to show localised chaotic behaviour at different positions in freely bubbling fluidised beds. In this chapter, localised chaotic behaviour is reported and related to the overall chaotic hydrodynamics of fluidised beds. Optical fibre signals have the advantage that they can be divided into two parts corresponding to the emulsion phase and the dense phase. Time series analysis of both these reconstructed signals helps to understand the chaotic structure o f bubbling fluidised beds. 3.2. Experimental Equipment and Probe Calibration A Plexiglas column was used to collect voidage signals. The column has an ID of 150 mm and an overall height of 1.0 m. A P M M A (polymethylmethyl acrylate) plate covered with fine screen was employed as the gas distributor with 109 holes of 3.2 mm in diameter uniformly distributed over the plate, giving an open area ratio of 3%. A n optical fibre probe was inserted at various heights along the column with its tip located at the axis of the column. Glass beads (provided by Manus Industries) of sizes between 60 to 80 U.S. mesh (i.e. 200 to 270 urn) and particle density 2470 kg/m 3 were used as the bed material. Broken or angular particles comprise not more than 3% by count. According to the powder classification of Geldart (1973), these materials are typical Class B particles. The fluidizing gas is compressed air at atmospheric temperature and pressure. A l l measurements were carried out at a superficial gas velocity of 0.2 m/s, a static bed height of 230 mm and an expanded bed depth of 310 mm. The minimum fluidisation velocity is 0.1 m/s, while the voidage at minimum fluidisation is 0.41. Reflection type optical fiber probes have been used previously in our laboratory (e.g. He., 1994; Zhou, 1994; Pianrosa, 1996; Issangya et. al., 1996). The principle o f operation was described by Qin and L iu (1982). The probe is composed of two bundles o f 0.015 mm diameter fibres arranged layer by layer, with alternate layers used for light projection and light reception. Visible light projected from the light source to the fluidised bed through half the fibres is reflected by particles and received by the receiving fibres. 41 The light signals are then converted to electrical signals by a photomultiplier. A n amplifier raises the resulting signal to the voltage range of 0-5 V after which it is recorded by a personal computer via a D A S - 8 AID converter using Labtech Notebook software. The accuracy of optical fibre probe data relies upon the precision of the calibration. Various methods have been utilised. Matsuno et. al. (1983) calibrated their probe in a stream of glass beads falling from a vibrating sieve and estimated the concentration of the particles passing in front o f the probe from the measured particle flux and particle velocity, assuming the latter to be equal to the single particle terminal falling velocity. This method is only valid for dilute systems. They found a linear relationship between voidage and the output of the fibre optical probe. Qin et. al.(1982) found similar results. However, Issangya and Bai (1996) found that there is a non-linear relationship for particles of the order of 100 urn in diameter or smaller. It has been reported that the linearity of the relationship is dependent on the size of the particles (Lischer and Lounge,1992; Yamzaki et. al., 1992; Werther et al., 1993), with the relationship being linear for particles larger than about 200 urn. The shape, colour and surface characteristics of particles play a role in this relationship also. Round and clear particles like glass beads are more likely to have a linear relationship, while amorphous particles are likely to show nonlinear relationship between particle concentration and reflected light intensity. In view of the different calibration results discussed above and the probable dependence of the signal intensity on both the fluid and solids properties, Issangya et al. (1996) developed a method to calibrate optical fibre probes. The equipment is shown in Figure 3.1. Particles are made to fall from an incipiently fluidised hopper through a short piece of tube in the centre of a punched plate distributor covered with fine wire mesh via 42 Equipment used for calibration of fibre optic/ momentum probe 500 mm 100 mm 100 mm 2 m 5 "OD 3" ID _screw (50 mm long) Air slide valve 1 slide valve 2 Balance screw 1/4" hole total 3 unformly distributed SLIDE VALVE 2" mnJL front view d3/16" thick Notes: 1. Al l materials are plexiglass except plate of slide valve 2. distance between the two valves depends on the size of pressure tap 3. dimensions can be modified upon discussion Figure 3.1 Optical fibre probe calibration equipment proposed by Issanga et. al. (1996) a 12 mm ED tube where they fall into a collection vessel. Once a steady flow is ensured, the signal from the optical fibre probe, inserted in the tube 1.5 m from the top, is recorded at 200 H z for a duration of 30 s. Then the two slide plates are closed simultaneously to 43 trap the small volume of solids surrounding the probe. Weighing these particles allows the corresponding solids concentration to be calculated. Because the glass beads fill the tube when the solids volumetric concentration is high, the manually controlled slide valves cannot be operated smoothly, promptly and simultaneously as needed for accurate measurement. Therefore calibration can only be achieved for intermediate and low solid concentrations. A n extra point at 0.41 voidage is achieved by inserting the optical probe into a loose packed bed. The resulting calibration is shown in Figure 3.2. From Figure 3.2, a linear relationship between particle concentration and the optical fibre probe signal is observed as a reasonable approximation for medium and low 2 >-V o i d a g e F i g u r e 3. 2 C a l i b r a t i o n c u r v e f o r o p t i c a l f i b e r p r o b e solid concentrations. From the above literature review, the linear relationship is expected to extend right to the high solids concentrations corresponding to minimum fluidisation. Hence a linear relationship is applied to the complete voidage range of interest (i.e. from S m f t O l ) 4 4 Intensity (V) Intensity (V) Intensity (V) ro co 4^  o i 9V Intensity (V) Intensity (V) Intensity (V) CQ C —s CD CO 4^ c 3 cr CD O o O OJ H 3' CD c+-T3 CO O ^ cr c/> CD CD CT) 3 c5" "2. I. CQ <=t- CD JQ =T C CD 0 CQ" £ 5 o -*» o ro o 3 ^ o o o ro o o o C O o o o 4^ o o o o o o CD O O o o o o I M ' I 1 I 1 +J~f 73 II o LP Intensity (V) Intensity (V) Intensity (V) 3.3 Experimental Voidage Distribution for a Freely Bubbling Fluidised Bed When the fluidised bed was operated at a superficial air velocity of 0.2 m/s, a characteristic bubbling regime was achieved. The static bed height is 230mm. The expanding bed height is 270mm. Small bubbles were generated from the distributor and i • i > r Radial posi t ion r/R Figure a 6 Radial voidage profi I es at U=0.2m/s for di f ferent bed height coalesced with each other to produce larger bubbles. The non-uniformity of bubble distribution can be quantified by optical fibre probe with simplicity, accuracy and low cost. A l l the information registered by optical fiber probe relates to very local voidage fluctuations. Figures 3.3, 3.4 and 3.5 are plots of optical fiber probe signals at different heights and radial positions in the bubbling fluidised bed. Figure 3.6 plots the local time-mean voidage from the optical fibre probe at different heights and radial positions. It clearly shows that the bubble distribution in this free bubbling fluidised bed is non-uniform. At 48 bottom, bubbles are distributed relatively uniformly along the radius. A s they coalescence, they migrate towards the axis of the fluidised bed as in earlier studies (e.g. Grace and Harrison, 1968) 3.4 Deterministic Analysis of Optical Fibre Signals Optical fibre probe signals were sampled at a frequency of 200 H z frequency for periods of 100 seconds. This sampling frequency and quantity of data are sufficient to yield reliable results, as discussed below. Although the optical fibre probe measuring systems is subject to a pulse noise which affects the measurements, the signals were sampled directly without any filter in order to improve the accuracy. Before performing any deterministic chaos analysis, however, numerical data smoothing was applied for noise reduction. 3.4.1. Analysis conditions In order to obtain reliable deterministic analysis results from optical fibre probe signals, the selection of parameters is critical. The key parameters are the embedding dimension, time delay, number of sample points and length of analysis vector. In principle, a reconstructed attractor can be obtained for any sufficiently large embedding dimension d and almost any time delay x. However, i f d is chosen too large, noise in the data tend to decrease the density of the points defining the attractors, making it harder to find replacement points. Noise, an infinite dimensional process, unlike the deterministic component of the data, fills each available phase dimension in a reconstructed attractor. Increasing d beyond what is minimally required has the effect of 49 unnecessarily increasing the level of contamination of the data and the computation time. Experience shows that the appropriate value is about twice the dimension of the attractor. This value becomes 8-15 for pressure signals obtained from fluidised beds that have a correlation dimension ranging from 3 to 7. Bai et. al.(1997) have shown that an appropriate value of embedding dimension can be selected as the value beyond which the determined correlation dimension attains a constant value. Some care is required in choosing the time delay for accurate estimation of the chaotic invariants. In general, excessively large time delays effectively result in uncorrelated points within a vector, thereby giving properties of stochastic time series. When the time delay is chosen too small, all subsequent points within a vector are highly correlated, and the attractor is stretched out along the diagonal in the pseudo-state space. In this respect, the time delay at which the autocorrelation function first attains a minimum or crosses zero can be considered an appropriate time delay for reconstructing the attractor. In this work, the time delay is chosen as the time at which the autocorrelation function decreases by 80%. This choice gives us similar results as the first zero-crossing. Theoretically, the number of data points required in the calculations is proportional to 10 d / 2 (where d is the dimension of the attractor). Experience has shown that 5,000 to 10,000 points are sufficient for fluidized bed signals (Bai et. al.,1997). In this work, the number of data points was selected as 10,000. The length of the vector should be large for the calculation. Experience has indicated that when the length of the vector is greater than 100, reliable results are achievable. The length in this work was chosen as 400. In order to show chaotic phenomena from the localised optical fibre probe in the fluidised bed, the optical fibre probe was placed at different positions. This allows 50 comparison of the information among these signals. The locations and designations are listed in Table 3.1. Table 3.1. Optical fibre probe sample locations in the free bubbling fluidised bed and their codes (U=0.2 m/s, H=270 mm) Z=height center 1/2 radius 3/4 radius 55 mm OL1 O L 2 OL3 110 mm OM1 O M 2 OM3 220 mm OH1 O H 2 OH3 Figure 3.7 plots the autocorrelation function for case O M 1 . The time delay is chosen as the time at which the autocorrelation function decreases by 80%, which is about 20. 51 Figure 3.8 Phase space plot for case OM1 A phase space plot for case OM1 is shown in Figure 3.8. The derivatives are assumed to be given by half the difference between the data points adjacent to each point divided by the s interval The strange attractor is revealed in the phase space plot. The correlation dimension for case OM1 is shown in Figure 3.9. When the embedding dimension increases from 23 to 27, the correlation dimension is virtually the same. This means that for the time delay 10, the embedding dimension of 23 is saturated. Lyapunov exponents o f optical fibre signals collected at different locations in a fluidised bed are listed in Table 3.2. Based on this result, in the high bubble frequency region, i.e. the bottom section of the bed, the Lyapunov exponent of the voidage signal is bigger than in the low bubble frequency region, found in the top section near the wall. In other words, the region with many bubbles is more chaotic than the region with few bubbles. On the other hand, in the region where the bubbles tend to be large, i.e. in the center of the top section, the Lyapunov exponent is close to the Lyapunov exponent 52 of voidage signals collected where bubbles tend to be small, i.e. the bottom section of the bed. This indicates that bubble size has only a weak influence on the Lyapunov exponent. Table 3.2. Lypunov exponents of optical fibre signals collected at different locations in a fluidised bed for U=0.2 m/s, H=270 mm Z=height, mm center 1/2 radius 3/4 radius 55 0.75 0.70 0.66 110 0.66 0.62 0.60 220 0.64 0.55 0.40 3.5 Deterministic chaos analysis of re-constructed optical fibre signals Double slope curves for the correlation dimension have been obtained in both pressure and voidage signals. Pressure signals can be depicted as superimposed signals. Voidage signals are collected by probes that alternatively see the bubble and dense phases Although pressure and voidage signals have different physical meanings, both originate from two-phase -bubble and particle phase. The voidage signal can be divided into two parts, one corresponding to the bubble phase, and the other corresponds to the dense phase. Figure 3.10 is autocorrelation function figure for superimposed signals of Lorenz attractor and Logistic Map, with the time delay chosen as 1. When the embedding dimension is increased from 1 to 5, 9, and 14, a double slope curve appears as illustrated in Figures 3.11, 3.12, and 3.13. This demonstrates that a double slope correlation dimension curve can be detected from the signal from two superimposed independent chaotic signals. For the gas fluidisation case, double slope correlation dimension curves possibly originate from the interaction of two phases (bubble and particulate). In order to further analyse the chaotic structure of voidage signal, the voidage is divided into two parts corresponding to the bubble phase and dense phase for dynamic analysis. The cut-off point is selected at voidage 0.42, which corresponds to a voltage 53 P9 ln(Cd(r)) 3 CD Q . CD DJ C/) DJ 3 T J CD CL Tl C CQ" C Q_ -I CD 0 —i I _ CO 11 - i 0 CD 0 O o O" C 0 CT DJ CT *— 1 " O 13 CQ O 5' O CD 13 CQ g_ —i (—t-* DJ o" C/) —h CT) o DJ «—h O DJ T J CT) O DJ DJ 3 ^ T J CT » 0 C Q ^ O Q~ 0 0 C 0 CT) £ <5 O 13 ^ DJ ro o o I N I .U c 0.8 - \ o \ 'Z 0.6 \ o \ § 0.4 M— 0.2 _ c ° 0.0 1 i ™ -0.2 • 0 1 2 3 4 CO b -0.4 -o -" -0.6 o 3 -0.8 < -1.0 T I M E D E L A Y F i g u r e 3. 1 0 A u t o c o r r e l a t i o n f u n c t i o n f o r s u p e r i m p o s e d s i g n a l s o f L o r e n z a t t r a c t o r a n d L o g m a p output of 4.1 V for the fibre probe equipment. The results are shown in Figures 3.14 and 3.15 for sampling point O M 1 . Figure 3.14 presents the correlation integrals for a reconstructed optical fiber probe signal, assuming that all dense phase voidages are equal to the voidage at minimum fluidisation. Only one slope is shown in Figure 3.14. Figure 3.15 plots the corresponding correlation integrals assuming that all bubble phase voidages are equal to 1. The correlation integrals curve is then seen to show jumps, related to the reassignment of a discrete voidage of 1 to all voidages greater than 0.42. Clearly the assumption commonly made in bubble models that there are only two voidages (that at minimum fluidization and unity) is oversimplified. 3.6 Summary In this chapter, it is demonstrated that voidage time series from a bubbling fluidised bed demonstrate chaotic properties. Furthermore, it is shown that the signal can be divided into two portions, one part corresponding to periods when the probe is immersed in the bubble phase and the other when it is in the dense phase. ^ 5 ln(Cd(r)) o CD N 03 T l CQ" C —\ CD CO O o ro Qj CD' O =r. -* o 3 - a •S eg « 3 o 0 3 £- 5) CD_ TJ 0) -CO —5 CT) c 0) TJ _ (/) CD ~1I 0 O § & I t CD QL C O C/> o" Z3 o o s. I CQ o CO O) 1—1—1—1—r~ 1 4^ 1 ro 1 4^ C O I ro • B CL Q_ II II CO J 1 L ro 8S ln(Cd(r)) CD N 0) CD O Tl CQ' C 0 CO C O o o — \ CD 0) "1—'—r-O CO CD i ro ro i — 1 — r i CO o 13 03 Q_ ro 3 "O eg ^ 00 O 3 o g S- a> cp_ -o ^ 0) o «2. ^ 3 "O CD O Z5 <P Q. 3' CD 00 o" =3 0) CD CQ —i 0 0 0) 00 c T3 CD 3 o oo CD CL 00 CQ" QD_ O O 13 00 00* I— Z3* CQ i CO I ro CL Q. II II J i l i I i I i I i r i i i i i CO DJ 3 CD CL C 13 C L CD ln(Cd(r)) CD < O CL DJ CO CD DJ 3 c CD Q . CL N" CD DJ — i-t-DJ o II -, CO 2 | "D CD 0_ C T C D" C 7 5' CQ O O 13 o. I—f-' o" 13 CT) DJ CT) CT) C T l CQ" c —1 0 CO 4^  O o cp_ DJ o' 13 i i I I ^. _ \ . v X | | O) 4^ NJ O 00 CD cn I i | i | • | i | • | i 4^  CO 0 CQ —» DJ_ CV) I NJ 13 CQ r-»-ZT DJ 5- ^ CQ -+i — i 0 JD C 0 13 O II NJ O O N Q_ 0 CO 0 "O FT 0 O O CT) <—i-—i C o <—t-0 CL O TJ o' OJ DJ CT) 0 < O CL OJ CQ 0 CO OJ 0 0 jQ C OJ CT 0 O cr 0 CO C Q ' 13 OJ NJ NJ O NJ cn i i CO NJ • • CL CL II II ro ro 4^ o I • I J i L NJ 09 ln(Cd(r)) K> o ro Chapter 4 Bubble Coalescence Simulation Based on Clift and Grace Model and Dynamic Behaviour Analysis 4.1 Introduction As outlined in Chapter 2, recent experimental studies suggest that the dynamics of gas fluidised beds are determined by deterministic chaos. In this chapter, we demonstrate that deterministic chaos in bubbling fluidised beds can arise from non-linear bubble-bubble interactions. Starting from the Clift-Grace bubble coalescence model, a three-dimensional bubbling fluidised bed model is constructed that exhibits characteristic features of deterministic chaos. 4.2 Modelling of Three-Dimensional Bubbling Fluidised Bed It has been found experimentally that an isolated bubble in a fluidised bed rises with nearly the same velocity as a large bubble of the same volume ascending through a liquid of low viscosity, the general relationship being the equation of Davies and Taylor (1950). uA=Kjg~^ (1) where UA~ the rise velocity for an isolated bubble, g= the acceleration due to gravity, a= the bubble radius of curvature at the nose 61 Equation (1) was originally derived for bubbles rising in liquids, assuming that the flow at the front of the bubble is incompressible and irrotational and that the pressure is constant over the bubble surface. The liquid at the nose of the bubble is stagnant relative to the bubble. The theoretical values for the constant K are 2/3 for three-dimensional and 0.5 for two-dimensional bubbles. There is experimental evidence for three-dimensional fluidised beds, that K is not constant (Rowe, 1971) but varies with material and particle size. The average value is around 0.95, and almost all values lie within 15% of this value for the materials and sizes investigated. Rowe suggested that the best correlation for K is ^ = 1 . 3 4 ^ + 0 . 3 1 (2) where voidage at minimum fluidisation. The condition of constant pressure at the bubble surface used by Davies and Taylor applies equally in fluidised beds and controls bubble shape and velocity, both for isolated and interacting bubbles. This condition has been used to develop a model for predicting the motion of interacting bubbles in a fluidised bed (Gi f t and Grace, 1970). In this derivation a single doublet in uniform flow was used to represent each bubble in the potential flow equations. This method does not give spherical or circular shapes for interacting bubbles, but it has been shown for a vertical row of bubbles in two and three dimensions that the deviations are small and qualitatively in agreement with the distortions observed experimentally, i.e. the trailing bubble and doublet become elongated. Thus it is not necessary to represent the bubbles by completely spherical or circular shapes in order to predict their movement. 62 Gi f t and Grace showed that each bubble velocity could be predicted as the vector sum of the rise velocity in isolation and the velocity at the nose caused by each other bubble as predicted by potential flow with each bubble represented by a doublet. This is the basis used to formulate the model which allows us to predict the motion of interacting bubbles: "The velocity o f a bubble in a fluidised bed may be approximated by adding to its rise velocity in isolation the velocity which the continuous phase would have at the position of the nose if the bubble were absent". The theoretical justification of the basic postulate is limited to the two-dimensional case, but the model was shown experimentally to also be applicable in three dimensions. Comparison with experimental results for selected bubble configurations, bubble pairs in vertical alignment, bubble pairs in oblique alignment and vertical bubbles (Gift and Grace, 1970, 1971, 1972) has shown excellent agreement with experimental data for both two- and three-dimensional bubbles. In using the model it is necessary to calculate the velocity of the particulate phase at the nose each bubble due to the motion of all other bubbles present. It has been shown that an isolated bubble can be treated approximately as a circular (2-D) or spherical (3-D) cap with an enclosed wake of particles completing the circle or sphere and travelling with the bubble. Rotational flow is assumed to be confined to the wake, with the motion outside the bubble and its wake well approximated by potential flow around a cylinder or sphere. A s in the derivation of the model, the potential flow approximation is used here to calculate velocities of the particulate phase due to interacting bubbles. The equations derived from the model and the representation of coalescence are treated in detail below. For the two-dimensional case, the predictions from the model are 63 compared with experimental results. The three-dimensional case is treated and predictions compared to experimental results from the literature. The Clift-Grace coalescence model was initially developed for a two-dimensional bed, and this treatment was then extended to three-dimensional beds by Johnsson (1973) and Johnsson et al (1974). The development follows the same pattern as in the two-dimensional case. The rise velocity for a single bubble in isolation and the interaction coefficients are different in three dimensions, but the logic concerning coalescing bubbles is the same in both cases and the control criteria for two dimensions can be used directly, a) Bubble Interaction i) The rise velocity of a bubble in isolation, u A , is calculated from the Davies-Taylor equation (1) using the theoretical isolated bubble three-dimensional value for K instead of equation (2), because equation (1) has already been shown to give good predictions of the motion of bubble pairs and bubble chains, while equation (2) was an average value for interacting bubbles, i.e. inappropriate for the detailed model: uA = 0.67VP^ (3) ii) The velocity of the particulate phase caused by the motion of all other bubbles, evaluated at the position of the nose of the bubble in question, but with that bubble absent, is calculated from potential f low theory assuming that each bubble is represented oy a doublet. The velocity is split into one vertical component q, and two horizontal velocity components, p and a.The velocity components for the particulate phase at the nose of bubble i caused by bubble j are: 64 9v=rVUJ+t(fVJ+l9WJ (4) P o = t U U J + S V V J + m V W J (5) crij=lijuj+mijvj+nijwj * (6) where qij =vertical particulate phase velocity at the nose of bubble i caused by bubble j Pi, =horizontal particulate phase velocity in the y-direction at the nose of bubble i caused by bubble j Oy =horizontal particulate phase velocity in the z-direction at the nose of bubble i caused by bubble j Uj ^vertical instantaneous velocity of bubble j Vj =horizontal instantaneous velocity of bubble j in the y-direction Wj=horizontal instantaneous velocity of bubble j in the z-direction Here x,y and z are Cartesian co-ordinates of the bubble centre, with x directed vertically upwards; r,j, Sy, n U i ty, 1^ m yare interaction coefficients depending on the size and position of bubbles i and j as defined in Table 4.1 The instantaneous velocity components for bubble i can now be calculated by summation: N N V 2>y ( 8 ) N W .=!>.> (9) 65 T A B L E 4.1 Interaction Coefficients for the Prediction of Bubble Mot ion in Three Dimensions (adapted from Johnsson, 1973) fy = a)(2X2 -Y2 -Z2)lAB stj - a*(2Y2 - X2 - Z2) I AB nij=a)(2Z2 -X2 -Y2)l AB tv = 3a)XYIAB lv =3a].XZ/AB mv=3a*YZ/AB where X = xt+ at - x. Z = z i - zj AB = 2(X2 +Y2+Z2f2 The bubble trajectories can then be found by integration of the three differential equations: dxi dyi dzt a- = "':*- = v' :*- = "'' <I0> Equations (7) to (10) cover the case where no wake entry has occurred. If the nose of bubble i has entered the wake of bubble j , we have simply: 66 Ui = UAi +Uj v, = vj (11) (12) (13) 4.2.1.Coalescence As in the two-dimensional case and as proposed by Clift and Grace (1972) , the following assumptions are made at the instant of coalescence: i) Bubble volumes are either additive upon coalescence, or the total volume increases by 10-15%; increases of this magnitude have been seen in practice (Grace and Venta, 1973). ii) Bubbles are close to vertical alignment when coalescing. iii) uAj is equal to the rise velocity of the composite bubble, and the interaction coefficients Sji nji tjj lji and nij; are all zero. iv) The time for completion of coalescence is 1.8 a/uAj-4.2.2 Numerical Solution a) General Method The numerical solution is carried out as in the two-dimensional case, the only difference being that N bubbles give 3N linear equations in three dimensions, rather than 2 N as in two dimensions. The solution is straightforward since the coefficent matrix has the largest elements along the main diagonal. Hence a simple iterative procedure can be used. The Gauss-Seidel method is used to calculate the velocities. The position of each bubble as a function of time is then obtained by integrating the 3 N differential equations simultaneously by the Kutta-Merson process. The integration step is chosen to be small enough to ensure that the check on bubble coalescence between steps is carried out at convenient intervals to prevent unrealistic situations, such as bubbles enclosing other bubbles. b) Wal l Effect Image bubbles were used to represent the walls of a square column. Two simulations were run with nine bubbles in the bed. In one case no image bubbles were included, while in the other image bubbles were included for all bubbles whose centres were less than two bubble radii from the walls: in this case five image bubbles were present. The largest effect was observed for the vertical velocity; one bubble was retarded slightly to 98.5% of the original velocity, while the effect on horizontal velocities was completely negligible. Consequently the wall effect on bubble velocities was assumed to be negligible and the image bubbles were omitted in the simulations. This was fortunate because it enables simulation in beds of circular cross-section also, without additional difficulties. c) Boundary Condition The boundary condition at the gas distributor (i.e. at the base of the bed) is usually taken as bubbles of uniform size, evenly-spaced in time, randomly distributed over the bed cross-section, but totally inside the bed, i.e. no bubble centre is closer than one bubble radius to the outer wall of the column. In the square bed the random distribution was obtained by picking both horizontal co-ordinates randomly. For the bed of circular cross-68 section the procedure was the same, but the position of each bubble was tested to see if it was inside the circle of radius R-rb0, where rb0 is the initial bubble radius. If not , it was discarded and a new bubble tried, d) Further Computational Details Bubbles were taken out of the bed at the top when their centres passed the level of the expanded bed height. If a bubble happened to pass through the wall during an integration step, it was put back into the bed after the step. This could happen because of the neglect of the wall in the simulation for high bubble concentrations. The horizontal positions and size of each bubble were recorded at given levels. This computer program is given in the Appendix. This program is derived from that of Johnsson (1973) for a Fortran M S Power Station compiler with an improved random number generator. 4.3. Signal Simulation The bubble information was used to determine the average voidage at selected positions in the bed. Observation "windows" were defined as horizontal circles. The instantaneous voidage fraction at each location is then determined after each time step by considering the array o f bubbles at that instant as predicted by numerical solution of the differential equations for the bubbles, as discussed above. The voidage fraction calculated in this manner led to predicted time series measurements for each window, corresponding to what a stationary observer or voidage sensor should see at that position. Obviously, the area of the "window" has some influence on the nature of the voidage fraction. It has been found empirically that the most useful windows have diameters which are 0.5 to 10 times the bubble diameter. 69 4.4. Chaos Analytical Method 4.4.1 Phase-Space Plots A phase-space plot is a graph in which the derivative X ' is plotted versus X at each data point in order to reveal the topology of the solution. For this purpose, the derivatives are assumed to be given by half the difference between the data points adjacent to each point divided by the X interval. Periodic data should appear as a closed curve on such a plot. With the special embedding method, chaotic data often appear in the form of a strange attractor having a fractal structure with a fractional dimension. 4.4.2 Lyapunov Exponent The Lyapunov exponent is a measure of the exponential rate at which nearby trajectories diverge in phase space. Chaotic orbits have at least one positive Lyapunov exponent, while for periodic orbits, all Lyapunov exponents are negative. The Lyapunov exponent is zero near bifurcation points. In general, there are as many exponents as there are dynamical equations. The most positive exponent is chosen as the Lyapunov exponent. 4.4.3. Autocorrelation Function The autocorrelation function, sometimes called the correlation function, is obtained by multiplying each X(t) by X(t-x) and summing the result over all data points. The sum is then plotted as a function of time delay x. This gives a measure o f how dependent data points are on their neighbours. The value of x at which the correlation function remains small is defined as the correlation time. Highly random data have no 70 correlation, and for these data the correlation function drops abruptly to zero, implying a small correlation time. Highly correlated data have a correlation amplitude that decreases only slowly with time. Chaotic data tend to show little correlation except when the time delay is small. 4.4.4 Correlation Dimension With each pass through the data, a new data point is taken, and a hyperdimensional sphere of embedding dimension d and radius r is centered on that point. The fraction of subsequent data points in the record within that sphere is then calculated for various values of r, and a plot is made of the log of this number versus the log of the radius. The correlation dimension is taken as the average slope of the cumulative curve over the middle one-quarter of the vertical scale, while the error is taken as half the difference between the maximum and minimum slope over the same range. 4.5 Sample Results And Discussion Because this model does not have any way of predicting bubble splitting, it is not possible to use it to simulate the bubble column used in the experiments reported in Chapter 3. Instead, the initial conditions were set as: The initial conditions were set as: initial bubble diameter = 10 mm bubble generation frequency = 100 H z column diameter = 500 mm static bed height = 1 m 71 time step = 0.05 s A "window" of diameter 10 mm was positioned with its center on the axis at height of 0.5 m above the distributor. The overall number of time steps of the computer run is 10,000, which required about 30 minutes to obtain results on a Pentium 100 H z computer. The simulated real time is 500 seconds. The program created time series data of bubble size, radial position and axial position. The voidage fraction is achieved by numerically calculating the area of three-dimensional bubbles in the fluidised bed cut by the horizontal window. The voidage of the cut area of the bubbles in the window is assumed to be l , while the remaining part of the window is assumed to have voidage of Smf. The average voidage fraction of the window is then assumed to be like the voidage signal obtained by an optical probe in the fluidised bed. The data are then used to remap the time series measurement into a phase or state. The phase space plot of the simulated voidage fraction signal is shown in Figure 4.1. It is clear that a strange attractor reveals in the phase space plot. The correlation function is shown in Figure 4.2. For time delay 0.05, embedding dimension 5 and 9, a single slope correlation integer is shown. 72 de/dt -5 -4 -3 -2 -1 ln(r) 0 1 2 Figure 4 2. Correlation integrals for simulated voidage fraction signal ( x=0.05) The correlation dimension for the time series is Dc=2.1 for the embedding dimension of 9. The largest Lyapunov exponent is equal to 0.182. A second set of conditions was also investigated: initial bubble diameter =15 n m bubble generation frequency = 100 H z column diameter = 500 mm static bed height = 1 m 74 time step = 0.05 s A n imaginary "window" of the same size with diameter 10 mm was again positioned centred on the axis at a height of 0.5 m above the distributor. 10,000 time steps of the computer run were taken, resulting in an overall simulated time of 500 seconds. The program converged and created time series data of bubble size and position. The resulting correlation function is shown in Figure 4.3. -5 -4 -3 -2 -1 0 1 2 | I 1 I 1 1 ' 1 > 1 1 1 1 1 1 1 I 1 1 1 1 1 i I i I -5 - 4 - 3 - 2 - 1 0 1 2 -l n ( r ) Figure 4 3. Correlation function figure for simulated voidage fraction signal ( x=0.05) The resulting correlation dimension for the time series is given by D c=2.3 while the embedding dimension is 9. The Lyapunov exponent is equal to 0.196. The third set of conditions was set as: initial bubble diameter = 10 mm bubble generation frequency = 120 H z 75 column diameter = 500 mm static bed height = 1 m time step = 0.05 s The same "window" was positioned at the same position as in the first two cases. Also the total number of time steps of the computer (10000), and the total simulated time (500 s) were as before. The program again converged and created a time series data of bubble size and position. Figure 4.4 gives the correlation function from the predicted signal for two values of the embedding dimension, d. • 6 - 4 - 3 - 2 - 1 0 1 2 "1 ' I 1 I 1 1 > 1 1 1 ' 1 1 1 I n ( r ) Figure 4 4. Correlation function figure for simulated voidage fraction signal (7=0.05) The correlation dimension for the time series is D c=2.8 while the embedding dimension is 9. The Lyapunov exponent is equal to 0.252. 76 The above three simulations, summarised in Table 4.2, demonstrate that dynamic chaos dominated this simulated 3-dimensional fluidised bed. B y comparing these three se of results, it can be seen that both bubble frequency and initial bubble size influence the chaotic characteristics o f the simulated voidage fraction signal, but bubble frequency appears to play a stronger role than the initial bubble size. .2 Summary of results of numerical simulation using three-dimensional bubble coalescence model. Simulation 1 Simulation 2 Simulation 3 Column diameter, mm 500 500 500 Static bed height, m 1 1 1 Initial bubble diameter, mm 10 15 10 Bubble generation frequency, H z 100 100 120 Time step, s 0.05 0.05 0.05 Figure showing correlation 4.2 4.3 4.4 function Correlation dimension, D 2.1 2.3 2.8 Lyapunov exponent, X 0.182 0.196 0.252 7 7 Chapter 5 Conclusions and Recommendations 5.1 Conclusion Recent research suggests that the dynamics o f gas fluidised beds are dominated by deterministic chaos. This thesis investigates the chaotic phenomena by examining experimental voidage signals from an optical fibre probe in a freely bubbling fluidised bed and simulated voidage fractions in a bubbling fluidised bed, predicted based on the three-dimensional version of the Clift and Grace bubble coalescence model It is shown that a bubbling fluidised bed can behave as a multi-fractal dynamic system. Two-slope correlation integral curves derived from experimental voidage signals possibly arise from superposition of the two curves corresponding to the bubble and dense phase. Voidage signals collected at different positions show that the bubble size and frequency influence the chaotic characteristics of the voidage signal. In the bottom section, the voidage signals generally are more chaotic than in the top section. In addition, near the axis of the top sections, voidage signals are less chaotic because of bubble coalescence. Simulated voidage fraction signals derived from the bubble coalescence also show chaotic characteristics, in qualitative agreement with experimental results. Both the experiments and numerical simulation indicate that the bubble frequency influences chaotic characteristics more strongly than the bubble size. Voidage signals from regions of high bubble frequency region are more chaotic than those in regions where there is 78 little coalescence. The initial size of bubbles appears to have only a weak influence on the signal's chaotic characteristics 5.2 Recommendations Though some progress in understanding the fundamentals of the chaotic structure of voidage signals in bubbling fluidised beds has been made in this study, more research is required both experimentally and theoretically to completely understand the chaotic behaviour of voidage signals. The influence of the non-linear relationship between voidage and experimental voltage optical fibre probe output for small particles on the derived signal chaotic characteristics is unknown. Research is needed to reach a conclusion as to whether the non-linear output influences chaotic index of signals, and i f it does have an influence, how strong the influence is. The multi-fractal chaotic structure assumption needs more direct proof. The numerical box counting method needs careful review to make sure that the double slope correlation integral curve is consistent and sufficient to support this assumption. In order to explore the elementary form of deterministic chaos in bubbling fluidised beds, further effort is needed to assume every point in a bubbling fluidised bed to consist of either void (bubble) or particulate phase. Imaginary optical fibre probes would then be positioned in the bed at different locations. A set of time series voidage values would then be generated. Because the output wi l l be an analogue signal, far more data points wi l l be needed to explore the dynamics of the simulated signals. Bubble splitting needs to be added to the Cliff-Grace model. 7 9 N O M E N C L A T U R E A constant, -A bed cross-sectional area, m 2 Cd(r) correlation integral corresponding to a specific embedding dimension, -C w (r) correlation integral corresponding to a cutoff parameter W, -d embedding dimension, -db radius of bubble, m D c correlation dimension of attractor, -Dcap capacity dimension of attractor, -Di information dimension of attractor, -f frequency, s"1 G B visible gas flow in the bubble phase, m 3/s K Kolmogorov entropy, -k number of data points of a subrecord, -M order of nearest neighbor or vector of time series, -m number of time steps, -N number of data points, -r distance, m r0 specific distance beyond which two nearby trajectories separate, -R range of cumulative time series, -t time, s tac correlation time, -tc characteristic time scale, -tet evolution time, -U superficial gas velocity, m/s Umf superficial gas velocity at minimum fluidisation, m/s U m b superficial gas velocity at minimum bubbling point, m/s Ubr rise velocity of isolated bubble, m/s X variable, -Y variable, -z axial coordinate, m X reconstructed vector, -5 nearest neighbour distance, -Sij distance between points Xt and X j , -9 dimensionless radial coordinate, -a standard deviation, -X time delay, -K largest Lyapunov exponent, -Xet Lyapunov exponent corresponding to evolution time t et, -y autocorrelation function, -Ymax maximum value of autocorrelation function, -REFERENCES Albano, A . M , Muench, J . , Schwartz, C , Mees, A . I. and Rapp, P. E., Singular-value decomposition and the Grassberger-Procaccia algorithm, Phys. Rev. A, 38, 3017-3026, 1988. Albano, A . M , Passamente, A . and Farrell, M . E., Using higher-order correlations to define an embedding window; Physica D, 54, 85-97 1992. Aleksic, Z . , 1991, Estimating the embedding dimension, Physica D, 52, 362-368. Badii, R. and Polit i, A . , 1985, Statistical description of chaotic attractors: the dimension function, / . Stat Phys., 40, 725-750. Bai , D., B i , EC. T. and Grace, J. R., Chaotic behavior of fluidized beds based on pressure and voidage fluctuations, A I C h E J. , 43, 1357-1361 (1997). Bai , D., Issangya, A . I., L im, K. S. and Grace, J. R., Gas fluidized beds: high-dimensional and multifractal dynamic systems? AIChE Annual Meeting, November 10-16, Chicago, 1996. Barnsley, M . , Fractals Everywhere, Academic Press, San Diego, C A , 1988. Ben Mizrachi, A . , Procaccia, I. and Grassberger, P., Characterization of experimental (noisy) strange attractors, Phys. Rev A, 29, 975-977, 1984. Bouillard, J. X . and Miller, A . L., Experimental investigations of chaotic hydrodynamic attractors in circulating fluidized beds, Powder Technoi, 79, 211-215., 1994. Clift R. and Grace J.R., 1970, Bubble interactions in fluidised beds, Client Eng. Prog. Symp. Ser., 66, No . 105, 14-27. Clift R. and Grace J R . , 1971, Coalescence of bubbles in fluidised beds, AIChE Symp. Ser., 67, No . 116, 23-33. Clift R. and Grace J.R., 1972, Coalescence of bubble chains in fluidised beds Trans. Inst. Chem. Eng, 50, No.4, 364-371. Clift R. and Grace J.R., Continuous bubbling and slugging , Chapter 3 in Fluidization, 2nd ed., J.F. Davidson, R. Clift and D. Harrison.,eds., Academic Press., London, 1985. Collins R., The Rise Velocity of Davidson's Fluidization Bubble, Chem, Eng. Sci., 20, 788(1965). 82 Davidson J.F. and Harrison D., Fluidized Particles, Cambridge Univ. Press.,Cambridge, 1963. Davidson J.F., Paul R.C. , Smith M . J.S., Duxburg H.A., The rise of bubbles in a fluidised beds, Trans. Inst. Chem. Eng., 37, 323 (1959). Davidson, J.F, Harrison D. and Guedes de Carvalho, J.R.F., On the liquid like behaviour of fluidised beds., Ann. Rev. Fluid Mech., 9,55 -86. Davies R . M . and Taylor G.I., Proc. Roy. Soc, A 2 0 0 , 375 (1950). Daw, C. S. and J. S. Halow, Characterization of voidage and pressure signals from fluidized beds using deterministic chaos theory, Proc. 11th International Fluidized Bed Combustion Conference, A S M E , 2 , 777-786, 1991. Daw, C. S. and Halow, J. S., Modeling deterministic chaos in gas-fluidized beds, AIChE Symp. Ser., 88, No.289, 61-69., 1992. Daw, C. S. and Halow, J. S., Evaluation and control of fluidization quality through chaotic time series analysis of pressure drop measurements, AIChE Symp. Ser., 89(296), 103-122., 1993. Daw, C. S., Lawkins, W. F., Downing, D. J. and Clapp, N. E., Chaotic characteristics of a complex gas-solids flow, Phys. Rev., A, 41, 1179-1181, 1990. Devaney R. L., A n Introduction to Chaotic Dynamical Systems, Addison- Wesley, 1989. Drahos, J . , Bradka, F. and Puncachar, M . , Fractal behaviour o f pressure fluctuations in a bubble column, Chem. Eng. Sci., 47, 4069-4075, 1992. Fan, L. T., Kang, Y . , Neogi, D. and Yashima, M . , 1993, Fractal analysis of fluidized particle behaviour in liquid-solid fluidized beds, AIChE J., 39, 513-517. Fan, L. T., Neogi, D., Yashima, M . and Nassar, R., 1990, Stochastic analysis of a three-phase fluidized bed: fractal approach, AIChE J., 36, 1529-1535. Farmer, J. D. and Ott, E., 1983, The dimension of chaotic attractors, Physica D, 7, 153-180. Feder, J . , 1988, Fractals, Plenum, New York. Franca, F., Acikgoz, M . , Lahey Jr, R. T. and Clausse, A . , 1991, The use of fractal techniques for flow regime identification, Int. J. Multiphase Flow, 17, 545-552. Fraser, A . M . and Swinney, H. L., 1986, Independent coordinates for strange attractors from mutual information, Phys. Rev. A., 33, 1134-1140. 8 3 Geldart, D., The size and frequency of bubble in two- and three- dimensioned gas fluidised beds, 1970, Powder Technology, 4, 41 Geldart, D., 1973, Types of gas fluidisation, Powder Technology, 7, 285-292. Glicksman. L. R., 1984, Scaling relationships for fluidized beds, Chem. Eng. Sci., 39, 1373-1379. Glicksman. L. R., 1988, Scaling relationships for fluidized beds, Chem, Eng. Sci., 43, 1419-1421. Glicksman, L. R., Hyre, M . R. and Farrell, P. A . , 1994, Dynamic similarity in fluidization, Int J. Multiphase Flow, 20, 331-386. Gogolek, P. E. and Grace, J. R., 1995, Fundamental hydrodynamics related to pressurized fluidized bed combustion, Prog. Energy Combust. Sci., 21, 419-451. Grace, J. R., 1984, Recent Advances in Engineering Analysis of Chemically Reacting Systems, L. K. Doraiswamy (Ed.), pp.237-235, Wiley Eastern, New Delhi. Grace, J. R., 1990, High-velocity fluidized bed reactors, Chem. Eng. Sci., 45, 1953-1966. Grace J. R. and Clift R. 1974, On the two-phase theory of fluidization, Chem. Eng. Sci., 29, 327-334, (1974). Grace, J.R. and Harrison, D,1968, The distribution of bubbles within a gas-fluidised bed, Instn, Chem. Engrs. Symp. Ser. No . 30, 105-113. Grace, J.R. and Venta, J (197#) Volumes Changes accompanying bubble splitting in fluidised beds, Can. J. Chem. Eng., 51, 110-111. Grassberger, P. and Procaccia, I., 1983a, Estimation of the Kolmogorov entropy from a chaotic signal, Phys. Rev., A, 28, 2591-2593. Grassberger, P. and Procaccia, I., 1983b, Measuring the strangeness of strange attractors, Physica D, 9, 189-208. Grassberger, P., Schreiber, T. and Schaffrath, C , 1991, Nonlinear time sequence analysis, Int J. Bifurcation and Chaos, 1, 521-547. Halow, J. S. and Daw, C. S., 1994, Characterizing fluidized-bed behaviour by decomposition of chaotic phase-space trajectories, AIChE Symp. Ser., 90(301), 69-91. Hammel, S. M . , 1990, A noise-reduction method for chaotic systems, Phys. Lett A, 148, 421-428. Harrison D. and. Leung L .S , The rate of rise of bubbles in fluidized beds, Trans. Inst Chem. Eng., 39, 409 (1961); 40, 146-151 (1962) 84 Hay, J. M . , Nelson, B. H. , Briens, C. L. and Bergougnou, M . A . , 1995, The calculation of the characteristics of a chaotic attractor in a gas-solid fluidized bed, Chem. Eng. Sci., 5 0 , 373-380. He, Y . , Hydrodynamic and scale-up studies of spouted beds, PhD Thesis, University of British Columbia, 1995. Hilborn, R. C , 1994, Chaos and nonlinear dynamics: an introduction for scientists and engineers, Oxford University Press, New York. Issangya, A . , Bai . , D., B i , H.-T., L im, K. S., Zhu, J . -X . and Grace, J. R , Axial solids holdup profiles in a high-density circulating fluidized bed, in Circulating Fluidized Bed, Technology V , M . Kwauk and J. L i , eds., Science Press, Beijing, pp.60-65, 1997. Izrar, B. and Lusseyran, F., 1993, Chaotic behaviour of an annular film of liquid unstabilized by an interfacial shear stress, Instabilities in Multiphase Flows (Edited by G. Gouesbet and A . Berlemont), pp. 1-15, Plenum Press, New York. Jackson R.,1963, Bubble velocity in fluidised beds Trans. Inst. Client Eng., 4 1 , 22-28. Johnsson, J.E. , Bubble distribution in fluidised bed, M.Eng. Thesis, McG i l l University, 1973. Johnsson, J .E. , Gi f t , R. and Grace J. R., 1974, Prediction of bubble distributions in a freely-bubbling two-dimensional fluidised beds, Instn. Chem, Engrs. Symp. Ser., No.38, 85. Kennel, M . B. and Isabelle, S., 1992, Method to distinguish possible chaos from colored noise and to determine embedding parameters, Phys. Rev. A, 4 6 , 3111-3118. Kennel, M . B., Brown, R. and Abarbanel, H. D. I., 1992, Determining embedding dimension for phase-space reconstruction using a geometrical construction, Phys. Rev. A, 4 5 , 3403-3411. Kostelich, E. J. and Schreiber, T., 1993, Noise reduction in chaotic time-series data: a survey of common methods, Phys. Rev. E, 4 8 , 1752-1763. Kostelich, E. J. and Swinney, H. L., 1987, Practical considerations in estimating dimension from times series data, in Chaos and Related Nonlinear Phenomena, Plenum, New York. Kostelich, E. J. and Yorke, J. A . , 1988, Noise reduction in dynamical systems, Phys. Rev. A, 3 8 , 1649-1652. Kostelich, E. J. and Yorke, J. A . , 1990 Noise reduction in dynamical systems, Physica D, 4 1 , 183-196. 85 Lacy, C. E. , Sheituch, M . and Dukler, A . E., 1991, Methods of deterministic chaos applied to the flow of thin wavy films, AIChE J., 37, 481-489. Liebert, W . and Schuster, H . G. , 1988, Proper choice of the time delay for the analysis of chaotic time series, Phys. Lett A, 142, 107-111. Lischer, D. J. and M . Y . Louge, Optical fiber measurements of particle concentration in dense suspensions: calibration and simulation, Appl. Optics, 31(24), 5106-5113 (1992). Luewisuthichat, W., Tsutsumi, A . and Yoshida, K., 1995, Fractal analysis of particle trajectories in three-phase systems, Trans. IChentE, 73 (Part A ) , 222-227'. Matsuno, Y . , H. Yamaguchi, T. Oka, H. Kage and K. Higashitani, The use of optic fiber probes for the measurement of dilute particle concentrations. Calibration and application to gas-fluidized bed carryover, Powder Technol., 36,215-221 (1983). Moon, F.C., 1992, Chaotic and Fractal Dynamics, John Wiley & Sons Inc., New York. Mull in, T., 1993, Nature of Chaos, Oxford Science Publications, Oxford. Murray J.D., On the mathematics of fluidization, / . Fluid Mech., 22, 57-80 (1965) Nguyen, X . T. and Leung, L. S. ,A note on bubble formation at an orfice in a fluidized Bed, Chem. Eng. Sci., 27, 1748-1750 (1972). Olofsen, E., Degoede, J. and Heijungs, R., 1992, A maximum likelihood approach to correlation dimension and entropy estimation, Bull. Math. Bio., 54, 45-58. Olowson P. A . and Almstedt A.E.,Hydrodynamics of a Bubbling Fluidized Bed, Client Eng. Sci., 47,357-366 (1992). Parker T. S. and Chua, L. O., 1989, Practical Numerical Algorithms for Chaotic Systems, Springer-Verleg, New York. Parlitz, U. , 1992, Identification of true and spurious Lyapunov exponents from time series, Int J. Bif. Chaos, 2, 155-165. Pence, D. V . , Beasley, D. E. and Riester, J. B., 1995, Deterministic chaotic behaviour of heat transfer in gas fluidized beds, J. of Heat Transfer, 117, 465-472. Pianarosa, D, Hydrodynamic studies of spouted and spout-fluid beds, M A S C Thesis, University of British Columbia, 1996 Rowe P.N. , in Fluidization, J.F. Davidson and D. Harrison, eds., 121, Academic Press, London, 1971 Rowe P.N, Partridge B A , Trans. Inst. Chem. Eng., 43, T157(1965) Sano, M . and Sawada, Y . , 1985, Measurement of the Lyapunov spectrum from chaotic time series, Phys. Rev. Lett., 55, 1082-1085. 86 Sato, S., Sano, M . and Sawada, Y . , 1987, Practical methods of measuring the generalized dimension and the largest Lyapunov exponent in high dimensional chaotic systems, Prog. Theor. Phys., 77, 1-5. Schouten, J. C. and van den Bleek, C. M . , 1991, Chaotic behaviour in a hydrodynamic model of a fluidized bed reactor, in Proc. 11th International Fluidized Bed Combustion Conference, pp.459-466, A S M E , New York. Schouten, J . C. and van den Bleek, C. M . , 1992, Chaotic hydrodynamics of fluidization: consequences for scaling and modeling of fluidized bed reactors, AIChE Symp. Ser., 88, No.289, 70-84. Schouten, J. C , Takens F. and van den Bleek, C. M . , 1994, Maximum-likelihood estimation of the entropy of an attractor, Phys. Rev. E, 49, 126-129. Schouten, J. C , van der Stappen, M . L. M . and van den Bleek, C. M . , 1992, Deterministic chaos analysis of gas-solids fluidization, in Fluidization VII (Ed. by O. E. Potter and D. J. Nicklin), pp. 103-111, Engineering Foundation, New York. Schreiber, T. and Grassberger, P., 1991, A simple noise reduction method for real data, Phys. Lett. A, 160, 411-418. Skrzycke, D. P., Nguyen, K. and Daw, C. S., 1993, Characterization of the fluidization behaviour of different solid types based on chaotic time series analysis of pressure signals, in Proc. of 12th International Fluidized Bed Combustion Conference, Vo l .1 , pp. 155-166, A S M E , New York. Stewart P.SB, Iso la ted Bubbles in a Fluidized Beds - Theory and Experiments, Trans. Inst. C h e m . Eng. , 46, T60-T66 (1968) Takens, F. 1981, in Lecture Notes in Mathematics, Vol.898, Springer, New York. Tarn, S. W. and Devine, M . K., 1991, Is there a strange attractor in a fluidized bed? in Measures of Complexity and Chaos (Ed. by N . B. Abraham, A . M . Albano, A. Passamante and P. E. Rapp), pp. 193-197, Plenum Press, New York. Toei R , 1967, in Proc. Int Symp. On Fluidization, A A . H Drinkenburg, ed., p.271-278, Netherlands Univ. Press, Amsterdam Theiler, J . , 1986, Spurious dimension from correlation algorithms applied to limited time-series data, Phys. Rev. A, 34, 2427-2432. Tsuji, Y . , Morikawa, A . and Shiomi, H., 1984, L D V measurement of an air-solid two phase flow in a vertical pipe, J. FluidMech., 139, 417. Valenzuela J.A. and Glicksman L.R., Gas flow distribution in a bubbling fluidized bed, Powder Tech., 44, 103-113 (1985) 87 van den Bleek C. M . and Schouten, J. C , 1993a, Can deterministic chaos create order in fluidized-bed scale-up? Chem. Eng. Sci., 48, 2367-2373. van den Bleek, C. M . and Schouten, J. C , 1993b, Deterministic chaos: a new tool in fluidized bed design and operation, Chem. Eng. J., 53, 75-87. van der Stappen , M . L. M . , Schouten, J. C , and van den Bleek, C. M . , 1993b, Application of deterministic chaos theory in understanding the fluid dynamic behaviour of gas-solids fluidization, AIChE Symp. Ser., 89, No.296, 91-102. van der Stappen, M . L. M . , Schouten, J. C. and van den Bleek, C. M . , 1993a, Deterministic chaos analysis of the dynamical behaviour of slugging and bubbling fluidized beds, in Proc. 12th International Fluidized Bed Combustion Conference, Vol .1 , pp. 129-140, A S M E , New York. van der Stappen, M . L. M . , Schouten, J. C. and van den Bleek, C. M . , 1994, Application of deterministic chaos analysis to pressure fluctuation measurements in a 0.96 m 2 C F B riser, in Circulating Fluidized Bed Technology IV (EA. by A . A . Avidan), pp.55-60, A IChE , New York. van der Stappen, M . L. M . , Schouten, J. C. and van den Bleek, C. M . , 1995, Chaotic hydrodynamics and scale-up of gas-solids fluidized beds, Fluidization VIII Preprints, pp.625-632. Werther, J . , Fluid mechanics of large-scale C F B units, in Circulating Fluidized Bed Technology IV, A. A. Avidan (ed.), 1-14 (1994). Wolf, A . , Swift, J. B., Swinney, H. L. and Vastano, J. A . , 1985, Determining Lyapunov exponents from a time series, Physica D, 16, 285-317. X ia , Y . , Zheng, C. and L i , H., 1992, Characterizing fast fluidization by optic output signals, Powder Technol., 72, 1-6. Yamazaki, H. , K. Tojo and K. Miyanani, Measurement of local concentration in a suspension by an optical method, Powder Technol., 70, 93-96 (1992). Yates, J .G. and Cheesman, D J , Voidage variations in the regions surrounding a rising bubble in a fluidized bed, AIChE Symp. Ser., No. 289, 88, 34-39 (1992). Yates, J .G. and Rowe, P.N. , The effect of pressure on the flow of gas in fluidized beds of fine particles, Chem. Eng. Sci., 38, 1935 -1945(1983). Zeng, X . , Eykholt, R. and Pielke, R. A . , 1991, Estimating the Lyapunov exponent from short time series of low precision, Phys. Rev. Lett., 66, 3229-3232. Zeng, X . , Pielke, R. A . and Eykholt, R., 1992, Extracting Lyapunov exponent from short times series of low precision, Mod Phys. Lett. B, 6, 55-75. 88 Appendix 3-DIMENSIONAL FLUIDISED BED SIMULATION: F o r t r a n Computer Code (adapted from Johnsson, 1973) C BUBBLE COALESCENCE SIMULATION IN FLUIDIZIED BED DIMENSION A(1000,1000), Y(1000), DY(1000), R(1000) COMMON /DEP/Y/YATC/A/GRAD/DY/RAD/R DIMENSION CLEV(6),NUM(6) COMMON /CLE/CLEV/NUM8/NUM C A: COEFFIENT MATRIX C Y: CARTECISIAN COORDINATES OF BUBBLE CENTER C DY: VELOCITY COMPONENTS C R: BUBBLE RADII C N: BUBBLES GIVEN IN 3-DIMENSION SIMULATION C o p e n ( u n i t = 2 , f i l e = ' s i m l . d a t ' , s t a t u s = 1 unknown 1) o p e n ( u n i t = 2 , f i l e = ' s i m l l . d a t ' , status='unknown') o p e n ( u n i t = 4 , f i l e = 1 s i m l 2 . d a t ' , status='unknown') open(unit=5, f i l e = 1 s i m l 3 . d a t 1 , status='unknown 1) o p e n ( u n i t = 3 , f i l e = 1 i n i l . d a t ' , s t a t u s = ' o l d ' ) WRITE(*,*)' INPUT ARGUMENTS FOR MERSON' READ(3,*) X1,DELX1,DX1,DXMIN1,T0LKM1 WRITE(*,*)' BED HEIGHT AND WIDTH 1 READ(3,*) BEDH,BEDW WRITE(*,*) 1 INPUT RECORD LEVEL' READ(3,*) (CLEV(K),K=l,6) 101 FORMAT(8F10.5) 1011 FORMAT(6F10.5) WRITE(*,*) 'INPUT BUBBLE SIZE AND FRQUENCY' • READ (3,*) RM,LOOP WRITE(*,*) 1 INPUT ARGUMENTS FOR SIDI' READ (3,*) ABC,NIK 102 FORMAT(F10.2,15) WRITE(6,201) X1,DELX1,DX1,DXMIN1,T0LKM1 201 FORMAT('ARGUMENTS FOR MERSON',/,' XI',F10.5,'DELX1=',Fl0.5 $ ,'DX1=',F10.5,'DXMIN1=',F10.5,'TOLKMNl=',F10.5) WRITE(6,301) BEDH,BEDW 301 FORMAT(2X,'BED HEIGHT,CM=',F6.2,'BED WIDTH,CM=',F6.2) WRITE(6,302) 302 FORMAT('HEIGHT OF RECORDING LEVELS') WRITE(6,303) (CLEV(K),K=l,2) 303 F0RMAT(1X,2F8.2) WRITE(6,304) RM,LOOP 304 FORMAT('ARGUMENTS FOR SUBROUTINE BEDBO',/,'INITIAL BUBBLE $ SIZE,CM=',F6.2,180,'FREQUENCY PARAMETER=',14) WRITE(6,305) ABC,NIK 305 FORMAT('ARGUMENTS FOR SUBROUTINE SIDI: ABC=',F10.2,/'NUMBER $ OF STEPS IN MAIN PROGRAM: NIK=', 15) IB=0 N=0 DO 5 1=1,6 5 NUM(I)=0 DO 40 1=1,NIK C BOUNDARY CONDITION AND BUBBLE SIZE AND POSITION AT X = 0 CALL BOUNDC(N,RM,BEDW,LOOP) 89 C RECORD BUBBLES PASSING GIVEN LEVELS CALL CLEVS(N,XI,ABC,RM) C WHEN A BUBBLE LEAVES THE BED, DECREASE N BY 3 CALL OUTBED(N,BEDH,RM) C FORCE OUT OF BED BUBBLE (THROUGH THE WALL ) BACK INTO THE BED CALL BUBOUT(N,BEDW,RM) C CONTROL OF COALESCENCE CALL COALES(N,X1,RM) CALL BUBOUT(N,BEDW, RM) NK=N CALL MERSON(XI,DELX1,DX1,DXMIN1,TOLKM1,NK,*100,*200,*300) 100 CONTINUE IB=IB+1 IF(IB.LT.5) GOTO 10 IB=0 CALL SKRIV(N,X1) GOTO 10 200 WRITE(2,450) XI 450 FORMAT(IX,'no convergence i n Merson INDEP. VAR.=',F10.5, $ ' v a l u e o f Y=') CALL SKRIV(NK,X1) GOTO 20 300 CONTINUE CALL SKRIV(NK,X1) GOTO 20 10 CONTINUE 40 CONTINUE 20 CONTINUE C C a l c u l a t e b u b b l e d i s t r i b u t i o n C CALL SIDI(XI,RM,BEDH,BEDW) STOP END SUBROUTINE OUTBED(N,BEDH,RM) C PURPOSE: When a bu b b l e l e a v e s t h e bed, t h e r e m a i n i n g bubble C a r e r e o r d e r e d and the number o f e q u a t i o n s i s d e c r e a s e d by 3 DIMENSION Y(1000),DY(1000),R(1000),INTG(6,500) COMMON /DEP/Y/GRAD/DY/RAD/R/INDIC/INTG NP=N DO 10 1=1,N,3 40 IF(Y(I).LT.BEDH/RM) GOTO 10 NP=NP-3 IF(I.GT.NP) GOTO 15 DO 20 J=I,NP K=J+3 Y(J)=Y(K) 20 DY(J)=DY(K) DO 25 J=1,NP,3 K=J+3 R(J)=R(K) DO 30 L=l,6 30 INTG(L,J)=INTG(L, K) 25 CONTINUE GOTO 40 15 CONTINUE 10 CONTINUE N=NP RETURN END 9 0 SUBROUTINE COALES(N,XI,RM) C Purpose: A f t e r c o a l e s c e n c e has o c c u r e d , the bubbles a r e C renumbered, t h e r e c o r d i n g i s c o n t r o l l e d , t h e number o f C e q u t i o n s i s de d u c t e d by 3 f o r each c o a l e s c e n c e . New p o s i t i o n s C and s i z e s o f bu b b l e s a r e c a l c u l a t e d . DIMENSION Y(1000),DY(1000), R(1000),INTG(6,500),ASTRID(6, 1000) $ ,AJAN(6, 1000),AKAREN(6, 1000),CLEV(6),NUM(6) ,XX(6) COMMON /DEP/Y/GRAD/DY/RAD/R/INDIC/INTG COMMON /AST/ASTRID/AJA/AJAN/AKA/AKAREN/CLE/CLEV/NUMB/NUM/XXX/XX NP=N DO 10 1=1,N,3 11=1+1 12=1+2 IF(R(I).LT.0.0001) GOTO 10 DO 20 J=1,N,3 80 I F (R(J).LT.0.0001) GOTO 20 IF(I.EQ.J) GOTO 20 J1=J+1 J2=J+2 C C o a l e s c e n c e c r i t e r i a I F ( A B S ( Y ( I ) + R ( I ) - Y ( J ) - R ( J ) ) . L T . 0 . 5 * ( R ( I ) + R ( J ) ) ) GOTO 90 I F ( A B S ( Y ( I ) - Y ( J ) ) . G T . 0 . 5 * ( R ( I ) + R ( J ) ) ) GOTO 20 90 I F ( ( Y ( I l ) - Y ( J l ) ) * * 2 + ( Y ( I 2 ) - Y ( J 2 ) ) * * 2 . G T . ( R ( I ) + R ( J ) ) * * 2 ) GOTO 20 C RECORD COALESCENCE BUBBLE IF ( ( Y ( J ) + R ( J ) ) . G T . ( Y ( I ) + R ( I ) ) ) GOTO 30 DO 15 K=l,6 IF((Y(I)+R(I)).LT.CLEV(K)/RM) GOTO 26 IF(I N T G ( K , J ) . E Q . l ) GOTO 15 IF(NUM(K).LT.1) XX(K)=X1 NUM(K)=NUM(K)+1 ASTRID(K,NUM(K))=Y(J1) AJAN(K,NUM(K))=Y(J2) AKAREN(K,NUM(K))=R(J) IF(K.NE.4) GOTO 15 WRITE(2,100) Y ( J 1 ) , Y ( J 2 ) , R ( J ) , X 1 100 FORMAT(4F10.4) 15 CONTINUE 30 DO 25 K=l,6 IF(Y{J)+R(J).LT.CLEV(K)/RM) GOTO 26 IF(INTG(K,I).EQ.1) GOTO 25 IF(NUM(K).LT.l) XX(K)=X1 NUM(K)=NUM(K)+1 ASTRID(K,NUM(K))=Y(I1) AJAN(K,NUM(K))=Y(12) AKAREN(K,NUM(K))=R(I) IF(K.NE.4) GOTO 25 WRITE(2,100) Y ( I 1 ) , Y ( I 2 ) , R ( I ) , X 1 25 CONTINUE 26 CONTINUE C New bub b l e s i z e s and p o s i t i o n s RK=R(I) RL=R(J) R(I) = (R(I)**3+R(J)**3)**0.333333*1.05 R(J)=0. IF(Y(I)+RK.GT.Y(J)+RL) GOTO 35 Y(I)=Y(J)+RL-R(I) GOTO 4 5 35 Y(I)=Y(I)+RK-R(I) 91 45 Y(I1)=(RK**3*Y(I1)+RL**3*Y(J1))/(RK**3+RL**3) Y(I2)=(RK**3*Y(I2)+RL**3*Y(J2))/(RK**3+RL**3) C Decrease number o f e q u a t i o n s by 3 NP=NP-3 IF(J.GT.NP) GOTO 55 DO 4 0 M=J,NP K=M+3 Y(M)=Y(K) 40 DY(M)=DY(K) DO 50 M=J,NP,3 K=M+3 R(M)=R(K) DO 60 L=l,6 60 INTG(L,M)=INTG(L, K) 50 CONTINUE 55 CONTINUE N3=NP+1 R(N3)=0. GOTO 80 2 0 CONTINUE 10 CONTINUE N=NP RETURN END SUBROUTINE CLEVS(N,XI,ABC,RM) C Purpose: T h i s s u b r o u t i n e r e c o r d b u b b les which have p a s s e d C g i v e n l e v e l s . INTG keeps t r a c k on bubbbles r e c o r d e d a l r e a d y . C I f INTG(1,4)=0, bubble 4 i s not r e c o r d e d a t l e v e l 1. C I f INTG(1,4)=1, b u b b l e 4 i s a l r e a d y r e c o r d e d a t l e v e l 1. DIMENSION Y(1000), R(1000),ASTRID(6,1000),AKAREN(6,1000), $ AJAN(6,1000),CLEV(6),NUM(6),XX(6),NUMX(6),INTG(6,500) COMMON /DEP/Y/RAD/R/AST/ASTRID/AJA/AJAN/AKA/AKAREN/CLE/CLEV/ $ NUMB/NUM/XXX/XX/NUMBX/NUMX/INDIC/INTG DO 10 1=1,N,3 DO 15 J=l,6 IF(Y(I)+R(I).LT.CLEV(J)/RM) GOTO 10 I F ( I N T G ( J , I ) . E Q . l ) GOTO 15 IF(NUM(J).LT.l) XX(J)=X1 NUM(J)=NUM(J)+1 K=I + 1 L=I+2 ASTRID(J,NUM(J))=Y(K) AJAN(J,NUM(J))=Y(L) AKAREN(J,NUM(J))=Y(L) AKAREN(J,NUM(J))=R(I) INTG(J,I)=1 IF(X1-ABC.LT.0.5) NUMX(J)=NUM(J) C NEXT RUN OF THE UPPER PART OF THE BED C P o s i t i o n s time and s i z e a r e punched when bubbles pass l e v e l 4 C i n the bed, they can be used as boundary c o n d i t i o n s f o r the IF(J.NE.4) GOTO 15 WRITE(2,100) Y ( K ) , Y ( L ) , R ( I ) , X I 100 FORMAT(4F10.4) 15 CONTINUE 10 CONTINUE RETURN END 92 SUBROUTINE BUBOUT(N,BEDW,RM) C Purpose: BUBBLES WHO PASSED THE WALL BY MISTAKE ARE PUT C BACK INTO THE BED DIMENSION Y(1000),R(1000) COMMON /DEP/Y/RAD/R DO 10 1=1,N,3 J=I + 1 K=I+2 IF(Y(J)**2+Y(K)**2.LT.(BEDW/(2.*RM)-1.)**2) GOTO 10 IF(Y(J).GT.O.) GOTO 15 15 Y(J)=Y(J)+0.1*R(I) 16 IF(Y(K).GT.O.) GOTO 17 Y(K)=Y(K)+0.1*R(I) GOTO 10 17 Y(K)=Y(K)-0.1*R(I) 10 CONTINUE RETURN END SUBROUTINE MART(N,T) C PURPOSE CALCULATION OF COEFFICIENT MATRIX A DIMENSION A(1000,1000), Y(1000), R(1000) COMMON /MATR/A/DEP/Y/RAD/R IF(T.GT.0.0001) GOTO 10 DO 5 1=1,150 5 A ( I , I ) = 1 . DO 6 1=1,147,3 J=I+1 K=I+2 A(I,J)=0.0 A(I,K)=0.0 A(J,I)=0.0 A(J,K)=0.0 A(K,I)=0.0 6 A(K,J)=0.0 10 NK=N+1 C L a s t column i n the m a t r i x DO 7 1=1,N,3 J=I+1 K=I+2 A(I,NK)=SQRT(R(I)) A(J,NK)=0.0 7 A(K,NK)=0.0 DO 15 1=1,N,3 11=1+1 12=1+2 DO 20 J=1,N,3 J1=J+1 J2=J+2 IF(I.EQ.J) GOTO 20 C C o n t r o l when bubble I e n t e r e d wake o f bubble J C = ( Y ( J ) - Y ( I ) ) * * 2 + ( Y ( J l ) - Y ( I l ) ) * * 2 + ( Y ( J 2 ) - Y ( I 2 ) ) * * 2 B=(R(J)+R(I))**2 IF(C.GT.B) GOTO 30 I F ( Y ( I ) . G T . Y ( J ) ) GOTO 21 C Set up c o e f f i c e n t m a t r i x when bubble I has e n t e r e d wake o f C bubble J A ( I , J ) = - l . • A(I,J1)=0. 9 3 A ( I , J 2 ) = 0 . A ( I I , J ) = 0 . A ( I 1 , J l ) = - 1 . A (I1,J2)=0. A ( I 2 , J ) = 0 . A ( I 2 , J 1 ) = 0 . A ( I 2 , J 2 ) = - l . DO 25 K=1,N,3 IF(K.EQ.I) GOTO 25 IF(K.EQ.J) GOTO 25 ; • L=K+1 M=K+2 A(I,K)=0. A ( I , L ) = 0 . A(I,M)=0. A(I1,K)=0. A(I1,L)=0. A(I1,M)=0. A(I2,K)=0. A(I2,L)=0. A(I2,M)=0. 25 CONTINUE GOTO 15 21 CONTINUE C s e t up c o e f f i c e n t s m a t r i x when bubble J has p e n e t r a t e d bubble I to C a c e r t a i n degree D=Y(I)-Y(J) E=R(I)+R(J)-SQRT(R(I)*R(J) ) IF(D.GT.E) GOTO 30 A(I,NK)=SQRT((R(I)**3+R(J)**3)**.3333) A ( I , J ) = 0 . A ( I , J l ) = 0 . A ( I , J 2 ) = 0 . A ( I 1 , J ) = 0 . A(I1,J1)=0. A(I1,J2)=0. A ( I 2 , J ) = 0 . A(I2,J1)=0. A(I2,J2)=0. GOTO 20 30 CONTINUE C Set up c o e f f i c e n t m a t r i x f o r the g e n e r a l case, where the bubble C a r e o u t s i d e t h e wake o f each o t h e r X I = Y ( I ) + R ( I ) - Y ( J ) YI=Y(I1)-Y(J1) ZI=Y(I2)-Y(J2) AB= (SQRT(XI**2+YI* * 2+ZI**2))**5* 2 . A ( I , J)=-1.*R(J)**3* ( 2.*XI**2-YI * * 2-ZI * * 2 )/AB A(I,J1)=-3.*R(J)**3*XI*YI/AB A ( I , J 2 ) = - 3 . *R(J)**3*XI*ZI/AB A ( I 1 , J ) = A ( I , J l ) A ( I I , J1)=-1.*R(J)**3*(2.*YI**2-XI**2-ZI**2) /AB A( I I , J 2 ) = - 3 . * R ( J ) * * 3 * Z I * Y I / A B A(12,J)=A(I,J2) A ( I 2 , J 1 ) = A ( I 1 , J 2 ) A ( I 2 , J 2 ) = - 1 . * R ( J ) * * 3 * ( 2 . * Z I * * 2 - X I * * 2 - Y I * * 2 ) / A B 20 CONTINUE 15 CONTINUE RETURN 9 4 END SUBROUTINE SIDI(TR,RM,BH,BW) purp o s e : C a l c u l a t e b u b b l e s i z e d i s t r i b u t i o n and s p a t i a l d i s t r i b u t i o n a t g i v e n l e v e l above t h e boundary c o n d i t i o n DIMENSION ASTRID(6, 1000),AKAREN(6, 1000),AJAN(6, 1000) , CLEV(6), NUM (6),XX(6),NUMX(6),BFQ(6),SD(10),SDC(10),SDA(11),SDP(10), BMV(6),FLOW(6),BMQ(6),FLOWQ(6),BMD(6) DIMENSION SM(6), NUMBT(5), BS(5),BSP(5),UMBT(5) DIMENSION AJAA(6,1000)/ ASTRAA(6,1000) COMMON /AST/ASTRID/AJA/AJAN/AKA/AKAREN/CLE/CLEV/NUMB/ NUM/XXX/XX/NUMBX/NUMX DO 10 1=1,6 NU=NUM(I) DO 5 J=1,NU AKAREN ( I , J) = AKAREN ( I , J) *RM ASTRID(I,J)=ASTRID(I, J)*RM AJAN ( I , J) =AJAN ( I , J) *RM TREAL=(TR-XX(I))/(0.67*SQRT(981.*RM)) C a l c u l a t e bubble f r e q u e n c y , BFQ(bubbles per second) BFQ(I)=NUM(I)/TREAL A r t h e m e t i c mean s i z e , d i a m e t e r SUM=0. DO 20 J = l , NU SUM=SUM+AKAREN(I, J ) * 2 . BMD(I)=SUM/NUM(I) A r t h e m e t i c mean s i z e , volume SUM1=0. DO 30 J=1,NU SUM1=SUM1+1.33333*3.1416*AKAREN(I,J)**3 BMV(I)=SUM1/NUM(I) V i s i b l e gas f l o w FLOW(I)=SUMl/(TREAL*BW**2*3.1416)*4. Frequency PER/cm**2 BMQ(I)=BFQ(I)/BW**2*4./3.1416 Gas flow, v o l u m e t r i c FLOWQ(I)=SUM1/TREAL C a l c u l a t e s i z e d i s t r i b u t i o n based on bubble d i a m e t e r DIT=BW/20.-0.000001 D=-0.000001 J=0 K=l DI=D+0.000001 D=DI+DIT SD(K)=0. DO 50 L=1,NU IF(AKAREN{I, L) .GT.D) GOTO 50 IF (AKAREN(I,L).LT.DI) GOTO 50 SD(K)=SD(K)+AKAREN(I,L)*2. J=J+1 CONTINUE K=K+1 IF(NUM(I).GT.J) GOTO 51 C a l c u l a t e b u bble s i z e d i s t r i b u t i o n p e r c e n t a g e k=k-l kp=k DO 60 L=1,K SDP (L)=SD(L)* 100./SUM CALCULATE CUMULATIVE SIZE DISTRIBUTION CUM=0. DO 70 L=1,K CUM=CUM+SDP(L) 70 SDC(L)=CUM K=K+1 C INTERVALS FOR SIZE DISTRIBUTION DO 80 L=1,K 80 SDA(L)=(L-1)*BW/10. K=K-1 C CALCULATE SPATIAL SIZE DISTRIBUTION, DIMENSIONLESS DISTANCE DO 90 L=1,NU AJAA(I,L)=AJAN(I,L)/BW*2. 90 ASTRAA(I,L)=ASTRID(I,L)/BW*2. C INTERVALS DO 95 L=l,6 95 SM(L)=(L-1)*0.2 C BUBBLE FLOW AND FREQUENCY DISTRIBUTION ON INTERVALS D=0. TAB=0. TUB=0. DO 61 K=l,5 TAL=0. NUMBT(K)=0. D=D+0.2 DI=D-0.2 DO 62 L=1,NU IF(ASTRAA(I,L)**2+AJAA(I,L)**2.GT.D**2) GOTO 62 IF (ASTRAA(I,L)**2+AJAA(I ,L)**2.LT.DI**2) GOTO 62 TAL=TAL+AKAREN(I,L)**3*1.33333*3 .1416 NUMBT(K)=NUMBT(K)+1 62 CONTINUE BS(K)=TAL/((D**2-DI**2)*(BW/2.)**2*TREAL) UMBT(K)=NUMBT(K)/(D**2-DI**2) TAB=TAB+BS(K) TUB=TUB+UMBT(K) 61 CONTINUE DO 63 K=l,5 BSP(K)=BS(K)/TAB*100. 63 UMBT(K)=UMBT(K)/TUB*100. C OUTPUT OF RESULTS WRITE(2,200) CLEV(I),BH,BW,NUM(I) 200 FORMAT('1','Height above t h e d i s t r i b u t o r cm" $ , ' d e v ' , f 8 . 2 / ' T o t a l bed h e i g h t , cm',5x, $ 'bedh=', F8.2/2x,'bed w i d t h cm','cw=', $ F8.2/ 'number o f b u b b l e s ' n * , 1 8 / / / ) WRITE(2,215) 215 FORMAT('0',' c a l c u l a t e d p o s i t i o n s and s i z e s ' ) WRITE(2,230) (ASTRTD(I,J) , J=1,NU) WRITE(2,230)(AJAN(I,J),J=1,NU) WRITE(2,230) (AKAREN(I,J) , J=1,NU) 230 FORMAT(IX,21F6.1) WRITE(2,265) BFQ(I),BMV(I) ,FLOWQ(I) 265 FORMAT('1 bubble f r e q u e n c y , BFQ=',f8.2,4x,'bubble mean s i z e cm** 3 $ BM=',F8.2,'bubble mean flow, cm**3/sec BF=',F8.2,//) WRITE(2,2 66) BMQ(I),BMD(I),FLOW(I) 2 66 FORMAT('bubble f r e q u e n c y , p e r cm**2', f8 . 3 , ' b u b b l e mean diameter, $ cm',F8.3,'bubble mean flow,cm/sec=',f8.3,//) WRITE(2,235) 96 FORMAT(IX,'size o f d i s t r i b u t i o n ' ) w r i t e ( 2 , 2 4 0 ) f o r m a t ( ' s i z e range cm=','bubble diam,cm', $ 'bubble diam p e r c e n t a g e ' , ' c u m u l a t i v e p e r c e n t a g e ' ) DO 97 L=1,KP M=L+1 WRITE(2,245) SDA(L),SDA(M),SD(L),SDP(L),SDC(L) FORMAT(F6.2, '- 1, F6.2, F8.2,F8.2, F8 . 2) WRITE(2,250) FORMAT(//'spatial d i s t r i b u t i o n , h a l f t h e bed',/) WRITE(2,255) FORMAT('position r a n g e ' , ' v i s . b u b b l e f l o w ' , $ 'bubble f l o w p e r c e n t a g e ' , ' b u b b l e f r e q u e n c y p e r c e n t a g e ' DO 96 L=l,5 I1=L+1 WRITE(2, 260) SM(L),SM(I1),BS(L),BSP(L),UMBT(L) FORMAT(F6.2, '- ' , F6 . 2, F8.2, F8.2,F8.2) CONTINUE RETURN END SUBROUTINE BOUNDC(N,RM,BEDW,LOOP) EXTERNAL RN1,RN2 DIMENSION Y(1000), R(1000), DY(1000) COMMON /DEP/Y/RAD/R/GRAD/DY DIMENSION INTG(6,500) COMMON /INDIC/INTG IF(N.EQ.O) RX=5 IF(N.EQ.O) RY=10 DO 10 J=l,LOOP N1=N+1 N2=N+2 N3=N+3 C 30 CALL RN0(RX) C CALL RNO(RY) C JX=JY 30 YFL=RN1(RX) XFL=RN2(RY) Y(N1)=0. Y(N2)=YFL*(BEDW/RM-2.)-BEDW/2./RM+1. Y(N3)=XFL*(BEDW/RM-2.)-BEDW/2./RM+1. I F (Y(N2)**2+Y(N3)**2.GT.(BEDW/2/RM-1)**2) GOTO 30 R(N1)=1. DY(N1)=SQRT(R(N1)) DY(N2)=0. DY(N3)=0. DO 20 L=l,6 20 INTG(L,N1)=0 10 N=N3 RETURN END SUBROUTINE MERSON(X,DELX,DX,DXMIN,TOLKM,N,*,*,*) DIMENSION Y(1000) ,YOLD(1000) ,FK(5, 1000),DY(1000) COMMON /DEP/Y/GRAD/DY TOLA=5.*TOLKM TOLB=TOLA/32. FINTS=DELX/DX+0.5 INTS=IFIX(FINTS) IF(INTS.LT.1) INTS=1 235 240 97 245 250 255 96 260 10 DX=DELX/INTS FMULT=DX/3. GOTO 4 C e r r o r check 1 IF(ERR.GT.TOLA) GOTO 20 IF(ERR.LT.TOLB) GOTO 21 C i n t e r g r a t i o n s a t i s f a c t o r y : c a l c u l a t e new p o i n t s 3 DO 2 1=1, N 2 Y(I)=YOLD(I)+0.5*FK(1,1)+2.0*FK(4,I)+0.5*FK(5,I) IF(INTS.EQ.1)RETURN 1 6 INTS=INTS-1 C p r e s e r v e c u r r e n t v a l u e 4 XOLD=X DO 5 1=1,N 5 YOLD(I)=Y(I) IHALF=0 GOTO 9 C e r r o r e x c e s s i v e : h a l v e s t e p 20 DX=0.5*DX IF(DX.LT.DXMIN) GOTO 19 INTS=INTS+INTS IHALF=1 GOTO 8 C s t e p l e n g t h t o o s m a l l : i n t e r g r a t i o n f a i l s 19 X=XOLD DO 23 1=1,N 23 Y(I)=YOLD(I) RETURN 2 C e r r o r s m a l l : s t e p l e n g t h may be i n c r e a s e d i f p o s s i b l e C check i f s t e p p r e v i o u s l y h a l v e d ( p r e v e n t s c y c l i n g ) 21 IF(IHALF.EQ.1) GOTO 3 C check i f INTS even IDUBLE=INTS/2 . IF(IDUBLE*2.EQ.INTS) GOTO 22 C not p o s s i b l e : INTS odd GOTO 3 C double s t e p l e n g t h 22 INTS=IDUBLE DX=2.*DX C go back t o l a s t p o i n t , i n t e r g r a t e w i t h new DX 8 FMULT=DX/3. DO 7 1=1,N 7 Y(I)=YOLD(I) X=XOLD C main i n t e g r a t i o n p r o c e s s s t a r t s here C advance X by DX 9 CALL DERIVS(X,N,*99) DO 18 IS=1,5 GOTO (31,30,32,33,30), IS 31 X=X+FMULT GOTO 30 32 X=X+FMULT*0.5 GOTO 30 33 X=XOLD+DX C update Y ( I ) 30 DO 10 1=1,N FK(IS,I)=FMULT*DY(I) GOTO(11,12,13,14,10),IS C p r e d i c t o r a t (X+DX/3.) 11 Y(I)=YOLD(I)+FK(l,I) 98 GOTO 10 C c o r r e c t o r f o r (X+DX/3.) 12 Y(I)=YOLD(I)+0.5*(FK(1,I)+FK(2,I)) GOTO 10 C advance t o (X+DX/2.) 13 Y(I)=YOLD(I)+0.375*FK(1,I)+1.125*FK(3,I) GOTO 10 C advance t o (X+DX) 14 Y(I)=YOLD(I)+1.5*FK(l,I)-4.5*FK(3,I)+6.*FKi<i,I) 10 CONTINUE IF(IS.EQ.5) GOTO 16 C e v a l u a t e d e r i v a t i v e CALL DERIVS(X,N,*99) GOTO 18 C on l a s t i n t e g r a t i o n , e v a l u a t e e r r o r 16 ERR=0.0 DO 17 1=1,N EI=ABS(FK(1,1)-4.5*FK(3,1)+4 .0*FK(4 ,1)-0.5*FK(5,1) ) IF(ERR.LT.EI). ERR=EI 17 CONTINUE 18 CONTINUE GOTO 1 9 9 RETURN 3 END SUBROUTINE DERIVS(XL,NL, *) DIMENSION A(1000,1000), DY(1000) COMMON/MATR/A/GRAD/DY INTEGER GUF C nmax i s maximal number o f i n t e r a t i o n NMAX=100 C emax i s maximal e r r o r i n convergence t e s t EMAX=0.05 C c o e f f i c e n t s from s u b r o u t i n e MART CALL MART(NL,XL) NL1=NL+1 C i t e r a t i o n b e g i n s DO 10 NM=1,NMAX GUF=1 DO 11 1=1,NL DYP=DY(I) DY(I)=A(I,NL1) C c a l c u l a t e new v a l u e s DO 12 J=1,NL IF(I.EQ.J) GOTO 12 IF(ABS(A(I,J)).GT.100000.) GOTO 20 IF(ABS(DY(J)).GT.10.) GOTO 21 DY(I)=D Y ( I ) - A ( I , J ) * D Y ( J ) 12 CONTINUE C convergence t e s t IF(ABS(DYP-DY(I) ) .LE.EMAX) GCTO 11 GUF=2 11 CONTINUE IF(GUF.EQ.2) GOTO 10 GOTO 13 10 CONTINUE C i f no convergence w r i t e c u r r e n t v a l u e s o f DY and stop WRITE'(2, 300)NM RETURN 1 20 WRITE(2, 400)1, J \ 99 RETURN 1 21 WRITE(2,500)J RETURN 1 13 CONTINUE 300 FORMAT('no convergence i n G a u s s - S e i d e l , number o f i t e r a t i o n s = ' , 1 5 $ , ' c u r r e n t v a l u e s o f Y,R and DY 1) 400 FORMAT(' A ( I , J ) i s t o o h i g h , I,J=',2I5) 500 FORMAT(' DY(J) i s too h i g h , J=',I5) RETURN END SUBROUTINE SKRIV(N,X1) DIMENSION Y(1000),DY(1000) ,R(1000) COMMON /DEP/Y/GRAD/DY/RAD/R C WRITE(2,400) XI 400 FORMAT('value o f i n d e p . v a r . T=' , F10.5,//) WRITE(2,410) (Y(K),K=1,N) 410 FORMAT(4X,3F10.5) WRITE(4,411) (R(K),K=1,N,3) 411 FORMAT(3X, F10.5) C WRITE(2,412) (DY(K),K=l,N) 412 FORMAT)'DY(I) 1, 3F10.5,3X,3F10.5) WRITE(5,413) N 413 FORMAT(IX, 13) RETURN END FUNCTION RN1(R) S=65536. U=2053. V=13849. M=R/S R=R-M*S R=U*R+V M=R/S R=R-M*S RN1=R/S RETURN END FUNCTION RN2(R) S=65536. U=2053. V=13849. M=R/S R=R-M*S R=U*R+V M=R/S R=R-M* S RN2=R/S RETURN END i 100 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items