Numerical Investigation of LiquidâSolid Countercurrent Fluidization in

Article pubs.acs.org/IECR

Numerical Investigation of Liquid−Solid Countercurrent Fluidization in a Solvent Extraction Column with an Expanding Tower Head Bo Liu,† Zhongyuan Li,*,†,‡,§ and Feng Xin† †

School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, China National Engineering Research Center for Distillation Technology, Tianjin 300072, China § Collaborative Innovation Center of Chemical Science and Engineering, Tianjin 300072, China ‡

ABSTRACT: A liquid−solid extraction column with an expanding tower head was constructed to study two-phase countercurrent fluidization processes. A composite solvent of pentane and hexane denoted as TU-A was used as the continuous phase, and glass beads were used for the dispersed phase. Numerical simulations were performed based on the particle trajectory model of the Eulerian−Lagrangian method to obtain the velocity field of the liquid and the concentration distribution of the solid particles inside the extraction column. The simulation method was previously confirmed by experiments. The operating zone and flooding area of the extraction column were determined by simulations. The results showed that the operating zone could be greatly increased by expanding the tower head, especially for small-particle systems. This study is very conducive to the industrial design and application of liquid−solid countercurrent extraction processes.

1. INTRODUCTION Solvent extraction is an effective technique for separating or removing some desirable species and pollutants from solid particles such as leachates, minerals, and contaminated soils. Accordingly, this technique has been widely used in many fields such as chemical engineering,1−3 metallurgy,4 pharmaceuticals,5 biochemistry,6,7 and soil remediation.8,9 Liquid−solid countercurrent fluidization is a general operation mode for the solvent extraction of solid particles because of its excellent mass- and heat-transfer performances. However, because of the inherent complexity of the liquid−solid countercurrent fluidization process, it still remains impossible to precisely describe the process by a theoretical approach. In recent years, the numerical simulation technique has become a powerful tool and an effective method for solving many complex systems, such as multiphase flow systems,10,11 and it has been extensively applied to various problems in science and engineering. Generally, modeling methods for fluidized beds can be classified in two categories: Eulerian− Eulerian (EE) and Eulerian−Lagrangian (EL) approaches.12 In the multiphase computational fluid dynamics (CFD) model, which is an EE approach, the fluid and the solid phases are considered to be continuous and fully interpenetrating, and the two phases are described in terms of separate sets of conservation equations with appropriate interaction terms representing the coupling between the phases. In the discrete phase method (DPM), which is an EL model, the threedimensional motion of each individual spherical particle is directly calculated from the forces acting on it, accounting for the interactions between the particle and the fluid phase. The fluid dynamic model of the continuous phase is based on the volume-averaged Navier−Stokes equations. In previous works, the EE method, especially the EE granular multiphase approach, has been widely applied to simulate the dynamical characteristics of gas−solid13−16 and liquid− solid17−19 fluidized beds. The EL method, as another widely © 2014 American Chemical Society

used approach, has also been extensively employed in simulations of fluidized beds.20,21 Li et al.22 proposed a new approach for predicting the characteristics of discrete phases of three-phase flow, in which the gas−liquid−solid flow in a fluidized bed was simulated by a combination of CFD with the DPM model and the volume tracking was performed by the volume-of-fluid (VOF) method. Chiesa et al.23 performed a computational study on the flow behavior of a laboratory-scale fluidized bed. The results obtained from the DPM model were qualitatively compared to those obtained from a multifluid CFD model. van der Hoef et al.24 classified previous models into five main categories according to whether a Lagrangian or Eulerian approach was used for the gas and/or particulate flow and indicated how the models could be solved numerically by finite-difference techniques. Some illustrative examples of a fluidized bed simulation were also given. Comparatively, the EL method is more frequently used in simulations of gas−solid fluidized beds. It is rarely employed for liquid−solid systems, primarily because the EL model is usually limited to relatively small hold-ups of solid particles for reasons of computational time.12 In general, liquid−solid fluidized beds can be classified as free type and constrained type depending on the operating mode. Obviously, the EL method is more accurate for small hold-ups of solid particles than the EE method in which the solid particles are considered as a continuous phase instead of as a practical discrete phase. This study was mainly focused on freely operated extraction columns under small hold-up of solid particles with a maximum of only 15%. In this case, the EL method was preferable for the liquid− solid fluidized system. Usually, to prevent the solid−liquid extraction tower from solid entrainment, the expanding tower Received: Revised: Accepted: Published: 285

March 25, December December December

2014 11, 2014 16, 2014 16, 2014 DOI: 10.1021/ie5012432 Ind. Eng. Chem. Res. 2015, 54, 285−295

Article

Industrial & Engineering Chemistry Research

Figure 1. Experimental flow scheme of countercurrent fluidization.

head is applied in practical processes. To fully understand the mechanism of the flow behavior inside an extraction tower with an expanding tower head, in this work, the DPM model was employed to simulate the liquid−solid two-phase flow inside the tower with an expanding tower head.

