J.-M Mencik

A wave finite element-based formulation for computing the forced response of structures involving rectangular flat shells

Keywords: wave finite elements, rectangular shells, component mode synthesis

The harmonic forced response of structures involving several non-coplanar rectangular flat shells is investigated using the Wave Finite Element method. Such flat shells are connected along parallel edges where external excitation sources as well as mechanical impedances are likely to occur. Also, they can be connected to one or several coupling elements whose shapes and dynamics can be complex. The dynamic behavior of the connected shells is described by means of numerical wave modes traveling towards and away from the coupling interfaces. Also, the coupling elements are modeled using the conventional finite element method. A finite element mesh tying procedure between shells having incompatible meshes is considered which uses Lagrange Multipliers for expressing the coupling conditions in wave-based form. A global wavebased matrix formulation is proposed for computing the amplitudes of the wave modes traveling along the shells. The resulting displacement solutions are obtained using a wave mode expansion procedure. The accuracy of the wave-based matrix formulation is highlighted in comparison with the conventional finite element method through three test cases of variable complexities. The relevance of the method for saving large CPU times is emphasized. Its efficiency is also highlighted in comparison with the component mode synthesis technique.

INTRODUCTION

The necessity to reduce the CPU times involved by the conventional Finite Element (FE) method appears crucial in many applications to describe the low-and mid-frequency (LF and MF) forced response of structures. The Wave Finite Element (WFE) method constitutes an efficient alternative to address this task. The method is confined to the study of waveguides whose FE models are periodic [START_REF] Mead | A general theory of harmonic wave propagation in linear periodic systems with multiple coupling[END_REF], i.e., which are composed of identical substructures. These waveguides can be connected either directly or by means of coupling elements whose local dynamics can be complex [START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF]. In many cases, elastic waveguides with uniform cross-sections are dealt with. Within the WFE framework, the kinematic and mechanical fields of such structures are assessed in terms of numerical wave modes traveling in positive and negative directions. For each waveguide, the wave modes are computed by considering the mass and stiffness matrices of one single FE substructure (see above). The number of those wave modes is generally very small compared to the number of DOFs required by a full FE model when the whole waveguide is meshed. This results in matrix formulations of small sizes for expressing the forced response of structures [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF]. Proposing a general procedure that makes the WFE method applicable to the study of complex structures involving an arbitrary number of waveguides, connected either directly or by means of complex coupling elements, is the motivation behind the proposed study. The case of rectangular flat shells which are connected along parallel edges involving arbitrary excitation sources and mechanical impedances, or which are connected with coupling elements, is specifically dealt with. Those coupled elastic systems seem to constitute a reasonable approximation for describing many industrial structures, e.g., those encountered in the automotive sector (see Figure 1).

The WFE method belongs to the class of deterministic wave-based approaches which can be used to assess the MF behavior of shell and plate-like structures. Among these approaches are the analytic Trefftz techniques like the Wave Based Method (WBM) [START_REF] Vanmaele | An efficient wave based prediction technique for dynamic plate bending problems with corner stress singularities[END_REF] and the Variational Theory of Complex Rays (VTCR) [START_REF] Ladevèze | Calculation of medium-frequency vibrations over a wide frequency range[END_REF] which make use of propagative plane waves (as well as evanescent waves) to describe the dynamic behavior of bounded domains. The study of rectangular flat shell assemblies has been particularly addressed by the WBM in ref. [START_REF] Vanmaele | On the use of the wave based method for the steady-state dynamic analysis of three-dimensional plate assemblies[END_REF]. In this framework, the wave amplitudes are computed by considering a variational formulation over the boundaries (including coupling interfaces) of the shells. Otherwise, enrichments techniques such as the Partition of Unity Method (PUM) [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF] and the Discontinuous Enrichment Method (DEM) [START_REF] Farhat | The discontinuous enrichment method[END_REF] are other approaches which consider analytic wave functions as well as conventional Lagrange polynomials as interpolation functions of finite elements; these techniques have been proved to be efficient to improve the convergence of FE models over the MF range. The feature of these analytic wavebased approaches is that they can theoretically be applied to structures having arbitrary shapes and boundary conditions, contrary to the WFE method where rectangular flat shells are considered (even though these shells can be connected, as explained above). Also, these analytic approaches are not subject to the discretization errors of the conventional FE method as exact wave solutions of the local governing equations of structures are considered. Thus they are less susceptible to numerical dispersion than the WFE method. As another drawback of the WFE method, when compared for instance with the WBM or the VTCR, is the requirement that the boundary conditions are to be periodic (or uniform) along the direction of waveguides and that the excitation sources are to be considered at the ends of waveguides. However, these aforementioned analytic approaches have to face several drawbacks in the sense that a large number of plane waves is usually required to achieve their convergence (the issue lies in the description of the boundary conditions) while ill-conditioned full matrix systems are involved. In contrast, the WFE method appears as an efficient means to address these issues. The fact that the WFE method uses a relatively small number of wave modes -which travel along one direction and which reflect the cross-section dynamics of waveguidesyields the sizes of the related matrix formulations (and thus the CPU times) to be usually small compared to those involved by the analytic approaches. The other feature of the WFE method is that it involves well-conditioned matrix formulations, as discussed below.

It is important to understand that the WFE method is in theory as accurate as the conventional FE method to describe the dynamic behavior of structures involving waveguides. This is true provided that the same FE models are used for both approaches (i.e., periodic meshes for the waveguides). The advantage of the WFE method is that it enables the CPU times to be largely decreased compared to the FE method. From this point of view, the WFE method constitutes an interesting alternative to the usual model reduction strategies like the Component Mode Synthesis / Craig-Bampton (CMS/CB) method [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] where fixed interface modes are involved to describe the behavior of waveguides. The drawback of the CMS/CB method is that a large number of fixed interface modes can be required to reach its convergence. This is explained because the boundary and coupling conditions at the ends of waveguides can be complex, but also because the lengths of waveguides can be large (meaning that a large number of fixed interface modes can be required to describe the spatial dynamics). Such issues are not involved within the WFE framework due to the fact that the wave bases are full, i.e., basis truncation errors are not involved. Also, the sizes of the wave bases do not depend on the lengths of waveguides.

Several works have been made to compute the forced response of single or coupled straight waveguides by means of the WFE method. A spectral element-based strategy that uses WFE wave modes for calculating the dynamic stiffness matrices of bounded waveguides has been proposed in ref. [START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF]. This procedure has been validated to describe the harmonic forced response of simply supported Kirchhoff Love plates excited by one punctual force. More recently, a WFE approach has been proposed which uses the reflection/transmission coefficients of the wave modes at waveguide boundaries [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF]. The interesting feature of this approach is that it is based on the computation of the wave amplitudes, instead of explicitly considering the displacements/rotations and forces/moments. This yields simple and well-conditioned matrix formulations to be considered to describe the LF and MF forced responses, which is explained because the variations of wave amplitudes along waveguides are governed by means of diagonal matrices * . This strategy has been successfully used to describe the forced response of several kinds of straight waveguides (beam-like structures, plates, multi-layered systems) whose left and right limits are subject to Neumann and Dirichlet conditions; also, it has been used to describe the dynamic behavior of elastic systems involving two beam-like structures connected to an elastic junction (namely, a coupling element) undergoing local resonances [START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF]. In ref. [START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF], a model reduction technique based on the CMS/CB method has been proposed to describe the dynamic behavior of the coupling element in terms of a few fixed interface modes. Other WFE strategies which use the same idea (i.e., computation of the wave amplitudes by consideration of the wave reflections/transmissions at waveguide boundaries) have been proposed in refs. [START_REF] Waki | Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides[END_REF][START_REF] Mead | The forced vibration of one-dimensional multi-coupled periodic structures: An application to finite element analysis[END_REF]. In ref. [START_REF] Waki | Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides[END_REF], a procedure has been proposed to assess the dynamic behavior of single waveguides excited by one punctual force. The procedure consists in splitting any waveguide into two subwaveguides and considering the reflection/transmission of the wave modes at the location of the punctual force. This last works has been extended in ref. [START_REF] Renno | On the forced response of waveguides using the wave and finite element method[END_REF] to describe the dynamic behavior of infinite excited plates subject to fluid loading.

To summarize, the WFE method appears as an efficient means to describe the forced responses of elastic waveguides. The method enables these forced responses to be computed by means of small-sized and well-conditioned matrix formulations. So far, the method has been mostly applied to single waveguides involving Neumann and Dirichlet boundary conditions, or coupled structures involving two waveguides. It seems that a general procedure that makes the WFE method applicable to an arbitrary number of waveguides, connected either directly or by means of one or several coupling elements of arbitrary shapes, has not been investigated yet. Also, proposing a WFE procedure which takes into account arbitrary excitation sources and mechanical impedances at waveguide interfaces appears as another open challenge, never investigated previously. Such issues are addressed within the present study considering the problem of several non-coplanar rectangular flat shells connected either directly, along parallel edges involving arbitrary excitation sources and mechanical impedances, or by means of coupling elements. If one restricts the proposed approach to the case of two coplanar shells connected along one edge involving one punctual force, this yields the plate problem already investigated in refs. [START_REF] Waki | Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides[END_REF][START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF]. In a more general view, the proposed approach aims at investigating the dynamic behavior of complex structures such as those depicted in Figure 1.

The rest of the paper is organized as follows. In Section 2, the basics of the WFE method are recalled regarding the description of wave modes traveling along structures. Numerical experiments are made to highlight the wave modes traveling along rectangular flat shells. Also, the wave-based coupling conditions between several shells connected to one elastic coupling element are recalled in accordance with past studies [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF][START_REF] Ichchou | Wave finite elements for low and mid-frequency description of coupled structures with damage[END_REF]. In Section 3, the wave-based coupling conditions between two connected non-coplanar rectangular flat shells are derived. The procedure enables external excitation sources as well as mechanical impedances to be taken into account over the coupling interface. Also, a tying formulation based on Lagrange Multipliers is proposed which addresses the issue of shells having incompatible FE meshes. In Section 4, the strategy for computing the forced response of structures involving several connected shells and coupling elements is investigated. A wave-based matrix formulation is proposed for computing the amplitudes of the wave modes in each shell. In Section 5, numerical experiments are made to highlight the accuracy of the WFE strategy. The forced response of structures is investigated considering the following test cases: (1) three non-coplanar shells, with a punctual force, whose coupling conditions involve lineic densities of translational and rotational springs (Figure 9); (2) two sets of two non-coplanar shells, with translational springs subject to vertical displacements, connected to one coupling element having a non-uniform curvature (Figure 12); (3) six non-coplanar shells, with punctual forces, connected to one cylinder having a conical head (Figure 17). The relevance of the WFE method, in terms of accuracy and CPU time savings, is emphasized in comparison with reference solutions issued from the FE method. Also, a comparison with the conventional CMS/CB method (test case 2) is carried out to highlight further on the efficiency of the WFE method.

