CFD modeling of turbulent reacting flow in a semi-batch stirred-tank reactor☆

Xiaoxia Duan ,Xin Feng ,*,Chao Yang ,*,Zaisha Mao

1 Key Laboratory of Green Process and Engineering,Institute of Process Engineering,Chinese Academy of Sciences,Beijing 100190,China

2 University of Chinese Academy of Sciences,Beijing 100049,China

1.Introduction

Semi-batch stirred reactors are often encountered in chemical and pharmaceutical industries;the particular reagent is added into the reactor with a feed-pipe.The scale of segregation of the feeding stream is reduced under the bulk inertial motion,followed with the stretch,deformation and engulfment of feed elements[1].At last,molecular diffusion reduces the intensity of segregation down to zero when the reagents achieve the molecular-scale homogeneity[2,3].For the mixing-sensitive multiple chemical reactions,if the characteristic reaction time is shorter than the time required for bringing homogeneity at molecular scale,micro-mixing can signi ficantly in fluence the product quality and selectivity[4].Mixing effects get more critical in reactor scale-up because the chemical reactions are scale-independent,while the local mixing time is sensitive to scale and position when the power per unit volume remains constant[2].Therefore,the studies of flow,mixing and chemical reaction that simultaneously occur in the stirred tank are significant to the design,optimization and scale-up of reactors.

In the CFD modeling of turbulent chemical reacting flow,the key problem is to select an appropriate micro-mixing model that can describe the local concentration information to close the chemical source term.Different mixing models have been developed and the most commonly used models include the engulfment model[5,6]and the probability density function methods,e.g. finite-node PDF[7,8],multipletime-scale turbulent mixer model[9–11],and DQMOM-IEM model[12].The DQMOM-IEM micro-mixing model is first presented by Fox[12]and has been successfully implemented to model the precipitation process[13–16]and parallel reaction scheme[17–19],combustion[20,21],etc.This method,compared with the transport-PDF method,is easily incorporated into existing CFD codes and does not introduce statistical errors[22].

Those successful applications validate that the DQMOM-IEM model can be used to simulate the reacting flow in the confined impinging jet reactor(CIJR)[14,15,17–19],plug- flow reactor[13],tubular reactor[23,24],etc.However,this model has not been used to model turbulent reacting flow in the stirred-tank reactors.The simulations of turbulent reacting flow are of great differences between semi-batch stirred tank and CIJR both in geometry and operating conditions,which are given as follows:(1)size of geometries:the diameters of stirred tank are often hundreds of times larger than the inlet diameters of CIJR that are a few millimeters[14,15,17–19].Furthermore,the interactions between the rotating impeller and stationary baffles produce complicated 3D- flow field,so the principal length scales and time scales of turbulent mixing process cover a wide range;(2)ratio between different feed streams:CIJR is carried out in a continuous operation,the flow ratesoftwo feed streams are often equal[17–19],butforthe semi-batch stirred-tank reactors,the feed volume is only 2%of the bulk flow which originally exists in the tank[5,10];Besides,the characteristics of feed ratio in fluence the initial concentrations if the quantities of all reagents are chemically equal[5,10];(3)feed time:the residence time of CIJR is far less than the feed times tfof semi-batch stirred tank.tfis often on the scale of several minutes which also makes the unsteady simulation time very long.All of these differences increase the difficulty ofsimulating turbulent reacting flow in stirred tanks.

The aim of this work is to employ the DQMOM-IEM model to predict the mixing effects on the course of parallel competing chemical reactions carried out in a semi-batch single-phase stirred-tank reactor.The simulation results are compared with the experimental data reported by Baldyga and Makowski[10].Some numerical strategies are adopted to solve the problems present in simulation of stirred tanks with complex geometry and differentoperations.In addition,the mixture fraction variance can be used to describe the degree of the local micromixing[25],so the state variable is used to mark the characteristics of reaction region in the semi-batch stirred-tank reactors.

