Mathilde Chevreuil

Anthony Nouy

Model Order Reduction based on Proper Generalized Decomposition for the Propagation of Uncertainties in Structural Dynamics

Keywords: Uncertainty propagation, Spectral stochastic methods, Structural dynamics, Model reduction, Proper Generalized Decomposition, Tensor product approximation, Separated representations

A priori model reduction methods based on separated representations are introduced for the prediction of the low frequency response of uncertain structures within a parametric stochastic framework. The Proper Generalized Decomposition method is used to construct a quasi-optimal separated representation of the random solution at some frequency samples. At each frequency, an accurate representation of the solution is obtained on reduced bases of spatial functions and stochastic functions. An extraction of the deterministic bases allows for the generation of a global reduced basis yielding a reduced order model of the uncertain structure which appears to be accurate on the whole frequency band under study and for all values of input random parameters. This strategy can be seen as an alternative to traditional constructions of reduced order models in structural dynamics in the presence of parametric uncertainties. This reduced order model can then be used for further analyses such as the computation of the response at unresolved frequencies or the computation of more accurate stochastic approximations at some frequencies of interest. The dynamic response being highly nonlinear with respect to the input random parameters, a second level of separation of variables is introduced for the representation of functions of multiple random parameters, thus allowing the introduction of very fine approximations in each parametric dimension even when dealing with high parametric dimension.

INTRODUCTION

The quantification of uncertainties and the prediction of their impact appear to be essential for the robust design of mechanical systems. Uncertainty in structural dynamics have predominantly been tackled in fields of application such as the optimisation of structures as regards their natural modes or the identification of models from modal type or more generally frequency type experimental data. In this context of structural dynamics, the random response can be severely disturbed by low variabilities on the model. The impact of random variabilities are traditionally estimated using sampling techniques such as Monte Carlo simulations. Nevertheless, sampling techniques turn out to be computationally expensive since they require numerous solutions of deterministic problems in order to estimate accurately statistical quantities of interest. For the characterization of the natural frequencies and mode shapes of uncertain dynamical systems, perturbation-based methods have extensively been used but give good results only for small variation ranges of the input parameters [START_REF] Kleiber | The Stochastic Finite Element Method. Basic Perturbation Technique and Computer Implementation[END_REF]. Some works have contributed to improve these methods (see e.g. [START_REF] Nair | An approximate solution scheme for the algebraic random eigenvalue problem[END_REF][START_REF] Adhikari | Random matrix eigenvalue problems in structural dynamics[END_REF]).

Alternative approaches have emerged for the modeling and propagation of uncertainties which adopt a functional point of view in stochastic analyses, where random quantities are systematically seen as functions of random parameters representing basic uncertainties. This parametric vision, combined with approximation theory and numerical analysis, brought to the so-called spectral stochastic methods [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF][START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF][START_REF] Matthies | Stochastic finite elements: Computational approaches to stochastic partial differential equations[END_REF][START_REF] Nouy | Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations[END_REF][START_REF] Xiu | Fast numerical methods for stochastic computations: a review[END_REF][START_REF] Le Maitre | Spectral methods for uncertainty quantification with applications to computational fluid dynamics[END_REF]. These methods provide an explicit representation of the random response in terms of basic random parameters, and as a by-product, they can be favorably used for parametric analyses such as optimisation, identification or sensitivity estimation. They also allow for fast post-processing to retrieve many quantities of interest (such as the envelope of the frequency response function, sensitivity with respect to input parameters...), sometimes unprecise or even unreachable with classical sampling techniques. Non-parametric methods shall also be mentioned among the strategies of uncertainty quantification in dynamic analysis. These latter methods enable to take the model error into account and have recently been associated to parametric modeling of the uncertainties [START_REF] Soize | A nonparametric model of random uncertainties on reduced matrix model in structural dynamics[END_REF][START_REF] Soize | Generalized probabilistic approach of uncertainties in computational dynamics using random matrices and polynomial chaos decompositions[END_REF].

This paper is devoted to the computation of the random frequency response of dynamical models in the low frequency band and in the presence of stochastic parametric uncertainties. The dynamic response of uncertain structures is usually a highly non regular functional of input random parameters. Therefore, when adopting a functional point of view, it proves necessary to use very fine approximation spaces for the accurate representation of this response (or eventually to use enriched polynomial approximation spaces [START_REF] Ghanem | Efficient characterization of the random eigenvalue problem in a polynomial chaos decomposition[END_REF][START_REF] Ghosh | Stochastic convergence acceleration through basis enrichment of polynomial chaos expansions[END_REF][START_REF] Nouy | extended stochastic finite element method for the numerical simulation of heterogenous materials with random material interfaces[END_REF]), which make classical spectral stochastic approaches inefficient for the computation of the frequency response over a whole frequency band, especially when dealing with complex mechanical models.

Model reduction techniques are classically introduced to make random dynamic analyses affordable. They can be substructuring or reduced basis type strategies. For the midfrequencies, adapted frequency-domain model reduction techniques have been developed (see e.g. [START_REF] Sarkar | Mid-frequency structural dynamics with parameter uncertainty[END_REF]). In the low frequency range, the usual modal superposition is classically used. Stochastic Reduced Basis Methods [START_REF] Nair | Stochastic reduced basis methods[END_REF][START_REF] Dohnal | Joint uncertainty propagation in linear structural dynamics using stochastic reduced basis methods[END_REF] have also provided a mean of searching the response on a reduced basis of random vectors associated with a stochastic Krylov subspace that need to be computed first and independently for each frequency of interest. Let us note that the resulting reduced order model is not conventional since it is a projection on a stochastic reduced basis.

In this article, we present a model order reduction strategy based on Proper Generalized Decomposition (PGD) for the prediction of the low frequency response of uncertain structures within a parametric stochastic framework. PGD methods have been recently introduced for the a priori construction of separated representations of the solution of variational problems defined in tensor product spaces [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids[END_REF][START_REF] Chinesta | Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models[END_REF][START_REF] Ladevèze | The latin multiscalecomputational method and the proper generalized decomposition[END_REF][START_REF] Falco | A Proper Generalized Decomposition for the solution of elliptic problems in abstract form by using a functional Eckart-Young approach[END_REF][START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]. They have been introduced in [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF] in the context of uncertainty propagation under the name Generalized Spectral Decomposition (GSD). These methods can be seen as generalizations of spectral decompositions (Karhunen Loève decompositions or Singular Value Decompositions) for the a priori construction of separated representations. In the present context, PGD methods are extended to the complex framework and are first used for the construction of a low rank separated representation of the random solution at some frequency samples:

u(x, ξ) ≈ m i=1 w i (x)λ i (ξ) (1)
with the w i (x) being complex-valued deterministic functions and the λ i (called stochastic functions) being scalar valued functions of random parameters ξ modeling the uncertainties . At each computed frequency, an accurate representation of the solution is obtained on a reduced basis of spatial functions and stochastic functions. An extraction of the deterministic bases then allows for the generation of a global reduced basis yielding a reduced order model of the uncertain structure which appears to be accurate on the whole frequency band under study and for all values of input random parameters. This strategy can be seen as an alternative to traditional constructions of reduced order models in structural dynamics in the presence of parametric uncertainties. The obtained reduced order model can then be used for further analyses such as the computation of the response at unresolved frequencies or the computation of more accurate stochastic approximations at some frequencies of interest.

As mentioned before, the dynamic response being a highly non regular function of input random parameters, a fine stochastic approximation is usually necessary to well represent the random response. Two strategies are finally presented to control the accuracy of the representation of the random response. The first one is an adaptive stochastic approximation strategy. The second one is a more tractable approach to handle multidimensional problems involving a large number of parameters ξ = (ξ 1 , • • • , ξ r ) and fine stochastic approximations in potentially all dimensions. It is based on the introduction of a second level of separation of variables for the representation of multiparametric functions:

λ(ξ) ≈ Z k=1 φ 0 k φ 1 k (ξ 1 ) • • • φ r k (ξ r ). (2)
In this article, we propose a construction of these representations with a multidimensional PGD method, introduced in [START_REF] Nouy | Proper Generalized Decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF] in the context of uncertainty quantification. We show the feasibility of this parametric splitting in the context of structural dynamic analyses and we illustrate the behavior of different definitions of PGD. Note that other definitions and algorithms have been recently proposed for the construction of such separated representations of multiparametric stochastic functions [START_REF] Doostan | A least-squares approximation of partial differential equations with highdimensional random inputs[END_REF][START_REF] Khoromskij | Tensor-structured galerkin approximation of parametric and stochastic elliptic pdes[END_REF].