Figure 1. Structures involving several connected rectangular flat shells subject to punctual force (F 0 ), springs (K, K ′ ) with imposed displacements (q 0 , q ′ 0 ), lineic density of translational and rotational springs (k t , k ′ t , kr): (a) structure involving three shells when one shell has one edge clamped; (b) structure involving four shells when two shells are connected to one coupling element.

WFE METHOD

Wave propagation along rectangular flat shells

The WFE method aims at numerically describing the waves traveling along periodic structures [START_REF] Mead | A general theory of harmonic wave propagation in linear periodic systems with multiple coupling[END_REF]. Such structures are called periodic in the sense that they can be described by means of identical substructures which are connected along a main axis x, referred to as the direction of propagation. Rectangular flat shells which are meshed periodically along their length x belong to that class of structures. In the present study, these shells are supposed to be elastic, dissipative (considering a loss factor η) and subject to harmonic disturbance under frequency ω/2π (ω being the pulsation). The basic assumptions, for those thin elastic structures, is that both bending and in-plane motions that include drilling degrees of freedom [START_REF] Hughes | Dynamic analysis and drilling degrees of freedom[END_REF], are taken into account. In this framework, these structures are meshed by means of 2D triangular flat shell elements with three nodes and six degrees of freedom (DOFs) per node (i.e., three displacements and three rotations) that incorporate both bending actions [START_REF] Argyris | The TRIC shell element: theoretical and numerical investigation[END_REF] as well as membrane actions with drilling DOFs [START_REF] Allman | Implementation of a flat facet shell finite element for applications in structural dynamics[END_REF]. A rectangular flat shell with a periodic FE mesh is shown in Figure 2. The related FE substructures have the same length d while their left and right boundaries (i.e., the edges coincident with the z-direction) contain the same number of DOFs, namely n.

The WFE method requires the mass and stiffness matrices of any FE substructure to be known; it uses a state vector representation [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF] for linking the kinematic/mechanical fields between the left (or right) boundaries of two adjacent substructures k and k -1. In the frequency domain, this relationship is expressed in terms of a 2n × 2n symplectic matrix S as [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF]

u (k) = Su (k-1) , (1)
where u = [q T ± F T ] T , q and F being the n × 1 vectors of displacements/rotations and forces/moments, respectively. The sign ahead of F in u results from the convention made for expressing the forces on the left or right boundaries of the substructures: in the present study, the convention -F (resp. F) is used to denote the left (resp. right) substructure boundaries. In Eq. ( 1), the matrix S is constructed from the condensed dynamic stiffness matrix of any of these substructures with regard to its left and right boundaries (cf. [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF] for further details).

The wave modes refer to the eigenvalues and eigenvectors of S which are denoted as {µ j } j and {Φ j } j , respectively. It must be noticed that S is expressed from the dynamic stiffness matrix of a substructure (see above), meaning that the wave modes depend on the frequency. According to Bloch's theorem [START_REF] Mencik | A model reduction strategy for computing the forced response of elastic waveguides using the wave finite element method[END_REF], the eigenvalues {µ j } j can be expressed as {e -iβj d } j , where {β j } j have the meaning of wave numbers. Also, the eigenvectors {Φ j } j , also known as wave shapes, relate the spatial distributions of the kinematic and mechanical fields over the width of the shell (z-direction). The fact that S is symplectic (see above) yields the wave modes {(µ j , Φ j )} j to be split into n incident and n reflected wave modes, i.e., n waves traveling towards and n waves traveling away from the right (or left) boundary of the shell. These incident and reflected wave modes are denoted as {(µ inc j , Φ inc j )} j and {(µ ref j , Φ ref j )} j ; they are usually defined so that |µ inc j | < 1 and |µ ref j | > 1 ∀j. † . Otherwise, the vectors of displacements/rotations q and forces/moments ±F, over any substructure boundary k (i.e., either a coupling interface between two consecutive substructures k -1 and k, or a limiting edge of the shell), can be expanded as [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF]

q (k) = Φ inc q Q inc(k) + Φ ref q Q ref(k) , ±F (k) = Φ inc F Q inc(k) + Φ ref F Q ref(k) , (2)
where

Φ inc q , Φ ref q , Φ inc F and Φ ref F
are square n × n matrices constituted from the displacement/rotation and force/moment components of the incident and reflected wave shapes; also, Q inc(k) and Q ref (k) are n × 1 vectors of wave amplitudes, whose variation along the shell is governed as [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF]

Q inc(k) = µQ inc(k-1) , Q ref(k) = µ -1 Q ref(k-1) , ( 3
)
where µ is a n × n matrix defined as µ = diag{µ inc j } j , such that ||µ|| 2 < 1 (||.|| 2 being the 2-norm) ‡ .

Numerical results

The wave modes of an aluminum rectangular shell are computed over a frequency band † Such a consideration follows from the fact that S is a symplectic matrix -i.e., its eigenvalues come in pairs as (µ, 1/µ) [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF] -while it is assumed that the shell is damped.

‡ The fact that ||µ|| 2 < 1 is readily proved since |µ inc j | < 1 ∀j.

[10 Hz, 300 Hz]. The characteristics of the structure are: density ρ = 2700 kg.m -3 , Young's modulus E = 70 × 10 9 Pa, Poisson's ratio ν = 0.3, loss factor η = 0.01, width (z-direction) L z = 1 m and thickness (y-direction) h = 5 × 10 -3 m. Within the WFE framework, the sought waves are those traveling along the length of the shell, namely the x-direction. The wave modes are computed by considering the FE model of a substructure of length d = 0.05 m with 21 nodes uniformly spread on its width, as shown in Figure 2 (here, the distance between two consecutive nodes is 0.05 m). The choice of this length 0.05 m between two consecutive nodes follows from the well known rule that consists in discretizing the wavelengths at least by means of 8 finite elements. To this end, one focuses on the wavelength of the bending wave traveling in an infinite equivalent Reissner-Mindlin plate, which is known analytically, at the maximum frequency considered (i.e., 300 Hz). Actually, it is assumed that this wavelength represents a minimum value regarding all the modes traveling along the shell. In other words, the fact to discretize this analytic wavelength at least by 8 elements yields the same conclusion regarding the other wavelengths. The FE model of the substructure is performed using the COMSOL Multiphysics R software. It is worth recalling that each triangular element incorporates both bending actions (with a shear correction factor κ = 5/6) and membrane actions with drilling DOFs. The fact that the left boundary (as well as the right boundary) contains n = 21 × 6 = 126 DOFs means that a same number of incident/reflected wave modes are to be computed. This is achieved by means of MATLAB R using the procedure proposed in ref. [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF] that consists in solving a generalized eigenproblem for the displacements/rotations only. The procedure requires one to get the mass and stiffness matrices of the substructure (and ultimately, the dynamic stiffness matrix condensed over the left and right boundaries) which is done by means of COMSOL Multiphysics R . The wave modes are computed at 581 discrete frequencies {f k } k uniformly spread on the frequency band [10 Hz, 300 Hz], considering an identical frequency step ∆f . Part of these wave modes are highlighted in Figure 3 regarding the dispersion curves -i.e., the real and imaginary parts of the wave numbers -over the frequency band. To track any wave mode r over the frequency domain, a correlation criterion among the wave shapes at two consecutive discrete frequencies f k and f k + ∆f is proposed as

max s |(Φ qs ) H f k +∆f (Φ qr ) f k | ||(Φ qs ) f k +∆f || 2 ||(Φ qr ) f k || 2 + |(Φ Fs ) H f k +∆f (Φ Fr ) f k | ||(Φ Fs ) f k +∆f || 2 ||(Φ Fr ) f k || 2 , (4)
where the superscript H denotes the conjugate transpose; Φ qr and Φ Fr (resp. Φ qs and Φ Fs ) represent the vectors of the displacement/rotation components and the force/moment components of the wave shape Φ r (resp. Φ s ), respectively. Clearly, once a wave mode has been numbered as r at the discrete frequency f k , the proposed criterion enables one to select among all the wave modes computed at the subsequent discrete frequency f k + ∆f the one which matches the best (from the point of view of the displacement/rotation components as well as the force/moment components) this wave mode r at the discrete frequency f k . Thus the procedure consists in renumbering the selected wave mode as r for the discrete frequency f k + ∆f , and so on. In other words, this criterion yields the wave modes to be clearly identified over the frequency domain (see ref. [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF] for further discussions).

The wave modes can be classified as propagating (i.e., the imaginary parts of the wave numbers are close to zero), evanescent (i.e., the real parts of the wave numbers are close to zero) or complex (i.e., the real and imaginary parts of the wave numbers are of the same order). For the sake of clarity, the complex wave modes are not shown in Figure 3. The wave shapes are generally disparate, as seen in Figure 4 where nine modes are considered at 10Hz (the arrow indicates the direction of wave propagation). Wave shapes (a), (b), (c) and (d) highlight the conventional longitudinal, flexural (in the zand y-directions) and torsional modes. Other wave shapes represent high order modes, part of them becoming propagating at certain frequencies (see Figure 3). It is worth emphasizing that the wave shapes are subject to changes as the frequency increases, as seen in Figure 5 at 300 Hz.

(a) (b) (c) (d) (e) (f) (g) (h) (i)

Wave-based coupling conditions between several shells connected to one coupling element

The issue of several waveguides connected to one elastic coupling element has been treated in depth in the literature (see for instance refs. [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF][START_REF] Mencik | Propagation multimodale dans les systèmes périodiques couplés (in French)[END_REF][START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF]). The related wave-based coupling conditions are recalled here for the sake of clarity. A kind of problem which is addressed within this topic is depicted in Figure 6, considering several flat shells having different orientations {(x i , y i , z i )} i ,

(a) (b) (c) (d) (e) (f) (g) (h) (i)
Figure 5. Several wave shapes (real parts of the displacement components) at 300 Hz.

connected to one coupling element of arbitrary shape. Within the WFE framework, the dynamic behavior of the shells is described by means of wave modes (see Section 2.1) while the local dynamics of the coupling element are assessed using the standard FE method. The procedure uses a mesh tying formulation that connects each shell with the coupling element. Lagrange Multipliers are invoked to address incompatible meshes between the shells and the coupling element. Here, these Lagrange Multipliers are discretized with respect to the FE mesh of the coupling element [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF]. For the sake of clarity, the present study will be restricted to the consideration of coupling elements which are free from excitation sources (the case of coupling elements undergoing imposed forces and prescribed displacements has been treated in depth in ref. [START_REF] Mencik | Propagation multimodale dans les systèmes périodiques couplés (in French)[END_REF] and may be considered without difficulty). According to ref. [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF], the wave-based coupling conditions between a set of R rectangular flat shells connected to one coupling element are expressed as

   Q ref 1 . . . Q ref R    = C c    Q inc 1 . . . Q inc R    , (5)
where

Q inc i and Q ref i
denote the vectors of wave amplitudes at the coupling interface, regarding incident/reflected wave modes traveling along each shell i (i = 1, . . . , R). The matrix C c is expressed as