2.Details of DQMOM-IEM Micro-mixing Model

DQMOM-IEM is a presumed probability density function model.This approach assumes that the joint composition PDF f?(ψ;x,t)can be approximately expressed as the sum of a finite number of multidimensional delta functions[26]:

where Neis the total number of environments,pn(x,t)is the probability or volume fraction of environment n,(x,t)is value of scalar α in environment n,Nsis the number of scalars,and ψ is the composition vector.

For the presumed PDF method,the non-reacting mixture fraction scalarand dimensionlessreaction progress variable are usefulto describe the mixing process between multiple non-premixed feed streams.

Considering the following second-order parallel competitive reaction system used in this work:

One stream contains reactant A with the initial concentration cA0,the other stream carries the pre-mixture B(cB0)and C(cC0)in the tank.

The mixture fraction is related to the local concentrations of reactants A,B and C,(i.e.,cA,cBand cC)de fined by[10]

The mixture fraction ranges from 0 to 1.ξ=1 is for the pure stream with concentration cA0.ξ=0 is for the pre-mixture inside the reactor.

The stoichiometric mixture fractions ξs1and ξs2are de fined respectively[12]

For the reaction system,the chemical species concentrations have to be expressed as the mixture fraction coupled with the reaction progress variable Y.The dimensionless reaction progress variables Y1and Y2for reactions(2)and(3)are de fined

According to the de finitions of Eq.(6),the reaction progress variables are normalized by the stoichiometric concentrations cA0ξs1and cA0ξs2.Furthermore,Y1and Y2are equal to 0 in the inlet and initial conditions.

Thus,chemical species concentrations can be expressed in terms of one mixture fraction variable ξ and two reaction progress variables(Y1and Y2).

As the micro-mixing testreactions,the firstreaction is instantaneous and the second reaction(side reaction)is a finite-rate reaction.When the feed stream of limited reactant A meets the other stream containing the premixed B and C,reactant A and B cannot coexist at the same spatial location.Thus,the first reaction progress variable Y1can be evaluated by the limiting value Y1∞as Eq.(8).

where 0≤Y2≤ξ/ξs2because the reaction progress variable is nonnegative.

If the two-environment DQMOM-IEM is employed to this reaction system, five transportation equations need to be solved.

1)The transport equation for the probability of environment 1(p1):

where the turbulent diffusivity ΓTis de fined as Eq.(10):

where Cμ=0.09,turbulent Schmidt number ScT=0.7,k is the turbulent kinetic energy and ε is the rate of energy dissipation.

Furthermore,the sum of probabilities is unity.Thus,the probability of the second environment can be calculated by p2=1-p1.

2)The transport equations for the weighted mixture fraction in environments 1(p1ξ1)and 2(p2ξ2):

where the firstand second termson the righthand ofthe transportequations are the micro-mixing terms and correction terms,respectively.The micro-mixing rate is de fined as:

where the micro-mixing rate constant C?is a function of the local Reynolds number,the detailed formulation can be referred to Liu and Fox[17].

3)The transport equations for the weighted reaction progress variable of the side chemical reaction in environments 1(p1Y21)and 2(p2Y22):

where the chemical source terms S2∞(ξ1,Y21)and S2∞(ξ2,Y22)are:

with

and the limiting value of the first chemical reaction progress variable in environments 1()and 2():

The reactant concentrations(cAn,cBn,cCn)and product concentrations(cRn,cSn)in the n th environments can be predicted by solving the transport equations of mixture fraction and reaction progress variables.Afterwards,the mean concentrations can be calculated using Eq.(21):

3.Numerical Simulation

3.1.Stirred tank and reaction system

As shown in Fig.1,the stirred tank used forsimulation in this work is the same as that used by Baldyga and Makowski[10].This vessel is equipped with a standard six-bladed Rushton turbine.

