# Seismic response of cylinder assemblies in axial flow

Roberto Capanna, Guillaume Ricciardi, Emmanuelle Sarrouy, Christophe Eloy

► To cite this version:

Roberto Capanna, Guillaume Ricciardi, Emmanuelle Sarrouy, Christophe Eloy. Seismic response of cylinder assemblies in axial flow. Journal of Fluid Mechanics, 2022, 943, 10.1017/jfm.2022.433. hal-03689213

HAL Id: hal-03689213

https://hal.science/hal-03689213

Submitted on 7 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Seismic response of cylinder assemblies in axial flow

Roberto Capanna¹,†, Guillaume Ricciardi¹, Emmanuelle Sarrouy² and Christophe Eloy³

¹CEA, DES, IRESNE, Department of Nuclear Technology, Cadarache, 13108 Saint-Paul-Lez-Durance, France

²Aix Marseille Univ, CNRS, Centrale Marseille, LMA UMR 7031, 13453, Marseille, France

³Aix Marseille Univ, CNRS, Centrale Marseille, IRPHE, 13453, Marseille, France

(Received 29 October 2021; revised 14 March 2022; accepted 27 April 2022)

Earthquakes are a great challenge for the safety of nuclear reactors. To address this challenge, we need to better understand how the reactor core responds to seismic forcing. The reactor core is made of fuel assemblies, which are themselves composed of flexible fuel rods immersed in a strong axial flow. This gives rise to strongly coupled fluid-structure interactions whose accurate modelling generally requires high computational costs. In this paper, we introduce a new model able to capture the mechanical response of the reactor core subjected to seismic forcing with low computational costs. This model is based on potential flow theory for the fluid part, and Euler-Bernoulli beam theory for the structural part, allowing us to predict the response to seismic forcing in presence of axial flow. The linear equations are solved in the Fourier space to decrease computational time. For validation purposes, first we use the proposed model to compute the response of a single cylinder in axial flow. We then implement a multiple-cylinder geometry made of four fuel assemblies, each made of 8 × 8 cylinders, corresponding to an experimental facility available at CEA. The comparison between numerical results and experiments shows good agreement. The model can predict correctly the added mass. It can also capture qualitatively the coupling between assemblies and the effect of confinement. This shows that a potential flow approach can give insight into the complex fluid-structure interactions within a nuclear reactor and, in particular, be used to predict the response to seismic forcing at low computational cost.

**Key words:** flow-structure interactions, slender-body theory

† Email address for correspondence: capanna@gwu.edu

© The Author(s), 2022. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
R. Capanna, G. Ricciardi, E. Sarrouy and C. Eloy

# 1. Introduction

One of the main concerns for the safety of nuclear power plants is represented by earthquakes. During an earthquake, the main risk is that fuel assemblies start to move and potentially touch each other or prevent the drop of the control rods used to cool the core. To better understand the motion of fuel assemblies subjected to seismic forcing, fluid-structure interaction models are needed. In this paper, we present such a model, which has been developed with the objective of gaining insight into the complex fluid-structure interactions at play inside a nuclear reactor, while remaining computationally efficient.

A reactor core of a pressurised water reactor (PWR) is made of fuel assemblies (between 150 and 250 depending on the power of the reactor). Each fuel assembly gathers about 100 fuel rods, stands between four and five metres high, has a cross-section of about 20 × 20 cm², and weighs about 800 kg. The fuel rods, which contain pellets of uranium, have a diameter of about 1 cm for 4 m in height; the space between each pair of fuel rods in the assembly is about 3 mm. Fuel rods are held together by grids (between 8 and 10 depending on the power of the reactor) distributed along the height of the fuel assembly. Springs and dimples are used in between the grid and the rods to avoid any drop of the fuel rods (De Mario & Street 1989).

The fuel assemblies are cooled by a strong axial flow. This water flow is upwards, with velocities of about 5 m s⁻¹ at 150 bar and 310 °C. The flow regime is fully turbulent, with a Reynolds number based on the rod diameter of $Re \approx 5 \times 10^5$. Note that even if the main flow is upwards, the root-mean-square average of the transverse component is between 5 % and 15 % of the vertical velocity.

The presence of the water flow gives rise to strongly coupled interactions between the fluid and the structure (Chen & Wambsganss 1972): the motion of the structure modifies the fluid flow, which itself exerts a force on the structure. A first attempt to describe fluid-structure interactions is to use the concepts of added mass and added damping. The added mass is the inertial mass added to a body because of the presence of a surrounding fluid. For simplicity, it can be viewed as a volume of fluid moving with the same velocity as the body, though in reality all fluid particles will be moving to various degrees (Buat 1779; Lamb 1895). The added mass $M_a$ of a non-deformable body moving at velocity $U$ through an unbounded fluid otherwise at rest can be defined such that the kinetic energy in the volume of fluid is $E_k = (M_a U^2)/2$, or equivalently,

$$M_a = \frac{\rho}{U^2} \int \|V\|^2 dv, \qquad (1.1)$$

where $V$ is the velocity field. A complete collection of added masses for different geometries and flow conditions can be found in Wendel (1956) or Brennen (1982), for instance.

While the added mass is due mainly to pressure forces exerted on the body, viscous forces and boundary layer separation give rise to drag and to an added damping effect. Taylor (1952) proposed a model for the damping force on a slender structure, which has been used widely (e.g. Païdoussis 1966a; Triantafyllou & Cheryssostomidis 1985; Gosselin & De Langre 2011; Singh, Michelin & De Langre 2012). In Taylor's model, the damping force is decomposed into two terms: a normal force akin to drag, and a longitudinal force depending on the body surface and on the incident angle.

Generally, the total fluid force exerted on the non-deformable body is assumed to be decomposed into its added mass component (proportional to the acceleration) and its drag component (proportional to the square of the velocity). This decomposition, known as the
Morison equation (Morison, Johnson & Schaaf 1950), can describe correctly experimental observations provided that drag and inertial coefficients are adjusted empirically (in general, they depend on the motion amplitude and frequency).

For deformable cylinders in axial flow, the concept of added mass does not apply directly. In that case, one can use slender body theory developed by Lighthill (1960a,b). This theory makes use of potential flow theory and arguments of momentum balance in slices of fluids along the slender body. The resulting normal force per unit length exerted by the fluid can be written in the limit of small displacement as

$$F(X, T) = -M_a (\partial_T + U \partial_X)^2 W, \qquad (1.2)$$

where $W(X, T)$ is the normal displacement, $U$ is the undisturbed axial flow velocity, and here, $M_a$ is the added mass per unit length ($M_a = \rho\pi A^2$ for a circular cylinder of radius $A$). In the case of a solid body, $W$ has no $X$ dependence and one recovers a force opposite to $M_a \partial_T^2 W$, as expected.

Based on the work of Lighthill (1960b), Païdoussis (1966a,b) studied theoretically and experimentally the dynamics of a flexible slender cylinder clamped at its leading edge, which is immersed in an axial flow. He found that above a critical flow velocity, a flutter instability appears. Later, Païdoussis (1973) and Païdoussis & Pettigrew (1979) extended this flutter stability analysis to confined geometries, by using the work of Classen (1972) and Chen & Wambsganss (1972) on added mass in confined geometries.

The stability analysis of a single cylinder in axial flow has been later generalised to clusters of cylinders (Païdoussis 1973, 1979; Païdoussis & Suss 1977; Païdoussis, Suss & Pustejovsky 1977). If each cylinder is free to move independently, then the fluid acts as a coupling medium and the motions of the cylinders are synchronised. Above a critical flow velocity, the coupled pinned-pinned cylinders lose stability by divergence (unstable mode of eigenfrequency zero). For higher flow velocities, the system may be subjected to several divergence and flutter instabilities simultaneously. In recent years, this work of Païdoussis has been extended (e.g. De Langre et al. 2007; Michelin & Smith 2009; Schouveiler & Eloy 2009) and studied numerically both for a single cylinder (De Ridder et al. 2013, 2015) and for a cluster of cylinders (De Ridder et al. 2017).

The flutter of a flexible cylinder in axial flow bears similarities with the flutter of an elastic plate often referred to as the flag instability. For a plate, two asymptotic limits can be studied: slender structures where Lighthill's slender body theory applies, and wide plates for which a two-dimensional approach is suited (Wu 2001). For intermediate aspect ratios, the pressure distribution on a flexible plate can be solved by projecting the problem in Fourier space (Guo & Païdoussis 2000; Eloy, Souillez & Schouveiler 2007; Eloy et al. 2008, 2010; Doaré, Mano & Bilbao Ludena 2011a; Doaré, Sauzade & Eloy 2011b). With this projection, the approaches of slender body approximation and large span approximation can be generalised to any aspect ratio. In this paper, we will follow a similar approach to describe the flow around assemblies of flexible cylinders found in a PWR.

A different approach was proposed by Ricciardi et al. (2009) using a porous medium method. It now becomes possible to model both the fluid and the structure dynamics of a whole core. Some local information is lost compared to a direct numerical simulation, such as the vibrations of individual rods, but fluid-mediated interactions between fuel assemblies can be modelled. This model shows good agreement with experimental results on the response of the whole core to external forcing, but computational cost remains high.

At present, the models describing the fluid-structure interactions within a reactor core can be divided into two families: (1) complex numerical models with high computational costs, which usually hinder the understanding of physical mechanisms; and (2) linear
R. Capanna, G. Ricciardi, E. Sarrouy and C. Eloy

