Unsteady forcing of turbulence by a randomly actuated impeller array

We investigate the unsteady forcing of turbulent flow in a well-stirred reactor using opposing arrays of pitched-blade impellers which randomly and independently reverse rotation. We systematically explore the dependence of the large-scale motions and the homogeneity and isotropy of the turbulence upon the forcing. We identify three dimensionless control parameters: the source fraction (the fraction of time spent in clockwise motion), the dimensionless forcing period and an impeller Reynolds number. We find the timescale of unsteady motion corresponds to the forcing period T, the average period of impeller reversal, independently of the impeller angular speed Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega$$\end{document} and source fraction. As in jet-stirred tanks, unsteady forcing substantially increases the unsteady kinetic energy, energy dissipation, integral length scale and Taylor microscale Reynolds number (Rλ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_\lambda$$\end{document}) and improves the homogeneity and isotropy of the flow, provided the source fraction is chosen optimally and the forcing period is sufficiently large (ΩT>103\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega T > 10^3$$\end{document}); impeller Reynolds number has a relatively small influence. The forcing period must be matched to angular speed: decreasing the forcing period below this threshold results in a less intense, more inhomogeneous turbulent flow. Spectra of two-point velocity increments demonstrate that unsteady energy injection is dominated by axial shear generated across impellers and becomes less prominent at smaller scales. However, even at Rλ≈354\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_\lambda \approx 354$$\end{document}, the signature of this unsteady forcing can still be detected in near-dissipation-range statistics. These observations provide insight into optimisation of forcing and the mechanism of energy transfer when using unsteady forcing to generate turbulence in confined vessels.


Introduction
For laboratory study of turbulent flows, it is often desirable to produce homogeneous, isotropic turbulence within a confined vessel. Multiple designs of varying complexity are available for this task. These include single-impeller stirred tanks (Levins and Glastonbury 1972;Sheng et al. 2000;Yoon et al. 2005;Steiros 2017), oscillating grids (Srdic et al. 1996;Mann et al. 1999), counter-rotating impellers (von-Kárman mixers) (Voth et al. 2002;Kuzzay et al. 2015;Worth 2010), fan or jet-stirred "bomb" type reactors (Hwang and Eaton 2004;De Jong et al. 2009;Chang et al. 2012;Dou et al. 2016;Bounoua et al. 2018) and randomly actuated jet arrays (Variano and Cowen 2008;Bellani and Variano 2014;Carter et al. 2016;Pérez-Alvarado et al. 2016;Esteban et al. 2019). The motivation behind these designs has been to generate turbulent flows which are homogeneous and isotropic within some, preferably large, fraction of the vessel. Additionally, it may be useful to: produce a relatively weak mean flow to facilitate Lagrangian measurements (Mann et al. 1999;Voth et al. 2002;Zimmermann et al. 2010); attain large Reynolds numbers whilst being able to measure dissipative scale motions (Worth and Nickels 2011;Lawson and Dawson 2015); or avoid damage of fragile, suspended particles (Myerson 2002;Hu et al. 2011).
For turbulent reacting flows, these devices constitute the canonical "well-stirred reactor". This idealisation allows the flow to be described by a handful of parameters and allows aggregate measurements of the reaction rate (e.g. by monitoring reaction products) to be interpreted as representative of the reaction rate throughout the entire vessel. However, in practice, the reaction may occur inhomogeneously, especially in regions of high shear (Myerson 2002). This issue means that the geometry of the agitation vessel must be accounted for in empirical correlations for particle-fluid mass transfer (Levins and Glastonbury 1972). Additionally, suspended particles may be fragile and become damaged under high shear, which increases with impeller size. This is often the case in bioreactors and crystallisation within stirred vessels (Myerson 2002;Hu et al. 2011), which leads to a trade-off between impeller size and mixing efficiency. These examples provide practical motivation for reactor designs which produce better flow homogeneity and mixing with smaller impellers.
Recently, designs based around arrays of randomly actuated jets have been shown to generate nearly homogeneous, isotropic turbulence over a substantial fraction of the measurement vessel (Variano and Cowen 2008;Bellani and Variano 2014;Carter et al. 2016;Pérez-Alvarado et al. 2016;Esteban et al. 2019). The turbulence intensity and mixing is optimised when only a small fraction ( ∼ 12.5% ) of the jets are actuated at any instant (Variano and Cowen 2008). By random, unsteady actuation of the jets, energy can be directly injected into large scale motions which promote mixing and result in improved homogeneity and higher Reynolds numbers. A limitation of these designs is that the jet velocity is fixed, which means the turbulence intensity cannot be easily varied in steady state operation. Similarly, unsteady forcing by impellers in von-Kárman and simply-stirred tanks can also achieve improved mixing (Roy and Acharya 2012;Woziwodzki 2011;Steiros 2017). In such designs, unsteadiness is introduced by changing the speed or direction of the impellers and turbulence intensity can be varied by changing the impeller speed.
Unsteady agitation by impellers, rather than intermittently actuated jets, differs in several important aspects. Unlike jet driven flows where the power input must be switched off to generate unsteadiness, impellers can provide continuous power input whilst reversing direction to create unsteadiness. The mechanism of energy injection differs, since both the linear and angular momentum imparted by the impellers varies in time. Moreover, in the jet driven designs discussed above, the jet velocity is fixed. In comparison, impeller driven flows naturally possess an additional degree of freedom, since the impeller speed can be varied. Therefore, the optimal forcing parameters of jet stirred flows (e.g. Variano and Cowen (2008); Pérez-Alvarado et al. (2016)) may differ in comparison to impeller stirred flows. Whilst random stirring algorithms have been partially explored in fan-stirred mixers (Zimmermann et al. 2010), a systematic investigation of the parameter space of unsteady forcing has not been conducted for impeller driven flows to date.
In this paper, we describe a new impeller stirred mixer configuration, which combines the advantages of large-scale mixing created by unsteady forcing with the planar symmetry of opposed jet array designs. The design consists of two opposing arrays of sixteen, independently driven impellers within a transparent, cuboidal tank. Using many (32) small impellers, rather than few large ones, we increase the number of degrees-of-freedom for the forcing and allow for lower tip velocities (and therefore lower shear near impellers) for a given power input. Unsteadiness is introduced by randomly reversing the direction of the impellers, whilst energy input and Reynolds number are controlled by impeller speed. We investigate the effect of three forcing parameters upon the mixing, homogeneity and isotropy of the turbulence. These are: the source fraction (the fraction of time spent in clockwise rotation), the characteristic timescale of the unsteady forcing and the speed of the impellers.
The paper is structured as follows. In §2, we describe: the opposed impeller array ( §2.1), the forcing parameter space studied ( §2.2), the PIV measurements performed ( §2.3) and their subsequent post-processing in to derive turbulence statistics ( §2.4). Our results are presented in §3. We begin by introducing the mean flow field in §3.1, then examine the timescale of the unsteady motions created by the forcing in §3.2. We then examine the optimisation of the forcing to maximise large scale energy injection, homogeneity and isotropy in §3.3 and §3.4, respectively. We document the influence of the forcing at different length scales of the flow in §3.5. We summarise our work in §4. A summary of nomenclature used throughout the paper is provided in appendix A.

