Hadrien Pinault

Elodie Arlaud

Etienne Balmes
email: balmes@sdtools.com

G Degrande

A general superelement generation strategy for piecewise periodic media

Keywords:

Structures composed of repetitions of multiple identical cells are common and have been the object of a large body of literature on waveguides, periodic media, and cyclic symmetry.

Starting from a cell model, possibly with a large number of Degrees of Freedom both inside the cell and on its edges, the objective of this paper is to propose a Ritz-Galerkin reduction procedure retaining the standard second order model form and periodicity properties of the original model, while controlling accuracy in terms of model bandwidth and reproduction of the forced response to applied loads. The procedure is classically decomposed in two phases: subspace learning and basis generation. For the learning phase, Wave Finite Element (WFE) and periodic computations are presented as alternatives. The latter are then preferred for their easier control on the reduced model accuracy using classical modal synthesis and the simple choice of few target wavelengths. For the basis generation phase, constraints needed to generate a periodic superelement are defined and numerical procedures to generate a basis verifying those constraints are proposed. The validity of the reduction is demonstrated for the case of a cell with random elastic properties presenting bandgaps, local modes and wavemode crossing. Accurate predictions of modes and damped forced response are given for both the infinite and finite cases, using frequency and time simulations. The proposed analysis illustrates tracking of waveshapes, evaluation of significant waves by computation of the forced response in the frequency/wavenumber domain and interpretation of the relation between the infinite and finite forced responses. The case of a railway track with edges and transitions between multiple periodic zones is finally used to illustrate scalability issues.

Introduction

Periodic structures, composed by the repetition of a reference cell, occur in many engineering applications and thus are represented by numerical models. There are two main types of periodic structures: periodicity around an axis (bladed disks, tires …), or along a line (beams, plates, meta-materials, train tracks …). To compute the response of such a structure to an excitation, full 3D-Finite Element (FE) models typically lead to an unacceptably high numerical cost, if the cell model needs to be refined and the number of repetitions is significant. Thus it is often useful to consider the structure as a waveguide [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF] and analyse its response in terms of waves propagating in the repetition direction.

A first group of methods considers the limit case of a structure invariant along the repetition direction. Such problems can be solved with semi-analytical finite element (SAFE) methods [START_REF] Marzani | A semi-analytical finite element formulation for modeling stress wave propagation in axisymmetric damped waveguides[END_REF][START_REF] Loveday | Semi-analytical finite element analysis of elastic waveguides subjected to axial loads[END_REF] or 2.5D approaches [START_REF] Yang | A 2.5d finite/infinite element approach for modelling visco-elastic bodies subjected to moving loads[END_REF][START_REF] Franois | A 2.5d coupled FE-BE methodology for the dynamic interaction between longitudinally invariant structures and a layered halfspace[END_REF][START_REF] Costa | Trackground vibrations induced by railway traffic: in-situ measurements and validation of a 2.5d FEM-BEM model[END_REF]. These approaches share the idea of considering the mechanical equations in the frequency/wave-number domain. For the trivial example of a beam of modulus E = 68.8 MPa, density 𝜌 = 1700 kg/m 3 and section A, excited in traction/compression, the space/time equations 𝜌A 𝜕 2 u 𝜕t 2 (x, t) -EA 𝜕 2 u 𝜕x 2 (x, t) = f (x, t)

(1) become in the frequency/wavenumber domain a rational transfer from force to displacement U(𝜅, 𝜔) F(𝜅, 𝜔)

= 1 -𝜌A𝜔 2 + EA𝜅 2 , (2)
where, at a given frequency 𝜔, free waves are solutions in the absence of excitation and thus have wavenumbers 𝜅 that verify -𝜌A𝜔 2 + EA𝜅 2 = 0, the so-called dispersion equation. Similarly, at a given wavenumber 𝜅, periodic modes are solutions in the absence of excitation and have poles that verify the same equation.

For structures that have an internal structure (blades, tracks, meta-materials …), or simply if a cell length is considered to bypass the need to develop specific elements in the SAFE methodology, space can no longer be considered on continuous line.

Responses are thus sampled in space and the Bloch or Floquet transforms (i.e. space DFT, for discrete Fourier transforms) are thus considered to analyse the dependence on wavenumbers, which are now defined on a real interval for infinite structures or on a discrete set for finite structures. Due to sampling, the spatial transform is now subject to aliasing: similarly to the timealiasing that occurs at

f s 2
when sampling a signal at a sampling frequency f s , a space-aliasing (or Bragg aliasing [START_REF] Boukadia | A Bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF]) occurs at

x 0 2 in a periodic medium of cell size x 0 . For a spatial step of x 0 = 0.4 m, Fig. 1 illustrates the transfer (2) in the frequency/wavenumber domain, which will be called 2D-DFT later on.

The main proposition of the paper is that periodic solutions can be used to build a cell superelement based on a Ritz-Galerkin reduction procedure with well-understood levers to control accuracy in terms of bandwidth and applied loads. Ritz-Galerkin model reduction methods can generally be decomposed in a subspace learning phase (static and eigen-solutions for mode synthesis [START_REF] Craig | A review of time-domain and frequency domain component mode synthesis methods[END_REF], snapshots for non-linear model reduction [START_REF] Kerschen | The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview[END_REF] …), followed by a basis generation algorithm where the independence of vectors is ensured and possibly other constraints are taken into account.

The assumption that the structure, or parts of it, is composed of a repetition of identical cells has led in the literature to two kinds of developments to efficiently evaluate the 2D-DFT.

A first category, called direct approaches in Ref. [START_REF] Boukadia | A Bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF], chooses to work at a given frequency (vertical line in Fig. 1). This has for example led to the Wave Finite Element (WFE) method (see e.g. Refs. [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF][START_REF] Mencik | A model reduction strategy for computing the forced response of elastic waveguides using the wave finite element method[END_REF][START_REF] Mencik | A wave-based model reduction technique for the description of the dynamic behavior of periodic structures involving arbitrary-shaped substructures and large-sized finite element models[END_REF][START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF][START_REF] Serra | Wave properties in poroelastic media using a wave finite element method[END_REF]). Transition matrices for displacements and loads defined on the left/right edges of a cell are defined. A numerical gain is achieved by diagonalizing the transition matrix to obtain wave definitions and simplify the expression of the successive powers of the transition matrix to describe propagation over an arbitrarily large number of cells. When considering problems with large interfaces, numerical cost can become significant. Consequently, model reduction has been considered in Refs. [START_REF] Mencik | A model reduction strategy for computing the forced response of elastic waveguides using the wave finite element method[END_REF][START_REF] Mencik | A wave-based model reduction technique for the description of the dynamic behavior of periodic structures involving arbitrary-shaped substructures and large-sized finite element models[END_REF][START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF] and strategies to preserve sparsity have been introduced in Ref. [START_REF] Collet | Floquet-Bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems[END_REF].

