Dimitri Goutaudier
email: dimitri.goutaudier@ensam.eu

Laurent Berthe

Francisco Chinesta

Proper Generalized Decomposition with time adaptive space separation for transient wave propagation problems in separable domains

Keywords: Proper Generalized Decomposition (PGD), Transient wave propagation, Time adaptive space separation, Separable domain, Scalar wave equation, Elastodynamics

Transient wave propagation problems may involve rich discretizations, both in space and in time, leading to computationally expensive simulations, even for simple spatial domains. The Proper Generalized Decomposition (PGD) is an attractive model order reduction technique to address this issue, especially when the spatial domain is separable. In this work, we propose a space separation with a time adaptive number of modes to efficiently c apture t ransient wave p ropagation i n s eparable domains. We combine standard time integration schemes with this original space separated representation for empowering standard procedures. The numerical behavior of the proposed method is explored through several 2D wave propagation problems involving radial waves, propagation on long time analyses, and wave conversions. We show that the PGD solution approximates its standard finite element solution counterpart with acceptable accuracy, while reducing the storage needs and the computation time (CPU time). Numerical results show that the CPU time per time step linearly increases when refining t he m esh, even with implicit time integration schemes, which is not the case with standard procedures.

Introduction

Transient wave propagation analyses are often necessary to fully understand the dynamic behavior of a mechanical system or the details of a mechanical process. This is notably the case for laser shock processes that consist in applying a sudden pressure, up to several GPa, locally on the surface of a material part [START_REF] Fabbro | Physics and applications of laser-shock processing[END_REF][START_REF] Bolis | Physical approach to adhesion testing using laser-driven shock waves[END_REF]. The outcomes of understanding these processes at the wave propagation level are important to optimize the laser parameters [START_REF] Ghrib | Generation of controlled delaminations in composites using symmetrical laser shock configuration[END_REF][START_REF] Ecault | A study of composite material damage induced by laser shock waves[END_REF] or to determine the mechanical properties of the medium [START_REF] Ehrhart | Development of a laser shock adhesion test for the assessment of weak adhesive bonded cfrp structures[END_REF][START_REF] Sagnard | Development of the symmetrical laser shock test for weak bond inspection[END_REF]. These laser shock processes are numerically explored through parametric studies on transient wave propagation problems in plate domains, by varying the loading parameters, the material properties or the mechanical behavior laws [START_REF] Ecault | Laser shock adhesion test numerical optimization for composite bonding assessment[END_REF][START_REF] Cuq-Lelandais | Spallation generated by femtosecond laser driven shocks in thin metallic targets[END_REF]. It would be beneficial to reduce the computation time and the storage requirements associated to these numerous simulations.

Many numerical tools have been developed over the past decades to solve transient wave propagation problems. The classical approach consists in using a time integration scheme with the Finite Element Method (FEM). It is however difficult t o o btain a ccurate s olutions w ith t his m ethod b ecause o f n umerical d ispersion a nd dissipation, resulting in non-physical wave velocities, period elongations and amplitude decays [START_REF] Deraemaeker | Dispersion and pollution of the fem solution for the Helmholtz equation in one, two and three dimensions[END_REF]. In addition, the error between the exact and the numerical solutions may severely increase as the waves travel in the medium, because of spurious numerical high frequency oscillations related to the Gibb's phenomenon. The numerical solution can straightforwardly be improved by refining t he m esh a nd b y c arefully s electing t he t ime s tep a ccording t o both the element size and the selected time integration scheme [START_REF] Bathe | Finite Element Procedures[END_REF]. Instead of refining t he m esh, w hich i ncreases the computation time, another approach is to select an appropriate time integration scheme with good dissipation and dispersion properties. For instance the Bathe method is an effective implicit time integration scheme for computing accurate transient wave propagation solutions [START_REF] Noh | Performance of an implicit time integration scheme in the analysis of wave propagations[END_REF][START_REF] Noh | The bathe time integration method with controllable spectral radius: The ρ∞-bathe method[END_REF]. The Spectral Element Method (SEM) has been developed to reduce the numerical dispersion with respect to the standard FEM [START_REF] Patera | A spectral element method for fluid dynamics: Laminar flow in a channel expansion[END_REF][START_REF] Seriani | Dispersion analysis of spectral element methods for elastic wave propagation[END_REF]. Instead of using standard finite elements, high-order Lagrangian-based finite e lements w ith s pecific in tegration po ints ar e us ed to bu ild a di agonal mass matrix. Consequently, the SEM can be effective for computing transient wave propagation solutions in explicit time integration. Similarly, Lagrangian-based finite elements enriched with harmonic functions have been proposed to significantly i mprove t he a ccuracy w ith r espect t o s tandard fi nite el ements [1 5]. Fi nally, ap proaches based on discontinuous finite e lement f ormulations h ave b een d eveloped t o s olve t ransient w ave p ropagation problems involving sharp gradients. For instance the Time Discontinuous Galerkin method (TDG) using space-time finite elements is effective for reducing high frequency spurious oscillations [START_REF] Hulbert | Space-time finite element methods for second-order hyperbolic equations[END_REF].

The computation time and the storage requirements of these methods significantly increase when rich discretizations are used to study the details of the wave propagation. Reduced order modeling is an active field o f research offering promising approaches to deal with this issue. In the present study we do not intend to propose a new model order reduction technique. We focus on the Proper Generalized Decomposition (PGD) which is an attractive technique when the geometry of the domain is separable, as for the laser shock simulations described above. The PGD has been developed to overcome the so-called curse of dimensionality when it comes to solve a multidimensional problem [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF]. This technique consists in computing on the fly a separated form representation of the solution, directly from the weak formulation of the problem. The main idea is to solve several problems of lower dimensions than the initial problem to alleviate the computational costs. When it comes to solve time-dependent partial differential equations, the PGD is usually implemented by considering the space-time separation presented in [START_REF] Ladevèze | The LATIN multiscale computational method and the proper generalized decomposition[END_REF]:

u(x, y, t) ≈ M ∑ i=1 S i (x, y)T i (t) (1)
where u represents a time-dependent scalar field defined over a 2D space domain. This space-time separation avoids using a time incremental solver, by analogy with the Proper Orthogonal Decomposition (POD). If the length of the sum M is small (between 10-100 for instance) then the time dependent solution can be stored with significantly less storage capacities. Unlike the POD that requires snapshots of the solution [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF], no prior knowledge on the solution needs to be known with the PGD. Over the last decade, the space-time separation has led to a plurality of demonstrations that the PGD computes accurate solutions faster than standard incremental solvers. This technique has been applied to a wide range of problems in different branches of physics: in rheology [START_REF] Ammar | Non incremental strategies based on separated representations: Applications in computational rheology[END_REF], thermal analysis [START_REF] Favoretto | Reduced order modeling via pgd for highly transient thermal evolutions in additive manufacturing[END_REF], kinetic theory of complex fluids [START_REF] Ammar | The proper generalized decomposition: a powerful tool for model reduction[END_REF], and even in quantum chemistry [START_REF] Ammar | Circumventing curse of dimensionality in the solution of highly multidimensional models encountered in quantum mechanics using meshfree finite sums decomposition[END_REF]. From a mathematical standpoint, several works [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[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] proved convergence results of the PGD algorithm with the space-time separation for problems involving elliptic operators. In elastodynamics, the space-time separation has been used to compute vibration responses up to the medium frequency domain [START_REF] Ladevèze | PGD in linear and nonlinear computational solid mechanics[END_REF][START_REF] Shirafkan | Quasistatic analysis of elastoplastic structures by the proper generalized decomposition in a space-time approach[END_REF]. However, difficulties have been experienced when the problem equations involve hyperbolic operators [START_REF] Ammar | The proper generalized decomposition: a powerful tool for model reduction[END_REF], which is typically the case with a wave equation. Therefore very few works have addressed the topic of wave propagation with the PGD. In [START_REF] Boucinha | Ideal minimal residual-based proper generalized decomposition for non-symmetric multi-field models -application to transient elastodynamics in space-time domain[END_REF][START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF] the authors presented a space-time separation within a TDG framework to remove spurious oscillations associated to transient wave propagation simulations. They applied their method to a 2D transient elastodynamic problem with a rectangular geometry, and they reported poor convergence results when using a standard PGD algorithm. To improve the numerical performances in terms of storage gain, the authors proposed a recompactification strategy to compute a lower rank approximation of the solution. However, this technique requires an experimented analyst to select the so-called auxiliary rank parameter of the method. In addition, difficulties are reported to compute the solution faster than with a standard TDG technique. As mentioned in their work, the computational cost associated to their method could be improved by considering a complete variables separation strategy, both in space and in time. However, to our knowledge, this approach has not been investigated. In this context, the aim of this paper is not proposing an algorithm making radically faster or better, the existing algorithms already perform well, and from that point of view the one here proposed reduces storage needs and speedup calculations. Here the main and crucial ingredient is the space separation, because nowadays engineered materials may involve thin coatings, geometries sometimes composed of extremely fine plies, and then when solutions are rich throughout the also rich domain thickness, 3D discretizations are sometimes unattainable. This is particularly the case for laser shock processes on composite laminate targets [START_REF] Ghrib | Generation of controlled delaminations in composites using symmetrical laser shock configuration[END_REF][START_REF] Ehrhart | Development of a laser shock adhesion test for the assessment of weak adhesive bonded cfrp structures[END_REF][START_REF] Ecault | Laser shock adhesion test numerical optimization for composite bonding assessment[END_REF]. In this study, we combine standard integrators with an original space separated representation for empowering standard procedures. The proposed methodology is here applied to 2D domains, with stringent loading configurations for the space separation, to analyze its numerical behavior. Therefore, the results do not reflect the entire potential of the method for reducing the CPU time and the storage needs for computationally intensive problems. Future research will be devoted to 3D transient elastodynamic problems in composite laminates.

The paper is organized as follows. In Section 2 we present a time adaptive space separation strategy. We discuss the main features of this approach and we highlight its limitations. We present an implementation of this formulation to the 2D scalar wave equation, but the same methodology can be employed for the 2D elastodynamics equations. Three commonly employed time integration schemes for wave propagation analyses are implemented: the central difference method (explicit, conditionally stable), the Newmark trapezoidal rule and the Bathe method (implicit, unconditionally stable). In Section 3 we explore the numerical behavior of the proposed PGD algorithm on several 2D problems involving radial waves, multiple reflections, propagation on long time analyses, and wave conversions. We compare the obtained results with the exact solution when available, and to numerical solutions computed with a standard finite element method using the same time integration scheme. Finally we discuss the numerical performances of the proposed approach in terms of storage gain and computation time.

Formulation of the method

In this section we introduce a new PGD formulation adapted to transient wave propagation problems in separable domains. The main idea is to perform a space separation at each time increment with an adaptive number of modes. Then we compare the advantages and the limitations of this approach with respect to PGD methods based on a conventional space-time separation. Finally we present the implementation of this formulation to the 2D scalar wave equation.

Space variables separation at each time increment

It is well known [START_REF] Tveito | Introduction to Partial Differential Equations: A Computational Approach[END_REF] that the solution of the 1D wave equation ü -c 2 u ′′ = 0 is under the form u(x, t) = f (xct) + g(x + ct). This result evidences that space and time variables are linked into phase variables x ± ct to describe the propagation of a disturbance. This observation could partially explain the difficulties encountered to implement standard PGD algorithms with a space-time separation for wave propagation problems [START_REF] Ammar | The proper generalized decomposition: a powerful tool for model reduction[END_REF]. Instead, we propose to implicitly separate some phase variables related to each spatial direction. This can be done by performing a standard space variables separation, but at each time increment t k , yielding in the 2D case to the following separated form representation:

u(x, y, t k ) ≈ M(t k ) ∑ i=1 φ i (x, t k )ψ i (y, t k ) (2)
The fields φ i (x, t) and ψ i (y, t) are referred as the ith PGD modes in xand y-directions, respectively. An important feature of (2) is the adaptive number of modes M(t) used to increase the enrichment of the numerical solution when needed, depending on the complexity of the wave field at a given instant. This number is expected to evolve as the waves travel within the medium, reflect at the boundaries, and exhibit sharp spatial variations, as discussed later in Section 2.2. This separated form representation, called time adaptive space separation in the following, is built up on the fly, directly from the considered wave equation, with a standard PGD algorithm based on a space separation at each time increment (see Section 2.4). Any time integration scheme compatible with the standard finite element method framework can be used. The PGD solution may then be computed faster than its standard FEM solution counterpart (see Section 3.1), by solving multiple 1D problems in space rather than a single 2D problem, and with less storage needs, as explained in Section 2.3.

This formulation is applicable when the spatial domain Ω is separable [START_REF] Quaranta | A new hybrid explicit/implicit in-plane-out-of-plane separated representation for the solution of dynamic problems defined in plate-like domains[END_REF]. In 2D, it limits the study to rectangular domains Ω = Ω x × Ω y , and in 3D to domains with an extrusion direction Ω = Ω x y × Ω z (plate, cylinder etc.). In fact, more complex geometries can be considered by using adequate geometrical transformations [START_REF] Ghnatios | Spurious-free interpolations for non-intrusive PGD-based parametric solutions: Application to composites forming processes[END_REF]. It is also important to mention that this space separation cannot exactly represent a wave front as simple as a circle with a finite number of terms in [START_REF] Bolis | Physical approach to adhesion testing using laser-driven shock waves[END_REF]. Hence, at first sight, radial waves seem difficult to be captured with this approach. Yet the numerical results demonstrate that the proposed method is actually effective in the context of numerical approximations for capturing complicated wave fields, including radial waves, with a moderate number of modes (typically 20 in the studied use cases). This can be explained by thinking of the proposed PGD algorithm as performing a Truncated Singular Values Decomposition (TSVD) of the wave field discretized in space at each time increment. It is well known that the TSVD provides accurate low rank approximations in many situations [START_REF] Trefethen | Numerical Linear Algebra[END_REF]. The results presented in Section 3 suggest that the proposed method is capable of inheriting this good property, although it is known that the PGD usually requires higher ranks M than the TSVD [START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF].

Time adaptive enrichment

In this work we emphasize the importance of using an adaptive number of modes M(t) through the simulation, instead of selecting a fixed rank M before the computation. The objective of this feature is to enrich the solution through time when needed, as the waves are initiated by an imposed loading, travel in the medium or interact with the boundaries.

First, similarly to the TSVD, it can be expected that the sharper the wave field is, for instance because of an imposed loading, the more modes will be needed to approximate the solution. This is illustrated on simple examples on Fig. 1a to c. M = 25 terms are necessary with the TSVD to approximate the surface displayed on Fig. 1a with a relative error of 10 -10 , while only M = 4 terms are sufficient to approximate the smoother surface displayed on Fig. 1b with the same relative error.

Second, the propagative behavior of the solution is also expected to affect the number of terms in the sum [START_REF] Bolis | Physical approach to adhesion testing using laser-driven shock waves[END_REF]. More precisely, it can be anticipated that M(t) will increase as the waves expand in the domain. This can be illustrated through the following fictive situation. Consider the diagonal propagation of a linear wave front in a square domain, see Fig. 1d, uniformly discretized in space. The coordinates of the nodes are (x i , y j ) and the discretized wave field at time step t k = k∆t is represented by the matrix with general coefficient U (k) i j = u(x i , y j , t k ). Suppose that the wave front is initiated at time t 1 , so that u(x 1 , y 1 , t 1 ) = 1, and that it diagonally propagates towards the bottom right corner (x n , y n ) until time t n , without interacting with the boundaries (ideal absorption). At each new time increment the wave front expands and reaches new nodes, so that the discretized wave field at any time t k ≤ t n is under the form:

U (k) = [ Ĩk 0 k,n-k 0 n-k,k 0 n-k,n-k ] ( 3
)
where Ĩk is a unit anti-diagonal square matrix of size k. It is straightforward to notice that Ĩk Ĩ T k = Ĩ T k Ĩk = I k , with I k the identity matrix of size k. It follows that the rank of the discretized wave field matrix linearly increases with time as the length of the wave front expands when traveling:

rank(U (k) ) = k for t k ≤ t n (4)
Of course, this fictive example is a worst case scenario in the sense that no low rank approximation can be used to accurately represent this extremely sharp wave front propagation (sudden 0 to 1 discontinuity at each time step). Indeed, all the singular values of U (k) i j are equal to 1 hence there is no truncation order characterized by a jump in the singular values, as observed with the examples presented in Fig. 1c. In practice, it is observed that M(t) generally increases as the waves expand, but up to acceptable values leading to a storage gain with respect to conventional methods.

In light of this discussion, it becomes clearer that an adaptive number of modes is necessary for transient wave propagation problems. On one hand, if the number of modes is not sufficient to describe the solution at a given instant, the numerical errors will accumulate through the simulation and will severely compromise the accuracy. On the other hand, if the number of modes is too high the computation of the PGD approximation will not be effective in terms of computation time and memory capacities.

Storage gain analysis

The storage gain, associated to a PGD formulation with respect to a standard time incremental finite element formulation, can be evaluated as the ratio of the amount of data to write the solution on disk. To do so, consider a rectangular spatial domain with N x and N y discretization points along each direction, and a number N t of output time steps to write the solution. With a time incremental finite element approach, N x N y nodal values of the wave field must be stored at each time step. Therefore the total amount of data to be stored is N x N y N t . On the other hand, with the space-time separation (1) the total amount of data to be stored is M(N x N y + N t ), and it is M(N x + N y )N t with the time adaptive space separation [START_REF] Bolis | Physical approach to adhesion testing using laser-driven shock waves[END_REF], where M is the mean number of modes over the selected output time steps. Three cases illustrate the storage gain associated to the space-time separation: storage gain with space-time separation

≈ ⎧ ⎪ ⎨ ⎪ ⎩ N /M if N x = N y = N t = N N t /M if N x = N y = N and N t ≪ N N 2 /M if N x = N y = N and N t ≫ N (5)
Regarding the time adaptive space separation presented in this paper, the number of time steps N t cancels out in any case, hence the storage gain is given by: storage gain with time adaptive space separation

= (N /2)/M if N x = N y = N (6)
An important storage gain can be obtained with both approaches if the PGD algorithm successfully computes a low rank approximation of the solution. Eq. [START_REF] Ehrhart | Development of a laser shock adhesion test for the assessment of weak adhesive bonded cfrp structures[END_REF] shows that the storage gain may increase with the space-time separation when considering longer time analyses, which is not the case with the time adaptive space separation. However, it is not guaranteed that a low rank approximation can be computed with the space-time separation for transient wave propagation problems, as indicated in Section 2.1 and as reported in [START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF]. On the contrary, the numerical results presented in Section 3 show that the time adaptive space separation effectively computes low rank approximations, and reveal that the mean number of modes M stabilizes when considering longer time analyses. It is also found that M stabilizes when refining the mesh, which indicates that more accurate solutions are computed with more storage gain.

Implementation on the 2D scalar wave equation

The propagation of a scalar wave in a 2D rectangular domain Ω with free contour Γ and null initial conditions is governed by the following equations [START_REF] Tveito | Introduction to Partial Differential Equations: A Computational Approach[END_REF]:

⎧ ⎪ ⎨ ⎪ ⎩ 1 c 2 ü -∆u = f for (x, y, t) ∈ Ω ×]0, T ] ∂u ∂n = 0 for (x, y, t) ∈ Γ ×]0, T ] u(x, y, 0) = 0, u(x, y, 0) = 0 for (x, y) ∈ Ω (7)
where u is the field variable, c is the wave velocity (assumed constant for simplicity), and f is a source term with a separated form representation f (x, y, t) = g(x, t)h(y, t).

