Matthieu Alexandre

Amine, Silouane Carlos

De part leurs propriétés d'amortissement, les matériaux viscoélastiques sont souvent employés dans l'industrie automobile. C'est le cas pour les pare-brises qui sont constitués de deux feuilles de verre séparées par une ou plusieurs couches de matériaux viscoélastiques. A l'origine, la présence de ces matériaux visait exclusivement à retenir les bris de verre en cas d'accidents. Avec le temps, on s'est aperçu qu'ils permettaient de réduire les vibrations du pare-brise et par conséquent, le niveau sonore dans l'habitacle. Suite à cela, le pare-brise est devenu un élément central dans la conception acoustique des véhicules.

Les propriétés mécaniques des matériaux viscoélastiques varient en fréquence, ce qui implique (i) d'adapter les méthodes numériques généralement utilisées en conception et (ii) de connaître les propriétés mécaniques à chaque pas de fréquence. En effet, les algorithmes classiques ne sont pas adaptés aux systèmes mécaniques composés de matériaux viscoélastiques. D'autres méthodes ont été développées mais peuvent être très coûteuses ou non adaptées à certains matériaux. Par ailleurs, les matériaux utilisés dans les parebrises font l'objet de nombreux brevets, et leurs propriétés sont souvent des données sensibles que l'on ne souhaite pas divulguer.

Cette thèse répond à ces deux problématiques en proposant une méthode originale permettant de réduire les modèles numériques de structures complexes avec des propriétés dépendantes de la fréquence. Les effets d'amortissements sont d'abord pris en compte par un modèle de Golla-Hughes-McTavish qui, en ajoutant des nouvelles variables, renvoie un système où les coefficients ne dépendent plus de la fréquence. Ensuite, un algorithme basé sur la Proper Orthogonal Decomposition permet de réduire efficacement la taille du modèle sans perte notable d'information. Cette procédure permet deux choses: un calcul très rapide de réponses en fréquence, et la création de super-éléments, modèles réduits destinés être connectés à d'autres modèles numériques. i

Numerical simulation is a powerful tool largely used in the design of many products. It allows to simulate the deformation of a structure, perform parametric studies or predict an acoustic field and so on. In doing so, many data can be obtained without the need to resort to any experiments, a trend every manufacturer try to follow.

In the automotive industry, the main function of planes, trains or cars is to transport people from one point to another. However, many physical phenomena make these structures vibrate and generate undesirable sounds into the passenger compartment. In order to avoid these noises, engineers perform vibroacoustic simulations with the aim of improving the acoustic comfort. All of these simulations are based on numerical models of the different elements constituting a vehicle. In this thesis we are specially interested in the windscreen which, due to its geometry, is mostly responsible for the noise in the passenger compartment.

Laminated structure

Unlike the other glass elements equipping a vehicle, a windscreen is composed of two sheets of glass bonded by one (Figure I.1a) or more (Figure I.1b) plastic interlayers often made of PolyVinyl Butyral (PVB). This viscoelastic material is aiming at holding the shards together when the glass breaks, improving safety and security in case of accidents. More than simple "wind-shields", these laminated glass are real safety features in vehicles. A vehicle is subject to many sources of vibrations [START_REF] Liang | The evaluation of vehicle interior impact noise inducing by speed bumps based on multi-features combination and support vector machine[END_REF] including engine, transmission and high-speed wind noises. Due to its small thickness and its large surface, the windscreen easily vibrates and then radiates sound into the passengers compartment. It is therefore necessary to control the amount of vibration in order to improve the acoustic comfort, hence the solutions developed by engineers. Among them, working on the damping properties can lead to great results and the PVB turns out to play a key role. Over the past decades, engineers have improved the mechanical properties of the PVB in order to maximize damping effects. They have also developed other polymers (acoustic PVB) resulting in more damped windscreens. Then, laminated glass also ensure acoustic performances beyond the safety aspect.

Windscreen modeling

Doing numerical simulations involves having numerical models, which often come from the Finite Element Method [START_REF] Taylor | The Finite Element Method : Its Basics and Fundamentals[END_REF][START_REF] Macneal | Finite Elements: Their deigne and Performances[END_REF][START_REF] Petyt | Introduction to finite element vibration analysis[END_REF]. When windscreens have to be modeled, difficulties occur because of the small thickness of the layers. As shown in Figure I.2, the elements are very stretched, which could lead to inaccurate models. This is specially true for windscreens since the viscoelastic core is much more flexible than glass and is more subject to shear deformation as represented in Figure I.3a. Many solutions allow to properly model the shear deformation, for instance the use of quadratic finite elements or the ZPST multilayered shell element developed by Sulmoni et al. [START_REF] Sulmoni | Modal validation of sandwich shell finite elements based on ap-order shear deformation theory including zigzag terms[END_REF] and adapted by Chazot et al. [START_REF] Chazot | Harmonic response computation of viscoelastic multilayered structures using a ZPST shell element[END_REF] to viscoelastic and poroelastic structures. Based on the p-order shear deformation theory, a zigzag function is added in order to better fit the displacement field into each layers (see Figure I.3). It has the benefit of only requiring to mesh the midplane of the geometry, since the thickness is integrated into the element formulation. However, in practice, for simplicity and for historical reasons, manufacturers still prefer classical 3D linear elements for windscreen modeling. The sheets of glass are essentially bending while the PVB interlayer is subject to shear deformation.

Reduced order model of windscreens

The numerical model of the standard windscreen assembled on the "Renault Zoé", represented in Figure I.2, is made of 1 268 linear 3D elements, 10 525 nodes for a total of 31 575 variables. A classical study consists in computing the displacement or acceleration at a given location on the structure, while another is subject to an external harmonic load. In doing so for a given set of frequencies, one computes the so called frequency response function. Then, thousands of equations must be solved several times, which becomes heavy when the number of frequencies is large. This motivates to reduce the number of equations to solve without compromising results accuracy, which is the aim of any model reduction technique.

A very popular method consists in representing the dynamics of a structure through a basis of modes found from the eigenvectors of the matrix system of equations. As this procedure has been largely used for the past decades, many efficient algorithms have been developed and are now the basis for model reduction methods implemented in software. Nevertheless, they only deal with materials with constant stiffness so difficulties arise when viscoelastic materials are considered. In fact, the stiffness of such materials depends on several parameters, including the frequency, which makes the classical modes superposition unfeasible as it is. In practice, several variants have been proposed and allow to build a reduced order model even though the stiffness is not constant but some drawbacks still occur: the algorithms are more complex, commercial software do not often provide such procedures, and finally the elastic modulus must be known at each frequency step, even with the reduced model.

Confidentiality issues about the PVB

The fact that the elastic modulus must be evaluated at each frequency prevents the manufacturers from sharing reduced models of windscreens since the mechanical properties related to the PVB are often confidential. Many patents about windscreens are registered and mainly concern the PVB. That is why the data related to the polymer are sensitive information that manufacturers do not want to share with partners. The present procedure consists in giving only general information as for instance the geometry, the kind of materials used and sometimes the finite element mesh. The user can then perform numerical simulations but only with a limited dataset. In doing so, every vibro-acoustic simulations considering windscreens models (reduced or not) are not always accurate and reliable.

This brings us to the aim of the thesis: build accurate reduced models of windscreens that do not reveal the confidential data related to the PVB and that facilitates the collaborations with other partners.

The following sections contain all relevant theoretical background forming the basis for the methodology proposed in this manuscript. Starting with a presentation of the windscreen structures, it then focuses on how to account for viscoelastic materials in numerical models. Next, a brief theory on dynamic substructuring methods is presented before ending on the POD-based model reduction techniques.

II. 1 Damping modeling

As Crandall said in the abstract of [START_REF] Crandall | The role of damping in vibration theory[END_REF], damping is responsible for the eventual decay of free vibrations and provides an explanation for the fact that the response of a vibratory system excited at resonance does not grow without limit. Indeed, a system with only mass and stiffness terms could eternally vibrates after being excited. Moreover, mechanical systems are subjected to resonance phenomenon which, when damping effects are ignored, leads to infinite displacements as shown in Figure II. 1. In order to avoid such behaviors, damping effects must be taken into account in systems. The most commune ways to model damping are gathered in the following paragraphs, while a more general study is available in [START_REF] Nashif | Vibration Damping[END_REF].

II. 1. 1 Viscous damping

Consider the matrix system in the frequency domain

-ω 2 M + K q = f , (II.1)
where M and K are the mass and stiffness matrices, f the load vector, q the displacement vector and ω the angular frequency. These equations can be decoupled using the eigenvalues λ j = ω 2 j and eigenvectors ϕ j , solutions of the eigenvalue problem Kϕ j = Mϕ j λ j .

(II.2)

Section 1 -Damping modeling

Storing the whole eigenvectors in the matrix Φ = [ , ..., ϕ j ... ], the transformation q = Φµ leads to -ω 2 Φ T MΦ + Φ T KΦ µ = Φ T f . (II. 3) with µ the new variable vector corresponding to the modal coefficients. When Φ is normalized such that Φ T MΦ = I, it also gives Φ T KΦ = Ω 2 = diag(ω 2 j ). Then, the previous matrix system is equivalent to the set of modal equations -ω 2 + ω 2 j µ j = ϕ T j f . (II. 4) It is clear that this undamped model leads to infinite modal coefficients when the angular frequency and the eigenvalue coincides, such that lim

ω→ω j |µ j | = lim ω→ω j ϕ T j f (-ω 2 + ω 2 j ) = ∞.
(II.5)

A modal viscous damping term can be added to each modal equations yielding -ω 2 + i ω2ζ j ω j + ω 2 j µ j = ϕ T j f (II. 6) where the modal damping factor ζ j must be measured or calculated in order to match the amplitude at resonance. Nevertheless, when no simulation data are available, ζ j is often set to a fixed arbitrary value, say 0.005, for each mode. A more sophisticated method consists in considering a viscous matrix D such that the equations of motions become

-ω 2 M + i ωD + K q = f , (II.7)
which, once the transformation q = Φµ applied, can be written as -ω 2 I + i ωΦ T DΦ + Ω 2 µ = Φ T f . (II.8)

Unfortunately, the matrix Φ does not diagonalize D, so that Φ T DΦ is a dense matrix, which makes the equations coupled. However, for lightly damped structures with sufficiently separated resonances, one can assume that only the diagonal of Φ T DΦ is essential so the modal equations become uncoupled again [START_REF] Imbert | Analyse des structures par éléments finis[END_REF], with approximated viscous terms defined by 2ζ j ω j ≈ ϕ T j Dϕ j .

(II.9)

One can also find the complex modes of equation (II. [START_REF] Benfield | Vibration Analysis of Structures by Component Mode Substitution[END_REF]) by first considering the equivalent first-order system i ωEx = Ax + F, (II. [START_REF] Berkooz | The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows[END_REF]) .11) As explained in [START_REF] Françoise | Meerbergen. The Quadratic Eigenvalue Problem[END_REF], N denotes any non-singular matrix and a natural choice is to take N = M in order to keep working with symmetric matrices. This linear system can be diagonalized after solving the associated eigenvalue problem Aχ j = Eχ j λ j .

