Turbulent flows in complex rod bundle geometries numerically predicted by the use of FEM and a basic turbulence model

Nuclear Engineering and Design 99 (1987) 351-363 North-Holland, Amsterdam

351

TURBULENT FLOWS IN COMPLEX ROD BUNDLE GEOMETRIES NUMERICALLY P R E D I C T E D BY T H E U S E O F F E M A N D A B A S I C T U R B U L E N C E M O D E L

H.G. K A I S E R a n d W. Z E G G E L Institut fftr Raumflugtechnik und Reaktortechnik, Technische Universitiit Braunschweig, PockelsstraBe 14, 3300 Braunschweig, Fed. Rep. Germany

1. Introduction Newer projects in nuclear reactor design tend to higher conversion ratios. An up-to-date PWR has a conversion ratio of approximately 0.53. Whereas APWRs are planned to have up to 0.95 and the breeder reactor is supposed to have a ratio better than 1.0. High conversion ratios necessitate tightly packed fuel rod lattices. Together with high burnup, necessary for economic efficiency, the slender fuel rods show a tendency to bend. The result of bent fuel rods is a distorted lattice. A further item which leads to irregular lattices are the tolerances in fuel rod assembly. The absolute values of tolerance which can be seen as fixed become relatively more important in tightly packed lattices. One aspect of the problem of distorted lattices is how the velocity distribution as well as the local heat transfer is affected. There are many possible cases for irregular cross sections of fuel elements. Not all of them can be measured. So there is a need for a computational tool which can manage the problem in a satisfactory way. For heat transfer calculations the velocity distribution has to be determined before or simultaneously with the temperature distribution, because the heat transfer depends strongly on the velocity distribution. If one considers the fluid properties independent of temperature, which is allowable in most technical cases, the velocity distribution may be calculated separately from the temperature distribution. For regular subchannels like the symmetric part of a triangular channel ((~) in fig. 1) several calculation procedures exist which give good approximations for momentum and heat transfer. Most of these procedures use very sophisticated turbulence models which always impose an iteration technique and some even demand the solution of additional differe.ntial equations other than the equations for momentum and heat transfer. Therefore most of these procedures are inefficient under the aspect of comput-

ing time and storage requirements to be used for larger solution domains, necessary to handle when no symmetry exists, i.e. in design tasks. There are only two computer codes known to the author which have been used to solve bigger solution domains than one subchannel. These codes are the VELASCO code developed by Eifler/Nijsing [1] and a F E M code from Slagter [2,3]. The idea of the VELASCO code is to assume a velocity profile perpendicular to the wall which depends mainly on the local wall shear stress. The momentum fluxes for a defined subzone under further assumption of secondary flow parallel to the wall and a global factor for anisotropic turbulence are summed. A complex geometry is put together from several subzones. An iteration technique makes the subzones fit together and determines the wall shear stress distribution. This code has one big advantage: It is extremely fast even for large solution domains. The main disadvantage is the inflexibility to turbulence model variations, because it does not solve the Reynolds equations directly.

c~nnetwoLL /

0 0 2 9 - 5 4 9 3 / 8 7 / $ 0 3 . 5 0 © Elsevier Science Publishers B.V. ( N o r t h - H o l l a n d Physics Publishing Division)

/ ~

/ ~

~ k

- -- - subzone . . . . divisionfineSDL . . . . . ymmetryline

iii!: ;!:rZW Fig. 1. Cross section of a w~l bounded seven rod bundle.

4

I

I

I

i

0

0

>

m

3

.

......:,:~,:!.!:~.~::~::~i::::ii:ill

.J

x

3

J ~

-:±

Ii x

-6

m

z/

..........

] i

_L

,

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

#]-
/

~ .

z~

H.G. Kaiser, W. Zeggel / Turbulent flows in complex rod bundle geometries

353

for momentum transfer, the same solution procedure and sometimes even the same domain discretisation can be used to solve it. So far Slagter has used two turbulence models namely the algebraic mixing length model after Meyder [4] and the one-equation model after Wolfshtein [5]. Both models do not account for secondary flow but consider anisotropic turbulent momentum transfer. The mixing length model as well as the one-equation model lead to a solution procedure demanding an iteration which is an essential disadvantage under the aspect of solving complex flow domains or heat transfer problems. The one equation model additionally suffers from the necessity to solve one extra differential equation. Therefore in this paper a new algebraic turbulence model is shown which is based on anisotropic eddy viscosities. This model does neither demand an iteration nor the solution of additional differential equations.

