Fluid flow and mixing in a novel intermittently rotating bioreactor for CAR-T cell therapy manufacturing

This work explores the mixing and fluid dynamics in a novel bioreactor currently used for the automated manufacturing of CAR-T cell therapy, which offers a single-dose cure for several forms of advanced blood cancer. The cylindrical bioreactor has a low aspect ratio and a free surface. Agitation is achieved by intermittent rotation of the entire vessel around its central axis. No engineering characterisation has been conducted to date for this system in a bioprocessing context. The study examines the fluid dynamics problems of spin-up from rest and spin-down to rest. Novel Particle Image Velocimetry and Planar Laser-Induced Fluorescence data is presented alongside reflective flakes results, shedding light on the different transient flow regions inside the intermittently rotating bioreactor, and determining the timescales of macro-and micromixing. The results presented can be used to design a custom rotation pattern of the bioreactor for improved mixing performance during the cell expansion step.


Introduction
Chimeric Antigen Receptor (CAR) T-cell therapy is a personalised cell-based gene therapy, which offers a single-dose curative treatment for certain types of late-stage blood cancer.Since the regulatory approval of the world's first commercial CAR-T cell therapy (Kymriah®, Novartis) just six years ago, the cell and gene therapy sector has experienced unprecedented growth in that direction, with six commercially available CAR-T cell therapies to date, and more than 650 candidates in ongoing clinical trials (Gustafson et al., 2022).However, the complex nature of these biological medicines results in numerous manufacturing challenges, high risk of process failure and ultimately, in their significant cost.
To address the manufacturing challenges and decrease cost of production, efforts have been made to design novel automated platforms tailored for CAR-T manufacturing.However, due to non-standard geometries and agitation strategies, process translation can prove challenging.This is the case for the CliniMACS Prodigy® (Miltenyi Biotec), one of the few commercially available automated cGMP-compliant platforms for CAR-T manufacturing, which can accommodate all unit operations (except cryopreservation) in one closed system.While the Prodigy is gaining popularity in both clinical and industrial settings (Mock et al., 2016;Zhu et al., 2018), fundamental fluid and mixing dynamics characterisation of its bioreactor, required for process scaling and optimisation, is still missing.
At the core of the Prodigy's functionality is the CentriCult™ chamber, a vertical low-aspect-ratio flat-bottom cylinder, serving both as a bioreactor, and as a centrifuge, with a radius of 6 cm and maximum working volume of 300 mL (75% volumetric efficiency) (Richter 2014).In contrast to the well-studied stirred-tank reactor (STR), the Cen-triCult™ is equipped with neither an impeller, nor a sparger.Instead, it rotates intermittently around its central axis to mix the cell broth and relies on surface aeration.The change of rotation rate of the chamber (start/stop) produces a transient secondary flow.The pattern of secondary flow depends on whether the rotation rate of the bioreactor is increased, decreased or reversed: spin-up, spin-down or spin-over, respectively.Each of the corresponding secondary flows has different implications for the cell expansion step in terms of mixing, oxygen transfer, suspension efficiency, shear stress and power input.
While the CentriCult™ has unique geometry, similarities exist with the Taylor-Couette bioreactor, also known as the Taylor vortex bioreactor, used for stem cell differentiation and tissue constructs (Chaudhuri and Al-Rubeai, 2005).It consists of two vertical coaxial cylinders with a narrow gap between them, filled with culture media.This bioreactor creates Taylor-Couette flow, which has long been examined from a fluid dynamics perspective (Greenspan, 1968).By selectively rotating either the outer or inner cylinder, while keeping the other stationary, different flow patterns can be achieved, the implications of which have more recently been addressed in the context of bioreactors.When the outer cylinder is rotated, laminar Couette flow is established, providing an easily controllable dissolved-oxygen environment, suitable for culturing cartilage constructs, as demonstrated by Saini andWick (2003, 2004).On the other hand, Curran andBlack (2004, 2005) have shown that the Taylor-vortex flow regimes, induced by rotating the inner cylinder, while keeping the outer stationary, offer improved cell suspension and oxygen mass transfer, albeit at the expense of higher shear stress.No such fluid dynamics and mixing characterisation has been found in the literature for the single intermittently rotating vertical cylinder bioreactor of interest to this work.Nevertheless, the problem of fluid adapting to the changing of rotation rate of a cylindrical enclosure has been addressed in other disciplines, such as ballistics, aerospace engineering and oceanography.For practical purposes, only the cases of spin-up from rest and spin-down to rest will be examined in the literature review, as they are most relevant to the operation of the bioreactor in question.
Nonlinear spin-up (from rest) in a cylinder was first studied analytically by Wedemeyer (1964), who extended the founding theory laid by Greenspan and Howard (1963) for linear spin-up (i.e.not from rest).He showed how the transient flow evolves for a vertical flat-bottom cylinder of radius, R, and height, H, completely filled with a fluid of kinematic viscosity, ν, which is impulsively accelerated from rest to a constant rotation rate, Ω (i.e.non-linear spin up).The Reynolds number is defined in Eq. ( 1) and the flow remains laminar for Re < 3 × 10 5 (Schlichting and Gersten, 2017).Within the first few radians of rotation, a quasi-steady viscous boundary layer develops at the (vertical) sidewall and at each (horizontal) endwall, while the interior of the fluid remains inviscid.
The boundary layer at each endwall, known as the Ekman layer, is of the von Karman type and, during spin up, expels fluid radially outward under the influence of centrifugal forces.To balance the radial outflow, a small axial influx is established at the centre of the cylinder from the inviscid interior into the Ekman layer.This convection, known as Ekman pumping, results in a geostrophic flow, characterised by a cylindrical shear front of the tangential velocity profile, dividing the volume of fluid into an outer rotating and an inner non-rotating domain.The shear front propagates from the cylinder wall to the axis of rotation in time.According to Wedemeyer's model for a completely filled cylinder, the radial location of the front in Ekman time, T, (Eq.( 2)) can be described by an exponential decay law.When the front reaches the axis of rotation, the fluid has essentially achieved solid body rotation (SBR) and spin-up is completed on the Ekman timescale.Experimental (Weidman 1976a;Savas 1985) and computational (Watkins and Hussey 1977;Kitchens 1980) studies are in good agreement with Wedemeyer's prediction.
Spin-up from rest in a partially filled cylinder was first addressed computationally by Goller and Ranov (1968), who considered an additional secondary circulation, generated by the gradually deforming free surface.Their computational model was in close agreement with experimental data for the transient free surface shape and was further advanced by Homicz and Gerber (1987) to include the cases when the free surface intersects one or both endwalls.Choi et al. (1989Choi et al. ( , 1991) ) have further validated these findings with free surface shape measurements, as well as Laser-Doppler Anemometry (LDA) measurements of the tangential velocity.They found that the propagation of the tangential velocity shear front follows an exponential decay as predicted in Wedemeyer's model for a completely filled cylinder, provided that the free surface intersects neither endwall.Van De Konijnenberg and Van Heijst (1995), however, reported that their experimentally determined free surface vorticity profiles deviate from those obtained through Wedemeyer's model, which no longer accurately estimates the Ekman pumping for spin-up from rest in a partially filled cylinder.
In the presence of a free surface, the Froude number, Fr, defined in Eq. ( 3), where g is the gravitational acceleration and H L is the liquid height at rest, becomes an additional scaling parameter to be considered.Greenspan and Howard (1963) first suggested that the spin-up process duration is inversely proportional to Fr.This was later experimentally confirmed by O'Donnell and Linden (1991) and Kloosterziel and Van Heijst (1992) who concluded that the vertical motion of the surface must be responsible for the increase in the characteristic timescale.The latter authors also showed analytically that no such Froude number effect is expected for linear spin-down.
Beside the Ekman layer on the endwall, an additional vertical boundary layer forms at the sidewall during spin up, i.e. the Stewartson layer, which reinforce the Ekman circulation by extending the radial outflow from the edge of the Ekman layer into an upflow along the sidewall (Greenspan and Howard, 1963;Van De Vooren, 1993).However, during spin-down to rest, this vertical boundary layer behaves in a vastly different manner and has been the main focus of study in the literature.Due to the curvature of the sidewall, the decelerating flow in its proximity becomes unstable and takes the form of horizontal counterrotating vortex pairs, known as Taylor-Görtler (T-G) vortices, which are at first identical to those formed in the Taylor-Couette bioreactor (Curran and Black, 2004).These occur when Re > Re crit , where Re crit varies in literature between 140 < Re crit < 320 (Neitzel and Davis 1981;Mathis and Neitzel 1985;Cui 2004;Kim and Choi 2006).Euteneuer (1968Euteneuer ( , 1972) ) first studied this flow experimentally in a completely filled transparent cylinder with high aspect ratio (H/R = 4).Using reflective flakes to visualise the vortex tubes on the sidewall, they showed that their onset time and wavelength are scalable with Re.These findings have been confirmed, again with reflective flakes, in other high aspect ratio configurations, H/R > 2 (Maxworthy, 1972;Mathis and Neitzel, 1985) and computationally in infinitely long cylinders (Neitzel and Davis 1981;Kaiser et al., 2019).Kim and Choi (2006) have also provided an analytical model which estimates the onset time and associated wavelength in an infinitely long cylinder.
Nevertheless, the behaviour of the sidewall boundary layer after the onset of the vortices, is still not well understood.Euteneuer (1972) observed that for the full range of Re studied (1000 < Re < 20,000), the primary vortex pairs grow linearly at first, after which they lose shape and begin continuously restructuring into fewer larger pairs.They also showed that for Re > 8000, the wavelength growth rate drastically increases compared to Re < 8000; a phenomenon they explained with the onset of three-dimensional turbulence in an otherwise laminar boundary layer.The transition to turbulence in the presence of T-G vortices, according to Saric (1994), occurs through the appearance of a secondary instability, which corrupts the structure of primary pairs.More recently this has been studied computationally by Kaiser et al. (2019) in an infinitely long cylinder (3000 < Re < 28,000), who observed a secondary instability in the form of 'hairpin' structures, which leads the flow to transition.Moreover, multiple studies suggest that if the boundary layer Reynolds number, Re δ , is sufficiently high, Tollmien-Schlichting (T-S) waves may also be present and interact with the T-G vortices, resulting in different mechanisms of premature transition to turbulence (Görtler, 1960;Bennett and Hall, 1988;Saric, 1994).
For spin-down in a finite cylinder, the boundary layer forming at the (flat) endwall is of the Bödewadt type and again induces Ekman circulation but in the opposite direction to spin-up (i.e., radial inflow along G.G. Atanasova et al. the endwall and an axial efflux from the Ekman layer into the interior).Faller (1991), who studied the flow over a rotating disc presented a comprehensive computational study alongside dye visualisation experiments, showing that both von Karman and Bödewadt boundary layers become unstable at Re crit = 70 and Re crit = 15, respectively, and each one can form two distinct types of stationary spiral vortex rolls, referred to as Type I and Type II waves.They identified multiple mechanisms of interaction between the waves, which can lead the Ekman layer to transition to turbulence.Bassom and Hall (1990) have analytically studied the interference of these waves with T-S waves, which may also be present in the Ekman layer and provide alternative mechanisms for the transition.Savas (1992) studied Type I and Type II waves experimentally in the context of a completely filled cylinder (H/R = 2) spinning down to rest.They visualised the instabilities with reflective flakes, and with LDA measurements of the tangential velocity component showed that the associated disturbances are not confined to the Ekman layer but perturbate into the interior.In this respect they behave differently than in the rotating disc apparatus of Faller.However, Savas (1992) only focused on disturbances due to Ekman layer instabilities and did not explore possible interactions with the sidewall boundary layer, where T-G vortices are simultaneously present.
To the authors' knowledge, there are no studies in the literature which explore the interaction between these two boundary layers in low aspect ratio cylinders (H/R < 1).Furthermore, no studies have been published for spin-down to rest in a partially filled cylinder, where the free surface deformation also induces a secondary motion.Unlike the pseudo-steady state flow in a stirred tank, which is controlled by the periodicity of the impeller blade passage, the flow in an intermittently rotating cylinder is transient and, depending on the stage of the spindown process, different regions with laminar transitional or turbulent flow can simultaneously exist within the reactor.Secondary motion originates at the walls of the vessel and spreads to the interior of the fluid in time.For the range of Re the CentriCult™ is operated at (Re > 37,000), the boundary layer flow is laminar but prone to a number of different types of instabilities which can interact and lead to transition to turbulence.The mechanisms of this transition are still not fully understood.The interaction between these regions and the premature transition to turbulence can be beneficial for mixing if the agitation strategy is well designed.However, no studies have been found in the literature, which address the mixing in an intermittently rotating partially filled vertical cylinder.
The aim of this work is to characterise the complex transient flows during spin-up from rest and spin-down to rest from solid body rotation in a partially filled low aspect ratio cylinder (H L /R = 0.5), and to clarify the impact of different agitation strategies on the flow in a bioprocessing context, specifically focusing on the interaction of micro-and macromixing.This was achieved by high spatial and temporal resolution measurements, Particle Image Velocimetry (PIV) and Planar Laser-Induced Fluorescence (PLIF), combined with standard methods based on reflective flakes (RF).The findings presented are applicable to vertical cylindrical bioreactors of similar aspect ratio to the one studied here (H L /R = 0.5-1), the operation of which resembles the CentriCult™, and are scalable with aspect ratio, Reynolds number and Froude number.Adopting an improved mixing strategy through the use of a customdesigned rotation pattern of the bioreactor will lead to improved oxygen transfer and therefore, to higher yields of the expansion step of the CAR-T manufacturing process and a reduction of costs.