stabilization tank. Meanwhile, the solvent containing few particles was filtered through a filter and transferred back to the solvent storage tank for recycling utilization by a diaphragm pump. Solid particles deposited in the column bottom were continuously sent out of the system by a discharging device and recovered in a solid dryer. In the experiments, the influence of the solid and liquid operating conditions on the solid hold-up fraction was obtained by measuring the relative pressure differences at different heights using a U-tube manometer or by directly removing the solid−liquid mixture inside the column from the sampling ports. The average particle diameters of the two types of glass beads used in this work were 285 and 170 μm, and the density was 2400 kg·m−3. The particle size distribution is shown in Figure 2. It can be seen in Figure 2 that the particle size distribution conformed to a Rosin−Rammler distribution25 describing a multiple-particle-size feedstock in FLUENT. Therefore, the particles in the extraction column were described with the Rosin−Rammler distribution in the numerical simulations

2. EXPERIMENTS 2.1. Experimental Flowchart and Methods. The experimental flow scheme of the solid−liquid countercurrent fluidization process is shown in Figure 1. The extraction column was constructed with an expanding tower head. The solid particles were fed to the extraction column filled with solvent by the solid feeder at a given flow rate, and the solvent in the storage tank was pumped to the feeding stage of the column bottom, flowed from bottom to top, and then contacted with the solid particles in countercurrent mode from the column top. The liquid overflowed the extraction column through the overflow hole, and most of the fine particles entrained in the overflow were intercepted by the 286

DOI: 10.1021/ie5012432 Ind. Eng. Chem. Res. 2015, 54, 285−295

Article


Figure 2. Particle diameter distribution of the glass beads. n

Yd = e−(d / d ̅)

Ud is the superficial velocity of the particle phase and can be calculated as

(1)

where Yd represents the mass (or volume) fraction of solid particles with diameters larger than d and d̅ is the average particle diameter. The distribution parameters of two types of glass beads are listed in Table 1. The fluid is a TU-A composite

Ud =

n

285 170

6.53 6.54

Up =

−3

solvent, with a density and viscosity of 710 kg·m and 0.45 cP, respectively. Experiments were carried out at a temperature of 298.15 K and atmospheric pressure. 2.2. Particle Slip Velocity and Hold-up Fraction. The slip velocity Uslip is an important parameter charactering the countercurrent fluidization process. It is the overall velocity of downward-moving solid particles relative to the upward-flowing fluid. The slip velocity is also not the velocity of a single particle relative to the flow field around the particle, but the average velocity of a large number of particles relative to the overall flow of the fluid26 Uslip =

Ud U0 + ϕ 1−ϕ

ϕ=

(5)

Δp (ρs − ρf )gH

(6)

Here, H is the effective height of the column, Δp is the relative pressure difference between the inlet and outlet of the column, ρs is the density of the solid phase, and ρf is the density of liquid phase.

(2)

3. MATHEMATICAL MODEL 3.1. Governing Equations. In this work, the particle trajectory model of the EL method was utilized to solve the flow field of the fluid phase by integrating the Navier−Stokes equations28 while neglecting the existence of the particle phase

Qf A

Ud ϕ

The total hold-up fraction of the column is another key parameter for a liquid−solid countercurrent extraction column. It can be calculated from the total relative pressure difference of the column, which can be attained by direct measurement or by numerical simulation27

Here, U0 is the superficial velocity of the fluid phase, which can be calculated from the flow rate of the fluid Qf and the crosssectional area A of the column: U0 =

(4)

where Wd is the mass flow rate of the particles, as well as the flow rate of the discrete phase into the system in the numerical simulation. The particle velocity Up is also a key parameter, defined as the overall velocity of solid particles relative to any fixed cross section in the column

Table 1. n Values of Two Types of Glass Beads mean diameter (μm)

Wd ρs A

(3) 287


Article

Industrial & Engineering Chemistry Research Table 2. Governing Equations for the Liquid and Solid Phases Governing Equations of the Continuous Phase Continuous equation ∂(ρf uj) ∂ρf + =0 ∂t ∂xj Momentum equation

∂(ρf ui) ∂t

+

∂(ρf uiuj) ∂xj

=−

∂peff ∂xi

+

∂u ⎞ ∂ ⎛⎜ ∂ ⎛ ∂uj ⎞ μeff i ⎟⎟ + ⎟, ⎜μ ⎜ ∂xj ⎝ ∂xj ⎠ ∂xj ⎝ t ∂xi ⎠

i = 1, 2, 3

Kinetic energy equation

∂(ρf k) ∂t

+

∂(ρf kuj) ∂xj

=

∂ ⎛⎜ ∂k ⎞⎟ αkμ + G k − ρε ∂xj ⎜⎝ eff ∂xj ⎟⎠

