Annals of Nuclear Energy 171 (2022) 109055

Contents lists available at ScienceDirect

Annals of Nuclear Energy

journal homepage: www.elsevier.com/locate/anucene

Numerical investigation on flow-induced vibration and heat transfer of

fuel rods in a small pitch-diameter ratio (P/D) channel
Yonghui Guo a, Xiaochang Li a,⇑, Ju Liu a, Guangliang Chen b, Ye Gao a
a
College of Aerospace and Civil Engineering, Harbin Engineering University, Harbin 150001, China
b
College of Nuclear Science and Technology, Harbin Engineering University, Harbin 150001, China

a r t i c l e i n f o a b s t r a c t

Article history: The Korea Atomic Energy Research Institute developed a dual-cooled annular fuel rod bundle structure to
Received 3 August 2021 enhance the core power of the nuclear reactor and improve its heat transfer efficiency. The structure also
Received in revised form 6 February 2022 reduced the pitch-diameter ratio (P/D) of the channel. In this study, a four-subchannel model is taken as
Accepted 28 February 2022
the research object, and the Reynolds stress model (RSM) is applied to calculate the flow-solid-thermal
Available online 7 March 2022
tri-physical field coupling for a typical bare rod beam channel. The accuracy of the numerical calculation
method used in this study was first verified by comparing it with the experimental data in the literature.
Keywords:
Then, the channel flow field under different P/D conditions was simulated, and the formation conditions
Flow-induced vibration
Gap vortex street
and influencing factors of the vortex-street structure were studied. Axial flow-induced vibrations of a
Heat transfer enhancement cylindrical structure in a channel were investigated using a weak coupling algorithm. Finally, the influ-
RSM model ence of the vortex-street structure on the heat transfer efficiency of the fluid in the channel is discussed.
The results show that the anisotropy-based RSM turbulence model can accurately capture the vortex-
street structure in small P/D channels. When P/D < 1.050, a stable vortex-street structure appeared in
the channel, which led to the stable micro-amplitude vibration of the cylindrical structure in the channel.
The presence of the vortex-street structure significantly improves the heat transfer efficiency of the fluid
in the channel.
Ó 2022 Elsevier Ltd. All rights reserved.

1. Introduction with the conventional fuel rod (As shown in Figure a), the dual-
cooled annular fuel rod bundle structure (As shown in Figure b)
The contradiction between the rapid growth of energy demand has two cooling channels, which greatly increases the contact area
and the energy shortage has become one of the main problems that between the fuel rods and coolant; however, the rod bundle chan-
the world is dealing with. Research shows that by 2030, the world’s nel is narrower owing to the increase in the rod bundle diameter.
energy consumption will increase by 35% to 49% compared with Fig. 1.2 shows the structural diagram of the rod bundle subchannel,
2010 (Galiana Gonzalez, 2010). The development and utilization where D is the diameter of the fuel rod, and P (pitch) represents the
of new clean energy will not only ease the increasingly serious distance between the centers of two adjacent fuel rods. The com-
energy crisis, but also reduce global warming to a certain extent. pactness of fuel rod bundles can be described by the pitch-
In the mid-20th century, the utilization of nuclear energy began diameter ratio (P/D). The literature (Guellouz and Tavoularis,
to be viewed in this context. At present, various types of nuclear 2000a) shows that when the P/D is small enough, a large-scale flow
reactors are widely used in electric power, heating, and power fluctuation occurs between rod bundle channels, called rod bundle
engineering. vortex network. The appearance of a vortex-street structure signif-
The Korea Atomic Energy Research Institute developed a dual- icantly improves the heat exchange efficiency between rod bundle
cooled annular fuel rod bundle structure to enhance the core channels, but large-scale flow pulsations will also cause the fuel
power of the nuclear reactor to improve its power and heat trans- rods to vibrate (De Ridder et al., 2016). Therefore, it is very impor-
fer efficiency (Chun et al., 2009; Shin et al., 2012). Fig. 1.1 shows a tant to study the vortex-street phenomena.
comparison between the structure of the dual-cooled annular fuel Tapucu et al. were the first to investigate vortex structures
rod bundle and conventional fuel rod (Lee et al., 2013). Compared (Tapucu and Merilo, 1977). They studied the transverse flow
between two parallel channels connected by a narrow transverse
⇑ Corresponding author. slot, and observed a sinusoidal pattern of fluctuations in both the
E-mail address: lixiaochang@hrbeu.edu.cn (X. Li). flow and pressure. Since then, a large number of experimental

https://doi.org/10.1016/j.anucene.2022.109055
0306-4549/Ó 2022 Elsevier Ltd. All rights reserved.
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

eet phenomena.

Coolant Cladding Cladding

Pellet

Coolant Pellet

Gap Gap
Coolant
(a) (b)

Fig. 1.1. Structural diagram of (a) the dual-cooled annular rod bundle channel and (b) conventional fuel rod.

studies have been conducted on this phenomenon (Guellouz and duced the experimental results most accurately. However, the
Tavoularis, 2000a,b; Möller, 1991; Möller, 1992). Meyer and results of unsteady Reynolds-averaged Navier–Stokes (URANS)
Rehme proposed the flow structure in a rod bundle in 1994 using simulation using the RSM are also reliable, and the flow field can
all of the information available at that time. They reported ‘‘Two be predicted with low computational cost. A large number of stud-
vortices driven by the higher velocities outside the gap and rotat- ies have proved that (Házi, 2005; Merzari et al., 2008; Horváth and
ing in opposite directions are transported with the velocity Uc (ap- Dressel, 2012; Chang and Tavoularis, 2007) owing to the strong
proximately the average of the velocity at the center and the edge anisotropic characteristics of the fluid in the rod bundle channel,
of the gap) axially transport within the gap” (Meyer et al., 1994). the RSM model can better predict the turbulent variables in the
The schematic is shown in Fig. 1.3. Meyer et al. provided a system- channel than the traditional two-equation model.
atic summary of previous studies and reviewed their main findings The flow-induced vibrations caused by vortex structures is
(Meyer, 2010). another topic of research interest. The micro-amplitude vibrations
However, the experimental measurements were costly, restric- of a cylinder in an axial flow are often caused by random pressure
tive, and the detailed laws of the flow field were difficult to study. fluctuations in the flow field around the structure. The causes of
With the rapid development of computer technology, numerical pressure fluctuations can be divided into two categories: the
calculations gradually gained mainstream acceptance in research. near-field component and far-field component of the pressure field
For example, Tavoularis (Chang and Tavoularis, 2005) numerically (De Santis and Shams, 2017). The near-field component is gener-
simulated the axial flow in a rectangular channel with a cylindrical ally caused by local pressure fluctuations in the turbulent bound-
rod by solving the unstable Reynolds mean Navier–Stokes equation ary layer, whereas the far-field component is attributed to
using the Reynolds stress model (RSM). It was observed that coher- disturbances in the flow field, such as upstream obstacles, pipe
ent structures existed in the channel as a pair of reverse vortices bending, or pump pulsation (Païdoussis and Gagnon, 1984). Thus,
and flowed along the axial direction. Tavoularis (Chang and the factors that induce structural vibrations in the flow field can
Tavoularis, 2012) carried out three-dimensional unsteady simula- be divided into two categories: the pressure fluctuations in the tur-
tions of developing turbulence in rectangular channels with a bulent boundary layer, and the disturbances caused by the macro-
cylindrical rod, and studied their sensitivity to boundary condi- scopic structure. It is well known that the URANS model calculates
tions and turbulence model selection. Among all the methods used the average of a set of random quantities of fluid flow; fluctuations
in the authors’ paper, the large-eddy simulation (LES) model repro- such as turbulent velocity and pressure are not calculated. There-
fore, the standard URANS equation cannot predict the near-field

P
Fig. 1.2. Schematic diagram of the rod bundle subchannel. Fig. 1.3. Flow model for the gap region.

2
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