Methodology
Experiments were conducted in a cylinder with radius, R = 49 mm, partially filled with Milli-Q® water.The working volume of liquid of 190 mL for all experiments was selected to match the aspect ratio, H L /R, of the CentriCult™ at maximum working volume (H L /R = 0.5).This corresponds to a fluid height, H L = 25 mm for experimental conditions investigated.In the experimental setup, shown in Fig. 1, the cylinder is submerged in a large square tank of water to minimise the impact of optical distortion by the curved cylinder wall.The axis of rotation is fixed with a ball bearing on each endwall, integrated into a square tank, as shown in Fig. 1.An electrical motor is connected on top via a belt mechanism.The motor position is offset from the central axis, to allow injection of dye while the cylinder is in rotation.The position of the laser and camera can be swapped to investigate different cross-sections of the flow (i.e., vertical or horizontal).A summary of the range of experimental conditions studied is provided in Table 1.
Preliminary validation experiments were carried out with reflective flakes, used to visualise the Taylor-Görtler vortices on the wall of the cylinder, using only natural light.The experimental technique is described in detail by Savas (1985) and has been used to study this flow on many instances as outlined in the literature review.The classical aluminium flakes were replaced with Silver Pearl Lustre Pigments (Cornelissen) in these experiments as they settle more slowly.An image obtained with this technique is shown in Fig. 2a, where the counterrotating vortex tubes are clearly visible as bright and dark bands on the cylindrical wall.
PLIF was carried out in the vertical plane which bisects the cylinder in half and to achieve a high resolution the field of view was restricted to half of the cross-section (r/R = 0-1, z/R = 0-0.5).It was utilised as a  visualisation technique to study the flow during the full duration of spinup and spin-down, and also to quantify the mixing at each Re.Beside the reference height H L = 25 mm, an additional, higher aspect ratio, H L /R = 1, was also studied with PLIF to further test mixing scalability.A 50 µL drop of 0.5 mM solution of rhodamine 6G (Sigma) in Milli-Q® water was injected at the centre of the bottom at rest, which matches the injection position (not necessarily the timing) in the CentriCult chamber.That particular injection position was selected for the experiments also because when the rotation starts, the dye is distributed into a thin sheet to both boundary layers by the Ekman pumping, and when it stops, both regions of the flow were illuminated simultaneously.The initial concentration of injected dye was selected to avoid image saturation at any instant during the mixing process.Reference images for each experiment were taken at the start before injection of the dye and at the end when the dye is fully mixed.Rhodamine fluoresces in green light (532 nm) and has an emission wavelength in the orange spectrum, which was isolated with a monochromatic bandpass filter (575 nm).The experiment temporal resolution varied for different Re between 1 ms (Re = 25,043) and 10 ms (Re = 3756), while the spatial resolution was kept fixed at 4.2 × 10 -5 m/px.An image obtained with this technique is shown in Fig. 2b, where the T-G vortex pairs are visible as mushroom-like structures stemming from the sidewall.Note that Fig. 2a and 2b show the same condition but at different instants of the spin-down process, as the instabilities take some time to propagate away from the wall and become visible in the cross-section.A complete video of the PLIF experiment at Re = 17500 is provided in the supplementary material (17500_vid.avi).2D time-resolved PIV was carried out in both the vertical and horizontal cross-sections.The flow on the horizontal cross-section was studied at a single height, z = 16 mm (z/H L = 0.64).Fig. 1 shows the set up for vertical plane measurements, but the position of the camera and laser would be swapped when looking at horizontal planes.Neutrally buoyant non-fluorescent Polyamide Seeding Particles (PSP) with mean diameter of 5 µm were used (Dantec Dynamics).The temporal resolution varied for the range of Re and different planes studied between 200 µs (Re = 17530, vertical plane) and 3 ms (Re = 3750, horizontal plane).The final 8 × 8 pixels spatial resolution after processing varied between 0.28 mm (3.5 × 10 -5 m/px, Re = 17530, vertical plane) and 0.73 mm (9.1 × 10 -5 m/px, horizontal plane).For the vertical cross-section, the field of view is r/R = 0-0.5,z/R = 0-0.5, and for the horizontal, it is r/R = 0-1, ϴ = − 90 • to 90 • , where a cylindrical coordinate system with origin at the centre of the bottom wall is considered and positive rotation is anti-clockwise, as seen from the top of the vessel.