Turbulent dissipation rate equation ∂(ρf εuj) ∂(ρf ε) εG ∂ ⎛⎜ ∂ε ⎞⎟ ε2 + = + C1ε k − C2*ερf αεμeff ⎜ ⎟ ∂t ∂xj ∂xj ⎝ ∂xj ⎠ k k peff is the effective pressure: peff = p + 2/3ρk μeff is the effective viscosity:29 d

ρ2 k εμ

( ) = 1.72

ν̂ ν 3̂ − 1 + Cν

dν ̂

ν̂ = μeff/μ Cν ≈ 100 μt is the turbulent viscosity: μt = μeff − μ Gk is the generation of turbulence kinetic energy: G k = μt

∂ui ⎛ ∂ui ⎜ ∂xj ⎝ ∂xj

+

∂uj ⎞ ⎟

∂xi ⎠

αk and αε are the inverse effective Prandtl numbers for k and ε

α − 1.3929 α0 − 1.3929

0.6321

α + 2.3929 α0 + 2.3929

0.3679

=

μ μeff

αk = αε = α, α0 = 1.0 C1ε and C2ε are model parameters: C1ε = 1.42, C2ε = 1.68

C2*ε = C2ε +

η=

Cμη3(1 − η /η0)

2SijSij k /ε ,

1 + βη3 Sij =

∂uj ⎞ 1 ⎛⎜ ∂ui ⎟ + ⎜ 2 ⎝ ∂xj ∂xi ⎟⎠

η0 = 4.38, β = 0.012, Cμ = 0.0845 Equations of Motion for Particles

du s = dt

∑ Fs +

ρs − ρf ρs

g

us is the particle velocity vector [(ρs − ρf)/ρs]g is the effective gravity per unit particle mass ∑Fs is the sum of all other forces per unit particle mass: ∑Fs = FD + Fsp + Fsa 18μ CDRe t (us ρs d p2 24

FD is the drag force: FD =

− u)

30

CD is the drag coefficient: b3Ret 24 CD = (1 + b1Ret b2) + Ret b4 + Ret b1 = exp(2.3288 − 6.4581φ + 2.4486φ2), b2 = 0.0964 + 0.5565φ, b3 = exp(4.905 − 13.8944φ + 18.4222φ2 − 10.2599φ3), b4 = exp(1.4681 + 12.2584φ − 20.7322φ2 + 15.8855φ3) Shape factor: φ = 1.0 Relative Reynolds number ρ |us − u| Ret = f dp μ Fluid velocity vector: u = u̅ + u′

ui′ = ζ ui′ 2 = ζ

2k , 3

i = 1, 2, 3 1

Fsp is the pressure gradient force: Fsp = − 6 πd p3∇p Fsa is the visual mass force: Fsa =

1 ρf d (u 2 ρs dt

− us)

288


Article

Industrial & Engineering Chemistry Research Table 2. continued Coupling between the Liquid and Solid Phases Momentum change of the control volume N

FCV = − ∑ (ṁ p, iΔt ∑ Fs) i=1

ṁ p,i is the mass flow rate of the ith particle, Δt is the time step of integration, and N is the number of particles passing through the control volume

and established boundary conditions. After that, the particle phase was introduced to calculate the forces exerted on a particle according to the flow field of the fluid phase previously obtained in the Lagrange coordinate. Successively, the tracks and velocities of the moving particles were obtained by integrating the differential equations describing the particlephase motion. The reactive force exerted on the fluid from the particles was embodied in the momentum equation of the fluid phase as a source term when the effect of the particles on the fluid phase was non-negligible. Then, the renewal flow field was obtained by solving the partial differential equations of the fluid phase with the source term. Finally, the renewal trajectories of the moving particles were obtained by further solving the differential equations of particle motion based on the renewal fluid field. These computational steps were repeated until convergence. The governing equations for the liquid and solid phases and their interactions are summarized in Table 2. 3.1.1. Modeling of the Turbulent Flow. For the system used in this work, when the system was stable, the statistical characteristics of a large number of micro groups and particles were mainly of interest, rather than the transient state of each particle or fluid micro group. In this case, the three-dimensional turbulent flow inside the fluidized column was determined by solving the steady-state Reynolds-average Navier−Stokes equation (RANS) with the proper turbulence model. Two turbulence models for the continuous phase combined with the consideration of the influence of fluid turbulence on the discrete phase were employed: the discrete random walk (DRW) model was used for predicting the motion of the particle. Specifically, the k−ε, k−ε + DRW, renormalization group k−ε (RNG-k−ε),28 and RNG-k−ε + DRW methods were selected to simulate the liquid−solid countercurrent fluidization process. By comparing the simulated results with the experiments, it was found that the RNG-k−ε + DRW method was the most precise among the four simulation approaches, as shown in Figure 3. When the solid−liquid countercurrent fluidization process neared stability, the flow regime of the continuous phase actually changed from transition flow to turbulent flow. In addition, considering the strongly disturbed influence of the solid particles on the flow field of the continuous phase, the continuous phase in the fluidization process was practically in a turbulent state. In this case, the RNG-k−ε model performed better because of its extensive flow-regime applicability (from transition flow to turbulent flow) and accuracy. Experiments also confirmed that the RNG-k−ε model for the continuous phase was appropriate, as shown in Figure 3. 3.1.2. Turbulent Dispersion of the Particles. The trajectory of each discrete-phase particle was predicted by integrating the force balance equation based on the Newtonian mechanics principle, which was derived in a Lagrangian reference frame. Forces acting on individual particles accounted for the momentum change rate per unit particle mass. Apart from the effective gravity, the viscous resistance, pressure gradient

