# Méthode globale d'équilibrage élémentaire pour la reconstruction de champs statiquement admissibles dans le cadre des éléments finis

Valentine Rey, Pierre Gosselet, Christian Rey

► To cite this version:

Valentine Rey, Pierre Gosselet, Christian Rey. Méthode globale d'équilibrage élémentaire pour la reconstruction de champs statiquement admissibles dans le cadre des éléments finis. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722332

HAL Id: hal-01722332

https://hal.science/hal-01722332v1

Submitted on 3 Mar 2018

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

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

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Méthode globale d'équilibrage élémentaire pour la reconstruction de champs statiquement admissibles dans le cadre des éléments finis

Valentine REY<fnref n="1" />, Pierre GOSSELET<fnref n="1" />, Christian REY<fnref n="1" />

**Résumé** — Ce papier traite des techniques de reconstruction de champs statiquement admissibles dans le cadre de l'estimation d'erreur *a posteriori*, le but étant de fournir une borne supérieure stricte de la norme énergétique de l'erreur de discrétisation. L'étape de reconstruction est souvent complexe et décisive puisque l'efficacité de l'estimateur dépend fortement de la qualité du champ reconstruit. Nous proposons une nouvelle technique de reconstruction se basant sur la même hypothèse que la méthode *Element Equilibration Technique* (EET) mais l'exploitant de manière globale, afin de décrire tout l'espace d'optimisation. Nous détaillons cette nouvelle technique et la comparons sur un exemple 2D aux estimateurs EET et *Flux-free* classiques.

**Mots clés** — méthode des éléments finis, vérification, champs statiquement admissibles

## 1 Introduction

Cet article traite des méthodes de reconstruction de champs statiquement admissibles (SA) à partir d'une résolution éléments finis en déplacement. Cette étape est nécessaire à l'obtention d'une borne majorante de l'erreur de discrétisation (distance entre la solution éléments finis et la solution exacte inconnue). Les méthodes basées sur l'utilisation de champs SA, comme l'erreur en relation de comportement [6, 8] ou la méthode des résidus équilibrés [2, 1], fournissent des estimateurs efficaces, pour des quantités globales ou locales [9] même dans des cadres non-linéaires [12] ou pour les méthodes de décomposition de domaine [15].

La construction d'un champ SA est une étape cruciale (la qualité de l'estimateur dépend du champ reconstruit) et complexe. Le champ SA optimal est celui qui minimise l'énergie complémentaire et peut être obtenu par une approche duale [5]. Cette approche étant extremement coûteuse, on lui préfère des techniques de post-traitement de la solution éléments finis comme la méthode EET [8, 11] et ses variantes [7, 16] ou la technique *Flux-free* [14, 3]. Elles impliquent un grand nombre de résolutions sur des éléments ou sur des star-patchs d'éléments et sont souvent délicates à mettre en oeuvre.

Nous revisitons la méthode EET en proposant une écriture globale de sa principale hypothèse : la condition de prolongement forte. Ainsi, nous parvenons à décrire l'ensemble des champs SA satisfaisant cette condition. Malgré la grande taille du système global à résoudre, la mise en oeuvre est facilitée et de coût faible grâce à l'exploitation des propriétés géométriques et topologiques du maillage [19].

Dans la Section 2, nous introduisons le problème type et la notion d'erreur en relation de comportement, puis rappelons les bases de la EET dans la Section 3. La Section 4 est consacrée à la présentation de la nouvelle méthode et les résultats sont présentés dans la Section 5, où la méthode est comparée à la EET classique et à la méthode *Flux-free*.

## 2 Problème de référence et notion d'estimation d'erreur

Considérons l'équilibre statique d'une structure occupant l'ouvert $\Omega \subset \mathbb{R}^d$ soumise à une force volumique $f$ dans $\Omega$, à un effort de traction $g$ sur $\partial_g \Omega$ et à un déplacement imposé $u_d$ sur la partie complémentaire du bord $\partial_u \Omega \neq \emptyset$. Nous nous plaçons dans l'hypothèse des petites perturbations et nous considérons le matériau linéaire élastique, caractérisé par le tenseur d'élasticité de Hooke $H$. $u$ est le champ de déplacement, $\varepsilon(u)$ la partie symétrique du gradient de $u$, $\sigma$ le tenseur de contraintes de Cauchy.

<footnote n="1">¹ LMT-Cachan, ENS-Cachan/CNRS/UPMC/Pres UniverSud Paris, prenom.nom@lmt.ens-cachan.fr</footnote>
Nous introduisons :

- Le sous-espace affine des champs cinématiquement admissibles (CA)

$$ CA = \{ \underline{u} \in (\mathbb{H}^1(\Omega))^d, \underline{u} = \underline{u}_d \text{ sur } \partial_u \Omega \} \quad (1) $$

