Numericalsimulation ofstirred tanks using a hybrid immersed-boundary method☆

Shengbin Di,Ji Xu ,QiChang ,3,WeiGe ,*

1 State Key Laboratory ofMultiphase Complex Systems,Institute ofProcess Engineering,Chinese Academy ofSciences,Beijing 100190,China

2 University of Chinese Academy ofSciences,Beijing 100049,China

3 Schoolof ChemicalEngineering and Technology,Tianjin University,Tianjin 300072,China

1.Introduction

Stirred tanks are widely used in chemicaland process industries,and hence deserve comprehensive experimentaland computationalinvestigation.Experimentalmeasurement techniques such as laser-doppler velocimetry(LDV)and particle image velocimetry(PIV),have signi ficantly contributed to our understanding of the flow structure in stirred tanks.However,most experimentalstudies have limitations,such as time and cost,scale-up dif ficulties,and the dif ficulty ofdirectly measuring some quantities like localvelocity especially for high-viscosity fluids[1].Computational fluid dynamics provides an alternative tool for unveiling the details of fluid flow and heat transfer in stirred tanks,but faces severalchallenges.The main dif ficulties lie in the relative motion between the impeller and vesseland the complex geometry of the impeller[2].To address these dif ficulties,severalnumericalmethods have been developed.

Conventionally,boundary-conforming nonorthogonalgridsare used to handle the complex geometries[3].In the so-called body-fitted approaches[4],the grid lines follow the boundaries exactly and the boundary conditions can therefore be implemented directly.However,the generation of such grids is both dif ficult and time consuming,and it has to be repeated or modi fied during the simulation for moving boundaries.Furthermore,the discretized governing equations for body-fitted grids are much more complex than those for Cartesian grids,which increase the computationalcost and decrease the numericalstability considerably.As a result,such a method is seldom directly applied to the simulation of stirred tanks.A non-inertial reference frame fixed to the moving body can be used to avoid mesh regeneration for moving boundaries,although the governing equations have to take a more complicated form.Multiple reference frame(MRF)[5]and sliding mesh(SM)models[6]are typicalexamples of this approach.They are among the most commonly used methods for the simulation ofstirred tanks and available in mainstream commercial codes[7–10].Both methods solve the governing equations in two domains:an inner domain containing the rotating impeller and an outer domain containing the stationary tank walland baf fles.In SM,the meshes are allowed to shear and slide to accommodate the relative motion at the domain interface,while in MRF,only a stationary solution of the flow field is provided for a given position of impeller with respect to the baf fles[11].The main disadvantages of SMare the explicit acceleration terms added to the momentum equation and the related numericalproblem atthe domain interface.For complex impellerdesigns and intermeshing extruders,both MRF and SM might be infeasible or computationally very demanding[12].

Immersed boundary(IB)methods,firstproposed by Peskin[13],are another category ofmethods for complex geometries thatcan be based on simple and regular Cartesian grids,and have received much attention recently[14–16].IB methods have higher computationalef ficiency when compared with body-fitted methods[16–19].In IB methods,the Cartesian grids do notnecessarily conformto the boundaries.The effects of the boundaries on the fluid are represented by additionalbody forces in the momentum equations.Although the number of grid points increases more quickly than in body-fitted methods as the Reynolds number increases,the per-grid-point-operation count is notably reduced owing to the absence of grid transformation[20].In recent years,IB methods have been applied in the simulation ofcomplex fluid flow in stirred tanks because of their ef ficiency and capacity to deal with moving boundaries directly[2,21,22].Unlike SMmodel,there is no constraint on the prescribed motion of the rigid bodies in IB methods[23].

Among various IB methods[15,16,24–26],two typicalapproaches are the bases ofthis study.In the firstapproach,the boundaries are represented by a set of Lagrangian points which are handled with direct forcing methods such as that proposed by Uhlmann[16].No-slip boundary conditions are assigned to these Lagrangian markers.The positions of the Lagrangian markers are updated in each time step if the boundaries are moving.As only the surfaces of the immersed bodies should be represented by the markers,the number of markers is relatively smallcompared with the totalnumber of Cartesian grid points.Each Lagrangian point is connected with a number of Cartesian grids around it.Mapping of the variables between Lagrangian and Eulerian locations is performed employing a regularized delta function.This technique effectively reduces non-physical force oscillations existing in moving-boundary simulations[27],however,it is not necessary for stationary-boundary problems.Furthermore,when the IBs lie close to the borders of the computational domain,where the separation is smaller than the support width of the delta function,it is dif ficult to find enough Cartesian grid points for velocity interpolation and force distribution when using this technique.

In the second approach,the body force is evaluated directly atthe IB using the desired velocity of the solid part.However,the fluid and solid velocities are smoothly interpolated using a solid volume fraction atthe boundary grid points.When the solid part is stationary,the solid fraction at the computational grid points does not change during the computation.As a result,the evaluation of the solid fraction needs to be carried out only once.The method of Kajishima et al.[28]belongs to this category.

Itis noted thatthe two approaches described above are more suitable for moving and stationary boundaries,respectively(for simplicity,they are referred as the moving and stationary approaches hereafter).Both moving and stationary boundaries,such as the boundaries ofrotating impellers and stationary tank walls,are present in a stirred tank simultaneously.However,most simulations use only one approach for all boundaries[22,29,30].Tyagi etal.[30]proposed a method to combine the advantages of the IB method to represent moving rigid geometries and the multi-block structured curvilinearmesh to representcomplex domains.However,it is stilla challenge to generate body-fitted mesh for complex geometries such as the helicalcoils inside a tank.To this regard,we willpropose in this work a hybrid IB technique in which the complex geometries are dealt with using different IB methods on a Cartesian grid.

