URANS simulation of the turbulent flow in tight lattice bundle

Abstract

The flow structure in tight lattice is still of great interest to nuclear industry. An accurate prediction of flow parameter in subchannels of tight lattice is likable. Unsteady Reynolds averaged Navier Stokes (URANS) is a promising approach to achieve this goal. The implementation of URANS approach will be validated by comparing computational results with the experimental data of Krauss. In this paper, the turbulent flow with different Reynolds number (5000–215000) and different pitch-to-diameter(P/D) (1.005–1.2) are simulated with computational fluid dynamics (CFD) code CFX12. The effects of the Reynolds number and the bundle geometry (P/D) on wall shear stress, turbulent kinetic energy, turbulent mixing and large scale coherent structure in tight lattice are analyzed in details. It is hoped that the present work will contribute to the understanding of these important flow phenomena and facilitate the prediction and design of rod bundles.

Keywords

tight rod bundle / flow structure / unsteady Reynolds averaged Navier Stokes (URANS)

Introduction

Tight lattice fuel assemblies have been proposed for advanced reactors. The fuel utilization will be enhanced by decreasing the pitch-to-diameter (P/D), i.e., less coolant volume fraction in the core which results in less moderation assures harder neutron energy spectrum and leads to higher conversion of ²³⁸U to ²³⁹Pu [1-3].

Early experimental observations on turbulent flow in rod bundle had been conducted in the early 1960s. With the development of measurement techniques, more experiments were presented in the late 1990 [4-9], which indicate that the turbulent flow in a rod bundle has completely different characteristics compared to that in a pipe. A high mixing in the gap region was observed, which was once explained by the secondary flow. However, later experiments prove that secondary flow is not the major factor for high mixing. The so called flow pulsation phenomenon was responsible for this high mixing. Vortices are transported in the longitudinal quasi-periodically with this oscillating flow. The interactions between transported vortices result in a gain in the momentum transfer and an increase in the mixing. These flow oscillations depend highly on the configuration of subchannels and the Reynolds number. This phenomenon presents a Reynolds threshold below which no actual oscillation is observed [10]. It was found that the pulsation frequency in a rod-wall gap decreases measurably as the gap size (width-to-depth ratio (W/D)) decreases in the range 1.015≤W/D≤1.250 [11]. Although a complete understanding of these oscillations has not yet been achieved, flow instability has mostly been accepted as the origin of these oscillations.

Up to the present, various approaches, such Reynolds averaged Navier-Stokes (RANS), unsteady RANS (URANS), large eddy simulation (LES), and direct numerical simulations (DNS) have been used in numerical investigations to study rod bundle flows [12-16]. Theses works show that RANS with isotropic turbulence models will miss the anisotropy in the rod bundle. Although the anisotropic turbulence models are capable of reproducing the turbulence-driven secondary flows in subchannels, the secondary flow is less dominant in the case of tightly packed geometries. The URANS approach captures the flow oscillation in the tight lattice rod bundle so that the accuracy in the prediction of averaged statistics is achieved since the wavelength of the oscillations has been grossly over-predicted. Of all these approaches, DNS is the most preferred because the global flow pulsation phenomena are still unknown or not clearly understood. LES has reproduced almost identical results with DNS for the turbulent flows in the rod bundle. Given the computational cost of the DNS and LES, URANS is a general practical approach capable of challenging arbitrary fuel rod-bundle design.

This paper focuses on the simulation of the turbulent flow inside different subchannels with RANS and URANS. The effect of the turbulence model on these simulations is investigated systemically. The features of the coherent structure in tight lattice with different P/D and different Reynolds number are studied in detail. This paper shows the effect of P/D and Reynolds number on turbulent flow in tight lattice. The simulation proves that there exists a critical P/D for specified Reynolds number, blow which the coherent structure disappear.

Numerical procedure

Numerical results were validated by the experiment of Krauss and Meyer [9]. Detailed information of experimental and numerical setup is listed in Table 1.

The experiment of Krauss and Meyers was conducted in a rod bundle of 37 parallel rods (O.D. 140 mm) arranged in a triangular array built in a hexagonal symmetric horizontal channel. The total length of the working section is L=11.50 m. The measurements of wall shear stress, axial velocity and turbulent intensity in tight lattice are performed in the experiment. The Reynolds stress model will be applied to the RANS&URANS simulation.