Figure 3. Comparisons between the simulation results with different turbulence models and experiments (glass beads with dp = 285 μm).

force, and visual mass force were considered in the present study. Calculation methods for these forces are summarized in Table 2. When the effect of turbulence was negligible, the particle motion equation could be integrated with the time-averaged flow velocity u̅ solved from the RANS equations for the continuous phase. When the effect of turbulence was nonnegligible, the particle motion equation had to be calculated considering the instantaneous fluctuating velocity u′ in a stochastic way. The DRW model was utilized to compute the moving trajectory of the particle in the present study. Two basic parameters the describing DRW model are a fluctuating velocity with a Gaussian random distribution and an eddy time scale. For the k−ε series model, the fluctuating velocity could be calculated as u′ = v′ = w′ = ζ(2k/3)1/2, where ζ follows the random number of a normal distribution. The time scale represented the eddy lifetime and was calculated as31 τε = 0.6

k ε

(7)

The time for a particle to pass through the eddy is 289

31


Article

Industrial & Engineering Chemistry Research Table 3. Boundary Conditions continuous phase inlet

|u| = outlet walls

Qf (πD2 /4)

,

discrete phase

k = 0.0384 |u|2 Re−1/4 ,

ε = Cμ3/4

surface injection and |us| = 0

3/2

k 0.07D

pgauge = 0 no slip boundaries

⎡ ⎛ Le ⎞⎤ τcross = −τ ln⎢1 − ⎜ ⎟⎥ ⎢⎣ ⎝ τ |u − us| ⎠⎥⎦

escape boundaries elastic rebound

Table 4. Geometric Parameters and Physical Properties parameter

(8)

Geometry diameter of the tower body diameter of the tower head height of the tower body height of the tower head average grid width grid count grid shape Particles shape factor density mean diameter Fluid density viscosity

where τcross is the time required for the particle to pass through the eddy, τ is the particle relaxation time, and Le is the eddy length scale. The particle interacted with the fluid-phase eddy over the smaller time calculated by eqs 7 and 8. When this time was reached, a new value of the instantaneous velocity was obtained by applying the new value of ζ. 3.1.3. Coupling between Phases. For the present solid− liquid countercurrent fluidization process, there is a significant density difference between the fluid and solid particles. Opposite forces acting on the fluid by the particles would have an obvious influence on the fluid motion. The momentum transfer from each particle to the continuous phase can be computed by examining the change in momentum of the particle as it passes through the fluid control volume. The sum of the momentum change for all particles passing through a control volume over each time step appears as a momentum sink in the momentum equation of the continuous phase. This process is called coupling between the discrete phase and the continuous phase. 3.2. Boundary Conditions. The fluid flowed into the cylinder column from its bottom, and the velocity inlet condition was used. The turbulent kinetic energy and its dissipation rate at the entrance were calculated as summarized in Table 3. At the top of the column, the outlet overflow was regarded as a free surface with static pressure zero. All of the wall surfaces were set as no-slip boundary conditions. The computed particles, which represent a group of real particles with the same sizes, velocities, and displacements, were injected into the computational domain from the plane of the upper extraction column, and the distance from the plane to the top of the column was 0.05 m. Both the top export and bottom import of the column were set as the particle escaping boundaries; that is, when a computed particle arrived at the two planes, the trajectory computation was terminated, and the particle was considered to be outside the computational domain. A particle coming into contact with the wall was considered to undergo elastic rebound. 3.3. Geometrical Discretization and Numerical Strategy. 3.3.1. Geometrical Discretization. The equipment geometry and the physical properties of fluid and particles for use in the numerical simulations are collected in Table 4. The computational domain was discretized into about 40000 hexahedral control volumes with 89 quadrangle cells in the cross section of a cylinder of equal proportions to the experimental equipment. 3.3.2. Numerical Strategy. The numerical simulation was run with the commercial CFD code ANSYS FLUENT 14.0. The continuous phase was solved by a steady-state segregated solver, and the SIMPLEC algorithm was employed for the coupling of pressure and velocity. The convection terms of the

value/comment 50 mm 100 mm 2100 mm 300 mm 5 mm about 40000 hexahedral 1.0 2400 kg·m−3 170/285 μm 710 kg·m−3 0.45 cP

