The prediction of an axisymetric turbulent jet by a three-equation model of turbulence

FOR FLUID DYNAMICS

TECHNICAL NOTE 114

THE PREDICTION OF AN AXISYMMETRIC TURBULENT JET BY A THREE-EQUATION MODEL OF TURBULENCE

S. BIRINGEN

JUL! 1975

TECBNISCHE HOGESCHOOL DEtH VLiEGTUIGSOUWKUNDE

BIBLIOTHEEK Kluyverweg 1 - DElfT

7 OKI 1975

~A~

-~O~-

RHODE SAINT GENESE BELGIUM

~VW

TECHNICAL NOTE 114

THE PREDICTION OF AN AXISYMMETRIC TURBULENT JET BY A THREE-EQUATION MODEL OF TURBULENCE

S. BIRINGEN

LIST OF SYMBOLS

1. INTRODUCTION

1.1 Literature survey

1.2 Statement of the problem

2. THE MATHEMATICAL PROBLEM 2.1 Governing equations 2.2 Closure assumptions

2.2.1 The turbulent dissipation

l 1 1

4

5

5

8 8

2.2.2 The diffusion terms 9

2.2.3 The pressure-velocity gradient correlations 11 2.2.4 The production term in the uv equations 13

2.2.5 The length scale equation 13

2.3 The governing equations in approximate form 15 2.3.1 Evaluation of the empirical constants 16

2.3.2 Calculation of bulk velocity . 17

3. THE CALCULATION PROCEDURE 19

3.1 The classification of the system of the governing equations

3.2 Boundary conditions and calculation on the symmetry line

3.3 The numerical scheme

4.

RESULTS AND DISCUSSION .

4.1 The accuracy of the calculation 4.2 Influence of the bulk veloeities

4.2.1 Qualitative observations 4.2.2 The evaluation of fl and f2 4.3 Comparison with experiment

4.3.1 Distribution of U, k, uv and the kinetic energy balance

4.3.2 The flow development

4.3.3 An evaluation of the method

19 21 23 28 28

29

29

31 32 32 34 36

REFERENCES APPENDICES I 11 111 IV V FIGURES

.

.

.

.

.

39

41

43

46

48

51

A dQ,d E(y)

t

f I (n ) f 2 (n ) h k k P L M m m n p P dQ dx Re R . . (y) I I SQ, t U,V,W u,v,W LIST OF SYMBOLS

coefficient matrix

(4x4)

appearing ln Eg. 34

constant appearing in the product ion of the uv eguation ( Eg. 23) a:: _{(C 5 - C4 )}

4-component vector appearing in Eg. 34 dissipation constant (Eg.

6)

diffusion constants (Egs. 22, 23, A.3-5)

constant appearing in the pressure-velocity gradient correlations (Eg. 13)

constant appearing in the pressure-velocity gradient correlations (Eg. 14)

constant appearing ln the product ion term of the uv eguation (Eg. 17)

constant appearing in the dissipation term of the L-eguation (Eg. 20)

nozzle diameter

initial line of errors (Eg. A.4-1)

4-component vector of dependent variables ih Eg. 34 empirical diffusion function for k

empirical diffusion function for uv step size in the cross stream direct ion turbulent kinetic energy, k ::

~

(u2_+v2_+;2) constant in the 2-eguation integral method integral length scale (Eq. 19)

momentum excess

index referrin~ to the cross stream direct ion ln Egs. 37-41

velocity ratio, m :: UI/Uj

index referring to the streamwise direct ion mean component of pressure

fluctuating component of the pressure rate of entrainment

Reynolds number, Re _ (k)1/2 L \!

2-point velocity correlation function, R .. (0)

=

2k

I I

slope of fl(n)

step size in the streamwise direct ion

mean components of the velocity in the x, y, and z directions respectively

fluctuating components of the velocity ln the x, y and z direction respectively

Vk,V_j x,y,z ö ö

n

e

K v 1 o m x

bulk convection velocities for k and uv respectively coordinate directions

effecti7c origin shift due to initial conditions reduced streamwise coordinate ~

=

x - xo

the minimum value of ~ (n )

coordinate along y-axis where the excess velocity lS half the centerline value

central difference operator turbulent energy dissipation

non dimensional cross stream coordinate,

. ( / 2) 1/2

momentum radlus , e:: M pUl

eddy viscosity coefficient (Eq. 1) velocity ratio, À :: UI/Ua

eigenvalues of A

constants appearing 1n formulae

40

and

41

kinematic viscosity

turbulent viscosity density

Subscripts

refers to the uniform outer stream value on the centerline

maximum value

Superscripts

n

=

1. Ö

pr1me refers to differentiation (except 1n Eq. 12) bar over velocity or velocity gradient correlating indicate time average, e.g.

u2 _ lT

fT

_u2 _dt o

The summation convention glves for repeated indices, e.g.

dUo l dx.

l

1. INTRODUCTION

1.1 Literature survey

The problem of turbulence closure takes its roots from the existence of fluctuating velocity correlations (Reynolds

stre~ses) in the time averaged equations of motion. The lack of a constitutive stress-strain relation ln turbulent flows neces-sitates the formulation of a model to describe the Reynolds stresses and obtain a closed set of equations (Refs. 1, 2 , 3 ) .

Earlier works, starting with that of Ref. 4 considered a direct analogy between laminar and turbulent flows by relating the turbulent stresses directly to the mean flow gradients to gl ve, e. g. ,