The reaction system for characterizing mixing efficiency is hydrolysis of ethyl chloroacetate in competition with neutralization of sodium hydroxide.The reaction rate constant k1=1.3 × 108m3·mol-1·s-1and k2=0.023 m3·mol-1·s-1at 20 °C.

Fig.1.Details of the stirred tank.

The NaOH solution with a concentration of 2000 mol·m-3is slowly fed into the stirred tank reactor which initially contains the premixture of HCl and CH2ClCOOC2H5,both with a concentration of 40 mol·m-3.The injection volume VA0is equal to 4.16 × 10-4m3,and the initial volume of stirred tank VBC,0is 0.02079 m3.Four different feed times are 8,10,15 and 20 min,respectively.The inner diameter of feedpipe is d=1 mm.The heightofthe feed point,located atthe middle plane between two adjacent baffles,is equal to the impeller clearance off the tank bottom and the radial position r=0.09 m.The segregation index XSis used to quantify the product distribution and calculated according to the volume-average concentrationafter the experiment[10].

whereandare equal and de fined as:

3.2.Simulation strategy

Numerical simulation is performed with the CFD commercial software ANSYS Fluent.The origin of three-dimensional coordinate is set at the center of the stirred tank bottom.The feed-pipe is created and the inlet is set as the velocity inlet boundary condition.The semibatch stirred tank is operated at a long feed time,so the velocity outlet is set at the bottom of the stirred tank to maintain the liquid level through draining out the equivalent amount working material.The computational domain which contains the whole tank geometry is divided into two sub-domains.The inner zone encompasses the rotating Rushton impeller and the multiple reference frames(MRF)approach is applied to model the motion of impeller.The remaining part of the tank is solved with the stationary reference frame equations.The boundariesbetween the rotating domain and the stationary domain are located at r=0.075 m and 0.06 m≤z≤0.14 m.The whole domain is discretized with tetrahedral and hexahedral grids using the GAMBIT mesh generation tool.The equisize skew is smaller than 0.75 to ensure good grid quality.The grids are re fined around the feed-pipe exit and the impeller where the variables have large gradient.

First,the flow field is solved with the Reynolds stress model(RSM)at steady-state.The RSM,which builds the model equations for the unclosed Reynolds stress and accounts for the anisotropy of turbulence,also belongs to RANS method.The previous simulation results[5–6]indicate that the RSM is relatively superior to the standard k–ε model on the prediction of turbulent flow field.Although LES approach better predicts the turbulent kinetic energy around the impeller than the RSM[27],the RSM method makes the computational cost more acceptable.For the sake of compromise between simulation accuracy and total computing time,the RSM approach is adopted to simulate the turbulent flow field in this work.The SIMPLEC algorithm is used to couple pressure and velocity.The second-order upwind scheme is adopted for the spatial discretization of transport equations.Then,the two-environment DQMOM-IEM micro-mixing model is loaded into the ANSYS Fluent through the user-de fined functions(UDFs).Based on the Fluent platform,the program written in the C programming language realizes the customization.The transport equations for five additional user-de fined scalars used in the DQMOM-IEM model are:

Those scalars are solved time-dependently with the second-order implicit time difference format.

The initialboundary conditions atthe velocity inletwhere reactant A is fed are:

and those for the computational domain inside the tank are:

Due to the long feed time,variable timestep approach is adopted to save the unsteady simulation time.The time step is set to 10-4s at the beginning of the simulation to overcome the numerical difficulties.

4.Results and Discussion

4.1.Veri fication of DQMOM-IEM procedure

Although the DQMOM-IEMmodelhas been used successfully for the simulation of turbulent reacting flow in the CIJR etc.,the procedure created in this work has to be veri fied before exploring the applicability of the DQMOM-IEM model in stirred tank reactor.The effect of mixing on the selectivity of parallel competitive reactions carried out in the CIJR is simulated with the two-environment DQMOM-IEMmodel.The detailed descriptions about geometry and operating parameters can be referred to Liu and Fox[17].The chemicalreaction sources added in the transport equationsofweighted reaction progressvariable and the computational boundaries are set according to the reaction system employed in CIJR[17].As shown in Figs.2 and 3,the predicted results about the mixture fraction variance 〈ξ'2〉and the segregation index X are in agreement with the experiment data[17]and the simulation results obtained with the same model by Liu and Fox[17],which confirm the validity of the DQMOM-IEM procedure.