The first step is to get the weak formulation of this problem. To do so, we multiply the first equation of ( 7) by an arbitrary test field u * and we use Green's first identity associated with the second equation of [START_REF] Ecault | Laser shock adhesion test numerical optimization for composite bonding assessment[END_REF], which yields:

∫ Ω u * ü + ∫ Ω c 2 ∇u * ∇u = ∫ Ω c 2 u * f (8)
The second step is to discretize equation [START_REF] Cuq-Lelandais | Spallation generated by femtosecond laser driven shocks in thin metallic targets[END_REF] in time with a time integration scheme. For the purpose of clarity we consider here the central difference method with a constant time step ∆t. The time increments are noted t k = k∆t.

Assuming that the solutions at times t k-1 and t k are known, the equation to be solved at t k+1 is:

∫ Ω u * u k+1 -2u k + u k-1 ∆t 2 + ∫ Ω c 2 ∇u * ∇u k = ∫ Ω c 2 u * f k (9)
Then a standard PGD algorithm with a space separation is used on Eq. ( 9) to build the separated form representation of the unknown field u k+1 . This latter is ultimately approximated under the form:

u k+1 (x, y) ≈ M k+1 ∑ i=1 φ k+1 i (x)ψ k+1 i (y) (10)
where a separate finite element discretization is used for each spatial direction to compute the PGD modes. More precisely, by noting P and Q the shape functions used for each spatial direction (linear 1D elements in this study), the finite element discretization in each spatial direction at time t k+1 reads:

φ k+1 i (x) ≈ P(x) T Φ k+1 i (11) ψ k+1 i (y) ≈ Q(y) T Ψ k+1 i (12)
where Φ k+1 i and Ψ k+1 i are the vectors containing the nodal values of the spatial modes φ i and ψ i at time t k+1 , respectively.

Iterative enrichment procedure

The PGD algorithm computes the approximation [START_REF] Bathe | Finite Element Procedures[END_REF] iteratively, by adding new pairs (φ k+1 i , ψ k+1 i ) to the solution within an enrichment loop. The number of modes M k+1 is automatically selected with an enrichment criterion. Each new pair is computed with a fixed point strategy within a stagnation loop. In practice these different sequences are as follows. Assuming that the first n -1 terms of the sum [START_REF] Bathe | Finite Element Procedures[END_REF] are known, the unknown field is enriched with a new pair:

u k+1 n = u k+1 n-1 + φψ (13)
The problem becomes to find a new pair (φ, ψ) satisfying the weak formulation. By injecting [START_REF] Patera | A spectral element method for fluid dynamics: Laminar flow in a channel expansion[END_REF] in [START_REF] Deraemaeker | Dispersion and pollution of the fem solution for the Helmholtz equation in one, two and three dimensions[END_REF], and by moving all the known terms to the right hand side of the equation, the new pair must satisfy:

∫ Ω u * φψ = - ∫ Ω u * u k+1 n-1 + 2 ∫ Ω u * u k - ∫ Ω u * u k-1 - ∫ Ω (c∆t) 2 ∇u * ∇u k + ∫ Ω (c∆t) 2 u * f k (14)
This nonlinear problem is solved with a fixed point method. First the field ψ is assumed known and the arbitrary test field is taken under the form u * = φ * ψ. Each integral over the plate domain is separated as a product of 1D integrals and the nodal approximations ( 11)-( 12) are used. It results in a 1D problem in φ consisting in the following linear system to be solved:

[ n A ∑ i=1 ( Ψ T A i 2 Ψ ) A i 1 ] Φ = n B ∑ i=1 B i 1 (B i 2 ) T Ψ (15)
where Φ and Ψ are the vectors of the nodal values of φ (unknown) and ψ (known), respectively. The operators involved in this equation are given in the Appendix. The same form of linear system is obtained when considering another time integration scheme or for the 2D elastodynamics equations. Then, φ is assumed known with the previously calculated approximation and the test field is taken under the form u * = φψ * . The equations are then manipulated in the exact same manner, leading to a 1D problem in ψ consisting in the following linear system to be solved:

[ n A ∑ i=1 ( Φ T A i 1 Φ ) A i 2 ] Ψ = n B ∑ i=1 B i 2 (B i 1 ) T Φ (16)
This alternative resolution in φ and in ψ is repeated until convergence. More precisely, by denoting (φ (s) , ψ (s) ) the pair calculated at the iteration s of the fixed point method, the stagnation criterion reads:

[ ∫ Ω |φ (s) ψ (s) -φ (s-1) ψ (s-1) | 2 ∫ Ω |φ (s-1) ψ (s-1) | 2 ] 1/2 < ϵ ( 17
)
where ϵ is a threshold value selected by the user. The enrichment loop over the number of pairs n is ultimately stopped when the following enrichment criterion is reached:

[ ∫ Ω |u k+1 n -u k+1 n-1 | 2 ∫ Ω |u k+1 1 | 2 ] 1/2 < ε ( 18
)
where ε is another threshold value selected by the user. This criterion leads to a different length of the sum [START_REF] Bathe | Finite Element Procedures[END_REF] depending on the complexity of the solution at the considered time increment. The convergence can be evaluated by using different error indicators [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] Nadal | A separated representation of an error indicator for the mesh refinement process under the proper generalized decomposition framework[END_REF][START_REF] Alfaro | An error estimator for real-time simulators based on model order réduction[END_REF]. However, in general, the one based on the stagnation works reasonably in most of the applications. The influence of the criteria ( 17)-( 18) on the accuracy and the numerical behavior of the PGD solution is investigated in Section 3.

Time integration schemes for wave propagation analyses

The implementation of the proposed method has been introduced in Section 2.4 with the central difference method, which is the second order accurate explicit scheme with the largest stability condition when combined with the finite element method [START_REF] Bathe | Finite Element Procedures[END_REF]. The advantage of this scheme, when combined with the finite element method, is that the stiffness matrix does not need to be factorized and the mass matrix can be lumped, which results in a diagonal linear system to be solved at each time step [START_REF] Wu | Lumped mass matrix in explicit finite element method for transient dynamics of elasticity[END_REF]. However, as reported in Section 3.5, this time integration scheme is found to be less relevant in the proposed PGD framework. Indeed, the left hand side operators ( 15)-( 16) are computed at each iteration within the enrichment and stagnation loops, while the left hand side operator (the mass matrix) is computed once and for all with the finite element method, before the time incremental resolution.

On the contrary, the proposed PGD framework becomes more efficient when implicit time integration schemes are employed, which involve consistent matrices to be inverted at each time step. The alternative resolution of 1D problems with the PGD may then become more efficient than a single resolution of a 2D problem with the FEM. Therefore, we propose to explore the performances of the proposed formulation with two widely used implicit time integration schemes for wave propagation analyses: the Newmark trapezoidal rule [START_REF] Bathe | Finite Element Procedures[END_REF] and the Bathe method [START_REF] Noh | Performance of an implicit time integration scheme in the analysis of wave propagations[END_REF]. These schemes are defined as follows:

• In the Newmark trapezoidal rule the equations used are:

uk+1 = uk + [ (1 -δ) ük + δ ük+1 ] ∆t (19)
u k+1 = u k + uk ∆t + [( 1 2 -α ) ük + α ük+1 ] ∆t 2 (20)
with δ = 0.5 and α = 0.25. • In the Bathe method the following equations are employed:

uk+1/2 = uk + ( ük + ük+1/2 ) ∆t 4 ( 21
)
u k+1/2 = u k + ( uk + uk+1/2 ) ∆t 4 (22) uk+1 = 1 ∆t u k - 4 ∆t u k+1/2 + 3 ∆t u k+1 (23) ük+1
= 1 ∆t uk - 4 ∆t uk+1/2 + 3 ∆t uk+1 (24)

Numerical results

In this section, we present the performances of the time adaptive space separation for several transient wave propagation problems. First, we consider a pre-stressed membrane with a concentrated impact loading in its center. We study the convergence properties of the PGD solution with comparisons to the exact solution and to a standard 2D finite element solution. Three different time integration schemes are tested. Then, we solve a water ripple propagation problem and we focus on long time analyses to describe the numerical behavior of the PGD solution. Finally, we study a 2D transient elastodynamic problem representative of a laser shock on a metallic target. We study in this last example the effect of varying the concentration of the loading, both in space and in time, on the performances of the proposed approach.