C c = -(A ref ) -1 A inc , (6)
where [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF] A

ref = T T D * T Ψ ref q + Ψ ref F , A inc = T T D * T Ψ inc q + Ψ inc F . (7)
Here, D * is the dynamic stiffness matrix of the coupling element condensed on the interface DOFs;

T is a block diagonal matrix with block components {B c i } i which reflect the use of Lagrange Multipliers to describe the coupling conditions between the shells and the coupling element; also,

Ψ inc q , Ψ ref q , Ψ inc F and Ψ ref F
are square matrices of same size expressed as

Ψ inc q =    L c 1 (Φ inc q ) 1 • • • 0 . . . . . . . . . 0 • • • L c R (Φ inc q ) R    , Ψ ref q =    L c 1 (Φ ref q ) 1 • • • 0 . . . . . . . . . 0 • • • L c R (Φ ref q ) R    , (8)
Ψ inc F =    L c 1 (Φ inc F ) 1 • • • 0 . . . . . . . . . 0 • • • L c R (Φ inc F ) R    , Ψ ref F =    L c 1 (Φ ref F ) 1 • • • 0 . . . . . . . . . 0 • • • L c R (Φ ref F ) R    ,
where, for each shell i (i = 1, . . . , R), L c i is a square n i × n i matrix (n i being the number of DOFs contained over the left/right boundary of the shell) which expresses the displacements/rotations and forces/moments of the shell in the coordinate system of the coupling element (x c , y c , z c ) (see Figure 6) [START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF].

Remark

When deriving Eq. ( 5), the shells are assumed to be connected to the coupling element over their right boundaries: in other words, the local frame axis x i of each shell i is assumed to point towards the coupling element (see Figure 6). Considering the left boundaries of some shells does not bring additional difficulties, however.

WAVE-BASED COUPLING CONDITIONS BETWEEN TWO CONNECTED SHELLS

Mesh tying formulation

Two structures involving one or several sets of rectangular flat shells, connected along parallel edges, are depicted in Figure 1. Within the WFE framework, the shells are meshed periodically in terms of identical FE substructures (see Section 2.1). The related FE models are illustrated in Figure 7 regarding two connected shells i and i + 1. These shells are described by means of two local reference frames R i = (x i , y i , z) and R i+1 = (x i+1 , y i+1 , z) whose orientations in the global reference frame -namely, R = (x, y, z) -are expressed around the z-direction in terms of two angles α i and α i+1 . The coupling interface -namely, Γ -represents the right boundary of shell i as well as the left boundary of shell i + 1. In the present study, it is supposed that the displacement/rotation field of the shells are continuous across Γ . Also, it is supposed that punctual and lineic excitation sources (forces/moments), as well as mechanical impedances (e.g., translational or rotational stiffnesses), are likely to occur over Γ . The substructures used to describe the two shells are enabled to have different lengths d i and d i+1 , as well as different FE discretizations. This yields different numbers of DOFs n i and n i+1 to describe the boundaries of shells i and i + 1. Here, it is assumed without loss of generality that shell i exhibits the finest mesh compared to shell i + 1 (see Figure 7). The issue is to transcribe the coupling conditions within the FE framework by considering incompatible meshes across the coupling interface Γ . Lagrange Multipliers are used to treat this issue. In this framework, a Lagrange Multiplier field is defined on Γ which is discretized using the trace of the mesh of either shell i or shell i + 1. In the present study, it is proposed to choose the finest mesh (shell i), the motivation behind this choice being that the spatial behavior of shell i can be correctly mapped on shell i + 1 by means of a sufficient number of constraints. It is worth pointing out that the consideration of the coarser mesh (shell i + 1), instead of the finest mesh, would have been problematic regarding zero-energy modes § . Considering such a discretization for the Lagrange Multiplier field yields the continuity of the displacement/rotation field across Γ to be expressed in the global reference frame R as

B i q Γ i -B i+1 q Γ i+1 = 0, (9)
where q Γ i and q Γ i+1 denote the vectors of displacements/rotations, resp. for shell i and shell i + 1, expressed on Γ ; also, B i and B i+1 are real matrices of the forms

B i = Γ N T p N q Γ i dS and B i+1 = Γ N T p N q Γ i+1 dS
(the superscript T denoting the matrix transpose); N p is the matrix of interpolation functions which discretize the Lagrange Multiplier field on Γ ; N q Γ i and N q Γ i+1 are the matrices of interpolation functions which discretize the displacement/rotation fields of shells i and i + 1 on Γ , respectively (see ref. [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF] for further details). If one assumes that N q Γ i and N p are equal yields the matrix B i to be square positive definite and thus invertible. According to Eq. ( 9), § Such an issue occurs considering for instance the case of a highly oscillating spatial behavior of side i with zero displacement points coincident with the DOFs of side i + 1, yielding a null displacement solution of shell i + 1.

this yields q Γ

i to be uniquely determined from the knowledge of q Γ i+1 as

q Γ i = Bq Γ i+1 , (10)
where B = B -1 i B i+1 . Also, the vectors of forces/moments expressed on Γ are linked as

F Γ i+1 = B T (-F Γ i + F ex -iωZq Γ i ), (11)
A justification on how Eq. ( 11) has been derived can be found in ref. [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF]. Here, F Γ i and F Γ i+1 are the vectors of forces/moments expressed on Γ , resp. for shell i and shell i + 1; also, F ex is a vector of excitation sources defined on Γ with respect to the mesh of shell i; ω is the pulsation; finally, Z is an impedance matrix defined with respect to the mesh of shell i, which reflects punctual or lineic springs (or masses, . . . ) on Γ . The fact that those springs can be excited by imposed displacements, instead of being clamped on one of their ends, can be taken into account without difficulty by expressing the vector F ex in a suitable way. To summarize, the coupling conditions between two non-coplanar rectangular flat shells i and i + 1 are expressed by means of Eqs. [START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF] and [START_REF] Waki | Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides[END_REF]. It is worth emphasizing that these relations are expressed in the global reference frame R = (x, y, z). Expressing these relations in the local reference frames R i = (x i , y i , z) and R i+1 = (x i+1 , y i+1 , z) of the shells enables the vectors of displacements/rotations and forces/moments to be expressed in terms of wave modes, as explained in Section 2 (Eq. ( 2)). Such a procedure is considered here for expressing the coupling conditions in wave-based form. For this purpose, it is proposed to introduce two rotation matrices L i and L i+1 , defined as

(q Γ i ) Ri = L i q Γ i , (F Γ i ) Ri = L i F Γ i , ( 12
) (q Γ i+1 ) Ri+1 = L i+1 q Γ i+1 , (F Γ i+1 ) Ri+1 = L i+1 F Γ i+1 . (13)
The matrices L i and L i+1 are constructed from the direction cosines of the local frames R i and R i+1 , respectively, in the global frame R. These matrices are real and orthogonal, i.e.,

L T i L i = L i L T i = I and L T i+1 L i+1 = L i+1 L T i+1 = I.
Introducing Eqs. ( 12) and ( 13) in Eqs. [START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF] and (11) while invoking the aforementioned orthogonality properties of L i and L i+1 enables the coupling conditions to be expressed as

(q Γ i ) Ri = L i BL T i+1 (q Γ i+1 ) Ri+1 , (14)
and

(F Γ i+1 ) Ri+1 = -L i BL T i+1 T (F Γ i ) Ri + L i iωZL T i (q Γ i ) Ri + BL T i+1 T F ex . (15)
Eqs. ( 14) and [START_REF] Ichchou | Wave finite elements for low and mid-frequency description of coupled structures with damage[END_REF] express the coupling conditions between the two shells, where the vectors of displacements/rotations and forces/moments are expressed in the local reference frames R i and R i+1 . Expanding these vectors in terms of the wave modes of the two shells yields the wave-based form of the coupling conditions. The procedure is detailed hereafter.

Wave-based coupling conditions

The wave mode expansion of the vectors of displacements/rotations and forces/moments follows from Eq. ( 2), where the superscripts inc and ref denote the wave modes which travel towards and away from the right boundary of any shell (see Figure 2). Regarding for instance shell i, the matrices of incident/reflected wave modes are written as

(Φ inc ) i = [(Φ inc q ) T i (Φ inc F ) T i ] T and (Φ ref ) i = [(Φ ref q ) T i (Φ ref F ) T i ],
while the related vectors of wave amplitudes are denoted as

Q inc(k) i and Q ref(k) i
. Should the left boundary of the shell (instead of its right boundary) be considered as a reference, the matrices of incident/reflected wave modes are denoted as (Φ inc⋆ ) i and (Φ ref⋆ ) i , while the vectors of wave amplitudes are denoted as

Q inc⋆(k) i and Q ref⋆(k) i
. These matrices and vectors are simply expressed as

(Φ inc⋆ ) i = (Φ ref ) i , (Φ ref⋆ ) i = (Φ inc ) i , Q inc⋆(k) i = Q ref(k) i , Q ref⋆(k) i = Q inc(k) i . ( 16
)
Such a convention involving the right or left boundary of the shell is introduced here as a means to clarify and simplify the subsequent developments made in the paper. The convention is highlighted in Figure 8. Invoking the aforementioned convention yields the wave mode expansion to be written for both shells i and i + 1 as

(q Γ i ) Ri = (Φ inc q ) i Q inc i + (Φ ref q ) i Q ref i , (F Γ i ) Ri = (Φ inc F ) i Q inc i + (Φ ref F ) i Q ref i , (17)
and

(q Γ i+1 ) Ri+1 = (Φ inc⋆ q ) i+1 Q inc⋆ i+1 + (Φ ref⋆ q ) i+1 Q ref⋆ i+1 , -(F Γ i+1 ) Ri+1 = (Φ inc⋆ F ) i+1 Q inc⋆ i+1 + (Φ ref⋆ F ) i+1 Q ref⋆ i+1 , (18)
where

(q Γ i ) Ri and (q Γ i+1 ) Ri+1 (resp. (F Γ i ) Ri and (F Γ i+1 ) Ri+1
) are the vectors of displacements/rotations (resp forces/moments) introduced in Section 3.1. The minus sign ahead of (F Γ i+1 ) Ri+1 in Eq. ( 18) follows from the convention mentioned below Eq. (1). Also,

Q inc i and Q ref i (resp. Q inc⋆ i+1 and Q ref⋆ i+1
) denote the vectors of wave amplitudes for shell i (resp. shell i + 1) at the coupling interface Γ . Introducing Eqs. ( 17) and [START_REF] Allman | Implementation of a flat facet shell finite element for applications in structural dynamics[END_REF] in Eqs. ( 14) and ( 15) leads to the following matrix system (see Appendix A for further details):

A ref Q ref i Q ref⋆ i+1 = -A inc Q inc i Q inc⋆ i+1 + 0 -(Φ ref⋆ F ) -1 i+1 BL T i+1 T F ex , (19)
where A ref and A inc are two (n i + n i+1 ) × (n i + n i+1 ) ¶ matrices expressed as