continuous governing equations were discretized spatially with the second-order-upwind scheme, whereas the viscous terms were discretized with the central difference scheme. All of the governing equations of the continuous phase were solved successively for 150 iterations between every two DPM iterations. The discrete phase was tracked in an unsteady way, computed particles were injected into the system each 0.1 s. After the continuous equations had been solved, the particle motion equation of each particle was integrated with the trapezoidal scheme numerically over the particle time step. The integration step was determined automatically by the software with a relative error tolerance of 10−5. The interaction between the discrete and continuous phases was considered as two-way coupling. The influence of fluid turbulence on the particle trajectory was simulated by the use of the DRW model. To increase the numerical stability, proper under-relaxation factors were used in this work, as detailed in Table 5. Table 5. Under-Relaxation Factors for the Simulation variable

under-relaxation factor

ui p k ε μt ∑Fs

0.4 0.3 0.6 0.6 1.0 0.5

The local solid hold-up ϕZ was calculated with the simulated result of the axis relative pressure gradient according to the equation ϕZ = 290

dp 1 (ρs − ρf )g dZ

(9) DOI: 10.1021/ie5012432 Ind. Eng. Chem. Res. 2015, 54, 285−295

Article

Industrial & Engineering Chemistry Research Then, the total solid hold-up was obtained by integrating this equation from the bottom to the top of the extraction column ϕ=

1 (ρs − ρf )g

H

Δp

∫Z=0 dp = (ρ − ρ )gH s

(10)

f

4. RESULTS AND DISCUSSION 4.1. Confirmation of the Simulation Method. At first, the experiments were carried out for the validation of

Figure 6. Velocity field of the fluid in the extraction column (dp = 285 μm, Qf = 0.3 m3·h−1, Wd = 3.16 × 104 kg·m−2·h−1).

Figure 4. Variation of the dimensionless slip velocity Uslip/Ut with solid hold-up fraction.

Table 6. Model Parameters of Different-Sized Particles in the Hold-up Fraction Function diameter (μm)

a

c0

c1

c2

285 170

1.697 × 10−2 −5.710 × 10−4

5.620 × 10−3 8.453 × 10−5

0.5904 0.6659

0.6020 0.2800

Figure 7. (Left) Particle concentration field and (right) fluid velocity vector field in the extraction column (dp = 285 μm, Qf = 0.3 m3·h−1, Wd = 3.16 × 104 kg·m−2·h−1).

from 0 to 0.25 m3·h−1 and from 5.09 × 102 to 4.07 × 104 kg· m−2·h−1, respectively, for 170-μm particles. The sediment processes of the particles were then simulated numerically under same operating ranges. To facilitate the description, the slip velocity was rendered dimensionless by Ut as Uslip/Ut, where Ut is the terminal sediment velocity of the particles. As a result, we obtained the “dimensionless slip velocity-particle hold-up fraction” profiles for different sizes of particles as shown in Figure 4. It can be clearly seen in Figure 4 that the simulation results agree well with the experimental results, indicating the reliability and accuracy of the simulation method used. In addition, it can also be seen in Figure 4 that the variation of the dimensionless slip velocity with the particle hold-up fraction was well-fitted by a power function empirical correlation with an average relative deviation of only 1.29% for dp = 285 μm and 0.89% for dp = 170 μm as compared to the experimental data. This function is

Figure 5. Velocity vector field of the fluid in the extraction column (dp = 285 μm, Qf = 0.3 m3·h−1, Wd = 3.16 × 104 kg·m−2·h−1).

simulation method. The total solid hold-up was measured with different liquid and solid flow rates ranging from 0 to 0.6 m3·h−1 and from 5.09 × 102 to 1.02 × 105 kg·m−2·h−1, respectively, for particles with a mean diameter of 285 μm and 291


Article


Figure 8. Particle concentration distribution on cross sections at different heights in the extraction column (dp = 285 μm, Qf = 0.3 m3·h−1, Wd = 3.16 × 104 kg·m−2·h−1).

Uslip Ut

=

ϕ+a c0 + c1(ϕ + a) + c 2(ϕ + a)2

mean diameter of 285 μm, Figure 5 shows the velocity vector distribution of the fluid in the extraction column with a liquid flow rate of 0.3 m3·h−1 and a solid flow rate of 3.16 × 104 kg· m−2·h−1. Figure 6 displays the velocity field of the fluid with the same flow rates. It can be seen that, when the particles and liquid flowed countercurrently under normal fluidizing conditions, due to the effect of particle sedimentation, some vortexes formed inside the column, which led to a downward flow near the wall and an upward flow of the fluid at the center of the column. The existence of the vortex could improve the mixing of the particles with the liquid and increase the masstransfer efficiency. Similar results were also obtained for particles with a mean diameter of 170 μm. 4.3. Concentration Distribution of the Particles. In view of the existence of the vortex in the flow field of the fluid, more particles were transferred toward the wall, resulting in a gradually increasing particle concentration along the radial direction, as shown in Figure 7.

