A tentative approach to the direct simulation of drag reduction by polymers

l~-~ian

Fluid Medmks

ELSEVIER

J. Non-Newtonian Fluid Mech., 60 (1995) 277 301

A tentative approach to the direct simulation of drag reduction by polymers Paolo Orlandi * Universiti?t di Roma, "La Sapienza', Diparttmento di Meccanica e Aeronautica, Via Eudossiana 16, 00184 Rome, lta(v

Received 28 April 1995: in revised form 20 July 1995

Abstract An efficient technique for drag reduction uses dilute solutions of a few p.p.m, of polymers. A possible reduction in drag of up to 80% is achieved. Several experimental observations have been made which tend to indicate that the polymers modify the turbulence structures within the buffer layer. Flow visualisations have shown that the changes consist of a weakening of the strength of the streamwise vortices. Existing literature reveals no attempts of numerical simulation of this phenomenon. In this paper an approach is pursued by using a constitutive equation which relates the elongation viscosity to the local properties of the flow. According to this model this viscosity is large in zones where the amount of strain rate is greater than the amount of vorticity, and is zero when the vorticity exceeds the strain rate. Simulations have been performed in a "minimal channel" to give good resolution with a limited number of grid points. The accuracy of the method is tested by comparison with the results of other techniques. For simulations with polymers, quantitative comparisons cannot be made, but the results reproduce the qualitative outputs of the experiments. The mean streamwise velocity is modified in the buffer layer; the peak of the streamwise turbulent intensity, in wall units, increases and its maximum moves far from the wall; and the vertical turbulent intensity is largely reduced in the wall region. An interesting outcome from both the simulation and the experiments is that the strength of the longitudinal vortices is reduced when the polymers are present. Ke)words: Channel flow; Dilute polymers; Direct simulation

* Tel: 39-6-44585878, Fax: 39-6-484854. e-mail [email protected] 0377-0257/95/$09.50 © 1995 ElsevierScience B.V. All rights reserved SSD1 0377-0257(95)01388-1

278

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272-301

I. Introduction

Viscous drag reduction is one of the most important techniques for reducing fuel consumption in transportation systems. Evidence of drag reduction, both by passive and active control, has been observed experimentally, but a clear explanation of the mechanism has not been given, mainly because the experimental observations cannot describe all the details of the flow in the wall region. In the literature there is a large number of papers devoted to the experimental study of the wall layer structures in drag reducing flows [1-3]. A common feature of these papers is the idea that drag reduction is due to modifications of the wall layer structures, particularly in the buffer region, the most active region in wall bounded flows. From flow visualisation, Tiederman et al. [3] concluded that the polymer solution inhibited the formation of low speed streaks and that, when these form, their spacing, in wall coordinates, increases with polymer concentration. In the last decade the direct simulation of wall-bounded flows has been very useful for obtaining a deeper knowledge of turbulence structures and of their role in wall friction. The importance of the bursting events, experimentally observed in Ref. [4], has been confirmed by direct simulation in Ref. [5]. The large amount of data generated by the direct simulation leads to a better understanding of the mechanisms that govern turbulence production, redistribution and dissipation. Attemps of simulate drag control have been made only recently since they require computational times much larger than those required for the plane channel. A first attempt to simulate the drag reduction by riblets can be found in Ref. [6]. Later, a more detailed simulation [7] was performed showing a drag reduction in accord with the experimental observation. In this case the numerical simulation requires a sophisticated numerical model to describe the geomertry of the riblets. In the case of drag reduction by a dilute polymer solution, numerical studies have not been performed to date, perhaps because of a lack of reliable constitutive equations relating the stresses to the strains. With the present paper a systematic study of the effects of polymers is initiated starting with a very simple constitutive equation, based on experimental observations. There are three important observations on which our model is based. The first is that the polymers do not interact with each other, with the consequence that the viscosity is proportional to the concentration. The second is that the polymers produce an effect on the bulk properties of the flow when they are in an extended state, while they are in a random coiled state, the bulk properties are those of the solvent alone. The third and most important observation is that the polymers must be in regions of the flow where the rate of strain prevails over the amount of vorticity to reach and maintain the extended configuration. With these observations, a very simple model for the elongational viscosity has been made by introducing a constant which is related to the concentration of the polymers. To be a realistic and invariant model, the relation between the stresses and the strains should be evaluated in the principal axis. This requires the evaluation of the eigenvalues and eigenvector of the tensor of rate of strain at each grid point. To overcome this large increase of computational time, the elongational viscosity has been introduced in the normal stresses referred to the

P. Orlandi /J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

279

coordinate axes, as a result the model is not frame invariant. However, this is a satisfactory approximation because, from the data of the direct simulation of the channel without polymers, it has been observed that the principal axes of the tensor of rate of strain in the buffer layer are approximately aligned with the coordinate axes. We consider the effect of the polymers in a "minimal channel" (see Ref. [8]); this is done to reduce the computational costs since the simulations of the full channel require large memory occupancy and hundreds of hours of CPU time on supercomputers. The simulation of the minimal channel unit [8] revealed that very good results are obtained at least as far as the wall structures are concerned. In Ref. [8] a pseudospectral code was employed, while in Ref. [6] a finite-difference method was used to perform similar simulations. Finite differences produced trubulent stress distribution in the wall region in good agreement with the results of the full channel simulation [5] and in the central region with the pseudospectral simulation of the minimal channel [8]. In the present paper the choice of the minimal channel has been done because it captures the complexity of the wall structures, but the lateral dimensions of the channel have been increased with respect to those in Ref. [8] to describe the growth of the spacing of the streaks when the polymers are present.

