Turbulent channel flow of generalized Newtonian fluids at a low Reynolds number

Abstract

Several studies concerning the turbulent pipe flow of generalized Newtonian (GN) fluids may be found in the literature, but not for channel flow, although that has been extensively studied for other types of non-Newtonian fluids, such as those with viscoelastic effects. Direct numerical simulations corresponding to statistically converged turbulent channel flow of GN fluids at a low frictional Reynolds number have been performed. The shear-dependent viscosity is introduced through the Carreau fluid model, and results corresponding to the Newtonian fluid case are compared to those of moderate shear-thickening and shear-thinning fluid behaviour. The different statistics studied reveal that shear-dependent fluid rheology appears mainly to affect the flow within the inner layer region and with shear-thinning behaviour; suppressing near-wall structures such as quasi-streamwise vortices and low-speed streaks, inhibiting turbulence generating events and leading to different drag reduction features. These include: enhancement of streamwise turbulence intensity and suppression of the other cross-sectional intensities, decrease of the Reynolds shear stress (leading to a lessening in turbulent production), decrease in energy redistribution between individual components of the Reynolds stress tensor through the velocity–pressure gradient term and overall increase in turbulence anisotropy at both small and large scales. In particular, it is noted that at the channel centre ‘rod-like’ turbulence states, a known low-Reynolds-number behaviour, are more clearly seen with shear-thinning fluid rheology.

1. Introduction

Shear-dominated wall-bounded turbulent flows such as in pipes, channels and boundary layers are of utmost importance. Approximately 25 % of the energy used in industry is destined to transport fluids through pipes and channels, or to propel vehicles in air or water, and approximately a quarter of that energy is irreversibly dissipated near walls (Jiménez 2013). Furthermore, in many industrial settings, such as bioreactors in biochemical plants or the drilling machines used in petroleum extraction, the transported fluid is non-Newtonian. A non-Newtonian fluid has non-uniform viscosity, which may depend on shear stress history and/or strain rate in a nonlinear manner, and it is typically classified within three main groups (see e.g. Irgens 2014): (i) time-independent fluids, in which fluid viscosity is not a function of time, (ii) time-dependent fluids and (iii) viscoelastic fluids consisting of materials with partial elastic recovery but also with viscous features.

The interest in non-Newtonian flows has increased since Toms (1948) reported frictional drag reductions in turbulent pipe flows due to a small amount of polymer additives. The first explanation of the drag reduction phenomenon may be attributed to Lumley (1969, 1973). He reasoned that the expansion of molecules, mainly within the buffer layer, leads to an increase in effective viscosity, the dampening of small eddies, the reduction of the Reynolds shear stress, the thickening of the sublayer and the consequent drag reduction. Results reported by Achia & Thompson (1977) also supported the idea about stretched molecules leading to less ‘bursting’ (Kline et al. 1967) and thus less turbulent kinetic energy production and subsequent drag reduction. The other explanation for frictional drag reduction, put forward by Tabor & de Gennes (1986), attributes the decrease in drag to the elasticity of the polymer additives and argues that the increase in effective viscosity is rather small and therefore negligible.

Many studies (e.g. Den Toonder et al. 1997; De Angelis, Casciola & Piva 2002; Min et al. 2003b; Ptasinski et al. 2003; Escudier, Nickson & Poole 2009; Shahmardi et al. 2019) have been committed to understand variations of drag in turbulent pipes, ducts and channels of polymeric solutions where viscoelastic effects are important. However, for a wide range of materials, the non-Newtonian rheology is mostly strain-rate-dependent and viscoelastic effects may be neglected (Rudman et al. 2004). Generalized Newtonian (GN) fluids are a class of time-independent, purely viscous, non-Newtonian fluids commonly encountered in numerous engineering and commercial applications, e.g. fluids in bioreactors displaying shear-thinning behaviour, drilling fluids, cosmetics or food products. In a GN fluid, the stress tensor due to viscous effects, ${\tau _{ij,vis}}$, is given by

(1.1)\begin{equation} {\tau_{ij,vis}} =2 \mu {\mathsf{S}}_{ij}, \end{equation}

where $\mu = \mu (\dot {\gamma })$ is the apparent dynamic viscosity, which solely depends on the strain rate $\dot {\gamma }=(2{\mathsf {S}}_{ij}{\mathsf {S}}_{ji})^{1/2}$, and ${\mathsf {S}}_{ij}$ is the strain-rate tensor. Note as well that certain materials flow like a GN fluid once a certain yield-stress value is exceeded. Such materials are called time-independent yield-stress or viscoplastic fluids.

Experimental studies about turbulent flows of GN fluids have primarily focused on friction factor measurements and one-point statistics. Metzner & Reed (1955) and Dodge & Metzner (1959) proposed correlations for the Fanning friction factor based on an alternative definition for the Reynolds number. Park et al. (1989) reported an increase in the mean velocity and axial turbulence intensities and a decrease in the tangential intensities, in pipe flows, caused by shear-thinning behaviour in GN fluids. Pinho & Whitelaw (1990) and Pereira & Pinho (1994) additionally reported a suppression of the radial turbulence intensities and a delay in transition from laminar to turbulent flow due to shear thinning. Rudman & Blackburn (2003) presented similar findings and also compared turbulent flows of GN fluids with those of yield-stress fluids. In Rudman & Blackburn (2003) it is recognized that even a small amount of elasticity may importantly affect the results.

Experimental issues such as the difficulty in removing unwanted plastic and viscoelastic effects from solutions with polymer additives and the necessity of specialized equipment to study non-Newtonian fluids, which are mostly not optically transparent (Gavrilov & Rudyak 2016a), have motivated a growing interest in numerical studies. Direct numerical simulation (DNS) is particularly appealing because it does not require closures for turbulence modelling and allows us to include/exclude plastic and elastic effects from the employed rheological model. Nonetheless, despite their industrial relevance and the advantages of direct simulations, few DNS studies about turbulent flows of GN fluids are available in the literature.

Rudman & Blackburn (2003), Rudman et al. (2004) and Rudman & Blackburn (2006) presented the first group of DNS studies for turbulent flows of GN and yield-stress fluids. In agreement with the experimental investigations for turbulent pipe flows, an increase in the mean axial velocity, delay in transition to turbulence, suppression of radial and tangential turbulence intensities and enhancement of the axial intensity with shear thinning were noted. Also, the decrease in the root-mean-square (r.m.s.) values of the axial vorticity fluctuations and reduced turbulence production with shear-thinning behaviour were documented. Gavrilov & Rudyak (2016a) showed similar results but in addition reported an increase in turbulent kinetic energy with increasing shear thinning and, motivated by the work of Escudier et al. (2009), studied for the first time large-scale anisotropy of a purely viscous GN fluid flow through anisotropy-invariant maps (Lumley & Newman 1977) of the Reynolds stress anisotropy tensor. Here the noted increase in anisotropy near the wall with shear thinning is attributed to a suppression of the mechanism of redistribution of fluctuation energy between individual components of the Reynolds stress tensor, but the corresponding budgets, showing such decrease in energy redistribution, are not actually presented. Also, small-scale anisotropy (see e.g. Antonia, Kim & Browne 1991; Yeung & Brasseur 1991), equally important for realizable turbulence, is not studied.

More recently, Singh, Rudman & Blackburn (2017b) considered the influence of increasing shear-thinning behaviour of GN fluids, in turbulent pipes flows, on first- and second-order statistics. In the same publication and likely motivated by the work of Ptasinski et al. (2003), the mean and turbulent kinetic energy budgets are presented for the first time. However, the individual budgets of the second moments of the velocity fluctuations are not shown. The budgets for the Reynolds stresses not only allow us to understand how the different terms contribute to the corresponding stresses, and to the overall turbulent kinetic energy, but also are necessary to directly evaluate closure models for turbulence. Singh et al. (2017b), based on the joint probability distributions of the axial and wall-normal velocity fluctuations at some wall-normal positions, suggested as well that shear-thinning rheology suppresses contributions from ‘sweeps’ and ‘ejections’ (Wallace, Eckelmann & Brodkey 1972) to the Reynolds shear stress. Nonetheless, in the publication, the cause for the variation with rheology in the contributions from those physical events is not explained.

Singh, Rudman & Blackburn (2017a), on the other hand, studied the effect of yield stress on a turbulent pipe flow of a GN fluid. Here, it is found that the effect of increasing the yield stress is similar to an increase in shear-thinning behaviour, with the important difference that the new stress, arising due to fluctuations in viscosity, increases as the pipe's core is approached. Subsequently, Singh, Rudman & Blackburn (2018) considered Reynolds-number effects on a turbulent pipe flow of a GN fluid. In this investigation, up to a moderate frictional Reynolds number, it is reported that rheological effects are still present. In the paper, it is also noted that the mean viscosity profile appears to be Reynolds-number-independent. Finally, with respect to recent DNS studies of turbulent flows of GN fluids, Zheng et al. (2019) compared finite-volume-scheme-based predictions obtained using OpenFOAM, for low-order statistics in a turbulent pipe flow, with high-order spectral-element DNS code results. The study reported that turbulence statistics predicted by OpenFOAM for shear-thinning fluids usually differ by less than 10 %.

In addition, direct simulations of GN fluid flow based on lattice Boltzmann methods have been performed as well, e.g. Gabbanelli, Drazer & Koplik (2005), Yoshino et al. (2007), Wang & Ho (2011) and more recently Chen & Shu (2020). However, those studies are generally limited to the laminar flow regime.

Other non-DNS numerical studies are based either on large-eddy simulations or on solving the time-averaged governing equations with some closure for the Reynolds stress tensor and the correlation term appearing with non-Newtonian rheology. Consider, for instance, Malin (1997), Cruz & Pinho (2003), Sungkorn, Derksen & Khinast (2012), Ohta & Miyashita (2014), Gnambode et al. (2015) and Gavrilov & Rudyak (2016b).

In this work, results from DNS of statistically converged turbulent channel flow corresponding to GN fluids, at a low Reynolds number, are presented and qualitatively compared to those of turbulent channel flow of viscoelastic fluids and another canonical flow of GN fluids such as turbulent pipe flow. In the results sections, aside from examining first- and second-order statistics of GN fluid flows, physical motions contributing to the turbulence production are considered in detail to comprehend how changes with rheology, in dominating fluctuations and their large intermittent values, cause variations in the shear stress budget and more specifically in the Reynolds shear stress. Also, all relevant non-zero Reynolds stress budgets for turbulent channel flow of GN fluids are presented for the first time, allowing us, for example, to better understand the decrease/increase in energy redistributed from streamwise fluctuations with shear thinning/shear thickening. Furthermore, not only large-scale but also small-scale anisotropy of turbulent GN fluid flow is appraised and, finally, at the end of the paper, the different reported drag-reducing features with shear-thinning behaviour are discussed in light of variations noted in the near-wall structures, i.e. quasi-streamwise vortices and velocity streaks.

2. Formulation

2.1. Governing equations and characteristic scales

Consider the equations, in a Cartesian coordinate system and index notation, governing mass and momentum conservation in the absence of external forces for a GN fluid, i.e.

(2.1)\begin{gather} \frac{\partial u_i}{\partial x_i} =0, \end{gather}

(2.2)\begin{gather} \frac{\partial u_i}{\partial t}+\frac{\partial (u_i u_j)}{\partial x_j} =\frac{1}{\rho}\frac{\partial}{\partial x_j}(-p\delta_{ij}+2\mu{\mathsf{S}}_{ij}), \end{gather}

