Maximum drag reduction simulation using rodlike polymers

Maximum drag reduction simulation using rodlike polymers

Department of Chemical Engineering, Delft University of Technology, Julianalaan 136, 2628 BL Delft, the Netherlands (Received 27 April 2012; revised manuscript received 21 September 2012; published 5 October 2012) Simulations of maximum drag reduction (MDR) in channel flow using constitutive equations for suspensions of noninteracting rods predict a few-fold larger turbulent kinetic energy than in experiments using rodlike polymers. These differences are attributed to the neglect of interactions between polymers in the simulations. Despite these inconsistencies the simulations correctly reproduce the essential features of MDR, with universal profiles of the mean flow and the shear stress budgets that do not depend on the polymer concentration.

DOI:10.1103/PhysRevE.86.046304 PACS number(s): 47.27.ek, 47.27.nd, 47.85.lb

I. INTRODUCTION

Polymer-induced drag reduction is the phenomenon wherein the friction factor of a turbulent boundary layer is reduced by the addition of long-chain, linear polymers; see, for example, Ref. [1]. The key property of a drag-reducing polymer is an extended shape. Rigid, rodlike polymers have this shape, while flexible polymers are randomly coiled in fluid at rest and need to be extended by sufficient fluid strain to become effective drag reducers [2].

Remarkably, the amount of drag reduction is bounded by an upper limit. When this maximum drag reduction (MDR) is reached, the flow does not change upon adding more polymers. In fact, MDR is not only independent of polymer concentration but also of polymer properties, such as chain length and flexibility; see, for example, Refs. [3–5].

In this paper we show that a numerical model of suspensions of noninteracting rods correctly reproduces the essential features of MDR: concentration-invariant profiles of mean flow and shear stress budgets. This means that the essential properties of MDR are recovered regardless of certain details of the polymer dynamics, such as the interactions between the polymers.

In this work we assess our MDR simulations by comparing the profiles of the mean flow, the Reynolds stress, and the turbulent kinetic energy (TKE) to experiments from the literature.

II. MATHEMATICAL MODEL A. Fluid equations of motion

We use direct numerical simulation (DNS) to compute turbulent channel flow of solutions of rodlike polymers of length l and diameter d, which are homogeneously dissolved at a volume concentration c in a solvent with a kinematic viscosity ν and a mass density ρ. The rods are inertialess, buoyantly free, and smaller than the Kolmogorov length scale but large enough such that the effects of Brownian motion can be neglected. Also, the rod aspect ratio r= l/d is assumed to be large enough, such that effects of a finite r can be ignored. The channel geometry is sketched in Fig.1. A fluid is driven by means of a constant pressure gradient−d/dx between two parallel, no-slip walls, separated by a distance D, in the y direction. The wall-parallel directions are referred to as streamwise (x) and spanwise (z). The velocity components in the streamwise, wall-normal, and spanwise directions are

denoted u, v, and w. Furthermore, subscripts x, y, and z or 1, 2, and 3 denote components of vectors and tensors. The dimensions of the channel are DLx and DLz in the x and z

directions. When denoted with the superscript +, a variable is scaled with density ρ, viscosity ν and friction velocity Uτ = [−(d/dx)(D/ρ)(1/4)]1/2.

We numerically integrate the incompressible, Navier-Stokes equations, supplemented by the divergence of the polymer stress tensor:

ρDu

Dt = ∇ · (−δ + 2μS + τ), (1a)

∇ · u = 0. (1b)

Here u is the fluid velocity vector, t is time, ∇ is the nabla operator, δ is the unit tensor, D/Dt = ∂/∂t + u · ∇ is the material derivative, S= 1₂(∇uT + ∇u) is the rate of strain tensor,  is the pressure, μ= νρ is the solvent dynamic viscosity, and τ is the polymer stress tensor, for which additional equations have to be solved.