2. Geometry Fig. 2b. Adopted model for eddy viscosity ep parallel to the wall in an EFC.

More sophisticated turbulence models like one or two equation models can not be implemented. The F E M code developed by Slagter is essentially a solution procedure for a differential equation like the momentum transport equation for fully developed flow. Consequently, it allows the use of any turbulence model that fits the Reynolds equations (provided the potentially necessary iteration converges). Since the type of differential equation for heat transfer is similar to that

12 10 08 o_

w06 O4 02~

CL ~0

O0 00

/

/

f

J

0.5

Y/~

10

Fig. 2c. Radial distribution of normalized eddy viscosity Ep//~p.

Geometric parameters play an important role in calculation of velocity distribution in rod bundles. Therefore definitions for occurring parameters and subchannels are given in the following section. All types of subchannels in hexagonal rod bundles with surrounding walls can be found in a seven rod arrangement (fig. 1). Three different types can be defined: the elementary flow cell (EFC) (no. 1), the wall channel (no. 2), the corner channel (no. 3). All subchannels are again subdivided into so called subzones by subzone division lines (SDL). These lines define either points of equal perpendicular distance to neighbouring walls or lines showing shortest or longest perpendicular distance of faceing walls. The SDL are good estimates for ridge or fall fines in the contour, which occur if the local velocity is plotted vertically over the solution domain (see fig. 10). In symmetrical subchannels of an infinite array the SDL are identical with the ridge and fall lines of the velocity contour and also identical with the geometric symmetry lines. Important geometric parameters in rod bundles are the pitch P, which is the distance between the center of two rods, the diameter of a rod D and the vertical distance W from a rod center to the nearest channel wall. Within one subzone, Y denotes the vertical distance of a certain point to the nearest wall and )3 the vertical distance from the wall to the next SDL (see fig. 2).

354

ft.G. Kaiser, 14/.Zeggel / Turbulent flows in coml;l~.~rod bundle geometrie~

3. Governing differential equations

4. Present algebraic turbulence model

The mean motion of turbulent flow can be described by the Reynolds equations for momentum transport and the equation for continuity. Considering fully developed longitudinal flow in rod bundles, the terms supposed to vary in the longitudinal direction x can be set to zero with the exception of the pressure gradient which is constant.

The development of a new algebraic turbulence model was intended to avoid iteration and the solution of additional differential equations to the momentum equation. The idea of anisotropic eddy viscosities has been adopted as a basis for the algebraic model with the main directions parallel and perpendicular to the nearest wall. In the major part of the solution domain the eddy viscosity es perpendicular to the wall is determined by the Prandtl/Kolmogorov expression [6].

0(~) 3x = 0 ,

3p ~ x =c°nst"

In the notation used, x is the coordinate parallel to the rod axis, p the time averaged pressure and u, v and w the time averaged velocities in the directions of x, y and z respectively. Measurements of v and w showed that the magnitude of secondary flow velocities is 1% or less of the bulk velocity. Their effect on the momentum transfer is neglected in the present turbulence model so that the goveming equations can be further reduced by setting v, w to zero. Employing Boussinesq's hypothesis of eddy viscosity the resulting equation for momentum transfer reads:

3

