Controlling secondary flow in Taylor-Couette turbulence through spanwise-varying roughness

Highly turbulent Taylor-Couette flow with spanwise-varying roughness is investigated experimentally and numerically (direct numerical simulations (DNS) with an immersed boundary method (IBM)) to determine the effects of the spacing and axial width $s$ of the spanwise varying roughness on the total drag and {on} the flow structures. We apply sandgrain roughness, in the form of alternating {rough and smooth} bands to the inner cylinder. Numerically, the Taylor number is $\mathcal{O}(10^9)$ and the roughness width is varied between $0.47\leq \tilde{s}=s/d \leq 1.23$, where $d$ is the gap width. Experimentally, we explore $\text{Ta}=\mathcal{O}(10^{12})$ and $0.61\leq \tilde s \leq 3.74$. For both approaches the radius ratio is fixed at $\eta=r_i/r_o = 0.716$, with $r_i$ and $r_o$ the radius of the inner and outer cylinder respectively. We present how the global transport properties and the local flow structures depend on the boundary conditions set by the roughness spacing $\tilde{s}$. Both numerically and experimentally, we find a maximum in the angular momentum transport as function of $\tilde s$. This can be atributed to the re-arrangement of the large-scale structures triggered by the presence of the rough stripes, leading to correspondingly large-scale turbulent vortices.


Introduction
In nearly all industrial applications and geophysical flows, turbulence is partly or completely wall bounded. In general, boundaries are not smooth but their surface is rather irregular and rough. Accordingly, such flows are extensively studied, although mainly under the approximation that the roughness is homogeneous (Jiménez 2004). Homogeneously rough surfaces have a characteristic length scale k that is much smaller than the largest wall-normal length scale δ. The effects of the roughness in these flows is believed to be confined to the immediate vicinity of the wall (i.e. the roughness sublayer), whereas in the outer, inertial layer, the flow only experiences the effective shear stress of the surface (Townsend's outer layer similarity, Townsend 1976). As such, the focus of many studies is to find functional relationships between the parameters describing the roughness geometry and the skin friction coefficient C f (Flack & Schultz 2010). In practice, however, flows are bounded by rough boundaries that not only vary on the scale of k, but also on a much larger scale s, with s = O(δ). Whereas these variations can occur either laterally (spanwise) or longitudinally (streamwise), we focus here only on the former. Such examples are found in shipping (i.e. the formation of stripes of bio-fouling on ship hulls (Schultz 2007)) and geophysical flows (e.g. the atmospheric flows over a spanwise-varying terrain (Ren & Wu 2011)).
Hitherto, the research is focused on the effects of spanwise-varying rough surfaces on canonical systems of wall-bounded turbulence, i.e. pipe (Koeltzsch, Dinkelacker & Grundmann 2002), boundary layer (Anderson et al. 2015) and channel flow (Chung, Monty & Hutchins 2018). The hallmark of flows over these surfaces is the presence of spanwise wall-normal secondary flows of size O(δ), with mean streamwise vorticity. Examples of studies where this has been observed are the works of Koeltzsch et al. (2002) on the effects of convergent and divergent grooves (reminiscent of shark skin) and the work by Wang & Cheng (2006) on spanwise-varying riverbeds. We note that earlier research dates back to the works of Hinze (1967Hinze ( , 1973 in the field of surface stress variations in duct flows. Following up on the work of Koeltzsch et al. (2002), Nugroho, Hutchins & Monty (2013) set out to perform a parametric study of the converging-diverging riblet surface in a zero pressure gradient boundary layer (BL). They found a thickening of the BL height above the converging regions, and vice versa above the diverging regions. Furthermore, the energy spectra showed an increased energy content of the larger scales.  performed stereo particle image velocimetry (PIV) in the spanwise wall-normal plane for the flow over a turbine blade replica and found spanwise variations of the order of δ in the mean velocity field. Within the same configuration, Mejia-Alvarez & Christensen (2013) identified regions of low momentum pathways (LMPs) and high momentum pathways (HMPs) in the instantaneous fields. Here, LMPs coincide with regions of enhanced turbulent kinetic energy (TKE) and Reynolds shear stress (RSS), and, rather remarkably, these regions do seem to occur at recessed roughness heights. Willingham et al. (2014) found very similar behaviour of the secondary flows for a much more regular surface geometry. Vanderwel & Ganapathisubramani (2015) found that only when s/δ 0.5 a secondary flow formation is observed. Here, s is the spacing between streamwise aligned Lego blocks. For s/δ 0.5 the secondary flows are confined to the roughness sublayer. Interestingly, contrary to the findings of Mejia-Alvarez & Christensen (2013), they found LMPs on top of their elevated blocks, and HMPs in between the roughness strips. Yang & Anderson (2017), however, found s/H 0.2, with H the channel Spanwise-varying roughness controls secondary flow in TC flow 883 A15-3 half-height, as the threshold for heterogeneous behaviour of the streamwise aligned pyramid elements.
It is an open question as to why the varying experiments in the literature show conflicting orientations of the secondary flow. We note that there are some differences, namely that the experiments in the paper by Vanderwel & Ganapathisubramani (2015) consist of small strips of elevated roughness, and wide regions of recessed roughness heights rather than alternating rough and smooth strips of equal width. One hypothesis concerns the location of the virtual origin of the rough part with respect to the smooth part. This could be either below (recessed) or above (protruding) the smooth surface, resulting in either an inflow (recessed) or an outflow (protruding) region above the roughness.
By carefully assessing the terms in the transport equation TKE, Anderson et al. (2015) found that spanwise variations of roughness lead to a local imbalance of production and dissipation of TKE, as already proposed by Hinze (1967). Since the secondary flows are driven by a spatial gradient of the RSS, they find that the mean secondary flows are Prandtl's secondary flow of the second kind (Bradshaw 1987). Medjnoun, Vanderwel & Ganapathisubramani (2018) observed a breakdown of outer layer similarity in the local profiles of the mean flow, turbulent intensity and the energy spectra, evidently induced by the presence of the secondary vortices. Finally, Chung et al. (2018) studied the influence of the spacing of idealized (i.e. no geometric induced disturbances to the flow) regions of low shear stress and high shear stress. They found that for s/δ 0.39 the notion of outer layer similarity is retained. Interestingly, for s/δ 6.28, they found a sign reversal of the isovels (stream velocity contour lines), with respect to the orientation of the secondary flows, that remain upwelling over low shear stress regions.
The aforementioned studies were all carried out in systems that lack two features which are intrinsic to many applications, namely the curvature in the streamwise direction (as in turbine blades), and the presence of strong secondary motions (as in the atmospheric boundary layer). Especially, the application for turbine blades is interesting, as also studied by . Here we come closer to this application by studying the turbulent flow over curved, rotating walls, present in the turbines.
1.1. Taylor-Couette flow A canonical system in which these two properties can be observed simultaneously is the Taylor-Couette (TC) flow. TC flow is the flow between two coaxially, independently rotating cylinders. Its geometry is characterized by two dimensionless parameters; the radius ratio η = r i /r o and the aspect ratio Γ = L/d. Here, r i and r o are the inner cylinder radius and outer cylinder radius respectively, L is the axial length of the inner cylinder and d = r o − r i the gap between the two cylinders. Since TC is a closed system, one can directly relate global and local quantities through exact mathematical relations (Eckhardt, Grossmann & Lohse 2007). The driving strength in TC flow is expressed in dimensionless form by the Taylor number where ω i and ω o are the inner and outer angular velocities of the cylinders, respectively, ν is the kinematic viscosity of the fluid andσ = ((1 + η)/(2 √ η)) 4 883 A15-4

D. Bakhuis and others
is the so-called geometric Prandtl number, in analogy to the Prandtl number in Rayleigh-Bénard convection (Eckhardt et al. 2007). Alternatively, when the outer cylinder is at rest (ω o = 0), the driving strength can also be expressed with a Reynolds number based on the inner scales Re i = r i ω i d/ν. This Reynolds number and Ta (ω o = 0), are related by Re i = (8η 2 /(1 + η) 3 ) √ Ta. In TC flow, the angular velocity flux J ω is radially conserved (Eckhardt et al. 2007), where the effects of the top and bottom plates are assumed negligible. Here, J ω = r 3 ( u r ω A,t − ν(∂/∂r) ω A,t ), where the brackets · A,t denote averaging over a cylindrical surface and time. The angular momentum flux for the case of laminar flow is J ω In this way, the response of the flow is quantified by the dimensionless Nusselt number (Nu ω ), which is also directly related to the torque T that is required to drive the cylinders at constant speed, i.e. (1.2) Here, ρ is the density of the working fluid. Alternatively, the torque of the system can be non-dimensionalized to form the friction coefficient C f = T /(ρLν 2 Re 2 i ), which is directly related to the Nusselt number, The inner friction velocity u τ ,i is also related to the torque by u τ ,i = T /(2πr 2 i ρL), which is used to scale quantities in the inner layer, indicated with the superscript '+'. The friction angular velocity ω τ equals u τ /r. Lastly, a frictional Reynolds number based on the inner scales can be defined as Re τ = u τ ,i d/(2ν).
Secondary flows are featured in TC flow, in the form of large-scale vortices with a mean streamwise vorticity component, the so-called turbulent Taylor vortices (TTVs). These structures are reminiscent of laminar Taylor vortices, which transition through a series of instabilities into turbulence once the flow becomes unstable (Taylor 1923). As noted by Chouippe et al. (2014), the axial wavelength λ/d of the TTVs, i.e. the distance between two rolls, is primarily a function of η and Re. When Re is large enough (O(10 6 )), the rolls are observed to persist in the system (Huisman et al. 2014). Here, multiple states for η = 0.716 can be observed in a certain regime of counterrotating cylinders, namely a ∈ [0.17, 0.51], where a = −ω o /ω i is their rotation ratio. These multiple states are characterized by a change in the number of rolls present in the system and, as a consequence, in their averaged axial wavelength (λ/d = 1.46 or λ/d = 1.96). These states, with the transition between them being strongly hysteretic, even at Re = O(10 6 ) (Huisman et al. 2014;van der Veen et al. 2016;Gul, Elsinga & Westerweel 2017), result in different torques for the same rotation rates, which reflects the importance of the large-scale structures (TTV) in transporting angular momentum. At pure inner cylinder rotation however (a = 0), no multiple states are found and the rolls are observed to be less coherent and stable. Finally, we note that the effect of the curvature of the cylinders is quantified by the radius ratio η, and it has a tremendous impact on the flow organization, as was reported by Ostilla-Mónico et al. (2014a,b). For a detailed review on turbulent Taylor-Couette flow we refer the reader to the review of Grossmann, Lohse & Sun (2016).
Roughness in a TC geometry has been studied in various ways: Cadot et al. (1997) and van den Berg et al. (2003) used obstacle roughness, in the form of axial riblets, to study the scaling of the angular momentum transport with the driving strength.
Spanwise-varying roughness controls secondary flow in TC flow 883 A15-5 Zhu et al. (2016) investigated the influence of grooves for large Ta (O(10 10 )), and found that at the tips of the grooves, plumes are preferentially ejected. In a more recent work, Zhu et al. (2018) found that, by using a similar configuration of rough walls as in van den Berg et al. (2003), the scaling Nu ω ∝ Ta 1/2 , where the logarithmic corrections of the ultimate regime vanish for asymptotically high Reynolds numbers, is obtained. This regime is referred to as the asymptotic ultimate regime (Grossmann et al. 2016;Zhu et al. 2018). They attribute this to a dominance of the pressure drag over the viscous drag on the cylinders. Very recently, Berghout et al. (2019) studied the influence of sandgrain roughness in TC flow, and found similarity of the roughness function with the same type of roughness in pipe flow (Nikuradse 1933). None of the TC papers described above reported an influence of the roughness variations in the axial direction, i.e. the spanwise direction.
In this paper we will fill this gap and study the effects of spanwise-varying roughness in highly turbulent TC flow with Ta up to O(10 12 ), for the case of pure inner cylinder rotation a = 0, where secondary flows are present in the form of TTVs. In particular, we focus on the effect of spanwise-varying roughness on the TTVs and, thus, on the global and local response of the flow. We introduce the roughness through a series of stripes which extend along the entire circumference of the inner cylinder (IC). This gives rise to a spanwise (axial) arrangement of roughness which we characterize with the widths of the roughness stripe. We conduct both experiments and direct numerical simulations (DNS) for various dimensionless stripe widthss = s/d, i.e. the width of the roughness stripe normalized with the gap width. For practical reasons, the roughness area is slightly larger than the smooth area, with a rough area coverage of 56 %.
The structure of the paper is as follows. In § 2, we introduce the experimental and numerical methods. In § 3.1 we show the local response of the flow due to the varying roughness arrangement. In § 3.2, we study its effect on the global quantities. In § 3.3, we link the global and local observations and explain the physical mechanism between the interaction of the rolls and the roughness. We finalize the paper in § 4 with some conclusions and an outlook to future work.

Numerical methods
The radius ratio η = 0.714, the same as in the experiments. The streamwise and spanwise lengths of the computational domain are set to match the minimum computational domain size as studied in Ostilla-Mónico et al. (2014c). As such, we simulate aspect ratios of 2.0 Γ 3.32 with periodic boundary conditions in the axial direction. The mean roughness height h r /d = 0.02 is small, so that the effective radius ratio is only slightly affected. We scale the roughness stripe such that the maximum roughness height, and thus the maximum blockage ratio, was max(h r ) = 0.1d. Depending ons, we cut out a portion of roughness from the scanned surface. The roughness was then mirrored and concatenated to obtain a streamwise homogeneous spanwise periodic stripe.
The sandpaper roughness was implemented in the code by an immersed boundary method (IBM) (Fadlun et al. 2000). In the IBM, the boundary conditions were enforced by adding a body force f to the Navier-Stokes equations. A regular, non-body fitting, mesh can thus be used, even though the rough boundary has a very complex geometry. We perform interpolation in the spatial direction preferential to the normal surface vector to transfer the boundary conditions to the momentum equations. The IBM has been validated previously (Fadlun et al. 2000;Iaccarino & Verzicco 2003;Stringano, Pascazio & Verzicco 2006;Zhu et al. 2016). A moving average over 10 × 10 points is employed to smooth the scan from measurement noise. Finally, we set the resolution based on the demands ( z + , r + i θ < 3), which is small enough to recover the smallest geometrical features of the surface.

Experimental apparatus with spanwise roughness
The experiments have been performed in the Twente Turbulent Taylor-Couette (T 3 C) facility as shown in figure 1(a) (details of the experimental facility can be found in van Gils et al. (2011)). The inner cylinder has a radius r i = 200 mm and the outer cylinder has a radius r o = 279.4 mm, such that the gap size is d = r o − r i = 79.4 mm, and the radius ratio η = r i /r o = 0.716. To minimize the effects of the endplate on the torque signal, the inner cylinder is partitioned into three sections and the torque is only measured at the middle section (L mid = 536 mm), shown as the highlighted section in figure 1a). The combined length of the cylinders is L = 927 mm, which leads to an aspect ratio Γ = L/d = 11.7. The outer cylinder (OC), is made of transparent acrylic which allows for optical access to the flow. The working fluid is demineralised water. We apply spanwise-varying roughness to the inner cylinder, which leads to patterns of homogeneously rough and hydrodynamically smooth bands in the spanwise direction (see figure 1a). The rough stripes are made of P36 ceramic industrial grade sandpaper and are fixed to the IC using double-sided adhesive tape (≈1 mm). In figure 2, we show the height scan of a roughness element using confocal microscopy. The scan revealed that the height (h r ) of the roughness is mostly within ±2σ (h r ) of the mean, giving a characteristic length scale k ≡ 4σ (h r ) = 695 µm (see figure 2b). More statistics of the roughness are shown in table 1. We fix the surface coverage  of the roughness at 56 % such that 0.56A i of the cylinder is rough, where A i = 2πr i L is the area of the entire IC. The middle section of the inner cylinder has a constant coverage of 56 %.

Global measurements: torque
We measured the torque T needed to drive the inner cylinder at constant angular velocity (the outer cylinder is kept at rest). For this we used a hollow flange reaction torque transducer connecting the drive shaft and the inner cylinder, sampled at 100 Hz. 883 A15-8 s = s/d is the normalized roughness width with smooth and rough indicating the fully smooth and fully rough cases respectively (i.e. 0 % or 100 % coverage of sandpaper roughness). N θ × N z × N r is the numerical resolution in the streamwise, spanwise and radial directions, respectively. Γ = L/d is the aspect ratio and L θ = r i 1 3 π is the constant streamwise length of the domain. ∆(ω) + is the downward shift of the angular velocity profile ω + . r + min is the minimum spacing in the wall-normal direction at the location of the maximum roughness height. r + max is the maximum spacing in the wall-normal direction. r + i θ = z + ≈ 2.7 (r + o θ ≈ 3.8) is the grid spacing in the streamwise and spanwise directions. In the DNS, the roughness height k + = 4σ (h r ) + = 130 ± 1 for all rough cases. t av /T is the averaging time needed to collect statistics, normalized with the bulk flow time scale T = d/(r i (ω i − ω o )). All experimental Nu ω are based on global torque measurements.

D. Bakhuis and others
We continuously measured the torque while quasi-statically ramping the frequency of the inner cylinder, f i , from 5 to 18 Hz, with increments of 0.002 Hz s −1 . This corresponds to Ta between 4 × 10 11 and Ta ≈ 6 × 10 12 . The system is temperature controlled such that all the experiments were performed at 21 ± 1 • C and all fluid flow properties are calculated at their actual temperature. Table 2 shows additional experimental parameters.

Local measurements: LDA and PIV
We performed a spanwise scan of the streamwise velocity with laser Doppler anemometry (LDA). The scan was performed at the middle of the gap,r = (r − r i )/d = 0.5, at a fixed Ta = 9.5 × 10 11 . Due to the axial symmetry, we scan the flow from mid-height to the bottom of the system. The flow was seeded using 5 µm-diameter polyamide particles with a density of 1030 kg m −3 that act as passive tracers . The laser beams went through the outer cylinder and were focused in the middle of the gap. We corrected for curvature effects (and refractive index effects) by numerically ray tracing the LDA beams, as was shown in Huisman, van Gils & Sun (2012). The spanwise extent of the LDA scans was 0 z/L 0.5. Particle image velocimetry measurements were performed at Ta = 9.5 × 10 11 (same as LDA) in the radial-azimuthal plane. The scan was performed for various heights and for alls. For the PIV measurements the flow was seeded with different particles that act as passive tracers, namely, fluorescent polymer particles (Dantec FPP-RhB-10) with diameters of 1-20 µm with a seeding density of ≈0.01 particles pixel −1 . These particles have an emission peak with a wavelength of ≈565 nm. We illuminated the particles with a Quantel Evergreen 145 532 nm, double pulsed laser. A cylindrical lens was used to create a light sheet of ≈1 mm thickness. The height of the laser sheet is altered using a custom traverse system. The images were captured with an Imager SCMOS (2560 × 2160 pixel) 16 bit camera with a Carl Zeiss 2.0/100 lens through a window in the bottom plate of the set-up and a mirror. The camera was operated in double frame mode with a frame rate f = 15 Hz which was much smaller than the inverse interframe time 1/ t, i.e. t 1/f . In order to enhance the particle contrast in the images, we added an Edmund High-Performance Longpass 550 nm filter to the camera lens. For everys, the spanwise extent of the experiments was different. This is done because -as will be shown later -the aspect ratios of the rolls change depending ons. For the smallests = 0.63 however, the spanwise resolution was δz/L ≈ 0.011 while for the largest values = 3.74, δz/L ≈ 0.022. Since we scan in the spanwise direction, the focus of the camera was changed accordingly. For each height, 1000 image pairs were recorded. The fields were resolved with a commercial PIV software (Davis 8.0) based on a multi-step method. The initial interrogation window size was set to 64 × 64 pixels and it decreased to 32 × 32 pixels for the last iteration, for increased accuracy. The fields are calculated in Cartesian coordinates, which we transformed to polar coordinates. The final result were the fields in the form u = u r (r, θ, t)ê r + u θ (r, θ, t)ê θ , where u r and u θ are the radial and streamwise velocity components which depend on the radius r, the azimuthal (streamwise) direction θ and time t. For an example of a typical azimuthal velocity field obtained from the experiments, see figure 3.

Laser Doppler anemometry
In order to get a first insight into the effect of the roughness on the flow, we performed spanwise scans of the streamwise flow velocity at mid-gap using LDA. Subsequently, we calculated the standard deviation of the streamwise velocity. The experiments will be compared with the simulations. Figure 4 shows an example of the snapshot of the instantaneous angular velocity from the simulations. In figure 5(a), we show the standard deviation of the streamwise velocity σ (u θ ) normalized by the inner cylinder velocity u i , as a function of the height, for variouss. Here, the spanwise coordinate is normalized using the cylinder gap width d such thatz = z/d. In the top row, we compare the results from direct numerical simulations and experiments. Since we employ different boundary conditions in the spanwise direction in DNS (periodic) and experiments (stationary solid plates), the comparison will be most reliable at midheight, where end effects are negligible. Therefore, we overlay the results from DNS (orange lines) over the results from experiments (blue lines) at matching locations with respect to the roughness stripes, atz ∈ (1, 6). In the bottom row we only assess the results from DNS and consequently plotz ∈ (0, Γ ). Figure 5(a) reveals that, for the case of the largest stripe width (s = 3.74), the smooth section has, on average, a value of σ (u θ )/u i ≈ 0.03, slightly larger than we Hatched regions indicate spanwise translated copies of the same data (including averages) -possible due to the periodic boundary condition in the axial direction -to allow for straightforward comparison. Imprints of the large secondary flows on the friction at the cylinder walls is observed, where impacting region experience a higher shear stress.
trend can be seen at the lower roughness section (z ≈ 0.93) of this case. However, this might be influenced by the lower bottom plate of the system. When looking at thes = 1.87 case, we see very similar, however more pronounced, dynamics. Streamwise velocity fluctuations are promoted in regions where the roughness is present, as suggested by the appearance of local peaks centred at the position of the rough stripes. This effect is further seen for the casess = 1.23 ands = 0.93, where we observe similar profiles. At their smooth areas however, we also observe enhanced velocity fluctuations, albeit less pronounced than the locations above the rough patches. This effect is not seen fors > 1.27. For the final case withs = 0.61 these trends seems to fade away and we see that σ (u θ ) becomes more spanwise independent, i.e. the peaks are less pronounced, and do not seem to follow the topology of the roughness stripes.
The results of the DNS, presented together with the experiments in figure 5(a) exhibit very similar behaviour. When normalized with u i , we find that the standard deviation of the streamwise velocity shows similar values as in the experiments. This is intriguing since the k/d values in the simulations are almost one order of magnitude larger than in the experiments. Above the rough stripes, we find enhanced σ (u θ )/u i , whereas over smooth stripes we find diminished σ (u θ )/u i . Fors = 0.61 however, the trends are somewhat different. There, we find enhanced σ (u θ )/u i above the smooth region in the DNS, which is explained by the recombination of plume ejection regions, forming one larger TTV above the smooth regions, see figure 7.
The findings presented in figure 5(a) show that the presence of the roughness affects the relative turbulence statistics in the bulk of the flow, far away from the roughness sublayer region (in contrast to homogeneous roughness in TC flow (Berghout et al. 2019)), reminiscent to what is found in studies of pipe and channel flows (Koeltzsch et al. 2002;Chung et al. 2018). Figure 5(b) shows the spanwise variations of the friction factor c f (z) (see § 1) on both the inner cylinder and the outer cylinder. The solid and dashed black lines represent c f (z) measured on the inner cylinder, and the solid (green) lines represent c f (z) on the outer cylinder. We average both in time and in the streamwise direction. Significant variations in c f (z) are observed which are linked to the orientation of the TTV. For impacting regions, i.e. plumes impacting on the inner cylinder boundary layer, the wall shear stress in the streamwise direction is enhanced. When plumes are ejected from the inner cylinder boundary layer, also known as the ejecting regions, the shear stress is reduced. This, once more, illustrates the relative strength of the secondary flow (TTV) and the mean flow. Fors = 0.93, the variations of the friction factor on the outer cylinder are even more pronounced, thus indicating that the strengths of the TTV are enhanced with enforcing the spanwise variations in the roughness. Fors = 0.47, the variations are not visible and the TTVs are severely weakened, but still present, see figure 7.

Particle imaging velocimetry and direct numerical simulations
To gain more insight into how the roughness alters the flow, we set out to measure the velocity field in the meridional plane using PIV at multiple heights. For the cases that required the highest axial resolution (s = 0.61 ands = 0.93) of the PIV measurements we employed δ z /d ≈ 0.13, where δ z is the spacing of laser sheets in z, so that every TTV roll pair was at least sampled 7 times. This is enough to resolve the wavelength of the structure. For all others, the sample resolution is higher than 7 times per wavelength. In figure 6, we show the temporal and streamwise averaged radial velocity component u r , normalized with u i , in the spanwise-wall-normal plane (z −r), where the radial coordinate is normalized such thatr = (r − r i )/d. Figure 6 shows that for the case ofs = 3.74, a very large structure can be seen, which consists of a large outflow region (positive u r ) aroundz = 5.84, while a large inflow region (negative u r ) is detected aroundz = 2.10. The situation is more pronounced for the cases ofs = 1.87,s = 1.23 ands = 0.93, where a clear roll-like structure (i.e. the TTV) can be observed. Note that the radial component in the flow changes sign along the spanwise direction as it should in the presence of a TTV. What is rather Spanwise-varying roughness controls secondary flow in TC flow 883 A15-13 remarkable is that the wavelength of the rolls λ changes for different values ofs. For the large structure ats = 3.74, the normalized wavelength isλ = λ/d ≈ 4.01. Ass decreases tos = 1.87,λ ≈ 1.49. Ats = 1.23, the wavelength decreases to a value ofλ ≈ 1.42. Ats = 0.93,λ ≈ 0.94, and finally for the smallest value ofs = 0.61, the wavelength increases slightly toλ = 1.09. We remind the reader of the work of Huisman et al. (2014), who revealed that for counter-rotation (−ω o /ω i ≈ 0.4), the average wavelength of the rolls could be eitherλ = 1.46 orλ = 1.96 depending on the state attained by the system. The current work shows that, by an appropriate choice ofs, the wavelength of the rolls can firstly, abandon its natural state, and secondly, it can be tuned within the rangeλ ∈ [0.94, 4.01]. The wavelengths described above were calculated by measuring the locations of two consecutive maxima and minima of u r t,θ,r bulk along z which are closest to mid-height. Here, the symbol · t,θ,r bulk denotes average over time, the streamwise direction and the bulk region, i.e.
In addition, we observe that outflow regions (flow away from the IC) are created in spanwise regions where the roughness is located; and conversely, inflow regions (flow towards the IC) are created in the smooth areas. Note that this orientation of the secondary flows is opposite to what is found in other canonical systems (e.g. pipe flow and channel flow (Willingham et al. 2014;Yang & Anderson 2017;Chung et al. 2018)), where one finds inflow regions above the rough stripes and outflow regions above the smooth stripes. Consistent with findings in boundary layers however, is the correlation between low momentum pathways and outflow regions (Vanderwel & Ganapathisubramani 2015). Indeed we find that LMPs -with respect to the IC and located on top of the roughness -are associated with outflow regions, similar to  Another interesting observation is that, since the driving is now from the BL rather than the bulk, the strength of the rolls changes depending on the value ofs, as evidenced by the magnitude of |u r |. In order to explore this feature in more detail, we quantify the strength of the rolls byũ r ≡ (u r /u i ) 2 t,θ,r bulk ,z λ as a function ofs. Here, the symbol · t,θ ,r bulk ,z λ denotes an average over time, the streamwise direction, the bulk region and the spanwise region that defines the wavelength of a single roll z λ . In figure 8(c), we showũ r as a function ofs, where we observe that the strength of the rolls increases with decreasings fors ∈ [0.93, 3.74]. However ats = 0.61 the trend is broken, and we observe thatũ r decreases with respect to the case of s = 0.93.
In order to obtain more insight into the mechanism(s) that lead to the variation of λ withs, we turn to DNS, albeit at a much lower Ta (≈1.0 × 10 9 ), and much higher roughness height (k/d ≈ 0.1). DNS at this high Re are very computationally expensive and we choose to simulate Γ 2 for which Ostilla-Mónico, Verzicco & Lohse (2015) showed that, at least for smooth walls, the aspect ratio is sufficient. Since very larges cases are not feasible for DNS, we focus on matching the exacts in the lower range. We will show that, despite the O(10 3 ) difference in Ta, the same observations found in the numerics are also found in the experiments.
First, we look at the streamwise velocity component. In figure 7, we plot the difference of the temporal and streamwise average of the angular velocity ω t,θ with respect to the temporal, streamwise and spanwise average of the angular velocity ω t,θ,z . This is done to emphasize the underlying organization of the TTVs. Here, we clearly observe that for alls, ejecting regions of angular velocity are originated in the rough stripes, similar to the preferential plume ejection sides at the tips of grooves in Zhu et al. (2016). These ejecting regions advect fluid from the roughness stripe on or at the inner cylinder towards the outer cylinder. These ejecting regions occur at each roughness stripe and as a consequence, an array of plume-like structures are formed along the spanwise direction. In TC flow (without roughness), plume-like structures are clear signatures of the presence of TTVs (Ostilla-Mónico et al. 2014b,c). A closer inspection of figure 7 reveals that for the largest value ofs = 1.23, the plumes have enough separation not to interact among them. Whens is decreased tos = 0.93, we observe that the plumes come closer and can, in fact, interact with each other. At the lowers = 0.61 however, the situation is rather different. The rough patches are too close to each other to create individual ejecting regions and therefore, merge to a larger collective plume. For thes = 0.61 case, we observe in both, experiments and DNS, that two ejecting regions, each located on top of a rough region, merge into a combined ejecting region. Whens is decreased to a value of 0.47, even three ejecting regions are combined to a single large plume. These combined plumes drive the TTV with a larger wavelength as compared to thes = 0.93 case, therefore decreasing the effective momentum transfer. The mechanisms leading to this optimum in transport remain elusive. We hypothesize that very lowλ constrains the roll structure too much, and ejecting and impacting plumes hinder each other, limiting the amount of transport. For highλ, the density of ejecting and impacting zones is too low, and the amount of transport is also limited.
The LDA, PIV and DNS explored in this section reveal that there is a mean effect of the spanwise-varying roughness on the large-scale secondary flows that exist in turbulent TC flow. We have seen thus far that the roughness pins the rolls, and that their wavelength and strength can be tuned depending of the choice ofs for a variety of Ta, and for a variety of roughness heights h. However, how does the flow respond globally, i.e. the angular momentum transport, to this change in morphology? This will be addressed in the following section.

Global response
We have performed torque measurements at an acquisition rate of 100 Hz at identical conditions, only varying the acceleration (ramping time). All measurements accelerated from 0 to 20 Hz at run times that varied between 30 minutes and 10 hours. The measured torque signal was binned and compared between various runs. If the Spanwise-varying roughness controls secondary flow in TC flow 883 A15-17 difference between the curves is smaller than the measurement error of the sensor (0.25 %), we call the runs indistinguishable. This is already the case for runs of 60 minutes.
The global response of the TC system can be expressed with Nu ω (1.2) or with the related friction coefficient C f . In figure 8(a), we show the compensated Nu ω as a function of the driving strength Ta, where a scaling of Nu ω ∝ Ta α , with α = 0.45 is observed for all thes explored. The shaded area indicates the error based on the standard deviation from three repeated measurements. The same data are represented as C f versus Re i in figure 8(d). In this figure, we also give best fits of the Prandtl friction law, defined as 1/ C f = a log 10 (Re i C f ) + b, for the fully smooth and fully rough cases. In the absence of roughness and within the same range of Ta, the scaling is found to be effectively Nu ω ∝ Ta 0.39 (Lathrop, Fineberg & Swinney 1992;van Gils et al. 2011;Paoletti & Lathrop 2011;Brauckmann & Eckhardt 2013;Huisman et al. 2014). In contrast, when both of the solid walls are made uniformly rough (i.e. pressure drag dominates), the scaling asymptotes to the ultimate regime predicted by Kraichnan, i.e. Nu ω ∝ Ta 0.5 (Kraichnan 1962;Zhu et al. 2018). In Zhu et al. (2018), the closest configuration to our study is the case of rough IC and smooth OC, for which an effective exponent α = 0.43 was found. We note that this exponent is slightly smaller than the ones observed in the current study. The reason behind this is currently unknown. We notice, however, that the roughness type in our study is rather different. In this study we use spanwise-varying sand grain roughness, while the roughness in Zhu et al. (2018) is made of rib obstacles and is oriented perpendicular to the streamwise direction.
In order to connect the observed dynamics of the TTVs with the global response, we plot in figure 8(b) the compensated Nusselt number Nu ω Ta −0.45 as a function of s for both the experiments and the numerics. We note that the exponent found for s = 1.87 (α = 0.44) is slightly smaller than α = 0.45. We notice rather remarkably, the appearance of a maximum arounds ≈ 0.93 for the experiments, ands = 0.61 for the DNS (shown in the inset of figure 8(b)). We attribute the appearance of this peak to the strengthening of the TTVs, which is caused by the variation ofs, and thus of λ. Explicitly, by lowerings, we can decrease the wavelength of the rolls, as seen in the inset of figure 8(c), thereby, bringing them closer together (see also § 3.1). As a consequence, the rolls are strengthened, which leads to an enhancement of the angular momentum transport; and thus, the peak arounds = 0.93. This peak is also visible in the normalized root mean square of the radial velocity,ũ r , plotted in figure 8(c). With decreasings,ũ r increases until an optimum is reached arounds = 0.93. Below the optimum,ũ r drops drastically to much lower values. The PIV measurements are carried out over a domain 0.3 < (r − r i )/d < 0.7, so that the value ofũ r in the boundary layers is not included. However, the DNS give us the value ofũ r throughout the entire domain, such that we include the values in the boundary layers. Considering thatũ r peaks in the BL, we can understand the higher values ofũ r for the DNS. Regardless of this, we find similar behaviours in the DNS and the experiments.
The optimum in transport is also observed by Huisman et al. (2014), although the mechanism leading to optimum transport there is quite different. While the rolls in their study are enhanced by counter-rotating the OC; in our case, the rolls are strengthened by forcingλ below their natural wavelength due to the right choice of the size of the spanwise-varying roughnesss. This is also supported by the observation that the magnitude of the radial velocity shows a maximum arounds = 0.93, as shown in figure 8(c). We note, however, that the torque is not measured throughout the entire spanwise length of the cylinders L = 927 mm, but in a smaller section of length L mid = 536 mm. As a result, the large structure identified previously for the case ofs = 3.74 (λ = 4.01), does not fit entirely in the measurement section (see the first panel of figure 6). The measured torque will be higher if part of the ejecting region (with respect to the IC) is excluded from the measurement zone, since there we experience a lower c f , see figure 5(b). The measured torque will consequently be lower if part of the impacting region (with respect to the IC) is excluded from the measurement zone, since there we experience a higher c f , see figure 5(b).
We also note that, in the case of the numerics, the position of the maximum is different than in the experiments. We attribute this to a combination of two effects. On the one hand, the DNS is performed at a lower Ta, which has an effect on the natural wavelength of the rolls, as was shown by Chouippe et al. (2014), who show that, for similar values of η, the wavelength of the rolls can decrease with decreasing Ta. On the other hand, the spanwise domain of the DNS is bounded by Γ ∈ [2.08, 3.32], which gives rise to limited box sizes. Thus, whens is varied, the rolls could suffer from an additional constraint due to the limited spanwise domain. In addition to this discrepancy, we also note that the scaling in the range of Ta for which the DNS are performed (≈1.0 × 10 9 ), is not known a priori. In the absence of a better choice, we compensate the numerical data using the same exponent as in the experiments ( figure 8(b)). However, we note that this exponent might be different due to the 2 decades of separation in Ta between the numerics and experiments, as was also shown by Zhu et al. (2018). We would like to emphasize, however, that in spite of these discrepancies, a maximum in angular momentum transport is observed for a givens in both the experiments and the numerics, which is solely a consequence of the varying spanwise wavelength of the TTV, dictated by the spanwise-varying roughness.

Velocity profiles
Having discussed the dynamics of the TTVs and the corresponding global response, in terms of the dimensionless torque, we now set out to study the streamwise, angular, velocity profiles (rather than the azimuthal profiles, as discussed in Grossmann, Lohse & Sun (2014) and Berghout et al. (2019)). To allow for a straightforward comparison between the respective velocity profiles, we run the DNS at constant friction Reynolds number Re τ = 690 ± 10. The profiles are then temporally, streamwise and spanwise averaged ω + = ω t,θ ,z /ω τ . The profiles still exhibit a logarithmic region when averaged over the entire spanwise coordinate range. Figure 7 shows however that the TTVs in the flow, following the spanwise-varying roughness, do not exhibit any outer similarity. Deviations of the streamwise and temporal averages from the mean logarithmic profiles are found up to ω + ≈ 2. For turbulent flows over rough walls, the streamwise velocity profiles retain their logarithmic form. However, the hallmark effect of rough walls is a downwards shift of this region (for any drag increasing surface), which can also be understood as an increase of the skin friction factor C f (Hama 1954). Figure 9(a) shows the angular velocity profiles ω + as a function of (y + − h + m ), where h + m is the virtual origin and equals the meltdown (i.e. mean) height of the rough surface and y + = (r − r i )/δ ν . We choose the meltdown height of the roughness over the full inner cylinder as the virtual origin. In figure 9(b) we show the velocity shift versus the wall-normal distance. The inset gives a vertical cut at y + = Re τ . Evident that also in this representation there is a maximum in the velocity shift. The position of this maximum (s = 0.61) is the same as the one obtained from the angular momentum transport (see § 3.2).
In figures 9(c) and 9(d) we employ conditional averaging of the angular velocity profiles over the smooth · · · smooth and rough · · · rough spanwise locations. The wall   (a) Angular velocity ω + profiles in the reference frame of the IC versus the wall-normal distance y + − h + m for variouss, where h + m is the virtual origin and equals the meltdown (i.e. mean) height h m /k of the rough surface and y + = (r − r i )/δ ν . The solid black line represents the uniformly rough case. (b) Angular velocity shift ω + as a function of y + − h + m for varyings. In the inset of (b), we show the angular velocity shift ω + versus the wall-normal distance y + − h + m . Here, we observe a maximum downward shift ( ω + = ω + s − ω + r > 0, where subscripts s and r represent the smooth and rough profiles respectively) of the angular velocity profile for the simulation where we cover the entire inner cylinder with sandpaper roughness (i.e. uniformly rough). (c) The angular velocity conditioned on the spanwise location: above smooth surface ( ω + smooth ) and rough surface ( ω + rough ). (d) The angular velocity shifts conditioned on the spanwise location: above smooth surface ( ω + smooth ) and rough surface ( ω + rough ). We conclude that the plumes originating from the roughness elements lead to enhanced mixing of streamwise momentum, and hence a downwards shift of the velocity profiles. Further away from the wall, the bulk is well mixed and the streamwise profiles above smooth and rough wall locations converge to a similar value.
shear stress in the viscous normalization is taken over the entire spanwise height L of the inner cylinder. We calculate h m also for the rough and smooth patches (where h m,smooth = 0). We already deduce from figure 7 that significant variations in the temporal and streamwise average of the velocity profiles are expected, at least close to the roughness. Indeed we find that for alls the region above the roughness is better mixed, due to the presence of plume-like structures originating from the rough surface. The angular velocity profiles is thus shifted downwards in comparison to the average over the entire IC. For the smooth wall conditioned profiles, we observe the opposite, such that the profiles lay higher.
The merging of plumes from different rough patches into a large-scale coherent TTV is also observed in the cross-over of ω + smooth and ω + rough fors = 0.61 in 883 A15-20 D. Bakhuis and others figure 9(d) at y + − h + m ≈ 210. Further into the bulk flow, turbulent processes mix out the inhomogeneous effects of rough wall attached plumes, and the angular velocity profiles converge to similar values. However, we note that even at y + − h + m = Re τ , ω + smooth and ω + rough differ to ≈0.5.

Conclusions and outlook
In conclusion, we have investigated, both numerically and experimentally, large Taylor number Taylor-Couette flow in the presence of spanwise-varying roughness, which consists of an arrangement of stripes of widths, that covers the entire circumference of the inner cylinder. In the experiments, the stripes were made from sandpaper, while in the numerics a confocal microscopy scan of the surface was implemented by means of the immersed boundary method.
Remarkably, we have found that by varyings in the ranges = [0.47, 3.74] we can alter the spanwise wavelength of the turbulent Taylor vortices within the rangẽ λ ∈ [0.94, 4.01], even if the roughness height was very low (k/d ≈ 0.01). This manipulation was observed to hold at Ta = O(10 9 ) and Ta = O(10 12 ).
In the experiments, the scaling of the Nusselt number with the driving strength was found to be effectively Nu ω ∝ Ta 0.45 for Ta ∈ [5 × 10 11 , 5 × 10 12 ]). The experiments and DNS also revealed that inflow regions (u r < 0) originated between the rough stripes, where the inner cylinder was hydrodynamically smooth (in contrast to secondary flows induced by spanwise-varying roughness in channel flow, where the orientation of the vortices is reversed (Chung et al. 2018)). Conversely, at the centre of the rough stripes, we observed the creation of outflow regions (u r > 0) which were accompanied by the promotion of streamwise velocity fluctuations σ (u θ ) at mid-gap. At these spanwise locations (centre of rough stripes), we observed, in both the numerics and experiments, the emission of plume-like structures, which are responsible for the creation and pinning of the rolls. Since the coverage of the roughness was fixed, we showed that by reducings, we can effectively bring these structures closer, and enhance the interaction of the rolls, as evidenced by the increment in |u r |. As a consequence of this interaction, the flow responded globally by inducing a maximum of angular momentum transport ats = 0.93 in the experiments, ands = 0.61 in the numerics.
The numerical simulations and experiments in this paper fall into the category of low momentum pathways above the roughness and high momentum pathways above the smooth stripes. Note that these are seen in the reference frame of the rotating inner cylinder. These findings agree with the hypothesis that points to a virtual origin effect, since the roughness is protruding and outflow regions are located above the roughness.
For increasing Ta, the relative strength of the turbulent Taylor vortices with respect to the mean flow decreases. Therefore, we find a weaker imprint of the vortices on the angular velocity profiles for PIV measurements at Ta = O(10 12 ), than we do for the DNS at Ta = O(10 9 ), for smooth walls. Consequently, we also find that the effects of spanwise-varying roughness on the flow are somewhat stronger at Ta = O(10 9 ) than for Ta = O(10 12 ).
We wish to stress that in this study the change in the morphology of the largescale structures is only due to the spanwise-varying roughness (of very low height) and not by a change of Γ or η, which opens the possibility of exploring different configurations in which the rolls can be tuned at such large turbulence levels.
Many questions arise from the aforementioned observations. Understanding the mechanisms leading to the merging of plume ejection regions, and accompanying Spanwise-varying roughness controls secondary flow in TC flow 883 A15-21 parameter boundaries at which this occurs, would lead to a further insight into the dynamics of the TTVs. Furthermore, it would be intriguing, in the spirit of Bakhuis et al. (2018), to study the influence of spanwise-varying regions of idealized high and low wall shear stress, without geometrical induced disturbances. It is an open question whether one could also alter λ, without the interaction of the plumes.