where streamwise, wall-normal and spanwise directions are denoted by ${\boldsymbol {x}}=(x_1,x_2,x_3)=(x,y,z)$, the corresponding instantaneous velocity field by ${\boldsymbol {u}} =(u_1,u_2,u_3)=(u,v,w)$ and the pressure by $p$. Here $t$ denotes time, $\rho$ is the density of the incompressible, isothermal GN fluid, $\delta _{ij}$ is the Kronecker delta and ${\mathsf {S}}_{ij}=(\partial u_i/\partial x_j + \partial u_j/\partial x_i)/2$ is the aforementioned strain-rate tensor.

For GN fluids, shear-thinning and shear-thickening behaviours may be reproduced through different rheological models (constitutive equations to relate apparent viscosity and strain rate) such as the power-law (PL) or the Carreau fluid models. In both models, certain parameters are to be specified based on experimental data obtained from a rheogram. Regarding particularities of such models, for example, PL is simpler but leads to non-physical results at large and low shear-rate values whilst the Carreau fluid model may be considered a truncated power-law introduced to avoid this issue. The apparent viscosity for the Carreau fluid model (see e.g. Irgens 2014) is given by

(2.3)\begin{equation} \mu = \mu_{\infty} + (\mu_0 - \mu_{\infty}) [1+(\varLambda\dot{\gamma})^{2}]^{(\alpha-1)/2}, \end{equation}

where $\mu _{\infty }$ and $\mu _0$ are the ‘infinite’ and ‘zero’ shear-rate viscosities, respectively, $\varLambda$ is a time constant and $\alpha$ is the flow index, which for shear thinning is to be less than unity and for shear thickening more than unity.

Considering a characteristic velocity $U_c$, viscosity $\mu _c$, length $L_c$, time $L_c/U_c$ and stress $\rho U_c^{2}$ for the flow, (2.1)–(2.3) can be rewritten in non-dimensional form as

(2.4)\begin{gather} \frac{\partial u_i}{\partial x_i} =0, \end{gather}

(2.5)\begin{gather} \frac{\partial u_i}{\partial t}+\frac{\partial (u_i u_j)}{\partial x_j} =\frac{\partial}{\partial x_j} \left(-p\delta_{ij}+\frac{2\beta}{ Re}{\mathsf{S}}_{ij}\right), \end{gather}

(2.6)\begin{gather} \beta = \frac{\mu_{\infty} }{\mu_{_c}} + \left(\frac{\mu_{0}}{\mu_c}-\frac{\mu_{\infty}}{\mu_{c}}\right) [1+(\varLambda\dot{\gamma})^{2}]^{(\alpha-1)/2}, \end{gather}

where $Re=\rho U_cL_c/\mu _c$ is the Reynolds number and $\beta =\mu /\mu _c$ is the viscosity ratio between the apparent fluid viscosity and the characteristic viscosity. Observe that, for simplicity, the same notation as in (2.1)–(2.3) has been used in (2.4)–(2.6).

For wall-bounded shear flows, typically $U_c=u_{\tau }$ and $L_c=h$ in the outer layer and $L_c = \mu _c/(\rho u_{\tau })$ in the inner layer; see for instance Pope (2000). Here $u_{\tau }$ and $h$ refer to the wall friction velocity and the outer length scale, respectively. Regarding the characteristic viscosity, its selection is not clear and is open to debate within the scientific community (see e.g. Rudman et al. 2004). Here, as in Pinho & Whitelaw (1990), Ptasinski et al. (2003) and Singh et al. (2017a), the nominal wall viscosity $\mu _w$ is taken as characteristic viscosity, i.e. $\mu _c = \mu _w$. The complete set of governing equations for a GN fluid is given by (2.4) and (2.5) and a constitutive equation for the apparent viscosity such as (2.6). It is worth mentioning that $\alpha = 1$, in the aforementioned constitutive equation, allows us to recover Newtonian fluid behaviour.

2.1.1. Averaged governing equations

Introducing the Reynolds decomposition, i.e. splitting the variables into an ensemble average $\bar {(~)}$ and a fluctuating component ${(~)}^{\prime }$, as $u_i = \bar {u}_i + u_i^{\prime },\ p=\bar {p} + p^{\prime },\ \beta = \bar {\beta } + \beta ^{\prime }$ and ${\mathsf {S}}_{ij} = \bar {{\mathsf {S}}}_{ij}+ {{\mathsf {S}}_{ij}}^{\prime }$, into (2.4) and (2.5) and taking the average leads to

(2.7)\begin{gather} \frac{\partial \bar{u}_i}{\partial x_i} =0, \end{gather}

(2.8)\begin{gather} Re\left[\frac{\partial \bar{u}_i}{\partial t}+\frac{\partial (\bar{u}_i \bar{u}_j)}{\partial x_j} \right]=- Re\frac{\partial \bar{p}}{\partial x_i} +\frac{\partial}{\partial x_j}\left(2\bar{\beta}\bar{{\mathsf{S}}}_{ij} +2\overline{\beta^{\prime}{\mathsf{S}}_{ij}^{\prime}}-\overline{u_i^{\prime} u_j^{\prime}}\right), \end{gather}

where, in comparison with the Reynolds-averaged Navier–Stokes equations for a Newtonian fluid, a new non-Newtonian term $(2\overline {\beta ^{\prime }{\mathsf {S}}_{ij}^{\prime }})$ arises. This new viscous stress, denoted as ‘turbulent viscous stress’ in Singh et al. (2017b), due to fluctuations in viscosity, is analogous to the recognized polymer stress in Ptasinski et al. (2001, 2003). Therefore, the total mean shear stress $\bar {\tau }=\bar {\tau }_{12}$ for a GN fluid is given by

(2.9)\begin{equation} \bar{\tau} = 2\bar{\beta}\bar{{\mathsf{S}}}_{12} +2\overline{\beta^{\prime}{\mathsf{S}}_{12}^{\prime}}-\overline{u_1^{\prime} u_2^{\prime}}. \end{equation}

Since in a GN fluid the viscosity depends on the velocity gradient, through the strain rate, its fluctuating part is not expected to vanish at the wall $(\cdot |_{w})$. Thus, $\bar {\tau }|_{w} = 2(\bar {\beta } \bar {{\mathsf {S}}}_{12}+\overline {\beta ^{\prime }{\mathsf {S}}_{12}^{\prime }})|_w$.

2.1.2. Reynolds stress budget equations

The transport equation for the correlation of the velocity fluctuations, $\overline {u_i^{\prime }u_k^{\prime }}$, corresponding to a GN fluid is deduced in a similar manner as it is deduced for a Newtonian fluid. Thus, the added products $(u_i^{\prime }\mathrm {D}(u_k^{\prime })/\mathrm {D}t +u_k^{\prime }\mathrm {D}(u_i^{\prime })/\mathrm {D}t)$ are time-averaged. Here $\mathrm {D}(~)/\mathrm {D}t = \partial (~)/\partial t + \bar {u}_j \partial (~)/\partial x_j$ is the material time derivative.

The transport equation for the velocity correlation, see for instance Pinho (2003) or appendix A, in non-dimensional form is given by

(2.10)\begin{align} Re\frac{\mathrm D }{\mathrm D t} \left(\overline{u_i^{\prime}u_k^{\prime}}\right) &= Re\left[-\overline{u_i^{\prime}u_j^{\prime}} \frac{\partial \overline{u_k}}{\partial x_j} -\overline{u_k^{\prime}u_j^{\prime}} \frac{\partial \overline{u_i}}{\partial x_j} -\frac{\partial}{\partial x_j} \left(\overline{u_i^{\prime}u_j^{\prime}u_k^{\prime}}\right)\right] + {\bar{\beta}}\frac{\partial^{2}}{\partial x_j \partial x_j}\left(\overline{u_i^{\prime}u_k^{\prime}}\right) \nonumber\\ &\quad- Re\overline{\left(u_i^{\prime}\frac{\partial p^{\prime} }{\partial x_k}+u_k^{\prime}\frac{\partial p^{\prime} }{\partial x_i}\right)} - 2\bar{\beta} \overline{\frac{\partial u_i^{\prime}}{\partial x_j}\frac{\partial u_k^{\prime}}{\partial x_j}}\notag\\ &\quad+ \left(\frac{\partial \bar{\beta}}{\partial x_j}\right)\left[\frac{\partial}{\partial x_j} \left(\overline{u_i^{\prime}u_k^{\prime}} \right)+ \overline{u_i^{\prime}\frac{\partial u_j^{\prime}}{\partial x_k}}+\overline{u_k^{\prime}\frac{\partial u_j^{\prime}}{\partial x_i}}\right]\nonumber\\ &\quad+\overline{\left(\frac{\partial \beta^{\prime}}{\partial x_j}\right)\left[\frac{\partial}{\partial x_j}\left(u_i^{\prime}u_k^{\prime}\right)+u_i^{\prime}\frac{\partial u_j^{\prime}}{\partial x_k}+u_k^{\prime}\frac{\partial u_j^{\prime}}{\partial x_i}\right]}+ \overline{\beta^{\prime}\frac{\partial^{2}}{\partial x_j \partial x_j}\left(u_i^{\prime}u_k^{\prime}\right)} \nonumber\\ &\quad -2\overline{\beta^{\prime}\frac{\partial u_i^{\prime}}{\partial x_j}\frac{\partial u_k^{\prime}}{\partial x_j}}+ \overline{\left(\frac{\partial \beta^{\prime}}{\partial x_j}\right)\left[u_i^{\prime}\left(\frac{\partial \overline{u_k}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_k}\right)+u_k^{\prime}\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)\right]} \notag\\&\quad+\overline{\beta^{\prime}\left[\frac{\partial^{2}}{\partial x_j \partial x_j}\left(u_i^{\prime}\overline{u_k} \right) +\frac{\partial^{2}}{\partial x_j \partial x_j}\left(\overline{u_i}u_k^{\prime} \right) \right]}\nonumber\\ &\quad- 2 \overline{\beta^{\prime}\left[\frac{\partial u_i^{\prime}}{\partial x_j}\frac{\partial \overline{u_k}}{\partial x_j}+\frac{\partial \overline{u_i}}{\partial x_j}\frac{\partial u_k^{\prime}}{\partial x_j} \right]} - \overline{\beta^{\prime}\left[\overline{u_k}\frac{\partial^{2} u_i^{\prime}}{\partial x_j \partial x_j} + \overline{u_i}\frac{\partial^{2} u_k^{\prime}}{\partial x_j \partial x_j} \right]}. \end{align}

As can be seen from (2.10), several new terms have arisen due to the non-Newtonian rheology. To facilitate the later discussion, the new terms are labelled as follows:

1. (i)

   (2.11)\begin{equation} P_{ik} = - Re\left[\overline{u_i^{\prime}u_j^{\prime}} \frac{\partial \overline{u_k}}{\partial x_j} +\overline{u_k^{\prime}u_j^{\prime}} \frac{\partial \overline{u_i}}{\partial x_j}\right] \quad \Longrightarrow \quad \text{total production rate}, \end{equation}

2. (ii)

   (2.12)\begin{equation} {TT}_{ik} =- Re\frac{\partial}{\partial x_j} \left(\overline{u_i^{\prime}u_j^{\prime}u_k^{\prime}}\right) \quad \Longrightarrow \quad \text{turbulent transport rate}, \end{equation}

3. (iii)

   (2.13)\begin{equation} {MD}_{ik} = \bar{\beta} \frac{\partial^{2}}{\partial x_j \partial x_j}\left(\overline{u_i^{\prime}u_k^{\prime}}\right) \quad \Longrightarrow \quad \text{mean viscous diffusion rate}, \end{equation}

4. (iv)

   (2.14)\begin{equation} {\varPi}_{ik} = - Re\overline{\left(u_i^{\prime}\frac{\partial p^{\prime} } {\partial x_k}+u_k^{\prime}\frac{\partial p^{\prime} }{\partial x_i}\right)} \quad \Longrightarrow \quad \text{velocity--pressure gradient transport rate}, \end{equation}