The computational domain consists of two sub-channels connected by a narrow gap. The boundary conditions include three couples (A, B, C) of periodic boundaries (Fig. 1) and non-slip walls. In the present case the computational length has been chosen to be equal to four times the average streamwise wavelength λ, obtained from the experiment (0.6 m).

The mesh presented in Fig. 1 has been used, for a total of more than 600000 meshes. Y+ is smaller than 20. In any case, the time step size has been ensured to satisfy

\mathrm{\Delta }t\ll {\displaystyle \frac{\mathrm{1}}{f}},

where f is the smallest dominant frequency that can be observed a posteriori as the simulation develops. Δt has been chosen equal to 10^-4 s.

Results and discussion

Figure 2 illustrates the comparison between experimental data and simulation results from URANS and RANS, where W is the axial velocity, W_b is the bulk velocity, y is the distance from the wall, y_max is the distance from the wall where the velocity reaches the maximum, τ is the wall shear stress, τ_m is the mean wall shear stress, k⁺ is the dimensionless turbulent intensity normalized by the average turbulent intensity, and \varphi is the azimuthal angle. The URANS simulation significantly improved the accuracy so that it is credible for prediction of the turbulent flow in tight lattice.

In order to investigate the effect of Reynolds number and P/D on turbulent flow in tight lattice, the flow behavior with Reynolds numbers ranging from 5000 to 215000 and P/D ranging from 1.001 to 1.2 is simulated in triangular array.

Wall shear stress

The homogeneity of the wall shear stress increases with the increase of P/D (Fig. 3(a)). The maximum wall shear stress appears in the widest flow region, while the maximum wall shear stress has obvious drift in DNS simulation [16]. The URANS simulations miss this monotonic trend which is also found in experiment and not fully understood yet. The homogeneity of the wall shear stress is poor for low Reynolds number (Fig. 3(b)). The wall shear stress is not sensitive to the Re for high Reynolds number.

Turbulent kinetic energy

The definition of the relative kinetic energy of none coherent structure is

kn{c}^{+}={\displaystyle \frac{u^{\mathrm{2}}+v^{\mathrm{2}}+w^{\mathrm{2}}}{\mathrm{2}{u}_{\tau }^{\mathrm{2}}}},

where u,v,w are axial, radial and azimuthal fluctuation velocity, respectively. u_{\tau } is shear velocity. The resolved velocity fluctuation is identified as coherent, and the solutions of the Reynolds stress equations are identified as non-coherent. The kinetic energy of the coherent velocity fluctuation is obtained by time averaging the sum of the squares of the resolved velocity component fluctuations as

k{c}^{+}={\displaystyle \frac{\overline{{\left(U_x-\overline{U_x}\right)}^{\mathrm{2}}}+\overline{{\left(U_y-\overline{U_y}\right)}^{\mathrm{2}}}+\overline{{\left(U_z-\overline{U_z}\right)}^{\mathrm{2}}}}{\mathrm{2}{u}_{\tau }^{\mathrm{2}}}},

where U_x,U_y,U_z is the transient axial, radial and azimuthal velocity. The total time-averaged turbulent kinetic energy per unit mass is determined as the sum of the above terms

k^{+}=k{c}^{+}+kn{c}^{+}.

