Strict upper bounds of the error in calculated outputs of interest for plasticity problems

This paper introduces a new computational technique for deriving guaranteed upper bounds of the error in calculated outputs of interest for plasticity problems under quasi-static conditions. This approach involves the resolution of some nonstandard additional problems, for which resolution techniques are given. All sources of error (spatial and temporal discretization, iteration stepping) are taken into account.


Introduction
Today more than ever, modeling and simulation are central to any mechanical engineering activity. A constant concern, both in industry and research, has always been the verification of models which, today, can attain very high levels of complexity. The novelty of the situation is that over the last twenty-five years truly quantitative tools for assessing the quality of a FE model have been put together into a branch of computational mechanics known as ''model verification''. Of course, the original continuum mechanics model remains the reference (Fig. 1). The state of the art on this topic is reviewed in [1][2][3][4][5].
Until 1990, only global error estimators, divided into three families introduced respectively by [30][31][32]were available. Since 1990, the evaluation of the quality of outputs of interest resulting from finite element analyses has become a key issue. This objective was beyond the reach of the early error estimators, which provided only global information that was totally insufficient for dimensioning purposes in mechanical design (where the dimensioning criteria involve local values of the stresses, displacements, stress intensity factors etc.). Among numerous relevant works concerning linear problems, the earliest were those of Peraire-Patera, Rannacher et al., Strouboulis-Babuška, Oden-Prudhomme, Ladevèze et al.; additional references can be found in [2][3][4][5]. The main idea which emerged then was that an output of interest can be expressed globally, enabling one to continue to use global error estimators, but an accurate error estimation is required for the finite element solution of the what is called the ''adjoint problem''. Extensions to the nonlinear case appeared in the early 90's with [6][7][8]; these approaches consisted in using linearization to revert back to the linear case during each time step.
Unfortunately, most of the estimates thus produced are not guaranteed to be upper or lower bounds, which is a very serious drawback as far as robust design is concerned. Consequently, the derivation of upper error bounds for the calculated values of outputs of interest has become one of today's pressing research challenges. Relatively few works have proposed answers, even in the linear case: [9][10][11][12]2,4,13,14].
Recent papers [15,16] have introduced a general method for deriving an upper bound of the error in a calculated output of interest for both linear and nonlinear time-dependent (quasi-static or dynamic) problems. This is probably the first published method leading to guaranteed upper error bounds for nonlinear problems and dynamic problems. Among the first developments, one can mention the error in pointwise quantities [17], the error in nonlinear goal-oriented quantities [18,19], and the error in outputs of interest in dynamics [20]. These methods are based on the concept of dissipation-type constitutive relation error (CRE) and on quasiexplicit techniques for the construction of the associated admissible FE solutions.
This paper introduces what is, to our knowledge, the first computable guaranteed error bounds ever published for goal-oriented quantities in plasticity. It is organized as follows. After this introduction, Section 2 presents the classical framework for plasticity problems and defines the concept of output of interest.
Section 3 discusses a version of the constitutive relation error which was already used in [11,21] and which involves admissible fields. This has the advantage of leading to a computable bound of the global discretization error. This is the first key point of our approach.
Section 4 introduces the second key point, which is the mirror problem and the associated basic identity which provides the difference between the exact value and the calculated value of the quantity of interest.
Section 5 presents the last key point (i.e. the bounding technique which combines the previous two points) and defines what we call the ''central problem'', for which an inexpensive calculation method is proposed in Section 5.
Finally, Section 6 presents a detailed example of the method. Initially, the structure being studied occupies a domain X bounded by @X. Let us assume small displacements, quasi-static loading and isothermal conditions, and let ½0; T be the time interval of interest. At any time t in ½0; T, the structure is placed in a given ''environment'' characterized by a displacement U d on a part @ u X of @X, a traction force density F d on @ f X (the part of @X which is complementary to @ u X), and a body force density f d within X.
The structure and its loading are illustrated in Fig. 1. The problem which describes the evolution of the structure over [0; T] is: Find the displacement field UðM; tÞ and the stress field rðM; tÞ, with t 2 ½0; T and M 2 X, which verify: the kinematic constraints: the equilibrium equations (principle of virtual work): and the constitutive relation: ðUÞ denotes the strain associated with the displacement ðUÞ ij ¼ 1 ; U ½0;T is the space containing the displace-ment field U defined over XÂ0; T½; S ½0;T is the space containing the stresses, also defined over XÂ0; T½; and U ad;0 is the vector space of the prescribed virtual velocities. Operator A, which is considered to be known and generally single-valued, characterizes the mechanical behavior of the material. Let U ½0;T ad denote the space of the displacement fields which verify the kinematic constraints (1).

