Two-way coupling of FENE dumbbells with a turbulent shear flow
Thomas Peters∗
Department of Physics, Philipps-Universität Marburg, D-35032 Marburg, Germany
Jörg Schumacher†
Department of Mechanical Engineering, Technische Universität Ilmenau, D-98684 Ilmenau, Germany
(Dated: August 3, 2021)
We present numerical studies for finitely extensible nonlinear elastic (FENE) dumbbells which are
dispersed in a turbulent plane shear flow at moderate Reynolds number. The polymer ensemble is
described on the mesoscopic level by a set of stochastic ordinary differential equations with Brow-
nian noise. The dynamics of the Newtonian solvent is determined by the Navier-Stokes equations.
Momentum transfer of the dumbbells with the solvent is implemented by an additional volume
forcing term in the Navier-Stokes equations, such that both components of the resulting viscoelas-
tic fluid are connected by a two-way coupling. The dynamics of the dumbbells is given then by
Newton’s second law of motion including small inertia effects. We investigate the dynamics of the
flow for different degrees of dumbbell elasticity and inertia, as given by Weissenberg and Stokes
numbers, respectively. For the parameters accessible in our study, the magnitude of the feedback
of the polymers on the macroscopic properties of turbulence remains small as quantified by the
global energy budget and the Reynolds stresses. A reduction of the turbulent drag by up to 20% is
observed for the larger particle inertia. The angular statistics of the dumbbells shows an increasing
alignment with the mean flow direction for both, increasing elasticity and inertia. This goes in line
with a growing asymmetry of the probability density function of the transverse derivative of the
streamwise turbulent velocity component. We find that dumbbells get stretched preferentially in
regions where vortex stretching or bi-axial strain dominate the local dynamics and topology of the
velocity gradient tensor.
PACS numbers: 47.27.ek, 83.10.Mj, 83.80.Rs
I. INTRODUCTION
When a few parts per million in weight of long-chained polymers are added to a turbulent fluid its properties
change drastically and a significant reduction of turbulent drag is observed. [1] Although the phenomenon is known
from pipe flow experiments for almost 60 years,[2, 3] a complete understanding is still lacking. One reason for
this circumstance is that the physical processes in a turbulent and dilute polymer solution cover several orders of
magnitude in space and time; in other words, we are faced with a real multiscale problem. [4, 5] In case of fully
developed turbulence, the integral scale L, which measures the extension of largest vortex structures in the flow,
exceeds the viscous Kolmogorov scale ηK, which stands for the extension of the smallest turbulent eddies, by a factor
of at least 1000. However, long-chained polymers barely exceed the viscous flow scale even in an almost stretched
state. Their equilibrium extension as given by the Flory radius R0 is usually by a factor of 100 smaller than ηK.[6]
In terms of time scales the situation differs slightly. The viscous Kolmogorov time τη can become smaller than the
slowest relaxation time τ of the macromolecules. Although macroscopic closures can rationalize some issues of drag
reduction [7], the challenging question remains of how the individual dynamics of numerous polymer chains, which is
present on sub-Kolmogorov and Kolmogorov scales, adds up to a macroscopic effect at scales r <∼ L as being observed
in several experiments. [8, 9, 10]
The description of dilute polymer solutions relies for most studies on one of the following two models: on one
side, macroscopic continuum models such as Oldroyd-B or FENE-P models [11, 12, 13, 14, 15] include the polymer
dynamics as an additional additive macroscopic stress field. Only the largest scales ℓ >∼ ηK of the viscoelastic fluid
are described in its full complexity. Numerical problems arise in connection with the pure hyperbolic character of
the equation of motion for the polymer stress field, such as the conservation of its positivity (see e.g. Ref. [16]
for a detailed discussion). In addition, the coarse graining to the macroscopic polymer stress can lead to deeper
conceptional difficulties, e.g., the failure of energy stability of viscoelastic flows, which is an important building block
for investigations of stability and upper bounds on the dissipation rate in Newtonian flows. [17] Further problems
∗ Present address: Institute for Theoretical Astrophysics, Ruprecht-Karls-Universität Heidelberg, D-69120 Heidelberg, Germany
† Corresponding author: joerg.schumacher@tu-ilmenau.de
http://arxiv.org/abs/0704.1255v1
arise for the macroscopic description of non-Newtonian fluids in the limits of very low and high frequencies, where
they should behave as Newtonian fluids and solids, respectively. [18, 19, 20]
On the other side, Brownian dynamics models [21, 22, 23, 24, 25] describe the polymer chain on a mesoscopic level
as overdamped coupled oscillators arranged in bead-spring chains. The models include complex conformations of the
macromolecules and screening effects due to the solvent such as hydrodynamic interaction.[26] The simplest of such
mesoscopic models for a polymer chain is a dumbbell where two beads are connected by a spring. The dynamics in
these models is on scales ℓ <∼ ηK. This means that the surrounding fluid is spatially smooth and either a steady [22],
a start-up shear flow [23], or a white-in-time random flow. [27] In a recent work by Davoudi and Schumacher[28],
numerical studies at the interface of both descriptions were conducted by combining Brownian dynamics simulations
(BDS) with direct numerical simulations of a turbulent Navier-Stokes shear flow. The simplest mesoscopic model
with a linear spring force - the Hookean dumbbell model - was taken there in order to study the stretching of the
dumbbell as a function of the outer shear rate and the elastic properties of the springs. However, a feedback of the
polymers on the shear flow was not included in their study.
In the following, we want to extend these investigations into two directions. Firstly, we will model the macro-
molecules more realistically as finitely extensible nonlinear elastic (FENE) dumbbells. Secondly, their feedback on
the shear flow is included via a two-way coupling. The effect of the FENE dumbbells on the statistical fluctuations of
the velocity and the velocity gradients will be studied. In addition, conformational properties of the dumbbells, such
as their extension and angular distribution with respect to the mean flow component, will be addressed. The polymer
feedback results in an additional forcing that has to be added to the right hand side of the Navier-Stokes equations
for the advecting Newtonian solvent similar to the case of two-phase flows with dispersed particles [29, 30, 31] or
bubbles.[32] We will keep the full dynamic equation of motion for the dumbbells, containing accelerations due to
elastic, friction and stochastic forces, and cannot neglect inertia. This step is necessary in order to describe the
momentum transfer of the dumbbells to the solvent as discussed in Ref.[33].
In contrast to the conventional BDS that neglect inertia effects from beginning, we will be left here with three
physical parameters: the Stokes number St for the particle inertia, the Weissenberg number Wi for the elastic
properties of the dumbbells, and the Reynolds number Re of the flow, respectively. The Reynolds number is defined
, (1)
with the characteristic (large-scale) velocity U , the characteristic length L (both are specified later in the text), and
the kinematic viscosity of the Newtonian solvent ν. The Weissenberg number Wi compares the characteristic dumbbell
relaxation time τ from a stretched to a coiled state with the characteristic time scale of the advecting flow, L/U , and
is given by
. (2)
The Stokes number St relates the particle response time to changes in the surrounding velocity, τst, with the charac-
teristic flow time scale. It follows to
. (3)
The physics of dispersed FENE dumbbells in a turbulent shear flow is thus described by three dimensionless numbers.
For a fixed Reynolds numbers Re, we can basically distinguish the following four limiting cases: (i) Wi ≫ 1, St ≫ 1;
(ii) Wi ≪ 1, St ≫ 1; (iii) Wi ≪ 1, St ≪ 1; (iv) Wi ≫ 1, St ≪ 1. Case (i) would stand for very heavy particles (or
dumbbells) which are stretched almost to their contour length. They will behave as dispersed rods. In case (ii), the
dumbbells would act as heavy spherical particles since they remain coiled in practical terms. The cases of interest for
dilute polymer solutions are (iii) and (iv), respectively. Inertia effects are then very small, [34] and the Weissenberg
number can vary from very small to large values implying an increasingly slower relaxation of the macromolecules
from a stretched non-equilibrium to a coiled equilibrium state in comparison to the characteristic flow variation time
scale. As we will discuss in the next section, the numerical treatment becomes challenging, on one hand due to the
finite extensibility, on the other hand due to the small Stokes numbers we are aiming at. The Stokes time τst sets
a small but finite time scale then, which can cause stiffness problems for an explicit integration algorithm. Despite
these efforts, our values for the Stokes number will still exceed the realistic magnitudes for polymer chains in solution
by orders of magnitude. Nevertheless, we think it is interesting and to some degree necessary to study the dumbbell
dynamics under these circumstances and to provide a systematic study of how a shear flow will be affected by the
presence of dispersed bead-spring chains with variable degree of inertia. This will shed some light on possible reasons
for drag reduction in our model.
The outline of the manuscript is as follows. In the next section the equations of motion, the two-way coupling
and the numerical scheme are presented. Afterwards, we discuss the results for the macroscopic energy balance as
well as for the Reynolds stresses. This is followed by studies of small-scale properties such as the statistics of the
extension and orientation of the dumbbells and of their impact on the fluctuations of velocity gradients. We conclude
with a discussion of our results and will give a brief outlook to extensions of the present work toward more realistic
parameter settings.
II. MODEL AND EQUATIONS
A. The Newtonian solvent
The Navier-Stokes equations that describe the dynamics of the three-dimensional incompressible Newtonian fluid
are solved by a pseudo-spectral method using a second-order predictor-corrector scheme for advancement in time.[28]
The equations of motion are
+ (u ·∇)u = −∇p+ ν∇2u+ f + fp , (4)
∇ · u = 0 , (5)
where u is the (total) velocity field, p the kinematic pressure field, f the volume forcing which sustains the turbulence,
and fp the feedback of the dumbbells (see section II C). The shear flow is modeled in a volume with free-slip boundary
conditions in the shear direction y and periodic boundaries in the streamwise and spanwise directions x and z. The
free-slip boundary conditions at y = 0, Ly are given by
uy = 0 ,
= 0 . (6)
Here, the total velocity field follows by a Reynolds (de)composition as a linear mean part with the constant shear
rate S and a turbulent fluctuating part
u = 〈u〉+ u′ = Syex + u′ . (7)
The notation 〈·〉 stands for the ensemble average, which will be a combination of volume and time averages for most
cases. The aspect ratio is Lx:Ly:Lz = 4π: 2: 2π. The characteristic length is the halfwidth of the slab, L = Ly/2.
Velocities are measured in units of the laminar flow profile U(y) = −
2 cos(πy/2)ex. We will take Ux(Ly/4) as the
characteristic velocity U (see also (1), (2), and (3)). The applied volume forcing sustains this laminar flow profile
and follows from (4) consequently to f(y) = −
2π2/(4ν) cos(πy/2)ex. Forcing amplitude and profile will remain
unchanged throughout this study. At sufficiently large Reynolds numbers this linearly stable laminar shear flow
becomes turbulent when a finite perturbation is applied.[35] The volume forcing f is then a permanent source of
kinetic energy injection into the shear flow which sustains turbulence in a statistically stationary state. Although
the steady forcing is of cosine shape, the resulting mean turbulent flow profile will be linear except for small layers
in the vicinity of both free-slip planes, where the boundary conditions have to be satisfied. Our mean profiles follow
to 〈ux(y)〉 ≃ S(y − 1) for y ∈ [0, 2] with S = 0.035 − 0.04 for Re = 800. This range of S-values remained nearly
unchanged for all parameter sets. In addition, 〈u′y〉 = 〈u′z〉 = 0. The shear flow can be considered therefore as being
nearly homogeneous.
The simulation program is run with two spectral resolutions. For Re = 400, a grid with 64 × 32 × 32 mesh
points was taken. For Re = 800, we took a grid with 128 × 32 × 64 points. The spectral resolution as given by
the product kmaxηK =
8πNx/(3Lx)ηK was 1.5 for the first case and 2.3 for the second. Here, ηK is the viscous
Kolmogorov scale and defined as ηK = ν
3/4/〈ε′〉1/4 with the mean turbulent energy dissipation rate 〈ε′〉, where
ε′(x, t) = (ν/2)(∂u′i/∂xj + ∂u
j/∂xi)
2 for i, j = x, y, z. Clearly, the spectral resolutions are not very large, but they
give us the opportunity to perform parametric studies in the three-dimensional space which is spanned by Re, Wi,
and St. Most of our following studies will be conducted for the better resolved case of Re = 800.
B. The FENE dumbbells
The smallest building block for the mesoscopic description of the polymer stretching can be accomplished by
considering dumbbells where two beads (that stand for several hundreds of monomers) are connected by a spring.
The entropic elastic force follows the Warner force law [11] and depends on the separation vector R(t) = x2(t)−x1(t)
that is spanned between both beads at positions x2(t) and x1(t), respectively. The force law is given by
Fel(R) =
1−R2/L20
, (8)
where L0 is the contour length of the dumbbells which cannot be exceeded. The spring constant is denoted byH . When
taking into account the elastic entropic force, hydrodynamic Stokes drag, and thermal noise, the second Newtonian
law for a FENE dumbbell written in relative coordinates R(t) and center-of-mass coordinates r(t) = (x1(t)+x2(t))/2
reads [27, 34]
ṙ = v , (9)
v̇ = −v + 1
(u1 + u2) +
ξr , (10)
Ṙ = V , (11)
V̇ = −V +∆u− 2HR
ζ (1−R2/L20)
ξR , (12)
where ∆u = u(x2, t) − u(x1, t) is the relative fluid velocity at the bead centers. The last terms in the velocity
equations, containing ξr and ξR, stand for vectors of thermal Gaussian noise with the properties
〈ξi(t)〉 = 0 , (13)
〈ξi(t)ξj(t′〉 = δijδ(t− t′) (14)
for i, j = x, y, z. The three components of each vectorial noise term are statistically independent stochastic processes.
Furthermore, the vectorial noise with respect to the center-of-mass velocity is statistically independent to that for the
relative velocity dynamics. The noise prevents the extension of a dumbbell to shrink below its equilibrium length
, (15)
with kB being the Boltzmann constant, T the temperature. Equation (15) follows from the equipartition theorem.
The contour length L0 = 10R0 is used throughout this study and R0 ≃ ηK. The relaxation time of the dumbbells is
given by [11]
, (16)
where
ζ = 6πρfνa (17)
is the Stokes drag coefficient of a spherical bead with radius a. The fluid mass density is ρf . Due to the current
resolution contraints the dumbbells will experience both the smooth and partly rough scales of the advecting flow.
Consequently, the velocity difference ∆u is kept in the equation and not approximated by the linearization ∆u ≈
(R · ∇)u as it is done in BDS where L0 ≪ ηK . For spatially smooth flows both expressions give the same results.
The equations (9) through (12) introduce the other two dimensionless parameters beside the Reynolds number Re,
the Weissenberg number Wi and the Stokes number St, respectively (see definitions (2) and (3)). The Stokes time τst
is the response time of an inertial particle which is required to speed up to the velocity of its local surrounding. A
zero Stokes time implies a behavior as a passive Lagrangian tracer. For beads, this time follows to τst = mb/ζ with
ζ as given above and consequently
τst =
. (18)
The density contrast ρp/ρf is to very good approximation unity[36], i.e. polymers are considered as neutrally buoyant.
In Ref. [28], we have compared the polymer relaxation time to the microscopic stretching time scale. This is given
by the inverse of the maximum Lyapunov exponent and is comparable to the microscopic time scale of the flow, the
Re = 400 Re = 800
Wi = 3 Wiη = 0.8 Wiη = 0.6
Wi = 20 Wiη = 5.1 Wiη = 4.3
Wi = 100 Wiη = 25.7 Wiη = 21.5
St = 5.0 × 10−4 Stη = 1.3× 10
−4 Stη = 1.1 × 10
St = 5.0 × 10−3 Stη = 1.3× 10
−3 Stη = 1.1 × 10
St = 5.0 × 10−2 Stη = 1.3× 10
−2 Stη = 1.1 × 10
St = 5.0 × 10−1 Stη = 1.3× 10
−1 Stη = 1.1 × 10
TABLE I: The Weissenberg and Stokes numbers rescaled by the Kolmogorov time τη of the flow. Wiη = τ/τη and Stη = τst/τη.
Note that τη is based on the pure Newtonian case. Only minor changes arise when polymers are added to the solvent.
Kolmogorov time τη =
ν/〈ε〉. Table 1 gives an overview of the values of St and Wi that have been used and of how
they translate into Stη and Wiη, respectively. We see that the Stokes numbers get as low as 10
−4 when measured in
viscous units, which is still orders of magnitude above the realistic estimates for dilute polymer solutions which are
about three to four order of magnitude below our minimal value.
In most cases, an ensemble of 6.3× 104 FENE dumbbells, i.e. 1.2× 105 beads, is advanced by a weak second-order
predictor-corrector scheme simultaneously with the flow equations.[21] The finite extensibility and the small Stokes
numbers require a semi-implicit time-stepping for some variables. In order to avoid a total length larger than L0,
we proceed in line with Ref. [21] and solve a cubic equation for R = |R| in the corrector step. Initially, the center
of mass of the dumbbells is seeded randomly in space with a uniform distribution and an initial extension of R0.
All Lagrangian interpolations were done with a trilinear scheme. Details on the numerical procedure are outlined in
appendix A.
In order to build a bridge to macroscopic simulations we provide an estimate for the contribution of the dumbbell
ensemble to the zero-shear viscosity. Following Ref. [21] it is defined as
ηp = ρpνp = npkBTτ , (19)
with the number density of dumbbells np. When applying (15) as well as definitions (16) and (17), and using ρf/ρp = 1
one gets
0νa (20)
with the solvent viscosity ν. The bead radius a is substituted by the Stokes time τst. Recalling the definitions for the
Kolmogorov length ηK = ν
3/4/〈ε′〉1/4 and for the Kolmogorov time τη =
ν/〈ε′〉, one ends with the relative viscosity
Stη . (21)
For the present simulations, one dumbbell is seeded per grid cell and therefore np ≈ 1/η3K. Additionally, R0 ≃ ηK.
Following table 1 for the runs at Re = 800, one gets ratios of s between between 0.1 for the smallest Stokes number
going up to 3 for the largest one. The latter value is rather large for polymer solutions. Values below unity are usually
taken, such as in DNS with the Oldroyd-B model.[14] Equation (21) is in this spirit consistent with the discussion
in the introductory part. Only the lower Stokes numbers result to values of s as taken for macroscopic DNS for
viscoelastic shear flows.
C. Two-way coupling
The back-reaction of the dumbbells on the fluid consists of contributions from the Stokes friction and the stochastic
noise term. In accordance with Newton’s third law, the force contribution from each of the two beads at positions xi
(i = 1, 2) follows to
Fi = −F (st)i − F
i = ζ(ẋi − u(xi))−
2kBTζ ξi . (22)
The force density generated by all FENE dumbbells results to
ρffp =
i δ(x− x
i ) , (23)
where Np is the number of dumbbells. The volume integral of (23) gives a force since the delta function carries the
dimension of an inverse volume due to
δ(x − x(j)i ) d
3x = 1. Consequently, the dimensionless form of the forcing
reads
i − u(x
i ))−
δ̃(x− x(j)i ) , (24)
where the bead volume follows to Vb = 4πa
3/3 = (4π/3)(9ντst/2)
3/2. The notation δ̃ is for the dimensionless delta
function. We have used again ρf/ρp ≈ 1. The force density has to be evaluated at space points that are between the
mesh vertices. Again the trilinear interpolation has to be used to evaluate the contributions of the point force to the
eight next neighboring mesh vertices.
III. LARGE-SCALE PROPERTIES
A. Energy balance
The first analysis step is the study of the effects of the two-way coupling on the macroscopic properties of turbulence.
Given the boundary conditions for our problem, eq. (4) results in the following balance for the total kinetic energy
E(t) = 1
|u|2 d3x with V = LxLyLz,
= −ν〈(∂ui/∂xj)2〉V + 〈u · f〉V + 〈u · fp〉V ,
= −ε(t) + εin(t)− εp(t) (25)
where 〈·〉V = 1V
· d3x is the short notation for the volume average. In case of statistical stationarity, one gets
d〈E〉t/dt = 0 and thus
〈εin〉 = 〈ε〉+ 〈εp〉 . (26)
Figure 1 shows the three mean rates as a function of the Stokes number for two Weissenberg numbers Wi = 20, 100.
The mean energy dissipation rate 〈ε〉 and the mean energy injection rate 〈εin〉 are of the same order of magnitude for
all cases. They remain nearly unchanged with respect to Weissenberg number, which indicates that the effect of the
dumbbell ensemble on the macroscopic flow properties is small. Nevertheless, one observes a slight increase of the
mean energy injection rate 〈εin〉 with respect to St going in line with a decrease of 〈ε〉 (see upper and mid panel of
Fig. 1). Recall that the energy injection rate will be maximal for the laminar case, i.e. for u ‖ f . The trend of the
data indicates that the streamwise flow component relaminarizes slightly with growing inertia. The lower panel of the
same figure shows the findings for the dissipation due to polymer stretching 〈εp〉. As an additional energy dissipation
mechanism, it consumes injected energy which goes into the elastic energy budget of the dumbbell ensemble. The
rate 〈εp〉 grows in magnitude with respect to both parameters, the Stokes and Weissenberg number. For Wi = 3, the
dumbbells are not significantly extended and no clear trend of 〈εp〉 with St could be observed. The dissipation rate
〈εp〉 is significantly smaller in comparison to the runs with larger Wi.
In order to estimate the maximum feedback of the dumbbells on the flow, we performed an “academic experiment”
for our system by tethering one of the two beads of a dumbbell at a fixed position. The dumbbells get then stretched
more efficiently and undergo strong conformational fluctuations. Figure 2 illustrates their dramatic effect on the total
kinetic energy. We compare the freely draining case with the tethered one and observe a significant decrease of the
kinetic energy. An inspection of the flow structures indicates that the turbulent fluctuations are supressed almost
completely. The flow becomes basically laminar. The magnitude of the feedback for freely draining dumbbells will
always remain significantly below this artifical limit with tethered dumbbells.
B. Reynolds stresses
Figure 3 shows the four non-vanishing components of the Reynolds stress tensor 〈u′iu′j〉/(2k) where k = 〈(u′i)2〉/2
is the turbulent kinetic energy (TKE). The moments are averages over the whole simulation volume for a sequence of
about 100 statistically independent snapshots of the time evolution of the shear flow. The results can be summarized
to the following trends. For the two smallest Stokes numbers, no dependence on the Weissenberg number is observed.
For St = 0.05 and 0.5, the mean streamwise fluctuations are enhanced while the remaining components of the Reynolds
stress tensor decrease as a function of Wi. This finding is in agreement with observations in a Kolmogorov flow by
Boffetta et al. [37]
Similar to the friction factor for a turbulent pipe [38], we can define a friction factor for the present flow where
the applied pressure gradient term has to be substituted by an amplitude of the static volume forcing profile f that
sustains the laminar cosine flow profile. Consequently,
〈ux(y = Ly)〉2
. (27)
Since f(y) = −
2π2/(4ν) cos(πy/2)ex, we take F = fx(y = Ly) =
2π2/(4ν). A similar definition was suggested for
a Kolmogorov flow which is also driven by a volume forcing.[37] Drag reduction by dispersed dumbbells would go in
line with a decrease of the dimensionless measure cf below the Newtonian value c
f . For the smallest Stokes number,
the ratio goes to about unity. The slight overshoot is attributed to the strong variations of the streamwise velocity
at the free-slip planes. Figure 4 indicates a reduction by 20%− 25% at St = 0.05, 0.5 and for the larger Weissenberg
numbers. The series with Wi = 3 gave cf ≃ cNf .
An important structural ingredient of shear flows are the asymmetric fluctuations of the three diagonal elements
of the Reynolds stress tensor. The streamwise fluctuations 〈(u′x)2〉 are spatially arranged in streamwise streaks which
interact with streamwise vortices in a so-called regeneration cycle of coherent structures. This cycle is sustained
by the non-normal amplification mechanism.[39, 40] The impact of long-chained polymers on the extension of the
streamwise streaks has been demonstrated in experiments [10] and numerical simulations.[41, 42] While streamwise
fluctuations were found to increase, the fluctuations in shear and spanwise directions decreased. This is in line with
our observations as discussed above. In Fig. 5, we show isolevels of the streamwise turbulent fluctuations for opposite
sign at Wi = 3, 20, 100. Although not very pronounced, a slight increase in the connectivity and extension of the
streamwise streaks can be observed with increasing Weissenberg number.
As we can see, the statistics of macroscopic turbulent properties is affected only slightly by the dispersed FENE-
dumbbells. Their impact increases with Weissenberg number as well as with Stokes number. In order to rule out
that particle inertia dominates the discussed trends of our studies, we considered the case of dispersed beads in the
same flow at the same Stokes numbers. This is achieved by switching off the elastic spring force, i.e. Fel = 0. The
Stokes friction force remained as the only force. The quantity fp models then the feedback of the particles on the
flow. We added the statistical means of injection and dissipation rates as a function of the Stokes number for this
case to Fig. 1. While the mean injection and mean dissipation rates are of the same magnitude, the dissipation due
to particle feedback is orders of magnitude smaller in comparison to the polymer feedback, except for the largest St.
In addition, we found no clear trends for the Reynolds stress components as a function of St.
IV. SMALL-SCALE PROPERTIES
A. Extensional and angular statistics of dumbbells
The finite extensibility of the dumbbells will affect the shape of the probability density function (PDF) of R, which
is supported on scales smaller than L0 only. Figure 6 reports our findings for p(R) for different Weissenberg and
Stokes numbers. For the lowest Weissenberg number, Wi = 3, the majority of the dumbbells remains at the extension
of about the Kolmogorov length ηK. This picture changes for larger values of Wi. At Wi = 100, the majority of the
ensemble is stretched to almost L0, which manifests in the sharp maximum at R <∼ L0. Qualitatively, the change of
the shapes of the PDFs with increasing Wi agrees well with experimental findings [43] and analytical studies [27, 44]
for the coil-stretch transition in random flows. The trends with the Stokes number remain small in all cases. However,
the data show that growing particle inertia suppresses the stretching to very extended molecules since the response
time of the molecules to the variation of the structures increases (see e.g. mid panel of Fig. 6).
As we have seen in the last section, the fluctuations of the turbulent velocity field in the shear flow vary strongly from
one space direction to another (see e.g. Fig. 3). The major contribution is contained in the streamwise component
〈(u′x)2〉 parallel to the direction of the mean turbulent flow. This suggests an investigation of the angular statistics of
the polymers since their stretching can be expected to become anisotropic as well. The following dumbbell coordinate
system will be used therefore throughout this text: Rx = R cosϕ cos θ, Ry = R sinϕ cos θ, and Rz = R sin θ, where
R is the distance between both beads. The notation differs from conventional spherical coordinates, but has the
advantage of giving perfect alignment with the outer mean flow direction for ϕ = θ = 0. ϕ is the azimuthal angle and
θ the polar angle. While the azimuthal angle always remains in the shear plane that is spanned by the streamwise
and shear directions, the polar angle θ 6= 0 indicates a dumbbell orientation out of this plane.
Davoudi and Schumacher [28] discussed the statistics of both angles as a function of the Weissenberg number for
passively advected Hookean dumbbells. The PDF of the polar angle was found to remain symmetric and to be less
sensitive with respect to variations of Wi. Our focus will be therefore on the statistics of the azimuthal angle ϕ
which can take values between −π/2 and π/2. The asymmetry between both quadrants is quantified by the following
measure for the PDF p(ϕ):
A(ϕ) = p(ϕ)− p(−ϕ) , (28)
with ϕ ∈ [0, π/2]. The measure A(ϕ) is plotted for two Weissenberg numbers in Fig. 7. A pronounced maximum
of A(ϕ) implies that the dumbbells are preferentially slightly tilted in the direction of shear, away from the mean
flow direction (see an illustration in Fig. 8). We find that with increasing Weissenberg number the asymmetry of the
angular distribution grows in magnitude. The same trend holds when the Stokes number grows at fixed Weissenberg
number. In each case, the graph of A(ϕ) shows an increasingly sharper maximum, which is shifted towards smaller
ϕ. Fluctuations of the dumbbells in the vicinity of ϕ = 0 are enhanced while the tails for very large ϕ are depleted.
Growing inertia amplifies this trend. Once the dumbbells are aligned along the mean flow they remain in this
orientation for longer periods of their evolution.
B. Velocity gradient statistics
Since the polymer dynamics takes place at the smallest scales of the turbulent flow, we study the impact of the
dumbbells on the small-scale statistical properties of the flow in the following. Recent experimental and numerical
studies in simple Newtonian shear flows indicate that in particular the statistics of the transverse derivative of the
streamwise turbulent velocity component ∂u′x/∂y is a sensitive measure for detecting deviations from local isotropy
in homogeneous or nearly homogeneous shear flows.[45, 46, 47] In a shear flow with a mean shear rate S > 0, one
expects a positive value for derivative skewness and other higher odd order moments which are defined as
M2n+1(∂u
x/∂y) =
〈(∂u′x/∂y)2n+1〉
〈(∂u′x/∂y)2〉n+1/2
. (29)
The derivative moments would be exactly zero in a perfectly isotropic flow. Their non-zero magnitudes indicate that
velocity gradient fluctuations of the streamwise component along the direction of the outer shear gradient are more
probable than the ones in the opposite direction. It can be expected that the asymmetry in the angular distribution,
which we discussed above, will have an impact on the statistics of exactly these gradient fluctuations. Figure 9
reports our findings for the PDF of the transverse derivative, which has been normalized by its root mean square
value for all cases. We observe in both figures a depletion of the left hand tail, which stands exactly for the velocity
gradient fluctuations opposite to the direction of the mean shear. The results suggest that the preferential orientation
fluctuations of the dumbbells at azimuthal angles ϕ > 0 go in line with a depletion of the negative tail of the PDF
of the transverse derivative. As sketched in Fig. 8, negative transverse gradients would be amplified by prefential
orientations with ϕ < 0 which correspond to the dumbbell colored in gray. The findings are consistent with our
observations on the ϕ-statistics. They can also be rationalized (but not explained) when considering the equation for
the Brownian dynamics of the FENE dumbbell [21]
= R · ∇u−
2τ(1 −R2/L20)
ξR . (30)
In the plane shear flow geometry the component Rx along the mean flow direction is of particular interest. Since we
are interested in stretched dumbbells with Rx > R0 and in Wi > 1 we neglect contributions from the spring force and
the noise for a moment. With the Reynolds decomposition (7) we get
Rx + ... , (31)
Rx + ... (32)
The important term is the first term on the r.h.s. of (31). The other three contributions will behave as noise terms.
Fluctuating gradients ∂u′x/∂y along Sey lead to a more rapid growth of Rx (for an angle ϕ > 0) and a prefered
alignment with the mean flow. This causes a more rapid decrease of Ry and consequently of Rx via (31). The dumbell
can be kicked afterwards again to larger ϕ values and transfers momentum to the flow which corresponds exactly to a
local patch of ∂u′x/∂y > 0 (see also Fig. (8)). Then Ry grows and this whole cycle starts anew. Small scale gradients
with the opposite sign diminish the total shear in the surrounding of the dumbbell and cause a less efficient stretching
and cycle. Clearly, this picture omits some important features such as the tumbling of the dumbbells.
The depletion of gradient fluctuations goes in line with experimental observations by Liberzon et al. [48, 49] The
authors found e.g. that the enstrophy production became anisotropic when polymers are added to the fluid. This
quantity is directly related to transverse gradient components discussed here.
C. Invariants of the velocity gradient tensor and dumbbell extension
The efficient stretching of the dumbbells is connected to particular local flow topologies. They are related to the
three eigenvalues λi of the velocity gradient tensor or the corresponding three velocity gradient tensor invariants,
which are denoted as I1, I2, and I3. The eigenvalues of the velocity gradient tensor ∂u
i/∂xj result as zeros of the
following third-order characteristic polynomial[50]
λ3 − I1λ2 + I2λ− I3 = 0 . (33)
For an incompressible flow [53],
I1 = λ1 + λ2 + λ3 = Tr
= 0 ,
I2 = λ1λ2 + λ2λ3 + λ3λ1 = −
I3 = λ1λ2λ3 = det
. (34)
The remaining coefficients of (33) are therefore I2 and I3, which span the I3 − I2 parameter plane. The scatter plots
for turbulent flows result in a typical skewed teardrop shape. With our definitions given above the following crude
classification scheme can be given. For I2 > 0, I3 > 0 vortex stretching is present corresponding to λ1 = a, λ2,3 =
−a ± ib (first quadrant); for I2 > 0, I3 < 0 vortex compression is present corresponding to λ1 = −a, λ2,3 = a ± ib
(second quadrant). The cases I3 < 0 are associated with bi-axial strain for I2 < 0 (third quadrant) corresponding to
λ1 = a, λ2 = b, λ3 = −(a + b) and with uniaxial strain at I2 > 0 (fourth quadrant) corresponding to λ1 = a, λ2 =
−b, λ3 = −(a−b). Constants a and b are larger than zero in all cases. Figure 10 relates the extension of the dumbbells
to the corresponding local velocity gradients in the I3−I2 plane (and consequently to the existing local flow topology).
The invariants of the velocity gradient were evaluated in the center of mass of each dumbbell. The typical teardrop
shape for the turbulence data in the parameter plane is detected.
Our findings can be summarized as follows. Strongly stretched dumbbells go in line with the largest excursions of
the gradients in the I3 − I2 plane. The longest dumbbells are found preferentially in regions where vortex stretching
or bi-axial strain dominate the local flow topology. The preferential stretching by bi-axial strain was discussed already
for the passive advection of FENE dumbbells in a minimal flow unit.[25] It corresponds to the scenario that different
parts of the dumbbell get pulled by counterstreaming streamwise streaks. The preferential extension close to vortex
stretching means that the polymers are pulled around streamwise vortices. This point was outlined in Ref. [42] on
the basis of an analysis of the energetics of viscoelastic turbulence. Here, we find both in a common description
based on the analysis of the full velocity gradient tensor, i.e. the symmetric strain tensor plus the anti-symmetric
vorticity tensor. We do also observe that the area of the teardrop shape shrinks with increasing Stokes number. This
indicates that the small-scale velocity gradients are supressed in magnitude, which goes in line with more limited
excursions across the I3 − I2 plane and a relaminarization of the turbulence. Again, this goes in line with very recent
experimental observations by Liberzon et al.[49]
V. SUMMARY AND DISCUSSION
The presented numerical studies aimed at connecting a macroscopic description for the Newtonian turbulent shear
flow to the mesoscopic description of an ensemble of FENE dumbbells which are advected in such flow. The momentum
transfer of the dumbbells with the fluid is implemented by an additional volume forcing in the Navier-Stokes equations.
In numerical terms, pseudospectral simulations for the solvent are coupled to a system of stochastic nonlinear ordinary
equations in order to model a viscoelastic fluid.
For the accessible parameters we found slight modifications of the macroscopic flow structures and mean statistical
properties only. This was demonstrated for the global energy balance and the mean components of the Reynolds
stress tensor. We conclude that dumbbell inertia effects are present, but remain subleading in comparison to the
elastic properties. For the present viscoelastic flow a drag reduction of up to 20% is achieved. The microscopic
properties of turbulence were found to be more sensitive with respect to the Weissenberg number. The statistics of
the azimuthal angle ϕ is consistent with former findings for elastic Hookean dumbbells. [28] A growing number of
dumbbells becomes increasingly aligned with the mean flow direction. The feedback of the FENE dumbbells on the
small-scale properties of turbulence is demonstrated for two gradient measures, the PDF of the transverse derivative of
the turbulent streamwise velocity component ∂u′x/∂y and the diminished scattering of the velocity gradient invariants
ampiltudes in the I3 − I2 plane with increasing Wi. The asymmetry of the PDF p(∂u′x/∂y) is found to increase with
increasing Wi. Furthermore, we determined that strongly stretched dumbbells can be found close to vortex stretching
or biaxial strain topologies of the advecting shear flow.
The present study should be considered as a first step for such class of hybrid models. One difference to the
situation in a dilute polymer solution is the relatively large Stokes number that had to be taken. Our dispersed
dumbbells behave in parts like deformable particles rather than polymer chains. Frequently, heavier quasi-particles
are used for the study of turbulence in particle-ladden flows.[31] Extensions of our investigations will have to go into
two directions. Firstly, it is desirable that larger spectral resolutions, like the ones in Ref. [28], are achieved. This
will require a fully parallel implementation of the current numerical scheme. Larger computational grids and higher
Reynolds numbers will give us the opportunity to decrease the ratio R0/ηK and to increase L0/R0 to more realistic
values. Secondly, eq. (21) implies the efforts that have to be taken in order to approach the situation in a polymer
solution. Decreasing values of R0 and St have to be compensated by np, e.g., a reduction of both – R0 and Stη –
by an order of magnitude requires an increase of the concentration (or number density) by a power of 5/2. Once
such operating point is reached, the time scale argument which is thought to be important for the drag reduction
effect, can also be studied.[1] Finally, a recent work by Vincenzi and co-workers [51] provides an interesting ansatz for
modelling the polymer dynamics. The authors studied a conformation-dependent Stokes drag coefficient that caused
a significant dynamical slow-down of the coil-stretch transition in steady elongational and random flows. The test of
these ideas in turbulent shear flows is still to be done.
Acknowledgments
This work was supported by the Deutsche Forschungsgemeinschaft (DFG) and the Deutscher Akademischer Aus-
tauschdienst (DAAD) within the German-French PROCOPE program. We thank for computing ressources on the
JUMP supercomputer at the John von Neumann Institute for Computing, Jülich (Germany). Further computations
have been conducted at the MARC cluster (Marburg) and the MaPaCC cluster (Ilmenau). Fruitful discussions with
F. de Lillo, B. Eckhardt, and D. Vincenzi are acknowledged.
APPENDIX A: SEMI-IMPLICIT INTEGRATION SCHEME FOR DUMBBELLS
The FENE dumbbells consist of two beads at positions x1(t) and x2(t) which are connected by a nonlinear elastic
spring. The velocities of the advecting flow at both beads are denoted by u1 and u2, respectively. Note that these
velocities coincide with ẋ1 and ẋ2, respectively, for St = 0 only. Since the beads are usually found between mesh
vertices, the values for u1 and u2 have to be determined by trilinear interpolation from the known velocity vectors
at the neighboring grid sites. The dynamical equations for the dumbbells are set up in relative and center-of-mass
coordinates. The relative coordinate (or separation) vector of the dumbbell is given by
R(t) = x2(t)− x1(t) . (A1)
The center-of-mass coordinate vector is given by
r(t) =
(x1(t) + x2(t)) . (A2)
The velocities which are assigned with the relative and center-of mass coordinates are denoted as V and v, respectively.
The Newtonian equations for the dynamics of the FENE dumbbells in dimensionless form, which follow then from
(9)-(12) with the definitions (3) and (2), are given by
= ṽ , (A3)
−ṽ + 1
(ũ1 + ũ2) +
, (A4)
= Ṽ , (A5)
−Ṽ + (ũ2 − ũ1)−
1− R̃2L2/L20
 . (A6)