4.2.Grid independence test

Fig.2.Profiles of〈ξ'2〉along the feedpipe axis predicted with differentmethods(Rej=400,dj=0.5 mm).

Fig.3.EffectofReynolds number Rej on the segregation index X(dj=0.5 mm,t r=4.8 ms).

Except the veri fication of the DQMOM-IEM procedure,it is also essential to analyze the sensitivity of the simulation results of turbulent flow field on the grid number.The turbulent flow field simulations are carried out with five different grids to achieve grid independence results.Grid 1 to grid 5 consist of 366381,601492,808938,1072953 and 1439368 cells,respectively.Fig.4 shows the axial profiles of dimensionless velocities components and turbulent kinetic energy at r=0.063 m.The results indicate that the turbulent kinetic energy is more difficult to achieve the grid independent in comparison with the velocities components.According to the numerical predictions,grid 4 is finally used for the simulations in this work.

4.3.Contour of volume fraction for environment

The ratio between reactor volume and feed volume is large for the semi-batch stirred-tank reactor.The feed stream consisting of reactant A slowly enters into the reactor from the feed inlet.Thus,the geometry and feeding characteristics of stirred tank determine the distribution of volume fraction of environment.

In Fig.5,contour plotofthe volume fraction ofenvironment1(p1)for the case N=214 r·min-1,tf=8 min is presented.At the beginning of the time-dependent simulation,the reactor is occupied by environment 2 except the velocity inlet.During the feeding period,environment 1 gradually enters into the reactor.The large gradient of p1decreases rapidly for a few seconds after the beginning of feeding and the system reaches the quasi-steady state due to the long feed time.As seen in Fig.5,environment 1 concentrates in the zone surrounding the feed-pipe exit and its volume fraction remains relatively small in the other region.

Fig.4.Comparison of simulation results or five different grids(N=214 r·min-1,r=0.063 m).(a)Dimensionless axial velocity.(b)Dimensionless radial velocity.(c)Dimensionless tangential velocity.(d)Dimensionless turbulent kinetic energy.

Fig.5.Contour plot of volume fraction of environment 1(N=214 r·min-1,t f=8 min).

4.4.Distributions of mixture fraction variance

By solving the transport equations of the probability of different environment(pn)and its weighted mixture fraction(pnξn),we can compute the mean mixture fraction 〈ξ〉and its variance 〈ξ′2〉according to Eqs.(29)and(30)[17]:

The mean mixture fraction 〈ξ〉is a non-reacting scalar which indicates the level of the macro-mixing and the mixture fraction variance〈ξ′2〉can be used as a measure for the small-scale segregation of the fluid[11].For the perfect mixing,the variance is zero.If the segregation is complete,it is equalto one.Fig.6 shows the segregation zone marked with different values of the variance ratio 〈ξ′2〉/〈ξ′2〉maxbetween 〈ξ′2〉and its maximum 〈ξ′2〉maxin the whole tank.With the increase of〈ξ′2〉/〈ξ′2〉max,the highlighted zone extremely shrinks and concentrates on the region close to the exit of feed-pipe.This result agrees with the conclusions by Vicum et al.[11]and Duan et al.[28].

4.5.Distributions of reactant concentrations