5. (v)

   (2.15)\begin{equation} {M\epsilon}_{ik}=2\bar{\beta} \overline{\frac{\partial u_i^{\prime}}{\partial x_j}\frac{\partial u_k^{\prime}}{\partial x_j}} \quad \Longrightarrow \quad \text{mean viscous dissipation rate}, \end{equation}

6. (vi)

   (2.16)\begin{align} {Mv}_{ik}&=\left(\frac{\partial \bar{\beta}}{\partial x_j}\right)\left[\frac{\partial}{\partial x_j} \left(\overline{u_i^{\prime}u_k^{\prime}}\right)+ \overline{u_i^{\prime}\frac{\partial u_j^{\prime}}{\partial x_k}}+\overline{u_k^{\prime}\frac{\partial u_j^{\prime}}{\partial x_i}}\right]\nonumber\\ &\Longrightarrow \quad \text{mean viscosity gradient transport rate}, \end{align}

7. (vii)

   (2.17)\begin{align} {Tv}_{ik} &=\overline{\left(\frac{\partial \beta^{\prime}}{\partial x_j}\right)\left[\frac{\partial}{\partial x_j}\left(u_i^{\prime}u_k^{\prime}\right)+u_i^{\prime}\frac{\partial u_j^{\prime}}{\partial x_k}+u_k^{\prime}\frac{\partial u_j^{\prime}}{\partial x_i}\right]} \nonumber\\ &\Longrightarrow \quad \text{turbulent viscosity gradient transport rate}, \end{align}

8. (viii)

   (2.18)\begin{equation} {TD}_{ik} = \overline{\beta^{\prime}\frac{\partial^{2}}{\partial x_j \partial x_j}\left(u_i^{\prime}u_k^{\prime}\right)} \quad \Longrightarrow \quad \text{turbulent viscous diffusion rate}, \end{equation}

9. (ix)

   (2.19)\begin{equation} {T\epsilon}_{ik} = 2\overline{\beta^{\prime}\frac{\partial u_i^{\prime}}{\partial x_j}\frac{\partial u_k^{\prime}}{\partial x_j}} \quad \Longrightarrow \quad \text{turbulent viscous dissipation rate}, \end{equation}

10. (x)

    (2.20)\begin{align} {\widetilde{Tv}}_{ik} &= \overline{\left(\frac{\partial \beta^{\prime}}{\partial x_j}\right)\left[u_i^{\prime}\left(\frac{\partial \overline{u_k}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_k}\right)+u_k^{\prime}\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)\right]} \nonumber\\ &\Longrightarrow \quad \text{turbulent viscosity gradient transport rate related to mean flow}, \end{align}

11. (xi)

    (2.21)\begin{align} {\widetilde{TD}}_{ik} &= \overline{\beta^{\prime}\left[\frac{\partial^{2}}{\partial x_j \partial x_j}\left(u_i^{\prime}\overline{u_k} \right) +\frac{\partial^{2}}{\partial x_j \partial x_j}\left(\overline{u_i}u_k^{\prime} \right)-\overline{u_k}\frac{\partial^{2} u_i^{\prime}}{\partial x_j \partial x_j} - \overline{u_i}\frac{\partial^{2} u_k^{\prime}}{\partial x_j \partial x_j} \right]}\nonumber\\ &\Longrightarrow \quad \text{turbulent viscous diffusion rate related to mean flow},\ \textrm{and} \end{align}

12. (xii)

    (2.22)\begin{align} {\widetilde{T\epsilon}}_{ik} &= 2\overline{\beta^{\prime}\left[\frac{\partial u_i^{\prime}}{\partial x_j}\frac{\partial \overline{u_k}}{\partial x_j}+\frac{\partial \overline{u_i}}{\partial x_j}\frac{\partial u_k^{\prime}}{\partial x_j} \right]}\nonumber\\ &\Longrightarrow \quad \text{turbulent viscous dissipation rate related to mean flow}. \end{align}

The overall equation is recast as

(2.23)\begin{equation} B_{ik} = P_{ik} + T_{ik} - \epsilon_{ik}, \end{equation}

where total dissipation $\epsilon _{ik}$ and transport $T_{ik}$ rates, contributing to budget $B_{ik} = Re\,\mathrm {D}(\overline {u_i^{\prime }u_k^{\prime }}) / \mathrm {D}t$, are given by

(2.24)\begin{equation} \epsilon_{ik}={M\epsilon}_{ik} + {T\epsilon}_{ik} + {\widetilde{T\epsilon}}_{ik} \end{equation}

and

(2.25)\begin{equation} T_{ik}= {TT}_{ik} + {MD}_{ik} + {\varPi}_{ik} + {Mv}_{ik} + {Tv}_{ik} + {TD}_{ik} +{\widetilde{Tv}}_{ik} + {\widetilde{TD}}_{ik}, \end{equation}

respectively. The transport equation for the turbulent kinetic energy is found by taking the summation of the diagonal components of (2.23) and dividing the resulting expression by $2$. Also, the transport equation for the velocity correlation $\overline {u_i^{\prime }u_k^{\prime }}$ corresponding to a Newtonian fluid may be recovered from (2.23) by considering constant viscosity.

It is worth pointing out that, in the following sections, the results are presented in inner (viscous) units unless otherwise specified. Consequently, all variables are non-dimensionalized with $\mu _c = \mu _w$, $L_c = (\mu _w/\rho )/u_{\tau }$ and $U_c = u_{\tau }$ and are identified as $\{~\}^{+}$ quantities.

2.2. Numerical set-up, computational domain and grid resolution

Direct numerical simulations of a statistically converged plane turbulent channel flow of GN fluids at a low frictional Reynolds number, $Re_{\tau } = \rho u_{\tau } h/\mu _w$, have been performed. A turbulent channel flow is statistically stationary and homogeneous in the spanwise and streamwise directions. Since the flow is restricted by the channel walls, non-slip and impermeability boundary conditions are imposed at the walls. In the homogeneous directions, periodic boundary conditions may be employed if the domain is large enough to contain the largest structures in the flow.

DNS requires resolution of all spatial and temporal scales within the flow. Typically, a domain is considered sufficiently large if two-point correlations of the turbulent fluctuations decay close to zero at a separation of half the period in the homogeneous directions (see Moin & Mahesh 1998). Regarding the grid resolution for DNS, it should be fine enough to capture the smallest scales in proximity to the wall. Generally, the resolution is $O(\eta )$, $\eta$ being the Kolmogorov length scale. In Kim, Moin & Moser (1987) and Moser, Kim & Mansour (1999), the grid resolution is considered adequate if there is an evident scale separation, i.e. energy density at high wavenumbers is several decades lower than the one at low wavenumbers and if no energy pile-up is happening at the smallest scales.

The numerical simulations have been carried out using a finite-volume method on a collocated grid. Central differencing is used for the spatial discretization whilst the Crank–Nicolson scheme is employed for the discretization in time. The numerical procedure is based on an implicit, two-time-step advancement technique where the Poisson equation for the pressure is solved with an efficient multigrid method (see Emvin 1997). For more details regarding the numerical procedure and the FORTRAN 77 code, called CALC-LES, see, for instance, Davidson & Peng (2003) and Davidson (2018). The GN fluid rheology has been incorporated into the code through the Carreau fluid model and the simulation parameters are summarized in table 1. A viscosity rheogram, $\beta$ versus $\dot {\gamma }^{+}$, is also shown in figure 1. For all GN fluid cases, a target $Re_{\tau } = 180$ is considered. For cases N180 and D180, a computational box as in Kim et al. (1987) is used. For case P180, a larger computational domain is employed since structures of larger size than those in the Newtonian case are expected. See, for instance, Rudman & Blackburn (2006) and Singh et al. (2017b).

Table 1. Parameters of the simulation. Here $L_x$ and $L_z$ are the periodic streamwise and spanwise dimensions of the computational domain, $L_y =2h$ is the domain in the wall-normal direction, $\Delta x^{+}$ and $\Delta z^{+}$ are the constant grid spacings in inner units corresponding to the streamwise and spanwise directions whilst $\Delta {y_{min}^{+}}$ and $\Delta {y_{max} ^{+}}$ are the minimum and maximum grid spacings in inner units corresponding to the wall-normal direction. P180 and D180 refer to the shear thinning or pseudo-plastic and shear thickening or dilatant fluid cases whilst N180 refers to the base Newtonian case.

Figure 1. Viscosity rheogram. The flow index $\alpha$ is equal to $0.8$ and to $1.2$ for cases P180 and D180, respectively, whilst $\varLambda (U_c/L_c)=0.1$ for both non-Newtonian cases. The profiles corresponding to P180 and D180 are identified by red and cyan colours, respectively.

Synthetic turbulence (Davidson 2007) has been used to initialize case N180, whilst the simulations for the non-Newtonian cases have been initialized using a flow field corresponding to the Newtonian case at the same target $Re_{\tau }$ and the initial transients have been discarded. Statistically steady state is considered to have been achieved once a linear profile for $\bar {\tau }$ is observed (Vinuesa et al. 2016). The initial transient time is over $200$ convective time units, and $150$ flow fields saved every $\approx 0.2$ eddy turnover times have been considered to compute the statistics.

With respect to the computational domain, for all cases, the two-point correlations between turbulent fluctuations, i.e. $R_{ij} = \overline {u_i^{\prime }({\boldsymbol {x}}, t ) u_j^{\prime }({\boldsymbol {x}}+{\boldsymbol {r}}, t )} / \, \overline {u_i^{\prime }({\boldsymbol {x}},t) u_j^{\prime }({\boldsymbol {x}},t)}$, where ${\boldsymbol {r}}$ is the separation vector between the two points, appear to be decreasing with increasing separation in the homogeneous directions and the employed domains are deemed adequate; see figure 2. Note that, since the two-point correlations are related to the integral length scales, as expected, more elongated/larger structures are present in case P180 and finer structures in case D180 when compared with the Newtonian base case. Such observation is further supported by the instantaneous flow structures shown in figure 3. As can be seen, longer low-speed streamwise velocity streaks are observed for the shear-thinning case. Figure 3 also suggests that, for all fluid cases, the turbulent flow regime is achieved despite the low frictional Reynolds number used for the simulations.

Figure 2. Normalized two-point correlations $R_{ij}$: streamwise separation at (a) $y^{+} \approx 5.07$ and (b) $y^{+} \approx 154.96$, spanwise separation at (c) $y^{+} \approx 5.07$ and (d) $y^{+} \approx 154.96$. Line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used to identify $R_{11},R_{22}$ and $R_{33}$, respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

Figure 3. Contours of instantaneous streamwise velocity fluctuations normalized by frictional velocity at $y^{+} \approx 5.07$ (a,c,e) and $y^{+} \approx 154.96$ (b,d,f): (a,b) P180, (c,d) N180 and (e,f) D180. White and black represent positive and negative fluctuations, respectively.

To check the grid resolution, pre-multiplied one-dimensional spectral energy densities based on the presented two-point correlations are considered; see figure 4. Here, for all cases, an evident energy drop-off is occurring as the wavenumber increases and the maximum value in the different spectral energy densities is expected close to the wavenumbers corresponding to the respective integral length scales. For instance, the maximum in the pre-multiplied energy spectrum based on the two-point correlation for $u^{\prime }$ is expected at a wavenumber corresponding to the longitudinal integral length scale. Note as well that the pre-multiplied spectra corresponding to case P180, when compared to the results from cases N180 and D180 for the same number of modes, are presenting lower maximum amplitudes due to the increase in the computational domain for that case.