-puv - pI(

au

-

ay

where K lS an eddy viscosity coefficient to be defined. A com-prehensive summary of attempts in this direct ion is given ln Ref. 1.

All formulations of this type are based on the assump-tion that the scale of turbulence (e .g., the integral length

scale) lS much smaller than the scale of the mean velocity varia-tion. This assumption, however, is incorrect because the scale of turbulence is generally of the same order of magnitude as e.g. the typical boundary layer thickness. An extensive discussion of this is given in Ref. 3.

Integral methods, first developed by von Karman (Ref.

5) constitute another type of approach applicable to turbulent boundary layers. ~hese methods mainly involve the integrated streamwise momentum equation and r equire assumptions relating the boundary layer thicknesses and the wall shear stress. Refer-ences

6

and

7

glve detailed discussions of these methods. The disadvantages with integr al methods are that they predict only the overall properties' of the flow and are liable to breakdown

1n rapidly developing flows (Ref.

8).

The more recent types of turbulence closure models, generally referred to as the field models, use transport equa-tions for the turbulence quantities. These are exact equaequa-tions that are derived from moments of the Navier-Stokes equations

(e.g., Ref. 1). The difficulty encountered in using these is the existence of higher order correlations introduced as new

un-knowns. The pioneering works in the field modelling are those of Chou (Ref.

9)

and Rotta (Ref. 10). The present day popularity of these methods is due to the availability of advanced digital computers which enable numerical solutions of nonlinear partial differentia1 equations.

Since it lS possible to generate transport equations for turbulence quantitites of any order, which contain unknowns of one higher order, it is a necessity to resolve for a level of closure. At the present the optimum level of closure is gene-rally accepted to be of second order. Transport equations for the Reynolds stresses are used and higher order correlations are

approximated in terms of the Reynolds stresses themselves. More-over, empirical inputs are required.

One of the ma1n problems s t i l l unresolved 1n the second order closure techniques is the choice of the necessary and ade-quate equations. Considering the turbulent kinetic energy equa-tion to be the only equation sufficiently documented by experiment Bradshaw et al. (Refs . 11, 12) have traa~formed i t into a rate equation for the shear stress. The basic assumption in this for-mulation is that profiles of all turbulence quantities at a given

streamwise location are related to the shear stress profile by

empirical functions. The outstanding feature lS the modelling

of the diffusion term by bulk convection alone. The shortcoming of this model l S that one would have to use different empirical

functions for different classes of flows, e.g., for boundary layers, duct flows, free shear flows, etc.

The original method(Ref. 12) was applied to three dimensional flows in referenCffi 13 and 14. An adapted version uS1ng the interaction hypothesis to deal with shear layers in which the shear stress changes sign was used to calculate sym-metrical duct flows 1n Ref. 15. Further, the same hypothesis was used in Ref. 16 to calculate two dimensional free shear flows, however, an empirical ~ate equation for the length scale was included in this calculation.

The interaction hypothesis mainly considers the flow as two separate interacting shear layers. These layers are assumed to interact only through the mean velocity profile. As pointed out in Ref. 16, although this approach is plausible 1n two dimensional flows, its applicability_t~the axisymmetric case lS questionable. For axisymmetric flows, one has to consider an infinite number of interactions. A through review of this method, its variations and further applications to flows with extra strain rates is given in Refs. 17 and 18.

The simult aneous use of the shear stress and the turbulent kinetic energy rate equations, which is generally referred to as the technique of two-equation models of tur-bulence seems to be more universal. A given set of empirical inputs is likely to be sufficient to represent a wider range of flows. Computational difficulties due to flow geometry, as found in the interaction hypothesis, are also eliminated.

Turbulence models 1n this direct ion have been de vel-oped and used by Daly and Harlow (Ref. 19), Donaldson (Ref. 20), Hanjalic and Launder (Ref. 21) and Rotta (Ref. 22). All have the common feature that turbulent diffusion is represented as a molecular transport phenomena. In reference 20, the dissipa-tion length scale of turbulence lS assumed to be a constant, whereas a .r.ate equation is used to represent it in Ref. 22. References 19 and 21 use explicit rate equations for the dis-sipation.

All of these methods use the basic formulation of Rotta (Ref. 10) to represent the pressure-velocity gradient correlations appearing in the shear stress transport equation. This term represents one major difficulty in the two-equation models of turbulence due to the impracticability of its measure-ment, and therefore to its direct quantitative evaluation.

References

8,

23, 24 and 25 contain the most recent developments in connection with two-equation models of tur-bulence and their applications to various flow situations.

1.2 Statement of the problem

The present work lS an attempt to investigate the applicability of several physical concepts in turbulence model-ling techniques. Three auxiliary partial differential equations are used to account for the turbulence terms appearing ln the time-averaged equations of motion. The diffusion terms ln the turbulent transport equations are modelled by using the concept of convective diffusion proposed earlier by Bradshaw et al.

(Ref. 12) and later used by Nash (Refs. 13, 14). Morel and

Torda (Ref. 16). The substantial difference of the present work from those is the use of the transport equation for the shear stress.

The resulting set of partial differential equations lS solved directly by an explicit finite difference scheme for the case of an axisymmetric jet ln a coflowing airstream. Physically this flow can be considered as a relatively complex turbulent flow which cannot be put in self-preserving form. This, together with the non zero velocities at the outer edges as weIl as the abseance of solid boundaries make i t an attractive case for testing models of turbulence.

2. THE MATHEMATICAL PROBLEM

In this section, equations that govern the flow are discussed 1n relation to various simplifying assumption. Trans-port equations for the turbulence quantities of interest are given. Assumptions required to close the set of equations are stated and empirical inputs are evaluated according to physical guidelines.

2.1 Governing equations

Equations of motion for incompressible turbulent flows 1n the abscence of body forces are given in Appendix I in

cylindrical coordinates. However, for free turbulent shear

flows, these equations can be simplified based on the following assumptions (Ref. 1):

a. Transverse component of the mean velocity lS much smaller than the streamwise component.

b. Transverse gradients of quantities are much larger than the corresponding streamwise gradients.

c. In the main flow direction, the pressure varles according to the pressure distribution outside the shear layer.

With these assumptions and considering rotational symmetry about the x-aX1S (Fig. 1) to imply W = 0;

:z

= 0, uw = vw = 0,

the equations of motion take the following form for high Reynolds number flows (Appendix I)

Continuity equation ~ + ~= 0 ax ay x-momentum u au +

v

au +

1.

a ax ay y ay y y-momentum uv + -1 dP1 P dx ( 1 )

=

°

where PI and UI are the pressure and velacity af the flaw aut-side the shear layer. Since generally u 2 ~ v2 everywhere in mast flaw fields, as a further appraximatian :x

(u

2

-v

2 )

=

0,

sa that in the case af flaws with canstant 'Outer stream vela-city, the shear stress uv, and the mean velacity campanents U

and Vare the anly unknawns. The immediate prablem therefare, is ta find a farmulatian far uv in 'Order ta clase the set af

equatians 1 and 2.

In the present study, the clasure will be abtained by using transpart equatians ta describe the turbulence. Trans-part equatians far the Reynalds stresses are given in Appendix

11 far a cylindrical caardinate system.

The summatian af the equatians far the different Reynalds narmal stresses yields a scalar equatian far the turbulent kinetic energy. Owing ta the cantinuity canditian, the pressure-velacity gradient carrelatians in the individual Reynalds narmal stress equatians add up ta zera. The eliminatian af these terms, which are the mast difficult terms ta madel, makes the kinetic energy equatian the simplest ane ta use. With the afare mentianed assumptians this equatian glves

U~+V~=

\_-

~

~

canvectian

-

au

- uv -

-ay

~v-···J praduc-tian 1 "a y (k + :Q.) v y oy P \.--.--- - -v- --.---- ..J diffusion di ffu s ian - € v1scaus dissipatian ( 4 )

where k

= ;

(u2+v2_{+w2 ), is the turbulent kinetic energy.}

Hawever, the use af this equatian requires a further assumptian ta relate uv ta k. Radi and Spalding (Ref. 26), Beckwith and Bushnell (Ref. 27), Ng and Spalding (Ref. 28) use the Kalmogarav-Prandtl farmulatian, which gives

uv _{= v T}

au

_ay

where v

This formulation is an improvement over pure eddy viscosity or m1x1ng length assumptions which connect the shear stress directly to the mean flow gradient. However, i t lS s t i l l contradictory to the shear stress transport equation 1n which the shear stress gradients (not the shear stress itself) are related to the

mean flow gradients (Ref. 17).

The first method which uses a transport equation without involving eddy viscosity is that of Bradshaw et al. (Ref. 12). They transform the turbulent kinetic energy trans-port equation (Eq. 4) into a rate equation for uv by using

-uv

=

ak

where a lS an empirical constant.

This assumption lS likely to break down 1n rapidly changing flows and is applicable only to boundary layers (Ref.

17); 1n symmetrical flows like jets and wakes, the turbulent kinetic energy has a high value on the axis of symmetry (e.g. the maximum value 1n measurements of Ref. 29), but the shear

stress is zero.

A more direct way of describing the shear stress lS to use its own transport equation. With the boundary layer type

assumptions mentioned before, this equation is identical for both two dimensional and axisymmetric flows and reads

u

\ ouv + V ouv

=

-

v2 oU ox oy dy . '---Y / \ .... _ y---' convection produc-tion - 2\!(~ ~) dX. dX. 1 1 \ ... - .... ._v _- -.. _..! viscous dissi-pat ion +

E

(~ +

lY)

p dy dX \ .. -.----v .. --... ---.J redistribu-tion d (v2 + E-) u dy p \ -.•... -.--••.•. V .. - - - -.. / diffusion

In this study Eq.

5

lS used to calculate the shear stress, and closure is obtained by approximating the higher order correla-tions in terms of the shear stress and the turbulent kinetic energy. This necessitates the simultaneous use of equations

4

and

5.

2.2 Closure assumptions

To close the set of equations 1, 2, 4 and 5, the

diffusion and viscous dissipation terms 1n Eq.

4,

the

produc-tion, diffusion and redistribution terms in Eq.

5

have to be

modelled. In this section the relevant empirical closure assump-tions are discussed.

To model the dissipation terms i t lS generally

ac-cepted to use the Kolmogorov hypothesis of locally isotropic

turbulence (Refs. 12, 20, 21, 22), i . e . , local isotropy of the

small scale motion will exist if the strain rates of small ed-dies are large compared to that of the mean flow. This lS

ob-tained at high Reynolds numbers, e . g. , Re

₌

(k)1/2 L _\) where L