Since there are no numerically feasible models available for τ that take interactions into account rigorously, we neglect these interactions in the present simulations, similarly as in previously reported simulations of (maximum) drag reduction; see, for example, Refs. [4,6]. Without interactions the expression forτ can be derived rigorously [7]:

τ = 2αμS :  pppp. (1c)

Here · · ·  is an average over the polymers, contained in a volume, which is centered at the point where the stress is to be determined and which is small compared to the Kolmogorov length scale, p is the polymer orientation unit vector, and α is the polymer concentration parameter, which measures the strength of the polymer stress relative to the Newtonian stress:

α≈ 0.1cr2. (1d)

Here we have ignored the logarithmic dependence of the numerical factor 0.1 on r. The scaling of α with cr2 _means that a polymer of length l makes a contribution to the stress, which is of the same order of magnitude as that of a sphere with diameter l. This explains that significant drag reduction effects have been observed using very small polymer concentrations, c∼ 10−5, provided that the aspect ratio is very large, r ∼ 104_; see, for example, Ref. [8].

As given by Eq.(1c), the polymer stress depends on the fourth-order moment pppp of the orientational distribution

x, u

z, w

y, v

−

D

L

_D

L

D

FIG. 1. Channel geometry. The coordinate axes are referred to as streamwise x, wall-normal y, and spanwise z. The velocity components in x, y, and z are denoted u, v, and w.

function. We approximate this quantity by solving the transport equation for the second-order moment pp [7]:

D p p

Dt − ∇u

T _{·  p p −  p p · ∇u}

= −2∇u :  pppp + β∇2_{ pp, (1e)} and expressing the fourth-order moment in terms of the second-order moment by means of the closure proposed by Cintra and Tucker [9]. The diffusive term β∇2_{ pp in Eq.(1e)} can be considered a model for the effect of the unresolved scales of pp, which are a few orders of magnitude smaller than the Kolmogorov scale [10]. In previous work we have shown that the fourth-order moment closure, as well as the subgrid model, do not introduce significant errors [10,11]. In order to guarantee a numerically stable solution, we choose β+ to range from 2 for small drag reduction up to 8 for maximum drag reduction. These values are comparable to previously used values in MDR simulations using flexible polymers [4,6]. Ptasinski et al. demonstrated that in this range the numerical solution for drag-reduced channel flow does not depend on the exact value of β [4].

As previously stated, we ignore the interactions between the polymers in the simulation. When classifying the interactions, it is customary to consider concentration regimes. In the dilute regime, the average distance s between the polymers is larger than the polymer length l, which means that hydrodynamic interactions can be neglected. Diluteness requires cr2_{ 1 and} this means that the polymer stress [Eq. (1c)] is negligible.

Therefore in this regime the flow is Newtonian, without any drag reduction. In the semidilute regime, the distance s is smaller than l but still larger than d, which means that hydrodynamic interactions are important while physical contacts are still rare. This regime corresponds to cr  1  cr2, which means that α 1, so that drag reduction can be expected. Finally, in the concentrated regime, s  d and the polymers are in constant physical contact with each other, which corresponds to cr  1.

We now estimate the number of hydrodynamic interactions per polymer∼cr2_{at MDR. To this end we use that the elastic} sublayer thickness δE+ at MDR equals the pipe radius [8].

Second, we use that for rodlike polymers [12]

δ_E+≈ 0.2α. (2)

We thus deduce that at MDR each polymer interacts with cr2∼ 102Reτ neighboring polymers, where we have used

that cr2≈ 10α and α ≈ 5δE+ and δ+E≈ Reτ/2. Typically, in

drag reduction experiments Reτ ∼ 102–103, so we find that at

MDR, each polymer interacts with cr2_{∼ 10}4_–105_neighboring polymers. It is noted that this estimate applies equally well to flexible polymers. This is because Eq. (2) also applies to flexible polymers, given that r in the expression for α [Eq.(1d)] is based on the average polymer extension near the wall [13].

These estimates suggest that interactions between polymers play an important role at MDR. We show, however, that the essential features of MDR are correctly simulated even when interactions between polymers are neglected in the numerical model.

