F Treyssède
email: fabien.treyssede@univ-eiffel.fr

A model reduction method for fast finite element analysis of continuously symmetric waveguides

Keywords: waveguide, reduction, finite element, mode, veering, leaky

A model reduction method for fast finite element analysis of continuously symmetric waveguides

Introduction

The knowledge of modal properties in elastic waveguides is of great interest for the dynamic analysis of thin and elongated structures (plates, bars, pipes, cables...). It provides the theoretical background for nondestructive evaluation, structural health monitoring, vibration and noise reduction, design of metamaterials, statistical energy analysis... Due to the multimodal and dispersive nature of guided waves, modeling tools are necessary to determine dispersion curves. These curves provide the properties of each wave mode, for instance given by phase velocity, group velocity or attenuation, as functions of frequency. In practice, accurate knowledge of dispersion curves allows the identification of modes that propagate, well or not, depending on the frequency range of interest.

The analysis of arbitrarily shaped three-dimensional structures requires numerical methods. In order to reduce the computational cost inherent to full three-dimensional models, specific finite element (FE) methods have been proposed to exploit the symmetry naturally arising in waveguide problems. A first approach, referred to as the Wave Finite Element (WFE) method, allows the consideration of periodic waveguides (i.e. obeying discrete symmetry). This approach uses Bloch theorem to reduce the problem to a single repetitive unit cell -see for instance Refs. [START_REF] Mead | Wave propagation in continuous periodic structures: Research contributions from southampton, 1964-1995[END_REF][START_REF] Gry | Dynamic modelling of railway track : A periodic model base on a generalized beam formulation[END_REF][START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF][START_REF] Mencik | A substructuring technique for finite element wave propagation in multi-layered systems[END_REF][START_REF] Collet | Floquet-bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems[END_REF][START_REF] Maurin | Bloch theorem with revised boundary conditions applied to glide, screw and rotational symmetric structures[END_REF].

For the study of uniform waveguides (i.e. obeying continuous symmetry), a second approach is usually preferred although the WFE method can also be applied. This approach consists in discretizing only the cross-section, which can be of arbitrary shape, while describing analytically the direction of wave propagation. This approach has been applied under various names in the literature: the extended Ritz technique, the thin layer method for stratified waveguides in geophysics, the strip-element method, the semi-analytical finite element (SAFE) method or the scaled boundary finite element method in ultrasonics (see e.g. Refs. [START_REF] Nelson | Vibrations and waves in laminated orthotropic circular cylinders[END_REF][START_REF] Barbosa | Perfectly matched layers in the thin layer method[END_REF][START_REF] Astaneh | Efficient computation of dispersion curves for multilayered waveguides and half-spaces[END_REF][START_REF] Liu | Strip element method to analyze wave scattering by cracks in anisotropic laminated plates[END_REF][START_REF] Hayashi | Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example[END_REF][START_REF] Bartoli | Modeling wave propagation in damped waveguides of arbitrary cross-section[END_REF][START_REF] Treyssède | Three-dimensional modeling of elastic guided waves excited by arbitrary sources in viscoelastic multilayered plates[END_REF][START_REF] Gravenkamp | The computation of dispersion relations for three-dimensional elastic waveguides using the scaled boundary finite element method[END_REF]). For conciseness, it will be referred to as the SAFE method throughout this paper. Since the problem is reduced on the cross-section (two-dimensional), the number of degrees of freedom (dofs) is significantly reduced compared to a three-dimensional approach.

Both the WFE and SAFE methods lead to a quadratic eigenvalue problem in the wavenumber space. To achieve a good resolution of dispersion curves, this type of eigenproblem has to be solved at numerous frequencies. Even with the SAFE method, this process is computationally costly when the problem size increases. This typically occurs when the FE mesh has to be refined in the high-frequency regime (i.e. when the wavelength becomes small compared to the cross-section length). Other situations leading to large-size problems are those involving open waveguides (i.e. waveguides unbounded in the transverse direction) which requires to mesh part of the unbounded exterior medium [START_REF] Castaings | Finite element model for waves guided along solid systems of arbitrary section coupled to infinite solid media[END_REF][START_REF] Treyssède | Spectral element computation of high frequency leaky modes in three-dimensional solid waveguides[END_REF], or those involving multi-component structures (see e.g. [START_REF] Mencik | A substructuring technique for finite element wave propagation in multi-layered systems[END_REF][START_REF] Treyssède | Dispersion curve veering of longitudinal guided waves propagating inside prestressed seven-wire strands[END_REF]).

For a fast and accurate analysis, it is hence desirable to further reduce the size of the modal problem. There are mainly two approaches to achieve this. The first approach consists in applying a component mode synthesis method. Component mode synthesis methods operate in the frequency domain and are well suited for the WFE method in order to reduce the three-dimensional model to its two-dimensional boundary dofs, where the Bloch periodicity is imposed [START_REF] Krattiger | Bloch mode synthesis : Ultrafast methodology for elastic band-structure calculations[END_REF][START_REF] Zhou | Wave finite element method based on reduced model for one-dimensional periodic structures[END_REF][START_REF] Palermo | Extended bloch mode synthesis: ultrafast method for the computation of complex band structures in phononic media[END_REF][START_REF] Krattiger | Generalized bloch mode synthesis for accelerated calculation of elastic band structures[END_REF][START_REF] Fan | Model reduction schemes for the wave and finite element method using the free modes of a unit cell[END_REF]. This approach is yet not relevant for the SAFE method since the continuous symmetry is accounted for in the operator itself, prior to the FE discretization, leading to a two-dimensional formulation.