lS an integral scale of t urbul enc e .

The eddy Slze can be interpreted in terms of the wave number spectrum of the turbulent energy such that i t will be decreasing with increasing wave number. The range of wave num-bers exhibiting local isotropy and including the dissipative e d die s i s r e fe r red t o a s t he" e q u i 1 i b r i um r a n ge" .

In this range, the major parameters are the viscosity, and the amount of energy cascading down the spectrum. The

eddies are in equilibrium because the energy transferred from

the large scale motion is equal to the energy dissipated by

the action of viscosity. Therefore, dissipation can be expressed

as a function of velocity and a length scale of the large scale

mot ion which is responsible for the rate of energy transfer from the mean motion to the turbulent motion.

The above arguments as weIl as dimensional considera-tions lead to E:

=

au.

2 ~(-~)

ax.

J + \) 2

-a

u.u. ~ J

=

ax.

dX. ~ J

=

_L ( 6 )

Also, s~nce ~n isotropic turbulence cross-correlations between velocity components must be identically zero,

2\)(~ ~)

=

0

ax. ax.

~ ~

where L ~s an integral scale of turbulence.

2.2.2 The diffusion terms

Diffusion terms ~n Eqs.

4

and

5

contain both third order velocity correlations and pressure-velocity correlations. It is generally accepted that the latter are identically zero (Refs. 12, 22, 23, 24) so that

E.E=

p ~ p

=

0

According to Townsend (Ref. 30), the existence of turbulent fluctuations over a broad spectrum of frequencies (or wave numbers) produces two different types of turbulent transport. One of these, the gradient diffusion, is due to the effect of the small-scale turbulence on the diffused quan-tity by transportating i t down the mean intensity gradient. In oTder that the gradient should have the same s~gn over a dif-fusing eddy, the scale of the difdif-fusing movements should be small with respect to the scale of variation of the intensity gradients. Mathematically, this can be interpreted as a Taylor series approximation of the diffused quantity retaining only the first order term with the condition that e.g.

where LT is a length pertaining to the scale of turbulence which 1S responsible for gradient diffusion. The turbulent transport by gradient diffusion is therefore conceptually similar to molecular diffusion.

In bulk convection, which is the turbulent transport due to the large scale mot ion, the rate of transport depends on the general distribution of the diffusing quantity. The scale of eddies which are responsible for this type of turbulent transport are comparable with the shear layer thickness.

The total turbulent transport of a quantity ~ will

be given as

( 9 )

where C and C2 are empirical constants, uT is a characteriitic

velocity of turbu~nce and (V~)i is a turbulent transport

veloci-ty which 1S characteristic of large scale motion. The relative

amounts of transport by gradient diffusion and bulk convection

depend on the quantity transported.

The most widely used models (Refs. 20, 21, 22) for the third order velocity correlations effectively assume the turbulent transport to be only of the gradient transport type,

such that in Eq.

9

C2

=

0. However, in the outer region af a

boundary layer flow (or across the whole flow field in a free shear flow) turbulent transport depends on large eddies. These are responsible for the intermittent behaviour of the shear layer edge and are also tied up with most of the turbulent

energy (Ref. 11). Therefore, 1n such flows Eq.

8

will not be

valid, S1nce LT ~ O.

The concept of bulk convection was first used by

Bradshaw et al. (Ref. 12) to model the turbulent kinetic energy

equation. They assumed in Eq. 9 C = 0, so that for ~i = k,

v Eq. 9 gives

u. = 1

(10)

The merits of this choice over the gradient-transport model are (Ref. 31)

a. Since differentiation through data points lS not needed V k lS easier to determine experimentally;

b. Only first order derivatives are involved.

In the present work, applicability of the concept of turbulent transport by bulk convection will be further investigated by applying it to the shear stress as well as the kinetic energy. This gives the following for the third order velocity fluctua-tions : vk

=

C2

V

_k k (11)

,

=

C₂ V. J uv

In Eq. 11 subscripts k and J are not indicei. They refer to the quantities k and uv.

,

C and C a r e constants to be determined from experiment. The choice of V

k and Vj will be explained later.