Figure 4. Pre-multiplied one-dimensional spectral energy density: against normalized streamwise wavenumber $k_xh$ at (a) $y^{+} \approx 5.07$ and (b) $y^{+} \approx 154.96$, and against normalized spanwise wavenumber $k_zh$ at (c) $y^{+} \approx 5.07$ and (d) $y^{+} \approx 154.96$. Line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used to identify results corresponding to $E_{11},E_{22}$ and $E_{33}$ spectral energy densities, respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

In addition, whilst considering the computational resolution, a length scale based on the mean viscosity and the total dissipation rate has been computed; see figure 5. For the Newtonian case, this is the Kolmogorov length scale. For cases P180 and D180, since Kolmogorov's first similarity hypothesis (see e.g. Pope 2000) is stated for constant viscosity, such length scale $\tilde {\eta } = (\bar {\mu }/\rho )^{3/4}/\epsilon _k^{1/4}$, where $\epsilon _k=(\epsilon _{11}+\epsilon _{22}+\epsilon _{33})/2$ is the total mean dissipation rate (see Bradshaw & Perot (1993) and Bradshaw (1995) for a note about true dissipation), is analogous to $\eta$ in the regions with minor variation in the mean viscosity profile, i.e. in the near-wall region and close to the channel's centre (see § 3.1). As seen from figure 5, for N180, $\Delta {y_{min} ^{+}}< \eta ^{+}$ in the near-wall region and $\Delta {y_{max} ^{+}}\approx O(\eta ^{+})$ at the channel's centre. Similar trends are noted for the non-Newtonian fluid cases when comparing the wall-normal grid resolution with $\tilde {\eta }^{+}$.

Figure 5. Kolmogorov's length scale $\eta ^{+}$ and analogous length scale ${\tilde {\eta }}^{+}$, for the non-Newtonian cases, against $y^{+}$. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

Based on the energy drop-off previously noted and on the ${\eta }^{+}$ and $\tilde {\eta }^{+}$ values, the computational resolution appears to be adequate for the different cases. Also, a verification with published data for the Newtonian case is presented in appendix B.

3. Results

3.1. Low-order statistics

The mean (averaged in time and in homogeneous directions) streamwise velocity profile ${\bar {u}}^{+}$ and its diagnosis function $y^{+} \mathrm {d} {\bar {u}}^{+} / \mathrm {d}y^{+}$ are presented in figure 6. In the same figure, the mean viscosity profile ${\bar {\mu }}^{+} = \bar {\beta }$ and an analogous quantity to the diagnosis but for the mean viscosity, $y^{+} \mathrm {d} {\bar {\beta }} / \mathrm {d}y^{+}$, are presented as well.

Figure 6. First-order statistics: (a) ${\bar {u}}^{+}$, (b) $y^{+} \mathrm {d} {\bar {u}}^{+} / \mathrm {d}y^{+}$, (c) $\bar {\beta }$ and (d) $y^{+} \mathrm {d}{\bar {\beta }}/ \mathrm {d}y^{+}$ versus $y^{+}$. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

To discuss variations against the wall-normal coordinate $y^{+}$, the classical flow-region subdivision based on case N180 is considered: there is an inner region, comprising a viscous sublayer $(y^{+} \lesssim 5 )$; a buffer region $(5\lesssim y^{+} \lesssim 55 )$; a quite limited – if it exists at all – log region $(55\lesssim y^{+} \lesssim 62 )$; and a remaining outer region. Within the viscous sublayer, for all cases, ${\bar {u}}^{+} \approx y^{+}$ and for $y^{+} \lesssim 10$, minor deviation for the non-Newtonian cases is noted in ${\bar {u}}^{+}$ and $\bar {\beta }$ when compared to case N180. For $y^{+} \gtrsim 10$, ${\bar {u}}^{+}$ increases with decreasing flow index $\alpha$, and a larger bulk velocity, implying a lower friction factor for a common driving pressure gradient, is observed. Also, as evidenced by the diagnosis function, the starting point of the log region and the corresponding $y^{+} \mathrm d {\bar {u}}^{+} / \mathrm {d}y^{+} \approx 1/\kappa$, where $\kappa$ is the von Kármán constant, increase with decreasing flow index.

Regarding $\bar {\beta }$, for $y^{+} \gtrsim 10$, it deviates rapidly from the corresponding constant Newtonian value, as expected. Since $\bar {\beta }$ appears to behave in a log manner within a certain $y^{+}$ range, an analogous diagnosis function based on $\bar {\beta }$ has been considered. On a first impression, such a log region in the mean viscosity $(y^{+} \mathrm {d}{\bar {\beta }}/\mathrm {d}y^{+} \approx \varGamma )$ indeed occurs but different slopes $\varGamma$ for cases P180 and D180 are observed; see figure 6(d).

For the outer layer region, the velocity defect profile ${\bar {u}}^{+} (y/h=1 )-{\bar {u}}^{+}$ shown in figure 7(a) is apprised. From the velocity defect profile, it is clear that larger mean centreline velocity values are observed with decreasing flow index. Also, it appears that, in the outer layer region $(y^{+} \gtrsim 62$ or $y/h \gtrsim 0.35)$, all curves collapse. Such behaviour suggests independence of the velocity defect profile in the outer region and that differences in the mean velocity profiles between the different cases (P180, N180 and D180) are most likely due to differences in the flow within the inner region. This observation has also been made by Singh et al. (2018).

Figure 7. (a) Velocity defect profile against $y/h$ and (b) low-order statistics corresponding to the strain rate against $y^{+}$; ‘—–’ mean strain rate ${\bar {\dot {\gamma }}}^{+}$ and ‘- - -’ r.m.s. values of strain-rate fluctuations r.m.s.$({\dot {\gamma }}^{+})$. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

Distribution in the $i$ direction ($i=1,2,3$) of the r.m.s. values of the velocity fluctuations and the vorticity fluctuations, i.e. r.m.s.$(u_i^{+})={\overline {(u_i^{\prime })^{2}}}^{1/2} / u_{\tau }$ and r.m.s.$(\omega _i^{+})={\overline {(\omega _i^{\prime })^{2}}}^{1/2} [\mu _w /(\rho u_{\tau }^{2})]$, respectively, is presented in figure 8. The r.m.s. values corresponding to fluctuations in viscosity r.m.s.$(\beta )={\overline {(\mu \prime )^{2}}}^{1/2} / \mu _w$ are also shown in figure 8. The r.m.s. values for the velocity fluctuations appear to be affected by the rheology mainly outside the viscous sublayer. The streamwise turbulence intensity increases with decreasing flow index whilst the opposite is observed in the spanwise and wall-normal turbulence intensities for the same trend. The wall-normal position at which a peak in the streamwise turbulence intensity appears moves slightly from the wall with decreasing flow index.

Figure 8. Second-order statistics: (a) turbulence intensities r.m.s.$(u_i^{+})$, (b) r.m.s. values of vorticity fluctuations r.m.s.$(\omega _i^{+})$ and (c) r.m.s. value of viscosity fluctuations r.m.s.$(\beta )$. In (a,b) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used to identify r.m.s. values based on $x,y$ and $z$ velocity and vorticity fluctuation components. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

With respect to the r.m.s. values of the vorticity fluctuations, in the near-wall region, there is an increase in the spanwise component whilst a decrease is observed for the other two components with decreasing flow index. Thus, we anticipate an increase in the magnitude of the mean viscous dissipation since, for homogeneous turbulence, the enstrophy is approximately equal to the ratio between mean viscous dissipation rate and mean kinematic viscosity. At a wall-normal position within the buffer region $(y^{+} \approx 10 )$, the opposite is then observed for the r.m.s. values of the spanwise vorticity fluctuation, i.e. the r.m.s. value decreases with shear thinning. Moreover, after such a point, the r.m.s. values of all vorticity fluctuations appear to decrease with decreasing flow index.

It is interesting to note that the presented results suggest an overall increase in the anisotropy of the velocity and vorticity correlation tensors, $\overline {u_i^{\prime } u_k ^{\prime }}$ and $\overline {\omega _i^{\prime } \omega _k ^{\prime }}$, with shear thinning (see § 3.5).

Regarding the r.m.s. values corresponding to fluctuations in viscosity, similar to the mean viscosity profile, the distribution appears to be approximately constant in the viscous sublayer and then it starts to increase with $y^{+}$ up to a certain value within the inner region. After such a wall-normal position, the r.m.s. values of the viscosity fluctuations start to decrease more rapidly. Here, the plateau in $\bar {\beta }$ and r.m.s.$(\beta )$ within the viscous sublayer is likely due to small variations in the mean strain rate and its r.m.s. values, respectively (see figure 7b). Note that, altogether, r.m.s.$(\beta )$ is larger for case P180, suggesting larger fluctuations from $\bar {\beta }$ with shear thinning. It is also noted that the peak in r.m.s.$(\beta )$ moves from the wall with shear thinning.

Finally, it is worth pointing out some qualitative similarities with channel flow of viscoelastic fluids, which is another type of drag-reducing fluid. Compared to the Newtonian case and as with shear thinning, there is a noticeable increase in the mean streamwise velocity with viscoelasticity. Also for $Re_{\tau } <1000$, a minimal – if any – presence of a log-law region is observed (Thais, Gatski & Mompean 2012). With respect to the turbulence intensities, compared to the Newtonian case and as with shear thinning, Ptasinski et al. (2003) and Min, Choi & Yoo (2003a) reported an enhancement of the streamwise turbulence intensity and a decrease of the wall-normal and spanwise turbulence intensities with viscoelasticity and in the case of small drag reduction (SDR) regime (see Warholic, Massah & Hanratty 1999). In contrast, the same publications reported a decrease in all turbulence intensities with viscoelasticity for the large drag reduction (LDR) regime. On the other hand, regarding low-order statistics reported for other canonical flows of GN fluids, it is worth mentioning that Singh et al. (2017b) showed similar trends with shear thinning for pipe flow. For $\alpha <1$, there is an increase in the mean axial velocity, mean viscosity and axial turbulence intensity whereas both radial and azimuthal turbulence intensities decrease.

3.2. Mean shear stress budget

For a statistically stationary fully developed flow of an incompressible GN fluid, the total mean shear stress $\bar {\tau }^{+}$ (see § 2.1.1) is given by

(3.1)\begin{equation} \bar{\tau}^{+} = \bar{\tau}^{+}_{{vis}} + \bar{\tau}^{+}_{{tur}} + \bar{\tau}^{+}_{{GN}}, \end{equation}

where $\bar {\tau }^{+}_{{vis}} = \bar {\beta } \mathrm {d}\bar {u}^{+}/\mathrm {d}y^{+}$, $\bar {\tau }^{+}_{{tur}}=-\overline {u^{\prime } v^{\prime }}/u_{\tau }^{2}$ and $\bar {\tau }^{+}_{{GN}} = \overline {2\beta ^{\prime } {{\mathsf {S}}_{12}^{\prime }}^{+}}$ are the viscous stress, the turbulent or Reynolds shear stress and the new stress due to fluctuation in viscosity, respectively. Owing to the constant pressure gradient driving the flow, the total mean shear stress is linear, i.e. $\bar {\tau }^{+} = 1-y^{+}/Re_{\tau }$.

