CFD SIMULATION OF THE MIXING AND DISPERSING OF FLOATING PARTICLES IN A VISCOUS SYSTEM

Based on the Gidaspow model, the distributions of velocity, turbulence intensity, and solid concentration in stirred vessels equipped with a down-pumping propeller (TXL), a six flat-blade disc turbine (Rushton), or a downpumping six 45° pitched-blade turbine (PBTD-6) in a viscous system were simulated. The power curve of the TXL propeller and the dimensionless solid concentrations of one sampling point in the vessel at different agitation speeds were obtained by simulation and experiment, which were in good agreement with each other. The results showed that: (1) both the tangential velocity and turbulence intensity on the liquid surface caused by a Rushton turbine were the highest of the three conditions at the same agitation speed; (2) the turbulence intensity on the azimuth of 90° behind the baffle near the shaft on the liquid surface was relatively larger than that in other regions; (3) the uniformity of solid concentration distribution in the stirred vessel equipped with a Rushton or PBTD-6 turbine was better than that with a TXL impeller at the same agitation speed.


INTRODUCTION
Mixing plays a vital role in occasions such as homogenization, emulsification, polymerization and fermentation.The quantity of stirred reactors is more than 85% of the total reactors in the production process of three synthetic materials, namely synthetic plastics, synthetic fiber, and synthetic rubber (Feng and Wang, 2010;Wang and Feng, 2000).In addition to experimental research, computational fluid dynamics (CFD) is gradually becoming an important method in the study of varieties of mixing conditions with the rapid development of the finite volume method and computer technology (Rajavathsacai et al., 2014).Compared with single-phase mixing, solidliquid mixing is more complicated because of the large density difference between phases.The critical speed, critical power, local velocity and concentration can be conveniently obtained through experiments.But it is hard to obtain the flow pattern and entire concentration distribution in the vessel due to the low efficiency and high cost (Tamburini et al., 2013).Based on this, it is essential to investigate the CFD simulation of solid-liquid mixing.
So far, the simulation studies of solid-liquid mixing focused more on the suspension process of sinking particles.Fan et al.(2005) compared the solid-liquid mixing of slender particles and spherical particles (ρ s = 1125 kg•m - 3 ) by simulation, and found that the shape of particles had little effect on the velocity field in the vessel equipped with a Rushton turbine.Kasat et al. (2008) simulated the solidliquid mixing of glass particles in water with a Rushton turbine.They concluded that the mixing time increased with the agitation speed until a peak value and then gradually decreased.It became a constant after getting over the just off-bottom suspension condition until a complete suspension condition appeared and then reduced slowly.Hosseini et al. (2010) investigated the effect of impeller type, impeller off-bottom clearance, agitation speed, particle size, and particle specific gravity on the mixing quality of sinking particles through numerical simulation and experimental method.Simulation results of the impeller torque and just suspended agitation speed agreed well with experimental values.Tamburini et al. (2011;2012) reviewed the simulation methods to estimate the just offbottom agitation speed of glass ballotini particles in water, and presented a universal CFD method.Liu and Barigou (2013) investigated the effect of solid concentration on the liquid velocity field by simulation.They found that the solid concentration of coarse glass particles had little effect on the liquid velocity field when it ranged below 0.1 g/g.As the solid concentration continued to increase, the liquid velocities near the impeller and along the wall of the vessel reduced significantly.
The mixing and dispersing of floating particles is as common as that of sinking particles in the process industry.But there is a large difference between these two processes.The latter process raises the sinking particles whose density is larger than the liquid phase; the former process draws down and disperses the low-density solids floating on the surface (Khazam and Kresta, 2008;Mersmann et al., 1998).The fluctuation of the free surface and eccentric vortex make the motion of floating particles more complex.However, literature on the numerical simulation of the mixing and dispersing of floating particles is relatively scarce.Waghmare et al. (2011) developed a semi-empirical correlation to predict the drawdown rate of floating particles through a conjunction of simulation and experimental measurement.The drawdown rate is correlated to the mean liquid velocity of the free surface.Li et al. (2014) simulated both the solid-liquid mixing of floating particles and sinking particles, which had a great difference in the suspension characteristics.The concentration of sinking particles decreased along with the height while that of floating particles was the opposite.Qiao et al. (2014) compared the effects of up-pumping and down-pumping impellers on the mixing and dispersing of polyethylene particles in water.Other scholars investigated the effects of different material parameters and/or stirred vessel structures on the floating particles mixing with the aid of experiments (Guida et al., 2009;Karcz and Mackiewicz, 2007;2009;Ozcan-Taskin, 2006;Wojtowicz, 2014).
The existing researches on the mixing and dispersing of floating particles used water as the liquid phase and lost sight of the influence of liquid viscosity.However, viscous materials are commonly adopted inthe actual production.This paper studies the mixing and dispersing of floating particles in a viscous system.Three common impellers were chosen to investigate the effect of impeller type on the solid-liquid mixing by numerical simulation.