The pressure-velocity gradient correlations (redis-tribution terms) in the individual normal Reynolds stress equa-tions add up to zero due to continuity. Therefore these terms do not contribute to the balance of the turbulent kinetic energy. Their effect has to be interpreted as a redistribution of energy among different velocity components. Since isotropic turbulence is a statistically more stable state, in non-isotropic turbulence these terms transfer energy from the higher to the lower energy levels to give the flow a "tendency to isotropy" (Ref. 1).

The general equation for these terms was first ob-tained by Chou (Ref.

9)

by taking the divergence of the Reynolds

point

Xo

away from the immediate vicinity of asolid boundary (in the case of a boundary layer) is

U

a

2

u u

,

,

[aum]'[au

i ]] au.

f

(au i )

2(:~~)

dvo~ P ~ 1

(ax~~x:l

(12 )

=

_~ +

p ax. ax. ax~ _aXj x

J J m

\ vOl _{v·- - -- --- - - - I} --v _______ __ -1

P

T PM

where a repeated index indicates summation; ~he terms without a superscript are evaluated at xo, and those with a superscript

'prime' at (XO+x).

Equation (12) shows that the pressure-velocity gradient correlations originate from two different physical mechanisms. The first part of the integral (P

T) consists of interactions of the velocity fluctuations between themselves and the second part (P_M) the interaction of the mean velocity gradients with velocity fluctuations. Assuming that P

T represents the tendency to isotropy in a sheared flow, Rotta (Ref. 10) has proposed the following form

where 0 .. ~s the Kronecker delta.

~J

The second part has been subject to controversy. Donaldson (Ref. 20), Mellor and Herring (Ref: 32) simply ne-glected it. Following Chou(Ref. 9), Hanjalic and Launder (Ref. 21) have assumed

au mi

P

M

=

(-~) ax _m a _~j

m~

where a . ~s a fourth order tensor which must satisfy the

condi-~J

tions of symm~try and of continuity. These conditions suffice

m~

to express a . ~n terms of combinations of the Reynolds stresses.

For flows with a mean velocity U

=

U(y), Rotta (Refs. 22, 23) has made an estimate of this term by developing the mean veloc~y in Taylor series. Retaining only the first term of the series (which accounts for about

60%

of the total) he pr'opo ses

It should be noted that in thin shear flows the model proposed in Ref. 21 (Eq. 14) reduces essentially to the form given by Eq. 15.

In the present study Eqs. 13 and 15 are used to model the redistribution term 1n Eq.

5

to give

l? (~

p ay +

~)

ax

=

L uv

(16)

where _{C 3} and C4 are constants to be determined experimentally.

Following Ref. 21 the normal stress terms will be assumed to be proportional to k. This gives

?

v-

au

ay

au

=

Cs k ay

where Cs lS a constant to be determined experimentally.

The closure assumptions glven by equations

6

and

16

rèqu1re an estimate of the dissipation length scale, L. In this study it is determined by a rate equation. Length scale rate equations can be obtained from the equations of motion by dif-ferentiation or integration with respect to a space variable (Ref. 31).

Differentiation yields a mean-square vorticity trans-port equation which can be related to the dissipation at high Reynolds numbers. This equation can be used directly to account for E appearing in the turbulent kinetic energy equation. Also

L, which will s t i l l be necessary to model the higher order cor-relations can be obtained from

L cc

3/2

(k)

E

(18)

Rodi (Ref. 33) has discussed the E equation. He pointed

out the necessity of grouping the generation by vortex stretching (which is equivalent to the spectral energy transfer to the

smallest eddies in the wave number space) and the viscous dis-sipation terms to obtain terms of the same order of magnitude ln the final form of the equation. However, all the contributions to these terms come from the higher wave number range. Therefore, the use of these quantities and consequently the E equation

itself te account for the length sc ale of large scale mot ion can be considered to be ill-conditioned (Ref. 32).

Integration applied to the rate equation for the two

-point velocity correlation u . A u. B gives transport equations 1 , 1 ,

for an integral length scale times u.u .. Rotta (Ref. 10) has 1 J

derived an equation for an integral length scale defined as

+ <X> 2kL

=

î

f

R .. (y) dy ₁₁ - <X> where R .. lS the 11 2k

=

R .. ( 0). The 11 lateral integral

two-point velocity correlation function and

3

factor

8

makes L to be identical with the scale in an isotropic turbulence field. With the high Reynolds number assumption and the thin shear flow approximations Rotta proposes a rate equation for the quantity kL (Refs. 22 and 23). In this equation the production, diffusion and dissipation terms have to be modelled (Appendix 111); the

. .

.

au.

.

In the present study the following approximations are made to put the kL equation in a tractable form :

a) only the first term in the series expansion for the product ion term is retained;

b) the dissipation term lS assumed to be proportional to the dissipation term in the kinetic energy transport equation; c) the diffusion term is modelled as convective diffusion. With these condideraiions Rotta's equation reduces to

DkL

Dt

= UakL + V akL = uv L au ay C" 2 Dk L Dt' ax ay an Noting that equat ion for

Y

the left-hand side consists _o f k DL _Dt +

DL can be obtained by subtracting L Dk

Dt Dt

from Eq. 20 and dividing by k (Appendix 111).

According to Bradshaw (Ref. 31) _{, ay}

iQ

which represents the effect of the mean strain rate on L can be neglected in

. ()1/2/ . . .

compar1son t o k L, Wh1Ch lS the effect of the fluctuat1ng strain rate. With these considerations and the requirement that V

L

=

Vk' the following rate equation for L is obtained

U aL _ax + V aL _ay

=

_- C ₂ V _{k ay -}aL

c

₁ (C _S-l ) (k) 1/2 ( 21 )

In this equation Cs lS a constant to be determined experimentally.

2.3 The governing equations 1n approximate form

Substitution of the closure assumptions into Eqs.

4

and

5

yields the following approximate equations for the turbu-lent kinetic energy and the shear stress

ak ak au C2 a (k)3/2 u + V =

-

uv _{Y Vkk}

-

_Cl (22) ax ay ay y ay L auv auv (C S-C4)k au

,

avJuv

₁₁

2 _uv ( 23 ) u + V =

-

C₂ CI C3(k) ax ay ay ay L

Equations (1), (2), (21), (22) and (23) now form a closed set of partial differential equations in which

x,y are the independent variables

U,V, uv, k, L are the dependent variables I

Cl' C_{2 ,}

Cz,

Ca, CLp Cs, Ca are constants to be determined

V

j , Vk are empirical informations to be given.

The primary objective in assigning values to the

constants appearlng 1n the approximate forms of the govern1ng

equations lS that they should be applicable to a large number

of different flows. Therefore simple types of flows are used to evaluate some of these constants.

Rotta (Ref.

6)

has noted that Uberoi's measurements

(Ref.

34)