On the other hand,one practicaldif ficulty in simulating complex flows in stirred tanks is the high computation cost.Consequently,we implement the IB method on a CPU–GPU hybrid architecture for massively parallel computing.Algorithms having good parallelism are therefore desirable.Since the fluid and IBs are discretized on two separate grids,their parallelization should be handled differently.The parallelization for Eulerian grids is easy because the grid points are fixed in time and are assigned to the same subdomain for the whole computation period[31].However,the Lagrangian markers move throughout the computationaldomain,so a given IB point may move across the subdomain border during one time step.Furthermore,the interactions between the fluid and IBs are another source of dif ficulty in the design of a distributed algorithm[32].To cope with these problems,we developed a partial velocity interpolation and force distribution scheme,which is described in Section 2.

To validate the proposed method,we firsttest its accuracy and robustness in two numericalexperiments,incorporating large eddy simulation(LES)which is found to be effective for turbulent flow in stirred tanks[33,34].Then,as an application,the discrete element method(DEM)is coupled with the presentmethod to carry outdirectnumerical simulation(DNS)of fluid–solid systems in stirred tanks.Satisfactory results reveal that complex flows in stirred tanks can be accomplished using hybrid IB methods on supercomputers.Itwillbecome a powerful toolfor optimaldesign of industrialstirred tanks if the method can be further developed.Finally,main conclusions are drawn with perspectives on wider applications of the present method.

2.Numerical Method

2.1.Governing equations

In the currentwork,we focus on incompressible fluid with constant properties.The governing equationsofthismethod describe the conservation ofmomentum and mass of the flow.They are established for the fluid phase and the solid objects(or boundaries)separately,and coupled through the force terms describing the interactions between them.

2.1.1.Fluid-phase equations

The governing equations for unsteady incompressible flow in the IB methods take the generalform[16,17]

whereρfandμfare the fluid density and dynamic viscosity,respectively,while u,and p are the velocity and pressure,respectively.They are basically the Navier–Stokes equations with additional body force terms f determined by the no-slip boundary conditions at the fluid–solid interface.f is non-zero only in the vicinity of fluid–solid interfaces.

2.1.2.Solid-phase equations

Only rigid bodies are considered in the present study.Given a prescribed motion,the desired velocity field within the solid region is[16]

where ut,ωand x0are the translationaland rotationalvelocity and center coordinates ofa rigid body.Otherwise,ifthe motion of the rigid solid domain is due to the fluid force,the motion is controlled by Newton's equations for linear and angular momenta ofa rigid body,

where vsandωsare the linear velocity at the mass center and the angular velocity of the rigid body,respectively.Fsand Tsrepresent the hydrodynamic force and torque acting on the rigid body due to the fluid.Gsand Nsare the externalforce and torque acting on the rigid body.msis the mass of the rigid body and Isis the moment ofinertia of the rigid body.

2.2.Large eddy simulation

In LES,the large-scale unsteady turbulent motions are resolved,while the effects of the smaller scale motions are modeled[35].In the previous studies[22,33],the well-known Smagorinsky model[36]was usually used:

whereand CSis the Smagorinsky constant.However,near the walls or the quiescent flow regions,the subgrid-scale stresses should vanish,which is notautomatically guaranteed in the standard Smagorinsky model[33].In this study,we adopted a Lagrangian dynamic subgrid-scale model[37]in which CSis de fined as:

where FLMand FMMcan be calculated by[37]

whererepresents filtering at scale 2Δand

where H(x)is the ramp function(H{x}=x if x≥0,and zero otherwise).

The turbulent kinetic energy k is for both resolved and unresolved scales,but Hartmann etal.[38]pointed out that the subgrid-scale kinetic energy in a midway baf fle plane is about 0.1%of that at the grid-scale near the impeller and less in the rest domains.As a result,only the resolved kinetic energy is considered in the present study

2.3.Determination of the source term f

According to the direct forcing method,the boundary body force f in Eq.(1)can be simply calculated as[26]

whereΔt is the time step,udis the desired velocity at an arbitrary point in the solid object,unis the velocity computed at step n,andΩis the domain in the vicinity of the IB.In general,the surface of the solid region where un+1=uddoes not coincide with the Cartesian grid points,especially when a staggered grid is applied,an interpolation procedure is thus needed.

As mentioned in the introduction,a hybrid technique for determining the source term f is proposed in this paper.That is,the scheme of Uhlmann[16]is used for moving boundaries and that of Kajishima et al.[28]is used for stationary boundaries.The details are described below.

2.4.IB method for moving boundaries(moving approach)

Uhlmann[16]proposed a direct forcing method where the IBs are represented by a set of Lagrangian markers.This method allows a smooth transfer between the Eulerian and Lagrangian representations while at the same time avoiding strong restrictions of the time step[16].The force term is evaluated at the Lagrangian markers to suppress the oscillations of the hydrodynamic force due to insuf ficientsmoothing when dealing with moving boundaries;i.e.,