The different contributions to the total mean shear rate are shown in figure 9. As can be seen, the mean viscous stress increases in the shear-thinning case. Since $\bar {\tau }^{+}_{{vis}}$ depends on $\bar {\beta }$ and $\mathrm {d}\bar {u}^{+}/\mathrm {d}y^{+}$, its increase can be discussed considering figure 6(b,c). For $y^{+} \lesssim 10$, there is minor variation in the mean viscosity and in the mean streamwise velocity gradient. Thus, the increase in $\bar {\tau }^{+}_{{vis}}$ for the shear-thinning case is attributed to a small increase in both quantities. For $y^{+} > 10$, outside the viscous sublayer, the increase in the mean viscous stress is mostly due to $\bar {\beta }$, which quite significantly increases in the shear-thinning case due to the higher ‘zero’ shear-rate viscosity (see viscosity rheogram, figure 1). In the near-wall region where the turbulent stress is close to zero, the decrease/increase in the viscous stress component due to shear thickening/shear thinning is compensated by the new stress component. Note that the sign of $\bar {\tau }^{+}_{{GN}}$ depends on the sign of $\alpha -1$. Outside the viscous sublayer, the new stress starts to decrease and the viscous stress is then compensated by the turbulent stress component. Thus, for example, an increase in the viscous stress due to shear thinning is compensated by a decrease in the Reynolds shear stress. It is also interesting to note that $\bar {\tau }^{+}_{{tur}}$, for all cases, appear to collapse in the outer region, suggesting that the cross-correlation $\overline {u^{\prime } v^{\prime }}/u_{\tau }^{2}$ is independent of the rheology in that region.

Figure 9. Mean shear stress budget. Line styles ‘—–’, ‘- - -’, ‘$\cdot \cdot \cdot$’ and ‘- $\cdot$ -’ are used to identify $\bar {\tau }^{+}$, $\bar {\tau }^{+}_{{vis}}$, $\bar {\tau }^{+}_{{tur}}$ and $\bar {\tau }^{+}_{{GN}}$, respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

Finally, it is worth contrasting the mean shear stress budget corresponding to the shear-thinning fluid with the budget of a drag-reducing polymer solution, i.e. a viscoelastic fluid. Compared to the Newtonian case and as with shear thinning, there is a significant overall reduction in the Reynolds shear stress with polymer additives (Min et al. 2003a; Ptasinski et al. 2003). In addition, Ptasinski et al. (2003) reported that the polymer stress, arising due to viscoelastic effects, is always positive but relatively small and mainly important near to the wall for SDR, whereas for LDR, this contribution is large and important across the whole channel. Note that, in contrast, the analogous $\bar {\tau }^{+}_{{GN}}$ is always negative for shear thinning and mainly important close to the wall. On the other hand, regarding Reynolds stress budgets reported for other canonical flows of GN fluids, Singh et al. (2017b) presented similar trends with shear thinning for pipe flow. For $\alpha <1$, the mean viscous stress slightly increases near the wall and more noticeably in the buffer layer region whereas the Reynolds shear stress decreases for all $y^{+}$ and a new negative stress arises in the total mean stress balance.

3.3. Quadrant analysis

To improve our understanding of the generation of Reynolds shear stress and the related production of turbulent kinetic energy, a quadrant analysis (Wallace et al. 1972) is carried out. Contributions to the cross-correlation $-\overline {u^{\prime } v^{\prime }}$ are classified according to the sign of the velocity fluctuations into four categories or quadrants: Q1 $(u^{\prime } > 0,\ v^{\prime } > 0)$, Q2 $(u^{\prime } < 0,\ v^{\prime } > 0)$, Q3 $(u^{\prime } < 0,\ v^{\prime } < 0)$ and Q4 $(u^{\prime } > 0,\ v^{\prime } < 0)$, and each of the quadrant motions is associated with a physical event: positive production of turbulent kinetic energy arises due to low-speed fluid moving from the wall (Q2 events) and high-speed fluid moving towards the wall (Q4 events). Such motions have been visualized, see for instance Kline et al. (1967) and Corino & Brodkey (1969), and denoted as ejection and sweep events, respectively. Q1 and Q3 motions, which correspond to high-speed fluid reflected outwards from the wall and low-speed fluid deflected towards the wall, account for negative production of turbulent kinetic energy, and in the absence of a better terminology may be called outwards and wallwards interactions (Wallace 2016).

The quadrant-conditioned contributions to the Reynolds shear stress are presented in figure 10 for the different cases. In the following discussion, changes in quadrant event contributions with shear thinning are also considered in light of variation with decreasing flow index in the skewness and flatness profiles corresponding to the velocity fluctuations. The profiles are shown in figure 11.

Figure 10. Quadrant contributions to $-\overline {u^{\prime } v^{\prime }}/u_{\tau }^{2}$: (a) non-fractional $-\overline {u^{\prime } v^{\prime }}_{q}/u_{\tau }^{2}$ and (b) fractional $\overline {u^{\prime } v^{\prime }}_{q} /\overline {u^{\prime } v^{\prime }}$ contributions. The subscript $q$ refers to quadrant-conditioned averages. Line styles ‘- - -’, ‘- $\cdot$ -’, ‘$\cdot \cdot \cdot$’ and ‘—–’ are used for Q1, Q2, Q3 and Q4 contributions, respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

Figure 11. Skewness $S$ and flatness $F$ profiles of velocity fluctuations: (a) $S(u^{+})=\overline {(u^{\prime })^{3}}/(\overline {u^{\prime } u^{\prime }})^{3/2}$, (b) $F(u^{+})=\overline {(u^{\prime })^{4}}/(\overline {u^{\prime } u^{\prime }})^{2}$, (c) $S(v^{+})=\overline {(v^{\prime })^{3}}/(\overline {v^{\prime } v^{\prime }})^{3/2}$, (d) $F(v^{+})=\overline {(v^{\prime })^{4}}/(\overline {v^{\prime } v^{\prime }})^{2}$, (e) $S(w^{+})=\overline {(w^{\prime })^{3}}/(\overline {w^{\prime } w^{\prime }})^{3/2}$ and (f) $F(w^{+})=\overline {(w^{\prime })^{4}}/(\overline {w^{\prime } w^{\prime }})^{2}$. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

From the non-fractional contributions to $-\overline {u^{\prime } v^{\prime }}/u_{\tau }^{2}$, a decrease in contributions from all quadrants’ events with shear thinning is observed. This trend is consistent with the perceived decrease in Reynolds shear stress with shear thinning seen in § 3.2. Also, from the non-fractional contributions, it is noted that the wall-normal position at which contributions from sweep and ejection events coincide has moved slightly from the wall for the shear-thinning case. A similar behaviour – with shear thinning – is then expected for the peaks corresponding to maximum production and maximum turbulent kinetic energy (see § 3.4).

Observing the fractional quadrant contributions, it is noticed that, for all fluid cases, sweep events appear to dominate in the very near-wall region whilst ejection events contribute more to $-\overline {u^{\prime } v^{\prime }}/u_{\tau }^{2}$ after the $y^{+}$ position where the contributions due to Q2 and Q4 motions are approximately the same. Within the viscous sublayer where quite large intermittency is present, a more pronounced increase in the fractional Q1 contribution is observed for the shear-thinning case. Such an increase, mainly compensated by fractional contributions due to Q4 events, is expected since the correlation between $u^{\prime } > 0$ and $v^{\prime } > 0$ appears to vary little with rheology (see non-fractional contribution profiles) and the overall Reynolds shear stress is decreasing with shear thinning even within such a very near-wall region. As reflected in the normalized skewness profiles, large positive streamwise and wall-normal velocity fluctuations appear to dominate more with shear thinning and, as shown in the non-fractional contributions, lead to a decrease in contributions from Q2 and Q4 events, which decreases the overall Reynolds shear stress.

Outside the viscous sublayer but before the wall-normal position where Q2 and Q4 contributions are approximately the same, there is little variation in the non-fractional contributions due to Q1 and Q3 events with rheology. The decrease in contribution from sweep and ejection events, together causing the decrease in $-\overline {u^{\prime } v^{\prime }}/u_{\tau }^{2}$, appears to be related to the appearance of more dominant large positive streamwise fluctuation and less dominant negative wall-normal fluctuations with shear thinning, as reflected in the normalized skewness profiles. After the point with equal contributions from Q2 and Q4 events but before a position $y^ + \approx 30\text {--}35$ where ejection events dominate due to large intermittent $v^{\prime } > 0$ values for all GN fluids, the decrease in Reynolds shear stress with shear thinning appears to be linked to a slight increase in dominant $u^{\prime } > 0$ and $v^{\prime } < 0$ signals.

Finally, after the position $y^ + \approx 30\text {--}35$, non-fractional contributions due to Q1, Q2 and Q3 events almost do not vary with rheology, and the observed decrease in $-\overline {u^{\prime } v^{\prime }}/u_{\tau }^{2}$ for the shear-thinning case is attributed to a decrease in the contribution from sweep events with shear thinning. Here, such behaviour is likely to be due to more dominant $u^{\prime } < 0$ and $v^{\prime } > 0$ signals with decreasing flow index, as reflected in the corresponding normalized skewness profiles.

3.4. Reynolds stress budgets

The Reynolds stress budgets, corresponding to the equations presented in § 2.1.2, are considered. The total production, transport and dissipation rates corresponding to budgets $B_{ik}^{+}$ of the relevant non-zero stresses in a fully developed turbulent channel flow, for the different GN fluid cases, are shown in figure 12. The turbulent kinetic energy budget $B_{k}^{+}$ and the turbulent kinetic energy profile $k^{+} = \overline {u_i^{\prime } u_i^{\prime }}/(2 u_{\tau }^{2})$ are presented in figure 13.

Figure 12. Reynolds stress budgets: (a) $B_{11}^{+}$, (b) $B_{12}^{+}$, (c) $B_{22}^{+}$ and (d) $B_{33}^{+}$. The line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used to identify $P_{ik}^{+}$, $T_{ik}^{+}$ and $-\epsilon _{ik}^{+}$, respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

Figure 13. Plots of (a) $B_{k}^{+}$ and (b) $k^{+}$. In panel (a), the line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used to identify $P_{k}^{+}$, $T_{k}^{+}$ and $-\epsilon _{k}^{+}$, respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

For all cases, and as expected, the interaction between $-\overline {u^{\prime } v^{\prime }}/u_{\tau } ^{2}$ and mean shear causes a production rate in budget $B_{11}^{+}$. Since the mean velocity gradients increase with shear thinning (see $\bar {u}^{+}$ profile in § 3.1), the decrease in $P_{11}^{+}$, in its near-peak region, is attributed to the observed decrease in Reynolds shear stress with decreasing flow index (see $\bar {\tau }_{{tur}}^{+}$ in § 3.2). Note as well that, as commented in § 3.3, the peak in production at $y^{+} \approx 12$ for the Newtonian case has moved slightly from the wall with shear thinning. In budget $B_{11}^{+}$, the energy from the production region is then distributed through $T_{11}^{+}$ towards and away from the wall. At the wall, a good portion of the transported energy from the maximum production region is then dissipated at a rate $\epsilon _{11}^{+}$. Everywhere $P_{11}^{+}$ is balanced by $T_{11}^{+} -\epsilon _{11}^{+}$. Note that, consistently with the decrease in the production rate with shear thinning, the amount of energy that is irreversibly dissipated decreases with decreasing flow index.

Since budgets $B_{22}^{+}$ and $B_{33}^{+}$ do not contain a production term, their source is energy being redistributed from budget $B_{11}^{+}$. Here, an apparent decrease in redistribution of energy from $B_{11}^{+}$ to $B_{22}^{+}$ and $B_{33}^{+}$, reflected in the decrease of $T_{11}^{+}$, $T_{22}^{+}$ and $T_{33}^{+}$, is occurring with decreasing flow index. This observation is consistent with the noticed increase in streamwise turbulence intensity and the decrease in the spanwise and wall-normal intensities with shear thinning, seen in § 3.1.

For budget $B_{12}^{+}$, the interaction between wall-normal turbulence intensities and mean shear yields production. Owing to the aforementioned decrease in the wall-normal intensities with decreasing flow index, there is a decrease in the production rate $P_{12}^{+}$ with shear thinning. This observation is consistent with the noted decrease in Reynolds shear stress with decreasing flow index. As for the $B_{11}^{+}$ budget, the source from the production region is transported towards and away from the wall through $T_{12}^{+}$. However, different from that budget, little of that source (even less so for the shear-thinning fluid case) is ‘dissipated’ at a rate $\epsilon _{12}^{+}$ and mostly is balanced by $T_{12}^{+}$.

