Predictions with low turbulent Reynolds number two-equation models of flow over a rotating cylinder in a quiescent fluid

Engineering Turbulence Modelling and Experiments 3 W. Rodi and G. Bergeles (Editors) © 1996 Elsevier Science B.V. All rights reserved.

101

Predictions with low turbulent Reynolds number two-equation models of flow over a rotating cylinder in a quiescent fluid A. Salhi a n d M. O n i r i D é p a r t e m e n t de P h y s i q u e , F a c u l t é des Sciences de T u n i s , 1060 T u n i s , T u n i s i a . Turbulent flow in a vicinity of a rotating cylinder in a quiescent fluid is examined via recently two-equation turbulence models. It is shown that the experimental d a t a of the mean velocity profile is reproduced by these models. Moreover, their predictions, which deviate from the usual two-dimensional logarithmic law, capture the logarithmic law proposed by Nakamura, Ueiki & Yamashita.

1.

INTRODUCTION

Curvature or rotation effects can strongly modify the structure of turbulence by stabilizing or destabilizing it. Since these effects occur in fields as diverse as engineering and geophysics, their predictions remain an important main objective of current turbulence research . As shown by Kristoffersen ci al [1] the low-Reynolds number second-order model of Launder &: Shima [2] (hereafter referred to as LS model) accurately captured rotation effects in the case of rotating channel flow. However, it yields results that are significantly different from those occurred in the case of a turbulent flow in a vicinity of a rotating cylinder in a quiescent fluid (see Andersson et al [3]). T h e purpose of this study, is to test the recently low-Reynolds two-equation models: the k — τ model of Speziale et al (4] (hereafter referred to as SAA model), and the k — c model of Rodi & Mansour [5] (hereafter referred to as RM model) in the case of this flow, which is a simplest kind of flow with curved streamlines and it satisfies a boundary layer approximation (see Townsend |6j). Nakamura ct al [7] have proposed a logorithmic law t h a t differs from that of the conventional two-dimensional turbulent layers. Even though, Andersson et al [3] have showed that the mean velocity profiles deviate from the conventional logarithmic law, they do not confirm the useful of the logarithmic law proposed by [7]. Moreover, Smith & Townsend [8|, who experimentally investigated a turbulent flow between two concentric cylinders with the inner one rotates, simply pointed out that near the rotating cylinder the ratio (r — R)/R, where R is the radius of the inner cylinder and r is the radial distance, is small and consequently the conventional logarithmic law and that one of [7] are not distinguishable. This point merits to be more clarified via the predictions of the SAA model and the RM model since they yield satisfactory results for the mean velocity profiles by comparison with the measurements of [3]. Also, we examine the radial distributions of the turbulence intensities given by these models in connection with the experimental d a t a of |3] and[7], and the predictions of the LS model given in [2j . T h e normal components of the Reynolds stress tensor are calculated using the non-linear k — € model of Speziale [9).

102

2.

GOVERNING EQUATIONS

2.1.

Curvature parameter

We will consider a full turbulent flow on a rotating cylinder in a quiescent fluid. A cylindrical coordinates (r, 0, z), where r is the radial distance from the axis (Oz) of the rotating cylinder and θ is the azinui thai angle, can be chosen to describe the flow motion. When the cylinder is very long the mean values depend only on the radial coordinate (see e.g. [6]). In this case, the mean velocity gradient VU is of the form / 0 (Vt/)

t i

0\

^

=

0

0

V ο

ο o/

where V is the peripheral component of the mean velocity. It is convenient to split the mean velocity gradient into its symétrie, 5, and antisymmetric, ω, parts (see e.g. Salhi et al [10]),

/ο

(£)« = s

1

1

0\

0

0

f . (a)« = s

U 0 oj S = rd(V/r)/0i\

\

0

1 +2Λ, 0\

-(1+2Λ.)

0

0

»

0

0/

where is the mean rate of deformation and R = (V/r)/S, is the ratio of the extra rate of strain V/r to the mean rate of deformation. The R parameter characterizes the curvature effects. If | R | < 1 the flow can be reasonably, approximated by the twodimensional flat-plate flow (see [7]). The particular value R = - 1 / 2 (i.e. the angular momentum is constant,V/r oc 1 / r ) characterizes the zero absolute vorticity case (i.e. ω = 0, see Cambon et al [11], Salhi h Cambon [12]). As pointed out by (7] the region where the angular momentum is constant, is not laminar but turbulent since the ratio of the RMS value in the peripheral direction to the mean velocity is constant. By using the displaced particle approach (see Tritton [13]) one can easily show that the Bradshaw number Β = 2R {\ + 2R ) is an indicator of the flow stability. The experiment of [7] reveals that V/r oc l / r where η > 0, so that Β < 0, the flow is then unstable or neutral. a

9

s

s

2

S

S

n

2.2.

Two-equation models to be adopted

In the background of two-equation models based on the eddy viscosity concept, we will adopt the nonlinear k — e model of Speziale [9] in conjunction with the SAA model and the RM model. In the case of the turbulent flow on rotating cylinder in a quiescent fluid, the Reynolds stress tensor components deduced from the nonlinear k — € model take the following forms, (1)

103

v

2

=

l

C ^Ç{y

k

+ 4

D

+ R,l

2

(2)

^=l ~l o^Y> k

(3)

C

uv = — / / 5 ,

(4)

t

where u,v and w are the fluctuating velocity components respectively, in the radial, peripheral and axial directions, and Co — 1.68 is a numerical constant (see [9]). Equations (1-4) will be solved in conjunction with the streamwise momentum equation,

(-ν**)*).

<*)