2. Physical model The most difficult task in the study of drag reduction in dilute solutions of polymers is constructing a satisfactory model to represent the effects of the polymers on the constitutive equations. The purpose of the present study is to make a tentative approach in this direction with a very simple anisotropic model for the viscosity. The model is based on experimental and theoretical observations of polymers in a dilute solution. A few p.p.m, of polymers in water changes the characteristics of the solution and at very low concentration, the first assumption is that the polymers do not interact with each other, giving an effect that is proportional to the concentration of the polymers. The theoretical analysis of this relation has been given in Ref. [9] for a suspension of force-free particles, where the correct way to do the ensemble averages is performed and the bulk stress owing to the presence of particles is derived. This theory of dilute suspension gives an average deviatoric stress owing to the particles which is proportional to the rate of strain of the flow. The quantity relating the stress to the rate of strain can be considered as a viscosity which depends on the kind of particles. The model in Refl [9], for a suspension of very slender rigid rods, shows that very high values of elongational viscosity, without any substantial increase in shearing viscosity, can be achieved. This situation is very similar to that encountered with polymers. When the polymers are in the random coiled state, it has been observed [10] that there is no change in the properties of the solution, whereas, when the polymers are in an extended configuration, there is a change in the properties of the solution that produces substantial drag reduction.

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

280

In order to reach their extended state, the polymers should be in a region of the flow where the strain is strong enough to maintain them oriented with the direction of the principal axes of the rate of strain. When the polymers are in a flow region where the vorticity is very intense, they will stay in the random coiled state and they will not affect the bulk properties of the solution. To construct a simple model which represents these observations at each grid position, the eigenvalues and the eigenvectors of the rate of strain tensor S o the mean square strain tensor S~jSji, and the mean square rotation rate RijRji must be calculated. With these quantities, a viscosity/~p, due to the effect of the polymers, is assumed to be tip = C 1

pp = O,

RijRji ~ s(jSj~i j,

Ro.Rj~ ~ SiiSji ,

Ro.Rj~ >~ so.SJ ~,

SgjSjs <~ ST.

siJsi i ~ ST,

(1) (2)

The quantities of these equations and in the remaining part of the paper are made dimensionless with respect to the centreline velocity of the Poiseuille flow UCL and the half channel dimension h. The Reynolds number thus is Re = UcLh/V. C is a constant that in the model is thought to be related to the polymers' concentration, and ST is related to the terminal relaxation time of the polymers. The polymers are elongate and give a contribution to the bulk properties of the solution only when the rate of strain is larger than ST. A study of the sensitivity of the model to changes in the constants C and ST has been performed. A substantial effect has been observed on the variation of C. It has been found that for a high value of C ( C = 20) the flow becomes laminar. However, the effect of ST is negligible when S T ~ 1. It is very difficult to relate the values of the constants to real polymer solutions, but I wish to emphasise that the aim of the present study is to perform qualitative comparisons with the experimental observation and if the results are realistic, to extract an explanation on the physical mechanism producing the drag reduction. In a dilute solution, each polymer contributes to the bulk stress of the solution. The relation between stress and rate of strain is then

:c~'~ = ~psk.

(3)

This equation is written in the coordinate system oriented with the principal axes of the rate of deformation tensor. This hypothesis is related to the physical observation that the polymers produce effects on the elongational viscosity and not on the shear viscosity. The stresses in the strain coordinates are related to the stresses in the physical coordinates through

vik V~q E~'~j = 2~ '~ 6kq,

(4)

where the V~ are the eigenvectors of the tensor of rate of strain. This calculation must be performed at each grid point. This is not a difficult task, but it requires a large amount of CPU time. To reduce the computational time the following relationship between the polymers stresses and the rate of strain is used:

:c~ = ~pS~j,

i =j,

(5)

P. Orlandi /J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

Y~' = 0,

281

i :~j.

(6)

This is not tensorial invariant but approximates Eq. (3) satisfactorily. In fact, the calculation of the eigenvectors of the tensor rate of strain in the turbulent channel at R~= 180 (R~= U~h/v), in all points where the condition [R~Rj,I~S,jS,, and S~jSj~ >~ Sr holds, shows that the probability distribution of the cosine of the angles of the eigenvectors of the largest positive eigenvalue (Fig. la), has a peak when it .20

.15h, i,~

J~

.10-

.p

a)

.::,'

.05-

.00

i

.0

.2

.4

.6

.8

1.0

cos(angle)

.20

.15~

.1o

i b)

:.7

i

.057j- .,-.

..

,,

.,. • :...

. o o

....................

.0

oo.,

: ..........

.2

-

.4

'

.6

.8

"

1.0

cos(angle) Fig. 1. P r o b a b i l i t y density function ( P D F ) o f the cosine o f the angle in the buffer layer a) o f eigenvectors of the p r i n c i p a l strain rate direction (positive eigenvalue) . . . . V~, V~, l'~, b) between vorticity vector a n d p r i n c i p a l strain rate directions positive, .. • negative, intermediate eigenvalues.

282

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272--301

is oriented at 45 ° with respect to the streamwise direction. The largest positive eigenvalue has been chosen because only stretching events characterised by two negative eigenvalues have been considered. These events are the most efficient to keep the polymers elongated. Also the orientation of the vorticity vector has been checked, and, as expected, Fig. l b shows that it is aligned with the intermediate eigenvalue. Eqs. (5) and (6) introduced in the momentum equations given 8~ + uj 8X j

8x, +-Re kSx~

8x~ 2#'"'

ax,j

(7)

