Multiresolution analysis of gas fluidization by Empirical Mode Decomposition and Recurrence Quantification Analysis.

The dynamics of gas fluidization is investigated by means of multiresolution analysis. Empirical mode decomposition (EMD) and the Hilbert-Huang transform approach applied to the signals of pressure fluctuation of the bed have been used for this purpose, operating in bubbling and slugging regimes. To elucidate the different components of the different scales, recurrence plots (RP) and recurrence quantification analysis (RQA) have been used. These techniques can distinguish the three different scales of gas fluidization: micro-scale, meso-scale and macro-scale, and classify every mode on its scale. Three modes from the EMD have been related to each dynamic component: particle interaction, local bubble dynamics and bed oscillation, showing evidence of this relationship. To show that the complexity of the modes matches with their characteristics, two measures have been computed: the apparent entropy and Lemped-Ziv complexity.


Introduction
Although gas fluidization has been industrially used for over a century, it remains a complex technique that still attracts researchers aiming to improve control and performance operation.It is recommended when good gas-solid contact is needed, improving mass and heat transfer.Also, its peculiar "fluid-like" dynamics greatly facilitates the handling and processing of solids in industrial processes.It is well known that particulate fluidization is desired for an optimal contact between the phases, but aggregative fluidization is quite common in industrial applications.Furthermore, for deep beds the bubbling fluidization can evolve to slugging fluidization, which is not convenient because then part of the gas bypasses the solid contact.
Gas fluidization is a complex dynamic system characterized by non-linearity and nonequilibrium.The complexity of its dynamics is due to the interaction of the phases involved at different scales in the heterogeneous flow structure.Complex systems are characterized by a multi-scale structure nature with respect to space and time, showing dissipative structures by non-linear and non-equilibrium interactions and exchanging energy, matter and information with their surroundings (Li, 2000;Li and Kwauk, 2003;Li et al., 2004).Also, gas fluidization is a typical dissipative structure consisting of a non-equilibrium system with particle and fluid self-organization.A considerable amount of the total input energy is dissipated to maintain the two phase heterogeneous structure.Nevertheless, all the phenomena that take place in gas fluidization are the result of the nonlinear interaction between the particles and the fluid with their own individual movement tendencies.The dissipative structure in fluidized systems has been found to show multi-scale characteristics and the sum of individual processes does not properly reflect the dynamics of the system.Therefore, different scales must be considered for a detailed analysis.The system can be structured into three basic scales: micro-scale (individual particle and fluid scale), meso-scale (cluster and dilute phase scale, or "bubble and emulsion") and macro-scale (effect of the equipment) (Li, 2000).This approach through different structures involving the fluidized bed is crucial to better understand the behavior of the bed dynamics and the influence of the different structures in the bed.
Among the techniques used to characterize the fluidization, the pressure fluctuations analysis is perhaps the most popular one because it is easy to implement and inexpensive, especially in industrial installations.Other techniques allow for more and better information, but their implementation is much more complex.In spite of the limited information provided by the pressure fluctuations of the fluidized bed, when properly threated, it may be a consistent and valuable source of information.It has even been considered to be the fingerprint of the system.The pressure fluctuations have been analyzed by means of time and frequency domain techniques (Fan et al., 1991;Johnsson et al., 2000;Llop and Jand, 2003;Sasic et al., 2007).However, given the complexity of multiphase flow and its nonlinear behavior, the analysis by nonlinear techniques has been introduced (van den Bleek and Schouten, 1993;Zijervelt et al., 1998, Johnsson et al., 2000, Llauró and Llop, 2006, Llop et al., 2012).
Several researchers have applied wavelet analysis to the experimental time series of pressure fluctuations for multiscale resolution (Lu and Li, 1999;Zhao andYang, 2003, Shou andLeu, 2005;Wu et al., 2007, Tahmasebpour, 2015).To a lesser extent, Empirical Mode Decomposition (EMD) of the signal has been used to extract intrinsic mode functions (IMF) with different frequencies, which can be related to the different scales.The signal analysis approach has attracted the attention in several research fields.Briongos et al. (2006) used the multiresolution analysis EMD to study the hydrodynamics of a gas-solid fluidized bed.
Wavelets can handle non-stationary signals due to the nature of wavelet functions and, although they are basically suited for linear signals, they have been used successfully in non-linear systems.EMD is suitable for nonlinear and non-stationary systems and makes it possible to simultaneously obtain the real time and the instantaneous frequency and can classify time or frequency dependent information with more accuracy.
In this work the multiscale resolution of bubbling and slugging regimes has been studied, decomposing the pressure fluctuations in the bed by EMD while, using the Hilbert-Huang transform method, the intrinsic frequency has been extracted.With the aim to analyze the different modes obtained with EMD, Recurrence Plots (RP) and Recurrence Quantification Analysis (RQA) have been used.The modes extracted have been related with the particle, bubbles and bulk structures in the fluidized bed.To analyze the structure and to be able to discuss the behavior of the modes generated with the EMD two measures of the complexity have been used: approximate entropy (ApEn) and Lempel-Ziv (L-Z) complexity, which are very useful parameters to characterize spatiotemporal patterns.

