Extended-PGD Model Reduction for Nonlinear Solid Mechanics Problems Involving Many Parameters

Reduced models and especially those based on Proper Generalized Decomposition (PGD) are decision-making tools which are about to revolutionize many domains. Unfortunately, their calculation remains problematic for problems involving many parameters, for which one can invoke the “curse of dimensionality”. The paper starts with the state-of-the-art for nonlinear problems involving stochastic parameters. Then, an answer to the challenge of many parameters is given in solid mechanics with the so-called “parameter-multiscale PGD”, which is based on the Saint-Venant principle.


Introduction
Numerical simulation has made a forceful entry into design and analysis offices. This revolution, which is anything but complete, has entered in a new stage, called simulation-driven "robust" design. It leads to a major scientific challenge: simulations should be performed in quasi real-time. The key is a new generation of reduced-order methods which comprises essentially the Proper Orthogonal Decomposition (POD), the Proper Generalized Decomposition (PGD) and the Reduced Basis Method (RB), the basis and recent developments of which are given in [6]. These are problems in which uncertainties or variations in parameters are to be taken into account, or problems with very high number of degrees of freedom, with multiple scales or interactions between several physics. These methods, together with the notions of offline and online calculations, also open the way to new approaches where simulation and analysis of structures can be carried out almost in real time.
The object of this work is the PGD, that was introduced in [13,14] for the treatment of nonlinear time-dependent problems of solid mechanics. Many developments have been made over the last thirty years: multi-scale, multiphysics, stochastic or non-stochastic parameters, acoustics, large displacements and deformations ...
In [18], one can find a synthesis of most of the developments in Cachan, where the method LATIN plays a central role. Currently, a number of tools have become mature and have been applied to industrial cases and then are competitors of classical computational methods (see book [7]). This is the case for the calculation of (visco)-plastic structures with less than ten parameters. The PGD not only makes it possible to construct reduced models that can be used in real time, but it also makes it possible to reduce drastically the calculation time in many situations. However, a major limitation is still the number of the parameters that can be involved (no more than 10). The paper starts with the state-of-the-art for nonlinear problems involving parameters which could be stochastic. It appears that the key to extend the PGD to a large number of parameters and then to cancel the "curse of dimensionality" is the same for linear and nonlinear problems.
Several attemps have been introduced for problems with a large number of parameters: Tucker tensors, Tensor Train tensors, Hierarchical Tucker tensors or more general tree-based Hierarchical Tucker tensors [3,8,9,12,24]. As they are generic a priori approximations, we believe that their improvement which is real, is limited. Here, in this paper we are following quite a different way. The idea is to notice that the operator underlying the problem treated is a binder between the parameters. Thus, at the level of the solution, the "curse of dimensionality" is no longer relevant. More, we go further with the parameter-multiscale PGD which takes its source in a fine analysis of the solution according to parameters, analysis built on the Principle of Saint-Venant. The basic idea is to introduce a two-level description of each parameter, the "macro" scale and the "micro" scale as one considers two scales for the space or the time. Also, to implement this vision, we have been led to use completely discontinuous approximations [18]. To carry out these idea, we use the Weak-Trefftz Discontinuous Method introduced in [16] and applied in [20] for the calculation of "medium frequency" phenomena. This paper is limited to fundamental aspects and the first numerical experiments for linear problems. More details can be founded in [25].

Nonlinear Problems with Parameters:
The State-of-the-art We have been working on ROM computation for 30 years with the so-called LATIN-PGD and what we are doing at the present time is the result of many works. PGD means Proper Generalized Decomposition and LATIN denotes the computational method which is nonincremental. The LATIN-PGD method was introduced in [13] for viscoplastic materials whose constitutive relations are described using a functional approach. Its extension to modern material descriptions involving internal variables, still for viscoplastic materials, was proposed in [14]. A number of mathematical properties regarding convergence and error indicators were proved in the book [15]. Overview could be found in [15,18]. Originally, PGD was called radial loading approximation, which, to us, meant a "mechanics" approximation in solid mechanics. LATIN-PGD leads today to a general and robust PGD computation technique which is based on an "abstract" reformulation of parametric nonlinear solid mechanics problems defined over a time-space-parameter domain. This work is the follow-up to [11,19,22,27]. Today, there are few other works except the works of Ryckelynck and its group done with POD [28]. However, things are changing and today there are more and more POD-approaches developed in relation with the homogenization technique FE2 [10,21,26]. Additional reduction or interpolation are then introduced to reduce their computation cost. This is done offline with the PGD computation. Our answer is the so-called Reference Point Method (RPM) [4].