A ref =     I ni -(Φ ref q ) -1 i L i BL T i+1 (Φ ref⋆ q ) i+1 -(Φ ref⋆ F ) -1 i+1 L i BL T i+1 T × (Φ ref F ) i + L i iωZL T i (Φ ref q ) i I ni+1     , (20)
and

A inc =     (Φ ref q ) -1 i (Φ inc q ) i -(Φ ref q ) -1 i L i BL T i+1 (Φ inc⋆ q ) i+1 -(Φ ref⋆ F ) -1 i+1 L i BL T i+1 T × (Φ inc F ) i + L i iωZL T i (Φ inc q ) i (Φ ref⋆ F ) -1 i+1 (Φ inc⋆ F ) i+1     . ( 21
)
In Eqs. [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF][START_REF] Mencik | A model reduction strategy for computing the forced response of elastic waveguides using the wave finite element method[END_REF][START_REF] Mencik | Propagation multimodale dans les systèmes périodiques couplés (in French)[END_REF],

(Φ ref q ) -1 i and (Φ ref⋆ F ) -1
i+1 are the inverses of the matrices (Φ ref q ) i and (Φ ref⋆ F ) i+1 , respectively. The existence of such inverses results from the fact that (Φ ref q ) i and (Φ ref⋆ F ) i+1 are full column rank (a proof of this statement is given in ref. [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF]). Invoking these inverses results in a better numerical conditioning of the matrix A ref . In fact, this circumvents the numerical issue involved when the matrix is partitioned into displacement/rotation and force/moment components [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF]. Also, each reflected or incident wave mode

Φ r = [Φ T qr Φ T Fr ] T is
) i+1 cannot lead to (L i BL T i+1 )(Φ ref⋆ q ) i+1 = 0, i.e., null(L i BL T i+1 ) = {0}, and (ii) ran(L i BL T i+1 ) cannot intersect the null space of (Φ ref q ) -1 i . Condition (i) is readily verified since L i BL T
i+1 is full column rank; this is proved taking into account that the Lagrange Multiplier field is discretized with regard to the finest mesh (side i); this would have failed otherwise (see comments at the beginning of Section 3.1). Also, condition (ii) is verified because (Φ ref q ) -1 i is full rank (the proof is given in ref. [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF]), i.e., null((Φ ref q ) -1 i ) = {0}. The proof that A ref 21 is full rank can be achieved on the same way as above, provided that the matrix

[(Φ ref F ) i + L i iωZL T i (Φ ref q ) i
] is full rank [START_REF] Mencik | On the low-and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation[END_REF]. This seems to be verified in general, taking into account that (Φ ref F ) i and (Φ ref q ) i are full rank matrices, except maybe in some very particular cases (e.g., when

L i iωZL T i (Φ ref q ) i = -(Φ ref F )
i ) which will be not considered here.

Thus, taking into account that the matrix A ref is invertible yields the solution of Eq. ( 19) to be expressed as

Q ref i Q ref⋆ i+1 = C Q inc i Q inc⋆ i+1 + F, (22)
where

C = -(A ref ) -1 A inc , F = (A ref ) -1 0 -(Φ ref⋆ F ) -1 i+1 BL T i+1 T F ex . ( 23
)
Eq. ( 22) expresses the wave-based coupling conditions between shells i and i + 1. Here, C is a square (n i + n i+1 ) × (n i + n i+1 ) matrix whose components are to be understood as the ¶ n i and n i+1 are the numbers of DOFs contained over the left/right boundaries of shells i and i + 1, respectively.

reflection/transmission coefficients of the wave modes (considering both shells i and i + 1) at the coupling interface Γ ; otherwise, F is a (n i + n i+1 ) × 1 vector that plays the role of excitation sources.

Remark: case when shell i + 1 has the finest mesh In this case, the matrix B is defined as B = B -1 i+1 B i . It can be shown without difficulty that the relationships ( 22) and (23) still hold provided that A ref and A inc are expressed as

A ref =     I ni -(Φ ref F ) -1 i L i+1 BL T i T × (Φ ref⋆ F ) i+1 -L i+1 iωZL T i+1 (Φ ref⋆ q ) i+1 -(Φ ref⋆ q ) -1 i+1 L i+1 BL T i (Φ ref q ) i I ni+1     , (24)
A inc =     (Φ ref F ) -1 i (Φ inc F ) i -(Φ ref F ) -1 i L i+1 BL T i T × (Φ inc⋆ F ) i+1 -L i+1 iωZL T i+1 (Φ inc⋆ q ) i+1 -(Φ ref⋆ q ) -1 i+1 L i+1 BL T i (Φ inc q ) i (Φ ref⋆ q ) -1 i+1 (Φ inc⋆ q ) i+1     . (25)
Also, the vector of excitation sources F is to be expressed as

F = (A ref ) -1 0 (Φ ref F ) -1 i BL T i T F ex . (26)

FORCED RESPONSE COMPUTATION

Preliminary comments and conventions

The harmonic forced response of structures involving several rectangular flat shells and coupling elements is investigated. Within the WFE framework, the computation of the forced response requires a set of wave-based coupling and boundary conditions to be expressed over the left/right boundaries of the shells. This yields a global wave-based matrix formulation to be considered which provides the wave amplitudes, and ultimately the kinematic and mechanical fields, for each shell. This strategy is depicted in the rest of the section.

As opening remarks of this study, it is worth emphasizing that any shell i is supposed to be composed from an integral number N i of identical substructures. This yields the number of substructure boundaries (i.e., the coupling interfaces between the substructures as well as the two limiting edges) involved along the shell to be N i + 1. As a convention, the substructures are numbered from the left to the right boundary of the shell. Also, the following notations are introduced:

Q inc i = Q inc(Ni+1) i , Q ref i = Q ref(Ni+1) i , Q inc⋆ i = Q inc⋆(1) i , Q ref⋆ i = Q ref⋆(1) i , (27)
where

{Q inc i , Q ref i }

Wave-based description of the global structure

The wave-based coupling conditions between two shells i and i + 1, connected along a common interface Γ i i+1 , have been derived in Section 3. They are expressed by Eq. ( 22), where the matrix C and the vector F are to be formulated using either Eqs. [START_REF] Mencik | A model reduction strategy for computing the forced response of elastic waveguides using the wave finite element method[END_REF][START_REF] Mencik | Propagation multimodale dans les systèmes périodiques couplés (in French)[END_REF] and (23) (case when shell i has the finest mesh), or Eqs. (24-26) (case when shell i + 1 has the finest mesh). The wave-based coupling conditions (22) are to be expressed for each interface Γ i i+1 considered. In block matrix form, this yields

Q ref i Q ref⋆ i+1 = C ii C ii+1 C ⋆ i+1i C ⋆ i+1i+1 Q inc i Q inc⋆ i+1 + F i F ⋆ i+1 , (28)
where C ii and C ⋆ i+1i+1 are square n i × n i and n i+1 × n i+1 matrices whose components denote the reflection coefficients of the wave modes, respectively for shells i and i + 1; C ii+1 and C ⋆ i+1i are rectangular n i × n i+1 and n i+1 × n i matrices whose components denote the transmission coefficients among the wave modes of the two shells; also, F i and F ⋆ i+1 are n i × 1 and n i+1 × 1 vectors of excitation sources (resp. for shells i and i + 1).

The wave-based coupling conditions between several shells connected to one coupling element have been derived in Section 2.2 (Eq. ( 5)). In the general case when several coupling elements are dealt with, the strategy consists in numbering as t sr (r = 1, . . . , R s ) a given shell connected to one coupling element s (s = 1, . . . , S). In block matrix form, this yields

   Q ref ts1 . . . Q ref tsR s    =    C c ts1ts1 • • • C c ts1tsR s . . . . . . . . . C c tsR s ts1 • • • C c tsR s tsR s       Q inc ts1 . . . Q inc tsR s    . (29)
When deriving this equation, the shells are implicitly assumed to be connected to the coupling element over their right boundaries. Considering their left boundaries does not bring difficulties, as mentioned below Eq. ( 8).

Apart from the coupling conditions, the left and right boundaries of any shell i may involve imposed forces/moments as well as imposed displacements/rotations, rather than coupling conditions, which also need to be considered. The related wave-based boundary conditions are expressed in Appendix B, considering imposed forces/moments and imposed displacements/rotations, as well as mechanical impedances. These boundary conditions are deduced from Eq. ( 11), taking into account that F Γ i+1 = 0. The resulting expressions, expressed at the left boundary and the right boundary of the shell, are

Q ref⋆ i = C ⋆ ii Q inc⋆ i + F ⋆ i , Q ref i = C ii Q inc i + F i , (30)
where C ⋆ ii and C ii are square n i × n i matrices expressing the reflection coefficients, while F ⋆ i and F i stand for the vectors of excitation sources (cf. Appendix B).

As a final comment, according to Eqs. ( 3) and ( 27), the vectors of wave amplitudes (expressed at the left and right boundaries of any shell i) are linked as

Q inc i = µ Ni i Q ref⋆ i , Q inc⋆ i = µ Ni i Q ref i , (31)
where µ i is the diagonal matrix of the wave mode parameters {(µ inc j ) i } j , which is such that ||µ i || 2 < 1 (cf. Section 2.1).

Matrix formulation

Considering the set of equations mentioned in Section 4.2, a global matrix formulation can be expressed which enables the computation of the wave amplitudes for each shell. To this aim, it is proposed to number as 1, 2, . . . all the rectangular flat shells involved to describe the whole structure. Two consecutive shells are numbered as i and i + 1. Also, it is proposed to denote as , where s (s = 1, . . . , S) relates any coupling element number. As a result, according to Eq. (31), Eqs. (28-30) yield the following system of equations:

Q ref⋆ i = C ⋆ ii µ Ni i Q ref i + C ⋆ ii-1 µ Ni-1 i-1 Q ref⋆ i-1 + F ⋆ i ∀i ∈ (O 1 ) L , (32)
Q ref i = C ii µ Ni i Q ref⋆ i + C ii+1 µ Ni+1 i+1 Q ref i+1 + F i ∀i ∈ (O 1 ) R , (33)
Q ref tsr = Rs q=1 C c tsrtsq µ Nt sq tsq Q ref⋆ tsq ∀t sr ∈ O s 2 r = 1, . . . , R s s = 1, . . . , S, (34)
Q ref⋆ i = C ⋆ ii µ Ni i Q ref + F ⋆ i ∀i ∈ (O 3 ) L , (35)
Q ref i = C ii µ Ni i Q ref⋆ i + F i ∀i ∈ (O 3 ) R . (36)
Invoking the aforementioned equations yields a matrix formulation of the form AQ = F, where Q and F are given by

