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.

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

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 (1969Lumley ( , 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, τ ij,vis , is given by where μ = μ(γ ) is the apparent dynamic viscosity, which solely depends on the strain ratė γ = (2S ij S ji ) 1/2 , and 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.
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.

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.
where streamwise, wall-normal and spanwise directions are denoted by x = (x 1 , x 2 , x 3 ) = (x, y, z), the corresponding instantaneous velocity field by u = (u 1 , u 2 , u 3 ) = (u, v, w) and the pressure by p. Here t denotes time, ρ is the density of the incompressible, isothermal GN fluid, δ ij is the Kronecker delta and S ij = (∂u i /∂ x j + ∂u j /∂ 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 where μ ∞ and μ 0 are the 'infinite' and 'zero' shear-rate viscosities, respectively, Λ is a time constant and α 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 μ c , length L c , time L c /U c and stress ρU 2 c for the flow, (2.1)-(2.3) can be rewritten in non-dimensional form as where Re = ρU c L c /μ c is the Reynolds number and β = μ/μ 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 τ and L c = h in the outer layer and L c = μ c /(ρu τ ) in the inner layer; see for instance Pope (2000). Here u τ 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 μ w is taken as characteristic viscosity, i.e. μ c = μ 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 α = 1, in the aforementioned constitutive equation, allows us to recover Newtonian fluid behaviour.

Averaged governing equations
Introducing the Reynolds decomposition, i.e. splitting the variables into an ensemble average( ) and a fluctuating component ( ) , as u i =ū i + u i , p =p + p , β =β + β and S ij =S ij + S ij , into (2.4) and (2.5) and taking the average leads to

7)
Re where, in comparison with the Reynolds-averaged Navier-Stokes equations for a Newtonian fluid, a new non-Newtonian term (2β S ij ) 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. (2001Ptasinski et al. ( , 2003. Therefore, the total mean shear stressτ =τ 12 for a GN fluid is given bȳ 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 (·| w ). Thus, τ | w = 2(βS 12 + β S 12 )| w .

Reynolds stress budget equations
The transport equation for the correlation of the velocity fluctuations, u i u k , corresponding to a GN fluid is deduced in a similar manner as it is deduced for a Newtonian fluid. Thus, the added products The transport equation for the velocity correlation, see for instance Pinho (2003) or appendix A, in non-dimensional form is given by (2.10) 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: (i) =⇒ turbulent viscosity gradient transport rate, (2.17) (viii) =⇒ turbulent viscosity gradient transport rate related to mean flow, (2.20) (xi) =⇒ turbulent viscous diffusion rate related to mean flow, and (2.21) (xii) =⇒ turbulent viscous dissipation rate related to mean flow. (2.22) The overall equation is recast as where total dissipation ik and transport T ik rates, contributing to budget B ik = Re D(u i u k )/ Dt, are given by and 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 u i u k 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 μ c = μ w , L c = (μ w /ρ)/u τ and U c = u τ and are identified as { } + quantities.

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 τ = ρu τ h/μ 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(η), η 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, 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, x + and z + are the constant grid spacings in inner units corresponding to the streamwise and spanwise directions whilst y + min and 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. 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, β versuṡ γ + , is also shown in figure 1. For all GN fluid cases, a target Re τ = 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). 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 τ and the initial transients have been discarded. Statistically steady state is considered to have been achieved once a linear profile forτ is observed (Vinuesa et al. 2016). The initial transient time is over 200 convective time units, and 150 flow fields saved every ≈ 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.
where r is the separation vector between the two points, appear to be decreasing with increasing '---' and '· · ·' 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. 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.
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 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. 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η = (μ/ρ) 3/4 / 1/4 k , where k = ( 11 + 22 + 33 )/2 is the total mean dissipation rate (see Bradshaw & Perot (1993) and Bradshaw (1995) for a note about true dissipation), is analogous to η 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, y + min < η + in the near-wall region and y + max ≈ O(η + ) at the channel's centre. Similar trends are noted for the non-Newtonian fluid cases when comparing the wall-normal grid resolution withη + .
Based on the energy drop-off previously noted and on the η + andη + 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.

Low-order statistics
The mean (averaged in time and in homogeneous directions) streamwise velocity profilē u + and its diagnosis function y + dū + /dy + are presented in figure 6. In the same figure, 10 1 5 × 10 -1 10 1 FIGURE 4. Pre-multiplied one-dimensional spectral energy density: against normalized streamwise wavenumber k x h at (a) y + ≈ 5.07 and (b) y + ≈ 154.96, and against normalized spanwise wavenumber k z h at (c) y + ≈ 5.07 and (d) y + ≈ 154.96. Line styles '--', '---' and '· · ·' 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.
the mean viscosity profileμ + =β and an analogous quantity to the diagnosis but for the mean viscosity, y + dβ/dy + , are presented as well.
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 + 5); a buffer region (5 y + 55); a quite limited -if it exists at all -log region (55 y + 62); and a remaining outer region. Within the viscous sublayer, for all cases,ū + ≈ y + and for y + 10, minor deviation for the non-Newtonian cases is noted inū + andβ when compared to case N180. For y + 10,ū + increases with decreasing flow index α, 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 + dū + /dy + ≈ 1/κ, where κ is the von Kármán constant, increase with decreasing flow index.
Regardingβ, for y + 10, it deviates rapidly from the corresponding constant Newtonian value, as expected. Sinceβ appears to behave in a log manner within a certain y + range, an analogous diagnosis function based onβ has been considered. On a first impression, such a log region in the mean viscosity ( y + dβ/dy + ≈ Γ ) indeed occurs but different slopes Γ for cases P180 and D180 are observed; see figure 6(d).
For the outer layer region, the velocity defect profileū + ( y/h = 1) −ū + 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 + 62 or y/h 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).
, respectively, is presented in figure 8. The r.m.s. values corresponding to fluctuations in viscosity r.m.s.(β) = (μ ) 2 1/2 /μ 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. 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 + ≈ 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, u i u k and ω i ω k , 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β and r.m.s.(β) 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.(β) is larger for case P180, suggesting larger fluctuations fromβ with shear thinning. It is also noted that the peak in r.m.s.(β) 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 τ < 1000, a minimal -if anypresence 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 α < 1, there is an increase in the mean axial velocity, mean viscosity and axial turbulence intensity whereas both radial and azimuthal turbulence intensities decrease.

Mean shear stress budget
For a statistically stationary fully developed flow of an incompressible GN fluid, the total mean shear stressτ + (see § 2.1.1) is given bȳ whereτ + vis =βdū + /dy + ,τ + tur = −u v /u 2 τ andτ + GN = 2β S 12 + 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.τ + = 1 − y + /Re τ .
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τ + vis depends onβ and dū + /dy + , its increase can be discussed considering figure 6(b,c). For y + 10, there is minor variation in the mean viscosity and in the mean streamwise velocity gradient. Thus, the increase inτ + 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β, 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τ + GN depends on the sign of α − 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τ + tur , for all cases, appear to collapse in the outer region, suggesting that the cross-correlation u v /u 2 τ is independent of the rheology in that region. 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 . Mean shear stress budget. Line styles '--', '---', '· · ·' and '-· -' are used to identifyτ + ,τ + vis ,τ + tur andτ + GN , respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively. 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τ + 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 α < 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.

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 −u v are classified according to the sign of the velocity fluctuations into four categories or quadrants: Q1 (u > 0, v > 0), Q2 (u < 0, v > 0), Q3 (u < 0, v < 0) and Q4 (u > 0, v < 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 The subscript q refers to quadrant-conditioned averages. Line styles '---', '-· -', '· · ·' 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. flow index in the skewness and flatness profiles corresponding to the velocity fluctuations. The profiles are shown in figure 11.
From the non-fractional contributions to −u v /u 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 −u v /u 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 > 0 and v > 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, . Skewness S and flatness F profiles of velocity fluctuations: (a) Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.
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 −u v /u 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 + ≈ 30-35 where ejection events dominate due to large intermittent v > 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 > 0 and v < 0 signals. The line styles '--', '---' and '· · ·' are used to identify P + ik , T + ik and − + ik , respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively.
Finally, after the position y + ≈ 30-35, non-fractional contributions due to Q1, Q2 and Q3 events almost do not vary with rheology, and the observed decrease in −u v /u 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 < 0 and v > 0 signals with decreasing flow index, as reflected in the corresponding normalized skewness profiles.

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 + = u i u i /(2u 2 τ ) are presented in figure 13. For all cases, and as expected, the interaction between −u v /u 2 τ and mean shear causes a production rate in budget B + 11 . Since the mean velocity gradients increase with shear · · ·' are used to identify P + k , T + k and − + k , respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively. thinning (seeū + 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τ + tur in § 3.2). Note as well that, as commented in § 3.3, the peak in production at y + ≈ 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 + 11 . Everywhere P + 11 is balanced by T + 11 − + 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 + 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 + 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 − + 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 u i u k , i.e. B + 22 and B + 33 , are presented in appendix C.
For the transport rate T + 11 , the traditional terms Π + 11 , TT + 11 and MD + 11 contribute the most. The velocity-pressure gradient term Π + 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 + 10, such an increase is attributed to an increase of diffusion of u u /u 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 TD + 11 . Since, in the very near-wall region, the fluctuations in viscosity are rather small (see the r.m.s.(β) profile in § 3.1), the diffusion of u ū/u 2 τ is considered the reason for the large magnitude in 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 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 Tv + 11 . For the total dissipation rate + 11 , M + 11 and T + 11 appear to contribute the most. The mean viscosity increases with decreasing flow index; however, at the wall and in the  Mv + 11 and Tv + 11 . Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively. very near-wall region, the mean viscosity is almost unaffected by rheology. Therefore, the increase of magnitude in M + 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 + 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 T + 11 increases in magnitude with decreasing flow index and mainly affects + 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 T + 11 is attributed to the correlation between mean streamwise velocity gradients and streamwise velocity fluctuation gradients. Note that the decrease in the magnitude of + 11 at the wall with shear thinning is mainly attributed to dissipation at rate T + 11 , since the magnitude of M + 11 is increasing with decreasing flow index. The other non-Newtonian dissipation rate T + 11 also causes the decrease in the magnitude of the total dissipation rate + 11 with shear thinning. Nonetheless, compared to T + 11 , it is smaller in magnitude. The dissipation rate T + 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 Π + 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 + 12 is quite low and it is mainly due to the mean viscous 'dissipation' rate M + 12 . The rate T + 12 presents the opposite sign to M + 12 for the same  non-Newtonian fluid case, whilst the rate T + 12 either increases or decreases the total rate + 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 + 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 + 11 , since the terms in budget B + 11 contribute the most to the turbulent kinetic energy budget. The profile of rate Π + 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 Π + 11 is reduced with viscoelasticity, in particular close to the region of maximum production. The paper also reported a significant reduction in rate Π + 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 α < 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.
In summary, some of the most relevant findings are as follows.
(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 . (ii) With shear thinning, the decrease in P + 11 is reflected by a decrease in the amount of energy redistributed through Π + 11 to budgets B + 22 and B + 33 . (iii) Such decrease in energy redistribution with decreasing flow index results in an increase of large-scale turbulence anisotropic behaviour (see figure 8a). (iv) The lessening of production with shear thinning also leads to a decrease in the total dissipation rate + 11 and consequently in + k . (v) Since the mean dissipation rates M + 11 and M + 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 T + 11 and T + k , respectively. (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 largeand small-scale motions, an invariant analysis is presented. Any symmetric tensor σ 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 a ik is then given by and since the anisotropy tensor is traceless (a ii = 0), the relevant non-zero invariants are − II a = 1 2 a ik a ki , III a = 1 3 a ik a kj a ji , (3.3a,b) which are denoted as the second and third invariants of 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 b ik , the Reynolds stress anisotropy tensor, defined by 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 respectively. Additionally, an anisotropy tensor e ik , for the mean viscous dissipation only, may be given by 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 a ik vanish and three-component isotropic turbulence (3C-IT) is found. If one diagonal component of a ik is void of the corresponding property (energy or enstrophy, for example), e.g. a 11 = −1/3, and the two other diagonal components are equal, e.g. a 22 = 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. a 11 = 2/3.
In the −II a versus III a plot, the origin is connected to the 1C and 2C-IT points through the relationship III a = ±2(−II a /3) 3/2 marking all cases of axisymmetric turbulence (i.e. two diagonal components of 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 ∂ū/∂ x and the symmetry axis corresponds to the spanwise z axis, then continuity and axisymmetry considerations lead to −0.5∂ū/∂ x = 0.5∂v/∂ y = ∂w/∂z. In a similar manner, the stress field will remain axisymmetric through the contraction/expansion, i.e. b 22 = b 33 = −b 11 /2 and thus the resulting invariants are II b = 3b 2 11 /4 and III b = b 3 11 /4. Eliminating b 11 then yields the previous relationship III b = ±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. v v = w w = 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. u u = 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 ζ 2 = −II/3 versus ξ 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.
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. v v ≈ ω 2 ω 2 ≈ 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 b ik , c ik and 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 u u /u 2 τ and the decrease of the wall-normal and spanwise components with shear thinning. with near-wall 2C states closer to the 1C point. In a similar manner, the increase in ω 3 ω 3 and in the magnitude of M 11 with shear thinning yield the observed movement of the data towards the 1C limit in the respective AIMs.
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 . Near-wall behaviour for correlation of velocity fluctuations. Line styles '--', '---', '· · ·' and '-· -' represent (u u 1/2 /u τ )/y + , 10(v v 1/2 /u τ )/y + 2 , (w w 1/2 /u τ )/y + and 350(−u v /u 2 τ )/y + 3 , respectively. Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively. the magnitude of + 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. b 12 = u v = 0 at the channel's centre).
As seen from the turbulence triangles, the data corresponding to b ik , c ik , d ik and 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 + 11 in the outer layer region. To complement the anisotropy analysis, the Lumley flatness (Lumley 1979) F = 1 + 27III + 9II and the axisymmetric parameter (Lee & Reynolds 1985) A = III/[2(−II/3) 3/2 ] are shown in figure 20. The Lumley flatness allows one to easily distinguish between 2C line states (F = 0) and 3C-IT (F = 1). Meanwhile, the axisymmetric parameter is a compact way to quantify axisymmetric modes: A = −1 corresponds to axisymmetric states where III < 0 and the correlation tensor (e.g. u i u k or ω i ω k ) has a smaller diagonal component than the other two equal ones; whilst 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 A = 1 and 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 σ ik under study, e.g. the Reynolds stress tensor or the vorticity correlation tensor.
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 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, A b , A d and A e remain relatively close to the A = 1 state (especially as the flow index decreases) although there is a localized decrease in the axisymmetric parameters, most noticeable for 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 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, 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 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 ζ versus ξ map and the slight dominance of the spanwise vorticity autocorrelation ω 3 ω 3 over ω 1 ω 1 and ω 2 ω 2 near the channel's centre as seen in § 3.1.

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 α is varied. Here, when a trend is associated to shear thinning α < 1, the opposite trend is associated to shear thickening α > 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 α < 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 + 10 and in consequence higher mean bulk velocity h + 0ū + dy + /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 α < 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 α < 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 α < 1. The budgets also show that, although the magnitude of the mean viscous dissipation M + k increases with shear thinning, the overall dissipation rate + 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.

Appendix B. Verification
The results corresponding to case N180 shown in figure 21 agree with those reported by Kim et al. (1987) and Moser et al. (1999). Contributions to the turbulent kinetic energy budget in figures 13(a) and 16(a-c) and the η + profile presented in figure 5 also appear to agree qualitatively with the profiles reported by Mansour et al. (1988) and Antonia et al. (1991), respectively.
Appendix C. Contributions to budgets B + 22 and B +

33
Contributions from the different terms in (2.23) to total transport and dissipation rates for budgets B + 22 and B + 33 are shown in figures 22 and 23. For transport rate T + 22 , the velocity-pressure gradient, the turbulent transport and the mean viscous diffusion rates contribute the most. The rate Π + 22 , which is the source of energy for budget B + 22 , decreases with decreasing flow index as expected. Similar to the production rate P + 11 in budget B + 11 , the . Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively. energy in Π + 22 is transported towards and away from the wall, mainly through the turbulent transport rate. Since there is less energy to be distributed with shear thinning, within the inner region, the magnitudes of TT + 22 and MD + 22 decrease with shear thinning as well. With respect to the non-Newtonian terms, the transport rate TD + 22 decreases with decreasing flow index whilst Mv + 22 and Tv + 22 increase. However, these terms are rather small and do not contribute significantly to the total transport rate.
The total dissipation rate + 22 balancing T + 22 is mainly due to the mean viscous dissipation rate M + 22 , which presents the same trends of Π + 22 with rheology, i.e. it decreases in magnitude with decreasing flow index. The other non-zero contribution to rate + 22 , although fairly low in magnitude, is T + 22 . This non-Newtonian dissipation rate increases in magnitude with shear thinning.
For transport rate T + 33 , mainly the pressure-velocity gradient term and the mean viscous diffusion (especially very close to the wall) are important. As for budget B + 22 , less energy is redistributed from budget B + 11 with decreasing flow index. The magnitudes of Π + 33 and consequently of TT + 33 and MD + 33 decrease with shear thinning. Other non-zero transport terms, associated with the non-Newtonian rheology, do not contribute much to rate T + 33 . As for T + 22 , the mean viscous dissipation rate contributes the most to the total rate T + 33 and decreases in magnitude with shear thinning since there is less energy that may be FIGURE 23. Contributions to T + 33 and − + 33 in budget B + 33 : (a) line styles '--' and '---' are used for Π + 33 and TT + 33 ; (b) line styles '--' and '· · ·' are used for MD + 33 and TD + 33 ; (c) line styles '--' and '· · ·' are used for −M + 33 and −T + 33 ; and (d) line styles '---' and '· · ·' are used for Mv + 33 and Tv + 33 . Profiles corresponding to P180, N180 and D180 are identified by red, black and cyan colours, respectively. dissipated. The other non-zero dissipation rate T + 33 is small in comparison to M + 33 and shows the opposite behaviour with decreasing flow index.