component of the pressure field. ter Hofstede et al. (2017), De selected as the research object, i.e., the four-subchannel model
Ridder et al. (2013), De Santis and Shams (2017), and De Santis combined with the periodic boundary conditions.
and Shams (2019) examined the feasibility of the URANS model Fig. 2.2 shows the schematic diagram of the four-subchannel
in fluid-structure interaction simulation involving turbulence. geometric model used in this study. The definitions of D and P
However, their results did not represent the actual behavior of are the same as those in the introduction, and L is the length of
structures observed in turbulence because there was no excitation the fuel rod. To determine the specific parameters of the fuel rods,
caused by pressure fluctuations. Because of the absence of external the geometric dimensions used in the study of De Santis and
excitation, the vibrations generated by the URANS model were not Shams (2019) were used, which are the same as those used in test
sustained. Therefore, we can only use the URANS model to simu- facilities built by SCK-CEN (De Pauw et al., 2016). The diameter of
late the fluid-structure coupling phenomenon caused by macro- the fuel rod is the same as that of the MYRRHA reactor (Baeten
scopic disturbances. et al., 2014); and its length is half of the actual value. According
The study of heat transfer efficiency in the rod beam channel is to Chang and Tavoularis (2012), the size of the computational
an important topic in reactor research. Yao et al. (1982) deter- domain length L significantly affects the formation of the vortex
mined the experimental correlation of the heat transfer coefficient structure and the accuracy of the numerical calculation results.
of the rod bundle channel based on the temperature measurement Although Guellouz and Tavoularis (2000a,b) used an experimental
results in a rod bundle channel with and without a grid. Holloway L of only 50 D, choosing the same L value in the experiment will
et al. (2004), Holloway et al. (2005), Holloway et al. (2008) con- cause significant errors because of the difficulty in providing the
ducted a series of experiments on the pressure drop and Nusselt real inlet and outlet conditions in the numerical calculation. The
number distribution laws for bare beam channel models with authors recommend that the length of the computational domain
and without a vane, and with different types of vanes. They be chosen to be 108D to allow the flow to approach full develop-
showed that the cross-mixed vane grid has a significant effect on ment at the exit of the channel. Integrating the above research
both the axial and circumferential heat transfer of the bar beam, results, the final D and L values were determined as 6.51 mm
and obtained the related heat transfer correlation equations. Over and 700 mm (107.5 D), respectively. To facilitate the subsequent
the past decades, the computational fluid dynamics (CFD) numer- analysis, this study defines the narrow-slit region between two col-
ical study of single-phase flow and heat transfer in rod bundle umns as the gap region, and the channel center region as the
channels has made significant progress. The application of CFD to subchannel.
predict the single-phase flow field of a reactor fuel assembly has
proved to be a mature and effective method. Liu et al. (2012) com-
2.2. Mesh setup
pared the effects of different turbulence models and near-wall
treatment methods on the numerical results from the perspective
The overall schematic diagram of the fluid domain mesh is
of convective heat transfer. By comparing the numerical results
shown in Fig. 2.3. A hexahedral mesh was used to partition the
and the experimental data, they showed that the SST k-x model
computational domain for the trade-off between computational
was more adaptable to different types of near-wall grids than the
accuracy and computational cost. Because the flow field calculated
k
e series and RSM model. In et al. (2004) used the anisotropic
in this study was strongly anisotropic, the RSM model must be
turbulence model to perform CFD numerical analysis on the flow
applied. There was a large velocity gradient near the wall, espe-
and heat transfer in a triangular bare bar bundle with different
cially in the gap region; hence, the two-layer wall model was used.
diameter-distance ratios. They showed that the anisotropic turbu-
The mesh was refined near the walls of the cylindrical structure to
lence model significantly improved the accuracy of the predicted
ensure that y+ was approximately 1. The local grid refinement dia-
time-averaged velocity and temperature in the channel compared
gram is shown in Fig. 2.4.
with the standard linear model.
The grid of the cylindrical structure is shown in Fig. 2.5. Accord-
In this study, the RSM model is applied to calculate the coupled
ing to the research of De Santis and Shams (2017), selecting the
flow-solid-thermal tri-physical field of a typical bare rod-bundle
type of structural grid element is more important than the number
channel. The effects of the P/D ratio on the formation of the
of cells. Therefore, in this study, quadratic elements were selected
vortex-street structure, the fluidic vibration of the fuel rod, and
for the calculation. To reduce the data exchange error at the cou-
the heat transfer efficiency of the channel flow field are investi-
pling interface, the numbers of structural elements and fluid ele-
gated. This paper is organized as follows. Section 2 presents the
ments in the axial direction were kept as consistent as possible.
geometric model and numerical methods. Section 3 verifies the
The number of cells in the solid domain cross section is 289, and
reliability of the numerical algorithm. Section 4 analyzes the calcu-
the number of axial cells is the same as that of the fluid domain.
lation results, and Section 5 concludes.
Detailed information about the meshing of the fluid domain can
be found in Section 3.3.

2. Geometric model and numerical methodology 2.3. Numerical algorithms

2.1. Geometric model Simulations were performed using the commercial CFD code
ANSYS 19.0 and ANSYS Fluent was used to solve the flow field
The size of an actual fuel assembly rod bundle is usually large; problems. Transient structural software for the Workbench plat-
for example, the typical fuel assembly of a pressurized water reac- form was used to solve the structural dynamics. Data exchange
tor is 17  17, as shown in Fig. 2.1. In published experimental stud- was implemented through the Workbench platform to ensure
ies, 3  3 or 5  5 fuel assemblies were adopted. If 17  17 or even equilibrium of forces and displacements at the fluid-solid coupling
5  5 rod bundle channels are directly taken as the research object, interface. To ensure the convergence of the iterations, the coupling
then the number of cells will be huge, which is not conducive to algorithm used in the calculations is the Quasi-Newton Stabiliza-
the research. Because the research object of this study is a bare tion Algorithm with a convergence residual of 0.001. The calcula-
channel, the flow field in the radial direction has obvious periodic- tions in this paper are first performed for the flow and
ity. To effectively reduce the computations, the smallest unit that temperature fields, and then the fluid-solid coupling is started
can reflect the flow field and heat transfer characteristics was when the flow and temperature fields reach steady state.
3
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Minimum repetitive element

Fig. 2.1. Minimum repetitive element of the fuel assembly with a hybrid wing lattice.

P
Gap Region
Subchannel

Fig. 2.2. Schematic diagram of four-channel geometric model.

Fig. 2.3. The overall schematic diagram of the fluid domain mesh.

According to the research of Chang and Tavoularis (2012), the approximately equal to 1. The results of literature (Liu et al.,
flow field with a vortex-street structure is strongly anisotropic; 2012) show that when the RSM turbulence model is used in com-
hence, the one-equation or two-equation model based on the bination with the y+ approximately equal to 1 grid model, the
Boussinesq hypothesis cannot provide reliable simulation results. calculated Nu number agrees well with the experimental results.
Although the LES model and RSM model can provide relatively The fuel rod boundary conditions were used with a constant heat
accurate numerical results, after a comprehensive consideration flow density (707,000 W/m2). Considering the small temperature
of the calculation costs and accuracy, the RSM-e model was difference between the inlet and outlet of the calculated condi-
selected as the turbulence model for the flow field simulation. tions, the fluid properties related to heat transfer are calculated
The ALE N–S have counted the effect of mesh motion in flow field, using the constant properties at the inlet temperature (310 K)
as the mesh motion is usually generated based on the common and pressure (0.1 MPa). Since the heat transfer mainly considers
movable boundaries between fluid and structure domains. the fluid temperature field and the heat transfer phenomenon
The two-layer model (called Enhanced Wall Treatment in Flu- between the fluid and the wall, and the boundary condition is a
ent) is used for the convective heat transfer model near the wall, constant heat flow density, the solids are not involved in the heat
and the first grid node of the wall is ensured to be located at y+ transfer calculation.

4
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Fig. 2.4. Local schematic of the fluid domain cross-sectional mesh.

Fig. 2.5. Grid distributions on the cylinder surface.