A second approach consists in employing a reduction basis to project the matrix operators of the eigenvalue problem. In Ref. [START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF], a reduced basis has been used from the least attenuated eigenvectors of the highest frequency of analysis. However, using a single-frequency basis for the whole band gradually alters the quality of the approximation as the frequency decreases. In Ref. [START_REF] Hussein | Reduced bloch mode expansion for periodic media band structure calculations[END_REF], a reduced basis has been proposed by concatenating the eigenvectors computed at several spectral points. In Ref. [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF], a reduced basis has been built from a set of positive-going propagative waves computed at the cut-on frequencies over the frequency band, separating the real and imaginary parts of eigenvectors. This kind of strategy has been improved in Ref. [START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF], collecting eigenvectors computed from a specific wavenumber sampling procedure, in order to cover the full complex spectrum. Inspired from Refs. [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF][START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF], a reduction basis strategy has also been applied in Ref. [START_REF] Palermo | A reduced bloch operator finite element method for fast calculation of elastic complex band structures[END_REF] for the fast calculation of complex band structures of bi-periodic media. As mentioned in these references, splitting eigenvectors into real and imaginary parts leads to a larger basis involving redundant vectors, and therefore to numerical instabilities. Filtering procedures based on modal assurance criterion or singular value decomposition must be post-processed to extract redundant vectors, and hereby, avoid ill-conditioning of the reduced problem [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF][START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF][START_REF] Palermo | A reduced bloch operator finite element method for fast calculation of elastic complex band structures[END_REF].

All the above-mentioned strategies have been implemented in the framework of WFE-like methods. The goal of this paper is to propose a model reduction strategy for the SAFE method, yielding fast determination of dispersion curves in continuously symmetric waveguides. Although the procedures described in Refs. [START_REF] Hussein | Reduced bloch mode expansion for periodic media band structure calculations[END_REF][START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF][START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF][START_REF] Palermo | A reduced bloch operator finite element method for fast calculation of elastic complex band structures[END_REF] could be adapted to the SAFE method, an alternative reduction strategy is proposed in the present paper. The reduced basis is built from complex eigenvectors computed at the maximum frequency of analysis enriched by cut-off eigenvectors. The resulting basis is then post-processed through a QR decomposition. As opposed to Refs. [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF][START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF][START_REF] Palermo | A reduced bloch operator finite element method for fast calculation of elastic complex band structures[END_REF], the reduction process does not require the solutions of full eigenvalue problems at several frequencies or several wavenumbers. Besides, the eigenvectors are not split into real and imaginary parts which avoids to control the quality of the reduced basis with user-defined tolerance criteria. The so-obtained reduction basis is rather simple to implement, yields fast and accurate computations for the full complex spectrum, including attenuated waves as will be shown in the numerical results.

The paper is organized as follows. Section 2 briefly recalls the SAFE method. Section 3 details the reduction strategy. In Sec. 4, numerical results are presented for three critical test cases. The first test case is that of a damped configuration over a large frequency range, given by a circular cylinder made of viscoelastic material. The second test case corresponds to a seven-wire strand. This test case involves strong curve veering phenomena, inherent to interwire contact mechanics. The last test case corresponds to an open waveguide involving a cylindrical core embedded into an unbounded elastic medium, truncated thanks to a perfectly matched layer (PML). In spectral theory, this kind of test case belongs to the class of open problems, which are known to be far more complex to handle than closed problems.

Background: semi-analytical finite element approach

This section briefly recalls the SAFE formulation to compute the wave modes of continuously symmetric waveguides. For more details, the reader may refer to Refs. [START_REF] Hayashi | Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example[END_REF][START_REF] Bartoli | Modeling wave propagation in damped waveguides of arbitrary cross-section[END_REF][START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF] for instance.

Let us consider a three-dimensional elastic waveguide of axis z. The cross-section coordinates are denoted as (x, y). A time Fourier transform is applied to the wave fields, with the following convention: g(ω) = +∞ -∞ g(t)e +jωt dt, where g(t) denotes an arbitrary time-dependent function. Under small strain and displacement assumptions, the elastodynamic equilibrium equations are:

∇ • σ(u) + ρω 2 u = f, (1)
where u is the displacement vector, σ(u) is the 2nd order stress tensor, ρ is the density and f is the body force vector (representing elastodynamic sources). The stress-strain relation is σ(u) = C : ǫ(u), where C is the 4th order constitutive tensor of linear elasticity. This tensor is complex in case of viscoelastic materials. The strain-displacement relation is given by ǫ(u) = (∇u + ∇u T )/2, where the superscript T denotes the transpose.

To separate the treatment of axial and cross-section coordinates, the strain-displacement relation is rewritten as follows:

ǫ = L S + L z ∂ ∂z u, (2)
where ǫ denotes the 6-dimensional strain vector using Voigt notation. The strain operator has been split into L S , containing all terms but derivatives with respect to the z-axis, and L z , the operator of z-derivatives.

A FE discretization is then applied to the cross-section coordinates. The displacement on each element is interpolated by u e (x, y, z, ω) = N e (x, y)U e (z, ω), where N e (x, y) is the matrix of two-dimensional interpolation functions and U e (z, ω) is the vector of nodal displacements of the element e. It can be shown that the cross-section discretization of the weak formulation associated with Eq. (1) yields the following global matrix system:

(K 1 -ω 2 M)U + (K 2 -K T 2 )U ,z -K 3 U ,zz = F. (3)
with the element matrices:

K e 1 = S e (L S N e ) T CL S N e dS , K e 2 = S e (L S N e ) T CL z N e dS , K e 3 = S e (L z N e ) T CL z N e dS , M e = S e
ρN eT N e dS , F e = S e

N eT fdS . ( 4
)
where C is the 6-by-6 constitutive matrix Voigt notation of the tensor C.

A spatial Fourier transform in the z-direction is finally applied to Eq. (3), which leads to the SAFE formulation:

K 1 -ω 2 M + jk(K 2 -K T 2 ) + k 2 K 3 Û = F, ( 5
)
where k is the axial wavenumber. The hat symbol denotes the spatial Fourier transform of wave fields, using the following convention: ĝ(k) = +∞ -∞ g(z)e -jkz dz, where g(z) is an arbitrary function of z. For the computation of wave modes, the source term in Eq. ( 5) is dropped ( F = 0). The source-free problem is a quadratic eigenvalue problem in terms of k. This quadratic eigenvalue problem can be transformed to a linear eigenproblem, more suitable for using efficient eigensolvers, as follows:

(A -kB)x = 0 (6)
with:

A = 0 I -(K 1 -ω 2 M) -i(K 2 -K T 2 ) , B = I 0 0 K 3 , x = Û k Û . (7)
Other linear forms are possible [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. Due to the FE nature of the discretization, the matrices involved in Eq. ( 5) are sparse. In this paper, the eigensystem (6) is solved with the ARPACK library [START_REF] Lehoucq | ARPACK User's Guide: Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods[END_REF] (eigs MATLAB implementation), appropriate for large sparse matrices and based on the implicitly restarted Arnoldi method. The calculation of dispersion curves requires solving the eigensystem repeatedly for numerous angular frequencies ω.

Due to the symmetry property of matrices K 1 , K 3 and M (see Eqs. 4), the solution of the eigenvalue problem associated with Eq. ( 5) consists of pairs of opposite wavenumbers. These pairs, denoted by (k m , U m ) and (k -m , U -m ) with the notation k -m = -k m (m = 1, . . . , M), represent M positive-going modes and M negative-going modes respectively.

The energy velocity v e of a given mode U m , which will be used for plotting the dispersion curves of this paper, can be post-processed as follows (see Ref. [START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF] for instance):

v e = 2ωIm U * m (K T 2 + ik m K 3 )U m Re U * m (K 1 + ω 2 M + ik m (K 2 -K T 2 ) + k 2 m K 3 )U m ( 8
)
where the (•) * operation represents the complex conjugate matrix transpose.

Model reduction strategy for fast dispersion analysis

Wave computation at the maximum frequency

The first step of the reduction process consists in computing the wave modes at ω max , the maximum frequency of analysis, corresponding to the following eigenvalue problem:

K 1 -ω 2 max M + jk(K 2 -K T 2 ) + k 2 K 3 U = 0. (9)
This quadratic eigenproblem is solved by using the same linear form as Eq. ( 6).

With the ARPACK library, a specified number of eigenvalues can be looked for around a user-defined shift. In this paper, this shift is set to zero. This amounts to compute the wavenumbers of smallest complex modulus, so that both positive-and negative-going waves are computed and that both real and complex waves are included in the solution.

Let us remind that for pure elastic waveguides (C real), both real and complex wavenumbers occur. Real wavenumbers correspond to pure propagating modes (Im k = 0, v e 0). Complex wavenumbers correspond to attenuated waves (Im k 0), which are non-propagating modes (v e = 0) and can be classified into evanescent (Re k = 0) or inhomogeneous (Re k 0) modes. For viscoelastic waveguides (C complex), all wavenumbers become complex but quasi-propagating waves can still be identified (these waves correspond to wavenumbers of small imaginary parts).

The number of eigenvalues specified with ARPACK, denoted as 2M ′ (M ′ positive-going plus M ′ negative-going waves), is chosen to get all the (quasi-)propagating modes occurring at the frequency ω max .

Among these 2M ′ modes, it is then desirable to filter out the most attenuated waves (typically, non-propagating waves) to reduce the size of the reduced problem. Hence, the reduced basis is constructed by selecting only the least attenuated modes, such that:

|Im k m | ≤ k ′′ max ( 10
)
where k ′′ max is a user-defined parameter. The resulting eigenvectors are collected into a matrix B max . This matrix includes both positive-going and negative-going waves. For clarity, let us denote B + and B -, the matrices of positivegoing and negative-going eigenvectors respectively:

B + = [U 1 U 2 . . . U M ], B -= [U -1 U -2 . . . U -M ] (11)
Then, B max is given by:

B max = [B -B + ] ( 12
)
The number of eigenvectors retained in B max will be denoted as 2M (M positive-going plus M negative-going waves).

To illustrate the construction of B max , a typical example of spectrum is shown in Fig. 1a for a viscoelastic waveguide. Note that another category of waveguides will be considered in this paper, corresponding to an open waveguide (i.e. buried). The eigenspectrum structure of open waveguides significantly differs from that of closed waveguides (i.e. in vacuum) and deserves a specific discussion (see Sec. 4.3).

Computation of the cut-off eigenvectors

In the second step, the cut-off wave modes are computed. By definition, these modes are such that k = 0 and can be obtained by solving the following linear eigenvalue problem: This problem is also solved with the ARPACK library (sparse matrices), with a specified number of modes denoted as I ′ . Let us denote the solutions (ω i , Γ i ).

(K 1 -ω 2 M)Γ = 0. ( 13
) (a) (b)
For pure elastic problems (lossless), the cut-off frequencies are real. As a truncating criterion, the eigenvectors whose frequency is greater than ω max will not be retained in the reduced basis. For lossy problems (e.g. viscoelastic), the matrices K 1 and/or M are complex, leading to complex cut-off eigenfrequencies ω i as well as complex eigenvectors Γ i . Similarly to the construction of B max , the most attenuated modes will be filtered out. Consequently, the selected modes for the reduction process are such that:

Re ω i ≤ ω max , |Im ω i | ≤ ω ′′ max ( 14
)
where ω ′′ max is a user-defined parameter. The cut-off eigenvectors are collected into a matrix B cut :

B cut = [Γ 1 Γ 2 . . . Γ I ] (15)
and the number of eigenvectors retained in B cut will be denoted as I.

To illustrate the construction of B cut , a typical example of spectrum is shown in Fig. 1b for a viscoelastic waveguide. As already mentioned, the case of open waveguides deserves a specific discussion and will be considered later in Sec. 4.3.

Projection matrix and reduction

As a third step, the eigenvalue problem can now be reduced. The reduction basis proposed in this paper corresponds to the concatenation of B max and B cut :

B = [B max B cut ] (16)
To avoid numerical ill-conditioning in the solution process of the eigenvalue problem, a QR decomposition of B is post-processed:

B = QR (17)
The projection matrix is then given by Q, so that the initial quadratic eigenvalue problem of Sec. 2 is reduced to:

K 1 -ω 2 M + jκ(K 2 -K T 2 ) + κ 2 K 3 ψ = 0. ( 18
)
where:

K 1 = Q T K 1 Q, K 2 = Q T K 2 Q, K 3 = Q T K 3 Q, M = Q T MQ (19)
Consequently, the quadratic eigenvalue problem has been reduced from N, the number of dofs in the FE model, to 2M + I, the number of modes retained in B. In practice, 2M + I is much lower than N. The displacement eigenvector can be post-processed from the reduced solution by:

U = Qψ (20)
Note that the solution obtained from the reduced system is denoted by (κ, U) to avoid confusion with the exact solution (k, U) of the full eigenproblem.

The solution to the reduced problem can be obtained with the same procedure as described in Sec. 2, i.e. by transformation to a linear problem solved by the ARPACK library. However, it should be noticed that the matrices involved in Eq. ( 19) are no longer sparse, so that other kinds of eigensolvers may be more efficient. In this paper, the algorithm described in Ref. [START_REF] Hammarling | An algorithm for the complete solution of quadratic eigenvalue problems[END_REF] is used instead (quadeig MATLAB implementation). Although the complete solution has to be computed with this algorithm (instead of a limited number of modes with ARPACK), the computation time is further reduced. Finally, the most attenuated solutions of the reduced problem are filtered out according to the same criterion as for B max , that is, only modes such that |Im κ m | ≤ k ′′ max are retained in the final solution.

Error evaluation

The model reduction strategy induces errors that have to be quantified compared to the exact solution of the full eigenproblem.

For a given wave mode at a given frequency, let us define the relative error on the wavenumber, complex in general, as:

ǫ k = κ m -k m k m ( 21
)
where | • | denotes the complex modulus. In practice, calculating this type of error for every mode at every frequency step is not an easy task because it requires to pair the solutions of the reduced model with those of the full model. Following Ref. [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF], the error can also be evaluated by injecting the reduced system solution (κ m , U m ) into the full eigenproblem and calculating the residual:

ǫ = ||D(κ m )U m || ||U m || (22)
where

D(κ) = K 1 -ω 2 M + jκ(K 2 -K T 2 )
+ κ 2 K 3 and || • || denotes the 2-norm of a vector. As opposed to ǫ k , the residual ǫ interestingly combines in a single value the discrepancies of both the eigenvalue (wavenumber) and the eigenvector (displacement mode shape), without knowledge of the exact solution.

Discussion

Several remarks can now be addressed. First, an efficient model reduction strategy should be able to provide the wave modes traveling in both positive and negative directions. In elastic waveguides, the negative-going eigenvectors cannot be generally determined from the positive-going ones (and reciprocally), except for particular cases involving isotropic configurations [START_REF] Treyssède | Three-dimensional modeling of elastic guided waves excited by arbitrary sources in viscoelastic multilayered plates[END_REF] or pure propagating modes under lossless assumption [START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF]. The eigenvectors in both directions are hence usually not related to each other in a straightforward manner and have to be computed by solving the eigenvalue problem itself. Note that the knowledge of wave modes in both directions is necessary to determine specific modal properties in a single direction, such as the group velocity in open waveguides [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF] or the excitability in damped waveguides [START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF][START_REF] Gallezot | A modal approach based on perfectly matched layers for the forced response of elastic open waveguides[END_REF] for the calculation of the forced response. This motivates the choice of retaining both positive-going and negative-going eigenvectors in the projection matrix B max .

Second, it can be remarked that the left projection of the reduced problem ( 18) is processed with the same reduced basis as the right projection (i.e. the matrix Q), although the full quadratic eigenvalue problem is not symmetric (K 2 -K T 2 is skew-symmetric). This choice is motivated by the structure of the eigenvalue problems. Let us write the eigenproblem for a left eigenvector V m associated with the eigenvalue k m (the adjoint problem):

V T m K 1 -ω 2 M + jk(K 2 -K T 2 ) + k 2 K 3 = 0. ( 23
)
Taking the transpose and using the symmetry of K 1 , K 3 and M leads to:

K 1 -ω 2 M + j(-k m )(K 2 -K T 2 ) + (-k m ) 2 K 3 V m = 0. (24)
The above result indeed shows that any left eigenvector V m associated with the eigenvalue k m can also be viewed as a right eigenvector, but traveling in the opposite direction -k m . This justifies that the left projection of the quadratic eigenproblem can be performed by using the same basis as the right projection in Eq. [START_REF] Zhou | Wave finite element method based on reduced model for one-dimensional periodic structures[END_REF]. Note also that the transpose of Q operates without complex conjugate. Because Q is complex, the left projection by Q T rather than by Q * allows the symmetry properties of matrices to be preserved (K 1 , K 3 and M remain symmetric after reduction). Consequently, the solution of the reduced problem [START_REF] Krattiger | Bloch mode synthesis : Ultrafast methodology for elastic band-structure calculations[END_REF] still delivers a symmetric spectrum consisting of pairs of opposite wavenumbers, as in the full problem. Furthermore, the matrix B is a complex basis. This choice differs from Refs. [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF][START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF][START_REF] Palermo | A reduced bloch operator finite element method for fast calculation of elastic complex band structures[END_REF] where a real basis is built by selecting the positive-going eigenvectors only and splitting them into real and imaginary parts, hence doubling the size of the projection vectors collection. These splitting strategies require user-defined tolerances to control the quality of the reduced basis and filter out almost redundant eigenvectors. In this paper, the construction of the reduced basis Q uses a QR decomposition, which does not require any external criterion. Besides, the strategy of Refs. [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF][START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF][START_REF] Palermo | A reduced bloch operator finite element method for fast calculation of elastic complex band structures[END_REF] requires the solutions of full eigenvalue problems at several frequencies or several wavenumbers. This can be costly for problems involving many sampling points. Conversely, the strategy proposed in this paper requires only one solution for the quadratic eigenvalue problem (the solution at ω max ) and one solution for the linear eigenvalue problem (the solution at k = 0).

As a final remark, the quality of the approximation of the reduced problem obviously depends on the subspace spanned by the projection matrix B. The enrichment of the projection matrix B + by B -and B cut is of special importance for the accuracy of the reduced eigenproblem. In the absence of mathematical proof (out of the scope of this paper), this will be illustrated by numerical tests later (see Sec. 4.2).

Summary of the reduction procedure

The reduction procedure can be summarized as follows:

• B max is computed with the ARPACK library for 2M ′ modes, where 2M ′ is chosen in order to include all the (quasi-)propagating modes,

• B max is then truncated to 2M ≤ 2M ′ modes by retaining eigenvectors such that |Im k m | ≤ k ′′ max , • B cut is computed with the ARPACK library for I ′ modes, where I ′ is chosen in order to exclude the solutions exceeding the maximum frequency of analysis, i.e. such that Re ω i > ω max ,

• B cut is truncated to I ≤ I ′ modes by retaining eigenvectors such that |Im ω i | ≤ ω ′′ max , • the reduced eigenvalue problem is fully solved iteratively (i.e. at each frequency step) with the algorithm of Ref. [START_REF] Hammarling | An algorithm for the complete solution of quadratic eigenvalue problems[END_REF],

• the solutions of the reduced problem such that |Im κ m | > k ′′ max are finally filtered out.

Results

This section compares the performance achieved by the reduction method with a full approach for three test cases. The reduction method follows the procedure summarized in Sec. 3.6. For the results of this paper, the user-defined parameters k ′′ max and ω ′′ max are set to k ′′ max a = ω ′′ max a/c s = 1, where a and c s will denote characteristic length and velocity depending on the test case. The FE meshes of this paper have been built from the free software Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. For conciseness, the dispersion curves will be usually post-processed for positive-going waves only (although it has been checked that negative-going waves accurately occur as opposite-going waves in the reduced model).

Viscoelastic cylindrical bar

As a first example, let us consider an isotropic homogeneous viscoelastic steel cylindrical bar of circular crosssection. The properties of steel are given in Table 1. With viscoelasticity, Young's modulus and Poissons ratio are complex and given by: with the complex bulk velocities cl and cs :

E = ρc 2 s 3c 2 l -4c 2 s c2 l -c2 s , ν = 1 2 c2 l -2c 2 s c2 l -c2 s ( 25
cl,s = c l,s 1 + i β l,s 2π -1 (26)
where β l and β s denote the bulk wave attenuations in Neper per wavelength. For this test case, the operators and the surface element dS involved in the SAFE matrices (see Eq. ( 4) are given by:

L S =                            ∂ ∂x 0 0 0 ∂ ∂y 0 0 0 0 ∂ ∂y ∂ ∂x 0 0 0 ∂ ∂x 0 0 ∂ ∂y                            , L z =                          0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0                          , dS = dxdy. ( 27
)
Owing to the viscoelastic properties, note that the SAFE matrices K 1 , K 2 and K 3 are complex.

The dispersion curves will be given for the normalized energy velocity v e /c s and the normalized attenuation Im(ka) as a function of the normalized frequency ωa/c s , where c s denotes the shear wave velocity of steel and a is the radius of the cylinder.

The FE mesh of the cross-section is shown in Fig. 2. Six-node triangles are used (quadratic interpolating functions), resulting in 4743 dofs. The averaged length l e of triangles satisfies the following meshing criterion: l e < λ/5, where λ = 2πc s /ω max is the smallest wavelength of the problem.

To compute B max and B cut , the number of modes have been set to 2M ′ = 350 and I ′ = 75. Discarding the most attenuated modes, i.e. such that |Im(ka)| > 1 or |Im ωa/c s | > 1, allows 2M = 152 and I = 75 modes to be retained in the bases (note that I = I ′ here because the imaginary parts of the cut-off eigenfrequencies are small owing to the weak viscoelasticity of steel). The quadratic eigenvalue problem is therefore reduced from 4743 dofs to 2M + I = 227.

Figure 3 plots the dispersion curves computed with the full model and with the reduced model up to a highfrequency regime, corresponding to a normalized frequency ω max a/c s = 10. The frequency range has been evenly divided into 100 steps. The results of the reduced model appear as superimposed to those of the full model over the whole frequency range. The attenuation in Fig. 3 is shown up to Im(ka) = 0.1, instead of 1, in order to focus on the quasi-propagating modes. Despite their relative weakness, the attenuation of these modes is also accurately determined by the reduction method. Note that the calculation of the energy velocity involves both the wavenumber and the displacement mode shape (see Eq. ( 8)). The agreement of results in Fig. 3a hence indicates that the eigenvectors are also accurately solved by the reduced model.

Figure 4 plots the wavenumber error ǫ k and the residual ǫ post-processed for all modes as a function of frequency. The wavenumber error remains lower than 0.1%, and the residual remains lower than 1%. In particular, it can be observed that both errors tend to stabilize or decrease as the frequency decreases (i.e. as it goes far away from ω max ).

Table 2 compares the computational time of the reduced model and the full model. The reduction method allows the computation time of the eigenvalue problem to be reduced by a factor 50. It can also be remarked that the time required to build the reduced system is roughly equal to the computation time of B max , that is to say only one frequency iteration of the full model. The other operations (calculation of B cut , QR factorization and projection of system) are significantly less time consuming than the calculation of B max . As already mentioned in Sec. 3.4, calculating the error ǫ k for every mode at every frequency step is not an easy task because it requires pairing the solution of both models. For convenience and since it follows the same trend as ǫ k , only the residual ǫ will be plotted in the next test cases.

Curve veering test case: seven-wire strand

The method is now applied to a more complex structure by considering a steel seven-wire strand as shown in Fig. 5a. The dynamic behavior of this kind of structure exhibits strong curve veering phenomena inherent to interwire contact mechanics and has been thoroughly studied in Refs. [START_REF] Treyssède | Dispersion curve veering of longitudinal guided waves propagating inside prestressed seven-wire strands[END_REF][START_REF] Treyssède | Investigation of the interwire energy transfer of elastic guided waves inside prestressed cables[END_REF].

Curve veering usually occurs in eigenvalue problems of weakly coupled systems [START_REF] Perkins | Comments on curve veering in eigenvalue problems[END_REF][START_REF] Pierre | Mode localization and eigenvalue loci veering in disordered structures[END_REF] and can be defined as the repulsion of two modal branches in the (k, ω) plane, suddenly veering away from each other instead of crossing. Problems involving curve veering can hence be considered as critical for reduction methods, since the reduction basis will be generally computed far away from the veering frequencies.

The material properties of steel are those of Table 1 but the viscoelasticity has been neglected for simplicity (β l,s = 0). The characteristic length a is chosen as the radius of the central wire. Because peripheral wires are not straight but helical, the analysis of guided waves requires a specific curvilinear coordinate system, called twisting coordinate system. A twisting system is a particular case of helical system, with zero curvature and a torsion τ = 2π/L (with L denoting the helix pitch of peripheral wires). In this twisting system, the cross-section of the whole structure is continuously symmetric along the z-axis so that guided waves truly exist [START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF].

The elastodynamic operators must be rewritten in the twisting coordinate system yielding expressions for L S and L z that depend on τ. These expressions are hence different from those expressed in Cartesian coordinates and given by Eq. ( 27). Furthermore, the strand is prestressed with a 0.6% elongation, which further modifies the operators. These operators are not recalled here for paper conciseness. The reader may refer to [START_REF] Treyssède | Dispersion curve veering of longitudinal guided waves propagating inside prestressed seven-wire strands[END_REF][START_REF] Frikha | Mechanical modeling of helical structures accounting for translational invariance. Part 1: Static behavior[END_REF][START_REF] Treyssède | Mechanical modeling of helical structures accounting for translational invariance. Part 2: Guided wave propagation under axial loads[END_REF] for details about the modeling of prestressed seven-wire strands. Although more involved, the formulation leads to a problem that keeps exactly the same form as in Sec. 2 (the reduction procedure described in Sec. 3 hence remains applicable).

Figure 5b shows the FE mesh of the cross-section. The finite elements are six-node triangles obeying the same meshing criterion as before (l e < λ/5), but the meshes are refined near interwire regions in order to accurately account for contact phenomena [START_REF] Treyssède | Dispersion curve veering of longitudinal guided waves propagating inside prestressed seven-wire strands[END_REF]. The radius of peripheral wires is set to 0.967a and τa = 0.0705.

The number of modes have been set to 2M ′ = 170 and I ′ = 39 for the computation of B max and B cut respectively. Discarding the attenuated modes (such that |Im(ka)| > 1 or |Im ωa/c s | > 1), 2M = 84 and I = 39 modes have been retained in the bases. Note that I = I ′ because the cut-off eigenfrequencies are real for this test case (no viscoelasticity has been introduced).

Figure 6 plots the dispersion curves computed with the full model and with the reduced model up to the normalized frequency ω max a/c s = 2. The frequency range has been evenly divided into 100 steps.

The results of both models coincide with each other, both for the energy velocity and the attenuation. As can be observed, seven-wire strands are complex structures involving many dispersion curves owing to the interwire coupling (a detailed description of these curves is out of the scope of this paper but can be found in Ref. [START_REF] Treyssède | Investigation of the interwire energy transfer of elastic guided waves inside prestressed cables[END_REF]). Figure 6a exhibits several velocity drops of modal branches as a function of frequency. These sudden changes of energy velocity, corresponding to curve veerings in the (k, ω) plane, are accurately captured by the reduced model. A typical example of veering occurs around ωa/c s =0.4 for the fastest mode, corresponding to the first compressionallike mode (the curve splits into two modes labeled L(0,1)a and L(0,1)b in Fig. 6a). This veering phenomenon has been thoroughly investigated in Ref. [START_REF] Treyssède | Dispersion curve veering of longitudinal guided waves propagating inside prestressed seven-wire strands[END_REF] and is due to the radial displacement constraint imposed between wires in the contact regions. For this test case, since the viscoelasticity has been neglected, non-attenuated modes correspond to pure propagating modes (Im k = 0, v e 0) while attenuated modes correspond to non-propagating modes (v e = 0). Although of less interest than propagating modes, the non-propagating modes are also accurately computed by the reduced model as shown by Fig. 6b.

Figure 7 plots the modal residual ǫ as a function of frequency. Similarly to the test case of Sec. 4.1, the residual remains lower than 1% and tend to stabilize as the frequency decreases. Table 3 compares the computational time of the reduced model and the full model. The reduction method allows the computation time of the eigenvalue problem to be reduced by a factor roughly equal to 80. The performance of the reduction method is even better than in Sec. 4.1 because of the mesh refinement required to accurately take into account interwire contact in the full model (here, the quadratic eigenvalue problem has been reduced from 12,369 dofs to only 2M + I = 123 modes). As in the previous test case, the time required to build the reduced system is roughly equal to the computation time of one frequency iteration of the full model.

To show the influence of cut-off eigenvectors in the quality of the reduction method, Fig. 8a shows the energy velocity dispersion curves computed with a reduced basis omitting B cut . The so-obtained reduced model completely fails to recover the results of the full model as the frequency decreases away from ω max .

Similarly, the influence of negative-going modes is illustrated in Fig. 8b, where the reduced basis has been built from B + and B cut but by omitting B -. Compared to the full model, significant discrepancies are found. Note that for this numerical test, the dispersion curves are given for positive-and negative-going modes: a closer look at the results shows that the so-obtained reduced model delivers an unsymmetric spectrum, that is to say, some negative-going modes are not accurately paired with positive-going ones. Note that this is not in contradiction with the remark of Sec. 3.5 about the symmetry of the spectrum, because such a symmetry is indeed only guaranteed for wavenumbers while Fig. 8b shows the energy velocity, which depends on the displacement mode shape. Interestingly, the asymmetry of the energy velocity spectrum can hence be related to the poor approximation of eigenvectors in the absence of the opposite-going wave subset B -. It can be concluded, from the comparison of Figs. 8a-b with Fig. 6a, that the enrichment of B + both by B -and B cut appears as an efficient reduction strategy for an accurate determination of dispersion curves over the whole frequency range of analysis.

Open waveguide test case: bar buried into a solid matrix

An open waveguide is now considered. The test case is taken from the paper of Castaings et al. [START_REF] Castaings | Finite element model for waves guided along solid systems of arbitrary section coupled to infinite solid media[END_REF]. It consists of a steel cylinder of radius a = 10 mm buried into an infinite concrete matrix. The material characteristics are given in Table 4. In the remainder, the characteristic velocity c s used to normalize the results will be the shear wave velocity of the core (steel).

For numerical purpose, the infinite surrounding medium is truncated by a finite PML of thickness h in order to absorb waves. A radial PML is introduced by analytic continuation [START_REF] Chew | A 3D perfectly matched medium from modified maxwell's equation with stretched coordinates[END_REF] into the complex transverse coordinate r:

r = r 0 γ(ζ)dζ ( 28
)
where γ is a user-defined complex-valued function of r in the absorbing region defined by

d < r ≤ d + h (d ≥ a).
For absorbing outgoing waves, the function γ must satisfy: γ(r)=1 for r ≤ d, and Im γ(r) > 0 for r > d. The userdefined parameters h and d denote the PML thickness and the radial position of the interface between the PML and the physical region respectively. In this paper, a parabolic profile independent of frequency is chosen for the radial PML function γ, expressed as follows:

γ(r) =            1 if r ≤ d 1 + 3( γ -1) r -d h 2 if r > d ( 29
)
where γ = 1 h d+h d γ(ζ)dζ denotes an averaged value of γ in the absorbing region. The introduction of a complex radial coordinate in the SAFE formulation yields expressions for L S and dS which depend on γ(r) and r(r), leading to complex matrices K 1 , K 2 , K 3 and M. For additional details on the SAFE-PML formulation, the reader may refer to Ref. [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF], where the present test case has also been used to validate the method.

Figure 9 shows the FE mesh of the cross-section. The finite elements used are six-node triangles satisfying the criterion l e < λ/5 (where λ is the shear wavelength of steel at the maximum frequency). The PML parameters are set as follows: d = a, h = a, γ = 2 + 4i. A Dirichlet condition is applied at the exterior boundary of the truncated section.

The number of modes have been set to 2M ′ = 600 and I ′ = 280 for the computation of B max and B cut respectively. Discarding modes such that |Im(ka)| > 1 or |Im ωa/c s | > 1 leads to 2M = 278 and I = 41 modes retained in the bases. The quadratic eigenvalue problem is hence reduced from 4971 dofs to 2M + I = 319 modes.

Let us emphasize that, even with a PML, the discretized eigenvalue problem keeps the same form as in Sec. 2. The reduction procedure described in Sec. 3 remains applicable and is hence quite general. However, as previously mentioned, this test case belongs to the class of open eigenvalue problems in spectral theory. Its eigenspectrum structure significantly differs from that of closed problems. Figure 10 shows the eigenspectrum of the buried bar computed for ω = ω max and that for k = 0.

From a theoretical point of view, the unbounded nature of the problem gives rise to continua of radiation modes (branch cuts). Numerically, these continua are discretized as finite sets due to the truncation of the PML to a finite thickness [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF][START_REF] Kim | The computation of resonances in open systems using a perfectly matched layer[END_REF]. In Fig. 10, the radiation modes can be recognized as clusters of poles emanating from special points (branch points), identified with black square symbols. These radiation modes, also called PML modes or Berenger modes, mainly oscillate inside the PML. Individually, a radiation mode cannot be considered as intrinsic to the physics since it depends on the user-defined PML parameters. Conversely, the so-called leaky modes indicated in Fig. 10 are intrinsic to the physics. These modes are revealed by the rotation of branch cuts, equal to -arg(d + h γ) in the complex plane [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF], due to the complex stretched coordinates of the PML. More details about the spectrum of open systems using perfectly matched layers can be found in Refs. [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF][START_REF] Kim | The computation of resonances in open systems using a perfectly matched layer[END_REF][START_REF] Gallezot | Contribution of leaky modes in the modal analysis of unbounded problems with perfectly matched layers[END_REF].

In practice, solving the eigenvalue problem leads to the computation of many radiation modes, which perturbs the visualization of dispersion curves of physical modes, namely the leaky modes (complex poles) and the trapped modes (real poles, if any) [START_REF] Nguyen | Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods[END_REF][START_REF] Treyssède | Finite element computation of trapped and leaky elastic waves in open stratified waveguides[END_REF]. Modal filtering must be post-processed to identify and remove the radiation modes. The filtering criterion used in this paper is based on the ratio of the imaginary part over the modulus of the kinetic energy [START_REF] Gallezot | A modal approach based on perfectly matched layers for the forced response of elastic open waveguides[END_REF]:

1 - Im E k |E k | > η min (30)
where E k = ω 2 4 U * m MU m is the kinetic energy, complex with the PML, and 0 < η min < 1 is a user-defined parameter. In this test case, the numerical results are shown after modal filtering with the parameter η min being set to 0.15.

Figure 11 plots the normalized dispersion curves of leaky modes obtained with the full model and with the reduced model, computed for a frequency range divided into 100 steps (for leaky modes, note that the energy velocity is integrated on the core only [START_REF] Bernard | Guided waves energy velocity in absorbing and non-absorbing plates[END_REF][START_REF] Pavlakovic | Disperse User's Manual[END_REF], so that the vectors and matrices involved in Eq. ( 8) must be understood as restricted to the dofs of the core [START_REF] Treyssède | Spectral element computation of high frequency leaky modes in three-dimensional solid waveguides[END_REF][START_REF] Gallezot | A modal approach based on perfectly matched layers for the forced response of elastic open waveguides[END_REF]). Figure 12 plots the residual of modes as a function of frequency. As in the previous test cases, the results of both models appear as superimposed, with a residual lower than 1%.

Table 5 compares the computational time of the reduced model and the full model. As previously noticed for the first two test cases, the time required to build the reduced system remains roughly equal to the computation time of one frequency iteration of the full model. The computation time of the eigenvalue problem is tremendously reduced, by a factor roughly equal to 90.

It is noteworthy that with the full model, a large number of PML modes is required in practice for determining only a limited number of leaky modes. The number of modes computed at each frequency step (2M ′ = 600 here) significantly increases the computation cost of eigensolvers [START_REF] Treyssède | Spectral element computation of high frequency leaky modes in three-dimensional solid waveguides[END_REF]. Most of the PML modes are yet highly attenuated and hence of no interest. With the reduction procedure, these modes are advantageously filtered out from the reduced basis B (yielding 2M + I = 319 modes for this test case). As a final remark, even if radiation modes are of little interest in practice, the reduction strategy proposed in this paper allows to capture them. As an example, Fig. 13 compares the eigenspectrum computed with the reduced model and the full model for a given frequency ω, different from the maximum frequency of analysis (here, ω = ω max /2). Although slight discrepancies can be observed for a few of them, the radiation modes turn out to be satisfyingly obtained by the reduced model also.

Conclusion

A reduction strategy has been proposed to fastly and accurately solved the quadratic eigenvalue problem inherent to the finite element analysis of continuously symmetric waveguides, including complex attenuated waves. The method is based on a specific modal reduced basis, built from the complex eigenmodes computed at the maximum frequency of analysis of the problem, including both positive-and negative-going complex waves, and enriched by cut-off eigenvectors. The resulting basis is further reduced by excluding the most attenuated modes and is simply postprocessed by a QR decomposition for well-conditioning of the reduced system. The reduction method is quite general and has been tested numerically by considering various critical test cases, considering high-frequency viscoelastic loss, strong curve veering, or leakage loss (buried waveguide). The numerical results show that the construction of the reduced model remains fast, of the order of a single frequency iteration of the full eigenvalue problem, and that the computational speed of the reduced eigenvalue problem can be increased by a factor greater than 50 compared to the full problem. Applying the reduction method appears to be particularly advantageous for solving problems involving mesh refinements, for instance required by high-frequency analyses or contact mechanics, as well as for solving open waveguide problems with perfectly matched layers. Although not considered in the present paper, it is noteworthy that the procedure could also be readily applied to periodic waveguides using the WFE method.

Figure 1 :

 1 Figure 1: Example of eigenspectrum for a viscoelastic waveguide. (a) Eigenspectrum for ω = ω max , (b) Eigenspectrum for k = 0 (note that Im ω i ≤ 0 according to the time Fourier transform convention adopted in this paper). Modes computed with the ARPACK library are inside the dashed red circle. The least attenuated modes retained in the reduced basis are inside the blue shaded region. The axes have been normalized (a and c s denote characteristic length and celerity respectively).

Figure 2 :

 2 Figure 2: FE mesh for the cylindrical bar (4743 dofs).

Figure 3 :

 3 Figure 3: Dispersion curves in a cylindrical bar for: (a) the energy velocity, (b) the attenuation. •: reduced model, •: full model. Black bullets are smaller than red ones in order to visualize the superimposition of both models.

Figure 4 :

 4 Figure 4: Error as a function of frequency for the cylindrical bar: (a) wavenumber error ǫ k , (b) residual ǫ.

Figure 5 :

 5 Figure 5: (a) Three-dimensional view of a seven-wire strand, (b) FE mesh of cross-section (12369 dofs).

Figure 6 :

 6 Figure 6: Dispersion curves in a seven-wire strand. Same legend as in Fig. 3.

Figure 7 :

 7 Figure 7: Residual ǫ as a function of frequency for the seven-wire strand.

Figure 8 :

 8 Figure 8: Numerical tests. Energy velocity dispersion curves of the seven-wire strand computed with a reduced basis: a) without B cut , (b) without B -(•: reduced model, •: full model).

Figure 9 :

 9 Figure 9: FE mesh for the buried bar (4971 dofs). Gray: steel cylindrical core, black: concrete matrix (PML).

Figure 10 :Figure 11 :

 1011 Figure 10: Eigenspectrum for the buried bar (open waveguide). (a) Eigenspectrum for ω = ω max (black squares: theoretical branch points, occurring at k = ±ω/c l∞ and k = ±ω/c s∞ , where c l∞ and c s∞ denote the bulk longitudinal and shear wave velocities of the external medium). (b) Eigenspectrum for k = 0 (black square: theoretical branch points, occurring at ω = 0). Modes computed with the ARPACK library are inside the dashed red circle. The least attenuated modes retained in the reduced basis are inside the blue shaded region.

Figure 12 :

 12 Figure 12: Residual ǫ as a function of frequency for the buried bar.

Figure 13 :

 13 Figure 13: Buried bar. Eigenspectrum computed at ω = ω max /2 (•: reduced model, •: full model). Encircled poles correspond to leaky modes (4 positive-and 4 negative-going leaky modes), the other poles correspond to PML modes.

Table 1 :

 1 Material properties for the viscoelastic cylindrical bar and the seven-wire strand.

	)

Table 2 :

 2

		Model size B max computation B cut computation QR factor. Projection EVP solution
	Full model	4743 dofs	-	-	-	-	19.8 s × 100
	Reduced model	227 dofs	17.4 s	1.34 s	0.046 s	1.74 s	0.391 s × 100

Cylindrical bar. Computational time required by the full model and the reduced model. The time for solving the eigenvalue problem (EVP) has been averaged over the frequency iterations (100 iterations).

Table 3 :

 3 Model size B max computation B cut computation QR factor. Projection EVP solution Seven-wire strand. Computational time required by the full model and the reduced model. The time for solving the eigenvalue problem (EVP) has been averaged over the frequency iterations (100 iterations).

	Full model	12,369 dofs	-	-	-	-	8.41 s × 100
	Reduced model	123 dofs	8.35 s	1.29 s	0.056 s	1.51 s	0.109 s × 100

Table 4 :

 4 of the elastodynamic equilibrium equations Material properties for the steel bar buried in concrete.

		Density	Longitudinal velocity Shear velocity
	Material ρ (kg/m 3 )	c l (m/s)	c s (m/s)
	Steel	7932	5960	3260
	Concrete	2300	4222.1	2637.5

Table 5 :

 5 Model size B max computation B cut computation QR factor. Projection EVP solution Buried bar. Computational time required by the full model and the reduced model. The time for solving the eigenvalue problem (EVP) has been averaged over the frequency iterations (100 iterations).

	Full model	4971 dofs	-	-	-	-	90.6 s × 100
	Reduced model	319 dofs	101.6 s	12.7 s	0.100 s	3.94 s	0.981 s × 100