The paper is organized as follows. Section 2 presents the dynamic model problem and the associated partial differential equations defined on a tensor product space (product of functions of the space variables and functions of the random parameters). Section 3 introduces the PGD method for the construction of a low-rank separated representation [START_REF] Adhikari | Random matrix eigenvalue problems in structural dynamics[END_REF] of the solution at a fixed frequency. Several definitions have been proposed for PGD [START_REF] Nouy | An extended stochastic finite element method for solving stochastic partial differential equations on random domains[END_REF][START_REF] Nouy | Generalized spectral decomposition method for stochastic non linear problems[END_REF][START_REF] Nouy | Proper Generalized Decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF]. Different variants of PGD are here introduced and the so called Galerkin PGD is retained and illustrated on two numerical examples. We then take an interest in the computation of the frequency response over a low frequency band in Section 4, where we propose a strategy for the construction of a global reduced basis of deterministic modes valid for the entire frequency band under study. This reduced basis appears to be relatively insensitive to the approximation at the parametric level. We then show that it can be obtained using a coarse stochastic approximation, and consequently with a reduced computational effort, and that it yields to the construction of an accurate reduced order model which can be reused for more accurate computations. On that occasion a third numerical example is inserted. Finally, Section 5 introduces the two possible strategies for the obtention of very accurate stochastic approximations: an adaptive stochastic approximation technique and a method for the construction of separated representations (2) of stochastic parametric functions. For the construction of such representations, we introduce and compare different multidimensional versions of PGD, namely Galerkin PGD or Minimal Residual PGD.

STOCHASTIC STRUCTURAL DYNAMICS PROBLEM

A model problem in structural vibration analysis

We consider a forced vibration problem of a damped elastic structure defined on a domain Ω. Γ D and Γ N denote the Dirichlet and Neumann boundaries, with Γ D ∩ Γ N = ∅ and Γ D ∪ Γ N = ∂Ω. In the first part of the paper, we consider the frequency response for a fixed frequency ω ∈ I, with I = (ω min , ω max ) being the low frequency band of interest. Equations of the forced vibration problem at frequency read

div(σ(u)) + f = -ω 2 ρu on Ω u = 0 on Γ D σ(u) • n = g on Γ N (3)
where f and g are volume and surface load densities, and where σ(u) is the stress field related to the displacement field u by the following constitutive relation:

σ(u) = (1 + iωη)K : ε(u)
where K is the Hooke's elasticity tensor and where η is the damping parameter, eventually frequency dependent. In the following, we omit the dependence on ω of functions of the frequency for the sake of clarity.

Uncertainty modeling

Uncertainties are supposed to be represented with a finite number of random variables ξ = (ξ 1 , . . . , ξ r ) which are random parameters of the continuous mechanical model. We let (Ξ, B, P ξ ) the associated probability space, where Ξ ⊂ R r is the range of ξ, B is the associated σ-algebra and P ξ is the probability measure of ξ. Material uncertainties are introduced by considering random parameterized material parameters K(ξ), ρ(ξ), η(ξ). Uncertainties on the forcing terms are also introduced by considering random loads f (ξ) and g(ξ).

Finite element model

A spatial-weak formulation of the vibration problem (3) reads

u(ξ) ∈ V, Re (a(u(ξ), v; ξ)) = Re( b(v; ξ)) ∀v ∈ V (4)
where Re(•) denotes the real part and where a and b are respectively sesquilinear and antilinear forms defined by

a(u, v; ξ) = Ω ε(v) : (1 + iωη(ξ))K(ξ) : ε(u)dx - Ω ω 2 ρ(ξ)v • udx (5) b(v; ξ) = Ω v • g(ξ)dx + Γ2 v • f (ξ)ds (6)
where v denotes the complex conjugate of v. V is the following complex Hilbert space:

V = {v ∈ (H 1 (Ω)) d ; v = 0 on Γ 1 } (7)
We now introduce a finite element approximation space

V N = {v N = N i=1 ϕ i v i ; v i ∈ C} ⊂ V. The Galerkin approximation u N ∈ V N of u is then defined by u N (ξ) ∈ V N , Re a(u N (ξ), v; ξ)) = Re b(v; ξ)) ∀v ∈ V N (8)
resulting in the following system of equations for the components

u(ξ) = (u 1 (ξ) . . . u N (ξ)) ∈ C N of approximation u N : Re v H A(ξ)u(ξ) = Re(v H b(ξ)) ∀v ∈ C N (9)
where v H = v T denotes the conjugate transpose of v ∈ C N and where A(ξ) is a random matrix defined by

A(ξ) = -ω 2 M(ξ) + iωC(ξ) + K(ξ)
with M, K and C being the so-called mass, stiffness and damping matrices respectively. These matrices can be random. In the following, we only consider the finite element system of equations [START_REF] Ghanem | Efficient characterization of the random eigenvalue problem in a polynomial chaos decomposition[END_REF]. Let us note that these equations may be associated with other types of mechanical linear models (beam, plate, ...) such as those introduced in the numerical examples, or also to other types of physical models.

Spectral stochastic methods

Under suitable regularity assumptions on random parameters of the continuous model, a weak solution of ( 9) is searched in L 2 P ξ (Ξ; C N ), the space of real valued square integrable functions defined on the probability space (Ξ, B, P ξ ) which can be identified with the tensor product space C N ⊗ S, with S = L 2 P ξ (Ξ; R). We then introduce the following weak formulation of (9):

u ∈ C N ⊗ S, v, Au = v, b ∀v ∈ C N ⊗ S (10) with •, • defined for u, v ∈ C N ⊗ S by v, u = Ξ Re(v(y) H u(y))dP ξ (y) := E Re v(ξ) H u(ξ) (11)
where E(•) denotes the mathematical expectation. Spectral stochastic methods consist in searching an approximation of the solution under the form of a decomposition u(ξ) = P α=1 u α H α (ξ), where {H α (ξ)} P α=1 designates a basis of an approximation space S P ⊂ S. This approximation can be classically defined with a Galerkin projection:

u ∈ C N ⊗ S P , v, Au = v, b ∀v ∈ C N ⊗ S P (12)
For every frequency, searching the approximation therefore implies the solution of [START_REF] Khoromskij | Tensor-structured galerkin approximation of parametric and stochastic elliptic pdes[END_REF], which is a system of N × P equations in the complex field. For dynamic problems, the number of degrees of freedom N can be large if we want to capture rather small scales features of the solution. Furthermore, the solution can be a highly non-linear function of the random variables ξ and it may require the use of very fine approximation spaces S P (such as finite element [START_REF] Deb | Solution of stochastic partial differential equations using galerkin finite element techniques[END_REF] or wavelets [START_REF] Le Maître | Uncertainty propagation using Wiener-Haar expansions[END_REF]). In this context, classical Krylov type iterative algorithms for the solution of ( 12) can become prohibitive.

PROPER GENERALIZED DECOMPOSITION

In this section, we present a Proper Generalized Decomposition method (PGD) for the solution of problem ( 12) at a given frequency ω. The idea is to exploit the tensor structure of the solution u ∈ C N ⊗ S P and to look for an optimal separated representation of u under the form

u m (ξ) = m i=1 w i λ i (ξ) (13)
where the w i ∈ C N are deterministic vectors and the λ i ∈ S P are stochastic functions. Separated representation u m is constructed without knowing the solution u but only the equation [START_REF] Khoromskij | Tensor-structured galerkin approximation of parametric and stochastic elliptic pdes[END_REF] it is solution of. In the context of uncertainty propagation, this method is also known as Generalized Spectral Decomposition method (GSD). In [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF][START_REF] Nouy | Proper Generalized Decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF], different definitions and algorithms were proposed for generating the approximate representations (13) of the solution. These definitions can be seen as generalizations of Karhunen-Loève expansion where optimality is defined with non usual criteria. Some of these definitions are here extended to the complex framework in order to solve the stochastic problem at each frequency.

Progressive PGD

Here, we introduce a progressive definition of the PGD which consists in constructing the decomposition [START_REF] Kleiber | The Stochastic Finite Element Method. Basic Perturbation Technique and Computer Implementation[END_REF] progressively, where a rank-one element w m λ m is added to a previously computed decomposition u m-1 to obtain a rank-m approximation u m = u m-1 + w m λ m .

-A(u m-1 + w m λ m ) 2 = min w∈C N ,λ∈SP b -A(u m-1 + wλ) 2 ( 14
)
where • is the norm associated with the inner product •, • defined by [START_REF] Ghosh | Stochastic convergence acceleration through basis enrichment of polynomial chaos expansions[END_REF]. Necessary conditions of optimality of a couple (w m , λ m ) are

w m = F m (λ m ) = arg min w∈C N A(wλ m ) -b + A(u m-1 ) 2 (15)
λ m = f m (w m ) = arg min λ∈SP A(w m λ) -b + A(u m-1 ) 2 (16)
Denoting r m = b -Au m-1 , equations ( 15) and ( 16) can be rewritten

w m λ m , A H Aw m λ m = w m λ m , A H r m ∀ w m ∈ C N (17)
w m λ m , A H Aw m λ m = w m λ m , A H r m ∀ λ m ∈ S P (18)
In fact, equations ( 17) and ( 18) are necessary but not sufficient conditions for the optimality for w m λ m , which belongs to the set of rank-one elements {wλ; w ∈ C N , λ ∈ S P } which is not a linear subspace (nor a convex subset) of C N ⊗ S P . An optimal couple (w m , λ m ) must thus verify ( 15) and ( 16) simultaneously, i.e.

w m = F m (λ m ) and λ m = f m (w m ) (19)
Remark 1. In this finite dimensional framework, random matrices A and A H are interpreted as linear operators from

C N ⊗ S P to C N ⊗ S P .
The downside of this definition of the PGD, though robust for non symmetric problems, is that it may lead to a slow convergence of the decomposition u m . Furthermore, tensor based methods capitalize on the separated representation of the operator and right-hand side. In this minimal residual formulation, the operator and right-hand are A H A and A H b and present a dramatically higher rank than the initial operator A and right-hand side b.

Progressive

Galerkin PGD Knowing u m-1 , the progressive Galerkin PGD consists in defining a new rank-one element w m λ m that verifies two Galerkin orthogonality criteria

w m λ m , Aw m λ m = w m λ m , b -Au m-1 ∀ w m ∈ C N (20)
w m λ m , Aw m λ m = w m λ m , b -Au m-1 ∀ λ m ∈ S P (21)
Such as for Minimal Residual PGD, we introduce mappings F m : S P → C N and f m : C N → S P , respectively defined by

w m λ m , AF(λ m )λ m = w m λ m , r m ∀ w m ∈ C N (22)
w m λ m , Aw m f m (w m ) = w m λ m , r m ∀ λ m ∈ S P (23)
with

r m = b -Au m-1 . A new couple (w m , λ m
) is thus defined as the optimal couple which verifies simultaneously

w m = F m (λ m ) and λ m = f m (w m ) ( 24
)
This definition is not associated with an optimality criteria. However, in the case of a selfadjoint operator A (which is not the case here), this problem could be classically reformulated as a minimization problem of a quadratic convex functional. Equations ( 22) and ( 23) would then be necessary (but not sufficient) conditions for optimality of a couple (w m , λ m ) and optimality of such a couple would be naturally defined from the underlying minimization problem. Although not based on an optimization criterion for general operators A, this PGD definition based on Galerkin orthogonality proves efficient to capture good rank-one approximations in many cases. The following interpretation as a pseudo eigenproblem allows to give a meaning to optimality.

Interpretation as a generalized spectral decomposition

In [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF], it is shown that problem [START_REF] Nair | Stochastic reduced basis methods[END_REF], and problem [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF] for the case of a self-adjoint operator A, can be interpreted as a pseudo eigenproblem where w m (respectively λ m ) is the dominant pseudo eigenfunction of the mapping

T m = F m • f m (respectively T * m = f m • F m ). Mappings T m and T *
m can be respectively interpreted as right and left (pseudo) correlation operators of the solution u, this solution being interpreted as an operator from C N to S P . The couple (w m , λ m ) can be further interpreted as a couple of right and left (pseudo) singular vectors of u, thus leading to a generalized singular value decomposition u m = m i=1 w i λ i . We can prove that in the case of a deterministic matrix A, these definitions exactly coincide with a classical definition of a singular value decomposition of u (Karhunen Loève decomposition), up to a change of basis in C N (the decomposition being optimal with respect to the metric on C N induced by A T A for the minimal residual PGD, and A for the Galerkin PGD). In the general case of a non self-adjoint operator A, problem [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF] which defines potentially optimal elements can still be interpreted as a pseudo eigenproblem. A couple (w m , λ m ) satisfying ( 24) can be interpreted as eigenfunctions of mappings T m and T * m . Optimality, however not clearly defined, could still be associated with the upper spectrum of these mappings, which would justify the behavior of the algorithm proposed below. For a detailed discussion on these interpretations, see [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF].

Algorithm

An alternated direction algorithm for solving [START_REF] Nair | Stochastic reduced basis methods[END_REF] or ( 24) is used to find the approximation of an optimal couple (w m , λ m ). Starting from an initial function λ 0 ∈ S P (chosen randomly), iterates {w k , λ k } k≥1 are recursively defined by