Figure 4 demonstrates the effect of P/D and Re on turbulent kinetic energy at 30° azimuthal angle. With the decrease of P/D, the effect of the coherent structures becomes more significant, which leads to greater contribution to the kinetic energy far away from the wall (Fig. 4(a)). The change law of k^{+} with the increase of Re is similar to that of kn{c}^{+} (Fig. 4(b). k^{+} becomes larger and closer to the wall when Re becomes larger. Furthermore, it is not sensitive to Re when Re reaches some critical value. This suggests that the increase rate of the coherent structure scale decreases with the increase of Reynolds number.

Turbulent mixing

For the present geometry, let us denote the two subchannels adjacent to the gap as 1 and 2, and assume that the bulk temperatures in the two subchannels are T_b1 and T_b2, respectively. Following Ref. [17], the convective heat transfer rate between these subchannels can be expressed as

Q=\rho C_p{w}_{\mathrm{e}\mathrm{f}\mathrm{f}}{\delta }_{\mathrm{12}}L(T_{\mathrm{b}\mathrm{1}}-T_{\mathrm{b}\mathrm{2}}),

where w_{\mathrm{e}\mathrm{f}\mathrm{f}} is an effective mixing velocity across the gap. It can be expressed in terms of an eddy viscosity {\upsilon }_{\mathrm{T}}, a mixing distance {\delta }_{\mathrm{12}} between the two subchannels, and an empirical mixing factor Y, which accounts for the subchannel shape, as

w_{\mathrm{e}\mathrm{f}\mathrm{f}}=Y{\displaystyle \frac{{\upsilon }_{\mathrm{T}}}{{\delta }_{\mathrm{12}}}}.

Further utilizing the empirical relationship in Ref. [17]

{\upsilon }_{\mathrm{T}}=\mathrm{0.0177}\upsilon \mathit{R}\mathit{e}\sqrt{f_{\mathrm{t}}},

where \upsilon is the kinematic viscosity of the fluid, f_{\mathrm{t}}=\mathrm{0.18}\mathit{R}{\mathit{e}}^{-\mathrm{0.2}} is the friction factor for smooth circular tubes, and specifying {\delta }_{\mathrm{12}} by geometrical reasoning. Moller [18] suggests to specify {\delta }_{\mathrm{12}} as the distance between the center of two subchannels (\sqrt{\mathrm{3}}P/\mathrm{3}) for tight lattice. Y needs to be determined from empirical information available. Rehme [17] proposes the empirical relationship as follow

Y={\displaystyle \frac{C}{\delta /D}}.

The empirical relationship predicts the increase of Y with the decrease of \delta /D. In this paper C=2.58 is obtained through Least square method. But the relationship means infinite mixing with zero gaps, which is in contradiction to common sense. Figure 5 depicts the effect of P/D on mixing factor.

When P/D>1.001, the numerical results agree well with the relationship. Apparently the relationship miss the disappearance of the mixing with P/D=1.001. No relationship is capable of predicting that the mixing is ignorable when P/D reaches some critical value.

Figure 6 shows the effect of P/D on mixing factor. Y is 41.61 from the relationship. The numerical results are within the error of 11.2%.

Coherent structures

The flow oscillation is caused by the coherent structure, as proposed by Jeong and Hussain [19]. The parameter Q is introduced to identify this coherent structure and is defined as the second invariant of the velocity gradient tensor:

Q=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial U_i}{\partial x_j}}{\displaystyle \frac{\partial U_j}{\partial x_i}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(S_{ij}{S}_{ij}+{\mathrm{\Omega }}_{ij}{\mathrm{\Omega }}_{ij}),

where the strain rate tensor and rotation rate tensor are defined as

S_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial U_i}{\partial x_j}}+{\displaystyle \frac{\partial U_j}{\partial x_i}}\right),

and

{\mathrm{\Omega }}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial U_i}{\partial x_j}}-{\displaystyle \frac{\partial U_j}{\partial x_i}}\right),

respectively. Positive values of Q indicate regions where vorticity overcomes strain. The present study fails to identify the coherent structure by Q. Therefore, a modified Q_m was introduce as