A second category, called inverse approaches in Ref. [START_REF] Boukadia | A Bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF], chooses to work at a given wavenumber (horizontal line in Fig. 1). This is the classical choice for cyclic symmetry computations [START_REF] Thomas | Standing waves in rotationally periodic structures[END_REF][START_REF] Sternchuss | Reduction of multi-stage disk models: application to an industrial roto[END_REF], where there is a finite number of spatial frequencies often called nodal diameters. For infinite waveguides, sample studies are [START_REF] Chebli | Dynamic response of high-speed ballasted railway tracks: 3D periodic model and in situ measurements[END_REF][START_REF] Arlaud | Receptance of railway tracks at low frequency: numerical and experimental approaches[END_REF]. The numerical gain comes from the separation of problems associated with each spatial harmonic, so that a full order periodic problem can be solved by considering a series of smaller problems, each related to an harmonic, and then performing the inverse spatial transform to rebuild the full solution. The forced response to a periodic load of a chosen wavenumber can be efficiently computed as shown in section 2.3. An advantage of this procedure is that, at a given wavelength, model reduction is a trivial application of modal synthesis where model accuracy is simply controlled by keeping modes within the band of interest and using a static correction for applied loads. This approach is addressed in section 2.

The work presented in this paper is a direct continuation of earlier cyclic [START_REF] Sternchuss | Reduction of multi-stage disk models: application to an industrial roto[END_REF] and periodic [START_REF] Arlaud | Receptance of railway tracks at low frequency: numerical and experimental approaches[END_REF] applications, which addressed parametric studies of cell properties (mistuning) and use of reduced models in non-linear transient simulations (train/track interaction). Boukadia et al. [START_REF] Boukadia | A Bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF] have later published a similar methodology, which essentially differs by the learning phase of the reduction and with illustrations focusing on the classical frequency domain responses. The lower cost of the strategy proposed here will thus be detailed in section 3 and innovation will be sought in the way to track waves in a dispersion diagram, to relate dispersion diagram and forced response and finally to relate forced response of an infinite repetition and forced response of a finite structure.

The proposed superelement construction procedure, detailed and illustrated in section 3 uses the fact that multi-model reduction [START_REF] Balmes | Parametric families of reduced finite element models. Theory and applications[END_REF][START_REF] Hammami | Numerical design and test on an assembled structure of a bolted joint with viscoelastic damping[END_REF] provides an efficient fixed basis reduction for the full spatial spectrum using periodic computations to build a modal synthesis basis at just a few learning wavenumbers. To prove that reduction is possible, subsection 3.2 illustrates that the subspace spanned by the periodic solutions does not grow significantly with the number of wavenumbers considered in a multi-wavelength resolution. Then, it is detailed how the desired end result is a periodic superelement, i.e. a model with standard second order equations considering mass, stiffness, damping matrices and Degrees Of Freedom (DOF) split as interior and matching left/right sets. From any starting learning subspace, a construction procedure is proposed to obtain such a superelement.

The example of a randomly-fluctuating, periodic beam is used to illustrate the validity of the proposed procedure. This example exhibits local modes (whose influence is limited to a small part of space), which interact with global modes (pressure and bending). Thus a significant part of wave energy is absorbed by the local modes: this phenomenon is called localization [START_REF] Castanier | Modeling and analysis of mistuned bladed disk vibration: status and emerging directions[END_REF]. Another property of this random medium is the apparition of frequency bandgaps, where energy cannot propagate, which can be used to design absorbing meta-materials [START_REF] Hussein | Metadamping: an emergent phenomenon in dissipative metamaterials[END_REF]. All those features are shown to be well reproduced with major numerical gains for an implementation using standard superelements in the Structural Dynamics Toolbox (SDT) [START_REF] Balmes | Structural Dynamics Toolbox (For Use with MATLAB)[END_REF].

The superelements resulting from the procedure of section 3 retain the standard form of a differential equation in time and can thus be used in a very general setting. This however does not address the problem of edges of finite structures or connections between finite periodic parts. Section 4 proposes to keep the reduction of a infinite medium for most of the structure and use interface models to represent the non-periodic areas close to the edges. The methodology is illustrated for the case of a finite heterogeneous beam, and of a railway track with a transition from a ballasted to a slab track thus illustrating the scalability and efficiency of the proposed methodology.

Forced response in the frequency/wavenumber domain

The purpose of this section is to highlight how to compute a forced response in the frequency domain at a given wavenumber (using periodic computations), contrary to the WFE method (in the wavenumber domain at a given frequency).

Direct and inverse spatial transforms

First the conventions for discrete spatial transforms are clarified. For any point 𝜉 in a reference cell Ω 0 , its geometric translations/rotations spaced by x 0 are defined by the series

( 𝜉 + nx 0 )
n∈E . Any repetition of Ω 0 is then denoted Ω n . An example of such a structure is presented in Fig. 2. Cases of interest are E = ℤ for infinite structures or E = [0 N] , N ∈ ℕ * for cyclic symmetry. There are several expressions of the Fourier, Bloch or Floquet transforms of field f depending on normalization. In the following, the DFT is used and defined by

F(𝜅, 𝜉) = ∑ n∈E f (𝜉 + nx 0 )e -2i𝜋𝜅n , ∀𝜅 ∈ [0 1] , ∀𝜉 ∈ Ω 0 , (3)
with 𝜅 the dimensionless wavenumber (or spatial frequency) expressed in cell -1 . The value of f at any point x = 𝜉 + nx 0 of space can be computed back with the inverse transform (IDFT), given by

f (𝜉 + nx 0 ) = ∫ 1 0 F(𝜅, 𝜉)e 2i𝜋𝜅n d𝜅. (4)
An important property of the DFT is the so-called (𝜅, x 0 )-Bloch periodicity relation, valid for any n ∈ E, and any periodic field of wavenumber 𝜅 0

F ( 𝜅 0 , 𝜉 + x 0 ) = e 2i𝜋𝜅 0 F ( 𝜅 0 , 𝜉 ) .
(

) 5
For more details on this definition and its relation with the DFT, refer to Ref. [START_REF] Allaire | The Bloch transform and applications[END_REF].

Equation of motion of a cell model

For a given cell Ω n , one assumes linear equations of motion derived from any FE software of the form