w k+1 = F m (λ k ) and λ k+1 = f m (w k+1 ), (25)
where mappings F m and f m are defined by equations ( 20) and ( 21) for the Galerkin PGD, and by equations ( 15) and ( 16) for the Minimal Residual PGD. Based on the interpretation of section 3.2, the above algorithm corresponds to a power-type algorithm to solve the pseudo eigenproblem which defines an optimal couple (w m , λ m ). The resulting construction of the separated representation (13) can then be interpreted as a power-type method with deflation for the computation of the generalized singular value decomposition † . In practice, only few iterations are sufficient to obtain a good estimation of an optimal couple (w m , λ m ).

Below, we detail the computational aspects of the application of mappings F m and f m in the case of Galerkin PGD. Note that applications of these mappings for Minimal Residual PGD can be simply obtained by replacing A by A H A and b by A H b.

3.3.1. Application of mapping F m For a given λ, the calculation of w = F m (λ), with F m defined in [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF], appears as a simple deterministic problem of size n (2n in the real field)

Re w H A λ w = Re w H r m λ ∀ w ∈ C N (26)
where, in the case of a Galerkin PGD, A λ = E(Aλλ) and r m λ = E(r m λ). Considering the real and imaginary parts of the quantities, equation ( 26) can be written under the form

E(A R λ 2 ) -E(A I λ 2 ) E(A I λ 2 ) E(A R λ 2 ) w R w I = E(r m R λ) E(r m I λ) (27)
where (•) R and (•) I denote the real and imaginary parts of a quantity (•). Concerning operator A, we have A R = Kω 2 M and A I = ωC.

Application of mapping f m

For a given w, the calculation of λ = f m (w), with f m defined in ( 23), appears as a simple stochastic algebraic equation approximated in S P (problem of size P )

E λA w λ = E λr m w ∀ λ ∈ S P (28)
where, in the case of a Galerkin PGD,

A w (ξ) = Re(w H A(ξ)w) = w R w I T A R (ξ) -A I (ξ) A I (ξ) A R (ξ) w R w I (29) r m w (ξ) = Re(w H r m (ξ)) = w R w I T r m R (ξ) r m I (ξ) (30)

Updated progressive PGD

In order to improve the convergence of the progressive PGD, an update of the stochastic functions can be performed [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF]. Let us suppose that a rank-m decomposition u m = m i=1 w i λ i has been obtained and let us introduce the notation

u m (ξ) = WΛ(ξ) (31)
where

W = (w 1 • • • w m ) ∈ C N ×m and Λ = (λ 1 • • • λ m ) T ∈ (S P
) m denote the sets of deterministic vectors and stochastic functions respectively. Given the deterministic basis W of this rank-m decomposition u m , stochastic functions Λ can be updated by solving

W Λ, AWΛ = W Λ, b ∀ Λ ∈ (S P ) m (32)
In other terms, it can be seen as a stochastic problem projected on the reduced basis W, W H AW and W H b being respectively the reduced operator and right hand side. More precisely, the obtained approximation u m = WΛ is a Galerkin projection of the solution u on the subspace span{w i } n i=1 ⊗ S P of C N ⊗ S P . The non-dimensional analysis considers a unitary mass density. The Young modulus E is chosen as a uniform random variable on (0.9, 1.1) and is expressed as a function of a uniform random variable ξ 1 on Ξ 1 = (0, 1) (E = 0.9 + 0.2ξ 1 ). We take the damping η = 0.005 2 ξ2 , with ξ 2 a uniform random variable on Ξ 2 = (0, 1). It corresponds for η to a log-uniform distribution on (0.005, 0.01). Sources of uncertainty are thus represented with 2 independent uniform random variables ξ = (ξ 1 , ξ 2 ) with values in Ξ = Ξ 1 ×Ξ 2 = (0, 1) 2 . Figure 2 shows some realizations of the frequency response function of the out of plane displacement of the upper right node of the structure. In figure 3 are plotted some realizations of the solution u of (9) for a given frequency. These two figures show that different values of input random parameters yield completely different dynamic behaviors of the structure, turning a resonance frequency to an anti-resonance frequency for instance. In order to well describe these irregularities with respect to parameters, we introduce an approximation space S P = S 1 P1 ⊗ S 2 P2 where S l P l are finite element spaces with 32 elements (uniform partition of Ξ 1 ) and polynomial degree 4 for l = 1 and with 1 element and polynomial degree 2 for l = 2. The dimension of the stochastic approximation space is P = 480. Convergence of PGD. We study the convergence of the progressive Galerkin PGD with respect to the rank m of the decomposition for three frequencies ω = 0.41 rad.s -1 , 0.67 rad.s -1 and 0.84 rad.s -1 . We choose as the reference solution the semi-discretized solution

u ref ∈ C N ⊗ S.
The convergence of the rank-m PGD u m is estimated with the relative error

ǫ m = u ref -u m u ref (33)
where

v 2 = E(Re(v H v))
. This norm is estimated with Monte-Carlo simulations with 1000 samplings which correspond to a converged estimation of the error indicator ǫ m . Figure 4 shows the relative error ǫ m with m for three frequencies. Only the first part of the curve before the plateau illustrates the convergence of the PGD, the plateau being the residual error due to the stochastic approximation. A good convergence is observed for each frequency, even for ω = 0.67 rad.s -1 and ω = 0.84 rad.s -1 , which are frequencies in the resonance zones of the structure. For these two frequencies, a relative error of the order of magnitude of the error due to stochastic approximation is rapidly reached with m ≈ 4. We observe an even faster convergence for ω = 0.41 rad.s -1 , which is a frequency outside of the resonance zones. In fact, away from the resonances, the response is much smoother with respect to ξ than in the resonance zones. Figure 5 illustrates the convergence with m of u m for one realization of input We now consider the response surfaces of the out of plane displacements of the top right node of the structure. Figure 6 (resp. 7) shows the response surfaces for the reference solution and for the rank-10 PGD approximation, for the frequency ω = 0.67rad.s -1 (resp. ω = 0.84rad.s -1 ). These figures illustrate the strong irregularity of the response with respect to the random parameters (more specifically with respect to ξ 1 ). The PGD is able to capture and well describe these irregularities with low rank decompositions. .

Figure 8 compares the displacement of the top right node obtained with the reference solution and the rank-10 PGD approximation for the three studied frequencies and for a fixed value of η = 0.005. As mentioned earlier, we can observe that the solution is much smoother for ω = 0.41 rad.s -1 (outside of the resonance zones) than for ω = 0.67 rad.s -1 and ω = 0.84 rad.s -1 (inside the resonance zones). Again, we observe a good concordance between the reference solution and the rank-10 PGD approximation.

Example 2

We now study the frequency response of the 3D structure presented in figure 9 with homogeneous Dirichlet boundary conditions on Γ 1 and a harmonic load density on Γ 2 . At the spatial level, a tetrahedral mesh composed of 9363 elements is used and the model contains N = 6186 degrees of freedom. We consider a unit mass density and a damping η = 0.001. The Young modulus E is chosen as a uniform random variable on (0.8, 1.2) (that is to say a variability of ±20% around the mean value) and is expressed as a function of a uniform random variable ξ with values in Ξ = (0, 1) (E = 1 + 0.4(ξ -1 2 )). At the stochastic level, a finite element type approximation is used with a partition of Ξ into 16 elements with polynomial degree 4 (P = 80). Figure 10 illustrates some realizations of the frequency response of a particular output degree of freedom of the 3D structure.