Experimental apparatus
The mixing tank and impeller array is shown in Fig. 1. It consists of a 45 L capacity transparent PMMA tank (internal dimensions W × H × H = 500 × 300 × 300mm ) with removable lid and N I = 32 independently driven pitchedblade impellers. These are spaced 75 mm apart on a 4 × 4 square grid and are mounted through the side-wall of the tank. The distance between the arrays is 428 mm . The tank is constructed from stainless steel, PMMA and PLA plastics in order to be non-reactive to strong acids and bases used in reacting flow experiments. The lid of the tank is removable; the tank is filled to a level slightly above the lid and careful insertion allows trapped air to be avoided. The working fluid (de-ionised water) can be recirculated through a builtin particle filtration and de-ionisation loop or drained into a storage tank below for retention. The water temperature was not regulated and varied between 15 − 18 • C . From this we estimate the mass density is 999 kg∕m 3 and kinematic viscosity varies between = 1.05 − 1.14 mm 2 ∕s.
The detail design of the pitched blade impellers is shown inset in Fig. 1. The impellers are 3D printed from PLA plastic and have six blades at a pitch of 45 • , thickness of 2 mm , a chord length of 14 mm and radius R = 35 mm . In preliminary tests, we trialled impellers with 90 • pitch, but found these to generate unsatisfactory low levels of turbulence. All impellers have the same handedness. A refinement of the design might consider alternating the handedness of impellers, in order to ensure no preference for handedness (e.g. through injection of helical motions). The impellers are directly driven by NEMA17 stepper motors (Ooznest, 1702HS133A) and are independently controlled using eight off-the-shelf Arduino-based stepper motor "CNC-Shields" each holding four A4988 stepper motor drivers. This allows for inexpensive, versatile and precise control of the forcing speed, set by the step rate of each motor. An important detail of the mechanical design is the sealing of the drive shafts which pass through the side walls. These are sealed with two nitrile rubber rotary shaft seals on either side of the side wall; a single shaft seal was insufficient to prevent leaking. Additionally, clearance between the motors and outer seal is designed to avoid fluid leaking into the motor body. The agitation speed is limited by the stall of the motors at high step rates largely due to friction from the shaft seals.

