On the techniques for constructing admissible stress fields in model verification: Performances on engineering examples

TL;DR

This paper compares three techniques for constructing admissible stress fields in finite element analysis, evaluating their accuracy, computational cost, and ease of implementation through industrial case studies.

Contribution

It provides a comprehensive comparison of EET, SPET, and EESPT methods, highlighting the hybrid EESPT as an effective compromise for engineering applications.

Findings

SPET is more accurate but computationally costlier.

EESPT offers a good balance between accuracy and practicality.

All methods show consistent convergence rates.

Abstract

Robust global/goal-oriented error estimation is used nowadays to control the approximate finite element solutions obtained from simulation. In the context of Computational Mechanics, the construction of admissible stress fields (\ie stress tensors which verify the equilibrium equations) is required to set up strict and guaranteed error bounds (using residual based error estimators) and plays an important role in the quality of the error estimates. This work focuses on the different procedures used in the calculation of admissible stress fields, which is a crucial and technically complicated point. The three main techniques that currently exist, called the element equilibration technique (EET), the star-patch equilibration technique (SPET), and the element equilibration + star-patch technique (EESPT), are investigated and compared with respect to three different criteria, namely the quality of associated error estimators, computational cost and easiness of practical implementation into commercial finite element codes. The numerical results which are presented focus on industrial problems; they highlight the main advantages and drawbacks of the different methods and show that the behavior of the three estimators, which have the same convergence rate as the exact global error, is consistent. Two- and three-dimensional experiments have been carried out in order to compare the performance and the computational cost of the three different approaches. The analysis of the results reveals that the SPET is more accurate than EET and EESPT methods, but the corresponding computational cost is higher. Overall, the numerical tests prove the interest of the hybrid method EESPT and show that it is a correct compromise between quality of the error estimate, practical implementation and computational cost.

Tables

Table 1. Table 1: Comparison of the error estimators given by the EET, the SPET, and the EESPT.

\topruleMethods	Estimate  $\theta$	Effectivity index  $\eta$	Normalized CPU time
\midruleEET	$812.999$	$2.3362$	$1.000$
SPET	$556.629$	$1.5995$	$4.218$
EESPT	$812.801$	$2.3357$	$1.156$
\bottomrule

Table 2. Table 2: Influence of the choice of the cost function involved in the EET and the EESPT on the quality of the error estimators.

\toprule		Estimate  $\theta$				Effectivity index  $\eta$
\ms      \msMethods		$J_0$	$J_1$	$J_2$		$J_0$	$J_1$	$J_2$
\midruleEET		$812.999$	$812.801$	$749.732$		$2.3362$	$2.3357$	$2.1544$
EESPT		$812.801$	$815.754$	$754.667$		$2.3357$	$2.3441$	$2.1685$
\bottomrule

Table 3. Table 3: Influence of the choice of the cost function involved in the EET and the EESPT on the bounds of the local effectivity indices.

\toprule		Bounds of the local effectivity indices  ${\eta }_E$
\ms \msMethods		$J_0$	$J_1$	$J_2$
\midruleEET		$0.675$ - $7.507$	$0.581$ - $7.375$	$0.639$ - $7.305$
EESPT		$0.581$ - $7.375$	$0.570$ - $7.209$	$0.610$ - $7.151$
\bottomrule

Table 4. Table 4: Comparison of the error estimators given by the EET, the SPET, and the EESPT.

\topruleMethods	Estimate  $\theta$	Effectivity index  $\eta$	Normalized CPU time
\midruleEET	$0.9227$	$5.4671$	$1.000$
SPET	$0.5928$	$3.5127$	$5.368$
EESPT	$0.9119$	$5.4031$	$1.009$
\bottomrule

Table 5. Table 7: Comparison of the error estimators given by the EET, the SPET, and the EESPT.

\topruleMethods	Estimate  $\theta$	Effectivity index  $\eta$	Normalized CPU time
\midruleEET	$812.771$	$6.001$	$1.000$
SPET	$604.311$	$4.462$	$5.017$
EESPT	$787.983$	$5.818$	$1.007$
\bottomrule

Table 6. Table 8: Influence of the choice of the cost function involved in the EET and the EESPT on the quality of the error estimators.

\toprule		Estimate  $\theta$				Effectivity index  $\eta$
\ms      \msMethods		$J_0$	$J_1$	$J_2$		$J_0$	$J_1$	$J_2$
\midruleEET		$812.771$	$787.989$	$782.711$		$6.001$	$5.818$	$5.779$
EESPT		$787.983$	$815.111$	$809.267$		$5.818$	$6.019$	$5.975$
\bottomrule

Table 7. Table 9: Influence of the choice of the cost function involved in the EET and the EESPT on the bounds of the local effectivity indices.

\toprule		Bounds of the local effectivity indices  ${\eta }_E$
\ms \msMethods		$J_0$	$J_1$	$J_2$
\midruleEET		$1.163$ - $58.900$	$1.201$ - $51.275$	$1.095$ - $59.823$
EESPT		$1.201$ - $51.274$	$0.948$ - $57.349$	$0.878$ - $62.153$
\bottomrule

Table 8. Table 10: Comparison of the error estimators given by the EET, the SPET, and the EESPT.

\topruleMethods	Estimate  $\theta$	Effectivity index  $\eta$	Normalized CPU time
\midruleEET	$\mathrm{58\hspace{0.17em}061}$	$15.071$	$1.000$
SPET	$\mathrm{18\hspace{0.17em}948}$	$4.918$	$4.828$
EESPT	$\mathrm{42\hspace{0.17em}078}$	$10.922$	$1.007$
\bottomrule

Table 9. Table 11: Influence of the choice of the cost function involved in the EET and the EESPT on the quality of the error estimators.

\toprule		Estimate  $\theta$				Effectivity index  $\eta$
\ms      \msMethods		$J_0$	$J_1$	$J_2$		$J_0$	$J_1$	$J_2$
\midruleEET		$\mathrm{58\hspace{0.17em}061}$	$\mathrm{42\hspace{0.17em}078}$	$\mathrm{42\hspace{0.17em}995}$		$15.071$	$10.922$	$11.160$
EESPT		$\mathrm{42\hspace{0.17em}078}$	$\mathrm{42\hspace{0.17em}985}$	$\mathrm{43\hspace{0.17em}574}$		$10.922$	$11.158$	$11.311$
\bottomrule

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,$$

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

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

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

