C Droz

C Zhou

M N Ichchou

J.-P Lainé

A hybrid wave-mode formulation for the vibro-acoustic analysis of 2D periodic structures

Keywords: Wave, finite element, periodic structure, reduced model, stiffened plate 1. Literature survey

In the framework of vibrational analysis of 2D periodic waveguides, Floquet-Bloch theorem is widely applied for the determination of wave dispersion characteristics. In this context, the Wave Finite Element Method (WFEM) combines Periodic Structure Theory (PST) with standard FE packages, enabling wave dispersion analysis of waveguides involving structurally realistic unit-cells. For such applications, the computational efficiency of the WFEM, depends on the choice of the formulation and can lead to numerical issues, worsen by extensive computational cost. This paper presents a coupled wave-mode approach for the determination of wave dispersion characteristics in structurally advanced periodic structures. It combines two scales of model order reduction. At the unit-cell's scale, Component Mode Synthesis (CMS) provides the displacement field associated with local resonances of the periodic structure, while the free wave propagation is considered using a spectral problem projection on a reduced set of shape functions associated with propagating waves, thus providing considerable reduction of the computational cost. An application is provided for a bi-directionally stiffened panel and the influence of reduction parameters is discussed, as well as the robustness of the numerical results.

higher orders of the theory including kinematic assumptions. However it requires high-order derivatives leading to an increased mathematical effort to obtain analytical solutions [START_REF] Sutyrin | Derivation of plate theory accounting asymptotically correct shear deformation[END_REF]. In [START_REF] Chronopoulos | Modeling the response of composite panels by a dynamic stiffness approach[END_REF], the author studied the dynamic response of a sandwich panel using a dynamic stiffness approach and several Equivalent Single Layer (ESL) theories using an homogenized orthotropic honeycomb core. The use of modal reduction can overcome some FEM methods limitations between low-and mid-frequency range, but reaches its limits for composite materials which exhibit high modal densities.

In the high frequency range, the mechanical energy mainly remains in the periodic cell and the diffusion between the structure subsystems is weak. Then, the vibro-acoustic behaviour is governed by the periodic cells dynamics and can be determined from a single periodic cell [START_REF] Ruzzene | Wave beaming effects in two-dimensional cellular structures[END_REF] using Bloch theorem. The averaged dynamic response of the structure under such excitation can be estimated using Statistical Energy Analysis (SEA) [START_REF] Lyon | Theory and application of statistical energy analysis[END_REF]. However, SEA method suffers significant limitations in the medium frequency range, considering it cannot provide the response of the longwavelength components.

In the medium-frequency range, the vibratory energy is partially localized in the periodic cell and propagates through the plate. In the last decade, numerical methods based on multi-modal propagation received a lot of attention. Among them, the Wave Finite Element (WFE) method [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF][START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF][START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF] is an effective tool for providing the wave dispersion characteristics of a wide range of waveguides. The main feature of the WFE method is to combine periodic structure theory with commercial finite element packages. Numerous applications were proposed for one-dimensional homogeneous and periodic waveguides [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF][START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF][START_REF] Ichchou | Multi-mode wave propagation in ribbed plates: Part i, wavenumber-space characteristics; part ii, predictions and comparisons[END_REF][START_REF] Houillon | Wave motion in thin walled structures[END_REF][START_REF] Waki | Free and forced vibrations of a tyre using a wave/finite element approach[END_REF][START_REF] Chronopoulos | Predicting the broadband response of a layered cone-cynlinder-cone shell[END_REF]. A formulation was introduced for two-dimensional waveguides [START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF] and applied to periodic structures in [START_REF] Collet | Floquet-bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems[END_REF][START_REF] Orrenius | Wave modelling in predictive vibro-acoustics: Applications to rail vehicles and aircraft[END_REF][START_REF] Renno | Calculating the forced response of two-dimensional homogeneous media using the wave and finite element method[END_REF]. Still, the above-mentioned methods are based on Floquet -Bloch theorem and face some drawbacks due to high computational costs when refined FEM meshes are required. Some numerical issues were investigated for example in [START_REF] Waki | Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides[END_REF] and a reduction method to compute the response of one-dimensional waveguides under harmonic excitations was proposed in [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] Renno | On the forced response of waveguides using the wave finite element method[END_REF]. However, there is still a lack of studies in the literature concerning model order reduction strategies for the WFE method in the case of 2D periodic waveguides.

As a consequence, the study is often restrained to lower frequencies due to coarsely meshed models or structural approximations. In [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF], the authors proposed a reduced strategy for studying propagating waves dispersion characteristics in structurally advanced or composite waveguides. Several formulations using Bloch's theorem can be found in the literature for 2D periodic waveguides [START_REF] Langley | A note on the force boundary conditions for two-dimensional periodic structures with corner freedoms[END_REF][START_REF] Akrout | Comportement dynamique déterministe et large bande des structures guidées[END_REF][START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF]. These methods require the resolution of an eigenproblem depending on the spatial variables used in the periodicity conditions, and possible degrees of freedom (DOF) condensation. One can discriminate these formulations between "unknown-frequency" and "unknown-wavenumber" resolutions. These different approaches are discussed in section 2. Since these formulations generally involve extensive computation when an important degree of discretization is required, the inners degrees of freedom condensation can benefit Component Mode Synthesis (CMS). In [START_REF] Cotoni | A statistical energy analysis subsystem formulation using finite element and periodic structure theory[END_REF], the authors successfully applied the CMS to a unknownfrequency WFE method. However, this formulation requires to choose arbitrary real phase constants and is not suited for damped structures. In order to provide a suitable numerical tool for wave analysis in structurally advanced periodic waveguides, this paper focuses on the prediction of propagating waves dispersion characteristics using reduced computation processes. The aim is to provide a suitable numerical tool for computing wave dispersion in two-dimensional periodic systems using a CMS approach for the inners cell DOF condensation, then reducing the spectral eigenproblem projecting the cell's cross-sectional transfer matrix on a reduced set of shape functions. It allows the propagating waves to be determined accurately by solving smaller eigenproblems, and enables low-cost broadband analysis of finely meshed two-dimensional periodic waveguides.

This paper is organized as follows. In section 2, some WFE formulations for 2D periodic waveguides are reviewed. The purpose of this section is to classify the methods, highlight their numerical issues and motivate the chosen "unknown-k" formulation. We also introduce the proposed developments presented in the following section. Section 3 describes the proposed reduction strategy. A CMS approach based on mode selection is first developed for the dynamic condensation required in the proposed WFE formulation. A modal analysis is performed on the fixed boundary unit cell; the main contributing modes of the periodic cell are retained to replace the inner cell's DOFs; a dynamic condensation of the modal participations allows to use classical 2D periodic WFE formulation used in [START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF]. Then a reduced formulation based on wave interpolation is provided to solve the spectral problem on a broadband frequency range. A subset of solutions is computed at waves cut-on frequencies; a shape basis is determined using a reduced number of propagating waves; an iterative process based on eigenvectors orthogonality is proposed to improve the reduced basis numerical stability; spectral problem is then projected on the reduced-size basis and easily solved for the whole frequency band. Numerical application is brought in section 4. Orthogonally stiffened panels are considered; classical WFE solutions are compared the proposed reduced formulation; the solutions sensitivity to both modal and wave selection is investigated and the computational gains are discussed.

Review of the WFE formulations for 2D periodic waveguides

Application of the periodicity relations

In this section, an overview of the WFE formulation is given for the general 2D periodic waveguide. These structures consists in a set of identical and regularly connected elements. A periodic cell is shown in Fig. 1, the degrees of freedom q are partitioned into edges (q B , q T , q L , q R ), corners (q 1 , q 2 , q 3 , q 4 ) and inner q I DOFs. The vector of nodal forces is defined similarly excepted for the internal DOFs where no external forces are applied: f I = 0. It is assumed that opposite edges and the four corners have the same number of DOFs. Apart from these conditions, the FEM discretization is arbitrary inside the periodic element. The relations between the element displacements can be expressed using Bloch's theorem:

q R = λ x q L , q T = λ y q B . (1)
for the edges, and

q 2 = λ x q 1 , q 3 = λ y q 1 , q 4 = λ x λ y q 1 . (2)
for the corners. The forces relations can be expressed similarly:

f R = -λ x f L , f T = -λ y f B . (3)
for the edges forces, and

f 2 = -λ x f 1 , f 3 = -λ y f 1 , f 4 = λ x λ y f 1 . (4)
for the corners forces. The propagation constants λ x = e -jk x d x and λ y = e -jk y d y in x and y direction, involve the wavenumbers k x and k y and the element size (d x × d y ). The values k x d x and k y d y are also called phase constants.

Similarly, the forces equilibrium leads to the following conditions

f B + λ -1 y f T = 0 , f L + λ -1 x f R = 0 (5)
for the edges, and

f 1 λ x λ y + f 2 λ y + f 3 λ x + f 4 = 0 (6)
for the corners. The cell's equation of motion can be used to form the spectral problem involving the wavenumbers k x , k y and the frequency. Denoting K and M the stiffness and mass matrices and ω the circular frequency, the governing equation in the element is written, assuming time-harmonic excitation:

(K -ω 2 M)q = f (7)
Different methods were provided to form a well-conditioned eigenproblem for the determination of wave dispersion characteristics in 2D periodic waveguides. The choice of a method strongly influences the problem size, the results usability and model order reduction possibilities. One can discriminate between a given (k x , k y ) formulation (see section 2.2) and a given (ω, k x ) formulation described in section 2.3.

Given (k x , k y ) formulation

The cell's DOFs can be ordered [START_REF] Langley | A note on the force boundary conditions for two-dimensional periodic structures with corner freedoms[END_REF] as follows:

     I 0 0 0 0 I 0 0 0 Iλ y 0 0 0 0 I 0 0 0 Iλ x 0 0 0 0 I 0 0 0 Iλ x 0 0 0 Iλ y 0 0 0 Iλ x λ y                                                       q I q B q L q 1               (9)
Substituting Eq.( 9) in the governing equation Eq.( 7) and pre-multiplying by R T yields:

R T (K -ω 2 M)Rq = 0 (10)
In practice, the eigenvalues ω i are found for each couple (λ x , λ y ) by specifying a range of values for k x and k y . Then, the propagating solutions can be sorted in order to plot phase constants surfaces illustrated in Fig. 2. One can note than Eq.( 10) requires the computation of a sizeable eigenproblem since the inner DOFs q I are involved.

A CMS approach was proposed in [START_REF] Cotoni | A statistical energy analysis subsystem formulation using finite element and periodic structure theory[END_REF] to reduced the computational expense associated with the phase constant surfaces determination described above. The method uses a subset fixed-interface modes of the periodic element to replace the interior DOFs. The inner DOFs q I are written as:

q I = Φp Φ -K -1 II K Ib q b ( 11
)
where Φ is a subset of the clamped boundary local modes, p Φ contains the associated modal participations, while the term K -1 II K Ib q b stands for the constraint modes where subscript b indicates the periodic element boundary DOFs. The number of modes retained is significantly lower than the number of inner DOFs, the reduced state vector q in Eq. [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF] being written: Then, the associated computation time is significantly reduced for each couple of (λ x , λ y ) specified. The same authors applied this procedure and obtained significant CPU time reduction in comparison with the standard approach. Notwithstanding, this formulation suffers some drawbacks. Firstly, the specified wavenumbers (k x , k y ) have to be assumed real hence restraining the study to undamped systems, since no information is given on the imaginary parts which are associated with propagating waves spatial attenuations. Secondly, the solutions provided here are not fitted for a further use in forced response analyses, where the frequency is a given parameter. Later formulations of the spectral problem were provided in section 2.3, allowing a dynamic condensation of the element inner DOFs.

