College of Energy Engineering,Zhejiang University,Hangzhou 310027,China

Keywords:Coaxial mixer Mixing Two-phase flow CFD

ABSTRACT A coaxial mixer consisting of an anchor and a Rushton turbine was selected as the research object,whose solid suspension characteristics were studied with the help of Computational Fluid Dynamics(CFD)method.Based on the Eulerian-Eulerian method and modified Brucato drag model,the just-suspension speed of impeller was predicted,and the simulation results were in good agreement with the experimental data.The quality of solid suspension under different rotation modes was also compared,and the results showed the coaxial mixer operating under co-rotation mode could get the best performance,and a larger anchor speed was beneficial to solid suspension by enhancing the turbulent intensity at the bottom.Compared with the anchor,the inner Rushton turbine played a dominant role in solid suspension due to its high rotational speed,whereas an extremely high inner impeller speed would make the uniformity of solid distributions become worse.Additionally,the effects of solid phase properties were also investigated,the results revealed that the higher the overall solid volume fraction and the smaller the Shields number,the worse the performance of solid suspension,meanwhile the solid suspension was more susceptible to solid density compared with particle diameter within the same Shields number gradient.

1.Introduction

The operation of solid-liquid mixing is widely applied in industrial production,such as catalysis,metallurgy,sewage treatment and dissolution.Solid suspension characteristics in solid-liquid stirred tanks have been widely studied[1-3].However,there are many limitations in experimental studies.For instance,flow patterns of solid and liquid phases are difficult to be obtained,and the experiment cost is higher compared with that of simulations.Therefore,CFD method is used to the research of solid-liquid mixing.

Establishing a reasonable mathematical model is critical to the numerical study,especially for the drag model which has great influence on the accuracy of simulation results.Common drag models include the Ergun model[4],the Wen and Yu model[5],the Gidaspow model[6]and the modified Brucato model[7].The Gidaspow model is established based on the Wen and Yu model and the Ergun model,which is used successfully in the numerical research of solid-liquid mixing[8-10].In addition,Brucato et al.[11]put forward a novel model whose drag coefficient was the function of particle diameter and Kolmogorov length scale,and the drag model could better reflect the influence of local turbulent state on drag force.Later on,Khopkar et al.[7]reduced the proportionality constant of the Brucato model to one-tenth of the original one,and found that better simulation results could be obtained.

The just-suspension impeller speed(Njs)is one of the important parameters for assessing the performance of solid suspension.Zwietering[1]first defined the Njsas the impeller speed at which no solid particles remained still at the bottom for more than 2 s.Researchers[1,12,13]have proposed many experimental methods for determining Njs,and the CFD method is also applied to predicting Njs[14-17].Hosseini et al.[17]proposed an approach based on the relationship between the average solid volume fractions on a horizontal plane where at the height of 1 mm from the bottom and the impeller speed.In this method,the impeller speed corresponding to the intersection point between tangent lines with the largest and smallest slopes is regarded as Njs.

Coaxial mixers are consisted of an outer and inner impeller for wallscarping and fast-dispersing respectively,and can be applied into various production occasions.Some researchers have studied coaxial mixers.Pakzad et al.[18]investigated an impeller consisting of the A200 and Scaba impellers,and they found that the combination of the anchor and the inner novel impeller could get lower power consumption and larger cavern sizes than other combinations.Furthermore,Pakzad et al.[19]also studied flow patterns of a coaxial mixer in non-Newtonian fluids numerically,and with the increase of No,there would be a secondary recirculation zone in the upper part of the tank.Kazemzadeh et al.[20-22]researched the effect of rheological parameters of non-Newtonian fluids on the mixing time and power consumption,and they discovered that consistency index and speed ratio had great impact on the mixing efficiency.Liu et al.[23]adopted the Eulerian-Eulerian method and the standard k-ε model to research the gas dispersion characteristics in viscous systems.By comparing the simulation results and experimental data of local gas volume fraction,they found that there was a good agreement.Hashemi et al.[24,25]researched the gas-liquid mixing with the help of electrical resistance tomography(ERT),and results showed that PBU-Anchor could get a better performance compared with the PBD-Anchor.However,most of current researches focused on the single-phase and gas-liquid flow in stirred tanks,and the numerical study of solid-liquid mixing using coaxial mixers was few.