In contrast to one or two equation models, the turbulent kinetic energy k = (u '2 + v'2 + w~2) (with u', v' and w' being local velocity fluctuations) and the mixing length l are not defined by differential equations but by algebraic expressions which depend on the geometry. So far the validity of these empirical expressions is limited to bundle geometries with 1.02 < ( P / D or W / D ) < !3. The algebraic formulation for the kinetic energy k seems allowable since most measured or calculated distributions for k show similar tendencies for different subchannel types. To illustrate the adopted form of distribution, fig. 2a shows a qualitative sketch of k ~ for an EFC in a triangular lattice, k + is the normalized turbulent kinetic energy k += k / u * 2 .

where # is the dynamic viscosity and p the fluid density. The eddy viscosity ~ is a fourth order symmetric tensor. Generally, eddy viscosities are defined perpendicular e~ and parallel ep to the nearest wall. In this coordinate system only the spur of the fourth order tensor differs from zero. If another coordinate system is chosen for calculation, the four occurring eddy viscosities in eq. (1) can be calculated by the following expressions:

where u* denotes the mean friction velocity of a subchannel. Perpendicular to the wall the distribution of k ~ can be approximated by a cubic equation which is defined by four boundary conditions: k~- defined by expr. (4a,b), (5a,c) and an interpolating cosine at (Y+ = Y . u * / ~ = 50), k ~ defined by expr. (4c,d), (5b,c) and an interpolating cosine at ( Y = 33) (k~,, + k ~ ) / 2

+ (k~-

k?s,)'0.15

e. . . . .= G cos24' + ep sin24',

(2a)

at the middle between the two points defined above,

COS24',

(2b)

3k-/3Y=O

(2c)

The values for the k ~ and k ~ vary with P / D or W / D , respectively, and the coordinate parallel to the rod or channel wall. Parallel to the wall the distribution of k shows different tendencies close to the wall and along the symmetry line. Close to the wall (fine AD) the maxi-

e. . . . = es sin24' + ep

e....... = e ..... = (e s - %) cos 4' sin 4',

where 4' denotes the angle of rotation from the coordinate system of definition to the system of calculation.

a t ( Y = ~).

355

H.G. Kaiser, W. Zeggel/ Turbulentflowsin complexrodbundlegeometries for 1.0 < (P/D or W/D) < 1.5;

mum of k + coincides with the maximum of the local velocity gradient perpendicular to the wall and has its minimum at the other end of the subzone. Along the symmetry line (BC) k + has a maximum in the middle of BC, where the velocity gradient parallel to that line has its maximum. The maxima at the wall kwmax+ and at the symmetry line k~4max a r e given by empirical algebraic expressions (4a,b,c,d) which have been deduced from measured data [9,17-22] (see fig. 3).

+

kwm~-

P/D

11.73. W/D + 17.53

+

kM . . . . k ~,lmin

for(P/D or W/D) > 1.26; + P/D kMmax = --17.05" W/D + 21.26

(4c)

k Mmax

+ k wmin

k Mmin

+

I= Y3~.3/k~v (4d) 1=)310.25

for( P/D or W/D) > 1.2; k ~Vmax k ~Vmin

+

k wma~

- 1.0

(5c)

The ratios of corresponding maxima and minima are also given by algebraic expressions (5a,b,c). Between the maximum and minimum of a line (AD or BC) a cosine is taken as an interpolation function. To calculate the eddy viscosity es perpendicular to the wall from expression (3) the mixing length I still needs to be defined. The following distribution was adopted which is a modified form of that used by Slagter [3] for a one-equation model.

for 1.0 < ( P / D or W/D) _< = 1.2; + kMmax = 0.8

(5b)

for(P/D or W/D) > 1.5.

for 1.0 < (P/D or W/D) < = 1.26; (4b)

P/D + 1.9 W/D

for 1.0 < (P/D or W/D) < 1.5; +

(4a)

+ k Wmax = 2.75

0.6"

+

for Y_< 0.25)3,

. ¢r 0.15sm(~(~-

0.8" P/D

(6a)

Y_

0.25))]

(6b)

for Y > 0,25)3.

W/D + 2.2

(5a) Most experimental data do not range closer to the wall x /8, 13/

t t

o/17/ [] /18/

k* t~ " 1 ~ Mmax I ~x ,

• /19/ + /20/

;t \0.8 0

k÷

t

~.o

t 6.0 k~Cmax

i.~

pig w/d~

p/d 11.73. w/d ÷ 17.53

2.0

•

Q

'•_

i.6

i.z

1.o

i.~.

p/d w/d~

i.6

x

1.t~

[]

Mmax

•

t 1.4

•

k ~ mox

÷

k Hmin

'1-

1.2

- -

.

-

1.0

10

1.2

.

kwmtn 1.2

1.t~

~16

pld ---. wld

p/d . 2.2

~

tc 1.0

•

121

1.2

.

•

1.t,

"~

pld --... w/d

Fig. 3. Turbulent kinetic energy measured data and deduced expressions.

i.6

356

H.G. Kaiser, W. Zeggel / Turbulent flows m comple v rod bundle geometr e~

than Y+ = 50. Therefore, the distribution of k + in the region Y+ < 50 cannot be described satisfactorily. This drawback was overcome by expressing the eddy viscosity with a formula developed for a rectangular channel with the assumption that two-dimensional effects so close to the wall can be neglected. The adopted formula developed by Reichardt [7] reads: e~ = vx [ Y+ - 11- t a n h ( Y + / l l ) ] .

