Gas-solid two-phase flow in the riser of circulating fluidized beds: mathematical modelling and numerical simulation

Gas-Solid Two-Phase Flow in the Riser of Circulating Fluidized Beds: Mathematical Modelling and Numerical Simulation

Luben Cabezas Gómez

[e-mail]

Fernando Eduardo Milioli

[e-mail]

Nucleo de Engenharia Térmica e Fluidos

Escola de Engenharia de São Carlos, USP

Av. Trabalhador São-carlense, N° 400, Centro

13566-590 São Carlos, SP. Brazil

A mathematical model is developed for gas-solids flows in circulating fluidized beds. An Eulerian formulation is followed based on the two-fluids model approach where both the fluid and the particulate phases are treated as a continuum. The physical modelling is discussed, including the formulation of boundary conditions and the description of the numerical methodology. Results of numerical simulation are presented and discussed. The model is validated through comparison to experiment, and simulation is performed to investigate the effects on the flow hydrodynamics of the solids viscosity.

Keywords: Two-fluid model, numerical modelling, two-phase gas-solid flow, circulating fluidized bed

Introduction

Gas-solids fluidization is widely applied in industry, including petroleum, chemical, metallurgical and energy industries (Kunii and Levenspiel, 1991). The largest applications of fluidized bed reactors occur in coal combustion for large scale thermoelectric power generation, and catalytic cracking of petroleum to produce gasoline and other fuels.

Circulating fluidized bed reactors are very competitive both technologically and economically for energy generation from coal combustion, and for gas fuel production from biomass gasification. The circulating fluidized bed combustors (CFBC) present several advantages such as fuel flexibility, wide range of possible operational conditions, high combustion efficiency, low emissions of both NO_x and SO₂, and good recirculation rates of solids through the riser. The CFBC are called the second generation of fluidized bed combustors (Tsuo, 1989).

In spite of the huge application, the CFBC technology faces some problematic features which need to be solved, such as the strong erosion on internal surfaces, and the blocking of fine particles in the riser. Also, the CFBC design results very complex in view of scaling-up difficulties, owing to the high sensibility of the flow regarding scale and operational conditions (Ding, 1990).

Hydrodynamic studies can significantly contribute for better understanding and thereby solving the above problems. The structure of multiphase flows in circulating fluidized bed columns is very complex, showing great variations on solids volumetric fraction through the riser, continuous formation and dissipation of clusters, and a high rate of recirculation of solids. Such flow conditions causes an intense superficial contact between gas and solids providing high reaction rates. The knowledge of the flow hydrodynamics is then of great importance, allowing relevant reactive and mass transport parameters to be determined.

The mathematical modelling of gas-solids fluidization processes represents an ancillary tool for minimizing the experimental efforts required for developing industrial plants. Experiment and prototype development are, of course, the main requirements for accurate engineering design in any industrial process. However, mathematical modelling and numerical simulation are in continuous development, contributing in a growing way for the better understanding of processes and physical phenomena, and thereby for design. Besides, mathematical models require experiment in order to be validated and, concerning fluidization, the required experiments involve complex measurements of difficult accomplishment. Therefore, mathematical modelling also represents an incentive for the development of new experimental methods and techniques.

In this work mathematical modelling is developed for the gas-solids flow through the riser of circulating fluidized beds using the two-fluids model approach. Problematic features related to this line of modelling include the treatment of pressure and viscosity of the solid phase, the formulation of reliable boundary conditions, and aspects related to physical instabilities of the flow as well as instabilities introduced by numerical procedures. In this work simulation is performed aiming to study the effects on the flow hydrodynamics of the solids viscosity.

Nomenclature

C Velocity of sound in gas, [m/s]
CD Drag coefficient for a particle on a particulate suspension
CDs Drag coefficient for a particle on an infinite medium
dp Particle mean diameter, [m]
D Gas phase residuous defined on Eq. (57)
f General function defined on Eq. (7)
Fdrag Stationary drag force on the interface, [N]
g Gravity acceleration, [m/s2]
G Solids phase elasticity modulus, [Pa]
Gs Solids mass velocity, [kg/(m2-s)]
Unitary tensor
Lp Particles mean free path, [m]
MkI Term for interface momentum transfer for phase k, [N/m3]
Generalised interface drag force for phase k, [N/m3]
Interface mass flux for phase k, [kg/(m2-s)]
External surface normal unitary vector for phase k
np Number of particles per unitary volume, [1/m3]
P Thermodynamic pressure, [Pa]
Pg Gas phase pressure, [Pa]
Ps Solids phase pressure, [Pa]
r Radius, point position vector, [m]
Rw Column internal radius, [m]
Re Reynolds number
t Time, [s]
T Temperature, [K]
Stress tensor for phase fase k, [Pa] Uk Mean velocity for phase k, [m/s]
Ug, Us Mean velocities of the gas and solids phases, respectively, [m/s]
Ur Relative velocity between the phases, [m/s]
U0 Slipping velocity, [m/s]
uI Interface velocity, [m/s]
uk Instantaneous velocity of phase k, [m/s]
Mean velocity fluctuation for phase k, [m/s]
V Volume, [m3]
vk Volume of phase k, [m3]
W Weighing factor, wall
x Transversal direction, [m]
Xk Phase density parameter
y Normal direction on Eqs. (49) and (51), [m]
z Axial direction, [m]

Greek symbols

ag Void fraction

ak Volumetric fraction of phase k

as Volumetric fraction of solids phase

b Interface drag function, [kg/m3-s]

dt Time interval, [s]

fk Volumetric source term

fI Volumetric source term at the interface

fs Particle sphericity

h Coefficient defined on Eq. (58)

lg, ls Volumetric viscosities of the gas and solids phases, respectively, [kg/m-s]

lk Volumetric viscosity of phase k, [kg/m-s]

