Application of the Two-Fluid Model with Kinetic Theory of Granular Flow in Liquid–Solid Fluidized Beds

Abstract

Numerical simulation of voidage distributions and bed expansions is carried out in a liquid–solid fluidized bed in the present work. Effects of drag force models as well as virtual mass force and lift force are studied in the prediction of particle flow characteristics; simulated results indicated that both virtual mass force and lift force could not be neglected in liquid–solid fluidized bed. Different superficial velocities of liquid phase are also studied to investigate the effects of operating conditions on the distribution of particle concentration and velocities. The coefficient of restitution varied from 0.6 to 0.99, and the effects of radial distribution function models on granular pressure and granular temperature are also studied. Different drag models exhibit various particle velocity distributions, while the Gibilaro drag model failed in predicting the liquid–solid drag to some extent in this study. A comprehensive simulation model was proposed for predicting the two-phase flow characteristics in the liquid–solid fluidized bed. Predicted axial void fraction agrees qualitatively and quantitatively well with the experimental results in the literature.

1. Introduction

The liquid–solid fluidized bed reactors are widely used in the pharmaceutical, chemical, food, petroleum, and many other industries; as a result, they have been the focus of much research. Liquid–solid fluidized beds exhibit a great variety of complex inhomogeneous flow structures. The origin and hierarchy of these of structures in liquid–solid fluidized beds have presented a challenge for both experimental and numerical research. A general understanding of their hydrodynamic behaviour is still under pursuit. One of the issues that make the accurate prediction of particulate flows to be difficult in liquid–solid flow system is the lack of accurate comprehensive simulation models and corresponded parameters.

Computational fluid dynamics (CFD) method can provide a series of fluid hydrodynamic information, which can hardly be obtained by modern measuring instruments. In the CFD modelling, the two-fluid model (TFM) assumes the liquid phase and solid phase as both continua and fully interpenetrating within each other. Among various attempts to formulate particulate flow stresses, the kinetic theory of granular flow (KTGF) is usually employed [1], which is an extension of the classical kinetic theory of gases to dense particulate flows. In this theory, the fluctuation energy of particles was described by introducing the concept of granular temperature. Thus, the particle flow behaviour can be predicted by TFM-KTGF model. A number of studies had shown the capability of the KTGF approach for modelling fluidized beds [2, 3, 4, 5, 6, 7].

In the TFM, the momentum transfer between fluid and particle phases is of the great significance for the momentum of both phases. An accurate closure law for fluid-particle interactions is highly required. Generally, the interaction terms in liquid–solid flow system include the drag force, the virtual mass force and the history force except that the pressure gradient and the gravity force. The momentum exchange is mainly represented by the drag force [8]. Hence, the drag force models are important in simulating the interphase momentum transfer between the liquid and solid phases. Traditionally, the drag force models are average-based in the literature [9, 10]. With the development of computational methods and instruments, numerous CFD models were applied to the simulation of dynamic processes in the liquid–solid circulating fluidized beds. Roy and Dudukovic [11] used TFM model combining with the KTGF to simulate the flow behaviours in liquid–solid circulating fluidized beds. Razzak [12] employed the KTGF based on TFM and simulated the particle viscosity and particle pressure, and a drag model proposed by Wen and Yu was adopted to calculate the interphase momentum exchange. Cheng and Zhu [13, 14] made a comprehensive study on the modelling and simulation of hydrodynamics in liquid–solid circulating fluidized beds using both similitude method and CFD technique.

Some other studies took the local inhomogeneity of the liquid–solid flow in a circulating fluidized bed into account to calculate the interphase momentum exchange. Liu et al. [15, 16, 17] proposed a multi-scale drag coefficient model that can show better capability of predicting distribution of particle concentration. Xie et al. [18] proposed a series of correlations of KTGF model for liquid–solid flow by calculating solid pressure and viscosity and found that the particle-particle interactions can affect suspension characteristics for large particle size and high solid loading systems. Rahaman et al. [19] tested and validated three established empirical drag law correlations used to explain momentum exchange between solid and liquid phases. It was found that Wen and Yu [9] and Gidaspow [1] drag law models showed greater predictive power in terms of pressure drop and voidage in the fluidized beds of multi-particle systems. Ozel et al. [20] compared the direct numerical simulation of a liquid–solid fluidized bed with experimental data and found that fluid velocity fluctuations were mainly driven by fluid–particle wake interactions (pseudo-turbulence) whereas the particle velocity fluctuations derive essentially from the large-scale flow motion (recirculation). However, there is still an absence of comprehensive evaluation of such models for a better numerical simulation.

In present study, a comprehensive investigation of the models, kinetic theory constitutive parameters as well as models are performed; their validity in predicting the liquid–solid fluidization is compared and evaluated.