(11)

where a, c0, c1, and c2 are equation parameters. For particles of different diameters, a, c0, c1, and c2 in eq 11 are different, as shown in Table 6. It can be clearly seen that the correlation curves also agree well with the simulation results. The average relative deviations between the correlation and the simulation are extremely small, at only 0.47% for dp = 285 μm and 0.11% for dp = 170 μm, indicating that, in practical applications, we can provide a fitted mathematical equation in terms of the simulation results to facilitate equipment design and operation. 4.2. Velocity Field of the Fluid. The flow field of the fluid and the concentration distribution of particles inside the extraction column are the two most important characteristics reflecting the fluidization state. Simulations were implemented using a group of liquid and solid velocities corresponding to the liquid and solid flow rates, respectively. For particles with a 292


Article


Figure 9. Simulated results of the operating zone without an expanding tower head.

Figure 10. Simulated results of the operating zone with an expanding tower head.

Figure 8 displays the particle concentration distribution on cross sections at different heights in the extraction column. It can be seen that the high-concentration zones of particles generally appeared near the wall of the column. The nonuniformity of the particles inside the column is mainly attributed to the uneven flow field of the fluid and the vortex. 4.4. Determination of the Operating Zone. 4.4.1. Operating Zone without the Expanding Tower Head. The escape of a computed particle from the top of the column indicates that liquid flooding takes place in the extraction column. In this work, the critical liquid flooding point is defined as the operating point when the escape rate equaled 1/100. The curve linking the critical flooding points is called the critical flooding point line. Generally, a larger escape rate implies more dramatic flooding. The operating zone was divided into two zones by the critical flooding point line: (I) the operating zone and (II) the liquid flooding zone. Figure 9 shows the simulated results of the operating zone without an expanding tower head, in which the operating zone lies above the flooding point line and the liquid flooding zone lies below the line. In Figure 9, the maxmum value, (Up/Ut)max, of the dimensionless quantity Up/Ut is actually the value of Uslip/Ut at U0 = 0. 4.4.2. Operating Zone with Expanding Tower Head. In practical applications, an expanding tower head is usually set in the top of the extraction column to abate the influence of fluctuations in the flow field on the solid−liquid fluidization operation and prevent possible particle entrainment at the column top. The expanding tower head obviously improved the

operating stability and reliability; however, the purely empirical method also often brings about great deviations from expectations. Obviously, a profound and full understanding of the effect of an expanding tower head on the solid−liquid fluidization process is indispensable and a prerequisite to the design of an adequate tower head and operation optimization. In this work, systematic numerical simulations on the effect of the tower head on the fluidization hydrodynamics of the extraction column were performed. Figure 10 presents the simulated results of the operating zones of the extraction column with an expanding tower head for two different particle diameters. Comparing Figure 9 with Figure 10, it can be clearly seen that the operating zones were markedly increased for both particle diameters. Moreover, the influence of the expanding tower head on the small-particle system was more notable than that on the large-particle system.

5. CONCLUSIONS The characteristics of liquid−solid countercurrent fluidization in an extraction column with an expanding tower head were numerically simulated. To validate the simulations, experiments were previously performed using TU-A composite solvent as the continuous phase and glass beads as the discrete phase. Both experiments and simulations were implemented at liquid and solid flow rates in the ranges of 0−0.6 m3·h−1 and (5.09 × 102)−(1.02 × 105) kg·m−2·h−1, respectively, for particles with a mean diameter 285 μm and 0−0.25 m3·h−1 and (5.09 × 102)− (4.07 × 104) kg·m−2·h−1, respectively, for 170-μm particles. The 293


Article


ui = velocity component in the ith direction, m·s−1 u′i = fluctuating velocity component in the ith direction, m· s−1 us = particle velocity vector, m·s−1 U0 = superficial velocity of the fluid phase, m·s−1 Ud = superficial velocity of the particle phase, m·s−1 Up = particle velocity, m·s−1 Uslip = slip velocity, m·s−1 Ut = terminal sediment velocity of the particles, m·s−1 Wd = particle flow rate, kg·m−2·h−1 Yd = mass (or volume) fraction of solid particles with diameters larger than d Z = tower height of a specific position, m

velocity field of the fluid and the concentration distribution of the solid particles were obtained. Simulations indicated that many vortexes can form inside the column as a result of particle sedimentation, leading to a downward flow of liquid near the wall and an upward flow at the center of the column. The operating zone and the flooding zone of the countercurrent fluidization were determined by a critical flooding line. Simulated results showed that the operating zones for particles with two different diameters were markedly increased because of the role of the expanding tower head and the influence of the expanding tower head on the small-particle system was more notable than that on the large-particle system.

■

AUTHOR INFORMATION

