The prediction of turbulent inclined flow in rod bundles

Computers & Fluids Vol. 14, No. 4, pp. 361-378, 1986

0045-7930/86 $3.00+ 0.00 Pergamon Journals Ltd

Printed in Great Britain

THE PREDICTION OF TURBULENT INCLINED FLOW IN ROD BUNDLES K. A. ANTONOPOULOS National Technical University of Athens, Mechanical Engineering Department, Thermal Section, 42 Patission Street, Athens 106 82, Greece (Received 1 November 1984;/n revised form 19 November 1985) AImtract--The present work describes a numerical procedure for the prediction of turbulent, fully developed inclined flow in rod bundles, including the limiting cases of purely axial and purely transverse flow. The procedure is based on a finite-difference solution of the momentum and continuity differential equations expressed in curvilinear--orthogonal co-ordinates. The effects of turbulence are simulated by the "k-~" turbulence model. Application of the procedure provides new information about the behaviour of turbulent inclined flow. The effects of inclination angle on the pressure-drop and on the other characteristics of the flow are quantified.

NOMENCLATURE A c

a~

c,, c~, c, c~ c,~ cp D Dh E G k P

eL

Pr Pw P

R% Re, Rec Re

SL Sr

So f,

~, V, W

6 Y y*

Greek symbols r °.dr

Ap~ Ap

K

0 p

cross-sectional area of a subehannel coefficient in finite-difference equations turbulence model constants axial pressure drop coefficient, c,. ffi ~J(0.5 p~,2) transverse pressure drop coefficient, c~ ffi ApJ(O.5pa2) pressure drop coefficient in the direction of an inclined flow, c~--Ap/(0.5pf '2) tube or rod diameter hydraulic diameter, D~ = 4(AJPw) turbulence model constant generation of kinetic energy of turbulence kinetic energy of turbulence metric coefficients associated with co-ordinates ~, ~, ~ respectively pitch of a rod bundle, defined as P = Pt = Pr in the case of the in-line square arrangement and as P = [(Pr/2)z+ p~.]o.s in the case of the equilateral triangular or the staggered square arrangement longitudinal pitch, defined as the distance between centre fines of adjacent transverse rows transverse pitch, defined as the centre-to-centre distance between adjacent rods in a transverse row wetted perimeter local pressure axial Reynolds number, Rea = p~D/p axial Reynolds number based on the hydraulic diameter, Re`--(P~Dh)/p transverse Reynolds number, Re~ = (p6D)/~ Reynolds number of an inclined flow, Re = (pFD)/I~ longitudinal pitch-to-diameter ratio, SL = PzJD transverse pitch-to-diameter ratio, S r = PT/D source term in the transport equations mean velocity in the direction of an inclined flow velocity components in the ~, r/and ~ directions respectively mean velocity through the minimum gap mean axial velocity distance of a grid node from the rod surface dimensionless distance of a grid node from the rod surface effective exchange coefficient for property @ pressure drop per row of rods in the transverse direction pressure drop per row of rods in the direction of an inclined flow volumetric rate of dissipation of turbulence kinetic energy rectilinear co-ordinate in the axial direction curvilinear-orthogonal co-ordinates on the cross-sectional plane turbulence model constant angle of inclination of the bulk velocity to the rod axes, O --- aretan (6/g') molecular viscosity effective viscosity fluid density 361

362

K . A . ANTONOPOULOS

0",~. 0", "t'~, To,,'t';~

turbulence model constants normal stresses shear stresses average wall shear stress around the periphery of a rod in the ~-direction general dependent variable which may stand for u, v, w, k, angular position of a point on rod surface measured from the rear of the rod

1. I N T R O D U C T I O N

The present paper is concerned with the numerical prediction of turbulent flow over regular arrays of parallel rods such as those illustrated in Fig. 1(a). Banks of rods or tubes feature in many types of industrial equipment including heat exchangers, boilers, condensers and nuclear reactors. Practical examples of tube-bank flows are usually two- and threedimensional, the simplest class being the fully developed "axial" or "parallel" flow (Fig. l a ) , as in, for example, a nuclear-reactor fuel assembly consisting of rod bundles with the coolant flowing axially in the spaces between the rods. A more complex two-dimensional tube band flow is the "transverse" flow or "cross-flow" (Fig. la) occurring in certain types of heat exchangers in which one stream is directed normal to the tube axes. In certain types of industrial equipment, including heat exchangers with baffles or other obstructions, the flow is inclined to the axes of the tubes (Fig. la). This is the case of primary interest in the present work. Previous numerical studies of turbulent flow in rod bundles have been confined, almost without exception, to the two limiting cases mentioned above, i.e. the purely axial and

'"" J

FLOW f INCUNED FLOW

