# Seismic response of cylinder assemblies in axial flow

Roberto Capanna1<EMAIL_ADDRESS>Guillaume Ricciardi1  Emmanuelle Sarrouy2
Christophe Eloy3  1 CEA, DES, IRESNE, Department of Nuclear Technology,
Cadarache, 13108, Saint-Paul-Lez-Durance, France 2Aix Marseille Univ, CNRS,
Centrale Marseille, LMA UMR 7031, Marseille, France 3Aix Marseille Univ, CNRS,
Centrale Marseille, IRPHE, Marseille, France

###### Abstract

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, we
first use the proposed model to compute the response of a single cylinder in
axial flow. We then implement a multiple cylinder geometry made of 4 fuel
assemblies, each made of $8\times 8$ cylinders, corresponding to an
experimental facility available at CEA. The comparison between numerical
results and experiments show good agreement. The model can correctly predict
the added mass. It can also qualitatively capture 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.

###### keywords:

##  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
effective.

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 meters
high, has a cross section of about $20\times 20\,$cm2 and weights about
$800\,$kg. The fuel rods, which contain the pellets of uranium, have a
diameter of about $1\,$cm for $4\,$m in height; the space between each fuel
rod 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.

The fuel assemblies are cooled down by a strong axial flow. This water flow is
upwards, with velocities of about $5\,$m s-1 at $150\,$ bar and
$310\,^{\circ}$C. The flow regime is fully turbulent with a Reynolds number
based on the rod diameter of $\Rey\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 & Wamnsganss, 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 (Lamb, 1895; Du Buat, 1786). The added mass $M_{a}$ of a non-
deformable body moving at a 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}=\frac{1}{2}M_{a}U^{2}$ or equivalently

$M_{a}=\frac{\rho}{U^{2}}\int\|\textbf{V}\|^{2}\,dv,$ (1)

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

While the added mass is mainly due 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 widely used (e.g., Païdoussis, 1966a;
Triantafyllou & Cheryssostomidis, 1985; Gosselin & De Langre, 2011; Singh et
al., 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 et al.,
1950), can correctly describe experimental observations providing drag and
inertial coefficients are empirically adjusted (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

$F(X,T)=-M_{a}\left(\partial_{T}+U\partial_{X}\right)^{2}W,$ (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^{2}_{T}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 Clasen (1972)
and Chen & Wamnsganss (1972) on added mass in confined geometries.

The stability analysis of a single cylinder in axial flow has been later
generalized to clusters of cylinders (Païdoussis, 1973; Païdoussis et al.,
1977; Païdoussis & Suss, 1977; Païdoussis, 1979). If each cylinder is free to
move independently, the fluid acts as a coupling medium and the motions of the
cylinders are synchronized. Above a critical flow velocity, the coupled
pinned-pinned cylinders loose 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; Schouveiler & Eloy, 2009; Michelin & Smith, 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 et al., 2007, 2008;
Doaré et al., 2011a; Eloy et al., 2010; Doaré et al., 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 PWR.

A different approach was proposed by Ricciardi et al. (2009) using a porous
medium approach. With this approach, it becomes possible to model both the
fluid and the structure dynamics of a whole core. Some local information are
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 physical mechanisms;
and (2) linear 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 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
1(a)). 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 same origin (Figure 1(a)).

(a) Real space.

(b) Fourier space.

Figure 1: Representation of confined pinned-pinned cylinder deformed under
axial flow in real space and in Fourier space.

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 1(a)).

We use the linearised Euler-Bernoulli beam equation (Bauchau & Craig, 2009) to
describe the dynamics of the cylinder deflection

$M_{s}\partial^{2}_{T}W+B\partial^{4}_{Z}W=F_{Y},$ (3)

where $M_{s}$ is the mass of the cylinder per unit length, $B=EI$ is the
bending rigidity ($E$ being the Young’s modulus and $I=\frac{1}{4}\pi A^{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 are mainly
originating 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)\,Ad\theta,$ (4)

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,

$\displaystyle 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}}.$

In dimensionless form, the linearised Euler-Bernoulli equation (3) becomes

$m\,\partial^{2}_{t}w+b\,\partial^{4}_{z}w=f_{y},$ (5)

where

$\displaystyle m=\frac{M_{s}}{\rho A^{2}},\qquad b=\frac{B}{\rho
U^{2}A^{4}},\qquad f_{y}=\frac{F_{Y}}{\rho U^{2}A}.$

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

$\textbf{V}=U\textbf{e}_{z}+\nabla\Phi,$ (6)

where $\textbf{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)

$\displaystyle\Delta{\phi}$ $\displaystyle=$ $\displaystyle 0,$ (7)
$\displaystyle\left.\frac{\partial{\phi}}{\partial{x}}\right|_{|{x}|={l}_{x}}$
$\displaystyle=$ $\displaystyle 0,$ (8)
$\displaystyle\left.\frac{\partial{\phi}}{\partial{y}}\right|_{|{y}|={l}_{y}}$
$\displaystyle=$ $\displaystyle 0,$ (9)
$\displaystyle\left.\frac{\partial{\phi}}{\partial{n}}\right|_{{r}=1}$
$\displaystyle=$
$\displaystyle-(\partial_{{t}}+\partial_{{z}}){w}(z,t)\sin\theta,$ (10)

where (8) and (9) come from the impermeability of the walls and (10) from the
impermeability on the cylinder wall.

Using the linearised unsteady Bernoulli equation, the dimensionless pressure
field can be linked to $\phi$ (Capanna, 2018)

${p}(x,y,z,t)=-(\partial_{{t}}+\partial_{{z}}){\phi}.$ (11)

Using this relation and applying the operator
$(\partial_{{t}}+\partial_{{z}})$ to the system (7–10) above yields

$\displaystyle\Delta{p}$ $\displaystyle=$ $\displaystyle 0,$ (12)
$\displaystyle\left.\frac{\partial{p}}{\partial{x}}\right|_{|{x}|={l}_{x}}$
$\displaystyle=$ $\displaystyle 0,$ (13)
$\displaystyle\left.\frac{\partial{p}}{\partial{y}}\right|_{|{y}|={l}_{y}}$
$\displaystyle=$ $\displaystyle 0,$ (14)
$\displaystyle\left.\frac{\partial{p}}{\partial{n}}\right|_{{r}=1}$
$\displaystyle=$
$\displaystyle(\partial_{{t}}+\partial_{{z}})^{2}{w}({z},{t})\sin\theta.$ (15)

It 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 airfoil theory.

###  2.2 Problem in Fourier space

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

$\hat{p}(x,y,k,t)=\mathcal{F}\left[p(x,y,z,t)\right],$ (16)

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

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

together with the inverse Fourier transform

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

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

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

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

$\displaystyle(\partial_{{x}}^{2}+\partial_{{y}}^{2})\hat{p}$ $\displaystyle=$
$\displaystyle{k}^{2}\hat{p},$ (20)
$\displaystyle\left.\frac{\partial\hat{p}}{\partial{x}}\right|_{|{x}|={l}_{x}}$
$\displaystyle=$ $\displaystyle 0,$ (21)
$\displaystyle\left.\frac{\partial\hat{p}}{\partial{y}}\right|_{|{y}|={l}_{y}}$
$\displaystyle=$ $\displaystyle 0,$ (22)
$\displaystyle\left.\frac{\partial\hat{p}}{\partial{n}}\right|_{{r}=1}$
$\displaystyle=$ $\displaystyle\hat{\gamma}({k},{t})\sin\theta,$ (23)