For ω = 1.67rad.s -1 in a resonance zone of the structure, a PGD approximation is computed. We focus on the displacement in the direction of the excitation of the output node (located in figure 9(a)). The semi-discretized solution, denoted u ref is taken as the reference solution. It is compared in figure 11 to the rank-4 PGD approximation u 4 . Once more, a very low rank PGD approximation is able to capture and well describe the irregularity if the response.

MODEL REDUCTION STRATEGY FOR BAND ANALYSIS

In this section, we take an interest to the computation of the frequency response over a certain frequency band. We propose a model reduction strategy based on PGD for the efficient solution of problem ( 12) over the whole frequency band. The PGD method introduced in section 4 can be seen as a strategy for the a priori construction of a reduced basis of functions which is optimal for the representation of the solution at a given frequency. However, these reduced bases can then be exploited for the generation of a global reduced order model which appears to be valid on the whole frequency band of interest.

Construction of a global reduced basis

Let {w i } M i=1 denote the global reduced basis and let introduce the associated matrix

W = (w 1 , . . . , w M ) ∈ C N ×M
For a new frequency ω, we proceed as follows. We compute a first approximation

u M = M i=1 w i λ i (ω, ξ) := WΛ(ω, ξ)
which is defined as the Galerkin projection of the solution u(ω, ξ) on the reduced approximation space associated with W. That means that stochastic functions Λ = (λ 1 , . . . , λ M ) T ∈ (S P ) M are defined as the solution of the following problem:

E( Λ T A W Λ) = E( Λ T b W ) ∀ Λ ∈ (S P ) M (34)
with A W (ω, ξ) = Re(W H A(ω, ξ)W) and b W (ω, ξ) = Re(W H b(ω, ξ)). The precision of the solution u M (ω, ξ) = WΛ(ω, ξ) is then evaluated by the following error indicator

ε res (u M ) = b -Au M b (35)
If the approximation u M does not satisfy a desired precision, the PGD method is used in order to build a higher rank decomposition

u M+K = u M + M+K i=M+1 w i λ i
with K such that u M+K satisfies the error criterion. The global reduced basis W is then enriched with the new computed vectors {w i } M+K i=M+1 , thus leading to the augmented global reduced basis W = (w 1 , . . . , w M , w M+1 , . . . , w M+K )

We finally denote by m(ω) = M + K the rank of the obtained approximation. Let us note that for a given desired precision, the required rank m(ω) depends on the way the frequency band is scanned. In the following examples, we choose to scan the frequencies by increasing values.

Example 1

We consider the example 1 described in section 3.5.1 and we choose a frequency band which contains the 10 first eigenmodes of the mean structure (mean value of parameters). We compute the response for 50 frequency samples {ω β } 50 β=1 uniformly distributed over the frequency band. For each frequency, we solve the problem with a precision ε res ≤ 10 -3 .

The problem is first solved for each frequency independently. Figure 12(a) shows the rank m(ω) obtained at each frequency. When a resonance appears, the representation of the solution clearly requires a higher rank decomposition (i.e. a richer basis of deterministic modes). A total of 265 deterministic modes have been computed to represent the random dynamic response when considering the 50 frequencies separately ( 50β=1 m(ω β ) = 265). We now consider the strategy using the global reduced basis. For the same required precision 10 -3 , figure 12(b) shows the number of deterministic modes added to the global basis at each frequency and the evolving size M of the reduced basis W. We recall that we scan the frequency band from the left to the right, which explains the increase of M with ω. We notice that enrichment of the global basis coincides with the resonances of the frequency response, for which new features which are not present in the global basis are detected in the solution. A total of only 30 deterministic modes are here computed for the prediction of the random dynamic response over the whole studied frequency band. The reuse of the reduced basis W enables to considerably reduce the computation cost.

Example 2

We proceed likewise with the second example described in section 3.5.2. The considered frequency band contains the first 10 eigenmodes of the 3D structure. We compute the response for 50 frequency samples {ω β } 50 β=1 uniformly distributed over the frequency band. For each frequency, we solve the problem with a precision ε res ≤ 10 -3 .

The problem is first solved for each frequency independently. Figure 13(a) shows the rank m(ω) obtained at each frequency. A total of 232 deterministic vectors are generated for the analysis of the response over the frequency band ( 50 β=1 m(ω β ) = 232). We now consider the strategy using the global reduced basis. For the same required precision 10 -3 , figure 13 Figure 14 shows some deterministic modes of the reduced basis W. We can easily have access to the envelope of the frequency response over the frequency band with a simple postprocessing. Figure 15 illustrates the very good quality of the envelope (quantiles of orders 0.1% and 99.9%) obtained with the fine stochastic approximation (see section 3.5.2). We can see that the Monte-Carlo simulations plotted on the same figure are inside the envelope. We have finally obtained a global reduced basis W with small dimension 32 that defines a reduced order model of an uncertain structure which gives an accurate representation of the dynamic response on the whole frequency band. It is further illustrated in the following section.

Robust model order reduction

In the previous section, we have introduced a way of constructing a global reduced deterministic basis W which is valid for a frequency band. This basis allows for the construction of a reduced order model which can be reused for computing non computed frequencies of the frequency band of interest and also for computing more accurate stochastic approximations of the solution. Actually, the global basis can be constructed with a reduced computational cost using a coarse stochastic approximation, and reused for computing a fine stochastic approximation of the associated reduced order model. More precisely, a fine approximation u M (ω, ξ) ∈ C N ⊗ S P on the frequency band can be computed in two steps:

• First, a global basis W is computed using a coarse stochastic approximation (P = P coarse ) and q frequency samples {ω β } q β=1 uniformly distributed over the frequency band I.

• Then, for any other frequency of interest in I, u M is calculated by just updating the stochastic functions Λ by solving equation ( 34) with a fine stochastic approximation (P = P f ine ).

As illustrated below, it appears that the global reduced basis W which is obtained by using a coarse stochastic approximation is very close to the one obtained by using directly a fine stochastic approximation. This can be interpreted as a property of robustness with respect to stochastic approximation.

Illustration on example 1

Let us consider the two-plate structure of example 1. In the following, I is the frequency band containing the first 10 eigenmodes. q = 50 frequency samples over I are used for the construction of the global deterministic bases. Figure 16 shows two global deterministic bases W which are obtained with the same required precision ε res ≤ 10 -3 and with the respective dimensions of stochastic approximation spaces P = 4 and P = 240. The angle between the subspaces spanned by the two global bases is 0.39rad (very close subspaces). Furthermore, we can note that the angle between the subspaces spanned by the modal basis made of the 10 first modes of the mean structure on the one hand and by W on the other hand is 4.10 10 -10 rad (resp. 2.04 10 -10 rad) with W = W P =240 (resp. W = W P =4 ), that is to say the supspace spanned by the modal basis of the mean structure is a subspace of that spanned by W. PGD expansions using a coarse stochastic approximation are thus able to capture a good approximation of the global reduced basis.

To support this observation, a separated representation u 30 = WΛ ∈ C N ⊗ S P f ine is now computed into two steps: we first generate the global basis W using a coarse stochastic approximation P = P coarse = 4. Then, we update for each frequency the stochastic functions Λ using the fine approximation P = P f ine = 240. Figure 17 shows the response surface of the approximation thus obtained before and after updating the stochastic functions with the fine stochastic approximation for three frequencies which belong to the set of frequencies used for the construction of W. The approximation u 30 matches perfectly the reference solution after the stochastic functions are updated using the fine stochastic approximation.

