Hydrodynamic and Heat Transfer Simulation of Fluidized Bed Using CFD

1. Introduction

The nuclear energy is suffering from the lack of public acceptance everywhere mainly due to the issues relating to reactor safety, economy and nuclear waste. The Fluidized Bed Nuclear Reactor (FBNR) concept has addressed these issues and tried to resolve such problems. The FBNR is small, modular and simple in design contributing to the economy of the reactor. It has inherent safety and passive cooling characteristics. Its spent fuel being small spherical elements may not be considered nuclear waste, and can be directly used as a source of radiation for applications in industry and agriculture resulting in reduced environmental impact [1].

With the increase of computational power, numerical simulation becomes an additional tool to predict the fluid dynamics and the heat transfer mechanism in multiphase flow. A numerical hydrodynamic and heat transfer model has been developed to simulate the gas fluidized bed. All of CFD, in one form or another, is based on the fundamental governing equations of fluid dynamics (continuity, momentum and energy equations). These equations speak physics. They are mathematical statements of three fundamental physical principals upon which all fluid dynamics is based: mass is conserved, Newton’s second law and energy is conserved [2].

This chapter aims to study a mathematical modeling and numerical simulation of the hydrodynamics and heat transfer processes in a two-dimensional gas fluidized bed with a vertical uniform gas velocity at the inlet. The velocity, volume fraction, temperature distribution for gas phase and particle phase are calculated. Also, gas pressure and a prediction of the average heat transfer coefficient are also studied.

Such a simulation technique allows performance evaluation for different bed input parameters, and can evolve into a tool for optimized design of fluidized beds for different industrial use.

The numerical setup consists of a two dimensional fluidized bed filled with particles.The cold gas enters to the bed to cool the hot particles. Based on conservation equations for both phases it is possible to predict particles and gas volume fractions, velocity distributions (for gas and particles), temperature distribution, heat transfer coefficient as well as gas pressure field.

2. General assumptions for the mathematical model

Fluidized beds are categorized as multiphase flow problems. There are currently two approaches to model multiphase flow problems as discussed in chapter two. The best overall balance between computational time and accuracy seems to be achieved by implementing an Eulerian-Eulerian approach. The following assumptions are introduced into the present analysis:

1. The bed is two-dimensional.

2. Eulerian-Eulerian approach is applied.

3. The gas has constant physical properties.

4. Uniform fluidization.

5. Constant input fluid flux.

6. There is no mass transfer or chemical reaction between the two phases.

7. All particles are spherical in shape with the same diameter,d_p.

8. The expanded bed region is considered in the analysis in addition to the out flowing gas, i.e. suspension and free board regions.

9. Viscous heat dissipation in the energy equation is negligible in comparison with conduction and convection.

3. Governing equations

Due to the high particle concentration in fluidized beds the particles interactions cannot be neglected. In fact the solid phase has similar properties as continuous fluid. Therefore, the Eulerian approach is an efficient method for the numerical simulation of fluidized beds.

A hydrodynamic and thermal model for the fluidized bed is developed based on schematic diagram shown in Figure (1). The principles of conservation of mass, momentum, and energy are used in the hydrodynamic and thermal models of fluidization. The general mass conservation equations and the separate phase momentum equations and energy equations (for each phase) for fluid–solids, nonreactive transient and two-phase flow will be discussed in the following sections.

Figure 1.

The schematic diagram of the present work model

3.1. Continuity equation

Particle phase:

$$\frac{\partial {\epsilon }_s}{\partial t}+\frac{\partial }{\partial x}({\epsilon }_s{u}_s)+\frac{\partial }{\partial y}({\epsilon }_s{v}_s)=0$$E1

Gas phase :

$$\frac{\partial {\epsilon }_g}{\partial t}+\frac{\partial }{\partial x}({\epsilon }_g{u}_g)+\frac{\partial }{\partial y}({\epsilon }_g{v}_g)=0$$E2

3.2. Volume fraction constraint

$${\epsilon }_g^{}\text{}+{\epsilon }_s^{}=1.0$$E3

3.3. Momentum equation

Particle phase:

$${\rho }_s{\epsilon }_s\frac{\partial u_s}{\partial t}+{\rho }_s{\epsilon }_s\frac{\partial }{\partial x}(u_s{u}_s)+{\rho }_s{\epsilon }_s\frac{\partial }{\partial y}(v_s{u}_s)=F_{sx}\text{“x -direction “}$$E4

$${\rho }_s{\epsilon }_s\frac{\partial v_s}{\partial t}+{\rho }_s{\epsilon }_s\frac{\partial }{\partial x}(u_s{v}_s)+{\rho }_s{\epsilon }_s\frac{\partial }{\partial y}(v_s{v}_s)=F_{sy}\text{“y -direction “}$$E5

The total force acting on particle phase is the sum of the net primary force and the force resulting from particle phase elasticity. The x and y components of forces acting on particle phase are as following:

$$\begin{array}{l}{F}_{sy}=C_d\frac{3{\epsilon }_s{\rho }_g(v_g-v_s)\left|v_g-v_s\right|}{4d_p}{(}^1-{\epsilon }_s{\rho }_{s}g\\ -{\epsilon }_s\frac{\partial p}{\partial y}-\frac{\partial {\epsilon }_s}{\partial y}(3.2g{d}_p{\epsilon }_s({\rho }_s-{\rho }_g))\end{array}$$E6

$$F_{sx}=C_d\frac{3{\epsilon }_s{\rho }_g(u_g-u_s)\left|u_g-u_s\right|}{4d_p}{(1-{\epsilon }_s)}^{-1.8}-{\epsilon }_s\frac{\partial p}{\partial x}$$E7

Gas phase :

$${\rho }_g{\epsilon }_g\frac{\partial u_g}{\partial t}+{\rho }_g{\epsilon }_g\frac{\partial }{\partial x}(u_g{u}_g)+{\rho }_g{\epsilon }_g\frac{\partial }{\partial y}(v_g{u}_g)=F_{gx}\text{“x -direction “}$$E8

$${\rho }_g{\epsilon }_g\frac{\partial v_g}{\partial t}+{\rho }_g{\epsilon }_g\frac{\partial }{\partial x}(u_g{v}_g)+{\rho }_g{\epsilon }_g\frac{\partial }{\partial y}(v_g{v}_g)=F_{gy}\text{“y -direction “}$$E9

The fluid phase forces are readily obtained from the particle phase relations for fluid – particle interaction (drag and pressure gradient force), which act in the opposite direction on the fluid, together with gravity. The x and y components of forces acting on gas phase are as following:

$$F_{gy}=-C_d\frac{3{\epsilon }_s{\rho }_g(v_g-v_s)\left|v_g-v_s\right|}{4d_p}{(1-{\epsilon }_s)}^{-1.8}-{\epsilon }_g{\rho }_{s}g-{\epsilon }_g\frac{\partial p}{\partial y}$$E10

$$F_{gx}=-F_{sx}=-C_d\frac{3{\epsilon }_s{\rho }_g(u_g-u_s)\left|u_g-u_s\right|}{4d_p}{(1-{\epsilon }_s)}^{-1.8}+{\epsilon }_s\frac{\partial p}{\partial x}$$E11

3.3.1. Relation between fluid and particle velocities

We assume that both particles and fluid are regarded as being incompressible. This was justified on the basis that only a gas phase is going to exhibit any significant compressibility, and the orders of magnitude differences in particle and fluid density for gas fluidization render quite insignificant the small change in gas density resulting from compression. This assumption led to the relation linking fluid and particle phase velocities at all location. By applying the overall mass balance, which is obtained by summing equations (1) and (2):

$$\frac{\partial }{\partial x}({\epsilon }_s{u}_s+{\epsilon }_g{u}_g)+\frac{\partial }{\partial y}({\epsilon }_s{v}_s+{\epsilon }_g{v}_g)=0$$E12

Equation (12) shows that the total flux (fluid plus particles) in x- direction and y-direction remains constant, equal to that of fluid entering the bed U_gin,V_gin. This result is a simple consequence of the particles and fluid being considered incompressible :

$$V_{gin}={\epsilon }_g{v}_g+{\epsilon }_s{v}_s$$E13

$$U_{gin}={\epsilon }_g{u}_g+{\epsilon }_s{u}_s$$E14

Equations (13) and (14) enable the fluid velocity variables to be expressed in terms of the particle velocity at all points in the bed.

3.3.2. Combined momentum equation

In this section the combined momentum equation is produced by combining the fluid and particle momentum equations (4), (5), (8) and (9) by elimination of the fluid pressure gradient, which appears in them. This yields the combined momentum equation:

x- direction:

[divide equation (4-4) by ε_s ] + [divide equation (4-8) by ε_g] which give us :

$${\rho }_s\left[\frac{\partial u_s}{\partial t}+\frac{\partial }{\partial x}(u_s{u}_s)+\frac{\partial }{\partial y}(v_s{u}_s)\right]+{\rho }_g\left[\frac{\partial u_s}{\partial t}+\frac{\partial }{\partial x}(u_g{u}_g)+\frac{\partial }{\partial y}(v_g{u}_g)\right]=0$$E15

y- direction:

[divide equation (4-3) by ε_s ] – [divide equation (4-7) by ε_g] which give us :

$${\rho }_s\left[\frac{\partial v_s}{\partial t}+\frac{\partial }{\partial x}(u_s{v}_s)+\frac{\partial }{\partial y}(v_s{v}_s)\right]$$

$$-{\rho }_g\left[\frac{\partial v_g}{\partial t}+\frac{\partial }{\partial x}(u_g{v}_g)+\frac{\partial }{\partial y}(v_g{v}_g)\right]=\frac{F_{sy}}{{\epsilon }_s}-\frac{F_{gy}}{{\epsilon }_g}$$E16

Equations (15) and (16) with the continuity equation for the two phase (1) and (2),now define the two phase system without account of fluid pressure variation.

3.3.3. Drag coefficient

An important constitutive relation in any multiphase flow model is the formula for the fluid-particle drag coefficient, which is may be expressed by the empirical Dallavalle relation as reported in [3] :

$$C_d={\left(0.63+\frac{4.8}{{\mathrm{Re}}_{}^{0.5}}\right)}^2$$E17

The particle Reynolds number, Re, based on particle diameter is given by :

$$\mathrm{Re}=\frac{{\epsilon }_g{\rho }_g{d}_p\left|\overrightarrow{U_r}\right|}{{\mu }_g}$$E18

3.3.4. Gas pressure drop

Figure (2) shows the relation between the total pressure drop across the bed "ΔP_B" and the input gas velocity [3].where:

$$\Delta P_B=\left({\rho }_g{\epsilon }_{g,mf}+{\rho }_s\left(1-{\epsilon }_{g,mf}\right)\right)g{H}_{mf}$$E19

The gas pressure at the entrance of the fluidized bed can be calculated from the following equation:

$$P_{g,in}=P_{atm}+{\rho }_{g}g\left(H-H_{mf}\right)+\left({\rho }_g{\epsilon }_{g,mf}+{\rho }_s\left(1-{\epsilon }_{g,mf}\right)\right)g{H}_{mf}$$E20

and the pressure drop at any position and head "h" can be calculated from:

$$\Delta P_g={\rho }_{sus}gh$$E21

Where $${\rho }_{sus}$$is the suspension density which calculated from:

$${\rho }_{sus}={\rho }_s{\epsilon }_s+{\rho }_g{\epsilon }_g$$E22

Figure 2.

Total pressure drop in fluidized bed

3.4. Energy equation

Particle phase:

$$\begin{array}{l}{\rho }_s{C}_{ps}\frac{\partial {\epsilon }_s{T}_s}{\partial t}+{\rho }_s{C}_{ps}\frac{\partial }{\partial x}({\epsilon }_s{u}_s{T}_s)+{\rho }_s{C}_{ps}\frac{\partial }{\partial y}({\epsilon }_s{v}_s{T}_s)=\frac{\partial }{\partial x}({\epsilon }_s{k}_s\frac{\partial T_s}{\partial x})\\ +\frac{\partial }{\partial y}({\epsilon }_s{k}_s\frac{\partial T_s}{\partial y})+h_v(T_g-T_s\text{})+{\epsilon }_s\stackrel{•}{q}\end{array}$$E23

GAS phase:

$$\begin{array}{l}{\rho }_g{C}_{pg}\frac{\partial {\epsilon }_g{T}_g}{\partial t}+{\rho }_g{C}_{pg}\frac{\partial }{\partial x}({\epsilon }_g{u}_g{T}_g)+{\rho }_g{C}_{pg}\frac{\partial }{\partial y}({\epsilon }_g{v}_g{T}_g)=\frac{\partial }{\partial x}({\epsilon }_g{k}_g\frac{\partial T_g}{\partial x})\\ +\frac{\partial }{\partial y}({\epsilon }_g{k}_g\frac{\partial T_g}{\partial y})+h_v(T_s-T_g\text{})\end{array}$$E24

3.4.1. Thermal conductivity values (k_g and k_s)

The thermal conductivities of the fluid phase and the solid phase (k_g and k_s) in the two fluid model formulation should be interpreted as effective transport coefficients which means that the corresponding microscopic (or absolute) coefficients k_g,o and k_s,o cannot be used. It can be represented in general as:

$$k_g=k_g(k_{g,o},k_{s,o},{\epsilon }_g,\text{ particle geometry)}$$E25

$$k_s=k_s(k_{g,o},k_{s,o},{\epsilon }_g,\text{ particle geometry)}$$E26

However, such a general formulation is not yet available for fluidized beds and approximate constitutive equations have to be used. These approximate equations have been obtained on modeling of the effective thermal conductivity k_b in packed beds. According to their model, the radial bed conductivity k_b consists of a contribution k_b,g due to the fluid phase only and a contribution due to a combination of the fluid phase and the solid phase[4].

$$k_b=k_{b,g}+k_{b,s}$$E27

where :

$$k_{bg}=\left(1-\sqrt[]{1-{\epsilon }_g}\right)k_{go}$$E28

$$k_{bs}=\sqrt{1-{\epsilon }_g}\left(\omega A+\left(1-\omega \right)\Gamma \right)k_{go}$$E29

$$\Gamma =\frac{2}{\left(1-\frac{B}{A}\right)}\left[\frac{A-1}{{\left(1-\frac{B}{A}\right)}^2}\frac{B}{A}Ln\left(\frac{A}{B}\right)-\frac{B-1}{\left(1-\frac{B}{A}\right)}-\frac{1}{2}\left(B+1\right)\right]$$E30

$$B=1.25{\left(\frac{1-{\epsilon }_g}{{\epsilon }_g}\right)}^{\frac{10}{9}}$$E31

$$A=\frac{k_{s,o}}{k_{g,o}}$$E32

$$\omega =7.26X{10}^{-3}$$E33

Thus the thermal conductivities of the fluid phase and the solid phase then are given by :

$$k_g=\frac{k_{b,g}}{{\epsilon }_g}$$E34

$$k_s=\frac{k_{b,s}}{{\epsilon }_s}$$E35

3.4.2. Interphase volumetric heat transfer coefficient h_v

The heat transfer coefficient is modelled using a correlation by Gunn as reported [5]. This correlation is applicable for gas voidage in the range of 0.35 to 1 and for Reynolds numbers up to Re = 10⁵, and gives the Nusselt number:

$$Nu=\frac{h_{gp}{d}_p}{k_{g,o}}$$E36

$$Nu=\left(7-10{\epsilon }_g+5{\epsilon }_g^2\right)\left(1+0.7{\mathrm{Re}}_p^{0.2}{\mathrm{Pr}}^{\frac{1}{3}}\right)$$

$$+\left(1.33-2.40{\epsilon }_g+1.20{\epsilon }_g^2\right){\mathrm{Re}}_p^{0.7}{\mathrm{Pr}}^{\frac{1}{3}}$$E37

where Reynolds number is defined by equation (18).The Prandlt number is defined by;

$$\mathrm{Pr}=\frac{C_{p,g}{\mu }_g}{k_{g,o}}$$E38

,and the overall heat transfer coefficient is evaluated from;

$$h_v=\frac{6\left(1-{\epsilon }_g\right)h_{gp}}{d_p}$$E39

4. Boundary and initial conditions

The system of conservation equations (1),(2), (3), (4), (5), (8), (9), (23), and (24), nine equations which are discussed in previous sections must be solved for the nine dependent variables: the gas-phase volume fraction ε_g, the particle-phase volume fraction ε_s, the gas pressure P_g, the gas velocity components u_g and v_g and the solids velocity components u_s and v_s in x-direction and y-direction, respectively,the gas temperature T_g and particle temperature T_s.We need appropriate boundary and initial conditions for the dependent variables listed above to solve the system of equations.

4.1. Boundary conditions

In this section the boundary conditions for the above governing equations, which relate to two dimensional fluidized bed with width "L" and height "H" to allow bed expansion typically i.e. the height of the bed is enough to prevent the particles being ejected from the bed. Boundary conditions are imposed as follow:

$$\begin{array}{l}x=0:v_g=v_s=u_g=u_s=0,\frac{\partial {\epsilon }_g}{\partial x}=\frac{\partial T_g}{\partial x}=\frac{\partial T_s}{\partial x}=0\\ x=\frac{L}{2}:\frac{\partial {\epsilon }_g}{\partial x}=\frac{\partial T_g}{\partial x}=\frac{\partial T_s}{\partial x}=\frac{\partial v_g}{\partial x}=\frac{\partial u_g}{\partial x}=\frac{\partial v_s}{\partial x}=\frac{\partial u_s}{\partial x}=0\left(symmetric\right)\\ y=0:v_g=V_{g,in},u_g=u_s=v_s=0,{\epsilon }_g=1,P_g=P_{g,in},T_g=T_{g,in},\frac{\partial T_s}{\partial y}=0\\ y=H:\frac{\partial v_g}{\partial y}=\frac{\partial u_g}{\partial y}=\frac{\partial T_g}{\partial y}=0,{\epsilon }_g=1,v_s=u_s=0,P_g=P_{atm}\end{array}$$

4.2. Initial conditions

For setting the initial conditions, the model is divided into two regions: the bed and the freeboard. For each of the regions specified above, an initial condition is specified.

Bed region

$${\epsilon }_g={\epsilon }_{g,mf},v_g=\frac{V_{g,mf}}{{\epsilon }_{g,mf}},u_g=v_s=u_s=0,T_g=T_{g,in},T_s=T_{s,in}$$

Freeboard region

$${\epsilon }_g=1,v_g=V_{g,mf},u_g=v_s=u_s=0,T_g=T_{g,in}$$

5. Finite difference approximation scheme

The conservation equations are transformed into difference equations by using a finite difference scheme.

5.1. Discretization of continuity equations

5.1.1. Particle phase continuity equation

The particle phase continuity equation,Eq.(1), is discretized at the node i; j in an explicit form as:

$$\begin{array}{l}{\left({\epsilon }_s\right)}_{i,j}^{n+1}={\left({\epsilon }_s\right)}_{i,j}^n-\frac{\Delta t}{\Delta x}{\left(u_s\right)}_{i,j}^n\{\begin{array}{c}{\left({\epsilon }_s\right)}_{i,j}^n-{\left({\epsilon }_s\right)}_{i-1,j}^n,if{\left(u_s\right)}_{i,j}^n\ge 0.0\\ {\left({\epsilon }_s\right)}_{i+1,j}^n-{\left({\epsilon }_s\right)}_{i,j}^n,if{\left(u_s\right)}_{i,j}^n<0.0\end{array}\\ -\frac{\Delta t}{\Delta y}{\left(v_s\right)}_{i,j}^n\{\begin{array}{c}{\left({\epsilon }_s\right)}_{i,j}^n-{\left({\epsilon }_s\right)}_{i,j-1}^n,if{\left(v_s\right)}_{i,j}^n\ge 0.0\\ {\left({\epsilon }_s\right)}_{i,j+1}^n-{\left({\epsilon }_s\right)}_{i,j}^n,if{\left(v_s\right)}_{i,j}^n<0.0\end{array}\\ \end{array}$$E40

The gas phase volume fraction is then calculated explicitly as:

$${\left({\epsilon }_g\right)}_{i,j}^{n+1}=1-{\left({\epsilon }_s\right)}_{i,j}^{n+1}$$E41

5.1.2. Gas phase continuity equation

The gas continuity equation residual, d_g, is discretized at i; j in a fully implicit way:

$$\begin{array}{l}{d}_g={\left({\epsilon }_g\right)}_{i,j}^{n+1}-{\left({\epsilon }_g\right)}_{i,j}^n+\frac{\Delta t}{\Delta x}{\left(u_g\right)}_{i,j}^{n+1}\{\begin{array}{c}{\left({\epsilon }_g\right)}_{i,j}^{n+1}-{\left({\epsilon }_g\right)}_{i-1,j}^{n+1},if{\left(u_g\right)}_{i,j}^{n+1}\ge 0.0\\ {\left({\epsilon }_g\right)}_{i+1,j}^{n+1}-{\left({\epsilon }_g\right)}_{i,j}^{n+1},if{\left(u_g\right)}_{i,j}^{n+1}<0.0\end{array}\\ +\frac{\Delta t}{\Delta y}{\left(v_g\right)}_{i,j}^{n+1}\{\begin{array}{c}{\left({\epsilon }_g\right)}_{i,j}^{n+1}-{\left({\epsilon }_g\right)}_{i,j-1}^{n+1},if{\left(v_g\right)}_{i,j}^{n+1}\ge 0.0\\ {\left({\epsilon }_g\right)}_{i,j+1}^{n+1}-{\left({\epsilon }_g\right)}_{i,j}^{n+1},if{\left(v_g\right)}_{i,j}^{n+1}<0.0\end{array}\end{array}$$E42

$$$$

5.2. Discretization of combined momentum equations

The combined momentum equations may after a time discretization be expressed in the following forms for x and y directions:

x –direction:

$$\begin{array}{l}{\rho }_s\left[\frac{{\left(u_s\right)}_{i,j}^{n+1}-{\left(u_s\right)}_{i,j}^n}{\Delta t}+\frac{{\left(u_s\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left(u_s\right)}_{i,j}^n-{\left(u_s\right)}_{i-1,j}^n,if{\left(u_s\right)}_{i,j}^n\ge 0.0\\ {\left(u_s\right)}_{i+1,j}^n-{\left(u_s\right)}_{i,j}^n,if{\left(u_s\right)}_{i,j}^n<0.0\end{array}\right]\\ +{\rho }_s\left[\begin{array}{l}\frac{{\left(v_s\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left(u_s\right)}_{i,j}^n-{\left(u_s\right)}_{i,j-1}^n,if{\left(v_s\right)}_{i,j}^n\ge 0.0\\ {\left(u_s\right)}_{i,j+1}^n-{\left(u_s\right)}_{i,j}^n,if{\left(v_s\right)}_{i,j}^n<0.0\end{array}\\ \end{array}\right]\\ +{\rho }_g\left[\begin{array}{l}\frac{{\left(u_g\right)}_{i,j}^{n+1}-{\left(u_g\right)}_{i,j}^n}{\Delta t}+\frac{{\left(u_g\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left(u_g\right)}_{i,j}^n-{\left(u_g\right)}_{i-1,j}^n,if{\left(u_g\right)}_{i,j}^n\ge 0.0\\ {\left(u_g\right)}_{i+1,j}^n-{\left(u_g\right)}_{i,j}^n,if{\left(u_g\right)}_{i,j}^n<0.0\end{array}\\ \end{array}\right]\\ +{\rho }_g\left[\frac{{\left(v_g\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left(u_g\right)}_{i,j}^n-{\left(u_g\right)}_{i,j-1}^n,if{\left(v_g\right)}_{i,j}^n\ge 0.0\\ {\left(u_g\right)}_{i,j+1}^n-{\left(u_g\right)}_{i,j}^n,if{\left(v_g\right)}_{i,j}^n<0.0\end{array}\right]=0\end{array}$$E43

y –direction:

$$\begin{array}{l}{\rho }_s\left[\frac{{\left(v_s\right)}_{i,j}^{n+1}-{\left(v_s\right)}_{i,j}^n}{\Delta t}+\frac{{\left(u_s\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left(v_s\right)}_{i,j}^n-{\left(v_s\right)}_{i-1,j}^n,if{\left(u_s\right)}_{i,j}^n\ge 0.0\\ {\left(v_s\right)}_{i+1,j}^n-{\left(v_s\right)}_{i,j}^n,if{\left(u_s\right)}_{i,j}^n<0.0\end{array}\right]\\ +{\rho }_s\left[\begin{array}{l}\frac{{\left(v_s\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left(v_s\right)}_{i,j}^n-{\left(v_s\right)}_{i,j-1}^n,if{\left(v_s\right)}_{i,j}^n\ge 0.0\\ {\left(v_s\right)}_{i,j+1}^n-{\left(v_s\right)}_{i,j}^n,if{\left(v_s\right)}_{i,j}^n<0.0\end{array}\\ \end{array}\right]\\ -{\rho }_g\left[\frac{{\left(v_g\right)}_{i,j}^{n+1}-{\left(v_g\right)}_{i,j}^n}{\Delta t}+\frac{{\left(u_g\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left(v_g\right)}_{i,j}^n-{\left(v_g\right)}_{i-1,j}^n,if{\left(u_g\right)}_{i,j}^n\ge 0.0\\ {\left(v_g\right)}_{i+1,j}^n-{\left(v_g\right)}_{i,j}^n,if{\left(u_g\right)}_{i,j}^n<0.0\end{array}\right]\\ -{\rho }_g\left[\frac{{\left(v_g\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left(v_g\right)}_{i,j}^n-{\left(v_g\right)}_{i,j-1}^n,if{\left(v_g\right)}_{i,j}^n\ge 0.0\\ {\left(v_g\right)}_{i,j+1}^n-{\left(v_g\right)}_{i,j}^n,if{\left(v_g\right)}_{i,j}^n<0.0\end{array}\right]=F_y\end{array}$$E44

where:

$$\begin{array}{l}{F}_y={\left(C_d\right)}_{i,j}^n\frac{3{\rho }_g\left({\left(v_g\right)}_{i,j}^{n+1}-{\left(v_s\right)}_{i,j}^{n+1}\right)\left|{\left(v_g\right)}_{i,j}^n-{\left(v_s\right)}_{i,j}^n\right|}{4d_p}{\left(1-{\left({\epsilon }_s\right)}_{i,j}^n\right)}^{-1.8}\\ -\frac{{\left({\epsilon }_s\right)}_{i,j}^n-{\left({\epsilon }_s\right)}_{i,j-1}^n}{\Delta y}\left(3.2g{d}_p\left({\rho }_s-{\rho }_g\right)\right)\\ +{\left(C_d\right)}_{i,j}^n\frac{3{\rho }_g{\left({\epsilon }_s\right)}_{i,j}^n\left({\left(v_g\right)}_{i,j}^{n+1}-{\left(v_s\right)}_{i,j}^{n+1}\right)\left|{\left(v_g\right)}_{i,j}^n-{\left(v_s\right)}_{i,j}^n\right|}{4d_p}{\left(1-{\left({\epsilon }_s\right)}_{i,j}^n\right)}^{-2.8}\end{array}$$E45

5.3. Discretization of energy equations

Particle phase:

$$\begin{array}{l}{\rho }_s{C}_{p,s}\frac{{\left({\epsilon }_s{T}_s\right)}_{i,j}^{n+1}-{\left({\epsilon }_s{T}_s\right)}_{i,j}^n}{\Delta t}\\ +\frac{{\rho }_s{C}_{p,s}{\left(u_s\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left({\epsilon }_s{T}_s\right)}_{i,j}^n-{\left({\epsilon }_s{T}_s\right)}_{i-1,j}^n,if{\left(u_s\right)}_{i,j}^n\ge 0.0\\ {\left({\epsilon }_s{T}_s\right)}_{i+1,j}^n-{\left({\epsilon }_s{T}_s\right)}_{i,j}^n,if{\left(u_s\right)}_{i,j}^n<0.0\end{array}\\ +\frac{{\rho }_s{C}_{p,s}{\left(v_s\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left({\epsilon }_s{T}_s\right)}_{i,j}^n-{\left({\epsilon }_s{T}_s\right)}_{i,j-1}^n,if{\left(v_s\right)}_{i,j}^n\ge 0.0\\ {\left({\epsilon }_s{T}_s\right)}_{i,j+1}^n-{\left({\epsilon }_s{T}_s\right)}_{i,j}^n,if{\left(v_s\right)}_{i,j}^n<0.0\end{array}\\ =\frac{{\left({\epsilon }_s{K}_s{T}_s\right)}_{i+1,j}^n-2{\left({\epsilon }_s{K}_s{T}_s\right)}_{i,j}^n+{\left({\epsilon }_s{K}_s{T}_s\right)}_{i-1,j}^n}{\Delta x^2}\\ +\frac{{\left({\epsilon }_s{K}_s{T}_s\right)}_{i,j+1}^n-2{\left({\epsilon }_s{K}_s{T}_s\right)}_{i,j}^n+{\left({\epsilon }_s{K}_s{T}_s\right)}_{i,j-1}^n}{\Delta y^2}\\ +{\left(h_v\right)}_{i,j}^{n+1}({\left(T_g\right)}_{i,j}^{n+1}-{\left(T_s\right)}_{i,j}^{n+1}\text{})+{(}_{{\epsilon }_s}^{\stackrel{•}{{\displaystyle q}}}\end{array}$$E46

Gas Phase:

$$\begin{array}{l}{\rho }_g{C}_{p,g}\frac{{\left({\epsilon }_g{T}_g\right)}_{i,j}^{n+1}-{\left({\epsilon }_g{T}_g\right)}_{i,j}^n}{\Delta t}\\ +\frac{{\rho }_g{C}_{p,g}{\left(u_g\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left({\epsilon }_g{T}_g\right)}_{i,j}^n-{\left({\epsilon }_g{T}_g\right)}_{i-1,j}^n,if{\left(u_g\right)}_{i,j}^n\ge 0.0\\ {\left({\epsilon }_g{T}_g\right)}_{i+1,j}^n-{\left({\epsilon }_g{T}_g\right)}_{i,j}^n,if{\left(u_g\right)}_{i,j}^n<0.0\end{array}\\ \\ +\frac{{\rho }_s{C}_{p,g}{\left(v_g\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left({\epsilon }_g{T}_g\right)}_{i,j}^n-{\left({\epsilon }_g{T}_g\right)}_{i,j-1}^n,if{\left(v_g\right)}_{i,j}^n\ge 0.0\\ {\left({\epsilon }_g{T}_g\right)}_{i,j+1}^n-{\left({\epsilon }_g{T}_g\right)}_{i,j}^n,if{\left(v_g\right)}_{i,j}^n<0.0\end{array}\\ =\frac{{\left({\epsilon }_g{K}_g{T}_g\right)}_{i+1,j}^n-2{\left({\epsilon }_g{K}_g{T}_g\right)}_{i,j}^n+{\left({\epsilon }_g{K}_g{T}_g\right)}_{i-1,j}^n}{\Delta x^2}\\ +\frac{{\left({\epsilon }_g{K}_g{T}_g\right)}_{i,j+1}^n-2{\left({\epsilon }_g{K}_g{T}_g\right)}_{i,j}^n+{\left({\epsilon }_g{K}_g{T}_g\right)}_{i,j-1}^n}{\Delta y^2}\\ +{\left(h_v\right)}_{i,j}^{n+1}({\left(T_s\right)}_{i,j}^{n+1}-{\left(T_g\right)}_{i,j}^{n+1}\text{})\end{array}$$E47

6. Dimensionless numbers

In this section we define a group of four dimensionless numbers which mainly affect and control fluidization field.

6.1. Archimedes number" Ar"

Archimedes Number is used in characterization of the fluidized state and is defined as follow:

$$Ar=\frac{g{d}_p^3{\rho }_g\left({\rho }_s-{\rho }_g\right)}{{\mu }_g^2}$$E48

6.2. Density number “De”

Which define as the density ratio:

$$De=\frac{{\rho }_g}{{\rho }_s}$$E49

6.3. Flow number “fl”

Which is defined as:

$$fl=\frac{V_{g,in}}{u_t}$$E50

6.4. Dimensionless gas velocity "$${\Omega }^{1/3}$$"

Which is defined as:

$${\Omega }^{1/3}={\left[\frac{{\rho }_g^2}{{\mu }_{g}g\left({\rho }_s-{\rho }_g\right)}\right]}^{1/3}{V}_{g,in}$$E51

6.5. Dimensionless time “$$\tau$$”

Which is defined as:

$$\tau =\frac{t{u}_t}{d_p}$$E52

$$$$

7. Methodology of solution

The method of solution used in the present work is described in details in this section. The procedures of solution are as following:

1. Input the following parameters :

• Total time of calculation

• Total nodes number in x and y direction

• Total height of the bed and free board

• Width of fluidized bed

• Initial height of the bed

• Properties of gas phase (μ_g,ρ_g,k_g,C_p,g )

• Properties of solid particles phase (ρ_s,k_s,C_p,s )

• Entering gas velocity

• Acceleration of gravity

• Particle diameter

• Entering gas temperature

• Initial particles temperature

1. Calculate the minimum fluidized velocity from the following equation [6] :

$$V_{g,mf}=33.7\left[{\left(1+3.59X{10}^{-5}Ar\right)}^{0.5}-1\right]\frac{{\mu }_g}{\left(d_p{\rho }_g\right)}$$E53

1. Calculate the gas void fraction at minimum fluidized velocity from the following equation [7]:

$${\in }_{g,mf}=\frac{1}{2.1}[0.4+{(\frac{200{\mu }_g{V}_{g,mf}}{d_p^2({\rho }_s-{\rho }_g)g})}^{1/3}]$$E54

1. Calculate the gas void fraction at entering gas velocity which given by:

$${\in }_{g,in}=\frac{1}{2.1}[0.4+{(\frac{200{\mu }_g{V}_{g,in}}{d_p^2({\rho }_s-{\rho }_g)g})}^{1/3}]$$E55

1. Calculate the particle terminal velocity from the following equation [2]:

$$u_t={\left[-3.809+{\left({3.809}^2+1.832A{r}^{0.5}\right)}^{0.5}\right]}^2\frac{{\mu }_g}{{\rho }_g{d}_p}$$E56

1. Specify stability of fluidized bed.

2. $$$$Determine the Δt, Δx and Δy

3. Determine the fluidized bed height at the entering velocity using the equation [5]:

$$H1=\frac{\left(1-{\epsilon }_{g,mf}\right)}{\left(1-{\epsilon }_{g,in}\right)}{H}_{mf}$$E57

1. Specify Initial and boundary conditions.

2. Call subroutine “cont” to solve particles phase continuity equations and evaluate the new time step particles volume fractions ($$\text{}{\epsilon }_s^{n+1}$$) consequently evaluate ($${\epsilon }_g^{n+1}\text{}$$) from equation (41 ).

Note: we use the excess solid volume correction [8] in a special subroutine:

The correction works out as a posteriori redistribution of the particle phase volume fraction in excess in each cell where:

$${\epsilon }_s\le {\epsilon }_{s,\max }$$E58

$${\epsilon }_{s,\max }=1-{\epsilon }_{g,mf}$$E59

and if $${\epsilon }_s>{\epsilon }_{s,\max }$$then:

$${\epsilon }_s^{ex}={\epsilon }_s-{\epsilon }_{s,\max }$$E60

The balance may be expressed in terms of particle volume fraction:

$${\epsilon }_{s,i,j}^{new}={\epsilon }_{s,i,j}^{old}-{\epsilon }_{s,i,j}^{ex}+\frac{{\epsilon }_{s,i+1,j}^{ex}}{4}+\frac{{\epsilon }_{s,i-1,j}^{ex}}{4}+\frac{{\epsilon }_{s,i,j+1}^{ex}}{4}+\frac{{\epsilon }_{s,i,j-1}^{ex}}{4}$$E61

Figure (3) shows the correction mechanism:

Figure 3.

The excess solid volume correction

1. Call subroutine “dirx” to solve momentum equation in x-direction and evaluate at the new time step x-components velocities ($$u_g^{n+1}\text{}and\text{}{u}_s^{n+1}$$). In this subroutine, equations (14) and (43) are solved to get values of x-component new step velocities ($$u_g^{n+1}\text{}and\text{}{u}_s^{n+1}$$). Equation (43) is reduced to the following form :

$$A_{11}{(u_s)}_{i,j}^{n+1}+A_{12}{(u_g)}_{i,j}^{n+1}=A_{13}$$E62

where:

$$A_{11}=\frac{{\rho }_s}{\Delta t}$$E63

$$A_{12}=\frac{{\rho }_g}{\Delta t}$$E64

$$\begin{array}{l}{A}_{13}=\frac{{\rho }_s}{\Delta t}{\left(u_s\right)}_{i,j}^n+\frac{{\rho }_g}{\Delta t}{\left(u_g\right)}_{i,j}^n\\ -\frac{{\left(u_s\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left(u_s\right)}_{i,j}^n-{\left(u_s\right)}_{i-1,j}^n,if{\left(u_s\right)}_{i,j}^n\ge 0.0\\ {\left(u_s\right)}_{i+1,j}^n-{\left(u_s\right)}_{i,j}^n,if{\left(u_s\right)}_{i,j}^n<0.0\end{array}\\ -\frac{{\left(v_s\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left(u_s\right)}_{i,j}^n-{\left(u_s\right)}_{i,j-1}^n,if{\left(v_s\right)}_{i,j}^n\ge 0.0\\ {\left(u_s\right)}_{i,j+1}^n-{\left(u_s\right)}_{i,j}^n,if{\left(v_s\right)}_{i,j}^n<0.0\end{array}\\ -\frac{{\left(u_g\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left(u_g\right)}_{i,j}^n-{\left(u_g\right)}_{i-1,j}^n,if{\left(u_g\right)}_{i,j}^n\ge 0.0\\ {\left(u_g\right)}_{i+1,j}^n-{\left(u_g\right)}_{i,j}^n,if{\left(u_g\right)}_{i,j}^n<0.0\end{array}\\ -\frac{{\left(v_g\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left(u_g\right)}_{i,j}^n-{\left(u_g\right)}_{i,j-1}^n,if{\left(v_g\right)}_{i,j}^n\ge 0.0\\ {\left(u_g\right)}_{i,j+1}^n-{\left(u_g\right)}_{i,j}^n,if{\left(v_g\right)}_{i,j}^n<0.0\end{array}\\ \end{array}$$E65

so equation (62) and equation (14) result the following system of equations:

$$\left(\begin{array}{cc}\begin{array}{c}{A}_{11}\\ {\left({\epsilon }_s\right)}_{i,j}^{n+1}\end{array}& \begin{array}{c}{A}_{12}\\ {\left({\epsilon }_g\right)}_{i,j}^{n+1}\end{array}\end{array}\right)\left(\begin{array}{c}{\left(u_s\right)}_{i,j}^{n+1}\\ {\left(u_g\right)}_{i,j}^{n+1}\end{array}\right)=\left(\begin{array}{c}{A}_{13}\\ U_{gin}\end{array}\right)$$

Which are solved to get ($$u_g^{n+1}\text{}and\text{}{u}_s^{n+1}$$)

1. Call subroutine “diry” to solve momentum equation in y-direction and evaluate at the new time step y-components velocities ($$v_g^{n+1}\text{}and\text{}{v}_s^{n+1}$$). In this subroutine, equations (13) and (44) are used to get values of ($$v_g^{n+1}\text{}and\text{}{v}_s^{n+1}$$). Equation (44) is reduced to the following form :

$$B_{11}{(v_s)}_{i,j}^{n+1}+B_{12}{(v_g)}_{i,j}^{n+1}=B_{13}$$E66

where:

$$B_{11}=\frac{{\rho }_s}{\Delta t}+\beta \left(1+\frac{{\left({\epsilon }_s\right)}_{i,j}^n}{\left(1-{\left({\epsilon }_s\right)}_{i,j}^n\right)}\right)$$E67

where:

$$\beta =C_d\frac{3{\left({\epsilon }_s\right)}_{i,j}^n{\rho }_g\left|{\left(v_g\right)}_{i,j}^n-{\left(v_s\right)}_{i,j}^n\right|}{4d_p}{(1-{\left({\epsilon }_s\right)}_{i,j}^n)}^{-1.8}$$E68

$$B_{12}=-\frac{{\rho }_g}{\Delta t}-\beta \left(1+\frac{{\left({\epsilon }_s\right)}_{i,j}^n}{\left(1-{\left({\epsilon }_s\right)}_{i,j}^n\right)}\right)$$E69

$$\begin{array}{l}{B}_{13}=\frac{{\rho }_s}{\Delta t}{\left(v_s\right)}_{i,j}^n+\frac{{\rho }_g}{\Delta t}{\left(v_g\right)}_{i,j}^n\\ -\frac{{\left(u_s\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left(v_s\right)}_{i,j}^n-{\left(v_s\right)}_{i-1,j}^n,if{\left(u_s\right)}_{i,j}^n\ge 0.0\\ {\left(v_s\right)}_{i+1,j}^n-{\left(v_s\right)}_{i,j}^n,if{\left(u_s\right)}_{i,j}^n<0.0\end{array}\\ -\frac{{\left(v_s\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left(v_s\right)}_{i,j}^n-{\left(v_s\right)}_{i,j-1}^n,if{\left(v_s\right)}_{i,j}^n\ge 0.0\\ {\left(v_s\right)}_{i,j+1}^n-{\left(v_s\right)}_{i,j}^n,if{\left(v_s\right)}_{i,j}^n<0.0\end{array}\\ -\frac{{\left(u_g\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left(v_g\right)}_{i,j}^n-{\left(v_g\right)}_{i-1,j}^n,if{\left(u_g\right)}_{i,j}^n\ge 0.0\\ {\left(v_g\right)}_{i+1,j}^n-{\left(v_g\right)}_{i,j}^n,if{\left(u_g\right)}_{i,j}^n<0.0\end{array}\\ -\frac{{\left(v_g\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left(v_g\right)}_{i,j}^n-{\left(v_g\right)}_{i,j-1}^n,if{\left(v_g\right)}_{i,j}^n\ge 0.0\\ {\left(v_g\right)}_{i,j+1}^n-{\left(v_g\right)}_{i,j}^n,if{\left(v_g\right)}_{i,j}^n<0.0\end{array}\\ -\frac{{\left({\epsilon }_s\right)}_{i,j}^n-{\left({\epsilon }_s\right)}_{i,j-1}^n}{dy}\left(3.2g{d}_p\left({\rho }_s-{\rho }_g\right)\right)\end{array}$$E70

So equation (66) and equation (13) result the following system of equations:

$$\left(\begin{array}{cc}\begin{array}{c}{B}_{11}\\ {\left({\epsilon }_s\right)}_{i,j}^{n+1}\end{array}& \begin{array}{c}{B}_{12}\\ {\left({\epsilon }_g\right)}_{i,j}^{n+1}\end{array}\end{array}\right)\left(\begin{array}{c}{\left(v_s\right)}_{i,j}^{n+1}\\ {\left(v_g\right)}_{i,j}^{n+1}\end{array}\right)=\left(\begin{array}{c}{B}_{13}\\ V_{gin}\end{array}\right)$$

Which are solved to get ($$v_g^{n+1}\text{}and\text{}{v}_s^{n+1}$$)

1. Call subroutine “temp” to solve energy equation and evaluate the new time step gas and solid particles temperatures ($$T_g^{n+1}\text{}and\text{}{T}_s^{n+1}$$).In this subroutine, equations (46) and (47) are used to get values of new time step temperatures for both phases ($$T_g^{n+1}\text{}and\text{}{T}_s^{n+1}$$). Equation (46) is reduced to the following form :

$$C_{11}{(T_s)}_{i,j}^{n+1}+C_{12}{(T_g)}_{i,j}^{n+1}=C_{13}$$E71

where:

$$C_{11}=\frac{{\rho }_s{C}_{p,s}}{\Delta t}\left({\epsilon }_s^{n+1}\right)+{\left(h_v\right)}_{i,j}^{n+1}$$E72

$$C_{12}={\left(h_v\right)}_{i,j}^{n+1}$$E73

$$\begin{array}{l}{C}_{13}=\frac{{\rho }_s{C}_{p,s}}{\Delta t}{\left({\epsilon }_s{T}_s\right)}_{i,j}^n\\ -\frac{{\rho }_s{C}_{p,s}{\left(u_s\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left({\epsilon }_s{T}_s\right)}_{i,j}^n-{\left({\epsilon }_s{T}_s\right)}_{i-1,j}^n,if{\left(u_s\right)}_{i,j}^n\ge 0.0\\ {\left({\epsilon }_s{T}_s\right)}_{i+1,j}^n-{\left({\epsilon }_s{T}_s\right)}_{i,j}^n,if{\left(u_s\right)}_{i,j}^n<0.0\end{array}\\ -\frac{{\rho }_s{C}_{p,s}{\left(v_s\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left({\epsilon }_s{T}_s\right)}_{i,j}^n-{\left({\epsilon }_s{T}_s\right)}_{i,j-1}^n,if{\left(v_s\right)}_{i,j}^n\ge 0.0\\ {\left({\epsilon }_s{T}_s\right)}_{i,j+1}^n-{\left({\epsilon }_s{T}_s\right)}_{i,j}^n,if{\left(v_s\right)}_{i,j}^n<0.0\end{array}\\ +\frac{{\left({\epsilon }_s{K}_s{T}_s\right)}_{i+1,j}^n-2{\left({\epsilon }_s{K}_s{T}_s\right)}_{i,j}^n+{\left({\epsilon }_s{K}_s{T}_s\right)}_{i-1,j}^n}{\Delta x^2}\\ +\frac{{\left({\epsilon }_s{K}_s{T}_s\right)}_{i,j+1}^n-2{\left({\epsilon }_s{K}_s{T}_s\right)}_{i,j}^n+{\left({\epsilon }_s{K}_s{T}_s\right)}_{i,j-1}^n}{\Delta y^2}+{(}_{{\epsilon }_s}^{\stackrel{•}{{\displaystyle q}}}\end{array}$$E74

Also equation (47) is reduced to the following form:

$$D_{11}{(T_s)}_{i,j}^{n+1}+D_{12}{(T_g)}_{i,j}^{n+1}=D_{13}$$E75

where:

$$D_{11}={\left(h_v\right)}_{i,j}^{n+1}$$E76

$$D_{12}=\frac{{\rho }_g{C}_{p,g}}{\Delta t}\left({\epsilon }_g^{n+1}\right)+{\left(h_v\right)}_{i,j}^{n+1}$$E77

$$\begin{array}{l}{D}_{13}={\rho }_g{C}_{p,g}\frac{{\left({\epsilon }_g{T}_g\right)}_{i,j}^n}{\Delta t}\\ -\frac{{\rho }_g{C}_{p,g}{\left(u_g\right)}_{i,j}^n}{\Delta x}\{\begin{array}{c}{\left({\epsilon }_g{T}_g\right)}_{i,j}^n-{\left({\epsilon }_g{T}_g\right)}_{i-1,j}^n,if{\left(u_g\right)}_{i,j}^n\ge 0.\\ {\left({\epsilon }_g{T}_g\right)}_{i+1,j}^n-{\left({\epsilon }_g{T}_g\right)}_{i,j}^n,if{\left(u_g\right)}_{i,j}^n<0.0\end{array}\\ -\frac{{\rho }_s{C}_{p,g}{\left(v_g\right)}_{i,j}^n}{\Delta y}\{\begin{array}{c}{\left({\epsilon }_g{T}_g\right)}_{i,j}^n-{\left({\epsilon }_g{T}_g\right)}_{i,j-1}^n,if{\left(v_g\right)}_{i,j}^n\ge 0.0\\ {\left({\epsilon }_g{T}_g\right)}_{i,j+1}^n-{\left({\epsilon }_g{T}_g\right)}_{i,j}^n,if{\left(v_g\right)}_{i,j}^n<0.0\end{array}\\ +\frac{{\left({\epsilon }_g{K}_g{T}_g\right)}_{i+1,j}^n-2{\left({\epsilon }_g{K}_g{T}_g\right)}_{i,j}^n+{\left({\epsilon }_g{K}_g{T}_g\right)}_{i-1,j}^n}{\Delta x^2}\\ +\frac{{\left({\epsilon }_g{K}_g{T}_g\right)}_{i,j+1}^n-2{\left({\epsilon }_g{K}_g{T}_g\right)}_{i,j}^n+{\left({\epsilon }_g{K}_g{T}_g\right)}_{i,j-1}^n}{\Delta y^2}\end{array}$$E78

so equation (71) and equation (75) result the following system of equations:

$$\left(\begin{array}{cc}\begin{array}{c}{C}_{11}\\ D_{11}\end{array}& \begin{array}{c}{C}_{12}\\ D_{12}\end{array}\end{array}\right)\left(\begin{array}{c}{\left(T_s\right)}_{i,j}^{n+1}\\ {\left(T_g\right)}_{i,j}^{n+1}\end{array}\right)=\left(\begin{array}{c}{C}_{13}\\ D_{13}\end{array}\right)$$

Which are solved to get ($$T_g^{n+1}\text{}and\text{}{T}_s^{n+1}$$)

1. Make gas residual check,which given from the equation (42):

• If $$\left|d_g\left(i,j\right)\right|\le \delta$$ go to step 15, where $$\delta$$ is a small positive value $$\delta$$=5X10^-3.

• Else adjust $${\epsilon }_g^{n+1}$$by :

$$\begin{array}{l}{\left({\epsilon }_g\right)}_{i,j}^{n+1}={\left({\epsilon }_g\right)}_{i,j}^n-\frac{dt}{\Delta x}{\left(u_g\right)}_{i,j}^{n+1}\{\begin{array}{c}{\left({\epsilon }_g\right)}_{i,j}^{n+1}-{\left({\epsilon }_g\right)}_{i-1,j}^{n+1},if{\left(u_g\right)}_{i,j}^{n+1}\ge 0.0\\ {\left({\epsilon }_g\right)}_{i+1,j}^{n+1}-{\left({\epsilon }_g\right)}_{i,j}^{n+1},if{\left(u_g\right)}_{i,j}^{n+1}<0.0\end{array}\\ -\frac{\Delta t}{\Delta y}{\left(v_g\right)}_{i,j}^{n+1}\{\begin{array}{c}{\left({\epsilon }_g\right)}_{i,j}^{n+1}-{\left({\epsilon }_g\right)}_{i,j-1}^{n+1},if{\left(v_g\right)}_{i,j}^{n+1}\ge 0.0\\ {\left({\epsilon }_g\right)}_{i,j+1}^{n+1}-{\left({\epsilon }_g\right)}_{i,j}^{n+1},if{\left(v_g\right)}_{i,j}^{n+1}<0.0\end{array}\\ \end{array}$$E79

and calculate $${\epsilon }_s^{n+1}$$ from equation (41).

• Go to step 11 and calculate new time step velocities

1. Calculate the gas pressure values.

2. Calculate dimensionless numbers.

3. End program.

8. Parametric study of the hydrodynamic and thermal results

The present work results show the effect of variation of several bed parameters such as particle diameter, terminal velocity of the particle, minimum fluidized velocity of the particle input gas velocity, fluidized material type and heat generation by particles on the hydrodynamic and thermal behavior of fluidized bed. In this section the effect of different parameters in hydrodynamic and thermal performance of fluidized bed is analyzed in detail.

8.1. Effect of particle diameter

Particle diameter is the most influential parameter in the overall fluidized bed performance. In view of that fact, the bed material in a fluidized bed is characterized by a wide range of particle diameter, so that the effect of particle diameter is analyzed in details. In this section particle diameter is changed from 100 μm to 1000 μm for sand as fluidized material to study the effect of particle diameter on the fluidization performance.

One of the important parameter of the fluidized bed study is the total pressure drop across the bed. Although it is constant after beginning of fluidization and equal to the weight of the bed approximately. But its value changes with change of particle diameter.The effect of change of particle diameter on total pressure drop is very important in design and cost of fluidized bed. Figure (4) shows that effect for sand particles of different diameters (100, 200, 300, 400, 500, 600, 700, 800, 900 and 1000 μm). It is clear that the pressure drop increase with increase of particle diameter, which means that small particle is betters in design of fluidized bed cost.

A key parameter in the Thermal analysis of fluidized bed is the average heat transfer coefficient. Figure (5) shows the variation of average heat transfer coefficient with the time for different particles diameter from 100 to 1000 μm sand particles.It is observed from the figure that particles with small size show higher values of average heat transfer coefficient. The particles with diameter 100 μm show higher for average heat transfer coefficient reaches to about 3 times of particles with diameter 1000 μm.

Figure 4.

Effect of particle diameter on total pressure drop across fluidized bed

Figure 5.

Effect of particle diameter on average heat transfer coefficient in fluidized bed

Particles diameter determines the type of particles on Geldart diagram, consequently the behavior of the fluidized bed.

Figures from (6) to (13) illustrate the effect of particles diameter on hydrodynamic and thermal behavior for particle type-B. These figures explain the high disturbance in different bed parameters (ε_g, ρ_sus,v_g,v_s,u_g,u_s,T_g,T_s) due to the bubbles formation. The disturbance decreases with increase of particles diameter. This may be due to come near D-type under uniform fluidization.

Fluidized bed is used for wide range of particle diameter. It is better to fluidize particle Dtype in spouted bed to decrease the pressure drop. However, using uniform fluidization for D type give a good and stable thermal behavior of fluidized bed. So in some applications like nuclear reactors the stability and safety is important than cost of pumping power.

Figures from (14) to (21) show the effect of change particle diameter on hydrodynamic and thermal behavior for particle D-type. For this particles type, the behavior of fluidized bed is more uniform in performance than B-type [be consistent with usage and define what the different types are]. Although D-type gives butter fluidization in spouted bed but it gives good performance under uniform fluidization with high pressure drop as shown in figure (4). This means an increase in pumping power and costs to achieve uniform fluidization.

Figure 6.

Effect of particle diameter on gas volume fraction (B-type)

Figure 7.

Effect of particle diameter on suspension density (B-type)

Figure 8.

Effect of particle diameter on vertical gas velocity (B-type)

Figure 9.

Effect of change particle diameter on vertical particle velocity (B-type)

Figure 10.

Effect of particle diameter on horizontal particle velocity (B-type)

Figure 11.

Effect of particle diameter on horizontal gas velocity (B-Type)

Figure 12.

Effect of particle diameter on gas temperature (B-type)

Figure 13.

Effect of particle diameter on particle temperature (B-type)

Figure 14.

Effect of particle diameter on gas volume fraction (D-type)

Figure 15.

Effect of particle diameter on suspension density (D-type)

Figure 16.

Effect of particle diameter on vertical gas velocity (D-type)

Figure 17.

Effect of change particle diameter on vertical particle velocity (D-Type)

Figure 18.

Effect of particle diameter on horizontal particle velocity (D-type)

Figure 19.

Effect of particle diameter on horizontal gas particle velocity (D-type)

Figure 20.

Effect of particle diameter on gas temperature (D-type)

Figure 21.

Effect of particle diameter on particle temperature (D-type)

8.2. Effect of input gas velocity

In this section the input gas velocity is changed from one to nine times minimum fluidized velocity for 500 μm sand particles to study the the effect of this parameter on fluidization behavior.

Input gas velocity has an effective role in thermal performance of the fluidized bed. In order to illustrate its effect, the relation between Nusselt number and flow number is described in figure (22). It is clear from this figure that with increase in flow numbers, the Nusselt number increases until reached an optimum flow number where the Nusselt number reaches its maximum value. After this optimum value the increase in flow number is associated with a decrease in Nusselt number. This decrease in Nusselt number may be due to the increase of input gas velocity toward the terminal velocity, consequently the bed goes to be empty bed.

Figure (23) shows the variation of Nusselt number with velocity number. Also the relation between Nusselt number and velocity number has the same trend as the Nusselt number with flow number. This confirms the result from figure (22). This means that there is an optimum input gas velocity to yield the best heat transfer characteristics. This velocity is the target of the fluidized bed designer. The value of this velocity depends on the fluidized gas properties, the fluidized material, particle diameter, and bed geometry.

Figure 22.

Variation of Nusselt number with flow number

Figure 23.

Variation of Nusselt number with velocity number

8.3. Effect of fluidized material type

Type of fluidized material controls hydrodynamic and thermal performance of fluidized bed.It affects on the different parameters of fluidization such as gas volume fraction, suspension density, gas velocity distribution and particle velocity distribution, gas phase temperature and particle phase temperature. In this section different types of materials such as sand, marble, lead, copper, aluminum and steel of particle diameters 1mm are used to study the effect of fluidized materials on the bed performance.

Figure (24) shows the change of gas volume fraction with time for different types of fluidized materials. The figure shows that the gas volume fraction of copper is the highest value followed by lead and steel. Gas volume fraction of sand, marble and aluminum are at the same level.

The change of suspension density with time for several types of fluidized materials is shown in Figure (25). It is clear that the highest suspension density lies with the material of the highest density.

Figures (26) to (29) show the effect of change of fluidized material type on horizontal and vertical velocities of gas and particle.

Figures (30) and (40) illustrate the effect of change of fluidized material type on particle and gas temperatures.

Figure 24.

Effect of fluidized material type on gas volume fraction

Figure 25.

Effect of fluidized material type on suspension density

Figure 26.

Effect of fluidized material type on vertical particle velocity

Figure 27.

Effect of fluidized material type on horizontal particle velocity

Figure 28.

Effect of fluidized material type on vertical gas velocity

Figure 29.

Effect of fluidized material type on horizontal gas velocity

Figure 30.

Effect of fluidized material type on particle temperature

Figure 31.

Effect of fluidized material type on gas temperature

8.4. Effect of heat generation by particles

The aluminum particles of 1 mm diameter are fluidized with different value of heat generated in particles (0, 500,1000,1500,2000 and 3000 Watt), the effect of change of heat generated in particles is studied. With the increase of heat generated by particles the gas phase temperature increases as shown in Figure (32). The value of the increase in gas temperature is approximately in range of 3 ^°C. Figure (33) shows that the particle phase temperature increases with the increase of heat generated by particles. The range of increase is about 45^°C. It is clear that the rate of increase in particle temperature is more than the rate of increase in gas temperature, consequently the temperature difference between the two phases increase. Figure (34) illustrates the relation between average heat transfer coefficient and heat generated by particles. The results of the present work shows that the average heat transfer coefficient dos not depend on heat generated by particles and all heat generated by particles converts to temperature difference between the two phases. This result agrees with that of reference [9].

8.5. Terminal velocity effect

Figure (35) shows the relation between terminal velocity and average heat transfer coefficient. It is clear from the figure that with the increase in terminal velocity the average heat transfer coefficient decreases.

8.6. Minimum fluidized velocity effect

Minimum fluidized velocity is the most important parameter in study of fluidization. This velocity distinguishes the fluidized bed from a packed bed and is an indicator that fluidization is occurred. Figure (36) shows the variation of average heat transfer coefficient with minimum fluidized velocity. The average heat transfer coefficient decreases with the increase of minimum fluidized velocity.

Figure 32.

Effect of particle heat generation on gas temperature

Figure 33.

Effect of particle heat generation in particle temperature

Figure 34.

Effect of particle heat generation on average heat transfer coefficient

Figure 35.

Effect of terminal velocity on average heat transfer coefficient

Figure 36.

Effect of Minimum Velocity on Average Heat Transfer Coefficient

Nomenclature

Ar	Archimedes number, $$Ar=\frac{g{d}_p^3{\rho }_g\left({\rho }_s-{\rho }_g\right)}{{\mu }_g^2}$$
Cd	Drag coefficient, $$C_d={\left(0.63+\frac{4.8}{{\mathrm{Re}}_{}^{0.5}}\right)}^2$$
Cp,g	Specific heat of fluidizing gas at constant pressure	J/kg.K
Cp,s	Specific heat of solid particles	J/kg.K
dg	Residue of the gas continuity equation	Kg/m^3 s
dp	Mean particle diameter	M
G	Acceleration due to gravity	m/s^2
Fgx	Total gas phase force in x direction per unit volume	N /m^3
Fgy	Total gas phase force in y direction per unit volume	N /m^3
Fsx	Total particle phase force in x direction per unit volume	N /m^3
Fsy	Total particle phase force in y direction per unit volume	N /m^3
fl	Flow number , $$fl=\frac{V_{g,in}}{u_t}$$
hgp	Heat transfer coefficient between gas phase and particle phase	W/m^2.K
hv	Volumetric heat transfer coefficient , $$h_v=\frac{6\left(1-{\epsilon }_g\right)h_{gp}}{d_p}$$	W/m^3.K
H	Total height of the bed and freeboard	M
Hmf	Minimum fluidized head of the bed	M
H1	Expansion head of bed at the input velocity	M
kg	Thermal conductivity of gas phase	W/m.K
ks	Thermal conductivity of particle phase	W/m.K
L	Width of the bed	M
Nu	Nusselt number based on particle diameter, (Nu= hgpdp/kg)
Pg	Gas pressure	Pa
Pr	Prandtl number , Pr= μgCpg/kg
$$\stackrel{•}{q}$$	Rate of heat generated within particle phase	W/m^3
ug	Gas phase velocity in x direction	m/s
Ugin	Input gas velocity to the bed in x direction	m/s
us	Particle phase velocity in x direction	m/s
ut	Particle terminal velocity	m/s
Ur	Relative velocity between two phases	m/s
Re	Reynolds number, $$\mathrm{Re}=\frac{{\epsilon }_g{\rho }_g{d}_p\left|\overrightarrow{U_r}\right|}{{\mu }_g}$$
S	Stability function
T	Time	s
Tg	Gas phase temperature	C^◦
Ts	Particle phase temperature	C^◦
TFM	Two fluid model
vg	gas phase velocity in y direction	m/s
Vg,mf	Gas minimum fluidized velocity	m/s
Vg,in	Input gas velocity to the bed in y direction	m/s
vs	Particle phase velocity in y direction	m/s

Greek Letters

ΔPB	Total Pressure drop across the bed	Pa
Δt	Time step	s
Δx	Length of cell in the computational grid	m
Δy	height of cell in the computational grid	m
∈g	Gas phase volume fraction
∈g,mf	Gas phase volume fraction at minimum fluidization
∈in	Gas volume fraction at the input velocity
∈s	Particle phase volume fraction
ρg	Density of gas phase	Kg/m^3
ρs	Density of particle phase	Kg/m^3
ρsus	Suspension density	Kg/m^3
μg	Viscosity of gas	Pa.s
δ	Small positive value = 5X10^-3