A coaxial mixer consisting of an anchor and a Rushton turbine was selected as the research object in the paper,and its solid suspension characteristics in viscous systems were investigated by means of CFD software—Fluent 15.0.

2.Study Object

The geometric model and corresponding dimension parameters are shown in Fig.1,which is consistent with our previous experimental research[26].The density and particle diameter of the solid phase are 2600 kg·m-3and 667 μm,respectively.The viscosity(μ)of the liquid phase is 0.38 Pa·s,and the corresponding density is 1327 kg·m-3.To verify the accuracy of the mathematical model,the Njsof the experimental research is compared with that of the numerical study.

3.Mathematical Model

3.1.Two-phase flow equations

The Eulerian-Eulerian method was selected to simulate multiphase flows.Compared with the Eulerian-Lagrangian method,the Eulerian-Eulerian method needs smaller computer resource and can be used in the conditions with high solid volume fraction,so it's used widely in the study of multiphase flows[9,27,28].In the Eulerian-Eulerian method,both liquid and solid phases are considered as the interpenetrating continua and have the same form of governing equations.In addition,the sum of all phases'volume fraction is equal to 1.The continuity equation and momentum equation are as follows:

Fig.1.Dimensional parameters of the stirred tank.

Continuity equation:

Momentum equation:

where subscripts i=l and i=s denote liquid phase and solid phase,respectively.Firepresents the effects of virtual mass force and lift force,etc.g and viare gravity acceleration and velocity,respectively.Rijis the interaction force between phase i and phase j,and p is the pressure shared by all phases.αiis the volume fraction,and it satisfies:

τiis stress-strain tensor,and it can be calculated as:

Here,μiand λiare the shear and bulk viscosity,respectively,andis the unit stress tensor.

3.2.Turbulence equations

In order to enclose the momentum equation,the general solution is to get the Reynolds stress by computing turbulent kinetic energy(k)and turbulent dissipation rate(ε).The standard k-ε mixture turbulence model was adopted to simulate flows in the stirred tank,and governing equations are written as:

where ρm,μm,vmand μt,mare the density,molecular viscosity,velocity and turbulent viscosity for the mixture respectively,μt,iis the turbulent viscosity of phase i,and they can be calculated by:

The production of turbulence kinetic energy(Gk,m)is defined as:

Here,C1ε,C2ε,Cμ,σkand σεare constants,whose values are 1.44,1.92,0.09,1.0 and 1.3,respectively.

3.3.Equation for interaction force

Interaction force is affected by friction,pressure and cohesion,and it satisfies:

In solid-liquid systems,interaction force can be expressed as:

where Ksl,vsand vlare the solid-liquid exchange coefficient,solid velocity and liquid velocity,respectively.

In the paper,the Gidaspow model and the modified Brucato model were compared to find a suitable drag model for the following calculation of solid-liquid mixing.The Gidaspow model is given by:

When αl＞0.8:

where dsis the particle diameter,μlrepresents the molecular viscosity of the liquid phase,and particle Reynolds number Repcan be expressed as:

In the modified Brucato model,Kslcan be calculated as:

where CDis the function of the particle diameter and the Kolmogorov length scale:

where λ and CDoare Kolmogorov length scale and drag coefficient in a still liquid,respectively.

Researchers[29]studied the influence of drag force,lift force and virtual mass force on the simulation results of solid volume fraction profiles,and results showed that the drag force was the most dominant,and other forces had little influence.In addition,many researches[7,8,27,30]without considering lift force and virtual mass force also obtained satisfactory results,therefore only the drag force was considered in this paper.

4.Numerical Simulation Methods

4.1.Geometry and discretization

Geometric model and mesh were generated using Gambit 2.4.6.The dimensions of geometric model are the same as that of the stirred tank shown in Fig.1.As the structure was complex,the geometric model was discretized by hybrid grids.Tetrahedral grids were generated in the region of the elliptical head and near the inner impeller,and hexahedral grids were generated in the other zones.

A mesh independence test was carried out to find a proper number of elements.Three models with element numbers of 595333,923092,and 1807931 were compared,and their maximum Equisize skew were 0.849,0.826,and 0.868,respectively.Radial positions at the height of 1 mm from the bottom and axial positions away from the wall of 100 mm were selected,and the distributions of local volume fractions(Φvl),solid radial velocity(Vr)and solid axial velocity(Va)in these positions are shown in Fig.2.