((:]1

~TRANSVERSE "~

FLOW

'..J

(',

•~--pL ---~

symmetry

(b)

6

~ F~ k.A

"rl

k.A F k. J

.,,k. J l

a rk ) (~

(c) oN

Fig. I. Flow in rod bundles. (a) Axial, transverse and inclined flow. (b) In-line and staggered arrangements. (c) Solution domain corresponding to staggered arrangments. (d) Curvilinearorthogonal computational mesh on the cross-sectional plane.

Turbulent inclinedflowin rod bundles

363

purely transverse flow. The work of Carajlescov and Todreas [i] is typical of the first case and that of Massey [2] and Le Feuvre [3] of the second. In the above studies, the governing differential equations were solved by finite-difference means with turbulence effects simulated by use of turbulence models of various kinds. Unlike the rather limited numerical studies, a large amount of experimental work has been reported but only for the limiting cases of axial and transverse flow. Typical is the Zhukanskas [4] collection of data. The inclined flow received very little attention and in those cases where experimental studies of this kind exist, the data provided refer to very special flow situations. Among the few correlations available for the inclined flow case, typical are those of the ESDU [5], which are based partly on Kazakievitch's [6] experiments and partly on theoretical considerations. The objective here is to provide an economical numerical procedure for predicting turbulent, fully developed inclined flow in rod bundles and to gain an insight into this virtually unexplored flow case by applying this procedure. The limiting cases of purely axial and purely transverse, fully developed flow are also examined for the purpose of evaluating the accuracy of the procedure developed by comparing results with existing data, since, as already mentioned, such data do not exist for the inclined case. The procedure developed is based on a finite difference solution of the continuity and momentum differential equations expressed in curvilinear-orthogonal co-ordinates, which allow the irregularly shaped solution domain to be mapped accurately and economically. Turbulence effects are simulated by the " k - e " turbulence model [7]. 2. CO-ORDINATE SYSTEM AND DIFFERENTIAL EQUATIONS The solution domain is of the general irregular shape shown in the example of Fig. ! (b) (cross-sectional view) or in Fig. l(c) (enlarged isometric view). A general orthogonalcurvilinear co-ordinate frame l - t / i s used to map the cross-sectional plane while the third co-ordinate ~ is taken to be rectilinear and is aligned with the axial direction, as showed in Fig. 1(c). Figure 1(d) shows an example of orthogonal-curvilinear computational mesh of co-ordinate lines, which maps the cross-sectional plane. The generation of the mesh for any particular shape of cross-section is performed by a numerical procedure described in [8]. Associated with the co-ordinates ¢ and r/are the spatially varying metric coefficients 1¢ and l~, which connect increments of ~ and r/to increments of physical distance. The values of l¢ and I~ are calculated by the mesh generation procedure [8]. The metric coefficient l¢, which corresponds to the axial co-ordinate ~, is specified to be unity. The transport equations governing turbulent three-dimensional flow in the co-ordinate system ¢, r/, ~ defined above, can be written in the following common form: I 0

1¢1,0¢ (putO

1 a

(p,,h*)+

(pw,)

o(_°°)

and the continuity equation may be written as ! 63

t,t, 63¢ (a,,t,) +

163

63 +

O,w) = o.

(2)

In the general transport equation, (1), q~ is the dependent variable and may stand for the velocity components u, v, w in the ¢, )7, ( directions respectively, for the kinetic energy of turbulence k and for the volumetric rate of its dissipation t. The latter two quantities pertain to the "k-~" turbulence model, employed to simulate turbulence transport effects, details of which may be found in [7]. The term F,.,~ in the equation is the "effective" exchange coefficient for property ¢, being the effective viscosity/~ for the momentum equations and the effective diffusivities Fk.~, F,.,fr for the k and t transport equations respectively, i.e. [7]

364

K.A.

ANTONOPOULOS

Table 2. Expressions for the source term Sd, in equation (I) S,p

Table I. Turbulence model constants Constant

Value

C. ok

0.09 1.00

o, ic Ci C2 E

K21[(C2

-

I dp - - - - + terms arising from co-ordinate curvature I: ?~ I dp - ~--~ + terms arising from co-ordinate curvature

u

r

C, )CI. 21

.'

0.4187 1.44 i .92 9.793

-~

dp

+ terms arising from co-ordinate curvature

k

G - p~

E

CiG'~-C2p ~

~2

F..ar F....~ F....er # + C, p k21E #ar

(3)

F,.,er = lt,rrlak r.<~ = mdcr..

(4)

=

=

=

=

(5)

where # is the molecular viscosity and C~, ak, a, are empirical coefficients whose values are given in Table 1. Lastly, the term S, in equation (1) stands for the source or sink of the property @, as shown in Table 2. Thus, when ¢ stands for the velocity components, So contains the appropriate pressure gradient term as well as terms arising from the co-ordinate curvature. When @ stands for k or E, So contains generation and dissipation terms as dictated by the turbulence model. The full expressions for S,, derived for each case, may be found in [8]. The term G, contained in the k and E source terms (Table 2) stands for the generation of k and is calculated as 1 , 2 2 G = 2--~ [~¢ + ~.. + r~¢ + 2(r~. + r ~,,. + ~.¢)1,

(6)

where es denote normal and shear stresses. The expressions derived for these in terms of the co-ordinates ~, q, ~ are given below:

(7) /1 ~v

"C,;¢ =

2

u cgl.k

(8)

0w

(9)

/.i¢ff-~-

[1 ,)v

10U

v OI~

t,~,t

td,~

u Ol¢

id.on

[1 t3w Ou) = "°"t,T, o7 + tel dw =

3.

BOUNDARY

(to) (11)

Or) +

J

02)

•

CONDITIONS

The boundaries of the solution domain (Fig. lc) include symmetry planes, inlet and outlet planes, and solid walls. On a symmetry plane, the velocity normal to this plane is zero as are the normal gradients of the remaining velocity components and of all other variables. The boundary conditions at inlet and outlet planes depend on the flow case to be simulated. These boundary conditions and the manner of imposing them, which is one

Turbulent inclined flow in rod bundles

365

of the novelties of the present study, are discussed in Section 5, where various flow cases are examined. The usual "non-slip" conditions apply at a wall. The various diffusive fluxes across the wall are calculated, for reasons of accuracy and economy [8], by use of "wall functions", which are based on the assumptions that in the region very near the wall: (a) the flow behaves as an one-dimensional Couette flow with low streamwise pressure gradient, and (b) production of k is approximately equal to its dissipation. Details about wall functions may be found in [7, 8]. A brief outline of the relations employed is given below. The nearest to the wall grid node P, located at distance y from the wall, is considered to lie in the fully turbulent region of the Couette flow, if the value of quantity y* at point P, defined as

=FpCTk"'Yl L

Y*

~

(13)

Jp

is greater than 11.63, while for y* ~< 11.63 the point P is taken to lie in the laminar sublayer. Thus, for the momentum equations, the wall shear stress component z,.¢ in the -direction at point P is calculated by use of the linear or the logarithmic law of the wall [9] as ~,.~= g

for y* ~< 11.63

(14)

P

l z..e L In(Ty,)"Je

for y*>11.63,

(15)

while the component ~.¢ in the ~-direction is calculated by relations similar to the above with u replaced by w. The values of constants C~, x and E are given in Table 1. For the k-transport equation, the source term Sk = G - p c (Table 2) at the points P nearest to the wall is modified by use of wall functions as follows: the generation term G [equation (6)] is calculated by use of the wall shear stress expressions given above, while the dissipation term pc is calculated as LoE]p = FPC~/4 L Yk3/2y* 1le

for y* ~< 11.63

Loc]e=[PC:/' k3/2ln(EY*)'/e _

for

y* > 11.63.

(16) (17)

For the E-transport equation, the value of c at the nearest to the wall points P is calculated as 314 312

Ep =

L

xy

]p

(18)

4. S O L U T I O N OF D I F F E R E N T I A L E Q U A T I O N S