of the decay of an isotropie field in the abscence of mean shear indicate that 2.6 < C3 < 3.0. A value of 2.8 was

used in Ref. 21, the same value is recommended also in Ref. 23.

Hanjalic and Launder (Ref. 21) have considered a

homo-geneous turbulent field in whichthe only non-zero velocity com-ponent U increases linearly with y. For this flow (since spatial

derivatives of time averaged turbulent quantities are identically

zero in a homogeneous turbulent field) equations (22) and (23)

glve

-

uv lQ

=

Cl (k)3/2 ay L ( 24) -(C S-C4) k lQ

=

_{CI C3 (k)} 1/2 uv ay L

Equations 24 and 25 glve

- 2

(CS- C4)

=

_{C 3} (uv) _k '" 0.1 C 3

For the dissipation constant Rotta (Ref. 10) has

estimated Cl

=

0.228 but his later calculations suggest a value

of 0.165. Further from his inductive treatment of the length

scale equation a value of Cs

=

0.8 is obtained.

Making use of the above discussion, ln the present study the following values are assigned to the constants

Cl

=

0.165 C3

=

2.8

(Cs-C4 ) '" 0.1 C3

Cs '" 0.8

I

C_{2 and} C_{2 are the only constants to be adjusted by computer}

optimization.

The bulk velocities for k and uv can be estimated by considering V

k and Vj to be simple functions of x, y (Ref. 12)

such that V k

=

C2 (uv ) 1/2 f 1 (y) m ) 1/2 ( 27 ) I (uv f 2 (y)

V.

₌

C₂ J m

In these expressions,

1

-

uv

11/2

- WhlCh

.

lS

.

the square

m

root of maXlmum uv - is chosen to be a velocity scale of the

large eddies . The functions fl(y) and f2(y) can be estimated

by correlating data for third order and second order velocity correlations. In the present study these functions have been

estimated from the data of Wygnanski and Fiedler (Ref . 29) as

shown in Fig. 2. It should be noted that the symmetry conditions for the distributions of the third and second order correlations

r-e-<fl:l: i:r e f 1 ( :y ) Y tl.

q

I.A ,,'12.

ba\:\.,

be ~e antisymmetFiQ aQd f~(y) to ~~ symmetFie.

- 1/2

The justification for the use of (uv) _m as a velo-city scale for the bulk convection can be made by an analysis analogous to that of Ref. 29. At the outer edge of the flow,

equations 22 and 23 reduce to "convection equals diffusion", so t hat

ok ok C2 Vk k 0

U + V =

-

_C2

-

V_k k

ox

oy y oy (28)

oUv

_V ouv

,

oV.uv J

U + = _C2

ox

oy oy

n

=

Y...

x

Considering that as y ~

00;

k ~ 0 and introducing these equations can be put into the following form

o

k (V-nU) - k Cl (V-nU) Cln on

=

-C2

~

on V k k

o

uv(V-n U ) -Cln uv

o

-(V-nU)

=

on , 0 -C V. 2 on J uv ( 30) ( 31 )

Equations 30 and 31 can be integrated over a small

range of n where uv and k ~ 0 to give

( nU -V )

₌

_{C2 Vk}

( n U-V)

₌

_{C2 V.}

,

J

At the edge of

rate of entrainment dQ _dx relation between dQ and

dx 10 the jet Bradshaw (uv ) to m

(nU-v) will be equal to the

has found an empirical

cor-glve

to correlate the bulk velocities and uv through Eqs. 32.

m

Another implication of Eq. 32 is that V_k is proportional to V. at the outer edges of the flow.

3. THE CALCULATION PROCEDURE

In this section the set of the governing equations is examined qualitatively from the standpoint of classification. An explicit numerical scheme is used to discretize the equations. Boundary conditions as well as conditions for convergence and stability are discussed.

3.1 The classification of the system of the governing equations

Every system of partial differential equations can be put into the following form (Refs . 35, 36)

-+ ~+ ay n L: i=l -+ -+ A.(f,x,y) 1 -+ -+ -+ -+ af ax. 1 -+ -+ -+ + B(f,x;y)= 0 '"2.le

where f, x, Bare vectors and A. are matrices. The system of 1

Eqs. 2, 21, 22, 23 with i = 1 can therefore be written as

-+ af -+ + A(f,x,y) ay -+ af -+-+ ax + B(f,x,y) = 0 ( 34)

where x and y refer to the coordinate directions shown in Fig.

-+ -+

1; f and Bare four-component vectors and A is a 4x4 matrix given by -+ f = -+ B

=

U k uv L 1 - x U V 1 uv A = x (Cs-C₄ )k U uv

avo

1/ ~ + C C (k) 2uv ay 1 3 L 1/2 CdCs-l)(k) 0 0 1 0 (V+C₂V k ) 0 0 (V+C;V.) 0 0 J 0 0 (V+C₂V k )

In any system of equations with the boundary layer type approximations, there will be at least one vertical characteristic (Ref. 12). In the present case, the equation along the vertical characteristic can be obtained by substi-tuting the continuity equation into the x-momentum equation and reads dV + dy 1 dU) U dy 1 Uy d(~ y)

=

0 dy

In Eq. 34 Slnce V appears only 1n the coefficient matrix A, it can be calculated separately from Eq. 35.

From the theory of characteristics (Refs. 36, 37), i t lS known that the type of the system of equations 34 will be determined according to the following matrix equation

(36 )

where I is the identity matrix and AT is the transpose of A. Equation 35 will define a hyperbolic system if the roots

À, ••• , À are real, it will de fine an elliptic system if all the

1 4

roots are complex and i t will define a hybrid elliptic-hyper-bolic system i f some roots are real and some roots are complex .

. In the present case the characteristic equation of the matrix AT, with a

=

Cs - C4' is

The roots of Eq. 37 (or the eigenvalues of the matrix A) are

À = (V+C₂V

k )/U (double characteristic) 1 ,2

(2V+C;V_J 2 1/2

J

À = ± ((C;V

J ) +4ak) /2U 3) 4

and 14ak" are always positive, all the

eharaeteristie direetions.

Under this eondition the system of Eqs. 34 lS eompletely

hyper-bolie.

The eharaeteristic directions reduce to the following

forms for the limiting cases of y ~ 0 (axis of symmetry) and

y ~ 00 (outer edge) (a) y ~ 0 À

=

0 1 , 2 À

=

± (ak) 1 2 / U 3 , 4

The double characteristie lS coincident with the symmetry line and the two others are inelined to the symmetry line at opposite angle s. k ~ 0 À 1 , 2 À 3,4

Here one characteristie coineides with the mean streamline and all the others are inclined at steeper angles due to the outward diffusion.

3.2 Boundary eonditions and calculation on the symmetry line

Equations 34 will de fine the particular problem only