It can be observed that when the element number increases from 923092 to 1807931,the profiles of solid volume fraction and velocity have the same tendency and change little.However,in radial positions,solid volume fraction and solid radial velocity distributions change greatly when the element number increases from 595333 to 923092.Therefore,the model with element number of 923092 was adopted,which could ensure adequate calculation accuracy and reduce computing time.

4.2.Boundary conditions

The Multiple Reference Frame(MRF)method was selected to simulate the rotation of impellers,considering its successful application in the steady calculation of coaxial mixers[23,31]and less calculation compared with the sliding mesh method.The solution domain was divided into four zones as shown in Fig.3.The inner and outer moving zones were set at the same rotational speed as Rushton turbine and anchor respectively,while other zones were kept stationary,and momentum between different zones could be exchanged through interfaces.The top of stirred tank was considered as the symmetry boundary,tank wall was defined as the stationary wall,and the surfaces of impellers and shaft were set as the rotational wall.

4.3.Other simulation strategies

The SIMPLE method was used to solve pressure-velocity coupling,and the standard wall function with no slip shear condition was selected.In order to mitigate the effect of false diffusion and improve accuracy of simulation results,the high order upwind scheme was adopted.In simulations,the convergence absolute criterion was set as 10-5,and steady calculation was performed.

5.Results and Discussion

5.1.Comparison of drag models

The drag model is critical to the accuracy of simulation results.The Gidaspow model and the modified Brucato model were compared to find a suitable one for the calculation of solid-liquid mixing.Due to lack of experimental data of coaxial mixers' solid-liquid flow field,a single-impeller mixer in solid-liquid mixing systems reported by Guha et al.[32]was selected as the research object.

The axial profiles of solid radial velocity at two radial locations are shown in Fig.4,and simulation results of the modified Brucato model are more accurate than that of the Gidaspow model,which is in agreement with the conclusion of Wadnerkar et al.[33].Therefore,the modified Brucato model was used in the following study.

Fig.2.Solid volume fraction and velocity profiles:(a)radial positions;(b)axial positions(Φv=4.3%,μ=0.38 Pa·s,No=25 r·min-1,Ni=230 r·min-1).

5.2.Prediction of just-suspension impeller speed

The method proposed by Hosseini et al.[17]was adopted to predict Njs.The average solid volume fraction(Φvp1)on a horizontal plane at the height of 1 mm from the bottom was monitored with increasing inner impeller speed(Ni),and the Φvp1-Nicurve is shown in Fig.5.

Fig.3.Division of the solution domain(1 outer stationary zone;2 outer moving zone;3 inner stationary zone;4 inner moving zone).

It can be seen that Nihas no obvious influence on Φvp1when Niis relatively low(＜50 r·min-1).However,with the increase of Ni,Φvp1declines sharply and tends to be stable finally.Tangent lines with the largest and smallest slopes are selected,and impeller speed corresponding to their intersection point is Njs.The simulation value and experimental data are 212 and 227 r·min-1,respectively,and the error is 6.6%.

In order to further validate the reliability of simulation values,Njswas predicted at different Noand Φv,and comparisons are shown in Figs.6 and 7.From figures,it can be observed that the curves of simulation values match well with the curves of experimental data,and the maximum error is 15.8%,which means the method for predicting Njsis applicable,and the mathematical model is accurate.

5.3.Effect of rotation modes

The most remarkable feature of coaxial mixers is that the inner and outer impellers can operate at different rotational directions and speeds,and the solid suspension performances under four rotation modes were compared based on the cloud height,solid volume fraction and solid velocity vector distributions.

The distributions of solid volume fraction and solid velocity vectors on the vertical plane are shown in Fig.8.Obviously,the distribution of solid particles is the most uniform under co-rotation mode and the worst under single-outer-impeller mode.The results can be illustrated from the circulation flows of solid phase.Under co-rotation mode,the range of circulation loops is very large,the upper recirculation zones are near the liquid level,and the lower circulation loops are almost full of the whole bottom.As a result,the coaxial mixer has stronger ability of circulation under co-rotation mode.However,under singleinner-impeller mode,the range of circulation loops is smaller than that under co-rotation mode,especially for the upper circulation loops which are important for the movement of solid phase toward the top,and the circulation loops are further compressed by the outer impeller under counter-rotation mode.Additionally,the intensity of circulation loops under single-outer-impeller mode is very weak due to its low impeller speed,causing solids accumulation at the bottom.