2. Liquid–solid two-fluid model

In the present work, an Eulerian multi-fluid model, which considers the conservation of mass and momentum for the solid and liquid phases, has been adopted. The kinetic theory of granular flow, which considers the conservation of solid fluctuation energy, has been used for closure. Conservation equations of mass and momentum of both phases result from the statistical average of instantaneous local transport equations. The governing equations are given below.

2.1. Governing equations

Both phases are continuous assuming a single liquid phase and single solid phase. The continuity and the momentum balance for both the phases are given. Interphase momentum transfer term includes the drag force, virtual mass force and lift force.

For simplification, we assume that (1) both liquid and solid phases are assumed to be isothermal; it is also assumed that there is no interphase mass transfer and (2) the solid phase is characterised by a spherical configuration with mean particle diameter and density. The continuity for phase i (i = l for liquid phase; i = s for solid phase):

$$\frac{\partial }{\partial t}\left({\epsilon }_i{\rho }_i\right)+\frac{\partial }{\partial x}\left({\epsilon }_i{\rho }_i{\mathit{u}}_i\right)=0$$E1

The momentum balance for the liquid phase is expressed by the Navier–Stokes equation, such equation is modified to include the drag force, virtual mass force and lift force to consider the momentum transfer between phases.

$$\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\epsilon }_l{\rho }_l{\mathit{u}}_l\right)+\nabla \cdot \left({\epsilon }_l{\rho }_l{\mathit{u}}_l{\mathit{u}}_l\right)={\epsilon }_l\nabla \cdot {\overline{\mathit{\tau }}}_l+{\epsilon }_l{\rho }_l\mathit{g}-{\epsilon }_l\nabla p-\beta \left({\mathit{u}}_l-{\mathit{u}}_s\right)-{\mathit{F}}_{\mathit{vm}}-{\mathit{F}}_{\mathit{lf}}$$E2

The stress tensor of liquid phase can be represented as:

$${\overline{\mathit{\tau }}}_l={\mu }_f\left[\nabla {\mathit{u}}_l+{\left(\nabla {\mathit{u}}_l\right)}^T\right]-\frac{2}{3}{\mu }_f\left(\nabla \cdot {\mathit{u}}_l\right)\overline{\mathit{I}}$$E3

where ${\mu }_f$ is combined with the laminar part and turbulent part and could be expressed as ${\mu }_f={\mu }_l+{\mu }_t$. The turbulent viscosity for the liquid phase is calculated as ${\mu }_t=c_{\mu }{\rho }_l{k}^2/\epsilon$. The liquid phase is described by a standard $k-\epsilon$ turbulence model.

The solid phase momentum balance is given as follows:

$$\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\epsilon }_s{\rho }_s{\mathit{u}}_s\right)+\nabla \cdot \left({\epsilon }_s{\rho }_s{\mathit{u}}_s{\mathit{u}}_s\right)=-{\epsilon }_s\nabla p-\nabla p_s+\nabla \cdot {\overline{\mathit{\tau }}}_s+{\epsilon }_s{\rho }_s\mathbf{g}+\beta \left({\mathit{u}}_l-{\mathit{u}}_s\right)+F_{\mathit{vm}}+F_{\mathit{lf}}$$E4

The solid phase stress tensor can be expressed in terms of the bulk solid viscosity ξ_s, and shear solid viscosity μ_s

$${\overline{\mathit{\tau }}}_s={\mu }_s\left\{\left[\nabla {\mathit{u}}_s+{\left(\nabla {\mathit{u}}_s\right)}^T\right]-\frac{2}{3}\left(\nabla \cdot {\mathit{u}}_s\right)\overline{\mathit{I}}\right\}+{\xi }_s\nabla \cdot {\mathit{u}}_s\overline{\mathit{I}}$$E5

2.2. Constitutive correlations

Analogical to the thermodynamic temperature for gases, the granular temperature θ = C²/3 was introduced as a measure for the energy of the fluctuating velocity of the particles [21]. The equation of conservation of solids fluctuating energy can be expressed as:

$$\frac{3}{2}\left[\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\epsilon }_s{\rho }_s\theta \right)+\nabla \cdot \left({\epsilon }_s{\rho }_s\theta \right){\mathit{u}}_s\right]=\left(-\nabla p_s\overline{\mathit{I}}+{\overline{\mathit{\tau }}}_s\right):\nabla {\mathit{u}}_s+\nabla \cdot \left(k_s\nabla \theta \right)-{\gamma }_s-3\mathit{\beta \theta }+D_{\mathit{ls}}$$E6

Its description is based on the kinetic theory of granular flow, where both the kinetic and the collisional influence are taken into account. The particle pressure can be calculated as follows:

$$p_s={\epsilon }_s{\rho }_s\theta \left[\frac{1}{1+\lambda /L}+2{\epsilon }_s{g}_0\left(1+e\right)\right]$$E7

g₀ is the radial distribution function at contact; it is employed to describe how density varies as a function of distance from a reference particle, which is a correlation factor that modifies the probability of collisions between grains when the solid granular phase becomes dense and can be regarded as a measure for the probability of inter-particle contact. Lun et al. [21] used the following equation for calculating radial distribution function.

$$g_0={\left[1-{\left({\epsilon }_s/{\epsilon }_{s,\text{max}}\right)}^{1/3}\right]}^{-1}$$E8

In this work, the equation proposed by Bagnold [22] is used, where ${\epsilon }_{s,\text{max}}$ is the maximum particle concentration at random packing.

Arastoopour et al. [23] also derived the similar forms of equation for calculating granular pressure.

$$g_0=\frac{1}{\left(1-\frac{{\epsilon }_s}{{\epsilon }_{s,\text{max}}}\right)}+\frac{3}{2}{d}_l\sum _{k=1}^N\frac{{\epsilon }_k}{d_k}$$E9

Ma and Ahmadi [24, 25] derived another form of radial distribution function; such form of the function is similar to that of the Arastoopour.

$$g_0=\frac{1+2.5{\epsilon }_s+4.59{\epsilon }_s^2+4.52{\epsilon }_s^3}{{\left(1-{\left(\frac{{\epsilon }_s}{{\epsilon }_{s,\text{max}}}\right)}^3\right)}^{0.687}}+\frac{1}{2}{d}_l\sum _{k=1}^N\frac{{\epsilon }_k}{d_k}$$E10

Lebowitz [26] derived the radial distribution function at contact for a mixture of hard spheres:

$$g_0=\frac{1}{\left(1-{\epsilon }_s\right)}+\frac{3\sum _{k=1}^N\frac{{\epsilon }_k}{d_k}}{{\left(1-{\epsilon }_s\right)}^2\left(d_l+d_k\right)}{d}_k{d}_l$$E11

The shear viscosity ${\mu }_s$ accounts for the tangential forces. It is capable of combining different inter-particle forces and using a momentum balance similar to that of a true continuous fluid. It is composed of three parts: ${\mu }_{s,\mathit{col}}$ represents the viscosity induced from the particle collisions, ${\mu }_{s,\mathit{kin}}$ represents the viscosity induced from particle fluctuations, and ${\mu }_{s,\mathit{fr}}$ from particle frictions; thus, solid viscosity can be expressed as:

$${\mu }_s={\mu }_{s,\mathit{col}}+{\mu }_{s,\mathit{kin}}+{\mu }_{s,\mathit{fr}}$$E12

$${\mu }_{s,\mathit{col}}=\frac{5{\rho }_s{d}_s}{96{\epsilon }_s}\sqrt{\frac{\theta }{\pi }}\left[\left(\frac{8{\epsilon }_s}{5\left(2-\eta \right)}\right)\left(1+\frac{8}{5}\eta \left(3\eta -2\right){\epsilon }_s{g}_0\right)+\frac{768}{25\pi }{\mathit{\eta \epsilon }}_s^2{g}_0\right]$$E13

$${\mu }_{s,\mathit{kin}}=\frac{5{\rho }_s{d}_s}{96{\epsilon }_s\eta \left(2-\eta \right)g_0}\sqrt{\frac{\theta }{\pi }}\left(\frac{1}{1+\lambda /L}\right)\left[1+\frac{8}{5}\eta \left(3\eta -2\right){\epsilon }_s{g}_0\right]$$E14

$${\mu }_{s,\mathit{fr}}=\frac{p_s\text{sin}\varphi }{2\sqrt{{\mathit{I}}_{2D}}}$$E15

where,

$$\eta =\left(1+e\right)/2$$E16

The bulk viscosity ${\xi }_s$ formulates the resistance of solid particles to compression and expansion. The following equation given by Ding and Gidaspow [2] is used in this work:

$${\xi }_s=\frac{4}{3}{\epsilon }_s^2{\rho }_s{d}_s{g}_0\left(1+e\right)\sqrt{\theta /\pi }$$E17

For the conductivity of granular energy k_s, following correlation is used:

$$k_s=\frac{25{\mathit{\pi \rho }}_s{d}_s\sqrt{\theta }}{128}\left[\left(\frac{1}{1+\lambda /L}\frac{8}{\eta g_0}+\frac{96{\epsilon }_s}{5}\right)\left(\frac{1+\frac{12}{5}{\eta }^2\left(4\eta -3\right){\epsilon }_s{g}_0}{41-33\eta }\right)+\frac{512{\mathit{\eta \epsilon }}_s^2{g}_0}{25\pi }\right]$$E18