Forcing
Like opposed jet arrays, the opposed impeller array provides a large parameter space to explore in terms of forcing. To simplify our investigation, we implement the "sunbathing" algorithm of Variano and Cowen (2008), which has been widely used in jet stirred facilities (Variano and Cowen 2008;Bellani and Variano 2014;Pérez-Alvarado et al. 2016;Carter et al. 2016). We drive the impellers at a fixed angular speed = 2 f I whilst randomly reversing the direction of motion. Each impeller switches between spinning clockwise, thrusting fluid towards the center of the tank for a time F , and then anti-clockwise, thrusting fluid towards the walls for a time R . We call the clockwise motion forward thrust and the anti-clockwise motion reverse thrust. Following Variano and Cowen (2008), the duration of forward F ∼ N( F , F ∕3) and reverse thrust R ∼ N( R , R ∕3) To describe this parameter space in dimensionless terms, we find it convenient to parametrise the forcing in terms of a forcing period T = F + R and a source fraction = F ∕T . The forcing period T represents the average cycle length of an individual impeller and is twice the average interval between impeller reversals. The source fraction represents the average fraction of time spent in forward thrust. Thus, for = 0 , the impellers are always in reverse thrust; for = 1 , the impellers are always in forward thrust; and at intermediate values the forcing is unsteady. In dimensionless terms, the forcing is characterised by the source fraction and dimensionless forcing period T . The geometry of the impellers and kinematic viscosity of the working fluid form an additional dimensionless parameter. We write this as an impeller Reynolds number Re I ≡ 2 R∕ based on the tip velocity of the impellers and is analogous to the jet Reynolds number in opposed jet arrays. Therefore, in contrast to opposed jet arrays where the exit velocity of the jet is fixed, impeller arrays have an extra degree of freedom: the dimensionless forcing period T.
This parametrisation takes the forcing period T to be the characteristic timescale of the forcing. There are other choices of characteristic timescale: the on-time F and the "temporary state" time s = T∕(2N I ) (i.e. the average duration over which a specific pattern of impeller motion is maintained) have also been suggested (Variano and Cowen 2008). We investigate the appropriate choice of characteristic timescale in §3.5.3. To non-dimensionalise flow field quantities, we choose the impeller radius R and angular speed to form characteristic length and timescales. This is because fluctuating velocities and dissipation rates in impeller driven flows are widely observed to scale with these two length and timescales at high Re I (see e.g. Voth et al. 2002;Yoon et al. 2005;Zimmermann et al. 2010). We confirm this scaling in §3.3.

PIV measurement
We conducted several PIV measurements to systematically explore this parameter space. The forcing parameters for each set of measurements are provided in Table 1. These are organised into three series: Series , in which the source fraction is varied and the the impeller speed f I = 4 Hz and forward time F = 8 s are fixed; Series , in which the impeller speed is varied and the forcing period T = 32 s and source fraction = 0.25 are fixed; and Series T, in which the forcing period is varied and the impeller speed f I = 4 Hz and source fraction = 0.25 are fixed. The effect of changing each dimensionless forcing parameter can therefore be isolated by comparing results across series. For example, data from Series and T can be combined to compare the flow at matched T but varying Re I (i.e. impeller speed). The flow at the mid-plane of the tank was measured using 2D PIV for each forcing condition. The coordinate system and field of view is illustrated in Fig. 1b. The field of view covers most of the tank, measuring 456 × 304 mm , but vectors can only be reliably obtained for a region −0.67 < x∕H < 0.67 and 0.48 < y∕H < 0.48 away from the impellers and walls, where H = 300 mm is the height of the tank.
For each case, 2000 image pairs were recorded at a sample rate of 0.5 Hz using an IMPERX Bobcat IGV-B6620 CCD camera equipped at a resolution of 6600 × 4400 pixels. The field of view was imaged using a Sigma 105 mm Macro lens at an f-number f # = 11 . Double-pulse laser illumination was provided by a Litron Bernoulli PIV laser 1 with a nominal pulse energy of 200 mJ per pulse at a wavelength of 532nm. 1 The beam was expanded into a laser sheet with a −10 mm focal length cylindrical lens and brought to a narrow waist near the centre of the tank with a 1000 mm focal length plano-convex spherical lens. Based on the reported beam parameter product of the laser, we estimate the laser sheet thickness to be 3.4 mm . We seed the flow with 56 m median diameter Nylon-12 seeding particles (Vestosint 2157) with mass density 1016 kg∕m 3 to achieve a seeding concentration of approximately 6 × 10 −3 particles per pixel. Based on our estimates of the Kolmogorov timescale , the particle Stokes number St does not exceed 3.3 × 10 −2 (Toschi and Bodenschatz 2009). The PIV time separation was chosen to provide a root-mean-square particle displacement at the centre of the tank of approximately 3 px , which represents a compromise between ensuring measurable displacements near the impellers (where velocity fluctuations are significantly larger) and minimising the influence of measurement noise for particles near the slower-moving centre.
PIV interrogation was performed using an in-house 2D PIV code. This is derived from a 3D PIV code developed by the author (Lawson 2016) and implements standard features such as multi-pass window deformation with (Gaussian) window weighting, Gaussian peak fitting and Lanczos image interpolation between steps (Raffel et al. 2018). The code has been benchmarked against 2D PIV interrogation using DaViS 8.22 (LaVision, GmbH) and found to produce comparable results in a fraction of the processing time. Images were preprocessed in three steps: levelisation by subtraction of the minimum background intensity level, intensity capping at two standard deviations in excess of the median intensity, and a low-pass Gaussian filter over a 3 × 3 pixel neighbourhood. Images were interrogated using a four-pass window deformation scheme, using Gaussianweighted square windows of 64 2 , 64 2 , 32 2 , 32 2 pixels, with a constant vector spacing x = 1.08 mm and window size W = 2.17 mm at the final pass. The interrogation window size at the final pass is provided in Kolmogorov units in Table 1. Between passes, invalid vectors are infilled using secondary or tertiary correlation peaks if these are spatially consistent or linearly interpolated from neighbours if these are inconsistent. Invalid vectors are detected using a normalised median test (Raffel et al. 2018).