Greek Letters

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China (Nos. 21306129 and 41201497) and the Municipal Natural Science Foundation of Tianjin (Nos. 12JCQNJC05300 and 11JCYBJC05400).

■

NOTATION A = cross-sectional area of the column, m2 b1, b2, b3, b4 = model constants of the drag coefficient C1ε, C2ε, C2ε * , Cν, Cμ = model constants of the RNG-k−ε model CD = drag coefficient D = diameter of the extraction column, m d̅ = average particle diameter, m dp = particle diameter, m FCV = momentum change of the control volume, N FD = drag force, N·kg−1 Fs = force vector per unit particle mass, N·kg−1 Fsp = pressure gradient force, N·kg−1 Fsa = visual mass force, N·kg−1 g = constant of gravitational acceleration, N·kg−1 Gk = generation of turbulence kinetic energy due to the mean velocity gradient, J·m−3·s−1 H = effective tower height, m k = turbulent kinetic energy, J·kg−1 Le = eddy length scale, m ṁ p,i = mass flow rate of the particles, kg·s−1 N = number of particles passing through the control volume p = pressure, Pa peff = effective pressure, Pa pgauge = gauge pressure, Pa Δp = relative pressure difference, Pa Qf = flow rate of the fluid, m3·h−1 Re = Reynolds number of the continuous phase, Re = ρfU0D/ μ Ret = relative Reynolds number of the particles, Ret = ρf |us − u|dp/μ Sij = strain rate tensor, s−1 t = time, s Δt = time step of integration, s u = fluid velocity vector, m·s−1 u̅ = fluid mean velocity vector, m·s−1 u′ = fluctuating velocity vector, m·s−1

■

αk = inverse effective Prandtl number for k αε = inverse effective Prandtl number for ε β = model constant of the RNG-k−ε model ε = turbulent dissipation rate, J·kg−1·s−1 ζ = random number of the normal distribution η0 = model constant of the RNG-k−ε model μ = dynamic viscosity, Pa·s μeff = effective viscosity, Pa·s μt = turbulent (or eddy) viscosity, Pa·s φ = particle shape factor ϕ = total solid hold-up of the column ϕZ = local solid hold-up ρf = density of the liquid phase, kg·m−3 ρs = density of the particles, kg·m−3 τ = particle relaxation time, s τε = eddy lifetime, s τcross = time required for a particle to pass through the eddy, s

REFERENCES

(1) Zannikos, F.; Lois, E.; Stournas, S. Desulfurization of Petroleum Fractions by Oxidation and Solvent Extraction. Fuel Process. Technol. 1995, 42, 35−45. (2) Visser, A. E.; Swatloski, R. P.; Rogers, R. D. pH-Dependent Partitioning in Room Temperature Ionic Liquids Provides a Link to Traditional Solvent Extraction Behavior. Green Chem. 2000, 2, 1−4. (3) Jeannot, M. A.; Cantwell, F. F. Mass Transfer Characteristics of Solvent Extraction into a Single Drop at the Tip of a Syringe Needle. Anal. Chem. 1997, 69, 235−239. (4) Zhang, Y. M.; Bao, S. X.; Liu, T.; Chen, T. J.; Huang, J. The Technology of Extracting Vanadium from Stone Coal in China: History, Current Status and Future Prospects. Hydrometallurgy 2011, 109, 116−124. (5) Fages, J.; Lochard, H.; Letourneau, J. J.; Sauceau, M. Particle Generation for Pharmaceutical Applications Using Supercritical Fluid Technology. Powder Technol. 2004, 141, 219−226. (6) Gailliot, F. P.; Gleason, C.; Wilson, J. J.; Zwarick, J. Fluidized Bed Adsorption for Whole Broth Extraction. Biotechnol. Prog. 1990, 6, 370−375. (7) Davison, B. H.; Thompson, J. E. Continuous Direct Solvent Extraction of Butanol in a Fermenting Fluidized-Bed Bioreactor with Immobilized Clostridium acetobutylicum. Appl. Biochem. Biotechnol. 1993, 39−40, 415−426. (8) Li, X. G.; Du, Y. L.; Wu, G. Z.; Li, Z. Y.; Li, H.; Sui, H. Solvent Extraction for Heavy Crude Oil Removal from Contaminated Soils. Chemosphere 2012, 88, 245−249. (9) Yang, Y. W.; Li, Z. Y.; Huang G. Q.; Li, X. G. Use of Solvent Extraction and Microbial Combined Method to Remediate Oilpolluted Soils in Oil Field Area. Presented at the 3rd International Conference on Soil Pollution and its Remediation, Nanjing, China, Oct 18−21, 2008.

294


Article

