On the skin friction coefficient in viscoelastic wall-bounded flows

International Journal of Heat and Fluid Flow 42 (2013) 49–67

Contents lists available at SciVerse ScienceDirect

International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff

On the skin friction coefficient in viscoelastic wall-bounded flows Kostas D. Housiadas a, Antony N. Beris b,⇑ a b

Department of Mathematics, University of the Aegean, Karlovassi 83200, Samos, Greece Department of Chemical & Biomolecular Engineering, University of Delaware, Newark, DE 19716, USA

a r t i c l e

i n f o

Article history: Received 25 February 2012 Received in revised form 13 November 2012 Accepted 16 November 2012 Available online 5 March 2013 Keywords: Drag reduction Skin friction coefficient Dilute polymer solutions Turbulent channel flow FENE-P model Giesekus model Oldroyd-B model

a b s t r a c t Analysis of the skin friction coefficient for wall bounded viscoelastic flows is performed by utilizing available direct numerical simulation (DNS) results for viscoelastic turbulent channel flow. The Oldroyd-B, FENE-P and Giesekus constitutive models are used. First, we analyze the friction coefficient in viscous, viscoelastic and inertial stress contributions, as these arise from suitable momentum balances, for the flow in channels and pipes. Following Fukagata et al. (Phys. Fluids, 14, p. L73, 2002) and Yu et al. (Int. J. Heat. Fluid Flow, 25, p. 961, 2004) these three contributions are evaluated averaging available numerical results, and presented for selected values of flow and rheological parameters. Second, based on DNS results, we develop a universal function for the relative drag reduction as a function of the friction Weissenberg number. This leads to a closed-form approximate expression for the inverse of the square root of the skin friction coefficient for viscoelastic turbulent pipe flow as a function of the friction Reynolds number involving two primary material parameters, and a secondary one which also depends on the flow. The primary parameters are the zero shear-rate elasticity number, El0, and the limiting value for the drag reduction at high Weissenberg number, LDR, while the secondary one is the relative wall viscosity, lw. The predictions reproduce both types A and B of drag reduction, as first introduced by Virk (Nature, 253, p. 109, 1975), corresponding to partially and fully extended polymer molecules, respectively. Comparison of the results for the skin friction coefficient against experimental data shows good agreement for low and moderate drag reduction which is the region covered by the simulations. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction Wall-bounded turbulent flows in channels and pipes are of great importance for many practical and engineering applications. A particularly interesting case is that of flows involving viscoelastic liquids, such as dilute and semi-dilute polymer solutions. Indeed, it is now well-known that under turbulent conditions, the addition of even minute (ppm range) quantities of a polymer in a low molecular weight solvent alters the flow dynamics resulting to a significant decrease of the pressure drop required to drive the flow (Virk, 1975a). This phenomenon is referred to as drag reduction by polymer additives and it has been the subject of numerous experimental and numerical studies (Gyr and Bewersdorff, 1995; White and Mungal, 2008; Beris and Housiadas, 2010; and references therein). Of particular interest, are the experimental observations that there is a maximum level of drag reduction that can be achieved by polymer additives (usually called ‘‘Virk’s maximum drag reduction (MDR) asymptote’’, Virk et al., 1967) as well as that they can be two types, A and B, of polymer drag reduction (Virk, 1975b; Virk and Wagger, 1990; Gyr and Bewersdorff, 1995). The former is asso⇑ Corresponding author. E-mail address: [email protected] (A.N. Beris). 0142-727X/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.ijheatfluidflow.2012.11.004

ciated with low and moderate levels of polymer extension by the flow while the latter corresponds to fully extended molecules. Also, of interest is the observation that there is usually a minimum flow deformation required for drag reduction, typically presented as a minimum Reynolds or Weissenberg number (Virk et al., 1967; Paterson and Abernathy, 1970). This experimental result has also been confirmed by DNS (Min et al., 2003; Housiadas and Beris, 2003). Significant progress has been made on the understanding of the mechanism underlying polymer-induced drag reduction. The two main theories that have been developed attempt to explain the phenomena based on either molecular (De Gennes, 1986; Tabor and De Gennes, 1986) or fluid-mechanical (Lumley, 1969; Seyer and Metzner, 1969) considerations. However, although we can say that the onset and intermediate drag reduction phenomena have been mostly understood, this is not the case regarding the high and maximum drag reduction regions, which are still an active area of investigation (Xi and Graham, 2012). More recently, direct numerical simulations (DNSs) of dilute and semi-dilute polymer solutions flows have been developed that further contributed on the understanding of the polymer-induced drag reduction phenomenon (Sureshkumar et al., 1997; DeAngelis et al., 2002; Min et al., 2003; Housiadas and Beris, 2003, 2004a; Yu and Kawaguchi, 2003, 2006; Dubief et al., 2004; Dimitropoulos

50

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

et al., 2006; Li et al., 2006a; Li et al., 2006b). DNS have been proved very useful because they provide significant information regarding not only the detailed structure of the flow but also that of the polymer molecules. Based on the DNS results, it appears that as the viscoelasticity in the flow increases the eddies formed under turbulent conditions are larger in size and take a longer time to develop, reducing the momentum transfer. It is also clear the increase of the width of the streaky structures proportionally to the drag reduction (Sureshkumar et al., 1997; Housiadas et al., 2005), which compares well with experiments (Luchik and Tiederman, 1988). However, no systematic effort has been made so far to study the quantitative influence of viscoelasticity on the skin friction coefficient. The latter is defined as the ratio of the total wall shear stress to a characteristic bulk kinetic energy of the fluid, and it is of importance for engineering applications. Thus, we utilize here our available DNS results using a variety of well-known constitutive equations such as the Oldroyd-B, Giesekus and FENE-P models, and present correlations among the model parameters and the skin friction coefficient. First, a systematic analysis of the skin friction coefficient on viscous, inertial and viscoelastic stress contributions is conducted, following the work of Fukagata et al. (2002) for Newtonian wall bounded flows and the work of Yu et al. (2004) for surfactant-induced drag reduction. The former authors derived the skin friction coefficient in terms of inertia, viscous and uniform blowing/suction on the wall contributions in order to evaluate the effect of the latter in controlling turbulence. Yu et al. (2004) performed a similar analysis for viscoelastic turbulent channel flow and analyzed the skin friction coefficient in terms of viscous, inertial, and viscoelastic contributions. They presented results using the Giesekus model only in a single case, where they showed that, although the addition of viscoelasticity introduces a new contribution to the skin friction coefficient, this term is much smaller than the associated reduction to the inertial contribution and, as a consequence, drag reduction results. Even more recently Tsukahara et al. (2011), presented additional results for drag-reduced surfactant solutions in turbulent channel flow. These authors used the Giesekus model as well, and presented the fractional contribution to the skin friction coefficient for five typical cases. They found that whereas the fractional viscoelastic contribution increases with increasing friction Weissenberg number, the inertial one decreases. Second, we represent the previously evaluated drag reduction results with a universal correlation as a function of the friction Weissenberg number. This curve only involves one scaling parameter, the limiting drag reduction at high Weissenberg numbers, LDR. This new material parameter can be evaluated from DNS or, in the absence of DNS, can be treated as a fitting parameter. Based on that, we develop a new parametric relationship for the skin friction coefficient for viscoelastic fluids most naturally in Prandtl– Karman coordinates, i.e. representing the inverse of the square root of the friction factor as function of the friction Reynolds number. Using a more accurate expression of the skin friction coefficient for low and intermediate Reynolds Newtonian turbulent pipe flow that we have obtained through integration of the best available expression originated from experimentally determined mean velocity data, we established an equivalent expression for dragreducing viscoelastic turbulent pipe flows. Comparison of the theoretical predictions with experimental data from the literature (Shetty and Solomon, 2009) shows good agreement provided that the relaxation time of the polymer, estimated from simple shear data, is suitably rescaled in order to better reflect the response of the polymer under extensional deformation encountered in turbulent flows. Last, we note that, as part of the definition of the friction or bulk Reynolds number as well as the friction Weissenberg number, the viscosity of the solution at the wall is required. This can be

approximated by the zero shear-rate value, except when significant shear thinning is present, in which case it can be estimated with a Carreau fitting of simple shear data. The rest of the paper is organized as follows. In Section 2 we present the basic definitions and equations together with the analysis of the skin friction coefficient for straight channels for laminar flow followed by an analysis for turbulent channel flow. In Section 3 we offer an analysis of the skin friction coefficient in terms of its viscous, inertial and viscoelastic contributions based to the averaged in space and time results of a selected set of DNS cases. In Section 4, we offer parametric relations of the drag reduction and the skin friction coefficient. In Section 5 we evaluate the skin friction coefficient for turbulent pipe flow based on mean velocity profiles for Newtonian fluids, as well as viscoelastic fluids at maximum drag reduction. In Section 6 we explore the possibility of a vanishing Reynolds stress at the maximum drag reduction state in both channel and pipe flows. Finally, in Section 7 we state our conclusions. 2. Analysis of the skin-friction coefficient For pipe and channel flows, the skin friction coefficient, or Fanning friction factor, Cf, is defined as the ratio of the wall stress, sw , applied by the flow to the wall, to an indicative kinetic energy:

sw ;  u2 q b 2

Cf  1

ð1Þ

where q⁄ is the constant mass density of the fluid, and ub is its bulk average velocity; here and in the following a superscript ⁄ is used to denote a dimensional quantity. The bulk velocity for pipe and channel flow is defined as:

ub 

2

Z

R



ux r  dr ; pipe flow; R2 0 Z h 1  u dy ; channel flow; ub   2h h x

ð2aÞ ð2bÞ

where R⁄ is the radius of the pipe, and h⁄ is the half channel-height. In Eq. (2), ux is the velocity component along the main flow direction. The skin friction coefficient is a function of the bulk Reynolds number. Due to the additional turbulent momentum transfer caused by eddies Cf is much higher in turbulent flows than in laminar flows. Note that in the literature, the Darcy friction factor, f, is also often used; the Darcy friction factor is simply four times the Fanning friction factor, f = 4Cf. In turbulent flows, it is also customary to define the friction velocity us :

us 

qffiffiffiffiffiffiffiffiffiffiffiffiffi sw =q :

ð3Þ

Combining definitions (1) and (3) gives:

ub ¼ us

sffiffiffiffiffi 2 : Cf

ð4Þ

The bulk and friction velocity are used to define the zero shearrate bulk and friction Reynolds numbers, respectively:

Reb0  Res0 

c‘ ub

m0 ‘ us

m0

;

;

ð5Þ ð6Þ

where m0  g0 =q is the zero shear-rate kinematic viscosity of the fluid, g0 is the zero shear-rate viscosity and ‘⁄ is a characteristic length scale; for pipes, ‘⁄ = R⁄, and for channels ‘⁄ = h⁄. The constant

51

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

c has been introduced for later convenience since for pipe flow c = 2, while various values are used for channel flow. By combining definitions (5) and (6) we have:

Reb0 ¼c Res0

sffiffiffiffiffi 2 ; Cf

ð7aÞ

or, alternatively

Cf ¼

Re2 2c2 s20 Reb0

ð7bÞ

In flow systems for which the shear viscosity is not constant we can also express the friction factor as C f ¼ 2c2 Re2s =Re2b , where Res and Reb are the wall friction and bulk Reynolds numbers, respectively. Those are defined analogously to Eq. (6) using the wall kinematic viscosity, mw  gw =q , instead of m0 (Housiadas and Beris, 2004b). Note that any quantity that is defined using a reference viscosity can be defined either with respect to the zero shear-rate value, g0 (and in this case a zero subscript is used) or in terms of its value at the wall, gw . The relationship between these two is conveniently described by the dimensionless wall viscosity lw  gw =g0 (Housiadas and Beris, 2004b; see also Section 4). For a general viscoelastic turbulent flow, the friction factor can be further analyzed in terms of its viscous, viscoelastic and inertial contributions. First, we derive the corresponding expressions for a straight channel under steady and laminar flow conditions (Section 2.1), and then under stationary and turbulent flow conditions (Section 2.2). The same analyses for steady and turbulent flows in a circular pipe is offered in Appendices A and B, respectively. 2.1. Steady laminar viscoelastic flow in a channel Let us use a Cartesian coordinate system x⁄y⁄, where x⁄ is the streamline coordinate and y⁄ is the distance in the perpendicular, shear, direction, extending from the centerline towards the wall, where h⁄ 6 y⁄ 6 h⁄. The linear momentum balance in the axial direction for steady, Poiseuille, flow is: 2

d ux dy

2

þ

dsxy 

dy

ð8Þ

;

where sxy is the extra stress tensor due to the viscoelastic contribution of the polymer in the flow, gs is the constant solvent viscosity,   p;x  dp =dx ; ux are the constant pressure gradient and velocity in the main flow direction, x⁄, respectively. Integrating Eq. (8) from the centerline, 0, to y⁄ leads to the following integral momentum balance: 

0 ¼ y p;x þ gs

dux   þ sxy : dy

ð9Þ ⁄

⁄

⁄

Integrating once more from y to the top channel wall, y = h , allows us to obtain an expression for the velocity profile as a function of y⁄:

ux ðy Þ ¼ 

p;x 2 1 ðh  y2 Þ þ  2gs gs

Z

h

y

sxy dy^ ;

ð10Þ

R þh ^ ¼ 0 and the symmetry of with the boundary constraint h sxy dy the velocity profile, ux ðy Þ ¼ ux ðy Þ, duly satisfied due to the antisymmetry of the shear stress, sxy ðy Þ ¼ sxy ðy Þ. In the Newtonian limit, sxy ¼ 0 and Eq. (10) reduces to the well-known Poiseuille profile. The bulk velocity can be found by inserting Eq. (10) in Eq. (2b):

ub ¼

2 h p;x  3 s

g

þ



sw   gs

   dux    ¼ p;x h :  þ sxy   dy  y ¼h