Postprocessing
Data are subsequently postprocessed to obtain statistics of the mean flow field and turbulence. The Reynolds decomposition is applied to decompose the instantaneous flow field y, z] and the associated Cartesian unit vectors are e 1 , e 2 , e 3 . An overline u denotes a temporal average over time t whilst a superscript ′ denotes temporal fluctuations. Since only two velocity components u 1 and u 2 are measured, the kinetic energy of the mean K = 1 2 (U 2 1 + 2U 2 2 ) and fluctuating k(x) = 1 2 (u �2 1 + 2u �2 2 ) components is approximated by assuming large scale axisymmetry, based on the four-fold discrete rotational symmetry of the tank around y = z = 0 . We expect this surrogate provides a reasonable approximation of the fluctuating kinetic energy near the symmetry axis of the tank, but may break down towards the edges. Temporal energy spectra E ij (x) of single-point velocity fluctuations are obtained by decomposing each timeseries of 2000 samples into eight windowed subseries of 500 samples each using a rectangular window function. Subsequent spatial averages of these quantities are obtained over the large field of view and centre of the tank and are denoted with angled brackets ⟨⋅⟩ c and ⟨⋅⟩ v . These regions are specified in §3.1. Vectors flagged as invalid are rejected from statistics.
The second order structure function at point x and scale and orientation r = rr . This definition allows us to recover the isotropic scaling properties of the turbulence from the anisotropic and inhomogeneous mean flow using the SO(3) decomposition (Knutsen et al. 2020;Arad et al. 1999). Using this decomposition, the longitudinal component expressed as D LL = D ijrirj is averaged over orientation r = [cos , sin cos , sin sin ] using, 13 Page 6 of 16 in which we have invoked axisymmetry about y = z = 0 . This may be justified given that D LL (x, r) is evaluated at positions x near the geometric centre of the tank, where axial symmetry is expected to hold. A similar operation is performed to obtain longitudinal autocorrelation function ⟨R LL ⟩̂r(x, r) from the two-point correlation tensor given by: Jong et al. (2009), we estimate the dissipation rate ⟨ ⟩ c near the centre of the tank using the compensated second order structure function method as where C 2 = 2.13 is the Kolmogorov constant (Sreenivasan 1995). This procedure is demonstrated in appendix B. The Taylor microscale Reynolds number R reported in Table 1 is based on the dissipation rate ⟨ ⟩ c and turbulent kinetic energy ⟨k⟩ c measured at the centre of the tank.
To obtain the spatial distribution of the time average kinetic energy dissipation rate (x) , we use the large eddy PIV method (Sheng et al. 2000;De Jong et al. 2009). In this, the kinetic energy dissipation rate is estimated using the Smagorinsky sub-filter stress model based on the velocity field filtered at a length scale chosen to be within the inertial subrange. Following the procedure detailed in appendix B, the filter scale is chosen = 6 x = 6.50 mm and the value of the Smagorinsky-Lilley constant is tuned so that estimates of the average dissipation rate obtained from the large eddy PIV method are consistent with estimates of the dissipation rate near the centre of the tank obtained from the second order structure function.

Mean flow and velocity fluctuations
To provide an overview of the flow and to initiate our discussion, we introduce the mean flow and the spatial distribution of velocity fluctuations generated under steady and unsteady forcing conditions. As an example, Fig. 2 presents the mean flow field at the mid-plane of the tank for varying source fraction = 0, 0.25 and 1 but constant angular velocity and mean on time F = 8 s . The flow field and turbulent kinetic energy distribution at = 0.125, 0.5 and 0.75 is qualitatively similar to that at = 0.25 . In all cases, the mean flow contains of a pair of large-scale, counter-rotating vortices, resulting in straining flow near center of the tank. This configuration resembles the mean flow induced by Ekman pumping in von-Kárman type mixers (Kuzzay et al. 2015). It persists even at = 1 , where the thrust of each impeller imparts a net momentum flux towards the centre of the tank. Qualitatively, we observe that the strength of this mean flow is weakened as the forcing is randomised and the intensity of unsteady fluctuations increase in magnitude. For steady forcing, we observe pronounced inhomogeneity in the energy of velocity fluctuations k near the impellers, with localised "hotspots" near regions of stronger mean strain. For unsteady forcing, the homogeneity of these velocity fluctuations is improved.
We now examine the velocity fluctuations in more detail. As an example, we consider the flow generated by forcing with = 0.25 , T = 32 s and f I = 4 Hz . In Fig. 3a, we present axial profiles of mean-square velocity fluctuations across the tank over the region −H∕8 < y < H∕8 (marked in orange in Fig. 2), in addition to the spatial average over the y direction in this region. For both axial and radial velocity components, velocity fluctuations are larger near the impellers, reaching a minimum value near the centre of the tank, and deviate by less than around 10% of the spatial average in this region. The extent of this homogeneous region and the relative magnitudes of each component depends upon the details of the forcing, which we shall examine in §3.4. Similarly, we show radial profiles of the velocity fluctuations over the region −H∕8 < x < H∕8 in Fig. 3b. Towards the walls, the radial velocity fluctuations decay to satisfy the no-penetration condition u 2 = 0 , whereas axial velocity fluctuations exhibit a peak near the wall. However, the kinetic energy of velocity fluctuations is uniform over a larger radial extent. This suggests that large-scale radial eddying motions are impeded by the wall and are deflected axially to preserve continuity, whilst approximately conserving their kinetic energy.