with E = D M N 0 , A = K 0 0 N , F = f 0 and x = q i ωq . ( II
(II.12)

Here, the eigenvalues are λ j = i ω j and the corresponding eigenvectors come in conjugate pairs with the particular structure χ j = ψ j λ j ψ j .

(II. [START_REF] Bouayed | Finite element analysis of the dynamic behavior of a laminated windscreen with frequency dependent viscoelastic core[END_REF])

The pencil (λ j , ψ j ) correspond to the eigenvalue and eigenvector of the second-order system such that λ 2 j M + λ j D + K ψ j = 0.

(II. 14) Then, there is a way to find the modes of a damped system. However, it requires to solve an eigenproblem twice bigger than the original one.

II. 1. 2 Structural damping

An other way to model damping effects is to consider a complex stiffness for each material employed in the system. It is then a matter of structural or hysteretic damping. After establishing the equations respecting q(t) = |q|e +i ωt , the matrix system can be written as

-ω 2 M + K + i H q = f , (II. [START_REF] Caughey | Classical Normal Modes in Damped Linear Dynamic Systems[END_REF] with H the structural damping matrix. Considering the complex symmetric (not Hermitian) stiffness matrix K g = K + i H, one can solve the eigenvalue problem K g ϕ j = Mϕ j λ j (II. [START_REF] Chahlaoui | Second-order balanced truncation[END_REF] leading to the complex eigenvectors ϕ j and complex eigenvalues λ j = ω j (1 + i η j ), ( II.17) where ω j = Re(λ j ) and η j = Im(λ j )/Re(λ j ) represents the loss factor of the j th mode. A particular case can be observed when only one material is used. In that case, the stiffness can be factored out such that K = E(1 + i η)K f (II. [START_REF] Craig | A review of time-domain and frequency-domain component mode synthesis method[END_REF] where the scalar E(1 + i η) is the complex constant stiffness used in the model. It can be shown that all the eigenvectors are real, whereas the eigenvalues are complex with the same loss factor, i. e. η j = η. This particular situation amounts to choose the identical modal damping coefficient for each eigenmodes.

II. 1. 3 Rayleigh damping

Under particular conditions, the matrix D can be diagonalized through the eigenvectors Φ of the undamped system. A damping model, known as the Rayleigh damping, proportional or classical damping consists in expressing the damping matrix D as a linear combination of mass and stiffness matrices so that

D = α 1 M + α 2 K,
(II. [START_REF] Roy | Coupling of Substructures for Dynamic Analyses[END_REF] with α 1 , α 2 two real scalars. In doing so, the transformation q = Φµ yields [START_REF] Crandall | The role of damping in vibration theory[END_REF] This method is attributed to Lord Rayleigh (1877) [START_REF] Rayleigh | The Theory of Sound[END_REF] and a general form has been developed by Caughey and O'Kelly in 1965 [14,[START_REF] Caughey | Classical Normal Modes in Damped Linear Dynamic Systems[END_REF] which shows that

-ω 2 I + i ω(α 1 I + α 2 Ω 2 ) + Ω 2 µ = Φ T f (II.
D = M N -1 k=0 α k M -1 K k (II.21)
is the necessary and sufficient condition to allows to diagonalize the damping matrix through the undamped modes. Rayleigh damping is then a special case of the so called Caughey series.

II. 1. 4 Linear viscoelastic damping

Since viscoelastic materials are very widespread in mechanical structures, numerical models have to correctly take into account their damping effects. To do so, the previous damping models turn out to be unsuited since the stiffness is now complex-valued and frequency-dependent, resulting in a complex and frequency-dependent stiffness matrix. In the Laplace domain, with s = i ω, the matrix system can be written

s 2 M + E(s) K q = f , (II.22)
where here a single material is considered with E(ω) its elastic modulus.

Let us briefly recall how viscoelastic material behaves under stress and strain solicitations. For more complete information, the reader is referred to [START_REF] Salençon | Viscoélasticité pour le calcul des structures[END_REF]. Two experiments can be carried out to determine the evolution of the elastic moduli (Young's, shear and bulk) in the time domain. The first one consists in applying a stress step at t = t 0 and observing the strain ϵ(t). The Figure II.2 shows that the strain increases gradually after a step at t = t 0 . One can define the function J describing the evolution of ϵ such that ϵ(t) = σ 0 J(t), ∀ t > t 0 .

(II. The opposite experiment can be done by applying a strain step and observing the evolution of the stress, as shown in Figure II.3. Here, a step is visible at t = t 0 and the stress disappears with time. The function G describes its evolution such that σ(t) = ϵ 0 G(t), ∀ t > t 0 .

(II. The instant responses, namely the steps at t = t 0 , are testament to the purely elastic part of the material and directly give the elastic modulus such as

J(t 0 ) = ϵ 0 σ 0 = 1 E 0 , G(t 0 ) = σ 0 ϵ 0 = E 0 , (II.25)
Thanks to the theory of distribution and the Stiltjes integrals, the constitutive stressstrain relation is given by

σ(t) = +∞ -∞ G(t -τ )dϵ(τ ). (II.26)
Considering σ(t < t 0 ) = 0, this integral may be replaced by

σ(t) = T t 0 G(t -τ ) ∂ϵ(τ ) ∂τ dτ. (II.27)
The function G(t) may be written as a sum of a constant and time-dependent part

G(t) = G ∞ + h(t), with G ∞ = lim t→∞ G(t), (II.28) leading to σ(t) = G ∞ ϵ(t) + T t 0 h(t -τ ) ∂ϵ(τ ) ∂τ dτ. (II.29) Applying the Laplace transform f (s) = Lf (t) = ∞ 0 f (t)e -st dt on the previous relation yields σ(s) = G ∞ + s h(s) ε(s). (II.30)
One can now define the complex frequency-dependent elastic modulus expressed in the frequency-domain

E(i ω) = E ′ (i ω) + i E ′′ (i ω) = G ∞ + s h(s) (II.31) with E ′ (i ω) = G ∞ +Re(s h(s)
) and E ′′ (i ω) = Im(s h(s)) are respectively called the storage modulus and the loss modulus. As any complex numbers, it can be recast as a magnitudephase form such that

E(i ω) = |E(i ω)| e i θ (II.32)
where

|E(i ω)| = E ′ (i ω) 2 + E ′′ (i ω) 2 and tan(θ) = E(i ω) ′′ E(i ω) ′ = η(i ω), (II.33) ωt ϵ(t) σ(t) θ
Figure II.4: Phase difference between stress and imposed sinusoidal strain on a viscoelastic material.

with η(i ω) the frequency-dependent loss factor. The angle can be put in evidence in plotting the stress and strain in the time further to a strain solicitation, as shown in Figure II. 4.

Only the frequency-dependency is taken into account in this manuscript, but it is important to state that viscoelastic materials depend on other parameters as pressure, pre-stress, amplitude of excitation, etc. And also temperature. Since polymers are often characterized once incorporated in their environment, engineers often study the evolution of the complex modulus with respect to the frequency and the temperature only, which are proving to be strongly linked. As shown in Figure II.5, the modulus varies in the same way with increasing frequency or decreasing temperature, and vice versa. From this statement, the principle of frequencytemperature superposition [START_REF] Williams | The Temperature Dependence of Relaxation Mechanisms in Amorphous Polymers and Other Glass-forming Liquids[END_REF] can be written as

ρT E(ω, T ) = ρ 0 T 0 E(α(T, T 0 ) ω) (II.34)
where the shift factor α(T, T 0 ) is generally determined with the William-Landel-Ferry law (WLF)

log 10 α(T, T 0 ) = -C 0 1 (T -T 0 ) C 0 2 + T -T 0 . (II.35)
The parameters C 0 1 and C 0 2 depend on the material and the temperature T 0 often taken as T 0 = T g , with T g the glass transition temperature.

II. 1. 5 Simple rheological models

Since viscoelastic materials combine purely elastic and purely viscous behaviors, some models are aiming at assemble springs and dampers (Figure II.6) in order to find a good representation of the elastic modulus. The simplest models are detailed in the following paragraphs . 4 Maxwell model. Represented in Figure II.7a, the Maxwell model [START_REF] Clerk | On the dynamical theory of gases[END_REF] is made of one spring of stiffness k and one damper of constant η in series. They are subject two the same stress σ but have their own strain denoted by ϵ s and ϵ d for the spring and the damper respectively. In series, the total strain is equal to the sum of the strain of each element and the stress is the same in each element. Then, it leads to the first two relations

k (a) Spring.
ϵ = ϵ s + ϵ d σ = σ s = σ d . (II.36)
The combination of these equations leads, in the frequency domain, to

σ = i ωkη k + i ωη , (II.37)
which gives a complex modulus written in a general form

E = K i ωτ 1 + i ωτ (II.38)
where K = k is the global stiffness and τ = η k the relaxation time.

Kelvin-Voigt model. Using a parallel association, the Kelvin-Voigt model [START_REF] Meyers | Mechanical Behavior of Materials[END_REF] assembles a spring and a damper in parallel as shown in Figure II.7b. This way, the strain is the same s in each element and the stress is the sum of the stress of each element. Then, .42) where τ = k η is the relaxation time, and K 0 and K ∞ are the modulus at low and high frequency.

σ = σ s + σ d ϵ = ϵ s = ϵ d . (II.
E = K 0 + K ∞ i ωτ 1 + i ωτ , ( II

II. 1. 6 Complex rheological models

The three previous models are often too simple to give a good representation of the real-life behaviors of viscoelastic materials. Other ones, more sophisticated, have been developed to alleviate this issue. The three principal are the generalized Maxwell [START_REF] Weichert | Ueber elastische Nachwirkung[END_REF], Anelastic Displacement Fields (ADF) [START_REF] Lesieutre | Finite element modeling of frequencydependent material damping using augmenting thermodynamic fields[END_REF], the Golla-Hughes-McTavish (GHM) [START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF] and fractional Zener [START_REF] Bagley | Fractional Calculus -A Different Approach to the Analysis of Viscoelastically Damped Structures[END_REF] models. They are all represented in the figures II.8, II.9, II.10 and II.11, while their associated complex modulus and their number of parameters are given in Table II. [START_REF] Bagley | Fractional Calculus -A Different Approach to the Analysis of Viscoelastically Damped Structures[END_REF].

II. 1. 7 Conclusions

This section introduced several approaches on how to take into account the damping properties of a mechanical system in equations. Four types of damping are presented, namely the viscous, structural, proportional and linear viscoelastic models. Since materials with constant mechanical properties could be well represented with a viscous, structural of proportional damping, special materials -as the polymer used in windscreens -require more complex models.

k 0 k 1 η 1 k 2 η 2 k n η n Figure II.8: Generalized Maxwell model. k 0 k 1 η 1 k 2 η 2 k n η n Figure II.9: Anaelastic Displacement Field (ADF) model. k 0 k 1 m 1 η 1 k 2 m 2 η 2 k n m n η n Figure II.10: Golla-Hughes-McTavish (GHM) model. k ∞ -k 0 τ α (k ∞ -k 0 ) k 0 Figure II.
(ω) = E 0 1 + n j=1 γ j i ωτ j 1 + i ωτ j 2n + 1 ADF E(ω) = E 0 1 + n j=1 ∆ j ω 2 + i ωΩ j ω 2 + Ω 2 j 2n + 1 GHM E(ω) = E 0 1 + n j=1 α j -ω 2 + i ω2ζ j ω j -ω 2 + i ω2ζ j ω j + ω 2 j 3n + 1 Fractional Zener E(ω) = E 0 + E ∞ (i ωτ ) α 1 + (i ωτ ) α 4
Table II.1: Complex moduli of the sophisticated rheological models.

II. 2 The Golla-Hughes-McTavish technique

The GHM technique, developed by David F. Golla, Peter C. Hughes and Donald J. McTavish [START_REF] Golla | Dynamics of Viscoelastic Structures-A Time-Domain, Finite Element Formulation[END_REF][START_REF] Mctavish | The mini-oscillator technique : A finite element method for the modelling of linear viscoelastic structures[END_REF][START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF], is used to take into account the damping effects in the motion equations of a system made of viscoelastic materials. Each complex and frequency-dependent moduli can be approximated by the n-order rational fraction

E(ω) = E 0 1 + n j=1 α j -ω 2 + i ω2ζ j ω j -ω 2 + i ω2ζ j ω j + ω 2 j .
(II. [START_REF] Meyer | Balancing and model reduction for secondorder form linear systems[END_REF] This expression aims at fitting the master curves of the real and imaginary parts of an elastic modulus obtained thanks to experimental measures from a Dynamic Mechanical Analyzer (DMA). For example, the real part, imaginary part and loss factor of a 1-order GHM model are illustrated in Figure II. [START_REF] Besselink | A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control[END_REF]. The storage modulus (real part) follows the same evolution as the one depicted in Figure II.5 showing the two constant regions separated by the glassy transition. For other materials, a 1-order GHM model could be insufficient so high-order models might be required.

f (Hz)

Loss factor

Re(E)

Im(E)

Figure II.12: Complex frequency-dependent stiffness k(ω) of a 1-order GHM model.

Once the material moduli are approximated, it is possible to recast the equations of motion having a complex frequency-dependent matrix as a real second-order matrix system in either the time or frequency domain. Nevertheless, there is a price to pay since new variables, called the dissipation coordinates, are introduced and increase with the order of the model. However, according some criteria (dimensions of the model, complexity of the material, level of accuracy desired, etc), numerical tools could still manage such systems due to the large computation power presently available.

II. 2. 1 One-dimensional spring-mass example

This paragraph considers the single variable spring-mass system drawn in Figure II.13a.

m k(ω) q f (a) Original spring-mass system. k α 1 k α 1 ω 2 1 k 2α 1 ζ 1 ω 1 k z 1 m q f (b) Equivalent GHM 1-order system.
Figure II.13: Study of a single variable spring-mass system modeled by a 1-order GHM model.

The equation of motion associated to such a system is

-ω 2 m + k(ω) q = f, (II.44)
where the frequency-dependent and complex stiffness k(ω) is modeled with a 1-order GHM mode, i. e. taking n = 1 in equation (II.43). The "equivalent" spring-mass-damper system is drawn in Figure II.13b. Then, the stiffness can be rewritten as

k(ω) = k 0 1 + α 1 -ω 2 + i ω2ζ 1 ω 1 -ω 2 + i ω2ζ 1 ω 1 + ω 2 1 ,
(II. [START_REF] Montier | Balanced Proper Orthogonal Decomposition Applied to Magnetoquasi-Static Problems Through a Stabilization Methodology[END_REF] where (k 0 , α to the numerator of the fraction, this expression can be recast as a sum of a constant and frequency-dependent terms, such as

k(ω) = k 0 1 + α 1 -k 0 α 1 ω 2 1 -ω 2 + i ω2ζ 1 ω 1 + ω 2 1 . (II.46)
At this stage, one can introduce the so called dissipation coordinate

z = γ 1 (ω)q, with γ 1 (ω) = ω 2 1 -ω 2 + i ω2ζ 1 ω 1 + ω 2 1 , (II.47)
which is a scalar depending on the physical displacement q, and rewrite the motion equation such that

-ω 2 m + k 0 1 + α 1 q -k 0 α 1 z = f, (II.48)
which forms one of the two main results that will lead to the final system. The second one comes from the relation z = γ 1 (ω)q that can be transformed as

-ω 2 + i ω2ζ 1 ω 1 + ω 2 1 z -ω 2 1 q = 0. (II.49)
Combining equations (II. [START_REF] Taylor | The Finite Element Method : Its Basics and Fundamentals[END_REF]) and (II.49), one can built the following second-order augmented matrix system

-ω 2 m 0 0 1 + i ω 0 0 0 2ζ 1 ω 1 + k 0 (1 + α 1 ) -α 1 k 0 -ω 2 1 ω 2 1 q z = f 0 , (II.50)
which contains only real and constant matrices. The so called GHM method consists then in using the GHM rheological model to fit damping experimental data and recast the original frequency-dependent system as a constant one, the price to pay being the creation of additional coordinates.

As explained in the previous section concerning the viscous damping, one can find the complex modes of this second order system. Taking for instance the parameters

m = 1, k 0 = 100, α 1 = 1, ζ 1 = 100 and ω 1 = 100, (II .51)
one can proceed to an eigenvalue analysis that yield the eigenvalues (λ j = i ω j ) and eigenmodes such as

Λ =         -2•10 4 -0.13 + 14.14i -0.13 -14.14i -0.25         , V =   2.5•10 -7 -0.03 -0.99i -0.03 + 0.99i 0.45 1 -0.04 + 2.5•10 -5 i -0.04 -2.5•10 -5 i 0.89   .
(II.52)

One can identify two elastic modes for j = 2 and 3, with complex and conjugate eigenvalues and eigenvectors. The corresponding deformation vector can be decomposed as a product such as, for the mode j = 2, where e i 14.14t is the oscillatory term and e -0.13t the damping term. The two other eigenvalues have only damping terms so the modes are not "dynamic" modes but only instantly damped deformations.

q 2 (ω 2 ) =    -0.03 -0.99i -0.04 + 2.5•10 -5 i    e i 14

II. 2. 2 Two-dimensional spring-mass example

The Figure II.14 depicts a system made of two mass connected through a complex frequency-dependent stiffness modeled by a 1-order GHM model.

m 1 m 2 k(ω) q 1 q 2 f 1 f 2
Figure II.14: Two-dimensions spring-mass system.

The associated matrix system can be written in the frequency domain as

-ω 2 M + k(ω)K q = f , (II .54)
where

M = m 1 0 0 m 2 , K = 1 -1 -1 1 , q = q 1 q 2 and f = f 1 f 2 (II.55)
are the mass matrix, the "factorized" stiffness matrix, the variable vector and the force vector respectively. One might introduce the dissipation coordinates stored this time in the vector z such that

z = z 1 z 2 = γ 1 (ω)q, with γ 1 (ω) = ω 2 1 -ω 2 + i ω2ζ 1 ω 1 + ω 2 1 .
(II.56)

The two relations needed to build the augmented system are

-ω 2 M + 1 + α 1 K 0 q -α 1 K 0 z 1 = f , (II.57a) -ω 2 + i ω2ζ 1 ω 1 + ω 2 1 z -ω 2 1 q = 0, (II.57b)
where K 0 = k 0 K. Following the same procedure as in the 1-dimension example, one could built an augmented system where the matrices are not symmetric. However, this issue can be alleviate in multiplying equation (II.57b) from the left by α 1

ω 2 1 K 0 , yielding -ω 2 M + 1 + α 1 K 0 q -α 1 K 0 z = f , (II.58a) -ω 2 α 1 ω 2 1 + i ω 2ζ 1 ω 1 + α 1 z -α 1 q = 0, (II.58b)
which leads to the second-order augmented matrix system

-ω M + i ω D + K v = f , (II.59)
where

M = M 0 0 α 1 ω 2 1 K 0 , D = 0 0 0 2α 1 ζ 1 ω 1 K 0 , K = (1 + α 1 )K 0 -α 1 K 0 -α 1 K 0 α 1 K 0 , v = q z and f = f 0 .
(II.60)

Here again, the matrices are real, constant and symmetric. Nonetheless, the matrix K 0 could contain rigid body modes which add singularities. To alleviate this issue, one must avoid the damping forces associated to these modes by first computing a spectral decomposition of K 0 and replacing it through its elastic modes only, such that

K 0 ←-VΛV T (II.61)
where Λ contains only the nonzeros eigenvalues of K 0 and V the associated eigenmodes.

Considering this decomposition in equations (II.58a) and (II.58b), the augmented matrices become

M = M 0 0 α 1 ω 2 1 Λ , D = 0 0 0 2α 1 ζ 1 ω 1 Λ , K = (1 + α 1 )K 0 -α 1 VΛ -α 1 ΛV T α 1 Λ (II.62)
with a modified variable vector

v = q V T z . (II.63)
Multiple variable system modeled by a n-order approximation. When the order of the GHM model is higher than one, the modulus takes the form (II. [START_REF] Meyer | Balancing and model reduction for secondorder form linear systems[END_REF]) and the matrices are extend to

M =        M 0 . . . 0 0 α 1 ω 2 1 K 0 0 . . . . . . 0 . . . 0 0 . . . 0 α k ω 2 k K 0        , D =       0 0 . . . 0 0 2α 1 ζ 1 ω 1 K 0 0 . . . . . . 0 . . . 0 0 . . . 0 2α k ζ k ω k K 0       , K =       K 0 (1 + k j=1 α j ) -α 1 K 0 . . . -α k K v -α 1 K 0 α 1 K 0 0 . . . . . . 0 . . . 0 -α k K 0 . . . 0 α k K 0       , f =          f 0 . . . 0          , v =          q z 1 . . . z k          . (II.64)
The same issue concerning the rigid-body modes can appear, so a spectral decomposition of K 0 must be considered, as explained previously. For a system with a single material modeled with a 1-order GHM expression, the number of coordinates is simply doubled. When a n-order GHM model is employed, the dimension of the augmented system is N (1 + n) with N the dimensions of the original system.

Application to real-life systems. Numerical simulations are often based on the Finite Element Method where the global mass and stiffness matrices stems from the assembly of element matrices and the displacement vector q corresponds to the nodal displacements.

The dissipation coordinates are local to each element. Nonetheless, the GHM method can be applied on each element, on a group of elements or directly on the global matrices. When a spectral decomposition is required to avoid damping forces related to the rigidbody modes, it is clear that manipulating element matrices are more suited than applying the method on the global matrices that could reach a million of variables. Where a highorder GHM model is needed to well approximate the modulus, the number of dissipation coordinate is not negligible. Nevertheless, this number can be reduced by truncating the highest modes of K 0 [25]. In some situations, it has been shown that more than half of the modes can be truncated, which considerably lightens the number of extra variables.

II. 2. 3 Conclusions

After having approximation the complex-valued elastic moduli through a particular rehological model, the GHM method allows to recast a frequency-dependent matrix system as a second-order one with larger but constant matrices. The presence of rigid body modes leads to singularities that can be avoided thanks to a spectral decomposition of the stiffness matrix. A problem arises for particular viscoelastic material. In fact, when a high-order GHM model is required, too many dissipation coordinates are incorporated so substantial numerical resources are required to deal with the final matrix system.

II. 3 Dynamic substructuring

This section gathers the most popular dynamic substructuring techniques and shows how they can be used to create a superelement of a windscreen. It first introduces the notion of superelement before presenting the well known procedures of Guyan, Craig-Bampton and MacNeal-Rubin. Then, it is shown how to adapt these methods in case of numerical models with frequency-dependent stiffness matrix.

II. 3. 1 Component Mode Synthesis

Developed in the 60's, the Component Modes Synthesis or modes synthesis aims at splitting a structure into several substructures, also called macro-or superelements, from which reduced order models are built and finally connected together. To illustrate this, let us consider a simple example depicted in Figure II.15 where a global structure is considered in Figure II.15a representing a mesh with 2D finite elements. After defining a "requires substantial numerical resources", 4 substructures can be identified, reduced through special procedures and connected to it as shown in Figure II.15b. This kind of decomposition was initially used to lighten the numerical effort needed to solve the global problem. In fact, instead of directly solving the equations related to the global structure, one can work separately on each substructure in order to reduce at best the dimensions of its associated matrix system. Once the subsystems are reduced, they are denoted as superelements and are connected to the host. One can then solve the global system containing much less variables, and, if needed, recover the original coordinates into each superelement. The areas shared between the host and the superelements, represented in a red color in Figure II.15, are called the interfaces or the boundaries. It contains the only variables that will communicate with other substructures. Several techniques have been developed for the past decades such as the Guyan reduction [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF] that uses only static terms so inertia and damped effects are generally poorly represented, which limits the method in the low-frequency analysis. The first modes synthesis method has been proposed by Hurty in 1965 [START_REF] Hurty | Dynamic analysis of structural systems using component modes[END_REF] and is based on normal modes with fixed interfaces. Later, Craig and Bampton developed in 1968 a similar method [START_REF] Roy | Coupling of Substructures for Dynamic Analyses[END_REF]. In the 70's, MacNeal [START_REF] Macneal | A hybrid method of component mode synthesis[END_REF] and Rubin [START_REF] Rubin | Improved Component-Mode Representation for Structural Dynamic Analysis[END_REF] introduced another formulation based on free interface normal modes. Other procedures have been developed over time, as those attributed to Benfield and Hruda [START_REF] Benfield | Vibration Analysis of Structures by Component Mode Substitution[END_REF] and Herting [START_REF] Herting | A general purpose, multi-stage, component modal synthesis method[END_REF]. See [START_REF] De Klerk | General Framework for Dynamic Substructuring: History, Review, and Classification of Techniques[END_REF] for a complete review of dynamic substructuring.

The efficiency of these methods is directly linked to the number of variables in the interface. If the interface consists in a large part of the structure, the reduction procedure will not be very effective, while if only a few coordinates define the interface, the ratio given by the projection will be substential. These methods also assume that the structure is lightly damped so the normal modes are representative of the real (damped) dynamic behavior. Several methods have been developed to tackle this issue and are reviewed for instance in [START_REF] De Klerk | General Framework for Dynamic Substructuring: History, Review, and Classification of Techniques[END_REF][START_REF] Craig | A review of time-domain and frequency-domain component mode synthesis method[END_REF]. The subsection 3. 5 shows how these methods can be modified in order to take the non negligible damping effects into account. In the style of E. Balmès and A. S. Plouin [3,4,[START_REF] Plouin | Pseudo-modal representations of large models with viscoelastic behavior[END_REF][START_REF] Plouin | Steel/Viscoelastic/Steel Sandwich Shells Computational Methods and Experimental Validations[END_REF], the idea is to enrich the basis made with the normal modes by adding correction or residual terms.

II. 3. 2 Guyan reduction

The Guyan reduction [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF] involves expressing a subset of coordinates as a function of the rest of the coordinates through a static projection. Consider the matrix system

-ω 2 M 11 M 12 M 21 M 22 + i ω D 11 D 12 D 21 D 22 + K 11 K 12 K 21 K 22 q 1 q 2 = f 0 , (II.65)
where the coordinates can be decomposed into two subsets. The first one are associated to the interface and can be subjected to exterior loads f . The second one represent the interior of the structure that does not communicate to other models and are free of any loads. Then, Guyan proposed to consider the stiffness terms of the second line of the equation to express q 2 as a function of q 1 such that

q 2 = -K -1 22 K 21 q 1 .
(II.66)

Finally, one can define a relationship between q 1 and the whole coordinates vector

q 1 q 2 = I -K -1 22 K 21 q 1 = T G q 1 , (II.67)
where the static transformation matrix T G allows projecting the initial matrix system to keep only the coordinates q 1 as variables such as

-ω 2 T -1 G MT G + i ωT -1 G DT G + T -1 G KT G q 1 = T -1 G f 0 . (II.68)
As any substructuring technique, this method is efficient when the number of variables to condensed is large. This transformation is a static projection since T G comes from a static relation (stiffness terms only) between the variables. Then, damping and inertia effects are generally not well taken into account. Other method, as the following ones, contains not only static but also dynamic terms.

II. 3. 3 Craig-Bampton

Certainly the most famous dynamic substructuring technique, the Craig-Bampton procedure [START_REF] Roy | Coupling of Substructures for Dynamic Analyses[END_REF] consists in expressing the interior variables as a contribution of static modes and fixed interface dynamic modes. The static modes contain the value of the interior variables when all the interface variables are considered as fixed, except one of them that is subjected to a unitary displacement. In doing so, one obtains as static modes as interface variables. Let q b and q i be the boundary (interface) and interior variables of the mechanical system

-ω 2 M bb M bi M ib M ii + K bb K bi K ib K ii q b q i = f 0 . (II.69)
The static terms of the second line allow to express the interior dofs q i as a function of q b thanks to the Guyan reduction

q i = -K -1 ii K ib q b , (II.70)
and, following the definition given previously, the static modes are defined by the columns of -K -1

ii K ib such that

Φ C = -K -1 ii K ib . (II.71)
Then, the dynamic modes correspond to the eigenvectors of the system considering the interface variables q b as null, so one must solve the eigenvalue problem

K ii ϕ N = λM ii ϕ N , (II.72)
related to the interior variable only. The eigenvalues are stored in the diagonal matrix Λ N and the dynamic modes are stored in the matrix Φ N , where the subscript N refers to the term normal modes, employed in the section 1. In practice, numerical simulations consider a finite frequency interval, generally 0, f max . A common rule consist in retaining the modes with eigenfrequencies smaller than αf max where α is usually taken as 1.5 or 2.

Once the two sets of modes computed, one expresses q i as a contribution of static and dynamic modes such that

q i = Φ C q b + Φ N δ, (II.73)
where δ is a vectors of modal coefficients. Then, the relationship between the initial and new variables is given by

q b q i = I 0 Φ C Φ N q b δ = T CB q b δ . (II.74)
This transformation can be seen as a augmented Guyan transformation since T CB is T G enriched with dynamic modes. This way, inertia and damping terms can be taken into account. A Graamschmidt procedure is applied to reorthogonalize the basis Φ C Φ N . The efficiency of this method is based on the assumption that a small number of modes can represent the dynamic behavior of the interior of the structure. Then, the large number of interior variables is reduced to a handful of modal coefficients.

II. 3. 4 MacNeal-Rubin

The MacNeal [START_REF] Macneal | A hybrid method of component mode synthesis[END_REF] and Rubin [START_REF] Rubin | Improved Component-Mode Representation for Structural Dynamic Analysis[END_REF] methods differs from the Craig-Bampton procedure by the nature of the modes. Here, a static modes is defined as the deformation of the whole structure due to a unitary load applied on one of the interface variables. These modes are stored in the matrix Φ A and verify

Φ A = K -1 I 0 . (II.75)
After, the dynamic modes contained in the matrix Φ N are computed considering a free interface so that the associated eigenvalue problem concerns the whole variables:

Kϕ = λMϕ. (II.76)
Then, the relationship between the initial and new variables is given by

q b q i = Φ A Φ N f b δ = T MN f b δ , (II.77)
where f b is the load vector related to the interface coordinates and δ the vectors of modal coefficients of the free interface normal modes.

Here, the static modes require solving a system of equations while in the Craig-Bampton procedure they are deduced from a matrix-vector product. Nonetheless, this step must be done only once so this offline cost could be seen as negligible regarding the final reduction ratio.

II. 3. 5 Substructuring with frequency-dependent stiffness matrix

The methods presented previously involve computing normal modes with fixed or free interface. These modes are the eigenvectors of the system with or without boundary conditions. When viscoelastic materials are employed, the stiffness matrix depends on the frequency so the equations of motion are written as

-ω 2 M + K(ω) q = f , (II.78)
or, splitting the stiffness matrix K(ω),

-ω 2 M + E(ω) Kv + K e q = f , (II.79)
where K e is the purely elastic stiffness matrix containing the information related to the constant elastic moduli, and E(ω) Kv = K v (ω) is the viscoelastic matrix related to a single frequency-dependent elastic modulus E(ω). Such a formulation leads to a nonlinear eigenvalue problem K(ω j )ϕ j = ω 2 j Mϕ j , (II.80)

where the eigenvalue ω 2 j and eigenvector ϕ j are complex-valued. This non-linearity makes the classical mode superposition unusable as it is. In practice, and on particular conditions, some techniques have been developed to alleviate this issue. A complete review made by Rouleau et al. [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF] The just quoted paper provides clear explanations on each methods so they will not be recalled in this manuscript. Using the transformation δ = Tq, the reduced system has the form

-ω 2 T H MT + T H K e T + E(ω)T H Kv T δ = T H f , (II.81)
where, even if the matrices are reduced, the user needs to know the value of the elastic modulus E(ω) at each frequency step. This prevents sharing models when confidentiality is the key issue, which is one of the main purposes of this thesis. Clearly, if some of the listed techniques are applied on a numerical model of windscreen, K e will represent the stiffness of the glass sheets and E(ω) Kv the stiffness related to the PVB, where E(ω) is the confidential data. It is shown in the following chapters that a POD-based model reduction method applied to a GHM system can lead to a reduced model containing all the information related to the damping properties of the PVB without the need to evaluate the modulus at each frequency step.

II. 3. 6 Superelement of a windshield

Let us finish this section with the application of reduction methods or substructuring techniques on numerical models of windscreens. Based on the illustration in Figure II.16, one can identify a windscreen connected to a car. The area shared between the two models defines the interface (or boundary), which is the surface where the polyurethane glue is applied. Then, one can split the global variable vectors of displacements q of the windscreen in two subsets : q b that gathers the variables of the interface, and q i that collects the rest of the variables. In order to illustrate a real life application, consider the finite element model of the windscreen assembled on the car "Renault Zoé" shown in Figure II. [START_REF] Chazot | Harmonic response computation of viscoelastic multilayered structures using a ZPST shell element[END_REF]. Three dimensional linear finite elements are used which generates a model with 7 200 elements, 9 788 nodes and 29 364 variables. Since the polyurethane bond is also represented, the interface can be identified and counts 728 variables, making q b ∈ C 728×1 . Then, the rest totalizes 28 636 variables, which gives q i ∈ C 28 636×1 . It is clear that if the 28 636 interior variables are expressed as a sum of a few modes, the reduction would be substantial.

Figure II.17: FE model of the windscreen of the car "Renault Zoé".

II. 3. 7 Conclusions

The principle of dynamic substructuring and the main methods have been introduced. The dynamic behavior of a structure can be expressed as a contribution of static and dynamic modes so the equations are rewritten in order to build the so called superelement. All of the reduced substructures can be assembled, i.e. connected to each others or to a host structure, which forms the global system. This allows much savings in terms of computation times and memory. Moreover, this is also helpful in industrial projects where different teams are involved since they can work on different substructures at the same time.

The dynamic modes are essentially the eigenvectors of the system. When viscoelastic materials are employed, the stiffness terms depend on the frequency so the dynamic modes basis can then no longer be determined as usual. However, for lightly damped structures, a basis made from eigenvectors computing with a constant elastic modulus could well represent the dynamics (MSE method). Nevertheless, such structures as windscreens are made of too high and too frequency-dependent materials so other techniques must be used. Several adaptations exist and provides great results [START_REF] Plouin | Steel/Viscoelastic/Steel Sandwich Shells Computational Methods and Experimental Validations[END_REF]3]. The point is the elastic modulus must be known at each frequency step, even in the reduced order models. When confidentiality is the key issue, this formulation could restrict the exchanges with other engineers.

II. 4 SVD-based model reduction methods

It has been shown in the previous section that a basis made of the eigenvectors of a matrix system could be used to build a reduced order model. When viscoelastic material are employed, the stiffness terms depend on the frequency and finding such a basis is no that simple. This section aims at presenting an other class of model reduction techniques, based on the Singular Value Decomposition (SVD).

This section first presents linear input-output systems following by the Balanced Truncation (BT), a famous SVD-based reduction method. This technique provides very good performance but can only be applied on small to medium size systems. In order to deal with higher dimensions, one may use the Balanced Proper Orthogonal Decomposition (BPOD) that finds an approximation to the reduced system obtained with BT using numerical algorithms adapted to large systems. The BPOD can be seen as Proper Orthogonal Decomposition (POD) respecting to a particular inner product.

The following paragraphs lie on several notions, especially the POD via the method of snapshots explained in Appendix A. The reader can also refer to Appendix B that gives some details on the algorithms used in BPOD.

II. 4. 1 LTI systems in state-space representation

The following subsections introduces the notion of state-space systems and some model reduction methods that could be applied on, as the Balanced Truncation (BT) and The Balanced Proper Orthogonal Decomposition (BPOD). Since this manuscript concerns vibro-acoustic simulations in the frequency domain, these methods will be introduced first in the time domain, and then in the frequency domain.

A Linear Time Invariant (LTI) system can be described in the time domain through the system of equations ẋ

(t) = Ax(t) + Gu(t) y(t) = Lx(t) (II.82)
where A ∈ R N ×N , G ∈ R N ×m and L ∈ R p×N are called the system matrices. The matrix A is called the state space matrix and G, L are called the input and output map respectively. Another matrix, called the direct transmission map could be introduced but is not considered in this work. The vectors x ∈ R N , u ∈ R m and y ∈ R p are called the state, input and output vectors of the system. Also, the first equation is called the state equation and the second is referred to the output equation. When a LTI is composed of only 1 input and 1 output, i.e. m = p = 1, it is called a Single Input Single Output (SISO) system, while when m, p ̸ = 1 it is qualified as a Multiple Input Multiple Output (MIMO) system.

m 1 m 2 k q 1 q 2 f 1 + 2f 2 f 1 -3f 2 + 2f 3 u =    f 1 f 2 f 3    , y =        q 1 q 2 q 1 -q 2 q 1 +q 2 2        Figure II.18:
Example of a MIMO system with 2 variables, 3 inputs and 4 outputs.

An example of a MIMO system is given in Figure II.18 where a spring-mass system is subjected to 3 inputs and 4 outputs. Note that the number of inputs and outputs can be greater than the number of sates (variables). Considering a state vector q = q 1 q 2 T , the input and output maps B and C are

f = Bu = 1 1 0 1 -3 2    f 1 f 2 f 3    and y = Cq =     1 0 0 1 1 -1 1 2 1 2     q 1 q 2 , (II.83)
where the product Bu is equivalent to a classical force vector f . Such systems are asymptotically stable when all the eigenvalues of A have a negative real part so the solution of the system tends to 0 as the time goes to infinity for an impulse response, which is always the case for passive and damped mechanical systems. A more general definition consists in qualifying the matrix A as Hurwitz-stable if its eigenvalues λ are all contained in the open left half of the complex plane, such as Re(λ) ∈ R -. We can also talk about Bounded Input Bounded Output (BIBO) stability. This more intuitive concept assumes that a system is stable if and only if, for a bounded input, the output is also bounded.

II. 4. 2 Second order models and generalized state-space form

Stiffness, damping and inertia terms yield second-order motion equations, which are directly linked to a time invariant second-order system

Mq(t) + D q(t) + Kq(t) = Bu(t) y(t) = Cq(t) , (II.84)
where q(t) is the state vector, M, D, K ∈ R N ×N are the mass, damping and stiffness matrices, B ∈ R N ×m the input map and C ∈ R p×N the output map. When it comes to vibration studies, accelerations could be interesting data to describe the dynamic behavior.

Then, it could be useful to rewrite the output equation with position and velocity output maps, respectively C p and C v , such that

y(t) = C p q(t) + C v q(t). (II.85)
Nevertheless, only position unknowns will be considered in the following so the systems will have the representation of Eq. (II.84). Some mathematical transformations can only be applied to first order systems, so our second order state-space model has to be turned back into a first order form. To do this, one must define a new state vector x(t) = q(t) T qT (t)

T and linearize the system as the first companion form

     N 0 0 M ẋ(t) = 0 N -K -D x(t) + 0 B u(t) y(t) = C 0 x(t) , (II.86)
or in the second companion form

     D M N 0 ẋ(t) = -K 0 0 N x(t) + B 0 u(t) y(t) = C 0 x(t) , (II.87)
where N ∈ R N ×N can be any non-singular matrix [START_REF] Françoise | Meerbergen. The Quadratic Eigenvalue Problem[END_REF] often chosen as the identity matrix, -K or M to make the linear system symmetric. The linearization returns a double size state-space in its generalized form

E ẋ(t) = Ax(t) + Gu(t) y(t) = Lx(t) , (II.88)
where E is called the descriptor matrix. It can be easily turned into a standard state-space by multiplying the state equation by E -1 from the left, which gives

ẋ(t) = E -1 Ax(t) + E -1 Gu(t) y(t) = Lx(t) (II.89) where E -1 A ∈ R 2N ×2N is the standard state matrix and E -1 G ∈ R 2N ×m is the standard input map.
The standard form can also be found by multiplying the second order state equation (II.84) by M -1 from the left side and choosing N = I. However, when M, D and K are sparse, seeking the standard form destroys the sparsity of the block matrices.

Working with the generalized instead of the standard form allows enjoying performances of sparse solvers. On top of that, when M, D and K are symmetric, choosing N = M makes E and A for the second companion form. When C = B T , the whole state-space is symmetric.

II. 4. 3 Transfer function and motivations for model reduction

The transfer function is defined as the relation between the inputs u and the outputs y. In the frequency domain, they are denoted by H 1 (ω) and H 2 (ω) for the first-and second-order forms respectively, such that

H 1 (ω) = L(i ωE -A) -1 G H 2 (ω) = C(-ω 2 M + +i ωD + K) -1 B.
(II.90) These two expressions link the same inputs and outputs such as H 1 (ω) = H 2 (ω) = y(ω)/u(ω). We conclude that it describes the same system, in an input/output meaning only. Applying the transformation x = Tx with T ∈ R 2N ×2N (or R N ×N for the second order form) a non-singular matrix modifies the matrices but does not affect the inputs and outputs vectors. Consequently, the transfer functions remain invariant under state-space transformations, such that

H 1 (ω) = L(i ωE -A) -1 G = LT(i ωT -1 ET -T -1 AT) -1 T -1 G, (II.91)
for the first-order form. A dynamic system can be represented with different state-space models called realizations of this system. With T ∈ R 2N ×r , the order of the new realization is r < 2N , which enables to reduce the size of the matrices without changing the transfer function, as illustrated in Figure II. [START_REF] Roy | Coupling of Substructures for Dynamic Analyses[END_REF]. This is the aim of model order reduction.

y = Ĉ ←-i ω Ê -Â -1 B u
Figure II.19: Illustration of two state-space systems, before and after reduction. Inputs and outputs remain the same while the system matrices are smaller.

In reality, making the matrices smaller modifies the system so the new transfer function is not exactly the same. The user will have to avoid working with r < r lim at risk of excessively degrade the quality of the reduced order model. The realization with r = r lim is called the minimal realization and r lim the McMillan degree of the system.

II. 4. 4 Controllability and observability Gramians

The notions of controllability and observability Gramians, used in control theory, are the basics of the Balanced Truncation (BT) reduction method. Considering the linear and stable generalized state-space

E ẋ(t) = Ax(t) + Gu(t) y(t) = Lx(t) (II.92)
one can define the controllability and observability Gramians P and Q such as

P = ∞ 0 R(t)R T (t)dt, with R(t) = e E -1 At E -1 G (II.93a) Q = ∞ 0 S(t)S T (t)dt, with S(t) = e (E T ) -1 A T t (E T ) -1 L T (II.93b)
The notation (•) T denotes the transpose. These matrices are usually computed by solving the associated generalized Lyapunov equations [62, 49]

APE T + EPA T = -GG T A T QE + E T QA = -L T L (II.94)
that become impossible to solve when it comes to very large models [6]. Two definitions of the Gramians can also be given in order to get a better comprehension of what they represent:

• The controllability Gramian P gives the degree of excitation of each state for a given input. Considering two different states x 1 and

x 2 with ||x 1 || = ||x 2 ||, if x T 1 Px 1 > x T 2
Px 2 , then the state x 1 is considered as more controllable than the state x 2 . It means that the input u needs less energy to drive the system from rest to x 1 than to x 2 .

• The observability Gramian Q gives the degree of contribution of each state for each output. For an initial state x 0 , and with no input, one has ||y|| 2 2 = x T 0 Qx 0 . States which excite larger output signals are denoted by more observable states. It is also possible to define the observability gramian as the controllability gramian of the adjoint system defined by

E T ẋ(t) = A T x(t) + L T u(t) y(t) = G T x(t) , (II.95)
which highlights their duality.

Furthermore, the Gramians depends on the coordinates in which they are expressed.

Moore [START_REF] Moore | Principal component analysis in linear systems: Controllability, observability, and model reduction[END_REF] shown that under a particular change of coordinates x = Tx, the Gramians become diagonal and equal as

P = Q = T -1 P T -1 T = T T QT =    σ 1 . . . σ 2N    , (II.96)
where σ i ⩾ σ i+1 ⩾ 0. The matrix T exists as long as the system is both controllable and observable, which is the case in structural dynamics most of the time. This involves that, in this new space, all the states have simultaneously the same degree of controllability and observability. These degrees are represented by the Hankel singular values (HSV) denoted by σ i such as the more controllable and observable states have a larger HSV than the less controllable and observable states. As for eigendecomposition, this transformation allows the user to keep the state with high HSV and delete those with low HSV. This is the aim of the Balanced Truncation model reduction presented in the following subsection.

II. 4. 5 Balanced Truncation

The Balanced Truncation (BT) technique consists first in rewriting the system in a new basis where the two Gramians are diagonal and equal, as explained in the previous subsection. Then, a truncation step allows to keep only the blocks of the matrices related to the most important states. This method was developed by Moore in 1981 in control theory [START_REF] Moore | Principal component analysis in linear systems: Controllability, observability, and model reduction[END_REF], but has also been used in other fields, such as structural dynamics [START_REF] Friswell | Reduced-Order Models of Structures with Viscoelastic Components[END_REF][START_REF] Benner | Improved Second-Order Balanced Truncation for Symmetric Systems[END_REF][START_REF] Benner | Balanced Truncation for Descriptor Systems with Many Terminals[END_REF].

The transformation matrix, denoted here by T, is generally composed of the appropriately scaled eigenvectors of the product PQ such that PQT = TΣ 2 with Σ = diag(σ 1 , . . . , σ 2N ) the matrix of HSV

P Q = Σ 2 T -1 P T -1 T T T QT = Σ 2 T -1 PQT = Σ 2 PQT = TΣ 2 .
(II.97)

Notice that the eigenvalues of PQ are the squares of the HSV. Once T (and T -1 ) is computed, the truncation is easy to realize since it consists in removing the last modes that correspond to the less controllable and observable states. In practice one must keep the first r columns of T and rows of T -1 .

Another algorithm can be used to find T. Let R and S ∈ R 2N ×2N be the Cholesky factors of the Gramians such as

P = R RT and Q = SS T . Then, compute the Singular Value Decomposition (SVD) of ST R such that ST R = U 1 U 2 Σ 1 Σ 2 V T 1 V T 2 ≈ U 1 Σ 1 V T 1 .
(II.98)

Hence the matrix T and its inverse T -1 appears to be

T = RV 1 Σ -1 2 1 and T -1 = Σ -1 2 1 U T 1 ST . (II.99)
Here, the truncation is done by keeping the first r Hankel singular values in Σ 1 and the associated modes in U 1 and V 1 .

Once T and T -1 found, one can use the transformation x = Tx and write the reduced state-space system Êx(t) = Âx(t) + Ĝu(t)

y(t) = Lx(t) , (II.100) where Ê = T -1 E T ∈ R r×r , Â = T -1 A T ∈ R r×r , Ĝ = T -1 G ∈ R r×m and L = L T ∈ R p×r .
To be considered as a Reduced Order Model (ROM) we must have r < 2N , or even r ≪ 2N for an effective ROM.

There are two major advantages for the balanced truncation method. First the stability of the original system is preserved, and then a priori error bound is available. After defining the two transfer functions H(ω) for the original system and Ĥ(ω) for the ROM such as

H(ω) = L i ωE -A -1 G, (II.101a) Ĥ(ω) = L i ω Ê -Â -1 Ĝ, (II.101b)
the error bound can be shown to be [START_REF] Moore | Principal component analysis in linear systems: Controllability, observability, and model reduction[END_REF]

∥H(ω) -Ĥ(ω)∥ ∞ ⩽ 2 2N i=r+1 σ i , (II.102)
and allow to determine the order of the reduced system a priori by choosing a prescribed error tolerance.

II. 4. 6 Balanced Proper Orthogonal Decomposition

We saw in the previous sections that the balancing transformation can be computed by solving an eigeinproblem or by performing a SVD on ST R where R and S are the Cholesky factors of the Gramians P and Q. The size of these matrices is 2N × 2N so computations become intractable for very large systems. The main idea of Balanced Proper Orthogonal Decomposition (BPOD) [START_REF] Willcox | Balanced Model Reduction via the Proper Orthogonal Decomposition[END_REF] is to compute an approximation to the Balanced realization, i.e. the reduced model coming from BT, from a low-rank approximation to the Gramians. The fact that exact Gramians never need to be computed constitutes a considerable benefit from a numerical point of view. Other techniques tackle the large scale issue and try to approximate the Gramians via special algorithms such as the Alternating Direction Implicit (ADI). It is used to find an approximation of the solution of the Lyapunov equations (II.94) [START_REF] Saak | Efficient Numerical Solution of Large Scale Algebraic Matrix Equations in PDE Controland Model Order Reduction[END_REF][START_REF] Benner | Improved Second-Order Balanced Truncation for Symmetric Systems[END_REF][START_REF] Benner | Balanced Truncation for Descriptor Systems with Many Terminals[END_REF]. Note that that ADI, as the BPOD, requires solving the system for several integration points (times or frequencies according the domain). This task could be numerically intense when systems with many variables are the subject of interest.

As explained in [START_REF] Rowley | Model Reduction For Fluids, Using Balanced Proper Orthogonal Decomposition[END_REF], BPOD can be presented as BT using the method of snapshots. The transformation matrix is found from an SVD, as developed previously, while the factors are not the Cholesky factors of P and Q but low-rank factors such that P ≈ RR T and Q ≈ SS T . It can be shown that BPOD has deep connections with Proper Orthogonal Decomposition (POD) while it is presented as an approximation to the BT. For a theoretical background on the POD via the method of snapshots, please refer to Appendix A.

Time domain formulation.

In order to approximate the Gramians, the integral expressions (IV.10a) and (IV.10b) can be considered on a finite time interval 0, T and turned to quadrature sums on J discrete times

P ≈ J j=1 R(t j )R T (t j )δ j Q ≈ J j=1 S(t j )S T (t j )δ j .
(II.103) where δ j are quadrature coefficients divided by 2π. These sums are the result of a numerical integration following a quadrature rule [START_REF] Davis | Methods of Numerical Integration[END_REF]. Note that different integration schemes could be considered for the two Gramians, which would lead to different integration points t j and weights δ j . For the sake of simplicity, the same rule is considered for P and Q.

The columns of R(t j ) and S(t j ) represent the state of the direct and dual systems, called snapshots, at discrete times t j . These data can be assembled into the matrices Then, one compute the SVD on S T R ∈ R pJ×mJ with mJ, pJ < 2N such as

R ∈ R 2N ×mJ and S ∈ R 2N ×pJ R = . . . R(t j ) δ j . . .
S T R = U 1 U 2 Σ 1 Σ 2 V T 1 V T 2 ≈ U 1 Σ 1 V T 1 . (II.106)
To be effective, the method of snapshots must be applied on state-space with few inputs and outputs such that mJ, pJ ⩽ 2N . This way, the size of S T R is much smaller than ST R and the SVD is performed very quickly. Then, the transformation matrix T and its inverse T -1 appear to be

T = RV 1 Σ -1 2 1 and T -1 = Σ -1 2 1 U T 1 S T . (II.107)
Frequency domain formulation. Since the BPOD is usually introduced in the time domain, this paragraph gives the formulation in the frequency domain. Let us start by rewriting the state-space system (II.92) in the frequency domain such that

i ωEx(ω) = Ax(ω) + Gu(ω) y(ω) = Lx(ω) . (II.108)
Note that the standard linear form (II.82) or the second-order form (II.84) could also be rewritten in a same way. Then, using an extension of Parseval's theorem, the integral expressions of the Gramians (IV.10a) and (IV.10b) given in the time domain can be rewritten in the frequency domain [START_REF] Karim Aghaee | Principle of frequency-domain balanced structure in linear systems and model reduction[END_REF] such that

P = 1 2π ∞ -∞ R(ω)R H (ω)dω, with R(ω) = i ωE -A -1 G, (II.109a) Q = 1 2π ∞ -∞ S(ω)S H (ω)dω, with S(ω) = -i ωE H -A H -1 L H , (II.109b)
where (•) H represents the conjugate transpose. As studies generally consider a bounded frequency interval, there no need to build Gramians valid on a infinite bandwidth. Instead, the Gramians can be defined only for a given frequency range ω a , ω b such as

P = 1 2π ω b ωa R(ω)R H (ω)dω, (II.110a) Q = 1 2π ω b ωa S(ω)S H (ω)dω. (II.110b)
Then, one can compute the snapshots, i.e. computing and storing R(ω j ) and S(ω j ) for J frequencies into the interval ω a , ω b , which lead to the approximation

P ≈ J j=1 R(ω j )R H (ω j )δ j Q ≈ J j=1 S(ω j )S H (ω j )δ j .
(II.111)

Finally, the transformation matrix T and its inverse T -1 are found the same way than using Cholesky factors.

Deep connections with POD. As Rowley elegantly explained in section 3.4 of [START_REF] Rowley | Model Reduction For Fluids, Using Balanced Proper Orthogonal Decomposition[END_REF], "the balanced truncation [using the method of snapshots, i.e. the BPOD] can be viewed as a bi-orthogonal decomposition, instead of the orthogonal decomposition given by POD.

Alternatively, balanced truncation may be viewed as a special case of POD, using a particular dataset (impulse responses), and using the observability Gramian as an inner product. Appendix B goes back over the explanations made by Rowley in [START_REF] Rowley | Model Reduction For Fluids, Using Balanced Proper Orthogonal Decomposition[END_REF], and reformulates it in the frequency domain.

II. 4. 7 Conclusions

This section has introduced some model reduction methods based on a the Singular Value Decomposition (SVD). After rewriting the system equations under a state-space form, one can reduce it using the Balanced Truncation (BT). In doing so, one can work with the same inputs and outputs while dealing with a smaller matrix system. With a good projection, numerical efforts could be substantially lightened.

When it comes to large scale systems, the BT procedure becomes unfeasible. Instead, one can compute an approximation to the BT in using the Balanced Proper Orthogonal Decomposition (BPOD). This mix between POD and snapshots allows to approximate the Gramians and make the BT applicable to large systems. Nevertheless, the BPOD requires solving the system for a set of discrete points (time or frequencies according the domain) since integrals are approximated by sums. Once all of these simulation stored, a SVD is performed on a correlation matrix where the dimensions are directly linked to the number of integration point and the number of inputs-outputs. When these parameters are large, the SVD becomes difficult to perform. Some scientist have worked on dealing with a large number of outputs, which tackle one of the special situations one can found. Nonetheless, one must compute all the numerical simulations for each inputs, which is still numerically intense if the number of input is large. This limit would be rapidly reached if a GHM system would be concerned. In fact, combining a large scale system of a viscoelastic structure, a high order GHM formulation, a large number of input-output and a large number of integration points, the BPOD would in turn become unfeasible.

II. 5 Conclusion and chapter outlines

This chapter has dealt with with four different topics, namely the damping modeling, the Golla-Hughes-McTavish technique, dynamic substructuring and SVD-based model reduction methods. We saw in the introduction that windscreens are made of glass and polymers with highly frequency-dependent damping properties. Due to this particular behavior, classical model reduction methods cannot be applied as they are. Then, adapted variants have been developed but involve evaluating the elastic modulus at each frequency step. For confidentiality aspects, another formulation of reduced equations where the elastic modulus does not appear can be more suited. This thesis tackles this issue and proposes an original procedure. Starting with a mechanical system with a frequencydependent stiffness matrix, a GHM model is found from experimental data from the viscoelastic material. Then, the GHM system is reduced using the BPOD. Finally, the reduced model is used as a superelement that can be connected to other numerical models. One originality of this work lies in the fact than the numerical costs related to the BPOD Section 5 -Conclusion and chapter outlines are substantially lightened, specially the computation of the snapshots. This manuscript is organized into 6 Chapters. Chapter I has introduced this work in first giving the context and motivations. It has also presented the windscreens in order to better understand the main issues it and why glass manufacturers attach importance to it.

Once the aim of the thesis is well defined, chapter II proposes a state of the art complemented by a theoretical background about the main concepts used in this work. As windscreens have damping properties, it has first presented the different ways to take into account damping in numerical models. It then put the emphasis on the Golla-Hughes-McTavish technique that allows to account for the frequency-dependency of an elastic modulus. Next, some recalls about the famous dynamic substructuring techniques have been given. The Craig-Bampton and the MacNeal methods have been detailed to understand the concept of superelement. Some methods allowing to deal with frequencydependent materials have also been presented. Finally, the concepts of state-space system, controllability and observability are explained in order to introduce the model reduction method used in this work, namely the Balanced Proper Orthogonal Decomposition, which has been, contrary to most papers, presented in the frequency domain.

The next chapters, following the state of the art, are mainly inspired by publications submitted and/or accepted during the thesis. The chapter III concerns the application of the Balanced Proper Orthogonal Decomposition on a matrix system with a frequencydependent stiffness matrix modeled with the Golla-Hughes-McTavish technique. The idea is to show that these two methods can be efficiently mixed, and what are the parameters to care about during the procedure. It puts the emphasis on the reduction step, considering the elastic modulus of the viscoelastic material perfectly approximated by the Golla-Hughes-MacTavish method. This way, the reduction of a transfer function of such systems through the Balanced Proper Orthogonal Decomposition is explained.

Chapter IV extends the procedure given in chapter III in order to build a superelement of a structure made of a frequency-dependent material. Here, the reduction steps is the same but the quality of the approximation of the elastic modulus is emphasized. Then, it is shown how to connect a superelement to another numerical model, just like a reduced order system of a windscreen connected to the numerical model of a car.

Then, chapter V shows how to recover a displacement field from a system reduced with the Balanced Proper Orthogonal Decomposition in order for example to estimate the acoustic radiation of the windscreen.

Finally, general conclusions and outlook are given in chapter VI.

Chapter III

The Balanced POD applied to GHM systems

III. 1 Motivations

The current chapter introduces an original procedure aiming at reducing input-output systems modeling mechanical structures made with highly frequency-dependent damped materials. It is mainly constituted of the paper [START_REF] Berthet | The balanced proper orthogonal decomposition applied to a class of frequency-dependent damped structures[END_REF] published during the thesis in the journal Mechanical Systems and Signal Processing. This publication shows how to reduce a matrix system transformed by the GHM method using a model reduction procedure based on the singular value decomposition, namely the Balanced POD. Several details are given on the numerical steps, and show that the GHM formulation is well suited to the POD. In fact, the POD via the method of snapshots requires solving the system for a set of integration points (times or frequencies according to the domain) and storing the results. Unfortunately, the GHM method could lead to very large systems, which also must be recast as a linear form in order to apply the POD. In doing so, the final dimension is such that it brings to very significant numerical costs from a memory and computation point of view, which makes the method unfeasible. Nevertheless, it turns out that the numerical effort can be substantially lightened since the snapshots of the final (and very large) system can be deduced from the solution of the original frequency-dependent FE matrix system, which can be calculated easily with any commercial software. This way, the method becomes applicable again.

After the paper, further works are presented about the formulation of the BPOD for second-order systems. The idea is to project directly the second-order GHM matrix system without the need to recast it as a linear one. After detailing the procedure, numerical examples are given to compare the effectiveness with the original BPOD.

Abstract

A model reduction technique aimed at computing efficiently Frequency Response Functions of damped structures is presented. The frequency-dependent complex moduli are approximated by a mini-oscillators model, known as the Golla-Hughes-MacTavish (GHM) model, which permits to recast the original problem as a more familiar second-order, constant-coefficient system of equations. The matrix system, although much larger, is then treated by application of the Balanced Proper Orthogonal Decomposition (BPOD) which aims at approximating the transfer function matrix, or equivalently the admittance matrix, connecting forces and displacements at a specified set of points of the vibrating structure. All the necessary ingredients of the reduction strategy as well as its efficiency measured in terms of data reduction, accuracy and computational cost are shown. Two illustrative examples of increasing complexity involving a clamped cantilever beam and a realistic windshield are presented. It is shown that the admittance matrix can be approximated by matrices of very small size which computation can be speeded up via diagonalization. It is concluded that the application of BPOD combined with the GHM decomposition of the frequency-dependent algebraic system proves extremely efficient for the modeling of vibrating structures made of different materials, either viscoelastic or purely elastic.

III. 2 Introduction

Viscoelastic materials [START_REF] Salençon | Viscoélasticité pour le calcul des structures[END_REF] are widely used for passive damping treatments to control vibrations of mechanical structures. Sandwich structures with constrained viscoelastic materials are employed for instance in automotive, specially for windscreens. Since this component can easily radiate sounds, it can affect the acoustic comfort of the passengers inside the vehicle. Consequently, vibroacoustic simulations [START_REF] Balmès | Model validation for heavily damped structures. Application to a windshield joint[END_REF][START_REF] Bouayed | Finite element analysis of the dynamic behavior of a laminated windscreen with frequency dependent viscoelastic core[END_REF] are necessary to limit this undesirable behaviour. These simulations are generally based on the use of Finite Element Method (FEM) which often leads to very large algebraic systems and this renders optimisations procedures usually needed at an earlier design stage very prohibitive if not impossible.

Viscoelastic materials are known to be frequency-and temperature-dependent. The master curves for the complex-valued shear modulus can be measured thanks to a Dynamic Mechanical Thermal Analyzer (DMTA) and a parametric mathematical expression can be derived using existing rheological models such as the generalized Maxwell or Kelvin-Voigt models [START_REF] Meyers | Mechanical Behavior of Materials[END_REF]. The frequency dependence of the FE matrix makes classical model reduction methods inoperative and more sophisticated techniques have emerged in the last decades to remedy this [START_REF] Plouin | Steel/Viscoelastic/Steel Sandwich Shells Computational Methods and Experimental Validations[END_REF][START_REF] Plouin | Pseudo-modal representations of large models with viscoelastic behavior[END_REF]. Although this a not a place for a complete study, we can refer to a recent review by Rouleau et al. [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF]. The modal-based reduction techniques discussed in the just quoted paper aims at describing the dynamical behavior of the whole structure and this is beneficial if for instance, one is interested in the identification of specific regions (a priori unknown) where the highest stresses or maximum displacements occur. There is another class of model reduction techniques which consists in approximating the inputoutput behavior, i.e. the transfer functions, of the original problem. These techniques which usually originate from the fields of numerical mathematics and systems and control can be extremely efficient as long as the internal behavior of the structure is of little interest [START_REF] Besselink | A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control[END_REF]. Among them, the Krylov subspace based model order reduction [2] and the Balance Truncation (BT) are probably the most famous and now widely used in a large range of applications in computational mechanics, electrical and control engineering. BT have been developed by Moore [START_REF] Moore | Principal component analysis in linear systems: Controllability, observability, and model reduction[END_REF] in 1981 for the control theory. Moore developed it for first-order state-space models and, later, second-order versions have been implemented [START_REF] Chahlaoui | Second-order balanced truncation[END_REF][START_REF] Reis | Balanced truncation model reduction of second-order systems[END_REF][START_REF] Benner | Improved Second-Order Balanced Truncation for Symmetric Systems[END_REF]. To compute the balanced realization (reduced model computing from the BT), one must compute the controllability and observability Gramians and find a transformation that makes them equal and diagonal so each state have the same degree of controllability and observability. The least controllable and observable transformed states may then be removed without altering the input-output behavior of the original system. Because BT, and also the Krylov subspace methods, have been developed for the treatment of linear or second order state-space systems, i.e. the matrix system must behave linearly or quadratically with frequency, they can not be employed due to the complicated frequency dependence of the original FE matrix involving viscoelastic materials unless some kind of linearization process is employed.

One technique, which has been investigated in the context of FE discretization of poroelastic materials [START_REF] Dazel | An Extension of Complexe Modes for the Resolution of Finite-Element Poroelastic Problems[END_REF] consists in using a Taylor expansion of the matrix coefficients but leads to tedious and rather heavy algebraic manipulations. Other techniques which have been very popular for the treatment of viscoelastic materials are the Golla-Hughes-MacTavish GHM [START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF]) and the Anaelastic Displacement Field (ADF [START_REF] Lesieutre | Finite element modeling of frequencydependent material damping using augmenting thermodynamic fields[END_REF]) methods. The idea is to consider a generic form for the complex-valued modulus which allows to recast the original matrix system into a low order algebraic system (with real-valued coeffi-Section 2 -Introduction cients) augmented with additional coordinates. This then allows classical model reduction by solving complex-valued eigenmodes and eigenvalues taking into account the damping properties of the material as shown by Friswell and Inman [START_REF] Friswell | Reduced-Order Models of Structures with Viscoelastic Components[END_REF]. However, this operation can be quickly computationally expensive. For instance, Vasques et al. [67] show that the 3M ISD112 material needs three sets of parameters in the GHM model, which means that the FE algebraic system associated with the viscoelastic material must be 4 times bigger, furthermore the second-order system must be linearized first [START_REF] Françoise | Meerbergen. The Quadratic Eigenvalue Problem[END_REF] in order to perform the eigendecomposition. This procedure can be avoided using BT if one is interested in computing the transfer functions between a small set of input-output variables. Friswell and Inman [START_REF] Friswell | Reduced-Order Models of Structures with Viscoelastic Components[END_REF] were probably the first authors to apply it on a GHM FE matrix system by finding the balanced realization using the Cholesky factors of the Gramians. More recently, Zhang et al. [START_REF] Zhang | Model reduction methods with application to frequency-dependent viscoelastic finite element model[END_REF] applied the technique to the modeling of bi-dimensional layered structures.

Controllability and observability Gramians are the unique positive solutions of the Lyapunov equations which can be computationally costly and is therefore limited to small or moderate size matrices (say up to a few thousands) [6]. In order to alleviate these limitations, low-rank factors of the Gramians can be found using algorithms based on the Alternating Direction Implicit (ADI) as shown in [START_REF] Saak | Efficient Numerical Solution of Large Scale Algebraic Matrix Equations in PDE Controland Model Order Reduction[END_REF]. In 2002, Willcox and Peraire [START_REF] Willcox | Balanced Model Reduction via the Proper Orthogonal Decomposition[END_REF] proposed another strategy in order to find, at relatively small cost, a good approximation of the balanced realization. Originally developed in the field of fluid dynamics, the authors of the article used the POD (Proper Orthogonal Decomposition) method of snapshots, developed earlier by Sirovich [START_REF] Sirovich | Analysis of turbulent flows by means of the empirical eigenfunctions[END_REF] in 1991, in order to compute the low-rank factors directly from the responses of the system. The method, called Balanced POD (BPOD) was originally devised for the description of time-evolution systems though alternate representation exist in the Laplace and Fourier domain [START_REF] Willcox | Balanced Model Reduction via the Proper Orthogonal Decomposition[END_REF]. It is the purpose of this work to show the applicability of this approach for solving efficiently the dynamical response of a vibrating structure made with viscoelastic materials. One direction of particular interest to us is the dynamic substructuring for laminated structures such as windshields by computing the frequency responses, or equivalently the admittance matrix, connecting forces and displacements at a specified set of points of the structure. To the authors' knowledge, the use of BPOD for the study of vibrational motions has not yet been explored and this paper aims to present all the necessary ingredients of the method as well as its efficiency measured in terms of data reduction, accuracy and computational cost.

The paper is organized as follows. The next section briefly reminds the GHM technique which aims at recasting the original frequency-dependent system of equations into a second order system with constant matrix coefficients. The BPOD method is then presented in section 4 which includes a discussion on the computational aspects of the method and in particular the benefit of the GHM decomposition is highlighted in this context. In the last section, the method presented here is applied to the vibration of a clamped laminated structure of rectangular shape. Numerical results show that Transfer Function Matrices can be constructed efficiently which allows extremely fast computation of frequency response functions.

III. 3 The GHM method

The equations of motion stemming from the finite element formulation for a structure made of a viscoelastic material can be written in the Laplace domain:

s 2 M v + K v (s) q v = F v , (III .1)
where s = i ω is the Laplace variable and ω is the angular frequency. Here, q v denotes the displacement vector, M v the mass matrix, K v (s) the frequency dependent stiffness matrix and F v the force vector. Since the shear modulus can be factored out such that K v (s) = G(s) Kv , Golla, Hughes and MacTavish [START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF] proposed a model to describe it as a series of mini-oscillator terms such that

G(s) = G 0 1 + k j=1 αj s 2 + 2 ζj ωj s s 2 + 2 ζj ωj s + ω2 j , (III.2)
where G 0 is the static modulus, and the k sets of real-valued parameters αj , ζj , ωj are found by curve fitting from experimental data. The number of parameters can be augmented in order to describe the modulus at best. Now, by adding k vectors z j of extra-coordinates called dissipation coordinates as

z j = ω2 j s 2 + 2 ζj ωj s + ω2 j q v , (III.3)
the equation (III.1) can be recast into a more familiar second-order, constant-coefficient form:

s 2 Mv + s Dv + Kv v v = Fv , (III .4)
where

Mv =        M v 0 . . . 0 0 α1 ω2 1 N 0 . . . . . . 0 . . . 0 0 . . . 0 αk ω2 k N        , Dv =        0 0 . . . 0
0 2 α1 ζ1 ω1 N 0 . . . . . . 0 . . . 0 0 . . . 0 2 αk ζk ωk N        , Kv =       K v (1 + k i=1 αj ) -α 1 K v . . . -α k K v -α 1 N α1 N 0 . . . . . . 0 . . . 0 -α k N . . . 0 αk N       , Fv =          F v 0 . . . 0          , v =          q v z 1 . . . z k          , (III.5)
and, to simplify the notation, we put

K v = K v (0) = G 0 Kv .
Here, and for the sake of generality, we intentionally introduce N which is in principle an arbitrary and invertible matrix. One particular choice is to take N = G 0 I. For symmetry reasons it is more judicious to choose N = K T v , though in this case, the presence of rigid body modes renders the matrix singular. These modes can be removed thanks to a spectral decomposition that makes the block matrices of Kv densely populated, which can be penalizing for large size matrices. III. 4 Balanced POD on a GHM model III. 4

. 1 State-space form and transfer function matrix

For the sake of generality, the method is presented for the case of a vibrating structure made of different materials, either viscoelastic or purely elastic. After finite element discretization and assembly of the global system, equations of motion have the general form

s 2 M + K(s) q = F. (III.6)
Let us define a vector of inputs u ∈ C p and a vector of outputs y ∈ C m . Here, inputs correspond to forces and outputs correspond to linear combinations of the states of the system. The input and output maps are defined such that F = Bu and y = Cq. By construction, the associated transfer function matrix is

H(s) = C s 2 M + K(s) -1 B. (III.7)
After applying the GHM method as described earlier, equation (III. 6) is recast as the second-order system

s 2 M + s D + K v = Bu y = Cv , (III.8)
where M, D and K ∈ R N ×N are matrices with real-valued coefficients. Here N stands for the total number of degrees of freedom including the extra-coordinates (III.3) and the new input/output maps are B ∈ R N ×p and C ∈ R m×N . The linearization of the system (III.8) consists in transforming the quadratic expression into a linear one (see for instance [START_REF] Françoise | Meerbergen. The Quadratic Eigenvalue Problem[END_REF]) so it can be rewritten as

sEx = Ax + Gu y = Lx , (III.9)
with, considering the new state vector

x = v T s v T T , E = D M M 0 , A = -K 0 0 M , G = B 0 , L = C 0 , (III.10)
where E ∈ R 2N ×2N is called the descriptor matrix, A ∈ R 2N ×2N the system matrix and G ∈ R 2N ×p , L ∈ R m×2N the input and output maps with 2N the final size of the linear state-space model. Note that the state-space representation is not unique [START_REF] Françoise | Meerbergen. The Quadratic Eigenvalue Problem[END_REF] and the choice made here is motivated by the need to have real-valued and symmetric matrices (A and E).

III. 4. 2 Balanced truncation using the method of snapshots

The aim of model reduction methods is to find a reduced model of order r where r ≪ 2N whilst keeping the input-output behaviour as close as possible to the original one. For the BPOD, we assume that the original system is stable which implies that E is non-singular and the spectrum of the pencil (A, E) is contained in the open left half complex plane, i.e. all the eigenvalues have a negative real part, which is the case here since we are dealing with a dissipative system. The reduced order model is computed from the controllability and observability Gramians respectively P and Q, positive-semidefinite matrices defined in the frequency-domain within the range of [ω a , ω b ] as

P = 1 2π ω b ωa R(ω)R * (ω)dω, (III.11a) Q = 1 2π ω b ωa S(ω)S * (ω)dω, (III.11b) with R(ω) = i ωE -A -1 G and S(ω) = -i ωE * -A * -1 L * respectively the solutions
of the direct and adjoint systems. Here, the notation (•) * denotes the conjugate transpose.

The BT technique consists first in rewriting the system in a new basis where the two Gramians are diagonal and equal as

P = Q = T -1 P T -1 * = T * QT = Σ =    σ 1 . . . σ 2N    (III.12)
with σ j ⩾ σ j+1 ⩾ 0. This involves that all states have simultaneously the same degree of controllability and observability. This means that the less controllable and observable states associated with small Hankel singular values σ j can be discarded. Classically, the transformation matrix T can be obtained in two different ways. The first one requires solving the eigenvalue problem PQT = TΣ 2 . The second way uses the Cholesky decomposition of the Gramians followed by a Singular Value Decomposition (SVD). This can however be computationally intensive due to the large size of the GHM matrix. In both cases, the construction of Gramians is usually accomplished by solving the Lyapunov equations [6]. Here again, this step can be computationally costly and is therefore limited to small or moderate size matrices. The computational burden can be alleviated by using the method of snapshots [START_REF] Sirovich | Analysis of turbulent flows by means of the empirical eigenfunctions[END_REF] allowing finding a low-rank approximation of the Gramians directly from the original definition (IV.10a) and (IV.10b). The factors R and S are constructed by approximating the integrals (IV.10a) and (IV.10b) using interpolatory quadrature rules, i.e.

P ≈ J j=1 R(ω j )R * (ω j )δ j = RR * , Q ≈ J j=1 S(ω j )S * (ω j )δ j = SS * , (III.13)
where R ∈ C 2N ×Jp and S ∈ C 2N ×Jm stacks in columns the so called snapshots for all the inputs and outputs R = . . . R(ω j ) δ j . . . j=1,...,J S = . . . S(ω j ) δ j . . . j=1,...,J .

(III.14)

We should mention in passing that there is a certain degree of freedom to select the integration parameters (ω j , δ j ) and, although the choice could be made in accordance with the complex-valued resonance frequencies of the vibrating structure, this does not necessarily provide the best choice as discussed later.

Once R and S have been computed, a SVD of the correlation matrix S * ER ∈ C mJ×pJ is performed1 2

S * ER = U 1 U 2 Σ 1 Σ 2 V * 1 V * 2 .
(III.15)

After keeping the first r Hankel singular values (HSV), one obtains a low-rank approximation

S * ER ≈ U 1 Σ 1 V * 1 (III.16)
from which a balancing transformation T and its inverse T -1 can be constructed

T = RV 1 Σ -1 2 1 and T -1 = Σ -1 2 1 U * 1 S * . (III.17)
Finally, the reduced system writes

sx = Âx + Ĝu y = Lx , (III.18)
where the hat symbol refers to the transformed state,

x = Tx, Â = T -1 A T ∈ C r×r , Ĝ = T -1 G ∈ C r×m , L = L T ∈ C p×r .
Thus, the size of the reduced model is directly related to the number of singular values retained during the truncation process in equation (III. [START_REF] Chahlaoui | Second-order balanced truncation[END_REF]). Note that, despite the approximation made in (III. [START_REF] Bouayed | Finite element analysis of the dynamic behavior of a laminated windscreen with frequency dependent viscoelastic core[END_REF]), the reduced order system (III.18) still preserves stability of the original model (which would be expected with the Balance Truncation) as explained in section 3.4 of [START_REF] Rowley | Model Reduction For Fluids, Using Balanced Proper Orthogonal Decomposition[END_REF].

III. 4. 3 Notes on computational aspects

The construction of the reduced model relies on two computational steps which are (i) the snapshots computation (III. [START_REF] Bouayed | Finite element analysis of the dynamic behavior of a laminated windscreen with frequency dependent viscoelastic core[END_REF]) and (ii) the singular value decomposition of the correlation matrix (III.15).

Step (i) involves computing simulations (snapshots) from the direct and adjoint systems. The computational burden can be reduced by observing that, in situations where G = L * , then the solution of the adjoint system is the conjugate of the solution of the direct system (recall that A and E are real-valued and symmetric matrices) which means that only solutions of the direct problem needs computing. Furthermore, the step can be substantially lightened by taking advantage of the GHM model via the definition of the extra-coordinates as shown in equation (III.3). Thus, snapshots can be calculated from the original system (III.1) which is much smaller and solved using standard multiple right-hand sides solvers. In fact, any commercial software could be used for this task regardless of the GHM model of the viscoelastic material.

The computational cost of Step (ii) relies strongly on the size of the correlation matrix (III.15) which is given by the product of the number of integration points by the number of inputs (or outputs). Thus the quality of the approximation (III. [START_REF] Bouayed | Finite element analysis of the dynamic behavior of a laminated windscreen with frequency dependent viscoelastic core[END_REF]) is crucial keeping in mind that the number of integration points must be kept as low as possible. This point is not well documented in the literature especially in the frequency domain although elements of discussion can be found in [START_REF] Montier | Balanced Proper Orthogonal Decomposition Applied to Magnetoquasi-Static Problems Through a Stabilization Methodology[END_REF][START_REF] Benner | Balanced Truncation for Descriptor Systems with Many Terminals[END_REF]. In the present work, different integration schemes are investigated to approximate the integrals (IV.10a) and (IV.10b). A first one, referred as pol, is to pick the integration points from the set of complex resonant frequencies, i.e. the poles, of the system which lie in the frequency range of interest. If eigenvalues could, in principle, be extracted from the generalized eigenvalue problem (III.9) (with u = 0), the size of the GHM matrix renders this operation too penalizing and it is preferable to consider the original frequency-dependent system (III.1) to determine the poles using an iterative scheme [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF][START_REF] Friswell | Reduced-Order Models of Structures with Viscoelastic Components[END_REF]. Another way, much simpler, is to consider a regular distribution of points, either linearly (lin) or logarithmically (log). Finally, standard Gaussian quadrature can be employed (gau). In this latter case, the associated weight δ i must be chosen accordingly whereas one simply put δ i = 1 for the other schemes just described. Note that adaptive schemes can also be employed in order to approximate the integrals in a iterative way [START_REF] Benner | Balanced Truncation for Descriptor Systems with Many Terminals[END_REF].

After constructing of the reduced order model, one can accelerate the calculation of the transfer function matrix by computing the eigendecomposition of the matrix  such that Λ = Ŵ *  V, (III. [START_REF] Roy | Coupling of Substructures for Dynamic Analyses[END_REF] where Λ is a diagonal matrix containing the eigenvalues and V and Ŵ the right and left eigenvectors normalized such that Ŵ * V = I. The reduced model has the alternative form

s μ = Λ μ + Ŵ * Ĝu y = L V μ , (III.20)
where μ are the modal coordinates. The transfer function matrix can be written as

Ĥ(s) = L V sI -Λ -1 Ŵ * Ĝ. (III.21)
It is clear that working in the eigenspace lightens the computational effort since the matrix to invert is diagonal while the new input and output maps Ŵ * Ĝ and L V have the same size than the original ones. Such an eigendecomposition is relatively cheap as the matrix size r generally never exceeds a few hundreds at most. All the aspects discussed are illustrated in the next section.

III. 5 Numerical examples

In order to illustrate the theory developed in the previous sections, here are presented two examples. First the method is tested on a simple sandwich beam with a viscoelastic core. Then the method is applied to a windshield to illustrate a real application.

III. 5. 1 Academic case : viscoelastic beam

We consider a clamped laminated structure of rectangular shape as shown in Figure III. 1. This example is representative of a windshield made of two layers of glass separated by a layer of viscoelastic material (here the 3M ISD112). The mechanical properties are reported in Table IV.1 and the associated GHM parameters are taken from [67], see Table III.2 and Figure III. 2. The finite element model is built using 3D quadratic hexaedrons elements with a total of 13 785 variables. The application of the GHM method with k = 3 terms in the model in equation (III.2) necessitates 19 197 dissipation coordinates so the second order system of equations (III.4) counts 32 982 variables which means that the state-space system reaches 65 964 variables. We are interested in computing the transfer functions between 18 inputs and outputs corresponding to the 3 displacements coordinates (u x , u y , u z ) located at 6 nodes of the FEM mesh as shown in Figure III. [START_REF] Bagley | Fractional Calculus -A Different Approach to the Analysis of Viscoelastically Damped Structures[END_REF]. Note that all computations are realized using Matlab, on a CPU i7-8665U, 1.90 GHz with 8 cores and 32 Go RAM.

Parameter

Glass III.2) for 10 Hz to 1000 Hz. Data taken from [67].

The first step consists in identifying the number of complex resonant frequencies which are present in the frequency range of interest, here [10Hz, 3000Hz]. Using an iterative procedure, 12 poles were identified and this gives the first integration scheme pol. The 3 other schemes lin, log and gau presented in section 4 are then generated with the same number of points in order to compare the effectiveness of the distribution only. The CPU cost for the computation of a single snapshot, shown in Table III.3 clearly shows the substantial saving permitted by the GHM model since only the inversion of the original matrix is needed as already discussed in the previous section.

Matrix system R(s), (s = i ω)

Size Time [ In order to assess the performances of the method, the accuracy of the reduced model is measured in terms of the error indicator:

ϵ(ω) = ||H(i ω) -Ĥ(i ω)|| 2 ||H(i ω)|| 2 . (III.22)
In the following study, we are interested in the effect of increasing the number J of integration points. Results are conveniently shown in Fig. III.3b, III.3d and III.3f where the four distributions are compared for different numbers of points. The first observation is that the integration scheme based on the location of poles (pol ) does not outperform the other schemes although it permits to obtain a reasonable estimate for the number of snapshots that need computing. As expected, the logarithmic distribution gives best results at low frequency and quickly deteriorates as frequency increases. Overall, the Gaussian quadrature appears to be the best choice as it guarantees a stable level of error over the whole frequency range. For each of these tests, the corresponding distribution of singular values (HSV) are shown in Fig. III.3a, III.3c and III.3e. All curves show similar behaviours and a plateau is reached at a certain threshold (nearly identical regardless of the integration scheme) above which singular values are too small to be of any interest and could therefore be discarded. In fact, a closer analysis reveals that the threshold corresponds to a condition number of around 10 -16 which is precisely what would be expected with standard double precision arithmetic. Finally we can notice that the threshold increases slightly with the number of quadrature points and this is a little more pronounced for the Gaussian quadrature. In order to see the effect of truncation in the low rank approximation (III. [START_REF] Chahlaoui | Second-order balanced truncation[END_REF]), four reduced order models are constructed using Gauss-Legendre quadrature with 23 points by keeping the first 10, 25, 75 and 212 Hankel singular values (the latter case corresponds to the threshold) as shown in Figure III. 4. As expected, truncating at threshold maintains a very high accuracy similar to the full model since all discarded values are negligible. Results show that retaining as little as 10 HSV is insufficient to provide acceptable accuracy while with 25 values, the error barely exceeds 1%. By increasing the number of singular values, errors decrease rapidly though faster convergence is reported in the field of Fluid Dynamics [START_REF] Rowley | Model Reduction For Fluids, Using Balanced Proper Orthogonal Decomposition[END_REF]. One reason for this probably lies in the resonant behaviour of the dynamic system, though dissipative, which is typical of vibrational structures. Now, it is instructive to measure the effect of the truncation on a specific frequency response (FRF).

To this end, we consider the two scenarios as depicted in Figure III.1: in the first one, we compute the out-of plane displacement at an edge of the cantilever ("Out" as indicated) due to a vertical and unitary force f 1 , in the second scenario we consider another force f 2 = -f 1 located at the opposite side so that torsional waves can be excited. Results which are all gathered in Figures III.5a, III.5b, III.5c and III.5d permit to identify the link between the global error indicator (IV.21) with the FRF. Results have nearly converged with r = 25 HSV and retaining more than 75 singular values would be superfluous from an engineering point of view. In Table III.4 are displayed the computational times required to calculate a sequence of FRFs (with 1000 frequencies). Clearly, depending on the formulation used, the gain can be very substantial showing a ratio of a few hundred thousands. The factorization of the transfer function matrix in its diagonal form allows fast computation of the FRF with a computational complexity that grows only mildly with the number of HSV retained. ) compared to the original system (III.6) ( ).

III. 5. 2 Application to a windshield FE model

In order to illustrate the performances of the method on a large size damped structure, a windshield FE model with free boundary conditions is investigated (see Figure III.6). The windshield is composed of two layers of glass material of about 1.2 m 2 , separated by a viscoelastic core. The structure is excited by a singular force applied at 2 nodes of the FE mesh (highlighted in red), while the displacement vector is computed at a specific node (highlighted in blue), thus 6 inputs and 3 outputs are considered in this illustrative example. The initial frequency-dependent model consists in 29 364 variables. The viscoelastic material is modeled with the GHM model taken from the previous example (see Table III.2 and Figure III.2), and this gives rise to 14 682 additional coordinates and an original system with 88 092 variables. A Gauss-Legendre quadrature rule with 50 points is employed to generate a reduced model in the frequency range [10, 1000 Hz]. The Table III.5 gathers the times related to the whole reduction procedure including the computation of the snapshots, the Singular Value Decomposition algorithm and so on. All these tasks are performed in 208 seconds (note parallelized computational procedures are used here) which is relatively cheap keeping in mind that this is performed only once. In Table III.6 are shown CPU times for the full model (III.6) as well as for three reduced models generated with 100, 50 and 25 singular values. Here again, reduced models remain much faster than the original one. Moreover, the offline computation times related to their constructions remain inferior to the time required to compute the FRF. Table III.6: FRF CPU times [s] for the windshield study with u = [1 0 0 1 1 0] T . Comparison between the original system (III.6) with dimensions 29 364 and the reduced models on 1000 frequency points.

The numerical precision can be checked in Figure III.7 that shows the relative error ϵ for the transfer functions matrix. As expected, the same level of error is observed for the specific frequency response (here the x-displacement) and for the input vector u = [1 0 0 1 1 0] T in Figure III.8.

III. 6 Conclusion

In this paper, we developed a model reduction technique in order to compute very efficiently Frequency Response Functions of vibrating structures with frequency-dependent damping due to the presence of viscoelastic materials. The method applies to the large class of materials for which the complex modulus master curves fit well with the minioscillators model given by equation (III.3). This specific form allows to recast the original problem into a more familiar second-order, constant-coefficient system of equations. After linearization, a state-space representation involving real-valued symmetric matrices is obtained. This representation is amenable to Balanced Truncation based model reduction techniques which consists in approximating the input-output behavior, i.e. the transfer functions matrix, between specific quantities of interest. The technique relies on the concept of controllability and observability Gramians which are difficult to compute and hardly applicable to large size systems. Following earlier work developed in the field of fluid dynamics, low-rank factors of the Gramians can be obtained directly from the responses of the direct and adjoint systems. One of the main originality of the present work lies in the way in which low-rank factors are calculated. One of the key ingredients is to utilize the mini-oscillator model equation ( III.3) which permits to restrict the calculation of the response from the inversion of the original frequency-dependent system (III.6). This is particularly advantageous as any commercial software could in principle be used for the task. The number of inversion needed, i.e. the snapshots, can be minimized by using an adequate distribution of integration points in equation (III.13). It is shown that choosing these points from the set of complex resonant frequencies, i.e. the poles, of the system which lie in the frequency range of interest does not necessarily provide best results, although it permits to identify a reliable estimate of the number of snapshots required to build the reduced model. A comparison of different integration schemes is carried out and it is observed that the classical Gaussian quadrature rule offers the best compromise. The method is applied to the construction of a transfer functions matrix between 18 inputs (forces) and outputs (displacements) for a clamped laminated structure of rectangular shape. The latter is modeled using Finite Element discretization with 13 785 variables. In this specific example, it is shown that the reduced model permits the calculation of the FRF very efficiently showing computational times reduced by 5 to 6 order of magnitudes. Such gain and efficiency is confirmed in the last section of the paper where the method is applied to the vibrational response of a realistic windshield.

In principle, the range of applicability of the method is only limited by available computer resources, recalling that the computational load mainly stems from the computation of snapshots and the SVD decomposition of the correlation matrix which size grows with the number of snapshots and the number of inputs/outputs in the model. Other aspects of the method, which are not dealt with in the present paper, such as the possibility in retrieving important data, such as displacements, stress and strain anywhere in the structure, remain to be explored.

One great advantage of this reduction technique is the explicit knowledge of the reduced system of the first order equation ( III.18) which permits an efficient dynamic substructuring whereby the reduced system is integrated in a host numerical model. One direction of particular interest for the authors is the development of super-elements for multi-layered structures and its utilization for vibro-acoustic predictions.

III. 7 Further work : Second-order BPOD

The previous paper shows how to reduce a GHM input-output system in its linear form. In some situations, it could be useful to preserve the original structure of the equations. In our case, it consists in directly projecting the GHM matrix system such that

-ω 2 M + i ω D + K v = Bu y = Ĉv , (III.23)
where v ∈ C r the new vector of coordinates, M, D and K ∈ C r×r the reduced mass, damping and stiffness matrices and B ∈ C r×m and Ĉ ∈ C m×r the reduced input and output maps respectively. Preserving the original structure allows a meaningful physical interpretation and generally provides more accurate approximations. Some researchers have adapted the Balanced Truncation (BT) to directly reduce the second-order system without the need to linearize it [START_REF] Reis | Balanced truncation model reduction of second-order systems[END_REF][START_REF] Benner | Improved Second-Order Balanced Truncation for Symmetric Systems[END_REF][START_REF] Chahlaoui | Second-order balanced truncation[END_REF]. This leads to the so called Second-Order Balanced Truncation (SO-BT). The Second-Order Balanced POD (SO-BPOD) presented in the following is largely based on the SO-BT. However, this section aims at giving a clear formulation in the frequency domain, which is not, to the author knowledge, available in the literature.

III. 7. 1 Numerical procedure.

The initial starting point is the particular structure of the state-vector. In fact, motivated by the need to work with a first-order system, the linearization process involves working with a new state vector x, twice bigger, with the following structure:

x = x p x v = x p i ωx p , (III.24)
where x p and x v are the position and velocity parts respectively. Knowing the block structure of the matrices, the position part correspond to the solution of the second-order system such that x p = v. This decomposition remains valid on the whole snapshots matrices such as

R = R p R v , S = S p S v . (III.25)
Then, the Gramians reveal their block structures .26) where P p , P v , Q p and Q v are called the position and velocity Gramians. They can also be decomposed as products of low-rank factors such as

P = P p P 1,2 P H 1,2 P v , Q = Q p Q 1,2 Q H 1,2 Q v , ( III
P p = R p R H p and P v = R v R H v
for the controllability, and Q p = S p S H p and Q v = S v S H v for the observability. A secondorder version of the Balanced POD can then be formed and consists in approximating the eigenvectors of P α Q β instead of PQ, with α, β ∈ {p, v}. This leads to four possible truncated SVD of rank r

S H β R α ≈ U αβ Σ αβ V H αβ . (III.27)
The associated four pairs of projection matrices are

Φ αβ = R α V αβ Σ -1 2 αβ and Ψ H αβ = Σ 1 2 αβ U H αβ S β (III.

III. 7. 2 Stability issues.

For BT, it has been proved that the balanced realization of a stable linear system is also stable. For the second-order forms, different algorithms have been proposed in the literature. It can be shown that for a symmetric and stable second-order state-space system, the SOBTpv [START_REF] Reis | Balanced truncation model reduction of second-order systems[END_REF] and the SOBTf v [START_REF] Meyer | Balancing and model reduction for secondorder form linear systems[END_REF] return a reduced systems also stable. However, for general systems, the methods are not stability-preserving.

The same question arises when using BPOD. For linear systems, the use of a special inner product makes the Galerkin projection of a stable system stability-preserving [START_REF] Rowley | Model Reduction For Fluids, Using Balanced Proper Orthogonal Decomposition[END_REF]. This assumption is valid when the original controllability Gramian Q is used as an inner product. For second-order forms, the SO-BPOD involves using Q p or Q v instead of Q such as the four sets

⟨R p , R p ⟩ Qp , ⟨R v , R v ⟩ Qp , ⟨R p , R p ⟩ Qv , and ⟨R v , R v ⟩ Qv (III.30)
are considered instead of ⟨R, R⟩ Q . Even though the question has not been studied here, it would be interesting to know if the SO-BPOD is also stability-preserving.

III. 7. 3 Numerical example.

Considering the second example of the previous paper, we propose to compare the original BPOD with the SO-BPOD in its four possible forms, i. In this case, the second-order form of the BPOD leads to better results than classical BPOD. One can see first that the singular values in Figure III.9a have a relatively faster convergence with the SO-BPOD. In fact, the threshold is reached after 100 HSV for the BPOD and around 90 for the SO-BPOD pp. As shown in Figure III.10, the four SO-BPOD algorithms provide globally better levels of error than the BPOD, especially the SO-BPOD pp. This algorithm can then built more accurate reduced models than BPOD for a same number of singular values kept. The drawbacks are the stability could be not preserved, and there is no prior way to know which second-order Gramians pair will give the better reduced model.

III. 8 Conclusions

In this chapter, an original model reduction procedure has been proposed. It allows to reduce a matrix system modeled by the Golla-Hughes-McTavish (GHM) technique thanks to a SVD-based algorithm, namely the Balanced Proper orthogonal Decomposition (BPOD). Starting with a linear input-output system with frequency-dependent stiffness terms, one obtains an accurate reduced order model with smaller and constant matrices.

Contrary to modal-based model reduction methods, the proposed procedure does not suffers from the large number of additional coordinates introduced by the GHM technique which is generally its major drawback. Then, the user can choose any high-order GHM model which well fits the experimental data of the viscoelastic material, before building the linear reduced model. Further works on the reduction of second-order systems have been conducted. Based on the Second-Order Balanced Truncation (SO-BT), the Second-Order Balanced Proper Orthogonal Decomposition (SO-BPOD) has been formulated and applied to the numerical model of a windshield. Four different algorithms are available to directly project the second-order GHM matrix system without the need to linearize it. The reduced model regains then its original form. Chapter IV

From BPOD-GHM to superelements

The Balanced POD in the frequency-domain

IV. 1 Motivations

The previous chapter has introduced an original methodology to build Reduced Order Models (ROMs) of mechanical systems containing viscoelastic materials. The idea was to apply the Balanced Proper Orthogonal Decomposition (BPOD) to a matrix system coming from the Golla-Hughes-McTavish (GHM) method. It has been shown that great results can be obtained and the numerical costs related the construction of the ROM can be substantially lightened.

This chapter aims at using the BPOD as a dynamic substructuring technique in order to create superelements of structures made with viscoelastic materials. The following paragraphs are mainly constituted of a preprint written at the end of the thesis. It starts by recalling the basics of the GHM and the BPOD methods and goes on the construction and the use of a superelement. Two examples are given in order to illustrate the performance of the procedure.

One of the main edges of this substructuring method is the form of the final matrix system. The presence of frequency-dependent stiffness involves using adapted modalbased methods [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF] that are used to build superelements [3]. The main procedures require computing at least the static modes and one basis of dynamic modes. We show here that our method can lead to ROMs with dimensions smaller than the number of interface variables. Smaller dimensions naturally lead to faster computations.

The second advantage is the formulation of the final matrix system. The adapted modal-based techniques yield frequency-dependent reduced stiffness matrix, which avoid sharing superelements if viscoelastic material properties are confidential. In our case, the matrices are constant and no longer require to evaluate the elastic modulus at each frequency step.

Dynamic substructuring for mechanical systems with frequency-dependent materials using a POD-based model reduction method

Alexandre Berthet, Emmanuel Perrey-Debain, Jean-Daniel Chazot and Sylvain Germès.

Preprint submitted to Journal of Sound and Vibration, 2023

Abstract

Based on the Balanced Truncation approach, a novel methodology for the construction of a super-element for the dynamic analysis of elastic structures made with viscoelastic materials is presented. Contrary to classical modal reduction techniques based on the Craig-Bampton or McNeal methods where the normal modes basis must be enriched to account for damping effects, the methodology presented here takes advantage of the Golla-Hughes-McTavish rheological model (GHM) before reducing the system via the Balanced Proper Orthogonal Decomposition (BPOD). Following the classical procedure used in dynamic substructuring techniques, Lagrange multipliers are employed to couple the reduced system to a host structure, which consists here of an elastic plate. Numerical examples show the efficiency of the super-element both in terms of computational time and data reduction. The specific features of the method, as opposed to usual modal-based reduction methods are that the construction of the reduced order model is relatively straightforward and can deliver very accurate results irrespective of the frequency behaviour of the elastic modulus of the viscoelastic layer. Furthermore it has the advantage to preserve confidentiality, i.e. the frequency dependence of the modulus of elasticity is not explicitly known.

IV. 2 Introduction

Dynamic substructuring is now a well established and powerful tool for engineers used to analyse the dynamic response of mechanical systems. The technique is particularly favored for the simulation of mechanical vibrations allowing different aspects (acoustics, fatigue/durability and comfort) to be studied. It consists in splitting the global structure into several substructures from which reduced models are built. In Finite Element software, these reduced models gave rise to the superelements [START_REF][END_REF]. The main advantages are (these classifications are taken from the very cited review paper on the subject [START_REF] De Klerk | General Framework for Dynamic Substructuring: History, Review, and Classification of Techniques[END_REF]): 1. it allows to tackle large and complex mechanical systems, 2. the algebraic set of equations associated with each substructure can be optimised and reduced, 3. the possibility of combining both modeled and experimentally identified components, 4. it allows sharing substructures from different project groups and allows companies to work on secure projects.

The basis of substructuring techniques comes from the work of Hurty [START_REF] Hurty | Dynamic analysis of structural systems using component modes[END_REF] in 1965. The Craig-Bampton [START_REF] Roy | Coupling of Substructures for Dynamic Analyses[END_REF] and MacNeal [START_REF] Macneal | A hybrid method of component mode synthesis[END_REF] methods are the two most famous procedures usually employed. The dynamics is represented through fixed and free interface modes respectively, which define the two main classes. Several variants have been developed, notably the dual Craig-Bampton [START_REF] Daniel | A dual Craig-Bampton method for dynamic substructuring[END_REF] or the mixed Craig-Bampton [68]. These methods fall in the wider class of Component-mode synthesis (CMS) methods which became particularly popular among the engineering community as a reduction method for finite element models. It is generally based on the modal analysis of the substructure, allowing to reduce the equations of motion from the physical to the truncated modal domain. Consequently, difficulties arise when structures made of viscoelastic materials [START_REF] Salençon | Viscoélasticité pour le calcul des structures[END_REF][START_REF] Nashif | Vibration Damping[END_REF] are considered since usual modal reduction algorithms are inapplicable as they are because of the frequency-dependence of the stiffness matrix. In order to tackle this, several alternative methods, gathered in a recent review [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF], have been developed in the last decades. The modal-based reduction techniques discussed in the just quoted paper have in common that they all relies on the concept of modal projection. Because of the non-linear nature of the eigenvalue problem (with respect to frequency), the computation of eigenmodes associated with more classical linear eigenvalue problems for a specific set of frequencies is in general favored. The benefit of each method, as shown in [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF], is measured both in terms of accuracy and computational cost. It is clear that the latter is expected to increase with the size of the finite element model and also depends strongly on the frequency behavior of the stiffness matrix. While these techniques aims at describing the dynamical behavior of the whole structure and this is beneficial if for instance, one is interested in the identification of specific regions, there is another class of model reduction techniques which consists in approximating the input-output behavior, i.e. the transfer functions, of the original problem. These techniques which usually originate from the fields of numerical mathematics and systems and control can be extremely efficient as long as the internal behavior of the structure is of little interest. Among them, the Balanced Truncation (BT) [START_REF] Moore | Principal component analysis in linear systems: Controllability, observability, and model reduction[END_REF] is now widely used in a large range of applications in computational mechanics [START_REF] Besselink | A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control[END_REF]. Recently, the authors of the present article have developed an original numerical method, based on the BT approach, aimed at computing efficiently Frequency Response Functions connecting forces and displacements at a specified set of points of a vibrating structure made with elastic and visco-elastic materials. The key ingredients of the method are (i) the approximation of the frequency-dependent complex moduli via a mini-oscillators model, known as the Golla-Hughes-McTavish (GHM) rheological model [START_REF] Golla | Dynamics of Viscoelastic Structures-A Time-Domain, Finite Element Formulation[END_REF][START_REF] Mctavish | The mini-oscillator technique : A finite element method for the modelling of linear viscoelastic structures[END_REF][START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF], which permits to recast the original problem as a more familiar second-order, constant-coefficient system of equations, (ii) the application of the Balanced Proper Orthogonal Decomposition (BPOD) [START_REF] Willcox | Balanced Model Reduction via the Proper Orthogonal Decomposition[END_REF] which aims at approximating the transfer function matrix.

It was shown that the method provides great performances, both in terms of data reduction and computational time, especially when a small set of input-output variables is considered. The aim of the present paper is to develop further the method in order to provide an efficient dynamic substructuring tool (or superelement) which can be applied to a large class of frequency-dependent damped vibrating structures. Besides the numerical performances, we note in passing that the technique has the advantage to preserve confidentiality, i.e. the frequency dependence of the modulus of elasticity is not explicitly known, in contrast with modal reduction techniques discussed above.

This paper is organized as follows. First, Sections 3 and 4 briefly recall the main theoretical ingredients which are at the heart of the GHM and BPOD techniques. In Section 5, the whole damped structure is modeled as a superelement whereby only physical variables (forces and displacements) on the interface with the host structure are considered in the reduced model. Finally, two numerical examples of increasing complexity are given in Section 6.

IV. 3 The Golla-Hughes-McTavish (GHM) approach

Following classical finite element discretisation, the equations of motion of a structure exclusively made with viscoelastic material can be written in the Laplace domain as

s 2 M v + K v (s) q v = F v , (IV.1)
where s = i ω is the Laplace variable and ω is the angular frequency. Here, q v denotes the displacement vector, M v the mass matrix, K v (s) the frequency dependent stiffness matrix and F v the force vector. Since the shear modulus can be factored out such that K v (s) = G(s)K v , Golla, Hughes and McTavish [START_REF] Golla | Dynamics of Viscoelastic Structures-A Time-Domain, Finite Element Formulation[END_REF][START_REF] Mctavish | The mini-oscillator technique : A finite element method for the modelling of linear viscoelastic structures[END_REF][START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF] introduced the so-called GHM model to describe it as a series of mini-oscillators as follows

G(s) = G 0 1 + k j=1 αj s 2 + 2 ζj ωj s s 2 + 2 ζj ωj s + ω2 j , (IV.2)
where G 0 is the static modulus, and the k sets of real-valued parameters (α j , ζj , ωj ) are found by curve fitting to experimental data. One can generate k vectors z j of extracoordinates called dissipation coordinates, all related to the displacement vector q v as z j = γ j q v , with γ j = ω2

j s 2 + 2 ζj ωj s + ω2 j . (IV.3)
And this yields the equivalent second order matrix system

s 2 Mv + s Dv + Kv v v = Fv , (IV.4)
with larger but real-valued and symmetric matrices (see the explicit form in [START_REF] Berthet | The balanced proper orthogonal decomposition applied to a class of frequency-dependent damped structures[END_REF]). Here, the vector v v = [q v z 1 . . . z k ] T contains both physical variables q v and the dissipation coordinates vectors z j .

IV. 4 The Balanced POD in the frequency-domain IV. 4

. 1 State-space representation

For the sake of generality, we consider a structure of arbitrary shape made of different materials, either elastic or viscoelastic. The equations of motion can be expressed, after discretisation with the finite element method, in the general form

-ω 2 M + K(ω) q = Bu y = Cq (IV.5)
where u ∈ C p and y ∈ C m are the input and output vectors whereas B ∈ R N ×p and C ∈ R m×N are the input and output maps. Note that the stiffness matrix contains both elastic and viscoelastic terms. After applying the GHM method, equation (IV.5) is recast as

-ω 2 M + i ω D + K v = Bu y = Cv . (IV .6)
In order to apply the BPOD technique, Eq. (IV.6) must be linearized [START_REF] Françoise | Meerbergen. The Quadratic Eigenvalue Problem[END_REF] to give

i ωEx = Ax + Gu y = Lx , (IV.7)
where matrices are now twice bigger with

E = D M M 0 , A = -K 0 0 M , G = B 0 , L = C 0 , (IV.8)
and x ∈ C 2N is the new state vector

x = v i ωv .
(IV.9)

Note that the state-space representation is not unique [START_REF] Françoise | Meerbergen. The Quadratic Eigenvalue Problem[END_REF] and the choice made here is motivated by the need to preserve symmetry and the sparsity of the block matrices, as opposed to the original work of Friswell [START_REF] Friswell | Reduced-Order Models of Structures with Viscoelastic Components[END_REF] where a standard linear state-space, i.e. the descriptor E is the identity matrix, was desired. Balanced Truncation and Balanced POD both rely on the concept of controllability and observability Gramians. These matrices allow measuring to what degree each state is excited by an input and to what degree each state excites future output respectively. They consist in two real, symmetric and positive semi-definite matrices, denoted by P and Q respectively. Using an extension of the Parseval's theorem [START_REF] Karim Aghaee | Principle of frequency-domain balanced structure in linear systems and model reduction[END_REF], it is shown that for a band-limited signal, i.e. ω a ≤ |ω| ≤ ω b , Gramians are obtained via the integrals1

P = 1 2π ω b ωa R(ω)R H (ω)dω, (IV.10a) Q = 1 2π ω b ωa S(ω)S H (ω)dω, (IV.10b)
where (•) H denotes the complex conjugate transpose, R(ω) = (i ωE -A) -1 G and S(ω) = (-i ωE -A) -1 L T are the snapshots matrices (in the frequency domain) for the direct and adjoint systems. As shown by the authors in [START_REF] Berthet | The balanced proper orthogonal decomposition applied to a class of frequency-dependent damped structures[END_REF], integrals can be fairly approximated using a standard Gauss-Legendre quadrature rule with Gauss points ω j and weights δ j giving

P ≈ 1 2π J j=1 R(ω j ) R H (ω j )δ j = RR H , (IV.11a) Q ≈ 1 2π J j=1 S(ω j ) S H (ω j )δ j = SS H , (IV.11b) where R = . . . (i ω j E -A) -1 G δ j /2π . . . , S = . . . (-i ω j E -A) -1 L T δ j /2π . . .
(IV.12) are called the snapshots matrices which store successive solutions of the direct and adjoint systems for a set of angular frequencies ω j .

IV. 4. 2 Balanced POD : an approximation to balanced truncation

As for any model reduction methods, the aim of the Balanced POD is to find a transformation x = Φx of rank r ≪ 2N that minimizes the error between the original and reduced spaces. Once the low-rank factors R and S of the Gramians are computed, one can apply the Singular Value Decomposition (SVD) on the correlation matrix Z = S H R and keep the first r components such that the following approximation can be made

Z = S H R = U 1 U 2 Σ 1 Σ 2 V H 1 V H 2 ≈ U 1 Σ 1 V H 1 , (IV.13)
where Σ 1 = diag(σ 1 , ..., σ r ) contains the r largest singular values and

U 1 , V 1 ∈ C 2N ×r .
The right and left projections of the Balanced POD Φ and Ψ are then given by IV. 4. 3 Notes on computational aspects (i). As shown in [START_REF] Berthet | The balanced proper orthogonal decomposition applied to a class of frequency-dependent damped structures[END_REF], performing the BPOD in the frequency-domain can yield very accurate results with only a few integration points. This is important since the number of points has a direct impact on the number of columns of the low-rank factors which controls the computational cost for the Singular Value Decomposition in eq. (IV.13). Nonetheless, when the number of inputs and outputs becomes large, the size of the correlation matrix increases accordingly and this may render the computational burden prohibitive. This encourages to use adapted techniques such as the randomized algorithm for low-rank matrix approximation developed by Yu et al. [START_REF] Yu | Efficient Randomized Algorithms for the Fixed-Precision Low-Rank Matrix Approximation[END_REF]. It consists in computing the SVD of a sketch (low-rank approximation) instead of dealing with the original matrix. In doing so, the numerical effort is lightened and only the first singular values and vectors corresponding to the main features of the original matrix are returned. This provides a way to automatically determine the number of modes to keep in the transformation (even though the user can select less modes as indicated).

Φ = RV 1 Σ -1 2 1 , Ψ H = Σ -1 2 1 U H 1 S H . (IV.
(ii). Numerical costs associated with the computation of snapshots can be substantially lightened by taking advantage of the particular structure of the state vector x. In Appendix A, it is shown that the solution of the large matrix system (IV.7) can be deduced from the solution of the original system (IV.5) which can be solved using any commercial software. Moreover, the particular structure of the correlation matrix Z = S H R can also be exploited as explained in Appendix B.

IV. 5 The BPOD as a dynamic substructuring technique

We propose in this section to use the BPOD as a dynamic substructuring technique applied to structures containing frequency-dependent viscoelastic materials. In the following, two systems of equations are considered labeled with indices 1 for the host, and 2 for the super-element. The first set of equations consists in a typical second-order system related to a mechanical structure with (real-valued) mass, damping and stiffness matrices and

-ω 2 M 1 + i ωD 1 + K 1 q 1 = B 1 u 1 + F 1 y 1 = B T 1 q 1 . (IV.16)
The second set of equations has the form of the reduced system equation (IV.15) (input and output maps are equal here):

i ω Ê2 -Â2 x2 = Ĝ2 u 2 y 2 = ĜH 2 x2
.

(IV.17)

The two structures share a common area, often called the interface, on which the two classical continuity conditions are satisfied. First, the displacements are equal, i.e. y 1 = y 2 , or in another form

B T 1 q 1 -ĜH 2 x2 = 0. (IV.18)
The sum of the efforts on the interface must be zero, i.e.

u 1 = -u 2 = u. (IV.19)
Once assembled, the whole global system has the following form

  -ω 2    M 1 0 0 0 0 0 0 0 0    + i ω    D 1 0 0 0 Ê2 0 0 0 0    +    K 1 0 B 1 0 -Â2 -Ĝ2 B T 1 -ĜH 2 0            q 1 x2 u      =      F 1 0 0      . (IV.20)
At this point, two important remarks should be made. First, in the MacNeal or Craig-Bampton methods, the projection basis is made of dynamic and static modes. As discussed in [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF] this basis must be enriched to correctly represent the damping effects and the accuracy strongly depends on the chosen method. Such questions do not arise with the Balanced POD since the only parameter that determines the dimension of the reduced system is the number of singular values retained. Second, the formulation of the superelement preserves confidentiality as the frequency-dependence of the elastic modulus is now hidden in the matrix Ê2 in equation (IV.17). In fact, this is not the case with modalbased methods mentioned above where the reduced system still contains the explicit dependence, see eq. ( 16) in [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF].

IV. 6 Numerical applications

This section deals with two applications. First, an academic case is studied in order to test the convergence of the method. A more realistic configuration is then considered.

IV. 6. 1 Cantilever sandwich beam

We consider the same structure as in [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF] which consists of a cantilever sandwich beam depicted in Figure IV. [START_REF] Bagley | Fractional Calculus -A Different Approach to the Analysis of Viscoelastically Damped Structures[END_REF]. It consists of two layers of steel separated by a damping layer made of viscoelastic material (Deltane 350). The corresponding material properties are given in Table IV. 1. In [57], the author compares the frequency response function (FRF) computed with several modal-based reduction methods. The FRF corresponds to the normal displacement at one point of the edge due to a normal force applied along the top edge (see arrows in Figure IV.1). Here, we propose to split the beam into two parts separated by the interface (red line): the super-element which is the biggest part of the structure and the host. Quadratic hexahedra finite elements are used which leads to 1218 and 11868 variables for the host and the substructure respectively. The elastic modulus of the viscoelastic material is approximated using the GHM mini-oscillator model by considering different order of approximation k = 1, 3, 5, 7 as shown in Figure IV. 2. Using 50 integration points in the BPOD method, 4 associated Reduced Order Models (ROMs) are generated. The SVD of the correlation matrix is performed using the algorithm svdsketch 2 . The number of singular values in Figure IV.3 gives the dimensions of the ROMs that are respectively 97, 156, 161 and 157. Figure IV.3 is instructive as it shows the convergence of the distribution of singular values with respect to order k. Clearly, convergence is reached as soon as k = 5 and one can anticipate a similar behaviour when computing the dynamic response of the super-element. One can also observe substantial data reduction since the BPOD reduces large matrices (IV.8) (reaching 100932 variables for k = 7) into a matrix system of very small dimension r = 157. The accuracy of the ROMs is evidenced in Figure IV.4 showing that the precision is strongly linked to the GHM model order k. The differences between the ROMs and the original system is measured using the same error criteria as employed by Rouleau et al. [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF]:

ϵ(ω) = ||q(ω) -q r (ω)|| 2 ||q(ω)|| 2 , (IV.21)
where q and q r are the displacement computed from the direct and reduced systems respectively, and || • || 2 is the ℓ 2 -norm. Plotting this criteria for the 4 ROMs returns error levels similar to the modal-based reduction methods investigated in [START_REF] Rouleau | A comparison of model reduction techniques based on modal projection for structures with frequency-dependent damping[END_REF]. Here again, convergence is reached as soon as k = 5 which was already reflected in Figure IV. 3. This suggests that the Singular Value Decomposition of the correlation matrix can provide a reasonable a priori indicator of the number of terms needed in the GHM model. For the sake of illustration, the different numerical costs are summarized in Table IV. 2. In this example, the construction of the reduced order model requires solving the system (IV.5) for 50 frequencies. The computational times related to these steps are about 400 s (note this corresponds to the execution of a Matlab code on workstation). More interestingly, the ratio in terms of computational time for the computation of the FRFs is only about 5 (3075 s with the original matrix and 560 s with the reduced model). These observations seem counter-intuitive if we have in mind that the overall number of variables has been divided by an order of magnitude. This could be explained by the special algebraic structure of the global matrix (some block matrices are full and complexvalued) whereas the full model preserves the sparsity associated with the finite element discretisation.

Models

Number of variables

Size of the global matrix system (IV. IV. 6. 2 Real-life application: a flat laminated plate of rectangular shape

A more realistic configuration, shown in Figure IV.5, is considered: it consists of a flat laminated plate of rectangular shape (1200mm × 1000m) with thickness 1.6mm × 0.8mm × 1.6mm, glued to a clamped elastic structure (2000mm × 1000mm × 20mm) through a layer of polyurethane (thickness 5mm) applied around the plate. The viscoelastic core is made of Deltane 350 and the constraint layers and the rest of the host are made of steel. The material properties are available in Table IV. [START_REF] Bagley | Fractional Calculus -A Different Approach to the Analysis of Viscoelastically Damped Structures[END_REF]. Quadratic hexahedral finite elements are used and this generates 38280 and 86019 variables for the host and the plate respectively. The interface corresponds here to the surface shared between the plate and the polyurethane, which totalizes 2520 variables. Reduced models of the laminated plate are built using 12 integration points (frequencies) in the range of [START_REF] Berkooz | The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows[END_REF]800 Hz].

In this scenario, the times for the computation of the frequency responses are shown in Table IV. 3. We may note that the reported values correspond to a reduced model with dimension r = 3010 which is automatically delivered by the SVD algorithm. In order to illustrate the effect of the truncation, two other reduced models are considered with

Models

Number of variables

Size of the global matrix system (IV. The two operational deformations before and after applying the reduction procedure (by taking r = 500) are represented in Figure IV.8. Note that the wavenumber of the flexural waves is higher for the laminated plate and this is due to the difference of thickness with the host structure. Clearly, the vibrational shapes are nearly identical, which is not the case where the number of singular values retained is too small (r = 150). Section B -Construction of the correlation matrix by z i = γ i Fq (see eq. (IV.3)), and the position vector R p (ω j ) can be written as

R p (ω j ) =        W R (ω j ) γ 1 (ω j )FW R (ω j ) . . . γ k (ω j )FW R (ω j )        . (IV.23) Here, W R (ω j ) = (-ω 2 j M + K(ω j )) -1 B (IV.24)
is the solution of the direct frequency-dependent system (IV.5) whose dimensions can be substantially smaller than the linearized system (IV.7). This signifies that the snapshot matrix R can be computed directly using any commercial software, for instance. Finally, by using the fact that matrices E and A are real-valued and symmetric, it is clear that

W S (ω j ) = (-ω 2 j M + K(ω j )) -1 C T . (IV.25)
Thus, by construction of the super-element, input and output maps are identical thus W S = W R and only the direct problem needs computing.

IV. B Construction of the correlation matrix

The particular structure of the snapshot matrices can also be exploited to compute the correlation matrix Z = S H R. Here again, memory can be saved since only the data sets W R = [ . . . W R (ω j ) . . . ] and W S = [ . . . W S (ω j ) . . . ] needs to be stored.

Considering the structure of R and S given in equation (IV.12), the correlation matrix Z has also a block structure as where the blocks Z ij can be decomposed using the position-velocity structure and

Z ij = S H (ω i )R(ω j ) = S H p (ω i ) S H v (ω i ) R p (ω j ) R v (ω j ) . (IV.27)
Since the derivative of a complex amounts to multiply it by i ω, it gives .28) This expression can also be simplified using equation (IV.3). In fact, after having decomposed S p (ω i ) and R p (ω j ) in the form of equation (IV.22), the matrix blocks Z ij can be written as

Z ij = S H p (ω i ) -i ω i S H p (ω i ) R p (ω j ) i ω j R p (ω j ) = 1 + ω i ω j S H p (ω i )R p (ω j ). ( IV
Z ij = 1 + ω i ω j W H S (ω i ) I + FF T k α=1 γ α (ω i )γ α (ω j ) W R (ω j ) (IV.29)
where it is reminded that only solutions of the original system (IV.5) needs computing and stored.

End of the paper.

IV. 8 Conclusion

A partial conclusion on this work has already been given in Section 7. Nonetheless, one may recall the main features.

In this chapter, the GHM/POD model reduction method has been employed to build superelements, which extends it to the dynamic substructuring field. Contrary to most modal based techniques, the form of the reduced system of equations does not require to know the value of the elastic modulus of the viscoelastic material at each frequency step. This is one of the main advantages of the method. Even though its construction can be computationally heavy, the superelement remains very accurate in providing very low level of errors. As a further work, it could be interesting to use the SO-BPOD instead of the classical BPOD. It is shown at the end of Chapter III that it could leads to more accurate models for the same number of singular values retained. Moreover, regaining the second-order structure could be more an advantage during the assembly process described in section 5. In fact, the first-order reduced system delivered by the BPOD must be connected to a second-order one. Intuitively, connecting systems with the same structure would lead to more accurate results.

V. 1 Introduction

Chapters III and IV provide the necessary ingredients in order to build a superelement of a mechanical structure made with viscoelastic materials. Once the Reduced Order Model (ROM) built, it is connected to another numerical model in order to perform global vibration studies. We focus here on vibro-acoustic simulations on vehicles where the ROM of a windscreen is connected to the car body of an automotive vehicle. In order to compute the acoustic pressure radiated inside the vehicle, the user must determine the vibrations of the surface in contact with the air cavity. This is a typical situation encountered in vehicle design. For instance, noises can emerge due to the road/tires contact creating structural vibrations that are propagated to the whole car body. When the windscreen is subject to structural excitation, it radiates noises in the passenger compartment. This phenomenon is clearly visible in the low frequency range of 0, 200Hz where acoustic modes of the cavity are coupled with structural modes of the windscreen. 1. In such a situation, several forces can be applied on the the car body (blue part) making each components of the vehicle vibrating. Among these components, windscreens are known to easily radiate sounds. These sounds can be numerically computed knowing the normal displacements over the surface in contact with the fluid. These data are easily obtained using classical FE models where the variables are the displacements at nodes. However, we consider here that the windscreen is reduced with the BPOD procedure. Then, the variables are no longer the displacements q but the generalized coordinates x, which requires to build a transformation q = Φx in order to recover the displacements of interest. Finally, the acoustic pressure resulting from the deformation of the interior (see Figure V.1a) face of the windscreen can be determined.

The quality of the acoustic field strongly depends on the quality of the displacements field over the windscreen. In fact, the normal displacements of the interior and exterior faces come from a reduced model. When a non adapted model reduction method is employed, the ROM can poorly approximate the original dynamic behavior. For BPOD, such a loss of quality arises when a too small number of POD modes are retained (see Chapter III). In this section, we propose to monitor the accuracy of the computed acoustic field for different "qualities" of ROM, i.e. for different numbers of singular values retained.

V. 2 Description of the models V. 2. 1 Problem and scheme

In this section, a description of the different vibroacoustic studies is given. The aim is to create an environment similar to the real-life situation of a windscreen glued to a car and radiating in the passenger compartment. To do so, let us consider the structure used in the second example of Chapter IV. As shown in Figure V.2, the superelement of the sandwich plate is connected to a 3-dimensional hard-walled acoustic cavity. Knowing the normal displacements of the wet surface Γ computed in vacuo, i.e. without considering the effect of the fluid on the structure, one can determine the acoustic field into the cavity Ω using an analytic procedure [START_REF] Guy | The transmission of sound through a cavitybacked finite plate[END_REF] presented in the following section. where p = p(x, y, z) is the value of the acoustic pressure at a coordinate (x, y, z), and k = ω/c 0 denotes the wave number with ω and c 0 the angular frequency and the speed of sound respectively. In addition, several boundary conditions must be respected. The first ones result from the presence of rigid walls:

∂p ∂n x=0 = ∂p ∂n x=Lx = ∂p ∂n y=0 = ∂p ∂n y=Ly = ∂p ∂n z=Lz = 0. (V.2)
Then, the speed continuity over the surface Γ must ensure

∂p ∂n Γ = -ρ 0 ω 2 q n , (V.3)
where ρ 0 is the air density and q n the normal displacement on the vibrating surface Γ.

Note that here the normal vector n denotes the z direction. The pressure p can be decomposed as a sum of weighted functions such that

p = l,m Φ lm ξ lm , (V.4)
where ξ lm is the coefficient of the function for the pair of parameter (l, m). In order to respect the boundary condition (V.3), each function is defined by

Φ lm (x, y, z) = cos(k l x) cos(k m y) cosh µ lm (L z -z) µ lm sinh µ lm L z , (V .5)
where

µ lm = k 2 l + k 2 m -k 2 , with k l = lπ L x and k m = mπ L y . (V.6)
Using this decomposition, one can rewrite equation (V.3) by multiplying it by Φ l ′ m ′ evaluated at z = 0 and integrated it over the wet surface Γ, which yields

Γ ∂ ∂z l,m ξ lm Φ lm z=0 Φ l ′ m ′ | z=0 dΓ = Γ -ρ 0 ω 2 q n Φ l ′ m ′ | z=0 dΓ. (V.7)
The linearity property of the integral and sum operators leads to

l,m ξ lm Γ ∂Φ lm ∂z z=0 Φ l ′ m ′ | z=0 dΓ = -ρ 0 ω 2 Γ q n Φ l ′ m ′ | z=0 dΓ, (V .8)
where

Φ l ′ m ′ | z=0 = cos(k l ′ x) cos(k m ′ y) cosh µ l ′ m ′ L z µ l ′ m ′ sinh µ l ′ m ′ L z , (V.9) and ∂Φ lm ∂z z=0 = cos(k l x) cos(k m y). (V.10)
Due to the orthogonality property of the functions Φ lm , the relation becomes [START_REF] Berthet | The balanced proper orthogonal decomposition applied to a class of frequency-dependent damped structures[END_REF] where I a (l, m) and I s (l, m, ω) are computed analytically and numerically respectively.

ξ lm Γ cos 2 (k l x) cos 2 (k m y)dΓ I a (l, m) = Γ -ρ 0 ω 2 q n cos(k l x) cos(k m y)dΓ I s (l, m, ω) , (V .
Analytic integral. Since the variables x and y are independent, the integral with square cosinus can be solved analytically such that

I a (l, m) = Γ cos 2 (k l x) cos 2 (k m y)dΓ = Lx 0 cos 2 (k l x)dx Ly 0 cos 2 (k m y)dy = α lm L x L y , (V .12)
where In doing so, one has

α lm =            1, if l = 0, m = 0, 1 2 , if l = 0, m > 0 or l > 0, m = 0, 1 4 , otherwise.
I s (l, m, ω) = -ρ 0 ω 2 Γ q n cos(k l x) cos(k m y)dΓ, (V.14)
and must solve the integral over the finite elements surfaces. Since the vibrating surface Γ coincides with the (⃗ x, ⃗ y) plane, the variable z is omitted. Moreover, the finite elements are defined in a reference coordinate system (u, v) as depicted in Figure V. 4. Then, the integral is evaluated in the reference coordinate system using a Gauss-Legendre quadrature rule, which allows to state

Γ e q e n (x, y) cos(k l x) cos(k m y)dΓ e ≈ g q e n (u g , v g ) cos(k l u g ) cos(k m v g ) det(J e g )δ g , (V.15)
where (u g , v g ) are the coordinates of the g th integration point expressed in the local system, det(J e g ) the Jacobian matrix and δ g the weight. The displacements q n (u g , v g ) of each element is obtained knowing two dataset: the displacements at each node of the element and the shape functions evaluated at the integration point. Final matrix system. The equation (V.11) must be solved for the whole pairs (l, m). The final matrix system for a single angular frequency ω can be written as

    . . . I a (l, m) . . .     I a        . . . ξ lm (ω) . . .        Ξ(ω) =     . . . I s (l, m, ω) . . .     I s (ω) . (V.16)
Then, the coefficients are computed from the simple diagonal matrix product

Ξ(ω) = I -1 a I s (ω), (V .17)
where I -1 a = diag(..., (α lm L x L y ) -1 , ...).

V. 3 Post-processing: vibroacoustic data

Before computing any acoustic data, one needs to determine the normal displacements over the interior or exterior vibrating surfaces. We consider here the model depicted in Figure V.2 and V.3 where the structure used in Chapter IV is connected to an acoustic cavity of 3.6 m 3 . The host is excited with a harmonic load while the sandwich plate is reduced through GHM/BPOD. This section starts in showing how to recover the displacements from the solution of the global system.

Instead of studying the whole vibroacoustic system or computing local data, we introduce two scalars that give global information on the elastic and acoustic domains. The first one is the mean square normal velocity v 2 n over the interior vibrating surface Γ that reports on the amount of deformation leading to sound pressure level inside the cavity. The second indicator is the average quadratic acoustic pressure p 2 over the acoustic domain Ω, which gives a report on the amount of sound in the cavity. Both indicators are computed along the frequency axis, and for different ROMs.

V. 3. 1 Recovering the displacements field over Γ

The solution of the equations of motion for the elastic domain leads to a vector of the form

   q 1 x2 u    , (V .18)
where q 1 are the displacements of the host, x2 are the generalized coordinates of the ROM of the sandwich plate and u the forces for each variable on the interface. In order to compute the acoustic field, one needs to evaluate the normal displacements over the vibrating surface in contact with the fluid. To do so, we employed the transformation used in the BPOD:

x 2 = Φx 2 , where x 2 = v 2 i ωv 2 with v 2 =          q 2 z 1 . . . z k          . (V.19)
Then, taking the appropriate lines of Φ in the matrix Φ allows to recover the required dataset directly from the generalized coordinates x2 , such as

q n = Φx 2 , (V.20)
where q n is a subset of q 2 containing the z-displacements for each node over Γ.

V. 3. 2 Mean square normal velocity over Γ

The normal velocity v n on the vibrating surface Γ is the derivative of the normal displacements q n . Then, the mean square velocity is defined by the integral V. 3. 3 Mean square acoustic pressure over Ω

v 2 n = 1 Γ Γ |i ωq n (x, y)| 2 dΓ. (V.
We introduce the mean square acoustic pressure

p 2 = 1 Ω Ω pp * dΩ, (V .24)
where p = p(x, y, z) is the acoustic pressure evaluated at a point in the cavity, the notation (•) * denotes the complex conjugate and L x , L y , L z are the x, y and z dimensions of the rectangular box. Using the eigendecomposition of p, the equation can be rewritten as

p 2 = 1 L x L y L z l,m Φ lm ξ lm l,m Φ lm ξ lm * dΩ. (V.25)
Due to the orthogonality of the functions Φ lm , the expression can be simplified as

p 2 = 1 L x L y L z l,m ξ lm ξ * lm Ω Φ lm Φ * lm dΩ. (V.26)
Based on the fact that the variables x, y and z are independent, the integrals over Ω can be analytically computed yielding

Ω Φ lm Φ * lm dΩ = Lx 0 cos 2 (k x x)dx Ly 0 cos 2 (k y y)dy Lz 0 ch lm (z)ch * lm (z)dz, (V.27) with the condensed form ch lm (z) = cosh µ lm (L z -z) µ lm sinh µ lm L z . (V.28)
Since the coefficient µ lm is either purely real or purely imaginary according the value of ω, these two cases must be distinguished. After calculations, it leads to

p 2 = 1 L z l,m |ξ lm | 2 α lm I z (µ lm )dz, (V .29)
where

I z (µ lm ) =            1 |µ lm | 2 sinh 2 (|µ lm |L z ) L z 2 + sinh(2|µ lm |L z ) 4|µ lm | , if µ 2 lm > 0, 1 |µ lm | 2 sin 2 (|µ lm |L z ) L z 2 + sin(2|µ lm |L z ) 4|µ lm | , if µ 2 lm < 0, (V.30)
where α lm is the same coefficient as defined in paragraph Analytic integral.

V. 4 Vibroacoustic study with different ROMs

This section gives numerical results such as the scalar indicators presented in section 3. Both mean square velocity and pressure are monitored along the frequency axis using different ROMs build from the same dataset of snapshots, but with a different number of singular values. As in any model reduction methods, projecting onto a small subspace leads to fast computations, but with a risk of losing information. Then, let us study the effect of the truncation of singular values on the vibroacoustic data in comparing the Full Order Model (FOM) with Reduced Order Models (ROMs) of dimensions 3010, 500 and 150 respectively.

V. 4. 1 Recovering the displacements

A first difference over the normal displacements obtained from the FOM and the ROMs respectively can be observed. In order to measure the quality of the recovering process, we introduce the following criterion

ϵ n = ||q FOM n -q ROM n || 2 ||q FOM n || 2 , (V.31)
which quantifies the error induced from the ROMs along the frequency axis. The Figure V.5 gathers this error criterion computed for the three ROMs. These results are not about the quality of the ROMs themselves, this is about their ability to recover the original variables, namely the physical displacements at each finite elements nodes. It can be seen that taking the 3010 singular values given by the algorithm svdsketch gives a level of error that is globally 1% or below. Considering 3010 singular values gives an error about 1% while the last ROM is not accurate enough to correctly recover the displacement field (about 100% of error). where p 0 = 2•10 -5 Pa is the reference acoustic pressure. It gives information about the global sound in the cavity along the frequency axis. The same remarks on the accuracy of the ROMs made for the mean square velocity can be done here. Moreover, more peaks are visible for this data. This is explained since this is not only due to the elastic resonances that are mainly responsible of the response in the low-frequency range, but also to the acoustic modes that provide numerous resonances for the mid-and high-frequencies.

V. 4. 4 Entire acoustic fields

It is also possible to compute the value of the acoustic pressure everywhere in the cavity. The field in the cavity is easily determined knowing the value of the coefficients ξ lm and the functions Φ lm (x, y, z) stored in vectors and matrices. Then, a simple matrixvector product gives the value of the pressure at any point in the cavity for a given frequency ω. The Figures V.8, V.9, V.10 and V.11 shows the normalized1 deformation of the whole structure and the corresponding acoustic field for a harmonic load at f = 400Hz. The scales are the same for the four figures. This allows to monitor the quality of the acoustic field according the displacement field of the vibrating surface in contact with the fluid. As expected, when the normal displacements of the surface Γ are sufficiently close to them obtained with the FOM (the represented in Figure V.8), the acoustic field is the same. This is confirmed with the third ROM that does not return the same acoustic data since the dynamic of the plate is poorly represented from a numerical point of view.

V. 5 Conclusion

This chapter has shown how to perform vibroacoustic simulations from a superelement created with GHM/BPOD. The model is made of two different domains: elastic (structure) and acoustic (cavity). It is representative of the real-life situation considering a windscreen glued to a car body and radiating sounds in the passenger compartment due to external vibrations.

First, the deformation of the structure is computed. At this stage, the dynamics of the superelement is described through a set of generalized coordinates that no longer have a physical meaning. The original variables (displacements) are recovered using the transformation matrix employed during the reduction procedure. An error criterion has been introduced in order to quantify the difference between the displacements obtained with the FOM and the ROMs. Then, the normal displacements over the wet surface are considered to compute the acoustic pressure in the cavity. Here, an analytic method is used to get the acoustic data, but numerical tools would be employed instead to model a real passenger compartment for instance.

We introduced two scalar indicators reporting global information over the elastic and acoustic domains respectively. These are the mean square normal velocity over the wet surface and the mean square acoustic pressure over the cavity. It allows to study the systems in a global way instead of monitoring data everywhere over the domains. Nonetheless, observing the pressure over the cavity is still possible and allows to visualize the difference between the FOM (model without reduction) and an inaccurate ROM (too few singular values).

This chapter has shown how to compute acoustic data from a vibrating structure reduced through the GHM/BPOD procedure. The results prove that the quality of the acoustic data are strongly linked to the quality of the reduced models. Then, the accuracy of any vibration or acoustic simulations rely on a single parameter, the appropriate choice of the number of singular values that determine both dimensions and precision of the ROM. the GHM/POD procedure can be used to build superelements of a structure made with viscoelastic materials. Here again, great results are obtained, even if an extra difficulty emerges. Examples in Chapter III shows that POD or BPOD can be easily and rapidly performed when only a few inputs or outputs are considered. However, when it comes to superelements with an interface made of a large number of variables, many precautions must be taken:

1. The number of snapshots must be optimized, because it directly determines the size of the matrix on which a Singular Value Decomposition (SVD) is performed. As the SVD is a very costly algorithm, taking this number as small as possible allows to lighten the numerical effort. Moreover, the time related to the computation of a single snapshots for large systems can be substantial, which is a second reason to minimize their number.

2. When working in the frequency domain, a rigorous application of the Fourier transform involves computing snapshots for positive and negative frequency-points. This doubles the size on the matrix on which the SVD is performed. Nevertheless, it is possible to consider only the positive values. This way, the SVD is performed on a complex but twice smaller matrix. The second example in Chapter IV shows that a sandwich structure with about 100 000 variables can be turned into a superelement of dimensions 500 only, while the interface counts 2520 variables. We must have in mind that this was made possible by considering and applying the 3 points cited above.

Finally, Chapter V proposes a contextualization where the second example of Chapter IV (sandwich plate glued to a steel structure) is extended. The aim is to recover the deformation over the plate and compute the acoustic pressure in a cavity. Once the equations of the structure-domain solved, the dynamics of the sandwich plate is described through generalized coordinates, so the first step is to recover the displacements (original coordinates) using the POD transformation matrix. Since then, the "quality" of recovering can be monitored thanks to an error criterion. Then, two indicators have been introduced in order to describe the structure-and acoustic-domains in a global way. Both are strongly linked to the quality of the superelement since the whole data comes out of the displacements field over the wet surface.

Future research perspectives

Further works associated to this study could take many directions. [START_REF] Bagley | Fractional Calculus -A Different Approach to the Analysis of Viscoelastically Damped Structures[END_REF]. The first suggestion would be to increase the scope of the GHM/POD procedure to parametric systems. This is motivated by the need to take into account more parameters related to the viscoelastic materials such as the temperature for instance, which also modifies the value of the stiffness.

2. After, one may work on how to build reduced models valid over a wide frequency range. The POD becomes costly when the range is wide because the number of snapshots could be high. An idea would be to split the range into several intervals, and to build a projection matrices for each of them. Once this is done, the collection of matrices could be concatenated and form the final projection basis valid over the whole range. In doing so, the number of snapshots is globally the same. However, one does not have to perform a single SVD on a large matrix, but several SVD on smaller matrices since every interval considers less snapshots.

3. It could be interesting to define a prior error criterion to the GHM/POD procedure.

Even if the svdsketch Matlab function directly sets the dimensions of the reduced model, it could be useful to know a priori the error between the original and reduced model without computing a posteriori the error on the transfer functions. 4. Finally, a work could be conducted on the numerical methods used to solve the global system connecting the superelement to an other numerical model. Even though the dimensions of the superelement matrices are relatively small, they are dense and complex, which could impact the performances of the solvers. Nevertheless, the global matrices always have the same block structure that could be suited for special solvers. Moreover, using the SO-BPOD instead of the BPOD would lead to reduced equations with the same structure of the host. This way, one can hope for better numerical stability and then better accuracy. The last lever would be to apply model reduction on the reduced system. The GHM/POD procedure has the advantage to build reduced systems with constant matrices, even if it represents the dynamics of a system with frequency-dependent properties. From this point, one may try to diagonalize the final system using modal-based algorithms for instance.

Appendices

Let U = u 1 , ..., u r and Λ = diag(λ 1 , ..., λ r ) be the first r eigenvectors and eigenvalues, the matrix of the POD modes is written as

Φ = RUΛ -1 2 . (A.8)
Since R T R is real symmetric and positive-definite matrix, its eigendecomposition coincides with its singular decomposition so the Singular Value Decomposition (SVD) procedure can be used instead of the eigendecomposition, which gives orthomnormal vectors.

The POD modes provides a basis that minimizes the error criteria (A.1). They are optimal in an energy sense since the first modes describes the dynamic concerning the more amount of energy that is directly related to the eigenvalues (or singular values) λ j such that (A.9)

Then, the larger λ j is, the more dynamically important the mode is. That is why only the first r modes are considered. This method is called the method of snapshots because it is based on the fact that data are given in discrete times where q(t j ) represents the state of the system at one of these times, which can be viewed as a "snapshots" of the system at t = t j . When the number of snapshots is smaller that the dimensions of the system, i.e. when J ≪ N , it is easier to solve the eigenvalue problem on R T R ∈ R J×J than RR T ∈ R N ×N . so Σ = diag(..., σ j , ...) and S H Φ are the eigenvalues and eigenvectors of ZZ H . Then, combining Eq. (B.7) and (B.8), one can write Finally, the left transformation matrix is

S H Φ = ZV 1 Σ -1 1 , S H Φ = S H RV 1 Σ -1 1 , (B.
Ψ H = ΦQ = RV 1 Σ -1 1 H SS H , = Σ -1 1 V H 1 R H SS H , = Σ -1 1 V H 1 S H R H S H , = Σ -1 1 V H 1 Z H S H , = Σ -1 1 V H 1 U 1 Σ 1 V H 1 H S H , = Σ -1 1 V H 1 V 1 Σ 1 U H 1 S

Figure I. 1 :

 1 Figure I.1: Laminated glass.

Figure I. 2 :

 2 Figure I.2: FE model of the windscreen of the car "Renault Zoé".

  (a) Deformation of a laminated glass. (b) The zigzag function used in the ZPST element.

Figure I. 3 :

 3 Figure I.3: Deformation of a laminated glass and approximation to the displacement. The sheets of glass are essentially bending while the PVB interlayer is subject to shear deformation.

Figure II. 1 :

 1 Figure II.1: Impact of damping on a single variable mechanical system.

Figure II. 2 :

 2 Figure II.2: Strain response to a stress step.

Figure II. 3 :

 3 Figure II.3: Stress response to a strain step.

Figure II. 5 :

 5 Figure II.5: Typical evolution of the real part of the elastic modulus.

Figure II. 6 :

 6 Figure II.6: Basic elements constituting rheological models.

Figure II. 7 :

 7 Figure II.7: Simplest rheological models.

  (a) Global structure.

Figure II. 15 :

 15 Figure II.15: Decomposition of a global structure into several substructures.

Figure II. 16 :

 16 Figure II.16: Connection of a numerical model of a windscreen on a numerical model of a car.

  j=1,...,J , S = . . . S(t j ) δ j . . . j=1,...,J , (II.104) that define the approximated Gramians such that P ≈ RR T , and Q ≈ SS T . (II.105)

2 OutFigure III. 1 :

 21 Figure III.1: Cantilever sandwich beam, 64 × 4 × 3 elements, 18 inputs and outputs (six nodes highlighted in the figure, and 3 components by node).

Figure III. 2 :

 2 Figure III.2: Real part and loss factor of the GHM model of 3M ISD112 (see TableIII.2) for 10 Hz to 1000 Hz. Data taken from[67].

Figure III. 3 :

 3 Figure III.3: Distribution of Hankel Singular Values (HSV) and relative error ϵ(ω) according the number of integration points for a 13 785 dofs finite element model with 19 197 GHM coordinated : lin ( ), log ( ), pol ( ) and gau ( ).

  Relative error ϵ for the transfer function.

Figure III. 4 : 2 M

 42 Figure III.4: Impact of the number of HSV : 10 ( ), 25 ( ), 75 ( ) and 212 ().

  Relative errors f 1 and f 2 .

Figure III. 5 :

 5 Figure III.5: FRFs for different number of HSV : 10 ( ), 25 ( ), 75 ( ) and 212 () compared to the original system (III.6) ( ).

Figure III. 6 : 5 :C(s 2 M

 652 Figure III.6: Finite element model of the windshield. The red and blue points correspond to the input and output nodes respectively.

  Relative error ϵ for the transfer function.

Figure III. 7 :

 7 Figure III.7: Hankel singular values and relative error ϵ between the full and reduced systems depending on the number of HSV retained: 100 ( ), 50 ( ) and 25 ( ).

  Relative error ϵ.

Figure III. 8 :

 8 Figure III.8: Frequency response function of the first output with u = [1 0 0 1 1 0] T as input vector for the windshield study. Comparison between three numbers of HSV : 100 ( ), 50 ( ) and 25 ( ) compared to the original system (III.6) with dimensions 29 364 ( ).

  [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF] that, thanks to the transformation v = Φ αβ v, reduce the second-order state-space (III.23) as -ω 2 Mαβ + i ω Dαβ + Kαβ v = Bαβ u y = Ĉαβ v (III.29) with M = Ψ H αβ MΦ αβ , D = Ψ H αβ DΦ αβ , K = Ψ H αβ KΦ αβ , B = Ψ H αβ B and Ĉ = CΦ αβ .

  e. the four possible pairs of Gramians: pp, pv, vp and vv. The first step is to plot the singular values, which are visible in Figure III.9. Then, choosing several fixed number of HSV to keep (100, 50 and 25), let us plot in Figure III.10 the FRF and associated errors found from the different ROMS.

  Truncation at r 1 = 100.

  Truncation at r 3 = 25.

Figure III. 9 :

 9 Figure III.9: Hankel singular values and errors on transfer function for the windshield study (see Figure III.6). Comparison between the BPOD and the four SO-BPOD for 3 different truncation thresholds (100, 50 and 25).

  Rel. err. for r 3 = 25.

Figure III. 10 :

 10 Figure III.10: Frequency responses and related errors for the windshield study (see Figure III.6). Comparison between the BPOD and the four SO-BPOD for 3 different truncation thresholds (100, 50 and 25).

. . 60

 60 Beginning of the paper .. . . . . . . . . . . . . . . . . . . . . . . 61

. 2 65 4. 3 65 5 66 6 67 6. 1 67 6. 2

 26536566671672 Balanced POD : an approximation to balanced truncation . . . . Notes on computational aspects . . . . . . . . . . . . . . . . . . . The BPOD as a dynamic substructuring technique . . . . . . Numericalapplications . . . . . . . . . . . . . . . . . . . . . . . Cantilever sandwich beam . . . . . . . . . . . . . . . . . . . .. . Real-life application: a flat laminated plate of rectangular shape 707 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . 74 A Fast computation of the snapshot matrix . . . . . . . . . . . . 74 B Construction of the correlation matrix . . . . . . . . . . . . . . 75 End of the paper . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

14 )

 14 Using the transformation x = Φx, the reduced system becomesi ω Êx = Âx + Ĝu y = Lx , (IV.15) with the reduced matrices Ê = Ψ H E Φ ∈ C r×r , Â = Ψ H A Φ ∈ C r×r , Ĝ = Ψ H G ∈ C r×p and L = LΦ ∈ C m×r .Though this is not place for a rigorous analysis, it is anticipated that the approximation equation (IV.13) should lead to a truncation error for the response function which is proportional to the sum of the discarded singular values.

Figure IV. 1 :

 1 Figure IV.1: Academic structure : sandwich beam split in two parts.

Figure IV. 2 :

 2 Figure IV.2: GHM models for Deltane 350 in the frequency range of [10, 800 Hz].

Figure IV. 3 :Figure IV. 4 :

 34 Figure IV.3: Singular values for the 4 GHM models of the beam given by the svdsketch Matlab function.

Figure IV. 5 :

 5 Figure IV.5: Laminated plate of rectangular shape connected to an elastic structure.

r

  = 500 and r = 150. These choices are identified in Figure IV.7. The superelement, i.e. the laminated plate, is used to compute the frequency responses due to a normal force exerted upon the host structure whereas the displacement is observed on the other side as depicted in Figure IV.5. In Figure IV.6 are shown the FRFs together with the ℓ 2 -norm error as defined earlier. Clearly, r = 500 is sufficient to provide very accurate results and keeping more singular values only have a marginal effect on the quality of the FRF.

Figure IV. 6 :

 6 FigureIV.6: FRF and relative error between FOM (Deltane 350) and ROM (BPOD on GHM order 5) for the plate.

Figure IV. 7 :

 7 Figure IV.7: Singular values for the laminated plate.

  (a) Original FE model of the plate of dimensions 86019. (b) Plate reduced order model of dimension 3010.(c) Plate reduced order model of dimension 500. (d) Plate reduced order model of dimension 150.

Figure IV. 8 :

 8 Figure IV.8: Deformation of the models at f = 400 Hz. The host structure has the same dynamic behaviour with or without reduction applied to the sandwich plate numerical model where 3010 or 500 singular values are retained, while taking only 150 is not sufficient and leads to inaccurate results.

  (a) Windscreen glued to the car body.

Figure V. 1 :

 1 Figure V.1: Recovering a displacement field from a reduced model of a windscreen.

  Figure V.2: 2D-slice of the vibroacoustic environment. The structure is the one used in Chapter IV where the bottom of the sandwich plate radiates in a 3D acoustic cavity.

(V. 13 )

 13 Figure V.3: 3D view of the vibroacoustic simulation environment.

Figure V. 4 :

 4 Figure V.4: Finite element mapping, relation between the local (left) and global (right) coordinate systems.

Figure V. 5 :

 5 Figure V.5: Relative errors ϵ n between the normal displacements obtained from the full models (without model reduction) and the three reduced models.

Figure V. 6

 6 Figure V.6 presents the averaged value of the square normal velocity over the vibrating surface in contact with the air in the cavity. The peaks indicate the resonances of the

Figure V. 6 :V. 4 . 3

 643 Figure V.6: Mean square normal velocity over the vibrating plate Γ.

Figure V. 7 :

 7 Figure V.7: Mean square acoustic pressure over the cavity Ω.

Figure V. 8 :

 8 Figure V.8: Inside system at f = 400Hz, FOM: sandwich plate with dimensions 86019.

Figure V. 9 :

 9 Figure V.9: Inside system at f = 400Hz, ROM 1: sandwich plate with dimensions 3010.

Figure V. 10 :

 10 Figure V.10: Inside system at f = 400Hz, ROM 2: sandwich plate with dimensions 500.

Figure V. 11 :

 11 Figure V.11: Inside system at f = 400Hz, ROM 3: sandwich plate with dimensions 150.

3 .

 3 When the two previous tricks are not sufficient, the SVD step is still a difficult stage. Nevertheless, since only the first singular values and vectors are of interest, one may use special algorithms that only returns a truncated SVD keeping the highest singular values, as the svdsketch Matlab function used in Chapter IV.

9 ) 1 (B. 10 )

 9110 which finally leads toΦ = RV 1 Σ -1that are the POD modes considering the new inner product in Eq. (B.3). The basis Φ = ..., ϕ j , ... is orthogonal respecting Eq. (B.3) such as⟨ϕ i , ϕ j ⟩ = ϕ H i Qϕ j . (B.11)However, by seeing ψ H j = ϕ H i Q as a left mode, the BPOD (POD with inner product (B.3)) can be seen as a biorthogonal projection: a orthogonal projection with different left and right transformation matrices. This projection is bi-orthogonal within the standard inner product, such that Ψ H Φ = I = ⟨Ψ, Φ⟩. (B.12)

  recalls the theory of the main methods:

	MSE : Modal Strain Energy
	IMSE : Iterative MSE
	ICE : Iterative Complex Eigensolution method
	MMSE : Modified MSE
	MM : Multi-Model approach
	MSE+C : MSE with first order correction
	MSE+R : MSE with displacement residuals

Table III

 III

				s]
	Linear	(sE -A) -1 G	65 964	675
	Quadratic	(s 2 M + s D + K) -1 B 32 982	312
	Original	(s 2 M + K(s)) -1 B	13 785	10.4

.3

: CPU times for the computation of a single snapshot. Results clearly show the interest in working with the original FE matrix system.

Table IV .

 IV 2: Calculation time of the FRF for the cantilever sandwich beam. Comparison between the FOM and the ROM with a GHM model of order 5. The offline cost of the GHM/BPOD procedure is about 400 s.

				Calculation time
			20)	(10 000 frequencies)
	FOM	n = 11868	13239	3074 s
	ROM	r = 161	1576	560 s

Table IV .

 IV 3: Calculation times of the FRF for the flat laminated plate. Comparison between the FOM and the ROMs with a GHM model of order 5 and with a different number of singular values retained.

				Calculation time
			20)	(10 frequencies)
	FOM	n = 86019	92619	1785 s
	ROM	r = 3010	9610	749 s
	ROM	r = 500	7600	739 s
	ROM	r = 150	7100	619 s

  [START_REF] Davis | Methods of Numerical Integration[END_REF] Using the linearity property of the integral, one can decomposed the mean square velocity as a sum on each finite element, such that

	v 2 n = ω 2	e	1 Γ e	Γ e	|q e n (x, y)| 2 dΓ e ,	(V.22)
	which is evaluated in the reference coordinate system
	v 2 n = ω 2	e	1 Γ e		g	|q e n (u g , v g )| 2 det(J e g )δ g .	(V.23)

  (B.[START_REF] Bouayed | Finite element analysis of the dynamic behavior of a laminated windscreen with frequency dependent viscoelastic core[END_REF] and then, after simplifications, the left transformation matrix of POD with inner product(B.3)) is Ψ H = U 1 S H . (B.14)It can be verifies that Ψ H Φ = I. Moreover, Φ and Ψ H are rescaling of the transformation matrices T and T -1 obtained within Cholesky factors in BT in Eq. (II.99). That allows to claim that BPOD modes are rescaling and approximated modes of BT modes using the method of snapshots.

H ,

The SVD can also be performed on S * R, which gives a reduced descriptor matrix Ê = T -1 ET. Applying the SVD on S * ER allows obtaining Ê = I.

In practice, the number of snapshots must be relatively small and in all cases it is expected that the size of the correlation matrix mJ and pJ should be much smaller than the original matrix size2N .

A rigorous application of the Fourier transform requires evaluating the integrals also along the negative real axis [-ω b , -ω a ], which can be achieved by considering the real-valued matrix [real(R) imag(R)] instead of the original complex-valued matrix R, as shown in[START_REF] Montier | Balanced Proper Orthogonal Decomposition Applied to Magnetoquasi-Static Problems Through a Stabilization Methodology[END_REF]. In practice, it was observed that integrals can be evaluated along the positive axis only (results have been checked to be identical) and this allows to perform the SVD on a smaller matrix.

The function svdsketch is a Matlab function based on the work of Yu[START_REF] Yu | Efficient Randomized Algorithms for the Fixed-Precision Low-Rank Matrix Approximation[END_REF]

The displacements have been normalized to set the absolute maximum to1.

Remerciements

Acknowledgements This work was carried out within the Chair for "intelligent and transparent surfaces for the automobile of the future", supported by the UTC, the UTC Foundation for Innovation and SAINT-GOBAIN. The authors would like to thank the European Regional Development Fund (ERDF) 2014/2020 and the Chair "intelligent and transparent surfaces for the automobile of the future" for the funding of this work.

End of the paper.

IV. 7 Conclusion

In this paper, an original dynamic substructuring technique dealing with structures made of viscoelastic materials is presented. The method applies to the large class of materials for which the complex modulus master curves fit well with the Golla-Hughes-McTavish mini-oscillators model. This specific form allows to recast the original problem, via linearization, into a conventional state-space representation which is amenable to Balanced Truncation model reduction techniques. In this regard, it is shown that the Balanced POD proves to be an excellent candidate for the numerical treatment as it considerably lightens the computational cost for the construction of the reduced basis which serves to approximate the displacement field of the structure. After construction, the reduced order model has an interesting feature that the state-space representation is preserved and this is exploited for the development of the super-element.

Two numerical examples of increasing complexity are studied in the last section of the paper. The first configuration, which consists of a cantilever sandwich beam made of two layers of steel separated by a viscoelastic material layer, serves to illustrate the effect of the number of terms retained in the GHM model on the accuracy of the computed solution (here the Frequency Response Function). In particular it is shown that the distribution of singular values of the correlation matrix can provide a reasonable indicator of convergence.

The second example is a more realistic configuration consisting of a flat laminated plate of rectangular shape connected to an elastic structure. The effect of the size of the reduced order model which is given by the number of singular values retained in the analysis is illustrated.

In both examples, it is shown that the super-element can be described with a number of variables which represents a very small fraction of the total number of degrees of freedom of the original system. The specific features of the method, as opposed to usual modalbased reduction methods are (i) the construction of the reduced order model is relatively straightforward and can deliver very accurate results irrespective of the frequency behavior of the elastic modulus of the viscoelastic layer, (ii) it preserves some confidentiality thus allowing companies to work on secure projects.

Appendix

IV. A Fast computation of the snapshot matrix

The particular structure of the state vector can be exploited in order to lighten the numerical effort required for the snapshots computation. First, it can observed that, by construction, vector R(ω j ) has the position (p) -velocity (v) structure

Second, if we let F the matrix defined so that Fq corresponds to the subset of q exclusively related to the viscoelastic material. Then, the vectors of dissipation coordinates are given

Chapter VI

Conclusion and outlooks

Summary and conclusions .

Summary and conclusions

The aim of this thesis is to propose a new model reduction method applied to numerical models of mechanical structures made with viscoelastic materials while respecting the two main criteria: (i) correctly taking into account the frequency-dependent damping properties and (ii) sharing reduced models without revealing the confidential data related to the viscoelastic material. This kind of situations rises in automotive industry where the glass manufacturer cannot directly share numerical models of windscreens without the risk of loosing confidential data related to the PVB.

The procedure presented in this manuscript can be decomposed in two major steps. First, the Golla-Hughes-McTavish (GHM) method is employed to properly consider the stiffness of the viscoelastic material. Then, a reduction procedure based on the Proper Orthogonal Decomposition (POD) efficiently reduces the dimensions of the system.

Chapter II has given a theoretical background on viscoelasticity, model reduction and dynamic substructuring techniques. Among this information, it emerges that reducing a numerical model with frequency-dependent stiffness -which is the case with viscoelastic materials -cannot be achieved using classical tools and requires particular strategies. The main disadvantage of these strategies lies in the form of the reduced equations where the stiffness of the viscoelastic material is clearly "visible". Even though it would still ensure point (i) (accuracy), it is not the case for point (ii) (confidentiality). In fact, sharing such reduced models involves giving sensitive data in the same time, which is precisely what manufacturers seek to avoid.

That is why the GHM/POD procedure detailed in Chapter III has been developed. The first step is to find the GHM model that properly fit the frequency-dependent experimental data. When it comes to particular materials with properties strongly linked to the frequency, a high-order GHM model is required, which leads to large augmented systems making modal-based model reduction methods intractable. Among the different reduction algorithms, the POD via the method of snapshots turns out to be suited to the GHM method. In fact, the projection matrix of the POD is found from a collection of snapshots, each one corresponding to an impulse response of the system. In the frequency domain, a snapshot is the solution of the system for a given frequency. Then, obtaining a set of snapshots requires solving a very large system many times, and storing the results. What we highlight in this work is that the solution of the very large system to reduce (GHM linearized) can be simply deduced from the solution of the initial system (FE system with frequency-dependent stiffness) which has not been enlarged by the GHM method. This way, the POD on a GHM system is feasible again even for a high-order GHM model.

Further works have been conducted about the reduction of the GHM system in its second-order form directly. In fact, the POD-based reduction procedure can only be applied on first-order systems, which involves linearizing the GHM system first before reducing it. The end of Chapter III has introduced the Second-Order Balanced Proper Orthogonal Decomposition (SO-BPOD) that applies POD directly on the second-order system, and returns a reduced order model also in a second-order form. This way, the classical form of mechanical systems of equations can be kept where the mass, damping and stiffness terms are clearly visible.

Chapter IV has been focused on the notion of superelement. It is explained how

Appendix A

The Proper Orthogonal Decomposition

As Berkooz written in [START_REF] Berkooz | The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows[END_REF], "the proper orthogonal decomposition is a procedure for extracting a basis for a modal decomposition from an ensemble of signals. Its power lies in the mathematical properties that suggest that it is the preferred basis to use in many circumstances". The origin of the Proper Orthogonal Decomposition (POD) is attributed to Lumley [START_REF] Lumley | The structure of inhomogeneous turbulence[END_REF] who has introduced it in the context of fluid dynamics in order to build reduced order models of fluids [START_REF] Lumley | Stochastic tools in turbulence[END_REF][START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF][START_REF] Holmes | Turbulence, coherent structures, dynamical systems and symmetry[END_REF]. The POD procedure has also been discovered in other disciplines under different names such as Karhunen-Loève method or Principal Component Analysis.

Giving the state vector q(t) ∈ R N with N the dimension of a generally large scale system, the idea of the POD procedure is to find a projection P of rank r ≪ N that minimizes the error

The projection matrix P comes from the first r eigenvectors of

where Zϕ j = ϕ j µ j is the N ×N eigenvalue problem to solve with µ j and ϕ j the eigenvalues and eigenvectors. Let Φ = ϕ 1 , ..., ϕ r be the first r eigenvectors of Z. Then, the basis P is given by

with ϕ j the so called POD modes. In practice, the vector q(t) is evaluated at a set of J discrete times, denoted by t j ∈ 0, T . In doing so, the integral can be approximated by the finite sum

or, in the matrix form, Z ≈ RR T , (A.5) with R = ..., q(t j ) δ j , ... ∈ R N ×J the matrix of snapshots and δ j the quadrature coefficients. In 1987, Sirovich [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF] founds from Eq. (A.5) that the N × N eigenvalue problem could be turned into a J × J eigenvalue problem, which lighten the numerical efforts when it comes to large systems. To do so, one might first solve

where R T R is now a J × J matrix, and λ j = µ j the same eigenvalues as with the previous problem. Once the vectors u j orthonormalized, the POD modes ϕ j are given by where Σ = diag(..., σ j , ...) contains the eigenvalues on its diagonal. The matrix Φ can be expressed from the spectral decomposition R H RV = VΣ R that yields

At this stage, two remarks can be observed : POD does not take into account the output equation of the state-space formulation, and the inner product used is the standard inner product as ⟨a, b⟩ = a H b. Nevertheless, the quantity x(ω) H Qx(ω) reflects the dynamical importance of the state x(ω), so computing the modes of ⟨R, R⟩ Q = R H QR rather than R H R would better reflect the dynamics of the system. We show in the following that making POD with ⟨a, b⟩ Q as inner product results in BT, and the BT modes (transformation matrix T defined in the previous subsections) are the BPOD modes normalized differently.

However, the singular modes are orthonormal such that

Then, it leads to the eigenvalue problem

List of Figures

List of Tables

List of Tables

Bibliography