where the uppercase letters denote the variables at the Lagrangian pointsand Udis the desired velocity atthe location on the fluid–solid interface and can be obtained from Eq.(3).The intermediate velocity~U can be evaluated by interpolation from its Eulerian counterpart~u,

where h is the spacing of the Eulerian grid points and~u can be computed from Eq.(1)without the body force f.δhis the discrete delta function usually constructed in the form[39]

where x,y,z and X,Y,Z are the components of x and X,respectively,andφis the one-dimensionaldiscrete delta function:

To obtain the forcing term in Eq.(1),the body force acting on Lagrangian markers should be distributed to the Cartesian grid points:

where Nmis the totalnumberofmarkers andΔVlis the volume assigned to the l th Lagrangian marker;the de finition refers to Uhlmann[16].

2.5.IB method for stationary boundaries(stationary approach)

Kajishima et al.[28]proposed an extremely simple scheme for modeling fluid–solid interactions.At the interface,the fluid and solid velocities are smoothly connected using the solid volume fraction at each Cartesian grid point as a weight factor:

whereαrepresents the solid volume fraction at the computationalgrid point,and~u and udare the velocities within the fluid and solid domains,respectively.For the stationary tank wall,ud=0,and the evaluation of the desired velocity at a Cartesian grid point can then be de fined as

Similarto the case for the directforcing method based on the no-slip boundary condition,the body force term at the boundary grid point should be

As the treated boundary is stationary,the volume fraction at the Cartesian grid points will not change during the computation.This means the volume fraction needs to be evaluated only once,usually before the simulation by a separate code.

Obviously,accurate determination of the solid volume fractionαis the key to this method.The discrete approximation approach[40]is employed in this study.The grid points lying suf ficiently close to the interior of the IB haveα=1 directly.Then,for grid points intersecting with the IBs,a series ofmarker points are generated uniformly in each of the three directions,with an intervalmuch smaller than the interval of Cartesian grid points h(e.g.,h/M,where M is an integer),and hence,the totalnumber ofmarker points within a grid cellis M3.A ray tracing technique[41]is then employed to countthe numberofmarkers within the solid region.Ifthis number is n,we immediately haveα=n/M3.Larger M gives higher accuracy but requires more computationaltime.However,as the computation is only carried outonce,large M is affordable in obtaining a fairly accurate estimation ofα.

Finally,in the context of the hybrid IB technique in this work,the body force at a Cartesian grid point can be represented as a sum of the two components for moving and stationary boundaries:

Naturally,the two terms are nonzero only at different times and positions except for some short and localized contact.

2.6.Parallelimplementation

The hybrid IB technique described above is implemented and parallelized on a supercomputer with both multi-core CPUs and manycore GPUs.The Navier–Strokes equations for the fluid flow are solved on a fixed Cartesian grid,which can be parallelized straightforwardly.In the parallelimplementation,the regularrectangle computationaldomain is divided into a number ofsub-domains,and each process(running on a CPU core and the corresponding GPU)is assigned to a sub-domain.The data exchange between sub-domains is communicated through interprocess data transfer in the industry-standard MPIprotocol[42].

However,the IBs crossing subdomain borders may pose challenges,which willbe tackled in this section.As the solid volume fraction ofeach grid pointis stored locally,no additionaldif ficulty is associated with the parallelization of the stationary approach.Meanwhile,for the moving approach,the moving Lagrangian markers at the borders between subdomains cause a variety of complications.The interpolation of the velocity on markers(Eq.(13))requires the exchange ofextra flow information at the subdomain borders.In general,to minimize the communication,tagging and searching of Lagrangian points and construction of the delta functions should be done ef ficiently[23].

The two steps for processing the IB markers,force distribution and velocity interpolation,both loop over the same indexes determined by the support of the regularized delta function,which has a width of four grid points in this study.Taking the two-dimensional case in Fig.1(a)for example,we see thatthe IB maker m interacts with the surrounding 4×4 Eulerian grid points(point 1 to point 16).Considering the staggered grid used,in the extreme case shown in Fig.1(a)when the Lagrangian markers are very close to the subdomain border,Eulerian grid points 1–8 are ghost grid points whose values should be transferred from the neighboring process.Therefore,in general,two layers of ghost grid points are needed to perform the interpolation and distribution steps while only one layer is needed ifa second order discretization scheme is applied.Therefore,two additionalcommunications are needed during the IB stage in this extreme case.

?Interpolation step:the data of the ghost Cartesian grid points(grid points 1–8)should be transferred from the adjacent process to compute the velocities of the Lagrangian markers(Eq.(13)).One extra layer ofghost grid points is needed when compared with the case of the fluid solver.

?Distribution step:the body force source of Lagrangian markers should be spread to the surrounding connecting Cartesian grid points(grid points 1–16).The velocity of the Cartesian grid points should be updated by adding the source terms to them.As a result,the calculated body force source term of Cartesian grid points 1–8 should be further transferred back to the neighboring process and added to the corresponding grid points.

Obviously,these two additionalcommunications greatly complicate the parallelization and increase the communication tremendously.