III. NUMERICAL METHODS AND PARAMETERS According to Eqs. (1), the simulated flow is governed by two dimensionless parameters: the polymer concentration parameter α and the frictional Reynolds number Reτ =

UτD/ν. We have carried out seven simulations, where we

used Reτ= 1000 and varied α from zero (Newtonian flow)

up to 1600 (maximum drag reduction). Additional numerical parameters are listed in Table I. For the Newtonian flow calculation, the streamwise and spanwise channel dimensions are 1000 and 500 wall units, which are sufficient to capture the near-wall vortical structures responsible for the friction factor. These structures increase in size with increasing drag reduction

TABLE I. Parameters used in the simulations. The markers correspond to Figs.3(a),4(a),5(a),6and7. Re is the bulk Reynolds number. α is the concentration parameter [Eq.(1d)]. Lxis the ratio of the streamwise and wall-normal

channel dimensions. Lzis the ratio of the spanwise and wall-normal channel dimensions. Nx, Ny, and Nzare the

number of grid points in streamwise, wall-normal, and spanwise directions. β is the artificial diffusivity in Eq.(1e). DR is drag reduction [Eq.(3)].

RUN Marker 10−3Re α Lx= 2Lz Nx = Nz Ny β+ DR 1 ◦ 18 0 1 96 384 2  24 50 2 96 192 2 0.42 3  26 100 3 96 192 4 0.53 4 32 300 6 96 96 4 0.68 5  39 600 15 192 96 7 0.77 6  42 800 18 192 96 8 0.82 7  43 1600 18 192 96 8 0.81

102 103 0.4 0.6 0.8

α

DR

FIG. 2. Simulated drag reduction [Eq.(3)] vs polymer concen-tration α [Eq.(1d)].

[14]. Therefore, we have used larger channel dimensions for the drag-reduced flow calculations. Extensive tests have shown that for the current parameters the results hardly depend on grid resolutions and domain sizes.

Spatial derivatives are computed with a Fourier basis for the homogeneous directions and a second-order, central, finite-differences scheme for the wall-normal direction. Time integration is achieved with the second-order, explicit Adams-Bashforth scheme. Conservation of mass is ensured using a standard projection method.

IV. RESULTS A. Drag reduction Drag reduction (DR) is defined as

DR= 1 − f

P

fN, (3)

where the friction factor f = D(d/dx)/(ρU2

b) for the

polymer solution fP _{and for the Newtonian fluid f}N _are

evaluated at the same frictional Reynolds number Reτ. Here

Ub is the bulk velocity. Figure 2 shows DR as a function

of the polymer concentration parameter α. Two regimes are observed. For α < 800 we are in the small drag reduction (SDR) regime, where DR scales as log(α). In Ref. [12] we rationalized this logarithmic scaling by showing that the elastic layer thickness δE+ scales linearly with α. For α > 800, DR

reaches a plateau value, which defines the maximum drag reduction (MDR) regime. At MDR the elastic layer spans the whole channel cross section.

This is a demonstration of a concentration invariant friction factor using direct numerical simulation. Besides concentra-tion invariance, we also investigate the profiles of mean flow, turbulent energy, and shear stress budgets in the next sections.

B. Mean flow

We compare our simulation results to the experimental data of Escudier et al. [5], who conducted channel flow experiments using aqueous solutions of xanthan gum (XG). Owing to the high rigidity of its backbone structure this polymer can be considered a rigid, rodlike particle. We consider three of their data sets, corresponding to Newtonian flow, SDR, and MDR. For these data sets, the polymer mass concentration and the frictional Reynolds number (cm; Reτ)

100 101 102 0 20 40 60

y

u

y

+

u

+

FIG. 3. (a) Simulated velocity profiles. The simulation param-eters are given in Table I. For clarity the markers are shown for only one tenth of the grid points. (b) Measured velocity profiles in Newtonian flow (circles) and drag-reduced flow using solutions of xanthan gum at concentration per weight cm= 6.7 × 10−4(squares)