( -𝜔 2 M + K ) q n (𝜔) = b f n (𝜔), (6)
with q n the DOF, 𝜔 the angular frequency, b the spatial position of the loads and f n their evolution in frequency. Assuming the structure to be periodic implies that any cell Ω n is simply obtained by translation (one then talks of periodic media), or rotation (for cyclic periodicity) of the reference cell Ω 0 . The matrices are thus independent of n, at least until piecewise periodicity is introduced in section 3. Without damping, the mass M and damping K matrices are real-valued. The development also holds for cases with viscous damping (i.e. i𝜔C matrix), hysteretic (i.e. K complex-valued) or viscoelastic (i.e. K(𝜔, T) complex-valued, depending on frequency and temperature [START_REF] Nashif | Vibration Damping, john wiley and sons Edition[END_REF]). To have more concise notations, the paper will use the dynamic stiffness matrix Z defined by

Z(𝜔) = -𝜔 2 M + K. (7)
DOF will be labelled I if they are inside their cell Ω n , R if they are on the right side, and L on the left side (considering a fixed axis direction). Interfaces DOF are common to adjacent cells and thus play a key role: right DOF of any cell Ω n are identical to the left DOF of the following cell Ω n+1 (see Fig. 2.b). This leads to the following continuity equation

q R n = q L n+1 , ∀n ∈ E. (8)
Eq. ( 6) in Ω n can then be rewritten as

⎡ ⎢ ⎢ ⎢ ⎣ Z LL Z LI Z LR Z IL Z II Z IR Z RL Z RI Z RR ⎤ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ q L n q I n q R n ⎫ ⎪ ⎬ ⎪ ⎭ = b ⎧ ⎪ ⎨ ⎪ ⎩ f L n f I n f R n ⎫ ⎪ ⎬ ⎪ ⎭ . ( 9
)
Cases that this paper seeks to address are those where the cell model is large. In the litterature, two main methods are used to reduce this problem. The first is known as dynamic condensation. If there are many internal DOF and no external forces on internal nodes (F I n = 0), the interior DOF can be eliminated [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF][START_REF] Duhamel | Finite Element Analysis of the Vibrations of Waveguides and Periodic Structures[END_REF] with the relation

q I = -Z -1 II ( Z IL q L + Z IR q R ) . ( 10
)
Eq. ( 9) can then be reduced to

[ ZLL ZLR ZRL ZRR ] { q L n q R n } = b { f L n f R n } , (11)
where

[ ZLL ZLR ZRL ZRR ] = [ Z LL -Z LI Z -1 II Z IL Z LR -Z LI Z -1 II Z IR Z RL -Z RI Z -1 II Z IL Z RR -Z RI Z -1 II Z IR ] . ( 12
)
If the cost of factoring Z II at multiple frequencies is large, a Craig-Bampton procedure [START_REF] Craig | A review of time-domain and frequency domain component mode synthesis methods[END_REF] can be considered as an alternative methodology. First, the modes 𝜑 I of the internal DOF are computed with fixed boundaries as the non-trivial solutions of

Z II 𝝋 I = 0. ( 13
)
Then the DOF are redefined with hybrid DOF containing: physical DOF on the interfaces q R , q L and a significantly smaller number of internal modal DOF p I . The initial internal DOF q I can then be recovered using the constant linear relation

q I = -K -1 II K IL q L + 𝝋 I p I -K -1 II K IR q R . ( 14
)
Nevertheless, for large interfaces (i.e. a high number of DOF in q L/R ), computations of the condensed dynamic stiffness Z will remain costly and the methodology proposed in this paper gives a solution.

Forced response in the frequency domain at a given wavenumber

To propose an alternative method to compute the 2D-DFT, one now starts doing the spatial DFT and thus chooses a wavenumber 𝜅 0 . Solutions that verify the periodicity Eq. ( 5) are then searched. This is a classical approach for a structure composed of cyclically repeated cells with constant properties. It is a well-known result (see Ref. [START_REF] Sternchuss | Reduction of multi-stage disk models: application to an industrial roto[END_REF] for example) that the transformation of a full 3D model to harmonic coordinates leads to a series of fully decoupled problems, associated with the generalized DOF related with each harmonic. This result, being derived from Fourier properties, can also be extended to an infinite periodic structure. Periodic solutions to be found verify the Eqs. of motion ( 6) on every cell and the continuity condition ( 8). This continuity is expressed in the frequency domain with Eq. ( 5), leading to [START_REF] Collet | Floquet-Bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems[END_REF] with c R (resp. c L ) an observation matrix which extracts the DOF on the right side (resp. left side), and Q(k 0 , 𝜔) being a Fourier transform in both time and space with 4 components at every DOF (combinations of real/imaginary parts in time and space).

c L Q(𝜅 0 , 𝜔) = e 2i𝜋𝜅 0 c R Q(𝜅 0 , 𝜔),
Starting with the time transform leads to the complex q(x, 𝜔). Doing the spatial transform leads to spatially real and imaginary parts of Q(k 0 , 𝜔). Thus Eq. ( 15) is rewritten in a matrix form

c(𝜅 0 ) • { Re ( Q(𝜅 0 , 𝜔) ) Im ( Q(𝜅 0 , 𝜔) ) } = 0, (16)
where the continuity matrix c, independent of frequency 𝜔, is

c(𝜅 0 ) = [ c L -cos(𝜅 0 )c R -sin(𝜅 0 )c R sin(𝜅 0 )c R c L -cos(𝜅 0 )c R ] . ( 17
)
In Eq. ( 16), Re and Im refer only to the spatial transform: thus Re(Q) and Im(Q) are complex-valued, containing real and imaginary part of the time transform.

A problem with a linear constraint can be solved using Lagrange multipliers, penalization or elimination. In the latter case, chosen in this work, a basis T c of the kernel of c is computed. More precisely the columns of matrix T c form a basis which spans the kernel and the abuse of language will be used later on as usually in the model reduction literature. The constrained problem for a given wavelength 𝜅 0 is finally written as

T T c [ Z 0 0 Z ] T c { Re(Q c ) Im(Q c ) } = T T c [ b 0 0 b ] { Re (F) Im (F) } , ( 18
)
where T c for every angular frequency 𝜔 can be very costly, but this is the classic problem of modal synthesis. It is well-known that a modal approximation with static correction (i.e. MacNeal reduction [START_REF] Craig | A review of time-domain and frequency domain component mode synthesis methods[END_REF]) is perfectly suited. That approximation has the form

Q = T c Q c and
{ Re ( Q c ) Im ( Q c ) } ≈ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ NM ∑ j=1 𝝋 j 𝝋 T j -𝜔 2 + 𝜔 2 j ⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟ modal approximation + K -1 c - NM ∑ j=1 𝝋 j 𝝋 T j 𝜔 2 j ⏟⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏟ static correction ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ T T c [ b 0 0 b ] { Re (F) Im (F) } , ( 19
)
where harmonic modes, with eigenvectors

( 𝝋 j (𝜅 0 )
) j and eigenfrequencies