the transport equation for the turbulent kinetic energy (k — ΰ7ΰϊ/2, in cartesian tensor form),

D , Dt

1 d (, rdr\

v dk\ a Or J x x

n 2

k

and the model transport equation for the dissipation rate (e = form),

i/UijUij

in cartesian tensor

or the model transport equation for the time scale τ = k/e (in this case the term e in (6) is replaced by k/τ),

Η";έ(<"£>£)+<'-«·'·>ϊ·>*«Κ> k dr dr

τ dr dr

-(l-C / / ). £ 2

2

3

(8)

Here, D(.)/Dt = d(.)/dt-\- i/.V(.), ν is the kinematic viscosity, / i , / , / are damping functions that tend toward unity away from the wall. In the both models / i = 1 and a = σ = σ = 1.3. The SAA model (or the RM model) uses C = 1.44 (C = 1.45), C = 1.83 ( C = 2). The term F in (7) corresponds to the model proposed by [5] for the gradient production; involving in the exact transport equation for 6. This model will be adopted to a cylindrical coordinate system. As in some k — e models, the RM model uses the dissipation € = e - D\ where D is an extra term, which will be specified later, to define the eddy viscosity i/ 2

3

k

£ l

tl

e 2

£

Τ

e 2

3

£

t9

u = t

^4^.

(9)

In contrast to €, the quantity ? tends to zero at the wall. Here, / is a damping function that reduces the eddy viscosity near the wall (see Patel et al [14]); Ο is a coefficient that depends on the ratio of production to dissipation of turbulence energy ( P / c ) but is typically taken to be 0.09 as in the local equilibrium layer where Ρ = e (see Rodi [15]). However, in the k — τ model v is defined without this modification since the time scale r tends to zero at the wall, μ

μ

%

v = CpfpkT. %

(10)

We will now adopt to a cylindrical coordinate system the damping functions, the extra term D and the model for Pf proposed by [5],

104

2.3. T h e e x p r e s s i o n s o f t h e q u a n t i t i e s D, F , / > ht /μ nate system 3

