Free Surface Turbulent Flow in an Unbaffled Stirred Tank : Detached Eddy Simulation and VOF Study

Stirred tanks are widely used in process industries for blending of singleand/or multiple-phase fluids, which often involves mass and/or heat transfer, chemical reactions, etc. Deep knowledge of flow field in the stirred tank is necessary to carry out these unit operations efficiently. As a result, extensive investigations of the flow in stirred tanks have been conducted over the past several decades. However, this does not mean that we have completely understood the hydrodynamics in stirred tanks. Generally speaking, the fluid flow in a stirred tank is highly complicated. It is more complex when the free-surface deformation is taken into account. Up to now, most published studies on the hydrodynamics in stirred tanks have been carried out in baffled stirred tanks. In this case, the liquid surface at the top of the stirred tank was commonly assumed to be flat (see for example, Alcamo et al., 2005; Armenante et al., 1997; Dong et al., 1994; Shekhar and Jayanti, 2002).1–4 But in fact, the liquid surfaces in baffled stirred tanks are wavy. There is surface macro-swell, i.e., elevation/depression of the liquid surface with time (Jahoda et al., 2011).5 In some process industries, such as animal cell cultures, micromixing (Aloi and Cherry, 1996; Rousseaux et al., 2001; Assirelli et al., 2008)6–8, the use of unbaffled stirred tanks may be desirable. In such cases, surface vortices are formed and the liquid surface can no longer be treated as flat. Accordingly, the free surface deformation must be modelled to obtain a more accurate result. The volume of fluid (VOF) method has been widely used for the modeling of fluid flow in stirred tanks. The first application was performed by Serra et al. (2001).9 They simulated the flows with wavy free-surface in a fully baffled stirred tank. Haque et al. (2006 and 2011)10–11 for the first time used an Eulerian-Eulerian multiphase flow model coupled with VOF to determine the gas-liquid interface in unbaffled vessels. The predicted liquid surface profiles, the mean flow characteristics and the turbulent kinetic energy were generally well predicted. However, some differences between the numerical predictions and experimental results could still be observed. They attributed this to the weaknesses of the RANS turbulence models they had used. Glover and Fitzpatrick (2007)12 and Torré et al. (2007)13 used the same approach for the modelling of vortex in an unbaffled stirred tank. They also pointed out Free Surface Turbulent Flow in an Unbaffled Stirred Tank: Detached Eddy Simulation and VOF Study


