A small amount of polymer additives can cause substantial reduction in the energy dissipation and friction loss of turbulent flow. The problem of polymer-induced drag reduction has attracted continuous attention over the seven decades since its discovery. However, changes in research paradigm and perspectives have triggered a wave of new advancements in the past decade. This review attempts to bring researchers of all levels, from beginners to experts, to the forefront of this area. It starts with a comprehensive coverage of fundamental knowledge and classical findings and theories. It then highlights several recent developments that bring fresh insights into long-standing problems. Open questions and ongoing debates are also discussed.

## I. INTRODUCTION

In 1883, Reynolds reported his historical experiments on the transition to turbulence in pipe flow:^{1} as fluid velocity increases, its trajectories in the flow change from a steady linear pattern (laminar flow) to a dynamic and sinuous one (turbulence). Dimensional analysis leads to the conclusion that for the same type of flow setup, the transition, as well as the whole flow behaviors, would depend solely on one nondimensional parameter, later known as the Reynolds number,

where *ρ* and *η* are the fluid density and viscosity (to be specific, dynamic viscosity as against the kinematic viscosity *η*/*ρ*), respectively, and *U* and *l* are the characteristic velocity and length scales of the flow, respectively. The laminar-turbulent (L-T) transition is marked by sharp changes in flow statistics. Most notably, the friction factor rises abruptly. In the turbulent regime, friction loss is significantly higher than that of a laminar flow when compared at the same Re. Techniques for friction drag reduction (DR) are thus of significant practical interest.^{2,3}

Polymer additives have long been known as highly potent drag-reducing agents. In 1948, Toms^{4} reported that dissolving a minute amount [*O*(10) wppm—parts per million by weight] of poly(methyl methacrylate) (PMMA) into monochlorobenzene can substantially reduce the friction drag, compared with that of the pure solvent, in a high-Re pipe flow. Similar observations were subsequently made in a wide variety of polymer-solvent pairs and, under certain circumstances, the percentage drop in friction drag can be as high as 80%.^{5,6} The ubiquity of DR across different chemical species shows the effect to be purely mechanical, caused by the coupling between polymer dynamics and turbulent flow motions rather than any specific chemical interaction. The most effective drag-reducing polymer molecules are linear long chains with flexible backbones,^{5–8} although rigid polymers are also known to cause DR.^{9,10}

The phenomenon of polymer-induced turbulent DR has been continuously studied for seven decades since its discovery, and findings were extensively reviewed in the literature. The classical review by Virk^{6} is a comprehensive and well-organized account of major experimental observations up to that time, especially flow statistics and their parameter dependence. Later availability of new research tools, in particular, direct numerical simulation (DNS)^{11} and particle image velocimetry (PIV),^{7,12} allowed the access to detailed flow and polymer stress fields, which has led to significant new discoveries in the past two decades. Many of those advances were covered in more recent reviews by Graham^{13} and by White and Mungal.^{14}

Despite a long history of research, this area has witnessed a wave of recent advances that pushed the boundaries of our knowledge. These developments, which mostly occurred over the past 10 years, were largely triggered by the shift of focus from ensemble flow statistics of turbulence to its dynamical heterogeneity, intermittency, and transitions between different flow states. Significant breakthroughs have been made in areas of long-standing difficulty, which is the primary motivation for this review. These newest developments will be covered in Sec. III. Established phenomenology and classical theories will still be reviewed in Sec. II to present a self-contained overview. Note that separation between Secs. II and III is not strictly chronological. A number of very recent and interesting contributions are covered in Sec. II for their better conceptual alignment with the more established framework of study.

Finally, as is the case in any other review, the current work is inevitably limited by the author’s scope of knowledge and expertise. Although every effort has been made to stay neutral and objective when describing previous findings and observations, interpretation and discussion naturally reflects my own opinion. The reader is advised to consult a number of earlier reviews for more balanced viewpoints.^{5,6,13–17}

## II. FUNDAMENTALS

This section starts with an introduction to the basic concepts, terminology, and theoretical foundation in polymer DR (Sec. II A). I attempt to take a pedagogical approach and write for the broadest readership possible, which I believe is necessary for the interdisciplinarity of this area. It is indeed an unusual marriage between two otherwise distant fields of physics, i.e., flow turbulence and polymer dynamics. Consequently, it is rather rare for a beginner to have prior exposure to both fields. After that, Sec. II B will summarize major transitions between flow regimes and phenomenological observations in each regime. Classical understanding and theories will be reviewed and discussed in Sec. II C.

### A. Basic concepts and background knowledge

#### 1. Flow geometry and setup

The most studied flow geometries are pipe and channel (Fig. 1) flows where the flow is driven by a pressure difference between the inlet and the outlet. The flow can be set up either with a fixed average pressure gradient (where the flow rate is allowed to vary) or with a fixed bulk average velocity *U*_{avg} (where the pressure drop is allowed to vary). (Hereinafter, bulk average refers to an average over the volume of the flow domain.) The latter seems to be more common in experiments while the former is more often seen in simulations, although with exceptions. When the bulk velocity is the controlled variable, it becomes the natural choice of the characteristic velocity *U* in the Re definition [Eq. (1)]. When the average pressure gradient is held constant, time-averaged velocity is available only *ex post facto*. Otherwise, a velocity scale can be deduced from the applied pressure gradient, e.g., using the laminar flow velocity generated from the same pressure gradient. For the characteristic length *l*, the pipe diameter or half (occasionally, full) channel gap height is the common choice.

Other flow setups often used include plane Couette (flow between two parallel plates driven by their relative velocity instead of pressure gradient),^{18} boundary layer (flow generated near the solid surface when a uniform bulk flow is moving over a stationary plate),^{19} and duct (pressure-driven flow in a straight conduit with uniform square or rectangular cross sections) flows.^{20}

#### 2. Definition of drag reduction

Friction loss is measured by using the Fanning friction factor *C*_{f} defined as^{21}

where *τ*_{w} is the average wall shear stress [subscript “w” indicates quantities at the wall(s)] and $(1/2)\rho Uavg2$ is the average fluid kinetic energy per volume. At steady states,

is a balance of the flow’s driving force (left-hand side—LHS), provided by the average pressure drop

and resistance due to wall friction (right-hand side—RHS). Areas of the flow cross section and the walls, *A*_{cross section} and *A*_{wall}, can be calculated from the flow geometry, which, combined with Eqs. (2) and (3), leads to the specific forms of the definition in pipe and channel flows as follows:

(*D* and *L* are the pipe diameter and length, respectively) and

(*H* and *L* are the channel half height and length, respectively). Note that Δ*P*/*L* is the average pressure gradient.

The level of drag reduction is quantified by

where *C*_{f,N} is the friction factor of a reference flow of a benchmark Newtonian fluid (“N” stands for “Newtonian”). The definition is straightforward when the solution is dilute enough that its shear viscosity is still constant and is indistinguishable from that of the pure solvent. Turbulent flow of the solvent, which presumably is a Newtonian fluid, in the same flow geometry and with the same pressure drop or flow rate, is a natural choice of the reference flow. When compared at the same flow rate (same *U*_{avg}), DR is reflected in the decrease in pressure drop. Likewise, when compared at the same Δ*P*, DR is reflected in the increase in flow rate (Fig. 1).

At higher polymer concentrations, the solution viscosity is no longer constant and decreases with shear rate $\gamma \u0307$, which is known as the shear-thinning effect. Its value at the zero-shear limit,

is usually considerably higher than the solvent viscosity. In such cases, it is common to define the benchmark fluid as a hypothetical Newtonian liquid having the same viscosity as the polymer solution shear viscosity $\eta (\gamma \u0307w)$, corresponding to the average shear rate measured at the walls $\gamma \u0307w$ (see, e.g., the work of Warholic, Massah, and Hanratty,^{22} and Ptasinski *et al.*^{23}). This is to ensure that any reduction of the friction factor caused by shear thinning, i.e., reduction of viscous shear stress resulting from the lowering viscosity, is not considered. Rather, DR% should only include the reduction of turbulent friction loss. Lumley^{5} insisted a more stringent criterion which uses the pure solvent, whose viscosity is lower than or at best equal to that of the polymer solution, as the benchmark fluid.

#### 3. Turbulent inner layer: Scales and flow characteristics

Flow statistics of wall turbulence are often reported in nondimensionalized quantities using the so-called viscous scales or turbulent inner scales.^{24} The practice is very common in the DR literature but can be confusing to beginners. At its heart, it is rooted in the separation of scales between turbulence in the bulk flow and that in near-wall layers.

For Newtonian flow, Re [Eq. (1)] is the only dimensionless parameter in the governing equations if *U* and *l*, sometimes called the outer scales, are used to nondimensionalize all flow quantities. When the flow is turbulent, boundary layers develop near the walls, within which it is expected that the most relevant scales are not those of the bulk flow (*U* and *l*), but scales derived from flow quantities at the wall, with wall shear stress

where

being the only choice available. Note that in this paper, spatial axis labels are assigned according to the commonly accepted convention in wall turbulence (also marked in Fig. 1):

*x*: streamwise direction, aligned with the mean flow;*y*: wall-normal direction, perpendicular to the wall; and*z*: spanwise direction, aligned with the vorticity of the laminar shear flow.

Here, *y*_{w} is the position of the wall, *U*_{m} is the mean velocity profile

and ⟨·⟩ indicates ensemble average, which, in wall turbulence, is averaged over homogeneous directions, *x* and *z*, and time *t*.

Combining *τ*_{w} with fluid properties *ρ* and *η*, the only velocity and length scales to be derived, i.e., the turbulent inner scales, are the friction velocity

and viscous length scale

The latter is also colloquially known as the “wall unit.” The friction Reynolds number

is defined based on friction velocity and the flow geometric length scale *l* (instead of *l*_{v}, which would give the same trivial value of 1 for all cases). Comparing Eqs. (13) and (14), it is clear that

i.e., the number of wall units in the flow domain size *l* (see Fig. 1). Combining Eqs. (12) and (14), we get $Re\tau =\rho \tau wl/\eta $. For flow under the constant-Δ*P* constraint, *τ*_{w} is constant [see Eq. (3)], and thus, Re_{τ} is constant and predetermined.

Quantities nondimensionalized with inner scales are marked with the superscript “+,” e.g.,

When discussing near-wall turbulence using the inner scales, it is also customary to use the wall coordinate in the wall-normal direction,

which measures the distance from the wall in wall units (Fig. 1).

At sufficiently high Re, many flow quantities scale with the inner scales up to |*y* − *y*_{w}| ∼ *O*(0.1)*l*. These regions are collectively called the turbulent inner layer and, correspondingly, the bulk of turbulence beyond the near-wall region is called the outer layer. The best known example of the success of the inner scaling is the von Kármán law of wall,^{24} which shows that for different flow geometries and for vastly different Re, the near-wall mean velocity profile of Newtonian turbulence follows a universal logarithmic relation when scaled in inner units:

where *A*^{+} and *B*^{+} are constants believed to be universal among Newtonian wall flows, with some small variations between different sources. Pope^{24} found that literature values are within 5% variation of *A*^{+} = 2.44 and *B*^{+} = 5.2. DNS data from the work of Kim, Moin, and Moser^{25} were best fit by *A*^{+} = 2.5 and *B*^{+} = 5.5. The logarithmic profile is found to hold over most of the near-wall region from *y*^{+} ≈ 30 up to |*y* − *y*_{w}| ∼ 0.3*l*, which is often called the “log-law layer.” In Fig. 2, a clear logarithmic dependence is observed in the water case at *y*^{+} ≳ 30. The polymeric cases will be discussed in Sec. II B 2.

In regions closest to the wall, viscous dissipation dominates and the mean velocity follows,

This region is called the “viscous sublayer,” and Eq. (19) is sketched in Fig. 2 as a curved solid line. Between *y*^{+} ≈ 5 and *y*^{+} ≈ 30 is a transition region called the “buffer layer.” In addition to flow statistics, inner scales are also found to describe flow structures well: e.g., the average spacing between near-wall velocity streaks (again, for Newtonian flow) is found to be around 100 wall units for various flow conditions.^{26–28}

#### 4. Experimental systems and parameters

Experimental parameters include the flow geometric size (pipe diameter or channel height), flow rate, and polymer solution properties. The latter include the polymer and solvent species, polymer molecular weight, and concentration. Although a variety of polymer-solution pairs have been tested especially in earlier years, after it was established that DR depends only on the mechanics of polymer molecules, aqueous solutions of flexible hydrophilic polymers, such as poly(ethylene oxide) (PEO) and polyacrylamide (PAM), have become the most widely used systems for fundamental research.^{5,6,8,29} More rigid bio-based polymers, in particular polysaccharides such as scleroglucan and xanthan gum (XG), are also often studied.^{30–33} Typical concentrations are in the range of *O*(10)–*O*(100) wppm. For most drag-reducing polymers, *O*(100) wppm is sufficient to observe a substantial increase in the zero-shear viscosity and clear shear-thinning. A shear viscosity vs shear rate curve needs to be measured for the proper calculation of DR%, which, as discussed in Sec. II A 2, requires the $\eta (\gamma \u0307w)$ value for each flow condition. Most of those molecules have very high molecular weights, *O*(10^{6}) Da, corresponding to over *O*(10^{4}) repeating units, or higher.

#### 5. DNS: Constitutive models and physical basis

##### a. Introduction: DNS and constitutive equations.

Ever since its first successful implementation by Sureshkumar, Beris, and Handler,^{11} DNS has become an indispensable tool in DR research which complements experiments with fully detailed and time-resolved flow field data as well as direct access to polymer conformation and stress fields. In DNS, the time-dependent Navier-Stokes (N-S) equation (momentum and mass balances) is numerically solved and the polymer stress contribution to momentum balance is modeled with a viscoelastic constitutive equation. Sureshkumar, Beris, and Handler^{11} adopted the FENE-P constitutive model (FENE-P stands for *f*initely *e*xtensible *n*onlinear *e*lastic dumbbell model with the *P*eterlin approximation),^{36} whereas Oldroyd-B and Giesekus models have also been occasionally used in later studies.^{37–40}

All three models treat polymer molecules as elastic dumbbells—two beads or force points connected by a spring force [Fig. 3(b)]. This effectively reduces the chain conformation into its end-to-end vector ** Q**, which describes the orientation and extension of the chain and ignores all internal degrees of freedom. (In this paper, boldface symbols indicate vectors or tensors.) Both Oldroyd-B and Giesekus models use the Hooke’s law for the spring force: i.e., the force is proportional to the chain extension |

*Q*|. (Hereinafter, |·| denotes the

*L*

^{2}-norm of a vector, i.e.,

.) This is clearly unrealistic at the limit of large deformation since polymer extension is bounded by its contour length, whereas a Hookean spring can be stretched infinitely. In a FENE dumbbell, the spring force depends nonlinearly on |*Q*|,

where *H*_{s} is the spring force constant and *Q*_{max} is the upper limit of |*Q*|: *Q*_{max} equals the chain contour length. It is easily verifiable that the force diverges as |*Q*| → *Q*_{max}, which ensures the upper-boundedness of |*Q*| [Fig. 3(b)].

The mechanical model behind FENE-P is clearly physically sounder and most closely represents the drag-reducing fluids of concern: i.e., dilute solutions of flexible polymers. Indeed, without the finite extensibility constraint, the Oldroyd-B model fails to predict shear-thinning and erroneously predicts infinite extensional viscosity at a finite extension rate. The Giesekus model incorporates interactions between Hookean dumbbells, making it a natural fit for more concentrated solutions. As a result, FENE-P has been by far the most widely adopted model in the DNS of viscoelastic turbulence. Meanwhile, limited comparisons between constitutive models found in the literature did not show any significant difference in their physical results.^{37,38} The following discussion will thus focus on the FENE-P model.

##### b. FENE-P: Basic concepts and assumptions.

Here, we briefly go over a few basic concepts in polymer physics in relation to dilute solutions,^{34} which will be mentioned repeatedly in this review. It is followed by a very brief introduction to the physical basis of the FENE-P model and its associated approximations.^{35,36}

**Ideal chain and solvent conditions**. An ideal chain model neglects nonbonded interactions, such as van der Waals (vdW) and electrostatic interactions, between chain segments.^{34} Obviously, real chain segments at least experience close-range steric repulsion and longer range vdW attraction. The latter is tunable through solvent effects. For a given polymer-solvent pair, the temperature at which the repulsion and attraction balance each other is called the θ-temperature *T*_{θ}. Polymer chains under this θ-solvent condition effectively behave like ideal chains. At *T* > *T*_{θ}, known as the good-solvent condition, the polymer-solvent interaction is more favorable, allowing polymer chains to swell and become more exposed to solvent molecules. This can be described as an *effective* reduction in the intersegmental attraction. Likewise, at *T* < *T*_{θ}, polymer chains contract, which is called the poor-solvent condition.

**Dilute solution**. A polymer solution is considered dilute when individual polymer chains are so far apart from one another that their interactions can be neglected. At equilibrium, polymer chains take a random coil conformation [see Fig. 3(a)], which is statistically most probable—i.e., having highest entropy. At the infinitely dilute limit, each coil is effectively isolated from the rest. With increasing polymer concentration *C*_{p}, the average distance between coils shortens. Upon a critical concentration level $Cp*$, this distance becomes shorter than the coil diameter and different coils start to penetrate into the space occupied by one another. This so-called overlapping concentration is traditionally regarded as the upper bound of the dilute regime.^{34} Recent experiments showed that, in both shear and extensional flows, the polymer relaxation time *λ*_{H}—i.e., the time scale for an extended polymer chain to retract to its equilibrium coil conformation—starts to increase with *C*_{p} well before $Cp*$ is reached.^{41,42} For polystyrene in both θ- and good-solvent conditions, concentration-dependence was found at *C*_{p} as low as $O(10\u22122)Cp*$. For a molecular weight of ≈7 × 10^{6} Da studied in the work of Del Giudice, Haward, and Shen,^{42} $Cp*$ was estimated at ≈5 × 10^{3} wppm in a θ-solvent and even lower in a good solvent. This indicates that under flow conditions, considerable interchain interactions may start to exist at *C*_{p} ≲ 100 wppm, which covers most DR systems. One plausible explanation is that, as polymer molecules are extended by the flow, each chain would span a larger volume than an equilibrium coil.^{42}

**Physical basis and assumptions of FENE-P**. FENE-P can be derived for polymer chains from a molecular mechanics approach by considering the free energy of stretching polymer chains in dilute solutions, with the following few assumptions and approximations:^{34–36}

The free energy of chain extension derived assuming ideal chain conformation leads to a force law

(*f*) in the form of an inverse Langevin function.*Q*The inverse Langevin force law can be approximated by a mathematically simpler empirical form of Eq. (21).

A mathematical Peterlin approximation is introduced to allow the closure of the constitutive model.

As such, FENE-P is an idealized and approximate model for drag-reducing polymer solutions. The validity of and errors from these simplifications depend on the specific polymer type and solvent conditions. There has not been much effort among DNS researchers in establishing the connection between FENE-P and realistic drag-reducing polymer solutions, which is one of the areas open for future research (see the discussion in Sec. IV). Existing DNS studies, in a way, can thus be viewed as describing the DR behaviors of a generic type of “FENE-P polymer” solutions.

##### c. Polymer conformation tensor.

Numerical solutions from DNS contain the time-resolved three-dimensional (3D) turbulent velocity $v$(** r**,

*t*) and pressure

*p*(

**,**

*r**t*) fields (

**denotes 3D spatial coordinates and**

*r**t*is time). Constitutive models such as FENE-P solve for the nondimensionalized polymer conformation tensor

field, from which the polymer stress contribution *τ*_{p} can be calculated. [Here, *k*_{B} is the Boltzmann constant, *T* is the absolute temperature, and ⟨·⟩, again, denotes an ensemble average—in the case of Eq. (22), it is the average over all individual polymer molecules at the local flow region of (** r**,

*t*). Symbol

**is sometimes used instead of**

*C***in the literature.] Note that the**

*α**xx*-,

*yy*-, and

*zz*-components of the ⟨

**⟩ tensor are $\u27e8Qx2\u27e9$, $\u27e8Qy2\u27e9$, and $\u27e8Qz2\u27e9$, respectively. Thus, the trace of the***QQ*

**tensor**

*α*is the nondimensional mean square end-to-end distance of polymer chains and a common measure of chain extension.

##### d. DNS parameters based on the FENE-P constitutive model.

Dimensional analysis of the governing equations leads to four nondimensional parameters (compared with only one in the case of Newtonian fluids) that fully define the system.

**Reynolds number**. Re is the same as that in Eq. (1), except that *ρ* and *η* are now the density and zero-shear viscosity *η*_{0} [Eq. (8)] of the polymer solution.

**Viscosity ratio**. *β* is defined as the ratio of the solvent viscosity to that of the solution

(subscripts “s” and “p” represent solvent and polymer contributions, respectively). Of course, *β* = 1 is the pure Newtonian fluid limit and *β* decreases with increasing polymer concentration. Typical DNS studies use *β* ≥ 0.9.^{11,43}

**Weissenberg number**. Wi is defined as the product of polymer relaxation time *λ*_{H} and a characteristic shear rate of the flow $\gamma \u0307$,

in which, for elastic dumbbell models,^{36}

and *ζ* is the friction coefficient of a single bead in the dumbbell: i.e., each bead experiences a viscous drag force from the surrounding solvent

which is proportional (and opposite) to its relative velocity—defined as the difference between bead velocity $v$_{b} and the solvent velocity at the bead location $v$(*r*_{b}). Note that $\gamma \u0307\u22121$ has the unit of time: Eq. (25) can be interpreted as

Higher Wi indicates slower relaxation compared with the flow deformation scale, resulting in stronger “memory effect” in the fluid and thus stronger elasticity. A purely viscous fluid would respond instantly to flow deformation and thus has Wi = 0. In DNS, the mean shear rate measured at the wall $\gamma \u0307w$ is the most common choice for the characteristic shear rate. With the advancement of numerical methods, Wi up to *O*(100) has been reported in latest studies.^{44,45}

**Finite extensibility parameter**. *b*, sometimes also denoted as *L*^{2},

enforces the finite extensibility of the FENE dumbbells—comparing Eqs. (23) and (29), it is clear that

On the other hand, since the equilibrium solution of the FENE-P equation is

where ** δ** is the identity tensor—i.e.,

the mean square end-to-end distance at equilibrium is

and the last “≈” relation is because, for flexible polymers, *b* ≫ *O*(1). Comparing Eqs. (33) and (29), one can write

i.e., $b/3$ is the ratio between the polymer extension when fully stretched and that at equilibrium. Earlier DNS studies have used *b* as low as 100, but *b* = *O*(10^{3})–*O*(10^{4}) are commonly seen in the later literature.^{11,43,44,46}

##### e. Connection with experimental parameters.

Let us now examine the relationship between model and experimental parameters for polymer solutions. Flow parameters are straightforward and thus not discussed. Note that in the definition of Wi [Eq. (25)], $\gamma \u0307$ is a flow parameter and only *λ*_{H} depends on the polymer solution used. The following discussion covers the effects of polymer solution properties—polymer and solvent species, polymer concentration, and molecular weight—on *β*, *λ*_{H}, and *b*.

**Effects of polymer concentration**. For a dilute solution in the strict sense (i.e., no interchain interactions), polymer concentration affects the governing equations only through the *β* parameter. Solving FENE-P for a simple shear flow, one can get the polymer contribution to viscosity at the $\gamma \u0307\u21920$ limit

(*M*_{w} is its molar mass, and $NAv$ is the Avogadro constant). Once the polymer and solvent species as well as polymer *M*_{w} are determined (i.e., *λ*_{H}, *b*, and *M*_{w} are fixed), *η*_{p} is proportional to polymer concentration *C*_{p}. Since

and for a very dilute solution (*C*_{p} → 0 limit), *η*_{p} ≪ *η*_{s}, we get 1 − *β* ≈ *η*_{p}/*η*_{s} ∝ *C*_{p}. Meanwhile, as discussed in Sec. II A 5 b, recent evidence suggested that at concentrations relevant to DR systems, interchain interactions may not be negligible at all. In this case, concentration would also affect *λ*_{H}: higher *C*_{p} → higher *λ*_{H}.

**Effects of polymer molecular weight**. Effects of *M*_{w} or, more accurately, chain length are illustrated with a scaling argument for polymer conformation.^{34,35} Without getting into the full complexity of solvent effects, the discussion here will be limited to the θ-solvent (ideal chain) condition as a simple demonstration. Effects of changing solvent conditions will be briefly discussed below but only at a qualitative level. In addition, the analysis here assumes true diluteness with no interchain interactions.

As sketched in Fig. 3(a), if we move along the contour of a flexible chain over a sufficiently long distance or arc length (long compared with the persistence length of the chain), in the absence of intersegmental interactions (ideal chain assumption), the orientation of the local segment would be decorrelated from that of the starting point. It is thus always possible to map the chain into a random walk or freely jointed chain (FJC) model in which each step covers a sufficiently large number of repeating units that directions of successive steps are uncorrelated. The step size *L*_{K} of this FJC is commonly termed the “Kuhn length.” For ideal chains, there is no energetic effect associated with changing chain conformation. Random coils are preferred at equilibrium solely because the conformational entropy would be lower with an increased chain extension. It can be shown that the “entropic force” pulling the chain ends together is equivalent to an elastic force with spring constant

which is the rationale for the elastic dumbbell model [Fig. 3(b)].

For a given polymer chemical species, the number of Kuhn steps *N*_{K} ∝ *M*_{w}. Note that *λ*_{H} is proportional to the ratio between *ζ* and *H*_{s} [Eq. (26)]. Other than *H*_{s}, *ζ* also depends on *N*_{k}, which follows the Zimm scaling

in a dilute solution with a θ-solvent. [The “∼” sign indicates that the two sides differ only by an *O*(1) constant prefactor.] Combining Eqs. (37) and (38) with Eq. (26),

For the *b* parameter, since

and, in a θ-solvent,

from Eq. (34), we have

Experimentally, parameters for the FENE-P model are obtained through fitting with rheological measurements and both scalings of $\lambda H\u221dMw3/2$ and *b* ∝ *M*_{w} have been observed under certain conditions.^{47,48}

**Effects of the solvent species**. The solvent affects the FENE-P model in two aspects. The first is the solvent-polymer interaction, which influences chain conformation. For instance, compared with the θ-solvent condition discussed above, a good solvent allows polymer molecules to expand at equilibrium, which, according to the Flory theory,^{34} changes the scaling of Eq. (41) to

and thus, the $Qmax2/\u27e8Qeq2\u27e9$ ratio must be reduced—so does *b*. Expansion of the polymer coil also changes the elastic force and thus *λ*_{H}. The second is solvent viscosity *η*_{s} which directly controls the friction coefficient *ζ* [Eq. (38)] and thus *λ*_{H}.

**Effects of the polymer species**. This is, of course, also included in the solvent-polymer interaction factor discussed above. In addition, changing chain mechanics directly affects the Kuhn length *L*_{K}. Increasing chain rigidity means more backbone carbon atoms must be represented by each Kuhn segment (higher *L*_{K}). When compared at the same contour length *Q*_{max} [Eq. (40)], it means that the total number of Kuhn segments *N*_{K} must drop. The combined results are higher *λ*_{H} and lower *b*. This is, of course, assuming that the chain is still flexible enough to be modeled by FENE-P.

##### f. Artificial diffusion (AD) in DNS.

Constitutive models for viscoelastic polymer fluids, such as FENE-P, do not have a diffusion term. Such purely convective partial differential equations are prone to numerical instability at high Wi. If pseudospectral methods are used for DNS, the only viable option for numerically stable solutions is to artificially introduce a diffusion term. In the case of FENE-P, the term

is added to the RHS of the *∂*** α**/

*∂t*dynamical equation. The Schmidt number Sc is defined as the ratio of kinematic viscosity

*η*/

*ρ*to polymer stress diffusivity—the artificial diffusivity of the transport of the

**tensor. Lower Sc corresponds to higher AD.**

*α*The practice of applying AD to viscoelastic DNS was first introduced by Sureshkumar and Beris^{49} who made the case that the numerical solution would converge to the accurate solution with mesh refinement, if the magnitude of AD is kept small and scales with the mesh size: 1/Sc ∼ *δx*^{2}/*δt* as *δx* ∼ *δt* → 0, where *δx* and *δt* are the characteristic mesh size and time step, respectively.^{11,37} A few finite-difference methods developed later can minimize and, in some cases, completely avoid the need of AD.^{39,50–54} Application of AD was prevalent in earlier years of DNS research and is still widespread among researchers. For the most part, the results represent the physical system reasonably well if the choice of Sc is carefully tested. Indeed, many meaningful insights were extracted from DNS studies using AD. The reader is referred to the work of Min, Yoo, and Choi^{50} and Yu and Kawaguchi^{39} (and to some extent, the work of Dubief *et al.*^{51}) for detailed comparative studies on the effects of AD.

The topic is brought up here because of the recent discovery of flow states where flow instabilities (not to be confused with numerical instabilities) are driven, at least in part, by polymer elasticity (which will be the focus of Sec. III B). Numerical schemes using AD are known to have difficulty with this particular class of flow states because AD smears the polymer stress field at regions with steep stress gradients, which are crucial for those elastic flow instabilities.^{55} For channel flow at Re_{τ} = 84.95, Sid, Terrapon, and Dubief^{45} reported that Sc < *O*(10) [typical pseudospectral methods require Sc ≤ *O*(0.1) for stability] would completely eradicate those flow instabilities in the solution and even higher Sc affects its accuracy. However, at a comparable Re_{τ}, Lopez, Choueiri, and Hof^{46} were able to capture those flow instabilities with Sc = 0.5 but only at a very high *b* = 40 000. Since AD is not a physically meaningful quantity, in this review, it will only be mentioned when the physical interpretation of results is expected to be affected.

#### 6. Vortex identification

Vortex is an instrumental concept for understanding turbulent dynamics. A vortex identification criterion would turn a fully resolved 3D instantaneous velocity field $v$(*x*, *y*, *z*) into a quantifiable and visualizable measure of vortex strength or intensity as well as its spatial distribution. This topic is not directly relevant to DR, but it is necessary to understand many flow visualization images from DNS.

At first glance, additional vortex identification criteria seem redundant as one would intuitively resort to the vorticity field

in which a streamwise (*x*-aligned) vortex will show as a region with large *ω*_{x} magnitude. The necessity for quantitative criteria beyond vorticity becomes clear in the example of a simple shear flow: $vx=\gamma \u0307y$ and $v$_{y} = $v$_{z} = 0, where $\omega z=\u2212\gamma \u0307$ (proportional to shear rate) even though there is no vortex at all. A vortex identification criterion must effectively differentiate swirling and shear flow motions.

The topic of vortex identification has been widely studied. Here, one of the most widely used method, the *Q*-criterion,^{56} is used as an illustrative example. From the velocity field, the rate of strain

and vorticity tensors

can be calculated, and the scalar identifier is defined as

where $\Vert \u22c5\Vert $ is the Frobenius tensor norm: e.g., $\Vert \Omega \Vert \u2261\u2211i\u2211j\Omega ij2$. In the simplest interpretation, Eq. (48) measures the difference between the magnitudes of fluid rotation $\Vert \Omega \Vert 2$ and strain $\Vert S\Vert 2$. The sign of *Q* reflects the local flow type—i.e.,

and the magnitude measures the strength of rotational or extensional motion. Vortices are defined as regions dominated by strong rotation with *Q* > *Q*_{threshold}. In flow field visualization, it is common to use the isosurfaces of *Q* = *Q*_{threshold} as a graphical representation of vortices. The choice of the threshold value *Q*_{threshold} depends on the flow field and is somewhat arbitrary. An approach for its systematic determination has been recently used in the work of Zhu *et al.*^{57} and Zhu and Xi,^{58} which was adapted from a similar approach for flow structure analysis by Lozano-Durán, Flores, and Jiménez.^{59}

There are quite a few other options for vortex identification. Nearly all of them, like the *Q*-criterion, turn the $v$ field into a scalar field that indicates the flow type and strength. For instance, the *λ*_{2}-criterion by Jeong and Hussain^{60} calculates the *λ*_{2} quantity from the velocity gradient tensor, in which *λ*_{2} < 0 corresponds to rotational flow regions. These criteria differ mathematically, but in complex turbulent flow fields, the results are practically equivalent. The reader is referred to the work of Chakraborty, Balachandar, and Adrian^{61} and Chen *et al.*^{62} for detailed comparisons.

### B. Phenomenology: Different stages of DR

The framework of DR phenomenology, in terms of the transitions between different flow regimes based on flow statistics, had been mostly established by the time Virk^{6} wrote his influential review. Meanwhile, since the late 1990s, direct access to turbulent flow structures and polymer conformation field enabled by PIV and DNS has fueled a wave of new observations, which greatly deepened our understanding of the dynamics within each regime. One key recent addition to the framework was the differentiation between low- and high-extent DR (LDR and HDR) first brought to light by the experiments of Warholic, Schmidt, and Hanratty.^{65} LDR and HDR were later shown to be two distinct flow stages driven by different DR mechanisms^{57} (see Sec. II B 2 d).

An overview of different regimes of DR behaviors in the Re-Wi parameter space is sketched in Fig. 4. Numerical simulations typically explore the parameter space along horizontal lines—i.e., increasing Wi with fixed Re, during which a series of transitions are typically observed, including the onset of DR, LDR-HDR transition, and maximum drag reduction (MDR). In relation to MDR, LDR and HDR are collectively called intermediate DR, which was indeed viewed as one homogeneous stage until the work of Warholic, Schmidt, and Hanratty.^{65} We will call the regime before the onset the preonset stage.

Experiments performed in the same flow apparatus and with the same polymer solution would see both Re and Wi increasing with flow rate [Eqs. (1) and (25)], but their ratio

would be constant. (Here, $\gamma \u0307=U/l$ is taken as the characteristic shear rate without loss of generality: for the same flow type, any choice of $\gamma \u0307\u221dU/l$.) In Fig. 4, such data points would follow an inclined line with zero-intercept (whose slope depends on polymer solution properties and flow geometric size), which is termed an experimental path by Li, Xi, and Graham^{63} and Li and Graham.^{66}

#### 1. Onset of DR

DR does not occur immediately upon introducing polymers. Rather, at a sufficiently low fluid elasticity (low Wi), flow statistics are indistinguishable from those of the Newtonian benchmark. For the given Re, there is a critical onset Wi, hereinafter denoted by Wi_{onset}, above which DR becomes discernible. The existence of this threshold is not surprising, considering that at the molecular level, polymer chains in a dilute solution are coiled at equilibrium and they unwind rather abruptly under increasing velocity gradient—the so-called “coil-stretch” (C-S) transition.^{67} In the case of FENE dumbbells in a uniaxial extensional flow, the C-S transition occurs when the Weissenberg number

reaches 1/2 ($\u03f5\u0307$ is the extension rate).^{36}

The rheological consequence of the C-S transition is substantial, including, e.g., a drastic increase in the extensional viscosity of the solution. It is thus expected that for DR, Wi_{onset} = *O*(1). For Wi defined based on the wall shear rate $\gamma \u0307w$, the onset is observed in the range of 5 ≲ Wi_{onset} ≲ 10 in DNS, with small variations between different Re^{38,43,44,57,68} and, possibly, different numerical settings.^{11,37,68} The number is higher than the expected *O*(1) magnitude because $\gamma \u0307w\u22121$ is not the time scale directly associated with DR. Indeed, dynamics within the viscous sublayer, which $\gamma \u0307w$ directly measures, is inconsequential as far as polymer-induced DR is concerned.^{69,70} Should we have a full grasp of the complex polymer-turbulence dynamics, a time scale of turbulent motion most relevant to DR, *τ*_{flow}, would be identified and the corresponding Wi = *λ*_{H}/*τ*_{flow} must be *O*(1). (Rigorously speaking, this should be called Deborah number De—see the De vs Wi discussion of Poole.^{71}) This *λ*_{H}/*τ*_{flow}, although it differs from the mean-flow definition $Wi=\lambda H\gamma \u0307w$ by one order of magnitude (estimated by comparing their onset magnitudes), is expected to be proportional to the latter.

At least in a region immediately after the onset, DR is primarily contributed by polymer effects in the buffer layer,^{6,13,14,57} where turbulence is dominated by streamwise vortices.^{28} Li, Sureshkumar, and Khomami^{44} noted that the root-mean-square (rms) streamwise vorticity fluctuations

–apostrophe denotes the fluctuating component, i.e.,

–in the buffer layer of Newtonian turbulence is $\omega x,rms\u2032=O(0.1)\gamma \u0307w$. Thus, Wi defined as *λ*_{H}*ω*_{x,rms} would be smaller than $\lambda H\gamma \u0307w$ by a factor of *O*(10). The onset threshold in the latter definition, as noted above, is 5–10, which leads to $(\lambda H\omega x,rms\u2032)onset=O(1)$. Of course, $\omega x,rms\u2032\u22121$ is only one of many plausible choices for *τ*_{flow}. How to choose this time scale and, ultimately, how to predict the onset depend on one’s interpretation of the DR mechanism, which itself is up for debate. Amore detailed discussion is deferred to Sec. II C.

#### 2. Intermediate DR: LDR vs HDR

##### a. Mean flow.

After the onset, DR increases with fluid elasticity (usually by increasing polymer concentration or molecular weight in experiments and increasing Wi or *b* in DNS). Figure 2 shows typical mean velocity profiles for various levels of DR% in channel flow experiments performed at a constant flow rate.^{22} With the bulk velocity *U*_{avg} fixed, Re_{τ} decreases with DR% from ≈1000 at the Newtonian (water) limit to ≲200 at highest DR levels. For the Newtonian (water) case, a pronounced log-law relation closely following the von Kármán asymptote [Eq. (18)] is found at *y*^{+} ≳ 30. With increasing polymer concentration, the profile is elevated. Since the profiles are scaled by the friction velocity [Eq. (12)], higher $Um+$ indicates lower *τ*_{w} and thus higher DR%. Up to DR% ≈ 30% [Fig. 2(a)], the increase in $Um+$ is caused by its higher slope in the buffer layer, which also thickens with DR%. Note that the lower limit of the log law layer increases from *y*^{+} ≈ 30 to ≈95 as DR% rises to 33%. The logarithmic relation itself still follows the same slope ($A+=ANewt.+=2.5$) but with higher intercepts (*B*^{+}) as a direct result of the velocity gain in the buffer layer. At DR% > 35% [Fig. 2(b)], the profile lifts up across the channel and the slope is substantially higher in regions where the von Kármán slope used to dominate (in Newtonian and LDR cases). Based on this clear difference in the shape of *U*_{m}(*y*^{+}), Warholic, Massah, and Hanratty^{22} divided DR into two regimes, LDR and HDR, roughly around the *DR*% ≈ 30%–35% line. This distinction was later confirmed in a large number of experimental and DNS studies of different flow geometries.^{12,23,33,43,44,52,57,72–75}

##### b. Fluctuations.

Fluctuating velocities in transverse (perpendicular to the mean flow—i.e., *y* and *z*) directions decrease with increasing DR%, and so does the Reynolds shear stress (RSS),^{11,22,23,33,43,57,73,75,76}

This is consistent with an intuitive general picture that, with DR, turbulent fluctuations are suppressed and turbulent motions are weakened, leading to less momentum redistribution to transverse directions and more efficient streamwise momentum transport.

Significance of *τ*_{R,xy} is evident once we write out the transport equation for turbulent kinetic energy (TKE),^{24,54,72}

where

is the TKE and *T*^{k} is the total TKE flux. The production rate of TKE

is generally positive and correlates directly with RSS—reduction of *τ*_{R,xy} thus suppresses turbulence generation. The $\u03f5vk$ term measures the rate at which TKE converts to heat via viscous dissipation. It is always positive ($\u2212\u03f5vk$ is negative—i.e., net loss of TKE), which is a reflection of the second law of thermodynamics. Finally, $\u03f5pk$ is the rate of TKE conversion to the elastic energy stored in stretched polymer chains, which, in theory, can be either positive (loss of TKE) or negative (gain of TKE).

Observations about streamwise velocity fluctuations $vx,rms\u2032$ are more conflicted. In the experimental $vx,rms\u2032+(y+)$ profiles of Warholic, Massah, and Hanratty,^{22} it was found that at LDR, with increasing DR%, $vx,rms\u2032+$ decreases within the viscous sublayer but increases moderately in buffer and log-law layers. The overall shape of the $vx,rms\u2032+$ profile still resembles that of Newtonian flow, which has a sharp peak near the wall. The peak position is at *y*^{+} ≈ 15 for Newtonian turbulence, and for LDR, it gradually shifts away from the wall but still stays in the buffer layer. At HDR, $vx,rms\u2032+$ decreases with increasing DR% across the entire flow domain. The overall profile takes a much flatter shape. A similar decline of $vx,rms\u2032+$ at HDR was found in DNS by Min, Choi, and Yoo^{72} and Dallas, Vassilicos, and Hewitt^{52} (see Fig. 5), but not in many other DNS studies^{23,43,57,68,73,76,77}—rather, their $vx,rms\u2032+$ continues with the same trend as LDR. Dallas, Vassilicos, and Hewitt^{52} attributed the failure of the other studies in capturing the sharp drop of $vx,rms\u2032+$ at HDR to the numerical artifact of using AD in the DNS (Sec. II A 5 f).

Remarkably, two recent experimental reports from the same group (Ghaemi and co-workers)^{33,75} using different polymers showed both types of behaviors. For PAM, which is a flexible polymer, the behavior is consistent with that found by Warholic, Massah, and Hanratty^{22} (and thus with that of Dallas, Vassilicos, and Hewitt^{52} and Min, Choi, and Yoo^{72}), but for XG, which is more rigid, the $vx,rms\u2032+(y+)$ profiles at HDR remain similar in shape as the Newtonian and LDR cases. Similar $vx,rms\u2032+(y+)$ behaviors can also be seen in XG solutions measured by Escudier, Nickson, and Poole.^{78} As discussed in Sec. II A 5 f, one known spurious effect of AD is its suppression of flow instabilities (or, loosely speaking, “turbulence”) that are driven, in part or in whole, by fluid elasticity.^{45,55} (By contrast, turbulence, in the conventional sense, is driven by fluid inertia and suppressed by polymers—see the more detailed discussion in Sec. III B.) In a way, DNS with AD can be viewed as a virtual experiment in which turbulent states that are elastic in nature are filtered. The lack of nonmonotonic $vx,rms\u2032+(y+)$ behaviors (as shown in Fig. 5) between LDR and HDR in those simulations, as well as in rigid polymer experiments, suggests that the drastic decrease and flattening of $vx,rms\u2032+(y+)$ could be associated with those elastic turbulent states, which may appear at HDR. Meanwhile, all other key features of HDR, including its characteristic behaviors of $Um+$ and $\u2212\u27e8vx\u2032+vy\u2032+\u27e9$ profiles as well as distinct flow structures (shown below), are still captured, indicating that turbulence driven by elasticity is not a necessary condition for the transition to HDR. Of course, this is so far only an educated guess—further research is needed for its validation.

Further complicating this issue, most DNS studies compared different DR% cases at the same pressure drop (thus same *τ*_{w} and Re_{τ}), whereas experiments such as the work of Warholic, Massah, and Hanratty,^{22} together with both aforementioned DNS cases reporting the $vx,rms\u2032+(y+)$ behaviors of Fig. 5,^{52,72} compared cases at the same flow rate (thus decreasing Re_{τ} with increasing DR%). This divide could also partially account for the differences in the observed $vx,rms\u2032+$ magnitudes: even for Newtonian flow, $vx,rms\u2032+$ goes down with decreasing Re_{τ}.^{79}

##### c. Flow structures.

Near wall turbulence is populated and, in many ways, sustained by flow structures with well-recognizable patterns [Fig. 6(a)]. The existence and importance of those so-called “coherent structures” are well documented in the literature.^{28,81–83} The best-known conceptual model involves streamwise vortices and velocity streaks. They are the characteristic structures in the buffer layer where TKE production is the highest.^{24,79,84,85} Those vortices align in the direction of the mean flow. Neighboring vortices often rotate in opposite directions, which creates stripes of either slower fluid near the wall being washed up or faster fluid closer to the bulk flushed down, forming low and high speed velocity streaks. (Such events are also referred to as “ejections” and “sweeps,” respectively, in the turbulence literature.^{86})

Availability of invariant solutions to the N-S equation enables the *a priori* study of isolated coherent structures that are otherwise braided into the complexity of full scale turbulence in DNS.^{87–91} They are fully nonlinear solutions corresponding to the “exact” forms of specific coherent structures, also known as exact coherent states (ECSs). The earliest-known form of ECS consists of a pair of staggered streamwise vortices separated by sinuous velocity streaks,^{89} which corresponds directly to buffer layer structures. (See Fig. 7. Note that at Re_{τ} = 44.2, the buffer-layer thickness is comparable to half channel height.)

Earlier attention on polymer DR effects was focused on the buffer layer, which was thought to be the primary region for DR. It was later known that this presumption is only valid for LDR where increase in the $Um+(y+)$ slope is limited to the buffer layer (Fig. 2). Indeed, until the study by Warholic, Massah, and Hanratty,^{22} the whole intermediate DR stage was thought to follow the LDR-type behavior. Based on this belief, Virk^{92} proposed his well-known three-layer model, in which DR only occurs in the buffer layer and both the viscous sublayer and log-law layer remain unaffected. Polymers cause turbulence to be suppressed in the buffer layer, leading to the enlargement of its thickness. Effects of polymers on buffer layer structures have thus been most extensively studied.

Earlier dye experiments were able to visualize velocity streak patterns. It was found that the average spacing between low-speed streaks increases with rising DR%, from 100 wall units in the Newtonian limit to over 200 wall units at high levels of DR.^{70,93} In DNS, the spanwise correlation length from velocity spatial autocorrelation functions was also found to increase with DR%.^{11,94} Later application of PIV allowed a more detailed view of velocity patterns, which showed that not only do the streaks grow wider, but they, especially at HDR, also extend along the flow direction for much longer distances without interruption and their contours become smoother and not as rugged.^{7,12} This was again widely confirmed in DNS^{39,57,73,80,95} [see Fig. 6(b)]. In particular, Li, Sureshkumar, and Khomami^{73} reported that the streamwise length scale of coherent structures can increase by over tenfold between LDR and HDR.

For vortices, DR is accompanied by the weakening of their strength and reduction of their density.^{43,51,57,73,80} The attenuated vortices dilate in their transverse scales, raising the wall-normal positions of vortex axis lines, where turbulent fluctuations are strongest, away from the wall. This is reflected in the outward shift of the peaks in *τ*_{R,xy}, velocity fluctuations, and vorticity profiles.^{12,22,33,43,52,57,72,75} Numerical investigation of viscoelastic ECS solutions rendered direct evidence that polymers weaken those coherent structures and eventually cause their extinction.^{63,66,96,97} Not surprisingly, DR is also observed in those solutions and the observed flow statistics (mean and fluctuating velocities) are consistent with the LDR behaviors discussed above.

The picture described so far targets the explanation of turbulence suppression and enlargement of the buffer layer, which is largely guided by the presumption of the Virk three-layer model. The model, despite its conceptual appeal (for its simplicity), is not consistent with later findings. An obvious miss is the HDR behaviors of the $Um+(y+)$ profiles. This suggests that structural insights summarized above are likely only accurate for LDR. Physical understanding of HDR requires the study of coherent structures at higher *y*^{+} which are, foreseeably, much more complex. Not much progress was made until fairly recently when new tools for the extraction and analysis of those structures have emerged.^{58,98} Details are discussed in Sec. III C.

##### d. Two stages of DR with two mechanisms.

Despite the common practice among researchers to use DR% = 30%–40% as an empirical threshold for the inception of HDR, sharp changes in flow statistics between LDR and HDR indicate that it is a qualitative transition between two stages of DR likely underlain by different mechanisms. There is thus no reason for it to be tied with any particular magnitude of DR%. In the absence of a presumed mechanism for the transition, the critical Wi_{LDR-HDR}, and thus corresponding DR%, of the transition can generally depend on Re, polymer solution properties, and flow geometry. (This opinion also seems to be shared by White, Dubief, and Klewicki.^{99}) Indeed, the LDR-HDR transition was observed in DNS at DR% as low as ≈15% at a low Re and in small flow domains.^{43}

A systematic investigation of the LDR-HDR transition using DNS was recently reported by Zhu *et al.*^{57} It concluded that the primary difference between LDR and HDR is the region, i.e., the span of the wall layer, influenced by DR effects. At LDR, DR effects are mostly contained in the buffer layer, although the layer does thicken with DR%. This stage of DR is relatively better understood, per the discussion above, and the Virk three-layer model still captures some essential elements. [The model’s depiction of buffer layer dynamics and its $Um+(y+)$ profile, however, is not accurate—see Sec. II B 3.] At HDR, DR effects are felt across nearly the entire domain, with the viscous sublayer being the only exception. This, of course, includes the key observation that the $Um+(y+)$ profile changes its shape in the (what used to be) log-law layer. In addition, DR effects are more directly measured by the RSS. In Fig. 8, $\tau R,xy+$ measured at a position within (*y*^{+} = 25) and one outside the buffer layer [*y* = 0.6*H* or *y*^{+} = 0.6Re_{τ}; recall Eq. (15)] is plotted against DR% for three different Re. The contrast is clear. Within the buffer layer [Fig. 8(a)], DR is continuous across the whole range of DR% and $\tau R,xy+$ drops consistently starting from the onset (DR% = 0). Outside the buffer layer [Fig. 8(b)], $\tau R,xy+$ nearly fully retains its Newtonian magnitude until DR% ≈ 20% where it starts to descend, indicating that turbulent dynamics at larger *y*^{+} remains minimally impacted until this threshold. Note that the transition point DR% ≈ 20% is, again, much lower than the commonly cited 30%–40% threshold for all three, albeit relatively low, Re tested. The conclusion that the RSS at higher *y*^{+} is only suppressed at HDR can be verified in the *τ*_{R,xy} profiles from various previous studies, using different experimental techniques or numerical algorithms, where the LDR profiles are suppressed only in the buffer layer and stay close to the Newtonian profile at higher *y*^{+}, but the HDR profiles are suppressed nearly anywhere except the *y*^{+} ≲ 5 region.^{12,22,52,76}

In addition to $Um+(y+)$ and $\tau R,xy+$, sharp transitions at higher *y*^{+} were also found, by the same study,^{57} in polymer shear stress and, more interestingly, energy spectra (which measure the distribution of TKE over different length scales). Changes in the latter showed that, at LDR, polymers suppress turbulent fluctuations and redistribute energy toward large scales at all *y*^{+} positions, while at HDR, there is a sharp increase in this effect at higher *y*^{+}, where the von Kármán law used to dominate. Disappearance of the inertia-dominated layer (i.e., the log-law layer—see Sec. II B 3) was also regarded by White, Dubief, and Klewicki^{99} as the sign of HDR based on a mean momentum balance analysis.

Comparing the flow structures of Newtonian/LDR [Fig. 6(a)] and HDR/MDR [Fig. 6(b)] flows, not only are the vortices attenuated in the latter case, but vortex distribution is also more localized, as the small number of remaining vortices tend to appear in conglomerates. Zhu *et al.*^{57} quantitatively analyzed the degree of localization at different regimes of DR and found that the localization starts at the LDR-HDR transition.

Sharp transitions in flow statistics and structures suggest the presence of two separate mechanisms for DR. The first sets in at the onset of DR and, by all accounts, is a general weakening of vortices whose effects are largely concentrated in the buffer layer. The second is triggered at the LDR-HDR transition, which suppresses turbulent momentum transport at higher *y*^{+} and extends DR effect across most of the domain height. The distinct changes in what used to be the log-law layer suggest an underlying shift in the coherent structure dynamics in that region. A plausible hypothesis was proposed by Zhu *et al.*^{57} and clear supporting evidence became available with a newest vortex tracking technique,^{58,98} which will be discussed in Sec. III C.

#### 3. Maximum drag reduction (MDR)

##### a. Basic observations.

Figure 2 shows DR% increases with increasing fluid elasticity (in their case by increasing polymer concentration). This trend is eventually bounded by an upper limit call the MDR asymptote—an “ultimate” flow regime where certain forms of turbulence still persist and the friction drag has converged to a level between the magnitudes of Newtonian turbulence and laminar flow. The phenomenon of MDR was first reported by Virk *et al.*^{69} It was subsequently found, contrary to intuition, that the MDR limit is universal—flows of different polymer solutions (changing polymer species, molecular weight, or concentration) or different geometric size *l* would converge to the same MDR asymptote.^{92} That is, DR% of various MDR flow states depends solely on their Re, even when they have different Wi, *β*, or *b* (using FENE-P parameters introduced in Sec. II A 5 d).

The Re-dependence can be completely wrapped into turbulent inner scales [note from, e.g., Eq. (15) that inner scales depend on Re] and the rescaled mean velocity profile in “+” units appears to be universal for MDR under different conditions. Virk and Baher^{100} and Virk^{92} used a logarithmic relation, same in form as Eq. (18) but with different constants, to fit various experimental MDR data. The resulting universal profile

has a markedly higher slope and seems to well approximate experimental profiles for most of the flow domain (i.e., no longer confined to a near-wall layer). Note that the two highest DR cases (DR% > 60%) in Fig. 2 are closely approaching the Virk profile (dotted-dashed line).

##### b. Validity of the logarithmic relation.

Equation (58) has long been seen as a gold standard for MDR, but its validity was recently challenged. Taking derivatives of both sides of Eq. (18), one obtains

which is the local slope value in the logarithmic relation. A true log law will show a region of nearly constant *A*^{+}. White, Dubief, and Klewicki^{101} calculated the log-law slopes of $Um+(y+)$ profiles from several recent experimental and DNS studies and concluded that, for HDR cases, a well-defined wall layer (*y*^{+} range) with a logarithmic relation is, in their word, “eradicated.” For cases with DR% ≥ 60% (both experimental and numerical, Re_{τ} ≲ 200), where the $Um+(y+)$ profile appears very close to the Virk MDR profile [Eq. (58)], the log-law slope *A*^{+} does not show any clear plateau, near 11.7 [as implied by Eq. (58)] or elsewhere. Rather, its magnitude rises up at low *y*^{+}, reaches its peak within 20 ≲ *y*^{+} ≲ 30 and then drops steeply. The peak magnitude is, nonetheless, comparable to (but not the same as) 11.7. Later experimental measurements in boundary layer flow by Elbing *et al.*^{74} led to largely similar results that for their MDR-like case (DR% = 64.8%), the $Um+$ profile does not show a logarithmic region anywhere near *A*^{+} = 11.7. For HDR (their DR% = 53% case), a roughly constant *A*^{+} region is found at *y*^{+} ≳ 200 with a magnitude higher than 2.5 (Newtonian level) but much lower than 11.7 (Virk level), which is consistent with the depiction of HDR by Warholic, Massah, and Hanratty^{22} (i.e., still logarithmic but higher slope). The missing logarithmic region in the work by White, Dubief, and Klewicki^{101} at HDR was likely caused by the relatively low Re analyzed there. Indeed, in a more recent study by the same authors^{99} (shown in Fig. 9), DNS channel flow data of Re_{τ} up to 1000 were included and, at HDR, a quasiflat region in the *A*^{+} profile can be spotted at *y*^{+} ≳ 200 (where 2.5 < *A*^{+} ≪ 11.7). Nevertheless, for MDR, conclusions of those studies are unanimous that a rigorously logarithmic profile with *A*^{+} anywhere close to 11.7 [Eq. (58)] is not found.

Recall that for Newtonian turbulence, the von Kármán law of wall was deduced from the following scaling argument:^{24}

- In the turbulent inner layer, the bulk flow rate and geometry are no longer relevant and the mean velocity gradient is determined by inner scales (Sec. II A 3) only as follows:[(60)$dUm+dy+=1y+A+(y+)$
*A*^{+}(*y*^{+}) is an undetermined nondimensional function universal to different flow conditions.]. At sufficiently high

*y*^{+}(well above the viscous sublayer but still within the inner layer), viscosity effects vanish and the flow is dominated by inertia—it is thus postulated that*A*^{+}is a constant [because*y*^{+}depends on viscosity through the definition of*l*_{v}—Eq. (13)].

Integration of Eq. (60) with a constant *A*^{+} leads to Eq. (18), and the numerical values of *A*^{+} and *B*^{+} are determined from fitting with experiments and DNS. Since $Um+(y+)$ profiles at MDR are nearly the same for different Re, i.e., inner scales are still very much relevant, the loss of a logarithmic layer at MDR is likely associated with the elimination of the inertia-dominated layer in near-wall turbulence.

The finding that a well-defined *A*^{+} = 11.7 logarithmic layer is not found at any *y*^{+} position and any stage of DR is in direct contrast to one of the key elements in Virk’s three-layer model, which postulates that the Virk MDR logarithmic profile [Eq. (58)] should appear even in the intermediate DR regime but only in a limited region. According to the model, before MDR, $Um+(y+)$ would follow the Virk profile within the enlarged buffer layer or, as termed by Virk, the “elastic sublayer” until its crossover to a logarithmic dependence with the von Kármán slope.^{6,92} The model not only misses the increased slope in the log law layer at HDR, but it is also clear from Fig. 9 that, at both LDR and HDR, the buffer layer profile does not follow the Virk logarithmic relation. Indeed, it appears that it is the reduction and elimination of the inertia-dominated log-law layer, rather than the generation and expansion of an elastic sublayer, that drives the transition to HDR and MDR.^{99}

##### c. Fundamental attributes of MDR.

Mechanistic understanding of MDR remains the most coveted target in this area, which is not surprising considering its mysterious nature. It is common practice among researchers to use the Virk profile [Eq. (58)] as an identifying trait of MDR: a flow state would be considered MDR if its $Um+(y+)$ profile matches the Virk profile. However, let us be reminded that Eq. (58) is based on empirical fitting to a somewhat arbitrarily chosen functional form: (1) there is no *a priori* physical argument for a logarithmic relationship other than a simple analogy to the von Kármán law [Eq. (18)]; (2) coefficients in Eq. (58) are empirically determined from experimental data. Latest examination of *A*^{+} reviewed above further shows that (1) a logarithmic relationship is probably inaccurate and (2) comparing $Um+(y+)$ profiles in semilog coordinates (such as Fig. 2) can be misleading: nonlogarithmic profiles are not sufficiently distinguished from logarithmic ones.

Therefore, research of MDR must go beyond this fixation on the Virk logarithmic profile of Eq. (58) and qualitative traits must be considered.^{57,102} In particular, any complete theory for MDR must consistently explains its three key attributes.

- Existence:
Why does there have to be an upper bound for DR at all? If polymers cause DR by suppressing turbulence, why can they not take the flow all the way to the laminar state? This question requires fundamental insights into the turbulent self-sustaining mechanism at MDR.

- Universality:
Why is the upper bound universal for different polymer solution properties and for different flow geometric dimensions? MDR occurs at the limit of high DR where polymer effects are strongest—it is counter-intuitive that this state, or at least the mean flow statistics thereat, is not affected by changing polymer properties.

- Magnitude:
Although Eq. (58) may not have used the most appropriate functional form, it is still undeniable that $Um+$ magnitudes at MDR are quantitatively close to the Virk profile. Why would DR converge to that particular magnitude? This question must be included because, as discussed below, there are possibly multiple flow states that address the previous two points and they cannot all be MDR.

A complete answer remains elusive. However, substantial progress has been made in the past decade, which will be reviewed in Sec. III.

### C. Mechanistic understanding

Most theoretical attempts at DR predated the discovery of the LDR-HDR transition.^{22} Their focus was thus on explaining why polymers cause DR and how to predict the onset of DR based on the supposed mechanism. Efforts were also made to address MDR but success has been limited. Because the LDR-HDR transition was not considered a qualitative transition involving different DR mechanisms until Zhu *et al.*,^{57} earlier theories covered in this section were all directed at the first DR mechanism: i.e., the one responsible for the onset of DR and LDR. Recent developments of new theoretical frameworks and research methodology have led to a series of mechanistic insights into HDR and MDR, which will be discussed in Sec. III.

#### 1. Polymer-turbulence interactions

At first glance, the notion of polymer addition causing reduced friction is counter-intuitive to most, as polymers are typically associated with higher viscosity. It should be clear by now that drag-reducing polymer solutions are often too dilute to have significant viscosity increase and, when they do, the definition of DR explicitly corrects for the viscosity change (Sec. II A 2). The term DR here is thus referring not to any change in the viscous shear stress but to the reduction of the extra stress attributed to turbulent motion—i.e., the Reynolds stress—as a consequence of turbulence-polymer interactions.

To cause DR, polymer molecules have to, one way or another, suppress turbulent motion. Availability of detailed polymer conformation (and, consequently, stress) field information—thanks to extensive DNS efforts in the past two decades—has provided direct evidence in this regard. Earlier experiments have shown, by injecting polymers to different wall positions, that DR becomes substantial when polymers reach the near-wall region covering the buffer layer and lower log-law layer.^{103,104} The buffer layer is also where turbulence production peaks in Newtonian flow and where significant reduction in RSS is observed when drag-reducing polymers are added. (For the latter, we now know it only applies to LDR.^{57}) Investigation of polymer effects on buffer-layer turbulence was thus the focus of earlier DNS studies.

It is now generally agreed that, within the buffer layer, polymers reduce turbulent intensity by opposing its dominant flow structure—streamwise vortices. Polymers alter fluid dynamics through an additional term in the momentum balance, which can be described as a polymer force. Multiple DNS studies have showed that, in transverse directions (i.e., in the plane of rotation of streamwise vortices), the polymer force counters velocity fluctuations. Direct inspection of flow field visualization images indicated that the effect is strongest immediately next to the vortices.^{51,94,105} Numerical computation of ECS solutions allowed those vortices to be isolated in a static form from the complex and dynamical backdrop of turbulent fluctuations. Investigation of polymer effects on ECS confirmed the same mechanism.^{66,97} A later study by Kim *et al.*^{106} constructed statistical representations of characteristic near-wall vortices through conditional sampling—extracting and averaging local flow structures that satisfy certain prescribed conditions.^{107} Their results reaffirmed that counter-rotating streamwise vortex pairs are the most significant structures in the buffer layer and the key conclusion of polymer forces counteracting those vortices, deduced previously from arbitrarily selected images in DNS, is statistically verifiable.

As streamwise vortices lessen in their intensity, they also dilate laterally, which leads to increased spacing between vortices and the lift of vortex axes away from the wall. Since velocity streaks are created between counter-rotating vortices, increased vortex separation is reflected in the enlargement of near-wall streak spacings (see Sec. II B 2 c). Meanwhile, upward shift of vortex axes is reflected in profile peaks of fluctuating quantities, such as $vx,rms\u2032$ and *τ*_{R,xy}, moving away from the wall (Fig. 5). It also directly accounts for the apparent thickening of the buffer layer.

Suppression of vortex motion reduces the strength of velocity streaks in between, as measured by their velocity fluctuation intensity. This explains the lowering RSS magnitude. To see this, note that low-speed streaks have negative $vx\u2032$ and positive $vy\u2032$, whereas high-speed streaks have positive $vx\u2032$ and negative $vy\u2032$—both contribute to higher −⟨$vx\u2032vy\u2032$⟩. Lowering RSS then contributes less to TKE production [see Eq. (57)], which leads to an overall flow with less fluctuations and more momentum retained in the mean flow. Thus far, a convincing depiction of DR mechanism, at least applicable to LDR, has arisen. The situation of HDR is different, where the dominant structures are more complex and DR is no longer solely attributable to an enlarged buffer layer. Study into this second stage of DR has been rather limited and occurred very recently, which will be discussed in Sec. III C 1.

#### 2. Classical theories: An introduction

Distilling a simple quantitative theory from the above micromechanical depiction of polymer-turbulence dynamics is not straightforward. There has been a long-standing debate between viscous and elastic interpretations of the polymer DR mechanism. Both theories are built on the common foundation of energy cascade in which turbulence is conceptualized to consist of a hierarchy of eddies of different sizes superimposed on one another in the same flow region. TKE is produced at the largest eddies, whose sizes are comparable to that of the flow domain, and “cascades” toward successively smaller eddies as larger ones erupt, until it reaches the lower end of the spectrum—the Kolmogorov scale—where the length scale is small enough for viscous dissipation to dominate. Experimental and numerical measurements have consistently shown, through energy spectra or Karhunen-Loève (KL) analysis, that small scale structures are most susceptible to suppression by polymers.^{12,57,68,108,109} Indeed, both viscous and elastic interpretations consider polymers to truncate or disrupt the energy cascade at a certain scale larger than the Kolmogorov scale.

##### a. Viscous mechanism.

The viscous theory was proposed by Lumley,^{5,15} which postulates that the energy cascade is truncated at a larger minimal length scale because of the viscosity increment caused by polymer additives. The smallest length scale in the energy cascade, i.e., the dissipative scale, of Newtonian turbulence is called the Kolmogorov length scale *l*_{d},

which depends on the fluid viscosity as well as density and *ϵ*—the rate of energy dissipation per unit mass of the fluid.^{24} Since drag-reducing polymers often have minimal impact on the shear viscosity of the solution, the higher viscosity responsible for shifting the turbulent dissipation length scale can only be interpreted as its extensional viscosity *η*^{ext}, which, for viscoelastic dilute polymer solutions, can be much higher than their shear viscosity *η*. Indeed, for drag-reducing polymer solutions, the Trouton ratio

can be as high as *O*(10^{4}).^{110,111}

The relevance of extensional viscosity is clear considering the kinematics of turbulence in the near-wall region, where strong but transient extensional motions are constantly being generated near vortices. Polymer molecules become locally aligned and highly stretched in those extension-intensive spots. Strong retractive forces of stretched polymer chains provide the resistance to extensional flow motion, which causes the suppression of vortex dynamics.

For dilute solutions of drag-reducing polymers in uniaxial extensional flow, as Wi^{ext} [Eq. (51)] increases, its Tr undergoes a sharp upsurge from the Newtonian value of 3 to several orders of magnitude higher, as a result of the abrupt C-S transition.^{48} For FENE-P, Wi^{ext} = 1/2 is the critical magnitude for this steep transition.^{36} At lower Wi, there is no appreciable change in *η*^{ext}, and thus, no DR is expected. Therefore, an inherent implication of the viscous mechanism is that the onset of DR is associated with a critical Wi independent of polymer concentration in so far as the solution is still dilute. The latter is because in a dilute solution the C-S transition of each individual polymer chain is not affected by the presence and state of other chains. For a given polymer-solvent pair and given molecular weight, the relaxation time *λ*_{H} is determined, which means that the onset of DR is determined by the flow time scale dropping below a critical value (i.e., *λ*_{H}/*τ*_{flow} higher than a critical value)—the so-called “time criterion.”^{5,15,112}

After the onset, as Wi or polymer concentration continues to increase, *η*^{ext} increases and the energy cascade is further truncated at larger scales as the dissipative scale increases. The effect will eventually be bounded—i.e., MDR is reached—when *l*_{d} becomes comparable to the largest scale of the flow which must be the geometric length scale *l*. This idea is in line with Virk’s three-layer model^{92} which predicts MDR as the limit where the enlarged buffer layer is restricted by the flow geometric size. Since this geometric constraint is independent of polymer properties, it offers a simple explanation for the universality of MDR, but the mechanism sustaining turbulence thereat remains unspecified.

##### b. Elastic mechanism.

The conceptual framework of the elastic theory for DR was constructed by de Gennes and co-worker^{113,114} through their scaling theory for polymer dynamics within the turbulent energy cascade. Sreenivasan and White^{115} further elucidated the theory and incorporated the effects of flow heterogeneity for DR prediction in pipe flow. Differences between the viscous and elastic interpretations of the DR mechanism originate from their underlying molecular assumptions. While the viscous mechanism assumes polymer chains to be highly stretched (past the C-S transition) to display a substantial extensional viscosity increase, de Gennes considered this scenario unlikely, especially immediately after the onset. He thus postulated that polymer chains are only partially stretched and the deformation is transient in nature. The level of deformation is not sufficient to cause a significant increase in extensional viscosity, but the elastic energy stored in individual stretched chains can add up. DR is expected when the elastic stress becomes comparable to or higher than Reynolds stress −*ρ*⟨$v$′$v$′⟩. (Equivalently, since this is a scaling argument, it could also be a comparison between elastic energy and TKE.)

As TKE flows toward smaller scales in the energy cascade, the corresponding time scale also decreases. There is thus a critical scale *r*^{*} above which polymers remain undeformed (i.e., *λ*_{H}/*τ*_{r*} equals a critical Wi or, rigorously, De magnitude). For *r* immediately below *r*^{*}, polymer chains are stretched but the elastic stress is small in compassion with the Reynolds stress. Polymers are progressively more stretched at smaller scales. At another critical scale *r*^{**}, the elastic stress starts to exceed the Reynolds stress and the energy cascade is presumed to be disrupted. The entire spectrum is divided into three regimes: at *r* > *r*^{*}, polymer chains are unaffected by turbulence; at *r*^{*} > *r* > *r*^{**}, polymers are stretched but not strongly enough to impact the flow; and at *r* < *r*^{**}, polymers disrupt the energy cascade.

Although *r*^{*} is determined by a time-scale criterion and thus independent of concentration, *r*^{**} is determined by stress or energy—quantities proportional to concentration. In the energy cascade, smaller length scales are associated with lower magnitudes of TKE or Reynolds stress. When either the polymer concentration is lower or polymer chains are less stretched, the accumulated elastic energy or stress from all chains is lower. It would thus intercept the turbulent cascade at a smaller *r*^{**}. If *r*^{**} < *l*_{d}—the smallest scale in the cascade—turbulence is unaffected at all scales. The onset of DR happens through increasing concentration or increasing Wi (which determines the level of stretching) until *r*^{**} is raised above *l*_{d}. It is thus obvious that a higher Wi_{onset} would be required when the polymer concentration *C*_{p} is lower. Compared with the viscous mechanism, this concentration-dependence of the onset is a distinguishing feature of the elastic mechanism.

Prediction of MDR from the elastic theory is not obvious. Sreenivasan and White^{115} fitted various experimental data points at the margin of MDR (where the flow just reaches MDR) to the theory and found that both *r*^{*} and *r*^{**} are in the order of magnitude of the flow geometric size *l*—i.e., even the largest eddies, which are least effective at stretching polymers, can generate sufficient elastic stress for themselves to be disrupted. The study, however, did not address the asymptotic convergence of the flow with further increasing elasticity.

#### 3. Discussion: Viscous vs elastic theories

The Lumley vs de Gennes debate remains one of the most important outstanding problems in this field. The gap between the fairly detailed knowledge, at least for LDR, of polymer-turbulence interactions (Sec. II C 1) and a conclusive predictive theory underscores the complexity of the DR phenomenon.

##### a. Experimental and numerical evidence.

Both theories have found supporting evidence but none seems sufficient for an unequivocal conclusion. Virk^{6} noted that for PEO solutions of different concentrations, the onset occurs at the same well-defined Re_{τ} and concluded that the wall-shear rate $\gamma \u0307w$ of the onset is independent of concentration for a given polymer (given *λ*_{H}), which is in line with the time criterion. Concentration independence was also reported by Berman^{116} for both PEO and PAM, who further adjusted polymer relaxation time *λ*_{H} by adding glycerol to change the solvent viscosity [see Eqs. (26) and (38)] and found that the onset wall shear rate $\gamma \u0307w$ is inversely proportional to *λ*_{H}: i.e., Wi_{onset} is constant. On the other hand, data showing sensitivity of the onset to concentration were also reported, most notably by Nadolink.^{117}

DNS seems naturally suited to test the concentration dependence of the onset as it avoids experimental complexities such as polymer polydispersity, concentration inhomogeneity, and degradation.^{13} So far, the only DNS study with a direct comparison between Wi_{onset} of different concentrations was the study by Xi and Graham,^{43} in which *β* = 0.97 and *β* = 0.99 (a roughly threefold difference in concentration) have nearly identical Wi_{onset}. The study used highly constrained domain size (motivation will be discussed in Sec. III A 1) and low Re (Re_{τ} = 84.85). Concentration dependence of Wi_{onset} at more realistic flow conditions has not been tested. Indeed, accurate determination of Wi_{onset} requires a number of extended simulation runs near the onset point. Long simulation time is required in each run for lower uncertainty in time-averaged flow quantities in order to discern potential changes in DR%. The overall requirement for computational resources is nontrivial.

As to the underlying molecular assumptions, DNS has consistently shown that at the onset, average tr(*α*) ≲ *O*(0.1)*b*, which is nowhere close to full extension.^{39,43,52,57,68,73} This is often quoted as evidence in support of de Gennes’s argument for an elastic mechanism. However, one must also consider the spatial variation and transient nature of turbulent flows—the extent of stretching can vary vastly between different regions. Indeed, Brownian dynamics simulation of FENE dumbbells in turbulent channel flow has shown that even at low Wi [*O*(1) based on $\gamma \u0307w$], a small portion of the chains is highly stretched in regions with strong local biaxial extension.^{118} Therefore, despite the low overall stretching near the onset, high extensional viscosity is still possible at certain spots as a result of the transient occurrence of strong local extension. In this sense, the extensional viscosity in the Lumley theory should be not interpreted as a domain average.

##### b. Further DNS evaluation.

DNS results have been interpreted in both theoretic frameworks. Sureshkumar, Beris, and Handler^{11} tested a very low finite extensibility parameter *b*, which suppresses extensional viscosity, and found no DR for Wi up to 50. This seemingly provides direct evidence in support of the viscous mechanism. However, it is unclear how a low *b* magnitude may interfere with an elastic mechanism. Indeed, the specific influence of polymer elasticity on turbulence is not even specified in the elastic theory.

A common practice in DNS is to inspect the TKE balance [Eq. (55)] and infer a mechanism based on profiles of different terms in the TKE budget. For the elastic conversion rate $\u03f5pk$ profiles, Min *et al.*,^{38} Min, Choi, and Yoo,^{72} and later Dallas, Vassilicos, and Hewitt^{52} have all observed a three-layer distribution at LDR: $\u03f5pk$ is positive (i.e., polymers suppress turbulence) in the viscous sublayer, becomes negative (i.e., elastic energy feeds turbulence) in a small region around *y*^{+} ∼ 10–15, and then turns positive again at higher *y*^{+}. The initial interpretation by those authors was that the polymer chains are stretched by the mean shear immediately near the wall and release the elastic energy back to the flow in the buffer layer. (Recent developments indicated that the second layer with negative $\u03f5pk$ is likely caused by a new state of turbulence; see Sec. III B.) At HDR and MDR, the third layer disappears and the second layer expands across the domain.

Treating the dynamical DR mechanism as a problem of 1D (*y*-direction) transport of averaged TKE budget terms is relatively simplistic. It neglects, once again, the spatial variations of flow patterns and polymer stretching. There is also no direct information on the dynamical sequence of events: e.g., the notion of polymers moving to the buffer layer after getting stretched at the wall is mostly a speculation. In addition, between negative and positive $\u03f5pk$, it is unclear which one is more relevant to DR. Therefore, the overall proposed scenario for an elastic mechanism based on TKE budget analysis is more of a presumption than a conclusion.

##### c. Relationship and distinctions.

The viscous and elastic theories are also not as diametrical as their names might suggest. It is not a debate of whether DR is caused by viscosity or elasticity. Rather, both mechanisms require a certain level of elasticity—a purely viscous Newtonian fluid would not cause DR under the viscous mechanism. One may even argue that they are partially reconcilable. Under the viscous mechanism, polymer feedback to turbulence is almost instant: polymer molecules are extensively stretched in extensional flows as they resist such deformation at the same time. During this process, TKE is converted into elastic energy stored in the stretched chains. This could also be one possible scenario for the generation of elastic energy or stress in the elastic theory since the latter does not specify the detailed mechanism for the polymer-turbulence interaction. Of course, the lack of details also allows the elastic theory to be more plastic. For instance, it may also be interpreted in terms of the “memory effect” of viscoelastic fluids—polymers store elastic energy at one place and release it at another (as in the case of the TKE budget analysis discussed above).

The fundamental division between these two theories is not on how polymers interact with turbulence, which is not explicitly specified in the elastic mechanism, but on how such interactions are felt by turbulence and cause DR. In the viscous mechanism, polymer effects are solitary and instant. As long as one local extensional flow spot is suppressed, DR, no matter how small, will occur. DR increases as more local flow regions become suppressed. Overall, this mechanism only requires localized C-S transition and thus depends solely on Wi (i.e., the time criterion). In the elastic mechanism, polymer effects are collective. DR is only observed when the total elastic energy exceeds a certain threshold. The onset thus depends on both Wi and concentration.

##### d. Explanation of MDR.

Neither theory offers a complete account for MDR. Both fall short of portraying a turbulence self-sustaining mechanism at MDR, which is essential for the question of MDR “existence” (Sec. II B 3 c). The viscous theory considers turbulence to be suppressed beyond the increased dissipative scale *l*_{d}, which leads to an enlarged viscous sublayer. At the limit of MDR, *l*_{d} becomes comparable to *l*, which would predict near complete laminarization. The elastic theory only compares the magnitudes of turbulent and elastic energies (or, equivalently, stresses) without specifying the specific dynamics. At *r* < *r*^{**}, it can be interpreted as either elastic energy quenching turbulence or replacing it as the new driving mechanism for instability. The latter possibility is raised here because of a series of recent studies of turbulent states driven by elasticity, which will be discussed in Sec. III B.

Warholic, Massah, and Hanratty^{22} observed in their channel flow experiments that *τ*_{R,xy} effectively vanishes at MDR. The deficit left behind in the shear stress balance must have been replaced by polymer shear stress. Since RSS is related to TKE production, it was inferred that at MDR, the conventional mechanism for turbulence generation is suppressed and turbulence must be sustained by elasticity. This finding is not supported by later experimental and numerical observations, all of which showed that although *τ*_{R,xy} at MDR is typically significantly lowered compared with its Newtonian magnitude, it remains finitely nonzero.^{23,52,72} Nevertheless, the magnitude of polymer shear stress exceeding that of RSS is at least consistent with the presumption of the elastic theory.

For “universality,” an explanation is readily available from the viscous mechanism in which MDR is considered a flow state where turbulence is suppressed and polymer effects are no longer important. It poses more challenge for the elastic theory since elastic energy inherently depends on polymer properties.

Finally, neither theory provides a quantitative prediction for the velocity “magnitude” at MDR. In a related development, L’vov, Procaccia, and co-workers have proposed a model for the mean velocity profile incorporating a number of empirical approximations.^{119,120} One notable approximation that allows the coupling of polymer effects into the mean flow prediction involves the definition of an “effective viscosity” $\eta peff(y)$ (which, inevitably, must vary with *y*) relating the polymer shear stress to the mean shear rate

and the assumption that it grows linearly with wall distance |*y* − *y*_{w}| outside the viscous sublayer. The model predicts DR that increases with the $\eta peff$ magnitude. Its solution at *τ*_{R,xy} = 0 gives a $Um+(y+)$ profile agreeing with the Virk logarithmic profile [Eq. (58)].

This development is sometimes characterized in the literature as a viscous theory for DR, which is not entirely accurate. The effective viscosity $\eta peff(y)$ is not a molecular viscosity in Lumley sense as it includes effects beyond molecular momentum transport. Indeed, it is an empirical model function lumping together the effects of polymer-turbulence dynamics, without specifying the nature of those interactions. It is thus more appropriate to be categorized as a phenomenological model than as a theory. Its MDR prediction builds on the assumption that turbulence is completely quenched (zero *τ*_{R,xy}), which is what the viscous theory would predict. On the other hand, polymer shear stress *τ*_{p,xy} (through the effective viscosity function) is high, which prevents the mean velocity from going back to the laminar value. The source of *τ*_{p,xy} is not specified. Thus, the model does not necessarily preclude the elastic mechanism which also predicts the polymer stress to exceed the Reynolds stress at *r* < *r*^{**}.

## III. RECENT ADVANCES

Section II is intended to provide an overview of the major phenomenology and fundamental understandings of DR that have been more established over the decades. There are evidently many outstanding gaps in those understandings, most notably in the HDR and MDR regimes, where progress has historically been limited. Several new developments and breakthroughs on these fronts have taken place over the past 10 years. These developments show significant breakaway from the established approaches in this field and, as such, bring forward new lines of thought in the mechanistic inquiry of DR. They will be reviewed in this section.

### A. A dynamical framework for DR research

In the classical approach of wall turbulence research, flow quantities are often averaged over homogeneous (streamwise and spanwise) directions and over time. The results are typically presented using profiles of mean velocity, rms velocity fluctuations, and Reynolds stress. Balances of shear stress and TKE are also commonly studied. Theoretical analysis would thus focus on turbulent transport between different wall layers and neglect both spatial intermittency and temporal intermittency. Earlier discussion of the viscous mechanism for DR has already shown some limitations of this perspective. In particular, it seems that the Lumley theory can better explain DNS results if the characteristic extensional viscosity is interpreted as the transient local magnitudes of *η*^{ext} at regions of strong deformation, rather than its ensemble average. A number of studies from the past decade have revealed significant new insights hidden in the dynamical intermittency, which prompted a wave of developments in DR research. A recent review by Graham^{16} specifically focused on intermittency, dynamical systems, and MDR, which covered many of the, of course, earlier developments in this wave.

#### 1. Minimal flow unit (MFU)

Coherent structures are spatially and temporally recurrent structural patterns in a turbulent flow field. Spatial patches showing such patterns, such as well-defined vortices and streaks, can be clearly identified in DNS and experimental flow visualization results. Those patterns do not strictly repeat themselves but share many similarities. To a first approximation, the flow field of near-wall turbulence can be viewed as an ensemble of individual structural units, each evolving through its own life cycle. Although no pairs of such units are identical at any given moment, all units are statistically equivalent.

These so-called minimal flow units (MFUs) can be numerically isolated by constraining the periodic simulation domain to the smallest size that still sustains turbulence. For Newtonian flow, Jiménez and Moin^{27} pioneered this idea and found that a MFU reflects the correct length scales of coherent structures: e.g., its spanwise domain size $Lz+\u2248100$ is comparable to the well-established characteristic near-wall streak spacing.^{26,28} At their relatively low Re_{τ} (100–200), flow structures captured in a MFU resemble that of the ECS (Fig. 7)—for a significant portion of time, the MFU contains one low-speed streak straddled by staggered streamwise vortices. Unlike ECS which is stationary, the MFU approach captures the dynamical life cycles of such structures. Remarkably, time-averaged flow statistics of a MFU also quantitatively agree with those of large-scale turbulence.

The approach was extended to viscoelastic turbulence more recently by Xi and Graham.^{43} It was found that the domain size of a MFU increases with increasing Wi and DR%, which resonates with the observation of increased streak spacings accompanying DR. At a low Re_{τ} = 84.85 and for Wi up to 29 tested (under constant *β* = 0.97 and *b* = 5000), viscoelastic MFUs are still dominated by streamwise vortices and streaks. However, the study reported strong intermittency in instantaneous wall shear rate $\gamma \u0307w$ magnitudes, which also seems to correlate with variations in turbulence strength. Especially at higher levels of DR%, instants with significantly lower wall shear rate are frequently observed, during which vortex strength, as measured by the *Q* quantity [Eq. (48)], experiences multifold decreases.

#### 2. Phenomenology: Active and hibernating turbulence

##### a. The concept.

These low-shear events are later known as “hibernating turbulence,” a term coined by Xi and Graham^{121} in their immediate follow-up study. They discovered that dynamics in a MFU undergoes clear transitions between distinct phases with strong (active turbulence) and weak (hibernating turbulence) turbulent activities (Fig. 10). By following the time series of common indicators of turbulent activities such as wall shear stress and *τ*_{R,xy}, hibernating turbulence can be identified as distinct intervals during which those quantities drop to uncharacteristically low magnitudes.

Although such events are more often observed at higher Wi and higher DR%, viscoelasticity is not a prerequisite. In Newtonian turbulence, the instantaneous wall shear stress or *τ*_{R,xy} can be seen fluctuating around its time-average value for most of the time, but occasionally [every *O*(1000)*l*/*U* on average at Re_{τ} = 84.85], the signal drops significantly as the flow enters the hibernating phase. It would bounce back to the normal level (active phase) after spending *O*(100)*l*/*U* in hibernation. Hibernating events are rather rare at the Newtonian limit where the time-average flow quantities almost entirely reflect characteristics of active turbulence. Its frequency stays low until well after the onset of DR but rises steadily after a higher critical Wi_{c}. For the only case so far where direct comparison has been made, Wi_{c} for the increase of hibernating frequency coincides with that of the LDR-HDR transition Wi_{LDR-HDR}, as observed in the work of Xi and Graham.^{43} At the highest Wi = 29 reported, the average duration of active intervals drops by one order of magnitude to *O*(100)*l*/*U*, while, interestingly, that of hibernating intervals remain nearly unchanged from the Newtonian limit.