( 𝜔 j (𝜅 0 ) ) j , are solutions of a standard eigenvalue problem Z c (𝜔 j ) { Re ( 𝝋 j ) Im ( 𝝋 j ) } = 0. ( 20
)
In Eq. ( 19), K c is defined similarly to Z c . In the following, for sake of readability, Z c , Q c and similar fields will be assimilated to Z, Q, etc. Damped problems can also be treated using appropriate reduction that will be detailed in section 3.4.

Though fully equivalent to WFE [START_REF] Mencik | A model reduction strategy for computing the forced response of elastic waveguides using the wave finite element method[END_REF][START_REF] Mencik | A wave-based model reduction technique for the description of the dynamic behavior of periodic structures involving arbitrary-shaped substructures and large-sized finite element models[END_REF][START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF][START_REF] Serra | Wave properties in poroelastic media using a wave finite element method[END_REF], this method however provides a clearer criterion to achieve model reduction. Since classical modal synthesis can be used, one simply needs to specify a bandwidth leading to a selection of NM modes and a static correction corresponding to the applied loads b. A multi-wavelength reduction strategy will be introduced in the next section to alleviate the cost of computing the response for a large number of wavelengths.

Superelement building using a multi-wavelength learning phase

The objective of this section is to propose a superelement generation procedure that will give a model with mass, damping, and stiffness properties as well as left, right and interior DOF, that is able to reproduce the 2D-DFT of an infinite periodic medium with a controlled accuracy. The extension to finite media will be addressed in section 4. The properties of the non-trivial heterogeneous medium chosen as example are detailed in section 3.1.

Section 3.2 first estimates a priori the reduction relevance by showing that the subspace spanned by the forced response is of low rank. The second step is to build a standard cell model with left, right and interior DOF and section 3.3 proposes an associated strategy starting from any subspace spanning the forced response. As a subspace generation strategy the multi-model reduction [START_REF] Balmes | Parametric families of reduced finite element models. Theory and applications[END_REF][START_REF] Hammami | Numerical design and test on an assembled structure of a bolted joint with viscoelastic damping[END_REF] is applied, using learning at a few wavelengths, and shown to be efficient.

Test case: 2D-DFT computation of an infinite periodic heterogeneous beam

To show the usefulness of the model reduction, a non-trivial case is chosen: an infinite periodic heterogeneous beam with a large cell model. The map of properties is taken from 3D, periodical, random fields, following a log-normal distribution. They are efficiently generated with a spectral method proposed by Shinozuka and Deodatis [START_REF] Shinozuka | Simulation of multi-dimensional Gaussian stochastic fields by spectral representation[END_REF]. For example, the density random field

𝜌 [x] is calculated from 𝜌 [x] = 𝔼 𝜌 (x) + 2 N rand ∑ n=1 𝜁 n √ 𝕊 𝜌 (nΔk)Δk cos ( nΔkx + 𝜙 n ) , (21)
where 𝔼 𝜌 and 𝕊 𝜌 are the expectation and autocorrelation of 𝜌, 𝜁 n standard normal fields, 𝜙 n uniform fields on [-𝜋 𝜋] and Δk the step in wavenumber. Bulk modulus

K = E 3(1-2𝜈)
and shear modulus G = E 2(1+𝜈) are simulated similarly. The correlation model is chosen gaussian, only determined by a correlation length 𝜆 corr identical in the three directions of space. This kind of structure can be used to model heterogeneous media, such as soil [START_REF] Gatti | Investigation of the earthquake ground motion coherence in heterogeneous non-linear soil deposits[END_REF] or railway ballast [START_REF] De Abreu Corra | Randomly-fluctuating Heterogeneous Continuum Model of a Granular Medium[END_REF]. The selected parameters are summarized in Table 1.

A 3D view of the isovalues of the density field is presented in Fig. 2. This representation highlights harder points (in red) included in a softer medium (transparent). Even with the heterogeneity, the periodicity is enforced in the reference cell to avoid interface discontinuities. However, in such a beam, many waves can propagate: 2 bending waves, 2 torsion waves, 1 pressure wave. To simplify the analysis, the beam is also supposed periodic in the transverse direction and lies on an elastic foundation with a constant stiffness (see Fig. 6.b). Then, only 1 pressure wave and 1 bending wave can propagate, and the problem of rigid modes is avoided.

Before analysing the forced response of a heterogeneous beam, it is interesting to consider the dispersion diagram, displaying the evolution of the modes 𝜔 k (𝜅 j ) for all k and j. Following the previously described method, a set of wavenumbers ( 𝜅 j ) j is fixed, and for each 𝜅 j , the associated modes 𝜔 k (𝜅 j ) are computed.

For a fully homogeneous beam (with the mean values of Table 1), the dispersion diagram is presented in Fig. 3.a. Two classes of modes are visible: bending modes (blue crosses) and pressure modes (red circles). To automatically identify mode classes, modal coordinates [START_REF] Bianchi | Using modal damping for full model transient analysis[END_REF] are used. Taking a reference wavelength 𝜅 0 (typically at 𝜅 ≈ 0), the periodic modes (20) are mass orthogonal. One can thus estimate modal amplitudes 𝛼 j of mode 𝜑 j using and use as marker that of the highest contributing reference mode.

𝛼 j = { Re ( 𝝋 j (𝜅 0 ) ) Im ( 𝝋 j (𝜅 0 ) ) } T [ M 0 0 M ] { Re ( Q(𝜅 0 , 𝜔) ) Im ( Q(𝜅 0 , 𝜔) ) } (22)
As expected, the evolution of the pressure mode frequency is quite linear, i.e. not dispersive, contrary to the bending mode, whose frequency has a quadratic evolution and is thus dispersive.

With the introduction of the heterogeneity, the dispersion diagram is modified, as shown in Fig. 3.b. When a periodic frequency reaches a wavenumber of 0.5 or 0, aliasing still occurs but a bandgap (range of frequencies where no waves propagate) is now present because the structure actually has a corresponding geometric wavelength. This behaviour is considered in attenuation studies using meta-materials [START_REF] Hussein | Metadamping: an emergent phenomenon in dissipative metamaterials[END_REF]. Also, the slopes of bending and pressure curves are also a bit steeper, which means that the associated wave speed is slower than in the homogenous case. This is a classical result of the homogeneization of random media [START_REF] Capdeville | 1-D non-periodic homogenization for the seismic wave equation[END_REF].

In the heterogeneous medium, the hypothesis of decoupled modes is also no longer verified: a gap is visible at the crossing bending and pressure waves. At frequencies within the gap, only evanescent waves exist.