The solution procedure of the governing differential equations (I) and (2) is based on a finite-volume formulation described in [10], which is suitably modified so as to incorporate the orthogonal-curvilinear co-ordinates. Accordingly, only a brief outline of the solution procedure is given below. Further details are available from [8]. The variables are stored at grid points formed by the intersections of the three families of co-ordinate lines. The velocities are taken to lie midway between the pressures which drive them. The finite-difference equations for each variable are derived by approximate integration of the differential equations over control volumes (Fig. ld) whose faces are formed by co-ordinate surfaces and whose centers are points at which the particular variable is stored. This preserves the conservation property of the parent equations. Such

366

K.A. ANTONOPOULOS

assumptions about the distributions of the variables between nodes are made in the process that the resulting differencing scheme is a hybrid of central/upwind differencing depending on the ratio of convective/diffusive flux coefficients. The resulting difference equations are of the form

n

n = E, W , N , S , D , U,

(19)

where q~ stands for u, v, w, k or c; the coefficients a express the combined effects of convection and diffusion; S~.e is the integrated source term; and the summation is over the neighbours E, W, N, S, D, U of the central node P (Fig. ld) of each control volume. The continuity equation (2) is approximated by central differences and is used in conjunction with the momentum finite-difference equations to derive a set of equations of the general form of (19) for a field of pressure perturbations which tend to drive the velocities towards the satisfaction of continuity. The sets of algebraic equations (19), resulting in the way outlined above, for each variable are solved sequentially and repeatedly. The velocities in each cycle of calculation are obtained from the momentum equations with the prevailing pressure field inserted. They are then adjusted, along with the pressure field, to accord with the solution of the pressure perturbation (i.e. continuity) equation. The equation sets for the variables k and E are then solved. The solution of the equation set for each variable is effected separately by a form of Alternating Direction Implicit (ADI) algorithm in conjunction with block pressure and velocity adjustments [8]. The solution is considered converged when the residual in each of the equations diminishes below a prescribed level.

5. FLOW CASES EXAMINED 5.1 Preliminary remarks

The method outlned in the previous section is in principle capable of calculating tube bank flows for arbitrary circumstances including those which give rise to a three-dimensional structure. In the present section the general equations (l) and (2) are simplified so as to simulate the following flow cases: (a) purely axial, fully developed flow (i.e. inclination angle O = arctan ( ~ / ~ ) = 0°); (b) purely transverse, fully developed flow (O = 90°); (c) fully developed, uniformly inclined flow (0° < O < 90°). The simplifications performed result in a considerable reduction in computer storage and time compared to that for the general case. Extracts are shown from a wide-ranging body of predictions for the cases mentioned above. In each case, laminar flow predictions are also given so as to allow evaluation of the numerical technique without the possible errors arising from the inadequacies of the turbulence model. Although the main interest of the present work lies in case (c), the first two cases were also included mainly for allowing the accuracy of the method to be evaluated by comparison with existing theoretical and experimental data, since no such information is available for the inclined flow case. Details about the grids employed (an example has already been shown in Fig. ld) for the numerical solution and about grid dependence tests may be found in [8]. The most important limitations imposed by the k ~ turbulence model in connection to the present application, are summarized below: (i) It is assumed that the flow within the tube bank is wholly turbulent and therefore transition effects cannot be simulated. The present predictions are therefore expected to be more reliable for the higher Reynolds numbers. This limitation is, however, less restricting in an inclined flow than in a purely transverse one because, for a given transverse Reynolds number, higher turbulence is expected in the former case than in the latter due

Turbulent inclinedfl~,~ in rod bundles

367

to extra turbulence generation associated with the axial velocity gradients in the ~ and r/ directions. (ii) Turbulence-driven secondary motions, which arise in purely axial flow [!] cannot be predicted. The present predictions in this case are, therefore, expected to be less accurate in respect of local quantities (i.e. local values of axial velocity, wall shear stresses and turbulence properties) than in respect of overall quantities (i.e. average drag) where the effects of secondary motions are to some extent averaged out. It should, however, be emphasized that in inclined flows, the transverse component of the flow will swamp the turbulence-driven secondary flows and therefore the k-f model is not precluded from producing realistic predictions of overall as well as of local features. (iii) The assumptions of low streamwise pressure gradient and of local equilibrium (i.e. production of k approximately equal to its dissipation) near the tube wall, which are invoked in deriving the wall functions are quite reasonable for a purely axial flow but less so for the case of transverse flow, where separation (and reattachment in the in-line tube arrangements) of flow occurs at the rod surface in many geometrical configurations.

5.2 Purely axial fully developedflow 5.2.1 Equations and boundary conditions In this case (inclination angle O = 0°), the general equations (1) and (2) are simplified as follows: (a) The condition of full development implies that all derivatives in respect to the axial co-ordinate ¢ vanish with the exception of the axial pressure gradient term ap/a~. (b) Transverse velocities u and v are set to zero. The problem is therefore described by the following equations: Transport equations for ~ = w, k, ~:

Since the axial co-ordinate ~ does not appear in the above equations (except for the axial pressure gradient term Op/O~contained in the source term of the w-momentum equation) the solution can be confined to a single cross-sectional plane ~-t 1with the term 0p/O~ either prescribed as an input or adjusted during the course of the calculation so that the desired mass flow rate may be obtained. The wall functions are used for imposing the boundary conditions on the rod surfaces AF and CD (Fig. 2a) while on all other boundaries of the solution domain, which are planes of symmetry, the normal gradients of w, k and c are set to zero.

5.2.2 Results Calculations for various rod arrangements have bcen performed for both laminar and turbulent flow. Comparisons of the results with existing theoretical and experimental data have shown that in the present limiting case of purely axial, fully developed flow, the accuracy of the method developed is satisfactory in the laminar as well as in the turbulent flow at high Reynolds numbers (Rea > 5 x 10~) within, of course, the limitations imposed by the neglect of the secondary motions. Typical examples of the results are shown in Fig. 2, discussed below. Figure 2(a) compares predicted dimensionless axial velocity contours w/~, for laminar flow through an equilateral triangular arrangement of pitch-to-diameter ratio P/D = 1.5, with the predictions of Meyder [1 I]. The agreement is very good. Predicted turbulent flow velocity contours for two values of Reynolds number are also shown in the same Fig. The dependence of the axial prcssure drop cocfllcient c~ on Reynolds number Reo is very accurately predicted: in the laminar case, the c~ Rea = const, relation has been verified by the present predictions and the values of this product are given in terms of P/D in Fig. 2(b) for both the equilateral-triangular and the square arrangements. The analytical solution of Sparrow and Loeffler [12], shown in the same figure, is in very good agreement with the present predictions. The predicted e~ for turbulent flow through an equilateral

