New bounding techniques for goal-oriented error estimation applied to linear problems

TL;DR

This paper introduces new bounding techniques based on Saint-Venant's principle and homothetic domains to improve the accuracy of goal-oriented error estimates in linear problems, demonstrated through numerical experiments.

Contribution

It presents a novel approach for error bounding in linear problems using properties of homothetic domains, enhancing the precision of error estimates.

Findings

Improved error bounds for linear problems.

Effective use of homothetic domains in error estimation.

Numerical experiments validate the proposed techniques.

Abstract

The paper deals with the accuracy of guaranteed error bounds on outputs of interest computed from approximate methods such as the finite element method. A considerable improvement is introduced for linear problems thanks to new bounding techniques based on Saint-Venant's principle. The main breakthrough of these optimized bounding techniques is the use of properties of homothetic domains which enables to cleverly derive guaranteed and accurate boundings of contributions to the global error estimate over a local region of the domain. Performances of these techniques are illustrated through several numerical experiments.

Tables

Table 1. Table 1: Calculated values of the quantities of interest

Quantity of interest	Exact value  $I_{ex}$	Approximate value  $I_h$
$I_1$	$0.0347969$	$0.0303747$
$I_2$	$-16.2939$	$-16.1596$
$I_3$	$2.8974$	$2.86161$

Table 2. Table 2: Values of constants h ℎ h and k 𝑘 k involved in both optimized improvements for different geometric shapes

Geometric shape	Constant  $h$		Constant  $k$
Dim 2	plane stress assumption	plane strain assumption
circle	$0.76923$	$0.7$	$2$
cracked circle with $\theta =0$	$0.79780$	$0.72122$	$2$
cracked circle with $\theta =\pi /6$	$0.80039$	$0.72315$	$2$
cracked circle with $\theta =\pi /3$	$0.80351$	$0.72546$	$2$
cracked circle with $\theta =\pi /2$	$0.80732$	$0.72829$	$2$
square	$0.85897$	$0.76667$	$2$
Dim 3
sphere	$0.53846$		$3$
parallelepiped	$0.64103$		$3$

Table 3. Table 3: Normalized bounds obtained using classical and optimized bounding techniques for quantity of interest I 2 subscript 𝐼 2 I_{2}

Number of layers	Classical bounds		Improved bounds 1		Improved bounds 2
of enriched nodes	${\overline{\xi }}_{\inf }$	${\overline{\xi }}_{\sup }$	${\overline{\chi }}_{\inf }$	${\overline{\chi }}_{\sup }$	${\overline{\zeta }}_{\inf }$	${\overline{\zeta }}_{\sup }$
$1$	$0.7168$	$1.2668$	$0.7291$	$1.2545$	$0.8512$	$1.1323$
$2$	$0.7567$	$1.2268$	$0.7387$	$1.2451$	$0.8543$	$1.1292$

Table 4. Table 4: Analytical values of constant h ℎ h for different geometric shapes

Geometric shape	Constant  $h$
Dim 2		plane stress assumption	plane strain assumption
unit circle	$\frac{2\mu +\lambda }{2\left(\mu +\lambda \right)}$	$\frac{1}{1+\nu }$	$1-\nu$
unit cracked circle111$\theta$ denotes the angle between the two lips of the crack.	$\frac{2\mu +\lambda }{2\left(\mu +\lambda \right)}+\frac{1}{3\left(2\pi -\theta \right)}\frac{\mu }{\mu +\lambda }$	$\frac{1}{1+\nu }+\frac{1}{3\left(2\pi -\theta \right)}\frac{1-\nu }{1+\nu }$	$1-\nu +\frac{1}{3\left(2\pi -\theta \right)}\left(1-2\nu \right)$
double unit square	$\frac{7\mu +3\lambda }{6\left(\mu +\lambda \right)}$	$\frac{7-\nu }{6\left(1+\nu \right)}$	$\frac{7-8\nu }{6}$
Dim 3
unit sphere	$\frac{2\mu +\lambda }{2\mu +3\lambda }$	$\frac{1-\nu }{1+\nu }$
double unit parallelepiped	$\frac{8\mu +3\lambda }{6\mu +9\lambda }$	$\frac{4-5\nu }{3\left(1+\nu \right)}$

Equations

$$\underline{u}\in\bm{\mathcal{U}};\quad\underline{u}_{\mathclose{}|\mathopen{}{\partial_{1}\Omega}}=\underline{U}_{d};\quad\depsil(\underline{u})=\frac{1}{2}\big{(}\nabla\underline{u}+\nabla^{T}\underline{u}\big{)};$$

