1. Introduction
Turbulence rapidly drains energy from a fluid flow. Therefore, in many applications – such as aviation, naval and road transport and industrial processes – a significant portion of supplied energy is lost to undesired and difficult-to-control turbulent motion. Minimising turbulent drag is thus a key strategy for improving energetic efficiency. A spectacular example of turbulent drag reduction is an effect first experimentally studied by B.A. Toms in the 1940s (Toms Reference Toms1949, Reference Toms1977). Adding a small amount of polymers to a turbulent pipe or channel flow can reduce the turbulent drag enormously, leading thereby to important savings of energy.
Despite the large number of investigations dedicated to the subject (see reviews of older contributions in Lumley Reference Lumley1973 and Virk Reference Virk1975 and more recent work in Procaccia, L’vov & Benzi Reference Procaccia, L’vov and Benzi2008, White & Mungal Reference White and Mungal2008; Xi Reference Xi2019 and Dubief, Terrapon & Hof Reference Dubief, Terrapon and Hof2023), the precise effect remains obscure. It is realised that the polymers interact with turbulent velocity fluctuations in a certain way, which allows us to reduce the overall energy dissipation and the momentum flux from the fluid towards the wall. Early efforts of De Gennes and Lumley pointed to possible mechanisms associated with the elastic properties of polymers (De Gennes Reference De Gennes1986; Tabor & De Gennes Reference Tabor and De Gennes1986) or the modification of the effective viscosity (Lumley Reference Lumley1973), an idea further explored, for instance, in L’vov et al. (Reference L’vov, Pomyalov, Procaccia and Tiberkevich2004) and Ryskin (Reference Ryskin1987b ). More recent studies have focused on spatio-temporal intermittency, or hibernating states associated with events of reduced drag reduction (Graham Reference Graham2004; Xi & Graham Reference Xi and Graham2012; Whalley et al. Reference Whalley, Park, Kushwaha, Dennis, Graham and Poole2017). Despite significant progress, none of these theories is entirely satisfactory or predictive.
Even though the precise mechanism behind turbulent drag reduction by polymers is unknown, the continued research efforts during eight decades have yielded a wealth of insights. Indeed, several key features of dilute polymer-containing flows near boundaries are now well established, and here we highlight some of these features.
Polymer-laden flows are visco-elastic, which adds to the difficulty of normal, viscous flows. Due to their elastic nature, the polymers need an extensional flow to get stretched. If the typical time scale of the flow is too large, the polymers will relax to their equilibrium coiled state. The dimensionless number which compares the elastic time scale with the flow time scale is called the Weissenberg number
$Wi$
. At low Weissenberg numbers, polymers have minimal influence on the flow.
As
$Wi$
increases beyond the coil–stretch transition (De Gennes Reference De Gennes1974; Watanabe & Gotoh Reference Watanabe and Gotoh2010), the drag reduction becomes more significant, reaching, for large enough values of
$Wi$
, an upper limit associated with a flow configuration known as the maximum drag reduction (MDR) state (Virk, Mickley & Smith Reference Virk, Mickley and Smith1970; Virk Reference Virk1975). Notably, there appears to be a functional relationship between
$Wi$
and the extent of drag reduction in parallel shear flows (Owolabi, Dennis & Poole Reference Owolabi, Dennis and Poole2017). For large values of the drag reduction, and up to the MDR state, the mean-velocity profile in the near-wall region transitions from the logarithmic profile – characteristic of turbulent flow in Newtonian fluids – to a steeper profile. Historically, this profile was approximated using an alternative logarithmic expression (Virk Reference Virk1975). However, more precise measurements have shown that the profile deviates from a logarithmic form (White, Dubief & Klewicki Reference White, Dubief and Klewicki2012), exhibiting a convex shape in log–linear representation (Ptasinski et al. Reference Ptasinski, Nieuwstadt, Van Den Brule and Hulsen2001; Owolabi et al. Reference Owolabi, Dennis and Poole2017).
In this study, we focus on the role of vortex-stretching reduction in drag reduction. The importance of vortex stretching has been mentioned in numerous studies over the years. As early as the works of Gadd (Reference Gadd1968) and Landahl (Reference Landahl1973), vortex-stretching suppression was identified as a significant factor. Subsequent studies (Sureshkumar, Beris & Handler Reference Sureshkumar, Beris and Handler1997; Yarin Reference Yarin1997) expanded on this idea, while experiments demonstrated that material-line stretching was reduced in polymer-laden flows (Liberzon et al. Reference Liberzon, Guala, Lüthi, Kinzelbach and Tsinober2005). Recent numerical simulations confirm that polymers attenuate vortex stretching (ur Rehman et al. Reference ur Rehman, Lee and Lee2022), and experiments have shown that this attenuation is central to the mechanism of polymer drag reduction (Warwaruk & Ghaemi Reference Warwaruk and Ghaemi2024).
Polymers reduce vortex stretching due to two key effects. First, it has been shown that rod-like passive particles or fibres align with the vorticity vector in turbulent flow (Pumir & Wilkinson Reference Pumir and Wilkinson2011; Ni, Ouellette & Voth Reference Ni, Ouellette and Voth2014), and this is expected to hold true for polymers as well. Second, the presence of polymers significantly increases the extensional viscosity (Metzner & Metzner Reference Metzner and Metzner1970; Hinch Reference Hinch1977; Lindner, Vermant & Bonn Reference Lindner, Vermant and Bonn2003). From these two observations – the alignment of polymers with vorticity and the increased extensional viscosity – it follows that polymers attenuate vortex stretching. This attenuation reduces drag since vortex stretching is one specific part of the nonlinear term, responsible for generating drag (Li et al. Reference Li, Fan, Modesti and Cheng2019).
The idea that drag reduction is linked to increased extensional viscosity is not new. Most theoretical studies acknowledge that this increase is a principal effect of polymers on visco-elastic flow. In purely extensional flows, this increase can be interpreted as an effective viscosity enhancement (Ryskin Reference Ryskin1987a ). When flows are not purely extensional, modifying the effective viscosity inhomogeneously provides a coarse model for the influence of polymers in wall-bounded flow, albeit without accounting for the difference between extensional and shear stresses. This approach underpins Lumley’s theory (Lumley Reference Lumley1973), which can be refined using local energy balances and dimensional arguments to estimate viscosity changes (L’vov et al. Reference L’vov, Pomyalov, Procaccia and Tiberkevich2004; Procaccia et al. Reference Procaccia, L’vov and Benzi2008). However, firstly, this does not explain observations in homogeneous shear flow (Robert et al. Reference Robert, Vaithianathan, Collins and Brasseur2010; Benzi & Ching Reference Benzi and Ching2018) where the viscosity should be statistically homogeneous and, secondly, turbulence is far from purely extensional. Incorporating a separate treatment of the two different stresses into the Navier–Stokes equations remains challenging. Only recently has progress been made on this for the case of two-dimensional flows (Oliveira Reference Oliveira2024; see also Poole Reference Poole2023).
The role of elasticity in polymer-laden turbulence is still a debated subject. At moderate and high Reynolds numbers it is the fully stretched polymers which are responsible for the drag reduction (Serafini et al. Reference Serafini, Battista, Gualtieri and Casciola2022), thereby suggesting that elastic effects are not necessary to explain drag reduction beyond the coil–stretch transition. However, around the MDR state, in particular at low Reynolds numbers, elastic instabilities might play an important role in maintaining a marginally unstable turbulent state (Choueiri, Lopez & Hof Reference Choueiri, Lopez and Hof2018; Dubief et al. Reference Dubief, Terrapon and Hof2023). We will come back to this point in the conclusion section. For the moment, we will ignore the effects of elasticity and consider flows beyond the coil–stretch transition.
In the present investigation, we avoid the complexities associated with altering the system’s viscosity or modifying the viscous stress tensor. Instead, we focus directly on the vortex-stretching mechanism at the level of the governing equations. This approach neither relies on purely phenomenological arguments nor involves detailed simulations of visco-elastic turbulence. The strength of our method is its ability to isolate a specific aspect of polymer–turbulence interaction for specific analysis. Such an approach necessarily omits certain features of the interaction (such as elastic instabilities (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013) and the effect of the weight distribution of the polymers (Brandfellner et al. Reference Brandfellner, Muratspahić, Bismarck and Müller2024; Serafini et al. Reference Serafini, Battista, Gualtieri and Casciola2025), etc.), but it will prove valuable since it allows us to derive an analytical prediction of the inertial velocity profile. We will then compare this prediction with specifically designed numerical experiments, before comparing it with laboratory experiments of polymer-laden turbulent pipe flow.
2. A law of the wall for turbulence with tamed vortex stretching
2.1. Derivation of the mean profile
To model the influence of reduced vortex stretching, we write the evolution equation of the fluid vorticity in the form
\begin{equation} \frac {\partial \boldsymbol \omega }{\partial t}+\boldsymbol u\boldsymbol{\cdot } \boldsymbol{\nabla }{\boldsymbol \omega} =(1-\lambda )\boldsymbol \omega \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol u -\lambda \boldsymbol{\nabla }P_\omega +\nu \Delta \boldsymbol \omega , \end{equation}
where the vorticity
$\boldsymbol \omega$
is the curl of the velocity
$\boldsymbol u$
,
$\nu$
the kinematic viscosity and
$P_\omega$
a Lagrange multiplier, ensuring
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol \omega =0$
(Wu, David & Bos Reference Wu, David and Bos2023)
$\lambda$
is here a control parameter which we will discuss below. A first property of this equation is that the left-hand side, representing the advection of the vorticity, is not modified. Intuitively, this represents the fact that the advection of a small fluid volume should be largely independent (at leading order) of what molecular particles, particles or fibres, the fluid particle contains. We hypothesise that this is different for the vortex-stretching term, since it is the term responsible for enstrophy generation, an effect dominated by extensional motion aligned with the vorticity. Indeed, the extensional stress is dramatically influenced by the presence of polymers in the flow.
In the case
$\lambda =0$
we retrieve the classical Navier–Stokes vorticity equation. The system (2.1) for
$\lambda =1$
was introduced in Bos (Reference Bos2021) and Wu & Bos (Reference Wu and Bos2021) and studied in statistically homogeneous flow. In the present investigation we consider (2.1) in wall-bounded flow as the simplest model for polymer-laden turbulence; we suppress vortex stretching, partially or completely, by allowing values
$0\leqslant \lambda \leqslant 1$
.
Before carrying out numerical experiments, we will first extend the classical arguments of Prandtl and von Kármán (von Kármán Reference von Kármán1930) to our system, in order to anticipate how the near-wall scaling is modified by the attenuation of vortex stretching. For this, we follow the same reasoning as von Kármán in his 1930 paper, which introduced the logarithmic law of the wall (or log law).
We focus on fully developed channel flow, with
$x,y,z$
the streamwise, wall-normal and spanwise directions, respectively. All average quantities are stationary and homogeneous in the streamwise and spanwise directions. The mean-velocity and vorticity fields are
\begin{align} \langle \boldsymbol u \rangle =U(y)\boldsymbol e_x\,\,\,\,\,\,\,\,\langle \boldsymbol \omega \rangle =-\frac {{\rm d}U(y)}{{\rm d}y} \boldsymbol e_z, \\[-18pt] \nonumber \end{align}
respectively. From (2.1) we derive the exact mean vorticity balance
$\partial _t \langle \omega _z \rangle$
\begin{eqnarray} 0=\lambda \frac {{\rm d}\langle \omega _z^{\prime}v' \rangle }{{\rm d} y}+(\lambda -1)\frac {d^2\langle u'v' \rangle }{{\rm d} y^2} +\nu \frac {d^3 U}{{\rm d} y^3}. \end{eqnarray}
Neglecting viscous effects as usual in the derivation of the log law, we proceed by modelling the turbulent fluxes of mean vorticity and momentum,
$\langle \omega _z^{\prime}v' \rangle$
and
$\langle u'v' \rangle$
, using a classical mixing-length model
\begin{eqnarray} \langle a'v' \rangle =\nu _T\frac {\partial \langle a\rangle }{\partial y}, \textrm { with } \nu _T=(\kappa y)^2 \left |\frac {\partial U}{\partial y}\right |, \end{eqnarray}
where
$a$
can be either the velocity or vorticity. This corresponds to a mixing length
$l= \kappa y$
proportional to the distance to the wall (von Kármán Reference von Kármán1930). In principle, different values can be used for the constant in the mixing-length estimate for vorticity or velocity (Hinze Reference Hinze1975), but this will not qualitatively change the outcome of this analysis. We obtain, since
${\rm d}\langle \omega _z\rangle /{\rm d}y=-{\rm d}^2U/{\rm d}y^2$
,
\begin{eqnarray} 0=\kappa ^2 \lambda y^2 \frac {{\rm d}^2U}{{\rm d}y^2}\frac {{\rm d}U}{{\rm d}y}+(1-\lambda )\kappa ^2 \frac {{\rm d}}{{\rm d}y}\left (y^2\frac {{\rm d}U}{{\rm d}y}\frac {{\rm d}U}{{\rm d}y}\right ). \end{eqnarray}
Equation (2.5) can be solved, yielding the solution
\begin{eqnarray} U^+(y^+)=\frac {1}{\kappa }\frac {{y^+}^{\varGamma }-1}{\varGamma }+U^+_{0}, \end{eqnarray}
with
\begin{equation} \varGamma =\lambda /(2-\lambda ) ,\end{equation}
where quantities are expressed in wall units (
$U^+=U/u_\tau$
, with
$u_\tau =(\boldsymbol{\nabla }\langle p\rangle h/\rho )^{1/2}$
,
$y^+=y u_\tau /\nu$
) with
$\rho$
the density. In the derivation, an integration constant is fixed, as is customary, by the observation that, in the inertial layer of Newtonian channel flow turbulence,
$|\langle u'v'\rangle |\approx u_\tau ^2$
(Pope Reference Pope2000).
According to this phenomenological analysis we thus find a power-law scaling with an exponent
$\lambda /(2-\lambda )$
for
$0\lt \lambda \leqslant 1$
. For
$\lambda =0$
we find as a special case the logarithmic profile. Indeed, in the limit
$\varGamma \downarrow 0$
, we find that
\begin{eqnarray} \lim _{\varGamma \rightarrow 0}U^+(y^+)=\frac {1}{\kappa }\ln (y^+)+U^+_0. \end{eqnarray}
The general
$\lambda$
-dependent mean-velocity profile (2.6) is the main analytical prediction which we will compare first with direct numerical simulations of (2.1) and subsequently with experimental results from literature on polymer turbulence.
2.2. Numerical experiments of turbulence with reduced vortex stretching
To assess our predictions, we carry out Direct Numerical Simulations (DNS) of (2.1) in plane channel flow. The walls are separated by a width
$2h$
, the width of the domain is
$\pi h$
and the length
$2\pi h$
. Boundary conditions are no slip at the wall and periodic in the streamwise and spanwise directions. The code is based on a classical Fourier–Chebyshev formulation with a resolution of
$N_x=256,N_y=192,N_z=192$
grid points. Further details can be found in Appendix A. Computations are carried out for
$R_\tau =395$
at a constant mass-flow rate. Simulations are run until a statistically steady state is obtained for the Navier–Stokes case
$\lambda =0$
, as measured by observing the wall shear stress. Starting from this steady state,
$\lambda$
is varied in the range
$\lambda \in [0,0.5]$
with steps of
$\delta \lambda =0.05$
. For simulations with
$\lambda \leqslant 0.35$
turbulent steady states are observed, while for
$\lambda \geqslant 0.4$
the flows relaminarise.
In figure 1 we show flow visualisations using the
$Q$
-criterion (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988) for
$Q\equiv \Delta p/2=1$
. Comparing the Navier–Stokes case (a) with the case
$\lambda =0.35$
(b), it is observed that, in the second case, the vortex intensity has reduced dramatically. Indeed, the vortex-stretching term is responsible for enstrophy production, and reducing it leads to the intuitive effect of less vortical activity.
Plane-Poiseuille flow driven by an imposed pressure gradient in the x-direction. Vortical motion is amplified in three-dimensional turbulence by vortex stretching. The flow visualisations show iso-surfaces of the
$Q$
-criterion for
$Q=1$
(see main text), coloured by the local streamwise velocity. (a) Traditional channel flow. (b) Channel flow with reduced vortex stretching (
$\lambda =0.35$
in (2.1)).
Comparison of analytical predictions with our DNS results. (a) Mean-velocity profiles in log–linear representation. The values of
$U_0^+$
are for
$\lambda \in [0,0.35]$
:
$5.1; 5.9; 6.6; 7.8; 8.3; 9.8; 9.9; 10.3$
. (b) The generalised indicator function
$\varXi =(y^+)^{1-\lambda /(2-\lambda )}{\rm d}U^+(y^+)/{\rm d}y^+$
, (2.10) for the different velocity profiles in (a).
In figure 2(a) we show mean-velocity profiles for the turbulent cases, determined by averaging during the statistically steady state. Also shown is our prediction, where
$\kappa =0.4$
and the only adjustable parameter is the velocity offset
$U^+_0(\lambda )$
, which is fitted to obtain a best fit in the region where the generalised indicator function (see below) is in the vicinity of
$1/\kappa$
. Even though the Reynolds number is rather small, the predicted profiles describe the data well. The relation between the values
$U^+_0(\lambda )$
and
$\lambda$
is approximately linear
\begin{equation} U^+_0(\lambda )=U^+_0(0)+\lambda b_U, \end{equation}
with
$U^+_0(0)=5.2\pm 0.2$
and
$b_U=15.8\pm 1$
.
In order to assess our predictions independently from the virtual origin
$U_0(\lambda )$
, we introduce a generalised indicator function
\begin{equation} \varXi [U,y]=y^{1-\varGamma }\frac {{\rm d}U(y)}{{\rm d}y} ,\end{equation}
where
$\varGamma =\lambda /(2-\lambda )$
. For the value
$\lambda =0$
this simplifies to the commonly used indicator function to verify the existence of a logarithmic profile (e.g. White et al. Reference White, Dubief and Klewicki2012; Laadhari Reference Laadhari2019). For
$\lambda \neq 0$
, this function should allow us to identify the power-law behaviour (2.6) by a plateau of value
$1/\kappa$
independently of the adjustable parameters other than the value of
$\lambda$
, which is fixed for each simulation.
Even though the Reynolds number is too low to have a substantial inertial profile, figure 2(b) shows that all 7 simulations coincide near the value
$1/\kappa \approx 2.5$
around
$y^+=100$
, supporting the idea that
$\kappa$
plays a central role, irrespective of the amount of vortex stretching that is suppressed in the governing equations.
3. Drag reduction and vortex stretching
3.1. Similarities between the modified Navier–Stokes equations and visco-elastic turbulence
The extended von Kármán phenomenology appears to capture the behaviour of (2.1), elucidating how reduced vortex stretching modifies the mean-velocity profile. The next question is: How relevant is (2.1) for modelling polymer turbulence?
Since (2.1) lacks derivation from a polymer stress constitutive relation (an important question for future research), its relevance must be inferred by comparing its properties with those of real flows. Before attempting such comparison, let us discuss certain similarities between polymer-laden turbulence and (2.1).
First, the steepened velocity profile in our system qualitatively resembles profiles in polymer turbulence. Although our approach does not model the coil–stretch transition (De Gennes Reference De Gennes1974), the continuous profile steepening aligns with high drag reduction regimes (Warholic, Massah & Hanratty Reference Warholic, Massah and Hanratty1999).
Second, our system stabilises at a marginally stable state for
$\lambda \approx 0.35$
, with relaminarisation at higher values. While polymer turbulence also reaches an asymptotic (MDR) state, the flow does not relaminarise at high Reynolds numbers. It is plausible that exceeding MDR suppresses vortex stretching to the point of relaminarisation, at which polymers relax to their coiled equilibrium state. This would cause the flow to revert to Newtonian-like properties, reintroducing instabilities and turbulence. Local spatio-temporal relaminarisation could correspond to hibernating turbulent states (Graham Reference Graham2004; Xi & Graham Reference Xi and Graham2012).
Third, drag reduction can occur for both elastic and rod-like polymers (Berman Reference Berman1978; Paschkewitz et al. Reference Paschkewitz, Dubief, Dimitropoulos, Shaqfeh and Moin2004; Japper-Jaafar, Escudier & Poole Reference Japper-Jaafar, Escudier and Poole2009), suggesting that polymer-induced drag reduction is not purely elastic in origin. Our approach, lacking relaxation time scales (inherent to elastic effects) is consistent with this observation.
Fourth, simulations of visco-elastic isotropic turbulence and turbulence without vortex stretching reveal shared features. Simulations employing the finite extensible nonlinear elastic model with Peterlin’s approximation show that, at high Weissenberg numbers, the kinetic energy spectrum transitions from a
$k^{-5/3}$
scaling to
$k^{-3}$
, with
$k$
the wavenumber (Valente et al. Reference Valente, Da Silva and Pinho2014, Reference Valente, da Silva and Pinho2016), as also observed in channel flows (Mitishita, Elfring & Frigaard Reference Mitishita, Elfring and Frigaard2023). This
$k^{-3}$
scaling, linked to enstrophy conservation, is similarly observed in turbulence without vortex stretching (Wu & Bos Reference Wu and Bos2022). Furthermore, advanced simulations (Watanabe & Gotoh Reference Watanabe and Gotoh2013) demonstrate that polymers suppress enstrophy production, promoting a conservative enstrophy cascade at high Reynolds numbers.
These similarities motivate further scrutiny of the connection between polymer drag reduction and (2.1). We therefore compare its predictions with experimental results.
Comparison of analytical predictions with experimental results. (a) Mean-velocity profile in pipe flow at MDR (Owolabi et al. Reference Owolabi, Dennis and Poole2017) compared with our power-law estimate. (b) Fanning friction factor in Prandtl–von Karman coordinates. Experimental results from smooth pipe measurements at MDR (Virk et al. Reference Virk, Mickley and Smith1970).
3.2. Comparison with experiments of visco-elastic shear flow
We compare our predictions with measures in pipe flow, which is, historically, the most investigated case. We first compare with the velocity profile in pipe flow at MDR measured experimentally at the University of Liverpool (Owolabi et al. Reference Owolabi, Dennis and Poole2017) in figure 3(a). Fitting our prediction (2.6) in the range
$100\leqslant y^+\leqslant 250$
, we obtain
$\varGamma \equiv \lambda /(2-\lambda )=0.349\pm 0.001$
and
$U^+_0=3.5\pm 0.13$
. Agreement at these scales is convincing and clearly superior to Virk’s approximation.
A larger amount of data are available on the friction factor, which is more easily measured in experiments. We will here compare with the pioneering compilation of results by Virk et al. (Reference Virk, Mickley and Smith1970). The friction factor is determined by measuring the pressure drop
$\delta p$
over a pipe of length
$L$
. Introducing the bulk velocity
$U_b$
(velocity averaged over the pipe cross-section of radius
$R$
) the friction factor is
\begin{equation} f=2\left (\frac {u_\tau }{U_b}\right )^2, \end{equation}
where
$u_\tau =\sqrt {R\delta p/(2\rho L)}$
. Experiments are then performed to measure the value of
$f$
as a function of the bulk Reynolds number,
$R_b=2R U_b/\nu$
. Results are traditionally plotted, in Prandtl–von Kármán representation,
$f^{-1/2}$
as a function of
$R_b f^{1/2}(=2\sqrt {2}R_\tau )$
.
We determine the bulk velocity by integrating the velocity profile over the pipe cross-section, yielding, after a change of variables
$r=R-y$
and
$y^+=((R-r)/R)R_\tau$
, the relation
\begin{equation} f^{-1/2}=\frac {\sqrt {2}}{ R_\tau }\int _0^{R_\tau }\left (1-\frac {y^+}{R_\tau }\right )U^+(y^+){\rm d}y^+. \end{equation}
We now use expression (2.6) for
$\lambda \neq 0$
, to obtain
\begin{equation} f^{-1/2}=\frac {2 \left (\left [R_b f^{1/2}/\sqrt {8}\right ]^{\varGamma }-1\right )/\varGamma -(\varGamma +3)}{\sqrt {2}\kappa (\varGamma +1)(\varGamma +2)}+\frac {U^+_0}{ \sqrt {2}} .\end{equation}
We show this expression in figure 3(b). Fitting the data we find
$\varGamma =0.348 \pm 0.003$
and
$U^+_0=5.6\pm 0.4$
. Also shown are the Newtonian expression
$f^{-1/2}(x)=4 \log (x)-0.4$
and Virk’s relation
$f^{-1/2}(x)=19 \log (x)-32.4$
.
What we can say is that, for large values of the Reynolds number, the agreement between experiments of the friction factor and our theory is at least as good as using Virk’s log law, with the difference that our fit is based on a theoretical expression for the velocity profile (2.6), valid for both Newtonian and visco-elastic flow, whereas Virk’s expression is an empirical fit to the MDR data. For the velocity profile for the scales
$100\lt y^+\lt 250$
the agreement with our theory is clearly superior (figure 3
a).
It is remarkable that we obtain, within error bars, exactly the same value for
$\varGamma$
in figures 3(a) and 3(b). This supports the validity of (2.6) to describe the mean profile at MDR. The slightly different value of
$U_0$
is not surprising since we obtain the estimate for
$U_b$
by integrating over the full pipe cross-section, although the power law does not hold for
$y^+$
near the wall and the centre. Another nice feature is that (3.3) holds also for the Newtonian case, when the value
$\varGamma =0$
is used.
A question we have not addressed is the link between the Weissenberg number and the vortex-stretching parameter
$\lambda$
. At present, we have no direct physical relation between
$\lambda$
and
$Wi$
derived from a visco-elastic constitutive relation. In order to obtain a tentative link, we compare the results of our system with the compilation of Owolabi et al. (Reference Owolabi, Dennis and Poole2017), who showed that the amount of drag reduction in different pipes and channels could be linked by a close-to-universal relation to
$Wi$
. The drag reduction, for a fixed mass-flow rate, is defined as
\begin{equation} \textit{DR}=\left (1-\left (\frac {u_\tau }{u_\tau ^{(0)}}\right )^2 \right ) ,\end{equation}
where
$u_\tau ^{(0)}$
is the friction velocity for the reference case. In figure 4 we show the experimental data points for the interval
$0\lt Wi\lt 6$
. We compare with the results of our DNS for different values of
$\lambda$
. To attempt a comparison and in the absence of a physical model, we try to compare the data using the linear relation
\begin{equation} Wi-Wi_{c}=a \lambda , \end{equation}
which is the simplest possible non-trivial relation, taking into account that the suppression of vortex stretching should increase with the Weissenberg number, at least for values between the coil–stretch transition (at the critical Weissenberg number
$Wi_c$
) and MDR. We use the value
$Wi_c=0.5$
and the value
$a=10$
, which seem to give a good qualitative agreement.
Drag reduction as a function of the Weissenberg number. Comparison of experimental results of Owolabi et al. (Reference Owolabi, Dennis and Poole2017) with our analytical predictions and numerical results for DNS with reduced vortex stretching. For the comparison we use a linear relation,
$Wi-Wi_c=10\lambda$
between the reduction of vortex stretching
$\lambda$
and
$Wi$
. The horizontal dashed line represents the MDR value (
$\textit{DR}=64\,\%$
) corresponding to the fit to the experimental results in Owolabi et al. (Reference Owolabi, Dennis and Poole2017).
To compare with the theory, we integrate the mean-velocity profile (2.6) over half of the channel to get an estimate of the bulk velocity
\begin{eqnarray} U_b&=&\frac {u_\tau }{h}\int _0^h U(y) {\rm d}y\nonumber \\ &=&\frac {u_\tau }{\kappa (\varGamma +1)}\left [\frac {\left (\dfrac {hu_\tau }{\nu }\right )^{\varGamma }-1}{\varGamma }-1\right ]+u_\tau U_0(\lambda ), \end{eqnarray}
yielding an implicit equation for
$u_\tau$
for fixed
$U_b$
, which can be numerically solved, using (2.9). The resulting values of
$u_\tau$
are used to compute the drag reduction, shown in figure 4. The tendency is the same as the DNS and the experiments. Quantitative disagreement with the DNS can be attributed to the assumption of using (2.6) over the complete channel, ignoring deviations near the wall and centre of the domain. Using this relation between
$\lambda$
and
$Wi$
, we find that the value
$\varGamma = 0.349$
we obtained in figure 3 corresponds to
$Wi_\textit{MDR}\approx 5.6$
, which is not inconsistent with the results in Owolabi et al. (Reference Owolabi, Dennis and Poole2017).
4. Conclusion and future issues
The complexity of polymer-laden turbulence is rooted in the multi-scale interaction of turbulence with polymer molecules (Koide & Goto Reference Koide and Goto2024) combined with the presence of walls, introducing both scale and position dependences. While advanced numerical experiments now provide access to detailed flow quantities in realistic geometries (Serafini et al. Reference Serafini, Battista, Gualtieri and Casciola2022), progress in understanding requires drastic simplifications of the flow physics to isolate specific features – an approach we have adopted.
Introducing the new system of governing equations (2.1) for the investigation of drag reduction by polymers opens several avenues for exploration. A first and important one is the assessment of the assumption that strong extensional viscosity primarily acts to damp the term
$\boldsymbol \omega \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol u$
in the vorticity dynamics. Such an assessment could start from the framework introduced in Oliveira (Reference Oliveira2024). Further comparisons should investigate the implications of (2.1) for the detailed inhomogeneous and anisotropic turbulence statistics (Escudier, Nickson & Poole Reference Escudier, Nickson and Poole2009).
Perhaps even more intriguingly, investigating the self-sustaining mechanism of wall turbulence (De Angelis et al. Reference De Angelis, Casciola and Piva2002; Waleffe Reference Waleffe1997) and the linear and nonlinear instability properties of shear flows with reduced vortex stretching, could provide theoretical insights into the nature of the MDR state. While at low Reynolds numbers certain results indicate that elastic instabilities might determine the precise asymptotic value of MDR (Dubief et al. (Reference Dubief, Terrapon and Hof2023) and references therein), this question is not settled at high Reynolds numbers. The MDR asymptote at higher Reynolds numbers might be determined by a subcritical transition (Morozov & van Saarloos Reference Morozov and van Saarloos2005, Reference Morozov and van Saarloos2007), and as such, its precise value is certainly not easy or impossible to determine by simple analytical reasoning.
In this context, an interesting observation in our study is that we observe laminarisation at relatively large values of the Reynolds number for sufficiently high values of
$\lambda$
, in the absence of elastic effects. Let us imagine what elastic effects would add in such a relaminarised state in an actual visco-elastic fluid. There, polymers would relax back to their equilibrium state, which would have two direct influences: first, the relaxation could reintroduce kinetic energy into the flow when the polymer stress changes sign, as observed in isotropic visco-elastic turbulence (Nguyen et al. Reference Nguyen, Delache, Simoëns, Bos and El Hajem2016; Valente, Da Silva & Pinho Reference Valente, Da Silva and Pinho2014) at small scales and high Weissenberg number, or in low Reynolds number shear flows due to elastic instabilities Choueiri et al. (Reference Choueiri, Lopez and Hof2018). This energy flux could possibly perturb the relaminarised state. Second, after relaxation, the extensional viscosity would recover its Newtonian value and the flow would become unstable again. Both effects would contribute to sustaining a marginally unstable state, which would act as a visco-elastic fluid at the MDR asymptote.
These arguments support the ideas that, firstly, (extensionally) viscous effects are the main actor in reaching the MDR state. Secondly, that elastic effects are essential to maintain a turbulent state once the MDR asymptote is reached. A more sophisticated model than (2.1) is needed to assess this reasoning and the development of such a model is left for further research.
Acknowledgements
We gratefully acknowledge discussions with R. Poole, F. Laadhari, B. Miquel, J. Xie and C. Xu. This work was supported by the National Natural Science Foundation of China (project approval nos. 12372214 and U2341231). The data that support the findings of this article are openly available.
For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical method
We carry out simulations of plane-Poiseuille flow. The half-width of the channel is denoted
$h$
, the computational domain has dimensions of
$L_x \times L_y \times L_z = 2\pi h \times 2h \times \pi h$
and employs
$N_x \times N_y \times N_z = 256 \times 192 \times 192$
grid points. Uniform grids are applied in the streamwise and spanwise directions, while a non-uniform grid is used in the wall-normal (y) direction, defined by
$y_j = \cos (\pi j / N_y)$
for
$j = 0, 1, \ldots , N_y$
. Note that in the numerical method
$y=0$
is the centreline of the channel. The initial mean-velocity profile is set as
$U(y) =(3F/4)(1-y^2)$
. By maintaining a constant mass-flow rate
$F=2$
, the flow achieves a Reynolds number of
$R_\tau = hu_\tau /\nu =395$
, with
$\nu$
the kinematic viscosity and
$u_\tau = \sqrt {\tau _w/\rho }$
the friction velocity. Here,
$\tau _w$
represents wall shear stress and
$\rho$
denotes density.
We conduct simulations with different values of
$\lambda$
, starting from a fully developed flow of conventional Navier–Stokes turbulence (
$\lambda =0$
). The numerical method for conventional turbulence has been validated in previous studies (Xu, Zhang & Nieuwstadt Reference Xu, Zhang and Nieuwstadt1996; Fang et al. Reference Fang, Shao, Bertoglio, Lu and Zhang2011; Chen & Fang Reference Chen and Fang2024). Spatial discretisation is conducted using the Fourier–Galerkin and Chebyshev–Gauss–Lobatto collocation methods: streamwise and spanwise directions are expanded using Fourier series, while the wall-normal direction utilises Lagrange interpolation polynomials at Chebyshev–Gauss–Lobatto collocation points.
Before explaining the time advancement, let us start by formulating the evolution equation for the velocity. Equation (A1) of the manuscript is formulated for the evolution of the vorticity, whereas the numerical method we use is based on the time advancement of the velocity. We therefore use the Biot–Savart operator to write
\begin{equation} \frac {\partial \boldsymbol{u}}{\partial t}=\unicode{x1D6E5}^{-1}\boldsymbol{\nabla }\times \left (\frac {\partial \boldsymbol{\omega }}{\partial t}\right ). \end{equation}
We now write two evolution equations. The first equation corresponds to Navier–Stokes turbulence
\begin{equation} \left \{ \begin{aligned} \frac {\partial \boldsymbol{u}}{\partial t} &= \boldsymbol{u} \times \boldsymbol{\omega }- \boldsymbol{\nabla }\varPi + \nu \nabla^{2} \boldsymbol{u} ,\\ \varPi &= \frac {\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{u}}{2} + \frac {p}{\rho }. \end{aligned} \right . \end{equation}
The second one corresponds to turbulence without vortex stretching
\begin{equation} \left \{ \begin{aligned} \frac {\partial \boldsymbol{u}}{\partial t}&=\unicode{x1D6E5}^{-1}\boldsymbol{\nabla }\times ((\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\omega })-\boldsymbol{\nabla }\varPi + \nu \nabla^{2} \boldsymbol{u} ,\\ \varPi &= \frac {\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{u}}{2} + \frac {p}{\rho }, \end{aligned} \right . \end{equation}
where we used (A1). Adding now the first and the second equation in proportions
$\lambda$
and
$1-\lambda$
, we obtain the evolution of the velocity with partially removed vortex stretching
\begin{equation} \left \{ \begin{aligned} \frac {\partial \boldsymbol{u}}{\partial t}&=\lambda \left (\unicode{x1D6E5}^{-1}\boldsymbol{\nabla }\times ((\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\omega })\right )+(1-\lambda )(\boldsymbol{u}\times \boldsymbol{\omega })\\ &\quad-\boldsymbol{\nabla }\varPi + \nu \nabla^{2} \boldsymbol{u} ,\\ \varPi &= \frac {\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{u}}{2} + \frac {p}{\rho }. \end{aligned} \right . \end{equation}
Now that we have a velocity-evolution equation we will describe the time advancement in detail.
The time-advancement scheme is based on a semi-implicit, three-step backward differentiation formula (that is,
$J_e=3$
) (Moin, Reynolds & Ferziger Reference Moin, Reynolds and Ferziger1978; Fan & Dong Reference Fan and Dong2022). During each time step, a third-order time-splitting method is employed, including three substeps that sequentially process the convective, pressure and diffusion terms, respectively. In the following parts, superscripts
$s$
,
$s+({1}/{3})$
,
$s+({2}/{3})$
and
$s+1$
will be used to denote the advancement from time step
$s$
to
$s+1$
, corresponding to the substeps, respectively.
In the present work without vortex stretching, the convective term in (A4) is advanced from
$s$
to
$s+({1}/{3})$
using two microsteps, by further introducing an instant
$s+({1}/{6})$
. The first microstep indicates the influence of the term
$\lambda (\unicode{x1D6E5}^{-1}\boldsymbol{\nabla }\times ((\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\omega }) )$
, written as
\begin{equation} \frac {\Delta \boldsymbol{u}^{s+\frac {1}{6}}-\sum _{q=0}^{J_e-1} \alpha _q \Delta \boldsymbol{u}^{s-q}}{\mathrm{d}t}=\sum _{q=0}^{J_e-1} \gamma _q \left [\lambda (\boldsymbol{\nabla }\times (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\omega })\right ]^{s-q}, \end{equation}
with no-slip boundary conditions
\begin{equation} \left \{\begin{aligned} &\!u_{y=\pm 1}=v_{y=\pm 1}=w_{y=\pm 1}=0 \end{aligned}. \right . \end{equation}
The second microstep indicates the influence of the term
$(1-\lambda )(\boldsymbol{u}\times \boldsymbol{\omega })$
\begin{equation} \frac {\boldsymbol{u}^{s+\frac {1}{3}}-\boldsymbol{u}^{s+\frac {1}{6}}}{\mathrm{d}t}=\sum _{q=0}^{J_e-1} \gamma _q \left [(1-\lambda )(\boldsymbol{u} \times \boldsymbol{\omega })\right ]^{s-q}, \end{equation}
with additional boundary conditions
\begin{equation} \frac {\boldsymbol{u}^{s+\frac {1}{3}}-\sum _{q=0}^{J_e-1} \alpha _q \boldsymbol{u}^{s-q}}{\mathrm{d}t}=\sum _{q=0}^{J_e-1} \gamma _q \left [(1-\lambda )(\boldsymbol{u} \times \boldsymbol{\omega })\right ]^{s-q}. \end{equation}
The calculation of the pressure term in (A4) (from
$s+({1}/{3})$
to
$s+ ({2}/{3})$
) advances as
\begin{equation} \left \{ \begin{aligned} \frac {\boldsymbol{u}^{s+\frac {2}{3}}-\boldsymbol{u}^{s+\frac {1}{3}}}{\mathrm{d}t}&=-\boldsymbol{\nabla }\varPi ^{s+\frac {2}{3}},\\ \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{s+\frac {2}{3}}&=0 ,\end{aligned} \right . \end{equation}
for the internal field, with the wall boundary condition for the pressure step modified as
\begin{align} \big(\boldsymbol{\nabla }\varPi ^{s+\frac {2}{3}}\big) \boldsymbol{\cdot }\boldsymbol{e_y}= \left (\sum _{q=0}^{J_e-1} \gamma _q \left [(1-\lambda )(\boldsymbol{u} \times \boldsymbol{\omega })\right ]^{s-q} +\nu \sum _{q=0}^{J_e-1} \gamma _q\left (-\boldsymbol{\nabla }\times \left (\boldsymbol{\nabla }\times \boldsymbol{u}^{s-q}\right )\right )\right ) \boldsymbol{\cdot }\boldsymbol{e_y}, \end{align}
where (A10) involves the unit normal vector
$\boldsymbol{e_y}$
. The coefficients in (A5), (A7), (A6), (A8) and (A10) remain consistent with conventional Navier–Stokes turbulence in the case
$\lambda =0$
.
The viscous term, advancing from time step
$s+({2}/{3})$
to
$s+1$
, is
\begin{equation} \frac {\gamma _3 \boldsymbol{u}^{s+1}-\boldsymbol{u}^{s+\frac {2}{3}}}{\mathrm{d}t}=\nu \nabla^{2} \boldsymbol{u}^{s+1}. \end{equation}
The coefficients in the semi-implicit scheme
$\alpha _0, \alpha _1, \alpha _2, \gamma _0, \gamma _1, \gamma _2$
and
$\gamma _3$
are reported in Karniadakis, Israeli & Orszag (Reference Karniadakis, Israeli and Orszag1991).
Open access