This system of equations, together with the continuity equation, describes the evolution of a turbulent flow inside a plane channel when initial and boundary conditions are given. Periodicity in the spanwise and streamwise directions and no-slip walls in the normal direction have been assumed. The size of the computational domain in the spanwise direction (L,-~ = 0.574~) captures about two streaks in the simulation of pure water and is adequate to reproduce the experimental observation that with the polymers the spacing between the steaks increases. In the review of Robinson [11], who analysed channel [5] and boundary layers data [12], it is reported that the streamwise structures have diameters of the order of d + = 20-40 and length of l + = 100-200. The + superscript indicates, as usual, quantities scaled in wall units, e.g. W + = L , U~/v, with U ~ = ( r w / p ) ,/2. If the channel width is W + = 200 and the length is L + = 1000, these structures can be described with a limited number of degrees of freedom. With the grid 32 x 97 x 64 the streamwise and the spanwise resolution is close to that in the full channel simulation [5], and so the vortical structures in the wall region are accurately represented. As shown in previous simulations of the "minimal channel" [6,8], the reduced dimensions of the channel width cause unreal profiles of the Reynolds stresses in the central region of the channel. In fact, the structures in the central part have dimensions comparable to the channel height, and the reduced lateral dimensions of the computational box do not allow the interaction of several of these large structures. A non-uniform spacing in the vertical direction was used, in consideration of the fact that the vortices in the wall region are smaller than those in the centre of the channel. The first direct simulations of turbulent channels were performed by pseudospectral methods with Chebishev polynomials in the vertical directions and Fourier expansions in the periodic directions. Recently, it has been shown that finite differences are as accurate as pseudospectral methods, if global conservation properties are preserved [6,13]. Finite difference schemes have the additional advantage of a very simple treatment of viscosity variations. The numerical scheme, described in Ref. [6], is based on an implicit scheme for the viscous terms, and an explicit scheme for the nonlinear terms. The viscosity describing the effect of the polymers (Eq. (1)) is evaluated at the previous time step. The fractional step developed in Ref. [14] is used for the time advancement. FFTs in the x3 and x~ directions are employed for the "pressure" equation yielding a tridiagonal matrix with the reduced wave numbers in the main diagonal. This non-iterative procedure

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1095)272-301

283

produces a solenoidal field within roundoff errors. The calculations were done both on a RISC IBM workstation and on the C R A Y YMP84 at N A S A Ames. The CPU time required on the R I S C 380 with a grid 32 × 97 x 64 is about 25 s for a full time step, which is composed of three substeps as described in Ref. [13]. Jim6nez and Moin [8] did a systematic analysis at each Reynolds number, on the minimal dimensions of the channel required to sustain turbulence. As in the full channel [5] the Reynolds number is Re = Uh/v =4200. At this Re and with the streamwise length 2~h, it was observed that in order to sustain turbulence with a fixed mean flow rate it is necessary to introduce a random perturbation with a large amplitude (~ = 0.25U) at t = 0. This perturbation is reduced in the wall region to avoid very small time steps during the initial transient.

3. Results and discussion The initial transient in a direct simulation is not physically interesting since it depends on the initial perturbation. However, the wall friction time history shows that it increases, reaches a maximum and then settles down oscillating around a value 0.002. By a non-dimensional time [ = t U c u h = 500, the turbulence has reached a statistical steady state. Fig. 2a shows the near-wall profiles of turbulence intensities in wall units compared with those in the full channel [5]. To check the dependence on grid resolution in the streamwise directions two grids were used, 32 x 97 × 64 and 64 x 97 x 64. The good agreement between these cases shows that the resolution and the numerical scheme are adequate to predict the rms turbulence profiles and thus to describe the dynamics of the structures in the near-wall region. In my opinion, only when the interactions among the vortical structures are represented will the correlations' profile agree with those measured in experiments. For a further check of the quality of the simulations Fig. 2b shows that the velocity profiles, in wall units, agree with those in the full channel. A comparison with the experimental results [15] at Re~ = 142 and Re~ = 208 has not been done since this comparison was performed in Ref. [5], with an excellent agreement between real and numerical experiments. The better results with the finer grid are due to a lower value of the wall stress. In fact, the profiles normalised by the centreline velocity do not change nor do the contour plots of ~ , , , the vorticity component producing high and low speed streaks. This comparison is a p r o o f that the 32 × 97 z 64 grid is adequate and all the results discussed in the remaining part of the paper have been obtained with this resolution. It could be asserted that the prediction of second-order correlation does not prove the quality of a direct simulation and that a better check involves the higher order statistics such as skewness S ( ~ ) = (u'/3 )/(u'i 2)3'2 and flatness F(Ui ) = (u'~4 )/ (u'~2) 2 so that whether the statistics are far from Gaussian (where S(U,)= 0 and F ( U ~ ) = 3 . 6 ) could be judged. Regions of high intermittency, characteristic of turbulent flows are those where high values of flatness occur. The analysis of skewness and flatness was also presented in Ref. [13] to prove the quality of the finite difference scheme. In that paper the authors claim that very good results

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

284

X× X×

2

x

a) 1

0

0

50 '

1O0

150

200

20

15

+_ 10b) 5-

.

10-I

.

.

.

.

.

.

.

10°

.

.

.

.

.

.

.

.

.

.

101

.

.

.

.

.

.

i

102

. . . . . . . .

103

Fig. 2. a) r.m.s., b) s t r e a m w i s e velocity profile, in wall c o o r d i n a t e s for the c h a n n e l flow at R~ = 180. x , Ref. [5]; - - , 64 x 97 x 64; . . . . , 32 × 97 x 64.

could be obtained only by high-order upwind schemes, moreover they have shown that, by centered schemes, the numerical dissipation produced results in poor agreement with those by full spectral simulations. In this paper the simulations were performed by centered schemes and the results are as good as those in [13]. In Fig. 3 a - c the skewness of the three velocity components are compared with those of the numerical simulation [5] and with those in the experiment of Ref. [16]. The present results, particularly for the streamwise components, agree with the experimental results. Both simulations predict the location where the skewness changes sign. It is

285

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995)272 301

"2_

".5/

-10

•

o

2o

t

o

O

~

O

x~

•

~

•

°

8o

•

100

1.0-

.5^

G

b)

.0 x -.5'

-1.C

2'o

40

6o

80

~00

2'0

4b

6o

~o

,oo

1.0

.5

~

-7,

.0

-.5 ¸

-1.£

b

Fig. 3. Skewness profiles of a) streamwise, b) normal, c) spanwise velocity components: x. Ref. [16]; O, Ref. [5]: - - - , present results.

286

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272-301