when the initial and boundary conditions are specified. Sinee the problem under investigation involves a symmetrical flow field, the line of symmetry is taken as the lnner boundary.

The numerical seheme employed in this study lS an ~xplicit

outer edge lS not necessary as in implicit schemes (e.g., Ref.

23). However~ Slnce the equations are solved in the physical

plane, an allowance must be made for the spread of the jet. Following Ref. 12, this is done by adding two new points in the transverse direct ion at each x-station. The values of the unknowns at these points are made to comply with the boundary

conditions specified at the outer edge. The edge of the jet is

taken as that point where the difference between two successive values of the mean velocity is less than a certain specified amount (Fig. 3)

The boundary conditions at the inner boundary are prescribed ak aL - , - uv ay ay' conditions d· t . .

au

accor lng 0 the symmetry condltlons such that ay'

are set to zero. At the outer edge the free stream are assumed to be U = Ul' k = 0, uv = 0 and L is calculated by linear extrapolation. The problem is completely defined once the initial conditions (profiles) are glven.

Due to the circular symmetry of the flow, there are -+

two terms in the vector B which will assume the indeterminate

form % at the inner boundary where y

=

o.

Limiting forms for

these terms can be evaluated by McLaurin series expansion,

such that

~(y)

=

uv(O) + y

:~v

(0) + ... Since uv(O)

=

0 lim

~(y)

=

y-+O Y Also auv ay ( 0) lim Vk(y)k(y) y-+O y aV k

=

k(O) ay (0) aV k ay (0) + ...

It has already been stated 1n part (2.3.2) that V

k lS an anti-symmetrical function of y. This adds the constraint that Vk(O)

=

O.

3.3 The numerical scheme

In section 3.1 i t was stated that the hyperbolicity of a system of equations depends on the existence of real characteristic directions. The equations can be reduced to or-dinary differential equations along these characteristic direc-tions. A convenient method, therefore, to solve hyperbolic equations is to make the appropriate choice of the coordinate system and i~tegrate the equations along the characteristic directions. This method, which is referred to as the method of characteristics, provides an accurate process of solving hyper-bolic partial differential equations especially when the initial data are discontinuous (Refs. 36, 38).

Finite difference methods, as opposed to the method of characteristics, offer greater flexibility in organizing the computations for evaluation on a digital computer. If relevant criteria for stability and convergence are obeyed,

the finite difference methods should be expected to glve results of comparable accuracy to the method of characteristics.

Finite difference schemes applicable to problems ln which there is a marching direct ion (x-direction in the present case) can be broadly classified as implicit and explicit schemes. An implicit scheme involves more than one grid point at the

advanced marching level and therefore the calculation of an unknown grid point value necessitates the solution of a set

of simultaneous equations. An explicit scheme, however, expresses one grid point at the advanced marching level in terms of known values atthe previous level. Although computationally easier, explicit methods for hyperbolic systems are prone to stability problems. Stability can be obtained only at the cost of small

A convenient class of explicit methods used for hyper-bol ic syst ems is t he t wo st ep Lax-Wendro ff type met hod. "Any method which can be interpreted as a second order Taylor series expanSl.on l.n terms of the time-like variabIe" is customarily referred to as a method of the Lax-Wendroff type (Ref. 39). This method mainly consists of using Lax's method to calculate the provisional values of the unknowns at the intermediate marching level. These values are then used to calculate the primary points by employing a midpoint leapfrog calculation

(Ref. 39). A typical computational molecule is shown in Fig. 4 . Using the same notation as in Fig.

4

such that the index m

refers to the cross stream direct ion and the index n refers to the streamwise (marching) direct ion, the two steps of this scheme can be written as·

n+l/2 f

=

iÏÎ+l/2 n+l f_

₌

m 1 n

n

_t af Ïi

"2(

fÏii.+l + f_} + _m -₂ (h}m+l/2 TI _{af n+l/2} f_ + t

(-a)-m x m

a

A formula relating the differential operator D

=

ax to the central difference operator 8 is given by ~ef. 40)

x D

₌

1 ( 8 h where

n

8 f_ = x m

n

ö3 f_

=

x m etc. 1 1 8 3 + l ' 3 2 8 5

-

...

) x _22'3! x 24'5! x

n

-n fiii+l / 2 fiii - I / 2 Ïl

n

n

f_{iii +l / 2 -} 3f

_m

_{+l / 2} + 3f

_m-

_{l / 2} -*

For the sake of simplicity of notation the vector f will be written without the overbar.

The original two-step Lax-Wendroff scheme employs only second order correct finite difference operators. Nash

(Ref. 14), however, reports the occurrence of instabilities ln his calculation of a three dimensional boundary layer using Bradshaw's turbulence model and employing the Lax-Wendroff

scheme. Stability was restored by using fourth order accurate central differences infue cross stream direction. He further reports that maximum precision was obtained when fourth order differences were used only in t he first step. Maximum stability was obtained when the prlmary points in the advanced marching level were calculated from the intermediate level only. Making use of formulas 37, 38, 39 and employing fourth order central differences in the first step, Eq. 34 can be put into the follo-wing difference form

n+l/2 fiii+l/2

=

and

n

n

n

;h Am+l/2 (f m+l - f m) +

n

n

n

n

n

~n + t Q A ( f 3f + 3 f f ) - -L B 2h

24

m+l/2 m+2 - m+l

m -

m-l 2 m+l/2 ( 40)

o

Ql n+l/2 n+l/2 n+l f_ m

=

(l-Ql) f

m

+ :2(fm+l / 2 - f m- l / 2 ) t n+l/2 0+1/2 0+1/2 n+l/2 (fm+l / 2 f

_m-

_{l / 2 )}

....

( 41 ) A_

-

-

t B_ h m m

In equations 40 and 41, Q and Ql are disposable con-stant s such that when Q

=

Ql

=

0, equations 40 and 41 will reduce to the original Lax-Wendroff type equations.

The hal f- spac e values (at ïii+l/2 and n+l/2) of A and

-

B can be evaluated by averaging values at data points such that n

1 n

n

n+l/2 1 fi+l/2 n+l/2

A =

(L

l+A_) A = _{(\ü+l/2+ Am-l/2)}

-+

and similarly for B.

It lS known that explicit finite difference

discreti-1

zations of hyperbolic equations will converge only if the numerical domain of dependence includes the exact domain of

.

5

l·t . . pn+l

·11

b

dependence. In F1g. lS assumed that pOlnt _{_ _m} _ Wl e calculated by using data !rom points

P~+l

and

P~-l

Then, the

. n + l . .

reg10n of dependence of p_ can be def1ned as the reglon on m

the data line which is bounded by the outermost characteristics through

P~.

If

'mI

and 'ÏÏ' are small the characteristics through