and cm= 1.5 × 10−3(triangles). Data taken from Ref. [5].

are approximately (0; 6× 102), (6.7× 10−4; 7× 102), and (1.5× 10−3; 5× 102_{), respectively.}

The measured velocity profiles are shown in Fig. 3(b). These profiles extend from the wall at y+= 0 to the channel centerline, which in+-units equals half the frictional Reynolds number. The Newtonian profile with cm= 0 (circles) follows

the well-known law of the wall: a linear viscous region and a logarithmic inertial region, with a 2.5 slope. In the polymer solution with cm= 6.7 × 10−4 (squares), the logarithmic

profile is shifted upward in parallel. In the more concentrated solution cm= 1.5 × 10−3(downward triangles) the flow is at

MDR, with a logarithmic profile with a 11.7 slope, known as Virk’s asymptote [8].

The computed velocity profiles are shown in Fig.3(a). With increasing concentration α the profiles change in accordance with the experimental data in Fig.3(b). The logarithmic profile in the inertial layer is shifted upward in parallel, due to the formation of an elastic layer between the viscous and inertial layers, that approaches Virk’s asymptote. When the elastic layer spans the whole channel cross section the flow is at MDR, which corresponds to our simulations with α= 800 (left-pointing triangles) and α= 1600 (downward triangles). Despite some slight variations in the computed slope, our simulated MDR profiles for α= 800 and α = 1600 are in reasonable agreement with Virk’s empirical asymptote, which is plotted with a dashed line in Fig. 3(a). The simulated velocity profiles for α= 0 (circles), 100 (squares), and

1600 (downward triangles) in Fig. 3(a) have similar levels of drag reduction as the measured profiles in Fig. 3(b). Therefore, the following comparisons are based on these three simulations:

It is noted that the Reynolds numbers in the experiments are different than in the corresponding simulations, which may cause some quantitative differences between the two. However, the experiments and the simulations are operated in the same flow regimes: Newtonian, SDR, and MDR. Therefore, the qualitative trends should be similar, and any qualitative differences would reflect shortcomings of the numerical model. We have excluded the fourth-order moment closure, the subgrid model, the domain sizes, and the grid resolutions as significant sources of error [10,11]. We therefore believe that any discrepancies would be due to the neglect of polymer interaction in the numerical model.

C. Turbulent kinetic energy

Next we consider the turbulent kinetic energy (TKE), mea-sured by the standard deviation of the fluid velocity vector urms. Simulation and experiment are compared in Fig.4, showing the profiles of urms for the streamwise and the wall-normal velocity components. Most of the energy is contained in the

100 101 102 0 2 4 6 (a)

y

+

u

+

_rms

y

+

u

+

_rms

FIG. 4. (a) Simulated standard deviation of streamwise fluid velocity (open markers) and wall-normal fluid velocity (filled markers). The simulation parameters are given in Table I. For clarity the markers are shown for only one tenth of the grid points. (b) Measured standard deviation of streamwise fluid velocity (open markers) and wall-normal fluid velocity (filled markers) in Newtonian flow (circles) and drag-reduced flow using solutions of xanthan gum at concentration per weight cm= 6.7 × 10−4 (squares) and

cm= 1.5 × 10−3(triangles). Data taken from Ref. [5].

streamwise velocity component. The anisotropy is largest near the wall, while in the channel center the TKE is more evenly distributed over the different directions. For Newtonian flow the kinetic energy peaks at y+≈ 10 and decreases toward the channel center. In drag-reduced flow, this peak shifts away from the wall. The Newtonian simulation (circles) is in perfect agreement with the corresponding experiment. The SDR simulation (squares) is in qualitative agreement with the corresponding experiment. In both simulation as experiment u_rms is increased and vrms is reduced as compared to the Newtonian flow. Quantitatively however, the SDR simulation predicts a larger urmsand a smaller vrms as compared to the experiment. However, these differences may be attributed to differences in Reτ and DR between simulation and