q =               p Φ q B q L q 1               (12)

Given (ω, k x ) formulation

In [START_REF] Akrout | Comportement dynamique déterministe et large bande des structures guidées[END_REF], the author proceeds a dynamic condensation of the inner DOFs, and gives the explicit coefficients X(ω, λ x ), Y(ω, λ x ), Z(ω, λ x ) so that the spectral problem Eq.10 can be written on the quadratic form:

(λ y X + Y + λ -1 y Z)           q L q B q 1           = 0 (13)
A given (ω, k y ) formulation can be obtained similarly by inverting the roles of k y and k x . Further investigations regarding the resolution of linear, quadratic, polynomial or transcendental eigenproblems were carried out in [START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF]. The complex solutions λ y are then sorted in order to retain mainly propagating waves. If the wavenumber k y is purely imaginary the wave is evanescent, if k y is complex the wave is propagating and | (k y )| is the spatial attenuation. Therefore, the wave dispersion characteristics can be determined on a given-frequency form, enabling the computation of forced responses from a wave-based method. An example of the results obtained by a given-(ω, k y ) formulation is described in Fig. 3 for a thick-layered plate. Since ω is given in this formulation, the group velocities can be determined to evaluate the waves directivity in 2D waveguides.

Group velocities 400kHz

Figure 3: Group velocities of a multi-layered orthotropic plate at 400 kHz, obtained using given-(ω, k y ) WFE formulation and plotted in polar coordinates.

Yet, the latter formulation requires a dynamic condensation of the inner DOFs for each frequency step. It may lead to extensive computation time when the periodic cell is investigated in higher frequencies, or contains a complex geometry and needs refined FEM. Besides, the condensed eigenproblem Eq.( 13) still involves numerous DOFs and remains a major drawback for analysing structurally advanced periodic waveguides. In the following section a hybrid formulation is proposed which combines the use of a modal reduction approach for the dynamic condensation with a wave shape frequency-interpolation.

Proposed model order reduction strategy

It should be noticed that the proposed formulation is a two-scale approach [START_REF] Langley | Some perspectives on wave-mode duality[END_REF], using hybrid CMS formulation in a wave approach context. The CMS, described in section 3.1 is performed at the cell's scale (mesoscopic scale), while a frequency-interpolation method based on propagating waves determination is detailed in section 3.2 to section 3.4. As a convention, the symbol Φ is used for modes and Ψ is used for waves.

A Craig-Bampton reduction for the cell's inner DOFs

The CMS procedure provides a significant reduction of the number of inner DOFs involved by replacing physical displacements by the cell's local modes. To this aim, the displacement vector q defined Eq.( 7) is partitioned into the inner cell's displacements q I , and the boundaries displacements:

q b = q T B q T T q T L q T R q T 1 q T 2 q T 3 q T 4 T (14)
Since there are no external forces applied on the cell's inner DOFs, f I = 0, Eq.( 7) can then be written:

K b K bI K Ib K I -ω 2 M b M bI M Ib M I q b q I = f b 0 (15)
Therefore, a reduced set of coordinates is required to replace the physical DOFs q I . This reduced Ritz basis involves vectors associated with static boundary modes Φ b and selected component modes Φ C so that the displacements can be written:

q b q I = I 0 Φ b Φ C q b P C = B q b P C (16)
where the vectors Φ b are obtained by static condensation

Φ b = K -1 I K Ib (17)
the component modes Φ C are selected among the eigenvectors Φ I of Eq.( 18)

K I -ω 2 M II Φ I = 0 (18)
and the inner DOFs displacements are replaced with a reduced set of modal participations P C . In Craig-Bampton procedure, modal selection is based on the lower resonance frequencies. This method can be extended in a wave approach context since the aim here is to capture the local deformed shape of the periodic cell. The underlying hypothesis is that when a wave propagates through the periodic waveguide, the displacements inside a unit-cell can be expanded on a subset of stationary modes. Therefore, one can choose [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] modes into the frequency range [0, 2 ω max ],

where ω max is the maximum frequency of interest for the wave dispersion analysis. Then, the stiffness and mass matrices can be written in the reduced set of coordinates using the projection matrix B Φ defined Eq.( 16):

K = B T Φ KB Φ , M = B T Φ MB Φ (19)
It should be noted that free modes could also be used for the CMS procedure. Yet, they do not improve the reduced model accuracy and require additional projection of the boundary degrees of freedom. As mentioned in section 2.3, the reduced matrices can be used to perform dynamic condensation in order to obtain a spectral problem of size

q T 1 q T L q T B T
. However the size of this state vector remains a major drawback when a large frequency band is considered. A wave-based reduced formulation for broadband analysis of 2D homogeneous or periodic waveguides is developed in the following.

Determination of a reduced frequency sample for the wave shape interpolation

The principle of the model order reduction (ROM) strategy proposed here is to reduce Eq.( 13) using a subset of the propagating waves solutions. The main assumption is that a waveguide cross-sectional deformed shape for a given wave type at a given frequency can be used to describe the motion of the periodic cell for the same wave in a large frequency range around this given frequency. It will be considered further in this paper for 2D periodic problems. In order to identify a suitable reduced subset Ω of the frequency range Ω, one can consider the specific behaviour of propagating waves in periodic waveguides. Indeed, in such structures the real wavenumbers remain between k = 0 and k = π/d (see [START_REF] Mead | Free wave propagation in two-dimensional periodic plates[END_REF]). The wave shapes are then retained at the mid-aliasing frequencies, defined by the wavenumber: k = π/(2d). Therefore, the frequency subset can be determined by using given (k x , k y ) formulation. From the reduced matrices defined Eq.( 19), the dispersion relation defined Eq.( 10) can be written using reduced modal participations:

Λ L (λ x , λ y ) K -ω 2 M Λ H L (λ x , λ y )               q 1 q L q B P C               = 0 ( 20
)
where Λ L is the definite positive Hermitian matrix detailed in Appendix A. In this given (k x , k y ) formulation, the propagating solutions are associated with the real eigenvalues of Eq.( 20). These solutions form the frequency subset Ω used in the following to span the reduced wave basis.

Resolution of the spectral problem on the frequency subset ∆ω

Considering matrices K and M defined section 3.1 and the given frequency ω ∈ Ω, the reduced dynamic stiffness defined as D = K -ω2 M can be ordered as follows:

Db DbI DIb DI q b P C = f b 0 (21)
Then, the modal participations P C can be condensed on the cell's boundaries so that:

D = Db -DbI D-1 I DIb (22)
The dynamic stiffness can then be ordered as in Eq.(B.1), so that the spectral problem can be written in terms of this condensed matrix:

Λ L (λ x , λ y ) D(ω) Λ R (λ x , λ y )           q 1 q L q B           = 0 (23)
In this case the problem is not linear since λ x is the unknown and (ω, λ y ) are given. Therefore, a quadratic eigenproblem in λ x can be formed using Eqs.(23,B.1):

          λ x           A 11 A 1L A 1B A L1 A LL A LB A B1 A BL A BB           +           B 11 B 1L B 1B B L1 B LL B LB B B1 B BL B BB           + 1 λ x           C 11 C 1L C 1B C L1 C LL C LB C B1 C BL C BB                               q 1 q L q B           = 0 (24)
where the coefficients A ij , B ij and C ij are detailed Appendix B. Then denoting A, B, C these matrices, the quadratic eigenproblem Eq.( 24) can be solved using the following linearisation:

-C 0 0 I q λ x q = λ x B A I 0 q λ x q (25)
For the sake of clarity, the spectral problem Eq.( 24) will be denoted: S λ x , λ y , ω Ψ = 0, where the solution Ψ is the wave shape on the condensed DOFs defined as:

Ψ =           q 1 q L q B           (26)
The solutions of Eq.( 24) are computed on the mid-aliasing frequency subset Ω defined in section 3.2. Accordingly, an appropriate reduced basis for the ROM can be spanned using solutions of the following system of equations:

                         S λ x , λ y , ω1 Ψ = 0 • • • S λ x , λ y , ωj Ψ = 0 ; ωj ∈ Ω • • • S λ x , λ y , ωn Ψ = 0 ( 27
)
where n is the number of elements in the subset Ω. For each given ωj , discrete values of k y are chosen in the range [-π/d y , π/d y ]. The values of λ y are derived from these choices of k y , with λ y = e -jd y k y . It is noteworthy that the CPU time required to solve the n eigenproblems Eq.( 27) is considerably lower than for solving Eq.( 24) on the whole frequency domain Ω. The reduced basis can be built using these eigenvectors. However, each line of Eq.( 27) yields 2 × N solutions, where N is the size of

Ψ = q T 1 , q T L , q T B T
, corresponding to the number of DOFs in the state vector. Retaining all of the 2 × N solutions Ψ at the n frequencies would lead to an extensive wave shape basis, while the only solutions of interest are associated with propagating waves. This issue is discussed below.

Building a reduced wave shape basis

The idea is to select propagating waves using their dispersion characteristics in order to build a reduced wave basis. A complex wavenumber k x = j ln(λ x ) is associated with each eigenvector Ψ, and its imaginary part (k) is the wave attenuation. The wave is evanescent when the wavenumber is purely imaginary, while real positive wavenumbers are associated with the least decaying positive-going waves. Yet, damped structures can produce highly decaying waves and the determination of propagating solutions can be delicate. In this study, a propagating wave is described by the following conditions:

             0 ≤ (k x ) ≤ π /d x Positive-going waves (k x ) ≤ 0 Physical solutions | (k x )| ≤ α (k x )
Least-decaying waves, α 0.01 [START_REF] Langley | Some perspectives on wave-mode duality[END_REF] We denote Ψ ( j) i the eigenvector associated with the i-th propagating eigenvector at the frequency ω j . Therefore, the solution subspace can be written:

Ψ = {Ψ ( j) 1 , . . . , Ψ ( j) p( j) } 1≤ j≤n (29)
where p( j) stands for the number of propagating waves retained at frequency ω j . We denote p = n i=1 p(i) the total number of vectors in Ψ. For a sake of clarity the j-th column of matrix Ψ is a vector, denoted Ψ( j). In practice, this set may involve redundant vectors and lead to an unnecessary large wave basis when a wave shape remains the same at several frequencies in Ω. Considering modal assurance criteria (MAC) as an indicator of correlation between two complex eigenvectors, the orthogonality relation ⊥ ε can be written:

ψ 1 ⊥ ε ψ 2 ⇔ (ψ T 1 ψ2 ) (ψ T 2 ψ1 ) (ψ T 1 ψ1 ) (ψ T 2 ψ2 ) ≤ ε (30)
The criteria ε ∈ [0 , 1] is chosen close to 1 and describes the tolerance for the orthogonality relation. This definition is extended to define the orthogonality between a vector and a set ψ 1 ⊥ ε {ψ}, when ψ 1 is orthogonal to all the vectors in {ψ}. An optimisation procedure is required to reduce the number of vectors in the wave basis. In order to create this subset Ψ of Ψ by elimination of correlated wave shapes and build a reduced wave basis providing good numerical stability, we propose the algorithm described Eq.( 31). First, the vector subset Ψ is built by recurrence using ⊥ ε orthogonality. Note that the choice of ε influences the size of the reduced basis and the ROM efficiency. Then, a simple way to improve numerical stability is to normalize the remaining vectors and dissociate their real and imaginary parts. Thus, the wave basis Γ is purely real and wellconditioned.

                                                                   Define Ψ(1) = Ψ(1) for j = 1 to j = p -1                  if Ψ( j + 1) ⊥ ε Ψ do Ψ = Ψ , Ψ( j + 1)