Timescale of unsteady motion
The examples presented so far have shown that unsteady forcing is effective at introducing energetic, large-scale motions which increase the magnitude of velocity fluctuations and improve their spatial homogeneity. What is the timescale of these motions? To answer this question, we compute spatially averaged temporal energy spectra ⟨E kk (x, )⟩ v of the velocity fluctuations averaged over the region shown in blue in Fig. 2. We interpret this statistic to be representative of the average kinetic energy over the entire tank. In Fig. 4a, we show spectra at fixed impeller speed and mean on time F , but varying source fraction . This has the effect of varying the forcing period T. When the forcing is unsteady, we observe a pronounced spectral peak at angular frequency T ≈ 2 , which corresponds to the average period of direction reversal. When the forcing is steady, we do not observe this peak. This spectral peak is considerably lower frequency than that associated with the average duration of each temporary state 2 ∕ s or the impeller actuation time 2 ∕ F . A similar spectral peak is not easily observed in spectra reported for jet stirred tank of Variano and Cowen (2008), but may be hidden amongst the statistical uncertainty of the data analysed. We conclude that for impeller stirred flows the characteristic timescale of the energy containing, unsteady motions is the forcing period T and not the mean forward time F or temporary state time s . To test the hypothesis that the timescale of the unsteady motion is the forcing period, we present the same spatially averaged spectra over a range of impeller frequency and forcing period spanning one order of magnitude in Fig. 4b . Across this range, we observe the same pronounced spectral peak at T ≈ 2 . Above this characteristic frequency, there is a decay in the spectral density of kinetic energy. Thus, the unsteady forcing injects energy into the flow at frequency 2 ∕T , which break down into more rapid, smaller scale motions.

Optimisation of large scale energy injection
To examine how the kinetic energy associated with the mean and fluctuating contributions varies under different forcing parameters, we compute spatial averages of the mean ⟨K⟩ v and fluctuating ⟨k⟩ v contributions to the average kinetic energy and energy dissipation rate ⟨ ⟩ v of the flow over the region marked in blue in Fig. 2. Again, we interpret this statistic to be representative of the average kinetic energy over the entire tank. Based on a dimensional analysis of the power input by the impellers, we expect ⟨k⟩ v ∕ 2 R 2 and ⟨ ⟩ v ∕ 3 R 2 to be functions of the dimensionless forcing period T , source fraction and impeller Reynolds number Re I = R 2 ∕ . These are shown in Fig. 5a as a function of the average source fraction , for fixed impeller speed and on time F . We observe that the majority of the kinetic energy is contained within unsteady motions, and that the fluctuating and total kinetic energy are maximised for a source fraction = 0.25 when the on-time F is held constant. The additional agitation introduced by unsteady forcing is profound, with the energy contained within unsteady motions roughly 2.6 times larger at = 0.25 than = 0 (steady clockwise rotation, reverse thrust). Moreover, unsteady forcing decreases the strength of the mean flow; at = 0.25 , the kinetic energy of fluctuations exceeds the mean flow by an order of magnitude. Similarly, the spatial average kinetic energy dissipation rate ⟨ ⟩ v is also maximised at = 0.25 . This demonstrates that the coupling of mechanical energy between the impellers and the flow is more effective under unsteady forcing conditions. We examine how these spatial average quantities vary with the forcing period and impeller speed for a fixed source fraction = 0.25 in figs. 5b and 5c. When the impeller speed is held constant, but the forcing period T is varied (orange symbols), we observe a weakly increasing trend in the spatial average fluctuating kinetic energy ⟨k⟩ v , but no discernible trend in ⟨K⟩ v or ⟨ ⟩ v . This corresponds to fixing the impeller Reynolds number Re I = R∕ , but varying the dimensionless forcing period T . A plateau appears to be reached at around T ≈ 10 3 . Similarly, when the forcing period T is held constant but the impeller speed is varied (green symbols), there is also a weakly increasing trend in the dimensionless kinetic energy ⟨k⟩ v ∕( 2 R 2 ) and no discernible trend in the dimensionless energy dissipation rate ⟨ ⟩ v ∕( 3 R 2 ) . Thus, we are able to compare data at matched dimensionless forcing period T but varying impeller Reynolds number. There is reasonable collapse at matched dimensionless forcing period, demonstrating that Reynolds number effects are relatively weak. We note that, in dimensional terms, ⟨k⟩ v and ⟨ ⟩ v vary by over two and three orders of magnitude respectively over the range of tested. We also note that, in preliminary tests with impellers with 90 • pitch, we also observed the same ⟨K⟩ v ∼ 2 scaling for the kinetic energy, which indicates the observation is robust to changes in impeller geometry. We conclude that the large-scale power input is an inviscid process, the power input scales with 3 R 2 and is a weak function of the dimensionless forcing period T and impeller Reynolds number, provided these are sufficiently large.