<figure><img src="image_4.png" /><figcaption>Figure 1. Representation of a confined pinned-pinned cylinder deformed under axial flow in (a) real space and (b) Fourier space.</figcaption></figure>

models with low computational costs, which do not capture important mechanisms such as fluid-mediated interactions between assemblies. In this paper, we introduce a new model based on potential flow theory, with the objective of providing an accurate, computationally effective modelling of fluid-structure interactions within a reactor core.

The paper is organised as follows. In § 2, the model is presented and the mathematical methodology used to solve the problem is described. The model is then applied to a single cylinder in § 3, and compared to results of the literature for validation purposes. In § 4, it is applied to a multiple-cylinder geometry, replicating the geometry of an experimental facility. Numerical results are compared to experimental data and discussed in § 5. Finally, in § 6, some conclusions are drawn.

## 2. Potential flow model

In this section, we describe the proposed model where the fluid flow is treated as potential and the structure as an Euler-Bernoulli beam. The model is introduced for a single cylinder in axial flow, but as we shall see, it can be generalised to multiple cylinders.

### 2.1. Problem statement

We consider a pinned-pinned cylinder of length $L$ and radius $A$ (figure 1a). We will use two coordinate systems: a Cartesian system $(X, Y, Z)$ with its origin at the bottom of the cylinder, and a cylindrical system $(R, \theta, Z)$ with the same origin (figure 1b).

The cylinder is immersed in a uniform axial flow with velocity $U$ and bounded in the $X$- and $Y$-directions by rigid walls at distances $L_x$ and $L_y$ from the cylinder axis. Without loss of generality, we consider that the cylinder deflects in the $Y$-direction, and the deflection is called $W(Z, T)$ (figure 1a).
We use the linearised Euler–Bernoulli beam equation (Bauchau & Craig 2009) to
describe the dynamics of the cylinder deflection

$$
M_s \partial_T^2 W + B \partial_Z^4 W = F_Y, \qquad (2.1)
$$

where $M_s$ is the mass of the cylinder per unit length, $B = EI$ is the bending rigidity (E
being Young’s modulus, and $I = (\pi A^4)/4$ the second moment of area), and $F_Y$ is the
force per unit length that the fluid exerts on the cylinder in the $Y$-direction.

We assume that the forces exerted by the fluid on the structure originate mainly from
pressure difference. Hence the force per unit length can be expressed as

$$
F_Y(Z, T) = - \int_{0}^{2\pi} P(A, \theta, Z, T) \sin(\theta) A \, d\theta, \quad (2.2)
$$

where $P(R, \theta, Z, T)$ is the pressure field in cylindrical coordinates.

The problem is made dimensionless using the cylinder radius $A$, the flow velocity $U$, and
the fluid density $\rho$ as characteristic length, speed and density. Dimensionless quantities are
designated with lowercase letters such that, for instance,

$$
r = \frac{R}{A}, \quad w = \frac{W}{A}, \quad z = \frac{Z}{A}, \quad l_x = \frac{L_x}{A}, \quad l_y = \frac{L_y}{A}, \quad t = \frac{UT}{A}, \quad p = \frac{P}{\rho U^2}. \quad (2.3a-g)
$$

In dimensionless form, the linearised Euler–Bernoulli equation (2.1) becomes

$$
m \partial_t^2 w + b \partial_z^4 w = f_y, \quad (2.4)
$$

where

$$
m = \frac{M_s}{\rho A^2}, \quad b = \frac{B}{\rho U^2 A^4}, \quad f_y = \frac{F_Y}{\rho U^2 A}. \qquad (2.5a-c)
$$

The flow is assumed to be potential, inviscid and incompressible, such that the flow
velocity is given by

$$
V = Ue_z + \nabla \Phi, \tag{2.6}
$$

where $e_z$ is the unit vector along $z$, and $\Phi(X, Y, Z, T)$ is the perturbation potential.
In dimensionless units, this potential is $\phi = \Phi/(UA)$. To find the flow around the
moving cylinder, one has to solve a Laplace problem with Neumann boundary conditions
(linearised for small displacements):

$$
\Delta\phi = 0, \tag{2.7}
$$

$$
\left. \frac{\partial \phi}{\partial x} \right|_{|x|=l_x} = 0, \qquad (2.8)
$$

$$
\left. \frac{\partial \phi}{\partial y} \right|_{|y|=l_y} = 0, \qquad (2.9)
$$

$$
\left. \frac{\partial \phi}{\partial n} \right|_{r=1} = -(\partial_t + \partial_z)w(z, t) \sin \theta, \quad (2.10)
$$

where (2.8) and (2.9) come from the impermeability of the walls, and (2.10) comes from
the impermeability on the cylinder wall.
Using the linearised unsteady Bernoulli equation, the dimensionless pressure field can be linked to $\phi$ (Capanna 2018) by

$$p(x, y, z, t) = -(\partial_t + \partial_z)\phi. \qquad (2.11)$$

Using this relation and applying the operator $(\partial_t + \partial_z)$ to the system (2.7)–(2.10) yields

$$\Delta p = 0, \qquad (2.12)$$

$$\left. \frac{\partial p}{\partial x} \right|_{|x|=l_x} = 0, \qquad (2.13)$$

$$\left. \frac{\partial p}{\partial y} \right|_{|y|=l_y} = 0, \qquad (2.14)$$

$$\left. \frac{\partial p}{\partial n} \right|_{r=1} = (\partial_t + \partial_z)^2 w(z, t) \sin \theta. \qquad (2.15)$$

This shows that the pressure field is also a solution to a Laplace equation with Neumann boundary conditions. This is why the term $P/\rho$ is called the acceleration potential in aerofoil theory.

## 2.2. Problem in Fourier space

To solve the set of equations (2.12)–(2.15), we will use the method proposed by Doaré et al. (2011b). It consists in expressing the Laplace problem in Fourier space along $z$ as

$$\hat{p}(x, y, k, t) = \mathcal{F}[p(x, y, z, t)], \qquad (2.16)$$

where $\hat{p}$ is the Fourier transform of $p$, with $\mathcal{F}[\cdot]$ defined as

$$\mathcal{F}[f(z)] = \frac{1}{2\pi} \int_{-\infty}^{+\infty} f(z) e^{-ikz} dz = \hat{f}(k), \qquad (2.17)$$

together with the inverse Fourier transform

$$\mathcal{F}^{-1}[\hat{f}(k)] = \int_{-\infty}^{+\infty} \hat{f}(k) e^{ikz} dk = f(z). \qquad (2.18)$$

The convolution product along $z$, denoted $\star$, is also introduced

$$f \star g = \int_{-\infty}^{+\infty} f(\zeta) g(z - \zeta) d\zeta = 2\pi \mathcal{F}^{-1}[\hat{f}\hat{g}]. \qquad (2.19)$$

In Fourier space, the three-dimensional Laplace problem is transformed into a two-dimensional Helmholtz problem:

$$(\partial_x^2 + \partial_y^2)\hat{p} = k^2\hat{p}, \qquad (2.20)$$

$$\left. \frac{\partial \hat{p}}{\partial x} \right|_{|x|=l_x} = 0, \qquad (2.21)$$

$$\left. \frac{\partial \hat{p}}{\partial y} \right|_{|y|=l_y} = 0, \qquad (2.22)$$

$$\left. \frac{\partial \hat{p}}{\partial n} \right|_{r=1} = \hat{\gamma}(k, t) \sin \theta, \qquad (2.23)$$
where $\hat{\gamma}(k, t)$ is the Fourier transform of the impermeability boundary condition, such that

$$ \hat{\gamma} = \mathcal{F}[(\partial_t + \partial_z)^2 w]. \qquad (2.24) $$

Equation (2.2) relates the force per unit length exerted by the fluid on the cylinder and the pressure field. Using dimensionless quantities and taking the Fourier transform of (2.2), one obtains the equivalent in Fourier space:

$$ \hat{f}_y(k, t) = - \int_0^{2\pi} \hat{p}(r = 1, \theta, k, t) \sin \theta \, d\theta, \qquad (2.25) $$

where $\hat{f}_y$ is the Fourier transform of $f_y$. Because the pressure $\hat{p}$ is a solution of the linear Helmholtz problem (2.20)–(2.23), it is proportional to $\hat{\gamma}(k, t)$. Hence (2.25) can be written as

$$ \hat{f}_y = -\hat{\mu}(k) \hat{\gamma}(k, t). \qquad (2.26) $$

where the function $\hat{\mu}(k)$ depends on $l_x$ and $l_y$ through the boundary conditions (2.21) and (2.22).

Equation (2.4) governing beam dynamics then becomes in Fourier space

$$ m \partial_t^2 \hat{w} + k^4 b \hat{w} = \hat{\mu}(k) [-(\partial_t^2 \hat{w} + 2ik \partial_t \hat{w} - k^2 \hat{w})], \qquad (2.27) $$

showing that $\hat{\mu}$ plays the role of an added mass in Fourier space and that added damping $2ik\hat{\mu}$ and added stiffness $-k^2\hat{\mu}$ are proportional to this added mass.