To simplify the above procedure,a partialvelocity interpolation and force distribution method is developed in this study,where the Lagrangian markers interact with only some of the 4×4 Eulerian grid points at the subdomain border.In this strategy,only one layerofghostgrid points is required and the interpolation and distribution steps are modi fied accordingly.It should be noted that only the Lagrangian markers assigned to the present process are computed.In Fig.1(b),for example,only Eulerian grid points 5–16 are used for the velocity interpolation and the body force of m isonly distributed to grid points9–16.Asa result,overlapping ofone ghostgrid pointis suf ficientand the body force does notneed to be transferred to the neighboring process any more.Two factors should be incorporated into the interpolation and distribution steps to ensure that the velocities are interpolated exactly and the body forces exerted by the Lagrangian markers are fully transferred to the Cartesian grid points[43,44]:

Fig.1.Schematic representation of the interpolation and distribution operations at the subdomain borders:(a)two layers ofghost Cartesian grid points,(b)one layer ofghostgrid points.

where fvand Flare de fined by

HereΩvis the Cartesian domain for the interpolation of the intermediate velocity atthe Lagrangian markers(in the extreme case ofFig.1(b),it includes grid points 5–16),and Ωfis the domain for the body force spreading(in the extreme case of Fig.1(b),it includes grid points 9–16).Naturally,when a Lagrangian marker lies more than two grid point widths from the subdomain border,ΩvandΩfwillconsist ofall 16 surrounding Eulerian grid points,and fv=1 and Fl=1.Applying this strategy,the communication becomes more ef ficient.This strategy is derived from and an improvement of the half-distributed forcing strategy for IB methods,which was demonstrated to be feasible in our previous publication[44].Furthermore,the forces and torques added to the fluid through the distribution procedures(Eq.(16))are not changed[44].From the discussions above,we conclude thatthe proposed partialinterpolation and distribution strategy willnotlead to additionalcommunications.As a result,the parallelization is much simpler and more ef ficient.

3.Validation

To evaluate the ef ficiency of the hybrid IB technique and the accuracy of the proposed parallelization strategy,we simulated the flow field in two stirred tanks with this method.The main results are reported in the follows.

3.1.A simple stirred tank

The method is firstapplied to a cylindricalunbaf fled stirred tank with eightstraightblades atthe mid-heightas an impeller.This tank has been studied experimentally[45]and numerically[2].Aschematic illustration of the stirred tank is presented in Fig.2,where the height of the tank H=10 cm.The impeller rotates at n=100 r·min-1,which leads to an impeller Reynolds number Re=ND2/ν=1024,in terms of the turbine diameter D,turbine frequency N and fluid kinematic viscosityν.Direct numericalsimulation is possible for the present simple stirred tank in which geometricaldimensions and rotation speed are both small.The rectangular computational domain is discretized using a uniform Cartesian grid of 256×256×256 grid points,which is finer than that used in the DNS by Verzicco etal.[2].The computationaldomain is divided into 2×2×2 subdomains in the parallelcomputing.No-slip boundary condition is imposed at the bottom of the tank,while free-slip condition is assigned at the top.Within the computationaldomain,noslip boundary conditions are enforced for the cylindricaltank walland the impeller using the hybrid IB technique described above.

Fig.2.Schematic illustration of the stirred tank(adapted from Dong etal.[45]and Verzicco et al.[2]).

We compare the computationalef ficiencies of the two IB methods described above in handling the stationary and moving IBs,as listed in Table 1.Note that,ifa moving boundary is handled by applying the stationary approach(Case 1),the grid volume fraction should be evaluated in each time step because the position of the boundary changes with time.In a stirred tank,the surface area of the rotating impeller is generally smallcompared with the stationary wall,and the totalnumber of Lagrangian markers is limited.Conversely,ifthe tank wallis also simulated with the moving approach(Case 2),the number of Lagrangian markers willincrease by approximately an order of magnitude(from 15611 to 132183 for the present stirred tank),and the computational time willincrease sharply.From Table 1,we see that the time required for the hybrid IB technique is much less than the traditionalapproach to a single IB for both stationary and moving boundaries.

Table 1 Comparison of the ef ficiency of the two IBmethods in the case ofa stirred tank.SAand MA denote the stationary approach and moving approach,respectively.Data are the time taken for 100 time steps.In Cases 1 and 3,a smallvalue of M(=3)is used when evaluating the grid solid fraction in the stationary approach

Furthermore,the averaged velocity components are compared with experimental data[45],as shown in Fig.3.Data are evaluated at two different axiallocations z/H=0.5 and 0.7.A satisfactory agreementis obtained for allthree velocity components in these three cases,except for smalldeviations in the axialvelocity at z/H=0.7 and in the radial velocity at z/H=0.5 in Case 1.This may be ascribed to the approximation of the solid volume fraction due to the small M used.

3.2.Turbulent flow in a Rushton stirred tank

To verify our method in more practicalapplications,we then simulated a tank stirred by a Rushton turbine,and LES is incorporated to describe the turbulent flow in it[21,22,33,46].The problem has been studied by Hartmann et al.[38]who presented the comparison of LES and RANS simulation with LDA measurements.

Fig.3.Comparison of the computed(lines)average velocity components with the experimentaldata from Dong etal.[45](symbols▲)along the radialdirection atdifferentaxiallocations z/H=0.5(a,b,c)and 0.7(d,e,f).u a,u r and u t represent axial,radialand tangentialvelocity components,respectively.

Fig.4.Cross-section of the tank driven by a Rushton turbine[33].

