Estimation of the validity domain of hyper-reduction approximations in generalized standard elastoviscoplasticity

Background: We propose an a posteriori estimator of the error of hyper-reduced predictions for elastoviscoplastic problems. For a given fixed mesh, this error estimator aims to forecast the validity domain in the parameter space, of hyper-reduction approximations. This error estimator evaluates if the simulation outputs generated by the hyper-reduced model represent a convenient approximation of the outputs that the finite element simulation would have predicted. We do not account for the approximation error related to the finite element approximation upon which the hyper-reduction approximation is introduced. Methods: We restrict our attention to generalized standard materials. Upon use of incremental variational principles, we propose an error in constitutive relation. This error is split into three terms including a tailored norm of the hyper-reduction approximation error. This error norm is defined by using the convexity of an incremental potential introduced to state the constitutive equations. The second term of the a posteriori error is related to the stress recovery technique that generates stresses fulfilling the finite element equilibrium equations. The last term is a coupling term between the hyperreduction approximation error at each time step and the errors committed before this time step. Unfortunately, this last term prevents error certification. In this paper, we restrict our attention to outputs extracted by a Lipschitz function of the displacements. Results: In the proposed numerical examples, we show very good preliminary results in predicting the validity domain of hyper-reduction approximations. The average computational time of the predictions obtained by hyper reduction, is accelerated by a factor of 6 compared to that of finite element simulations. This speed-up incorporates the computational time devoted to the error estimation. Conclusions: The numerical implementation of the proposed error estimator is straightforward. It does not require the computation of the incremental potential. In the numerical results, the estimated validity domain of hyper-reduced approximations is inside the reference validity domain. This paper is a first attempt for a posteriori error estimation of hyper-reduction approximations.