Finally, an additional mode is present around 200 Hz (pink diamonds). The slope of the corresponding curve is very high, which means that the associated wave speed is almost zero. This particular mode will be named local, as the corresponding wave energy has a very low speed and consequently stays confined in a small portion of space. To clarify the distinction between these local and global modes, the shapes of few modes are now presented. On Fig. 4, bending and pressure modes are plotted on 5 consecutive cells (the selected modes are labelled in Fig. 3.b). They can be named global because their mode shape is not restricted to a small zone of the cell, but spread over the whole structure. Fig. 5 then illustrates three points around the local mode interaction. Far from the interaction, the global mode (Fig. 5.a) and the local mode (Fig. 5.c) have a distinct nature. The global mode is mainly pressure and the mode shape is spread over all cells. On the contrary, the local mode has a shape mainly limited to one portion of the cell, and the repetition is only due to the periodicity of the structure. Closer to the interaction (Fig. 5.b), the shape is mixed and difficult to interpret.

These remarks still apply to any sampling of the heterogeneous medium. However, the slopes of the curves, the width of the bandgaps, and the frequencies of the local modes will vary.

Finally, the forced response of the heterogeneous beam is computed, with the method developed in section 2. The beam (periodic along the x-axis) is excited by a uniform shear force at x = 0 in the vertical direction, as presented in Fig. 6.b. The mean vertical acceleration of the section x = 0 is computed and the 2D-DFT of the acceleration is shown in Fig. 6c, with colours proportional to amplitude that is visible in the 3D view (Fig. 6d). The superimposed black curves are taken from the dispersion diagram in Fig. 3.b. Fig. 6.a presents the time-FFT (Fast Fourier Transform) of the acceleration. Placed above the 2D-DFT it eases the understanding of the correspondence between the two plots: indeed, the FFT can be seen as the integral of the 2D-DFT along the wavenumber axis. Thus, bandgaps at aliasing points or mode crossings are visible as anti-resonances on the FFT (Fig. 6.a) or valleys on the 2D-DFT (Fig. 6c andd).

As expected, the bending modes are particularly excited, i.e. the amplitude of the 2D-DFT is maximum on the bending wave curve. The presence of bandgaps is clearly visible as the existence of valleys in the response. At the crossing with the pressure mode at 82 Hz, a part of the energy is shared between pressure and bending. Also, the local mode is excited with a similar amplitude and thus cannot be neglected.

In the 3D view of the 2D-DFT in Fig. 6d, irregular lines are visible at low frequencies and wavenumbers. These are a wellknown numerical artefact, due to the discrete nature of computation which typically misses the maximum of the response for low damping values and insufficient frequency spacing.

The computation of this forced response may require huge computing resources if the cell model is large and furthermore if it has many interface DOF. The following subsections will describe the method proposed to reduce this cost.

A priori evaluation of reduction relevance

The main purpose of the proposed reduction method is to build a superelement of slices based on a second level of Ritz Galerkin method (the FE method already being such an approximation). The goal is then to find a constant matrix T (i.e. inde-pendent of 𝜅 and 𝜔) whose columns span the Ritz subspace and reduced (or generalized) coordinates Q r (𝜅, 𝜔) such that the Ritz approximation of the response is of the form Q(𝜅, 𝜔) = TQ r (𝜅, 𝜔) . [START_REF] Hussein | Metadamping: an emergent phenomenon in dissipative metamaterials[END_REF] Note that T is of size N DOF × N r , where N DOF is the number of DOF in Ω 0 and N r the number of DOF in the reduced model. Model reduction is useful only if N r << N DOF is achievable.

To estimate an accurate N r , a Singular Value Decomposition (SVD) can be used with, as usual in mechanics [START_REF] Balmes | Optimal Ritz vectors for component mode synthesis using the singular value decomposition[END_REF], the strain energy chosen as a norm. A first set of N vectors containing the eigenvectors

( 𝝋 k ( 𝜅 j ))
1≤k≤N , solutions of Eq. ( 20) obtained at N 𝜅 given wavenumbers

( 𝜅 j ) 1≤j≤N 𝜅
, is stored as columns of matrix T per ,

T per = [ 𝝋 1 ( 𝜅 1 ) • • • 𝝋 N ( 𝜅 1 ) • • • 𝝋 1 ( 𝜅 N 𝜅 ) • • • 𝝋 N ( 𝜅 N 𝜅 )] . ( 24
)
Then, using the stiffness matrix K, the strain energy in the modal basis  is computed by

 = T T per KT per . ( 25
)
Finally the classical SVD algorithm is applied to  leading to an estimate of N r as the number of significant singular values. The associated singular vectors can then be used as columns of T.

To have an overview of the global behaviour of the track, 20 wavenumbers are chosen linearly between 2 and 100 cells and the associated periodic modes are kept. In the frequency band of interest, Fig. 3.b shows 3 modes (pressure, bending and local). Thus, only the first 4 modes are kept for each wavelength. The decrease of strain energy rate linked to the subspace of the heterogeneous beam is displayed in Fig. 7 normalized by the predominant shape. It appears that the decrease rate is very fast, because only N r = 6 shapes are necessary to represent more than 99% of the energy, compared to the 20 × 4 = 80 initial vectors. Moreover, a clear drop appears after the 6th shape. 20 samples from the same random properties are presented in grey, together with the mean value (in red). This proves that reduction is relevant in this case.

The case of high wavenumbers (or short wavelengths close to the half spectrum minimum of two cells) is a limit. The shape within the cell then becomes very sensitive to wavenumber. This corresponds to the usual limit of FE models, where one needs at least a few elements per wavelength. It is however interesting to note that local modes within the element are still quite well reproduced as will be illustrated later.

Superelement generation process for a periodic waveguide

The next section will discuss the learning phase where simple computations are used to generate a subspace, associated with the span of a 𝚯 matrix, that is close to containing the desired responses. This section first focuses on showing that for any 𝚯, a clear procedure exists to generate a matrix T whose columns form a basis and verify some geometry constraints that make it suitable for use as a the reduction subspace of a superelement.

The matrix T is searched of the form

⎧ ⎪ ⎨ ⎪ ⎩ q L q I q R ⎫ ⎪ ⎬ ⎪ ⎭ N DOF = [ T L T I T R ] q r = ⎡ ⎢ ⎢ ⎢ ⎣ T LL 0 0 T LI T II T RI 0 0 T RR ⎤ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ q Lr q Ir q Rr ⎫ ⎪ ⎬ ⎪ ⎭ N R , ( 26
)
and verifying the following properties.

• the continuity between two slices for reduced DOF equality q n+1 Rr = q n Lr is assured by taking T RR = T LL , • displacements at the left interface are generated only by T L , so T IL = T LR = 0, • displacements at the right interface are generated only by T R , so T RL = T IR = 0. Starting from the matrix 𝚯, the internal DOF are separated from interfaces ones. This operation is possible because the resulting space contains the initial subspace:

span(𝚯) ⊂ span ⎡ ⎢ ⎢ ⎢ ⎣ 𝚯 L 0 0 0 0 𝚯 I 0 𝚯 R 0 ⎤ ⎥ ⎥ ⎥ ⎦ . ( 27
)
As before, the notation span(𝚯) means the linear span of the set of vectors associated with the columns of the matrix 𝚯.

First, the matrix blocks T RR = T LL are constructed. A unique subspace is considered for left and right interfaces span

[ 𝚯 L 𝚯 R ]
. Then T RR is obtained by orthogonalization of the matrix

[ 𝚯 L 𝚯 R ]
. If classical iterative Gram-Schmidt method is a robust solution, the generation of a complete subspace is not compulsory to get an accurate solution. A partial truncation using SVD can be used instead. Here a truncation threshold was set to 10 -4 , but this value is still an open question.

Once the blocks T LL = T RR related to the interface displacements are built, the interior blocks have to be built to obtain representative displacements inside the slice. Most of the time, a simple interpolation between the interfaces is not possible, because there are huge heterogeneities inside a slice of a waveguide. Then, a classical static expansion method (also called Schur complement or Guyan reduction [START_REF] Craig | A review of time-domain and frequency domain component mode synthesis methods[END_REF]) is used. Then T LI (and similarly T RI ) is computed by

T LI = -K -1 II K LI T LL , T RI = -K -1 II K RI T RR . ( 28
)
Finally only T II still has to be constructed. It corresponds to displacements generated by 𝚯 I where both interfaces are not excited, which were not included neither in T IL nor T IR . In applications shown here, this is achieved by the following linear combination

T II = 𝚯 I -T L T + LL 𝚯 L -T R T + RR 𝚯 R , ( 29
)
where T + LL is the pseudo-inverse of T LL . But alternate strategies to build T II should be considered. While this work is derived from Refs. [START_REF] Sternchuss | Reduction of multi-stage disk models: application to an industrial roto[END_REF][START_REF] Arlaud | Receptance of railway tracks at low frequency: numerical and experimental approaches[END_REF], a comparison to the independent [START_REF] Boukadia | A Bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF] is useful. The static expansion of edge motion ( 28) is here explicit, when it was implicit due to the use of initial Craig-Bampton reduction of the cell. For the model considered in Table 1, the Craig-Bampton basis would be 3 GB and require a costly 2883 static evaluations. Thus bypassing for this first reduction is interesting for large interface models. For interior DOF, Ref. [START_REF] Boukadia | A Bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF] chooses to keep all fixed interface modes of the Craig-Bampton reduction. The proposal of ( 29) is to use the a basis of the learned subspace, which gives an easier control of model accuracy than choosing a high cutoff at 4 × 𝜔 max .

Subspace learning phase: a multi-wavelength strategy

In this subsection, the method to generate a matrix 𝚯 whose columns span a subspace close to containing the solution is discussed. This basis is obtained by learning, that is solving some reference problems for target wavenumbers and possibly viscoelastic damping or design parameters, and assuming that the solution for a range of parameters is well represented within the generated subspace.

The chosen reference problems are the 2D-DFT for a small set of wavenumbers

( 𝜅 j ) 1≤j≤N 𝜅
and frequency band

[ 0 𝜔 max ]
. For any such wavelength 𝜅 j , the modes

( 𝝋 k (𝜅 j )
) k in the frequency band of interest (i.e. such as 𝜔 k (𝜅 j ) < 𝜔 max ) are computed with the very classical Lanczos algorithm, and stored as columns of a matrices 𝚽(𝜅 j ). The modes are then classically complemented with the static response to applied loads, obtained by solving Eq. ( 18) for 𝜔 = 0 and 𝜅 j , because it is important to allow loads localized in space. This static response will be then referred by the notation

( K -1 b ) (𝜅 j ).
If damping or any parametric model variation [START_REF] Balmes | Parametric families of reduced finite element models. Theory and applications[END_REF][START_REF] Hammami | Numerical design and test on an assembled structure of a bolted joint with viscoelastic damping[END_REF] is considered, a first order correction K -1 Im(K)𝚽 or multi-model reduction with at least two values of the viscoelastic modulus E v1 and E v2 can also be introduced. Rouleau and al. [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF] clearly demonstrated that this approach is among the best to generate reduced damped models. Then the wanted matrix 𝚯 is obtained as the concatenation of all these components

𝚯 = [ 𝚽(𝜅 j ) ( K -1 b ) (𝜅 j ) ( K -1 Im (K) 𝚽 ) (𝜅 j ) ] 1≤j≤N N 𝜅 . ( 30
)
Here the single adjustment parameter is the selection of the wavelength set as originally proposed in Ref. [START_REF] Sternchuss | Reduction of multi-stage disk models: application to an industrial roto[END_REF]. For the same objective [START_REF] Boukadia | A Bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF], starts by a costly Craig-Bampton reduction with a cutoff frequency for fixed interface modes at a high 4 × f max .

The direct solution at 𝜔 max provides information about both modes and response to loads. Then wavenumber sampling is performed as done here, but with a much larger number of wavenumbers. Indeed, a classical result of parametric model reduction [START_REF] Balmes | Parametric families of reduced finite element models. Theory and applications[END_REF] is that modal crossing is well represented with very few learning points. Indeed, in most cases, two are sufficient. In the present study, wavenumber sampling is a form of multi-model reduction and two points will indeed be shown to be sufficient. Finally, residue iterations [START_REF] Bobillot | Iterative computation of modal sensitivities[END_REF] provide a generic mechanism to extend the learning subspace by checking the validity of reduced solutions. This could easily be applied here at additional wavenumbers.

To analyse the accuracy of this reduction methodology, it is applied to the case of the heterogeneous beam described in section 3.1. The simplest possible choice it to keep a single wavenumber 𝜅. The result obtained for 𝜅 = 0.5 (corresponding to 2 cells) is presented in Fig. 8.a for the FFT and in Fig. 9.a for the 2D-DFT. Similarly to the non-reduced case, FFT and 2D-DFT are placed above each other for easier comparison. The dashed blue curve in Fig. 8.a is the full model FFT of Fig. 6c. It appears that this choice leads to a very wrong response, although the modes for this particular wavenumber are correct (in the red rectangle). This means that the corresponding subspace is too poor to represent the full behaviour of the heterogeneous beam.

With a lower wavenumber (e.g. 𝜅 = 0.2, corresponding to 5 cells), the result is closer to the reference, but the amplitudes are not satisfying. This consideration shows that high wavenumbers are the limit of the strategy, because the stress discontinuity due to the discretisation becomes too high in comparison to the typical wavelength.