lp Average distance among particles, [m]

mg, ms Dynamic viscosities of the gas and solids phases, respectively, [kg/m-s]
mk Dynamic viscosity of phase k, [kg/m-s]
rg, rs Densities of the gas and solids phases, respectively, [kg/m3]
rk Density of phase k, [kg/m3]
rp Particulate density, [kg/m3]
Viscous stress tensor for gas and solids phases, respectively, [Pa]
Viscous stress tensor for phase k, [Pa]
yk Conserved intensive property

Subscripts

g Gas phase

I Interface

k Phase index (1-gas, 2-solids)

s Solids phase

p Particulate

r Relative

1, 2 Gas and solids phase, respectively

Superscripts

cin kinetics

col Collision

gas Gas

w Weighing

® Vector

Tensor

‘ Fluctuation

General Lines of Mathematical Modelling in Circulating Fluidization

Harris and Davidson apud Pugsley and Berruti (1996) present a classification of models for simulating circulating fluidized beds. According to the authors there are three types of mathematical models:

1. Models that predict the axial variation of the density of solids and disregard its radial variations.

2. Models that predict the radial variation of the density of solids and the high average slipping velocities, accounting for two or more regions of different flow characteristics (such as annular plug flow or annular flow with clusters).

3. Models that apply the fundamental conservative equations of the fluid dynamics for predicting the two-phase gas-solids flow.

The first two types of models are mostly used for design, mainly for investigating the effects on the process of geometry and operational conditions. These models can easily include chemical kinetics for simulating the performance of reactors (Pugsley and Berruty, 1996).

The models of the third type, as for example a two-fluids model, are more suitable for research allowing, for instance, the behaviour of flow local structures and the effects of local geometry to be studied (Pugsley and Berruty, 1996). Figure 1 presents actual lines and tendencies for the third type mathematical modelling. There are two major tendencies for modelling, following a treatment either Eulerian for both phases (Eulerian formulation) or Eulerian for the fluid phase and Lagrangian for the particulate phase (Eulerian-Lagrangian formulation).

In the Eulerian formulation both phases are treated as a continuum. Each phase is modelled separately in terms of a system of conservative equations for mass, momentum and energy. The conservative equations present terms accounting for interface interaction, which are related to mass, momentum and energy exchanges through the interface. In their traditional formulation, those models require the dynamic viscosity of the solid phase to be specified. A constant value for the solids viscosity can be obtained through momentum balances applied to experimental data (Miller and Gidaspow, 1992). The traditional Eulerian formulation has been extensively applied to fluidization (Gidaspow and Ettehadieh, 1983, Syamlal and Gidaspow, 1985, Gidaspow, 1986, Bouillard et al., 1989 among others). A new approach has been developed by a number of researchers for dealing with solids phase viscosity (Jenkins and Savage, 1983, Lun et al., 1984, Jenkins and Richman, 1985). This is the kinetic theory of granular flow (KTGF), which is based on the kinetic theory of dense gases (Chapman and Cowling, 1970). Bagnold (1954) apud Gidaspow (1994) is generally credited with starting the kinetic theory approach of granular flow (Gidaspow, 1994). The KTGF is based on an analogy between the flow of granular materials and the motion of gas molecules (Peirano, 1996). In spite of allowing the direct determination of the solids viscosity, the KTGF implies in more complex numerical procedures and higher computation times. There are two different main procedures for applying the KTGF, either considering (Ma and Ahmadi, 1988, Balzer et al., 1996) or not (Ding and Gidaspow, 1990) the presence of interstitial gas among the particles. Gidaspow (1994) was the first to present a detailed theoretical derivation of KTGF and to consider its application to particulate flows. Peirano and Leckner (1998) presented a detailed description concerning mathematical modelling of gas-solids turbulent flows in circulating fluidized beds using a KTGF. Among the research works applying the KTGF to fluidization modelling are Ding and Gidaspow (1990), Boemer et al. (1995), Kim and Arastoopour (1995) and Samuelsberg and Hjertager (1996).

In Eulerian-Lagrangean formulation only the gas phase is assumed as a continuum. This approach allows a better understanding of the particle-particle interactions. However, the Eulerian-Lagrangean formulation requires a complete set of equations (Newton’s second and third laws) to be written for each particle in the flow field, difficulting its application to fluidization in view of computational limitations. Otherwise it is a very useful tool for the development of new rheological models for fluidized suspensions, and to enhance the formulation of closing laws required by two-fluids Eulerian models. An explanation of general procedures regarding Eulerian-Lagrangian formulation is presented by Kuipers and van Swaaij (1997).

Theoretical Formulation of the Two-Fluids Model

The formal approach for the formulation of two-fluids Eulerian models for two-phase flows is well described in the literature (Ishii, 1975, Delhaye and Achard, 1976 and 1977, Delhaye, 1981, Drew, 1983) using different averaging procedures. Enwald et al. (1996) applied the procedure for the two-phases gas-solids flows as shown in Figure 2. Initially, mass and momentum integral balances are performed on a fixed control volume comprising both phases. The theorems of Leibniz and Gauss are then applied reducing the integral balances into two types of equations: the local instantaneous conservative equations for each phase, and the local instantaneous jump equations which account for interface interactions. Following, the Euler averaging procedure is applied to the local instantaneous equations, thereby eliminating the local instantaneous fluctuations characteristic of the flow field. As a result new variables arise on the conservative equations related to terms of interface interaction. Then, closing laws are applied in order to model unknown terms of the conservative equations. Finally, initial and boundary conditions are formulated.

Local Instantaneous Equations

The local instantaneous equations are the fundamental basis for multiphase models derived from averaging procedures. In fact, when the different phases separated by interfaces are assumed as continua, the local instantaneous formulation assumes a character of extreme mathematical rigorousness. Therefore, it is advisable that all models for multiphase flows should be derived from the local instantaneous formulation associated to some averaging procedure. Besides, the local instantaneous equations allow the direct modelling of separated flow processes (Ishii, 1975).