The global reduced basis of dimension 30 has been computed using 50 frequency samples (and a coarse stochastic approximation). Using this reduced basis, the approximation u 30 = WΛ is now calculated by updating the stochastic functions Λ for the 50 frequency samples used for the construction of W on the one hand and for 200 other frequency samples on the other hand. The envelopes of the frequency response (quantiles of orders 0.1% and 99.9%) are then post-processed and plotted in figure 18. The envelope obtained in figure 18(b) is of very good quality even though the frequencies for which the approximation u 30 is computed are not those used in the construction of W. The envelope is notably well computed for the low frequencies ω < 0.1rad.s -1 . In conclusion, it is possible to compute the global reduced basis of deterministic modes W with a coarse stochastic approximation and a coarse frequency sampling, and then to update a global basis of stochastic functions using a suited refined stochastic approximation (the global basis W can be enriched if needed). This construction of the reduced basis enables considerable computational savings and has proved robust.

Illustration on example 2

Similar observations can be done on the 3D structure. We define I the frequency band containing the 10 first eigen modes of the mean structure. q = 50 frequency samples are used for the construction of W with a coarse stochastic approximation P = P coarse = 12. The stochastic functions are then updated with a fine stochastic approximation P = P f ine = 80. The parametric response surface thus obtained matches very well the reference solution once the stochastic functions are updated with the fine stochastic approximation, as illustrated in figure 19.

Illustration on example 3

Let us introduce a third example which consists of a 3D elastic structure presented in figure 20. Homogeneous Dirichlet conditions are imposed at the base of the arch and the structure is excited on its left vertical face with a normal harmonic load density. The structure has four cylindrical inclusions with a random Young modulus for each of them. At the spatial level, a tetrahedral mesh composed of 8159 elements is used and the model contains N = 5346 degrees of freedom. We choose the mass density ρ = 1 and the damping coefficient η = 0.01. We denote by E 1 to E 4 the Young moduli of the four inclusions distributed clockwise looking at figure 20(a), E 1 corresponding to the Young modulus of the bottom left inclusion. E j , with j = {1, 2, 3, 4}, are uniform random variables on [0.95, 1.05]. The sources of uncertainties are then represented with 4 independent uniform random variables ξ = {ξ 1 , ξ 2 , ξ 3 , ξ 4 } with values in Ξ = × 4 j=1 Ξ j = (0, 1) 4 and E j = 1 + 0.1(ξ j -1 2 ). Figure 21 shows samples of the frequency response function of a degree of freedom over the frequency band containing the first 9 eigen modes of the mean structure.

Here we take an interest in the narrow frequency band I = [0.25rad.s -1 , 0.32rad.s -1 ] that contains the eigen modes 6 to 9 of the mean structure and for which variabilities of the input parameters are of important incidence on the frequency response. We construct a global basis W over this frequency band with a precision ε res ≤ 10 -3 . It is generated using q = 10 increasing frequency samples uniformly spread over I and a coarse stochastic approximation consisting in a degree 1 polynomial approximation in each dimension (P j = 2), yielding a total dimension of P = 16. We finally obtain a reduced basis W = W P =16 of 75 deterministic modes. This reduced basis is compared with the modal bases V s made of the eigen modes 6 to 9 of the deterministic structures associated with different realizations of the input parameters ξ s by computing the angles between the subspaces spanned by W and V s . Using 100 samplings we find a maximum angle of 1.82 10 -6 rad. We can conclude that the subspace spanned by the global basis W contains almost surely all the subspaces spanned by the modal bases V s of the deterministic structures. This reduced basis W thus generated proves to be an accurate basis for the computation of the frequency response for a frequency ω ∈ I as will be seen in section 5.2.2.

STRATEGIES FOR FINE STOCHASTIC APPROXIMATIONS

Since the dynamic response may happen to be highly nonlinear with respect to the random variables, the dynamic analysis requires a fine stochastic approximation. This points out the need for strategies for fine stochastic approximation to improve the accuracy of the representation of the random response. We here propose and discuss two strategies for fine stochastic approximation.

Adaptive stochastic approximation

Various strategies of adaptive approximation can be considered. The strategy used here is a simple h-adaptive finite element approximation. Let {Ξ i } N i=1 be a partition of Ξ into N non-overlapping elements. The error estimation per element on the partition {Ξ i } N i=1 of Ξ is evaluated based on two solutions computed with different polynomial degrees. Let us denote by u p M and u p+q M the rank M approximations computed with a previously constructed deterministic basis W of M modes and an update of the stochastic functions with polynomial degree respectively p and p + q, with q > 0. Considering u p+q M as the reference solution, the global error estimation on Ξ is defined by

ǫ = u p+q M -u p M L 2 P ξ (Ξ;C n ) u p+q M L 2 P ξ (Ξ;C n ) (36)
and the error estimation per element is defined by

ǫ i = u p+q M -u p M L 2 P ξ (Ξ i ;C n ) √ P i u p+q M L 2 P ξ (Ξ;C n ) ( 37
)
where

v 2 L 2 P ξ (Ξ i ;C n ) = Ξ i
Re(v H v)dP ξ and P i is the probability measure of the ith element. It comes out that providing ǫ i < e ensures that ǫ < e. Thus an element is splitted if the error ǫ i on the element is higher than a refinement tolerance e. Conversely two adjacent elements i and i + 1 are merged if

P i ǫ 2 i + P i+1 ǫ 2 i+1 P i + P i+1 < e γ 2 with γ > 1 ( 38
)
where P i ǫ 2 i represents the error contribution of element i. The adaptive strategy is illustrated on example 1 for which the most sensitive stochastic dimension is the first dimension. For this dimension, we set p = 4 and q = 1. Let γ = 10. The resulting partition of Ξ 1 for a relative precision of e = 10 -3 and the response surface in terms of ξ 1 of the out of plane displacement are illustrated in figure 22 for frequencies ω = 0.41Hz, ω = 0.67Hz and ω = 0.84Hz.

Adaptive strategies provide appropriate stochastic approximation spaces as regards the response surfaces. Naturally, bigger elements are obtained when the response is smooth and smaller elements are needed where the response presents irregularities, which leads to a high dimension of the stochastic approximation space. The authors refer to [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic diffential equations[END_REF][START_REF] Wan | Error control in multi-element generalized polynomial chaos method for elliptic problems with random coefficients[END_REF][START_REF] Mathelin | Dual-based a posteriori error estimate for stochastic finite element methods[END_REF][START_REF] Almeida | Solution verification, goal-oriented adaptive methods for stochastic advection-diffusion problems[END_REF] for an efficient assessment of the error estimation and some improvements of the adaptive methods. However, these adaptive strategies become untractable when increasing the stochastic dimension even for a moderate stochastic dimension. It then proves necessary to propose alternative approaches.

Separated representation in the parametric space

In this section, we leave aside the concerns about the choice of the most adapted stochastic approximation. The motivation is to use a straightforward and sufficiently fine stochastic approximation. So far, the stochastic algebraic equations (32) (resp. ( 28)) have been solved using a classical spectral stochastic approximations which requires the solution of a system of respectively M × P (resp. P ) equations, with a possibly very high P . Here, we exploit the tensor product structure of the stochastic space S P = S 1 P1 ⊗ . . . ⊗ S r Pr and we introduce a method for the a priori construction of a separated representation of the stochastic functions Λ ∈ R M ⊗ S P under the form

Λ(ξ) ≈ Λ Z (ξ) = Z k=1 φ 0 k φ 1 k (ξ 1 ) . . . φ r k (ξ r ), (39)
with φ 0 k ∈ R M and φ j k ∈ S j Pj for j = 1 . . . r, and where Z denotes the rank of the decomposition. This type of representation allows for high dimensional stochastic problems to be dealt with (see [START_REF] Nouy | Proper Generalized Decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF]). Note that it can be interpreted as a generalized spectral expansion