Q =      Q ref⋆ 1 Q ref 1 Q ref⋆ 2 . . .      , F =      F ⋆ 1 F 1 F ⋆ 2 . . .      . (37)
Here, A and F are, respectively, a 2 i n i × 2 i n i matrix and a 2 i n i × 1 vector whose exact expressions depend on the kind of structures investigated. These will be specified in the next section considering different test cases. The solution of the matrix formulation provides the wave amplitudes as Q = A -1 F (this is true provided that A is invertible). The matrix formulation is to be computed at every discrete frequency considered within the frequency band studied. The determination of the vectors of displacements/rotations q (k) i and forces/moments F (k) i at any substructure boundary k (k = 1, . . . , N i + 1) along any shell i follows from Eq. ( 2), where the vectors

Q inc(k) i and Q ref(k) i
are to be expressed as (cf. Eqs. ( 3) and ( 27)):

Q inc(k) i = µ k-1 i Q ref⋆ i , Q ref(k) i = µ Ni-(k-1) i Q ref i k = 1, . . . , N i + 1. (38)
Again, it is important to understand that the solution provided by the WFE method -i.e., by solving the aforementioned matrix system AQ = F and using the wave mode expansion (2) -gives in theory exactly the same result as the conventional FE method (see Section 1). This is true provided that both FE and WFE approaches are based on the same FE model for the whole structure, i.e., with periodic meshes for the shells. Should the related FE meshes be different, FE and WFE solutions may be slightly different.

Remark: CPU times

The WFE method requires a matrix A (see above) of size 2 i n i × 2 i n i (2 i n i being the total number of DOFs involved over the left and right boundaries of the shells) to be inverted at several discrete frequencies for computing the forced response of structures. In case when coupling elements are considered, it also requires the dynamic stiffness matrices of these elements to be condensed on their interfaces DOFs at each discrete frequency: such a procedure is achieved by inverting particular block components of these matrices [START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF]. In comparison, the classic FE method requires a dynamic stiffness matrix with a larger size -i.e., which includes all the DOFs of the shells and eventual coupling elements -to be inverted at each frequency. Taking into account that the number of DOFs used to describe the shells is usually much greater compared to those contained over their left and right boundaries, this yields the size of FE-based matrices to be considerably large compared to those involved by the WFE method. Considering that the number of matrix inversions is usually large (i.e., many discrete frequencies can be considered for expressing the dynamic behavior of the structure), this yields the CPU times involved by the WFE method to be very small compared to those of FE.

Numerical conditioning

The wave-based matrix formulation involved by Eqs. (32-36) is expressed as AQ = F whose solution is Q = A -1 F. The matrix formulation is well-conditioned provided that the condition number of A is small enough. This feature is highlighted as follows. The matrix A involves identity matrices as diagonal block terms, while its off-diagonal block terms are of the form -Cµ N (N denoting the number of substructures in shells). The feature of this matrix form is that ||µ|| 2 < 1 (see below Eq. ( 3)), meaning that A tends to the identity matrix when N increases. In other words, the matrix is likely to be inverted without difficulty when a sufficient number of substructures are dealt with. The second feature is the use of µ N in -Cµ N which results in a filtering effect for high order wave modes -i.e., those for which the parameters {µ j } j are close-to-zero -that do not contribute for expressing the forced response of the structure. Those non-contributing modes may be understood as spurious solutions of the WFE method which can be sources of numerical problems. The fact that such modes can be filtered through the consideration of matrix terms of the form -Cµ N in A is explained as follows. Let us assume that the matrices C and µ are partitioned as

C = C E 1 E 2 C , µ = µ 0 0 µ , (39)
where the matrices E 1 and E 2 , as well as the matrix µ, have close-to zero components. Here, µ relates the diagonal matrix of the wave parameters associated to the aforementioned spurious modes, while E 1 and E 2 reflect coupling matrices for linking those spurious modes to the other wave modes.

The fact that E 1 and E 2 have close-to-zero components is explained since the spurious modes are decoupled from the modes which contribute for expressing the forced response of the structure (this is understood since the spurious modes do not contribute to the forced response, by definition). As a result, right multiplying the matrix C by µ N gives

Cµ N = C µ E 1 µ E 2 µ C µ ≈ C µ 0 0 0 . ( 40
)
Here, it is assumed that the matrix C µ can be neglected compared to C µ, because the matrix µ has close-to zero components (see above). To summarize, the strategy consisting in expressing the matrix formulation AQ = F by means of Eqs. (32-36) yields the influence of high order spurious wave modes to be considerably lowered when computing the forced responses. This makes the proposed matrix formulation relatively well-conditioned.

NUMERICAL EXPERIMENTS

Introduction

The relevance of the WFE strategy, as proposed in the last section, is discussed there for computing the forced response of different kinds of structures, namely:

1. Three connected shells, with a punctual force, subject to lineic densities of translational and rotational springs (see Figure 9). 2. Two sets of two connected shells, with translational springs subject to vertical displacements, connected to one coupling element having a non-uniform curvature (see Figure 12). This kind of structure may be understood as coarsely representing a part of a train structure. 3. Six shells, with punctual forces, connected to one cylindrical structure having a conical head (see Figure 17). This kind of structure may be understood as coarsely representing an aircraft structure.

The first test case (Figure 9) is quite simple. Its investigation aims at validating the strategy depicted in Section 3 which describes the wave-based coupling conditions between a set of noncoplanar shells that exhibit lineic impedances and punctual forces over their coupling interfaces. The other test cases (Figures 12 and17) appear more complex as they involve coupling elements whose shapes are not as simple as rectangular flat shells. Their investigation aims at validating the strategy depicted in Section 4 which mixes wave-based coupling conditions, between connected shells and between shells and coupling elements, in a whole matrix formulation.

For each test case, the frequency response function (say, the quadratic velocity ) is computed at some measurement point using the WFE strategy described in Section 4. In this framework, the flat shells are meshed periodically (see Section 2.1) while incompatible meshes can be considered over the coupling interfaces (see Sections 3 and 2.2). The criterion for meshing these shells is to discretize the wavelength of the bending wave traveling in an infinite equivalent Reissner-Mindlin plate at least by means of 8 elements at the maximum frequency considered (see Section 2.1).

Within the WFE framework, a matrix formulation of the form AQ = F is considered (see Section 4.3) which is solved using MATLAB R . As a preliminary step, the mass/stiffness matrices of several FE substructures and coupling elements are to be obtained. This is done by means of COMSOL Multiphysics R . The mass/stiffness matrices of the substructures are used to compute the wave modes (see Section 2.1), while the mass/stiffness matrices of the coupling elements are used to compute the matrices {C c ij } mentioned in Section 2.2 and Eq. (34). For each test case, the WFE solution is compared with a reference FE solution provided by COMSOL Multiphysics R when the whole structure is finely meshed using 2D shell elements of arbitrary shapes. Considering the reference solution, the FE meshes are supposed to be compatible across the coupling interfaces. Regarding test cases 2 and 3, a Component Mode Synthesis / Craig-Bampton (CMS/CB) procedure is investigated within the WFE framework which uses a few fixed interface modes of the coupling elements to compute the matrices {C c ij } [START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF]. The relevance of this CMS-based WFE procedure, in terms of accuracy and CPU time saving, is discussed compared to the FE method. Additional experiments are carried out regarding the test case 2 to compare the WFE method with the conventional CMS/CB method when all the shells, as well as the coupling element, are modeled in terms of fixed interface modes.

Test case 1: three connected shells with lineic densities of translational and rotational springs, and a punctual force

The frequency response of a structure whose FE model is depicted in Figure 9 is investigated using the WFE method over a frequency band [10 Hz, 300 Hz]. The global structure is composed

Here, the quadratic velocity is defined as the norm of the velocity vector.

of three rectangular shells of same width 1 m (z-direction) while their respective lengths are 0.9 m, √ 2 × 0.6 m and 0.7 m. Both shells 1 and 3 are oriented along the horizontal direction, while shell 2 is rotated from an angle α = 45 o around the z-direction. Both shells 1 and 3 are made of aluminum material whose characteristics are given in Section 2.1 (it is worth recalling that their thickness is 5 × 10 -3 m). In contrast, shell 2 is made of steel material (density ρ = 7800 kg.m -3 , Young's modulus E = 210 × 10 9 Pa, Poisson's ratio ν = 0.3, loss factor η = 0.01) with a thickness 10 × 10 -3 m. Shell 1 is subject to a transverse punctual force (vector F 0 ) of magnitude 10 N at location (x = 0.4 m, z = 0.3 m). Shells 1 and 2 are coupled over a common interface with a uniform lineic density of translational springs k t = 10 7 N.m -2 acting in the y-direction. Also, shells 2 and 3 are coupled over a common interface which is supposed to be free from external excitation sources and mechanical impedances. Apart from the coupling conditions, the left edge of shell 1 is subject to a uniform lineic density of rotational springs k r = 5 × 10 3 N acting in the z-direction, while the right edge of shell 3 is clamped. The reference FE model of the structure, as provided by COMSOL Multiphysics R using 2D triangular shell elements of arbitrary shapes, is depicted in Figure 9. Regarding this reference FE model, the shell meshes are compatible across the coupling interfaces; also, the maximum element size is 0.04 m, which leads to 17, 304 DOFs. This yields the wavelengths to be described at least by means of 10 elements (14 elements for shell 2) at the maximum frequency considered (i.e., 300 Hz), which is completely satisfactory from the point of view of numerical dispersion. Within the WFE framework, the shells are meshed periodically, as depicted in Figure 10. Also, the first shell is to be split into two sub-shells 1 and 2 of respective lengths 0.4 m and 0.5 m, whose common interface Γ 12 involves the punctual force (F 0 ). The need to split this shell appears as a requirement of the WFE procedure since excitation sources can only be considered on the left or right boundaries of waveguides. In this framework, the other shells are numbered as shell 10). The FE discretization of each shell is chosen so that the rule of 8 elements per wavelength is satisfied at 300 Hz (see Section 2.1). The FE meshes turn out to be incompatible across Γ 23 and Γ 34 , i.e., between shells 2 and 3 and between shells 3 and 4. The need to consider different meshes is explained since shell 3 exhibits a bending rigidity which is high compared to other shells, meaning that it can be discretized using a coarser mesh with a view to reducing CPU times. To address these incompatible meshes, two Lagrange Multipliers fields are respectively introduced on Γ 23 and Γ 34 . These Lagrange Multipliers fields are discretized with respect to the sides with the finest meshes (see Section 3.1), using the same interpolation functions as the corresponding shells. In the present case, these sides appear linked to shells 2 and 4. Also, the lineic densities of springs (interface Γ 23 and left edge of shell 1) are discretized in the same way as shell 2 and shell 1, respectively. As a result, the impedance matrices turn out to be of the form k N T q Γ N q Γ dS (N q Γ being the matrix of interpolation functions used to discretize the displacement fields of the shells over their left or right boundary). Over the shell boundaries, the displacement fields are supposed to be discretized by means of linear interpolation functions, which in the present case are chosen as Lagrange polynomials.

           I n1 -C ⋆ 11 µ N1 1 0 0 • • • 0 0 -C 11 µ N1 1 I n1 0 -C 12 µ N2 2 • • • 0 0 -C ⋆ 21 µ N1 1 0 I n2 -C ⋆ 22 µ N2 2 • • • 0 0 0 0 -C 22 µ N2 2 I n2 • • • 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 • • • I n4 -C ⋆ 44 µ N4 4 0 0 0 0 • • • -C 44 µ N4 4 I n4                      Q ref⋆ 1 Q ref 1 Q ref⋆ 2 Q ref 2 . . . Q ref⋆ 4 Q ref 4           =           0 F 1 F ⋆ 2 0 . . . 0 0           , (41)