The term “hibernating” turbulence was thus intended as a metaphoric reference to its most notable dynamical trait that the flow enters a state of weak activities for an extended time interval but will eventually “revive” by itself to normal activity. Indeed, flow fields during hibernation show drastically weaker vortex structures. At moment (iii) of Fig. 10, *Q* values are so low across the domain that no isosurface can be found at the prescribed level, whereas clear large vortices are found at active turbulence—moment (v)—using the same *Q* level. A velocity streak is still observed, but its strength is markedly lower. In comparison, streaks in active turbulence [moment (v)] induce stronger distortion on the velocity isosurface and are also clearly wavier in the flow (*x*) direction.

##### b. Proposed link to MDR.

Figure 10 shows the instantaneous mean velocity profiles for moments before (i), during (ii–iv), and after (v) a representative hibernating interval. [For instantaneous profiles, inner units—see Eqs. (12) and (13)—are defined using the wall shear stress of that moment at the nearest wall, with no average between walls or over time. They are thus marked with “*” instead of “+.”] Time average profiles for Newtonian and viscoelastic (Wi = 29) cases are also shown: the former follows closely the von Kármán log law at *y*^{+} > 30 and the latter is lifted with a raised slope. Profiles of active turbulence [(i) and (v)] imitate Newtonian and LDR profiles—they appear roughly parallel to the von-Kármán asymptote in the log-law layer, although with larger fluctuations than time-average profiles. Profiles deep in hibernation [(iii) and (iv)], strikingly, appear close to MDR. Moment (ii) is transitional and thus less raised.