$${\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({\hat{\underline{u}}}_h,{\widehat{\text{dsigma}}}_h)$$

$${\int }_E\left({\widehat{\text{dsigma}}}_h-{\text{dsigma}}_h\right)\underline{\nabla }{\phi }_i\mathrm{d}\Omega =\underline{0}\forall E\in \mathcal{E},\forall i\in {\mathcal{N}}_E^{\mathcal{E}},$$

$${\widehat{\text{dsigma}}}_{h|E}{\underline{n}}_E={\eta }_E{\hat{\underline{F}}}_h\text{on}\partial E,$$

$${\eta }_E{\hat{\underline{F}}}_h={\underline{F}}_d\text{on}\partial E\subset {\partial }_2\Omega$$

$${\int }_{\partial E}{\eta }_E{\hat{\underline{F}}}_h\cdot {\underline{u}}_s\mathrm{d}S+{\int }_E{\underline{f}}_d\cdot {\underline{u}}_s\mathrm{d}\Omega =0\forall {\underline{u}}_s\in {\mathcal{U}}_{S|E},$$

$${\int }_{\partial E}{\eta }_E{\hat{\underline{F}}}_h{\phi }_i\mathrm{d}S={\int }_E\left({\text{dsigma}}_h\underline{\nabla }{\phi }_i-{\underline{f}}_d{\phi }_i\right)\mathrm{d}\Omega ={\underline{\mathcal{Q}}}_E^{(i)}\forall E\in \mathcal{E},\forall i\in {\mathcal{N}}_E^{\mathcal{E}},$$

$$\displaystyle\int_{\partial E}\eta_{E}\>\hat{\underline{F}}_{h}\cdot\underline{u}_{h}^{\ast}\mathop{}\mathopen{}\mathrm{d}S=\int_{E}\left(\operatorname{Tr}\big{[}\dsigma_{h}\>\depsil(\underline{u}_{h}^{\ast})\big{]}-\underline{f}_{d}\cdot\underline{u}_{h}^{\ast}\right)\mathop{}\mathopen{}\mathrm{d}\Omega\quad\forall\>\underline{u}_{h}^{\ast}\in\bm{\mathcal{U}}^{p}_{h\mathclose{}|\mathopen{}{E}},$$

$$\sum _{\Gamma \in {\mathcal{J}}_i^{\mathcal{N}}\cap \partial E}{\eta }_E{\hat{\underline{b}}}_{|\Gamma }^{(i)}={\underline{\mathcal{Q}}}_E^{(i)}\forall E\in {\mathcal{E}}_i^{\mathcal{N}},$$

$${\hat{\underline{b}}}_{|\Gamma }^{(i)}={\int }_{\Gamma }{\hat{\underline{F}}}_h{\phi }_i\mathrm{d}S.$$

$$\sum _{E\in {\mathcal{E}}_i^{\mathcal{N}}}{\underline{\mathcal{Q}}}_E^{(i)}=\underline{0}\text{for an internal node}i\in \Omega ,$$

$$\sum _{E\in {\mathcal{E}}_i^{\mathcal{N}}}{\underline{\mathcal{Q}}}_E^{(i)}={\int }_{\Gamma \in {\mathcal{J}}_i^{\mathcal{N}}\cap {\partial }_2\Omega }{\underline{F}}_d{\phi }_i\mathrm{d}S\text{for a node}i\subset {\partial }_2\Omega .$$

$$\sum _{E\in {\mathcal{E}}_i^{\mathcal{N}}}{\underline{\mathcal{Q}}}_E^{(i)}={\int }_{\Omega }\left({\text{dsigma}}_h\underline{\nabla }{\phi }_i-{\underline{f}}_d{\phi }_i\right)\mathrm{d}\Omega ,$$

$$J({\hat{\underline{b}}}^{(i)})=\frac{1}{2}\sum _{\Gamma \in {\mathcal{J}}_i^{\mathcal{N}}}({\hat{\underline{b}}}^{(i)}-{\underline{b}}^{(i)}{)}_{|\Gamma }^2,$$

$${\underline{b}}_{|\Gamma }^{(i)}=\frac{1}{2}{\int }_{\Gamma }(\frac{1}{{\eta }_E}{\text{dsigma}}_{h|E}{\underline{n}}_E+\frac{1}{{\eta }_{E^{\prime }}}{\text{dsigma}}_{h|E^{\prime }}{\underline{n}}_{E^{\prime }}){\phi }_i\mathrm{d}S;$$

$${\underline{b}}_{|\Gamma }^{(i)}={\int }_{\Gamma }\frac{1}{{\eta }_E}{\text{dsigma}}_{h|E}{\underline{n}}_E{\phi }_i\mathrm{d}S;$$

$${\underline{b}}_{|\Gamma }^{(i)}={\int }_{\Gamma }\frac{1}{{\eta }_E}{\underline{F}}_d{\phi }_i\mathrm{d}S.$$

$${\hat{\underline{F}}}_{h|\Gamma }=\sum _{j\in {\mathcal{N}}_{\Gamma }^{\mathcal{J}}}{\hat{\underline{F}}}_j^{\Gamma }{\phi }_{j|\Gamma }.$$

$${\widehat{\text{dsigma}}}_{h|E}\in {\mathcal{S}}_{{\hat{\underline{F}}}_h}⟺\{\begin{array}{rl}& {\widehat{\text{dsigma}}}_{h|E}\in \mathcal{S}\\ & \underline{\mathrm{div}}{\widehat{\text{dsigma}}}_{h|E}+{\underline{f}}_d=\underline{0}\text{in}E\\ & {\widehat{\text{dsigma}}}_{h|E}{\underline{n}}_E={\eta }_E{\hat{\underline{F}}}_h\text{on}\partial E\end{array}$$

$$\displaystyle\int_{\Omega}\operatorname{Tr}\big{[}\textbf{\em K}\>\depsil(\underline{e}_{h})\>\depsil(\underline{u}^{\ast})\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega$$

$$={\mathcal{R}}_h({\underline{u}}^{\ast })\forall {\underline{u}}^{\ast }\in {\mathcal{U}}_0,$$

$${\lambda }_i({\underline{x}}_j)={\delta }_{ij}\forall j\in \mathcal{I};\sum _{i\in \mathcal{I}}{\lambda }_i(\underline{x})=1\text{for}\underline{x}\in \Omega .$$

$$\displaystyle\int_{\Omega}\operatorname{Tr}\big{[}\textbf{\em K}\>\depsil(\underline{e}_{h})\>\depsil(\underline{u}^{\ast})\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega=\sum_{i\in\mathcal{I}}\mathcal{R}_{h}(\lambda_{i}\>\underline{u}^{\ast})\quad\forall\>\underline{u}^{\ast}\in\bm{\mathcal{U}}_{0}.$$

$$\displaystyle\int_{\Omega_{i}}\operatorname{Tr}\big{[}\textbf{\em K}\>\depsil(\underline{e}_{i})\>\depsil(\underline{u}^{\ast})\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega=\mathcal{R}_{h}(\lambda_{i}\>\underline{u}^{\ast})\quad\forall\>\underline{u}^{\ast}\in\bm{\mathcal{U}}_{0\mathclose{}|\mathopen{}{\Omega_{i}}},$$

$${\mathcal{R}}_h({\underline{u}}_h^{\ast })=0\forall {\underline{u}}_h^{\ast }\in {\mathcal{U}}_{h,0}^p,$$

$$\displaystyle\int_{\Omega_{i}}\operatorname{Tr}\big{[}\textbf{\em K}\>\depsil(\underline{e}_{i})\>\depsil(\underline{u}^{\ast})\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega=\mathcal{R}_{h}(\lambda_{i}\>(\underline{u}^{\ast}-\Pi\underline{u}^{\ast}))\quad\forall\>\underline{u}^{\ast}\in\bm{\mathcal{U}}_{0\mathclose{}|\mathopen{}{\Omega_{i}}}.$$

$$\Pi {\underline{u}}_s={\underline{u}}_s\forall {\underline{u}}_s\in {\mathcal{U}}_{S,0|{\Omega }_i}⟹{\mathcal{R}}_h({\lambda }_i({\underline{u}}_s-\Pi {\underline{u}}_s))=0\forall {\underline{u}}_s\in {\mathcal{U}}_{S,0|{\Omega }_i},$$

Full text

The numerical results which are presented focus on industrial problems; they highlight the main advantages and drawbacks of the different methods and show that the behavior of the three estimators, which have the same convergence rate as the exact global error, is consistent. Two- and three-dimensional experiments have been carried out in order to compare the performance and the computational cost of the three different approaches. The analysis of the results reveals that the SPET is more accurate than EET and EESPT methods, but the corresponding computational cost is higher. Overall, the numerical tests prove the interest of the hybrid method EESPT and show that it is a correct compromise between quality of the error estimate, practical implementation and computational cost. Furthermore, the influence of the cost function involved in the EET and the EESPT is studied in order to optimize the estimators.

\NME

16002810

\runningheads

F. Pled, L. Chamoin and P. Ladevèze On techniques for constructing admissible stress fields in model verification

\corraddr

P. Ladevèze, LMT-Cachan (ENS-Cachan/CNRS/Paris 6 University), 61 Avenue du Président Wilson, 94235 CACHAN Cedex

1 INTRODUCTION

Assessing the global/goal-oriented discretization error, henceforth known as model verification, has become a major challenge and a topical issue for both industrial applications and academic research. The widespread availability of computer hardware and numerical tools has contributed to a recent upsurge in the development of virtual prototyping. The modeling of physical problems requires the use of an initial mathematical model, which is considered as the reference to build a discretized model whose calculation is performed by numerical methods suited to computing tools. One of the most widespread approximation methods is the Finite Element Method (FEM) providing a numerical approximated solution of the a priori unknown exact solution of the reference model.

Within the finite element framework, the first works dealing with verification date back from the late 1970s [1, 2, 3, 4] and provide a global estimation of the discretization error which allows the global quality of a FE calculation to be quantified [5, 6, 7, 8]. For over fifteen years, various methods had been introduced to control the numerical quality of specific quantities of interest, such as local stresses, displacement values; they lead to bounds on functional outputs which are relevant information for design purposes in Mechanics [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. However, among those global/goal-oriented error estimators built from different methods, few lead to guaranteed and relevant bounds, which is a serious drawback in the domain of robust design. Only methods based on the construction of statically admissible stress fields actually yield guaranteed bounds of the global discretization error and allow the treatment of a wide range of mechanical problems (visco-elasticity, visco-plasticity, transient dynamics, …).

Various techniques have been developed for the construction of such stress fields; the first technique, proposed by Fraeijs de Veubeke, is based on a dual formulation, which corresponds to the best approach as regards the quality of the error estimator [22, 23]. However, despite its significant performances, this method is very expensive and intractable in practice as it requires the calculation of another global approximated solution of the reference problem; therefore its implementation is not suited for current FE codes. Three other techniques are well-suited to error estimation as they provide a statically admissible stress field from the FE solution: the element equilibration technique (EET) [2, 24, 25, 26, 27]; the star-patch equilibration technique (SPET) [28, 29, 30, 31, 32, 33, 34, 35, 36]; the element equilibration + star-patch technique (EESPT) [37, 38].

First, the EET, also called standard technique, introduces equilibrated tractions along finite element edges, leading to a quasi-explicit calculation of admissible stress fields at the element scale. The construction of tractions, in equilibrium with the external loading, is based on properties of the finite element stress field through an energy relation called the prolongation condition. Despite its affordable computational cost, the procedure for defining such tractions is in general difficult to implement into existing finite element codes. On the contrary, the SPET, also called the flux-free technique, precludes the main drawbacks of the EET as it circumvents the need of flux-equilibration. The principle of the SPET is to exploit a partition of unity in order to define self-equilibrated local problems (it avoids the calculation of equilibrated tractions) and results in a simple implementation. However, this technique requires the fine solution of local problems on sets of elements, also called patches or stars, involving a large number of degrees of freedom (particularly for 3D applications) and, therefore, leading to the calculation of admissible stress fields at high computational cost. Eventually, the EESPT is a hybrid technique which takes advantage of the ingredients of both EET and SPET methods. Indeed, this technique combines the prolongation condition used in the EET and the partition of unity involved in the SPET in order to construct equilibrated tractions over element edges. This hybrid procedure comes down to solving simple, local, independent problems defined over patches of elements, then solving problems at the element scale at reasonable computational cost.

This paper aims at comparing these three techniques for constructing statically admissible stress fields in terms of quality of the associated error estimate, practical implementation and computational cost. The performance of error estimators is investigated by comparing the estimate with the energy norm of the exact error. The numerical experiments are carried out on industrial cases. Thus, this work is in line with previous studies even though it goes further than that presented in [38], which compared the three error estimation methods on two-dimensional academic examples, and provides more details on implementation. A good match is observed between the estimates provided by the three techniques and the energy norm of the reference error. Furthermore, different cost functions involved in the minimization step of both EET and EESPT methods are considered and analyzed to assess their influence on the corresponding error estimates.

The remainder of the paper is structured as follows: after this introduction,

2

7 defines the reference problem and its discretized form in order to introduce the basic concept of admissible solution;

3

8 deals with the main principles of the three existing methods for constructing admissible stress fields, while

4

12 is devoted to a detailed analysis of their practical implementation; several two- and three-dimensional numerical results are presented in

5

13 to investigate the advantages and drawbacks of those techniques on complex industrial structures; finally,

6

15 suggests several perspectives on which future research should focus.

7 BASICS ON ERROR ESTIMATION AND ADMISSIBLE SOLUTIONS

7.1 Presentation of the reference problem

Let us consider an open bounded domain $\Omega$. Its boundary $\partial\Omega$ is divided into two disjoints, complementary parts $\partial_{1}\Omega\,(\neq\emptyset)$ and $\partial_{2}\Omega$ such that $\overline{\partial_{1}\Omega\cup\partial_{2}\Omega}=\partial\Omega$, $\partial_{1}\Omega\cap\partial_{2}\Omega=\emptyset$. The problem is assumed to be linear under the assumptions of elastic material, small perturbations and the evolution in time is considered to be quasi-static and isothermal. The mechanical structure is subjected to a loading represented by (see Figure 1):

• •

a prescribed displacement $\underline{U}_{d}$ on part $\partial_{1}\Omega$;

• •

a traction force density $\underline{F}_{d}$ on part $\partial_{2}\Omega$;

• •

a body force field $\underline{f}_{d}$ in $\Omega$.

The reference problem consists of searching a displacement/stress pair $(\underline{u}(\underline{M}),\dsigma(M)),M\in\Omega$, which verifies:

• •

the kinematic conditions:

[TABLE]

• •

the equilibrium equations:

[TABLE]

• •

the constitutive relation:

[TABLE]

where $\bm{\mathcal{U}}=[\mathcal{H}^{1}(\Omega)]^{3}$ and $\bm{\mathcal{S}}=\left\{\dsigma/\dsigma=\dsigma^{T},\dsigma\in[\mathcal{L}^{2}(\Omega)]^{6}\right\}$ are the spaces of finite-energy solutions; $\mathcal{H}^{1}(\Omega)$ stands for the standard Sobolev space of square integrable functions and first derivatives, while $\mathcal{L}^{2}(\Omega)$ refers to the space of square integrable functions. $\bm{\mathcal{U}}_{0}=\left\{\underline{u}^{\ast}\in\bm{\mathcal{U}},\underline{u}^{\ast}_{\mathclose{}|\mathopen{}{\partial_{1}\Omega}}=\underline{0}\right\}$ represents the vectorial space associated to $\bm{\mathcal{U}}$, and K is the Hooke elasticity tensor.

In the following, the exact solution of the reference problem is designated by $(\underline{u},\dsigma)$.

7.2 Presentation of the finite element problem

In most practical cases, the exact solution $(\underline{u},\dsigma)$ cannot be obtained explicitly. Typically, numerical methods are employed to achieve an approximation of $(\underline{u},\dsigma)$. One of the most widespread approximation methods is the finite element method (FEM) which provides a numerical solution $(\underline{u}_{h},\dsigma_{h})$ lying in the discretized space $\bm{\mathcal{U}}_{h}\!\times\!\bm{\mathcal{S}}_{h}\subset\bm{\mathcal{U}}\!\times\!\bm{\mathcal{S}}$. This is defined from piecewise continuous polynomial displacement shape functions associated with a spatial discretization (finite element space mesh $\mathcal{M}_{h}$) of the domain $\Omega$. The prescribed displacement field $\underline{U}_{d}$ is assumed to be compatible with the FE discretization. Thus, the finite element problem consists of finding a displacement/stress pair $(\underline{u}_{h}(\underline{M}),\dsigma_{h}(\underline{M})),\underline{M}\in\Omega$, which verifies:

• •

the kinematic conditions:

[TABLE]

• •

the equilibrium equations:

[TABLE]

• •

the constitutive relation:

[TABLE]

where $\bm{\mathcal{U}}_{h,0}=\bm{\mathcal{U}}_{h}\cap\bm{\mathcal{U}}_{0}$.

In the displacement-type finite element framework, the FE solution $(\underline{u}_{h},\dsigma_{h})$ satisfies both kinematic conditions (1) and constitutive relation (3) of the reference problem, but not equilibrium equations (2).

On the one hand, let us define the discretization error $\underline{e}_{h}=\underline{u}-\underline{u}_{h}$, also called the exact error or true error, corresponding to the difference between the exact displacement solution and the FE one; this error enables to control and assess the numerical quality of the FE solution $(\underline{u}_{h},\dsigma_{h})$. Usually, it is measured in terms of a suitable norm, such as the 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}, which leads to a global discretization error $\left\lVert\underline{e}_{h}\right\rVert_{u,\Omega}$. On the order hand, local errors $\Delta I=I(\underline{u})-I(\underline{u}_{h})$ can be defined if one seeks to evaluate quantities of interest, i.e. functional outputs $I(\underline{u})$ of the solution.

7.3 Construction of an admissible solution

In the framework of error estimation, a major interest concerns methods leading to strict and robust (global or local) error bounds for various types of mechanical problems. For this purpose, the construction of admissible solutions is required to obtain such reliable bounds [8, 2]. Furthermore, the quality of the yielded admissible solution can have a significant impact on the accuracy of the corresponding error estimate. An admissible pair, denoted $(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$, should verify (1) and (2), which amounts to saying that $\hat{\underline{u}}_{h}$ and $\hat{\dsigma}_{h}$ are kinematically and statically admissible, respectively. The strategy to set up an admissible solution of the reference problem relies on the post-processing of the data of both reference problem and FE solution $(\underline{u}_{h},\dsigma_{h})$ at hand. When using a displacement-type FE method, a natural way to reconstruct $\hat{\underline{u}}_{h}$ is to exploit straightforwardly the FE displacement field $\underline{u}_{h}$ as this one already satisfies the kinematic conditions. Thus, one generally chooses $\hat{\underline{u}}_{h}=\underline{u}_{h}$ except in the case of incompressible materials, for which a special method is employed [39]. Conversely, the reconstruction of an admissible stress field $\hat{\dsigma}_{h}$ is more complicated as the FE stress field $\dsigma_{h}$ does not verify the equilibrium equations (it does on the FE sense only). This stage represents the key difficulty of the construction of an admissible pair. Various techniques enable to construct an admissible stress field, the best one with respect to the quality of the error estimator being based on a dual formulation of the reference problem. However, it relies on the FE calculation of two completely distinct formulations, which represents a serious drawback in terms of computing time and implementation [22, 23]. Three concurrent and well-suited techniques will be described and discussed in details in the next section.

Eventually, an appraisal of the quality of the new approximate solution $(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ is provided by the measure $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}$ of the error on the constitutive relation (3), 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}. The constitutive relation error (CRE) is linked to the global discretization error through the well-known Prager-Synge theorem [40] implying that the constitutive relation error is a guaranteed upper bound of the global discretization error:

[TABLE]

Therefore, the so-called constitutive relation error $e_{\mathrm{cre}}(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$ represents an error estimate of the global discretization error $\left\lVert\underline{e}_{h}\right\rVert_{u,\Omega}$; it can be used for goal-oriented error estimation as well [13, 41, 19, 18, 20, 21]. Now, let us recall the basic ideas of the different techniques for constructing admissible stress fields in order to highlight the pros and cons of each method. These ideas have been more deeply presented in [38], particularly as regards the EESPT.

8 PRINCIPLES OF THE DIFFERENT TECHNIQUES FOR CONSTRUCTING ADMISSIBLE STRESS FIELDS

8.1 Notations

Let us define $\mathcal{E}$, $\mathcal{N}$, $\mathcal{I}$ and $\mathcal{J}$ the set of elements, nodes, vertices and edges of the FE mesh $\mathcal{M}_{h}$, respectively. $\mathcal{E}_{i}^{\mathcal{N}}\subset\mathcal{E}$, $\mathcal{E}_{i}^{\mathcal{I}}\subset\mathcal{E}$ and $\mathcal{E}_{\Gamma}^{\mathcal{J}}\subset\mathcal{E}$ stand for the set of elements connected to node $i$, vertex $i$ and edge $\Gamma$, respectively. $\mathcal{J}_{i}^{\mathcal{N}}\subset\mathcal{J}$ and $\mathcal{J}_{i}^{\mathcal{I}}\subset\mathcal{J}$ represent the set of edges connected to node $i$ and vertex $i$, respectively. $\mathcal{N}_{E}^{\mathcal{E}}\subset\mathcal{N}$ and $\mathcal{N}_{\Gamma}^{\mathcal{J}}\subset\mathcal{N}$ denote the set of nodes associated with element $E$ and edge $\Gamma$, respectively. $\mathcal{I}_{E}^{\mathcal{E}}\subset\mathcal{I}$ and $\mathcal{I}_{\Gamma}^{\mathcal{J}}\subset\mathcal{I}$ represent the set of vertices connected to element $E$ and edge $\Gamma$. Moreover, $\underline{x}_{i}$ designates the position of vertex $i$ in the FE mesh $\mathcal{M}_{h}$. Eventually, it is assumed that the FE displacement field $\underline{u}_{h}$ belongs to $\bm{\mathcal{U}}^{p}_{h}$, where $\bm{\mathcal{U}}^{p}_{h}$ corresponds to the FE interpolation space of maximum degree $p$. $\mathcal{U}^{p}_{h}$ denotes its one-dimensional correspondent.

8.2 The element equilibration technique (EET) - standard method

8.2.1 Principle of the construction

The first technique, called the element equilibration technique (EET), was introduced by Ladevèze [1, 2]. The principle is to exploit the FE properties of the stress field $\dsigma_{h}$ through an energy relation, called prolongation condition, which is the key point of the method. That relation links the searched admissible stress field $\hat{\dsigma}_{h}$ to the FE stress field $\dsigma_{h}$ under the form:

[TABLE]

where $\varphi_{i}\in\mathcal{U}^{p}_{h}$ stands for the FE shape function associated with node $i$. This prolongation condition is a physically sound relation as it imposes that the unknown admissible stress field $\hat{\dsigma}_{h}$ provides the same work as the FE stress field $\dsigma_{h}$ for each element $E$ of the FE mesh and for all FE displacement fields.

Classically, the EET is a quasi-explicit technique with a two-stage procedure:

1. (i)

construction of tractions $\hat{\underline{F}}_{h}$ in equilibrium with the external loading $(\underline{F}_{d},\underline{f}_{d})$ on element boundaries $\partial E$ of the spatial mesh $\mathcal{M}_{h}$; 2. (ii)

calculation of an admissible stress field $\hat{\dsigma}_{h}$ in equilibrium with these equilibrated tractions $\hat{\underline{F}}_{h}$ and body force field $\underline{f}_{d}$ at the element level.

8.2.2 First stage: construction of tractions $\hat{\underline{F}}_{h}$

Tractions $\hat{\underline{F}}_{h}$ are intended to represent the stress vectors $\hat{\dsigma}_{h\mathclose{}|\mathopen{}{E}}\>\underline{n}_{E}$ on edges $\partial E$ of element $E\in\mathcal{E}$:

[TABLE]

where $\underline{n}_{E}$ is the unit outward normal vector to element $E$ and $\eta_{E}=\pm 1$ are functions ensuring the continuity of the stress vector across $\partial E$.

Furthermore, these tractions $\hat{\underline{F}}_{h}$ are built in equilibrium with the external loading $(\underline{F}_{d},\underline{f}_{d})$:

[TABLE]

where $\bm{\mathcal{U}}_{S\mathclose{}|\mathopen{}{E}}$ is the set of rigid body displacement fields over element $E$.

The procedure for calculating tractions $\hat{\underline{F}}_{h}$ on the element boundaries of the spatial mesh $\mathcal{M}_{h}$ is quasi-explicit. First, tractions $\hat{\underline{F}}_{h}$ are set equal to the traction force density $\underline{F}_{d}$ on the subdomain $\partial_{2}\Omega$: $\eta_{E}\>\hat{\underline{F}}_{h}=\underline{F}_{d}$ on $\partial E\subset\partial_{2}\Omega$. Second, by applying integration by parts (Green’s theorem) before using definition (9) of tractions $\hat{\underline{F}}_{h}$ as well as the strong form of the problem associated to equilibrium equations (2) satisfied by $\hat{\dsigma}_{h}$, the prolongation condition (8) can be recasted in the following form:

[TABLE]

Remark 1

The derivation of (12) also assumes some regularity for $\hat{\dsigma}_{h\mathclose{}|\mathopen{}{E}}$ over each element $E\in\mathcal{E}$: either $\hat{\dsigma}_{h\mathclose{}|\mathopen{}{E}}$ is continuous inside $E$ or, if a discontinuity occurs along a path $\Gamma\in E$, the normal component $\hat{\dsigma}_{h\mathclose{}|\mathopen{}{\Gamma}}\>\underline{n}_{\Gamma}$ along $\Gamma$ is continuous.

It is worthy noticing that (12) enforces the equilibrium condition (11) of tractions $\hat{\underline{F}}_{h}$ with the body force field $\underline{f}_{d}$ over each element $E\in\mathcal{E}$, as (12) is equivalent to

[TABLE]

and $\bm{\mathcal{U}}_{S\mathclose{}|\mathopen{}{E}}\subset\bm{\mathcal{U}}^{p}_{h\mathclose{}|\mathopen{}{E}}$.

Thus, the prolongation condition (12) defines local problems $\mathcal{P}_{i}^{\mathcal{N}}$ associated with each node $i\in\mathcal{N}$. Problem $\mathcal{P}_{i}^{\mathcal{N}}$ associated with node $i\in\mathcal{N}$ is a linear system that reads:

[TABLE]

where the unknown quantity $\displaystyle\hat{\underline{b}}_{\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ over the edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{N}}$ is the projection of traction $\hat{\underline{F}}_{h\mathclose{}|\mathopen{}{\Gamma}}$ over the FE shape function $\varphi_{i}$:

[TABLE]

Solvability of local problems $\mathcal{P}_{i}^{\mathcal{N}}$ (14) depends on the type of node $i$ considered. Indeed, existence of a solution for system (14) most of the time requires the verification of a compatibility condition, that reads:

[TABLE]

Noticing that:

[TABLE]

the compatibility conditions (16) and (17) are ensured by the FE equilibrium equations (5) satisfied by the FE stress field $\dsigma_{h}$, as the vectorial equation (18) is linked after some straightforward computations to the scalar equilibrium equations (5).

If local problems $\mathcal{P}_{i}^{\mathcal{N}}$ (14) can admit several solutions, the uniqueness of the solution is ensured by performing the least-squares minimization of a cost function of the form [8]:

[TABLE]

where the known quantity $\displaystyle\underline{b}_{\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ represents the projection of the FE stress vector $\dsigma_{h}\>\underline{n}$, over the edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{N}}$, on the FE shape function $\varphi_{i}$:

• •

for an internal edge $\Gamma$ connected to elements $E$ and $E^{\prime}$,

[TABLE]

• •

for an edge $\Gamma\subset\partial_{1}\Omega$ connected to element $E$,

[TABLE]

• •

for an edge $\Gamma\subset\partial_{2}\Omega$ connected to element $E$,

[TABLE]

Other cost functions will be considered in

9

13 to study the influence of its choice on the quality of the error estimate.

Eventually, densities $\hat{\underline{F}}_{h}$ along edges $\Gamma\in\mathcal{J}$ are recovered from $\displaystyle\hat{\underline{b}}^{(i)}$, using the same interpolation degree as the FE displacement field $\underline{u}_{h}$. In other words, one seeks $\displaystyle\hat{\underline{F}}_{h\mathclose{}|\mathopen{}{\Gamma}}\in\bm{\mathcal{U}}^{p}_{h\mathclose{}|\mathopen{}{\Gamma}}$ under the form:

[TABLE]

Therefore, the evaluation of tractions $\displaystyle\hat{\underline{F}}_{h\mathclose{}|\mathopen{}{\Gamma}}$ along each edge $\Gamma\in\mathcal{J}$ requires the solution of a set of linear local problems at the edge level involving the projections $\displaystyle\hat{\underline{b}}_{\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ associated to each node $i\in\mathcal{N}_{\Gamma}^{\mathcal{J}}$.

9.0.1 Second stage: calculation of an admissible stress field $\hat{\dsigma}_{h}$ at the element scale

Once the flux-equilibration procedure has been performed, one seeks the local restriction $\hat{\dsigma}_{h\mathclose{}|\mathopen{}{E}}$ of an admissible stress field $\hat{\dsigma}_{h}$ to each element $E\in\mathcal{E}$ as the solution of an equilibrium local problem $\mathcal{P}_{E}^{\mathcal{E}}$ where the previously calculated tractions $\displaystyle\hat{\underline{F}}_{h}$ play the role of the external loading over the element boundaries:

[TABLE]

It is important to recall that the construction of tractions $\displaystyle\hat{\underline{F}}_{h}$ leads to equilibrated local problems $\mathcal{P}_{E}^{\mathcal{E}}$ (21) ensuring their solvability, since $\hat{\dsigma}_{h\mathclose{}|\mathopen{}{E}}$ is assumed to have no discontinuity inside each element $E\in\mathcal{E}$. In practice, a dual formulation (displacement-type FEM) of local static problems $\mathcal{P}_{E}^{\mathcal{E}}$ (21) combined with a $p$-refinement or $h$-refinement of the FEM is used to obtain a numerical approximation of $\hat{\dsigma}_{h\mathclose{}|\mathopen{}{E}}$. In other words, the discretization of each element $E\in\mathcal{E}$ consists of considering either a single element $E$ along with an interpolation of degree $p+k$, where $p$ denotes the degree of the FE analysis and $k$ an additional degree, or a subdivision of element $E$ along with an interpolation of degree $p$. The $p$-refinement technique, introduced by Babus̆ka and Strouboulis [42], has a higher convergence rate than the $h$-refinement one and circumvents the need to generate a new refined mesh. Using an extra-degree $k=3$ yields a very good and accurate approximation of a statically admissible stress field, thus leading to an error estimate of good quality [43]. The $p$-refinement technique could also be combined with a $h$-refinement technique. The interested reader can refer to [8] for more information.

Remark 2

The stress field $\hat{\dsigma}_{h\mathclose{}|\mathopen{}{E}}$ obtained using a $p$-refinement technique is not strictly admissible in the sense that local problems $\mathcal{P}_{E}^{\mathcal{E}}$ (21) are not solved exactly; such a stress field leads to only asymptotic error bounds with respect to a refined solution.

Remark 3

Local problems $\mathcal{P}_{E}^{\mathcal{E}}$ (21) can however be solved exactly: a first procedure consists of constructing an admissible stress field analytically; this construction requires a subsplitting of elements, the statically admissible stress field being searched using a piecewise polynomial basis [8, 25]. This analytical construction yields strictly admissible stress fields provided that body force field $\underline{f}_{d}$ is a polynomial with a degree compatible with that of $\hat{\dsigma}_{h\mathclose{}|\mathopen{}{E}}$. Another procedure consists of using a stress-type FEM. Both procedures yield admissible stress fields and therefore strict error bounds with respect to the exact solution.

The EET has a very attractive feature, namely its affordable computational cost, as local problems $\mathcal{P}_{i}^{\mathcal{N}}$ (14) and $\mathcal{P}_{E}^{\mathcal{E}}$ (21) are both reasonably costly. Despite this advantage, the flux-equilibration procedure needed to properly set boundary conditions for local problems $\mathcal{P}_{E}^{\mathcal{E}}$ (21) is complex to implement, which strongly limits its practical use into existing FE codes. Nevertheless, it has already been implemented into the industrial software SAMCEF [44].

9.1 The star-patch equilibration technique (SPET) - flux-free method

The second method, called the star-patch equilibration technique (SPET), was formerly developed for fluid mechanics [29, 30, 31] and has been revisited and adapted to solid mechanics by Parés, Díez and Huerta [32] under the name flux-free technique. This technique does not require any flux recovery or flux splitting technique as there is no need to perform flux equilibration; this constitutes a serious advantage for implementation. The novelty of this technique is the introduction of a partition of unity which allows to boil down to self-equilibrated local problems defined on sets of elements, also called patches or stars.

9.1.1 Definition of self-equilibrated local problems over patches of elements

The global problem defining the discretization error $\underline{e}_{h}=\underline{u}-\underline{u}_{h}$ can be obtained from the equilibrium equations (2) replacing $\underline{u}$ by $\underline{e}_{h}+\underline{u}_{h}$. It reads:

Find $\underline{e}_{h}\in\bm{\mathcal{U}}_{0}$ such that:

[TABLE]

where $\mathcal{R}_{h}$ stands for the weak residual functional associated with the FE approximation $\underline{u}_{h}$. This residual equation conveys the failure to comply with the FE equilibrium equations (5).

The main idea consists of introducing at this stage the partition of unity defined by the linear FE shape functions $\lambda_{i}\in\mathcal{U}_{h}^{1}$ based on vertices $i\in\mathcal{I}$, which are such that:

[TABLE]

Global problem (22) can thus be reformulated as:

[TABLE]

Let us define the patch $\Omega_{i}$ as the support of $\lambda_{i}$; in other words, $\Omega_{i}$ designates the support of the set $\mathcal{E}_{i}^{\mathcal{I}}$ of elements connected to vertex $i\in\mathcal{I}$ (see Figure 2).

Therefore, the SPET consists of solving a set of local problems in each patch $\Omega_{i}$:

Find $\underline{e}_{i}\in\bm{\mathcal{U}}_{0\mathclose{}|\mathopen{}{\Omega_{i}}}$ such that:

[TABLE]

where $\bm{\mathcal{U}}_{0\mathclose{}|\mathopen{}{\Omega_{i}}}$ is the restriction of $\bm{\mathcal{U}}_{0}$ to patch $\Omega_{i}$: $\bm{\mathcal{U}}_{0\mathclose{}|\mathopen{}{\Omega_{i}}}=\left\{\underline{u}^{\ast}\in\bm{\mathcal{U}}_{\mathclose{}|\mathopen{}{\Omega_{i}}},\underline{u}^{\ast}_{\mathclose{}|\mathopen{}{\partial_{1}\Omega\cap\Omega_{i}}}=\underline{0}\right\}$.

Now, let us recall the key property of the Galerkin approach, which is the Galerkin orthogonality:

[TABLE]

where $\bm{\mathcal{U}}_{h,0}^{p}=\bm{\mathcal{U}}_{h}^{p}\cap\bm{\mathcal{U}}_{0}$.

Let us note that the weak residual $\mathcal{R}_{h}(\lambda_{i}\>\bullet)$ does not allow to ensure the solvability of local problems (24). Indeed, existence of a solution of (24) is guaranteed if and only if the following compatibility condition holds: $\mathcal{R}_{h}(\lambda_{i}\>\underline{u}_{s})=0\quad\forall\>\underline{u}_{s}\in\bm{\mathcal{U}}_{S,0\mathclose{}|\mathopen{}{\Omega_{i}}}$, where $\bm{\mathcal{U}}_{S,0\mathclose{}|\mathopen{}{\Omega_{i}}}$ is the restriction of $\bm{\mathcal{U}}_{S,0}$ to patch $\Omega_{i}$, i.e. the set of rigid body displacement fields defined over $\Omega_{i}$ which vanish on $\partial_{1}\Omega\cap\Omega_{i}$. However, this condition is not verified for an interpolation degree $p=1$, as $\lambda_{i}\>\underline{u}_{s}\in\bm{\mathcal{U}}_{h,0}^{2}\not\subset\bm{\mathcal{U}}_{h,0}^{1}$.

In [32], the projector $\Pi\colon\bm{\mathcal{U}}_{0}\to\bm{\mathcal{U}}_{h,0}^{1}$ whose image is in the space of linear FE shape functions was defined. By using Galerkin orthogonality (25), the right-hand side term of (24) can be replaced by $\mathcal{R}_{h}(\lambda_{i}\>(\underline{u}^{\ast}-\Pi\underline{u}^{\ast}))$ leading to a new set of local problems $\mathcal{P}_{i}^{\mathcal{I}}$ over each patch $\Omega_{i}$ in the form:

Find $\underline{e}_{i}\in\bm{\mathcal{U}}_{0\mathclose{}|\mathopen{}{\Omega_{i}}}$ such that:

[TABLE]

By observing that:

[TABLE]

local problems $\mathcal{P}_{i}^{\mathcal{I}}$ (26) are well-posed and solvable, the projector being the key ingredient to ensure their solvability, especially for an interpolation degree $p=1$. Let us note that for scalar (i.e. one-dimensional) problems and for two- or three-dimensional mechanical problems with high-order interpolation degree (at least quadratic), it is not necessary to introduce the projector $\Pi$ into the right-hand side term of (24) to achieve the solvability of local problems $\mathcal{P}_{i}^{\mathcal{I}}$ (26).

Therefore, local problems (26) defined at the patch scale are self-equilibrated so that their solution does not require the construction of tractions or flux jumps along element boundaries and there is no need to perform any flux equilibration to achieve equilibrium, which is a very interesting feature for implementation purposes. In practice, the solution of local problems (26) is computed using a $p$-refinement or $h$-refinement over patches $\Omega_{i}$. The calculation is classically performed using the original FE mesh $\mathcal{M}_{h}$ with a $p+3$ discretization over each patch $\Omega_{i}$ in order to obtain a fairly good approximation of the solution of (26).

9.1.2 Construction of both an admissible stress field and a global estimate

Summing (26) for all vertices $i\in\mathcal{I}$ leads to a relation between the discretization error $\underline{e}_{h}$ and numerical solutions $\underline{e}_{i}$:

[TABLE]

Consequently, numerical solutions $\underline{e}_{i}$ obtained from (26) allow to define a global error estimate, which is an overestimation of the energy norm of the discretization error:

[TABLE]

with \displaystyle\left\lVert\bullet\right\rVert_{u,E}=\left(\int_{E}\operatorname{Tr}\big{[}\textbf{\em K}\>\depsil(\bullet_{\mathclose{}|\mathopen{}{E}})\>\depsil(\bullet_{\mathclose{}|\mathopen{}{E}})\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega\right)^{1/2}.

Equivalently, one can define an admissible stress field over each element $E\in\mathcal{E}$:

[TABLE]

Finally, the SPET evades the main disadvantage of the EET, because it does not require the construction of equilibrated tractions making its implementation less cumbersome. On the other hand, local problems $\mathcal{P}_{i}^{\mathcal{I}}$ (26) are defined at the patch scale and, therefore, involve a large number of degrees of freedom, particularly in 3D, compared with local problems $\mathcal{P}_{E}^{\mathcal{E}}$ (21) defined at the element scale for the EET. Thus, the calculation of an admissible stress field using the SPET may result in simply defined but costly computations.

Remark 4

Local problems $\mathcal{P}_{i}^{\mathcal{I}}$ (26) defined on patches can be solved using a dual approach, that is a stress-type FEM [35].

9.2 The element equilibration and star-patch technique (EESPT) - hybrid method

9.2.1 Principle of the construction

This new technique, based on recent works of Ladevèze et al [38], is a hybrid method insofar as it takes advantage of the ingredients of both EET and SPET methods. As for the EET, the procedure to construct an admissible stress field involves two main stages:

1. (i)

construction of tractions $\hat{\underline{F}}_{h}$ in equilibrium with the external loading $(\underline{F}_{d},\underline{f}_{d})$ on element boundaries $\partial E$ of the spatial mesh $\mathcal{M}_{h}$; 2. (ii)

construction of an admissible stress field $\hat{\dsigma}_{h}$ solution of the static local problem $\mathcal{P}_{E}^{\mathcal{E}}$ (21) over each element $E\in\mathcal{E}$.

The second stage is similar to that of the EET. As previously explained, those local problems $\mathcal{P}_{E}^{\mathcal{E}}$ are solved numerically by using a dual formulation of (21) and the FEM through a refined mesh ($p$-refinement); one can refer to

10

9.0.1 for details. The main thrust relates to the first stage, which aims at constructing equilibrated tractions, in order to facilitate its practical implementation. The principle is to exploit not only the prolongation condition used in the EET, but also the partition of unity involved in the SPET.

10.0.1 Definition of a set of local problems over patches of elements

Prolongation condition (8), in which the FE shape functions $\varphi_{i}$ belong to $\mathcal{U}^{p}_{h}$, can be reformulated in the global form:

[TABLE]

where $\bm{\mathcal{V}}_{h}^{p}$ designates the space of polynomial functions which are continuous over each element $E\in\mathcal{E}$ but possibly discontinuous across the inter-element edges.

One can replace the vectorial space $\bm{\mathcal{V}}_{h}^{p}$ of degree $p$ by the vectorial space $\bm{\mathcal{V}}_{h}^{q}$ of degree $q$ such that $1\leqslant q\leqslant p$, since this one suffices to ensure the equilibrium condition (11). Let us note that using $q<p$ is not equivalent to enforcing prolongation condition (8). For the sake of simplicity and affordable computational cost, one limits to the case $q=1$, but the method can be generalized for functions $\underline{v}_{h}^{\ast}\in\bm{\mathcal{V}}_{h}^{q}$ of degree $q$ with $1\leqslant q\leqslant p$.

By considering linear FE shape functions $\varphi_{i}\in\mathcal{U}_{h}^{1}$ only, (31) becomes:

[TABLE]

Using the weak equilibrium between $\hat{\dsigma}_{h}$, $\hat{\underline{F}}_{h}$ and $\underline{f}_{d}$ gives:

[TABLE]

Then, incorporating the partition of unity defined by the linear FE shape functions $\lambda_{i}\in\mathcal{U}_{h}^{1}$ into (33) leads to:

[TABLE]

In order to offer better flexibility, we will consider the following set of local problems $\mathcal{P}_{i}^{\mathcal{I}}$ introduced in [38] and defined over each patch $\Omega_{i}$ by observing that $\Omega_{i}$ is the support of $\lambda_{i}$:

Find $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ such that:

[TABLE]

where $\displaystyle\mathcal{Q}_{\Omega_{i}}(\lambda_{i}\>\underline{v}_{h}^{\ast})=\int_{\Omega_{i}}\left(\operatorname{Tr}\left[\dsigma_{h}\>\depsil(\lambda_{i}\>\underline{v}_{h}^{\ast})\right]-\underline{f}_{d}\cdot\lambda_{i}\>\underline{v}_{h}^{\ast}\right)\mathop{}\mathopen{}\mathrm{d}\Omega$.

Let us note that the sum of local problems (35) (or (34)) for all vertices $i\in\mathcal{I}$ leads to equation (33), which shows that global equilibrium is ensured.

By noticing that solutions $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ of (35) are nonzero only along the edges connected to vertex $i$, i.e. $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$, local problems $\mathcal{P}_{i}^{\mathcal{I}}$ (35) can be rewritten as:

Find $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ such that:

[TABLE]

Both sets $\mathcal{J}_{i}^{\mathcal{I}}$ and $\mathcal{E}_{i}^{\mathcal{I}}$ are represented in Figure 3.

Consequently, once local problems (36) have been solved, tractions $\hat{\underline{F}}_{h}$ are recovered on each edge $\Gamma\in\mathcal{J}$ by summing the solutions $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ of local problems (36) for all the vertices connected to the edge $\Gamma$, i.e. for $i\in\mathcal{I}_{\Gamma}^{\mathcal{J}}$:

[TABLE]

In practice, local quantities $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ are sought along each edge $\Gamma\in\mathcal{J}$ in the space of piecewise continuous polynomial functions of degree $p$ defined on the element edges. Indeed, more precisely, $\displaystyle\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}\in\bm{\mathcal{U}}^{p}_{h\mathclose{}|\mathopen{}{\Gamma}}$ in order to recover $\hat{\underline{F}}_{h}$ with the same interpolation degree as the FE displacement field $\underline{u}_{h}\in\bm{\mathcal{U}}_{h}^{p}$.

Furthermore, let us recall that conditions $\eta_{E}\>\hat{\underline{F}}_{h}=\underline{F}_{d}$ have to be enforced over edges $\Gamma\subset\partial_{2}\Omega$.

10.0.2 Solvability and uniqueness of solution $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ for local problem $\mathcal{P}_{i}^{\mathcal{I}}$

First, let us address the question of solvability of local problems $\mathcal{P}_{i}^{\mathcal{I}}$. A detailed analysis of local problems $\mathcal{P}_{i}^{\mathcal{I}}$ over $\Omega_{i}$ is carried out in [38]. Let us recall the main points. By considering the set ${\bar{\bm{\mathcal{U}}}_{h,0\mathclose{}|\mathopen{}{\Omega_{i}}}^{1}}$ of polynomial functions $\underline{v}_{h}^{\ast}\in\bm{\mathcal{V}}_{h}^{1}$ which are continuous across edges $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$ and which vanish on $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}\cap\partial\Omega$, i.e. such that $\underline{v}_{h}^{\ast}=\underline{0}$ along $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}\cap\partial\Omega$, the existence of a solution of problem (36) is ensured for $p\geqslant 2$. However, in the case $p=1$, local problem $\mathcal{P}_{i}^{\mathcal{I}}$ over $\Omega_{i}$ is replaced by:

Find $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ such that:

[TABLE]

where $\displaystyle\mathcal{Q}_{\Omega_{i}}(\lambda_{i}\>(\underline{v}_{h}^{\ast}-\Pi^{i}\underline{v}_{h}^{\ast}))=\mathcal{Q}_{\Omega_{i}}(\lambda_{i}\>\underline{v}_{h}^{\ast}(\underline{x}_{i}))=\int_{\Omega_{i}}\left(\operatorname{Tr}\left[\dsigma_{h}\>\depsil(\lambda_{i}\>\underline{v}_{h}^{\ast}(\underline{x}_{i}))\right]-\underline{f}_{d}\cdot\lambda_{i}\>\underline{v}_{h}^{\ast}(\underline{x}_{i})\right)\mathop{}\mathopen{}\mathrm{d}\Omega$. Operator $\Pi^{i}\colon\bm{\mathcal{V}}_{h}^{1}\to\bm{\mathcal{V}}_{h,i}^{1}$ designates the projector onto the space $\bm{\mathcal{V}}_{h,i}^{1}$ of functions $\underline{v}_{h}^{\ast}\in\bm{\mathcal{V}}_{h}^{1}$ defined over $\Omega_{i}$ which vanish at vertex $i$. Mathematically, it reads $\Pi^{i}\underline{v}_{h\mathclose{}|\mathopen{}{E}}^{\ast}=\underline{v}_{h\mathclose{}|\mathopen{}{E}}^{\ast}-\underline{v}_{h\mathclose{}|\mathopen{}{E}}^{\ast}(\underline{x}_{i})\quad\forall\>E\in\mathcal{E}_{i}^{\mathcal{I}},\ \forall\>\underline{v}_{h}^{\ast}\in\bm{\mathcal{V}}_{h}^{1}$. The solvability of local problem $\mathcal{P}_{i}^{\mathcal{I}}$ defined over patch $\Omega_{i}$ by (38) is ensured and those local problems introduced for $p=1$ implicitly require that (see [38]):

[TABLE]

Remark 5

In the case $q=p$, $\underline{v}_{h}^{\ast}\in\bm{\mathcal{V}}_{h}^{p}$ and local problems (36) are not solvable. Therefore, the modification introduced in (38) should be used to ensure solvability.

Now, let us focus on the uniqueness of solution $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ of (36) (or (38)). This one is obtained by performing a least-squares minimization step over each patch $\Omega_{i}$ involving a cost function whose form is similar to (19):

[TABLE]

The known quantity $\lambda_{i}\>\underline{F}_{h\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ involves the projection $\Pi_{\dsigma_{h}}^{\Gamma}$ of the FE stress field $\dsigma_{h}$ over edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$:

• •

for an internal edge $\Gamma$ connected to elements $E$ and $E^{\prime}$,

[TABLE]

• •

for an edge $\Gamma\subset\partial_{1}\Omega$ connected to element $E$,

[TABLE]

• •

for an edge $\Gamma\subset\partial_{2}\Omega$ connected to element $E$,

[TABLE]

Then, $\lambda_{i}\>\underline{F}_{h\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ is defined for each edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$ as:

• •

for $p\geqslant 2$, $\lambda_{i}\>\underline{F}_{h\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ is the projection of $\lambda_{i}\>\Pi_{\dsigma_{h}}^{\Gamma}$ over the space of polynomial functions of degree $p$;

• •

for $p=1$, due to implicit conditions (39), $\lambda_{i}\>\underline{F}_{h\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ is such that:

[TABLE]

It is important to point out that, with respect to the expressions of $\lambda_{i}\>\underline{F}_{h}^{(i)}$, those quantities are sought along each edge $\Gamma\in\mathcal{J}$ in the same space as $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$, namely the space of piecewise continuous polynomial functions of degree $p$.

Let us recall that several cost functions will be defined and investigated in

11

13.

Computational aspects of the various methods and a detailed analysis of their implementation are addressed in the next section. This section particularly concerns the computation of the local tractions in the EET and EESPT methods.

12 PRACTICAL IMPLEMENTATION OF THE DIFFERENT METHODS

12.1 Practical implementation of the EET

Local problem $\mathcal{P}_{i}^{\mathcal{N}}$ (14) associated with a node $i\in\mathcal{N}$ takes the matrix form:

[TABLE]

where $\hat{\underline{b}}^{(i)}$ is the unknown generalized vector which contains the combination of vectors $\hat{\underline{b}}_{\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ for every edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{N}}$, and $\underline{Q}^{(i)}$ is the generalized vector which is the combination of vectors $\underline{\mathcal{Q}}_{E}^{(i)}$ for every element $E\in\mathcal{E}_{i}^{\mathcal{N}}$.

Let us recall that the solvability of (42), for both internal node $i\in\Omega$ and node $i\in\partial_{2}\Omega$, requires the verification of a compatibility condition resulting from the FE equilibrium. For such systems, the vectorial equation coming from one, and only one, element $E\in\mathcal{E}_{i}^{\mathcal{N}}$ is removed to fix the kernel of matrix $\mathbb{B}^{(i)}$. Thus, system (42) is changed to:

[TABLE]

Let $n_{ind}^{(i)}$ be the number of independent equations coming from the prolongation condition, i.e. the number of rows of matrix $\tilde{\mathbb{B}}^{(i)}$ involved in (43). Let $n_{unk}^{(i)}$ be the number of unknowns of system (43), i.e. the size of vector $\hat{\underline{b}}^{(i)}$.

Furthermore, for all nodes $i\in\partial_{2}\Omega$, conditions $\displaystyle\hat{\underline{b}}_{\mathclose{}|\mathopen{}{\Gamma}}^{(i)}=\int_{\Gamma}\underline{F}_{d}\>\varphi_{i}\mathop{}\mathopen{}\mathrm{d}S$ over each edge $\Gamma\subset\partial_{2}\Omega$ are enforced by Lagrange multipliers. Those additional conditions read:

[TABLE]

Let $n_{enf}^{(i)}$ be the number of enforced equations of system (44), i.e. the number of rows of matrix $\mathbb{C}^{(i)}$.

Besides, one performs a minimization step if and only if the following minimization condition holds: $n_{unk}^{(i)}>n_{enf}^{(i)}+n_{ind}^{(i)}$. In this case, cost function (19) can be rewritten mathematically as:

[TABLE]

where the generalized vector $\underline{b}^{(i)}$ is the combination of known vectors $\underline{b}_{\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ for every edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{N}}$ and $\mathbb{M}^{(i)}$ is a diagonal matrix.

Finally, the problem to be solved for each node $i\in\mathcal{N}$ depends on the minimization condition:

• •

if $n_{unk}^{(i)}>n_{enf}^{(i)}+n_{ind}^{(i)}$, it consists of minimizing cost function (45) under both constraints (44) enforced over $\partial_{2}\Omega$ and constraints (43) coming from prolongation condition (8). Thus, introducing the Lagrangian:

[TABLE]

the system to be solved takes the matrix form:

[TABLE]

where $\underline{\Lambda}_{\mathbb{C}}^{(i)}$ and $\underline{\Lambda}_{\mathbb{B}}^{(i)}$ correspond to the vectors of Lagrange multipliers associated with constraints (44) and (43), respectively;

• •

else, it consists of solving explicitly the following system:

[TABLE]

It is worthy noticing that, if $n_{unk}^{(i)}<n_{enf}^{(i)}+n_{ind}^{(i)}$, one gives priority to enforced constraints (44) involving traction force density $\underline{F}_{d}$. Thus, in this case, matrix $\tilde{\mathbb{B}}^{(i)}$ and vector $\tilde{\underline{Q}}^{(i)}$ are truncated so that the number of constraints coming from (43) (i.e. resulting from the prolongation condition) is reduced to $n_{unk}^{(i)}-n_{enf}^{(i)}$.

Over each edge $\Gamma\in\mathcal{J}$, traction $\hat{\underline{F}}_{h\mathclose{}|\mathopen{}{\Gamma}}$ reads:

[TABLE]

where $\hat{\underline{F}}^{\Gamma}$ is the unknown vector of components of $\hat{\underline{F}}_{h\mathclose{}|\mathopen{}{\Gamma}}$ over $\Gamma$, and $[\varphi_{\mathclose{}|\mathopen{}{\Gamma}}]^{T}$ is the vector of FE shape functions of degree $p$ over $\Gamma$.

Finally, once projections $\hat{\underline{b}}^{(i)}$ have been calculated, tractions $\hat{\underline{F}}_{h\mathclose{}|\mathopen{}{\Gamma}}$ are recovered by solving a linear local system over each edge $\Gamma\in\mathcal{J}$ in the form:

[TABLE]

where the generalized vector $\hat{\underline{b}}^{\Gamma}$ is the combination of vectors $\hat{\underline{b}}_{\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ for each node $i\in\mathcal{N}_{\Gamma}^{\mathcal{J}}$, and $\mathbb{K}^{\Gamma}$ is similar to a mass matrix associated with edge $\Gamma\in\mathcal{J}$.

12.2 Practical implementation of the SPET

In order to simply remove the first (linear) part $\Pi\underline{u}^{\ast}$ of $\underline{u}^{\ast}$ involved in the r.h.s. term of (26), quantities $\underline{e}_{i}$ and $\underline{u}^{\ast}$ are assumed to belong to the FE interpolation space defined from hierarchical shape functions of degree $p+3$ over each element $E\in\mathcal{E}_{i}^{\mathcal{I}}$, instead of Lagrange shape functions. Furthermore, solutions $\underline{e}_{i}$ of (26) have to verify the following kinematic conditions $\underline{e}_{i\mathclose{}|\mathopen{}{\partial_{1}\Omega\cap\Omega_{i}}}=\underline{0}$ (i.e. $\underline{e}_{i}\in\bm{\mathcal{U}}_{0\mathclose{}|\mathopen{}{\Omega_{i}}}$); those kinematic constraints are enforced by substitution in (26). Therefore, the procedure used for the calculation of solutions $\underline{e}_{i}$ of (26) is fairly simple to implement as the use of FE hierarchical shape functions enables to retain only the high-order part of $\underline{u}^{\ast}$, which apparently seems to be the main difficulty associated with this method.

12.3 Practical implementation of the EESPT

Searched quantity $\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ and known quantity $\lambda_{i}\>\underline{F}_{h}^{(i)}$ are discretized over edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$ in the form:

[TABLE]

where $[\varphi_{\mathclose{}|\mathopen{}{\Gamma}}]$ is the matrix of FE shape functions of degree $p$ over $\Gamma$ and $\hat{\underline{f}}_{h,\Gamma}^{(i)}$ (respectively $\underline{f}_{h,\Gamma}^{(i)}$) is the vector of components of $\lambda_{i}\>\hat{\underline{F}}_{h\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ (respectively $\lambda_{i}\>\underline{F}_{h\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$).

Besides, quantity $\underline{v}_{h}^{\ast}$ is discretized over element $E\in\mathcal{E}_{i}^{\mathcal{I}}$ which reads:

[TABLE]

where $[\lambda_{\mathclose{}|\mathopen{}{E}}]$ is the matrix of linear FE shape functions over $E$ and $\underline{X}_{h,E}^{(i)}$ is the vector of components of $\underline{v}_{h\mathclose{}|\mathopen{}{E}}^{\ast}$.

Local problems $\mathcal{P}_{i}^{\mathcal{I}}$ (36) (or (38)) then take the matrix form:

[TABLE]

where $\hat{\underline{f}}_{h}^{(i)}$ and $\underline{X}_{h}^{(i)}$ are the generalized vectors which contain the combination of vectors $\hat{\underline{f}}_{h,\Gamma}^{(i)}$ for every edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$ and vectors $\underline{X}_{h,E}^{(i)}$ for every element $E\in\mathcal{E}_{i}^{\mathcal{I}}$, respectively. $\mathbb{A}^{(i)}$ can be seen as the combination of mass matrices associated with every edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$.

Let us recall that conditions $\lambda_{i}\>\hat{\underline{F}}_{h\mathclose{}|\mathopen{}{\Gamma}}^{(i)}=\lambda_{i}\>\underline{F}_{h\mathclose{}|\mathopen{}{\Gamma}}^{(i)}$ along edges $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}\cap\partial_{2}\Omega$ must be verified for vertex $i\in\partial_{2}\Omega$. Those equalities read $\hat{\underline{f}}_{h,\Gamma}^{(i)}=\underline{f}_{h,\Gamma}^{(i)}$ over $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}\cap\partial_{2}\Omega$ and can be enforced using penalty terms in cost function (40). Thus, cost function (40) takes the matrix form:

[TABLE]

where the generalized vector $\underline{f}_{h}^{(i)}$ is the combination of known vectors $\underline{f}_{h,\Gamma}^{(i)}$ for every edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$ and $\mathbb{P}^{(i)}$ is a diagonal penalty matrix.

• •

In the case $p\geqslant 2$, known vector $\underline{f}_{h,\Gamma}^{(i)}$ over edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$ is calculated explicitly as it is the projection of $\lambda_{i}\>\Pi_{\dsigma_{h}}^{\Gamma}$ over the space of polynomial functions of degree $p$, where $\Pi_{\dsigma_{h}}^{\Gamma}$ depends on the FE stress field $\dsigma_{h}$ and the traction force density $\underline{F}_{d}$.

• •

In the case $p=1$, conditions (41) required for the calculation of quantity $\underline{f}_{h,\Gamma}^{(i)}$ over edge $\Gamma\in\mathcal{J}_{i}^{\mathcal{I}}$ mathematically read under the form:

[TABLE]

where $\dbeta_{\Gamma,j}^{(i)}$ represents the projection of $[\varphi_{\mathclose{}|\mathopen{}{\Gamma}}]$ over edge $\Gamma$ and over linear FE shape function $\lambda_{j}$. Grouping those conditions together, quantities $\underline{f}_{h,\Gamma}^{(i)}\quad\forall\>i\in\mathcal{I}_{\Gamma}^{\mathcal{J}}$ are recovered by solving a linear local system over edge $\Gamma\in\mathcal{J}$:

[TABLE]

where $\underline{f}_{h}^{\Gamma}$ denotes the generalized vector which is the combination of vectors $\underline{f}_{h,\Gamma}^{(i)}\quad\forall\>i\in\mathcal{I}_{\Gamma}^{\mathcal{J}}$ and $\mathbb{B}^{\Gamma}$ involves matrices $\dbeta_{\Gamma,j}^{(i)}\quad\forall\>j\in\mathcal{I}_{\Gamma}^{\mathcal{J}}\setminus\left\{i\right\},\ \forall\>i\in\mathcal{I}_{\Gamma}^{\mathcal{J}}$ and identity matrix $\mathbb{I}_{d}$.

Now, let us recall that kernel $\mathbf{Ker}(\mathbb{A}^{(i)})$ is the set of vectors $\underline{X}_{h}^{(i)}$ coming from continuous functions defined over patch $\Omega_{i}$ which vanish on $\mathcal{J}_{i}^{\mathcal{I}}\cap\partial\Omega$, i.e. $\underline{v}_{h}^{\ast}\in{\bar{\bm{\mathcal{U}}}_{h,0\mathclose{}|\mathopen{}{\Omega_{i}}}^{1}}$. In practice, this kernel is automatically reduced to $\underline{X}_{h}^{(i)}=\underline{0}$ by fixing some degrees of freedom over patch $\Omega_{i}$. This automatic procedure constitutes a great advantage of the method. Figure 4 displays free and fixed degrees of freedom for both an internal vertex $i\in\Omega\setminus\partial\Omega$ and a vertex $i\in\partial\Omega$.

Once the procedure to fix the kernel of matrix $\mathbb{A}^{(i)}$ has been performed, local problem (53) becomes:

[TABLE]

Finally, the problem to be solved over patch $\Omega_{i}$ consists of minimizing cost function (54) under constraint (57) coming from the global prolongation condition (31). Thus, introducing the Lagrangian:

[TABLE]

where $\underline{\tilde{X}}_{h}^{(i)}$ represents the vector of Lagrange multipliers associated with constraints (57), the system to be solved takes the matrix form:

[TABLE]

Finally, referring to section 10.0.1, one recovers tractions $\hat{\underline{F}}_{h}$ along edges $\Gamma\in\mathcal{J}$ from quantities $\hat{\underline{f}}_{h,\Gamma}^{(i)}$ for all vertices $i\in\mathcal{I}_{\Gamma}^{\mathcal{J}}$:

[TABLE]

As for the EET, an important and noteworthy point related to the procedure for constructing tractions $\hat{\underline{F}}_{h}=\displaystyle\sum_{i\in\mathcal{I}}\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}$ concerns its explicit and non-intrusive nature.

13 NUMERICAL RESULTS

All the two- and three-dimensional numerical experiments presented in this section are based on finite elements with Lagrange shape functions for both EET and EESPT methods, and hierarchical shape functions for the SPET. All the resolutions of local problems described in

14

8 at the element scale for both the EET and EESPT or at the patch scale for the SPET are performed by using a $p+k$ discretization with $k=3$ ($p$-refinement). The value of the penalty terms involved in the EESPT for enforcing the Neumann conditions on the edges which lie on $\partial_{2}\Omega$ is set to $10^{5}$. The selected examples are commonly used for industrial applications:

• •

the first case is a weight sensor whose sensor capacity ranges from $1$ kg to $2.5$ kg. This two-dimensional numerical test is a plane-stress linear elastic problem;

• •

the second case is an open hole specimen, which is represented by a three-dimensional plate containing a central hole. This structure is generally employed to study bolted joints in construction and machine design;

• •

the third case is a dome-shaped closure head, forming a part of a nuclear reactor vessel closure assembly;

• •

the fourth case is the hub of the main rotor of the NH90 helicopter developed by the Eurocopter company. It acts as a coupling sleeve between the helicopter frame and the rotor system.

For all the structures considered, the material is chosen to be isotropic, homogeneous, linear, and elastic with Young’s modulus $E=1$ and Poisson’s ratio $\nu=0.3$.

The three techniques used to calculate admissible stress fields are analyzed in terms of quality of the error estimators, computational cost and simplicity of practical implementation. The corresponding error estimates are denoted by $\theta=e_{\mathrm{cre}}(\hat{\underline{u}}_{h},\hat{\dsigma}_{h})$. The considered error estimates are also compared to the reference error, which usually replaces the unknown exact discretization error for practical purposes. Indeed, the exact error is estimated using a reference error lying in a discretized space much refined with respect to $\bm{\mathcal{U}}_{h}$; this reference error is thus computed from an “overkill” solution. Let us note that the same notation $\left\lVert\underline{e}_{h}\right\rVert_{u,\Omega}$ will be used to denote the energy norm of the reference error and that of the true or exact discretization error. The quality of the error estimators is measured with the usual effectivity index $\eta_{\bullet}$, which is the ratio between any error estimate $\theta_{\bullet}$ and the energy norm of the reference error $\left\lVert\underline{e}_{h}\right\rVert_{u,\Omega}$:

[TABLE]

the subscript $\bullet$ corresponding to the three estimators studied, i.e. “EET”, “SPET”  and “EESPT”. Thus, the accuracy of the error estimator is given by the value of $\eta$: $\eta>1$ indicates that $\theta$ is an overestimation of the reference error, whereas $\eta<1$ demonstrates that the former underestimates the latter.

14.1 Weight sensor

First, let us consider the structure of Figure 5, which contains two symmetric holes. The structure is clamped along the bottom-right side and subjected to a unit force density $\underline{f}=-\underline{y}$ along the top-left side. The other remaining boundaries are traction-free. The mesh used for the calculation consists of $11\,807$ linear triangular elements and $6\,320$ nodes (i.e. $12\,640$ degrees of freedom), see Figure 5. It has been heuristically adapted by refining the elements in the regions of the two holes where stresses are larger (indeed, only the high-stress zones are interesting for design purposes in mechanical design). The reference mesh used to compute the overall reference solution contains $3\,326\,963$ linear triangular elements and $1\,668\,711$ nodes (i.e. $3\,337\,422$ d.o.f.).

14.1.1 Comparison of the three error estimators

The local minimization step for both the EET and EESPT is performed using the cost function $J_{0}$ (19), which does not take the element size into account. A detailed analysis of the influence of the choice of the cost function is carried out in the next part. The FE stress field is shown in Figure 6 and the admissible stress fields obtained from the three techniques are shown in Figure 7. More precisely, Figure 6 and Figure 7 display the magnitudes \displaystyle\sqrt{\operatorname{Tr}\big{[}\dsigma_{h}\>\dsigma_{h}\big{]}} and \displaystyle\sqrt{\operatorname{Tr}\big{[}\hat{\dsigma}_{h}\>\hat{\dsigma}_{h}\big{]}} of the FE and admissible stress fields, respectively.

The direct calculation of $\left\lVert\underline{e}_{h}\right\rVert_{u,\Omega}$ is computationally unaffordable. As the structure is only subjected to kinematic conditions imposing zero displacement field on $\partial_{1}\Omega$, the exact value of the energy norm of the reference error has been calculated using:

[TABLE]

The required CPU time is about $2$ s, which is very convenient if one seeks to assess the energy norm of the reference error only. Indeed this computation is very low compared to that needed to compute the local contributions to the energy norm of the reference error whose corresponding computational cost reaches about $15$ hours, since the finite element and reference meshes are not being nested and the projection procedure is very time consuming. However, in order to have access to the spatial distribution of the element-by-element contributions to $\left\lVert\underline{e}_{h}\right\rVert_{u,\Omega}$, the calculation of $\left\lVert\underline{e}_{h}\right\rVert_{u,E}\quad\forall\>E\in\mathcal{E}$, with \left\lVert\bullet\right\rVert_{u,E}=(\int_{E}\operatorname{Tr}\big{[}\textbf{\em K}\>\depsil(\bullet)\>\depsil(\bullet)\big{]}\mathop{}\mathopen{}\mathrm{d}\Omega)^{1/2}, has been performed.

The performance of the error estimators is summarized in Table 1. First, one verifies that the global effectivity index is higher than $1$ for the three methods. Therefore the three error estimates $\theta_{\bullet}$ are upper bounds with respect to the energy norm of the reference error $\left\lVert\underline{e}_{h}\right\rVert_{u,\Omega}$. The behaviors of the error estimators are consistent. Moreover, the results reveal that the error estimate obtained with the SPET is more accurate than those given by the EET and EESPT: the global effectivity index $\eta$ is closer to $1$ for the SPET. The computational costs required by the EET and the EESPT are quite similar, whereas the CPU time needed for the SPET is four to five times higher.

Spatial distribution of local contributions to the energy norm of the reference error is shown in Figure 8, whereas that of local contributions to the error estimates for the three methods are displayed in Figure 9. Local contributions in each element of the FE mesh to the error estimates provide local indicators of the local energy norm of the reference error, that could be useful for a remeshing strategy in an adaptive procedure. From a qualitative viewpoint, the estimated error distribution is in good agreement with the reference error distribution. The main contributions of the error are concentrated around the top and bottom regions of the holes, i.e. located in the highly-loaded regions.

One defines a local effectivity index as the ratio between the element-by-element contribution to the error estimate and the one to the energy norm of the reference error. Spatial distribution of the local effectivity indices is shown in Figure 10. Let us recall that contributions to the effectivity index are computed with respect to the reference error. The local effectivity indices range between $0.675$ and $7.507$ for the EET; $0.586$ and $3.802$ for the SPET; $0.581$ and $7.375$ for the EESPT. One can observe that a quite large set of elements has a local effectivity index inferior to $1$, even in low-stress regions.

14.1.2 Influence of the choice of the cost function involved in the EET and the EESPT

Quantities ${\Delta\underline{q}^{(i)}}_{\mathclose{}|\mathopen{}{\Gamma}}$ involved in the cost function are $(\hat{\underline{b}}^{(i)}-\underline{b}^{(i)})_{\mathclose{}|\mathopen{}{\Gamma}}$ and $(\lambda_{i}\>\hat{\underline{F}}_{h}^{(i)}-\lambda_{i}\>\underline{F}_{h}^{(i)})_{\mathclose{}|\mathopen{}{\Gamma}}$ for the EET and the EESPT, respectively. Several cost functions have been considered for the least-squares minimization step over $\mathcal{J}_{i}^{\square}$, where $\square$ stands for $\mathcal{N}$ (resp. $\mathcal{I}$) in the case of the EET (resp. EESPT):

1. (i)

$J_{0}(\underline{q}^{(i)})=\frac{1}{2}\displaystyle\sum_{\Gamma\in\mathcal{J}_{i}^{\square}}(\Delta\underline{q}^{(i)}\cdot\Delta\underline{q}^{(i)})_{\mathclose{}|\mathopen{}{\Gamma}}$; 2. (ii)

$J_{1}(\underline{q}^{(i)})=\frac{1}{2}\displaystyle\sum_{\Gamma\in\mathcal{J}_{i}^{\square}}\frac{1}{l_{\Gamma}^{2}}\>(\Delta\underline{q}^{(i)}\cdot\Delta\underline{q}^{(i)})_{\mathclose{}|\mathopen{}{\Gamma}}$ which takes the size $l_{\Gamma}$ of each edge $\Gamma\in\mathcal{J}_{i}^{\square}$ into account; 3. (iii)

$J_{2}(\underline{q}^{(i)})=\frac{1}{2}\displaystyle\sum_{\Gamma\in\mathcal{J}_{i}^{\square}}\frac{1}{l_{\Gamma}^{2}}\>\frac{1+\nu}{E}\left(\frac{1-2\nu}{1-\nu}(\Delta\underline{q}^{(i)}\cdot\underline{n})_{\mathclose{}|\mathopen{}{\Gamma}}^{2}+2\>(\Pi\Delta\underline{q}^{(i)}\cdot\Pi\Delta\underline{q}^{(i)})_{\mathclose{}|\mathopen{}{\Gamma}}\right)$ which corresponds to the density of elastic energy stored in each edge $\Gamma\in\mathcal{J}_{i}^{\square}$. It has been explicitly obtained by assuming that the plane part of the strain field associated with the statically admissible stress field is continuous across the inter-element edges. The interested reader can refer to [45] for details.

Their influence on the quality of the two error estimators EET and EESPT has been studied; results are shown in Table 14.1.2. They reveal that the more physical information the cost function contains, the more accurate the yielded error estimator is. Indeed the global effectivity index experiences a $8\%$ decrease when the cost function $J_{2}$ is used instead of $J_{0}$. The bounds of the local effectivity indices calculated using each cost function are shown in Table 14.1.2 for both EET and EESPT methods.