Results and discussion
Firstly, the core flow during spin-up from rest is addressed in Section 3.1, followed by the core flow during spin-down to rest from solid body rotation (SBR) in Section 3.2.1.The data is used to establish the timescales and scalability of each process for different Re.Then, in Section 3.2.2, the study focuses on the vertical boundary layer during spin-down and the generation of T-G vortices and turbulence.RF results are compared to data available in the literature for a completely filled high aspect ratio and novel PIV data is presented to characterise the vorticity intensity and decay timescale of T-G vortices.Finally, in Section 3.3, PLIF is utilised to capture the full length of spin-up and spin-down processes and to study the macro-and micromixing in the vessel.

Spin-up from rest
Profiles of the instantaneous tangential velocity, U θ , along the radius, are presented in Fig. 3a for different number of revolutions after the start of rotation for Re = 3756, Fr = 0.02, Ω′ = 100π s − 2 .The tangential component of the velocity was obtained from PIV measurements in the horizontal plane, at height, z = 16 mm (z/H L = 0.64).The computational model of Goller and Ranov (1968) for impulsive spin-up from rest in a partially filled cylinder separates the flow into two distinct regions: the Ekman boundary layer flow, and the flow in the inviscid interior, far from the Ekman layer.The free surface deformation also induces a secondary motion, in addition to the Ekman pumping.As a consequence the horizontal plane in these experiments was positioned away from both the Ekman layer and the free surface (z/H L = 0.64), where the tangential component is independent of the vertical coordinate and is thus more indicative of the core flow.Each profile is comprised of more than 6000 points, sorted by radial location in the horizontal plane.Data for solid body rotation (Nt = ∞) was obtained from a separate experiment, after the fluid had been allowed to spin-up for a sufficiently long time.To validate the state of solid body rotation, a straight line was fit through the tangential velocity profile, and the initial gradient showed deviation from the theoretical profile no larger than 0.8% for all experimental conditions investigated.
Initially, only fluid close to the wall is in rotation, while the core is at rest.As spin-up time advances, the vertical boundary between rotating and non-rotating fluid, characteristic of the geostrophic flow, propagates from the wall towards the central axis, until the entire fluid rotates as a solid body.In order to find the radial position, r*, of that propagating shear front, an 8th-degree polynomial was fit through each instantaneous tangential velocity profile.The shear front radial coordinate, r*, was taken when the normalised tangential velocity fell below 0.05, as indicated by the dashed line in Fig. 3a.
The position of the propagating front of tangential velocity, indicative of the progress of spin-up in time, is shown in Fig. 3b for different Re.The faster the chamber is spinning, the more revolutions it takes for the fluid to reach the steady state of solid body rotation, even though the absolute time is shorter (see Table 2).Nevertheless, when the time is normalised with the Ekman timescale (Eq.( 2)), see Fig. 3c, the spin-up process occurs on the same convective timescale for all Re, in agreement with Wedemeyer's analytical model for a completely filled cylinder.According to the model the position of the propagating front follows the exponential decay law of Eq. ( 4).However, the decay rate coefficient, b, is not known a priori.Greenspan (1968) reported that for spin-up which is not from rest, in an impulsively accelerated completely filled cylinder, the value of b can be 0.886, 1 or 1.1.In the practical case of interest, the cylinder is only partially filled, and it is almost impulsively accelerated from rest; a case for which decay rate coefficients are not reported in the literature.For a partially filled cylinder configuration Choi et al. (1991) suggest calculating the Ekman time from the model of a completely filled cylinder with a height H = 2H L , where H L is the liquid height.This is to account for the absence of the top Ekman layer, and therefore assuming a core flow at the free surface height.The values of the decay rate coefficient, b, for the experimental runs shown in Fig. 3c are included in Table 2.
In Fig. 3c it is also seen that the characteristic decay timescale, T d (i.e. 1/b), of the shear front radial coordinate becomes longer at larger Froude numbers (higher Re).This is due to the increasing impact of the free surface deformation on the flow, in line with the literature.The surface remains flat at solid body rotation at Re = 3756, however at Re = 25043, it is a steep paraboloid.The gradual deformation of the surface during spin-up generates a secondary flow, in addition to the Ekman pumping by the bottom wall, as outlined by Goller and Ranov (1968).Greenspan and Howard (1963) found analytically that for linear spin-up (i.e.not from rest), the decay timescale expressed in Ekman time increases linearly with Fr.This was later confirmed experimentally by O'Donnell and Linden (1991) for 0.38 < Fr < 15.31.However, no reports have been found in the in the literature for how that relationship might change for spin-up from rest (i.e.nonlinear spin-up), which is the case of the current work.The spin-up completion Ekman time, T SU95 (t SU95 in dimensional units) was taken when each decay curve falls below 0.05, meaning spin-up is 95% complete (99.75% of the volume of   Greenspan and Howard (1963).It is worth to mention that despite the increase in the Ekman spin-up time, the dimensional completion times decrease with Re.
3.2.Spin-down to rest from solid body rotation