For obtaining the local instantaneous conservative equations and jump conditions, integral balances are performed for a variable y _k on an Eulerian control volume for both phases (k=1 for gas and k=2 for solids). The phases are separated by an interface of area A_I(t) which moves at a velocity . The theorems of Gauss and Leibniz are then applied, giving rise to the following equations:

Continuity local instantaneous equation:

Continuity jump condition at the interface:

Momentum local instantaneous equation:

Momentum jump condition at the interface:

where r k and are the density and velocity of phase k, is the normal unitary vector external to the interface of the volume occupied by phase k, is the stress tensor, is the gravity acceleration and is the mass flow rate at the interface, defined as:

Averaging Procedure and Average Equations

The averaging procedure eliminates terms of the local instantaneous formulations which result from fluctuations characteristic of multiphase flows, and the governing equations become expressed in terms of mean parameters and properties. So, only the macroscopic parameters of the flow are modelled, which are the features of major interest in industry (Ishii, 1975). Multiphase flow simulation is performed by solving the resulting average conservative equations together with closure laws and initial and boundary conditions.

The procedure followed to obtain the average equations comprise two steps. Firstly, the local instantaneous equations are multiplied by the phase density function X_K, a general average operator satisfying the axioms of Reynolds. In the second step the averages of products of dependent variables are transformed into products of averages by applying both Reynolds decomposition and weighted averaging. The phase density function is given by:

Three types of averaging operators are of common use, i.e. time, volume and ensemble averaging operators. The last is more representative for unsteady heterogeneous processes such as gas-solids flows. However, the evaluation of results obtained through statistical operators becomes almost impossible due to the nonrepetibility of experiments. Nevertheless, in multiphase flows it is common to assume the ergodicity hypothesis, which states that the three operators are equivalent for steady homogeneous flows. By enforcing this hypothesis its is possible to compare, for example, numerical results obtained through the application of ensemble average operator to time averaged experimental results. A detailed description of operators and averaging procedures can be found in Ishii (1975), Delhaye and Achard (1976, 1977) and Delhaye (1981).

When applying decomposition of Reynolds new terms are generated involving fluctuations. For a generical variable f it follows that:

where represents the average weighed value of f and f ' the fluctuation around this average. In general, the average weighed value of f is defined as:

where W is a weighing factor. There are two different weighing procedures, either through the phase density function, or through the product of the phase density function to the density of the phase. They are named, respectively, phase average and mass weighed average or Favre average, and are given by:

By applying the above described averaging procedures on the conservative local instantaneous equations, the average conservative equations result:

Average continuity equation:

Average momentum equation:

where represents the mean velocity of the phase k, and represents the fluctuation of the stress tensor for the phase k, . MkI accounts for the momentum transfer to the phase k at the interface, and is given by:

A detailed derivation of the above equations can be found in Enwald et al. (1996).

Closure Laws

The closure laws may be classified in three types:

1. Constitutive laws: account for physical properties of the phases by using experimental data and axioms. Constitutive laws are required for modelling the stress tensor, the dynamic viscosity, the volumetric viscosity and the pressure for all phases.

2. Transfer laws: describe different interface interactions between the phases through empirical correlations. Transfer laws are required for modelling mass, momentum and energy transfer between the phases at the interface.

3. Topological laws: describe the spatial distribution of a specific flow variable.

The literature presents a great variety of proposals for closure laws. In fact, a suitable set of closure laws should be specified taking into account the particular physical phenomena to be described and a mathematical model to be developed. Detailed discussions on this issue can be found in Arnold et al. (1990), Truesdell and Toupin (1960) and Aris (1962). The closure laws commonly adopted in the classic two-fluids model for gas-solids flows are presented next.

Constitutive Laws

Stress tensor:

Despite it is well known that the stress tensor for both phases should be a function of the void fraction, the spatial derivatives of and memory effects, there is still no available general formulation providing reliable values for materials constants. In fact, the studies related to the rheology of fluidized powders have not conducted yet to the proposition of an unified rheological model (Kuipers et al., 1992). In view of this the viscous stress tensor for both phases, , is modelled under the assumption of Newtonian fluid having in view the hypothesis of Stokes as considered in Aris (1962) and White (1992), i.e.:

The stress tensor is given by:

Dynamic viscosity:

There are two procedures for determining the dynamic viscosity of the solids phase, m_s (Gidaspow, 1994 and Kuipers and van Swaaij, 1997). Those are:

1. To assume ms as a constant, considering the solid phase as a Newtonian fluid.

2. To calculate ms from the kinetic theory of granular flows (KTGF).

In the first procedure ms is calculated from experimental data. Integrating the momentum equation for the mixture in the axial direction over the cross section, Tsuo (1989) and Tsuo and Gidaspow (1990) obtained:

where R_w is the internal radius of the circulating fluidized bed. Turbulence of the solids phase is disregarded. The calculation methodology as well as required experimental data are described in detail by Miller and Gidaspow (1992), who also presented an example of calculation and a table of compiled results.

The literature presents several formulations for the dynamic viscosity of the solids phase derived from the KTGF. All of these formulations are based on the composition of kinetic and collisional effects. Boemer et al. (1995) present a discussion on this subject.

Bulk viscosity:

The bulk viscosity, lk, is determined having in view the hypothesis of Stokes, as follows:

Similar to ms, the bulk viscosity for the solids phase can also be obtained form the KTGF. Boemer et al. (1995) present a discussion on this issue.

Gas pressure:

The gas pressure is the thermodynamic pressure, i.e.:

Solids pressure:

The physical concept of solids or particulate pressure, Ps, is difficult to define. Boemer et al. (1995) defined the pressure of the solids phase as the normal force per unit area exerted over the solids phase due to interactions among particles. Campbell and Wang (1991) defined Ps as the force per unit area exerted over the surface of the solids phase, reflecting the total transport of momentum which can be attributed to the motion of the particles and their interactions.

The literature reports two different general procedures to determine the pressure of the solids phase. Similarly to the dynamic viscosity of the solids phase, Ps can also be obtained by:

1. Applying the KTGF.

2. Applying the traditional procedure which is based on particle properties and the local void fraction a_g.

The first procedure is described in Gidaspow (1994), Boemer et al. (1995) and Enwald et al. (1996) among others.

The traditional procedure is followed in this work. It is fully developed in Gidaspow (1994), Enwald et al. (1996) and Cabezas- Gómez (1999) among others. According to this procedure, three different effects account for the solids pressure: the transport of momentum due to fluctuations on particles velocity, ; the collisions among particles, ; and the contribution of the gas phase pressure, (Enwald et al. 1996). Neglecting the effect due to the fluctuations on particles velocity, the gradient of the solids pressure results:

where G(a_g) is a coefficient of interaction among particles also known as the elasticity modulus for the particulate phase, which is obtained from experiment. Table 1 presents literature correlations for G(a_g). According to Gidaspow and Ettehadieh (1983) Equation (20) gives significant results when the void fraction is lower than its value at minimum fluidizing conditions. This correlation contributes for making more stable the numerical system of equations, since it turns into real the imaginary characteristic directions. According to Gidaspow (1986) some calculation is required to adjust G(a_g), as done in Syamlal (1985), in order to prevent the void fraction of reaching too low values causing particle concentrations higher than the maximum possible. In spite of the fact that Equations (20-27) give very different predictions for G(a_g), this may affect very little the solutions for the time averaged pressure, velocities and phase volumetric fractions (Enwald et al., 1996). According to Enwald et al. (1996) this is a subject for further research.

Transfer Laws

In this work only the interface momentum transfer, M_kI, is considered. It is modelled through empirical correlations, where the jump conditions impose the necessary restrictions. Enwald et al. (1996) present the formulation of M_kI in detail. Following the classical theory M_kI is considered to be approximately equal to the so called generalised drag force per unit volume exerted at the interface of a suspension with particles of mean diameter d_p, represented by . This force comprises various effects, but only the effect due to stationary drag at the interface is usually considered. The stationary drag force at the interface, , may be determined through two different procedures:

1. From correlations for the drag coefficient over a single particle in a suspension of particles, C_D.

2. From the pressure drop per unit length through a suspension of particles, DP/L.

The literature presents several models for determining using both the above procedures. Whatever the model considered, the aim is to express the drag force as a function of the dependent variables of the conservative equations. In general is defined in function of the relative velocity between the phases, , and the momentum exchange coefficient at the interface, b , also known as drag function. MkI is given by:

where n_p is the number of spherical particles per unit volume. The calculation of the drag function depends on the adopted procedure. Under the first procedure b results:

where may be either the relative velocity or the drift velocity , depending on which velocity C_D is related to. The literature available expressions for the drag coefficient are valid for single isolated particles, C_Ds. In order to account for this Equation (29) includes a function dependent on the gas volumetric fraction, f(a_g), so that the influence of other particles of the suspension on C_D is considered. The drift velocity is defined as:

Following the second procedure for determining the stationary drag force at the interface, the drag function results:

The above expression results from the application of the momentum conservation equation for the gas phase under stationary flow condition, disregarding acceleration and neglecting attrition against the walls and the gravity effect (Gidaspow and Ettehadieh, 1983).

Equations (29) and (31) produce partially general results, depending on the correlations used to calculate C_Ds, , f(a_g) and DP/L. Cabezas- Gómez (1999) presents a compilation of several correlations used in the literature for modelling b.

The Closed General PDE System for a Gas-Solids Flow

The closed general PDE system for modelling gas-solids flows considering the adopted classical formulation procedure results:

Continuity equation

Momentum equation

The above equations where obtained under the following general hypothesis: both phases are continua; laminar and isothermal flow for both phases; no mass transfer at the interface; the particulate phase is homogeneously characterized by an average diameter; unreactive flow.

Hydrodynamic Models A and B from IIT/ANL

The IIT/ANL group (Illinois Institute of Technology/Argonne National Laboratory) has developed modelling of gas-solids flows by applying the above developed procedures, which they called Model A. From this formulation they derived a second model which they called Model B. Table 2 presents the equations of both models. In Model A the pressure gradient arises in the momentum equations for both phases. In Model B, owing to some mathematical manipulations, the pressure gradient arises only in the momentum equation for the gas phase, and the expression for the drag function b is modified to account for both the principle of Archimedes and the usual relation of minimum fluidization given by Kunii and Levenspiel (1991).

Model B was proposed in view of the conditional numerical stability of Model A owing to the occurrence of characteristic directions which become complex (Bouillard et al., 1989). It is recognised that Model B represents the transformation of the conditionally stable equations of Model A into a well-posed set of equations for initial value unidimensional problems. However, there is no evidence in the literature to support that the well-posed character extends to multidimensional problems (Enwald et al., 1996). The procedure for deriving Model B can be found in Bouillard et al. (1989) and Gidaspow (1994).

In the equations of Table 2, concerning the gas (g) and solids (s) phases, and are mean velocities (m/s), r_g and r_s are densities (kg/m³), a_g and a_s are volumetric fractions, and and are viscous stress tensors (Pa). Also, P is the thermodynamic gas pressure (Pa), is the acceleration of gravity (m/s²), G is the modulus of elasticity of the solids (N/m²), b is the drag function at the interface (kg/m²-s), C_Ds is the drag coefficient for a single particle in an infinite medium, Re_p is the number of Reynolds based on the particle diameter, f_s is the sphericity of the particles, m stands for dynamic viscosity (kg/m-s), R_g is the ideal gas constant (kJ/kg-K) and t is the time (s).