The stirred tank is encapsulated by a standard cylindrical vessel 150 mm in diameter(T=150 mm)with four equispaced baf fles.The geometry is sketched in Fig.4.The density of the working fluid(ρ)is 1039 kg·m-3and the dynamic viscosity(μ)is 15.9 mPa·s.The rotating speed is 2672 r·min-1,resulting in an impeller Reynolds number Re=7300.This Reynolds number is not within the fully turbulent regime,however,the flow is fully turbulent near the turbine[29].A uniform Cartesian grid of 5123is used in the computation and divided into 43subdomains.Starting from the quiescent condition,the stirred tank flow requires approximately 100 turbine revolutions to reach a statically steady state[29].The start-up was firstcomputed on the 2563grid for about100 revolutions and then interpolated onto the 5123grid for another 15 revolutions.Then,the flow variables were collected from the last 10 revolutions.In the computation,the timestep was selected according to the tip speed Utipand the grid width Δh,i.e.,ΔtUtip/Δh=0.2.The computation requires approximately 1122 timesteps for a full revolution on the grid of2563and 2244 timesteps on the grid of 5123.

To compare the present LES results with experimentaland numericaldata from other publications,the phase averaged variables are evaluated.Fig.5 shows the experimental,RANS and LES results of the phase averaged flow field and their associated kinetic energy.Fig.5(d)is the present LES results,Fig.5(c)is the LES results using a standard Smagorinsky subgrid-scale model with constant CS=0.1[33].The main features of the flow field are the two large circulation loops on both the top and bottom of the impeller.The working fluid is pumped radially towards the tank wall,partly goes upward above the impeller and partly goes downward below the impeller.Finally,the fluid flows into the impeller region along the axis.In this figure,the RANS and LES for both subgrid-scale models show good agreementto the experiment.The magnitudes of the kinetic energy of the simulation results are comparable with that of the experiment.As can be seen,the kinetic energy in the flow jet area is much higher than the rest of the tank.

A quantitative comparison of phase-averaged radialand tangential velocity components atthree radiallocations is shown in Fig.6.A satisfactory agreement is obtained between the simulation and the experiment except a smalldeviation for radialvelocity at r/T=0.317 and the tangential velocity at location r/T=0.183.All the simulation results overestimate the tangential velocity at location r/T=0.183,but the present LES resultis better than thatof RANS and LES(s).Fromthe comparison,we can see that the results of RANS are obviously worse than those from LES,especially for the tangentialvelocity at r/T=0.25 and r/T=0.317.LES with different subgrid-scale models get similar results.Although the advantage is not obvious,the results from the present LES is a bit better in general.For the present LES,the radialvelocity at r/T=0.183 and the tangentialvelocity at r/T=0.317 are almostidenticalto the experimentaldata,while bigger deviation is seen for the LES with standard subgrid-scale model(LES(s)).

Fig.5.Phase averaged plot of the velocity vector field and turbulent kinetic energy.(a)Experiment,phase averaged turbulent kinetic energy are only measured in the impeller region,(b)RANS,(c)LES with Smagorinsky model(LES(s)),(d)Present LES model.(a)–(c)are taken from Hartmann etal.[38].

Fig.6.Phase averaged axialpro file of the radial(left)and tangential(right)velocity components atthree differentradiallocations r/T=0.183(top),r/T=0.25(middle)and r/T=0.317(bottom).

Fig.7 shows the axialpro file of the turbulentkinetic energy atthree different radiallocations.The deviation of the kinetic energy between the simulation and the experiment is much larger than that ofvelocity.Underestimation is observed for allthe simulations at the impeller tip(i.e.r/T=0.183).RANS predictions signi ficantly underestimate the kinetic energy at allthe three locations,mesh re finements are not expected to improve substantially the prediction[10].LES with the standard subgrid-scale model overestimate at locations away from the impeller region(i.e.r/T=0.25 and r/T=0.317).While the peak values of the kinetic energy predicted by the present LES are comparable with that of the experiment at locations away from the impeller region(i.e.r/T=0.25 and r/T=0.317),but the width of the pro file is smaller.

From this study,we can find that the predictions of the present LESIB method compare satisfactorily in generalagainst the experimental data,and it is at least as accurate as other LES-IB methods.It suggests that this method can be promising for the simulation of stirred tanks atlarge scale with complex geometries and athigh Reynolds numbers.

4.Application of the Method to Particle-resolved Simulation of Stirred Tanks

Solid suspensions in liquid are widely encountered in technical systems.There is no universalapproach for modeling and simulation ofsuch problems,the appropriate method depends on the flow regime and on the required details and accuracy[47],butlimited by the computationalcost,Eulerian–Eulerian(two fluid model,TFM)and Eulerian–Lagrangian(discrete particle model,DPM)methods are used most often[48,49].However,empirical correlations are needed to account for unresolved parts of the suspensions physics in these two models and we cannotobtain the complex flow field around the solid particles.

Discrete element method(DEM),originally developed by Cundall and Strack[50],is a soft particle modelwhich allows overlap between the contacting particles[50,51].A linear spring-dashpot force modelis used in the implementation.In this study,IBMis coupled with the particle simulator to handle the solids suspension inside stirred tanks.The motion of particles and the particle–particle collisions are treated by DEM[51],while no-slip conditions at the surfaces of the particles are assigned by using IBM.This method is usually called particle-resolved simulation[49],where the Eulerian grid width is approximately smaller by an order ofmagnitude.In this method,no correlations are required,both the particle–particle and particle–fluid interactions are modeled in a realistic way.There is obviously a computationalpenalty for resolving down to finer scales and the physicalsize of the domain thatcan be simulated gets limited.As a result,smaller stirred tanks compared to laboratory scale are studied here.