difficult to explain why the present profiles of S ( U 1 ) (Fig. 3a) agrees better with the measured one than that of the full simulation, in regions far from the walls, where the width of the channel is not large enough to reproduce the motion of the largest scales without the unphysical effects of periodicity condition. Fig. 3b shows that the position of S(U2)= 0 in the present simulation is closer to the wall than in the experiment. The reason, as pointed out in Ref. [13], is due to the lack of resolution in the vertical direction. In fact Rai and Moin [13] performed two simulations, the first with 101 non-equally spaced points in the vertical direction and the second with 129 as in the spectral simulation. With the second grid, the S ( U 2 ) was in perfect agreement with that of the spectral simulation. Fig. 3c shows that the spanwise component S(U3) has a smaller vlaue than the other two components in agreement with the other results. The agreement with the full channel simulation can be observed also for the flatness (Fig. 4a c). The numerical simulations show in Fig. 4b that the flatness F(U2) of the vertical component increases near the wall, while the experiments show that it decreases approaching the wall. Kim et al. [5] claim that the the measurements are suspect near the wall. The present results show smaller values of F(U2) than those for the full channel and this is due to the lack of resolution in x2. Satisfactory results for the spanwise flatness component F(U3) are also found near the wall. I do not wish to claim that the present results are as accurate as those by a pseudospectral simulation. My intent was to show that, contrary to the usual belief, finite difference simulations of the "minimal channel", with a centered scheme for the nonlinear terms, produce satisfactory results for studying the turbulence physics of the wall region in channel flows. Thus this method is adequate for performing simulations with non-isotropic viscosity, varying with the local flow field, which can be very cumbersome with a pseudospectral method. However pseudospectral methods are used in large eddy simulations where an eddy viscosity is introduced [17], but they require smaller time steps for the larger Courant-Friedrichs-Lewy restrictions than in finite difference. The flow field at t = 520 has been taken as initial condition for the cases with polymers. Calculations have been performed without any additive (C = 0) and with two different values for the constant C in Eq. (I). Since the present model is based on the viscosity in Eq. (1), it is worthwhile to analyse how that quantity is distributed in the turbulent channel. Fig. 5a and b show, in four equally spaced vertical planes, contour plots of 1-(RisRji/So.Sji ) (solid lines), with the restrictions in Eq. (1), superimposed respectively to o9~ and ~ox2, where positive vorticity has dashed contours and negative vorticity has dotted contours. The result is that, in the buffer region and mainly where o~x, is concentrated, the viscosity of the present model is active. The highest values of the viscosity of the polymers are located in regions of high gradient of ~ox2. At C = 20 the simulation predicts that the turbulent flow becomes laminar with a parabolic velocity profile. For C = 10 the simulation produces a drag reduction of 27.7% and for C = 5 of 23.4%, these values have been obtained by the time history of the mean pressure gradient, required to keep a constant flow rate given in Fig. 6. In fact the average pressure gradients, performed

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

287

15

10a)

5-

0

0

iO

40

~+

60

8'0

100

b) 5

x

0

•

•

20

4'0

60

80

100

20

40

60

8'0

100

15 ¸

10

a.

b

Fig. 4. Flatness profiles of a) streamwise, b) normal, c) spanwise velocity components: ×. Ref. [16]; Q, Ref. [5]; - - - , present results.

P. Orlandi I J. Non-Newtonian Fluid Mech. 60 (1995) 272-301

288

t, .,',' ,,

',

"

,;r~:~..-.,

~

.-"~I--,'~

.. ~--.:- .... .:~:I

'~. "-

"~

• -'..

'--"

Q

~"--~ ..T,' ". "X-",'

'-"

', '~ :

,,

t-~:" ',',I' ..... " >o~kl ",i~ .-" .... ,'//' ,'" ):....,i c~.,!.

~,..~-'..::,.

'

£,.L.~I

'~

, ,'. ~

~-_o~',,,,,,

~' ,.~.~- .,, :::;~.

",,~,::--,,::~

; "--.~.'~1:~

...':,?-,,i!1

,

,

.-" ~ ¢ . 1 ' ,

v,,.,

: ~^ %;'1 ,, ,.. ~,<,,111

,~ ,. -,-- A . , M

',.

~,~.

"-v-"

",'-.:',~,.,-..

~1

.':'A

~:,.

,..,..--~'..' "...<~,,,',..=.~I

~.:~._. ~.

....

,,--

,~. ,,-~ ...'

. . . .

;>

~- ~.

._.: ~.~,,

0 I I I I

, >~-~,, ~--.,

..-

i ,;~= --~,, b.T.'.:

I~,.',

,, ",,

~ii

,

,

,,

/

..

"L

, .... :..~,~ .o

,

I;~=-- .'. ,' " : - , "-'--" ," I~,~:":~.~'v" ' : . . . . .

F~,:;--, ~,,~ i "-" *"

,. ,2.~..,,~

'-... ~.._:*~.'~

,~.~ " . ~ 0 - -'.

, " ,'-

-.~.

t

@

-I

-,

. .::.9.

I!.',, ~ , ' ], -':-, " ~" ~"~'" " .

", ,

K'~: *''-; . . . . . . ~ " , ,:-'-"A

"':', % "',. "'-'q

~,,.~;-.

..... .. ,:~..~...

~F~.'~..~ "' ."::--'Z-";.-'.

c~

.E_

"",~

¢) ,,,'h

',::i::~.-~,..... ,, , :~ -,, ,..-

..-:-~ -i!!i)l

, ,.. ,.-

. ',~'~,~]

"-' :',--

".L~

i ~ -~-,-,,

~-~.'. ~.~ ~'

."

.

C~

~,~"

/ ,~

'

"',

t'k~.~*-.'"- " .~ .----::-_~

r.,

~-.--

/. _: "-,":~

~

I', '5': ',',,

~

~,,~

....

..~/- !.~.~

t

~-.,, '.,

,t

.

~ '~k'~,l

'i"

g

:.%--_~*~1

,--

/.--:,', ,,.%.-7-"£1

I

