1. Introduction
Turbulence is a state of erratic, chaotic and unpredictable fluid motion at very high Reynolds numbers – the ratio of the fluid’s inertial forces to its viscous forces (Frisch Reference Frisch1995; Pope Reference Pope2001). Statistically stationary, homogeneous and isotropic turbulent flows represent the most basic example of the fundamental problem of turbulence, where the statistical characteristics of velocity fluctuations can be probed (Ishihara, Gotoh & Kaneda Reference Ishihara, Gotoh and Kaneda2009). Typically, energy is injected into the flow at large length scales, whereas viscous processes play a significant role in dissipating energy from the flow at the small scales (Alexakis & Biferale Reference Alexakis and Biferale2018). Scaling exponents are universal in the intermediate range of scales, meaning that they are independent of the turbulence generation process, with the energy spectra decaying as
$k^{-5/3}$
, with
$k$
being the wavenumber, as predicted by Kolmogorov’s dimensional arguments (Kolmogorov Reference Kolmogorov1941).
Both in nature and in industry, turbulent flows are frequently multiphase, meaning that they are packed with particles, they may consist of fluid mixes or contain additives like polymers, and complex dynamics can arise, such as fluid elasticity (Benzi & Ching Reference Benzi and Ching2018), yield stress (Balmforth, Frigaard & Ovarlez Reference Balmforth, Frigaard and Ovarlez2014) or shear-dependent viscosity (Larson & Desai Reference Larson and Desai2015). In recent years, a lot of research has been done on fluids that include long-chain polymers, since adding high molecular weight polymers to a turbulent pipe flow at small concentrations can significantly lower the drag (White & Mungal Reference White and Mungal2008). Fluid elasticity can also lead to flow instabilities (Shaqfeh Reference Shaqfeh1996) and turbulence (Steinberg Reference Steinberg2021), even at very small Reynolds numbers. The disordered fluid motion that results from elastic instabilities at low or vanishing Reynolds numbers exhibits characteristics resembling those of classical Newtonian turbulence, such as a power-law scaling in the velocity spectrum, and intermittency (Datta et al. Reference Datta2022; Singh et al. Reference Singh, Perlekar, Mitra and Rosti2024). At all Reynolds numbers, however, the energy cascade from large to small scales is significantly altered by the presence of polymers, with the turbulent kinetic energy spectrum exhibiting a multiscale behaviour (Valente, Da Silva & Pinho Reference Valente, Da Silva and Pinho2014; Rosti, Perlekar & Mitra Reference Rosti, Perlekar and Mitra2023), and the cascade process being altered not only quantitatively but also qualitatively, thus strongly departing from Kolmogorov theory.
A less studied effect of polymers is to provide a shear-dependent viscosity, i.e. with a nonlinear relation between the shear stress and the shear rate, which is often difficult to decouple from the fluid elasticity. When the shear stress increases more than linearly with the shear rate, the fluid is called dilatant or shear-thickening, whereas in the case of the opposite behaviour, i.e. when the shear stress increases less than linearly with the shear rate, the fluid is called pseudo-plastic or shear-thinning. In this work, we focus on fluids whose rheological response does not depend explicitly on time, but depends only on the present shear rate, often called generalised Newtonian fluids (Bird Reference Bird1976). Recent results (Rosti et al. Reference Rosti, Perlekar and Mitra2023; Soligo & Rosti Reference Soligo and Rosti2023; Amor et al. Reference Amor, Soligo, Mazzino and Rosti2024) have shown a decay of the energy spectrum of such fluids similar to that predicted by Kolmogorov for Newtonian fluids, but a clear understanding is still missing, as well as a theory to explain this behaviour, which are the focus of the present work. Thus after introducing the governing equation, the fluid model and the set-up considered in § 2, we analyse the results obtained by direct numerical simulations in § 3, and summarise our findings in § 4.
2. Methodology
We consider an incompressible three-dimensional turbulent flow field, governed by the generalised Navier–Stokes equations. In an inertial, Cartesian frame of reference, the equations of momentum and mass conservation for the incompressible flow read as
\begin{equation} \rho \left ( \frac {\partial u_i}{\partial t} + \frac {\partial u_i u_j}{\partial x_j} \right ) = -\frac {\partial p}{\partial x_i} + \frac {\partial 2 \unicode{x03BC} \mathcal{S}_{ij}}{\partial x_j} + \rho f_i^T \quad\textrm {and}\quad \frac {\partial u_i}{\partial x_i} = 0, \end{equation}
where
$u_i$
is the fluid velocity field,
$p$
is the pressure,
$\mathcal{S}_{ij}= ({1}/{2}) ( {\partial u_i}/{\partial x_j} + {\partial u_j}/{\partial x_i} )$
is the strain rate tensor,
$f_i^T$
is the force used to sustain the turbulent flow, and
$\rho$
and
$\mu$
are the density and dynamic viscosity of the fluid (being their ratio
$\nu =\mu /\rho$
, the kinematic viscosity). Note that we are using the Einstein sum rule for repeated indices. The forcing is used to generate and sustain a fully turbulent flow with quasi-homogeneous, isotropic and stationary statistics; in particular, turbulence is sustained using the Arnold–Beltrami–Childress (ABC) cellular-flow forcing, which is a combination of sinusoids with wavelength equal to the domain size and amplitude
$F^T$
, kept the same for all cases. The equations of motion are solved in a triperiodic cubic box of size
$L$
, with periodic boundary conditions applied in all three Cartesian directions, discretised with
$1024$
grid points in each direction to ensure that all the scales down to the smallest dissipative ones are properly solved, i.e.
$\eta \lesssim \Delta x$
, where
$\Delta x$
is the grid spacing, and
$\eta$
is the Kolmogorov scale of the Newtonian simulation. While this criterion may not be directly applicable to the shear-thinning fluids, we ensure the correctness of the results by using a grid that is suitable for a Reynolds number larger than the one reached here (see e.g. Olivieri, Cannon & Rosti (Reference Olivieri, Cannon and Rosti2022) and Cannon, Olivieri & Rosti (Reference Cannon, Olivieri and Rosti2024), who reached
$Re_\lambda \approx 433$
).
In the current study, the fluid is non-Newtonian, and we focus on the simple inelastic power-law models, where the local viscosity of the fluid is a function of the sole local value of shear rate, i.e.
$\mu =K \dot {\gamma }^{n-1}$
, where
$n$
is the flow power index, and
$K$
is the fluid consistency index. A Newtonian behaviour is recovered when
$n=1$
, while values of the flow index above and below unity,
$n\gt 1$
and
$n\lt 1$
, denote shear-thickening and shear-thinning fluids, respectively. The consistency index
$K$
measures how strongly the fluid responds to the imposed deformation rate. In the previous relation, the local shear rate
$\dot {\gamma }$
is the second invariant of the strain-rate tensor
$\mathcal{S}_{ij}$
and is computed by its dyadic product, i.e.
$\dot {\gamma }=\sqrt {2 \mathcal{S}_{ij} \mathcal{S}_{ij}}$
. The viscosity of a power-law shear-thinning fluid becomes infinite for null shear rate; to overcome this, the Carreau fluid model is used instead, in which the local viscosity is computed as
\begin{equation} \mu =\mu _\infty + \left ( \mu _0 - \mu _\infty \right ) \left [ 1 + \left ( \mathcal{K} \dot {\gamma } \right )^2 \right ]^{(n-1)/2}. \end{equation}
In this equation,
$\mu _0$
and
$\mu _\infty$
indicate the lower and upper limits of fluid viscosity at zero and infinite shear rates, and the flow index
$n$
and time constant
$\mathcal{K}$
have similar interpretations to the material properties of the power-law fluid.
We use the in-house Fortran code Fujin to simulate the flow. The fluid equations are solved numerically on a staggered grid using a second-order finite difference code, with pressure points located at the cell centre, and velocity components at the cell faces. Equations (2.1) are advanced in time by a fully explicit fractional step method, based on the second-order Adams–Bashforth method; finally, the Poisson equation is solved by fast Fourier transform.
2.1. Set-up
We consider various fluids with shear dependent viscosity, and report their rheology in figure 1(
$a$
). In particular, we consider a Newtonian fluid with uniform viscosity
$\mu _0$
, and two shear-thinning fluids with shear dependent viscosity
$\mu$
defined by (2.2). These have a fixed ratio of zero to infinity viscosity
$\mu _0/\mu _\infty$
equal to
$10^4$
, and varying power indices
$n$
equal to
$0.7$
and
$0.4$
. Finally, the consistency index
$\mathcal{K}$
is chosen to ensure a low probability of being in the plateaus of viscosity at
$\dot {\gamma } \to 0$
and
$\infty$
. All the fluids are tested in a number of simulations with different values of Reynolds number
$Re_{ABC}$
, defined based on the amplitude of the forcing and on the zero-shear viscosity
$\mu _0$
. The parameters are chosen to achieve in the purely Newtonian case a Taylor Reynolds number
$Re_\lambda = \rho u_{\textit {rms}} \lambda /\langle \mu \rangle \approx 300$
,
$200$
and
$100$
, where
$u_{\textit {rms}} = \sqrt {\langle u_i u_i \rangle /3}$
is the root mean square of the fluid velocity, and
$\lambda =u_{\textit {rms}} \sqrt {15 \langle \mu \rangle / ( \rho \langle \varepsilon \rangle ) }$
is the Taylor length scale, with
$\langle \varepsilon \rangle = \langle 2 \nu \mathcal{S}_{ij} \mathcal{S}_{ij} \rangle$
the mean turbulent dissipation rate, with
$\langle \cdot \rangle$
indicating average over time, space and ensemble. In order to show that the results are general and independent of the values of
$\mu _0/\mu _\infty$
and
$\mathcal{K}$
, we perform three additional simulations, one where we set
$\mu _\infty =0$
, and another two where we increase
$\mathcal{K}$
. All the simulation parameters are summarised in table 1. We allow the flow to reach a statistically steady state, and measure the statistics presented in the rest of the paper, averaged in time over
$40$
snapshots collected over approximately
$15$
large-eddy turnover times
$L/u_{rms}$
.
Figure 1. (
$a$
) Fluid viscosity as a function of the shear rate. (
$b$
) Probability density function (PDF) of the fluid viscosity (main) and shear rate (inset) in the turbulent flows. The dark grey, blue and brown colours distinguish the different power indices
$n=1$
,
$0.7$
and
$0.4$
, respectively. The solid, dashed and dotted lines distinguish the different Reynolds numbers
$Re_\lambda \approx 300$
,
$200$
and
$100$
. The circles show the case with
$\mu _\infty =0$
, while the cyan and gold triangles indicate the case with different
$\mathcal{K}$
and
$n=0.7$
and
$0.4$
.
Table 1. Main input and output parameters of the cases investigated in the present study, together with the colours and line styles used throughout the work.
3. Results
We start the analysis of how turbulence is affected by the shear-dependent viscosity by looking at the changes of the Taylor Reynolds number in table 1. Interestingly, we notice that the Taylor Reynolds number does not change significantly for the different fluids, provided that the averaged viscosity is used in the definitions of the Reynolds number and in the Taylor microscale, while the local viscosity is used when averaging the turbulent dissipation rate. This result is significant, hinting towards adapting Kolmogorov theory to the generalised Newtonian fluids, as will be discussed later on. First, however, we show that the viscosity is indeed changing significantly throughout the domain by looking at its probability distribution function in figure 1(
$b$
), with table 1 reporting the first three moments of the distributions. At a fixed Reynolds number, the viscosity spans a monotonically increasing range of values as
$n$
is reduced, while the distribution shrinks when the Reynolds number is reduced at constant
$n$
. When
$n \to 1$
, the mean viscosity
$\langle \mu \rangle$
tends to the classical Newtonian viscosity. Also, the distributions never reach the values of the infinite viscosity
$\mu _\infty$
, and very rarely approach the zero one
$\mu _0$
; we can thus conclude that what is observed hereafter is not influenced by the values of
$\mu _0$
and
$\mu _\infty$
, which is true for large
$\mu _0/\mu _\infty$
. It is thus representative of the broader family of power-law fluids. In the opposite limit when
$\mu _0/\mu _\infty$
is tending to unity, the change of viscosity is rather limited, whatever the value of
$n$
, thus the Newtonian results are expected to hold again.
Figure 2. Instantaneous visualisation of the vorticity magnitude on a plane in the middle of the cubic box for (a) Newtonian fluids with
$n=1$
, and (b) shear-thinning fluids with
$n=0.4$
. Colours go from white to blue from zero to maximum vorticity.
Figure 3. Energy spectra (
$a$
) for different power indices
$n$
at
$Re_\lambda \approx 300$
(different colours), and for different Reynolds numbers
$Re_\lambda$
with
$n=0.4$
(different line styles). The latter have been shifted downwards for clarity. (
$b$
) Energy spectra for different Reynolds numbers and power indices, normalised with the Kolmogorov scale
$\eta$
and turbulent dissipation rate
$\varepsilon$
.
The variable viscosity profoundly alters the flow, especially the velocity gradients; see e.g. figure 2. To understand how the variable viscosity is affecting the flow at all scales, we start by looking at the energy spectra
$E(k)$
, reported in figure 3(
$a$
) for different power indices
$n$
and different Reynolds numbers. For every power index
$n$
, the spectra exhibit an inertial range where the energy decays according to the classical Kolmogorov scaling
$k^{-5/3}$
, followed by a sharp decrease in energy within the dissipation range. As
$n$
is reduced, we observe an extension of the inertial range of scales, and a shallower decay of energy at the smaller scales. Instead, as
$Re_\lambda$
is reduced, the typical reduction of the inertial range and expansion of the dissipation range are observed. Thus the variable viscosity keeps unaltered the large scales of the flow, and brings higher level of energy at the small scales, with an extension of the inertial range and a weaker energy decay in the dissipation range.
While finding the Kolmogorov scaling in the inertial range of scales of a variable viscosity fluid may seem strange, its appearance can be understood by rethinking the Kolmogorov hypotheses in our framework of variable viscosity. The hypothesis of local isotropy at the smallest scale remains unchanged, while the first and second similarity hypotheses need to be revisited. Indeed, both involve definitions of viscosity and dissipation rate, which are then used to define the Kolmogorov scales, which are uniquely defined in classical turbulence while not so in our case of variable viscosity. As done before when computing the Taylor Reynolds number in table 1, here we take the mean viscosity
$\langle \mu \rangle$
and the mean turbulent dissipation rate defined as
$\langle \varepsilon \rangle = \langle 2 \nu \mathcal{S}_{ij} \mathcal{S}_{ij} \rangle$
. With these two quantities, we can then extend Kolmogorov’s theory in a straightforward manner, leading to the definition of the Kolmogorov scale
$\eta = ( \langle \mu \rangle ^3/\langle \varepsilon \rangle )^{1/4}$
, which allows us to separate the small scales of motion
$r \ll \eta$
, where the statistics are supposed to be universal functions of
$\langle \mu \rangle$
and
$\langle \varepsilon \rangle$
(first similarity hypothesis), from the intermediate ones
$r \gg \eta$
, where the statistics are independent of
$\langle \mu \rangle$
and uniquely determined by
$\langle \varepsilon \rangle$
(second similarity hypothesis). By simple dimensional argument, we can then rescale the energy spectra as done in figure 3(
$b$
), where we confirm that all the cases with different
$n$
and
$Re_\lambda$
indeed collapse onto a master curve in the inertial range of scales. However, this is not the case in the dissipative region, where the spectra remain segregated by the different
$n$
, while collapsing for the different
$Re_\lambda$
. This is because in our case of variable viscosity, the smallest scales should be expected to be universal functions of
$\langle \mu \rangle$
,
$\langle \varepsilon \rangle$
and
$n$
, where the latter is an additional non-dimensional parameter; the same does not apply for the intermediate range of scales, where the effects of viscosity are not relevant. Thus we can expect the Kolmogorov theory to hold true at the large and intermediate scales for
$r\gg \eta$
, while at the smallest scale, a dependence on
$n$
can remain. In the following, we will thus verify this, and study the dependence of the small scales on
$n$
.
Figure 4. (
$a$
) Energy balance for different power indices
$n$
at
$Re_\lambda \approx 300$
. For
$n=0.4$
, we also show the dissipation term decomposed into mean (upward triangles) and fluctuating (downward triangles) parts. (
$b$
) Energy balance for different Reynolds number and power index, normalised with the Kolmogorv scale
$\eta$
. The ordinates in both (a) and (b) are divided by
$\langle \varepsilon \rangle$
.
To gain a more detailed insight, we look at the scale-by-scale energy transfer balance
\begin{equation} \mathcal{P} \!\left ( k \right ) + \Pi \!\left ( k \right ) + \mathcal{D} \!\left ( k \right ) = \langle \varepsilon \rangle , \end{equation}
where
$\mathcal{P} ( k )$
is the production term associated with the external forcing,
$\Pi ( k )$
is the energy flux associated with the nonlinear convective term, and
$\mathcal{D} ( k )$
is the viscous dissipation, defined as
$\mathcal{P} ( \kappa ) = \int _\kappa ^\infty ( \hat {f} \, \hat {u}^\star + \hat {f}^\star \, \hat {u} )/2\, {\rm d}\kappa$
,
$\Pi ( \kappa ) = \int _\kappa ^\infty - ( \hat {G} \, \hat {u}^\star + \hat {G}^\star \, \hat {u} )/2\, {\rm d}\kappa$
and
$\mathcal{D} ( \kappa ) = \int _0^\kappa ( \hat {D} \, \hat {u}^\star + \hat {D}^\star \, \hat {u} )\, {\rm d}\kappa$
. Here, a hat denotes the Fourier transform operator, a superscript
$^\star$
denotes complex conjugate, and
$\hat {G}$
and
$\hat {D}$
are the Fourier transforms of the nonlinear and viscous terms (Pope Reference Pope2001). Note that the production term and the flux are obtained by integrating from
$\kappa$
to
$\infty$
, while the dissipation term is integrated from
$0$
to
$\kappa$
to obtain a positive quantity that matches
$\langle \varepsilon \rangle$
; also, due to the variable viscosity, the usual relation
$\hat {D} ( \kappa ) = 2 \nu k^2\, E ( \kappa )$
no longer holds.
In figure 4(
$a$
), we show the energy fluxes and dissipation for different power indices
$n$
at large Reynolds number. The energy is injected into the system by the forcing trough
$\mathcal{P} ( k )$
, which is not null only at the largest scale
$k=k_{\textit{max}}$
(not shown); energy is then carried from large to small scales by the advective term
$\Pi (k )$
, where it is ultimately dissipated by the viscous term
$\mathcal{D} ( k )$
. The picture remains qualitatively the same when
$n$
is reduced, except for a stronger advective term
$\Pi (k )$
that dominates over a larger range of scales, confining the dissipation to even smaller scales. Here,
$\mathcal{D} (k )$
is defined with the local viscosity; we can decompose this term into a term coming from the mean viscosity
$\langle \mu \rangle$
(upward triangles), which is formally the same as in classical turbulence, and the rest which arises from its fluctuation
$\mu - \langle \mu \rangle$
(downward triangles). The decomposition suggests that the fluctuations of viscosity provide a negative dissipation, or in other words, are a production mechanism of turbulent kinetic energy acting at the small scales, which is balanced on average by the excess of dissipation coming from the mean viscosity. Finally, figure 4(
$b$
) shows again the energy balance for all
$n$
and
$Re_\lambda$
, normalised with the Kolmogorov scaling. An excellent collapse of the curves is found.
Figure 5. (
$a$
) Structure function
$S_p ( r )$
for different power indices
$n$
at
$Re_\lambda \approx 320$
. The structure functions of orders
$2$
,
$4$
and
$6$
are shifted vertically for visual clarity. (
$b$
) Structure function
$S_p ( r )$
compensated with the expected scaling for
$r \rightarrow 0$
. (
$c$
) Compensated third-order structure function
$S_3 ( r )$
. (
$d$
) Extended self-similarity of structure functions, with the prediction based on the refined similarity hypothesis. The inset shows the multifractal spectra of the dissipation rate, with the symbols representing experimental data for a Newtonian fluid, taken from Meneveau & Sreenivasan (Reference Meneveau and Sreenivasan1991).
The changes observed in the energy spectra can be appreciated also when looking at the structure functions in the real space; the longitudinal structure functions are defined as
$S_p ( r ) = \langle [ \delta u ( r ) ]^p \rangle$
, where
$p$
is the order of the moment, and
$\delta u ( r ) = u ( x+r ) - u ( x )$
is the velocity increment across a length
$r$
. For the single-phase case,
$S_p \sim r^p$
at small scales, and
$S_p \sim r^{p/3}$
in the inertial range, as predicted by the Kolmogorov theory (Kolmogorov Reference Kolmogorov1941), with some deviations due to the flow intermittency (Ishihara et al. Reference Ishihara, Gotoh and Kaneda2009), caused by extreme events that are localised in space and time that break the Kolmogorov similarity hypothesis. These extreme events correspond to the large tails in the velocity increment distributions, and make significant contributions to the high-order moments.
Figure 5(
$a$
) shows the structure functions of orders
$2$
,
$4$
and
$6$
versus the separation
$r$
, for different values of the power index
$n$
at large Reynolds number. We observe that the structure functions collapse well for different
$n$
in the inertial range of scale, while they fan out at smaller scales in the dissipation range, consistently with the observation from the spectra. Assuming that the velocity is a differentiable function, we can expect
$S_p ( r ) \sim r^p$
for
$r\rightarrow 0$
, and following Schumacher, Sreenivasan & Yakhot (Reference Schumacher, Sreenivasan and Yakhot2007), we plot in figure 5(
$b$
) the same structure functions, compensated with their expected scaling for the viscous range
$r^p$
, such that the function must approach a constant for reducing
$r$
; this is evident for all the Newtonian structure functions, while for the shear-thinning fluid it is clear only for
$S_2$
and less for the other orders, indicating that even smaller
$r$
are needed to achieve the predicted scalings compared to a Newtonian fluid. Next, we assess the validity of the most important relationship valid in classical turbulence, the celebrated
$4/5$
law by Kolmogorov:
$S_3 = -{4}/{5} \langle \varepsilon \rangle r$
. Verifying its validity is the most direct proof of the possibility of extending Kolmogorov theory to shear-thinning fluids and on the correctness of using
$\langle \varepsilon \rangle$
. The results of our simulations are shown in figure 5(
$c$
), where the third-order structure function has been compensated with the expected value. The expected power law and the proper negative sign, which is connected to a direct energy cascade, are both visible, with an extension of the range where the scaling holds with
$n$
. Finally, we assess the flow intermittency using the extended self-similarity form (Benzi et al. Reference Benzi, Ciliberto, Tripiccione, Baudet, Massaioli and Succi1993), by plotting
$S_6$
against
$S_2$
, as shown in figure 5(
$d$
). In the limit case, where extreme events do not occur, the
$S_6 \sim S_2^{6/2}$
power law holds, while the deviation from this behaviour is a measure of the flow intermittency. It should be noted that there are already deviations from the theoretical prediction in the single-phase case, with the correction suggested by Kolmogorov (Reference Kolmogorov1962) offering a good prediction of the data. Moving to the effect of the variable viscosity, figure 5(
$d$
) shows that a reduction of the power-law index generally leads to a larger deviation from the expected scaling at the small scales; the deviation from the single phase is monotonic with
$n$
, and starts at larger
$r$
as
$n$
decreases, leading to a stronger flow intermittency.
Examining the intermittency of the dissipation in space is an alternative method of examining the impact of extreme events. Using the procedure outlined by Frisch (Reference Frisch1995), we average the dissipation within a spherical region of radius
$l$
to get
$\varepsilon _l$
. We pick a range of moments
$-6 \le q \le 6$
,
$\varepsilon _l^p$
, and average them over various locations and times. At large
$l$
,
$\varepsilon _l$
is by definition equivalent to the bulk value
$\langle \varepsilon \rangle$
, while the two differ for small
$l$
, due to the localised regions of high dissipation in the fluid. We assume that
$\langle \varepsilon _l^q \rangle \sim l^{\tau _q}$
and then compute
$\tau _q$
through fitting; finally, we obtain the multifractal spectra
$\mathcal{F} ( \alpha )$
using a Legendre transformation:
$\alpha = {\rm d}\tau / {\rm d}q +1$
and
$\mathcal{F} = q ( \alpha -1 ) - \tau _q +3$
. The inset of figure 5(
$d$
) shows the multifractal spectra for all cases at large Reynolds number, and we observe that all the multifractal spectra have peaks at
$\alpha \approx 1$
and
$\mathcal{F} \approx 3$
, showing a background of space-filling dissipation; however, the presence of tails in the spectra indicates that the dissipation fields are not self-similar. No significant difference is found for the various
$n$
, indicating that the added intermittency observed in the velocity field is compensated by the variable viscosity in the dissipation rate.
Figure 6.
$(a)$
Probability density function of the alignment of the vorticity unit vector
$\hat {\boldsymbol{\omega }}$
with the eigenvectors
$\hat {\boldsymbol{s}}_1$
(solid line),
$\hat {\boldsymbol{s}}_2$
(dashed line) and
$\hat {\boldsymbol{s}}_3$
(dotted line) of the strain-rate tensor.
$(b)$
Joint histograms of
$\mathcal{Q}$
and
$\mathcal{R}$
. The black curves show where the discriminant of the polynomial equation is zero, i.e.
$27 \mathcal{R}^2/4 + \mathcal{Q}^3 = 0$
. Data are for
$Re_\lambda \approx 300$
.
To look more closely at the local flow altered by the variable viscosity, we compute the alignment of the unit-length eigenvectors
$\hat {\boldsymbol{s}}_i$
of the strain-rate tensor
$\mathcal{S}_{ij}$
with the vorticity
$\boldsymbol{\omega } = \epsilon _{ijk} (\partial u_k/\partial x_j)\, \hat {\boldsymbol{e}}_i$
(Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987), where
$\hat {\boldsymbol{e}}_i$
are the Cartesian unit vectors, and
$\epsilon _{ijk}$
is the Levi–Civita symbol. For an incompressible fluid, the three eigenvalues sum to zero, so the largest eigenvalue is never negative, and its eigenvector
$\hat {\boldsymbol{s}}_1$
corresponds to the direction of elongation in the flow. Similarly, the smallest eigenvalue is negative, and its eigenvector
$\hat {\boldsymbol{s}}_3$
corresponds to the direction of compression in the flow. Finally,
$\hat {\boldsymbol{s}}_2$
completes the tri-orthogonal system. Figure 6(
$a$
) shows the probability density functions of the cosine of the angle between the vorticity
$\hat {\boldsymbol{\omega }}$
and the eigenvectors
$\hat {\boldsymbol{s}}_i$
. There is a strong alignment of vorticity with the intermediate eigenvector
$\hat {\boldsymbol{s}}_2$
, which is frequently attributed to the axial stretching of vortices (Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987), with the maximum probability at
$|\hat {\boldsymbol{\omega }} \boldsymbol{\cdot} \hat {\boldsymbol{s}}_2 | = 1$
; the third eigenvector is largely perpendicular to the vorticity, producing a peak at
$\hat {\boldsymbol{\omega }} \boldsymbol{\cdot} \hat {\boldsymbol{s}}_3 = 0$
, while the first eigenvector
$\hat {\boldsymbol{s}}_1$
exhibits very little correlation with
$\hat {\boldsymbol{\omega }}$
. When the fluid is shear-thinning and
$n$
is reduced, the same pictures remain, with an extremely weak tendency of vorticity aligning even more with the intermediate eigenvector
$\hat {\boldsymbol{s}}_2$
, and becoming more perpendicular to the last eigenvector
$\hat {\boldsymbol{s}}_3$
.
The three principle invariants of the velocity gradient tensor
$\partial u_j/\partial x_i$
can be used to fully characterise the local topology of a flow (Cheng Reference Cheng1996). The first invariant is zero due to incompressibility, the second invariant
$\mathcal{Q} = ({1}/{4})\omega _i\omega _i- ({1}/{2})\mathcal{S}_{ij}\mathcal{S}_{ij}$
represents the balance between strain and vorticity, and the third invariant
$\mathcal{R} = ({1}/{4})\omega _i \mathcal{S}_{i j} \omega _j-({1}/{3}) \mathcal{S}_{i j} \mathcal{S}_{j k} \mathcal{S}_{k i}$
represents the balance between strain and vorticity production. The roots of the polynomial equation
$\Lambda ^3 + \mathcal{Q} \Lambda + \mathcal{R} = 0$
are the eigenvalues of
$\partial u_j/\partial x_i$
. Figure 6(
$b$
) shows the joint probability distributions for all of our flows at large Reynolds number. We also plot a line where the discriminant is zero: below this line, all three eigenvalues are real, and strain dominates the flow, whereas above this line,
$\partial u_j/\partial x_i$
has one real and two complex eigenvalues, and vortices dominate the flow. The top left and bottom right quadrants feature tails in the distribution, which represent stretched vortices and areas where the flow compresses along a single axis (Cheng Reference Cheng1996). The shape remains overall unchanged for the non-Newtonian fluids, but with the probability falling off much more sharply when
$n$
is reduced, thus indicating weaker velocity gradients, compensated by the local variation of the viscosity.
4. Conclusions
In this work we have extended the Kolmogorov theory to shear-thinning fluids with variable viscosity, with the main quantities emerging being
$\langle \mu \rangle$
,
$\langle \varepsilon \rangle$
and
$n$
, and with the Kolmogorov theory for a Newtonian fluid recovered when
$n=1$
. By analysing data from direct numerical simulations, we have shown that the variable viscosity brings no significant changes to the large and intermediate scales of the flow, while affecting only the smallest scale. We have shown that these fluids exhibit a universal inertial range with power-law spectrum scaling
$k^{-5/3}$
, extending in range as the level of shear-thinning is increased, and caused by a strengthening of the nonlinear energy flux. As a consequence, the dissipation range shrinks. Remarkably, small-scale universality is lost, with a remaining dependency on the power-law index
$n$
. The modified flows exhibit enhanced intermittency in the velocity field, with local flow structures having similar directionality to the Newtonian case, but with weaker intensity, the two effects being compensated by the variation of the viscosity. The results hold for different
$\mu _0/\mu _\infty$
,
$\mathcal{K}$
and
$n$
, as long as the viscosity falls in the purely power-law region of the rheology. In particular, the effect of
$\mathcal{K}$
, while significant in changing the viscosity quantitatively, can be collapsed into a change in the actual Reynolds number experienced by the flow.
The results discussed in this work are relevant when studying polymeric flows, where often elasticity and shear-thinning are not decoupled. Since shear-thinning promotes the nonlinear energy flux, while viscoelasticity weakens it, complex effects can arise when studying polymeric flows at both large and small Reynolds numbers: (i) shear-thinning acts against elasticity at large Reynolds numbers, potentially hiding the elasto-inertial range of scales (Rosti et al. Reference Rosti, Perlekar and Mitra2023); (ii) shear-thinning can promote instabilities at low Reynolds numbers (see e.g. for a jet Soligo & Rosti Reference Soligo and Rosti2023), which are, however, not elastic, with the resulting turbulent flow which may be inertia-dominated, notwithstanding the low Reynolds number (Amor et al. Reference Amor, Soligo, Mazzino and Rosti2024).
Acknowledgements
The author acknowledges the computer time provided by the Scientific Computing section of Research Support Division at OIST and by HPCI, under the Research Project grants hp210269, hp220099, hp230018, hp250021 and hp250035.
Funding
The research was supported by the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding to M.E.R. from the Cabinet Office, Government of Japan. M.E.R. acknowledges funding from the Japan Society for the Promotion of Science (JSPS), grants 24K17210 and 24K00810.
Declaration of interests
The authors report no conflict of interest.
Open access