Fig.4.Axial profiles of solid radial velocity:(a)r/R=0.33;(b)r/R=0.5.

Fig.5.Φvp1varies with Ni(Φv=4.3%,μ=0.38 Pa·s,No=25 r·min-1).

Fig.6.Comparisons between the simulation and experimental results at different No(Φv=4.3%,μ=0.38 Pa·s).

The cloud height is defined as the height of the interface between solid-rich and clear liquid zones.In simulations,the cloud height can be obtained through the iso-surface with a threshold,and the overall solid volume fraction is selected as the threshold as suggested by Kasat et al.[30].From Fig.9,it can be observed that the relative cloud height(Cloud height/Liquid level height)is the highest under corotation mode and lowest under single-outer-impeller mode.In addition,the relative cloud height under single-inner-impeller mode is far higher than that under single-outer-impeller mode,and the reason is that the inner impeller as the fast-dispersing impeller has a higher rotational speed than the outer impeller used for wall-scraping.

In summary,the coaxial mixer can obtain the best performance of solid suspension under co-rotation mode,which was adopted in the following study.Additionally,by comparing the performance under single-inner-impeller and single-outer-impeller modes,it can be concluded that the inner impeller plays a more important role in solid suspension than the outer anchor.

Fig.7.Comparisons between the simulation and experimental results at different Φv(No=25 r·min-1,μ=0.38 Pa·s).

Fig.8.Solid volume fraction and solid velocity vector distributions under four rotation modes:(a)co-rotation mode;(b)single-inner-impeller mode;(c)counter-rotation mode;(d)single-outer-impeller mode(Φv=4.3%,μ=0.38 Pa·s).

Fig.9.Cloud height under different rotation modes(Φv=4.3%,μ=0.38 Pa·s).

5.4.Effect of outer impeller speed

The simulations with outer impeller speed ranging from 10 to 25 r·min-1were carried out,and the outer impeller speed was low considering the energy consumption of the outer impeller would be great at a high rotational speed out of its large impeller diameter.

Under different outer impeller speed,solid volume fraction and turbulent kinetic energy distributions on the horizontal plane located 1 mm from the bottom are shown in Fig.10.The results reveal that with the increase of No,both the area of high solid concentration zones and the maximum value of solid volume fraction decrease,besides,the turbulent kinetic energy also increases significantly.Over the entire range of No,it can be found there is the lowest turbulent kinetic energy in the center,and this is because the radial distance of the center position is the shortest,which leads to a lower flow velocity.By comparing solid volume fraction and turbulent kinetic energy distributions at the same No,it can be found that the area with high solid volume fraction corresponds to the region with low turbulent kinetic energy.Therefore,it can be concluded that the outer impeller can improve the suspension performance by enhancing the turbulence intensity at the bottom since its lower blade is very close to the tank bottom.

The solid velocity distributions on the vertical line located 100 mm away from the tank wall are shown in Fig.11.It indicates that the solid velocity is controlled by the inner impeller and has the maximum value at the height of the inner impeller,besides,the maximum values are almost the same out of the same inner impeller speed in the four cases.However,solid velocity in the zone away from the inner impeller is significantly affected by No,especially for the region above the inner impeller,and the larger No,the larger solid velocity.

5.5.Effect of inner impeller speed

In solid-liquid mixing of coaxial mixers,the inner impeller plays an important role.Therefore,simulations were carried out to investigate the effect of inner impeller speed on solid suspension,especially at high rotational speed.

The solid volume fraction distributions on the vertical plane are demonstrated in Fig.12,and local solid volume fractions in the same line as Fig.11 are shown in Fig.13.It's obvious that with the increase of Ni,the local solid volume fraction at the bottom reduces,while increases at the upper part of the tank.Additionally,solid distribution on the vertical line is very uniform except the bottom and top of the tank,even when Niis low.As the speed continues to increase,the solid distributions change slowly when Niis more than 500 r·min-1.