The balance of the turbulent kinetic energy budget $B_{k}^{+}$ is based on the terms appearing in $B_{11}^{+}$, $B_{22}^{+}$ and $B_{33}^{+}$, where budget $B_{11}^{+}$ with turbulent production dominates over the others. Thus, the balance of $B_{k}^{+}$ and the corresponding observations are similar to those made for budget $B_{11}^{+}$. It is interesting to note that most variations in the terms of $B_{k}^{+}$ with rheology appear to be restricted to the inner layer region. Also, with respect to the total dissipation rate, since from the r.m.s. values of the vorticity fluctuations (see § 3.1) an increase in the magnitude of the mean viscous dissipation rate is expected with shear thinning, the overall decrease in the magnitude of $\epsilon _{k}^{+}$ appears to be due to the non-Newtonian terms contributing to it.

Regarding $k^{+}$, as expected, there is an increase with shear thinning and, similar to the production rate, its peak has moved slightly from the wall with decreasing flow index. The increase in turbulent kinetic energy may be explained by considering the deficit in redistribution of energy from budget $B_{11}^{+}$ to $B_{22}^{+}$ and $B_{33}^{+}$ with shear thinning, since it causes the observed increase in anisotropy between the turbulence intensities. The profile of $k^{+}$ may also be understood while examining the total production and dissipation rates. Within the region where production exceeds dissipation and the transport rate $T_{k}^{+}$ is negative, there is an increase in how much $P_k^{+}$ exceeds $-\epsilon _k^{+}$ with decreasing flow index. This is reflected by an increase in the turbulent kinetic energy profile within the same region.

Contributions from the different terms in (2.23) to total transport and dissipation rates for budgets $B_{11}^{+}$, $B_{12}^{+}$ and $B_{k}^{+}$ are shown in figures 14–16. Contributions to budgets corresponding to the other diagonal components of $\overline {u_i^{\prime } u_k^{\prime }}$, i.e. $B_{22}^{+}$ and $B_{33}^{+}$, are presented in appendix C.

Figure 14. Contributions to $T_{11}^{+}$ and $-\epsilon _{11}^{+}$ in budget $B_{11}^{+}$: (a) line styles ‘—–’ and ‘- - -’ are used for $\varPi _{11}^{+}$ and ${TT}_{11}^{+}$; (b) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used for ${MD}_{11}^{+}$, ${\widetilde {TD}_{11}}^{+}$ and ${TD}_{11}^{+}$; (c) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used for $-{M\epsilon }_{11}^{+}$, $-{\widetilde {T\epsilon }}_{11}^{+}$ and $-{T\epsilon }_{11}^{+}$; and (d) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used for ${\widetilde {Tv}}_{11}^{+}$, ${Mv}_{11}^{+}$ and ${Tv}_{11}^{+}$. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

For the transport rate $T_{11}^{+}$, the traditional terms $\varPi _{11}^{+}$, ${TT}_{11}^{+}$ and ${MD}_{11}^{+}$ contribute the most. The velocity–pressure gradient term $\varPi _{11}^{+}$ notably decreases in magnitude with decreasing flow index. This behaviour is to be expected since redistribution of energy to budgets $B_{22}^{+}$ and $B_{33}^{+}$ decreases with shear thinning. The mean viscous diffusion ${MD}_{11}^{+}$, which constitutes the largest contribution to $T_{11}^{+}$, appears to be mainly affected by rheology within the viscous sublayer and increases with decreasing flow index. Since the mean viscosity varies little for $y^{+} \lesssim 10$, such an increase is attributed to an increase of diffusion of $\overline {u^{\prime } u^{\prime }}/u_{\tau } ^{2}$ with shear thinning. The turbulent transport ${TT}_{11}^{+}$, on the other hand, is mainly affected outside the viscous sublayer but within the inner layer region by rheology. Here ${TT}_{11}^{+}$ mostly decreases in magnitude with decreasing flow index.

The non-Newtonian transport term mainly affecting $T_{11}^{+}$ within the viscous sublayer is ${\widetilde {TD}}_{11}^{+}$. Since, in the very near-wall region, the fluctuations in viscosity are rather small (see the r.m.s.$(\beta )$ profile in § 3.1), the diffusion of $u^{\prime }\bar {u}/u_{\tau }^{2}$ is considered the reason for the large magnitude in ${\widetilde {TD}}_{11}^{+}$ within this region. Also, the turbulent viscous diffusion rate related to mean flow decreases the total transport rate $T_{11}^{+}$ with shear thinning. Other transport terms such as ${TD}_{11}^{+}$, ${Tv}_{11}^{+}$ and ${Mv}_{11}^{+}$ related to fluctuation in viscosity, their gradients and gradients in the mean viscosity either decrease or increase with rheology at different wall-normal positions but almost do not affect the total transport rate. The remaining transport term ${\widetilde {Tv}}_{11}^{+}$ increases with decreasing flow index and appears to be mainly relevant in the inner region outside the viscous sublayer. In this region both gradients of fluctuations in viscosity and advected gradients of the mean streamwise velocity are expected to contribute to ${\widetilde {Tv}}_{11}^{+}$.

For the total dissipation rate $\epsilon _{11}^{+}$, ${M\epsilon }_{11}^{+}$ and ${\widetilde {T\epsilon }_{11}}^{+}$ appear to contribute the most. The mean viscosity increases with decreasing flow index; however, at the wall and in the very near-wall region, the mean viscosity is almost unaffected by rheology. Therefore, the increase of magnitude in ${M\epsilon }_{11}^{+}$ with shear thinning is mainly attributed to an increase in the autocorrelation of streamwise velocity fluctuation gradients. This observation is consistent with the increase in the r.m.s. values of the spanwise vorticity fluctuations at the wall and very near it with shear thinning (see § 3.1). Outside the viscous sublayer and mainly for $y^{+} \gtrsim 10$, the mean viscosity starts to increase in a more noticeable manner with decreasing flow index, and contributes more to the increase in magnitude of the mean viscous dissipation rate within that region.

The non-Newtonian dissipation rate ${\widetilde {T\epsilon }_{11}}^{+}$ increases in magnitude with decreasing flow index and mainly affects $\epsilon _{11}^{+}$ in the inner region, specially within the viscous sublayer. In the vicinity of the wall, the fluctuations in viscosity are fairly small and the large magnitude observed for ${\widetilde {T\epsilon }}_{11}^{+}$ is attributed to the correlation between mean streamwise velocity gradients and streamwise velocity fluctuation gradients. Note that the decrease in the magnitude of $\epsilon _{11}^{+}$ at the wall with shear thinning is mainly attributed to dissipation at rate ${\widetilde {T\epsilon }_{11}}^{+}$, since the magnitude of ${M\epsilon }_{11}^{+}$ is increasing with decreasing flow index. The other non-Newtonian dissipation rate ${T\epsilon }_{11}^{+}$ also causes the decrease in the magnitude of the total dissipation rate $\epsilon _{11}^{+}$ with shear thinning. Nonetheless, compared to ${\widetilde {T\epsilon }_{11}}^{+}$, it is smaller in magnitude. The dissipation rate ${T\epsilon }_{11}^{+}$ arises due to interactions between fluctuations in viscosity and squared streamwise velocity fluctuation gradients.

For the transport rate ${T}_{12}^{+}$, the important terms appear to be the velocity–pressure gradient term and the turbulent transport rate. Most of the production ${P}_{12}^{+}$ is balanced by $\varPi _{12}^{+}$ and ${TT}_{12}^{+}$. Within the inner region, the mean viscous diffusion rate appears to contribute in no significant manner to the rate ${T}_{12}^{+}$ and, even less so, other non-zero transport terms associated to the non-Newtonian rheology.

The total ‘dissipation’ rate $\epsilon _{12}^{+}$ is quite low and it is mainly due to the mean viscous ‘dissipation’ rate ${M\epsilon }_{12}^{+}$. The rate ${T\epsilon }_{12}^{+}$ presents the opposite sign to ${M\epsilon }_{12}^{+}$ for the same non-Newtonian fluid case, whilst the rate ${\widetilde {T\epsilon }_{12}}^{+}$ either increases or decreases the total rate $\epsilon _{12}^{+}$, depending on which part of the inner region is being considered. Nonetheless, none of the ‘dissipation’ rates appearing with non-Newtonian rheology contribute considerably to the total rate $\epsilon _{12}^{+}$.

Regarding the terms contributing to the total transport and dissipation rates in budget $B_{k}^{+}$, the observed trends with rheology are the same as those noted for $T_{11}^{+}$ and $\epsilon _{11}^{+}$, since the terms in budget $B_{11}^{+}$ contribute the most to the turbulent kinetic energy budget. The profile of rate $\varPi _{k}^{+}$ is, of course, different from the profile of the velocity–pressure gradient term in budget $B_{11}^{+}$, but the trend of a decrease in magnitude with shear thinning, due to less energy being redistributed, is the same.

Also, it is worth mentioning some similarities found with the budgets reported for another type of drag-reducing fluid, namely a viscoelastic fluid. For the channel flow of a viscoelastic fluid, Dimitropoulos et al. (2001) showed that the most significant changes in budget $B_{k}^{+}$ relative to Newtonian values are observed in the turbulence production, mean viscous dissipation and turbulence transport rates. Compared to the Newtonian case and as with shear thinning, there is a decrease in the production term and in the absolute peak values for the turbulent transport term with viscoelasticity. In addition, opposite to the behaviour seen for shear thinning, Dimitropoulos et al. (2001) reported a decrease in the mean viscous dissipation with shear thinning. The paper also shows that viscoelastic terms, contributing to $B_{k}^{+}$, are relatively small when compared to the production, mean diffusion or mean viscous dissipation rates. Regarding the individual budgets, Dimitropoulos et al. (2001) revealed that viscoelastic effects are most pronounced in the velocity–pressure gradient term. Compared to the Newtonian case and as with shear thinning, in the budget $B_{11}^{+}$, the rate $\varPi _{11}^{+}$ is reduced with viscoelasticity, in particular close to the region of maximum production. The paper also reported a significant reduction in rate $\varPi _{12}^{+}$ of budget $B_{12}^{+}$. Note that figure 15(a) shows a similar behaviour with shear thinning. On the other hand, regarding budgets reported for other canonical flows of GN fluids, Singh et al. (2017b) presented the budget $B_{k}^{+}$ for pipe flow and reported similar trends with shear thinning. For $\alpha <1$, there is a decrease of turbulence production in the buffer layer region, whereas the magnitude of the turbulence transport and mean viscous dissipation rates increase in the near-wall region. Also for shear thinning, in pipe flow, the magnitude of the total dissipation rate is decreased due to contributions from the non-Newtonian terms. It is worth commenting that an apparent inner region dependence on the flow index for the budget $B_{k}^{+}$ is also reported by Singh et al. (2017b) for pipe flow.

Figure 15. Contributions to $T_{12}^{+}$ and $-\epsilon _{12}^{+}$ in budget $B_{12}^{+}$: (a) line styles ‘—–’ and ‘- - -’ are used for $\varPi _{12}^{+}$ and ${TT}_{12}^{+}$; (b) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used for ${MD}_{12}^{+}$, ${\widetilde {TD}_{12}}^{+}$ and ${TD}_{12}^{+}$; (c) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used for $-{M\epsilon }_{12}^{+}$, $-{\widetilde {T\epsilon }}_{12}^{+}$ and $-{T\epsilon }_{12}^{+}$; and (d) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used for ${\widetilde {Tv}}_{12}^{+}$, ${Mv}_{12}^{+}$ and ${Tv}_{12}^{+}$. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