Core flow
In this section, analogous analysis of the core flow is carried out on spin-down to rest from solid body rotation (Ω′ = − 100π s − 2 ), while the boundary layer flow is examined later, in Section 3.2.2.
Instantaneous tangential velocity profiles are shown in Fig. 4a for the same values of Re and Fr as spin-up.While fluid close to the wall decelerates and is brought to rest, the core remains in solid body rotation and its angular velocity decays in time.
To find the rate of decay of angular momentum of the fluid core, a straight line was fit to the tangential velocity profile at each timestep, from the axis of rotation to r/R = 0.75.The gradient of that line represents the normalised angular velocity of the fluid core, which retracts as the sidewall boundary layer expands.As will be shown in Section 3.3, the sidewall boundary layer remains confined to r/R > 0.75 for the low aspect ratio configuration with a free surface, investigated in this work.
The decay of the angular velocity, ω, of the inner inviscid core with number of revolutions is shown in Fig. 4b for different Re.Similarly to spin-up from rest, the faster the chamber was spinning, the more revolutions it takes to complete the process.
Nevertheless, the four curves fully scale when the decay is plotted against the Ekman time, indicating no dependence on Fr, unlike spin-up from rest.This is in agreement with the findings of Kloosterziel and Van Heijst (1992), who modelled this analytically for linear spin-down (not to rest).In the early stages of the process, the spin-down decay curve is steeper than the analogous curve for spin-up, which indicates that angular momentum is lost faster than it is built during spin-up.The curve does, however, appear, to flatten in the later stages of the process and to roughly match the spin-up curves at large Ekman times.
A single exponential decay can no longer describe the entire time range of data because of the steeper decay initially, and the flattening of the curve in the later stages.Therefore, in order to model the decay of angular velocity in the later stages more accurately and find the spindown completion times, t SD95 and T SD95 , an exponential decay curve, modelled according to Eq. ( 5), was fit through all the data sets for T > 0.8 (see in Fig. 4c).
Values of a = 0.59, b = 1.60 were found, and using those, the dimensionless spin-down time, T SD95 was calculated, at which the value of the angular velocity falls below 5% of that at solid body rotation.This is presented in Table 3, along the dimensional spin-down times at each Re, denoted as t SD95 .In accordance with the analytical models by Benton (1973) for spin-up and spin-down between infinite coaxially rotating plates, non-linear spin-down is overall slower than non-linear spin-up, despite the early stage of rapid decay in angular velocity.The characteristic Ekman times for spin-up and spin-down between infinite plates, reported by Benton (1973), however, are much longer than those for fluid enclosed in a cylinder.