For the following, we omit the tilde symbol for the dimensionless quantities. The predictor values of the center-of-mass
vector r and the distance vector R are calculated by an explicit Euler step whereas the corresponding velocities are
treated by an implicit Euler step, giving
r∗ = rl +∆tvl , (A7)
St + ∆t
St vl +
(ul1 + u
2)∆t+
, (A8)
R∗ = Rl +∆tV l , (A9)
V ∗ =
St + ∆t
StV l + (ul2 − ul1)∆t−
2Wi (1− (Rl)2L2/L20)
. (A10)
The corrector step for the center-of-mass and distance vectors is given as
rl+1 = rl +
(v∗ + vl)∆t (A11)
l+1 =
St + ∆t
Stv∗ +
(u∗1 + u
2)∆t+ Stv
(ul1 + u
(A12)
l+1 = Rl +
(V l + V l+1)∆t (A13)
V l+1 =
St + ∆t
StV ∗ + (u∗2 − u∗1)∆t+ StV l + (ul2 − ul1)∆t−
2Wi (1− (Rl+1)2L2/L20)
∆t− R
2Wi (1− (Rl)2L2/L20)
∆W l] . (A14)
Note that the corrector step for the distance vector is semi-implicit in the velocity in order to avoid stiffness of the
equation system at small Stokes numbers. The corrector step for the distance velocity V has to be semi-implicit in
the separation vector R due to the finite extensibility of the dumbbells.[21] When inserting (A14) into (A13) one gets
(∆t)2
8Wi (St + ∆t) (1− (Rl+1)2L2/L20)
Rl+1 = A , (A15)
where the abbrevation A contains terms only which are known. By taking the norm of (A15) one ends up with a
cubic polynomial for Rl+1. The formula for the “casus irreducibilis” of three real solutions of the polynomial goes
back to F. Viète [52] and yields directly the unique solution for R = |R| between 0 and L0. From (A15) follows now
Rl+1 = Rl+1
. (A16)
This value is inserted into (A14) which completes the corrector step.
[1] J. L. Lumley, “Drag reduction by additives,” Annu. Rev. Fluid Mech. 1, 367 (1969).
[2] B. A. Toms, “Observations on the flow of linear polymer solutions through straight tubes at large Reynolds numbers,” in
Proceedings of the International Congress on Rheology (Holland 1948), Vol. 2, pp. 135-141 North-Holland, Amsterdam
(1949).
[3] P. S. Virk, “Drag reduction fundamentals,” AIChE J. 21, 625 (1975).
[4] G. H. McKinley and T. Sridhar, “Filament-stretching rheometry of complex fluids,” Annu. Rev. Fluid Mech. 34, 375
(2002).
[5] R. G. Larson, “The rheology of dilute solutions of polymers: Progress and problems,” J. Rheol. 49, 1 (2005).
[6] M. Doi, Introduction to polymer physics, Cambridge University Press, Cambridge, 1996.
[7] R. Benzi, E. De Angelis, V. S. L‘vov, I.Procaccia, and V. Tiberkevich, “Maximum drag reduction asymptotes and the
cross-over to the Newtonian plug,” J. Fluid Mech. 551, 185 (2006).
[8] M. D. Warholic, H. Massah, and T. J. Hanratty, “Influence of drag-reducing polymers on turbulence: effects of Reynolds
number, concentration and mixing,” Exp. Fluids 27, 461 (1999).
[9] K. R. Sreenivasan and C. M. White, “The onset of drag reduction by dilute polymer additives and the maximum drag
reduction asymptote,” J. Fluid Mech. 409, 149 (2000).
[10] C. M. White, V. S. R. Somandepalli, and M. G. Mungal, “The turbulence structure of drag reduced boundary layer flow,”
Exp. Fluids 36, 62 (2004).
[11] R. B. Bird, R. C. Amstrong, and O. Hassager, Dynamics of polymeric liquids, John Wiley & Sons, New York, 1987.
[12] R. Sureshkumar, A. N. Beris, and R. A. Handler, “Direct numerical simulation of the turbulent channel flow of a polymer
solution,” Phys. Fluids 9, 743 (1997).
[13] P. Ilg, E. De Angelis, I. V. Karlin, C. M. Casciola, and S. Succi, “Polymer dynamics in wall turbulent flow,” Europhys.
Lett. 58, 616 (2002).
[14] B. Eckhardt, J. Kronjäger, and J. Schumacher, “Stretching of polymers in a turbulent enviroment,” Comp. Physics Comm.
147, 538 (2002).
[15] C. D. Dimitropolous, Y. Dubief, E. S. G. Shaqfeh, P. Moin, and S. K. Lele, “Direct numerical simulation of polymer-induced
drag reduction in a turbulent shear flow,” Phys. Fluids 17, 011705 (2005).
[16] T. Vaithianathan and L. R. Collins, “Numerical approach to simulating turbulent flow of a viscoelastic polymer solution,”
J. Comp. Phys. 187, 1 (2003).
[17] C. R. Doering, B. Eckhardt, and J. Schumacher, “Failure of energy stability in Oldroyd-B fluids at arbitrarily low Reynolds
numbers,” J. Non-Newtonian Fluid Mech. 135, 92 (2006).
[18] H. Temmen, H. Pleiner, M. Liu, and H. R. Brand, “Convective nonlinearity in non-Newtonian fluids,” Phys. Rev. Lett.
84, 3228 (2000).
[19] A. N. Beris, M. D. Graham, I. Karlin, and H.-C. Öttinger, “Comment on Convective nonlinearity in non-Newtonian fluids,”
Phys. Rev. Lett. 86, 744 (2001).
[20] H. Temmen, H. Pleiner, M. Liu, and H. R. Brand, “Temmen et al. reply,” Phys. Rev. Lett. 86, 745 (2001).
[21] H.-C. Öttinger, Stochastic processes in polymeric fluids, Springer Verlag, Berlin, 1996.
[22] A. Celani, A. Puliafito, and K. Turitsyn, “Polymers in linear shear flow: a numerical study,” Europhys. Lett. 70, 464
(2005).
[23] J. S. Hur, E. S. G. Shaqfeh, H. P. Babcock, D. E. Smith, and S. Chu, “Dynamics of dilute and semidilute DNA solutions
in the start-up of shear flow,” J. Rheol. 45, 421 (2001).
[24] P. A. Stone and M. D. Graham, “Polymer dynamics in a model of the turbulent buffer layer,” Phys. Fluids 15, 1247 (2003).
[25] V. E. Terrapon, Y. Dubief, P. Moin, E. S. G. Shaqfeh, and S. K. Lele, “Simulated polymer stretch in a turbulent flow
using Brownian dynamics,” J. Fluid Mech. 504 61 (2004).
[26] B. H. Zimm, “Dynamics of polymers in dilute solution - viscoelasticity, birefringence and dielectric loss,” J. Chem. Phys.
24, 269 (1956).
[27] A. Celani, S. Musacchio, and D. Vincenzi, “Polymer transport in random flow,” J. Stat. Phys. 118, 531 (2005).
[28] J. Davoudi and J. Schumacher,“Stretching of polymers around the Kolmogorov scale in a turbulent shear flow,” Phys.
Fluids 18, 025103 (2006).
[29] K. D. Squires and J. K. Eaton, “Particle response and turbulence modification in isotropic turbulence,” Phys. Fluids A 2,
1191 (1990).
[30] S. Elghobashi and G. C. Truesdell, “On the two-way interaction between homogeneous turbulence and dispersed solid
particles. I: Turbulence modification,” Phys. Fluids A 5, 1790 (1993).
[31] T. Bosse, L. Kleiser, and E. Meiburg, “Small particles in homogeneous turbulence: Settling velocity enhancement by
two-way coupling,” Phys. Fluids 18, 027102 (2006).
[32] I. M. Mazzitelli, D. Lohse, and F. Toschi, “On the relevance of the lift force in bubbly turbulence,” J. Fluid Mech. 488,
283 (2003).
[33] P. Ahlrichs and B. Dünweg, “Simulation of a single polymer chain in solution by combining lattice Boltzmann and molecular
dynamics,” J. Chem. Phys. 111, 8225 (1999).
[34] J. D. Schieber and H.-C. Öttinger, “The effect of bead inertia on the Rouse model,” J. Chem. Phys. 89, 6972 (1988).
[35] J. Schumacher and B. Eckhardt, “Evolution of turbulent spots in a plane shear flow,” Phys. Rev. E 63, 046307 (2001).
[36] P. J. Flory, Principles of Polymer Chemistry, Cornell University Press, Ithaca, 1953.
[37] G. Boffetta, A. Celani, and A. Mazzino, “Drag reduction in the turbulent Kolmogorov flow”, Phys. Rev. E 71, 036307
(2005).
[38] H. T. Schlichting, Boundary Layer Theory, McGraw Hill, New York, 1979.
[39] F. Waleffe, “On a self-sustaining process in shear flows,” Phys. Fluids 9, 883 (1997).
[40] S. Grossmann, “The onset of shear flow turbulence,” Rev. Mod. Phys. 72 , 603 (2000).
[41] P. A. Stone, A. Roy, R. G. Larson, F. Waleffe, and M. D. Graham, “Polymer drag reduction in exact coherent structures
of plane shear flow,” Phys. Fluids 16, 3470 (2004).
[42] Y. Dubief, C. M. White, V. Terrapon, E. S. G. Shaqfeh, P. Moin, and S. K. Lele, “On the coherent drag-reducing and
turbulence-enhancing behaviour of polymers in wall-flows,” J. Fluid Mech. 514, 271 (2004).
[43] S. Gerashchenko, C. Chevallard, and V. Steinberg, “Single polymer dynamics: coil-stretch transition in a random flow,”
Europhys. Lett. 71 221 (2005).
[44] M. Chertkov, I. Kolokolov, V. Lebedev, and K. Turitsyn, “Polymer statistics in a random flow with mean shear,” J. Fluid
Mech. 531, 251 (2005).
[45] Z. Warhaft, “Turbulence in nature and laboratory,” Proc. Nat. Acad. Sci. 99, 2481 (2002).
[46] J. Schumacher, K. R. Sreenivasan, and P. K. Yeung, “Derivative moments in turbulent shear flows,” Phys. Fluids 15, 84
(2003).
[47] L. Biferale and I. Procaccia, “Anisotropy in turbulent flows and in turbulent transport,” Phys. Rep. 414, 43 (2005).
[48] A. Liberzon, M. Guala, B. Lüthi, W. Kinzelbach, and A. Tsinober, “Turbulence in dilute polymer solutions,” Phys. Fluids
17, 031707 (2005).
[49] A. Liberzon, M. Guala, W. Kinzelbach, and A. Tsinober, “On the kinetic energy production and dissipation in dilute
polymer solutions,” Phys. Fluids 18, 125101 (2006).
[50] P. A. Davidson, Turbulence, Oxford University Press, Oxford, 2004.
[51] A. Celani, A. Puliafito, and D. Vincenzi, “Dynamical slowdown of polymers in laminar and random flows,” Phys. Rev.
Lett. 97, 118301 (2006).
[52] T. Needham, Visual Complex Analysis, Oxford University Press, Oxford, 1997.
[53] We will use the definitions as given in Ref. [50] (p. 264) with I1 = P , I2 = Q, and I3 = R
PSfrag replacements
Wi = 20
Wi = 20
Wi = 100
Wi = 100
Fel = 0
Fel = 0
FIG. 1: Mean dissipation and injection rates as a function of the Stokes and Weissenberg numbers. Upper picture: mean energy
injection rate due to shear flow forcing 〈εin〉. Mid panel: mean energy dissipation rate 〈ε〉. Lower panel: mean dissipation rate
which arises from the coupling to the dumbbell ensemble 〈εp〉. The Reynolds number is Re = 800. The case with Fel = 0 is for
the case without spring force and stands for a shear flow with dispersed inertial particles.
PSfrag replacements
Wi = 20, St = 0.05
tethered
0 5 10 15 20
FIG. 2: Comparison of the kinetic energy for two cases at Wi = 20 and St = 0.05: for the tethered case one bead of each
dumbbell is fixed at a grid site while the second bead can fluctuate. The freely draining case is the usual situation which allows
the free motion of the dumbbells through the turbulent flow volume. The Reynolds number is Re = 400.
PSfrag replacements
St = 0.0005
St = 0.005
St = 0.05
St = 0.5
1.0 · 10−4
2.0 · 10−4
3.0 · 10−4
4.0 · 10−4
5.0 · 10−4
2.5 · 10−1
2.0 · 10−1
1.5 · 10−1
1.0 · 10−1
0.5 · 10−1
9.0 · 10−2
8.0 · 10−2
7.0 · 10−2
6.0 · 10−2
5.0 · 10−2
4.0 · 10−2
7.0 · 10−1
8.0 · 10−1
9.0 · 10−1
9.0 · 10−1
8.5 · 10−1
8.0 · 10−1
7.5 · 10−1
7.0 · 10−1
FIG. 3: Reynolds stresses 〈u′iu
j〉 normalized by the turbulent kinetic energy k = 〈(u
2〉/2 as a function of the Weissenberg
and Stokes numbers. From top to bottom: streamwise fluctuations, fluctuations in shear direction, spanwise fluctuations, and
shear stress. The Reynolds number is Re = 800.
10−2 10−1
1.2 Wi = 20
Wi = 100
FIG. 4: Ratio of friction factors as a function of St for the largest Wi. The friction factor for the fluid with the dispersed
dumbbells is cf (cf. Eq. (27)). The quantity c
f is the friction factor of the Newtonian fluid.
FIG. 5: Isosurface plot of the fluctuations of the streamwise turbulent velocity component u′x. The snapshots are for Re = 800
and St = 0.0005. The isolevels are for ±0.04 in each case.
Wi = 3
St = 0.0005
St = 0.005
St = 0.05
St = 0.5
0.0 0.2 0.4 0.6 0.8 1.0
Wi = 20
St = 0.0005
St = 0.005
St = 0.05
St = 0.5
0.0 0.2 0.4 0.6 0.8 1.0
Wi = 100
St = 0.0005
St = 0.005
St = 0.05
St = 0.5
0.0 0.2 0.4 0.6 0.8 1.0
FIG. 6: Probability density function (PDF) of the extension R normalized by the contour length L0. Three different Weissenberg
numbers are shown. The Stokes numbers of the data are indicated in the legend. Data are for Re = 800.
PSfrag replacements
Wi = 3
Wi = 100
−3π/8
St = 0.0005
St = 0.0005
St = 0.005
St = 0.005
St = 0.05
St = 0.05
St = 0.5
St = 0.5
FIG. 7: Asymmetry of the probability density function (PDF) of the azimuthal angle ϕ. It is defined as A(ϕ) = p(ϕ)− p(−ϕ).
The upper panel shows the data for Wi = 3 and four different Stokes numbers. The lower panel shows the data for Wi = 100
and four different Stokes numbers. The analysis is for Re = 800.
FIG. 8: Sketch of the orientation of a dumbbell in the turbulent shear flow. The mean turbulent flow profile is indicated.
The dark-colored dumbbell stands for the preferentially oriented one while the gray-colored orientation is less probable. This
orientation asymmetry leads to the asymmetry in the angular distribution as given in Fig. 7.
19PSfrag replacements
x/(∂yu
x)rms
x/(∂yu
x)rms
Newtonian
Newtonian
Wi = 3
Wi = 20
Wi = 20
Wi = 100
Wi = 100
St = 0.0005
St = 0.5
FIG. 9: Probability density function (PDF) of the transverse velocity gradient of the streamwise turbulent fluctuations, ∂u′x/∂y.
The Newtonian case is compared with the two larger values of the Weissenberg number at St = 5× 10−4 and 0.5, respectively.
The data are for Re = 800.
FIG. 10: Relation between the extension of the dumbbells and the local velocity gradient at the center of mass of the dumbbells.
The local flow topology that is related to the velocity gradient is quantified by the second and third invariants I2 and I3 (see
eqns. (34) for the definition). Quadrant I stands for vortex stretching, II for vortex compression, III for bi-axial strain, and IV
for uniaxial strain, respectively. The gray color coding of the bins for 0 < R/L0 < 0.25, 0.25 ≤ R/L0 < 0.5, 0.5 ≤ R/L0 < 0.75,
0.75 ≤ R/L0 ≤ 1 is indicated by the legend for each figure. Data are for Re = 800.
	Introduction
	Model and equations
	The Newtonian solvent
	The FENE dumbbells
	Two-way coupling
	Large-scale properties
	Energy balance
	Reynolds stresses
	Small-scale properties
	Extensional and angular statistics of dumbbells
	Velocity gradient statistics
	Invariants of the velocity gradient tensor and dumbbell extension
	Summary and discussion
	Acknowledgments
	Semi-implicit integration scheme for dumbbells
	References

  We present numerical studies for finitely extensible nonlinear elastic (FENE)
dumbbells which are dispersed in a turbulent plane shear flow at moderate
Reynolds number. The polymer ensemble is described on the mesoscopic level by a
set of stochastic ordinary differential equations with Brownian noise. The
dynamics of the Newtonian solvent is determined by the Navier-Stokes equations.
Momentum transfer of the dumbbells with the solvent is implemented by an
additional volume forcing term in the Navier-Stokes equations, such that both
components of the resulting viscoelastic fluid are connected by a two-way
coupling. The dynamics of the dumbbells is given then by Newton's second law of
motion including small inertia effects. We investigate the dynamics of the flow
for different degrees of dumbbell elasticity and inertia, as given by
Weissenberg and Stokes numbers, respectively. For the parameters accessible in
our study, the magnitude of the feedback of the polymers on the macroscopic
properties of turbulence remains small as quantified by the global energy
budget and the Reynolds stresses. A reduction of the turbulent drag by up to
20% is observed for the larger particle inertia. The angular statistics of the
dumbbells shows an increasing alignment with the mean flow direction for both,
increasing elasticity and inertia. This goes in line with a growing asymmetry
of the probability density function of the transverse derivative of the
streamwise turbulent velocity component. We find that dumbbells get stretched
referentially in regions where vortex stretching or bi-axial strain dominate
the local dynamics and topology of the velocity gradient tensor.

Introduction to polymer physics, Cambridge University Press, Cambridge, 1996.
[7] R. Benzi, E. De Angelis, V. S. L‘vov, I.Procaccia, and V. Tiberkevich, “Maximum drag reduction asymptotes and the
cross-over to the Newtonian plug,” J. Fluid Mech. 551, 185 (2006).
[8] M. D. Warholic, H. Massah, and T. J. Hanratty, “Influence of drag-reducing polymers on turbulence: effects of Reynolds
number, concentration and mixing,” Exp. Fluids 27, 461 (1999).
[9] K. R. Sreenivasan and C. M. White, “The onset of drag reduction by dilute polymer additives and the maximum drag
reduction asymptote,” J. Fluid Mech. 409, 149 (2000).
[10] C. M. White, V. S. R. Somandepalli, and M. G. Mungal, “The turbulence structure of drag reduced boundary layer flow,”
Exp. Fluids 36, 62 (2004).
[11] R. B. Bird, R. C. Amstrong, and O. Hassager, Dynamics of polymeric liquids, John Wiley & Sons, New York, 1987.
[12] R. Sureshkumar, A. N. Beris, and R. A. Handler, “Direct numerical simulation of the turbulent channel flow of a polymer
solution,” Phys. Fluids 9, 743 (1997).
[13] P. Ilg, E. De Angelis, I. V. Karlin, C. M. Casciola, and S. Succi, “Polymer dynamics in wall turbulent flow,” Europhys.
Lett. 58, 616 (2002).
[14] B. Eckhardt, J. Kronjäger, and J. Schumacher, “Stretching of polymers in a turbulent enviroment,” Comp. Physics Comm.
147, 538 (2002).
[15] C. D. Dimitropolous, Y. Dubief, E. S. G. Shaqfeh, P. Moin, and S. K. Lele, “Direct numerical simulation of polymer-induced
drag reduction in a turbulent shear flow,” Phys. Fluids 17, 011705 (2005).
[16] T. Vaithianathan and L. R. Collins, “Numerical approach to simulating turbulent flow of a viscoelastic polymer solution,”
J. Comp. Phys. 187, 1 (2003).
[17] C. R. Doering, B. Eckhardt, and J. Schumacher, “Failure of energy stability in Oldroyd-B fluids at arbitrarily low Reynolds
numbers,” J. Non-Newtonian Fluid Mech. 135, 92 (2006).
[18] H. Temmen, H. Pleiner, M. Liu, and H. R. Brand, “Convective nonlinearity in non-Newtonian fluids,” Phys. Rev. Lett.
84, 3228 (2000).
[19] A. N. Beris, M. D. Graham, I. Karlin, and H.-C. Öttinger, “Comment on Convective nonlinearity in non-Newtonian fluids,”
Phys. Rev. Lett. 86, 744 (2001).
[20] H. Temmen, H. Pleiner, M. Liu, and H. R. Brand, “Temmen et al. reply,” Phys. Rev. Lett. 86, 745 (2001).
[21] H.-C. Öttinger, Stochastic processes in polymeric fluids, Springer Verlag, Berlin, 1996.
[22] A. Celani, A. Puliafito, and K. Turitsyn, “Polymers in linear shear flow: a numerical study,” Europhys. Lett. 70, 464
(2005).
[23] J. S. Hur, E. S. G. Shaqfeh, H. P. Babcock, D. E. Smith, and S. Chu, “Dynamics of dilute and semidilute DNA solutions
in the start-up of shear flow,” J. Rheol. 45, 421 (2001).
[24] P. A. Stone and M. D. Graham, “Polymer dynamics in a model of the turbulent buffer layer,” Phys. Fluids 15, 1247 (2003).
[25] V. E. Terrapon, Y. Dubief, P. Moin, E. S. G. Shaqfeh, and S. K. Lele, “Simulated polymer stretch in a turbulent flow
using Brownian dynamics,” J. Fluid Mech. 504 61 (2004).
[26] B. H. Zimm, “Dynamics of polymers in dilute solution - viscoelasticity, birefringence and dielectric loss,” J. Chem. Phys.
24, 269 (1956).
[27] A. Celani, S. Musacchio, and D. Vincenzi, “Polymer transport in random flow,” J. Stat. Phys. 118, 531 (2005).
[28] J. Davoudi and J. Schumacher,“Stretching of polymers around the Kolmogorov scale in a turbulent shear flow,” Phys.
Fluids 18, 025103 (2006).
[29] K. D. Squires and J. K. Eaton, “Particle response and turbulence modification in isotropic turbulence,” Phys. Fluids A 2,
1191 (1990).
[30] S. Elghobashi and G. C. Truesdell, “On the two-way interaction between homogeneous turbulence and dispersed solid
particles. I: Turbulence modification,” Phys. Fluids A 5, 1790 (1993).
[31] T. Bosse, L. Kleiser, and E. Meiburg, “Small particles in homogeneous turbulence: Settling velocity enhancement by
two-way coupling,” Phys. Fluids 18, 027102 (2006).
[32] I. M. Mazzitelli, D. Lohse, and F. Toschi, “On the relevance of the lift force in bubbly turbulence,” J. Fluid Mech. 488,
283 (2003).
[33] P. Ahlrichs and B. Dünweg, “Simulation of a single polymer chain in solution by combining lattice Boltzmann and molecular
dynamics,” J. Chem. Phys. 111, 8225 (1999).
[34] J. D. Schieber and H.-C. Öttinger, “The effect of bead inertia on the Rouse model,” J. Chem. Phys. 89, 6972 (1988).
[35] J. Schumacher and B. Eckhardt, “Evolution of turbulent spots in a plane shear flow,” Phys. Rev. E 63, 046307 (2001).
[36] P. J. Flory, Principles of Polymer Chemistry, Cornell University Press, Ithaca, 1953.
[37] G. Boffetta, A. Celani, and A. Mazzino, “Drag reduction in the turbulent Kolmogorov flow”, Phys. Rev. E 71, 036307
(2005).
[38] H. T. Schlichting, Boundary Layer Theory, McGraw Hill, New York, 1979.
[39] F. Waleffe, “On a self-sustaining process in shear flows,” Phys. Fluids 9, 883 (1997).
[40] S. Grossmann, “The onset of shear flow turbulence,” Rev. Mod. Phys. 72 , 603 (2000).
[41] P. A. Stone, A. Roy, R. G. Larson, F. Waleffe, and M. D. Graham, “Polymer drag reduction in exact coherent structures
of plane shear flow,” Phys. Fluids 16, 3470 (2004).
[42] Y. Dubief, C. M. White, V. Terrapon, E. S. G. Shaqfeh, P. Moin, and S. K. Lele, “On the coherent drag-reducing and
turbulence-enhancing behaviour of polymers in wall-flows,” J. Fluid Mech. 514, 271 (2004).
[43] S. Gerashchenko, C. Chevallard, and V. Steinberg, “Single polymer dynamics: coil-stretch transition in a random flow,”
Europhys. Lett. 71 221 (2005).
[44] M. Chertkov, I. Kolokolov, V. Lebedev, and K. Turitsyn, “Polymer statistics in a random flow with mean shear,” J. Fluid
Mech. 531, 251 (2005).
[45] Z. Warhaft, “Turbulence in nature and laboratory,” Proc. Nat. Acad. Sci. 99, 2481 (2002).
[46] J. Schumacher, K. R. Sreenivasan, and P. K. Yeung, “Derivative moments in turbulent shear flows,” Phys. Fluids 15, 84
(2003).
[47] L. Biferale and I. Procaccia, “Anisotropy in turbulent flows and in turbulent transport,” Phys. Rep. 414, 43 (2005).
[48] A. Liberzon, M. Guala, B. Lüthi, W. Kinzelbach, and A. Tsinober, “Turbulence in dilute polymer solutions,” Phys. Fluids
17, 031707 (2005).
[49] A. Liberzon, M. Guala, W. Kinzelbach, and A. Tsinober, “On the kinetic energy production and dissipation in dilute
polymer solutions,” Phys. Fluids 18, 125101 (2006).
[50] P. A. Davidson, Turbulence, Oxford University Press, Oxford, 2004.
[51] A. Celani, A. Puliafito, and D. Vincenzi, “Dynamical slowdown of polymers in laminar and random flows,” Phys. Rev.
Lett. 97, 118301 (2006).
[52] T. Needham, Visual Complex Analysis, Oxford University Press, Oxford, 1997.
[53] We will use the definitions as given in Ref. [50] (p. 264) with I1 = P , I2 = Q, and I3 = R
PSfrag replacements
Wi = 20
Wi = 20
Wi = 100
Wi = 100
Fel = 0
Fel = 0
FIG. 1: Mean dissipation and injection rates as a function of the Stokes and Weissenberg numbers. Upper picture: mean energy
injection rate due to shear flow forcing 〈εin〉. Mid panel: mean energy dissipation rate 〈ε〉. Lower panel: mean dissipation rate
which arises from the coupling to the dumbbell ensemble 〈εp〉. The Reynolds number is Re = 800. The case with Fel = 0 is for
the case without spring force and stands for a shear flow with dispersed inertial particles.
PSfrag replacements
Wi = 20, St = 0.05
tethered
0 5 10 15 20
FIG. 2: Comparison of the kinetic energy for two cases at Wi = 20 and St = 0.05: for the tethered case one bead of each
dumbbell is fixed at a grid site while the second bead can fluctuate. The freely draining case is the usual situation which allows
the free motion of the dumbbells through the turbulent flow volume. The Reynolds number is Re = 400.
PSfrag replacements
St = 0.0005
St = 0.005
St = 0.05
St = 0.5
1.0 · 10−4
2.0 · 10−4
3.0 · 10−4
4.0 · 10−4
5.0 · 10−4
2.5 · 10−1
2.0 · 10−1
1.5 · 10−1
1.0 · 10−1
0.5 · 10−1
9.0 · 10−2
8.0 · 10−2
7.0 · 10−2
6.0 · 10−2
5.0 · 10−2
4.0 · 10−2
7.0 · 10−1
8.0 · 10−1
9.0 · 10−1
9.0 · 10−1
8.5 · 10−1
8.0 · 10−1
7.5 · 10−1
7.0 · 10−1
FIG. 3: Reynolds stresses 〈u′iu
j〉 normalized by the turbulent kinetic energy k = 〈(u
2〉/2 as a function of the Weissenberg
and Stokes numbers. From top to bottom: streamwise fluctuations, fluctuations in shear direction, spanwise fluctuations, and
shear stress. The Reynolds number is Re = 800.
10−2 10−1
1.2 Wi = 20
Wi = 100
FIG. 4: Ratio of friction factors as a function of St for the largest Wi. The friction factor for the fluid with the dispersed
dumbbells is cf (cf. Eq. (27)). The quantity c
f is the friction factor of the Newtonian fluid.
FIG. 5: Isosurface plot of the fluctuations of the streamwise turbulent velocity component u′x. The snapshots are for Re = 800
and St = 0.0005. The isolevels are for ±0.04 in each case.
Wi = 3
St = 0.0005
St = 0.005
St = 0.05
St = 0.5
0.0 0.2 0.4 0.6 0.8 1.0
Wi = 20
St = 0.0005
St = 0.005
St = 0.05
St = 0.5
0.0 0.2 0.4 0.6 0.8 1.0
Wi = 100
St = 0.0005
St = 0.005
St = 0.05
St = 0.5
0.0 0.2 0.4 0.6 0.8 1.0
FIG. 6: Probability density function (PDF) of the extension R normalized by the contour length L0. Three different Weissenberg
numbers are shown. The Stokes numbers of the data are indicated in the legend. Data are for Re = 800.
PSfrag replacements
Wi = 3
Wi = 100
−3π/8
St = 0.0005
St = 0.0005
St = 0.005
St = 0.005
St = 0.05
St = 0.05
St = 0.5
St = 0.5
FIG. 7: Asymmetry of the probability density function (PDF) of the azimuthal angle ϕ. It is defined as A(ϕ) = p(ϕ)− p(−ϕ).
The upper panel shows the data for Wi = 3 and four different Stokes numbers. The lower panel shows the data for Wi = 100
and four different Stokes numbers. The analysis is for Re = 800.
FIG. 8: Sketch of the orientation of a dumbbell in the turbulent shear flow. The mean turbulent flow profile is indicated.
The dark-colored dumbbell stands for the preferentially oriented one while the gray-colored orientation is less probable. This
orientation asymmetry leads to the asymmetry in the angular distribution as given in Fig. 7.
19PSfrag replacements
x/(∂yu
x)rms
x/(∂yu
x)rms
Newtonian
Newtonian
Wi = 3
Wi = 20
Wi = 20
Wi = 100
Wi = 100
St = 0.0005
St = 0.5
FIG. 9: Probability density function (PDF) of the transverse velocity gradient of the streamwise turbulent fluctuations, ∂u′x/∂y.
The Newtonian case is compared with the two larger values of the Weissenberg number at St = 5× 10−4 and 0.5, respectively.
The data are for Re = 800.
FIG. 10: Relation between the extension of the dumbbells and the local velocity gradient at the center of mass of the dumbbells.
The local flow topology that is related to the velocity gradient is quantified by the second and third invariants I2 and I3 (see
eqns. (34) for the definition). Quadrant I stands for vortex stretching, II for vortex compression, III for bi-axial strain, and IV
for uniaxial strain, respectively. The gray color coding of the bins for 0 < R/L0 < 0.25, 0.25 ≤ R/L0 < 0.5, 0.5 ≤ R/L0 < 0.75,
0.75 ≤ R/L0 ≤ 1 is indicated by the legend for each figure. Data are for Re = 800.
	Introduction
	Model and equations
	The Newtonian solvent
	The FENE dumbbells
	Two-way coupling
	Large-scale properties
	Energy balance
	Reynolds stresses
	Small-scale properties
	Extensional and angular statistics of dumbbells
	Velocity gradient statistics
	Invariants of the velocity gradient tensor and dumbbell extension
	Summary and discussion
	Acknowledgments
	Semi-implicit integration scheme for dumbbells
	References
