Strong polymer-turbulence interactions in viscoelastic turbulent channel flow
Vassilios Dallas, J. Christos Vassilicos, Geoffrey F. Hewitt
aa r X i v : . [ phy s i c s . f l u - dyn ] N ov Strong polymer-turbulence interactions in viscoelastic turbulent channel flow
V. Dallas ∗ and J. C. Vassilicos † Institute for Mathematical Sciences, Imperial College, London, SW7 2PG, UK andDepartment of Aeronautics, Imperial College, London, SW7 2AZ, UK
G. F. Hewitt
Department of Chemical Engineering and Chemical Technology, Imperial College, London, SW7 2AZ, UK
This paper is focused on the fundamental mechanism(s) of viscoelastic turbulence that lead topolymer induced turbulent drag reduction phenomenon. A great challenge in this problem is thecomputation of viscoelastic turbulent flows, since the understanding of polymer physics is restrictedto mechanical models. An effective state-of-the-art numerical method to solve the governing equa-tion for polymers modelled as non-linear springs, without using any artificial assumptions as usual,was implemented here for the first time on a three-dimensional channel flow geometry. The capa-bility of this algorithm to capture the strong polymer-turbulence dynamical interactions is depictedon the results, which are much closer qualitatively to experimental observations. This allowed amore detailed study of the polymer-turbulence interactions, which yields an enhanced picture on amechanism resulting from the polymer-turbulence energy transfers.
I. INTRODUCTION
A few parts per million by weight polymer in a wall-bounded turbulent flow are enough to reduce the forcenecessary to drive the flow through a channel by a fac-tor of up to 70%, as was discovered by Toms [1] whileperforming experiments on the degradation of polymers.Turbulence is a multiscale phenomenon with a vast spec-trum of spatial scales and therefore a very large num-ber of degrees of freedom. Due to the fact that eventhe maximum polymer molecule end-to-end distance L p is much less than the Kolmogorov viscous scale η , onemight anticipate that the small size polymers can onlyaffect sub-Kolmogorov scale processes and that scales oflength ℓ > η would remain unaffected. Surprisingly, thedynamics of the small polymer chains are able to funda-mentally modify the large scale structures and statistics,as observed in the drag reduction (DR) phenomenon [2].Polymer drag reduction in wall-bounded turbulentflows induces higher mean velocities, implying deviationsfrom the classical phenomenology of hydrodynamic wall-bounded turbulence and hence from the von K´arm´an law U + = 1 κ ln y + + B (1)where ‘+’ denotes normalisation with the friction veloc-ity u τ and the viscous length scale δ ν ≡ ν/u τ with ν being the fluid’s kinematic viscosity. Moreover, κ in Eq.(1) is the von K´arm´an coefficient [3, 4], which is usu-ally considered to be a constant taking the value 0 . B ≃ . et al. [5] distinguished be-tween polymer induced drag reduced flows at low drag re-duction (LDR) and high drag reduction (HDR) regimes, ∗ Electronic address: [email protected] † Electronic address: [email protected] based on the statistical trends of the turbulent velocityfield. When | DR | .
40% (LDR), the mean velocity pro-file is a log-law parallel to the von K´arm´an law Eq. (1)with a higher value of B , i.e. larger mean velocity. How-ever, for 40% < | DR | .
60% (HDR), the slope of thelog-region increases until it reaches the empirical maxi-mum drag reduction (MDR) asymptote U + = 1 κ v ln y + + B v (2)where κ − ≃ . B v ≃ et al. [6, 7] and confirmedexperimentally in channel flow by Warholic et al. [5].Overall, the mean velocity profile is bounded betweenthe von K´arm´an law Eq. (1) and the MDR law Eq. (2),the latter being independent of the Newtonian solvent,the polymer characteristics and the flow geometry.Polymer induced drag reduction has been known formore than sixty years and has attracted attention bothfrom the fundamental and applied perspective. However,no generally accepted theory has been provided to ex-plain adequately the phenomenon. Such a theory shouldprovide an explanation of the drag reduction onset, aswell as the MDR law and its universality, which playsa significant fundamental role in understanding the phe-nomenon. Several theoretical concepts have been pro-posed but all have been subjected to criticism. The pro-posed theories fall mainly into two categories, that ofviscous [8, 9] and that of elastic effects [10–12].Recent progress in DNS of viscoelastic turbulence hasbegun to elucidate some of the dynamical interactionsbetween polymers and turbulence, which are responsiblefor drag reduction. The aim of this study is to investigatethe polymer dynamics, their influence on flow quantitiesand the phenomenology of drag reduction in the variousregimes through DNS of viscoelastic turbulent channelflow using the finite extensible non-linear elastic modelwith the Peterlin linearisation (FENE-P) [15], the mostwidely used coarse-grained model in such studies.The paper is organised as follows. The necessary de-tails on the DNS of viscoelastic turbulent channel floware provided in section II. We analyse various viscoelas-tic turbulent statistics in section III for all the drag re-duction regimes achieved in this study with a state-of-the-art numerical approach, which we have adapted towall-bounded flows [13], aimed at capturing discontinu-ities in the polymer field. This approach is described insome detail in appendix 1. Specifically, the effects ofpolymer extensibility and Reynolds number are brieflyconsidered, whereas the statistics of mean velocity, fluc-tuating velocities and vorticities are examined in depthdemonstrating that our computations are qualitativelycloser to experimental observations than previous numer-ical studies. Section IV presents the conformation tensorstatistics and the scaling of polymer stress tensor com-ponents at the high Weissenberg number limit, whichassists in a new asymptotic result for the shear stressbalance (see section V). Finally, the polymer-turbulenceinteractions are studied in section VI through the energybalance. A refined and extended picture of a conceptualmodel for drag reduction based on viscoelastic dissipa-tion is proposed in section VII before summing up ourmost important results (see section VIII). II. DNS OF VISCOELASTIC TURBULENTCHANNEL FLOW
The enormous number of degrees of freedom of eachcoil means that polymers are an extraordinarily complexsystem, whose dynamics depend on the conformationsof the polymer molecules, i.e. orientation and degreeof stretching of a coil. The study of detailed motionsof this complex system and their relations to the non-equilibrium properties would be prohibitive. Only af-ter elimination of the fast relaxation processes of localmotions in favour of stochastic noise, is it possible tostudy the dynamics of longer relaxation time scales [14],such as the end-to-end conformation, that are responsiblefor many physical properties of polymers in fluids, suchas viscoelastic turbulence and polymer drag reduction.Thus, coarse-grained mechanical models, such as bead-rod-spring models, are very crucial in DNS of viscoelasticturbulence.The computationally demanding Navier-Stokes equa-tions in three-dimensions makes a Lagrangian approachfor the polymer equally prohibitive and also limits poly-mer models to simple representations. A successful modelfor DNS studies of turbulent drag reduction is the FENE-P model in the Eulerian frame of reference, representing aconformation field of polymer macromolecules that havebeen modelled as non-linear bead-spring dumbbells [15].The standard approach to numerically solve the FENE-Pmodel and its slight variations [16, 17] add an artificialdiffusion term in the conformation field equation to avoid the loss of strict positive definiteness (SPD) of the con-formation tensor and subsequently numerical breakdowncaused by the hyperbolic nature of the FENE-P model(see Eq. (5)).Jin and Collins [18] stress the fact that much finer gridresolutions are required to fully resolve the polymer fieldthan the velocity and pressure fields. Indeed, the hyper-bolicity of the FENE-P model admits near discontinuitiesin the conformation and polymer stress fields [19]. Qual-itatively similar problems occur with shock waves andtheir full resolution in gas dynamic compressible flows,which is not practical using finer grids. In this case,high resolution numerical schemes such as slope-limiterand Godunov-type methods [20] have proved successfulat capturing the shock waves by accurately reproducingthe Rankine-Hugoniot conditions across the discontinuityto ensure the correct propagation speed.Motivated by these schemes, Vaithianathan et al. [21]adapted the second-order hyperbolic solver by Kurganovand Tadmor [22], which guarantees that a positive scalarremains positive over all space, to satisfy the SPD prop-erty of the conformation tensor and therefore avoid lossof evolution. Vaithianathan et al. [21] further demon-strated that this scheme dissipates less elastic energythan methods based on artificial diffusion, resulting instronger polymer-turbulence interactions. Moreover, ac-cording to the most recent review on the subject [2], thereare a lot of divergent and misleading results because ofthis artificial term introduced in the governing equations.For these reasons a modification of this shock-capturingscheme was developed in this present study to complywith non-periodic boundary conditions (see appendix 1).
A. Governing equations
The dimensionless incompressible Navier-Stokes equa-tions for a viscoelastic fluid take the form ∇ · u = 0 ∂ t u + ( u · ∇ ) u = − ∇ p + β Re c ∆ u + ∇ · σ (3)where β ≡ µ s /µ is the ratio of the solvent viscosity µ s to the total zero-shear-rate viscosity of the solution µ ,Re c ≡ U c δ/ν is the Reynolds number based on U c ≡ U b with U b the bulk velocity of the flow kept constant in timeand the channel’s half-width δ . The extra force in Eq.(3) arises due to polymers and the polymer stress tensorfor the FENE-P dumbbells is defined by the Kramersexpression σ = 1 − β Re c We c ( f ( tr C ) C − I ) (4)where We c ≡ τ p U c /δ with τ p the polymer relaxation timescale, f ( tr C ) ≡ L p − L p − tr C is the Peterlin function [23] and C ≡ h QQ i is the conformation tensor, which is definedas the dyadic product of the end-to-end vector Q of adumbbell that specifies its configuration. The Peterlinfunction prevents the dumbbell to reach its maximumextensibility, i.e. tr C ≤ L p , since as tr C → L p the forcerequired for further extension approached infinity. Notethat C and L p are made dimensionless by the equilib-rium length scale p k B T /H , where k B is the Boltzmannconstant, T is the solution temperature and H is theHookean spring constant and they have been normalisedsuch that the equibrium condition is C eq = I . Then, theconformation tensor is governed by the FENE-P model ∂ t C +( u · ∇ ) C − C · ∇ u − ∇ u ⊤ · C = − c ( f ( tr C ) C − I )(5)where the left hand side is the material derivative fora tensor field preserving its Galilean invariance and theright hand side represents deviation from the isotropicequilibrium due to Warner’s finite extensible non-linearelastic spring-force law [24].The elastic potential energy per unit volume E p storedby FENE-P dumbbells can now be specified by taking theintegral of Warner’s spring-force law over the end-to-endvector and after some algebra we obtain E p = 12 (1 − β )Re c We c ( L p −
3) ln( f ( tr C )) + E p (6)where E p is a constant reference energy at equilibrium.After that, taking the time derivative of the elastic po-tential energy ∂ t E p = 12 (1 − β )Re c We c ( L p −
3) 1 f ∂f∂C ii ∂C ii ∂t = 12 (1 − β )Re c We c f ∂C ii ∂t , (7)using the trace of Eq. (5), viz. ∂C ii ∂t = 2 C ik ∂ k u i − c ( f ( C kk ) C ii − δ ii ) (8)and similarly for the ∇ E p , we can derive the follow-ing balance equation for the elastic potential energy ofFENE-P dumbbells ∂ t E p + u · ∇ E p = σ · ∇ u − c f ( tr C ) tr σ (9)where E p is produced by σ · ∇ u , dissipated by c f ( tr C ) tr σ and transported by u · ∇ E p . B. Numerical parameters and procedures
Incompressible viscoelastic turbulence in a channel wassimulated in a rectangular geometry by numerically solv-ing the non-dimensional Eqs. (3)-(5) in Cartesian co-ordinates. After obtaining the new update of the con-formation tensor from the FENE-P model using themethod described in appendix 1, Eqs. (3)-(5) are numer-ically integrated with a fractional step method using a second-order Adams-Bashworth/Trapezoidal scheme (seeappendix 2). The fractional step method projects thevelocity field to a divergence free velocity field and thePoisson pressure equation is solved in Fourier space witha staggered grid for the pressure field [25]. The staggeredgrid for the pressure was used for numerical stabilitypurposes as well as the skew-symmetric implementationof the non-linear term in Eqs. (3). Spatial derivativesare estimated using sixth-order compact finite-differenceschemes [26]. The grid stretching technique used in theinhomogeneous wall-normal direction maps an equallyspaced co-ordinate in the computational space to a non-equally spaced co-ordinate in the physical space, in or-der to be able to use Fourier transforms [25, 27]. Fur-ther details of our numerical method are provided in[13]. Moreover, a validation of the algorithm just for theNavier-Stokes equations for turbulent channel flow can befound in [25], where this methodology was compared withspectral and second-order finite-difference schemes show-ing the necessity of spectral-like accuracy of the compacthigh-order schemes against second-order finite-differencesin turbulence computations.To simulate incompressible channel flow turbulence weapplied periodic boundary conditions for u ≡ ( u, v, w ) inthe x and z homogeneous directions and no-slip boundaryconditions u = 0 at the walls. The mean flow is in the x direction, i.e. h u i = ( h u ( y ) i , , h i in this paperdenotes averages in x , z spatial directions and time. Thebulk velocity U b in the x direction was kept constantfor all computations at all times by adjusting the meanpressure gradient − d h p i / d x at each time step. The choiceof U b in the computations for the Newtonian fluid is madebased on Dean’s formula Re τ ≃ . / c [28, 29] fora required Re τ ≡ u τ δν , where u τ is the friction velocityfor Newtonian fluid flow, i.e. β = 1 (see N cases in TableI).The procedure used for the computation of the vis-coelastic turbulent channel flows of Table I is the follow-ing. First, DNS of the Newtonian fluid, i.e. β = 1, wereperformed for the various Reynolds numbers until theyreached a steady state. Then, the initial conditions forthe viscoelastic DNS were these turbulent Newtonian ve-locity fields as well as the stationary analytical solution ofthe FENE-P model, given a steady unidirectional shearflow u = ( U ( y ) , , C ij tensor components C = 1 f ( C kk ) c f ( C kk ) (cid:18) d U d y (cid:19) ! C = We c f ( C kk ) d U d yC = C = 0 C = C = 1 f ( C kk ) f ( C kk ) = 23 cosh φ φ = cosh − (cid:0) Ω + 1 (cid:1) , Ω = √ c L p d U d y and dd y U = − y − assuming that U ( y ) = 0 . − ( y − ) ∀ y ∈ [0 ,
2] is a close approximation to the averaged velocityprofile of a Newtonian fully developed turbulent chan-nel flow at moderate Reynolds numbers [30]. Initially,the governing equations were integrated uncoupled, i.e. β = 1, until the conformation tensor achieved a station-ary state. From then on the fully coupled system of equa-tions, i.e. β = 1, was marched far in time, while u and C statistics were monitored for several successive timeintegrals until a fully developed steady state is reached,which satisfies the total shear stress balance across thechannel, viz. β Re c d h u i d y − h u ′ v ′ i + h σ i = u τ (cid:16) − yδ (cid:17) (11)where h σ i = − β Re c We c h L p − L p − C kk C i is the mean poly- mer shear stress. Finally, after reaching a statisticallysteady state, statistics were collected for several decadesof through-flow time scales L x /U b . In addition, existingturbulent velocity and conformation tensor fields wererestarted for computations where We c or L p was modi-fied. In these cases, the flow undergoes a transient time,where again sufficient statistics were collected after reach-ing a stationary state.According to Eqs. (3)-(5), the four dimensionlessgroups that can fully characterise the velocity and theconformation tensor fields are We c , L p , β and Re c , andthey are tabulated below. ?he reasons behind the choiceof the particular parameter values is outlined below. Therationale here follows the thorough parametric study byLi et al. [31]. TABLE I: Parameters for the DNS of viscoelastic turbulent channel flow. The friction Weissenberg number is defined byWe τ ≡ τ p u τ /ν . LDR cases: A, B, D2, I, J; HDR cases: C, D, D1, E, F, G, K; MDR case: H. Case We c We τ L p β Re c Re τ L x × L y × L z N x × N y × N z DR ( % ) N1 - - - 1 2750 123.8 6 . πδ × δ × . πδ × ×
100 0N2 - - - 1 4250 181 4 . πδ × δ × πδ × ×
100 0N3 - - - 1 10400 392.6 2 πδ × δ × . πδ × ×
100 0A 2 15.4 120 0.9 4250 167.7 4 . πδ × δ × πδ × ×
100 -14.2B 4 30.8 120 0.9 4250 147.3 4 . πδ × δ × πδ × ×
100 -33.8C 7 54 120 0.9 4250 121.8 4 . πδ × δ × πδ × ×
100 -54.7D 9 69.4 120 0.9 4250 118.3 4 . πδ × δ × πδ × ×
100 -57.3D1 9 69.4 60 0.9 4250 124.7 4 . πδ × δ × πδ × ×
100 -52.5D2 9 69.4 30 0.9 4250 150.3 4 . πδ × δ × πδ × ×
100 -31E 11 84.8 120 0.9 4250 113.3 4 . πδ × δ × πδ × ×
100 -60.8F 13 100.2 120 0.9 4250 112.4 4 . πδ × δ × πδ × ×
100 -61.4G 15 115.6 120 0.9 4250 111.4 4 . πδ × δ × πδ × ×
100 -62.1H 17 131 120 0.9 4250 107.8 8 πδ × δ × πδ × ×
100 -64.5I 2 29.6 120 0.9 10400 323.3 2 πδ × δ × . πδ × ×
100 -32.2J 4 22.3 120 0.9 2750 106.9 6 . πδ × δ × . πδ × ×
100 -25.4K 7 39 120 0.9 2750 91.1 6 . πδ × δ × . πδ × ×
100 -45.9
Drag reduction effects are expected to be stronger athigh Weissenberg numbers. In fact, higher levels of per-centage drag reduction at MDR have also been measuredfor higher Reynolds numbers [7], showing the Reynoldsnumber dependence on drag reduction amplitude. There-fore, in this work, an extensive parametric study has beencarried out by mainly varying We c for the computation-ally affordable Re c = 4250 to determine the impact ofpolymer dynamics on the extent of drag reduction. TheReynolds numbers considered here, Re c = 2750 , c and Re c . The chosen L p = b + 3 valuesare representative of real polymer molecule lengths whichcan be related through b ≈ N C / ( σ sf N ) where N C is thenumber of carbon atoms in the backbone of the polymermacromolecule, σ sf is an empirical steric factor, N is thenumber of monomers and b is the finite dumbbell exten-sibility and is a large number [14]. Note that in the limit b → ∞ , the Hookean spring-force law is recovered, whichgoverns a linear spring.Low β values were used in most prior DNS to achievehigh levels of drag reduction, in view of the attenuation ofthe polymer-turbulence interactions due to the additionalartificial diffusion term in the FENE-P model and theirmoderate Reynolds numbers, usually Re τ ≡ δ/δ ν ≤ β = 0 . β values may lead tosignificant shear-thinning [11] unlike in experiments ofpolymer drag reduction. The fact that the numericalscheme applied in our study for the FENE-P model isexpected to capture the strong polymer-turbulence in-teractions allows the value of β , which is inversely pro-portional to the polymer concentration, to be high, i.e. β = 0 .
9, representative of dilute polymer solutions usedin experiments.The box sizes L x × L y × L z , where subscripts indi-cate the three Cartesian co-ordinates, were chosen withreference to the systematic study by Li et al. [31] ofhow the domain size influences the numerical accuracy.Specifically, they point out that long boxes are requiredin DNS of polymer drag reduction, particularly in thestreamwise direction because of longer streamwise corre-lations at higher percentage DR, as opposed to the min-imal flow unit [33] used in many earlier works. Differ-ent grid resolutions N x × N y × N z were tested for con-vergence. In particular, the following set of resolutions128 × ×
64, 200 × ×
100 and 256 × ×
128 were triedfor Re τ ≃
180 with the two latter giving identical meanvelocity profiles and not significantly different rms veloc-ity and vorticity profiles. Similarly, grid sensitivity testswere carried out for the other Re τ cases. Eventually, theresolutions for each Newtonian fluid computation werevalidated against previously published databases for thecorresponding Re τ cases [34–36]. Note that if the res-olutions for Newtonian turbulent computations are ad-equately resolving the flow scales, then the same reso-lutions are sufficient for viscoelastic turbulent computa-tions, since the size of vortex filaments in these flowsincreases while their number decreases as drag reduces[2].For a given resolution, viscoelastic computations re-quire approximately 4 times more memory and 2 timesmore CPU time per time step compared to the Newto-nian case. The time step ∆ t used in viscoelastic compu-tations is typically a factor of 5 smaller than that usedin the Newtonian cases due to the stricter CFL conditionof the present numerical method for the FENE-P model(see Eq. (A.9) in the appendix and [26] for more detailson the time step constraint using compact schemes). Ul-timately, the viscoelastic computations require approxi-mately 10 times more CPU resources than the Newtoniancomputations for a given computational time period. the shear stress increases slower than linear σ ∝ S III. VISCOELASTIC TURBULENCESTATISTICSA. Polymer drag reduction
Since the computations are performed with a con-stant flow rate by adjusting the pressure gradient, DRis manifested via a decrease in skin friction, i.e. lowerRe τ = δ/δ ν values as drag reduces. Here, we define per-centage drag reduction as a negative quantityDR ≡ − d h p i d x − (cid:16) − d h p i d x (cid:17) (cid:12)(cid:12) − d h p i d x (cid:12)(cid:12) × u τ − u τ (cid:12)(cid:12) u τ (cid:12)(cid:12) × (cid:18) Re τ Re τ (cid:19) − ! × u τ = − δρ d h p i d x . Variables with and without subscript0 in Eq. (12) refer to Newtonian and viscoelastic fluidflow, respectively. Note that for a direct comparison be-tween the various cases with different skin friction wechoose our plots to be presented in terms of y/δ insteadof y/δ ν because δ is the same for all our computations,whereas δ ν increases with drag reduction.Figure 1 depicts the capability of the current numer-ical scheme for the FENE-P model to enable strongerpolymer-turbulence interactions than artificial diffusionmethods. Higher values of percentage drag reduction asfunction of Weissenberg number are obtained comparingwith earlier DNS studies (see Fig. 1b in [37]) withoutthe need for low β values [32]. These DR values extend We c = 1 MDR limit A B C D D1 D2 E F G H DR ( % ) We τ o FIG. 1: (Colour online) Variation of percentage drag reduc-tion with Weissenberg number. Any departure from the Newtonian behaviour, i.e. σ ij ∝ S ij ,with some constant of proportionality independent of the rate ofstrain, could be called non-Newtonian. throughout the drag reduction regimes. The MDR limitis approached in this case at | DR | ≃
65% because ofthe moderate Re c in our computations. Even so, thisamount of drag reduction falls within the MDR regime,based on the classification of drag reduction by Warholic et al. [5], allowing to study the MDR dynamics of thepolymer molecules and their effects on the flow in thisasymptotic state. B. Effects of polymer extensibility and Reynoldsnumber
The effects of maximum dumbbell extensibility isbriefly considered for three different extensibilities L p =30 ,
60 and 120 but the same We c and Re c (see D casesin Table I). Figure 1 shows that the extent of drag re-duction is amplified by longer polymer chains consistentwith other DNS studies [31, 38]. This effect is relatedto the fact that the average actual length of the dumb-bells, represented by the trace of the conformation tensor h C kk i , increases further for larger L p according to Fig.2a, inducing stronger influence of the polymers on theflow. The percentage increase, however, of the polymersextension is less for larger FENE-P dumbbells (see Fig.2b), suggesting that large polymer molecules could be lesssusceptible to chain scission degradation, which causesloss of drag reduction in experiments [2]. The near-wallturbulence dynamics play an important role for all threecases, as most of the stretching happens near the wall,where the highest fluctuating strain rates are expected.Eventually, the largest maximum length, i.e. L p = 120,was used for the rest of the computations considered inthis work in order to explore the polymer dynamics ateffective drag reductions, which are interesting not onlyfundamentally but also in many real life applications.Based on DNS with artificial diffusion methodology,Housiadas and Beris [39] claim that the extent of dragreduction is rather insensitive to Reynolds numbers inthe range between 125 ≤ Re τ ≤
590 for LDR flows.On the other hand, avoiding the use of artificial diffusionin our study, the Reynolds number dependence on dragreduction for cases with identical We c values but differentReynolds numbers, i.e. Re c = 2750 , c at all instances (see Table I). This Reynolds numberdependence is further depicted in the polymer dynamicsthrough the profiles of h C kk i /L p , which amplify closer tothe wall due to more intense strain rates in this region athigher Re c and collapse towards the centre of the channel(see Fig. 3). The disparate behaviour of h C kk i withrespect to y/δ due to the Reynolds number dependenceis anticipated by the broader spectra of flow time scalesthat are encountered at higher Re c by the dumbbells withfixed relaxation time scale. The fact that the currentDNS could capture the Reynolds number dependence on y/ δ < C kk > Case DCase D1Case D2 y/ δ < C kk > / L p FIG. 2: (Colour online) Effect of maximum dumbbell length.Plots of (a) average actual dumbbell extensibility h C kk i and(b) percentage average dumbbell extensibility h C kk i /L p asfunctions of y/δ . Note: case D ( L p = 120); case D1 ( L p = 60);case D2 ( L p = 30). drag reduction and polymer dynamics, emphasises oncemore the strong polymer-turbulence interactions that canbe captured by the present numerical approach even atlow levels of drag reduction.It is essential to note at this point that the interme-diate dynamics between the von K´arm´an and the MDRlaw, i.e. the LDR and HDR regimes, are non-universalbecause they depend on polymer concentration, chemicalcharacteristics of polymers, Reynolds number, etc. [7, 9].Here, this is illustrated by the maximum dumbbell lengthand Reynolds number dependencies of the polymer dy-namics in Figs. 2 and 3, respectively. However, at theMDR limit, which is achieved at We c ≫ c ≫ y/ δ < C kk > / L p Case ACase BCase CCase ICase JCase K
FIG. 3: (Colour online) Effect of Reynolds number on per-centage average dumbbell extensibility as function of y/δ .Identical symbols correspond to cases with the same We c val-ues. Note: Compare case A (We c = 2, Re c = 4250) with caseI (We c = 2, Re c = 10400); case B (We c = 4, Re c = 4250) withcase J (We c = 4, Re c = 2750); case C (We c = 7, Re c = 4250)with case K (We c = 7, Re c = 2750). C. Mean and fluctuating velocity statistics
The picture of drag reduction can be analysed in fur-ther detail with the statistics of the turbulent velocityfield introduced in Figs. 4 and 5. The distinct differ-ences in the statistical trends of the turbulent velocityfield between the LDR and HDR regime, that have beenobserved experimentally, are clearly identified in theseresults. For clarity, a few indicative cases from the dataof Table I have been chosen for plotting, representing theLDR, HDR and MDR regimes for different Weissenbergnumbers at Re c = 4250.According to Fig. 4 and noting that β = 0 . y + .
10 to the linear variation U + = β − y + , which can be deduced for viscoelastic flows, byrewriting Eq. (11) in viscous scales β d U + d y + − h u ′ v ′ i u τ + h σ i u τ = 1 − y + Re τ (13)and neglecting the normalised Reynolds and mean poly-mer shear stress in the viscous sublayer y + → c values at the same Re c . The profileof the Newtonian case N2 is in agreement with the vonK´arm´an law Eq. (1), which does not hold for viscoelasticturbulent flows. Specifically, the curves of cases A andB (LDR regime) are shifted upwards with higher valuesof the intercept constant B , i.e. parallel to the profile y + U + Case N2Case ACase BCase DCase GCase H
FIG. 4: (Colour online) Mean velocity profiles versus y + forthe LDR, HDR and MDR regimes. – · –: U + = y + , - - -: U + = . log y + + 6 . · · · : U + = . log y + −
17. Note: case N2(DR = 0%); case A (DR = − . − . − . − . − . of the Newtonian flow (see Fig. 4), increasing DR. Thispicture is consistent with the phenomenological descrip-tion by Lumley [8, 40], where the upward shift of theinertial sublayer can be interpreted as a thickening ofthe buffer or elastic layer for viscoelastic flows, which isequivalent to drag reduction. HDR cases D and G exhibitdifferent statistical behaviour than LDR flows with theslope of the log-region increasing until the MDR asymp-tote is reached by case H. Overall, the same behaviouracross the extent of drag reduction in viscoelastic tur-bulent flows have been seen in several experimental andnumerical results [2].Different statistical trends between low and high dragreduction have also been observed experimentally [5, 41]for the rms streamwise velocity fluctuations normalisedwith u τ . Figure 5a illustrates the growth of the peakin the profile of u ′ + for LDR case A and B at lowWe c and a notable decrease for the rest of the casesat HDR/MDR with high We c values. The peaks movemonotonically away from the wall throughout the dragreduction regimes indicating the thickening of the elas-tic layer, which is compatible with the behaviour of themean velocity profile.Note that this is the first time that a DNS computationcan so distinctly attain this behaviour. This is attributedto the accurate shock-capturing numerical scheme we ap-plied for the FENE-P model in this study. It has tobe mentioned however that there have been three earlierstudies [32, 37, 42], which use the artificial diffusion algo-rithms for FENE-P and showed similar but not as cleartrends for u ′ + in a DNS of viscoelastic turbulent channelflow. In fact, Min et al. [37] reached the HDR/MDRregime at roughly | DR | ≃ et al. [32] had to y/ δ u , + Case N2Case ACase BCase DCase GCase H y/ δ v , + y/ δ w , + FIG. 5: (Colour online) Rms velocity components for theLDR, HDR and MDR regimes. (a) Streamwise u ′ + , (b) wall-normal v ′ + and (c) spanwise w ′ + profiles versus y/δ . Note:case N2 (DR = 0%); case A (DR = − . − . − . − . − . use β = 0 . β values, have not been able to ob-tain this transition effect on the statistics of u ′ + betweenthe drag reduction regimes.Finally, the wall-normal v ′ + and spanwise w ′ + rms ve-locity fluctuations in Figs. 5b and 5c, respectively, arecontinuously attenuated while DR is enhanced by in-creasing the polymer relaxation time scale. Again, themonotonic displacement of their peaks towards the cen-tre of the channel as drag reduction amplifies is consistentwith that of the mean velocity profile and with experi-mental and other numerical studies [2]. D. Fluctuating vorticity statistics
The rms statistics of the fluctuating vorticity field nor-malised by viscous scales, i.e. ω ′ + ≡ ω ′ δ ν /u τ , are pre-sented in Fig. 6 for representative cases from Table I atvarious levels of drag reduction. The streamwise vortic-ity fluctuations ω ′ x + demonstrate a persistent attenuationalong the normalised distance y/δ as drag reduction en-hances due to the increase of We c (see Fig. 6a). In thenear-wall region y/δ < . ω ′ x + provides evidence for a dragreduction mechanism based on the suppression of thenear-wall counter-rotating steamwise vortices [46, 47],which underpin considerable amount of the turbulenceproduction [48].The wall-normal rms vorticity is zero at the wall dueto the no-slip boundary condition and reaches its peakwithin the buffer layer (see Fig. 6b). The intensity of ω ′ y + is reduced for all levels of drag reduction accordingto Fig. 6b, with the position of the near-wall peaks mov-ing towards the centre of the channel as We c becomeslarger, representing once more the thickening of the elas-tic layer in a consistent way. Most of the inhibition of ω ′ y + happens near the wall and slightly towards the cen-tre of the channel only for the HDR/MDR cases G andH, i.e. for | DR | > ω ′ z + ,where the spanwise vorticity fluctuations decrease inthe near-wall region y/δ . . u ′ + between the LDRand HDR/MDR regimes (see Fig. 5a) plus the contin-uous drop of v ′ + (see Fig. 5b) in viscoelastic drag re- y/ δ ω x , + Case N2Case ACase BCase DCase GCase H y/ δ ω y , + y/ δ ω z , + FIG. 6: (Colour online) Rms vorticity profiles for the LDR,HDR and MDR regimes. (a) Streamwise ω ′ x + , (b) wall-normal ω ′ y + and (c) spanwise ω ′ z + profiles versus y/δ . Note: case N2(DR = 0%); case A (DR = − . − . − . − . − . duced flows. As a final note, ω ′ z + > ω ′ x + > ω ′ y + inthe viscous sublayer, i.e. y/δ < .
05 for all cases and ω ′ z + ≃ ω ′ x + ≃ ω ′ y + in the inertial and outer layer for theNewtonian case N2. However, ω ′ z + > ω ′ y + > ω ′ x + awayfrom the wall when drag reduces for viscoelastic flows,which manifests the dominance of small scale anisotropyin the inertial and outer layer at HDR and MDR. IV. CONFORMATION AND POLYMER STRESSTENSOR
Before looking at the mean momentum and energybalance, the study of the conformation tensor field isessential to get an understanding of the polymer dy-namics in support of the results provided by this newnumerical method for the FENE-P model in turbulentchannel flow. The symmetries in the flow geometry de-termine properties of tensor components in the averagesense [49]. In the current DNS of turbulent channelflow, statistics are independent of the z direction andthe flow is also statistically invariant under reflectionsof the z co-ordinate axis. Therefore, for the probabilitydensity function f ( Q ; x , t ) of a vector Q , these two con-ditions imply ∂f /∂z = 0 and f ( Q , Q , Q ; x, y, z, t ) = f ( Q , Q , − Q ; x, y, − z, t ). Then, at z = 0 reflectionalsymmetry suggests that h Q i = −h Q i ⇒ h Q i = 0 andsimilarly for h Q Q i = h Q Q i = 0. So, in this case themean conformation tensor reduces to h C ij i = h C i h C i h C i h C i
00 0 h C i (14)where the non-zero components scaled with L p arepresented in Figs. 7 and 8 with respect to y/δ for cases at various drag reduction regimes (see Ta-ble I). The zero components in our study havebeen found to be zero within the precision accuracy.Turbulent channel flow is also statistically symmet-ric about the plane y = δ . Therefore, this re-flectional symmetry imposes f ( Q , Q , Q ; x, y, z, t ) = f ( Q , − Q , Q ; x, − y, z, t ), which implies that the nor-mal components of h C ij i are even functions and h C i isan odd function comparable to the Reynolds stress tensorcomponents.The normalised trace of the mean conformation tensor h C kk i /L p is plotted in Fig. 7a together with h C i /L p .Notice that the dominant contribution in the trace comesfrom h C i , i.e. h C kk i ≃ h C i at all Weissenberg num-bers, reflecting on average a strong preferential orien-tation of the stretched dumbbells along the streamwisedirection. The fact that h C i ≫ h C i ≃ h C i > h C i denotes the strong anisotropic behaviour of the meanconformation tensor caused by the mean shear in turbu-lent channel flow. This anisotropy influences the statis-tics of the fluctuating velocity field particularly at smallscales, as was mentioned in section III D. The curves of h C i /L p and consequently of h C kk i /L p constantly risewith most of the stretching happening close to the walland growing towards the centre of the channel, since0 y/ δ < C kk > / L p , < C > / L p Case ACase BCase DCase GCase H −3 y/ δ < C > / L p FIG. 7: (Colour online) Profiles of (a) h C kk i /L p (line-symbols) and h C i /L p (solid lines), (b) h C i /L p as functionsof y/δ for the LDR, HDR and MDR regimes. Note: case A(DR = − . − . − . − . − . higher values of polymer time scale are influenced from awider spectrum of flow time scales. A local minimum anda maximum emerge in the near-wall region y/δ < .
2, in-duced by the streamwise vortices [42, 50]. These peaksmove apart from each other and away from the wall forhigher We c values. Figure 7a also shows that the ampli-tudes of these peaks seem inversely proportional to thepeak amplitudes of ω ′ x + as drag reduces (see also Fig.6a).Moreover, as We c increases the profiles of h C i /L p and h C i /L p amplify, reaching their peaks at not muchdifferent y/δ for each We c case (see Figs. 7b and 8b).In particular, the values of h C i /L p at the wall are de-pendent on the polymer relaxation time scale unlike for h C i /L p . On the other hand, the values of h C i /L p depend on Weissenberg number at y = δ in contrast to h C i /L p , which is zero for all cases because of the sym-metry mentioned earlier. The behaviour of h C i /L p in −3 y/ δ < C > / L p Case ACase BCase DCase GCase H −3 y/ δ < C > / L p FIG. 8: (Colour online) Profiles of (a) h C i /L p and (b) h C i /L p as functions of y/δ for the LDR, HDR and MDRregimes. Note: case A (DR = − . − . − . − . − . Fig. 8a is more peculiar with respect to We c , with theprofiles increasing within the LDR regime and attenuatefor HDR and MDR cases, in a similar manner to u ′ + (seeFig. 5a). Its peak values are achieved closer to the coreof the channel in comparison to the rest of the confor-mation tensor components. This points out the differentflow time scales that are important for h C i , exemplify-ing the complex dynamics of the polymers, even in thissimple mechanical model.It is interesting to mention that the components of h C ij i have different asymptotic rates of convergence to-wards the limit We c → ∞ . It is known that for We c ≫ h C kk i ≤ L p and subse-quently in our case h C i . L p (see Fig. 7a), where thisupper bound is far from achieved in our computations.This result demonstrates that highly stretched polymersare not required for the manifestation of drag reduc-tion or even of the MDR asymptote, as de Gennes [51]claims against Lumley’s [8] assumption of a coil-stretch1transition, i.e. highly stretched polymer molecules, forthe enhancement of intrinsic viscosity. The components h C i /L p and h C i /L p seem to have almost reached theirasymptotic limit with the MDR case H according to Figs.7b and 8b, respectively. Finally, h C i /L p has not yetconverged to its limit, decreasing with a slow rate to-wards very small values for high We c . In fact, it hasbeen argued theoretically that h C i → et al. [53] consider the FENE-P model integrated overthe x , z spatial directions and time, assuming statisticalstationarity and homogeneity in x and z h u ∂ x C ij i = h C ik ∂ x k u j i + h C jk ∂ x k u i i− c h f ( C kk ) C ij − δ ij i . (15)Then, taking the Reynolds decomposition of the velocityfield u i = h u i i + u ′ i , one obtains1We c h f ( C kk ) C ij − δ ij i = h C ik i ∂ x k h u j i + h C jk i ∂ x k h u i i + Q ij (16)where Q ij = h C ik ∂ x k u ′ j i + h C jk ∂ x k u ′ i i − h u ′ ∂ x C ij i .Therefore, the average polymer stress tensor defined byEqs. (4) takes the form h σ ij i == 1 − β Re c h C i ∂ h u i ∂x + Q h C i ∂ h u i ∂x + Q Q h C i ∂ h u i ∂x + Q Q Q Q Q Q . (17)Now, the important assumption at the limit of a localWeissenberg number We S ≡ τ p dd y h u i → ∞ is that Q and Q can be neglected, considering the polymers tobe stiff, i.e. C ij → h C ij i , mostly in the main stretchingdirections and the correlations between fluctuating con-formation tensor and velocity fields in the other Cartesiandirections to remain minimal at this limit. In this case,as a result h σ i = A − β Re c h C i ∂ x h u i (18) h σ i = A − β Re c h C i ∂ x h u i (19)where A and A are expected to be independent of y and equal to 1 at some intermediate region in the flow asWe S ≫
1. This hypothesis is checked in Fig. 9 againstvarious viscoelastic DNS from Table I.Figure 9a shows clearly that A tends to a constantand reaches 1 in the region 0 . y/δ . . S values justifying that Q can be neglected for HDR and y/ δ A ≡ < f C − > / ( W e c < C > d < u > / d y ) y/ δ A ≡ < f C > / ( W e c < C > d < u > / d y ) Case ACase BCase DCase GCase H
FIG. 9: (Colour online) Scalings of the compensated polymerstress components (a) A ≡ h σ i / (cid:16) − β Re c h C i d h u i d y (cid:17) and (b) A ≡ h σ i / (cid:16) − β Re c h C i d h u i d y (cid:17) with respect to y/δ . Note: caseA (DR = − . − . − . − . − . MDR cases. Note that A deviates from 1 towards thecentre of the channel because We S becomes small in thisregion. A is approximately independent of y in someintermediate region for almost all cases and appears totend towards 1 as We S increases (see Fig. 9b). However,the polymer relaxation time scales used in this study arenot sufficiently high for A →
1. So, for our DNS re-sults the polymer shear stress can be considered to be h σ i ∝ − β Re c h C i ∂ x h u i in a range 0 . . y/δ . . h C i is the component in-volved in the MDR dynamics, bearing in mind that h C i ≫ h C i ≃ h C i > h C i . In the end, both Figs.7b and 9a confirm the claims that h C i has reached itsasymptotic limit within the Weissenberg numbers con-sidered at this particular Reynolds number in this study,unlike h C i (see Figs. 8a and 9b).2 V. SHEAR STRESS BALANCE
The balance of total shear stress Eq. (11) is con-sidered in this section. The total shear stress in vis-coelastic turbulent channel flow contains the shear stress βν dd y h u i coming from the mean flow, the Reynolds stress −h u ′ v ′ i rising from turbulence and the mean polymershear stress h σ i due to polymers in the flow, whichis also referred to as the Reynolds stress deficit since ν dd y h u i − h u ′ v ′ i 6 = u τ (1 − y/δ ). In addition, the gradientof the Reynolds shear stress rises due to the non-linearterm in the Navier-Stokes equations and it is related tothe lift force (Magnus effect) experienced by a vortex lineexposed to a velocity u [54]. The Reynolds shear stressand the rotational form u × ω of the non-linear termin Eq. (3) are related through the following equation ∂∂y h u ′ v ′ i = h v ′ ω ′ z i − h w ′ ω ′ y i . It is true that tangles ofvery intense and slender vortex filaments can exist downto very fine scales creating an energy drain on the meanflow. These tangles can be seen as one form of inter-mittency in turbulent flows [55]. However, it is unclearhow the intermittency of the vorticify field ω affects theaverages h v ′ ω ′ z i and h w ′ ω ′ y i and thereby Reynolds shearstress via the relation ∂∂y h u ′ v ′ i = h v ′ ω ′ z i − h w ′ ω ′ y i .The viscous stress of the solvent, the Reynolds shearstress and the mean polymer shear stress normalised withviscous scales are presented in Fig. 10 at different levelsof percentage DR for cases with the same Re c from TableI. At the wall, the no-slip boundary condition enforces −h u ′ v ′ i (cid:12)(cid:12) y =0 = 0. The wall shear stress is governed by90% viscous as well as 10% polymer contribution for allviscoelastic cases as opposed to the Newtonian case N2.Viscosity is the dominant parameter in the near-wall re-gion but becomes more influential in the outer regionsas drag reduction enhances. This is clear from Fig. 10awhere β dd y + U + increases monotonically towards the cen-tre of the channel as We c increases. Viscoelastic effectsbecome also more significant towards the centre of thechannel for higher We c cases. Reynolds shear stress, onthe other hand, is constantly decorrelated at higher per-centage DR with its peak shifting away from the wall.It is interesting to see that for the HDR case D −h u ′ v ′ i and h σ i are comparable and as MDR is approached thepolymer shear stress plays an increasingly fundamentalrole in sustaining turbulence due to the vast attenua-tion of the Reynolds shear stress at these finite Reynoldsnumber computations. This becomes apparent in thenext section by analysing the turbulent kinetic energybudget.Notice that Reynolds shear stress remains finite atMDR confirming the experimental measurements byPtasinski et al. [41] against the complete depletion of −h u ′ v ′ i reported by Warholic et al. [5] and their sub-sequent claim that turbulence is sustained entirely bypolymer stresses. What can be said theoretically on thiscontroversy is the following. Consider first the limit ofWe S → ∞ , where A → y/ δ β d U + / d y + Case N2Case ACase BCase DCase GCase H y/ δ −< u ’ v ’ > / u τ y/ δ < σ > / u τ FIG. 10: (Colour online) Profiles of (a) viscous shear stress,(b) Reynolds shear stress and (c) mean polymer shear stressversus y/δ for the LDR, HDR and MDR regimes. Note: caseN2 (DR = 0%); case A (DR = − . − . − . − . − . Reynolds numbers, as Fig. 9b suggested. Then, the to-tal shear stress balance Eq. (11) can be rewritten using3Eq. (19) ν ( β + (1 − β ) h C i ) d h u i d y − h u ′ v ′ i ≃ u τ (cid:16) − yδ (cid:17) (20)where ν eff ( y ) ≡ ν ( β + (1 − β ) h C i ) is an effective vis-cosity similar to the one encountered in Lumley’s phe-nomenology [9, 40]. Now, when We S ≫ h C i becomes minimal based on theoretical claims by[9, 52] and observational indications in this study. Then,for high enough Reynolds number along the universalMDR asymptotic line, i.e. taking first the infinite Weis-senberg number limit and then the infinite Reynoldsnumber limit, one might expect an intermediate region δ ν ≪ y ≪ δ of approximately constant Reynolds shearstress, i.e. −h u ′ v ′ i /u τ →
1, implied by Eq. (20) whentaking the limits of y/δ → y/δ ν → ∞ with the rea-sonable assumption that νβ dd y h u i → y ≫ δ ν . Thisstatement suggests that the classical way of turbulenceproduction does not vanish in the infinite Weissenbergand Reynolds number limit.Ultimately, the conjecture here is that h σ i becomesnegligible in the stress balance Eq. (11) when carefullytaking the double Weissenberg and Reynolds numberlimit in the right order, so that we are on the universalMDR asymptotic line. This, however, does not indicatethat drag reduction is depleted, it rather suggests thatthe MDR asymptote could be entirely determined by theenergetics in these infinite limits. Nevertheless, polymersplay a crucial role in the dynamics at MDR and this willbe explored further in the next section. VI. POLYMER-TURBULENCE DYNAMICALINTERACTIONS
The balance equation for the turbulent kinetic energyof a viscoelastic fluid provides further insight into the dy-namical interactions between polymers and turbulence.Assuming statistical stationarity and homogeneity in x and z directions for the mean turbulent kinetic energybalance and integrating over the y direction, we obtain Z P d y = Z ε N d y + Z ε P d y (21)with no contribution from the transport terms due to theno-slip boundary condition, using the divergence theo-rem. The turbulence production by Reynolds shear stressis denoted here by P ≡ −h u ′ v ′ i dd y h u i , the viscous dissi-pation rate ε N ≡ νβ h s ij s ij i and the viscoelastic dissipa-tion rate ε P ≡ h σ ′ ij ∂ x j u ′ i i , which arises due to fluctuatingpolymer stresses. Note that ε P has a dual nature, i.e. itcan serve either as dissipation or production dependingon the signs of the polymer stress fluctuations and thatof the fluctuating velocity gradients.Figure 11 presents each term of Eq. (21) normalisedby δ ν /u τ with respect to We τ for all cases from Ta-ble I at Re c = 4250. An asymptotic behaviour to a N2 A B C D E F G H We τ o ∫ P δ ν / u τ d y , ∫ ε N δ ν / u τ d y , ∫ ε P δ ν / u τ d y ∫ P dy ∫ ε N dy ∫ ε P dy FIG. 11: (Colour online) Terms of the y -integrated turbulentenergy balance with respect to We τ . marginal flow state can be observed by increasing thepolymer relaxation time scale with a vast attenuation oc-curring in the total production and viscous dissipation,while viscoelastic dissipation grows mildly in the LDRregime and constantly decays within HDR and MDR.Overall, R ε P d y becomes pivotal in the dynamics of theflow relative to R P d y and R ε N d y for HDR and MDRflows. Most importantly R ε P d y < τ val-ues according to Fig. 11, in agreement with experimen-tal measurements [41], implying that polymers somehowcan sustain turbulence by producing turbulent kineticenergy. Notice, that in this plot both dissipations arepresented as positive quantities and this was done onpurpose to emphasise the interplay between productionand viscous dissipation from LDR to HDR. It is notewor-thy that R P d y > R ε N d y for LDR cases A and B but R P d y < R ε N d y for HDR cases and gets even smalleras drag reduction approaches its maximum limit. Thisobservation hints that polymer dynamics get somehowinvolved in the production of turbulent kinetic energy sothat turbulence does not die out at HDR and MDR.Lets now look in more detail at the profiles of P , ε N and ε P scaled by δ ν /u τ with respect to normalised dis-tance from the wall y/δ for representative cases at variouslevels of drag reduction from Table I (see Fig. 12). Dissi-pation represents drain of energy, hence, ε N and ε P havebeen plotted here as negative quantities.The production of turbulent energy by Reynoldsstresses, which is continuously reduced over the extend ofdrag reduction as a function of We c , serves to exchangekinetic energy between the mean flow and the turbulence.The local peak of P is reached within the buffer layerand in fact for Newtonian flows we can easily show thatthe maximum production occurs where −h u ′ v ′ i = ν dd y h u i and P max δ ν /u τ < [49]. The peak turbulence produc-tion within the LDR regime also occurs at the intersec-tion point of viscous and Reynolds shear stress (compare4 y/ δ P δ ν / u τ Case N2Case ACase BCase DCase GCase H y/ δ − ε N δ ν / u τ y/ δ − ε P δ ν / u τ FIG. 12: (Colour online) Profiles of (a) turbulence production P δ ν /u τ , (b) viscous dissipation − ε N δ ν /u τ and (c) viscoelasticdissipation − ε P δ ν /u τ for the LDR, HDR and MDR. Note:case N2 (DR = 0%); case A (DR = − . − . − . − . − . Figs. 10a and 10b with Fig. 12a), which shifts awayfrom the wall as We c increases, indicating the thickening of the elastic layer. However, for HDR and MDR cases P max δ ν /u τ is within 0 . < y/δ . .
3, where the maxi-mum Reynolds stress roughly appears, without followingthe −h u ′ v ′ i = βν dd y h u i intersection point, which does noteven exist for cases G and H (see Figs. 10a and 10b).Viscous dissipation exhibits monotonic attenuation asdrag reduces for higher values of We c with the maximumdissipation arising at the wall for the Newtonian case N2and the LDR cases A and B (see Fig. 12b). Althoughthe kinetic energy is zero at the wall since u ′ | y =0 = 0 im-posed by the no-slip boundary conditions, the fluctuatingstrain rate and consequently ε N is non-zero. At high per-centage DR, we surprisingly observe that the highest fluc-tuating strain rates are encountered away from the wallproviding a completely different picture of the near-walldissipation dynamics. The local kink in the buffer/elasticlayer, which arises due to intense activity in this region,exists at corresponding y/δ with P max δ ν /u τ for all casesconsidered in Fig. 12b and becomes a global minimumfor the HDR and MDR cases, dominating the profiles ofviscous dissipation.The profiles of viscoelastic dissipation obey a charac-teristic transitional trend similar to what has been al-ready observed for u ′ + (see Fig. 5a) and h C i (see Fig.8a) from LDR to HDR regime, as We c increases. Indetail, the curves of LDR cases A and B shift down-wards increasing viscoelastic dissipation but those of theHDR/MDR cases move upwards enhancing the positivenature of − ε P δ ν /u τ . The dual nature of ε P is clearly de-picted in Fig. 12c with polymers dissipating and produc-ing turbulent kinetic energy in different regions, whichdepend on the polymer relaxation time scale at a givenReynolds number. A Reynolds number dependence ofthese regions is expected owing to the effect of differentflow time scales on dumbbells with a particular relaxationtime scale. Figure 13 compares cases of identical Weis-senberg numbers and different Reynolds numbers (seeTable I), illustrating a weaker Re c dependence on vis-coelastic dissipation in comparison to the stronger We c dependence in Fig. 12c, particularly at HDR and MDR.The part of the total dissipation that occurs in the threeregions defined by the profile of viscoelastic dissipation inFig. 12c can be estimated based on the profiles in Figs.12b and 12c. Approximately 15% to 25% of the totaldissipation takes place in the first region, 25% to 60% inthe second region and 20% to 70% in the third region. Inother words, the majority of the total dissipation occursaway from the wall.Now, considering each component of the correlationmatrix ε P ≡ h σ ′ ij ∂ x j u ′ i i , where summation applies overthe indices i and j , we can observe that components with i = 2 , i = 1 components according to Fig. 14 whichis very similar to Fig. 12c. The qualitative featuresof ε P are clearly captured by h σ ′ j ∂ x j u ′ i , simplifying theunderpining dynamics of viscoelastic dissipation. How-ever, to be precise, ε P is neither exactly approximatenor proportional to h σ ′ j ∂ x j u ′ i . Note that the positive5 y/ δ − ε P δ ν / u τ Case ACase BCase CCase ICase JCase K
FIG. 13: (Colour online) Effect of Reynolds number on vis-coelastic dissipation as function of y/δ . Identical symbolscorrespond to cases with the same We c values. Compare caseA (We c = 2, Re c = 4250) with case I (We c = 2, Re c = 10400);case B (We c = 4, Re c = 4250) with case J (We c = 4,Re c = 2750); case C (We c = 7, Re c = 4250) with case K(We c = 7, Re c = 2750). y/ δ − < σ , j ∂ x j u , > δ ν / u τ Case ACase BCase DCase GCase H
FIG. 14: (Colour online) −h σ ′ j ∂ x j u ′ i δ ν /u τ as function of y/δ for the LDR, HDR and MDR regimes. Note: case A (DR = − . − . − . − . − . nature of ε P is caused by the correlations −h σ ′ ∂ x u ′ i and −h σ ′ ∂ x u ′ i (see Figs. 15a and 15c). The rest of thecomponents are negative for all cases considered here anddecrease monotonically as We c increases like −h σ ′ ∂ x u ′ i in Fig. 15b. The only exception though is −h σ ′ ∂ x u ′ i ,which also exhibits a dual trend, negligible however incomparison to the components presented in Fig. 15.Finally, the correlations in Figs. 15a and 15c are alsoresponsible for the transitional behaviour of viscoelasticdissipation profiles from LDR to HDR discussed earlier.The current picture of the dual nature of ε P was first −3 y/ δ − < σ , ∂ x u , > δ ν / u τ Case ACase BCase DCase GCase H y/ δ − < σ , ∂ x u , > δ ν / u τ −3 y/ δ − < σ , ∂ x u , > δ ν / u τ FIG. 15: (Colour online) Profiles of viscoelastic dis-sipation components for the LDR, HDR and MDRregimes. (a) −h σ ′ ∂ x u ′ i δ ν /u τ , (b) −h σ ′ ∂ y u ′ i δ ν /u τ and(c) −h σ ′ ∂ z u ′ i δ ν /u τ as function of y/δ . Note: case N2(DR = 0%); case A (DR = − . − . − . − . − . predicted by Min et al. [56] at low Weissenberg num-6bers, adding an artificial diffusion term to numericallysolve the FENE-P model. However, the present DNSare the first to capture so clearly these regions through-out the drag reduction regimes, predicting the appropri-ate dynamics at corresponding We c values. Once more,this is attributed to our numerical approach applied herefor the FENE-P model that is able to capture strongerpolymer-turbulence interactions than algorithms basedon artificial diffusion. There are even results using theartificial diffusion methodologies that erroneously predictpolymers never feeding energy back to the flow [9, 32].Hence, in view of the current distinctly transparent ob-servations a conceptual model for the mechanism of dragreduction is deduced in the next section. VII. DRAG REDUCTION MECHANISM
The recent review on polymer drag reduction by Whiteand Munghal [2] reports that the numerical evidence issomewhat conflicting regarding the flow regions wherepolymers extend and contract. In this study, these re-gions can be identified by applying the Reynolds decom-positions u i = h u i i + u ′ i and σ ij = h σ ij i + σ ′ ij to Eq. (9),following the spirit of [32, 56]. Then, we can notice that h σ ′ ij ∂ x j u ′ i i appears as a production term due to turbu-lence for the mean polymer elastic energy. Now, fromthe definition of polymer elastic energy Eq. (6), it isevident that h E p i ∝ ln h f ( C kk ) i . So, the FENE-P dumb-bells are stretched when − ε P δ ν /u τ < y from the wall.According to Figs. 12c and 13 there are three mainregions in the profiles of viscoelastic dissipation − ε P δ ν /u τ < , ≤ y/δ < δ (We c , Re c ) > , δ (We c , Re c ) ≤ y/δ ≤ δ (We c , Re c ) < , δ (We c , Re c ) < y/δ ≤ . (22)The first region is at the proximity of the wall, wherepolymers unravel because of the high mean shear, con-sistent with other studies [42, 56–58], storing elas-tic potential energy. The range of this region has aweak dependence on Weissenberg and Reynolds numberwith its upper bound being within the viscous sublayer δ (We c , Re c ) . .
05 for all We c and Re c cases considered.The second region is the most interesting since poly-mers release energy back to the flow, contracting towardstheir equilibrium length, as they are convected away fromthe wall by the near-wall vortical motions. The mani-festation of turbulence production by polymers can beinterpreted in terms of the correlation of the polymerswith the local fluctuating strain rates and their persis-tence in this region. In particular, −h σ ′ ij ∂ x j u ′ i i revealsthat −h σ ′ ∂ x u ′ i as well as −h σ ′ ∂ z u ′ i are responsible for the contraction of the dumbbells and consequently forthe release of the stored elastic energy, since they arepositively correlated in this region away from the wall(see Figs. 15a and 15c). This region exists in an in-termediate y/δ range, whose upper bound δ (We c , Re c )is strongly dependent on We c and less on Re c values.As drag reduction amplifies for larger polymer relaxationtime scales this positive region expands to a wider y/δ range, which dominates the nature of − ε P δ ν /u τ at MDR(see Fig. 12c).Finally, polymers transported away from the wall getalso negatively correlated with the persistent fluctuatingstrain rates (see Fig. 15b) and are extended in the region δ (We c , Re c ) < y/δ ≤
1, which is a sink for turbulentkinetic energy, prevailing the LDR flows. However, thisregion is diminished for HDR and MDR flows (see Fig.12c) due to the interplay between the productive anddissipative inherent features of ε P , which mainly dependon the polymer relaxation time scale and the existence ofintense velocity fluctuations that are able to stretch thepolymer molecules.The phenomenology of the proposed mechanism sharesmany similarities with various conceptual models of ear-lier works [2]. In this study, the basic idea is that thetransport of the elastic potential energy, stored by poly-mers near the wall, is mainly associated with the polymerrelaxation time scale. The latter determines the distribu-tion of energy away from the wall and as a consequencethe near-wall turbulence dynamics weaken. Up to thispoint, the mechanism agrees with the interpretation ofMin et al. [56], which is essentially confirmed by thepresent illustrative computations. However, the noveltyhere is that this mechanism is valid for higher We c valuesand levels of percentage DR in contrast to Min et al. [37],who claim that it is not valid for HDR/MDR flows bas-ing their arguments on their debatable numerical results(see also section III C).In addition, the refinement of the proposed con-ceptual mechanism resides on the reduction of ε P to h σ ′ j ∂ x j u ′ i and even more on the correlations h σ ′ ∂ x u ′ i and h σ ′ ∂ z u ′ i , which are responsible for the turbulenceproduction by polymer coils. The existence of a thirddissipative region away from the wall is also emphasisedin this mechanism, where polymers, after their contrac-tion, are now stretched by the intense fluctuating veloc-ity field. This outer region dominates the viscoelasticdissipative dynamics of the LDR regime and diminishesasymptotically as We c increases but it never disappears.Ultimately, this picture along with the anisotropy intro-duced into the components of turbulent kinetic energy,i.e. E = ( u ′ + v ′ + w ′ ), comprise the drag reductionmechanism deduced in this study. VIII. CONCLUSION
This paper is devoted to the polymer dynamics in vis-coelastic turbulent channel flow and their effects on the7flow, reproducing turbulent drag reduction by DNS us-ing a state-of-the-art numerical scheme in wall-boundedflows to solve the FENE-P model. The potential of thismethodology to capture the strong polymer-turbulencedynamical interactions allowed β values to remain high,more representative of dilute polymer solutions used inexperiments. Even then, higher percentage DR values areobtained for given We c than previous numerical studies.The effects of L p and Re c on the results support theclaims for non-universality of the dynamics for interme-diate levels of DR between the von K´arm´an and theMDR law. The universal MDR asymptote, on the otherhand, is reached in this study under the combination ofhigh polymer extensibility L p with high enough elasticitygiven by large values of We c at a given moderate Re c .The experimentally observed distinct differences in thestatistical trends of the turbulent velocity field, particu-larly for u ′ + (see Fig. 5a), are clearly identified withthe current numerical approach in comparison with othersimulations, most of which do not even approach such acharacteristic trend. Overall, the peaks of the statisticalprofiles of velocity and vorticity fluctuations shift awayfrom the wall as DR increases, in agreement with otherexperimental and numerical studies, indicating the thick-ening of the buffer layer. At the same time, νβ dd y h u i in-creases towards the centre of the channel for higher We c ,denoting the importance of viscosity away from the wallat these moderate Reynolds number DNS.Lumley’s phenomenology Lumley [8] on the manifes-tation of drag reduction is based on the conjecture ofcoil-stretch transition, i.e. exponential full uncoiling ofpolymer molecules, for the build-up of intrinsic viscos-ity. However, our numerical results illustrate that theonset of drag reduction and even the MDR asymptoticstate can be reached while h C kk i ≪ L p with L p largeenough, in agreement with the initial claim by Taborand de Gennes Tabor and de Gennes [10] that even highspace-time strain rate fluctuations near the wall can onlypartially stretch polymer coils. We also showed that thepercentage polymer extension is less but the actual ex-tension is more for larger L p , amplifying DR. Thus, largepolymer coils that do not reach their critical full exten-sibility should be of interest to experimental investiga-tions on scission degradation of polymer chains and dragreduction effectiveness. Such macromolecules would beless vulnerable to rupture avoiding the loss of the dragreduction effect. Besides, they should be able to stretchsubstantially to make a stronger impact on turbulent ac-tivity and consequently enhance percentage drag reduc-tion.The analysis of the conformation tensor field providesgreat insight into the polymer dynamics and their influ-ence on the flow. The dominant anisotropic behaviourof the mean conformation tensor, i.e. h C i ≫ h C i ≃h C i > h C i , due to the mean shear in viscoelastic tur-bulent channel flow, influences the anisotropy of the fluc-tuating flow field. The anisotropy in the HDR and MDRregimes is depicted at the small scales of our DNS outside the buffer layer and towards the centre of the channel by ω ′ z + > ω ′ y + > ω ′ x + .Different asymptotic rates of convergence are observedfor the conformation tensor components towards the limitof infinite Weissenberg number demonstrating the com-plex polymer dynamics even in this simplified dumbbellmodel. In the limit We S → ∞ polymers are consid-ered stiff, i.e. C ij → h C ij i , mostly in the main di-rections of elongation and the correlations of the fluc-tuating conformation tensor and velocity fields in theother directions are assumed to remain minimal at thislimit. Therefore, h σ i = A − β Re c h C i dd y h u i and h σ i = A − β Re c h C i dd y h u i , with A → A → A → A . A on the other hand isabout contant in the range 0 . . y/δ . . S increases.The following theoretical view was stated in this pa-per with regards to the controversy over the existence ornot of Reynolds shear stress at the MDR limit, whichis of fundamental importance to the dynamics of tur-bulence production at this limit. It is conjectured thatat the MDR limit h σ i is negligible. This was basedon the idea mentioned above about the stiffness of poly-mers at We S → ∞ plus the assumption that h C i be-comes negligible at the same limit. Then, it is supposedthat this behaviour is also valid under both the infiniteWeissenberg and Reynolds number limits by taking care-fully these limits, so that we go along the universal MDRasymptotic line. Hence, one might expect an interme-diate region δ ν ≪ y ≪ δ of approximately constantReynolds shear stress, i.e. −h u ′ v ′ i /u τ →
1, implied bythe balance of shear stresses when taking the limits of y/δ → y/δ ν → ∞ with the reasonable assumptionthat νβ dd y h u i → y ≫ δ ν . In summary, the classicalturbulence generation by −h u ′ v ′ i seems to survive at theMDR limit, based on the above assumptions.Polymer-turbulence dynamical interactions were ex-pressed through viscoelastic dissipation ε P ≡ h σ ′ ij ∂ x j u ′ i i ,which can either dissipate or produce turbulent kineticenergy. For HDR and MDR flows, R ε P d y becomes vitalin the flow dynamics in proportion to R P d y and R ε N d y due to the vast inhibition of Reynolds shear stress andfluctuating strain rates, respectively. In particular, adifferent view of the near-wall dissipation dynamics isshown for HDR/MDR flows, with the maximum dissipa-tion arising away from the wall. It is intriguing to notethat ε P follows a transitional pattern from LDR to HDRregime (see Fig. 12c) similar to u ′ + (see Fig. 5a) and h C i (see Fig. 8a). This characteristic behaviour is alsoreproduced on average in R ε P d y , where its dissipativefeature enhances in the LDR regime but attenuates forHDR/MDR flows, with the productive nature dominat-ing for high percentage drag reduction. Thus, polymersget somehow involved in the production dynamics of tur-bulent kinetic energy.8In view of the current viscoelastic DNS the followingconceptual picture of drag reduction is deduced, whichis an extension to and refinement of the mechanism pro-posed by Min et al. [56]. Polymers in the near-wall re-gion extract energy from the flow due to the uncoilingcaused by the mean shear and release some portion of thisstored elastic energy back to the flow by contracting asthey move away from the wall. This transport of energydepends on Weissenberg number which determines thedistribution of energy away from the wall. Ultimately,this process undermines the dynamics of near-wall tur-bulence. Note that polymers also unravel due to velocityfluctuations, as they move towards the core region of theflow, extracting again energy from the flow. This mech-anism appears to be valid for all drag reduction regimeswith the dissipative and productive elements of viscoelas-tic dissipation competing in different parts of the flow fordifferent levels of DR. We also observe that correlation h σ ′ j ∂ x j u ′ i is able to resemble the dynamics of ε P andspecifically that h σ ′ ∂ x u ′ i and h σ ′ ∂ z u ′ i are the correla-tions responsible for the production of turbulent kineticenergy by polymers.So far, in the limited context of the FENE-P modeland at moderate Reynolds number DNS, the proposedphenomenology agrees with the majority of experimentaland numerical data, where dampening of near-wall turbu-lence has long been speculated with various analyses andinterpretations. Here, however, the transfer of energyfrom the flow to the polymers, its redistribution by thelatter in the flow field and the prevalence of anisotropyover the components of E ≡ h| u ′ | i in the three Carte-sian directions is suggested as a possible cause of dragreduction. Acknowledgments
We are grateful to S. Laizet for providing the Navier-Stokes solver and to Halliburton for the financial support.We would also like to thank J. G. Brasseur, L. R. Collinsand T. Vaithianathan for useful discussions.
Appendix: Computational method1. Numerical method for the FENE-P model
The numerical scheme adapted here for non-periodicboundary conditions was initially developed by Vaithi-anathan et al.
Vaithianathan et al. [21] for periodic do-mains. The main idea behind the high-resolution cen-tral schemes employed here is the use of higher-orderreconstructions, which enable the decrease of numericaldissipation so as to achieve higher resolution of shocks.In essence, they employ more precise information ofthe local propagation speeds. A key advantage of cen-tral schemes is that one avoids the intricate and time-consuming characteristic decompositions based on ap- proximate Riemann solvers [20]. This is because theseparticular schemes realise the approximate solution interms of its cell averages integrated over the Riemannfan (see Fig. 16). i+1/2,j,k C + i+1/2,j,k CH i+1/2,j,k − FIG. 16: Central differencing approach – staggered integra-tion over a local Riemann fan denoted by the dashed-doubledotted lines.
Considering the discretisation of the convection term ofthe FENE-P model only in the x -direction, using the re-construction illustrated in Fig. 16, the following second-order discretisation is obtained ∂ C ni,j,k ∂x = 1∆ x ( H ni +1 / ,j,k − H ni − / ,j,k ) (A.1)where H ni +1 / ,j,k = 12 u i +1 / ,j,k ( C + i +1 / ,j,k + C − i +1 / ,j,k ) − | u i +1 / ,j,k | ( C + i +1 / ,j,k − C − i +1 / ,j,k ) (A.2)with C ± i +1 / ,j,k = C ni +1 / ± / ,j,k ∓ ∆ x · ∂ C ∂x (cid:12)(cid:12)(cid:12)(cid:12) ni +1 / ± / ,j,k (A.3)and ∂ C ∂x (cid:12)(cid:12)(cid:12)(cid:12) ni,j,k = x ( C ni +1 ,j,k − C ni,j,k ) x ( C ni,j,k − C ni − ,j,k ) x ( C ni +1 ,j,k − C ni − ,j,k ) . (A.4)Similarly, Eqs. (A.2)-(A.4) can be rewritten for H ni − / ,j,k . The appropriate choice of the derivative dis-cretisation in Eq. (A.4) limits the slope so that the SPDproperty for C is satisfied. The SPD criterion for thischoice is that all the eigenvalues of the conformation ten-sor should be positive, viz. λ i > C ) >
0, is not sufficient9to guarantee the SPD property for the tensor [59]. In casenone of the options in Eq. (A.4) satisfy the criterion, thenthe derivative is set to zero reducing the scheme to firstorder locally in space. The proof for C being SPD us-ing this numerical scheme can be found in Vaithianathanet al. [21]. The eigenvalues of the conformation tensor inthis implementation are computed using Cardano’s an-alytical solution [60] for the cubic characteristic polyno-mial avoiding any complicated and time-consuming linearalgebra matrix decompositions and inversions for just a3 × et al. Vaithianathan et al. [61], since they only consideredperiodic boundary conditions. So, the implementation ofthe numerical method near the walls of the channel wasmodified for this study considering ghost nodes beyondthe wall boundaries to keep the original formulation un-altered, preserving in that way the second-order accuracyat the boundaries. The values at the ghost nodes werelinearly extrapolated from the interior solution [20], i.e. C ni,j +1 ,k = C ni,j,k + ( C ni,j,k − C ni,j − ,k )= 2 C ni,j,k − C ni,j − ,k . (A.5)The time advancement is done simply using the for-ward Euler update, treating implicitly the third, the forthterm on the left hand side and the right hand side of Eq.(5) due to the potential finite extensibility of the poly-mer. Hence, the fully discretised form of the FENE-Pmodel is C n +1 i,j,k = C ni,j,k − ∆ t ∆ x ( H ni +1 / ,j,k − H ni − / ,j,k ) − ∆ t ∆ y j ( H ni,j +1 / ,k − H ni,j − / ,k ) − ∆ t ∆ z ( H ni,j,k +1 / − H ni,j,k − / )+ ∆ t ( C n +1 i,j,k ∇ u ni,j,k + ∇ u n ⊤ i,j,k C n +1 i,j,k ) − ∆ t (cid:18) c f ( C n +1 i,j,k ) C n +1 i,j,k − I (cid:19) (A.6)with C ni,j,k = 16 ( C − i +1 / ,j,k + C + i − / ,j,k + C − i,j +1 / ,k + C + i,j − / ,k + C − i,j,k +1 / + C + i,j,k − / ) (A.7) so that the convection term and the explicit term comingfrom the time derivative can be assembled in a convexsum C ∗ = C ni,j,k + ∂ C ni,j,k ∂ x = N X l =1 s l C l (A.8)where all coefficients s l ≥ P Nl =1 s l = 1, with C ∗ being SPD if the matrices C l are SPD, ensuring thefinite extensibility of the dumbbell, i.e. the trace of theconformation tensor is bounded tr C = λ + λ + λ ≤ L P [21]. The following CFL condition needs to be satisfiedfor the coefficients s l to be non-negativeCFL = max (cid:26) | u | ∆ x , | v | ∆ y min , | w | ∆ z (cid:27) · ∆ t <
16 (A.9)and it also determines the time step ∆ t . Note that thisCFL condition is more strict than the one for compactfinite differences [26] used for Newtonian turbulence com-putations.The numerical solution of Eq. (A.6) is carried out byfirst rewriting it in a Sylvester-Lyapunov form [62], sep-arating the implicit and explicit terms, i.e. A ⊤ X + XA = B ⇒ ( I ⊗ A ⊤ + A ⊤ ⊗ I ) x = b (A.10)where A ≡ [1 + f ( C n +1 i,j,k ) ∆ t We c ] I − ∆ t ∇ u ni,j,k , X ≡ C n +1 i,j,k and B ≡ C ∗ + ∆ t We c I are 3 × I ⊗ A ⊤ + A ⊤ ⊗ I )is a 9 × x ≡ vec( X ), b ≡ vec( B ) are9 × × ×
2. Time advancement
After obtaining the new update of the conformationtensor C n +1 i,j,k , the two-step, i.e. three time-level, second-order Adams-Bashforth/Trapezoidal scheme is used forthe time integration of Eqs. (3) through the followingprojection method [64] u ∗ − u n ∆ t = 12 (3 F n − F n − ) + 12 ( P ∗ n +1 + P n ) (A.11) u n +1 − u ∗ ∆ t = − ∇ ˜ p n +1 (A.12)where F = −
12 [ ∇ ( u ⊗ u ) + ( u · ∇ ) u ] + 1Re c ∆ u (A.13)0and P = 1 − β Re c We c ∇ · L p − L p − tr C C − I ! (A.14)with ˜ p n +1 = 1∆ t Z t n +1 t n p dt. (A.15)The incompressibility condition ∇ · u n +1 = 0 is verifiedby solving the Poisson equation ∇ · ∇ ˜ p n +1 = ∇ · u ∗ ∆ t (A.16) which is done in Fourier space [25]. It is well known thatthese multistep methods are not self-starting and requirea single-step method to provide the first time level [64,65]. In this study, explicit Euler was chosen for just thefirst iteration of these computations, viz. u n = u n − +∆ t F n − . [1] B. A. Toms (North-Holland, Amsterdam, 1948), vol. 2 of Proceedings 1st International Congress on Rheology , pp.135–141.[2] C. M. White and M. G. Mungal, Ann. Rev. Fluid Mech. , 235 (2008).[3] H. M. Nagib and K. A. Chauhan, Phys. Fluids , 101518(2008).[4] V. Dallas, J. C. Vassilicos, and G. F. Hewitt, Phys. Rev.E , 046306 (2009).[5] M. D. Warholic, H. Massah, and T. J. Hanratty, Exp.Fluids , 461 (1999).[6] P. S. Virk, E. W. Merrill, H. S. Mickley, K. A. Smith, andE. L. Mollo-Christensen, J. Fluid Mech. , 305 (1967).[7] P. S. Virk, AIChE Journal , 625 (1975).[8] J. L. Lumley, Ann. Rev. Fluid Mech. , 367 (1969).[9] I. Procaccia, V. S. L’vov, and R. Benzi, Rev. Mod. Phys. , 225 (pages 23) (2008).[10] M. Tabor and P. G. de Gennes, Europhys. Lett. pp. 519– 522 (1986).[11] D. D. Joseph, Fluid dynamics of viscoelastic liquids (Springer-Verlag, 1990).[12] K. R. Sreenivasan and C. M. White, J. Fluid Mech. ,149 (2000).[13] V. Dallas, Ph.D. thesis, Imperial College London (2010).[14] H. C. ¨Ottinger,
Stochastic processes in polymeric fluids (Springer Berlin, 1996).[15] R. B. Bird, C. F. Curtis, R. C. Armstrong, and O. Has-sager,
Dynamics of Polymeric Liquids , vol. 2: KineticTheory (John Wiley and Sons, 1987).[16] R. Sureshkumar and A. N. Beris, J. Non-NewtonianFluid Mech. , 53 (1995).[17] T. Min, J. Y. Yoo, and H. Choi, J. Non-Newtonian FluidMech. , 27 (2001).[18] S. Jin and L. R. Collins, New Journal of Physics , 360(2007).[19] D. Joseph and J. Saut, J. Non-Newtonian Fluid Mech. , 117 (1986).[20] R. J. LeVeque, Finite volume methods for hyperbolicproblems (Cambridge University Press, 2002).[21] T. Vaithianathan, A. Robert, J. G. Brasseur, and L. R.Collins, J. Non-Newtonian Fluid Mech. , 3 (2006).[22] A. Kurganov and E. Tadmor, J. Comput. Phys. , 241(2000).[23] A. Peterlin, Polymer , 257 (1961). [24] H. R. J. Warner, Ind. Eng. Chem. Fundam. , 379(1972).[25] S. Laizet and E. Lamballais, J. Comput. Phys. , 5989(2009).[26] S. K. Lele, J. Comput. Phys. , 16 (1992).[27] A. B. Cain, J. H. Ferziger, and W. C. Reynolds, J. Com-put. Phys. , 272 (1984).[28] R. B. Dean, ASME J. Fluids Eng. , 215 (1978).[29] M. Lesieur, Turbulence in Fluids (Kluwer Academic Pub-lishers, 1997).[30] P. Moin and J. Kim, J. Comput. Phys. , 381 (1980).[31] C.-F. Li, R. Sureshkumar, and B. Khomami, J. Non-Newtonian Fluid Mech. , 23 (2006).[32] P. K. Ptasinski, B. J. Boersma, F. T. M. Nieuwstadt,M. A. Hulsen, H. A. A. Van den Brule, and J. C. R.Hunt, J. Fluid Mech. , 251 (2003).[33] J. Jim´enez and P. Moin, J. Fluid Mech. , 213 (1991).[34] R. D. Moser, J. Kim, and N. N. Man-sour, Phys. Fluids , 943 (1999), URL http://turbulence.ices.utexas.edu .[35] K. Iwamoto, Y. Suzuki, and N. Kasagi, Int.J. Heat Fluid Flow , 678 (2002), URL .[36] Z. W. Hu, C. L. Morfey, and N. D. Sand-ham, AIAA Journal , 1541 (2006), URL .[37] T. Min, H. Choi, and J. Y. Yoo, J. Fluid Mech. , 91(2003).[38] C. D. Dimitropoulos, R. Sureshkumar, and A. N. Beris,J. Non-Newtonian Fluid Mech. , 433 (1998).[39] K. D. Housiadas and A. N. Beris, Phys. Fluids , 2369(2003).[40] J. L. Lumley, J. Polym. Sci. Macromol. Rev. , 263(1973).[41] P. K. Ptasinski, F. T. M. Nieuwstadt, B. H. A. A. van denBrule, and M. A. Hulsen, Flow, Turbul. Combust. ,159 (2001).[42] Y. Dubief, M. C. White, V. E. Terrapon, E. S. G.Shaqfeh, P. Moin, and S. K. Lele, J. Fluid Mech. ,271 (2004).[43] R. A. Handler, K. D. Housiadas, and A. N. Beris, Int. J.Numer. Methods Fluids , 1339 (2006).[44] J. Kim, P. Moin, and R. Moser, J. Fluid Mech. , 133(1987). [45] R. Sureshkumar, A. N. Beris, and R. A. Handler, Phys.Fluids , 743 (1997).[46] K. Kim, C.-F. Li, R. Sureshkumar, S. Balachandar, andR. J. Adrian, J. Fluid Mech. , 281 (2007).[47] K. Kim, R. J. Adrian, S. Balachandar, and R. Sureshku-mar, Phys. Rev. Lett. , 134504 (2008).[48] H. T. Kim, S. J. Kline, and W. C. Reynolds, J. FluidMech. , 133 (1971).[49] S. B. Pope, Turbulent Flows (Cambridge UniversityPress, 2000).[50] C. D. Dimitropoulos, Y. Dubief, E. S. G. Shaqfeh,P. Moin, and S. K. Lele, Phys. Fluids , 011705 (2005).[51] P. G. De Gennes, Introduction to polymer dynamics (Cambridge University Press, 1990).[52] V. S. L’vov, A. Pomyalov, I. Procaccia, and V. Tiberke-vich, Phys. Rev. E , 016305 (2005).[53] R. Benzi, E. De Angelis, V. S. L’vov, I. Procaccia, andV. Tiberkevich, J. Fluid Mech. , 185 (2006).[54] H. Tennekes and J. Lumley, A first course in turbulence (MIT press, 1972).[55] U. Frisch,
Turbulence: the legacy of A. N. Kolmogorov (Cambridge university press, 1995).[56] T. Min, J. Y. Yoo, H. Choi, and D. D. Joseph, J. FluidMech. , 213 (2003). [57] H. Massah and T. J. Hanratty, J. Fluid Mech. , 67(1997).[58] V. E. Terrapon, Y. Dubief, P. Moin, E. S. G. Shaqfeh,and S. K. Lele, J. Fluid Mech. , 61 (2004).[59] G. Strang,
Linear algebra and its applications (Thomp-son Learning, Inc., 1988).[60] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.Flannery,
Numerical recipes in Fortran 77 (CambridgeUniversity Press, 1996).[61] T. Vaithianathan, J. G. Brasseur, and L. R. Collins,
Per-sonal communication (Cornell University, USA, 2007).[62] K. B. Petersen and M. S. Pedersen,
The Matrix Cookbook (Technical University of Denmark, 2008).[63] J. Dennis and R. Schnabel,
Numerical methods for uncon-strained optimization and nonlinear equations (Prentice-Hall, 1983).[64] R. Peyret,
Spectral methods for incompressible viscousflow (Springer-Verlag, 2002).[65] R. J. LeVeque,