Initial and Boundary Conditions

The initial and boundary conditions commonly used in the classical modelling approach are formulated below for both phases in a circulating fluidized bed. It is important to note that the boundary conditions at the walls for the solids phase are not trivial. Basically, velocities must be established at the walls. For the gas phase the non-slip condition is assumed, i.e.:

For the solids phase the condition defining the tangential velocity at the walls oscillates between non-slip and free-slip conditions. A possible assumption is that the fine particles stick to the walls while the coarse particles roll over its surface. According to Soo (1983) apud Tsuo and Gidaspow (1990) this assumption give rise to:

where Lp is the free mean path of the particles, which is given by:

Ding and Gidaspow (1990) followed a model proposed by Eldighity et al. (1977), and assumed that, since the particles diameter is usually higher than the characteristic roughness of the wall, the particles slide in a partial way. They obtained:

where l_p is a sliding parameter which represents the mean distance among the particles, and is given by:

In Equations (49) and (51) y represents the direction normal to the wall. The KTGF can also be used in order to evaluate the boundary conditions at the wall (Jhonson and Jakcson, 1987). A general discussion regarding initial and boundary conditions, including those at the entrance and at the exit of a circulating fluidized bed can be found in Tsuo (1989), Gidaspow (1994), Peirano and Lekcner (1998) and Cabezas- Gómez (1999).

Numerical Solution Methodology

The literature presents several computational codes used for the solution of classical two-phases gas-solids models (see Kuipers and van Swaaij, 1997). This work uses the MULTIFIX code developed by Syamlal (1985). It is an extension of the previous K-FIX code, initially developed for gas-liquid flows (Rivard and Torrey, 1977 apud Syamlal, 1985), and later on adapted to deal with gas-solids flows (Ettehadieh, 1982). K-FIX is based on a numerical method developed by Harlow and Amsden (1975), which is an extension of the implicit continuous-fluid Eulerian technique (ICE) developed by Harlow & Amsden (1971).

A modified version of MULTIFIX is used following Tsuo (1989). The VSSCG procedure is used for modelling the gas phase stress tensor, while the similar procedure VSSCS is used for the solids phase stress tensor (Ding, 1990). Equations (51) and (52) are used for the calculation of momentum boundary conditions for the solids phase at the walls.

Following the procedure, the hydrodynamic conservative equations are discretized in finite differences. The continuity equations are discretized implicitly, while the momentum equations are discretized according to an explicit-implicit technique. In the momentum equations, the terms related to gas pressure, solids pressure and interface momentum exchange are discretized implicitly, while the remaining terms are discretized explicitly. A fixed cell 2D computational grid is used. The scalar variables are set at the center of the cell while the vector variables are placed at the boundaries of the cell.

The discrete continuity equations for the gas and the solids phases result:

The discrete momentum equations for the gas and the solids phases result:

A convergence criterious is defined by writing the gas phase continuity equation as:

where is the residue of the mass conservation equation for the gas phase, which must be zero for absolute convergence. The subscripts r and z as well as i and j represent respectively the radial and the axial directions. The superscripts n and n+1 represent the values at the respective time levels t and t + Dt.

The overall iterative calculation procedure is as follows:

1. The convective explicit terms and are determined.

2. The drag coefficients are explicitly determined.

3. Velocities are calculated for the level n+1 from Equations (55) and (56).

4. The coefficient is calculated as:

where is the velocity of sound in air. The coefficient represents the derivatives of the residue related to the pressure, and is analytically obtained from the momentum conservation equations.

5. Pressure is then updated. This is performed through the domain in each computational cell until the convergence criterious is met, or until the number of iterations exceeds an inner maximum limit allowed. At the end of such computational sweep, if a pressure adjustment is necessary in any of the cells, the procedure is repeated until simultaneous convergence in all the cells is obtained. The number of iterations, however, is restricted by an outer iteration limit.

The iterative calculation procedure for pressure is as follows:

1. The residue is determined from Equation (57). If the convergence criterious for is met the calculation is finished, otherwise go to step 2.

2. The pressure is updated. The initial iteration for pressure is performed using the Newton point iteration method in each cell. This procedure is continued till changes sign. After this the next pressure correction is done using the secant method. Thereafter a combination of the secant method and the bisection method is used (Syamlal, 1985, Ding, 1990).

3. The velocities are calculated from Equations (55) and (56)

4. The mass fluxes of both phases are determined.

5. The volumetric fraction of solids is calculated from Equation (54)

6. a_g is found from the following relation:

7. Back to step 1.

Numerical Simulations

Tsuo (1989) developed mathematical modelling using the two-fluids model for studying heat transfer to surfaces, cluster formation and the structure of the annular flow through the riser in circulating fluidised beds. The author simulated some installations, including an IIT transport system which had already been studied before by Luo (1987). This particular IIT installation is also simulated in this work for studying the effects on the flow hydrodynamics of the viscosity of the solids phase. The mathematical modelling developed follows model B of IIT/ANL in a 2D cylindrical coordinate system.

The concerning IIT installation is described in Luo (1987). Figure 3 presents the initial and boundary conditions taken into account in this work, which are the same considered by Tsuo (1989). The column geometry was considered symmetrical, allowing the assumption of 2D cylindrical coordinate system and symmetry condition along the axis.

Figures 4 to 8 show the results of simulation. Figures 4 to 6 show the results for a fixed value of the solids phase viscosity ms = 0.509 kg/(ms), which was determined by Tsuo (1989) using Eq. (16) for the same conditions assumed in the present simulations. Figures 7 and 8 show the results for different values of m_s allowing an analysis of its effects on the flow hydrodynamics.

Results of Simulation for a Fixed ms