et notons $CA^0$ l'espace vectoriel associé.

- Le sous-espace affine des champs statiquement admissibles (SA)

$$ SA = \left\{ \underline{\tau} \in (L^2(\Omega))^d \right\}_{\text{sym}} ; \forall \underline{v} \in CA^0, \int_{\Omega} \underline{\tau} : \underline{\varepsilon}(\underline{v}) d\Omega = \int_{\Omega} \underline{f} \cdot \underline{v} d\Omega + \int_{\partial_g \Omega} \underline{g} \cdot \underline{v} dS \right\} \quad (2) $$

- L'erreur en relation de comportement

$$ e_{CR_{\Omega}}(\underline{u}, \underline{\sigma}) = \|\underline{\sigma} - \mathbb{H} : \underline{\varepsilon}(\underline{u})\|_{\mathbb{H}^{-1}, \Omega} \text{ avec } \|\underline{x}\|_{\mathbb{H}^{-1}, \Omega} = \sqrt{\int_{\Omega} (\underline{x} : \mathbb{H}^{-1} : \underline{x}) d\Omega} \quad (3) $$

Le problème mécanique posé sur $\Omega$ s'écrit :

$$ \text{Trouver } (\underline{u}_{ex}, \underline{\sigma}_{ex}) \in CA \times SA \text{ tel que } e_{CR_{\Omega}}(\underline{u}_{ex}, \underline{\sigma}_{ex}) = 0 $$

La solution de ce problème, dite solution exacte, existe et est unique.

## 2.1 Approximation éléments finis

Soit $\Omega_H$ un maillage par des triangles de $\bar{\Omega}$. On note $\mathcal{T}$ l'ensemble des éléments de $\Omega_H$, $\mathcal{E}$ l'ensemble des arêtes et $\mathcal{V}$ l'ensemble des sommets. Notons $|\mathcal{A}|$ le cardinal d'un ensemble $\mathcal{A}$. Introduisons l'opérateur de bord $\partial$ qui extrait les arêtes d'un élément et les sommets d'une arête.

Au domaine $\Omega_H$, nous associons le sous-espace de dimension finie $CA_H$ de $CA$. La méthode des éléments finis consiste à chercher $\underline{u}_H \in CA_H$ et $\underline{\sigma}_H = \mathbb{H} : \underline{\varepsilon}(\underline{u}_H)$ tels que :

$$ \forall \underline{v}_H \in CA_H^0 \int_{\Omega_H} \underline{\sigma}_H : \underline{\varepsilon}(\underline{v}_H) d\Omega = \int_{\Omega_H} \underline{f} \cdot \underline{v}_H d\Omega + \int_{\partial_g \Omega_H} \underline{g} \cdot \underline{v}_H dS \quad (4) $$

Après avoir introduit la matrice $\varphi_H$ des fonctions de forme et le vecteur des inconnues nodales $\mathbf{u}$ tel que $\underline{u}_H = \varphi_H \mathbf{u}$, nous obtenons le système linéaire :

$$ \begin{pmatrix} \mathbf{K}_{rr} & \mathbf{K}_{rd} \\ \mathbf{K}_{dr} & \mathbf{K}_{dd} \end{pmatrix} \begin{pmatrix} \mathbf{u}_r \\ \mathbf{u}_d \end{pmatrix} = \begin{pmatrix} \mathbf{f}_r \\ \mathbf{f}_d \end{pmatrix} + \begin{pmatrix} 0 \\ \lambda_d \end{pmatrix} \quad (5) $$

où $\mathbf{K}$ est la matrice de rigidité du domaine $\Omega_H$ et $\mathbf{f}$ est le vecteur des efforts généralisés ; l'indice $d$ est associé aux degrés de liberté de Dirichlet et l'indice $r$ représente les degrés de liberté tel que les inconnues soient $\mathbf{u}_r$ et $\lambda_d$ où le vecteur $\lambda_d$ représente les réactions nodales.

## 2.2 L'erreur en relation de comportement

La mesure de la qualité de la solution via l'erreur en relation de comportement consiste à construire des champs admissibles à partir des champs obtenus par la méthode des éléments finis et mesurer la non-vérification de la relation de comportement :