Formulation in the thermodynamics framework
Using some global notations, the reference problem (1)-(3) can be rewritten in the framework of classical thermodynamics with internal variables. Let us introduce the generalized quantities, corresponding to three sets of variables: where p and e denote the inelastic strain and the elastic strain respectively; x is a vector of additional internal variables (e.g. hardening variables) and y the associated force vector (with x and y having the same dimension). One has ¼ e þ p , and the dissipation bilinear form is: Admissibility conditions Let us introduce the classical kinematic and static admissibility conditions (K-admissibility and S-admissibility): -_ e is said to be a K-admissible solution if it exists U such that -s is said to be an S-admissible solution if r 2 S ½0;T ; 8t 20; T½8U Ã 2 U ad;0 Constitutive relations The constitutive relations are divided into two parts: -the state equations: e e ¼ KðsÞ; -the state evolution equations: _ e p ¼ BðsÞ, where K and B denote material operators. K is assumed to be linear, symmetric and positive: this is the case for most elastic-(visco) plastic materials, generally after a change of internal variables (see [21]). This is called the ''normal'' material description. Let us observe that what we are doing here could be easily extended to the case where K derives from a convex potential. B can be nonlinear and multivalued, as in plasticity. However, we assume that we are dealing with a usual plastic models (where H-assumption holds). Thus, we have the following properties: B is a monotonous mutlivalued operator the state equations are separated; with r¼K e and y ¼ Hx HÀassumption: An important special case is the family of standard materials whose state evolution laws can be expressed using two potentials u and u Ã which are dual convex functions such that, for ðt; MÞ 2 ½0; T Â X. One can obtain: 8ð_ e p ; sÞ 2 S ½0;T u Ã ðsÞ þ uð_ e p Þ À s Á _ e p P 0; u Ã ðsÞ þ uð_ e p Þ À s Á _ e p ¼ 0 () _ e p ¼ BðsÞ; and BðsÞ ¼ @ s u Ã : The reference problem Using the generalized quantities, the reference problem becomes This is a rewritten form of Eqs. (1)- (3), where e and e e are auxiliary variables defined as functions of e p and s. Example: plasticity with isotropic hardening In plasticity with isotropic hardening, the additional internal variable is the accumulated plastic strain p (a scalar quantity, here p ¼ x). The associated force is the threshold R (also a scalar quantity, here R ¼ y). Then, the state equations are: and the evolution law B is multivalued and can be written as: This is a standard model whose corresponding potentials u and u Ã can be found in [4]: where w is the indicator function associated with the convex domain f _ e p 2 e ½0;T such that Tr p Â Ã ¼ 0 and k p k 6 _ pg, and w Ã is the indicator function associated with the convex domain fs 2 s ½0;T such that kr D k 6 Rg.