The objective now will be to calculate the added mass $\hat{\mu}(k)$. In practice, to calculate it, we can assume $\hat{\gamma} = 1$, compute numerically the solution $\hat{p}$ of the Helmholtz problem (2.20)–(2.23) for different values of $k$, $l_x$ and $l_y$, and finally use (2.25) and (2.26) to compute $\hat{\mu}(k)$.

### 2.3. Modal decomposition and range of interest for k

Considering an external forcing at a given frequency $\Omega$, the displacement of the cylinder can be written in dimensionless form as

$$ w(z, t) = \sum_{j=1}^{\infty} q_j(t) \chi_j(z), \qquad (2.28) $$

where $\chi_j(z)$ are the beam eigenmodes, and $q_j(t) = \eta_j e^{i\omega t}$ are the generalised coordinates, with $\eta_j$ the modal amplitudes (which do not depend on time), and $\omega = \Omega A/U$ the dimensionless forcing frequency.

To determine the interesting range for $k$ for which $\hat{\mu}$ should be computed, let us examine the eigenmodes $\chi_j(z)$ and their Fourier transform $\hat{\chi}_j(k)$ for pinned-pinned beams. Eigenmodes are given by (Blevins 2015)

$$ \chi_j(z) = \sqrt{2} \sin(k_j z), \quad \text{with } k_j = j\pi/l, \qquad (2.29) $$

where $l = L/A$ is the dimensionless cylinder length. Using this normalisation, the scalar product of two eigenmodes is

$$ \langle \chi_i, \chi_j \rangle = \frac{1}{l} \int_0^l \chi_i(z) \chi_j(z) \, dz = \delta_{ij}. \qquad (2.30) $$

As figure 2(b) shows, $k$ values greater than $10k_j$ have contributions at least 100 times smaller than low $k$ values in Fourier space. Hence they can be truncated.
R. Capanna, G. Ricciardi, E. Sarrouy and C. Eloy

<figure><img src="image_5.png" /><figcaption>Figure 2. (a) Eigenmodes for the first three modes of a pinned-pinned beam, and (b) their Fourier transforms.</figcaption></figure>

Classical values for PWR beams are $A \in [10^{-3}, 10^{-2}]$ m and $L \in [1, 10]$ m. Hence classical slenderness ratios $l$ belong to $[10^2, 10^4]$, and dimensionless wavenumbers are $k_j = j\pi/l \in [3 \times 10^{-4}, 0.3]$ when considering the first 10 modes. This implies that the interesting range for which $\hat{\mu}$ should be evaluated is roughly $k \in [10^{-4}, 10]$.

## 2.4. Equation of motion in generalised coordinates

Taking the Fourier transform of the modal decomposition (2.28), inserting it into the beam equation (2.27) in Fourier space, performing an inverse Fourier transform, and taking the scalar product (2.30) with the eigenmodes $\chi_j$, gives a set of differential equations for the generalised coordinates $q_i(t)$ that can be written in vector form as

$$ M\ddot{q} + Kq = -(M_a\ddot{q} + C_a\dot{q} + K_aq), \quad (2.31) $$

where $q(t) = [q_1(t), q_2(t), \dots]^T$ is the vector of generalised coordinates, the left-hand side corresponds to the unforced beam equation, and the right-hand side corresponds to the forcing by pressure forces. Here, we have omitted the external forcing term representing the seismic forcing.

The matrix $M = ml$, with $l$ the identity matrix, is the mass matrix, and $K = b \text{diag}(\{k_j^4\})$ is the stiffness matrix (where $\text{diag}(\{k_j^4\})$ means the diagonal matrix with $k_1^4, k_2^4$, etc. on the diagonal). Equating the left-hand side of (2.31) to zero allows us to recover the eigenmodes and the eigenfrequencies $\omega_j = \sqrt{bk_j^4/m}$ of the beam in vacuum.

The matrices $M_a, C_a$ and $K_a$ correspond to the added mass, added damping and added stiffness matrices, respectively. Their coefficients can be calculated as

$$ (M_a)_{ij} = \langle \chi_i, \mathcal{F}^{-1}[\hat{\mu}\hat{\chi}_j] \rangle = \frac{1}{2\pi} \langle \chi_i, \mu \star \chi_j \rangle, \quad (2.32a) $$

$$ (C_a)_{ij} = \langle \chi_i, \mathcal{F}^{-1}[2ik\hat{\mu}\hat{\chi}_j] \rangle = \frac{1}{2\pi} \langle \chi_i, 2\mu \star \chi_j' \rangle, \quad (2.32b) $$

$$ (K_a)_{ij} = \langle \chi_i, \mathcal{F}^{-1}[-k^2 \hat{\mu} \hat{\chi}_j] \rangle = \frac{1}{2\pi} \langle \chi_i, \mu \star \chi_j'' \rangle, \quad (2.32c) $$

where $\chi_j'$ and $\chi_j''$ respectively denote first and second derivatives of $\chi_j$, $\mu = \mathcal{F}^{-1}[\hat{\mu}]$, and we have used the property (2.19) of the convolution product.
We see here that knowledge of $\hat{\mu}(k)$ or $\mu(z)$ is enough to compute the matrices $M_a, C_a$ and $K_a$. In addition, these matrices depend only on the geometry of the problem (i.e. the wall distances $l_x$ and $l_y$) and not on the forcing. Finally, note that the Helmholtz problem (2.20) depends only on $k^2$, which means that $\hat{\mu}(k)$ and $\mu(z)$ are real and even functions. The coefficients of the matrices $M_a, C_a$ and $K_a$ are thus all real, as expected.

### 3. Single cylinder

In this section, the numerical calculations performed for the single-cylinder geometry using the code implemented on FreeFEM++ (Hecht 2012a,b) will be presented. FreeFEM++ requires the weak formulation of problems. Problem (2.20)–(2.23) is then implemented as

$$ \int_{\partial C} v \sin(\theta) \, dS - \int_{D} (\nabla \hat{p} \cdot \nabla v + k^2 \hat{p}v) \, dV = 0, \quad \forall v, \quad (3.1) $$

where $\partial C$ denotes the rod’s boundaries, and $D$ denotes the fluid domain (see Appendix A).

#### 3.1. Meshing

As described previously, the numerical resolution of the problem takes advantage of the Fourier transform. This allows us to solve several two-dimensional problems, instead of solving a three-dimensional problem.

In this subsection, we present numerical solutions of the Helmholtz problem described by (2.20)–(2.23). Our objective is to compute the value of the function $\hat{\mu}(k; l_x, l_y)$ for different values of $k, l_x$ and $l_y$. Thanks to the symmetries of the problem, the equations may be solved on a quarter of the domain only (figure 3).

Before showing the results of these numerical computations, some considerations have to be made on the discretisation of the domain. A convergence study has been performed in order to find the optimal meshing for different combinations of wall distance and wavenumber. Typically, the size of the meshes must be decreased when the wavenumber $k$ is increased. In addition, the refinement of meshes will increase as the enclosure size decreases and when getting closer to the cylinder.

Explored values are $k \in [10^{-5}, 10^3]$ and $l_x, l_y \in [1.25, 50]$. The variations of the number of cells in the mesh with $l_x$ and $l_y$ do not depend on $k$: number of cells is maximum for $l_x = l_y = 5$. This maximum is equal to 864 for $k \le 1$, and increases quadratically for $k > 1$. The FreeFEM++ code used in this study is included in Appendix A. Raw data can be found in the csv file attached to this paper; it contains the 2448 combinations of $l_x, l_y$ and $k$ values, and the resulting $\hat{\mu}$ values.

#### 3.2. General trends

Figure 4 shows the dependence of the added mass $\hat{\mu}$ on the wavenumber $k$ for different values of the wall distances $l_x$ and $l_y$. It shows that the added mass does not change substantially for $k \lesssim 10^{-2}$ or $k \gtrsim 10^2$, for the enclosure sizes considered. This was expected because the only scales in the problem are the radius $a = 1$ and the enclosure lengths $l_x$ and $l_y$, which are of order 1. Hence for $l_x$ and $l_y$ fixed, most variations of the added mass are expected to occur near $k \approx 1$. Note that the slenderness ratio of cylinders in a reactor core assembly is very large, and the range of relevant wavenumbers in this
R. Capanna, G. Ricciardi, E. Sarrouy and C. Eloy

<figure><img src="image_6.png" /><figcaption>Figure 3. Helmholtz problem represented in the FEM domain: limit conditions and example mesh for $k = 0.01$, $l_x = 1.5$ and $l_y = 2$.</figcaption></figure>

context is $k \in [3 \times 10^{-4}, 0.3]$, as stressed already in § 2.3. For the sake of clarity, the raw data used for the plots in figure 4 are published online joint to this paper.

### 3.3. Slender body limit

In the asymptotic limit of large wavelengths ($k \ll 1$) and large enclosure sizes ($l_x, l_y \gg 1$), we expect to recover the result of slender body theory (Lighthill 1960b) given by (1.2), which in dimensionless form is

$$f_y = -m_a (\partial_t + \partial_z)^2 w, \qquad (3.2)$$

with $m_a = \pi$, the added mass per unit length of a circular cylinder in dimensionless units. By analogy with (2.26), one thus expects

$$\lim_{k \to 0} \hat{\mu} = \pi, \quad \text{for } l_x, l_y \gg 1, \qquad (3.3)$$