Figure 16. Contributions to $T_{k}^{+}$ and $-\epsilon _{k}^{+}$ in budget $B_{k}^{+}$: (a) line styles ‘—–’ and ‘- - -’ are used for $\varPi _{k}^{+}$ and ${TT}_{k}^{+}$; (b) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used for ${MD}_{k}^{+}$, ${\widetilde {TD}}_{k}^{+}$ and ${TD}_{k}^{+}$; (c) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used for $-{M\epsilon }_{k}^{+}$, $-{\widetilde {T\epsilon }}_{k}^{+}$ and $-{T\epsilon }_{k}^{+}$; and (d) line styles ‘—–’, ‘- - -’ and ‘$\cdot \cdot \cdot$’ are used for ${\widetilde {Tv}}_{k}^{+}$, ${Mv}_{k}^{+}$ and ${Tv}_{k}^{+}$. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

In summary, some of the most relevant findings are as follows.

1. (i) Rheological variations affecting the Reynolds shear stress and the mean shear lead to shear-dependent changes in the production rate $P_{11}^{+}$ and therefore also in $P_{k}^{+}$.

2. (ii) With shear thinning, the decrease in $P_{11}^{+}$ is reflected by a decrease in the amount of energy redistributed through $\varPi _{11}^{+}$ to budgets $B_{22}^{+}$ and $B_{33}^{+}$.

3. (iii) Such decrease in energy redistribution with decreasing flow index results in an increase of large-scale turbulence anisotropic behaviour (see figure 8a).

4. (iv) The lessening of production with shear thinning also leads to a decrease in the total dissipation rate $\epsilon _{11}^{+}$ and consequently in $\epsilon _{k}^{+}$.

5. (v) Since the mean dissipation rates ${M\epsilon }_{11}^{+}$ and ${M\epsilon }_{k}^{+}$ actually increase with decreasing flow index, the noted decrease in the corresponding total dissipation is attributed to the non-Newtonian contributions, in particular to ${\widetilde {T\epsilon }_{11}}^{+}$ and ${\widetilde {T\epsilon }_{k}}^{+}$, respectively.

6. (vi) Finally, note that, for all budgets, the terms associated with non-Newtonian rheology appear to be mainly important within the inner layer region.

3.5. Invariant analysis

All turbulent flows of practical interest are anisotropic, and drag-reducing fluids are known to have an even higher degree of anisotropy compared to Newtonian fluids (Escudier et al. 2009). The analysis of invariants corresponding to relevant tensors are typically performed to determine their degree of anisotropy and to identify realizable states of turbulence. To study the anisotropic behaviour of turbulence in GN fluids, for both large- and small-scale motions, an invariant analysis is presented. Any symmetric tensor $\sigma _{ik}$, such as the strain-rate tensor or the Reynolds stress tensor, may be decomposed into an isotropic and a traceless deviatoric part. The non-dimensional anisotropy tensor ${\mathsf {a}}_{ik}$ is then given by

(3.2)\begin{equation} {\mathsf{a}}_{ik} = \frac{\sigma_{ik}}{\sigma_{jj}}-\frac{1}{3} \delta_{ik} \end{equation}

and since the anisotropy tensor is traceless $({\mathsf {a}}_{ii}=0 )$, the relevant non-zero invariants are

(3.3a,b)\begin{equation} -II_a = \tfrac{1}{2} {\mathsf{a}}_{ik} {\mathsf{a}}_{ki}, \quad III_a = \tfrac{1}{3} {\mathsf{a}}_{ik} {\mathsf{a}}_{kj} {\mathsf{a}}_{ji}, \end{equation}

which are denoted as the second and third invariants of ${\mathsf {a}}_{ik}$, respectively.

Lumley & Newman (1977) introduced an anisotropy invariant map (AIM) consisting of a $- II$ versus $III$ plot to investigate the anisotropy of the Reynolds stress tensor. Such a map, nowadays called the Lumley triangle, was originally based on the invariants of ${\mathsf {b}}_{ik}$, the Reynolds stress anisotropy tensor, defined by

(3.4)\begin{equation} {\mathsf{b}}_{ik} = \frac{\,\overline{u_i^{\prime} u_k^{\prime}}\,}{\overline{u_j^{\prime} u_j^{\prime}}}-\frac{1}{3} \delta_{ik}. \end{equation}

Analogous AIMs have been used to study small-scale anisotropy by considering quantities such as the dissipation rate and the vorticity correlation tensor, see for instance Mansour, Kim & Moin (1988), Antonia et al. (1991) and Barri & Andersson (2010). Here, the anisotropy tensors corresponding to the vorticity correlation and the total dissipation rate are given by

(3.5)\begin{equation} {\mathsf{c}}_{ik} = \frac{\,\overline{\omega_i^{\prime} \omega_k^{\prime}}\,}{\overline{\omega_j^{\prime} \omega_j^{\prime}}}-\frac{1}{3} \delta_{ik} \end{equation}

and

(3.6)\begin{equation} {\mathsf{d}}_{ik} = \frac{\epsilon_{ik}}{\epsilon_{jj}}-\frac{1}{3} \delta_{ik}, \end{equation}

respectively. Additionally, an anisotropy tensor ${\mathsf {e}}_{ik}$, for the mean viscous dissipation only, may be given by

(3.7)\begin{equation} {\mathsf{e}}_{ik} = \frac{{M\epsilon}_{ik}}{{M\epsilon}_{jj}}-\frac{1}{3} \delta_{ik}. \end{equation}

In all the mentioned AIMs, realizable turbulent states are constrained within certain limiting values; see figure 17. At the origin, $-II_{a} = III_{a} =0$, all elements of a non-dimensional anisotropy tensor ${\mathsf {a}}_{ik}$ vanish and three-component isotropic turbulence (3C-IT) is found. If one diagonal component of ${\mathsf {a}}_{ik}$ is void of the corresponding property (energy or enstrophy, for example), e.g. ${\mathsf {a}}_{11} = -1/3$, and the two other diagonal components are equal, e.g. ${\mathsf {a}}_{22} = {\mathsf {a}}_{33} = 1/6$, then two-component isotropic turbulence (2C-IT) is encountered. Finally, one-component (1C) turbulence corresponds to a situation where all of the corresponding property is within a diagonal component, e.g. ${\mathsf {a}}_{11} = 2/3$.

Figure 17. Limiting states in anisotropy invariant map.

In the $-II_{a}$ versus $III_{a}$ plot, the origin is connected to the 1C and 2C-IT points through the relationship $III_a = \pm 2 (-II_a /3 )^{3/2}$ marking all cases of axisymmetric turbulence (i.e. two diagonal components of ${\mathsf {a}}_{ik}$ are equal and all off-diagonal components have vanished). The 2C-IT and 1C points, meanwhile, are connected by the line $-II_a = 3III_a +1/9$ where two-component (2C) states reside. With the aim of adding physical context to the limiting states in the Lumley triangle, consider, for example, axisymmetric deformation created by passing initially isotropic turbulence through a hypothetical slip-free axisymmetric nozzle or diffuser (see the same example as in Hanjalić & Launder (2011)). If the streamwise strain rate is $\partial \bar {u}/\partial x$ and the symmetry axis corresponds to the spanwise $z$ axis, then continuity and axisymmetry considerations lead to $-0.5\partial \bar {u}/\partial x = 0.5\partial \bar {v}/\partial y =\partial \bar {w}/\partial z$. In a similar manner, the stress field will remain axisymmetric through the contraction/expansion, i.e. ${\mathsf {b}}_{22}={\mathsf {b}}_{33}=-{\mathsf {b}}_{11}/2$ and thus the resulting invariants are ${II}_{b} =3{\mathsf {b}}_{11}^{2}/4$ and ${III}_{b} ={\mathsf {b}}_{11}^{3}/4$. Eliminating ${\mathsf {b}}_{11}$ then yields the previous relationship $III_b = \pm 2 (-II_b /3 )^{3/2}$; marking all cases of axisymmetric turbulence. The positive region (right-hand side in Lumley triangle) corresponds to axisymmetric expansion, which in the limit results in the 1C state, i.e. $\overline {v^{\prime } v^{\prime }}=\overline {w^{\prime } w^{\prime }}=0$. On the other hand, the negative region corresponds to axisymmetric contraction and in the extreme limit results in the 2C isotropic turbulence state, i.e. $\overline {u^{\prime } u^{\prime }}=0$.

The Lumley triangles for the different defined anisotropy tensors are presented in figure 18. In the same figure, the corresponding turbulence triangles (Choi & Lumley 2001) are shown as well. The map $\zeta ^{2} = -II/3$ versus $\xi ^{3} = III/2$ is used to emphasize the region in the proximity of the isotropic state. As noted by Emory & Iaccarino (2014), the Lumley triangle appears to provide more insight while studying states near the 1C and 2C limiting states whilst the turbulence triangle stretches the lower left quadrant of the Lumley triangle and focuses on the regions near the 2C-IT and 3C-IT limits.

Figure 18. AIMs: (a) Lumley triangle for ${\mathsf {b}}_{ik}$, (b) turbulence triangle for ${\mathsf {b}}_{ik}$, (c) Lumley triangle for ${\mathsf {c}}_{ik}$, (d) turbulence triangle for ${\mathsf {c}}_{ik}$, (e) Lumley triangle for ${\mathsf {d}}_{ik}$, (f) turbulence triangle for ${\mathsf {d}}_{ik}$, (g) Lumley triangle for ${\mathsf {e}}_{ik}$ and (h) turbulence triangle for ${\mathsf {e}}_{ik}$. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

For all GN fluids, the considered invariants vary from the 2C state limit in the vicinity of the wall to nearly isotropic at the channel's centre. Near the wall, the wall-normal diagonal component of the considered symmetric tensors are negligible, e.g. $\overline {v^{\prime } v^{\prime }} \approx \overline {\omega _2^{\prime } \omega _2^{\prime }} \approx 0$, and, as $y^{+}$ increases, the data transition from the upper boundary to the right-hand boundary of the Lumley triangle. With shear thinning, in comparison to the Newtonian case, a larger maximum value for $-II$ is observed on the anisotropy invariant maps corresponding to ${\mathsf {b}}_{ik}$, ${\mathsf {c}}_{ik}$ and ${\mathsf {d}}_{ik}$. In the AIMs, the movement of the data towards the 1C limit with shear thinning is attributed to the deficit in energy redistribution from budget $B_{11}^{+}$ to budgets $B_{22}^{+}$ and $B_{33}^{+}$, mentioned in § 3.4. In the case of the Reynolds stress tensor, for instance, near-wall behaviour (see figure 19) displays the expected increase in the streamwise component $\overline {u^{\prime } u^{\prime }}/u_{\tau }^{2}$ and the decrease of the wall-normal and spanwise components with shear thinning. Such increase of $\overline {u^{\prime } u^{\prime }}/u_{\tau }^{2}$, with decreasing flow index, is reflected in the Lumley triangle with near-wall 2C states closer to the 1C point. In a similar manner, the increase in $\overline {\omega _3^{\prime } \omega _3^{\prime }}$ and in the magnitude of ${M\epsilon }_{11}$ with shear thinning yield the observed movement of the data towards the 1C limit in the respective AIMs.

Figure 19. Near-wall behaviour for correlation of velocity fluctuations. Line styles ‘—–’, ‘- - -’, ‘$\cdot \cdot \cdot$’ and ‘- $\cdot$ -’ represent $(\overline {u^{\prime } u^{\prime }}^{1/2}/u_{\tau })/y^{+}$, $10(\overline {v^{\prime } v^{\prime }}^{1/2}/u_{\tau })/{y^{+}}^{2}$, $(\overline {w^{\prime } w^{\prime }}^{1/2}/u_{\tau })/y^{+}$ and $350(-\overline {u^{\prime } v^{\prime }}/u_{\tau }^{2})/{y^{+}}^{3}$, respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