experiment.

For MDR (downward triangles) the agreement is poor, with large differences in urms. In contrast to the experiment, the MDR simulation shows a substantial increase in urms as compared to the Newtonian flow, reaching a maximum of approximately 5.5, being two times as large as the experimental value. Also, the peak of the urmsprofile is shifted further outwards, as compared to the experiment. Previous MDR simulations showed a similarly overpredicted urms[4,6].

D. Shear stress

Finally we consider the various components of the shear stress balance. This balance is obtained by integrating the x

0 0.1 0.2 0.3 0.4 0 0.2 0.4 0.6 0.8 1 (a)

y/D

y/D

FIG. 5. (a) Simulated Reynolds shear stress−uv+(open mark-ers) and polymer shear stress τxy+ (filled markers). The simulation

parameters are given in TableI. For clarity the markers are shown for only one tenth of the grid points. (b) Measured Reynolds shear stress−uv+(open markers) in drag-reduced flow using solutions of xanthan gum at concentration per weight cm= 6.7 × 10−4(squares)

component of the Reynolds averaged Navier-Stokes equations over the y direction. In nondimensional form the result reads 1−2y D = du+ dy+ − uv++ τxy +_{. (4)}

The left-hand side represents the driving pressure gradient force, which balances the shear stress on the right-hand side. The shear stress is divided into the viscous stress du+/dy+, the Reynolds stress−uv+, and the polymer stress τxy+. Drag

reduction is equivalent to an increase of the nondimensional velocity gradient du+/dy+. According to Eq. (4) this is realized when the polymer stress is overwhelmed by the reduction of the Reynolds stress. In Fig.5(a) we show the computed Reynolds stresses and polymer stresses. These results are compared to the experimental data, which is given in Fig.5(b). For both SDR (squares) and MDR (downward triangles), there is good agreement between the simulations and the experiments. It is especially noteworthy that despite the large differences for TKE, there is very good agreement for Reynolds stress at MDR.

We rewrite the shear stress balance [Eq.(4)] into

1−2y D = du+ dy+(1+ ν + E+ νP+), (5) 100 101 102 10−4 10−3 10−2 10−1 100 101 102

y

y

FIG. 6. (a) Simulated eddy viscosity νE+ (black lines). The

molecular viscosity ν+= 1 is shown with the horizontal line and the polymer viscosity νP+is shown with the grey lines, for reference.

(b) Simulated polymer viscosity ν_P+ (black lines). The molecular viscosity ν+= 1 is shown with the horizontal line and the eddy viscosity ν_E+ is shown with the grey lines, for reference. The simulation parameters are given in TableI. For clarity the markers are shown for only one tenth of the grid points.

100 101 102 0 10 20 30

y

k

FIG. 7. Simulated turbulent kinetic energy. The simulation pa-rameters are given in TableI. For clarity the markers are shown for only one tenth of the grid points.

where we introduce the eddy viscosity νE= −uv(du/dy)−1

and the polymer viscosity νP = τxy(du/dy)−1. Figure 6(a)

shows that in the viscous layer (y+<10), the eddy viscosity νEis smaller than ν and it scales as νE∼ y3. For larger y, we

observe the well-known inertial layer scaling νE∼ y. With

the introduction of the polymers, νE decreases and the inner

scaling of νE∼ y3 extends outward, which is accompanied

by the emergence of the elastic layer. In the viscous and elastic layers, the eddy viscosity is observed to decrease as α−1. At large concentrations (α= 800 and α = 1600), the νE profiles collapse on a single curve, which reflects MDR

conditions. Figure6(b)shows that the polymer viscosity scales as νP ∼ y, and for small enough y the data for different α

roughly collapse on a single curve. The linear scaling of νP

and the collapse of νP and νE at MDR are consistent with

the universal, logarithmic velocity profile due to Virk [8]. We thus find that the Reynolds shear stress is independent of α at MDR, as well as the polymer shear stress and the viscous shear stress.