and this is indeed what is found numerically when $l_x \gtrsim 10$, $l_x \gtrsim 10$ and $k \lesssim 10^{-2}$ (figure 4c). This validates our numerical results for small wavenumbers $k$, in the limit of negligible enclosure sizes.

### 3.4. Comparison with cylinder in annular enclosure

In the previous section, the numerical computations of the added mass have been validated in the limit of long wavelengths and large enclosure size by comparison with the slender body theory of Lighthill (1960b). In this section, we perform further comparisons with the theoretical results of Chen (1985) on a single cylinder confined in a concentric
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 4. Value of $\hat{\mu}$ as function of the wavenumber $k$ for different values of the wall distances $l_x$ and $l_y$. Dots denote computed values, and the solid line shows the fit (3.6), valid for $k \ge 10$. Wall distances are: (a) $l_x \le 5$, $l_y \le 5$; (b) $l_x \le 5$, $l_y > 5$; (c) $l_x > 5$, $l_y > 5$; (d) $l_x > 5$, $l_y \le 5$.</figcaption></figure>

annular chamber. In a potential flow, the added mass of such a cylinder is

$$ \hat{\mu} = \pi \frac{a_{ext}^2 + 1}{a_{ext}^2 - 1}, \qquad (3.4) $$

where $a_{ext}$ is the radius of the annular enclosure in dimensionless units.

We will compare this prediction with our numerical calculations of a cylinder in a square enclosure ($l_x = l_y$). To perform this comparison, we need to relate the wall distances $l_x = l_y$ to the annular radius $a_{ext}$. One possibility is to use the hydraulic diameter equivalence, which leads to an annular radius

$$ a_{hydr} = 1 + \frac{4l_x^2 - \pi}{4l_x + \pi}. \qquad (3.5) $$

Alternatively, one can simply use the inscribed circle, which gives an annular radius $a_{insc} = l_x$. Finally, one can look for the best linear fit with the numerical data, which leads to $a_{fit} = 1.0775 l_x$.

Figure 5 shows a comparison between our computations of the added mass for a square enclosure ($l_x = l_y$ and $k \ll 1$) and the added mass predicted by (3.4) for an annular enclosure. This comparison shows that the effect of the square enclosure is equivalent to an annular enclosure with radius $a_{fit} = 1.0775 l_x$. The hydraulic diameter equivalence gives a good estimate of the effect of the enclosure. This validates our numerical calculations in a square enclosure geometry in the limit $k \ll 1$.
R. Capanna, G. Ricciardi, E. Sarrouy and C. Eloy

<figure><img src="image_11.png" /><figcaption>Figure 5. Comparison between the computed added mass of a single cylinder in a square enclosure and the theoretical predictions of Chen (1985) for a concentric annular enclosure. (a) Geometry of the square enclosure ($l_x = l_y = 1.25$) superimposed with the three proposed annular enclosure radii: $a_{hydr}$, $a_{insc}$ and $a_{fit}$. (b) Comparison between our numerical results (dots) for $l_x = l_y$ and the prediction given by (3.4) for the three proposed radii (solid lines).</figcaption></figure>

### 3.5. Confinement effects

Figure 4 shows that $\hat{\mu}$ does change significantly for small values of $k$, but still depends on the wall distances $l_x$ and $l_y$ when $l_x \lesssim 5$ or $l_y \lesssim 5$. Two extreme cases can be considered: (1) a small value of the wavenumber $k = 10^{-5}$; and (2) a large value $k = 10$.

Figure 6 shows how the added mass depends on both $l_x$ and $l_y$ for the two extreme wavenumbers considered ($k = 10^{-5}$ and $k = 10$). For $k = 10^{-5}$, the enclosure length along the x-direction has much more influence than its counterpart along the y-direction (figure 6a). For $k = 10$, it is the opposite. Note, however, that these effects have different magnitudes: for $k = 10^{-5}$, the value of $\hat{\mu}$ roughly triples between $l_x = 10$ and $l_x = 1.25$, while for $k = 10$, $\hat{\mu}$ varies by only 15% between $l_y = 10$ and $l_y = 1.25$.

Figure 4 thus shows that for large wavenumbers, the dependency on the enclosure size is weak. Assuming that $\hat{\mu}$ is independent of $l_x$ and $l_y$ for large $k$ values, we can find a power-law approximation of $\hat{\mu}$ displayed as black solid line in figure 4:

$$ \hat{\mu} \approx 2.5 k^{-1}, \quad k \ge 10. \qquad (3.6) $$

For small wavenumbers, although the cylinder can move only in the y-direction, it is the wall distance along $x$ that has the strongest influence. To better understand this paradoxical result, we plot the velocity field in the Fourier space (figure 7). For $k = 10^{-5}$, a strong flow in the y-direction (opposite to the cylinder displacement) is observed along the wall $x = l_x$ around the cylinder diameter ($y \approx 0$), as shown in figure 7(a). For large wavenumbers, however, no such flow is observed (figure 7b). By virtue of mass conservation, when the wall distance along $x$ is reduced, this flow increases. This is why there is a strong dependency of the fluid-structure interaction force on the x-direction confinements for small wavenumbers.

Finally, for all the wavenumbers that have been considered and for all enclosure sizes, the added mass $\hat{\mu}$ appears to increase as the enclosure lengths decrease
Seismic response of assemblies

<figure><img src="image_12.png" /><figcaption>Figure 6. Confinement effects for (a) small wavenumbers and (b) large wavenumbers.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 7. Velocity and pressure fields in Fourier space for both small and large wavenumbers.</figcaption></figure>

(in both directions). This observation is in agreement with the literature on channel flows by Chen & Wambsganss (1972), Chen (1985) and Païdoussis & Pettigrew (1979).

The calculations performed for a single-cylinder geometry allowed us to prove the reliability of the proposed simplified model. The introduction of the potential flow theory and the use of Fourier transforms are key to relate directly the cylinder displacement and the resulting pressure force.

This approach leads to fast calculations for any kind of geometry. In the next section, it will be applied to a group of cylinders and compared to experimental results in the same geometry.
<figure><img src="image_14.png" /><figcaption>Figure 8. Scheme of the ICARE experimental facility. LDV, laser Doppler velocimetry.</figcaption></figure>

## 4. Assemblies of cylinders

This section is dedicated to a geometry with multiple cylinders. This geometry is built to represent the four surrogate fuel bundles of the experimental facility ICARE that we used to validate our results.

### 4.1. ICARE experimental apparatus

The ICARE experimental facility is made of a closed water loop powered by a centrifugal pump, the test section, a compensation tank and a heat exchanger (figure 8). The test section is placed in vertical position with a square section measuring 22.5 cm × 22.5 cm, and hosts up to four fuel assemblies arranged in a 2 × 2 lattice. The length of the fuel assemblies is $L = 2.57$ m. Each of them constitutes a square lattice of 8 × 8 rods, in which there are 60 rods simulating fuel rods made of acrylic, and four stainless steel guide tubes. The guide tubes are empty inside, and they are welded to five metallic spacer grids along the length of the assembly. The guide tubes have a structural function since they give rigidity to the assembly and they hold together fuel rods. Each assembly has a section measuring 10.1 × 10.1 cm², and the 60 rods constituting it have diameter $2A = 9$ mm with pin pitch $P/(2A) = 1.39$. The top and bottom of the assembly are fixed rigidly to the guide tubes, and they are fixed to the test section through the lower and upper support plates.

A hydraulic actuator applies dynamic forcing to one of the four fuel assemblies (assembly[1] in figure 8). This forcing in one dimension aims to simulate an earthquake. The actuator is screwed into one of the grids of the assembly, and a force sensor is installed in between the actuator and the stem.

The test section is equipped with 24 linear variable differential transformer (LVDT) position sensors to measure the displacement of each grid in two directions (figure 9). For more details on the ICARE experiment, please refer to Capanna *et al.* (2019).
Seismic response of assemblies

<figure><img src="image_15.png" /><figcaption>Figure 9. Scheme of the displacement sensors on the ICARE test section.</figcaption></figure>

## 4.2. Problem geometry and meshing

The ICARE geometry is implemented in the code to allow for comparison with experiments with large ($C = 8$ mm) or small ($C = 4$ mm) confinements (figures 10(a) and 10(b), respectively). Confinement size denotes the length of the gap between different assemblies and between the assemblies and the walls. For both confinement sizes, the distance in between the rods of the same assembly remains the same.

Assemblies are numbered as depicted in figure 10(a). As in the experiments, the rods of assembly[1] are forced along the x-direction, and they all exhibit the same displacement $w_x^{[1]}(z, t)$. Those rods are assumed not to move along the y-direction, and the rods in the other assemblies are assumed not to move.

Based on § 3, the Helmholtz problem to solve is

$$ \Delta \hat{p} = k^2 \hat{p}, \quad \text{in the domain,} \qquad (4.1a) $$

$$ \left. \frac{\partial \hat{p}}{\partial n} \right|_{\partial C_w} = 0, \quad \text{on the walls,} \qquad (4.1b) $$

$$ \left. \frac{\partial \hat{p}}{\partial n} \right|_{\partial C^{[1]}} = \hat{\gamma}_x^{[1]}(k, t) \cos(\theta), \quad \text{for assembly [1],} \qquad (4.1c) $$

