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

TL;DR

This study uses direct numerical simulations to investigate how sand grain roughness affects turbulence in Taylor--Couette flow, revealing that roughness causes a downward shift in the logarithmic layer similar to pipe flow, across a range of flow conditions.

Contribution

The paper introduces an improved simulation model for monodisperse sand grain roughness in Taylor--Couette flow and analyzes its impact on turbulence and skin-friction drag.

Findings

Roughness causes a downward shift in the logarithmic layer.

Shift behavior is similar to that observed in pipe flow.

Effect consistent across different Taylor numbers.

Abstract

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.

Figures

Tables

Table 1. Table 1 : Input parameters, numerical resolution, and global response of the simulations. We run 4 4 4 sets of simulations, which are separated by an empty horizontal line. Within every set we keep T ​ a 𝑇 𝑎 Ta constant. n θ × n z subscript 𝑛 𝜃 subscript 𝑛 𝑧 n_{\theta}\times n_{z} gives the number of ellipsoids in the streamwise ( θ 𝜃 \theta ) and spanwise ( z 𝑧 z ) directions, respectively. N e ​ l ​ l subscript 𝑁 𝑒 𝑙 𝑙 N_{ell} expresses the resolution ( N θ × N z subscript 𝑁 𝜃 subscript 𝑁 𝑧 N_{\theta}\times N_{z} ) per elementary building block of size 2 ​ k × 2 ​ k 2 𝑘 2 𝑘 2k\times 2k . N θ × N z × N r subscript 𝑁 𝜃 subscript 𝑁 𝑧 subscript 𝑁 𝑟 N_{\theta}\times N_{z}\times N_{r} presents the total resolution of the computational domain. C f subscript 𝐶 𝑓 C_{f} is the friction factor. N ​ u ω 𝑁 subscript 𝑢 𝜔 Nu_{\omega} is the dimensionless torque. k s + = 1.33 ​ k + superscript subscript 𝑘 𝑠 1.33 superscript 𝑘 k_{s}^{+}=1.33k^{+} is the equivalent sandgrain roughness height in viscous units, where k + superscript 𝑘 k^{+} is the size of the sandgrains (ellipsoids with axes 2 ​ k × 1.4 ​ k × k 2 𝑘 1.4 𝑘 𝑘 2k\times 1.4k\times k , see figure 2 ), also in viscous units. R ​ e τ 𝑅 subscript 𝑒 𝜏 Re_{\tau} is the friction Reynolds number, R ​ e τ = ( d / 2 ) ​ u τ / ν 𝑅 subscript 𝑒 𝜏 𝑑 2 subscript 𝑢 𝜏 𝜈 Re_{\tau}=(d/2)u_{\tau}/\nu . r i + ​ Δ ​ θ superscript subscript 𝑟 𝑖 Δ 𝜃 r_{i}^{+}\Delta\theta is the grid spacing in the streamwise and spanwise directions in viscous units, and Δ ​ r + Δ superscript 𝑟 \Delta r^{+} is the minimal grid spacing, at the maximum roughness height, in the wall-normal direction in viscous units.