cipal has been used. The idea of the variational method in the FEM is to find a functional 1 which has the equation (2) as its Euler-Lagrange equation.

I 0p 0u ~u +2(v + e,.~,:l~Ty ~7 + u 3-7x [! d G.

(7)

A v o n Karaman constant of g = 0.38 was found to give good results. The eddy viscosity ep parallel to the wall was formulated directly without using the Prandtl/Kolmogorov expression. If the eddy viscosities ep calculated by Rehrne [8-11] from measured data are correct, a very irregular distribution for the mixing length would have to be taken to relate Ep to the turbulent kinetic energy k. The modeled distribution for ep is shown in a three dimensional sketch for an EFC in fig. 2b. In the narrow gap between two rods ep has high values, while at the other end of the subzone ep is low. In radial direction ep has a not quite symmetric hat-like profile (fig 2c) to approximate the natural distribution. Neelen [14] deduced an expression for ep/(33, u*) as a function of 33 from measured data. He correlated the radially averaged values of ep measured by Rehme [8,9] ep = exp[18.53 e 9.,~/R+o.m8) _ 0.833/R - 1.3]. u'33

(8) An overly high momentum transfer in peripheral direction resulting from this correlation was compensated by taking the values from eq. (8) as maxima (confer fig. 2b). An important aspect of the developed turbulence model is the fact, that in all cases the friction velocity u*, which is used in several expressions, is the mean value of the subchannel considered. It is estimated by the help of the friction factor for a circular tube with the same hydraulic diameter as the subchannet considered.

5. Solution procedure The finite element method (FEM) was used as numerical solution procedure because of its flexibility to geometric variation of the solution domain. To derive the necessary element equations representing the governing differential equation (2) a variational prin-

(9)

This functional must be rendered stationary by finding a distribution for u ( y , z) which makes the first variation of 1 vanish. In the FEM the solution is elementwise approximated by linear, parabolic or higher order expressions which depend on discrete values of the unknown variable at the knots of each element. For this investigation, a six knot triangle with parabolic form function was chosen because of its geometric flexibility and good accuracy expressing strongly bent solutions. The solution within one element is described by u(y, z)=(N,).(u,),

i=1,6,

(10)

where N, are the form functions and ui the unknown velocities at the knots i. If expression (10) is used in (9), the functional has a stationary value when the first derivative, with respect to the unknown velocities ui, vanishes. This leads to six equations which must be formulated in the element routine of the FEM code. By superposition of all element equations of the solution domain a set of simultaneous equations is received. The solution of this set of equations is the approximation for the solution of eq. (2). The FEM calculations have been carried out with the help of the FEMAP kit system developed by Harbord [15]. 5.1. Domain discretization

The solution of a FEM calculation should not depend on the domain discretization (grid). This means the grid must be sufficiently fine. The necessary discretization for a wall channel hs been investigated under the aspect of computer time and storage economy. It was determined that in the radial direction three zones are to be defined in which a certain number of elements are sufficient to produce results that do not significantly change when the grid is refined. The zones and the necessary number of elements are: I: II: III:

0 _<_Y~ _< 5 0 < Y+ <_ 50 50 < Y ~ < ( Y = 33)

2 elements radially, 6 elements radially, 4 / 7 elements radially.

H.G. Kaiser, W. Zeggel / Turbulentflows in complex rod bundle geometries

357

Zone

I+1I III

Fig. 4. (a) Demonstration grid for a wall channel, (b) realistic grid for a wall channel.

Since the radial width of the third zone varies with the azimuth, the radial discretization is also varied (see fig. 4). In the peripheral direction the 90 degree sector of the wall channel along the rod is divided into 100 elements in the first zone, 52 elements in the second zone and 28 elements in the third zone. The discretization parallel to the straight wall shows a similar structure. The number of elements there is always half of that for the 90 degree rod section. Fig. 4a shows a grid where all zones can be distinguished. This grid is not realistic and only used for demonstration. A correct grid for a Reynolds number of 123000 is given in fig. 4b. It consists of 2061 elements and 4338 knots. The radial zones I and II resembling the viscous sublayer and the transition zones have shrunk to a thick line.

6. Results To prove the applicability of the developed solution procedure to rod bundle geometries several subchannel types for which measured data were available have been numerically investigated. The intention was to show that the three types of subchannels occurring in rod bundle geometries can be adequately simulated. Asymmetric flow situations resulting from distorted geometry or combination of wall and central channels have been

studied by comparing calculated solutions with measured data from Rehme [11,12].

6.1. Elementary flow cell (EFC) The experiments of Bartzis/Todreas [16] and Carajilescov/Todreas [17] have been chosen as reference for an example of an EFC. Fig. 5a shows isovel

~

122

12 1 20

\\\\\ \

\\\\\\

Um

"~~105 0.80.595180

08 0.9

I0

Fig. 5. (a, left) Isovel plot of measured data [17] for an EFC, P/D = 1.123, Re = 27000, (b, right) present prediction.

H.G. Kaiser, W. Zeggel / Turbulent flows in complex rod bundle geomet,ve~

358

1()5

}