Yr+l m . . . .

p_ ean be approx1mated by stralght llnes, and the

characterlS-~ . . ÏÏ fi+l

tlC d1rectlons through p_ can be used to represent those of P_

m . m

Convergence will be insured only if (e.g. Refs.

36, 37)

This condition, which is referred to as the Courant-Friedrichs-Lewy condition (CFL-condition), can be shown to be coincident with the stability requirements for Eqs.

40, 41,

e.g. by using the von Neumann stability analysis (Appendix IV). The CFL-con-dition for stability is necessary but not sufficient; however, i t must be obeyed as a criterion for choosing the step size in the marching direction.

The accuracy of a finite difference scheme can be estimated from the local discretization error. In Appendix V i t lS shown that the explicit formulation of Eq.

34

glven by Eqs.

40

and

41

is second order accurate in both x and y.

The V-velocity is calculated from Eq. 35 by uSlng a simple predictor-corrector scheme so that

1

A finite difference scheme lS convergent if at a fixed point p(x,y) the difference between the exact solution and the

n

finite difference solution, Z_(x,y) + 0 as

m,n

+ 0 (e.g., Ref. m

dV

=

V (y.) + f:,y (y . )

1 dy 1

( 43 )

4. RESULTS AND DISCUSSION

In this section varlOUS factors influencing the accuracy of the calculations are discussed. Effects of the diffusion func-tions on the solution are investigated. The numerical results are compared with experimental results and finally a critical evaluation of the method is made.

4.1 The accuracy of the calculation

The finite difference form of equations 34, which lS glven by Eq. 40 and 41 was programmed in Fortran IV and the calculations were performed on a Mitra 15 computer. The details of the program are given in Ref. 41.

The accuracy of calculations was checked by changing the step Slze ln the cross-stream direction. No large differences were observed ln the results due to this change. As another

necessary check, the excess momentum of the jet was calculated at each station by integrating the calculated mean velocity profile by using a simple trapezoidal rule. The momentum was

conserved within a few per cent over a distance of about 200 ~.

This error was assumed to be negligible so that no special

techniques were used to correct this momentum unbalance.

Input profiles of the unknowns were glven with constant step Slze in the y-direction with a typical value of 0.10 which

corresponds approximately to 25 points on the profiles.

It lS weIl known that the solutions of hyperbolic equations depend on the initial data. This necessitates the pro-vision of accurate experimental data for all the unknown quan-tities .. In the present case the inputs were obtained from the measurements reported in Ref. 42. The lateral integral length

scale was assumed to be proportional to the longitudinal integra1 length scale. Although this cannot be physically justified, it

propagation of irregularities in the distributions of the measured quantities into the computation domain, all the

pro-files were smoothed. To check the dependence of the solution on the particular input conditions, calculations were repeated using two different inputs, at ~

=

10 and 20.

The step S1ze 1n the x-direct ion is g1ven by the step S1ze 1n the y-direction divided by the largest character-istic angle. Owing to the non-zero velocities at the outer edges (almost everywhere UI ~ O.lUO) the characteristic angles were always less than one. This ensured relatively large step sizes in the marching direction. Further, to ensure stability the step S1ze in the marc hing direct ion was multiplied by a safety factor smaller than one.

4.2 Influence of the bulk velocities

An overall estimate of the functions fl and f2, which describe the bulk velocities of k and uv respectively (Eq. 2T)

was made according to their distributions obtained by correlating the data of Ref. 29. These functions were both prescribed in

terms of n

=

Z.

The following additional conditions had to be ê

considered for numerical and physical consistency:

a) fl(n) must be zero at the centerline to meet the condition of symmetry;

b) at the outer edge, the characteristics have to preserve their respective centerline orientations such that

(V+C~V.)/u J

will be the outermost characteristic. This condition 1S necessary to avoid intersection of characteristics of the same family, e.g. of

è

and ~

=

(V+C_{2 Vk )/U.}

c) the data of Ref. 29 show that the diffusion term in the uv equation is zero on the centerline. This term is given by

auv 2 a ( uv avo auv)

=

c'

uv V.

₌

c'

-2. + V. ay 2 ay J 2 ay J ay

Since at the centerline 0, auv

"

0, to have diffusion

uv

=

zero

ay

v./

_{J y=O}

=

0 or _f(O)

₌

°

The edge of the jet coincides with the outermost characteristic. Therefore the value of À and consequently the

4 val ue 0 f V.

J is the decisive factor in the diffusivity of the turbulent field.

Numerous trial and error modifications to the shapes of f i and f 2 finally yielded the following :

(a) In order to obtain increasing kinetic energy level with the streamwise coordinate on the centerlin~ fi had to attain negative values. Although this lS not a direct outcome of the correlation measurements of Ref. 29, i t is compatible with the basic trans-port principle which says that energy should be transtrans-ported away from regions of maximum energy.

(b) The zero crossing of fi was found to have a significant effect on the calculated profiles of the kinetic energy and the mean velocity. As pointed out in Ref. 16 a sufficient zero

crossing was found to coincide with the position of maximum shear.

(c) The shape of the shear stress profile was found to be lnsen-sitive to f2. However, near the edges of the flow it had to be adjusted to allow for the rate of spread of the jet.

Considering all this, the functions shown ln Fig.

6

are found to be adequate to describe the diffusion of kinetic energy and the shear stress.

The calculations were performed by setting the

fun-damenta~ constantsXto their values specified in section 2.3.1

and adjusting the diffusion functions for the case of m

=

0.1.

.

,

For th1S C2 and C2 were set equal to one so that the only con-stants to be optimized were yO and the slope St (Fig.

6).

Although order of magnitude values for these constants could be obtained from the data of Ref. 29, the precise evaluation had to be made by numerical experimentation. The value of yO was found to have a very significant effect on the distributions of the kinetic energy. The mean velocity and uv profiles were observed to be much less sensitive to yO' The results of three different calculations with Yo = -0.15, Yo = -0.20 and yO = -0.30 are shown in Fig.

7.

In these calculations the value of all the other constants were kept the same. For Yo

=

-0.30, the kinetic energy profile produces an overshoot towards the

center-x

line which increases with

d .

For YO

=

-0.15, an opposite trend to the above is observed. In this case, the value of the kinetic energy on the centerline decreased with downstream distance up

x

to about 50

d'

For yo = -0.20, the distributions of the kinetic energy retained their original shapes; moreover, the centerline values of the kinetic energy were observed to increase, as

expected, as early as a few diameters in front of the initial data line. Considering these, Yo =-0.20 was chosen as an ade-quate value.

This type of dependence of the solution on the diffu-S10n functions (or constants) ha s also been observed by Rotta (Ref. 22). In his calculations for the plane wake flow, diffu-sion was modelled as gradient transport type (section 2.2.2). For certain values of the diffusion constant, Rotta reports the occurrence of saddle-shape4 mean velocity profiles as well as significant decrease of the kinetic energy on tpe centerline.