Optimisation of large scale homogeneity and isotropy
We now consider the optimisation of the homogeneity and isotropy of the turbulence generated by unsteady forcing. A simple measure of the homogeneity is to compare averages of quantities near the centre of the tank, denoted ⟨.⟩ c , to the corresponding spatial average over the entire tank, denoted ⟨.⟩ v . Because our field of view is limited, we take the region outlined in blue in Fig. 2 to represent the behaviour over the entire tank, whereas the centre region is outlined in orange. For a perfectly homogeneous flow, we would expect the ratio ⟨.⟩ c ∕⟨.⟩ v of any given quantity to be one. We present such ratios of fluctuating kinetic energy and energy dissipation rate in Fig. 6 as a function of the dimensionless forcing period. For fixed source fraction = 0.25 , the data collapse remarkably well across forcing period and impeller speed. Furthermore, the fluctuating kinetic energy is more uniformly distributed than the energy dissipation. As the dimensionless forcing period is increased, the flow tends to become more homogeneous, up to a threshold T ≈ 10 3 beyond which the homogeneity is not further improved. This is consistent with the plateau in the optimisation of the overall kinetic energy in Fig. 5b. The effect of source fraction is more difficult to discern from our data, since in these experiments the dimensionless forcing period was not held constant and has a strong influence upon homogeneity. This parameter appears to have a smaller influence than the dimensionless forcing period, since data over varying 0.125 ≤ ≤ 0.75 largely collapse upon the data at = 0.25 . However, we note that the best homogeneity is achieved at a source fraction = 0.125 , which also corresponds to the largest dimensionless forcing period tested. To examine the homogeneity in more detail, we compare spatially averaged profiles of the magnitude of the fluctuating kinetic energy and energy dissipation relative to their value near the centre of the tank in Fig. 7a under forcing conditions with matched dimensionless forcing period T . The spatial average is obtained over |y| < H∕8 and represents the flow along the symmetry axis the tank. We observe excellent collapse across profiles of fluctuating kinetic energy and energy dissipation with matched dimensionless forcing period and source fraction. Furthermore, we confirm that the homogeneity is markedly improved by increasing the dimensionless forcing period and that the best homogeneity is attained at a source fraction = 0.125 . For this case, the fluctuating kinetic energy is uniform to within 10% of the centre value over a region measuring 69% of the spacing between arrays ( 297mm ), whilst the energy dissipation is uniform to within 10% over a region 49% of the array spacing ( ∼ 210 mm ). This homogeneous region covers a substantially larger fraction of the vessel than achieved by other jet stirred apparatuses to date: using the same metric, Carter et al. (2016) report a central homogeneous region spanning 27% of the distance between jet arrays. Similarly, Bellani and Variano (2014) find the central homogeneous region between spans 11.7% of the distance between their jet arrays.
In Fig. 8, we examine the ratio of the magnitude of the velocity averaged over the centre of the tank, as a function of both the source fraction and dimensionless forcing period. Large scale isotropy implies 1 ∕ 2 = 1 and deviations from unity mark increasingly pronounced anisotropy. We observe that, for all forcing conditions tested, the axial velocity fluctuations are always weaker than radial velocity fluctuations at the centre of the tank, with 1 ∕ 2 ≈ 0.84 at large T . This stands in contrast to jet stirred arrays, where the axial velocity fluctuations are typically larger than the radial fluctuations, e.g. Bellani and Variano (2014) report 1 ∕ 2 = 1.05 at the centre with forcing conditions optimised for homogeneity, whilst Carter et al. (2016) report 1 ∕ 2 = 1.28 − 1.88 over a wide range of forcing conditions. Notably, Carter et al. (2016) found the large scale motions became increasingly anisotropic as the forcing period is increased. In contrast, the dependence we observe upon the forcing parameters is quite weak. The source fraction appears to have the greatest influence upon the isotropy, with 1 ∕ 2 varying between 0.54 and 0.86 at = 0 and 0.125. The forcing period has a weaker influence; nonetheless, across varying impeller speeds and forcing period, the data collapse well upon the dimensionless forcing period T and the large-scale isotropy is slightly improved by increasing this quantity.

Influence of forcing across scales
Thus far, we have focused on the largest scales of motion. In this section, we examine the influence of forcing at different length scales.

Velocity increment statistics
To examine the energy of turbulent motions at different scales, we consider the second order structure function. Figure 9 shows the shell average second order structure function ⟨D LL ⟩̂r(x 0 , r) obtained for a point x 0 near the geometric centre of the tank, under steady and unsteady forcing conditions with varying source fraction but fixed forcing period T = 64 s and impeller speed f I = 4 Hz . Error bars show the range of variation over different orientations of the separation vector r . At larger scales r > 0.1H , the introduction of unsteady forcing increases the range over which the expected inertial range scaling D LL ∼ r 2∕3 1∕2 c evaluated at the centre of the tank, as a function of the forcing period T Fig. 9 Shell averaged second order longitudinal velocity structure function ⟨D LL ⟩̂r(x 0 , r) at varying source fraction but fixed forcing period T = 64 s and impeller speed f I = 4 Hz . Error bars show maximum and minimum variation over orientation . Dashed lines indicate expected small-and inertial-range scaling holds. There is relatively little change in the shell average at smaller scales: remarkably, the shell averaged second order structure function collapses at = 0, 0.125, 0.25 . However, the isotropy of the turbulence is improved across all scales: with steady forcing ( = 0, 1 ), the second order structure function exhibits large variations as the orientation of the separation vector r is varied; with unsteady forcing ( = 0.125, 0.25 ), these variations are markedly reduced. The anisotropy remains strongest at the largest scales of motion r ≈ H , where the influence of the boundary of the side walls becomes apparent.
To quantify the influence of the forcing period and impeller speed upon the isotropy of the turbulence across scale, we examine the ratio of longitudinal velocity increments ⟨D 11 (x, re 1 )⟩ c ∕⟨D 22 (x, re 2 )⟩ c obtained over a region near the centre of the tank. Data for Series-and T are plotted as a function of dimensionless forcing period in Fig. 10 at scales r∕ 0 = 1, 1∕4, 1∕16 where 0 = 69.3 mm . Across all scales and forcing conditions, axial velocity increments D 11 (x, re 1 ) tend to be smaller in magnitude than radial increments D 22 (x, re 2 ) , consistent with the behaviour of large scale velocity fluctuations seen in §3.4. When parametrised by the dimensionless forcing period T , the data exhibit good agreement at all scales. As the forcing period or impeller speed is increased, the isotropy of the turbulence is improved, with the largest improvements seen at the largest scales. Beyond T ≈ 10 3 , there is little improvement in the isotropy. This is consistent with our observations on the optimisation of large scale homogeneity in §3.4 and reinforces the point that the forcing period must be sufficiently large (compared to the impeller speed) to attain the best approximation of homogeneous and isotropic turbulence.

