CFD Modelling of Fluidized Bed with Immersed Tubes

Fluidized beds are widely used in combustion and chemical industries. The immersed tubes are usually used for enhancement of heat transfer or control of temperature in fluidized beds. By his turn, tubes subjected to the solid particle impact may suffer severe erosion wear. Many investigations have been devoted to erosion in tubes immersed in fluidized beds on the various influencing factors (cf. Lyczkowski and Bouillard, 2002). As pointed by Achim et al. (2002), the factors can be classified as particle characteristics, mechanical design and operating conditions.


Introduction
Fluidized beds are widely used in combustion and chemical industries. The immersed tubes are usually used for enhancement of heat transfer or control of temperature in fluidized beds. By his turn, tubes subjected to the solid particle impact may suffer severe erosion wear. Many investigations have been devoted to erosion in tubes immersed in fluidized beds on the various influencing factors (cf. Lyczkowski and Bouillard, 2002). As pointed by Achim et al. (2002), the factors can be classified as particle characteristics, mechanical design and operating conditions. Some previous experimental studies have focused on bubble and particle behaviors (Kobayashi et al., 2000, Ozawa et al., 2002, tube attrition, erosion or wastage (Bouillard and Lyczkowski, 1991;Lee and Wang, 1995;Fan et al., 1998;Wiman, 1994), heat transfer (Wong andSeville, 2006, Wiman andAlmstedt, 1997) and gas flow regimes (Wang et. al, 2002).
Previous numerical studies were also performed using different CFD codes. Recently He et al. (2009He et al. ( , 2004, using the K-FIX code adapted to body fitted coordinates investigated the hydrodynamics of bubbling fluidized beds with one to four immersed tubes. The erosion rates predicted using the monolayer kinetic energy dissipation model were compared against the experimental values of Wiman (1994) for the two tube arrangement. The numerical values were three magnitudes lower than the experimental ones. Also employing a eulerian-eulerian model and the GEMINI numerical code, Almstedt (2000, 1999) performed numerical computations and comparison against experimental results (Enwald et. al., 1999). As reported for those authors, fairly good qualitative agreement between the experimental and numerical erosion results were obtained, and the contributions to the erosion from the different fluid dynamics phenomena near the tube were identified.
In the present study, we revisit the phenomena of the immersed tubes in a gas fluidized bed with different immersed tube arrangements employing the eulerian-eulerian two fluid model and the MFIX code. The purposes of the numerical simulations are to compare and explore some effects not previously investigated in the above-mentioned references.

Two fluid and erosion model
The mathematical model is based on the assumption that the phases can be mathematically described as interpenetrating continua; the point variables are averaged over a region that is www.intechopen.com large compared with the particle spacing but much smaller than the flow domain (see Anderson, 1967). A short summary of the equations solved by the numerical code in this study are presented next. Refer to  and Syamlal et al. (1993) for more detailment.
The continuity equations for the fluid and solid phase are given by : In the previous equations  f ,  s ,  f ,  s , f v  and s v  are the volumetric fraction, density and velocity field for the fluid and solids phases.
The momentum equations for the fluid and solid phases are given by: (4) f S s S are the stress tensors for the fluid and solid phase. It is assumed newtonian behavior for the fluid and solid phases, i.e., In the above equation P, ,  are the pressure, bulk and dynamic viscosity, respectively.
In addition, the solid phase behavior is divided between a plastic regime (also named as slow shearing frictional regime) and a viscous regime (also named as rapidly shearing regime). The constitutive relations for the plastic regime are related to the soil mechanics theory. Here they are representated as: In the above equation * ε is the packed bed void fraction and  is the angle of internal friction.
A detailing of functions f 1 to f 4 and f 9 can be obtained in Benyahia (2008).
On the other hand, the viscous regime behavior for the solid phase is ruled by two gas kinetic theory related parameters (e,  The solid stress model outlined by Eqs. (6) and (7) will be quoted here as the standard model. Additionally, a general formulation for the solids phase stress tensor that admits a transition between the two regimes is given by : According to Pannala et al.(2009), two diferent formulations for the weighting parameter "" can be employed : In the above equation the void fraction range  and the shape factor  are smaller values less than unity. It must be emphasized that when  goes to zero and  equals to unity, the "switch" model as proposed by Syamlal et al. (1993) based on the Schaeffer (1987) can be recovered. The models based on eqs. (9a) and (9b) will be referred in the numerical simulations as BLEND S and BLEND T, respectivelly.
On the other hand, the Srivastava and Sundaresan (2003), also called "Princeton model", can be placed on the basis of Eq. (9) Also in equations (4) and (5) fs I  is the momentum interaction term between the solid and fluid phases, given by There is a number of correlations for the drag coefficient  (Eqs. 11 to 16). The first of the correlations for the drag coefficient is based on Wen and Yu (1966) work. The Gidaspow drag coefficient is a combination between the Wen Yu correlation and the correlation from Ergun (1952). The Gidaspow blended drag correlation allows controlling the transition from the Wen and Yu, and Ergun based correlations. In this correlation the  blending function was originally proposed by Lathowers and Bellan (2000) and the value of parameter C controls the degree of transition. From Eq. (14), the correlation proposed by Syamlal and O'Brien (1993) carries the advantage of adjustable parameters C 1 and d 1 for different minimum fluidization conditions. The correlations given in Eq. (15)  Re 1000 For closing the model, a transport equation for the granular energy  provides a way of determine the pressure and viscosity for the solid phase during the viscous regime. Equation (17) is a transport equation for the granular energy . Its solution provides a way of determine the pressure and viscosity for the solid phase during the viscous regime. The terms  s  and  gs are the granular energy conductivity, dissipation and exchange, respectively.
  ss s ss In the algebraic approach, instead solving the full equation (17) , the granular energy is obtained by equating the first term on the right hand side with the dissipation term.
The model where Eqs. (5) to (8) and (17) are solved is the kinetic theory model, termed here as KTGF. Conversely, in the constant solids viscosity model (CVM) the solids pressure is defined as in Eq. (6) and the solids viscosity in either plastic and viscous regimes is set equal to a constant.

www.intechopen.com
For erosion calculations in this work we use the monolayer energy dissipation model (Lyczkowski and Bouillard, 2002). In that model the kinetic energy dissipation rate for the solids phase in the vicinity of stationary immersed surfaces is related to erosion rate in m/s by multiplication with an appropriate constant. This constant is function of surface hardness, elasticity of collision and diameter of particles hitting the surface. The kinetic energy dissipation rate  s in W/m 3 for the solids phase is given by : The hydrodynamic model is solved using the finite volume approach with discretization on a staggered grid. A second order accurate discretization scheme was used and superbee scheme was adopted for discretization of the convective fluxes at cell faces for all equations in this work. With the governing equations discretized, a sequential iterative solver is used to calculate the field variables at each time step. The main numerical algorithm is an extension of SIMPLE. Modifications to this algorithm in MFIX include a partial elimination algorithm to reduce the strong coupling between the two phases due to the interphase transfer terms. Also, MFIX makes use of a solids volume fraction correction step instead of a solids pressure correction step which is thought to assist convergence in loosely packed regions. Finally, an adaptive time step is used to minimize computation time. See Syamlal (1998) for more details. The immersed obstacles were implemented using the cut-cell technique available in the code (Dietiker, 2009) The numerical runs were based on experiments of Wiman (1994Wiman ( , 1997 in an air pressurized bed with horizontal tubes for four different tube-bank geometries. The T2 and I4 are portrayed in Fig. 1 whereas the S4 and S4D tubes geometry are given in Figure 3. Figure 2 details the domain for the numerical simulations and circumferential angle coordinate for erosion measurement. In Fig. 2a is detailed the mesh for the T2 arrangement This nonuniform stretched grid, follows Cebeci et al. (2005) approach, has fine spacing close to the surface of the tubes and coarse spacing away from the surface. For all the other arrangements the mesh was uniform as depicted in Fig. 2b. The mesh employed for the bed without tube, and the I4 arrangement was 60  260 cells, and for the S4 and S4D arrangement 60  340. The bed was operated at room temperature (24 C) at pressures between 0.1 and 1.6 MPa and at two different excess velocities: U f1 -U mf = 0.2 m/s and U f1 -U mf = 0.6 m/s. Here, U f1 is the superficial fluidization velocity based on the free bed crosssection. The voidage at minimum fluidization was 0.46 and the minimum fluidization velocities was 0.42 , 0.31 and 0.18 m/s, for 0.1, 0.4 and 1.6 MPa pressures, correspondingly. The particle diameter and density were 700 m and 2600 kg/m 3 . The bubble parameters www.intechopen.com obtained from simulation were based on methodology described in Almstedt (1987), using numerical probes in the domain, centered at (0.15, 0.55) and separated from 15 mm. The target tube for erosion measurements is the one centered at (0.18, 0.55) for all the tube arrangements. More information about the experiments can be accessed from Wiman (1994Wiman ( , 1997. In this work, the parameters for controlling the numerical solution (e.g., under-relaxation, sweep direction, linear equation solvers, number of iterations, residual tolerances) were kept www.intechopen.com as their default code values. Moreover, for setting up the mathematical model, when not otherwise specified the code default values were used. The computer used in the numerical simulations was a PC with OpenSuse linux and Intel Quad Core processor. The simulation time was 20 s.
For generating the numerical results, e employed the parameters listed above, referred here as baseline simulation. In addition, for the baseline simulation we employed the Syamlal-O´Brien drag model, the standard solid stress model, and slip and non-slip condition for solid and gas phase, correspondingly. The previous set of models will be referred in the result's section as baseline simulation models. (4) shows the bubble splitting mechanism taking place and the characteristic time scale for the bubble passage. After the bubble passage, the solid wake has higher solid velocity magnitude around the obstacle. The solids upward movement following the bubble wake and wall downward movement is kept unchanged. Generally, for beds without internal obstacles, bubble size increases with bed height, particle size and superficial velocity. Analysis of CFD results shows that the presence of tubes was found to alter such general trends for bubble growth.  Also, analysis of transient simulated results shows the highest values of kinetic energy dissipation rate occur when the particle fraction suddenly changes from a low to a high value, which corresponds to the tube being hit by the wake of a bubble. Figure 6 presents the result of kinetic energy dissipation for different solid stress models. As shown the peak values of the dissipation rate are close to 180 degrees for the model with (baseline) and without the blending function discussed in section 2. On the other hand, the maximum value predicted by the constant solids viscosity model (NO KTGF) locates around 180 degrees, and it is superior to the predicted by the solids kinetic energy theory.   Figures 8 and 9 shows a comparison between the numerically predicted time averaged values of the kinetic energy dissipation rate, using the baseline simulation models discussed in section 2, and those based on the experimental values for the T2 arrangement in the work of Gustavson and Almstedt (2000). A comparison of the numerical and experimental results for the two different operational pressures (c.f . Figs 8 and 9 ) shows some degree of discordance. However, the lack of minutely agreement with experimental values, the results can be compared for recognizing similar drifts. For instance, from the Fig 8, the higher simulated values of the dissipation rate occurring in the bottom position of the tube, i.e. for  < 180 o , is in agreement with the experimental counterpart. Similarly, the increase of dissipation rate with increased operational pressure for the simulated results is in agreement with the experimental values. Regarding the erosion and baseline models used for the simulation, some remarks towards better agreement with experimental values can also be done. According to the monolayer erosion model and its discussion above Eq. (19) some degree of uncertainty is associated to the multiplying constant, as the exact value of elasticity of collision is not known. It is also expected, that adjustments in the baseline simulation models, such as those outlined in Figs. (5) to (7) would produce better agreement. Figure 10 is a sampling plot showing the instantaneous gas volumetric fraction fields for different tube arrangements. Analysis of Fig. (10) shows the influence of immersed obstacles on the bubble splitting mechanism taking place and the bubble passage pattern. Above the tube bank, the bubble appears to grow to size similar to the without tube geometry. For the geometry with tubes the bubble encompasses the obstacles but not at the full width of the bed. The interaction is stronger for the denser tube geometry. Figures 11 and 12 shows a comparison between the numerically predicted values of bubble frequency, using the baseline simulation models discussed in section 2, and those based on the experimental measurements from the works of Almstedt (1987) and Wiman (1995). As shown in the Fig. 10, the calculated values of Nb are underestimated at higher pressures, while at low pressures, there is a quite good agreement between calculated and experimental results. This conclusion holds true both for the I4 and for the S4D tube arrangement. As in the experimental results there are no noticeable differences for the two tubes arrangements. Figures 13 and 14 shows a comparison for the bubble frequency for the bed without tube and with the S4 arrangement. As shown in Fig. 13 for the bed without tubes the trend points to a frequency agreement between 0.4 and 0.6 MPa. Up this range the values differences increases as pressures increases up to 1.6 MPa. For the S4 arrangement the trend is similar to the I4 and S4D arrangement. The experimental results depicted in Figs. 11 to 14 suggest that the frequency increases with pressure both with and without tubes. For numerical results this holds true only for low pressures, i.e., 0.1 and 0.4 MPa. The numerical results suggest a maxima occurring between 0.4 and 1.6 MPa. On the other hand, the numerical results corroborate the experimental trend that the mean frequency is higher for the bed with tubes than for the freely bubbling bed.

Conclusion
In this work was investigated numerically the hydrodynamics of two dimensional beds with immersed tubes. The simulations were based on an experimental bed with different tube bank geometries, operating pressures and gas excess velocities. The objective of this study was two fold: explore and investigate some effects not previously explored in the literature, to verify the feasibility of the MFIX code for such a kind of study.
The simulation's results were framed in terms of averaged solids kinetic energy dissipation rate and bubble parameters (frequency and mean velocity). A comparison between the numerical results for frequency and the experiments shows good agreement for low pressures. Also, some similar drifts were identified, e.g., greater frequency for greater excess velocities. By his turn, the bubble velocity numerical results agrees better with experiments for high pressures. For the bed with tubes, a maxima of bubble velocity for intermediate pressure is identifiable, this could not be verified for the numerical results. Our results points to significant influences in the predicted dissipation rates and consequently, in the erosion rate, when employing different drag models. The dissipation rate is also influenced by either the use of blending functions for the transition between the plastic and viscous regime of solids flow or the use of a constant viscosity model. The results are key sensitive to the slip condition for the gas phase in the surface of the tube. In the case of a free slip condition for the gas phase the lowest values of dissipation are obtained By his turn, for the solids kinetic dissipation energy, the range of angles that gives maximum values are shifted in relation to the experiments. Finally, remarks towards better agreement with experimental values can also be done for the energy dissipation model. Specifically, according to the monolayer erosion model and its discussion above Eq. (19) some degree of uncertainty is associated to the multiplying constant, as the exact value of elasticity of collision is not known.