PHYSICAL MODEL AND EXPERIMENTAL METHODOLOGY
As shown in Figure 1, a stirred vessel of 380 mm inner diameter (T) and 456 mm liquid level height (H) with a standard elliptical head was adopted.A single baffle of 38 mm breadth (B), 10 mm thickness (d), and 7.6 mm clearance to the vessel wall (c = T/50) was fitted in the vessel.A down-pumping propeller (TXL), a six flatblade disc turbine (Rushton), or a down-pumping six 45° pitched-blade turbine (PBTD-6), as shown in Figure 2, corresponding to axial flow, radial flow and mixed flow impellers, respectively, were adopted in the investigation.The three impellers have the same diameter (D) of 200 mm, agitation speed (N) of 150 rpm, and impeller submergence of 95 mm (S = T/4).The vessel structure and agitation speed were determined on account of previous experiments (Chen, 2015).Malt syrup and polyethylene particles were selected as the liquid and solid phases, whose physical parameters are listed in Table 1.
To certify the reliability of simulations, the corresponding experiments were conducted with the same structure, agitation speed, and working medium with simulations.The agitation speed (N) was adjusted by a frequency converter.The shaft torque (M) and system viscosity (μ) were measured by a TQ-660 torque sensor and digital viscometer, respectively.The power number (N P ) and Reynolds number (Re) were calculated by equations (1) and (2).

(
) ( ) The solid concentrations were sampled by a peristaltic pump.

Governing equations
The common models for the simulation of solid-liquid mixing include the Euler-Lagrange model and the Euler-Euler model.The Euler-Lagrange model regards the liquid phase and solid phase as the continuous phase and dispersed phase, respectively.The continuous phase is treated in an Eulerian framework, and the motion of the dispersed phase (1) (2)    is simulated by solving the force balance of each particle.
The approach based on this model can obtain more details on the trajectory of the particles when the solid concentration is low.However, with increasing volume fraction of the dispersed phase, the interaction between the two phases increases.It will consume a large sum of computer resources by tracing a large number of particles (Ochieng et al., 2010;Shah et al., 2015).To reduce the computational cost, most researches adopt the Euler-Euler model in the simulation of solid-liquid mixing (Fan et al., 2005;Kasat et al., 2008;Li et al., 2014;Liu et al., 2013;Hosseini et al., 2010;Qiao et al., 2014;Tamburini et al., 2011).The Euler-Euler model is also known as the twofluid model, which regards both the solid and liquid phases as continuous phases.It is widely used for solid-liquid, gas-liquid and gas-liquid-solid mixing simulation for its own advantages.The simulation results obtained by this method agree well with the corresponding experimental data (Hosseini et al., 2010).Hence, the Euler-Euler model was adopted to simulate the solid-liquid mixing of floating particles in this paper.
In the current simulation, the influence of temperature and reaction was neglected.The governing equations include the continuity equation and momentum conservation equation.

Continuity equation ( ) (
) where i denotes the phase, i = s for the solid phase and i = l for the liquid phase; α i , ρ i , and i v  are the volume fraction, density, and velocity vector for each phase, respectively; p is the pressure for all the phases;  ̿ is the unit stress tensor; i F  is the external body force; ̿ � is the stress-strain tensor; ij R  is the interaction force between phases; μ i and λ i are the shear viscosity and bulk viscosity for each phase (Jiang and Huang, 2010).
For the solid-liquid two-phase flow, especially in the viscous system, the drag force between phases plays a major role compared with the force between particles (Jiang and Huang, 2010).Hence, equation ( 6) was taken to calculate the interaction force between phases.