Initially,environment 1 contains pure reactant A.It is injected into the stirred-tank reactor from the beginning of time-dependent micromixing simulation.Environment 2 is originally filled with the mixture of reactants B and C.The interaction between environment 1 and environment 2 causes the compositional change,so the chemical reactions occur in both environments.Fig.7 shows the concentration distributions of reactants B and C in environments 1 and 2(at N=214 r·min-1,tf=8 min),respectively.The micro-mixing test reactions are fastand the feed position is in the discharge region ofthe Rushton impeller.Thus,reactants B and C are consumed almost in the vicinity of the feed-pipe exit.As shown in Fig.7,the consumption zone of reactant C is more concentrated than that of reactant B on the region where the value of the mixture fraction variance is high.This also indicates that the side reaction mainly occurs in the large segregation zone[28].Furthermore,the reaction zone extends to the tank walldue to the strong radial discharge flow.The characteristics of concentration distributions re flect the flow pattern generated by the standard Rushton turbine.

4.6.Effect of agitation speed and feed time on segregation index

The effect of agitation speed on the reaction process for two cases(feed time tf=20 min and tf=15 min,respectively)are reported in Figs.8 and 9 in terms of the segregation index XSwhich represents the selectivity of by-product.For the parallel competitive reactions given by Eqs.(22)and(23),the acid–base neutralization reaction occurs only and the yield of the undesired product tends to zero in the well mixed condition,because there is a great disparity between the reaction rate constants k1and k2.However,in the complete segregation,the initial reactant concentrations determine the distribution of products.According to experimental parameters listed in Section 3.1,the segregation index is 0.5 that is also the upper limit.

The segregation indexes predicted by the DQMOM-IEM micromixing model are compared with the experimental data and the numerical simulation results in literature[10].Model I is the multipletime-scale turbulent mixing model and Model III neglects the concentration fluctuations.As shown in Figs.8 and 9,the predicted segregation index decreases with the increase of agitation speed because the micro-mixing performance is intensi fied by the enhanced turbulence level.

There is apparent difference between the predicted results with Model III and micro-mixing model(either DQMOM-IEM model or multiple-time-scale turbulentmixing model).This indicates that the effect of micro-mixing on the course of parallel competing chemical reactions is significant.The suitable closure models have to be used to calculate the non-linear chemical source term and describe the subgrid mixing.

In addition,the predicted results with Model I are really in good agreementwith the experimentdata as well.Itis a presumed PDF method and the beta PDF is employed to represent the mixture fraction PDF.For the non-premixed turbulent flows with multiple feed streams,the mixture fraction PDF is poorly approximated by a beta PDF[12].In this work,the semi-batch stirred tank is predicted and there is only one feed stream.For this binary mixing,the mixture fraction PDF is wellapproximated by a beta PDF[12].In orderto selectthe mostappropriate modelto simulate differentcases,the comprehensive comparison aboutthe predicted ability between the DQMOM-IEMmodeland multitime-scale turbulent mixing model should be conducted in more CFD studies.

Altering the feed time at N=214 r·min-1,the predicted results are plotted in Fig.10 and have the same trend with experimentdata.For the fixed volume of reactant A,the feed velocity is inversely proportional to the feed time.Atthe relatively low feed velocity,the feed stream is more easily dispersed into the mainstream which is bene ficial to maximize the yield of the desired product.

The visualization of reaction plume at the end of the feeding is also presented with the mean concentration of reactant A in Fig.11.As it can be seen,the volume of zone marked with 〈cA〉＞0.0001 mol·m-3at the feed time tf=20 min is also smaller than it at tf=8 min and the reaction zone is more close to the feed point with the increase of feed time.These visualized results also prove that the injected reactant A can be more easily dispersed and consumed by the reactants B and C at low feed velocity.

Fig.7.Contour plots of reactant concentrations in different environments(N=214 r·min-1,t f=8 min,in the y=0 plane).(a)Concentration of reactant B in environment 1.(b)Concentration of reactant C in environment 1.(c)Concentration of reactant B in environment 2.(d)Concentration of reactant C in environment 2.

Fig.8.Effect of agitation speed on X S at t f=20 min.

5.Conclusions

Fig.9.Effect of agitation speed on X S at t f=15 min.

Fig.10.Effect of feed time on X S at N=214 r·min-1.