(

a cylindrical coordi­

m

2

At the wall, the exact transport equation of the turbulent kinetic energy is reduced, 1 0 , Ok' _ —(//r — )

c(r = R)

(11)

Equation (11) can be used to define the extra term D since ? takes a zero value at the wall. However, due to the fact that (11) can lead to considerable numerical problem (see [14]), several authors have then proposed for D an alternative form. For the developed channel flow, Rodi & Mansour [5] have used the following equation, D = 2v

oy

where y is the cartesian coordinate perpendicular to the wall. We propose to transpose the latter equation as

.Or

D = 2//[4-v A?l , J /

since equations (11) and (12) yield the same result at the wall as can be shown from a near wall asymptotic analysis. Note that the curvature effects are not appreciable close to the wall. For shear-layer flows Rodi & Mansour have proposed for P ; which according to the DNS data, becomes negative near the wall; the following model, 3

c

k, >-\ c

(13)

where U is the mean velocity component in the direction Ox and ( . ) = d(.)/dy. For the particular values C? = 0.5 and C\ = 0.006 the RM model well fit the DNS data (see [5]). If the scheme used by Rodi & Mansour for the modeling of P , is adopted to a cylindrical coordinate system and the correlation ΰν^ is neglected for a first approximation, one obtains the following equation, y

3

P = 2 C*(S,)>

+ vC\{\

3

(

m

+

(4+i?,)G'f

+ R,){S,

+ -S)—k, r e

T

+

* £ ± M l ± l τ

s >

m

-^ (.C ).r--K'/ ).r.S' . 2

(14)

2

1

1

For simplicity, we will not use model (14), which, with respect to (13), has a complicated form. The following simplified form is used, Pf = 2iw Cl{S f t

+ uCl^S k .

tr

>r

(15)


Near the wall (where R < 1) the second term in the rigtli-hand side of (15) is negative {k r > 0 S < 0, S > 0) as the term involving C' in (13). The damping functions /2,/3 and can be adopted to a cylindrical coordinate system without difficulty. The SAA model uses the following expressions, 9

3

t

y

/ , = (! +

tr

3.45

)ia7i/i(r /70), y/TÏ7t

f

+

2

=

1 -exp(-r /4.9) \ / +

3

= l-^exp(-(/2

/6) ), 2

e t

where r = (u R/2i/)(l - (R /r )), R = kr/u i s th e turbulen t Reynold s numbe r an d u the frictio n velocity . While , th e R M mode l uses , +

2

2

T

et

/„ = 1 - exp(-0.0002r + - 0 . 0 0 0 6 5 r ) , f = ( l - 0 . 2 2 e x p [ - ( 7 2 +2

2

/6) ]^, 2

e <

T

is

105

It should be mentioned that the above expression for is due to [16].

3.

RESULTS A N D D I S C U S S I O N S

3.1.

Numerical scheme

The governing differential equations (5),(6) and (7), in which the advection term is neglected (D(.)/Dt % d(.)/dt), were solved by the control-volme method. This is performed via a fully implicit method where a second order, space centered differencing and a first order time diffe­ rencing is used (see [17]). Among others, Ekandcr [18] has used this numerical method, which is discussed in some detail in his report, to predict via the standard k — ( model and the algebraic stress model the rotation and curvature effects in the cases of plane rotating channel flow and circular channel flow. We have used the report [18] to carry out the numerical calculations. In addition, we have focused attention for integrating source terms S over the control volume and for the linearisation of the average cS" = S -f £,,Ψ by ensuring that S and S are respectively positive and negative (see [17]). The choice of an appropriate grid depends on the Reynolds number, R = V Rju, where V is the linear velocity of the rotating cylinder. 0 4 — 85 grid points across the flow domain were used in our calculations. The first point of the nonuniform grid chosen away from the wall at about r+ = 0.5 (or r = 2.5) when R = 3.9 10 (or R = 0.33 1 0 ) . For the choice of the nonuniform grid and the initial profiles of V, k and c- or τ we have used c

c

p

e

+

5

e

p

p

5

e

- the universal velocity distribution law of [7], which will be specified later, - the near wall asymptotic analysis, - the reduced form of the momentum equation in the local equilibrium layer (—r uv = 2

R u , P = e.) 2

2

T

The boundary conditions are, - At the wall: V = V ,k = 0,c = D,r = 0, p

- Far from the wall: V(r = r ) = 0,k(r = r ) = 0 , φ · = r ) = 0 , r ( r = r„) = r ( r = r„_i), n

n

n

where ΊΙ denotes the number of grid points. A number of tests have been undertaken, in particular, comparisons of numerical results with the exact solution for V in the case of laminar flow between rotating cylinders where the inner one rotates, have been made by setting the eddy viscosity equal to zero. Without introduction of relaxation and when the time-step size takes the value At = 1 0 ~ a fully converging solution requires a number of iterations greater than 20000. This iterations number gives rise to a maximum relative error for V, k,( or r about 10" 12, (computations were effected on the HP 9000 station). We found that the relation (v -f v )r S = Constant, can be also used as a useful way of testing the accuracy of the numerical solution in each grid point. 3

2

t

106

3.2.

Comparison with experimental

data

Figure 1 shows the mean velocity profiles (V/V versus (r — R)/R) yielded by the models: SAA, RM and LS for R = 10 . It appears that the predictions of the both two-equation models agree with the experimental data of [3|. Accordingly, the relevance of the logarithmic law proposed by [7] can be checked. p

r>

e

F i g u r e 1. M e a n v e l o c i t y profiles (V/V R M model; " " L S model;

·

v e r s u s r/R - 1) for R = 10 .

S A A model;

5

p

e

Experimental data of [3j.

Among others, | 3 | , [7| and [10] have applied their mean flow data to the usual logarithmic law,

™

+

Ξ

^

= ό > ' -

+

5

<

(

1

6

)

where r+ = (R - r)u /v. Their comparisons clearly show that the data do not fit this later law. This fact has motivated |7| and |20] to introduce the dimensionless variable η such as (άη/dr) = ( U 7 2 ) / ( I > T ) , and derive, after integration of the streamwise momentum equation, the following law, T

3

3

t

VN

+

ξξ

V - -V — = 1.824 In r + 9.8. r

(17)

+

u

T

As shown by figures 2, giving the variation of VA versus r+ for R — 0.33 10 and R — 10 , the both two-equation models yield results that exhibit the destabilized curvature effects as the experimental data of |3], which deviate from the law (16). Pointing out that [3] have used the correlation proposed by Theodorsen Sz Regier [21] for the skin friction coefficient, Cf = 2ul/V , +

5

e

5

e

2

p

= 4.071og (fl /cV) - 0.6, 10

e>

(18)

107 to obtain the wall friction velocity since their data are considerably scattered and not conclusive (see fig. 3). For R > 0.55 10 the both two-equation models yield results for Cj that are close to those calculated from (18). This can explain the fact that the models better reproduce the experimental mean velocity profile when R = 10 (fig. 2-b). 5

%

5

e

ιο·

ίο

ίο

1

ίο

2

ιο·

3

F i g u r e 2. M e a n v e l o c i t y profiles ((V - V)/u R = 10 , (b). S A A model; R M model; of [31 p

T

5

t

10

10

1

v e r s u s r+ for R LS m o d e l ; ·

10*

2

= 0.33 ΙΟ , (a) a n d Experimental data 5

e

Figure 4 gives the variation of VN versus r for R = 3.910 . As can be seen there is no expected deviation of numerical results from the law (17). Moreover, the differences observed between computations and (17), become négligeable if the wall friction velocity deduced from the correlation proposed by |7], +

+

5

e

-4f

= 3.091og (/t /cV) + 4.79. l0

(19)

fV

is used instead of the predicted value for u . T

F i g u r e 3 . D r a g coefficient Cj v e r s u s 2R

e

108

Since the k—€ model cannot reflected the sensitivity of the turbulent shear stress to curvature an empirical ad-hoc modification of Ο is proposed. On the assumption of local equilibrium of turbulence Ρ = € the right-hand sides of (l)-(4) are respectively set equal to the algebraic expressions of u , v , w and uv formulated by using the algebraic stress model of Rodi [15]. This scheme implies € = 1/(1 4- fi Β), where fi is a constant taken to be 1.83, recall that Β is the Bradshaw number. It appears that this had-hoc modification does not modify in appreciable manner the radial variations of the mean velocity (see fig. 4) and the turbulence intensities. μ

2

2

2

μ

F i g u r e 4. M e a n v e l o c i t y profiles (\V - (R/r)V]/u v e r s u s r+ w i t h ( c i r c l e s ) a n d wi­ t h o u t (line) m o d i f i c a t i o n of Ο for R = 3.9 10 . ( a ) : S A A m o d e l ; ( b ) : R M m o d e l . p

T

5

μ

e

F i g u r e 5. T u r b u l e n c e i n t e n s i t y profiles for R = 10 . S A A model; R M mo­ del; • • • • L S m o d e l ; · E x p e r i m e n t a l d a t a of [3]. ( a ) : radial c o m p o n e n t (VvP/Vp); ( b ) : p e r i p h e r a l c o m p o n e n t (\Z7/V ); ( c ) : axial c o m p o n e n t (y/u7/V ). 5

e

p

p

109

By comparing the predicted turbulence intensities to the experimental data of [3] for R = 0.3310 and R = 10 . It appears that the two-equation models predicted better the radial distributions of the turbulence intensities when R = 10 . This can be tied to the fact that for the former value, the flow relatively exhibits a complete turbulent character as shown by the intermittency factor variation (see fig. 7 in [3]). The radial variation of the turbulence intensities for R = 10 given in figure 5 indicates clearly that the two-equation models are more close to the experimental data than the predictions of the LS model. According to [3], the dificiency of the LS model is due to the modeling of the pressure-strain correlation. On the other hand, the both models qualitatively reproduce the experimental data of [7] when r/R > 1.075. However, for r/R. < 1.075 (or ln[\ - (R/r) ] < - 2 ) there is a significant deviation as shown by figure 6 giving the radial evolution of the dimensionless turbulent kinetic energy k+ = (R/r) (2k/u ) for R = 3.9 10 . e

5

5

e

5

e

5

e

2

2

2

5

T

(r/R

F i g u r e 5.

e

- 1 )

Hl-RVr ) 2

F i g u r e 6. R a d i a l v a r i a t i o n o f k

+

=

(R/r) (2k/u ) v e r s u s 1 Ι Ι ( 1 - Λ / Γ ) for R = 3.9 10 . 2

2

2

2

T

5

e

•

4.

S A A model; R M model; E x p e r i m e n t a l d a t a o f [7].

CONCLUDING REMARKS

The turbulent flow in a vicinity of a rotating cylinder in a quiescent fluid was examined using two low-Reynolds number models: the SAA model and the RM model. More attention is focused on the adoption to a cylindrical coordinate system of the different terms in the governing equations. Comparisons with avalable experimental data have been made. In contrast to the LS model, the RM and SAA models yield satisfactory results for the mean velocity profile. Accordingly, the relevance of the universal logarithmic law proposed by [7] is checked. It appears that the both two-equation models qualitatively reproduce this law. The ad-hoc modification for the eddy viscosity coefficient involving the Bradshaw number is introduced in order to take into account the curvature effects. It is shown that this correction does not lead to a significant difference with respect to the case where Ο is constant. μ

110

Even though there are no great difference between the predictions of the two-equation models, the RM model remains more performent than the SAA model. This is due to the values of the constant C ; since far from the wall the damping functions tend toward unity and is négligeable. t 2

REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23.

R. Kristofferesen, P.J. Nilsen and II.I. Andersson, In. Eng. Mod. Exp. Elsevier (1990). B.E. Launder and N. Shima, AIAA J.27 (1989) 1319. H.I. Andersson, B. Johansson, L. Lofdahl and P.J. Nilsen 8th Symp. T.S.F. (1991) 30-1. C.G. Speziale, R. Abid, E.G. Anderson E.G., NASA CR-182068 (1990). W. Rodi, and N.N. Mansour, J. Fluid Mech. 250 (1993) 509. A.A. Townsend, The structure of turbulent shear flow, Cambridge Univ. Press (1976). I. Nakamura, Y. Ueiki and S. Yamashita, 4th Symp. on T.S.F, (1983) 2.21. G.P. Smith, and A.A. Townsend, J. Fluid Mech. 123 (1982) 187. C.G Speziale, J. Fluid Mech. 187 (1987) 450. A. Salhi, T. Lili and J.F. Sini, Phys. Fluids 5 (1993) 2014. G. Cambon, J.P. Benoit, L. Shao and L. Jacquin, J. Fluid Mech. 287 (1994) 175. A. Salhi and C. Cambon, Revisiting rotating shear flow at high shear rate, submitted to J. Fluid Mech. (1005). D.J. Tritton, J. Fluid Mech. 241 (1992) 503. V.C. Patel, W. Rodi and G. Scheuerer, AIAA J. 23 (1985) 1308. W. Rodi, Ph.D. Thesis, London (1972). T.H. Shih and N.N. Mansour, In. Eng. Mod. Exp. Elsevier (1990). S.V. Patankar, Numerical Heat Transfer and Fluid Flow, M.G. Hill (1980). II. Ekander, TRITA-MEK-88-06, The Royal Institute of Tech. Sweden (1988). N. Kasagi and S. Ilirata, ASME Pub. 75WA/IIT (1975). Y. Furuya, Y. Nakamura and S. Yamashita, Mem. Fac. Eng. Nagoya Univ., (1978) 30-1. T. Theodorsen, A. Regicr, NACA report No. 793 (1944). B.E. Launder, G.J. R.eece and W. Rodi, J. Fluid Mech., 68 (1975) 537. C.G. Speziale, N. Mac Giolla Mhuiris, J. Fluid Mech., 209 (1989) 591.