The rate of dissipation of fluctuation kinetic energy due to particle collisions is expressed as,

$${\gamma }_s=3\left(1-e^2\right){\epsilon }_s^2{\rho }_s{g}_o\theta \left(\frac{4}{d_s}\sqrt{\frac{\theta }{\pi }}-\nabla \cdot {\mathit{u}}_s\right)$$E19

The last term D_ls in Eq. (6) is the rate of energy dissipation per unit volume caused by the transfer of liquid phase fluctuations to the particle phase fluctuations. In this study, the value of D_ls is calculated with Koch’s expression [27] as follows:

$$D_{\mathit{ls}}=\frac{d_s{\rho }_s}{4\sqrt{\mathit{\pi \theta }}{g}_o}{\left(\frac{18{\mu }_l}{d^2{\rho }_s}\right)}^2{\left|{\mathit{u}}_l-{\mathit{u}}_s\right|}^2$$E20

2.3. Drag model

The inter-phase drag term in the liquid and solid phase momentum equations is expressed as $\beta \left({\mathit{u}}_l-{\mathit{u}}_s\right)$, the product of the interphase momentum exchange coefficient $\beta$ and the slip velocity (relative velocity of the solid phase to the liquid phase). The inter-phase drag coefficient $\beta$ can be expressed by the correlations given by Gidaspow [1]. This correlation is a combination of the works of Ergun [28] and Wen and Yu [9]; the formulation presented by Ergun [28] is used at the porosity less than and equal to 0.8 where the suspension is dense, whereas the formulation provided by Wen and Yu [9] is used for the porosity greater than 0.8, where the suspension is regarded as dilute.

$${\beta }_{\mathit{\text{Ergun}}}=150\frac{{\left(1-{\epsilon }_l\right)}^2{\mu }_l}{{\left({\epsilon }_l{d}_s\right)}^2}+1.75\frac{{\rho }_l\left(1-{\epsilon }_l\right)\left|{\mathit{u}}_l-{\mathit{u}}_s\right|}{{\epsilon }_l{d}_s}{\epsilon }_l\le 0.8$$E21

$${\beta }_{\mathit{Wen}\&\mathit{Yu}}=\frac{3}{4}{C}_d\frac{{\rho }_l\left(1-{\epsilon }_l\right)\left|{\mathit{u}}_l-{\mathit{u}}_s\right|}{d_s}{\epsilon }_l^{-2.65}{\epsilon }_l>0.8$$E22