Figure 4 shows the time variation of the solids volumetric fraction through riser, where darker shades of grey mean lower solids volumetric fraction. The figure shows the formation and dissipation of clusters. The first cluster is formed at 1.5 seconds. The results are qualitatively similar to those obtained by Tsuo (1989) using a 2D Cartesian coordinate system for the same installation and conditions. Nevertheless, contrary to the results of Tsuo (1989), who did not assumed symmetry conditions, the solids accumulate mainly along the axis of symmetry. Otherwise, Neri and Gidaspow (2000) obtained the same effect obsereved in this work when simulating the IIT CFB, shown in Miller and Gidapow (1992), using the symmetry boundary condition at center of the riser. The axis of symmetry acts as a mirror, not allowing for any mass transfer from one wall to the other (Neri and Gidaspow, 2000).

Figure 5 shows the time variation of the solids phase mass velocity G_s in the riser 5.3 meters above the entrance. G_s is obtained by integrating the local solids phase mass velocity over the cross section of the riser. The graph shows the classic behaviour of circulating fluidized beds, including low frequency oscillations. Tsuo (1989) reported bed oscillations of about 0.2 Hertz. The results of simulation are in good agreement with the empirical data of Luo (1987). A tendency to a well established stationary regime is not observed, which is in agreement to practice.

Figure 6 shows radial profiles of the local axial velocities for both phases at 3.4 meters above the entrance. Simulated profiles averaged over different time intervals are compared with the experimental data of Luo (1987). For the gas phase the predictions present great deviations from experimental data. For both phases the profiles show an unreal behaviour close to the axis, which is possibly due to the adoption of a 2D cylindrical coordinate system and symmetry around the axis. This suggests that the symmetry boundary condition is not adequate for vertical gas-solids flows in circulating fluidized beds. It is seen that despite different time intervals were considered for averaging the velocity profiles, the results show similar tendencies. This must be emphasized since the flow under consideration is characterized by great fluctuations of the velocity fields, and in fact never reaches the condition of stationary flow.

Results of Simulation for Various ms

Results of simulation are generated for various ms keeping all the remaining conditions as defined in Figure 3. Values of m s were arbitrarily chosen around the experimental value, 0.509 kg/(ms), established by Tsuo (1989).

Figure 7 shows the time variation of the solids phase mass velocity Gs in the riser 5.3 meters above the entrance, for various solids phase viscosities. It is seem that the solids phase viscosity has little effect over Gs. In fact this behaviour was expected in view of the continuity enforcement. The predictions are in good agreement with the empirical data of Luo (1987) no matter the value of ms. The low frequency oscillations characteristic of circulating beds are observed. No tendency is observed towards a well established stationary regime, which was expected owing to the unstable nature of the flow.

Despite not affecting significantly the solids phase mass velocity, the solids phase viscosity considerably affects velocity profiles of both phases. Figure 8 shows radial profiles of the local axial velocities for both phases 3.4 meters above the entrance, for various solids phase viscosities, in comparison to experiment (Luo, 1987). The predicted profiles were averaged over the time interval from 10 to 48 seconds for all solids phase viscosities, except for m_s = 0.509 kg/(ms) where the time interval was from 10 to 40 seconds (In this case the procedure diverged for times greater than 42 seconds). The comparison of the profiles for different m_s shows relative deviations of up to 60% for both phases. As m_s increases the velocity profiles become increasingly flat, showing a cross sectional diffusion effect of m_s. This effect does not seem to be a consequence of numerical diffusion since the spreading of the profile for m_s = 0.75 kg/(ms) is less than expected (it would be flatter in comparison with the profile for m_s = 0.509 kg/(ms) if numerical diffusion was the reason). The results are particularly unsatisfactory for m_s = 0.25 kg/(ms). In this case the influence of the symmetry boundary condition at the axis seems to become more intense.

It is clear from Figure 8 that the solids phase viscosity considerably affects the hydrodynamic predictions. Therefore, a more precise theoretical methodology is required for accurately determining this parameter, such as the kinetic theory of the granular flow.

General Remarks

A number of aspects related to gas-solids two-fluids modelling require further development if accuracy is to be achieved. Assumed hypothesis such as spherical particles, particulate phase homogeneously characterized by an average diameter, and treated as Newtonian continuum fluid are certainly cause of inaccurate predictions. Besides, most of the information required for establishing reliable closure laws are based upon both theories which need enhancement and experimental data which are scattered and limited in view of experimental techniques shortcomings.

The solids phase pressure is modelled through empirical correlations which give predictions quite different from each other. This seems to be a consequence of the fact that the parameter is introduced in the conservative equations for numerical stability reasons, as shown by Gidaspow (1986). The interface momentum transfer is modelled considering only the effect of stationary drag. According to Enwald et al. (1996) the literature proposed correlations for the drag function must be used with caution since most of them are derived from experiments in homogeneous liquid-solids fluidized beds.

Turbulence effects may be accounted for in two-fluids modelling through Reynolds stress terms. According to Enwald et al. (1996), in bubbling fluidized beds which are characterised by high solids concentrations, the particulate inertia attenuates the gas phase turbulence. However, for circulating fluidized beds which are mostly characterised by regions of low solids concentrations, turbulence becomes significant.

The prescription of entrance boundary conditions for both phases and wall boundary conditions for the solids phase are not straightforward (Enwald et al., 1996). One-dimensional plug flow is prescribed for both phases at the column entrance even though it is known that the entrance flow conditions are in fact dependent upon the particular distributor plate project. The assumed geometry and symmetry are also sources of inaccuracy. Gas-solids flows in circulating fluidized beds are characterized by three-dimensional effects due to non-uniform shapes and sizes distributions of the particles, asymmetric solids feeding, asymmetric geometry at the exit of the riser, and the presence of solids separators such as cyclones when the whole installation is modelled.