$$\dsigma\in\bm{\mathcal{S}};\quad\forall\>\underline{u}^{\ast}\in\bm{\mathcal{U}}_{0},\quad\int_{\Omega}\operatorname{Tr}\big{[}\dsigma\>\depsil(\underline{u}^{\ast})\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega=\int_{\Omega}\underline{f}_{d}\cdot\underline{u}^{\ast}\mathop{}\mathopen{}\mathrm{d}\Omega+\int_{\partial_{2}\Omega}\underline{F}_{d}\cdot\underline{u}^{\ast}\mathop{}\mathopen{}\mathrm{d}S;$$

$$\dsigma(M)=\textbf{\em K}\>\depsil\big{(}\underline{u}(M)\big{)}\quad\forall\>M\in\Omega,$$

$$\displaystyle I=\mathcal{L}(\underline{u})=\int_{\Omega}\left(\operatorname{Tr}\big{[}\tilde{\dsigma}_{\Sigma}\>\depsil(\underline{u})\big{]}+\tilde{\underline{f}}_{\Sigma}\cdot\underline{u}\right)\mathop{}\mathopen{}\mathrm{d}\Omega,$$

$$\tilde{\underline{u}}\in {\mathcal{U}}_0;$$

$$\tilde{\dsigma}\in\bm{\mathcal{S}};\quad\forall\>\underline{u}^{\ast}\in\bm{\mathcal{U}}_{0},\quad\int_{\Omega}\operatorname{Tr}\big{[}\tilde{\dsigma}\>\depsil(\underline{u}^{\ast})\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega=\mathcal{L}(\underline{u}^{\ast})=\int_{\Omega}\left(\operatorname{Tr}\big{[}\tilde{\dsigma}_{\Sigma}\>\depsil(\underline{u}^{\ast})\big{]}+\tilde{\underline{f}}_{\Sigma}\cdot\underline{u}^{\ast}\right)\mathop{}\mathopen{}\mathrm{d}\Omega;$$

$$\tilde{\dsigma}(M)=\textbf{\em K}\>\depsil\big{(}\tilde{\underline{u}}(M)\big{)}\quad\forall\>M\in\Omega.$$

$$\displaystyle I_{ex}-I_{h}=\mathcal{L}(\underline{u})-\mathcal{L}(\underline{u}_{h})=\mathcal{L}(\underline{e}_{h})=\int_{\Omega}\operatorname{Tr}\big{[}\textbf{\em K}\>\depsil(\tilde{\underline{u}})\>\depsil(\underline{e}_{h})\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega=\mathcal{R}_{h}(\tilde{\underline{u}}).$$

$$\displaystyle\mathcal{L}(\underline{u})=\int_{\Omega}\operatorname{Tr}\big{[}\depsil(\underline{u})\>\textbf{\em K}\>\depsil(\tilde{\underline{u}})\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega=\mathcal{F}(\tilde{\underline{u}}),$$

$$\mathcal{F}(\tilde{\underline{u}})={\int }_{\Omega }{\underline{f}}_d\cdot \tilde{\underline{u}}\mathrm{d}\Omega +{\int }_{{\partial }_2\Omega }{\underline{F}}_d\cdot \tilde{\underline{u}}\mathrm{d}S.$$

$${\Vert {\underline{e}}_h\Vert }_{u,\Omega }^2={\Vert \underline{u}-{\hat{\underline{u}}}_h\Vert }_{u,\Omega }^2⩽{\Vert \underline{u}-{\hat{\underline{u}}}_h\Vert }_{u,\Omega }^2+{\Vert \text{dsigma}-{\widehat{\text{dsigma}}}_h\Vert }_{\text{dsigma},\Omega }^2=e_{\mathrm{cre}}^2,$$

$${\widehat{\text{dsigma}}}_h^m=\frac{1}{2}\left({\widehat{\text{dsigma}}}_h+\text{emK}\text{depsil}({\hat{\underline{u}}}_h)\right),$$

$${\Vert \text{dsigma}-{\widehat{\text{dsigma}}}_h^m\Vert }_{\text{dsigma},\Omega }=\frac{1}{2}{e}_{\mathrm{cre}}.$$

$$I_{ex}-I_h-I_{hh}$$

$$={⟨\text{dsigma}-{\widehat{\text{dsigma}}}_h^m,{\widehat{\tilde{\text{dsigma}}}}_h-\text{emK}\text{depsil}({\widehat{\tilde{\underline{u}}}}_h)⟩}_{\text{dsigma},\Omega },$$

$$I_{hh}$$

$$={⟨{\widehat{\tilde{\text{dsigma}}}}_h^m,{\widehat{\text{dsigma}}}_h-\text{emK}\text{depsil}({\hat{\underline{u}}}_h)⟩}_{\text{dsigma},\Omega }+\mathcal{L}({\hat{\underline{u}}}_h-{\underline{u}}_h),$$

$$\left|I_{ex}-I_h-I_{hh}\right|⩽\frac{1}{2}{e}_{\mathrm{cre}}{\tilde{e}}_{\mathrm{cre}}.$$

$${\xi }_{\inf }⩽I_{ex}⩽{\xi }_{\sup },$$

$${\xi }_{\inf }=I_h+I_{hh}-\frac{1}{2}{e}_{\mathrm{cre}}{\tilde{e}}_{\mathrm{cre}};$$

$${\xi }_{\sup }=I_h+I_{hh}+\frac{1}{2}{e}_{\mathrm{cre}}{\tilde{e}}_{\mathrm{cre}}.$$

$${\omega }_{\lambda }=\mathcal{H}[O;\lambda ]({\omega }_1)$$

$${\int }_{{\omega }_{\lambda }}\bullet \mathrm{d}\Omega$$

$$\frac{\mathrm{d}}{\mathrm{d}\lambda }\left[{\int }_{{\omega }_{\lambda }}\bullet \mathrm{d}\Omega \right]$$

$$\mathrm{d}S=\{\begin{array}{ll}\mathrm{d}s& \text{in 2D};\\ a(s_1,s_2)\mathrm{d}{s}_1\mathrm{d}{s}_2& \text{in 3D},\end{array}$$

$${\underline{v}}_{\lambda }=\{\begin{array}{ll}\frac{\lambda }{\bar{\lambda }}{\underline{v}}_{\bar{\lambda }}(s)\text{parameterized by}(\lambda ,s)& \text{in 2D};\\ \frac{\lambda }{\bar{\lambda }}{\underline{v}}_{\bar{\lambda }}(s_1,s_2)\text{parameterized by}(\lambda ,s_1,s_2)& \text{in 3D},\end{array}$$

$${\int }_{{\omega }_{\lambda }}\bullet \mathrm{d}\Omega$$

$${\int }_{\partial {\omega }_{\lambda }}\bullet \mathrm{d}S$$

$$n=\{\begin{array}{ll}1& \text{in 2D};\\ 2& \text{in 3D}.\end{array}$$

$${⟨\bullet ,\circ ⟩}_{u,{\omega }_{\lambda }}$$

$${⟨\bullet ,\circ ⟩}_{\text{dsigma},{\omega }_{\lambda }}$$

$${\left(\bullet ,\circ \right)}_{u,\partial {\omega }_{\lambda }}$$

$${\left(\bullet ,\circ \right)}_{\text{dsigma},\partial {\omega }_{\lambda }}$$

$$\frac{\mathrm{d}}{\mathrm{d}\lambda }\left[{\Vert \bullet \Vert }_{u,{\omega }_{\lambda }}^2\right]$$

Full text

1 Introduction

In the context of finite element model verification, research works currently focus on the evaluation of the numerical quality of specific quantities of practical interest. Such worthwhile methods, known as goal-oriented error estimation methods, have been emerging for about twenty years [1, 2, 3, 4, 5, 6, 7, 8] and have been recently extended to a wide range of mechanical problems [9, 10, 11, 12, 13, 14, 15, 16]. Nevertheless, among all of them, only a few actually lead to robust and relevant estimates ensuring the recovery of strict and high-quality error bounds.

In order to achieve robust goal-oriented error estimation, a general method, initially introduced in [13], has been prone to considerable developments. It is based on classical and powerful tools, such as the concept of constitutive relation error (CRE), and more recently handbook techniques [12, 17] and projection procedures [16]. First, methods based on the CRE enable to set up a guaranteed global error estimate through the construction of an admissible pair, which constitutes the key technical point. Various techniques enable to construct such an admissible solution from the prescribed data of the reference problem and the approximate finite element (FE) solution (see [18, 19, 20, 21, 22, 23, 24, 25] for further information). Then, a measure of the non-verification of the constitutive relation by this admissible pair (with respect to an energy norm) leads to an upper bound of the global discretization error. The goal-oriented error estimator considered in this work is based on extraction (or adjoint-based) techniques and involves the solution of an auxiliary problem, often referred to as dual (or adjoint in the linear case, or mirror in the non-linear case [13]) problem. An accurate solution of this auxiliary problem is required to achieve sharp local error bounds. A natural but intrusive way to properly solve this problem merely consists of performing a local space(-time) refinement of the mesh being considered. Alternatively, handbook techniques, initially proposed in [26], rely on a local enrichment of the solution of the auxiliary problem introduced through a partition of unity method (PUM). The enrichment functions correspond to (quasi-)exact local solutions of the auxiliary problem calculated either analytically or numerically in (semi-)infinite domains. This enrichment is particularly well-suited to handle truly pointwise quantities of interest without resorting to any regularization of the quantity of interest being studied. In practice, this technique preserves the non-intrusive nature offered by the FE framework as it circumvents the need to perform local refinement of the auxiliary problem to provide precise error bounds. Lastly, projection procedures allow the treatment of non-linear quantities of interest without having recourse to any classical linearization techniques and keeping intact the strict nature of the resulting bounds.

However, among those various tools, only the use of handbook techniques leads to the derivation of outstandingly accurate local error bounds, provided that appropriate enrichment functions are available in a library of pre-computed handbooks functions. Indeed, without additional techniques, the classical bounding technique may achieve low-quality guaranteed error bounds even with an extremely refined mesh. The reason why such error estimates are crude comes from the basic bounding technique, i.e. the Cauchy-Schwarz inequality, which is only efficient if terms involved in this inequality are close to be collinear, i.e. when the main contributions to the discretization errors related to both reference and auxiliary problems are located in the same regions. In other words, when the zone of interest does not coincide with the most concentrated error regions associated to the reference problem, the classical bounding technique yields inaccurate and sometimes useless bounds.

The present work sheds some light on this crucial question. In this paper, we revisit the classical bounding technique based on global error estimation methods applied to both reference and auxiliary problems, and introduce two new bounding techniques to alleviate problems related to the classical bounding technique. These improved techniques are based on both classical and innovative tools; they lean on Saint-Venant’s principle and are therefore restricted to linear problems. Basic extraction techniques and specific homotheticity properties are employed to get guaranteed and relevant bounds of better quality than the classical bounding technique. More precisely, the main idea consists of considering separately the zone of interest and the remainder of the structure. Both enhanced techniques use homothetic domains properties to cleverly derive accurate bounds over a local region surrounding the zone of interest. In this work, we examine only the case of linear quantities of interest with respect to the displacement field associated to linear elasticity problems. Finally, these new basic bounding techniques are combined with handbook techniques to solve in a very accurate manner the adjoint problem, thus resulting in efficient error bounds. 2D applications are performed and show the important resulting gain on the quality of the computed bounds using the two new basic bounding techniques.

The paper is organized as follows: after this introduction,

2

6 presents both reference and adjoint problems and defines the discretization error;

3

7 recalls basics on goal-oriented error estimation using the concept of admissible solutions and associated constitutive relation error; Sections 8 and 9 provide a detailed description of the first and second improved bounding techniques, respectively, both based on homothetic domains properties; several numerical experiments are presented in

4

15 to demonstrate the effectivity and relevance of the two proposed techniques; eventually,

5

17 draws some conclusions and may provide road maps to future developments.

6 Reference and adjoint problems

6.1 Reference problem and discretization error

Let us consider a mechanical structure occupying an open bounded domain $\Omega$, with Lipschitz boundary $\partial\Omega$. The prescribed loading acting on $\Omega$ consists of: a displacement field $\underline{U}_{d}$ on part $\partial_{1}\Omega\subset\partial\Omega$ $(\partial_{1}\Omega\neq\emptyset)$; a traction force density $\underline{F}_{d}$ on the complementary part $\partial_{2}\Omega$ of $\partial\Omega$ such that $\overline{\partial_{1}\Omega\cup\partial_{2}\Omega}=\partial\Omega$, $\partial_{1}\Omega\cap\partial_{2}\Omega=\emptyset$; a body force field $\underline{f}_{d}$ within $\Omega$. Structure $\Omega$ is assumed to be made of a material with isotropic, homogeneous, linear and elastic behavior characterized by Hooke’s tensor K. Assuming a quasi-static loading as well as isothermal and small perturbations state, the reference problem which describes the behavior of the structure consists of finding a displacement/stress pair $(\underline{u},\dsigma)$ in the space domain $\Omega$, which verifies:

• •

the kinematic conditions:

[TABLE]

• •

the weak form of equilibrium equations:

[TABLE]

• •

the constitutive relation:

[TABLE]

where $\depsil(\underline{u})$ represents the classical linearized strain tensor corresponding to the symmetric part of the gradient of displacement field $\underline{u}$. Affine spaces $\bm{\mathcal{U}}=\left\{\underline{u}\in[\mathcal{H}^{1}(\Omega)]^{3}\right\}$ and $\bm{\mathcal{S}}=\left\{\dsigma\in\mathcal{M}_{s}(3)\cap[\mathcal{L}^{2}(\Omega)]^{6}\right\}$ guarantee the existence of finite-energy solutions, $\mathcal{M}_{s}(n)$ representing the space of symmetric square matrices of order $n$. Lastly, $\bm{\mathcal{U}}_{0}\subset\bm{\mathcal{U}}$ denotes the vectorial space associated to $\bm{\mathcal{U}}$, i.e. containing the functions subjected to homogeneous kinematic boundary conditions over $\partial_{1}\Omega$.

In practical applications, the exact solution of the reference problem, hereafter denoted $(\underline{u},\dsigma)$, remains usually out of reach and only an approximate solution, referred to as $(\underline{u}_{h},\dsigma_{h})$, can be obtained through numerical approximation methods (such as the finite element method (FEM) associated with a space mesh $\mathcal{M}_{h}$ mapping $\Omega$). Such a numerical approximation is searched in a discretized space $\bm{\mathcal{U}}_{h}\!\times\!\bm{\mathcal{S}}_{h}\subset\bm{\mathcal{U}}\!\times\!\bm{\mathcal{S}}$. A displacement-type FEM leads to a displacement field $\underline{u}_{h}$ verifying kinematic constraints (1) and a stress field $\dsigma_{h}$ computed a posteriori from constitutive relation (3).

The resulting discretization error, denoted $\underline{e}_{h}=\underline{u}-\underline{u}_{h}$, can be evaluated by means of:

• •

a global measure, such as the classical energy norm \left\lVert\bullet\right\rVert_{u,\Omega}=\left(\int_{\Omega}\operatorname{Tr}\big{[}\textbf{\em K}\>\depsil(\bullet)\>\depsil(\bullet)\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega\right)^{1/2}, providing a global discretization error $e^{\mathrm{glob}}_{h}=\left\lVert\underline{e}_{h}\right\rVert_{u,\Omega}$.

• •

a local measure defined with respect to a specified output of interest $I(\underline{u})$ of the problem, providing a local error $e^{\mathrm{loc}}_{h}=I(\underline{u})-I(\underline{u}_{h})$. Under the assumption of a linear quantity of interest with respect to displacement $\underline{u}$, it merely reads: $e^{\mathrm{loc}}_{h}=I(\underline{e}_{h})$.

6.2 Adjoint problem

The quantity of interest, hereafter denoted $I$, is a goal-oriented output, such as the mean value of a stress component over a local region or the displacement value at a specific point, for instance. These meaningful quantities of engineering interest are usually defined by means of extraction techniques [2, 7, 27], i.e. by expressing the local quantity $I$ being considered in the global form involving global extraction operators, also called extractors. In this work, for the sake of simplicity, the quantity of interest is represented as a linear functional $\mathcal{L}$ of displacement field $\underline{u}$ on a finite support under the following global form:

[TABLE]

where so-called extractors $\tilde{\dsigma}_{\Sigma}$ and $\tilde{\underline{f}}_{\Sigma}$, known analytically, can be mechanically viewed as a prestress field and a body force field, respectively. These extractors can be defined explicitly or implicitly depending on the selected output of interest. Let us note that the extraction technique provides a natural framework to handle a wide range of local quantities. In the following, let $I_{ex}=\mathcal{L}(\underline{u})$ and $I_{h}=\mathcal{L}(\underline{u}_{h})$ be the unknown exact value of the quantity of interest $I$ being studied and its approximate value obtained through the FEM, respectively.

Once the quantity of interest has been put into such a global form, the classical approach then consists of introducing an auxiliary problem, also called adjoint problem, which is similar to the reference problem, except that the external mechanical loading $(\underline{F}_{d},\underline{f}_{d})$ is replaced by the extractors on the one hand, and the non-homogeneous Dirichlet boundary conditions are changed to homogeneous kinematic constraints on the other hand. The adjoint problem consists of finding a displacement/stress pair $(\tilde{\underline{u}},\tilde{\dsigma})$, in the space domain $\Omega$, which verifies:

• •

the kinematic conditions:

[TABLE]

• •

the weak form of equilibrium equations:

[TABLE]

• •

the constitutive relation:

[TABLE]

Under the assumption of a linear functional $\mathcal{L}$, the following equality holds:

[TABLE]

The solution $\tilde{\underline{u}}$ of the adjoint problem can thus be viewed as an influence function [5] indicating how the weak residual functional $\mathcal{R}_{h}$ (with respect to the discretization error $\underline{e}_{h}$) affects the local discretization error $\mathcal{L}(\underline{e}_{h})$ (with respect to the specific measure $\mathcal{L}$). Let us note that solutions $\underline{u}$ and $\tilde{\underline{u}}$ of reference and adjoint problems, respectively, are mutually adjoint to each other [8] insofar as

[TABLE]

where

[TABLE]

As for the reference problem, the exact solution $(\tilde{\underline{u}},\tilde{\dsigma})$ of the adjoint problem remains out of reach in most practical applications, and one can only obtain an approximate solution, denoted $(\tilde{\underline{u}}_{h},\tilde{\dsigma}_{h})$. This last solution lies in discretized FE spaces associated with a space mesh $\tilde{\mathcal{M}}_{h}$, mapping the physical domain $\Omega$, regardless of the FE mesh $\mathcal{M}_{h}$ used to solve the reference problem.

7 Basics on goal-oriented error estimation based on constitutive relation error

We review here the classical procedure to get guaranteed local error bounds on functional outputs which constitute valuable and relevant information in standard engineering practice.

7.1 Constitutive relation error

In verification research activities, setting up robust error estimation methods has become an overriding concern. The construction of what is called an admissible pair is currently an essential and crucial step in order to obtain guaranteed error bounds. An admissible pair $(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ verifies all the equations of the reference problem, apart from constitutive relation (3). On the one hand, a kinematically admissible displacement field is generally obtained by merely taking $\hat{\underline{u}}_{h}$ equal to $\underline{u}_{h}$ (apart from the case of incompressible materials, see [28]). On the other hand, the derivation of a statically admissible stress field can be achieved by using various balance techniques suitable to error estimation [18, 19, 20, 21, 22, 23, 24, 25]. Such an admissible stress field $\hat{\dsigma}_{h}$ can be recovered from the data and the FE stress field $\dsigma_{h}$ alone. Starting from an admissible solution $(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ provided by one of the existing techniques, one can measure the global residual on constitutive relation (3), called the constitutive relation error (CRE) and hereafter referred to as $e_{\mathrm{cre}}\equiv e_{\mathrm{cre}}(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})=\left\lVert\hat{\dsigma}_{h}-\textbf{\em K}\>\depsil(\hat{\underline{u}}_{h})\right\rVert_{\dsigma,\Omega}$, with \left\lVert\bullet\right\rVert_{\dsigma,\Omega}=\left(\int_{\Omega}\operatorname{Tr}\big{[}\bullet\>\textbf{\em K}^{-1}\>\bullet\big{]}\>\mathop{}\mathopen{}\mathrm{d}\Omega\right)^{1/2}. Computing the CRE $e_{\mathrm{cre}}(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ provides a guaranteed upper bound of the global discretization error $\left\lVert\underline{e}_{h}\right\rVert_{u,\Omega}$, as the well-known Prager-Synge hypercircle theorem [29] leads to the following bounding inequality:

[TABLE]

which conveys the guaranteed nature of the CRE $e_{\mathrm{cre}}$.

Introducing the average admissible field:

[TABLE]

one can directly deduce another fundamental relation, called the Prager-Synge’s equality:

[TABLE]

Equations (11) and (13) are key relations to derive guaranteed error bounds in both global and local robust error estimation methods.

In the same way as for the reference problem, an admissible solution of the adjoint problem, hereafter referred to as $(\hat{\tilde{\underline{u}}}_{h},\hat{\tilde{\dsigma}}_{h})$, can be derived from one of the existing equilibration techniques. Then, the associated CRE $\tilde{e}_{\mathrm{cre}}\equiv e_{\mathrm{cre}}(\hat{\tilde{\underline{u}}}_{h},\hat{\tilde{\dsigma}}_{h})$ of the adjoint problem can be computed leading to a global estimate of the discretization error $\tilde{\underline{e}}_{h}=\tilde{\underline{u}}-\tilde{\underline{u}}_{h}$ of the adjoint problem.

Now, let us focus on the main principles of the classical bounding technique involved in goal-oriented error estimation method based on extraction techniques and CRE.

7.2 Basic identity and classical bounding technique

The definition of the quantity of interest $I$ recast in the global form (4) and properties of both admissible solutions $(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ and $(\hat{\tilde{\underline{u}}}_{h},\hat{\tilde{\dsigma}}_{h})$ lead to the following basic identity:

[TABLE]

with $\displaystyle\hat{\dsigma}_{h}^{m}=\frac{1}{2}(\hat{\dsigma}_{h}+\textbf{\em K}\>\depsil(\hat{\underline{u}}_{h}))$; $\left\langle\bullet,\circ\right\rangle_{\dsigma,\Omega}$ denotes an energetic inner product defined on the stress field space over $\Omega$. $I_{hh}$ can be viewed as a computable correction term involving known quantities of both reference and adjoint problems:

[TABLE]

with $\displaystyle\hat{\tilde{\dsigma}}^{m}_{h}=\frac{1}{2}(\hat{\tilde{\dsigma}}_{h}+\textbf{\em K}\>\depsil(\hat{\tilde{\underline{u}}}_{h}))$, and leading to a new approximate solution $I_{h}+I_{hh}$ of the exact value $I_{ex}$ of the quantity of interest. A complete and detailed proof of this basic identity can be found in [13, 16].

The fundamental equality (14), which does not require any orthogonality property of the FE solutions (contrary to pioneering techniques [2, 5, 8]) and allows to decouple discretizations of reference and adjoint problems, is the cornerstone of the classical bounding technique as well as the improved ones described in section 8. Besides, this bounding technique could conceivably be extended to problems solved using numerical approximation methods different from the FEM.

Subsequently, the classical bounding procedure merely consists of applying the Cauchy-Schwarz inequality to (14) with respect to inner product $\left\langle\bullet,\circ\right\rangle_{\dsigma,\Omega}$ and then using Prager-Synge’s equality (13). This yields:

[TABLE]

Subsequently, the derivation of strict lower and upper bounds $(\xi_{\mathrm{inf}},\xi_{\mathrm{sup}})$ of $I_{ex}$ (or, equivalently, of the local error $I_{ex}-I_{h}$) can be achieved straightforwardly, just having a global error estimation procedure at hand:

[TABLE]

with

[TABLE]

Besides, owing to the independent natures of spatial discretizations associated to reference and adjoint problems, a convenient way to achieve accurate and sharp bounds of $I_{ex}$ is to perform a local space refinement of the adjoint mesh $\tilde{\mathcal{M}}_{h}$ alone around the zone of interest in order to properly solve the adjoint problem while keeping a reasonable computational cost. In most practical applications, the discretization error related to the adjoint problem is concentrated in the vicinity of the zone of interest, whereas that related to the reference problem may be scattered around zones which present some singularities or other error sources. However, when the error related to the reference problem is mostly located outside and far from the zone of interest, the classical bounding technique may yield large and low-quality local error bounds and thus makes useless bounding result (17). This is the point that we are revisiting here.

The proposed bounding techniques we present in the two following sections are intended to circumvent this serious drawback inherent to the classical technique in order to get sharp local error bounds.

8 First improved bounding technique

8.1 Definitions and preliminaries

Let us consider a reference subdomain, denoted $\omega_{1}$ and included in $\Omega$, defined by a point $O$ and a geometric shape. The set of homothetic domains $\omega_{\lambda}$ is defined as:

[TABLE]

where $\mathcal{H}[O;\lambda]$ stands for the homothetic transformation centered in point $O$, called homothetic center, and parameterized by a nonzero positive number $\lambda\in\mathopen{[}0\mathclose{}\mathpunct{},\lambda_{\max}\mathclose{]}$, also called magnification ratio, scale factor or similitude ratio, such that $\omega_{\lambda}\subset\Omega$, i.e. $\omega_{\lambda}$ is a subset of $\Omega$ (see Figure 1). The geometric shape defining the set of homothetic domains $\omega_{\lambda}$ is arbitrary, but in practice these physical domains are assumed to be basic, such as a circle or a rectangle in 2D, and a sphere or a rectangular cuboid (also called rectangular parallelepiped or right rectangular prism) in 3D, for instance.

By considering the parameterization $(\lambda,s)$ (resp. $(\lambda,s_{1},s_{2})$) of a given domain $\omega_{\lambda}\subset\Omega$ in 2D (resp. 3D), where $s$ denotes the curvilinear abscissa along boundary $\partial\omega_{\lambda}$ (see Figure 1), the following equalities hold:

[TABLE]

where $\mathop{}\mathopen{}\mathrm{d}S$ is defined in a generic manner as:

[TABLE]

and depends on the geometric shape defining the set of homothetic domains.

For a given pair $(\omega_{\lambda},\omega_{\bar{\lambda}})$ of homothetic domains included in $\Omega$, represented in Figure 1 and parameterized by $(\lambda,\bar{\lambda})$, such that $\omega_{\lambda}\subset\omega_{\bar{\lambda}}\subset\Omega$, i.e. $\lambda\in\mathopen{]}0\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$, the position $\underline{v}_{\lambda}$ of a point $M_{\lambda}$ along $\partial\omega_{\lambda}$ can be defined from the position $\underline{v}_{\bar{\lambda}}$ of the corresponding point $M_{\bar{\lambda}}$ along $\partial\omega_{\bar{\lambda}}$ by the following relation:

[TABLE]

where $s$ (resp. $s_{1}$ and $s_{2}$) represent the curvilinear abscissa along boundary $\partial\omega_{\bar{\lambda}}$ in 2D (resp. 3D).

Such a parameterization leads to the following relations:

[TABLE]

where $n$ denotes a positive integer defined as:

[TABLE]

Let us define inner products and associated norms over a given homothetic domain $\omega_{\lambda}\subset\Omega$:

[TABLE]

Similarly, let us also define inner products and associated norms over boundary $\partial\omega_{\lambda}$ of a given homothetic domain $\omega_{\lambda}$, such that $\omega_{\lambda}\subset\omega_{\bar{\lambda}}\subset\Omega$:

[TABLE]

One can easily prove that the following equalities hold:

[TABLE]

where $n$ is defined by (27).

Proof.

Starting from definition (28) of $\left\lVert\bullet\right\rVert^{2}_{u,\omega_{\lambda}}$, then using relation (25) and definition (30) of $\left\lvert\bullet\right\rvert^{2}_{u,\partial\omega_{\bar{\lambda}}}$, one gets directly:

[TABLE]

Eventually, differentiation of (34) with respect to variable $\lambda$ completes the proof of relation (32). Similarly, the derivation of relation (33) can be proved in a straightforward manner. ∎

Definition 1**.**

Measures $e_{\mathrm{cre},\lambda}(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ and $e_{\mathrm{cre},\lambda}(\hat{\tilde{\underline{u}}}_{h},\hat{\tilde{\dsigma}}_{h})$ of the non-verification of the constitutive relations related to both reference and adjoint problems in part $\omega_{\lambda}\subset\Omega$ are defined by:

[TABLE]

and

[TABLE]

respectively.

Similarly, measures $e_{\mathrm{cre},\setminus\lambda}(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ and $e_{\mathrm{cre},\setminus\lambda}(\hat{\tilde{\underline{u}}}_{h},\hat{\tilde{\dsigma}}_{h})$ of the non-verification of the constitutive relations related to both reference and adjoint problems in the complementary part $\Omega\setminus\omega_{\lambda}$ are defined by:

[TABLE]

and

[TABLE]

respectively.

8.2 General idea

First, let us recall that the quantity $q$ to bound for building an upper error bound is (see (14)):

[TABLE]

where $\hat{\dsigma}_{h}^{m}$ and $\hat{\tilde{\dsigma}}_{h}-\textbf{\em K}\>\depsil(\hat{\tilde{\underline{u}}}_{h})$ are given quantities and $\dsigma$ is the unknown exact stress solution of the reference problem.

Let us consider a subdomain $\omega_{\lambda}$ of domain $\Omega$, whose complementary part is denoted by $\Omega\setminus\omega_{\lambda}$ in the following. Then, $q$ can be split into two parts:

[TABLE]

where

[TABLE]

and

[TABLE]

If quantity $\hat{\tilde{\dsigma}}_{h}-\textbf{\em K}\>\depsil(\hat{\tilde{\underline{u}}}_{h})$ is concentrated over $\omega_{\lambda}$, i.e. if subdomain $\omega_{\lambda}$ surrounds the zone of interest $\omega$, part $q_{\setminus\lambda}$ can be merely bounded as follows:

[TABLE]

$\tilde{e}_{\mathrm{cre},\setminus\lambda}$ being a relatively small computable term. It follows that the main contribution to the error comes from $q_{\lambda}$. Consequently, quantity $q_{\lambda}$ has to be correctly bounded.

In order to derive accurate bounds for part $q_{\lambda}$, the discretization error $\underline{e}_{h}=\underline{u}-\hat{\underline{u}}_{h}$ on $\omega_{\lambda}$ is split into:

• •

a local error, denoted $\underline{u}_{1}$;

• •

a pollution error, denoted $\underline{u}_{2}$,

which are solutions of the following two local problems defined on $\omega_{\lambda}$, referred to as $(\mathcal{P}_{1})$ and $(\mathcal{P}_{2})$, respectively:

1. $\bullet$

Problem $(\mathcal{P}_{1})$ consists of searching $(\underline{u}_{1},\dsigma_{1})\in\bm{\mathcal{U}}\!\times\!\bm{\mathcal{S}}$ such that:

• $\circ$

$\underline{u}_{1}=\underline{0}\quad\text{on}\ \partial\omega_{\lambda}$

• $\circ$

$\underline{\operatorname{div}}(\dsigma_{1})=\underline{\operatorname{div}}(\textbf{\em K}\>\depsil(\underline{u}-\hat{\underline{u}}_{h}))$

• $\circ$

$\dsigma_{1}=\textbf{\em K}\>\depsil(\underline{u}_{1})$ 2. $\bullet$

Problem $(\mathcal{P}_{2})$ consists of searching $(\underline{u}_{2},\dsigma_{2})\in\bm{\mathcal{U}}\!\times\!\bm{\mathcal{S}}$ such that:

• $\circ$

$\underline{u}_{2}=\underline{u}-\hat{\underline{u}}_{h}\quad\text{on}\ \partial\omega_{\lambda}$

• $\circ$

$\underline{\operatorname{div}}(\dsigma_{2})=\underline{0}$

• $\circ$

$\dsigma_{2}=\textbf{\em K}\>\depsil(\underline{u}_{2})$

Decomposition $\underline{u}-\hat{\underline{u}}_{h}=\underline{u}_{1}+\underline{u}_{2}$ is the starting point for deriving the main technical result presented in the next section. It is worthy noticing that there is no need to perform the extra-resolutions of local problems $(\mathcal{P}_{1})$ and $(\mathcal{P}_{2})$ to get the main technical and final bounding results (49) and (53) introduced in Sections 8.3 and 8.4, respectively. This point is discussed and proved in Appendix A.

8.3 Main technical result

Let us consider the space $\bm{\mathcal{V}}$ of functions satisfying equilibrium conditions:

[TABLE]

and let us introduce the Steklov constant, or Steklov eigenvalue, $h$ defined in [30] as:

[TABLE]

with

[TABLE]

where subdomain $\omega_{1}$ denotes the homothetic domain $\omega_{\lambda}$ associated to a constant parameter $\lambda=1$. Then, for any homothetic domain $\omega_{\lambda}\subset\Omega$ parametrized by $\lambda>0$, one can derive a relation involving the product of constant $h$ and parameter $\lambda$:

[TABLE]

with

[TABLE]

Proposition 1**.**

Let $(\omega_{\lambda},\omega_{\bar{\lambda}})$ be a pair of homothetic domains such that $\lambda\in\mathopen{]}0\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$, i.e. $\omega_{\lambda}\subset\omega_{\bar{\lambda}}$. The following key inequality holds:

[TABLE]

where

[TABLE]

The proof of Proposition 1 is given in Appendix A.

Let us note that, using Prager-Synge’s equality (13), unknown term $\left\lVert\dsigma-\hat{\dsigma}_{h}\right\rVert_{\dsigma,\omega_{\bar{\lambda}}}$ involved in the right-hand side term of fundamental inequality (49) in Proposition 1 is readily bounded as:

[TABLE]

It follows that fundamental result (49) can be rewritten in terms of perfectly known quantities as:

[TABLE]

Remark 1**.**

Quantity $\left\lVert\dsigma-\hat{\dsigma}_{h}\right\rVert_{\dsigma,\omega_{\bar{\lambda}}}$ can be straightforwardly bounded without introducing $\hat{\dsigma}^{m}_{h}$ by simply using $\left\lVert\dsigma-\hat{\dsigma}_{h}\right\rVert^{2}_{\dsigma,\omega_{\bar{\lambda}}}\leqslant\left\lVert\dsigma-\hat{\dsigma}_{h}\right\rVert^{2}_{\dsigma,\Omega}\leqslant e^{2}_{\mathrm{cre}}$. Nevertheless, the corresponding bound, namely $e^{2}_{\mathrm{cre}}$, is less accurate than the one given by (51), namely $1/4\left(e_{\mathrm{cre}}+e_{\mathrm{cre},\bar{\lambda}}\right)^{2}$.

8.4 Final bounding result

Proposition 2**.**

The final improved bounding result reads:

[TABLE]

where

[TABLE]

and

[TABLE]

are all computable from the calculated approximate solutions of both reference and adjoint problems. $I_{h}+I_{hh}+I_{hhh}$ can be viewed as a new approximate solution of the exact value $I_{ex}$ of the quantity of interest.

Proof.

First, let us address the question of bounding of part $q_{\setminus\lambda}$. Starting from inequation (43) and using Prager-Synge’s equality (13), quantity $q_{\setminus\lambda}$ can be bounded as follows:

[TABLE]

Second, let us now handle the question of bounding of part $q_{\lambda}$. Equation (41) can be rewritten as follows:

[TABLE]

where $I_{hhh}$ is a calculable known term defined in (55).

Applying the Cauchy-Schwarz inequality to (58) with respect to scalar product $\left\langle\bullet,\circ\right\rangle_{\dsigma,\omega_{\lambda}}$, one has:

[TABLE]

Then, introducing the key inequality (52) coming from Proposition 1 into (59) leads to the following bounding result:

[TABLE]

Finally, using both inequalities (56) and (60), one gets:

[TABLE]

which completes the proof of Proposition 6. ∎

Thus, this improved technique provides the following guaranteed lower and upper bounds $(\chi_{\mathrm{inf}},\chi_{\mathrm{sup}})$ of $I_{ex}$:

[TABLE]

with

[TABLE]

Remark 2**.**

These bounds depend on both parameters $\lambda$ and $\bar{\lambda}$. In order to get a practical minimizer, one seeks to reduce ratio $\displaystyle\frac{\lambda}{\bar{\lambda}}$ as much as possible by choosing:

• •

the smallest parameter $\lambda$ such that domain $\omega_{\lambda}$ surrounds the zone of interest $\omega$;

• •

the largest parameter $\bar{\lambda}$ such that domain $\omega_{\bar{\lambda}}$ remains a homothetic mapping of $\omega_{\lambda}$ (preserving its geometric shape) contained in $\Omega$,

and leading to sharp error bounds.

Remark 3**.**

*The Steklov constant $h$ could be easily computed, $\bm{\mathcal{S}}_{1}(\underline{v})$ being a Rayleigh quotient associated with a symmetric eigenproblem. By considering a material with isotropic, homogeneous, linear and elastic behavior, in the two-dimensional (three-dimensional, respectively) case of a unit circle111a unit circle (sphere, respectively) corresponds to a circular (spherical, respectively) domain of radius one. (unit sphere111a unit circle (sphere, respectively) corresponds to a circular (spherical, respectively) domain of radius one., respectively) $\omega_{1}$, it has been shown numerically that the maximum of $\bm{\mathcal{S}}_{1}$ is reached for $\underline{v}=\underline{OM}$, where $O$ is the homothetic center and $M\in\partial\omega_{1}$; it follows that $\depsil(\underline{v})=\mathbb{I}_{d}$. Let us note that in the two-dimensional case of a unit cracked circle and of a double unit square222a double unit square (parallelepiped, respectively) corresponds to a squared (parallelepiped, respectively) domain of side length two. as well as in the three-dimensional case of a double unit parallelepiped222a double unit square (parallelepiped, respectively) corresponds to a squared (parallelepiped, respectively) domain of side length two., the same maximum eigenfunction for $\bm{\mathcal{S}}_{1}$ has been obtained numerically. Analytical expressions of constant $h$ for various shape domains are reported in Table 4 (see Appendix C). Note that, in the particular shape domains we considered, constant $h$ only depends on Poisson’s ratio $\nu$. Besides, lower bounds of this constant for a circular shape domain are given in [31].

9 Second improved bounding technique

9.1 General idea

As for the first improved bounding method presented in

10

8, let us consider a subdomain $\omega_{\bar{\lambda}}\subset\Omega$, with complementary part $\Omega\setminus\omega_{\bar{\lambda}}$. Quantity $q$ previously defined in (39) can be decomposed as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

Similarly to the previous improved technique, quantity $q_{\setminus\bar{\lambda}}$ can be easily bounded as follows:

[TABLE]

whereas the bounding of $q_{\bar{\lambda}}$ differs widely. As previously mentioned in

11

8.2, the main contribution to the error derives from $q_{\bar{\lambda}}$.

In order to build sharp bounds for part $q_{\bar{\lambda}}$, let us introduce the following local problem defined on $\omega_{\bar{\lambda}}\subset\Omega$, referred to as $(\mathcal{P}^{h}_{ex})$, which consists of searching $(\underline{u}^{h}_{ex},\dsigma^{h}_{ex})\in\bm{\mathcal{U}}\!\times\!\bm{\mathcal{S}}$ such that:

• $\circ$

$\underline{u}^{h}_{ex}=\hat{\underline{u}}_{h}\quad\text{on}\ \partial\omega_{\bar{\lambda}}$

• $\circ$

$\underline{\operatorname{div}}(\dsigma^{h}_{ex})+\underline{f}_{d}=\underline{0}$

• $\circ$

$\dsigma^{h}_{ex}=\textbf{\em K}\>\depsil(\underline{u}^{h}_{ex})$

Similarly to local problems $(\mathcal{P}_{1})$ and $(\mathcal{P}_{2})$ previously defined in

12

8.2 for the first improved bounding technique, one can mention that the final bounding result (85) introduced in

13

13.3 does not require the solution of local problem $(\mathcal{P}^{h}_{ex})$, as the present improved technique circumvents the need to perform any additional resolution of $(\mathcal{P}^{h}_{ex})$.

It follows that the discretization error $\underline{e}_{h}=\underline{u}-\hat{\underline{u}}_{h}$ and $\dsigma-\hat{\dsigma}_{h}$ on $\omega_{\bar{\lambda}}$ could be rewritten:

[TABLE]

and

[TABLE]

respectively.

Then, quantity $q_{\bar{\lambda}}$ can be decomposed into two parts:

[TABLE]

where

[TABLE]

and

[TABLE]

It will be shown that second part $q_{\bar{\lambda},2}$ could be easily bounded, while a sharp bound of first part $q_{\bar{\lambda},1}$ will be derived thanks to Saint-Venant’s principle.

13.1 Bounding of part $q_{\bar{\lambda},2}$

First, let us address the question of bounding of part $q_{\bar{\lambda},2}$. Applying the Cauchy-Schwarz inequality to (73) with respect to scalar product $\left\langle\bullet,\circ\right\rangle_{\dsigma,\omega_{\bar{\lambda}}}$ leads to:

[TABLE]

Proposition 3**.**

Let $(\underline{u}^{h}_{ex},\dsigma^{h}_{ex})\in\bm{\mathcal{U}}\!\times\!\bm{\mathcal{S}}$ be the exact solution of problem $(\mathcal{P}^{h}_{ex})$ defined over $\omega_{\bar{\lambda}}$. Given an admissible approximate solution $(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ of the reference problem, the following equality holds:

[TABLE]

Let us note that replacing $\dsigma^{h}_{ex}$ by the exact solution $\dsigma$ of the reference problem in Proposition 3 leads generally to a incorrect result except over the whole domain $\Omega$, where Prager-Synge’s equality (13) holds.

Proof.

Noticing that the restriction of $(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ to $\omega_{\bar{\lambda}}$ is also an admissible solution of problem $(\mathcal{P}^{h}_{ex})$ defined over $\omega_{\bar{\lambda}}$, the well-known Prager-Synge hypercircle theorem [29] leads to the following bounding inequality:

[TABLE]

which conveys the guaranteed nature of the CRE $e_{\mathrm{cre},\bar{\lambda}}$ on $\omega_{\bar{\lambda}}$ for problem $(\mathcal{P}^{h}_{ex})$.

Introducing the average admissible field $\hat{\dsigma}^{m}_{h}$ defined by (12) completes the proof of Proposition 3. ∎

Finally, incorporating result (75) of Proposition 3 into (74), one gets:

[TABLE]

From now on, in order to get accurate bounds for part $q_{\bar{\lambda},1}$, let us introduce the main technical result.

13.2 Main technical result

Let us consider the space $\bm{\mathcal{V}}$ introduced in (44) and let us define the following dimensionless constant:

[TABLE]

with

[TABLE]

for a given homothetic domain $\omega_{\bar{\lambda}}\subset\Omega$ associated to parameter $\bar{\lambda}>0$. Then, for any domain $\omega_{\lambda}$ homothetic to $\omega_{\bar{\lambda}}$ such that $\omega_{\lambda}\subset\omega_{\bar{\lambda}}\subset\Omega$, i.e. for any $\lambda\in\mathopen{]}0\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$, one can derive a relation involving the product of constant $k$ and ratio $\displaystyle\frac{\bar{\lambda}}{\lambda}$:

[TABLE]

with

[TABLE]

Remark 4**.**

This problem is connected to the description of a 2D or 3D homogeneous domain as an abstract beam for which a semi-group could be defined [32]. However, we do not know any name for this constant which is strictly positive for star-shaped domains such that $\underline{v}_{\bar{\lambda}}\cdot\underline{n}>0$.

Proposition 4**.**

Let $(\omega_{\lambda},\omega_{\bar{\lambda}})$ be a pair of homothetic domains such that $\lambda\in\mathopen{]}0\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$, i.e. $\omega_{\lambda}\subset\omega_{\bar{\lambda}}\subset\Omega$. The following key inequality holds:

[TABLE]

The proof of Proposition 4 is given in Appendix B.

The bounding of part $q_{\bar{\lambda},1}$ involves quantity $\underline{u}-\underline{u}^{h}_{ex}$ (resp. $\dsigma-\dsigma^{h}_{ex}$) of decomposition (69) (resp. (70)) of the discretization error. Recalling that $\underline{u}-\underline{u}^{h}_{ex}\in\bm{\mathcal{V}}$, applying fundamental result (82) of Proposition 4 to $\underline{u}-\underline{u}^{h}_{ex}$ leads to:

[TABLE]

or, equivalently:

[TABLE]

13.3 Final bounding result

Proposition 5**.**

The final improved bounding result reads:

[TABLE]

where

[TABLE]

and

[TABLE]

are fully calculable from the calculated approximate solution of adjoint problem.

Proof.

Let us first handle the question of bounding of part $q_{\bar{\lambda},1}$, which is actually the key point of this new technique. Using relation (25) for domain $\omega_{\bar{\lambda}}$ itself and applying the Cauchy-Schwarz inequality with respect to scalar product $\left(\bullet,\circ\right)_{\dsigma,\partial\omega_{\bar{\lambda}}}$, one has:

[TABLE]

and

[TABLE]

In order to derive an upper bound of (89), let us introduce a function $\mu(\lambda)>0$:

[TABLE]

Subsequently, applying result (33) to $\dsigma-\dsigma^{h}_{ex}$ in the left-hand side term of (90) leads to:

[TABLE]

Integrating the first term of the left-hand side of (91) by parts, one obtains:

[TABLE]

Assuming that $\displaystyle\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d}\lambda}\left[\mu(\lambda)\right]\leqslant 0\quad\forall\>\lambda\in\mathopen{]}0\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$, using fundamental result (84) coming from Proposition 4 leads to:

[TABLE]

Integrating the second term of the right-hand side of (93) by parts, one obtains:

[TABLE]

thus resulting in the following inequality:

[TABLE]

The optimal function $\mu$ is the one which minimizes the right-hand side of (95):

[TABLE]

It reads:

[TABLE]

under the assumption that $\mu$ is a strictly positive, monotonically decreasing function of $\lambda$ on the interval $\mathopen{]}0\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$, which implies the following condition:

[TABLE]

Remark 5**.**

Note that, in the two- and three-dimensional cases of circular and spherical domains, respectively, condition (98) holds (see Remark 7).

Then, replacing function $\mu(\lambda)$ by expression (97) into (95), one gets:

[TABLE]

where $\tilde{\theta}_{\bar{\lambda}}$ is defined by (86). Finally, collecting bounding results (99) and (77) for $q_{\bar{\lambda},1}$ and $q_{\bar{\lambda},2}$, respectively, $q_{\bar{\lambda}}$ can be bounded as follows:

[TABLE]

Corollary 1**.**

Let $(\underline{u}^{h}_{ex},\dsigma^{h}_{ex})\in\bm{\mathcal{U}}\!\times\!\bm{\mathcal{S}}$ be the exact solution of problem $(\mathcal{P}^{h}_{ex})$ and $(\underline{u},\dsigma)\in\bm{\mathcal{U}}\!\times\!\bm{\mathcal{S}}$ the one of the reference problem. Given an admissible approximate solution $(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ of the reference problem, the following equality holds:

[TABLE]

Proof.

Using result (75) of Proposition 3, one obtains:

[TABLE]

with

[TABLE]

and

[TABLE]

Eventually, we end up the proof of Corollary 1 with:

[TABLE]

∎

It follows from result (101) of Corollary 1:

[TABLE]

it can be rewritten in the following form:

[TABLE]

thereby getting back to the bounding (65) of quantity $q$ by using bounding results (107) and (68) for $q_{\bar{\lambda}}$ and $q_{\setminus\bar{\lambda}}$, respectively, one obtains:

[TABLE]

Let us now introduce a scalar $\nu>0$ in order to shrewdly regroup $\left\lVert\dsigma-\hat{\dsigma}_{h}^{m}\right\rVert_{\dsigma,\omega_{\bar{\lambda}}}$ and $\left\lVert\dsigma-\hat{\dsigma}_{h}^{m}\right\rVert_{\dsigma,\Omega\setminus\omega_{\bar{\lambda}}}$, which are such that:

[TABLE]

one has:

[TABLE]

The optimal scalar $\nu$ is the one which minimizes the right-hand side of (110):

[TABLE]

It reads:

[TABLE]

which is a strictly positive scalar.

Then, changing scalar $\nu$ by expression (112) into inequation (110), one gets:

[TABLE]

which proves result (85) of Proposition 5. ∎

Thus, this improved technique provides the following guaranteed lower and upper bounds $(\zeta_{\mathrm{inf}},\zeta_{\mathrm{sup}})$ of $I_{ex}$:

[TABLE]

with

[TABLE]

Remark 6**.**

These bounds involve only one parameter $\bar{\lambda}$, while bounds $(\chi_{\mathrm{inf}},\chi_{\mathrm{sup}})$ defined by (63) and (64) in

14

8.3* depend on two parameters $(\lambda,\bar{\lambda})$. In practice, one can determine an optimum value $\bar{\lambda}_{\mathrm{opt}}$ for parameter $\bar{\lambda}$ that minimizes upper bound (113) of $q$ in order to obtain very precise bounds $(\zeta_{\mathrm{inf}},\zeta_{\mathrm{sup}})$ of $I_{ex}$. Besides, practically subdomain $\omega_{\bar{\lambda}}$ should recover the zone where the solution of adjoint problem has stiff gradients. *

Remark 7**.**

*The constant $k$ could be easily computed, $\bm{\mathcal{R}}_{\bar{\lambda}}(\underline{v})$ being a Rayleigh quotient associated with a symmetric eigenproblem. In the two-dimensional (three-dimensional, respectively) case of a circle111circle (sphere, respectively) $\omega_{\bar{\lambda}}$ corresponds to a circular (spherical, respectively) domain of radius $\bar{\lambda}$. (sphere111circle (sphere, respectively) $\omega_{\bar{\lambda}}$ corresponds to a circular (spherical, respectively) domain of radius $\bar{\lambda}$., respectively) $\omega_{\bar{\lambda}}$, it has been shown numerically that the minimum of $\bm{\mathcal{R}}_{\bar{\lambda}}$ is reached for $\underline{v}=\underline{OM}$, where $O$ is the homothetic center and $M\in\partial\omega_{\bar{\lambda}}$; furthermore, for a circular and spherical geometric shapes, it has been shown in [33] (see Appendix D); it follows that $\depsil(\underline{v})=\mathbb{I}_{d}$. Let us note that in the two-dimensional case of a cracked circle and of a double square222double square (parallelepiped, respectively) $\omega_{\bar{\lambda}}$ corresponds to a squared (parallelepiped, respectively) domain of side length $2\bar{\lambda}$., as well as in the three-dimensional case of a double unit parallelepiped222double square (parallelepiped, respectively) $\omega_{\bar{\lambda}}$ corresponds to a squared (parallelepiped, respectively) domain of side length $2\bar{\lambda}$., the same minimum eigenfunction for $\bm{\mathcal{R}}_{\bar{\lambda}}$ has been obtained numerically.

For all aforementioned two- (three-, respectively) dimensional cases, constant $k$ is equal to $2$ ($3$, respectively). Note that, in these particular shape domains, constant $k$ does not depend on any material parameter.

Remark 8**.**

In the present work, we restrict to only one subdomain, but the results could be easily extended to several ones. Multiple subdomains could be associated to different zones where the local CRE is relatively high either for adjoint or reference problem.

15 Numerical results

All numerical experiments have been performed assuming that the material remains isotropic, homogeneous, linear and elastic with Young’s modulus $E=1$ and Poisson’s ratio $\nu=0.3$. Furthermore, the two-dimensional examples are assumed to satisfy the plane-stress approximation. The balance technique used to derive a statically admissible stress field is the element equilibration technique (EET) combined with a $p$-refinement technique consisting of a $p+k$ discretization, $p$ being the FE interpolation degree and $k$ an additional degree equal to $3$ (see principles of the different techniques for constructing admissible stress fields in [24] for more information).

Performances of the proposed bounding techniques are illustrated through a two-dimensional cracked structure, already considered in [19, 23, 25].

The results obtained for classical bounding technique as well as first and second improved variants are presented in terms of the normalized bounds $(\bar{\xi}_{\mathrm{inf}},\bar{\xi}_{\mathrm{sup}})$, $(\bar{\chi}_{\mathrm{inf}},\bar{\chi}_{\mathrm{sup}})$, $(\bar{\zeta}_{\mathrm{inf}},\bar{\zeta}_{\mathrm{sup}})$, respectively, defined by:

[TABLE]

15.1 Presentation of the 2D cracked structure

Let us consider the two-dimensional structure shown in Figure 2, which presents two round cavities. A homogeneous Dirichlet boundary condition is imposed to the bigger circular hole, whereas a unit internal constant pressure $p_{0}$ is applied to the smaller one. Furthermore, the top-left edge is subjected to a unit normal traction force density $\underline{t}=+\underline{n}$. Besides, a single edge crack emanates from the bottom of the smaller cavity. The two lips of this crack as well as the remaining sides are traction-free boundaries. The FE mesh consists of $7\,751$ linear triangular elements and $4\,122$ nodes (i.e. $8\,244$ d.o.f.), see Figure 2. It has been adaptively refined in the vicinity of the crack tip. The reference mesh used to compute an “overkill” solution and to define a “quasi-exact” value, denoted $I_{ex}$ for convenience, of the quantity of interest is built up by dividing each element into $256$ elements; thereby, it is made of $1\,984\,256$ linear triangular elements and $996\,080$ nodes (i.e. $1\,992\,160$ d.o.f.).

The quantities of interest being considered in this work are:

• •

the average value in a local zone $\omega\subset\Omega$ of the component $(.)_{xx}$ of the stress field $\dsigma$:

[TABLE]

where extraction domain $\omega$ corresponds to a finite element $E$ of FE mesh $\mathcal{M}_{h}$ illustrated in Figure 3 and $\left\lvert\omega\right\rvert$ represents its measure;

• •

the pointwise value of the component $(.)_{x}$ of the displacement field $\underline{u}$ at a point $P$:

[TABLE]

where point $P$ coincides with a node of FE mesh $\mathcal{M}_{h}$ illustrated in Figure 3, the corresponding extraction domain $\omega$ being pointwise;

• •

the stress intensity factor $K_{\mathrm{I}}$ involved in the finite-energy analytical asymptotic expression of the stress field $\dsigma$ in the vicinity of the crack tip [10, 15] and classically used in crack propagation criteria for linear elastic fracture mechanics (LEFM) problems:

[TABLE]

All the local quantities being considered are linear functions of displacement field $\underline{u}$ associated to reference problem. Exact values $I_{ex}$ and approximate values $I_{h}$ obtained for the three considered quantities of interest are reported in Table 1.

In the case of pointwise quantity of interest $I_{2}$, for which the loading of the adjoint problem is a pointwise force defined using a Dirac-type function applied to a node of FE mesh $\mathcal{M}_{h}$, a natural regularization would consist in treating the pointwise force as a nodal force $\tilde{\underline{f}}_{\Sigma}=\delta(\underline{x}_{P})\tilde{\underline{F}}_{\Sigma}$, where $\tilde{\underline{F}}_{\Sigma}=\displaystyle\begin{bmatrix}1\\ 0\end{bmatrix}$ and $\underline{x}_{P}$ denotes the position of point $P$. Let us note that this approach can be employed only in the case of pointwise values of the displacement field related to the position of a node $P$ of the FE mesh and leads to a coarse approximate solution of the adjoint problem. Another classical way to handle truly pointwise quantities of interest is to have recourse to the use of mollification, also called mollifying process, replacing the initial quantity of interest by a local weighted average value [5]. A last approach, henceforth known as handbook techniques [12, 17, 16], consists in introducing a local enrichment of the solution of the adjoint problem particularly well-suited to pointwise error estimation without performing any regularization of the pointwise quantity of interest being considered. The singularities involved in the loading of the adjoint problem are captured explicitly in a non-intrusive way introducing adapted enrichment functions, also called handbook functions and denoted $(\tilde{\underline{u}}^{hand},\tilde{\dsigma}^{hand})$, through the PUM (satisfied by the linear FE shape functions) allowing to ensure displacement compatibility. These functions represent local (quasi-)exact solutions of the adjoint problem. In the context of pointwise error estimation, they correspond to the well-known and possibly infinite-energy Green functions describing the singular solutions of the adjoint problem over an infinite (or semi-infinite) domain. Several examples of such functions and a detailed description of handbook techniques can be found in [12, 17, 16]. The Green functions can be determined analytically using approaches based on strain nuclei combined with the image method [34, 35, 36, 37, 38]; for quantity of interest $I_{2}$, handbook functions are given in Figure 4.

Enrichment is applied via the PUM to a region $\Omega^{PUM}$ of domain $\Omega$, which is decomposed into two non-overlapping subregions: a zone $\Omega_{1}^{PUM}$ containing the zone of interest $\omega$ over which the quantity of interest is defined and a complementary zone $\Omega_{2}^{PUM}$ surrounding $\Omega_{1}^{PUM}$. In the case of quantity of interest $I_{2}$, an enrichment by the PUM involving only one or two layers of nodes is sufficient to capture the local high gradients of the exact solution $\tilde{u}$ of the adjoint problem; the definition of subregions $\Omega_{1}^{PUM}$ (containing $6$ (resp. $24$) elements and $7$ (resp. $19$) enriched nodes) and $\Omega_{2}^{PUM}$ (containing $18$ (resp. $30$) elements) in region $\Omega^{PUM}$ is given in Figure 5 by considering a single layer (resp. two layers) of nodes involved in the enrichment.

Finally, the global solution $(\tilde{\underline{u}},\tilde{\dsigma})$ of the adjoint problem is composed of an enrichment part $(\tilde{\underline{u}}^{hand}_{PUM},\tilde{\dsigma}^{hand}_{PUM})$ introduced a priori either analytically or numerically and a residual part $(\tilde{\underline{u}}^{res},\tilde{\dsigma}^{res})$ computed a posteriori numerically by using the FEM. The new adjoint problem consists in finding the residual pair $(\tilde{\underline{u}}^{res},\tilde{\dsigma}^{res})$. The corresponding force vector involved in the equilibrium equations comes down to a traction load $-\tilde{\dsigma}^{hand}\>\underline{n}_{12}$ over $\partial\Omega_{1}^{PUM}$, where $\underline{n}_{12}$ is the outgoing normal vector to $\Omega_{1}^{PUM}$, and a prestress $-\tilde{\dsigma}^{hand}_{PUM}$ over $\Omega_{2}^{PUM}$ by using properties of the handbook functions. Let us note that the loading of the adjoint problem becomes finite and smooth and that a correct and satisfactory approximate solution can be merely obtained by using the same spatial discretization as the one employed for the reference problem. Besides, observing that enrichment part $(\tilde{\underline{u}}^{hand}_{PUM},\tilde{\dsigma}^{hand}_{PUM})$ satisfies the constitutive relation over $\Omega$, the CRE $e_{\mathrm{cre}}(\hat{\tilde{\underline{u}}}_{h},\hat{\tilde{\dsigma}}_{h})$ of initial approximate adjoint problem is exactly equal to the CRE $e_{\mathrm{cre}}(\hat{\tilde{\underline{u}}}^{res}_{h},\hat{\tilde{\dsigma}}^{res}_{h})$ of residual approximate adjoint problem.

In the case of quantity of interest $I_{3}$, open-mode stress intensity factor $K_{\mathrm{I}}$ (i.e. associated with the Mode-I loading) is defined from an extension of the contour integral method proposed in [39]. The interested reader can refer to [39, 1, 10, 40, 15] for a detailed description of this technique. The calculation of stress intensity factor $K_{\mathrm{I}}$ comes down to evaluating a surface integral over an arbitrary crown $\Omega_{c}$ surrounding the crack tip:

[TABLE]

where $(\underline{v}_{\mathrm{I}},\dsigma_{\mathrm{I}})$ corresponds to a singular (infinite-energy) analytical solution of the elastic problem defined near the crack tip and $\phi$ is an arbitrary continuously derivable scalar function that is equal to [math] (resp. $1$) along external boundary $\Gamma_{e}$ (resp. internal boundary $\Gamma_{i}$) of crown $\Omega_{c}$. One usually considers a circular annulus $\Omega_{c}$ centered around the crack tip and a linear continuous function $\phi$ [10, 40, 15], but other polynomial functions $\phi$ could be considered [1]. Here, internal part $\Gamma_{i}$ and external part $\Gamma_{e}$ of boundary $\partial\Omega_{c}$ are assumed to be circles centered at the crack tip with radii $R_{i}=6$ and $R_{e}=8$, respectively; the definition of crown $\Omega_{c}$ is given in Figure 6.

Function $\phi$ is chosen to be linear with radius $r$ and is defined as:

[TABLE]

Let us mention that taking a polynomial function $\phi$ with a degree at most compatible with that of the FE analysis is appropriate, since it can be correctly represented in FE framework. The effectiveness and accuracy of the classical bounding technique to derive lower and upper bounds for open- and shear-mode stress intensity factors with a standard FEM have been demonstrated on two-dimensional benchmark problems in [10].

The loading of the adjoint problems involves the following extractors:

• •

a uniform prestress field $\tilde{\dsigma}_{\Sigma}=\textbf{\em K}\>\tilde{\depsil}_{\Sigma}$ over element $E$, where $\tilde{\depsil}_{\Sigma}=\displaystyle\frac{1}{\left\lvert E\right\rvert}\begin{bmatrix}1&0\\ 0&0\end{bmatrix}$ and $\left\lvert E\right\rvert$ denotes the measure of element $E$, for quantity of interest $I_{1}$;

• •

a traction load $\tilde{\underline{F}}_{\Sigma}=-\tilde{\dsigma}^{hand}\>\underline{n}_{12}$ over $\partial\Omega_{1}^{PUM}$ and a prestress field $\tilde{\dsigma}_{\Sigma}=-\tilde{\dsigma}^{hand}_{PUM}$ over $\Omega_{2}^{PUM}$ for quantity of interest $I_{2}$;

• •

a prestress field $\tilde{\dsigma}_{\Sigma}=\textbf{\em K}\>\depsil(\phi\>\underline{v}_{\mathrm{I}})-\phi\>\dsigma_{\mathrm{I}}$ and a body force field $\tilde{\underline{f}}_{\Sigma}=-\dsigma_{\mathrm{I}}\>\nabla\phi$ over crown $\Omega_{c}$ for quantity of interest $I_{3}$.

In this work, circular shape domains are considered for quantities of interest $I_{1}$ and $I_{2}$, while cracked circular shape domains are used for quantity of interest $I_{3}$, but other geometric shapes could have been investigated, the main technical aspect being the calculation of the constants $h$ and $k$ involved in the first and second improved techniques, respectively. Values of constants $h$ and $k$ have been calculated analytically and computed numerically for different geometric shapes by considering an isotropic, homogeneous, linear and elastic material with Poisson’s ratio $\nu=0.3$. Results are given in Table 2.

15.2 Average value of a field over a local zone: case of quantity of interest $I_{1}$

Maps of the local contributions to the global error estimates for both reference and adjoint problems are displayed in Figure 7. The adjoint mesh has been slightly refined around the zone of interest $\omega$ (corresponding to element $E$), and contains $8\,973$ linear triangular elements and $4\,733$ nodes (i.e. $9\,466$ d.o.f.). The main contributions to the error estimate associated to reference problem are by a majority located near the crack tip, while that associated to adjoint problem are concentrated around the zone of interest $\omega$. Therefore, the error estimates for both reference and adjoint problems are localized in disjoint regions.

The values of parameters $\lambda$ and $\bar{\lambda}$ involved in the first improved technique are set to $2r$ and $14r$, respectively, where $r$ corresponds to the radius of the circle circumscribed by element $E$. The value of parameter $\bar{\lambda}_{\mathrm{opt}}$ involved in the second improved technique is set to $9r$, which enables to achieve the sharpest bounds for quantity of interest $I_{1}$. The corresponding subdomains $\omega_{\lambda}$, $\omega_{\bar{\lambda}}$ and $\omega_{\bar{\lambda}_{\mathrm{opt}}}$ are illustrated in Figure 8.

Figure 9 shows the evolutions of the normalized bounds on $I_{ex}$ for quantity of interest $I_{1}$ as functions of the number of elements $\tilde{N}_{e}$ of adjoint problem for the classical bounding technique as well as the two improved ones. The adjoint mesh has been locally refined near the zone of interest $\omega$, since the loading and the contributions to the global error estimate of the adjoint problem are highly localized in this region. One can see a slight improvement in the bounds obtained with the first improved technique compared to the classical one. As regards the second improved technique, a very clear improvement is observed allowing to achieve sharp local error bounds without refining too much the adjoint problem, thus keeping an affordable computing time. Similar results can be obtained for the average values in element $E$ of the other components $(.)_{xy}$ and $(.)_{yy}$ of stress field $\dsigma$, as well as for the ones of the different components of strain field $\depsil$.

Figure 10 represents the evolutions of the normalized exact value $I_{ex}/I_{ex}$ of local quantity $I_{1}$, its normalized approximate value $I_{h}/I_{ex}$ obtained through the FEM and its new normalized approximate value $(I_{h}+I_{hh})/I_{ex}$ as functions of the number of elements $\tilde{N}_{e}$ of adjoint problem. One can see that $I_{h}+I_{hh}$ corresponds to an approximation of better quality with respect to $I_{ex}$ compared to $I_{h}$.

15.3 Pointwise value of a field: case of quantity of interest $I_{2}$

In order to properly solve the adjoint problem and to capture the singularities of its solution, specific pre-calculated handbook functions are introduced locally in the neighborhood of the loading of the adjoint problem, i.e. in the vicinity of point $P$. Let us recall that handbook techniques enable to achieve accurate local error bounds by using exactly the same spatial mesh for both adjoint and reference problems. The graph represented in Figure 11 shows the evolution of the normalized bounds $(\bar{\xi}_{\mathrm{sup}},\bar{\xi}_{\mathrm{inf}})$ of $I_{ex}$, obtained using the classical bounding technique coupled with handbook techniques, with respect to the number of enriched nodes. The corresponding number of layers of enriched nodes ranges from $1$ to $7$. By using the same spatial discretization for both residual adjoint and reference problems, meaningful and effective bounds can be obtained with an enrichment via the PUM involving only few nodes. In the following, we restrict to an enrichment applied only to one or two layers of nodes without performing any refinement of the residual adjoint problem.

Maps of the local contributions to the global error estimates for (residual) adjoint problem with both enrichments are displayed in Figure 12. The main contributions to the error estimate associated to residual adjoint problem are concentrated in zone $\Omega_{2}^{PUM}$, which constitutes the support of the loading of residual adjoint problem. Similarly to the previous case presented in

16

15.2, the error estimates for both reference and residual adjoint problems are located in distant regions.

For the enrichment involving a single layer of nodes, the values of parameters $\lambda$ and $\bar{\lambda}$ involved in the first improved technique are set to $1.7r$ and $7r$, respectively, where $r$ corresponds to the distance between point $P$ (defining the pointwise zone of interest) and its farthest neighbor node. The value of parameter $\bar{\lambda}_{\mathrm{opt}}$ involved in the second improved technique is set to $4.4r$, which enables to achieve the sharpest bounds for quantity of interest $I_{2}$. The corresponding subdomains $\omega_{\lambda}$, $\omega_{\bar{\lambda}}$ and $\omega_{\bar{\lambda}_{\mathrm{opt}}}$ are illustrated in Figure 13.

For the enrichment involving two layers of nodes, the values of parameters $\lambda$ and $\bar{\lambda}$ involved in the first improved technique are set to $2.5r$ and $7r$, while the value of parameter $\bar{\lambda}_{\mathrm{opt}}$ involved in the second improved technique is set to $4.4r$, which enables to achieve the sharpest bounds for quantity of interest $I_{2}$. The corresponding subdomains $\omega_{\lambda}$, $\omega_{\bar{\lambda}}$ and $\omega_{\bar{\lambda}_{\mathrm{opt}}}$ are illustrated in Figure 14.

The normalized bounds obtained using classical and improved bounding techniques combined with handbook techniques for an enrichment by the PUM involving either one or two layers of nodes are summarized in Table 3.

The analysis of the results reveals that bounds $(\bar{\zeta}_{\mathrm{inf}},\bar{\zeta}_{\mathrm{sup}})$ obtained using the second improved technique are more accurate than the ones $(\bar{\xi}_{\mathrm{inf}},\bar{\xi}_{\mathrm{sup}})$ obtained using the classical technique as well as the ones $(\bar{\chi}_{\mathrm{inf}},\bar{\chi}_{\mathrm{sup}})$ obtained using the first improved technique for both considered enrichments. It is worth noticing that bounds $(\bar{\chi}_{\mathrm{inf}},\bar{\chi}_{\mathrm{sup}})$ obtained using the first improved technique are coarser than the ones $(\bar{\xi}_{\mathrm{inf}},\bar{\xi}_{\mathrm{sup}})$ obtained using the classical technique for the second enrichment. Indeed, performances of the first improved technique strongly depend on ratio $\displaystyle\frac{\lambda}{\bar{\lambda}}$; the larger the enriched zone is, the higher ratio $\displaystyle\frac{\lambda}{\bar{\lambda}}$ is.

Similar results can be obtained for the value at point $P$ of the other component $(.)_{y}$ of displacement field $\underline{u}$.

16.1 Extracted value of a field: case of quantity of interest $I_{3}$

Map of the local contributions to the global error estimate for adjoint problem is depicted in Figure 15. The adjoint mesh density has been slightly increased toward the crown $\Omega_{c}$; the adjoint mesh is made of $10\,699$ linear triangular elements and $5\,597$ nodes (i.e. $11\,194$ d.o.f.). The highest contributions to the error estimate associated to adjoint problem are localized along the boundary of crown $\Omega_{c}$. Therefore, the error estimates for both reference and adjoint problems are localized in close regions, contrary to previous cases presented in Sections 15.2 and 15.3.

The values of parameters $\lambda$ and $\bar{\lambda}$ involved in the first improved technique are set to $1.2R_{e}$ and $2.1R_{e}$, respectively, where $R_{e}$ corresponds to the radius of external circle $\Gamma_{e}$. The value of parameter $\bar{\lambda}_{\mathrm{opt}}$ involved in the second improved technique is set to $1.6R_{e}$, which enables to achieve the sharpest bounds for quantity of interest $I_{3}$. The corresponding subdomains $\omega_{\lambda}$, $\omega_{\bar{\lambda}}$ and $\omega_{\bar{\lambda}_{\mathrm{opt}}}$ are illustrated in Figure 16.

Figure 17 shows the evolutions of the normalized bounds on $I_{ex}$ for quantity of interest $I_{3}$ as functions of the number of elements $\tilde{N}_{e}$ of adjoint problem for the classical bounding technique as well as the two improved ones. A local refinement of adjoint mesh has been performed adaptively near the zone of interest $\omega$, especially along internal and external boundaries of crown $\Omega_{c}$, since the contributions to the global error estimate of the adjoint problem are mainly located in this region. One can see a moderate decrease in the bounds obtained with the first improved technique compared to the classical one., while the second improved technique and the classical one give similar results as regards the accuracy of the bounds. Consequently, the improvements presented in this work may lead to useless bounds in cases where the major part of estimated error related to both reference and adjoint problems are localized in close regions. Similar results can be obtained for the shear-mode stress intensity factor $K_{\mathrm{II}}$.

Figure 18 represents the evolutions of the normalized exact value $I_{ex}/I_{ex}$ of local quantity $I_{3}$, its normalized approximate value $I_{h}/I_{ex}$ obtained through the FEM and its new normalized approximate value $(I_{h}+I_{hh})/I_{ex}$ as functions of the number of elements $\tilde{N}_{e}$ of adjoint problem. One can observe once again that $I_{h}+I_{hh}$ corresponds to an approximation of better quality with respect to $I_{ex}$ compared to $I_{h}$.

17 Conclusion and prospects

In this paper, we introduced two new approaches related to the general framework of robust goal-oriented error estimation dealing with extraction techniques. These techniques are based on mathematical tools which are not classical in model verification. Several numerical experiments clearly demonstrate the efficiency of these methods to produce strict and relevant bounds on the errors in linear local quantities of interest compared to the classical bounding technique, especially when the discretization error related to the reference error is not concentrated in the local zone of interest. Nevertheless, the second proposed technique seems to achieve sharper local error estimates than the first one.

Finally, such powerful methods may open up opportunities and help widen the field of robust goal-oriented error estimation methods. Both techniques could be easily extended to other quantities of interest but are restricted to linear problems, i.e. cases where Saint-Venant’s principle is well established. The extension to broader classes of mechanical problems such as time-dependent non linear problems is thus an open question. Besides, coupling these new improved methods with handbook techniques offers several potentially fruitful directions for future research. This crucial issue is still under active consideration.

Appendix A. Proof of the main technical result of Proposition 1

First, let us recall that the CRE in subdomain $\omega_{\lambda}$, included in $\Omega$, defined by (35) can be recast in the following form:

[TABLE]

Let us now consider the following decomposition of the discretization error $\underline{e}_{h}=\underline{u}-\hat{\underline{u}}_{h}$ on subdomain $\omega_{\lambda}$:

[TABLE]

where the local error $\underline{u}_{1}$ and the pollution error $\underline{u}_{2}$ are solutions of local problems $(\mathcal{P}_{1})$ and $(\mathcal{P}_{2})$, respectively, defined on $\omega_{\lambda}$ and introduced in

18

8.2.

Such a decomposition verifies

[TABLE]

and

[TABLE]

where ${(\bullet)}_{\mathrm{sym}}$ represents the symmetric part of matrix $\bullet$.

Then, relation (126) becomes:

[TABLE]

By observing that $\underline{u}_{2}\in\bm{\mathcal{V}}$, where space $\bm{\mathcal{V}}$ has been defined in (44), let us now introduce the Steklov constant $h\lambda$ defined in (47), which leads to:

[TABLE]

As a result, the following inequality can be deduced from (133) and (134):

[TABLE]

Besides, using equality (22) for quantity \operatorname{Tr}\big{[}(\dsigma-\hat{\dsigma}_{h})\>\textbf{\em K}^{-1}\>(\dsigma-\hat{\dsigma}_{h})\big{]}, one obtains:

[TABLE]

Incorporating (136) into (135) and rearranging the terms, inequality (135) can be rewritten as follows:

[TABLE]

where $\displaystyle\alpha_{\lambda}\equiv\alpha_{\lambda}(\hat{\underline{u}}_{h},\hat{\dsigma}_{h},\underline{u}_{1},\underline{u}_{2})=e^{2}_{\mathrm{cre},\lambda}-\left\lVert\underline{u}_{1}\right\rVert^{2}_{u,\omega_{\lambda}}-\left[\left\lVert\underline{u}_{2}\right\rVert_{u,\omega_{\lambda}}-\sqrt{h\lambda}\left\lVert\dsigma-\hat{\dsigma}_{h}\right\rVert_{\dsigma,\partial\omega_{\lambda}}\right]^{2}$. Let us note that first order ordinary differential inequation (137) verified by $\left\lVert\dsigma-\hat{\dsigma}_{h}\right\rVert^{2}_{\dsigma,\omega_{\lambda}}$ can be reformulated as:

[TABLE]

where $\displaystyle f\colon\lambda\mapsto\lambda^{-1/h}$ is a strictly positive function for any $\lambda>0$. Finally, a worthwhile relation between $\left\lVert\dsigma-\hat{\dsigma}_{h}\right\rVert^{2}_{\dsigma,\omega_{\lambda}}$ and $\left\lVert\dsigma-\hat{\dsigma}_{h}\right\rVert^{2}_{\dsigma,\omega_{\bar{\lambda}}}$ can be derived integrating inequality (138) over $\mathopen{[}\lambda\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$. It reads:

[TABLE]

where $\displaystyle\beta_{\lambda,\bar{\lambda}}\equiv\beta_{\lambda,\bar{\lambda}}(\hat{\underline{u}}_{h},\hat{\dsigma}_{h},\underline{u}_{1},\underline{u}_{2})=\int_{\lambda^{\prime}=\lambda}^{\bar{\lambda}}\left[\left(\frac{\lambda^{\prime}}{\lambda}\right)^{-1/h}\frac{1}{h\lambda^{\prime}}\alpha_{\lambda^{\prime}}\right]\mathop{}\mathopen{}\mathrm{d}\lambda^{\prime}$ is a function involving displacement solutions $\underline{u}_{1}$ and $\underline{u}_{2}$ of problems $(\mathcal{P}_{1})$ and $(\mathcal{P}_{2})$, respectively. In practice, inequality $\displaystyle\alpha_{\lambda^{\prime}}\leqslant e^{2}_{\mathrm{cre},\lambda^{\prime}}\quad\forall\>\lambda^{\prime}\in\mathopen{[}\lambda\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$ allows to work around the extra resolutions of problems $(\mathcal{P}_{1})$ and $(\mathcal{P}_{2})$ and leads to the following fundamental result of Proposition 1:

[TABLE]

where $\displaystyle\gamma_{\lambda,\bar{\lambda}}$ is completely defined by (50) in terms of computable known quantities. For practical purposes, one-dimensional numerical integration methods, such as the simple trapezoidal rule or the classical Gauss-Legendre quadrature rule, can be employed to get an accurate approximation of $\gamma_{\lambda,\bar{\lambda}}$.

Remark 9**.**

Other expressions for $\gamma_{\lambda,\bar{\lambda}}$ can be derived directly by integration by parts; it readily reads:

[TABLE]

However, implementation of this latter expression is more cumbersome compared to the former one, since it requires the integration of the first derivative of function $e^{2}_{\mathrm{cre},\lambda^{\prime}}$ instead of function $e^{2}_{\mathrm{cre},\lambda^{\prime}}$ itself. Nevertheless, applying relation (33) to quantity $\hat{\dsigma}_{h}-\textbf{\em K}\>\depsil(\hat{\underline{u}}_{h})$ for any homothetic domain $\omega_{\lambda^{\prime}}$ such that $\lambda^{\prime}\in\mathopen{[}\lambda\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$, one obtains:

[TABLE]

consequently, using relation (25) for both domains $\omega_{\lambda}$ and $\omega_{\bar{\lambda}}$, one has:

[TABLE]

where the last term of the right-hand side:

[TABLE]

can be viewed as a weighted constitutive relation error in $\omega_{\bar{\lambda}}\setminus\omega_{\lambda}$. Even though this last relation (143) does not call for one-dimensional numerical integration methods contrary to (50) and (141), it requires the numerical evaluation of a definite integral of a rational function over $\omega_{\bar{\lambda}}\setminus\omega_{\lambda}$.

Eventually, in order to perform an accurate calculation of function $\gamma_{\lambda,\bar{\lambda}}$, expression (50) computed by means of a basic trapezoidal integration method with a large number of integration points is preferred among the different aforementioned expressions.

Appendix B. Proof of the main technical result of Proposition 4

Let us consider $\underline{v}\in\bm{\mathcal{V}}$. One can deduce the following inequality from (80):

[TABLE]

Then, using (26), one has: $\displaystyle\left\lvert\underline{v}\right\rvert^{2}_{u,\partial\omega_{\lambda}}={\left(\frac{\lambda}{\displaystyle\bar{\lambda}}\right)}^{n}\left\lvert\underline{v}\right\rvert^{2}_{u,\partial\omega_{\bar{\lambda}}}$ and (146) becomes:

[TABLE]

According to result (32), one has: $\displaystyle{\left(\frac{\lambda}{\displaystyle\bar{\lambda}}\right)}^{n}\frac{1}{\displaystyle\bar{\lambda}}\left\lvert\underline{v}\right\rvert^{2}_{u,\partial\omega_{\bar{\lambda}}}=\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d}\lambda}\left[\left\lVert\underline{v}\right\rVert^{2}_{u,\omega_{\lambda}}\right]$; consequently, the following first order ordinary differential inequation satisfied by $\left\lVert\underline{v}\right\rVert^{2}_{u,\omega_{\lambda}}$ holds:

[TABLE]

which can be recast as:

[TABLE]

where $\displaystyle g\colon\lambda\mapsto\lambda^{-k}$ is a strictly positive function for any $\lambda>0$; as a consequence, one gets:

[TABLE]

Finally, a fundamental result connecting $\left\lVert\underline{v}\right\rVert^{2}_{u,\omega_{\lambda}}$ and $\left\lVert\underline{v}\right\rVert^{2}_{u,\omega_{\bar{\lambda}}}$ can be derived integrating inequality (150) over $\mathopen{[}\lambda\mathclose{}\mathpunct{},\bar{\lambda}\mathclose{]}$. It reads:

[TABLE]

which completes the proof of Proposition 4.

Appendix C. Analytical values of constant $h$ involved in the first improved technique

The values of constant $h$, expressed in terms of Lamé’s cœfficients $(\lambda,\mu)$ in column two and in terms of Poisson’s ratio $\nu$ in the last columns (under both plane stress and plane strain assumptions for two-dimensional geometric shapes), are given in Table 4.

Appendix D. Analytical computation of constant $k$ involved in the second improved technique for circular and spherical shape domains

Let us consider the space $\bm{\mathcal{V}}$ of functions satisfying homogeneous equilibrium equations expressed in displacements, also known as Lamé-Navier equations. Trefftz (or T-) functions are homogeneous solutions of the governing differential equations inside the domain (corresponding to Lamé-Navier equations in our case). Trefftz-type approaches, which consist in using a set of linearly independent solutions of a differential equation, were initially introduced by Trefftz [41]. T-functions were firstly employed as basis functions in Trefftz methods [42]. Nowadays, they are classically used as interpolation functions to define Trefftz-type finite elements in (hybrid) FEM, Boundary Element Method (BEM) also called Boundary Integral Equation Method (BIEM); they can be found in the form of polynomials, Legendre, harmonic, exponential, Bessel, Hankel, Kelvin (also called singular Kupradze), Boussinesq functions, depending on the governing equations [43, 44, 45]. Some special purpose T-functions enable to satisfy not only the governing equations, but also take into account special boundary conditions a priori. The interested reader can refer to [46, 43] for a complete description of basic sets of T-functions associated to Laplace, Helmholtz and biharmonic equations in 2D and 3D problems, both for bounded and unbounded domains. Combinations of those classical T-functions allow the derivation of T-functions associated to a broad class of problems, such as 2D and 3D elasticity or Mindlin-type plates.

Let us consider a circular (resp. spherical) shape domain $\omega_{\bar{\lambda}}$ of radius $\bar{\lambda}$ in 2D (resp. 3D). The set of T-functions associated to Lamé-Navier equations constitutes a basis of functions belonging to space $\bm{\mathcal{V}}$ on $\omega_{\bar{\lambda}}$ [33]. Thus, any homogeneous solution $\underline{v}\in\bm{\mathcal{V}}$ can be defined by linear combination of a set of Trefftz functions. Suitable Trefftz functions corresponding to solutions of Lamé-Navier equations can be found in cartesian coordinates system in [47] and in polar coordinates system in [48]. T-functions associated to Lamé-Navier equations can be expressed in cartesian coordinate system $(x,y)$ in 2D (resp. $(x,y,z)$ in 3D) as polynomials [47]. Transforming cartesian coordinates into polar coordinates $(r,\theta)$ in 2D (resp. spherical coordinates $(r,\theta,\phi)$ in 3D), displacement T-functions $\underline{T}$ can be formulated as regular harmonic polynomials of the form:

[TABLE]

where functions $f_{2}$ and $g_{2}$ (resp. $f_{3}$, $g_{3}$ and $h_{3}$) depend on material parameters $\nu$, $E$ as well as polar coordinate $\theta$ in 2D (resp. spherical coordinates $\theta,\phi$ in 3D) only. Let us consider T-function $\underline{T}^{(m)}$ of maximum degree $m$ in polar (resp. spherical) coordinate $r$, or, equivalently, in cartesian coordinates $(x,y)$ (resp. $(x,y,z)$) in 2D (resp. 3D), with $m\geqslant 1$ (i.e. discarding rigid body motions). Corresponding stress T-function $\textbf{\em K}\>\depsil(\underline{T}^{(m)})$ can be derived in a polynomial form of degree $m-1$; it follows that energy e^{(m)}\equiv\operatorname{Tr}\big{[}\depsil(\underline{T}^{(m)})\>\textbf{\em K}\>\depsil(\underline{T}^{(m)})\big{]} can be put in the following polynomial form:

[TABLE]

where function $\Psi_{2}$ (resp. $\Psi_{3}$) depends on material parameters $\nu$, $E$ as well as polar coordinate $\theta$ in 2D (resp. spherical coordinates $\theta,\phi$ in 3D) only. Thus, after some straightforward computations, corresponding Rayleigh quotient $\bm{\mathcal{R}}_{\bar{\lambda}}(\underline{T}^{(m)})$ reads as:

[TABLE]

Therefore, any displacement T-function $\underline{T}$ associated to non-vanishing strain satisfies $\displaystyle\bm{\mathcal{R}}_{\bar{\lambda}}(\underline{T})\geqslant 2$ in 2D and $\displaystyle\bm{\mathcal{R}}_{\bar{\lambda}}(\underline{T})\geqslant 3$ in 3D. Besides, let us note that T-functions are also orthogonal with respect to both inner products $\left(\bullet,\circ\right)_{u,\partial\omega_{\bar{\lambda}}}$ and $\left\langle\bullet,\circ\right\rangle_{u,\omega_{\bar{\lambda}}}$ related to Rayleigh quotient $\bm{\mathcal{R}}_{\bar{\lambda}}$.

Finally, given that any function $\underline{v}\in\bm{\mathcal{V}}$ (i.e. any homogeneous solution of Lamé-Navier equations) can be approximated as a linear combination of T-functions, one can show that:

[TABLE]