Case	$Ta$	$k/d$	$n_{\theta }\times n_z$	$N_{ell}$	$N_{\theta }\times N_z\times N_r$	$C_f$	$N{u}_{\omega }$	$k_s^{+}$	$R{e}_{\tau }$	$r_i^{+}\Delta \theta$	$\Delta r^{+}$
AS	$1.0\times {10}^7$	–	–	–	$280\times 240\times 256$	$0.161$	$6.41$	$0.00$	$81.9$	$1.51$	$0.26$
A1	$1.0\times {10}^7$	$0.022$	$59\times 48$	$64$	$472\times 384\times 256$	$0.170$	$6.79$	$4.96$	$84.3$	$1.10$	$0.22$
A2	$1.0\times {10}^7$	$0.038$	$34\times 28$	$64$	$272\times 224\times 300$	$0.182$	$7.31$	$8.93$	$87.4$	$1.60$	$0.20$
A3	$1.0\times {10}^7$	$0.055$	$24\times 20$	$144$	$288\times 240\times 300$	$0.191$	$7.63$	$12.93$	$89.3$	$1.55$	$0.22$
A4	$1.0\times {10}^7$	$0.073$	$18\times 15$	$256$	$288\times 240\times 400$	$0.208$	$8.30$	$18.00$	$93.2$	$1.62$	$0.20$
BS	$5.0\times {10}^7$	–	–	–	$280\times 240\times 448$	$0.098$	$8.78$	$0.00$	$143.3$	$2.67$	$0.18$
B1	$5.0\times {10}^7$	$0.022$	$60\times 48$	$64$	$720\times 576\times 600$	$0.107$	$9.58$	$8.86$	$149.7$	$1.11$	$0.17$
B2	$5.0\times {10}^7$	$0.038$	$34\times 28$	$64$	$272\times 224\times 600$	$0.124$	$11.06$	$16.43$	$160.9$	$3.08$	$0.19$
B3	$5.0\times {10}^7$	$0.055$	$24\times 20$	$144$	$288\times 240\times 600$	$0.136$	$12.12$	$24.49$	$168.4$	$3.06$	$0.20$
BY	$5.0\times {10}^7$	$0.055$	$24\times 20$	$144$	$288\times 240\times 600$	$0.136$	$12.20$	$24.57$	$169.0$	$3.06$	$0.20$
B4	$5.0\times {10}^7$	$0.073$	$18\times 15$	$256$	$288\times 240\times 600$	$0.148$	$13.24$	$33.99$	$176.0$	$3.20$	$0.22$
B5	$5.0\times {10}^7$	$0.087$	$14\times 12$	$400$	$280\times 240\times 600$	$0.155$	$13.84$	$41.71$	$180.0$	$3.37$	$0.23$
CS	$5.0\times {10}^8$	–	–	–	$512\times 512\times 640$	$0.060$	$16.94$	$0.00$	$354.0$	$3.62$	$0.23$
C1	$5.0\times {10}^8$	$0.019$	$68\times 56$	$144$	$816\times 672\times 800$	$0.076$	$21.48$	$20.39$	$398.6$	$2.56$	$0.29$
C2	$5.0\times {10}^8$	$0.026$	$50\times 40$	$256$	$800\times 640\times 800$	$0.084$	$23.66$	$29.11$	$418.4$	$2.74$	$0.26$
C3	$5.0\times {10}^8$	$0.034$	$38\times 32$	$256$	$608\times 512\times 800$	$0.091$	$25.77$	$39.95$	$436.6$	$3.76$	$0.23$
C4	$5.0\times {10}^8$	$0.041$	$32\times 26$	$256$	$512\times 416\times 800$	$0.096$	$26.98$	$48.55$	$446.7$	$4.57$	$0.28$
CX	$5.0\times {10}^8$	$0.041$	$32\times 26$	$256$	$512\times 416\times 1000$	$0.096$	$27.01$	$48.66$	$447.0$	$4.57$	$0.20$
C5	$5.0\times {10}^8$	$0.047$	$28\times 23$	$324$	$504\times 414\times 800$	$0.101$	$28.50$	$56.98$	$459.2$	$4.77$	$0.35$
DS	$1.0\times {10}^9$	–	–	–	$512\times 512\times 640$	$0.054$	$21.70$	$0.00$	$476.5$	$4.87$	$0.30$
D1	$1.0\times {10}^9$	$0.026$	$50\times 40$	$256$	$800\times 640\times 1000$	$0.080$	$31.81$	$40.13$	$576.8$	$3.78$	$0.27$
D2	$1.0\times {10}^9$	$0.034$	$38\times 32$	$400$	$760\times 640\times 1000$	$0.086$	$34.20$	$54.74$	$598.2$	$4.12$	$0.29$
D3	$1.0\times {10}^9$	$0.041$	$32\times 26$	$484$	$704\times 572\times 1200$	$0.089$	$35.75$	$66.47$	$611.6$	$4.55$	$0.23$
D4	$1.0\times {10}^9$	$0.047$	$28\times 23$	$784$	$784\times 644\times 1000$	$0.094$	$37.45$	$77.66$	$625.9$	$4.18$	$0.45$
D5	$1.0\times {10}^9$	$0.055$	$24\times 20$	$1024$	$768\times 640\times 1000$	$0.097$	$38.54$	$91.95$	$635.0$	$4.33$	$0.46$