Industrial & Engineering Chemistry Research (10) Kasat, G. R.; Khopkar, A. R.; Ranade, V. V.; Pandit, A. B. CFD Simulation of Liquid-Phase Mixing in Solid−Liquid Stirred Reactor. Chem. Eng. Sci. 2008, 63, 3877−3885. (11) Colagrossi, A.; Landrini, M. Numerical Simulation of Interfacial Flows by Smoothed Particle Hydrodynamics. J. Comput. Phys. 2003, 191, 448−475. (12) Reddy, R. K.; Joshi, J. B. CFD Modeling of Solid−Liquid Fluidized Beds of Mono and Binary Particle Mixtures. Chem. Eng. Sci. 2009, 64, 3641−3658. (13) Gera, D.; Gautam, M.; Tsuji, Y.; Kawaguchi, T.; Tanaka, T. Computer Simulation of Bubbles in Large-Particle Fluidized Beds. Powder Technol. 1998, 98, 38−47. (14) Benyahia, S.; Arastoopour, H.; Knowlton, M. T.; Massah, H. Simulation of Particles and Gas Flow Behavior in the Riser Section of a Circulating Fluidized Bed Using the Kinetic Theory Approach for the Particulate Phase. Powder Technol. 2000, 112, 24−33. (15) Chen, Z.; Gibilaro, L. G.; Foscolo, P. U. Two-Dimensional Voidage Waves in Fluidized Beds. Ind. Eng. Chem. Res. 1999, 38, 610− 620. (16) Taghipour, F.; Ellis, N.; Wong, N. Experimental and Computational Study of Gas−Solid Fluidized Bed Hydrodynamics. Chem. Eng. Sci. 2005, 60, 6857−6867. (17) Cheng, Y.; Zhu, J. CFD Modeling and Simulation of Hydrodynamics in Liquid−Solid Circulating Fluidized Beds. Can. J. Chem. Eng. 2005, 83, 177−185. (18) Roy, S.; Dudukovic, M. P. Flow Mapping and Modeling of Liquid−Solid Risers. Ind. Eng. Chem. Res. 2001, 40, 5440−5454. (19) Jena, H. M.; Nair, R.; Chidambaram, A. CFD Simulation of Bed Voidage Behavior of a Liquid-Solid. Presented at the National Seminar on Recent Advances in Chemical Engineering (RACE 2010), Gunupur, Orissa, India, Jan 30−31, 2010. (20) Tryggvason, G.; Bunner, B.; Esmaeeli, A.; Juric, D.; Al-Rawahi, N.; Tauber, W.; Han, J.; Nas, S.; Jan, Y. J. A Front-Tracking Method for the Computations of Multiphase Flow. J. Comput. Phys. 2001, 169, 708−759. (21) Feng, Z. G.; Michaelides, E. E. Proteus: A Direct Forcing Method in the Simulations of Particulate Flows. J. Comput. Phys. 2005, 202, 20−51. (22) Li, Y.; Zhang, J.; Fan, L. S. Numerical Simulation of Gas− Liquid−Solid Fluidization Systems Using a Combined CFD-VOFDPM Method: Bubble Wake Behavior. Chem. Eng. Sci. 1999, 54, 5101−5107. (23) Chiesa, M.; Mathiesen, V.; Melheim, J. A.; Halvorsen, B. Numerical Simulation of Particulate Flow by the Eulerian−Lagrangian and the Eulerian−Eulerian Approach with Application to a Fluidized Bed. Comput. Chem. Eng. 2005, 29, 291−304. (24) van der Hoef, M. A.; van Sint Annaland, M.; Deen, N. G.; Kuipers, J. A. M. Numerical Simulation of Dense Gas−Solid Fluidized Beds: A Multiscale Modeling Strategy. Annu. Rev. Fluid Mech. 2008, 40, 47−70. (25) Brown, W. K.; Wohletz, K. H. Derivation of the Weibull Distribution Based on Physical Principles and Its Connection to the Rosin−Rammler and Lognormal Distributions. J. Appl. Phys. 1995, 78, 2758−2763. (26) Wallis, G. B. One-Dimensional Two-Phase Flow; McGraw-Hill: New York, 1969. (27) Lapidus, L.; Elgin, J. C. Mechanics of Vertical-Moving Fluidized Systems. AIChE J. 1957, 31, 63−68. (28) Hinze, J. O. Turbulence; McGraw-Hill: New York, 1975. (29) Orszag, S. A.; Yakhot, V.; Flannery, W. S. Renormalization Group Modeling and Turbulence Simulations. Presented at the International Conference on Near-Wall Turbulent Flows, Tempe, AZ, Mar 15−17, 1993. (30) Haider, A.; Levenspiel, O. Drag Coefficient and Terminal Velocity of Spherical and Nonspherical Particles. Powder Technol. 1989, 58, 63−70. (31) Gosman, A. D.; Loannides, E. Aspects of Computer Simulation of Liquid-Fuelled Combustors. J. Energy 1983, 7, 482−490.

295


Numerical Investigation of LiquidâSolid Countercurrent Fluidization in

Recommend Documents

Numerical Investigation of LiquidâSolid Countercurrent Fluidization in