Opposite to the mean dissipation, the data of the total dissipation in the vicinity of the wall do not approach states closer to the 1C limit with decreasing flow index. Owing to the non-Newtonian terms contributing to the total dissipation rate, there is a decrease in the magnitude of $\epsilon _{11}^{+}$ with shear thinning (as seen in § 3.4) causing the observed behaviour in the corresponding Lumley triangle.

For all GN fluids, after the location of the maximum in the different AIMs, the transition to anisotropic states occurs closer to the axisymmetric limit where $III>0$ takes place because one diagonal component is larger than the other two. Consider, for example, the streamwise component of the Reynolds stress in comparison with the other lateral components (see turbulence intensity profiles in § 3.1). Note that an actual axisymmetric state only occurs once the off-diagonal components of the corresponding anisotropy tensor are zero (e.g. ${\mathsf {b}}_{12} = \overline {u^{\prime } v^{\prime }}=0$ at the channel's centre).

As seen from the turbulence triangles, the data corresponding to ${\mathsf {b}}_{ik},\ {\mathsf {c}}_{ik},\ {\mathsf {d}}_{ik}$ and ${\mathsf {e}}_{ik}$ move closer to the axisymmetric limiting states with shear thinning. This observation is consistent with the decrease in the off-diagonal components corresponding to the different anisotropy tensors with shear thinning. See, for instance, the variation of the Reynolds shear stress with shear thinning in § 3.2. It is also observed that there is no appreciable difference in the turbulence triangles describing the mean viscous dissipation rate anisotropy and the total dissipation rate anisotropy in proximity to the axisymmetric limit. This behaviour is consistent with the observed decrease in the non-Newtonian terms contributing to $\epsilon _{11}^{+}$ in the outer layer region.

To complement the anisotropy analysis, the Lumley flatness (Lumley 1979) $\mathcal {F}=1+27III+9II$ and the axisymmetric parameter (Lee & Reynolds 1985) $\mathcal {A}=III/[ 2(-II/3)^{3/2}]$ are shown in figure 20. The Lumley flatness allows one to easily distinguish between 2C line states $(\mathcal {F}=0)$ and 3C-IT $(\mathcal {F}=1)$. Meanwhile, the axisymmetric parameter is a compact way to quantify axisymmetric modes: $\mathcal {A}=-1$ corresponds to axisymmetric states where $III<0$ and the correlation tensor (e.g. $\overline {u_i^{\prime } u_k^{\prime }}$ or $\overline {\omega _i^{\prime } \omega _k^{\prime }}$) has a smaller diagonal component than the other two equal ones; whilst $\mathcal {A}=1$ corresponds to axisymmetric states where $III>0$ and the correlation tensor has one larger diagonal component than the other two equal ones. Here terminology such as ‘rod-like’ or ‘cigar-shaped’ turbulence and ‘disk-like’ or ‘pancake-shaped’ turbulence is used while considering $\mathcal {A}=1$ and $\mathcal {A}=-1$ states of the vorticity correlation anisotropy tensor, respectively. Alternative terminology is omitted here to avoid misunderstanding; see, for instance, Simonsen & Krogstad (2005) where it is clarified that notation used for the axisymmetric states refers to the shape of the symmetrical second-order tensor with zero trace $\sigma _{ik}$ under study, e.g. the Reynolds stress tensor or the vorticity correlation tensor.

Figure 20. Lumley's flatness and axisymmetric parameter: (a) line styles ‘—–’, ‘- - -’, ‘$\cdot \cdot \cdot$’ and ‘- $\cdot$ -’ are used for $\mathcal {F}_b$, $\mathcal {F}_c$, $\mathcal {F}_d$ and $\mathcal {F}_e$; and (b) line styles ‘—–’, ‘- - -’, ‘$\cdot \cdot \cdot$’ and ‘- $\cdot$ -’ are used for $\mathcal {A}_b$, $\mathcal {A}_c$ (right axis), $\mathcal {A}_d$ and $\mathcal {A}_e$. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.

With the exception of the very near-wall behaviour in the total dissipation rate, the overall increase in anisotropy with decreasing flow index at the channel's centre and elsewhere is reflected by larger $-II$ values in the AIMs and consequently lower Lumley flatness. For all GN fluids, the profiles corresponding to the axisymmetric parameters show that, close to the end of the viscous sublayer, i.e. before the buffer layer region, the state $\mathcal {A}=1$ is approached. This observation is consistent with the noticed transition from the upper boundary to the right-hand boundary in the AIMs. After the maximum value of $-II$ is achieved on the respective maps, $\mathcal {A}_b$, $\mathcal {A}_d$ and $\mathcal {A}_e$ remain relatively close to the $\mathcal {A}=1$ state (especially as the flow index decreases) although there is a localized decrease in the axisymmetric parameters, most noticeable for $\mathcal {A}_b$, in the outer layer. This behaviour is observed on the respective turbulence triangles as a sudden but short deviation from the near-right boundary region towards the centre of the map, more clearly seen for increasing flow index.

With respect to the $\mathcal {A}_c$ scalar, an apparent change from a cigar-shaped axisymmetric state, at the end of the viscous sublayer, to a disk-like state $(III_c<0,\ \mathcal {A}_c=-1)$, within the buffer region, is noticed. The wall-normal position at which the pancake-shaped axisymmetric mode is observed appears to move slightly away from the wall with shear thinning. After such point, it is difficult to point out a trend in $\mathcal {A}_c$ since oscillations at different $y^{+}$ positions are noted. Nonetheless, states closer to disk-like turbulence appear to be more frequent, especially with decreasing flow index. Interestingly enough, close to the channel's centre, rod-like states are suddenly approached. This observation is consistent with the presented $\zeta$ versus $\xi$ map and the slight dominance of the spanwise vorticity autocorrelation $\overline {\omega _3^{\prime } \omega _3^{\prime }}$ over $\overline {\omega _1^{\prime } \omega _1^{\prime }}$ and $\overline {\omega _2^{\prime } \omega _2^{\prime }}$ near the channel's centre as seen in § 3.1.

4. Discussion and final remarks

Direct numerical simulations for statistically converged turbulent channel flow of GN fluids at a low Reynolds number have been performed. A GN fluid presents time-independent rheology and is free of plastic effects. In the simulation, GN fluid rheology has been incorporated through a relatively simple constitutive equation, the Carreau fluid model. To investigate the difference between Newtonian and shear-dependent fluid behaviour, the flow index $\alpha$ is varied. Here, when a trend is associated to shear thinning $\alpha <1$, the opposite trend is associated to shear thickening $\alpha >1$. Note that the selected rheological model is likely to be of no consequence since, for the same rheology characterization at large strain rates, quite similar statistics are expected even if a different rheological model is implemented (see e.g. Singh, Rudman & Blackburn 2016).

Through different statistics and analyses, it is found that shear-dependent fluid rheology seems to affect the channel flow mainly within the inner layer region. As we move further away from the viscous sublayer, the monotonic increase in the apparent fluid viscosity for $\alpha <1$ leads to drag reduction. For a constant driving pressure gradient, a decrease in the friction factor with decreasing flow index is reflected by an increase in the mean streamwise velocity for $y^{+} \gtrsim 10$ and in consequence higher mean bulk velocity $\int _{0}^{h^{+}}\bar {u}^{+}\,\mathrm {d}y^{+}/h^{+}$ and flow rate. The diagnosis function also reveals that the limited log region starts further away from the wall and with a slightly larger slope with shear-thinning fluid rheology. This thickening of the buffer layer is consistent with the observed increase in streamwise turbulence intensity and the shift in its peak value for $\alpha <1$.

The previous well-known drag-reducing-related characteristics are attributed to changes in the near-wall structures within the buffer layer region for shear-thinning fluid rheology. Quasi-streamwise vortices (rolls) are suppressed; as figure 8(b) shows, shear-thinning effects reduce the intensity of the streamwise vorticity fluctuations and also move the location of the local minimum and maximum values towards the channel's centre. These two extrema locations correspond to the average locations of the edge and centre, respectively, of the near-wall rolls (Moser & Moin 1984). Hence, the quasi-streamwise vortical structures not only decrease in intensity but also grow in size and depart from the wall.

The suppression of the near-wall streamwise rolls is accompanied by variations in the high- and low-speed fluid alternating regions (streaks) in the spanwise direction. It is recognized that there is an interaction between the rolls and the mean spanwise vorticity which induces low-speed streamwise streaks (see e.g. § 4.2.6 in Davidson 2015). These streaks eventually begin to oscillate and lift away from the wall during the so-called burst process. Each burst contains one or more ejections of low-speed fluid resulting from the same streak instability (Luchik & Tiederman 1987). Since burst or ejection events occur in a quasi-periodic manner, particular attention is paid to statistical quantities such as the average time between bursts and average spanwise spacing of the near-wall streaks. As seen in figure 2(c), for $\alpha <1$, the location of the minimum in the two-point correlation of the streamwise fluctuating velocity occurs at a larger spanwise separation, and a similar trend (not shown here) was found at other wall-normal positions within the buffer layer region. Note that the location of this minimum is the mean distance between a high-speed and a low-speed streak; thus the average streak spacing is twice such distance (Moser & Moin 1984). Since, for shear-thinning fluid rheology, a larger average streak spacing is found, a larger average time between bursts is expected and, in consequence, there is an inhibition of the related turbulence generating event. As revealed by the quadrant analysis, the decrease in Reynolds shear stress – leading to a reduction in the production of turbulent kinetic energy – arises due to a decrease from all positive production events but, within the buffer layer, especially due to a diminishing in contributions from ejection events.

Another common characteristic of drag-reducing flow, such as the decrease in the spanwise and wall-normal turbulence intensities, is explained from the Reynolds stress budgets. In connection to the decrease in turbulence production with shear thinning, there is a decrease in the velocity–pressure gradient terms, which are commonly split into a pressure–transport term and an energy redistributive pressure–strain rate term. Thus, since the velocity–pressure gradient terms play a dominant role in the redistribution of energy from streamwise to wall-normal and spanwise directions, their decrease leads to the observed trends for the turbulent intensities and consequently to an increase in anisotropy at the largest scales with $\alpha <1$. The budgets also show that, although the magnitude of the mean viscous dissipation ${M\epsilon }_{k}^{+}$ increases with shear thinning, the overall dissipation rate ${\epsilon }_{k}^{+}$ decreases in magnitude since there is less energy available for irreversible dissipation at the wall; hence the importance of the non-Newtonian terms, in particular the turbulent viscous dissipation rate related to the mean flow, whilst studying the total dissipation in GN fluids.

Turbulence anisotropy variations at both large and small scales due to the non-Newtonian rheology are studied through anisotropy invariant maps. The presented analysis reveals an overall increase in anisotropy for the Reynolds stress, vorticity correlation and mean viscous dissipation with shear thinning. In contrast to the mean viscous dissipation, in the near-wall region, the anisotropy of the total dissipation rate decreases with shear thinning due to non-Newtonian effects. On the other hand, regarding the data closer to the axisymmetric limits, as one moves from the wall towards the channel's centre, states closer to the axisymmetric limit where $III>0$ are seen for the Reynolds stress budgets and the dissipation rates. In the case of vorticity correlation, states between rod-like and disk-like turbulence are seen in the direction of the channel's centre. Here, states closer to disk-like or pancake-shaped turbulence appear to be preferred, especially with decreasing flow index. Finally, at the channel's centre, rod-like turbulence is approached. This observation is consistent with the slight dominance of the spanwise vorticity correlation over the other two components for all GN fluid cases. Anisotropy at the smallest scales is a recognized behaviour at low Reynolds numbers (Andersson, Zhao & Variano 2015) and appears to be even more noticeable with shear thinning.