The output of interest
The output of interest a is a goal-oriented quantity, such as a mean stress value, or any internal variable component defined over an element or a set of elements and over a time interval. True local quantities in both time and space could also be considered. Such an output of interest can be expressed globally as whose extractors are ðd _ R ; dr R Þ over ½0; T Â X and where ''ex'' denotes the exact solution of the conceptual structural model of (9).
Here, for the sake of simplicity, nonlinear functions of the exact solution are not considered, but the corresponding extensions would involve no serious difficulty as long as the output of interest is local in space (see [17]).
Using generalized quantities, a ex is equal to The notation d indicates that the extractors must be interpreted as finite but small quantities relative to the solution of the conceptual structural problem. However, no linearization is involved.
Let us note that (11) can lead to outputs of interest which are related to the inelastic strain. For example, one has: Example. The following example is adopted in the remainder of the paper to illustrate the proposed error bounding method. As an illustration let us consider, under the plane stress assumption, the two-dimensional structure of Fig. 2 fixed along its lower left edge and subjected to a loading F d (t) along its upper right edge. We used the following characteristics for the problem: with the material parameters: The calculation was carried out with a spatial mesh composed of 8; 035 linear quadrilateral elements. The time interval ½0; T was divided into 20 steps. Fig. 2(d) shows the resulting approximate accumulated plastic strain at time T obtained with the FEM and a backward Euler scheme. The output of interest was defined as the mean plastic strain over domain x in the e y e y direction between times t 1 and t 2 . x was localized in the lower left corner, and t 1 and t 2 were chosen at the end of ½0; T (t 1 ¼ 0:85:T and t 2 ¼ T) as shown in Fig. 3. This output of interest can be expressed as Description of the reference problem (b), its geometry (a), its loading history (c), and the approximate finite element accumulated plastic strain p h at t ¼ T (d).
Considering its support, the output of interest can also be written as Remark 1. Any output of interest can be associated with an extractor.
Remark 2. If necessary in order to get a relatively small extractor, one introduces a sizing coefficient k; then, the goal-oriented quantity becomes:

The dissipation constitutive relation error and its upper error bounds
After reviewing the dissipation constitutive relation error, let us now derive upper bounds of the discretization error. This is the first of the three key points leading to guaranteed bounds for plasticity problems.

The admissible solution associated with the calculated solution of the reference problem
In this section, our objective is to associate an admissible solution with the calculated solution and with the data. In the dissipation error framework, ''admissible'' means that ð _ e p;h ;ŝ h Þ is such that: The approach used is well known and the construction is only outlined here; further details can be found in [4].
Having defined an admissible solution as one that verifies all the equations except the evolution laws, let us assume that the calculated solution was obtained using the FEM. Thus, at discrete times t m belonging to ½0; T, one knows ½ _ e h ; s h t ; t 2 ½0; t 1 ; . . . ; t n ¼ T: One assumes that ð _ e h ; s h Þ verifies the kinematic constraints and the equilibrium equations at t 2 ½0; t 1 ; . . . ; t n ¼ T in the FE sense. Considering that the behavior of the data during each time step is linear, one can extend this FE solution, originally defined at discrete times, across the whole time-space domain. This leads to ð _ e h ; s h Þ 2 S ½0;T which verifies the kinematic constraints and the equilibrium equations in the FE sense at any time t 2 ½0; T.
In order to get an associated admissible solution, let us begin with the same technique as that used in elasticity to define a displacement-stress pair ðÛ h ;r h Þ which verifies a set of kinematic constraints, equilibrium equations and initial conditions over ½0; T Â X. (Let us note that in the case of elastic-(visco) plastic behavior, with the constraint Tr _ p Â Ã ¼ 0, the previous displacement must be modified so that Tr _ p h i ¼ 0; for additional information on the application to the 3D case, see [22].) The additional internal variables ðx h ;ŷ h Þ, which must verify the state equations, can be obtained easily through the resolution of local problems which are related to the minimization of the dissipation-type constitutive relation error. Finally, one obtains an admissible solution ð _ e p ;ŝÞ 2 S ½0;T of the reference model. Further details can be found in [4,11]. Let us note that similar admissible solutions could also be derived using recent techniques such as [23][24][25][26].

The dissipation error
Compared to [15], we introduce a different version of the dissipation constitutive relation error in order to minimize the importance of the additional internal variables. Then: Hence: E T CRE ð_ e p ; sÞ ¼ 0: Ã