Quantitative agreement with the Virk asymptote should not be generalized considering the intermittent nature of the dynamics—no two hibernating events are identical and the $Um*(y*)$ profile rises to different levels at different instances of hibernation. Nevertheless, qualitative similarity to MDR is hard to ignore. Regardless of the magnitude, the profile shape is unequivocally different between active and hibernating phases with the latter always showing a more lifted silhouette characteristic of HDR and MDR. In a later study, Xi and Graham^{122} showed that flow statistics between active and hibernating phases are statistically differentiable. Other than mean velocity, Reynolds stress (*τ*_{R,xy} and rms velocity fluctuations) profiles are significantly more suppressed at hibernation across nearly the entire domain, which is, again, similar to MDR. Finally, the distinct patterns of turbulent structures observed during hibernation, including weak vortices and straightened streaks, are also consistent with flow visualization at HDR and MDR.^{7,39,57,73,80,95}

Above all, the clearest connection is the insensitivity of hibernating turbulence to polymer properties, which is obviously reminiscent of the universality of MDR. As mentioned above, the average duration of an active interval drops for one order of magnitude from the Newtonian limit to Wi = 29. Within the same range, the average duration of a hibernating interval is nearly invariant with increasing Wi. Xi and Graham^{121} performed a numerical experiment in which polymer stress is suddenly turned off (i.e., switching to Newtonian dynamics) at different points during the cycle in Fig. 10. It was found that at moment (i), which is immediately before the system pivots toward hibernation, removing polymer stress prevents upcoming hibernation altogether. However, once hibernation has started, turning off polymer stress does not stop the flow from slipping deeper into hibernation. Transition between these two behaviors is rather sharp—there appears to be a turning point after which the flow is on its path to hibernation and polymer stress is no longer relevant. Therefore, polymers can suppress active turbulence, shorten its duration, and allow hibernation to occur more frequently, but hibernating turbulence itself is largely unaffected by polymers.