Standard 2D FEM solution counterpart of the PGD solution

The objective of this numerical study is to demonstrate that the proposed PGD formulation computes as accurate transient wave propagation solutions as the standard finite element approach, but faster and with less storage needs. In order to assess the influence of the space separation, we systematically compare the PGD solution to its socalled FEM solution counterpart. As described in Section 2.4, we use a finite element discretization for each spatial direction with 2-node linear elements to compute the PGD solution. Therefore, the FEM solution counterpart is computed with standard 4-node quadrilateral elements, with the same total number of degrees of freedom. As an example, the finite element discretization associated to the scalar wave equation [START_REF] Ecault | Laser shock adhesion test numerical optimization for composite bonding assessment[END_REF] with f = 0 is given by [START_REF] Bathe | Finite Element Procedures[END_REF]:

M Ü + c 2 K U = 0 ( 25
)
where U is the vector containing the nodal values of the solution, M and K are the mass and stiffness matrices, and the contribution to these matrices of element (e) with domain Ω (e) is:

M (e) = ∫ Ω (e)
(H (e) ) T H (e) ( 26) where H (e) is the element interpolation matrix. The matrices M and K are computed from the elementary matrices M (e) and K (e) with a standard summation process. Eq. ( 25) is then discretized in time with the same time integration scheme used to compute the PGD solution.

H (e) = ∫ Ω (e) (∇ H (e) ) T (∇ H (e) ) (27)

Impacted pre-stressed membrane

We consider a pre-stressed membrane impacted in its center, see configuration on Fig. 2, for which the transverse displacement is governed by the following set of equations:

⎧ ⎪ ⎨ ⎪ ⎩ 1 c 2 ü -∆u = f (0, 0, t)
for (x, y, t) ∈ Ω ×]0, T ] u(x, y, t) = 0 for (x, y, t) ∈ Γ ×]0, T ] u(x, y, 0) = 0, u(x, y, 0) = 0 for (x, y) ∈ Ω [START_REF] Shirafkan | Quasistatic analysis of elastoplastic structures by the proper generalized decomposition in a space-time approach[END_REF] with f (0, 0, t) = 4(1 -(2t -1) 2 )H (1 -t) where H (t) is the unit step function, and we select a unit wave velocity c = 1 for convenience. The duration of the time analysis is set to T = 13, namely before the first reflection, to compare the numerical solutions with the exact one given in [START_REF] Yue | Dispersion-reducing finite elements for transient acoustics[END_REF]. Only the domain [0, 15] × [0, 15] is considered due to symmetry. A uniform discretization of the domain is used so that ∆x = ∆y.

This transient wave propagation problem has been extensively studied in the work of Bathe et al. [START_REF] Noh | Performance of an implicit time integration scheme in the analysis of wave propagations[END_REF][START_REF] Ham | A finite element method enriched for wave propagation problems[END_REF][START_REF] Noh | An explicit time integration scheme for the analysis of wave propagations[END_REF][START_REF] Kim | Transient implicit wave propagation dynamics with overlapping finite elements[END_REF] to investigate the numerical behavior of various time integration schemes, and especially to demonstrate the good filtering properties of the Bathe method (see Section 2.6). The sudden impact, concentrated at the center, generates radial waves with sharp displacement variations, resulting in spurious numerical oscillations. This use case is even more stringent within the PGD framework because of the low rank approximation challenge evoked in Section 2.2 to capture circular structures with a space separation.

We compare the numerical results with three different time integration schemes when refining the mesh: the central difference method (CDM), the Newmark trapezoidal rule (TR) and Bathe method. The time step ∆t is constant through the simulation and selected according to each time integration scheme to reduce at best the dissipation and dispersion errors [START_REF] Noh | Performance of an implicit time integration scheme in the analysis of wave propagations[END_REF]. This is done by setting the so-called CFL number defined by CFL = c∆t/∆x to the following recommended values: CFL = 1 for CDM, CFL = 0.65 for TR, and CFL = 1 for Bathe method.

Fig. 3a presents the numerical results obtained with the finite element method. As reported in [START_REF] Noh | An explicit time integration scheme for the analysis of wave propagations[END_REF] the central difference method is not effective on this problem for computing accurate solutions: numerous spurious oscillations are visible. As expected, the accuracy is even worse when considering oblique propagation angles with respect to

Table 1

Relative error with respect to the exact solution and discrepancy between PGD/FEM solutions. Indicated values are computed with a 2-norm on the nodal points the 251 by 251 mesh at the end of the simulation (t = 13).

Error FEM/exact

Error PGD/exact Discrepancy PGD/FEM CDM (CFL = 0.7) 7.3% 6.7% 5.7% TR (CFL = 0.65) 5.8% 5.9% 2.0% Bathe (CFL = 1) 3.1% 3.3% 0.53% the mesh orientation. The implicit TR scheme performs better than CDM but spurious oscillations are still visible, especially at the propagation angle θ = 45 • . The Bathe method removes these spurious oscillations, even with the coarse 101 by 101 mesh. However, it should be noted that when similar time increments are used, implicit schemes are more expensive than explicit ones in terms of memory and computation time. Also, spurious oscillations can be damped with explicit schemes with specific methods, such as the bulk viscosity method and the Tchamwa-Wielgosz method [START_REF] Maheo | Numerical damping of spurious oscillations: a comparison between the bulk viscosity method and the explicit dissipative tchamwa-Wielgosz scheme[END_REF]. Fig. 3b presents the numerical results obtained with the proposed PGD formulation. The stagnation and enrichment criteria ( 17)-( 18) are set equal to ϵ = ε = 10 -3 . The FEM recommended values of the CFL are used for the trapezoidal rule (0.65) and Bathe method (1), but the PGD solution does not converge with the central difference method if the CFL number is not reduced (down to 0.7 with the finest mesh). Interestingly, with the finest mesh it is observed that the spurious oscillations are rather concentrated at low propagation angles with the PGD while it is the opposite with the FEM. Fig. 4 compares the PGD, FEM and exact solutions at the propagation angles θ = 0 • and θ = 45 • . The PGD solution provides the same level of error than its FEM counterpart (see Table 1). For instance, with the 251 by 251 mesh, and with Bathe method, the relative error with respect to the exact solution at the end of the simulation is around 3% for both the PGD and the FEM solutions. The discrepancy between the PGD solution and its FEM counterpart is less than 1%.

Fig. 5a displays the evolution of the number of modes M(t) of the PGD solutions with the different time integration schemes with the 251 by 251 mesh. For each time integration scheme, M(t) monotonically increases, but up to different values (57 with CDM, 75 with TR and 32 with Bathe method). This observation agrees with the fictive situation studied in Section 2.2: M(t) increases as the waves expand in the domain. It is also observed that the less the numerical solution exhibits spurious oscillations, the lesser number of modes are required to represent the solution with a separated form. Consequently, the Bathe method provides the lowest number of modes in the PGD solution. Fig. 5b displays the evolution of the mean number of modes M when the mesh is refined with the three time integration schemes. It is observed that M stabilizes with Bathe method when the mesh is fine enough (below 20), indicating that a limited number of modes is actually sufficient to represent the solution. It follows that more accurate solutions, obtained by refining the mesh, are computed with a higher storage gain. The storage gain reaches 6.7 with Bathe method with the 251 by 251 mesh (see Fig. 5c).

Water ripple propagation

We consider here the propagation of a ripple in a rectangular pond of water and we focus on the numerical behavior of the PGD solution on long time analyses, with multiple wave reflections on the boundaries. This situation is described by imposing an initial displacement profile (see Fig. 6) and by considering free boundaries:

⎧ ⎪ ⎨ ⎪ ⎩ 1 c 2 ü -∆u = 0 for (x, y, t) ∈ Ω ×]0, T ] ∂u ∂n = 0 for (x, y, t) ∈ Γ ×]0, T ] u(x, y, 0) = e -2(x 2 +y 2 ) , u(x, y, 0) = 0 for (x, y) ∈ Ω (29)
The spatial domain is the same as with the previous use case, the wave velocity is also set to c = 1, but we consider a time analysis duration of T = 6L/c = 180, giving rise to multiple wave reflections and crossings. A preliminary mesh convergence study showed that a 201 by 201 mesh with the Newmark trapezoidal rule (CFL = 0.65) is sufficient to provide numerical results with acceptable accuracy. This configuration does not generate spurious oscillations as with the previous use case. Fig. 7 compares the wave fields computed by the PGD and its FEM counterpart at time t = 13, namely just before the first reflection. The solution consists in a donut-shaped ripple (positive front followed by a negative front) propagating from the center of the pond to the boundaries and reflecting multiple times. Fig. 8a shows the evolution of the number of modes M(t) through the simulation depending on the mesh refinement. The parameters ϵ = ε = 10 -3 are used in the PGD algorithm. For the two coarser meshes, M(t) linearly increases from 0 to 23-25 and then stabilizes around 25. This observation suggests that the maximum number of modes is rather linked to the waves expansion in the domain than to the complexity of the wave field, since numerous wave reflections and crossings occur during this simulation (see Fig. 9). However, it is observed that the finest the mesh is, the sooner the M(t) curve starts increasing again (see Fig. 8a). This numerical behavior is removed by reducing the enrichment criterion ε, as depicted on Fig. 8b andc. However reducing ε too much increases the mean number of modes, which reduces the storage gain. In summary, it can be concluded that around 30 modes are in essence sufficient to represent the solution on the complete time interval, regardless how fine the mesh is (even if it means reducing the enrichment criterion ε).