Conclusions

The predicted time variation of the solids volumetric fraction through the riser of circulating fluidized beds are indicative of a quite complex hydrodynamic behaviour. There are steep variations on both radial and axial solids concentration profiles, together with frequent formation and dissociation of clusters flowing both upwards and downwards.

It is observed a very significant influence of the solids phase viscosity on the predictions. It is clear that this parameter should be rigorously modelled through a suitable theory which can accurately characterize the fluids interactions and the physical phenomena involved. The kinetic theory of the granular flow is a promising theory for this purpose.

The predicted time variation of the solids phase mass velocity is in good agreement with experiment. However, the predicted time averaged radial profiles of the local axial velocities for both phases greatly deviate from experiment. This is possibly due to model assumptions and boundary conditions (specifically the symmetry at the center of the column), which impose conditions quite different from practice. On the other hand, the uncertainty on the experimental measurements which are characteristic of gas-solids multiphase flows must be taken into account.

Results of simulation from literature and presented here strongly suggest that many features related to two-fluids Eulerian mathematical modelling applied to circulating fluidized beds need improvement. Special attention should be directed towards interface momentum transfer, turbulence in both phases, pressure and viscosity of the solids phase, and realistic boundary conditions. On the numerical side it becomes evident that more refined grids are required which unfortunately becomes prohibitive in view of the current stage of development of computers.

Acknowledgments

The authors wish to acknowledge the support of FAPESP (process 98/13812-1) and CNPq for this work.

Harlow, F. H. and Amsden, A. A., 1971, "A numerical fluid dynamics calculation method for all flow speeds", Journal of Computational Physics, v.8, p.197-213.

Harlow, F.H. and Amsden, A. A., 1975, "Numerical calculations of multiphase fluid flow", Journal of Computational Physics, v.17, p.19-52.

Peirano, E. and Leckner, B., 1998, "Fundamentals of turbulent gas-solid flows applied to circulating fluidized bed combustion", Progress in Energy and Combustion Science, v.24, n.4, p.259-296.

Manuscript received: June 2000. Technical Editor: Aristeu S. Neto.