This is explained by the substantially reduced polymer extension during hibernation, as a result of its weaker turbulent intensity. Indeed, during hibernating intervals, nearly all polymer stretching can be attributed to the mean shear of the flow and polymer stretching in transverse directions, which can only come from turbulence, becomes negligible.^{122} Overall, the level of polymer extension anticorrelates with the level of DR: moments with higher DR (lower drag) show less polymer stretching. This would be counter-intuitive from the classical ensemble-average perspective of turbulence, where higher DR is associated with higher polymer stress (required for suppressing turbulence), and can only be interpreted in the dynamical framework of active-hibernating intermittency.

##### c. Further evidence.

Intermittency between active and hibernating turbulence has been reported in a number of later studies from different researchers. Pereira *et al.*^{123} reported a similar investigation of Couette flow where the key observations are all well consistent with the channel flow findings of Xi and Graham.^{121,122} Similar findings were also reported in pipe flow by Lopez, Choueiri, and Hof^{46} for their low to moderate Wi regime (their highest Wi results are in a different regime to be discussed in Sec. III B). The existence of hibernating turbulence in Newtonian flow was repeatedly confirmed,^{124–127} mostly at Re close to the L-T transition, where these states are more frequently visited. Conditionally averaged flow statistics and structures in those hibernating intervals again clearly resemble MDR observations.