$$ \left. \frac{\partial \hat{p}}{\partial n} \right|_{\partial C^{[n]}} = 0, \quad \text{for other assemblies,} \qquad (4.1d) $$
<figure><img src="image_16.png" /><figcaption>Figure 10. Meshes for the two configurations simulating the ICARE set-up. These two geometries differ only by the value of the confinement C, as pictured. (a) Four assemblies in large confinement. (b) Four assemblies in small confinement.</figcaption></figure>

where the superscripts [n] refer to the assembly number, $\partial C^{[n]}$ denotes the boundaries of rods of assembly [n] (with $\theta$ the local polar angle for each rod), and

$$ \hat{\gamma}_{x}^{[1]} = \mathcal{F}[(\partial_t + \partial_z)^2 w_{x}^{[1]}]. \quad (4.2) $$

As in § 3, the solution $\hat{p}$ of the Helmholtz problem (4.1) is proportional to $\hat{\gamma}_x^{[1]}$, and one can define added masses $\hat{\mu}_x^{[n]}$ and $\hat{\mu}_y^{[n]}$ from the forces exerted on each assembly:

$$ \hat{f}_{x}^{[n]}(k, t) = - \int_{\partial C^{[n]}} \hat{p} \cos(\theta) \, d\theta = - \hat{\mu}_{x}^{[n]} \hat{\gamma}_{x}^{[1]}, \quad (4.3a) $$

$$ \hat{f}_{y}^{[n]}(k, t) = - \int_{\partial C^{[n]}} \hat{p} \sin(\theta) \, d\theta = - \hat{\mu}_{y}^{[n]} \hat{\gamma}_{x}^{[1]}. \quad (4.3b) $$

In practice, $\hat{\mu}_x^{[n]}$ and $\hat{\mu}_y^{[n]}$ can be calculated by assuming $\hat{\gamma}_x^{[1]} = 1$ and numerically solving the Helmholtz problem (4.1) for $\hat{p}$. The convergence of the calculations depends on the mesh size. We thus conducted a convergence study and found that convergence is ensured when the walls are divided into at least 150 sections, and the cylinder borders into 50 sections. The calculations presented here have been performed with 200 divisions on the external walls and 75 divisions on the cylinder borders, which corresponds to approximately 500 000 cells in the domain and a computing time of about 60 s.

Figure 11 shows the results of these computations for $k \in [10^{-5}, 10]$. In Appendix B, we detail these results for individual rods. As in the single-cylinder case, the added masses do not change much for extreme values of k, i.e. $k \lesssim 10^{-2}$ or $k \gtrsim 1$, except for assembly [1], the forced assembly. For large k, the added masses converge to zero.

### 4.3. Equations of motion in real space

The functions $\hat{\mu}_x^{[n]}(k)$ and $\hat{\mu}_y^{[n]}(k)$ can be interpolated using piecewise sixth-order polynomials, which can then be used to calculate $\mu_x^{[n]}(z)$ and $\mu_y^{[n]}(z)$ with an inverse
<figure><img src="image_17.png" /><figcaption>Figure 11. Added masses $\hat{\mu}_x^{[n]}$ (blue solid line) and $\hat{\mu}_y^{[n]}$ (red dashed line) on each assembly (in Fourier space) for different values of $k$.</figcaption></figure>

Fourier transform (Appendix C). By analogy with the calculations for a single cylinder done in § 2.4, we can define added mass, added damping and added stiffness matrices for the assemblies too. There are 24 such matrices: 3 kinds (mass **M**, damping **C**, stiffness **K**), 2 axes ($x$ and $y$, subscript $\circ$) and 4 assemblies, noted with the superscript $[n]$. Similarly to (2.32a-c), these matrices are expressed using the added masses $\mu_{\circ}^{[n]}(z)$:

$$ (\mathbf{M}_{\circ}^{[n]})_{ij} = \langle \chi_i, \mathcal{F}^{-1}[\hat{\mu}_{\circ}^{[n]} \hat{\chi}_j] \rangle = \frac{1}{2\pi} \langle \chi_i, \mu_{\circ}^{[n]} \star \chi_j \rangle, \quad (4.4a) $$

$$ (\mathbf{C}_{\circ}^{[n]})_{ij} = \langle \chi_i, \mathcal{F}^{-1}[2ik\hat{\mu}_{\circ}^{[n]} \hat{\chi}_j] \rangle = \frac{1}{2\pi} \langle \chi_i, 2\mu_{\circ}^{[n]} \star \chi_j' \rangle, \quad (4.4b) $$

$$ (\mathbf{K}_{\circ}^{[n]})_{ij} = \langle \chi_i, \mathcal{F}^{-1}[-k^2 \hat{\mu}_{\circ}^{[n]} \hat{\chi}_j] \rangle = \frac{1}{2\pi} \langle \chi_i, \mu_{\circ}^{[n]} \star \chi_j'' \rangle, \quad (4.4c) $$

where $\chi_i(z)$ is the $i$th beam eigenmode defined in (2.29).

Because it is forced at mid-height, we can assume that assembly [1] moves mainly in its first beam eigenmode, i.e. $w_x^{[1]}(z, t) = q_x^{[1]}(t) \chi_1(z)$. We can assume further that the added mass and added stiffness matrices are diagonal, and that the added damping is negligible (Appendix C). With these hypotheses, the equation of motion for the $x$-component of assembly [1] is

$$ (m_s + m_x^{[1]})\ddot{q}_x^{[1]} + c_v \dot{q}_x^{[1]} + (k_s + k_x^{[1]})q_x^{[1]} = f_e e^{i\omega t}, \quad (4.5) $$

where $m_s = 64m$ is the dimensionless mass of one assembly per unit length, $k_s = 64bk_1^4$ is its stiffness, $c_v$ accounts for the damping induced by the viscosity effect, $m_x^{[1]}$ and $k_x^{[1]}$ denote the component (1, 1) of the added matrices $\mathbf{M}_x^{[1]}$ and $\mathbf{K}_x^{[1]}$, and $f_e e^{i\omega t}$ is the external
<table><caption>Table 1. Added mass matrix coefficients for small and large confinements.</caption><thead><tr><th>Confinement</th><th><i>m</i><sub>x</sub><sup>[1]</sup></th><th><i>m</i><sub>x</sub><sup>[2]</sup></th><th><i>m</i><sub>x</sub><sup>[3]</sup></th><th><i>m</i><sub>x</sub><sup>[4]</sup></th><th><i>m</i><sub>y</sub><sup>[1]</sup></th><th><i>m</i><sub>y</sub><sup>[2]</sup></th><th><i>m</i><sub>y</sub><sup>[3]</sup></th><th><i>m</i><sub>y</sub><sup>[4]</sup></th></tr></thead><tbody><tr><td>Large</td><td>350</td><td>-22.1</td><td>72.7</td><td>18.6</td><td>-17.0</td><td>20.5</td><td>-20.5</td><td>26.2</td></tr><tr><td>Small</td><td>384</td><td>-24.0</td><td>79.0</td><td>21.3</td><td>-22.4</td><td>25.3</td><td>-25.3</td><td>29.3</td></tr></tbody></table>

forcing projected onto the first eigenmode, with $f_e$ a scalar. The $m_x^{[1]}$ and $k_x^{[1]}$ coefficients represent the added mass and stiffness encountered by an assembly that moves in its first mode in either the x-direction or the y-direction. Here, we consider only a weak coupling between the assemblies, such that we can neglect the hydrodynamic forces due to the motion of the other assemblies on assembly [1] (the validity of this assumption will be assessed below). The solution of (4.5) is simply

$$q_x^{[1]} = h_x^{[1]} f_e e^{i\omega t}, \quad (4.6)$$

with the transfer function

$$h_x^{[1]}(\omega) = \frac{1}{-\omega^2(m_s + m_x^{[1]}) + i\omega c_v + k_s + k_x^{[1]}}. \quad (4.7)$$

It is possible to estimate the motion of assemblies [2], [3] and [4] by considering the superimposition of the two following cases: still rods undergoing fluid forces induced by assembly [1] on the one hand, and moving rods undergoing fluid forces generated by their own movement while other assemblies stand still:

$$(m_s + m_x^{[1]}) \ddot{q}_o^{[n]} + c_v \dot{q}_o^{[n]} + (k_s + k_x^{[1]}) q_o^{[n]} = -(-\omega^2 m_o^{[n]} + k_o^{[n]}) q_x^{[1]}, \quad (4.8)$$

where the coefficients $m_x^{[1]}$, $c_v$ and $k_x^{[1]}$ appear on the left-hand side because, by symmetry, they also correspond to the pressure force induced by the motion $q_o^{[n]}$ on the assembly [n] itself. The right-hand side represents the action of assembly [1] motion via the fluid.

From (4.4c), we see that the added stiffness is $k_o^{[n]} = -k_1^2 m_o^{[n]}$, with $k_1 = \pi/l$ (since $\chi_i'' = -k_i^2 \chi_i$). We will now assume that $\omega \gg k_1$, which is true here because typical values of $\omega$ are around 0.3, and $l \approx 571$. This means that the added stiffness term $k_o^{[n]}$ can be neglected compared to the added mass term $-\omega^2 m_o^{[n]}$ in (4.5) and (4.8). The values of the added masses $m_o^{[n]}$ are given in table 1. From this table, we see that the ratio of added masses $m_o^{[n]}/m_x^{[1]}$ is at most 21% (for $n > 1$). This justifies neglecting the hydrodynamic forces due to the motion of the assemblies [n] onto assembly [1], which can be estimated as the square of this ratio and is thus of order 4% at most.