$$ \left\{ \begin{array}{l} \underline{u}_H \in CA_H \\ \underline{\sigma}_H = \mathbb{H} : \underline{\varepsilon}(\underline{u}_H) \\ e_{CR_{\Omega_H}}(\underline{u}_H, \underline{\sigma}_H) = 0 \end{array} \right. \rightarrow \left\{ \begin{array}{l} \hat{\underline{u}} \in CA \\ \hat{\underline{\sigma}} \in SA \\ e_{CR_{\Omega_H}}(\hat{\underline{u}}, \hat{\underline{\sigma}}) \ge 0 \end{array} \right. \quad (6) $$

$e_{CR_{\Omega_H}}(\hat{\underline{u}}, \hat{\underline{\sigma}})$ fournit une borne stricte de la norme énergétique de l'erreur en contrainte et en déplacement [6, 10].
Ici, nous considérons que la construction de $\hat{\underline{u}} \in CA_H$ est aisée puisque $CA_H \subset CA$ (nous choisissons donc $\hat{\underline{u}} = \underline{u}_H$). Construire le champ $\hat{\underline{\sigma}}$ est bien plus complexe. Plusieurs techniques existent. Citons par exemple la méthode d'équilibrage des éléments [8, 11] (EET) détaillée dans la Section 3, la méthode Flux-free [13] et la méthode d'équilibrage élémentaire par star-patchs (EESPT) [7, 17].

Une fois $\hat{\underline{\sigma}} \in SA$ construit, le calcul de l'erreur est effectué sur chaque élément $T$ de la manière suivante :

$$e_{\text{CR}_T}(\underline{u}_H, \hat{\underline{\sigma}}) = \frac{1}{2} \int_T \hat{\underline{\sigma}} : \mathbb{H}^{-1} : \hat{\underline{\sigma}} d\Omega + \frac{1}{2} \int_T \underline{\varepsilon}(\underline{u}_H) : \mathbb{H} : \underline{\varepsilon}(\underline{u}_H) d\Omega - \int_T \hat{\underline{\sigma}} : \underline{\varepsilon}(\underline{u}_H) \quad (7)$$

# 3 Méthode EET classique

La méthode EET se décompose en 2 étapes : premièrement, il s'agit de construire des efforts de tractions équilibrés $\hat{\underline{F}}^{\Gamma}$ sur chaque arête $(\Gamma)$ des éléments. Les champs SA sont ensuite calculés indépendamment sur chaque élément $(T)$ avec comme conditions de Neumann $\underline{\sigma}_{\underline{T}} \cdot \underline{n}^{\Gamma} = \delta_{\underline{T}}^{\Gamma} \hat{\underline{F}}^{\Gamma}$ (le rôle de $\delta_{\underline{T}}^{\Gamma} = \pm 1$ est explicité dans l'équation (10)).

Afin de simplifier l'estimation de l'erreur (7), l'hypothèse de condition de prolongement forte est formulée :

$$\forall T \in \mathcal{T}, \forall \underline{v} \in CA_H, \int_T (\underline{\hat{\sigma}} - \underline{\sigma}_H) : \underline{\varepsilon}(\underline{v}) d\Omega = 0 \quad (8)$$

qui traduit que le champ SA reconstruit développe le même travail que la solution éléments finis dans l'espace éléments finis.

Cette condition permet d'éliminer le terme couplé dans (7). Par conséquent, seule l'énergie complémentaire du champ SA est estimée et plus cette énergie est faible, meilleure est l'estimation.

## 3.1 Efforts de traction équilibrés

Sur chaque arête $\Gamma$, nous devons contruire des tractions $\hat{\underline{F}}^{\Gamma}$ satisfaisant la condition de prolongement forte.

$$\left\{
\begin{array}{l}
\underline{\hat{\sigma}}_{\underline{T}} \cdot \underline{n}^{\Gamma} = \delta_{\underline{T}}^{\Gamma} \hat{\underline{F}}^{\Gamma}, \quad \forall T \in \mathcal{T}, \forall \Gamma \in \partial T \\
\underline{\hat{F}}^{\Gamma} = \underline{g} \quad \text{sur } \partial_g \Omega_H \\
0 = \int_T (\underline{\hat{\sigma}} - \underline{\sigma}_H) : \underline{\varepsilon}(\underline{v}) d\Omega, \quad \forall \underline{v} \in CA_H,
\end{array}
\right.
\quad (9)$$

Le coefficient $\delta_{\underline{T}}^{\Gamma} = \pm 1$ est utilisé pour assurer le principe d'action-réaction de deux éléments voisins

$$\left\{
\begin{array}{l}
\forall \Gamma \in \partial \Omega_H, T = \partial^{-1}\Gamma, \delta_{\underline{T}}^{\Gamma} = 1 \\
\forall \Gamma \in \mathcal{E} \setminus \partial \Omega_H, \{T, T'\} = \partial^{-1}\Gamma, \delta_{\underline{T}}^{\Gamma} + \delta_{\underline{T}'}^{\Gamma} = 0
\end{array}
\right.
\quad (10)$$

Puisque les mouvements de corps rigides appartiennent à l'espace $CA_H$, la condition de prolongement forte implique l'équilibre des éléments vis-à-vis des corps rigides. L'alternative de Fredholm est vérifiée et $\underline{\hat{\sigma}}$ est bien défini.

En développant la condition de prolongement et en utilisant l'équilibre éléments finis, nous obtenons :

$$\forall T \in \mathcal{T}, \sum_{\Gamma \in \partial T} \delta_{\underline{T}}^{\Gamma} \int_{\Gamma} \underline{\hat{F}}^{\Gamma} \cdot \underline{v} d\Gamma = \int_{T} (\underline{\underline{\sigma}}_{\underline{H}} : \underline{\underline{\varepsilon}}(\underline{v}) - \underline{f} \cdot \underline{v}) d\Omega = \underline{\dot{R}}_{\underline{T}}(\underline{v}) \quad (11)$$

## 3.2 Construction des efforts de traction équilibrés

Dans la version orginale de la EET [8], chaque traction $\hat{\underline{F}}^{\Gamma}$ est supposée linéaire. L'équation (11) est alors testée sur chaque fonction de forme afin que l'évaluation soit limitée à leur support, nommé star-patch. On note $\omega_N$ le support de la fonction de forme $\phi_H^N$ associée au noeud $N \in \mathcal{V}$. Nous pouvons calculer le vecteur $\underline{\hat{W}}_{\underline{N}}^{\Gamma}$ :

$$\underline{\hat{W}}_{\underline{N}}^{\Gamma} = \int_{\Gamma} \left( \underline{\hat{F}}^{\Gamma} \cdot \underline{\varepsilon}_{\underline{x}} \phi_{\underline{H}}^{N} \right) d\Gamma \quad (12)$$
Nous réécrivons alors la condition de prolongement forte comme plusieurs petits problèmes découpés. La nature du système linéaire associé au star-patch dépend de la position du noeud. Pour les noeuds internes et ceux sur le contour de Dirichlet, le système est sous-déterminé si bien qu'il est possible d'optimiser les efforts de traction calculés.

Une fois ces systèmes résolus, il faut calculer les tractions à partir des travaux $\hat{W}_N^\Gamma$. En faisant l'hy-
pothèse que les tractions se décomposent sur l'espace éléments finis, il ne reste alors qu'à inverser des
matrices de masse (2×2).

3.3 Calcul du champ SA sur chaque élément

Une fois les tractions ($\hat{F}^{\Gamma}$) déterminées, celles-ci sont appliquées comme conditions de Neumann pour chaque élément. Pour des cas simples de chargement volumique, le champ de contrainte élémentaire exact peut être obtenu en effectuant une résolution polynômiale. Plus généralement, une solution numérique satisfaisante est obtenue avec une approximation éléments finis plus précise en se plaçant dans un espace $CA_p^T$ plus riche (si cet espace n'est pas assez riche, l'énergie est sous-estimée et la borne obtenue n'est plus stricte).

$$
\left\{
\begin{array}{l}
\displaystyle \int_T \underline{\underline{\sigma}}_p : \underline{\underline{\epsilon}}(\underline{\Phi}) d\Omega - \left( \int_T \underline{\underline{f}} \cdot \underline{\underline{\Phi}} d\Omega + \sum_{\Gamma \in \partial T} \int_{\Gamma} \delta_T^\Gamma \hat{F}^\Gamma \cdot \underline{\underline{\Phi}} dS \right) = 0, \forall \underline{\underline{\Phi}} \in CA_p^T \\[1em]
\underline{\underline{\sigma}}_p = \mathbb{H} : \underline{\underline{\epsilon}}(\underline{\mu}_p), \quad \underline{\underline{\mu}}_p \in CA_p^T
\end{array}
\right.
\qquad (13)
$$

4 Nouvelle technique de reconstruction de champs SA

La nouvelle technique que nous proposons repose sur l'interprétation de la condition de prolonge-
ment forte à l'échelle globale donnant la possibilité d'une optimisation globale. L'acronyme de cette
méthode est STARFLEET pour star-free and lazy element equilibration technique où nous insistons sur
la facilité d'implémentation et sur le rôle des star-patchs qui aparaissent naturellement dans le dévelop-
pement de la méthode.

Soit $\mathcal{E}$ l'ensemble des arêtes internes (celles séparant deux éléments), $\mathcal{E}_{\partial}$ les arêtes externes, $\mathcal{V}$ l'ensemble des sommets internes (qui ne sont pas sur la peau du domaine $\Omega_H$). Rappelons que, afin d'obtenir un maillage valide pour l'approximation éléments finis, un élément ne peut avoir 2 arêtes sur la peau du domaine [4], donc $|\mathcal{E}_{\partial}| = |\mathcal{V}_{\partial}|$.

4.1 Elimination des arêtes externes

La première étape de cette méthode est d'incorporer à l'équation des résidus (11) l'information
connue sur le bord du domaine. Pour un champ test $v_H \in CA_H :

- Si l'élément T est interne alors $R_T(v_H) = \mathring{R}_T(v_H)$

- Si l'élément T possède une arête de Neumann $\Gamma_g : R_T(v_H) = \mathring{R}_T(v_H) - \int_{\Gamma_g} g \cdot v_H d\Gamma$

- si l'élément T possède une arête de Dirichlet $\Gamma_u : R_T(v_H) = \mathring{R}_T(v_H) - \sum_{V \in \partial \Gamma_u} \alpha_V^{\Gamma_d} \lambda_d(V) \cdot v_H(V)$, où $V$ représente les sommets de $\Gamma_D$. $\lambda_d(V)$ est associé à la valeur de réaction nodale, connue de la résolution éléments finis, et $v_H(V)$ est la valeur du champ test au sommet $(\alpha_V^\Gamma)$ étant une partition de l'effort nodal de réaction entre les éléments qu'il impacte).

4.2 Equation de prolongement forte à l'échelle globale

Il est nécessaire que les coefficients $(\delta_T^\Gamma)$ soient cohérents avec l'orientation des éléments et arêtes.
En 2D, cela revient à choisir, pour chaque arête, un sommet initial et un sommet final et pour chaque
élément, un sens de rotation. Soit $\Delta$ la matrice rectangulaire de $|\mathcal{T}|$ lignes et $|\mathcal{E}|$ colonnes de coefficients
$(\delta_T^\Gamma)$, appelée matrice d'incidence. Il s'agit d'une matrice très creuse et booléenne signée. Chaque ligne
contient au plus 3 coefficients non nuls (2 pour une ligne associée à un élément du bord), et chaque
colonne contient 2 coefficients non nuls.
Pour un champ $\underline{v} \in CA_H$ donné, notons $\hat{\mathbf{W}}(\underline{v})$ le vecteur de taille $|\mathcal{E}|$ des travaux des tractions par le champ $\underline{v}$ et notons $\mathbf{R}(\underline{v})$ le vecteur de taille $|\mathcal{T}|$ :

La condition de prolongement forte écrite pour ces champs tests est :

$$ \forall \underline{v}_H \in CA_H, \Delta \hat{\mathbf{W}}(\underline{v}_H) = \mathbf{R}(\underline{v}_H) \quad (14) $$

Dans la méthode EET classique, les fonctions test sont les fonctions de forme. Dans notre approche, la base test est la base canonique $\mathcal{B} = (1\underline{\epsilon}_x, x\underline{\epsilon}_x, y\underline{\epsilon}_x, 1\underline{\epsilon}_y, x\underline{\epsilon}_y, y\underline{\epsilon}_y)$. Puisque la condition de prolongement est écrite indépendamment sur chaque élément, cette écriture de la condition est équivalente à celle de la EET.

Satisfaire l'équation (14) n'est pas suffisant. En effet, les travaux doivent être linéaires vis-à-vis des fonctions test :

$$ \forall (\underline{v}_1, \underline{v}_2) \in CA_H \times CA_H, \forall \alpha \in \mathbb{R}, \hat{\mathbf{W}}(\underline{v}_1 + \alpha \underline{v}_2) = \hat{\mathbf{W}}(\underline{v}_1) + \alpha \hat{\mathbf{W}}(\underline{v}_2) $$

Cette propriété implique que les travaux sur les arêtes ne sont pas indépendants. Les coordonnées spatiales sont liées : pour chaque arête $\Gamma$, il existe des coefficients $(a_\Gamma, b_\Gamma, c_\Gamma)$ tels que,

$$ (x, y) \in \Gamma \iff a_{\Gamma}x + b_{\Gamma}y + c_{\Gamma} = 0 \quad (15) $$

qui peuvent être choisis tels que $a_\Gamma^2 + b_\Gamma^2 = |\Gamma|^2$. Ainsi, dans chaque direction $\underline{\epsilon}$,

$$ a_\Gamma W^\Gamma(x\underline{\epsilon}) + b_\Gamma W^\Gamma(y\underline{\epsilon}) + c_\Gamma W^\Gamma(\underline{\epsilon}) = W^\Gamma((a_\Gamma x + b_\Gamma y + c_\Gamma)\underline{\epsilon}) = 0 $$

Soit $\mathbf{a}, \mathbf{b}, \mathbf{c}$ les matrices diagonales respectives des coefficients $a_\Gamma, b_\Gamma$ and $c_\Gamma$.

Puisque les composantes $(\underline{\epsilon}_x, \underline{\epsilon}_y)$ sont indépendantes, nous restreignons notre analyse au problème scolaire (thermique) et la condition de prolongement s'écrit :

$$ \mathbf{G} \begin{pmatrix} \hat{\mathbf{W}}(1) \\ \hat{\mathbf{W}}(x) \\ \hat{\mathbf{W}}(y) \end{pmatrix} = \begin{pmatrix} \mathbf{R}(1) \\ \mathbf{R}(x) \\ \mathbf{R}(y) \\ 0 \end{pmatrix} \quad \text{avec } \mathbf{G} = \begin{pmatrix} \Delta & & \\ & \Delta & \\ & & \Delta \\ \mathbf{c} & \mathbf{a} & \mathbf{b} \end{pmatrix} \quad (16) $$

## 4.3 Etude du système

Nous détaillons à présent les noyaux des différentes matrices introduites afin de s'assurer que le problème est bien posé et caractériser l'espace d'optimisation.

### 4.3.1 Noyaux de $\Delta$

$\Delta$ possède un unique vecteur dans son noyau à gauche qui est $\mathbf{i} = (1 \dots 1)^T$ si tous les éléments ont la même orientation. En effet, une arête interne est vue par deux éléments : une fois positivement, une fois négativement. Le rang de $\Delta$ étant $|\mathcal{T}| - 1$, la dimension du noyau à droite est $|\mathcal{E}| - |\mathcal{T}| + 1$.

Rappelons l'identité d'Euler :

$$ |\mathcal{T}| - |\mathcal{E}| + |\mathcal{V}| = c - h \quad (17) $$

où $c$ est le nombre de composantes connexes et $h$ le nombre de trous. Par la suite, nous aurons $c = 1$.

Rappelons également que $|\mathcal{V}_\partial| = |\mathcal{E}_\partial|$ et que $|\mathcal{V}_\partial| + |\hat{\mathcal{V}}| = |\mathcal{V}|$. La dimension du noyau à droite de $\Delta$ est $|\hat{\mathcal{V}}| + h$

- Chaque noeud interne définit un star-patch. En combinant les arêtes irradiant le noeud, les éléments sont comptés deux fois avec des signes opposés.

- Les arêtes internes irradiant les noeuds du contour du trou forment un chemin fermé entourant le trou.

Cette base peut être facilement obtenue via une analyse de graphe (ou en utilisant la topologie algébrique [18]). Par la suite, $N_V$ est le vecteur du noyau associé au star-patch du noeud V et $N_h$ celui associé au trou h.
### 4.3.2 Noyau à droite de G

Considérons un star-patch de noeud central V de coordonnées $(x_V, y_V)$, et $\Gamma$ une branche d'extremité $A^\Gamma = (x_A, y_A)$. Les coefficients de la branche sont proportionnels à :

$$a_{\Gamma} = (y_{A} - y_{V}), b_{\Gamma} = -(x_{A} - x_{V}), c_{\Gamma} = y_{V}(x_{A} - x_{V}) - x_{V}(y_{A} - y_{V}) \quad (18)$$

Le centre étant un noeud sommet de l'arête, nous pouvons écrire $a_\Gamma x_V + b_\Gamma y_V + c_\Gamma = 0$. Alors

$$\left\{ \begin{array}{l} W_{I}^{\Gamma} = N_{V} \\ W_{x}^{\Gamma} = x_{V}N_{V} \\ W_{y}^{\Gamma} = y_{V}N_{V} \end{array} \right. \qquad (19)$$

appartient au noyau de G. Chaque noeud interne est associé à un vecteur du noyau et aucun autre vecteur ne peut appartenir à ce noyau. Notons ce noyau $\mathcal{N}_G = (\mathbf{N}_V^T, \mathbf{x}_V \mathbf{N}_V^T, \mathbf{y}_V \mathbf{N}_V^T)^T$ où $\mathbf{x}_V$ et $\mathbf{y}_V$ sont les matrices diagonales des coordonnées des noeuds internes.

### 4.3.3 Noyau à gauche de G

Décrivons le noyau à gauche de G afin de s'assurer que le problème est bien posé. En effet, l'alternative de Fredholm doit être vérifiée. Ainsi, nous devons garantir que le second membre du système est orthogonal au noyau à gauche de l'opérateur. Le théorème du rang donne :

$$\dim(\ker(\mathbf{G}^T)) = 3|\mathcal{T}| + |\mathcal{E}| - (3|\mathcal{E}| - |\mathcal{V}|) = |\mathcal{V}| + |\mathcal{E}_\partial| + 3|\mathcal{T}| - (2|\mathcal{E}| - |\mathcal{E}_\partial|) = |\mathcal{V}|$$

L'existence du noyau à gauche est directement liée à la possibilité de construire une fonction chapeau (fonction de forme) pour chaque star patch (fonction continue, linéaire par élément, nulle sur le bord).

Multiplier à gauche le second membre de (16) par un vecteur du noyau à gauche correspond au calcul du résidu pour une fonction de forme sur le star-patch, ce qui donne 0 par l'orthogonalité de Galerkin. Le problème est donc bien posé.

## 4.4 Résolution efficace

La connaissance des propriétés de ce système très creux nous permet d’envisager une résolution efficace dont la complexité est faible. Nous avons exhibé une base $\mathcal{N}_G$ du noyau de G si bien que toutes les solutions sont de la forme $\hat{\mathbf{W}} = \hat{\mathbf{W}}_0 + \mathcal{N}_G\gamma$ où $\gamma$ est un vecteur de taille $|\mathcal{V}|$ utilisé pour optimiser $\hat{\mathbf{W}}$ vis-à-vis d’un critère bien choisi. La matrice $\Delta$ étant la description algébrique du maillage, il existe une forme normale de Smith :

$\exists(\mathbf{U}, \mathbf{V})$ matrices booléennes signées inversibles d'inverses booléennes signées telles que

$$\mathbf{U}\mathbf{\Delta}\mathbf{V} = \begin{pmatrix} 1 & & & 0 & \dots & 0 \\ & \ddots & & \vdots & & \vdots \\ & & 1 & 0 & & 0 \\ & \dots & & 0 & \dots & 0 \end{pmatrix} \quad \text{avec} \quad \begin{array}{l} \mathbf{U} = [\tilde{\mathbf{U}}^T, \mathbf{i}]^T \\ \mathbf{V} = [\tilde{\mathbf{V}}, \mathbf{N}] \end{array} \qquad (20)$$

La factorisation est alors facilement réalisable par un pivot de Gauss total. La matrice $\mathbf{\Delta}$ étant très creuse, la complexité de cette opération est linéaire. Une fois la factorisation de Smith effectuée, il est possible d'utiliser une factorisation LU avec pivot afin de préserver la structure très creuse de la matrice.

## 4.5 Critères de choix de l'élément du noyau

Nous avons déterminé une solution $\hat{\mathbf{W}}_0$ et la base $\mathcal{N}_G$ du noyau de G. Quel que soit $\gamma$, $\hat{\mathbf{W}} = \hat{\mathbf{W}}_0 + \mathcal{N}_G\gamma$ est le travail d'un champ SA satisfaisant la condition de prolongement. Cela laisse une opportunité d'optimisation dans un espace de taille $|\mathcal{V}|$ pour un problème scalaire.
Un premier choix peut être de se rapprocher des travaux développés par la solution éléments finis, au sens d'une norme M. Les résultats seront donnés pour $M = I$ et l'estimateur basé sur ce choix est noté $\text{STARFLEET}_{\|\cdot\|_2}$. Nous avons montré dans [19] qu'une norme bien choisie conduit à la minimisation faite classiquement dans la EET. Il est aussi possible de rendre optimale l'estimation de l'erreur en minimisant l'énergie complémentaire. En contrepartie, le coût de calcul associé à ce choix (indiqué par la suite $\text{STARFLEET}_{\|\cdot\|_{\text{erdc}}}$) est plus élevé que pour les précédents.

## 4.6 Résolution précise des problèmes élémentaires

Pour résoudre chaque problème élémentaire, il est nécessaire de construire une base plus riche que la base éléments finis. Nous choisissons d'augmenter le degré de la base canonique. Ainsi, $u_p$ et $v_p$ sont développés sur la base $x, y, xy, x^2, y^2, \dots$. Les travaux correspondant au premier ordre ont été déterminés précédemment. Pour déterminer les travaux d'ordre supérieurs, nous faisons l'hypothèse que ces travaux sont issus d'une distribution d'effort linéaire pour chaque arête.

# 5 Résultats

Considérons la structure suivante (Figure 1) dans le cadre statique, isotrope ($E = 1$ et $\nu = 0.3$) et linéaire, encastrée à sa base et soumise à un effort de traction sur le bord supérieur droit. La Figure 2 donne l'évolution des estimateurs en fonction du pas du maillage $h$. La résolution élémentaire a été effectuée en utilisant des fonctions de forme de degré $p+2$. On donne aussi la pente correspondant à $h^{2/3}$ (ce qui correspond à la singularité du problème).

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Conditions aux limites et chargement</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Résolution (p+2)</figcaption></figure>

Nous observons la nouvelle méthode de reconstruction mène à un estimateur aux performances comparables à la méthode EET. L'estimateur $\text{STARFLEET}_{\|\cdot\|_2}$ se comporte moins bien quand le maillage est fin. Cela est sans doute causé par la singularité du problème. Le critère de minimisation de l'énergie complémentaire, plus coûteux, donne de meilleurs résultats que les deux autres approches. Notons que l'estimation fournie reste proche de celle provenant de la EET.

# 6 Conclusion

Nous avons présenté une nouvelle méthode de reconstruction de champs statiquement admissibles pour l'estimation d'erreur *a posteriori*. Cette méthode reprend l'ingrédient principal de la EET qui est la condition de prolongement forte et la réécrit de manière globale, ce qui offre plusieurs possibilités d'optimisation. En exploitant les propriétés géométriques et topologiques du maillage, cette technique apporte une compréhension géométrique de la nature des systèmes issus de la condition de prolongement et permet une mise en oeuvre peu coûteuse profitant du caractère très creux de la matrice d'incidence.
Les résultats montrent que les bornes d’erreur obtenues et les temps de calcul sont très comparables à ceux de la technique EET. Le critère se basant sur la minimisation de l’estimateur permet d’étudier la pertinence de la condition de prolongement et montre que la EET en fait une exploitation tout à fait correcte.

La méthode s'étend aux éléments quadrangulaires, aux éléments d'ordre 2 et aux problèmes 3D pour lesquelles la mise en oeuvre semble très simplifiée par rapport aux techniques existantes. De plus, la méthode STARFLEET avec minimisation de l'énergie complémentaire est un estimateur efficace pour certains cas pathologiques (maillages avec triangles très applatis).

Références

[1] I. Babuška and A. Miller. *A feedback finite element method with a posteriori error estimation. part 1 : The finite element method and some basic properties of the a posteriori error estimator.* Computer Methods in Applied Mechanics and Engineering, 61 :1–40, 1987.

[2] I. Babuška and W. C. Rheinboldt. *A-posteriori error estimates for the finite element method*. International Journal for Numerical Methods in Engineering, 12(10) :1597–1615, 1978.

[3] R. Cottereau, P. Díez, and A. Huerta. *Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method*. Computational Mechanics, 44(4) :533–547, 2009.

[4] P.L. George. *Automatic mesh generation. Applications to finite element methods*. Wiley, 1991.

[5] Martin Kempeneers, Jean-François Debongnie, and Pierre Beckers. *Pure equilibrium tetrahedral finite elements for global error estimation by dual analysis*. International Journal for Numerical Methods in Engineering, 81(4) :513–536, 2009.

[6] P. Ladevèze. *Comparaison de modèles de milieux continus*. PhD thesis, Université P. et M. Curie, 1975.

[7] P. Ladevèze, L. Chamoin, and E. Florentin. *A new non-intrusive technique for the construction of admissible stress fields in model verification*. Computer Methods in Applied Mechanics and Engineering, 199 :766–777, 2010.

[8] P. Ladevèze and D. Leguillon. *Error estimate procedure in the finite element method and application*. SIAM Journal of Numerical Analysis, 20(3) :485–509, 1983.

[9] P. Ladevèze and N. Moës. *Adaptive control for finite element analysis in plasticity*. Computers and Structures, 73(1-5) :45–60, 1999.

[10] P. Ladevèze and J. P. Pelle. *Mastering Calculations in Linear and Nonlinear Mechanics*. Springer, New York, 2004.

[11] P. Ladevèze and P. Rougeot. *New advances on a posteriori error on constitutive relation in finite element analysis*. Computer Methods in Applied Mechanics and Engineering, 150(1-4) :239–249, 1997.

[12] F. Louf, J. P. Combe, and J. P. Pelle. *Constitutive error estimator for the control of contact problems involving friction*. Computers and Structures, 81(18-19) :1759–1772, 2003.

[13] N. Parés, P. Díez, and A. Huerta. *Subdomain-based flux-free a posteriori error estimators*. Computer Methods in Applied Mechanics and Engineering, 195(4-6) :297–323, 2006.

[14] Núria Parés, Pedro Díez, and Antonio Huerta. *Exact Bounds for Linear Outputs of the Advection-Diffusion-Reaction Equation Using Flux-Free Error Estimates*. SIAM J. Sci. Comput., 31(4) :3064–3089, 2009.

[15] A. Parret-Fréaud, C. Rey, P. Gosselet, and F. Feyel. *Fast estimation of discretization error for FE problems solved by domain decomposition*. Computer Methods in Applied Mechanics and Engineering, 199(49-52) :3315–3323, 2010.

[16] F. Pled, L. Chamoin, and P. Ladevèze. *On the techniques for constructing admissible stress fields in model verification : Performances on engineering examples*. International Journal for Numerical Methods in Engineering, 88(5) :409–441, 2011.

[17] F. Pled, L. Chamoin, and P. Ladevèze. *An enhanced method with local energy minimization for the robust a posteriori construction of equilibrated stress fields in finite element analyses*. Computational Mechanics, 49(3) :357–378, 2012.

[18] F. Rapetti, F. Dubois, and A. Bossavit. *Discrete vector potentials for nonsimply connected three-dimensional domains*. SIAM J. NUMER. ANAL., 41(4) :1505–1527, 2003.

[19] V. Rey, P. Gosselet, and C. Rey. *Study of the strong prolongation equation for the construction of statically admissible stress fields : implementation and optimization*. Soumis, 2013.