Fig.7.Phase averaged axialpro file of the kinetic energy of the random velocity fluctuations at three different radiallocations r/T=0.183,r/T=0.25 and r/T=0.317.

It should be mentioned that the lubrication force is also considered in the collisions to account for condition when particles get in very close proximity[52].During the collision,a thin lubrication layer is formed between the particles and the fluid is squeezed out ofthis gap when the particles approach and is pushed back into the gap during rebound.Viscous forces hence become important in the approaching and rebounding process and can lead to sizeable dissipation.This feature mightbe resolved by an adaptive localgrid re finement[53],butthis usually results in substantially increased computation time.An appropriate kind of“subgrid-scale model”accounting for the film between the surfaces should be implemented when it is too thin to be resolved by the Eulerian grid.This contribution can be quantitatively described by[49,54]

where n is the unit vector from particle i to particle j,Δuij=ui-uj,s0=0.2R0and s1=2×10-4R0.

The physical parameters in the computation are identical to the experiments[55]and the tank size is smaller than that in experiment in consideration of the computation cost.A schematic diagram of the stirred tank is shown in Fig.8,and the numericaland physicalparameters are shown in Table 2.

Three cases with different geometric con figurations are carried out in this study,the setup are listed in Table 3.The average solid volume fractions in these three cases are the same,?=10.4%.The no-slip boundary condition is enforced on all surfaces including the impeller,the tank wall,the baf fles and the particles using the hybrid method of this study.It is dif ficult to judge the relative location between the particles and the impeller,for the convenience to deal with particle-impeller collisions,the impeller is represented by the combination of particles with diameter of d=1 mm.As a result,there are two sets of representations of impeller:combination of particles during particle-impeller collision;and Lagrangian markers during fluid–solid interactions.The particles are initially located at both the upside and the underneath of the impeller,as shown in Fig.9.

The computationalgrid in this simulation is 512×512×512,the resolution of the present grid is not fine enough for DNS,as a result LES with wall-layer model is applied in the simulations.A resolution of16 grid spacing over a particle diameter is used to resolve the flow around a solid sphere.CFL is around 0.2 and it takes about 4177 timesteps for a fullrevolution.

Fig.8.Schematic illustration of the lab-scale stirred tank.Left:side view,right:top view.D=24 mm,b=9.17 mm.

Table 2 The numericaland physicalparameters in the simulation ofsolid suspension in a lab-scale stirred tank

Table 3 The setup of the three cases in the simulation ofsolid suspension in a stirred tank

We investigated the startup of the particle suspension in the three cases.As shown in Fig.9,we can conclude from the qualitative analysis that PBTU is more conducive to particle suspension than PBTD,and more uniform particle distribution in the radial direction is obtained when baf fles are used.This difference is shown more clearly in Fig.10 with the distribution of solids concentration,and more quantitatively in Fig.11 using a globaluniformity index(ξ)[56]where ?iis the solid fraction in i th grid.The index is conveniently defined so thatasξ→0 the uniformity of the suspension isatitsminimum;and whenξ→1 the solids are uniformly distributed within the vessel.

Fig.9.Evolution of the solids suspension process:instantaneous realizations ofparticles locations at four different moments as indicated for the three cases.

Fig.10.Time-averaged solids volume fraction in the verticalmiddle plane(left)and in the horizontalplane with z/H=0.1(right).

Fig.11.Time variation ofglobaluniformity index in the three cases.

Fig.12 shows the phase averaged velocity of the liquid phase,solid phase and the solid velocity between the two.The averaged solid velocity is calculated byUsing cylindricalcoordinates,the radial–axial2D velocity maps were obtained by azimuthally averaging the 3D velocity data and projecting them onto the 2D radial–axialplane.Whilst being two-dimensional,these flow patterns for both the fluid and the solid are qualitatively consistent with the experimentalresults of Guida et al.[55],which can be considered as a veri fication of the present results.

5.Conclusions

An ef ficientmethod for the DNS/LES ofcomplex flow in stirred tanks on fixed Cartesian grids is developed in the framework of the finite difference method.The moving and stationary boundaries are handled respectively by the direct forcing method using Lagrangian markers and a method based on solid volume fractions.This hybrid technique is more ef ficient than using a single method applied to both stationary and moving boundaries.

Fig.12.Phase averaged velocity of the liquid phase,solid phase and the difference between the two.

Meanwhile,the method is furtherimplemented on a supercomputer with CPU–GPU hybrid architecture.The existence of the IB complicates the parallelization.To minimize the amount ofinformation that needs to be communicated through the inter-subdomains,an easy partialvelocity interpolation and force distribution strategy is developed for the parallelization of the IB,which greatly reduces the communication and simpli fies the implementation procedure.

The accuracy and robustness of the present method were validated in two relatively simple simulations:fluid flow in a simple stirred tank and turbulent flow in a Rushton tank.Both validation cases demonstrated,respectively,the ef ficiency of the hybrid IB technique and the capability of the combined LES and IBM implementation.Furthermore,hybrid IB technique is coupled with DEMto study the particle suspension inside the stirred tanks.In this method,no correlations are required;both the particle–particle and particle–fluid interactions are modeled in a realistic way.As a result,this method can be applied to study the fundamental mechanism of solid suspensions in complex flows.