Table 2. Table 2 : Parameters describing the surface geometry. From left to right; Case number in accordance with table 1 . The roughness height is measured with respect to the inner cylinder location r = r 0 𝑟 subscript 𝑟 0 r=r_{0} . k ¯ / k m ​ a ​ x ¯ 𝑘 subscript 𝑘 𝑚 𝑎 𝑥 \bar{k}/k_{max} is the mean roughness height respective to the maximum roughness height. k ¯ ¯ 𝑘 \bar{k} is the first moment of the roughness height distribution; 1 S ​ ∫ S k ​ ( z , θ ) ​ 𝑑 S 1 𝑆 subscript 𝑆 𝑘 𝑧 𝜃 differential-d 𝑆 \frac{1}{S}\int_{S}k(z,\theta)dS . k s ​ t ​ d subscript 𝑘 𝑠 𝑡 𝑑 k_{std} is the standard deviation of the surface height distribution ( 1 S ​ ∫ S ( k ​ ( z , θ ) − k ¯ ) 2 ​ 𝑑 S ) 1 / 2 superscript 1 𝑆 subscript 𝑆 superscript 𝑘 𝑧 𝜃 ¯ 𝑘 2 differential-d 𝑆 1 2 (\frac{1}{S}\int_{S}(k(z,\theta)-\bar{k})^{2}dS)^{1/2} , Sk is the skewness 1 k s ​ t ​ d 3 ​ S ​ ∫ S ( k ​ ( z , θ ) − k ¯ ) 3 ​ 𝑑 S 1 superscript subscript 𝑘 𝑠 𝑡 𝑑 3 𝑆 subscript 𝑆 superscript 𝑘 𝑧 𝜃 ¯ 𝑘 3 differential-d 𝑆 \frac{1}{k_{std}^{3}S}\int_{S}(k(z,\theta)-\bar{k})^{3}dS and Ku the kurtosis, 1 k s ​ t ​ d 4 ​ S ​ ∫ S ( k ​ ( z , θ ) − k ¯ ) 4 ​ 𝑑 S 1 superscript subscript 𝑘 𝑠 𝑡 𝑑 4 𝑆 subscript 𝑆 superscript 𝑘 𝑧 𝜃 ¯ 𝑘 4 differential-d 𝑆 \frac{1}{k_{std}^{4}S}\int_{S}(k(z,\theta)-\bar{k})^{4}dS . k r ​ m ​ s subscript 𝑘 𝑟 𝑚 𝑠 k_{rms} is the root mean square of the mean height ( 1 S ​ ∫ S ( k ​ ( z , θ ) ) 2 ​ 𝑑 S ) 1 / 2 superscript 1 𝑆 subscript 𝑆 superscript 𝑘 𝑧 𝜃 2 differential-d 𝑆 1 2 (\frac{1}{S}\int_{S}(k(z,\theta))^{2}dS)^{1/2} . k p subscript 𝑘 𝑝 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, E ​ S = 1 L θ ​ L z ​ ∫ L θ ∫ L z | ∂ k ​ ( z , θ ) ∂ ( r i ​ θ ) | ​ 𝑑 z ​ d ​ ( r i ​ θ ) 𝐸 𝑆 1 subscript 𝐿 𝜃 subscript 𝐿 𝑧 subscript subscript 𝐿 𝜃 subscript subscript 𝐿 𝑧 𝑘 𝑧 𝜃 subscript 𝑟 𝑖 𝜃 differential-d 𝑧 𝑑 subscript 𝑟 𝑖 𝜃 ES=\frac{1}{L_{\theta}L_{z}}\int_{L_{\theta}}\int_{L_{z}}|\frac{\partial k(z,\theta)}{\partial(r_{i}\theta)}|dzd(r_{i}\theta) as introduced by Napoli et al. ( 2008 ) . It can be shown that the ES parameter is twice the often used solidity parameter λ 𝜆 \lambda (Jimenez, 2004 ) . α 𝛼 \alpha is the surface gradient in the streamwise direction, ∂ k ​ ( z , θ ) ∂ θ 𝑘 𝑧 𝜃 𝜃 \frac{\partial k(z,\theta)}{\partial\theta} , and α r ​ m ​ s subscript 𝛼 𝑟 𝑚 𝑠 \alpha_{rms} is the root mean square of this distribution. S f subscript 𝑆 𝑓 S_{f} is the total frontal projected area of the roughness (in the streamwise direction). S s subscript 𝑆 𝑠 S_{s} is the total windward wetted area of the roughness. Λ Λ \Lambda is the density parameter which relates the latter two parameters by Λ = ( S S f ) ​ ( S f S s ) − 1.6 Λ 𝑆 subscript 𝑆 𝑓 superscript subscript 𝑆 𝑓 subscript 𝑆 𝑠 1.6 \Lambda=(\frac{S}{S_{f}})(\frac{S_{f}}{S_{s}})^{-1.6} where S 𝑆 S is the total area of the surface without roughness (Sigal & Danberg, 1990 ) . The fractal dimension of surface B ​ 2 𝐵 2 B2 is calculated from the slope of the power density spectrum versus the wavenumber, C ∝ q − 4.65 proportional-to 𝐶 superscript 𝑞 4.65 C\propto q^{-4.65} , for details we refer to (Persson et al. , 2005 ) . Here we find a Hurst exponent H 𝐻 H of 1.32 1.32 1.32 and a fractal dimension D 𝐷 D of 1.68 1.68 1.68 . Such we conclude that the surface can at no scale be considered fractal.