end if end for define the normalized Γ i = Ψ(i)

Ψ(i) dissociate complex vectors: Γ i into: (Γ i ) (Γ i ) build the wave basis Γ = [Γ 1 , . . . , Γ r ] (31)
The state vector can be written in terms of wave participations:

          q 1 q L q B           =             Γ (1) 0 0 0 Γ (L) 0 0 0 Γ (B)                              = BΨ w =                                    Γ (1) Γ (2) 0 Γ (3) Γ (4) Γ (L) Γ (B) 0 Γ (R) Γ (T )                                   
where 4) , Γ (L) = Γ (R) and Γ (B) = Γ (T ) . It is noteworthy that wave participations are defined separately for each of the periodic cell's components. Geometrically, it means that three different shape basis are defined for the corners and both pairs of opposite edges, and each of these components are provided with independent coefficients. Since the coefficients in Appendix B are linear with respect to matrix D, one can build reduced dynamic stiffness BT Ψ D BΨ . Then denoting:

Γ (1) = Γ (2) = Γ (3) = Γ (
B Ψ = BΨ 0 0 I (34)
and using Eq.( 19), one can write reduced mass and stiffness matrices leading to a reduced CPU computation for the dynamic condensation Eq.( 22) and the spectral eigenproblem Eq.( 24):

KΨ = B T Ψ B T Φ K B Φ B Ψ , MΨ = B T Ψ B T Φ M B Φ B Ψ (35)
To summarize, a modal basis is defined using Eqs. [START_REF] Chronopoulos | Predicting the broadband response of a layered cone-cynlinder-cone shell[END_REF][START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF][START_REF] Collet | Floquet-bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems[END_REF][START_REF] Orrenius | Wave modelling in predictive vibro-acoustics: Applications to rail vehicles and aircraft[END_REF] in order to replace the inner DOFs of the periodic cell with a reduced set modal participations so that the dynamic condensation Eq.( 22) can be easily done on the frequency range Ω. Then, the spectral problem Eq.( 24) is solved on a small frequency subset Ω, defined using Eq. [START_REF] Renno | Calculating the forced response of two-dimensional homogeneous media using the wave and finite element method[END_REF]. The solutions of system Eq.( 27) associated with positive-going propagating waves are selected using the conditions Eq.( 28), and a reduced set of eigenvectors is built using the algorithm provided Eq. [START_REF] Orrenius | Calculation of wave propagation in rib-stiffened plate structures[END_REF]. Therefore, the wave participations are defined using Eq.( 33) and reduced mass and stiffness matrices are provided using Eqs. (34,35). Finally, a reduced size spectral problem Eq.( 24) can be solved on the whole frequency range Ω, while the physical solutions can be derived using Eqs. [START_REF] Chronopoulos | Predicting the broadband response of a layered cone-cynlinder-cone shell[END_REF]33).

Numerical results: orthogonally stiffened panels

Problem description

In this section, the formulation developed in the paper is used to provide wave dispersion characteristics (wavenumbers, group velocities, wave modes) in bidirectionally stiffened panels. The wave propagation in such structures is extensively studied in the literature [START_REF] Orrenius | Calculation of wave propagation in rib-stiffened plate structures[END_REF][START_REF] Legault | Vibro-acoustic response of orthogonally stiffened panels: the effects of finite dimensions[END_REF]. The mode selection and wave-based interpolation parameters introduced in section 3 are used here to study the stability and accuracy of the proposed method. A flat panel reinforced by stiffeners with rectangular cross-sections is considered in this example. The unit-cell of the periodic plate is shown in Fig. 4. The skin is an aluminium plate of thickness h = 0.8 mm and the stiffeners are rectangular aluminium beams of width l 1 = l 2 = 10 mm and height h 1 = 20 mm and h 2 = 10 mm. A structural loss factor η = 0.5% is introduced and the spacing between the stiffeners is L 1 = 0.08 m in the y-direction and L 2 = 0.12 m in the x-direction. The The FE model built using ANSYS 14.0 Software involves 600 4-nodes SHELL63 (DKT) elements for the skins, and 50 rectangular BEAM44 elements with cubic interpolation for the stiffeners. Therefore, the periodic cell has a total of 3906 DOFs, involving 3306 DOFs inside the periodic cell and 600 DOFs on the cell's boundaries. Since the stiffened panel exhibits two localized resonance at 840 Hz and 2040 Hz, the frequency range studied in this application is [0 -2500] Hz.

h l 1 l 2 h 1 h 2 L 1 L 2 z y x

Reduced WFEM solution

The proposed reduction strategy is applied to the model described above and compared to the standard formulation of the WFEM. Five stationary modes are employed to describe the inner cell's displacements, according to the CMS procedure described section 3.1 and 15 waves are determined using the procedure Eq.( 31) to build the reduced solution subspace. Real wavenumbers associated with the propagating waves in the x-direction are presented in Fig. 5 using the standard WFEM and the proposed formulation. Vertical lines, where the full eigenproblems are solved denote the mid-aliasing frequencies. Note that all of intersections between the propagating waves and the mid-aliasing frequency are correctly predicted. A very good correlation is observed between the standard WFEM formulation and the reduced model, involving modal condensation and wave-based interpolation. The maximal error on the wavenumbers is err = 0.8% for the flexural wave. The modal basis is built using a subset of stationary modes between 0 and 2500 Hz, and the wave basis is computed using a correlation criteria ε = 0.99 described Eq.( 30), and provides a reduction of the quadratic eigenproblem from 600 boundary DOF's to 120 wave participations.

WFEM

Proposed formulation Wave basis

Figure 5: Real wavenumbers in the y-direction (for k x = 0) using the standard WFE method and the proposed CMS -wave interpolation strategy A significant reduction of the computational cost is observed for the proposed method. CPU time using Matlab 7.14 running on Linux with two AMD 64×2 dual core 6000+ processors are given in Table 1. The CMS procedure requires an additional modal analysis, but provides an important reduction (95.5%) of the time required for the dynamic condensation. On the other hand, the wave based reduction strategy also provides a considerable reduction (95.7%) of the time required to solve the quadratic eigenproblems, in comparison with the standard formulation involving the edges unit-cell nodal displacements. Therefore, the proposed formulation provides a total of 80% reduction of the computational cost. Noteworthy that a higher frequency sampling will increase the wave-based reduction efficiency. The aforementioned reduction strategy is then applied in the 2D k-space, and compared to the given -(ω, k y ) standard WFE formulation, presented in section 2.3. The wavenumbers k x and k y are shown in Fig. 6 at 1500 Hz, and exhibit a very good correlation between the standard WFEM and the proposed method. The flexural ('o'), shear (' ') and longitudinal waves ('+') are accurately predicted by the reduced formulation ('O'). Besides, in Fig. 7 the localized resonances, associated with important variations of the directivity and appearance of angular bandgaps for the flexural wave, are predicted as well at 2400 Hz. At this frequency, the reduced model exhibits a good accuracy, although k x is given for all the values of the first Brillouin zone. Therefore, the k-space can be accurately determined in broadband frequency range and provide the required 2D wave dispersion characteristics of the periodic structure.

Sensitivity of the CMS to the modal basis

The accuracy of the proposed ROM strategy relies essentially on two parameters: the number of modes retained for the CMS procedure and the number of wave shapes used to span the solution subspace. For the CMS procedure, the fixed-interface modes are usually retained between 0 and two times the maximal frequency. An insufficient number of modes is expected to produce significant errors in the determination of high-order wave shapes, since the displacement field inside the periodic cell can no longer be described. In Fig. 8, the dispersion curves are compared are shown in Fig. 9. The second basis (CMS 1500 Hz, '×') requires the two modes below 1500 Hz and accurately describes the first band-gap, but fails to predict the second band-gap at 2000 Hz. Similarly, the third basis (CMS 0 Hz, '+') is limited to the the static constrained modes. Therefore the model is then limited to the low frequencies and fails to predict the band-gaps and high-order flexural waves. However it can be noticed that the longitudinal and shear waves are in good agreement in the three configurations. It can be explained since their shapes remain the same on the frequency range and can be described using the static constrained modes. The shape of the longitudinal wave at 1000 Hz and the shear wave at 750 Hz are shown respectively in Fig. 10

Sensitivity of the method to wave-based reduction

The aforementioned wave-based interpolation involves a parameter ε introduced in Eq. [START_REF] Mead | Free wave propagation in two-dimensional periodic plates[END_REF], to build the reduced set of shape vectors using a chosen subspace of the spectral problem solutions. The effect of this wave basis refinement is shown in Fig.12. In the proposed example, a choice of ε = 0.85 would considerably reduce the number of vectors retained in the shape basis but leads to an important error for the second and third orders of the flexural wave. The third order flexural wave shape is shown at 2500 Hz in Fig. 13. Noteworthy, this displacement field cannot be described using the first and second orders flexural wave shapes. Oppositely, ε = 1 would mean that all of the propagating wave solutions are retained at each interpolation frequency. Therefore, it can be seen that the wave correlation parameter strongly influences the wave selection and the accuracy of the reduced model.

Concluding remarks

The effectiveness of wave dispersion analysis for predicting the vibro-acoustic behaviour of large-scaled 2D waveguides led to numerous research work in the last decade. The Wave Finite Element Method (WFEM) is suitable tool for such purposes, but suffers numerical drawbacks throughout structurally advanced structures, industrial applications or for analysing localized wave propagation in broadband frequency ranges. In this paper, a formulation of the WFEM is described for 2D homogeneous or periodic waveguides, requiring large-scaled, structurally advanced or finely meshed finite element descriptions. This formulation is based on a two-scale model order reduction combining modal synthesis in the periodic cell with an innovative wave-based interpolation strategy, involving a projection of the spectral problem on a solution subspace associated with propagating wave shapes. The proposed formulation was applied for the k-space analysis of a bi-directionally stiffened panel on the frequency range [0 -2500Hz]. It provided very good results and a significant reduction of the computational cost. The stability of the method was investigated considering two key parameters of the wave -mode selection procedure: First, the influence of the fixed-interface modes used in the CMS to predict local resonances effects (important for stiffened panels) was highlighted. Noteworthy, these modes are essential to describe the inner cells displacement field, thereby the vibro-acoustic characteristics

Y (normalized)

Figure 11: Deformed shape of the periodic cell for the shear wave at 750 Hz. Blue ' ': undeformed ; red ' ': deformed shape of the panel. Therefore, the authors recommend to select the modes between 0 and 2 × ω max , where ω max is the upper limit of the frequency range considered. Secondly, the effects of wave correlation criteria ε were investigated to study the method's robustness. The proposed wave-based interpolation was successfully applied to reduce the spectral problem arising from Bloch's theorem, using propagating wave shapes participations on the cell's edges instead of the nodal displacements. A wave correlation criteria ε of 0.99 was chosen and provided accurate dispersion curves and significant reduction of the computational time. Besides, since most of the computational cost is associated to the propagating waves determination, an iterative procedure can be employed to enrich the wave basis to reach a required level of accuracy. It is hoped that it will provide an efficient numerical tool, encompassing numerical difficulties associated with the spectral analysis of industrial of structurally advanced 2D, continuous or periodic waveguides. An efficient prediction of local resonance effects and high-order propagating waves will provide useful results for the acoustic analysis of periodic panels.