Introduction
Stirred tanks are widely used in process industries for blending of single-and/or multiple-phase fluids, which often involves mass and/or heat transfer, chemical reactions, etc. Deep knowledge of flow field in the stirred tank is necessary to carry out these unit operations efficiently.As a result, extensive investigations of the flow in stirred tanks have been conducted over the past several decades.However, this does not mean that we have completely understood the hydrodynamics in stirred tanks.Generally speaking, the fluid flow in a stirred tank is highly complicated.It is more complex when the free-surface deformation is taken into account.
Up to now, most published studies on the hydrodynamics in stirred tanks have been carried out in baffled stirred tanks.2][3][4] But in fact, the liquid surfaces in baffled stirred tanks are wavy.There is surface macro-swell, i.e., elevation/depression of the liquid surface with time (Jahoda et al., 2011). 5n some process industries, such as animal cell cultures, micromixing (Aloi and Cherry, 1996  [6][7][8] , the use of unbaffled stirred tanks may be desirable.In such cases, surface vortices are formed and the liquid surface can no longer be treated as flat.Accordingly, the free surface deformation must be modelled to obtain a more accurate result. The volume of fluid (VOF) method has been widely used for the modeling of fluid flow in stirred tanks.The first application was performed by Serra  et al. (2001). 9They simulated the flows with wavy free-surface in a fully baffled stirred tank.Haque et  al. (2006 and 2011)  [10][11] for the first time used an Eulerian-Eulerian multiphase flow model coupled with VOF to determine the gas-liquid interface in unbaffled vessels.The predicted liquid surface profiles, the mean flow characteristics and the turbulent kinetic energy were generally well predicted.However, some differences between the numerical predictions and experimental results could still be observed.They attributed this to the weaknesses of the RANS turbulence models they had used.Glover and Fitzpatrick (2007) 12 and Torré et al. (2007) 13 used the same approach for the modelling of vortex in an unbaffled stirred tank.They also pointed out that large eddy simulation (LES), detached eddy simulation (DES), scale adaptive simulation (SAS) or Reynolds stress model (RSM) that can improve the calculation accuracy of turbulence quantities of the flow fields should be adopted.The combination of RSM and VOF was first employed by Mahmud  et al. (2009)  14 to simulate the free-surface turbulent flow in an unbaffled reactor agitated by a cylindrical magnetic stirrer.The predicted free-surface shape was in good agreement with measurements, but the vortex depth was under-predicted.They did not give the reason, and we think this may be due to the weakness of the RSM model they had used.In a recent study, VOF together with the LES model was used for prediction of liquid flow and free liquid surface motion in a stirred tank by Jahoda et al.  (2011). 5The predicted results were found to be in good agreement with an experimental investigation by a conductivity method.The VOF model was also used by Motamedvaziri and Armenante (2012) 15 to model the air-water interface in stirred vessels for different liquid fill ratios.The vortex formation and air entrainment phenomenon were well predicted.
In the present work, numerical simulations of the free-surface hydrodynamics in an unbaffled dished-bottom stirred tank was performed by using the combination of the VOF multiphase model and the DES model.In order to validate the numerical methods, the predictions were compared with data reported by Haque et al. (2011)  11 .DES combines the advantages of RANS and LES.It can resolve the flow field at a low cost and high accuracy.[18][19]

Stirred vessel configuration
As shown in Fig. 1, the stirred tank is an unbaffled, dished-bottom, cylindrical tank agitated by a standard Rushton impeller.The tank diameter is where C DES = 0 65 .is the model empirical constant, and D is the largest dimension of the grid cell.This modification of the S-A model changes the interpretation of the model substantially.In regions close to the wall, where d C < DES ∆ , it behaves as a RANS model.Away from the wall, where d C > DES ∆ , it behaves in a Smagorinsky-like manner and is changed to the LES model.
The governing equation of DES model can be given as follows: The methodology proposed here has already been used in our previous work.For the purpose of simplicity, this is not given here.For more information, please refer to Yang et   Modeling of the free surface turbulent flow was performed by using the VOF method.The mass and momentum conservation equations of the VOF model are given as: where t is time, u is velocity, p is pressure, g is gravitational acceleration, F is the force acted on the fluid element, ρ is the volume-averaged density of the fluids in the tank, μ is the volume-averaged viscosity: where α a and α w are volume fraction of air and water, and can be determined by solving a continuity equation for one (or more) of the phases.
For the i-th phase, this equation has the following form: where α i is α a or α w and: In a given computational cell, the volume fraction of unity represents pure water and zero represents pure air.Accordingly, the air-water interface can be determined by identifying the cells where the volume fraction is 0 1 < < α i .The continuum surface force (CSF) model was used to calculate the surface tension.The geometric reconstruction scheme based on piecewise linear interpolation was used for the reconstruction of the interface.

Computational details Computational domain and grid arrangement
Although initial water level in the vessel was 1 T, the height of the computation domain was extended up to 1.5 T to capture the air-water interface.The sliding mesh method was used for modelling the relative rotation between the moving impeller and the stationary tank.An inner rotating cylindrical volume centered on the impeller (rotor region) and an outer stationary zone containing the rest of the tank (stator region) were defined.The origin of coordinate system coincides with the center of the stirred tank bottom.The boundary of the rotor region was positioned at r = 60 mm and 0 026 .m 0.076 m < < z (where z is the axial distance from the bottom of the tank).
The grid was prepared with the pre-processor Gambit 2.3.The computational domain was discretized using a non-uniformly distributed unstructured mesh.The maximum skewness of the volume mesh was less than 0.84.Flow in the impeller region is the most complicated.In order to resolve the flow field in this region accurately, more nodes were generated.Along the impeller width, 38 nodes were assigned with the minimal grid length equalling 0.5 mm.In addition, grid independence tests were conducted to examine the effect of mesh size on the numerical predictions.Three grids composed of approximately 1,100,000, 1,408,000 and 1,601,000 cells were generated and referred to as coarse, medium, and fine grid, respectively.All the grids were generated according to the guidance of Spalart (2001&2009) 20-21 with y + ranged from 1 to 20: where u T is the friction velocity, y is the distance from the nearest wall, and ν is the kinematic viscosity.The sensitivity of the simulation to the grid resolution was verified by comparing local velocities of the flow fields calculated on the three grids.It was found that the velocities obtained from the simulation using the fine grid were nearly identical to those obtained using the medium grid, and accordingly, the medium grid was used to perform the numerical simulation.

Initial and boundary conditions
The initial velocities of the fluids in the stirred tank were assumed to be zero.The upper 1/3 region of the computational domain was filled with compressible air.The remaining regions were patched with water with a volume fraction of one.The initial interface between air and water was assumed to be flat.No-slip wall boundary conditions were applied to all solid surfaces.The impeller rotates with the rotor and for the shaft in the stator region, it rotates with the same speed as the impeller.The zero-shear boundary condition is used at the top of the computation domain in the gas phase.

Modeling method
The ANSYS FLUENT 13.0 solver was used to carry out the numerical simulations.The pressure-based Navier-Stokes algorithm was used for the solution of the model with implicit solver formulation where the absolute velocity formulation was adopted.To ensure smooth and better convergence, initially the k-ε computation was performed until the steady state flow field was obtained.Subsequently, the result of the steady-state computation was used as the initial solution to perform the unsteady DES computation.
For the k-ε simulation, the standard wall functions were used to solve the near-wall flow.The SIMPLE algorithm was performed to couple the velocities and pressure terms.The continuity equation, momentum conservation equation, and k-ε equation were all discretized using the second order upwind scheme.For the DES computation, the PISO algorithm was adopted to couple the velocities and pressure terms.The first order implicit scheme was used for time discretization.An explicit VOF scheme with the implicit body force formulation was selected with the Geo-Reconstruct scheme to capture the interface.The body force weighted scheme was employed to compute the pressure.The bounded central differencing scheme was adopted for spatial discretization of momentum and the modified turbulent viscosity equation.For the VOF multiphase model coupled with level-set method, the compressible air was set as the primary phase, and water and tracer were set as the secondary phase.This treatment can improve the solution stability.The surface tension coefficient σ for the interface between air and water was specified as a constant which equaled 0.072 N m -1 .Gravity was defined as 9.81 m s -2 , along the negative Z direction.The specified operating density was set as that of air, which was 1.225 g m -3 .
The time step is crucial to the transient DES modeling.The time step must be small enough to capture the flow features induced by the motion of the impeller.Furthermore, it also must be considered with the grid to ensure a stable and converged solution.In the present study, the time step was set as follows: At the initial stage of the computation, the fixed time stepping method with a time step of 1• 10 -6 s was chosen for the VOF calculation.The solution was considered to be fully converged when the normalized residuals for all variables were less than 1• 10 -4 .Then the variable time stepping method was chosen to accelerate the calculation.The time-step change factor was limited in the range of 0.5-2 to avoid the sudden increase in time step size.The global Courant number was controlled below 2.

Results and discussion
Profiles of the liquid surface Typical time variations of the numerically predicted water surfaces, which are presented in the form of water iso-surface volume fractions, are shown in Fig. 2. For comparison, the predictions based on the DES and k-ε model were both presented.The influence of the value of the water volume fraction used to determine the free surface shape has been studied by Torré et al. (2007). 22They pointed out that small gradient of water volume fraction profiles around the air/water interface is sufficient to give a difference in the free surface representation.A water volume fraction of 0.9 was suggested to stand for the existence of water surface, and this value was adopted here.
In the absence of experimentally measured free surface profiles, the numerically predicted water surfaces were compared with the flow patterns reported by Haque et al. (2011). 11It has been generally recognized that the numerically determined depths of the free surfaces in concentrically unbaffled tanks based on RANS models, such as k-ε model, RST model and/or SST model, were all found to be under-predicted.This may be attributed to the drawbacks of these models.However, the free surface profiles predicted by these models qualitatively agree well with that estimated using a correlation given by Markopoulos and Kontogeorgaki (1995). 23In Fig. 2, we can see that the water surface profiles are unsteady due to the transient nature of the free-surface hydrodynamics in the stirred tank.The results obtained by using k-ε model are consistent with the results of Haque et al. (2011) 11 , which indicates that our simulation strategies are reliable.By comparison, the profiles predicted by the DES model are deeper, which are thought to be much closer to the experimental results.In this sense, the combination of DES model and VOF technique, and the associated numerical strategies adopted in this study, can predict the transient free liquid surface in stirred tanks with considerable accuracy.

Instantaneous quantities
The mean quantities given in the following sections were all obtained from a pseudo-steady simulation.Here, attainment of the pseudo-steady condition is demonstrated by plotting the time-dependent instantaneous axial, radial, and tangential velocity values acquired at r = 0.421D, θ = 30° and z = 0.061 m.For clarity of representation, a time interval of 20 impeller revolutions, which corresponds to 6 s, was tracked.The results are shown in Fig. 3.It can be observed that the turbulent fluctuations are well es-tablished after about 5 impeller revolutions.Besides, the instantaneous velocities oscillate with time.This is named as the macro-instability phenomena.Frequency analysis was applied to the velocity recordings based on fast Fourier transforms (FFT).The frequency spectra are presented in Fig. 4. From these plots, a main frequency component of 0.11 Hz can be observed.

Mean velocity distributions
In order to make a quantitative analysis, mean velocity profiles in the vertical plane θ = 0° are given below.For comparison purposes, the results obtained by Haque et al. (2011) 11 are also presented.In these plots, the radial position r is normalized by the radius R of the stirred tank, and the mean velocity is normalized by the blade tip velocity u tip .DES computation can only provide instantaneous velocities and it is necessary to carry out post-processing.The DES results were averaged only after the convergence of the instantaneous velocity was ob-tained.Subsequently, the mean velocity components were angular-averaged to compare with the LDV results.
The mean axial (u z ), radial (u r ), and tangential (u θ ) velocity profiles at three axial heights, z = 24, 41, and 61 mm, are shown in Figures 5-7, respectively.The first and second measurement positions are below the impeller and the last position is above the impeller.In Fig. 5, the numerically predicted radial profiles of the axial velocity are depicted and compared with the experimental results.Good agreement between the DES and LDV results of the axial velocity can be clearly observed from these plots.When compared with the k-ε results, the DES predictions are more close to the LDV data, especially in the regions below the impeller.At all locations, the average difference between DES and LDV results is within 10 %.Therefore, it could be said that the methodology here applied has yielded satisfactory agreement between experimental and simulated values.The radial velocity is shown as radial profile at three axial levels in Fig. 6.Generally speaking, the DES results agree well with the LDV data.Although they were under-predicted, the velocity profiles tendency is much the same.also simulated the mean velocities with SST and RST model.They found that the radial velocities were all significantly under-predicted.In addition, it should be pointed out that the magnitudes of the radial velocities are much smaller than the tangential velocities, and therefore the accuracy may be debatable.Fig. 7 reports radial profiles of the tangential velocity at the selected three axial levels.It is evident that of all the velocities, the tangential velocity is much larger over the bulk of the flow field.This is true for the fluid flow in unbaffled stirred tank.For this component, excellent agreement between the numerical results (both DES and k-ε) and LDV data can be achieved.By comparison, DES gives more accurate predictions.The maximum difference is no more than 5 %.

Turbulent kinetic energy distributions
In incompressible flows, the instantaneous velocity can be decomposed into three parts: a time average component, a deterministic periodic oscillation component, and a turbulent fluctuation component: Accordingly, a periodic kinetic energy and a turbulent kinetic energy can be determined.In the reference paper of Haque et al. (2011) 11 , the turbulent kinetic energy was obtained from the phase-averaged measurement of rms velocities.Following the same approach, the turbulent fluctuating component of the velocities was used to compute the turbulent kinetic energy using the following expression: The quantitative comparisons between the numerically predicted and experimentally measured turbulent kinetic energies are shown in Fig. 8.It is  clear that DES predictions agree well with the LDV results.Best agreement can be found at z = 24 mm, which is in the region below the impeller.At another two positions close to the impeller, i.e., z = 41 mm and 61 mm, the agreement is not as good.
Overall, the turbulent kinetic energies are somewhat under-predicted.In most locations, the error is within 10 %.However, at z = 61 mm and r/R = 1/3, the experimental result was found to be nearly 4 times the numerical prediction.It needs to be stated that, the experimentally determined turbulent kinetic energy consists of a coherent part due to the periodic fluctuations related with the impeller blade passage frequency, and a random part induced by random fluctuations associated with turbulence.Unfortunately, the random part of the turbulent kinetic energy was not extracted in the experimental study.This may partly account for the under-predictions of numerical results with the experimental data.By contrast, the k-ε model can only predict the distribution tendency.At most positions, turbulent kinetic energies were all under-predicted.This finding indicates that DES works better than k-ε in the prediction of turbulent kinetic energy in the unbaffled stirred tank.

Conclusions
The free-surface turbulent flow in an unbaffled dished-bottom stirred tank was modelled using a combination of DES and VOF method.The mean velocity and turbulent kinetic energy were determined and compared with the LDV measurements and k-ε predictions reported by Haque et al. (2011)  11 .DES predictions of the free surface profiles were compared with the k-ε simulations.In general, it is confirmed that the DES model gives better predictions than the k-ε model.Combining DES with VOF models produces predicted vortex shapes "closer to reality".In terms of angular-averaged instantaneous velocity, the mean velocity and turbulent kinetic energy generated by DES based on Spalart-Allmaras model can be well predicted.More specifically, the mean radial velocity was not predicted with satisfactory accuracy.The error is about 50 %.This may be attributed to the insufficient grid resolution and the inaccurate measurements of the radial velocities, which has been pointed out by Haque et al. in their study.By comparison, the errors are much smaller for the mean axial and tangential velocity, which are no more than 10 % and 5 %, respectively.Overall, turbulent kinetic energies are also well predicted with the error less than 10 %. 21306105) is gratefully acknowledged.
The turbulent flow field was modelled by the DES model.It is formulated by replacing the distance function d in the Spalart-Allmaras (S-A) model with a modified distance function given below:

F i g . 1 -
Schematic diagram of the stirred systemPrediction of free-surface flow
Haque et al. (2011) 11 pointed out that their measurements of the radial velocities are much larger than those reported by others (see for example, Alcamo et al., 2005; Torré et al., 2007) 1,22 .This may be one reason why such a large error was detected.Besides the k-ε model, Haque et al. (2011)

F i g . 3 -
Time dependence of the instantaneous (a) axial, (b) radial, and (c) tangential velocities acquired at r = 0.421D, θ = 30° and z = 0.061 m F i g .4 -Frequency spectra of (a) axial, (b) radial, and (c) tangential velocity recordings acquired at r = 0.421D, θ = 30° and z = 0.061 m F i g .5 -Comparison between the numerically predicted and experimentally measured radial profiles of the mean axial velocity at three axial heights: z = 24, 41 and 61 mm.The k-ε and LDV data were adopted from Haque et al. (2011) 11

F i g . 6 -
Comparison between the numerically predicted and experimentally measured radial profiles of the mean radial velocity.Same conditions as shown in Fig.5.F i g .7 -Comparison between the numerically predicted and experimentally measured radial profiles of the mean tangential velocity.Same conditions as shown inFig. 5.

F i g . 8 -
Comparison between the numerically predicted and experimentally measured radial profiles of the turbulent kinetic energies.Same conditions as shown in Fig.5.

N 1 u
o m e n c l a t u r e C -impeller off-bottom clearance, m C DES -empirical constant of DES model, dimensionless d -distance function, m  d -modified distance function, m D -impeller diameter, m F -force, Pa g -gravitational force, m s -2 H -liquid level, m k -turbulent kinetic energy, m 2 s -2 N -impeller rotational speed s -1 p -pressure, Pa r -radial coordinate, m Re -Reynolds number, dimensionless t -time, s T -vessel diameter, m u r -mean radial velocity, m s -1 u tip -impeller blade tip speed, m s -1 u z -mean axial velocity, m s -1 u θ -mean tangential velocity, m s -1 u -velocity m s --time average velocity, m s -1  u -deterministic periodic oscillation velocity, m s -1 ′ u -turbulent fluctuation velocity, m s -1 u T -friction velocity, m s -1 y -distance from node to nearest wall, m z -axial coordinate, m G r e e k s y m b o l s α -volume fraction, % Δ -largest dimension of the grid cell, m ε -turbulent dissipation rate, m 2 s -3 θ -tangential coordinate, °μ -dynamic viscosity, Pa s ν -kinematic viscosity, m 2 s -1 ρ -density, kg m -3 σ -surface tension, N m -1 S u b s c r i p t s a -air w -water