To obtain a more precise response, two wavenumbers are now selected: 𝜅 ∈ {0.01, 0.2}, corresponding to wavelengths of 100 and 5 cells. Comparing to the reference in Fig. 6c, the results after reduction (Fig. 8.b and Fig. 9.b) are excellent until 180 Hz. Of course, the two selected wavelengths (in red rectangles) are still perfectly correct as expected for a Galerkin procedure, but the response is also accurate for the other wavenumbers. The choice of these two wavelengths (long and small) seems sufficient to represent the global response of this periodic structure. Since accuracy is already excellent, selecting more wave numbers does not seem necessary. Furthermore, the exact values chosen do not affect the quality of results.

To illustrate the performance of the strategy, the forced response of the random beam has been computed with a refined mesh to be sure to catch the heterogeneity properties (h = 1 cm). Consequently, the number of interface DOF is quite high, as detailed in Table 2. With this strategy, the number of DOF is reduced by 3 orders of magnitude, and so is the computation time and the memory required.

Error evaluation for the reduced model

Most studies on periodic damped media focus on the prediction of wavemodes with complex wavenumbers at a given frequency. Since our focus is on forced response for models that have a time domain formulation, computing these values require unnecessary developments and, on cases with large cell models and extended interfaces, the cost of exact computations of these wavenumbers would be large. Furthermore, introducing an alternative to the classical approach may be inspiring for others. Thus as for the learning phase, we consider that propagating waves generate peaks in the 2D-DFT curve in Fig. 9.b. Rather than evaluating the error on the propagating part of a wave at a given frequency, one will here use the error on the frequency of a periodic mode at a given wavenumber. To evaluate the attenuation accuracy, one could use the damping of the periodic peak, but will prefer the amplitude at the real frequency. As explained in section 3.4, only 2 wavenumbers have been selected as the 2D-DFT provides a satisfying shape. First the frequencies computed in the dispersion diagram of the reduced heterogeneous beam (black lines on Fig. 9.b) are compared to the reference frequencies of the full model (see Fig. 3.b). Fig. 10.a displays the relative error. For clarity, the markers and their position in the frequency axis are the same as in Fig. 3.b. It can be seen that the frequencies predicted by the reduced model are very satisfying as the relative error is always lower than 10 -6 %, which is negligible. For the amplitude, in Fig. 10.b, a few higher frequency errors are greater, around 7% but these correspond to poorly excited modes (blue/green resonances in Fig. 9.b). Providing an error estimation based on the number of wavenumbers selected in the learning phase (part 2.4) does not seem necessary, because only 2 are satisfying in this case. On all models already tested by the authors, it was never necessary to consider more than 3 wavenumbers. However, a more precise selection methodology is a perspective, to ensure that this methodology covers other application cases.

Extension to piecewise periodic or finite structures

Earlier sections focused on strictly periodic structures. However, the purpose of this work is more general. As periodicity is generally only achieved for a limited length, accounting for edge effects or transitions between multiple periodic parts is of interest for many applications.

Unlike in the strictly periodic case, evanescent waves are now present on edges of the finite model, or close to the interface between periodic sections. The learning phase leading to the selection of shapes and construction of a superelement assumes that such edge effects do not have a major influence within the area where the reduction is valid. As done in Ref. [START_REF] Mencik | A wave-based model reduction technique for the description of the dynamic behavior of periodic structures involving arbitrary-shaped substructures and large-sized finite element models[END_REF] for example, the proposed method consists in keeping the superelements obtained from the solutions for an infinite medium for the interior of the finite model, and avoiding reduction in areas close to edges where such effects can occur. This is very similar to the Saint Venant principle, which states that the response of a beam excited by a load, far enough from the application point, does not depend on how the load was applied, but only on the resulting force and moment. Two cases are adressed in this section: first a finite heterogenous beam, then a piecewise periodic railway track model.

Transfer computation on a finite heterogeneous beam

To demonstrate the validity of the procedure, a time domain transient is computed with a classical Newmark scheme for a finite structure subjected to a point load. A 20 cells finite beam with the heterogeneous properties discussed in section 3 is excited with a point force on the center of the left face (x = 0) and observed at the other extremity. The 18 interior cells are reduced with the same method as for the infinite beam. The 2 extreme cells at the finite boundaries are however non reduced, and coupled to the reduced cells with Multiple Point Constraints (MPC). A time Fourier transform is used to generate the response shown in Fig. 11. Rayleigh damping is used to obtain a case with reasonable levels at resonances. The maximum relative error on amplitude for the response shown is 3 10 -5 . Similarly, the modal frequencies are predicted with an error below 0.1% except for two modes of the local band which are still below 2%.

The learning phase is as earlier performed using 0.2 cell -1 (5 cells) and 0.01 cell -1 (100 cells) modes and a static correction. As expected for a multi-model reduction, the result is fairly independent of the choice of wavelengths. Static correction for the applied load is however necessary.

Piecewise periodic railway track model for transition zones

To provide a more challenging illustration, the case of the transition zone of Chauconin on the East-European High Speed Line between Paris and Strasbourg is considered. The transition is between a ballasted track (i.e. which contains a layer of a gravel material called ballast) and a slab track (i.e. ballastless). The track structure is modified with under-ballast mats in order to adjust the track stiffness between the ballasted and the slab track. As shown in Fig. 12.a, the geometry of the site is approximated by the combination of 5 periodic structures (1 for the ballasted track, 3 for the zone with an inclined under-ballast mat, 1 with an even mat of different stiffness). A portion of 105 m of track is modelled with finite elements. This model in practice has 2000 DOF for each interface, and 15,000 interior DOF for the reference cell in Fig. 12.b, resulting in more than 6,000,000 DOF in total.

The role of this interface is to ensure a smooth continuity between the 2 periodic structures A and B: it must be sufficiently large to allow the representation of local effects. The same consideration can be applied to finite structures: the edges where evanescent waves are also present can be dealt with a similar procedure. Here, the two extreme edges of the track are also modelled with an interface.

Defining a proper numerical strategy to evaluate the width needed for interfaces is still an open question. For the considered application, it is proposed to consider the static deflection of the track at the pass-by of a train. The width of the interfaces will thus be increased until the static response is perfectly smooth.

In Fig. 13, the reduced model response in the soil at the pass-by of a moving load at 80 m/s (288 km/h) for an interface width of 4 cells (2.4 m). The interfaces are shown as white areas. In this space/time representation, the interface fulfils its role as there is hardly no reflection.

In this more complex structure, each periodic section can be treated separately as presented before. For the reduction basis, three wavelengths of 2, 5 and 50 were found to be necessary for good accuracy in the superelement building [START_REF] Arlaud | Receptance of railway tracks at low frequency: numerical and experimental approaches[END_REF]. The only remaining points of attention are the interfaces, defined as sections of the structure between 2 periodic sections A and B, with strictly positive width. As done for cells, nodes inside this interface are split between those common with structure A, those common with structure B and inner ones. It is assumed that the superelements on each structure A and B are already constructed, with a reduced basis of the form [START_REF] Nashif | Vibration Damping, john wiley and sons Edition[END_REF]. Then one can keep the edge superelement shape T A LL from the structure A on the left, and similarly the superelement shape from the structure B on the right. This leads to an interface superelement using a basis of the form

