slugging flow behavior in a cylindrical fluidized bed: ECT measurements and two-fluid simulations

Journal Pre-proofs Bubbling/Slugging Flow Behavior in a Cylindrical Fluidized Bed: ECT Measurements and Two-Fluid Simulations Brajesh K. Singh, Shantanu Roy, Vivek V. Buwa PII: DOI: Reference:

S1385-8947(19)32532-X https://doi.org/10.1016/j.cej.2019.123120 CEJ 123120

To appear in:

Chemical Engineering Journal

Received Date: Revised Date: Accepted Date:

22 January 2019 1 October 2019 9 October 2019

Please cite this article as: B.K. Singh, S. Roy, V.V. Buwa, Bubbling/Slugging Flow Behavior in a Cylindrical Fluidized Bed: ECT Measurements and Two-Fluid Simulations, Chemical Engineering Journal (2019), doi: https:// doi.org/10.1016/j.cej.2019.123120

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier B.V.

CEJ-D-19-00977 (Revised Manuscript) Bubbling/Slugging Flow Behavior in a Cylindrical Fluidized Bed: ECT Measurements and Two-Fluid Simulations Brajesh K. Singh, Shantanu Roy, Vivek V. Buwa* Department of Chemical Engineering, Indian Institute of Technology Delhi, New Delhi 110016, India (* corresponding author: [email protected])

Abstract While the Eulerian two-fluid models (TFM) are now reasonably well established to predict the “time–averaged” flow behavior of gas–solid flow in fluidized beds, quantitative predictions of inherently unsteady local solid-phase distribution and bubbling/slugging behavior using the TFM continues to be a challenge. In the present work, we have investigated time-evolution of solid volume fraction (s) distribution and bubbling/slugging behavior in a laboratory-scale fluidized bed using the Electrical Capacitance Tomography (ECT) measurements and TFM simulations. The ECT was used to measure time- and spatially- resolved local s distribution and these measurements were used further to characterize the bubble chord length fluctuations, bubbling frequency and bubble size distributions. The effects of different models used to calculate gas-solid momentum exchange through the drag force were investigated. The results showed that the predicted time-averaged s distribution differed marginally and were in a satisfactory agreement with the corresponding measurements. However, the predicted bubbling frequency and bubble size distribution, in particular, the formation of large bubbles/slugs, predicted using different drag models were considerably different and were not even in a qualitative agreement with the measurements. The modification proposed to the drag model in the present work led to an improvement in the predicted bubbling behavior and showed a qualitative agreement with the measurements. Further, simulations performed to investigate the effect of solid-phase frictional viscosity (s,fr) showed that an increase in the solid-phase frictional pressure leads 1

to an increase in s,fr. Such an increase in s,fr leads to compaction of the solid phase and to the formation of large bubbles/slugs, in line with what is observed experimentally. The predicted bubbling frequency and bubble size distribution were in a quantitative agreement with those measured using the ECT for different gas velocities. Keywords: Bubbling, Fluidized bed, Dynamics, Electrical Capacitance Tomography, Eulerian Two-fluid model, Frictional viscosity

2

1

Introduction Owing to excellent mixing, heat– and mass–transport characteristics; gas–solid fluidized

beds are extensively used to carry out gas–solid mixing and/or reactions in chemical process industry (e.g. coal and biomass combustion, gasification, etc.). Overall reaction rate, heat and mass transfer between the phases, and eventually the performance of the gas–solids fluidized bed often depend on the inherently unsteady (in time and space) bubbling/slugging behavior (e.g. bubble size distribution and bubbling frequency) and solid-phase distribution [1]. While the Eulerian two-fluid models (TFM) are now reasonably well established to predict the “time–averaged” flow behavior of gas–solid flow in fluidized beds, their ability to predict the time-evolution of solid volume fraction or bubbling/slugging quantitatively still remains to be established. Therefore, it is important to establish the capability of the TFM to predict the time-evolution of phase distribution by comparing with their corresponding time- and spatially-resolved phase distribution measurements. The Eulerian two-fluid model is widely used to simulate gas-solid flow in fluidized beds. In most of the previous studies, the predictions of the TFM were validated using the timeaveraged gas velocities or time-averaged solid volume fraction distributions along the radial or axial direction (e.g. [2–4]). Also, a comparison between the measured local/crosssectional averaged phase volume fraction fluctuations or pressure fluctuations with their corresponding simulations have been reported (e.g. [5–8]). In a few cases, the TFM was validated with time-averaged parameters and further it was used to study the dynamics of gas–solid flow (e.g. [9–11]), assuming that the validated model with time-averaged parameters can also provide accurate predictions of the dynamics of the flow. The literature reports cited above actually show that this is not the case i.e., the model validation of timeaveraged profiles does not imply predictability of the dynamics. For example, the predictions of the TFM were compared with the local solid-volume fraction measured using optic fiber 3

probe and pressure fluctuation measurements. It was found that the frequency and amplitude of pressure fluctuations agreed well with the predictions of the TFM, but differed when compared with the frequency and amplitude of measured solids volume fraction fluctuations [8]. Attempts have been made to compare the average bubble size predicted using the TFM simulations with the average bubble size calculated using different correlations (e.g. [12]). Owing to significant differences in the calculated average bubble size using different correlations (for the same fluid properties and bed geometry), the ability of the TFM to predict the bubble size distribution could not be established. Therefore, to quantify the abilities of the TFM to predict bubbling behavior, it is important to carry out one to one comparison between measured and predicted bubbling behavior. In one of the studies [13], the spatial distribution of the solid phase measured using ultra-fast X-ray tomography was used to calculate bubble size and bubble rise velocity. Also, the predictions of the TFM were compared with their corresponding measurements, and a significant deviation was observed between the predicted and measured bubble size distribution for most of the cases [13]. Such differences in the measured and simulated bubbling characteristics (size or frequency of the bubbles) can be attributed to the models used for interfacial (drag force) and viscous forces (solid stresses) in the TFM. This is a topic that demands a thorough study to establish the effect of these forces on the bubbling behavior. The calculation of gas-solid momentum exchange, which is accounted primarily through the drag force, is important for the prediction of gas-solid flow. Recently, Gao et al. [2] categorized the drag force models into heterogeneous (e.g. EMMS [14], Sarkar et al. [15], Radl et al. [16], etc.) and homogeneous (e.g. Gidaspow [17], Gao et al. [18], etc.) drag force on the basis of the formulation of the models (for details see Gao et al. [2]) and showed a comparison between the predictions of different drag force models [2]. The predicted timeaveraged axial and radial voidage distribution was found to be different for different 4