The unstable boundary layer
This section focuses on the generation of instabilities in the vertical boundary layer during spin-down to rest and examines their  implications for localised mixing (micromixing).
Reflective flakes were used to visualise the vortices as they form at the wall and their onset wavelength and time were measured (Fig. 5a and 5b, respectively).Results from this study are compared to computational results (empty markers), experimental data (filled markers), and the theoretical prediction (dashed line) found in literature.
The decrease of the onset wavelength, λ c , with increasing Re in Fig. 5a is robustly predictable and the experimental data points from this study follow nicely the theoretical model of Eq. ( 6) provided by Kim and Choi (2006), where A = 28.56.
A second higher dashed line is plotted in Fig. 5a with magnitude 2 × λ c .This is to compare with the early work by Euteneuer (1968), who first used reflective flakes to visualise Taylor-Görtler vortices in a rotating cylinder and defined the onset wavelength as the combined distance of one bright and one dark band (as seen in Fig. 2a), which is equivalent to the diameter of a vortex pair.On the other hand, Maxworthy (1972) and Mathis and Neitzel (1985) employed a photodiode array to differentiate between bright and dark bands and calculated the wavelength (of a pair) as the distance between two adjacent peaks in the photodiode voltage.Since their results for the wavelength are exactly two times higher than those of Euteneuer (1968), Kim and Choi (2006) assumed that with the photodiode technique, two counter-rotating vortices are imaged as a single vortex, and thus the photodiode results should be scaled down by a factor of 2.
The effect of gradual (as opposed to impulsive) angular deceleration of the chamber, Ω' = − 2π rad/s 2 , on the onset wavelength was also studied using reflective flakes.This data is shown in Fig. 5a as magenta markers for three values of Re.The measured wavelength deviates from the impulsive reference line (dashed line) and visibly increases as a result of the slower stop of the wall.This is in line with Weidman (1976b), who reported gradual as opposed to impulsive deceleration of the cylinder leads to an increase in the critical Reynolds.For example, at Re = 25043, λ c /R increases from 0.0314 for an impulsive stop to 0.0399 for the slowest deceleration.A λ c /R of 0.0399 would be achieved with impulsive stop at Re of 19150, according to the theoretical model of Kim and Choi (2006).This means that a higher deceleration (faster stop) of the chamber would be more beneficial for mixing as this is linked to smaller and faster-rotating local eddies which increase the local rate of viscous dissipation of kinetic energy and therefore enhance micromixing.
The onset time of the vortices, t c , plotted in Fig. 5b, is also very scalable with Re and data available in the literature widely agrees.At Re > 10 3 the theoretical correlation provided by Kim and Choi (2006) and shown in Eq. ( 7), may be used to find the dimensional onset time of instabilities for any flat-bottomed cylinder (A = 9.4).The higher the Re number, the sooner the vortices appear.At very high Re, the instabilities appear instantaneously after the stop.However, experiments with reflective flakes (data not included) showed the onset of instabilities can occur before the cylinder has reached complete stop, should the deceleration be sufficiently low and the Reynolds number sufficiently high.
The reflective flakes technique can provide a quick and easy validation of the experimental setup; however, it provides only visual information of the structures close to the wall during the early phase of spin down and cannot be used to gain further insight into the flow dynamics and mixing in the interior.All subsequent analysis presented is derived from time-resolved vertical-plane PIV and PLIF data, which, to the authors' knowledge, has not been presented in the literature for this flow before.2019), the exact moment of transition to turbulence occurs at the moment of 'knick', described by Euteneuer (1972).This corresponds to the moment the   original coherent structures of vortex pairs lose shape and begin rapidly reshuffling, as their characteristic wavelength begins growing at a much higher rate.This can be seen in Fig. 6, at Re = 15026, for example.The first frame, at T = 0.021, illustrates the last moment the primary structured pairs are seen, while in the next timeframe (T = 0.039), the vorticity structures observed are random and the thickness of the boundary layer has grown significantly in proportion to the increasing wavelength.Therefore, the moment of 'knick' has passed.
In order to gain further insight into the intensity of the instabilities, the stage of transition to turbulence, and the subsequent decay of turbulence, points which contain vortices must be differentiated from experimental noise.Several vortex identification methods were tested (Q-criterion, Omega method, Q-R method), and the λ 2 technique was found most effective.This method relies on the fact that a pressure minimum is present at the core of a vortex, and this can be identified if both eigenvalues of the matrix S 2 + A 2 are negative, where S and A are the symmetric and antisymmetric components of the velocity gradient tensor (Jeong and Hussain, 1995).In the current work, a threshold of ) 2 was selected to identify the vortices and to filter out erroneous local gradients due to noise in the instantaneous velocity fields.
The space-average normalised instantaneous vorticity, ωθ /Ω, of all points in the plane, for which λ 2 is smaller than the threshold, is plotted against the Ekman time in Fig. 7a.This allows a characterisation of the strength of the T-G vortices in time.Each vorticity curve peaks shortly after the stop of rotation, slightly earlier with increasing Re.The corresponding vorticity contour maps are shown in Fig. 6 at T = 0.078, T = 0.039 and T = 0.021 from low to high Re, and they represent the last moment the primary structured pairs are observed, and the instant when the boundary layer thickness begins increasing at a faster rate.It therefore becomes evident that the peak of average vorticity in the vertical plane coincides with the moment of 'knick', described by Euteneuer (1972) and Kaiser et al (2019) and marks the onset of turbulence in the sidewall boundary layer.The time of peak average vorticity, t onset , is included in Fig. 7b for the range of Re investigated (orange square markers), while the corresponding numerical values can be found in Table 4.It is also worth noting that the dimensional peak average vorticity in the vertical plane, ωθ (s − 1 ), was found to vary linearly with the Reynolds number in the range of Re studied (i.e.ωθ = 0.005Re − 8.176).
The decay of average vorticity in the plane occurs on a timescale much shorter than the convective Ekman timescale, on which the rotation of the core decays.This can be visualised better if the average vorticity in the plane is multiplied by the area of the cross-section pertinent to that vorticity, A ωθ .The resulting curves are shown in Fig. 7b for the range of Re considered.As energy dissipates and the vorticity decays, more points in the plane fall outside the λ 2 threshold.Eventually, each curve will reach 0, which means any remaining vorticity in the plane can no longer be differentiated from experimental noise.
As can easily be seen in Fig. 7b, the total vorticity drops to/close to 0 by Ekman time 0.5 for the entire range of Re, while around that Ekman time, the fluid core still has about 30% of its original angular velocity (Fig. 4c).To find the timescale of vorticity decay and compare it to the Ekman spin-down time, an exponential decay was fit to each curve in Fig. 7b from the peak onwards.The vorticity decay time, T decay , was calculated as the instant at which 99% of the total vorticity has been dissipated (i.e.time corresponding to 99% of the area under the curve being covered).This is given for each Re in Table 4 with the corresponding dimensional time, t decay , alongside turbulence onset times, found from the corresponding peak of vorticity.By the time 99% of the vorticity has decayed below detectability the boundary layer begins relaminarizing and therefore its potential for micromixing due to the T-G instabilities has been largely exhausted.At that point, rotation can be restarted in the same direction to conserve the remaining rotation of the core and begin another spin-up.The remaining angular momentum of the core, ω/Ω, at T = T decay , which would be conserved, is found from Fig. 4c, and is included in Table 4.The last column in Table 4 denotes the percentage of the spin down time over which the TG instabilities are  active and can effectively enhance micromixing.From this table a trend becomes apparent: the higher the Re, the longer it takes in Ekman time for the total vorticity to decay, therefore, the portion of spin-down time, during which micromixing in the vertical boundary layer is enhanced due to the turbulent T-G instabilities, is larger.This indicates that spindown becomes progressively more efficient at mixing with increasing Re, not only as a result of the higher intensity of the instabilities, but also due the longer time taken to dissipate the turbulence of these instabilities.