Λ Z (ξ) = Z k=1 φ 0 k H k (ξ) (40)
with basis functions H k (ξ) = φ 1 k (ξ 1 ) . . . φ r k (ξ r ) which are not selected a priori but constructed in an optimal way.

PGD method

For the presentation of the methodology, we consider the solution of problem [START_REF] Soize | Generalized probabilistic approach of uncertainties in computational dynamics using random matrices and polynomial chaos decompositions[END_REF], equivalently rewritten as follows:

Λ ∈ R M ⊗ S P , Λ, BΛ = Λ, d ∀ Λ ∈ R M ⊗ S P (41)
where random matrix B(ξ) and random vector d(ξ) are defined by

B(ξ) = Re(W H A(ξ)W) ∈ R M×M , d(ξ) = Re(W H b) ∈ R M Remark 2.
Note that the solution of equation ( 28) can be also recasted in the same format by considering

B(ξ) = Re(w H A(ξ)w) ∈ R and d(ξ) = Re(w H r m ) ∈ R.
Minimal residual PGD. We first introduce a progressive definition of separated representation (39) based on a Minimal Residual PGD. We suppose that a rank-(Z -1) separated representation Λ Z-1 has been computed. A new rank-one element

φ 0 Z φ 1 Z (ξ 1 ) . . . φ r Z (ξ r ) is then defined by the optimization problem B(Λ Z-1 + φ 0 Z φ 1 Z . . . φ r Z ) -d 2 = min φ 0 ,φ 1 ,...,φ r B(Λ Z-1 + φ 0 φ 1 . . . φ r ) -d 2 (42)
The following r + 1 coupled equations are necessary conditions of optimality of the new rank one element φ 0 Z φ 1 Z . . . φ r Z :

< φ 0 φ 1 Z . . . φ r Z , B T Bφ 0 Z φ 1 Z . . . φ r Z >=< φ 0 φ 1 Z . . . φ r Z , B T d >, ∀ φ 0 ∈ R M < φ 0 Z φ 1 . . . φ r Z , B T Bφ 0 Z φ 1 Z . . . φ r Z >=< φ 0 Z φ 1 . . . φ r Z , B T d >, ∀ φ 1 ∈ S 1 P1 • • • < φ 0 Z φ 1 Z . . . φ r , B T Bφ 0 Z φ 1 Z . . . φ r Z >=< φ 0 Z φ 1 Z . . . φ r Z , B T d >, ∀ φ r ∈ S r

Pr

An alternated minimization algorithm is then used to construct a rank-one element solution of (42), consisting in minimizing successively on φ 0 , φ 1 , ... φ d .

Galerkin PGD. A definition based on Galerkin orthogonality criteria can also be introduced, by searching a new rank one element φ 0 Z φ 1 Z . . . φ r Z that satisfies simultaneously:

< φ 0 φ 1 Z . . . φ r Z , Bφ 0 Z φ 1 Z . . . φ r Z >=< φ 0 φ 1 Z . . . φ r Z , d >, ∀ φ 0 ∈ R M < φ 0 Z φ 1 . . . φ r Z , Bφ 0 Z φ 1 Z . . . φ r Z >=< φ 0 Z φ 1 . . . φ r Z , d >, ∀ φ 1 ∈ S 1 P1 • • • < φ 0 Z φ 1 Z . . . φ r , Bφ 0 Z φ 1 Z . . . φ r Z >=< φ 0 Z φ 1 Z . . . φ r Z , d >, ∀ φ r ∈ S r Pr Preprint
A new rank-one element is then computed by applying an iterative algorithm for solving the above system of equations (see [START_REF] Nouy | Proper Generalized Decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF] for details). Such as for the PGD in 2 dimensions presented in section 3, optimality of a rank-one element can be clearly defined only for symmetric problems associated with an optimization problem. This definition however proves efficient for many applications. Below, this definition is compared with the Minimal Residual definition on our problem of interest.

Remark 3. In [START_REF] Nouy | Proper Generalized Decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF], the author introduces additional update steps that can be carried out in some or all of the dimensions in order to improve the convergence of the separated representation. In the following example, we perform one updating step along each dimension (see [START_REF] Nouy | Proper Generalized Decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF] for details).

Illustration

Example 1. The methodology is applied to example 1 reusing the global basis W = W P =4 constructed in section 4.2.1. The set of stochastic functions Λ Z is constructed with the definition above where the approximation spaces are those defined in section 3.5.1. We then search Λ

Z (ξ 1 , ξ 2 ) ∈ R M ⊗ S 1 P1 ⊗ S 2
P2 with M = 30, P 1 = 160 and P 2 = 3. We study the convergence of Λ Z with respect to the rank Z. In the same way as in section 3.5.1, the convergence is illustrated with the relative error between u Z = WΛ Z and the semi-discretized solution

u ref ∈ C N ⊗ S ǫ Z = u ref -u Z u ref (43)
Figure 23 depicts the relative error ǫ Z with Z for the three frequencies, ω = 0.41 rad.s -1 , ω = 0.67 rad.s -1 and ω = 0.84 rad.s -1 , using a Galerkin PGD and a Minimal Residual PGD for the definition of the separated representation Λ Z . Only the Minimal Residual PGD proves robust for all three of the frequencies. The relative error reaches the plateau (observed before in figure 4) representing the error due to the stochastic approximation with Z = 22. The progressive Galerkin PGD happens to finally converge only for relatively smooth parametric responses but fails for the construction of a separated representation of a highly nonlinear parametric response. In figure 24, are plotted the response surfaces of the displacement u Z with Z = 5, 10 and 30 obtained with the Minimal Residual PGD. They have to be compared to those of figures 6 and 7.

Example 3. The strategy is now applied on example 3 where the global basis W = W P =16 constructed in section 4.2.3 is reused. The set of stochastic functions Λ Z is constructed with the Minimal Residual PGD. We introduce the approximation spaces S j Pj of finite element type with 40 elements (uniform partition of Ξ j ) and polynomial degree 4 for j = {1, 2, 3, 4} and

Λ Z is searched in R M ⊗ S 1 P1 ⊗ S 2 P2 ⊗ S 3 P3 ⊗ S 4
P4 with M = 75 and P 1 = ... = P 4 = 200. Λ Z is computed for the frequency which corresponds to the seventh eigen frequency of the mean structure (ω ≈ 0.27 rad.s -1 ) and which is not in the set of q = 10 frequencies used for the construction of W. we have ε res (Λ 50 ) = 1.7 10 -3 . These results support the observation made in section 4.2 of the ability to construct a good quality reduced basis W over a frequency band I that can be reused for the computation of frequency responses in I. and ω = 0.84 rad.s -1 (right figures), and for an increasing rank Z from top to bottom figures