Combining Eqs. (11) and (12), we find

:

0 ¼ p;x þ gs

where the dimensional viscoelastic contribution is defined as R þh  ðlamÞ IV   12 h y sxy dy . Moreover, the magnitude of the wall shear  stress, sw , results by evaluating Eq. (9) at the top wall, y⁄ = h⁄:

1  2h gs

Z

þh

h



y sxy dy ¼

2 h p;x  3 s

g



ðlamÞ IV  h s

g

;

ð11Þ

sw ¼

ð12Þ

sw in terms of ub :

3gs ub 3IðlamÞ þ V 2 :  h h

ð13Þ

⁄  Using q u2 b ; ub and h as the characteristic scales for the pressure,  velocity, and distance, respectively, the relation sw ¼ p;x h becomes p,x = Cf/2. Furthermore, non-dimensionalizing the polymer  extra-stress by gp0 ub =h , the dimensionless velocity becomes:

ux ðyÞ ¼ C f

Reb0 1  b0 ð1  y2 Þ þ 4cb0 b0

Z

1

sxy0 dy:

ð14Þ

y

Integrating Eq. (14) with respect to y and using the total mass balR þ1 ance, i.e. 1 ux ðyÞdy ¼ 2, the skin friction coefficient is obtained:

Cf ¼

6cb0 6cð1  b0 Þ ðlamÞ þ IV0 ; Reb0 Reb0

ð15Þ

where

the dimensionless viscoelastic contribution is R þ1 ðlamÞ IV0   12 1 ysxy0 dy. For a pure Newtonian fluid b0 = 1, and Eq. (15) reduces to the known expression Cf = 6c/Reb0 for smooth channels. Note that here and in the following we keep c as a parameter because different values of it are used in the literature. For an Oldroyd-B constitutive model sxy0 ¼ dux =dy  u0x ðyÞ, ðlamÞ

which results in IV0 ¼ 1 and Eq. (15) gives Cf = 6c/Reb0 for any b0 and Wes0. Therefore, when the constitutive model does not exhibit shear thinning, the resulted skin friction coefficient remains at its Newtonian value. In this case, the velocity given by Eq. (14) reduces to the quadratic Poiseuille profile, ux = 3(1  y2)/2. On the other hand, when the constitutive model allows for a ðlamÞ shear-thinning effect, IV 0 is always less than unity. For instance, consider the steady, laminar flow of a viscoelastic fluid which obeys to the FENE-P constitutive model. Simplifying the original equations we end up to the following system of equations:

8 9   C f Reb0 0 > > y þ u0x ðyÞ b0 þ 1b ¼ 0 > > 2c f > > P > > < =  2 3 2 0 d fP  fP  2 Weux ðyÞ ¼ 0 >; > > > > > > > R1 : ; u ðyÞ dy ¼ 1 x 0

ð16Þ

d ¼ cWes0 =ðReb0 LÞ is a modified Weissenberg number, where We  Wes0  k u2 s =m0 is the zero shear-rate friction Weissenberg number, L is the maximum extensibility parameter of the polymer molecules, fP = fP(y) is known as the Peterlin function, and the viscoelastic extra-stress for this model is given by sxy0 ¼ u0x ðyÞ=fP . The first equation in (16) results from the momentum balance in the streamwise direction, the second one from the definition of the Peterlin function and the constitutive model, and the last is the total mass balance, with the symmetry with respect to the midplane having also been used. It is worth mentioning, as revealed by Eq. (16), that for the FENE-P model the independent parameters are CfReb0, b0 and d Therefore, there are no independent contributions of the fricWe. tion Weissenberg number, Wes0, and the maximum extensibility parameter, L. In the limiting case for which the fluid is pure viscoelastic, i.e. for b0 = 0, we can solve Eq. (16) analytically to obtain:

52

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

ðlamÞ

ðlamÞ

7 sinhð5x=6Þ þ 5 sinhð7x=6Þ ; pffiffiffi d3 2835 6 We

d 2 Þ: x  cosh ð1 þ 243 We 1

ðlamÞ

lam Iv 0 ðb0 Þ  ð1  b0 ÞIlam v 0 ðb0 ¼ 0Þ þ b0 Iv 0 ðb0 ¼ 1Þ

1.01

0

β0=0.1

β0=0.2

1.3

solid lines: channel flow dotted lines: pipe flow

β0=0.3 0

1

2

3

4

5

6

7

^

We = c Weτ0 / ( L Reb0 ) Fig. 2. The ratio of the skin friction coefficient for a viscoelastic fluid to the skin friction coefficient of a Newtonian fluid, in pipes and channels (laminar case) (a) for values of the viscosity ratio b0 close to unity, (b) for small values of the viscosity ratio b0.

deviations which grow as b0 decreases from b0 = 1 (the Newtonian limit) to b0 = 0 (for pure polymer). ðNewtÞ In Fig. 2, the ratio C f =C f ¼ C f Reb =ð6cÞ for channel flow, and

1.0

solid lines: pure polymer ( β0=0) dashed lines: infinitely dilution limit ( β0=1)

ðNewtÞ

the ratio C f =C f

¼ C f Reb =16 for pipe flow (see Appendix A) are

d for different b0 values. It is seen that offered as a function of We the ratio is always above unity with the maximum deviation increasing with decreasing b0 up to about 45% for a pure polymer (b0 = 0). However, for dilute polymer solutions, i.e for b0 > 0.9, Fig. 2a shows that the maximum deviation from unity is less than 1%. Therefore, under those conditions the bulk Reynolds number defined according to the wall viscosity approximately captures the shear-thinning effects and the Newtonian expressions, Cf = 6c/Reb for channel flow and Cf = 16/Reb, respectively, for pipe flow, can be utilized as rough first approximations for the skin friction coefficient.

0.6

Iv0

(lam)

5

1.4

R þ1

where b  gs =gw ¼ b0 =lw ; Reb ¼ Reb0 =lw , and   12 1 ysxy dy involves the polymer extra-stress sxy scaled by gpw ub =‘ , where gpw is the wall polymer viscosity. For a shear-thinning fluid, it is interesting to evaluate Cf based on the wall properties, following Eq. (20). Although in this case it turns out that Cf, given in terms of Reb, has a relationship close to the Newtonian case, Cf  6c/Reb, this is only exact in the Newtonian limit; there are considerable

4

β0=0

1.5

1.0

ðlamÞ IV

3

1.6

1.1

ð20Þ

2

We = c Weτ0 / ( L Reb0 )

For viscoelastic models which show a dependence of the shear viscosity on the shear rate Eq. (15) can also be expressed in terms of the wall properties:

6cb 6cð1  bÞ ðlamÞ Cf ¼ þ IV ; Reb Reb

1

^

1.2

channel

0.2

β0=0.9

1.00

ð19Þ

0.4

β0=0.8

ð18Þ

The results by expressions (17) and (18) are presented in Fig. 1, along with the corresponding results for pipe flow (see Eqs. (A.10) ðlamÞ and (A.11) in Appendix A). It is seen that Iv 0 ðb0 ¼ 0Þ and ðlamÞ d They Iv 0 ðb0 ¼ 1Þ decrease monotonically with increasing We. are also very similar in magnitude which means that for constant d the dependence of IðlamÞ on b0 is very weak. Indeed, although We v0 the velocity profile undergoes considerable changes from plug flow-like at b0 = 0 to fully parabolic at b0 = 1, the contribution of the extra-stress to the skin friction coefficient remains approximately constant. It appears that as the shear thinning increases and the velocity profile becomes steeper near the wall, the enhanced stress contributions are compensated by much smaller contributions near the centerline where the velocity profile flattens. Thus, for b0 – 0,1 and for all practical purposes, a linear interpolation between the two limiting expressions can be used:

0.8

dotted lines: pipe flow

1.02

(Newt)

In the infinitely dilution limit, for which R þ1 ðlamÞ b0 ! 1; Iv 0 ¼  12 1 ysxy0 dy can also be evaluated analytically:

Iv 0 ðb0 ¼ 1Þ ¼

β0=0.7

ð17Þ

:

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 d 27 We d þ 10 þ 729 We d2 22=3 9 We

solid lines: channel flow

1.03

Cf / C f

¼

( ) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2=3 1=3 d d 2 27 We þ 10 þ 729 We  10

(Newt)

5

1=3

Cf / Cf

Iv 0 ðb0 ¼ 0Þ

pipe

0.0 0

1

2

3

4

5

6

^ We = c Weτ0 / ( L Reb0 ) Fig. 1. The viscoelastic contribution to the skin friction coefficient for laminar flow in channels and pipes.

2.2. Stationary turbulent viscoelastic flow in a channel For fully developed turbulent flow in a channel the Reynolds averaged governing equations are derived by first decomposing each dependent variable, f⁄, as follows:

53

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

f  ¼ f  ðy Þ þ f 0 ðx ; y ; z ; t Þ;

ð21Þ

where the overbar denotes time and spatial (in all the periodic directions) averaging and, by definition, f 0 ¼ 0. Because of the continuity equation, it is easily shown that uy ¼ 0, i.e. the average velocity component in the wall-normal direction is zero. The nontrivial component of the momentum balance is:

q

2   ds  d    d u ;x þ gs 2x þ xy  ux uy ¼ p  ; dy dy dy

ð22Þ

Integrating twice Eq. (22) with respect to y⁄ we get the velocity profile in the main flow direction:

 x ðy Þ ¼ u

;x 2 p 1 ðh  y2 Þ þ  2gs gs

Z

h

y

sxy dy^ 

q gs

Z

h

y

^ : ux uy dy

ð23Þ

Eq. (23) can be written in dimensionless form by using the characteristic scales reported for the laminar flow:

ux ðyÞ ¼ ðp;x Þ 

Z

Reb0 1  b0 ð1  y2 Þ þ 2cb0 b0

Z

1

sxy0 dy^ 

y

Reb0 cb0 ð24Þ

y

Then, Eq. (24) is integrated with respect to y between 1 and 1, and R þ1 using 1 ux dy ¼ 2; C f ¼ p;x =2, and rearranging the result gives the final expression for the skin-friction factor:

Cf ¼

6cb0 6cð1  b0 Þ þ IV0 þ 3IR ; Reb0 Reb0

ð25Þ

R þ1 R þ1 xy0 y dy and IR ¼  1 where IV0 ¼  12 1 s ðux uy Þy dy. Note that we can also reformulate Eq. (25), by analogy to Eq. (15), using the wall viscosity as reference viscosity simply by replacing b0, Reb0, IV0 with b, Reb, IV:

Cf ¼

6cb 6cð1  bÞ þ IV þ 3IR : Reb Reb

ð26Þ

This procedure is similar to that used by Fukagata et al. (2002) for Newtonian wall bounded flows and by Yu and Kawaguchi (2003); note that these authors use explicitly c = 2. Eq. (26) shows that the total skin friction coefficient is decomposed in three individual contributions: a laminar (or viscous) contribution that is inversely proportional to the zero shear-rate bulk Reynolds number, a Reynolds turbulent contribution, 3IR, which is present for both Newtonian and viscoelastic flows, and a purely viscoelastic turbulent contribution, which is proportional to the quantity IV. Note that the expressions (25) and (26) have been derived by using as characteristic velocity the bulk velocity ub , which is a known quantity for experiments conducted with constant flux. When the experiments are conducted with constant pressure gradient, the suitable characteristic velocity is the friction velocity, us , Eq. (3), while the characteristic length is the same as before. In this case, the corresponding expression for the skin friction coefficient in channel flow is:

sffiffiffiffiffi   2 Res0 3 3ð1  b0 Þ ¼ 1  IR;s  IV0;s ; Cf 2 Res0 3b0

ð27Þ

where IV0,s and IR,s have been made dimensionless using us as the characteristic velocity. They are connected with IV0 and IR as:

Reb0 IV0;s ¼ IV0 ¼ cRes0 IR;s ¼

Re2b0 c2 Re2s0

IR ¼

sffiffiffiffiffi 2 IV0 Cf

2 IR : Cf

3. Analysis of the skin friction coefficient using DNS results Based on Eq. (26) we have calculated the three contributions to the skin friction coefficient in turbulent channel flow using previous DNS results for selected cases (Housiadas and Beris, 2003, 2004a; Housiadas et al., 2005). The corresponding flow and rheological parameters are reported in Table 1. In addition to nine viscoelastic cases, we have also included for comparison the results of four Newtonian simulations. In the third to the last column in Table 1, we present the zero shear-rate elasticity parameter which is defined as follows:

El0 

1

^: ux uy dy

drag. They also explicitly state that the most important quantities for accurate predictions of the skin friction coefficients and the resulted drag reduction are specific weighted averages of the Reyxy , respectively. nolds stress and the viscoelastic stress, ux uy and s As such, they can be readily calculated from the existing simulations data.

ð28aÞ ð28bÞ

The above expressions are useful because they analyze directly the effect of the Reynolds and the viscoelastic stress on the frictional

k m0 h

2

¼

Wes0

ð29Þ

Re2s0

The advantage of using El0, instead of Wes0, is that it does not involve the flow velocity, and therefore, is a material parameter, albeit not pure as it also involves the flow geometry, h⁄. Physically it represents the dimensionless ratio between the zero shear-rate kinematic viscosity and the ratio of the square of the channel half width to the polymeric relaxation time. In the second to the last column, we give the maximum Trouton ratio which is defined as:

ge ; We!1 g 0

Trm  lim

ð30Þ