~

These constants which are not related to the diffusion will be referred to as the fundament al constant.

The value of s2 was adjusted to obtain an acceptable variation of 0 with x. An optimal value for this was found to be 0.9. Derivative discontinuities on fl at n

=

1 due to dif-ferent slopes at

n

< 1 and

n

> 1 (Fig. 6) had no apparent effect on any of the profiles.

The same values for yO and s~ which were selected for the case of m

=

0.1 were used to calculate the other two cases, namely, m

=

0.2 and m

=

0.3.

4.3 Comparison with experiment

In this section the results of the computations are compared with experiments reported in Ref. 42.

4.3.1 Qi~t~i~u~i~n~ ~f_U~ ~,_u~ and ~h~ ~i~e~i~ ~n~rgy_b~l~n~e

Profiles of mean excess velocity, kinetic energy and shear stress across the jet are shown in Figs. 8a to 16a. For

-

x

all the velocity ratios and the whole range of ~' mean velocity profiles attain self similarity. Agreement between predictions and the experimental distributions of the excess mean velocity are very close.

Calculated profiles of the Reynolds shear stress are also ln very good agreement with the experimental distributions. For almost all the cases, the maximum shear as weIl as its

location are weIl predicted. The only notabIe discrepancies are at

~

= 53 for m = 0.2 and Rd = 65 for m = 0.3. The larger

~ ~

values of the predicted maxima suggest a quicker rise in the

-x

calculated shear stress level with

d'

than indicated by experiment.

o

For m

=

0.1 (Figs. 8a, 9a and 10a), the calculated kinetic energy, k, profiles are generally in good agreement with experiment. CenteEline values are closely predicted. How-ever, with increasing ~, the calculated maxima become higher

o

energy to diffuse towards the centerline with increasing ~ is

d~

displayed.

For m

=

0.2 (Figs. lla, 12a and 13a) the centerline values of k are reasonably weIl predicted all along the flow. The difference between the predicted and experimental values, however, increase with ~. Also, the locations of the maxima

do shift towards the centerline.

The predicted centerline values of k for m

=

0.3 (Figs. 14, 15 and 16) are considerably higher than the experi-mental results. The locations of the maxima show a marked devia-

-~

tion towards the centerline with increasing ~ . For ~

=

182

do d~

(Fig. 16a), the shape of the calculated distribution lS con-siderably different than the experimental profile:

In figures 8b to 16b, the calculated kinetic energy balances are presented. As a means for comparison, the energy balance for a round jet (Ref. 29) and for a round wake (Ref. 43) are given in Figs. 17 and 18. The most pronounced differ-ence of the computed energy balance from those for the jet and the wake is the convection term. The computed convection term has a smaller contribution to the energy balance. It should be noted that this term is given by

u

~ + V

ak

ax

ay

so that the transversal distribution of k will have a large effect on its convection. In the data of Ref. 29 the location of maX1mum kinetic energy coincides with the axis of symmetry

ak

so that

_ay

is always negative. The experimental results reported 1n Ref. 42 which have been used as the initial data in the

present calculations display regions of positive

~~

around the

axis of symmetry.

In almost all the cases the requirement that "diffu-sion should integrate to zero" is closely satisfied. The ratio

of production to dissipation of energy shows a slight decrease with increasing UI/Ua. However, even at the large st UI/Ua (cor-responding to m = 0.3, ~ = 182) this ratio stays much larger

d

than that found in the round wake.

In figur3s 19 - 22 the flow development is investigated variations of

d~

(Fig. 19) is reasonably For m

=

0.1 and m

=

0.2, experiment al

. x T

as a funct~on of ~. he weIl predicted for all m.

cS

x

results indicate a more linear dependenee of ~ on ~.

The agreement between the experimental and computed

_,.

variations of the mean velocity decay with

:0

is fairly good

for m

=

0.1. The discrepancy between predictions and experimental results increases with increasing m. This could be

partly attributed to the momentum unbalance present in the experiments to have pronounced effects at large m.

The predicted variations of

uVfu/U~

are ~n close agreement with experiment. Only

with x (Fig. 21) dO

for m

=

0.3, the increase in the computed value is more rapid than indicated by experiment. For

x

d > 60 the predicted uv for m

=

0.3 decreases

f<# 0 m

slightly with d~ towards an asymptotic value ~f 0.03. The agree-ment between the predicted variations of ka/Ua are generally in f~ agreement with experiment (Fig. 22). For m

=

0.3,

0/°0 _{the pred~cted} . _{values are h~gher} . t _{han ~nd~cated} . . _by exper-iment

In most works on fully developed free turbulent flows it ~s generally assumed that some distance downstream of the nozzle the flow development will depend only on the variables,

(x-xO)/8. 8 is the excess momentum of the jet given by

where M _

00

2wp

J

U(U-UI) dy

o

the initial conditions at the nozzle are assumed to affect only xO' which is the location of the virtual origin.

In figures 23 -26 the flow development lS investigated

as a function of (x-xo)/6. The predictions are compared with an integral method using a two-equation model of turbulence

(Ref.

44).

This method uses Rotta's gradient transport type

formulation for shear diffusion. In reference

44

it was found

to give~best agreement with

o

·th x h th d·ff .

e

Wl

e

w en e l USlon

experimental variations of UI/Ua and

function kp was set to 1.

The mean flow development (Fig. 23 and 24) lS predicted

reasonably well by both the present finite difference computa-tions (FDC) and the integral method. It should be noted that the experiment as well

haviour of tand UI/Ua

as the predictions show a universal be-~

with x

6

Maximum values of uv and centerline values of k are

.

it

shown as a functlon of

e

ln Figs. 25 and 26. The experimental

r~sul~s dis~lay a non-~niversal behaviour for both uV_m and k

wlth - for - < 20 and - < 30 respectively. This is reproduced

e

6 6 ~

x

in the results of the FDC as well. For large

e

a tendency to a

universal behaviour is observed ln both the experimental results

and the results of the FDC. For m

=

0.3, the predicted and

experimental asymptotic values of uv are very close. The results

m ~

x

of the integral method are erroneous especially at large 6

Figure 27 shows the variations of the centerline

H

of

d~

with

d~

.

d~

varies as astrong function of :0

~nd

is

dependent on m. Figure 28 shows centerline values of

6

as a

.

ie

. .

L . f

values also

func-tlon of 6 The expected lncrease ln

6

lS observed only or

m

=

0.3 for

~

> 30. This anomalous behaviour of

~

should be

regarded as an explanation for the sudden increase in uv and k

3t

(Figs. 23, 24) at - ~ 30.