Comparison between RANS and 3D-PTV measurements of Newtonian and non-Newtonian fluid flows in a stirred vessel in the transitional regime

(cid:1) 3D-PTV measurements and RANS of transitional ﬂows in a lab-scale stirred vessel have been cross-validated. (cid:1) The turbulence model had little effects on the main ﬂow. (cid:1) The discretisation scheme had signiﬁcant effects on the turbulent quantities. (cid:1) Non-Newtonian rheology exacerbated the difﬁculties associated to RANS of transitional ﬂows.


Introduction
Liquid mixing in mechanically agitated vessels is widely encountered in the formulation and manufacture of complex liquid products across a wide range of industries.Examples of such products include paints and inks, foods, home and personal care products catalyst intermediates and battery materials.This unit operation has received significant research attention, with flow measurement techniques such as Laser Doppler Anemometry (LDA), Particle Image Velocimetry (PIV) and Planar Laser Induced Fluorescence (PLIF), enabling quantitative analysis of mixing time, flow patterns, turnover rates and turbulent fluctuations (Galletti et al., 2003(Galletti et al., , 2004;;Ducci andYianneskis, 2005, 2007;Simmons et al., 2007;Chung et al., 2009;Gabriele et al., 2009;Story et al., 2018;Taghavi and Moghaddas, 2019).
Despite the fact that mixing under transitional flow conditions is very common in industry, especially for products whose rheology is evolving during manufacture, most experimental studies in the literature focus on fully laminar or fully turbulent flows.Recently, the authors of this paper applied 3D Particle Tracking Velocimetry (3D-PTV) to obtain Lagrangian trajectories of flow tracers in a lab-scale vessel at Re ¼ 12; 000 (Romano et al., 2021a), where Re ¼ qD 2 N=l is the impeller Reynolds number, q is the fluid density, l is the viscosity, D is the impeller diameter and N is the rotational speed.They used a polynomial interpolation of the Lagrangian coordinates with respect to time to filter the experimental errors and found, based on the definition of goodness of fit, that a second order polynomial was better than a first or third order.This means that not only the instantaneous velocity, but also the Lagrangian acceleration is required for an accurate description of the tracer trajectories.In addition, higher order derivatives are not necessary.They also extended their study to examine the agitation of Newtonian and non-Newtonian fluids in a vessel operated in the transitional flow regime (Romano et al., 2021b).The analysis of the distributions of Lagrangian acceleration and shear rate within a control volume around the impeller revealed that their average values were strongly correlated over many combinations of Reynolds number and fluid rheology.
Computational Fluid Dynamics (CFD) is a powerful tool for simulating fluid flow applications, including mixing in agitated vessels, and has some significant advantages compared to experiments.Accurate measurements of shear rate, energy dissipation rate and phase dispersion may be unfeasible or too expensive.Numerical simulations can provide this information at high spatial and temporal resolution and the effects of relevant variables can be isolated easily.For instance, the fluid rheological parameters can be changed one at a time.However, numerical simulations always represent a compromise between accuracy and associated computational cost, and validation against experimental data is always necessary.The numerical solution is influenced by many simplifying assumptions, the most impactful being the turbulence modelling approach, the boundary conditions at the impeller and the numerical discretisation scheme.These are discussed below.
CFD techniques are classified into Direct Numerical Simulations (DNS), Large Eddy Simulations (LES) or Reynolds-Averaged Navier-Stokes (RANS) simulations, according to how the turbulent fluctuations of simulated quantities are handled.DNS solve the instantaneous flow equations down to the dissipative scales.This requires extremely fine computational grids and limits DNS to low Re operations (Re 10 3 ).To date, only a few authors have conducted DNS of stirred tanks (Bartels et al., 2002;Sbrizzai et al., 2006;Gillissen and Van den Akker, 2012;Bas ßbug ˘et al., 2017;Tamburini et al., 2018).In LES, the contribution of the large anisotropic eddies to momentum transfer is computed explicitly, while the effect of the isotropic small scales is modelled.This is a compromise between DNS resolutions and RANS affordability.LES of stirred tanks provide better predictions of turbulence quantities than RANS (Yeoh et al., 2004;Alcamo et al., 2005;Delafosse et al., 2008Delafosse et al., , 2009;;Zadghaffari et al., 2010;Kysela et al., 2017) and represent the current best practice.However, they remain unpractical in industrial environments with limited computational budgets.
The RANS approach is the most widespread method for simulating turbulent fluid mixing (Coroneo et al., 2011).The equations for the mean flow are obtained by Reynolds-averaging all terms in the mass and momentum equations.This introduces the need to model the effects of turbulence on the mean flow.Many RANS models have been developed, each one with intrinsic benefits and limitations.Eddy-viscosity models, such as the standard k-e (Launder and Spalding, 1974), the k-x (Wilcox, 1988) and their variants, are based on Boussinesq's hypothesis and cannot predict turbulence anisotropic effects.However, experimental evidence and numerical studies indicate that turbulence in agitated vessels is anisotropic (Derksen et al., 1999;Galletti et al., 2004b;Yeoh et al., 2004;Ducci and Yianneskis, 2005;Escudié and Liné, 2006;Haque et al., 2011;Zamiri and Chung, 2018).As a result, eddyviscosity models applied to stirred tank flows tend to predict the mean velocity accurately, but not the quantities related to turbulence and dissipation (Bartels et al., 2000;Alexopoulos et al., 2002;Bartels et al., 2002;Aubin et al., 2004;Singh et al., 2011;Zhang et al., 2017).
RANS models can also differ in their near-wall formulation and these wall models are at least as important as the turbulence model for the bulk flow (Menter, 2009).For instance, the Shear Stress Transport (SST) model (Menter, 1994) combines the robustness of k-e model in the bulk with the near-wall accuracy of k-x by simply applying a blending function.In some cases, the SST model has been proved superior to k-e for stirred vessels simulations (Singh et al., 2011;Xiao et al., 2014;Lane, 2017;Tamburini et al., 2021).However, as the full resolution of the flow in the boundary layer requires extremely fine meshes, it is common to use wall functions in this region.These are empirical correlations that model the velocity profile based on the law-of-the-wall (Launder and Spalding, 1974) rather than solving it directly from the flow equations.
The impeller rotation with respect to the stationary elements, such as baffles, can be handled with the sliding mesh (SM) or the multiple frames of reference (MFR) approach.The vessel is divided in an outer region, which includes the stationary elements, and an inner cylinder around the impeller.The SM technique (Luo et al., 1993) can be used in transient simulations to resolve the impeller-baffle interactions and pseudo-turbulence.At each time step, the inner region rotates and the solution is recomputed according to the new impeller position.In the MFR approach (Luo et al., 1994), the flow equations in the inner region are solved in a rotating frame of reference, so that both the inner and outer flows are assumed steady in their own reference frame.Although the MFR approach does not capture the impeller-baffle interactions, many researchers have reported no significant differences in the mean velocities and flow patterns obtained between MFR and SM methods (Aubin et al., 2004;Patel et al., 2015;Lane, 2017).However, turbulent quantities are sensitive to the impeller model, the SM being more accurate (Aubin et al., 2004).
The discretisation scheme is the method for interpolating cell face-values from the values at the neighbouring cell centres.First order schemes are numerically stable but also prone to numerical diffusion.Higher order schemes are less sensitive to numerical diffusion, but they may suffer convergence issues.Few authors have investigated the impact of the discretisation scheme on the predicted flows in stirred tanks.In most cases, little effect was observed on the mean velocities, but a first order scheme exacerbated the errors in the turbulent quantities and power number (Aubin et al., 2004;Deglon and Meyer, 2006;Coroneo et al., 2011).Although a myriad of numerical studies of agitated vessel flows are available in the literature, most focus on fully turbulent or laminar conditions or involve Newtonian fluids.There is still a dearth of CFD studies involving non-Newtonian fluids in the transitional flow regime, despite such operations being extremely common in industry.Therefore, the best practice for simulation of these systems are not yet clearly defined.In this paper, steady RANS simulations of an agitated vessel equipped with a Rushton turbine have been conducted in the transitional regime.The effects of the turbulence model and discretisation scheme have been addressed based on two Newtonian fluids of different viscosity.Then, four non-Newtonian fluids, two with shear thinning power law rheology and two with a Herschel Bulkley rheology, have been studied.The numerical solutions have been validated against 3D-PTV measurements obtained by Romano et al. (2021b).Lagrangian acceleration measurements obtained through PTV and shear rate data simulated with CFD have been interpolated in a Eulerian grid to test the spatial correlation between the two quantities within the vessel.

Equipment set-up and flow conditions
The vessel under study had diameter T ¼ 180 mm and height H ¼ T. It was stirred with a Rushton disc turbine of diameter D ¼ T=3 placed at a height C ¼ T=3 from the bottom.The impeller blades were D=4 in width and D=5 in height.The vessel was equipped with four baffles of width w ¼ T=10.The thickness of impeller blades, disc and baffles was 2 mm.The setup is identical that used in Romano et al. (2021b).
The six working fluids investigated were aqueous solutions of glycerol (Sigma-Aldrich, 75 % and 80 % wt.concentration), carboxymethylcellulose sodium salt (medium viscosity CMC, Sigma-Aldrich, 0.5 % and 1.0 % wt.) and crosslinked polyacrylic acid polymer (Carbopol 940, Lubrizol plc, 0.1 % wt. at pH 5.5 and 0.15 % wt. at pH 5.0) at room temperature.The rheological parameters of the liquids are given in Table 1.The glycerol solutions have constant viscosity, l, given the definition of a Newtonian fluid.The aqueous CMC solutions possessed a shear-thinning flow rheology; the apparent viscosity, l a , is a function of the local shear rate, _ c, according to the power-law model, where K is the fluid consistency index and n is the power law exponent.The Carbopol 940 solutions showed a yield stress, s 0 , and their viscosity was expressed with a piecewise Herschel-Bulkley model: In eq. ( 2), l 0 is the zero-shear viscosity, which is used in the CFD solver to avoid the numerical singularity at _ c ¼ 0 (Xiao et al., 2014;Patel et al., 2015).Its value should be large to ensure that the fluid flow is negligible when s s 0 .In the simulations presented in this paper, the zero shear viscosity was set conservatively as viscosity and effective shear rate according to the Metzner-Otto's rule (Metzner and Otto, 1957), and k s ¼ 11:5 is the typical value for a Rushton turbine.

CFD framework
All simulations described below have been carried out in Open-Foam Ò on a cluster node with 32 CPUs working at 2.20 GHz.The impeller rotation has been modelled using the MRF approach and the simulations have been conducted in steady-state mode.Since experimental Eulerian information for cross-validation was available in the form of averaged PTV data, transient simulations were unnecessary.Given the fully baffled configuration of the vessel, the liquid free surface has been modelled as flat, imposing zero-stress and zero-flux conditions.Wall functions have been used with the standard k-e and RNG turbulence models.With the SST model, the flow equations have been solved in the boundary layer.This approach has been validated retrospectively by checking that the normalised wall unit was y þ < 1.This condition means that the grid was sufficiently fine to resolve the viscous sub-layer.
All turbulence model constants were set to their default values.The second order discretisation scheme was coupled with a leastsquare gradient scheme, a cell-based gradient limiter to avoid under-and over-shooting of the turbulent quantities and 3 nonorthogonality correctors.The solutions were considered numerically converged when the torque monitored at the impeller had become independent from further iterations and the sum of normalised residuals had dropped below 10 -5 , with most of them being in effect less than 10 -6 .

PTV velocity interpolation
CFD and PTV velocity data have been compared in terms of azimuthally averaged velocity at fixed positions in the vessel.This required the interpolation of the PTV Lagrangian velocity data on a 2D Eulerian grid.Each PTV velocity data point calculated is an average over many observations within a small voxel.The r-y vertical plane of the tank has been divided into 30 Â 60 ¼ 1; 800 voxels, each one identified by the indices i and j and having width dr ¼ 3 mm and height dy ¼ 3 mm.If the index k indicates the kth velocity event in a given voxel, the ensemble average over all n velocity events, i.e.
is slightly biased in favour of low velocity values.That is because (i) detectability and traceability of the PTV tracers are negatively correlated with their velocity and (ii) the residence time of a particle in the voxel is inversely proportional to its velocity, hence low-velocity particles are observed many times in the voxel while high-velocity particles are seen comparatively few times.For this reason, an alternative method has been used, based on the following steps.The velocity events observed in the voxel were grouped in individual trajectories and a first average was calculated for each trajectory (4).Then, the final value was obtained by averaging over all the trajectories entering that voxel (5).
The index l identifies a specific trajectory, mði; j; lÞ is the number of events for a specific trajectory in that voxel, and nði; jÞ is the number of trajectories entering that voxel.This method is conceptually similar to a weighted average based on residence time, as proposed by Chiti et al. (2011) for single-tracer PEPT data in stirred tanks and it mitigates the bias in the interpolated velocity.At low Reynolds numbers (Re 310), the differences between the ''trajectory-based" and the ''ensemble average" method were negligible.At higher Reynolds numbers (Re ¼ 943 and Re ¼ 1; 086), the differences were appreciable in the impeller jet and in the two recirculation loops, reaching a maximum of $ 5 % of the tip speed ($16 % of the local fluid velocity).
2.4.Newtonian simulations 2.4.1.Mid transitional regime A first Newtonian case study has been defined to conduct a mesh independence test and to assess the effects of turbulence model and discretisation scheme under mid-transitional flow conditions (Re ¼ 943).The vessel has been discretised in three meshes of approximately 1.9, 4.2 and 6.8 million cells.In each one, a cylinder of diameter 2 T=3 and height T=3 has been defined around the impeller to apply the MRF method.For each mesh, approximately 40 % of the total number of cells was concentrated in this region to resolve the steep gradients.Mesh refinement was also present in the proximity of the tank walls and baffles.Maximum nonorthogonality was about 65 in all three meshes, and average non-orthogonality was about 6.
The flow has been simulated with the three grids using the standard k-e model, standard wall functions and a first-order upwind (UW) discretisation scheme.The solutions have been compared in terms of azimuthally averaged velocity profiles and power number.This has been calculated with two methods.The first involved the integration of the energy dissipation rate in the whole volume, that is a weighted sum based on the cell volumes V j : The second was based on the torque, K (N m), experienced by the moving parts, hence on the stresses acting on the surface of impeller and shaft: Ideally, the two equations should provide the same result.However, it is well known that RANS simulations tend to underestimate the energy dissipation rate in agitated vessels, even with fine computational grids (Bartels et al., 2000;Alexopoulos et al., 2002;Bartels et al., 2002;Aubin et al., 2004;Singh et al., 2011;Lane, 2017).Instead, the torque on solid walls is mainly determined by the wall model.Consequently, the values of power number obtained with the two methods do not match.The differences depend on the mesh size and particularly on the discretisation scheme (Coroneo et al., 2011).For this reason, an additional simulation has been conducted using the intermediate mesh and a second order upwind scheme (2UW).
In order to assess the effects of turbulence modelling, the flow has also been simulated on the intermediate mesh with the k-e RNG and k-x SST models.Some statistics of the mesh wall unit obtained with the SST simulation are given in Table 2 to demonstrate the viability of the boundary layer resolution.

Low transitional regime
In principle, one turbulence model may perform better than others only within a certain range of Reynolds number.Moreover, stirred tank flows at moderate Reynolds numbers are typically simulated using a laminar solver, even if the overall flow regime is transitional and some turbulence is produced by the impeller (Zalc et al., 2001;Kelly and Gigas, 2003;Ameur and Bouzit, 2012;Patel et al., 2015;Mendoza et al., 2018).However, the predictions worsen rapidly as the Reynolds number increases.For this reason, the laminar, standard k-e, RNG k-e, and k-x SST models have been applied to the flow of a glycerol-water solution (80 %-20 % wt.) in the low-transitional regime (Re ¼ 86), to address how RANS perform at low Re.The simulations have been conducted on the intermediate mesh of 4.2 million cells and the 2UW scheme has been used.The velocity data have been compared and validated against PTV measurements.The turbulent quantities obtained through RANS must be interpreted carefully.In fact, the flow could be locally laminar in some zones of the vessel.In the laminar regions, turbulent fluctuations are not present, while RANS models would predict very low values of k and e, which should not be trusted.Instead, it is possible to verify the predictions in those regions where turbulence occurs.In fact, the integral of e should return the expected power consumption of the mixer geometry.

Non-Newtonian simulations
In light of the findings obtained with the two Newtonian case studies, the non-Newtonian simulations have been conducted using the standard k-e model and the 2UW scheme on the mesh of 4.2 million cells.The Carbopol solutions formed caverns of moving fluid around the impeller, while the fluid outside remained static due to the fluid yield stress.Cavern formation has detrimental effects on fluid mixing.The phenomenon has been investigated in several studies and some models have been proposed to predict the cavern shape and size (Solomon et al., 1981;Elson et al., 1986;Amanullah et al., 1998;Wilkens et al., 2005;Hui et al., 2009).The cavern boundary is commonly determined as the isosurface (or isoline, with 2D data) were the fluid velocity equals 1 % of the tip speed.In the present paper, the caverns have been identified based on the azimuthally averaged 2D velocity.

Test of the correlation between shear rate and Lagrangian acceleration
The Lagrangian acceleration is a measure of the force acting on an infinitesimal element of fluid which moves with the local and instantaneous flow.It is responsible for stretching, contraction and curvature of the fluid element.The statistics of Lagrangian acceleration are of great importance in turbulence research, as they are related to the intermittent bursts of vorticity and strain that characterize turbulent flows (Zeff et al., 2003).Large values of acceleration and energy dissipation are strongly correlated (Reynolds et al., 2005).
A unique feature of Lagrangian measurement techniques, including PTV, is that the Lagrangian acceleration along the tracer trajectories can be easily obtained through finite differentiation of the coordinates with respect to time.Conversely, the calculation of Eulerian derivatives, such as shear rate, is quite difficult.CFD can be used to fill this gap, as the shear rate can be directly calculated in each cell as.
where S is the symmetric part of the velocity gradient.impeller region was highly correlated with the average shear rate over a range of Reynolds numbers between 73 and 1,257.Since stirred tanks operated at moderate Reynolds numbers are often designed and scaled-up based on the effective shear rate in the impeller region, this finding highlighted the importance of the Lagrangian acceleration measurements that can be obtained through PTV.
In the present paper, the simultaneous availability of acceleration and shear rate gives the capability to test whether the local values of the two quantities are also spatially correlated.The data from the two Newtonian case studies have been used for this analysis, as they provided the best agreement between simulations and experiments.The shear rate and the acceleration data have been interpolated on a 2D Eulerian grid of 1,800 voxels (the same used for the PTV velocity data) to obtain azimuthally averaged values.Then, since the two quantities ranged within several orders of magnitude and a power law relationship was expected, the Pearson correlation coefficient (PCC) between Log 10 h _ ci and Log 10 hai has been calculated.The correlation was tested in the whole vessel, based on all the 1,800 voxels, and in a cylinder around the impeller, defined as.
which comprised of 200 voxels.

Mesh independence
The comparison of the azimuthally averaged velocity profiles obtained with the three numerical grids showed no significant differences; radial velocity is displayed in Fig. 1.The use of a second order discretisation scheme instead of a first order led to small differences, mainly localised in the impeller discharge flow.The jet was slightly wider and the maximum velocity slightly lower in the case of the 2UW simulation.
The power numbers are listed in Table 3.With the k-e model and the first order UW scheme, the torque on the moving parts was essentially mesh independent, but the volume averaged energy dissipation rate was significantly underestimated and kept increasing slightly as the mesh was refined.However, such differences were much smaller than those due to the discretisation scheme.With the 2UW scheme, the power number obtained from torque was very close to the expected value of about 3.5 for a Rushton turbine operated at Re % 10 3 (Distelhoff et al., 1995), while the resolved energy dissipation was about 80 % of the expected value.This level of resolution was in line with the figures reported by other authors (Coroneo et al., 2011;Lane, 2017) and has been considered satisfactory.For this reason, all the simulations in the remainder of the paper have been run on the intermediate mesh, coupled with a 2UW discretisation scheme.

Effects of the turbulence model
The choice of turbulence model had very little effect on the azimuthally averaged velocity fields, with the profiles matching quite well with the PTV experimental data (Fig. 2).Appreciable differences were only present at r ¼ R=3 at the impeller tip, and in the lower recirculation loop, particularly for the tangential component.PTV measurements underestimated the velocity very close to the impeller (up to À23%), likely because the PTV tracer traceability decreases with their velocity.
Good agreement was also found in the 2D maps of azimuthally averaged velocity magnitude and flow streamlines.Note that even the smallest differences in the Eulerian velocity components can lead to significant discrepancies in the resulting streamlines.Fig. 3 shows the case of the k-e simulation.The other two models  gave very similar results.The shadowed area with white dashed border represents the projection of the baffles.In particular, the size and position of the two main recirculation loops have been predicted well.Some differences can be observed upper region of the vessel, and particularly in the flow separation region close to the wall at y % 0:8 H.Here the fluid flows through an adverse pressure gradient and a third circulation loop is formed.Romano et al. (2021b) showed that there was a large population of Lagrangian trajectories pointing downwards in the area downstream of the baffles at y % 0:8 H.The fluid would then meet the upper recirculation loop, gaining tangential velocity.However, due to the low mean velocity (< 5% U tip ), its relative variance was large.The result is that, in the upper part of the vessel, the experimental streamlines appear much more chaotic than the simulated ones.Flow sep-aration was also present at the edges between the wall and the bottom, as captured by both PTV and CFD.Below the impeller, the experimental velocity magnitude was much higher than the simulated counterpart.Most of the difference was due to the tangential component.
The choice of turbulence model had appreciable effects on the spatial distributions of the turbulent quantities, namely the energy dissipation rate, e (m 2 s À3 ), and the turbulent kinetic energy, k (m 2 s À2 ).All models predicted e % 1 D 2 N 3 in the impeller discharge jet, e % 10 À2 D 2 N 3 in the bulk, and e % 10 2 D 2 N 3 around the impeller blades (Fig. 4).Moderate values of energy dissipation were also present around the baffles and close to the wall at y % H=3, due to the impeller jet impinging these surfaces.These orders of magnitude agree with those reported in the literature, although most published works concerning the energy dissipation in stirred tanks focus on fully turbulent conditions.Soos et al. (2013) conducted LES of a flow driven by a Rushton turbine at Re ¼ 12; 000 and reported phase-resolved profiles of e spanning between 10 D 2 N 3 and 130 D 2 N 3 .It is well-known that high values of turbulent kinetic energy and dissipation can be found in the trailing vortices behind the impeller blades (Derksen and Van den Akker, 1999;Sharp and Adrian, 2001;Escudié et al., 2004;Delafosse et al., 2009).Fig. 5 shows the isosurfaces corresponding to e ¼ 2 D 2 N 3 for the standard k-e model and e ¼ 1 D 2 N 3 for the RNG and SST models, highlighting the shape of the trailing vortices.A higher threshold was necessary for the k-e simulation because of the higher dissipation predicted in the whole impeller discharge flow.The standard k-e model provided also the highest integral of energy dissipation in the whole tank.This was directly reflected in the predicted power number, given in Table 3. Being much closer to the expected value, the standard k-e model outperformed the RNG and SST models.
The spatial distributions of the turbulent kinetic energy (Fig. 6) were qualitatively similar to those of the dissipation rate, with high values close to the impeller blades (k % 10 À1 U tip 2 ), moderate values in the discharge jet (k % 10 À2 U tip 2 ) and low values in the bulk (k % 10 À3 U tip 2 ).There were some differences between the three turbulence models.The standard k-e model provided the highest values in the discharge region (k ¼ 0:037 U tip 2 at y ¼ H=3 and r ¼ R=2).In the SST simulation, the velocity gradients and the turbulent kinetic energy close to the wall have been determined through direct integration of the flow equations within the boundary layer.Values as low as k % 10 À5 U tip 2 can be observed in the regions of flow separation close to the peripheral wall at y % 0:8 H and to the vessel bottom.Although the k-e and RNG models also predicted an adverse pressure gradient in those regions, this was not reflected strongly in the levels of turbulent kinetic energy.This may be ascribed to the use of wall functions.

Newtonian simulation, low transitional regime
At Re ¼ 86, the azimuthally averaged velocity profiles obtained with the three turbulence models were virtually identical and agreed very well with the PTV measurements (Fig. 7).The laminar simulation also provided extremely similar results.That was because the ratio of the turbulent viscosity over the fluid viscosity was small (l T =l % 10 À2 in the impeller jet and much lower in the bulk).This suggests that the flow was essentially laminar in most parts of the tank.Transitional flows in agitated vessels are notoriously Reynolds number dependent and, as expected, all three normalised velocity components were different from those at Re ¼ 943.
The distributions of normalised energy dissipation rate and turbulent kinetic energy were also very different from those in the mid transitional case study.It is well-known that transitional flows may show turbulent behaviour close to the impeller and laminar behaviour in the bulk (Gabelle et al., 2013).The proportion between the turbulent and laminar regions changes with the Rey-     ).Close to the baffles and peripheral wall, due to the jet impinging on the solid surfaces, e % 1 D 2 N 3 and k % 10 À3 U tip 2 .The turbulence levels in the jet were smaller than in the previous case (e % 10 À3 D 2 N 3 , k % 10 À4 U tip 2 ).As previously explained, in a purely laminar flow turbulent fluctuations are not present, and thus k and e cannot be defined.Therefore, the very small values of k and e obtained in the bulk flow are likely to be unphysical and should be ignored.At Re ¼ 86, the expected power number for this geometry is Po % 3:5 (Distelhoff et al., 1995).The simulations resolved most of the energy dissipated in the vessel, while the torque on the moving parts was extremely close to the expected value (Table 4).Again, the k-e model gave the best predictions, although the differences with the other models were small.