Laser shock on a metallic target

We now consider a 2D transient elastodynamic problem similar to a laser shock on a metallic target. The objectives are to capture wave conversion phenomena and to study the influence of the space-time concentration of the pressure loading on the numerical behavior of the PGD solution.

The displacement field has two components u(x, y, t) and v(x, y, t) gathered in the displacement field u = (u, v). The studied configuration, illustrated on Fig. 10, consists in a rectangular homogeneous linear elastic media under plane strain deformation and small perturbation, subjected to a disymmetric pressure loading on the portion Γ of the top edge. Due to symmetry, only half of the target is considered by imposing a null displacement v = 0 along the middle axis of the real target. The displacement field of this problem is governed by the following set of equations:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ρ ü + div ( σ ) = 0 for (x, y, t) ∈ Ω ×]0, T ] ∂u ∂n = p(0, y, t) for (x, y, t) ∈ Γ ×]0, T ] v(x, 0, t) = 0 for (x, t) ∈ [0, H ]×]0, T ] u(x, y, 0) = 0, u(x, y, 0) = 0 for (x, y) ∈ Ω (30)
where σ is the Cauchy's stress tensor. The pressure loading is under the separated form p(0, y, t) = f (t)g(y) with the time and spatial profiles depicted on Fig. 10. The simulation parameters are presented in Table 2. The PGD operators of this problem are derived with the methodology presented in the Appendix.

Fig. 11 shows snapshots of the displacement field computed with the proposed PGD method with Newmark trapezoidal rule (CFL = 0.65). It is also observed with this use case that M(t) increases as the two components of the waves fill the domain, and then M(t) stabilizes (below 20). Fig. 12 compares snapshots of u(x, 0, t) and u(0, y, t) computed with the PGD and with its FEM counterpart. The discrepancy between the PGD and the FEM solutions is larger near x = H . Our interpretation is that the free boundary condition at x = H doubles the displacement magnitude, leading to sharper variations of the displacement field that are more difficult to be captured by the PGD solution.

The applied loading in x-direction generates a 2D displacement field propagating in the medium. As expected, the applied pressure pushes Γ downwards, which generates a compressive longitudinal wave propagating in xdirection and a transverse wave propagating in y-direction (left column). The speed of the longitudinal and transverse waves can be roughly evaluated by dividing the distance traveled during a given time interval ∆T . The uncertainty on the velocity estimation is given by 2∆x/∆T , where ∆x is the spatial discretization and 2∆x represents the maximum error on the distance approximation with this technique. With the selected parameters (∆T = 360∆t, ∆x = ∆y = 0.1mm) the confidence interval on the velocity estimations is ±114 m/s. It gives 6200 ± 114 m/s for the longitudinal wave propagating along the x-axis, which agrees well with the theoretical value c L = 6198 m/s. And it gives 3100 ± 114 m/s for the transverse wave propagating along the y-axis, which is also consistent with the theoretical value c T = 3122 m/s.

Then the compressive wave propagating along the x-axis progressively initiates waves propagating in y direction. There are two contributions to this effect. First, the initial disturbance progressively expands diagonally (compare the propagation directions indicated by the red arrows at t = 585 ns and t = 1559 ns). This is due to the conversion of the longitudinal wave propagating along the x-axis into transverse waves propagating along the y-axis. Second, the compression front propagating in x-direction locally pushes the material in y-direction, just as a compressed elastic bar clamped at y = 0 and free at y = r would extend its length (Poisson effect). As depicted on Fig. 12 (top left image), this wave conversion results in a progressive decay of the amplitude of the initial compressive wave propagating along the x-axis.

We now investigate the influence of the pressure loading concentration in space or in time. To do so, we consider pressure loadings with a symmetric half-sine time profile. Fig. 13a shows the evolution of the mean number of modes M by varying the ratio T s /(H/c L ) with r/W = 0.3. The shorter the impact duration is, the higher the mean number of modes. This is because the displacement wave fronts become sharper as the duration of the impact decreases. It results that the storage gain increases as the impact duration increases. Fig. 13b shows the same data by varying the ratio r/W with T s /(H/c L ) = 1. The influence of the space concentration is less noticeable than the time concentration, but it is observed that the mean number of modes increases as r diminishes. This is because the waves become more radial if the impact is concentrated in space hence more difficult to be captured with the space variables separation. With T s /(H/c L ) = 1 and r/W = 0.3 the mean number of modes is 8.6 and the corresponding storage gain is 7.

Computation time per time step

Although the proposed method has been implemented on 2D use cases, with stringent loading configurations for the space separation, the following results are promising for addressing computationally intensive 3D problems in future work. We compare hereunder the calculation time to compute the PGD solutions and their FEM counterparts. The simulations are run with Matlab R2018a on a laptop with a single core Intel i7-5500U 2.40 GHz CPU and 8 GB RAM. We focus on the CPU time per time step in the time incremental resolution. The CPU gain is computed as the ratio between the CPU time per time step to compute the FEM solution and the PGD solution with k = 10.

Table 3 summarizes the CPU times per time step obtained for all the use cases by varying the maximum number of iterations in the stagnation loop (k = 10, 50, ∞). The relative error with respect to the exact solution can be computed for the first use case. It can be seen that k = 10 iterations are sufficient to compute numerical solutions as accurate as with k = ∞. Hence the proposed PGD method speeds up the calculations without decreasing the accuracy, and leads to lighter solutions in terms of storage (see Fig. 5c).

Fig. 14 compares the evolution of the CPU time per time step for the first use case depending on the time integration scheme. With the PGD, it is observed that for all the time integration schemes the CPU time per time step linearly increases with the number of degrees of freedom (NDOFs), which is not the case with the FEM. The exception is with the explicit time integration scheme, namely with the central difference method, as explained in Section 2.6. The PGD is more efficient when using implicit time integration schemes for which consistent matrices are inverted. The CPU time per time step is reduced by 5 with Bathe method and the finest mesh (63,001 DOFs).

Two aspects contribute to the reduction of the CPU time per time step. First, the proposed method lowers the complexity of the resolution by solving multiple 1D problems in space rather than a single 2D problem. Second, the space separation strategy reduces the memory consumption. Indeed, the solutions at previous time steps, involved in the right hand side of Eq. ( 15) due to the time discretization, are stored under separated form, which are lighter than full 2D representations if the number of modes is small enough. In addition, when using implicit time schemes, the contribution of the stiffness matrix to the left hand side operator is also stored under separated form, within the operators B i 1 and B i 2 in Eq. ( 15), of sizes N x × M j and N y × M j , respectively, with j = k -1, k -2.

Conclusion

In this work, we proposed a time adaptive space separation strategy to solve transient wave propagation problems in separable domains. The separated form representation of the solution is computed on the fly with a PGD algorithm ran at each time increment. The numerical results show that the PGD solution is enriched when needed, as the waves are generated by the loading and travel in the domain. The number of modes in the PGD solution stabilizes when long time analyses are considered. The time integration scheme should be carefully selected to avoid spurious high frequency oscillations, leading to an increased number of modes. For all the studied use cases, either on the 2D scalar wave equation or on the 2D elastodynamic equations, the proposed method computes the solution faster and with less storage capacities than the standard finite element method. Besides, the number of modes stabilizes when refining the mesh, indicating that the storage gain increases when richer solutions are computed.

Future work will be conducted in the following directions to improve the performances. First, it could be interesting to use the spectral element method instead of the finite element method for the spatial discretization. The spatial accuracy of the PGD solution would be significantly improved by increasing the order of the interpolants, leading to less spurious oscillations, hence to a lower number of modes. Therefore, the memory consumption and the storage needs would be reduced even more. Besides, the spectral elements could be defined on the Gauss-Lobatto-Legendre (GLL) points to diagonalize the left hand side operator that needs to be inverted at each time step, within the enrichment and stagnation loops. The complexity of the resolution would be reduced and so would be the computation time. Second, although the proposed method has been implemented on 2D use cases only, the numerical behavior is promising for future applications on 3D transient wave propagation problems. Indeed, the crucial ingredient is the space separation, which could speed up calculations associated to rich solutions on nowadays engineered materials, involving thin geometries sometimes composed of extremely fine plies. Either a 2D/1D or a 1D/1D/1D space separation strategy could be considered. Finally, the proposed method could be extended to parametric studies on wave propagation problems, to circumvent the so-called curse of dimensionality when it comes to solve high dimension problems. To be consistent with the proposed time incremental approach, a possible implementation could be to run a standard parametric PGD algorithm at each time step.

A.1. Left hand side operators

Let us note P = ∫ Ω x P P T and Q = ∫ Ω y Q Q T . Now let us introduce the finite element discretization ( 11)-( 12) in the left hand side of Eq. ( 14):

∫ Ω u * φψ = ∫ Ω x φ * φ × ∫ Ω y ψψ = ( (Φ * ) T [∫ Ω x P P T ] Φ ) × ( Ψ T [ ∫ Ω y Q Q T ] Ψ ) = (Φ * ) T ( Ψ T QΨ ) (PΦ) (31)

A.2. Right hand side operators