Notations
With the assumption of small perturbations, let us consider the quasi-static and isothermal evolution of a structure defined over the time-space domain [0, T] × and depending of the parameters ∈ . This structure is subjected to prescribed body forces f d , traction forces F d over a part 2 of the boundary, and displacements u d over the complementary part 1 (see Fig. 1). All the data and the material characteristics depend on the parameters ∈ . The state of the structure is defined by the set of the fields = (̇p,Ẋ, , Y) (where the dot notation□ denotes the time derivative), in which: • p designates the inelastic part of the strain field which corresponds to the displacement field u, uncoupled into an elastic part e such that p = − e ; X designates the remaining internal variables; • designates the Cauchy stress field and Y the set of variables conjugate of X (Y and X have the same dimension). X could be hardening variables, damage variables, chemical variables, . . .
• generalized quantities over the parameter-time-space domain.
All these quantities are defined over the parameter-time-space domain × [0, T] × and assumed to be sufficiently regular. For the sake of simplicity, the displacement u alone is assumed to have a nonzero initial value, denoted u 0 . Introducing the following notations for the primal fields: and for the dual fields: The mechanical dissipation rate for the entire structure is: where ⋅ denotes the contraction adapted to the tensorial nature of X and Y. Notation • denotes the contraction operator for generalized quantities. Let us introduce the following fundamental bilinear "dissipation" form: along with and , the spaces of the fieldṡp and which are compatible with (4). These spaces enable us to define = × , the space in which the state = (̇p, ) of the structure is being sought.

The State Equations over × [ , T] ×
Following [15], a normal formulation with internal state variables is used to represent the behavior of the material. If denotes the mass density of the material, from the free energy Ψ( e , X) with the usual uncoupling assumptions, the state law yields: where the Hooke's tensor and the constant, symmetric and positive definite tensor are material characteristics.

× [ , T] ×
The state evolution laws can be written: where is a positive operator which is also for most viscoplastic models maximal monotone [15]. Let us introduce now the space ,0 the associated vectorial space. The compatibility equation can be written as: (7) It follows that the stress = ( (u) − p ) can be written: where is a linear given operator and d is a prestress depending on the data. Introducing the generalized stress, the admissibility conditions can be written as: where is a linear symmetric positive operator. Finally, the problem to solve, which is defined over Consequently, one has to solve a first order differential equation with an initial condition; The operators and as well as the right-hand-side member d depends on the parameter belonging to the parameter set .

The LATIN Solver
The LATIN method is an iterative strategy which differs from classical incremental or step-by-step techniques in that, at each iteration, it produces an approximation of the complete structural response over the whole loading history being considered. A review of the state-of-the-art and more recent extensions could be found in [18]. The LATIN method is designed as a mechanics-based computational strategy whose aim is to achieve the best possible performance level for solid mechanics problems. Consequently, this alternative approach is rooted in some remarkable properties which are verified by most of the models encountered in structural mechanics.
The LATIN method operates here over the parameter-time-space domain × [0, T] × , and its first principle (P1) consists in separating the difficulties. Thus, the equations are divided into: • a set of linear equations which can be global in the space variables: the equilibrium and compatibility equations, the state equations; • a set of equations which are local in the space variables but can be nonlinear: the state evolution laws.
The reformulation (11) of the reference problem enters into this framework because: • is a linear operator; • is at least local in space variables.
In the geometric representation given Fig. 2, and represent the solutions of the first and second set respectively, being associated to the free energy and with the dissipation. The LATIN second principle (P2) is also very natural. It consists in solving the two sets of equations alternatively until practical convergence. In order to do that, one uses search directions given as parameters of the LATIN method. One possible choice (Newton search direction) consists of the tangent direction and its conjugate direction (see Fig. 3).
Local stage principle at iteration n + -Find̂n +1∕2 = (̂̇p ,n+1∕2 ,̂n +1∕2 ) ∈ such that: Fig. 2 The geometric representation associated to the reformulation of the reference problem The search direction + is a LATIN parameter. Practically, one takes a linear positive operator which is local both on time and space variables. This local stage is very suitable for parallel computing.
Linear stage principle at iteration n + -Find n+1 = (̇p ,n+1 , n+1 ) ∈ such that: The search direction − is a LATIN parameter. This is a linear positive operator which is local both on time and space variables. It is associated to the material operator . One has to solve a first order linear differential equation with an initial condition, the operator being non-explicit.
Remark In practice, − is chosen close to the tangent to the manifold at the point n+1∕2 = (̂̇p ,n+1∕2 ,̂n +1∕2 ). For + , one takes or − . The convergence of the iterative process has been proved in the case of non-softening materials and contacts without friction [15]. The distance between two successive approximations gives a good and easily computed error indicator. Let us also note that one often uses an additional relaxation with a coefficient equal to 0.8.