An upper error bound
The following theorem leads to an upper bound of global quantities involving the exact solution. This is the first key point of our approach.
Theorem 1. One has: with the initial conditions p Ã ¼ 0 and r Ã ¼ O at t ¼ 0: Fig. 4 shows the dissipation relation error of the approximate solution of the previous example. One can observe that the error is localized mainly in the vicinity of the geometric singularity and that its time evolution is quite homogeneous.

The basics of upper error bounds for goal-oriented verification
In this section, we derive the method which enables the calculation of strict upper bounds for goal-oriented verification in (visco) plastic problems. This development follows [15] with slight modifications associated with our choice of the dissipation error. The method is based on Theorem 1, the mirror problem and a duality property.

The reference mirror problem and the associated identity
The mirror problem, introduced in [15], is associated with the extractor and replaces the adjoint problem, with which it coincides in the linear case. Its main feature, inherited from the adjoint problem, is that it enables one to build an identity which gives the error in the quantity of interest, i.e. the difference between its exact value and its calculated value. The mirror problem is similar to the initial reference problem, except that time goes backward (s ¼ T À t) and modified evolution laws apply. This mirror problem can be expressed in terms of d-quantities as: Find the solution ðd _ e p ; dsÞ 2 S ½0;T such that The mirror problem can be solved using the same technique and the same discretization as for the initial problem. However, for extractors which are localized in space or in time, it may be preferable to use a discretization which is more refined locally. Let The interest of the mirror problem resides in the following identity, which was proven in [16], in which the mirror problem is interpreted as a perturbation problem.
where ð_ e p;ex ; s ex Þ is the exact solution of the reference problem; ð _ e p;hŝh Þ is an admissible solution of the reference problem; ðdê j p;h ; dŝ h Þ is an admissible solution of the mirror problem.
All the quantities on the right-hand side of (20) are known, except for the exact solution ( _ e p;ex ; s ex ) of the initial reference problem. In this case, one uses the identity given by Property 2, which is a direct consequence of (20). This Property 2 and the mirror problem are the second keypoint of our approach to get guaranteed error bounds.

dtdX:
Proof. The difference between two expressions of a ex À a h is

The upper error bound
This is the third and last keypoint of our approach. We describe here the method to get error upper bound for an output of interest a ex . Using Theorem 1, let us introduce: One also has: In order to get an upper bound, one begins by multiplying the previous expression of the error by a positive factor g 2 R þ , and then by adding the inequality (22): Then, the computable bound is a ex À a h 6 min gP0 1 jgj Taking Àg 2 R þ , one obtains a lower bound. Finally: where À h ¼ min In order to get an upper bound, one must solve a temporal problem at every point of the domain. An important contribution of this paper is that it proposes an inexpensive technique for the resolution of this problem which, from now on, we will call the ''central problem''.

The technical point: the resolution of the central problem
In the previous section, we introduced a method for the determination of strict upper bound of the error in calculated outputs of interest. This method requires the resolution of the maximization problem of (25), called the central problem. In the linear framework, as previously illustrated in [16] in visco-elasticity, this maximization problem can be solved explicitly using the global measure of the dissipation CRE. The framework of nonlinear materials, which is the subject of this paper, requires the introduction of a specific numerical method. The purpose of this section is first to study this central problem, its parameters and their magnitude, and then to propose a numerical method for its resolution, in the framework of standard plasticity model with isotropic hardening.

The functional of the central problem in plasticity with isotropic hardening
This functional is defined explicitly and divided into three parts which will be considered separately.
One has: Property 3. h ðgÞ has the following upper bound: where: Let us note that min kr Ã D k6kR Ã l 1 ð _ p Ã ; R Ã ; kr Ã D kÞ È É is equal to: In order to get l 2 ð _ p Ã ; R Ã Þ, one adds the quantity: This quantity is equal to zero because dR hjT ¼ 0 and DR Ã j0 ¼ 0:

Analysis of the order of magnitude of functional '
This analysis is presented in Appendix A. The different contributions to the error are studied, and the main conclusions are that in order to obtain sharp bounds. the extractors and the D terms must be small to get sharp bounds, which can always be achieved by an appropriate choice of the sizing parameter k; many terms of the functional can be neglected.