In this work,the two-environment DQMOM-IEM model coupled with CFD is implemented to predict and quantify the effects of mixing on the distributions of products in a semi-batch stirredtank reactor.The segregation index of parallel competitive chemical reactions decreases with the increase of feed time.The higher the agitation speed is,the lower the yield of by-product.In the semibatch stirred tank,the reaction zone mainly concentrates on the region close to the feed-pipe exit where the mixture fraction variance is large.The concentration distributions of reagents and the visualization of reaction plume can be presented with the twoenvironment DQMOM-IEM micro-mixing model,which cannot be achieved by the empirical and theoretical models solved in the Lagrangian framework.Although the simulation results agree with the experimental data,it is essential to investigate the predictive behavior of DQMOM-IEM micro-mixing model used for stirred tanks in other geometry parameters and operating conditions in the further works.

Nomenclature

C impeller off-bottom clearance,m

ciconcentration of specie i,mol·m-3

cinconcentration of specie i in environment n,mol·m-3

〈ci〉 mean concentration of specie i,mol·m-3

D impeller diameter,m

djdiameter of the impinging jets,m

H height of stirred tank,m

k turbulent kinetic energy,m2·s-2

N agitation speed,r·min-1

p probability

RejReynolds number in the CIJR

r radial distance from Z-axis,m

T stirred tank diameter,m

t time,s

tffeed time,s

trcharacteristic reaction time,ms

X,XSsegregation index

x spatial position vector

Y reaction progress

z axial distance from tank bottom,m

γ micro-mixing rate,s-1

ε rate of energy dissipation,m2·s-3

ν kinematic viscosity,m2·s-1

ξ mixture fraction

ξ1mixture fraction in environment 1

ξ2mixture fraction in environment 2

Fig.11.Visualization of reaction plume with 〈c A〉＞0.0001 mol·m-3 at the end of the feeding.(a)N=214 r·min-1,t f=8 min.(b)N=214 r·min-1,t f=20 min.

Subscripts

A,B,C,R,S chemical species

j direction

0 initial value

References

[1]J.Baldyga,J.R.Bourne,A fluid mechanical approach to turbulent mixing and chemical reaction part II micromixing in the light of turbulence theory,Chem.Eng.Commun.28(4–6)(1984)243–258.

[2]E.L.Paul,V.A.Atiemo-Obeng,S.M.Kresta,Handbook of Industrial Mixing:Science and Practice,John Wiley&Sons,USA,2004.

[3]P.V.Danckwerts,The effect of incomplete mixing on homogeneous reactions,Chem.Eng.Sci.8(1–2)(1958)93–102.

[4]J.C.Cheng,X.Feng,D.Cheng,C.Yang,Retrospect and perspective of micro-mixing studies in stirred tanks,Chin.J.Chem.Eng.20(1)(2012)178–190.

[5]O.Akiti,P.M.Armenante,Experimentally-validated micromixing-based CFD model for fed-batch stirred-tank reactors,AIChE J.50(3)(2004)566–577.

[6]Y.Han,J.J.Wang,X.P.Gu,L.F.Feng,Numerical simulation on micromixing of viscous lf uids in a stirred-tank reactor,Chem.Eng.Sci.74(0)(2012)9–17.

[7]D.L.Marchisio,A.A.Barresi,CFD simulation of mixing and reaction:the relevance of the micro-mixing model,Chem.Eng.Sci.58(16)(2003)3579–3587.

[8]Z.Wang,Q.H.Zhang,C.Yang,Z.S.Mao,X.Q.Shen,Simulation of barium sulfate precipitation using CFD and FM-PDF modeling in a continuous stirred tank,Chem.Eng.Technol.30(12)(2007)1642–1649.

[9]J.Baldyga,M.Henczka,L.Makowski,Effect of mixing on parallel chemical reactions in a continuous- flow stirred-tank reactor,Chem.Eng.Res.Des.79(8)(2001)895–900.