Mixing study
Complementary to the PIV experiments in Section 3.2, PLIF was employed to gain insight into the later stages of spin-down to rest, the macro-and micromixing dynamics and the role of the Ekman boundary layer.
Fig. 8 illustrates PLIF data at different stages of spin-down to rest at Re = 7513, Fr = 0.1 (Ω′ = 100π).As mentioned in Section 2, a 50 µL drop of concentrated rhodamine solution was injected at the centre of the bottom at rest.The dye was then distributed in a very concentrated thin film close to the wall by the Ekman convection during spin-up.The free surface at this Froude number remains practically flat and is at the top boundary of each subfigure; the left boundary is the axis of rotation, and the right boundary is the sidewall.The colormap shows the concentration of dye in each pixel normalised by the initial concentration in each experiment (i.e. at the moment of injection, C = 1).At the beginning of spin-down, when the structured primary T-G vortices are first seen in the vertical cross-section, the dye entrained by the vortices is still around 50% of the initial concentration, even after the long spin-up time (ref.to Fig. 8a).
Fig. 8a shows the moment right before the 'knick' when the average vorticity in the vertical plane peaks (from PIV) and the T-G vortices are last seen as structured pairs.At the 'knick' the flow in the sidewall boundary layer transitions to turbulence, as discussed in Section 3.2.Fig. 8b-8e (T = 0.13-0.48)capture the complex mixing dynamics during the period of turbulence decay.Throughout the mixing process three regions can be distinguished which interact with each other: the sidewall boundary layer (right side), where the T-G instabilities decay, the endwall boundary layer (bottom), which is dominated by the (clockwise) Ekman convection, and the inviscid interior (top-left side).Mixing takes place in the boundary layers, while the inviscid interior is a mixing dead zone.In agreement with PIV data (ref.to Fig. 6), the sidewall boundary layer expands radially towards the axis of rotation.At the same time the endwall layer is expanding towards the free surface as a result of the slow axial motion generated by the Ekman convection.This is well captured by the thin horizontal strip of dye, outlining the interface between the endwall layer and the inviscid interior, which is seen to travel towards the surface (from one plot to the next).At this Re, the interface, outlined by the strip of dye, is nearly flat (i.e.uniform with radial distance), and only a slight undulation can be seen close to the interface with the sidewall boundary layer.The thin film of dye is also seen to curve axially down past the sidewall layer due to the Ekman circulation loop slowly closing clockwise.The three distinct regions can be seen in Fig. 8f for the last time.At that point, the remaining vorticity has completely dissipated based on PIV results discussed above, and radial and axial velocity components in the sidewall layer are very low in comparison to the tangential one.In Fig. 8g, the dye close to the sidewall, which has been mixed by the T-G vortices, is slowly drawn into the Ekman layer.At this instant the Ekman circulation takes over the corner of the cylinder, which was previously within the sidewall BL.At that point, the sidewall boundary layer relaminarises, and the convective motion (axially down) begins to dominate the flow in proximity of the lower sidewall.A single slow surface vortex remains visible in the top-left side of Fig. 8g-8h, which eventually also decays, and the dye entrained in it is drawn into the Ekman convection loop (Fig. 8i).The pattern of dye freezes when the fluid comes to complete rest (Fig. 8j), which happens in very long Ekman times (T > 6), long after all velocity components, including the tangential one, have decayed below detectability.At this point, any remaining locally concentrated dye (i.e. as seen at the surface in Fig. 8j) mixes only through diffusion.A complete video a of the PLIF experiment at Re = 7500 including early and late stages of mixing is provided in the supplementary material (7500_vid.avi).
In addition, Fig. 9 shows the same instant in Ekman time, T = 0.41 for Re = 3756 (Fig. 9a) and Re = 25043 (Fig. 9b); the lowest and highest Re considered.The horizontal dashed line marks the thickness of the endwall layer, δ E , and the vertical dashed line indicates the thickness of the sidewall boundary layer, δ S .In Fig. 9a, the interface of the endwall layer with the inviscid interior, illuminated by the bright strip of dye, is relatively uniform in radius, however, clear oscillations are present in Fig. 9b near the interface with the sidewall layer.These perturbations are a result of the high intensity turbulence present at high Re in the sidewall boundary layer, which affects the region of the endwall layer in its proximity.
Analysing how the thickness of each boundary layer changes in time can shed more light on their interaction and the overall mixing dynamics.In order to find δ E , points at the interface of the endwall layer with the inviscid interior were selected with an image processing technique based on the local brightness gradient, but only in the range of r/R < 0.25, where that interface remains smooth and nearly flat for the entire range of Re considered.Results for the endwall boundary layer are shown in Fig. 10a.The sidewall boundary layer thickness, on the other hand, varies with the axial coordinate, z, due to the onset of turbulence.As a consequence points at the interface were selected across the entire liquid height, where the border of the sidewall and endwall boundary layers was clearly visible.In the early stages of the spin-down process, the vertical boundary layer grows at a very fast rate and a higher temporal resolution was used in the PLIF measurements than for later stages of spin-down.The variation in time of the Ekman and sidewall boundary layers are shown in Fig. 10a and 10b, respectively.Fig. 10a shows the thickness of the endwall boundary layer, normalised by the liquid height, in Ekman time for the range of Re and aspect ratios studied.To compare the rate of expansion of the endwall boundary layer to the loss of angular velocity of the fluid core during spin-down, a reference dashed curve corresponding to 1 -ω/Ω is also shown in Fig. 10a.For comparison, the vorticity decay timescale is also included as two separate grey bands in the background.The thinner band (on the left) close to T = 0, indicates the window in which the average vorticity in the vertical plane reaches its peak, and corresponds to the 'knick' and the onset of turbulence due to T-G instabilities.The wider band (on the right), at 0.35 < T < 1.1 indicates the window in which 99% of the total vorticity of the T-G instabilities has decayed for the range of Re, and the flow relaminarises.The hue in the grey bands darkens with increasing Re (i.e., full decay of vorticity occurs over a longer period at higher Re, cf.Table 4).
It is immediately seen that the Ekman layer growth rate in Ekman time is the same across the range of Re, and that it matches closely the rate of decay of rotation of the fluid core, illustrated by the dashed curve.The timescale of convective mixing is therefore the spin-down timescale already discussed in Section 3.2.1.Based on Fig. 4c, a non-dimensional spin-down time, T SU95 , of 2.68 was calculated as the Ekman spin-down time for all Re, at which the angular velocity of the core has decayed to 5% of its initial value during solid body rotation.At the same Ekman time, the thickness of the endwall boundary layer has increased to 95% of the liquid height, therefore the inviscid interior, which does not promote any mixing, takes up < 5% of the volume of liquid.Compared to the Ekman timescale, on which rotation decays and convective mixing is completed, the vorticity decay of the TG instabilities in Ekman time is significantly shorter and is dependent on Re (cf.Section 3.2.2).At the lowest rotation speed (Re = 3456), the TG instabilities vorticity has fully decayed at T decay = 0.35, which corresponds to a point on the lefthand side of the wider grey band.At that time, the endwall layer thickness is still in the stage of rapid growth.On the other hand, when a greater  The initial steep increase in δ S occurs within the thinner grey band, which marks the onset of turbulence.This is in line with the observations of Euteneuer (1972) and Kaiser et al (2019), who reported that the characteristic T-G wavelength growth rate, and therefore the sidewall boundary layer thickness, suddenly increase at the moment of 'knick'.Shortly after, as seen in Fig. 10b, the sidewall boundary layer reaches a nearly constant maximum thickness, which is maintained while the vorticity of the T-G instabilities decays.After the turbulence has decayed, the sidewall boundary layer relaminarises, and a single slow surface vortex remains visible, while the Ekman circulation extends all the way to the sidewall and takes over the corner, previously occupied by the ).The wavelength of the single surface vortex that remains was selected as the characteristic δ S for later Ekman times after the turbulence had decayed, as it is the only remnant of the original T-G instabilities, while the rest of the fluid near the sidewall is entrained by the Ekman convection.Similarly to the T-G instabilities and measurements of their vorticity in the vertical plane, the surface vortex persists for longer Ekman times at higher Re.Finally, at very long Ekman times, the surface vortex decays as well, and the curves drop to 0 because the dye circulated in the vortex is drawn into the main Ekman circulation loop (visualised in Fig. 8h-8j).
To quantify the total level of mixing, the normalised instantaneous brightness, b, was calculated for each pixel in time, as shown in Eq. ( 9).M is the value of the instantaneous brightness of a pixel, which can vary between 0 and 65535.M 0 is the value of the pixel brightness before the injection of dye, and M 1 refers to the final brightness when the dye is completely mixed.
A pixel is qualified as mixed, when |b| < 0.4, which corresponds to a +/− 40% deviation from the final normalised value of brightness for that pixel when the dye is completely mixed.To put this threshold in perspective, instantaneous values of |b| in the order of thousands were measured for each experiment shortly after the injection, when the dye is concentrated in a very small volume.The maximum |b| recorded was 3500, which translates to a maximum instantaneous deviation from the final pixel brightness of 350000%.The fraction of mixed pixels from the total number of pixels in the cross-section, γ M , represents the instantaneous level of mixedness.The fluid is mixed at mixing time, T M (t M in dimensional units), when 95% of the pixels in the cross-section are considered mixed, (i.e.γ M = 0.95).Fig. 11a shows how mixing is achieved in Ekman time for four different rotation speeds (15 rpm-100 rpm) and two aspect ratios (H L /R = 0.5 and 1).Each curve is the average of three experimental repeats.24,000 images in a sequence were analysed to produce the individual mixing curve for each experimental repeat.At higher Re, the variability was smaller and mixing characteristics were more reproducible.The standard deviation, σ, was between 0.02 and 0.06 for the range of Re studied.The green vertical dashed line marks the moment of start of rotation for all Re.The injection occurs just before that, and the dye is then spread in a thin layer over the wall during rotation.T = 0 is the moment the cylinder wall comes to a complete stop (Ω′ = − 100π).The red vertical dashed line corresponds to T SD95 when 95% of the angular velocity of the core has decayed, and the endwall boundary layer has expanded to 95% of the liquid height.The horizontal dashed line shows the threshold at which the fluid is taken to be mixed, γ M = 0.95.They grey bands illustrate the onset of turbulence (left) and the decay of turbulence (right) across the range of Re considered, where darker hues indicate higher Re.
The plot is divided into mixing during spin-up (T < 0) and mixing during spin-down (T > 0).During spin-up, each mixing curve plateaus at a constant low value of γ M between 0.05 and 0.1 around the spin-up Ekman time, -T SU95 (found in Table 2), when the fluid reaches solid body rotation.As expected, little micromixing occurs while the fluid is gaining angular momentum.At the moment of stop of rotation, the fraction of mixed pixels begins rapidly increasing at different rates for different Re.The rapid increase in the level of mixing slows down around the time vorticity decays at that Re, which is indicated by the wider grey band.At lower Re, vorticity decays earlier, therefore the mixing rate slows down earlier (in the lighter end of the grey band) than for higher Re (which happens in the darker end of the grey band).By T SD95 = 2.68 (red dashed line), 95% of the rotation of the core has decayed, and the endwall boundary layer has expanded to 95% of the liquid height.Spindown is therefore 95% complete, however, the mixing curves show that only the highest rotation speed (Re = 25043) has resulted in >95% total level of mixing.At the lowest speed studied (Re = 3756), the total level of mixedness is below 60% at the time by which the fluid has practically come to rest.Any further mixing from that point onwards occurs mainly through diffusion and its rate is very slow, as can be seen right of the red dashed line.
The nature of spin-down to rest from solid body rotation dictates that the core of fluid spins down as a solid body (the rotation of which decays in Ekman time).It thus represents a mixing dead zone, as outlined earlier in Fig. 8, which slowly shrinks on the Ekman timescale.The inviscid interior is only fully recirculated at T SD95 , as seen in Fig. 10a.Therefore the convective Ekman timescale controls macromixing during spin-down to rest from solid body rotation.This is confirmed in Fig. 11a, where the rate of mixing at Re = 25043 is seen to match the rate of expansion of the endwall layer (Fig. 10a) and the decay rate of the angular momentum of the interior fluid, (Fig. 4c).For this condition the mixing curve nearly passes through the intersection point of the red vertical dashed line, corresponding to the end of spin down, and the horizontal dashed line, showing to the fully mixed condition, i.e.T M = T SD95 .This is further confirmed by the mixing curve at Re = 25043 for the higher aspect ratio, H L /R = 1, which superimposes to the corresponding curve for the smaller aspect ratio, H L /R = 0.5.
These findings can be used to estimate the dimensional macromixing time, t M , from the macromixing Ekman time, T M , in any cylinder for any Re which is high enough, so that the rate of micromixing due to turbulent instabilities is higher than the rate of convective mixing due to the Ekman circulation.For example, from Fig. 11a it is evident that the minimum Re from the range tested satisfying this condition is 25043.The measured non dimensional and dimensional mixing times for Re = 25043 from Fig. 11a are T M = 2.34 and t M = 36 sec, respectively.These are very close to the nondimensional and dimensional spin down times T SD95 = 2.68 and t SD95 = 41 sec confirming that for such high Re, T SD95 can provide a reasonable approximation of T M .On the other hand, if micromixing in the boundary layers due to turbulence is slower than convective mixing due to the Ekman pumping (e.g.Fig. 9a vs Fig. 9b), then the boundary layer will not be homogenised before it reaches the free surface, at which point the fluid has run out of momentum.In this case, final local homogenisation is controlled by diffusion and T m → ∞.
In order to compare the rate of micromixing in the boundary layer regions, which varies with Re in Ekman time, to the rate of convective mixing due to the Ekman circulation, which is constant for the range of Re studied, the effect of the inviscid core mixing dead zone must be excluded when estimating the mixing time.For that reason, a second mixedness index based on a smaller area of the vertical cross-section including only the bottom half of the endwall layer (and the corner shared with the sidewall boundary layer) was further analysed.This area is delimited horizontally by the full radius, r/R = 0-1, and vertically as shown in Eq. ( 10), to match the endwall boundary layer thickness increase in time.
The level of mixing in the boundary layer is quantified by the fraction of mixed pixels from the total number of pixels within the area defined in Eq. ( 10), γ m .Similarly to Fig. 11a, in Fig. 11b, a horizontal dashed line, a red vertical dashed line and grey bands are used to visualise the mixed condition, the spin down time and the T-G vortices decay time, respectively.The solid black curve represents the convective (macromixing) timescale on which rotation of the core decays, as the endwall boundary layer expands and the inviscid interior mixing dead zone is depleted by the Ekman convection.
Fig. 11b illustrates the competition between micromixing due to the turbulent instabilities and macromixing due to the Ekman convection.At high rotation rates (e.g.Re = 17530 and Re = 25043), the scale and intensity of turbulence are sufficiently high, and the resulting rate of micromixing is higher than the rate of convective mixing due to the Ekman circulation, illustrated by the solid black line.In that case fluid is mixed locally in the boundary layers before convection is completed.Micromixing times of T m = 0.41 for Re = 25043 (t m = 6.32 s for H L /R = 0.5 and t m = 12.6 s for H L /R = 1) and at T m = 0.80 (t m = 14.74 s) for Re = 17530 were measured.Fig. 9b illustrates the exact instant in time when fluid in the area of the boundary layer examined qualifies as mixed at Re = 25043.The mixing curves for the higher and lower aspect ratio at Re = 25043 almost overlap and result in identical Ekman mixing time, which shows that the level of micromixing is independent of aspect ratio in the range examined.On the other hand, if the micromixing rate is lower than the rate of convective mixing (e.g.Re = 3756), the boundary layers fail to mix the fluid at a fine scale as spin-down progresses (i.e.Fig. 9a), and total mixing will only be achieved on the diffusive timescale, long after macromixing is complete and the fluid has come to rest.This is particularly well illustrated by the curve at Re = 7513, where at first, the rate of micromixing is just high enough to match that of convection.However, in later stages it slows downs and the level of mixedness plateaus at around 80%.This happens around the middle of the wider grey band, when vorticity in the vertical plane decays, followed by relaminarisation of the flow.In the absence of the turbulent instabilities, the rate of micromixing is surpassed by that of the Ekman convection, and the fluid comes to rest before complete mixing is achieved.