T L q T B q T R q T T P T C T = q b P C (A.1)
Then, using the periodicity relations defined section 2.1, the displacement vector involving modal participations can be written:

q b P C = Λ R               q 1 q L q B P C               (A.2) with Λ R = Λ R 0 0 I C (A.3)
where I C is an identity matrix of the size of P C , and Λ R contains the propagation constants from the periodicity conditions:

Λ R =                                    I s 0 0 λ x I s 0 0 λ y I s 0 0 λ x λ y 0 0 0 I 0 0 0 I 0 λ x I 0 0 0 λ y I                                    (A.4)
Note that the vector P C is kept in the spectral eigenproblem since the frequency is unknown and no condensation of the inner DOFs can be done. Similarly, the matrix Λ L being defined as:

Λ L f b = 0 (A.5)
it yields:

Λ L = Λ L 0 0 I C (A.6)
where

Λ L = Λ H R =           I λ -1 x I λ -1 y I λ -1 x λ -1 y I 0 0 0 0 0 0 0 0 I 0 λ -1 x I 0 0 0 0 0 0 I 0 λ -1 y I           (A.7)
Therefore, the eigenproblem providing mid-aliasing frequency values: k x = π 2d x , retained for the reduced solution subspace computation, can be written, with λ y given and λ x = exp(-j π 2 ) =j as:

Λ L (λ x , λ y )B T Φ K -ω 2 j M B Φ Λ R (λ x , λ y ) Φ = 0 (A.8)
Appendix B. Formulation of the spectral eigenproblem for the 2D periodic waveguide

The condensed dynamic stiffness defined Eq.( 22) can be ordered as follows:

                                  
                                  
Using the matrices Λ L and Λ R defined in Appendix A, the product Λ L (λ x , λ y ) D(ω) Λ R (λ x , λ y ) can be developed using the decomposition proposed for D. Therefore, the coefficients A, B and C are constant in a given (ω, λ y ) formulation. Eq.( 23) becomes a quadratic eigenvalue problem in λ x , written:

          λ x           A 11 A 1L A 1B A L1 A LL A LB A B1 A BL A BB           +           B 11 B 1L B 1B B L1 B LL B LB B B1 B BL B BB           + 1 λ x           C 11 C 1L C 1B C L1 C LL C LB C B1 C BL C BB                               q 1 q L q B           = 0 (B.2)
where λ x coefficients are:

Figure 1 :

 1 Figure 1: Illustration of the 2D waveguide periodic cell.

Figure 2 :

 2 Figure 2: Phase constant surfaces of a bi-directionally stiffened panel, obtained using given-(k x , k y ) WFE formulation.

Figure 4 :

 4 Figure 4: Illustration of the stiffened panel's periodic cell

Figure 6 :

 6 Figure 6: Comparison of the k-space solutions at 1500 Hz using standard WFEM and the proposed formulation (black 'o').

1 )Figure 7 :

 17 Figure 7: Comparison of the k-space solutions at 2400 Hz using standard WFEM and the proposed formulation (black 'o').

Figure 8 :Figure 9 :

 89 Figure 8: Influence of the modal basis used in the CMS procedure. Real part of wavenumbers in the y-direction (k x = 0). using three different mode selection for the modal synthesis. The first basis (CMS 2500 Hz, 'o') involves the four modes until 2500 Hz. It show good agreement with the full model, including the frequency band-gaps for the flexural waves at 850 Hz and 2000 Hz. The wave shapes associated with both local resonances at 800 Hz and 2000 Hz

Figure 10 :

 10 Figure 10: Deformed shape of the periodic cell for the longitudinal wave at 1000 Hz. Blue ' ': undeformed ; red ' ': deformed shape

Figure 12 :Figure 13 :

 1213 Figure 12: Influence of the wave selection criteria on the reduced model accuracy. The wavenumbers k y are determined for a given k x = 0.

  q

D 11 D

 11 12 D 13 D 14 D 1L D 1B D 1R D 1T D 21 D 22 D 23 D 24 D 2L D 2B D 2R D 2T D 31 D 32 D 33 D 34 D 3L D 3B D 3R D 3T D 41 D 42 D 43 D 44 D 4L D 4B D 4R D 4T D L1 D L2 D L3 D L4 D LL D LB D LR D LT D B1 D B2 D B3 D B4 D BL D BB D BR D BT D R1 D R2 D R3 D R4 D RL D RB D RR D RT D T1 D T2 D T3 D T4 D TL D TB D TR D TT

A 11 =

 11 D 12 + D 34 + D 32 λ -1 y + D 14 λ y (B.3a)A 1L = D 1R + D 3R λ -1 y (B.3b) A L1 = D L2 + D L4 λ y (B.3c) A LL = D LR (B.3d) A B1 = D B2 + D T4 + D T2 λ -1 y + D B4 λ y (B.3e) A BL = D BR + D TR λ -1 y (B.3f) A 1B = A LB = A BB = 0 (B.3g)the constant coefficients are:B 11 = D 11 + D 22 + D 33 + D 44 + (D 31 + D 42 )λ -1 y + (D 13 + D 24 )λ y (B.4a) B 1L = D 1L + D 2R + (D 3L + D 4R )λ -1 y (B.4b) B 1B = D 1B + D 3T + D 3B λ -1 y + D 1T λ y (B.4c) B L1 = D L1 + D R2 + D L3 λ y + D R4 λ y (B.4d) B LL = D LL + D RR (B.4e) B LB = D LB + D LT λ y (B.4f) B B1 = D B1 + D T3 + D T1 λ -1 y + D B3 λ y (B.4g) B BL = D BL + D TL λ -1 y (B.4h) B BB = D BB + D TT + D TB λ -1 y + D BT λ y (B.4i)and the coefficients in 1 λ x are:C 11 = D 21 + D 43 + D 41 λ -1 y + D 23 λ y (B.5a) C 1L = D 2L + D 4L λ -1 y (B.5b) C 1B = D 2B + D 4T + D 4B λ -1 y + D 2T λ y (B.5c) C L1 = D R1 + D R3 λ y (B.5d) C LL = D RL (B.5e) C LB = D RB + D RT λ y (B.5f) C B1 = C BL = C BB = 0 (B.5g)

Table 1 :

 1 Stiffened plate: CPU time for the standard and reduced WFE formulations.

	Time per iteration (s)	WFEM CMS-WFEM Interp-WFEM Proposed method
	Modal Analysis (×1)		258		258
	Dynamic Condens.(×500)	1.85	0.082	1.85	0.082
	Building wave basis (×1)			17	8.16
	Spectral problem (×500)	1.55	1.55	0.066	0.066
	Total Time	1700	1074 (-37%)	975 (-43%)	340 (-80%)

Acknowledgements

The authors would like to gratefully acknowledge Airbus Helicopters for their financial support.

Appendix A. Determination of a reduced frequency subset using a fixed (k x , k y ) formulation

The following convention is used:


Member, IEEE Peng Xu
email: xupeng@cqupt.edu.cn

Senior Member, IEEE, Gaofeng Pan, Senior Member, IEEE Gaojie Chen
email: gaojie.chen@leicester.ac.uk

Fellow, IEEE Marco Di Renzo
email: marco.direnzo@centralesupelec.fr

Ergodic Secrecy Rate of RIS-Assisted Communication Systems in the Presence of Discrete Phase Shifts and Multiple Eavesdroppers

Keywords: Reconfigurable intelligent surface, discrete phase shifts, multiple eavesdroppers, ergodic secrecy rate

This letter investigates the ergodic secrecy rate (ESR) of a reconfigurable intelligent surface (RIS)-assisted communication system in the presence of discrete phase shifts and multiple eavesdroppers (Eves). In particular, a closed-form approximation of the ESR is derived for both non-colluding and colluding Eves. The analytical results are shown to be accurate when the number of reflecting elements of the RIS N is large. Asymptotic analysis is provided to investigate the impact of N on the ESR, and it is proved that the ESR scales with log 2 N for both non-colluding and colluding Eves. Numerical results are provided to verify the analytical results and the obtained scaling laws.

I. INTRODUCTION

Reconfigurable intelligent surfaces (RISs) utilize a large number of passive reflecting elements to customize wireless communication environments [START_REF] Di Renzo | Smart radio environments empowered by reconfigurable AI meta-surfaces: An idea whose time has come[END_REF]- [START_REF] Di Renzo | Reconfigurable intelligent surfaces vs. relaying: Differences, similarities, and performance comparison[END_REF]. Due to the low-cost, high energy-efficiency and full-duplex advantages, RISs are regarded as a promising technology for next-generation wireless communications and hence have recently received significant academic and industrial attention [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]- [START_REF] Pan | Intelligent reflecting surface aided mimo broadcasting for simultaneous wireless information and power transfer[END_REF].

RISs have various potential applications in wireless communications, which include the design of secure wireless systems based on the concept of physical layer security (e.g., [START_REF] Cui | Secure wireless communication via intelligent reflecting surface[END_REF]- [START_REF] Yang | Secrecy performance analysis of RIS-aided wireless communication systems[END_REF]). In [START_REF] Cui | Secure wireless communication via intelligent reflecting surface[END_REF]- [START_REF] Yu | Robust and secure wireless communications via intelligent reflecting surfaces[END_REF], the authors investigated optimization problems to jointly design the beamforming vectors and phase shifts at the transmitter and RIS, respectively. In general, there exist two objectives for the design of the phase shifts at the RIS: (i) to strengthen the legitimate channels by co-phasing the reflected signals with the signal directly received from the transmitter; and (ii) to suppress the eavesdropping channels by setting the reflected signals at the eavesdroppers (Eves) to be in opposite phase with respect to the signal from the transmitter. The key idea behind the optimization problems in existing works [START_REF] Cui | Secure wireless communication via intelligent reflecting surface[END_REF]- [START_REF] Yu | Robust and secure wireless communications via intelligent reflecting surfaces[END_REF] lies in achieving a favorable trade-off between these two design objectives, which requires the knowledge of the instantaneous eavesdropping channel state information (CSI) at the transmitter and RIS. However, the instantaneous eavesdropping CSI is difficult to obtain in practice, since the Eves are usually passive and do not actively communicate with other nodes. Motivated by this consideration, the authors of [START_REF] Wang | Intelligent reflecting surfaces assisted secure transmission without eavesdropper's CSI[END_REF] and [START_REF] Yang | Secrecy performance analysis of RIS-aided wireless communication systems[END_REF] considered RIS-assisted secrecy communications without assuming the knowledge of the instantaneous eavesdropping CSI.

Different from these existing works, this letter investigates the ergodic secrecy rate (ESR) of RIS-assisted systems in the presence of discrete phase shifts and multiple Eves. In particular, by approximately characterizing the distribution of the received signal-to-noise-ratios (SNRs) at the Eves, we obtain a closed-form approximation of the ESR for both noncolluding and colluding Eves. The analysis of the ESR, in fact, is essentially different from the analysis of the ergodic rate without security constraints [START_REF] Li | Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors[END_REF], since the phase shifts at the RIS lead to a different impact on the intended receiver and Eves. Moreover, based on passive beamforming, the received SNRs at the destination and Eves depend on phase quantization errors and cascaded channels, that are different from those in massive multiple-input multiple-output (MIMO) systems. In order to provide insights, asymptotic analysis is also provided, which shows that the ESR scales with log N for both non-colluding and colluding Eves. Numerical results are illustrated to verify that the analytical results are accurate for large values of N .