Resolution of the central problem
Using the solutions ðr h ; _ h Þ; ðr h ; _ h Þ and ðdr h ; d j h Þ, the problem can be written as a minimization problem over Q ½0;T : : Coefficients I; J and K are known and are small, which leads to the following local equations: Qð0; R Ã Þ 6 Q ðþ1; R Ã Þ < 1: Remark 5. Again in the case in which Q is piecewise linear, one starts with Q ð0; R h Þ, then introduces the solution as a sum of Dirac distributions Q ð _ p Ã 1 ; R Ã 1 Þ, and iterates if necessary.

First illustration
The reference problem and the output of interest are the same as in the previous example. The output of interest takes the value a h ¼ 4:19 10 À3 : Using the notations introduced in (12), one can define an upper bound of the local error denoted Da h : Table 1 The values of the lower and upper bounds as functions of parameter k (k > 0 or k < 0).  ja h þ a hh À a ex j 6 Da h : This bounding can be applied directly to the exact value of the output of interest, leading to its extremum: Let us introduce the nondimensional error measure: Table 1 shows the evolution of the bounds for different values of parameters (g; k). First, one can observe the significant role played by the sign of parameter k. This can be explained by considering the stress direction. The zone of interest is located in a bending zone. In particular, one can extract the range of each component of the stress field in that zone: In practice, based on this behavior, we implemented a procedure which enables the optimum parameters to be calculated inexpensively. Usually, only five computations are necessary to obtain an optimum parameter a À or a þ .
Remark 8. In practice, the properties of the mirror problem are similar to those of the adjoint problem, its equivalent in the linear case. The method we presented here had been previously implemented in the case of linear problems, for which the bounds can be calculated explicitly by means of the dissipation error. This has already be done for viscoelasticity [27,17,18,28] and for dynamics [20,29]. In the linear case, it was shown that the more accurate the resolution of the adjoint problem, the tighter the bounds. We showed that the same property also holds in the nonlinear case.

A proposal for improvement and conclusion
By analogy with the approach of [27], one could introduce a weighting function aðtÞ. Indeed, the temporal evolution of the global term introduced in the bounding procedure (see (23)) is quite homogeneous. Conversely, the behavior of the quantity which is associated with the error in the output of interest is local in time. If one introduces a function aðtÞ in order to balance the global quantities, the bounding relation becomes The choice of a piecewise constant aðtÞ leads to a new associated central problem which is similar to the previous one. The resolution of this problem could be obtained through an extension of the geometric approach. This function aðtÞ must be such that the evolution of the associated global term is comparable with that of the local term. We chose a simple piecewise constant function whose two values a 1 and a 2 are inspired by the quadratic behavior of the bounds, as in Fig. 8.
With this definition of aðtÞ, the bounding results were: a h þ a hh ¼ 4:07 Â 10 À3 and Da % h ¼ 27:31%: The work presented in this paper is based on the general method introduced in [15,16] for the calculation of strict error bounds of outputs of interest in quasi-static or dynamic problems involving any elastic-(visco) plastic material under small-displacement and convexity assumptions. Applications of this theoretical approach had been previously presented for viscoelastic materials [27,20].  This work is a first attempt at using this general framework to derive guaranteed error bounds for plasticity problems. This proposed method was found to be promising. It could be improved [17,18,20] by refining the calculation of the mirror problem in order to obtain more accurate bounds. Finally, it could be extend to nonlinear quantity of interest, such as the maximum Von Mises equivalent stress or the cumulated plastic strain. For linear problems including viscoelastic ones such a result has already been obtained introducing Dirac distributions for the extractors. For the second pattern, the plastic strain can be expressed as dp ¼ sup s 0 6s kP r D;h ðdr Djs 0 Þ þ Oð 0 Þk; 0 n o , which leads to small quantities.
In the general case, the D terms are often small. Although the mirror problem is applicable to any type of extractor, it is usually better to work with small quantities. In this case, the choice of parameter k enables this property to be always satisfied.