Tamano, Graham, and Morinishi^{19} reported a DNS study of viscoelastic boundary layer flow, where the streamwise direction is no longer spatially homogeneous and temporal intermittency is reflected in the downstream spatial variation of flow quantities. They observed that regions with higher DR show lower polymer stretching and vice versa, which is consistent with the same anticorrelation in the temporal dynamics of Xi and Graham:^{122} in particular, note the strong similarity between Fig. 13 from the work of Tamano, Graham, and Morinishi^{19} and Fig. 30 from the work of Xi and Graham.^{122} They also performed a numerical test similar to that in the work of Xi and Graham^{128} by turning off polymer stress after a certain downstream position with low drag and found that the low drag state can survive a substantial distance further downstream without polymer stress, which reaffirms the conclusion that polymer stress is not important within the hibernating phase.

##### d. Connecting MFU dynamics with large-scale turbulence.

Following the discussion in Sec. III A 1, MFU can be mapped to realistic large-scale flows by invoking the ergodicity hypothesis of turbulence,^{129} which postulates that both temporal and spatial statistics are equivalent to ensemble statistics. In this sense, temporal intermittency in a MFU would map to spatial intermittency in a larger domain. A single structural unit can experience active and hibernating phases at different times. If a larger flow domain is viewed as a jigsaw assembled by such units, at any moment, some units will be caught at the active phase and some in hibernation. Their proportions should reflect the time share of each phase in a single unit. Of course, this would be neglecting any correlation between individual units, which is a nontrivial assumption (see Sec. III C).

This view was adopted by Kushwaha, Park, and Graham^{124} and Wang, Shekar, and Graham^{130} for Newtonian and viscoelastic flows, respectively. They divided DNS flow fields of an extended domain into active and hibernating patches based on magnitudes of local velocity gradient in the buffer layer. Conditionally averaged flow statistics of these two categories compare well with MFU results. Conditional average was also applied to large-scale experimental measurements based on local wall shear stress,^{125,126} which again confirmed the same picture. Kushwaha, Park, and Graham^{124} showed that conditional averages of active and hibernating turbulence based on temporal intermittency and spatial intermittency analysis are consistently alike.

#### 3. DR from a dynamical systems perspective

##### a. Solutions in the state space.

From a dynamical system perspective, turbulent dynamics obtained from a MFU can be viewed as a chaotic dynamical trajectory in a state space governed by numerous invariant solutions such as ECS [see Fig. 11(a) for a simple case].^{91,131,132} Other than the laminar state, none of those solutions are linearly stable. All known ECS solutions are saddle points in the state space, having both stable and unstable directions (stable and unstable manifolds). (An illustrative visualization of how invariant solutions and their unstable manifolds guide the dynamics of turbulence in a MFU was shown in Fig. 9 from the work of Gibson, Halcrow, and Cvitanotić.^{133}) For Newtonian turbulence close to the L-T transition (low Re), dynamics in the “kernel” of strong turbulence is dominated by the so-called upper-branch (UB) ECS solutions,^{133–135} while hibernating intervals appear to be intermittent escapes from the kernel.

Since hibernating turbulence (1) shows distinctly weak turbulent activities and (2) exists in both Newtonian and viscoelastic flows, for its origin, one would first look for solutions representing weak turbulent states in Newtonian flow. In canonical flow types, transition to turbulence follows a nonlinear mechanism and requires finite-amplitude disturbances. (For channel flow, a linear instability of the laminar state does occur at Re_{τ} ≈ 107.4, which leads to a new solution called the Tollmien-Schlichting or T-S wave.^{136,137} Since Re_{τ,crit} ≈ 44.7 for the L-T transition is much lower,^{138} a nonlinear mechanism is the only pathway for the transition at Re_{τ} ≲ 107.4. Without further complicating the discussion, the linear instability scenario will not be discussed which does not exist in pipe or Couette flows.) Therefore, the state space has *at least* [additional solution(s) may exist; see Sec. III B] two basins of attraction: laminar and turbulent regions. The ridge separating them can thus be considered the most marginal or weakest form of turbulence. This so-called “edge of chaos” (EoC) can be numerically computed, noting that trajectories above and below the ridge divert to opposite basins, through repeated bisections.^{139,140} For the relatively low Re studied so far, dynamics along the EoC converges to an asymptotic solution called the edge state (ES) [Fig. 11(a)], which often appears as a quasiperiodic [Fig. 11(b)] or nonperiodic oscillatory (with varying period length) trajectory itself.^{102,141} The ES has only one unstable direction,^{139} which means for an *N*-dimensional system (*N* being the number of degrees of freedom), its stable manifold spans an (*N* − 1)-dimensional subspace and forms the separating surface between the two basins—i.e., the EoC. Lower-branch (LB) ECS solutions with one-dimensional unstable manifolds have been found and are believed to be the key constituents of the ES.^{142,143}

Xi and Graham^{128} computed the ES for viscoelastic channel flow whose flow structure and mean velocity profile turn out to be remarkably similar to hibernating turbulence. Most importantly, mean flow of the ES was found to be insensitive to viscoelasticity—from Newtonian flow to Wi = 28, the profiles seamlessly overlap with one another. The collapsed profile at Re_{τ} = 84.85 happens to align with the Virk profile [Eq. (58)], although this agreement cannot be generalized to higher Re.^{102} This insensitivity, same as the case of hibernating turbulence, is again caused by the lack of polymer stretching in this weak turbulent state, where the flow kinematics is dominated by shear rather than extension or rotation.^{102}

##### b. Active-hibernating-bursting (AHB) cycles.

The proposed dynamical scenario is sketched in Fig. 11(a). For Newtonian flow, a turbulent trajectory would spend most time sampling a well-defined kernel region of strong (active) turbulence formed by UB solutions. Occasionally, it may hit certain points of exit in the kernel region and embark on excursion toward the laminar state, only to be blocked by the EoC. The trajectory will then travel along the EoC until it reaches the vicinity of the ES, where it will be pivoted toward the unstable direction of the latter and start its return journey. During the excursion, the flow spends substantial time in regions with very weak turbulence, which is reflected as a hibernating interval in DNS.