In summary,the presentmethod can be used as a promising toolfor simulating turbulent flows in stirred tanks with complex geometries.

Nomenclature

[1]E.S.Szalai,P.Arratia,K.Johnson,F.J.Muzzio,Mixing analysis in a tank stirred with Ekato Intermig? impellers,Chem.Eng.Sci.59(2004)3793–3805.

[2]R.Verzicco,M.Fatica,G.Iaccarino,P.Orlandi,Flow in an impeller-stirred tank using an immersed-boundary method,AIChE J.50(2004)1109–1118.

[3]J.H.Ferziger,M.Peric,Computationalmethods for fluid dynamics,third ed.Springer,Berlin;New York,2002.

[4]W.J.Gordon,C.A.Hall,Construction of curvilinear co-ordinate systems and applications to mesh generation,Int.J.Numer.Methods Eng.7(1973)461–477.

[5]J.Y.Luo,R.I.Issa,A.D.Gosman,Prediction ofimpeller induced flows in mixing vessels using multiple frames of reference,Eighth European conference on mixing 1994,pp.549–556.

[6]J.Y.Luo,A.D.Gosman,R.I.Issa,J.C.Middleton,M.K.Fitzgerald,Full flow-field computation of mixing in baf fled stirred vessels,Chem.Eng.Res.Des.71(1993)342–344.

[7]D.Deglon,C.Meyer,CFD modelling of stirred tanks:Numerical considerations,Miner.Eng.19(2006)1059–1068.

[8]F.Kerdouss,A.Bannari,P.Proulx,CFD modeling ofgas dispersion and bubble size in a double turbine stirred tank,Chem.Eng.Sci.61(2006)3313–3322.

[9]M.Vakili,M.Esfahany,CFD analysis of turbulence in a baf fled stirred tank,a threecompartment model,Chem.Eng.Sci.64(2009)351–362.

[10]K.Ng,N.Fentiman,K.Lee,M.Yianneskis,Assessment of sliding mesh CFD predictions and LDA measurements of the flow in a tank stirred by a Rushton impeller,Chem.Eng.Res.Des.76(1998)737–747.

[11]M.Sommerfeld,S.Decker,State of the art and future trends in CFD simulation of stirred vesselhydrodynamics,Chem.Eng.Technol.27(2004)215–224.

[12]F.Thibault,P.A.Tanguy,Power-draw analysis ofa coaxialmixer with Newtonian and non-Newtonian fluids in the laminar regime,Chem.Eng.Sci.57(2002)3861–3872.

[13]C.S.Peskin,Flow patterns around heart valves:A numericalmethod,J.Comput.Phys.10(1972)252–271.

[14]J.Kim,An immersed-boundary finite-volume method for simulations of flow in complex geometries,J.Comput.Phys.171(2001)132–150.

[15]E.A.Fadlun,R.Verzicco,P.Orlandi,J.Mohd-Yusof,Combined immersed-boundary finite-difference methods for three-dimensionalcomplex flow simulations,J.Comput.Phys.161(2000)35–60.

[16]M.Uhlmann,An immersed boundary method with direct forcing for the simulation of particulate flows,J.Comput.Phys.209(2005)448–476.

[17]Z.-G.Feng,E.E.Michaelides,Heattransfer in particulate flows with Direct Numerical Simulation(DNS),Int.J.Heat Mass Transf.52(2009)777–786.

[18]N.G.Deen,S.H.L.Kriebitzsch,M.A.van der Hoef,J.A.M.Kuipers,Direct numerical simulation of flow and heat transfer in dense fluid–particle systems,Chem.Eng.Sci.81(2012)329–344.

[19]B.E.Grif fith,R.D.Hornung,D.M.Mcqueen,C.S.Peskin,An adaptive,formally second order accurate version of the immersed boundary method,J.Comput.Phys.223(1)(2010)10–49.

[20]R.Mittal,G.Iaccarino,Immersed boundary methods,Annu.Rev.Fluid Mech.37(2005)239–261.

[21]J.G.M.Eggels,Direct and large-eddy simulation of turbulent fluid flow using the lattice-Boltzmann scheme,Int.J.Heat Fluid Flow 17(1996)307–323.

[22]J.Derksen,H.E.A.Van den Akker,Large eddy simulations on the flow driven by a Rushton turbine,AIChE J.45(1999)209–221.

[23]F.Sbrizzai,V.Lavezzo,R.Verzicco,M.Campolo,A.Soldati,Direct numericalsimulation ofturbulent particle dispersion in an unbaf fled stirred-tank reactor,Chem.Eng.Sci.61(9)(2006)2843–2851.

[24]D.Goldstein,R.Handler,L.Sirovich,Modeling a no-slip flow boundary with an externalforce-field,J.Comput.Phys.105(1993)354–366.

[25]T.Kajishima,S.Takiguchi,H.Hamasaki,Y.Miyake,Turbulence structure ofparticleladen flow in a verticalplane channeldue to vortex shedding,JSME Int.J.Ser.B Fluids Therm.Eng.44(2001)526–535.

[26]Mohd-Yusof,Combined immersed boundary B spline methods for simulations of flow in complex geometries,Annu.Res.Briefs(1997)317–327.

[27]X.Yang,X.Zhang,Z.Li,G.-W.He,A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations,J.Comput.Phys.228(2009)7821–7836.