where We  e_  k with e_  being the constant extensional rate for which the extensional viscosity ge is measured. The maximum Trouton ratio Trm depends physically on the polymer concentration, which for the models used here is proportional to 1  b0, and the molecular characteristics. It actually describes the maximum degree of extensional thickening, i.e. the increased resistance to the flow observed in uniaxial extensional flow. For a Newtonian fluid Trm = 3, for the FENE-P model Trm = 3b0 + 2L2(1  b0) where L is the molecular maximum extensibility parameter, for the Giesekus model Trm = 3b0 + 2(1  b0)/aG where aG is the Giesekus mobility parameter, while it is unbounded for the Oldroyd-B model. Trm, as independent of the flow, is a pure material parameter. Finally, in the last column of Table 1, we offer the dimensionless wall viscosity lw so that one can evaluate the dimensionless groups in terms of the wall properties; see the discussion below Eq. (20). For an alternative estimation of lw see the subsequent section. For completion, we also mention here (not shown in Table 1, though) about the ratio of the second to the first normal stress difference, N2/N1. For simple shear flow, this is zero for the FENE-P and Oldroyd-B fluid models, whereas it is equal to aG/2 for the Giesekus model. In Table 2 we show the skin friction coefficient as well as its analysis in the three contributions discussed in Section 2.2; see Eq. (25). We report the skin friction coefficient evaluated in two different ways: (a) according to the averaged equation, Eq. (25), and (b) according to its definition Eq. (7b). Comparison between these two values of Cf, shows agreement within 2% demonstrating the consistency of the calculations. Also, in Table 2, we show the drag reduction which is obtained from:

DR ¼ 1 

1

ðv iscÞ

Reb

l2w ReðNewtÞ b

!2n

ðv iscÞ

¼1 Res0

Reb

!2n ;

ðNewtÞ

Reb

Res

ð31Þ

54

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

Table 1 Flow and rheological parameters for the selected DNS cases used in the present work. Case

MODEL

Res0

Reb0/c

Wes0

b0

El0  104

Trm

lw

1 2 3 4 5 6 7 8 9 10 11 12 13

FENE-P, L = 30 FENE-P, L = 30 FENE-P, L = 30 FENE-P, L = 30 FENE-P, L = 30 FENE-P, L = 30 GIESEKUS, aG = 1/900 GIESEKUS, aG = 1/900 OLDROYD-B NEWTONIAN NEWTONIAN NEWTONIAN NEWTONIAN

125 180 395 590 590 590 182.7 182.7 180 125 180 395 590

2270 3450 8500 12150 13500 14470 3730 4270 3480 1840 2790 6890 10905

50 50 50 20 50 100 50 125 25 0 0 0 0

0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 1.0 1.0 1.0 1.0

32.0 15.4 3.21 0.575 1.44 2.87 15.4 38.6 7.72 0 0 0 0

182.7 182.7 182.7 182.7 182.7 182.7 182.7 182.7 1 3 3 3 3

0.944 0.944 0.944 0.969 0.944 0.930 0.932 0.911 1.0 1.0 1.0 1.0 1.0

Table 2 The skin friction coefficient and its contributions according to Eq. (25) for the cases reported in Table 1. In parenthesis the percentage of each contribution is given. Case

DR%

Cf  103 Eq. (25)

Cf  103 Eq. (7b)

Contribution 6cb0 Reb0

1 2 3 4 5 6 7 8 9 10 11 12 13

29.8 29.9 29.6 16.3 29.8 37.6 38.6 51.2 32 0 0 0 0

where n  1.14775 and

6.23 5.53 4.35 5.00 3.75 3.49 4.59 3.39 5.47 9.16 8.08 6.51 5.80

6.07 5.45 4.32 4.72 3.84 3.38 4.66 3.56 5.35 9.23 8.32 6.57 5.85



0  lw ¼ b0 = 1 þ 1b s Res0 xy



following and

extending (in the second equality) the theoretical development by Housiadas and Beris (2004b). The first equality in Eq. (31) is used when the numerical simulations are performed under constant pressure drop conditions for given zero shear-rate friction Reynolds number, i.e. for constant Res0, for both Newtonian and viscoelastic fluids. In this case, the achieved drag reduction is evaluated indirectly through the increase in the mean velocity of the viscoelastic fluid, which results in different bulk Reynolds numbers between the Newtonian and the viscoelastic flow. The second equality is used in Eq. (35) below when we derive an expression for the bulk Reynolds number, Reb, involving the friction Reynolds number, Res. From an analysis of the DNS results given in Table 2, several important conclusions can be readily extracted: (a) First we mention that the viscous contribution is given explicitly, inversely proportional to the bulk Reynolds number. (b) As far as the Newtonian cases are concerned (cases 10–13), the inertial contribution shows an interesting behavior with respect to an increasing bulk Reynolds number. It initially increases and then decreases. This behavior of the inertial contribution is encountered in the viscoelastic cases too, providing that the other rheological parameters remain fixed (cases 1, 2, 3, and 5). In contrast, there, the viscoelastic contribution 6c(1  b0)IV0/Reb0 decreases monotonically with increasing friction Reynolds number. The combination of these two effects leads to an approximately constant drag reduction as the Reynolds number increases as already noted in a previous publication (Housiadas and Beris,

 103

2.38 (38%) 1.57 (28%) 0.635 (15%) 0.444 (9%) 0.401 (11%) 0.376 (10%) 1.45 (31%) 1.23 (37%) 1.55 (28%) 3.26 (36%) 2.15 (26%) 0.871 (13%) 0.550 (9.5%)

3IR  103

6cð1b0 Þ Reb0 IV0

3.22 3.40 3.36 4.41 3.10 2.78 2.52 1.28 3.29 5.89 5.96 5.64 5.25

0.62 (10%) 0.56 (10%) 0.35 (8%) 0.15 (3%) 0.25 (7%) 0.34 (10%) 0.63 (14%) 0.78 (23%) 0.63 (11%) 0 (0%) 0 (0%) 0 (0%) 0 (0%)

(52%) (62%) (77%) (88%) (82%) (80%) (55%) (40%) (61%) (64%) (74%) (87%) (90.5%)

 103

IV0 2.35 3.22 4.96 3.04 5.63 8.20 3.92 5.55 3.65 0 0 0 0

2003). Percentage-wise, the inertial contribution in both the Newtonian and the viscoelastic cases monotonically increases with increasing Reb0 eventually dominating the Cf, albeit in the viscoelastic case is always less than the Newtonian. (c) Keeping the friction Reynolds number constant and increasing the friction Weissenberg number (cases 4, 5 and 6 for the FENE-P model, and cases 7 and 8 for the Giesekus model), the inertial contribution decreases whereas the viscoelastic contribution increases. This appears to hold independently of the constitutive model. However, the sum of these two decreases with increasing Wes0, consistent with the increasing drag reduction. The Wes0 number, along with the maximum Trouton ratio, is the most critical parameters for drag reduction, as also mentioned in our previous works (Housiadas and Beris, 2003, 2004a, 2005; Beris and Housiadas, 2010). It should be mentioned here that Tsukahara et al. (2011) observed a percentage increase of the viscoelastic contribution, as we did, as well, however, they saw that this contribution slightly decreases in absolute value with increasing Weissenberg number; instead we observed a slight increase. (d) The rheological parameters of the FENE-P and Giesekus models are so chosen so that the extensional characteristics, represented by the maximum Trouton ratio, shown in Table 1, are the same. In consequence, in cases 2 and 7, where both the friction Reynolds and Weissenberg numbers are also the same, the bulk Reynolds number is similar (and almost the same drag reduction is obtained as also noted by Housiadas and Beris, 2004a) resulting to similar contributions to the skin friction coefficient. There are some

55

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

secondary effects with the Giesekus model being slightly more drag-reducing, which have been attributed to the non-zero second normal stress difference in simple shear flow. These effects lead to a slightly higher bulk Reynolds number which also translates to a slightly reduced inertia contribution and a slightly increased viscoelastic one. On the other hand, comparing cases 4 and 9 which correspond to similar friction Weissenberg numbers (20 and 25, respectively) but with very different Trouton ratios (180 and infinity, respectively) one sees that the drag reduction has been doubled in the second case. (e) The Giesekus model gives the lowest skin friction coefficient (highest drag reduction) in case 6. This corresponds to the lowest inertial contribution and the maximum viscoelastic contribution. However, the decrease of the Reynolds contribution is higher than the increase of the viscoelastic contribution and this is why a larger drag reduction is observed. (f) The quantity IV0 is always larger than the laminar limit, i.e. ðlamÞ IV 0 > IV0 ; in fact IV0 > 1 in all viscoelastic cases considered here. One is reminded (see Section 2.1) that in the laminar ðlamÞ case IV0  IV 0 6 1 with the equality holding only in the case of a non-shear-thinning fluid. In fact, it is characteristic to the shear-thinning behavior of the models used here that ðlamÞ d reaching zero as We d!1 IV 0 decreases steeply with We (see Fig. 2). In contrast, for all the cases considered here under turbulent conditions, IV0 > 1, increasing significantly with increasing either Reynolds or Weissenberg number. This appears to be related to the extensional thickening behavior associated with these models as well as the extensional character of turbulent flow. In Table 3 we show the inverse of the square root of the skin friction coefficient as well as its three contributions discussed in Section 2.2; see Eq. (27). In parallel to Table 2, we report the skin friction coefficient evaluated in two different ways: (a) according to the averaged equation, Eq. (27), and (b) according to its definition Eq. (7b). Again, comparison between these two values of Cf, shows agreement demonstrating the consistency of the calculations. The last two columns in Table 3 can be obtained dividing the third and the second to the last columns of Table 2 by the skin friction coefficient which is evaluated using Eq. (25). Due to this dependency, the information contained in Table 3 is not really new. However, it may be helpful in interpreting the results. For example, because of the nature of Eq. (27), it is obvious that the sum of the two last columns in Table 3, is strictly bounded from above by unity. Their magnitudes are easier to compare against the pure viscous contribution which, as shown by Eq. (27), is unity. It is worth

Table 3 The inverse square root of the skin friction coefficient and its contributions according to Eq. (30) for the cases reported in Table 1. pffiffiffiffiffi C f Eq. (27)

Case

DR%

1=

1 2 3 4 5 6 7 8 9 10 11 12 13

29.8 29.9 29.6 16.3 29.8 37.6 38.6 51.2 32 0 0 0 0

12.55859 13.38347 15.21979 13.59741 16.48171 16.38135 14.78918 18.49465 13.35791 10.51783 11.13168 12.44226 13.1871

1=

pffiffiffiffiffi C f Eq. (7b)

12.6694 13.44737 15.16196 14.14214 16.32993 16.92728 14.76025 17.17513 13.52092 10.44846 11.12485 12.39394 13.13064

3/2IR,s

3ð1b0 ÞIV 0;s Res0

0.51685 0.61483 0.77241 0.882 0.82667 0.79656 0.54902 0.37758 0.60146 0.64301 0.73762 0.86636 0.90517

0.09952 0.10127 0.08046 0.03 0.06667 0.09742 0.13725 0.23009 0.11517 0 0 0 0

noticing that some of the effects described above, are clearer here. For instance, the inertial contribution in the Newtonian case monotonically increases with increasing friction Reynolds number. Last, we should mention that Eq. (27) and the data presented in Table 3, provide an alternative of representing friction data, to that of Eq. (25), and the corresponding data in Table 2. The former simply use the friction Reynolds number while the latter the bulk Reynolds number as the key flow parameter. Both have been used in the literature, and are reported here for completion. In conclusion we see that when viscoelasticity becomes important all contributions to the skin friction coefficient, viscous, inertia and viscoelastic contributions, are significant. Therefore, in order to be more quantitative on the skin friction coefficient, as well as to allow a comparison against experimental data, we need to use approximate relations for the various dimensionless parameters. This is the topic of the next section. 4. Development of a parametric relation for the skin friction coefficient In order to reconcile theoretical predictions to experimental data, we develop here an approximate relation for the skin friction coefficient for viscoelastic turbulent channel flow. This relationship is most naturally developed in Prandtl–von Karman coordinates, i.e. for the inverse of the square root of the skin friction coefficient with respect to the friction Reynolds number. As a motivation for this analysis, we reproduce in Fig. 3 the data offered by Paterson and Abernathy (1970). They have been collected on various concentrations of Polyox WSR-301 aqueous solutions in 0.248 inches internal diameter pipe and they are reported as Moody friction factor vs a bulk Reynolds number calculated based on the diameter of the pipe. Those data span different concentrations for the same system in the same pipe, thus better enabling us to infer trends among them. Note that similar data on turbulent pipe flow have also been collected for a number of different polymer solutions by other investigators; see for instance, Escudier et al. (1999), and especially Shetty and Solomon (2009), who also reported data for various concentrations of Polyox WSR-301 aqueous solutions in pipe flow, also reproduced here in the insert in Fig. 4. However, the effect of viscoelasticity is most clearly seen in the Prandtl–von Karman coordinates as shown in the main Fig. 4, also taken from Shetty and Solomon. To enable comparison to the data we first need to connect the material parameters that characterize the experimental system to the dimensionless parameters used in the simulation. In the process, we have to realize the limitations: The original experimental data were obtained with a broad molecular weight system with no detailed modeling of its rheology from which the relevant relaxation time and extensional behavior could be extracted. The situation is better with the latest data of Shetty and Solomon (2009) albeit, still there, there is no sufficient information to a priori fit all the model parameters. Thus, at this stage, the best we can offer is a comparison against the parametric predictions obtained numerically as we vary the values of key characteristic material dimensionless groups. The most pertinent of these is a dimensionless relaxation time characterizing memory effects within the material, the zero shear-rate elasticity parameter:

El0 ¼ k m0 =h

2

¼ Wes0 =Res20 ¼

1

lw

Wes =Re2s :

ð32Þ

An estimate of this parameter requires an evaluation of the effective relaxation time. Note that this may be different form the relaxation time governing the shear behavior at laminar conditions (Clasen et al., 2006). Furthermore, for associating polymers,

56

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

Fig. 3. The Darcy friction factor, f, vs the bulk mean Reynolds number, Re = Reb(c = 2), determined experimentally for various concentrations of Polyox WSR-301 aqueous solutions (reproduced by permission from Paterson and Abernathy, JFM, 43(4), 689–710, 1970).

