CFD simulations and experimental of homogenisation curves and mixing time in stirred liquids

Chemical Engineering Science60(2005)2427–2437/locate/cesCFDsimulations and experimental validation of homogenisation curves and mixing time in stirred Newtonian and pseudoplastic liquids Giuseppina Montante a,Michal Moštˇe k b,Milan Jahoda b,Franco Magelli a,∗a Department of Chemical,Mining and Environmental Engineering,University of Bologna,Viale Risorgimento2,40136Bologna,Italyb Department of Chemical Engineering,Institute of Chemical Technology,16628Praha,Czech RepublicReceived in revised form15October2004;accepted1November2004Available online22January2005AbstractThe aim of this work is to identify a proper Computational Fluid Dynamics(CFD)simulation strategy for the calculation of mixing time in stirred tanks forfluids characterised by either Newtonian or non-Newtonian rheological behaviour.Two different baffled tanks are considered,both stirred with multiple down-pumping45◦pitched blade turbines(PBT):a small one for extensive,preliminary simulations conducted with water and aimed at determining a generic procedure suitable with low-viscosity,Newtonian liquids;and a bigger one for testing the results with moderately viscous liquids and non-Newtonianfluids.The influence on the predicted mixing time of the tracer injection position,turbulence modelling and numerical parameters is investigated.The results of the simulations are compared with original experimental data and good agreement is obtained in all cases,once the proper mathematical models,solution procedure and relevant parameters are established.᭧2004Elsevier Ltd.All rights reserved.Keywords:Mixing time;Computationalfluid dynamics;Multiple PBT impellers;Liquid rheological behaviour1.IntroductionMechanical mixing operations are largely spread in,and often crucial,for a great variety of processes of the chemi-cal and process industry as well as in biotechnology,envi-ronmental remediation,etc.Several parameters that depict equipment behaviour are traditionally used;among these, mixing time provides a simple and a powerful means for the evaluation of mixers effectiveness.It is defined as the time needed to get a given homogeneity level starting from a non-equilibrium condition and it is typically measured after a small amount of a tracer had been injected into the stirred liquid.Phenomenological theories exist as for its mechanism (Nienow,1997)and correlations have been proposed for its prediction(Grenville and Tilton,1996;Nienow,1998;Nere et al.,2003).Over the years,several experimental techniques ∗Corresponding author.Tel.:+390512093147;fax:+39051581200.E-mail address:franco.magelli@mail.ing.unibo.it(F.Magelli). 0009-2509/$-see front matter᭧2004Elsevier Ltd.All rights reserved. doi:10.1016/j.ces.2004.11.020have been developed and employed for mixing time mea-surement.Nere et al.(2003)have recently reviewed such techniques,pointing out their main features and limitations. Some of these drawbacks such as the difficult applicability in the case of non-transparent vessels and severe experimental conditions,the confinement of information to a small num-ber of locations inside the stirred vessel and the high cost for big installations could be overcome with application of reliable CFDsimulations.The capability of CFDtools to forecast the mixing time is implicitly considered as an obvious achievement of these methods and encouraging results have been obtained so far (e.g.,Sahu et al.,1999;Bujalski et al.,2002a).However, some drawbacks are apparent in the studies that have been published in the last few years.Indeed,most of the simu-lations were performed using“black box”methods for the impeller description(Ranade et al.,1991;Patwardhan and Joshi,1999;Sahu et al.,1999;Lunden et al.,1995;Patil et al.,2001).When fully predictive simulation strategies were employed,discrepancies in the order of100%have been2428G.Montante et al./Chemical Engineering Science 60(2005)2427–2437often found between the predicted and experimental val-ues (e.g.,Jaworski et al.,2000;Bujalski et al.,2002b ).In some cases,the predicted mixing times were claimed to be severely affected by the precise point of tracer addition—a feature which was not observed experimentally.Indeed,the results obtained from the various authors do not provide a consistent description of the influence of mean flow and tur-bulence on mixing time predictions.In same cases,scarce sensitivity to eddy diffusivity is stated (Ranade et al.,1991;Patwardhan and Joshi,1999),while strong dependency on the same parameter is shown in other cases (Lunden et al.,1995).Also,most authors assert that the mean flow field is crucial for the correct determination of tracer distribu-tion,but sometimes good mixing time predictions have been obtained by using poorly simulated velocity fields (Nere et al.,2003).Moreover,mixing time predictions for non-Newtonian fluids have never been performed.Overall,a fully1st 2nd 3rd 4th (a)(b)Fig.1.Geometrical configuration of the experimental vessels:(a)T29vessel;(b)T48vessel.predictive strategy that can be confidently applied to any stirred vessel geometrical configuration has still to be iden-tified.The purpose of this paper is to further investigate the ca-pability of two of the currently available commercial CFD codes (namely,CFX and Fluent)to correctly forecast both the homogenisation process and mixing time in baffled stirred tanks and to identify a computational procedure that can be confidently applied to the case of Newtonian and non-Newtonian pseudoplastic liquids.The results of the simulations are compared with original experimental data.2.Stirred vessels and experimental conditionsThe stirred vessels considered in this investigation were two baffled tall cylindrical tanks.The smaller vessel wasG.Montante et al./Chemical Engineering Science60(2005)2427–24372429 T=29cm in diameter,H=4T height and was stirred withfour,evenly spaced6-bladed45◦PBTs(D=T/3,b/D=0.14,C1=T/3);the bigger one was T=48cm,H=3T andwas stirred with two different sets of impellers:three4-bladed45◦PBTs(D=0.40T,b/D=0.20,C1=T/2)andthree6-bladed45◦PBTs(D=T/3,b/D=0.14,C1=T/2).Experiments were conducted with water,aqueous solu-tions of polyvinylpyrrolidone(PVP,Newtonian behaviourwith viscosity equal to6and24mPa s)and aqueous shear-thinning solutions of Carbomer(TM)(pseudoplastic be-haviour),at selected impeller speeds.For these last liquids,the average apparent viscosity, aa,was calculated throughthe Metzner–Otto methodaa=k s n−1a with a=KN,(1)where k s is the consistency index, a the average shear rateand n theflow behaviour index equal to0.68in the inves-tigated conditions.For the constant K,the value of11wasadopted,as recommended for high-speed impellers.Underthese assumption the average apparent viscosity varied from25to42mPa s depending on solution concentration and im-peller speed;the modified rotational Reynolds number de-fined asRe= ND2aa(2)was in the range2500–10,000.The mixing time was evaluated from transient concentra-tion(homogenisation curves)of an electrolytic tracer mea-sured by conductivity probes.The experiments were per-formed by injecting a small amount of the process liquid saturated with KCl at the top of the vessel and measuring the resulting concentration evolution with time at selected positions in the vessel:one single position close to the ves-sel bottom in the case of the bigger vessel and four different elevations in the case of the smaller ually,the time necessary to reach95%homogeneity,the so-called t95,is referred to by researchers and practitioners and is adopted in this study too.Also,the whole homogenisation curve will be considered for stricter evaluation of the CFDsimulation procedures.Further details about the geometrical characteristics of the two stirred vessels and of the tracer injection and measure-ment positions are reported in Fig.1(a)and(b).3.CFD simulationsSince mixing time is a single parameter resulting from the treatment of transient experiment recordings,a suitable procedure had to be devised for the simulations.In order to reproduce the experiments as closely as possible,the fully developed single-phaseflowfield was calculatedfirst;later, the time evolution of the added tracer concentration super-imposed on theflowfield was determined.The simulation approach adopted for calculating the liquid flowfield was based on the use of(i)the transient“Sliding Grid”(SG)method,implemented in the commercial CFD code CFX-4and(ii)the steady Multiple Reference Frame method(MRF),implemented in Fluent6.The grid employed for the simulations consisted of432,560cells for the smaller vessel,up to1,000,000cells for the bigger one.In all cases, the whole2 azimuthal extent of the stirred vessel was considered.The Reynolds averaged Navier–Stokes(RANS) equations were selected for the modelling.The standard k– model was used for dealing with fully developed turbulence, while the k– model based on the Renormalization Group theory,known as RNG k– model,was chosen for the early turbulence regime,as it is more suitable for dealing with low Reynolds numberflows(Van den Akker,2000).For the non-Newtonian liquid simulations,the laminar contribution to the effective viscosity,to be used in the conservation and turbulent equations,was assumed to be the average,apparent viscosity value, aa,based on the Metzner–Otto method.The use of this simplified approach (Montante and Magelli,2003),instead of the rigorous rhe-ological equation for the laminar viscosity at each location of the grid,was obliged by the commercial codes that allow adopting only constant viscosity values in conjunction with the turbulent models.After checking whether the solution of theflow equations was reached,a small amount of tracer was added close to the vessel top(in the same location as the injection point of the experiments).The dynamic distribution of the tracer inside the stirred volume was then calculated by solving a Reynolds-averaged time-dependent scalar transport equa-tion,based on the assumption that the tracer is distributed in the vessel by convection and diffusion(Montante and Magelli,2004)jj t+∇ U =∇D m∇ −tt∇,(3)where is the tracer volumetric fraction,U is the mean velocity vector, is thefluid density,D m is the molecular diffusivity, t is the turbulent viscosity and t is the turbulent Schmidt number.As for the two parameters of Eq.(3),the molecular dif-fusivity was always assumed to be equal to10−9m2/s,a typical value for liquids.Actually,it was not a critical value since the contribution of molecular diffusion to the overall tracer dispersion process is negligible.Instead,more than one value was adopted for the turbulent Schmidt number. It is worth noting that one of the factors determining the magnitude of the turbulent contribution to tracer dispersion is parameter t,i.e.,the ratio between the square turbulent kinetic energy,k2,and its dissipation, ,(times an empirical constant,C ).The uncertainty on the t value reflects on the more suitable turbulent Schmidt number to adopt.In the past,different t values have been selected by the various authors,ranging from1(e.g.,Ranade et al.,1991;Patward-han and Joshi,1999)to0.2(He et al.,1999),and different2430G.Montante et al./Chemical Engineering Science60(2005)2427–2437 recommendations have been given depending on the partic-ular application(Yimer et al.,2002).Therefore,its influenceon the predictions was also investigated in this work:thedefault value of0.70has been used down to0.10.An ad-ditional discussion on the turbulent Schmidt number can befound in Montante and Magelli(2004).The standard time step adopted for the simulations wasequal to0.1s,i.e.,small enough as compared with the ex-perimental data acquisition rate(4Hz);its influence on thesimulated tracer concentration was checked by performingadditional simulations with either one order of magnitudesmaller and one order of magnitude greater time steps(fur-ther details are given below).The initial condition was of zero tracer concentration in the whole vessel volume,ex-cept for the injection region(either a large portion of the uppermost plane or a small number of cells,see below), where a volume fraction of tracer equal to one was consid-ered.As for the boundary conditions at wall,zeroflux was imposed.The solution of the resulting set of equations is a three-dimensional time-dependent map of tracer concentra-tion.For straightforward analysis and comparison with the experimental results,the time evolution of the tracer con-centration was monitored at the same position of the exper-imental probes.4.Results and discussion4.1.T29simulationsThe above-mentioned smaller vesselfilled with water and stirred with four6-bladed PBTs was simulated with the pur-pose of analysing the capability of the CFDmodels to cor-rectly forecast the mixing time value as well as the tracer concentration distribution in a Newtonian liquid during the time interval from the injection to its(almost)homogeneous distribution inside the vessel volume.The liquidflowfield produced by the stirring action at the impeller speed N=5s−1was simulatedfirst;Fluent 6was used,with MRF as the impeller simulation strategy and the standard k– turbulence model.Even though the se-lected computational strategies for single-phase stirred tanks are consolidated and their advantages and drawbacks are well known,as they have already been widely discussed in a number of previous works(e.g.Brucato et al.,1998), the evaluation of the accuracy of these results is very im-portant.Indeed,the correctness of the velocityfield plays a very important role in the prediction of the temporal and spatial distribution of a passive scalar(Ranade et al.,1991; Patwardhan and Joshi,1999;Bujalski et al.,2002b).The experimental mean axial velocity data collected by Ranade and Joshi(1989)in a single-impeller vessel were adopted for performing a quantitative evaluation of ourflow field predictions.This choice was imposed by the lack of experimental velocity data in analogous multiple impeller vessels;it can be considered acceptable because the distance-0.6-0.4-0.20.00.20.40.60.00.20.40.60.8 1.02r/TU/V tipFig.2.Experimental and simulated mean axial velocity profiles along a radius below the impellers.T=29cm vessel stirred with four6-bladed PBTs.Symbols:experimental data(Ranade and Joshi,1989);lines:simu-lations,each line corresponding to the profile calculated for each impeller.between two consecutive impellers was bigger than1.5D, thus allowing to assume that each impeller behaves inde-pendently(Mishra et al.,1994;Bittorf and Kresta,2000; Montante and Magelli,2004).In Fig.2,the experimental axial velocity profile nor-malised with the impeller tip speed,V tip,is reported along with those obtained from the simulations.The data refer to a radius midway between two baffles,1cm below each impeller(which corresponds to the experimental position z/T=0.27).As can be observed,the computed velocity pro-files below each impeller are very similar,almost superim-posed to each other,thus confirming that each impeller be-haves independently of the others;all computed profiles are in good agreement with the experimental data.The power number,as computed from the torque on the baffles,and the pumping numbers for each impeller were equal to1.7and 0.81,respectively.In the literature,experimental values of power number ranging from1.47(Ranade and Joshi,1989) to1.60(Ibrahim and Nienow,1995)and pumping number from0.77(Hiraoka et al.,2001)to0.97(Foˇr t,1986)were found for the6-bladed45◦PBTs,the differences being likely due to different experimental techniques used for obtaining these data as well as to small differences in turbine geomet-rical details.Therefore,the simulated values can be consid-ered of reasonable accuracy and,overall,the simulated liquid mean velocityfield was confirmed to be accurate enough. The N p value computed from the turbulent kinetic energy dissipation rate, ,was equal to1.06,i.e.,smaller by about 34%than that computed from torque on baffles and under-estimated with respect to the experimental one:this fact is consistent with what usually happens with the k– turbulent model because of inherent deficiency of this modelling ap-proach(Ng and Yianneskis,2000).At this stage,the experimental and computed transient dimensionless tracer concentration values can be compared.G.Montante et al./Chemical Engineering Science 60(2005)2427–2437243100.511.5050100150200t (s)c * (--)00.511.5050100150200t (s)c * (--)00.511.550100150200t (s)c * (--)2468050100150200t (s)c * (--)(a)(b)(c)(d)Fig.3.Experimental and simulated dimensionless tracer concentration vs.time.T =29cm vessel stirred with four 6-bladed PBTs.Liquid:water,N =5s −1.Symbols:experimental data;lines:Fluent results.(a)1st probe (z/H =0.06);(b)2nd probe (z/H =0.31);(c)3rd probe (z/H =0.56);(d)4th probe (z/H =0.81).The preliminary simulations refer to the single case N =5s −1and =1mPa s.In the simulations,the tracer was assumed to be injected at the vessel top into the previously calculated,fully developed flow field.At the injection instant,the tracer was distributed homogenously in the horizontal section close to the vessel top,covering the azimuthal extension of /4.The tracer concentration calculated at each time step was recorded in small volumes corresponding to the positions of the four experimental probes.The solution of Eq.(3)was obtained with Fluent 6.For the numerical solution of the passive scalar transport equation,a first order discretization scheme was adopted,a time step of 0.1s was chosen and for each time step a number of internal iterations ensuring very low (10−7)and constant values of mass residuals were performed.As for the turbulent Schmidt number,the default value of 0.70was taken.In all the cases,the concentration curves were normalised with respect to the concentration value corresponding to the fully homogeneous condition.The predicted concentration dynamic curves obtained at four different elevations along the vessel height are compared with the corresponding measured curves in Fig.3(a)–(d).As it is observed,overall,the calculated tracer profiles predict a (slightly)slower dynamics with respect to the experimental evidence.In particular,the response of the lowest probe,plotted in Fig.3(a),is not entirely satisfac-tory:a clearly slower tracer distribution is calculated with respect to the experiment.The corresponding calculated mixing time,t 95,was overestimated by 15%with respect tothe experimental value of 97s.The second probe from the vessel base gives a response (Fig.3(b))that is reproduced much more satisfactorily by the simulation;nevertheless,the predicted mixing time is higher than the experimen-tal one by 11.3%.The trend of the tracer concentration recorded at the third probe is similar for the measured and the simulated curves (Fig.3(c)),though the experimental curve is rather erratic.As in the case of transient solids dis-tribution (Nocentini et al.,2002),at heights corresponding to z/H =0.3–0.4the curves are more critical to evaluate and less reliable for parameter determination.Finally,the simulated tracer concentration curve corresponding to the probe closer to the vessel top (Fig.3(d))shows an overall slower distribution with respect to the experimental curve,though the scattering of the experimental data makes it quite difficult to evaluate the quality of the simulation results for the first 10s of tracer dispersion.In this case,the simulated mixing time is 39.7%higher than the experimental one.Overall,all these simulations provide a reasonable predic-tion of the qualitative time evolution of the tracer dispersion inside the vessel,but the mixing time is affected by an error ranging from 11%up to 40%,depending on the measuring position.It is also apparent that the whole homogenisation curve gives a more complete picture of the tracer dispersion process than the mere mixing time value.In order to investigate if different simulation choices would result in more accurate mixing time forecasts,the dependency of the predictions on a number of factors was2432G.Montante et al./Chemical Engineering Science 60(2005)2427–243701234204060t (s)c * (--)parison of dimensionless tracer concentration at the uppermost probe vs.time,for various time steps.T =29cm vessel stirred with four 6-bladed PBTs.Liquid:water,N =5s −1;fine line: t =0.01s;thick line: t =0.1s;dotted line t =1s.ly,the numerical settings,the value of the turbulent Schmidt number,the injection and probe po-sitions.It is worth noting that the position of the injection was found to affect the simulations (Bujalski et al.,2002a ),while the lack of knowledge of the precise probe location was mentioned as a possible error source in mixing time predictions (Patil et al.,2001).Grid independency of the results had already been verified for the liquid flow field predictions and,therefore,it was not considered at this step.The additional simulations performed for investigating the numerical issues have shown that the effect of the discretiza-tion scheme order is negligible:concentration profiles iden-tical to those obtained with a first-order scheme were ob-tained with a second order scheme.The adoption of a suit-able time step is,on the contrary,more important.In Fig.4the calculated curves corresponding to the highest probe are shown,as obtained with time steps of 0.01,0.1and 1s;only the first 60s are considered in order to magnify the differences among the curves.It can be seen that the time step selected originally (0.1s)was small enough:in fact,one order of magnitude smaller time step produces practically coincident profiles,while with t =1s the calculated con-centration curves are slightly different.This fact is not sur-prising,considering that a data acquisition rate of 4Hz was adopted for the experiments.(The results at lower heights are practically unaffected by the time step).It is also worth observing that the predicted value of t 95is nearly the same also when differences can be observed in the tracer concen-tration at the beginning of its dispersion inside the vessel.Therefore,the proper selection of the time step size seems to mainly affect the capability of catching the time evolution of the tracer dispersion rather than the global parameter.The influence of the tracer injection point was observed by changing its axial position with respect to the previous case,0.511.550100150200t (s)c * (--)σt =0.70σt =0.202468050100150200t (s)c * (--)σt =0.70σt =0.2(a)(b)parison of dimensionless tracer concentration vs.time obtained using different turbulent Schmidt numbers.T =29cm vessel stirred with four 6-bladed PBTs.Liquid:water,N =5s −1.Symbols:experimental data;thick line: t =0.70(same as in Fig.3);fine line t =0.20.(a)1st probe;(b)4th probe.i.e.,by moving it down by 1cm.No changes in concentration profiles have been observed at any of the four probes loca-tion.The same result was obtained once the tracer was added in a very small region close to the vessel top,correspond-ing to a single computational cell,rather than to a whole /4sector.The exact location of the probes was found not to be critical.In particular,several probe positions around the standard one were investigated:2cm above,2cm below,2cm closer to the vessel centreline and 1cm to the opposite direction (closer to the tank wall).In all the cases,only neg-ligible differences were observed in both the dimensionless concentration curves and the resulting mixing time (±0.7s around the “standard”position).The most critical parameter for tracer concentration pre-diction was found to be the value of the turbulent Schmidt number,which was varied from 0.7to 0.1.As can be ob-served in Fig.5(a)–(b),changing the t value from 0.7to 0.2G.Montante et al./Chemical Engineering Science 60(2005)2427–24372433-0.6-0.4-0.20.00.20.40.60.00.20.40.60.8 1.02r/TU /V t i pFig.6.Same comparison as in Fig.2.Simulations performed with standard settings and simplified simulation tools provided by a commercial code.Symbols:experimental data (Ranade and Joshi,1989);lines:simulations,each line corresponding to the profile calculated for each impeller.produces clearly different profiles.Those characterised by t =0.2interpret the experimental curves in a much better way than those with t =0.7:for the lower probes the exper-imental and the simulated curves are practically coincident.The corresponding errors in the predicted t 95values are 7.8%and 11.7%for the lower and upper probe,respectively,that are much lower than those mentioned above—predicted by using the standard value for the turbulent Schmidt number.It is worth noting at this point that homogenisation curves and mixing times have also been determined for several ro-tational speeds with one of the standard tools provided by some commercially available codes specifically developed for stirred vessels,through which this kind of simulations are claimed to be performed without the need of particu-lar expertise.Standard settings were used.As expected,grid size,which is one of the main items to select,has a severe influence on the predictions accuracy.However,in all the cases studied,the value of mixing time was severely over-estimated,by about 40%with a finer grid of about 400,000cells and by about 110%with a coarse grid of 73,000cells.These data are to be compared with an error in the order of 10%for all the conditions as obtained with the procedure outlined above.Scrutiny of the results shows that the main reason of such discrepancy between experimental and sim-ulated data may be found in severe inaccuracy in the flow field forecast when using automatically generated computa-tional grids and standardised solution procedures.This fact can be appreciated in Fig.6,where the same comparison as reported in Fig.2is shown for the radial profiles calculated with the finer of the two grids mentioned above:apparently,the agreement is much worse than that discussed previously and,in our opinion,this result is mainly due to the worse grid quality and,perhaps,to a poorer convergence criterion.4.2.T48simulationsThe mixing simulations relevant to the bigger vessel were performed considering two geometrically different sets of PBTs,various liquids characterised by Newtonian and non-Newtonian rheological behaviour,several impeller rotational speeds,both in the early and fully turbulent regimes.At-tention was devoted to the construction of the computa-tional grid.As in the case of other vessel,the reliability of the computed liquid flow field was evaluated in terms of power,N p ,and flow numbers,N Q ,that are reported along with the experimental data in Table 1.In all the cases,the power number computed from the torque on the baf-fles was practically equal to 1.5for both the 6-bladed and the 4-bladed PBTs (range of variation:1.47–1.50).Once again,the N p value computed from the turbulent dissipa-tion rate, ,was underestimated,as usually happens when the k – turbulence model is used:by about 50%for the 4-bladed PBTs and 25%for the 6-bladed PBTs.The sim-ulated flow number for each impeller,N Q ,was equal to about 0.85for the 4-bladed PBTs and 0.80for the 6-bladed PBTs.Both the computed dimensionless numbers are in good agreement with the experimental data.For calculat-ing the homogenisation curves,the same procedure as men-tioned above was adopted;the time step of 0.1s was cur-rently used (additional simulations with t =0.01s pro-duced results that were practically identical).In view of the results obtained with the T =29cm vessel,the dependency of the simulations on t was extensively tested in all these cases.Unless stated differently,the calculations were per-formed with the Fluent code.Only homogenisation curves obtained with the bottom probe were determined with this vessel.The vessel equipped with three 6-bladed PBTs is con-sidered first.Three different experimental conditions were considered:water stirred at 6s −1,a PVP solution of vis-cosity =23.7mPa s stirred at 7s −1,and a non-Newtonian liquid characterised by an average apparent viscosity of 14.5mPa s and stirred at 6s −1.As a picture of the tracer distribution with time,the contours of the tracer concen-tration for the last case at different agitation times in a vertical plane mid-way between two baffles is shown in Fig.7.As can be seen,the triple impeller configuration is characterised by gradual spreading and homogenisation of the tracer from the injection zone (top)to the other zones and in particular toward the vessel bottom.Limited radial gradients are noted especially on planes close to this position.The experimental and simulated tracer concentration curves obtained with t varying from 0.1to 0.7are shown for the three mentioned systems in Figs.8–10.The best agreement between experiments and simulations was ob-tained with the turbulent Schmidt number equal to 0.1.A very good agreement was found also for the aqueous so-lution of PVP and for that of Carbomer (Figs.9and 10),while for water a slightly worse match can be observed。