Under these hypotheses, we can calculate the displacement ratio at the forcing frequency $\omega$ as

$$\frac{q_o^{[n]}}{q_x^{[1]}} \approx - \frac{-\omega^2 m_o^{[n]}}{-\omega^2(m_s + m_x^{[1]}) + i\omega c_v + k_s}, \quad (4.9)$$

which becomes, at the resonant frequency $\omega_{res}$,

$$\frac{q_o^{[n]}}{q_x^{[1]}} \approx -i\omega_{res} \frac{m_o^{[n]}}{c_v}, \quad \text{with} \quad -\omega_{res}^2(m_s + m_x^{[1]}) + k_s = 0. \quad (4.10)$$

Note that this equation is applicable only to $n \neq 1$.
## 4.4. Comparison with a periodic assembly of cylinders

The model presented above allows us to compute the added mass of an assembly of cylinders organised in a square array. A qualitative comparison can be made with an infinite assembly made of a periodic arrangement of cylinders (Pettigrew & Taylor 2003). In this case, by analogy with (3.4), the added mass of each cylinder is expressed as

$$ \hat{\mu} = \pi \frac{a_{\text{equiv}}^2 + 1}{a_{\text{equiv}}^2 - 1}, \qquad (4.11) $$

where $a_{\text{equiv}}$ is a measure of the equivalent confinement that can be approximated as

$$ a_{\text{equiv}} = (0.535 + 0.14p) p, \qquad (4.12) $$

with $p = P/A$, the period along $x$ and $y$ of the arrangement.

Using the value $P/A = 2.78$ of the ICARE assemblies, (4.11) corresponds to an added mass per cylinder $\hat{\mu} = 4.26$. Considering an assembly made of 64 cylinders, the total added mass is $m_x = 274$. This value can be compared to the value $m_x^{[1]} = 350$ obtained in our simulations (table 1). The two values are of the same order, but in the ICARE geometry, the presence of other assemblies and of the walls tends to increase the added mass.

# 5. Comparison with experimental results

This section compares results obtained using the simplified model and results from a set of experiments conducted on the ICARE facility described in § 4.1. Experiments are performed imposing a sine sweep ranging from 0 to 10 Hz with constant amplitude 1 mm and axial fluid velocity $1 \text{ m s}^{-1}$ when water is present. To catch the transfer function of a system using a swept sinusoidal excitation, the sweep rate needs to be small enough to avoid any transient phenomena and contamination of different harmonics in the system. For these experiments, a sweep rate of $0.05 \text{ Hz s}^{-1}$ is applied, respecting the international standard ISO-7626 indications. The frequency range embraces assembly [1] first and second mode in water.

## 5.1. Added mass

Let us consider the transfer function $H_{\text{exp}}^{[1]}$ between the displacement $W_x^{[1]}$ along $x$ of the grid with vertical position $Z_w$, and the force $F_g$ applied by the actuator on the grid with position $Z_f$. Since the data are acquired with a swept sinusoidal excitation, the transfer functions are calculated using the cross-correlation product, and are defined as

$$ H_{\text{exp}}^{[1]} = \left\| \frac{W_x^{[1]} \star W_x^{[1]}}{W_x^{[1]} \star F_x} \right\|, \qquad (5.1) $$

where $W_x^{[1]} \star W_x^{[1]}$ is the spectral autocorrelation of the displacement signal of the first assembly along the $x$-direction, and $W_x^{[1]} \star F_x$ is the spectral cross-correlation between the displacement signal of the first assembly and the force imposed on the first assembly. Cross-correlation functions are calculated using an Hamming windowing filter (5000 points per windows) to clean up the noise from the signal.
The theoretical transfer function modulus for a one-degree-of-freedom system is given by

$$H_{theo}^{[1]}(\Omega) = \frac{1}{\sqrt{(K_e - \Omega^2 M_e)^2 + \Omega^2 C_e^2}}, \qquad (5.2)$$

where $M_e$, $C_e$ and $K_e$ are the coefficients of mass, damping and stiffness identified, and $\Omega$ is the circular frequency.

Assembly [1] is excited around its first mode; to identify its dynamical properties, one has to find $M_e$, $C_e$ and $K_e$ that make $H_{theo}^{[1]}$ fit the experimental data $H_{exp}^{[1]}$. This can be achieved easily by a conjugate gradient optimisation method.

Using an Euler-Bernoulli beam model for the rods as described in (2.1), one can write the virtual work of assembly [1] as

$$\begin{aligned} M_t \int_0^L \partial_T^2 W_x^{[1]} \delta W \, dZ + C_t \int_0^L \partial_T W_x^{[1]} \delta W \, dZ + B_t \int_0^L \partial_Z^2 W_x^{[1]} \partial_Z^2 \delta W \, dZ \\ = \int_0^L F_g \delta Z_f \, dZ, \end{aligned} \qquad (5.3)$$

with

$$M_t = M_s + M_f, \quad B_t = B + B_f, \qquad (5.4a,b)$$

where $M_f$ and $B_f$ are respectively the participation of the fluid into mass and stiffness, $C_t$ is a linear damping coefficient, and the sum accounts for the 64 rods inside the assembly.

When the structure vibrates on its jth natural mode,

$$W_x^{[1]}(Z, T) = N_j(T) X_j(Z), \qquad (5.5)$$

with $X_j(Z) = \chi_j(Z/A)$, one can rewrite the energy of the system as

$$M_t \int_0^L X_j^2 \, dZ \, \partial_T^2 N_j + C_t \int_0^L X_j^2 \, dZ \, \partial_T N_j + B_t \int_0^L (\partial_Z^2 X_j)^2 \, dZ \, N_j = F_g X_j(Z_f). \qquad (5.6)$$

Noting that the experimental signal is

$$W_g^{[1]} = N_1(T) X_1(Z_w), \qquad (5.7)$$

one can relate the mass per unit of length to the identified coefficient

$$M_t = \frac{X_1(Z_f) X_1(Z_w)}{\int_0^L X_1^2 \, dZ} M_e. \qquad (5.8)$$

Measurements take into account the effect of the fluid added mass $M_f$ and the mass of the structure $M_s$, therefore to isolate the fluid added mass, the participation of the structure is subtracted, based on experimental results obtained in air. Finally, one then obtains the results listed in table 2 where added mass values are provided for the whole assembly,

$$M_x^{[1]} = M_f. \qquad (5.9)$$

To get a comparison with results from numerical simulation obtained in § 4, the following transformation is applied

$$M_x^{[1]} = \rho A^2 m_x^{[1]}, \qquad (5.10)$$

where $\rho = 1000 \text{ kg m}^{-3}$ is the water density.
<table><caption>Table 2. Values of added mass $M_x^{[1]}$.</caption><thead><tr><th>M<sub>x</sub><sup>[1]</sup></th><th>Experiment</th><th>Simulation</th></tr></thead><tbody><tr><td>Large confinement</td><td>4.06 kg m<sup>-1</sup></td><td>7.09 kg m<sup>-1</sup></td></tr><tr><td>Narrow confinement</td><td>10.36 kg m<sup>-1</sup></td><td>7.78 kg m<sup>-1</sup></td></tr></tbody></table>

Table 2 shows the values of added mass for large and narrow confinement given by experiments and simulations. Both agree on the increase of added mass when the confinement gets narrower. Experiments give a large increase, more than double, whereas simulations predict an increase of about 10%. Although the order of magnitude is respected, the differences observed could cast doubt on the validity of the model. To analyse these results properly, one has to keep in mind that the methodology used to obtain the experimental results is subject to uncertainties coming from the accuracy of measurement, the optimisation process, and the fact that we try to identify constants of a linear one-degree-of-freedom system from a complex nonlinear system. As a result, one could estimate that the added mass coefficient given by the simulations is a reasonable estimation.

The added stiffness is orders of magnitude smaller than the ‘in air’ structural stiffness: the added stiffness estimated from simulations in § 4 is about 4.7 N m², while the measured flexural stiffness of a fuel assembly is about 2300 N m². As a consequence, we cannot extract the added stiffness from the total stiffness measured experimentally.

## 5.2. Coupling

To characterise the coupling between the fuel assemblies induced by the fluid, let us consider the transfer function between the imposed displacement and the displacement of the other fuel assemblies at the third grid level in the two transverse directions. Similarly to (5.1), transfer functions are calculated as

$$H_{\circ}^{[n]} = \left\| \frac{W_{\circ}^{[n]} \star W_{\circ}^{[n]}}{W_{\circ}^{[n]} \star W_{x}^{[1]}} \right\|, \quad \circ \in \{x, y\}, \qquad (5.11)$$

where $W^{[n]}$ is the displacement signal of the n-th assembly along the $\circ$ direction.

For each transfer function, the maximum value – which occurs at the first natural frequency – is accounted for. These experimental results are compared to the displacement ratio given by (4.10), with the added mass coefficients $m_{\circ}^{[n]}$ given in table 1, and $c_v$ and $\omega_{res}$ estimated from identification of dynamical parameters using (5.6) as follows:

$$C_v = \frac{X_1(Z_f) X_1(Z_w)}{\int_0^L X_1^2 dZ} C_e, \quad c_v = \frac{C_v}{\rho UA}, \quad \omega_{res} = \frac{A}{U} \sqrt{\frac{K_e}{M_e}}. \qquad (5.12a-c)$$