( )
where K ij is the momentum exchange coefficient between phases; f is the function of drag force; τ s is the particle relaxation time; d s is the diameter of solid particles.
The common drag models used to calculate the momentum exchange coefficient K ij include the Syamlal-O'Brien model, the Wen-Yu model, the Gidaspow model, and so on.The Gidaspow model is the combination of the Wen-Yu model and the Ergun equation (Jiang and Huang, 2010).Based on the previous numerical tests and pertinent literature, the Gidaspow model was taken to calculate the momentum exchange coefficient (Chen, 2015).When α l is greater than 0.8, the momentum exchange is calculated by equation ( 9) as below: (3) Here the drag coefficient C D is calculated by equation ( 10).
( ) The relative Reynolds number Re s is defined as below: When α l does not exceed 0.8, the momentum exchange coefficient is calculated by equation ( 12): Solid phase pressure estimated for the granular flows (whose volume fraction is lower than the maximum allowable value) used to calculate the term p s is composed of a kinetic term and a second term as follows (Jiang and Huang, 2010): The shear viscosity μ s and bulk viscosity λ s of the solid phase in equation ( 5) can be calculated as follows: ,col ,kin s s s where Θ s is the granular temperature; e ss is the restitution coefficient for particle collision; g 0,ss is the radial distribution function; μ s,col and μ s,kin are the collisional and kinetic viscosity of the solid phase.The granular temperature conservation equation based on the kinetic theories takes the following form (Jiang and Huang, 2010): where �− �  ̿ + ̿ � �∇ ⃑ � is the generation of energy by the solid stress tensor; ( ) The term ls Φ is represented as follow: Comparing with the other two extensions of the standard k-ε turbulence model (dispersed and per-phase), the mixture k-ε model has less computational demand and better representation of the solid distribution (Qiao et al., 2014).Hence, it was adopted as the turbulence model for simulation in this paper.The equations are listed as below.
where k is the turbulent kinetic energy; ε is the turbulent kinetic energy dissipation rate; μ t,m is the turbulent viscosity; G k,m is the generation of turbulence in the mixture.The mixture density (ρ m ) and mixture velocity ( ) are calculated by equations ( 23) to (24) (Jiang and Huang, 2010).

Simulation method
The commercial CFD software FLUENT was utilized for the numerical simulation of the solid-liquid mixing of floating particles.The multiple reference frame method (MRF) was applied to calculate the flow field, which divided the whole vessel into moving area and static area.The moving area consisted of a central cylinder region with the entire impeller area, which used a rotating reference frame with the agitation speed.Other regions were defined as the static area with static reference frame.No-slip boundary conditions with the standard wall function were imposed on the solid walls for the liquid phase.Partialslip boundary condition (Johnson and Jackson, 1987) was considered for the solid phase.The solid walls included the vessel wall, the surface of the baffle, impeller and shaft.The symmetry boundary condition was imposed on the free liquid surface.The interface boundary condition was imposed on the interface between the moving and static areas.
The second order upwind scheme was used for the momentum, turbulent kinetic energy, and turbulent dissipation rate.The phase coupled SIMPLE algorithm was adopted for the pressure-velocity coupling.The time step was set to 0.001 s.The residual criterion for convergence of the flow field and concentration field were 10 -6 and 10 - 8 , respectively.The addition of a second phase employed the method of patch in initialization.The floating particles with an initial volume fraction of 0.05 L/L were uniformly distributed in the liquid phase.

Mesh generation
The models were meshed by the CFD pre-processing software GAMBIT.The moving area around the impeller and the static elliptical head region were meshed with the unstructured tetrahedral element because of their irregularity.Other regions were meshed with the structured hexahedral element.Local grids around the impeller were refined in addition.The specific mesh is shown in Figure 3.
Establishing adequacy of computational grids is important for the numerical simulation.The appropriate mesh quantity will save the computer resource, improve the computing speed, and ensure the accuracy.Based on the grid independence test, we increased the number of grids by the factor of 2 until the turbulence intensity stop changing.The appropriate grid numbers of the three stirred vessels are listed in Table 2.

Effect of boundary condition for the solid phase
Specularity coefficient (Φ) is an empirical parameter to qualify the particle-wall collisions in partial-slip boundary condition, whose range is from 0 to 1 (Li et al., 2010).The smaller coefficient corresponds to a smoother wall with less friction.It is hard to obtain the specularity coefficient through an experimental method.Data simulated on the sampling line L shown in Figure 1 in a stirred vessel equipped with a TXL impeller under four different boundary conditions for the solid phase, were compared in Figure 4.The first condition corresponded to the no-slip boundary condition, while the other three corresponded to partialslip boundary conditions with specularity coefficient of 0, 0.5, and 1.The impeller speed, liquid viscosity, and mean solid concentration were 150 rpm, 75.3 mPa•s and 0.05 L/L, respectively.It can be concluded that the effect of specularity coefficient on the liquid velocity, turbulence intensity, and solid concentration can be neglected in this study.Compared with the other region, the area near the solid wall was very small.The motions of the solid particles were controlled by the flow pattern of the liquid phase.Hence, the effect of the boundary condition for the solid phase was quite limited.In fact, most research on the simulation of solid-liquid mixing adopted the noslip boundary condition on the wall for the solid phase instead of the partial-slip boundary condition (Fan et al., 2005;Hosseini et al., 2010;Qiao et al., 2014;Tamburini et al., 2012;Waghmare et al., 2011).In the remaining sections of this paper, the partial-slip boundary condition with a specularity coefficient of 0 was adopted as the wall condition for the solid phase.