We proceed similarly with each term of the right hand side of Eq. ( 14). The following property is used to present under compact form the terms associated to u ) ( QY k+1 n-1

) T Ψ (33) +2 ∫ Ω u * u k = (Φ * ) T ( 2PX k ) ( QY k ) T Ψ (34)
-

∫ Ω u * u k-1 = (Φ * ) T ( -PX k-1 ) ( QY k-1 ) T Ψ (35)
+ ∫ Ω (c∆t) 2 u * f k = (Φ * ) T ( (c∆t) 2 ∫ Ω x Pg k ) ( ∫ Ω y Qh k ) T Ψ (36)
Now let us introduce P ,x x = ∫ Ω x P ′ (P ′ ) T and Q ,yy = ∫ Ω y Q ′ (Q ′ ) T . Then the remaining term of the right hand side is: [START_REF] Wu | Lumped mass matrix in explicit finite element method for transient dynamics of elasticity[END_REF]

- ∫ Ω (c∆t) 2 ∇u * ∇u k = (Φ * ) T ( -(c∆t) 2 P ,x x X k ) ( QY k ) T Ψ (37)
+ (Φ * ) T ( -(c∆t) 2 PX k ) ( Q ,yy Y k ) T Ψ (

Fig. 1 .

 1 Fig. 1. Surfaces with (a) sharp or (b) smooth variations (discretization with 75 points in each direction). (c) Relative error between the TSVD approximations and the exact surfaces depending on the truncation order. (d) Linear wave front diagonally propagating in a square domain with absorbing boundaries.

Fig. 2 .

 2 Fig. 2. Pre-stressed membrane with concentrated impact force. Numerical results are displayed on blue area.

Fig. 3 .

 3 Fig. 3. Transverse displacement field computed with (a) FEM and (b) PGD at the end of the simulation. Mesh convergence depending on selected time integration scheme. CDM: Central Difference Method. TR: Newmark Trapezoidal Rule.

Fig. 4 .

 4 Fig. 4. Comparison of the numerical solutions (251 by 251 mesh) with the exact solution at different propagation angles (θ = 0 • and θ = 45 • ) at the end of the simulation (t = 13).

Fig. 5 .

 5 Fig. 5. Evolution of (a) the number of modes M(t) through the simulation, (b) the mean number of modes M depending on the mesh, (c) the storage gain as defined by formula (6) depending on the mesh.

Fig. 6 .

 6 Fig. 6. Initial displacement profile (left) imposed at the center of the domain (right). Numerical results are displayed on the blue area.

Fig. 7 .

 7 Fig. 7. Transverse displacement field at t = 13 computed with the FEM (left) and the PGD (right) by using the trapezoidal rule (CFL = 0.65) and a 201 by 201 mesh.

Fig. 8 .

 8 Fig. 8. Evolution of the number of modes M(t) through the simulation depending on (a) the mesh refinement and (b-c) the enrichment criterion.Table2Simulation parameters.

Fig. 9 .

 9 Fig. 9. Snapshots of the transverse displacement field computed with the PGD by using the trapezoidal rule (CFL = 65) and the 201 by 201 mesh.

Fig. 10 .

 10 Fig. 10. Laser shock on a metallic target. Half of the domain is computed by imposing v = 0 on the left edge (corresponding to the middle axis of the real target).

Fig. 11 .

 11 Fig. 11. Snapshots of the 2D displacement field (Newmark trapezoidal rule with CFL = 0.65). The red arrows represent the propagation direction of the wave fronts. Left column: component u. Right column: component v.

Fig. 12 .

 12 Fig. 12. Snapshots of u(x, 0, t) (left) and u(0, y, t) (right) before the first reflection at x = H (top) and after (bottom).

Fig. 13 .

 13 Fig. 13. Evolution of the mean number of modes and storage gain depending on the loading concentration (a) in time or (b) in space.

Fig. 14 .

 14 Fig. 14. CPU time per time step for the simulations on the first use case (see Section 3.2). PGD solutions are computed with ε = 10 -3 and a limited number of iterations in the stagnation loop (k = 10).

k+1 n- 1 ,

 1 u k and u k-1 : i ] = a T (AB) (BD) T c (32) where B = [b 1 |...|b M ] and D = [d 1 |...|d M ]. Let us use the notation X k = contributions to the right hand side are:-∫ Ω u * u k+1 n-1 = (Φ * ) T ( -PX k+1 n-1

Table 2

 2 Simulation parameters.

	E [GPa]:	70
	ν:	0.33
	ρ [kg/m 3 ]:	2700
	H [mm]:	10
	W [mm]:	15
	T s [µs]:	1.32
	r [mm]:	2
	T [µs]:	3.23
	∆t [ns]:	4.87
	N x :	100
	N y :	150
	ϵ:	10 -4
	ε:	10 -4

Table 3

 3 CPU times required to solve the 2D transient wave propagation problems. For the first use case the percentage under parenthesis indicates the relative error with respect to the exact solution at the end of the simulation. The Bathe method is used for the membrane problem with CFL=1, the Newmark trapezoidal rule is used for the water ripple and laser shock problems with CFL=0.65.

	Use case	NDOFs FEM	PGD			CPU gain a
			CPU time per time step (s)	CPU time per time step (s)
			(number of increments)
				k = 10	k = 50	k = ∞
	Membrane	63,001	1.1 (109) (<4%)	0.22 (<4%)	0.26 (<4%)	0.30 (<4%)	5.0
	Water ripple 40,401	0.31 (1847)	0.14	0.19	0.26	2.2
	Laser shock	30,000	0.46 (306)	0.14	0.16	0.17	3.3
	a

Acknowledgments

This work was supported by a government funding under the project FUI, France MONARQUE.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix

We derive herein the PGD operators associated to the 2D wave propagation problem presented in Section 2.4. The same methodology can be used for deriving the operators associated to the elastodynamic problem studied in Section 3.4. More details regarding the manipulation of a multi-component displacement field u = (u, v) within a PGD framework can be found in [START_REF] Quaranta | A new hybrid explicit/implicit in-plane-out-of-plane separated representation for the solution of dynamic problems defined in plate-like domains[END_REF].


Laurent Bessie

Gurvan Madec

Florent Lyard

Global tidal residual mean circulation: Does it affect a climate OGCM?

A global description of the ocean tidal residual mean circulation (TRMC) is provided and its effect on a global ocean climate model is evaluated. The TRMC is computed from a global barotropic tidal model based on pure hydrodynamical equations and forced by the major tidal constituents. It thus includes both Stokes drift and nonlinear tidal dynamics. The TRMC map exhibits remarkable basin scale structures in the Indonesian archipelago-Australian region and on the Malvinas shelf. Regions of largest transport are characterized by sites of high tidal amplitudes, abrupt topography and western boundaries. Yet transport magnitudes remain at a few tenths of a Sverdrup. This TRMC is introduced as an external forcing in an ocean climate model. At equilibrium, the thermohaline structure and the poleward heat transport of the ocean are almost insensitive to the TRMC. It is concluded that the TRMC is negligible for climate purposes in an ocean climate model.

Introduction

[2] Tides play a major role in the ocean overturning circulation [START_REF] Munk | Abyssal recipes II: Energetics of tidal and wind mixing[END_REF]]. They generate internal waves which dissipate in the interior ocean and this dissipation induces vertical mixing that contributes to sustain the thermohaline circulation [START_REF] Wunsch | Vertical mixing, energy, and the general circulation of the oceans[END_REF]. This process has been incorporated in several climate Ocean General Circulation Model (OGCM) [e.g., [START_REF] Simmons | Tidally driven mixing in a numerical model of the ocean general circulation[END_REF]L. Bessie `res et al., Improved tidally driven mixing in a numerical model of the ocean general circulation, submitted to Ocean Modelling, 2008, hereinafter referred to as Bessie `res et al., submitted manuscript, 2008]. In all cases, the introduction of tidally driven mixing leads to an improvement of the thermohaline structure and the circulation of the simulated ocean. Nevertheless tidal mixing is not the only aspect that can influence the ocean. Though quasi-periodic, tides also generate a mean circulation, the tidal residual mean circulation (TRMC). This circulation is due to both Stokes transport and non-linear tidal dynamic [START_REF] Longuet-Higgins | On the transport of mass by time-varying currents[END_REF][START_REF] Zimmerman | On the Euler-Lagrangian transformation and the Stokes drift in the presence of oscillatory and residual currents[END_REF]. The TRMC is extremely difficult to extract from observations, as it is almost impossible to isolate it from the mean wind and thermohaline driven circulations. Nevertheless, TRMC can be evaluated using numerical models. In coastal oceanography, its importance is well known, in particular its effect on the transport of sediments [START_REF] Hoitink | Flow asymmetry associated with astronomical tides: Implications for the residual transport of sediment[END_REF][START_REF] Me ´nesguen | Observation and modelling of a natural retention structures in the English Channel[END_REF]]. At regional scales, several studies have shown that the TRMC can be the main contributor to the mean circulation [e.g., Le [START_REF] Provost | Tidal spectroscopy of the English Channel with a numerical model[END_REF][START_REF] Sentchev | Residual circulation in the English Channel as a dynamically consistent synthesis of shore-based observations of sea level and currents[END_REF]. Furthermore, [START_REF] Schiller | Effects of explicit tidal forcing in an OGCM on the water-mass structure and circulation in the Indonesian Throughflow region[END_REF] has shown that the TRMC contributes significantly to an inter-basin exchange. He found that it is responsible for a transport of $0.7 Sv from the Pacific Ocean to the Indian Ocean, i.e. 7% of the total transport of the Indonesian Throughflow (ITF). To our knowledge the effects of the TRMC at global scale has never been addressed, and even more the global TRMC itself has never been described. The aim of this paper is to provide such a global description and to test its effect on a global ocean climate model.

[3] Section 2 provides a description of the TRMC at a global scale, computed from a barotropic tidal model. In section 3 the TRMC is introduced in a climate OGCM as an external forcing and the ocean response is investigated. [START_REF] Carre `re | Modelling the barotropic response of the global ocean to atmospheric wind and pressure forcing: Comparisons with observations[END_REF][START_REF] Lyard | Modelling the global ocean tides: Modern insights from FES2004[END_REF]. The finite element space discretisation (triangular P1 Â P1 elements) allows us to adapt the model resolution with respect to the local dynamical constraints (e.g. classical wavelength criterion) and coastal geometry (e.g. local slope criterion). The resolution varies from $10 km along the coasts and in steep topography or shallow water areas, to 200 km in the deep ocean. The domain includes all marginal seas as well as tides under Weddell and Ross Sea ice-sheets. Tidal dissipation is introduced with both a quadratic bottom friction and an internal wave drag dissipation parameterization. The simulation has been forced by the tidal potential of the major tidal constituents, namely M 2 , S 2 , N 2 , K 1 and O 1 . This set of components explains more than 90% of the observed barotropic tidal signal.

Global Tidal Residual Mean Circulation

[5] We have chosen not to use data assimilation in order to maintain the consistency of the pure hydrodynamical equations. This is important as the signal we want to extract (TRMC) is weak compared to the tides themselves. Even without assimilation, the results show a very good global accuracy [START_REF] Lyard | Modelling the global ocean tides: Modern insights from FES2004[END_REF]. Rms differences of 4.0 cm (M 2 ) and 1.2 cm (K 1 ) are found when compared to the best altimetry derived global models, such as GOT00 or FES2004 [START_REF] Lyard | Modelling the global ocean tides: Modern insights from FES2004[END_REF]. The residual mean transport has been computed as the averaged instantaneous model transport over a one year period. The resulting 2D vector field has been decomposed into a potential function and a barotropic stream function (BSF) associated with the divergent and rotational parts, respectively [START_REF] Lynch | Deducing the wind from vorticity and divergence[END_REF]. The transport reconstructed from the divergent potential function is almost zero everywhere, demonstrating that a one year average is long enough to reduce the aliasing effect associated with the presence of several tidal components to a negligible level. The TRMC is thus entirely described by the BSF. Then, the BSF has been interpolated onto a global geographic mesh of 1/8°i.e. a resolution comparable to the smallest finite element used. This choice of a BSF formulation has been motivated by its capacity to preserve transports, especially when interpolating on an OGCM grid with different resolution and bathymetry.

[6] The BSF associated with TRMC is shown in Figure 1a. Its intensity ranges from À0.54 Sv to 0.66 Sv (Positive and negative transports respectively match the clockwise and counter-clockwise sense of currents). Areas of largest transports are located around the Philippines (0.4 Sv), Australia (0.2 Sv) and over the Malvinas shelf (0.3 Sv). All these locations are characterized by: (1) sites of large amplitude M2, K1 and M4 constituents (>75 cm), ( 2) abrupt topography such as continental slopes and (3) western boundaries. In addition, barotropic flows along boundaries induce recirculation gyres as they have to loop back. The global map reveals a set of remarkable cells which circulation loop around these locations. Weaker but still significant cells appear along the Kurill and Aleutian Islands and north of the undersea ice shelf of Weddell Sea. A clockwise cell is also centred around the Salomon and Vanuatu Islands which fills the whole South Pacific Ocean area, but it is quite anecdotal in terms of transport. Barotropic currents associated with the TRMC range between 0 and 14 cm.s À1 , with a typical value of 1 cm.s À1 . The maximal values that exceed 10 cm.s À1 are found in two narrow and shallow places: Magellan straits and Tatar straits (between Sakhalin Island and Russia).

[7] Two of the three largest TRMC cells are located between the Pacific and Indian Oceans, a key position for water mass exchange (cf. Figure 1b). The northern cell is centred on the Philippine and Celebes archipelagos (the clockwise residual mean transport is: from 0.2 Sv up to 0.4 Sv). In this area the largest residual mean transport is found along the Philippine trench (0.4 Sv) with a peak in magnitude in the Luzon Strait (0.5 Sv), which represents 10% of the eastward volume transport estimated by [START_REF] Tian | Observation of Luzon Strait transport[END_REF] in this strait. The transport along the Philippine trench flows southward down to the Celebes Sea where the circulation splits into two branches. A westward branch (0.2 Sv) loops back and re-circulates northward. The second branch (0.2 Sv) follows the eastern route of South Pacific subthermocline waters via the Molucca (0.1 Sv), Halmahera and Seram (0.1 Sv) Seas down to the Banda Sea. At Flores Strait the circulation loops back and re-circulates in the Makassar Strait (0.2 Sv) up to the Luzon strait. The southern cell is centred around the Australian continent with a westward extension into the Indian Ocean (the induced counter-clockwise residual mean transport is: À0.2 Sv). This circulation splits into two branches before entering the ITF. A northern branch (0.1 Sv) circulates north of Papua New Guinea before supplying Lombok and Ombai Straits. The southern branch supplies Timor Passage via Torres Strait (0.1 Sv). Finally, the overall flow pattern in the ITF is only due to the southern cell. This ultimately leads to a transport of $0.25 Sv between Australian northern coasts and the archipelago of Indonesia. The total mean transport in the ITF is $10 Sv [START_REF] Gordon | Oceanography of the Indonesian Seas and their throughflow[END_REF]; therefore the TRMC in the ITF represents about 2.5% of the data-estimated transport, which is lower than the $7% estimated by [START_REF] Schiller | Effects of explicit tidal forcing in an OGCM on the water-mass structure and circulation in the Indonesian Throughflow region[END_REF]. The higher resolution of the tidal T-UGO model (15 km coastal resolution) compared to [START_REF] Schiller | Effects of explicit tidal forcing in an OGCM on the water-mass structure and circulation in the Indonesian Throughflow region[END_REF] model ($50 km) provides a more detailed description of the Indonesian archipelago region. This is particularly apparent in the area between Borneo and the Philippines where the narrow and shallow Strait of Mindoro is open, north of Palawan Passage. This strait also carries the return branch of the tidal residual mean transport (0.4 Sv) associated with the Philippine and Celebes archipelagos. Therefore, in this study the spatial structure of the TRMC is basically different than the one of [START_REF] Schiller | Effects of explicit tidal forcing in an OGCM on the water-mass structure and circulation in the Indonesian Throughflow region[END_REF].

[8] The third main cell of TRMC has a semi-basin extension in the South West Atlantic Ocean. The patterns of the circulation draw a clockwise loop around the Malvinas shelf and the Scotia Ridge, and are correlated with the topography features. The flow is northward, offshore from the Patagonian coast up to the Santos Plateau where it loops back and re-circulates southward along the shelf break (its residual mean transport is: +0.3 Sv) (cf. Figure 1c). Areas with the highest magnitude are found near the Malvinas shelf break (0.3 Sv) and around the Scotia Ridge (0.4 Sv). A peak in magnitude is found downstream of Drake Passage (0.55 Sv). Yet, the tidal residual mean transport remains two to three orders of magnitude smaller than the transport associated with the Western boundary current and the Antarctic Circumpolar Current (136.7 ± 7.8 Sv) [START_REF] Cunningham | Transport and variability of the Antarctic Circumpolar Current in Drake Passage[END_REF].

[9] These results are extremely difficult to validate, as it is almost impossible to identify the TRMC in observations of the mean wind and thermohaline driven circulations [Le [START_REF] Provost | Tidal spectroscopy of the English Channel with a numerical model[END_REF]]. In addition, simple time averaging of the velocity data obtained at a fixed mooring station might lead to a spurious material transport [START_REF] Katsumata | Frictionless generation of a tidal Eulerian residual flow over a sill in a narrow channel[END_REF]. TRMC maps computed from barotropic tidal models are also rare in the literature, except for coastal areas [e.g., Me ´nesguen and Gohin, 2006] or those showing only the Stokes drift [START_REF] Sentchev | Residual circulation in the English Channel as a dynamically consistent synthesis of shore-based observations of sea level and currents[END_REF]. However, an indirect validation can be made for the two components of the TRMC: the Stokes drift and the nonlinear tidal dynamic component. The Stokes drift component is well represented by the model since its accuracy is directly correlated with the accuracy of tidal constituents. The non-linear tidal dynamic component should be also well reproduced. Indeed, the T-UGOm solution for the quarter-diurnal tidal constituents is in good agreement with observations. A rms difference of 0.3 cm is found when compared to bottom pressure records deriving from BPR data base. Since the quarter-diurnal tides are produced by the same non-linear processes as the residual mean current, we have some confidence in the present results.

[10] It must be emphasized that the described TRMC does not take into account of the transport associated with the propagation of internal waves. Nevertheless, such a contribution may be considered as negligible for two reasons. The former is that a residual mean transport associated with the propagation of internal waves would require non-linear effects, the origin of which is different from the one of wave breaking and is much more localized (e.g., wave-topography reflexions or wave-wave interactions). So, internal waves do not necessarily develop such non-linearities, making such a transport hypothetical. The latter reason is that, if such a transport would happen, the fraction of tidal energy that would be dedicated to it would not range more than from about 20 to 30% of the barotropic tidal energy. Indeed, the fraction of tidal energy dedicated to internal waves ranges from about 1/4 to 1/3 of the total tidal energy that is dissipated in the ocean [START_REF] Munk | Abyssal recipes II: Energetics of tidal and wind mixing[END_REF], furthermore only 2/3 of this fraction is likely to radiate away from the generation site [START_REF] St | The role of internal tides in mixing the deep ocean[END_REF]. Therefore, we assume in this study that the baroclinic contribution to TRMC is negligible. Yet, it must be noticed that the baroclinic tidal contribution associated with internal wave breaking has been taken into account in the OGCM via a tidally driven mixing parameterization (Bessie `res et al., submitted manuscript, 2008).