Integral length scale
To identify the effect of our forcing upon the length scale of the largest turbulent motions produced, we examine the dependence of the integral length scale of the turbulence upon the forcing in Fig. 11. This is obtained from the shell averaged longitudinal autocorrelation function ⟨R LL (x 0 , r)⟩̂r evaluated at a point x 0 near the centre of the tank using where r max corresponds to the extent of the field of view (0.96H) or the first zero crossing, whichever is smaller. This definition has the property that the integrand approaches 1 at r → 0 and is consistent with the definition of the integral length scale in homogeneous, isotropic turbulence (Pope 2000). The dashed lines show the integral length scale obtained for steady forcing, whereas markers show the integral length scale at various unsteady forcing conditions. We first remark that the introduction of unsteady forcing considerably increases the length scale over which velocity fluctuations are correlated within the tank. We observe that, for sufficiently large T ≥ 400 , the integral length scale becomes independent of the forcing period and source fraction and saturates at a value L int ∕H ≈ 0.21 . This is comparable to both the impeller diameter 0.233H and the spacing between impellers 0.25H. Saturation at this length scale is consistent with the notion that the large-scale motions are generated by the unsteady shear produced by adjacent impellers. Since the impeller spacing is fixed in our experiment, we are unable to examine whether the integral length scale would be further increased by increasing the impeller spacing. However, we note that in the jet-stirred turbulent flow of Carter et al.  Table 1 (2016), increasing the inter-jet spacing was not found to be effective in increasing the integral length scale.

Signature of forcing across length scales
To what extent does the signature of the large scale forcing present in the velocity increment statistics persist down scale? To answer this question, we examine (spatially averaged) spectra of two-point velocity increments u (x, re , t) , defined where is the temporal Fourier transform of the two-point velocity increment. These are shown in Fig. 12 for forcing at = 0.25 with f I = 4 Hz and T = 32 s for a range of scales r∕ 0 = 1, 1∕2, … 1∕32 , where 0 = 69.3 mm ≈ 1.07L int . At the largest scale r = , there is a clear spectral peak associated with the forcing at T ≈ 2 . This peak is largest in magnitude for axial velocity increments D 11 and is most prominent for increments measured in the radial direction, i.e. D 11 (re 2 , ) . In comparison, the peak is less pronounced in the radial component of velocity increments (i.e. D 22 ). We conclude that the forcing predominantly injects energy as large scale shear generated by the differential velocity between the jets produced by adjacent impellers. The signature of this forcing persists down scale and, for this case, can still be distinguished in spectra of transverse increments even at r = 0 ∕32 , which corresponds to 14 . A similar phenomenon can be observed in measurements of oscillating grid turbulence with artificial large scale intermittency, where fluctuating large scale energy injection has been shown to leave signatures in dissipation range statistics (Chien et al. 2013).

A comment on the effect of the forcing
One possible explanation for the improved homogeneity, isotropy and maximised kinetic energy associated with unsteady forcing at large T is to imagine the mean flow as a superposition of deterministic, but randomly chosen, quasi-steady temporary states, which average out to be statistically homogeneous and isotropic. This notion is applicable in the hypothetical limit of large T, where the temporary state time s = T∕2N I is also large with respect to the turbulent integral timescale L int ∕k 1∕2 . Two pieces of evidence suggest this is not the case here. Firstly, the integral timescale is the same order of magnitude as the temporary state duration (around ∼ 1 to 10 times the temporary state duration). Thus, each temporary state is far from equilibrium with respect to the turbulence and does not satisfy the quasi steady condition. Secondly, the energy spectra presented in figs. 4 and 12 demonstrate a roll off in spectral energy density beyond T ≈ 2 for the largest scale motions, indicating the initiation of a turbulent energy cascade at this timescale. This suggests the impeller reversal mechanism naturally provides a way to ensure the temporary states are far from equilibrium.