368

K.A. ANTONOPOULOS

tot, t03 tOO

B

09/,

0.81

t06

1.10

138t02 n ~ C EQUILATERAL-TRIANGULARP/D=1.5

.L."'.,"". :", "".",L-:"-. _ -.,.

| I - - PRESENTPREDICTIONS'I ~ ~ " MEYDERD1]PREDICT~ LAMINAR \--PRESENTP 'D Cr S. R, =I.0S,aOS

..,

.... ~

ESENT PREDIO'IONS Reh~=l.2/.xl06

1 (°)

F 50

E

LAMINAR FLOW SPARROWAND LOEE0 2] SENT PREDICTIONS

TURBULENT FLOW EOUILATERAI.-TRIANGULAR P/D=t866 -t- PRESENT PREDICTIONS 10"1 = COURTAUDET AL [13] EXPERIM. S.

(J

o 30 20

10-2

10

(b)

1.0

SOJ~E

-

i

I

1.5

2.0

~ P/D

(c) 2.5

104

105

106

Reh

Fig. 2. Purelyaxial fullydevelopedflow.(a) Dimensionlessaxial velocitycontours w/~.,.(b) Product cr~ Reo as a function of pitch-to-diameter ratio. (c) Pressure drop coefficientas a function of Reynolds number.

triangular arrangement of P/D = 1.866, as a function of the hydraulic diameter-based Reynolds number Re] is successfully compared in Fig. 2(c) with the experimental data of Courtaud et ai. [13] for the same arrangement.

5.3 Purely transverse, fully developed flow 5.3.1 Equations and boundary conditions In this case (O = 90°), the terms in the general equations (!) and (2) relating to the axial co-ordinate ~ vanish and the axial velocity component w is set to zero. The problem is therefore described by the following equations: Transport equations for ¢, = u, v, k, E:

lJ a¢ C o n t i n u i t y equation:

1 d

z;l

ul)

(o .

I

(22)

Since no variations of flow conditions exist with axial direction ~, the solution of the above equations is confined to a single cross-sectional plane ~-~/. For the simulation of the fully developed character of the flow, two different procedures have been devised: the "Explicit" and the "Implicit Periodic Procedure". The former extends the computational mesh by one line EF and then transfers the profiles of the dependent variables from line DG to BA and from CJ to EF (Fig. 3a) after each iteration. The latter, which is much more economical in respect to computer time, is described below. It will be recalled that in the ADI method (Section 4) the equations are solved along grid lines, taking values at neighbouring lines as known. This is done because the coefficient

Turbulent inclined flow in rod bundles

B C

369

D E "al.l al. 2 a2.1 a2.2 a2.3 a3.2 a3.3 a3./,

(a)

z

al, 1 al, 2 a2,2 a2,3 a3,3 a3,/.

I

(c)

an.-1,n-2 an-l.n-1 an-l. n| ~n.~l an.,, J

an.1 b) MATRIX M1

H

al.n a2.n a3.n an-l. n-1 an-1. an,n

.2 a3,3

(d)

MA'mlX ~2

al'n ]

a3, n 1 " a2'n an-1. n-1 an-1.n an.n

MATRIX M3

Fig. 3. Explicit and implicit periodic procedures for purely transverse fully developed flow. (a) Grid arrangement. (b) Matrix of coefficients before first-stage eliminations. (c) Matrix of coefficients after first-stage eliminations. (d) Matrix of coefficientsafter second-stage eliminations.

matrix for each line is then tri-diagonal and can therefore, be solved directly and efficiently by a Tri-Diagonal Matrix Algorithm (TDMA). The present Implicit Periodic Procedure differs from the Explicit one in that, first, the grid is not artificially extended and, second, the cyclically repeating condition at inlet and outlet is implicitly incorporated into the ADI solution of the momentum, k and e (but not pressure) equations, at the stage when the ADI passes are made along grid lines intersecting inlet and outlet boundaries, as follows: the grid nodes along these grid lines are numbered in a cyclic fashion, as depicted in Fig. 3(a), and the finite-difference equations for each node i are written in the way suggested by equation (19) as a~O i. ~-ar.l+*l -

*

Or+l + a~r_. iOr-i + (remaining terms).

(23)

Therefore, the coefficient matrix M, of the finite-difference equations corresponding to a grid line, has the form shown in Fig. 3(b), i.e. the usual tri-diagonal structure is destroyed by the elements a~., and a,.t, resulting from the cyclically repeating condition at inlet and outlet boundaries. A simple method of solving this equation set, for which the TDMA is no longer applicable, has been devised. This method consists of a two stage elimination of unknowns in the equation set of each grid line, i.e. (a) First stage: with reference to matrix Mi (Fig. 3b), • I is eliminated by combining the first and second equations; the resulting equation is then combined with the third to eliminate ¢2, etc. The final coefficient matrix M, is of the form shown in Fig. 3(c). (b) Second stage: with reference to matrix M,, O, is eliminated by combining the first and last equations; the resulting equation is combined with the second equation to elimination 02; the resulting equation is then combined with the third to eliminate • 3, etc. The final coefficient matrix M~ is of the form shown in Fig. 3(d) and therefore can be solved by a simple back-substitution formula. The wall functions are used for imposing the boundary conditions on the rod surfaces AJI and FGH (Fig. 3a) while on the upper and lower boundaries of the solution domain BCDE and IH, respectively, which are planes of symmetry, v-velocity is zero as are normal gradients of u, k, E. 5.3.2 Results

Typical examples of the comparisons made between the present predictions and existing experimental data are displayed in Fig. 4, discussed below.

370

K.A. ANTONOPOULOS I0

1

(b) 0.1 .AMINAR FLOW in-line 1.25xl.25 • PRESENT PREDICTIONS l - - BERGEUNEl" AL [14] EXP. 0.1 i t 100 101 102

Y~

Re c

103

TURBULENT FLOW in-line 13xl.3 • PRESENT PREDICTIONS STAS.ANO SAM. [15.16]EXR 0.01 K i 104 10 5 106 Re c