Non-Newtonian simulations
With the two power-law fluids, the velocity magnitude predicted by the simulations was higher than the experimental values in every location of the vessel, particularly in the impeller jet.Although the differences in the individual velocity components were not excessively large, they had dramatic effects on the resulting mean flow patterns (Fig. 9).In particular, the position, shape and size of the two main recirculation loops were quite different in both cases.The causes for such discrepancies are not clear and need further investigation.
The maps of the normalised, azimuthally averaged energy dissipation rate and turbulent kinetic energy of the flow at Re ¼ 1; 086   are illustrated in Fig. 10.A direct comparison with results obtained for the Newtonian fluid agitated at Re ¼ 943 (Figs. 4 and 6) revealed a strong similarity between the two flows.Although the two Reynolds numbers are very similar, it seems unlikely that the fluid rheology had so little impact on the distribution of the turbulent quantities.However, the power-law exponent of the CMC solution was n ¼ 0:7, not much smaller than unity, and the fluid apparent viscosity at the Metzner-Otto shear rate was very close to that for the Newtonian glycerol solution (4.4 Â 10 À2 Pa s and 4.6 Â 10 À2 Pa s, respectively).In order to isolate the role of fluid rheology on the turbulence distribution in the vessel, future CFD simulations should involve non-Newtonian fluids with a much smaller value of power law index.The simulations of the two yield stress fluids were in reasonable agreement with the experimental data.Compared to PTV, CFD slightly overestimated the cavern extension (Fig. 11).In this figure, hUi indicates the azimuthally averaged magnitude of the 2D velocity, i.e. the combination of radial and axial components, as obtained through CFD, hUi CFD , and through PTV, hUi PTV .Roughly, CFD predicted U CFD ¼ 2% U tip where U PTV ¼ 1% U tip was measured.
In the literature, much larger discrepancies have been reported.For instance, Xiao et al. (2014) compared six CFD turbulence models in the range 7 Re 163 and found that the k-e model overpredicted the cavern size quite significantly, particularly at Re 70.
In conclusion, non-Newtonian rheology seems to amplify the difficulties associated to the RANS of stirred tanks operated in the transitional regime.