to higher shear thinning), a proper evaluation of the wall viscosity is needed (Tsukahara et al. 2011). For the development of an approximate relation for the skin friction coefficient we use the insight gained from previous DNS work. One of the most important conclusions from our previous works (Housiadas and Beris, 2003, 2004a, 2005, 2006; Housiadas et al., 2005) and others (Yu and Kawaguchi, 2006)—see also the review (Beris and Housiadas, 2010)—is that the drag reduction remains unchanged when the friction Reynolds number is varied provided that Wes and all the other material parameters are kept the same. Indeed, as it turns out, all the previous results for DR, obtained by DNS, can be fit with respect to the friction Weissenberg number, Wes, simply by using a single new parameter, LDR. The latter is the asymptotic limit of drag reduction at large Weissenberg numbers, which is used to scale DR (see Fig. 5). A convenient form for this fitting is the following:

8 > < 0; DR ¼ 1 LDR > : Fig. 4. Prandtl–von Karman plots for different PEO WSR-301 concentrations ranging from 1 to 25 ppm. The lower curve is the Prandtl–von Karman (PK) law for a Newtonian solvent and the upper curve is the maximum drag reduction asymptote for polymer turbulent drag reduction. The inset plot is a Cf vs Reb (c = 2) plot for the different PEO WSR-301concentrations studied above (reproduced by permission from Shetty and Solomon, Polymer, 50, 261–270, 2009).

Wes < WeðonsetÞ s 1þexp



2

ðonsetÞ Wes We s DWes

ð33Þ

 ; Wes P WeðonsetÞ s

1.2

1.0

0.8

DR / LDR

filament formation and aggregation may be taking place as the polymer molecules are getting extended by the turbulent flow following experimental evidence (Baik et al., 2005). This may also affect the effective relaxation time of the material; see also the discussion at the end of this section. The second parameter that is proposed here is the limiting drag reduction, LDR, observed at high Weissenberg numbers. This is connected to the rheological parameters of the constitutive model that are responsible for the extensional behavior of the material. It can be determined a priori from the DNS results for a given rheological model, or it can be fit to experimental data a posteriori. The El0 and LDR parameters are anticipated to be the most critical material parameters for drag reduction based on all previous experimental and, more recently, numerical investigations (Housiadas and Beris, 2004a; Tsukahara et al., 2011). In addition, especially with more concentrated systems (and therefore more prone

0.6

Boltzman fitting curve Giesekus, α=1/900 Oldroyd-B FENE-P, L=30 FENE-P, L=20 FENE-P, L=10

0.4

0.2

0.0 0

20

40

60

80

100

120

140

160

180

200

Weτ Fig. 5. The universal fitting curve for the ratio of drag reduction along with the results obtained by DNS of viscoelastic turbulent channel flow.

57

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

where WeðonsetÞ represents the Weissenberg number for the onset of s drag reduction, WeðonsetÞ  6, and D Wes  25 as they arise from the s fitting of the existent DNS results. The fitting curve, along with the DNS results for the drag reduction, are shown in Fig. 5. The LDR parameter corresponding to the various constitutive models and their parameter values are reported in Table 4. Note that the LDR parameter in general increases with increasing Trouton ratio; however, other parameters may also affect LDR (contrast, for instance, columns 3 and 4 in Table 4). An explanation is warranted here regarding the data used in developing our universal drag reduction curve. Ideally, one would have liked to use all available literature data to generate/validate the universal fit proposed in Fig. 5, namely from DeAngelis et al. (2002), Min et al. (2003), Yu and Kawaguchi (2006), Dubief et al. (2004), Li et al. (2006a,b), Dallas et al. (2010) and others. However, quantitative details in turbulent friction predictions, on which the universal curve critically relies upon, are affected by a number of numerical issues, such as the computational domain size and resolution, the numerical method and numerical parameters, such as the artificial diffusion coefficient and the time of integration to collect the data. Different groups using different methods, and different values of all the above parameters, have generated different accuracy results that are not directly comparable. Thus, we have opted here to rely on our own data only, since, as all of them have been generated using the same method and most the same of the numerical parameters, they represent the most self-consistent results available over a wide range of material parameters that can lead to meaningful comparisons. Next, we use the relation

Reb ¼ cCRens ;

Parameters

Housiadas and Beris (2004b)

Dean (1978)

Pope (2000)

C n

7.2001 1.14775

7.320676 1.14286

7.71496 1.13636

in developing an expression for friction factor that can be used for a wide range of Reynolds number values for turbulent pipe flows a Colebrook-type relation is more appropriate as it naturally emerges from the integration of the mean velocity profile (Colebrook, 1937; Pope, 2000; Bird et al., 2002; Fox et al., 2004). More specifically, we propose using the equation developed in Section 5 as it arises from the integration of a mean velocity profile developed as an approximation of experimental data, following the suggestion of Bird et al. (2002). This relation, developed in Eq. (58), and in connection with the effect on the drag reduction seen in Eqs. (35) and (36) above, leads to the following general relationship applicable for viscoelastic pipe flow (c = 2): !  qffiffiffiffiffi 1 1 162:3 1586 pffiffiffiffiffi ¼ pffiffiffiffiffi þ 1:7678  ln Reb C f  0:60  : C f ð1  DRÞn~=2 Reb C f Re2b C f ð37aÞ

Rearranging Eq. (37a) gives: ! pffiffiffi  pffiffiffi  2 2Res 162:3 1586 p ffiffiffi Reb ¼ 1:7678  ln 2 2 Re þ ;  0:60  s ~ 2 2Res 8Re2s ð1  DRÞn=2

ð37bÞ

ð34Þ

which is established for turbulent channel flow in the Newtonian limit (Dean, 1978; Pope, 2000; Housiadas and Beris, 2004b), and the second definition for the drag reduction that appears in Eq. (34), in order to express the bulk mean Reynolds number in a viscoelastic turbulent channel flow as follows:

Reb ¼

Table 5 The parameters C and n, developed for turbulent channel flows, used in Eqs. (34)– (36).

cC ð1  DRÞ

n=2

Rens ;

ð35Þ

In the above equation, C and n are fitting parameters to experimental data or DNS Newtonian results; see Table 5 adapted from Housiadas and Beris (2004b). From Eq. (35), and the definition of the skin friction coefficient C f ¼ 2c2 Re2s =Re2b (see Section 2) we can also get the following relationship for the skin friction coefficient for turbulent channel flow:

1 CRen1 s pffiffiffiffiffi ¼ pffiffiffi : Cf 2ð1  DRÞn=2

ð36Þ

Similarly, for turbulent pipe flow, we can also use analogues of Eqs. (35) and (36) with simply the parameters C, n replaced by the e; n ~ applicable for pipe flow; see Seccorresponding parameters C tion 5 for possible sets of values. However, as also mentioned in Section 5, power-law Blasius-type relations developed for Newtonian turbulent flows, such as those described above for a given e; n ~ , have a relatively limited range of Reyset of coefficients C, n or C nolds number values validity. For Newtonian turbulent channel flow, as the application exploited here involved a limited range of Reynolds numbers, that was not that much of an issue. However,

where Eq. (62), further below, can be used to get a very good ~ which is a weak function of approximation for the coefficient n the bulk Reynolds number, giving (for c = 2)

~ ¼1þ n

1:085 6:535 : þ lnðReb Þ ðlnðReb ÞÞ2

ð38Þ

~  1:18 as an initial guess allows for a rapid convergence. Use of n In Eqs. (35)–(38) the relative wall viscosity lw appears implicitly in the definition of the friction and bulk Reynolds numbers as well as in the drag reduction, DR, since this is given by Eq. (33) as a function of the friction Weissenberg number that depends, from Eq. (32), on the zero shear-rate elasticity number, the friction Reynolds number and the relative wall viscosity. As shown in Housiadas and Beris (2003, 2004a,b) lw can be approxiðssÞ mated to a very good extent by the steady shear viscosity, lw , at the mean wall shear rate, rendered dimensionless with respect to its zero shear-rate value. This only depends on the friction Weissenberg number and the material rheological parameters of the model, but not on the bulk Reynolds number. For instance, for ðssÞ the FENE-P model lw can be found as:

lðssÞ w  b0 1  b0

¼

1 þ 2 cosh



3

1 1 cosh ð1 3

þ 27ðWes =LÞ2 Þ

:

ð39Þ

Note also that, in the limiting case of zero viscous solvent contribution, i.e. b0 = 0, Eq. (39) can also be recast into the much simpler form

Table 4 The asymptotic DR at large We, LDR, the maximum Trouton ratio, TRm, and the mc Carreau parameter for the various constitutive models used here. Model

FENE-P (L = 10, b0 = 0.9)

FENE-P (L = 20, b0 = 0.9)

FENE-P (L = 30, b0 = 0.9)

Giesekus (aG = 1/900, b0 = 0.9)

Oldroyd-B (b0 = 0.9)

LDR TRm mc

0.15 22.7 0.472

0.32 82.7 0.472

0.40 182.7 0.47

0.51 182.7 0.212

0.64 1 1

58

lðssÞ w ¼

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

1 1 þ 2ðWes0 =LÞ2

and, therefore, the simulations cannot predict this kind of behavior (which is beyond the scope of the current work). In the following Figs. 7–10, we switch from pffiffiffiffiffi Cf vs Reb plots, pffiffiffiffiffito Prandtl–von Karman plots, p i.e. in plots 1= C vs log Cf , Re f b   10 pffiffiffi ffiffiffiffiffi noting that the quantity Reb C f ¼ c 2Res is essentially the friction Reynolds number. This is a more suitable way to discuss and interpret drag reduction phenomena as the differences in friction factor values are much more accentuated. For comparison purposes, the predictions corresponding to the maximum drag reduction (Virk) asymptote (Virk, 1975; Gyr and Bewersdorff, 1995) are also shown in those figures as derived based on the integration of the best approximation of experimental data available for this type of flow described by Eq. (66) in Section 5. First, in Fig. 7 we examine the effect of the shear-thinning behavior between three different cases: a non-shear thinning case, obtained in the limit of a very dilute polymer solution (i.e. for b0 = 1), a regular dilute solution, corresponding to b0 = 0.9, value that is most often used in DNS simulations of viscoelastic turbulent flow, such as the cases reported in this work, and for illustration of an extreme shear thinning case, a pure polymer case (i.e. for b0 = 0). All the other material parameters being the same (such as the zero shear-rate Elasticity, El0 and the Limiting Drag Reduction, LDR), the observed differences are solely due to differences in the relative wall viscosity lw. From Fig. 7,pitffiffiffiffiffiis clear pffiffiffithat for constant friction Reynolds number values p Reffiffiffiffiffi C f ¼ c 2Res , shear-thinning reb sults to lower values of 1= C f , and that the curves shift away from the MDR asymptote. However, in order to see a significant effect, the shear thinning has to be quite pronounced, certainly more than what corresponds to typical dilute polymer solution solvent viscosity ratio values (i.e. for b0 = 0.9) which case produces almost indistinguishable results from the non-shear-thinning limiting case of (i.e. for b0 = 1). It is also seen that the effect of shear thinning is smaller at larger values of the LDR parameter. Nevertheless, even in cases where shear thinning is important to quantitatively affect the results, their qualitative character of the predictions is not altered. In Fig. 8, the effect of the LDR parameter is presented in a Prandtl–von Karman plot, for fixed elasticity number, El0 = 8  104. The wall viscosity has been calculated using the FENE-P constitutive model with L = 30, b0 = 0.9. It is observed that when the LDR increases the curves move closer to the MDR asymptote, as expected. More specifically, we observe that the initial slope of the curves, as this develops away from the Newtonian limit, is the

ð40Þ

;

which is explicit in terms of the zero shear-rate friction Weissenberg number, Wes0; one is reminded here that Wes0 = Weslw. ðssÞ Similarly, for the Giesekus model lw is a function of Wes and the mobility parameter aG. Note that the friction Weissenberg number Wes is directly related to the zero shear-rate elasticity parameter, El0, as given by Eq. (32). The procedure to generate predictions for the skin friction coefficient for given flow, fluid and geometry has as follows. For a given fluid and geometry, i.e. for fixed material/geometry parameters El0, b0 and L (for FENE-P), or aG (for Giesekus), and for a given flow as ðssÞ described by Wes, one first determines lw (using, for example, Eq. (39) for the FENE-P model). Then, one determines Res from Eq. (32), and DR from Eq. (33). Finally, from DR and Res we find Reb, through Eq. (36), and from there, the skin friction coefficient as C f ¼ 2c2 Re2s =Re2b . In the case of pipe flow, this process does involve ~ is a function of the bulk an iterative procedure as the exponent n Reynolds number, as Eq. (38) shows. However, as mentioned above, this dependence is very week, and for all practical purposes ~  1:18 represents a very good guess that will require the value n only one or two iterations for convergence. The results for Cf for turbulent channel flow are shown in Fig. 6,  as a function of the bulk Reynolds number Reb  2h ub =mw , for various constitutive models. It is clear that the behavior seen in Fig. 3, i.e. for the experimental data, is also observed in Fig. 6, i.e. for the theoretical predictions. As the concentration and, correspondingly, the maximum extensional viscosity and therefore the LDR increase, we do see a shifting of the curves from the Newtonian to the MDR limit; for the highest concentration (highest extensional effect in the theoretical predictions) we see an approach to the MDR, which shows that the theoretical data are in the proper range of parameters to fit experimental results. On the other hand, the elasticity parameter, which depends primarily on the molecular characteristics of the polymer, controls the onset of drag reduction; the higher El0 the earlier the onset. In comparison with the experimental data shown in Fig. 3, we need to mention here that the non-monotonic behavior of those experimental curves as a function of the Reynolds number is to be attributed to the degradation of the macromolecules, which usually takes place at high Reynolds numbers. A mechanism for degradation has not been included in the equations

0.01 0.009 0.008 0.007

-4

El0=2x10 Newtonian -4

El0=0.5x10

-4

El0=8x10

Newtonian (pipes) -4

El0=0.08x10

0.006 0.005

Cf

0.004 0.003

Laminar Cf=12/Reb

0.002 MDR (pipes) Solid lines (+): FENE-P, L=10 Solid lines (.): FENE-P, L=20 Solid lines (x): FENE-P, L=30

squares: Oldroyd-B triangles: Giesekus, α=1/900

1E-3 4

5

10

10

Reb Fig. 6. The skin friction coefficient (or Fanning friction factor) as a function of the bulk Reynolds number, in a channel; c = 2.

59

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

25

β0=0

MDR 20

Laminar

1/Cf

0.5

LDR =0.7 β0=0

15

LDR =0.4

10

Newtonian -4

El0=8x10

5 2.2

2.4

2.6

2.8

3.0

Fig. 9. Schematic of types A and B drag reduction (reproduced by permission from Gyr and Bewersdorff (1995)).

3.2

3.4

0.5 log10 (RebCf )

25 Fig. 7. The effect of shear thinning in a Prandtl–von Karman plot for pipe flow at LDR = 0.4 and 0.7, where Reb is calculated with c = 2. The lines with dots correspond to constant viscosity fluids (lw = 1). The shear thinning is evaluated with the FENEP model using L = 30 and b0 = 0.9 (lines with dashes). Results for the FENE-P model with L = 30, b0 = 0, i.e. for pure polymer, are also shown (solid lines). The laminar case, Newtonian turbulence and maximum drag reduction are included for comparison.

0.5

1/Cf

0.7

MDR

0.5

1.8x10

Laminar

-4

0.8x10 0.2x10

15

0.6

Laminar

Newtonian 10

0.5

1/Cf

-4

20

-4

25

20

-4

5.0x10

MDR

0.4

LDR=0.6

15 5 2.2

2.4

2.6

2.8

3.0

3.2

3.4

0.5

log10 (RebCf ) 10

Fig. 10. The effect of the zero shear-rate Elasticity number, El0, in a Prandtl–von Karman plot for pipe flow. Results are shown for El0  104 = 0.2,0.8,1.8 and 5. Reb is calculated with c = 2 and the shear thinning is evaluated with the FENE-P model using L = 30 and b0 = 0.9. Results for the laminar case, Newtonian turbulence and maximum drag reduction are also included for comparison.

Newtonian -4

El0=8x10

5 2.2

2.4

2.6

2.8

3.0

3.2

3.4

0.5 log10 (RebCf )

Fig. 8. The effect of the limiting drag reduction parameter, LDR, in a Prandtl–von Karman plot for pipe flow. Results are shown for LDR = 0.4, 0.5, 0.6 and 0.7. Reb is calculated with c = 2 and the shear thinning is evaluated with the FENE-P model using L = 30, b0 = 0.9. Results for the laminar case, Newtonian turbulence, and maximum drag reduction are also included for comparison.

one which most sensitively depends on LDR. The gradually increasing slope, with increasing LDR (which for a given polymer additive increases along with the concentration) while the onset remains the same, is also tantalizing close to the experimental findings of Shetty and Solomon (2009)) obtained with different PEO WSR301 concentrations—see Fig. 4, taken from that work.