The relationship between the cloud height and inner impeller speed is demonstrated in Fig.14.When the inner impeller stops running,the cloud height is very low only with the rotation of the outer impeller,and the cloud height grows rapidly as the inner impeller begins to rotate.According to the 90%cloud height criterion[34],the impeller speed can be regarded as Njswhen the relative cloud height reaches 0.9.However,the critical value determined by the criterion is smaller than the result obtained through the method which is discussed in Section 5.2,and the latter is closer to the experimental value.The reason is that the cloud height of coaxial mixers is higher than that of singleshaft impellers at the critical solid suspension condition out of the rotation of the outer impeller and the removal of baffles.

To further characterize the effect of Nion solid suspension,the standard deviation(σ)of Φvlon the vertical plane was calculated by Eq.(23):

The smaller the σ,the more uniform the solid volume fraction distributions.As demonstrated in Fig.15,the relationship of σ and Niis similar to the relationship shown in Fig.5.However,along with the Niincreasing,σ will rise slowly when Niis extremely high(＞800 r·min-1),which means the uniformity of solid distributions gets worse.Therefore,the extremely high inner impeller speed can't lead to a better performance of solid suspension,and similar results are also found in single-shaft impeller systems[35,36].When Nireaches Njsof the situation,the value of σ is still great.The reason is that when the liquid level height is high,although the just off-bottom solid suspension condition has been reached at the bottom,there still exists a larger difference of local solid volume fraction between the top and bottom of the tank.As a result,various parameters should be considered to judge the quality of solid suspension,such as cloud height and local solid volume fraction.

5.6.Effect of overall solid volume fraction

Under different overall solid volume fractions,the solid phase distributions and their standard deviation on the vertical plane are demonstrated in Fig.16.Obviously,the uniformity of solid phase distributions decreases with the increase of Φvat same impeller speed.This is because the higher the Φv,the greater the apparent density of phases,and more difficult for the impellers to promote flows of solid and liquid phases,causing non-uniform distributions of solids.

The average solid volume fraction on the horizontal plane was simulated,and the ratio of the average solid volume fraction to solids packing limit was calculated,whose profile is shown in Fig.17.At the same impeller speed,the value of maximum solid concentration on the plane is greater with higher overall solid volume fraction,and the ratio increases rapidly,which means more power consumption is needed to reach the just off-bottom suspension condition.

The variation of cloud height with solid holdup is shown in Fig.18,and there is no obvious change in cloud height under different overall solid volume fractions.The simulated results indicate that the cloud height is insensitive to solid concentration when the impeller speed is fixed.However,it should be noted that although the impeller speed is the same,the energy input is larger at high solid concentrations out of the increasing apparent density.

5.7.Effect of particle diameter and density

The properties of solid phase including solid diameter and density have a significant impact on solid suspension,and the Shields number(θ)which represents the ratio of the force suspending particles to the gravity of particle in flows was selected as the research parameter.

Fig.10.Solid volume fraction and turbulent kinetic energy distributions at different No(Φv=4.3%,μ=0.38 Pa·s,Ni=230 r·min-1).

Fig.11.Distributions of solid velocity in the vertical line at different No(Φv=4.3%,μ=0.38 Pa·s,Ni=230 r·min-1).

The equation[37]for calculation is as follows:

Fig.13.Vertical profiles of local solid volume fraction at different Ni(Φv=4.3%,μ=0.38 Pa·s,No=25 r·min-1).

where N and D represent rotational speed and impeller diameter respectively,and the inner impeller diameter and speed are selected for the calculation.

Fig.14.Cloud height at different Ni(Φv=4.3%,μ=0.38 Pa·s,No=25 r·min-1).

Under different solid diameters and densities,the relations between the Shields number and solid suspension performance are shown in Figs.19 and 20.The results show that at the same rotational speed,average solid volume fraction on horizontal plane and the standard deviation of local solid volume fraction on the vertical plane all decrease with the increase of the Shields number.From Eq.(24),it can be found that the smaller Shields number corresponds to the larger particle diameter and density,which leads to a greater deposition rate.Therefore,the smaller Shields number is detrimental to solid suspension.In addition,within the same Shields number gradient,the descent velocity of average solid volume fraction and standard deviation under different solid densities is greater than that under different particle diameters,which means that within the same Shields number gradient,the influence of solid density on suspension is greater compared with that of particle diameter.