Empirical Mode Decomposition (EMD)
Wavelet analysis has been used for the decomposition of the pressure fluctuation signal in different levels of resolution, and related to the three scales associated to the fluidization dynamics: micro-scale, meso-scale and macro-scale.Zhao and Yang (2003) studied the fractal behavior of resulting levels from Daubechies wavelets with the Hurst exponent.Lu and Li (1999) used the discrete wavelet transform to analyze the pressure signal in a bubbling fluidized bed and related the bubble size with the average peak value of level 4 and the bubble frequency with the peak frequency of this level.Tahmasebpour et al. (2015)  The algorithm of extraction proposed by Huang (Huang et al. 1998) generates upper and lower smooth envelopes enclosing the signal.These envelopes are generated by the identification of all local extrema, which are connected by cubic spline lines.A new function is obtained by subtracting the running mean of the envelope from the original data signal.If this function has the same number of zero-crossing points and extrema, the first IMF is obtained, which contains the highest frequency oscillations in the signal.
Otherwise, the process must continue until an acceptable tolerance is reached.To extract the following IMF, the previous IMF is subtracted from the original signal.The difference will be treated like the original data and the process is applied again until the above mentioned condition is fulfilled.The process of finding the several modes is carried out until the last mode (the residue) is found.The original signal is the sum of the different modes generated, Where C i is every mode extracted and r n is the residual part of the signal.
Once the modes have been extracted, a second process must be done.The instantaneous frequency is computed by applying the Hilbert transform to every mode so that the time-frequency distribution of the signal energy is obtained.Each mode function C i (t) is associated with its Hilbert Spectral Analysis H i (t): and the combination of C i (t) and H i (t) gives the analytical signal which can be expressed as: is the amplitude of the signal and  i (t) is the phase of the oscillation mode "i".Hence, the original time series, neglecting the residual part, can be expressed as: Re meaning the real part.The amplitude   t A i and the phase  i (t) times series can be computed by: The instantaneous frequency (f i ) can be obtained by differentiating the phase angle: For each mode, the Hilbert spectrum can be defined as the square amplitude The spectrum provides an intuitive visualization of the instantaneous frequencies of the signal in the time scale, showing where the energy is concentrated in time and frequency space.This method has been proved as an efficient procedure to distinguish the signal trend from its small scale fluctuations.