where the expressions of matrices C ij , C ⋆ ij (i, j = 1, . . . , 4) and vectors F 1 and F ⋆ 2 follow from Section 3 and Appendix B. Also, the vectors F 1 and F ⋆ 2 relate the punctual force (F 0 ) on Γ 12 . Solving this matrix formulation while invoking the wave mode expansion (2) provides the displacements/rotations, as well as the internal forces/moments, at any location within the shells. Particularly, the quadratic velocity of shell 1, at the location (x = 0.3 m, z = 0.3 m), is computed over the frequency band [10 Hz, 300 Hz] using 581 discrete frequencies with identical frequency steps 0.5 Hz. The WFE solution is compared with the reference FE solution when the full FE model depicted in Figure 9 is solved. The results are shown in Figure 11(a). Also, the relative errors between the FE and WFE solutions are computed. Here, both the norms of the real and imaginary parts of the velocity vector, averaged over 15 frequency bands of same width * * , are investigated (see Figure 11(b)). In other words, the magnitude as well as the phase of the velocity vector are highlighted. It is shown that the WFE solution perfectly agrees with the reference FE solution, even at high frequencies when the frequency response function is rather complex. Here, the relative errors are less than 10% over the whole frequency band. This particularly means that the WFE-based matrix formulation (41) can be used without difficulty even when lineic impedances are dealt with, without introducing numerical ill-conditioned problems (e.g., spurious resonances, among others).

CPU times and numerical conditioning

The accuracy of the WFE method has been emphasized above compared to the classic FE method. Its efficiency compared the FE method lies in the fact that the CPU times are considerably reduced.

In fact, the WFE method requires a matrix A of size 2 i n i × 2 i n i (= 684 × 684 in the present case) to be inverted at each discrete frequency considered (i.e., 581 frequencies). In contrast, the FE method requires a matrix of size 17, 304 × 17, 304 to be inverted at the same frequencies. Here, it takes almost 980 s to compute the forced response of the structure with the WFE method (with MATLAB R ) against 6, 560 s with FE (with COMSOL Multiphysics R ), considering a FE model with an element size 0.05 m -i.e., which is comparable with the WFE model regarding shells 1, 2 and 4 -using an Intel R Core TM 2 Duo processor. It must be noticed that even if the CPU times required to calculate the wave modes are taken into account -i.e., 890 s to compute the wave modes of shells 1, 2 and 4 (here, these shells involve identical substructures, meaning that they exhibit the same wave modes) and 330 s to compute the wave modes of shell 3 -, the global computational cost of the WFE method remains relatively small, i.e., 66% less than FE. It must also be noticed that the wave modes are to be computed once and for all, regardless of the kind of boundary and coupling conditions considered over the left or right edges of the shells. This interesting feature of the WFE method enables the forced response of the global structure to be computed many times with very small CPU times (i.e., without the need to re-compute the wave modes) to investigate different kinds of boundary and coupling conditions (parametric analysis).

Otherwise, the well-conditioning of WFE matrix formulations has been emphasized in Section 4.4. The WFE solution can be claimed to be less sensitive to small rounding errors of input data provided that the condition number of A (see Eq. ( 41)) -κ(A) = ||A|| 2 ||A -1 || 2 -is small enough. In the present case, κ(A) < 10 6 over the whole frequency band which is completely satisfactory from the point of view of m-digit floating point arithmetics, i.e., the relative error for expressing the structure behavior is expressed as 10 α-m , where 10 -m is the machine precision (m = 32 for many softwares) while κ(A) ≤ 10 α (α < 6 in the present case). It is worth emphasizing that the fact that κ(A) is small, although due to the particular form of the matrix A (see Section 4.4), is also explained because the structure is damped: considering no damping will have as effect to make A singular at the resonance frequencies of the structure, leading to infinite responses (which is of course unrealistic).

Test case 2: two sets of two connected shells, with translational springs subject to vertical displacements, connected to one coupling element having a non-uniform curvature

The forced response of the structure depicted in Figure 12 is investigated using the WFE method over a frequency band [1 Hz, 100 Hz]. Such a structure might be viewed as coarsely representing a part of a train structure. It is composed of four rectangular flat shells of same width 2 m (z-direction), i.e., two vertically oriented shells of length 2.6 m connected to two 45 o -oriented shells of length 0.5 m. The vertically oriented flat shells are distant from 2.4 m (x-direction) and connected to one coupling element which represents a curved shell with a non-uniform curvature † † . The four rectangular flat shells exhibit a same thickness of 7 × 10 -3 m and are made of steel material (see Section 5.2 for the material characteristics). Also, the coupling element exhibits a thickness of 4 × 10 -3 m and is made in aluminum (see Section 2.1 for the material characteristics). Each bottom corner of the 45 o -oriented flat shells is connected to three translational springs k t = 10 7 N.m -1 in the three x-, yand z-directions. The bottom ends of the vertical springs (y-direction) are subject to vertical displacements (vector ±q 0 of magnitude 10 -3 m) acting in opposite directions (see Figure 12). The reference FE model of the structure is depicted in Figure 12. Here, a FE mesh involving a maximum element size of 0.08 m is considered which is continuous across the coupling interfaces, leading to 29, 670 DOFs. Such a FE mesh enables the wavelengths to be discretized at least by 8 elements (see Section 2.1).

Within the WFE framework, the flat shells -namely, shells 1, 2, 3, 4 -are meshed periodically, as depicted in Figure 13. Each flat shell is described by means of identical FE substructures of length 0.1 m whose left/right boundaries involve 21 nodes (see Figure 13). As a matter of rule, these nodes are uniformly spread on the substructure boundaries. A number of N 2 = N 3 = 26 substructures are used to describe shells 2 and 3, while N 1 = N 4 = 5 substructures are used to describe shells 1 and 4. The maximum element size involved by these FE substructures is 0.1 m which appears coherent regarding the rule of thumb of 8 elements (at least) per wavelength (see above). Here, mesh compatibility is assumed between the shells. On the other hand, the coupling element is meshed using shell elements of arbitrary shapes with a maximum element size of 0.08 m, disregarding the FE mesh used to describe the connected shells. Thus the FE meshes between shells 2, 3 and the coupling element turn out to be incompatible. Such an issue is addressed by means of the strategy depicted in Section 2.2 (see also ref. [START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF]). Within the WFE framework, the harmonic forced response of the structure is computed by solving a wave-based matrix formulation of the form AQ = F (see Section 4.3) which is based on Eqs. (32-36). In the present case, this matrix formulation is expressed as

           I n1 -C ⋆ 11 µ N1 1 0 0 0 0 • • • -C 11 µ N1 1 I n1 0 -C 12 µ N2 2 0 0 • • • -C ⋆ 21 µ N1 1 0 I n2 -C ⋆ 22 µ N2 2 0 0 • • • 0 0 -C c 22 µ N2 2 I n2 -C c 23 µ N3 3 0 • • • 0 0 0 0 I n3 -C ⋆ 33 µ N3 3 • • • 0 0 -C c 32 µ N2 2 0 -C c 33 µ N3 3 I n3 • • • . . . . . . . . . . . . . . . . . . . . .                      Q ref⋆ 1 Q ref 1 Q ref⋆ 2 Q ref 2 Q ref⋆ 3 Q ref 3 . . .           =           F ⋆ 1 0 0 0 0 0 . . .           , (42)
where the expressions of {C ij }, {C ⋆ ij }, {C c ij } follow from Section 3, Section 2.2 and Appendix B. Otherwise, the vectors F ⋆ 1 and F ⋆ 4 (not mentioned in Eq. ( 42)) relate the consideration of imposed displacements (vector ±q 0 ) at the bottom corners of shells 1 and 4. Solving the matrix formulation AQ = F provides the amplitudes of the wave modes traveling along the shells, and further on, the displacements/rotations and internal forces/moments. Here, 126 incident/reflected wave modes are considered to describe the behavior of each shell. The quadratic velocity of shell 2 (measurement point: medium of the width, 1 m above the bottom edge) is computed over the frequency band [1 Hz, 100 Hz] considering 397 discrete frequencies with identical frequency steps. To reduce further the CPU times associated with the computation of {C c ij }, a Component Mode Synthesis / Craig-Bampton (CMS/CB) procedure is used as detailed in ref. [START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF]. In this framework, the dynamic stiffness matrix of the coupling element is described by means of a small number of fixed interface modes. Static modes are also invoked (see ref. [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF]). In the present case, 40 fixed interface modes are used. The resulting CMS/CB based WFE solution is compared with the reference FE solution, when the full FE model depicted in Figure 12 is solved. The results are shown in Figure 14(a). Also, the relative errors between the FE and WFE solutions, averaged over 15 frequency bands of same width (see last test case), are computed. Again, both the norms of the real and imaginary parts of the velocity vector are highlighted (see Figure 14(b)). It is shown that the WFE solution perfectly agrees with the reference FE solution. In particular, the relative errors are less than 10% over the whole frequency band. Again, this highlights the efficiency of the WFE strategy. It can be used without difficulty even when coupling elements, whose characteristics can be strongly different from those of the shells, are dealt with. Regarding the CPU times, it takes almost 12, 650 s with FE to compute the forced response of the structure. On the other hand, the WFE procedure requires: (i) the computation of wave modes (660 s); (ii) the computation of the matrices {C c ij } based on 40 fixed interface modes of the coupling element (180 s); (iii) the computation of the matrix formulation AQ = F (800 s). As a whole, this leads to 1640 s, i.e., 87% CPU time saving compared to the FE method. This highlights the relevance of the WFE method.

Comparisons with the Craig-Bampton Method

A comparison between the WFE method and the CMS/CB method [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] is proposed to predict the dynamic behavior of the whole structure. Within the CMS/CB framework, it is proposed to describe each shell, as well as the coupling element, by means of static modes and fixed interface modes.