• Aris, R., 1962, "Vector, tensors, and the basic equations of fluid mechanics", Dover Publications Inc., New York.
• Arnold, G.S., 1990, "An assessment of multiphase flow models using the second law of thermodynamics", International Journal of Multiphase Flow, v.16, n.3, p.481-494.
• Balzer, G., Simonin, O., Boëlle, A. and Lavieville, J., 1996, "A unified modelling approach for the numerical prediction of dilute and dense gas-solid two-phase flow", In: Kwauk, M., Li, J., eds. Circulating Fluidized Bed Technology IV, Science Press, Beijing, p. 432-439.
• Boemer, A., Qi, H., Renz, U., Vasquez, S. and Boysan, F., 1995, "Eulerian computation of fluidized hydrodynamics – A comparison of physical models", Proc.13th Int. Conf. on Fluidized Bed Combustion, v.2, p.775-787.
• Bouillard, J.X, 1986, "Hydrodynamics of sedimentation, fluidization, and erosion", Ph.D. Thesis, Illinois Institute of Technology, Chicago.
• Bouillard, J.X., 1989, "Porosity distributions in a fluidized bed with an immersed obstacle", AIChE Journal, v.35, n.6, p.908-922.
• Cabezas- Gómez, L., 1999, "An study of mathematical modelling and numerical simulation of gas-solids two-phase flows", MS Dissertation, Escola de Engenharia de São Carlos, USP, Brasil.
• Campbell, C.S. and Wang, D.G., 1991, "Particle pressures in gas-fluidized beds", Journal of Fluid Mechanics, v.227, p.495-508.
• Chapman, S. and Cowling, T.G., 1970, "The mathematical theory of non-uniform gases", 3rd ed., Cambrige Univ. Press, Cambrige, U.K.
• Delhaye, J.M. and Achard, J.L., 1976, "On the averaging operators introduced in two-phase flow modeling", Specialists’ Meeting on Transient Two-Phase Flow, August 3-4, Toronto.
• Delhaye, J.M. and Achard, J.L., 1977, "On the use of averaging operators in two-phase flow modeling", Thermal and hidraulic aspects of nuclear reactor safety, 1: Light water reactors, ASME winter meeting.
• Delhaye, J.M., 1981, "Thermohydraulics of two-phase systems for industrial design and nuclear engineering", ed. J. M. Delhaye.
• Ding, J., 1990, "A fluidization model using kinetic theory of granular flow", Ph.D. Thesis, Illinois Institute of Technology, Chicago.
• Ding, J. and Gidaspow, D., 1990, "A bubbling model using kinetic theory of granular flow", AIChE Journal, v.36, n.4, p.523-538.
• Drew, D.A., 1983, "Mathematical modeling of two-phase flow", Annual Review of Fluid Mechanics, v.15, p.261-291.
• Eldighidy, S.M., Chen, R.Y. and Comparin, R.A., 1977, "Deposition of suspensions in the entrance of a channel", Journal of Fluids Engineering – Transactions of the ASME, v., p. 365-370.
• Enwald, H., Peirano, E., and Almstedt, A.-E., 1996, "Eulerian two-phase flow theory applied to fluidization", International Journal of Multiphase Flow, v.22, p.21-66, Supplement.
• Ergun, S., 1952, "Fluid flow through packed columns", Chemical Engineering Progress, v.48, n.2, p.89-94.
• Ettehadieh, B., 1982, "Hydrodynamic analysis of gas-solids fluidized beds", Ph.D. Thesis, Illinois Institute of Technology, Chicago.
• Ettehadieh, B., Gidaspow, D. and Lyczkowski, R.W., 1984, "Hydrodynamics of fluidization in a semicircular bed with a jet", AIChE Journal, v.30, n.4, p.529-536.
• Gidaspow, D. and Ettehadieh, B., 1983, "Fluidization in two-dimensional beds with a jet. 2. Hydrodynamics modeling", I & EC Fundam., v.22, p.193-201.
• Gidaspow, D., 1986, "Hydrodynamics of fluidization and heat transfer: supercomputing modelling", Applied Mechanics Review, v.39, n.1, p.1-23.
• Gidaspow, D., Shih, Y.T., Bouillard, J. and Wasan, D., 1989, "Hydrodynamics of a lamella electrosettler", AIChE Journal, v.35, n.5, p.715-724.
• Gidaspow, D., 1994, "Multiphase flow and fluidization: continuum and kinetic theory descriptions", Academic Press Inc., Bostom.
• Ishii, M., 1975, "Thermo-fluid Dynamic Theory of Two-phase Flow", Eyrolles, Paris.
• Jenkins, J.T. and Savage, S.B., 1983, "A theory for the rapid flow of identical, smooth, nearly elastic, spherical particles, Journal of Fluid Mechanics", v.130, p.187-102.
• Jenkins, J.T. and Richman, M.W., 1985, "Grad’s 13 moment system for a dense gas of inelastic spheres", Archive for Rational Mechanics and Analysis, v.87, n.4, p.355-377.
• Johnson, J.T. and Jakcson, R., 1987, "Frictional-collisional constitutive relations for granular materials with application to plane shearing", Journal of Fluid Mechanics, v.176, p.67-93.
• Kim, H.S. and Arastoopour, H., 1995, "Simulation of FCC particles flow behavior in a CFB using modified kinetic theory", Canadiam Journal of Chemical Engineering, v.73, p.603-611.
• Kuipers, J.A.M., van Duin, K.J., van Beckum, F.P.H. and van Swaaij, W.P.M., 1992, "A numerical model of gas fluidized beds", Chemical Engineering Science, v.47, n.8, p.1913-1924.
• Kuipers, J.A.M. and van Swaaij, W.P.M., 1997, "Application of computational fluid dynamics to chemical reaction engineering", Reviews in Chemical Engineering, v.13, n.3.
• Kunii, D. and Levenspiel, O., 1991, "Fluidization Engineering", 2nd Ed., Butterworth-Heineman, Boston.
• Lun, C.K.K., Savage, S.B., Jeffrey, D.J. and Chepurniy, N., 1984, "Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general flowfield, Journal of Fluid Mechanics, v.140, p.223-256.
• Luo, K. M., 1987, "Dilute, Dense-Phase and Maximum Solids-Gas Transport", Ph.D. Thesis, Illinois Institute of Technology, Chicago.
• Ma, D. and Ahmadi, G., 1988, "A kinetic model for rapid granular flows of nearly elastic particles including interstitial fluid effects", Powder Technology, v.56, p.191-207.
• Massoudi, M., Rajagopal, K.R., Ekmann, J.M. and Mathur, M.P., 1992, "Remarks on the modelling of fluidized systems", AIChE Journal, v.38, n.3, p.471-472.
• Miller, A. and Gidaspow, D., 1992, "Dense, vertical gas-solid flow in a pipe", AIChE Journal, v.38, n.11, p.1801-1815.
• Neri, A. and Gidaspow, D., 2000, "Riser hydrodynamics: Simulation using kinetic theory", AIChE Journal, v.46, n.1, p.52-67.
• Nieuwland, J.J., van Sint Annaland, M., Kuipers, J.A.M. and van Swaaij, J.P.M., 1996, "Hydrodynamics modelling of gas/particle flows in riser reactors", AIChE Journal, v.42, n.6, p.1569-1582.
• Peirano, E., 1996, "The Eulerian/Eulerian formulation applied to gas-solid flows", Report A96-218, Department of Energy Conversion, Chalmers University of Technology, Sweden.
• Pugsley, T.S. and Berruti, F., 1996, "A predictive model for circulating fluidized bed risers", Powder Technology, v.89, p.57-69.
• Rowe, P.N., 1961, "Drag forces in a hydraulic model of a fluidised bed – part II", Transactions of the Institut of Chemicals Engineers, v.39, p.175-180.
• Samuelsberg, A. and Hjertager, B.H., 1996, "An experimental and numerical study of flow patterns in a circulating fluidized bed reactor", International Journal of Multiphase Flow, v.22, n.3, p.575-591.
• Shih, Y.T., Gidaspow, D. and Wasan, D.T., 1987, "Hydrodynamics of electro-fluidization: separation of pyrites from coal", AIChE Journal, v.33, p.1322-1333.
• Syamlal, M., 1985, "Multiphase hydrodynamics of gas-solids flow", Ph.D. Thesis, Illinois Institute of Technology, Chicago.
• Syamlal, M. and Gidaspow, D., 1985, "Hydrodynamics of fluidization: prediction of wall to bed heat transfer coefficients", AIChE Journal, v.31, n.1, p.127-135.
• Syamlal, M. and O’Brien, T.J., 1988, "Simulation of granular layer inversion in liquid fluidized beds", International Journal of Multiphase Flow, v.14, n.4, p.473-481.
• Truesdel, C. and Toupin, R., 1960, "Classical field theories", In: Handbuch der Physik, Springer, Berlin.
• Tsuo, Y.P., 1989, "Computation of flow regimes in circulating fluidized beds", Ph.D. Thesis, Illinois Institute of Technology, Chicago.
• Tsuo, Y.P. and Gidaspow, D., 1990, "Computation of flow patterns in circulating fluidized beds", AIChE Journal, v.36, n.6, p.885-896.
• White, F.M., 1991, "Viscous fluid flow", 2nd ed. McGraw-Hill Inc., NY.

Publication Dates

• Publication in this collection
  24 Oct 2001
• Date of issue
  2001

History

• Received
  June 2000