,2 : .~ ',,'" -.-:.) . "- , ~'.,,., --., .~.~ ~'~.'t.-,...,~., ............ . ,,p~.,,-

0

k

~)

',~~

r-~.,'," ;,'" "--'X/; i '-;

i~..~--~,;

"~,. ;

..-,,,

~,,¢.,e,..._...

~",~,-,

";'-

"-';:';',':~1

.....~.;., ¢-,

I

,,,~ ~ .?. , q

,,'!:

,-'.."

'..~ . ~ 0

~,:::~..::::.~ ,..-;--?~

,H

m

d~

~

I~ '~ ~..~.,-:"

".,

,'. i~ ,,ilt

'.. ".,~'~,~

'q~-'-:i'-.~_']

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995)272 301

2~9

for t > 580, for C equal to 0, 5 and 10 respectively are --0.002102, -0.001609 and -0.001518. The case without additives (C = 0) has been performed to compare the flow structures with and without polymers. It is difficult to relate the values of the constant C to the amount of polymer concentration and to the type of polymers used in the experiments [10,18,19] but, as mentioned before, in this study the intent is not to reproduce a real experiment, first because the experiments were not at this low Reynolds number, and second because Eq. (1) is not an exact model of polymers dissolved in water. The aim was to analyse the tendency of the flow to reduce the drag when the concentration increases. Lumley [20] defined drag reduction as the reduction of skin friction in turbulent flow below that of the solvent alone. With this definition, only when the flow remains turbulent do polymers produce drag reduction. Indication of whether the flow is still turbulent is given by the profiles of the turbulent statistics and of the streamwise velocity in wall units. The mean pressure gradient time history (Fig. 6) shows that for all values of C, a different statistical steady state is reached. The sharp decrease at t = 520, the time when the polymers are thought to initiate their effect, is caused by the viscosity increase. The most interesting observation is that the time intermittency is reduced with the increase of polymers' concentration: this effect is related to the large modifications of the wall structures as shown later. In a review article Tiederman [21], analysing experimental results in pipe and in large aspect ratio channel flows, asserted that the interpretation of the data has improved but a full understanding of dilute polymer flows has not yet been reached. He summarised the main aspect of the flow, which I wish to recall in order to see which of them the present simplified model is able to reproduce. Regarding the mean velocity profile, Tiederman [21] claims that, for both Newtonian and dilute -.0010

-.0015

.... - ..........

ii