Notation: C and Z denote the complex domain and integer set, respectively; we denote [1 : M ] {1, . . . , M }, where M is a positive integer, and [x] + max{0, x}; CN denotes the complex Gaussian distributions; E[•] denotes the expectation of a random variable; log(•) and ln(•) denote the base-two and natural logarithms, respectively; and κ is Euler's constant.

II. SYSTEM MODEL AND PRELIMINARIES

We consider an RIS-assisted secure communication system with a source (S), an RIS (R) with N reconfigurable elements, a destination (D) and K Eves (E k , ∀k ∈ [1 : K]). The reconfigurable elements of the RIS are arranged in a uniform array of tiny antennas spaced half of the wavelength apart. All nodes are assumed to be equipped with a single antenna 1

. The channels S → D, S → E k , S → R, R → D and R → E k are denoted by h SD ∈ C, h SE k ∈ C, h SR ∈ C N ×1 , h RD ∈ C N ×1 and h k ∈ C N ×1
, respectively. These channels are modeled as

h SD = g SD d -α 2 SD , h SE k = g SE k d -α 2 SE k , [h SR ] n = g SR,n d -α 2 SR , [h RD ] n = g RD,n d -α 2 RD and [h k ] n =
1 RIS-aided transmission has several applications when multiple antennas are not available at either the transmitter or the receiver, e.g., in device-todevice communications [START_REF] Di Renzo | Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead[END_REF].

g k,n d -α 2 k
, where g SD , g SE k , g SR,n , g RD,n , g k,n ∼ CN (0, 1) denote the small-scale fading 2

, d SD , d SE k , d SR , d RD and d k denote the distances S → D, S → E k , S → R, R → D and R → E k , respectively, ∀n ∈ [1 : N ], k ∈ [1 : K],
and α is the pass-loss exponent. Then, the received signal at D and E k can be written as

yD = √ P ηh T SR ΦhRD + hSD xS + nD, (1)
yE k = √ P ηh T SR Φh k + hSE k xS + nE k , (2)
respectively, where x S is the transmitted signal, E(|x S | 2 ) = 1, P is the transmit power, n D and n E k ∼ CN (0, δ 2 ) are the additive white Gaussian noises at D and E k , respectively, η ∈ (0, 1] is the amplitude reflection coefficient, Φ diag(e jφ1 , . . . , e jφN ) and φ n ∈ [0, 2π) is the phase shift of the nth element of the RIS. We assume that the RIS does not have access to the instantaneous eavesdropping CSI, so that it cannot design φ n in order to suppress the received SNRs at the Eves. However, the RIS is assumed to know the instantaneous legitimate CSI. Under these assumptions, the optimal value of φ n that maximizes the received SNR at D is φ * n = θSD -θSR,n -θRD,n, where θ SD , θ SR,n and θ RD,n denote the phases of g SD , g SR,n and g RD,n , respectively. Due to hardware limitations, φ n can only take a number of discrete values. In particular, the set of discrete phase shifts is denoted by

F 0, 2π 2 b , . . . , (2 b -1)2π 2 b
, where b denotes the number of quantization bits. Accordingly, we set

φ n = f 1 (φ * n ), where the function f 1 (φ * n ) maps φ * n to the nearest point in F , i.e., f1(φ * n ) = φi, if |φ * n -φi| ≤ |φ * n -φj |, φi, φj ∈ F, ∀j = i. (3)
Therefore, the phase quantization error is

Θ n = f 1 (φ * n ) -φ * n , which is uniformly distributed in -π 2 b , π 2
b , similar to [START_REF] Li | Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors[END_REF], [START_REF] Xu | Reconfigurable intelligent surfaces assisted communications with discrete phase shifts: How many quantization levels are required to achieve full diversity?[END_REF]. According to (1) and ( 2), the received SNRs at D and E k can be formulated, respectively, as follows

γD = ρ |hSD| + η N n=1 |[hSR]n[hRD]n| e jΘn 2 = ρ d -α 2 SD |gSD| + ηd -α 2 SR d -α 2 RD N n=1 |gSR,ngRD,n|e jΘn 2 , (4)
γE k = ρ hSE k + η N n=1 |[hSR]n[h k ]n| e jψ k,n 2 = ρ d -α 2 SE k gSE k + ηd -α 2 SR d -α 2 k N n=1 |gSR,ng k,n | e jψ k,n 2 , (5)
where

ψ k,n f 2 (φ * n , θ SR,n ) + θ k,n , θ k,n is the phase of g k,n and the function f 2 (φ * n , θ SR,n ) is defined as follows f2(φ * n , θSR,n) f1(φ * n ) + θSR,n. (6)
The ESR 3 can be expressed as follows

Rs = [RD -RE] + , (7)
where

R D = E γD [log(1 + γ D )]
and R E denote the ergodic rates from S to D and from S to the Eves, respectively. Given {Θ n } N n=1 , an approximated expression of R D can be found in [15, Eq. ( 13)]. By averaging over {Θ} N n=1 , R D can be calculated as shown in [START_REF] Pan | Intelligent reflecting surface aided mimo broadcasting for simultaneous wireless information and power transfer[END_REF] at the top of the next page, where

A 1 η 2 d -α SR d -α RD , A 2 √ πη2 b 4 d -α 2 SD d -α 2 SR d -α 2 RD sin π 2 b and A 3 η 2 2 2b 32 d -α SR d -α RD 1-cos 2π 2 b . Remark 1:
The analysis of the ESR for the considered RISassisted system relies only on the knowledge of the statistical eavesdropping CSI, which can be obtained by using several methods, such as those used in [START_REF] Shang | Unmanned aerial vehicle meets vehicle-to-everything in secure communications[END_REF].

In the following sections, R E is calculated for non-colluding and colluding Eves, respectively.

III. NON-COLLUDING EVES

In the non-colluding case, R E can be expressed as follows

RE = max k∈[1:K] RE k , (9)
where

R E k E γE k [log(1 + γ E k )].
In order to derive R E k , the distribution of γ E k in (5) needs to be computed.

A. Distribution of γ E k

Before deriving the distribution of γ E k , we introduce the following lemma.

Lemma

1: The phase ψ k,n , k ∈ [1 : K], n ∈ [1 : N ], in (5) has the following properties: a) ψ k,n is uniformly distributed in [0, 2π); b) ψ k,n is independent of f 2 (φ * n , θ SR,n ) defined in (6); c) ψ k,i is independent of ψ k,j , ∀i = j, i, j ∈ [1 : N ]. Proof 1: See Appendix A.
Based on Lemma 1, the distribution of γ E k in ( 5) is provided in the following lemma.

Lemma 2: When N is large, γ E k can be approximated with an exponential random variable with mean Lemma 2] and the fact that {ψ k,n } N n=1 are independent and identically distributed uniform random variables in [0, 2π) as proved in Lemma 1, [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF], the proof follows.

λ E k = ρ d -α SE k + N B k , where B k η 2 d -α SR d -α k . Proof 2: Define G k N n=1 |g SR,n g k,n | e jψ k,n , k ∈ [1 : K]. Based on [20,
G k ∼ CN (0, N ) as N → ∞. Furthermore, since g SE k is independent of G k , we have d -α 2 SE k g SE k +ηd -α 2 SR d -α 2 k G k ∼ CN 0, d -α SE k +N B k , as N → ∞. Recalling that γ E k = ρ d -α 2 SE k g SE k + √ B k G k 2 in (
Remark 2: The authors of [START_REF] Li | Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors[END_REF] approximated R D based on the fact that Jensen's inequality is tight, rather than an upper bound, if Var[γD] E 2 [γD] → 0, as N → ∞. However, Lemma 2 shows that Jensen's inequality cannot be applied for approximating R E , since 3 We assume that the RIS appropriately customizes the wireless channel but we consider that the distribution of the signal transmitted by S is always Gaussian. It is worth mentioning that the information-theoretic characterization of RIS-assisted transmission and the calculation of the optimal input distribution in the presence of an RIS is an open issue that is currently under active research [START_REF] Karasik | Beyond max-snr: Joint encoding for reconfigurable intelligent surfaces[END_REF]. This is, however, beyond the scope of this letter.

Var[γE k ] E 2 [γE k ] → 1, as N → ∞.
RD ≈ log 1+ρ N η 2 d -α SR d -α RD +d -α SD + π 3 2 η 4 d -α 2 SD d -α 2 SR d -α 2 RD N n=1 EΘ n [cos Θn]+ π 2 η 2 8 d -α SR d -α RD N-1 i=1 N k=i+1 EΘ i ,Θ k [cos(Θ k -Θi)] = log 1 + ρN A1 + ρd -α SD + ρN A2 + ρN (N -1)A3 , (8)

B. Ergodic Secrecy Rate

The ESR for non-colluding Eves is summarized in the following theorem.

Theorem 1: When N is large, the ESR for non-colluding Eves can be expressed as follows

Rs ≈ RD + 1 ln 2 max k∈[1:K] e 1 λ E k Ei - 1 λE k + , (10)
where Ei(x) -

∞ -x e -t
t dt, x < 0, is the exponential integral function [START_REF] Gradshteyn | Table of Integrals, Series and Products[END_REF]Eq. 8.211].

Proof 3: When N is large, based on Lemma 2, R E k in (9) can be approximated as follows

R E k ≈ ∞ 0 log(1 + x) 1 λ E k e -x λ E k dx = 1 ln 2 ∞ 0 e -x λ E k 1 + x dx = - e 1 λ E k ln 2 Ei - 1 λ E k , (11)
where the last equality is based on [START_REF] Gradshteyn | Table of Integrals, Series and Products[END_REF]Eq. 3.352.4]. Combining ( 7), ( 9) and ( 11), the theorem is proved.

C. Asymptotic Analysis

To obtain insights from the obtained ESR, its asymptotic behavior is analyzed in the following corollary.

R D ≈ log ρA 3 N 2 1 ρA 3 N 2 + A 1 A 3 N + d -α SD A 3 N 2 + A 2 -A 3 A 3 N +1 → 2 log N + log ρ + log A 3 , as N → ∞. (12)
In addition, R E k in ( 11) can be further expressed as follows

R E k (a) ≈ e 1 λ E k ln 2 -κ + ln(λ E k ) + ∞ i=1 (-1) i+1 i • i! • λ i E k (b) → log(λ E k ) - κ ln 2 (c) = log (N ρB k ) + log 1 + d -α SE k N B k - κ ln 2 → log N + log ρ + log B k - κ ln 2 , as N → ∞. ( 13
)
where (a) is based on [START_REF] Guan | Intelligent reflecting surface assisted secrecy communication: Is artificial noise helpful or not?[END_REF] 7), ( 9), ( 12) and ( 13), the corollary follows.

Remark 3: Compared with the scaling law 2 log N for nonsecrecy transmission with discrete phase shifts [START_REF] Li | Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors[END_REF], Corollary 1 shows that the ESR scales with log N .

IV. COLLUDING EVES

When the Eves are colluding, they can combine their received signals for information interception. Based on [START_REF] Pinto | Secure communication in stochastic wireless networks-part II: Maximum rate and collusion[END_REF], R E in (7) can be expressed as follows