Q_{\mathrm{m}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(C_{\mathrm{q}}{S}_{ij}{S}_{ij}+{\mathrm{\Omega }}_{ij}{\mathrm{\Omega }}_{ij}),

where C_q<1 is an empirical factor which reduce the weight of strain effect. By selecting C_q=0.55 [14] and Q_m=1, it becomes possible to identify the surface of the coherent structure. Figure 7 presents the coherent structure identified in the tight lattice (Re=28500). It is observed that the coherent structure appears in pairs on either side of the gap.

The contour of Q_m on plane YOZ with different P/D values is exhibited in Fig. 8. With the increase in P/D, the scale and the configuration of the coherent structure become smaller and more irregular. The coherent structure is not obvious at P/D=1.2.

The contour of Q_m on plane YOZ with different Reynolds numbers is given in Fig. 9. It is observed that the scale of the coherent structure appears in pairs and increases with the increase in Reynolds number. The coherent structure exists even with low Reynolds number (5000). It was pointed out that there exists threshold value of Reynolds number, below which no oscillation occurs [20]. This phenomenon is also validated in this numerical study.

Conclusions

The validity of the methodology used in this study is based on the experimental data of Krauss and Meyer [9].With the Reynolds number ranging from 5000 to 215000 and P/D from 1.001 to 1.2, the effect of the Reynolds number and the bundle geometry on the flow oscillation is investigated.

In very tight lattice (P/D<1.1), the effect of P/D on the wall shear stress and turbulent kinetic energy is significant due to the dramatic variation of the amplitude and the frequency of the flow oscillation in the gap region.

This paper verifies the inverse ratio between mixing factor Y and geometric factor before P/D is above critical value.

For the fixed geometry (P/D), the flow parameter. i.e., stream wise velocity, wall shears stress, turbulent kinetic energy is not sensitive to the Re when Re is higher than some value (9600 in this paper).

The scale of the coherent structure increases when the Reynolds number increases or the P/D decreases. There exists a critical P/D for specified Reynolds number, blow which the coherent structure disappears.

References

[1] Oldekop W, Berger H D, Zeggel W. General features of advanced pressurized water reactors with improved fuel utilization. Nuclear Technology, 1982, 59(2): 212-227

[2] Uchikawa S, Okubo T, Kugo T, Akie H, Nakano Y, Onuki A, Iwamura T. Investigations on innovative water reactor for flexible fuel cycle (LFWR). In: Proceedings of GLOBAL, Tsukuba, Japan, 2005, Paper No. 358

[3] Cheng X, Liu X J, Yang Y H. A mixed core for supercritical water-cooled reactors. Nuclear Engineering and Technology, 2008, 40(1): 1-10

[4] Rehme K. Simple method of predicting friction factors of turbulent flow in non-circular channels. International Journal of Heat and Mass Transfer, 1973, 16(5): 933-950

[5] Trupp A C, Azad R S. The structure of turbulent flow in triangular array rod bundles. Nuclear Engineering and Design, 1975, 32(1): 47-84

[6] Trippe G, Weinberg D. Non-isotropic eddy viscosities in turbulent flow through rod bundles. In: Kakac S, Spalding D B, eds. Turbulent Forced Convection in Channels and Bundles, Vol. 1. Washington: Hemisphere Publishing Corporation, 1979

[7] Seale W J. Turbulent diffusion of heat between connected flow passages Part 1: Outline of problem and experimental investigation. Nuclear Engineering and Design, 1979, 54(2): 183-195

[8] Rehme K. The structure of turbulent flow through rod bundles. Nuclear Engineering and Design, 1987, 99(1): 141-154

[9] Krauss T, Meyer T. Experimental investigation of turbulent transport of momentum and energy in a heated rod bundle. Nuclear Engineering and Design, 1998, 180(3): 185-206

[10] Meyer L, Rehme K. Large-scale turbulence phenomena in compound rectangular channels. Experimental Thermal and Fluid Science, 1994, 8(4): 286-304

[11] Baratto F, Bailey S C C, Tavoularis S. Measurements of frequencies and spatial correlations of coherent structures in rod bundle flows. Nuclear Engineering and Design, 2006, 236(17): 1830-1837

[12] Baglietto E, Ninokata H. Turbulence models evaluation for heat transfer simulation in tight lattice fuel bundles. In: Proceedings of the 10th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-10), Seoul, Korea, 2003

[13] Chang D, Tavoularis S. Unsteady numerical simulations of turbulence and coherent structures in axial flow near a narrow gap. ASME Journal of Fluids Engineering, 2005, 127(3): 458-466

[14] Chang D, Tavoularis S. Numerical simulation of turbulent flow in a 37-rod bundle. Nuclear Engineering and Design, 2007, 237(6): 575-590

[15] Merzari E, Ninokata H, Baglietto E. Large eddy simulation of the vortex street between rectangular channels. In: NTHAS 5, Jeju Island, Korea, 2006

[16] Ninokata H, Merzari E, Khakim A. Analysis of low Reynolds number turbulent flow phenomena in nuclear fuel pin subassemblies of tight lattice configuration. Nuclear Engineering and Design, 2009, 239(5): 855-866

[17] Rehme K. The structure of turbulence in rod bundles and the implications on natural mixing between the subchannels. International Journal of Heat and Mass Transfer, 1992, 35(2): 567-581

[18] Moller S V. Single-phase turbulent mixing in rod bundles. Experimental Thermal and Fluid Science, 1992, 5(1): 26-33

[19] Jeong J, Hussain F. On the identification of a vortex. Journal of Fluid Mechanics, 1995, 285(1): 69-94

[20] Lexmond A S, Mudde R F, van der Haagen T H J J. Visualization of the vortex street and characterization of the cross flow in the gap between two subchannels. In: Proceedings of the 11th Nureth, Avignon, France, 2005

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg