Direct numerical simulations of Taylor--Couette turbulence: the effect of sand grain roughness

Progress in roughness research, mapping any given roughness geometry to its fluid dynamic behaviour, has been hampered by the lack of accurate and direct measurements of skin-friction drag, especially in open systems. The Taylor--Couette (TC) system has the benefit of being a closed system, but its potential for characterizing irregular, realistic, 3-D roughness has not been previously considered in depth. Here, we present direct numerical simulations (DNSs) of TC turbulence with sand grain roughness mounted on the inner cylinder. The model proposed by Scotti (\textit{Phys. Fluids}, vol. 18, 031701, 2006) has been improved to simulate a random rough surface of monodisperse sand grains, which is characterized by the equivalent sand grain height $k_s$. Taylor numbers range from $Ta = 1.0\times 10^7$(corresponding to $Re_\tau = 82$) to $Ta = 1.0\times 10^9$($Re_\tau = 635$). We focus on the influence of the roughness height $k_s^+$ in the transitionally rough regime, through simulations of TC with rough surfaces, ranging from $k_s^+=5$ up to $k_s^+ = 92$, where the superscript `$+$' indicates non-dimensionalization in viscous units. We find that the downwards shift of the logarithmic layer, due to transitionally rough sand grains exhibits remarkably similar behavior to that of the Nikuradse (\textit{VDI-Forschungsheft} 361, 1933) data of sand grain roughness in pipe flow, regardless of the Taylor number dependent constants of the logarithmic layer.


Introduction
Many turbulent flows in nature and industry are bounded by rough boundaries. Examples are the surface of the planet with respect to geophysical flows or fouling on ship hulls with respect to open waters. Such rough boundaries strongly influence the total drag, with often adverse consequences in the form of higher transport costs. Therefore, it becomes of paramount importance to understand the physics behind such changes in drag, ultimately leading to better informed management of the problem. One key recurring question concerns the influence of the roughness topology on the drag coefficient.
Seminal work by Nikuradse (1933) investigated the influence of the height of closely packed, monodisperse, sand grains in pipe flow. This work has become one of the pillars in the field. Later, a vast amount of research was carried out to study the influence of roughness on the canonical systems of turbulence -pipe, channel, and boundary layer flows -aiming for a better understanding of the roughness effects on turbulent flows Ligrani & Moffat (1986); Schultz & Flack (2007); Chan et al. (2015); Busse et al. (2017); also see Jimenez (2004), and Schultz & Flack (2010) for comprehensive reviews.
Next to pipe flow, Taylor-Couette (TC) flow -the flow between two coaxial, independently rotating cylinders -is another canonical system in turbulence . Closely related to its 'twin' of Rayleigh-Bénard (RB) turbulence (Busse 2012;Eckhardt et al. 2007), it serves as an ideal system to study the interaction between boundary layer and bulk flow. Very long spatial transients, as found in open systems, are bypassed by the circumferential restrictions. Since the domain is closed, global balances can be easily derived and monitored, giving room for extensive comparison between theory, experiments and simulations. Further, the streamwise curvature effects find many applications in industry. For these reasons, we set out to investigate the effects of roughness on the turbulent fluid flow in the TC system.
Over the last century, much work has been carried out with the aim of understanding the effect of the roughness topology on fluid flow. One of the consequences of roughness is the change of the wall drag, which can be expressed as a shift of the mean streamwise velocity profile ∆U + ≡ (U s − U r )/u τ , where ∆U + is known as the Hama roughness function (Hama 1954) and U s , U r are the smooth-wall and the rough-wall mean streamwise velocities, respectively. Clauser (1954) and Hama (1954) already observed that roughness effects are confined to the inner region of the boundary layer. Namely that the velocity deficit occurs in the log layer. Shockling et al. (2006) and Schultz & Flack (2007) later found strong support for this notion. The shift in the mean streamwise velocity profile U (y) then becomes (Pope 2000) u + (y + ) = 1 κ log y + + A − ∆u + (1.1) As usual, the superscript '+' indicates a scaling in viscous units (e.g. length y + = yu τ /ν and velocity u + = u/u τ ) and u τ is the friction velocity, u τ = τ w /ρ with τ w being the total tangential stress at the wall, and ρ the fluid density. Note that in this representation, the skin-friction coefficient C f is related to the friction velocity by C f = 2(u τ /U 0 ) 2 , where U 0 is the centerline velocity (Pope 2000). It has been found that for TC turbulence κ and A are not constant anymore, but are functions of the inner cylinder Reynolds number Re i at least until Re i = 10 6 (Huisman et al. 2013). Therefore, for TC we here suggest the generalization; u + (y + ) = f 1 (Re i ) log(y + ) + f 2 (Re i ) − ∆u + (1. 2) with f 1 (Re i ) and f 2 (Re i ) being unknown functions. The questions now are: i) How does ∆u + depend on the parameters that characterize the surface geometry. ii) Can ∆u + be generalized to other flows. Although many parameters influence the Hama roughness function ∆u + (Schlichting 1936;Musker 1980;Napoli et al. 2008;Chan et al. 2015;MacDonald et al. 2016), the most relevant parameter is the characteristic height of the roughness k + . In a regime in which the pressure forces dominate the drag force, any surface can be collapsed onto the Nikuradse data by rescaling the roughness height to the so-called 'equivalent sandgrain roughness height' k + s . Nikuradse (1933) found that three regimes of the characteristic roughness height k + s can be differentiated with respect to the effect of roughness (Flack et al. 2012). For k + s 5, the rough wall appears to be hydrodynamically smooth and the roughness function ∆u + goes to zero. For k + s 70, the wall drag scales quadratically with the fluid velocity and the friction factor C f is independent of the Reynolds number, indicating that hydrodynamic pressure drag (also called form drag) on the roughness dominates the total drag. The transitionally rough regime is in between these two regimes. Where in the fully rough regime, a surface is fully determined by k + s to give a collapse onto the fully rough asymptote (Schultz & Flack 2010), in the transitionally rough regime, any unique surface corresponds to a unique roughness function. This can be attributed to the delicate interplay between pressure drag, viscous drag, and the weakening of the so-called viscous generation cycle (Jimenez 2004).
An intriguing feature of the data from Nikuradse (1933) is at k + s ≈ 5, where roughness effects suddenly result in an inflectional increase of ∆u + , as compared to the gradual increase of the roughness function found by Colebrook (1939) who extracted an empirical relationship from many industrial surfaces (Shockling et al. 2006). Later, this inflectional behavior was also observed for tightly packed spheres (Ligrani & Moffat 1986), honed surfaces (Shockling et al. 2006), and grit-blasted surfaces (Thakkar et al. 2018). Chan-Braun et al. (2011) had too few points to find the inflectional behavior; however, their two simulations of monodisperse spheres in regular arrangement collapsed on the Nikuradse curve. In the Moody (Moody 1944) representation, this inflectional behavior manifests itself as a dip in the friction factor C f , leading to a significantly lower drag coefficient (≈ 20% (Bradshaw 2000)) in comparison to the monotonic behavior of Colebrook (1939), on which the original Moody diagram is based. Although it is proposed that the inflectional behavior has to do with the monodispersity of the roughness leading to a critical Reynolds number at which the elements become active (Jimenez 2004), recent simulations by Thakkar et al. (2018) for a polydisperse surface (containing irregularities with a range of sizes) also show this inflectional behavior. In a broader sense, the DNS by Thakkar et al. (2018) are interesting since they show, for the first time, a surface that very closely and quantitatively resembles the Nikuradse (1933) roughness function in all regimes, the authors found k + s = 0.87k rms . Regarding TC flow, only a few studies have looked at the effect of roughness (Cadot et al. 1997;van den Berg et al. 2003). Recently, the effect of regular roughness on TC turbulence has also been investigated by means of DNS (Zhu et al. 2016(Zhu et al. , 2017(Zhu et al. , 2018b. Zhu et al. (2016) looked at the effect of axisymmetric grooves on the torque scaling, boundary layer thickness, and plume ejections. They find that enhanced plume ejection from the roughness tips can lead to an ultimate torque scaling exponent of N u ∝ T a 0.5 , although for higher T a the exponent saturates back to the ultimate scaling effective exponent of 0.38. Zhu et al. (2017) then simulate transverse bar roughness elements on the inner cylinder to disentangle the separate effects of viscosity and pressure, and find that the ultimate torque scaling exponent of N u ∝ T a 0.5 is only possible when the pressure forces dominate at the rough boundary (Zhu et al. 2018b).
In contrast to the above mentioned previous work, in which the roughness consisted of well-defined transverse bars with constant distance and heights (Zhu et al. 2017(Zhu et al. , 2018b, in this research we set out to investigate the effect of irregular, monodisperse roughness, resembling the sand grain roughness reported by Nikuradse (1933). We model the roughness as randomly rotated and semi-randomly translated ellipsoids of constant Figure 1: Illustration of the Taylor-Couette setup. (a) Taylor-Couette setup with inner cylinder sand grain roughness. ω i is the inner cylinder angular rotation rate, r i is the inner cylinder radius, r o is the outer cylinder radius, and d is the gap width. (b) A more detailed and accurate zoom of the sand grain roughness that is modeled. The outer cylinder is stationary and smooth. volume and aspect ratio, based on the roughness model (subgrid-scale) of Scotti (2006). Previously, a fully resolved version of the model by Scotti (2006) was also used for largeeddy simulations in channel flow (Yuan & Piomelli 2014). Taylor numbers in our DNS range from T a = 1.0×10 7 (Re τ = 82) to T a = 1.0×10 9 (Re τ = 635); therefore, we capture both classical (laminar-like boundary layers) and ultimate (turbulent boundary layers) regimes (Ostilla-Mónico et al. 2014;Grossmann et al. 2016;Zhu et al. 2018b). Moreover, whereas previous research on roughness in TC flow focussed on the torque scaling, we now look at the effects of the roughness height on the Hama roughness function ∆U + in the transitionally rough and fully rough regimes, ranging from k + s = 5 to k + s = 92. This manuscript is structured as follows. In §2 & §3, we elaborate on the TC setup, the roughness model and the numerical procedure. In §4.1, we study the velocity profiles and present the effects of the roughness height on the Hama roughness function. In §4.2, we present the global response of the system. In §4.3 we study the flow structures. In §4.4 the fluctuations close to the roughness are studied and in §4.5 we present radial profiles of various other quantities. Finally, in section §5 we draw our conclusions and propose future research directions.

Taylor-Couette flow
The TC setup, as shown in figure 1, comprises independently co-or counter-rotating concentric cylinders around their vertical axis. The flow, driven by the shear on both of the cylinders, fills the gap between the cylinders. r i is the inner cylinder radius, r o is the outer cylinder radius, and the radius ratio is defined as η = r i /r o . For this research, we set η ≈ 0.714, in accordance with the experimental T 3 C setup (Huisman et al. 2013), and previous simulations (Zhu et al. 2017). Γ = L/d is the aspect ratio, where L is the height of the cylinders, and d = r o − r i = 0.4r i is the gap width. Here, Γ ≈ 2 such that one pair of Taylor vortices fits in the gap, and the mean flow statistics become independent of the aspect ratio (Ostilla-Mónico et al. 2015). In the azimuthal direction we employ a rotational symmetry of order 6 to save on computational expense such that the streamwise aspect ratio of our simulations becomes L θ /d = (r i 2π/6)/d = 2.62. Brauckmann & Eckhardt (2013) and Ostilla-Mónico et al. (2015) found that this reduction of the streamwise extent does not effect the mean flow statistics. This gives L θ /(d/2) = 5.24 and L/(d/2) ≈ 4.0.
To maintain convenient boundary conditions, we solve the Navier-Stokes (NS) equations in a reference frame rotating with the inner cylinder (ω i e z ) The NS equations in that reference frame become T a 1/2 = Re −1 where f (η) is the so-called 'geometric Prandtl' number (Eckhardt et al. 2007): here f (0.714) ≈ 1.23 and T a is the Taylor number, which is a measure of the driving strength of the system, Note that the pressure in the equations above represents the 'reduced pressure' that incorporates the centrifugal term;P = p − ω 2 i d 2r2 2r 2 i |ωi−ωo| 2 e r with p = ρr 2 i |ω i − ω o | 2 p and p is the physical pressure. It is directly clear that the centrifugal force in TC flow is analogous to the gravitational force in RB flow (Eckhardt et al. 2007). The final term on the right-hand side of eq. (2.1) and (2.2) gives the Coriolis force, with Ro −1 being the inverse Rossby number (2.7) Analogous to RB flow, the global response of TC flow can be expressed in terms of a Nusselt number. In the former, the Nusselt number expresses the dimensionless conserved heat flux, whereas in the latter the Nusselt number expresses the dimensionless conserved angular velocity flux J ω , calculated by: with the laminar flux given by J ω lam = 2νr 2 i r 2 o ω1−ω2 r 2 2 −r 2 1 where ν is the kinematic viscosity and . z,θ,t indicates averaging over the spatial directions z, θ and time t. For incompressible flows, it can be derived from the NS equations that J ω is conserved in the radial direction, ∂ r J ω = 0 (Eckhardt et al. 2007). In both cases, the values are made dimensionless by their respective laminar, conducting, values. For TC Nusselt becomes: The angular velocity flux J ω can easily be written in terms of the torque T on any of the cylinders: J ω = T (2πLρ) −1 with ρ being the fluid density. Consequently, the shear stress on the inner cylinder τ w,i is related to the angular velocity flux by τ w,i = ρJ ω r 2 i . Since part of our endeavour is to compare the effects of sand grain roughness on TC turbulence with the effect of sand grain roughness in other canonical systems (e.g. pipe flow), where the use of N u ω is not conventional, we choose to also present the global response in terms of the friction factor C f . Here we follow Lathrop et al. (1992), and and Re i is the inner cylinder Reynolds number Re i = riωid ν . The translation between N u ω and C f is straightforward, namely , which is different from eq. (2.10) by a factor which depends on the radius ratio η (Lathrop et al. 1992).

Numerical procedure
3.1. Roughness model Figure 2 exhibits the setup of the 'virtual' sand grain roughness model that is used in this research. The inner cylinder is divided up into square tiles of size 2k × 2k, each tile containing exactly 1 ellipsoid. The height L is slightly varied (0.85d ± 0.03d) to ensure that an integer amount of tiles fits into the domain. Unlike in the original model by Scotti (2006), we also introduce a random translation of the center of the ellipsoid by applying r i ∆θ and ∆z, where r i ∆θ, and ∆z are random uniform translations from the center of the 2k × 2k tile. This random translation allows for the surface to be more irregular and as such to relate more closely to a realistic sand grain surface. As also introduced by Scotti (2006), we employ a constant translation of the center of the ellipsoid in the radial direction, with ∆r = −0.5k from r = r i . It can be seen from figure 1(a) that part of the cylinder (≈ 15%) is not covered by rough elements.

Numerical method
The NS equations are spatially discretized by using a central second-order finitedifference scheme and solved in cylindrical coordinates by means of a semi-implicit procedure (Verzicco & Orlandi 1996;van der Poel et al. 2015a). The staggered grid is homogeneous in both the spanwise and streamwise directions (the axial and azimuthal directions respectively).
The wall-normal grid consists of a double cosine (Chebychev-type) grid stretching. Below the maximum roughness height, we employ a cosine stretching such that the  Table 2. (d) Schematic of an ellipsoidal building block of the rough surface. Every rectangular tile of size 2k × 2k contains exactly one ellipsoid. M indicates the center of the tile. The radii l 1 = 2.0k, l 2 = 1.4k, l 3 = 1.0k are kept constant for each ellipsoid to maintain a monodisperse rough surface. Randomness is ensured by giving the ellipsoid 5 degrees of freedom; 2 translational shifts of the center of the ellipsoid from the center of the tile M , (∆z and r i ∆θ) and 3 rotational degrees of freedom (α 1 , α 2 , α 3 ) from (r, θ, z) to (l 1 , l 2 , l 3 ). We also employ a constant translation of the center of the ellipsoid in the radial direction, with ∆r = −0.5k from r = r i . maximum grid spacing is always smaller than 0.5 times the viscous length scale. In the bulk of the fluid, we employ a second stretching, such that the maximum grid spacing in the bulk is approximately 1.5 times the viscous length scale. The minimum grid spacing is thus located at the position of the maximum roughness height, where we expect the highest shear stress, see table 1 for the exact values.
Time advancement is performed by using a fractional-step third-order Runge-Kutta scheme in combination with a Crank-Nicolson scheme for the implicit terms. The Courant-Friedrichs-Lewy (CFL) ( u∆t ∆x < 0.8) time-step constraint for the non-linear terms is enforced to ensure stability.
Sand grain roughness is implemented in the code by an immersed boundary method (IBM) (Fadlun et al. 2000). In the IBM, the boundary conditions are enforced by adding a body force f to the momentum equations (2.1) -(2.3). 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 preferentially to the normal surface vector to transfer the boundary conditions to the momentum equations. The IB has been validated previously (Fadlun et al. 2000;Iaccarino & Verzicco 2003;Stringano et al. 2006;Zhu et al. 2016Zhu et al. , 2017Zhu et al. , 2018b. Simulations are run until they become statistically stationary, such that Nusselt N u ω number remains constant to within ≈ 1%, which requirest ≈ 200 . Thereafter, we gather statistics until the results converge, which requirest ≈ 50. The resolution constraints of the domain are typically derived from the literature and are based on the grid spacing in '+' (viscous) units. Grid independence checks of the time averaged statistics are performed by ensuring that N u ω remains constant with increasing grid resolution in all directions and presented along with the results in Table 1.   Table 1: Input parameters, numerical resolution, and global response of the simulations. We run 4 sets of simulations, which are separated by an empty horizontal line. Within every set we keep T a constant. n θ × n z gives the number of ellipsoids in the streamwise (θ) and spanwise (z) directions, respectively. N ell expresses the resolution (N θ × N z ) per elementary building block of size 2k × 2k. N θ × N z × N r presents the total resolution of the computational domain. C f is the friction factor. N u ω is the dimensionless torque. k + s = 1.33k + is the equivalent sandgrain roughness height in viscous units, where k + is the size of the sandgrains (ellipsoids with axes 2k × 1.4k × k, see figure 2), also in viscous units. Re τ is the friction Reynolds number, Re τ = (d/2)u τ /ν. r + i ∆θ is the grid spacing in the streamwise and spanwise directions in viscous units, and ∆r + is the minimal grid spacing, at the maximum roughness height, in the wall-normal direction in viscous units. Figure 3 (left column) presents the streamwise (i.e. azimuthal) velocity profiles u + = (u(r)−u(r i )) z,θ,t /u τ (solid ) and angular velocity profiles ω + = u(r i )−r i u(r)/r z,θ,t /u τ (dashed ) versus the wall normal distance y + = r + − h + m , where . z,θ,t indicates averaging over the streamwise and spanwise directions and in time, and h + m is the mean roughness  s . Close overlap with the Nikuradse curve is observed in the transitionally rough regime. The overlap is slightly better for the angular velocity shift, for which we also obtain k + s = 1.33k + . The solid blue line represents the fully rough asymptote; ∆u + = 2.44 log(k + s ) + 5.2 − 8.5. The green lines represent the fully rough asymptotes obtained from the simulations, with κ u , κ ω , A u and A ω extracted from figure 3(g). The spread between statistically similar surfaces, with similar mean and maximum heights, is indicated by the vertical bar.

Roughness function
height. Every row corresponds to simulations at constant rotation rate of the inner cylinder (Taylor number), and increasing roughness height.
In line with the previous observations of Huisman et al. (2013) and Ostilla-Mónico et al. (2014), we also find that the logarithmic profiles of the streamwise velocity u + in smooth-wall Taylor-Couette do not fit the κ = 0.4 slope, as found in other wall bounded flows (e.g. pipe, boundary layer, channel) -for similar values of the friction Reynolds number Re τ . However, this asymptotic value is experimentally observed at very high shear rates of T a = O(10 12 ) (Huisman et al. 2013), much higher than can be obtained by the present DNS. The logarithmic profiles of angular velocity ω + have a slope that is closer to the κ = 0.4 asymptote (Ostilla-Mónico et al. 2015), especially for the higher Ta figure 3(g, h). Here we investigate the effects of roughness on both u + and ω + .
For rough wall simulations, the logarithmic region shifts downwards -a hallmark effect of a drag increasing surface. Figure 3 (right column) present the shifts, where ∆u + = u + s − u + r and ∆ω + = ω + s −ω + r . For lower T a, there is a small but observable difference between ∆u + and ∆ω + , see figure 3(b), whereas for the higher T a, this difference diminishes, see figure 3(f, h). Figure 4(a, b) presents the shift of the streamwise and angular velocity profiles, respectively, versus the equivalent roughness height k + s , for all T a. Care is taken to ensure overlap for varying T a and similar k + s , to study the T a dependence of ∆u + and ∆ω + . However, despite the varying T a numbers, all data collapse onto a single curve, with some scatter. Note that scatter is expected due to the randomness of the surfaces, which are reproduced for every simulation. To obtain an estimate of the expected scatter, we run two simulations with statistically similar surfaces. These are indicated by B3 and BY in table 1 and the velocity profiles are found in figures 3(c, d). We find a difference between the two cases of 0.2∆u + , 0.2∆ω + .
A comparison with the findings of Nikuradse (1933) can be carried out by scaling the Figure 5: (a) Azimuthal velocity u + (solid ) and the angular velocity ω + (dashed ) versus the wall normal distance y/k s , where y = (r − h m ) and h m is the mean roughness height. The three simulations with the highest roughness are plotted (D3, D4 and D5 respectively, see Table 1) to convey collapse of the profiles for the fully rough cases. (b) The Nikuradse constantB versus the equivalent sand grain roughness height k + s for both the azimuthal velocity (squares) and the angular velocity (diamonds). Horizontal black line atB = 6.0 gives the asymptotic value that is observed for fully rough behavior.
fully rough regime to obtain k + s = Ck + , where C is a constant that depends on the surface topology. In figure 5(a) we plot the velocity profiles versus (r − h m )/k s for the highest roughness (D3,D4 and D5 respectively, see Table 1). Excellent collapse of the D4 and D5 profiles indicates that those simulations are indeed fully rough. In this fully rough regime viscosity can be neglected (y k δ ν ), whereas the velocity profile is also independent of the system outer length scales (y d) i.e. the overlap argument (Pope 2000). The gradient of the velocity profile becomes; d U dy = uτ y Φ(y/k), where Φ(y/k) is a universal function that will go to 1/κ. Integration then gives; where B is a constant and y = r − h m .B is the Nikuradse constant. The roughness function in the fully rough regime, (i.e. the fully rough asymptote), is obtained by subtracting (4.1) from the smooth wall profile u + s = 1/κ log (y + ) + A and rescaling it to overlap with the Nikuradse data: In figure 4(a), the blue solid line is the fully rough asymptote, with κ, A,B as found in pipe flow (Pope 2000). The green solid line is the fully rough asymptote as obtained from our simulations. κ −1 u = 1.22 and A u = 8.0 are taken from the fit of the smooth wall simulation at identical T a as the fully rough cases (figure 3g). The fits are in the domain y + = [150, 500], as there the slope becomes approximately constant (figure 3h).B is plotted in figure 5(b), where we find thatB ≈ 6.0 for the fully rough cases. The mismatch of the slopes in the fully rough regime makes a rescaling to find k + s a priori impossible -a statement that we wish to emphasize. However, to proceed with the comparison of the transitionally rough cases in TC and pipe flow, we choose to rescale the fully rough cases (D4 and D5) with the Nikuradse fully rough asymptote in figure 4. We find that k + s = 1.33k + and very close collapse of our data with the Nikuradse data. Figure 6: Profiles of (a) the friction factor C f , and (b) the Nusselt number N u ω versus the roughness height. k/d is the roughness height k relative to the gap width d.
In parallel, we analyse the behavior of ∆ω + versus k + s , shown in figure 4(b). Again, the blue solid line represents the fully rough asymptote of Nikuradse. The green solid line is the fully rough asymptote obtained from fits (y + = [150, 500]) of the smooth wall angular velocity profile at identical T a as the fully rough cases, see figure 3(g). We find κ −1 ω = 2.17 (A ω = 3.7), close to the asymptotic value κ −1 = 2.44. Although the differences are marginally, ∆ω + fits to the Nikuradse data slightly better than ∆u + (note that also here the rescaling is with the Nikuradse data, k s = 1.33k + ). However, the major difference is the closeness of the fully rough asymptotes.
These results suggest that the near-wall effects of transitionally rough sand grains (and other rough surfaces) in TC flow are similar to the effects of transitionally rough sand grains (and other rough surfaces) in pipe flow (and other canonical systems). Further, we find that these transitionally rough effects are independent of slope of the logarithmic law, whereas in the fully rough regime, they, a priori, depend on this slope. This is confirmed with similar values of ∆u + at similar k + s , for varying T a, see figure 4. Also the similarity between ∆u + and ∆ω + in the transitionally rough regime confirms this, whereas the fully rough asymptotes are very dissimilar. We like to point out that simulations (or experiments) at high enough T a (= 10 12 (Huisman et al. 2013)), where κ = 0.4, then are expected to also give a collapse to the Nikuradse data in the fully rough regime. Figure 6(a) presents the friction factor C f versus the dimensionless roughness height k/d for varying Re i . For increasing T a numbers, the friction factor decreases (e.g. figure  7.23 in Pope (2000)), an indication that pressure is not yet dominant. For constant T a, as expected, the friction factor increases for increasing roughness height. Figure 6(b) then presents the global response in terms of N u ω . We observe an increase in N u ω for increasing T a, corresponding to the increased transport of the angular velocity that is due to the increased turbulent mixing. Higher roughness leads to increased J ω as compared to the smooth wall at the same T a, which also relates to a higher intensity of the turbulent mixing (the r 3 ( u r ω z,θ,t ) term of equation (2.8)) and more plumes ejecting from the boundary layer radially outwards (Zhu et al. 2017), on which we will elaborate in §4.3. By assuming a logarithmic profile, and integrating this profile over the entire gap (thereby neglecting the contributions of the viscous sublayer), we arrive at an implicit equation for the friction factor C f , namely the celebrated Prandtl's friction law:

Global response
(4.3) Figure 7 shows the friction factor C f versus the inner cylinder Reynolds number Re i , for both smooth wall as the rough wall data. An upward shift of the friction factor for increasing roughness height is consistent with what is observed for sand grains in pipe flow (Nikuradse 1933) and recently also for tranverse ribs in Taylor-Couette flow (Zhu et al. 2018b). Note that this upward shift, is directly related to the downward shift of the mean streamwise velocity profile (the roughness function). Since the friction factor and the Nusselt number are related, as expressed in equation (2.10), we expect the Nusselt number to increase, for increasing roughness height. This is confirmed in figure 7(b), where we plot the Nu number versus the Ta number. The number of simulations with constant k/d is limited, and we vary the Ta number over 2 decades only. However, we clearly observe that the asymptotic, ultimate scaling of N u ω ∝ T a 0.5 , as found for fully rough transverse ribs in (Zhu et al. 2018b), is not reached. This is expected, as for Ta numbers the system is still in the transitionally rough TC regime, where vicous forces still strongly contribute to N u ω , and only the inner cylinder is covered with roughness.

Flow structures
To obtain a qualitative understanding of the effect of inner cylinder roughness on the turbulent flow in the gap, we present two series of snapshots of the streamwise azimuthal velocity u θ (r, θ, z). Figure 8(a − c) exhibit the snapshots for T a = 5.0 × 10 7 . It is known, and observed here, that for this Taylor number the coherence length of the dominant flow structures becomes smaller than the gap width d, and turbulence develops in the bulk ). On the other hand, the boundary layers remain predominantly laminar and as such the regime is referred to as the 'classical regime' of TC turbulence. A divergent colormap is chosen to highlight the turbulent structures in the bulk. A snapshot for the smooth inner cylinder simulation is presented in figure 8(a). At z/d ≈ 0.3, one observes an ejecting structure (plume) that detaches from the inner cylinder laminar boundary layer at the location of an adverse pressure gradient. Locally, where this ejecting plume detaches, the flow will be different (i.e. more chaotic) to that in the other parts Ultimate turbulent state Contour fields of the instanteneous azimutal velocity profile u θ (r, θ, z, t) for T a = 1.0 × 10 9 in the meridional plane. (d) Smooth inner cylinder (DS) with many plumes ejecting, considerably more chaotic than in (a). (e) Inner cylinder wall roughness (indicated in grey) k/d = 0.035 (D2) and (f ) inner wall roughness k/d = 0.055 (D5). For the rough cases, we observe more plumes ejecting and more mixing in the bulk, leading to enhanced radial transport of the angular velocity J ω , expressed in a higher N uω. Figure 8 of the boundary layer. As such, one also expects the local variables (e.g. skin friction, turbulence intensity) to be different. Later we will therefore also employ local averaging, to investigate the spatial differences in the flow associated with this structure (Zhu et al. (2018a); van der Poel et al. (2015b); Ostilla-Mónico et al. (2016)). The ejecting and impacting (located at z/d ≈ 1.3) plumes have very strong radial velocity components u r . From the first term on the right-hand side of (2.8), r 3 ( u r ω z,θ,t ), we then directly see that they strongly contribute to N u ω . This brings us to the remaining (figure 8 b and c) snapshots. Many more, small, plumes are seen to eject from the inner cylinder. The roughness there promotes the detachment of ejecting structures and in that way contributes to a higher N u ω . An increase in the level of turbulence, as suggested by the increased level of turbulence dissipation, is quantitatively reflected by a decrease in the Kolmogrov scale (η = (ν 3 / ) 1/4 ), namely η/d = 7.1 × 10 −3 for the smooth wall case BS and η/d = 6.5 × 10 −3 for the highest roughness case B5. Note that the decrease in the Kolmogorov scale η is only small, since η/d ∝ ( d 4 /ν 3 ) −1/4 . For TC flow, the volume averaged dissipation rate is related to the angular velocity transport N u ω with: = ν 3 d −4 σ −2 (N u ω − 1)T a + lam , where lam is the laminar volume averaged dissipation rate, d is the gap width of the setup and σ is a geometric parameter (Eckhardt et al. 2007). As such, we see that η/d ∝ N u −1/4 ω only. Figure 8(d − f ) presents snapshots of a flow in the ultimate turbulent state at T a = 1.0 × 10 9 . Although less pronounced than for T a = 5.0 × 10 7 , we still observe distinct ejecting and impacting regions, indicating the survival of the turbulent Taylor rolls. A similar rationale as applied above, to the classical turbulence case, can also be used to explain the enhancement of the Nusselt number for rough inner cylinders in the ultimate turbulent state. In fact, we can also observe more intense plumes for the highest roughness (D5, figure 8f ), in comparison to a lower roughness case (D2, figure 8e). Note that here we do not observe the stable vortex formation in between roughness elements and the associated ejection of plumes from sharp peaks, as was reported by Zhu et al. (2016) for grooved cylinders, for similar Taylor numbers and roughness heights. The increase in the turbulence level is also quantitatively confirmed by a decrease in the Kolmogorov scale here, η/d = 2.7 × 10 −3 for the smooth wall case DS and η/d = 2.1 × 10 −3 for the highest roughness case D5.

Roughness sublayer
The existence of Taylor roll structures is already anticipated in the snapshots of the instantaneous flow in figure 8, from which we observe the ejecting and impacting plume regions. Contour plots of the time and azimuthally averaged azimuthal velocity field, as presented in figure 9, confirm this. Note that the Taylor roll is spatially fixed, allowing for convenient averaging over impacting (solid line), shearing (dashed line) and ejecting (dashed dotted line) regions, a method that we also have employed in RB flow (van der Poel et al. 2015b). For an increasing roughness height, the white region (representing u θ θ ≈ 0.5) shifts radially outwards and the azimuthal velocity in the bulk increases. This process previously has been seen in Zhu et al. (2018b), where it is referred to as the bulk velocity 'getting slaved to' to the velocity of a cylinder covered with roughness, reflecting the enhanced drag on the rough side.
A better impression of the local fluid flow disturbances induced by the roughness elements, is obtained from the time-averaged azimuthal velocity u θ , a contour of which is shown in figure 9. We zoom in on only a few roughness elements and overlay the contour with isolines of u + θ . It can quickly be seen that local disturbances are limited to a region of only a few times the roughness height, above which the isolines relax to approximate u θ (r, θ, z, t) = u θ θ (r, z) +ũ θ (r, θ, z) + u θ (r, θ, z, t) (4.4) where u θ (r, θ, z, t) = u θ (r, θ, z, t) − u θ (r, θ, z) is the temporal fluctuation andũ θ (r, θ, z) = u θ (r, θ, z) − u θ θ (r, z) is the component that is strongly related to the roughness induced disturbances and therefore termed the form-induced (or dispersive) velocity fluctuation (for brevity now written asũ). Note that by applying the triple decomposition to the full NS equations, one will recover the related form-induced stress tensor ũ iũj θ . However, here we only discussũ. The root mean square of the form-induced fluctuations √ũ 2 + at various heights above the roughness is given in figure 11(a). The horizontal axis corresponds to the roughness elements shown in figure 9. Already for r/k = 4.5 (with k being the roughness height) we find it hard to detect spatial fluctuations along θ. For r/k = 12.0, the line is barely distinguishable from the horizontal axis. If we average the lines over the azimuthal direction, we obtain the behavior of the dispersive fluctations as a function of the wall normal distance; see figure 11(b). We plot the lines for three respective axial locations, that is the impacting, sheared and ejecting regions, and find little variance in the wall normal extent of the form-induced fluctuations, for the varying locations. The vertical black line represent the maximum roughness height, below wich we average only over fluid regions. It is convenient to introduce a height at which dispersive fluctuations vanish, and this height is commonly known as the 'roughness sublayer height' h r . In this research, we define y + = h + r where ũ +2 θ,z = 0.01 ū + θ and find h r = 3.70k (h r = 2.78k s ). This value agrees well with a roughness sublayer height of 2 k s 5, as typically found in other canonical flows (Pokrajac et al. 2007). Figure 10: Contour of the time-averaged azimuthal velocity u + θ , zoomed in on only a few roughness elements, indicated in grey, for T a = 1.0 × 10 9 (D2). Isolines of u + θ overlay the contour. On the vertical axis, we display the wall normal coordinate normalized by the gap width d (left) and in viscous units (right). On the horizontal axis, we display the azimuthal coordinate, normalized by the gap width d (below ) and in viscous units (above). The arrow indicates the direction of the mean flow.

Radial profiles
We plot the time and azimuthally averaged azimuthal velocity for T a = 5.0 × 10 8 in figure 12(a). The roughness covers the inner cylinder, i.e. (r − r i )/d 0.07. Below the roughness crest, the velocity is averaged over the entire domain, including the solid rough regions. For increasing roughness height, see the legend for viscous roughness heights, the azimuthal velocity in the bulk increases. Figure 12(b) then presents the corresponding azimuthal velocity profiles for constant k/d = 0.060 ± 0.002 and varying Taylor number. Figure 13(a) shows the double-averaged radial profiles (u θ ) + rms = u 2 θ θ,z − u θ 2 θ,z 1/2 θ,z /u τ for constant Taylor number T a = 1.0 × 10 9 . For the smooth wall case, the peak of the rms of the azimuthal velocity is located at u θ,rms | max ≈ 10y + . This agrees with the smooth case in (?) for TC flow, but is closer to the solid boundary than that is found in channel flow; u θ,rms | max ≈ 12y + . We observe a slight decrease in the viscous scaled turbulence intensity (u θ ) + rms for increasing roughness heights, close to the inner cylinder. Further outside the profiles, the profiles collapse. This is once more a sign that the effects of roughness are only confined to the fuid flow close to the wall.
It has already been mentioned that the angular velocity flux J ω (y + ) is conserved in the radial direction, see section 4.2. However, the individual components -i.e. the viscous J ω ν (y + ) and Reynolds stress J ω urω (y + ) terms -are functions of the wall normal  (radial) coordinate. The profiles of these terms are shown in figure 13(b). We normalize all terms with max(J ω (y + )). The viscous stress terms are presented as solid lines and the Reynolds stress terms as dashed lines. The blue lines represent the smooth wall case. Very close to the inner cyliner (y + = 1), J ω ≈ J ω ν ; this is expected, since there the gradient of the mean streamwise velocity is maximum and the wall normal component Figure 13: (a) Root mean square of the mean streamwise velocity (u θ ) + rms = u 2 θ θ,z − u θ 2 θ,z 1/2 θ,z /u τ versus the radius (r − r i )/d for T a = 1.0 × 10 9 . We observe an overall decrease of the viscous scaled turbulence intensity for increasing roughness height k + s . (b) The viscous J ω ν (y + ) ≡ −νr 3 ∂ r ω z,θ,t Solid line and Reynolds stress J ω urω (y + ) ≡ − r 3 ν ∂ r ω z,θ,t (dashed line) terms of the conserved angular velocity flux J ω (r + ) versus the wall normal distance y + for T a = 1.0×10 9 . The sum of the individual terms represents the total conserved angular velocity flux radially outwards. The black horizontal line at J ω = 1.0 is the sum of the two blue lines for the smooth wall case. It can be seen that J ω is conserved above the maximum roughness height (indicated by cross markers).
of the velocity vector goes to zero. On the far right of figure 13(b), the gradient of the mean velocity approaches zero in the bulk, and the correlation function between the wall normal and angular velocity goes to a maximum; as such J ω ≈ J ω urω . The black line represents the sum of the viscous and Reynold stress terms, for the smooth (blue) case, and is independent of y + . The situation for the rough cylinder cases is more complex. For increasing roughness height, the viscous stress reaches a maximum below the maximum roughness height. This can be explained by the recirculation zones behind the roughness elements. Then, above the roughness, the viscous stress goes to zero in the bulk. For the Reynolds stress terms the increase is monotonic, and similar to the smooth wall case. However, we observe a steeper increase for the rough cases. Note that under the maximum roughness height, the viscous and Reynolds stress terms do not add up to unity, since the rough surface acts as a radial dependent momentum source term to the NS equations there.

Summary & conclusions
We have performed Direct Numerical Simulations of turbulent Taylor-Couette flow with inner cylinder rotation and inner cylinder sand grain roughness. The Taylor number ranges from T a = 1.0 × 10 7 (Re τ = 82) up to T a = 1.0 × 10 9 (Re τ = 635), covering thereby both the classical and the ultimate regimes of turbulent TC flow. In particular, we studied the effects of the roughness height on the fluid flow in the transitionally rough and fully rough regimes, with the equivalent sand grain roughness height ranging from k + s = 5 to k + s = 92. We modelled the sand grains as randomly rotated and translated ellipsoids of constant size and shape (monodisperse), similar to the model proposed by Scotti (2006). The surface was implemented into a second order finite difference code by means of the Immersed Boundary Method.
We confirm an increase in the dimensionless torque, expressed as the Nusselt number, for increasing roughness height. This is attributed to the enhanced boundary layer detachment, resulting in plume ejection regions, which we observed in snapshots of the azimuthal velocity field. The plumes contain a strong radial velocity component, and as such contribute strongly to the Reynolds stress term of the angular velocity flux. This mechanism is analogous to what was found for the Nu increase for grooved TC turbulence by Zhu et al. (2016).
To quantify the degree of roughness induced disturbances to the velocity field, we measured the dispersive fluctuations of the azimuthal velocity component,ũ θ = u θ t − u θ t,θ . This dispersive term was obtained by means of a triple decomposition to the azimuthal velocity. We defined the height of the roughness sublayer h r there where the dispersive fluctuations become very small, such that ũ +2 θ,z = 0.01 ū + θ and found the height of the roughness sublayer to be h r = 2.78k s . This height of the roughness sublayer compares well with values found for other canonical systems (Pokrajac et al. 2007).
The hallmark of turbulent flows over rough walls is the shift of the logarithmic streamwise velocity profile ∆u + θ . The shift is a function of any parameters describing the roughness topology. Here we focused in particular on the effect of the sand grain size k and the roughness function becomes ∆u + θ = f (k). It was shown in Huisman et al. (2013) that the constants of the logarithmic law are not constant in the Taylor number range of our simulations. Hence we proposed the generalization; u + θ = f 1 (Re i ) log(y + ) + f 2 (Re i ) − ∆u + θ . We performed simulations at four Ta and various roughness heights and ensured that the k + s range for the various Ta numbers overlaps. First, we concluded that the velocity shift is independent of Ta, despite the Ta dependence of the constant in the logarithmic layer. As such, all simulations collapse onto a single curve. Second, we saw a strong overlap between the roughness function calculated from our DNS in TC flow and the seminal work by Nikuradse (1933) on monodisperse sand grains in pipe flow, in the transitionally rough regime. We found k + s = 1.33k. Only for very low k + s values, close to the hydrodynamically smooth regime, we found that the simulations slightly differ from the Nikuradse curve.
It is remarkable that the Hama roughness function appears to be universal for similar surfaces in such different systems. Note in particular that we have a streamwise curvature and strong secondary motions (Taylor rolls), which were absent in the pipe flow experiments of Nikuradse. As such, our findings point towards a universal behavior of the roughness function for very different fluid flow systems. However, many more comparison studies, of identical rough surfaces and varying fluid flow systems, are needed to confirm this notion.  measured with respect to the inner cylinder location r = r 0 .k/k max is the mean roughness height respective to the maximum roughness height.k is the first moment of the roughness height distribution; 1 S S k(z, θ)dS. k std is the standard deviation of the surface height distribution ( 1 S S (k(z, θ) −k) 2 dS) 1/2 , Sk is the skewness 1 k 3 std S S (k(z, θ) −k) 3 dS and Ku the kurtosis, 1 k 4 std S S (k(z, θ) −k) 4 dS. k rms is the root mean square of the mean height ( 1 S S (k(z, θ)) 2 dS) 1/2 . k p is the mean peak height. The peaks are obtained by a peak finding algorithm ('find peak local ' in van der Walt et al. (2014)) with zero threshold and a minimum spacing of the peaks of 4 grid nodes. ES is the effective slope, ES = 1 L θ Lz L θ Lz | ∂k(z,θ) ∂(riθ) |dzd(r i θ) as introduced by Napoli et al. (2008). It can be shown that the ES parameter is twice the often used solidity parameter λ (Jimenez 2004). α is the surface gradient in the streamwise direction, ∂k(z,θ) ∂θ , and α rms is the root mean square of this distribution. S f is the total frontal projected area of the roughness (in the streamwise direction). S s is the total windward wetted area of the roughness. Λ is the density parameter which relates the latter two parameters by Λ = ( S S f )( S f Ss ) −1.6 where S is the total area of the surface without roughness (Sigal & Danberg 1990). The fractal dimension of surface B2 is calculated from the slope of the power density spectrum versus the wavenumber, C ∝ q −4.65 , for details we refer to (Persson et al. 2005). Here we find a Hurst exponent H of 1.32 and a fractal dimension D of 1.68. Such we conclude that the surface can at no scale be considered fractal.