"'''"'''"'

:iii

•

''4

,""

/ 1

/l

V

-.0020

-.0025

500

, 600

Fig. 6. M e a n p r e s s u r e g r a d i e n t t i m e h i s t o r y f o r dilute p o l y m e r s ( C = 10).

' 700 solvent,

800

dilute p o l y m e r s ( C = 5), "

290

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

25

20-

""

•

ii: ,,*" ,-''"X... x • ,oo" .....'" ~,o X X.X ..--

..q" X

...o,( .............""

15 +~

I0 . ......-''"

5

,......,. y'"'"''" r

10 -]

'

I

[

'

'

'

'1

l0 °

.

.

.

.

.

'

I l l

10 ]

.

.

.

.

.

'

;

11

,

,

1

l0 z

Fig. 7. Streamwise velocity profiles in wall coordinates ×, solvent, • polymers Ref. [19]; - .... polymers (present results C = 5), • - - (l/k) In x2~ + 5.0.

,

103 solvent,

polymers' profiles, the UT = x~ relationship holds, and that the viscous layer extension does not change. The buffer region is extended in drag-reducing flows yielding a thicker wall layer. The velocity profiles in wall units, measured in Ref. [19] in pure water and in a dilute solution of 2.1 p.p.m, of polymers (Separan AP-273), compared with the present results are shown in Fig. 7. The trend of the measurements to have the same thickness of the viscous layer and a thicker buffer layer is reproduced. The quantitative differences are due to two reasons: the first is that the Reynolds number in the experiment is three times larger than in the simulation, the second is that the simulation was done in a channel of reduced width which could contribute to the discrepancies at large x~. As mentioned before, the aim of the present simulation is mainly devoted to understanding and representing the wall region. The laser velocimeter measurements of turbulent intensities described in Ref. [19], scaled with inner variables, showed that the peak value of u']+ occurs at larger values of x + and is broader in flows with polymers than with respect of the Newtonian case. Fig. 8a shows that the trend of the measurements is predicted by the simulation; however, also in the Newtonian case, the simulation produces a faster decrease of u']+ , in the log region, than that in the experiment. The reason should be related to the different Reynolds numbers; in fact, as shown in Fig. 3, the present results agreed with those of the full simulation [5] that were in a perfect agreement with those in the experiment of Ref. [15]. Regarding the different profiles of u'~ in polymers' solution, Tiederman [21] claims that "this increase must be due in part to the higher mean velocities in the outer region of these drag-reducing flows". In my opinion, this is not a clear explanation. However, from the simulation contour plots of u~, ,2 in a X l - X 3 plane, (e.g. at x~- = 15) for the solvent and for

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

291

the polymers, are easily obtained. Fig. 8b and c show that in the polymers' solution the levels are reduced with respect to solvent and that the peaks are located in a reduced number of zones. These numerical visualisations prove that

2-

t

~

--..

•

•

1-

a/

20

40

60

80

100

b)

Xl

Fig. 8. a) Streamwise velocity r.m.s, profiles in wall coordinates • solvent, • polymers Ref. [19]: solvent, polymers (present results C = 5); contour plots of u'~2 in x~ x3 planes at x~ ~ 15 b} solvent, c) polymers C = 5.

292

P. Orlandi /J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

1.5~

].0 ........

.5

.0

a)

0

,

i

r

i

20

40

60

80

100

b)

o

z1

Fig. 9. Normal velocity r.m.s, profiles in wall coordinates • solvent, • polymers Ref. [19]; - - solvent, - - - polymers (present results C = 5); contour plots of u~2 in x~-x3 planes at x + ~ 15 b) solvent, c) polymers C = 5.

the increase o f u'l + , is due to the decrease o f the wall stress. T h e s i m u l a t i o n s by the same x l - x 3 c o n t o u r plots for u22 (Fig. 9b a n d c) show that, in the p o l y m e r s ' solution, there is a very large r e d u c t i o n o f the vertical velocity fluctuations. This r e d u c t i o n is n o t c o m p e n s a t e d b y the r e d u c t i o n in wall stress with the p o l y m e r s a n d thus the u~+ profile with the p o l y m e r s has values smaller t h a n those in p u r e water, c o n f i r m i n g the m e a s u r e m e n t s [19] in Fig. 9a. F r o m these results it is clear t h a t the r e d u c e d activity in the wall region is m a i n l y due to the decrease o f the intensity a n d the n u m b e r o f the streamwise vortices as is discussed in Section 4.

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

293

The behaviour of the Reynolds stress from the experimental data is more controversial. In fact, in Refs. [22,23] it was found that the viscous stress and the Reynolds stress do not add to yield the expected linear variation, while in Ref. [19] a linear relationship was found. Willmarth et al. [22] conjectured that the reason was due to the inefficient mixing of their solution. Bewersdoff and Berman [23], however, suggested that the deficit could be due to the fact that the viscosity of the polymers contributes as a shear viscosity. The results of the present simulation, shown in Fig. 10, predict the same trend of the measurements in Ref. [19], confirming that the modelling of an elongational viscosity and not of a shear viscosity is satisfactory to reproduce the polymers drag reduction. Quadrant analysis of the u] and u" (Fig. 11), normalised by the corresponding r.m.s. fluctuations, shows that events in quadrant 2 and in quadrant 4 are those contributing to the Reynolds stress. From the experimental data, Tiederman [21] observed a rotation of the principal axes that were more aligned with the laboratory axes. This rotation produced lower u2 and a lower level for (u'lu2). The results of the simulation at three different distances from the wall are shown in Fig. 11. As t\~r the solvent the higher intensity of the negative Reynolds stress in the wall region is due to a reduced number of events in quadrants 1 and 3 compared to quadrants 2 and 4. The reduction of the Reynolds stress in the polymers solution is caused by the strong reduction of u2 fluctuations. From Fig. 11, this reduction is not directly observed, but considering that the fluctuations are normalised by the r.m.s, value at the corresponding distance from the wall, this statement is understood. The decoupling of the u'l and u" fluctuations is caused by a reduced amplitude of ~,), ,.

1.0 '... %'"'"....,.... •

•

'~t ....

"'""'.., •

-{:

.... "'"....

.5-

V

.0

I -.5

Fig. 10. a) Reynolds stress profiles in wall coordinates • solvent, • polymers Ref. [19]: polymers (present results C = 5), " total stress.

solvent,

294

P, Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272 301 5,

I

"L',." ": "i • • .. :~.....o~.. • • • ~, :'.7~'.~ "'.. ". ~'~."~ ~i." " .

0"

:,.

• ".i'..,..'. .................... ~.~~.,~..~:~:.-.-

~.',~ ...

................

~)

•

•

"':

";":.:~'~.C'_

"° t

:.'

. . . . .~, .: .,

..... ..'...

..~

.. -~.. ,.

•

•

..

•

..

•

-51 -5

-5

•

,

! ~.,.

•

... • .~

•

:,..~..~ ~:,'. ~;.;.,:~.. "~.". . . . . . . . . . . . . . . . . . . . . .

-...

•

~;~;~s~

°l

"'~

b)

~.'4:

• .;.2'~

.

,£.~.. :.; .. ! . ].:

-5 I

-5"

-5

5

5

~

•. •

..~ ~.'" "~:3'.~

•

•

"

. ~:.~. ~ . . . ~

:..

: . ::,2.~~::. E

.

..:"-~gi;2) ?t,s.,: ..

I

.................... ~'2

,¢,1

0

. •vj.°

~...

o. .................. =!~.~

.' "

•

'

.'

' ...

-:'?.

t ", .-.j. ".

. :.~:." ~-:~'~.'.-G',,:.

•

.

......?.g.?:.,'~..:~. -,.. ....

•

,::....................

o

c)

.

' ~;~'i;:.

• ~...

.....

J

.

i

1 .5

0 t t Ul/Ulrms

5

-5

6

ui/ui..,

Fig. 11. Quadrant analysis of u'~ and u~ for solvents (left) and for polymers (right) at a) x~ =5, b) x + = 15, c) x~ = 50. The axes are normalized by the corresponding r.m.s, fluctuations.

T h e wall s t r u c t u r e s c a n be visualised b y c o n t o u r plots in vertical cross sections o f s t r e a m w i s e v o r t i c i t y a n d velocity• F r o m the l a t t e r q u a n t i t y , ejection m e c h a n i s m a n d the f o r m a t i o n o f i n t e n s e shears is o b s e r v e d , while f r o m the c o n t o u r plots o f ~ox,, the role o f the l o n g i t u d i n a l vortices o n the s t r e a k f o r m a t i o n is d e d u c e d . I n Figs. 12 a n d 13 c o n t o u r plots in f o u r e q u a l l y s p a c e d vertical sections o f s t r e a m w i s e voriticity a n d velocity are s h o w n for the s o l v e n t a n d for the p o l y m e r s ' s o l u t i o n . W i t h the

295

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272-301

6/

,

,,,_a

j ".y_, &

E © [,

f~' .,..;

~ . S i '~

k tN

--&

x

~

!?:'l

¸¸

©

E

7, --& ©

e-i

?.x

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272-301

296

o -I

0

c~

C~

© e~

Dv

e~

I

r~

E

©

o o L~

u.,

~G0 ~

~0

P. Orlandi J. Non-Newtonian Fhdd Mech. 60 (1995) 272 301

297

solvent alone, there is a strong correlation between the low and high speed streaks and ~ , . This correlation has been predicted also by a quasi 2D model [24], and it was demonstrated that the high speed streaks contribute for the major part to the wall stress. While a large portion of the wall channel is occupied by the high speed streaks, on the other hand the high turbulence level is located in the low speed streaks and mainly between the longitudinal vortices. Fig. 12a shows that when a longitudinal vortex is near the wall, it induces vorticity of opposite sign at the wall: this could be the mechanism creating new vortices that, convected far from the wall are responsible for the turbulence generation. However it should be remembered that vortex stretching and vortex tilting are responsible for the increase of the vorticity field. Fig. 13a indicates that the central part of the channel is less active. In this region the large structures are more convoluted, as observed in a large number of flow visualisations. The results of this simulation agree with the previous observations in the full channel [5] and in the minimal channel [8]. The contours of streamwise vorticity, in the presence of polymers, (Fig. 12b) show that the effects of the elongational viscosity, mainly in the buffer region, reduce the number of streamwise w)rtices and their intensity with respect to the Newtonian flows. These weaker vortices create low and high speed streaks of reduced intensity, as shown by the smoother distribution of streamwise velocity causing a reduction both in mean wall stress and in turbulent fluctuations. From this instantaneous picture it follows that the spacing between the streaks increases, as experimentally observed [2]. The streaks are satisfactorily detected by contour

a)

~

'

-

~

/

_

,

~:~.~,,.¢ .....

.~!!!!~;¢::Z~ '~ ~

~

~

;gl

b)

...................

Fig. 14.. Contour plots of normal vorticity in x, .v~ planes at x~ = 15 for a) solvent, b) polymers positive. """ negative values.

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

298

plots in xj - x3 planes of co~2 as shown in Fig. 14. The patterns in the two cases are very different, for the solvent alone Fig. 14a shows that the streaks are closer than in the flow with polymers (Fig. 14b). In the first case the streaks are more intermittent than in the former case; this means that the elongational viscosity reduces the number of streamwise vortices leading to a reduction of opposite sign vorticity generated at the wall. This process is related to the regeneration process of streamwise vortices, since in the dilute solution of polymers the turbulent activity is thus reduced. The reduction of intermittency in the dilute solution is confirmed by the reduction of a factor of two of the flatness in the wall region. Robinson [11] analysing the

1.0

.5-

..

.0

.

.--"~

..................

.............S.?.?.......~....................................... a)

.......~i::S ........................................ .....

-.5

0

.~-°°

2'5

5'0

7'5

100

1.0

.5-

.........

.0i

IIIIIIU

-.5 0

2'5

~t

5'0

,...|..

b)

7'5

Fig. 15. C o r r e l a t i o n velocity profiles at x + = 10 for a) solvent, b) polymers, vertical, ' " streamwise components.

100 spanwise,

P. Orlandi / J. Non-Newtonian Fluid Mech. 60 (1995) 272 301

299

boundary layers data [12], observed that the high speed streaks are shorter than the low speed streaks, the same behaviour occurs with the polymers. From contour plots similar to those in Fig. 14, not reported here, it has been observed that the regions where co,. is concentrated are shorter than those where there is ~,),,, and this holds for solvent and polymers solutions. The quantitative spacing of the streaks is obtained by the spanwise correlations of the three velocities components. Fig. 15 shows the correlations for the solvent and the polymers at x~ = 10, the most representative distance to evaluate the streaks spacing. The correlations show that with polymers the correlation R~t reaches a value of 140, which is in a good agreement with the relation given in Ref. [2]: 2 ÷ = 1.9DR + 99.7 with DR = 23.4 the percent drag reduction in this simulation. Tiederman [21] in his discussion section claims that the mechanism which leads to the drag reduction was not determined in a conclusive way, but that it was important to note that the dilute polymer flows were not laminar. Moreover he affirmed that the turbulence structure was modified in various ways but that all the features of the Newtonian wall structure remained. With the present results 1 hope to have shed some light on the drag reduction mechanism. In fact, for the Newtonian flows [24] it was observed that the key mechanism of the wall stress is the presence of streamwise vortices. Then for an optimum control leading e.g. to drag reduction it is necessary to reduce the number and the strength of these vortices. The dilution of polymers was modelled by an elongational viscosity, a function of the local properties of the flow. This viscosity reduces the number and strength of these vortices leading to the drag reduction and leaving the flow turbulent, as stated in Ref. [21].

4. Conclusions

The present simulation in a minimal channel predicts the near wall layer structures when the solvent alone is considered: the results for the second and higher order statistics agree with those measured and with those calculated in a full channel. The minimal channel unit has been used to perform a tentative approach to study drag reduction when polymers are added to the solvent. The polymers' solution requires a constitutive equation for the stress strain-rate relation. In this preliminary study an heuristic model that does not fulfill the requirements of rational continuum mechanics has been used. The model is based on the idea that the polymers are active when elongated, while they do not act when they are in the random coiled state. The state of the polymers depends on the local enstrophy and on the square of the strain-rate tensor S~S;,. When the enstrophy prevails the polymers are in the random coiled state and they do not produce any elongational viscosity, while the elongational viscosity is high when S~/Si, reaches a certain threshold level and when it prevails over the local enstrophy. A simple linear relationship between normal stress and strain allows prediction of the main features, experimentally observed, in dilute polymers' solutions. That is, the mean velocities and turbulent quantities, scaled in the wall variables, do not change in the

300

P. Orlundi / J. Non-Newtoniun Fluid Mech. 60 (1995) 272 301

viscous layer while they change in the buffer layer. The streamwise turbulent stress in wall units does not increase appreciably and moves its maximum towards the central part of the channel. The more interesting outcome of the present direct simulation is that, as from flow visualisations [2], the longitudinal vortical structures decrease their strength and increase their spacing when the polymer concentration is increased. Quantitative comparisons were not done because at present it is difficult to quantify the effect of the polymer concentration on the viscosity reduction, and because the computers do not allow the simulation at the high R e numbers of the experiments. The first improvement on the model is to introduce more sophisticated models for the polymers as e.g. the Oldroyd type B model, which Lumley [1] proved to be equivalent to the dumbbell model. In this case the direct simulation becomes very onerous since it is necessary to solve a system of partial differential equation for the six components of the stresses. Because the streamwise vortices play a key role for drag reduction, and that these structures are very elongated in the streamwise direction, the Oldroyd type B model could be introduced in a 2D simulation, representative of a x2 x3 cut of a full simulation. This quasi 2D approximation has been attempted [25], showing that the resolution of equations of the polymers stresses are difficult. A further direction could be to eliminate from the upper convected Maxwell model, the material time derivative while keeping the nonlinearities. This approach was used in Ref. [26] and it seems to be promising.

Acknowledgements This work started with an idea of J. Lumley while he was visiting the engineering faculty in Roma, and I am greatly indebted to him. I am especially indebted to the fruitful discussions with J. Jimdnez and P. Moin about the minimal channel and with K. Shariff on the construction of the polymers model. Finally I would like to thank L. Piomeli who revised the first draft of the paper. I would like to thank one of the referees who sent Ref. [26] to me, which will be useful for my future work. The research was supported by Ministero Pubblica Istruzione and the Center for Turbulence Research. The computers of N A S A Ames Research Center were used.

References [1] J.L. Lumley, Drag reduction in two phase and polymer flows, Phys. Fluids, 20 (1977) $64 $71. [2] D.K. Oldaker and W.G. Tiederman, Spatial structures of the viscous sublayer in drag-reducing channel flows, Phys. Fluids, 20 (1977) S133-S144. [3] W.G. Tiederman, T.S. Luchik and D.G. Bogard, Wall-layer structure and drag reduction, J. Fluid. Mech., 156 (1985) 419-437. [4] H.T. Kim, S.J. Kline and W.C. Reynolds, The production of turbulence near a smooth wall in turbulent boundary layers, J. Fluid. Mech., 50 (1971) 133 160. [5] J. Kim, P. Moin and R. Moser, Turbulence statistics in fully developed channel flow at low Reynolds number, J. Fluid. Mech., 177 (1987) 133 166.

P. Orlandi / J. Non-New:onian Fluid Mech. 60 (1995)272 301

3111

[6] P. Orlandi, A numerical method for direct simulation of turbulence in complex geometries, Annual Research Briefs 1989, Center for Turbulence Research, 1990, pp. 215 229. [7] H. Choi, P. Moin and J. Kim, Direct numerical simulation of turbulent flow over riblets, J. Fluid. Mech., 255 (1993) 503 539. [8] J. Jim~nez and P. Moin, The minimal flow unit in near wall turbulence. J. Fluid. Mech., 225 (1991 213 240. [9] G.K. Batchelor, The stress system in u suspension of force-free particles, J. Fluid. Mech.. 41 (I 97[)t 545 570. [10] P.S. Virk and D.L. Wagger, Aspects of mechanisms in type B drag reduction, in A. Gyr ted.), Proc. Second IUTAM Symp. on Structure of Turbulence and Drag Reduction, Springer-Verlag. Berlin, 1990, pp, 201 213, [I 1] S.K. Robinson, Coherent motions in the turbulent boundary layer, Ann. Rev. Fluid Mech.. (1991) 610 639. [12] P.R. Spalart, Direct numerical simulation of a turbulent boundary layer up to R,, = 1410. J. Fluid Mech. 187 (1988) 61 98. [13] M.M. Rai and P. Moin, Direct simulations of turbulent flow using finite-difference schemes, J. Comp. Phys.. 96 (1991) 15 53. [14] J. Kim and P. Moin, Application of a tractional step method to incompressible Navier Stokes equations, J. Comp. Phys., 59 (1985) 308 323. [15] H. Eckelmann, The structure of the viscous sublayer and the adjacent wall region in a turbulent channel flow, J. Fluid. Mech., 65 (1974) 439-.459. [16] R.S. Barlow and J.P. Johnston, Structures of turbulent boundary layers on a concave surface. Rep. M D-47. Department of Mechanical Engineering, Stanford University, Stanford CA, 1985. [17] U. Piomelli, P. Moin and J.H. Ferziger, Model consistency in a large eddy simulation of turbulent channel flows, Phys. Fluids, 31 (1988) 1884 1891. [18] B. Gumpert and C.K. Yong, The influence of polymer additives on the coherent structures of turbuelent channel flow, in A. Gyr (Ed.), Proc. Second IUTAM Symp. on Structure of Turbulence and [)rag Reduction, Springer-Verlag, Berlin, 1990, pp. 223 232. [19] T.S. Luchik and W.G. Tiederman, Turbulent structure in low-concentration drag-reducing channel flows, J. Fluid. Mech., 190 (1988) 241 263. [20] J.L. Lumley, Drag reduction by additives, Ann. Rev. Fluid Mech., 1 (1969) 367 383. [21] W.G. Yiederman, The effect of dilute polymer solutions on viscous drag and turbulence structure, in A. Gyr (Ed.), Proc. Second IUTAM Symp. on Structure of Turbulence and Drag Reduction, Springer-Verlag, Berlin, 1990, pp. 187 200. [22] W.W. Willmarth, T. Wei and C.O. Lee. Laser anemometer measurements of Reynolds stress in a turbulent channel flow with drag reducing polymer additives, Phys. Fluids, 30 (1987) 933 935. [23] H.W. Bewersdorff and N.S. Berman, The influence of flow-induced non-Newtonian fluid properties on turbulent drag reduction, Rheoh Acta, 27 (1988) 130 136. [24] P. Orlandi and J. Jim~nez, On the generation of turbulent wall friction, Phys. Fluids. 6 (1994) 634 641. [25] P. Orlandi, G.M. Homsy and J. Azaiez, A viscoelastic model for direct simulation of polymers drag reduction, C.T.R. Proc. Summer Program 1992, Stanford University, Stanford, CA, pp. 165 174. [26] M. Renardy, On the mechanism of drag reduction, J. of Non-Newtonian Fluid Mech., 59 (1995:) 93 101.