heterogeneous drag models. Also, the predictions using the homogeneous drag models were not in agreement with the corresponding measurements [2]. The homogeneous drag models are basically the combinations of different sub-models that are used for different regions in the bed. These regions are divided using different values of solid/gas volume fraction and one sub-model for each region is used. This was done due to significant variations in solid volume fraction in different parts of bed and use of any single sub-model for the entire range of solid/gas volume fraction was not sufficient to provide accurate results [19]. However, the division of the regions in the bed, where different sub-models were applied for calculation of gas-solid momentum exchange, was carried out rather on an ad-hoc basis. In addition to the drag force models, the solid viscosity is known to influence the bubbling behavior in the gas-solid fluidized bed (e.g. [20–22] ). The effects of different solid viscosities on the predicted bubble size, bubble detachment time, and time-averaged phase distribution were analyzed by different researchers (e.g. [20–22]). However, such studies were not validated with the measured bubble size distribution. It should also be noted that most importantly the studies were carried on the pseudo-3D rectangular beds. Due to the significant influence of the wall on the flow of gas and solids, the bubbling characteristics in pseudo-3D beds could be different than that of the beds with large depth (e.g. cylindrical beds). Therefore, it is important to investigate the effect of solid phase viscosity on the bubbling behavior in cylindrical beds. Several time–resolved measurement techniques have been developed and used for the measurement of local pressure fluctuations, local phase fraction fluctuations (optical fiber probe (e.g. [23,24]), capacitance probe (e.g. [25]), and spatial distribution of phase fraction fluctuations (e.g. ultrafast X–ray tomography and Electrical Capacitance Tomography (ECT)) in gas-solid fluidized beds. The measurements of local and instantaneous variations of the pressure or the volume fraction in the bed have been widely studied to derive the 5

information of local dynamics of the flow. While such local measurements can be used to understand the dynamics, flow regimes, etc., these local measurements do not provide complete information on bubble size distribution and its time-evolution. In order to obtain complete information, the spatial- and time-resolved measurements of phase volume fraction distribution are required. This can be achieved only with the time-resolved tomographic techniques (e.g. ultra-fast X-ray tomography [26], or ECT [27]). ECT is a non–intrusive and time-resolved tomographic technique which can be used for laboratory-scale as well as for industrial-scale fluidized beds. In addition to the timeaveraged solid volume fraction distributions (e.g. [28–30]), the ECT has been also used to characterize the bubbling behavior by measuring cross-sectional averaged solid volume fraction fluctuations (e.g. [30–32]). While such measurements of cross-sectional averaged solid volume fraction fluctuations are used for quantification of flow regimes and dynamic characteristics of the flow, it may not represent the dynamics of the flow accurately due to the spatial averaging over the bed cross-section. Also, the time-evolution of phase distribution has been studied, but without the quantification of the bubble size and bubbling frequency (e.g. [33]). There are only a few literature reports on the characterization of bubbling behavior performed using the time- and spatially-resolved solid volume fraction measurements (ECT [27] and ultra-fast X-ray tomography [13,26]). To establish the capability of the TFM to predict the bubbling behavior, there is a need for detailed information on time- and spatially-resolved measurement data and therefore the ECT measurements are performed in the present work. From the brief analysis of relevant literature presented above, it can be seen that there are only a few reports that compare the predicted instantaneous local or cross-sectional averaged flow properties with that of the measurements. The effect of drag models on the bubble size distribution was studied but without their corresponding validation with the measurements. 6

Further, the effect of different solid viscosities on bubble size was studied mostly on pseudo3D rectangular fluidized beds with their qualitative comparison using photographic images. To the best of author’s knowledge, quantitative predictions of bubbling/slugging behavior (in terms of bubbling frequency and bubble size distribution) using Euler two-fluid model still continues to be a challenge.

Figure 1. Schematic of the experimental set–up (A: compressor, B: pressure regulator, C: rotameter, D: pin valve, E: distributor, F, and G: ECT sensors, H: fluidized bed and I: ECT data acquisition system.

The objectives of the present work are (i) to investigate the ability of the TFM to predict time-evolution of phase distribution and the bubbling behavior, in terms of the bubbling frequency and bubble size distribution, in a cylindrical fluidized bed, (ii) to verify the predictions using the time- and spatially- resolved local s distribution, chord length fluctuations, bubbling frequency and bubble size distributions measured using the ECT and (iii) to comprehend the effects of gas-solid drag force and solid-phase viscosity on formation of large bubbles/slugs.

7

2

ECT Measurements Experiments were performed in a semi-batch cylindrical fluidized bed made of PMMA

(Db = 9 cm and Hb = 150 cm) using the ECT (see Figure 1). Glass beads (ds = 430 µm) of Geldart Group-B (Umf = 0.13 m/s) and air under ambient conditions were used as the working fluids. For all the experiments, the static bed height was kept at 22 cm and the superficial gas velocity (UG) was varied from 0.9 m/s (~7×Umf) to 2.6 m/s (20×Umf), using a uniform distributor. The uniform distributor consisted of 400 holes of 1 mm ID arranged in a square pitch of 4 mm over the entire column cross–section, which provided a fractional open area of 4.9%. A commercially available ECT measurement instrument (model: m3c, Industrial Tomography Systems plc, UK) was used to measure the permittivity distribution over a cross-section of the fluidized bed, at different axial locations, viz. z = 8.9, and 24 cm from the gas distributor. The ECT hardware was comprised of three components: sensor, data acquisition system (DAS) and image reconstruction software (see Figure 1). The sensors are comprised of 12 electrodes, mounted on the periphery at an equal interval (see ECT sensor schematic shown as an inset in Figure 1). The electrodes were connected to the DAS using coaxial cables and the DAS was connected to a computer via a USB interface. A commercial software m3c was used to provide initial data (injection current, permittivity, data acquisition time etc.) for every experiment. The adjacent measurement protocol was used in this study in which the electrical current was applied to an electrode (called source electrode) and the voltage was measured from rest of the electrodes (called detector electrodes). The voltage measurements were performed by taking all electrodes one by one as the source electrode and rest as detector electrodes. An injection current of 12 mA and frequency of 9,600 Hz were used in all experiments. There were 66 voltage measurements for one frame. Further, these voltage measurements were used for image reconstruction algorithm which 8

provided the permittivity distribution for each frame that consisted of 812 voxels on the sensor cross-section. The measurements carried out using the ECT technique for unary gas-solids flow have been benchmarked quantitatively by other measurement techniques (e.g. time-resolved X– ray [28], optic fiber probe [34], and MRI [35]), and a good agreement was reported. In some of the studies performed using the ECT, the offline iterative linear back projection (ILBP) method was used for the image reconstruction [28,34] and the results were compared with the time-resolved optic fiber probe [34] and X-ray tomography [28] and a good agreement was reported. In one of the recent study by Agarwal et al. [27], different image reconstruction algorithms (Landweber, NN-MOIRT, and LBP) were compared and LBP was found suitable for the reconstruction of tomograms. Therefore, the LBP algorithm was used in the present work for image reconstruction. In the present work, the measurements were carried out to investigate the effect of the UG on the time–evolution of solid volume fraction (𝛼𝑠) distribution on the bed cross-section. The time–evolution of bubble chord lengths was calculated at different axial locations from the tomograms and was analyzed to quantify the bubbling frequencies and the bubble chord length (size) distribution. The results are discussed in Section 4 together with the simulations results. The spatial resolution of the tomograms is estimated to be 5% of the column diameter and the data was recorded at 33 frames per second. There are 812 voxels in a single tomogram and permittivity value of each voxel represents a spatially average 𝛼𝑠 over the volume of one voxel (~0.45 × 0.45 × 3 cm3). The measured permittivity values were then converted to the solid volume fractions 𝛼𝑠 using the Maxwell equation [36] as,

m 

 g [2 g   s  2 g   g   s ]

(1)

[2 g   s   g   g   s ]

9

For air (εg) and solid (εs) system, εg = 1 s  1 

2 1   m    s 1   m 

(2)

2 1   s    m 1   s 

where αs is the solid volume fraction and εm is the measured permittivity of the mixture using the ECT. Table 1. Drag models used in the present work Gidaspow [17] If αg < 0.80, 1  150 If αg ≥ 0.80,

2 

 s 1   g   g  g d s2

 1.75

 g s U s  U g ds

(2)

3 CD  s g  g U s  U g  g2.65 4 ds

where, CD 

(1)

(3)

24 [1  0.15( g Re s )0.687 ] for Re s  1000  g Re s

 0.44 Re s 

for Re s  1000

g ds U s  U g g

Switch function given by Lu et al. [37] was used for smooth and rapid transition from Ergun[38] to Wen & Yu [39] drag equations

i 

tan1[1501.75( g  i )]



(4)

 0.5

where, 1 correspond to  1 = 0.8

K gs  (1  1 ) 1  1 2

(5)

Gao et al.[18] If αg < 0.80, 1  150

 s 1   g  g  g ds2

 1.75

g s U s  U g ds

gs g U s U g 5 If 0.80 < αg ≤ 0.933, 2  CD 72 ds (1   g )0.293 If 0.933 < αg ≤ 0.99 ,

If 0.99 < αg ≤ 1.0 ,

3  0.75CD

 4  0.75CD

(6)

 g s U s  U g  g2.65

(7)

(8)

ds

 g s U s  U g ds 10

(9)

CD 

24 [1  0.15(Re s )0.687 ] Re s

for Re s  1000

 0.44 Re s 

i 

(10)

for Re s  1000

 g g ds U s  U g g

tan1[1501.75( g  i )]



 0.5

where, 1 ,  2 and  3 correspond to  1 = 0.8,  2 = 0.933 and  3 = 0.99 respectively

K gs  (1  1 ) 1  1 (1  2 )  2  2 (1  3 ) 3  3  4 

(11)

Hong et al. [40] If αg ≤ εmf , K gs  150

 s 1   g  g  g ds2

 1.75

ds

If εmf < αg ≤ 0.5, Hd  183.4 1138.6 g  2355.2 g 1624.4 g

(13)

If 0.5 < αg ≤ 0.99, Hd  0.715  2.84 g  2.96 g

(14)

2

2

If αg > 0.99,

3

Hd 1

 3 CD  K gs  H d   s g  g U s  U g  g2.65   4 ds 

3

(12)

g s U s  U g

(15)

Computational model A two-fluid Eulerian CFD model was used to simulate unsteady gas-solid flow in the

cylindrical fluidized bed with the geometry and dimensions as considered for the bed used in the measurements. The volume- and time-averaged mass (Eqs. (S1) – (S2)) and momentum (Eqs. (S3) – (S4)) conservation equations as given in the Supplementary Table 1, are solved. The gas-phase and solid-phase shear stresses are modeled using the Eq (S5) and Eq. (S6), respectively. Also, Kgs term of the momentum conservation equations (Eqs (S3) and (S4) in the supplementary material) is calculated using different models given in Table 1 and using the modified drag model that is explained in the Section 4.2.1. The Kinetic Theory of Granular Flow (KTGF) model (Eq. (S8) in the supplementary material) was used to calculate granular temperature which was further used in the solid viscosity models to calculate solid stresses (kinetic and collisional) from equations given in Supplementary material (see Eqs. 11

(S14) – (S15)). The solid frictional viscosity was modeled using two different solid pressure equations and these frictional viscosity models were named as KTGF and Johnson-Jackson, respectively as given in the Supplementary material (see Eqs. (S16) – (S17)). Table 2. Boundary conditions, closure models, and simulation parameters Boundary conditions: Gas inlet Velocity inlet Outlet Pressure outlet (atmospheric pressure) Wall Gas No-slip Solid Johnson-Jackson KTGF closure models: Solid pressure ( p s ) Lun et al. [41] Solid bulk viscosity (  s )

Lun et al.[41]

Collisional viscosity (  s ,col )

Gidaspow et al. [42]

Kinetic viscosity ( s,kin )

Gidaspow et al. [42]

Frictional viscosity ( s, fr )

Schaeffer [43]

Granular diffusion ( k s )

Gidaspow et al. [42]

Granular dissipation (   s )

Lun et al. [41]

Radial distribution ( g o , ss )

Syamlal et al. [44]

Frictional pressure ( Ps , f )

KTGF-based [41]

Simulations parameter: 𝜌𝑠 𝑑𝑠 𝜌𝑔 𝜇𝑔 𝐻𝑏 𝐷𝑏 𝐻𝑠 𝑈𝑔 𝜑 𝑒𝑠𝑠 𝛼𝑠,𝑚𝑎𝑥 ∆𝑡

2500 kg/m3 430 µm 1.225 kg/m3 1.7894×10-5 Pa.s 1.35 m 0.09 m 0.22 m 0.9 – 2.6 m/s 0.01 0.9 0.65 1×10-4 s

12

The TFM model equations were solved using the commercial software “Ansys Fluent 14.0.0” for the cylindrical domain using closure models, solution parameters, and boundary conditions as given in Table 2. The phase coupled SIMPLE algorithm, which is an extension of the SIMPLE algorithm was used for the pressure-velocity coupling. For unsteady formulation, the 2nd order implicit scheme was used for discretization of time-derivatives and the QUICK scheme was used for discretization of convective terms. Each simulation was performed for a time period of 160 s whereas the time-averaged of solid volume fraction was obtained between 10 s and 35 s. A small time step of 1×10-4 was used to avoid numerical instability. The convergence criteria were set such that the residuals for all flow variables fall less than 1×10-4. In the simulations, solid phase (glass beads) was patched up to a height of 22 cm from the distributor which was the same as that used in the experiments. The initial velocity of the solid phase was kept at zero. The user defined functions (UDFs) were used for implementing different closures equations (Table 1) in the CFD model. Unless mentioned otherwise, all the simulations were performed using the closure models and the simulations parameters given in Table 2.

3.1 Effect of grid resolution The main focus of the present work is to characterize the local dynamics (using time varying spatial distribution of solid volume fraction) of gas-solid flow and therefore it is important to examine the dependency of local 𝛼𝑠 fluctuations using different time steps and grid resolutions. The simulations were performed using different number of cells (Ncells) = 11000 (coarse grid), 22,000 (medium grid), 46,000 (fine grid), and 92,000 (very fine grid). The time-averaged 𝛼𝑠 along the axial (see Figure S1(a)) and radial direction (see Figure S1 (b)) of the bed were compared. The time-averaged distribution was marginally different for coarse, medium, fine and very fine grids. Further, to examine the effect of the grid resolution

13

on the predictions of time-evolution of 𝛼𝑠, the power spectra of local 𝛼𝑠 fluctuations were compared for the medium grid, fine grid, and the very fine grid as shown in Figure S1(c), and marginal change in the frequency distribution was seen for fine and very fine grid, therefore the fine grid was selected for further studies in this work.

(a) (b) Figure 2. Comparison of measured and predicted time-averaged (a) particle dispersion height and (b) radial distribution of time-averaged 𝛼𝑠 at UG of 0.9 m/s and z = 8.9 cm using different drag models

4

Results and discussion

4.1 Comparison of measured and predicted time-averaged solid volume fraction

A comparison of measured and predicted dispersion (expanded bed) heights and timeaveraged 𝛼𝑠 distribution using different drag models is shown in Figure 2. The maximum and minimum dispersion heights were measured using the high-speed imaging at all UG and the spread bars are used to show the maximum and minimum dispersion heights. It can be seen that the spread bars increase with an increase in UG, due to increased fluctuations in the dispersion height at higher gas flow rates. The predicted dispersion height was in a good agreement at lower UG (< 1.3 m/s) using all the drag models, but at higher UG ( ≥ 1.3 m/s), considerable differences in the dispersion heights were observed for different drag models. This could be due to a different magnitude of drag force (as discussed in Section 4.2.1), 14

which primarily depends on the regions divided based on the local value of 𝛼𝑠 and the choice of drag models used in those regions. Further, the time-averaged 𝛼𝑠 was measured along the diameter of the bed at z = 8.9 cm (from the gas distributor plate) using the ECT and is compared with predicted time-averaged 𝛼𝑠 as shown in Figure 2(b). The time-averaged 𝛼𝑠 profiles predicted using different drag models were found to be marginally different than that of the measured profile, except that using the model of the Hong et al. [40].

(a) (b) (c) (d) (e) Figure 3. Time-evolution of 𝛼𝑠 distribution measured using (a) ECT and predicted using different drag models of (b) Gidaspow [17], (c) Gao et al. [18], (d) Hong et al. [40], and (e) modified model used in the present work (UG = 0.9 m/s, z = 8.9 cm).

As observed in the present work, the time-averaged 𝛼𝑠 distribution can be predicted reasonably well. This also agrees with the literature reports (e.g. [2]). However, as discussed earlier in the Section 1, it is important to predict the dynamic characteristics (e.g. timeevolution of local 𝛼𝑠 distributions, bubbling behavior, bubble size distribution) and compare the predictions with the literature. Such information is not available in the literature. These results are discussed in the following sections.

15

(a) (b) (c) (d) (e) Figure 4. Time-evolution of 𝛼𝑠 distribution measured using (a) ECT and predicted using different drag models of (b) Gidaspow [17], (c) Gao et al. [18], (d) Hong et al. [40], and (e) modified model used in the present work (UG = 1.3 m/s, z = 8.9 cm).

(a) (b) (c) (d) (e) Figure 5. Time-evolution of 𝛼𝑠 distribution measured using (a) ECT and predicted using different drag models (b) Gidaspow [17], (c) Gao et al. [18], (d) Hong et al. [40], and (e) modified model used in the present work (UG = 2.6 m/s, z = 8.9 cm).

16

17

4.2 Comparison of measured and predicted time-evolution of solid volume fraction The time-evolution of 𝛼𝑠 was measured using the ECT on the cross-section of the fluidized bed and the time-evolution of 𝛼𝑠 distribution along the diameter of the bed at z = 8.9 cm was used for further analysis. Also, the predicted time-evolution of 𝛼𝑠 distribution along the diameter was calculated and compared with the corresponding measurements (see Figure 3 - 6). The measured 𝛼𝑠 distribution showed mostly larger gas bubbles passing through the center of the bed. Also, due to smaller diameter of the bed, and flow conditions, single bubbling regime was seen at UG =0.9 m/s as shown in Figure 3(a). The time-evolution of predicted 𝛼𝑠 distributions along the bed diameter using the drag models Gidaspow [17], Gao et al. [18] and Hong et al. [40] are shown in Figure 3(b), (c), and (d), respectively. Similarly, the corresponding predicted 𝛼𝑠 distributions using aforementioned drag models at higher UG (1.3 m/s and 2.6 m/s) are shown in Figure 4(b)-(d) and 6(b)-(d), respectively. The bubbles predicted by the drag models of Gidaspow [17], Gao et al. [18] and Hong et al. [40] were smaller in size in comparison to the measured gas bubbles.

4.2.1 Modified drag model The observed differences in the time-evolution of 𝛼𝑠 distribution along the bed diameter (shown in Figure 3 – 5) can be attributed to the differences in the magnitude of drag forces predicted by the different drag models. As mentioned earlier, multiple regimes with varying value of 𝛼𝑠 can co-exist in the bed at a particular gas velocity [18], such co–existence of the multiple regimes necessitated the usage of multiple drag models with different classification criteria in terms of local 𝛼𝑠. However, the classification criteria on the use of different drag models in different regimes was done on an ad–hoc basis. For example, Gidaspow [17] used Ergun [38] equation for 𝛼𝑔 < 0.8 and the drag model proposed by Wen and Yu [39] for 𝛼𝑔 ≥ 0.8 (see Table 1). Though the validity of the Ergun equation was limited in the range of 0.4 18

< 𝛼𝑔 < 0.65 [45–47], it was used for 𝛼𝑔 upto 0.8. The possible reason behind the use of Ergun equation for 𝛼𝑔 upto 0.8 may be due to the fact that the Wen and Yu’s drag model is developed for single particle and it is inappropriate to use it for the dense phase. Following Gidaspow’s model that uses Ergun equation to calculate drag force up to 𝛼𝑔 = 0.8, the models those were developed further (e.g. [18,48]) use the Ergun equation up to 𝛼𝑔 = 0.8 and different drag models for 𝛼𝑔 > 0.8.

(a)

(b)

(c) (d) Figure 6. Effect of local slip velocities and gas volume fraction on the magnitude of Kgs for (a) Gidaspow [17], (b) Gao et al. [18], (c) Hong et al. [40], and (d) modified (present work) drag models.

From the preliminary calculations, it is found that the magnitude of drag force calculated using the Ergun equation is several times higher than that estimated by Wen and Yu’s model (see Figure 6(a)). Also, the use of Ergun for 𝛼𝑔 up to 0.8 implies larger regions in the bed with higher drag force (see Figure 6(a) and (b)), which may increase significantly 19

with the increase in the gas velocity and decrease in particle size. This could be one of the possible reasons for over-prediction of the dispersion height at higher gas velocities using Gidaspow’s [17] drag model and also using the models based on the Gidaspow’s [17] models. Table 3. Modified drag model used in the present work If αg ≤ 0.65 , 1  150

 s 1   g   g  g d s2

 1.75

 g s U s  U g ds

1.8  17.3   g s U s  U g  g If 0.65 < αg ≤ 0.94 ,  2  0.01  0.336    Re  ds g s  

If 0.94 < αg ≤ 1.0 ,  3  0.75CD

CD 

 g s U s  U g  g2.65

i 

(17)

(18)

ds

24 [1  0.15( g Re s )0.687 ] for Re s  1000  g Re s

 0.44

(16)

(19)

for Re s  1000

tan 1[150 1.75( g   i )]



(20)

 0.5

where, 1 and  2 correspond to 1 = 0.65 and  2 = 0.94, respectively

K gs  (1  1 )  1  1{(1   2 )  2   2  3 }

(21)

In the present work, Ergun [38] equation is used for the dense phase (𝛼𝑔 < 0.65), and the drag equation proposed by Gibilaro et al. [49] were used for sub–dense phase (0.65 < 𝛼𝑔 ≤ 0.94) with a modified scaling factor for the glass beads. In the literature, the scaling factor was obtained for FCC particles by comparing the bed expansion and time-averaged 𝛼𝑔 distribution (e.g. [50]). Since the scaling factor was not available for the glass beads, it was obtained by comparing the predicted bed expansion and time-averaged 𝛼𝑠 distribution with their corresponding measurements, the scaling factor for the present case was found to be 0.01. For the dilute phase (0.94 < 𝛼𝑔 ≤ 1.0), where single-particle drag model is applicable, the drag model of Wen and Yu [39] was used. In order to provide a smooth transition

20

between different drag models, a switch function [51] is used. The proposed drag model (this work) is summarized in the Table 3. The modified drag model showed the formation of larger gas bubbles together with smaller bubbles. The time-evolution of 𝛼𝑠 distribution predicted using the modified drag model showed an improved agreement with the measured time-evolution of 𝛼𝑠 distribution (see Figure 3). Similar improvements were observed over the other drag models at higher gas velocities (see Figure 4 and 5). The observed improvements in the predicted 𝛼𝑠 distribution using the modified drag model (see Figure 3(e) – 6(e)) can be attributed to the differences in the magnitude of drag forces. As it was already shown in Figure 6, the use of different drag models with varying dividing criteria results into different magnitudes of the drag force. To quantify such variations in the magnitude of drag force caused by the use of different drag models and dividing criteria, the fluctuations in the cumulative magnitude of drag force (sum of drag force in all the cells of the domain) and their averaged values (obtained using different drag models) at UG = 2.6 m/s are shown in Figure 7. The dispersion height was found to be directly proportional to the magnitude of drag force (see Figure 2(a) and Figure 7). For example, the average cumulative drag force is approximately same (~ 0.5 kg/m2s2) for the Gao’s drag model and the modified drag model and as a result the corresponding dispersion height are also approximately same (see Figure 2(a) and Figure 7). However, an over prediction in the average cumulative drag force in the case of Gidaspow’s drag model led to a considerable increase in the dispersion height and similarly a small under prediction in the average cumulative drag force in the case of Hong’s model led to a considerable decrease in the dispersion height.

21

Figure 7. Time-evolution of the cumulative drag force per unit volume in the domain predicted using different drag models at UG = 2.6 m/s.

However, the variations in bubble size and the number of bubbles per unit time shown in Figure 3 – 5 were difficult to correlate with the average cumulative drag force. The drag force acting in the cells covered by a bubble and the cells which are around that bubble can affect the bubble size. The drag force depends on the discrimination criteria and the drag sub-models used in those regions (see Supplementary Figures S2). The larger the region in the bed in which Ergun equation (e.g. 𝛼𝑔< 0.8) is applied, higher is the drag force magnitude in those regions. This reduces the formation of large bubbles by increasing upward momentum in a large region of the bed. Also, the use of Ergun equation in a large region reduces the fluctuations in the drag force which can directly influence the 𝛼𝑠 distribution in the bed (Figure 7 and Supplementary Figure S2).

22

Figure 8. Time-evolution of the cumulative total solid viscosity in the domain predicted using different drag models and KTGF based model at UG = 2.6 m/s.

Further, it is interesting to note that the use of different drag models in turn also affect the total solid viscosity (which is a sum of kinetic, collisional and frictional viscosities) in the domain as shown in Figure 8. In most of the earlier works on Eulerian two-fluid simulations reported in the literature on gas-solid fluidization, attempts were made to improve the predictions by using different drag models or by modifying the available drag models. It was shown in a few studies for pseudo-3D beds, without rigorous validation with their corresponding measurements, that the solid viscosity can also affect the bubble size and bubbling frequency [20–22]. Therefore, it is important to analyze the effect of solid viscosities on predicted time-evolution of 𝛼𝑠 distribution. From Figure 8, it can be seen that the average of cumulative solid viscosity is smaller for drag models of Gidaspow [17] and Gao et al. [18] in comparison to the drag models of Hong et al. [40] and the modified drag model. With smaller values of the solid viscosity in the cells around the bubbles, formation of bubbles is expected to be faster and the size of bubbles to be smaller in comparison to the system with relatively large viscosity values (e.g. compare Figure 5 and Figure 8), also see 23

the instantaneous total solid viscosity distribution along the diameter of the bed shown in Supplementary Figure S3. From the result presented above, it can be seen that the solid viscosity, in addition to the drag force contribute to the disagreement seen in the predicted results with the measurements. The drag force can also affect the solid viscosities and therefore the bubbling characteristics. Therefore the effect of solid viscosities on the bubbling characteristics is studied and the results are discussed in Section 4.4.

4.3 Comparison of measured and predicted bubbling characteristics 4.3.1 Bubbling frequency The time-evolution of 𝛼𝑠 distribution along the diameter of the bed was used to calculate the bubble chord-length fluctuations. Since the 𝛼𝑠 distribution along the diameter of the bed (or on bed cross-section) was measured using the ECT, the measurements of vertical chord lengths (in the axial direction) are not available (unlike in two-plane or volumetric ECT). At lower location (z = 8.9 cm) in the bed, with an increase in the UG from 0.9 to 1.3 m/s, the bubbling frequency was found to increase as shown in Figure 9(i)(a) and (b). With further increase in the amount of gas (UG = 2.6 m/s), the newly formed bubbles coalesce and lead to the formation of bigger bubbles or slugs. Also at upper location (z = 24 cm) in the bed, similar observations were seen (see Figure 9(ii)). The bubbling flow regimes observed at UG = 0.9 and 1.3 m/s (see Figure 3(a) and Figure 4(a)) agreed well with the flow regime map provided by Kunii and Levenspiel [52].

24

(i)

(ii)

Figure 9. Measured time-evolution of bubbles chord length fluctuations for UG of (a) 0.9 m/s, (b) 1.3 m/s, and (c) 2.6 m/s at (i) z = 8.9 cm and (ii) z = 24 cm.

Further, the power spectra of bubble chord length fluctuations (see Figure 9(i)) are shown in Figure 10. Since a distribution of frequencies around the dominant frequency was observed (see Figure 10(a) and (b)), the flow of bubbles was aperiodic at UG = 0.9 m/s and 1.3 m/s. The bubbling frequency at 1.3 m/s was higher than the bubbling frequency at 0.9 m/s. With further increase in the UG, bubbles were found to coalesce and also to break, thus leading to the presence of smaller as well as larger bubbles. Therefore, the power spectra did not show any particular dominant frequency (Figure 10(c)).

25

(a)

(b)

(c) Figure 10. Power spectra of measured bubble-chord length fluctuations at (a) 0.9 m/s, (b) 1.3 m/s and (c) 2.6 m/s ( z = 8.9 cm).

The power spectra of simulated time-evolution of bubble chord length fluctuations using different drag models at UG = 0.9 m/s are shown in Figure 11. Unlike the measured power spectra (at UG = 0.9 m/s see Figure 10(a)), the predicted power spectra showed multiple dominant frequencies for all the drag models except the model of Hong et al. [40]. Also, the predicted dominant frequency values (0.1 Hz – 7 Hz) were significantly higher for the drag models of Gidaspow [17] and Gao et al. [18] than the measurements (2.0 Hz). This was due to the formation of a large number of small bubbles of different sizes in the simulations in comparison to that in the measurements (see Figure 3). However, the higher frequencies were less prominent in the cases of Hong et al. [40] and the modified drag models as compared to the other drag models (see Figure 11). As discussed earlier, the drag models that lead to higher drag force and smaller solid viscosity, lead to the formation of 26

smaller bubbles with higher bubbling frequency. Therefore, the drag models of Gidaspow [17] and Gao et al. [18] showed higher bubbling frequency in comparison to that by Hong et al. [40] and the modified drag model.

Gao et al. [18]

Gidaspow [17]

Hong et al. [40]

Modified (present work)

Figure 11. Power spectra of simulated time-evolution of bubble-chord length for different drag models (UG = 0.9 m/s, z = 8.9 cm).

4.3.2 Bubble size distribution The measured time-evolution of bubble chord length crossing the bed diameter (at z = 8.9 and 24 cm) was calculated from the time-evolution of 𝛼𝑠 distribution along the diameter. The maximum chord length of a particular bubble was identified and considered as the horizontal bubble diameter (dh). The dh distribution was obtained from the measurements performed for ~150 s and it was independent of time. The effect of UG on the measured dh distribution is shown in Figure 12. With an increase in UG, it was seen that the dh was increased and the maximum value of dh was seen to reach to 9 cm, which is the internal diameter of the bed. At the lower location of the bed (z = 8.9 cm), a wider bubble size 27

distribution was observed and the average dh was found to increase with the increase in the UG (see Figure 12(i)). At the upper location of the bed (z = 24 cm), large bubbles or slugs were observed as shown in Figure 12(ii).

Figure 12. Measured bubble chord length distribution at different UG (a) 0.9 m/s, (b) 1.3 m/s and (c) 2.6 m/s at axial locations (i) 8.9 cm and (ii) 24 cm from the distributor.

While we have measured the horizontal bubble diameter in the present work, these measurements are different than that was reported in the literature using the local and crosssectional averaged measurements. The local measurements report the vertical chord length of the bubbles and the probe location will have a significant impact on the measured chord lengths due to the chaotic bubbling over the cross-section of the bed. On the other hand, in the cross-sectional averaged fluctuations that provide the horizontal size of bubbles, the cross-sectional averaging leads to loss of information on the actual behavior of bubbles. The measured bubble size distribution in the present work was further used to evaluate the effect of drag models and the solid phase viscosities on the bubble size distribution.

28

Figure 13. Predicted bubble chord length distribution using different drag models at z = 8.9 cm from the distributor and UG = 0.9 m/s.

The simulated dh distribution was calculated (following the same methodology as used for measurements) using the predicted time-evolution of solid volume fraction for different drag models (Figure 13 and Figure 14). The simulated dh distribution was also obtained from the simulated data for ~150 s and such distribution did not change further with time. The dh distribution calculated at UG = 0.9 m/s using different drag models showed a marginal difference (see Figure 13) and agreed well with the measured dh distribution except for larger bubble sizes (dh > ~ 7.0 cm) (see Figure 12(i)(a)). However, the presence of larger bubbles (dh > ~ 7.0 cm) was seen in the case of the modified drag model. The Sauter mean cord length (l32) was calculated as l32   li3 i

l

2 i

where is li is the chord length of ith

i

bubble and the summation is performed over the number of bubbles considered in an experiment or a simulation. The predicted Sauter mean chord lengths using the Gidaspow [17], Gao et al.[18], Hong et al.[40], and modified drag models were 4.66, 5.05, 5.10 and 29

5.89 cm, respectively. Among all the various drag models considered in the present work, the Sauter mean chord length predicted using the modified drag model (=5.89 cm) agreed well with the corresponding measured Sauter mean chord length (=6.46 cm).

Figure 14. Predicted bubble chord length distribution using different drag models at 2.6 m/s (z = 8.9 cm).

At higher UG ( = 2.6 m/s), the predicted dh distribution, in particular the % of larger bubble (dh > ~ 7.0 cm), was considerably different than that of the measured size for all the drag models except that predicted using the modified drag model (see Figure 14). The Sauter mean chord lengths predicted using Gao et al.[18], Hong et al.[40], and modified drag models were 6.03, 6.77 and 7.45 cm, respectively. Again, the Sauter mean chord length predicted using the modified drag model (=7.45 cm) agreed well with the measured Sauter mean chord length (=7.64 cm). At UG = 2.6 m/s due to large drag force in the bed which led to higher dispersion height, solid phase was highly dispersed and no well-defined bubbles were formed (see Figure 5(b)) for the drag model of Gidaspow [17]. Therefore, bubble size distribution could not be calculated for this case. The distribution of dh predicted using 30

modified drag model was in good agreement with the measurements due to the improved calculation of the drag force and therefore the solid viscosity in comparison to the other drag models.

4.4 Effect of solid viscosities on the bubbling behavior From the preliminary simulation studies, it was seen that among the kinematic, collisional, and frictional viscosities of the solid phase, the effect of frictional viscosity on the bubbling behavior was found to be more dominant than the dynamic (kinematic and collisional viscosities) viscosity. There are a few models available for the calculation of solid frictional viscosity (𝜇𝑠,𝑓𝑟), e.g. Srivastava and Sundaresan [20] and Shaeffer [43]. In both the models, the 𝜇𝑠,𝑓𝑟 is directly proportional to the normal component of the frictional stress, which is also called frictional pressure (𝑃𝑠,𝑓). Depending on the choice of 𝜇𝑠,𝑓𝑟 model (Srivastava and Sundaresan [20], Shaeffer [43], etc.) and the 𝑃𝑠,𝑓 models (KTGF, Johnson and Jackson [53], and Syamlal et al. [44]), different combinations for the solid frictional stress model are used in the literature. In the present work, the Shaeffer [43] model was used to calculate the 𝜇𝑠,𝑓𝑟as given in Eq. (22).

s, fr 

Ps, f sin  2 I2D

(22)

In all the simulations reported in this manuscript so far, the 𝑃𝑠,𝑓 was calculated using the KTGF based model as given by Eq. (23). Ps , f   s  s  s  2  s 1  ess   s2 g 0, ss  s

(23)

The KTGF based model used for calculations of 𝑃𝑠,𝑓 (referred as KTGF model to calculate 𝜇𝑠,𝑓𝑟) leads to a lower magnitude of 𝑃𝑠,𝑓 and therefore smaller values of 𝜇𝑠,𝑓𝑟, which results in the formation of smaller bubbles with higher bubbling frequency. In further simulations, 31

the model given by Johnson and Jackson [53] (see Eq. (24)) for 𝑃𝑠,𝑓 (referred as JohnsonJackson model to calculate 𝜇𝑠,𝑓𝑟) was used, which predict higher values of 𝜇𝑠,𝑓𝑟

Ps , f

if  s   s , f ,min

0  r    s   s , f ,min  F s   s ,max   s 

if  s   s , f ,min

(24)

where 𝐹 = 0.1𝛼𝑠, 𝑟 = 2, and 𝑠 = 5 are constants. The value of 𝑃𝑠,𝑓 depends on the 𝛼𝑠,𝑚𝑖𝑛and 𝛼𝑠,𝑚𝑎𝑥 and often different values of 𝛼𝑠,𝑚𝑖𝑛 (in the range of 0.6 – 0.63) and 𝛼𝑠,𝑚𝑎𝑥 (in the range of 0.6 – 0.65) are used in the study reported in the literature [20,22,54]. Also, it was reported that with increase in the 𝑃𝑠,𝑓, the voidage around the bubble increased and therefore the size of the bubble decreased (due to gas leakage) despite of having large value of 𝜇𝑠,𝑓𝑟 [21]. Whereas these predictions of bubbling behavior were compared with the photographic images qualitatively, the quantitative comparison in terms of measured and predicted bubble size distribution and bubbling frequency is not attempted before for cylindrical fluidized beds.

4.4.1 Time-evolution of phase distribution The effect of viscosity on the time-evolution of simulated 𝛼𝑠 distribution along the diameter of the bed with the corresponding measured time-evolution of 𝛼𝑠 distribution at UG = 0.9 m/s is shown in Figure 15. Two solid frictional viscosity models viz. KTGF (Eqs. 22 and 23) and Johnson-Jackson (Eqs. 22 and 24) was used and their corresponding predicted 𝛼𝑠 distribution is shown in Figure 15(ii)(a) and 15(iii)(a), respectively. With the use of Johnson-Jackson to calculate 𝜇𝑠,𝑓𝑟, the bubble size increased and bubbling frequency decreased in comparison to the KTGF based 𝜇𝑠,𝑓𝑟 model (see Figure 15(ii)(a) and 15(iii)(a)). This was seen due to the increase solid viscosity which led to the increase in the 32

compactness of the solid phase around the bubble (see Figure 15(ii)(b) and 15(iii)(b)). Further, the simulations were performed using the Johnson-Jackson model to predict the time evolution of 𝛼𝑠 at different UG (1.3 m/s and 2.6 m/s) and the results are shown in Supplementary Figure S4. The agreement between the measured and simulated timeevolution of 𝛼𝑠 was found to improve significantly over KTGF based 𝜇𝑠,𝑓𝑟 model (see Supplementary Figure S4).

(a) (b) (a) (b) (i) (ii) (iii) Figure 15: (i) Measured time-evolution of 𝛼𝑠 distribution along the diameter of the bed using ECT, predicted time-evolution of (a) 𝛼𝑠 distribution and (b) 𝜇𝑠 predicted using 𝜇𝑠,𝑓𝑟 models based on (ii) KTGF and (iii) Johnson-Jackson frictional viscosity models (UG = 0.9 m/s, z = 8.9 cm).

In addition to the time-evolution of 𝛼𝑠 distribution, the power spectra of bubble chord length fluctuations were also compared with that of the measurements as shown in Figure 16. The higher and widely distributed bubbling frequencies values (> 4 Hz) obtained by KTGF based 𝜇𝑠,𝑓𝑟 model (see Figure 16(b)) were drastically reduced with the increase in the 𝜇𝑠,𝑓𝑟 by the use of Johnson-Jackson model (see Figure 16(c)), and the latter showed a good agreement with the measured power spectra (see Figure 16(a)). Further, the predicted distribution of bubbling frequency and the dominant frequency ~1.8 Hz was found to be in a 33

good quantitative agreement with the measured bubbling frequency distribution and the corresponding dominant frequency ~2 Hz (see Figure 16(a) and (c)). Similarly, the power spectra of predicted bubble chord length fluctuations using the Johnson-Jackson based 𝜇𝑠,𝑓𝑟 model at higher UG were also calculated (see Supplementary Figure S5). At UG = 1.3 m/s, the predicted dominant frequency (~3 Hz) is found to be in a good agreement with the measured dominant frequency (see Supplementary Figure S5(a)). With further increase in the UG, the flow was found to be chaotic and therefore multiple frequencies were seen and were in a satisfactory agreement with the measurements (as shown in Supplementary Figure S5(b)).

(a)

(b) (c) Figure 16: Power spectra of time-evolution of bubble chord lengths (a) measured using ECT, predicted using (b) KTGF based 𝜇𝑠,𝑓𝑟 model and (c) Johnson-Jackson based 𝜇𝑠,𝑓𝑟 model (UG = 0.9 m/s, z = 8.9 cm).

34

(a)

(b) (c) Figure 17: Bubble chord length distribution (a) measured using ECT, predicted using (b) KTGF based model for 𝜇𝑠,𝑓𝑟 and (c) Johnson-Jackson based model for 𝜇𝑠,𝑓𝑟 (UG = 0.9 m/s, z = 8.9 cm).

4.4.2 Bubble size distribution Further, the bubble size distributions were calculated from the predicted timeevolution of bubble chord length using different models for 𝜇𝑠,𝑓𝑟 and are shown in Figure 17. As explained earlier, KTGF based 𝜇𝑠,𝑓𝑟 model lead to the formation of increased number of smaller bubbles in comparison to Johnson-Jackson based 𝜇𝑠,𝑓𝑟 model (see Figure 17). The bubble size distribution predicted using Johnson-Jackson based 𝜇𝑠,𝑓𝑟 model (Sauter mean chord length = 6.42 cm) was found in a good agreement with the corresponding measurements (Sauter mean chord length = 6.46 cm). Similarly, at higher UG = 1.3 m/s and 2.6 m/s, the Johnson-Jackson based 𝜇𝑠,𝑓𝑟 model led to satisfactory predictions of bubble size 35

distribution (see Supplementary Figure S6), in particular the formation of large bubbles. The bubbles size was found to increase due to increase in the solid viscosity arising from the change in the frictional pressure models (KTGF v/s Johnson-Jackson) as shown in Supplementary Figure S7(a) for UG = 0.9 m/s and in Supplementary Figure S7(b) for UG = 2.6 m/s. From the above discussion, it was seen that the solid viscosities play an important role in the bubbling behavior of the gas-solid flow in fluidized bed. Further, it is important to estimate appropriate values of 𝛼𝑠,𝑚𝑖𝑛 and 𝛼𝑠,𝑚𝑎𝑥 with a scientific basis and the particleresolved simulations can be very useful for this purpose.

5

Conclusions In the present work, we have investigated the time-evolution of local s distribution and

bubbling/slugging behavior in a laboratory-scale fluidized bed (Db = 9 cm) using the Electrical Capacitance Tomography (ECT) measurements and TFM simulations. Glass beads (ds = 430 µm) of Geldart Group-B were used and the superficial gas velocity was varied from 0.9 m/s (~7×Umf) to 2.6 m/s (20×Umf), using a uniform distributor. The ECT was used to measure time- and spatially- resolved local s distribution and these measurements were used to characterize the bubble chord length fluctuations, bubbling frequency and bubble size (horizontal diameter) distributions. The effects of different models, e.g. involving multi-regime models (depending on the local s values) and the EMMS-based models, used to calculate gas-solid momentum exchange through the drag force were investigated. The results showed that the predicted time-averaged s distribution differed marginally and was in a satisfactory agreement with the corresponding measurements, except the EMMS-based model of Hong et al. [40]. However, the predicted bubbling frequency and bubble size distribution using different drag 36

models, in particular, the formation of large bubbles/slugs, were considerably different and were not even in a qualitative agreement with the measurements. A modified drag model that relies on the use of Ergun [38] equation for the dense solid phase (𝛼𝑔 < 0.65) and other drag models for sub-dense and dilute solid phases with modified scaling and switching functions was used in the present work. This led to a reduction in the regions in the bed in which higher drag force was applied (by the use of Ergun equation), allowed the formation of larger bubbles/slugs and led to an improvement in the predicted bubbling behavior. It was also found that the application of different drag models in turn also affects the total solid viscosity and therefore the bubbling behavior. The posterior analysis showed that among the kinematic, collisional, and frictional viscosities of the solid phase, the effect of frictional viscosity on the bubbling behavior was found to be more dominant than the kinematic and collisional viscosities. Further, simulations performed to investigate the effect of solid-phase frictional viscosity (s,fr) showed that an increase in the solid-phase frictional pressure leads to an increase in s,fr. Such an increase in s,fr leads to compaction of the solid phase and to the formation of large bubbles/slugs, those are observed experimentally. The predicted bubbling frequency and bubble size distribution were in a quantitative agreement with those measured using the ECT for different gas velocities. While we have shown that the bubbling frequency and bubble size distribution can be predicted quantitatively well by controlling the s,fr for Geldart group B particles and for superficial velocities in the range of ~7Umf to 20Umf, considered in the present work, further work is required to extend these findings to particles of other Geldart regimes and also to understand the formulation of s,fr models through particle-resolved simulations.

37

Acknowledgments: The authors gratefully acknowledge the financial support provided by IIT Delhi for setting up the ECT system. One of the authors (Brajesh K. Singh) acknowledges the financial assistance provided through the “High Impact Research (HIR)” program funded by IIT Delhi (MI00808).

Nomenclature 𝐶𝐷

drag coefficient, -

𝐷𝑏

fluidized bed (internal) diameter, cm

𝑑𝑠

solid particle diameter, μm

𝑒𝑠𝑠

particle restitution co–efficient, -

𝑔

gravitational acceleration, m/s2

𝑔0,𝑠𝑠

radial distribution function, -

𝐻𝑏

bed height, cm

𝐻𝑑

heterogeneity index, -

𝐻𝑠

static bed height, cm

𝐼2𝐷

second invariant of the strain rate tensor, 1/s2

𝐼

identity tensor, -

𝐾𝑔𝑠

interphase momentum exchange co–efficient, kg/(m3∙s)

li

chord length, cm

l32

Sauter mean chord length, cm

𝑝𝑠

solid pressure, Pa

𝑃𝑠,𝑓

solid frictional pressure, Pa

𝑅𝑒𝑠

particles Reynolds number, -

𝑡

flow time, s

𝑈

velocity, m/s

𝑈𝐺

superficial gas velocity, m/s

Greek Letters α

instantaneous phase volume fraction, -

𝛼𝑖

gas volume fraction value for dividing criteria in the switch function, -

𝛼𝑠,𝑚𝑎𝑥 𝛼𝑠,𝑚𝑖𝑛

maximum solid phase volume fraction, minimum solid phase volume fraction, -

38

𝜇

viscosity, Pa∙s

ρ

density, kg/m3

𝜏

shear stress, N/m2

𝜀

relative permittivity, -

𝛽

interphase momentum exchange coefficient, kg/(m3s)

𝜑

specularity co–efficient, -

𝜆𝑠

solid bulk viscosity, Pa.s

𝜅𝜃,𝑠

granular diffusion co-efficient, kg/(m.s)

𝛾𝜃,𝑠

granular dissipation, W/m3



granular temperature, m2/s2



switch function, degree



angle of internal friction, degree

Subscripts col

collisional

fr

friction

𝑔

gas

kin

kinetic

𝑚

mixture

𝑚𝑓

minimum fluidization

𝑟

relative

𝑠

solid phase

𝑠𝑓

solid phase friction

Acronyms CFD

Computational Fluid Dynamics

ECT

Electrical Capacitance Tomography

ILBP

Iterative Linear Back Projection

KTGF

Kinetic Theory of Granular Flow

TFM

Two-fluid Model

References 39

[1]

S. Kaart, J.C. Schouten, C.M. Van Den Bleek, Improving conversion and selectivity of catalytic reactions in bubbling gas-solid fluidized bed reactors by control of the nonlinear bubble dynamics, Catal. Today. 48 (1999) 185–194.

[2]

X. Gao, T. Li, A. Sarkar, L. Lu, W.A. Rogers, Development and validation of an enhanced filtered drag model for simulating gas-solid fluidization of Geldart A particles in all flow regimes, Chem. Eng. Sci. 184 (2018) 33–51.

[3]

S. Hu, X. Liu, N. Zhang, J. Li, W. Ge, W. Wang, Quantifying cluster dynamics to improve EMMS drag law and radial heterogeneity description in coupling with gassolid two-fluid method, Chem. Eng. J. 307 (2017) 326–338.

[4]

X. Gao, C. Wu, Y.W. Cheng, L.J. Wang, X. Li, Experimental and numerical investigation of solid behavior in a gas-solid turbulent fluidized bed, Powder Technol. 228 (2012) 1–13.

[5]

F. Taghipour, N. Ellis, C. Wong, Experimental and computational study of gas-solid fluidized bed hydrodynamics, Chem. Eng. Sci. 60 (2005) 6857–6867.

[6]

A. Mahecha-Botero, T. Li, F. Haseidl, A. Nguyen, J.R. Grace, Experimental and computational fluid dynamic study of the change of volumetric flow in fluidized-bed reactors, Chem. Eng. Sci. 106 (2014) 231–241.

[7]

A. Busciglio, G. Vella, G. Micale, L. Rizzuti, Analysis of the bubbling behaviour of 2D gas solid fluidized beds. Part II. Comparison between experiments and numerical simulations via Digital Image Analysis Technique, Chem. Eng. J. 148 (2009) 145– 163.

[8]

M. De Vega, A. Acosta-iborra, C. Sobrino, F. Herna, Experimental and computational study on the bubble behavior in a 3-D fluidized bed, Chem. Eng. Sci. 66 (2011) 3499– 3512.

[9]

X. Gao, C. Wu, Y. wei Cheng, L. jun Wang, X. Li, Experimental and numerical investigation of solid behavior in a gas-solid turbulent fluidized bed, Powder Technol. 228 (2012) 1–13.

[10]

X. Lv, H. Li, Q. Zhu, Simulation of gas-solid flow in 2D/3D bubbling fluidized beds by combining the two-fluid model with structure-based drag model, Chem. Eng. J. 236 (2014) 149–157.

[11]

M. Deza, F. Battaglia, A CFD Study of Pressure Fluctuations to Determine Fluidization Regimes in Gas–Solid Beds, J. Fluids Eng. 135 (2013) 101301.

[12]

J. Li, B. Yang, CFD simulation of bubbling fluidized beds using a local-structuredependent drag model, Chem. Eng. J. (2017).

[13]

V. Verma, J.T. Padding, N.G. Deen, J.A.M.H. Kuipers, F. Barthel, M. Bieberle, M. Wagner, U. Hampel, Bubble dynamics in a 3-D gas-solid fluidized bed using ultrafast electron beam X-ray tomography and two-fluid model, AIChE J. (2014).

[14]

H. Chen, S. Gu, H. Li, Simulation gas-solid flow in the downer with new structurebased drag model, Powder Technol. 323 (2018) 163–175.

[15]

A. Sarkar, F.E. Milioli, S. Ozarkar, T. Li, X. Sun, S. Sundaresan, Filtered sub-grid constitutive models for fluidized gas-particle flows constructed from 3-D simulations, Chem. Eng. Sci. 152 (2016) 443–456. 40

[16]

S. Radl, S. Sundaresan, A drag model for filtered Euler-Lagrange simulations of clustered gas-particle suspensions, Chem. Eng. Sci. 117 (2014) 416–425.

[17]

D. Gidaspow, Multiphase flow and fluidization: continuum and kinetic theory descriptions, Academic Press, San Diego, 1994.

[18]

J. Gao, X. Lan, Y. Fan, J. Chang, G. Wang, C. Lu, C. Xu, CFD modeling and validation of the turbulent fluidized bed of FCC particles, AIChE J. 55 (2009) 1680– 1694.

[19]

J. Gao, X. Lan, Y. Fan, J. Chang, G. Wang, C. Lu, C. Xu, CFD modeling and validation of the turbulent fluidized bed of FCC particles, AIChE J. 55 (2009) 1680– 1694.

[20]

A. Srivastava, S. Sundaresan, Analysis of a frictional-kinetic model for gas-particle flow, Powder Technol. 129 (2003) 72–85.

[21]

D.J. Patil, M. Van Sint Annaland, J.A.M. Kuipers, Critical comparison of hydrodynamic models for gas-solid fluidized beds-Part II: Freely bubbling gas-solid fluidized beds, Chem. Eng. Sci. 60 (2005) 73–84.

[22]

S. Cloete, S.T. Johansen, A. Zaabout, M. van Sint Annaland, F. Gallucci, S. Amini, The effect of frictional pressure, geometry and wall friction on the modelling of a pseudo-2D bubbling fluidised bed reactor, Powder Technol. 283 (2015) 85–102.

[23]

P. Sauriol, H. Cui, J. Chaouki, Gas-solid structure in the vicinity of a sparger nozzle in a fluidized bed, Powder Technol. 228 (2012) 131–140.

[24]

H. Cui, N. Mostoufi, J. Chaouki, Characterization of dynamic gas-solid distribution in fluidized beds, Chem. Eng. J. 79 (2000) 133–143.

[25]

N. Ellis, H.T. Bi, C.J. Lim, J.R. Grace, Hydrodynamics of turbulent fluidized beds of different diameters, Powder Technol. 141 (2004) 124–136.

[26]

J. Ma, J.R. van Ommen, D. Liu, R.F. Mudde, X. Chen, E.C. Wagner, C. Liang, Fluidization dynamics of cohesive Geldart B particles. Part I: X-ray tomography analysis, Chem. Eng. J. (2018) 0–1.

[27]

V. Agrawal, Y.H. Shinde, M.T. Shah, R.P. Utikar, V.K. Pareek, J.B. Joshi, Estimation of Bubble Properties in Bubbling Fluidized Bed Using ECVT Measurements, Ind. Eng. Chem. Res. 57 (2018) 8319–8333.

[28]

C. Rautenbach, R.F. Mudde, X. Yang, M.C. Melaaen, B.M. Halvorsen, A comparative study between electrical capacitance tomography and time-resolved Xraytomography, Flow Meas. Instrum. 30 (2013) 34–44.

[29]

K. Dubrawski, S. Tebianian, H.T. Bi, J. Chaouki, N. Ellis, R. Gerspacher, R. Jafari, A. Kantzas, C. Lim, G.S. Patience, T. Pugsley, M.Z. Qi, J.X. Zhu, J.R. Grace, Traveling column for comparison of invasive and non-invasive fluidization voidage measurement techniques, Powder Technol. 235 (2013) 203–220.

[30]

G. Qiu, J. Ye, H. Wang, W. Yang, Investigation of flow hydrodynamics and regime transition in a gas-solids fluidized bed with different riser diameters, Chem. Eng. Sci. 116 (2014) 195–207.

[31]

B. Du, W. Warsito, L. Fan, Bed nonhomogeneity in turbulent gas‐solid fluidization, 41

AIChE J. 49 (2003) 1109–1126. [32]

X. Li, A.J. Jaworski, X. Mao, Comparative study of two non-intrusive measurement methods for bubbling gas-solids fluidized beds: Electrical capacitance tomography and pressure fluctuations, Flow Meas. Instrum. 62 (2018) 255–268.

[33]

C.E. Agu, L.A. Tokheim, M. Eikeland, B.M.E. Moldestad, Determination of onset of bubbling and slugging in a fluidized bed using a dual-plane electrical capacitance tomography system, Chem. Eng. J. 328 (2017) 997–1008.

[34]

T. Pugsley, H. Tanfara, S. Malcus, H. Cui, J. Chaouki, C. Winters, Verification of fluidized bed electrical capacitance tomography measurements with a fibre optic probe, Chem. Eng. Sci. 58 (2003) 3923–3934.

[35]

T.C. Chandrasekera, A. Wang, D.J. Holland, Q. Marashdeh, M. Pore, F. Wang, A.J. Sederman, L.S. Fan, L.F. Gladden, J.S. Dennis, A comparison of magnetic resonance imaging and electrical capacitance tomography: An air jet through a bed of particles, Powder Technol. 227 (2012) 86–95.

[36]

Maxwell J C, A Treatise on Electricity and Magnetism, 1881.

[37]

L. Huilin, D. Gidaspow, Hydrodynamics of binary yuidization in a riser: CFD simulation using two granular temperatures, Chem. Eng. Sci. 58 (2003) 3777–3792.

[38]

S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (1952) 89–94.

[39]

C.Y. Wen, Y.H. Yu, Mechanics of fluidization, Chem. Eng. Prog. Symp Ser. 62 (1966) 100–110.

[40]

K. Hong, Z. Shi, W. Wang, J. Li, A structure-dependent multi-fluid model (SFM) for heterogeneous gas-solid flow, Chem. Eng. Sci. 99 (2013) 191–202.

[41]

C.K.K. Lun, S.B. Savage, D.J. Jeffrey, N. Chepurniy, Kinetic theories for granular flow: inelastic particle in couette flow and slightly inelastic particle in a general flow field, J. Fluid Mech. 140 (1984) 223–256.

[42]

D. Gidaspow, R. Bezburuah, J. Ding, Hydrodynamics of circulating fluidized beds: kinetic theory approach, Fluid. VII, Proc. 7th Eng. Found. Conf. Fluid. (1991) 75–82.

[43]

D.G. Schaeffer, Instability in the evolution equations describing incompressible granular flow, J. Differ. Equ. 66 (1987) 19–50.

[44]

T.J. Syamlal, M., Rogers, W., O’Brien, Mfix documentation: volume I, theory guide. Technical Report DOE/METC-9411004, NTIS/DE9400087, 1993.

[45]

C. K. Zoltani, Flow Resistance in Packed and Fluidized Beds: An Assessment of Current Practice, U.S. Army Ballist. Res. Lab. Maryl. USA. Technical (1992).

[46]

S. Ergun, Fluid Flow Through Packed Columns, Fluid Flow Through Packed Columns. 48 (2). (1952) 89–94.

[47]

I.F. Macdonald, M.S. El-Sayed, K. Mow, F.A.L. Dullien, Flow through Porous Media - The {Ergun} Equation Revisited, Ind. Chem. Eng. Fundam. 18 (1979) 199–208.

[48]

Y. Zhang, J.M. Reese, The drag force in two-fluid models of gas-solid flows, Chem. Eng. Sci. 58 (2003) 1641–1644.

[49]

L.G. Gibilaro, R.D.I. Felice, P.U. Foscolo, Generalized Friction Factor and Drag 42

Coefficient Correlations for Fluid Particle Interactions, Chem. Eng. Sci. 40 (1985) 1817–1823. [50]

T. McKeen, T. Pugsley, Simulation and experimental validation of a freely bubbling bed of FCC catalyst, Powder Technol. 129 (2003) 139–152.

[51]

L. Huilin, D. Gidaspow, Hydrodynamics of binary uidization in a riser : CFD simulation using two granular temperatures, Energy. 58 (2003) 3777–3792.

[52]

D. Kunii, L. Octave, Fluidization engineering, Elsevier, 2013.

[53]

P.C. Johnson, R. Jackson, Friction collisional constitutive relations for granular materials, with application to plane shearing, J. Fluid Mech. 176 (1987) 67–93.

[54]

A. Passalacqua, L. Marmo, A critical comparison of frictional stress models applied to the simulation of bubbling fluidized beds, Chem. Eng. Sci. 64 (2009) 2795–2806.

43

44

Bubbling/slugging flow behavior in a cylindrical fluidized bed characterized using ECT Local αs distribution, bubbling frequency and bubble size distribution predicted using different drag models vary significantly. Increase in solid frictional viscosity leads to compaction of solid phase and formation of large bubbles/slugs Quantitative predictions of bubbling/slugging behavior using modified drag and frictional viscosity

45