Case	$\overline{k}/k_{max}$	$\overline{k}/k$	$k_{std}/k$	Sk	Ku	$k_{rms}/k$	$k_p/k$	ES	${\alpha }_{rms}$	$S_f/r_i^2$	$S_s/r_i^2$	$\Lambda$
A1	$0.38$	$0.57$	$0.35$	$-0.06$	$-0.77$	$0.67$	$1.11$	$0.44$	$0.17$	$0.15$	$0.21$	$9$
A2	$0.38$	$0.57$	$0.35$	$-0.08$	$-0.78$	$0.67$	$1.11$	$0.44$	$0.17$	$0.16$	$0.21$	$9$
A3	$0.37$	$0.56$	$0.34$	$-0.09$	$-0.76$	$0.66$	$1.42$	$0.46$	$0.11$	$0.17$	$0.24$	$9$
A4	$0.37$	$0.56$	$0.36$	$-0.11$	$-0.87$	$0.67$	$1.05$	$0.48$	$0.12$	$0.18$	$0.28$	$10$
B1	$0.38$	$0.57$	$0.35$	$-0.06$	$-0.77$	$0.67$	$1.11$	$0.44$	$0.17$	$0.15$	$0.21$	$9$
B2	$0.38$	$0.57$	$0.35$	$-0.08$	$-0.78$	$0.67$	$1.12$	$0.44$	$0.17$	$0.16$	$0.21$	$9$
B3	$0.37$	$0.56$	$0.34$	$-0.09$	$-0.76$	$0.66$	$1.04$	$0.46$	$0.11$	$0.17$	$0.24$	$9$
BY	$0.38$	$0.57$	$0.35$	$-0.07$	$-0.75$	$0.67$	$1.07$	$0.47$	$0.15$	$0.17$	$0.25$	$9$
B4	$0.37$	$0.56$	$0.36$	$-0.11$	$-0.87$	$0.67$	$1.05$	$0.48$	$0.12$	$0.18$	$0.28$	$10$
B5	$0.37$	$0.56$	$0.35$	$-0.13$	$-0.86$	$0.66$	$1.02$	$0.44$	$0.11$	$0.17$	$0.26$	$9$
C1	$0.38$	$0.57$	$0.35$	$-0.08$	$-0.81$	$0.67$	$1.07$	$0.47$	$0.13$	$0.17$	$0.24$	$9$
C2	$0.37$	$0.56$	$0.35$	$-0.09$	$-0.85$	$0.66$	$1.03$	$0.47$	$0.11$	$0.16$	$0.24$	$10$
C3	$0.37$	$0.56$	$0.36$	$-0.11$	$-0.85$	$0.67$	$1.04$	$0.48$	$0.07$	$0.18$	$0.29$	$10$
C4	$0.37$	$0.56$	$0.35$	$-0.10$	$-0.80$	$0.66$	$1.04$	$0.47$	$0.12$	$0.17$	$0.25$	$10$
C5	$0.37$	$0.56$	$0.35$	$-0.10$	$-0.82$	$0.67$	$1.04$	$0.48$	$0.12$	$0.17$	$0.26$	$10$
D1	$0.37$	$0.56$	$0.35$	$-0.09$	$-0.85$	$0.66$	$1.03$	$0.47$	$0.11$	$0.16$	$0.24$	$10$
D2	$0.37$	$0.56$	$0.36$	$-0.11$	$-0.85$	$0.67$	$1.03$	$0.48$	$0.06$	$0.19$	$0.30$	$10$
D3	$0.37$	$0.56$	$0.35$	$-0.11$	$-0.82$	$0.66$	$1.01$	$0.48$	$0.08$	$0.17$	$0.27$	$11$
D4	$0.37$	$0.56$	$0.36$	$-0.10$	$-0.83$	$0.67$	$1.02$	$0.49$	$0.08$	$0.18$	$0.29$	$11$
D5	$0.37$	$0.56$	$0.35$	$-0.12$	$-0.85$	$0.67$	$1.01$	$0.50$	$0.08$	$0.19$	$0.32$	$10$

Equations

$$u^{+}(y^{+})=\frac{1}{\kappa }\log y^{+}+A-\Delta u^{+}$$

$$u^{+}(y^{+})=f_1(R{e}_i)\log (y^{+})+f_2(R{e}_i)-\Delta u^{+}$$

$${\partial }_{\hat{t}}{\hat{u}}_r+\hat{\mathbf{u}}\cdot \hat{\nabla }{\hat{u}}_r-\frac{{\hat{u}}_{\theta }^2}{\hat{r}}=-{\partial }_{\hat{r}}\hat{P}+\frac{f(\eta )}{T{a}^{1/2}}({\hat{\nabla }}^2{\hat{u}}_r-\frac{{\hat{u}}_r}{{\hat{r}}^2}-\frac{2}{{\hat{r}}^2}{\partial }_{\theta }{\hat{u}}_{\theta })-R{o}^{-1}{\hat{u}}_{\theta }$$

$${\partial }_{\hat{t}}{\hat{u}}_{\theta }+\hat{\mathbf{u}}\cdot \hat{\nabla }{\hat{u}}_{\theta }-\frac{{\hat{u}}_{\theta }{\hat{u}}_r}{\hat{r}}=-\frac{1}{\hat{r}}{\partial }_{\hat{\theta }}\hat{P}+\frac{f(\eta )}{T{a}^{1/2}}({\hat{\nabla }}^2{\hat{u}}_{\theta }-\frac{{\hat{u}}_{\theta }}{{\hat{r}}^2}+\frac{2}{{\hat{r}}^2}{\partial }_{\theta }{\hat{u}}_r)+R{o}^{-1}{\hat{u}}_r$$

$${\partial }_{\hat{t}}{\hat{u}}_z+\hat{\mathbf{u}}\cdot \hat{\nabla }{\hat{u}}_z=-{\partial }_{\hat{z}}\hat{P}+\frac{f(\eta )}{T{a}^{1/2}}({\hat{\nabla }}^2{\hat{u}}_z)$$

$$\hat{\nabla }\cdot \hat{\mathbf{u}}=0$$

$$f(\eta )=\frac{(1+\eta {)}^3}{8{\eta }^2}$$

$$Ta=\frac{(1+\eta {)}^4}{64{\eta }^2}\frac{(r_o-r_i{)}^2(r_i+r_o{)}^2({\omega }_i-{\omega }_o{)}^2}{{\nu }^2}.$$

$$R{o}^{-1}=\frac{2{\omega }_{i}d}{r_i|{\omega }_i-{\omega }_o|}.$$