Climate OGCM Response

[11] In order to explore the effects of the TRMC on the large-scale ocean circulation, we have used a climate OGCM derived from NEMO ocean modelling framework [START_REF] Madec | OPA8.1 OGCM reference manual[END_REF]. The model configuration is ORCA2-LIM: on a horizontal 2°Â 2°Mercator tripolar grid; 31 vertical levels (with a 10 m resolution in the upper 100 m, increasing to 500 m near the bottom). The physics and the forcing are described by [START_REF] Iudicone | The global conveyor belt in a Southern Ocean perspective[END_REF]. The physics includes a sea-ice model; a free surface formulation; lateral mixing evaluated along isoneutral surfaces; an eddy-induced velocity; vertical mixing computed by a turbulent kinetic energy closure scheme; a diffusive bottom boundary layer; convection processes parameterized by locally enhancing the vertical diffusivity; and a synchronous time step of 5400 sec. In addition, the model uses a partial-step representation of the bottom topography and tidally driven bottom intensified mixing (Bessie `res et al., submitted manuscript, 2008). At the surface the model is forced by heat and fresh water fluxes, using bulk formulae and monthly climatologies of wind stress, air temperature, humidity, cloud cover and precipitations. A restoring term to sea surface salinity is also added.

[12] As for the eddy induced velocities [START_REF] Gent | Isopycnal mixing in ocean mixing circulation[END_REF], the TRMC is introduced in the model as an external forcing, so that its associated currents are added to the eulerian velocities on the tracer advection. Interpolation consistent with the OGCM grid is required. Directly imposing with tidal residual mean transports would have led to spurious divergent quantities when interpolating, since the OGCM grid has a different spatial resolution and topography. This problem vanishes when the BSF is interpolated. Indeed, potentials are easily conserved when interpolating. As a result, extrema of the previously discussed gyre structures have not been altered, except those around the Kurill and Aleutian Islands whose spatial scales are too close to the model grid spacing. Attention has been paid to preserve the initial potential for each island accounted for in the climate OGCM. In particular, the ITF transport associated with the TRMC is exactly the same in ORCA2-LIM as in T-UGOm. The barotropic currents associated with the TRMC are then reconstructed in the model from the BSF via numerical derivations. In order to quantify the climate OGCM response to TRMC forcing, two experiments have been performed, one with the TRMC forcing (RES) and the other one without (NoRES). An integration over 1500 years was sufficient to reach a steady state i.e. at which time the trends in globally averaged potential temperature and salinity were respectively less than 0.02°C and 0.001 psu per millennia.

[13] The meridional overturning stream function (MSF) diagnosed for RES experiment (Figure 2a) is in agreement with what is expected with this configuration [START_REF] Iudicone | The global conveyor belt in a Southern Ocean perspective[END_REF]Bessie `res et al., submitted manuscript, 2008]. The TRMC barotropic forcing is associated with an overturning circulation which is represented as MSF in Figure 2c. This resulting overturning forcing has the same order of magnitude as the BSF (O(0.1) Sv). Its spatial structure shows a set of meridional cells whose location and regional extension are correlated with the main gyres of the TRMC. The OGCM response to the TRMC forcing is presented in term of MSF difference between RES and NoRES (Figure 2b). Its order of magnitude (O(0.1) Sv) is the same as that of the overturning forcing (Figure 2c). The spatial distribution of the ocean response is not associated with an intensification of the meridional overturning cells. Instead, these structures seem to be located in maximal gradient areas, where small displacements induce strong differences. Therefore, the weak ocean response indicates that the TRMC has only a slight effect on the global meridional overturning circulation.