1

8.

(c)

(J

1

(d)

r.) t~

¢n

c~

c; 0.1

0.1

in-line 2x2 • PRESENT PREDICTIONS - - ESDU [ 5 ] CORREL. 0.01 " - ZHU KAU~; KAS [4] COI~REL.

1°~

1°5

z°6 Rec

TURBULENT FLOW stoggered 13xl • PRESENT PREDICTIONS

msrAS.ANOSA~ [lsj EXP. 1°7

0.011

' 105

106 Rec

1°7

Fig. 4. Pressure drop coefficient as a function of Reynolds number for purely transverse fully developed flow in various arrangements.

The predicted pressure-drop coefficient c~ as a function of Reynolds number Re, for laminar flow through an in-line arrangement of transverse and longitudinal pitch-todiameter ratios Sr x S, = 1.25 x 1.25 is compared successfully with the measurements of Bergelin et al. [14] in Fig. 4(a). Similar predictions for turbulent flow through an in-line arrangement of Sr x SL = i.3 x 1.3 are compared in Fig. 4(b) with the measurements of Stasiulevicius and Samoshka [15, 16]. The agreement is satisfactory at high Re,. but at the lower values of Rec it is less so, owing to the inapplicability of the k-E model in this Reynolds-number region. The slight underestimation of the data may be due to entrance effects, since the experiments were performed in a bank with only seven rows of tubes. The predicted c~ for turbulent flow through an in-line arrangement of Sr x SL = 2 x 2 is compared successfully (especially at high Re,.) in Fig. 4(c) with empirical correlations for the same arrangement, taken from Zhukauskas [4] and ESDU [5]. Predictions for turbulent flow through a staggered arrangement of Sr x SL = i.3 x 1 are compared in Fig. 4(d) with the measurements of Stasiulevicius and Samoshka [15]. The agreement here is less satisfactory than that for the in-line arrangements displayed above. Extensive comparisons have shown that although in the laminar case, pressure-drop characteristics of both in-line and staggered arrangements are equally well predicted, in turbulent flow the predictions concerning in-line arrangements are more accurate than those for the staggered ones. The explanation lies in the wall functions employed. These, as mentioned in Section 3, are based on the assumption of low pressure gradient along the rod surface. This assumption is closer to the reality in the case of in-line than in the staggered arrangements. For this reason, in the remaining of the study only in-line arrangements will be considered. It is important to note in the above discussed Figs 4(b), (c) and (d) that the pressure drop coefficient c~ is predicted to reach a practically constant value at high Re,. This is the value which can (and has been) predicted accurately enough by the turbulence model employed. The minimum which in reality appears in the c~ vs Re,. curve of the widely spaced arrangements (but not in the closely spaced ones, as suggested by Zhukauskas's [4] collection of data) owing to the transition of the laminar boundary layer on the rod surface to a turbulent one, cannot be predicted by the present turbulence model.

Turbulent inclined flow in rod bundles

LAMINARFLOW , Rec=10 •>

1.009

371

TURBULENT FLOW,Re(;=IO5

Lo..

1.0

°_1

1.25x1.25 -(zoo16\ / / - PRESENT ((1) 1 v I---LE FEUVRE[3]

TURBULENT FLOW, Rec=lO4 L ~ 0 67

(b) / '" I

l OOj

I-- p S_ENT_ I--- SSE¥ [Z]

O,OO

(c) I TURBULENT FLOW, Re¢'-lO6

t25 xl.25

~

TURBULENT FLOW, Re¢=lO6

in- l i n e ' ~ . . . . ~ ~

2x=

\ (e)~

\

Fig. 5. Predicted flow patterns for purely transverse, fully developed flow through various in-line arrangements at various Reynolds numbers.

Typical examples of predicted flow patterns through in-line arrangements at various Reynolds numbers are displayed in Fig. 5(a)-(e). Figures 5(a) and (b) compare the present predictions for laminar (Re¢-- 10) and turbulent (Re~ = 10:) flow with the predictions of Le Feuvre [3] and Massey [2] respectively. Agreement with both authors is good. The predicted influence of Re¢ on the flow pattern is illustrated in Figs 5(a), (c) and (d), which correspond to the same rod bundle (in-line, Sr x SL = 1.25 X !.25) at Re~ = 10, 104 and 106 respectively. The predicted effect of rod spacing is shown by comparing Figs 5(d) and (e), which refer to the same Reynods number (Re¢= 106) but different spacings (Sr x SL = 1.25 x 1.25 and 2.0 x 2.0): in the widely spaced arrangement, the recirculation zone does not bridge the gap between the tubes. Predicted contours of the kinetic energy of turbulence k (normalized by the mean velocity fi ) for the arrangement of Fig. 5(b) are shown in Fig. 6(a). The maximum k occurs near the reattachment point where generation is greatest due to the steep velocity gradients which exist there. This peak value is surrounded by an area of high k which extends along the tube surface to just beyond the outflow plane. Throughout this area steep velocity gradients occur along the wall, which contribute to the high turbulence levels. Along the vertical bisector of the cross-section, k is lowest in the gap between the tubes (where the flow is constrained and therefore smooth changes occur) and increases towards the shear layer separating the vortex from the main flow. Within the latter and turbulence energy decays towards the symmetry plane where the velocity is practically uniform and generation is therefore small. The predicted effect of rod spacing on the level of k in the bundle is illustrated in Fig. 6(b), where the distribution of k has been plotted along the centre-line between the longitudinal rows of various in-line arrangements at Re,. = i0 s. It is clear that the

372

K. A. ANTONOPOULOS

in-fine

2.06xl.38,

Rec:105

/ ooo

,.~°°3 .~~ in-line, Rec:lO 5

12s,1.2s

(]

(b)

x~o

0.1 0.2 0.3 (1/, 0.5 0.6 0.7 Q8 0.9 1.0

x/sL

Fig. 6. Purely transverse, fully developed, turbulent flow. (a) Predicted contours of kinetic energy of turbulence k/~ 2. Co) Predicted distribution ofk/~ ~along the centre-line between the longitudinal rows of various in-line arrangements.

turbulence level is lower for the widely spaced bank and practically constant along the centre-line, because diffusion levelled out the non-uniformities produced by the tubes at which the turbulence is generated. By contrast, for the closely spaced arrangement the level of turbulence is higher, since more turbulence is generated per unit volume, and k varies along the centre-line, the highest values being at locations corresponding to the closest proximity to tubes and the lowest corresponding to the gap between the latter.

5.4 Fully developed inclinedflow 5.4.1 Equations and boundary conditions. In this case (0° < O < 90 °) it is assumed that the bulk velocity vector makes everywhere the same angle with the rods and the flow is fully developed. This implies that there is no variation of flow conditions with the axial direction ~, except for pressure. Hence all 0/0~ derivatives in the general equations (1) and (2) vanish with the exception of the axial pressure gradient term Op/O~contained in the source term of the axial momentum equation. Therefore, the transport equations take the form of equation (21) for • = u, v, w, k, e and the continuity relationship takes the form of equation (22). It is noteworthy that here the transverse momentum equations [i.e. equation (21) for • = u and • = v] do not contain the axial velocity component w. Therefore, in the laminar case, they are decoupled from the axial velocity field, i.e. the latter has no effect on the transverse motion. However, this is not the case with turbulent flow because of an indirect coupling introduced via P,fr: This quantity is calculated from equation (3) in terms of k and ~, the transport equations of which [i.e. equation (21) for 4, = k and • = ¢] contain the axial velocity w in the source terms S, and S,, as deduced from Table 2 and equations (6)-(12). Since the axial co-ordinate ~ does not appear in the governing equations (except for the axial pressure gradient term Op/O~), the solution of all equations may be confined to a single cross-sectional plane l-r/, with the term Op/O~ either prescribed as an input or adjusted during the solution procedure to give the desired axial mass flow rate.

Turbulent inclined flow in rod bundles

373

LAMINAR FLOW, in-line 1.25xl.25 •

>

1.00

,

(Q) TRASVERSE \ ~ / /

.

COM~,~ /_0.0,S/ R,-IO0

W/~V =1.0 1 / 1 . ~

1:o

1./, I

2.15 L-

j

.o:o

COMPONENT ~ \ ~ /

o5

1.6 .-J ~.-4-1.4

~

0:51", Rec: 100

03

Fig. 7. Laminar, fully developed inclined flow. (a) Predicted transverse flow pattern. (b) Predicted dimensionless axial velocity contours w/~ for purely axial and for inclined fully developed flow.

The boundary conditions at inlet BA and outlet DG (Fig. 3a) are imposed in the same manner as in the case of purely transverse, fully developed flow (Section 5.3.1), i.e. by use of the Explicit or the Implicit Periodic Procedure. The wall functions are employed for imposing the boundary conditions on the rod surfaces AJI and H G F (Fig. 3a), while on the upper and lower boundaries BCDE and IH respectively, which are planes of symmetry, v-velocity is zero as are the normal gradients of u, w, k and c. 5.4.2 Results (a) Flow pattern. First, an example of laminar flow calculations is displayed in Fig. 7, which refers to an in-line arrangement of Sr x SL = 1.25 x 1.25. The values of the axial Re° and the transverse Rec Reynolds numbers are 81 and 100 respectively and therefore the inclination angle of the mean velocity vector to the rod axes is 8 = arctan (6/@) = arctan (ReJRe°) -- 51 °. The transverse flow field, shown in Fig. 7(a), is independent of the axial flow component because, as mentioned earlier, the transverse momentum equations do not contain the axial velocity w. On the contrary, the axial flow field is strongly influenced by the transverse motion, since velocities u and v feature in the axial momentum equation [i.e. equation (21) for • = w]. This influence is illustrated in Fig. 7(b), where dimensionless axial velocity contours w/~v have been plotted for both inclined (O = 51 °, Rec = 100) and purely axial (O = 0 °, Re,. = 0) flow. It is evident that in the latter case the w/W contours are symmetrical about the vertical bisector of the cross-section with the maximum dimensionless velocity lying at the point furthest removed from the rod surfaces. The effect of imposing transverse flow is to destroy the symmetry of the axial motion, i.e. the isovels are "convected" towards the outlet of the cross-section as is the location of the maximum dimensionless velocity. In the case of turbulent flow, although the axial velocity component w is not explicitly contained in the transverse momentum equations, its influence is indirectly introduced into these via #el, as mentioned earlier. This influence is shown in Fig. 8(a), where transverse flow patterns are presented for Re~ fixed at 106 and for two values of inclination angle O, i.e. O = 90 ° (purely transverse flow) and O = 17.5° (Re, = 3.17 x 106). It is noteworthy, in this figure, that the recirculation bubble between the rods becomes slightly larger in size for decreasing O with consequent upward deflection of the external stream lines. The explanation is that for the inclined case (O = 17.5 °) the level of turbulence energy near C.AF

14 "~-'D

374

K. A. ANTONOPOULOS TURBULENT

FLOW.in line2.06xl.38

/

o.o_

(O) T

. _

.

~

.

J

,o=3 x'06.

•..~.~tlO5~'~,,~. . . . . . . .

\~

.

0 Re~0,Rec--106

COMINeNT \ {('0

Ib) AXIAL

.

:~L030

/7

./'/__e=o o, Rea=3xl05,Rec=0

COMPONENT~ ~ / _ _ _

e=l?ff.R~=3.2x106. Rec=106

Fig. 8. Turbulent, fully developed inclined flow. (a) Predicted transverse flow patterns for purely transverse and for inclined fully developed flow. (b) Predicted dimensionless velocity contours w/P for purely axial and for inclined fully developed flow.

the wall increases due to the extra generation associated with the axial velocity gradients in the ~ and ~ directions: this results in higher wall shear stresses thus causing slightly earlier separation. Also noteworthy is the increase of the recirculation strength from -0.0303 at O = 9 0 ° to -0.0387 at O = 17.5°. It is interesting that despite the large difference in the inclination angle O, the differences in the flow pattern outside the recirculation zone are not more than about 3%. The effect of the transverse flow on the axial velocity distribution is expected to be strong, as was the case with laminar flow. The predicted axial velocity contours w/fv for the same case of inclined flow discussed above (O = 17.5°, Rec = 106, Reo = 3.17 x 106) are shown in Fig. 8(b) together with the predicted w/f¢ contours for a purely axial flow (O = 0°). It is clear that the isovels for O = 0 ° are symmetrical about the vertical bisector of the cross-section while the O = 17.5° they are "convected" towards the outlet. The influence of the recirculation zone produced by the transverse flow (Fig. 8a) is clear in the w-velocity contours, which are biased in the direction of the recirculation. Also noteworthy is the slight change of slope of the isovels near the separation point. Another example of turbulent flow calculations is displayed in Figs 9(a) and (b), which correspond to an in-line arrangement of ST x SL ffi 1.25 x 1.25. The first of these figures shows the predicted transverse flow patterns for Rec fixed at l06 and O ---90 ° (purely transverse flow) and 8 = 77.4 °. Here, the distortion of the transverse flow pattern by the axial flow is less than that for the arrangement of Fig. 8(a), but this behaviour is expected because of the smaller difference of inclination angles in the present case. The corresponding axial velocity contours are shown in Fig. 9(b): the recirculation zone of the transverse motion exerts a stronger influence here than that shown in Fig. 8(b) because of the higher inclination angle. Thus, the isovels outside the recirculation zone are convected towards the outlet, as usual, but the internal isovels remain within the recirculation zone, i.e. the axial motion is effectively separated into two distinct zones. The inclined flow cases of Figs 8 and 9 discussed above, serve to justify the statement made earlier, i.e. that the turbulence-driven secondary motions of the axial flow have negligibly small effect in an inclined flow. This can be predicted in advance simply from

375

Turbulent inclined flow in rod bundles

TURBULENT FLOW, In line 1.25xI.25

t00

tRANsverse \ ~ ' ~ 7 ~ / - - e ~ ° . ~ ' e °. r'~'~°6 (O) CO~I;'O'I~" ' ~ / - - O=,'/.'p.Rea=2~x105. Re¢=106

(b)

[I1o.,,. o=-'.

\4 ,2 h

'.¢:'0s

s

Fig. 9. Turbulent, fully developed inclined flow. (a) Predicted transverse flow patterns for purely transverse and for inclined fully developed flow. (b) Predicted dimensionless axial velocity contours

wl~.

knowledge of the ratio fi/~. Thus, in these inclined flows [where O = arctan (~/~,) = 17.5° and O = arctan (~/~,)= 77.4 ° respectively], the mean transverse velocity ~ is 32% and 448% of the mean axial velocity ~, respectively, while the turbulence driven secondary velocities are generally not more than 2% as found by [1]. It is therefore clear that the latter are swamped by the main transverse motion. (b) Kinetic energy of turbulence. Figures 10(a) and (b) compare the predicted variation of the kinetic energy of turbulence k (normalized by the mean velocity ~2) in terms of the

in line 2.06xi.38 O=90=,Reo=0, Rec=I06

0.101 TURB,FLOW,

I

- - - e=n¢'. Re~a2x~.

0.06~ ", (10/,

r,.c=~ 6

• . . . . . . . . . . . . . . . .

- .....

0.O2 ~( C[7~: 90o 00

I

I

!

i

0.08

0.10

0.02 (10/,

0.06

I

I

I

1102

0.06

0.08

0.10

0.12

~flD

0.06 I:J

0.0~

0.02 0 0

0.0~

0.12

Y/D

Fig. 10. Predicted variation of kinetic energy of turbulence k/O' as a function of the dimensionless distance y/D from the wall for a purely transverse and for an inclined fully developed flow. (a) At angular position 0 = 90°- (b) At angular position 0 = 60'%

376

K.A. ANTONOPOULOS 5.0

LAM. FLOW, in line 1.25xl.25 ~ PRESENT PREDICTIONS [12]FOR ~c-O

(~ u

101

TURB.FLOW. in line 1.25xl.25

1.0 UI O.5

II

ioo-//

0.1L..d.

O.I! 10

50

(a)

0.1

Rec=106

",.,,\\\\

TI,.RaFLOW, in line

100

500

Recz

o'

,

, 105

¢b)

1"01 ~ ' " ~ TURS.FLOW,inIine 1.25xl.25 _ <..L " ~ . ~ " [ wI~S. PRI~ ,105<~R_e<~105

1.25xl.25

"°/

ok (c)

°u(']'

061

\-....

°'t

0.01

0.00t

,

,

1o5

,

i

lO5

105 , I 10 5 Reh 10 7

,

:

I°7

90 80 (d )

\'"

?0

60

50

40

30 20 10 e (degrees ]

0

Fig. 11. Fully developed inclined flow. (a) Predicted axial pressure-drop coefficient as a function of the axial Reynolds number, for various values of the transverse Reynolds number. (b) Predicted transverse pressure-drop coefficient as a function of the axial Reynolds number, for various values of the transverse Reynolds number. (c) Predicted axial pressure-drop coefficient as a function of the axial Reynolds number, for various values of the transverse Reynolds number. (d) Ratio .4 of total to transverse (at Re~ = Re) pressure-drop coefficients as a function of the inclination angle.

dimensionless distance y/D from the wall at two angular positions for a purely transverse ( 8 = 9 0 °, Rccffi 106) and for an inclined flow (O = 17.5 °, Re,.= 106, Re, = 3.17 x 106) through the in-line arrangement of Sr x SL = 2.06 X !.38 discussed earlier. It is clear that at a fixed Rec the level of k is higher in the inclined than in the purely transverse flow due to extra turbulence generation associated with the axial velocity gradients in the ~ and r/ directions. The augmentation of k is greatest near the wall, as expected. It is also worth noting that for both cases (i.e. the purely transverse ahd the inclined case) the highest turbulence energy is observed at the location corresponding to the top of the tubes (~ --90 °) because steeper velocity gradients occur there due to the narrow passage, but also because this location is closer to the reattachmcnt point where k attains high values, as already shown in Fig. 6(a). (c) Pressure-drop coe~icient. In the laminar flow, the transverse pressure-drop coefficient c~ is a function of the transverse Reynolds number Re,. and not of the axial one Re,, since the transverse flow field is independent of the axial motion, as mentioned earlier. The predicted dependence of cw on Rec has already been shown in Fig. 4(a). On the contrary, the axial pressure-drop coefficient c~ depends on both Re,. and Re,, as shown in Fig. I I (a), which corresponds to an in-line arrangement of Sr x SL -- 1.25 x 1.25. It is seen in this figure that for constant value of Re,., c~ varies inversely with Re, as in purely axial fully developed flow. The effect of transverse flow is to increase uniformly Cp~ by an amount which is a function of Re~. The minimum cp~ is observed at Re, = 0 (at which the present predictions agree very well with the analytical solution of Sparrow and Loeffler [I 2]), thus

Turbulent inclined flow in rod bundles

377

showing that the axial component of the flow :is distributed in purely axial flow in such a way as to provide the least resistance. Turbulent flow predictions for the same arrangement discussed above are displayed in Figs l l(b)-(d): The predicted effect of the transverse Re~ and the (hydraulic diameterbased) axial Re~ Reynolds numbers on the transverse pressure drop coefficient c~ is small in the Reynolds number ranges examined, as shown in Fig. 1 l(b). The predicted axial pressure-drop coefficient c~ as a function of Re~ for various values of Re~ is shown in Fig. l l(c). It is evident that the minimum c~ is obtained at Re,. = 0. This behaviour suggests that, as was the case with the laminar flow, the axial velocity is distributed at Re~ = 0 in such a way as to provide the least resistance. The distortion produced by the transverse flow causes the axial resistance to increase. The predicted influence of inclination angle O on the pressure drop coefficient cp in the direction of the inclined flow is compared below with the recommendations of ESDU [5]. These recommendations, based mainly on theoretical considerations, suggest that cp, at a value Re of the inclined flow Reynolds number, can be calculated as cv = A c w ,

(24)

where cp~ is the pressure-drop coefficient in a purely transverse flow at Re~ = Re and ,4 postulated to be a function of O only (for a specified geometry) for all Reynolds numbers Re > 103. The recommended form of this function is given in Fig. 1 l(d), which corresponds to the same arrangement discussed above. In the same figure, the present predictions for 105 < Re < 106 have also been plotted. It is seen that cp increases with O, the maximum cv being obtained at O = 90 ° and the minimum at O = 0 °. The agreement of the present predictions with the ESDU [5] recommendations is good for the higher values of O, while increasing disagreement is observed with decreasing O. It is interesting to note that the present predictions agree with the ESDU suggestions in that the predicted A vs O curve is independent of Re. 6. C O N C L U S I O N

The prediction method developed allows calculations of flow in rod bundles to be made for a wider range of circumstances than has hitherto been possible. The numerical technique has been tested successfully in the laminar case, where no turbulence model uncertainties were present. In the turbulent case, despite the limitations imposed by the turbulence model, the method gave encouragingly accurate results at high Reynolds numbers, as was proved by comparisons with existing data in the limiting cases of purely axial and purely transverse flow. Although comparisons in the case of fully developed inclined flow were not possible owing to absence of related information, it is believed that the present predictions give a valuable insight into this almost unexplored flow case. Among the main findings is that although the effect of the axial flow component on the transverse motion is very small, the latter strongly influences the axial velocity distribution and the axial pressure-drop coefficient. Thus, at a fixed transverse Reynolds number Rec, the axial pressure drop coefficient c,. decreases with increasing axial Reynolds number Re.h in about the same way as in a purely axial flow. The minimum c,. is obtained at Rec -- 0 thus suggesting that the axial velocity is distributed in a purely axial flow in such a way as to provide the least resistance. At constant R ~ , the effect of transverse flow is to increase c,. by an amount which is mainly a function of Re~ and to a lesser extent to Reh.. The predicted effect of inclination angle O on the pressure-drop coefficient cp in the direction of the inclined flow is that the latter increases with O at a constant Re, the maximum cp being obtained in a purely transverse flow (O = 90 °) and the minimum in a purely axial flow (O -- 0°). Also, at a constant O there is practically no variation of cp with Re at the higher Reynolds numbers. REFERENCES I. P. Carajlescov and N. E. Todreas, Experimental and analytical study of axial turbulent flows in an interior subchannel of a bare rod bundle. ASME Winter Meeting, 75-WA/HT-SI (1975).

378

K . A . ANTONOPOULOS

2. T. H. Massey, The prediction of flow and heat transfer in banks of tubes in cross-flow. Ph.D. Thesis, Central Electricity Research Laboratories, Leatherhead, Surrey (19767. 3. R. F. Le Feuvre, Laminar and turbulent forced convection processes through in-line tube banks. Ph.D. thesis, Imperial College of Science and Technology 0973). 4. A. Zhukauskas, Heat transfer from tubes in cross-flow. Adr. Heat Trans. 8, 93-160 (19727. 5. ESDU, Pressure loss during cross-flow of fluids with heat transfer over plain tube banks without baffIes. Engineering Sciences Data Item No 74040 (19747. 6. F. P. Kazakevich, Influence of the angle of approach of a gas stream on the aerodynamic resistance of a tube bundle. Izv. Faes. Teglotekh. Inst. Imeni F. E. Dzerzhinskogo, No. 8 (1952). 7. B. E. Launder and D. B. Spalding, Mathematical Models of Turbulence. Academic Press, New York 0972). 8. K. A. Antonopoulos, Prediction of flow and heat transfer in rod bundles. Ph.D. thesis, Imperial College of Science and Technology, London (19797. 9. H. Schlichting, Boundary Layer Theory. McGraw-Hill, New York 0968). 10. L. S. Caretto, A. D. Gosman, S. V. Patankar and D. B. Spalding, Two calculation procedures for steady, three-dimensional flows with recirculation. In Proc. 3rd Int. Conf. Numerical Methods in Fluid Mechanics, Vol. 2, pp. 60--68. Springer Verlag, Berlin 0972). I I. R. Meyder, Solving the conservation equations in fuel rod bundles exposed to parallel flow by means of curvilinear--orthogonal co-ordinates. J. comput. Phys. 17, 53-67 0975). 12. E. M. Sparrow and A. L. Loeffler, Longitudinal laminar flow between cylinders arranged in regular array. A.LCh.E.J. 5, 325-330. (1959), 13. M. Courtaud, R. Ricque and M. Martinet, Etude de laertes de charge dans des conduites circulaire contenant un faisceau de barreaux. Chem. Engng. Sci. 21, 881-893 0966). 14. O. P. Bergelin, G. A. Brown and S. C. Doberstein, Heat transfer and fluid friction during flow across banks of tubes, part IV. Tram A S M E 74, 953-960 09527, 15. J. K. Stasiulevicius and P. S. Samoshka, Lienmos TSR MA Darbai, ser. B4(35), 77, 83 0963). 16. J. K. Stasiulevicius and P. S. Samoshka, Lietuvos TSR MA Darbai, ser B4(55), 133 0968).