Here, the shells and the coupling element are meshed in the same way as the WFE modeling (cf. Figure 13). Also, 2 × 126 static modes and 126 fixed interface modes are used to describe the dynamic behavior of each shell. In other words, the number of fixed interface modes is chosen so that it matches the number of incident/reflected wave modes used by the WFE method. Also, 40 fixed interface modes are used for the coupling element, i.e., in the same way as the WFE method. Using the CMS/CB method yields the quadratic velocity (the measurement point has been specified above) to be calculated as shown in Figure 15(a). Also, the relative error (between CMS/CB and FE solutions) is shown in Figure 15(b) regarding the real part of the velocity vector; the result is compared with the relative error involved by the WFE method. Although the frequency response function seems in good agreement with the reference solution, further insights of the relative error reveal a drawback of the CMS/CB method for predicting the vibratory levels. This is particularly true around the two main resonances of the structures (over [40 Hz, 60 Hz]) when the error goes up to 20%, as opposed to the WFE method where the error remains below 10%. Also, the CMS/CB solution appears less accurate than WFE at high frequencies. Otherwise, the condition numbers ‡ ‡ of both CMS/CB and WFE matrix formulations have been plotted in Figure 16 to highlight further on the accuracy of the WFE method. It is shown that the condition number of the WFE matrix formulation is almost 500 times smaller compared to the CMS/CB method when averaged over the whole frequency band. This highlights the efficiency of the WFE method, in terms of numerical conditioning, compared to the CMS/CB method. The feature of the WFE matrix formulation follows from the particular form of the matrix A which has been explained in depth in Section 4.3. It could be emphasized that the lack of accuracy of the CMS/CB method can be solved in theory by increasing the number of fixed interface modes for the shells, which could however has the consequence of increasing the condition number of the related matrix formulation (this is understood since the size of the matrix formulation is enlarged).

Regarding CPU times, it takes around 500 s with the CMS/CB method to compute the forced response of the structure, against 1640 s with the WFE method (see above). Those CPU times are reached using an Intel R Core TM 2 Duo processor, and MATLAB R . Within the CMS/CB framework, the static modes and fixed interface modes are computed twice regarding shells 1 and 2 (indeed, the fixes interface and static modes of shells 3 and 4 can be simply deduced from those of shells 1 and 2). In contrast, the WFE method makes use of the same wave basis for all the shells. The fact that the same wave basis can be used for all the shells is explained since the wave modes do not depend on the lengths and boundary conditions of waveguides. This interesting feature demonstrates the capability of the WFE method to involve CPU times that can be reduced further, regarding structures involving many shells.

Test case 3: six shells, with punctual forces, connected to one cylinder having a conical head

The forced response of the structure depicted in Figure 17 Within the WFE framework, the flat shells are meshed periodically, as depicted in Figure 18 considering the half of the structure (see above). Here, the shell involving the punctual force is split into two sub-shells -namely, shells 1 and 2. The need to split this shell into two sub-shells has been explained in Section 5.2. Each shell is described by means of identical FE substructures of length 0.05 m: the substructures of shells 1 and 2 involve 11 nodes over their left or right boundaries, compared to 7 nodes regarding shells 3 and 4 (see Figure 18). Those nodes are supposed to be uniformly spread on the substructure boundaries. The maximum element size involved by these FE substructures is 0.05 m which appears similar to the FE model depicted in Figure 17. Also, the coupling element is meshed using 2D shell elements of arbitrary shapes with a maximum element size of 0.05 m (say, in a similar way as in Figure 17). The connection between this arbitrary mesh and the periodic meshes used to discretize the shells is carried out by means of the strategy depicted in Section 2.2 (see also ref. [START_REF] Mencik | Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies[END_REF]). Otherwise, mesh compatibility is assumed between the shells. Within the WFE framework, the harmonic forced response of the structure is computed by solving a wave-based matrix formulation of the form AQ = F (see Section 4.3) which is based on Eqs. (32-36). In the present case, this matrix formulation is expressed as

           . . . . . . . . . . . . . . . . . . . . . • • • I n2 -C ⋆ 22 µ N2 2 0 0 0 0 • • • -C c 22 µ N2 2 I n2 -C c 23 µ N3 3 0 -C c 24 µ N4 4 0 • • • 0 0 I n3 -C ⋆ 33 µ N3 3 0 0 • • • -C c 32 µ N2 2 0 -C c 33 µ N3 3 I n3 -C c 34 µ N4 4 0 • • • 0 0 0 0 I n4 -C ⋆ 44 µ N4 4 • • • -C c 42 µ N2 2 0 -C c 43 µ N3 3 0 -C c 44 µ N4 4 I n4                      . . . Q ref⋆ 2 Q ref 2 Q ref⋆ 3 Q ref 3 Q ref⋆ 4 Q ref 4           =           . . . F ⋆ 2 0 0 0 0 0           , (43)
where the matrices {C ij }, {C ⋆ ij }, {C c ij } are expressed by considering Section 3, Section 2.2 and Appendix B. The vector F ⋆ 2 relates the punctual force (vector F 0 ) at the interface between shells 1 and 2. Here, 66 incident/reflected wave modes for shells 1, 2 and 42 incident/reflected wave modes for shells 3, 4 are considered to describe the dynamic behavior of the whole structure. Also, the dynamic stiffness matrix of the coupling element is modeled using the CMS/CB procedure described previously (last test case), by means of 20 fixed interface modes. The resulting WFE solution is computed over the frequency band [10 Hz, 300 Hz] using 581 discrete frequencies with identical frequency steps. The WFE solution is compared with the reference FE solution when the full FE model (half of the structure) depicted in Figure 17 is solved. The results are shown in Figure 19(a) regarding the frequency response function. Also, the relative errors between FE and WFE solutions, averaged over 15 frequency bands of same width composing the whole frequency band, are shown (see Figure 19(b)). Again, both the real and imaginary parts of the velocity vector are considered. It is shown that the WFE solution perfectly agrees with the reference FE solution. In particular, the relative errors appears less than 8% over the whole frequency band. The condition number of the WFE matrix formulation -i.e., κ(A) -is almost 10 4 when averaged over the whole frequency band. Again, this highlights the well-conditioning of WFE matrix formulations.

Regarding CPU times, it takes almost 13, 000 s with FE to compute the forced response of the structure, against 380 s with the WFE method (this includes the computation of wave modes, the computation of the matrices {C c ij } based on 20 fixed interface modes of the coupling element, the computation of the matrix formulation AQ = F). This leads to 97% CPU times savings compared to the FE method. Again, this gives credit to the WFE method as an efficient means to describe the forced response of structures that can be complex, even at high frequencies.

CONCLUSION

In this paper, a wave finite element based strategy has been proposed for computing the forced response of structures involving several connected rectangular flat shells. Within the WFE framework, the shells are connected along their left or right boundaries; also, each shell is meshed periodically by means of identical FE substructures. The shell dynamic behavior is described by means of numerical wave modes traveling towards and away from the coupling interfaces. In the present study, those rectangular flat shells are enabled to be oriented in different ways as well as they can be connected along coupling interfaces where mesh compatibility is not necessarily assumed. Also, the coupling interfaces are supposed to involve several kinds of external excitation sources (i.e., punctual and lineic forces/moments) as well as mechanical impedances (e.g., lineic densities of springs). Finally, the shells are supposed to be connected to one or several coupling elements whose shapes, as well as dynamic behavior, can be relatively complex.

A mesh tying formulation based on Lagrange Multipliers, which describes the coupling conditions between two connected shells or between several shells connected to one coupling element, has been proposed and adapted to the WFE framework. The resulting wave-based coupling conditions have been formulated. A WFE strategy has been proposed to assess the forced response of structures involving several connected shells and coupling elements. A CMS/CB procedure has been used to describe the dynamic behavior of the coupling elements by means of a small number of fixed interface modes. The WFE strategy is based on a global wave-based matrix formulation for computing the amplitudes of the wave modes traveling along the shells. The determination of the displacements/rotations and forces/moments within the shells follows from a wave mode expansion procedure. One feature of the proposed wave-based matrix formulation is that it is well conditioned, which particularly means that it can be used without difficulty even when a large number of shells are dealt with. As a second feature, it enables the CPU times to be considerably reduced compared to the conventional FE method. The relevance of this WFE strategy has been clearly highlighted considering the forced response of three kinds of structures involving several connected flat shells and coupling elements of variable complexities. Also, its efficiency in terms of accuracy has been highlighted in comparison with the conventional CMS/CB method.

APPENDIX A. DERIVATION OF MATRICES A ref AND A inc

Inserting Eqs. ( 17) and [START_REF] Allman | Implementation of a flat facet shell finite element for applications in structural dynamics[END_REF] in Eqs. ( 14) and [START_REF] Ichchou | Wave finite elements for low and mid-frequency description of coupled structures with damage[END_REF] gives

(Φ ref q ) i Q ref i -L i BL T i+1 (Φ ref⋆ q ) i+1 Q ref⋆ i+1 (A-1) = -(Φ inc q ) i Q inc i -L i BL T i+1 (Φ inc⋆ q ) i+1 Q inc⋆ i+1 ,
and

-(Φ ref⋆ F ) i+1 Q ref⋆ i+1 + L i BL T i+1 T (Φ ref F ) i Q ref i + L i iωZL T i (Φ ref q ) i Q ref i = --(Φ inc⋆ F ) i+1 Q inc⋆ i+1 + L i BL T i+1 T (Φ inc F ) i Q inc i + L i iωZL T i (Φ inc q ) i Q inc i + BL T i+1 T F ex . (A-2)
Left multiplying Eq. (A-1) by (Φ ref q ) -1 i leads to