EIO0 i

Table 1 Measured (f,,e~) and predicted (fp,c) friction factors for the subchannels considered

o- fd--

Z

]

o exp dataBartz,s

0%,t

!

--

pres

!

i

0

SO

100

15 0

Comer subchannel

1.148 1.148 123000 0.01751 0.01751

1.071 1.071 59700 0.02606 0.02031

prediction_~

090 085

Wall subchannel

200 250 300 waU coordinate ¢' (degree)

Fig. 5c. Wall shear stress distribution of an EFC measured by Bartzis [16], compared to present prediction.

plots of measured and fig. 5b calculated distributions of U/Um with Um being the mean velocity. The calculated and measured friction factors are depicted in table 1. The distribution of the wall shear stress ~w/Tw,m is given in fig. 5c. The agreement of calculated and measured data is satisfactory. The deviation between predicted and measured values can be minimized by the use of a more complicated version of the eddy viscosity model [23]. 6.2. Wall channel

Fig. 6a shows isovel plots of a wall channel ( P / D = W / D = 1.148) measured by Rehme [9]. The calculated

P/D W/D Re f,,¢.~ fpr~

.10, Omax /Umean

Re

. . . .

--t23 000

o7

~111t k._

I~'

97

large

] 071 ! .026 54600 0.01921 0.01755

1.071 1.118 105000 0.01766 0.0174

'~\'\ h "," "

%\\\ '

small

distribution in fig. 6b agrees well with the measured one. The distribution of the wall shear stress related to the mean value of the entire subchannel is shown in fig. 6c. The measured and calculated data agree qualitatively and quantitatively well with the exception of a zone from 0 to 30 degrees at the rod wall and close to the rod to duct gap. The obvious asymmetric experimental results are a consequence of a non fully developed flow situation as Rehme recently reported [24]. Comparing the computed distributions of VELASCO [9] and Slagter [3] to the present predictions demonstrates that an adequate algebraic turbulence model based on anisotropy of eddy viscosities ignoring secondary flow can approximate rod bundle flow as well as more complex models. The calculated friction factor for this channel can be found in table 1.

I /~\~,,\\",,\ \~X_- 1121 122 "X~" \\"'

I

Asymmetric subchannel

--: =

1.23

1 20

%~1 ~

12088s o.,o

0.90

Fig. 6. (a, left) Isovel plot of measured data [9] for a wall channel, P / D = 1.148, Re = 123 000, (b, fight) present prediction.

H.G. Kaiser, W. Zeggel / Turbulent flows in complex rod bundle geometries 6.3. Corner channel

1.2

Since no measured data of sharp edged comer channels with an angle of 120 degrees (hexagonal rod bundle) was available to the author, the 90 degree corner ( P / D = W / D = 1.071) measured by Rehme [13], fig. 7, was used to compare measured and calculated data. Fig. 7 again shows isovel plots of calculation and measurement as well as distributions of the wall shear stress for the channel wall. The data predicted by VELASCO

1.1 E 3~

1.0 0.9

0.a~-

359

+ i 20

40

i

i

I

I

I

I_.~

60 80 100 120 WoLf Coordinate x/ram; ~o/deg

Fig. 6c. Comparison of measured and predicted data of wall shear stress.

Fig.(Tb)

1.2

E

1.0 + °'~//A

L__ 0,8

)

exp. x-direck exp.y-direct. ----

0.6

VELASCO pres. predict.

L

I

i

2'0

40

Fig.(Tc)

/ " \

i

. _ ~i I L L' 80 [ram]

60 x,y

0.93

0.9 0.8

1.1 '1 Fig. 7. (a) Isovel plot of present prediction for a 90 degree comer channel P / D = W / D = 1.071, Re = 59700, (b) isovel plot of measured velocity distribution, (c) wall shear stress variation along the channel wall; measured and predicted data.