Simulation reliability verification
The power curve of the TXL impeller in a system of 27.3 mPa•s viscosity simulated with the Gidaspow model was compared with the experimental result in Figure 5.The maximum and minimum errors between simulation and experiment are 18.0% and 8.59%, respectively.
Furthermore, Figure 6 shows the dimensionless solid concentrations (C/C 0 ) of the sampling point P1 shown in Figure 1 obtained by simulation and experiment.The dimensionless solid concentrations with the agitation speed lower than 80 rpm are close to 0. Hence, the data errors are relatively large.Leaving out the first three experimental data, the maximum error between simulation and experiment was 24.3% at an agitation speed of 90 rpm.The minimum error was 1.97% at an agitation speed of 110 rpm.This verifies the Gidaspow model and other settings are suitable for the simulation of solid-liquid mixing with floating particles, which can also provide a basis and support for the follow-up simulation.

Distribution of liquid velocity field
The distributions of liquid velocity field in a system of 75.3 mPa•s viscosity in the vessel equipped with a TXL impeller, a Rushton turbine, and a PBTD-6 turbine were investigated by simulation and shown in Figure 7.In addition, liquid velocities on a radius in the three vessels are compared in Figure 8.As shown in Figure 1, the coordinate origin was located at the center of the elliptical head; the z-axis was upward along the symmetry axis of the stirred vessel.The radius is 1 mm under the liquid surface where z = 360 mm shown in Figure 7.
Figure 7 shows that the radial flow was developed around the impeller when the Rushton turbine was adopted.The fluid was divided into two streams near the vessel wall: one stream flowed upward, converged to the center in a spiral way near the surface to generate a surface vortex, and flowed down to the turbine; another stream flowed downward, converged to the center at the bottom of the vessel, and flowed upward along the axis.When a PBTD-6 turbine was adopted, the liquid flowed with a downward bias under the effect of blades, and was divided into two streams upward and downward near the wall.
As shown in Figure 8, (1) the average axial, radial, and tangential liquid velocities on the radius near the surface in the vessel equipped with a TXL impeller were smaller than those in the vessels equipped withthe other two impellers at the agitation speed of 150 rpm.This means that higher agitation speed was required to draw down and disperse the solids with a TXL impeller.(2) The maximum tangential liquid velocities near the surface of the vessels equipped with a TXL impeller, a Rushton turbine, or a PBTD-6 turbine are 0.366 m/s, 0.777 m/s, and 0.433 m/s, respectively.The tangential liquid velocity in the vessel equipped with a Rushton turbine was nearly twice as much as those corresponding to the other two impellers.For the Rushton turbine, the normal direction of the blade was the same with the linear velocity direction of the impeller.Hence, the force of the blade acting on the liquid was stronger than that of the TXL impeller and PBTD-6 turbine.