PGD Computation
Let us recall that the problem is defined over the parameter-time-space domain Linear stage at iteration n + -Let us introduce corrections: 7 where n+1 = (̇p ,n , n ) has been computed at iteration n. The problem to solve over × [0, T] × at iteration n + 1 is then: The main idea is to interpret the search direction as a linear constitutive relation, the operator − being local both on time and space variables and positive definite as the Hooke tensor. Consequently, one introduces the associated constitutive relation error which will be minimized: and with = (̇p, ) ∈ . The problem (15) becomes: Find ∈ minimizing ∈ ↦ R( ) ∈ ℝ with the constrains =̇p and p = 0 at t = 0 One only prescribes that: with i (0) = 0 (initial condition), i ∈ L 2 ( ), i (t) ∈ L 2 [0, T] and i ( ) ∈ L 2 ( ). It follows, using admissibility conditions, that: where i (x) are computed solving several elasticity problems. For the sake of simplicity, let us consider now that does not depend on t and belonging to × [0, T] × . The general case does not involve serious difficulties, the -constraint being satisfied in a mean sense.
A greedy algorithm with updating is used to solve the minimization problem. More details can be founded in [17,18]. For few parameters, it could be advantageous to describe point-by-point the parameter space using the remarkable property of the LATIN method: the initialization of the iterative process can be any function defined over [0, T] × [27].
Local stage at iteration n + -One has the solve a problem which is local over the parameter-time-space domain for which "hyper-reduction" techniques are also welcome. A very popular technique is the Empirical Interpolation Method (EIM) [2] and its discrete version named DEIM [5]. The Hyperreduction method [28] makes the most of a restricted subdomain of the space domain. Here we use the Reference Point Method (RPM) developed in [4]. Let us note that such method can be also used for solving the linear stage.

An Engineering Illustration
To illustrate the use of the technique to deal with parametrized problems, we consider an example issued from [22] and which is freely inspired from a blade of the Vulcain engine of the Ariane 5 launcher. The geometry, boundary conditions and mesh are presented on Fig. 4. A four-sinusoidal-cycles displacement with is prescribed on the lower part. The total number of DOFs is 141,500 and the time interval is discretized using 120 time steps. The material coefficients used for the Marquis-Chaboche elastic-viscoplastic material are typical of a Titanium TA6V material at 500 • K. The parametric study is defined by Table 1. It concerns concerns the influence of the loading amplitude, of the limit stress and the power coefficient in the evolution law. The range of variation of each parameter was discretized into 10 values, leading to 1000 different problems. Figure 5 shows the total agreement of the results obtained with the LATIN-PGD compared to the one obtained with ABAQUS.  To illustrate the performances of LATIN-PGD, we compare compare the CPU times: • about 50 days (estimated time) are necessary to complete the 1000 resolution with ABAQUS; • less than 17 h are necessary to complete the 1000 resolution with the multiple runs algorithm.
The gain is a about 70 using the muliple runs algorithm, but can achieve more than 700 when using also the RPM strategy. The important point is that, once this parametric study has been performed, a reduced-model of the nonlinear problem is built. For example, some stochastic studies can be performed very easily. Assuming a probabilistic distribution of the three parameters, a probabilistic distribution of the quantity of interest can be computed in quasi-real time. Let us consider now that the first parameter is fuzzy and the two others follow a normal law. Then, it is easy to compute the probability law related to the maximum Mises stress over [0, T] × (see Fig. 6).

Model Problem
The basic principles of parameters-multiscale-PGD are given in [17]. One considers an elastic media which occupies the domain divided into sub-domains or elements E , E ∈ . The parameters E ∈ ℝ q E are associated to the rigidity of the volume E . For the sake of simplicity, we consider that E is a scalar; it belongs to [−1∕2; 1∕2]. Let us introduce ≡ { E ⧵ E ∈ }, the corresponding space being . The problem to solve can be written as: where is a linear positive definite operator depending on . F d is a given loading which could also depend on .
Numerical illustration with the standard PGD-The standard parameter-PGD has been introduced by [1]. Stochastic framework has been considered by [23]. An overview is given in [6]. Figure 7 shows an illustration of the model problem. This is a cube submitted to a uniaxial traction displacement, the opposite face being clamped. The parameters could be interpreted as damage intensity. Convergence curves for the classical PGD are given Fig. 8 for 8, 27 and 64 parameters; they show that convergence cannot be obtained when the number of parameters is more than 30.