With increasing LDR the slope increases, while the subsequent transition to a secondary, less steep, slope occurs at higher pffiffiffiffiffi 1= C f values. This is tantalizing close to the mechanisms A and B of drag reduction first suggested by Virk (1975b) and reported more extensively later by Gyr and Bewersdorff (1995); see Fig. 9 with A corresponding to the primary, higher, slope and B corresponding to the secondary, smaller, slope of the curves. In fact, we conjecture from the data that the MDR curve represents simply the loci where the transition of the primary to the secondary slope is taking place for high LDR values. It is also instructive based on the modeling of drag reduction, Eq. (36), to correlate type A behavior with partially extended molecules (the steep slope in the DR curve) and type B with fully extended ones (the upper saturated

60

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

part of the curve). Therefore our theoretical predictions seem to be in full agreement with both the letter, i.e. regarding the steepness of the slope of the curve, as well as the spirit, i.e. regarding the underlying physical reasoning, of the proposed A and B mechanisms. The A mechanism, corresponding to the initial slope, is associated with a partial extension of the molecules whereas the B mechanism, corresponding to the final slope, with fully extended ones, exactly as proposed by Virk (1975b), Virk and Wagger (1990), and Virk et al. (1997). In Fig. 10, the effect of El0 in a Prandtl–von Karman plot is shown. It is seen that as El0 decreases, the predictions develop a longer Newtonian tail with the results departing from the Newtonian asymptote at larger Reynolds numbers. The departing point corresponds to the drag reduction onset which takes place at higher Reynolds numbers for smaller elasticity values. Inversely, as El0 increases significant departures from the Newtonian curve are observed at lower Reynolds numbers, with the curves shifting towards the MDR asymptote. It is also interesting that the maximum slope of the El0 = 5  104 curve is even higher than that of the MDR asymptote in agreement with the trends schematically shown for curves 1 and 3 in Fig. 9. Note that the same holds for even higher values of El0. The elasticity parameter, therefore, controls the onset of drag reduction with respect to the bulk mean Reynolds number. Regarding the evaluation of the wall viscosity ratio, which is important in order to assess the shear-thinning effects, we note here that either we have to solve the corresponding constitutive equations as described above, or, alternatively, we can use experimental viscometric data. The latter can be conveniently represented by a Carreau viscosity model (Bird et al., 1987):

n



gp ¼ gp0 1 þ kCarreau c_ w

2 omC21

;

ð41Þ

where kCarreau is a fitting parameter with dimensions of time, and mc is a dimensionless exponent. This parallels the Cross model fitting, previously used by Escudier et al. (1999) to model the viscosity of a polymer solution, studied in turbulent pipe flow. In dimensionless form Eq. (41) becomes:

lðssÞ w  b0 1  b0

mC 1 2

¼ f1 þ ðXWes Þ2 g

;

ð42Þ

where X is the ratio of kCarreau to the relaxation time of the model k⁄. Note that X is not necessarily close to unity. For example, for the FENE-P and Giesekus models, very good Carreau viscosity fits can pffiffiffi be obtained using X ¼ kCarreau =kFENEP ¼ 7=L and X ¼ kCarreau = ffiffiffiffiffi p kGiesekus ¼ 2 aG , and the exponents mc shown in Table 4. In the last two expressions for X; kFENEP and kGiesekus are the relaxation times for the FENE-P and Giesekus constitutive models, respectively, that are used to fit the available experimental data with these models. The fact that X is generally much different than unity (and in particular X  1 under conditions of high Trouton ratios) brings forward the danger of estimating the dominant polymer relaxation time from simple shear data through a Carreau fitting. Moreover, as data extracted for dilute polymer solutions in both shear and extensional transient flows seem to indicate (Clasen et al., 2006), it appears that the relaxation time extracted from simple shear viscometry may be significant less to that needed in extension. This may be further complicated if, as it is usual the case for associated polymers, we have the development and aggregations of filaments due to the turbulent flow (Baik et al., 2005) further increasing the effective relaxation time. For instance, as a further corroboration on this issue, we report the work by Shetty and Solomon (2009) who found that the critical Weissenberg number for the onset of drag reduction, based on shear relaxation time, is WeðonsetÞ  1:6 while from our DNS works we inferred here that s

WeðonsetÞ  6. Independently, Kalashnikov (1998) determined that s  5:6 for any viscoelastic flow near smooth walls coming WeðonsetÞ s much closer to our predictions. Similarly, Li et al. (2006a,b) reported that their simulations indicated that the onset friction Weissenberg number lies in the window 5 < WeðonsetÞ < 10. s Before concluding this section, it is worth mentioning that we can use empirical expressions for the mean velocity boundary layer profile, in order to evaluate approximately, either the Newtonian or the MDR skin friction coefficient data. This analysis is shown in the subsequent section. A side benefit of this capability is the potential to evaluate expressions for the inertia contribution to the skin friction coefficient, IR, in the limiting case of pure Newtonian turbulent flow and for the viscoelastic contribution to the skin friction coefficient, IV, of MDR flows. 5. The skin friction coefficient for turbulent pipe flow based on the mean velocity profile From Eq. (1) it is straightforward to evaluate the skin friction coefficient from expressions supplied for the mean velocity profile. It is instructive to pursue this correlation in order to explore the connection between the contributions that the skin friction coefficient consist of, such as those shown in Eq. (30), and popular approximate expressions (experimentally determined and theoretically justified) for the mean velocity profile. 5.1. The mean velocity profile for Newtonian turbulent flow It is well known that in a Newtonian turbulent pipe flow, the mean dimensionless velocity profile can be best approximately expressed as follows (Lin et al., 1953; Bird et al., 2002). Starting with the averaged axial momentum equation with respect to the two periodic spatial directions, h, z and time, t for a Newtonian system (see also Eq. (B.2)), we get:

q

   1 d       1 d  dux  : r u u þ g r ¼  p     ;x s  x r r  dr r dr dr

ð43Þ

Weighting the above equation by r⁄, integrating along the radial direction from 0 to r⁄, and dividing the final result by r⁄, gives

r 2

q ux ur ¼  p;x þ gs

 x du  : dr

ð44Þ

   Using q u2 s ; us , and m =us as characteristic scales for the pressure, velocity, and distance, respectively, Eq. (44) can be rewritten in dimensionless (plus) turbulent quantities as:

uþx uþr ¼ 

þ rþ þ du ;xþ þ þx : p 2 dr

ð45Þ

With those characteristic scales the relation sw ¼ p;x R =2, which results from the total force balance along the main flow direction, becomes p,x = 2/Res and Eq. (45) simplifies to

  þ yþ du þ þx ; uþx uþr ¼ 1  Res dr

ð46Þ

where we also used the following identity for the expressing the radial coordinate, r+ in terms of the distance from the wall, y+, defined as

yþ  Res  r þ :

ð47Þ

According to (Lin et al., 1953) the following approximation can be used in representing the eddy viscosity, e+, of the Reynolds stress contribution to the shear stress in the viscous sublayer next to the wall, i.e. for r⁄  R⁄:

61

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

eþ  uþx uþr =

 þ 3 þx du y ; þ ¼ 14:5 dr

0 6 yþ 6 5:

ð48Þ

When we substitute this expression in Eq. (46) the following expression can be obtained for the dimensionless (in plus units) derivative of the mean axial velocity with respect to y+:

  yþ 1  Re  þx  þx du du s  yþ 3 ; þ ¼  þ ¼ dy dr 1 þ 14:5

0 6 yþ 6 5:

ð49Þ

This expression can be integrated analytically from the wall (where due to no-slip conditions the velocity is zero), following [Lin et al., 1953], in order to arrive to the following closed form expression for the velocity in the viscous sublayer, close to the wall:

0 1 0     yþ 1 þ 14:5 14:5 14:5 B C pffiffiffi B 14:5  þx ¼ þ 3 1  ln @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u A @ 1þ  yþ 2 3 Res Res yþ þ 14:5 1  14:5 !! ! 2yþ 1 p ; 0 6 yþ 6 5: þ ð50Þ  tan1 14:5pffiffiffi 6 3 Alternatively, since the corrections from the straight profile are small, Eq. (50) can be used to evaluate the first and higher derivatives at the wall location. Those, together with the no-slip wall condition, can then be used to form a Taylor expansion from the wall for the velocity in the sublayer, following what was done in Eq. (5.3-13) in Bird et al. (2002):

 þx u

þ

¼y

 3 ! 1 þ 1 yþ 1 y  ; 2Res 4 14:5

0 6 yþ 6 5:

ð51Þ

In the viscous sublayer, 0 6 y+ 6 5, Eqs. (50) and (51) give results that differ at most 0.05%. Thus, Eq. (51) is preferred because is much simpler than Eq. (50). For the buffer layer the approximation suggested by Lin et al. (1953) is also used. Consequently, the overall velocity profile is given as

 þ ðyþ Þ ¼ u

 8  þ 3  > yþ y þ þ þ 1 > 0 6 yþ 6 5 > < u1 ðy Þ  y 1  2Res  4 14:5 uþ2 ðyþ Þ > > > : þ þ u3 ðy Þ

2

þ

þ

 j lnðy þ aÞ þ b

5 6 y 6 30

 j1 lnðyþ Þ þ B

yþ > 30

8   þ 3   þ 4  þ 6  > y y y 1 1 þ 1 : þ  3  j lnðy Þ þ B; yþ > yþc u þ

þ

ð54Þ

yþ c

where is determined so that the velocity is continuous and where, for better accuracy, due to the extended nature of the viscous sublayer, a higher order Taylor expansion was used to approxþ imate the exact solution provided by Eq. (50) for u 1 . Note that this expansion exhibits at most 2 parts per thousand error over that extended viscous sublayer region justifying its use in subsequent calculations. 5.2. Analysis of the skin friction coefficient in turbulent flow In the following, we develop an expression for Cf in turbulent pipe flows, corresponding to the mean velocity profiles given by Eq. (52) or (54). More specifically, Eq. (52), or (54), is used in conjunction with the definitions of the skin friction factor, Eq. (7), and the friction and bulk Reynolds numbers, Eqs. (5)–(7), in order to derive a relation between Cf and Reb, i.e. Cf = Cf(Reb), for Newtonian or maximum drag reduction viscoelastic turbulent pipe flows, respectively. For turbulent flow in a pipe, the ratio ub =us can be described as an integral over the distance from the wall y+, defined through Eq. (47), as:

Z

ub 2 ¼ us Re2s

Res

þ r þ dr þ ¼ u

0

Z

2 2

Res

Res

 þ ðRes  yþ Þ dyþ : u

ð55Þ

0

For a Newtonian fluid, substituting Eq. (52) in Eq. (55) gives:

ð52Þ

