1. Introduction
Turbulent drag reduction is a common goal for many different applications, such as pipelines, military vehicles and commercial transport. For many such situations, it is necessary to understand how turbulence produces drag in order to mitigate it as much as possible. It is well known that wall-bounded turbulence has characteristic coherent structures (Fiedler Reference Fiedler1988) that contribute to Reynolds stresses, which in turn increase the skin-friction drag (Bernard, Thomas & Handler Reference Bernard, Thomas and Handler1993). One particularly effective way of disrupting these drag-producing structures is through the addition of polymers.
Since Toms (Reference Toms1948) demonstrated that dilute polymer solutions can significantly reduce drag, the field has seen substantial advancements through hundreds of experiments and numerical simulations. Lumley (Reference Lumley1969, Reference Lumley1973) synthesised many experimental results to characterise polymer drag reduction and proposed that the key factor for effective drag reduction was the relaxation time of polymer molecules relative to the local flow time scales, leading to the well-known Lumley criterion. His proposed mechanism for drag reduction was then based on a local effective viscosity of polymer chains which changed due to their elongation in the flow. Virk (Reference Virk1975) introduced the concept of a universal maximum drag-reduction (MDR) asymptote which establishes a theoretical ideal performance of any given polymer solution, although Choueiri, Lopez & Hof (Reference Choueiri, Lopez and Hof2018) recently demonstrated that, under certain circumstances, this limit can be exceeded with an effective relaminarization of the flow. Later, Tabor & de Gennes (Reference Tabor and de Gennes1986) proposed an alternative, energy-based mechanism due to the elastic effects of the long-chain polymers interfering with turbulent motions, motivated by the observation that drag reduction could occur in flows beyond near-wall regions. However, this elastic energy theory was developed for homogeneous isotropic turbulence, and its applicability to drag reduction in wall-bounded flows is unclear even with efforts by de Gennes (Reference de Gennes1986) to address this. More recently, experiments have shown a clear connection between maximum viscoelastic drag reduction and elasto-inertial turbulence (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013), in which small amounts of polymer introduce flow instabilities that induce transition at a lower Reynolds number than that for equivalent Newtonian flow (see Dubief, Terrapon & Hof Reference Dubief, Terrapon and Hof2023). In this state, the flow falls into the MDR regime and lends credence to the two-way energy transfer of de Genne’s theory (Dubief, Terrapon & Soria Reference Dubief, Terrapon and Soria2013). Recent simulations by Serafini et al. (Reference Serafini, Battista, Gualtieri and Casciola2025), however, emphasise the role of polymer stretching and the associated extensional viscosity in the onset of drag reduction, consistent with Lumley’s theory. Sreenivasan & White (Reference Sreenivasan and White2000) explored the connections between elastic energy-based models and predictions from effective viscosity frameworks, noting that, under some flow conditions, the two approaches can yield similar scaling. The apparent evidence supporting both theories and their limitations in universally describing polymer drag reduction have led to the development of nonlinear unified theories that integrate both viscous and elastic contributions (Procaccia, L’vov & Benzi Reference Procaccia, L’vov and Benzi2008). Nevertheless, the original theories of Lumley and de Gennes continue to guide many experiments and simulations leading to further characterisation of polymer drag reduction mechanisms (Benzi & Ching Reference Benzi and Ching2018). Further discussion on the evidence and limitations of the two theories can be found in the review of Xi (Reference Xi2019).
Regardless of the precise mechanism responsible for the observed drag reduction by polymer additives, the effects are best understood within the context of a polymer model. For this numerical study, we choose to use the finitely extensible, nonlinear elastic dumbbell model with the Peterlin approximation (FENE-P) model, which was developed by Wedgewood & Bird (Reference Wedgewood and Bird1988) and, despite its simplified description of the polymers, has proven to be the best choice of polymer model in most cases (White & Mungal Reference White and Mungal2008). Under the FENE-P model, a polymer molecule is modelled as a pair of masses connected by a spring, and the macroscopic dynamics of the flow are determined by the evolution of the end-to-end vector connecting the two masses, averaged across all polymer molecules, represented by the so-called conformation tensor,
$C_{\textit{ij}}$
, which describes the orientation and stretching of the molecules. As computational power and technology have increased, however, it has become more practical for some to use a coupled Lagrangian simulation for individual polymer molecules with other fluid solvers. Mortimer & Fairweather (Reference Mortimer and Fairweather2022) used this kind of approach in direct numerical simulations (DNSs) to study the relationship between polymer extension and turbulent boundary layer features. They confirmed that the extension of the polymer is greatest in the buffer layer since this is where the largest streamwise velocity gradients are found. Serafini et al. (Reference Serafini, Battista, Gualtieri and Casciola2022) demonstrated that this sophisticated hybrid method could more accurately simulate polymer drag reduction in pipe flows beyond the capabilities of standard DNSs with the FENE-P model. Serafini et al. (Reference Serafini, Battista, Gualtieri and Casciola2023, Reference Serafini, Battista, Gualtieri and Casciola2024) further used this method to examine how polymer solutions affect the kinetic energy budget in a wall-bounded turbulent flow, comparing Lagrangian modelling for polymer molecules with the more common FENE-P model, and they showed that the FENE-P model tends to overestimate drag reduction. However, the FENE-P model is still widely used and the current study offers a useful framework under this model with which to study non-uniform polymers and their effects on wall-bounded turbulence.
The ability to simulate non-uniform polymer distributions is not only important for more realistic modelling of inhomogeneous polymer solutions, but also allows us to address the operational and environmental challenges of polymer drag reduction. This naturally motivates an investigation into targeted application methods that maximise drag reduction per unit of polymer introduced. While polymer drag reduction has proven successful in closed systems such as pipe flow, it is not as practical for many external flows like those of marine vehicles (White & Mungal Reference White and Mungal2008). With such external flows, polymer cannot be easily recycled into the flow, and thus there are concerns with both the cost of supplying enough polymer additives to be continually effective and the ecological impact of leaving polymer additives in bodies of water. Both of these concerns can be addressed in part by works such as that of Rajappan & McKinley (Reference Rajappan and McKinley2019) and Xie et al. (Reference Xie, Jiang, Meng, Li, Wen and Hu2023), which are focused on producing cheap, biodegradable alternatives to typical polymers. However, it is still necessary to understand how to efficiently use polymer additives to target vital regions of drag-producing structures.
While simulations and models of inhomogeneous polymer solutions have been around for some time (Mavrantzas & Beris Reference Mavrantzas and Beris1992), the idea of a targeted approach to polymer drag reduction is relatively new. Larson (Reference Larson2003) modelled polymer injection into an external boundary layer and discussed physical considerations for polymer effectiveness and the problem of turbulent transport removing material from the near-wall region where it is more effective. The author’s conclusions, however, were that injectors would need to be spaced at least 4–10 cm apart in the streamwise direction to achieve adequate coverage without losing material to the outer region of the boundary layer. Such a stringent requirement makes standard injection methods infeasible. The most relevant work on the subject of practically targeting specific regions in the flow was first done by Yang et al. (Reference Yang, Park, Chun, Ceccio, Perlin and Lee2014) in which they mixed a powder of poly(ethylene oxide) (PEO) with an anti-fouling paint. When applied to a ship’s hull, the shear induced by the motion of the ship ablated the coating and the PEO powder mixed with the water to produce a dilute polymer solution which reduced skin-friction drag by up to 33 % more than a typical anti-fouling paint alone. Recent studies have explored this idea further experimentally (Rowin et al. Reference Rowin, Asha, Narain and Ghaemi2021) and numerically (Du et al. Reference Du, Mutsuda, Wasada and Nakashima2023). This type of approach avoids the issue of having to inject polymer discretely from the ship’s hull via thousands of injectors, but it does not address where the polymer should go once it is in the boundary layer to achieve optimum drag reduction. One key drawback, however, is that it requires consistent reapplication of paint since it is stripped away during the voyage; according to the estimated polymer mass flux required for sufficient drag reduction, the paint would also need to contain at least 1 kg of PEO per square meter for roughly 30 h of travel (Larson Reference Larson2003), which is neither economical nor environmentally friendly. We are therefore motivated to explore alternative or complementary polymer addition methods which are informed by existing literature.
It is generally understood that polymers are most effective in the buffer layer where they oppose the streamwise vortices dominating the flow in that region. The experiments of McComb & Rabie (Reference McComb and Rabie1982a , Reference McComb and Rabieb ) demonstrated that injecting polymer directly into the buffer layer of a turbulent pipe flow induces a quicker onset of drag reduction and may even overshoot the asymptotic drag-reduction value obtained with a uniformly distributed polymer solution or injection at the centreline. A later DNS study of Kim et al. (Reference Kim, Li, Sureshkumar, Balachandar and Adrian2007) verified that counter-rotating streamwise vortex pairs are the most significant structures in the buffer layer and polymer interactions with these structures indeed contribute greatly to the thickening of the buffer layer and subsequent drag reduction. Additionally, it has been observed that the counteracting force of a polymer on a vortex is strongest immediately next to the vortices associated with regions of strong extensional flow (De Angelis, Casciola & Piva Reference De Angelis, Casciola and Piva2002; Sibilla & Baron Reference Sibilla and Baron2002; Terrapon et al. Reference Terrapon, Dubief, Moin, Shaqfeh and Lele2004; Dubief et al. Reference Dubief, Terrapon, White, Shaqfeh, Moin and Lele2005). A recent study by Han et al. (Reference Han, Kelly, Goldstein, Handler, Suryanarayanan and Burtsev2025) studied how small amounts of polymer affect the enstrophy of a simplified two-dimensional vortex and found a strong connection between enstrophy reduction and strain rate. These conclusions lead naturally to potential targeting schemes by which we can more effectively and efficiently reduce drag via polymer additives.
One practical method we explore is to emit polymer from simple particles, which can move differently than fluid elements, notably crossing streamlines. Wang & Maxey (Reference Wang and Maxey1993) showed that bubbles tend to accumulate significantly in the vortices produced by homogeneous isotropic turbulence, which means that light particles may naturally target coherent structures in complex flows. Other types of particles may also prove useful in targeting drag-producing structures and help elucidate the fundamental physics of turbulence. The literature suggests possible targets such as the buffer layer, high-strain regions or vortex structures themselves, which we explore in this work by using different types of particles and injection methods in a DNS study. In § 2, we introduce the governing equations for a viscoelastic fluid with non-uniform polymer distribution and the equations of motion for particles. In § 3.2, we first use a uniformly distributed polymer solution to verify that our code and set-up are accurate and provide a baseline for comparing the efficacy of the different targeting methods, including particle-based methods and highly controlled theoretical targeting to explore more fundamental results of polymer drag reduction. We then show how the targeting methods perform compared with an equivalent polymer ocean.
2. Methodology
2.1. Viscoelastic fluid equations
The continuity and momentum equations for an incompressible viscoelastic fluid are given by
\begin{equation} \dfrac {\partial u_i}{\partial x_i}=0 ,\end{equation}
and
\begin{equation} \dfrac {\partial u_i}{\partial t} + u_{\!j} \dfrac {\partial u_i}{\partial x_{\!j}} = -\frac {1}{\rho } \dfrac {\partial p}{\partial x_i} + \nu _s \frac {\partial ^2 u_i}{\partial x_{\!j} \partial x_{\!j}} + \frac {\partial \tau _{\textit{ij}}^p}{\partial x_{\!j}} + f_i, \end{equation}
where
$u_i$
are the fluid velocities,
$x_i$
are the coordinates,
$t$
is time,
$p$
is pressure,
$\rho$
is the mass density of the fluid,
$\nu _s$
is the kinematic viscosity of the solvent,
$\tau _{\textit{ij}}^p$
is the polymer stress tensor (per unit density) and
$f_i$
is the constant uniform pressure gradient which drives the flow. The model assumes that the polymer solution is dilute, and that the polymer stress, which appears in the last term of (2.2), is given according to the FENE-P model as follows:
\begin{equation} \tau _{\textit{ij}}^p = \frac {(\nu _0 - \nu _s)}{\lambda } [f(C_{\textit{kk}}) C_{\textit{ij}} - \delta _{\textit{ij}}]. \end{equation}
Here,
$\nu _0$
is the kinematic viscosity of the solution,
$\lambda$
is the polymer relaxation time,
$C_{\textit{ij}}$
is the polymer conformation tensor,
$f(C_{\textit{kk}})$
is the Peterlin function (defined in (2.6)) and
$\delta _{\textit{ij}}$
is the Kronecker delta. Note that, when there is no polymer added, the solvent and solution viscosities are identical and (2.2) is thus reduced to the Navier–Stokes equations. The conformation tensor evolves according to
\begin{equation} \dfrac {\partial C_{\textit{ij}}}{\partial t} + u_k \dfrac {\partial C_{\textit{ij}}}{\partial x_k} = C_{ik} \dfrac {\partial u_{\!j}}{\partial x_k} + C_{kj} \dfrac {\partial u_i}{\partial x_k} - \frac {1}{\lambda }[f(C_{\textit{kk}}) C_{\textit{ij}} - \delta _{\textit{ij}}] + \alpha _{p} \frac {\partial ^2 C_{\textit{ij}}}{\partial x_k \partial x_k}. \end{equation}
The last term represents an artificial diffusion of the conformation, which is added for numerical stability. Here,
$\alpha _{\!p}$
is a numerical diffusion coefficient whose only purpose is to stabilise the numerical scheme used for solving (2.4). This artificial diffusion term is necessary for spectral and pseudo-spectral methods using the FENE-P model, but it also fundamentally changes the detailed physics which the model represents (Xi Reference Xi2019; Dubief et al. Reference Dubief, Terrapon and Hof2023). Recent work by Beneitez, Page & Kerswell (Reference Beneitez, Page and Kerswell2023) and Couchman et al. (Reference Couchman, Beneitez, Page and Kerswell2024) showed that this diffusive term can even be the source of a physical instability that is not present in the diffusion-free system. However, the work of Sureshkumar & Beris (Reference Sureshkumar and Beris1995) has justified the use of
$Sc \equiv \nu _s/\alpha _{\!p} \lt 1$
, and under these conditions, the model still captures the drag-reduction mechanism of polymer additives quite well (Dubief et al. Reference Dubief, Terrapon and Hof2023). Generally, the best choice is the largest
$Sc$
that allows for a numerically stable simulation, which we determined through extensive testing. Several approaches to reduce the influence of this numerical diffusivity without sacrificing numerical stability have been proposed (e.g. Min, Yoo & Choi Reference Min, Yoo and Choi2001; Fattal & Kupferman Reference Fattal and Kupferman2005; Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019), but they all generally lead to higher computational costs. In this work, we define the Weissenberg number as
\begin{equation} {\textit{Wi}}_{\tau } \equiv \frac {\lambda u_\tau ^2}{\nu _s}, \end{equation}
where
$u_\tau$
is the friction velocity of the Newtonian turbulent channel. The Peterlin function
\begin{equation} f\!(C_{\textit{kk}}) = \frac {L_{\textit{max}}^2 - 3}{L_{\textit{max}}^2 - C_{\textit{kk}}}, \end{equation}
sets a limit on polymer stretching based on a specified maximum molecule length,
$L_{\textit{max}}$
, given in molecular diameters, which we choose to be
$L_{\textit{max}}=100$
.
For targeting polymer,
$\nu _0$
becomes a function of concentration,
$\gamma$
, whose evolution equation is given by
\begin{equation} \dfrac {\partial \gamma }{\partial t} + u_{\!j} \dfrac {\partial \gamma }{\partial x_{\!j}} = \alpha \dfrac {\partial ^2 \gamma }{\partial x_{\!j} \partial x_{\!j}} + s_{\gamma } \end{equation}
where
$\gamma$
is the concentration of polymer in parts per million (PPM) by weight,
$\alpha$
is the diffusion coefficient of polymer in the solvent and
$s_{\gamma }$
is a source term which depends on the targeting method chosen. Like
$\alpha _{\!p}$
, the value of
$\alpha$
is restricted by the numerical scheme used to solve (2.7). While it corresponds to the physical mass diffusivity of polymer in a solvent, the value is not necessarily representative of a particular solution’s physical mass diffusivity. Instead,
$\alpha$
is set to be equal to the numerical diffusion coefficient to ensure that spatial diffusion and numerical diffusion happen at the same rate. Equation (2.7) is a simplified description of polymer motion and inhomogeneity in the solution and notably does not account for other factors introduced through kinetic theory (see, e.g. Beris & Mavrantzas Reference Beris and Mavrantzas1994), but it is still descriptive of a dilute polymer solution in which molecules generally move with the flow and remains computationally tractable. In the FENE-P and Oldroyd-B models, the kinematic viscosity ratio,
$\beta \equiv \nu _s/\nu _0$
, is a function of concentration and is often used to describe how much polymer is added to the flow. In this study, we use an exponential model, based on the so-called Martin’s viscosity–concentration relationship (Martin Reference Martin1942; Weissberg, Simha & Rothman Reference Weissberg, Simha and Rothman1951), to calculate
$\beta$
based on scalar concentration
\begin{equation} \beta = \exp (-A_p \gamma ). \end{equation}
This particular relationship uses
$A_p = 2.418 \times 10^{-3}$
PPM
$^{-1}$
and is calculated to correspond to PEO dissolved in water based on the experimental work of Nsom & Latrache (Reference Nsom and Latrache2018).
2.2. Particle tracking equations
For two targeting cases presented, we use very small rigid particles that are moved based on the local flow according to a slightly modified version of the Maxey–Riley equations (Maxey & Riley Reference Maxey and Riley1983) with Faxén corrections
\begin{align} \dfrac {\text{d} {\boldsymbol{v}}}{\text{d} t} &= \varPi _{\!f} \bigg [ C_{\!D} \frac {3}{8a} \bigg |{\boldsymbol{u}} - {\boldsymbol{v}} + \frac {a^2}{6}{\nabla} ^2 {\boldsymbol{u}}\bigg | \bigg ({\boldsymbol{u}} - {\boldsymbol{v}} + \frac {a^2}{6}{\nabla} ^2 {\boldsymbol{u}}\bigg ) + \frac {3}{2} \dfrac {\text{D} {\boldsymbol{u}}}{\text{D} t} + \dfrac {\text{d} }{\text{d} t} \bigg (\frac {3a^2}{20} {\nabla} ^2 {\boldsymbol{u}}\bigg ) \nonumber \\&\quad + C_{\!L} \bigg ({\boldsymbol{u}} - {\boldsymbol{v}} + \frac {a^2}{6}{\nabla} ^2 {\boldsymbol{u}}\bigg ) \times \boldsymbol{\omega } \bigg ], \end{align}
where
$\boldsymbol{v}$
is the particle velocity,
$\boldsymbol{u}$
is the fluid velocity at the particle position
\begin{equation} \varPi _{\!f} \equiv \frac {\rho _{\!f}}{\rho _{\!p} + \dfrac {1}{2}\rho _{\!f}} ,\end{equation}
is a density ratio between the fluid (
$\rho _{\!f}$
) and particle (
$\rho _{\!p}$
);
$a$
is the particle radius;
$C_{\!D}$
is the drag coefficient of a rigid sphere based on the empirical relation first proposed by Schiller & Naumann (Reference Schiller and Naumann1935) and confirmed by Magnaudet, Rivero & Fabre (Reference Magnaudet, Rivero and Fabre1995)
\begin{equation} C_{\!D} = \frac {24}{\textit{Re}}\left(1 - 0.15 {\textit{Re}}^{-0.687}\right)\!, \quad {\textit{Re}} \leqslant 300, \end{equation}
where
${\textit{Re}} = 2a|{\boldsymbol{u}} - {\boldsymbol{v}}|/\nu _s$
is the particle Reynolds number;
$C_{\!L}$
is the lift coefficient for a sphere, taken to be
$1/2$
based on the work of Auton (Reference Auton1987);
$\mathrm{D}/\mathrm{D}t$
is the time derivative following a fluid element; and
$\boldsymbol{\omega }$
is the fluid vorticity at the centre of the particle. Two different particles are used in this study, referred to as bubbles and heavy particles. The bubble is a simulated particle with
$\rho _{\!p}/\rho _{\!f} = 0$
and the heavy particle is one with
$\rho _{\!p}/\rho _{\!f} = 10$
. Since these particles will carry some polymer inside them, their densities will ultimately depend on the amount and composition of polymer material, but for the purpose of studying potential targeting methods, we only examine these two limiting cases of an effective bubble and heavy particle in a flow of water. Both types of particles have a Stokes number of 0.95, given by
\begin{equation} \textit{St} \equiv \frac {2a^2}{9\nu _s \varPi _{\!f} \tau _\eta } ,\end{equation}
where
$\tau _\eta$
is the volume-averaged Kolmogorov time scale of the steady-state Newtonian turbulence. This Stokes number regime represents the largest possible bubbles we can simulate with (2.9) on our grid, allowing for distinctly different particle behaviour than tracers. In particular, we expect bubbles to be generally attracted to the vortex cores, and for the heavy particles to generally accumulate in regions of higher strain (Kelly et al. Reference Kelly, Goldstein, Suryanarayanan, Botto Tornielli and Handler2021; Wang & Maxey Reference Wang and Maxey1993).
The notable difference between the full Maxey–Riley equations and (2.9) is the lack of the buoyancy and Basset history terms. We are interested only in potential ways to transport polymer over a relatively short period of time, and we do not expect gravitational settling of particles to be a useful feature to capture. The Basset history term, while not always negligible (Prasath, Vasan & Govindarajan Reference Prasath, Vasan and Govindarajan2019), can be safely ignored under certain circumstances, such as when particles are heavy relative to the carrier fluid (Michaelides Reference Michaelides2006). For bubbles, the history term may be more important, but its overall impact on the overall drag-reduction results is expected to be minor since its effect on vortex entrapment has been shown to be minor (Tanga & Provenzale Reference Tanga and Provenzale1994). The most significant qualitative difference would be the reduction of preferential concentration (Daitche Reference Daitche2015), but this is already weak in our case. The particles used in this study are assumed to have negligible impact on the flow prior to releasing polymer and are therefore only one-way coupled to (2.2) due to their small size and low volume fraction (Marchioli Reference Marchioli2017; Zhai, Fairweather & Colombo Reference Zhai, Fairweather and Colombo2020). A further consideration would be to alter the drag coefficient based on viscoelastic effects of the polymer, which has been considered recently with high fidelity DNS (Faroughi et al. Reference Faroughi, Fernandes, Miguel Nóbrega and McKinley2020) and machine learning models (Faroughi, Roriz & Fernandes Reference Faroughi, Roriz and Fernandes2022). The correction factors for the Newtonian drag coefficient (taken to be
$C_{\!D} = 24/Re_p$
), are relatively small for a wide range of parameters, within 30 %. Certain extreme parameter cases can lead to correction factors up to approximately 6 or 7, but even that amount of correction has little impact on the overall drag reduction results presented in § 3.2, mostly because the bulk of particle motion happens before the particles are in the presence of any viscoelastic fluid.
2.3. Numerical method
To solve (2.1)–(2.7), we use a pseudo-spectral channel flow DNS code based on Kim, Moin & Moser (Reference Kim, Moin and Moser1987), Handler, Hendricks & Leighton (Reference Handler, Hendricks and Leighton1989) and Sureshkumar, Beris & Handler (Reference Sureshkumar, Beris and Handler1997). It uses a rectangular domain with evenly spaced grid points in the
$x$
- and
$z$
-directions and a Gauss–Lobatto grid in the
$y$
-direction. The
$x$
- and
$z$
-directions are periodic, while in the
$y$
-direction, no-slip boundary conditions are imposed for (2.1) and (2.2), and zero-flux conditions are imposed for (2.4) and (2.7). Equations (2.1)–(2.7) are solved by splitting them into linear and nonlinear components. The linear terms are solved in spectral space with a Crank–Nicholson scheme, and the nonlinear terms are solved in physical space using a second-order Adams–Bashforth integrator.
Unlike (2.1)–(2.7), (2.9) is only coupled one-way and therefore it can be solved independently. We solve for particle trajectories in physical space using a fourth-order Runge-Kutta integrator, with the particles programmed to bounce elastically off the top and bottom walls of the channel. The velocity, vorticity and velocity gradients are solved by DNS and then interpolated at each particle location to be used in (2.9). The interpolation is trilinear, which is sufficiently accurate when the grid is already well resolved for turbulence.
We simulate a fully turbulent three-dimensional channel flow at
${\textit{Re}}_{\tau } = u_\tau \delta /\nu _s = 180$
, with channel half-height,
$\delta$
, on a domain of size
$L_x \times L_y \times L_z =4\pi\! \delta \times 2\delta \times 2\pi\! \delta$
. The pseudo-spectral discretisation uses Fourier expansions in the periodic streamwise (
$x$
) and spanwise (
$z$
) directions, and Chebyshev polynomials in the wall-normal (
$y$
) direction, with a grid of
$N_x \times N_y \times N_z = 512 \times 129 \times 128$
. The
$3/2$
dealiasing rule was used in the
$x$
and
$z$
directions. Time is scaled by a large eddy turnover time, given by
$t^* \equiv {\textit{tu}}_\tau /\delta$
, with a time step of
$\Delta t^* = 0.9 \times 10^{-4}$
. The flow, which was initially a laminar Poiseuille flow, was tripped to turbulence using an obstruction generated via immersed boundary forces; the flow then develops into statistically stationary Newtonian turbulence with a constant pressure gradient. Eight snapshots of the steady-state turbulence, each separated by a few correlation times, are used as initial conditions for the various polymer additions, and those runs are ensemble averaged for each case. All non-dimensional parameters used in this study are listed and defined in table 1.
3. Results and discussion
3.1. Validation: polymer ocean
The motivation for using targeted polymer is that, for placing the same amount of polymer in a given domain, concentrating it in important regions will hopefully exhibit a greater effect overall than if the polymer were distributed uniformly (which we refer to as a polymer ocean). Therefore, identifying these regions and focusing polymer concentration there should boost the efficiency of polymer drag reduction. As a validation case, we show a polymer ocean with
$\beta = 0.9$
,
${\textit{Wi}}_{\tau } = 32.4$
and
$Sc = 0.7$
, which demonstrates that our model and code is reasonably accurate for predicting drag reduction.
Table 1. Important dimensionless quantities used throughout this study. Note that, since we enforce
$\alpha = \alpha _{\!p}$
, the Schmidt number can be defined equivalently with either term.
Qualitatively, we know that vortices are disrupted by the polymer and we see a greatly reduced vortex field after polymer is added and has affected the flow (Xi Reference Xi2009). Figure 1 demonstrates this qualitative effect of polymer on vortex structures in the channel for our polymer ocean case with
$\beta = 0.9$
. At the initial state of Newtonian turbulence (figure 1
a), there are many vortex structures identified by the Q-criterion. At the steady state of polymeric turbulence (figure 1
b), nearly all of the structures at that contour level have vanished, leaving only more dispersed and weaker vortex structures.
Figure 1. Vortex structures in (a) Newtonian turbulence and (b) viscoelastic turbulence. The vortex structures are drawn as iso-surfaces of the Q-criterion at 1000 s
$^{-1}$
, coloured by streamwise velocity. The snapshot in (a) is the initial turbulent state of a later time in (b) after polymer is added uniformly to the channel.
More importantly, however, is the extent to which drag is reduced. For a channel flow with constant pressure gradient, drag reduction is calculated based on the skin-friction coefficient,
$C_{\!f}$
, as
\begin{equation} \textrm{DR} = \frac {C_{\!f}^{N} - C_{\!f}}{C_{\!f}^{N}} \times 100\,\% ,\end{equation}
where
$C_{\!f}^{N}$
is the skin-friction coefficient of the Newtonian channel and
$C_{\!f}$
is the skin-friction coefficient of the channel with polymer additives. The skin-friction coefficient can be calculated as
\begin{equation} C_{\!f} = \frac {\tau _w}{\dfrac {1}{2}\rho U_b^2} ,\end{equation}
where
$\tau _w$
is the wall shear stress (WSS),
$\rho$
is the mass density of the fluid and
$U_b$
is the bulk velocity of the channel. For a viscoelastic fluid, we must account not only for the shear stress of the fluid, but also the shear stress of the polymer, so
\begin{equation} \tau _w \equiv \left .\left(\mu _s \dfrac {\partial u}{\partial y} + \tau ^p_{\textit{xy}}\right)\right|_{y = \pm \delta }\! ,\end{equation}
where
$\mu _{s}$
is the dynamic viscosity of the solvent, and
$\tau ^p_{\textit{xy}}$
is the polymer shear stress. Thus, in general, the reduction in drag is based on total shear stress at the wall and bulk velocity, both of which are captured in figure 2.
Figure 2. (a) Bulk velocity of the polymer ocean relative to Newtonian when polymer is added abruptly and uniformly at
$t^* = 0$
. (b) Total WSS normalised by average Newtonian WSS.
We ran eight independent simulations of the polymer ocean, and for each one, the polymer is turned on at time
$t^* = 0$
, when the channel flow is statistically stationary Newtonian turbulence. The results presented throughout this study correspond to the ensemble-averaged values of all eight simulations, further averaged in space and/or time as specified. Each simulation used a unique initial state that was uncorrelated with any other. The flow was then allowed to evolve according to (2.1)–(2.4) until it reached a quasi-steady state where the disrupted turbulence regeneration cycle reaches a new equilibrium. At this point, we average the flow over an extended time to compare with the Newtonian case to determine drag reduction. In the long-time quasi-steady state of polymeric turbulence, the balance of the driving pressure gradient and the total shear stress on the walls necessitates an equilibrium equivalent to that of the Newtonian turbulence. Thus for long times, the drag reduction may be calculated as
\begin{equation} \textrm{DR} = \bigg [ 1 - \bigg (\frac {U_b^N}{U_b}\bigg )^2 \bigg ] \times 100\,\% . \end{equation}
The average bulk velocity ratio over
$t^*=70$
to
$t^*=90$
is 1.246, which translates to a drag reduction of
$35.4\,\%$
$(\pm 0.5\,\%)$
. This value agrees well with several reports in the literature, including the DNS study by Min, Yoo & Choi (Reference Min, Yoo and Choi2003) where a channel flow with
${\textit{Re}}_\tau = 135$
,
${\textit{Wi}}_{\tau } = 32.4$
and
$\beta = 0.9$
can be inferred to achieve approximately 40 % DR; and a similar case by Dubief et al. (Reference Dubief, White, Terrapon, Shaqfeh, Moin and Lele2004) with
${\textit{Re}}_{\tau } = 300$
,
${\textit{Wi}}_\tau = 36$
and
$\beta = 0.9$
achieved 35 % DR. Li, Sureshkumar & Khomami (Reference Li, Sureshkumar and Khomami2006) even provides a DR % scaling equation, which for our
${\textit{Re}}_\tau$
,
${\textit{Wi}}_\tau$
and
$L_{\textit{max}}$
estimates 33.9 % DR. Additionally, we have used the pseudo-spectral partial differential equation solver, Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), to independently set up a turbulent channel with uniform polymer added instantaneously with identical features to the current problem. We found that an average steady-state bulk velocity ratio settled around 1.27, giving a reduction in drag of 38.2 %, which is reasonable agreement with our result. The moderate DR % value puts the current simulations in an intermediate DR regime. Low drag-reduction (LDR) flows tend to have DR % values around 35 % and below, whereas high drag-reduction (HDR) flows tend to have much higher values, generally above 38 % (Kim et al. Reference Kim, Li, Sureshkumar, Balachandar and Adrian2007; Xi Reference Xi2019). However, other characteristics of the flow align well with the LDR regime, as shown later.
In figure 2(b), we see that there is initially a large spike the WSS corresponding to the stretching of polymer molecules and a corresponding increase in polymer shear stress, which appears as a small decrease in relative mass flux. Soon after, however, the turbulence regeneration cycle starts breaking down and we see a sharp decrease in wall shear stress and a corresponding increase in the mass flux. At some time around
$t^* = 50$
, the flow starts to stabilise in a new viscoelastic turbulence state where the WSS and mass flux fluctuate about a steady average value. Notably, the WSS fluctuations are significantly larger than the Newtonian fluctuations, most likely due to the polymer effects increasing dynamical complexity of the system (Graham & Floryan Reference Graham and Floryan2021).
Figure 3 shows the mean velocity profile averaged over the portion of the simulation when the flow is at a steady state. The Newtonian profile agrees very well with the DNS statistics of
${\textit{Re}}_\tau = 180$
channel flow of Kim et al. (Reference Kim, Moin and Moser1987), and matches the log law with the equation
$U^+(y^+) = 2.5 \ln y^+ + 5.5$
. In the polymer ocean profile, we observe the characteristic thickening of the buffer layer and lifting up of the log layer profile (Lumley Reference Lumley1973), but the flow profile does not reach the MDR asymptote given by
$u^+ = 11.7 \ln y^+ - 17$
(Virk Reference Virk1975), which is expected for LDR flows.
Figure 3. Mean velocity profile,
$u^+$
vs
$y^+$
, for a polymer ocean (blue) with
$\beta =0.9$
and
${\textit{Wi}}_\tau = 32.4$
compared with a Newtonian flow (red). The mean flow profile is averaged from
$t^* = 70$
to
$t^* = 90$
.
Figure 4. Shear stress profiles for a polymer ocean (
$-$
) with
$\beta =0.9$
and
${\textit{Wi}}_\tau = 32.4$
and a Newtonian flow (
${-}{-}$
). The black lines represent the total shear stress profile, the blue lines represent the Reynolds shear stress, the green lines represent the mean shear stress and the red lines represent the polymer shear stress.
Figure 5. The r.m.s. fluctuating velocity profiles for a polymer ocean (
$-$
) with
$\beta =0.9$
and
${\textit{Wi}}_\tau = 32.4$
, and a Newtonian flow (
${-}{-}$
).
Figure 4 shows the average shear stress profile, whose magnitude is averaged over both halves of the channel and time averaged from
$t^* = 70$
to
$t^* = 90$
. The total shear stress is computed from the combination of the mean shear stress, Reynolds shear stress and polymer shear stress
\begin{equation} T_{\textit{xy}} = \mu \dfrac {\text{d} u}{\text{d} y} - \rho \overline {u'\!v'} + \tau _{\textit{xy}}^p. \end{equation}
For the Newtonian flow,
$\tau _{\textit{xy}}^p\equiv 0$
. At a steady state, we see that the total shear stress profiles are nearly identical, but the components change significantly. The Reynolds shear stress (RSS) profile significantly decreases across the channel, especially in the buffer layer, and the peak shifts away from the wall, as expected (Min et al. Reference Min, Yoo and Choi2003; Wang, Shekar & Graham Reference Wang, Shekar and Graham2017; Serafini et al. Reference Serafini, Battista, Gualtieri and Casciola2023). The mean shear stress for the polymer ocean is smaller than the Newtonian case near the wall but larger as we look away from the wall, which seems to be aligned with the mean flow profile in figure 3; the polymer shear stress produces a force that opposes the flow strongly near the walls, which widens the viscous sublayer slightly, reducing the slope prior to
$y^+=10$
. Beyond this point, the flow drastically increases in the thickened buffer layer, which produces a velocity profile with a larger slope from
$y^+=10$
to the centre of the channel. The polymer stress is largest near the wall, where the velocity gradients there cause maximal polymer stretching and consequently increase the polymer shear stress (Mortimer & Fairweather Reference Mortimer and Fairweather2022).
Figure 5 shows the root-mean-square (r.m.s.) fluctuating velocity profiles of the polymer ocean and Newtonian flows. The streamwise velocity fluctuations increase moderately away from the wall with a peak shifted to higher in the buffer layer, while the spanwise and wall-normal velocity fluctuations decrease moderately throughout the channel. These profiles are consistent with the LDR regime (Xi Reference Xi2019), although a similar trend can be seen in HDR flows with rigid polymers (Mohammadtabar, Sanders & Ghaemi Reference Mohammadtabar, Sanders and Ghaemi2017).
3.2. Targeted polymer
We examine five targeted injection cases: two in which physical particles act as sources of polymer, two in which polymer is artificially injected to target regions in or around vortex structures and one to simulate a case of polymer being injected directly into the buffer layer, near the walls. These cases start out distinctly different from one another and the polymer ocean but eventually each case will end at the same state due to diffusion and mixing; that is, all targeting cases revert to a polymer ocean after sufficient time has elapsed. Therefore, the key time period to examine is early times relative to the eventual asymptotic limit of drag reduction. This is also relevant for physical applications with a finite length of wall over which the viscoelastic flow may not reach a steady state by the time it moves beyond the end of the wall, or for some other limiting factor of time such as polymer degradation due to scission.
Figure 6. Initial positioning of particles near the wall with vortex structures shown by iso-surfaces of Q-criterion at 1000 s
$^{-1}$
coloured by streamwise velocity. Only a portion of the channel is shown, but the particles are distributed the same across the entirety of both walls of the channel.
Since the duration of the release is long compared with the time it takes for the polymer to affect the flow, the targeting cases are compared with a polymer ocean which is gradually increasing in concentration for a more fair comparison. The gradual polymer ocean case increases in concentration uniformly and steadily from a state of
$\beta = 1$
at
$t^* = 0$
until it reaches a concentration such that
$\beta = 0.95$
, which corresponds to a dilute polymer concentration of approximately 21 PPM by weight. This is roughly half the concentration of
$\beta =0.9$
(approximately 43 PPM by weight) used in the polymer ocean validation study. The rate at which polymer concentration increases is set such that it corresponds to the time it takes for any given targeting method to reach the final average concentration, which happens around
$t^* = 7.9$
.
Just as in § 3.1, eight independent simulations were run and subsequently ensemble averaged to produce the results presented for these new polymer addition methods. These different methods are compared up to
$\overline {\beta } = 0.95$
with
${\textit{Wi}}_\tau = 32.4$
, where
$\overline {\beta }$
refers to the spatially averaged value of
$\beta$
of the final state of polymer, although this is somewhat larger than what is typically used in polymer drag-reduction studies. For example, the simulations of Sureshkumar et al. (Reference Sureshkumar, Beris and Handler1997), Min et al. (Reference Min, Yoo and Choi2003), Dubief et al. (Reference Dubief, White, Terrapon, Shaqfeh, Moin and Lele2004) and Mortimer & Fairweather (Reference Mortimer and Fairweather2022) use a value of
$\beta = 0.9$
under the FENE-P or Oldroyd-B models, and most of the experiments reported in Lumley (Reference Lumley1969), Virk (Reference Virk1975) and Xi (Reference Xi2019) use polymer concentrations of
$O(10)$
to
$O(100)$
PPM. Using a lower concentration is necessary to avoid numerical instabilities for the targeting methods which are much more sensitive than a polymer ocean due to gradients in concentration and stretching, but the current study also assumes a monodisperse polymer solution. Much of the prior experimental work deals with polydisperse polymer solutions, which have a wide range of molecular weight polymer chains, and a higher concentration is necessary to compensate for this polydispersity (Brandfellner et al. Reference Brandfellner, Muratspahić, Bismarck and Müller2024; Serafini et al. Reference Serafini, Battista, Gualtieri and Casciola2025). The Schmidt number is increased from 0.7 in the validation study to 1.0 here in order to provide more robust and distinct results between targeting cases. If the Schmidt number is too small, then the polymer will diffuse quickly and render the targeting methods indistinct after only a short time.
Figure 7. Visualisation of targeting methods. (a) Bubbles (black dots) releasing polymer (green contours); (b) injection of polymer targeting regions of high Q-criterion; and (c) near-wall injection of polymer. Only a quarter of the channel is shown (
$[0,2\pi\! \delta ] \times [-\delta ,0] \times [0,2\pi\! \delta ]$
). Vortex structures are shown as iso-surfaces of the Q-criterion at 1000 s
$^{-1}$
, coloured by streamwise velocity.
Figure 6 shows the initial seeding of particles near the walls for the particle targeting methods. Bubbles and heavy particles are placed near both walls at
$y^+=18$
, then move passively in a Newtonian flow for a certain duration to allow time for preferential accumulation. These targeting methods will be referred to by their respective particle types. The particles then begin releasing polymer at a constant rate until a spatially averaged concentration is reached such that
$\beta = 0.95$
. A snapshot of this process for a case of bubbles is shown in figure 7(a). If we distribute this amount of polymer across all particles, this would give an actual density ratio of bubbles of
$\rho _{\!p}/\rho _{\!f} = 0.037$
and for heavy particles,
$\rho _{\!p}/\rho _{\!f} = 6.1$
, assuming PEO is the only material encased in a virtually massless shell. The size of particles differs based on density such that they both have the same Stokes number. The near-wall injection method is similar to the particle methods in that the polymer is injected from point sources located at the initial particle positions, however, the location of the point sources remains fixed with time. This results in the polymer being emitted essentially as a sheet across the domain, which is distinctly different than the experiments of McComb & Rabie (Reference McComb and Rabie1982a
,Reference McComb and Rabie
b
), although the height of injection is similar. The artificial targeted injection methods are functionally the same and are implemented at each time step as follows:
-
(i) calculate the field in physical space;
-
(ii) identify all grid points with discrete local maxima in the field;
-
(iii) sort the local maxima by field value;
-
(iv) insert polymer at a fixed rate in the top
$N$
local maxima, skipping regions where polymer concentration has already reached a specified maximum;
where ‘field’ refers to either the Q-criterion or the sum of squares of the axial strain rates
\begin{equation} \sum _{k=1}^3 S_{\textit{kk}}^2 = \bigg (\dfrac {\partial u_1}{\partial x_1}\bigg )^2 + \bigg (\dfrac {\partial u_2}{\partial x_2}\bigg )^2 + \bigg (\dfrac {\partial u_3}{\partial x_3}\bigg )^2 .\end{equation}
The goal of this algorithm is to put polymer in regions of the highest value of some field variable but distribute it smoothly in small volumes around several local maxima without pouring too much polymer into one particular region. For further explanation of the targeting algorithm, see Appendix A. These targeting methods are referred to as high-Q targeting and high-strain targeting. Extensive testing revealed that choosing only the high normal strain rate regions as targets was more effective at identifying regions around vortices than choosing regions of high total strain (characterised by the Frobenius norm of the total strain rate tensor) or large negative Q-criterion. Targeting these more general regions resulted in most of the polymer being injected near the walls and was virtually the same as the near-wall injection method. Therefore, the term high-strain targeting refers only to the case of targeting high axial strains. This targeting method is partially inspired by Terrapon et al. (Reference Terrapon, Dubief, Moin, Shaqfeh and Lele2004), who found that polymer stretching is associated with regions of strong biaxial extensional flow, which is encompassed in (3.6). An example of high-Q targeting is depicted in figure 7(b). For consistency with the particle targeting, we chose to target
$N = 1024$
cells and set an arbitrary limit of local
$\beta = 0.8$
to be the threshold for skipping a targeted cell candidate. For example, if a targeted cell already has a
$\beta$
value at or below 0.8, we do not add more polymer there and choose a different cell to target. Local maxima are chosen in order to avoid too much overlap in a region of particularly large field values and to spread out the targeting to affect as much of the domain as possible. Similarly, the threshold is in place to limit the local concentration of polymer which could lead to both numerical and physical issues. Since the FENE-P model is only accurate for dilute polymer solutions, unchecked polymer addition may incorrectly represent highly elastic regions, and large gradients in polymer concentration can lead to numerical instabilities. In most cases, however, this limit was merely a precaution and did not affect the results of the polymer addition.
All targeting cases inject polymer into the domain at the same rate via a Gaussian distribution of polymer centred at either the particle location or the calculated injection cell. The width of the Gaussian is chosen to be slightly larger than the particle diameter to closely mimic material oozing from a particulate carrier or dissolving from the surface. Due to the proximity of some injection locations to the wall, the profiles may be cut off and the injection rate at that particle or grid cell is thus reduced. Therefore, the specific times for the global polymer concentration to reach the cutoff threshold is slightly different from case to case, but only by a few time steps. The injection rate of polymer is chosen such that the time it takes to reach the global concentration threshold is long compared with flow time scales.
We can examine how the polymer distribution changes based on the various targeting methods by averaging the concentration across planes parallel to the wall. Using the average concentration, we can compute a local,
$y$
-dependent
$\beta$
value using (2.8). Figure 8 shows how this distribution of polymer aligns with the distribution of vortical structures, educed by the swirling strength criterion,
$\lambda _{ci}$
, for its convenient property of being entirely non-negative, and positive
$\lambda _{ci}$
values are kinematically identical to positive Q-criterion values (Jeong & Hussain Reference Jeong and Hussain1995). The distributions are shown for values at
$t^* = 8$
, so this is well after the vortex structures have been disrupted and diminished and it is immediately after polymer addition stops. We notice in figure 8(a) that the high-strain and high-Q targeting methods produce similar polymer distribution shapes to that of the
$\lambda _{ci}$
distribution, as expected, and that the near-wall injection has much higher concentration at and near the walls. The two particle cases are shown in figure 8(b) at the same time, and we notice that the heavy particles exhibit virtually no targeting, being very closely aligned with the horizontal polymer ocean profile. The bubbles are only slightly better but with more polymer being injected in the centre of the channel as the bubbles tend to be driven away from the walls by the turbulent motions of the flow.
Figure 8. Average swirl criterion across the channel (left axis, blue) compared with polymer distribution (right axis, orange) at
$t^* = 8$
for (a) high-strain targeting and near-wall injection and (b) bubbles and heavy particles compared with the polymer ocean.
Since the targeting methods are compared at earlier times, the drag reduction will be driven both by the shear stress at the wall as well as the bulk velocity of the flow. The relative mass flux (equivalent to the relative bulk velocity for identical geometry) is shown in figure 9(a), and it is illuminating for how the flow is changing. At long times we will see the total WSS equilibrate at the Newtonian value, but the mass flux will generally increase until it reaches a steady level. At early time, we see that for all cases of polymer addition the mass flux does not increase until approximately
$t^* = 1$
. It takes time for the polymer to affect drag-producing structures and for this in turn to be manifested as an increase in the mass flux of the channel. The most dramatic increase in mass flux is after this point, and it subsequently increases more slowly until it eventually levels off around a steady state. For the current runs, this happens around
$t^* = 50$
, just as for the polymer ocean with
$\beta = 0.9$
, which is well after times shown in figure 9. At
$\beta = 0.95$
, the high-strain and high-Q targeting cases increase more rapidly throughout the simulation and reach the steady state sooner than the polymer ocean, suggesting that targeting is in fact more effective throughout the drag-reduction process even though all the polymer addition cases eventually end up at the same level. The high-strain and high-Q targeting methods have larger increases in mass flux throughout the time shown than the other methods. The particle-based methods are virtually the same as the polymer ocean case, and examining the distribution of polymer during the simulation showed that neither of these particulate cases exhibited significant clustering of particles in vortical regions; rather, they were largely dispersed throughout the domain. Notably, the near-wall injection method exhibits distinctly lower mass flux than the other methods, which is due to the increased polymer stresses very close to the wall. The near-wall injection has the highest polymer concentration at the walls by far and therefore behaves differently than the other methods, as we see in figure 9.
Figure 9. (a) Relative mass flux and (b) WSS of polymer addition methods for
$\overline {\beta } = 0.95$
.
The ensemble-averaged magnitude of total shear stress at the walls,
$\langle \tau _w \rangle$
, is shown in figure 9(b), normalised by the time-averaged Newtonian WSS,
$\overline {\tau }_{w,N}$
. We see a brief period where polymer is initially stretched and the shear stress remains near the Newtonian case, but then we see a drastic decrease in the total WSS, most significantly by the two structure targeting injection methods. The delay for reducing the wall shear stress corresponds to the delay of the increasing mass flux, happening around
$t^* = 1$
. The near-wall injection exhibits a much larger spike in WSS than the other five cases, which also aligns with the increased polymer stress and high concentration near the walls. The high-strain and high-Q targeting methods have the most significant WSS reduction, but all the cases tend slowly upward until they ultimately fluctuate around the Newtonian WSS value due to the force balance with the constant pressure gradient (at around
$t^*=50$
, not shown in figure 9
b).
Figure 10. Mean velocity profile,
$u^+$
vs
$y^+$
, for high-strain targeting, polymer ocean, near-wall injection and Newtonian at
$\beta = 0.95$
averaged over
$t^* = 8$
to
$t^* = 12$
.
In figure 10, the mean flow profile is plotted for the high-strain targeting, polymer ocean and near-wall injection cases along with the Newtonian case. Since the flow profiles are so close together, only these three polymer cases are shown. The velocity profiles are scaled by the Newtonian inner scaling variables, and they are generated from the mean velocity profile averaged from
$t^* = 8$
to
$t^* = 12$
. This time period is after all polymer has been added and when the mass flux increase starts to diminish. Here, we already see the characteristic thickening of the buffer layer and a lifting up of the profile in the log layer for both cases with polymer additives, although there is not much distinction between the different configurations. It should be noted that this mean profile does not reflect the steady state, since, as shown in figure 2(a), the mass flow continues to develop until well after the initial addition of polymer.
Figure 11. Shear stress components for high-strain targeting (
$-$
), polymer ocean (
${-}{-}$
), near-wall injection (
$-\boldsymbol{\cdot }-$
) and Newtonian (
$\boldsymbol{\cdots}$
) at
$\beta = 0.95$
averaged over
$t^* = 8$
to
$t^* = 12$
.
In figure 11, the transient shear stress profiles are shown, again only for the high-strain targeting, polymer ocean and near-wall injection cases, compared with the Newtonian flow. Just like in § 3.1, we see a significant reduction in RSS throughout the domain, especially in the buffer layer, which agrees with prior observations (Zhu et al. Reference Zhu, Schrobsdorff, Schneider and Xi2018; Xi Reference Xi2019), although in the averaged transient state, the shear stress profiles are uniformly lower than what is shown for the
$\beta =0.9$
polymer ocean case in figure 4. The near-wall injection case generally has the highest RSS profile of all the polymer additive methods, which corresponds to the lowest DR % across the cases; however, the polymer ocean case demonstrates the most RSS reduction despite the high-strain targeting method outperforming it. Interestingly, the mean shear for the near-wall injection case is slightly lower across the domain, but the polymeric shear stress component is much higher. This suggests that the higher concentration near the walls leads to a much more significant effect of polymer stretching, and since the region where the polymer concentration is largest is not aligned well with the location of most of the drag-producing structures, the mass flux is lower, the WSS is higher and the overall drag reduction is much less effective than for the other methods. Figure 12 shows the r.m.s. fluctuating velocity profiles for the same methods, and all three polymer methods are very similar. However, it is interesting to note that these profiles are all similar to that of the polymer ocean with
$\beta = 0.9$
in figure 5, indicating that the fluctuating velocities are changed very rapidly by the polymer addition regardless of how it is added.
Figure 12. The r.m.s. fluctuating velocity profiles for high-strain targeting (
$-$
), polymer ocean (
${-}{-}$
), near-wall injection (
$-\boldsymbol{\cdot }-$
) and Newtonian (
$\boldsymbol{\cdots}$
) at
$\beta = 0.95$
averaged over
$t^* = 8$
to
$t^* = 12$
.
The mass flux and WSS together determine the overall drag reduction, but a unique challenge for this study is capturing the drag-reduction performance for a transient flow. To quantify this in a holistic way, we use a cumulative measure of the drag reduction through a continuous average
\begin{equation} \overline {\textrm{DR}}(t) = \frac {1}{t} \int _{0}^{t} \bigg (\frac {\overline {C}_{\!f}^{N} - C_{\!f}(\tau )}{\overline {C}_{\!f}^{N}}\bigg ) \, {\rm d}\tau \times 100 \,\% ,\end{equation}
Table 2. Cumulative drag-reduction comparison for various targeting methods of polymer addition at early times.
Figure 13. Continuously averaged drag reduction over time for all targeting cases at
$\overline {\beta } = 0.95$
.
Figure 14. Mean swirling strength, spatially averaged over the entire channel, for all polymer addition methods compared with a Newtonian channel.
where
$\overline {C}_{\!f}^{N}$
is the mean skin-friction coefficient for the Newtonian flow, averaged over the entire simulation. It can also be interpreted as the normalised relative impulse due to polymer additives which corresponds to a relative change in momentum of the fluid due to drag reduction (see Appendix B). This time-averaged drag reduction provides a smoother, more representative picture of the overall performance of the polymer additives over the entire transient period. The results reported in table 2 correspond to the drag-reduction values at a few key points during the transient state, and the drag reduction over the entire duration is shown in figure 13. Based on the trends in this figure, we can distinguish four stages of drag reduction. The first stage is the brief oscillation from
$t^* = 0$
to
$t^* \approx 1$
where the addition of polymer is small and has not disrupted the turbulence regeneration process. In this window, the DR plot is very sensitive to small changes and the initial discrepancies between the instantaneous start condition vs the time-averaged values used for normalisation. After this point, we see the onset of drag reduction, which rapidly increases as more polymer is added in the channel and it starts to significantly disrupt drag-producing structures. This stage, from
$t^* \approx 1$
to
$t^* \approx 5$
is where we see the primary differences in the effects of the targeting methods. The injections targeting high-strain regions and Q-criterion enter this regime more quickly than the other cases and also have the largest slopes, suggesting that both of these methods are more effective at disrupting the drag-producing structures in the flow. This region is also where the bulk of vortex disruption happens, as is evident in figure 14, which shows the spatially averaged swirl criterion over time. By
$t^* = 4$
–5, the mean swirl throughout the domain has settled on a particular value for all polymer addition cases, which means that beyond this point the structures are no longer being significantly disrupted despite more polymer being added. While there are vortex structures which persist in the viscoelastic turbulence, they are different in strength and shape (Wang et al. Reference Wang, Shekar and Graham2017), so the targeting methods may not be as effective beyond this point with further distinction between targeting methods being determined primarily by local polymer concentration prior to this point. The third stage, around
$t^* = 5$
to
$t^* = 50$
, is where we start to see diminishing returns, and the slope of the cumulative drag reduction levels off. Once the polymer alters the structure of the turbulence, the relative effectiveness will naturally decrease since the structures it affects have diminished and the flow is settling into a new equilibrium. After the polymer is no longer injected into the domain, it can only move under the influence of advection and diffusion, ultimately leading to uniformly distributed polymer at long times. Finally, the fluid reaches a new quasi-steady state at some asymptotic value of drag reduction, and all targeting cases become indistinguishable (at some point around and after
$t^* = 50$
, not shown here but may be inferred from figure 2).
Table 2 shows the values for the targeting methods at three particular times: when the drag reduction is finishing its initial rapid increase (
$t^* = 4$
), when the polymer is finished being added (
$t^* = 8$
) and when the drag reduction is diminishing (
$t^*=12$
). The long-term drag reduction settles around 32 % for all cases at long times, so the values reported in table 2 capture the overall effect quite well. Across the board, the near-wall injection performs the worst, which is surprising given the success of McComb & Rabie (Reference McComb and Rabie1982a
,
Reference McComb and Rabieb
). However, it must be noted that the simulated method is distinctly different from their experiments in a few different ways. Aside from the difference in geometry (pipe flow vs channel flow), the experimental set-up used a single injector whereas the current simulations flood the near-wall regions from distributed point sources. Additionally, their experiments involved much higher concentrations of polymer and they achieved up to 69 % DR, which is a HDR flow where the dynamics may be different. The effective range for their injection was lower bounded at
$y^+=15$
, which is right around the centre of polymer injection in the present simulations, leading to significant local polymer concentrations even closer to the wall. The result for the current simulations is a much higher near-wall concentration and increased WSS due to its significantly larger polymeric shear stress. However, this generally only sets the performance back initially by delaying the onset of drag reduction; we see that the increase in DR throughout the simulation for the near-wall injection is on par with all the other cases. On the other hand, both high-Q and high-strain targeting methods perform consistently above the rest until overall performance starts diminishing, when diffusion and mixing make the targeting methods less distinct. Together, these results show a significant disparity between a ‘poor’ targeting method (injecting directly near or at the walls) and a ‘good’ targeting method (injecting in or around the vortical structures), which provides strong evidence that, even in the LDR regime, placement of polymer can be an effective way to improve drag-reduction performance.
4. Conclusions and future work
In this study, we simulated various methods of targeted polymer injection into a Newtonian turbulent channel to examine how distribution of polymer affects drag reduction. At long times, the polymer will be uniformly distributed throughout the domain and there is no significant difference between any method of targeting. In the transient state of polymer drag reduction, however, we see a consistent trend of targeting methods reducing the drag more effectively than a polymer ocean. Using light and heavy particles as carriers for polymer to be released into the flow appeared to produce similar results as a polymer ocean due to modest accumulation. The most effective methods proved to be targeting coherent vortex structures in the flow, both by injecting polymer at local maxima of the Q-criterion (targeting the centres of those structures) and injecting polymer at local maxima of the normal strain rate (targeting regions which tend to be just outside of the vortex structures), which increased drag reduction up to 2.5 % compared with the polymer ocean. The near-wall injection displayed the worst performance since this method produced very high polymer shear stresses at the beginning of the simulation and little effect directly on the primary drag-producing structures. These results provide a theoretical foundation for improved efficiency of polymer drag reduction when selectively adding polymer to certain regions of the flow. Conversely, the poorer performance of near-wall injection implies that there are ‘good’ and ‘poor’ methods of targeted polymer addition, which demonstrates that a uniformly mixed polymer solution is not necessarily ideal for polymer-based drag reduction and offers insight into the fundamental mechanism of polymer drag reduction which apparently hinges on interactions with coherent fluid structures. Indeed, between the best and worst cases, we saw a gap of roughly 4 %–5 % drag reduction. We acknowledge that this is a modest difference within the context of polymer drag reduction which can practically see up to 80 % DR in certain cases (Lumley Reference Lumley1969; Virk Reference Virk1975). We also note that the targeting method used to achieve the highest possible drag reduction is not physically viable with current technological limitations. Nevertheless, we believe that further exploration of this approach is warranted to determine its viability for large-scale polymer drag reduction for marine vehicles with historically limiting ecological and economical concerns.
Since targeted polymer drag reduction has relatively little emphasis in the literature, there are still many avenues to explore to further refine our understanding and develop new targeting techniques. For example, the simulated injection of polymer does not account for changes in mass or momentum; while it can be argued that these should be negligible for such dilute concentrations, it is still a simplifying assumption that may have an impact on physical reproductions of these techniques. Similarly, the lack of thermodynamic closure presents a problem for the accuracy of the current model and may be improved to better simulate reality. This study, while promising, may not fully demonstrate the benefits of targeted polymer drag reduction due to model and simulation limitations, nor does it thoroughly examine the fluctuating polymer-induced forces and torques at play. We chose a few specific ways to target based on hypotheses and our understanding of the mechanisms involved in polymer drag reduction, but there may be various other configurations with larger reductions in drag , or even larger differences in DR between polymer addition methods. The current framework could also serve as a useful platform for more comprehensive studies of polymer–turbulence interactions. In particular, targeting strategies like high-Q and high-strain injection naturally isolate zones of strong velocity gradients, where polymer stress gradients are likely to produce localised body forces. Future work could leverage these targeted methods to explore how such forces influence vortex suppression and the buffer-layer dynamics. Future development of targeting methods should test with much larger Weissenberg numbers for physical relevance as well as removing the dependence on
${\textit{Re}}_{\tau }$
(Serafini et al. Reference Serafini, Battista, Gualtieri and Casciola2024). However, the FENE-P model is not suited for that regime and the simulations would therefore require a more robust and physically representative polymer model such as the direct simulation of polymer molecules in a hybrid approach. This would also provide more freedom to lower
$\overline {\beta }$
and increase the Schmidt number without numerical instabilities plaguing the simulations, allowing us to understand how these results translate to a HDR regime. Lastly, the dependence of drag reduction on the release rate for the targeting methods has only been modestly explored due to numerical stability issues, and further testing may demonstrate greater performance improvements than what have been shown in this study.
Funding
The authors gratefully acknowledge funding from the NSF, awards 1905288 and 1904953.
Declaration of interests
The authors report no conflict of interest.
Appendix A.
The structure targeting algorithm was constructed carefully to be comparable to particle addition methods and to ensure it reliably targeted structures throughout the time of polymer addition, regardless of the duration. The algorithm is presented below as a pseudocode example of targeting large values of the Q-criterion. In the case of high-strain targeting, the only difference is the first function call, which computes whatever field value we are interested in. The details of certain subroutines are omitted for readability, although their purposes should be clear. The test for local minimum is a simple search across grid points. For a well-resolved simulation, the grid cell sizes are much smaller than the structures, so this should sufficiently target certain flow regions. The Qmin array is not a consistent size since the total number of local minima in the domain may vary with each time step. It is only used to sort the index vectors using a shell-sorting algorithm. The vectors Qx, Qy and Qz contain the indices of the corresponding local minima, and these are used to subsequently reference specific points in the overall field. The algorithm starts at the maximum field value and first checks if the local value of
$\beta$
is higher than 0.8. If it is lower than this threshold value, that particular local maximum is skipped to prevent overloading the polymer at certain regions. From testing, the choice of this value had little effect on the overall results and is more of a precaution than a necessity.
Appendix B.
The balance of forces integrated over the volume of the domain is given by
\begin{equation} m \dfrac {\text{D} U_b}{\text{D} t} = F_{\!P} + F_{\!f} , \end{equation}
where
$m$
is the mass of fluid in the domain,
$F_{\!P}$
is the force due to pressure and
$F_{\!f}$
is the friction force due to shear stress. For a viscoelastic fluid, the latter term will encompass both fluid shear stress and polymer shear stress. If we normalise (B1) by the dynamic pressure
$({1}/{2})\rho U_b^2$
, we get the force balance in terms of the skin-friction coefficient,
$C_{\!f}$
\begin{equation} \frac {2 m}{\rho U_b^2} \dfrac {\text{D} U_b}{\text{D} t} = \frac {2 F_{\!P}}{\rho U_b^2} + C_{\!f}. \end{equation}
Integrating over time
$t$
, we get the impulse–momentum relationship
\begin{equation} \Delta \tilde {p} = I_{\!P} + I_{\!f}, \end{equation}
where
$\Delta \tilde {p}$
is the change in normalised momentum,
$I_{\!P}$
is the impulse of the normalised pressure force and
$I_{\!f}$
is the impulse of the normalised friction force. Therefore, the drag reduction at time
$t$
may be interpreted as the relative change in normalised momentum due to the change in the friction force, expressed via the impulse
\begin{equation} \textrm{DR}(t) = \frac {I_{\!f}^{\!N} - I_{\!f}^{\!p}(t)}{I_{\!f}^{\!N}} \times 100\,\%, \end{equation}
where the superscripts
$N$
and
$p$
refer to the Newtonian case and polymer case, respectively. Writing this out in terms of
$C_{\!f}$
gives us
\begin{equation} \textrm{DR}(t) = \frac {\int _{0}^{t} (\overline {C}_{\!f}^N-C_{\!f}(\tau )) \, {\rm d}\tau }{\int _{0}^{t} \overline {C}_{\!f}^N \, {\rm d}\tau } \times 100\,\%. \end{equation}
Since we are comparing the polymer case to the average steady-state Newtonian turbulence,
$\overline {C}_{\!f}^N$
is a constant, and the denominator in (B5) becomes
$t\overline {C}_{\!f}^N$
. The expression may then be rewritten as
\begin{equation} \textrm{DR}(t) = \frac {1}{t} \int _{0}^{t} \frac {\overline {C}_{\!f}^N - C_{\!f}(\tau )}{\overline {C}_{\!f}^N} \, {\rm d}\tau \times 100\,\%, \end{equation}
which is identical to the continuously averaged drag reduction presented in (3.7).
Open access