Conclusion
We have presented a novel design of well stirred reactor which uses opposing arrays of randomly actuated impellers to generate homogeneous, isotropic turbulence under confinement. This design, which uses off-the-shelf, inexpensive components has a natural advantage in comparison to previous approaches using randomly actuated jets: the impeller speed may be used to vary the intensity of the turbulence. In comparison to other fan-stirred mixer designs, this approach achieves significant improvement in the large scale homogeneity by introducing large scale motions through unsteady stirring, rather than increasing the impeller size, which may prove advantageous for the study of homogeneous reacting flows and shear-sensitive particles.
We have conducted a systematic characterisation of the large-scale turbulent flow generated under random forcing using the "sunbathing" algorithm (Variano and Cowen 2008). This characterisation may inform the optimisation of forcing in other fan-stirred mixers and forms a basis for future studies of turbulent, particle laden and reacting flows using this design. We find that unsteady forcing substantially increases the fluctuating kinetic energy, energy dissipation rate and integral length scale and decreases the strength of the mean flow. The forcing induces low-frequency energy containing motions, which subsequently decay turbulently. Their timescale grows with the forcing period T rather than the forward time F and the largest integral length scale obtained is comparable to the impeller spacing or diameter. Dimensional analysis reveals that three parameters are important: the dimensionless forcing period T , the source fraction and the impeller Reynolds number Re I .
We have characterised the effect of varying these parameters in terms of the overall turbulent kinetic energy, energy dissipation and the homogeneity and isotropy of the turbulence generated. Our most significant observation is that the timescale of the forcing must be appropriately matched to the impeller speed: in dimensionless terms, we require T > 10 3 to maximise the fluctuating kinetic energy, energy dissipation and homogeneity and isotropy of the turbulence across the vessel. This can be achieved either by increasing the forcing period T or impeller speed . For a fixed forward time F , the overall kinetic energy and energy dissipation are maximised by unsteady forcing at a source fraction ≈ 0.25 . However, over the parameter space tested, the best homogeneity of velocity fluctuations and kinetic energy dissipation is achieved for = 0.125 . The impeller Reynolds number Re I is found to have comparatively little influence.
We have examined the mechanism of energy injection by the forcing. Spectra of two point velocity increments suggest that the most energetic motions are due large scale shear, generated by the differential velocity between the jets of adjacent impellers. The saturation of the integral length scale of motion at a length scale comparable to the impeller spacing is consistent with this hypothesis, but not conclusive. The signature of this forcing is found to persist across scales. Even at a Taylor microscale Reynolds number R ≈ 354 , this signature can be observed in spectra of near inertial-and dissipation-range velocity increments.

Nomenclature
The nomenclature used within the manuscript is summarised for convenience in Table 2.

Measurement of kinetic energy dissipation rate
In this section, we illustrate the application of the methods outlined in §2.4 to obtain estimates of the kinetic energy dissipation rate. In Fig. 13a, we present the compensated, second order longitudinal structure function 1 r D 3∕2 LL obtained for Series-, after orientation and spatial averaging over the centre of the tank. This quantity is maximised in the inertial range and defines our estimate of the dissipation rate ⟨ ⟩ c near the centre of the tank according to (2). We use the large eddy PIV method to obtain the spatial distribution of the mean kinetic energy dissipation rate (x) (Sheng et al. 2000;De Jong et al. 2009). In this, the kinetic energy dissipation rate is approximated as (x) = − ijsij , where ij is the modelled sub-filter stress tensor and s ij = 1 2 ( ũ i ∕ x j + ũ j ∕ x i ) is the resolved strain tensor obtained by filtering the measured velocity field with a box filter at scale . Using the Smagorinsky model for the sub-filter stress, the dissipation rate is therefore obtained as where C s is the Smagorinsky-Lilley constant. Since only four of nine components of the resolved strain tensor are available, we approximate the tensor norm s ijsij as following De Jong et al. (2009), which has identical mean to s ijsij under isotropy assumptions. This method takes the filter scale to be within the inertial subrange.
Our estimate of the dissipation rate at the centre of the tank ⟨ ⟩ c obtained from the second order structure function (2) is subsequently used to tune the Smagorinsky-Lilley constant for large eddy PIV estimation of the dissipation field (x) using (6). Under this approximation, it is necessary to select a filter scale where the resultant estimate of the dissipation rate is independent of the (6) (x) = 2C 2 s 2 (s ijsij ) 3∕2 (7) s ijsij = 2s 2 11 + 2s 2 22 + 2s 11s22 + 6s 2 12 filter scale. Figure 13b demonstrates that a filter scale = 6 x = 6.50 mm , the dissipation estimate becomes scale-invariant. By subsequently equating our estimate of the dissipation rate at the centre of the tank ⟨ ⟩ c to the corresponding large-eddy PIV estimate, we obtain a tuned value of C 2 s = ⟨ ⟩ c ∕(2 2 � (s ijsij ) 3∕2 � c at the chosen filter scale. Across all our data, there is little variation in this quantity, with values of C s ranging between 0.104 and 0.117.
Author Contributions JML designed the study, built the apparatus, conducted experimental work and drafted the manuscript. BG provided supervision and edited the manuscript.
Funding This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 846648.
Availability of data and material pertinent data for this paper will be made available at doi: https:// doi. org/ 10. 5258/ SOTON/ D2056 upon publication.

Financial interests
The authors declare that they have no financial interests.

Non-financial interests BG is an associate editor for Experiments in Fluids.
Code availability The in-house PIV interrogation software used in this study will be made available at https:// git. soton. ac. uk/ jml1g 18/ pivto ols.

Ethics approval Not applicable.
Consent to participate Not applicable.