Experimental values give $c_v = 69.2$ and $\omega_{res} = 0.12$.

Figure 12 compares graphically the displacement ratios as observed in the experiments with the ratios predicted by (4.10) with these values of $c_v$ and $\omega_{res}$. Experiments show a ratio of about 10% for most of the fuel assemblies in the two directions except for fuel assembly [3], which shows a coupling of about 20% in the direction of forcing.
<figure><img src="image_18.png" /><figcaption>Figure 12. Displacement ratios obtained by experimental data (a) and approximated from numerical simulations (b) for four fuel assemblies in large confinement at first natural frequency.</figcaption></figure>

Although simulations seem to underestimate the values of coupling, they give a good prediction of the order of magnitude. Moreover, the general pattern is well reproduced, with a significantly more important coupling for fuel assembly [3] in the direction of forcing. Therefore, the methodology developed seems to be appropriate to simulate the coupling between fuel assemblies.

## 6. Conclusions

In this paper, a new model for fluid-structure interaction of PWR fuel assemblies is proposed. The model used the well-known potential flow theory, and the equations are solved by using a Fourier approach. This approach leads to the important result of relating the fluid-structure interaction forces directly to the displacement of the structure itself. This result is an important achievement since it allows us to reduce computational time drastically, avoiding the necessity to solve fluid equations. As a drawback, it should be considered that the assumption of potential flow leads us to neglect completely the fuel grids, which introduce vorticity in the flow.

Calculations with a potential flow model are first performed for a single-cylinder geometry, and results have been compared with reference works in the literature. The model showed perfect agreement with the slender body model, and it also shows good predictions with respect to the confinement size. The model has been validated, thus the mathematical approach used to solve the equations has been demonstrated to be consistent and reliable. The model was thus improved for multiple-cylinder geometries.

Finally, in order to validate the multiple-cylinder model, calculations simulating the ICARE geometry are discussed and compared with experimental results. Modal parameters are identified as functions of the flow rates and compared with experimental ones. The model fits reasonably experimental data for the added mass, and gives a good estimation of the coupling between fuel assemblies.

Many perspectives are opened as a result of this work. As discussed in this paper, the model that has been implemented does not account for viscous forces. An improvement that would be of fundamental importance is to introduce with empirical factors the viscous
Seismic response of assemblies

forces in the model. One possible way to introduce viscosity in the model is to consider
the axial water flow as the sum of two flows: a potential bulk flow and a viscous flow that
has a defined empirical distribution. In this way, the viscous forces would be taken into
account with an empirical formulation, and it would be extremely interesting to assess the
effects of such improvements, especially on the added damping coefficient estimation.

**Supplementary material.** Supplementary material is available at https://doi.org/10.1017/jfm.2022.433.

**Funding.** The authors are grateful for the financial support of Electricité de France (EDF) and FRAMATOME.

**Declaration of interests.** The authors report no conflict of interest.

Author ORCIDs.

Roberto Capanna https://orcid.org/0000-0001-9580-5140;

Christophe Eloy https://orcid.org/0000-0003-4114-7263.

Appendix A. FreeFEM++ code for single cylinder

1
real k = 0.01;
real lx = 1.5;
real ly = 2;
real A = 1.;

6 // Geometry

Cxb(t=0, 1){x=A+(lx-A)*t^2); y=0*t; label=1;}
border Cyr(t=0, 1){x=lx+0*t; y=tan(atan(ly/lx)*t)/tan(atan(ly/lx))*ly; label=2;}
border Cxt(t=0, 1){x=tan(atan(lx/ly)*(1-t))/tan(atan(lx/ly))*lx; y=ly+0*t; label=2;}
border Cyl(t=0, 1){x=0*t; y=A+(ly-A)*(1-t)^2); label=2;}

11 border Beam(t=0, pi/2){x=A*cos(t); y=A*sin(t); label=3;}

// Mesh

real ds = 1/k/10;
int NbA = max(ceil(pi/2*A/ds), 21.);
int Nbx = ceil(atan(lx/ly)/(pi/2)*(NbA))+1;
int Nby = ceil(atan(ly/lx)/(pi/2)*(NbA))+1;
int NbxmA = ceil((lx-A)/lx*NbA)+1;
int NbymA = ceil((ly-A)/ly*NbA)+1;

mesh cMesh = buildmesh(Cxb(NbxmA) + Cyr(Nby) + Cxt(Nbx) + Cyl(NbymA) + Beam(-NbA));

21

// FE

fespace Vh(cMesh, P2);
Vh u, v;

26

// Resolution

solve pFourier(u,v)
    = - int2d(cMesh)(u*v*k^2 + dx(u)*dx(v) + dy(u)*dy(v))
        + int1d(cMesh,3)(y/A*v)
        + on(1,u=0);

31

// mu computation

real mu = 0;
mu = -4*int1d(cMesh,3)(-u*y/A);

**Appendix B. Qualitative analysis**

After solving numerically the Helmholtz problem (4.1) for $\hat{\gamma}_x^{[1]} = 1$, the resulting forces on each rod $q$ of each assembly $n$ can be evaluated in both the $x$- and $y$-directions.
R. Capanna, G. Ricciardi, E. Sarrouy and C. Eloy

<figure><img src="image_19.png" /><figcaption>Figure 13. Added mass $\hat{\mu}_x^{[n,q]}$ along $x$ for each rod for four values of $k$, as labelled. Solid circles note positive values, hollow ones negative values. The confinement is 8 mm.</figcaption></figure>

From this, we can calculate the associated added masses $\hat{\mu}_x^{[n,q]}$ and $\hat{\mu}_y^{[n,q]}$:

$$ \hat{\mu}_x^{[n,q]} = \int_{\partial C^{[n,q]}} \hat{p} \cos(\theta) \, d\theta, \qquad (B1a) $$

$$ \hat{\mu}_y^{[n,q]} = \int_{\partial C^{[n,q]}} \hat{p} \sin(\theta) \, d\theta, \qquad (B1b) $$

with $\partial C^{[n,q]}$ the boundary of rod $q$ in assembly $n$. Note that $\hat{\mu}_x^{[n,q]}$ and $\hat{\mu}_y^{[n,q]}$ are unchanged when $k$ is changed into $-k$ since $k$ always appears squared in the problem.

These added masses are illustrated in figures 13 and 14 for the 8 mm confinement.

Analogy with the single-cylinder case explains the positive values observed on assembly [1] for $\hat{\mu}_x^{[n,q]}$. Confinement size is about the rod diameter, which may be compared to $l_x = l_y = 3$ in §§ 3.3 and 3.4 for which $\hat{\mu} > \pi$ when $k \ll 1$. This is consistent with the values observed for assembly [1] (even though the rods are not confined properly). An interesting result is that rods in assemblies [2], [3] and [4] experience very small forces. As noted in § 3.5, this phenomenon is all the more true as $k$ becomes large: pressure perturbations remain closer to the moving rods when $k$ is large, and interactions with far away rods are almost negligible.
Seismic response of assemblies

<figure><img src="image_20.png" /><figcaption>Figure 14. Same as figure 13 for the y-component of the added mass, $\hat{\mu}_{y}^{[n,q]}$.</figcaption></figure>

Appendix C. Calculation of pressure forces in physical space

In § 4.3, we explained how to compute numerically the functions $\hat{\mu}_{\circ}^{[n]}(k)$, which are then interpolated with piecewise sixth-order polynomials. These interpolations can be used to compute the added masses $\mu_{y}^{[n]}(z)$ with an inverse Fourier transform (figure 15).

The functions $\mu_{\circ}^{[n]}(z)$ are needed in (4.4a-c) to compute the added mass, added damping and added stiffness matrices. Note that the functions $\mu_{\circ}^{[n]}(z)$ are non-zero on an interval $\Delta z \lesssim 100$ (figure 15). This means that on a z-scale of order $l = L/A \approx 571$, these functions are almost proportional to a Dirac function, such that the convolution products between $\mu_{\circ}^{[n]}$ and $\chi_j$, $\chi'_j$ or $\chi''_j$ appearing in (4.4a-c) can be approximated by simple product with $\int_{-\infty}^{\infty} \mu_{\circ}^{[n]} dz$.