Background
Reduced-Order models aim at reducing the computational time required to obtain solutions of Partial Differential Equations (PDE) which are physically-based and parameter dependent. They reduce the computational complexity of optimization procedures or parametric analyses [1][2][3].
Reduced-Order Models are very useful to model complex mechanical phenomena, when parametric studies are mandatory to setup convenient nonlinear constitutive equations or boundary conditions. In such situations, many open questions about physical assumptions arise and the scientist's intuition alone cannot guarantee that convenient simplified numerical approximations are chosen. Here, we consider hyper-reduced (HR) simulations, involving both a reduced-basis approximation and reduced-coordinate estimation by using a mesh restricted to Reduced Integration Domain (RID) [4,5]. The introduction of the RID is crucial for elastoviscoplastic models, because in many practical cases, no speed-up is achieved if the mesh is not restricted to the RID. In such a framework, error estimation provides a valuable algorithmic approach to check if hyperreduced simulations have been performed in the validity domain of the hyper-reduced model, and if they allow a physical understanding of the simulation outputs.
HR simulations may not be performed in faster time than predictions using metamodeling, response surface methods, or virtual charts [6]. But, unlike these methods, Reduced-Order Models involve equations related to physical principles (i.e. the balance of momentum), physical constitutive equations and data, validated by experiments. Such an approach enables certification of outputs of interest [7], or error estimation of predictions obtained by using Reduced-Order Models [7][8][9][10][11]. Moreover, error estimation is mandatory in many model reduction of parametric Partial Differential Equations, when one needs to evaluate, in the parameter space, a trust region or a validity domain related to Reduced-Basis accuracy [12][13][14]. Error estimation is also mandatory to perform a priori model reduction by using greedy algorithms [15,16] or the A Priori Hyper-reduction method [4].
Unfortunately, in nonlinear mechanics of materials, parametric PDEs are rarely affine in parameter, and it is quite difficult to reduce equations and data to an approximate affine form as proposed in [15] by using the Empirical Interpolation Method [17,18]. The novelty of this paper is the use of an incremental variational principle [19][20][21] to develop an a posteriori estimator of errors for HR predictions in case of generalized standard elastoviscoplastic models. As usual, the finite element (FE) solution fulfills incremental equations in time, but no differential equations related to the continuous time variable. The incremental variational principle enables to preserve formulations discrete in time when comparing the HR solution to the FE solution. This approach enables to state a relationship between a tailored norm of the hyper-reduction approximation error and the a posteriori estimator of this error. Here, we term hyper-reduction approximation error (HRAE) the difference between the FE solution of the PDE and its estimation by the hyper-reduced model. The reference solution is the FE solution, not the exact solution of the PDE. In the sequel, the proposed tailored norm of the displacement is denoted by u u , where u is the predicted displacement over the domain . The error estimator is denoted by η. We follow the theory of constitutive relation error (CRE) proposed by P. Ladevèze [22]. This kind of a posteriori error has been developed for a large set of constitutive equations [9,[23][24][25][26]. Although constitutive relation error (CRE) can handle all approximation errors due to time discretization, FE approximation, and eventually a separated representation of the displacements [27], we restrict our attention to approximation errors due to the additional assumptions incorporated in the numerical model by the hyper-reduced formulation. Here, the CRE η requires the construction of stresses fulfilling FE equilibrium equations.
In the numerical part of this paper, we apply the error estimation to forecast the validity domain of hyper-reduced predictions of simulation outputs. Adjoint sensitivity analyses have been introduced for certification of outputs in [8,28]. In case of time dependent problems, the adjoint equation or the mirror problem [25] is solved backward in time from T to 0. Here, we circumvent such a complex numerical solution by restricting our attention to Lipschitz functions of the displacements for simulation outputs, similarly to [29]. The output function is denoted by s(u) ∈ R N s . The Lipschitz continuity of the output s(u) reads: where V h is the convenient Hilbert space related to the FE approximation of the PDE, · F is the Frobenius norm on R N s . The elastoviscoplastic problem is both discretized in space and in time by the finite element method and the forward Euler scheme. The solution of this FE problem at discrete time t n ∈ {t o , . . . , t N t } is denoted by u n FE : where (ξ i ) N ξ i=1 are the FE shape functions, x the space coordinate, μ is a vector of parameters and D is a given parameter space.The hyper-reduction approximation of this FE solution is denoted by u n HR : where (ψ k ) N ψ k=1 is a given reduced-basis, (γ n k ) N ψ k=1 are the reduced coordinates evaluated by using equilibrium equations set-up on the RID denoted by Z , Z ⊂ , by following the hyper-reduction method [5]. The HRAE is denoted by e n , e n = u n FE − u n HR . Because of the time discretization, we cannot consider the approximation error e n separately for each time index n. We must account for the error at all discrete times simultaneously. Therefore, the upper bound of interest related to the Lipschitz continuity reads: The reference validity domain of the hyper-reduced model is denoted by D V . Its definition reads: where D is a given tolerance. The proposed numerical error estimation aims to forecast a reliable approximation of D V by the recourse to hyper-reduced simulations, a pos-teriori estimations of errors and a FE simulation used to generate reduced-bases. The approximate validity domain is denoted by D V : where c η is a substitute for the Lipschitz constant c o , η is a substitute for the tailored norm of the HRAE, (u n HR ) N t n=1 is the sequence of hyper-reduced predictions over the full time interval and ( σ n ) N t n=1 is a sequence of stresses fulfilling the FE equilibrium equation at each discrete time t n , where σ n is an affine function of the stress estimated from u n HR . In the sequel, both sequences (u n HR ) N t n=1 and ( σ n ) N t n=1 are denoted by u HR and σ respectively. Various Reduced-Basis construction methods can be incorporated into the hyperreduced model. For the sake of simplicity, in the following numerical examples, the Reduced-Basis is provided by a Proper Orthogonal Decomposition [30][31][32] of the FE solution related to a given parameter value μ 1 ∈ D. The POD modes are orthonormal solutions of the following mimization problem: where ·, · is the L 2 scalar product defined on and · L 2 ( ) is the related norm. The hyper-reduced model is set-up in order to provide an accurate approximate solution for μ = μ 1 : c η η u HR ·; μ 1 , σ ·; μ 1 HR , where HR is a given tolerance. Hence μ 1 ∈ D V . Moreover c η fulfills the following constraint: N t n=1 s u n FE ·; μ 1 − s u n HR ·; μ 1 2 F c η η u HR ·; μ 1 , σ ·; μ 1 (9) In this paper μ 1 is arbitrary chosen. Nevertheless, the knowledge of an approximate validity domain aims to validate this choice or to adapt it. In the sequel, for clarity, we do not mention the vector of parameters μ in the notations.

Methods
This section is organized as follows. The generalized standard constitutive equations are presented first. They are based on an incremental variational principle proposed in [33]. Then, we introduce the hyper-reduced incremental problem to be solved. The following section is devoted to the construction of the FE-equilibrated stress within the framework of the hyper-reduced modeling. The next sections introduce the error estimator, the tailored norm related to the CRE and the partition of the CRE. We finish by remarks on the numerical implementation of the proposed approach.

Incremental potential for the constitutive equations
The constitutive laws are described by using an incremental potential in the framework of the irreversible thermodynamic processes. A priori error estimators and incremental variational formulations were introduced in [33] for mechanical problems of bodies undergoing large dynamic deformations. Extensions of this approach were proposed in [20,21] for effective response predictions of heterogeneous materials. The strain history is taken into account by using internal variables denoted by α. These variables are the lump sum of the history of material changes. This approach has its roots in the works by Biot [34], Ziegler [35], Germain [36] or Halpen and Nguyen [37] and has proven its ability to cover a broad spectrum of models in viscoelasticity, viscoplasticity, plasticity and also continuum damage mechanics. The FE solution is approximated by an hyperreduced predictions denoted by u n HR , α n HR , σ n HR respectively for the displacements, the internal variables and the Cauchy stress, at discrete time t n . For the sake of simplicity, we denote by ε n HR the strain tensor ε(u n HR ), in the framework of the infinitesimal strain theory.
According to the framework of the incremental variational formulations, the constitutive law is defined by a condensed incremental potential, denoted by w (ε n HR ), such that: In practice, we do not implement the computation of the partial derivative of w with respect to the component of the strain tensor. Nevertheless, Equation (10) is fulfilled up to the computer precision. The internal variables are solutions of the following equation [21]: where w(ε, α) is the free energy of the material, and ϕ(α) is its dissipation potential [36]. The two potentials w and ϕ are convex functions of their arguments (ε, α) andα respectively, according to the theory of generalized standard materials [36,37]. Examples of constitutive laws can be found in [38]. A detailed example is given in the last section of this paper. The convexity of w is proved in [21] under the assumption that w and ϕ are convex functions. In the sequel, the explicit knowledge of the condensed incremental potential is not required for the mathematical formulation of the error estimator.
By recourse to the Legendre transformation, the dual potential of w , denoted by w , reads: We restrict our attention to convex functions w , hence: Therefore: The definition of the partial derivatives ∂w ∂ε and ∂w ∂σ is extended to fulfill the following equations: Hence, by construction, σ n HR and ε n HR fulfill the following equation:

Hyper-reduced setting
The continuous medium is occupying a domain . The boundary ∂ of is denoted by We denote by V the matrix form of empirical modes defined on the FE basis: By following the hyper-reduction method [39], we generate the RID Z by using the empirical modes (ψ k ) N ψ k=1 and the strain modes (ε(ψ k )) N ψ k=1 . When the RID is known, we determine the set of available FE residual entries when the stress estimation is restricted to Z . This set is denoted by Z such that: Therefore, we can define truncated test functions insuring a weak form of the balance of momentum restricted to the domain Z , since these test functions are set to zero over \ Z . These truncated test functions are denoted by (ψ Z k ) N ψ k=1 such that: The statement of the hyper-reduced incremental problem is the following. Given a parameter vector μ, find an estimation of the reduced coordinates γ n such that u n HR ∈ u c + V ROM fulfills the constitutive equations and the principle of virtual work: where : is the contracted product for second-order tensors. We must emphasize the fact that Equation (24) gives access to the displacement in all the domain although the equi-librium is set only on Z . Hence, when the hyper-reduced solution is known, the stress σ n HR can be estimated on the full domain by using Equation (10). The FE equations are obtained by substituting in equations (24) to (26), the following items: However, in Equation (27), the incremental potential related to the hyper-reduced prediction is not the incremental potential of the FE prediction, because α n−1 FE and α n−1 HR may differ. So, σ n FE fulfills the constitutive equation, but it may not be a derivative of w . Therefore, we introduce a correction term in stress, denoted by δσ n , such that: where ε n FE = ε(u n FE ) and δσ n account for the variation of the internal variables due to the difference between ε i FE n−1 i=1 and ε i HR n−1 i=1 at time instants before t n . The convexity of w and φ insures that the FE solution is unique, if no rigid mode is available.

Dual reduced-subspace
In this section, we introduce the technique of construction of the stress fields σ n fulfilling the finite element equilibrium equation at each discrete time t n . These stress fields are used in the next sections to build a CRE. The FE equilibrium equation is a linear equation upon stresses denoted by σ n . It defines an affine space such that σ n ∈ σ n N + S h , where σ n N is a particular solution of Neumann boundary conditions (i.e. in the linear elastic case for instance), and S h is the following vector space: The following linear elastic problem gives us σ n N : where E o ∈ R + is an abritary constant. By the recourse to the stress predicted by the FE simulation for μ = μ 1 , we obtain the snapshots of stress fields σ n FE (·; μ 1 ) − σ n N N t n=1 that span a subspace of S h : Similarly to the approach proposed in [10] for elasticity, we introduce a dual reducedsubspace denoted by S ROM ⊂ S h . S ROM is generated by the usual POD applied to σ n FE (·; μ 1 ) − σ n N N t n=1 . Its dimension is denoted by N σ ψ , and it is such that N σ ψ ≤ N t .
Hence, the hyper-reduced predictions of the stress tensor can be projected into the space of admissible stresses fulfilling the FE equilibrium equation. This projection reads: The above minimization problem is a global L 2 projection of the stress σ n HR onto the reduced basis that span S ROM . The dual basis being a POD basis, it is orthonormal with respect to the L 2 scalar product. So the computational complexity of the stresses projection scales linearly with the number of Gauss points of the mesh and N σ ψ . In practice, this computational complexity is negligible compared to the complexity of the evaluation of σ n HR , by using the constitutive equations.

The constitutive relation error and its partition
The reference solution being incremental in time, no continuous formulation in time is considered for the error estimator. The error estimation proposed in [33] for incremental variational formulation, is related to an asymptotic convergence assumption in order to get an upper bound of the approximation error related to the FE discretization. This upper bound depends on a constant related of the weak form of the partial differential equations.
Here, we propose to apply the Constitutive Relation Error method proposed in [22] to estimate the constant c η without assuming an asymptotic convergence of the HRAE. Following the method proposed in [22], a constitutive relation error, denoted by η, can be introduced as follows: where σ n is defined by Equation (34). The following properties hold: The first property (36) comes from the Legendre transformation (14). The property (37) comes from Equation (19). The proof of Property (38) is the following. If η(u n HR , σ n ) = 0 ∀ n ∈ {1, . . . , N t } then (u n HR , σ n , α n HR ) fulfills the constitutive equations, the Dirichlet conditions and the FE equilibrium equation. As (u n HR , σ n , α n HR ) fulfills all the equations of the FE problem, it is a solution of the FE problem and e n = 0.
We propose for e n u an Hilbert norm parametrized by the approximate solution u n HR and the exact solution u n FE . This parametrized norm is defined by: where ε(·) is the symmetric part of the gradient and G a symmetric positive-definite fourth-order tensor. If the identity tensor is substituted for G, we obtain a usual H 1 ( ) norm. We assume that there is no rigid mode, neither in the FE solution nor in the reduced basis.
The following mathematical developments aim to establish a relation between η(u HR , ∂w ∂ε (ε FE )) and the tailored norm e n u . We choose G by the recourse to the convexity of the incremental potential w . According to the Legendre transformation, we have: Hence, we can remove the contribution of the dual potential w * in η u HR , ∂w ∂ε (ε FE ) : Let us define the scalar function f (λ): We assume that f is C 2 on [ 0, 1]. The application of the Taylor Lagrange formula at order 2 leads to: The function evaluation and its derivatives read: Therefore, by choosing: we obtain the following property: This property is an intermediate result before establishing the relationship between e n 2 u and the error estimator. The incremental potential being strictly convex, G is definite positive. A schematic view of N t n=1 e n 2 u is shown on Figure 1. Similarly, a norm on stresses can be defined to measure the distance between the admissible stress σ n and σ n FE − δσ n , by the recourse to the dual potential w such that: where H is a positive symmetric fourth order tensor. To establish the relation between H and the Hessian matrix of the potential w , we consider η(u FE , σ ). Intermediate stresses (S n (λ)) N t n=1 are introduced such that: By the recourse to the Taylor Lagrange formula, we obtain the following property: Property. The proposed constitutive relation error fulfills the following partition: The last term of the sum is a coupling term between the error e n and the HRAE committed before the discrete time t n . Since we can't certify that this term is positive, η is not an upper bound of the HRAE.
Thanks to the definition of δσ by Equation (28), we have the following properties: We start proving the partition property by considering the following sum of terms: In A, the contributions of w and w are canceled when considering the definition of η, and we obtain: −ε n HR : σ n + ε n HR : σ n FE − δσ n + ε n FE : σ n − ε n FE : σ n FE − δσ n d Then: where ε n FE − ε n HR = ε(e n ). Therefore, according to properties (53), (55) and (56), the following partition is fulfilled: ε(e n ) : δσ n d Since e n ∈ V h and σ n FE − σ n ∈ S h , then: This ends the proof.
Property. If e n = 0, for n ∈ {1, . . . , N t } then u n HR = u n FE , α n HR = α n FE , δσ n = 0 and: In the general case, the closer σ n to σ n FE − δσ n the better the error estimation. The proposed error estimator incorporates errors related to the projection of the stress onto the dual reduced basis.

Numerical implementation of the CRE
The numerical implementation of the incremental potential and its dual is not required for the error estimation. Let's consider the following function of a scalar coordinate λ: Then: Therefore, the following property holds: Here, ε n (λ) is the strain tensor given by the constitutive equation upon the stress τ (λ) at time t n and the internal variables α n−1 HR at time t n−1 . Hence, the numerical estimation of ε n (λ) can be performed by the usual implementation of the constitutive equations of generalized standard materials [40][41][42].
In the following numerical simulations, we have estimated the integral on λ by the value of first derivative of g for λ = 1.

Numerical estimation of the constant c η
In practice, we have more modes than necessary to achieve accurate predictions for μ = μ 1 , because the POD of (u n FE (·, μ 1 )) N t n=1 is accurate up to the numerical precision. The total number of available modes is denoted by N ψ such that N ψ ≤ min(N ξ , N t ) and N ψ ≤ N ψ . This gives access to several hyper-reduced models by restricting the number of modes involved in the reduced basis related to the displacements. The constant c η is setup by using these approximate hyper-reduced solutions and the known stresses predicted by the FE model: min η u HR ·; μ 1 , σ FE ·; μ 1 , η u HR ·; μ 1 , σ ·; μ 1 (66) Hence the constraint (9) is fulfilled. Let us denote by N c ψ the number of displacement modes for which the maximum in Equation (66) is reached. In our opinion, if we expect that the error estimator behaves like an upper bound, we should not take N c ψ modes to generate the final HR model. In the following example, we are setting N ψ = N c ψ + 1.

Results and discussion
We apply the estimation of the validity domain to a beam composed with a layer of tin alloy (SnX) on a layer of copper (Cu). This is a cantilever beam shown in Figure 2. The thickness of the tin is 0.155 mm, the thickness of the copper in 0.111 mm, the width and the length of the beam are respectively 2. mm and 12. mm. These are the dimensions of an experimental specimen used at the Centre des Matériaux. The simulation output is the difference between an experimental transverse displacement at point A, shown in Figure 2, and its prediction by the mechanical model: where y is the transverse axis and u n exp is a given experimental measurement of the transverse displacement of the point A. The copper is assumed to be elastic and isotropic. Its Young modulus is 63000. MPa and its Poisson coefficient is 0.3. The tin alloy is isotropic and elastoviscoplastic with a linear kinematic hardening. The potentials of the tin alloy read: The FE simulation used to generate the reduced basis is performed on an arbitrary sampling point μ 1 = (20. MPa, 30000. MPa). The POD of the FE solution gives N ψ = 9 displacement modes with non-zero eigenvalue. It also provides N σ ψ = 48 admissible stress modes. We denote by (ψ σ k ) N σ ψ k=1 the statically admissible modes related to the POD of (σ n FE (·; μ 1 ) − σ n N ) N t n=1 . Some of these modes are shown in Figure 3. Figure 4 shows that the proposed error estimator incorporates errors related to the projection of the stress onto the dual reduced basis. In this example, if the number of modes of the dual reduced-basis is too small (i.e. smaller than 24), the error estimator does not account for the improvement of u HR when the reduced-basis (ψ k ) N ψ k=1 has more than 3 modes. The convenient choice of N σ ψ is set-up by considering its influence on the smallest value of η for large values of N ψ . Here, N σ ψ = 29 makes this influence negligible for all values of N ψ . When N σ ψ has been fixed, we can estimate c η according to Equation (66). Figure 5 shows the curve related to the HRAE and the error estimator c η η obtained for the sampling point μ = μ 1 , when varying the number of displacement modes for the maximization problem (66).
For the proposed example, inside the plastic zone, the stress field through the thickness of the beam is very sensitive to modifications of the parameters h and C. Here, the validity domain of the hyper-reduction approximations is defined by D = 0.05 2 N t n=1 (u n exp ) 2 . It is shown in Figure 6.  step. The speed-up obtained with this hyper-reduced model, including the error estimation, is 6. Without error estimation the speed-up is 8. Hence, the cost of error estimation is quite reasonable and the computational-complexity reduction is preserved. Here, the FE model is very simple. The bigger the FE model, the larger we expect the speed-up to be. The estimation of the validity domain is shown on the right of Figure 6. It is quite restrictive. This means that most of the points in D V are in the reference validity domain D V .

Figure 6
Error prediction on outputs. On the left, in blue, the reference validity domain D V over the parameter space D, on the right the approximate one D V . The white star is located at the sampling point μ 1 . The black dashed line is the limit of the reference validity domain defined by an error on the outputs lower than 5%.
Here, for small values of the parameter h, the prediction of the validity domain is too conservative. Such a situation can appear when the average influence of a parameter on the solution is larger than its influence on the output, because the error estimator is linked to the HRAE in an average sense.
We have tried the proposed numerical approach in case of N ψ = N c ψ (here N c ψ = 6), by using the same RID. The estimated validity domain and the boundary of the reference validity domain for the related HR model are shown in Figure 7. This modification has significant influence on the numerical results. Once again, all the points in D V are in the reference validity domain.
We have performed additional numerical simulations by substituting a Dirichlet boundary condition to the pointwise load. In this case, the time evolution of the load is:  ones because the reduced bases are different. They have a better accuracy with respect to N ψ = 7 and N σ ψ = 29. We can notice that the higher h and C, the smaller the plastic strains in the beam. Since the hyper-reduced model is accurate for low plastic strains, the validity domain covers the high values of h and C.