Despite this, the turbulence structure at MDR is not concentration invariant. This is understood from Eq.(1c)as follows. When α is increased, the turbulent velocity gradient ∇u as well as the polymer orientation p must adjust in order to produce a constant τxy= α∇u : pppxpy. In our simulations,

this adjustments is accompanied by an increase in the turbulent kinetic energy k= 1₂u· u, which is plotted in Fig.7.

V. DISCUSSION AND CONCLUSION

Direct numerical simulation has been used to study DR in turbulent channel flow induced by rodlike polymers. At present there is no computationally feasible, rigorous model for polymer entanglements in turbulent flows. Therefore, we ignore these interactions in the simulations, similar to previous DR simulations reported in the literature [4,6].

We compared the DNS to experimental data from the liter-ature, concerning the profiles of the mean flow, the Reynolds, stress and the turbulent kinetic energy. At SDR the trends are in good agreement, while at MDR the simulations produce TKE values which are substantially larger than the experimentally found values. These discrepancies are due to the neglect of interactions between the polymers. A similar situation holds for MDR simulation using flexible chains [4,6] that also

predicts a structurally larger TKE than experimental values [3–5].

Despite this, our MDR simulations correctly reproduce the measured profiles of the mean flow, as well as the concentration invariance of these profiles. Our simulations predict polymer viscosity profiles and eddy viscosity profiles that collapse on asymptotic curves for large polymer concentrations. These results indicates that not only the mean flow but also the shear stress budgets are universal functions at MDR.

We conclude that, even in the absence of polymer interac-tions, our simulations contain sufficient physics to correctly reproduce the mean flow and the shear stress at MDR. These findings are in line with the general notion that MDR is independent of the polymer type and concentration. The details

of the constitutive equations, which depend for instance on the nature of the polymer interactions, have a negligible effect on the essential properties of MDR, while they do have an effect on secondary properties, such as the TKE.

Finally, it is noted that experimental data of Reynolds stress and TKE at MDR are sparse and we recommend further experimental research to characterize these quantities.

ACKNOWLEDGMENT

The research has been partially supported through the PETROMAKS program funded by the Research Council of Norway.

[1] C. M. White and M. G. Mungal,Ann. Rev. Fluid Mech. 40, 235 (2008).

[2] P. Virk,Nature (London) 253, 109 (1975).

[3] M. D. Warholic, H. Massah, and T. J. Hanratty,Exps. Fluids 27, 461 (1999).

[4] P. K. Ptasinski, B. J. Boersma, F. T. M. Nieuwstadt, M. A. Hulsen, B. H. A. van den Brule, and J. C. R. Hunt,J. Fluid Mech. 490, 251 (2003).

[5] M. Escudier, A. Nickson, and R. Poole,J. Non-Newtonian Fluid Mech. 160, 2 (2009).

[6] C. F. Li, R. Sureshkumar, and B. Khomami,J. Non-Newtonian Fluid Mech. 139, 177 (2006).

[7] M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Clarendon, Oxford, 1986).

[8] P. Virk,J. Fluid Mech. 45, 417 (1971).

[9] J. S. Cintra and C. L. Tucker, J. Rheol. 39, 1095 (1995).

[10] J. J. J. Gillissen, B. J. Boersma, P. H. Mortensen, and H. I. Andersson,Phys. Fluids 19, 115107 (2007).

[11] J. J. J. Gillissen, B. J. Boersma, P. H. Mortensen, and H. I. Andersson, Phys. Fluids 19, 035102 (2007).

[12] J. J. J. Gillissen, B. J. Boersma, P. H. Mortensen, and H. I. Andersson,J. Fluid Mech. 602, 209 (2008).

[13] J. J. J. Gillissen, Phys. Rev. E 78, 046311 (2008).

[14] C. F. Li, R. Sureshkumar, and B. Khomami,J. Non-Newtonian Fluid Mech. 140, 23 (2006).