CONCLUSION

In this paper, we have first shown that the PGD method enables to automatically generate a quasi-optimal separated representation of the frequency response of a structure on a reduced basis of deterministic modes and stochastic functions. A global reduced basis of deterministic modes valid for a low frequency range can be constructed progressively with a reuse and enrichment of bases extracted from PGD expansions at some frequency samples. This global basis appears to be robust with respect to stochastic approximation and it allows the construction of an accurate reduced order model that can be reused for further analyses.

In particular, it can be used for computing the dynamic response at unresolved frequency samples. It can also be used for computing more accurate stochastic approximations, by using adaptive stochastic approximation methods or a further separated representation method for the representation of multiparametric stochastic functions. This latter method enables to tackle multidimensional parametric problems with fine approximation in each parametric dimension. For the construction of these separated representations, a Minimal Residual PGD method has been introduced. It has proved robust but it yields to a rather slow convergence of separated representations. This points out the need for further investigations and the possible introduction of more efficient definitions of PGD for the present context of structural dynamics.

3. 1 . 1 .

 11 Progressive Minimal Residual PGD Problem (12) can be reformulated as an optimization problem by introducing a minimal residual formulation. Knowing u m-1 , a new couple w m λ m can be defined by b

  Description of the example, stochastic modeling and approximation. The method is applied on a slightly damped elastic structure made of two plates represented in figure1. It is subjected to homogeneous Dirichlet boundary conditions on the part Γ 1 of the boundary and to a harmonic loading on the part Γ 2 of the boundary. The complementary part of the boundary is a free boundary. At the space level, a finite element approximation is used, where the mesh is composed of 1778 DKT plate elements and the discrete deterministic model contains N = 5556 degrees of freedom.

Figure 1 .

 1 Figure 1. Description of example 1: elastic plate structure under harmonic bending load. Geometry and boundary conditions (a) and finite element mesh (b)

Figure 2 .Figure 3 .

 23 Figure 2. Samples of the frequency response function of the modulus of the out of plane displacement of the upper right node of the two-plate structure

Figure 4 .

 4 Figure 4. Convergence with m of um for three frequencies ω = 0.41 rad.s -1 , 0.67 rad.s -1 and 0.84 rad.s -1 . Error indicator ǫ m estimated with Monte-Carlo simulations parameters and for different frequencies.

Figure 5 .

 5 Figure 5. Convergence with m of PGD approximation um. From left to right, the solutions are plotted for three different frequencies and a given sample of random parameters. The top line represents the reference solution and the following lines represent rank-m PGD approximations for ranks m = 1, m = 2 and m = 10

u 10 Figure 6 .u 10 Figure 7 .

 106107 Figure 6. Response surface of the modulus of the out of plane displacement |u| of the upper right node of the two-plate structure for ω = 0.67 rad.s -1.

1 Figure 8 .

 18 Figure 8. Modulus of the out of plane displacement of the upper right node of the two-plate structure as a function of random variable E and for a fixed value of η = 0.005. Comparison of the rank-10 PGD u10 with the reference solution u ref for different frequencies ω.

Figure 9 .

 9 Figure 9. Description of example 2. Geometry and boundary conditions (a) and finite element mesh (b)

Figure 10 .Figure 11 .

 1011 Figure 10. Samples of the frequency response function of the output degree of freedom of the 3D structure

.

Figure 12 .

 12 Figure 12. Example 1. (a) rank m(ω) of the PGD approximation for each of the 50 independently computed frequencies and (b) progression of the rank M of the decomposition by reusing the deterministic basis over 50 computed frequencies, and number of added modes at each frequency (frequencies are computed from the left to the right)

Figure 13 .

 13 Figure 13. Example 2. (a) rank m(ω) of the PGD approximation for each of the 50 independently computed frequencies and (b) progression of the rank M of the decomposition by reusing the deterministic basis over 50 computed frequencies ω, and number of added modes at each frequency (frequencies are computed from the left to the right)

Figure 14 .

 14 Figure 14. Some deterministic modes of the global reduced basis W

Figure 15 .

 15 Figure 15. Envelope of the frequency response of the 3D structure: quantiles of orders 0.1% and 99.9%. Envelope of PGD approximation associated with a residual error of 10 -3 and samples of the reference solution.

  .

Figure 16 .

 16 Figure 16. Global deterministic bases W of the two plate structure obtained with two different stochastic approximations: (a) P = 4 and (b) P = 240

1 Figure 17 .

 117 Figure 17. Modulus of the out of plane displacement of the upper right node of the two-plate structure as a function of random variable E and for a fixed value of η = 0.005. Comparison between the reference solution u ref and the PGD approximation using the global basis W for frequencies (a)ω = 0.41 rad.s -1 , (b)ω = 0.67 rad.s -1 and (c)ω = 0.84 rad.s -1 . W is computed with P = 4 and with 50 frequency samples.

( a )

 a Update of the stochastic functions on 50 frequency samples (b) Update of the stochastic functions on 200 frequency samples

Figure 18 .

 18 Figure 18. Envelopes (quantiles of orders 0.1% and 99.9%) of the frequency response of the modulus of the output degree of freedom for the approximation u30 = WΛ. W is first computed with 50 frequency samples and using a coarse stochastic approximation, and Λ is updated at the same 50 frequency samples (a) or at 200 different frequency samples (b), using a fine stochastic approximation.Black curves represent samples of the reference response surface.

Figure 19 .

 19 Figure[START_REF] Nair | Stochastic reduced basis methods[END_REF]. Modulus of the output degree of freedom of the 3D structure as a function of random variable E(ξ). Comparison between the reference solution u ref and the approximation using the global basis W for the frequency ω = 1.67 rad.s -1 . W is first computed with P = Pcoarse = 12 and with 50 frequency samples.

Figure 20 .Figure 21 .

 2021 Figure 20. Description of example 3: mesh and boundary conditions of a 3D elastic arch structure.

  (a) ω = 0.41 rad.s -1 , P = 30 (b) ω = 0.67 rad.s -1 , P = 150 (c) ω = 0.84 rad.s -1 , P = 75

Figure 22 .

 22 Figure 22. Comparison between the h-adapted approximate solution and the reference solution u ref (top figures) and final partition of Ξ1 and error ǫi per element (bottom figures).

1 ,Figure 23 .

 123 Figure 23. Convergence with Z of uZ . Error indicator ǫ Z estimated with Monte-Carlo simulations

  ω = 0.67 rad.s -1 , Z = 5 ω = 0.84 rad.s -1 , Z = 5 ω = 0.67 rad.s -1 , ω = 0.84 rad.s -1 , ω = 0.67 rad.s -1 , ω = 0.84 rad.s -1 , Z = 30

Figure 24 .

 24 Figure 24. Response surface of the modulus of the out of plane displacement of the upper right node of the two-plate structure (for the separated approximation uZ) for ω = 0.67 rad.s -1 (left figures)and ω = 0.84 rad.s -1 (right figures), and for an increasing rank Z from top to bottom figures

† Note that for the case of a deterministic operator A, when the obtained separated representation coincides with a classical SVD of u, the proposed algorithm is exactly a power algorithm with deflation for the computation of this SVD.

Contract/grant sponsor: French National Research Agency (ANR); contract/grant number: TYCHE (ANR-2010-BLAN-0904)