Signs of such cycles can be spotted in earlier DNS studies. For instance, Min, Choi, and Yoo^{72} obtained their viscoelastic solutions by suddenly switching on viscoelasticity from a Newtonian solution. After leaving the Newtonian kernel, all solutions (different Wi) can be seen approaching the same low-drag (drag coefficient close to but higher than the laminar magnitude) asymptotic state and staying in its vicinity for *O*(10^{2})*l*/*U*, before bouncing backward to their corresponding viscoelastic kernels. Interestingly, the simulation domain used in that study is much larger than MFUs where dynamical cycles involving active and hibernating states were later more systematically studied.^{102,121,122,128}

The connection of LB ECSs (which, as discussed, are believed to dominate the ES) to hibernation is supported by their direct comparison with conditionally averaged hibernating events in extended flow domains from both experiments and DNS.^{125,126} Furthermore, Pereira, Thompson, and Mompean^{127} and Kushwaha, Park, and Graham^{124} found that the frequency of hibernation in Newtonian channel flow increases with decreasing Re, which is consistent with the common belief that active turbulence and hibernating turbulence are organized around UB and LB ECSs, respectively—solutions on the two branches are closer, and thus, the transition frequency becomes higher—as Re decreases toward the bifurcation point. With increasing Re, hibernation becomes rare, which only becomes unmasked when Wi is higher.^{130}

Returning through the unstable manifold of the ES often leads to a strong overshoot of turbulent intensity—a bursting event—before it settles back to the active kernel [Fig. 11(a)]. The term “bursting,” although prevalent in the turbulence literature, is not unequivocally defined and interpretation varies depending on the context. Here, it refers to events marked with a strong sharp increase in TKE. These events are strongly nonequilibrium in nature (as opposed to quasi-steady-state solutions like ECS). Their importance in the turbulence self-sustaining cycle was pointed out by Jiménez *et al.*^{134} In time series showing hibernating intervals, a strong burst of TKE, as well as other quantities associated with turbulence intensity, usually follows hibernation.^{121,122} Pereira, Thompson, and Mompean^{127} showed similar observations. They made the differentiation between “strong” and “moderate” active states, which, in our terminology, roughly correspond to bursting and active phases. Kushwaha, Park, and Graham^{124} conditionally averaged instances of hibernating intervals by aligning them at their end points. The result clearly showed that, statistically, hibernating intervals are often followed by strong overshoots of wall shear stress—i.e., “bursting.”

Zhu *et al.*^{54} systematically studied the nature of bursting events and showed that trajectories following the unstable manifold of the ES *en route* to the turbulent kernel demonstrate well-defined bursting behaviors. It starts with quick intensification and lift up of near-wall vortices, which results in a sharp rise of the RSS [and thus TKE production—see Eq. (57)]. These primary vortices then undergo swift rupture, generating countless much smaller but intense vortex pieces as their debris, which is accompanied by TKE quickly reaching its peak. Viscous dissipation then stabilizes the flow and regulates the fluctuations into larger vortices typical of the kernel of active turbulence. Park, Shekar, and Graham^{143} also showed that strong bursting events in a MFU are closely aligned with the unstable manifolds of a class of LB solutions likely sitting on the ES.^{135} The scenario depicted above resonates well with the earlier KL analysis of MFU solutions by Webber, Handler, and Sirovich.^{144} They showed that at the early stage of a bursting event (they used the term “entropy event”), the representational entropy between modes reaches minimum, which indicates that TKE is contained in a few large-scale modes—i.e., the primary vortices. A quick rise in entropy soon follows, which corresponds to the later rupture of primary vortices and redistribution of TKE between a wide range of scales.

DNS trajectories of active-hibernating-bursting (AHB) cycles, projected to a two-dimensional (2D) plane, are shown in Fig. 11(b).^{102} A numerically computed ES solution is also shown, which appears as a quasiperiodic orbit located between the laminar state and turbulent kernel. Starting from the ES and with minimal disturbances [see the solid blue line in Fig. 11(b)], instabilities would grow and lead to a strong bursting event (overshoot) before the flow decays to the turbulent kernel. Afterward, occasional escapes and excursions, i.e., hibernating events, are observed, during which the flow visits regions close to the ES. Bursting is typically observed after those visits, which is reflected in increased $\tau R,xy*$, but not necessarily in $A25*$. [$A25*$ is the logarithmic slope calculated, according to Eq. (59), from instantaneous $Um*(y*)$ profiles at *y*^{*} = 25, which is a mean flow quantity not directly tied to instantaneous turbulent intensity.]

A more detailed state-space visualization can be found in Fig. 11 from the work of Park and Graham^{135} and Fig. 9 from the work of Park, Shekar, and Graham.^{143} Although those authors did not directly compute the ES, they showed that hibernating turbulence corresponds to intermittent visits to the vicinity of a class of LB solutions believed to be part of the ES, all having only one unstable dimension. Park, Shekar, and Graham^{143} also noted that although hibernation is a necessary precursor to strong bursts, not every hibernating interval is followed by clear bursting. This is not surprising considering the intermittent nature of the dynamics: not every hibernation trip reaches close proximity to the ES. The returning trajectory is only expected to follow the unstable manifold closely, which leads to strong bursting, if it gets very close to the ES during hibernation.

##### c. Theory for MDR and HDR.

We are now ready to discuss the effects of polymers on this AHB dynamical cycle from which a new theory for explaining MDR will arise. AHB cycles exist in both Newtonian and viscoelastic flows, but at high Wi, the time scale for active turbulence is severely reduced and the frequency of hibernation, whose drag is substantially lower, is raised. Although this could also increase the frequency of bursting, whose drag is even higher than active turbulence, Zhu *et al.*^{54} showed that polymers can suppress the strongest form of bursting and reroute the dynamics to avoid vortex rupture and intense fluctuations. The net outcome is that both time-averaged flow quantities and statistically dominant flow structures will be increasingly represented by hibernation, which, as discussed in Sec. III A 2, shows key characteristics of MDR. It was thus postulated that at sufficiently high Wi, active intervals will be very short so that hibernation becomes predominant. This is when the overall flow converges to MDR.

A simple yet quantitative mathematical model was derived, by Xi and Graham,^{122} based on this conceptual picture, assuming that (1) at Wi higher than a critical value Wi_{c}, polymer chains are persistently stretched in active turbulence where polymer stress grows exponentially with time; (2) finite extensibility of polymer chains is neglected; (3) the polymer stretching rate at active turbulence is independent of Wi; and (4) active turbulence is suppressed once the nondimensional polymer stress reaches a threshold value *S*_{T}. It leads to a simple expression for the duration of active intervals

which fits perfectly with DNS data at the high-Wi end. (At lower Wi, *T*_{act} is independent of Wi.) For Re_{τ} = 84.85 used in that study, the fitting yields Wi_{c} = 17.36, which is a stunning agreement with the observed value of Wi_{c} ≈ 18 (where *T*_{act} starts to decline). At the limit of Wi → ∞, the model predicts *T*_{act} ≈ 127*l*/*U*, significantly below the average duration of hibernation *T*_{hib} ≈ 200*l*/*U*, which affirms the postulation that MDR is a flow state where hibernation becomes statistically dominant.

The overall picture is depicted in Fig. 12. Intermittent AHB transition cycles occur in both Newtonian and viscoelastic flows, only that in the former they are rare and active turbulence lasts for long periods. After the onset of DR, polymers are stretched at regions of strong flow motion (e.g., local extensional spots), which leads to an overall weakening of vortex structures at active turbulence (and, presumably, during bursting as well), following the mechanism described in Sec. II C 1. The stretching is likely transient in nature and polymer stress never grows to *S*_{T}, the level required to disrupt active intervals. At Wi > Wi_{c} (>Wi_{onset}), polymer stretching in active turbulence becomes persistent and polymer stress keeps ramping up until active turbulence is quenched, which triggers early exit to hibernation. During hibernation, polymer chains recoil and polymer effects are minimal, but, because of the unstable nature of hibernation (both ES and LB solutions have unstable directions), instability will still grow, which is often followed by a bursting event before the start of a new active interval.

Further increasing Wi continues to reduce *T*_{act} and increase the statistical weight of hibernating turbulence (which remains mostly unaffected). At MDR, some form of strong turbulence (bursting or active turbulence) must still exist and appear intermittently because polymers do not stabilize hibernating turbulence. Indeed, Zhu *et al.*^{54} found that as the flow leaves the ES, polymer effects remain small during the initial growth of instability (intensification and lift up of primary vortices) and only become significant in later stages of bursting (rupture of primary vortices), which occurs ≈50 − 70*l*/*U* after the initial instability. This is the same order of magnitude as the ≈127*l*/*U* prediction of Eq. (64) at the Wi → ∞ limit. The prediction is higher by a factor of ∼2, which is not unexpected considering the assumptions behind Eq. (64), especially the neglect of finite extensibility. It is thus inferred that bursting/active phases are terminated almost immediately after polymers are stretched. Thus, the only form of strong turbulence remaining at MDR would be the budding stage of instability before polymer stress is significant.

The theory essentially links MDR to a class of fundamentally Newtonian flow states associated with the L-T transition, in which polymer elasticity is not important. This is supported by the observations of White, Dubief, and Klewicki^{99} who compared the mean velocity and RSS profiles between MDR and transitional states in Newtonian boundary layer flow (from the work of Wu and Moin^{145}) and found remarkable similarities. The main difference, as they noted, is that at MDR, such transient states are stabilized. According to the current theory, this “stabilization” is more accurately described as an intrinsically intermittent process—viscoelasticity does not make hibernation more stable *per se* but keeps driving the flow back to hibernation to make it statistically more persistent.

Compared with classical attempts, this theory takes a unique dynamical perspective and explains the transition to MDR in terms of a shifting dynamical balance between flow states. By offering a consistent explanation for two of the three key attributes of MDR, as outlined in Sec. II B 3, it represents a major step forward toward solving this mystery. For the “existence” attribute, the theory postulates that MDR is a form of turbulence dominated by weak streaks and vortices (hibernation) with intermittent short-lived eruptions of stronger activities (bursting/active turbulence). The flow does not return to the laminar state because of the insensitivity of the ES and EoS to polymer effects—polymers cannot remove the barrier between turbulent and laminar basins. For the “universality” attribute, the explanation is obvious: hibernating turbulence, which presumably dominates MDR, is a class of fundamentally Newtonian flow states. It does not rely on viscoelasticity to exist and does not seem to be affected by the latter either. Polymers only unmask those states by compressing active intervals. The question regarding the “magnitude” of MDR remains open, which we will discuss shortly below.

The transition from vortex weakening by transiently stretched polymers (first DR mechanism at Wi_{onset} < Wi < Wi_{c}) to the suppression of active turbulence by persistently stretched polymers (second DR mechanism at Wi > Wi_{c}) also offers an explanation for the LDR-HDR transition, where a second mechanism is needed to explain the sharp differences between the two regimes. So far, the coincidence between Wi_{c} and the LDR-HDR transition has only been shown for one parameter setting (Re, $\beta $, and *b* combination). More research is still required to establish this link.

#### 4. Limitations and remaining gaps

Despite its great promise, this theoretical framework based on ABH cycles still has a few gaps to fill. Developments in more recent years (after the study by Xi and Graham^{121} in 2010) have revealed encouraging new directions, which will be discussed in Secs. III B and III C. The answers, however, are far from complete.

##### a. Quantitative magnitude of MDR.

The theory falls short of offering a mean flow prediction for MDR and cannot fully explain the quantitative magnitude of the Virk profile [Eq. (58)]. Although Virk-like profiles are often observed both during hibernation and at the ES, it is not always generalizable—other instances of hibernation or ESs at other parameters may differ from the asymptote. The intrinsic difficulty again stems from the dynamical nature of the theory, in which MDR is determined by not necessarily one state but the statistical ensemble of many. Since the ES is not altered by polymer elasticity,^{102,128} it is reasonable to anticipate that, near the ES, there is a band of the state space within which turbulence is too weak to be influenced by polymers. Mean flow at MDR, based on this theory, would be determined by the average of those states weighted by their sampling frequency during hibernation.

##### b. Asymptotic convergence.

Direct numerical evidence for the dynamical theory was only obtained for limited Wi. For Xi and Graham^{121,122} with their parameter settings, it is up to Wi = 29. At higher Wi, their DNS trajectories return to the laminar state after finite time periods, which is inconsistent with the expectation of MDR as an asymptotic state of self-sustaining turbulence. The theory builds on the immutability of the EoC and ES—these solutions are required to block any laminarization attempt as they are asymptotically and dynamically approached by the flow trajectory. This immutability is postulated based on the insensitivity of ES mean flow quantities to viscoelasticity. The observed laminarization at higher Wi apparently suggests that the L-T boundary (EoC and ES) becomes somehow penetrable. Although this does not necessarily contradict the immutability of the ES mean velocity, it is at odds with the predicted asymptotic dynamics at MDR.

There are two very different, but possibly coexisting, explanations for the premature laminarization. One is the small domain size used in those MFU studies. Wang *et al.*^{146} later reported that the maximum Wi for sustained turbulence is limited by the domain size—in particular, as $Lx+$ increases from 360 to 1000, the highest turbulent Wi increases from 29 to 90 (same Re_{τ}, *β*, and *b* as those in the work of Xi and Graham^{121}). Intermittency between active and hibernating states, similar to that in MFU, was also observed in their larger domains. Nevertheless, the asymptotic behavior of MDR was still not truly captured—although their Wi = 90 case has DR% close to 60% and its $Um+(y+)$ profile is comparable to the Virk profile of Eq. (58), further increasing Wi would still cause laminarization. Indeed, in other studies with larger domain sizes ($Lx+$ up to 4000), laminarization was also observed at sufficiently high Wi.^{57,68,73} The other possible explanation emerged very recently. It was found that under similar flow parameter settings, a new class of turbulencelike instabilities, driven at least in part by polymer elasticity, exist at high Wi.^{45,147} As discussed in Sec. II A 5 f, turbulence driven by elasticity is very difficult to capture in DNS when AD is used,^{45,55} which was the case in all above-quoted studies. In this explanation, the eventual extinction of lasting turbulence, in the conventional sense, at high Wi is real, but the newly emerged “turbulent” states would keep the flow from laminarization. Further discussion is deferred to Sec. III B, after those states are introduced.

##### c. Dynamics at higher Re.

The state-space depiction illustrated in Fig. 11(a) is based on numerical invariant solutions and DNS at fairly low Re, typically Re ∼ *O*(10^{3}) or Re_{τ} ≲ *O*(10^{2}) and not very far above Re_{crit} for the L-T transition. Extension of this scenario to higher Re, where experiments are typically conducted (Re_{τ} ≳ 10^{3}), is nontrivial. The complexity of the state space is expected to grow explosively with increasing Re because of the rapid bifurcation of existing solutions and emergence of new ones.^{135,141,148} In addition, most known ECS solutions depict streak-vortex structures (Fig. 7) representative of turbulence in the buffer layer, which, at higher Re, accounts only for a very small portion of the flow domain. Structures at higher *y*^{+} are vastly more complex,^{83,149} most of which are not included in the current framework of dynamical systems. (Recently, Shekar and Graham^{150} reported the first class of ECSs resembling hairpin vortices, which are involved in log-law layer dynamics.) Understanding polymer effects on those structures is particularly relevant to the study of HDR and MDR, where, as reviewed above, DR reduction effects are felt well beyond the buffer layer.

Overall, the state-space picture of higher Re and higher *y*^{+} is rather murky at the moment. It is unclear whether any simple dynamical theory can still formed in a similar bottom-up fashion—from invariant solutions and MFU trajectories. An alternative approach is to extract and analyze flow structures *a posteriori* from the DNS data of realistic flow conditions. Recent efforts on that front will be discussed in Sec. III C 2.

##### d. Dynamics in extended flow domains.

MFU allows direct access to the temporal evolution of isolated structures but neglects correlations between different structural units. Dynamics in a MFU can be mapped to large-scale turbulence citing ergodicity and spatial-temporal equivalency, as did by Wang, Shekar, and Graham,^{130} which, however, still precludes interactions between structures. For turbulence in larger flow domains, the generation, growth, spreading, and demise of each unit structure are inevitably coupled with the dynamics of other structures nearby and, possibly, far away. For example, Lopez, Choueiri, and Hof^{46} reported that with increasing length of their simulation domain (pipe geometry), although a transition from temporal intermittency to spatial intermittency was indeed observed, the latter cannot be accurately depicted as a statistical ensemble of different states in the temporal trajectory. Rather, with sufficient domain length, highly localized turbulent patches would form, which is a clear sign for interstructural correlation. In a separate example, Zhu *et al.*^{57} found that spatial localization of vortex structures is one of the key characteristics of HDR, which again would not occur without interaction between structures. A detailed analysis of such dynamics in the complex backdrop of turbulence appears daunting, for which latest development in vortex analysis methodology, discussed in Sec. III C 2, can be particularly instrumental.

### B. Elasticity as a driving force for turbulence

Conventional turbulence, e.g., as in a Newtonian fluid, is driven by inertial effects. Adding drag-reducing polymers damps those turbulent motions and coerces the flow to a state of less disorderness and, consequently, less drag. This has been the storyline dominating DR research for decades and has also been the narrative of this review thus far. Other than inertia, in polymer solutions, fluid elasticity provides an additional source of nonlinearity in system dynamics, which, in principle, could also drive (instead of damping) flow instabilities. Such elastic instabilities are more likely to occur as Wi increases. For this reason, it has been long speculated that MDR is formed by flow instabilities that are more elastic in nature. The report of vanishing RSS (and thus large polymer shear stress) at MDR by Warholic, Massah, and Hanratty^{22} further fueled such speculations. On the other hand, DR is still a high-Re phenomenon—it is equally possible that stronger conventional turbulence driven by fluid inertia prevails over any potential elastic instability. Indeed, solid evidence for the relevance of such instabilities, driven at least in part by elasticity, to DR only emerged very recently: the first and so far only established case is the so-called elastoinertial turbulence (EIT) reported by Samanta *et al.*^{147}

#### 1. Background

##### a. Historical context of elastic instability research.

It has been widely known that viscoelastic polymer fluids can display instabilities that are driven either purely or in part by elasticity. Purely elastic instabilities can occur at the limit of vanishing Re (i.e., the inertia-less limit).^{151–155} Many earlier studies focused on polymer melts for rheometry and polymer processing applications. There has also been substantial research of dilute solutions of long-chain (i.e., drag-reducing) polymers in, e.g., microfluidics, where, because of their small geometric dimensions, Re is too low to generate turbulence and elastic instabilities are explored for the purpose of promoting fluid mixing. The best understood type of elastic instabilities occur in flows with curved streamlines, such as Taylor-Couette flow—shear flow in the annular space between two coaxial cylinders with relative rotation.^{156} The mechanism of this so-called “hoop-stress” instability involves the coupling between polymer normal stress and streamline curvature.^{157,158} Flow instabilities typically manifest as secondary and often oscillatory flow patterns, which cause increased drag compared with stable laminar flow. In a rotating parallel-plate setup, Groisman and Steinberg^{159} observed a particular type of instability showing chaotic flow patterns. This instability was shown to exist at arbitrarily low Re, indicating its purely elastic nature. The term “elastic turbulence” was coined to describe this turbulencelike instability that does not rely on inertia.

Elastic instabilities also occur in flow around a stagnation point, where the streamlines make sharp turns, such as cross-slot flow.^{160,161} One mechanism for those instabilities is the coupling between the incoming convection of polymer stress and fluctuations in the width of the so-called “birefringent” strand—a thin sheet of fluids carrying highly stretched polymers, which initiates from the stagnation point but extends far downstream.^{162}

##### b. Dual origin for flow instabilities.

Mathematically, flow instabilities stem from the nonlinear relationship between the stress and rate of strain, which is found in both fluid inertia (measured by Re) and elasticity (measured by Wi). Combining these two effects can also give rise to new types of instabilities. Such inertio-elastic instabilities (IEIs) were found in a variety of flow geometries. In a microfluidic planar contraction-expansion flow, Rodd *et al.*^{163} observed that the instability secondary flow pattern varies between different Re-Wi combinations: instabilities occurring at finite Re (i.e., IEI) are notably different from those at the purely elastic (vanishing Re) limit (both at high Wi). In a cross-slot flow, Burshtein *et al.*^{164} showed that instability around the stagnation point morphs from an inertia-dominated form to an elasticity-dominated form by raising Wi at constant finite Re.

According to their physical origin, flow instabilities in viscoelastic fluids can be categorized into the following three major types:

- elasticity-driven instability (EDI):
often called purely elastic or inertia-less instabilities in the literature, which are driven

*entirely*by elasticity and can occur at the Re → 0 limit; - inertia-driven instability (IDI):
instabilities that would also occur in a Newtonian flow at sufficiently high Re; and

- inertio-elastic instability (IEI):
instabilities where both inertia and elasticity are essential and both Re and Wi need to reach critical threshold values.

These different instabilities were systematically explored in the Re-Wi parameter space by Rodd *et al.*^{163} (see Fig. 15 in their work) and Burshtein *et al.*^{164} (see Fig. 13 in their work). A generic overview of possible flow states at different Re-Wi regimes is illustrated in Table I.

Re > Re_{crit} | Inertia-driven | IDT with DR | Elastoinertial |

turbulence (IDT) | turbulence (EIT) | ||

O(1) ≲ Re < Re_{crit} | Laminar flow/inertia-driven | Laminar flow/IDI/inertio-elastic | IEI/EIT |

inability (IDI) | instability (IEI) | ||

Re ≪ O(1) | Laminar flow | Laminar flow/elasticity-driven | EDI/elasticity-driven |

instability (EDI) | turbulence (EDT) | ||

Wi ≲ O(1) | Wi ≳ O(1) | Wi ≫ O(1) |

Re > Re_{crit} | Inertia-driven | IDT with DR | Elastoinertial |

turbulence (IDT) | turbulence (EIT) | ||

O(1) ≲ Re < Re_{crit} | Laminar flow/inertia-driven | Laminar flow/IDI/inertio-elastic | IEI/EIT |

inability (IDI) | instability (IEI) | ||

Re ≪ O(1) | Laminar flow | Laminar flow/elasticity-driven | EDI/elasticity-driven |

instability (EDI) | turbulence (EDT) | ||

Wi ≲ O(1) | Wi ≳ O(1) | Wi ≫ O(1) |

##### c. Elastic instabilities in parallel shear flows.

Canonical turbulent flow types are all parallel shear flows whose mean flow has straight streamlines. Instability mechanisms reviewed above clearly do not apply. Indeed, common viscoelastic parallel shear flows, such as Couette, channel, and pipe flows, are known to be linearly stable at the inertia-less (or purely elastic, i.e., Re → 0) limit.^{155,165} Attention has thus been turned to possible nonlinear mechanisms for instability.

A linear instability can be triggered with infinitesimal disturbances, whereas a nonlinear instability requires a finite-amplitude disturbance—one that exceeds a certain threshold magnitude. A classic example for the latter is the L-T transition in Newtonian fluids: for pipe flow, it typically occurs at Re_{crit} ≈ 2100 (defined based on bulk velocity) with finite-amplitude disturbances but can be delayed to much higher Re if disturbances are well controlled and minimized.^{21} It has been speculated that elastic instability in parallel shear flow could follow a similar scenario along the Wi axis.^{155} Evidence for such a nonlinear transition to purely elastic [Re ≲ *O*(10^{−2})] instability was shown by Pan *et al.*^{166} in microfluidic (≈100 *μ*m) open channel flow.

DR flow systems can no longer be considered inertia-less. At finite Re, Zhang *et al.*^{167} performed a nonmodal linear stability analysis, which studies the transient growth of specific modes of disturbances via linear mechanisms. They found that a streak mode disturbance becomes transiently amplified with increasing Wi, although it does still eventually decay over time. The mechanism, however, is probably unrelated with any of the later found IEIs, as the amplification is caused by stronger turbulence production while conversion to polymer elastic energy still acts as a suppressing force for the growth of TKE [see Eq. (55)]. Several types of IEIs in viscoelastic pipe and channel flows were reported in latest studies, and at least some of them are linear in nature.^{147,150,168} In particular, the so-called *elastoinertial turbulence* (*EIT*), discovered and named by Samanta *et al.*,^{147} has brought significant developments to the understanding of DR, which is the focus of this section.

##### d. Further notes on the terminology.

Fundamentally, EIT is a type of flow instability that is driven at least in part by fluid elasticity and likely of inertio-elastic nature. The latter statement is because so far EIT has only been reported for finite Re, but without a clear mechanistic understanding, the possibility of it being purely elastic cannot as yet be ruled out. Samanta *et al.*^{147} invoked the word “turbulence” to describe this instability likely to reflect its chaotic flow patterns, much like “turbulence” as in “elastic turbulence” in the work of Groisman and Steinberg.^{159} In addition, in the context of DR, EIT is viewed as a new state of turbulence that emerges at high levels of polymer elasticity. Some may frown upon such a loose reference to the term “turbulence,” especially if they subscribe to a narrower concept of turbulence which must satisfy a certain set of criteria, such as the seven characteristics put forth by Tennekes and Lumley.^{169} Without getting into the philosophical debate about what counts as turbulence, the term EIT is used here simply to follow its wide adoption in the recent literature.

On the other hand, one may also encounter the term “elastoinertial turbulence” being used for flow states of entirely different natures. For example, Gillissen^{170} used it to describe a 2D viscoelastic decaying homogeneous isotropic turbulence, even though the fluctuations are fueled solely by inertia and elasticity plays a suppressing role. To clarify, in this review, as well as in much of the DR literature, EIT refers to a specific category of instabilities that are (1) driven by both inertia and elasticity, (2) showing turbulencelike chaotic flow patterns, and (3) self-sustaining in time. In DNS, the first criterion can be unequivocally determined by examining the TKE balance equation (55): elasticity becomes a driving-force for instabilities if and only if $\u03f5pk<0$ (i.e., $\u2212\u03f5pk>0$).

Likewise, EDIs showing turbulencelike chaotic flow patterns will be termed *elasticity-driven turbulence* (*EDT*), whereas conventional turbulence sustained by inertial instability will be referred to as *inertia-driven turbulence* (*IDT*). Many in the literature called the latter “Newtonian turbulence” because the turbulence generation mechanism stays the same, at least qualitatively, as that in Newtonian flow, only that in viscoelastic flows, turbulence is attenuated by DR polymers. That term is a bit ambiguous and even misleading since it may be misconstrued as turbulence of Newtonian fluids only. In Table I, EDT and EIT are marked at high-Wi regimes as a representative scenario. In reality, whether and when do these turbulencelike instabilities occur depends on the specific flow type.

#### 2. Phenomenology

##### a. Transitions between different types of turbulence.

Samanta *et al.*^{147} measured the Re_{crit} for the L-T transition in Newtonian (water) and viscoelastic (aqueous solution of PAM with the same molecular weight and varying concentration *C*_{p} up to 500 wppm) fluids in pipe flow. At the same pipe diameter, they found that Re_{crit} initially increases with *C*_{p}—i.e., drag-reducing polymers delay the transition, which is intuitive considering their role in resisting turbulent motion. At higher *C*_{p}, however, Re_{crit} starts to decrease and eventually drops well below the value of Newtonian flow—the so-called “early turbulence” phenomenon. Both delayed^{171–173} and early^{174,175} transitions have been previously reported in various studies, which seemed conflicting at first. Observing both behaviors in the same system (same polymer-solvent pair and same flow setup) with varying *C*_{p} suggests that these are two coexisting transition pathways determined by the level of fluid elasticity.

Choueiri, Lopez, and Hof ^{64} later conducted a more comprehensive investigation of flow behaviors in a Re-*C*_{p} parameter space, as shown in Fig. 13. The boundary between laminar (white) and turbulent (shaded) regions, i.e., Re_{crit}, clearly varies nonmonotonically with *C*_{p}, with a peak at *C*_{p,crit} ≈ 25 wppm. For *C*_{p} < *C*_{p,crit}, increasing Re at constant *C*_{p} would see the occurrence of localized turbulent patches, the so-called “puffs” and “slugs,”^{176} before IDT (the authors labeled it as “Newtonian turbulence” in the figure) fills the entire domain. This transition sequence is identical to that in Newtonian flow except that Re_{crit} becomes higher with increasing *C*_{p}. For *C*_{p} > *C*_{p,crit}, the transition directly enters, without discernible structural localization, a new type of turbulence, which is presumed to be EIT.

If we instead explore the parameter space horizontally (i.e., along constant Re lines), for Re between ≈2300 (Newtonian Re_{crit}) and ≈3600 (highest point of the L-T borderline, found at *C*_{p,crit}), as *C*_{p} increases, the flow would laminarize first before reentering the turbulence zone. Choueiri, Lopez, and Hof ^{64} also performed a test where *C*_{p} was ramped up at a slow pace with Re controlled at 3150. The Newtonian limit shows well-developed IDT with a friction factor following the Blasius correlation.^{24} Within IDT, the friction factor drops sharply with increasing *C*_{p}. After *C*_{p} reaches the approximate range of 15 wppm to 20 wppm, the flow becomes intermittent and consists of extended quiescent regions separated by occasional bursts of turbulent activities. For a window around *C*_{p} ≈ 20 wppm, the friction factor matches that of the Virk MDR asymptote. This stage is strongly reminiscent of the high intermittency between hibernating turbulence and active/bursting phases reported by Xi and Graham^{121,122} and Zhu *et al.*^{54} (see Secs. III A 2 and III A 3). Like those earlier simulations, the flow relaminarizes at *C*_{p} ≳ 20 wppm. Turbulence, however, returns at *C*_{p} ≳ 45 wppm. The friction factor later converges again to the Virk MDR asymptote at *C*_{p} ≈ 60 wppm and stays nearly constant with further increasing *C*_{p}.

Although the two regimes of *C*_{p} ≈ 20 wppm and *C*_{p} ≳ 60 wppm share nearly the same friction factor, and thus DR%, and both agree with the Virk asymptote, the fact that they are separated by a laminar window clearly indicates two distinctly different stages of turbulence. Spatiotemporal patterns of streamwise velocity also appear different between them. Although both show smooth elongated streaks, the former (IDT with intermittency) also shows sporadic bursts of strong turbulence whilst the latter (presumed EIT) does not. Since the latter stage is the asymptotic limit of high elasticity, the authors proposed, same as Samanta *et al.*,^{147} that EIT is the ultimate MDR state.

Lopez, Choueiri, and Hof ^{46} performed DNS of pipe flow to complement the above experiments.^{64} The study focused on one Re and explored different regimes by adjusting Wi. Several simulation domain lengths were tested: in all cases, IDT is quenched and the flow laminarizes at sufficiently high Wi. Consistent with earlier studies,^{43,57,146} the observed laminarization is always preceded by a regime of high intermittency in IDT. The intermittency is temporal in shorter domains and resembles the AHB cycles in MFUs,^{54,121,122} but given sufficient domain length, it becomes spatial and takes the form of localized turbulence separated by laminarlike regions. The localized structures appear first as slugs and, at higher Wi, change to puffs before laminarization. This is strongly reminiscent of the L-T transition where the same group of structures appear but in an opposite order: with increasing Re, it is laminar →puffs →slugs →space-filling IDT. Events leading up to laminarization with increasing Wi were thus dubbed “reverse transition” by the authors.^{46,64} In another agreement with the work of Xi and Graham,^{43} DR converges to an asymptotic level in a small window of Wi immediately before laminarization. This asymptotic DR level gets higher as the domain length increases, and in the longest domain reported (50*D*), where slugs and puffs occur, it agrees with the Virk MDR. At a very high finite extensibility parameter *b* = 40 000, the study reported a second stage of turbulence with clearly different flow structures, which was believed to be EIT. The same as experiments at comparable Re, it occurs after a window of laminarization. (AD was used in the study which may affect the existence boundary of EIT—see Sec. II A 5 f.)

Back to the parameter space (Fig. 13), with increasing Re, both the window of laminarization and that of turbulence localization (puffs and slugs) shrink and eventually vanish, after which the transition from IDT to EIT and the ultimate convergence to MDR become a continuous process, as typically observed in the literature. The boundary of the MDR regime, which the authors presumed to be complete EIT, can be determined based on the friction factor magnitude. However, for the regime immediately preceding MDR, which is labeled “EIT dominant” in Fig. 13, its underlying turbulent dynamics is harder to determine using only experimental means. Choueiri, Lopez, and Hof ^{64} also referred to it as a “coexistence phase,” which is probable. Recent DNS by Zhu *et al.*^{54} reported the observation of EIT-like structures in a thin near-wall layer while the rest of the channel is populated by typical IDT structures.

##### b. Characteristics of flow and polymer conformation fields.

The most compelling evidence for EIT being a different state of turbulence is its distinctive flow structures observed in DNS, which was first reported in the work of Samanta *et al.*^{147} and confirmed by a number of studies.^{45,46,177–179} A typical image is shown in Fig. 14. In contrast to IDT (see Figs. 6, 7, and 10), where vortices in the buffer layer align in the streamwise direction (their downstream head may lift up), vortices at EIT align spanwise. They appear in two characteristic sizes with larger rolls separated by thinner threads. The polymer conformation field features distinct thin tilted sheets of high polymer extension [large tr(*α*)]—bright stripes in the color contours shown on the side panel of Fig. 14. This also differs from IDT where structures of polymer conformation align closely with vortices.^{51,66,106}

Recent DNS of Sid, Terrapon, and Dubief ^{45} found flow states with similar spanwise vortices and tilted tr(*α*) sheets in an *xy*-2D channel (i.e., no spanwise dependence), which strongly indicates that the underlying instability for EIT is intrinsically 2D. EIT in 3D DNS appears more like the 2D instability as Wi increases, presumably because 3D structures of IDT, especially those of active turbulence, become increasingly suppressed at higher Wi. Importance of 2D instability in EIT was further confirmed by Shekar *et al.*^{179} who found its velocity spectrum to be dominated by spanwise-independent modes.

#### 3. Origin of instability

##### a. Elasticity for instability.

It is clear that EIT is a new type of turbulence following a different self-sustaining mechanism. Evidence for an elastic mechanism is abundant. Both experiments and DNS have shown that EIT can exist at Re well below the Re_{crit} for the Newtonian L-T transition^{64,147,178}—to see this, compare the L-T borderline in Fig. 13 between the Newtonian and high *C*_{p} limits. Early occurrence of turbulence indicates that an additional source for flow instability must at least supplement, and possibly override, the classical inertia-driven mechanism, which in viscoelastic fluids can only come from elasticity. In addition, Samanta *et al.*^{147} tested the laminar-EIT transition with different pipe diameters [for the same polymer solution and same pipe diameter, Re and Wi are directly correlated—see Eq. (50)] and found that the transition occurs at the same critical Wi for different Re.

Direct measure of the instability driving mechanism is given in the TKE balance [Eq. (55)], which can only be evaluated in DNS. Turbulence generation from inertial and elastic mechanisms is measured by the production $Pk$ and elastic conversion $\u03f5pk$ terms, respectively. In IDT, polymers suppress TKE and thus $\u2212\u03f5pk<0$. For elasticity to become a driving force for turbulence $\u2212\u03f5pk$ must turn positive. In this sense, EIT was found in earlier DNS studies by Min *et al.*,^{38} Min, Choi, and Yoo,^{72} and Dallas, Vassilicos, and Hewitt.^{52} As discussed above (Sec. II C 3 b), positive $\u2212\u03f5pk$ appears in a thin slab around *y*^{+} ≈ 10–15 at LDR, which expands across most of the domain at MDR and at least some HDR cases. Positive $\u2212\u03f5pk$ was also observed in the transient development of an oblique mode disturbance in the nonmodal linear stability analysis of Zhang *et al.*^{167}

The experimental discovery of EIT prompted a more systematic study by Dubief, Terrapon, and Soria,^{178} where $\u2212\u03f5pk$ at EIT (as identified by its characteristic spanwise vortices) was found to be positive for all *y*^{+} ≳ 8. $P\u2009k$ is still positive, but its magnitude depends on Wi. For Re_{τ} ≈ 130, $P\u2009k$ is comparable to $\u2212\u03f5pk$ at Wi = 96, but at Wi = 720, it becomes insignificant, suggesting a diminishing role of inertia in EIT with increasing Wi. Similarly, at Re_{τ} = 84.85, Sid, Terrapon, and Dubief^{45} noted that even within in the broadly defined EIT regime, flow structures in their 3D simulation evolve with Wi. At Wi = 40, streamwise vortex structures typical of IDT appear intermittently and can often been seen blended with the spanwise rolls of EIT. They then claimed that those structures vanish in their Wi = 100 case where the 3D EIT solution more closely resembles the 2D one.

This brings up the possibility that the so-called EIT may be a combination of an underlying 2D instability, which could also be driven purely by elasticity (i.e., EDT), and an ensemble of intermittently occurring IDT states. The latter become less important at high Wi. Note that Samanta *et al.*^{147} proposed the term EIT prior to any knowledge of its instability mechanism. It was chosen mainly because inertial effects were presumed to be relevant at the high Re where they are found. A purely elastic mechanism, which of course awaits further validation, would not contradict that original study.

##### b. Linear vs nonlinear instability.

The laminar-IDT transition is nonlinear and Re_{crit,IDT} would increase when disturbances are reduced. In the Newtonian pipe flow of Samanta *et al.*,^{147} removing the imposed external perturbation delays the transition from Re_{crit,IDT} ≈ 2000 to ≈6500. Meanwhile, for the laminar-EIT transition, which is dominant at *C*_{p} > *C*_{p,crit}, they showed that the unperturbed and perturbed transitions occur at the same Re_{crit,EIT}. This clearly indicates a linear instability mechanism whose Re_{crit,EIT} is unaffected by the disturbance magnitude. At *C*_{p} < *C*_{p,crit}, both transition pathways can coexist in a certain range of *C*_{p}. Reducing disturbances would cause the laminar-IDT border to retreat to higher Re_{crit,IDT} and expose more of the laminar-EIT border (see Fig. 13). Indeed, Samanta *et al.*^{147} did find that at one lower *C*_{p}, where the perturbed transition is delayed $Recrit, pperturb>Recrit,Newt.perturb$ as expected from the laminar-IDT pathway, removing the external perturbation exposes a new transition point $Recrit, punpert.$ and $Recrit, pperturb<Recrit, punpert.\u226aRecrit,Newt.unpert.$. This $Recrit, punpert.$ is presumably the laminar-EIT border Re_{crit,EIT} normally hidden behind the laminar-IDT border Re_{crit,IDT} when sufficient disturbances are present.

In DNS, however, finite-amplitude disturbance is required to trigger EIT from the laminar state.^{46,178,179} Recent evidence also indicates a dependence on the specific form of disturbance.^{177} Both observations contradict the experimental conclusion that EIT is a linear instability.

Two latest stability analysis studies also pointed toward opposite directions. The first, by Garg *et al.*,^{168} found a linear instability in viscoelastic pipe flow at finitely large Re [≳*O*(10^{2}) and depends on polymer properties]. It was a 2D linear stability analysis in the longitudinal plane (spanned by axial and radial directions), the same plane on which 2D EIT solutions are found (spanned by streamwise and wall-normal directions). Unlike the more commonly studied Re → 0 limit, which as discussed above (Sec. III B 1) is linearly stable for viscoelastic parallel shear flow, the study found that at higher Re, a single unstable mode emerges in the laminar state. The instability consists of a streamwise array of spanwise rolls and neighboring vortices rotate in opposite directions, which resembles the train of spanwise vortices found in EIT. The structure is, however, localized in a layer adjacent to the pipe axis, whereas EIT structures are known to originate from the walls. A similar instability was also reported in channel flow.^{180}

The second, by Shekar *et al.*,^{179} performed linear stability analysis in viscoelastic channel flow at the dominant wavenumbers of the EIT solution. The laminar state was found to be linearly stable. However, the least stable eigenmode appears strikingly similar to the corresponding Fourier mode (same wavenumbers) of EIT. Slow decay of this mode makes it susceptible to nonlinear amplification of finite-amplitude disturbances. The most-amplified mode again closely resembles the same Fourier mode of EIT. The study further connected EIT to the Tollmien-Schlichting (T-S) wave, a self-sustaining TW solution to the N-S equation in channel flow, which, for Newtonian flow, stems from a linear instability of the laminar state at Re_{τ} ≈ 107.4.^{136,137,181} The bifurcation is subcritical and the solution extends to much lower Re (down to Re_{τ} ≈ 75 at the wavenumbers studied by Shekar *et al.*^{179}), where finite-amplitude disturbance would be required for it to be triggered (i.e., nonlinear instability). For comparison, however, EIT has been found at Re_{τ} as low as 40.^{178}

The study computed fully nonlinear T-S wave solutions for viscoelastic channel flow, where sheets of high polymer extension, a signature structure of EIT (Fig. 14), were observed. They are generated by near-wall stagnation points in the T-S wave, in a manner similar to the formation of a birefringent strand in cross-slot flow.^{162} This offers an explanation for an otherwise peculiar structural feature of EIT. Unlike the linear instability of Garg *et al.*,^{168} polymer sheets in T-S wave are localized near the walls. However, Shekar *et al.*^{179} were not able to find direct overlap between T-S wave and EIT in the parameter space—the latter exists at much higher Wi than the former. Therefore, any specific dynamical connection between the two, if existent, likely involves additional nonlinear mechanisms. Generalization to other flow types is also a challenge. For example, T-S wave is not found in pipe flow which is believed to be linearly stable for all Re.^{136,176}

#### 4. Discussion: Implications for MDR

Discovery of EIT not only provides direct proof for the existence of instabilities in parallel shear flow that are (fully or partially) elastic in nature, but it also reconciles the seemingly conflicting observations of both delayed turbulence and early turbulence in drag-reducing polymer solutions. Furthermore, it brings new perspectives and raises new questions for the MDR problem, especially in light of the major knowledge gaps in the dynamical systems framework (Sec. III A 4). Considering this latest development, three possible scenarios are discussed here for the asymptotic limit of high polymer elasticity.

##### a. Scenario I: EIT is the single form of instability underlying MDR.

This is the most straightforward possibility and also seems to be the opinion held by most researchers. Experimental and numerical findings reviewed above seem to demonstrate that, at least at the relatively moderate Re examined, EIT is the ultimate state of the high polymer elasticity limit. This is the conclusion from the work of Choueiri, Lopez, and Hof,^{64} as shown in the parameter space of Fig. 13 and also supported by the diminishing presence of IDT structures with increasing Wi.^{45} The dynamical AHB cycle, in its originally proposed form described in Sec. III A 3, was clearly also observed in those studies but only for a limited window of *C*_{p} or Wi, which eventually gives way to EIT—a self-sustaining process itself.^{46,64} This addresses the question of the “existence” of MDR (see three attributes of MDR in Sec. II B 3 c). In addition, mean velocity of EIT seems to follow the Virk MDR asymptote,^{46,64,147,178} in accordance with the “magnitude” attribute of MDR.

However, it would appear inexplicable for a flow instability, which relies on elasticity, to be independent of polymer solution properties—the attribute of “universality” is a conspicuous gap to fill. The scenario will also have difficulty with the experimental observation that solutions of rigid polymers, which, under FENE-P, would be described by a very low *b* parameter [Eq. (34)] and presumably cannot support EIT, are bounded by the same Virk MDR asymptote.^{182,183}

##### b. Scenario II: EIT is a phenomenological reflection of a new form of the AHB cycle.

This scenario interprets the term EIT in the existing literature not as a reference to a single form of flow instability but as an imprecise umbrella term encompassing many flow states with varying extent of elastic effects. Indeed, DNS results do show the presence of different amounts of IDT structures blended with spanwise rolls in 3D flow states generically described as EIT.^{45,177,178} It is highly possible that the complex appearance of 3D EIT simply reflects a statistical ensemble containing both states dominated by IDT and those belonging to a “pure” form of EIT—the latter is likely the corresponding 2D EIT.^{45,168,179} The 2D form could also be purely elastic (i.e., EDT) rather than elastoinertial, considering that the magnitude of average TKE production $P\u2009k$ in 3D EIT diminishes with increasing Wi.^{178} In this case, the observed inertial characteristics would be solely attributed to the intermittent occurrences of IDT.

In this scenario, intermittent AHB cycles are likely still involved, except that their specific dynamics has to change at higher Wi. The original theory, discussed in Sec. III A 3 and depicted in Fig. 11, would still be valid but only up to a certain Wi. At sufficiently high Wi, the EoC, which is originally in charge of shielding the flow from laminarization, becomes penetrable. EIT in its pure form (or EDT if it is purely elastic) takes over to keep turbulence self-sustaining (i.e., prevent laminarization). In 3D flow, IDT still grows intermittently to initiate the active phases of turbulence, which would be quickly suppressed by polymers. The destruction of EoC and emergence of EIT/EDT do not always coincide, which would explain the appearance of a laminar window at lower Re (Fig. 13).

The scenario’s reliance on intermittency does not necessarily contradict the seeming convergence to an “EIT” state observed by Choueiri, Lopez, and Hof.^{64} (Indeed, flow structures observed in both experiments and DNS do not strictly converge with *C*_{p} or Wi.^{46,64}) The observed “converged EIT” state could be an evolving dynamical cycle between IDT and EIT/EDT, only that variations between states are difficult to detect from streamwise velocity patterns, which are relied on in experiments for flow state determination. The lack of apparent bursting in velocity patterns also does not preclude the existence of active IDT, since under polymer stress the flow can bypass the strongest form of bursting and energy eruption *en route* to active turbulence.^{54}

For the attributes of MDR, this scenario again offers explanation for “existence” but with a self-sustaining mechanism that is more dynamical in nature. The challenge with “universality” still exists, similar to the previous scenario, but there appear to be more avenues for its explanation. For example, because the intermittent occurrence of active turbulence keeps the flow from fully reaching the pure form of EIT/EDT, its mean flow is not the same as that of the latter but determined by the dynamical balance between EIT/EDT and IDT. Both EIT/EDT and IDT statistics could change with increasing Wi. How do their dynamical average stays constant at the Virk level would be the key to addressing both the “universality” and “magnitude” of MDR.

##### c. Scenario III: EIT is only important at low Re and/or in small domains.

This scenario assumes that IDT would not be terminated by polymer effects and replaced by EIT, if higher Re or larger flow domains are used. It will only be suppressed to a weaker form in which polymer effects are minimal—i.e., the theory for MDR based on AHB cycles still holds. Compared with previous scenarios, here, EIT is no longer a part of the self-sustaining process of MDR, which thus easily addresses its “universality.”

Existing knowledge of EIT and its related transition scenario was based on studies at fairly low Re (in turbulence standards) where the total thickness of the wall-normal layer, from the wall to the channel center, equals Re_{τ} ≲ *O*(100) wall units. This limits IDT to primarily buffer-layer structures, such as streamwise vortices, and only some lower log-law layer structures, while leaving out structures of the upper log-law and outer layers. Complexity of turbulent structures and dynamics increases with *y*^{+}.^{28,83} It is thus possible that structures at higher *y*^{+} are not affected by polymers in the same way as low-Re IDT and their self-sustaining cycles are not completely disrupted by polymer elasticity.

From a scaling argument, since the log-law layer extends from ≈30 wall units to ∼*O*(0.1)*l* (Sec. II A 3), it is not fully developed and the inner and outer layers are not completely separate unless *l*/*l*_{v} = Re_{τ} ≫ *O*(100). Indeed, in Newtonian flow, the mean velocity profile does not even fully collapse on to the von Kármán law [Eq. (18)] until Re_{τ} ≳ 400.^{58,79} At Re_{τ} > 5000, a secondary peak would arise at higher *y*^{+} in the $vx,rms\u2032+$ profile, reflecting the emergence of new structures not captured at lower Re_{τ}.^{184} Large-scale coherent structures at high Re_{τ} have attracted much interest among turbulence researchers in recent years. One example is the so-called very-large-scale motions (VLSMs), multiple packets of strong turbulent structures organized into higher-order patterns, which are typically studied at the Re_{τ} ∼ *O*(10^{3}) regime.^{185} Another example is the tall attached structures studied by Lozano-Durán, Flores, and Jiménez,^{59} which originate from the buffer layer but can extend across the entire channel height at Re_{τ} ≈ 2000.