RE = E {γ E k } K k=1 log 1 + K k=1 γE k . ( 14
)
Since common random variables {g SR,n } N n=1 are present in every γ E k as shown in [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF], {γ E k } K k=1 are correlated random variables. However, the following lemma shows that such correlation is negligible for large values of N .

Lemma

3: γ Ei is independent of γ Ej if N → ∞, i, j ∈ [1 : K], i = j. Proof 5: Let us define HE k d -α 2 SE k gSE k + √ B k N n=1 |gSR,ng k,n | e jψ k,n . Thus, γ E k = ρ|H E k | 2 . According to Lemma 1-a), E e jψ k,n = 0 and E [HE k ] = 0, ∀k ∈ [1 : K]. Moreover, since ψ i,n is independent of ψ j,m if i = j or n = m, we have E[HE i HE j ] = 0, ∀i = j.
Therefore, the covariance of H Ei and H Ej is zero. Furthermore, based on Lemma 2, {H k } K k=1 are uncorrelated complex Gaussian variables if N → ∞, and hence {H k } K k=1 are independent of each other. This completes the proof.

A. Ergodic Secrecy Rate

Based on Lemma 2 and Lemma 3, the ESR is provided in the following theorem.

Theorem 2: When N is large, and λ Ei = λ Ej , ∀i = j, i, j ∈ [1 : N ], the ESR for colluding Eves can be approximated as follows

Rs ≈   RD + 1 ln 2 K i=1 e 1 λ E i Ei - 1 λE i K j=1,j =i λE i λE i -λE j   + . ( 15
)
Proof 6: Based on Lemma 2 and Lemma 3, if λ Ei = λ Ej , ∀i = j, K k=1 γ E k has the following probability density function (PDF) [START_REF] Ross | Introduction to probability models[END_REF]:

f K k=1 γ E k (x) = K i=1 1 λE i e -x λ E i K j=1,j =i λE i λE i -λE j . ( 16
)
Combining ( 11), ( 14) and ( 16), we obtain

RE ≈ - 1 ln 2 K i=1 e 1 λ E i Ei - 1 λE i K j=1,j =i λE i λE i -λE j . ( 17
)
Recalling ( 7), the theorem is proved. Remark 4: Theorem 2 corresponds to the case that the Eves lie in different locations, so that {γ E k } K k=1 have different means. When the Eves are clustered relatively closely together, {γ E k } K k=1 have the same (or a very similar) mean. In this case, the ESR can be analyzed in a similar way, whose details are not provided due to space limitations.

B. Asymptotic Analysis

The asymptotic behavior of the obtained ESR for colluding Eves is provided in the following corollary.

Corollary 2:

As N → ∞, R s → log N + log A 3 + κ ln 2 - K i=1 log B i K j=1,j =i

Bi

Bi-Bj , which implies that the ESR for colluding Eves also scales with log N .

Proof 7: From ( 11), ( 13) and ( 17), we have

R E → K i=1 log N +log ρ- κ ln 2 +log B i K j=1,j =i λ Ei λ Ei -λ Ej (a) → K i=1 log N +log ρ- κ ln 2 +log B i K j=1,j =i B i B i -B j (b) = log N +log ρ- κ ln 2 + K i=1 log B i K j=1,j =i B i B i -B j , (18)
where (a) holds since λ

E k = ρN B k 1 + d -α SE k N B k → ρN B k
as N → ∞, and (b) is based on the fact that

K i=1 K j=1,j =i

Bi

Bi-Bj = 1, as proved in [START_REF] Ross | Introduction to probability models[END_REF]Chapter 5]. Combining [START_REF] Pan | Multicell mimo communications relying on intelligent reflecting surfaces[END_REF], ( 12), ( 14) and ( 18), the proof follows. Remark 5: Comparing Corollaries 1 and 2, we evince that only the last terms for the asymptotic ESR are different, i.e., max k∈[1:K] log B k and

K i=1 log B i K j=1,j =i

Bi

Bi-Bj for noncolluding and colluding Eves, respectively. In addition, the ESRs for both non-colluding and colluding Eves have the same scaling law, i.e., log N .

C. Large Number of Eves

To obtain more insights from the obtained ESR, we provide a simplified expression of R s for large values of K in the following corollary.

Corollary 3: When N and K → ∞, the ESR for colluding Eves can be approximated as follows

Rs ≈ RD -log 1 + K k=1 λE k + . ( 19
)
Proof 8: According to [START_REF] Sanayei | Opportunistic beamforming with limited feedback[END_REF]Theorem 4], R E in ( 14) can be approximated as RE ≈ log 1 + E K k=1 γ k . This completes the proof.

Based on Lemma 3, Var K k=1 γ E k → K k=1 λ 2 E k as N → ∞. Thus, Var[ K k=1 γ k ] (E[ K k=1 γ k ]) 2 → K k=1 λ 2 E k ( K k=1 λE k ) 2 → 0 as N and K → ∞.
V. NUMERICAL RESULTS In this section, numerical results are provided to verify the analytical results stated in the theorems and corollaries. For Fig. 1 shows the impact of the number of reconfigurable elements N on the ESR, when the number of Eves is K = 5. We observe that the approximated analytical results in Theorems 1 and 2 match well with Monte Carlo simulations almost for all values of N . This is because the ESR is zero for small values of N , and starts to be positive for large values of N , i.e., N ≥ 46 in the considered case. We also observe that the ESR increases with N . For example, the ESRs are about 3.7 bps/Hz and 2.5 bps/Hz for non-colluding and colluding Eves, respectively, if N = 100. In addition, the analytical results obtained in Corollaries 1 and 2 asymptotically approach the simulations as N becomes sufficiently large, which confirms the scaling laws. The setup with non-colluding Eves provides a larger secrecy rate since only the "best" Eve determines the ESR. There exists a constant gap of about 1 bps/Hz between the ESRs for non-colluding and colluding Eves if N exceeds 10 4 .

In Fig. 2, the ESR as a function of the number of Eves K for N = 250 is shown. The figure confirms the findings in Corollary 3 for colluding Eves, and we observe that the approximation in [START_REF] Shang | Unmanned aerial vehicle meets vehicle-to-everything in secure communications[END_REF] becomes tighter as K increases. The ESR for non-colluding Eves is less affected by K. We observe, in particular, that there exists an ESR floor of about 5 bps/Hz for large values of K. This is because the ESR for noncolluding Eves is determined by the nearest Eve to the source.

In the considered simulation setup, the nearest Eve is located at around (0, -20) m, when K is large.

VI. CONCLUSIONS

This letter investigated the ESR of an RIS-assisted communication system in the presence of discrete phase shifts and multiple Eves. We obtained an approximated closed-form expression of the ESR and unveiled that the ESR scales with log N in the presence of both non-colluding and colluding Eves. An interesting future direction is to analyze the ESR of RIS-assisted systems in the presence of multi-antenna transmitters. For example, the recent research works in [START_REF] Qian | Beamforming through reconfigurable intelligent surfaces in single-user MIMO systems: SNR distribution and scaling laws in the presence of channel fading and phase noise[END_REF], [START_REF] Zappone | Overhead-aware design of reconfigurable intelligent surfaces in smart radio environments[END_REF] could be generalized in order to take into account security constraints.

APPENDIX A PROOF OF LEMMA 1

Let us denote X 1,n f 2 (φ * n , θ SR,n ), X 2,n θ k,n and Y n ψ k,n . Accordingly, Y n = X 1,n + X 2,n as shown in Section II. We note that the random phases X 1,n , X 2,n and Y n have a uniform circular distribution.

Given X 1,n = x 1 , ∀x 1 ∈ [0, 2π),

Y n = x 1 + X 2,n is uni- formly distributed in [x 1 , x 1 + 2π) = [x 1 , 2π) ∪ [2π, x 1 + 2π).
Since Y n has a circular uniform distribution, [2π, x 1 + 2π) is equivalent to [0, x 1 ). Thus, Y n is uniformly distributed in [0, 2π), which proves Lemma 1-a).

Since X 1,n and X 2,n are independent, their joint PDF is We can construct the following Jacobian matrix

Thus, the joint PDF of X 1,n and Y n can be written as

which implies that Y n is independent of X 1,n . Thus, Lemma 1-b) is proved. For ∀i = j and i, j ∈ [1 : N ], we have Y i = X 1,i + X 2,i and Y j = X 1,j + X 2,j . Although the same random phase θ SD is present in both Y i and Y j as shown in Section II, Y i is still independent of Y j , due to the following two facts: (i) Y i and Y j are independent of X 1,i and X 1,j , respectively, according to Lemma 1-b); (ii) X 2,i is independent of X 2,j . Therefore, Lemma 1-c) is proved.

Corollary 1 :ln 2 -

 12 As N → ∞, R s → log N + log A 3 + κ max k∈[1:K] log B k , which implies that the ESR for noncolluding Eves scales with log N .Proof 4: From (8), we have

Fig. 1 .

 1 Fig. 1. Ergodic secrecy rate vs. N , for K = 5.

  illustrative purposes, we set α = 3, b = 3 bits, P = 20 dBm, σ 2 = -96 dBm and η = 0.8. In addition, S, R and D are located at (0, 0) m, (100, 0) m and (90, 20) m, respectively. As for the Eves, E k is located at 90k K , -20 m, k ∈ [1 : K]. For the considered simulation setup, the ESR is zero in the absence of RIS. In this case, in fact, the ESR for non-colluding Eves can be expressed as Rs = f3(dSD) -max k∈[1:K] f3(dSE k ) -α /(1 + t)dt, x > 0. Therefore, R s = 0 since f 3 (x) is a decreasing function of x and d SD ≥ d SE k in the considered network configuration, ∀k ∈ [1 : K].

fX 1 ,

 1 n ,X 2,n (x1, x2) = fX 1,n (x1)fX 2,n (x2) = 1 2π fX 1,n (x1).(20)

JX 1 , 1 ∂x 1

 111 n ,Yn (x1, x2) = ∂x

fX 1 ,

 1 n ,Yn (x1, y) = fX 1,n ,X 2,n (x1, x2) det(JX 1,n ,Yn (x1, x2)) = 1 2π fX 1,n(x1)= fX 1,n (x1)fY n (y),

  Fig. 2. Ergodic secrecy rate vs. K, for N = 250.

		8
						Non-colluding Eves, simulation result
						Non-colluding Eves, approx. result (Theorem 1)
		7				Colluding Eves, simulation result
						Colluding Eves, approx. result (Theorem 2)
						Colluding Eves, approx. result (Corollary 3)
	Ergodic secrecy rate (bps/Hz)	3 4 5 6
		2
		1
		2	4	6	8	10 12 14 16 18 20 22 24 26 28 30
						Number of Eves K

Since the elements of the RIS are spaced half of the wavelength apart and we assume that the location of the RIS cannot be optimized to ensure strong line-of-sight links, we have, as a first approximation similar to[START_REF] Yang | Secrecy performance analysis of RIS-aided wireless communication systems[END_REF]-[START_REF] Qian | Beamforming through reconfigurable intelligent surfaces in single-user MIMO systems: SNR distribution and scaling laws in the presence of channel fading and phase noise[END_REF], that the channels can be modeled as independent and identically distributed, and follow a Rayleigh distribution.


C Droz

Z Zergoune

R Boukadia

O Bareille

M N Ichchou

Vibro-acoustic optimisation of sandwich panels using the wave/finite element method

Keywords: Honeycomb, cellular core, finite element, sandwich, wave, acoustic, design