[14] Quantitative arguments based on TRMC forcing analysis show that the meridional heat transport (MHT) should be negligible. Firstly, barotropic flows along boundaries induce re-circulation gyres as they have to loop back. Gyres are known to be less efficient in conveying heat than the meridional overturning circulation. In the case of the Gulf-Stream, the circulation induced along the western boundary is about O(40) Sv near the surface and the associated heat transport is approximately O(0.1) petawatt (PW). Since the TRMC induced re-circulation gyres are about O(0.1) Sv over the whole water column, the TRMC induced MHT can not be higher than O(0.01) PW. Secondly, the overturning circulation due to the North Atlantic deep waters induces a O(10) Sv meridional mass transport associated with a O(1) PW MHT [START_REF] Ganachaud | Large-scale ocean heat and freshwater transports during the World Ocean Circulation Experiment[END_REF]]. If we consider that the MHT is proportional to the strength of the overturning circulation times the temperature difference across the pycnocline, the O(0.1) Sv meridional mass transport induced by the TRMC should not lead to a MHT higher than O(0.01) PW. Consequently, neither the gyre circulations induced by the TRMC nor the overturning ones are likely to modify the MHT beyond a hundredth of a PW. After computation, the difference in the MHT between the RES and NoRES cases is smaller than a hundredth of a PW, as expected (not shown).

[15] Largest changes in heat content are found in the Atlantic Ocean, south of Greenland, where a difference of À0.2°C is found in the first hundred meters between RES and NoRES. This results from a O(10) km horizontal displacement in a strong temperature gradient area, which does not affect deep water mass formation in Labrador Sea. Almost everywhere else, differences of heat content in the first hundred meters remain below 0.05°C of magnitude. For example, the TRMC transport increase of 0.25 Sv in the ITF leads to a slight raise of the equatorial thermocline by about 0.5 m in the Pacific Ocean, which is proportional to the raise occurring when the ITF is opened/closed in the model [START_REF] Gorgues | Biogeochemical impact of tropical instability waves in the equatorial Pacific[END_REF].

Conclusion

[16] The residual mean circulation induced by tides (TRMC) has been computed at a global scale and its effect on an OGCM has been investigated. The computation of the TRMC has been obtained from an accurate barotropic tidal model. The circulation spreads over large-scale gyres, whose intensities range from À0.54 Sv to 0.66 Sv, which can be large compared to coastal circulation but remains negligible for mean large-scale circulation. These gyres generate an overturning circulation of weak amplitude ($0.2 Sv), which is introduced as an external forcing in an ocean climate model. The effects of this forcing on the global meridional overturning circulation are negligible. Therefore, the TRMC can be neglected in current ocean components of earth system modelling for climate purposes. Concerning the particular case of the ITF, it reinforces the idea that tidally driven mixing remains the major way for tides to affect water-masses properties [START_REF] Koch-Larrouy | On the transformation of Pacific Water into Indonesian Throughflow Water by internal tidal mixing[END_REF].

[17] Nevertheless, local non-linearities may still lead to regional modifications of the ocean state. In this regional scope, a higher resolution TRMC would better account for the bathymetry variations, which are central for resolving the patterns of the TRMC at a local scale. Although the present version of the TRMC remains sufficiently accurate for large-scale studies, it could be improved for coastal and regional studies.

[18] Acknowledgments. This work was supported by the CNRS, the DGA, Mercator-Ocean (project 100017), MERSEA (project SIP3-CT-2003-502885), and IDRIS. The authors would like to thank Rosemary Morrow and Sebastien Masson for fruitful discussions. They express their gratitude in memory of Christian Le Provost, who initiated the study.

[ 4 ]

 4 In order to compute the barotropic TRMC at a global scale, the Toulouse Unstructured Grid Ocean model (T-UGOm) has been used to generate the global ocean tides. The external mode of T-UGOm is identical to MOG2D-G non-linear hydrodynamical barotropic model [Carre `re and

Figure 1 .

 1 Figure 1. (a) Global barotropic stream function (BSF) associated with the tidal residual mean transport. Units in Sverdrup. Positive and negative transports respectively match the clockwise and counter-clockwise sense of currents. Zooms have been added (b) in the Indonesian archipelago and (c) in the Malvinas shelf area, where TRMC induced transports are the most significant.

Figure 2 .

 2 Figure 2. Meridional overturning stream functions (MSFs) for the global ocean from the ocean general circulation model ORCA2-LIM, after 1500 years of run. Units in Sverdrup. (a) Numerical experiment RES including the forcing induced by the tidal residual mean circulation (TRMC). The contour interval is 2 Sv. (b) Differences between the experiments RES and NoRES (NoRES is a control run without TRMC). The contour interval is 0.1 Sv. (c) Meridional overturning circulation induced by the TRMC alone. The contour interval is 0.1 Sv.