Q ref i -(Φ ref q ) -1 i L i BL T i+1 (Φ ref⋆ q ) i+1 Q ref⋆ i+1 = -(Φ ref q ) -1 i (Φ inc q ) i Q inc i -(Φ ref q ) -1 i L i BL T i+1 (Φ inc⋆ q ) i+1 Q inc⋆ i+1 , (A-3)
while left multiplying Eq. (A-2) by -(Φ ref⋆ F ) -1 i+1 leads to • Case when two vectors of forces/moments F ex = F 0 and F ex = F ′ 0 are respectively considered on the left and right boundaries. Expressing these vectors in the local frame R i of the shell, i.e., (F (1) i ) Ri = L i F 0 and (F (Ni+1) i

Q ref⋆ i+1 -(Φ ref⋆ F ) -1 i+1 L i BL T i+1 T (Φ ref F ) i Q ref i + L i iωZL T i (Φ ref q ) i Q ref i = -(Φ ref⋆ F ) -1 i+1 (Φ inc⋆ F ) i+1 Q inc⋆ i+1 -(Φ ref⋆ F ) -1 i+1 L i BL T i+1 T × (Φ inc F ) i Q inc i + L i iωZL T i (Φ inc q ) i Q inc i -(Φ ref⋆
) Ri = L i F ′ 0 , while invoking the wave mode expansion (2) leads to

(Φ inc⋆ F ) i Q inc⋆ i + (Φ ref⋆ F ) i Q ref⋆ i = -L i F 0 , (B-1) (Φ inc F ) i Q inc i + (Φ ref F ) i Q ref i = L i F ′ 0 , (B-2)
where the minus sign ahead of L i F 0 is introduced since the left boundary of the shell is of concern (see comment below Eq. ( 1)). As a result, Eqs. (B-1) and (B-2) lead to

Q ref⋆ i = C ⋆ ii Q inc⋆ i + F ⋆ i (left boundary) and Q ref i = C ii Q inc i + F i (right boundary),
where

C ⋆ ii = -(Φ ref⋆ F ) -1 i (Φ inc⋆ F ) i , C ii = -(Φ ref F ) -1 i (Φ inc F ) i , F ⋆ i = -(Φ ref⋆ F ) -1 i L i F 0 , F i = (Φ ref F ) -1 i L i F ′ 0 .
(B-3)

• Case when two vectors of displacements/rotations q 0 and q ′ 0 are respectively considered on the left and right boundaries. Expressing these vectors in the local frame R i , i.e., (q (1) i ) Ri = L i q 0 and (q (Ni+1) i ) Ri = L i q ′ 0 , while invoking the wave mode expansion (2) leads to

(Φ inc⋆ q ) i Q inc⋆ i + (Φ ref⋆ q ) i Q ref⋆ i = L i q 0 ,
(B-4)

(Φ inc q ) i Q inc i + (Φ ref q ) i Q ref i = L i q ′ 0 . (B-5)
Again, this yields

Q ref⋆ i = C ⋆ ii Q inc⋆ i + F ⋆ i (left boundary) and Q ref i = C ii Q inc i
+ F i (right boundary), where

C ⋆ ii = -(Φ ref⋆ q ) -1 i (Φ inc⋆ q ) i , C ii = -(Φ ref q ) -1 i (Φ inc q ) i , F ⋆ i = (Φ ref⋆ q
) -1 i L i q 0 , F i = (Φ ref q ) -1 i L i q ′ 0 .

(B-6)

• Case when two matrices of mechanical impedance Z = Z 0 and Z = Z ′ 0 and two vectors of forces/moments F ex = F 0 and F ex = F ′ 0 are respectively considered on the left and right boundaries. In that case, the boundary conditions are deduced from Eq. ( 11) with F Γ i+1 = 0. According to Eq. ( 12), this gives

-(F i ) Ri + L i F ex -L i iωZL T i (q i ) Ri = 0. (B-7)
Invoking the wave mode expansion (2) leads to

(Φ inc⋆ F ) i Q inc⋆ i + (Φ ref⋆ F ) i Q ref⋆ i = (L i iωZ 0 L T i ) (Φ inc⋆ q ) i Q inc⋆ i + (Φ ref⋆ q ) i Q ref⋆ i -L i F 0 , (B-8) (Φ inc F ) i Q inc i + (Φ ref F ) i Q ref i = -(L i iωZ ′ 0 L T i ) (Φ inc q ) i Q inc i + (Φ ref q ) i Q ref i + L i F ′ 0 . (B-9)
Again, this yields

Q ref⋆ i = C ⋆ ii Q inc⋆ i + F ⋆ i (left boundary) and Q ref i = C ii Q inc i
+ F i (right boundary), where

C ⋆ ii = -(Φ ref⋆ F ) i -(L i iωZ 0 L T i )(Φ ref⋆ q ) i -1 (Φ inc⋆ F ) i -(L i iωZ 0 L T i )(Φ inc⋆ q ) i , C ii = -(Φ ref F ) i + (L i iωZ ′ 0 L T i )(Φ ref q ) i -1 (Φ inc F ) i + (L i iωZ ′ 0 L T i )(Φ inc q ) i , F ⋆ i = -(Φ ref⋆ F ) i -(L i iωZ 0 L T i )(Φ ref⋆ q ) i -1 L i F 0 , F i = (Φ ref F ) i + (L i iωZ ′ 0 L T i )(Φ ref q ) i -1 L i F ′ 0 .
(B-10)

Figure 2 .

 2 Figure 2. Finite element model of a rectangular flat shell and representative substructure of length d (the sought waves are those traveling in the x-direction).

Figure 3 .

 3 Figure 3. Several dispersion curves for the wave modes traveling along the rectangular flat shell depicted in Figure 2: (--) wave modes whose shapes are depicted in Figures 4 and 5.

Figure 4 .

 4 Figure 4. Several wave shapes (real parts of the displacement components) at 10 Hz.

Figure 6 .

 6 Figure 6. Several non-coplanar flat shells connected to one coupling element having an arbitrary shape.

Figure 7 .

 7 Figure 7. Two coupled shells i and i + 1 having different FE meshes and different orientations α i and α i+1 (around the z-direction), connected along one common interface Γ involving excitation sources (e.g., punctual force F 0 ) and mechanical impedances.

Figure 8 .

 8 Figure 8. WFE-based description of the coupled shells depicted in Figure 7.

  are the vectors of wave amplitudes expressed at the right boundary of the shell (i.e., considering the substructure boundary N i + 1), while {Q inc⋆ i , Q ref⋆ i } are the vectors of wave amplitudes expressed at the left boundary (i.e., considering the substructure boundary 1).

Figure 9 .

 9 Figure 9. FE model of a structure involving three non-coplanar flat shells with lineic densities of translational and rotational springs (k t and kr), a punctual force (F 0 ) and a clamped edge.

3 and shell 4 .

 4 Shells 1, 2 and 4 are discretized by means of identical substructures (see Figure 10) of length d = 0.05 m and containing 21 nodes (i.e., 126 DOFs) uniformly spread on each left or right boundary. In this case, shells 1, 2 and 4 are respectively composed of N 1 = 8, N 2 = 10 and N 4 = 14 substructures. In contrast, shell 3 is composed of N 3 = 12 substructures of length d 3 = √ 2 × 0.05 m, which are discretized by means of 15 nodes (i.e., 90 DOFs) over each left or right boundary (see Figure

Figure 10 .

 10 Figure 10. WFE-based description of the structure depicted in Figure 9, involving four shells.

Figure 11 .

 11 Figure 11. (a) Frequency response function of the structure depicted in Figure 9: (--) FE reference solution; (• • •) WFE solution. (b) Relative errors of the velocity vector: (-•-) real part; (-o --) imaginary part.

Figure 12 .

 12 Figure 12. FE model of a structure involving two sets of two non-coplanar flat shells, with translational springs subject to vertical displacements (±q 0 ), connected to one coupling element having a non-uniform curvature.

Figure 13 .

 13 Figure 13. WFE-based description of the structure depicted in Figure 12, involving four flat shells (for each shell, the direction of wave propagation is mentioned by means of oscillating arrows).

Figure 14 .

 14 Figure 14. (a) Frequency response function of the structure depicted in Figure 12: (--) FE reference solution; (• • •) WFE solution. (b) Relative errors of the velocity vector: (-•-) real part; (-o --) imaginary part.

Figure 15 .

 15 Figure 15. (a) Frequency response function of the structure depicted in Figure 12: (--) FE reference solution; (• • •) CMS/CB solution (40 fixed interface modes for the coupling element; 126 fixed interface modes for each shell). (b) Relative errors for the real part of the velocity vector: (-•-) CMS/CB; (-•-) WFE.

Figure 16 .

 16 Figure 16. Condition number of the matrix formulation: (--) CMS/CB; (--) WFE.

  is investigated using the WFE method over a frequency band [10 Hz, 300 Hz]. Such a structure might be viewed as coarsely representing a part of an aircraft structure. It is composed of six rectangular flat shells involving two horizontally oriented shells of length 3.5 m and width 0.5 m, and four 45 o -oriented shells of length 0.5 m and width 0.3 m. These flat shells are connected to one coupling element representing a cylinder of length 3 m and radius 0.25 m whose ends are respectively connected to one conical head of length 0.5 m and a disk-shaped flat cap. The shells and the coupling element exhibit the same thickness 5 × 10 -3 m and are made of steel material (see Section 5.2 for the material characteristics). Each horizontally oriented flat shell is excited by a vertical punctual force (vector F 0 ) of magnitude 10 N acting at a location distant of 1.5 m from the coupling element, at the middle of the shell width. The FE model of the structure is depicted in Figure 17. It involves 2D triangular shell elements of arbitrary shapes for the flat shells as well as the coupling element. Here, a FE mesh involving a maximum element size of 0.05 m for the flat shells and the coupling element is considered which is continuous across the coupling interfaces. The global problem exhibits a plane symmetry (plane (y, z)) with regard to the main axis of the coupling element, meaning that only half of the structure can be studied. The number of DOFs involved by the related FE model is 21, 726. Such a FE model is considered for computing the forced response of the structure taken as a reference FE solution.

Figure 17 .

 17 Figure 17. FE model of a structure involving six non-coplanar flat shells, with punctual forces (F 0 ), coupled with one cylinder having a conical head.

Figure 18 .

 18 Figure 18. WFE-based description of the structure depicted in Figure 17, involving four flat shells (for each shell, the direction of wave propagation is mentioned by means of oscillating arrows).

Figure 19 .

 19 Figure 19. (a) Frequency response function of the structure depicted in Figure 17: (--) FE reference solution; (• • •) WFE solution. (b) Relative errors of the velocity vector: (-•-) real part; (-o --) imaginary part.

F ) -1 i+1 BL T i+1 T

 i+1 F ex . (A-4) Expressing, in matrix form, Eqs. (A-3) and (A-4) leads to Eq. (19) where A ref and A inc are given by Eqs. (20) and (21). APPENDIX B. EXPRESSIONS OF MATRICES C ⋆ ii , C ii , AND VECTORS F ⋆ i AND F i OVER THE LEFT/RIGHT BOUNDARIES WHERE COUPLING CONDITIONS DO NOT APPLY The wave-based boundary conditions on the left and right boundaries of any shell i are expressed as follows:

  normalized with respect to its euclidean norm (Φ H qr Φ qr + Φ H Fr Φ Fr ) 1/2 to improve further on the conditioning of A ref . According to Eq. (20), A ref is formulated in such a way it involves identity matrices as diagonal block terms. It is invertible provided that det(I ni+1 -A ref 21 A ref 12 ) = 0 (A ref 21 and A ref 12 being the bottom left and top right off-diagonal block terms of the matrix). Such an assumption appears to be satisfied provided that the matrix Im{A ref 21 A ref 12 } is full rank (where Im denotes the imaginary part). This question can be viewed as proving that A ref 21 and A ref 12 are full rank (it is assumed that these conditions are sufficient to prove that Im{A ref 21 A ref 12 } is full rank, taking into account that the matrices are complex). Regarding for instance A ref

[START_REF] Mead | The forced vibration of one-dimensional multi-coupled periodic structures: An application to finite element analysis[END_REF]

, this requires one to verify that (i) a non-zero matrix (Φ ref⋆ q

  • (O 1 ) L (resp. (O 1 ) R ) the set of integers which correspond to shells that are connected to other shells over their left (resp. right) boundaries, • O 2 the set of integers which correspond to shells that are connected to coupling elements, • (O 3 ) L (resp. (O 3 ) R ) the set of integers which correspond to shells whose left (resp. right) boundaries do not involve coupling conditions.Regarding the coupling elements, the set O 2 can be partitioned as O 2 = ∪ s O s 2

* * The objective behind this averaging process is to lower the influence of slight shifts among the resonance frequencies (between the WFE and FE solutions) which are sources of high errors although they are of minor importance

† † The curvature of the coupling element is actually based on a cubic Bézier function.