This paper investigates the use of a wave-based method in the framework of structural optimisation of composite panels involving advanced components. The wave/finite element method (WFEM) is used to evaluate the influence of a core's geometry on the transition frequency of a sandwich panel involving composite skins. This transition occurs is a sandwich panel when the transverse shear stiffness has a significant influence on the flexural motion, compared to the bending stiffness. It follows that the modal density and the acoustic radiation will considerably increase above this frequency. The periodic waveguide is modelled at the mesoscopic scale using a 3D finite element model of the unit-cell. Therefore this method does not require an homogenisation of the core based on Gibson and Ashby formulations to provide the wave dispersion characteristics. Although the cellular cores compared in this study share the same mass-to-stiffness ratio, a significant alteration of the transition frequency and modal density can be observed compared to honeycomb cores. A periodic octagonal core is designed, providing up to 70 % increase of the transition frequency and a significant reduction of the modal density.

Introduction

Sandwich panels are extensively used in automotive, naval and aerospace industry. These structures have high stiffness-to-weight ratios and are usually made of a moderately thick, lightweight core surrounded by glass or carbonfibre composite skins. The core can be made of polymer foam or more complex structures involving aluminium or resin honeycomb constructions. Although the material and geometrical properties of the panel are usually designed to provide specific stiffness and density characteristics, structural optimisation also concerns acoustic radiation efficiency in order to improve the acoustic comfort. Therefore, there is an increasing need for reliable optimisation tools for design engineers, providing fast vibro-acoustic evaluation of large-scaled structures involving composite components.

In this context, the knowledge of accurate wave dispersion characteristics in two laminated orthotropic skins connected by 3D cellular honeycomb core is a key information for the prediction of the acoustic transmission parameters.

Numerical methods for analysing the vibrational behaviour of complex composite or periodic panels in a broadband frequency range were extensively investigated in the last decade. Honeycomb sandwich panels were usually modelled using classical laminate plate (CPLT) of Reissner-Mindlin (FSDT) theories involving three homogeneous layers. An extensive number of analytical formula for analysing wave dispersion characteristics in two-dimensional waveguides can be found in the literature (see Reddy [START_REF] Reddy | Mechanics of laminated composite plates -Theory and analysis[END_REF]), starting with the asymptotic model for analysing wave dispersion in thick symmetric sandwich panels developed by Kurtze and Watters [START_REF] Kurtze | New wall design for high transmission loss or high damping[END_REF], then later in the classic work of Wilkinson [START_REF] Wilkinson | Modal densities of certain shallow structural elements[END_REF] and Erickson [START_REF] Erickson | Modal densities of sandwich panels: theory and experiment[END_REF]. The prediction of wave propagation characteristics of structurally advanced structures using analytical models has been a subject of intense research over the past years. Renji et al. [START_REF] Renji | Modal density of composite honeycomb sandwich panels[END_REF] introduced a new analytical formulation providing a core transverse shear stiffness and bending matrices in both directions for orthotropic honeycomb sandwich panels. Ghinet and Atalla [START_REF] Ghinet | Modeling thick composite laminate and sandwich structures with linear viscoelastic damping[END_REF] also used an analytical model based on Discrete Laminate Theory (DLT) to compute the dispersion characteristics of layered structures and obtain their acoustic parameters.

Recently, Guillaumie [START_REF] Guillaumie | Vibroacoustic flexural properties of symmetric honeycomb sandwich panels with composite faces[END_REF] proposed an analytical solution for the eigenmodes and modal densities of symmetric sandwich panels involving composite skins and honeycomb core. This model was found to be effective in the lowand medium-frequency range for moderately thick sandwich structures and was confirmed by experiments. Florence and Renji [START_REF] Florence | Modal density of thin composite cylindrical shells[END_REF] also derived explicit solutions for the modal density of composite cylinders.

However, these analytical solutions require a homogenisation of the hexagonal cells of the core based on Gibson and Ashby formulations [START_REF] Gibson | Cellular solids: structure and properties[END_REF]. These analytical approaches therefore suffer some drawbacks in the medium-frequency range, since they cannot handle arbitrarily shaped 3D cells for the core. Besides, these models are often inaccurate when compared to experimental characterizations in higher frequency range, especially when dealing with complex industrial products involving thick, multi-layered constructions made of anisotropic glass or carbon fibres skins.

Therefore, numerical methods were recently developed to perform wave analysis in complex periodic waveguides defined using classical finite element packages (FEM). Among others, the wave finite element method (WFEM) [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF][START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF] has been successfully applied to predict the radiation efficiency of numerous thick layered structures. Nevertheless, these models involve homogeneous orthotropic layers for the honeycomb construction, and may become inaccurate when the wavelengths are reaching the periodic cells dimensions. Although the design optimisation of honeycomb cells has received a lot of attention recently [START_REF] Alkhader | Low frequency acoustic characteristics of periodic honeycomb cellular cores: The effect of relative density and strain fields[END_REF][START_REF] Catapano | A multi-scale approach for the optimum design of sandwich plates with honeycomb core. part ii: the optimisation strategy[END_REF], the dispersion characteristics of sandwich panels involving realistic composite skins and arbitrarily shaped cores was not investigated to the author's knowledge. This paper focuses on the influence of the meso-scale parameters of the core's construction on the vibro-acoustic response of the sandwich panel. A special interest is given to the modal density, a typical vibro-acoustic parameter and the transition frequency, associated with the passage from a behaviour where the flexural wave is ruled by the skins stiffness to a behaviour where the core's transverse shear parameters governs the flexural motion. First, finite element models are developed for the periodic cells. The WFEM, which combines FEM and the Periodic Structures Theory (PST) is applied to evaluate the wave dispersion characteristics. Since numerous degrees of freedom are usually involved to ensure accuracy in a broadband frequency range, reduced order modelling (ROM) was recently investigated for periodic [START_REF] Casadei | Multiscale finite element analysis of wave propagation in periodic solids[END_REF][START_REF] Zhou | Multi-scale modelling for two-dimensional periodic structures using a combined mode/wave based approach[END_REF] and large-scaled [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF] waveguides. A hybrid wave-mode ROM is therefore applied to provide suitable computational efficiency. When dealing wit hexagonal cores, the equivalent orthotropic material can be described using Gibson and Ashby formulations, and the analytical expression of the wavenumber can be derived for the symmetric sandwich. The possibility to use a regression on the wavenumber function is therefore examined for various core's configurations, in order to derive an analytical expression of the transition frequency from the equivalent bending and shear parameters. Finally, the effects of several geometrical parameters of the core on the modal density and the transition frequency are investigated.

Wave propagation in 2D periodic structures

The wave finite element method

Numerous formulations of the WFEM were given in the literature. The method is briefly reviewed hereby. Consider a periodic cell of a 2D waveguide of dimensions d x and d y in x and y directions (see Fig. 1), modelled using a FE software. The mass M, stiffness K and damping C matrices can be extracted and the dynamic stiffness matrix is written at the circular frequency ω: D = K + jωC -ω 2 M. The degrees of freedom q (DOF) on the edges, corner and centre of the cell can be reordered as follows: 3

q = q T I q T B q T T q T L q T R q T 1 q T 2 q T 3 q T 4 T . (1)
It should be noted that apart from the mesh compatibility between the opposite edges of the unit-cell, the discretization is arbitrary inside the periodic element. Then using the periodicity relations and Bloch's theorem [START_REF] Bloch | Über die quantenmechanik der elektronen in kristallgittern[END_REF], a reduced state vector can be defined using the propagation constants λ x and λ y in x-and y-directions:

                                                                q I q B q T q L q R q 1 q 2 q 3 q 4            I 0 0 0 0 I 0 0 0 Iλ y 0 0 0 0 I 0 0 0 Iλ x 0 0 0 0 I 0 0 0 Iλ x 0 0 0 Iλ y 0 0 0 Iλ x λ y                          q I q B q L q 1                          = Rq . (2)
Then, denoting k x , k y the wavenumbers in x and y directions, the propagation constants λ x = e -jk x d x and λ y = e -jk y d y are the solutions of the problem:

R T (λ x , λ y )DR(λ x , λ y )q = 0. (3)
These solutions describe the propagating and evanescent waves in the periodic waveguide. Each solution (k x , k y , ω)

is associated with an eigenvector Ψ describing the deformed shape of a given wave. The propagating and evanescent waves can be discriminated from the imaginary parts, since ( (k x ), (k y )) describes wave attenuation in x-and ydirection, respectively. In dissipative waveguides, all the waves are decaying and a customised selection procedure, usually based on the value of | (k)| (k) can be employed.

Resolution and post-processing

It should be noted that Eq.( 3) is a transcendental eigenvalue problem whose resolution is not straightforward when internal DOFs are involved. In order to turn this problem into a classical eigenvalue problem, two parameters are usually given (the circular frequency ω and the wavenumber in one direction: k x or k y ). The directional wavenumber k(θ, ω) can be obtained from the discrete (k x (ω), k y (ω)) solutions. An accurate description of the 2D dispersion curves (also referred as k-space) hence requires an important number of resolutions of Eq.( 3).

In order to enable fast design optimisation, through broadband wave analyses of finely meshed periodic honeycomb constructions, a reduced formulation of the WFEM (see [START_REF] Droz | A hybrid wave-mode formulation for the vibro-acoustic analysis of 2d periodic structures[END_REF]) is used in the numerical applications. It employs a projection of the state vector q on solution subset involving a reduced number of wave solutions and component modes of the periodic cell.

The flexural waves can then be discriminated using a wave-matching algorithm. At each frequency step, the wave shape is compared to the whole set of propagating waves using a MAC criterion. It expresses the correlation between two vectors Φ i and Φ j by: MAC =

(Φ t i Φ j )(Φ t j Φ i ) (Φ t i Φ i )(Φ t j Φ j ) . (4)
Therefore, the wavetypes can be sorted and the flexural wave is identified by the transverse vertical motion of the unit-cell. Although this procedure is not always necessary in lower frequencies it becomes essential for the identification of the flexural wave in higher frequencies, where numerous high-order (or structural) waves can be subjected to veering effects and conversions.

Generally, refined angular sampling for k(θ, ω) is necessary for orthotropic or anisotropic structures, which may exhibit high directivity. Therefore, the modal density can be written using classical Statistical Energy Analysis (SEA) expression [START_REF] Lesueur | Rayonnement acoustique des structures: vibroacoustique, interactions fluide-structure[END_REF]:

n(ω) = π 0 n(θ, ω)dθ, (5)
where n(θ, ω) denotes the angular distribution of the modal density:

n(θ, ω) = A 2π 2 k(θ, ω) |c g (θ, ω)| . (6)
A is the surface area and c g (θ, ω) = ∂ω ∂k is the group velocity in the direction θ. An equivalent wavenumber k f (ω) is defined so that the modal density can be written:

n(ω) = A 2π k f (ω) ∂ω ∂ k f (ω) . (7)
This expression enables the comparison of scalar solutions of the dispersion relation while ensuring the accurate computation of the modal density in the 2D structure, since the directivity is taken into account. The averaged wavenumber can be expressed in terms of the angular wavenumber:

k f (ω) = 1 2π 2π 0 k 2 f (θ, ω)dθ. (8)
The modal density can be also determined by the area integration of the wavenumber space (see Langley [START_REF] Langley | The modal density and mode count of thin cylinders and curved panels[END_REF]). In case of an elliptic orthotropy along x-and y-directions, the equivalent wavenumber can be defined so that the k-space areas remain the same:

πk x k y = π k f 2 . It yields k f = k x k y .
For a sake of clarity, the equivalent flexural wavenumber < k f (ω) > and group velocity ∂ω

∂ k f (ω) are simply
denoted k and c g in the following.

Evaluation of the transition frequency

In sandwich structures, one can discriminate between three frequency domains associated with different kinetic energy distributions. In low frequency, the structure has a classical isotropic plate behaviour. The asymptotic wavenumber modulus k( f → 0) and modal density n( f → 0) are written:

k( f → 0) = m s D 1/4 √ ω ; n( f → 0) = A 2 m s D (9)
where A, D, m s denote the area, equivalent bending rigidity and mass per unit area of the plate, respectively.

In a medium frequency (usually acoustic) range, the transverse vertical shear rigidity S , increases and influences the dynamical behaviour. The asymptotic solutions for the orthotropic plate are:

k( f ) M.F. ω m s S ; n( f ) M.F. A m s S ω (10)
In terms of waves, higher frequencies are associated with fully localised propagation, meaning that flexural waves can propagate independently in both skins. The transition frequency considered here refers specifically to the passage from the isotropic plate behaviour to a motion governed by the out-of-plane shear rigidity. In Guillaumie [START_REF] Guillaumie | Vibroacoustic flexural properties of symmetric honeycomb sandwich panels with composite faces[END_REF], this frequency is calculated as the intersection between the asymptotic wavenumbers: k( f → 0) and k( f ) M.F. . The influence of a core's shear rigidity on the modal density and transition frequency of a sandwich plate was well analysed by Renji et al. [START_REF] Renji | Modal density of composite honeycomb sandwich panels[END_REF]. Above the transition frequency the modal density increases linearly with frequency and produces higher levels of acoustic radiation. In the following, the possibility to shift the transition frequency by using different cores geometries is investigated. The influence of a core's geometry on the transition frequency and modal density is the concern of this paper.

Wave dispersion characteristics in arbitrary shaped sandwich panel

Consequently, the mechanical properties and thicknesses of the skins, core and Nomex paper are identical in all the models. The sandwich panels developed below only differ in the geometry of their cores and the dimensions d x , d y of their unit-cells.

The cores are created from the parametrized models shown in Figure 2. The rectangular boxes describe the periodic cell's dimensions, while the inner lines describe the top-view position of the Nomex paper in the core. A parametrized hexagonal core is proposed in Figure 2.a. The core construction is designed using the three points M, P 1 , In Figure 2b, two additional points P 1 and P 2 are used for the construction of octagonal cores. Although this model can be easily generalised to create more complex shapes, it is reminded that the computational cost associated with optimisation procedures depends on the number of design parameters.

P
The finite element model of five unit-cells obtained using the first parametrized model is shown in Figure 3.

These models involve linear block elements having 3 degrees of freedom per node. Case (a) is a classical hexagonal model (honeycomb). Case (b) is a rectangular, Case (c) is triangular, Case (d) is a pentagonal core and Case (e) is a

parallelogram. An finely meshed octagonal core was also designed in Figure 4. It should be noted that all the models are converged in the considered frequency range, but refined meshes such as in Figure 4 would be necessary in higher frequencies to describe the local dynamics at the cell's scale.

( Besides, the geometry of the core has a noticeable influence on the Nomex distribution, hence on the surface density of the sandwich panel. For this reason, the dimensions of the cellular core is chosen so that the overall bending stiffness-to-mass ratio remains constant in the sandwich. It means that the modal density is expected to be identical for all the considered models at low frequency. The size of the unit-cells, d x and d y are detailed in Table 2.

Case (a) (b) (c) (d) (e) (f)

(f)

Figure 4: 3D Finite element model of the octagonal construction.

Effects on the transition frequency

The wave dispersion characteristics of the periodic waveguides defined by unit-cells (a)-(f) are obtained using the method described in Section 2. The flexural wavenumbers k are compared in Figure 5a. It is seen that the core geometry has a significant influence on the wavenumbers of the sandwich panel. More precisely, the flexural behaviours remain similar due to the constant bending rigidity and differ above the transition frequency. A 95% difference can then be observed between the (c) and (f) wavenumbers at 1800 Hz.

The transition frequency is located at the intersection between the asymptotic shear and bending curves. The method is shown in Figure 5b for the core geometry (b), where the bending k b and shear k s wavenumbers are both derived from the WFEM solution. It should be noted that an accurate bending wavenumber requires a refined frequency sampling in the lower frequencies. For the considered models, it can be seen that a clear bending behaviour is observed below 200 Hz, while the shear behaviour occurs above 600 Hz except for case (f).

For the honeycomb core (a), an exact analytical expression can be derived for the transition frequency. The dispersion relation in the sandwich can be written in terms of the mass per unit area m s , the equivalent transverse shear stiffness S and bending rigidity D of the plate [START_REF] Guillaumie | Vibroacoustic flexural properties of symmetric honeycomb sandwich panels with composite faces[END_REF]:

k 2 = m S 2S          ω 2 + ω ω 2 + 4S 2 m S D          (11)
Therefore, equalizing the asymptotic solutions defined Eqs.( 9) and (10) yields: The transition frequency ω can therefore be expressed analytically in terms of the equivalent bending, shear and mass parameters.

ω T = S √ m S D (12)
For other cores' geometries, it should be noted that the transition frequency can be only approximated using Eqs.( 9) and [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF], since Gibson and Ashby homogenisation is no longer valid. The equivalent parameters mS , D and S are obtained from Eq.( 11) and the numerical solution k t (ω) evaluated using non-linear least-square regression. The error on the regression function is written:

ε reg (ω) = k t (ω) -f reg (ω, D, mS , S ) k t (ω) (13)
The values of ε reg are compared in Figure 6. The regression appears in very good agreement with the WFEM solutions, except for Case (c) which exhibits up to 3% error in the considered frequency range. On the other hand Case (d) and (f) fit very well with the honeycomb model described in Eq. [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF]. It should be noted that the regression accuracy may depends on the chosen frequency range for f reg .

To summarize, the transition frequency can be determined from the asymptotic wavenumbers k s and k b derived from the WFEM solutions. In order to reduce the sensitivity to sampling for the determination of k b (ω → 0), the The modal density in the sandwich panel is calculated from the propagating wavenumbers and group velocities as mentioned above. The WFEM solutions are compared in Figure 8 for a surface S=1 m 2 , although it is reminded that the transition frequency does not depend on the panel's surface. The scaling procedure to preserve the mass-tostiffness ratios is validated since the modal densities n(ω → 0) are the same for the six models. It can be seen that n(ω) critically increases above ω T for each model. It should be noted that the influence of the core's geometry on the overall flexural behaviour will increase if the ratio between the core h T and the skins h s increases.

Concluding remarks

This paper is concerned with the design of sandwich panels with improved vibro-acoustic parameters. Although honeycomb constructions provide excellent compressional strength to density ratio, other types of periodic cores can be used to increase the transition frequency and reduce the modal density in a given frequency range. In this work, the WFEM was employed to compare various types of core constructions for a given sandwich panel.

This method does not suffer the limitations associated with core homogenisation and provides accurate wave dispersion characteristics with a reasonable computational cost. Additionally, the method can be further used to compute bandgaps and local resonances of the periodic cells which have a significant influence on the acoustic parameters in higher frequency (see Song et al. [START_REF] Song | Reduction of the sound transmission of a periodic sandwich plate using the stop band concept[END_REF]). It was shown that a 70 % increase of the transition frequency can be obtained by using an octagonal core geometry, while maintaining the same bending stiffness-to-mass ratio, as well as the skin and core materials. As a consequence of the increased equivalent shear modulus, the octagonal core also produces significantly higher group velocity and reduced modal density.

This transition frequency can be used as a relevant output for optimisation procedures oriented toward improving vibro-acoustic parameters. It should be noted that the proposed application is restrained to 2D geometric parameters (the core is unchanged along its thickness), since a reduced number of design parameters is required for sensitivity of optimisation procedures. An interesting development would be to understand how 3D cores could be designed in order to significantly postpone the bending to shear transition to higher frequencies.

Figure 1 :

 1 Figure 1: Top-view of the modelled periodic element with its edge (B, T, R, L), corner (1, 2, 3, 4) and inner (I) degrees of freedom.

3. 1 .

 1 Material properties and geometry of the core The periodic panel considered here is a symmetric sandwich plate involving composite faces made of carbon reinforced epoxy resin. The core is made of Nomex paper (aramid fibre and phenolic resin) of thickness 76.2 µm assumed isotropic. The skins are made with two plies of transverse isotropic, 0.25 mm-thick, unidirectional composite (UD) with the staking sequence [0 o , 90 o ]. The thickness of the core is h T = 25 mm and the material properties are detailed in Table 1. Material Density (kg.m -3 ) Young modulus (Pa) Shear modulus (Pa) Poisson ratio Nomex core ρ c = 1240 E c = 5.5 × 10 9 G c = 2.7 × 10 9 ν c = 0.33 E 1 = 133.6 × 10 9 G 23 = 2.6 × 10 9 ν 23 = 0.49 UD laminate ρ s = 1550 E 2 = 7.7 × 10 9 G 13 = 3.1 × 10 9 ν 13 = 0.29 E 3 = 7.7 × 10 9 G 12 = 3.1 × 10 9 ν 12 = 0.29

  2 and the corner parameters A x , A y , B x and B y . The other parameters are imposed by the compatibility requirement between the unit-cells. It can be seen that triangular, squares and pentagonal geometries can be obtained by setting coincident points M = P 1 , P 1 = P 2 or A x = B x and A y = B y .

Figure 2 :

 2 Figure 2: Top view of the 2D parametrized cores. Hexagonal (a) and octagonal (b) constructions.

Figure 3 :

 3 Figure 3: 3D Finite element models of the considered periodic cells.

Figure 5 :

 5 Figure 5: Flexural wavenumbers for the 5 different core configurations (a). Definition of the transition frequency : intersection between Case (b) asymptotic shear and bending waves (b).

Figure 7 :

 7 Figure 7: Averaged energy velocity in the sandwich waveguide.

Figure 8 :

 8 Figure 8: Modal density in the sandwich panels with different cores.

Table 1 :

 1 Mechanical properties of Nomex paper and UD fibre.

Table 2 :

 2 Size of the rectangular unit-cells in the six cases.

	d x (mm) 6.93	8	4.29	6.93 13.85 20
	d y (mm)	12	8	14.48 10.53	8	20

equivalent shear, bending and mass parameters can be obtained through a regression model. An analytical expression can therefore be derived for ω T with a good accuracy, except for case (c).

Influence on vibro-acoustic parameters

The transition can be observed on the averaged group velocity c g , as the passage from a square to linear function.

The velocities are compared in Figure 7. It can be seen that the cores have a significant influence on the energy velocity in the panel, and the waves are still dispersive above the transition frequency. Noteworthy, the maximal group velocity is obtained around the transition, meaning that reduced wave attenuations are expected in this frequency range.

The transition frequencies are detailed in Table 3 and compared with the honeycomb model (a). Case (c), associated with triangular construction has the lowest shear rigidity. The rectangular and parallelogram cores, cases (b) and (e) are associated with reduced transition frequencies compared to the hexagonal (a) core. The pentagonal core (d), produces a slight increase (9.11 %) of the transition frequency compared with honeycomb. On the other hand, the octagonal core (f) provides a considerable increase (70 %) of the transition frequency. Therefore, this geometry can postpone the transition frequency while maintaining the stiffness-to-mass ratio.