Anthony Giacoma

David Dureisseix

Anthony Gravouil

Michel Rochette

Toward an optimal a priori reduced basis strategy for frictional contact problems with LATIN solver

Keywords: 74M10, 74M15, 65M60, 74S05 Frictional contact, LATIN method, Reduced Order Model, Proper Generalized Decomposition (PGD)

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Finite element simulations involving frictional contact dare the computational mechanics community. In spite of sophisticated algorithms, such simulations can lead to strong numerical difficulties and large computational cost. Contact mechanics is a broad topic involving several geometrical, optimization and numerical aspects relying on strong mathematical foundations. A great coverage of computational contact mechanics with all its basic ingredients is available in [START_REF] Yastrebov | Numerical methods in contact mechanics[END_REF]: from contact detection to implementation of contact algorithms in a finite element software.

Frictional contact leads to a strong non-linear problem on the contacting interface. Indeed, such phenomenon is characterized by non-smooth constitutive laws (for both normal and tangential components). From various computational fields and over the last decade, numerous numerical methods were proposed to solve such problems. Generally speaking, contact problems can be formulated as an optimization problem subjected to constraints. Then, the Lagrange multiplier method can be used. For frictionless cases, such a formulation can be casted into a Linear Complementarity Problem and solved with dedicated approaches (Lemke's algorithm, projected successive over-relaxation, active-set methods...). For frictional cases, [START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF] proposes a method based on a discretization of the Coulomb's cone affecting toughly the precision of the solution. Augmented Lagrangian approaches are the most widespread formulations and many methods were proposed to solve both frictionless and frictional problems with these approaches [START_REF] Fortin | Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems[END_REF][START_REF] Simo | An augmented lagrangian treatment of contact problems involving friction[END_REF][START_REF] Laursen | Formulation and treatment of frictional contact problems using finite elements[END_REF][START_REF] Wriggers | Finite element algorithms for contact problems[END_REF][START_REF] Heinstein | An algorithm for the matrix-free solution of quasistatic frictional contact problems[END_REF][START_REF] Alart | A scalable multiscale latin method adapted to nonsmooth discrete media[END_REF][START_REF] Bussetta | The adapted augmented lagrangian method: a new method for the resolution of the mechanical frictional contact problem[END_REF][START_REF] Pietrzak | Large deformation frictional contact mechanics: continuum formulation and augmented lagrangian treatment[END_REF]. Non-linear Gauss-Seidel solvers are also developed and applied to problems involving multi-contact bodies (rigid or flexible) [START_REF] Jourdan | A gauss-seidel like algorithm to solve frictional contact problems[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF], but they are not usually the most efficient way to solve large systems due to a somehow low convergence rate. Gradient methods [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF][START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF][START_REF] Renouf | Conjugate gradient type algorithms for frictional multi-contact problems: applications to granular materials[END_REF][START_REF] Dostál | Minimizing quadratic functions subject to bound constraints with the rate of convergence and finite termination[END_REF] and generalized Newton's methods [START_REF] Alart | Méthode de newton généralisée en mécanique du contact[END_REF][START_REF] Li | Smoothing newton and quasi-newton methods for mixed complementarity problems[END_REF][START_REF] Curnier | A generalized newton method for contact problems with friction[END_REF] were also adapted to contact problems to reproduce performances of these well-known techniques. The non-smooth behavior of the frictional contact generally prevents to provide a convergence proof of these algorithms. Finally, the bipotential method [START_REF] Saxcé | The bipotential method: a constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF][START_REF] Joli | Uzawa and newton algorithms to solve frictional contact problems within the bi-potential framework[END_REF] provides a formulation and a suited algorithm allowing also some computer cost reduction. In [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF], main results on space discretization scheme and non-linear solver for the solution of frictional contact problems are provided.

Usually and even within a quasi-static context, all these non-linear solvers may lead to a prohibitive computing time for industrial purposes. To address this computational issue, several acceleration methods can be used. First, an acceleration strategy based on multigrid methods were proposed in [START_REF] Alart | Solution of frictional contact problems using ilu and coarse/fine preconditioners[END_REF][START_REF] Lebon | Multigrid methods for unilateral contact problems with friction[END_REF]. It consists in computing cheaply corrections on coarser discretizations (grids) to accelerate the convergence rate of a non-linear solver (e.g. non-smooth Newton method) also called smoother on the finest description of the problem. Second, domain decomposition techniques (FETI-based) were also proposed for frictional contact problems [START_REF] Dostál | A scalable tfeti algorithm for two-dimensional multibody contact problems with friction[END_REF][START_REF] Dostál | A theoretically supported scalable tfeti algorithm for the solution of multibody 3d contact problems with friction[END_REF]. Such techniques aims at partitioning bodies into nonoverlapping sub-domains processed independently. Then, interesting scalability property can be gained. On the other hand and over the last decade, model reduction techniques surge in enthusiasm within the computational mechanics field. For both linear and non-linear problems, impressive accelerations can be obtained. Nevertheless, for frictional contact mechanics, the use of such reduced basis methods seem to be difficult due the non-smoothness of the constitutive laws. As proposed and claimed in [START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF], using a space-time approach allows

g f ext u p Ω ∂1Ω ∂2Ω ∂3Ω t 0 t1 • • • tk • • • T Figure 1: Problem set.
to circumvent this difficulty and benefiting from computational gain of such reduced basis approaches.

Classically, two different approaches are arising from model reduction techniques. Both consist in solving a given large scale reference problem in a reduced size subspace. The difference lies on the method providing the so-called Reduced Basis (RB). The first approach is called a posteriori method and consists firstly in designing a Reduced Order Model (ROM) [START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF][START_REF] Amsallem | Interpolation method for adapting reducedorder models and application to aeroelasticity[END_REF][START_REF] Amsallem | An online method for interpolating linear parametric reduced-order models[END_REF][START_REF] Carlberg | The gnat method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows[END_REF][START_REF] Amsallem | An online method for interpolating linear parametric reduced-order models[END_REF][START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF]. This is achieved by computing a RB with solution snapshots using e.g. Proper Orthogonal Decomposition (POD) or Singular Value Decomposition (SVD). Then, the reference problem and its governing equations are projected into the subspace engendered by a truncation of this RB. The ROM is built during an offline phase and may be expensive. The online phase consists in computing approximated solutions of the ROM. Such a technique aims to provide real-time and on-board computations and address efficiently multiparametric problems. Nevertheless, accuracy of the ROM depends strongly on the relevancy of the initial RB [START_REF] Heyberger | A rational strategy for the resolution of parametrized problems in the PGD framework[END_REF][START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF].

The second approach is the a priori method [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Ryckelynck | On the a priori model reduction: Overview and recent developments[END_REF]. Contrary to the previous one, it does not require solution snapshots. The RB generating the approximation for the solution is computed along the resolution process. The best-known a priori method is the Proper Generalized Decomposition (PGD) [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Ladevèze | The LATIN multiscale computational method and the proper generalized decomposition[END_REF][START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF]. It consists in searching the solution in a separable form reducing the computation cost. Moreover, the expensive offline phase is avoided.

In practice, a posteriori and a priori methods are complementary. The first approach approximates efficiently a given reference problem. The second one provides an efficient tool to adapt or enrich a given RB on-the-fly to generate suited approximations for the solution [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: An adaptive model order reduction for highly nonlinear mechanical problems[END_REF][START_REF] Heyberger | A rational strategy for the resolution of parametrized problems in the PGD framework[END_REF][START_REF] Galland | An adaptive model reduction approach for 3D fatigue crack growth in small scale yielding conditions[END_REF][START_REF] Galland | A global model reduction approach for 3d fatigue crack growth with confined plasticity[END_REF].

In the present article, the authors propose an a priori approach close to the LATIN-PGD approach [START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF][START_REF] Relun | A model reduction technique based on the pgd for elastic-viscoplastic computational analysis[END_REF][START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF] allowing to address frictional contact problems. This new non-linear solver uses ingredients from both the Large Time INcrement method (LATIN) [START_REF] Ladevèze | Nonlinear Computational Structural Methods: New Approaches and Non-Incremental Methods of Calculation[END_REF] and the PGD.

Reference problem

A flexible body is contacting a rigid basement (figure 1). We define Ω, the domain occupied by this solid whose boundary is ∂Ω at time t ∈ [0, T ] ⊂ R. This boundary can be split into three complementary parts:

• ∂ 1 Ω: the part of the boundary where displacements u p are prescribed;

• ∂ 2 Ω: the part of the boundary where external loads f ext are prescribed;

• ∂ 3 Ω: the part of the boundary where contact conditions may occur (potential contact interface).

Strong formulation

Linear assumptions are assessed (small displacements, homogeneous-isotropiclinear-elastic material, small sliding). External loads are time-dependent but a quasi-static regime is assumed. The problem can be set as follows: find the displacement field and the Cauchy's stress field denoted respectively by u(x, t) and σ(x, t) satisfying:

• kinematic admissibility:

u(x, t) ∈ U [0,T ] and U [0,T ] = u ∈ H 1 (Ω) with u | ∂1Ω×[0, T ] = u p (1)
the strain-displacement (in small perturbations) relationship reads:

∀(x, t) ∈ Ω × [0, T ] : ε(u) = 1 2 ∇u + ∇ T u (2)
and the trace of the displacement field on the contact interface is

v = u | ∂3Ω×[0, T ] .
• static admissibility: σ(x, t) is balanced with external force and contact forces λ, i.e.

     ∀(x, t) ∈ ∂ 2 Ω × [0, T ] : σ n = f ext ∀(x, t) ∈ ∂ 3 Ω × [0, T ] : σ n = λ ∀(x, t) ∈ Ω × [0, T ] : div σ = 0 (3)
where n is the local outward normal vector.

• constitutive laws:

-Hooke law for elasticity:

∀(x, t) ∈ Ω × [0, T ] : σ = Kε(u) ( 4
)
where K is the linear elasticity operator.

-Coulomb frictional law for contact. With v and λ, the traces of displacement field and contact force field (both normal and tangential) on the contacting interface ∂ 3 Ω, these conditions can be written formally as in [START_REF] Dolbow | An extended finite element method for modeling crack growth with frictional contact[END_REF] as follows:

R(v, λ) = 0 (5)
The fields at the contacting interface ∂ 3 Ω can be split into a normal and a tangential part as:

∀(x, t) ∈ ∂ 3 Ω × [0, T ] : v(x, t) = u N n + u T λ(x, t) = λ N n + λ T (6)
According to Signorini conditions and Coulomb law on the displacement field and contact force field, operator R governs the contact behavior at the interface

∀(x, t) ∈ ∂ 3 Ω × [0, T ]:
• Normal contact with Signorini conditions given in [START_REF] Signorini | Questioni di elasticità non linearizzata e semilinearizzata[END_REF],

u N = u • n g Non-penetration condition λ N = (σ n) • n 0 Compressive contact force λ N (g -u N ) = 0
Complementarity condition [START_REF] Amsallem | An online method for interpolating linear parametric reduced-order models[END_REF] where g is the initial gap.

• Tangential contact with Coulomb friction law [START_REF] Coulomb | Théorie des machines simples[END_REF]:

λ T 2 = -µλ N if sliding: δu T 2 = 0 and δu T = -ρλ T with ρ > 0 λ T 2 -µλ N if sticking: δu T 2 = 0 (8)
where µ ∈ [0 , +∞[ is the friction coefficient and δu T is the increment of tangential displacement.

Weak formulation

We define the following set of homogeneous kinematically admissible functions as:

u(x, t) ∈ U [0,T ] 0 with U [0,T ] 0 = u ∈ H 1 (Ω) and u | ∂1Ω×[0, T ] = 0 (9)
Equilibrium ( 3) is equivalent to the following integral formulation with compatibility of contact force field and displacement field to contact conditions [START_REF] Amsallem | An online method for interpolating linear parametric reduced-order models[END_REF] and [START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF]:

∀u * ∈ U [0,T ] 0 : Ω σ : ε(u * )dV - ∂2Ω f ext • u * dS - ∂3Ω λ • u * dS = 0 (10)

Semi-discretized and discretized weak formulations

Using Hooke law for σ = Kε and the finite element approximation for displacement field yield to the so-called equilibrium equation with respect to frictional contact conditions at each contacting node:

∀u * (t) : T 0 u * T (t)Ku(t) dt = T 0 u * T (t) f ext (t) + f ctc (t) dt with R = 0 (11)
K is the stiffness matrix, f ext (t) and f ctc (t) (corresponding to λ(t)) are generalized forces and B is a boolean matrix mapping the global vector of nodal values to the values on contact boundary nodes. All in all, we have to find displacement field u(t) and contact force field f ctc (t) verifying [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF] for all weighting functions u * (t).

The time interval is discretized into a sequence of time steps t 1 k m describing a regular time stepping (i.e. t k+1 = t k + ∆t). Then, the time-space discretized reference problem is defined as follows:

∀t ∈ [t 1 , t m ] : Ku = f ext + f ctc with R = 0 ( 12
)
This system characterized by stiffness matrix K has a size n × n (number of degrees of freedom) and has to be verified for m time steps. The boolean mapping operator B is such that:

v = Bu = 0 1 0 0 0 1    ū u T u N    and f ctc = B T λ =   0 0 1 0 0 1   λ T λ N (13)
Relationship 12 defines a general non-linear problem including a non-smooth constitutive law which can be extended to other behaviors without much change in the formulation (cohesive contact, etc.).

Reduced basis and adaptation

Reduced basis methods are attractive techniques to tackle computational cost and data storage issues. They are related to numerous techniques (principal component analysis, POD...) and are also precious tools for model reduction techniques.

A posteriori approach and the LATIN/FAS principle

Given a reference problem, a posteriori strategies consist in using prior knowledge about the solution to compute inexpensively other close solutions (e.g. multiparametric approach). This knowledge is obtained from some solution snapshots (or coarse meshes, analytic solutions...) and are compressed into a small basis using for instance SVD or POD [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF]). Taking into account the most dominant directions, a RB is obtained. Then the equations of the problem can be Galerkin-projected into the RB providing the ROM for the considered problem. Efficiency of the ROM depends strongly on the relevancy of the RB. Thus, a compromise between precision and size of the ROM has to be reached.

In [START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF], analysis of RB provided by SVD are carried out for frictional contact problems and the LATIN/FAS strategy is proposed. It is exemplified that the RB has a multiscale content. First modes (the most energetic ones) are global at the overall structural scale, whereas higher order modes are localized to the contacting zone. The LATIN/FAS method proceeds as the non-linear FAS (Full Approximation Scheme) multigrid solver considering the RB as a coarse "grid". This strategy is embedded into the LATIN non-linear solver considered as a smoother. Coarse corrections are computed on the RB reducing error at the structural scale. Thus, the computational effort is concentrated to solve the localized frictional behavior. A consequent computational speedup is obtained and this strategy is efficient for multiparametric studies. Nevertheless, and in order to reach high levels of accuracy, the RB may not be able to produce relevant corrections. So, additional vectors are computed and added along the iterations to the RB (enrichment). Consequently, the main challenge of the LATIN/FAS is to provide a suited RB to compute relevant coarse corrections for iterated solution. Moreover, its suggested enrichment strategy may lead to expensive computations.

Hereinafter, the authors propose a robust PGD method consisting in solving the considered problem with an appropriate and evolving RB along iterations. Contrary to the LATIN/FAS method, the suggested a priori method does not require prior knowledge to design a suited RB which is computed on-the-fly and possibly from scratch. Furthermore, the strategy employed to achieve enrichment of the RB is far less expensive than the one proposed for the LATIN/FAS. In addition, it will be shown that the designed on-the-fly reduced basis has a quasi-optimal property. In other word, the computed RB is close to the best that one can compute for iterated solutions. Furthermore, as the PGD method is employed, a low rank approximation for the solution has to be handled. Consequently, a great computational work and huge memory requirements can be spared.

A priori LATIN-PGD method

Solution fields defined in [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF] for all time steps can be collected into a n-by-m snapshot matrix. For example, a snapshot matrix of the discretized displacement field is:

U =      u 1 (t 1 ) u 1 (t 2 ) • • • u 1 (t m ) u 2 (t 1 ) u 2 (t 2 ) • • • u 2 (t m ) . . . . . . . . . . . . u n (t 1 ) u n (t 2 ) • • • u n (t m )      = u(t 1 ) u(t 2 ) • • • u(t m ) (14)
The snapshot matrix ( 14) can be approximated using a separated form (or a low-rank approximation) as follows:

Ǔ = V 1 Φ T 1 + V 2 Φ T 2 + • • • + V p Φ T p u(t) = V 1 φ 1 (t) + V 2 φ 2 (t) + • • • + V p φ p (t) (15)
The vectors V k (respectively Φ k ) of size n (respectively m) are associated to spatial modes (respectively temporal modes φ k (t)). The number of pairs of space-time modes (or dyads) is denoted by p. Consequently, a basis can be computed (e.g. a spatial basis) to approximate Ǔ. Such a basis can be considered as a RB for U. In the following, we denote by V k p 1 the p-tuple of space vectors and by Φ k p 1 the p-tuple of time vectors. If p < 1, p-tuples is the empty set. The amount of data to store for ( 14) is of order n × m whereas for the low-rank form [START_REF] Boucard | A nonincremental approach for large displacement problems[END_REF] it is p(n + m). If p < mn m+n , then low-rank form is advantageous. Generally speaking, p has to be chosen such that the low-rank form [START_REF] Boucard | A nonincremental approach for large displacement problems[END_REF] generates an accurate approximation for ( 14) keeping an attractive size for the expansion [START_REF] Boucard | A nonincremental approach for large displacement problems[END_REF].

The PGD method is introduced in [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF]. It consists in building [START_REF] Boucard | A nonincremental approach for large displacement problems[END_REF] with progressive enrichment stages (i.e. given an approximation Ǔ, new vectors V and Φ are sought and added to the expansion) making Ǔ more accurate and close to U. This way and due to the fact that only vectors have to be computed and stored, computational work and memory usage is drastically reduced. Interested readers can also refer to [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF] for exhaustive reviews and details.

Optimality issues and downsizing

Intrinsically, the PGD method does not ensure the optimality of the separated form (15) (i.e. it is a non-unique form and the optimal one is such that p is minimal). According to [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF] and considering the Frobenius norm, the best separated form is provided by the SVD (the so-called optimality property).

To reach a certain level of accuracy for Ǔ, PGD process may need numerous enrichment stages leading to large values for p (with potentially memory overflow and computational cost issues). So for the sake of numerical efficiency, the size p of (15) has to be controlled during the PGD process.

To enforce optimality property, a brute force solution is to perform a full SVD computation on Ǔ after each enrichment step [START_REF] Golub | Singular value decomposition and least squares solutions[END_REF][START_REF] Baker | Low-rank incremental methods for computing dominant singular subspaces[END_REF]. A less expansive method consists in deploying SVD updating techniques [START_REF] Bunch | Updating the singular value decomposition[END_REF][START_REF] Brand | Fast low-rank modifications of the thin singular value decomposition[END_REF]. Nevertheless, all of those lies on burdensome methods to compute precisely the SVD of the approximation Ǔ. This high computational effort is not worthwhile knowing that Ǔ may be an inaccurate relatively to the solution of the reference problem.

This lead to assume that optimality property is useless and a soft downsizing approach suffices (i.e. maintain the quality of Ǔ while decreasing p). To design a suited algorithm, in a first step, an orthogonality property for V k p 1 is considered (i.e. V T i V j = δ ij with δ ij the Kronecker symbol). Then, the amplitude of each mode is attributed to the corresponding time vector Φ k p 1 . To preserve this orthonormality property at each enrichment stage, one can use algorithm 1.

Algorithm 1: Orthogonal enrichment.

Input:

Preexisting p-tuples:

V k p 1 , Φ p 1 p 1 with V T i V j = δ ij New pair of modes: V, Φ Output: Enriched (p + 1)-tuples: V k p+1 1 , Φ k p+1 1 with V T i V j = δ ij 1 for k = 1 to p do 2 Projection on existing space modes: α = V T k V 3 Update corresponding time modes: Φ k = Φ k + αΦ 4 Subtract projected component: V ← V -αV k 5 end 6 Norm of remaining space mode: β = V 2 7 Normalize space mode: V p+1 ← V/β 8 Amplify time mode: Φ p+1 ← βΦ
Notwithstanding orthonormality for V k p 1 , numerical experiments show that the inflation of the PGD basis is not stemmed, and a redundancy also occurs for the time modes Φ k p 1 . To compress time vectors, the authors propose the following dyad to dyad projection strategy. In order to illustrate it, we consider an application for an approximation Ǔ written with two dyads:

Ǔ = V 1 Φ T 1 + V 2 Φ T 2 with V i V j = δ ij (16)
Equation ( 16) can be rewritten considering that Φ 2 has a redundant component with Φ 1 (i.e. Φ 2 = αΦ 1 + Φ2 and ΦT 2 Φ 1 = 0) as follows:

Ǔ = (V 1 + αV 2 )Φ T 1 + V 2 ΦT 2 = V1 Φ T 1 + V 2 ΦT 2 (17)
Rewriting [START_REF] Brand | Incremental singular value decomposition of uncertain data with missing values[END_REF] does not respect orthogonality for space vectors (i.e. V 2 = β V1 + V2 and VT 1 V2 = 0). Applying a orthonormalization yields:

Ǔ = V1 (Φ T 1 + β ΦT 2 ) + V2 ΦT 2 = V1 ΦT 1 + V2 ΦT 2 = V1 V1 V1 ΦT 1 + V2 V2 V2 ΦT 2 (18)
This process does not change Ǔ and preserves orthogonality property for space vectors but does not ensure it for time vectors. The first projection step [START_REF] Brand | Incremental singular value decomposition of uncertain data with missing values[END_REF] allows to compress redundancy between Φ 1 and Φ 2 into a new vector Φ1 but looses the orthogonality for space vectors, therefore justifying step [START_REF] Brand | Fast low-rank modifications of the thin singular value decomposition[END_REF]. Incidentally, enforcing this property makes Φ1 and Φ2 no longer orthogonal. This process can be repeated until having Φ1 and Φ2 orthogonal (i.e. having produced the SVD of Ǔ).

Algorithm 2 generalizes this dyad to dyad process for larger expansions for Ǔ. The process [START_REF] Brand | Incremental singular value decomposition of uncertain data with missing values[END_REF], [START_REF] Brand | Fast low-rank modifications of the thin singular value decomposition[END_REF] is applied in a certain sorting of modes, selected with time vectors amplitudes (i.e. the dyad with the weakest amplitude is projected on the one with the largest amplitude).

Algorithm 2: Downsizing stage

Input: Existing p-tuples: V k p 1 , Φ k p 1 such that V T i V j = δ ij Max iterations: ξ max Cut-off amplitude: ε Output: Downsized q-tuples: V k q 1 , Φ k q 1 such that V T i V j = δ ij 1 for ξ = 1 to ξ max do 2 Index p-tuples such that: Φ 1 2 • • • Φ p 2 3 for i = p down to 2 do 4 for j = 1 up to i -1 do 5 Projection of time mode: α = Φ T j Φ i /(Φ T j Φ j ) 6
Update corresponding space mode:

V j ← V j + αV i 7 Subtract projected component: Φ i ← Φ i -αΦ j 8
Projection of space mode:

β = V T j V i /(V T j V j ) 9
Update corresponding time mode:

V i ← V i -βV j 10
Subtract projected component:

Φ j ← Φ j + βΦ i 11
Norm of the j th space mode:

ζ = V j 2 12
Normalization:

V j ← V j /ζ and Φ j ← ζΦ j 13 end 14
Norm of the i th space mode:

γ = V i 2 15
Normalization:

V i ← V i /γ
is c/ξ max = n( 5 2 p 2 -1 2 p -2) + m( 5 2 p 2 -3 2 p -1) + 3 2 p 2 -1 2 p -1.
A cut-off amplitude ε can be fixed relatively to a desired accuracy for Ǔ. Embedded in the PGD process, the downsizing algorithm 2 stems efficiently the inflation of the expansion Ǔ. In practice, a very few iterations ξ max (1 or 2 iterations) is sufficient. More iterations allow to make Ǔ converge to its SVD. To compare with classical and widespread techniques: Lanczos thin-SVD method has a O(nmp 2 ) complexity [START_REF] Golub | Matrix computations[END_REF] and fastest SVD updating methods have a complexity O(nmp) [START_REF] Brand | Incremental singular value decomposition of uncertain data with missing values[END_REF].

Formulation of the LATIN-PGD for frictional contact problems

One of the major ingredients of the proposed strategy is the LATIN method [START_REF] Ladevèze | Nonlinear Computational Structural Methods: New Approaches and Non-Incremental Methods of Calculation[END_REF]. This innovative non-incremental approach is well known for its ability to solve efficiently various non-linear time dependent problems such as frictional contact problems [START_REF] Champaney | Large scale applications on parallel computers of a mixed domain decomposition method[END_REF][START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF][START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF], large displacement [START_REF] Boucard | A nonincremental approach for large displacement problems[END_REF], non-linear materials [START_REF] Boisse | A new approach in non-linear mechanics: The large time increment method[END_REF][START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF][START_REF] Relun | A model reduction technique based on the pgd for elastic-viscoplastic computational analysis[END_REF], transient dynamics [START_REF] Ladevèze | A modular approach to 3-D impact computation with frictional contact[END_REF][START_REF] Lemoussu | A 3D shock computational strategy for real assembly and shock attenuator[END_REF][START_REF] Odièvre | A parallel, multiscale domain decomposition method for the transient dynamic analysis of assemblies with friction[END_REF]. The non-incremental LATIN method was initially proposed as a commitment of three principles:

(P1) Separation of the linear and non-linear behaviors. A denotes the set of solutions s = v, λ satisfying linear constitutive law, kinematic admissibility and static admissibility. These are defined on the whole space-time domain Ω × [0, T ]. Γ denotes the set of solutions ŝ = v, λ verifying frictional contact conditions and are defined locally at the contacting interface and on the whole time interval

∂ 3 Ω × [0, T ]. The solution of the problem is s ∈ A ∩ Γ.
(P2) A two-staged iterative algorithm. The solution of the problem is searched with the construction of two sequences of approximations belonging alternatively to A and Γ. At the i th iteration, the local stage consists in finding ŝi = vi , λi ∈ Γ with a search direction ŝi -

s i-1 = vi -v i-1 , λi - λ i-1 ∈ E + . Note that s i-1 = vi-1 , λi-1 is known from the previous iteration.
Then, the global stage consists in finding

s i = v i , λ i ∈ A with another search direction s i -ŝi = v i -vi , λ i -λi ∈ E -. Note that ŝi = vi , λi is known from the previous local stage.
(P3) Radial approximation or space-time separation. Unknowns are represented as a sum of products between a space function and a time function. This supplementary constraint makes the problem over-determined. In order to enforce the constraints in (P1), the search direction equation E - is then verified in a weak sense. It was initially stated to limit memory usage.

These principles are illustrated on figure 2. For certain cases and for sake of simplicity, the LATIN method can be formulated without the space-time separation stated in (P3). In this case several similarities can be stated with augmented Lagrangian methods [START_REF] Alart | Using Nonsmooth Analysis for Numerical Simulation of Contact Mechanics[END_REF].

Local stage Given a solution s = v, λ verifying the internal balance, kinematic admissibility and static admissibility, the local stage is an updating stage occurring at the contacting interface. At each node, non-linear frictional contact conditions and search direction equation have to be verified. The nodal search direction equations for both normal and tangential components reads:

(E + ) : ĉN = λN -λ N = k N vN -v N ĉT = λT -λ T = k T vT -v T (19)
k N and k T are parameters of the method defined at each node. It is similar to the penalty parameter in an augmented Lagrangian formulation. Different values for the normal and tangential problem can be chosen. These values influence only the convergence rate and optimal values are related to mechanical properties of the studied body [START_REF] Ladevèze | Nonlinear Computational Structural Methods: New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF][START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Trollé | Optimization of a stabilized X-FEM formulation for frictional cracks[END_REF]. According to [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF], close-to-optimum values are searched as follows: Normal components

k N = EL c and k T = EL c ( 20
) Γ×[0, T ] A×[0, T ] s s i = v i , λ i ŝi = vi , λi s i-1 = v i-1 , λ i-1 ŝi-1 = vi-1 , λi-1 s 0 = v 0 , λ 0 (a) Γ×[0, T ] A×[0, T ] s i = v i , λ i s i-1 = v i-1 , λ i-1 s s 0 = v 0 , λ 0 ŝi = vi , λi ŝi-1 = vi-1 , λi-1 ∆s i = (∆u i , ∆λ i ) (b)
λN = c N - vN = g- 1 k N c N + Tangential components Sliding: c T 2 > µ| λN | Sticking: c T 2 µ| λN | t = c T / c T 2 λT = µ| λN |t vk T = vk-1 T - 1 k T c T 2 -µ| λN | t λT = c T vk T = vk-1 T Table 1: Solutions of the local stage.
with E the Young's modulus and L c contact characteristic dimensions. k N,T have the dimension of a stiffness. Search direction equations can be equivalently rewritten introducing the gap parameter (for the normal component) and the displacement of the previous time step (for the tangential component, with a superscript k -1). This yields to the following equations:

c N = ĉN + k N g = λ N -k N (v N -g) = λN -k N (v N -g) c T = ĉT + k T vk-1 T = λ T -k T (v k T -vk-1 T ) = λT -k T (v k T -vk-1 T ) (21)
v k T refers to the nodal displacement on the contacting interface in the tangential plane at time step t k . Then, contact conditions can be easily applied on each term of c N and c T .

Explicit solutions of the local stage are given in table 1 and have to be computed for each node belonging to the contacting interface and for each time step. Hence, the solution of the local stage ŝ = v, λ verifies strictly frictional contact conditions and gives a prediction/correction given a solution s = v, λ verifying internal balance. Moreover, ŝ is computed explicitly.

Global stage Given a solution ŝ = v, λ issued from the previous local stage and verifying frictional contact conditions. The global stage consists in finding a space-time solution defined on the whole space domain and time interval verifying the linear constitutive law, kinematic admissibility, static admissibility. We introduce an additional search direction equation to introduce the local stage frictional contact prediction. This stage consists in finding a displacement field and a contact force field s = v, λ defined over ∂ 3 Ω × [0, T ]. These fields have to verify the admissibility equations:

Ku = f ext + B T λ v = Bu (22)
v is the trace of the displacement field over ∂ 3 Ω × [0, T ] whereas u is defined over the whole space-time domain. According to the second principle of the LATIN method, the search direction has to be verified as well. Let us write search direction equations for the whole contacting interface (with I d the identity matrix):

(E -) : λ N λ T -λN λT = k N I d 0 0 k T I d k vN vT - v N v T (23)
All in all, taking into account only the first and the second principle of the LATIN method leads to solve the following linear system for each time step:

     [K + B T kB]u = f ext + B T λ + kv v = Bu λ = λ + k v -v (24)
Note that the operator [K + B T kB] is symmetric definite positive and remains constant along the iterations. K is the stiffness matrix, taking into account the Dirichlet boundary conditions [START_REF] Alart | Méthode de newton généralisée en mécanique du contact[END_REF]. If these last conditions are not sufficient, K is only semi-definite positive. If the genuine problem is well-posed, the active and/or sticking contact conditions for the solution are sufficient to prevent any undefined rigid body motion. The full matrix [K + B T kB] is the previous stiffness matrix completed by a Robin-like boundary condition (with stiffness k) on the whole potential contact area. This last term is therefore sufficient to get a regular left hand side. Nevertheless, one has to take care of its condition number to solve efficiently such a linear system.

Initialization and stopping criterion Assuming that s ∈ A, so does the initial guess s 0 . Consequently, the iterative process is initialized with the linear contactless elastic solution:

s 0 = u 0 , λ 0 such that Ku 0 = f ext λ 0 = 0 (25)
We also use the convergence indicator [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF] introduced in [START_REF] Ribeaucourt | A new fatigue frictional contact crack propagation model with the coupled X-FEM/LATIN method[END_REF]. Iterations are stopped once a threshold is reached: I < ε.

I = max s N -ŝN 2 ∞ s N 2 ∞ + ŝN 2 ∞ , s T -ŝT 2 ∞ s T 2 ∞ + ŝT 2 ∞ with s N,T 2
∞ = max x,t 1 k N,T λ 2 N,T + k N,T u 2 N,T (26)
where k N,T is the search direction parameter. Note that subscripts N and T refer respectively to normal and tangential components at a given node belonging to the contacting interface. Using maximum norms, the convergence criterion ( 26) is a very tough criterion. A comparison with an other far less restrictive indicator is available in [START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF].

Global stage with space-time separation

The third principle consists in searching the solution of the reference problem (both displacement field and contact force field) under a separated form. Using the finite element framework, it reads:

u(t) = u 0 (t) + p k=1 V k φ k (t) λ(t) = λ 0 (t) + p k=1 L k ψ k (t) (27)
Vectors V k (resp. L k ) are space vectors associated to the displacement field defined over Ω (resp. contact force field defined over ∂ 3 Ω) and functions φ k (t) (resp. ψ k (t)) are time functions associated to V k (resp. L k ) and defined on [0, T ]. Vectors u 0 and λ 0 are initial guesses (which can be also provided as low rank approximations). The approximation ( 27) can be written equivalently using a matrix notation as follows:

u = u 0 + p k=1 V k Φ T k λ = λ 0 + p k=1 L k Ψ T k (28)
First, we reformulate equations of the global stage with a corrective scheme. At a given iteration (indexed by i), the correction ∆s = ∆u, ∆λ is computed and added to the solution of the previous global stage s i-1 , so that the quality of the updated solution s = s i-1 + ∆s improves. Subtracting appropriate quantities to equations [START_REF] Champaney | Large scale applications on parallel computers of a mixed domain decomposition method[END_REF] yields to the following problem to solve:

     K∆u = B T ∆λ ∆v = B∆u ∆λ + k∆v -res sd = 0 (29)
where

res sd = λ -λ i-1 + k(v -v i-1
) is known at this stage.

Second, we introduce the separated representation for corrective increments:

∆u = Vφ(t) on Ω × [0, T ] ∆λ = Lψ(t) on ∂ 3 Ω × [0, T ] (30)
Substituting representations [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF] in the internal balance equation of ( 29) yields to a relationship between space and time functions:

∀t ∈ [0, T ] : KVφ(t) = B T Lψ(t) ⇔ KV = B T L φ(t) = ψ(t) (31)
The problem (29) including representation [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF]

with W = L + kBV. The auxiliary space unknown W defined on ∂ 3 Ω is introduced to take into account the linear relationship (31) between V and L ensuring the internal balance. To compute them, the following relationships can be used:

[K + B T kB]V = B T W L = W -kBV (33)

Find time functions with fixed space functions

In this stage, one looks for solving global stage defined in [START_REF] Curnier | A generalized newton method for contact problems with friction[END_REF] with given space functions W k p 1 provided by an initial guess or by the previous iterations. Corrections φ k p 1 for time functions are sought according the following optimization problem:

φ k = argmin φ W k φ -res sd F (34)
A Galerkin method is applied with the trial function ω = W k φ (t) and yields to:

φ k = W T k res sd (35)
Consequently, updated time modes are

(φ k ) p 1 ← (φ k + φ k ) p 1 .
The corresponding pseudo-code is given in algorithm 3.

φ k ← φ k + W T k res sd 3 end
Because space functions are fixed, convergence may be slow and accuracy for the solution depends the given space basis. Usually, this stage is only deployed to update time functions and does not suffice to solve the considered problem. To make the space basis more relevant, extra space modes have to be computed and added.

Enrichment stage

The enrichment stage consists in generating a new pair of space-time functions extending the basis and making the iterated solution more accurate. So, the following optimization problem for the global stage defined on [START_REF] Curnier | A generalized newton method for contact problems with friction[END_REF] has to be solved:

W, φ(t) = argmin W φ(t) W φ(t) -res sd F ( 36
)
Based on a space-time weak form modeling method as [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF], a Galerkin method is applied with the trial function ω = W φ(t) + Wφ (t) and yields to:

T 0 ω Wφ(t) -res sd dt = 0 ( 37
)
The non-linear system ( 37) is solved with a robust iterative strategy [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF]. It proceeds in two steps, repeated until convergence. The first step consists in computing vector W knowing φ(t) from the previous step. Then, the second step consists in updating φ(t) knowing W from the first step. The process is stopped when W and φ(t) are no more significantly updated.

First step Space vector W is updated knowing φ(t) from the previous step:

W = T 0 res sd φ(t) dt T 0 φ 2 (t) dt ( 38
)
Second step Time function φ(t) is searched knowing W from the first step:

φ(t) = W T res sd (W T W) (39)
In practice, a very few iterations are needed and a normalization condition is applied for one of the two vectors. The pseudo-code of this method is given in the algorithm 4. Once the new pair is found, space and time basis are enriched using algorithm 1.

Preliminary stage

In [START_REF] Boisse | A new approach in non-linear mechanics: The large time increment method[END_REF], the preliminary stage consists in updating time functions as described in 4.1. Then, the enrichment stage is used only if the convergence criterion stagnates (e.g. if the error indicator decreases by a factor less than 10%). We propose herein to enhance the preliminary stage by updating space modes as well with corrections W k p 1 for fixed time modes. This new updating phase is defined from the global stage equation [START_REF] Curnier | A generalized newton method for contact problems with friction[END_REF] and the following optimization problem:

W k = argmin W Wφ k (t) -res sd F (40)
A Galerkin method with the trial function ω = W k φ k (t) is used and yields to:

W k = T 0 res sd φ k (t) dt T 0 φ 2 k (t) dt (41)
Because of orthonormality condition, space modes W k p 1 can not be updated directly and each correction W k has to be projected in the spatial basis. The corresponding pseudo-code is given in algorithm 5. Output: p-tuple of updated space modes: W k p 1

1 for i = 1 to p do 2 Space mode correction: W i = T 0 res sd φ i (t) dt T 0 φ 2 i (t) dt 3
Update: res sd ← res sd + W i φ i (t) Update corresponding time modes:

Φ j = Φ j + αΦ i 7 Subtract projected component: W i ← W i -αW j 8 end 9
Correction of space mode:

W i ← W i + W i 10
Norm of space mode:

β = W i 2 11
Normalization:

W i ← W i /β 12 Amplification: Φ i ← βΦ i 13 end
After the two updating stages for time and space modes, the enrichment stage is eventually used according to the decreasing of the convergence indicator. These stages fulfill the global stage and provide a kinematically and statically admissible solution. In addition, and once a global stage is solved, a downsizing stage defined with algorithm 2 is used in order to control the basis size.

The quasi-optimal LATIN-PGD algorithm for frictional contact

Two-dimensional numerical application

The proposed LATIN-PGD method is applied to the two dimensional frictional model illustrated in figure 3 with numerical properties given in Table 2. It is a flexible indenter subject to prescribed displacements. These are designed such that boundary ∂ 3 Ω is contacting a rigid part with successively different regimes (active contact or not, sliding, sticking) within the studied time interval.

To solve the initial reference problem, the LATIN non-linear solver is used.

A posteriori analysis of the reference solution

As in [START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF], a posteriori analyses can be carried out to exemplify the reducibility of the space-time reference solution (u r , λ r ) and its multiscale content. Moreover, a SVD basis of reference vectors W r k P 1 can be extracted from the following singular values decomposition:

W r 1 v r 1 (t) + • • • + W r K v r K (t) = svd (λ r -λ 0 ) + kB(u r -u 0 ) (42)
Associated vectors V r k P 1 and L r k P 1 correcting the elastic solution can be also computed and an approximated solution (ũ r , λr ) can be rebuilt taking into account only P modes. The larger is P , the more accurate is the solution (if P = K the solution is exactly recovered). On figure 4 different error indicators are provided to illustrate the reducibility of the solution. For instance, if one takes into account 10 modes, the error on the displacement field is 2.7 × 10 -4 , on the contact force field 4.8 × 10 -2 and the LATIN error measure is 3.3 × 10 -4 . All in all, depending on the field of interest, very few modes are sufficient to generate the solution. Usually, the displacement field needs less modes than the contact force field to be recovered.

A few space modes for both displacement and contact force field and associated time modes are given on figure 5. The same observations as in [START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF] can be made as far as the scale separability of a frictional contact problem is concerned: First modes (the most energetic ones) depict a global scale of the solution (scale of the structure and spread over the whole time interval) whereas higher order modes bring corrections localized at the contacting interface and at specific time steps.

Application of the quasi-optimal LATIN-PGD method

In this section, performances are compared between the LATIN method (without separated representation) and the LATIN-PGD (see figure 6). To reach the same level of precision, the LATIN-PGD requires less iterations than the LATIN method. Moreover, note that the complexity of one LATIN-PGD iteration is far less than one of the LATIN method. So the computation cost is reduced drastically using the LATIN-PGD.

A very interesting point is the evolution of the PGD basis. As the LATIN-PGD converges to the solution of the reference problem, the PGD basis has to span the same subspace as the one given by [START_REF] Fortin | Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems[END_REF]. Remarkably, the size of the PGD basis grows up until reaching, after a certain number of iterations, the size of the reference SVD basis.

The behavior of the PGD basis along iterations can be illustrated using Modal Assurance Criterion (MAC) diagrams [START_REF] Allemang | The modal assurance criterion-twenty years of use and abuse[END_REF]. Given two sets of vectors of same dimension (X i ) p 1 and (Y i ) q 1 , the MAC matrix M is defined as follows:

(M) 1 i p

1 j q = |X T i Y j | 2 X i 2 Y j 2 ∈ [0, 1] (43)
The coefficient M ij measures the correlation between modes X i and modes Y j .

If M ij = 1, then X i and Y j are colinear (highly correlated). On the contrary, M ij = 0 means that X i and Y j are orthogonal (uncorrelated). In practice MAC matrix is plotted using gray scale or 3D histogram to visualize correlation between a given set of modes and a reference one. Figure 7 reports the MAC matrices between the first 30 reference SVD modes and iterated PGD ones. The main observations can be pointed out:

• The PGD modes are clearly correlated to the corresponding SVD modes (according their order and depending on the precision of the solution).

• The PGD modes are computed hierarchically: structural modes are found at first.

• The larger the number of iterations is, the more accurate the solution is and the more correlations between PGD modes and SVD modes appear.

Thanks to the downsizing stage, iterated solutions are provided in a quasioptimal separated form (i.e. close to the SVD). To make it even more close to the reference SVD basis, the parameter ξ max can be set to a higher value. But and as it was explained in section 3.3, the strict optimality property is superfluous.

This observation emphasizes and exemplify a multiscale approach within the LATIN-PGD. The problem is firstly solved on the structural scales. Then, the problem is solved locally at the contacting interface. This strengthens the approach developed in [START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF] and fosters multiscale methods.

Initialization using coarse computations

In this section, a pre-computation method is proposed to initialize the PGD basis to accelerate convergence of the LATIN-PGD. It consists in finding global scale PGD modes (both in space and time) using a coarse space-time discretization (see figure 8). Indeed, the global scale behavior of the structure can be tracked using a coarse mesh and larger time steps. This pre-computation (less expensive than the reference model) allows to give a first guess of PGD modes. For the present example, we propose to use the coarse model described in the figure 9 and in the table 3.

Once the coarse problem is solved (i.e. the solution s C = u C , λ C ) is known), the displacement field is interpolated (using [START_REF] Dureisseix | Information transfer between incompatible finite element meshes: application to coupled thermo-viscoelasticity[END_REF] or a trilinear interpolation ...) on the reference space-time discretization providing ũ. Then, a first guess of the PGD basis can be deduced easily using the following formula: Reference space mode order PGD space mode order Table 3: Simulation parameters for the coarse model. Only first modes of (44) are considered (global modes) because higher order ones may be irrelevant (dependence on the discretization and the interpolation method). This initialization (see figure 10) provides interesting features. First, the initial solution is more accurate than the elastic solution (roughly 3 decades are gained). Second, the error of the initialized LATIN-PGD decreases rapidly during first iterations until recovering the same convergence rate as the uninitialized LATIN-PGD (pre-computed coarse modes are rapidly corrected on the fine discretization). In practice, to reach a very high-level of accuracy (says 10 -4 ), this initialization allows to divide by 2 the number of iterations (global modes are already found and few iterations are needed to capture the local behavior) relatively to the standard LATIN-PGD and by 4 relatively to the LATIN method. Nevertheless for higher levels of accuracy, more computational work is needed to describe precisely the local behavior diminishing relatively the gain of this initialization method. Of course, the gain of the initialization depends on the quality of the displacement projector from coarse to fine mesh.

W Φ T = svd BK(ũ -u 0 ) + kB(ũ -u 0 ) (44)
t Ω(t 0 ) Ω(t f ) Ω(t k ) W(t k ) (a) Reference space-time discretiza- tions t Ω(t C 0 ) Ω(t C f ) Ω(t C k ) W C (t k ) (b) Coarse space-time discretizations

Three-dimensional numerical application

The previous example can be easily transposed to a three dimensional case as described in the figure 11 and in the table 4.

Same a posteriori analyses can be carried out for the solution to exhibit the multiscale content of the SVD basis (see figures 12 and 13). As the previous two-dimensional problem, a scale separability phenomenon can be pointed out as far as spatial and time modes are concerned. The LATIN method, the LATIN-PGD and the LATIN-PGD initialized with 10 coarse space-time modes are compared on figure 14. LATIN-PGD overtakes again the LATIN method in terms of iterations. In addition, a LATIN-PGD iteration is less costly than one of the LATIN method. The initialization of the PGD basis allows to gain roughly one decade which is limited in comparison to the two dimensional previous example. To achieve it and for sake of simplicity, the authors use a collocation projection which is simpler but less accurate than [START_REF] Dureisseix | Information transfer between incompatible finite element meshes: application to coupled thermo-viscoelasticity[END_REF]. Nevertheless, the LATIN-PGD shows a great potential to solve 3D frictional problems well-known for their toughness.

Conclusion and perspectives

This article provides the formulation of the quasi-optimal LATIN-PGD solver dedicated to frictional contact problems. This new a priori reduced basis strategy converges faster compared to classic LATIN solver and requires less computational effort and memory usage. Moreover, the quasi-optimality property is a striking observation. The solution is found progressively according to its optimal modes provided by the SVD. This is achieved by a enhanced preliminary stage increasing the convergence rate and an efficient downsizing stage awarding the suggested approach a quasi-optimality property for PGD basis. At convergence, the solution of the considered problem is provided into a separated space-time form close to the SVD. Numerical applications illustrate the ability of the suggested approach to solve two-dimensional and three-dimensional problems. In addition and to take advantage of the scale separability of the space-time basis, an initialization method was proposed. Global and structural modes of the solution are computed on coarse space-time discretizations. Then, they are reused to provide an initial guess for the space-time basis on reference discretizations sparing a consequent computational effort.

In this paper, the tackled problem involves a single flexible body contacting a rigid basement for the sake of simplicity. Nonetheless, multibody system can be easily implemented by considering appropriate frictional contact conditions between a couple of flexible bodies. So, the local stage has to be generalized to these new conditions. But, the overall of the approach remains the same. The global stage aims still at solving internal balance with PGD space-time representation but different approaches are possible. The first consists in tackling globally all solids with a single family of space-time functions wherein whole unknowns are assembled. The second consists in tackling each solid individually with a dedicated family of space-time functions. This last approach is suited for parallel computing but harder to implement.

Furthermore, the LATIN method provides also a framework to take into account other kind of non-linearities such as finite strains, non-linear material such as plasticity accounting the space-time separation. Then, it would be possible to take into account those supplementary non-linearities in the presented work with the developed method.

This strategy can be easily extended to the multiparametric framework as introduced in [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF]. For that purpose, a parameter of the reference problem (e.g. frictional coefficient) may considered as a variable of the problem. Then, the approximated displacement field solution is still sought as a separated form solution u(x, t; µ) = u(x, t, µ) = u 0 (x, t, µ) + i V(x)φ(t)f (µ). Nevertheless, downsizing approaches and compression techniques such as PARAFAC [START_REF] Bro | Parafac. tutorial and applications[END_REF][START_REF] Bro | A new efficient method for determining the number of components in parafac models[END_REF] or promising HOSVD [START_REF] Luo | Are tensor decomposition solutions unique? on the global convergence hosvd and parafac algorithms[END_REF][START_REF] Kolda | Tensor decompositions and applications[END_REF][START_REF] Lathauwer | A multilinear singular value decomposition[END_REF][START_REF] Ma | Dynamic updating and downdating matrix svd and tensor hosvd for adaptive indexing and retrieval of motion trajectories[END_REF] remain challenging.

Figure 2 :

 2 Figure 2: Illustration of LATIN non-linear iterative solver. (a) using (P1) and (P2) principles -(b) using (P1), (P2) and (P3) principles: dashed arrows symbolize exact search directions whereas waving line are approximated search directions.

Algorithm 3 : 1 p 1 1 1

 3111 Time modes update. Input: res sd p-tuple of existing space modes: W k p -tuple of existing time modes: φ k p Output: p-tuple of updated time modes: φ k p for k = 1 to p do 2 Update time modes:

Algorithm 4 :1 5 Space mode amplitude: α = W 2 6 7

 4527 Enrichment stage. Input: res sd Max subiterations: ξ max Output: New space vector: W New time mode: φ(t) Initialize φ(t) (e.g. φ(t) = 1) 2 for ξ = 1 to ξ max do 3 Compute space mode: W = T 0 res sd φ(t) dt T 0 φ 2 (t) dt 4 Compute time mode: φ(t) = W T res sd (W T W) Normalize space mode: W = W/α Amplify time mode: φ(t) = αφ(t) 8 end

Algorithm 5 : 1 p

 51 Space modes update Input: res sd p-tuple of existing space modes: W k p -tuple of existing time modes: φ k p 1

4 for j = 1 to p do 5

 45 Projection of the correction: α = W T j W i 6

Algorithm 6 : 1 1 1 = 4 ** Local stage ** 5 7 * 8 if p 1 then 9 Compute 10 Time modes update → φ k p 1 11 1 13 end 14 * 17 Enrichment stage → {W, φ} 18 Orthogonal enrichment → { W k p+1 1 , φ k p+1 1 } 19 Increment basis size p ← p + 1 20 end 21 ** Downsizing stage ** 22 if p 1 then 23 Downsizing1 , φ k q 1 } 24 Update basis size p = q 25 end 26 ** Convergence check ** 27 Save

 61145789101114171811191212223124252627 Quasi-optimal LATIN algorithm for frictional contact Data: Precisions ε,εmax Enriching criterion 0 θ Compute linear elastic solution Ku0 = fext and λ0 = 0 2 Initialization of the PGD basis p = 0, W k p ∅ and φ k p 1 = ∅ 3 while I > ε do Build interface solution → (v, λ) 6 Solve local stage → (v, λ) * Extended preliminary stage ** res sd = ( λ + kû) -(λ0 + kBu0) + p k=1 W k φ k (t) Update res sd = ( λ + kû) -(λ0 + kBu0) + p k=1 W k φ k (t) 12 Space modes update → W k p * Enrichment stage ** 15 if (1 -I/I old < θ) then 16 Compute res sd = ( λ + kû) -(λ0 + kBu0) + p k=1 W k φ k (t) stage with ξmax = 1 and ε = εmax → { W k q previous convergence criterion I old ← I 28 Compute convergence criterion I 29 end

Figure 3 :Figure 4 :

 34 Figure 3: (a) Model -(b) Prescribed displacement evolution -(c) Mesh -(d) Shape of the deformed body at t 150 .

Figure 5 :

 5 Figure 5: Space modes V 1,3,7 (displacement field), L 1,3,7 (contact force field) computed from W 1,3,7 and time modes φ 1,3,7 (t).

  basis Evolution of the size of PGD basis versus iterations LATIN-PGD (b)

Figure 6 :

 6 Figure 6: (a) Error given by the formula (26) versus iterations -(b) Size of the PGD basis versus iterations.

Figure 7 :

 7 Figure 7: MAC diagrams after a certain number of iterations for both space modes and time modes.

Figure 8 :

 8 Figure 8: A cheap pre-computation on coarse discretizations can be used to catch global time-pace modes to initialize LATIN-PGD on the fine reference discretizations.

Figure 9 :

 9 Figure 9: Coarse mesh used for pre-computation.

Figure 10 :Figure 11 :

 1011 Figure 10: Convergence plots with an initial guess for the PGD basis.

Figure 12 :

 12 Figure 12: Errors of rebuilt solution taking into account a given number of modes for the 3D problem.

Figure 13 :Figure 14 :

 1314 Figure 13: Space-time modes V 1,3,5 (displacement field), L 1,3,5 (contact force field) computed from W 1,3,5 and φ 1,3,5 (t) for the 3D problem.

Table 4 :

 4 Simulation parameters.

			Fine model Coarse model
	Number of time steps	m	200	50
	Number of dofs	n	12099	7920
	Number of contacting nodes	n C	385	172
	Young modulus	E	10 6 Pa	10 6 Pa
	Poisson ratio	ν	0.3	0.3
	Friction coefficient	µ	0.15	0.15
	Final time	T	1 s	1 s
	Number of reduced SVD modes		160	44
	Search direction	k T	3000 N/m	3000 N/m
	Search direction	k N	1000 N/m	1000 N/m

The whole pseudo-code of the quasi-optimal LATIN-PGD is given in algorithm 6. Note that, the cut-off amplitude criterion of the downsizing stage is choosen in relation to the LATIN-PGD convergence indicator.

Acknowledgment

We gratefully acknowledge the French National Association for Research and Technology (ANRT, CIFRE grant number 2011/0705).


H A Gueniche

P A Glaude
email: glaude@ensic.inpl-nancy.fr

R Fournet

F Battin-Leclerc

Experimental study of the structure of a rich premixed 1,3-butadiene/CH 4 /O 2 /Ar flame

The structure of a laminar rich premixed 1,3-C 4 H 6 /CH 4 /O 2 /Ar flame have been investigated. 1,3-Butadiene, methane, oxygen and argon mole fractions are 0.033; 0.2073; 0.3315, and 0.4280, respectively, for an equivalent ratio of 1.80. The flame has been stabilized on a burner at a pressure of 6.7 kPa (50 Torr). The concentration profiles of stable species were measured by gas chromatography after sampling with a quartz probe. Quantified species included carbon monoxide and dioxide, methane, oxygen, hydrogen, ethane, ethylene, acetylene, propyne, allene, propene, cyclopropane, 1,3-butadiene, butenes, 1-butyne, vinylacetylene, diacetylene, C 5 compounds, benzene, and toluene. The temperature was measured thanks to a thermocouple in PtRh (6%)-PtRh (30%) settled inside the enclosure and ranged from 900 K close to the burner up to 2100 K.

Introduction

Soots and polyaromatic hydrocarbons (PAH) present in the exhaust gas of engines and combustion devices represent a high part of urban pollution. Their formation in combustion involves small unsaturated hydrocarbons for which the chemistry is still very uncertain. Different ways are proposed for the formation and the oxidation of the first aromatic compounds, involving the reactions of C 2 (acetylene), C 3 or C 4 unsaturated species [1][2][3][4][5].

Allene and propyne flames have been previously investigated [START_REF] Gueniche | Proceedings of theEuropean Combustion Meeting[END_REF] to study the determinant role of the propargyl radicals in benzene formation. Among the other unsaturated soot precursors implied in the aromatic ring formation, 1,3-butadiene is of primary importance as a primary product in hydrocarbon oxidation and because of its structure involving two conjugated unsaturated bonds. Furthermore, 1,3-butadiene is an interesting molecule in studying the formation of benzene, due to the possibility of branching between C 2 , C 3 , and C 4 chains. The oxidation of butadiene has been already studied in some oxidation conditions: shock tube [START_REF] Fournet | [END_REF], flow reactor [8], and premixed flat flame [1]. The previous work in flame involved a mixture oxygen/1,3-butadiene/argon 67.5/29.5/3.0 corresponding to an equivalence ratio of 2.4 under a pressure of 2.67 kPa (20 Torr), in near-sooting conditions.

The purpose of the present work is to experimentally investigate the structure of a premixed flat 1,3-C 4 H 6 /CH 4 /O 2 /Ar flame. These results will be use in a further extend to improve the mechanism previously developed in our laboratory for the reactions of C 3 -C 4 unsaturated hydrocarbons (propyne, allene, 1-butyne, 2-butyne, 1,3-butadiene) [START_REF] Fournet | [END_REF]9].

Experimental procedure

The experiments were performed using an apparatus newly developed in our laboratory to study temperature and stable species profiles in a laminar premixed flat flame at low pressure. The body of the flat flame matrix burner, provided by McKenna Products, was made of stainless steel, with an outer diameter of 120 mm and a height of 60 mm (without gas/water connectors). The burner, used in this investigation, was built with a bronze disk (95% copper, 5% tin). The porous plate (60 mm diameter) for flame stabilization was water cooled (water temperature: 60°C) with a cooling coil sintered into the plate. The burner could be operated with an annular co-flow of argon to favor the stabilization of the flame.

This burner was housed in a water-cooled vacuum chamber evacuated by two primary pumps and maintained at 50 Torr (6.7 kPa) by a regulation valve. This chamber was equipped of four quartz windows for an optical access, a probe for samples taking and a thermocouple for temperature measurements. The burner could be vertically translated, while the housing and its equipments were kept fixed. A sighting telescope measured the position of the burner relative to the probe or the thermocouple with an accuracy of 0.01 mm. The flame was lighted on using an electrical discharge.

Gas flow rates were regulated by RDM 280 Alphagaz and Bronkhorst (El-Flow) mass flow regulators. 1,3-Butadiene had a 99.5 % purity, methane a 99.95 % purity and both were supplied by Alphagaz -L'Air Liquide. Oxygen (99.5% pure) and argon (99.995% pure) were supplied by Messer.

Temperature profiles were obtained using a PtRh (6%)-PtRh (30%) type B thermocouple (diameter 100 µm). The thermocouple wire is sustained horizontally by a fork and crosses the flame to avoid conduction. The junction is put at the centre of the burner. In the case of the flame perturbed by the probe, the junction of the thermocouple is located at two probe's inlet diameter or so from the probe's tip, corresponding to about 200 µm. The thermocouple was coated with an inert layer of BeO-Y 2 O 3 to prevent catalytic effects [10]. The ceramic was done by damping in a hot solution of Y 2 (CO 3 ) 3 (93% mass.) and BeO (7% mass.) followed by a drying in a Meker burner flame. This process was reiterated until the whole metal was covered. Radiative heat losses were corrected using the electric compensation method [11]. The temperature measured by the thermocouple had been first calibrated as a function of the heating electrical power in a high vacuum. The real temperature of the gas was obtained thereafter at the intersection of the calibration curve and the experimental one measured in the flame.

Stable species profiles were determined by gas chromatography. Gas samples were directly obtained by connecting the quartz micro probe to a volume, which was previously evacuated by a turbo molecular pump down to 10 -7 kPa and which was then filled up to a pressure of 10 Torr (1.3 kPa). The diameter of the extremity of the probe was about 100 µm, while the upper part diameter was 6 mm, with a angle of ~40°. The sudden decrease of the pressure from the flame gases to the inside of the probe insured all reactions to be frozen.

Chromatographs with a Carbosphere packed column and helium or argon as carrier gas were used to analyse O 2 , H compounds, benzene, and toluene) were analysed on a haysep packed column by flame ionisation detection (FID) and nitrogen as gas carrier gas. Water was not directly measured, but its mole fraction was determined by mass balance on O atom. The change in the mole fractions due to the dilution by the reaction products was taken into account from a first simulation. Some small oxygenated compounds were detected but were not quantitatively analysed.

Results and Discussion

A methane flame doped with 1,3-butadiene have been investigated. The gas flow rate was 3.32 l/min and the mixtures contained 20.7% methane, 3.3% 1,3-butadiene, 33.2% oxygen and 42.8% argon corresponding to an equivalent ratio of 1.80. C 4 compounds represented then 16% of the amount of methane. Figure 1 presents the temperature profiles obtained for the flame. The lower temperatures were measured at 0.04 cm from the burner and were around 600 K without the probe, 440 K with the probe. The higher temperatures were reached between 0.7 and 0.9 cm above the burner and were around 2150 K without the probe. The temperature decreased thereafter because of the heat losses.

Figure 2 to 11 present the profiles of the reactants and of the products of the reaction in function of the height above the burner. Figure 2 displays the profiles of the fuels, methane and 1,3butadiene, and figure 3 those of oxygen, hydrogen and water that was calculated from the mass balance. Carbon oxides are presented on figure 4. Figure 5 shows the profiles of C 2 products, figures 6 and 7 those of C 3 products, figures 8 and 9 those of C 4 products. Some C 5 products, such as methyl butadiene or cyclopentadiene, have been measured. They have not been displayed here since they have not been calibrated yet. Figure 10 displays the profile of benzene and figure 11 that of toluene that is the heaviest analyzed compound.

The consumption of both fuels occurs in the first stage of the flame, but 1,3-butadiene is consumed a little bit closer to the burner, at 5 mm height while some methane remains up to 7 mm. The higher reactivity of the butadiene can be explained by the weaker C-H bonds on the tertiary Catoms of this reactant (99.9 kcal/mol [9]) that allows easier H-abstractions yielding a resonance stabilized i-C 4 H 5 (CH 2 =CH-(•)C=CH 2 ) radical:

C 4 H 6 + R• → RH + i-C 4 H 5 Furthermore the two double bonds involved in this fuel allow the addition of small radicals, especially H atoms yielding a resonance stabilized butenyl free radical.

C 4 H 6 + H → C 4 H 7 The shape of the profile of oxygen is very similar to that of methane (figure 3). Final carboncontaining products are CO and CO 2 (figure 4). Since the mixtures are rather rich, ϕ = 1.80, CO remains the major carbon oxide in the post-flame. Hydrogen is also one of the major products with a mole fraction reaching 0.2 in the end gas. This species is produced by the H-abstractions by Hatom on the fuels. The low amount of oxygenated free radicals caused by the high equivalence ratio does not permit to oxidize thereafter H 2 to water.

The main stable intermediates are C 2 molecules, which are mainly yielded by the combustion of methane (figure 5). Ethane and ethylene are produced promptly and reach their higher amount close to the burner. In this flame, the peak of ethylene is located a little bit closer to the burner than that of ethane, while that of acetylene is located further at 0.6 cm. These molecules are produced in the combustion of methane via the combination of methyl radicals that is yielded by a H-abstraction from the reactant. CH 3 + CH 3 → C 2 H 6 The adduct can also decompose directly CH 3 + CH 3 → C 2 H 5 + H Ethane reacts by H-abstraction and yields also C 2 H 5 radicals, that decompose quickly into ethylene C 2 H 5 → C 2 H 4 + H The prompt formation of ethylene, which peak is reached before that of ethane, can be due to reaction of 1,3-butadiene, for instance by addition with H-atoms followed by the decomposition of the adduct:

C 4 H 6 + H → C 2 H 4 + C 2 H 3 A minor amount of ethylene react by addition with small radicals, but the main channel is the Habstraction by H or OH yielding C 2 H 3 . These radicals oxidize then with O 2 or decompose into acetylene.

C 2 H 3 → C 2 H 2 + H Since this last molecule is the most stable in flame conditions, its amount remains higher than that of other C 2 species. In the rich mixtures studied here, the amount of acetylene decreases slowly in the post flame.

Among the C 3 products, the most important is propene. It can come from both oxidation of methane and butadiene. C 3 H 6 can be produced by the termination of methyl radicals with C 2 H 3 or its addition on C 2 H 4 followed by the decomposition of the adduct CH 3 + C 2 H 4 → C 3 H 6 +H C 3 H 6 can also be produced by addition of methyl radicals on butadiene followed by the decomposition of the adduct:

CH 3 + C 4 H 6 → C 3 H 6 + C 2 H 3
In this case, it is mostly resulting of a coupling between the two fuels since methyl radical are yielded by a H-abstraction on methane.

Cyclopropane must be mainly formed by the molecular isomerization of propene.

Propyne and allene are also present among the stable intermediates in the flame (figure 7). They can be produced by many pathways involving C 2 , C 3 and C 4 compounds. A first route yielding propyne is the addition of methyl radicals on acetylene followed by the decomposition of the adduct:

C 2 H 2 + CH 3 → p-C 3 H 4 +H The H-abstractions on propene produced allyl radicals that decompose also to allene or propyne. The H-abstractions on 1,3-butadiene yield resonance stabilized i-C 4 H 5 that can decompose to allene and a methyl radical. The amounts of allene and propyne are linked together by the molecular isomerization that is easier in the high temperature conditions.

The major C 4 products are 1-butyne, 2-butene (figure 8), vinylacetylene and diacetylene (figure 9), which are produced in the very first stage of the combustion. Even if some routes from the C 2 species are involved in their formation, they must be yielded mainly in the combustion of 1,3butadiene. 2-Butene can be produced by the termination reaction between resonance stabilized i-C 4 H 5 and H-atom while 1-butyne can be produced by successive molecular isomerizations. 1,3butadiene isomerized to 1,2-butadiene that produced then 1-butyne. 1,2-Butadiene has not been quantified experimentally since it could not be separated from the reactant 1,3-butadiene by chromatography.

Vinylacetylene is produced directly from 1,3-butadiene by a H-abstraction yielding i-C 4 H 5 followed by a decomposition by β-scission: The heaviest and most unsaturated compounds quantified here are benzene (figure 10) and toluene (figure 11). In these conditions, benzene can be produced by both the C 3 and C 4 -pathways, as proposed by many authors [2][3][4]. A first route is the combination of propargyl radicals that are produced by allene and propyne: C 6 H 5 + CH 3 → C 7 H 8 Some addition reaction C 5 and C 2 compounds followed by internal rearrangements could also be involved in the formation of toluene.

C 4 H 6 + R• → RH + i-
C 3 H 3 + C 3 H 3 → C 6 H 6 C 3 H 3 + C 3 H 3 → C 6 H

Conclusion

This paper presents new experimental results for a rich premixed laminar flat flame of methane and 1,3-butadiene. Profiles of temperature and of stable species have been measured for many products and in peculiar many unsaturated compounds that are soot precursors. The first aromatic rings, benzene and toluene have been measured and in these conditions, their formation seems due to both C 3 and C 4 routes.

Future work will include a kinetic modeling of these data to improve the reaction mechanisms of small soot precursors responsible for the first aromatic rings formation. Flames doped with other unsaturated compounds such as cyclopentene will be also investigated.

Figure 1 :

 1 Figure 1: Experimental temperature profiles for the 1,3-butadiene flames with and without probe perturbation, measured with a coated type B thermocouple.

Figure 2 :

 2 Figure 2: Mole fraction of the fuels.

Figure 3 :

 3 Figure 3: Mole fraction of oxygen, hydrogen and water

Figure 4 :

 4 Figure 4: Mole fraction of carbon oxides.

Figure 5 :

 5 Figure 5: Mole fraction of ethane, ethylene and acetylene.

Figure 6 :

 6 Figure 6: Mole fraction of propane and cyclopropane.

Figure 7 :

 7 Figure 7: Mole fraction of allene (propadiene) and propyne.

Figure 8 :

 8 Figure 8: Mole fraction of 1-butyne and 2-butene.

Figure 9 :

 9 Figure 9: Mole fraction of vinylacetylene and diacetylene.

Figure 10 :

 10 Figure 10: Mole fraction of benzene.

Figure 11 :

 11 Figure 11: Mole fraction of toluene.

  2 , CO and CO 2 by thermal conductivity detection and CH 4 , C 2 H 2 , C 2 H 4 , C 2 H 6 by FID. Heavier hydrocarbons (a-C 3 H 4 , p-C 3 H 4 , c-C 3 H 6 , C 3 H 6 , C 4 H 2 , C 4 H 4 , 1,3-C 4 H 6 , 1-C 4 H 6, 2-C 4 H 8 , C 5

  5 + H The adduct can stabilized or decompose into phenyl radical and H-atom. Another channel involves the C 4 compounds deriving from 1,3-butadiene. C 4 H 5 radical can add on acetylene yielding benzene plus H: C 4 H 5 + C 2 H 2 → C 6 H 6 + H C 2 H 3 can add too on C 4 H 4 , the adduct react then by internal addition to form a benzene ring: C 2 H 3 + C 4 H 4 → C 6 H 6 + H Toluene can be produced by the addition of methyl on benzene followed by the elimination of a H atom: C 6 H 6 + CH 3 → C 7 H 8 + H It can be produced too by the termination step between phenyl and methyl radicals:

Acknowledgements

Financial support of this work by the European Union within the "SAFEKINEX" Project EVG1-CT-2002-00072 is gratefully acknowledged.