Distribution of turbulence intensity
High turbulence intensity is beneficial to increase the rate of drawdown and dispersing for solids in the stirred vessel.The distributions of turbulence intensity caused by the TXL impeller, Rushton turbine, and PBTD-6 turbine were investigated in the system with 75.3 mPa•s viscosity, which are shown in Figure 9.
As seen in Figure 9, (1) the turbulence intensity around the impeller was larger than those of the other regions in the vessel.The turbulence intensity on the blade tip of the impeller was the largest.This is because the linear velocity ofthe blade tip was the highest.The trailing vortexes generated behind the blades would also enhance the turbulence of fluid.( 2  azimuth of 90° behind the baffle near the agitation shaft on the surface was relatively large for all three impellers in the vessel equipped with a single baffle.As affected by the baffle, the tangential flow near the wall flowed to the center radially.It was mixed with the tangential flow near the center and generated the eccentric vortex.The position of the vortex observed in experiment (Liu et al., 2016), as shown in Figure 10, was basically the same where the largest intensity appeared in the simulation (Chen, 2015).
(3) The turbulence intensity near the surface in the vessel equipped with a Rushton turbine was larger than those in the vessels equipped with other impellers.As shown in Figure 7 and Figure 8, the liquid velocity near the surface in the vessel equipped with a Rushton turbine was relatively larger.Higher velocity resulted in stronger turbulence intensity.

Distribution of solid concentration field
The distribution of solid concentration field is an important index to evaluate the performance of solid-liquid mixing.Usually, it is hard to obtain the entire concentration distribution in the vessel by experiment.However, it can be visually displayed by simulation, which is beneficial to enhance the understanding of solid-liquid mixing.The distributions of solid concentration field in vessels equipped with a TXL impeller, a Rushton turbine, and a PBTD-6 turbine were investigated by simulation.
Figure 11 shows that, (1) particles tend to accumulate on the surface around the agitation shaft in the three vessels during steady-state conditions.As seen in Figure 7, the liquid on the surface converged to the center in a spiral way with particles carried.So particles concentrated to the center of the liquid surface.(2) The uniformity of the vessel equipped with a TXL impeller was the worst of the three at the same agitation speed of 150 rpm.The literature (Chen, 2015) indicates that the just drawdown speed (N jd ) of the TXL impeller was higher than those of the Rushton and the PBTD-6 turbines.The working speed of 150 rpm was much closer to the critical speed of the Rushton turbine.The just drawdown speed of the TXL impeller in these conditions was nearly 160 rpm; consequently, it was hard to draw down and disperse the floating particles on the surface at the agitation speed of 150 rpm which was lower than the just drawdown speed.The simulation result shows that the TXL impeller needs a higher speed than the other two impellers to reach a just drawdown state, which was in consonance with the experimental results.
(3) The distribution of particles in the stirred vessel was not symmetric and had a certain offset relative to the axis.The structural asymmetry of the stirred vessel with a single baffle resulted in an eccentric vortex and flow asymmetry in the vessel.It led to an asymmetry of solid concentration in the vessel.
Solid concentrations on the sampling line L shown in Figure 1 and Figure 11 for the three impellers are compared in Figure 12.The degrees of homogeneity for the three vessels were represented by the standard deviation (σ) defined as equation ( 25) (Bohnet and Niesmak, 1979).
The standard deviations for the three vessels equipped with a TXL impeller, a Rushton turbine, or a PBTD turbine were 0.242, 0.0275, and 0.0334, respectively.It can be concluded that, at the same agitation speed of 150 rpm, the homogeneity in the axial direction in the vessel equipped with a Rushton turbine or a PBTD-6 turbine was much better than that in the vessel equipped with a TXL.

CONCLUSIONS
The distributions of velocity, turbulence intensity, and solid concentration in the stirred vessel equipped with a  single baffle and a TXL impeller, or a Rushton turbine, or a PBTD-6 turbine were compared by simulation with the help of the Gidaspow model.Three impellers of 200 mm diameter, 95 mm immersion depth and 150 rpm agitation speed were adopted in the simulations.The results showed that: (1) The power curve of the TXL impellerand the dimensionless solid concentration of the point P1 were obtained by simulation with the Gidaspow model and experiment.It can be confirmed that the Gidaspow drag model is applicable to simulate the solid-liquid mixing process of floating particles in a viscous system by comparing the simulation results andexperimental data.
(2) It is beneficial to generate an eccentric vortex on the azimuth of 90° behind the baffle near the agitation shaft on the surface when the stirred vessel was equipped with a single baffle.The turbulence intensity of this position was found to be relatively larger according to the simulations.
(3) The agitation speed of 150 rpm adopted inthe simulation was much closer to the just drawdown speed of the Rushton turbine.It did not reach the critical state of the TXL impeller.The uniformity of the vessel with a TXL impeller was the worst of the three.The critical speed of the TXL impeller was higher than those of the Rushton turbine and PBTD-6 turbine.

Figure 1 .
Figure 1.Dimensions of the stirred vessel: (a) Top view; (b) Front view; L-Sampling line; P-Sampling point.
is the diffusion of energy; s k Θ is the diffusion coefficient; s γ Θ is the collisional dissipation of energy; ls Φ is the energy exchange between the fluid phase l and the solid phase s.The term s γ Θ is represented byLun et al. (1984):

Figure 12 .
Figure 12.Comparison of the solid concentration on the sampling line L shown in Figure 1 corresponding to different impellers (N = 150 rpm, μ = 75.3mPa•s, C 0 = 0.05 L/L).
Abbreviation CFD computational fluid dynamics MRF multiple reference frame PBTD-6 down-pumping six 45° pitched-blade turbine Rushton six flat-blade disc turbine TXL down-pumping propeller

Table 1 .
Physical parameters of the materials.