2.1. Governing equations

In the present work, the incompressible, fully turbulent flow of PL fluids in the channel geometry is considered. The governing equations are the continuity and Navier–Stokes equations, with the general viscous stress term:

(2.1) \begin{align}&\qquad\qquad\quad \boldsymbol{\nabla \cdot} \boldsymbol{{u}} = 0, \end{align}

(2.2) \begin{align}& \frac {\partial \boldsymbol{{u}}}{\partial t} + \boldsymbol{{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{{u}}= - \frac {1}{\rho }\boldsymbol{\nabla} p+\frac {1}{\rho } \boldsymbol{\nabla \cdot} \boldsymbol{\tau }, \end{align}

where $\boldsymbol{{u}}$ , $\rho$ , $p$ and $\boldsymbol{\tau }$ denote respectively the velocity vector, density, pressure and the viscous stress tensor. To maintain a constant mass flow rate, a fixed pressure drop along the streamwise direction has been also imposed on the system.

2.2. Rheology

Equation (2.2) is complemented by the constitutive model of a GNF, i.e. $\boldsymbol{\tau } = 2 \mu (\dot {\gamma }) \unicode{x1D64E}$ , where $\mu$ is the viscosity function, $\unicode{x1D64E}={1}/{2}(\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} ^{\text{T}} \boldsymbol{u})$ is the rate-of-strain tensor and $\dot {\gamma }\equiv \sqrt {2\unicode{x1D64E}:\unicode{x1D64E}}$ is the deformation rate, which is associated with the second invariant of $\unicode{x1D64E}$ . The power-law model is given as a viscosity function of the following form:

(2.3) \begin{equation} \mu ( \dot {\gamma } ) = K \dot {\gamma }^{n-1}, \end{equation}

where $K$ is the consistency parameter and $n$ is the power-law index, which measures the degree to which a fluid is shear-thickening $(n \gt 1)$ or shear-thinning $(n \lt 1)$ . Needless to say, when $n = 1$ , the PL model is reduced to Newtonian rheology. The simulations carried out here consider both the Newtonian and PL rheologies.

2.3. Scaling

When the fluid viscosity varies in both space and time, the appropriate reference viscosity to be used in the definition of the Reynolds number is not obvious. Numerous choices are available, including zero shear-rate (Beris & Dimitropoulos Reference Beris and Dimitropoulos1999), average, mean centreline and mean wall viscosity. Using the latter, the Reynolds number would assume the following expression:

(2.4) \begin{equation} Re_w = \frac {V h}{\nu _w}, \end{equation}

where $h$ and $\nu _w$ stand for half the channel height and the wall kinematic viscosity, which is given by

(2.5) \begin{equation} \nu _w = \frac {\tau _w}{\rho } \left ( \frac {K}{\tau _w} \right )^{1/n}. \end{equation}

In many works such as Rudman et al. (Reference Rudman, Graham, Blackburn and Pullum2002, Reference Rudman, Blackburn, Graham and Pullum2004), Rudman & Blackburn (Reference Rudman and Blackburn2006), Guang et al. (Reference Guang, Rudman, Chryss, Slatter and Bhattacharya2011) and Rudman & Blackburn (Reference Rudman and Blackburn2012), the Reynolds number based on (2.4) and (2.5) was employed. However, as we will show in § 3.2, this choice is not adequate since it does not reflect the desirable underlying effects.

Choosing a characteristic viscosity of the form $K(\textit{V/h})^{n-1}$ leads to the following power-law expression for the Reynolds number:

(2.6) \begin{equation} Re_n = \frac {\rho V^{2-n}h^n}{K}. \end{equation}

It is worth observing that Metzner & Reed (Reference Metzner and Reed1955) developed an expression for the Reynolds number for power-law fluids in a tube by matching the friction factor in the laminar case, where a factor of $8\,[n/(6n+2) ]$ was found to apply. In the case of the plane channel, a similar analysis would lead to a factor of $6\,[n/(4n+2) ]$ (Thompson & Oishi Reference Thompson and Oishi2021). However, Anbarlooei et al. (Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018) have shown that (2.6) is more effective in expressing the friction in turbulent flows of power-law fluids because more universal coefficients were obtained, in contrast to coefficients which are more dependent on the flow exponent ( $n$ ). This matter will be further discussed in § 3.2.

As well as in the Newtonian case, the friction velocity ( $u_\tau =\sqrt {\tau _w / \rho }$ ) has been used for the velocity scaling $U^+=U/u_\tau$ . By considering the momentum equation in the viscous layer, one can find the viscous length scale for the PL fluid as

(2.7) \begin{equation} \delta _{\nu } = u_\tau \left ( \frac {K}{\tau _w} \right )^{1/n}, \end{equation}

which recovers $U^+=y^+=y/\delta _\nu$ near the wall. Based on (2.7), the friction Reynolds number $Re_{\tau }$ is defined as

(2.8) \begin{equation} Re_{\tau } = \frac {h}{\delta _\nu } = \left ( \frac {\rho u_\tau ^{2-n} h^n}{K} \right )^{1/n}. \end{equation}

2.4. Numerical method

The unsteady, incompressible, three-dimensional equations of motion were solved with the open source code OpenFOAM (Greenshields Reference Greenshields2015). The standard solver nonNewtonianIcoFoam used in the simulations, with some minor changes, is based on the finite-volume method of Ferziger & Peric (Reference Ferziger and Peric1999) and the predictor–corrector PISO algorithm of Issa (Reference Issa1986). The second-order set-up in space and time was used for all simulations. The preconditioned bi-conjugate gradient method with incomplete-Cholesky preconditioner (DIC) was used to solve the linear systems with a local accuracy of $10^{-7}$ for all dependent variables at each time step. A Courant–Friedrichs–Lewy (CFL) number equal to $0.5$ was used through all simulations.

Although it is common to use higher-order methods (i.e. spectral methods) for DNS investigations, OpenFOAM has demonstrated its applicability in different DNS studies (e.g. Komen et al. Reference Komen, Shams, Camilo and Koren2014; Vo et al. Reference Vo, Kronenburg, Stein and Hawkes2016; Zheng et al. Reference Zheng, Rudman, Singh and Kuang2019; Karahan et al. Reference Karahan, Ranjan and Aidun2023). It has been observed before (Sprague Reference Sprague2010) that to achieve similar accuracy as spectral methods, OpenFOAM requires at least eight times more grid points (two times more in each direction). However, using the finite volume approach allows us to study more complex geometries. This advantage suppresses the former shortcoming in the simulation of the flow with rough surfaces as in the present work.

2.5. Geometries, grid and domain size

A sketch illustrating the geometry of the channel, the flow orientation and the coordinate system is presented in figure 2. In the present notation, the dimensions of the channel are denoted by ( $L_x h, 2h, L_z h$ ), which have been discretised with a grid of $N_x \times N_y \times N_z$ points. The grid points were equally spaced in the streamwise ( $x$ ) and spanwise ( $z$ ) directions, while in the normal direction to the walls ( $y$ ), distributed with a simple geometrically expanding rate.

Figure 2.

Coordinate system for the plane channel geometry.

Figure 3 shows the rod-roughened channel used to study the effect of roughness on the turbulent flow of the power-law fluids. Square rods with a height of $1.7\,\%$ of the channel size $2h$ were used to make both walls of the channel rough. The pitch-to-height ratio for the roughness elements was $8$ , which corresponds to a ‘k-type’ roughness. This configuration is the same as that used by Ashrafian et al. (Reference Ashrafian, Andersson and Manhart2004); Krogstad et al. (Reference Krogstad, Andersson, Bakken and Ashrafian2005) to study the turbulent flow of Newtonian fluid over rough surfaces. The computational grid has been divided into two parts: the near-wall region and the central region. The latter one has been meshed with the same structure as the smooth channel flow, with equally spaced grid points in the streamwise and spanwise directions, and gradually expanding points in the normal direction. However, in the near-wall region, and up to three heights of the roughness element from the wall, equally spaced grid points were used. The computational mesh is also illustrated in figure 3.

Figure 3.

Rod-roughened channel and its computational grid.

Assuming that the features that are smaller than Kolmogorov’s micro-scale $\eta$ will be dissipated by the viscosity, it is an established practice to distribute grid points considering the local variation of $\eta$ (Coleman & Sandberg Reference Coleman and Sandberg2010). Anbarlooei et al. (Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018) extended Kolmogorov’s micro-scales to GNF. Considering that these scales depend only on the turbulent dissipation $\varepsilon$ and on the rheology, it is possible to obtain the following relations for the power-law Kolmogorov’s length and velocity scales (Anbarlooei et al. Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018):

(2.9) \begin{align}& \rho ^3 \varepsilon ^2 \eta ^{2(n+1)} - K^3 \varepsilon ^n = 0, \end{align}

(2.10) \begin{align}&\quad \rho u_\eta ^{2(n+1)} - K \varepsilon ^n = 0. \end{align}

Equations (2.9) and (2.10) reduce to

(2.11) \begin{equation} \eta \sim \left ( \frac {K^3}{\rho ^3 \varepsilon ^{2-n}} \right )^{\frac {1}{2(n+1)}} ; \qquad u_\eta \sim \left ( \frac {K \varepsilon ^n}{\rho } \right )^{\frac {1}{2(n+1)}}. \end{equation}

Therefore, setting $n=1$ recovers the Newtonian Kolmogorov’s micro-scales. Assuming $\varepsilon \sim V_L^3/L$ , it is possible to obtain $\eta /L$ and $u_\eta /V_L$ , where $V_L$ and $L$ are the large (energetic) velocity and length scales. For the PL fluid, this will lead to

(2.12) \begin{equation} \frac {\eta }{L} \sim Re_{n,L}^{-\frac {3}{2(n+1)}} ; \qquad \frac {u_\eta }{V_L} \sim Re_{n,L}^{-\frac {1}{2(n+1)}}, \end{equation}

where $Re_{n,L}$ is the power-law Reynolds number based on $L$ and $V_L$ .

This extension of K41 allows us to check the grid requirements mentioned before for the GNF. This issue has not been adequately addressed in previous DNS studies of non-Newtonian fluids, on both viscoplastic and viscoelastic ones (i.e. Rudman et al. Reference Rudman, Blackburn, Graham and Pullum2004; Rudman & Blackburn Reference Rudman and Blackburn2006; Singh et al. Reference Singh, Rudman and Blackburn2017a ; Zheng et al. Reference Zheng, Rudman, Singh and Kuang2019), using the simulations of the turbulent power-law fluid, showed that by decreasing the index $n$ , the intensity of the structures decreases as well. Therefore, it has been inferred that Kolmogorov’s length scale will also increase. However, examining (2.12) makes clear that, while keeping the Reynolds number constant, decreasing $n$ results in finer micro-scales. Hence, one can use (2.9) to verify the selected grid spacing. Based on this argument, and also on the requirement of finer grids in OpenFOAM for obtaining comparable accuracy as spectral methods, in this work, more grid points were used to simulate GNF than the Newtonian simulations. This is in contrast to the grid used by Zheng et al. (Reference Zheng, Rudman, Singh and Kuang2019), where the adapted resolution is not suitable even for Newtonian simulations. We believe that one of the reasons for the reported discrepancy of Zheng et al. (Reference Zheng, Rudman, Singh and Kuang2019) between OpenFOAM simulations and spectral (Semtex) ones is this inferior grid resolution.

All but one of the simulations in the current study were carried out for $Re_{\tau } = 395$ , with the same computational grid as that used by Moser, Kim & Mansour (Reference Moser, Kim and Mansour1999) for the Newtonian case ( $256 \times 192 \times 192$ ). Finer grids were also examined for the Newtonian and power-law simulations with $n=0.5$ . In these cases, it was inferred that grid dependence is negligible at this resolution for the first and second moments. However, it is also evident that by decreasing the power-law exponent, the discrepancy between coarse and finer grids grows. In addition, following the guidelines of Coleman & Sandberg (Reference Coleman and Sandberg2010), it was assured that in each simulation, the computational grid has at least one grid point before $y^+=1$ and ten points before $y^+=10$ .

Previous numerical experiments (Anbarlooei, Cruz & Freire Reference Anbarlooei, Cruz and Freire2014) have shown that by decreasing the power-law index $n$ , the required length of the domain to obtain statistically converged results increases. Anbarlooei et al. (Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018) found that a streamwise domain size twice as large as that usually employed in the Newtonian simulations is necessary for the power-law fluids with indices down to $n = 0.5$ . In this regard, the simulations in the present study were carried out on a domain of size $(4 \pi h \times 2h \times \pi h)$ . As a result, two times more grid points were used in the streamwise direction when compared with the Newtonian simulations of Moser et al. (Reference Moser, Kim and Mansour1999). These domain lengths were also verified using two-point velocity fluctuations correlations for streamwise and spanwise directions (Anbarlooei et al. Reference Anbarlooei, Cruz and Freire2014). To summarise, grids with at least $1.88 \times 10^7$ cells and up to $1.5 \times 10^8$ cells (grid dependence study) were used in the present work. For comparison, the maximum number of cells used by Zheng et al. (Reference Zheng, Rudman, Singh and Kuang2019) is $8 \times 10^6$ .

2.6. Fully developed flow and averaging process

The fully developed Newtonian turbulent flow at the same $Re_{\tau }$ was used to initialise the velocity field for the non-Newtonian simulations, while the flow index $n$ was gradually changed from 1 to the target value.

To ensure the fully developed flow, the averaging process was started after at least 30 recirculations of the flow inside the domain, i.e. ${t V}/{L_x h} \gt 30$ . Two other criteria were also adopted. The averaging process was only started when the fluctuation in the bulk velocity was less than 1 % of its average over 10 recirculations of the flow, and also the difference between the friction factor obtained by the simulation and that calculated by the theory developed by Anbarlooei et al. (Reference Anbarlooei, Cruz, Ramos and Freire2015b , Reference Anbarlooei, Cruz, Ramos, Santos and Freire2017, Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018) was less than 5 %. Then, the turbulence statistics were gathered for another 100 recirculations of the flow in the domain.

The summary of the DNS studies carried out for the present work is provided in table 1.

Table 1.

Simulation parameters for the smooth and rough channel flows of power-law fluids studied in the present work. $n$ is the power-law index.

3.1. Power-law simulations

The first set of results presented here is obtained for the PL rheology with $n=2.0, \, 1.0, \, 0.75 \text{ and } 0.5$ , for the smooth wall and $Re_\tau = 395$ . The instantaneous velocity contours near the wall ( $y^+ \sim 15$ ) are shown in the left part of figure 5. These contours show that keeping $Re_{\tau }$ fixed, decreasing $n$ results in larger (longer) structures near the wall. The same behaviour has been observed in the pipe geometry (Rudman et al. Reference Rudman, Blackburn, Graham and Pullum2004; Rudman & Blackburn Reference Rudman and Blackburn2006; Singh et al. Reference Singh, Rudman and Blackburn2017b ). This indeed indicates that longer domains are necessary to capture these structures and supports our previous observation (Anbarlooei et al. Reference Anbarlooei, Cruz and Freire2013, Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018). However, in the shear-thickening case ( $n=2.0$ ), the random and small-scale fluctuations (structures) are visible.

The right-hand side of figure 5 shows the instantaneous velocity contours in a plane normal to the walls. Here, we can notice a lower turbulent intensity in the shear-thinning cases than in the Newtonian simulation. For the shear-thickening case, the trend is inverse. These trends are similar to the effects of changing the Reynolds number in Newtonian flows and are aligned to what was reported by Pereira, Thompson & Mompean (Reference Pereira, Thompson and Mompean2019) in the context of viscoelastic fluids. Speaking qualitatively, shear-thinning fluids exhibit signatures of low-Reynolds-number Newtonian flows, while the shear-thickening case appears equivalent to higher Reynolds number flows.

Figure 6 shows the turbulent structures in half of the channel for the PL simulations, obtained using the Q-criterion. To ensure the same threshold value for every case, $Q$ was non-dimensionalised by $u_\tau / \delta _\nu$ . Following the ideas of Pereira et al. (Reference Pereira, Mompean, Thais and Thompson2017), we show elliptical-like ( $Q\lt 0$ ), as well as hyperbolic-like ( $Q\gt 0$ ) structures. The iso-surfaces in figure 6(a) are for the values of $Q = 0.1$ and $Q = -0.1$ in figure 6(b), the colours indicate the velocity magnitude. A movie demonstrating the contours of the PL cases of figure 6(a) in the progression of the flow can be accessed in the Supplementary material, which also demonstrates the iso-surfaces with a larger threshold of $Q = 1.0$ .

From this perspective, we also note that by decreasing $n$ , elliptical and hyperbolic structures become more sparse, larger and longer, which resembles the qualitative tendency obtained by decreasing the Reynolds number in Newtonian turbulence. This fact was also observed by Pereira et al. (Reference Pereira, Thompson and Mompean2019) with viscoelastic DNS of the FENE-P model. There, the increase in elasticity and shear-thinning behaviours led to a turbulence state comparable with the Newtonian turbulence in the transition range. Since FENE-P fluids are shear-thinning-extensional-thickening materials, the qualitative alignment of the turbulence of power-law shear-thinning fluids (where elasticity is absent) with the viscoelastic turbulence indicates that viscous effects can be predominant in some observed phenomena in viscoelastic drag-reducing turbulent flows.

Figure 6.

Three-dimensional structures represent iso-surfaces of vortical regions defined by the Q-criterion. For the four simulations, the values of Q were normalised by $u_\tau / \delta _\nu$ . The contours in panels (a) and (b) are for the threshold of $Q = 0.1$ and $Q = -0.1$ , respectively. The colours indicate the velocity magnitude. A movie demonstrating the iso-surfaces in the progression of the flow can be accessed in the Supplementary material.

This topic warrants further investigation and will be the focus of our future work. However, from figure 5, it is evident that although the sizes of the structures vary with different values of the power-law index, they exhibit a certain scale similarity with the inverse of the corresponding Kolmogorov scale. This observation suggests that, for the current settings ( $Re_\tau = 395$ and $y^+ = 15$ ), the average transverse number of structures ( $N_s$ ) is proportional to the dimensionless Kolmogorov scale for each power-law index, as $N_s \propto \eta$ . This phenomenon is illustrated in figure 7, where the sizes of the contours shown in figure 5 are scaled according to the inverse of the Kolmogorov scale for each power-law index as

(3.1) \begin{align} L_{n=1} &= L_{n=2} \times \frac {\eta _{n=1}}{\eta _{n=2}} \approx 0.53, \nonumber\\[4pt] L_{n=0.75} &= L_{n=2} \times \frac {\eta _{n=0.75}}{\eta _{n=2}} \approx 0.36, \\[4pt] \nonumber L_{n=0.5} &= L_{n=2} \times \frac {\eta _{n=0.5}}{\eta _{n=2}} \approx 0.21. \end{align}

These figures show that the structures for different values of $n$ look qualitatively similar and confirm that the structure sizes remain consistent when scaled according to the inverse of the Kolmogorov length, indicating that simulations using $n=2.0$ (or even higher) might effectively capture certain features of turbulent flows with reduced computational effort.

Figure 7.

Velocity contours for the Newtonian and power-law simulations at $y^+ = 15$ , scaled according to the inverse of the Kolmogorov length.

Although such a strongly shear-thickening fluid studied here ( $n=2.0$ ) does not have any known physical counterpart in practice, it can provide a means to understand what is happening in shear-thinning-extensional-thickening viscoelastic fluids, in the case where turbulence induces an extensional motion. In addition, based on the previous observations drawn here, it could be used as a toy model to study the effect of high Reynolds turbulent effects with less computational effort. As (2.11) indicates, Kolmogorov’s length scale of shear-thickening fluids is larger than the Newtonian case. For example, for the case of $n=2$ and $Re_n=5000$ , the simulation needs at least $1/8$ of the grid points in each direction in comparison with a Newtonian simulation at the same $Re_n=5000$ . This would ideally result in at least $8^4$ times less computational load in the isotropic turbulence simulations. This topic needs further investigation and is the subject of our future work.

The mean axial velocity profiles of the simulations are presented in figure 8. In the near-wall region (viscous sub-layer), the power-law behaviour is in agreement with the classical linear velocity profile ( $U^+=y^+$ ). Away from the wall and in the logarithmic region, decreasing the PL index ( $n$ ) results in a higher mean velocity, indicating less friction and higher maximum velocity values at the channel centre. This behaviour is already expected as the flow rate increases in shear-thinning fluids (drag reduction). In contrast, the flow rate decreases in the shear-thickening case, resulting in lower mean velocities. We notice that when $n=0.75$ , the inclination in the log-law layer has not changed significantly with respect to the Newtonian case, which is a fingerprint of a low drag reduction (LDR) subregime. However, for $n=0.5$ , the shift in the MVP is accompanied by a mild change in the inclination, which can be a sign of a high drag reduction (HDR) effect. The MVP obtained with $n=2.0$ shows a departure from the Newtonian behaviour with a higher inclination and a lower discrepancy closer to the centre of the channel.

Next, we examine the functional form of the MVP. Figure 9(a) inspects the existence of a logarithmic region in the profile. As observed, even if one considers the presence of this region at such a low Reynolds number, its extent is very limited. Also, by decreasing $n$ , the Von Karmán constant decreases too; however, the shear-thickening case does not follow this trend. Figure 9(b) examines the presence of a power-law-type region ( $U^+=\alpha y^{+ \beta }$ ) in the velocity profile. As the fluid becomes more shear-thinning, the power-law relation performs better and is applicable to a larger portion of the velocity profile. In § 3.3, these velocity profiles are studied in more detail.

Figure 8.

Mean axial velocity profiles from the four smooth channel simulations with $Re_\tau = 395$ .

Figure 9.

Diagnostic quantities for the log-law $y^+ {\textrm d}U^+/{\textrm d}y^+$ and the power-law $y^+/U^+ {\textrm d}U^+/{\textrm d}y^+$ velocity profiles. The constant value in panel (a) represents the von Karmán constant in the log-law, while that in panel (b) represents the exponent in the power-law velocity profile.

The averaged turbulence intensities $\overline {u'u'}$ , $\overline {v'v'}$ and $\overline {w'w'}$ , Reynolds shear stress $\overline {u'v'}$ , turbulent kinetic energy $k$ , and viscous shear stress $\overline {\tau _{xy}}$ are plotted in figure 10. Here, the Reynolds average operator is denoted by the upper bar $\overline {q}$ , where $q$ is the averaged quantity. Figure 10(a) shows that decreasing $n$ increases the streamwise turbulent intensity $\overline {u'u'}$ . In contrast, the normal $\overline {v'v'}$ (figure 10 b) and spanwise $\overline {w'w'}$ (figure 10 c) intensities, and the Reynolds shear stress $\overline {u'v'}$ (figure 10 d) become higher by increasing $n$ , despite the strange behaviour of the case $n=2.0$ for the spanwise intensity. Also, it can be noted from figure 10(e) that the turbulent kinetic energy $k$ increases by decreasing $n$ . These trends have been reported before for the pipe geometry (Rudman et al. Reference Rudman, Blackburn, Graham and Pullum2004; Singh et al. Reference Singh, Rudman and Blackburn2017b ). Table 1 shows that flows with lower indices, for constant $Re_\tau$ , have lower generalised Reynolds numbers. For the Newtonian case, it has been observed that a decrease in the Reynolds number provokes a decrease in the maximum of the turbulent intensities and the turbulent shear stress (Moser et al. Reference Moser, Kim and Mansour1999). Despite the streamwise intensity in figure 10(a), other intensities and Reynolds shear stress behave the same here. It is important to note that the statistics in this figure are scaled using $u_\tau$ . As the shear-thinning effect increases, the friction velocity decreases. Consequently, although turbulent intensities decrease with decreasing $ n$ , this effect is not apparent under this scaling. However, if one were to scale, for example, $\overline {u'u'}$ using the bulk velocity, it would become evident that the shear-thinning effect reduces the turbulent intensity, leading to lower energy dissipation (see figure 11). It is noted in figure 10(e) that, as the fluid becomes more shear-thinning, the peaks in $k$ are also slightly shifted away from the wall. Overall, the same trends have been observed in the numerical (Dubief et al. Reference Dubief, White, Terrapon, Shaqfeh, Moin and Lele2004; Pereira et al. Reference Pereira, Mompean, Thais and Thompson2017) and experimental (Pinho & Whitelaw Reference Pinho and Whitelaw1990; Warholic, Massah & Hanratty Reference Warholic, Massah and Hanratty1999) viscoelastic context, where increasing the elasticity weakened the Reynolds shear stress and the stream-normal intensities, while the streamwise intensity was strengthened. In that sense, decreasing $n$ seems analogous to increasing the elasticity on viscoelastic fluids, and the net effect of the changes in the turbulent intensities heightened the peaks of $k$ , even though $Re_n$ is decreased. These effects have also been reported for pipe flows of PL fluids by Singh et al. (Reference Singh, Rudman and Blackburn2017b , Reference Singh, Rudman and Blackburn2018) and in the channel flow of Carreau fluids by Arosemena et al. (Reference Arosemena, Andersson and Solsvik2020). A connection between the effects of increasing elasticity and increasing the shear-thinning degree was discussed by Rudman et al. (Reference Rudman, Blackburn, Graham and Pullum2004), although the authors did not report the increase in the streamwise fluctuations, contradicting the observations of the present work and of Singh et al. (Reference Singh, Rudman and Blackburn2017b , Reference Singh, Rudman and Blackburn2018). Arosemena et al. (Reference Arosemena, Andersson and Solsvik2020) also pointed to the similarities to viscoelastic flows in the LDR regime.

Figure 10.

Averaged turbulence intensities $\overline {u'u'}$ , $\overline {v'v'}$ , $\overline {w'w'}$ , Reynolds shear stress $\overline {u'v'}$ , turbulent kinetic energy $k$ and viscous shear stress $\overline {\tau _{xy}}$ for the smooth channel cases with $Re_\tau = 395$ .

Figure 11.

Averaged streamwise turbulence intensity $\overline {u'u'}$ scaled by bulk velocity for the smooth channel cases with $Re_\tau = 395$ .

The Reynolds averaged momentum equation in streamwise direction is

(3.2) \begin{equation} \frac {d}{{\textrm d} y^+ } \overline {u'v'^+} - \frac {d}{{\textrm d} y^+ } \overline {\tau _{xy}^+} = \frac {1}{h^+}, \end{equation}

where the first term on the left represents the Reynolds shear force (the derivative of stress) and the second is the average viscous shear force. The superscript $+$ indicates scaling with proper wall quantities ( $\tau _w$ and $\delta _\nu$ ). By integrating this equation, one finds the relation between Reynolds and viscous stresses as

(3.3) \begin{equation} \overline {\tau _{xy}^+} - \overline {u'v'^+} = 1 - \frac {y^+}{h^+}. \end{equation}

It is common practice (Pinho Reference Pinho2003; Singh et al. Reference Singh, Rudman and Blackburn2017b ) to decompose viscosity into the average and the fluctuating part as $\mu = \overline {\mu } + \mu '$ and, therefore, it is possible to write the viscous stress as

(3.4) \begin{equation} \overline {\boldsymbol{\tau }} = 2 \overline {\mu } \overline {\unicode{x1D64E}} + 2 \overline{\mu^{\prime} {\unicode{x1D64E}}^{\prime}}. \end{equation}

The latter term in (3.4) vanishes in the Newtonian case. As observed by Anbarlooei et al. (Reference Anbarlooei, Cruz and Freire2013) and Singh et al. (Reference Singh, Rudman and Blackburn2017b ), it is non-negligible for GNF, especially in the near-wall region. It is important to notice that this decomposition does not reduce to wall quantities by approaching the boundaries ( $\overline {\mu } \neq \mu _w$ and $2 \overline {\mu } \overline {\unicode{x1D64E}} \neq \tau _w$ ) (Singh et al. Reference Singh, Rudman and Blackburn2017b ). Because of this, the length and velocity scales have been modified by Singh et al. (Reference Singh, Rudman and Blackburn2017b ) to overcome the related difficulties. However, it is more natural to work with stress in non-Newtonian fluids rather than the viscosity. Therefore, in the remainder of the present work, the decomposition $\tau _{ij} = \overline {\tau _{ij}} + {\tau ^{\prime}_{ij}}$ is preferred to the other forms and no further decomposition is applied to the stress tensor.

For the Newtonian case where $\overline {\tau _{xy}^+} = {\textrm d}U^+/{\textrm d}y^+$ , the viscous term rapidly approaches zero after the maximum of the Reynolds stress. However, as illustrated in figures 10(d) and 10( f), decreasing the PL index $n$ results in a ‘residual viscous stress’, which remains substantial until very close to the centre of the channel. In other words, the viscous stresses in shear-thinning fluids, in contrast to the Newtonian or shear-thickening fluids, have a much larger domain effect. The consequences of this effect on friction prediction are also discussed in § 3.2.

Figure 12.

Ratio of the viscous stress to the Reynolds stress gradients $F_V / F_R$ for the smooth channel simulations with $Re_\tau = 395$ .

It is interesting to study the effect of viscosity from the perspective of the viscous forces $F_V={\textrm d} \overline {\tau _{xy}^+} / {\textrm d}y^+$ and the Reynolds forces $F_R={\textrm d} \overline {u'v'^+} /{\textrm d}y^+$ . For this purpose, the four-layer structure of Wei et al. (Reference Wei, Fife, Klewicki and Mcmurtry2005) has been obtained for the current simulations. As discussed by Wei et al. (Reference Wei, Fife, Klewicki and Mcmurtry2005) for the Newtonian turbulent flow, one can divide the domain into: (i) the first layer next to the wall, where viscous force is dominant and balances the pressure force; (ii) the second layer where viscous and Reynolds forces are in balance; (iii) the third layer where the viscous force again dominates the flow; and (iv) the last layer in which the turbulence (Reynolds shear force) dominates the flow. Figure 12 shows the ratio of the viscous to the Reynolds forces ( $F_V/F_R$ ) for different PL indices. As evident, the first layer is the same as in the Newtonian case. The length of the second layer increases while $n$ decreases. In the third layer, the change in the position of the maximum Reynolds stress by changing $n$ is visible. The most noticeable and important change happens in the last layer, where it seems to be divided into two layers. In the first part, the turbulence dominates the flow as in the fourth layer of the Newtonian case. However, near the centre, viscous force becomes important again and cannot be neglected. To our knowledge, this behaviour has not been reported before. The end positions of each layer in the four-layer structure are listed in table 2. Note that because the fourth layer’s end is at the centre of the channel for all cases, this information has been omitted from table 2.

Table 2.

End of each layer in the four-layer structure of Wei et al. (Reference Wei, Fife, Klewicki and Mcmurtry2005). The end of the fourth layer is at $y/h = 1.0$ for all cases.

Figure 13 shows the average viscosity and shear rate for the current simulations. As expected, in figure 13(a), where the average viscosity is scaled with the viscosity at the centre of the channel ( $\overline {\mu }_c$ ), the shear-thinning fluid viscosity decreases as the wall is approached while the opposite is true for the shear-thickening fluid. Figure 13(c) shows the average shear rate scaled by the wall shear rate. As can be noticed, the shear rate drops faster near the centre for the more shear-thinning case. Also, it seems that it approaches the wall shear rate near the wall. However, close examination shows that by decreasing $n$ , there will be a growing difference between the average shear rate and the wall value ( ${\approx}3 \,\%$ at maximum). The same difference has been reported by Rudman et al. (Reference Rudman, Blackburn, Graham and Pullum2004), where authors propose a physical explanation for that. However, we believe this difference is a numerical artefact and should be eliminated using finer grids. Figure 13(b) compares the average viscosity and the viscosity obtained using the average shear rate and (2.3). There is once again a match near the wall with a growing difference towards the centre of the channel in the shear-thinning cases. Finally, in figure 13(d), the shear component $\overline {S_{xy}}$ of the mean strain rate tensor is compared with the strain rate. As evident, while $\overline {S_{xy}}$ vanishes at the centre, the mean shear rate never drops to zero. This is a very important observation since the absence of shear at the centre of the channel/pipe is commonly expected.

Figure 13.

Average viscosity and strain rate ratios of the DNS with $Re_\tau = 395$ . (a) Ratio between average viscosity $\overline {\mu }$ and the average channel centreline viscosity $\overline {\mu }_c$ . (b) Ratio between the average viscosity and the average shear rate viscosity $K\overline {\dot {\gamma }}^{n-1}$ . (c) Ratio between the average shear rate and the wall shear rate $\dot {\gamma }_w$ . (d) Ratio between the shear component of the mean strain rate tensor $\overline {S_{xy}}$ normalised by the average shear rate $\dot {\overline {\gamma }}$ .

To finalise this section, the effect of an increment in the $Re_\tau$ on the turbulent flow of PL fluid is investigated in figure 14. In this figure, the simulation for the case $n=0.5$ has been repeated for $Re_\tau =590$ . As it can be seen in figure 14(a), both velocity profiles collapse pretty well in the viscous, buffer and log regions. This result and the results presented by Rudman et al. (Reference Rudman, Graham, Blackburn and Pullum2002, Reference Rudman, Blackburn, Graham and Pullum2004), Singh, Rudman & Blackburn (Reference Singh, Rudman and Blackburn2018) and Zheng et al. (Reference Zheng, Rudman, Singh and Kuang2019) clearly show that there is a common log (or power) region for a fixed power-law index (as expected). This contradicts the results obtained using the LES method developed by Gnambode et al. (Reference Gnambode, Orlandi, Ould-Rouiss and Nicolas2015), where the increase of the Reynolds number broke such overlap. This indicates this LES method, although modified to consider the non-Newtonian effects, is not adequate, missing to properly consider some effects. However, in the LES approach of Inagaki (Reference Inagaki2021), increasing $Re$ caused the velocity profile to better adhere to the log law, the same as observed by Singh et al. (Reference Singh, Rudman and Blackburn2018) and Zheng et al. (Reference Zheng, Rudman, Singh and Kuang2019).

Turbulent intensities in figures 14(b), 14(c) and 14(d) behave the same as in the Newtonian case when the Reynolds number is increased. The Reynolds stress in figure 14(e) shows that the effect of viscosity in the middle section (residual viscous effects) decreases for the higher value of the Reynolds. This is also clearly visible in figure 14(d) for the average viscous stress. We expect that the presence of this residual viscous effect will be insignificant in higher Reynolds numbers.

Figure 14.

Velocity $U^+$ , turbulence intensities $\overline {u'u'}$ , $\overline {v'v'}$ and $\overline {w'w'}$ , Reynolds shear stress $\overline {u'v'}$ and viscous shear stress $\overline {\tau _{xy}}$ profiles for the power-law smooth channel case for $Re_\tau = 395$ and $Re_\tau = 590$ with $n=0.5$ .

3.2. Friction

As stated in § 1, our goal is to verify and develop a phenomenological model for power-law fluids. In this section, the friction factor computed from DNS is compared to the predictions of (1.2). For the PL fluid, the friction equation could be obtained by combining (1.2) with Kolmogorov’s velocity scale of (2.12), which results in

(3.5) \begin{equation} f = 0.059 \, Re_n^{-1/2(n+1)}. \end{equation}

The top part of table 3 contains the comparison of the predictions of (3.5) and the results of our DNS. As demonstrated, the theory developed by Anbarlooei et al. (Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018) predicts the friction factor of PL fluids with very good accuracy. However, we notice an increment in the error of (3.5) when the power-law index $n$ is decreased. Several reasons could be the cause. Based on the experimental evidence for the pipe flow and in the range of available data, the proportionality constant in (1.2) seems to be rheology independent (Anbarlooei et al. Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018). However, this may not be the case for very low values of $n$ . This discrepancy may be attributed to the residual viscous stress observed in figure 10( f).

Table 3.

Comparison of the DNS friction factor, $f_{DNS}$ , for turbulent channel flow of power-law fluids with predictions of (3.5), $f_p$ . Both the results of the current work and the reference results by Singh et al. (Reference Singh, Rudman and Blackburn2017b , Reference Singh, Rudman and Blackburn2018) have been analysed.

However, the discrepancy in the shear-thickening case could have an alternative explanation. Recently, Anbarlooei, Ramos & Cruz (Reference Anbarlooei, Ramos and Cruz2022) demonstrated that at sufficiently high Reynolds numbers, the mechanism of momentum transfer from the core flow to the viscous layer changes. In this extreme Reynolds regime, attached eddies, rather than Kolmogorov-sized eddies in the Blasius region, are responsible for the momentum transfer. This results in the following friction equation for Newtonian turbulent flows:

(3.6) \begin{equation} f = C_k Re_\tau ^{-1/6}, \end{equation}

where $C_k$ is the proportionally constant. The transition from the Blasius region to the extreme Reynolds regime begins around $Re = 10^5$ in the Newtonian case. Our preliminary investigation suggests that as $n$ increases, this transition occurs at lower Reynolds numbers. Given that the Kolmogorov scale is much larger for $n=2$ than in the Newtonian case, the attached eddies more quickly assume the role of Kolmogorov eddies, becoming the primary contributors to the momentum transfer process to the wall. If this explanation holds and considering the validity of (3.6) for the extreme shear thickening case, the error in friction factor prediction reduces to less than $5\,\%$ . We believe that a DNS study of shear-thickening cases ( $n=2$ ) could offer valuable insights into understanding this phase transition in the momentum transfer mechanism and enhance our understanding of extreme Reynolds number Newtonian turbulent flows.

As evident from figure 10( f), by decreasing $n$ , the average viscous stress in the vicinity of the centre of the channel starts to increase, which contradicts one of the main assumptions of the phenomenological models commonly employed (Gioia & Chakraborty Reference Gioia and Chakraborty2006; Gioia et al. Reference Gioia, Guttenberg, Goldenfeld and Chakraborty2010; Anbarlooei et al. Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018). However, it is expected that the viscous effects decrease by increasing the Reynolds number and, as a result, the viscous stress should also diminish. Therefore, the predictions of a more robust model should converge to the observations (DNS or experiments) for higher Reynolds numbers. Considering that most experimental data available have higher Reynolds numbers than current DNS, comparisons with Anbarlooei et al. (Reference Anbarlooei, Cruz, Ramos and Freire2015b , Reference Anbarlooei, Cruz, Ramos, Santos and Freire2017, Reference Anbarlooei, Cruz, Ramos, Santos and Freire2018) show very good agreement between (3.5) and experiments. Considering the slight change in the Reynolds number and the associated decrease in error, the result for $Re_\tau = 590$ partially supports this conjecture.

The same analysis for the DNS results of Singh et al. (Reference Singh, Rudman and Blackburn2017b , Reference Singh, Rudman and Blackburn2018) for the pipe geometry is shown in the lower part of table 3, which demonstrates the accuracy of our friction formula. Here, $0.079$ is assumed as the proportionality constant. When the residual viscous stress is present, especially for the cases where $n=0.6$ , the same trend discussed before is again visible.

It is also interesting to notice the behaviour of $Re_{\eta _w}$ in table 3 as an indicator of the bulk changes in the flow. Based on the results of the previous section, it was observed that the state of the turbulence was dramatically changed when $n$ decreased. Therefore, one expects to see the reflection of such changes in the Reynolds number defined based on the bulk properties. However, while $Re_n$ clearly shows a significant change, $Re_{\eta _w}$ remains almost indifferent. Based on this result, we believe that $Re_{\eta _w}$ is not an adequate quantity for turbulent flows of GNF.

3.3. Mean velocity profiles

In this section, the mean velocity profile of the PL simulations is compared with the logarithmic laws available in the literature and also with a new power-law profile developed here.

The MVP of PL fluids has been investigated before. Dodge & Metzner (Reference Dodge and Metzner1959) used measurements of frictional pressure-loss in pipes to propose a mean velocity profile of the following form:

(3.7) \begin{align} U^{+} &= \frac {5.66}{n^{0.75}} \log (y^+)-\frac {5.66}{n^{1.2}} \nonumber\\[6pt] &\quad +\,\frac {3.475}{n^{0.75}}\left [ 1.960 + 0.815\,n-1.628\,n\log \left (3+\frac {1}{n} \right ) \right ]. \end{align}

The slope of the log-term was assumed to vary with $n^{-0.75}$ , while further work by Clapp (Reference Clapp1961) showed no dependence on $n$ as

(3.8) \begin{equation} U^+ = A \, \ln (y^+) + B = 2.8 \, \ln (y^+) + \frac {3.8}{n}. \end{equation}

This equation was based on experiments conducted with fluids with PL indices between $0.7$ and $0.8$ and for Reynolds numbers up to $5 \times 10^4$ . Anbarlooei et al. (Reference Anbarlooei, Cruz and Freire2013) further adjusted (3.8) after comparisons to DNS and also to ensure it recovered the Newtonian log-law when $n=1$ , resulting in $A=2.5$ and $B=5.0 / n$ . Finally, using available DNS and a wide range of experimental data, Anbarlooei et al. (Reference Anbarlooei, Cruz, Ramos and Freire2015a ) refined (3.8) leading to

(3.9) \begin{equation} U^+ = 2.5 \ln (y^+) + (5.0 - 5.44 \ln (n) ). \end{equation}

Figure 15.

Mean velocity profiles of shear-thinning power-law fluids compared with the logarithmic laws of (3.9) and (3.8) on the left and the power-law velocity profile of (3.14) on the right.

Figure 15(a,c,e) shows the comparison of the obtained DNS mean velocity profiles with (3.9) and (3.8) with modified coefficients. The former equation represents the log region more accurately than Clapp’s equation, especially for the lower $n$ . The same conclusion has been reported by Zheng et al. (Reference Zheng, Rudman, Singh and Kuang2019).

One can obtain an expression for the mean velocity profile from the friction factor, by assuming that for sufficiently large $Re$ , the following incomplete similarity relation holds: $ U(y)=u_\tau \tilde {\Phi } (y/\delta _{\nu })$ for any $\delta _{\nu } \ll y \lt \delta$ . This is in contrast to Prandtl’s friction equation which is obtained by integrating the logarithmic velocity profile. Indeed, assuming the above mentioned incomplete similarity, one can find such an approximation for most of the flow’s domain, as

(3.10) \begin{equation} \frac {V}{u_\tau } = \left ( \frac {f}{2} \right )^{-1/2} = \frac {1}{h} \int _0^h \tilde {\Phi }(y/\delta _{\nu }) {\textrm d}y = \frac {1}{Re_\tau } \int _0^{Re_\tau } \tilde {\Phi }(y^+) {\textrm d}y^+. \end{equation}

Equation (3.10) can be rewritten as

(3.11) \begin{equation} \frac {Re_{\tau }}{\sqrt {C_{\tau }Re_{\tau }^{-\frac {2n}{5n+2}}}} = \int _{0}^{Re_{\tau }} U^+ (y^+) {\textrm d}y^+, \end{equation}

where here the proportionality constant $C_{\tau }$ is different from the constant in (3.5). The velocity profile can be obtained from (3.11) by changing $Re_{\tau }$ for $y^+$ and applying the derivative $ (\partial /\partial y^+ )$ on both sides to obtain

(3.12) \begin{equation} \frac {\partial }{\partial y^+} \frac {y^+}{\sqrt {f(y^+)}} = U^+ (y^+). \end{equation}

Finally, from the (3.12) and calculating the derivative, we have

(3.13) \begin{equation} U^+ = \frac {2(1+3n)(y^+)^{\frac {n}{2+5n}}}{(2+5n)\sqrt {C_{\tau }(n)}}. \end{equation}

This implies that the MVP for PL fluids satisfies

(3.14) \begin{equation} U^+=\frac {u}{u_\tau }\left (y\right ) = \tilde {\Phi }\Big(\frac {y}{\delta _{\nu }}\Big) = \beta {y^+}^{\frac {n}{5n+2}}. \end{equation}

If $n=1$ , we recover Blasius’ power law $U^+ (y^+)=\beta _{n=1} \,{y^+}^{(1/7)}$ . We remark that this power-law approximation can be entirely compatible with the log-law since they can be valid in different regions of the flow’s domain, and they can also match in a large region of the flow. One can easily find the relation between $C_f$ and $C_\tau$ as

(3.15) \begin{equation}\begin{aligned} C_{\tau ,{channel}} &= ( C_{f,{channel}} )^{\frac {4(n+1)}{5n+2}} 2^{\frac {n-2}{5n+2}}, \\[2pt] C_{\tau ,{pipe}} &=( C_{f,{pipe}} )^{\frac {4(n+1)}{5n+2}} 2^{-\frac {n+2}{5n+2}}.\end{aligned} \end{equation}

Having $C_f$ in (3.5), it is possible to use (3.10) and (3.15) to find $\beta$ as

(3.16) \begin{align} \begin{aligned} \beta _{{channel}} &= \left ( \frac { C_{\tau ,{channel}}}{2}\right )^{-\frac {1}{2}} \left ( 1 + \frac {n}{5n+2} \right ), \\[2pt] \beta _{{pipe}} &= \frac {1}{2} \left ( \frac { C_{\tau ,{pipe}}}{2}\right )^{-\frac {1}{2}} \left ( 1 + \frac {n}{5n+2} \right ) \left ( 2 + \frac {n}{5n+2} \right ). \end{aligned} \end{align}

Using (3.16) for $n=1$ , one finds $\beta _{{channel}}=8.5581$ and $\beta _{{pipe}}=8.5684$ . In addition to the surprisingly universal behaviour ( $c_{{channel}}\approx c_{{pipe}}$ ), these coefficients are in very good agreement with the empirical values provided by McKeon et al. (Reference McKeon, Li, Jiang, Morrison and Smits2004).

Figure 15(b,d,f) represents both the DNS and the power-law profiles obtained in this section. For each PL index, all three profiles of (3.14) for $n=1.0$ , $0.75$ and $0.5$ are plotted. As observed, the derived power-law velocity profile works very well for $n=1.0$ and $n=0.75$ , while it is close to the $n=0.5$ curve.

It must be emphasised here that although there is no fitted constant in this equation, the result still agrees well with the simulations. Also, it must be considered that (3.14) is obtained with the assumption of infinite $Re_\tau$ , while all of our DNS studies were done for quite low Reynolds numbers. In figure 16, the effect of increasing Reynolds number on the velocity profile is studied. Figure 16(b) shows the Newtonian DNS for different Reynolds numbers from Lee & Moser (Reference Lee and Moser2015). In this figure, the blue dashed lines represent (3.14) for the Newtonian case (power of $1/7$ ) and the green dash-dotted lines represent the best power-law fit for the lowest $Re_\tau$ . It is clear that by increasing $Re_\tau$ , the MVPs start to collapse on the power-law $1/7$ profile, while the lowest one does not show such a property. It is well known that for high enough Reynolds numbers, the Newtonian MVP will follow ${y^+}^{(1/7)}$ for a large portion of the channel. The same trend can be seen for the PL cases in figure 16(a). As evident, increasing $Re_\tau$ causes the higher Reynolds MVP to deviate from the best fit for the lowest into the so-called asymptotic curve (for $n=0.5$ , asymptotic will be $y^{1/9}$ ). So, it is logical to expect a better agreement at higher $Re_\tau$ .

Figure 16.

Comparison between the asymptotic power-law mean velocity profile of (3.14) and best fit power-law profile: (a) $n=0.5$ simulation; (b) Newtonian simulations of Lee & Moser (Reference Lee and Moser2015).

Given the lack of experimental data, studies and geometry mismatch in non-Newtonian turbulent flows, it is not easy to perform a direct comparison. However, the data published by Pinho & Whitelaw (Reference Pinho and Whitelaw1990) could be partially used for comparison with the results of our current work in smooth channel conditions. This reference investigates the turbulent flow of CMC solutions at varying concentrations in water, specifically within a pipe geometry. For low CMC concentrations, we can reasonably assume that viscoelastic effects are negligible. For instance, for a $0.2\,\%$ CMC solution, the rheology can be approximated using a power-law model with ( $n=0.75, K = 0.044$ ). Notably, one of the fully developed turbulent Reynolds numbers reported in the aforementioned study corresponds $Re_\tau\,{\approx}\,400$ .

Figure 17 presents a comparison of the experimental mean velocity profile with the power-law MVP developed in our study, (3.14). As shown, our power-law MVP captures a significant portion of the flow domain effectively. In this figure, the value of $\beta$ is determined using (3.16) without any fine-tuning to align the MVP with experimental data.

Figure 17.

Comparison of experimental and analytical mean velocity profiles for turbulent flow of a non-Newtonian power-law-type fluid at $Re_\tau \approx 400$ . Experimental mean velocity profile from Pinho & Whitelaw (Reference Pinho and Whitelaw1990) compared with the analytical profile from (3.14), with no fine-tuning of proportionality constants.

4.1. Verification of the simulation strategy

To verify the simulation strategy, our Newtonian results were compared with the DNS results of Ashrafian et al. (Reference Ashrafian, Andersson and Manhart2004). Figure 18 demonstrates the mean streamwise velocity profile obtained in the present work along with the reference profile. To compare, the velocity profile of the turbulent smooth channel at the same $Re_\tau$ is also shown in figure 18. The rough velocity profile is shifted downward. This is an expected behaviour since von Karman’s constant is not affected by the wall roughness, as reported by White (Reference White1991). As observed in figure 18, the present work’s simulation agrees well with the reference. The velocity shift, or roughness function, $\Delta U^+$ , is equal to 7.2, in agreement with the value of 7.1 reported by Ashrafian et al. (Reference Ashrafian, Andersson and Manhart2004).

Figure 18.

Streamwise mean velocity profiles for the Newtonian simulations of the smooth and rough channel of the present work, and reference rough channel by Ashrafian et al. (Reference Ashrafian, Andersson and Manhart2004).

The streamlines associated with the average velocities are demonstrated in figure 19, where a recirculation region between the roughness elements can be observed. Downstream of the rough element, the height of the recirculation is comparable to the rough element. Following the flow direction, we see a neck in the recirculation bubble after which the recirculation height increases again, upstream to the next rough element, but achieving a lower height. Comparing the Newtonian and the $n=0.5$ cases of figure 19, we can notice that the rheology timidly affects the recirculation area sizes and configuration. No separation is observed, while near the wall, there is counter flow, and directed flow occurs immediately above the rough element. In figure 19, four streamwise locations in the channel are highlighted. Location I is at the maximum stream function value, Location II is at the neck of the recirculation, Location III is at a point upstream of the rough element where the recirculation height is more pronounced and Location IV is at the roughness crest. The profiles of the turbulent stresses will be demonstrated along these four locations.

Figure 19.

Recirculations inside the cavities of the rough channel for the Newtonian and $n=0.5$ simulations, at $Re_\tau = 395$ . Locations I, II, III and IV are downstream, neck, upstream and at the roughness crest, respectively.

Figure 20.

Profiles of the Reynolds shear stress $\overline {u'v'}$ and of the streamwise $\overline {u'u'}$ , normal $\overline {v'v'}$ and spanwise $\overline {w'w'}$ turbulence intensities in the Newtonian rough channel simulation with $Re_\tau =395$ . Reference DNS results of Ashrafian et al. (Reference Ashrafian, Andersson and Manhart2004) and smooth channel results are displayed for comparison.

The averaged turbulence intensities and the Reynolds shear stress at location I of figure 19 are demonstrated in figure 20. The profiles of the reference DNS results of Ashrafian et al. (Reference Ashrafian, Andersson and Manhart2004) are also presented. The match between the MVPs of the present and reference works is also observed in the turbulent stresses.

4.2. Insights from the phenomenological model

This study began by reviewing the phenomenological model, which we believe explains the behaviour of the friction factor in both smooth and rough cases. As outlined in the introduction, after the transition to turbulence, the friction-Reynolds-number domain can be categorised into three main regions: (i) the Blasius region; (ii) the extreme Reynolds number region; and (iii) the fully rough regime. Naturally, there are transition intervals between these regions.

As discussed in the previous section, the cases corresponding to $n=1.0$ , $0.75$ and $0.5$ in the smooth channel fall within the Blasius region, while the $n=2.0$ case in the smooth channel lies in the extreme Reynolds number range. Both of these regions pertain to smooth-wall turbulence.

According to the phenomenological model, as the Reynolds number increases, the size of the dominant eddies responsible for momentum transfer from the core flow to the wall decreases. Eventually, these eddies scale with the roughness element size (or the spacing between them). Beyond this point, further increasing the Reynolds number no longer affects the dominant eddy size, leading to a plateau in the friction factor. These regions are illustrated in figure 21(a), where, for the sake of demonstration, we include experimental data from Nikuradse’s seminal experiments (sand-roughened pipe). This serves to highlight the key transitions in friction behaviour across different flow regimes.

Figure 21(b) presents our results for $n = 1.0$ in both smooth and rough channels ( $Re_\tau\,{=}\,395$ ). Two sets of experimental data are included in this figure. The first set, represented by triangles, corresponds to experimental data obtained in a smooth channel for a Newtonian fluid by Yimprasert et al. (Reference Yimprasert, Kvick, Alfredsson and Matsubara2021). The second set, shown as an open blue circle, refers to experimental data from Krogstad et al. (Reference Krogstad, Andersson, Bakken and Ashrafian2005) in a rough channel with the same geometry studied here, but at higher Reynolds number, $Re_\tau = 600$ . Additionally, we have included the laminar solution given by

(4.1) \begin{equation} f = \frac {2 \left ( 2 + \frac {1}{n}\right )}{Re_n} \end{equation}

and the Blasius equation (3.5) for comparison.

Figure 21.

Friction factor versus Reynolds number. (a) Experimental data on (Newtonian, sand-roughened pipe) from Nikuradse. (b) Experimental data (smooth channel) from Yimprasert et al. (Reference Yimprasert, Kvick, Alfredsson and Matsubara2021), compared with DNS results and experimental data (riblet-rough channel) from Krogstad et al. (Reference Krogstad, Andersson, Bakken and Ashrafian2005) ( $n=1.0$ ). (c) DNS results for $n = 0.75$ and (d) DNS results for $n = 0.5$ . In all panels, the solid line represents the laminar friction profile, the dash-dotted line corresponds to the Blasius correlation from (3.5) and the dotted line indicates the fully rough friction factor plateau for the Newtonian case.

As shown, our DNS results in the smooth channel aligns with the Blasius curve, as expected. However, the DNS results in the rough channel deviates from the Blasius curve and approaches the fully rough regime (see table 4 for the values). In fact, the experimental friction factor from Krogstad et al. (Reference Krogstad, Andersson, Bakken and Ashrafian2005) closely matches our DNS results, suggesting that the flow in the rough DNS is either in or very close to the fully rough regime. Therefore, the friction factor from this experiment is assumed to represent the constant friction factor in the fully rough regime, as indicated by the dotted line in the figure.

Table 4.

Comparison of the friction factor for the rough channel, $f_{{DNS}}^{R}$ , as obtained from DNS, with the predictions for the smooth channel $f_{{Pred}}^{S}$ derived from (3.5). $\Delta f$ highlights the increase in friction attributed to surface roughness.

To the best of our knowledge, experimental data (at least for the channel geometry) are not available for the non-Newtonian case. Therefore, figures 21(c) and 21(d) present only our DNS results. These figures also include the constant friction factor (fully rough regime) from the Newtonian case for comparison. Additionally, figure 21(d) presents our DNS results for the smooth channel at $Re_\tau = 595$ .

Again, for both cases ( $n = 0.75$ and $n = 0.5$ ), the smooth channel results follow the Blasius equation (3.5), while the rough simulations exhibit a sharp increase in the friction factor. It is important to note that for both the smooth and rough channel cases, $Re_\tau$ is held constant. As a result, the roughness effect is reflected in a reduction in the flow rate within the channel. As shown in table 4, a roughness size equivalent to $3.4\,\%$ of the channel height causes the friction factor to more than double.

Returning to our phenomenological model, the key question is whether the flow in the last two cases has reached the fully rough regime. If it has, the size of the dominant eddies should scale with the gap between the rough elements. If not, the dominant eddies near the wall should be larger than this gap. To address this question, it is necessary to analyse the near-wall vortical structures.

4.3. Near-wall vortical structures

Figure 22 illustrates the vortical structures near the wall for both smooth and rough channels for different rheologies. To identify these structures, we use the method proposed by Jeong & Hussain (Reference Jeong and Hussain1995), which is based on the eigenvalues of the symmetric tensor $\unicode{x1D64E}^2+\Omega ^2$ . Here, $\unicode{x1D64E}$ is the rate-of-strain tensor and $\Omega$ is the vorticity tensor. According to this approach, the vortex core geometry is visualised as an iso-surface of $\lambda _2$ , the median eigenvalue of $\unicode{x1D64E}^2+\Omega ^2$ . The $\lambda _2$ vortex criterion states that a vortex exists in a region where $\lambda _2 \lt 0$ . In such a case, the flow is dominated by rotational effects rather than stretching or shearing.

Figure 22.

Vortical structures for $\lambda _2=-0.03$ at $y=1.5k$ plane. From top to bottom, $n=1.0, 0.75, 0.5$ ; and from left to right, smooth and rough channel. All cases for $Re_{\tau }=395$ . The iso-surfaces are shown coloured with dimensionless velocity $u/V$ . Flow direction from left to right.

Figure 22 presents the iso-surfaces of $\lambda _2 = -0.03$ up to $y = 1.5 k / h$ , i.e. slightly above the height of the rough element. These surfaces are colour-coded based on the streamwise velocity, normalised by the bulk velocity. Panels (a,c,e) correspond to the smooth channel case, while panels (b,d,f) represent the rough channel. From top to bottom, the cases correspond to $n = 1.0$ , $0.75$ and $0.5$ .

As shown in this figure, in all smooth channel cases, the majority of the vortical structures are aligned in the streamwise direction. In contrast, the presence of roughness elements significantly disrupts the near-wall structures. For example, when $n=1.0$ (figure 22 b), the structures in the rough case are predominantly smaller, as they are disrupted by the roughness elements. We observe that vortex structures form between the riblets, and their sizes scale with the spacing between them. To provide a clearer visualisation, a zoomed-in view of the rough channel structures is presented in figure 23. This figure highlights only the structures trapped between the rough elements. For $ n = 1.0$ (figure 23 a), the structures completely fill the gap and their size is clearly proportional to the spacing between the roughness elements. This, combined with our phenomenological model, supports our claim from the previous section that the $ n = 1.0$ rough case is in the fully rough regime. Returning to figures 22(a) and 22(b), the velocities of the structures near the wall decrease substantially in the rough case. This reduction in vortex velocity is reflected in a shift in the velocity profile and an overall increase in the friction factor in the rough channel.

A similar behaviour is observed for the shear-thinning cases (see figure 22). Again, the long, streamwise-aligned structures in the smooth channel are transformed into structures predominantly confined between the roughness elements. However, as the shear-thinning effect increases, these structures occupy a progressively smaller portion of the riblet spacing (see figures 23 b and 23 c). As seen in figure 23(b), for $ n = 0.75$ , the structures still occupy a significant portion of the gap. Interestingly, in figure 21(c), the friction factor of this shear-thinning case is close to the value obtained for the fully rough Newtonian case. This may suggest that friction in the fully rough regime, where dominant structures near the wall occupy all the gap, is independent of rheology. However, additional simulations are required to confirm or refute this hypothesis.

Figure 22(c) further illustrates how roughness disrupts the near-wall structures, fragmenting them into smaller vortices. Additionally, as shown in figure 23(c) for the $ n = 0.5$ case, the structures are significantly smaller than the gap between the roughness elements. According to our phenomenological model, smaller vortices correspond to lower friction factors. This trend is evident in figure 21(d), where the friction factor for the rough case with $ n = 0.5$ is lower than that of the fully rough Newtonian case. It is also noteworthy that as $n$ decreases, the structures retain a larger portion of their streamwise velocity. This can be attributed to the drag reduction resulting from the shear-thinning effect.

According to the phenomenological theory discussed in the introduction, as the Reynolds number increases, the dominant eddy size – responsible for momentum transfer towards the wall – eventually matches the size of the eddies trapped between the roughness elements. Beyond this point, further increasing the Reynolds number no longer affects the dominant eddy size, resulting in a constant friction factor (the fully rough friction plateau in figure 1 or figure 21 a). The results for the $ n = 0.5$ rough channel may suggest that shear-thinning fluids exhibit different behaviour at higher Reynolds numbers compared with their Newtonian counterparts. It is possible that organised roughness, such as that investigated in the present study, behaves differently at high Reynolds numbers than the sand-roughened tubes used in Nikuradse’s experiments. Experimental evidence for the latter case, as reported by Santos et al. (Reference Santos, Anbarlooei, Cruz and Ramos2018), suggests that after reaching a plateau in friction, a further increase in Reynolds number leads to a deviation from this plateau, with friction beginning to decrease as the Reynolds number continues to rise. To refine the phenomenological model introduced at the beginning of this work and accurately explain the observed behaviour, additional DNS studies at higher Reynolds numbers are required.

The present results suggest the existence of two distinct mechanisms contributing to drag enhancement induced by rough wall elements. The first is an increase in momentum transfer, driven by the generation of vortical structures at the upstream upper corner of each roughness element. The second mechanism is associated with form drag, which arises from the bluff-body-like resistance of the elements themselves.

Both mechanisms are influenced by the geometric arrangement of the roughness elements – specifically, their height, their length and the spacing between them. As shown in figure 22, the roughness-induced disruption of vortex structures enhances turbulence and leads to an increase in turbulent shear stress.

The form drag mechanism, in contrast, involves resistance similar to that experienced by flow past a bluff body. In this case, the intensity of the drag force is governed by the size of vortices trapped in the cavities between adjacent roughness elements, which originate near the downstream upper corner. As illustrated in figure 23, form drag is strongest for $n=1.0$ and decreases as the flow behaviour index $n$ is reduced.

Finally, figure 25 presents the spatial distribution of local shear stress. Together with figure 23, it supports the proposed correlation between the size of the trapped vortex structures and the local friction factor, reinforcing the dual-mechanism interpretation of roughness-induced drag enhancement.

Figure 23.

Vortical structures closer view between two rough elements. $\lambda _2=-0.03$ , $y=0.99k$ plane. From left to right, $n=1.0, 0.75, 0.5$ . All cases for $Re_{\tau }=395$ . The iso-surfaces are shown coloured with dimensionless velocity $u/V$ .

4.4. Turbulence statistics

The velocity profiles of the rough wall PL simulations are presented in figures 24(a) and 24(b) for the cases with $n=0.75$ and $n=0.5$ , respectively. The corresponding smooth wall profiles were also included in both figures to demonstrate the velocity shifts. The behaviour of the MVPs in both cases is similar to the Newtonian one; they were shifted with similar $\Delta U^+$ values of $7.6$ in figure 24(a) and $7.8$ in figure 24(b). This indicates that the power-law index of the velocity profiles is not significantly affected by roughness.

Figure 24.

Shear-thinning velocity profiles for the rough and smooth channels with $Re_\tau = 395$ . The logarithmic (3.9) and power-law (3.14) profiles are also displayed. (a) Power-law index of $n = 0.75$ . (b) Power-law index of $n = 0.5$ .

The theoretical logarithmic and power-law profiles of (3.9) and (3.14) are also presented in figure 24. In the $n=0.75$ simulation, the proposed equation reproduces well the numerical results, indicating that the power-law expression (3.14) appropriately describes the velocity profile at least as well as the log-law of (3.9). The results for $n=0.5$ show that neither (3.9) nor (3.14) can accurately predict the velocity profile, with the former performing slightly better. Also, as discussed in § 3.1, the power-law profile (3.14) adheres better to higher Reynolds numbers, while the rough simulations have a modest Reynolds of $Re_\tau = 395$ .

The protruding shape of the roughness elements causes a combination of form drag and skin friction on the channel wall. The local friction profile $f(x) = 2 \tau _w(x)/\rho V^2$ along the bottom of the channel and the crest of the rods is presented in figure 25 for the shear-thinning and Newtonian cases, where $\lambda$ is the distance from the beginning of one gap until the beginning of the subsequent gap (which also coincides with the end of the rough element). In the bottom of the cavities, the recirculations observed in figure 19 induce slightly negative local friction. In this region, the $f(x)$ profiles were considerably close for the three values of the power-law index considered. The saddle point of location II causes a local peak at $x/\lambda \approx 0.65$ , where $f(x)$ becomes less negative. However, on the top of the crests, the friction profiles were dependent on the PL index $n$ . As expected, a considerable drag reduction is observed as the fluid becomes more shear-thinning.

Figure 25.

Local friction $f(x)$ profiles for the shear-thinning and Newtonian rough simulations with $Re_\tau = 395$ . $\lambda$ is the distance between the rods.

The Reynolds shear stress and the turbulence intensities profiles are shown in figure 26 at the four locations of figure 19 for the simulation with $n=0.75$ . Except for the streamwise intensity $\overline {u'u'}$ , the roughness effects are only present at up to approximately $10\,\%$ of the channel’s half-height $h$ . Figure 26(b) demonstrates that for $\overline {u'u'}$ , the roughness affects a longer segment, causing the largest deviation from the smooth case. At up to approximately $20\,\%$ of $h$ , $\overline {u'u'}$ is considerably below the smooth case along all of the four locations. The wall-normal $\overline {v'v'}$ and spanwise $\overline {w'w'}$ intensities presented peaks at locations II and III. At locations I and IV, a decrease in the magnitude of the shear stress $\overline {u'v'}$ was observed, while an increase occurred at location III. Out of the rough layer in which roughness considerably alters the profiles ( $y/h \lessapprox 0.2$ ), a close match with the smooth profiles was observed.

Figure 26.

Profiles of the Reynolds shear stress $\overline {u'v'}$ and of the streamwise $\overline {u'u'}$ , normal $\overline {v'v'}$ and spanwise $\overline {w'w'}$ turbulence intensities in the $n = 0.75$ rough channel simulation with $Re_\tau =395$ .

The results of the $n=0.75$ case did not deviate much from the Newtonian scenario. Ashrafian et al. (Reference Ashrafian, Andersson and Manhart2004) reported a discrepancy between smooth and rough turbulent stresses that covered approximately $17\,\%$ of $h$ , with the largest differences also observed on $\overline {u'u'}$ . However, the difference between roughened and smooth channels was much more pronounced in the more shear-thinning case, with $n=0.5$ . Figure 27 reveals considerable discrepancies between smooth and rough turbulent stresses along the channel height. In contrast to what was observed for the $n=0.75$ case in figure 26, the effects of roughness were not restricted to a rough sub-layer of $y/h \leqslant 0.2$ and the whole profile was changed. More discrete differences were observed on $\overline {u'v'}$ (figure 27 a), where the same trends along locations I–IV were observed and the remainder of the profile was slightly shifted from the smooth one. Again, the highest deviation from the smooth case was observed on the streamwise intensity $\overline {u'u'}$ . Figure 27(b) demonstrates that in the rough sub-layer, the smooth profile achieved much higher values, while in the remainder of the channel, the rough case was considerably above. The overall behaviours at locations I–IV for $\overline {v'v'}$ and $\overline {w'w'}$ in figures 27(c) and 27(d) were the same as observed for $n=0.75$ but with attenuated magnitudes. It is observed in figure 27(c) that $\overline {v'v'}$ behaves similarly to $\overline {u'u'}$ , but with more discrete differences between rough and smooth profiles and shifted peak positions. In the spanwise intensity $\overline {w'w'}$ (figure 27 d), a higher peak on the rough channel is observed, followed by a large segment where the smooth profile is above the rough. Close to the centre of the channel, $\overline {v'v'}$ and $\overline {w'w'}$ are slightly below the smooth profile.

Figure 27.

Profiles of the Reynolds shear stress $\overline {u'v'}$ , and of the streamwise $\overline {u'u'}$ , normal $\overline {v'v'}$ and spanwise $\overline {w'w'}$ turbulence intensities in the $n = 0.5$ rough channel simulation with $Re_\tau =395$ .