Conclusion
The transient flow during spin-up from rest and spin-down to rest to/ from solid body rotation (SBR) in a cylinder has been characterised by several studies, however to the best of the author's knowledge, none have addressed the flow from a mixing perspective.While the first appearance of the Taylor-Görtler (T-G) instabilities during spin-down is well characterised in the literature, very few studies focused on the later stages of spin-down, after the transition to turbulence, which are of particular interest to this work due to their implications for mixing.The present study offers novel experimental data and provides insight into the interior flow, boundary layer interaction and mixing dynamics during spin-up and spin-down in an intermittently rotating cylinder bioreactor with a small aspect ratio and a free surface.
In this work it was found that during both spin-up and spin-down, the fluid is divided into boundary layer regions (endwall and sidewall), where mixing takes place, and the inviscid interior fluid, which is a no mixing zone that is slowly depleted by the Ekman convection.Little overall mixing is generated by the stable convective (geostrophic) flow during spin-up, captured in the schematics in Fig. 12a.However, shortly after the stop of rotation, a steep increase in the overall level of homogeneity was observed for all Re.This was related to the transition of the sidewall boundary layer to turbulence under the influence of the T-G instabilities.Their growth, transition to turbulence and dissipation are the driving factors of mixing in the expanding boundary layer regions.In later stages of spin-down, the vorticity decays (earlier for lower Re), the flow relaminarises, and the remaining driving factor of mixing is the slow Ekman convection, which has little benefit for local small-scale homogenisation.The intensity and duration of turbulence were found to increase with Re.For Re = 3700, the portion of spin-down time during which the micromixing is enhanced by turbulence, is around 9%, with this number increasing to 35% at Re = 20,000.
Since micromixing is confined to the sidewall and bottom boundary layers, which thicken in time, their interaction during spin-down is of particular interest to this work (see schematics in Fig. 12b).It was found that the endwall (Bödewadt type) boundary layer expands towards the free surface in Ekman time, and its thickness is inversely proportional to the decaying angular velocity of the inviscid fluid core.On the other hand, the turbulent sidewall boundary layer expands (radially towards the axis of rotation) at a much higher rate than the endwall layer at first.However, its expansion is limited to about 0.25R, where the thickness plateaus while the vorticity decays independently of the Ekman convection.It is noteworthy that during that stage, at high Re, the region of the endwall layer close to the interface with the sidewall layer is seen to be impacted by the turbulent oscillations, which additionally enhances mixing in the endwall layer.Once turbulence has decayed completely, the Ekman circulation is seen to extend all the way to the sidewall, where the T-G vortices were previously identified, as the endwall boundary layer reaches the free surface (the inviscid core is depleted) and the fluid comes to complete rest.This happens on the spin-down Ekman timescale, which was found to be independent of Re, and an Ekman spin-down time of T SD95 = 2.68 was calculated for all rotation speeds.
It thus becomes evident that the mixing process in the intermittently rotating bioreactor is the result of a fine balance between local fine-scale (micro) mixing in the boundary layers, generated by the T-G instabilities, which is dependent on Re, and the large-scale convection, associated with the Ekman pumping, which is independent of Re.If the Reynolds number is not sufficiently high, the intensity and duration of the T-G instabilities are insufficient to mix the fluid locally, while the endwall layer expands and recirculates it.Therefore, the resulting rate of micromixing is lower than the rate of convective mixing, which results in extremely long macromixing times governed by the diffusive timescale.On the other hand, if the Reynolds is high enough (i.e.Re > 17,000 from the conditions tested), fluid in the boundary layers is mixed locally first and when the Ekman convection is completed, macromixing is also complete.In that case, the spin-down Ekman time, T SD95 = 2.68, can provide a reasonable approximation of the Ekman macromixing time, T M , from which the dimensional macromixing time can be calculated for a cylinder with a similar aspect ratio to those studied here.
Finally, it is worth to mention that operating the bioreactor at very high Re will result in an improvement in the micromixing due to the higher intensity and duration of turbulence, however the macromixing Ekman time, T M , will remain the same since the inviscid interior is only depleted on the Ekman timescale.Thus, it becomes evident that there must also be an upper boundary for the optimal Re at which the bioreactor is to be operated, beyond which shear stress and power input increase, but micromixing can no longer improve, and macromixing is still limited by the slow depletion of the inviscid interior by the Ekman convection.In order to avoid this limitation altogether, the rotation pattern of the bioreactor must be designed in such a way to generate flow which disrupts the inviscid interior early on, so that the enhanced mixing due to turbulence is not confined to the boundary layers.This will be the focus of future work, where spin-down, not from solid body rotation, and spin-over (reversal of the direction of rotation) will be examined.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 3 .
Fig. 3. a) Instantaneous tangential velocity radial profiles during spin-up form rest, z/R = 0.64 (Re = 3756); b) Radial position of the propagating front of tangential velocity as a function of number of revolutions for different Re; c) Position of the propagating front as a function of Ekman time for different Re.