Correlation between shear rate and Lagrangian acceleration
The analysis of the data interpolated on the vertical plane of the vessel revealed a strong correlation between the shear rate, obtained from CFD, and the Lagrangian acceleration, measured with PTV (Fig. 12).In the mid transitional regime (Re ¼ 943), the Pearson's coefficient of correlation between the two quantities was PCC ¼ 0:657.The underlying scaling rule was a power law fitting of the form hai=a tip ¼

Conclusions
The transitional flows of Newtonian and non-Newtonian fluids in a lab-scale vessel stirred with a Rushton turbine have been simulated through steady RANS.Such simulations are very common, despite their limitations, because of their computational affordability.The numerical results have been validated against PTV experiments.To the best of the authors' knowledge, this is the first scientific contribution presenting a comparison between Finite Volume Method (FVM) and PTV data in stirred tanks.The main findings are summarised below.Recently, a paper discussing the comparison between Lattice-Boltzmann simulations and PTV measurements has been published by Hofmann et al. (2022).
The mean flow of a Newtonian fluid agitated at Re ¼ 943 was independent from the use of a first or second order discretisation scheme.However, realistic estimations of the turbulent quantities could only be obtained with a 2UW scheme.This has been used for all subsequent simulations.The turbulence model had little effects on the mean flow.The velocity profiles predicted by the k-e, RNG and SST models matched well the PTV measurements.Appreciable differences were only present close to the impeller tip.The cause was likely the low traceability of high velocity PTV tracers close to the impeller.The choice of the turbulence model slightly affected the distributions of the turbulent quantities.The standard k-e model provided the power number closest to the expected value.About 80 % of the total energy dissipation was resolved.
With the Newtonian fluid at Re ¼ 86, the normalised mean velocity, energy dissipation and turbulent kinetic energy were very different from the previous case at Re ¼ 943.The turbulence offset in the bulk flow was significantly lower.Some authors have reported superior performance of the RNG and SST models over the standard k-e in simulating stirred tank flows.In this work, the opposite has been observed both in the mid and in the low transitional flow regime.This might be a consequence to the MRF approach to the impeller rotation.With shear thinning fluids agitated at Re ¼ 310 and Re ¼ 1; 086, CFD predicted a higher mean velocity in most of the tank, compared to PTV.The flow patterns and the shape of the recirculation loops did not correspond to the measured ones.The specific reason for such discrepancies is unclear.
In the case of Herschel-Bulkley fluids, the cavern boundaries predicted by the k-e model at Re ¼ 133 and Re ¼ 220 were in reasonable agreement with those obtained through PTV measurements.
Based on the Newtonian fluid data at Re ¼ 943 and Re ¼ 86, the measured Lagrangian acceleration and simulated shear rate were correlated in space.This is a further confirm that Lagrangian acceleration measurements obtained by means of PTV can be used to infer the local shear rate intensity.The underlying scaling rule was a power law of the form hai $ hc Romano et al. (2021b) have shown that, for different fluids (both Newtonian and non-Newtonian), the average Lagrangian acceleration in the
M.G.Romano, F. Alberini, L. Liu et al.Chemical Engineering Science 267 (2023) 118294 nolds number and the fluid properties.Here, turbulence was only present in the impeller jet and in those regions were the flow interacted with solid walls.For brevity, only the k-e solution is displayed in Fig. 8. High values of e and k were only present in a small region around the impeller blades (e % 10 2 D 2 N 3 , k % 10 À1 U tip 2

.
The exponent b was 0.66 and 0.88 in the two cases.More flows should be analysed to assess

Fig. 11 .
Fig. 11.Boundaries of the caverns formed by the two yield stress fluids, as simulated with CFD and measured with PTV, at a) Re ¼ 133 and b) Re ¼ 220.

Table 1
Rheological properties of the fluids and flow conditions.

Table 2
Average and maximum y þ obtained with the SST simulation on a mesh of 4.2 million cells (Newtonian fluid, Re ¼ 943).

Table 3
Power number in function of the mesh refinement, discretisation scheme and turbulence model (Newtonian fluid, Re ¼ 943).

Table 4
Power numbers predicted from different turbulent models (Newtonian fluid, Re ¼ 86).