360

ft.G. Kaiser, W. Zeggel / Turbulent flows in conq~h'r rod bundle .~eome',~t~',

code [13] do not range to the very edge of the corner since it is not possible to calculate sharp edges with this code. The overall agreement is good including the friction factor in table 1. 6.4. Asymmetric flow situation If one connects two channels with different hydraulic diameters the flow in the wider channel usually drags along that in the smaller one. This is the case at the intersection of an E F C and a wall channel or central channels of radially distorted lattices. To investigate the effect of m o m e n t u m transport via subchannel boundaries experiments have been done by Rehme [11,12]. Here, the four rods in a rectangular channel are moved out of the channel's centerline. The

isovel plots of measured fig, 8a and calculated fig. 8b data are qualitatively the same. Quantitatively however, they differ more than the symmetric channels discussed above. The same is true for the distribution of the wai! shear stress. The results from the VEI,ASCO code show larger discrepancies in the smaller subchanncl. The friction factor for the two subchannels can be taken from table 1, 6.5. Set;en rod bundle All regular cases discussed so far can be found in a symmetric seven rod bundle. For easier calculation the symmetry can be used so that only a 30 degree sector has to be discretisized. In fig. 9a the isovel plot of the normalized velocity u/urn is shown. It is easy to see

IIIIIIIIIIIIIII Fig. 8. (a) Isovel plot of the measured velocity distribution [11,12] in two interconnected subchannels with different hydraulic diameters P / D = 1.071; (W/D)1 = 1.026, (W/D)2 = 1.118; Re = 79 800, (b) prediction of velocity distribution of the channel shown in fig. 8a, (c) predicted and measured wall shear stress distribution for the 180 degree rod sector.

361

H.G. Kaiser, W. Zeggel / Turbulent flows in complex rod bundle geometries

0.95 0.9

0.85 0.9

O.B

Fig. 9a. Predicted isovel plot u~ u m for a symmetric 30 degree sector of a wall-bounded 7-rod bundle.

vertical parameter is shown in fig. 10 for the whole rod bundle. This graphic representation helps for easy understanding of the situation and clearly shows tendencies which are more difficult to see in other types of presentation.

that the large wall channel drags along the smaller corner and central channels. The wall shear stress distribution for an outer rod is given in fig. 9b. A three dimensional plot with the normalized velocity u / u m as

1.5

-

EFC

wait subchonnet

=

corner subchonnet

1.3

.J

1.1

E

09

J

....._

0.7

0.5 180

150

120 =

90 60 watt coordinate ~P

30

Fig. 9b. Predicted wall shear stress "g,,/'rw,' at the perimeter of the outer rod.

I

362

tt.G. Kaiser, 1+: Zeggel / Turbulentflows in comph-r rod bundle geomeu'ic~

Fig. 10. Three dimensional plot of the velocity a/um distribution for a wall bounded 7-rod bundle (P/D=I.1, W/D =1.1, Re = 100000).

7. Conclusions

Acknowledgement

From the results shown in this paper, one can see that fully developed turbulent fluid flow, if limited to rod bundle geometries, can be simulated using an algebraic turbulence model which provides the possibility of straight forward calculation [25]. The model is based on the idea of anisotropic eddy viscosities neglecting the effect of secondary flow. In contrast to an one equation model like that used by Slagter [3], no iteration is necessary which saves about 75% computing time if one assumes three iteration cycles to be necessary for a proper solution, Additionally, the solution of a second differential equation, for example the transport equation for the kinetic turbulence energy, can be avoided by the use of the empirical correlations, The developed simulations procedure for turbulent fluid flow is a good basis for further investigations of heated rod bundles. In addition to the use of this empirical turbulence model it is possible to apply more complex turbulence models like one or two equation models with the developed FEM code for theoretical investigations on turbulence structure.

The sponsorship of the Deutsche Forschungsgemeinschaft has to be acknowledged gratefully.

Nomenclature dh fp~. f, .... k k k ~v k~ p u, v, w u',v',w' •lli x,y,z D G I N P W

equivalent hydraulic diameter, predicted friction factor, measured friction factor, turbulent kinetic energy, normalized kinetic energy, = k / ( u * ) 2 k + close to the wall (Y~ = 50), k + at the subzone division line SDL (Y = ~), time averaged pressure, time averaged velocities in direction of X,T~Z, velocity fluctuations in direction of x , y , z , mean value of the time averaged velocity for one subchannel, Cartesian coordinates, rod diameter, domain, functional, form function, rod pitch, distance rod-wall (fig. 1),