Fig. 4 .
Fig. 4. a) Instantaneous tangential velocity radial profiles during spin-down to rest from sbr,z/R = 0.64 (Re = 3756); b) Decay of the angular velocity of the interior core of fluid vs. number of revolutions for different Re; c) Decay of the angular velocity vs Ekman time for different Re.

Fig. 6
shows an example of the instantaneous normalised tangential vorticity in the vertical cross section, ω ϴ /Ω, at different Ekman times for three different Re.In agreement with RF experiments, the onset of T-G vortices occurs earlier for increasing Re.The different-sized primary pairs of counter-rotating T-G vortices can be seen at T = 0.078, T = 0.039 and T = 0.021 from low to high Re.The instabilities first appear smaller and are seen to briefly grow before they lose shape and the local flow transitions to turbulence (e.g.: T = 0.059 for Re = 10017 and T = 0.039 for Re = 15026).According to Kaiser et.al. (

Fig. 5 .
Fig. 5. a) Wavelength of Taylor-Görtler instabilities during spin-down to rest from solid body rotation; b) onset time of Taylor-Görtler instabilities during spin-down to rest from solid body rotation.

Fig. 6 .
Fig. 6.Instantaneous normalised vorticity in the vertical plane obtained with PIV; x-axis shows the Ekman time of the snapshot, and the y-axis shows the Re.The field of view is r/R = 0.5-1 and z/R = 0-0.5.

Fig. 7 .
Fig. 7. a) Mean vorticity in the cross-section in Ekman time for different Re; b) Total vorticity in the cross section as a function of Ekman time for different Re; All curves are obtained independently for two experimental repeats, and only the average is plotted.

Fig. 10 .
Fig. 10.a) Average axial thickness of the endwall boundary layer during spindown to rest as a fraction of the liquid height against Ekman time; b) average radial thickness of the sidewall boundary layer during spin-down to rest as a fraction of the liquid height in Ekman time.

Fig. 11 .
Fig. 11.a) Fraction of mixed pixels from the total vertical cross-section, showing the rate of macromixing; b) fraction of mixed pixels from the boundary layer cross-section, showing the different rates of micromixing.

Fig. 12 .
Fig. 12. Schematics of the flow during spin-up from rest and spin-down from rest, (a) and (b) respectively.Dashed symbols show an arbitrarily later instant in time.

Table 1
Range of operating conditions studied and experimental methods used.

Table 2
Characteristic spin-up time for different Re and Fr.
N [rps]Re[-] Fr [-] b [-] T d [-] T SU95 [-] t SU95 [s] G.G. Atanasova et al.fluid is in rotation).The decay rate coefficient, b, the corresponding Ekman decay timescale, T d (i.e.1/b), and spin-up completion times, t SU95 and T SU95 , are reported in Table2for the different combinations of Re and Fr considered.The Ekman decay timescale, T d , and spin-up completion times are no longer seen to increase linearly with Fr, and longer spin-up Ekman times were measured than those expected from the linear correlation proposed by

Table 3
Characteristic spin-down time for different Re and Fr.

Table 4
Summary of spin-down from solid body rotation for different Re.