[10]J.Baldyga,L.Makowski,CFD modeling of mixing effects on the course of parallel chemical reactions carried out in a stirred tank,Chem.Eng.Technol.27(3)(2004)225–231.

[11]L.Vicum,S.Ottiger,M.Mazzotti,L.Makowski,J.Baldyga,Multi-scale modeling of a reactive mixing process in a semibatch stirred tank,Chem.Eng.Sci.59(8–9)(2004)1767–1781.

[12]R.O.Fox,Computational Models for Turbulent Reacting Flows,Cambridge University Press,Cambridge,2003.

[13]L.Wang,R.O.Fox,Comparison of micromixing models for CFD simulation of nanoparticle formation,AIChE J.50(9)(2004)2217–2232.

[14]E.Gavi,L.Rivautella,D.L.Marchisio,M.Vanni,A.A.Barresi,G.Baldi,CFD modelling of nano-particle precipitation in confined impinging jet reactors,Chem.Eng.Res.Des.85(5)(2007)735–744.

[15]E.Gavi,D.L.Marchisio,A.A.Barresi,On the importance of mixing for the production of nanoparticles,J.Dispers.Sci.Technol.29(4)(2008)548–554.

[16]J.Sierra-Pallares,D.L.Marchisio,M.T.Parra-Santos,J.García-Serna,F.Castro,M.J.Cocero,A computational fluid dynamics study of supercritical antisolvent precipitation:mixing effects on particle size,AIChE J.58(2)(2012)385–398.

[17]Y.Liu,R.O.Fox,CFD predictions for chemical processing in a confined impinging-jets reactor,AIChE J.52(2)(2006)731–744.

[18]E.Gavi,D.L.Marchisio,A.A.Barresi,CFD modelling and scale-up of confined impinging jet reactors,Chem.Eng.Sci.62(8)(2007)2228–2241.

[19]D.L.Marchisio,Large eddy simulation of mixing and reaction in a confined impinging jets reactor,Comput.Chem.Eng.33(2)(2009)408–420.

[20]V.Raman,H.Pitsch,R.O.Fox,Eulerian transported probability density function subfilter model for large-eddy simulations of turbulent combustion,Combust.Theor.Model.10(3)(2006)439–548.

[21]R.Yadav,A.Kushari,V.Eswaran,A.K.Verma,A numerical investigation of the Eulerian PDF transport approach for modeling of turbulent non-premixed pilot stabilized flames,Combust.Flame 160(3)(2013)618–634.

[22]J.Lee,Y.Kim,Transported PDF approach and direct-quadrature method of moment for modeling turbulent piloted jet flames,J.Mech.Sci.Technol.25(1)(2011)3259–3265.

[23]J.Akroyd,A.J.Smith,L.R.McGlashan,M.Kraft,Comparison of the stochastic fields method and DQMOM-IEM as turbulent reaction closures,Chem.Eng.Sci.65(20)(2010)5429–5441.

[24]J.Akroyd,A.J.Smith,L.R.McGlashan,M.Kraft,Numerical investigation of DQMOMIEM as a turbulent reaction closure,Chem.Eng.Sci.65(6)(2010)1915–1924.

[25]R.O.Fox,CFD models for analysis and design of chemical reactors,Adv.Chem.Eng.31(2006)231–305.

[26]R.O.Fox,On the relationship between Lagrangian micro-mixing models and computational fluid dynamics,Chem.Eng.Process.37(6)(1998)521–535.

[27]B.N.Murthy,J.B.Joshi,Assessment of standard k–ε,RSM and LES turbulence models in a baffled stirred vessel agitated by various impeller designs,Chem.Eng.Sci.63(22)(2008)5468–5495.

[28]X.X.Duan,X.Feng,C.Yang,Z.S.Mao,Numerical simulation of micro-mixing in stirred reactors using the engulfment model coupled with CFD,Chem.Eng.Sci.140(2)(2016)179–188.