6.Conclusions

The solid suspension characteristics of the coaxial mixer consisting of a Rushton turbine and an anchor were studied with the help of the CFD method.The Eulerian-Eulerian method,the standard k-ε mixture turbulence model and the modified Brucato model were adopted to simulate the solid-liquid mixing in the stirred tank.The main conclusions are as follows:

Fig.15.σ of solid volume fraction distributions varies with Ni(Φv=4.3%,μ=0.38 Pa·s,No=25 r·min-1).

Fig.16.Effect of Φvon standard deviation(μ=0.38 Pa·s,Ni=230 r·min-1,No=25 r·min-1).

1)Compared with the Gidaspow model,the modified Brucato model was found to be more accurate and suitable for the simulation of solid-liquid mixing.With the modified Brucato drag mode,the just-suspension impeller speed was predicted,and there was a good agreement between simulation results and experimental data.

2)Under co-rotation mode,the coaxial mixer could achieve the best solid suspension performance,and a large outer impeller speed was beneficial to solid suspension by enhancing the turbulent intensity at the tank bottom.

3)The inner Rushton turbine played a dominant role in solid suspension,and the increasing inner impeller speed could quickly improve the performance of solid suspension,whereas an extremely high inner impeller speed would make the uniformity of solid distributions become worse.

4)The high overall solid volume fraction and small Shields number were detrimental to solid suspension.Besides,it was found that the cloud height was insensitive to solid concentration at the fixed rotational speed,and the solid suspension performance was more susceptible to solid density compared with particle diameter within the same Shields number gradient.

Fig.17.Effect of Φvon horizontal average solid volume fraction(μ=0.38 Pa·s,Ni=230 r·min-1,No=25 r·min-1).

Fig.18.Effect of Φvon cloud height(μ=0.38 Pa·s,Ni=230 r·min-1,No=25 r·min-1).

Fig.20.Effect of Shields number on standard deviation:Group a ds=0.667 mm,ρs=2200-4000 kg·m-3;Group b ρs=2600 kg·m-3,ds=0.5 mm-1.2 mm(Φv=4.3%,μ=0.38 Pa·s,Ni=230 r·min-1,No=25 r·min-1).

Nomenclature

CDDrag coefficient

CDoDrag coefficient in a still liquid

dsParticle diameter,m

Gk,mProduction of turbulence kinetic energy,kg·s-3·m-1

g Gravity acceleration,m·s-2

H Height of liquid level,mm

h Height of monitoring point,mm

KslSolid-fluid exchange coefficient

k Turbulent kinetic energy,m2·s-2

N Impeller speed,r·min-1

NiInner impeller speed,r·min-1

NjsJust-suspension impeller speed,r·min-1

NoOuter impeller speed,r·min-1

p Pressure shared by all phases,Pa

RijInteraction force between phase i and phase j,N·m-3

RepParticle Reynolds number

r Radial position,mm

Fig.19.Effect of Shields number on horizontal average solid volume fraction:Group a ds=0.667 mm,ρs=2200-4000 kg·m-3;Group b ρs=2600 kg·m-3,ds=0.5 mm-1.2 mm(Φv=4.3%,μ=0.38 Pa·s,Ni=230 r·min-1,No=25 r·min-1).

VaSolid axial velocity,m·s-1

VrSolid radial velocity,m·s-1

VtipImpeller tip speed,m·s-1

viVelocity of phase i,m·s-1

z Axial position,mm

αiVolume fraction of phase i

ε Turbulent dissipation rate,m2·s-3

θ Shields number

λ Kolmogorov length scale,m

λiBulk viscosity of phase i,Pa·s

μ Liquid viscosity,Pa·s

μiShear viscosity of phase i,Pa·s

μt,iTurbulent viscosity of phase i,Pa·s

ρiDensity of phase i,kg·m-3

σ Standard deviation

ΦvOverall solid volume fraction,%

ΦvlLocal solid volume fraction,%

Φvp1Average solid volume fraction on a horizontal plane located at the height of 1 mm from the bottom,%

Subscripts

i,j Phase number

l Liquid phase

s Solid phase

m Mixture phase

t Turbulent