$$J^{\omega }=r^3(⟨u_r\omega {⟩}_{z,\theta ,t}-\nu {\partial }_r⟨\omega {⟩}_{z,\theta ,t})$$

$$N{u}_{\omega }=\frac{J^{\omega }}{J_{lam}^{\omega }}$$

$$C_f\equiv \frac{2\pi {\tau }_{w,i}}{\rho d^2|{\omega }_i-{\omega }_o{|}^2}=2\pi N{u}_{\omega }{J}_{lam}^{\omega }(\nu R{e}_i{)}^{-2}$$

$$u_r^{+}=\frac{1}{\kappa }\log (y/k)+B=\frac{1}{\kappa }\log (y/k_s)+\tilde{B}$$

$$\Delta u^{+}=\frac{1}{\kappa }\log (k_s^{+})+A-\tilde{B}.$$

$$(C_f/2{)}^{-1/2}=C_1\log ((C_f/2{)}^{1/2}R{e}_i)+C_2.$$

$$u_{\theta }(r,\theta ,z,t)=⟨\overline{u_{\theta }}{⟩}_{\theta }(r,z)+{\tilde{u}}_{\theta }(r,\theta ,z)+u_{\theta }^{\prime }(r,\theta ,z,t)$$

Full text

Here, we present direct numerical simulations (DNSs) of TC turbulence with sand grain roughness mounted on the inner cylinder. The model proposed by Scotti (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 analyse the global response of the system, expressed both by the dimensionless angular velocity transport $Nu_{\omega}$ and by the friction factor $C_{f}$. An increase in friction with increasing roughness height is accompanied with enhanced plume ejection from the inner cylinder. Subsequently, we investigate the local response of the fluid flow over the rough surface. The equivalent sand grain roughness $k_{s}^{+}$ is calculated to be $1.33k$, where $k$ is the size of the sand grains. 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 (VDI-Forschungsheft 361, 1933) data of sand grain roughness in pipe flow, regardless of the Taylor number dependent constants of the logarithmic layer.

1 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 (Grossmann et al., 2016). 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 $\Delta U^{+}\equiv(U_{s}-U_{r})/u_{\tau}$, where $\Delta 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)

[TABLE]

As usual, the superscript ‘+’ indicates a scaling in viscous units (e.g. length $y^{+}=yu_{\tau}/\nu$ and velocity $u^{+}=u/u_{\tau}$) and $u_{\tau}$ is the friction velocity, $u_{\tau}=\sqrt[]{\tau_{w}/\rho}$ with $\tau_{w}$ being the total tangential stress at the wall, and $\rho$ 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_{\tau}/U_{0})^{2}$, where $U_{0}$ is the centerline velocity (Pope, 2000). It has been found that for TC turbulence $\kappa$ 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;

[TABLE]

with $f_{1}(Re_{i})$ and $f_{2}(Re_{i})$ being unknown functions. The questions now are: $i)$ How does $\Delta u^{+}$ depend on the parameters that characterize the surface geometry. $ii)$ Can $\Delta u^{+}$ be generalized to other flows.