T = ⎡ ⎢ ⎢ ⎢ ⎣ T A LL 0 0 -K -1 II T A LL T II -K -1 II T B LL 0 0 T B LL ⎤ ⎥ ⎥ ⎥ ⎦ . ( 31
)
The first difference with the periodic case is that left and right shapes T A LL and T B LL are distinct. The second is in the handling of interior motion T II . Since the number of interior DOF may be quite large, a reduction procedure is necessary. The proposed strategy is to use the Craig-Bampton approach [START_REF] Craig | A review of time-domain and frequency domain component mode synthesis methods[END_REF]. Indeed the resulting cost is acceptable because the number of interface DOF (that is to say the number of independent columns in T LL , or the number of generalized interface DOF in the superelements) is small. Now that the interfaces are validated, this track model can be used to compute forced response of the track, at the pass-by of trains. In Table 3, some data are given to illustrate the interest of the reduction. The computation of the reduction basis takes 23min and requires 800 MB of memory. Then, each transient computation only takes 3min and 106 MB, to be compared to a full 3D computation, that would take more than 3 h (indirectly estimated time) and 96 GB of memory to store the output. Furthermore, the interest of this reduction is even higher for a parametric study of a structure: if the same basis can be kept for each test of the parametric space, the gain is really dramatic.

Conclusion

This paper proposed a superelement building strategy using a Ritz-Galerkin procedure to construct a reduced model, in standard second order form, starting from any learning subspace. Compared to the independently derived but very similar procedure [START_REF] Boukadia | A Bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF], the proposed learning phase only uses periodic solves for the full model, thus bypassing the need for a costly Craig-Bampton reduction. Furthermore, the reference to multi-model reduction [START_REF] Balmes | Parametric families of reduced finite element models. Theory and applications[END_REF] also motivates the fact that no more than two wavelengths are needed, thus alleviating the need for a costly wavenumber sampling strategy and the use of WFE at the maximum frequency of interest. The only example known to the authors where three wavelengths of 2, 5 and 50 were found to be necessary a track structure [START_REF] Arlaud | Receptance of railway tracks at low frequency: numerical and experimental approaches[END_REF], where reproduction of the so called pin-pin mode requires the addition of the wavelength corresponding to 2 cells. Automating the evaluation of errors for intermediate wavelengths using a residue iteration process [START_REF] Bobillot | Iterative computation of modal sensitivities[END_REF] is a clear perspective, as is the development of more systematic criteria for edge truncation using the singular value decomposition. To account for applied loads, viscoelastic damping [START_REF] Hammami | Numerical design and test on an assembled structure of a bolted joint with viscoelastic damping[END_REF], and even possibly radiation in infinite media, standard parametric model reduction tools are applicable. Finally, obtaining a reduced cell model that can be used in both frequency and time domains seems a definite advantage over earlier work.

A random beam example combined difficulties of presenting both local and global modes, Bragg bandgaps and wave-mode crossing. The model reduced with only two wavelengths, a small (5 cells) and a long (100 cells), gave very accurate predictions of frequencies and forced response of infinite and finite structures. The illustrations sought to highlight the fact that computing forced responses in the frequency/wavenumber domain is a useful alternative to the classical consideration of complex eigen wavenumbers and a wave shape tracking procedure was introduced. The second industrial case of a transition zone in a railway track was used to discuss extension to piecewise periodic models and performance issues. Here both the computation time and the required memory were lowered by 1-2 orders of magnitude. This example however left open questions about the size of non reduced areas at interfaces, which will need to be addressed.

A final perspective would be an extension cases with periodicity in more than one direction.

Fig. 1 .

 1 Fig. 1. Transfer of a beam in traction/compression in the frequency/wavenumber domain (the colour scale is the logarithm of the amplitude of the 2D-DFT). The theoretical dispersion curve is superimposed in black. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 2 .

 2 Fig. 2. Isovalues of density in a periodic random waveguide (80% of the maximal density in blue and 99% in red). See section 3.1 for properties. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

  the spatial distribution of loads b is replicated for the spatial real and imaginary contributions of the applied forces F. Inverting the constrained dynamic stiffness matrix Z c = T T c [

Fig. 3 .

 3 Fig. 3. Comparison of dispersion diagrams of infinite periodic beams.

Fig. 4 .

 4 Fig. 4. Examples of bending and pressure modes of a heterogeneous beam.

Fig. 5 .

 5 Fig. 5. Interaction of the pressure mode with the local mode.

Fig. 6 .

 6 Fig. 6. Views of the non-reduced 2D-DFT computation of the forced response in bending of an infinite periodic heterogeneous beam. The black lines correspond to the dispersion curves presented in Fig. 3.b.

Fig. 7 .

 7 Fig. 7. Strain energy rate for 20 samples of randomly-fluctuating heterogeneous beams. The mean value is shown in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 8 .Fig. 9 .

 89 Fig. 8. Comparaison between the FFT of the forced response with model reduction (in plain orange line) and without reduction (in dashed blue line). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 10 .

 10 Fig. 10. Evaluation of the error for the reduction of the heterogeneous beam model.

Fig. 11 .

 11 Fig. 11. Transfer for a 20 cells heterogeneous beam. Plain lines represents the reduced model and triangle markers the modes of the structure.

Fig. 12 .

 12 Fig. 12. FE model of the transition zone of Chauconin [37].

Fig. 13 .

 13 Fig. 13. Space/Time diagram of the vertical deflection of the reduced model of the track for a moving load at 80 m/s.

Table 1

 1 Properties of the heterogeneous beam.

	Property	Mean value	Standard deviation
	Bulk modulus K	120 MPa	361MPa
	Shear modulus G	24.5MPa	73.4MPa
	Density 𝜌	1700kg/m 3	5100kg/m 3
	Correlation length 𝜆 corr	0.1m	-
	Length reference cell x 0	0.6m	-
	Section surface A	0.25 × 0.25m 2	-
	Minimal mesh size h	0.04m	-

Table 2

 2 Comparison of DOF number between the initial and reduced model.

	DOF number	Initial model	Reduced model	Gain
	Interface (L/R)	2883	49	58×
	Interior (I)	141,267	41	3445×
	Total	147,033	139	1057×

Table 3

 3 Comparaison of numerical parameters in a computation of transition zone.

		Initial model	Reduced model	Gain
	DOF number (total)	6,000,000	7125	842×
	Computation time	≫ 3 h (each)	23min (red.) + 3min (each)	≫ 7×
	Memory required

>96 GB (each) 800 MB (red.) + 106 MB (each) 108×