H.G. Kaiser, W. Zeggel / Turbulentflows in complex rod bundle geometries Y

distance of a certain point perpendicular to the nearest wall (fig. 2a), eddy viscosity, e perpendicular to the wall, e parallel to the wall, radially averaged eddy viscosity % / 0 3 . u*) (eq. (8)), element of the fourth order tensor, of eddy viscosity, of eddy viscosity, of eddy viscosity, von K a r m a n constant, kinematic viscosity, fluid density, wall shear stress, mean wall shear stress of a subchannel, angle.

8 Es Ep

~p 8xyxy xzxz /?xyxz /?xzxy /¢ P

p % 'rw, m

Abbreviations EFC SDL

elementary flow cell (central channel), subzone division line.

References [1]

[2]

[3]

[4]

[5]

W. Eifler and R. Nijsing, VELASCO - Velocity field in asymmetric rod configuration, Report EUR-4950e (1973). W. Slagter, Finite element analysis for turbulent flows of incompressible fluids in fuel rod bundles, Nucl. Sci. Engrg. 66 (1978) 84-92. W. Slagter, Finite element solution of axial turbulent flow in a bare rod bundle using a one-equation turbulence model, Nucl. Sci. Engrg. 82 (1982) 243-259. R. Meyder, Bestimmung des turbulenten Geschwindigkeits- und Temperaturfeldes in Stabbilndeln mit Hilfe von krummlinig orthogonalen Koordinaten, KfK (Kernforschungszentrum Karlsruhe) 2029 (1974). M. Wolfshtein, The velocity and temperature distribution in one-dimensional flow with turbulence augmenta-

363

tion and pressure gradient, Int. J. Heat Mass Transfer 12 (1969). [6] B.E. Launder and D.B. Spalding, Mathematical Models of Turbulence, (Academic Press, London, New York, 1972). [7] H. Reichardt, Vollst~_qdige Darstellung der turbulenten Geschwindigkeitsverteilung in glatten Leitungen, ZAMM, Bd. 31, Nr. 7 (1951) 208-219. [8,9] K. Rehme, KfK 2983 (1980); KfK 2617 (1978). [10,11] K. Rehme, KfK 3177 (1981); KfK 3318 (1982). [12,13] K. Rehme, KfK 3324 (1982); KfK 2512 (1977). [14] N. Neelen, Zur Thermohydraulik enger hexagonaler Brennstabbi~ndel, unpublished paper, K 8005, IfRR, Techn. Univ. Braunschweig (1980). [15] R. Harbord, Ein Programmpaket zur Ermittlung von N~herungsRSsungen f'flr statische und dynamische Probleme mit Hilfe der FEM, Berlin (1983). [16] J.G. Bartzis and N.E. Todreas, Hydrodynamic behavior of a bare rod bundle, Depart. of Nucl. Engrg., MIT, C00-2245-48TR (1977). [17] P. Carajilescov and N.E. Todreas, Experimental and analytical study of axial turbulent flows in an interior subchannel of a bare rod bundle, J. Heat Transfer 98 (1976) 262-268. [18] P. Kjellstr~Sm,Studies of turbulent flow parallel to a rod bundle of triangular array, Aktiebolaget Atomenergi, AE-487 (1974). [19] A.C. Trupp and R.S. Azad, The structure of turbulent flow in triangular array rod bundles, Nucl. Engrg. Des. 32 (1975) 47-84. [20] J.D. Hooper, Developed single phase turbulent flow through a square-pitch rod cluster, Nucl. Engrg. Des. 60 (1980) 365-379. [21,22] K. Rehme, KfK 2983 (1980); KfK 3069 (1980). [23] W. Zeggel and N. Neelen, Validation of wall parallel eddy viscosity formulation, paper presented at IAHR 5th Specialists Meeting on Liquid Metal Thermal Hydraulics, 1986, Grenoble, France. [24] K. Rehme, KfK 4027 (1986). [25] H.G. Kaiser, NRherungsRSsungen for den Impulstransport bei turbulenter StrSmung in engen Stabgitterbiindeln nach der Methode der Finiten Elemente, Thesis, Techn. Univ. Braunschweig (1986).