Although many parameters influence the Hama roughness function $\Delta 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}^{+}\lesssim 5$, the rough wall appears to be hydrodynamically smooth and the roughness function $\Delta u^{+}$ goes to zero. For $k_{s}^{+}\gtrsim 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}^{+}\approx 5$, where roughness effects suddenly result in an inflectional increase of $\Delta 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 ($\approx 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, 2017, 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 $Nu\propto Ta^{0.5}$, although for higher $Ta$ 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 $Nu\propto Ta^{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, 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 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 large-eddy simulations in channel flow (Yuan & Piomelli, 2014). Taylor numbers in our DNS range from $Ta=1.0\times 10^{7}$($Re_{\tau}=82$) to $Ta=1.0\times 10^{9}$($Re_{\tau}=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 $\Delta 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.

2 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 $\eta=r_{i}/r_{o}$. For this research, we set $\eta\approx 0.714$, in accordance with the experimental T3C setup (Huisman et al., 2013), and previous simulations (Zhu et al., 2017). $\Gamma=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, $\Gamma\approx 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_{\theta}/d=(r_{i}2\pi/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_{\theta}/(d/2)=5.24$ and $L/(d/2)\approx 4.0$.

To maintain convenient boundary conditions, we solve the Navier–Stokes (NS) equations in a reference frame rotating with the inner cylinder ($\omega_{i}\textbf{e}_{z}$) The NS equations in that reference frame become

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with no-slip boundary conditions $u_{\theta}|_{r=r_{i}}=0$, $u_{\theta}|_{r=r_{o}}=r_{o}(\omega_{o}-\omega_{i})$. Equation (6) expresses the incompressible restriction. Hatted symbols indicate the respective dimensionless variables, with $\mathbf{u}=r_{i}|\omega_{i}-\omega_{o}|\hat{\mathbf{u}}$, $r=d\hat{r}$ and $t=\frac{d}{r_{i}|\omega_{i}-\omega_{o}|}\hat{t}$. We define $\frac{f(\eta)}{Ta^{1/2}}=Re^{-1}$ where $f(\eta)$ is the so-called ‘geometric Prandtl’ number (Eckhardt et al., 2007):

[TABLE]

here $f(0.714)\approx 1.23$ and $Ta$ is the Taylor number, which is a measure of the driving strength of the system,

[TABLE]

Note that the pressure in the equations above represents the ‘reduced pressure’ that incorporates the centrifugal term; $\hat{P}=p^{\prime}-\frac{\omega_{i}^{2}d^{2}\hat{r}^{2}}{2r^{2}_{i}|\omega_{i}-\omega_{o}|^{2}}\mathbf{e}_{r}$ with $p=\rho r_{i}^{2}|\omega_{i}-\omega_{o}|^{2}p^{\prime}$ 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. (3) and (4) gives the Coriolis force, with $Ro^{-1}$ being the inverse Rossby number

[TABLE]

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^{\omega}$, calculated by:

[TABLE]

with the laminar flux given by $J^{\omega}_{lam}=2\nu r_{i}^{2}r_{o}^{2}\frac{\omega_{1}-\omega_{2}}{r_{2}^{2}-r_{1}^{2}}$ where $\nu$ is the kinematic viscosity and $\langle.\rangle_{z,\theta,t}$ indicates averaging over the spatial directions $z$, $\theta$ and time $t$. For incompressible flows, it can be derived from the NS equations that $J^{\omega}$ is conserved in the radial direction, $\partial_{r}J^{\omega}=0$ (Eckhardt et al., 2007). In both cases, the values are made dimensionless by their respective laminar, conducting, values. For TC Nusselt becomes:

[TABLE]

The angular velocity flux $J^{\omega}$ can easily be written in terms of the torque $\mathcal{T}$ on any of the cylinders: $J^{\omega}=\mathcal{T}(2\pi L\rho)^{-1}$ with $\rho$ being the fluid density. Consequently, the shear stress on the inner cylinder $\tau_{w,i}$ is related to the angular velocity flux by $\tau_{w,i}=\frac{\rho J^{\omega}}{r_{i}^{2}}$.

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 $Nu_{\omega}$ 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 define $C_{f}\equiv G/Re_{i}^{2}$, where $G$ is the dimensionless torque $G=\mathcal{T}/(\rho\nu^{2}L)$ and $Re_{i}$ is the inner cylinder Reynolds number $Re_{i}=\frac{r_{i}\omega_{i}d}{\nu}$. The translation between $Nu_{\omega}$ and $C_{f}$ is straightforward, namely

[TABLE]

Note that one can also define $C_{f}\equiv 2\tau_{w}/\rho(r_{i}\omega_{i})^{2}=\frac{(1-\eta)^{2}}{\pi\eta^{2}}\frac{G}{Re_{i}^{2}}$, which is different from eq. (12) by a factor which depends on the radius ratio $\eta$ (Lathrop et al., 1992).

3 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\times 2k$, each tile containing exactly $1$ ellipsoid. The height $L$ is slightly varied ($0.85d\pm 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}\Delta\theta$ and $\Delta z$, where $r_{i}\Delta\theta$, and $\Delta z$ are random uniform translations from the center of the $2k\times 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 $\Delta r=-0.5k$ from $r=r_{i}$. It can be seen from figure 1($a$) that part of the cylinder ($\approx 15\%$) is not covered by rough elements.

3.2 Numerical method

The NS equations are spatially discretized by using a central second-order finite-difference 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 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) ($\frac{u\Delta t}{\Delta 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 $\mathbf{f}$ to the momentum equations (3) - (5). 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., 2016, 2017, 2018b).

Simulations are run until they become statistically stationary, such that Nusselt $Nu_{\omega}$ number remains constant to within $\approx 1\%$, which requires $\hat{t}\approx 200$ . Thereafter, we gather statistics until the results converge, which requires $\hat{t}\approx 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 $Nu_{\omega}$ remains constant with increasing grid resolution in all directions and presented along with the results in Table 1.

4 Results

4.1 Roughness function

Figure 3 (left column) presents the streamwise (i.e. azimuthal) velocity profiles $u^{+}=\langle(u(r)-u(r_{i}))\rangle_{z,\theta,t}/u_{\tau}$ (solid) and angular velocity profiles $\omega^{+}=\langle u(r_{i})-r_{i}u(r)/r\rangle_{z,\theta,t}/u_{\tau}$ (dashed) versus the wall normal distance $y^{+}=r^{+}-h_{m}^{+}$, where $\langle.\rangle_{z,\theta,t}$ indicates averaging over the streamwise and spanwise directions and in time, and $h_{m}^{+}$ is the mean roughness 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 $\kappa=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_{\tau}$. However, this asymptotic value is experimentally observed at very high shear rates of $Ta=O(10^{12})$ (Huisman et al., 2013), much higher than can be obtained by the present DNS. The logarithmic profiles of angular velocity $\omega^{+}$ have a slope that is closer to the $\kappa=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 $\omega^{+}$.

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 $\Delta u^{+}=u^{+}_{s}-u^{+}_{r}$ and $\Delta\omega^{+}=\omega^{+}_{s}-\omega^{+}_{r}$. For lower $Ta$, there is a small but observable difference between $\Delta u^{+}$ and $\Delta\omega^{+}$, see figure 3($b$), whereas for the higher $Ta$, 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 $Ta$. Care is taken to ensure overlap for varying $Ta$ and similar $k_{s}^{+}$, to study the $Ta$ dependence of $\Delta u^{+}$ and $\Delta\omega^{+}$. However, despite the varying $Ta$ 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 $\lesssim 0.2\Delta u^{+},0.2\Delta\omega^{+}$.

A comparison with the findings of Nikuradse (1933) can be carried out by scaling the 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\gg k\gg\delta_{\nu}$), whereas the velocity profile is also independent of the system outer length scales ($y\ll d$) i.e. the overlap argument (Pope, 2000). The gradient of the velocity profile becomes; $\frac{d\langle U\rangle}{dy}=\frac{u_{\tau}}{y}\Phi(y/k)$, where $\Phi(y/k)$ is a universal function that will go to $1/\kappa$. Integration then gives;

[TABLE]

where $B$ is a constant and $y=r-h_{m}$. $\tilde{B}$ is the Nikuradse constant. The roughness function in the fully rough regime, (i.e. the fully rough asymptote), is obtained by subtracting (13) from the smooth wall profile $u_{s}^{+}=1/\kappa\log{(y^{+})}+A$ and rescaling it to overlap with the Nikuradse data:

[TABLE]

In figure 4($a$), the blue solid line is the fully rough asymptote, with $\kappa,A,\tilde{B}$ as found in pipe flow (Pope, 2000). The green solid line is the fully rough asymptote as obtained from our simulations. $\kappa_{u}^{-1}=1.22$ and $A_{u}=8.0$ are taken from the fit of the smooth wall simulation at identical $Ta$ as the fully rough cases (figure 3$g$). The fits are in the domain $y^{+}=[150,500]$, as there the slope becomes approximately constant (figure 3$h$). $\tilde{B}$ is plotted in figure 5($b$), where we find that $\tilde{B}\approx 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.

In parallel, we analyse the behavior of $\Delta\omega^{+}$ 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 $Ta$ as the fully rough cases, see figure 3($g$). We find $\kappa_{\omega}^{-1}=2.17$ ($A_{\omega}=3.7$), close to the asymptotic value $\kappa^{-1}=2.44$. Although the differences are marginally, $\Delta\omega^{+}$ fits to the Nikuradse data slightly better than $\Delta 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 $\Delta u^{+}$ at similar $k_{s}^{+}$, for varying $Ta$, see figure 4. Also the similarity between $\Delta u^{+}$ and $\Delta\omega^{+}$ 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 $Ta$ ($=10^{12}$ (Huisman et al., 2013)), where $\kappa=0.4$, then are expected to also give a collapse to the Nikuradse data in the fully rough regime.

4.2 Global response

Figure 6($a$) presents the friction factor $C_{f}$ versus the dimensionless roughness height $k/d$ for varying $Re_{i}$. For increasing $Ta$ numbers, the friction factor decreases (e.g. figure 7.23 in Pope (2000)), an indication that pressure is not yet dominant. For constant $Ta$, as expected, the friction factor increases for increasing roughness height. Figure 6($b$) then presents the global response in terms of $Nu_{\omega}$. We observe an increase in $Nu_{\omega}$ for increasing $Ta$, corresponding to the increased transport of the angular velocity that is due to the increased turbulent mixing. Higher roughness leads to increased $J^{\omega}$ as compared to the smooth wall at the same $Ta$, which also relates to a higher intensity of the turbulent mixing (the $r^{3}(\langle u_{r}\omega\rangle_{z,\theta,t})$ term of equation (10)) 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:

[TABLE]

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 (12), 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 $Nu_{\omega}\propto Ta^{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 $Nu_{\omega}$, and only the inner cylinder is covered with roughness.

4.3 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_{\theta}(r,\theta,z)$. Figure 8($a-c$) exhibit the snapshots for $Ta=5.0\times 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 (Grossmann et al., 2016). 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\approx 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 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\approx 1.3$) plumes have very strong radial velocity components $u_{r}$. From the first term on the right-hand side of (10), $r^{3}(\langle u_{r}\omega\rangle_{z,\theta,t})$, we then directly see that they strongly contribute to $Nu_{\omega}$. 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 $Nu_{\omega}$. 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 ($\eta=(\nu^{3}/\epsilon)^{1/4}$), namely $\eta/d=7.1\times 10^{-3}$ for the smooth wall case BS and $\eta/d=6.5\times 10^{-3}$ for the highest roughness case B5. Note that the decrease in the Kolmogorov scale $\eta$ is only small, since $\eta/d\propto(\epsilon d^{4}/\nu^{3})^{-1/4}$. For TC flow, the volume averaged dissipation rate $\epsilon$ is related to the angular velocity transport $Nu_{\omega}$ with: $\epsilon=\nu^{3}d^{-4}\sigma^{-2}(Nu_{\omega}-1)Ta+\epsilon_{lam}$, where $\epsilon_{lam}$ is the laminar volume averaged dissipation rate, $d$ is the gap width of the setup and $\sigma$ is a geometric parameter (Eckhardt et al., 2007). As such, we see that $\eta/d\propto Nu_{\omega}^{-1/4}$ only.

Figure 8($d-f$) presents snapshots of a flow in the ultimate turbulent state at $Ta=1.0\times 10^{9}$ (Grossmann et al., 2016). Although less pronounced than for $Ta=5.0\times 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 8$f$), in comparison to a lower roughness case (D2, figure 8$e$). 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, $\eta/d=2.7\times 10^{-3}$ for the smooth wall case DS and $\eta/d=2.1\times 10^{-3}$ for the highest roughness case D5.

4.4 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 $\langle\overline{u_{\theta}}\rangle_{\theta}\approx 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 $\overline{u_{\theta}}$, 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_{\theta}^{+}$. 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 horizontal lines. Small recirculation zones (closed isolines) are detected in open regions behind high roughness elements. To quantify the degree of roughness-induced velocity disturbances, we apply a triple decomposition to the instanteneous azimuthal velocity $u_{\theta}(r,\theta,z,t)$ (Pokrajac et al., 2007):

[TABLE]

where $u^{\prime}_{\theta}(r,\theta,z,t)=u_{\theta}(r,\theta,z,t)-\overline{u_{\theta}}(r,\theta,z)$ is the temporal fluctuation and $\tilde{u}_{\theta}(r,\theta,z)=\overline{u_{\theta}}(r,\theta,z)-\langle\overline{u_{\theta}}\rangle_{\theta}(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 $\tilde{u}$). Note that by applying the triple decomposition to the full NS equations, one will recover the related form-induced stress tensor $\langle\overline{\tilde{u}_{i}\tilde{u}_{j}}\rangle_{\theta}$. However, here we only discuss $\tilde{u}$. The root mean square of the form-induced fluctuations $\sqrt{\tilde{u}^{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 $\theta$. 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 $\sqrt{\langle\tilde{u}^{+2}\rangle_{\theta,z}}=0.01\langle\bar{u}^{+}\rangle_{\theta}$ and find $h_{r}=3.70k$ ($h_{r}=2.78k_{s}$). This value agrees well with a roughness sublayer height of $2\lesssim k_{s}\lesssim 5$, as typically found in other canonical flows (Pokrajac et al., 2007).

4.5 Radial profiles

We plot the time and azimuthally averaged azimuthal velocity for $Ta=5.0\times 10^{8}$ in figure 12($a$). The roughness covers the inner cylinder, i.e. $(r-r_{i})/d\leq 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\pm 0.002$ and varying Taylor number.

Figure 13($a$) shows the double-averaged radial profiles $(u_{\theta})_{rms}^{+}=\langle\langle\overline{u^{2}_{\theta}}\rangle_{\theta,z}-\langle\overline{u_{\theta}}\rangle_{\theta,z}^{2}\rangle_{\theta,z}^{1/2}/u_{\tau}$ for constant Taylor number $Ta=1.0\times 10^{9}$. For the smooth wall case, the peak of the rms of the azimuthal velocity is located at $u_{\theta,rms}|_{max}\approx 10y^{+}$. This agrees with the smooth case in (ostilla16) for TC flow, but is closer to the solid boundary than that is found in channel flow; $u_{\theta,rms}|_{max}\approx 12y^{+}$. We observe a slight decrease in the viscous scaled turbulence intensity $(u_{\theta})_{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^{\omega}(y^{+})$ is conserved in the radial direction, see section 4.2. However, the individual components – i.e. the viscous $J_{\nu}^{\omega}(y^{+})$ and Reynolds stress $J_{u_{r}\omega}^{\omega}(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^{\omega}(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^{\omega}\approx J_{\nu}^{\omega}$; this is expected, since there the gradient of the mean streamwise velocity is maximum and the wall normal component 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^{\omega}\approx J_{u_{r}\omega}^{\omega}$. 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.

5 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 $Ta=1.0\times 10^{7}$ ($Re_{\tau}=82$) up to $Ta=1.0\times 10^{9}$ ($Re_{\tau}=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, $\tilde{u}_{\theta}=\langle u_{\theta}\rangle_{t}-\langle u_{\theta}\rangle_{t,\theta}$. 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 $\sqrt{\langle\tilde{u}^{+2}\rangle_{\theta,z}}=0.01\langle\bar{u}^{+}\rangle_{\theta}$ 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 $\Delta u_{\theta}^{+}$. 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 $\Delta u_{\theta}^{+}=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_{\theta}^{+}=f_{1}(Re_{i})\log(y^{+})+f_{2}(Re_{i})-\Delta u_{\theta}^{+}$.

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.

Acknowledgements

We wish to express our gratitude to Jelle Will for his help with various illustrations. This project is funded by the Priority Programme SPP 1881 Turbulent Superstructures of the Deutsche Forschungsgemeinschaft. We also acknowledge PRACE for awarding us access to Marconi, based in Italy at CINECA under PRACE project number 2016143351. We also acknowledge PRACE for awarding us access to MareNostrum based in Italy at the Barcelona Supercomputing Center (BSC) under PRACE project number 2017174146. This work was partly carried out on the national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organization for Dutch education and research.

Appendix A