The Key Idea
One has to solve: ( )X( ) = F d Let us introduce 0 = (0) and let us suppose for the sake of simplicity that F d does not depend on . One has: Let be X 0 = 0 −1 F d ; the solution can be written as: where X 1 is -linear and X 2 -quadratic. It follows that with only few terms, one gets a good approximation and then, the curse of dimensionality disappears; the operator links the parameters. Let us go further: on each element, one has: where AE∈ is the finite element assembly operator. Introducing the operator E , giving the restriction V E of a spatial vector V on the subdomain E through the It follows: The term E Z E is associated to a self-equilibrated loading. Therefore, from the Principle of Saint-Venant, the solution is localized in the neighborhood of the element E, essentially over the elements sharing a common point with E noted C E . Let be: which can be seen as negligible over the complement of C E . One has: and X 1 ( ) = − ∑ E∈ E Z 1,E ; X 1 is then -linear and Z 1,E is localized over the neighborhood C E of E. A similar property holds for X 2 : It follows from the Saint-Venant principle that the quadratic term of X 2 ( ) is such that Z 2,EE ′ is in practice localized over a close neighborhood of C E . The second term is linear with respect to each parameter; more, the E -contribution concerns essentially the neighborhood of C E .
The parameter-multiscale PGD that we propose is based on these remarks. Thus we introduce two scales, micro and macro, to describe the parameter space . In this way, we propose as approximation: where f M and g m are respectively "macro" and "micro" functions. C E denotes a chosen neighborhood of E defining the "micro" impact in space of the parameter E . One can take the elements having a common point with E. Here, we will choose a linear discretization of the "macro" functions f M E ′′ and thus consider only their value on two points, E ′′ = {±1∕2}.
However, there is a difficulty: X( ) is discontinuous from an element to another. It is removed using the so-called Weak-Trefftz Discontinuous Galerkin method (WTDG) proposed in [19] and extended to the quasistatic loadings in [18]. This approach is introduced in the next paragraph.

The WTDG Method
The domain is still split into elements or subdomains E , E ∈ on which the Hooke tensor E is constant but depends on E , belonging to [−1∕2, 1∕2]. The classical WTDG is modified here by adding a regularization term which assure the positivity of the operator, even if the rigidity vanishes over one or several subdomains.
The admissible space associated to E ∈ is denoted: U h, E,ad and the associated vectorial space: U h, E,0 . The problem to solve is then: = 0 (23) with = K (U). EE ′ is the common boundary of two adjacent subdomain E and E ′ ; EE is the common boundary of E and . One consider here approximations such that the interior equation is satisfied in an average sense, i.e. in resultant and moment over E . Two elements WP1 and WP2 can be easily associated to the classical elements P1 and P2. With these elements, one gets a coercivity property leading to the unicity of the solution [17].
On each element E, the WTDG leads to a contribution: ∑ The symmetric part is E and its value for = 0 is 0 E . The generalized force given by the WTDG is then:

Basic Operators
• Computation of E (R E ) ≡Z E ∈ ℝ n R is a known residual, the contribution to the subdomains E ∈ being R E . One computes here a search direction in space using as conditioner the symmetric operator 0 E . Let be ( E ( ),Z E ) ∈ ℝ × ℝ n ; the corresponding space is × ℝ n . One defines: The -integration can be done using the macro description.
• Computation of the extensionX = A E∈ a EZ E with a = (R,Z) FromZ E , E ∈ , ones defines a space search direction over . One has: which is equivalent to find the maximum of the Rayleigh quotient: which is relatively easy to compute.
• Computation of the "macro" functions f M = (R,X) One minimize the residual: and then where The classical PGD technique is used with for example two complete iterations over the parameter space.
• Computation of the "micro" functions {g EE ′ , E ′ ∈ C E } = E (R,X, f M ) Let us consider that the initial residual is: One introduces a new approximation which is equal to f MX over the complementary part of C E : This residual minimization problem which is here a small problem is solved classically.

Computational Strategy
The computational strategy uses error indicators as: The optimal computational strategy is under study. More details can be found in [17,25].

A First Illustration
The parameter-multiscale PGD has been implemented on a monodimensional example. A displacement is imposed at the extremity of a cantilever beam and each element has an independent Young modulus. A solution for a random set of parameters can be seen Fig. 9. This one-dimension example respects the Saint-Venant principle in stress, but not in displacement, requiring the addition of a rigid-body displacement associated to each micro function. This specificity will disappear in a two or three dimensional space.
First, we have computed the errors associated to the different approximations following the development (21) in the case where the variations of E are ±50%. Let us introduce the following norms and errors which are E-independent for the studied problem: Fig. 9 Displacement in a beam, each element has an independent and random Young Modulus One gets: • 0-order approximation: = 0.5 = 0.29 • 1-order approximation: = 0.25 = 0.10 • 2-order approximation: = 0.11 = 0.042 Then, the parameter-multiscale PGD has been implemented with the following strategy on a problem involving 64 parameters. For each element, one starts the computation with a micro-function corrected by a macro-one. The error at this stage is already small: Let us note that such an approximation gives the same error for any number of parameters. After introducing macro-functions, the error continues to decrease but slowly.

Conclusion
The parameter-multiscale PGD seems to be a very promising way to built reduced order model for problems involving a large number of parameters. Further work will be devoted to the derivation of verification tools and to the extension to nonlinear problems such as viscoplastic ones.