The fluid-solid coupling algorithm can be divided into direct the selection of the grid cell length Dz, please see the subsequent
solution algorithm and partitioned iterative solution algorithm section on grid correlation verification. The final time step chosen
according to the different ways of solving the control equations. in this paper is 2.5E
5s.
Compared with the direct solution algorithm, the partitioned iter- The boundary conditions of the dynamic calculation domain are
ative solution algorithm can make maximum use of the existing shown in Fig. 2.8. The fixed-constraint boundary conditions were
mature solvers and has a wider application range. Therefore, this applied to the two ends of the cylinder. The cylinder surface was
paper finally adopts the partitioned iterative solution algorithm. set as the fluid-structure interaction interface. The Newmark algo-
Fig. 2.6 shows a schematic diagram of the boundary conditions rithm based on time integration was used to solve the structural
of the fluid domain. The left side is the velocity inlet boundary, the dynamics. The Newmark method is a method of numerical integra-
right side is the pressure outlet boundary, and all the fuel rod walls tion used to solve differential equations. It is widely used in
are no-slip walls. Because the selected four-subchannel model is numerical evaluation of the dynamic response of structures and
the smallest repeating unit, the side boundaries in the figure were solids such as in finite element analysis to model dynamic systems.
set as periodic boundary conditions. The specific periodic boundary The time step used for the calculation is consistent with the fluid
conditions are shown in Fig. 2.7. Detailed information on the peri- domain. The material used for the dynamic calculation was struc-
odic boundary conditions can be found in the literature (Li and tural steel with a density qs of 7850 kg/m3, Young’s modulus of
Gao, 2014). The Semi-Implicit Pressure Linked Equations Consis- 200 GPa, and Poisson’s ratio ts of 0.3.
tent (SIMPLEC) algorithm was used for the pressure and velocity
coupling. The pressure term was discretized in a second-order
scheme; the second-order upwind scheme was used to discretize 3. Numerical algorithm validation
the momentum, turbulent kinetic energy, turbulent dissipation
rate, and Reynolds stress. The second-order implicit scheme was 3.1. Validation of numerical algorithm for flow field
used to discretize the time term.
The fluid is water (20 °C) with a density q of 998.2 kg m3 and a The numerical example in Chang and Tavoularis (2012) was
kinetic viscosity l of 0.001003 m2/s. The inlet velocity U b of water simulated to verify the reliability of the numerical algorithm in this
is 5 m/s, and the hydraulic diameter is 0.0035 m. The calculation study. The working conditions calculated by Chang et al. are the
time step was recommended by Chang and Tavoularis (2012), that
P1 P2
is, the Courant–Friedrichs–Lewy value (U b Dt=Dz) is equal to 0.15,
which ensures the stable convergence of the computations. For

P3 P3
Y

Wall X

Outlet
P4 P4
Inlet

Periodic P1 P2

Fig. 2.6. Boundary conditions of the fluid domain. Fig. 2.7. Distribution diagram of periodic boundary conditions.

5
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Fixed

Fluid-solid interface

Fixed

Fig. 2.8. Boundary conditions of the dynamic calculation domain.

fpD
St (3-1)
Ub

2D
D
W
y

x
3D

Fig. 3.1. Schematic diagram of the geometric parameters of the cross-section of the computational model (left) and schematic diagram of the mesh (right).

Fig. 3.2. Comparison of numerical results (right) and experimental results (left) of the axial velocityV z .

experimental results of Guellouz and Tavoularis (2000a). The geo- in Eq. (3-1), where is the frequency of the transverse velocity, is
metric model and boundary conditions used in the calculations are the diameter of the cylinder, and is the axial velocity. represents
the same as those in the original paper. The geometric model and the axial velocity of flow in the gap centre. represents the average
mesh of the computational domain are shown in Fig. 3.1. The cal- wavelength. It was found that the value and position of the maxi-
culation results of this paper are compared with those of Chang mum axial velocity predicted by the numerical method were very
and Tavoularis (2012) paper using RSM model. close to the experimental results. The velocity contour distribution
Fig. 3.2 shows the comparison between the numerical results of at the corner and gap was also consistent. Table 3.1 presents the
this study and the flow field of the experimental results in comparison of the present numerical results with those of Chang
Guellouz et al. The definition of the Strouhal number () is shown et al. and the experimental results of Guellouz et al. The errors of

6
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Table 3.1
Comparison of numerical results in this study and results in literature.

St U c =U b k=D
Guellouz (Experimental) 0.17 0.78 4.20
Chang (Numerical) 0.20 0.68 3.40
Current numerical results 0.21 0.65 3.71
Error with experimental results 23.5% 16.7% 11.7%
Error with numerical results 5.00% 4.41% 9.12%

the numerical results in this study and those of Chang et al. were
relatively close, i.e., within 10%. There were some differences
between the numerical and experimental results, with the maxi-
mum error of 23.5%. However, considering many factors such as
the gap between the model used in the numerical calculations
and the experimental model, and the fact that the authors did
not give specific entrance parameters, as well as the error of the
numerical calculation itself, the errors are considered to be within
acceptable limits. We believe that the comparison between the
present calculation results and the experimental results is consis-
tent, and the numerical algorithm is reliable. Fig. 3.3. Distribution curve of the Nusselt number.

f pD 3.2. Validation of fluid–structure coupling algorithm

