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

Pierre Ladevèze, Benoît Blaysat, Éric Florentin

► To cite this version:

Pierre Ladevèze, Benoît Blaysat, Éric Florentin. Strict upper bounds of the error in calculated outputs of interest for plasticity problems. Computer Methods in Applied Mechanics and Engineering, 2012, 245-246, pp.194-205. 10.1016/j.cma.2012.07.009. hal-01647872

HAL Id: hal-01647872

https://hal.science/hal-01647872

Submitted on 13 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Strict upper bounds of the error in calculated outputs of interest for plasticity problems

Pierre Ladevèze<sup>a,b,*</sup>, Benoît Blaysat<sup>a</sup>, Éric Florentin<sup>a</sup>

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.

## 1. 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–5].

Until 1990, only global error estimators, divided into three families introduced respectively by [30–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–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–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–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 quasi-explicit 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

<footnote><sup>a</sup>LMT-Cachan (ENS-Cachan/CNRS/UPMC/UniverSud Paris Pres), 61 Avenue du Président Wilson, 94230 Cachan, France</footnote>

<footnote><sup>b</sup>EADS Foundation Chair of Advanced Computational Structural Mechanics</footnote>

<footnote>* Corresponding author at: LMT-Cachan (ENS-Cachan/CNRS/UPMC/UniverSud Paris Pres), 61 Avenue du Président Wilson, 94230 Cachan, France.
E-mail address: pierre.ladeveze@lmt.ens-cachan.fr (P. Ladevèze).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic representation of the environment (i.e. the prescribed conditions).</figcaption></figure>

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.

## 2. The reference problem and its outputs of interest

### 2.1. The conceptual structural model

#### 2.1.1. Standard formulation

Initially, the structure being studied occupies a domain $\Omega$ bounded by $\partial\Omega$. 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 $\partial_u\Omega$ of $\partial\Omega$, a traction force density $F_d$ on $\partial_f\Omega$ (the part of $\partial\Omega$ which is complementary to $\partial_u\Omega$), and a body force density $f_d$ within $\Omega$. 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 $\sigma(M, t)$, with $t \in [0; T]$ and $M \in \Omega$, which verify:

* the kinematic constraints:

$$ U \in \mathcal{U}^{[0:T]}; \quad U_{|\partial_u \Omega} = U_d \text{ over } [0, T]; \qquad (1) $$

* the equilibrium equations (principle of virtual work):

$$ \sigma \in S^{[0:T]}; \quad \forall t \in [0, T], \ T[ \ \forall U^* \in \mathcal{U}_{ad,0} - \int_{\Omega} \{\text{Tr}[\sigma\epsilon(U^*)]\}d\Omega \\ + \int_{\Omega} \{f_d \cdot U^*\}d\Omega + \int_{\partial_f\Omega} \{F_d \cdot U^*\}dS = 0 \qquad (2) $$

* and the constitutive relation:

$$ \forall t \in [0; T] \quad \forall M \in \Omega \quad \sigma|_t = A(\epsilon(\dot{U})|_t; \tau \le t). \qquad (3) $$

$\epsilon(U)$ denotes the strain associated with the displacement $(\epsilon(U)|_{ij} = \frac{1}{2}(U_{ij} + U_{j,i}); U^{[0:T]})$ is the space containing the displacement

field $U$ defined over $\Omega \times [0, T]; S^{[0:T]}$ is the space containing the stresses, also defined over $\Omega \times [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_{ad}^{[0:T]}$ denote the space of the displacement fields which verify the kinematic constraints (1).

#### 2.1.2. 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:

$$ s = \{ \boldsymbol{\sigma}; \mathbf{y} \} \quad \dot{\boldsymbol{\epsilon}}_p = \{ \dot{\boldsymbol{\epsilon}}_p; -\dot{\mathbf{x}} \} \quad \dot{\boldsymbol{\epsilon}}_e = \{ \dot{\boldsymbol{\epsilon}}_e; \dot{\mathbf{x}} \}, \qquad (4) $$

where $\epsilon_p$ and $\epsilon_e$ denote the inelastic strain and the elastic strain respectively; $\mathbf{x}$ is a vector of additional internal variables (e.g. hardening variables) and $\mathbf{y}$ the associated force vector (with $\mathbf{x}$ and $\mathbf{y}$ having the same dimension). One has $\epsilon = \epsilon_e + \epsilon_p$, and the dissipation bilinear form is:

$$ d : e^{[0:T]} \times s^{[0:T]} \to \mathfrak{R}, \\ (\dot{\epsilon}_p, s) \to \int_0^T \{s \cdot \dot{\epsilon}_p\}dt = \int_0^T \{\text{Tr}[\boldsymbol{\sigma}\dot{\epsilon}_p] - \mathbf{y} \cdot \dot{\mathbf{x}}\}dt \qquad (5) $$

* Admissibility conditions

Let us introduce the classical kinematic and static admissibility conditions (K-admissibility and S-admissibility):

- $\dot{\mathbf{e}}$ is said to be a K-admissible solution if it exists $U$ such that

$$ \dot{\mathbf{e}} = \{ \dot{\epsilon}(U); \mathbf{0} \}, \quad \text{where} \quad U \in \mathcal{U}^{[0:T]}; \quad U_{|\partial_u \Omega} = U_d \text{ on } [0, T], \qquad (6) $$

- $\mathbf{s}$ is said to be an S-admissible solution if

$$ \left\{ \begin{array}{l} \sigma \in S^{[0:T]}; \forall t \in [0, T] \ \forall U^* \in \mathcal{U}_{ad,0} \\ - \int_{\Omega} \{\text{Tr}[\sigma\epsilon(U^*)]\}d\Omega + \int_{\Omega} \{f_d \cdot U^*\}d\Omega + \int_{\partial_f\Omega} \{F_d \cdot U^*\}dS = 0 \end{array} \right. \qquad (7) $$

* Constitutive relations

The constitutive relations are divided into two parts:

- the state equations: $e_e = \Lambda(s)$;

- the state evolution equations: $\dot{\epsilon}_p = \mathbf{B}(s)$,

where $\mathbf{A}$ and $\mathbf{B}$ denote material operators.

$\mathbf{A}$ 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 $\mathbf{A}$ derives from a convex potential.

$\mathbf{B}$ can be nonlinear and multivalued, as in plasticity.

However, we assume that we are dealing with a usual plastic models (where $\mathcal{H}$-assumption holds). Thus, we have the following properties:

$$ \left. \begin{array}{l} \bullet \mathbf{B} \text{ is a monotonous multivalued operator} \\ \bullet \text{the state equations are separated, with} \\ \qquad \boldsymbol{\sigma} = \mathbf{K}\boldsymbol{\epsilon}_e \text{ and } \mathbf{y} = \mathbf{H}\mathbf{x} \\ \text{where } \mathbf{K} \text{ and } \mathbf{H} \text{ are positive definite material operators.} \end{array} \right\} \mathcal{H}\text{-assumption.} $$

An important special case is the family of standard materials whose state evolution laws can be expressed using two potentials $\varphi$ and $\varphi^*$ which are dual convex functions such that, for $(t, M) \in [0, T] \times \Omega$. One can obtain:
$$
\begin{gathered}
\forall (\dot{e}_p, s) \in S^{[0,T]} \quad \varphi^*(s) + \varphi(\dot{e}_p) - s \cdot \dot{e}_p \ge 0, \\
\varphi^*(s) + \varphi(\dot{e}_p) - s \cdot \dot{e}_p = 0 \iff \dot{e}_p = \mathbf{B}(s), \\
\text{and } \mathbf{B}(s) = \partial_s \varphi^*.
\end{gathered}
\tag{8}
$$

• The reference problem

Using the generalized quantities, the reference problem
becomes

Find $(\dot{e}_p, s) \in S^{[0,T]}$ such that

• $\dot{e} = \dot{e}_e + \dot{e}_p$ K-admissible,

• $(\dot{e}, s)$ S-admissible,

• $e_e = \Lambda(s)$ (state equations),

• $\dot{e}_p = \mathbf{B}(s)$ (state evolution laws),

• $s = 0, e = 0, \underline{U} = 0, \dot{\underline{U}} = 0$ at $t = 0$ (initial conditions). \hfill (9)

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 = \mathbf{x}$). The associated force is the threshold $R$ (also a scalar
quantity, here $R = \mathbf{y}$). Then, the state equations are:

$R = k \cdot p + R_0$ with $k \ge 0$ and $R_0 \ge 0$;

$\boldsymbol{\sigma} = \mathbf{K} \boldsymbol{\epsilon}_e$

and the evolution law **B** is multivalued and can be written as:

$$
\begin{align}
\{ \dot{e}_p; -\dot{p} \} &= \dot{p} \{ \frac{\sigma_D}{R}; -1 \} && \text{with } \dot{p} \ge 0, \nonumber \\
f(\boldsymbol{\sigma}, R) &= \| \boldsymbol{\sigma}_D \| - R \le 0 && \text{and } \dot{p} \cdot f(\boldsymbol{\sigma}, R) = 0. \tag{10}
\end{align}
$$

This is a standard model whose corresponding potentials $\varphi$ and $\varphi^*$ can be found in [4]:

$$
\varphi(\dot{e}_p) = R_0 \dot{p} + \psi,
$$

$$
\varphi^*(s) = \psi^*,
$$

where $\psi$ is the indicator function associated with the convex domain $\{\dot{e}_p \in e^{[0,T]}$ such that $\mathrm{Tr}[\epsilon_p] = 0$ and $\|\epsilon_p\| \le \dot{p}\}$, and $\psi^*$ is the indicator function associated with the convex domain $\{s \in s^{[0,T]}$ such that $\|\boldsymbol{\sigma}_D\| \le R\}$.

2.2. The output of interest

The output of interest $\alpha$ 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

$$
\alpha_{ex} = \int_{\Omega} \int_{0}^{T} \{ \mathrm{Tr} [ \boldsymbol{\sigma}_{ex} \delta \dot{\boldsymbol{e}}_{\Sigma} ] + \mathrm{Tr}[ \dot{\epsilon}_{ex} \delta \bar{\boldsymbol{\sigma}}_{\Sigma} ] \} dt d\Omega, \quad (11)
$$

whose extractors are $(\delta\dot{\epsilon}_{\Sigma}, \delta\bar{\sigma}_{\Sigma})$ over $[0; T] \times \Omega$ 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, $\alpha_{ex}$ is equal to

$$
\begin{align}
\alpha_{ex} &= \int_{\Omega} \int_{0}^{T} \{ s_{ex} \cdot \delta \dot{e}_{\Sigma} + \dot{e}_{ex} \cdot \delta \bar{s}_{\Sigma} \} dt d\Omega \quad \text{with} \nonumber \\
\delta \dot{e}_{\Sigma} &= \{ \delta \dot{e}_{\Sigma}; \mathbf{0} \} \quad \text{and} \quad \delta \bar{s}_{\Sigma} = \{ \delta \bar{\sigma}_{\Sigma}; \mathbf{0} \}. \tag{12}
\end{align}
$$

The notation $\delta$ indicates that the extractors must be interpreted
as finite but small quantities relative to the solution of the concep-
tual 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:

$$
\begin{gather*}
\alpha_{ex} = \int_{\Omega} \int_{0}^{T} \{ \mathrm{Tr}[\epsilon_{p,ex} \delta \tilde{\mathbf{A}}_{\Sigma}] \} dt d\Omega, \\
\text{with } \delta \dot{\epsilon}_{\Sigma} = -\mathbf{K} \delta \tilde{\mathbf{A}}_{\Sigma}, \delta \bar{\sigma}_{\Sigma} = \int_{t}^{T} \{ \delta \tilde{\mathbf{A}}_{\Sigma} \} d\tau.
\end{gather*}
$$

**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:

$$
L = 100 \text{ mm}; \quad T = 1 \text{ s}; \quad F_0 = 0, 15 \text{ N mm}^{-2},
$$

with the material parameters:

$$
R_0 = 1 \text{ MPa}; \quad k = 8.16 \text{ MPa}; \quad E = 244.95 \text{ MPa}; \quad \nu = 0.3.
$$

The calculation was carried out with a spatial mesh composed of
8,035 linear quadrilateral elements. The time interval $[0; T]$ was di-
vided into 20 steps. Fig. 2(d) shows the resulting approximate accu-
mulated 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 $\omega$ in the $e_y \otimes e_y$ direction between times $t_1$ and $t_2$. $\omega$ 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 \cdot T$ and $t_2 = T$) as shown in Fig. 3.

This output of interest can be expressed as

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** 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).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Definition of the zone of interest in space and in time. $\omega$ is an $l \times l$ square, with $l = L/5$ (left); the associated extractors are shown on the right.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

$$ \alpha_{ex} = \frac{1}{\text{mes}(\omega)(t_2 - t_1)} \int_{\omega} \int_{t_1}^{t_2} \{ \epsilon_{p,ex,yy} \} dt d\Omega \quad \text{where } \omega \subset \Omega \text{ and} \\ 0 \le t_1 < t_2 \le T. $$

Considering its support, the output of interest can also be written as

$$ \alpha_{ex} = \int_{\Omega} \int_{0}^{T} \{ \text{Tr}[\epsilon_{p,ex} \delta \tilde{\mathbf{A}}_{\Sigma}] \} dt d\Omega \quad \text{with} \\ \delta \tilde{\mathbf{A}}_{\Sigma} = \begin{cases} \frac{\epsilon_y \otimes \epsilon_y}{\text{mes}(\omega)(t_2 - t_1)}, & \text{for } \forall (\mathbf{M}, t) \in \omega \times [t_1; t_2] \\ \emptyset, & \text{otherwise} \end{cases} \quad (13) $$

**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 $\lambda$; then, the goal-oriented quantity becomes:

$$ \lambda \alpha_{ex} = \int_{\Omega} \int_{0}^{T} \left\{ \text{Tr}\left[\boldsymbol{\sigma}_{ex}\lambda \cdot \dot{\boldsymbol{\delta}\hat{\mathbf{\epsilon}}}_{\Sigma}\right] + \text{Tr}\left[\hat{\mathbf{e}}_{ex}\lambda \cdot \dot{\boldsymbol{\delta}\tilde{\mathbf{\sigma}}}_{\Sigma}\right] \right\} dt d\Omega. \quad (14) $$

### 3. 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.

#### 3.1. 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 $(\dot{\hat{e}}_{p,h}, \hat{s}_h)$ is such that:

* $\hat{e}_h$: K-admissible, with $\hat{e}_h = \hat{e}_{e,h} + \hat{e}_{p,h}$;

* $\hat{s}_h$: S-admissible;

* $\hat{s}_h = \Lambda \hat{e}_{e,h}$.

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

$$ [\hat{e}_h, s_h]_t, \quad t \in [0, t_1, \dots, t_n = T]. $$

One assumes that $(\hat{e}_h, s_h)$ verifies the kinematic constraints and the equilibrium equations at $t \in [0, t_1, \dots, 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 $(\hat{e}_h, s_h) \in S^{[0:T]}$ which verifies the kinematic constraints and the equilibrium equations in the FE sense at any time $t \in [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 $(\hat{\mathbf{J}}_h, \hat{\boldsymbol{\sigma}}_h)$ which verifies a set of kinematic constraints, equilibrium equations and initial conditions over $[0; T] \times \Omega$. (Let us note that in the case of elastic-(visco) plastic behavior, with the constraint Tr[$\hat{\mathbf{e}}_p$] = 0, the previous displacement must be modified so that Tr[$\hat{\mathbf{e}}_p$] = 0; for additional information on the application to the 3D case, see [22].) The additional internal variables $(\hat{\mathbf{x}}_h, \hat{\mathbf{y}}_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 $(\hat{\mathbf{e}}_p, \hat{s}_h) \in 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–26].

#### 3.2. 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.

**Definition 1.** The local (in terms of space) dissipation CRE $E_{CRE}^T$ is

$$ E_{CRE}^{T}(\dot{\hat{e}}_p, s, \underline{M}) = \max_{(\hat{e}_p^*, s^*) \in \Gamma^{[0:T]}} \left\{ - \int_{0}^{T} \int_{0}^{t} \left\{ \text{Tr}\left[ (\boldsymbol{\sigma}(t, \underline{M}) - \boldsymbol{\sigma}^*(t, \underline{M})) \times (\dot{\hat{e}}_p(t, \underline{M}) - \dot{\hat{e}}_p^*(t, \underline{M})) \right] \right\} \frac{d\tau dt}{T} \right\}, $$

where $\Gamma^{[0:T]}$ is the set of the pairs $(\hat{e}_p^*, s^*) \in S^{[0:T]}$ which verify the initial condition as well as the behavior.

$$ (\hat{e}_p^*, s^*) \in \Gamma^{[0:T]} \iff \left\{ \begin{array}{l} (\hat{e}_p^*, s^*) \in S^{[0:T]}; \\ \text{the evolution law } \dot{\hat{e}}_p^* = \mathbf{B}(s^*); \\ \text{the state law } e_e = \Lambda(s); \\ \text{the initial conditions } e_p^* = 0 \text{ and } s^* = 0 \text{ at } t=0. \end{array} \right. $$

**Definition 2.** The global dissipation error $E_{CRE}^{T, \Omega}$ is
$$E_{\text{CRE}}^{T}(\dot{e}_p, s) = \int_{\Omega} \{E_{\text{CRE}}^{T}(\dot{e}_p, s, \underline{M})\} d\Omega.$$

**Property 1.** For $(\dot{e}_p, s) \in S^{[0:T]}$, one has:

(i) $E_{\text{CRE}}^{T}(\dot{e}_p, s) \ge 0$;

(ii) Under the $\mathcal{H}$-assumption, $E_{\text{CRE}}^{T}(\dot{e}_p, s) = 0$ if $(\dot{e}_p, s) \in \Gamma^{[0:T]}$.

**Proof.** (i) is straightforward.

For (ii), let $(\dot{e}_p, s)$ and $(\dot{e}_p^*, s^*)$ be two elements of $\Gamma^{[0:T]}$; since $\mathbf{B}$ is monotonous, one has:

$$ 
\begin{aligned}
q &= \int_{0}^{T} \int_{0}^{t} \{ (\dot{e}_p^* - \dot{e}_p)(s^* - s) \} \frac{d\tau dt}{T} \ge 0 \\
&= \int_{0}^{T} \int_{0}^{t} \{ \text{Tr}[(\dot{e}_p^* - \dot{e}_p)(\boldsymbol{\sigma}^* - \boldsymbol{\sigma})] - (\dot{\mathbf{x}}^* - \dot{\mathbf{x}})(\mathbf{y}^* - \mathbf{y}) \} \frac{d\tau dt}{T} \ge 0.
\end{aligned}
\quad (15) $$

Then:

$$ 
\begin{aligned}
& - \int_{0}^{T} \int_{0}^{t} \{ \text{Tr}[(\dot{e}_p^* - \dot{e}_p)(\boldsymbol{\sigma}^* - \boldsymbol{\sigma})] \} \frac{d\tau dt}{T} \\
& \le -\frac{1}{2} \int_{0}^{T} \{ (\mathbf{x}^* - \mathbf{x})(\mathbf{y}^* - \mathbf{y}) \} \frac{dt}{T} \le 0.
\end{aligned}
\quad (16) $$

Hence:

$$E_{\text{CRE}}^{T}(\dot{e}_p, s) = 0. \quad \square$$

### 3.3. 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:

$$ 
\begin{aligned}
& (i) - \int_{\Omega} \int_{0}^{T} \int_{0}^{t} \{ \text{Tr}[(\boldsymbol{\sigma}_{\text{ex}} - \hat{\boldsymbol{\sigma}}_{h})(\dot{\boldsymbol{e}}_{p,\text{ex}} - \dot{\boldsymbol{e}}_{p,h})] \} \frac{d\tau dt d\Omega}{T} \\
& = \frac{1}{2} \int_{\Omega} \int_{0}^{T} \{ \text{Tr}[(\boldsymbol{\sigma}_{\text{ex}} - \hat{\boldsymbol{\sigma}}_{h})\mathbf{K}^{-1}(\boldsymbol{\sigma}_{\text{ex}} - \hat{\boldsymbol{\sigma}}_{h})] \} \frac{d\tau dt d\Omega}{T}
\end{aligned}
$$

$$ (ii) - \int_{\Omega} \int_{0}^{T} \int_{0}^{t} \{ \text{Tr}[(\boldsymbol{\sigma}_{\text{ex}} - \hat{\boldsymbol{\sigma}}_{h})(\dot{\boldsymbol{e}}_{p,\text{ex}} - \dot{\boldsymbol{e}}_{p,h})] \} \frac{d\tau dt d\Omega}{T} \le E_{\text{CRE}}^{T}(\dot{e}_{p,h}, \hat{s}_{h}) $$

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** The dissipation relation error $\epsilon$ associated with the approximate solution: map (left) and time evolution (right).</figcaption></figure>

**Proof.**

(i) follows from the admissibility of $(\dot{e}_{p,h}, \hat{s}_h)$ and $(\dot{e}_{p,\text{ex}}, \hat{s}_{\text{ex}})$;

(ii) is a consequence of Definition 2. $\square$

### 3.4. Example: plasticity with isotropic hardening

$$ 
\Gamma^{[0:T]} = 
\begin{cases}
\dot{e}_p^* = \{\dot{e}_p^*; -\dot{p}^*\}, & s^* = \{\boldsymbol{\sigma}^*; R^*\} \ (\dot{e}_p^*, s^*) \in S^{[0:T]}; \\
\dot{e}_p^* = \dot{p} \frac{\boldsymbol{\sigma}_D^*}{\|\boldsymbol{\sigma}_D^*\|}; \dot{p}^* \le 0; \|\boldsymbol{\sigma}_D^*\| - R^* \le 0 \text{ and } \dot{p}^* (\|\boldsymbol{\sigma}_D^*\| - R^*) = 0 \\
R^* = k \cdot p^* + R_0; \\
\text{with the initial conditions } p^* = 0 \text{ and } \boldsymbol{\sigma}^* = \mathbf{0} \text{ at } t=0.
\end{cases}
\quad (17) $$

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.

## 4. 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.

### 4.1. 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 ($\tau = T - t$) and modified evolution laws apply. This mirror problem can be expressed in terms of $\delta$-quantities as:

$$ 
\text{Find the solution } (\delta\dot{\hat{e}}_p, \delta\hat{s}) \in S^{[0:T]} \text{ such that } 
\begin{cases}
\delta\dot{\hat{e}} = \delta\dot{\hat{e}}_e + \delta\dot{\hat{e}}_p; (\delta\dot{\hat{e}} - \delta\dot{\hat{e}}_\Sigma) \text{ K-admissible to zero} \\
(\delta\hat{s} - \delta\hat{s}_\Sigma) \text{ S-admissible to zero} \\
\delta\hat{e}_e = \Lambda(\delta\hat{s}) \text{ (state equations)} \\
\bullet \quad \delta\hat{e}_p = \tilde{\mathbf{B}}(\delta\hat{s}) \text{ (evolution laws)} \\
\delta = 0, \delta\hat{s} = 0, \delta\underline{U} = 0 \text{ (initial conditions at } \tau = 0)
\end{cases}
\quad (18) $$
<figure><img src="image_7.png" /><figcaption>Fig. 5. The approximate finite element solution of the mirror problem in terms of the accumulated plastic strain ($\delta p_h$): solution map at time $t = 0$ (left) and time evolution at a Gauss point located in the lower left corner (right).</figcaption></figure>

From operator $\mathbf{B}$ of the conceptual model, $\tilde{\mathbf{B}}$ can be taken as:

$$ \delta \hat{\mathbf{e}}_p = \tilde{\mathbf{B}}(\delta \tilde{\mathbf{S}}) \equiv \mathbf{B}(s_h + \delta \tilde{\mathbf{S}}) - \mathbf{B}(s_h), \quad (19) $$

where “■” denotes the reverse-time derivative ($\hat{f} = \frac{\partial f}{\partial t}$) and $s_h(\tau)$ is the calculated solution of the reference problem. With this definition of the central problem, the extractors are defined taking into account the retrograde aspect, i.e. $\delta \dot{\tilde{\mathbf{e}}}_{\Sigma}$ and $\delta \tilde{\mathbf{S}}_{\Sigma}$ are functions of $\tau = (T - t)$.

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 $\langle \delta \hat{\mathbf{e}}_{p,h}, \delta \tilde{\mathbf{S}}_h \rangle$ be the calculated solution of the mirror problem and

$\langle \delta \hat{\mathbf{e}}_{p,h}, \delta \tilde{\mathbf{S}}_h \rangle$ the associated admissible solution.

Fig. 5 shows some properties of the calculated solution of the output of interest considered.

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.

$$ (\alpha_h - \alpha_{ex}) = \int_{\Omega} \int_{0}^{T} \left\{ (\dot{\mathbf{e}}_{p,ex} - \dot{\hat{\mathbf{e}}}_{p,h}) \cdot \delta \hat{\mathbf{S}}_h - (s_{ex} - s_h) \delta \hat{\mathbf{e}}_{p,h} \right\} dt d\Omega, \quad (20) $$

where

* $(\dot{\mathbf{e}}_{p,ex}, s_{ex})$ is the exact solution of the reference problem;

* $(\dot{\mathbf{e}}_{p,h} \delta \tilde{\mathbf{S}}_h)$ is an admissible solution of the reference problem;

* $(\delta \hat{\mathbf{e}}_{p,h}, \delta \tilde{\mathbf{S}}_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 $(\dot{\mathbf{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.

**Property 2.**

$$ (\alpha_h - \alpha_{ex}) = \int_{\Omega} \int_{0}^{T} \left\{ \mathrm{Tr} \left[ (\dot{\mathbf{e}}_{p,ex} - \dot{\hat{\mathbf{e}}}_{p,h}) \cdot \delta \hat{\mathbf{\sigma}}_h - (\boldsymbol{\sigma}_{ex} - \boldsymbol{\sigma}_h) \delta \hat{\mathbf{e}}_{p,h} \right] \right\} dt d\Omega. $$

**Proof.** The difference between two expressions of $\alpha_{ex} - \alpha_h$ is

$$ \begin{align*} z &\equiv \int_{\Omega} \int_{0}^{T} \left\{ -(\dot{\mathbf{x}}_{ex} - \dot{\hat{\mathbf{x}}}_{h}) \cdot \delta \hat{\mathbf{y}}_{h} + (\mathbf{y}_{ex} - \mathbf{y}_{h}) \delta \hat{\mathbf{x}}_{h} \right\} dt d\Omega \\ &= \int_{\Omega} \left\{ - \left[ (\mathbf{x}_{ex} - \hat{\mathbf{x}}_{h}) \cdot \delta \hat{\mathbf{y}}_{h} \right]_{0}^{T} \right\} d\Omega = 0 \quad \text{as} \quad \delta \hat{\mathbf{y}}|_{t=T} \\ &= 0 \quad \text{and} \quad (\mathbf{x}_{ex} - \hat{\mathbf{x}}_{h})|_{t=0} = 0. \quad \square \end{align*} $$

## 4.2. 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 $\alpha_{ex}$. Using Theorem 1, let us introduce:

$$ \gamma(\dot{\mathbf{e}}_{p,ex}, s_{ex}) = - \int_{0}^{t} \{ \mathrm{Tr}[(\boldsymbol{\sigma}_{ex} - \hat{\boldsymbol{\sigma}}_{h})(\dot{\mathbf{e}}_{p,ex} - \dot{\mathbf{e}}_{p}, h)] \} \frac{d\tau}{T}. \quad (21) $$

One also has:

$$ 0 \le \int_{\Omega} \int_{0}^{T} \{ \gamma(\dot{\mathbf{e}}_{p,ex}, s_{ex}) \} + E_{CRE}^{T, \Omega}. \quad (22) $$

In order to get an upper bound, one begins by multiplying the previous expression of the error by a positive factor $\eta \in \mathbb{R}^{+}$, and then by adding the inequality (22):

$$ \alpha_{ex} - \alpha_h \le \frac{1}{\eta} \left[ \int_{\Omega} \left\{ \eta \int_{0}^{T} \left\{ \mathrm{Tr} \left[ (\dot{\mathbf{e}}_{p,ex} - \dot{\hat{\mathbf{e}}}_{p,h}) \cdot \delta \hat{\boldsymbol{\sigma}}_{h} - (\boldsymbol{\sigma}_{ex} - \boldsymbol{\sigma}_{h}) \delta \hat{\mathbf{e}}_{p,h} \right] \right. \right. \right. \\ \left. \left. - \gamma(\dot{\mathbf{e}}_{p,ex}, s_{ex}) \right\} dt \right\} d\Omega + E_{CRE}^{T, \Omega} \right]. \quad (23) $$

Then, the computable bound is

$$ \alpha_{ex} - \alpha_h \le \min_{\eta \ge 0} \underbrace{\left\{ \frac{1}{|\eta|} \left[ \int_{\Omega} \{g_h(\eta)\} d\Omega + E_{CRE}^{T,\Omega} \right] \right\}}_{\epsilon_h(\eta)} = \epsilon_h^+, \quad (24) $$

where

$$ \forall \underline{M} \in \Omega, \quad g_h(\eta) = \max_{(\dot{\mathbf{e}}_{p,h}^*, s^*) \in \Gamma^{[0,T]}} \left\{ \eta \int_{0}^{T} \left\{ \mathrm{Tr} \left[ (\dot{\mathbf{e}}_{p,ex} - \dot{\hat{\mathbf{e}}}_{p,h}) \cdot \delta \hat{\boldsymbol{\sigma}}_{h} - (\boldsymbol{\sigma}_{ex} - \boldsymbol{\sigma}_{h}) \delta \hat{\mathbf{e}}_{p,h} \right] \right\} dt - \int_{0}^{T} \left\{ \gamma(\dot{\mathbf{e}}_{p,h}^*, s^*) \right\} dt \right\}. \quad (25) $$

Taking $-\eta \in \mathbb{R}^{+}$, one obtains a lower bound. Finally:

$$ |\alpha_{ex} - \alpha_h - \alpha_{hh}| \le \frac{1}{2}(\epsilon_h^+ + \epsilon_h^-), \quad (26) $$

where

$$ \begin{align*} \epsilon_h^- &= \min_{\eta \ge 0} \left\{ \frac{1}{\eta} \left[ \int_{\Omega} \{g_h(-\eta)\} d\Omega + E_{CRE}^{T,\Omega} \right] \right\} \\ \alpha_{hh} &= \frac{1}{2}(\epsilon_h^+ - \epsilon_h^-). \end{align*} $$

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”.
## 5. 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.

### 5.1. 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.** $\epsilon_h(\eta)$ has the following upper bound:

$$ \epsilon_h(\eta) \le - \min_{(\hat{p}^*, R^*) \in Q^{(0,T)}} \left\{ \int_0^T \int_0^T \{\ell(\hat{p}^*, R^*)\} \frac{d\tau dt}{T} \right\}, $$

where:

$$ \bullet \ Q^{(0:T)} \equiv \{(\hat{p}, R) \mid p \in L^2([0; T]), \hat{p} \ge 0, R = R_0 + k \cdot p\}. $$

$$ \bullet \ \ell(\hat{p}^*, R^*) = \ell_1(\hat{p}^*, R^*) + \ell_2(\hat{p}^*, R^*) + \ell_3(R^*). $$

$$ \bullet \ \ell_1(\hat{p}^*, R^*) = -\|R^*\| \cdot \left( \dot{\hat{\sigma}}_{p,h} + \frac{\hat{\sigma}_h}{2\mu \cdot (T-\tau)} + \Delta \dot{\hat{\epsilon}}_p - \delta \hat{\epsilon}_{p,h} \right) \\ + \dot{p}^* \cdot \left( \hat{\sigma}_h + \Delta \sigma + \delta \hat{\sigma}_h \right) + R^* \cdot \left( \dot{p}_h + \frac{P_{\sigma_{D,h}}(\hat{\sigma}_{D,h})}{2\mu \cdot (T-\tau)} + P_{\sigma_{D,h}}(\Delta \dot{\hat{\epsilon}}_p) - P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) \right) \\ + \dot{p}^* \cdot (\|\hat{\sigma}_{D,h}\| + P_{\sigma_{D,h}}(\Delta \sigma) + P_{\sigma_{D,h}}(\delta \hat{\sigma}_h)). $$

$$ \bullet \ \ell_2(\hat{p}^*, R^*) = R^* \cdot \dot{p}^* + \operatorname{Tr}[\hat{\sigma}_h \cdot \dot{\hat{\epsilon}}_{p,h}] - \dot{p}^* \cdot \left( \|\hat{\sigma}_{D,h}\| + P_{\sigma_{D,h}}(\Delta \sigma) + P_{\sigma_{D,h}}(\delta \hat{\sigma}_h) \right) - R^* \cdot \dot{p}_h + P_{\sigma_{D,h}}(\Delta \dot{\hat{\epsilon}}_p) - P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) \\ + \operatorname{Tr}[\dot{\hat{\epsilon}}_{p,h} \cdot \delta \hat{\sigma}_h - \hat{\sigma}_h \cdot \delta \hat{\epsilon}_{p,h}] + \Delta \dot{p}^* \cdot \delta \hat{\underline{R}}_h - \Delta R^* \cdot \delta \hat{\underline{p}}_h. $$

$$ \bullet \ \ell_3(R^*) = \frac{(R^* - \|\hat{\sigma}_D^*\|)(R_1 - \|\hat{\sigma}_D^*\|)}{2\mu \cdot (T-\tau)} + \frac{(\|\hat{\sigma}_D^*\| - \|\hat{\sigma}_{D,h}\|)^2}{4\mu \cdot (T-\tau)} \quad \text{with} \quad \|\hat{\sigma}_D^*\| = \min\{R^*, R_1\} $$

using the notations:

$$ \bullet \quad \left\{ \begin{array}{l} \Delta \sigma = \hat{\sigma}_h - \hat{\sigma}_h \\ \Delta \dot{\hat{\epsilon}}_p = \dot{\hat{\epsilon}}_{p,h} - \dot{\hat{\epsilon}}_{p,h} \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{l} \Delta R^* = R^* - R_h \\ \Delta \dot{p}^* = \dot{p}^* - \dot{p}_h \end{array} \right. $$

$$ \bullet \quad \forall \mathbf{T} \in S^{(0:T)} \|\mathbf{T}\| \neq 0, \quad P_{\mathbf{T}}(\square) = \frac{\operatorname{Tr}(\square, \mathbf{T})}{\sqrt{\operatorname{Tr}(\square, \mathbf{T})}}, \quad \text{and} \quad \mathbf{P}_{\mathbf{T}}^{\perp}(\square) = \square - P_{\mathbf{T}}(\square), \quad \mathbf{T}^{\perp} = \mathbf{T} $$

$$ \bullet \quad R_1 = \|\hat{\sigma}_{D,h}\| + 2\mu \cdot (T - \tau) \cdot \dot{p}_h + P_{\sigma_{D,h}}(\Delta \dot{\hat{\epsilon}}) - P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) \cdot \blacksquare $$

**Proof.**

$$ \begin{aligned} l(\hat{p}^*, R^*, \|\hat{\sigma}_D^*\|) ={}& R^* \dot{p}^* + \hat{\sigma}_h : \dot{\hat{\epsilon}}_{p,h} - \hat{\sigma}^* : \dot{\hat{\epsilon}}_{p,h} - \hat{\sigma}_h : \dot{\hat{\epsilon}}_p^* - \delta \hat{\sigma}_h^* \\ & : \dot{\hat{\epsilon}}_p^* + \sigma^* : \delta \hat{\epsilon}_{p,h} + \delta \hat{\sigma}_h : \dot{\hat{\epsilon}}_{p,h} - \hat{\sigma}_h^* \\ & : \delta \hat{\epsilon}_{p,h} - \frac{\hat{\sigma}_D^* : \hat{\sigma}_{D,h}}{2\mu \cdot (T-\tau)} + \frac{\|\hat{\sigma}_D^*\|^2 + \|\hat{\sigma}_{D,h}\|^2}{4\mu \cdot (T-\tau)}. \end{aligned} \quad (27) $$

By minimizing the terms, which involves a tensor direction, one gets:

$$ \begin{aligned} l(\dot{p}^*, R^*, \|\hat{\sigma}_D^*\|) ={}& R^* \dot{p}^* + \hat{\sigma}_h : \dot{\hat{\epsilon}}_{p,h} \\ & - \|\hat{\sigma}_D^*\| \left( \dot{\hat{\epsilon}}_{p,h} + \frac{\hat{\sigma}_{D,h}}{2\mu \cdot (T-\tau)} - \delta \hat{\epsilon}_{p,h} \right) + \dot{p}^* (\hat{\sigma}_{D,h} + \delta \hat{\sigma}_{D,h}) \\ & + \delta \hat{\sigma}_h : \dot{\hat{\epsilon}}_{p,h} - \hat{\sigma}_h : \delta \hat{\epsilon}_{p,h} \\ & + \frac{\|\hat{\sigma}_D^*\|^2 + \|\hat{\sigma}_{D,h}\|^2}{4\mu \cdot (T-\tau)}, \end{aligned} \quad (28) $$

$$ \begin{aligned} ={}& l_1(\hat{p}^*, R^*, \|\hat{\sigma}_D^*\|) + R^* \dot{p}^* + \hat{\sigma}_h : \dot{\hat{\epsilon}}_{p,h} + \delta \hat{\sigma}_h : \dot{\hat{\epsilon}}_{p,h} - \hat{\sigma}_h : \delta \hat{\epsilon}_{p,h} \\ & + \frac{\|\hat{\sigma}_D^*\|^2 + \|\hat{\sigma}_{D,h}\|^2}{4\mu \cdot (T-\tau)} - \|\hat{\sigma}_D^*\| \dot{p}_h + P_{\sigma_{D,h}}(\Delta \dot{\hat{\epsilon}}_p) + \frac{P_{\sigma_{D,h}}(\hat{\sigma}_{D,h})}{2\mu \cdot (T-\tau)} - P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) \\ & - \dot{p}^* (\|\hat{\sigma}_{D,h}\| + P_{\sigma_{D,h}}(\Delta \sigma) + P_{\sigma_{D,h}}(\delta \hat{\sigma}_{D,h})). \end{aligned} \quad (29) $$

Let us note that $\min_{\|\hat{\sigma}_D^*\| \le R^*} \{l_1(\hat{p}^*, R^*, \|\hat{\sigma}_D^*\|)\}$ is equal to:

$$ l_1(\hat{p}^*, R^*) = l_1(\hat{p}^*, R^*, R^*), $$

$$ \begin{aligned} l(\hat{p}^*, R^*, \|\hat{\sigma}_D^*\|) ={}& l_1(\hat{p}^*, R^*) + R^* \dot{p}^* + \hat{\sigma}_h : \dot{\hat{\epsilon}}_{p,h} + \delta \hat{\sigma}_h : \dot{\hat{\epsilon}}_{p,h} - \hat{\sigma}_h : \delta \hat{\epsilon}_{p,h} \\ & \blacksquare : \delta \hat{\epsilon}_{p,h} - R^* \left( \dot{p}_h + P_{\sigma_{D,h}}(\Delta \dot{\hat{\epsilon}}_p) - P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) \right) \\ & - \dot{p}^* (\|\hat{\sigma}_{D,h}\| + P_{\sigma_{D,h}}(\Delta \sigma) + P_{\sigma_{D,h}}(\delta \hat{\sigma}_{D,h}) \\ & + l_3(R^*, \|\hat{\sigma}_D^*\|) \end{aligned} \quad (30) $$

$$ \text{with } l_3(R^*, \|\hat{\sigma}_D^*\|) = (R^* - \|\hat{\sigma}_D^*\|) \dot{p}_h + P_{\sigma_{D,h}}(\Delta \dot{\hat{\epsilon}}_p) - P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) \\ + \frac{(\|\hat{\sigma}_D^*\| - \|\hat{\sigma}_{D,h}\|)^2}{4\mu \cdot (T-\tau)} + \frac{\|\hat{\sigma}_D^*\| (\|\hat{\sigma}_h\| - P_{\sigma_{D,h}}(\hat{\sigma}_h))}{2\mu \cdot (T-\tau)}. $$

Consequently,

$$ l_3(R^*) = \min_{\|\hat{\sigma}_D^*\| \le R^*} \{l_3(R^*, \|\hat{\sigma}_D^*\|)\} $$

can be obtained for

$$ \|\hat{\sigma}_D^*\| = R^* $$

$$ \text{or } \|\hat{\sigma}_D^*\| = \|\hat{\sigma}_{D,h}\| + 2\mu \cdot (T - \tau) \cdot \dot{p}_h + P_{\sigma_{D,h}}(\Delta \dot{\hat{\epsilon}}_p) - P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) = R_1. $$

Then,

$$ l_3(R^*) = \frac{(R^* - R_1)(R_1 - \|\hat{\sigma}_{D,h}\|)}{2\mu \cdot (T - \tau)} + \frac{(R_1 - \|\hat{\sigma}_{D,h}\|)^2}{4\mu \cdot (T - \tau)} \quad \text{if } R_1 \le R^* \\ \text{and } l_3(R^*) = \frac{(R^* - \|\hat{\sigma}_{D,h}\|)^2}{4\mu \cdot (T - \tau)} \quad \text{if } R_1 > R^*. \quad (31) $$

In order to get $l_2(\hat{p}^*, R^*)$, one adds the quantity:

$$ z = \int_0^T \int_0^T \left\{ \delta \hat{\underline{R}}_h \Delta \dot{p}^* - \delta \hat{\underline{p}}_h \Delta R^* \right\} \frac{d\tau dt}{T} \quad (32) $$

$$ \text{with } \Delta \dot{p}^* = \dot{p}^* - \dot{p}_h \quad \delta R^* = R^* - R_h. $$

This quantity is equal to zero because $\delta \hat{\underline{R}}_h|_T = 0$ and $\Delta R^*_0 = 0$:

$$ z = \int_0^T \left\{ \left( \frac{\delta \hat{\underline{R}}_h \Delta R^*}{k} \right)_t \right\} dt = 0. \quad \square \quad (33) $$

### 5.2. Analysis of the order of magnitude of functional $\ell$

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 $\Delta$ terms must be small to get sharp bounds, which can always be achieved by an appropriate choice of the sizing parameter $\lambda$;
<figure><img src="image_8.png" /><figcaption>Fig. 6. The geometric approach used for the resolution of the central problem.</figcaption></figure>

<table><caption>**Table 1**
The values of the lower and upper bounds as functions of parameter $\lambda$ ($\lambda > 0$ or $\lambda < 0$).</caption><thead><tr><th>λ</th><th>max (<sup>Δσ<sub>u</sub></sup>&frasl;<sub>σ<sub>h,h</sub></sub>), max (<sup>Δε<sub>u,u</sub></sup>&frasl;<sub>ε<sub>h,h</sub></sub>))</th><th>Lower bound of -η (×10<sup>-3</sup>)</th><th>Upper bound of η (×10<sup>-3</sup>)</th><th>α<sub>h</sub> + α<sub>hh</sub> (×10<sup>-3</sup>)</th><th>Δα<sup>η</sup><sub>h</sub> (%)</th></tr></thead><tbody><tr><td>10<sup>2</sup></td><td>(1.9610<sup>-1</sup>, 1.31)</td><td>(-1.9710<sup>-1</sup>, -1.66)</td><td>(2.1810<sup>-2</sup>, 3.78)</td><td>5.25</td><td>51.80</td></tr><tr><td>10<sup>1</sup></td><td>(1.9610<sup>-2</sup>, 1.3110<sup>-1</sup>)</td><td>(-2.18, -1.46)</td><td>(6.5910<sup>-1</sup>, 3.53)</td><td>5.22</td><td>47.73</td></tr><tr><td>10<sup>0</sup></td><td>(1.9610<sup>-3</sup>, 1.3110<sup>-2</sup>)</td><td>(-2.2010<sup>1</sup>, -1.44)</td><td>(6.64, 3.51)</td><td>5.22</td><td>47.41</td></tr><tr><td>10<sup>-1</sup></td><td>(1.9610<sup>-4</sup>, 1.3110<sup>-3</sup>)</td><td>(-2.2010<sup>2</sup>, -1.44)</td><td>(6.6510<sup>1</sup>, 3.50)</td><td>5.22</td><td>47.34</td></tr><tr><td>10<sup>-2</sup></td><td>(1.9610<sup>-5</sup>, 1.3110<sup>-4</sup>)</td><td>(-2.2010<sup>3</sup>, -1.44)</td><td>(6.6610<sup>2</sup>, 3.50)</td><td>5.21</td><td>47.33</td></tr><tr><td>10<sup>-3</sup></td><td>(1.9610<sup>-6</sup>, 1.3110<sup>-5</sup>)</td><td>(-2.2010<sup>4</sup>, -1.44)</td><td>(6.6610<sup>3</sup>, 3.50)</td><td>5.21</td><td>47.33</td></tr><tr><td>-10<sup>2</sup></td><td>(1.9610<sup>-1</sup>, 1.31)</td><td>(6.3010<sup>-1</sup>, -4.15)</td><td>(-5.0410<sup>-2</sup>, 6.09)</td><td>5.16</td><td>99.22</td></tr><tr><td>-10<sup>0</sup></td><td>(1.9610<sup>-3</sup>, 1.3110<sup>-2</sup>)</td><td>(6.3010<sup>1</sup>, -4.15)</td><td>(-5.04, 6.09)</td><td>5.16</td><td>99.22</td></tr><tr><td>-10<sup>-2</sup></td><td>(1.9610<sup>-5</sup>, 1.3110<sup>-4</sup>)</td><td>(6.3010<sup>3</sup>, -4.15)</td><td>(-5.0410<sup>2</sup>, 6.09)</td><td>5.16</td><td>99.22</td></tr></tbody></table>

• many terms of the functional can be neglected.

### 5.3. Resolution of the central problem

Using the solutions $(\hat{\sigma}_h, \dot{\hat{\epsilon}}_h)$, $(\sigma_h, \dot{\epsilon}_h)$ and $(\delta\hat{\sigma}_h, \delta\dot{\epsilon}_h)$, the problem can be written as a minimization problem over $Q^{[0:T]}$:

$$-r = \min_{(\hat{p}^*, R^*) \in Q^{[0:T]}} \left\{ \int_0^T \int_0^t \{ \Delta R^* \Delta \hat{p}^* + \Delta \hat{p}^* \cdot \mathcal{I} + \Delta R^* \cdot \mathcal{J} + \mathcal{K} \} \frac{d\tau dt}{T} \right\}.$$

Coefficients $\mathcal{I}$, $\mathcal{J}$ and $\mathcal{K}$ are known and are small, which leads to the following local equations:

$$\left\{
\begin{array}{l}
\displaystyle \int_t^T \{ \Delta R^* \frac{d\tau}{T} \} + Q \ge 0, \\
\displaystyle \dot{p}^* \ge 0, \\
\displaystyle \left( \int_t^T \{ \Delta R^* \} \frac{d\tau}{T} + Q \right) \cdot \dot{p}^* = 0, \\
\displaystyle R^* = R_0 + k \cdot p^*, \\
\displaystyle p^*(t=0) = 0,
\end{array}
\right.
\quad (34)$$

where $Q = \frac{T-t}{T} \mathcal{I} + \int_t^T \{ \frac{T-t}{T} k \cdot \mathcal{J} \} d\tau$.

In order to solve this problem inexpensively, let us build the solution geometrically. First, let us define $\bar{Z}_{|t}$ as:

$$\bar{Z}_{|t} = \int_t^T \{k \cdot p_h\} d\tau - Q_{|t}.$$

Function $\int_t^T \{k \cdot p^*\}$ is the convex hull, which is such that:

$$R^* - R_0 = 0 \quad \text{at} \quad t = T,$$

$$\frac{d}{dt}[R^* - R_0] = 0 \quad \text{at} \quad t = 0. \quad (35)$$

Obviously, the corresponding curve closely matches that which corresponds to $\int_t^T \{k \cdot p^*\} d\tau$. This approach is illustrated in Fig. 6.

**Remark 3.** If the discretization of $Q$ is piecewise linear within $[0; T] = \cup_{1 \le i \le n} \{[t_{i-1}, t_i]\}$, the $\dot{p}^*$ solution is a combination of Dirac distributions corresponding to the times $t_i$ ($1 \le i \le n$) of the discretization.

**Remark 4.** It is easy to take the complete expression $l(\dot{p}^*, R^*)$ into account. In this case, $Q$ depends on $(\dot{p}^*, R^*)$ and one has

$$Q(0, R^*) \le Q(+\infty, R^*) < \infty.$$

**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(\dot{p}_1^*, R_1^*)$, and iterates if necessary.

## 6. 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

$$\alpha_h = 4.19 \ 10^{-3}.$$

Using the notations introduced in (12), one can define an upper bound of the local error denoted $\Delta\alpha_h$:
<figure><img src="image_9.png" /><figcaption>Fig. 7. The quadratic behavior of the upper and lower bounds; the calculated values are represented by diamonds and the solid lines are quadratic interpolations of $\epsilon_h(\eta, \lambda = 1)$.</figcaption></figure>

$$|\alpha_h + \alpha_{hh} - \alpha_{ex}| \le \Delta\alpha_h.$$

This bounding can be applied directly to the exact value of the output of interest, leading to its extremum:

$$\alpha_{ex} = (\alpha_h + \alpha_{hh}) + \Delta\alpha_h.$$

Let us introduce the nondimensional error measure:

$$\Delta\alpha_h^{cr} = \frac{\Delta\alpha_h}{\alpha_h + \alpha_{hh}}.$$

Table 1 shows the evolution of the bounds for different values of parameters $(\eta, \lambda)$. First, one can observe the significant role played by the sign of parameter $\lambda$. 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:

$$-0.26 \le \sigma_{h,xx} \le 0.84; \quad 1.10 \le \sigma_{h,yy} \le 1.74 \quad \text{and} \quad -0.06 \le \sigma_{h,xy} \le 0.39.$$

The prestress extractor in the zone of interest takes the values:

$$\begin{cases} \lambda = 1: & \delta\tilde{\sigma}_{\Sigma,h,xx} = 0, \quad \delta\tilde{\sigma}_{\Sigma,h,yy} = 0.25 \quad \text{and} \quad \delta\tilde{\sigma}_{\Sigma,h,xy} = 0, \\ \lambda = -1: & \delta\tilde{\sigma}_{\Sigma,h,xx} = 0, \quad \delta\tilde{\sigma}_{\Sigma,h,yy} = -0.25 \quad \text{and} \quad \delta\tilde{\sigma}_{\Sigma,h,xy} = 0. \end{cases}$$

**Remark 6.** On this example (Table 1), for negative values of parameter $\lambda$, the behavior of the mirror problem is elastic and corresponding bounds are not sharp. For positive values of parameter $\lambda$, the best values of the parameter are chosen to optimize the bounds.

**Remark 7.** Parameters $\eta$ and $\eta'$ were chosen by means of a quadratic fitting method. Near the optimum values, $\epsilon_{h,\lambda \text{cst}}$ has quadratic behavior, as illustrated in Fig. 7. In this figure, the lower and upper bounds were plotted using a constant extraction parameter $\lambda$ equal to 1. The calculated points are represented by diamonds and a quadratic interpolation of $\epsilon_h(\eta, \lambda = 1)$ is also shown. 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 $\alpha^-$ or $\alpha^+$.

**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 been done for viscoelasticity [27,17,18,28] and for dynamics [20,29]. In the linear case, it was shown that the more

<figure><img src="image_10.png" /><figcaption>Fig. 8. The parameters of $a(t)$ and the values chosen for each bound.</figcaption></figure>

accurate the resolution of the adjoint problem, the tighter the bounds. We showed that the same property also holds in the nonlinear case.

## 7. 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

$$\lambda(\alpha_{ex} - \alpha_h) \le \frac{1}{\eta} \int_{\Omega} \left\{ \max_{(\dot{\epsilon}_p, \dot{\sigma}) \in I^{0,T}} \left\{ \int_0^T \left\{ \eta \cdot \text{Tr} \left[ (\dot{\epsilon}_p^* - \dot{\epsilon}_{p,h}) \delta\hat{\sigma}_h - (\sigma^* - \hat{\sigma}_h) \delta\hat{\epsilon}_{p,h} \right] \right\} \right\} dt - \int_0^T \left\{ a(t) \left( \int_0^t \left\{ \text{Tr}[(\sigma^* - \sigma_h^*)(\dot{\epsilon}_p^* - \dot{\epsilon}_{p,h})] \right\} \frac{d\tau}{T} \right. \right. \\ \left. \left. + \frac{1}{2T} \text{Tr}[(\sigma^* - \sigma_h)\mathbf{K}^{-1}(\sigma^* - \sigma_h)] \right) \right\} dt \right\} dV. \quad (36)$$

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:

$$\alpha_h + \alpha_{hh} = 4.07 \times 10^{-3} \quad \text{and} \quad \Delta\alpha_h^{cr} = 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 extended 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.

## Appendix A. Analysis of the functional of the central problem

This appendix focuses on some technical points relative to the central problem. The analysis relies on the decomposition of functional $ℓ$ into several parts $ℓ_1, ℓ_2$ and $ℓ_3$. These terms are then studied separately. The study is done in the very interesting case where the $δ$ terms $Δ$ are small.

### A.1. Analysis of term $l_1(\dot{p}^*, R^*)$

$$ \ell_1(\dot{p}^*, R^*) = - \left\| \begin{aligned} & R^* \cdot \dot{\epsilon}_{p,h} + \frac{\hat{\sigma}_{D,h}}{2\mu \cdot (T-\tau)} + \Delta \dot{\epsilon}_p - \delta \hat{\epsilon}_{p,h} \\ & + \dot{p}^* \cdot (\sigma_{D,h} + \Delta \sigma_D + \delta \hat{\sigma}_{D,h}) \\ & + R^* \cdot \dot{p}_h + \frac{P_{\sigma_{D,h}}(\hat{\sigma}_{D,h})}{2\mu \cdot (T-\tau)} + P_{\sigma_{D,h}}(\Delta \dot{\epsilon}_p) - P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) \\ & + \dot{p}^* \cdot (\|\sigma_{D,h}\| + P_{\sigma_{D,h}}(\Delta \sigma) + P_{\sigma_{D,h}}(\delta \hat{\sigma}_h)) \end{aligned} \right\|^2 $$

$ℓ_1(\dot{p}^*, R^*)$ can also be written as

$$ \ell_1(\dot{p}^*, R^*) = -\frac{N}{D}, \quad (A.1) $$

where

$$ \left\{ \begin{aligned} N &= \left\| R^* \left( \dot{\epsilon}_{p,h} + \frac{\sigma_{D,h}}{2\mu \cdot (T-\tau)} + Q \right) + \dot{p}^* \cdot (\sigma_{D,h} + R) \right\|^2 \\ &\quad - \left( R^* \left( \dot{p}_h + \frac{\|\sigma_{D,h}\|}{2\mu \cdot (T-\tau)} + P_{\sigma_{D,h}}(Q) \right) + \dot{p}^* \cdot (\|\sigma_{D,h}\| + R) \right)^2, \\ D &= \left\| R^* \left( \dot{\epsilon}_{p,h} + \frac{\sigma_{D,h}}{2\mu \cdot (T-\tau)} + Q \right) + \dot{p}^* \cdot (\sigma_{D,h} + R) \right\|^2 \\ &\quad + \left( R^* \left( \dot{p}_h + \frac{\|\sigma_{D,h}\|}{2\mu \cdot (T-\tau)} + P_{\sigma_{D,h}}(Q) \right) + \dot{p}^* \cdot (\|\sigma_{D,h}\| + R) \right). \end{aligned} \right. $$

The numerator $N$ is equal to

$$ N = (R^* \cdot \mathbb{P}_{\sigma_{D,h}}^{\perp}(Q) + \dot{p}^* \cdot \mathbb{P}_{\sigma_{D,h}}^{\perp}(R))^2, $$

which is equivalent to

$$ N = R^* \cdot \mathbb{P}_{\sigma_{D,h}}^{\perp}(\Delta \dot{\epsilon}_p) + \frac{\mathbb{P}_{\sigma_{D,h}}^{\perp}(\Delta \sigma_D)}{2\mu \cdot (T-\tau)} + \dot{p}^* \mathbb{P}_{\sigma_{D,h}}^{\perp}(\Delta \sigma_D) \\ + \left[ \dot{p}^* \mathbb{P}_{\sigma_{D,h}}^{\perp}(\delta \hat{\sigma}_h) - R^* \cdot \mathbb{P}_{\sigma_{D,h}}^{\perp}(\delta \hat{\epsilon}_{p,h}) \right]^2 \quad (A.2) $$

Using (A.15), one gets:

$$ \begin{aligned} & \dot{p}^* \cdot \mathbb{P}_{\sigma_{D,h}}^{\perp}(\delta \hat{\sigma}_h) - R^* \cdot \mathbb{P}_{\sigma_{D,h}}^{\perp}(\delta \hat{\epsilon}_{p,h}) \\ & \sim \Delta \dot{p}^* \cdot \mathbb{P}_{\sigma_{D,h}}^{\perp}(\delta \hat{\sigma}_h) - \Delta R^* \cdot \mathbb{P}_{\sigma_{D,h}}^{\perp}(\delta \hat{\epsilon}_{p,h}). \end{aligned} \quad (A.3) $$

Therefore, $ℓ_1(\dot{p}^*, R^*)$ is very small. The FE solution $ℓ_1(\dot{p}_h, R_h)$ gives a very good approximation of its value. This is what we use in practice, even though it would not be difficult to use the exact term $ℓ_1(\dot{p}^*, R^*)$.

### A.2. Analysis of the central problem when the $δ$ terms equal zero

The central problem consists in minimizing $-r(\dot{p}^*, R^*)$ over $Q^{[0:T]}$ with

$$ \begin{aligned} \ell(\dot{p}^*, R^*) = & \ell_1(\dot{p}^*, R^*) + \ell_3(R^*) + R^* \cdot \dot{p}^* + \hat{\sigma}_h : \dot{\epsilon}_{p,h} \\ & - R^* \left( \dot{p}_h + P_{\sigma_{D,h}}(\Delta \dot{\epsilon}_p) - (\|\sigma_D\| + P_{\sigma_{D,h}}(\Delta \sigma)) \cdot \dot{p}^* \right) = \ell_1(\dot{p}^*, R^*) \\ & + \ell_3(R^*) + \Delta \sigma_D : \Delta \dot{\epsilon}_p - (R_h - \|\sigma_{D,h}\|) : \dot{\epsilon}_{p,h} + \Delta R^* \cdot \Delta \dot{p}^* \\ & + \underbrace{(R_h - \|\sigma_D, h\|) : \dot{p}^* - \Delta \dot{p}^* P_{\sigma_{D,h}}(\Delta \sigma_D) - \Delta R^* P_{\sigma_{D,h}}(\Delta \dot{\epsilon}_p)}_{\ge 0} \end{aligned} \quad (A.4) $$

Term $ℓ_3$, which represents an elastic energy, is small. Over $[0; T]$, there are two main possibilities:

(i) $\dot{p}^* > 0$ and $R^* = \|\sigma_D^*\|$; then, $ℓ_3$ is positive.

(ii) $\dot{p}^* = 0$ and $R_1 = \|\sigma_{D,h}\| + O(\epsilon(\|\sigma_{D,h}\|))$. (A.5)

In the latter case of (A.5), $ℓ_3$ can be written as:

$$ \begin{aligned} \ell_3 = & (\Delta R^* + (R_h - \|\sigma_D\|)) \cdot \frac{R_1 - \|\sigma_{D,h}\|}{2\mu \cdot (T-\tau)} + \frac{(\|\sigma_D\| - \|\sigma_{D,h}\|)^2}{4\mu \cdot (T-\tau)} \\ = & \Delta R^* \frac{\|\sigma_D\| - \|\hat{\sigma}_{D,h}\|}{2\mu \cdot (T-\tau)} + \frac{(R_h - \|\sigma_D\|)(R_1 - \|\hat{\sigma}_{D,h}\|)}{2\mu \cdot (T-\tau)} \\ & + \frac{(\|\sigma_D\| - \|\hat{\sigma}_{D,h}\|)^2}{4\mu \cdot (T-\tau)}. \end{aligned} \quad (A.6) $$

The positive part of $ℓ(\dot{p}^*, R^*)$ corresponds to a negative contribution to the error. Then, $ℓ$ leads to

$$ \begin{aligned} \Delta \dot{p} &= O(\epsilon), \\ \Delta R &= O(\epsilon) \end{aligned} \quad (A.7) $$

and, therefore, to a value of $\bar{r}$ which is $O(\epsilon^2)$.

### A.3. Analysis of the central problem when the $Δ$ terms equal zero

In this case, the central problem consists in minimizing $-\bar{r}(\dot{p}^*, R^*)$ over $Q^{[0:T]}$ with

$$ \begin{aligned} \ell(\dot{p}^*, R^*) = & \ell_1(\dot{p}^*, R^*) + \ell_3(R^*) + R^* \cdot \dot{p}^* + \|\sigma_{D,h}\| \dot{p}_h - \|\sigma_{D,h}\| \\ & - R^* \cdot \dot{p}_h - \dot{p}^* \cdot P_{\sigma_{D,h}}(\delta \hat{\sigma}_{D,h}) + R^* \cdot P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) \\ & + \dot{p}^* \cdot P_{\sigma_{D,h}}(\delta \hat{\sigma}_{D,h}) - \|\sigma_{D,h}\| P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) + \Delta \dot{p}^* \delta \hat{R} \\ & - \Delta R^* \delta \hat{p}_h = \ell_1(\dot{p}^*, R^*) + \ell_3(R^*) + \Delta R^* \Delta \dot{p}^* + (R_h - \|\sigma_{D,h}\|) \dot{p}^* \\ & + \Delta \dot{p}^* (\delta \hat{R}_h - P_{\sigma_{D,h}}(\delta \hat{\sigma}_{D,h})) + \Delta R^* (P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) - \delta \hat{p}_h) \\ & + (R_h - \|\sigma_{D,h}\|) P_{\sigma_{D,h}}(\delta \hat{\epsilon}_p). \end{aligned} \quad (A.8) $$

The analysis of $l_3$ is similar to that of Paragraph A.2. One can also note that $P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) \sim \delta \hat{p}_h > 0$ and $P_{\sigma_{D,h}}(\delta \hat{\epsilon}_{p,h}) - \delta \hat{p}_h = O(\epsilon')$. Let us analyze the term

$$ \omega = (R_h - \|\sigma_{D,h}\|) \dot{p}^* + \dot{p}^* (\delta \hat{R}_h - P_{\sigma_{D,h}}(\delta \hat{\sigma}_h)) - \dot{p}_h (\delta \hat{R}_h - P_{\sigma_{D,h}}(\delta \hat{\sigma}_{D,h})) $$

for the two main patterns

$$ \bullet (R_h - \|\sigma_{D,h}\|) > 0 \text{ and } \dot{p}_h = 0 $$

$$ \omega = (R_h - \|\sigma_{D,h}\| + \delta \hat{R}_h - P_{\sigma_{D,h}}(\delta \hat{\sigma}_{D,h})) \dot{p}^* \quad (A.9) $$

Since the perturbation $\delta \hat{R}_h - P_{\sigma_{D,h}}(\delta \hat{\sigma}_{D,h})$ is small, $\omega$ is positive over $[0; T]$.

Therefore, in this case, $r = O(\epsilon^2)$.

$$ \bullet (R_h - \|\sigma_{D,h}\|) = 0 $$
$$ \omega = \frac{(R_h - \|\boldsymbol{\sigma}_{D,h}\| + \delta\hat{R}_h - P_{\sigma_{D,h}}(\delta\hat{\boldsymbol{\sigma}}_{D,h}))\dot{p}^* - \dot{p}_h(\delta\hat{R}_h}{\ge 0} - P_{\sigma_{D,h}}(\delta\hat{\boldsymbol{\sigma}}_{D,h})). \qquad (A.10) $$

Consequently, $(\Delta\dot{p}^*, \Delta R^*) \sim O(\epsilon')$. Then, $\omega = \Delta\dot{p}^*(\delta\hat{R}_h - P_{\sigma_{D,h}}(\delta\hat{\boldsymbol{\sigma}}_{D,h}))$. This term is still small and is $O(\epsilon'(P_{\sigma_{D,h}}(\delta\hat{\boldsymbol{\sigma}}_{D,h})^2))$.

### A.4. Orders of magnitude of the $\Delta$ and $\delta$ terms

**$\Delta$ terms** Since they represent an error from the reference problem, the $\Delta$ terms are intrinsically very small quantities. One can also write:

$$ \begin{aligned} \Delta\boldsymbol{\sigma} &\sim O(\epsilon) \cdot \|\boldsymbol{\sigma}_{D,h}\|, \\ \Delta\dot{\epsilon}_p &\sim O(\epsilon) \cdot \frac{\|\boldsymbol{\sigma}_{D,h}\|}{k \cdot T}. \end{aligned} \qquad (A.11) $$

If the behavior is elastic (i.e. $\|\boldsymbol{\sigma}_{D,h}\| - R_h < 0$), the second expression becomes

$$ \Delta\dot{\epsilon}_p \sim O(\epsilon) \cdot \frac{\|\boldsymbol{\sigma}_{D,h}\|}{2\mu T}. $$

**$\delta$ terms** The $\delta$ terms require a more detailed analysis. By construction, the mirror solution based on the two extractors is defined as a perturbation. Parameter $\lambda$ is what enables this property to be verified. In practice, one chooses $\lambda$ such that each component of $\delta\hat{\boldsymbol{\sigma}}_{\Sigma}$ (respectively $\delta\bar{\boldsymbol{\epsilon}}_{\Sigma}$) is about $10^2$ times smaller than $\boldsymbol{\sigma}_h$ (respectively $\dot{\epsilon}_h$). Therefore, $\delta\boldsymbol{\sigma}$ is relatively small:

$$ \delta\boldsymbol{\sigma} = O(\epsilon') = O(\sigma_h). $$

Furthermore, the $\delta$ terms inherently satisfy the constitutive relation $\dot{B}$. Therefore, one has:

$$ \delta \boldsymbol{\epsilon}_p = (\dot{p}_h + \delta \dot{p}) \cdot \frac{\boldsymbol{\sigma}_{D,h} + \delta \boldsymbol{\sigma}_D}{\|\boldsymbol{\sigma}_{D,h} + \delta \boldsymbol{\sigma}_D\|} - \dot{\epsilon}_{p,h} \qquad (A.12) $$

$$ \begin{aligned} & \text{with } \delta \dot{p} (\|\boldsymbol{\sigma}_{D,h} + \delta \boldsymbol{\sigma}_D\| - (R_h + \delta R)) = 0 \\ & \delta R_{|\tau} = \sup_{\tau' \le \tau} \{ \|\boldsymbol{\sigma}_{D,h|T-\tau'} + \delta \boldsymbol{\sigma}_{D|\tau'}\| - R_{h|T-\tau'}, 0 \}. \end{aligned} $$

Consequently:

$$ \delta \boldsymbol{\epsilon}_p = \delta \dot{p} \frac{\boldsymbol{\sigma}_{D,h} + \delta \boldsymbol{\sigma}_D}{\|\boldsymbol{\sigma}_{D,h} + \delta \boldsymbol{\sigma}_D\|} + \dot{p}_h \left( \frac{\boldsymbol{\sigma}_{D,h} + \delta \boldsymbol{\sigma}_D}{\|\boldsymbol{\sigma}_{D,h} + \delta \boldsymbol{\sigma}_D\|} - \frac{\boldsymbol{\sigma}_{D,h}}{\|\boldsymbol{\sigma}_{D,h}\|} \right), \qquad (A.13) $$

which leads to:

$$ P_{\sigma_{D,h}}(\delta \boldsymbol{\epsilon}_p) = \delta \dot{p} + O(\epsilon'), \qquad (A.14) $$

$$ \mathbb{P}^{\perp}_{\sigma_{D,h}}(\delta \boldsymbol{\epsilon}_p) = \frac{\dot{p}_h}{\|\boldsymbol{\sigma}_{D,h}\|} \mathbb{P}^{\perp}_{\sigma_{D,h}}(\delta \boldsymbol{\sigma}_D) + O(\epsilon'), \qquad (A.15) $$

$$ \text{and } \delta R_{|\tau} = \sup_{\tau' \le \tau} \{ \|\boldsymbol{\sigma}_{D,h|T-\tau'}\| - R_{h|T-\tau'} + P_{\sigma_{D,h}}(\delta\boldsymbol{\sigma}_{D|\tau'}) + O(\epsilon'), 0 \}. \qquad (A.16) $$

Moreover, with regard to $(\dot{p}_h, R_h)$, one can identify two main patterns over $[0; T]$:

$$ \bullet R_h - \|\boldsymbol{\sigma}_{D,h}\| \gg O(\epsilon \cdot \|\boldsymbol{\sigma}_{D,h}\|, \|\delta\boldsymbol{\sigma}_D\|) \quad \text{and} \quad R_h - \|\boldsymbol{\sigma}_{D,h}\| > 0; \quad \dot{p}_h = 0, $$

$$ \bullet R_h - \|\boldsymbol{\sigma}_{D,h}\| = 0; \quad \dot{p}_h \gg O\left(\frac{\epsilon \cdot \|\boldsymbol{\sigma}_{D,h}\|}{T \cdot k}, \|\delta \dot{\boldsymbol{\epsilon}}_p\|\right). \qquad (A.17) $$

The portion of $[0; T]$ where the FE solution is not described by either of the two patterns is very small. Then, one has:

$$ \begin{aligned} \bullet \delta \boldsymbol{\epsilon}_p &= 0 \quad \text{for } R_h - \|\boldsymbol{\sigma}_{D,h}\| > 0, \\ \bullet P_{\sigma_{D,h}}(\delta \boldsymbol{\epsilon}_p) - \delta \dot{p} &= O\left(\epsilon' \frac{\delta \boldsymbol{\sigma}_D}{2\mu}\right) \quad \text{for } R_h - \|\boldsymbol{\sigma}_{D,h}\| = 0. \end{aligned} \qquad (A.18) $$

For the second pattern, the plastic strain can be expressed as $\delta p = \sup_{\tau' \le \tau} \{ \|P_{\sigma_{D,h}}(\delta\boldsymbol{\sigma}_{D|\tau'}) + O(\epsilon')\|, 0 \}$, which leads to small quantities.

In the general case, the $\Delta$ 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 $\lambda$ enables this property to be always satisfied.

### References

[1] P. Ladevèze, J. Oden (Eds.), Advances in Adaptive Computational Methods in Mechanics, Elsevier, 1998.

[2] I. Babuška, T. Strouboulis, The Finite Element Method and its Reliability, Oxford Science Publications, 2001.

[3] E. Stein, E. Ramm, E. Rank, R. Rannacher, K. Schweizerhof, E. Stein, W. Wendland, G. Wittum, P. Wriggers, W. Wunderlich, Error-controlled Adaptive Finite Elements in Solid Mechanics, Wiley, 2003.

[4] P. Ladevèze, J. Pelle, Mastering Calculations in Linear and Nonlinear Mechanics, Springer, NY, 2004.

[5] N. Wiberg, P. Diez, Comput. Methods Appl. Mech. Engrg. Special issue, 2006.

[6] K. Eriksson, D. Estep, P. Hansbo, C. Johnson, Introduction to adaptive methods for differential equations, Acta Numer. 4 (1995) 105-158.

[7] R. Becker, R. Rannacher, A feed-back approach to error control in finite element methods: basic analysis and examples, East-West J. Numer. Math. 4 (1996) 237-264.

[8] R. Rannacher, F.T. Suttmeier, A posteriori error control in finite element methods via duality techniques: application to perfect plasticity, Comput. Mech. 21 (1998) 123-133.

[9] F. Cirak, E. Ramm, A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem, Comput. Methods Appl. Mech. Engrg. 156 (1998) 351-362.

[10] S. Prudhomme, J. Oden, On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors, Comput. Methods Appl. Mech. Engrg. 176 (1999) 313-331.

[11] P. Ladevèze, P. Rougeot, P. Blanchard, J. Moreau, Local error estimators for finite element linear analysis, Comput. Methods Appl. Mech. Engrg. 176 (1999) 231-246.

[12] T. Strouboulis, I. Babuška, D.K. Datta, K. Copps, S.K. Gangaraj, A posteriori estimation and adaptive control of the error in the quantity of interest. part i: A posteriori estimation of the error in the von mises stress and the stress intensity factor, Comput. Methods Appl. Mech. Engrg. (2000) 261-294.

[13] N. Parés, J. Bonet, A. Huerta, J. Peraire, The computation of bounds for linear-functional outputs of weak solutions to the two-dimensional elasticity equations, Comput. Methods Appl. Mech. Engrg. (2006) 406-429.

[14] J. Peraire, A. Patera, Bounds for linear-functional outputs of coercive partial differential equations: local indicators and adaptive refinement, in: P. Ladevèze, J. Oden (Eds.), Advances in Adaptive Computational Methods in Mechanics, Studies in Applied Mechanics, vol. 47, Elsevier, 1998, pp. 199-216.

[15] P. Ladevèze, Upper error bounds on calculated outputs of interest for linear and non linear structural problems, C. R. Académie des Sci. - Méc. Paris 334 (2006) 399-407.

[16] P. Ladevèze, Strict upper error bounds on computed outputs of interest in computational structural mechanics, Comput. Mech. 42 (2008) 271-286.

[17] L. Chamoin, P. Ladevèze, A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems, Comput. Methods Appl. Mech. Engrg. 197 (2008) 994-1014.

[18] L. Chamoin, P. Ladevèze, Strict and practical bounds through a non-intrusive and goal-oriented error estimation method for linear viscoelasticity problems, Finite Elem. Anal. Design 45 (2009) 251-262.

[19] J. Panetier, P. Ladevèze, F. Louf, Strict bounds for computed stress intensity factors, Comput. Struct. 87 (2009) 1015-1021. Computational Structures Technology.

[20] P. Ladevèze, J. Waeytens, Model verification in dynamics through strict upper error bounds, Comput. Methods Appl. Mech. Engrg. 198 (2009) 1775-1784.

[21] P. Ladevèze, N. Moës, A posteriori constitutive relation error estimators for nonlinear finite element analysis and adaptive control, in: P. Ladevèze, J. Oden (Eds.), Advances in Adaptive Computational Methods in Mechanics, Studies in Applied Mechanics, vol. 47, Elsevier, 1998, pp. 231-256.

[22] P. Ladevèze, P. Marin, J. Pelle, J. Gastine, Accuracy and optimal meshes in finite element computation for nearly incompressible materials, Comput. Methods Appl. Mech. Engrg. 94 (1992) 303-315.

[23] R. Cottereau, P. Diez, A. Huerta, Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method, Comput. Mech. 44 (2009) 533-547.

[24] N. Parés, H. Santos, P. Díez, Guaranteed energy error bounds for the poisson equation using a flux-free approach: solving the local problems in subdomains, Int. J. Numer. Methods Engrg. 79 (2009) 1203-1244.
[25] P. Ladevèze, L. Chamoin, É. Florentin, A new non-intrusive technique for the construction of admissible stress fields in model verification, Comput. Methods Appl. Mech. Engrg. 199 (2009) 766-777.

[26] F. Pled, L. Chamoin, P. Ladevèze, On the techniques for constructing admissible stress fields in model verification: performances on engineering examples, Int. J. Numer. Methods Engrg. (2011).

[27] L. Chamoin, P. Ladevèze, Bounds on history-dependent or independent local quantities in viscoelasticity problems solved by approximate methods, Int. J. Numer. Methods Engrg. 71 (2007) 1387-1411.

[28] P. Ladevèze, L. Chamoin, Calculation of strict error bounds for finite element approximations of non-linear pointwise quantities of interest, Int. J. Numer. Methods Engrg. 84 (2010) 1638-1664.

[29] J. Waeytens, L. Chamoin, P. Ladevèze, Guaranteed error bounds on pointwise quantities of interest for transient viscodynamics problems, Comput. Mech. (2011) 1-17.

[30] P. Ladevèze, (1975) Comparaison de modéle de mileux continus. Thése de doctorat Université de Pierre et Marie Curie.

[31] I. Babuška, et W. Rheinboldt, Error estimates for adaptative finite element computation, SIAM Journal on Numeral Analysis 15 (4) (1978) 736-754.

[32] O.C. Zienkiewicz, et J.Z. Zhu, A simple error estimator and adaptative procedure for practical engineering analysis, International Journal for Numerical Methods in Engineering. 24 (1987) 337-357.