where $\hat{\gamma}(k,t)$ is the Fourier transform of the impermeability
boundary condition, such that

$\hat{\gamma}=\mathcal{F}\left[(\partial_{{t}}+\partial_{{z}})^{2}{w}\right].$
(24)

Equation (4) 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 (4), 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,$
(25)

where $\hat{f}_{y}$ is the Fourier transform of $f_{y}$. The pressure
$\hat{p}$ being solution of the linear Helmholtz problem (20–23), its solution
is proportional to $\hat{\gamma}(k,t)$. Hence, (25) can be written

$\hat{f}_{y}=-\hat{\mu}(k)\hat{\gamma}(k,t).$ (26)

where the function $\hat{\mu}(k)$ depends on $l_{x}$ and $l_{y}$ through the
boundary conditions (21) and (22).

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

$m\partial^{2}_{t}\hat{w}+k^{4}b\hat{w}=\hat{\mu}(k)\left[-(\partial^{2}_{t}\hat{w}+2ik\partial_{t}\hat{w}-k^{2}\hat{w})\right],$
(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 (20–23) for different values
of $k$, $l_{x}$ and $l_{y}$, and finally use (25) and (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),$ (28)

where $\chi_{j}(z)$ are the beam eigenmodes and $q_{j}(t)=\eta_{j}e^{i\omega
t}$ 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\mbox{with }k_{j}=j\pi/l,$ (29)

with $l=L/A$ 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}.$
(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.

Figure 2: (a) Eigenmodes for the first 3 modes of a pinned-pinned beam and (b)
their Fourier transform.

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 $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 (28), inserting it
into the beam equation (27) in Fourier space, performing an inverse Fourier
transform, and taking the scalar product (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

$\mathsfbi{M}\,\mathbf{\ddot{q}}+\mathsfbi{K}\,\mathbf{q}=-\left(\mathsfbi{M_{a}}\,\mathbf{\ddot{q}}+\mathsfbi{C_{a}}\,\mathbf{\dot{q}}+\mathsfbi{K_{a}}\,\mathbf{q}\right),$
(31)

where $\mathbf{q}(t)=[q_{1}(t),q_{2}(t),\cdots]^{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 $\mathsfbi{M}=m\mathsfbi{I}$, with $\mathsfbi{I}$ the identity
matrix, is the mass matrix, and
$\mathsfbi{K}=b\,\mathsfbi{diag}(\\{k_{j}^{4}\\})$ is the stiffness matrix
(where $\mathsfbi{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 (31) to zero allows to recover the eigenmodes and the eigenfrequencies
$\omega_{j}=\sqrt{bk_{j}^{4}/m}$ of the beam in vacuum.

The matrices $\mathsfbi{M_{a}}$, $\mathsfbi{C_{a}}$, and $\mathsfbi{K_{a}}$
correspond to the added mass, added damping and added stiffness matrices
respectively. Their coefficients can be calculated as

$\displaystyle\left(\mathsfbi{M_{a}}\right)_{ij}$ $\displaystyle=$
$\displaystyle\langle\chi_{i},\mathcal{F}^{-1}\left[\hat{\mu}\hat{\chi}_{j}\right]\rangle=\frac{1}{2\pi}\langle\chi_{i},\mu\star\chi_{j}\rangle,$
(32a) $\displaystyle\left(\mathsfbi{C_{a}}\right)_{ij}$ $\displaystyle=$
$\displaystyle\langle\chi_{i},\mathcal{F}^{-1}\left[2ik\hat{\mu}\hat{\chi}_{j}\right]\rangle=\frac{1}{2\pi}\langle\chi_{i},2\mu\star\chi^{\prime}_{j}\rangle,$
(32b) $\displaystyle\left(\mathsfbi{K_{a}}\right)_{ij}$ $\displaystyle=$
$\displaystyle\langle\chi_{i},\mathcal{F}^{-1}\left[-k^{2}\hat{\mu}\hat{\chi}_{j}\right]\rangle=\frac{1}{2\pi}\langle\chi_{i},\mu\star\chi^{\prime\prime}_{j}\rangle,$
(32c)

where $\chi_{j}^{\prime}$ and $\chi_{j}^{\prime\prime}$ respectively denote
first and second derivatives of $\chi_{j}$,
$\mu=\mathcal{F}^{-1}\left[\hat{\mu}\right]$ and we have used the property
(19) of the convolution product.

We see here that the knowledge of $\hat{\mu}(k)$ or $\mu(z)$ is enough to
compute the matrices $\mathsfbi{M_{a}}$, $\mathsfbi{C_{a}}$, and
$\mathsfbi{K_{a}}$. In addition, these matrices only depend on the geometry of
the problem (i.e. the confinements $l_{x}$ and $l_{y}$) and not on the
forcing. Finally, note that the Helmholtz problem (20) only depends on
$k^{2}$, which means that $\hat{\mu}(k)$ and $\mu(z)$ are real and even
functions. The coefficients of the matrices $\mathsfbi{M_{a}}$,
$\mathsfbi{C_{a}}$, and $\mathsfbi{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
(20)-(23) is then implemented as follows

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

where $\partial\mathcal{C}$ denotes the rod’s border and $\mathcal{D}$ the
fluid domain (see Appendix A).

###  3.1 Meshing

As previously described, the approach for the numerical resolution of the
mathematical problem is to take advantage of the Fourier transform. This
allows us to solve several 2D problems instead of solving a 3D problem.

$x$
---
$y$
---
$l_{x}$
---
$l_{y}$
---

$1$
---

$\theta$
---
$\hat{p}=0$
---
$\displaystyle\frac{\partial\hat{p}}{\partial n}=0$
---
$\displaystyle\frac{\partial\hat{p}}{\partial n}=\sin(\theta)$
---
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$.

In this section, we present numerical solution of the Helmholtz problem
described by Equations (20-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 can 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 confinement size and wave number. Typically, the
size of the meshes must be decreased when the wave number $k$ is increased. In
addition, the refinement of meshes will increase as the confinement size
decreases and when getting closer to the cylinder.

Explored values are $k\in[10^{-5},10^{3}]$, $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\leq 1$ and increases quadratically for $k>1$. FreeFEM++
code is included in Appendix A.

###  3.2 Slender body limit

In the asymptotic limit of long wavelengths ($k\ll 1$) and negligible
confinement ($l_{x},l_{y}\gg 1$), we expect to recover the result of slender
body theory Lighthill (1960b) given by (2), which writes in dimensionless form

$f_{y}=-m_{a}(\partial_{{t}}+\partial_{{z}})^{2}w,$ (33)

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

$\lim_{k\to 0}\hat{\mu}=\pi,\quad\mbox{for }l_{x},l_{y}\gg 1.$ (34)

Figure 4 shows the behaviour of $\hat{\mu}$ as function of the wavenumber $k$
for different ranges of confinement sizes. The added mass does not change
substantially for wavenumbers below $10^{-2}$ or above $10^{2}$, whatever the
confinement. Moreover, for confinements larger than two diameters
($l_{x},l_{y}>5$), the added mass converges towards $\pi$, as expected from
slender body theory (Figure 4c).

Figure 4: Value of $\hat{\mu}$ as function of the wavenumber $k$ for several
confinements ranges. Dots denote computed values and solid line for $k>1$
denotes the asymptotic extrapolation.

###  3.3 Confinement Effects

Figure 5: Confinement effects for small and big wavenumbers.

Figure 4 shows that $\hat{\mu}$ does not evolve much for small values of $k$,
but still depends on the confinement size as soon as the unconfined hypothesis
cannot be made (that is $l_{x}$ or $l_{y}<5$). Two extreme cases are
considered: (1) very small values of wavenumbers ($k=10^{-5}$); and (2) large
wavenumbers ($k=10$).

Figure 5 shows how the added mass depend on both $x$ and $y$ confinement sizes
for the two extreme wavenumbers considered ($k=10^{-5}$ and $k=10$). For
$k=10^{-5}$, the confinement along the $x$-direction has much more influence
than its counterpart along the $y$-direction (Figure 5a). For $k=10$, this 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 only by 15% between
$l_{y}=10$ and $l_{y}=1.25$.

Figure 4 thus shows that, for large wavenumbers, the dependency on the
confinement is weak. Hence, we can assume that $\hat{\mu}$ is independent of
$l_{x}$ and $l_{y}$, for large $k$ values. Based on this approximation, we can
construct a power-law approximation of $\hat{\mu}$ displayed as black line in
Figure 4

$\hat{\mu}\approx 2.5\,k^{-1},\ k\geq 10.$ (35)

For small wavenumbers, although the cylinder can only move in the $y$
direction, it is the confinement along $x$ that has the strongest influence.
To better understand this paradoxical result, we plot the velocity field in
the Fourier space (Figure 6). 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 showed in Figure 6a.
For large wavenumbers however, no such flow is observed (Figure 6b). By virtue
of mass conservation, when the confinement 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.

Figure 6: Velocity and pressure fields in Fourier space for both small and
large wavenumbers.

Finally, for all the wavenumbers that have been considered and for all
confinements, the added mass $\hat{\mu}$ appears to increase as the
confinement size decreases (in both directions). This observation is in
agreement with the literature works on channel flows made by Chen & Wamnsganss
(1972); Chen (1985); 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 transform are key to directly
relate the cylinder displacement and the pressure resulting force. This
approach leads to fast calculations for any kind of geometries. In the next
section, it will be applied to a group of cylinders and compared to
experimental results in the same geometry.

##  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 7). The test section is placed in vertical position with a square
section of 22.5 cm x 22.5 cm and hosts up to 4 fuel assemblies arranged in a 2
x 2 lattice. The length of the fuel assemblies is $L=2.57$ m. Each of them is
constituted of a squared lattice of $8\times 8$ rods, in which there are 60
rods simulating fuel rods made of acrylic and 4 stainless steel guide tubes.
The guide tubes are empty inside, and they are welded to 5 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 of $10.1\times 10.1\,$cm2, and the 60 rods
constituting it have a diameter of $2A=9$ mm with a pin pitch $P/(2A)=1.33$.
The top and bottom of the assembly are rigidly fixed to the guide tubes, and
they are fixed to the test section through the lower support plate and the
upper support plate.

Figure 7: Scheme of the ICARE experimental facility.

A hydraulic actuator is used to force one of the four fuel assemblies in one
direction; it is attached to one of the grids of the fuel assembly through a
screw, 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 8). For more details on the ICARE experiment, please refer
to Capanna et al. (2019).

Figure 8: Scheme of the displacement sensors on the ICARE test section.

###  4.2 Problem geometry and meshing

The ICARE geometry is implemented in the code to allow for comparison with
experiments with small (4 mm) or large (8 mm) confinements (see Figures 9(a),
9(b)). Confinement size denotes the smallest distance between rods border and
casing walls as well as the smallest distance between rods from different
assemblies.

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

Based on §3, the Helmholtz problem to solve is

$\displaystyle\Delta\hat{p}$ $\displaystyle=k^{2}\hat{p},$
$\displaystyle\text{in the domain},$ (36a)
$\displaystyle\left.\frac{\partial\hat{p}}{\partial
n}\right|_{\partial\text{Cas.}}$ $\displaystyle=0,$ $\displaystyle\text{on the
walls},$ (36b) $\displaystyle\left.\frac{\partial\hat{p}}{\partial
n}\right|_{\partial\mathcal{C}^{[1]}}$
$\displaystyle=\hat{\gamma}^{[1]}_{x}(k,t)\cos(\theta),$
$\displaystyle\text{for assembly [1]},$ (36c)
$\displaystyle\left.\frac{\partial\hat{p}}{\partial
n}\right|_{\partial\mathcal{C}^{[n]}}$ $\displaystyle=0,$
$\displaystyle\text{for other assemblies},$ (36d)

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

$\hat{\gamma}^{[1]}_{x}=\mathcal{F}\left[(\partial_{t}+\partial_{z})^{2}w^{[1]}_{x}\right].$
(37)

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

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

(a) 4 assemblies in large confinement.

(b) 4 assemblies in small confinement.

Figure 9: Mesh for the two configurations simulating the ICARE set-up.

In practice, $\hat{\mu}_{x}^{[n]}$ and $\hat{\mu}_{y}^{[n]}$ can be calculated
by assuming $\hat{\gamma}^{[1]}_{x}=1$ and numerically solving the Helmholtz
problem (36) 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 at least into $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 10: 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$.

Figure 10 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 6th-order polynomials, which can then be used to
calculate ${\mu}_{x}^{[n]}(z)$ and ${\mu}_{y}^{[n]}(z)$ with an inverse
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 $\mathsfbi{M}$, damping $\mathsfbi{C}$, stiffness $\mathsfbi{K}$), 2
axes ($x$ and $y$, noted $\circ$) and 4 assemblies noted with the superscript
$[n]$. Similarly to (32a–c), these matrices are expressed from the added
masses ${\mu}_{\circ}^{[n]}(z)$

$\displaystyle\left(\mathsfbi{M}_{\circ}^{[n]}\right)_{ij}$ $\displaystyle=$
$\displaystyle\langle\chi_{i},\mathcal{F}^{-1}\left[\hat{\mu}_{\circ}^{[n]}\hat{\chi}_{j}\right]\rangle=\frac{1}{2\pi}\langle\chi_{i},\mu_{\circ}^{[n]}\star\chi_{j}\rangle,$
(39a) $\displaystyle\left(\mathsfbi{C}_{\circ}^{[n]}\right)_{ij}$
$\displaystyle=$
$\displaystyle\langle\chi_{i},\mathcal{F}^{-1}\left[2ik\hat{\mu}_{\circ}^{[n]}\hat{\chi}_{j}\right]\rangle=\frac{1}{2\pi}\langle\chi_{i},2\mu_{\circ}^{[n]}\star\chi^{\prime}_{j}\rangle,$
(39b) $\displaystyle\left(\mathsfbi{K}_{\circ}^{[n]}\right)_{ij}$
$\displaystyle=$
$\displaystyle\langle\chi_{i},\mathcal{F}^{-1}\left[-k^{2}\hat{\mu}_{\circ}^{[n]}\hat{\chi}_{j}\right]\rangle=\frac{1}{2\pi}\langle\chi_{i},\mu_{\circ}^{[n]}\star\chi^{\prime\prime}_{j}\rangle,$
(39c)

where $\chi_{i}(z)$ is the $i$-th beam eigenmode defined in (29).

Because it is forced at mid-height, we can assume that the assembly [1] moves
mainly on its first beam eigenmode, i.e.
$w_{x}^{[1]}(z,t)=q_{x}^{[1]}(t)\chi_{1}(z)$. We can further assume 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 the assembly [1] is

$\left(m_{s}+m_{x}^{[1]}\right)\ddot{q}_{x}^{[1]}+c_{v}\dot{q}_{x}^{[1]}+\left(k_{s}+k_{x}^{[1]}\right)q_{x}^{[1]}=f_{e}e^{i\omega
t},$ (40)

where $m_{s}=64m$ is the dimensionless mass of one assembly per unit length,
$k_{s}=64bk_{1}^{4}$ its stiffness, $c_{v}$ account for the damping induced by
the viscosity effect, $m_{x}^{[1]}$ and $k_{x}^{[1]}$ note the component
$(1,1)$ of the added matrices $\mathsfbi{M}_{x}^{[1]}$ and
$\mathsfbi{K}_{x}^{[1]}$ and $f_{e}e^{i\omega t}$ is the external forcing
projected onto the first eigenmode with $f_{e}$ a scalar. $m_{x}^{[1]}$ and
$k_{x}^{[1]}$ coefficients represent added mass and stiffness encountered by
an assembly which moves on its first mode in either $x$ or $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 (40) is simply

$q_{x}^{[1]}=h_{x}^{[1]}f_{e}e^{i\omega t},$ (41)

with the transfer function

$h_{x}^{[1]}(\omega)=\frac{1}{-\omega^{2}\left(m_{s}+m_{x}^{[1]}\right)+i\omega
c_{v}+k_{s}+k_{x}^{[1]}}.$ (42)

It is possible to estimate the motion of assemblies [2] to [4] by considering
the superimposition of the two following cases: still rods undergoing fluids
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:

$\left(m_{s}+m_{x}^{[1]}\right)\ddot{q}_{\circ}^{[n]}+c_{v}\dot{q}_{\circ}^{[n]}+\left(k_{s}+k_{x}^{[1]}\right)q_{\circ}^{[n]}=-\left(-\omega^{2}m_{\circ}^{[n]}+k_{\circ}^{[n]}\right)q_{x}^{[1]},$
(43)

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_{\circ}^{[n]}$ on the assembly $[n]$ itself.
The right-hand side represents action of assembly [1] motion via fluid.

From (39c), we see that the added stiffness is
$k_{\circ}^{[n]}=-k_{1}^{2}m_{\circ}^{[n]}$, with $k_{1}=\pi/l$ (since
$\chi_{i}^{\prime\prime}=-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_{\circ}^{[n]}$ can be neglected compared to the added mass term
$-\omega^{2}m_{\circ}^{[n]}$ in (40) and (43). The values of the added masses
$m_{\circ}^{[n]}$ are given in table 1. From this table, we see that the ratio
of added masses, $m_{\circ}^{[n]}/m_{x}^{[1]}$ is at most 21% (for $n>1$).
This justifies to neglect 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 thus of order 4% at most.

Table 1: Added mass matrix coefficients for small and large confinement confinement | $m_{x}^{[1]}$ | $m_{x}^{[2]}$ | $m_{x}^{[3]}$ | $m_{x}^{[4]}$ | $m_{y}^{[1]}$ | $m_{y}^{[2]}$ | $m_{y}^{[3]}$ | $m_{y}^{[4]}$
---|---|---|---|---|---|---|---|---
large | 350 | -22.1 | 72.7 | 18.6 | -17.0 | 20.5 | -20.5 | 26.2
small | 384 | -24.0 | 79.0 | 21.3 | -22.4 | 25.3 | -25.3 | 29.3

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

$\frac{q_{\circ}^{[n]}}{q_{x}^{[1]}}\approx-\frac{-\omega^{2}m_{\circ}^{[n]}}{-\omega^{2}\left(m_{s}+m_{x}^{[1]}\right)+i\omega
c_{v}+k_{s}},$ (44)

which become at the resonant frequency $\omega_{\mathrm{res}}$

$\frac{q_{\circ}^{[n]}}{q_{x}^{[1]}}\approx-i\omega_{\mathrm{res}}\frac{m_{\circ}^{[n]}}{c_{v}},\quad\mbox{with
}-\omega_{\mathrm{res}}^{2}\left(m_{s}+m_{x}^{[1]}\right)+k_{s}=0,\quad$ (45)

note that the latter equation is applicable only for $n\neq$1.

##  5 Comparison with experimental results

This section compares results obtained using the simplified model and results
from an experimental campaign on ICARE facility described in §4.1. Experiments
are performed imposing a sine sweep ranging from 0 to 10 Hz with a constant
amplitude of 1 mm and an axial fluid velocity of 1 m/s when water is present.
To catch the transfer function of a system using a swept sinus excitation, the
sweep velocity needs to be small enough to avoid any transient phenomena and
contamination of different harmonics in the system. For these experiment a
sweep velocity of 0.05 Hz 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^{[1]}_{exp}$ between the displacement
$W^{[1]}_{x}$ along $x$ of the grid with vertical position $Z_{w}$, and the
force $F_{g}$ applied by the actuator on grid with position $Z_{f}$. Since the
data are acquired with a swept sinus excitation, the transfer functions are
calculated using the cross correlation product, and are defined as

$\displaystyle H^{[1]}_{exp}=\left\|\frac{W^{[1]}_{x}\star
W^{[1]}_{x}}{W^{[1]}_{x}\star F_{x}}\right\|$ (46)

where $W^{[1]}_{x}\star W^{1]}_{x}$ is the spectral autocorrelation of the
displacement signal of the $1$-st assembly along the $x$ direction and
$W^{[1]}_{x}\star F_{x}$ is the spectral cross correlation between the
displacement signal of the $1$-st 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 :

$\displaystyle
H^{[1]}_{theo}(\Omega)=\frac{1}{\sqrt{(K_{e}-\Omega^{2}M_{e})^{2}+\Omega^{2}C_{e}^{2}}},$
(47)

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^{[1]}_{theo}$ fit the experimental data $H^{[1]}_{exp}$. This can be easily
achieved by a conjugate gradient optimization method.

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

$\begin{split}M_{t}\int^{L}_{0}\partial^{2}_{T}W^{[1]}_{x}\delta
WdZ+C_{t}\int^{L}_{0}\partial_{T}W^{[1]}_{x}\delta
WdZ+B_{t}\int^{L}_{0}\partial^{2}_{Z}W^{[1]}_{x}\partial^{2}_{Z}\delta
WdZ\qquad\\\ \qquad=\int_{0}^{L}F_{g}\delta_{Z_{f}}\delta W\,dZ,\end{split}$
(48)

with

$M_{t}=M_{s}+M_{f},\quad B_{t}=B+B_{f},$ (49)

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 $j$-th natural mode

$W^{[1]}_{x}(Z,T)=N_{j}(T)X_{j}(Z),$ (50)

with $X_{j}(Z)=\chi_{j}(Z/A)$, one can rewrite the energy of the system :

$M_{t}\int^{L}_{0}X_{j}^{2}dZ\,\partial^{2}_{T}N_{j}+C_{t}\int^{L}_{0}X_{j}^{2}dZ\,\partial_{T}N_{j}+B_{t}\int^{L}_{0}(\partial^{2}_{Z}X_{j})^{2}dZ\,N_{j}=F_{g}X_{j}(Z_{f}).$
(51)

Noting that the experimental signal is

$W^{[1]}_{g}=N_{1}(T)X_{1}(Z_{w}),$ (52)

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^{L}_{0}X_{1}^{2}dZ}M_{e}.$ (53)

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 without water. Finally, one then obtains the results listed in Table
2 where added mass are provided for the whole assembly:

$M^{[1]}_{x}=M_{f}.$ (54)

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

$M^{[1]}_{x}=\rho A^{2}\,m_{x}^{[1]},$ (55)

with $\rho=1000\,$kg m-3 the water density.

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 gives a large increase, more than
double when simulation predict an increase of about 10%. Although the order of
magnitude is respected the differences observed could question the model. To
properly analyse these results, one has to keep in mind that the methodology
used to obtained the experimental results is subjected to uncertainties coming
from the accuracy of measurement, the optimization 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.

Table 2: Values of added mass $M^{[1]}_{x}$ $M^{[1]}_{x}$ |  experiment |  simulation
---|---|---
Large confinement |  4.06 kg/m |  7.09 kg/m
Narrow confinement |  10.36 kg/m |  7.78 kg/m

No experimental identification can be made on added stiffness as it is very
low compared to structure initial stiffness. To give orders of magnitude, the
measured stiffness of a fuel assembly $B_{t}=2300$ N.m2 while considering the
estimation of added stiffness from numerical simulation in the large
confinement case obtained in §4, $k^{[1]}_{x}\rho U^{2}A^{4}/(k_{1}^{4})=4.7$
N.m2. Experimental results are not accurate enough to measure this less than
0.5% added stiffness effect.

###  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 (46) transfer functions are calculated as

$\displaystyle H^{[n]}_{\circ}=\left\|\frac{W^{[n]}_{\circ}\star
W^{[n]}_{\circ}}{W^{[n]}_{\circ}\star W^{[1]}_{x}}\right\|,\
\circ\in\\{x,y\\},$ (56)

with $W^{[n]}$ being 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
displacement ratio given by (45), with the addedd mass coefficients
$m_{\circ}^{[n]}$ given in table 1 and $c_{v}$ and $\omega_{\mathrm{res}}$
estimated from identification of dynamical parameters using (51) as follows

$C_{v}=\frac{X_{1}(Z_{f})X_{1}(Z_{w})}{\int^{L}_{0}X_{1}^{2}dZ}C_{e},\quad
c_{v}=\frac{C_{v}}{\rho
UA},\quad\omega_{\mathrm{res}}=\frac{A}{U}\sqrt{\frac{K_{e}}{M_{e}}}.$ (57)

Experimental values give $c_{v}=69.2$ and $\omega_{\mathrm{res}}=0.12$.

Figure 11 compare graphically the displacement ratios as observed in the
experiments with the ratios predicted by (45) with these values of $c_{v}$ and
$\omega_{\mathrm{res}}$. Experiments show a ratio of about 10 % for most of
the fuel assemblies in the two directions except for the fuel assembly [3]
which shows a coupling about 20 % in the direction of forcing. Although
simulations seems to underestimate the values of coupling, they give a good
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.

Figure 11: Displacement ratios obtained by experimental data (left) and
approximated from numerical simulations (right) for 4 fuel assemblies in large
confinement at first natural frequency.

##  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 equation 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 to drastically reduce computational
time, avoiding to solve fluid equations. As a drawback it should be considered
that the assumption of potential flow leads to the complete negligence of the
fuel grids, which introduce vorticity in the flow.

Calculations with potential flow model are first performed for single cylinder
geometry and results have been compared with reference works in literature.
The model showed perfect agreement with the Slender Body Model and it also
shows good predictions respect to the confinement size. This validation has
proven the model, and thus the mathematical approach used to solve the
equations, 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 reasonably fits 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 which 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 forces in the model. One possible way in
order 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 which has a
defined empirical distribution. In this way the viscous forces are 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.

## Acknowledgement

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

## Declaration of Interests

The authors report no conflict of interest.

## References

  * Bauchau & Craig (2009) Bauchau, O. A. & Craig, J. I. 2009 Structural Analysis: With Applications to Aerospace Structures. Springer.
  * Blevins (2015) Blevins, R.D. 2015 Formulas for Dynamics, Acoustics and Vibration. John Wiley & Sons, Ltd.
  * Brennen (1982) Brennen, C.E. 1982 A review of added mass and fluid inertial forces. Technical Report Naval Civil Engineering Laboratory .
  * Capanna (2018) Capanna, R. 2018 Modelling of fluid structure interaction by potential flow theory in a PWR under seismic excitation. PhD Thesis .
  * Capanna et al. (2019) Capanna, R., Ricciardi, G., Eloy, C. & Sarrouy, E. 2019 Experimental study of fluid structure interaction on fuel assemblies on the icare experimental facility. Nuclear Engineering and Design 352.
  * Chen (1985) Chen, S.S. 1985 Flow-induced vibration of circular cylindrical structures. Argonne National Lab ANL 85-81.
  * Chen & Wamnsganss (1972) Chen, S.S. & Wamnsganss, M.W. 1972 Parallel-flow-induced vibration of fuel rods. Nuclear Engineering and Design 18 (2), 253–278.
  * Clasen (1972) Clasen, P. 1972 Die ’hydrodynamische masse’ eines rohrs in einem rohrbundel. Forshung im IngenieurwesenNuclear Engineering and Design 38 (2), 33–37.
  * De Langre et al. (2007) De Langre, E., Paaïdoussis, M.P., Doaré, O. & Modarres Sadeghi, Y. 2007 Flutter of long flexible cylinders in axial flow. Journal of Fluid Mechanics 571, 371–389.
  * De Ridder et al. (2013) De Ridder, J., Degroote, J., Van Tichelen, K. & Vierendeels, J. 2013 Modal characteristics of a flexible cylinder in turbulent axial flow from numerical simulations. Journal of Fluids and Structures 43, 110–123.
  * De Ridder et al. (2017) 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. Journal of Fluids and Structures 74, 90–110.
  * De Ridder et al. (2015) 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. Journal of Fluids and Structures 55, 139–154.
  * Doaré et al. (2011a) Doaré, O., Mano, D. & Bilbao Ludena, J. C. 2011a Effect of spanwise confinement on flag flutter: experimental measurements. Physics of Fluids, American Institute of Physics 23, 11704.
  * Doaré et al. (2011b) Doaré, O., Sauzade, M. & Eloy, C. 2011b Flutter of an elastic plate in a channel flow: confinement and finite-size effects. Journal of Fluids and Structures 27 (1), 76–88.
  * Du Buat (1786) Du Buat, P.L.G. 1786 Principles d’hydraulique. Imprimerie du Monsieur .
  * Eloy et al. (2010) Eloy, C., Doaré, O., Duchemin, L. & L., Schouveiler 2010 A unified introduction to fluid mechanics of flying and swimming at high reynolds number. Experimental Mechanics 50 (9), 1361–1366.
  * Eloy et al. (2008) Eloy, C., Lagrange, R., Souilliez, C. & Schouveiler, L. 2008 Aeroelastic instability of cantilevered flexible plates in uniform flow. Journal of Fluid Mechanics 611, 97–106.
  * Eloy et al. (2007) Eloy, C., Souillez, C. & Schouveiler, L. 2007 Flutter of a rectangulare plate. Journal of Fluids and Structures 23 (6), 904–919.
  * Gosselin & De Langre (2011) Gosselin, F.P. & De Langre, E. 2011 Drag reduction by reconfiguration of a poroeslatic system. Journal of Fluids and Structures 27 (7), 1111–1123.
  * Guo & Païdoussis (2000) Guo, C.Q. & Païdoussis, M. P. 2000 Stability of rectangular plates with free side-edges in two-dimensional inviscid channel flow. Journal of Applied Mechanics 67 (1), 171–176.
  * Hecht (2012a) Hecht, F. 2012a Freefem++ v.3.19-1. Universite Pierre et Marie Curie. http://www.freefem.org/ff++/ftp/freefem++doc.pdf .
  * Hecht (2012b) Hecht, F. 2012b New development in FreeFem++. Journal of Numerical Mathematics 20 (3-4), 251–265.
  * Lamb (1895) Lamb, H. 1895 Hydrodynamics. Cambridge University Press .
  * Lighthill (1960a) Lighthill, M.J. 1960a Mathematics and Aeronautics. The Journal of the Royal Aeronautical Society 64 (595), 375–394.
  * Lighthill (1960b) Lighthill, M.J. 1960b Note on the swimming of slender fish. Journal of Fluid Mechanics 9 (2), 305–317.
  * Michelin & Smith (2009) Michelin, S. & Smith, S.L. 2009 Linear stability analysis of coupled parallel flexible plates in an axial flow. Journal of Fluids and Structures 25 (7), 1136–1157.
  * Morison et al. (1950) Morison, J.R., Johnson, J.W. & Schaaf, S.A. 1950 The force distribution exerted by surface waves on piles. Journal of Petroleum Technology 2 (5), 149–154.
  * Païdoussis (1966a) Païdoussis, M. P. 1966a Dynamics of slender cylinders in axial flow: Part 1 Theory. Journal of Fluid Mechanics 26 (4), 717–736.
  * Païdoussis (1966b) Païdoussis, M. P. 1966b Dynamics of slender cylinders in axial flow: Part 2 Experiments. Journal of Fluid Mechanics 26 (4), 737–751.
  * Païdoussis (1973) Païdoussis, M. P. 1973 Dynamics of cylindrical structures subjected to axial flow. Journal of Sound and Vibration 29 (3), 365–385.
  * Païdoussis (1979) Païdoussis, M. P. 1979 The dynamics of clusters of flexible cylinders in axial flow: theory and experiments. Journal of Sound and Vibration 65 (3), 391–417.
  * Païdoussis & Pettigrew (1979) Païdoussis, M. P. & Pettigrew, M.J. 1979 Dynamics of flexible cylinders in axisymetrically confined axial flow. Journal of Applied Mechanics 46 (1), 37–44.
  * Païdoussis & Suss (1977) Païdoussis, M. P. & Suss, S. 1977 Stability of a cluster of flexible cylinders in bounded axial flow. Journal of Applied Mechanics 44 (3), 401–408.
  * Païdoussis et al. (1977) Païdoussis, M. P., Suss, S. & Pustejovsky, M. 1977 Free vibration of clusters of cylinders in liquid-filled channels. Journal of Sound and Vibration 55 (3), 443–459.
  * Ricciardi et al. (2009) Ricciardi, G., Bellizzi, S., Collard, B. & B., Cochelin 2009 Modelling Pressurized Water Reactor cores in terms of porous media. Journal of Fluids and Structures 25 (1), 112–133.
  * Schouveiler & Eloy (2009) Schouveiler, L. & Eloy, C. 2009 Coupled flutter of parallel plates. Physics of Fluids 21 (8).
  * Singh et al. (2012) Singh, K., Michelin, S. & De Langre, E. 2012 Energy harvesting from fluid-elastic instabilities of a cylinder. Journal of Fluids and Structures 30, 159–172.
  * Taylor (1952) Taylor, G.I. 1952 Analysis of the swimming of long and narrow animals. Proceedings of the Royal Society 214 (1117), 158–183.
  * Triantafyllou & Cheryssostomidis (1985) Triantafyllou, G.I. & Cheryssostomidis, C. 1985 Stability of a string in axial flow. ASME Journal of Energy Ressources Technology 107 (421), 421–425.
  * Wendel (1956) Wendel, K. 1956 Hydrodynamic masses and hydrodynamic moments of inertia. TMB Translation 260.
  * Wu (2001) Wu, T.Y. 2001 Mathematical biofluiddynamics and mechanophysiology of fish locomotion. Mathematical Methods in the Applied Sciences 24 (17-18), 1464–1541.

##  Appendix A FreeFEM++ code for one cylinder

⬇

1real k = 0.01;

real lx = 1.5;

real ly = 2;

real A = 1.;

6// Geometry

border 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;}

11border 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.);

16int 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 (36) 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. From this,
we can calculate the associated added masses $\hat{\mu}_{x}^{[n,q]}$ and
$\hat{\mu}_{y}^{[n,q]}$

$\displaystyle\hat{\mu}_{x}^{[n,q]}$
$\displaystyle=\int_{\partial\mathcal{C}^{[n,q]}}\hat{p}\cos(\theta)\,d\theta,$
(58a) $\displaystyle\hat{\mu}_{y}^{[n,q]}$
$\displaystyle=\int_{\partial\mathcal{C}^{[n,q]}}\hat{p}\sin(\theta)\,d\theta,$
(58b)

with $\partial\mathcal{C}^{[n,q]}$ the border of the 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 12 and 13 for the 8 mm
confinement.

Figure 12: Added mass along $x$ $\hat{\mu}_{x}^{[n,q]}$ for each rod for four
values of $k$ as labelled. Full circles note positive values, hollow ones
negative values. The confinement is 8 mm. Figure 13: Same as Figure 12 for the
$y$-component of the added mass: $\hat{\mu}_{y}^{[n,q]}$.

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

##  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 6th-
order polynomials. These interpolations can be used to compute the added
masses ${\mu}_{\circ}^{[n]}(z)$ with an inverse Fourier transform (Figure 14).

Figure 14: 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).

The functions ${\mu}_{\circ}^{[n]}(z)$ are needed in (39a–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 14). This means that on a $z$-scale of order $l=L/A\approx 571$,
these functions are almost proportional to a Dirac, such that the convolution
products between $\mu_{\circ}^{[n]}$ and $\chi_{j}$, $\chi_{j}^{\prime}$, or
$\chi_{j}^{\prime\prime}$ appearing in (39a–c) can be approximated by simple
product with $\int_{-\infty}^{\infty}\mu_{\circ}^{[n]}\,dz$.

As a consequence, the added mass matrix $\mathsfbi{M}_{\circ}^{[n]}$ and the
added stiffness matrix $\mathsfbi{K}_{\circ}^{[n]}$ are almost diagonal (since
$\langle\chi_{i},\chi_{j}\rangle=\delta_{ij}$) and the diagonal terms of the
added damping matrix $\mathsfbi{C}_{\circ}^{[n]}$ are almost zero (since
$\langle\chi_{i},\chi_{i}^{\prime}\rangle=0$). This can be verified by
computing the first $3\times 3$ terms of the matrices $\mathsfbi{M}_{x}^{[n]}$
and $\mathsfbi{C}_{x}^{[n]}$ for $n=1$ and $4$ (we chose the assembly $[4]$
because it is the assembly for which the extension of ${\mu}_{x}^{[n]}(z)$ is
widest, see figure 14)

$\displaystyle\mathsfbi{M}_{x}^{[1]}$ $\displaystyle=$
$\displaystyle\left(\begin{tabular}[]{ccc}350&0&0.0182\\\ 0&347&0\\\
0.0182&0&342\end{tabular}\right),$ (59d) $\displaystyle\mathsfbi{M}_{x}^{[4]}$
$\displaystyle=$ $\displaystyle\left(\begin{tabular}[]{ccc}18.6&0&0.0882\\\
0&17.4&0\\\ 0.0882&0&15.4\end{tabular}\right),$ (59h)
$\displaystyle\mathsfbi{C}_{x}^{[1]}$ $\displaystyle=$
$\displaystyle\left(\begin{tabular}[]{ccc}0&-3.26&0\\\ 3.26&0&-5.79\\\
0&5.79&0\end{tabular}\right),$ (59l) $\displaystyle\mathsfbi{C}_{x}^{[4]}$
$\displaystyle=$ $\displaystyle\left(\begin{tabular}[]{ccc}0&-0.167&0\\\
0.167&0&-0.272\\\ 0&0.272&0\end{tabular}\right).$ (59p)

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


can we use the insights from our results to design defenses that are valid in the classical PAC-learning model? As we discussed in Section <ref>, this might be challenging because certifying randomness between entities that are fully classical seems unreachable in our setup. But it is possible that there exists a third model, neither PAC nor quantum PAC, in which it is possible to obtain a result of similar flavor. Moreover, we can turn the question around and ask for a reduction of the problem of certifiable randomness to the adversarial robustness problem. If such a reduction exists then this strongly implies that in order to tackle adversarial robustness we should first solve the problem of certifiable randomness. Finally, our results are valid for the Hellinger distance. It is natural to ask if our reduction can be generalized to other similarity metrics and if such a potential generalization might be advantageous.

§ ACKNOWLEDGMENT

We would like to thank Thomas Vidick, whose superb lecture notes [Vidick, 2021] have aided us substantially in learning about the recent progress in quantum protocols.

[Arjovsky, 2021]
M. Arjovsky.
Out of distribution generalization in machine learning.
CoRR, abs/2103.02667, 2021.
URL <https://arxiv.org/abs/2103.02667>.

[Babai et al., 1990]
L. Babai, L. Fortnow, and C. Lund.
Nondeterministic exponential time has two-prover interactive
In Proceedings of the 31st Annual Symposium on Foundations of
Computer Science, SFCS '90, page 16–25 vol.1, USA, 1990. IEEE Computer
ISBN 081862082X.
URL <https://doi.org/10.1109/FSCS.1990.89520>.

[Biamonte and Love, 2008]
J. D. Biamonte and P. J. Love.
Realizable hamiltonians for universal adiabatic quantum computers.
Phys. Rev. A, 78:0 012352, Jul 2008.
URL <https://link.aps.org/doi/10.1103/PhysRevA.78.012352>.

[Brakerski et al., 2021]
Z. Brakerski, P. Christiano, U. Mahadev, U. Vazirani, and T. Vidick.
A cryptographic test of quantumness and certifiable randomness from a
single quantum device.
J. ACM, 680 (5), Aug. 2021.
ISSN 0004-5411.
URL <https://doi.org/10.1145/3441309>.

[Bshouty and Jackson, 1995]
N. H. Bshouty and J. C. Jackson.
Learning dnf over the uniform distribution using a quantum example
In Proceedings of the Eighth Annual Conference on Computational
Learning Theory, COLT '95, page 118–127, New York, NY, USA, 1995.
Association for Computing Machinery.
ISBN 0897917235.
URL <https://doi.org/10.1145/225298.225312>.

[Bubeck et al., 2018]
S. Bubeck, Y. T. Lee, E. Price, and I. P. Razenshteyn.
Adversarial examples from cryptographic pseudo-random generators.
CoRR, abs/1811.06418, 2018.
URL <http://arxiv.org/abs/1811.06418>.

[Bubeck et al., 2019]
S. Bubeck, Y. T. Lee, E. Price, and I. Razenshteyn.
Adversarial examples from computational constraints.
In K. Chaudhuri and R. Salakhutdinov, editors, Proceedings of
the 36th International Conference on Machine Learning, volume 97 of
Proceedings of Machine Learning Research, pages 831–840, Long Beach,
California, USA, 09–15 Jun 2019. PMLR.
URL <http://proceedings.mlr.press/v97/bubeck19a.html>.

[Chu, 1965]
J. T. Chu.
Optimal decision functions for computer character recognition.
J. ACM, 120 (2):0 213–226, Apr. 1965.
ISSN 0004-5411.
URL <https://doi.org/10.1145/321264.321271>.

[Colbeck, 2009]
R. Colbeck.
Quantum and relativistic protocols for secure multi-party
arXiv: Quantum Physics, 2009.

[Einstein et al., 1935]
A. Einstein, B. Podolsky, and N. Rosen.
Can quantum-mechanical description of physical reality be considered
Phys. Rev., 47:0 777–780, May 1935.
URL <https://link.aps.org/doi/10.1103/PhysRev.47.777>.

[Engstrom et al., 2019]
L. Engstrom, B. Tran, D. Tsipras, L. Schmidt, and A. Madry.
A rotation and a translation suffice: Fooling CNNs with simple
transformations, 2019.
URL <https://openreview.net/forum?id=BJfvknCqFQ>.

[Fitzsimons et al., 2018]
J. F. Fitzsimons, M. Hajdu ššek, and T. Morimae.
Post hoc verification of quantum computation.
Phys. Rev. Lett., 120:0 040501, Jan 2018.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.120.040501>.

[Gluch and Urbanke, 2020]
G. Gluch and R. L. Urbanke.
Query complexity of adversarial attacks.
CoRR, abs/2010.01039, 2020.
URL <https://arxiv.org/abs/2010.01039>.

[Gluch and Urbanke, 2021]
G. Gluch and R. L. Urbanke.
Adversarial robustness: What fools you makes you stronger.
CoRR, abs/2102.05475, 2021.
URL <https://arxiv.org/abs/2102.05475>.

[Goldwasser et al., 2020]
S. Goldwasser, A. T. Kalai, Y. Kalai, and O. Montasser.
Beyond perturbations: Learning guarantees with arbitrary adversarial
test examples.
In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin,
editors, Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, 2020.

[Ji et al., 2020]
Z. Ji, A. Natarajan, T. Vidick, J. Wright, and H. Yuen.
CoRR, abs/2001.04383, 2020.
URL <https://arxiv.org/abs/2001.04383>.

[Kitaev et al., 2002]
A. Y. Kitaev, A. H. Shen, and M. N. Vyalyi.
Classical and Quantum Computation.
American Mathematical Society, USA, 2002.
ISBN 0821832298.

[Macchiavello et al., 1997]
C. E. Macchiavello, B. Y. R. Cleve, A. Ekert, and C. Macchiavello.
Quantum algorithms revisited.
In Proceedings of the Royal Society of London A, pages
339–354, 1997.

[Mahadev, 2018]
U. Mahadev.
Classical verification of quantum computations.
2018 IEEE 59th Annual Symposium on Foundations of Computer
Science (FOCS), pages 259–267, 2018.

[Montasser et al., 2021]
O. Montasser, S. Hanneke, and N. Srebro.
Adversarially robust learning with unknown perturbation sets.
In M. Belkin and S. Kpotufe, editors, Proceedings of Thirty
Fourth Conference on Learning Theory, volume 134 of Proceedings of
Machine Learning Research, pages 3452–3482. PMLR, 15–19 Aug 2021.
URL <https://proceedings.mlr.press/v134/montasser21a.html>.

[Nguyen et al., 2015]
A. M. Nguyen, J. Yosinski, and J. Clune.
Deep neural networks are easily fooled: High confidence predictions
for unrecognizable images.
In CVPR, pages 427–436. IEEE Computer Society, 2015.
ISBN 978-1-4673-6964-0.

[Nielsen and Chuang, 2011]
M. A. Nielsen and I. L. Chuang.
Quantum Computation and Quantum Information: 10th Anniversary
Cambridge University Press, USA, 10th edition, 2011.
ISBN 1107002176.

[Quionero-Candela et al., 2009]
J. Quionero-Candela, M. Sugiyama, A. Schwaighofer, and N. D. Lawrence.
Dataset Shift in Machine Learning.
The MIT Press, 2009.
ISBN 0262170051.

[Raghunathan et al., 2018]
A. Raghunathan, J. Steinhardt, and P. Liang.
Certified defenses against adversarial examples.
In 6th International Conference on Learning Representations,
ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track
Proceedings, 2018.
URL <https://openreview.net/forum?id=Bys4ob-Rb>.

[Recht et al., 2019]
B. Recht, R. Roelofs, L. Schmidt, and V. Shankar.
Do ImageNet classifiers generalize to ImageNet?
In K. Chaudhuri and R. Salakhutdinov, editors, Proceedings of
the 36th International Conference on Machine Learning, volume 97 of
Proceedings of Machine Learning Research, pages 5389–5400. PMLR,
09–15 Jun 2019.
URL <https://proceedings.mlr.press/v97/recht19a.html>.

[Reichardt et al., 2013]
B. W. Reichardt, F. Unger, and U. Vazirani.
A classical leash for a quantum system: Command of quantum systems
via rigidity of chsh games.
In Proceedings of the 4th Conference on Innovations in
Theoretical Computer Science, ITCS '13, page 321–322, New York, NY, USA,
2013. Association for Computing Machinery.
ISBN 9781450318594.
URL <https://doi.org/10.1145/2422436.2422473>.

[Szegedy et al., 2014]
C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, and
R. Fergus.
Intriguing properties of neural networks.
In International Conference on Learning Representations, 2014.
URL <http://arxiv.org/abs/1312.6199>.

[Taori et al., 2020]
R. Taori, A. Dave, V. Shankar, N. Carlini, B. Recht, and L. Schmidt.
Measuring robustness to natural distribution shifts in image
In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin,
editors, Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, 2020.

[Tomamichel, 2015]
M. Tomamichel.
Quantum Information Processing with Finite Resources:
Mathematical Foundations.
Springer Publishing Company, Incorporated, 1st edition, 2015.
ISBN 3319218905.

[Tramer and Boneh, 2019]
F. Tramer and D. Boneh.
Adversarial training and robustness for multiple perturbations.
In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural
Information Processing Systems, volume 32, pages 5866–5876. Curran
Associates, Inc., 2019.

[Tsipras et al., 2019]
D. Tsipras, S. Santurkar, L. Engstrom, A. Turner, and A. Madry.
Robustness may be at odds with accuracy.
In International Conference on Learning Representations, 2019.
URL <https://openreview.net/forum?id=SyxAb30cY7>.

[Unruh, 2016]
D. Unruh.
Computationally binding quantum commitments.
In M. Fischlin and J. Coron, editors, Advances in Cryptology -
EUROCRYPT 2016 - 35th Annual International Conference on the Theory and
Applications of Cryptographic Techniques, Vienna, Austria, May 8-12, 2016,
Proceedings, Part II, volume 9666 of Lecture Notes in Computer
Science, pages 497–527. Springer, 2016.
URL <https://doi.org/10.1007/978-3-662-49896-5_18>.

[Vidick, 2021]
T. Vidick.
Interactions with quantum devices.
<http://users.cms.caltech.edu/ vidick/teaching/fsmp/fsmp.pdf>,
Accessed: 2021-08-30.

[Wong and Kolter, 2018]
E. Wong and J. Z. Kolter.
Provable defenses against adversarial examples via the convex outer
adversarial polytope.
In Proceedings of the 35th International Conference on Machine
Learning, ICML 2018, Stockholmsmässan, Stockholm, Sweden, July 10-15,
2018, pages 5283–5292, 2018.
URL <http://proceedings.mlr.press/v80/wong18a.html>.

Verifier    Prover

p := 0, S := ∅

i ∈{1,…N}  where  N = O (1/η^2 log(1/δ))

|ψ_^A⟩  = (ℂ^2)^⊗n

= 0:

Set  b_1, b_2, …, b_n  to be the measurement of

|ψ_^A⟩ in the $Z$ basis

S := S ∪{(b_j)_j∈{1 ,…, n}}

= 1:

Set  p_i  to be the measurement of the out register of

|0⟩_|ψ_^A⟩_|0^⊗n⟩_ in the Z basis

p := p + p_i/N

p < 1 - 2η^2

Set  x  to be an element of  S  chosen uniformly at random

The interactive protocol for the model where the verifier has access to a quantum computer and the prover doesn't need to act in an i.i.d. fashion.