St ¼ ð3-1Þ
Ub
De Santis and Shams (2019) simulated the fluid-induced vibra-
Because the analysis of convective heat transfer is included in
tion of a single cylinder in an axial flow. In the study by De Santis
this study, the reliability of the results of the heat transfer model
et al. fluid-structure interaction (FSI) simulations were performed
must be verified. The verification of convective heat transfer in a
for a single cylindrical rod in axial turbulent flow of water, analyz-
bare rod bundle is based on convective heat transfer experiments
ing the effect of the working fluid and the type of rod. The numer-
by El-Genk et al. (1993). The experiments were performed with a
ical example of De Santis was selected to verify the fluid-structure
square arrangement of rod bundle channels, where P/D = 1.38,
coupling algorithm in this study. All the geometric parameters
the medium was liquid water, the Reynolds number was 10000,
used in the validation were consistent with those reported in the
and the Prandtl number was 4.53. The experimental data for the
literature. The flow-solid coupling algorithm used in this paper is
relevant operating conditions are given in the literature as well
consistent with that in Section 2.3.
as the experimental correlation equations for the Nusselt number
The working condition simulated by De Santis in the paper is a
at different P/D ratios and Prandtl numbers. The geometric model
circular rod structure subjected to impulse loads to calculate the
and parameters used in the validation process are consistent with
vibration displacement response of the rod in water. The impact
Elgenk’s study, where the P/D was 1.194 and the Reynolds number
load is approximated by the following smooth function:
was 48000. The numerical results were validated using the follow-
ing experimental correlations: (
a  e
t2 =2r2 ; if 0  t  2:5  10
3 ;
dF ¼
0; otherwise ;
NuF;T ¼ CRe0:8 Pr0:333 ; C ¼ 0:028ðP=DÞ
0:006 ð3-2Þ
With r ¼ 0:0007 and load factor a chosen (after some prelimi-
nary tests) such that the maximum displacement of the center
where NuF;T is the fully developed Nusselt number; the velocity point of the rod is similar in the selected cases. The impulsive force
inlet boundary condition was given at the inlet; and the inlet tem- is applied both in positive y-direction and negative x-direction.
perature was taken as the temperature corresponding to the Prandtl Because no specific load factor a was given in the literature, the
number of 4.53 at the system pressure. The outlet was set as the final load factor was determined as 0.17, according to the selection
pressure outlet. All the walls were set as no-slip walls, and the principle of the load factor in the literature. Details on the method
two groups of side boundaries were set as periodic boundaries. of determining the load factor can be found in the paper by De
A comparison between the results of the numerical calculations Santis and Shams (2019).
and the above empirical formulation is shown in Fig. 3.3. Because Figs. 3.4 and 3.5 show the displacement of the cylinder center
the empirical formula is based on the measurement results under position over time. The current numerical calculation results and
the premise that both the flow and heat transfer are fully devel- the literature results are basically the same in terms of amplitude
oped, and the detailed inlet conditions have not been given in and vibration frequency, but the error gradually increased after
the literature, it is difficult to avoid the influence of the inlet effect 0.3 s, owing the time integration algorithm used for the vibration
in the calculation. As can be seen in the figure, the current numer- displacement solution; hence, the error accumulated with the
ical results underwent a stage of drastic changes at the inlet, owing increase in time. Meanwhile, because the specific load parameters
to the velocity inlet boundary conditions that caused the inlet to go were not given in the literature, the parameters used in this study
through a fully developed stage. At locations away from the inlet, were obtained by trial calculations according to the instructions in
the fluid and heat transfer reached a fully developed state, and the literature; thus, the above errors are considered to be accept-
the calculated Nusselt number gradually approximated the result able. The first, second, and third mode vibration frequencies of
of the empirical formulation. Therefore, it is considered that the the cylinder and the comparison with the results predicted by De
current model used for heat transfer calculations can effectively Santis are presented in Table 3.2. The results indicate that the
respond to real physical phenomena. errors between the numerical results and literature results are all
7
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

16

12 Santis(2019)
Current Numerical Results

Y-Displacement(mm)
8

-4

-8

-12

-16
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
Time(s)
Fig. 3.4. Comparison of current numerical results and literature results: Y-displacement.

16

12 Santis(2019)
Current Numerical Results
X-Displacement(mm

-4

-8

-12

-16
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
Time(s)
Fig. 3.5. Comparison of current numerical results and literature results: X-displacement.

Table 3.2 Y
Comparison of current numerical results and literature results: frequency.

f 1 ðHzÞ f 2 ðHzÞ f 3 ðHzÞ

De Santis and Shams (2019) 45.08 145.54 295.93
Current numerical results 44.51 146.64 302.88
Error 1.26% 0.76% 2.35% X

R
less than 3%. Considering that the load parameters and mesh were
different from those in the literature, such an error is considered
acceptable. Therefore, the fluid-structure coupling algorithm
adopted in this study is reliable.

3.3. Grid independence verification Fig. 3.6. Position of monitoring point R.

Reasonable grid division is very important for the numerical

computation. Too few cells will lead to low computational accu-

racy, while too many cells will waste computational resources. (L/Vz) was monitored, and the time-averaged axial velocity V z
Grid independence verification was carried out for the examples within a period T was obtained. The calculated results are shown
used in this study. in Table 3.3. The grid-independence verification includes the axial
To facilitate the analysis of the influence of grid parameters on grid independence verification and XY-plane grid independence
the flow field, the monitoring point R (
P/2, 0, 80D), as shown in verification. Different computing grids were generated by control-
Fig. 3.6, was selected. After the flow field stabilized, the change ling the number of grid nodes in different directions. First, the
in the axial velocity Vz of point R in the calculation period T number of grid nodes in the XY-plane was controlled to 70  70,
8
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Table 3.3 while the number of grid nodes in the Z -direction was increased
Variation of time-averaged axial velocity with grid nodes. from 300 to 450. The results show that the difference between

Mesh 03 and Mesh 04 is negligible. To minimize the calculation
Mesh XYZ V z (m/s)
cost while maintaining the accuracy, the number of grid nodes in
Mesh 01 70*70*300 2.28471
the Z-direction was maintained at 400 in the following calculation,
Mesh 02 70*70*350 3.06085
Mesh 03 70*70*400 3.13747 and the independence of the XY-plane grid was further verified.
Mesh 04 70*70*450 3.11841 When the number of grid nodes in the XY-direction increased from
Mesh 05 60*60*400 3.18038 60  60 to 70  70, the time-averaged axial velocity decreased by
Mesh 06 90*90*400 3.12321
1.35% relative to Mesh 05. When it continued to increase to
Mesh 07 100*100*400 3.10687
90  90, the time-averaged axial velocity decreased by 0.45%
(<1%) relative to Mesh 03. Therefore, Mesh 03 was selected as
the final grid scheme.

4. Results and discussion

Y
LX LZ3 LZ1
4.1. Effect of parameter P/D on the flow field

The discussion in this section is divided into three parts: the

influence of the P/D ratio on the formation of the vortex structure,
X
the influence of the P/D ratio on the heat transfer efficiency of the
fluid in the channel, and the phenomenon of fluidic vibration of the
circular tube in the channel. The geometric model and numerical
calculation method were described in Section 2. The effects of
LZ2 three operating conditions: P/D = 1.050, 1.075, and 1.100 on the
above phenomena were analyzed separately. P/D < 1.050 was not
LY studied because the gap between fuel rods is too small under this
condition. Deformation may cause the fuel rods to come into con-
Fig. 4.1. Schematic diagram of the reference lines. tact with each other, resulting in fuel rod failure. Therefore, the
extremely narrow case is no longer a practical geometry for

0.20
0.15 Line LZ1
Line LZ2
0.10
Line LZ3
Vx/Ubulk

0.05
0.00
-0.05
-0.10
-0.15
-0.20
0 12 24 36 48 60 72 84 96 108
Z/D

(a)

0.20
0.15 Line LZ1
Line LZ2
0.10 Line LZ3
Vy/Ubulk

0.05
0.00
-0.05
-0.10
-0.15
-0.20
0 12 24 36 48 60 72 84 96 108
Z/D

(b)
Fig. 4.2. Distributions of the (a) horizontal velocity V x and (b) lateral velocity V y along the axis direction at P/D = 1.050.

9
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

0.0012
Line LZ1
0.0008 Line LZ2
Line LZ3
0.0004
Vx/Ubulk

0.0000

-0.0004

-0.0008

-0.0012
0 12 24 36 48 60 72 84 96 108
Z/D

(a)

0.0012
Line LZ1
0.0008 Line LZ2
Line LZ3
0.0004
Vy/Ubulk

0.0000

-0.0004

-0.0008

-0.0012
0 12 24 36 48 60 72 84 96 108
Z/D

(b)
Fig. 4.3. Distributions of the (a) horizontal velocity V x and (b) lateral velocity V y along the axis direction at P/D = 1.075.

nuclear fuel applications (Lee et al., 2013). When P/D  1.100, the pared with the fluctuations of the horizontal velocity on line LZ1,
vortex-street phenomenon disappears, which is also not discussed the amplitude decreased by 93% because the position of line LZ3
in this paper. is far away from the gap region, resulting in a weakened tendency
For the convenience of the following discussion, we first define of the fluid to move to the sides. The horizontal and lateral veloc-
the reference lines to be used, as shown in Fig. 4.1. We define line ities do not fluctuate significantly on LZ2 and LZ1 because the fluid
(0,P/2,Z) as LZ1, line (
P/2,0,Z) as LZ2, line (
P/2, P/2,Z) as LZ3, line is blocked by the circular pipe wall and it is difficult to flow to both
(X,P/2,81D) as LX, and line (
P/2,Y, 81D) as LY. sides in the corresponding directions.
Since all the calculations in this paper are transient, all the When P/D = 1.075, the horizontal and lateral velocities also fluc-
results are selected after the flow field reaches stability. It is veri- tuated on LZ1 and LZ2, respectively. However, the fluctuations
fied that the flow field can be stabilized after 8 cycles of calculation started with a lag compared with the condition of P/D = 1.050.
(the cycles mentioned here are in terms of channel length/inlet The fluctuations occurred only when Z/D > 60, with the maximum
velocity). Therefore, all subsequent results are selected after 8 fluctuation amplitude of 0.0009 U bulk , which is only 0.5% of the
computational cycles. average fluctuation amplitude under the condition of P/
Figs. 4.2–4.4 show the distribution curves of the horizontal D = 1.050. This indicates that when P/D increased to 1.075, the ten-
velocity V x and lateral velocity V y along lines LZ1–LZ3 at the work- dency of the fluid to flow to both sides in the slit zone diminishes
ing conditions P/D = 1.050–1.100. In the case of P/D = 1.050, when significantly and only very weak fluctuations exist. When P/D con-
Z/D > 48, regular fluctuations of the horizontal velocity appeared tinued to increase to 1.100, it can be seen from Fig. 4.4 that the
on line LZ1, with an average amplitude of 0.159 U bulk and average fluctuations of horizontal and lateral velocities on lines LZ1–LZ3
wavelength of 2.83D. The wavelength mentioned here is the axial completely disappeared. There was a small fluctuation only at
distance between two consecutive peaks. Regular fluctuations of the exit position, which indicates that the tendency of the fluid
the lateral velocity also appeared on line LZ2, with an average to flow to both sides in the gap region can be neglected at this
amplitude of 0.175 U bulk and average wavelength of 2.83D. These operating condition.
indicate a strong tendency for the fluid in the interstitial zone to The axial velocity distribution on lines LZ1–LZ3 at different P/D
flow to both sides of the channel under current conditions. Owing values is shown in Fig. 4.5. The axial velocity decreased with the
to the geometric symmetry, the fluctuations in lateral and trans- decrease in P/D when Z/D < 60 on line LZ1. As shown in Fig. 4.5
verse velocities are very similar in magnitude and wavelength for (a), when P/D = 1.100, the axial velocity was stable at 0.73U bulk ;
LZ1 and LZ2. It should be noted that the horizontal and lateral when P/D decreased to 1.050, the stable axial velocity decreased
velocities showed obvious fluctuations on line LZ3, with an average to 0.53U bulk , which is 27.4% less than that at P/D = 1.100. This is
amplitude of 0.012 U bulk and average wavelength of 2.96D. Com- attributed to the fact that as the P/D value decreased, the gap

10
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

0.0012
Line LZ1
0.0008 Line LZ2
Line LZ3
0.0004
Vx/Ubulk

0.0000

-0.0004

-0.0008

-0.0012
0 12 24 36 48 60 72 84 96 108
Z/D

(a)

0.0012

0.0008 Line LZ1

Line LZ2
0.0004 Line LZ3
Vy/Ubulk

0.0000

-0.0004

-0.0008

-0.0012
0 12 24 36 48 60 72 84 96 108
Z/D

(b)
Fig. 4.4. Distributions of the (a) horizontal velocity V x and (b) lateral velocity V y along the axis direction at P/D = 1.100.

between the circular pipes narrowed, decreasing the axial flow as shown in Fig. 4.1. The value Z = 81D was selected for the analysis
area of the fluid and impeding the flow in the axial direction. When because a stable vortex-street structure formed at this position,
Z/D > 60, the axial velocity at P/D = 1.100 and 1.075 did not change. and was far away from the exit position to avoid the influence of
At P/D = 1.050, the axial velocity began to fluctuate and rise, and the exit boundary.
tended to be in a stable fluctuating state when Z/D > 72. This is Fig. 4.6 shows the distribution of the time-averaged axial veloc-
because a stable vortex-street was formed in the area where Z/ ity on line LX. It can be seen that the magnitude of is positively cor-
D > 60, which increased the axial flow velocity owing to the mixing related with the cross-sectional area of the channel. The larger the
effect of the vortex-street structure. Moreover, owing to the sym- channel cross-sectional area, the greater the axial velocity. The two
metry of the geometric structure, the phenomena on line LZ2 were maximum values of appeared on both sides of the cylinder, i.e.,
basically the same as those on line LZ1. because this position is the largest the circulation area. The mini-
The phenomenon on line LZ3 was the opposite of the results on mum value of appeared at the position x = 0, which is the position
the above two lines: when Z/D < 60, the axial velocity increased as with the smallest circulation area. It is worth noting that the time-
P/D decreased. As shown in Fig. 4.5(c), when P/D = 1.100, the axial averaged axial velocities obtained at x = 0 were basically the same
velocity was stable at 1.27U bulk ; when P/D decreased to 1.050, the for the two operating conditions of P/D = 1.050 and 1.100. This is
stable axial velocity increased to 1.30U bulk , which is 2.4% higher because of the appearance of the vortex-street structure at the
than that at P/D = 1.100. This is because with the decrease in P/ gap region position when P/D = 1.050, which increased the axial
D, the axial velocity in the gap region decreased, leading to the velocity. In other words, the appearance of the vortex-street struc-
diversion of fluid to both sides. However, the increase in velocity ture accelerates the axial flow velocity in the gap region.
was not significant because of the large area of the subchannel. Fig. 4.7 shows the distribution of the time-averaged lateral
Similarly, the presence of the vortex-street structure caused signif- velocity on line LX. When P/D = 1.050, the time-averaged lateral
icant fluctuations of the axial velocity in the region of Z/D > 60 at P/ velocity is stable at approximately 0 in the region where
0.1 < X/
D = 1.050. D < 0.1. This is because the circulation area is the smallest in this
The above analyses were conducted on the distribution of tran- region, and the fluid was impeded from moving in the Y-
sient values, but because the vortex structure is an unstable flow direction. In the region where 0.1 < X/D <
0.1, the time-
field, the time-averaged values were more responsive to the char- averaged lateral velocity first decreased to
0.0012, then increased
acteristics of the flow field. The following is an analysis of the dis- rapidly, reached the maximum value of 0.0018 at , decreased to
tribution of the time-averaged values of the axial velocity,
0.0005, and finally stabilized when it increased to approximately
horizontal velocity, and lateral velocity on line LX, which is defined 0. This could be explained by the initial stability of the fluid in the

11
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

1.2
P/D:1.050
P/D:1.075
1.0 P/D:1.100
Vz/Ubulk

0.8

0.6

0.4
0 12 24 36 48 60 72 84 96 108
Z/D

(a)

1.2
P/D:1.050
P/D:1.075
1.0 P/D:1.100
Vz/Ubulk

0.8

0.6

0.4
0 12 24 36 48 60 72 84 96 108
Z/D

(b)

1.4

1.3
Vz/Ubulk

1.2

P/D:1.050
1.1 P/D:1.075
P/D:1.100

1.0
0 12 24 36 48 60 72 84 96 108
Z/D

(c)
V
Fig. 4.5. Distributions of the axial velocity V z along lines LZ1–LZ3 for P/D = 1.050–1.100: (a) the axial velocity on line LZ1, (b) the axial velocity on line LZ2, and (c) the axial
velocity on line LZ3.

gap region, followed by appearance of a pair of opposite vortices values of in the region of , i.e., the extreme value of 0.0054 at
that caused the fluid to diffuse to both sides. When P/D increased and the minimum value of
0.0054 at . The time-averaged hori-
to 1.075, the time-averaged lateral velocity remained stable at zontal velocity was greater than 0 in the region of , and less than
approximately 0 without significant fluctuations owing to the 0 in the region of . The main reason for this velocity distribution
absence of a vortex-street structure. is the presence of secondary flow in the subchannel when the value
Fig. 4.8 shows the distribution of the time-averaged horizontal of P/D is large. When P/D = 1.050, there were four extremes in the
velocity on line LX. When P/D  1.075, there were two extreme range of for the time-averaged horizontal velocity, two of which

12
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

1.6 disappeared and only very small fluctuations existed near the out-
P/D:1.050 let. When P/D further increased to 1.100, all the vortices
P/D:1.075 disappeared.
1.4 P/D:1.100 The above analysis shows that when P/D = 1.050, a stable
vortex-street structure appeared at the narrow-slit channel, caus-
ing the fluid in the gap region to diffuse to both sides. When P/D
Time-average(Vz)/Ubulk

1.2 increased to 1.075, the vortex-street structure disappeared.

4.2. Effect of parameter P/D on structural vibration

1.0
The appearance of the vortex-street structure not only intensi-
fied the turbulent mixing but also triggered fluid-induced vibra-
0.8
tions of the cylinder. This vibration differs from the general
turbulent pulsation-excited vibration, which is due to the motion
of macroscopic-scale vortices. It causes fluctuations in the pressure
0.6
and velocity of the flow field around the cylinder, thus inducing the
vibration of the cylinder. The computational setup of the flow field
0.4 is consistent with that described above, and the RSM model was
-1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 adopted. To facilitate the analysis of structural displacements, the
X/D cylinder was divided into ten equal parts, and nine points were
taken to monitor the displacements; the specific geometric nodes
Fig. 4.6. Distribution of the time-averaged axial velocity on line LX. are schematically shown in Fig. 4.10.
The central point of the cylinder, i.e., point-5, was taken as the
monitoring point. The displacements of the node over time under
0.0025 different working conditions were given, as shown in Figs. 4.11–
P/D:1.050
4.13. Fig. 4.11 shows the displacement curve of the midpoint of
0.0020
P/D:1.075 the cylinder in different directions with time under the working
P/D:1.100 condition P/D = 1.050. Owing to the presence of the vortex-street
0.0015
structure, the vortex motion caused the fluid to transfer the
Time-average(Vy)/Ubulk

0.0010 momentum in the Z-direction to the X- and Y-directions, resulting

in a stable micro-amplitude vibration of the cylinder under the
0.0005 current working condition. When P/D increased to 1.075, as shown
in Fig. 4.12, the vibration displacement of point-5 showed an obvi-
0.0000 ous attenuation trend; that is, a stable vibration could not be main-
tained because the stable vortex-street structure did not appear
-0.0005
under the current working condition. Despite some fluctuations
in the flow field that induced the vibration of the structure, the
-0.0010
vibration quickly decayed to the equilibrium position owing to
-0.0015 the damping effect of the axial fluid. When P/D continued to
increase to 1.100, the fluctuation phenomenon disappeared; hence,
-0.0020 the vibration of the structure was very weak and quickly attenu-
-1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 ated to the equilibrium position.
X/D To better analyze the influence of P/D on the vibration of the
cylindrical structure, the displacements of the nine nodes shown
Fig. 4.7. Distribution of the time-averaged lateral velocity on line LX.
in Fig. 4.10 were monitored, and the RMS values were calculated.
The results are shown in Fig. 4.14. When P/D = 1.050, the RMS val-
ues of the displacements of all the nodes in the X- and Y-directions
occurred at the same position as those at P/D = 1.075. The other were grater than 0. The maximum value appeared at Z = 0.6 L, and
two extreme values appeared in the region of , which indicate that the maximum RMS values of the displacement in the Y- and X-
in the region of , there is a tendency for the fluid to flow to both directions were 0.5348 lm and 0.4954 lm, respectively. In the
sides due to the presence of the vortex structure; the movement case of P/D = 1.075, the maximum RMS value of the displacement
of the vortex causes the fluid to move to both sides. In the outer occurred at Z = 0.5 L. The maximum RMS value of the displacement
region of the cylinder, the horizontal velocity distribution was in the Y-direction was 0.0855 lm, which decreased by 84.01%
the same for all three operating conditions. compared with that at P/D = 1.050. The maximum RMS value of
To show the influence of the P/D value on the flow field more displacement in the X-direction was 0.0513 lm, which decreased
intuitively, the contour clouds of the axial velocity are given below. by 89.64% compared with that at P/D = 1.050. When P/D further
The plane where Y = P/2 was selected, and the results are shown in increased to 1.100, the maximum RMS value of the displacement
Fig. 4.9. When P/D = 1.050, a regular vortex-street phenomenon also appeared at Z = 0.5 L. When P/D is equal to 1.050, the trigger
appeared in the region of Z > 60D. In the region of for the maximum displacement is located downstream due to the

0:5D < X < 0:5D, a pair of vortices with opposite rotational direc- presence of the vortex structure. When P/D is>1.075, the vortex
tions moved along the axial direction, and the axial velocity on structure disappears and the main trigger point is not the vortex,
both sides was significantly greater than that of the one in the mid- so the maximum displacement is the antinode of the first mode.
dle. The presence of such a velocity gradient leads to the formation The maximum RMS value of the displacement in the Y-direction
of a vortex-street structure. At the same time, owing to the move- was 0.0009 lm, which was 99.82% less than that at P/D = 1.050.
ment of the vortex, the fluid in the middle spread to both sides. The maximum RMS value of the displacement in the X-direction
When P/D increased to 1.075, the regular vortex-street structure was 0.0454 lm, which was 91.51% less than that at of P/
13
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

al velocity distribution was the same for all three operating conditio
0.008
P/D:1.050
0.006 P/D:1.075
P/D:1.100

0.004

Time-average(Vx)/Ubulk
0.002

0.000

-0.002

-0.004

-0.006

-0.008
-1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25
X/D

Fig. 4.8. Distribution of the time-averaged horizontal velocity on line LX.

(a)

(b)

(c)

Fig. 4.9. Contour clouds of the axial velocity: (a) P/D = 1.050, (b) P/D = 1.075, and (c) P/D = 1.100.

D = 1.050. It can be seen that as P/D increased, the vortex-street vortex-street structure can induce stable micro-amplitude vibra-
structure disappeared, leading to the disappearance of the vibra- tions. This result is in agreement with the conclusion of the study
tion of the cylindrical structure. Therefore, we believe that the by Dolfen et al. (Dolfen et al., 2019).

14
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

p p g y g
Fixed

Fluid-solid interface

Fixed
Point-10

Point-8 Point-9
Point-6 Point-7
Point-5
Point-4
Point-3
Point-0 Point-2
Point-1

Fig. 4.10. Schematic diagram of the monitored nodes.

2.0

1.5

1.0
X--Dissplacement( m)

0.5

0.0

-0.5

-1.0

-1.5

-2.0
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
T(s)

(a)

1.5

1.0
ement( m)

0.5

0.0
place

-0.5
Disp
Y-D

-1.0

-1.5

-2.0
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
T(s)

(b)

Fig. 4.11. Vibration displacement curve of the cylinder midpoint at P/D = 1.050: (a) displacement in the X-direction and (b) displacement in the Y-direction.

15
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Fig. 4.12. Vibration displacement curve of the cylinder midpoint at P/D = 1.075: (a) displacement in the X-direction and (b) displacement in the Y-direction.

4.3. Effect of parameter P/D on heat exchange that at the position of Z = 84D, the temperature for P/D = 1.050
was the lowest because when P/D decreased to 1.050, a signifi-
The vortex-street structure not only causes structural vibrations cant vortex-street structure appeared in the region of Z > 60D.
but also reinforces the convective heat exchange of the fluid in the Due to the mixing effect of the vortex, the fluid in the middle
channel owing to the mixing effect triggered by the vortex move- of the gap region has a tendency to flow to both sides, which is
ment. The effect of P/D on the fluid heat transfer in the channel is beneficial to the fluid heat transfer, leading to the decrease in
analyzed as follows. In Fig. 4.15, the distribution curves of the tem- the fluid temperature. Owing to the symmetrical effect of the
perature T on lines LZ1–LZ3 are presented for P/D = 1.050–1.100, geometry, the temperature distribution on line LZ2, shown in
where the initial temperatureT Initial = 310 K. Fig. 4.15(a) shows Fig. 4.15(b), is consistent with the results in Fig. 4.15(a); no
the distribution of temperature T on line LZ1. In the region of detailed analysis was carried out.
Z < 60D, the temperature on the line increased with the decrease Fig. 4.15(c) shows the distribution of temperature T on line LZ3.
in P/D. This is attributed to the fact that the distance between It can be seen that when P/D = 1.050, the temperature did not show
two cylinders in the channel decreased with the decrease in P/D. a downward trend in the area of Z > 60D, but continued to rise
According to the previous analysis, the fluid on both sides tended slowly along line LZ3. This is because line LZ3 is located in the sub-
to move toward the middle, which is not conducive to the heat channel region with a large fluid flow area; hence, the temperature
exchange between fluids, thus causing the temperature to rise. In variation trends of the three working conditions were the same,
the region of Z > 60D, the temperature of the fluid on the line con- and the difference between temperature values was small.
tinued to increase as P/D decreased from 1.100 to 1.075. However, The presence of the vortex-street structure caused fluctuations
when P/D further decreased to 1.050, the temperature of the fluid at transient values of temperature T. To better analyze the temper-
on the line showed a significant fluctuating downward trend and ature variations under different working conditions, the distribu-
stabilized in the region after Z > 72D. tion of the time-averaged values of temperature T on line LX is
To facilitate the analysis, the position Z = 84D was chosen as shown in Fig. 4.16. The temperature in the three working condi-
the reference point because the temperature for P/D = 1.050 tions show a distribution of high in the middle and low on both
tended to be stable at this point, and the point is far away from sides. This is because the circulation area in the middle position
the outlet position to avoid the boundary effect. It can be seen was the smallest and not conducive to the heat dissipation of the

16
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Fig. 4.13. Vibration displacement curve of the cylinder midpoint at P/D = 1.100: (a) displacement in the X-direction and (b) displacement in the Y-direction.

fluid. Away from the gap region, the circulation area gradually high temperature at the narrow slit because of the large circulation
increased, and the fluid temperature gradually decreased. In the area in the gap region. However, the temperature was lower in the
region of
0:25D < Z < 0:25D, there was a maximum value of the subchannel, and the temperature distribution was not uniform.
time-averaged temperature under the three working conditions. When P/D decreased to 1.075, owing to the decreased flow area
When P/D = 1.100, the maximum value of the time-averaged tem- in the gap region, an obvious high-temperature zone appeared in
perature was 1.047T Initial ; when P/D decreased to 1.075, the maxi- the section after Z > 0.42 m, and the same problem of lower tem-
mum value of the time-averaged temperature increased to perature and uneven temperature distribution in the subchannel
1.065T Initial . However, when P/D continued to decrease to 1.050, existed. When P/D further decreased to 1.050, the high-
the maximum value of the time-averaged temperature decreased temperature zone appeared earlier in the cross-section at
to 1.039T Initial . The reason for this phenomenon is consistent with Z = 0.14 m because the circulation area in the gap region was too
the above analysis. At the same time, it is worth noting that in this small, and the high-temperature zone became more obvious with
region, when P/D decreased to 1.050, the time-averaged tempera- the increase in Z. However, the high-temperature zone disap-
ture distribution became more uniform and the temperature gradi- peared at the section after Z > 0.56 m, and the temperature at
ent decreased significantly. This indicates that the emergence of the subchannel position increased owing to the appearance of
the vortex-street structure completely mixed the fluids, enhanced the vortex-street structure, which accelerated the fluid mixing
the heat transfer between fluids, and reduced the temperature gra- and caused the heat in the gap region to diffuse to both sides, mak-
dient. In the region outside of
0:25D < Z < 0:25D, the tempera- ing the temperature distribution more uniform.
ture distribution was essentially constant when P/D decreased To quantitatively analyze the influence of the vortex-street
from 1.100 to 1.075. However, when P/D decreased to 1.050, the structure on the heat transfer effect, the dimensionless parame-
temperature increased significantly, which indicates that the heat ter, the Nusselt number (Liu et al., 2012), was introduced in this
in the gap region was transported to both sides owing to the mix- study. The Nusselt number is the dimensionless temperature gra-
ing effect of the vortex-street structure, which caused the temper- dient on the heating wall that represents the convective heat
ature on both sides to increase significantly. transfer ability. Take the temperature of a section on the wall
Fig. 4.17 presents the time-averaged temperature distribution and calculate the average value to gettw . Similarly take the tem-
clouds on a total of six cross-sections from Z = 0 m to 0.7 m, respec- perature value of the fluid domain at the section and calculate
tively. It can be seen that when P/D = 1.100, there was no obvious the average value to gett b .

17
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Fig. 4.14. RMS values of displacements of all nodes: (a) P/D = 1.050 (b) P/D = 1.075, and (c) P/D = 1.100.

18
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

1.12
P/D:1.050
1.10 P/D:1.075
P/D:1.100
1.08
T/TInitial

1.06

1.04

1.02

1.00
0 12 24 36 48 60 72 84 96 108
Z/D

(a)

1.12
P/D:1.050
1.10 P/D:1.075
P/D:1.100
1.08
T/TInitial

1.06

1.04

1.02

1.00
0 12 24 36 48 60 72 84 96 108
Z/D

(b)
1.05
P/D:1.050
1.04 P/D:1.075
P/D:1.100
T/TInitial

1.03

1.02

1.01

1.00
0 12 24 36 48 60 72 84 96 108
Z/D

(c)
Fig. 4.15. Temperature distributions on lines LZ1–LZ3 for P/D = 1.050–1.100: (a) temperature on line LZ1, (b) temperature on line LZ2, and (c) temperature on line LZ3.

hz D h qw
Nu ¼ ; hz ¼ ð3-1Þ The Nusselt number distribution in the axial direction of the
K tw
tb
channel is shown in Fig. 4.18. When P/D = 1.100, the Nusselt num-
where: ber gradually decreased with the increase in Z. This is because of
the large temperature difference between the wall surface and
hz —local convective heat transfer coefficient, W/(m2K); the fluid at the inlet, which is favorable to heat transfer from the
Dh —hydraulic diameter, m; wall surface to the fluid. As Z increased, the fluid temperature grad-
K—heat conductivity of fluids, W/(mK); ually increased, and the temperature difference between the wall
qw —wall heat flux, W/m2; surface and fluid gradually decreased, leading to the gradual
tw —average wall temperature, K; decrease in the Nusselt number. When P/D decreased to 1.075,
tb —temperature of the main region, and the average of the sec- the trend of the Nusselt number was basically the same as that
tion is taken, K. at P/D = 1.100; however, the overall Nusselt number was smaller

19
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

1.08 than that at P/D = 1.100. When P/D was further reduced to 1.050,
P/D:1.050 the Nusselt number first decreased in the region of Z < 60D, then
1.07 P/D:1.075 increased in the region of Z > 60D, and stabilized near an equilib-
P/D:1.100 rium value. Taking the position of Z = 84D as a reference, the Nus-
1.06 selt number was 90.76 when P/D = 1.100, and decreased to 84.86
when P/D decreased to 1.075, which is a 6.51% decrease relative
Time-average(T)/TInitial

to that at P/D = 1.100. When P/D was further reduced to 1.050,

1.05
the Nusselt number increased to 94.33, which is a 3.93% increase
compared with that at P/D = 1.100. It can be observed that the
1.04
presence of the vortex-street structure can significantly enhance
the heat transfer in the channel.
1.03

1.02
5. Conclusions
1.01
In this study, the RSM-e model was used to simulate a typical
four-subchannel flow field containing a single cylinder. The forma-
1.00
-1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 tion conditions of the vortex-street structure were investigated,
X/D and the fluid-induced vibration phenomenon of the cylinder in
the channel was studied using a weak coupling algorithm. Finally,
Fig. 4.16. Distribution of time-averaged temperature on line LX. the influence of the vortex-street structure on the heat transfer

Z=0.70
Z=0.56
Z=0.42
Z=0.28

Z=0.14

Z=0.00

(a)

Z=0.70
Z=0.56
Z=0.42
Z=0.28

Z=0.14

Z=0.00

(b)
Fig. 4.17. Time-averaged temperature clouds on each section. (a) P/D = 1.100, (b) P/D = 1.075, and (c) P/D = 1.050.

20
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Z=0.70
Z=0.56
Z=0.42
Z=0.28

Z=0.14

Z=0.00

(c)

Fig. 4.17 (continued)

150 When P/D = 1.050, the vortex-street structure inside the chan-
nel caused a stable micro-amplitude vibration of the cylinder.
P/D:1.050
140 When P/D increased, the vortex-street structure disappeared
P/D:1.075
P/D:1.100 along with the vibration of the cylinder.
130 4. The vortex-street structure intensified the heat transfer in the
flow field, resulting in a significant increase of the Nusselt num-
120 ber in the channel. This is because of the fact that the presence
of the vortex-street structure accelerated the mixing of the flu-
110 ids and caused the heat at the narrow slit to diffuse to both
Nu

sides, which strengthened the heat transfer of the fluids in the

100 channel and made the temperature distribution more uniform.

90 CRediT authorship contribution statement

80 Yonghui Guo: Conceptualization, Methodology, Software, Writ-

ing – review & editing. Xiaochang Li: Data curation, Writing – orig-
70 inal draft, Funding acquisition. Ju Liu: Formal analysis. Guangliang
Chen: Investigation, Funding acquisition. Ye Gao: Supervision.
60
0 12 24 36 48 60 72 84 96 108
Declaration of Competing Interest
Z/D
The authors declare that they have no known competing finan-
Fig. 4.18. Nusselt number distribution in the axial direction of the channel.
cial interests or personal relationships that could have appeared
to influence the work reported in this paper.
efficiency of the fluid in the channel was discussed. The conclu-
sions are as follows: Acknowledgements

1. When P/D is very small (1.050 in this study), a vortex-street The authors are grateful for the support of the National Natural
structure appears in the gap region of the channel. This is Science Foundation of China (Project No. 11805044, China), Science
because the velocity of the fluid in the gap region is less than and Technology on Reactor System Design Technology Laboratory
that of the fluid in the channels on both sides, and the shear funded project (HT-KFKT-02-2020007, China), National Natural
stress generated by the velocity gradient causes the vortex- Science Foundation of China (Project No. 51909045, China).
street structure to appear. When P/D increased to a certain
value (1.075 in this study), the vortex-street structure in the References
gap region of the channel disappeared.
2. The vortex-street structure is a pair of vortices that rotate in Baeten, P., Schyns, M., Fernandez, R., De Bruyn, D., Van den Eynde, G., 2014.
MYRRHA: A multipurpose nuclear research facility. EPJ Web of Conferences 79,
opposite directions and move along the axial direction. The 03001.
presence of vortex motion in the gap region led to a tendency Chang, D., Tavoularis, S., 2005. Unsteady Numerical Simulations of Turbulence and
for the fluid in the gap region to move to both sides. Coherent Structures in Axial Flow Near a Narrow Gap. J. Fluids Eng. 127 (3),
458–466.
3. The presence of the vortex-street structure caused a steady Chang, D., Tavoularis, S., 2007. Numerical simulation of turbulent flow in a 37-rod
micro-amplitude vibration of the cylinder inside the channel. bundle. Nucl. Eng. Des. 237 (6), 575–590.

21
Y. Guo, X. Li, J. Liu et al. Annals of Nuclear Energy 171 (2022) 109055

Chang, D., Tavoularis, S., 2012. Numerical simulations of developing flow and vortex Holloway, M.V., Beasley, D.E., Conner, M.E., 2008. Single-phase convective heat
street in a rectangular channel with a cylindrical core. Nucl. Eng. Des. 243, 176–199. transfer in rod bundles. Nucl. Eng. Des. 238 (4), 848–858.
Chun, T.H., Shin, C.W., In, W.K., Lee, K.H., Bae, K.H., Song, K.W., 2009. A Potential of Horváth, Á., Dressel, B., 2012. Numerical simulations of square arrayed rod bundles.
Dual Cooled Annular Fuel for OPR-1000 Power up-rate. In: Proceedings of the Nucl. Eng. Des. 247, 168–182.
Water Reactor Fuel Performance Meeting – WRFPM / Top Fuel 2009, France. p. In, W.-K., Shin, C.-H., Oh, D.-S., Chun, T.-H., 2004. Numerical analysis of the
268. turbulent flow and heat transfer in a heated rod bundle. Nucl. Eng. Technol. 36
De Pauw, B., Lamberti, A., Ertveldt, J., Rezayat, A., van Tichelen, K., Vanlanduit, S., (2), 153–164.
Berghmans, F., 2016. Vibration Monitoring Using Fiber Optic Sensors in a Lead- Lee, C.Y., Shin, C.H., In, W.K., 2013. Effect of gap width on turbulent mixing of
Bismuth Eutectic Cooled Nuclear Fuel Assembly. Sensors 16 (4), 571. https:// parallel flow in a square channel with a cylindrical rod. Exp. Therm Fluid Sci. 47,
doi.org/10.3390/s16040571. 98–107.
De Ridder, J., Degroote, J., Van Tichelen, K., Schuurmans, P., Vierendeels, J., 2013. Li, X., Gao, Y.e., 2014. Methods of simulating large-scale rod bundle and application to
Modal characteristics of a flexible cylinder in turbulent axial flow from a 1717 fuel assembly with mixing vane spacer grid. Nucl. Eng. Des. 267, 10–22.
numerical simulations. J. Fluids Struct. 43, 110–123. Liu, C.C., Ferng, Y.M., Shih, C.K., 2012. CFD evaluation of turbulence models for flow
De Ridder, J., Van Tichelen, K., Degroote, J., Vierendeels, J., 2016. Vortex-induced simulation of the fuel rod bundle with a spacer assembly. Appl. Therm. Eng. 40,
vibrations by axial flow in a bundle of cylinders. In: 11th International 389–396.
Conference on Flow-Induced Vibration, pp. 1–8. Merzari, E., Ninokata, H., Baglietto, E., 2008. Numerical simulation of flows in tight-
De Santis, D., Shams, A., 2017. Numerical modeling of flow induced vibration of lattice fuel bundles. Nucl. Eng. Design 238 (7), 1703–1719.
nuclear fuel rods. Nucl. Eng. Des. 320, 44–56. Meyer, L., 2010. From discovery to recognition of periodic large scale vortices in rod
De Santis, D., Shams, A., 2019. Analysis of flow induced vibrations and static bundles as source of natural mixing between subchannels—a review. Nucl. Eng.
deformations of fuel rods considering the effects of wire spacers and working Des. 240 (6), 1575–1588.
fluids. J. Fluids Struct. 84, 440–465. Meyer, L., Rehme, K.J.E.T., Science, F., 1994. Large-scale turbulence phenomena in
Dolfen, H., Bertocchi, F., Rohde, M., Degroote, J., 2019. Vibrations in a 7-rod bundle compound rectangular channels, vol. 8, no. 4, pp. 286–304.
subject to axial flow: Simulations and experiments. Nucl. Eng. Des. 353, 110227. Möller, S.V., 1991. On phenomena of turbulent flow through rod bundles. Exp.
https://doi.org/10.1016/j.nucengdes.2019.110227. Therm Fluid Sci. 4 (1), 25–35.
El-Genk, M.S., Su, B., Guo, Z., 1993. Experimental studies of forced, combined and Möller, S.V., 1992. Single-phase turbulent mixing in rod bundles. Exp. Therm Fluid
natural convection of water in vertical nine-rod bundles with a square lattice. Sci. 5 (1), 26–33.
Int. J. Heat Mass Transf. 36 (9), 2359–2374. Païdoussis, M.P., Gagnon, J.O., 1984. Experiments on vibration of clusters of
Galiana Gonzalez, B., 2010. Material requirements for a thorium based nuclear fuel. cylinders in axial flow: Modal and spectral characteristics. J. Sound Vib. 96 (3),
Guellouz, M.S., Tavoularis, S., 2000a. The structure of turbulent flow in a rectangular 341–352.
channel containing a cylindrical rod – Part 1: Reynolds-averaged Shin, C.-H., Chun, T.-H., Oh, D.-S., In, W.-K., 2012. Thermal hydraulic performance
measurements. Exp. Therm Fluid Sci. 23 (1-2), 59–73. assessment of dual-cooled annular nuclear fuel for OPR-1000. Nucl. Eng. Des.
Guellouz, M.S., Tavoularis, S., 2000b. The structure of turbulent flow in a rectangular 243, 291–300.
channel containing a cylindrical rod - Part 2: phase-averaged measurements. Tapucu, A., Merilo, M., 1977. Studies on diversion cross-flow between two parallel
Exp. Therm Fluid Sci. 23 (1-2), 75–91. channels communicating by a lateral slot. II: Axial pressure variations. Nucl.
Házi, G., 2005. On turbulence models for rod bundle flow computations. Ann. Nucl. Eng. Des. 42 (2), 307–318.
Energy 32 (7), 755–761. ter Hofstede, E., Kottapalli, S., Shams, A., 2017. Numerical prediction of flow induced
Holloway, M.V., McClusky, H.L., Beasley, D.E., Conner, M.E., 2004. The Effect of vibrations in nuclear reactor applications. Nucl. Eng. Des. 319, 81–90.
Support Grid Features on Local, Single-Phase Heat Transfer Measurements in Yao, S.C., Hochreiter, L.E., Leech, W.J., 1982. Heat-Transfer Augmentation in Rod
Rod Bundles. J. Heat Transfer 126 (1), 43–53. Bundles Near Grid Spacers. J. Heat Transfer 104 (1), 76–81.
Holloway, M.V., Conover, T.A., McClusky, H.L., Beasley, D.E., Conner, M.E., 2005. The
Effect of Support Grid Design on Azimuthal Variation in Heat Transfer
Coefficient for Rod Bundles. J. Heat Transfer 127 (6), 598–605.

22