where j, and B are the law-of-the-wall turbulent layer constants,  þ  u   =us ; yþ  y = m =us and the coefficients a, b are so deteru mined as to ensure continuity of the velocity profile, with 5, 30, loosely determining the extent of the intermediate, buffer, sublayer, respectively. For example, the original pipe-flow measurements for Newtonian turbulence, by Nikuradse in 1932 suggested that j  0.4 and B  5.5 (for which a = 0.1934 and b = 2.9707, when Res = 180) while later data correlations by Coles and Hirst (1968) showed that j  0.41 and B  5.0 (B  5.2 is used by Pope, 2000, p. 291) are improved values (for which a = 0.4178 and b = 3.3631, when Res = 180). When, more recent, higher Reynolds numbers experimental data are used (Zaragola and Smits, 1997), a slightly better fit is obtained with j = 0.436 and B = 6.13 (Pope, 2000 p. 291). In general, the values of the coefficients a,b are given as a function of j, B as

j   6X ; X  exp uþ3 ð30Þ  uþ1 ð5Þ 2 X1 pffiffiffiffiffiffi ! 2 30 þ B: b ¼ ln j ð30 þ aÞ

It is of interest that a similar expression to the law of the wall has also been established by Virk (1975a) even for the viscoelastic case and actually characterizes the limit of ‘‘maximum drag reduction’’. The corresponding parameters are j  1/11.7 and B  17. In this case though, it is customary to use that expression all the way from its meeting of the viscous sublayer (Gyr and Bewersdorff, 1995) (which, in this case and for Res = 180, occurs at a critical plus distance yþ c  8:3). Thus, the mean velocity profile in the viscoelastic case can be described then as follows:

a¼5

ð53Þ

In the following, we use the values of the first set as we are interested to represent the friction factor behavior for primarily low to intermediate Reynolds numbers, where those values do a better job.

ub 2 ¼ 2 us Res Z þ

Z

5 0

Res

10

 þ1 ðRes  yþ Þ dyþ þ u

 þ3 ðRes  yþ Þ dyþ u

Z

10

5

 þ2 ðRes  yþ Þ dyþ u



0" #y¼5 1 2 @ y2 y3 y5 A   ¼ Res 2 6Res 20ð14:5Þ3 y¼0

0" #y¼5 1 2 @ y3 y4 y6 A   2  3 8Res 24ð14:5Þ3 Res y¼0

 y¼10 ! 2 a ðy þ aÞðlnðy þ aÞ  1Þ þ by þ 1þ Res Res j y¼5 !

 2 y¼10 2 1 1 y ðy þ aÞ2 ðlnðy þ aÞ  Þ þ b  2 2 j 2 y¼5 Res

y¼Res ! 2 1 þ yðlnðyÞ  1Þ þ By Res j y¼10

y¼Res ! 1 1 2 1 1 2  2 y ðlnðyÞ  Þ þ By  lnðRes Þ 2 j j Res y¼10   3 þ B þ X; 2j

 2

ð56Þ

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

where X is so defined above as to represent the viscous sublayer and buffer correction. Note that this correction is negligible at large Reynolds numbers being proportional on the first and higher powers of the inverse of the Reynolds number, as evaluated from the above equation, substituting for the parameters a, b with the expressions provided from Eq. (53). Eq. (56) with the aid of Eqs. (5)–(7) gives (taking into account that for pipe flow c = 2):

ð57Þ

For j = 0.4, B = 5.5, c = 2, Eq. (57) gives:

pffiffiffiffiffi! Reb C f 1 2:5 1:75 X pffiffiffi pffiffiffiffiffi ¼ pffiffiffi  ln þ pffiffiffi þ pffiffiffi Cf 2 2 2 2 2  qffiffiffiffiffi 162:3 1586 pffiffiffiffiffi þ ;  1:7678  ln Reb C f  0:60  Reb C f Re2b C f

ð61Þ

ð58Þ

where the parametric expression for the viscous and buffer sublayer corrections has been found using a linear regression of the exact value of X,pffiffiffibased pffiffiffiffiffi on Eq. (56), in terms of a polynomial of 1=Res  2 2= C f Reb for the given parameters j, B (c = 2), accurate within the range of interest to less than 0.05% relative error. Without the sublayer and buffer layer corrections, Eq. (58) is the Colebrook expression for the friction factor (Fox et al., 2004), first obtained by Colebrook (1939) using a similar approach and neglecting the viscous and buffer sublayer effects. However, as pointed out by Bird et al. (2002), if one neglects the viscous and buffer layer effects, and one integrates the law of the wall all the way from y+ = 1, one needs to ‘‘fudge’’ a p bitffiffiffi the parameters of the pffiffiffi so-derived Colebrook relation from 2:5= 2 to 2:45= 2 and from pffiffiffi pffiffiffi 1:75= 2 to 2:0= 2 to get better agreement with experimental friction data:

pffiffiffiffiffi! Reb C f 1 2:45 2:0 pffiffiffi pffiffiffiffiffi  pffiffiffi  ln þ pffiffiffi Cf 2 2 2 2  qffiffiffiffiffi  1:7324  ln Reb C f  0:387:

ð59Þ

However, it is of interest that there is a comment in the same reference (Bird et al., 2002, p. 167) suggesting that this fudging would probably be not necessary if the integration was performed correctly taking into account the full velocity information over all the layers. Interestingly enough this was not carried out in that reference. Indeed, when our predictions from Eq. (58) are compared against those of Eq. (59) they come very close, especially at higher Reynolds numbers, as shown in Fig. 11. This makes sense as the more rigorous expression derived here is expected to be of more value (and more accurate) at lower Reynolds numbers, where the importance of the viscous and buffers layers is higher (which is also the reason this formula is used here). A very similar expression to Eq. (59) with coefficients 1.737 and 0.40, respectively is also offered by Pope (2000) as Eq. (7.98). An alternative expression to Eq. (58) is also useful, even if more approximate, representing the friction factor Cf as a direct function of the bulk Reynolds number, Reb. This is the role of the power-law Blasius expression

e n () Reb ¼ c CRe s ~

¼ 2c

2ð~ n1Þ ~ n

~1Þ 2ðn ~ n

;

~ ¼1þ n

1:085 6:535 ; þ lnðReb Þ ðlnðReb ÞÞ2

ð62Þ

Eq. (62) was derived originally from experimental data by Zaragola et al. (1997). Indeed, when we used the predictions derived from the integration of the mean velocity profile, Eq. (58), to perform a linear regression within the range of interest for our work (3000 < Reb < 50000) we get as best fit, for c = 2, the following relationship:

Cf ¼

0:6010 0:3162 Reb () 4



~ ¼ 1:1878 n e ¼ 5:2989 ; 3000 < Reb < 50000: C ð63Þ

The qualities of both fits, obtained with Eqs. (61) and (63) above can be judged from Fig. 11 where they are compared against the more rigorous fits corresponding to Eqs. (58) and (59). As we can

14

13

12

11 Eq.(58) Eq.(59) Eq.(61) Eq.(63)

10

9

8 2.4

2.6

2.8

3.0

3.2

3.4

3.6

3.8

0.5

log10 (Re b Cf )

 2n~  n~ e ðn~1Þ CReb 2 Reb  ¼ ðn~1Þ () C f Cf cRes c

e 2n~ Re C b

However, as discussed in great length in Pope (2000, Section 7.3.4) power law relationships such as the above are inappropriate and when used provide only approximate relationships valid for a limited range of Reynolds number. In fact Pope proposed a relationship for the power low coefficient that best fits data as (Pope, 2000, Eq. (7.158)) (for c = 2)

0.5

pffiffiffi pffiffiffiffiffi! Reb C f 2 1 3 pffiffiffi pffiffiffiffiffi ¼ ln þB þ X: 2j Cf j 2 2

erate different frictional relations (this is most prominent in the viscous laminar limit, where in channel flow we have Cf = 12/Reb as opposed to Cf = 16/Reb that holds for cylindrical Poiseuille flow). Similarly, the friction factor expressions valid under turbulent flow conditions are different corresponding to higher friction values for cylindrical tube flow than for channel flow, albeit in this case not by as much as in the laminar case, maybe to a higher value by about 5% instead of 33%; see Fig. 6. The corresponding to Eq. (60) parameters that are offered as fits of experimental data by Pope are (Pope, 2000, p. 312) (for c = 2) 8 ~ ¼ 8=7  1:1429 n < 0:3164 1=4  pffiffi 4=7 Reb ; 104 < Reb < 105 : Cf ¼ () e ¼ 8 42 :C 4  6:9925 0:3164

1/Cf

62

ð60Þ

similar to that used for turbulent channel flow—see Table 5 and discussion around it. However, it is well known that the geometrical difference between channel and Poiseuille flow is sufficient to gen-

pffiffiffiffiffi pffiffiffiffiffi  Fig. 11. Comparison of the Newtonian expressions for 1= C f vs log10 Reb C f , when c = 2. Black solid line: Eq. (58) obtained in this work through formal integration of the full mean velocity profile. Red line (with stars): Eq. (59) Colebrook correlation suggested by Bird et al. (2002). Blue line (with dots): Eq. (61), the approximate power-law Blasius relation for high Reynolds numbers. Magenta line (with squares) Eq. (63): the approximate power-law Blasius relation for low Reynolds numbers. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

63

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

see from there, the fit offered by Pope agrees very well with the general relationship offered by Bird et al., which is expected to be valid at the higher Reynolds numbers, whereas, the fit that we got comes closer to our relationship, valid for lower Reynolds num~ coefficients, correbers. It is of interest also to note, that the n sponding to the two power-law fits described by Eqs. (61) and (63), also shown in Fig. 11, are indeed variable and seem to be in agreement with the general rule, Eq. (62), offered by Pope. Indeed, ~ roughly corresponds to the optimum in the first case, the value of n value at the upper limit of validity Reb = 105, whereas in the latter case, the value of n⁄ found falls well in between the suggested optimum values corresponding to the limits of validity of that relation: ~ ð50000Þ < 1:1878 < 1:2375 ¼ n ~ ð3000Þ. Thus, Eq. (62) 1:1561 ¼ n seems to provide a fairly valid estimate for the (variable) power~ and justifies its use for the establishment of the law exponent n proposed formula for friction factor accounting for the effects of drag reduction—see Eq. (37) and others connected to it. A similar development can be followed for the viscoelastic MDR case. Substituting into Eq. (55) the approximation for the mean velocity profile provided by Eq. (54) gives: ! Z Res Z yc ub 2 þ þ þ þ þ þ   u u ¼ ðRe  y Þ dy þ ðRe  y Þ dy s s 1 3 us Re2s 0 yc 0" #y¼yc 1 2 @ y2 y3 y5 y6 y8 A  þ þ  ¼ 2 6Res 20ð14:5Þ3 30ð14:5Þ3 Res 80ð14:5Þ6 Res y¼0 0" #y¼yc 1 2 @ y3 y4 y6 y7 y9 A   2  þ þ 3 8Res 24ð14:5Þ3 35ð14:5Þ3 Res 90ð14:5Þ6 Res y¼0

y¼Res ! 2 1 þ yðlnðyÞ  1Þ þ By j Res y¼yc

y¼Res ! 1 1 2 1 y ðlnðyÞ  Þ þ By2  2 j 2 Res y¼yc   1 3 þ X;  lnðRes Þ þ B  j 2j

ð64Þ where X is so defined above as to represent the viscous sublayer correction. As in the Newtonian case, note that this correction is negligible at large Reynolds numbers being proportional on the first and higher powers of the inverse of the Reynolds number, as evaluated from the above equation. Eq. (64) with the aid of Eqs. (5)–(7) gives (taking into account that for pipe flow c = 2):

pffiffiffi pffiffiffiffiffi!   Reb C f 1 1 3 2 pffiffiffi pffiffiffiffiffi ¼ ln þB þ X: j 2 Cf j 2 2

ð65Þ

For the maximum drag reduction parameters, i.e. for j  1/11.7 and B  17, and for c = 2, Eq. (65) gives:

pffiffiffiffiffi! Reb C f 1 11:7 46:715 X pffiffiffi pffiffiffiffiffi ¼ pffiffiffi  ln  pffiffiffi þ pffiffiffi Cf 2 2 2 2 2  qffiffiffiffiffi 265:35 2129 pffiffiffiffiffi  ;  8:273  ln Reb C f  33:03 þ Reb C f Re2b C f

1 11:7 28:1 X pffiffiffiffiffi ¼ pffiffiffi  lnðRes Þ  pffiffiffi þ pffiffiffi Cf 2 2 2  8:273  lnðRes Þ  20:3 þ

6. Analysis of the potential contributions to Cf of MDR From the above relations for the skin friction coefficient in channel flow, and those derived in Appendix B for pipe flow, one can deduce parametric dependences for the skin friction coefficient as a function of the bulk Reynolds number. Here, we further explore this information in deducing at certain limiting cases the consequences of these relations to the viscoelastic and inertial contributions, IV and IR. In particular, Eqs. (26) and (B.6) can be solved in terms of the Reynolds stress contribution:

(1 IR ¼

3

16 C f  3Re ðb þ ð1  bÞIV Þ; pipe flow b

1 C 3 f

ð68Þ

2c  Re ðb þ ð1  bÞIV Þ; channel flow b

It is interesting to study Eq. (68) in the two limiting cases of a Newtonian turbulent flow (for which b = 1) and a viscoelastic turbulent flow under the condition of maximum drag reduction. In the first case, Eq. (68) gives: ðNewtÞ IR

¼

8 < 1 C ðNewtÞ  3

f

: 1 C ðNewtÞ 3 f

16 3Reb

; pipe flow

2c  Re ; channel flow

ð69Þ

:

b

ðNewtÞ

However, the relation between C f and Reb can be considered known; see Eq. (34) for channel flow, and Eq. (58) for pipe flow. This allows for a unique evaluation of the Reynolds stress contribuðNewtÞ tion for any given bulk Reynolds number. Results for IR as function of Reb0, in logarithmic scale, are shown in Fig. 12. In order to allow for a better comparison between channel and pipe flow, we have used twice the hydraulic radius in both cases, i.e. c = 2 for pipe flow and c = 4 for channel flow. For small Reynolds numbers, it is seen that the Reynolds stress contribution is quite different between the two cases. Furthermore, it increases substantially with the bulk Reynolds number, reaching a maximum at Reb  3500 for pipe flow, and at Reb  104 for channel flow. For large Reynolds numbers, beyond Reb  104, the results become almost identical, decreasing very slowly with increasing Reynolds number. Clearly, at high Reynolds numbers the effect of the geometry disappears with the hydraulic radius becoming a unified way of defining the bulk Reynolds number. In the second case, Eq. (68) gives: ðMDRÞ

where the parametric expression for the viscous sublayer correction in the second line has been found using a linear regression of the exact value X,ffiffiffiffiffibased on Eq. (64), in terms of a polynomial of pffiffiffiofp 1=Res  2 2= C f Reb for the given parameters, accurate within the range of interest to better than 0.5 parts per thousand relative error. A similar approach can also be followed in evaluating an approximate relationship for the channel flow. It leads to the following relationship:

ð67Þ

where the parametric expression for the viscous sublayer correction in the second line has been found using a linear regression of the exact value of X, based on an equivalent expression to Eq. (64), in terms of a polynomial of 1/Res for the given parameters, accurate within the range of interest to better than 0.1% relative error.

IR ð66Þ

47:0 64:2  ; Res Re2s

8   ðMDRÞ 16 > < 13 C ðMDRÞ  3Re b þ ð1  bÞIV ; pipe flow f b   ¼ : ðMDRÞ ðMDRÞ > 2c : 13 C f  Re b þ ð1  bÞIV ; channel flow

ð70Þ

b

ðMDRÞ

Even if the relationship between C f

and Reb can be also conðMDRÞ

sidered known, Eq. (70) shows that for any given Reb and b; IR

is

ðMDRÞ still a function of IV . One is reminded here the constraints ðlamÞ ðNewtÞ ðlamÞ ðMDRÞ 6 IR 6 I R and IV 6 IV . Therefore, the state of max0 ¼ IR

imum drag reduction may not be unique since for any value of the Reynolds stress contribution in the window between zero and the pure Newtonian value, there is a corresponding value for the viscoelastic turbulent contribution that results to maximum drag reduction.

64

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

approximately. Therefore, we believe that the zero Reynolds stress which is occasionally observed at the high or MDR regions cannot be sustained at high bulk Reynolds numbers.

7. Conclusions

IR

(Newt)

1E-3

Channel flow Pipe flow

1E-4

β0=0.9 3.0

3.5

4.0

4.5

5.0

5.5

6.0

log10 (Reb ) Fig. 12. The Reynolds stress contribution to the skin friction coefficient as a function of the bulk Reynolds number in logarithmic scale, Eq. (69). The skin friction coefficient has been calculated according to Eq. (58) for pipe flow (c = 2), and Eq. (34) for channel flow (c = 4).

It is worth mentioning here that some researchers report a zero ðMDRÞ Reynolds stress close or at maximum drag reduction, i.e. IR ¼0 (Warholic et al., 1999; Vlachogiannis et al., 2003; Ptasinski et al., ðMDRÞ 2003). In this case, we can solve Eq. (70) in terms of IV :



ðMDRÞ  IV 

ðMDRÞ

IR

¼0

8 1 ðMDRÞ Reb  b; 1 < 16 C f ¼ 1  b : 1 C ðMDRÞ Reb  b; f

6c

  ðMDRÞ

ðMDRÞ 

The results for IV

IR

¼0

pipe flow channel flow

:

ð71Þ

are shown as a function of Reb in ðMDRÞ

Fig. 13, forb = 0.9. The skin friction coefficient C f

is evaluated

using Eqs. (66) and (67) for pipe and channel flow, respectively. It is seen that under maximum drag reduction conditions, the Reynolds stress can remain zero only when extremely large values of ðMDRÞ

IV are developed. Actually, considering that there is shear thinning, b ? 1 as Reb ? 1 and the above estimates are considerably lower to the required values. This is a very unlikely situation except, perhaps, for Reynolds numbers in the window 104–105,

4

10

3

10

2

10

IV

(MDR)

Channel flow Pipe flow

1

10

β0=0.9 0

10

4

5

6

7

log10 (Reb0 ) Fig. 13. The viscoelastic stress contribution to the skin friction coefficient as a function of the bulk Reynolds number in logarithmic scale, under maximum drag reduction conditions, Eq. (71). The skin friction coefficient has been calculated according to Eq. (66) for pipe flow (c = 2), and Eq. (67) for channel flow (c = 4).

In this work, we have focused on the predictions for the skin friction coefficient and its contributions due to inertia, viscous and viscoelastic forces, for both laminar and turbulent wall bounded viscoelastic flows. For laminar flow conditions we have elucidated the role of the polymer shear thinning and showed that the use of the effective wall viscosity only partially can account for that effect. For the analysis under turbulent conditions we have used previous results obtained by direct numerical simulations of viscoelastic turbulent channel flow. A variety of well-known differential constitutive equations, such as the FENE-P, Giesekus and Oldroyd-B constitutive models, has been used. The results revealed that when viscoelasticity becomes important all contributions to the skin friction coefficient are significant. Among other conclusions, it is worth mentioning that, for constant Trouton ratio and friction Reynolds number, as the Weissenberg number increases the inertial contribution decreases and the viscoelastic one increases, eventually dominating the skin friction coefficient. For viscoelastic turbulent flow in a channel, we have also derived a universal fitting curve to the DNS results for the ratio of the achieved drag reduction to the asymptotic limit of drag reduction for large Weissenberg numbers, LDR, as a function of only one parameter, the friction Weissenberg number. By employing that curve and by using a previously established relation between the friction and the mean bulk Reynolds number for a Newtonian fluid, we have derived a new expression between the skin friction coefficient and the mean bulk Reynolds number. This involves two primary parameters, the zero shear-rate elasticity number, El0, which is a dimensionless effective relaxation time and controls the onset of drag reduction, and the limiting drag reduction, LDR, which represents the asymptotic degree of drag reduction at high Weissenberg numbers and controls the steepness of the skin friction curves. It also involves, as a secondary parameter, the relative wall viscosity lw that describes shear-thinning effects and it can be very well approximated by shear viscometric data. Using experimentally determined fits to the logarithmic region of the mean velocity and a formal integration we have established expressions for the inverse of the square root if the skin friction coefficient in terms of the friction Reynolds number for both the Newtonian and MDR limits of turbulent pipe flows. We have used the first, together with the universal drag reduction curve that we established from our DNS in order to generate predictions in Prandtl–von Karman coordinates for viscoelastic drag reducing turbulent cylindrical pipe flows. Our results reproduce both types A and B of drag reduction, corresponding to partially and fully extended polymer molecules, respectively, as first described by Virk (1975b). Furthermore, the predictions compare well with the experimental data by Shetty and Solomon (2009) provided a rescaling of the reported relaxation time in this work is performed. This is consistent with previous works suggesting that for associative polymers the effective relaxation time of a polymer solution may be different in extension than in a shear-dominated flow. Also, using the derived expressions for the skin friction coefficients, we have explored the maximum drag reduction state and we have shown that it may not be unique since for any value of the Reynolds stress between zero and the Newtonian value, there is always a corresponding value for the viscoelastic contribution which can result in maximum drag reduction. In the limit of zero Reynolds stress, as some researchers argue about being the case under maximum drag reduction, we have shown that the

65

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

viscoelastic contribution becomes extremely large at high bulk mean Reynolds numbers. Thus we believe that this is unlikely expect, perhaps, at small bulk Reynolds numbers. Acknowledgments K.D.H. would like to thank the Department of Mathematics, at the University of the Aegean, Greece, for the sabbatical leave of absence in the spring semester of the academic year 2011–2012, during which part of this work was conducted. Appendix A. Steady laminar viscoelastic flow in a pipe In pipe flow, the analysis for the skin friction factor proceeds along similar lines to that of channel described in subsection 2.1. Using a cylindrical coordinate system x⁄r⁄, where x⁄ is the axial coordinate and r⁄ is the radial coordinate, the linear momentum balance in the axial direction for steady, Poiseuille flow including the extra stress tensor due to the viscoelastic contribution of the polymer in the flow, is:

0 ¼ p;x þ gs

  1 d 1 d   dux þ   r srx : r    r dr r dr dr

ðA:1Þ

⁄

Multiplying by r , integrating with respect to the radial coordinate from the centerline, r⁄ = 0, up to r⁄ and dividing by r⁄ leads to the following integral momentum balance: 

1 du 0 ¼  r  p;x þ gs x þ srx : 2 dr

ðA:2Þ

Integrating Eq. (A.2) with respect to the radial coordinate, from r⁄ to the pipe wall R⁄, gives an expression for the velocity profile as a function of the radial distance:

ux ðr  Þ

p;x 2 1 ¼ ðR  r2 Þ þ  4gs gs

Z

ub ¼

g



ðlamÞ IV R2 s

g

16b0 16ð1  b0 Þ ðlamÞ þ IV0 ; Reb0 Reb0

s

r

 ^ rx dr :

ðA:3Þ

ðA:4Þ

;

ðlamÞ

and c = 2 has been used (as usual for pipe flow) in the definition of the bulk Reynolds number, Eq. (5). For a pure Newtonian fluid b0 = 1 and Eq. (A.8) reduces to the well-known expression Cf = 16/ Reb0 for smooth pipes. ðlamÞ It is also interesting to evaluate IV0 for different viscoelastic models. In the limiting case of the Oldroyd-B model, for which srx0 = dux/dr, it is easy to verify that IðlamÞ ¼ 1. Thus, Eq. (A.8) reV0 duces to Cf = 16/Reb0 for any b0 and Wes0, and Eq. (A.7) gives back the quadratic Poiseuille profile, ux = 2(1  r2), as it should. In general, in the laminar case, for a model that does not show a shearthinning behavior, the skin friction coefficient remains the same as in the Newtonian case. However, when the constitutive model allows for a shear-thinðlamÞ ning effect, IV0 is always less than one. For instance, if we consider the steady, laminar flow of a viscoelastic fluid which obeys to the FENE-P constitutive model (Housiadas and Beris, 2003), the following system of equations results:

8 9   C f Reb0 x 0 > r þ du b0 þ 1b ¼ 0> > > 4 dr f > > P > > < =  2 du 3 2 d x ¼0 ;  f  2 We f P > P > dr > > > > > > R1 : ; 2 0 ux ðrÞrdr ¼ 1



   dux   þ s ¼ p;x R =2:  rx  dr r  ¼R

By combining Eqs. (A.4) and (A.5), ub :

sw ¼

ðA:5Þ

sw can be expressed in terms of

4gs ub 4IðlamÞ þ V 3 : R R

ðA:6Þ

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2=3 d þ 1 þ 144 We d2 12 We 1 ðlamÞ

Iv 0 ðb0 ¼ 0Þ ¼

ux ðrÞ ¼ C f

Reb0 1  b0 ð1  r2 Þ þ 8b0 b0

Z

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 : d 12 We d þ 1 þ 144 We d2 8 We

ðA:10Þ

Eq. (A.9) can also be solved in the infinitely dilution limit, i.e. for b0 ? 1: ðlamÞ

Iv 0 ðb0 ¼ 1Þ ¼

4   sinh ðx=6Þ h d 2 þ 20 cosh x 14 þ 1728 We d 4 3 259270 We   2x 1 d 2 Þ: ðA:11Þ ; x  cosh ð1 þ 243 We þ11 cosh 3

Expressions (A.10) and (A.11) are visualized in Fig. 1; for more details see Section 2.1. For viscoelastic models for which the shear viscosity varies with shear rate, Eq. (A.8) may also be given in terms of the wall values as:

Cf ¼

16b 16ð1  bÞ ðlamÞ þ IV Reb Reb ðlamÞ

⁄  Using q u2 b ; ub and R as the characteristic scales for the pressure, velocity, and distance, respectively, the relation sw ¼ p;x R =2 gives Cf = p,x. Furthermore, when gp0 ub =R is used to scale the polymer extra shear stress, srx , where gp0 is the zero shear-rate viscosity of the polymer, Eq. (A.3) becomes:

ðA:9Þ

d and fP are the same as for the channel flow, and where We srx0 ¼ f1P dudrRx . Solving Eq. (A.9) we can find fP ; dudrx and then 1 ðlamÞ x 2 IV0 ¼  0 f1P du r dr. In the limiting case for which the fluid is pure dr viscoelastic, i.e. for b0 = 0, we can solve these equations analytically to obtain:

where the dimensional viscoelastic contribution is defined as RR  ðlamÞ IV   0 srx r2 dr . Furthermore, evaluating Eq. (A.2) at r⁄ = R⁄ produces the shear stress at the wall:

sw   gs

ðA:8Þ

where the dimensionless viscoelastic contribution is IV0     R1 s r2 dr; b0 is the viscosity ratio, b0  gs =g0 ¼ gs = gs þ gp0 , 0 rx0

R

In the Newtonian limit, Eq. (A.3) reduces to the well-known Poiseuille profile. Inserting Eq. (A.3) in (3) and performing the integration gives:

p;x R2  8 s

Cf ¼

ðA:12Þ

R1

where IV   0 srx r 2 dr. For a shear-thinning fluid, Cf based on the wall properties is shown in Fig. 2; see Section 2.1 for further discussion on the results. Appendix B. Stationary turbulent viscoelastic flow in a pipe

1

srx0 d^r:

ðA:7Þ

r

By integrating Eq. (A.7) with respect to r, from 0 to 1 and using R1 2 0 ux ðrÞrdr ¼ 1 (the total mass balance) the skin friction coefficient is obtained:

For fully developed turbulent flow in a pipe the Reynolds averaged governing equations are derived by first decomposing each dependent variable, f⁄, as follows:

f  ¼ f  ðr  Þ þ f 0 ðx ; r  ; u ; t  Þ;

ðB:1Þ

66

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67

where the overbar denotes time and spatial (in all the periodic directions) averaging and, by definition, f 0 ¼ 0. Because of the continuity equation, it is easily shown that ur ¼ 0, i.e. the average velocity component in the wall-normal direction is zero. The nontrivial component of the momentum balance is:

   1 d 1 d du 1 q   ðr ux ur Þ ¼ p;x þ gs   r x þ  r dr r dr r dr d   rx ;   r s dr

for accurate predictions of the skin friction coefficients and the resulted drag reduction are specific weighted averages of the Reyxy , respectively. nolds stress and the viscoelastic stress, ux ur and s References



ðB:2Þ

Integrating twice Eq. (B.2) with respect to r⁄, we get the velocity profile in the main flow direction:

 x ðr  Þ ¼ u

;x 2 p 1 ðR  r2 Þ þ  4gs gs

Z

R

r

srx d^r 

q gs

Z

R

r

ux ur d^r ;

ðB:3Þ

Eq. (B.3) can be written in dimensionless form by using the characteristic scales reported for the laminar flow:

Reb0 1  b0 ð1  r 2 Þ þ 4cb0 b0

ux ðrÞ ¼ ðp;x Þ 

Z

Z

1

r

srx0 d^r 

Reb0 cb0

1

ux ur d^r ; pipe flow;

ðB:4Þ

r

Eq. (B.4) is multiplied by r and is integrated with respect to r beR1 tween 0 and 1. Using 2 0 ux rdr ¼ 1; C f ¼ p;x and rearranging the result gives the final expression for the skin friction factor:

Cf ¼

16b0 16ð1  b0 Þ þ IV0 þ 3IR ; Reb0 Reb0

ðB:5Þ

R1 R rx0 r 2 dr and IR ¼  01 ðux ur Þr2 dr. Eq. (B.5) can also where IV0 ¼  0 s be reformulated using the wall viscosity as reference viscosity simply by replacing b0, Reb0, IV0 with b, Reb, IV:

Cf ¼

16b 16ð1  bÞ þ IV þ 3IR : Reb Reb

ðB:6Þ

This procedure is similar to that used by Fukagata et al. (2002) for Newtonian wall bounded flows and by Yu and Kawaguchi (2004). Eqs. (B.5) and (B.6) show that the total skin friction coefficient is decomposed in three individual contributions: a laminar (or viscous) contribution that is inversely proportional to the zero shear-rate bulk Reynolds number, a Reynolds turbulent contribution, 3IR, which is present for both Newtonian and viscoelastic flows, and a purely viscoelastic turbulent contribution, which is proportional to the quantity IV. Note that expressions Eqs. (B.5) and (B.6) have been derived by using as characteristic velocity the bulk velocity ub , which is a known quantity for experiments conducted with constant flux. When the experiments are conducted with constant pressure gradient, the suitable characteristic velocity is the friction velocity, us , Eq. (3), while the characteristic length is the same as before. In this case, the corresponding expressions for the skin friction coefficient in pipe flow is:

sffiffiffiffiffi   2 Res0 3 2ð1  b0 Þ ¼ 1  IR;s  IV0;s Cf 8 Res0 4b0

ðB:7Þ

where IV0,s and IR,s have been made dimensionless using us as the characteristic velocity. They are connected with IV0 and IR as:

Reb0 IV0;s ¼ IV0 ¼ cRes0 IR;s ¼

Re2b0

IR c2 Re2s0

¼

sffiffiffiffiffi 2 IV0 Cf

2 IR Cf

ðB:8aÞ ðB:8bÞ

The above expressions are useful because they analyze directly the effect of the Reynolds and the viscoelastic stress on the frictional drag. They also explicitly state that the most important quantities

Baik, S., Vlachogiannis, M., Hanratty, T.J., 2005. Use of particle image velocimetry to study heterogeneous drag reduction. Exp. Fluids 39, 637–650. Beris, A.N., Housiadas, K.D., 2010. Computational Viscoelastic Fluid Mechanics and Numerical Studies of Turbulent Flows of Dilute Polymer Solutions, contributed chapter in ‘‘Modeling and simulation in Polymers’’. In: Gujrati, Leonov (Eds.). Wiley-VCH. Bird, R.B., Armstrong, R.C., Hassager, O., 1987, . Dynamics of Polymeric Fluids, second ed., vol. 1. John Wiley and Sons, New York. Bird, R.B., Stewart, W.E., Lightfoot, E.N., 2002. Transport Phenomena. second ed.. John Wiley and Sons, New York. Clasen, C., Plog, J.P., Kulicke, W.-M., Owens, M., Macosko, C., Scriven, L.E., Verani, M., McKinley, G.H., 2006. How dilute are dilute solutions in extensional flows? J. Rheol. 50, 849–881. Colebrook, C.F., 1939. Turbulent flow in pipes, with particular reference to the transition region between smooth and rough pipe laws. J. Inst. Civil Eng. 11 (4), 133–156. Coles, D.E., Hirst, E.A., 1968. Computation of turbulent boundary layers–1968 AFORS-IFP Stanford conference. In: Proc. 1968 Conf., vol. 2, Stanford Univ., Stanford, California, USA. Dallas, V., Vassilicos, J.C., Hewitt, G.F., 2010. Strong polymer-turbulence interactions in viscoelastic turbulent channel flow. Phys. Rev. E 82, 066303-1-066303-19. DeAngelis, E., Casciola, C.M., Piva, R., 2002. DNS of wall turbulence: dilute polymers and self-sustaining mechanisms. Comput. Fluids 31, 495–507. De Gennes, P.G., 1986. Towards a scaling theory of drag reduction. Physica D 140, 9– 25. Dean, R., 1978. Reynolds number dependence of skin friction and other bulk flow quantities in-two dimensional rectangular duct flow. Trans. ASME I: J. Fluids Eng. 100, 215. Dimitropoulos, D.D., Dubief, Y., Shaqfeh, E.S.G., Moin, P., 2006. Direct numerical simulations of polymer-induced drag reduction in turbulent boundary layer flow of inhomogeneous polymer solutions. J. Fluid Mech. 566, 153–166. Dubief, Y., White, C.M., Terrapon, V.E., Shakfeh, E.S.G., Moin, P., Lele, S.K., 2004. On the coherent drag-reducing and turbulence-enhancing behaviour of polymers in wall flows. J. Fluid Mech. 514, 271–280. Escudier, M.P., Presti, F., Smith, S., 1999. Drag reduction in the turbulent pipe flow of polymers. J. Non-Newton. Fluid. Mech. 81, 197–213. Fox, R.W., McDonald, A.T., Pritchard, P.J., 2004. Introduction to Fluid Mechanics. John Wiley and Sons, New York. Fukagata, K., Iwamoto, K., Kasagi, N., 2002. Contribution of Reynolds stress distribution to the skin friction in wall-bounded flows. Phys. Fluids 14, L73. Gyr, A., Bewersdorff, H.-W., 1995. Drag Reduction of Turbulent Flows by Additives. Kluwer Academic Publishers.. Housiadas, K.D., Beris, A.N., 2003. Polymer-induced drag reduction: effects of the variations in elasticity and inertia in turbulent viscoelastic channel flow. Phys. Fluids 15, 2369–2384. Housiadas, K.D., Beris, A.N., 2004a. An efficient fully implicit spectral scheme for DNS of turbulent viscoelastic channel flow. J. Non-Newton. Fluid. Mech. 122, 243–262. Housiadas, K.D., Beris, A.N., 2004b. Characteristic scales and drag reduction evaluation in turbulent channel flow of non-constant viscosity viscoelastic fluids. Phys. Fluids 16, 1581–1586. Housiadas, K.D., Beris, A.N., 2005. Direct numerical simulations of viscoelastic turbulent channel flows at high drag reduction. Korea–Australia Rheol. J. 17, 131–140. Housiadas, K.D., Beris, A.N., 2006. Extensional behavior influence on viscoelastic turbulent channel flow. J. Non-Newton. Fluid Mech. 140, 41–56. Housiadas, K.D., Beris, A.N., Handler, R.A., 2005. Viscoelastic effects on higher order statistics and on coherent structures in turbulent channel flow. Phys. Fluids 17, 035106. Kalashnikov, V.N., 1998. Dynamical similarity and dimensionless relations for turbulent drag reduction by polymer additives. J. Non-Newton. Fluid Mech. 75, 209–230. Li, C.F., Sureshkumar, R., Khomami, B., 2006a. Influence of rheological parameters on polymer-induced turbulent drag reduction. J. Non-Newton. Fluid Mech. 140, 23–40. Li, C.F., Gupta, V.K., Sureshkumar, R., Khomami, B., 2006b. Turbulent channel flow of dilute polymeric solutions: drag reduction scaling and an eddy viscosity model. J. Non-Newton. Fluid Mech. 139, 177–189. Lin, C.S., Moulten, R.W., Putnam, G.L., 1953. Mass transfer between solid wall and fluid streams-mechanism and eddy distribution relationships in turbulent flow. Ind. Eng. Chem. 45, 636–640. Luchik, T.S., Tiederman, W.G., 1988. Turbulent structure in low-concentration dragreducing channel flows. J. Fluid Mech. 190, 241–263. Lumley, J.L., 1969. Drag reduction by additives. Ann. Rev. Fluid Mech. 1, 367–384. Min, T., Yoo, J.Y., Choi, H., Joseph, D.D., 2003. Drag reduction by polymer additives in a turbulent channel flow. J. Fluid Mech. 486, 213–238. Paterson, R.W., Abernathy, F.H., 1970. Turbulent flow drag reduction and degradation with dilute polymer solutions. J. Fluid. Mech. 43, 689–710.

K.D. Housiadas, A.N. Beris / International Journal of Heat and Fluid Flow 42 (2013) 49–67 Pope, S.B., 2000. Turbulent Flows. Cambridge University Press, Cambridge. Ptasinski, P.K., Boersma, B.J., Nieuwstadt, F.T.M., Hulsen, M.A., Van den Brule, B.H.A.A., Hunt, J.C.R., 2003. Turbulent channel flow near maximum drag reduction: simulations, experiments and mechanisms. J. Fluid Mech. 490, 251–291. Seyer, A., Metzner, A.B., 1969. Turbulence phenomena in drag-reducing systems. A.I.Ch.E.J. 15, 426–434. Shetty, A.M., Solomon, M.J., 2009. Aggregation in dilute solutions of high molar mass poly(ethylene) oxide and its effect on polymer turbulent drag reduction. Polymer 50, 261–270. Sureshkumar, R., Beris, A.N., Handler, R.A., 1997. Direct numerical simulation of the turbulent channel flow of a polymer solution. Phys. Fluids 9, 743–755. Tabor, M., De Gennes, P.G., 1986. A cascade theory of drag reduction. Europhys. Lett. 2, 519–522. Tsukahara, T., Ishigami, T., Yu, B., Kawaguchi, Y., 2011. DNS study on viscoelastic effect in drag-reduced turbulent channel flow. J. Turbul. 12, 1–25. Virk, P.S., Merrill, E.W., Mickley, H.S., Smith, K.A., Molo-Christensen, E.L., 1967. The Toms phenomenon: turbulent pipe flow of dilute polymer solutions. J. Fluid Mech. 30, 305–328. Virk, P.S., 1975a. Drag reduction fundamentals. A.I.Ch.E.J. 21, 625–656. Virk, P.S., 1975b. Drag reduction by collapsed and extended polyelectrolytes. Nature 253, 109–110. Virk, P.S., Sherman, D.C., Wagger, D.L., 1997. Additive equivalence during turbulent drag reduction. A.I.Ch.E.J. 43, 3257–3259. Virk, P.S., Wagger, D.L., 1990. Aspects of mechanisms in type B drag reduction. In: Gyr, A. (Ed.), Structure of Turbulence and Drag Reduction, IUTAM Symp. Zurich, 1989. Springer Verl., pp, 201–213.

67

Vlachogiannis, M., Liberatore, M.W., McHugh, A.J., Hanratty, T.J., 2003. Effectiveness of a drag reducing polymer: relation to molecular weight distribution and structuring. Phys. Fluids 15, 3786–3794. Warholic, M.D., Massah, H., Hanratty, T.J., 1999. Influence of drag-reducing polymers on turbulence: effects of Reynolds number, concentration and mixing. Exp. Fluids 27, 461–472. White, C.M., Mungal, M.G., 2008. Mechanisms and prediction of turbulent drag reduction with polymer additives. Annu. Rev. Fluid Mech. 40, 235–256. Xi, L., Graham, M.D., 2012. Dynamics on the laminar-turbulent boundary and the origin of the maximum drag reduction asymptote. Phys. Rev. Lett. 108, 028301. Yu, B., Kawaguchi, Y., 2003. Effect of Weissenberg number on the flow structure: DNS study of drag reducing flow with surfactant additives. Int. J. Heat Fluid. Trans. 24, 491–499. Yu, B., Kawaguchi, Y., 2006. Parametric study of surfactant-induced drag reduction by DNS. Int. J. Heat Fluid Flow 27, 887–894. Yu, B., Li, F.C., Kawaguchi, Y., 2004. Numerical and experimental investigation on turbulent structures in a drag reducing flow with surfactant additives. Int. J. Heat Fluid Flow 25, 961–974. Zaragola, M.V., Smits, A.J., 1997. Scaling of the mean velocity profile for turbulent pipe flow. Phys. Rev. Lett. 78, 239–242. Zaragola, M.V., Perry, A.E., Smits, A.J., 1997. Log laws or power laws: the scaling in the overlap region. Phys. Fluids 9, 2094–2100.