Although the TKE production rate peaks in the buffer layer for all Re, at higher Re_{τ}, the log-law and outer layers are larger in volume. It was thus estimated that at Re_{τ} ≈ 4200 their accumulated contribution to TKE production exceeds that of the buffer layer.^{85} Whether EIT can still overtake IDT in high-Re regimes where vastly different coherent structures are at play is anything but certain.

Domain size constraints, especially in the streamwise direction, could also be a factor. We know that, in the streamwise direction, the velocity correlation length and minimal domain size required for self-sustaining IDT increase rapidly with Wi and DR%.^{46,73,146} It is well possible that in a highly extended domain beyond the current computational power, IDT, in a form with strong spatial intermittency (i.e., puffs separated by large laminarlike regions^{46}), can persist to higher Wi. The flow may eventually laminarize if the “reverse transition” hypothesis^{46,64} is proven accurate. However, as long as this ultimate Wi for laminarization is too high to be practically relevant, spatially localized turbulence would appear as the asymptotic state. Obviously, domain size is only meaningful in simulation, which is however the only way the nature of instability (IDT, EIT, or EDT) can be unambiguously determined. Experimentally, it is inferred from flow patterns, which is how Fig. 13 was obtained.^{64} For this reason, the above discussion of the domain size factor is of general relevance.

