# Row of fuel assemblies analyses under seismic loading: modelling and experimental validation

Guillaume Ricciardi, Sergio Bellizzi, Bruno Collard, Bruno Cochelin

► To cite this version:

Guillaume Ricciardi, Sergio Bellizzi, Bruno Collard, Bruno Cochelin. Row of fuel assemblies analyses under seismic loading: modelling and experimental validation. Nuclear Engineering and Design, 2009, 239 (12), pp.2692-2704. 10.1016/j.nucengdes.2009.08.029. hal-00461501

HAL Id: hal-00461501

https://hal.science/hal-00461501v1

Submitted on 28 Mar 2023

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 CC BY 4.0 - Attribution - International License
# Row of fuel assemblies analysis under seismic loading: Modelling and experimental validation

Guillaume Ricciardi<sup>a, b, *</sup>, Sergio Bellizzi<sup>b</sup>, Bruno Collard<sup>a</sup>, Bruno Cochelin<sup>b</sup>

The aim of this study was to develop a numerical model for predicting the impact behaviour at fuel assembly level of a whole reactor core under seismic loading conditions. This model was based on a porous medium approach accounting for the dynamics of both the fluid and structure, which interact. The fluid is studied in the whole reactor core domain and each fuel assembly is modelled in the form of a deformable porous medium with a nonlinear constitutive law. The contact between fuel assemblies is modelled in the form of elastic stops, so that the impact forces can be assessed. Simulations were performed to predict the dynamics of a six fuel assemblies row immersed in stagnant water and the whole apparatus was placed on a shaking table mimicking seismic loading conditions. The maximum values of the impact forces predicted by the model were in good agreement with the experimental data. A Proper Orthogonal Decomposition analysis was performed on the numerical data to analyse the mechanical behaviour of the fluid and structure more closely.

## 1. Introduction

Earthquakes can irreversibly damage nuclear power plants especially in the core, where the fuel assemblies containing enriched uranium have to be particularly resistant. Before building a nuclear power plant, it is necessary to make sure that the core will resist the worst possible earthquake conditions liable to occur at the site. Engineers need special tools for designing and maintaining reactor cores. Therefore when Pressurized Water Reactors (PWRs) are subjected to seismic loading, the spacer grids strike each other, and safety measures are required to prevent the spacer grids from buckling. The model presented here can be used to predict the impact forces generated between grids when the core is subjected to seismic loading. Since direct simulations taking the complex geometry, all the nonlinear phenomena, and the fluid-structure interactions into account would result in too many degrees of freedom, simplified models have been developed. Some of these models are briefly reviewed below.

Fuel assemblies have frequently been modelled in terms of a single beam (Viallet et al., 2003; Rigaureau, 1997) which undergoes

fluid effects via added mass and damping. These simple models can be used to simulate a row of fuel assemblies, and the large number of simulations required due to the statistical character of seismic loading are possible in this case. Some authors have presented multi-beam models for fuel assemblies (Ben Jedida, 1993; Fontaine and Politopoulos, 2000) in which friction laws were used to model the rod/grid connections. These models are in good agreement with "in air" experiments, but the friction problems arising are difficult to solve numerically. Broc et al. (2003) proposed a model with two degrees of freedom for a fuel assembly, taking the fuel assembly coupling into account, and then developed a model for the whole reactor core. Pisapia (2004) developed a single degree of freedom nonlinear empirical model for a fuel assembly giving good agreement with the experimental data. All these models, except for that by Pisapia (2004), take the fluid effects into account in the form of an added mass and a linear added damping term. Païdoussis (1966, 2003) used a more complex expression for the fluid forces acting on a fuel assembly, in which the velocity and the relative direction of the flow with respect to the fuel assembly were accounted for. This model, which was used by Chen and Wambsganss (1970, 1972); Beaud (1997) and Pomîrleanu (2005), has given much better results than the simplified models. However, these complex models do not account for the perturbation flow through the fuel assemblies, although some authors have modelled the fluid flow and the structure using homogenization methods (Jacquelin et al., 1998; Zhang, 1998). In these fluid/structure models, the coupling between fuel assemblies is provided by the fluid flow.

<footnote><sup>a</sup> CEA CADARACHE DEN/DTN/STRI/LHC, 13108 Saint-Paul-Lez-Durance Cedex, France</footnote>

<footnote><sup>b</sup> LMA CNRS, 31 chemin Joseph Aiguier, 13402 Marseille Cedex 20, France</footnote>

<footnote>* Corresponding author at: CEA CADARACHE DEN/DTN/STRI/LHC, 13108 Saint-Paul-Lez-Durance Cedex, France.

E-mail addresses: ricciardi@crans.org (G. Ricciardi), bellizzi@lma.cnrs-mrs.fr (S. Bellizzi), bruno.collard@cea.fr (B. Collard), bruno.cochelin@egim-mrs.fr (B. Cochelin).</footnote>
## Nomenclature

<table><tr><td>N<sub>FA</sub></td><td>number of fuel assemblies</td></tr><tr><td>N<sub>g</sub></td><td>number of spacer grids</td></tr><tr><td>Ω<sub>c</sub></td><td>internal domain of the core</td></tr><tr><td>Ω<sub>s</sub></td><td>domain of the core occupied by the structure</td></tr><tr><td>Ω<sub>f</sub></td><td>domain of the core occupied by the fluid</td></tr><tr><td>Ω<sub>FAi</sub></td><td>domain of the core occupied by the fuel assembly, denoted i</td></tr><tr><td>I<sub>ΩFAi</sub></td><td>indicator function of Ω<sub>FAi</sub></td></tr><tr><td>V<sub>eq</sub></td><td>velocity of the equivalent fluid</td></tr><tr><td>P<sub>eq</sub></td><td>pressure of the equivalent fluid</td></tr><tr><td>u<sub>eqi</sub></td><td>displacement of the fuel assembly, denoted i</td></tr><tr><td>θ<sub>eqi</sub></td><td>rotation of the fuel assembly, denoted i</td></tr><tr><td>u<sub>baffle</sub></td><td>displacement imposed on the baffle</td></tr><tr><td>v<sub>baffle</sub></td><td>velocity imposed on the baffle</td></tr><tr><td>L<sub>cx</sub></td><td>dimension of Ω<sub>c</sub> in the vertical direction</td></tr><tr><td>L<sub>cz</sub></td><td>dimension of Ω<sub>c</sub> in the horizontal direction</td></tr><tr><td>g<sub>p</sub></td><td>gap between two spacer grids</td></tr><tr><td>a</td><td>distance between two rod centers</td></tr><tr><td>β</td><td>porosity</td></tr><tr><td>ρ<sub>eq</sub></td><td>density of the equivalent fluid</td></tr><tr><td>ρ</td><td>density of the fluid</td></tr><tr><td>μ<sub>Teq</sub></td><td>viscosity of the equivalent fluid</td></tr><tr><td>μ</td><td>viscosity of the fluid</td></tr><tr><td>m<sub>FA</sub></td><td>mass per unit length of a fuel assembly</td></tr><tr><td>I<sub>FA</sub></td><td>inertial moment per unit length of a fuel assembly</td></tr><tr><td>m<sub>f</sub></td><td>added mass coefficient</td></tr><tr><td>C<sub>D</sub></td><td>fluid damping coefficient</td></tr><tr><td>S<sub>FA</sub></td><td>area of the cross-section of a fuel assembly</td></tr><tr><td>G<sub>1</sub></td><td>shear stiffness coefficient (linear part)</td></tr><tr><td>G<sub>2</sub></td><td>shear stiffness coefficient (nonlinear part)</td></tr><tr><td>E<sub>1</sub></td><td>bending stiffness coefficient (linear part)</td></tr><tr><td>E<sub>2</sub></td><td>bending stiffness coefficient (nonlinear part)</td></tr><tr><td>μ<sub>G</sub></td><td>structural shear damping</td></tr><tr><td>μ<sub>E</sub></td><td>structural bending damping</td></tr><tr><td>l</td><td>quadratic moment of a fuel assembly</td></tr><tr><td>k<sub>impact</sub></td><td>impact stiffness between two spacer grids</td></tr><tr><td>v<sub>k+1</sub></td><td>discretized fluid velocity unknowns</td></tr><tr><td>P<sub>k+1</sub></td><td>discretized fluid pressure unknowns</td></tr><tr><td>u<sub>k+1</sub></td><td>discretized structure unknowns</td></tr><tr><td>Δt</td><td>time step</td></tr><tr><td>POD</td><td>Proper Orthogonal Decomposition</td></tr><tr><td>Φ<sub>j</sub></td><td>jth Proper Orthogonal Mode (POM)</td></tr><tr><td>α<sub>j</sub>(t)</td><td>Proper Orthogonal Component (POC) of the jth mode</td></tr><tr><td>λ<sub>i</sub></td><td>Proper Orthogonal Value (POV) of the jth mode</td></tr><tr><td>PWR</td><td>Pressurized Water Reactor</td></tr></table>

Ricciardi et al. (2009) developed a model for the whole reactor core using a fluid-structure coupling method and the porous media theory and taking fuel assemblies to resemble porous Timoshenko beams showing nonlinear visco-elastic behaviour. In this model, the beams are subjected to a vertical flow and the fluid equations are space averaged. The model proposed in the present paper is based on the latter model by Ricciardi et al. (2009), and an impact term is added, but since the water in the present study is still, the fluid-structure coupling force is different from that used in Ricciardi et al. (2009). The validity of the predictions obtained with this model was tested under seismic loading conditions. Since engineers need to have the values of the impact forces exerted on the fuel assembly grids to design reactor cores, the simulated data were compared with experimental values. This paper is organised as follows. In Section 2, the porous media modelling process is presented.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Row of N<sub>FA</sub> fuel assemblies.</figcaption></figure>

In Section 3, the numerical method is briefly described. In Section 4, a numerical example involving four fuel assemblies is presented. In Section 5, the results obtained are compared with the experimental data obtained on a row of six fuel assemblies immersed in still water. In the last section, a Proper Orthogonal Decomposition (POD) analysis is performed on numerical results.

## 2. Model formulation

### 2.1. Description of the problem

Let us consider a 2D core comprising N<sub>FA</sub> fuel assemblies arranged in a line and submerged in still water (see Fig. 1). The internal domain of the core is denoted Ω<sub>c</sub>, and its dimensions in the vertical (e<sub>x</sub>) and horizontal (e<sub>z</sub>) directions are L<sub>cx</sub> and L<sub>cz</sub> (i.e. Ω<sub>c</sub> = [0, L<sub>cz</sub>] × [0, L<sub>cx</sub>]). All the fuel assemblies are identical and L<sub>cx</sub> is also its height. A non-penetration condition by the fluid is adopted, on the front and back faces of the domain. A fuel assembly is held in place with N<sub>g</sub> spacer grids and the distance (or gap) in the e<sub>z</sub> direction between two adjacent spacer grids is denoted g<sub>p</sub>. The gap g<sub>p</sub> is almost equal to the distance between any two adjacent rods in a fuel assembly.

The baffle of the core is assumed to be rigid and the fuel assemblies are clamped at their ends to the top and bottom of the baffle. The whole core is made to vibrate, which imposes motion on the rigid frontier. This study focuses on motions simulating earthquake conditions.

### 2.2. Assumptions and methods

The procedure used to model the complex fluid-structure interaction problem under investigation involves making the following assumptions:

H1: The fluid is viscous, incompressible and Newtonian.

H2: Gravity effects are neglected.

H3: The rod section does not undergo any deformation.

H4: The distance between two rods in the same fuel assembly remains constant.

H5: The turbulent kinetic energy is negligible in comparison with the turbulent diffusion processes.

The model presented here was initially based on a vertical flow (Ricciardi et al., 2009), which meant that the gravity effects were negligible in comparison with the inertial and viscous forces. In the
<figure><img src="image_3.png" /><figcaption>Fig. 2. Porous modelling and method.</figcaption></figure>

present study, since there is no vertical flow, the assumption H2 is no longer necessarily valid, and this point will be discussed in Section 2.6. The presence of grids means that the distance between any two rods in a fuel assembly will remain constant until an impact occurs between two grids, in which case the grids may be deformed and the distance between rods may change, whereas the distance between two adjacent rods belonging to two different fuel assemblies will not remain constant as the two fuel assemblies move. This will also be discussed in Section 2.6.

The procedure used here, which is summarized in Fig. 2 is based on a porous medium approach taking $\Omega_c = \Omega_f \cup \Omega_s$ where $\Omega_f$ denotes the domain occupied by the fluid, $\Omega_s$ denotes the domain occupied by the structure, and the porosity $\beta$ is defined by the ratio between the fluid volume and the total volume. This approach yields an equivalent fluid model and an equivalent structure model, both of which are defined in the whole domain $\Omega_c$. Motion equations for the equivalent fluid and the equivalent structure are first drawn up separately. In the case of the fluid part, equations for the total fluid flow through the rod bundle are obtained by applying a local spatial averaging procedure to the Navier-Stokes equations written with an Arbitrary Lagrangian-Eulerian approach. The resulting equivalent fluid is characterized by the unknowns $\mathbf{V}_{eq}$ and $P_{eq}$, which are defined in the whole domain $\Omega_c$. Structure related effects on the fluid are accounted for by a body force, which is also defined in the whole domain. The structure equations are also space averaged, but each fuel assembly is modelled individually as an equivalent structure satisfying a Timoshenko beam model (Timoshenko and Goodier, 1970) with nonlinear visco-elastic behaviour. With each fuel assembly, the unknowns are reduced to the displacement $\mathbf{u}_{eq_i}$ and the rotation $\theta_{eq_i}$, and the unknowns associated with the equivalent structure are obtained by collecting these unknowns with $i = 1, \dots, N_{FA}$. Fluid-related effects on the structure are accounted for by adding a body force, which is defined in the whole domain. Lastly, fluid-related effects on the structure, $\mathbf{F}_{fluid \to structure}$, and structure-related effects on the fluid, $\mathbf{F}_{structure \to fluid}$, have opposite signs and are obtained by modelling the fluid forces acting on a rod subjected to a cross flow.

## 2.3. Coupled model formulation

Let us consider the following subdivision of the whole domain $\Omega_c = \bigcup_{i=1}^{N_{FA}} \Omega_{FA_i}$ where $\forall i, \forall j, i \neq j, \Omega_{FA_i} \cap \Omega_{FA_j} = \emptyset$ and each domain $\Omega_{FA_i}$ includes one fuel assembly (see Fig. 3). Using the results established in Ricciardi et al. (2009) and neglecting the geometric nonlinear effects, the equations of motion governing the equivalent fluid and the equivalent structure mentioned in the previous section are given by the following partial differential

equations in the domain $\Omega_c$:

$$ \rho_{eq} \frac{\partial \mathbf{V}_{eq}}{\partial t} + \rho_{eq} \mathrm{div} \, \mathbf{V}_{eq} \otimes \mathbf{V}_{eq} = -\nabla P_{eq} + \mu_{\mathrm{Teq}} \Delta \mathbf{V}_{eq} + 2 \rho_{eq} \frac{\partial \mathbf{u}_{eq}}{\partial t} \\ \cdot \nabla \mathbf{V}_{eq} - \rho_{eq} \mathbf{V}_{eq} \cdot \nabla \frac{\partial \mathbf{u}_{eq}}{\partial t} + \mathbf{F}_{\mathrm{structure} \to \mathrm{fluid}}, \quad (1) $$

$$ \mathrm{div} \, \mathbf{V}_{eq} = 0, \quad (2) $$

$$ m_{\mathrm{FA}} \frac{\partial^2 \mathbf{u}_{\mathrm{eq}}}{\partial t^2} = \frac{\partial \mathbf{T}}{\partial \mathbf{x}} + \mathbf{F}_{\mathrm{fluid} \to \mathrm{FA}} + \mathbf{F}_{\mathrm{impact}}, \quad (3) $$

$$ I_{\mathrm{FA}} \frac{\partial^2 \theta_{\mathrm{eq}}}{\partial t^2} = \frac{\partial \mathbf{M}}{\partial \mathbf{x}} + \mathbf{e}_{\mathrm{x}} \wedge \mathbf{T}, \quad (4) $$

where $\mathbf{u}_{eq} = \sum_{i=1}^{N_{FA}} u_{eq_i} \mathbf{I}_{\Omega_{FA_i}}$ and $\theta_{eq} = \sum_{i=1}^{N_{FA}} \theta_{eq_i} \mathbf{I}_{\Omega_{FA_i}}$ ($\mathbf{I}_{\Omega_{FA_i}}$ denotes the indicator function of $\Omega_{FA_i}$). $\mathbf{V}_{eq}$, $P_{eq}$, $u_{eq_i}$ and $\theta_{eq_i}$ are the unknowns of the problem.

Eqs. (1) and (2) describe the behaviour of the equivalent fluid. These equations were obtained by applying spatial averaging procedure to the Navier-Stokes equations introducing the equivalent fluid density $\rho_{eq} = \beta\rho$ and the equivalent fluid viscosity $\mu_{\mathrm{Teq}} = \beta\mu$ in terms of the density $\rho$ and the viscosity $\mu$ of the fluid and the porosity $\beta$. The body force $\mathbf{F}_{structure \to fluid}$ is obtained by averaging the stress tensor at the fluid-structure interface. It will be specified in the next section.

Eqs. (3) and (4) describe the behaviour of the equivalent structure. These equations are based on a Timoshenko beam model showing nonlinear visco-elastic behaviour, where $m_{FA}$ denotes the mass per unit length of a fuel assembly, $I_{FA}$ denotes the inertial moment per unit length of a fuel assembly, $\mathbf{M}$ denotes the bending moment and $\mathbf{T}$ denotes the shear force. The term $\mathbf{F}_{fluid \to FA}$ characterizing the fluid-fuel assembly interactions, the term $\mathbf{F}_{impact}$ characterizing the interaction between the fuel assemblies and the terms $\mathbf{T}$ and $\mathbf{M}$ will be described in the next section. Eqs. (1)-(4) are coupled, first because, as classically occurs, Eqs. (3) and (4) include terms depending on the behaviour of the fluid and Eqs. (1) and (2) include terms depending on the behaviour of the structure, but also because all the unknowns are defined in the whole domain $\Omega_c$.

The equations of motion are completely defined based on the boundary conditions. In the case of the structure part, the boundary conditions take into account the clamped state of the fuel assemblies and the imposed motion of the core baffle:

$$ \forall i = 1, \dots, N_{\mathrm{FA}}, \quad \left\{ \begin{array}{l} \mathbf{u}_{eq_i}(0) = \mathbf{u}_{\mathrm{baffle}}, \quad \mathbf{u}_{eq_i}(L_{cx}) = \mathbf{u}_{\mathrm{baffle}}, \\ \theta_{eq_i}(0) = 0, \quad \theta_{eq_i}(L_{cx}) = 0, \end{array} \right. \quad (5) $$

<figure><img src="image_4.png" /><figcaption>Fig. 3. Modelling six fuel assemblies.</figcaption></figure>
where **u**<sub>baffle</sub> denotes the displacement imposed on the core baffle. Note that the imposed motion is not restricted to the horizontal direction **e**<sub>z</sub>.

In the case of the fluid part, a non-penetration condition of the fluid is adopted, and the normal velocity of the fluid is therefore equal to the normal velocity of the baffle:

$$
\forall z, \quad \mathbf{V}_{\text{eq}}(0, z) \cdot \mathbf{e}_x = \mathbf{v}_{\text{baffle}} \cdot \mathbf{e}_x, \quad \mathbf{V}_{\text{eq}}(L_{\text{CX}}, z) \cdot \mathbf{e}_x = \mathbf{v}_{\text{baffle}} \cdot \mathbf{e}_x, \quad (6)
$$

$$
\forall x, \quad \mathbf{V}_{\text{eq}}(x, 0) \cdot \mathbf{e}_z = \mathbf{v}_{\text{baffle}} \cdot \mathbf{e}_z, \quad \mathbf{V}_{\text{eq}}(x, L_{\text{CZ}}) \cdot \mathbf{e}_z = \mathbf{v}_{\text{baffle}} \cdot \mathbf{e}_z, \quad (7)
$$

$$
\text{where } \mathbf{v}_{\text{baffle}} = \partial \mathbf{u}_{\text{baffle}} / \partial t.
$$

2.4. Behaviour modelling

We now briefly describe the models used to account for
the fluid-structure coupling forces, the nonlinear visco-elastic
behaviour of the fuel assemblies and the impacts between fuel
assemblies and between the fuel assemblies and the baffle.

2.4.1. Fluid-structure coupling forces

In line with the local spatial averaging procedure used, the body
force **F**<sub>structure→fluid</sub>, which characterizes the effects of the fluid on
the structure, is modelled empirically using the fluid forces acting
on one rod. It is proposed here to use a different expression for the
fluid-structure coupling forces from that adopted by Ricciardi et
al. (2009), where the flow was vertical, in order to account for the
special case of a cross flow. In line with Morison et al. (1950), the
fluid forces acting on a rod are decomposed into an inertial term
and a drag term as follows:

$$
\mathbf{F}_I = -m_f \frac{\partial}{\partial t} \left( \frac{\partial u_{eqz}}{\partial t} - V_{eqz} \right) \mathbf{e}_z, \quad (8)
$$

$$
\mathbf{F}_D = -C_D \left| \frac{\partial u_{eqz}}{\partial t} - V_{eqz} \right| \left( \frac{\partial u_{eqz}}{\partial t} - V_{eqz} \right) \mathbf{e}_z, \quad (9)
$$

giving,

$$
\mathbf{F}_{\text{structure}\to\text{fluid}} = \frac{1}{a^2} (\mathbf{F}_I + \mathbf{F}_D), \quad (10)
$$

where *a* is the distance between two rod centers. The model is
completely defined by taking the numerical values of the constants
*m*<sub>*f*</sub> and *C*<sub>*D*</sub>, which depend on the geometry, the roughness and the
casing (the gap).

In line with the kinematics of the beam model, the fluid forces
acting on a fuel assembly are the opposite of the body forces acting
on the fluid integrated over the fuel assembly cross-section:

$$
\mathbf{F}_{\text{fluid}\to\text{FA}} = \int_{S_{\text{FA}}} \mathbf{F}_{\text{fluid}\to\text{structure}} \, \mathrm{d}S, \quad (11)
$$

where $S_{FA}$ is the area of the fuel assembly cross-section.

2.4.2. Fuel assembly behaviour

We now focus on describing the mechanical behaviour of a fuel
assembly in terms of an equivalent Timoshenko beam model. In
a fuel assembly, the interactions between spacer grids and rods
give rise to complex behaviour (including contact and friction pro-
cesses) which makes it difficult (as well as being beyond the scope
of this study) to establish the overall behaviour analytically from
the averaged structure equations. It is therefore proposed to model
the overall nonlinear visco-elastic behaviour of the fuel assembly
empirically. Since the results obtained by Pisapia et al. (2003) using
methods based on the quadratic stiffness and damping were in good
agreement with the experimental data, a quadratic law was used
to account for the stiffness, whereas the damping was assumed to
be linear because the structural damping is small in comparison
with the fluid damping (Collard et al., 2003), which means that the

<figure><img src="image_5.png" /><figcaption>Fig. 4. Impact model for four fuel assemblies with three grids.</figcaption></figure>

nonlinear part of the structural damping is negligible. Hence, the
shear force and the bending moment are written as follows:

$$
\mathbf{T} = \left( G_1 S_{\text{FA}} \left( \frac{\partial u_{\text{eqz}}}{\partial x} + \theta_{\text{eqy}} \right) + G_2 S_{\text{FA}} \left| \frac{\partial u_{\text{eqz}}}{\partial x} + \theta_{\text{eqy}} \right| \left( \frac{\partial u_{\text{eqz}}}{\partial x} + \theta_{\text{eqy}} \right) \right) \times \mathbf{e}_{\text{z}} + \mu_{\text{G}} S_{\text{FA}} \frac{\partial}{\partial t} \left( \frac{\partial u_{\text{eqz}}}{\partial x} + \theta_{\text{eqy}} \right) \mathbf{e}_{\text{z}}, \quad (12)
$$

$$
\mathbf{M} = \left( E_1 I \frac{\partial \theta_{eqy}}{\partial x} + E_2 I \left| \frac{\partial \theta_{eqy}}{\partial x} \right| \frac{\partial \theta_{eqy}}{\partial x} + \mu_E I \frac{\partial^2 \theta_{eqy}}{\partial t \partial x} \right) \mathbf{e}_y, \quad (13)
$$

where $G_1$ and $E_1$ denote the linear stiffness coefficients, $G_2$, and $E_2$ denote the nonlinear stiffness coefficients, $\mu_G$ and $\mu_E$ are the (linear) structural damping coefficients, and $l$ is the quadratic moment.

2.4.3. Contact force model

When two fuel assemblies collide, the impact surface is located
at the grids, and buckling of the grids will occur before two rods
strike each other. We therefore model the contact force between
two fuel assemblies or between a fuel assembly and the baffle by
adding elastic stops with gaps in the **e**<sub>z</sub> direction at each grid posi-
tion (Fig. 4). For example, in the case of two adjacent fuel assemblies
denoted *l* and *k*, the impact force at the *i*th grid is given by

$$
\mathbf{F}_{\text{impact}} = \begin{cases} \frac{k_{\text{impact}}}{2} (u_{\text{eqz}}^k(x_{\text{gi}}) - u_{\text{eqz}}^l(x_{\text{gi}}) - g) \mathbf{e}_z & \text{if } u_{\text{eqz}}^k(x_{\text{gi}}) - u_{\text{eqz}}^l(x_{\text{gi}}) - g < 0, \\ 0 & \text{if } u_{\text{eqz}}^k(x_{\text{gi}}) - u_{\text{eqz}}^l(x_{\text{gi}}) - g > 0, \end{cases} \quad (14)
$$

where $x_{\text{gi}}$ is the coordinate in direction **e**<sub>x</sub> of the *i*-th grid. The impact force is assumed to be either a linear elastic force, with spring constant $k_{\text{impact}}$ when grid–grid contact occurs, or zero otherwise. The impact force is also assumed to have only horizontal components, i.e., the grid–grid friction effects are neglected.

The impact force between a fuel assembly and the baffle is mod-
elled in the same way except that the baffle is assumed to be rigid
enough not to be deformed.

2.5. Discussion

The equivalent (or global) model for the behaviour of a 2D core
with fuel assemblies arranged in line presented here is defined
in the whole reactor core domain in that both fluid and struc-
ture are defined in the whole domain. The geometrically complex
physical fluid-structure problem (in each fuel assembly, a large
number of rods are linked by numerous contacts with friction) has
been modelled using a more simplified geometry (with a parallel
epipedic shape). The global model thus avoids the large number
of degrees of freedom involved in directly simulating the complete fluid and structure problem numerically, and both the fluid and structure dynamics of a whole core can be simulated. However, some of the local interactions occurring are lost, such as the vibrations of rods into spacer grids. These effects are modelled by assuming each beam to simply show nonlinear visco-elastic behaviour. However, the equivalent model accounts for the interactions between fuel assemblies, via the fluid effects and the elastic stops inserted between adjacent fuel assemblies. It is therefore reasonable to expect that the effects of an external excitation (such as an earthquake or a plane crash) on the impact forces between fuel assemblies will be accurately simulated.

The mathematical model is completely defined here in terms of the numerical values of all the parameters. The parameters actually required are of two kinds, namely those which are defined by the geometry and the physical characteristics of the materials involved in the problem, and those introduced empirically into the model, which have to be determined on the basis of experimental data. The first set of parameters includes $m_{FA}$, $I_{FA}$ and $S_{FA}$ (for the structure) and $\rho_{eq}$, $\mu_{Teq}$ (for the fluid). The second kind of parameters, which define the behavioural laws of the model, can be divided in two sets, including $G_1$, $G_2$, $E_1$, $E_2$, $\mu_G$, $\mu_E$ and $k_{impact}$ for the first set, and $m_f$ and $C_D$ for the second one. The first set relates only to the behaviour of the structure in air, and the second set defines the fluid-structure coupling forces. The second kind of parameters depend on the geometry of the structure, its roughness, the material, and the casing ($g_p$), but not on the fluid flow conditions or the type of excitation.

## 2.6. Limitations

Experiments were carried out on a row of fuel assemblies submerged in still water, which were subjected to in-line seismic loading. These conditions were obviously not very realistic, but experiments in which vertical flow and 3D seismic loading are exerted on several "scale one" fuel assemblies would require expensive equipment. However, the experimental data used here sufficed to confirm the validity of the model.

The model presented here was originally designed to account for a vertical flow, and its validity was confirmed under these conditions (Ricciardi et al., 2009). The aim of the present study was to test the predictions of the model in the case of seismic loading, and unfortunately the experimental seismic loading set-up available did not include a vertical flow. The fuel assemblies were therefore subjected to a cross flow when the baffles moved in response to the seismic loading, and we had to adapt the model to account for the specificities of a cross flow. The effects of gravity were assumed to be negligible in comparison with the inertial and viscous terms corresponding to vertical flow conditions, but since there was no vertical flow in the present study, there was no obvious reason for neglecting the gravity effects: this assumption may be one of the limitations of this model. We must remember that the final aim of this model was to simulate the mechanical behaviour of a reactor core under operating conditions, which involve a vertical flow.

We assumed that the distance between two rods remains constant (H4), but this is no longer true at the interfaces between fuel assemblies, or when the spacer grids deform under the influence of impact. This is a strong assumption which induces the most important limitation of the model, because changing in the distance between rods may result in fluid forces that are not modelled. Accounting for these fluid forces would be a good way of improving the present model.

Another limiting factor is the fact that the torsional rotations of the fuel assemblies are not accounted for in the model. This torsion has an effect on the impact stiffness, as the spacer grids of two adjacent fuel assemblies will impact each other with a relative

angle so that the contact will occur at the corner of the grids. The impact stiffness will therefore be very different from that involves between two perfectly parallel grids. Moreover, the compressibility of the fluid in the contact region is neither accounted for. However, these limitations will have little effects, Collard and Vallory (2001) observed experimentally that the grids were never parallel. Another way of improving the model could be to include a more complex impact model in order to account for the compressibility of the fluid and the torsion of the fuel assemblies.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Fluid-structure mesh used to discretize one fuel assembly.</figcaption></figure>

## 3. Numerical model

In line with classical procedures, the present numerical model was obtained by performing spatial discretization on the mathematical model (1)–(4) using a 2D finite element method, and temporal discretization on the resulting differential equations.

For the spatial discretization, a mixed "velocity-pressure" element with 9 nodes and 22 degrees of freedom was developed to approximate the fluid part. This fluid element is a quadrilateral defined with eight geometrical nodes, resulting in a quadratic element. The velocity field is interpolated into each element using all nine nodes, whereas the pressure field interpolation involves only the four vertices maintaining the spatial continuity of the pressure field. A companion element was also used to approximate the displacement field characterizing the behaviour of the structure. Recalling that the equivalent structure consists of a set of juxtaposed Timoshenko beams, the classical three-node beam element was used for this purpose. The three nodes are located on the mean line of the beam and the degrees of freedom at each node are the displacement and the rotation. The spatial mesh of a single fuel assembly is shown in Fig. 5 as an example. The finite element mesh covers the whole domain of the fuel assembly. The mesh is constructed beginning with the beam mesh, and then superimposing on each beam element (S1 in Fig. 5 for instance) several fluid elements (F1, F2, F3, F4 and F5 in Fig. 5 for instance) to cover the whole full assembly domain in $e_z$-direction. The three nodes of the beam element coincide with the three middle nodes of the central fluid element (F3 in Fig. 5 for instance). Note that the number of degrees of freedom is larger in the case of the fluid than the structure.

In the time discretization procedure, two different classical schemes were used to discretize the fluid and structure equations: the fluid equations are time discretized with an Uzawa scheme (see Langtangen et al., 2002) and the structure equations are discretized
<table><caption>**Table 1**

Geometrical and material coefficients of the structure and fluid.</caption><thead><tr><th>m<sub>FA</sub></th><th>l<sub>FA</sub></th><th>S<sub>FA</sub></th><th>ρ<sub>eq</sub></th><th>μ<sub>Teq</sub></th></tr></thead><tbody><tr><td>1000 kg/m</td><td>0.83 kg m</td><td>0.04 m<sup>2</sup></td><td>532.44 kg/m<sup>3</sup></td><td>0.1 kg m<sup>2</sup>/s</td></tr></tbody></table>

<table><caption>**Table 2**

Coefficients related to the behaviour of the structure in air.</caption><thead><tr><th>G<sub>1</sub></th><th>G<sub>2</sub></th><th>E<sub>1</sub></th><th>E<sub>2</sub></th><th>μ<sub>G</sub></th><th>μ<sub>E</sub></th><th>k<sub>impact</sub></th></tr></thead><tbody><tr><td>2e7 Pa</td><td>-2e8 Pa</td><td>5e8 Pa</td><td>-6e9 Pa</td><td>2e4 Pa s</td><td>1e6 Pa s</td><td>2e7 N/m</td></tr></tbody></table>

<table><caption>**Table 3**

Coefficients related to the fluid-structure coupling forces.</caption><thead><tr><th>m<sub>t</sub></th><th>C<sub>D</sub></th></tr></thead><tbody><tr><td>0.3 kg/m<sup>3</sup></td><td>2e3 kg m</td></tr></tbody></table>

with a Newmark scheme (see Krenk, 2006). At each time step, the following nonlinear algebraic system has to be solved:

$$ \begin{pmatrix} A_F(\mathbf{v}_{k+1}, \mathbf{u}_{k+1}) & \alpha \Delta t B & H_F(\mathbf{v}_{k+1}) \\ B^T & 0 & 0 \\ H_S(\mathbf{v}_{k+1}) & 0 & A_S(\mathbf{v}_{k+1}, \mathbf{u}_{k+1}) \end{pmatrix} \begin{pmatrix} \mathbf{v}_{k+1} \\ \mathbf{p}_{k+1} \\ \mathbf{u}_{k+1} \end{pmatrix} = \begin{pmatrix} G_k \\ 0 \\ L_k \end{pmatrix}, \quad (15) $$

where $\mathbf{v}_{k+1}$ is the discrete unknown fluid velocity vector, $\mathbf{p}_{k+1}$ is the discrete unknown fluid pressure vector, $\mathbf{u}_{k+1}$ is the discrete unknown structure vector, and $k$ refers to the $k$th time step $\Delta t$. Spatial and temporal discretization give the linear and nonlinear matrices: $A_F(\mathbf{v}_{k+1}, \mathbf{u}_{k+1})$, $B$, $H_F(\mathbf{v}_{k+1})$, $H_S(\mathbf{v}_{k+1})$, $A_S(\mathbf{v}_{k+1}, \mathbf{u}_{k+1})$, $G_k$, and $L_k$. The nonlinear system is solved with Newton's method.

## 4. A numerical example

In the following application, four fuel assemblies ($N_{FA} = 4$) arranged in line and submerged in still water are considered. The characteristic dimensions of the structure under study are: $L_{CX} = 4$ m, $g_p = 2$ mm and $a = 10$ mm. The numerical values of the model parameters are defined in Tables 1-3. The displacement of the baffle is taken to be equal to a sweep sinus (from 0 Hz to 4 Hz) of the form:

$$ \mathbf{u}_{\text{baffle}}(t) = a_0 \sin \left( \frac{\pi(f_s - f_l)}{T_s} t^2 + 2\pi f_l t \right) \mathbf{e}_z, \quad (16) $$

where $a_0 = 1$ mm, $f_l = 0$ Hz, $f_s = 4$ Hz, $T_s = 20$ s ($T_s$ is the time duration of the sweep).

The simulations were performed using the numerical model obtained by performing spatial discretization on the mathematical model with six elements in the vertical direction and one fluid element in the horizontal direction in the case of each fuel assembly. The time step was $\Delta t = 12.5$ ms.

Fig. 6 shows the simulated relative displacements of the fourth spacer grid (from the bottom) of each fuel assembly versus time. All the fuel assemblies move in phase. They contact each other at the same time and the impacts occur when the velocity imposed on the baffle makes a change of sign. The displacements can increase until they reach the limits imposed by the baffle, which is assumed to be rigid. For instance, the relative displacement of the first fuel assembly is bounded with the lower bound $-g_p$ ($-2$ mm) and the upper bound $4g_p$ ($8$ mm), the second assembly with the lower bound $-2g_p$ ($-4$ mm) and the upper bound $3g_p$ ($6$ mm), and so on. It can be noted that some displacement peaks can cross the bounds, which means that the grids modelled by elastic stops are compressed as the result of the impact. Fig. 7 shows the evolution of the impact

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Relative displacements of the fourth grid of each fuel assembly (m) versus time (s).</figcaption></figure>

force between the fourth grid in the first fuel assembly and the baffle versus time. When an impact occurs, the maximum impact force is located at the beginning of the contact. The impact force reaches its maximum value when the natural frequency of the fuel assemblies is excited by the imposed displacement $u_{\text{baffle}}$. It can be seen from Fig. 8 that $u_{\text{baffle}}$ also excites the frequencies from 0 to 4 Hz, and that the FFT of the second fuel assembly shows a peak at its natural frequency.

## 5. Validation on a row of fuel assemblies

In this section, the validity of 2D numerical model is confirmed using experimental data obtained with an experimental apparatus involving “scale one” fuel assemblies.

### 5.1. Experimental apparatus

The experimental apparatus used here consisted of a row of six “scale one” fuel assemblies ($N_{FA} = 6$) immersed in still water (Fig. 9). The six fuel assemblies were confined and the whole apparatus, including the fuel assemblies and confinement baffles, was placed on a shaking table simulating seismic loading. Each fuel

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Reaction force (or impact force) (daN) versus time (s) between the wall and the fourth grid in the first fuel assembly.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Magnitude of the FFT (m) of $u_{\text{baffle}}$ (left), and that of the relative displacements of the fourth grid in the second fuel assembly (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 11.** Experimental apparatus, top view.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Experimental apparatus.</figcaption></figure>

assembly was held in place with eight spacer grids ($N_g = 8$). Displacement sensors were placed at grid numbers 2–7 in the first and last fuel assemblies. In the case of the other fuel assemblies, one displacement sensor was placed on the fourth spacer grid. All the displacement sensors were fixed to the baffle and used to measure the relative displacement between the grids and baffle in the $e_z$ direction. Force sensors were placed on the baffle opposite the spacer grids to measure the impact force in the $e_z$ direction. An acceleration sensor was also fixed to the shaking table, in the $e_z$ direction.

Tests were performed with and without water, under simulated seismic loading conditions at several levels selected between 0.05 g and 0.4 g. The duration of each test was 20 s. Fig. 10 shows the magnitude of the Fourier transform of the acceleration trajectory measured on the shaking table in response to simulated seismic loading with a 0.4 g excitation level. The seismic loading energy was concentrated in the [0.03, 10] Hz frequency range with a maximum at around 3 Hz and a bandwidth of 5 Hz.

In water, tests were performed with two gap sizes of 1.5 mm and 2 mm between the adjacent full assemblies (Fig. 11), whereas in air, tests were performed only with a 2 mm gap.

## 5.2. Coefficient identification

The model includes several coefficients that have to be estimated. As mentioned in Section 2.5, the coefficients were divided

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Magnitude of the Fourier transform of the acceleration (m s⁻², left), velocity (m s⁻¹, middle) and displacement (m, right) of the shaking table versus frequency (Hz), with a 0.4-g excitation level.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Relative displacements of the fourth grid in each fuel assembly (m) versus time (s): comparison between numerical simulation (continuous line) and experimental data (dashed line), in air with a 0.18-g excitation level.</figcaption></figure>

into three sets and each set of coefficients was estimated separately. The numerical values of the first set of coefficients ($m_{FA}$, $I_{FA}$, $S_{FA}$, $\rho_{eq}$, $\mu_{Teq}$) were calculated analytically from the geometrical and material properties of the experimental apparatus.

The numerical values of the second set of coefficients ($G_1$, $G_2$, $E_1$, $E_2$, $\mu_G$, $\mu_E$, $k_{impact}$) were identified from the experimental data. With this set of coefficients including only structural parameters, only data obtained in air experiments were used. The identification

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Relative displacements of the fourth grid in each fuel assembly (m) versus time (s): comparison between numerical simulation (continuous line) and experimental data (dashed line), in water with a 0.28 g excitation level.</figcaption></figure>
process was based on the minimization of the cost function:

$$ \sum_{j=1}^{N_{\text{sensors}}} \sum_{i=1}^{N_s} \frac{\left| d_{\text{exp}}^{j}(t_i) - d_{\text{sim}}^{j}(t_i, G_1, G_2, E_1, E_2, \mu_G, \mu_E, k_{\text{impact}}) \right|}{\left| d_{\text{exp}}^{j}(t_i) \right|} \quad (17) $$

where $d_{\text{sim}}^j$ denotes the simulated displacement at the location of the sensor j, $d_{\text{exp}}^j$ denotes the experimental displacement measured by the sensor j, $N_{\text{sensors}}$ denotes the number of displacement sensors and $N_s$ denotes the number of samples. The estimation procedure was carried out using the data obtained at 0.18 g excitation level (comparisons between simulation and experiments are shown in Fig. 12). This set of coefficients (obtained by optimizing the air data) was used to perform the water simulations.

The last set of coefficients ($m_f, C_D$) characterizing only the fluid-structure coupling forces was first predetermined by performing pressure measurements on rod bundles subjected to fluid flow, and then empirically adjusted to the 0.4 g excitation data obtained with the gap $g_p = 2$ mm.

For confidential reasons, the numerical values of the coefficients used in the following simulations are not disclosed in this paper.

## 5.3. Some validation results

All the simulations were performed with six structural elements in the vertical direction and one fluid element in the horizontal direction with each fuel assembly, giving a total number of 36 fluid elements and 36 beam elements (6 for each fuel assembly). The time step was $\Delta t = 12.5$ ms. Recalling that the model parameters were calculated on the basis of the displacement data obtained in the air configuration and the water configuration with a gap $g_p$ equal to 2 mm and an excitation level equal to 0.4 g, the displacement data and the force data obtained at low excitation levels were used to check the validity of the model. Here we compare only the experimental and simulation data in the water configuration.

Fig. 13 shows the measured and simulated relative displacements of the fourth spacer grid of each fuel assembly with a gap $g_p = 2$ mm and an excitation level equal to 0.28 g. Good agreements was found to exist between the measurements and the numerical results. More specifically, the numerical model was more accurate in the case of the fuel assemblies located at the center of the row (full assembly numbers 3 and 4) than with those located at the ends. The simulation displacement was always “in phase” with the measured displacement and differences are clearly visible on extrema values in the case of fuel assembly numbers 1 and 6, possibility due to the limitation introduced by assumption H4, as discussed above in Section 2.6.

Fig. 14 shows the measured and simulated impact force of the fourth spacer grid in the first fuel assembly with a gap $g_p$ equal to

<figure><img src="image_15.png" /><figcaption>Fig. 14. Impact force (daN) between the casing and the fourth fuel assemblies versus time in water with a 2-mm gap and with a 0.28 g excitation level: numerical simulation (continuous line) and experimental data (dashed line).</figcaption></figure>

2 mm and an excitation level equal to 0.28 g. The numerical model did not predict the evolution of the impact force between the fuel assemblies and the baffle as accurately as the evolution of the relative displacement between full assemblies. Some small impacts were not detected. The large amplitude impact forces were predicted by the numerical model, with a small lag time. We focus now on the ability of the numerical model to predict the maximum impact force occurring between the baffle and the fuel assemblies. The maximum impact force predicted by the numerical simulations is compared in Fig. 15 with the values measured on the experimental apparatus. First of all, as observed experimentally, the maximum impact force simulated increased with the excitation level. In addition, with the gap size $g_p = 2$ mm, the model predictions were satisfactory at all the excitation levels tested (see Fig. 15(a)). Note that this gap configuration was used to estimate the numerical values of the model parameters. In the case of the gap $g_p = 1.5$ mm, the maximum impact force was always underestimated by the model (see Fig. 15(b)), which gives non-conservative estimates and is therefore not suitable from the regulatory standpoint. This gap dependency was not shown by the model, possibly due to the limitations of the model discussed in Section 2.6. The model could be improved by adding an impact term accounting for the compressibility of the fluid and the effects of fuel assembly torsion. However, it is also possible to change the coefficients in order to make the numerical results fit the experimental data at $g_p = 1.5$ mm. We have to remember that since the coefficients were identified with the $g_p = 2$ mm data, from the regulatory standpoint, the same procedure should be used with the $g_p = 1.5$ mm data,

<figure><img src="image_16.png" /><figcaption>Fig. 15. Maximum impact force between the casing and the fourth fuel assemblies: comparison between numerical simulation and experimental data at various excitation levels, in water, (a) with a 2-mm gap and (b) with a 1.5-mm gap.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 16. POM in the fluid-structure field, from left to right first, second, third and fourth mode, from top to bottom 0.1 g, 0.28 g and 0.4 g excitation levels.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. POM in the fluid field, from left to right first, second, third and fourth mode, from top to bottom 0.1 g, 0.28 g and 0.4 g excitation levels.</figcaption></figure>

or data obtained with any gap size, to obtain accurate numerical simulations.

As we have seen, the present model predicts the fuel assembly dynamics reasonably accurately, but it is necessary to perform experiments to obtain the coefficients used in the model. However the aim of the model is not to simulate six fuel assemblies but a whole core. The coefficients can therefore be obtained by performing low cost experiments on a few fuel assemblies with an appropriate configuration (gap), and using the data obtained to simulate other fuel assemblies with various operating conditions. For instance, various configurations combining end of life and

beginning of life fuel assemblies can be analysed with this model under various excitation conditions.

## 6. POD analysis of the dynamics of the fuel assembly row model

In this section, a numerical study on the dynamics behaviour of the fuel assembly row model is presented. The tool used here was the Proper Orthogonal Decomposition (POD) (or Karhunen-Loève decomposition) method, which provides an optimum basis for modelling fields. It has been used successfully to analyse com-
<table><caption>**Table 4**
Energy captured by POMs at various excitation levels.</caption><thead><tr><th rowspan="2">Excitation level</th><th colspan="3">Fluid-structure field</th><th colspan="3">Fluid field</th><th colspan="3">Structure field</th></tr><tr><th>0.1 g</th><th>0.28 g</th><th>0.4 g</th><th>0.1 g</th><th>0.28 g</th><th>0.4 g</th><th>0.1 g</th><th>0.28 g</th><th>0.4 g</th></tr></thead><tbody><tr><td>λ<sub>1</sub></td><td>92.05%</td><td>63.75%</td><td>46.35%</td><td>92.27%</td><td>63.86%</td><td>46.41%</td><td>99.07%</td><td>96.39%</td><td>95.13%</td></tr><tr><td>λ<sub>2</sub></td><td>5.96%</td><td>31.47%</td><td>44.99%</td><td>5.98%</td><td>31.53%</td><td>45.04%</td><td>0.44%</td><td>1.68%</td><td>2.60%</td></tr><tr><td>λ<sub>3</sub></td><td>1.36%</td><td>1.41%</td><td>4.36%</td><td>1.24%</td><td>1.41%</td><td>4.37%</td><td>0.02%</td><td>1.37%</td><td>1.60%</td></tr><tr><td>λ<sub>4</sub></td><td>0.21%</td><td>0.83%</td><td>1.02%</td><td>0.14%</td><td>0.82%</td><td>1.02%</td><td>0.01%</td><td>0.23%</td><td>0.25%</td></tr><tr><td>Σ<sup>4</sup><sub>i=1</sub> λ<sub>i</sub></td><td>99.58%</td><td>97.46%</td><td>96.72%</td><td>99.63%</td><td>97.62%</td><td>96.84%</td><td>99.54%</td><td>99.67%</td><td>99.58%</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Fig. 18.** POM in the structure field, from left to right first, second, third and fourth mode, from top to bottom 0.1 g, 0.28 g and 0.4 g excitation levels.</figcaption></figure>

plex spatial-temporal fields (Païdoussis et al., 2005) and for model reduction purposes (Graham and Kevrekidis, 1996). The POD is used here to extract spatial information from a set of simulated fluid and structure fields. A detailed description of the method can be found in Kerschen et al. (2005).

In this study, the POD method is applied to the spatially discretized field $W(t) = (w(x_1, t)^T, \dots, w(x_{N_n}, t)^T)^T$, where $N_n$ denotes the number of nodes. Starting with the time-averaged correlation matrix:

$$ \mathbf{R} = \frac{1}{N_s} \sum_{i=1}^{N_s} \bar{\mathbf{W}}(t_i) \bar{\mathbf{W}}(t_i)^T, \qquad (18) $$

where $\bar{\mathbf{W}}(t_i) = \mathbf{W}(t_i) - \sum_{j=1}^{N_s} \bar{\mathbf{W}}(t_j)$, and solving the eigenproblem $\mathbf{R}\mathbf{\Phi} = \lambda\mathbf{\Phi}$, the field $\mathbf{W}$ can be decomposed into

$$ \bar{\mathbf{W}}(t) = \sum_{j} \alpha_{j}(t) \mathbf{\Phi}_{j}, \qquad (19) $$

where the eigenvectors $\mathbf{\Phi}_j$ have been numbered based on the associated eigenvalues $\lambda_j$ in decreasing order. Classically, $\mathbf{\Phi}_j$ denotes the jth Proper Orthogonal Mode (POM), $\alpha_j(t)$ the associated Proper Orthogonal Component (POC) and $\lambda_j$ the associated Proper Orthogonal Value (POV). Due to the symmetry of $\mathbf{R}$, the POMs are orthogonal and the POVs (with $\lambda_j > 0$) characterize the energy of the modal component and satisfy $||\mathbf{W}||^2 = \sum_j \lambda_j$. It should be

noted that POMs have no physical significance apart from the fact that they give a good idea of the spatial behaviour of the system. However, after making some assumptions about the system (the linearity, proportional viscous damping, etc.), the POMs have been found to converge with the physical modes of the system (Bellizzi and Sampaio, 2006).

All the results presented in this section were obtained using the mathematical model described and assessed in Section 5.

The POD analysis is applied here to three field configurations: the complete fluid-structure field (i.e. $\mathbf{w}(x, t) = (v(x, t)^T, u(x, t)^T)^T$), the fluid part of the field (i.e. $\mathbf{w}(x, t) = v(x, t)$) and the structure part of the field (i.e. $\mathbf{w}(x, t) = u(x, t)$). Fig. 16 (respectively Figs. 17 and 18) shows the first four POMs obtained in each field configuration at three excitation level values (0.1 g, 0.28 g, 0.4 g). The corresponding POVs are given in Table 4. Fig. 19 shows the evolution of the first four POCs versus the frequency at the highest excitation level tested.

First it is worth noting that in the configurations including structural variables, the structural components of the POMs differ significantly from the linear normal modes of the underlying linear structural system: the set of resonance frequencies obtained for the underlying linear structural system coincides with the set of resonance frequencies of a single beam. The order of multiplicity of each resonance frequency is equal to the number of fuel assemblies (6 in case of the present experimental apparatus), and at each natural resonance frequency, the normal modes are deduced from the
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Magnitude of the FFT of the first four POCs (from left to right $\alpha_1(t)$, $\alpha_2(t)$, $\alpha_3(t)$, $\alpha_4(t)$) versus frequency (Hz), in the fluid-structure field (top), the fluid field (middle), and the structure field (bottom), at a high excitation level (0.4g).</figcaption></figure>

normal mode of the single beam combined with five beams at rest. Hence, the structural components of the POMs (see Figs. 16 and 18) result from the bent shapes of the first and third natural modes of a single beam. The second natural mode is not observed due to the symmetrical property of the system. With a given POM, the bent shapes can be all in phase (see for example the first line in Fig. 16 or the first column in Fig. 18). The POMs with in phase bent shapes are associated with motions where simultaneous large amplitude impacts between fuel assemblies, and/or between the baffle and fuel assemblies can occur.

The structural components (Fig. 16) of the POMs obtained from the fluid-structure field differ from the POMs obtained from the structure field only(Fig. 18). Conversely, the fluid components of the POMs obtained from the fluid-structure field are very similar to those obtained from the fluid field (Fig. 17). The same conclusion can be drawn from a POC analysis (see Fig. 19). This finding is not surprising, because of the difference between the energy levels of the fluid field and the structure field, and it shows that the spatial coherence of the fluid field depends strongly on the skeleton of the structure. The energy distributions over the POCs are similar in both

the fluid field and the fluid-structure field, but they differ considerably in the structure field (see Table 4), where the first mode is always predominated at all the excitation levels tested (where at least 95% of the energy is captured). This was not so in the case of the fluid field, where the percentage energy corresponding to the first mode decreased as the excitation level increased and the second mode captured at least 45% of the energy at high excitation levels. Lastly, the modelling power of the POM basis was investigated. Two bases were tested. The first basis corresponds to the POMs obtained from the fluid-structure field, whereas the second basis was a combination between the POMs obtained from the fluid field and those obtained from the structure field. The first basis was found to give satisfactory fluid variables with a small number of POMs. To obtain a similar level of accuracy in the structure variables, however a large number of POMs are required (see Fig. 20). The second basis gives satisfactory fluid and structure variables with a small number of POMs for (see Fig. 21). This analysis shows that for order reduced modelling purposes, the decomposition should be carried out on the fluid field and the structure field separately in order to obtain satisfactory results with a small number of modal components.

<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Displacement (m) of the fourth grid in the first fuel assembly versus time (s) at a 0.4g excitation level; model with 10, 50 and 100 POMs obtained from the fluid-structure field.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 21. Displacement (m) of the fourth grid in the first fuel assembly versus time (s) at a 0.4 g excitation level; model with 1, 3 and 5 POM obtained from the structure field.</figcaption></figure>

## 7. Conclusion

It was proposed here to model a nuclear reactor core based on a porous medium approach accounting for the dynamics of both the fluid and structure, and including the nonlinear behaviour of the structure. This model can be used to simulate the dynamics behaviour of a core subjected to seismic loading, and therefore to estimate the impact forces generated between the grids of the fuel assemblies. Numerical simulations were performed and compared with experimental data obtained on a row of six "scale one" fuel assemblies immersed in water, in response to a 0.4 g seismic excitation. The displacements and the impact forces predicted by the numerical model were in good agreement with the experimental data. A POD analysis was also performed on the numerical results obtained. This analysis shows that for order reduced modelling purposes, the decomposition should be carried out on the fluid field and the structure field separately in order to obtain satisfactory results with a small number of modal components.

## References

* Beaud, F., (ICONE5-2290) 1997. An analytical model for the prediction of fluidelastic forces in a rod bundle subjected to axial flow: theory, experimental validation and application to PWR fuel assemblies. In: 5th International Conference on Nuclear Engineering Nice, France.

* Bellizzi, S., Sampaio, R., 2006. POMs analysis obtained from Karhunen-Loève expansion for randomly vibrating systems. Journal of Sound and Vibration 297, 774-793.

* Ben Jedida, A., 1993. Etude du comportement sismique d'un coeur de réacteur Nucléaire de type REP. Thèse de doctorat. Université Paris VI.

* Broc, D., Queval, J.C., Viallet, E., 2003. Seismic behaviour of a PWR reactor core: fluid structure interaction effects. In: Transactions of the 17th International Conference on Structural Mechanics in Reactor Technology (SMIRT 17) Prague, Czech Republic.

* Chen, S.S., Wambsganss, M.W., 1972. Parallel-flow induced vibrations of fuel rods. Nuclear Engineering and Design 18, 253-278.

* Chen, S.S., Wambsganss, M.W., 1970. Response of a flexible rod to near-field flow. In: Noise Proceedings of the Conference on Flow-induced Vibrations in Reactors System Components, ANL-7685.

* Collard, B., Pisapia, S., Bellizzi, S., Witters, F., 2003. PWR fuel assembly modal testing and analysis. In: Symposium of Flow-induced Vibration ASME PVP Conference, Cleveland, OH, USA.

* Collard, B., Vallory, J., 2001. Impact forces on a core shroud of an exited PWR fuel assembly. In: 9th International Conference on Nuclear Engineering Nice Acropolis, France.

* Fontaine, B., Politopoulos, I., 2000. A non linear model for the PWR fuel assembly seismic analysis. Nuclear Engineering and Design 195, 321-329.

* Graham, M.D., Kevrekidis, I.G., 1996. Alternative approaches to the Karhunen-Love decomposition for model reduction and data analysis. Computers & Chemical Engineering 20 (5), 495-506.

* Jacquelin, E., Lainé, J.P., Trollat, C., Jézéquel, L., 1998. Modelling the behaviour of a PWR core by a homogenization technique. Computer Methods in Applied Mechanics and Engineering 155, 1-13.

* Kerschen, G., Golinval, J.C., Vakakis, A., Bergman, L., 2005. The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview, Nonlinear Dynamics. Special Issue on Reduced Order Models: Methods and Applications 41, 141-170.

* Krenk, S., 2006. Energy conservation in Newmark based time integration algorithms. Computer Methods in Applied Mechanics and Engineering 195, 6110-6124.

* Langtangen, H.P., Mardal, K.A., Winther, R., 2002. Numerical methods for incompressible viscous flow. Advances in Water Resources 25, 1125-1146.

* Morison, J.R., O'Brien, M.P., Johnson, J.W., Schaaf, S.A., 1950. The forces exerted by surface waves on piles. Petroleum Transactions, AIME 189, 149-157.

* Païdoussis, M.P., Sarkar, A., Semler, C., 2005. A horizontal fluid-conveying cantilever: spatial coherent structures, beam modes and jumps in stability diagram. Journal of Sound and Vibration 280, 141-157.

* Païdoussis, M.P., 2003. Fluid-Structure Interactions: Slender Structures and Axial Flow, vol. 2. Elsevier Academic Press, London.

* Païdoussis, M.P., 1966. Dynamics of flexible slender cylinders in axial flow. Part 1. Theory. Journal of Fluids Mechanics 26, 717-736.

* Pisapia, S., 2004. Etude du comportement vibratoire non-linéaire d'un assemblage combustible de réacteur à eau pressurisée. Thèse de doctorat de l'Université de la Méditerranée Aix-Marseille.

* Pisapia, S., Collard, B., Bellizzi, S., Mori, V., 2003. Modal testing and identification of a PWR fuel assembly. In: Transactions of the 17th International Conference on Structural Mechanics in Reactor Technology (SMIRT 17) Prague, Czech Republic.

* Pomirleanu, R.O., 2005. Spectral response to harmonic excitation of rods in a confined nuclear fuel mini-bundle. In: Proceedings of the PVP 2005: ASME Pressure Vessels and Piping Division Conference, Denver Colorado Paper 71486.

* Ricciardi, G., Bellizzi, S., Collard, B., Cochelin, B., January 2009. Modelling Pressurized Water Reactor cores in terms of porous media. Journal of Fluids and Structures 25 (1), 112-133.

* Rigaudeau, J., (ICONE5-2568) 1997. Grid modelling and strength criterion in the lateral response of PWR fuel assemblies under accident conditions. In: 5th International Conference on Nuclear Engineering Nice, France.

* Timoshenko, S.P., Goodier, J.N., 1970. Theory of Elasticity, third ed. McGraw-Hill, New York.

* Viallet, E., Bolsee, G., Ladouceur, B., Goubin, T., Rigaudeau, J., 2003. Validation of PWR core seismic models with shaking table tests on interacting scale fuel assemblies. In: Transactions of the 17th International Conference on Structural Mechanics in Reactor Technology (SMIRT 17) Prague, Czech Republic.

* Zhang, R.J., 1998. Structural homogenized analysis for a nuclear reactor core. Nuclear Engineering and Design 183, 151-156.