As a consequence, the added mass matrix $M_{\circ}^{[n]}$ and the added stiffness matrix $K_{\circ}^{[n]}$ are almost diagonal (since $\langle \chi_i, \chi_j \rangle = \delta_{ij}$) and the diagonal terms of the added damping matrix $C_{\circ}^{[n]}$ are almost zero (since $\langle \chi_i, \chi'_i \rangle = 0$). This can be verified by computing the first $3 \times 3$ terms of the matrices $M_{\circ}^{[n]}$ and $C_{\circ}^{[n]}$ for $n = 1$ and $4$ (we chose assembly [4] because it is the assembly for which the extension of $\mu_{\circ}^{[n]}(z)$ is widest; see figure 15):

$$ M_{x}^{[1]} = \begin{pmatrix} 350 & 0 & 0.0182 \\ 0 & 347 & 0 \\ 0.0182 & 0 & 342 \end{pmatrix}, \qquad (C1a) $$
R. Capanna, G. Ricciardi, E. Sarrouy and C. Eloy

<figure><img src="image_21.png" /><figcaption>Figure 15. Added mass $\mu_x^{[n]}(z)$ (blue solid line) and $\mu_y^{[n]}(z)$ (red dashed line) as a function of $z$ for the different assemblies (large confinement).</figcaption></figure>

$$ \mathbf{M}_x^{[4]} = \begin{pmatrix} 18.6 & 0 & 0.0882 \\ 0 & 17.4 & 0 \\ 0.0882 & 0 & 15.4 \end{pmatrix}, \qquad (C1b) $$

$$ \mathbf{C}_x^{[1]} = \begin{pmatrix} 0 & -3.26 & 0 \\ 3.26 & 0 & -5.79 \\ 0 & 5.79 & 0 \end{pmatrix}, \qquad (C1c) $$

$$ \mathbf{C}_x^{[4]} = \begin{pmatrix} 0 & -0.167 & 0 \\ 0.167 & 0 & -0.272 \\ 0 & 0.272 & 0 \end{pmatrix}. \qquad (C1d) $$

Note that the matrices $\mathbf{K}_\circ^{[n]}$ satisfy $(\mathbf{K}_\circ^{[n]})_{ij} = -k_j^2 (\mathbf{M}_\circ^{[n]})_{ij}$, such that they share the same characteristic of being almost diagonal with the matrices $\mathbf{M}_\circ^{[n]}$.

REFERENCES

BAUCHAU, O.A. & CRAIG, J.I. 2009 *Structural Analysis: With Applications to Aerospace Structures*. Springer.

BLEVINS, R.D. 2015 *Formulas for Dynamics, Acoustics and Vibration*. John Wiley & Sons, Ltd.

BRENNEN, C.E. 1982 A review of added mass and fluid inertial forces. *Tech. Rep.* Naval Civil Engineering Laboratory, Sierra Madre, CA, USA.

BUAT, D. 1779 *Principles d'hydraulique*. Monsieur le Chevalier Du Buat.

CAPANNA, R. 2018 Modelling of fluid structure interaction by potential flow theory in a PWR under seismic excitation. PhD thesis, Ecole Centrale Marseille, France.

CAPANNA, R., RICCIARDI, G., ELOY, C. & SARROUY, E. 2019 Experimental study of fluid structure interaction on fuel assemblies on the ICARE experimental facility. *Nucl. Engng Des.* **352**, 110146.
CHEN, S.S. 1985 Flow-induced vibration of circular cylindrical structures. *Tech. Rep. No.* ANL-85-51.
Argonne National Lab. (ANL), Argonne, IL (United States).

CHEN, S.S. & WAMBSGANSS, M.W. 1972 Parallel-flow-induced vibration of fuel rods. *Nucl. Engng Des.* **18** (2), 253–278.

CLASSEN, P. 1972 Die 'hydrodynamische masse' eines rohrs in einem rohrbundel. *Forsch. Ing. Wes.* **38** (2), 33–37.

DE LANGRE, E., PAÏDOUSSIS, M.P., DOARÉ, O. & MODARRES SADEGHI, Y. 2007 Flutter of long flexible cylinders in axial flow. *J. Fluid Mech.* **571**, 371–389.

DE MARIO, E. & STREET, B.D. 1989 Nuclear fuel rod grid spring and dimple structures. US Patent 4,803,043.

DE RIDDER, J., DEGROOTE, J., VAN TICHELEN, K. & VIERENDEELS, J. 2013 Modal characteristics of a flexible cylinder in turbulent axial flow from numerical simulations. *J. Fluids Struct.* **43**, 110–123.

DE RIDDER, J., DEGROOTE, J., VAN TICHELEN, K. & VIERENDEELS, J. 2017 Predicting modal characteristics of a cluster of cylinders in axial flow: from potential flow solutions to coupled CFD-CSM calculations. *J. Fluids Struct.* **74**, 90–110.

DE RIDDER, J., DOARÉ, O., DEGROOTE, J., VAN TICHELEN, K., SCHUURMANS, J. & VIERENDEELS, J. 2015 Simulating the fluid forces and fluid-elastic instabilities of a clamped-clamped cylinder in a turbulent flow. *J. Fluids Struct.* **55**, 139–154.

DOARÉ, O., MANO, D. & BILBAO LUDENA, J.C. 2011a Effect of spanwise confinement on flag flutter: experimental measurements. *Phys. Fluids* **23**, 11704.

DOARÉ, O., SAUZADE, M. & ELOY, C. 2011b Flutter of an elastic plate in a channel flow: confinement and finite-size effects. *J. Fluids Struct.* **27** (1), 76–88.

ELOY, C., DOARÉ, O., DUCHEMIN, L. & SCHOUVEILER, L. 2010 A unified introduction to fluid mechanics of flying and swimming at high Reynolds number. *Exp. Mech.* **50** (9), 1361–1366.

ELOY, C., LAGRANGE, R., SOUILLEZ, C. & SCHOUVEILER, L. 2008 Aeroelastic instability of cantilevered flexible plates in uniform flow. *J. Fluid Mech.* **611**, 97–106.

ELOY, C., SOUILLEZ, C. & SCHOUVEILER, L. 2007 Flutter of a rectangular plate. *J. Fluids Struct.* **23** (6), 904–919.

GOSSELIN, F.P. & DE LANGRE, E. 2011 Drag reduction by reconfiguration of a poroelastic system. *J. Fluids Struct.* **27** (7), 1111–1123.

GUO, C.Q. & PAÏDOUSSIS, M.P. 2000 Stability of rectangular plates with free side-edges in two-dimensional inviscid channel flow. *Trans. ASME J. Appl. Mech.* **67** (1), 171–176.

HECHT, F. 2012a Freefem++ v.3.19-1. Université Pierre et Marie Curie. http://www.freefem.org/ff++/ftp/freefem++doc.pdf.

HECHT, F. 2012b New development in FreeFem++. *J. Numer. Maths* **20** (3–4), 251–265.

LAMB, H. 1895 *Hydrodynamics*. Cambridge University Press.

LIGHTHILL, M.J. 1960a Mathematics and aeronautics. *J. R. Aero. Soc.* **64** (595), 375–394.

LIGHTHILL, M.J. 1960b Note on the swimming of slender fish. *J. Fluid Mech.* **9** (2), 305–317.

MICHELIN, S. & SMITH, S.L. 2009 Linear stability analysis of coupled parallel flexible plates in an axial flow. *J. Fluids Struct.* **25** (7), 1136–1157.

MORISON, J.R., JOHNSON, J.W. & SCHAAF, S.A. 1950 The force distribution exerted by surface waves on piles. *J. Petrol. Tech.* **2** (5), 149–154.

PAÏDOUSSIS, M.P. 1966a Dynamics of slender cylinders in axial flow Part 1. Theory. *J. Fluid Mech.* **26** (4), 717–736.

PAÏDOUSSIS, M.P. 1966b Dynamics of slender cylinders in axial flow Part 2. Experiments. *J. Fluid Mech.* **26** (4), 737–751.

PAÏDOUSSIS, M.P. 1973 Dynamics of cylindrical structures subjected to axial flow. *J. Sound Vib.* **29** (3), 365–385.

PAÏDOUSSIS, M.P. 1979 The dynamics of clusters of flexible cylinders in axial flow: theory and experiments. *J. Sound Vib.* **65** (3), 391–417.

PAÏDOUSSIS, M.P. & PETTIGREW, M.J. 1979 Dynamics of flexible cylinders in axisymmetrically confined axial flow. *Trans. ASME J. Appl. Mech.* **46** (1), 37–44.

PAÏDOUSSIS, M.P. & SUSS, S. 1977 Stability of a cluster of flexible cylinders in bounded axial flow. *Trans. ASME J. Appl. Mech.* **44** (3), 401–408.

PAÏDOUSSIS, M.P., SUSS, S. & PUSTEJOVSKY, M. 1977 Free vibration of clusters of cylinders in liquid-filled channels. *J. Sound Vib.* **55** (3), 443–459.

PETTIGREW, M.J. & TAYLOR, C.E. 2003 Vibration analysis of shell-and-tube heat exchangers: an overview – Part 1: flow, damping, fluidelastic instability. *J. Fluids Struct.* **18** (5), 469–483.

RICCIARDI, G., BELLIZZI, S., COLLARD, B. & COCHELIN, B. 2009 Modelling pressurized water reactor cores in terms of porous media. *J. Fluids Struct.* **25** (1), 112–133.
SCHOUVEILER, L. & ELOY, C. 2009 Coupled flutter of parallel plates. *Phys. Fluids* **21** (8), 081703.

SINGH, K., MICHELIN, S. & DE LANGRE, E. 2012 Energy harvesting from fluid-elastic instabilities of a cylinder. *J. Fluids Struct.* **30**, 159–172.

TAYLOR, G.I. 1952 Analysis of the swimming of long and narrow animals. *Proc. R. Soc. Lond. A* **214** (1117), 158–183.

TRIANTAFYLLOU, G.I. & CHERYSSOSTOMIDIS, C. 1985 Stability of a string in axial flow. *J. Energy Resour. Technol.* **107** (421), 421–425.

WENDEL, K. 1956 Hydrodynamic masses and hydrodynamic moments of inertia. *TMB Translation* **260**.

WU, T.Y. 2001 Mathematical biofluidynamics and mechanophysiology of fish locomotion. *Math. Meth. Appl. Sci.* **24** (17–18), 1464–1541.