Recurrence Quantification Analysis
In the phase space the time evolution of the variables of a dynamic system can be represented with a trajectory.If the trajectory is attracted to a region of the space, this is called an attractor.The trajectory can be derived from the differential equations which define the system but, when the variable under study comes from experimental data included in a time series, the attractor can be reconstructed from the experimental data and copies delayed in time as postulated in the "embedding theorem" (Takens, 1981).Several invariant parameters can characterize the attractor, probably the most popular are: the Kolmogorov entropy, the correlation dimension and the largest Lyapunov exponent.But the attractor can also be characterized by means of Recurrence Plots (RP) and Recurrence Quantification Analisys (RQA) (Marwan et al., 2002, Tahmasebpour et al., 2015;Llop et al., 2015).Many dynamic systems, especially the nonlinear systems, exhibit recurrence behavior.The recurrence of states takes place in the system phase space and can be the source of relevant information (Marwan et al. 2007).The RP is the tool that measures and visualizes recurrences of a trajectory in phase space.However, for attractor reconstruction the appropriate embedding dimension 'm' and time delay '' must be chosen, according to the embedding theorem.
With Recurrence Plots (RP's) it is possible to analyze the time series without having to care too much about these parameters, since it can be considered that both take the value of one (March et al., 2005 andZbilut andWebber, 2006).Moreover, if the dimension is higher than three, obviously, the state space trajectory is very difficult to visualize.In contrast, the phase state trajectory can be visualized in a two dimension plot by the Recurrence Plot because this one is independent of the trajectory dimension (Eckmann et al. 1987).The recurrence occurs when a point of the trajectory repeats itself.The term repetition means that a point of the trajectory is close enough to another one within a distance suitably selected in advance.
The RP of a signal is generated from the reconstructed attractor computing the matrix, N is the number of states in the space state, i y  and are two different points of the space trajectory,  is a threshold distance of neighborhood, • represents the norm and  is the Heaviside function     . Once generated the N N  matrix with zeros and ones, the two dimensional graphical representation of R ij can be plotted by assigning a white dot to the value zero and a black dot to the value one.The graphical representation is named Recurrence Plot (RP).The qualitative patterns of RP give information about some characteristics like homogeneity, periodicity, drifting, disruption, etc. (Marwan et al. 2002).Furthermore, the RP's can give information about the influence of the micro-scale and macro-scale in the fluidization dynamics (Babaei et al. 2012, Llop et al. 2015).Although the plots yield very useful information, there must be a quantifiable criterion to detect a transition in the system dynamics.Zbilut and Webber (1992) and Webber and Zbilut (1994) introduced a methodology to quantify the plots named Recurrence Quantification Analysis (RQA) extended later by Marwan et al. (2002).RQA is a set of parameters conceived to describe the complex structure of the RP.These parameters quantify the complexity of the system and their computation is based on recurrence point density and on diagonal and vertical lines in the structure of the plot.Although the embedding dimension and the delay are not critical parameters to reconstruct the attractor, the suitable choice of a threshold distance () is crucial.For too small values of , no recurrences are identified and no information is extracted.
Conversely, if  is too large the consecutive points of the trajectory will be considered a recurrence.Several criteria have been proposed for the threshold distance choice (Marwan et al. 2007).One of the simplest criteria is to choose  in the scaling region of the recurrence point density parameter (RR) of RQA (Webber and Zbilut, 2005).
Nevertheless, the choice of  strongly depends on the considered system under study (Marwan et al. 2007).This has been the approach used in this contribution.
The RQA parameters used in this work have been: the recurrence rate (RR), the determinism (DET), the average length of diagonal lines (L), the entropy (ENT), the laminarity (LAM), the trapping time (TT), the recurrence time of first type (RET1) and the recurrence time of second type (RET2).A description of these parameters is skipped here but can be found in Marwan et al. (2007) or, for a similar case, in Llop et al. (2015).

Complexity of Lempel and Ziv
The complexity can be measured by different parameters.Some of them come from chaos theory such as Kolmogorov entropy, Lyapunov exponents or fractal dimension.
Nevertheless, these parameters are difficult to estimate and require a long computation time and a large amount of data.Lempel and Ziv (1976) introduced a measure of the complexity which has proved an excellent tool to characterize spatiotemporal patterns (Kaspar and Shuster, 1987).The main feature of the Lempel-Ziv complexity (LZ) is that it contains the notion of complexity in the deterministic sense as well as in a statistical sense and can be computed more easily with little computation time while being unaffected by the external noise.
The LZ complexity is analyzed by transforming the signal into a sequence P having a finite number of symbols.This sequence is monitored and when a new subsequence of consecutive characters is found, a counter of complexity increases in a unit thus obtaining c(n).Essentially the algorithm is as follows: Being the sequence  

If
one.The next step is to update the sequences S and Q so that . All these steps are repeated until the last character r = n is reached.In a last step, a normalization of the counter c(n) is applied dividing it by the factor . The objective of this normalization is to obtain a value of complexity independent of the length of the time series.The time series has been transformed into a symbolic sequence composed of zeros and ones using four symbols, which is enough to reconstruct the time series (Wang et al. 2010).

Approximate entropy
Pincus developed the approximate entropy (ApEn) as a measure of complexity for time series of experimental data, applicable to noisy and medium-sized datasets by modifying an exact regularity statistic, the Kolmogorov-Sinai entropy (Pincus, 1991).
ApEn measures the logarithmic probability that the series of patterns that are close remain close in the following incremental comparisons.ApEn is a value that reflects the predictability of future values in a time series based on previous values, with larger values corresponding to more complexity in the data.
Being the time series of pressure fluctuations with dimension N m is the length of window to be compared.For each "i" and being "r" the tolerance for accepting matches, The Approximate Entropy, then, is where  defines the criterion of similarity between patterns that represents the noise filter level.The larger the value of ApEn, the higher the complexity of the tested time series.

Experimental
The experimentation was carried out in two different cylindrical columns of diameters (D c ) 30 mm and 48 mm, made of PMM to see inside.Bronze particles of two sizes were used (99 m and 211 m), with a density of 8770 kg/m 3 and the bed heights (H s ) were from 80 mm to 108 mm of settled particles.A 2 mm thick perforated plate was used for the distributor with a nylon mesh placed above to avoid the weeping of particles when the flow was cut.Raschig rings were placed in the calming section below the distributor, in order to increase the uniformity of the gas distribution.The fluidization gas was air at ambient pressure and room temperature, driven by a compressor.The fluctuations due to the compressor were kept to a minimum since the air was driven through an intermediate tank accumulating the sufficient amount of air at 6 bar with a narrow pressure regulation.The fluidizing gas passed through a pressure reducing valve to place it at working pressure.The inlet air flow rate was measured under ambient condition using a rotameter manifold and controlled by a needle valve.
The differential pressure fluctuations in the bed were measured by a probe of 3 mm inner diameter, vertically inserted within the bed, close to the distributor.The probe was connected to the first connection of a piezoresistive differential pressure transducer with a response time of 5 ms (200 Hz).The second connection of the transducer was open to the freeboard.The analogic signal generated by the transducer was digitalized in a 12 bits A/D converter and stored in a PC with the adequate software to be analyzed later.For every run (every gas velocity selected) a time series of 8192 data with 100 Hz of sampling frequency was obtained.For the chaotic analysis of time series large amounts of data are recommended.Johnsson et al. (2000) used 65000 data, but other authors have obtained reliable results with less data; Llop et al. 2012 used 8192 data.On the other hand, the data length is not critical for the analysis with RP and the length of the time series used in this work is sufficient (Tahmasebpour et al., 2015;Llop et al., 2015).For the frequency sampling (Llop et al. 2015) in RP analysis 100 Hz are also enough.The velocity of incipient fluidization was obtained plotting the pressure loss through the bed versus the gas velocity while gradually reducing it from its maximum.

Results and discussion
The EMD has been computed for the time series of pressure fluctuations obtained experimentally in the fluidized bed for the two sizes of particles used and for bubbling and slugging regimes at several velocities.Fig. 1 shows an example of the decomposition of the original signal in several IMF and the residual component.For every IMF the percentage of energy contribution and the mean instantaneous frequency have been computed.Also, the RP has been generated for every IMF and the corresponding RQA has been calculated.With this information it has been possible to elucidate which IMF are included in the different scales of the fluidization dynamics in bubbling and slugging regimes.

Energy and frequency distribution
The phenomenon that predominates in the micro-scale is the interaction between the particles and between the particles and the bed wall.It is qualitatively obvious that high frequency and low energy are associated to it.The meso-scale reflects the bubble and  the bubble-particle interaction (the main cause of pressure fluctuations).Hence, the frequency decreases while the energy significantly increases.The macro-scale mainly reflects the equipment contribution and is characterized by low frequency and low energy.
The energy distribution for every IMF obtained with the decomposition is shown in Fig. 2. For bubbling regime the energy increases quickly from the lower mode to the third, which has the maximum energy, and afterwards quickly decreases to very low values.
As the higher contribution to the fluidization energy is due to the bubble dynamics, this third mode is, most probably, strongly related to the gas phase.Regarding to the frequency, it is very high for the first mode and quickly decreases with the mode number.In the slugging regime (Fig. 2 (c)), the energy is very low for the first three modes and near constant.It is larger from the fourth to the eight modes with a maximum at the fifth mode.This maximum happens at lower frequencies than in the bubbling regime.It is quite likely that the maximum is due to the bubble dynamics.Slugging produces large bubbles, or a train of large bubbles that rise through the bed, explode on the surface and impel the solid particles upstream, which afterwards return to the bed.The slugging dynamics causes a strong periodic signal in the pressure fluctuations with a lower frequency than in the bubbling regime, because the probe detects fewer bubbles in the same period of time.For the bubbling regime the frequency of maximum energy is 3.2 Hz (Fig. 2 (b)) while for slugging regime it is 1.7 Hz (Fig. 2 (d)).The bubble frequency for slugging regime is lower because the slugs fill the column diameter and hence, for the same flow of gas, the frequency of bubble passage is lower.A similar behavior has been verified for particle diameter of 96 m (Fig. 3  Comparing the micro-scale in the bubbling and the slugging regimes it must be observed that for slugging regime LBA is bigger.It is so because in the slugging regime there is an alternating flow between a slug and a dense phase portion.As the particles are confined in this phase, the interaction is more intense as shown in the RP plot, the In the meso-scale there is no big difference in LWA or LBA between the two regimes.
Nevertheless, it has been noticed that its distribution is less uniform in the case of slugging regime with zones including more concentration of LBA.This can also be attributed to the alternation of slugs and dense phase.

Recurrence quantification analysis
4.3.1.Bubbling regime Fig. 5 shows the evolution of the RQA parameters with the levels generated by EMD decomposition of the pressure fluctuations in the bubbling regime.As can be seen in Fig. 5(a) the RR evolution is not affected when changing the gas velocity.Also, the evolution of RR is similar for the two threshold values tested.
On the other hand, it is worth noting that the RR decreases from level C 1 to C 2 and increases from C 3 to C 10 .The increase is done with two different slopes.The first of them is rather high except for the relative gas velocity u/u mf = 3.35.The second slope, slightly lower, includes the levels from C 7 to C 10 .This evolution suggests different structures in these level intervals.
As stated above, micro-scale bed dynamics is originated mainly by particle interaction and particle-wall interaction.This movement is done with a certain periodicity so that the interactions take place with high frequency.As seen in   The RR value of C 2 is lower, near to the maximum of energy percentage (C 3 ).
Considering also Fig. 5(a), the RR of C 2 considerably differs from that of C 1 , but in decreasing order, approaching to the value of the C 3 .The combination of high frequency and low energy suggests that C 1 and C 2 are related to the micro-scale.
Meso-scale is related to the bubble phase and the interaction of bubbles with particles, which are the main contributor to the fluidization energy.For this reason, it can be suggested that the meso-scale is located in the zone of medium frequency and high energy (Fig. 2 and Fig. 3).Considering Fig. 5 (a), the RR from C 3 to C 6 follows the same trend, but from C 7 to C 10 the trend slightly changes, decreasing the slope.The level interval from C 7 to C 10 has both very low frequency and energy, which points out to equipment interaction.The change of structure of this level is backed by the turnaround of the RR value (Fig. 5 (a)).
The DET has been plotted in Fig. 5 (b) and shows a remarkable different evolution depending on the threshold radiuses.For the higher value  =0.05, the DET is saturated (constant) from C 4 to the end, while for •=0.005 it is saturated from C 7 onwards, where the trend change appears.This behavior is similar to the evolution observed in LAM (Fig. 5 (e)).This suggests that the bigger value of  is not appropriate to properly discriminate the different structures of the levels.On the other hand, as shown in Figs. 5 (c), (d) and (f), the behavior of the average length of diagonal lines (L), entropy (ENT) and trapping time (TT) are very similar.The fact that there is a common or similar behavior with the lower  for all the RQA investigated parameters induces to consider that this value of radius threshold is the most appropriate.
The RR, DET and LAM evolution shows a similar behavior in the region of micro-scale and meso-scale, and they can be used as indicators of the transition between these regions.Nevertheless, the DET and LAM values become saturated when IMF is C 7 and superior.This happens because the macro-scale (C 7 to C 10 ) has a strong and highly deterministic periodic component, which gives a maximum value of one.This observation is similar to those of Tahmasebpour et al. (2015), which found the maximum value of determinism for the higher scales, although in their case the signal was decomposed by wavelets.As can be seen in Fig. 5, all RQA parameters can group each level on its scale when =0.005 is used.

Slugging regime
The evolution of the RQA parameters for the slugging regime has been plotted in Fig. 6.For consistency, the chosen threshold radius has been the same as for the bubbling regime ( =0.005).As can be seen, all parameters can characterize the different regions corresponding to the three studied scales.
The differentiation with L and TT parameters is done with excellent precision (Figs.6(c) and 6 (f)).The other parameters can separate the three regions, although with less precision, in the intersection or confluence of two different scales, especially to clearly identify the border between micro-scale and meso-scale.However, the increase in energy that can be seen in Figs.2(c) and 3(c) confirms that there is a change in the structure of level.Assuming that the maximum energy corresponds to the bubble contribution, this region is constituted of meso-scale.The levels C 1 , C 2 and C 3 have very low energy.In the next levels there is an increase of the energy to a maximum in C 11 to the macro-scale.Thus, the classification is reasonably suitable despite some uncertainty at the scales confluence.
The RQA parameters evolution for slugging regime is similar to that for bubbling regime except for a significant difference: the maximum fluidization energy is displaced to higher level scales: from C 3 for bubbling regime to C 5 for slugging regime, which means that the bubble phase decreases its frequency as shown in Figs.2(d) and 3(d).
For both bubbling and slugging regimes the determinism (DET) decreases in the microscale, increases in the meso-scale and remains constant in the macro-scale.The determinism is related to the predictability of the system and has low values for stochastic systems and high values for predictable systems.
As shown in Fig. 1 (a) for bubbling regime, the level is much more periodic in the macro-scale (C 8 ) than in the meso-scale (C 4 ) i.e. more deterministic and, because of that, the determinism is 1 for macro-scale.All in all, comparing the DET of the detail levels (IMF) for both regimes, the values of the slugging regime are higher because this regime has a higher periodic component.
The entropy measures the complexity of the system.One would expect the ENT to decrease for macro-scale, the more periodic system.But the ENT increases with the periodicity, an apparent contradiction.This effect was reported by Webber and Zbilut (2005) who argued that the truncation of the diagonal lines in the RP influences the computation.As will be seen below, other measures of complexity (ApEn and LZ) will give results consistent with this argumentation.Even if the evolution of entropy is not logical, its increase in this region is also an indicator to distinguish the scales.

Intrinsic mode function analysis
The pressure fluctuations in the fluidized bed are the result of the contribution of different causes.The bubbles passing near the pressure measurement probe cause a variable pressure field and, consequently, pressure fluctuations.The pressure increases when the bubbles come near to the probe and decreases when they move away.Furthermore, the bubbles eruption on the bed surface originates surface waves; the sloshing motion of the surface and the inertial forces associated to the dense phase cause the bed oscillation with its natural frequency.The pressure fluctuations in the

Mode
bulk of the bed are originated by this natural bed oscillation (Baskakov et al.,1986).The contribution to the pressure fluctuations due to the particle dynamics (particle-particle interaction and particle-wall interaction), despite its low energy, must also be taken into consideration (Briongos et al., 2006).Other minor causes of pressure fluctuations are bubble formation, bubble coalescence and bubble disintegration.
Thus, the pressure fluctuations on the bed are the result of these individual contributions associated with the local and the global dynamics.Some of these oscillations would have different frequencies because of its different nature.EMD generates several IMF with different mean frequencies.It is likely, then, that some IMF can be associated to an individual cause of the pressure fluctuations.This makes reasonable to investigate the possible relationship of the modes obtained by EMD from the pressure fluctuations with the particle dynamics, the bulk dynamics and the bubble dynamics.Briongos et al. (2006) considered three detail levels of interest for engineering purposes in the analysis of the slugging fluidization dynamics.The first of them is the particle dynamics component: particle-particle, particle-wall interactions and coherent structures.The second level, or bulk component, is due to the natural oscillation of the bed.The bed moves up and down due to the gas circulation and to the elastic movement of the particles in the emulsion phase.The third level corresponds to the bubble dynamics, which has great importance in the overall dynamics of bubbling and slugging fluidized beds.A very significant amount of the dissipated energy in the fluidization corresponds to the bubble phase.Hence, the behavior of the bubble phase is of key importance for the design of fluidized beds.This is the reason for countless studies on this subject for decades (Cranfield and Geldart, 1974;Mori and Wen, 1975;Darton et al., 1977;Hilligardt and Werther, 1986;Choi et al., 1988;Choi et al., 1998).
To prove the connection of the IMF with these pressure fluctuations contribution, some models that describe their behavior has been used.The model of Baskakov et al. (Baskakov et al., 1986) has been used to predict the natural frequency of the bed (f bs ), with the following equation, The model of Choi et al. (Choi et al., 1988;1998) has been used to estimate the bubble frequency in the bed for the bubbling regime.With this model it is possible to predict the bubble frequency from the two phase theory and the bubble growth theory.The bubble frequency (f b ) can be estimated by the following expression: Where N f is the number-based bubble flow rate (s -1 ), d bf the frontal bubble diameter and A the bed transversal area.N f can be calculated by, being d b the bubble diameter.The bubble growth used in this theory is, For slugging regime the model of bubble frequency (f s ) proposed by De Luca et al.

 
It can be seen in Fig. 2 that the maximum energy corresponds to level C 3 for bubbling regime and to C 5 for slugging regime.It is coherent to assume that the maximum of the fluidization energy is associated to the local bubble dynamics.Briongos et al. (2006) justified that the level C 5 , obtained by EMD, was related to the local bubble dynamics for slugging regime.
The mean frequency evolution of every level versus the excess of gas velocity from incipient fluidization (u-u mf ) for the two particle diameters and regimes investigated is shown in Figs.7 and 8.It can be seen that the frequency values for every level do not change significantly with the gas velocity, but this is due to the narrow range of values investigated.
The plots for the bubbling regime (Figs.As can be seen in Fig. 9 (a) and (c) the C 1 component has more complexity (higher value of ApEn) than the original signal for all the relative velocities that have been studied.At first, the entropy does not change so much with the velocity of the gas but, for higher velocities, it considerably grows.This evolution is due to the increase in complexity originated by the increase of the inter-particle and wall-particle collisions at higher velocities.
The level C 3 and the original signal have similar complexity.But the original signal has a slightly higher complexity than the level C 3 surely due to the influence of the disordered particle motion on the bed, since it brings together all the involved phenomena, including particle interaction.
Figs. 9 (b) and (d) show the evolution of the ApEn for both particle diameters investigated and for slugging regime.The plot for d p =211 m shows that the component C 1 increases its complexity at higher velocity due to the intense particle motion.The level C 5 is almost constant with a trend to slightly decrease when the gas velocity grows, which is more pronounced for particles of 96 m of diameter.A decrease in the entropy would be justified by an increase in slug size when the gas velocity grows.The slug occupies more space in the column and there is less possibility of mixing with the dense phase.However, the entropies are lower in the case of the slugging regime for all component levels, since the dynamics for this regime is much more ordered (lower complexity).For this regime the level C 5 related to the local bubble component also has lower entropy than the original signal because the original signal includes the influence of the disordered movement of the particles when the slugs burst.On the other hand, the level C 4 related to the bulk component has an ApEn something higher than the original signal and also than the level component C 5 .This happens because this level only reflects the bulk dynamics of the bed and hence it is more complex than the dynamics of the bubbles (C 5 ) or than the combination of all phenomena including the particleparticle interaction (original signal).
A similar behavior can be seen in Fig. 10 where the LZ complexity has been plotted versus the difference of gas velocity.This parameter shows the clearest trends, especially for the slugging regime.The level component C 1 , as in the case of approximate entropy, shows the highest entropy while C 3 , the lowest.The original signal has an intermediate complexity between those levels.This is obvious due to the motion of the particles and the lesser complexity of the bubble motion component (C 3 ).
As the original signal includes both phenomena, it reflects the behavior of both contributions.For 211 m particles, de C 4 level (related to the bulk component) has an LZ complexity slightly higher than the original signal and the C 5 level.This can be ascribed to the increase in particle motion, something also indicated by ApEn.In contrast, for 96 m particles the C 4 level has slightly lower complexity than the original, this may be due to the lower energy of particle interaction.But in any case they have a very similar value.
The hypothesis that the detail level C 1 from EMD decomposition is related to particle dynamics is coherent analyzing the mutual information of these level scales.The that is it is more dissipative.

Conclusions
The structure of the gas solid fluidized bed has been analyzed by means of the multiresolution approach using the empirical mode decomposition (EMD) for bubbling and slugging regimes.It has been shown that the Recurrence Quantification Analysis (RQA) is an excellent tool to differentiate the scales into which the structure of gas fluidization can be classified.The parameters computed by RQA are able to distinguish three different zones which correspond to the micro-scale, meso-scale and macroscale.
The analysis of IMF extracted from the pressure fluctuations of the bed by EMD reflects the characteristics of the bubbling and the slugging regimes.The average frequency of the more relevant intrinsic mode functions related to local bubble component is lower for slugging regime than for bubbling regime.For the bubble regime the range average frequency is from 3.5 Hz to 5.5 Hz and for slugging regime from 1.5 Hz to 2.5 Hz.
Analyzing both the energy distribution and the RQA parameters of every intrinsic mode function computed (IMF), it can be concluded that C 1 and C 2 belong to micro-scale, C 3 to C 6 are included in meso-scale and from C 7 above belong to macro-scale for bubbling regime, while for slugging regime C 1 to C 3 belong to micro-scale, C 4 to C 8 are included in micro-scale and from C 9 above belong to macro-scale.
Comparing the evolution of the Recurrence plots (RP) the LWA pattern increases in the sense of increasing the level of the scale.That is to say, the level related to the local bubble dynamics (belonging to meso-scale) has considerably more LWA than the level related to particle interaction (belonging to micro-scale).Also, the level representative of macro-scale has more LWA than the level related to the local bubble dynamics.
Conversely, the LBA has an inverse behavior.
It has been also observed that the levels with more energy of fluidization are reasonably related to local bubble dynamics.Nevertheless, the level C 3 which has the highest energy value for bubbling regime, is probably the combination of two causes: the local bubble dynamics and the bed oscillation.Moreover, the level C 4 is related to bed oscillation and C 5 to the local bubble dynamics for slugging regime.
The behavior of the C 1 level shows that it is related to the particle-particle interaction and particle-wall interaction for both the bubbling and slugging regime.
Both the approximate entropy and the Lempel-Ziv complexity values are coherent with the evidence of the relation of the levels obtained by EMD with the dynamics of the main phenomena associated to gas-solid fluidization.
decomposed the pressure signals by means of Daubechies wavelets in three sub-signals representing the three scales of fluidized bed dynamics and analyzed them by means of Recurrence Plots and Recurrence Quantification Analysis.Both the Fast Fourier Transform (FFT) and the Wavelet Transform (WT) can analyze nonstationary signals but have limited accuracy when used to classify time or frequency dependent information.The Hilbert-Huang Transform (HHT) is a timefrequency analysis method developed for the analysis of non-stationary and non-linear time series introduced by Huang et al. (1998), particularly suited for nonlinear processes.The result is a combination of an empirical approach with a theoretical tool, which has been successfully used in several fields of research like meteorology, seismology, multiphase flow, etc.The HHT is based on the Empirical Mode Decomposition (EMD), which decomposes the signal in several oscillatory modes, named Intrinsic Mode Functions (IMF).The EMD is based on the sequential extraction of energy associated with the intrinsic mode functions of the signal, from finer temporal scales (high frequency modes) to coarser ones (low frequency modes).

.
SQ is the SQ sequence removing the last character.Being, Being r the index of the symbol analyzed.

Fig. 1 .
Fig. 1.EMD decomposition of pressure fluctuations for d p =211 m particles at u/u mf =3.35(a) for bubbling regime and (b) for slugging regime.

Fig. 2 .
Fig. 2. Energy and mean frequency distribution for the different modes decomposed for d p =211•m particles (a) and (b) for bubbling regime (Dc=48mm and Hs=108 mm), and (c) and (d) for slugging regime (Dc=30 mm and Hs=80mm).

Fig. 3 .
Fig. 2. (continued) more important.The interaction of particles in bubbling regime is less intense because their movement extends all along the bed.
Fig 2(b), the levels C 1 and C 2 have the biggest frequencies, suggesting that both belong to the micro-scale.On the other hand, C 1 shows low percentage of energy, which may agree with the collision between particles.

C 5 ,Fig. 5 .
Fig. 5. Evolution of the parameters of RQA for the different IMF computed for bubbling regime and particles of d p =211 m.

Fig. 6 .
Fig. 6.RQA evolution for the different modes computed for slugging regime, particles of