$$C_d=\left\{\begin{array}{l}\frac{24}{\mathit{Re}}\left(1+0.15{\mathit{Re}}^{0.687}\right)\text{Re}\le 1000\\ 0.44\text{Re}>1000\end{array}\right.$$E23

However, the transition proposed by Gidaspow [1] results in the drag law discontinuous in solid concentration even if it is continuous in Reynolds number. As a matter of fact, the drag force is a continuous function for both Reynolds number and solid concentration, the abrupt change in drag at bed voidage equals to 0.8 and can also cause numerical instabilities [29], and therefore the continuous forms of the drag law is in demand to correctly simulate liquid–solid flows.

To avoid discontinuity of these two correlations, a switch function φ is introduced by Huilin and Gidaspow [7] to give a smooth transition from the dilute regime to the dense regime.

$$\phi =\frac{\text{arctan}\left[150\times 1.75\left(0.2-{\epsilon }_s\right)\right]}{\pi }+0.5$$E24

Thus, the interphase momentum transfer coefficient becomes

$$\beta =\left(1-\phi \right){\beta }_{\mathit{\text{Ergun}}}+{\mathit{\phi \beta }}_{\mathit{Wen}\&\mathit{Yu}}$$E25

Syamlal et al. [10] proposed a correlation to calculate the momentum transfer coefficient based on the experimental data of particle sedimentation from Garside and Al-Dibouni [30]; the following equation is usually adopted:

$$\beta =\frac{3}{4}\frac{\left(1-{\epsilon }_{\mathit{l}}\right){\epsilon }_{\mathit{l}}{\rho }_{\mathit{l}}\left|{\mathit{u}}_{\mathit{l}}-{\mathit{u}}_{\mathit{p}}\right|}{V_r^2{d}_{\mathit{p}}}{\left(0.63+4.8\sqrt{\frac{V_r}{{\mathit{Re}}_{\mathit{p}}}}\right)}^2$$E26

2.4. Virtual mass force and lift force

In gas fluidization, virtual mass force and lift force could be neglected with respect to the much smaller density of gas phase, compared with solid phase. However, in liquid–solid fluidization, these forces could not be ignored due to the fact that liquid density is usually in the same order of magnitudes as the solids.

The virtual mass force can be expressed with Ishii and Mishima [31] correlation:

$${\mathit{F}}_{\mathit{vm}}={\epsilon }_s{\rho }_l{C}_{V,d}\left(\frac{\mathrm{d}{\mathit{u}}_l}{\mathrm{d}t}-\frac{\mathrm{d}{\mathit{u}}_s}{\mathrm{d}t}\right)$$E27

where $C_{V,d}=0.5$.

And Drew et al. [32] expression is used to calculate the lift force:

$${\mathit{F}}_{\mathit{lf}}={\epsilon }_s{\rho }_l{C}_{L,d}\left({\mathit{u}}_s-{\mathit{u}}_l\right)\times \left(\nabla \times {\mathit{u}}_l\right)$$E28

where $C_{L,d}=0.5$.

2.5. Boundary conditions

The aforementioned governing equations are numerically solved with appropriate boundary and initial conditions. Initially, solid particles are packed in the bed with a fixed solid concentration, and there are no motions for both the liquid and solid phases in the bed. At the inlet, the liquid velocity is constant with the concentration of unity. The particle velocity and the granular temperature are set to be zero. At the top of the bed, Neumann boundary conditions are applied to both the liquid and solid phases, and the liquid pressure is 1 atm.

The no-slip condition is set to the liquid phase at the wall. For the solid phase, the normal velocity is also set to be zero. The following boundary equations are applied for the tangential velocity and granular temperature of solid particles at the wall [33]:

$$u_{t,w}=-\frac{6{\mu }_s{\epsilon }_{s,\text{max}}}{{\mathit{\pi \rho }}_s{\epsilon }_s{g}_o\sqrt{3\theta }}\frac{\mathit{\partial }{u}_{s,w}}{\mathit{\partial n}}$$E29

$${\theta }_w=-\frac{k_{\mathit{se}}\theta }{e_w}\frac{{\mathit{\partial \theta }}_w}{\mathit{\partial n}}+\frac{\sqrt{3}{\mathit{\pi \rho }}_s{\epsilon }_s{u}_s{g}_o{\theta }^{3/2}}{6{\epsilon }_{s,\text{max}}{e}_w}$$E30

This simulation is carried out with the CFD codes and incorporates with kinetic theory of granular flow. To solve difference equations obtained from the differential equations, the second-order Total Variation Diminishing method (TVD) scheme is used. In all simulations, the time step is set to be constant of 1.0 × 10⁻³ s, such time step have been validated to be suitable for the simulation of liquid solid flows in such a fluidization system. Time-average distributions of flow variables are computed covering a period of 100 s corresponding to 2–3 weeks of computational time on a personal computer, and the last 40 s results are selected to make the average.

The solid and liquid phases treated as fully interpenetrating continua based on the extended granular flow theory; thus, the TFM is used to simulate the two-phase flows by FLUENT 6.3 based on the finite volume method. Since three-dimensional simulations are currently out of reach practically with the consideration of computation power, the simulations are carried out in two-dimensional rectangular Cartesian coordinates by ignoring front and rear wall effects.

Detailed parameter values for the simulation as well as in the experiments are reported in Table 1.

Parameters	Value	Parameters	Value
Bed height (m)	2.04	Particle diameter (m)	3 × 10^−3
Bed width (m)	5.1 × 10^−2	Particle density (kg/m^3)	2500
Initial bed height (m)	0.35	Terminal velocity (m/s)	0.318
Superficial liquid velocity (m/s)	0.093, 0.175, 0.247	Restitution coefficient of particles	0.99, 0.9, 0.8, 0.6
Maximum packing concentration	0.63	Restitution coefficient of wall	0.99, 0.9, 0.8, 0.6
Liquid viscosity (Pa-s)	1.0 × 10^−3	Liquid density (kg/m^3)	998

Table 1.

Numerical simulation parameters.

3. Results and discussion

3.1. Effect of drag force models

Figure 1 compares the volume fraction and the axial solid velocity against dimensionless radial position, using different drag models proposed by Gibilaro et al. [34], Huilin Gidaspow [7], Symalal-Obrien [10], and Wen-Yu [9], respectively. Experimental data are reported by Limtrakul et al. [35] From Figure 1, we can see that simulations by using Huilin-Gidaspow [7] model and Syamlal-Obrien [10] model are consistent with experimental data in axial solid velocity. The others deviate too much from the experimental data. Comparing the volume fraction between Huilin-Gidaspow [7] model and Syamlal-Obrien [10] model, the latter are much more different from experimental data than the former at the dimensionless radial position between 0.6 and 0.9. Hence, we chose Huilin-Gidaspow [7] drag model for further study.

Figure 1.

Influence of different drag models on (a) volume fraction; (b) axial Solid Velocity at a superficial liquid velocity of 0.07 m/s.

3.2. Effect of virtual mass

Figure 2 shows various volume fraction against dimensionless radial position at three different bed heights of H = 0.27, 0.43 and 0.60 m. In this simulation, the superficial liquid velocity is 0.1 m/s. It is obvious that volume fraction without virtual mass is higher than volume fraction with virtual mass in the whole column. Figure 3 shows various axial solid velocity against dimensionless radial position at three different bed heights of H = 0.27, 0.43 and 0.60 m. It is significant to notice that axial solid velocity without virtual mass is higher than solid velocity with virtual mass in the central region, lower near the wall. When the particle moves with accelerated motion in the fluid, it will accelerate the surrounding fluid in return. Due to inertia of fluid, fluid will give particle a reactive force. At this time, the reactive force will be greater than inertia force of particle itself, as if particle quality is increased, the bed voidage will also increase. This explains why volume fraction with virtual mass has a lower value. From Figure 3, we can see that particles go up in the central region and fall down close to the wall. Reactive forces from particles hinder the movements of fluid at the centre and reinforce the movements near the wall. So the value of axial solid velocity is lower at the central region and higher close to the wall when considering virtual mass.

Figure 2.

Influence of virtual mass on volume fraction at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

Figure 3.

Influence of virtual mass on axial solid velocity at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

3.3. Effect of lift

Figure 4 shows various volume fractions against dimensionless radial position at three different bed heights of H = 0.27, 0.43 and 0.60 m at a superficial liquid velocity is 0.1 m/s. In the simulation, lift coefficient between particles and liquid is 0.5. From Figure 7, it is clear that volume fraction with lift has a little higher of volume fraction than that without lift at the central region. The phenomenon is opposite near the wall, where the volume fraction of particles is lower when the lift force is considered. In general, solid volume fractions in both conditions have not much difference, not like the effect of virtual mass, and they have the similar trend. From Figure 5, we can see that axial solid velocity with lift is almost higher than the axial solid velocity without lift in the whole radial direction at bed height of 0.27 m. And axial solid velocity without lift is higher than solid velocity with virtual mass in the central region, lower near the wall at bed heights of 0.43 and 0.60 m.

Figure 4.

Influence of lift on volume fraction at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

Figure 5.

Influence of lift on axial solid velocity at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

3.4. Effect of liquid velocity

We can see that the volume fraction decrease with the increase in superficial liquid velocity obviously in Figure 6(a). This is because at higher liquid velocity, it will enhance bed expansion and the bed voidage will increase, leading to decreased volume fraction. The Figure 6(b) shows that axial solid velocity increases with the superficial liquid velocity changing from 0.07 to 0.10 m/s while the axial solid velocity at 0.13 m/s is lower than the above two velocities. Normally, the increase in superficial liquid velocity increases the energy input to the system, leading to enhanced solid motion. The deviation of axial solid velocity at 0.13 m/s may be due to serious fluctuation of energy at a high velocity.

Figure 6.

Influence of superficial liquid velocity on (a) volume fraction and (b) axial solid velocity.

3.5. Comparison of wall effects

In this work, the Johnson and Jackson [33] wall boundary condition with different specularity coefficients are used. The volume fraction and axial solid velocity are plotted in dimensionless radial position at different bed heights of H = 0.27, 0.43 and 0.60 m above the inlet. In Figure 7, between the centre and wall (about at the dimensionless radial position of 0.8), volume fraction reaches its maximum. The volume fractions are comparatively close to each other for all values of e at central region. The volume fractions increase at central region and decrease near the wall from 0.8 to 1. From Figure 8, we can see that axial solid velocities are similar for all values of e at bed height of 0.27 m but obviously different from each other at bed height of 0.43 m. It indicates that the behaviour of fluid is very complicated and difficult to have similar characters at middle of column. Besides, the trend of axial velocity by e = 1 seriously deviates from others in central region at bed height of 0.60 m. When volume fraction reaches its maximum, axial solid velocity reaches its minimum. And it is the coupling effect of upward fluid and downward fluid that gathers the particles together near the wall, and the volume fraction reach its maximum.

Figure 7.

Wall effects on volume fraction at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

Figure 8.

Wall effects on axial solid velocity at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

3.6. Comparison of restitution coefficient

In the particle collision process, particles collide and rebound with energy dissipation during the contact. The ratio between relative velocities after and before collision is the restitution coefficient. In the gas–solid flow system, restitution coefficient is usually a constant, while in the liquid–solid flows, such restitution coefficient is considered varying with the change of Stokes number of particles. In the discrete particle model, the effect of Stokes number on restitution coefficient could be well considered with a varying value; however, in the TFM, we take it as different constant values, so the effect of the relative velocity variation due to the viscosity of fluid could be taken into consideration. In this work, the restitution coefficient varies from 0.7 to 1.0; the effect of different restitution coefficients on particle volume fraction distribution is listed in Figure 9 for different heights of the bed. Different restitution coefficient can generate different particle volume fraction distribution along the radial direction. Particles with restitution coefficient of 1.0 generate the highest volume fraction in the bed. For other restitution coefficient values, the value of 0.7 generates a comparatively higher particle volume fraction distribution due to energy dissipation during colliding process and makes particle difficult to be transported by the fluid. Figure 10 shows the effect of different restitution coefficient on axial solid velocity distribution along radial position. The biggest restitution coefficient of 1.0 corresponds to the lowest axial solid velocity distribution at lower part of the bed as Figure 10(a) and (b) shows. Since particles at the lower part has higher volume fraction, particles tend to have more opportunities to collide, thus along with the velocity distribution, this tends to generate more energy dissipation. At the higher part of the bed, particles near the walls have higher volume fraction and lower velocity distribution due to the frictional resistance of the walls.

Figure 9.

Restitution coefficient on volume fraction at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

Figure 10.

Restitution coefficient of axial solid velocity at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

3.7. Comparison of radial distribution

Figure 11 shows the radial distribution of volume fraction at three different bed heights under prediction of TFM-KTGF model with different radial distribution function formations. The Syamlal-O’Brien [10] radial distribution function presents the lowest volume fraction near the walls at three different bed heights, and the Lun [21] radial distribution function model provides the highest volume fraction at most of the radial directions for higher heights in the bed. It is clear that the distribution function can generate obvious different solid volume fraction distribution; however, the distribution tendencies are the same.

Figure 11.

Radial distribution of volume fraction at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

Figure 12 shows the radial distribution of axial solid velocity at three different bed heights under prediction of TFM-KTGF model for different radial distribution function formations. The Syamlal-O’Brien radial distribution function presents the lowest axial solid velocity near the walls at three different bed heights. Particle axial velocity distribution shows a decrease of solid velocity from the centre to the near-wall region and then increase near the walls of the bed. Such distribution pattern agrees well with the previous experimental results.

Figure 12.

Radial distribution of axial solid velocity at three different bed heights of (a) H = 0.27 m; (b) H = 0.43 m; (c) H = 0.60 m at a superficial liquid velocity of 0.10 m/s.

3.8. Fluctuations at different liquid velocities

Fluctuation of solid phase mean volume fraction at three different liquid velocities of 0.07, 0.10 and 0.13 m/s at the bed height of 0.27 m between 80 and 100 s is shown in Figure 13. As one can find from Figure 13 (a) to (c), the amplitude of the fluctuation increased, and the frequency also increased with the increase of liquid velocity. Figure 14 shows the fluctuations of mean solid velocity at three different liquid velocities of (a) v = 0.07 m/s; (b) v = 0.10 m/s; (c) v = 0.13 m/s at the bed height of 0.27 m. The trend for the effect of liquid velocity on solid volume fraction fluctuations can also be found in Figure 14 for the fluctuations of the mean solid velocity.

Figure 13.

Fluctuation of solid volume fraction at three different liquid velocities of (a) v = 0.07 m/s; (b) v = 0.10 m/s; (c) v = 0.13 m/s at the bed height of 0.27 m.

Figure 14.

Fluctuation of solid velocity at three different liquid velocities of (a) v = 0.07 m/s; (b) v = 0.10 m/s; (c) v = 0.13 m/s at the bed height of 0.27 m.

3.9. Analysis of granular parameters

Granular pressure distribution as with increasing solid volume fraction for different granular pressure models is shown in Figure 15(a); granular pressure increases with the increase of solid volume fraction to a maximum value for lower volume fractions and then decreases with the increase of solid volume fraction for higher solid volume fractions. Such a distribution is because that the granular pressure is highly related to particle collision, for lower solid volume fraction, increase of solid volume fraction will generate more chance for particles to collide with each other, thus result in a higher granular pressure, while such collision reaches its maximum, the granular pressure will decrease with the increase of solid volume fraction due to collision mechanism is being hindered by more particles per unit volume, and quasi-static contact will play a more important role, thus granular pressure decrease at such solid volume fractions. It is obvious that all the granular pressure models can predict such distribution of granular pressure; however, the quantitative prediction of granular pressure distribution differs for these models. The model proposed by Lun et al. [21] get the highest value while the model of Syamlal and O’Brien [10] get the lowest simulation result. The coefficient of restitution is the ratio of the final to initial velocity difference between two particles after they collide, where one indicates a perfect elastic collision. When it is assumed that the collision is elastic, the granular pressure distribution is totally different from that inelastic collisions where the coefficient of restitution is less than one. As a result, when taking the numerical simulation, it is usually inaccurate to assume that such a liquid–solid system is elastic, since it will result in a falsehood granular pressure distribution.

Figure 15.

Granular pressure distribution for (a) different granular pressure models, (b) different coefficient of restitution.

Granular temperature, $\theta =\frac{1}{3}\left(⟨v_x^2⟩+⟨v_y^2⟩+⟨v_z^2⟩\right)$, is the mean value of the squares of fluctuating velocities at three directions. As indicated previously, with the increase of solid volume fraction, particle collision possibilities increase thus result in a higher granular temperature, while the granular temperature reaches its maximum, more particles per unit volume will bring about a quasi-static contact of particles, thus granular temperature decreases at such solid volume fractions. From Figure 16(a), we can obtain that all of the models for granular pressure predict similar trend of such granular temperature distribution, and the Lun’s [21] model provides the highest value for most of the solid volume fraction. When we change the coefficient of restitution coefficient, the elastic collision assumption will result in the highest granular temperature and the trend for such distribution is also unreasonable, since particle collisions in liquid solid flow systems are inelastic, more dissipation will be generated due to fluid drag, fluctuation, etc. As a result, it should be sensitive and careful to select the coefficient of restitution according to the flow system.

Figure 16.

Granular temperature distribution for (a) different granular temperature models, (b) different coefficient of restitution.

4. Conclusions

The TFM combined with the kinetic theory of granular flow (KTGF) is employed to investigate the hydrodynamics of particles in gas–solid as well as liquid–solid fluidized bed. A variety of models and parameters including drag models, granular pressure models, coefficient of restitution are selected when carrying out such numerical simulation. The effect of such models along with the selection of the values for such parameters is comprehensively studied in this work. Numerical investigation of the particle concentration distribution in a liquid–solid fluidized bed is carried out to study the effects of drag force models as well as virtual mass force and lift force in predicting of particle flow characteristics. Different density ratios of solid/liquid, liquid viscosity as well as superficial velocities of liquid phase are also studied to investigate the effects of operating conditions on the distribution of particle concentration.

The predicted axial particle concentration shows a nearly uniform distribution throughout the bed for the investigated particles. Different drag models exhibit various particle velocity distributions indicating that selection of the drag models should be careful. Virtual mass force and lift force should be considered due to the low solid–fluid density ratio. Distribution of granular pressure and granular temperature indicates that elastic assumption for liquid–solid fluidized bed in improper and more energy dissipation due to fluid interstitial effect should be taken into account.

Nomenclature

ε	concentration of each phases (−)
ρ	density (kg/m3)
p	thermodynamic pressure (N)
τ¯l	viscous stress tensor of liquid phase (Pa)
μt	turbulent viscosity for the liquid phase (N·s/m2)
k	turbulent kinetic energy (m2/s2)
ε	dissipation rate of turbulent kinetic energy (m2/s2)
μs	shear solid viscosity (Pa·s)
C	particle fluctuating velocity (m/s)
e	coefficient of restitution for particle-particle collisions (−)
L	characteristic length scale (m)
εs,max	maximum particle concentration at random packing (−)
μs,kin	the viscosity induced from particle fluctuations(N·s/m2)
ϕ	internal friction angle (°)
ks	conductivity of granular energy
Dls	the rate of energy dissipation per unit volume
us	velocity vector of solid (m/s)
Fvm	virtual mass force (N)
ut,w	tangential velocity at the wall (N·s/m2)
u	velocity vector (m/s)
g	gravity acceleration (N/kg)
β	inter-phase drag coefficient (−)
μf	viscosity of liquid phase (N·s/m2)
μl	laminar viscosity for the liquid phase (N·s/m2)
τ¯s	stress tensor of solid phase (Pa)
ξs	bulk solid viscosity (Pa·s)
θ	granular temperature (k)
ps	solid pressure (Pa)
λ	mean free path of particles (m)
g0	radial distribution function at contact (−)
μs,col	the viscosity induced from the particle collisions (N·s/m2)
μs,fr	the viscosity induced from particle frictions (N·s/m2)
I2D	the second invariant of deviatoric stress tensor (Pa)
γs	rate of dissipation of fluctuation kinetic energy due to particle collisions
ul	velocity vector of liquid(m/s)
φ	switch function (−)
Flf	lift force (N)
θw	granular temperature of solid particles at the wall (k)