### C. Dynamical analysis of coherent structures in extended domains

Recent progress reviewed in Secs. III A and III B builds on a shift of focus from ensemble average quantities to spatiotemporal intermittency. Numerical simulation played a pivotal role in those developments. MFU-based DNS is the most common approach for tracking the temporal evolution of individual structures, which is complemented by the direct computation of invariant solutions (ECS). Although many DNS studies reviewed in Secs. III A and III B were not performed in domains that are strictly minimized, they are still sufficiently constrained that dynamical evolution of structures is reflected in the time series of domain-average quantities. For DNS in more extended domains, patches of different flow states can be analyzed by invoking the equivalency between temporal intermittency and spatial intermittency, as in the work of Wang, Shekar, and Graham.^{130} This naïve approach views the dynamical life cycles of turbulent structures as statistical sampling of individual flow states in an ensemble, without any correlation or interaction between structures. Such a correlation becomes increasingly important at higher Wi, which is evident from the localization of vortices and their organization into slugs and puffs.^{46,57} This dynamics is excluded in MFU but fully captured in DNS in extended domains. Effective extraction of such information from complex turbulent flow fields is a nontrivial challenge.

In the area of DR, two major approaches were most used in the literature to analyze coherent structures in large-scale turbulent flow fields. One is KL analysis or proper-orthogonal decomposition (POD),^{24,95,108,186} which decomposes velocity fields into orthogonal spatial modes—a basis set for velocity fields. None of the modes directly represent any real occurring coherent structure, but the leading modes can be viewed as the most important constituents (in terms of the amount of kinetic energy contained) of statistically representative structures. Wang *et al.*^{146} extended the analysis to also include the polymer conformation field in the decomposition. Another approach is conditional sampling, which finds the statistical representation of local flow fields around a predefined event of interest. The event can be the occurrence of streamwise vortices, in which case patches of flow fields around qualifying vortices are aligned at their center axes and statistically averaged.^{109,187} It can also be high streak intensity, such as in the work of Kim *et al.*^{106} in which patches are selected based on the detection of strong eruption ($v$_{x}′ < 0, $v$_{y}′ > 0) events making largest contributions to the RSS [Eq. (54)]. Compared with KL analysis, conditional sampling is a local approach with no direct tie to the imposed simulation domain. Its results thus connect more closely with individual coherent structures. However, its outcome is inevitably influenced by the subjectivity in choosing the selection condition for the event.

Those earlier approaches generate statistical representations of coherent structures and neglect their individuality. Recent attention to spatiotemporal intermittency calls for the capability to extract the real instances of individual structures and trace their temporal evolution in complex flow fields. The first step, i.e., extraction of structural instances, is static and recently made possible by a method called VATIP (vortex axis tracking by iterative propagation).^{58} This section takes the mechanism of LDR-HDR transition as a case study for illustrating the application of VATIP in vortex conformation analysis. Future work of extending VATIP for dynamical analysis (the second step listed above) is also briefly discussed.

#### 1. Case study: Vortex regeneration mechanism for HDR

Distinctions between LDR and HDR were thoroughly studied by Zhu *et al.*^{57} and discussed in Sec. II B 2 d. The LDR-HDR transition is marked by sharp changes in a variety of flow statistics and believed to have its roots in the localization of vortex structures. Understanding the origin of localization requires the analysis of individual vortex objects and, when possible, their dynamics and interactions.

Zhu *et al.*^{57} proposed a mechanistic explanation citing the conceptual framework of vortex regeneration mechanisms summarized by Schoppa and Hussain.^{188} In this framework, new turbulent vortices are continually generated from existing ones following two distinct pathways. The first, the streak-instability pathway, is an extension of the classical conceptual model for the turbulent self-sustaining process by Waleffe.^{189} Counter-rotating streamwise rolls (straight vortices) generate low speed streaks in between. Adding streamwise sinuous perturbations can trigger the so-called “streak breakdown”^{190}—strong instabilities leading to streak intensification and vortex growth. Those vortices have a strong tendency to lift up—a process in which the downstream end of a vortex pivots away from the wall and extends vertically into higher *y*^{+} regions. Lifted vortex heads are often swung sideways by local transverse flows and many form the so-called “hairpins”—Ω-shaped vortices whose downstream head, the arch in Ω, is lifted up while both streamwise legs extend upstream and stay closer to the wall.^{81,191} Rapid lifting and growth of vortices ends with their sudden eruption—a bursting phase (Sec. III A 3 b), which was studied in detail by Zhu *et al.*^{54} The resulting fragments and intense fluctuations can spread and trigger new streak instability elsewhere to start the next cycles. An example displaying part of the process is shown in Fig. 15(a). A vortex tube marked with “A” merges with a nearby vortex to form a crescent at *T* = 8, which lifts up and grows into a full hairpin (*T* = 36) and then bursts (*T* = 40) into fragments (not shown), leaving only one leg in its residue (*T* = 48).

Zhu *et al.*^{57} proposed that strong polymer effects can block this cycle by suppressing vortex lift up. In a separate study, Zhu *et al.*^{54} also found that polymers can subdue bursting and suppress the production of high-intensity fluctuations. Without lifting and bursting, vortices are stabilized in the streamwise direction, which allows them to be extensively elongated, as clearly shown in the conditional eddies at HDR obtained by Kim *et al.*^{106} It would also disrupt the streak-instability pathway and expose the second so-called “parent-offspring” pathway to dominate regeneration cycles. In the latter, as illustrated in Fig. 15(b), vortices are successively generated in close proximity to form a chain—new “offspring” vortices are generated at the tip of an existing “parent,” in the shear layer between the parent and the wall. Obviously, this mechanism is intrinsically more local than the other. Both pathways exist in Newtonian flow, but at high Wi, the parent-offspring one becomes more prominent, leading to the apparent clustering and localization of vortices. Applying the VATIP algorithm (discussed below), Zhu and Xi^{98} showed that the suppression of vortex lift-up coincides with the LDR-HDR transition, which presents a compelling depiction of the second DR mechanism required for its explanation (see Sec. II B 2 d).

This explanation is compatible with the earlier theory based on AHB cycles (Sec. III A 3) as vortex lift up is part of the bursting and active phases and its suppression leads to higher presence of hibernation. It, however, goes one step further and explains the vortex localization mechanism, which requires interactions between vortices not considered in the AHB framework.

#### 2. Vortex conformation analysis by VATIP

The above mechanism was proposed based on direct anecdotal inspection of visualized flow field images (Fig. 15), which is too often relied on in turbulence research. Reliable conclusions cannot be reached without the objective and quantitative analysis of vortex conformation and dynamics. Although methods for vortex identification, i.e., derivation of quantitative indicators for vortex regions—such as *Q*, have been extensively studied (Sec. II A 6), they only tell if a region belongs to *any* vortex without identifying which regions constitute the same well-defined vortex, which is thus not adequate for vortex conformation analysis.

Each vortex has an axis of rotation which fully describes its topological shape and instantaneous conformation. An axis-line is formed by connecting the center of rotation—the axis-point—on each cross section of the vortex tube, which can be defined as the 2D maximum of *Q* (or any other indicator of vortex strength) on the plane. For linear streamwise vortices, grouping axis-points into axis-lines is relatively straightforward through a so-called “cone-detective” procedure by Jeong *et al.*^{192} Axis-lines so extracted provide the reference points for aligning streamwise vortices, which attracted much attention in earlier turbulence research, for conditional sampling.^{109,187,192}

As discussed in Sec. II C 1, existing studies on polymer-vortex interaction were mostly limited to streamwise vortices with few exceptions.^{53,106,193} Understanding how complex vortices found at higher *y*^{+}, such hairpins, are affected by polymers is essential for understanding DR at higher Re. As shown in Fig. 16, even at a fairly low Re_{τ} = 172.31, hairpin vortices are prevalent especially in Newtonian and LDR cases. Kim *et al.*^{193} extracted representative hairpin vortices using conditional sampling and studied their transient development with DNS. Their approach is akin to *in vitro* experiments in the sense that model hairpins are studied in isolation from their “living” environment. The new VATIP algorithm, developed by Zhu and Xi,^{58} is capable of extracting axis-lines of general 3D vortices, including hairpins, from near-wall turbulence “*in vivo*.”

VATIP borrows the cone-detective idea from Jeong *et al.*^{192} but extends the search to all three dimensions to capture vortices with nonlinear and branched shapes. Without getting into the details, the algorithm appears like the video game “snake,” where growing axis-lines keep looking for new axis-points to “swallow” (i.e., connect to). If no more eligible axis-points are found in their growing direction, the propagating ends turn to the next search direction. Iteration between different search directions continues until none of the axis-lines can grow further. The method reliably captures all known types of vortices generated by a no-slip wall, including both symmetric and asymmetric hairpins. This is clearly shown in Fig. 16 where axis-lines extracted by VATIP (colored dots) faithfully reproduce the contour shapes of all vortices (gray tubes) from Newtonian flow to HDR. Extracted axis-lines can be classified based on their topologies into categories such as streamwise vortices, hairpins, and branches. (It, however, may not work as well with structures having no direct interactions with the wall, which typically occur at much higher *y*^{+} than those studied so far for DR^{82,149}—see detailed discussion in the work of Zhu and Xi.^{58})

Application of VATIP to DNS data at different regimes of DR allows a comprehensive investigation of polymer effects on vortex conformation (Fig. 16).^{98} The results show that although streamwise vortices can be found either lying flat or lifted up, hairpins and hairpinlike vortices are almost always lifted, which supports the above postulation that hairpins are formed from the lift up of streamwise vortices. At LDR, vortex intensity is reduced, but their distribution pattern and conformation statistics remain largely unchanged from Newtonian flow. At HDR, the number of lifted vortices declines sharply. Suppression of vortex lift up interrupts the generation of hairpins and hairpinlike vortices. It also reduces the turbulent momentum flux between buffer (streamwise vortices) and log-law (lifted and hairpin vortices) layers, which would explain the changing flow statistics in the log-law layer between LDR and HDR. Overall, the results are consistent with the proposed mechanism in Sec. III C 1.

The current analysis is static. It shows, within any frozen instant of turbulence, what is the state and conformation of each vortex, without establishing the temporal connection between vortices found at different instants. Expansion of the VATIP framework to cover temporal dynamics is foreseeable. For instance, a recent method for temporal coherent structure analysis by Lozano-Durán and Jiménez^{194} may be adapted to the axis-lines extracted by VATIP. This would open the door for the detailed, quantitative, and *in situ* analysis of vortex life-time dynamics in realistic extended flow domains.

## IV. SUMMARY AND OUTLOOK

For over 70 years, turbulent drag reduction by polymers has remained an active area of research, as new approaches and new discoveries continue to emerge. The past decade (2010s) has witnessed a surge of interesting developments in the fundamental understanding of DR. Many of them have challenged the established line of thought and brought significant progress in answering some of the most puzzling questions—in particular, the nature of MDR. A central theme among those latest advances is to go beyond the ensemble statistics of fluctuating turbulent flow fields and study the dynamical evolution and intermittency of coherent flow structures.

The current review covers both classical contributions and recent developments in the fundamental inquiries into polymer DR. The ultimate target is the mechanistic understanding of transitions between different turbulence regimes, and the status of current research is summarized briefly as follows:

- Onset of DR:
This problem is equivalent to understanding the mechanism of polymer DR that leads to LDR. Great strides have been made in deciphering the interactions between polymers and turbulent flow structures, thanks to the advancement of numerical simulations over the past two decades. It is well established now that polymers can counteract and suppress streamwise vortices—dominant structures in turbulent buffer layer where the onset is believed to occur (Sec. II C 1). However, the debate of whether such effects should be wrapped theoretically as a viscous (Lumley) or elastic (de Gennes) mechanism of polymers remains unsettled.

- HDR:
Historically, HDR has received much less attention because of its later discovery but also because it is often presumed to be simply a precursor of MDR. There is certainly justification for this argument. Most notably, HDR is marked by the expansion of DR effects from a limited region (buffer layer) to nearly the entire channel.

^{57,99}Increasing length scale of drag-reduced turbulence to the scale of flow geometry is also central to several theories for MDR.^{5,15,92,115}On the other hand, if the nature of MDR is indeed at least partially elastic (Sec. III B), as believed by many, the fact that the LDR-HDR transition can be observed without EIT (or EDT)^{57}would indicate the decoupling between HDR and MDR. This transition reflects the start of a second DR mechanism, in relation to the first mechanism at the onset of DR. This new mechanism has been explained in terms of changing vortex regeneration mechanisms in recent studies.^{57,98}Despite some strong evidence from vortex conformation analysis (Sec. III C), further investigation, such as temporal analysis of vortex lifetime cycles, is still needed. - MDR:
This long-standing mystery finally seems to have opened up a few cracks. A complete theory must consistently address three attributes of MDR: existence (self-sustaining mechanism), universality (insensitivity to changing polymer solutions and domain size), and magnitude (empirical Virk profile). Two different but interconnected streams of developments occurred in the past decade. The first is the theoretical framework based on AHB cycles (Sec. III A 3). It offers a straightforward explanation for “universality.” As for “existence,” a plausible self-sustaining mechanism is proposed but so far cannot be fully reconciled with numerical observations (Sec. III A 4). The second builds on the discovery of new flow states driven by polymer elasticity—EIT (or possibility EDT). For the “existence” attribute, it offers a self-sustaining mechanism backed by experimental evidence. However, how to address the “universality” attribute remains an open question. Both streams fall short of a quantitative theory for the “magnitude” of MDR.

In addition to the outstanding questions summarized above, many research opportunities lie in areas that have been under-explored so far.

### A. High Re

Existing understanding and knowledge are mostly based on numerical simulations at fairly low Re. There have been some limited attempts at DNS at higher Re: the highest, to my knowledge, is Re_{τ} = 1000 in the work of Thais, Gatski, and Mompean^{76} and Pereira *et al.*^{195} However, in-depth analysis of coherent structures, flow states, and dynamical intermittency has been performed mostly at Re_{τ} ∼ *O*(100), where the state space is relatively simple. As discussed in Sec. III B 4 c, new families of coherent structures would emerge at higher Re, which brings new turbulent self-sustaining dynamics into the picture. Those complex structures are unlikely to follow the same interaction mechanism with polymers as streamwise vortices (described in Sec. II C 1). Both frameworks for MDR, based on AHB cycles and EIT, respectively, need to be reexamined for higher Re. The newly developed VATIP algorithm, for its capability of processing complex 3D vortex structures, could be instrumental (Sec. III C 2).

### B. Connecting numerical models with realistic polymer solutions

Numerical investigation has been overly reliant on “FENE-P polymers,” which (recall Sec. II A 5 b) is only an idealized model with substantial simplifications. There has not been much effort in quantitatively connecting simulation models with experimental drag-reducing solutions. Understanding which polymer-solvent combinations are more accurately modeled by FENE-P and how do model parameters map to realistic experimental systems is essential for making chemically specific predictions from DNS, which are of practical interest for the selection and design of drag reducers. It is also important for interpreting simulation results through the practical lens. For instance, Wi = 500 and *b* = 10^{4} map to what molecular weight and shear rate for a given polymer? What parameter ranges are practically relevant? Discussion in Sec. II A 5 e provides the theoretical basis for the determination of materials properties in FENE-P. In practice, parameters are more often obtained by fitting with experimental rheology data,^{47} and parameters from the best fit do not always match theoretical expectations.^{196}

### C. DR by rigid polymers

Comparison between flexible and rigid polymers is intriguing. Although both seem to be bounded by the same MDR asymptote, the transition patterns before reaching MDR are totally different. Based on this difference, Virk dubbed them type A (flexible) and type B (rigid) drag reducers.^{182,197} Streamwise velocity fluctuation patterns from flexible and rigid polymers also appear to differ^{33,75} (Sec. II B 2 b). Since FENE-P is a model for flexible (type A) polymers, insight into DR by rigid (type B) polymers is limited. Meanwhile, rigid polymers also provide a reference system where significant DR is achievable, but elastic instabilities are not expected, which can help us separate the roles of different mechanisms in DR.

### D. DR by surfactants

Surfactants are also highly potent drag reducers. They assemble into cylindrical “wormlike” micelles—long chain structures held together not by covalent bonds but through intermolecular (vdW and electrostatic) interactions. Wormlike micelles can be viewed as “living” polymers: they would break in regions with strong turbulence but can reassemble once external strain is removed. Experiments have shown that surfactant DR systems can exceed the Virk MDR asymptote,^{198,199} which indicates a different mechanism of DR, at least in that regime. This could be related to the “living” nature of the microstructure, or to the formation of higher-order assemblies such as shear-induced structures.^{200,201} Fundamental understanding is rather limited.

### E. Flow-induced scission of polymer chains

Long-chain polymers are subject to mechanical degradation in strong turbulence, which is responsible for the gradual loss of their drag-reducing capacity.^{202} This effect is not considered in current numerical models. Quantitative prediction of chain scission is, however, of practical significance. Any model for chain scission would be similar to those for surfactants as living polymers, less the reassembly part. One idea to combat flow-induced scission is to use supramolecular chemistry,^{203} where smaller polymeric building blocks are jointed together through noncovalent interactions to form much longer chains. The noncovalent bonding sites would break in strong flow, but, like surfactants, they can reconnect afterward.

## ACKNOWLEDGMENTS

Much gratitude is due to many co-workers (mentors, colleagues, students) with whom I have collaborated in this area or had enlightening discussions, in particular, Michael D. Graham (U Wisconsin-Madison), who initially led me into this area with patient guidance and visionary insights, and Lu Zhu (McMaster U.), whose recent contributions have injected new energy to this area of research in my group. Many thanks should also go to John F. Gibson (U. New Hampshire) for sharing his Newtonian DNS code ChannelFlow,^{204} which benefited our research as well as that of many others, to Hecke Schrobsdorff (Max Planck Inst.), Tobias Kreilos, and Tobias M. Schneider (EPFL) for helping with the parallelization of our in-house DNS code, and to Javier Jiménez (Polytechnic U. Madrid) for discussion that inspired our interest in vortex analysis, which led to VATIP. National Science Foundation (NSF) Grant No. NSF PHY11-25915 partially supported my stay at the Kavli Institute for Theoretical Physics (KITP) at UC Santa Barbara. Financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through its Discovery Grant Program (No. RGPIN-2014-04903), the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET: www.sharcnet.ca), and the computing resources allocation awarded by Compute/Calcul Canada are all acknowledged.

## REFERENCES

_{τ}=590