[28]T.Kajishima,S.Takiguchi,H.Hamasaki,Y.Miyake,Turbulence structure ofparticleladen flow in a verticalPlane Channeldue to vortex shedding,JSME Int.J.Ser.B Fluids Therm.Eng.44(2001)526–535.

[29]J.J.J.Gillissen,H.E.A.Van den Akker,Direct numerical simulation of the turbulent flow in a baf fled tank driven by a Rushton turbine,AIChE J.58(2012)3878–3890.

[30]M.Tyagi,S.Roy,A.D.H.Iii,S.Acharya,Simulation of laminar and turbulent impeller stirred tanks using immersed boundary method and large eddy simulation technique in multi-block curvilinear geometries,Chem.Eng.Sci.62(2007)1351–1363.

[31]J.K.Wiens,J.M.Stockie,An ef ficient parallelimmersed boundary algorithm using a pseudo-compressible fluid solver,J.Comput.Phys.281(2013)917–941.

[32]E.Givelberg,K.Yelick,Distributed immersed boundary simulation in titanium,SIAM J.Sci.Comput.28(2006)1361–1378.

[33]H.Hartmann,J.J.Derksen,C.Montavon,J.Pearson,I.S.Hamill,H.E.A.van den Akker,Assessmentoflarge eddy and RANS stirred tank simulations by means ofLDA,Chem.Eng.Sci.59(2004)2419–2432.

[34]R.M.Jones,A.D.Harvey,S.Acharya,Two-equation turbulence modeling for impeller stirred tanks,J.Fluids Eng.123(2001)640.

[35]S.B.Pope,Turbulent flows,Cambridge University Press,Cambridge;New York,2000.

[36]J.Smagorinsky,Generalcirculation experiments with the primitive equations,I.The basic experiment,Mon.Weather Rev.91(1963)99–164.

[37]C.Meneveau,T.S.Lund,W.H.Cabot,A Lagrangian dynamic subgrid-scale modelof turbulence,J.Fluid Mech.319(1996)353–385.

[38]H.Hartmann,J.J.Derksen,C.Montavon,J.Pearson,I.S.Hamill,H.E.A.V.D.Akker,Assessment of large eddy and RANS stirred tank simulations by means of LDA,Chem.Eng.Sci.59(2419–2432)(2004).

[39]C.S.Peskin,The immersed boundary method,Acta Numer.11(2003)479–517.

[40]B.K.Cook,A numericalframework for the direct simulation of solid–fluid systems,Massachusetts Institute of Technology,2001.

[41]J.O'Rourke,Computationalgeometry in C,second ed.Cambridge University Press,Cambridge;New York,1998.

[42]W.Gropp,E.Lusk,Using MPI-2:a problem-based approach,Recent advances in parallelvirtualmachine and message passing Interface,Proceedings,3666 2005,p.8.

[43]C.Ji,A.Munjiza,J.J.R.Williams,A noveliterative direct-forcing immersed boundary method and its finite volume applications,J.Comput.Phys.231(2012)1797–1821.

[44]S.Di,W.Ge,Simulation of dynamic fluid–solid interactions with an improved direct-forcing immersed boundary method,Particuology 18(2015)22–34.

[45]L.Dong,S.T.Johansen,T.A.Engh,Flow induced by an impeller in an unbaf fled tank—I.Experimental,Chem.Eng.Sci.49(1994)549–560.

[46]A.Bakker,L.M.Oshinowo,Modelling of turbulence in stirred vessels using large Eddy simulation,Chem.Eng.Res.Des.82(2004)1169–1178.

[47]J.Derksen,Highly resolved simulations of solids suspension in a smallmixing tank,AICHE J.58(2012)3266–3278.

[48]J.Li,W.Ge,W.Wei,X.L.Ning Yang,L.Wang,X.He,X.Wang,J.Wang,M.Kwauk,From multiscale modeling to meso-science:a chemical engineering perspectiveprinciples,modeling,simulation and applications,Springer,London,2013 476.

[49]M.A.van der Hoef,M.van Sint Annaland,N.G.Deen,J.A.M.Kuipers,Numericalsimulation ofdense gas–solid fluidized beds:Amultiscale modeling strategy,Annu.Rev.Fluid Mech.40(2008)47–70.

[50]P.A.Cundall,O.D.L.Strack,A discrete numerical-model for granular assemblies,Geotechnique 29(1979)47–65.

[51]J.Xu,H.Qi,X.Fang,L.Lu,W.Ge,X.Wang,M.Xu,F.Chen,X.He,J.Li,Quasi-real-time simulation of rotating drum using discrete element method with parallel GPU computing,Particuology 9(2011)446–450.

[52]T.Kempe,J.Fr?hlich,Collision modelling for the interface-resolved simulation of sphericalparticles in viscous fluids,J.Fluid Mech.709(2012)445–489.

[53]H.H.Hu,Directsimulation of flows ofsolid–liquid mixtures,Int.J.Multiphase Flow 22(1996)335–352.

[54]T.Kempe,J.Fr?hlich,An improved immersed boundary method with direct forcing for the simulation of particle laden flows,J.Comput.Phys.231(2012)3663–3684.

[55]A.Guida,A.W.Nienow,M.Barigou,PEPT measurements of solid–liquid flow field and spatial phase distribution in concentrated monodisperse stirred suspensions,Chem.Eng.Sci.65(2010)1905–1914.