Thèse Pour

Obtenir Le Grade De

Docteur De L'université De Montpellier

Daniele Di Pietro

Alexandre Ern

Roland Becker

Françoise Krasucki

Vielen Dank Maria

Markus Michi

Julia Und

Daniele Di

Sylvie Granet

Sous la direction de

Keywords: Algorithmes adaptatifs, estimation d'erreur a posteriori, reconstruction équilibrée de tenseur de contraintes, espace d Adaptive algorithms, a posteriori error estimation, equilibrated stress tensor reconstruction, Arnold-Winther nite element space, Arnold-Falk-Winther nite element space, poroplasticity

A lg o rith mes ad ap tatifs p o u r la p o ro méc an iq u e et la p o ro -p lastic ité

Résumé Dans cette thèse nous développons des estimations d'erreur a posteriori par équilibrage de ux pour la poro-mécanique et la poro-plasticité. En se basant sur ces estimations, nous proposons des algorithmes adaptatifs pour la résolution numérique de problèmes en mécanique des sols. Le premier chapitre traite des problèmes en poro-élasticité linéaire. Nous obtenons une borne garantie sur l'erreur en utilisant des reconstructions équilibrées et Hpdivq-conformes de la vitesse de Darcy et du tenseur de contraintes mécaniques. Nous appliquons cette estimation dans un algorithme adaptif pour équilibrer les composantes de l'erreur provenant de la discrétisation en espace et en temps pour des simulations en deux dimensions. La contribution principale du chapitre porte sur la reconstruction symétrique du tenseur de contraintes. Dans le deuxième chapitre nous proposons une deuxième technique de reconstruction du tenseur de contraintes dans le cadre de l'élasticité nonlinéaire. En imposant la symétrie faiblement, cette technique améliore les temps de calcul et facilite l'implémentation. Nous démontrons l'éfcacité locale et globale des estimateurs obtenus avec cette reconstruction pour une grande classe de lois en hyperélasticité. En ajoutant un estimateur de l'erreur de linéarisation, nous introduisons des critères d'arrêt adaptatifs pour le solveur de linéarisation.

Le troisième chapitre est consacré à l'application industrielle des résultats obtenus. Nous appliquons un algorithme adaptatif à des problèmes poro-mécaniques en trois dimensions avec des lois de comportement mécanique élasto-plastiques.

Contexte industriel

En France, l'Agence Nationale pour la Gestion des Déchets Radioactifs (ANDRA 1 ) est chargée de concevoir un centre de stockage souterrain pour des déchets dont le niveau de radioactivité et la durée de vie sont élevés. Ces déchets proviennent principalement de l'industrie éléctronucléaire. Si sa création est autorisée, ce centre industriel de stockage géologique (Cigéo 2 , cf.

Figure 1.1) sera construit à Bure, dans l'est de la France, où se situe depuis 2007 un laboratoire de recherche à cet eet. Il sera nancé par les producteurs des déchets concernés, donc principalement par EDF, qui entretient le parc nucléaire français composé en 2017 de 58 réacteurs.

EDF R&D participe donc à la recherche pour garantir la sûreté, la faisabilité et l'optimisation des installations de Cigéo.

Le contexte industriel de cette thèse est la simulation du creusement de tunnels en 3D. Il s'agit d'un problème dont la complexité vient à la fois du couplage multi-physique et des géométries complexes (voir la Figure 1.2 pour un exemple). Ces simulations sont importantes pour estimer l'endommagement du sous-sol, qui pourrait favoriser le passage de radionucléides. Le modèle poro-mécanique décrivant le comportement de la roche destinée à recevoir les déchets radioactifs prend en compte deux phénomènes physiques résultants de la structure du sol : d'une part, le comportement mécanique du squelette solide déformable ; d'autre part, le comportement hydraulique de l'eau occupant son espace poreux interstitiel. Ces deux phénomènes doivent être considérés de façon couplée, car la déformation du squelette peut inuencer la pression de l'eau et, réciproquement, l'écoulement souterrain induit des déformations du squelette. On se place donc ici dans le cadre de la poro-mécanique [37].

Il est à ce jour impossible de réaliser des études poro-mécaniques avec une loi plastique sur une structure 3D complexe dans des temps de calcul raisonnables. Les deux problèmes principaux qu'EDF R&D a rencontrés pour ces simulations, et qui ont motivé cette thèse, sont les suivants :

• Le couplage des deux phénomènes entraine un problème à deux variables primales le déplacement du squelette mécanique et la pression de l'eau qui ont des unités et des ordres de grandeur diérents. Même si dans le modèle considéré la non linéarité du système à résoudre ne provient que de la partie mécanique, la présence de la pression rend problématique la dénition d'un critère de convergence purement algébrique pour le solveur de Newton.

1 http://www.andra.fr/ 2 http://www.cigeo.com/ 3 http://www.code-aster.org 4 http://www.code-aster.org/outils/homard/index.fr.html Algorithmes adaptatifs et estimation d'erreur a posteriori L'objectif de cette thèse est de proposer une solution pour les deux problèmes ci-dessus en développant des algorithmes de résolution adaptative basés sur des estimateurs d'erreur a posteriori obtenus par équilibrage de ux. Le fonctionnement d'un tel algorithme est illustré dans la Figure 1.3 : l'idée derrière le terme adaptatif est de concentrer l'eort de calcul là où l'erreur est grande, par exemple : en utilisant des techniques de remaillage pour obtenir une distribution homogène de l'erreur de discrétisation spatiale ; en adaptant le pas de temps pour équilibrer les contributions spatiales et temporelles de l'erreur de discrétisation ; ou bien en arrêtant les solveurs itératifs linéaires et non linéaires dès que leurs contributions à l'erreur globale deviennent négligeables. Dans le contexte de la modélisation de reservoirs de pétrole, le développement et l'implémentation d'un tel algorithme a permis de reduire considérablement les temps de calcul [114].

Il est clair qu'une estimation able des composantes de l'erreur citées plus haut (erreur en espace, en temps, de linéarisation, etc.) est indispensable pour pouvoir les équilibrer. Contrairement aux estimations d'erreur a priori qui sont utilisées pour contrôler la convergence d'une méthode de discrétisation (en général en fonction de la régularité de la solution analytique), les estimations a posteriori donnent une borne calculable (donc notamment indépendante de la solution analytique) de l'erreur. Cette borne dépend de la solution discrète et se calcule en conséquence après (a posteriori ) la résolution du problème discret. Il existe de nombreuses techniques d'estimation d'erreur a posteriori, mais, à notre connaissance, aucune n'a été appliquée à des problèmes poro-mécaniques non linéaires. Dans cette thèse nous avons choisi de nous baser sur la technique introduite dans [28,38] pour le problème de Poisson et de l'étendre à des problèmes poro-mécaniques. Cette méthode utilise des reconstructions équili-brées et Hpdivq-conformes de ux, où -en fonction du problème physique considéré -le ux satisfait une loi de type Fick et est en équilibre avec les forces extérieures. Cette technique ore plusieurs avantages qui répondent aux problèmes ci-dessus :

• L'expression de l'erreur en termes de ux permet de distinguer et, plus particulièrement, de comparer les diérentes composantes de l'erreur dans l'esprit de [55]. Elle permet notamment de comparer les contributions provenant des solveurs itératifs à l'erreur de discrétisation, ce qui est indispensable pour la dénition de critères d'arrêt adaptatifs pour ces solveurs. L'expression en terme de ux permet également de distinguer les deux phénomènes physiques du problème, et donc notamment de ne considérer que la partie mécanique qui induit la non linéarité quand on dénit ces critères d'arrêt.

• • Les reconstructions des ux ne dépendent pas de la relation entre le ux et la variable primale. L'implémentation de la méthode est donc indépendante de la loi de comportement, et se prête en conséquence à l'application dans des codes orant un grand choix de telles lois, comme c'est le cas dans Code_Aster.

Développements techniques eectués dans cette thèse

Dans les problèmes poro-mécaniques considérés dans cette thèse apparaissent deux ux : la vitesse de l'eau pour la partie hydraulique et le tenseur de contraintes pour la partie mécanique.

Pour obtenir une reconstruction équilibrée du ux hydraulique nous nous basons sur [28,38].

La contribution principale de cette thèse est de développer des méthodes de reconstruction du tenseur de contraintes en prenant en compte la symétrie de ce tenseur. Pour des calculs en 2D nous proposons deux techniques diérentes, en imposant la symétrie fortement ou faiblement. Une comparaison sur un problème en élasticité linéaire montre que les deux méthodes produisent des résultats similaires. Pour les simulations en 3D, nous nous focalisons sur la deuxième variante, qui est plus facile à implémenter et moins coûteuse en termes de temps de calcul.

Nous avons progressivement intégré le calcul de ces reconstructions de ux discrets et des estimateurs d'erreur dans le logiciel industriel Code_Aster, et avons obtenu les résultats suivants en réponse aux problèmes qui ont motivé cette thèse :

• Critère de convergence : Comme mentionné ci-dessus, l'estimation de l'erreur en termes de ux permet de séparer la partie hydraulique et la partie mécanique du problème. Nous pouvons ainsi dénir des estimateurs de l'erreur de discrétisation pour les deux parties, et de l'erreur de linéarisation pour la partie mécanique (puisque nous ne considérons que des lois hydrauliques de Darcy linéaires). En arrêtant le solveur de Newton dès que l'erreur provenant de la linéarisation est négligeable par rapport à celui de la discrétisation de la partie mécanique, nous obtenons un critère de convergence able qui prend en compte la physique du problème et qui, de plus, permet d'éviter des itérations inutiles.

• Ecacité et remaillage : Pour des problèmes en hyperélasticité, nous démontrons que l'estimation d'erreur a posteriori développée pendant cette thèse est ecace. De plus, nous avons vérié sur des tests numériques avec des solutions analytiques que les estimateurs réètent bien la distribution de l'erreur de discrétisation sur le maillage.

Ces propriétés sont importantes pour fournir à HOMARD des estimateurs adaptés au problème de poro-mécanique non linéaire, et pour garantir que les maillages soient bien adaptés au problème à chaque pas de calcul. En ranant les maillages en fonction de la distribution des estimateurs, nous avons obtenu des ordres de convergence (en termes d'erreur en fonction du nombre de degrés de liberté) plus élevés qu'avec un ranement uniforme.

• Application directe à diérentes lois de comportement mécanique : Dans Code_Aster, nous avons implémenté l'outil de reconstruction locale de façon à ce qu'il puisse être appliqué à toutes les lois de comportement utilisées au sein du formalisme poro-mécanique.

Structure du manuscrit

La structure du manuscrit repose sur l'ordre chronologique des étapes de la thèse. Ω un domaine polyédrique de R d , d t2, 3u, occupé par le milieu poreux déformable, et soit p0, t F q avec t F ¡ 0 l'intervalle de temps considéré.

LE PROBLÈME DE BIOT

Equilibre mécanique

Nous commençons par considérer la partie mécanique. Dans toute la suite, nous allons nous restreindre au cas des petites déformations sans forcément le rappeler à chaque fois. Etant donné une force volumique f : Ω ¢ p0, t F q Ñ R d (souvent f ρ ref F m , où ρ ref est la masse volumique homogénéisée et F m la force de gravité), la première équation du modèle de Biot décrit l'équilibre entre le tenseur de contraintes σ : Ω ¢ p0, t F q Ñ R d¢d sym tτ R d¢d ; τ T τu et cette force volumique :

∇ ¤ σ f 0.

(1.1)

Pour décrire le comportement mécanique d'un milieu poreux, l'hypothèse de Terzaghi [110] permet de décomposer le tenseur de contraintes comme la somme du tenseur de contraintes eectives σ I induites par les déformations ε : Ω ¢ p0, t F q Ñ R d¢d sym de la structure solide, et du tenseur de contraintes de pression σ p induites par la pression p : Ω ¢ p0, t F q Ñ R du uide :

σ σ I σ p .
De ce fait, σ est aussi appelé tenseur de contraintes totales. Les contraintes de pression dépendent linéairement et de manière isotrope de la pression. En introduisant le coecient de Biot-Willis b ¡ 0, elles s'écrivent σ p ¡bpI d .

Les contraintes eectives, quant à elles, sont exprimées en fonction de la loi de comportement mécanique. Dans le cas le plus simple de l'élasticité linéaire elles correspondent au comportement élastique formulé par la loi de Hooke σ I 2µε λ trpεqI, ε 1 2 p∇u ∇u T q,

(1.2) où u : Ω ¢ p0, t F q Ñ R d est le déplacement du squelette solide, et les paramètres de Lamé µ ¡ 0 et λ ¥ 0 décrivent les propriétés mécaniques du matériau.

Dans le Chapitre 3 nous considérons également deux autres lois de comportement décrivant des problèmes hyperélastiques : le modèle de Hencky-Mises et un modèle d'endommagement isotrope réversible. Les deux sont des variantes non linéaires de (1.2) ; dans le modèle de Hencky-Mises, les paramètres de Lamé ne sont plus constants, mais des fonctions scalaires de la partie déviatorique des déformations devpεq trpε 2 q ¡ 1 d trpεq 2 . Dans le modèle d'endommagement, on introduit une fonction d'endommagement D : R d¢d sym Ñ R et on multiplie l'expression de σ I dans (1.2) par 1 ¡ Dpεq.

Enn, nous traitons des lois élasto-plastiques dans le Chapitre 4. Ces lois combinent deux comportements diérents : la réaction élastique d'un matériau tant que les contraintes n'atteignent pas un certain critère ; et une déformation plastique, donc irréversible, si le critère est atteint. Nous développons des expressions pour la loi de comportement en fonction de ce cri-tère de charge en prenant également en compte des eets d'écrouissage, qui apparaissent si des déformations plastiques préalables ont modié ce critère par une modication de la structure microscopique du matériau. Nous considérons deux lois de comportement mécanique utilisées dans la mécanique de sols : d'une part la loi de Drucker-Prager [49], et d'autre part la loi L&K développée au sein d'EDF [46].

Conservation de la masse du uide

La deuxième équation du modèle de Biot exprime la conservation de la masse m : Ω¢p0, t F q Ñ R du uide en postulant l'équilibre entre la variation de la masse et le ux massique Φ :

Ω ¢ p0, t F q Ñ R d : f t m ∇ ¤ Φ 0.
(1.3)

La masse du uide est le produit de la masse volumique ρ du uide et de la porosité lagrangienne ϕ du squelette, qui correspond, dans le cas saturé, à la fraction du volume occupé par le uide, donc au volume du uide. La variation de la masse volumique du uide peut être exprimée en fonction de sa pression par f t ρ ρ f t p K w , où K w est la compressibilité du uide. Le changement du volume des pores est relié à la variation de la déformation du squelette par la relation suivante :

f t ϕ bf t ∇ ¤ u.

Pour la variation de la masse, on obtient donc f t m f t pρϕq ϕf t ρ ρf t ϕ ρpc 0 f t p bf t ∇ ¤ uq,

(1.4) où le paramètre c 0 ϕK ¡1 w mesure la quantité de uide qui peut être forcée dans l'espace poreux. Si le uide considéré est incompressible, on a K ¡1 w 0, et donc la variation de sa masse ne dépend pas de la pression, mais uniquement de la variation du volume des pores qu'il occupe. Nous supposons ici que le ux hydraulique φ : Ω ¢ p0, t F q Ñ R d suit la loi de Darcy avec une conductivité hydraulique κ : Ω Ñ R donnée. On peut donc l'écrire φppq ¡κ∇p.

La conductivité hydraulique est donnée par κ K

int K ¡1 visc , où K int : Ω Ñ R ¦ et K visc R ¦
dénotent respectivement la perméabilité intrinsèque et la viscosité du uide. Pour le ux massique on obtient donc Φ ρκp¡∇p ρF m q.

(1.5) Dans cette section nous décrivons de façon plus détaillée la stratégie de reconstruction de ux.

Nous expliquons ensuite comment ces développements ont été implémentés dans Code_Aster.

L'évolution dans l'implémentation, qui sera de plus en plus intégrée dans la structure de

Code_Aster au fur et à mesure des chapitres de la thèse, est explicitée dans chacune des sections correspondantes de l'introduction.

La stratégie de reconstruction des ux

Lorsqu'on discrétise en espace le problème (1.7) par la méthode des éléments nis de Lagrange (H 1 -conformes), les variables primales discrètes vérient de nouveau u h H 1 pΩq et p h H 1 pΩq à chaque étape de calcul dans l'algorithme de la Figure 1.3, tandis que les ux σpεpu h q, p h q et φpp h q calculés en fonction de

Introduction

Biot's poro-elasticity problem was originally proposed by von Terzaghi [110] and Biot [19] to describe the hydro-mechanical coupling between the displacement eld u of a linearly elastic, porous material and the pressure p of an incompressible, viscous uid saturating its pores.

Let Ω R 2 be the simply connected polygonal region occupied by the porous material, and let t F ¡ 0 denote the simulation time. For the sake of simplicity, we assume that the material is clamped at its impermeable boundary, and x the BiotWillis coecient equal to 1. We also assume that the deformation of the material is much slower than the ow rate, so that the problem can be considered in quasi-static form. Then, the displacement eld u : Ω ¢ p0, t F q Ñ R 2 and the pore pressure p : Ω ¢ p0,

t F q Ñ R are determined by ¡∇ ¤ σpuq ∇p f in Ω ¢ p0, t F q, (2.1a) f t p∇ ¤ u c 0 pq ¡ ∇ ¤ pκ∇pq g in Ω ¢ p0, t F q, (2.1b) u 0 on fΩ ¢ p0, t F q, (2.1c) κ∇p ¤ n Ω 0 on fΩ ¢ p0, t F q, (2.1d) up¤, 0q u 0 in Ω, (2.1e)
where f denotes the volumetric body force acting on the material, g a volumetric uid source (which, if c 0 0, is assumed to verify the compatibility condition ³ Ω gpx, tqdx 0 for each t p0, t F q), and the eective stress tensor σ is linked to the strain tensor through Hooke's law σpuq 2µ puq λ trp puqqI 2 , puq 1 2 p∇u ∇u T q, (2.2)

where I 2 is the two-dimensional identity matrix. The Lamé parameters λ and µ, describing the mechanical properties of the material, are assumed such that µ ¡ 0 and λ 2 3 µ ¡ 0 uniformly in Ω. The scalar eld κ : Ω Ñ R describes the mobility of the uid and we assume that there exist positive real numbers κ S and κ U such that κ S ¤ κ ¤ κ U a.e. in Ω. When the specic storage coecient c 0 is zero, we enforce uniqueness of the pore-pressure in (2.1) by further requiring that » Ω pp¤, tqdx 0 in p0, t F q.

(2.3)

On the other hand, for c 0 ¡ 0, we complement (2.1) by the following initial condition on the pressure:

pp¤, 0q p 0 in Ω.

(2.4)

In practice, even when c 0 0, the initial velocity eld u 0 in (2.1e) is usually obtained by rst setting p 0 equal to the solution of a hydrostatic computation and then calculating u 0 by solving (2.1a) with p p 0 (cf. Remark 2.1 in [85]). The well-posedness of Biot's consolidation problem has been analyzed in [98,111]. A suitable approximation method consists of using TaylorHood H 1 -conforming nite elements in space (using piecewise polynomials of order k ¥ 1 for the pressure and of order pk 1q for the displacement) and a backward Euler scheme in time. The corresponding a priori error analysis can be found in [7779]. This discretization strategy is adopted in the Code_Aster 1 software, which is used for the numerical examples presented in this work. Several other discretization methods have been studied in the literature, among which we cite, in particular, the fully coupled algorithm of [20], where the Hybrid High-Order method of [39] is used for the elasticity operator, while the weighted discontinuous Galerkin method of [40] is used for the Darcy operator.

The two governing equations (2.1a) and (2.1b) express, respectively, the conservation of mechanical momentum and uid mass. In particular, the Darcy velocity φppq : ¡κ∇p and the total stress tensor θpu, pq : σpuq ¡ pI 2 have continuous normal component across any interface in the domain Ω, and the divergence of these elds is locally in equilibrium with the sources (and the accumulation terms) in any control volume. It is well known that the use of H 1 -conforming nite elements does not lead to discrete uxes φpp n h q and θpu n h , p n h q (where pu n h , p n h q denotes the discrete solution at a given discrete time t n ) that satisfy the discrete counterpart of the above properties across mesh interfaces and in mesh cells. The rst contribution of this work is to ll this gap by reconstructing equilibrated uxes from local mixed nite element solves on cell patches around mesh vertices. The Darcy velocity reconstruction uses, as in [28,38,56], RaviartThomas mixed nite elements [92] on cell patches around vertices of the original mesh. The construction we propose for the total stress tensor is, to our knowledge, 1 http://web-code-aster.org 2.2. SETTING novel and is based on the use of the ArnoldWinther mixed nite element [10], again on the same vertex-based cell patches. This construction provides, in particular, a symmetric total stress tensor. The Darcy velocity and the total stress tensor are reconstructed at each discrete time, they have continuous normal component across any mesh interface, and their divergence is locally in equilibrium with the sources (averaged over the time interval) in any mesh cell.

In steady-state linear elasticity, element-wise (as opposed to patch-wise) reconstructions of equilibrated tractions from local Neumann problems can be found in [4,36,73,84], whereas direct prescription of the degrees of freedom in the ArnoldWinther nite element space is considered in [82].

The second contribution of this work is to perform an a posteriori error analysis of Biot's poro-elasticity problem using the above reconstructed uxes to compute the error indicators. Equilibrated-ux a posteriori error estimates for poro-elasticity appear to be a novel topic (residual-based error estimates can be found, e.g., in [52,76]). Equilibrated-ux a posteriori error estimates oer several advantages. On the one hand, error upper bounds are obtained with fully computable constants. The idea can be traced back to [89] and was advanced amongst others by [START_REF]A posteriori error estimate[END_REF]28,38,54,56,69,71,75,93]. Another interesting property is the polynomial-degree robustness proved recently for the Poisson problem in [27,56]. A third attractive feature introduced in [55] is to distinguish among various error components, e.g., discretization, linearization, and algebraic solver error components, and to equilibrate adaptively these components in the iterative solution of nonlinear problems. This idea was applied to multi-phase, multi-components (possibly non isothermal) Darcy ows in [4143]. For simplicity, we consider in the present work a global error measure which lends itself naturally to the development of equilibrated-ux error estimators, and dened as the dual energy-norm of the residual of the weak formulation. This paper is organized as follows. In Section 2.2, we introduce the weak and discrete formulations of Biot's poro-elasticity problem (2.1), along with some useful notation and preliminary results. In Section 2.3, we present the equilibrated reconstruction for the Darcy velocity and the total stress tensor. In Section 2.4, we derive a fully computable upper bound on the residual dual norm. We then distinguish two dierent error sources in the upper bound, namely the spatial and the temporal discretization, and we propose an algorithm adapting the mesh and the time step so as to equilibrate these error sources. Finally, we show numerical results in Section 2.5.

Setting

In this section we introduce some notation, the weak formulation, and the discrete solution of problem (2.1).

Weak formulation

We denote by L 2 pΩq, L 2 pΩq and L 2 pΩq the spaces composed of square-integrable functions taking values in R, R 2 and R 2¢2 respectively, and by p¤, ¤q and ¤ the corresponding inner product and norm. We also let L 2 0 pΩq : tq L 2 pΩq | pq, 1q 0u. H 1 pΩq stands for the Sobolev space composed of L 2 pΩq functions with weak gradients in L 2 pΩq and H 1 0 pΩq for its zero-trace subspace. Hpdiv, Ωq and Hpdiv, Ωq denote the spaces composed of L 2 pΩq and L 2 pΩq functions with weak divergence in L 2 pΩq and L 2 pΩq, respectively, H s pdiv, Ωq the subspace of Hpdiv, Ωq composed of symmetric-valued tensors, and H 0 pdiv, Ωq : tϕ Hpdiv, Ωq | ϕ ¤ n Ω 0 on fΩu.

We assume henceforth, for the sake of simplicity, that the volumetric body force f and the uid source g lie in L 2 p0, t F ; L 2 pΩqq and L 2 p0, t F ; L 2 0 pΩqq, respectively. In order to write a weak formulation of this poro-elastic problem, we dene

U : H 1 0 pΩq, P : H 1 pΩq, (2.5)
where in the case c 0 0 we require additionally that P H 1 pΩq L 2 0 pΩq, and introduce the following Bochner spaces:

X : L 2 p0, t F ; U q ¢ L 2 p0, t F ; P q, (2.6a)
Y : H 1 p0, t F ; U q ¢ H 1 p0, t F ; P q.

(2.6b)

Let u, v U and p, q P . We dene the bilinear forms apu, vq : pσpuq, pvqq,

bpv, qq : ¡pq, ∇ ¤ vq,

(2.7b) cpp, qq : pc 0 p, qq, (2.7c)
dpp, qq : pκ∇p, ∇qq.

(2.7d)

Then, we consider the following weak formulation: nd pu, pq Y , verifying the initial condition (2.1e) with u 0 H 1 0 pΩq and (2.4) with p 0 H 1 pΩq if c 0 ¡ 0, and such that, for a.e. t p0, t F q, apuptq, vq bpv, pptqq pfptq, vq dv U, (2.8a) ¡bpf t uptq, qq cpf t pptq, qq dppptq, qq pgptq, qq dq P.

(2.8b)

The well-posedness of Biot's consolidation problem in slightly dierent weak formulations is shown in [98,111]. The uniqueness of the solution to (2.8) can be shown by energy arguments.

Assuming the existence of the solution in Y , we denote by σpuq the resulting eective stress tensor, by θpu, pq σpuq ¡pI 2 the total stress tensor and by φppq ¡κ∇p the Darcy velocity.

SETTING

They verify the following properties: θpu, pq L 2 p0, t F ; H s pdiv, Ωqq, ¡∇ ¤ θpu, pq f,

(2.9a) φppq L 2 p0, t F ; H 0 pdiv, Ωqq, ∇ ¤ φppq g ¡ f t p∇ ¤ u c 0 pq. (2.9b)

Discrete setting

For the time discretization, we consider a sequence of discrete times pt n q 0¤n¤N such that t i t j whenever i j, t 0 0, and t N t F . For each 1 ¤ n ¤ N , let I n : pt n¡1 , t n q and τ n : t n ¡ t n¡1 . For a space-time function v, we denote v n : vp¤, t n q and dene the backward dierencing operator f n t v τ ¡1 n pv n ¡ v n¡1 q. At each time step 1 ¤ n ¤ N , the space discretization is based on a conforming triangulation T n h of Ω, i.e. a set of closed triangles with union equal to Ω and such that, for any distinct T 1 , T 2 T n h , the set T 1 T 2 is either a common edge, a vertex or the empty set. We assume that T n h veries the minimum angle condition, i.e., there exists α min ¡ 0 uniform with respect to all considered meshes such that the minimum angle α T of each triangle T T n h satises α T ¥ α min . The set of vertices of the mesh is denoted by V n h ; it is decomposed into interior vertices V n,int h and boundary vertices V n,ext h . For any subdomain ω Ω we denote V n ω the set of vertices in ω. For all a V n h , T n a is the patch of elements sharing the vertex a, and ω a the corresponding open subset of Ω. For all T T n h , V n T denotes the set of vertices of T , h T its diameter and n T its unit outward normal vector.

For all n N and all k N, we denote by P k pTq the space of bivariate polynomials in T T n h of total degree at most k and by P k pT n h q tϕ L 2 pΩq | ϕ| T P k pTq dT T n h u the corresponding broken space over T n h .

The following Poincaré's inequality holds for all T T n h :

v ¡ Π 0,T v T ¤ C P,T h T ∇v T dv H 1 pTq, (2.10)
where Π 0,T : L 1 pTq Ñ P 0 pTq is such that ³ T pv ¡ Π 0,T vqdx 0 and C P,T 1{π owing to the convexity of the mesh elements (see e.g. [15]). Let

RM : tb cpx 2 , ¡x 1 q T | b R 2 , c Ru (2.11)
denote the space of rigid body motions. We have the following Korn's inequality, again valid for all T T n h :

∇pv ¡ Π RM,T vq T ¤ C K,T pvq T dv H 1 pTq, (2.12)
where Π RM,T :

H 1 pTq Ñ RM is such that ³ T pv¡Π RM,T vqdx 0 and ³ T rotpv ¡Π RM,T vqdx 0 (with rotpvq : f x 1 v 2 ¡ f x 2 v 1 ),
T T n h , v ¡ Π RM,T v T ¤ h T π c 2 sinpα T {4q pvq T dv H 1 pTq.
(2.13)

Discrete problem

We will focus on the conforming TaylorHood nite element method using for each time step

1 ¤ n ¤ N the spaces U n h : P k 1 pT n h q U, P n h : P k pT n h q P, (2.14)
with k ¥ 1. This method was rst proposed in [102] for incompressible ows and is known to provide stable pore pressure approximations (cf.

[79] and [96,113]) and is a classical choice for the discretization of poro-mechanical problems by conforming nite elements.

Assumption 2.1 (Piecewise-constant-in-time source terms). For simplicity of exposition, we assume henceforth that the functions f and g are constant-in-time on each time interval I n and denote f n : f | In and g n : g| In .

Using the TaylorHood nite element spaces and a backward Euler scheme to march in time, the discrete problem reads: given u 0

h and, if c 0 ¡ 0, p 0 h , nd pu n h , p n h q U n h ¢ P n h , for all 1 ¤ n ¤ N , such that apu n h , v h q bpv h , p n h q pf n , v h q dv h U n h , (2.15a)
¡bpf t u n hτ , q h q cpf t p n hτ , q h q dpp n h , q h q pg n , q h q dq h P n h ,

where we denote by u hτ , p hτ the discrete space-time functions which are continuous and piecewise ane in time, and such that, for each 0 ¤ n ¤ N , pu hτ , p hτ qp¤,

t n q pu n h , p n h q, so that f t u n hτ : f t u hτ | In τ ¡1 n pu n h ¡ u n¡1 h q and f t p n hτ : f t p hτ | In τ ¡1 n pp n h ¡ p n¡1 h q.

Quasi-static ux reconstructions

In contrast to (2.9), we have in general θpu n h , p n h q H s pdiv, Ωq and φpp n h q H 0 pdiv, Ωq. In this section, we restore these properties by reconstructing Hpdivq-conforming discrete uxes.

These reconstructions are devised locally on patches of elements around mesh vertices. We rst present the reconstructions in an abstract setting; then we apply these reconstructions to Biot's poro-elasticity problem. Since the time variable is irrelevant in devising the reconstructions, we drop the index n in the abstract presentation.

QUASI-STATIC FLUX RECONSTRUCTIONS

pigure PFI ! he mixed virt!homs (nite element for l 0 @leftA nd l 1 @rightA

Darcy velocity

We reconstruct the Darcy velocity using mixed RaviartThomas nite elements of order l ¥ 0.

For each element T T h , the local RaviartThomas polynomial spaces are dened by W T : P l pTq xP l pTq,

(2.16a)

Q T P l pTq.
(2.16b) Figure 2.1 shows the corresponding degrees of freedom for l 0 and l 1. For each vertex a V h , the mixed RaviartThomas nite element spaces on the patch domain ω a are then dened as

W a h : tv h Hpdiv, ω a q | v h | T W T dT T a u, (2.17a)
Qa h : tq h L 2 pω a q | q h | T Q T dT T a u.
(2.17b)

We need to consider the following subspaces associated with the setting where a zero normal component is enforced on the velocity:

W a h : tv h W a h | v h ¤ n ωa 0 on fω a u, (2.18a) Q a h : tq h Qa h | pq h , 1q ωa 0u.
(2.18b)

The distribution of the degrees of freedom of functions in W a h is presented in Figure 2.2. Note that we enforce the zero normal condition also on patches associated with boundary vertices since a zero normal Darcy velocity is prescribed in the exact problem; see Remark 2.4 for other types of boundary conditions.

Construction 2.2 (Darcy velocity φ h ).

For each a V h , let γ a L 2 pω a q be such that pγ a , 1q ωa 0, and let Γ a L 2 pω a q. Consider the following constrained minimization problem:

ϕ a h argmin w h W a h , ∇¤w h Π Q a h γa w h ¡ Γ a ωa , ( 2
Q a h such that ϕ a h , w h ¨ωa ¡ ps a h , ∇ ¤ w h q ωa pΓ a , w h q ωa dw h W a h , (2.21a)
∇ ¤ ϕ a h , q h ¨ωa pγ a , q h q ωa dq h Q a h .

(

pγ h ¡ ∇ ¤ φ h , qq T 0 dq Q T dT T h . (2.22)
Remark 2.4 (Other boundary conditions). Suppose that we are given a partition of the boundary as fΩ fΩ N,P fΩ D,P (the subsets fΩ N,P and fΩ D,P are conventionally closed in fΩ, i.e., fΩ N,P fΩ D,P is the common boundary of the two subsets) and that an inhomogeneous Neumann condition is enforced on the ux φppq on fΩ N,P (and a Dirichlet condition is enforced on p in fΩ D,P ). Assume that the mesh is tted to the boundary partition, so that any mesh edge on the boundary belongs to either fΩ N,P or fΩ D,P . As detailed in [47], this situation can be accommodated in Construction 2.2 up to minor modications for all a V ext h (the construction is unmodied for all a V int h ). For the ux, we consider for the trial and test spaces, respectively,

W a h,N : tw h W a h | w h ¤ n ωa | fωazfΩ 0, w h ¤ n ωa | fωafΩ N,P Φ a,N u, W a h,0 : tw h W a h | w h ¤ n ωa | fωazfΩ 0, w h ¤ n ωa | fωafΩ N,P 0u,
with Φ a,N related to the Neumann condition, whereas we set Q a h : Qa h if a lies on some edge in fΩ D,P and Q a h as in (2.18b) otherwise.

QUASI-STATIC FLUX RECONSTRUCTIONS

pigure PFQ ! ilement digrms for the pir pΣT , VT q in the ses m 1 @leftA nd m 2 @rightA

Total stress tensor

We reconstruct the total stress tensor using mixed ArnoldWinther nite elements of order m ¥ 1. One advantage of using these elements is that the reconstructed stress tensor is symmetric. For each element T T h , the local ArnoldWinther polynomial spaces are dened by

Σ T : P s,m 1 pTq tτ P s,m 2 pTq | ∇ ¤ τ 0u (2.23a) tτ P s,m 2 pTq | ∇ ¤ τ P m pTqu, V T : P m pTq, (2.23b)
where P s,m pTq denotes the subspace of P m pTq composed of symmetric-valued tensors.

For each vertex a V h , the mixed ArnoldWinther nite element spaces on the patch domain ω a are dened as

Σa h : tτ h H s pdiv, ω a q | τ h | T Σ T dT T a u, (2.24a) Ṽ a h : tv h L 2 pω a q | v h | T V T dT T a u.
(2.24b) Figure 2.3 shows the corresponding degrees of freedom in the cases m 1 and m 2. The dimension of V T is pm 1qpm 2q, and it is shown in [10] that dimpΣ T q p3m 2 17m 28q{2. For the lowest-order case m 1, the 24 degrees of freedom in Σ T are

• The values of the three components of the (symmetric) stress tensor at each vertex of the triangle (9 dofs);

• The values of the moments of degree 0 and 1 of the normal components of the stress tensor on each edge (12 dofs);

• The value of the moment of degree 0 of each component of the stress tensor on the triangle (3 dofs).

We need to consider subspaces where a zero normal component is enforced on the stress tensor.

Since the boundary condition in the exact problem prescribes the displacement and not the normal stress, we distinguish the case whether a is an interior vertex or a boundary vertex a fΩ pigure PFR ! he degrees of freedom of the spe Σ a h in the se m 1 on pth for a V int h @leftA nd

a V ext h @rightA (see Remark 2.7 for other types of boundary conditions). For a V int h , we set Σ a h : tτ h Σa h | τ h n ωa 0 on fω a , τ h pbq 0 db V ωa fω a u, (2.25a) V a h : tv h Ṽ a h | pv h , zq ωa 0 dz RM u, (2.25b)
and for a V ext h , we set

Σ a h : tτ h Σa h | τ h n ωa 0 on fω a zfΩ, τ h pbq 0 db V ωa pfω a zfΩqu, (2.26a)
V a h : Ṽ a h .

(

Note that, as argued in [10], the nodal degrees of freedom on fω a are set to zero if the vertex separates two edges where the normal stress is enforced to be zero. The distribution of the degrees of freedom in Σ a h is presented in Figure 2.4.

Construction 2.5 (Total stress tensor θ h ). For each a V h , let λ a L 2 pω a q be such that pλ a , zq ωa 0 for all z RM and all a V int h , and let Λ a L 2 pω a q. Consider the following constrained minimization problem:

ϑ a h argmin τ h Σ a h , ∇¤τ h Π V a h λa τ h ¡ Λ a ωa , (2.27)
where Π V a h denotes the L 2 -orthogonal projection onto V a h . Then, extending ϑ a h by zero outside ω a , set θ h :

aV h ϑ a h .
(2.28)

The condition on λ a for all a V int h ensures that the constrained minimization problem (2.27) is well-posed. This problem is classically solved by nding ϑ a

h Σ a h and s a h V a h such that ϑ a h , τ h ¨ωa ¡ s a h , ∇ ¤ τ h ¨ωa Λ a , τ h ¨ωa dτ h Σ a h , (2.29a) ∇ ¤ ϑ a h , v h ¨ωa pλ a , v h q ωa dv h V a h . (2.29b)
This problem is well-posed (see [10]), and we obtain the following result:

Lemma 2.6 (Properties of θ h ). Let θ h be prescribed by Construction 2.5. Then, θ h H s pdiv, Ωq, and letting λ h L 2 pΩq be dened such that λ h | T °aV T λ a for all T T h , we have

pλ h ¡ ∇ ¤ θ h , vq T 0 dv V T dT T h .
(2.30)

Proof. All the elds ϑ a h are in H s pdiv, ω a q and satisfy appropriate zero normal conditions so that their zero-extension to Ω is in H s pdiv, Ωq. Hence, θ h H s pdiv, Ωq. Let us prove (2.30). Let a V int h . Since pλ a , zq ωa 0 for all z RM , we infer that (2.29b) actually holds for all

v h Ṽ a h . The same holds true if a V ext h by denition of V a h . Hence, ∇ ¤ ϑ a h , v h ¨ωa pλ a , v h q ωa for all v h Ṽ a
h and all a V h . Since Ṽ a h is composed of piecewise polynomials that can be chosen independently in each cell T T a , we conclude that (2.30) holds.

Remark 2.7 (Other boundary conditions). In the spirit of Remark 2.4 with the boundary partition fΩ fΩ N,U fΩ D,U , the minor modications of Construction 2.5 are as follows for all a V ext h (the construction is unmodied for all a V int h ): For the stress tensor, we consider for the trial and test spaces, respectively,

Σ a h,N : tτ h Σa h | τ h n ωa | fωazfΩ 0, τ h n ωa | fωafΩ N,U Θ a,N , τ h pbq 0 db V ωa pfω a zfΩq, τ h pbq θ a,N db V ωa pfΩzfΩ D,U qu, Σ a h,0 : tτ h Σa h | τ h n ωa | fωazfΩ 0, τ h n ωa | fωafΩ N,U 0, τ h pbq 0 db V ωa pfω a zfΩq, τ h pbq 0 db V ωa pfΩzfΩ D,U qu,
with Θ a,N and θ a,N related to the Neumann condition, whereas we set V a h as in (2.26b) if a lies on some edge in fΩ D,U and V a h as in (2.25b) otherwise.

Remark 2.8 (Extension to 3D

). The extension of Construction 2.5 to three dimensions hinges on the existence of mixed nite element spaces producing three dimensional, Hpdivqconforming, symmetric tensors. These were introduced in [8], but are complex to implement and require signicant computational eort, due to the high number of degrees of freedom per element (162 for the stress tensor).

Application to Biot's poro-elasticity problem

In this section, we apply the above constructions to the discrete Biot poro-elasticity problem (2.15). The reconstructed Darcy velocity and total stress tensor are space-time functions that are piecewise constant in time, i.e. these functions are calculated at every time step. We use Constructions 2.2 and 2.5 where we now specify the data γ a , Γ a , λ a , and Λ a for each 1 ¤ n ¤ N . For this purpose, we consider for any mesh vertex a V n h , the piecewise ane hat" function ψ a P 1 pT h q H 1 pΩq supported in ω a , which takes the value 1 at the vertex a and zero at the vertices lying on the boundary of ω a ; cf.

a V n T ψ a | T 1 dT T n h .
(2.31) Construction 2.9 (Darcy velocity and total stress reconstructions

). Let 1 ¤ n ¤ N . Dene for all a V n h , γ a ψ a g n ¡ ψ a f t p∇ ¤ u hτ c 0 p hτ q n ∇ψ a ¤ φpp n h q, Γ a ψ a φpp n h q, (2.32a) λ a ¡ψ a f n θpu n h , p n h q∇ψ a , Λ a ψ a θpu n h , p n h q, (2.32b)
where we recall that φpp n h q ¡κ∇p n h and θpu n h , p n h q σpu n h q ¡ p n h I 2 . Then dene φ n h H 0 pdiv, Ωq and θ n h H s pdiv, Ωq using Constructions 2.2 and 2.5, respectively, with l tk¡1, ku and m k, where k is the degree of the used TaylorHood element in (2.15). Remark 2.10 (Choice of l and m in Construction 2.9). The polynomial degree k in the Taylor Hood nite element method (2.14) corresponds to degree k for the pressure and degree k 1 for the displacement, implying polynomial degree k ¡ 1 for φpp h q and k for θpu h , p h q. Thus, it seems reasonable that the polynomial degree for the reconstruction of the velocity is one lower than for the reconstruction of the total stress tensor. It has been shown in [27,56] that k ¡1 or k are suitable choices for the velocity reconstruction, and for k 1 we can observe the expected convergence rates for l 0 and m 1 in the numerical test of Section 2.5.2. Lemma 2.11 (Darcy velocity and total stress reconstructions). Construction 2.9 is well dened, and the following holds:

p¡∇ ¤ θ n h , zq T pf n , zq T dT T n h , dz RM, (2.33a) p∇ ¤ φ n h , 1q T pg n ¡ f t p∇ ¤ u hτ c 0 p hτ q n , 1q T dT T n h .
(2.33b)

Proof. Construction 2.2 is well dened provided pγ a , 1q ωa 0 holds for all a V h , and this follows by taking q h ψ a in (2.15b) (this is possible since hat functions are contained in the discrete space for the pressure). Construction 2.5 is well dened provided pλ a , zq ωa 0 for all z RM and all a V n,int h , and this follows by taking v h ψ a z in (2.15a) (recall that k ¥ 1 in 2.4. A POSTERIORI ERROR ANALYSIS AND SPACE-TIME ADAPTIVITY (2.14), so that this choice is legitimate) and using that pθpu n h , p n h q, pψ a zqq pθpu n h , p n h q, ∇pψ a zq T q pθpu n h , p n h q, zp∇ψ a q T q pθpu n h , p n h q, ψ a p∇zq T q pθpu n h , p n h q, zp∇ψ a q T q.

Finally, the properties on the divergence follow from Lemmas 2.6 and 2.3 and the partition of unity (2.31) which implies that °aV n T ∇ψ a | T 0.

Remark 2.12 (Other boundary conditions). If inhomogeneous Neumann boundary conditions

φppq ¤ n Ω Φ N or θpu, pqn Ω Θ N are imposed on fΩ N,P and fΩ N,U , respectively, we modify the constructions using Remarks 2.4 and 2.7. In particular, in Remark

2.4, Φ a,N is the L 2 - projection of ψ a Φ N onto W a h ¤ n Ω , whereas in Remark 2.7, Θ a,N is the L 2 -projection of ψ a Θ N onto Σa h n Ω .

A posteriori error analysis and space-time adaptivity

In this section, we derive an a posteriori error estimate at every time step n based on the quasiequilibrated ux reconstructions of Section 2.3.3 for Biot's poro-elasticity problem. Using these estimators, we devise an adaptive algorithm including the adaptive choice of the mesh size and of the time step.

A posteriori error estimate

To derive the a posteriori error estimate, we consider the residual of the weak formulation (2.8). To combine the two equations in (2.8) into a single residual, the two equations must be written using the same physical units. Therefore, we introduce a reference time scale t and a reference length scale l which we will use as scaling parameters together with the Young modulus E µp3λ 2µqpλ µq ¡1 Bppu, pq, pv, qqqptq pf, vqptq t pg, qqptq dpv, qq X.

(2.35)

For any pair pu hτ , p hτ q Y , we dene the residual Rpu hτ , p hτ q X I of (2.35) as xRpu hτ , p hτ q, pv, qqy X I ,X :

» t F 0 Bppu ¡ u hτ , p ¡ p hτ q, pv, qqqptqdt.
Its dual norm is dened as Rpu hτ , p hτ q X I : sup pv,qqX, pv,qq X 1 xRpu hτ , p hτ q, pv, qqy X I ,X , with pv, qq 2 X :

» t F 0 pE pvq q 2 pl ∇q q 2 dt.

We rst derive a local-in-time a posteriori error estimate. Let the time step 1 ¤ n ¤ N be xed. Let us set

X n : L 2 pI n ; U q ¢ L 2 pI n ; P q and Y n : H 1 pI n ; U q ¢ H 1 pI n ; P q, and dene the norm ¤ Xn on X n in the same way as ¤ X on X.

e n P pqq : t » In f t p∇ ¤ u c 0 pq ¡ f t p∇ ¤ u hτ c 0 p hτ q, q ¨nptq ¡ φppq ¡ φpp hτ q, ∇q ¨ptqdt, (2.37b)
where we recall that I n pt n¡1 , t n q and that both u hτ and p hτ are continuous, piecewise ane functions in time.

For all 1 ¤ n ¤ N , let φ n h and θ n h be the constant in time elds over I n dened by Construction 2.9. For all T T n h , we dene the residual estimators η n R,T,U , η n R,T,P by η n R,T,U :

h T π c 2 sinpα T {4q E ¡1 f n ∇ ¤ θ n h T , (2.38a)
η n R,T,P :

h T π t l g n ¡ f n t p∇ ¤ u hτ c 0 p hτ q ¡ ∇ ¤ φ n h T , (2.38b)
where α T denotes the minimum angle of the triangle T , and the ux estimators η n F,T,U ptq,

η n F,
q, |T 1 ptq| § § § § § Ţ T n h pf n ∇ ¤ θ n h , pv ¡ Π K,T vqptqq T § § § § § ¤ Ţ T n h η n R,T,U E pvqptq T ,
where we have used (2.33a) to insert Π K,T v inside the integral and (2.13) to conclude. For the second term, using the CauchySchwarz inequality readily yields

|T 2 ptq| ¤ Ţ T n h θ n h ¡ θpu hτ , p hτ qptq T pvqptq T Ţ T n h η n F,T,U ptqE pvqptq T .
Inserting these results into (2.41) and applying the Cauchy-Schwarz inequality yields

|e U pvq| ¤ ¤ ¥ » In Ţ T n h pη n R,T,U η n F,T,U ptqq 2 dt 1 {2 ¢ ¢» In pE pvqptq q 2 dt 1 {2
.

Proceeding in a similar way for e n P using (2.33b) and Poincaré's inequality (2.10) in place of (2.33a) and (2.13), respectively, we obtain

|e P pqq| ¤ ¤ ¥ » In Ţ T n h pη n R,T,P η n F,T,P ptqq 2 dt 1 {2 ¢ ¢» In pl ∇qptq q 2 dt 1 {2
.

Combining these results and using again the Cauchy-Schwarz inequality yields

|e U pvq e P pqq| ¤ ¤ ¥ » In Ţ T n h pη n R,T,U η n F,T,U ptqq 2 pη n R,T,P ptq η n F,T,P ptqq 2 dt 1 {2
¢ pv, qq Xn , and passing to the supremum concludes the proof.

Remark 2.14 (Data oscillation). Lemmas 2.3 and 2.6, and the mixed nite element space

property ∇ ¤ Σ n h pTq V n h pTq and ∇ ¤ W n h pTq Q n h pTq for any T T n h imply η n R,T,U h T π c 2 sinpα T {4q E ¡1 f n ¡ Π V n h pTq f T , η n R,T,P h T π t l g n ¡ f n t p∇ ¤ u hτ c 0 p hτ q ¡ Π Q n h pTq pg n ¡ f n t p∇ ¤ u hτ c 0 p hτ qq T .
For the sake of convenience, we assumed the source terms f and g to be piecewise constant in time.

η n N,T,U ȩE N,U T h T p2h e C t q 1 {2 E sinpα T {4q|T| 1 {2 θ n h n Ω ¡ θ N e , η n N,T,P ȩE N,P T t h T ph e C t q 1 {2 l |T| 1 {2 φ n h ¤ n Ω ¡ φ N e ,
where C t 0.77708, and (2.40) now reads

e n ¤ ¤ ¥ » In Ţ T n h tpη n R,T,U η n F,T,U ptq η n N,T,U q 2 pη n R,T,P η n F,T,P ptq η n N,T,P q 2 udt 1 {2
.

To dene a global-in-time a posteriori error estimate, we additionally dene the initial condition errors η IC,T,U and η IC,T,P by setting η IC,T,U :

¢ 1 2 E ¡1 t σpu 0 ¡ u hτ p¤, 0qq, pu 0 ¡ u hτ p¤, 0qq ¨T 1 {2 , (2.42a)
η IC,T,P :

¢ 1 2 E ¡1 pt q 2 c 0 ppp 0 ¡ p hτ p¤, 0qq, p 0 ¡ p hτ p¤, 0qq T 1 {2 , (2.42b)
and we set

e IC η IC ¤ ¥ Ţ T 0 h η 2 IC,T,U η 2 IC,T,P 1 {2 . (2.43)
We dene the global error as e : Rpu hτ , p hτ q X I e IC .

(

e ¤ ¤ ¥ N ņ1 » In Ţ T n h tpη n R,T,U η n F,T,U ptqq 2 pη n R,T,P η n F,T,P ptqq 2 udt 1 {2 η IC .
(2.45)

Proof. For each 1 ¤ n ¤ N , let ξ n X be the Riesz-representative of J n : X n Ñ R with J n pv, qq

³

In Bppu ¡ u hτ , p ¡ p hτ q, pv, qqqdt. Then the function ξ X dened by ξ |In : ξ n will be the Riesz-representative of J :

X Ñ R with pv, qq Þ Ñ ³ t F 0 Bppu ¡ u hτ , p ¡ p hτ q, pv, qqqdt, so that Rpu hτ , p hτ q 2 X I J 2 X I ξ 2 X N ņ1 ξ n 2 Xn N ņ1 J n 2 X I n N ņ1 pe n q 2 .
(2.46)

Inserting this result into (2.44) and applying Theorem 2.13 concludes the proof.

Distinguishing the space and time error components

The goal of this section is to elaborate the error estimate (2.40) so as to distinguish the error components resulting from the spatial and the temporal discretization. This is essential for the development of Algorithm 2.18 below, where the space mesh and the time step are chosen adaptively. Therefore, we add and subtract the discrete uxes in the ux estimators (2.39) and apply the triangle inequality. We obtain, for all T T n h , the following local spatial and temporal discretization error estimators :

η n sp,T,U : η n R,T,U E ¡1 θ n h ¡ θpu n h , p n h q T , (2.47a)
η n sp,T,P :

η n R,T,P t l φ n h ¡ φpp n h q T , (2.47b) η n tm,T,U ptq : E ¡1 θpu n h , p n h q ¡ θpu hτ , p hτ qptq T , (2.47c)
η n tm,T,P ptq : t l φpp n h q ¡ φpp hτ qptq T .

(2.47d)

For each of these local estimators we can dene a global version by setting

η n ,tU,Pu : ¤ ¥ 2 » In Ţ T n h ¡ η n ,T,tU,Pu ptq © 2 dt 1 {2 . (2.48)
Inserting them into (2.40) and applying the triangle inequality yields the following result.

Theorem 2.17 (A posteriori error estimate distinguishing the error components). Let 1 ¤ n ¤ N . Let pu, pq be the weak solution of (2.8) and let pu hτ , p hτ q Y n be the discrete solution of (2.15). Let θ n h , φ n h be the equilibrated uxes of Construction 2.9. (2.49)

Adaptive algorithm

Based on the error estimate of Theorem 2.17, we propose an adaptive algorithm where the mesh size and time step are locally adapted. The idea is to compare the estimators for the two error sources with each other in order to concentrate the computational eort on reducing the dominant one. Thus, both the spatial mesh and the time step are adjusted until space and time discretization contribute nearly equally to the overall error. For this purpose, let Γ tm ¡ 1 ¡ γ tm ¡ 0 be user-given weights and crit n , for all 0 ¤ n ¤ N , a chosen threshold that the error on the time interval I n should not exceed. For each of the considered error sources we dene the corresponding estimator as η : η ,U η ,P , so that (2.49) becomes e n ¤ η n sp η n tm .

(2.50) Algorithm 2.18 (Adaptive algorithm).

1. Initialisation (a) Choose an initial triangulation T 0 h , an initial time step τ 0 , and set t 0 : 0 (b) Initial mesh adaptation loop 1 {2 crit 0 .

(2.52)

In order to keep computational costs in the algorithm low, the initial mesh and time step should be chosen in a way that they match the criteria crit 0 and crit 1 . This can be achieved by performing only one time step before running the whole computation, and by modifying the initial discretization if they do not.

Numerical results

In this section we illustrate numerically our theoretical results on four test cases. For all tests we use the TaylorHood nite elements (2.14) with k 1 and Construction 2.9 with l 0 and m 1. In the rst two test cases, analytical solutions are known; the rst one is a purely elastic, stationary problem and the second one a Biot's poro-elasticity problem. We analyze the convergence rates of the error estimators and compare them to those of an energy-type norm of the analytical error. The third test is the quarter ve-spot problem, where we compare the results of the adaptive algorithm to a standard solution with xed mesh and time steps.

In the fourth test, the excavation of two parallel tunnels is simulated. It shows an industrial application of the error estimators used for remeshing and again compares the performance of Algorithm 2.18 to a standard resolution.

Purely mechanical analytical test

For this stationary, purely mechanical test we consider the mode I loading of a cracked plate, corresponding to pure tension at the top and the bottom applied at the innity. Following [112],

an analytical solution around the crack tip is given by upr, θq

1 ν E c 2π c r £ ¡ cosp θ 2 qp3 ¡ 4ν ¡ cospθqq sinp θ 2 qp3 ¡ 4ν ¡ cospθqq , (2.53)
pigure PFU ! irror estimtion @leftA nd nlytil error @rightA on n initil mesh nd fter three mesh re(nements leading to a singularity of the stress tensor at the crack tip. For our test, we restrain ourselves to the domain Ω p¡ 1 2 , 1 2 q ¢ p¡ 1 2 , 1 2 q with a straight crack from p¡ 1 2 , 0q (cf. Figure 2.6), and impose the analytical solution (2.53) as Dirichlet boundary condition on fΩ and the crack faces to obtain the discrete solution u h . The Young modulus and the Poisson ratio are set to E 1 and ν 0, leading to the Lamé parameters µ 0.5 and λ 0. Since in this purely mechanical test case there is no need for nondimensionalization, we omit the scaling factor E ¡1 in the error estimators. Figure 2.7 compares the distribution of the error estimators and the analytical error measured in the energy norm u ¡ u h en apu ¡ u h , u ¡ u h q 1 {2 . Besides detecting the dominating error at the crack tip due to the singularity of σpuq, the error estimators reect the distribution of the analytical error in the whole domain, as can be seen in the lower panel for the ner mesh.

Poro-elastic analytical test

Let Ω p0, 1q ¢ p0, 1q.

Quarter ve-spot problem

In this standard conguration considered in petroleum engineering, the injection of water at the center of a square domain and the production at the four corners is simulated on a quarter of the domain. In our test, this quarter is a square of 100m side length, divided into two parts with dierent mobilities; a circle around the injection point of radius 50m with κ 8 ¤ 10 ¡9 m 2 Pa ¡1 s ¡1 , and κ 10 ¡9 m 2 Pa ¡1 s ¡1 in the rest of the domain. The Young modulus and the Poisson ratio are given by E 10 9 Pa, ν 0.3, and we set c 0 0. The initial state is given by θ 0 0, φ 0 0 and p 0 10 5 Pa. le PFQ ! he three omputtions in our test for the qurter (veEspot prolem pigure PFW ! snitil mesh @leftA nd meshes t the end of the (rst @enterA nd the seond @rightA exvtion in the dptive lgorithm

We use Algorithm 2.18 to perform space-time adaptivity. We start with an initial mesh of 10,638 vertices and with an initial time step of τ 0 12h. For the space-time error balancing, we set γ tm 0.8 and Γ tm 1.3 and x the error limit for each time step to crit n 0.005τ n .

We compare the performance of the adaptive algorithm to two static computations (i.e. with xed meshes and time steps), one, called equivalent, where the discretization is chosen in a way to have approximately the same number of space-time unknowns as in the adaptive algorithm, and one where the discretization is very ne, so its solution can be taken as a reference solution. Table 2.3 compares the number of space-time unknowns and performed iterations (i.e. the number of time steps, counting repetitions in the adaptive algorithm), and the values of the error estimators of the three computations.

The left graphic in Figure 2.8 shows the discrete pressure along the diagonal going from the bottom left to the top right of the domain (as indicated in the right graphic) at three dierent times obtained by the static computations (solid and dotted lines) and the adaptive algorithm (dashed lines). The loosely dotted vertical line marks the edge between the two parts of Ω with dierent permeabilities. At each of these times, the discrete solution of the adaptive algorithm is closer to the reference solution than the equivalent computation using a xed mesh and time step. At the last time step, all the results get closer as the solution converges in time to a constant state.

Excavation damage test

In the context of the conception of a radioactive waste repository site, the excavation of tunnels destined to contain waste packages is numerically simulated. The domain Ω is a 80m ¢ 60m uxes are given by θ 0 θ ref and φ 0 0, while the initial pressure is p 0 . The parameters describing the rock are the Young modulus E 5800MPa, the Poisson ratio ν 0.3, the specic storage coecient c 0 0, and the hydraulic mobility κ 10 ¡13 m 2 Pa ¡1 s ¡1 . For the nondimensionalization of the problem, we used, along with E, the parameters t 1h and l 100m.

The performance of Algorithm 2.18 is tested on four dierent initial meshes with crit n 7 ¤ 10 ¡3 τ n for the coarsest one and, with the mesh getting ner, crit n 4¤10 ¡3 τ n , crit n 2¤10 ¡3 τ n and crit n 1 ¤ 10 ¡3 τ n . In all the calculations we x γ tm 0.8, Γ tm 1.5 and τ 0 3.9d. Figure 2.9 illustrates the evolution of the second coarsest mesh with crit n 4 ¤10 ¡3 τ n . During the rst excavation, the renement takes only place around the left tunnel, whereas the area around the right tunnel is only rened after the beginning of the second excavation. The calculations resulting from the adaptive algorithm are compared to calculations with xed meshes and time steps. Each of these meshes is slightly ner around the tunnels than in the rest of Ω, and the time steps are chosen in a way that η sp η tm . Figure 2.10 compares the spatial discretization error estimators at the nal time t F of the static algorithm to those of the adaptive algorithm, which are much more evenly distributed over the domain. Furthermore, the left graphic in Figure 2.11 shows that in our test, the use of the adaptive algorithm reduces the number of space-time-unknowns for a similar value of the error estimator.

In the right graphic of (or both). Thus, only the square shaped points in the graphic contribute to the overall error estimate. In the consolidation phase between the two excavations (from t 17.4 days to t 29 days), the mesh is slightly coarsened and the time step considerably increased, since the dominating error source in this phase is the spatial discretization.

Conclusion

The analytical test cases show that the distribution and convergence rates of our error estimators reect those of the analytical error. The eciency of Algorithm 2.18 has been illustrated in industrial tests, where the number of space-time unknowns is considerably decreased for a comparable overall error estimate. We also observe that the price for computing the ux reconstructions can be substantially reduced by pre-processing, a task that is fully parallelizable.

As shown in the rst test, the stress reconstruction and a posteriori estimate presented in this work are directly applicable to pure linear elasticity problems. The second test shows that the presented error estimate also delivers sharp bounds (as reected by moderate eectivity indices) of more accessible error measures computed using energy-type norms. In the third and fourth tests, comparing the proportions of the estimators for the hydraulic and the mechanical parts reects the physical properties of the problem: in the quarter ve-spot test, the dominating estimators are those for the hydraulic part; for the excavation damage test, they are approximately of the same order of magnitude, with the mechanical estimator dominating in regions of stress concentration.

Abstract

We consider hyperelastic problems and their numerical solution using a conforming nite element discretization and some iterative linearization algorithm. For these problems, we present equilibrated, weakly symmetric, Hpdivq-conforming stress tensor reconstructions, obtained from local problems on patches around vertices using the ArnoldFalkWinther nite element spaces. We distinguish two stress reconstructions, one for the discrete stress and one representing the linearization error.

The reconstructions are independent of the mechanical behavior law. Based on these stress tensor reconstructions we derive an a posteriori error estimate distinguishing the discretization, linearization, and quadrature error estimates, and propose an adaptive algorithm balancing these dierent error sources. We prove the eciency of the estimate, and conrm it on a numerical test with analytical solution for the linear elasticity problem. We then apply the adaptive algorithm to a more applicationoriented test, considering the Hencky-Mises and an isotropic damage model.

Introduction

In this work we develop equilibrated Hpdivq-conforming stress tensor reconstructions for a class of (linear and) nonlinear elasticity problems in the small deformation regime. Based on these reconstructions, we can derive an a posteriori error estimate distinguishing the discretization and linearization errors for conforming discretizations of the problem.

Let Ω R d , d t2, 3u, be a bounded, simply connected polyhedron, which is occupied by a body subjected to a volumetric force eld f L 2 pΩq. For the sake of simplicity, we assume that the body is xed on its boundary fΩ. The nonlinear elasticity problem consists in nding a vector-valued displacement eld u :

Ω Ñ R d solving ¡∇ ¤ σp∇ s uq f in Ω, (3.1a)
u 0 on fΩ,

where ∇ s u 1 2 pp∇uq T ∇uq denotes the symmetric gradient and expresses the strain tensor associated to u. In these applications, the solution is often approximated using H 1 -conforming nite elements. For nonlinear mechanical behavior laws, the resulting discrete nonlinear equation can then be solved using an iterative linearization algorithm yielding at each iteration a linear algebraic system to be solved, until the residual of the nonlinear equation lies under a predened threshold.

In this paper we develop an a posteriori error estimate allowing to distinguish between the CHAPTER 3. NONLINEAR ELASTICITY PROBLEM 51 error stemming from the linearization of the problem and the one due to its discretization, as proposed in [55] for nonlinear diusion problems. Thanks to this distinction we can, at each iteration, compare these two error contributions and stop the linearization algorithm once its contribution is negligible compared to the discretization error.

The a posteriori error estimate is based on equilibrated stress reconstructions. It is well known that, in contrast to the analytical solution, the discrete stress tensor resulting from the conforming nite element method does not have continuous normal components across mesh interfaces, and that its divergence is not locally in equilibrium with the source term f on mesh elements. In this paper we consider the stress tensor reconstruction proposed in [94] for linear elasticity to restore these two properties. This reconstruction uses the ArnoldFalk (except for the constant in the upper bound). Therefore, its implementation is independent and directly applicable to these laws, which makes the method convenient for FEM softwares in solid mechanics, which often provide a large choice of behavior laws. In addition, equilibrated error estimates were proven to be polynomial-degree robust for several linear problems in 2D, This paper is organized as follows. In Section 3.2 we rst formulate the assumptions on the stress-strain function σ and provide three examples of models used in the engineering practice.

We then introduce the weak and the discrete formulations of problem (3.1) and its linearization, along with some useful notation. In Section 3.3 we present the equilibrated stress tensor reconstructions, rst assuming that we solve the nonlinear discrete problem exactly, and then, based on this rst reconstruction, distinguish its discrete and its linearization error part at each iteration of a linearization solver. In Section 3.4 we derive the a posteriori error estimate, again rst assuming the exact solution of the discrete problem and then distinguishing the different error components. We then propose an algorithm equilibrating the error sources using adaptive stopping criteria for the linearization and adaptive remeshing. We nally show the eciency of the error estimate. In Section 3.5 we evaluate the performance of the estimates for the three behavior laws given as examples on numerical test cases.

Setting

In this section we will give three examples of hyperelastic behavior laws, before writing the weak and the considered discrete formulation of problem (3.1).

Continuous setting

Assumption 3.1 (Stress-strain relation). We assume that the symmetric stress tensor σ :

R d¢d sym Ñ R d¢d sym is continuous on R d¢d
sym and that σp0q 0. Moreover, we assume that there exist real numbers C gro , C mon p0, Vq such that, for all τ, η R d¢d sym ,

|σpτq| d¢d ¤ C gro |τ| d¢d , (growth) (3.2a) ¡ σpτq ¡ σpηq © : ¡ τ ¡ η © ¥ C 2 mon |τ ¡ η| 2 d¢d , (strong monotonicity) (3.2b)
where τ : η : trpτ T ηq with trpτq : °d i1 τ ii , and |τ| 2 d¢d τ : τ.

We next discuss a number of meaningful stress-strain relations for hyperelastic materials that Before presenting the variational formulation of problem (3.1), some useful notations are introduced. For X Ω, we respectively denote by p¤, ¤q X and ¤ X the standard inner product and norm in L 2 pXq, with the convention that the subscript is omitted whenever X Ω. The same notation is used in the vector-and tensor-valued cases. H 1 pΩq and Hpdiv, Ωq stand for the Sobolev spaces composed of vector-valued L 2 pΩq functions with weak gradient in L 2 pΩq, and tensor-valued L 2 pΩq functions with weak divergence in L 2 pΩq, respectively. Multiplying equation (3.1a) by a test function v H 1 0 pΩq and integrating by parts one has pσp∇ s uq, ∇ s vq pf, vq. Given f L 2 pΩq, nd u H 1 0 pΩq s.t., dv H 1 0 pΩq, apu, vq pf, vq.

(3.13) From (3.13) it is clear that the analytical stress tensor σp∇ s uq lies in the space H s pdiv, Ωq : tτ L 2 pΩq | ∇ ¤ τ L 2 pΩq and τ is symmetricu.

Discrete setting

The discretization (3.13) is based on a conforming triangulation T h of Ω, i.e. a set of closed triangles or tetrahedra with union Ω and such that, for any distinct T 1 , T 2 T h , the set T 1 T 2 is either a common edge, a vertex, the empty set or, if d 3, a common face. We assume that T h veries the minimum angle condition, i.e., there exists α min ¡ 0 uniform with respect to all considered meshes such that the minimum angle α T of each T T h satises α T ¥ α min .

The set of vertices of the mesh is denoted by V h ; it is decomposed into interior vertices V int h and boundary vertices V ext h . For all a V h , T a is the patch of elements sharing the vertex a, ω a the corresponding open subdomain in Ω and V a the set of vertices in ω a . For all T T h , V T denotes the set of vertices of T , h T its diameter and n T its unit outward normal vector. For all p N and all T T h , we denote by P p pTq the space of d-variate polynomials in T of total degree at most p and by P p pT h q tϕ L 2 pΩq | ϕ |T P p pTq dT T h u the corresponding broken space over T h . In the same way we denote by P p pTq and P p pTq, respectively, the space of vector-valued and tensor-valued polynomials of total degree p over T , and by P p pT h q and P p pT h q the corresponding broken spaces over T h .

In this work we will focus on conforming discretizations of problem (3.11) of polynomial degree p ¥ 2 to avoid numerical locking, cf [107]. The discrete formulation reads: nd u h H 1 0 pΩq P p pT h q such that dv h H 1 0 pΩq P p pT h q, apu h , v h q pf, v h q.

(3.14

)
This problem is usually solved using some iterative linearization algorithm dening at each iteration k ¥ 1 a linear approximation σ k¡1 of σ. Then the linearized formulation reads: nd u k h H 1 0 pΩq P p pT h q such that dv h H 1 0 pΩq P p pT h q, pσ k¡1 p∇ s u k h q, ∇ s v h q pf, v h q.

(3.15)

For the Newton algorithm the linearized stress tensor is dened as

σ k¡1 p∇ s u k h q : fσpτq fτ | τ∇s u k¡1 h ∇ s pu k h ¡ u k¡1 h q σp∇ s u k¡1 h q.
(3.16)

Equilibrated stress reconstruction

In general, the discrete stress tensor σp∇ s u h q resulting from (3.14) does not lie in Hpdiv, Ωq and thus cannot verify the equilibrium equation (3.1a). In this section we will reconstruct from σp∇ s u h q a discrete stress tensor σ h satisfying these properties. Based on this reconstruction, we then devise two equilibrated stress tensors representing the discrete stress and the linearization error respectively, which will be useful for the distinction of error components in the a posteriori error estimate of Section 3.4.2. Σ T : P q pTq, V T : P q¡1 pTq, Λ T : tµ P q¡1 pTq | µ ¡µ T u.

Patchwise construction in the

For q 2, the degrees of freedom are displayed in Figure 3.1. On a patch ω a the global space Σ h pω a q is the subspace of Hpdiv, ω a q composed of functions belonging piecewise to Σ T . The spaces V h pω a q and Λ h pω a q consist of functions lying piecewise in V T and Λ T respectively, with no continuity conditions between two elements.

Let now q : p. On each patch we need to consider subspaces where a zero normal component is enforced on the stress tensor on the boundary of the patch, so that the sum of the local solutions will have continous normal component across any mesh face inside Ω. Since the boundary condition in the exact problem prescribes the displacement and not the normal stress, we distinguish the case whether a is an interior vertex or a boundary vertex. If a V int h we set

Σ a h : tτ h Σ h pω a q | τ h n ωa 0 on fω a u, (3.18a) V a h : tv h V h pω a q | pv h , zq ωa 0 dz RM d u, (3.18b) Λ a h : Λ h pω a q, ( 3.18c)
where

RM 2 : tb cpx 2 , ¡x 1 q T | b R 2 , c Ru and RM 3 : tb a ¢ x | b R 3 , a
V a h : V h pω a q, ( 3.18e)
Λ a h : Λ h pω a q.

(3.18f )

For each vertex a V h we dene its hat function ψ a P 1 pT h q as the piecewise linear function taking value one at the vertex a and zero on all other mesh vertices.

Construction 3.5 (Stress tensor reconstruction). Let u h solve (3.14). For each

a V h nd pσ a h , r a h , λ a h q Σ a h ¢ V a h ¢ Λ a h such that for all pτ h , v h , µ h q Σ a h ¢ V a h ¢ Λ a h ,
pσ a h , τ h q ωa pr a h , ∇ ¤ τ h q ωa pλ a h , τ h q ωa pψ a σp∇ s u h q, τ h q ωa , (

p∇ ¤ σ a h , v h q ωa p¡ψ a f σp∇ s u h q∇ψ a , v h q ωa , (3.19b) pσ a h , µ h q ωa 0. Lemma 3.6 (Properties of σ h ). Let σ h be prescribed by Construction 3.5. Then σ h Hpdiv, Ωq, and for all T T h , the following holds: pf ∇ ¤ σ h , vq T 0 dv V T dT T h .

(3.21)

Proof. All the elds σ a h are in Hpdiv, ω a q and satisfy appropriate zero normal conditions so that their zero-extension to Ω is in Hpdiv, Ωq. Hence, σ h Hpdiv, Ωq. Let us prove (3.21). Since (3.20) holds for all a V int h , we infer that (3.19b) is actually true for all v h V h pω a q. The same holds if a V ext h by denition of V a h . Since V h pω a q is composed of piecewise polynomials that can be chosen independently in each cell T T a , and using σ h | T °aV T σ a h | T and the partition of unity °aV T ψ a 1, we infer that pf ∇ ¤ σ h , vq T 0 for all v V T and all T T h .

Discretization and linearization error stress reconstructions

Let now, for k ¥ 1, u k h solve (3.15). We will construct two dierent equilibrated Hpdivqconforming stress tensors. The rst one, σ k h,disc , represents as above the discrete stress tensor σp∇ s u k h q, for which we will have to modify Construction 3.5, because the Neumann compatibility condition (3.20) is not satised anymore. The second stress tensor σ k h,lin will be a measure for the linearization error and approximate σ k¡1 p∇ s u k h q ¡ σp∇ s u k h q. The matrix resulting from the left side of (3.19) will stay unchanged and we will only modify the source terms.

We denote by σp∇ s u k h q the L 2 -orthogonal projection of σp∇ s u k h q onto P p¡1 pT h q such that pσp∇ s u k h q ¡ σp∇ s u k h q, τ h q 0 for any τ h P p¡1 pT h q. The so obtained problem reads:

nd pσ a h , r a h , λ a h q Σ a h ¢ V a h ¢ Λ a h such that for all pτ h , v h , µ h q Σ a h ¢ V a h ¢ Λ a h ,
pσ a h , τ h q ωa pr a h , ∇ ¤ τ h q ωa pλ a h , τ h q ωa pψ a σp∇ s u k h q, τ h q ωa , p∇ ¤ σ a h , v h q ωa p¡ψ a f σp∇ s u k h q∇ψ a ¡ y k disc , v h q ωa , pσ a h , µ h q ωa 0.

Σ a h ¢ V a h ¢ Λ a h such that for all pτ h , v h , µ h q Σ a h ¢ V a h ¢ Λ a h ,
pσ a h , τ h q ωa pr a h , ∇ ¤ τ h q ωa pλ a h , τ h q ωa pψ a pσ k¡1 p∇ s u k h q ¡ σp∇ s u k h qq, τ h q ωa , p∇ ¤ σ a h , v h q ωa ppσ k¡1 p∇ s u k h q ¡ σp∇ s u k h qq∇ψ a y k disc , v h q ωa , pσ a h , µ h q ωa 0.

Then set σ k h,lin : °aV h σ a h .

Notice that the role of y k disc is to guarantee that, for interior vertices, the source terms in Constructions 3.7 and 3.8 satisfy the Neumann compatibility conditions p¡ψ a f σp∇ s u k h q∇ψ a ¡ y k disc , zq ωa 0 dz RM d , ppσ k¡1 p∇ s u k h q ¡ σp∇ s u k h qq∇ψ a y k disc , zq ωa 0 dz RM d .

Lemma 3.9 (Properties of the discretization and linearization error stress reconstructions).

Let σ k h,disc and σ k h,lin be prescribed by Constructions 3.7 and 3.8. Then it holds

1. σ k h,disc , σ k h,lin Hpdiv, Ωq, 2. pf ∇ ¤ pσ k h,disc σ k h,lin q, vq T 0 dv V T dT T h , 3.
As the Newton solver converges, σ k h,lin Ñ 0.

Proof. The proof is similar to the proof of Lemma 3.6. The rst property is again satised due to the denition of Σ a h . In order to show that the second property holds, we add the two equations (3.19b) obtained for each of the constructions. The right hand side of this sum then reads p¡ψ a f σ k¡1 p∇ s u k h q∇ψ a , v h q ωa . Once again we can, for any z RM d , take ψ a z as a test function in (3.15) to show that this term is zero if v h RM d , and so the equation holds for all v h V h pω a q. Then we proceed as in the proof of Lemma 3.6.

A posteriori error estimate and adaptive algorithm

In this section we rst derive an upper bound on the error between the analytical solution of (3.13) and the solution u h of (3.14), in which we then identify and distinguish the discretization and linearization error components at each Newton iteration for the solution u k h of (3.15).

Based on this distinction, we present an adaptive algorithm stopping the Newton iterations once the linearization error estimate is dominated by the estimate of the discretization error.

Finally, in a more theoretical part, we show the eectivity of the error estimate.

Guaranteed upper bound

We measure the error in the energy norm v 2 en : apv, vq pσp∇ s vq, ∇ s vq, (3.25) for which we obtain the properties and u h the discrete solution of (3.14). Let σ h be the stress tensor dened in Construction 3.5.

C 2 mon C ¡2 K ∇v 2 ¤ v 2 en ¤ C gro ∇ s v 2 , ( 3
Then, Proof of Theorem 3.10. We start by bounding the energy norm of the error by the dual norm of the residual of the weak formulation (3.13). Using (3.26), (3.2b), the linearity of a in its second argument, and (3.13) we obtain

u ¡ u h en ¤ c 2C gro C ¡3 mon £ Ţ T h h T π f ∇ ¤ σ h T σ h ¡ σp∇ s u h q T ¨2 1 {2 . ( 3
u ¡ u h 2 en ¤ C gro ∇ s pu ¡ u h q 2 ¤ C gro C ¡2 mon |apu, u ¡ u h q ¡ apu h , u ¡ u h q| C gro C ¡2 mon ∇pu ¡ u h q § § § § § a £ u, u ¡ u h ∇pu ¡ u h q ¡ a £ u h , u ¡ u h ∇pu ¡ u h q § § § § § ¤ C gro C ¡3 mon C K u ¡ u h en sup vH 1 0 pΩq, ∇v 1 tapu, vq ¡ apu h , vqu C gro C ¡3 mon C K u ¡ u h en sup vH 1 0 pΩq, ∇v 1
tpf, vq ¡ pσp∇ s u h q, ∇ s vqu.

and thus

u ¡ u h en ¤ C gro C ¡3 mon C K sup vH 1 0 pΩq, ∇v 1
tpf, vq ¡ pσp∇ s u h q, ∇ s vqu.

(3.28)

Note that, due to the symmetry of σ we can replace ∇ s v by ∇v in the second term inside the supremum. Now x v H 1 0 pΩq, such that ∇v 1. Since σ h Hpdiv, Ωq, we can insert p∇ ¤ σ h , vq pσ h , ∇vq 0 into the term inside the supremum and obtain pf, vq ¡ pσp∇ s u h q, ∇vq pf ∇ ¤ σ h , vq pσ h ¡ σp∇ s u h q, ∇vq.

pf ∇ ¤ σ h , vq § § ¤ § § § Ţ T h pf ∇ ¤ σ h , v ¡ Π 0 T vq T § § § ¤ Ţ T h h T π f ∇ ¤ σ h T ∇v T , (3.30)
whereas the CauchySchwarz inequality applied to the second term directly yields § § pσ h ¡ σp∇ s u h q, ∇vq § § ¤ Ţ T h σ h ¡ σp∇ s u h q T ∇v T .

Inserting these results in (3.28) and again applying the CauchySchwarz inequality yields the result.

Distinguishing the dierent error components

The goal of this section is to elaborate the error estimate (3.27) so as to distinguish dierent error components using the equilibrated stress tensors of Constructions 3.7 and 3.8. This distinction is essential for the development of Algorithm 3.14, where the mesh and the stopping criteria for the iterative solver are chosen adaptively.

Notice that in Theorem 3.10 we don't necessarily need σ h to be the stress tensor obtained in Construction 3.5. We only need it to satisfy two properties: First, equation (3.29) requires σ h to lie in Hpdiv, Ωq. Second, in order to be able to apply the Poincaré inequality in (3.30), σ h has to satisfy the local equilibrium relation pf ¡ ∇ ¤ σ h , vq T 0 dv P 0 pTq dT T h . Theorem 3.12 (A posteriori error estimate distinguishing dierent error sources). Let u be the analytical solution of (3.13), u k h the discrete solution of (3.15), and σ h :

σ k h,disc σ k h,lin . Then, u ¡ u k h en ¤ c 2C gro C ¡3 mon η k disc η k lin η k quad η k osc ¨, (3.32)
where the local discretization, linearization, quadrature and oscillation error estimators on each T T h are dened as

η k disc,T : σ k h,disc ¡ σp∇ s u k h q T , (3.33a)
η k lin,T : σ k h,lin T ,

(3.33b) η k quad,T : σp∇ s u k h q ¡ σp∇ s u k h q T , (3.33c)
η k osc,T :

h T π f ¡ Π p¡1 T f T , (3.33d)
with Π p¡1 T denoting the L 2 -projection onto P p¡1 pTq, and for each error source the global estimator is given by

η k ¤ : ¡ 4 Ţ T h pη k ¤,T q 2 © 1 {2 . ( 3
u¡u k h en ¤ c 2C gro C ¡3 mon £ Ţ T h h T π f ∇ ¤ pσ k h,disc σ k h,lin q T σ k h,lin σ k h,disc ¡ σp∇ s u k h q T ¨2 1 {2
.

Applying the second property of Lemma 3.9 in the rst term yields the oscillation error estimator. In the second term we add and substract σp∇ s u k h q and apply the triangle inequality to

obtain u ¡ u k h en ¤ c 2C gro C ¡3 mon £ Ţ T h η k disc,T η k lin,T η k quad,T η k osc,T ¨2 1 {2
.

Owing to (3.34), the previous bound yields the conclusion.

Remark 3.13 (Quadrature error). In practice, the projection σp∇ s u k h q of σp∇ s u k h q onto P p¡1 pT h q for a general nonlinear stress-strain relation cannot be computed exactly. The quadrature error estimator η k quad,T measures the quality of this projection.

Adaptive algorithm

Based on the error estimate of Theorem 3.12, we propose an adaptive algorithm where the mesh size is locally adapted, and a dynamic stopping criterion is used for the linearization iterations. The idea is to compare the estimators for the dierent error sources with each other in order to concentrate the computational eort on reducing the dominant one. For this purpose, let γ lin , γ quad p0, 1q, be user-given weights and Γ ¡ 0 a chosen threshold that the error should not exceed.

Algorithm 3.14 (Adaptive algorithm).

Mesh adaptation loop

1. Choose an initial function u 0 h H 1 0 pΩq P p pT h q and set k : 1 2. Set the initial quadrature precision ν : 2p (exactness for polynomials up to degree ν)

η k quad,T ¤ γ quad pη k disc,T η k lin,T η k osc,T q dT T h , (3.37a) η k lin,T ¤ γ lin pη k disc,T η k osc,T q dT T h , (3.37b)
where it is also possible to dene local weights γ lin,T and γ quad,T for each element. These local stopping criteria are necessary to establish the local eciency of the error estimators in the following section, whereas the global criteria are only sucient to prove global eciency.

Local and global eciency

Let us start by introducing some additional notation used in this section. For a given element T T h , the set T T collects the elements sharing at least a vertex with T . The set F h contains all faces (if d 2 we will, for simplicity, refer to the edges as faces) of the mesh and is decomposed into boundary faces F ext h and interfaces F int h . For a vertex a V h , we denote by F a the faces containing a, and by F ωa the faces of the patch ω a . For any T T h the set F T contains the faces of T , whereas F T T collects all faces sharing at least a vertex with T and we denote F int

T T F T T F int h .
In what follows we let a À b stand for a ¤ Cb with a generic constant C, which is independent of the mesh size, the domain Ω and the stress-strain relation, but that can depend on the shape regularity of the mesh family tT h u h and on the polynomial degree p.

To prove eciency, we will use a posteriori error estimators of residual type. Following [104,105] we dene for X Ω the functional R X : H 1 pXq Ñ H ¡1 pXq such that, for all v H 1 pXq, w H 1 0 pXq, xR X pvq, wy X : pσp∇ s vq, ∇ s wq X ¡ pf, wq X .

Ţ I T T h 2 T I ∇ ¤ σp∇ s u k h q Π p T f 2 T I F F int T T h F σp∇ s u k h qn F 2 F , pη k S,T q 2 : Ţ I T T h 2 T I ∇ ¤ pσp∇ s u k h q ¡ σp∇ s u k h qq 2 T I F F int T T h F pσp∇ s u k h q ¡ σp∇ s u k h qqn F 2 F .
(3.38)

The quantity η k S,T obviously measures the quality of the approximation of σp∇ s u k h q by σp∇ s u k h q and can be estimated explicitly. The following result is shown in [104, Section 3.3]. Denoting by η k osc,T T :

2 2 °T I T T pη k osc,T I q 2 @ 1 {2 , it holds η k U,T À R T T pu k h q H ¡1 pT T q η k S,T η k osc,T T . (3.39)
In order to bound the dual norm of the residual, we need an additional assumption on the stress-strain relation which, in particular, implies the growth assumption (3.2a).

Assumption 3.15 (Stress-strain relation II

). There exists a real number C Lip p0, Vq such that, for all τ, η R d¢d and to the fact that ¡∇ ¤ σp∇ s uq f L 2 pT T q, using the CauchySchwarz inequality and the Lipschitz continuity (3.40) of σ, it is inferred that

R T T pu k h q H ¡1 pT T q : sup wH 1 0 pT T q, w H 1 0 pT T q ¤1 pσp∇ s u k h q, ∇ s wq T T ¡ pf, wq T T sup wH 1 0 pT T q, w H 1 0 pT T q ¤1 pσp∇ s u k h q ¡ σp∇ s uq, ∇ s wq T T ¤ sup wH 1 0 pT T q, w H 1 0 pT T q ¤1 σp∇ s u k h q ¡ σp∇ s uq T T ∇ s w T T ¤ C Lip ∇ s pu ¡ u k h q T T .
Thus, by (3.39), the previous bound, and the strong monotonicity (3.2b) it holds

η k U,T À C Lip C ¡1 mon u ¡ u k h en,T T η k S,T η k osc,T T .
(3.41) Theorem 3.16 (Local eciency). Let u H 1 0 pΩq be the solution of (3.13), u k h H 1 0 pΩq P p pT h q be arbitrary and σ k h,disc and σ k h,lin dened by Constructions 3.7 and 3.8. Let the local stopping criteria (3.37) be veried. Then it holds for all T T h ,

η k disc,T η k lin,T η k quad,T η k osc,T À C Lip C ¡1 mon u ¡ u k h en,T T η k S,T η k osc,T T . (3.42)
It is well known that there exist nonconforming nite element methods which are equivalent to mixed nite element methods using the BrezziDouglasMarini spaces (see e.g. [6]). Following the ideas of [50, 55, 62] and references therein, we use these spaces to prove Theorem 3.16. We will denote by M h pω a q the extension to vector valued functions of the nonconforming space introduced in [6] on a patch ω a . Recall that Σ h pω a q is the subspace of Hpdiv, Ωq containing tensor-valued piecewise polynomials of degree at most p.

For d 2, the space M T on a triangle T T h is given by M T :

5 tv P p 2 pTq | v |F rP p 1 pFq dF F T u if p is even, tv P p 2 pTq | v |F P p pFq Pp 2 pFq dF F T u if p is odd, (3.43)
where Pp 2 pFq is the L 2 pFqorthogonal complement of P p 2 pFq in P p 1 pFq. The degrees of freedom are given by the moments up to degree pp ¡ 1q inside each T T h and up to degree p on each edge F F h . On a patch ω a this means that M h pω a q contains vector-valued functions lying piecewise in M T such that p m h , v h q F 0 dF F a zF ext h dv h P p pFq. (3.44) We will denote by M a h the subspace of M h pω a q with functions m h verifying pm h , zq ωa 0 dz RM d ,

(3.45) if a V int h , and pm h , v h q F 0 dF F a F ext h dv h P p pFq, (3.46) if a V ext h .
We will use the space M a h together with Proposition 3.17 to prove Theorem 3.16. For Proposi- tion 3.17 we introduce two equivalent formulations of Construction 3.7 based on the following spaces

ΣT : tτ Σ T | pτ, µq T 0 dµ Λ T u, (3.47) Σh pω a q : tτ h L 2 pΩq | τ h ΣT dT T a u, (3.48) Σa h : Σ a h Σh pω a q tτ h Σ a h | pτ h , µ h q ωa 0 dµ h Λ a h u, (3.49) L a h : tl h P p pF ωa q | l h 0 on fω a if a V int h , l h 0 on fω a zfΩ if a V ext h u, (3.50)
where F ωa collects the faces of the patch. The rst equivalent formulation of Construction 3.7

A POSTERIORI ERROR ESTIMATE AND ADAPTIVE ALGORITHM

consists in nding σ a h Σa h and r a h V a h such that for all pτ h , v h q Σa h ¢ V a h pσ a h , τ h q ωa pr a h , ∇ ¤ τ h q ωa pψ a σp∇ s u k h q, τ h q ωa ,

(3.51a) p∇ ¤ σ a h , v h q ωa p¡ψ a f σp∇ s u k h q∇ψ a ¡ y k disc , v h q ωa . (3.51b)
The second formulation is the rst step when hybridizing the mixed problem (3.51). Following

[6] it consists in using the broken space Σh pω a q instead of Σa h and imposing the continuity of the normal stress components by Lagrange multipliers. Its solution is pσ a h , r a h , l a h q Σh pω a q ¢

V a h ¢ L a h such that for all pτ h , v h , l h q Σh pω a q ¢ V a h ¢ L a h pσ a h , τ h q ωa Ţ Ta pr a h , ∇ ¤ τ h q T ¡ F Fω a pl a h , τ h n F F q F pψ a σp∇ s u k h q, τ h q ωa , (3.52a) Ţ Ta p∇ ¤ σ a h , v h q T p¡ψ a f σp∇ s u k h q∇ψ a ¡ y k disc , v h q ωa , (3.52b) ¡ F Fω a p σ a h n F F , l h q F 0, (3.52c)
where we denote by n T F the outward normal vector of T on F and by n F the normal vector of F with an arbitrary, but xed direction. In particular, (3.52a) can be reformulated as

pσ a h ¡ ψ a σp∇ s u k h q, τ T q T pr a h , ∇ ¤ τ T q T F F T pl a h , τ T n T F q F dτ T ΣT dT T a . (3.53)
Proposition 3.17. Let a V h and let pσ a h , r a h , l a h q Σh pω a q ¢ V a h ¢ L a h be dened by (3.52).

Let ra h be a vector-valued function verifying for all T T a and for all F F T ,

ra h|T M T , (3.54a) Π L F ra h|F l a h|F , (3.54b) Π V T ra h|T r a h|T , (3.54c)
where Π L F and Π V T denote, respectively, the L 2 -projections on L F P p pFq and V T P p¡1 pTq.

Then ra h M a h .

Proof. From dimpV T q 3dimpL F q p 2 p 3p2p 2q p 2 7p 6 dimpM T q we infer that problem (

η k disc,T σ k h,disc ¡ σp∇ s u k h q T À η k U,T η k osc,T T . (3.56)
We dene ra h by (3.54). Then using the fact that ra h M a h by Proposition 3.17 and [109, Lemma 5.4], stating that the dual norm on M h is an upper bound for the H 1 -seminorm, we obtain

σ a h ¡ ψ a σp∇ s u k h q ωa ¤ ∇r a h ωa À sup m h M a h , ∇m h 1 pσ a h ¡ ψ a σp∇ s u k h q, ∇m h q ωa . (3.57) Now x m h M a h such that ∇m h ωa 1. Then, by (3.44), it follows pσ a h ¡ ψ a σp∇ s u k h q, ∇m h q ωa Ţ Ta pσ a h ¡ ψ a σp∇ s u k h q, ∇m h q T ¡ Ţ Ta p∇ ¤ σ a h ¡ ∇ ¤ pψ a σp∇ s u k h q, m h q T loooooooooooooooooooooooomoooooooooooooooooooooooon :T 1 F Fa p ψ a σp∇ s u k h qn F , m h q F looooooooooooooooomooooooooooooooooon :T 2 .
Using (3.51b) (which, as in the proof of Lemma 3.9, is valid for all v h V h pω a q) and the fact that for all T T a and τ Σ T it holds p∇ ¤ τ, m h q T p∇ ¤ τ, Π V T m h q T , due to the property ∇ ¤ Σ T V T , we can write for the rst term

T 1 ¡ Ţ Ta p¡ψ a f σp∇ s u k h q∇ψ a ¡ ∇ ¤ pψ a σp∇ s u k h qq, Π V T m h q T ¡ Ţ Ta pψ a pf ∇ ¤ σp∇ s u k h qq, Π V T m h q T ¡ Ţ Ta pΠ p T f ∇ ¤ σp∇ s u k h q, ψ a Π V T m h q T ¤ £ Ţ Ta h 2 T ψ a pΠ p T f ∇ ¤ pσp∇ s u k h qqq 2 T 1 {2 £ Ţ Ta h ¡2 T m h 2 T 1 {2 À £ Ţ Ta h 2 T Π p T f ∇ ¤ σp∇ s u k h q 2 T ψ a 2 L V pTq 1 {2 ∇m h ωa ,
where we used the Cauchy-Schwarz, the discrete Poincaré inequality of [108, Theorem 8.1] 3.5. NUMERICAL RESULTS together with (3.45) if a V int h and the discrete Friedrichs inequality of [108, Theorem 5.4] together with (3.46) if a V ext h , and ψ a L V pTq 1. For the second term we proceed in a similar way, using the discrete trace inequality m h F À h ¡1 {2 F m h T , and obtain

T 2 F F int a pψ a σp∇ s u k h qn F , m h q F ¤ ¤ ¥ F F int a h F ψ a σp∇ s u k h qn F 2 F 1 {2 ¤ ¥ F F int a h ¡1 F m h 2 F 1 {2 À ¤ ¥ F F int a h F σp∇ s u k h qn F 2 F 1 {2 ∇m h ωa .
Inserting these results into (3.57) yields (3.56).

From the local stopping criteria (3.37), the denition of the discretization error estimator (3.33a) and the local approximation property (3.56) it follows that

η k disc,T η k lin,T η k quad,T À η k disc,T σ k h,disc ¡ σp∇ s u k h q T À η k U,T η k osc,T T .
Then (3.41) yields the result.

Theorem 3.18 (Global eciency). Let u H 1 0 pΩq be the solution of (3.13), u k h H 1 0 pΩq P p pT h q be arbitrary and σ k h,disc and σ k h,lin dened by Constructions 3.7 and 3.8. Let the stopping criteria (3.35) and (3.36) be veried. Then it holds

η k disc η k lin η k quad η k osc À C Lip C ¡1 mon u ¡ u k h en η k S η k osc . ( 3
η k disc η k lin η k quad η k osc À η k disc η k osc À ¡ 2 Ţ T h pη k U,T T η k osc,T T q 2 © 1 {2 η k osc À η k U η k osc .
(3.59)

Then, using again (3.41) we obtain the result.

Numerical results

In this section we illustrate numerically our results on two test cases, both performed with the Code_Aster 1 software, which uses conforming nite elements of degree p 2. Our intention is, rst, to show the relevance of the discretization error estimators used as mesh renement pigure QFP ! vEshped domin with liner elstiity modelF histriution of the error estimtors @topA nd the nlytil error @ottomA for the initil mesh @leftA nd fter three @middleA nd six @rightA dptive mesh re(nementsF

L-shaped domain

Following [4,66,82], we consider the L-shaped domain Ω p¡1, 1q 2 zpr0, 1s ¢ r¡1, 0sq, where for the linear elasticity case an analytical solution is given by upr, θq

1 2µ r α £ cospαθq ¡ cosppα ¡ 2qθq A sinpαθq sinppα ¡ 2qθq
, with the parameters µ 1.0, λ 5.0, α 0.6, A 31{9.

Linear elasticity model

We compute the analytical error and its estimate on two series of unstructered meshes. Starting with the same initial mesh, we use uniform mesh renement for the rst one and adaptive renement based on the error estimate for the second series. estimates and errors for each mesh, as well as their eectivity index corresponding to the ratio of the estimate to the error. We obtain eectivity indices close to one, showing that the estimated error value lies close to the actual one, what we can also observe in the graphics on the left. As expected, the adaptively rened mesh series has a higher convergence rate, with corresponding error an order of magnitude lower for 10 3 elements.

HenckyMises model

For the HenckyMises model we choose the Lamé functions μpρq : a bp1 ρ 2 q ¡1 {2 , λpρq : κ ¡ The results are shown in Figure 3.4. In the left graphic we observe that the linearization error estimate in the adaptive case is much higher than in the one without adaptive stopping criterion. We see that this does not aect the discretization error estimator. The table shows the number of performed Newton iterations for both cases. The algorithm using the adaptive stopping criterion is more ecient. In the right graphic we compare the global distretization error estimate on two series of meshes, one rened uniformly and the other one adaptively, based on the local discretization error estimators. Again the convergence rate is higher for the adaptively rened mesh series.

Notched specimen plate

In our second test we use the two nonlinear models of Examples 3.3 and 3.4 on a more application-oriented test. The idea is to set a special sample geometry yielding to a model discrimination test, namely dierent physical results for dierent models. We simulate the uniform traction of a notched specimen under plain strain assumption (cf. Figure 3.5).

The notch is meant to favor strain localization phenomenon. We consider a domain Ω p0, 10mq ¢ p¡10m, 10mqztx R 2 | xm ¡ p0, 11mq T ¤ 2mu, we take f 0, and we prescribe a displacement on the boundary leading to the following Dirichlet conditions:

u x 0m if x 0m, u y ¡1.1 ¤ 10 ¡3 m if y ¡10m, u y 1.1 ¤ 10 ¡3 m if y 10m.
In many applications, the information about the material properties are obtained in uniaxial corresponding to the Lamé parameters µ 3 26 ¤ 10 9 Pa and λ 9 52 ¤ 10 9 Pa. For both stressstrain relations we apply Algorithm 3.14 with γ lin 0.1. We rst compare the results to a computation on a very ne mesh to evaluate the remeshing based on the discretization error estimators. Secondly, we perform adaptive remeshing based on these estimators but without applying the adaptive stopping of the Newton iterations and compare the two series of meshes.

As in Section 3.5.1, we verify if the reduced number of iterations impacts the discretization error.

Figure 3.6 shows the result of the rst part of the test. In each of the four images the left specimen corresponds to the HenckyMises and the right to the isotropic damage model. To illustrate the dierence of the two models, the top left picture shows the trace of the strain tensor. This scalar value is a good indicator for both models, representing locally the relative volume increase which could correspond to either a damage or shear band localization zone.

In the top right picture we see the distribution of the discretization error estimators in the reference computation (209,375 elements), whereas the distribution of the estimators on the sixth adaptively rened mesh is shown in the bottom right picture (60,618 elements for Hencky Mises, 55,718 elements for the damage model). The corresponding meshes and the initial mesh for the adaptive algorithm are displayed in the bottom left of the gure. To ensure a good discretization of the notch after repeated mesh renement, the initial mesh cannot be too corse in this curved area. We observe that the adaptively rened meshes match the distribution of the discretization error estimators on the uniform mesh, and that the estimators are more pigure QFT ! xothed speimen plteD omprison etween renky!wises @left in eh pitureA nd dmge model @rightAF Top left: trp∇s u h qF Top right: η disc on (ne mesh @no dptive re(nementAF Bottom left: meshes fter six dptive re(nementsF Bottom middle: initil meshF Bottom right: η disc on the dptively re(ned meshesF evenly distributed on these meshes.

The results of the second part of the test are illustrated in Figure 3.7. As for the L-shape test, we observe that the reduced number of Newton iterations does not aect the discretization error estimate, nor the overall error estimate which is dominated by the discretization error estimate if the Newton algorithm is stopped.

Conclusions

In this work we have developed an a posteriori error estimate for a wide class of hyperelastic problems. The estimate is based on stress tensor reconstructions and thus independent of the stress-strain relation, except for two constants. In a nite element software providing dierent mechanical behavior laws it can be directly applied to any of these laws. The assumptions we make on the stress-strain relation are only used to obtain the equivalence of the energy norm and the dual norm of the residual of the weak formulation. Using the latter as error The main goal of engineering and research in solid mechanics is to predict this deformation, which depends on several factors, such as the geometry and the physical properties of the material, and of course the way the forces are applied. The rst step for this prediction is the denition of variables describing the intrinsic rigidity of the material and the observation of their reaction during tests simulating dierent loadings on material specimen of simple geometries. Let Ω R 3 be the considered body and p0, t F q with t F ¡ 0 the time interval. By dening in each point of the structure a local mechanical state, it is possible to distinguish the contribution of the geometry to the deformation from the one due to material properties.

Under the hypothesis of small deformations, this local state is typically characterized by the stress σ : Ω ¢ p0, t F q Ñ R 3¢3 sym and the strain ε : Ω ¢ p0, t F q Ñ R 3¢3 sym . The evaluation of the test results then yields a model describing the material properties as the relation between these quantities or their time derivatives 9

σ and 9 ε. The second step is to apply this model to a specic case of interest, involving in general more complicated geometries and loading terms than in the tests. In general, this can not be done analytically, so the challenge of the second step is the numerical resolution of the resulting problem, taking into account the material behaviour described by the model. Nowadays, some well-established models for groups of materials exist, involving material-specic parameters.

Step one then often consists in identifying the parameters corresponding to the considered material through a series of experiments.

In this chapter we start by reviewing a general elasto-plastic stress-strain relation in Section 4.1, which we then, in Section 4.2, consider in the context of poro-mechanical problems.

Section 4.3 is dedicated to the numerical solution of the obtained problem. In Section 4.4 we discuss an error measure and its estimation using the equilibrated ux reconstruction techniques elaborated in the previous chapters. In Section 4.5 we present an adaptive algorithm based on the error estimate. We then discuss some mechanical models used in the modeling of soils and rocks in Section 4.6. Finally, in Section 4.7 we assess the error estimate on an analytical test case and apply the adaptive algorithm to the simulation of a tunnel excavation in three space dimensions. Let us recall that elasticity is based on the assumption that every deformation is reversible, meaning that a body always returns to its original size and shape when loadings are removed.

A typical example for elastic behavior is a spring. On the other hand, we speak of plastic behaviour, whenever the deformation in response to an applied loading is non reversible and the shape of the material stays changed after the loading is removed, as for example with a piece of modeling clay. Some materials only behave elastically until the stress reaches a certain threshold and show plastic behaviour once the threshold is exceeded -like a spring that always returns in its originial position, unless it is pulled too hard. These materials are called elasto-plastic and the threshold expressing the transition between elastic and plastic behaviour is called yield and plays a fundamental role in models describing these materials.

An important property for the modelization of elasto-plastic materials is that the strain tensor ε is composed of an elastic and a plastic part, i.e.

ε ε e ε p ,

The yield function

To modelize the material behavior, we need a criterion telling us if the reaction at a point x in the material to a stress σ will be elastic or plastic deformation. This criterion is expressed as a convex function F : R 3¢3 sym Ñ R, called the yield function, dening the yield surface tσ; F pσq 0u in the space spanned by the components of σ. The reaction of a material to the stress σ is dened by the position of σ with regard to the yield surface of the material :

• if σ lies inside the yield surface, i.e. F pσq 0: the material is elastic at σ. In particular this means that ε p 0 and thus that stress and strain satisfy Hooke's law.

• if the point σ lies on the yield surface, i.e. F pσq 0, and the direction of 9 σ points inwards the yield surface: the material reaction is still elastic. Again, stress and strain satisfy Hooke's law.

• if the point σ lies on the yield surface and the direction of 9 σ points outwards the yield surface: the material is plastic at σ.

Hardening and softening

It is quite intuitive to understand that once a body has been irreversibly deformed, its behavior might change. This is the case if the energy absorbed by the material during the deformation process modies the internal structure of the material. The resulting eect is called hardening or softening, depending on whether the deformed material will support more or less stress than the original one before reaching the yield criterion. This means that the value of F for a given σ increases or decreases with ongoing plastic deformation and consequently that the convex elasticity zone tσ; F pσq 0u decreases or increases respectively. To account for this change of the yield surface, we introduce the parameter α : Ω ¢ p0, t F q Ñ R in the function F . This parameter, called accumulated plastic strain, describes the length of the deformation path until the time t, and thus is a way of tracking the strain history in the time interval p0, t F q.

In general, geomaterials, which are the topic of this chapter, show pure softening behavior, meaning that after every plastic deformation the resulting yield surface tpσ, αq; F pσ, αq 0u lies entirely inside the previous one. Equivalently, we speek of pure hardening behavior when the yield surface lies entirely outside the previous surfaces after any plastic deformation. At this point, it is already worth noticing that only elasto-plasticity problems describing purely hardening materials have a unique solution. The reason can be explained in a simplied way using Figure 4.1, showing two basic stress-strain relations for a one-dimensional problem; left for a hardening material, and right for a softening one. Both materials show an elastic behavior until σ rst reaches σ crit and plastic behavior if ε increases from that point on. At every point of the domain, the couple pε, σq lies on the curve indicated in Figure 4.1. Let us assume that the material satises the equilibrium ∇ ¤ σ 0, which -in one dimension -means that σ is constant over the whole domain, and that the domain is under traction with an increasing displacement u D . Then, since σ is constant, we see that in the case of the left picture, ε will be constant over the domain too. In the right picture on the other hand, once σ reaches the value σ crit , each point can either take the left or the right branch of the curve, as long as the integral over ε is equal to u I . The middle picture shows what happens when the traction then decreases to zero: only the elastic deformation ε e is reversible and obeys Hooke's law σ Eε e , where E ¡ 0 is Young's modulus, whereas the plastic deformation remains. When -as illustrated in the right picture -the traction increases again, the material resists to less stress than the undeformed one.

The elasto-plastic behavior law

We recall that the mechanical behavior describing the relation between stress and strain is essential for the resolution of a mechanical problem (where external forces are known and the resulting displacement is sought), as the stress satises the equilibrium with the external forces, and the strain is the symmetric gradient of the displacement. To derive this stressstrain relation taking into account the yield criterion, we assume that the yield function F : R 3¢3 sym ¢ R Ñ R satises the following properties:

• F is a convex and piecewise analytical function of pσ, αq • The point p0, 0q does not lie on the yield surface, i.e. F p0, 0q 0 • F is dierentiable in all points pσ, αq satisfying F pσ, αq 0

In plasticity, the behavior law is usually expressed as the relation between the time derivatives of the stress and the strain. For a given time t p0, t F q, we distinguish between the elastic behavior, which takes place if pσ, αq is such that F pσ, αq 0 or if F pσ, αq 0 and pfF{fσq : 9 σ ¤ 0 (i.e. pσ, αq lies inside or will move towards the inside of the yield surface), and the plastic behavior occuring if F pσ, αq 0 and pfF{fσq : 9 σ ¡ 0. In the elastic case we know that σ and ε satisfy Hooke's law, and that the plastic increment is zero:

9 σ D : 9 ε, (4.2a) 9 p 0, (4.2b)
where D is a symmetric, positive denite fourth-order tensor, and a : b produces a second order tensor by pa : bq ij : °3 k,l1 a ijkl b kl . Similarly we dene pb : a q kl : °3 i,j1 b ij a ijkl .

As mentioned above, in the plastic case we split the strain into two parts. On the one hand

ε e 1 ε σ σ crit,0 ε 1 E ε σ σ crit,0 ε p 1 ε e 1 E E ε e 1 ε σ σ crit,0 σ crit,1 ε p 1 E
pigure RFP ! lsti deformtion with softening 4.1. ELASTO-PLASTICITY the elastic strain ε e still satisfying Hooke's law (4.2a), and on the other hand the plastic strain ε p ε ¡ ε e . In order to write the mechanical behavior law, we eliminate this plastic strain, using the normality assumption, also called Prandtl-Reuss ow rule. It states that, for every t p0, t F q, there exists λ ¥ 0, such that 9 ε p λ fF fσ pσ, αq and 9 α ¡λ fF fα pσ, αq, (

meaning that, on the one hand, the direction of the plastic ow is equal to the outward normal to the yield surface in the stress space, and that, on the other hand, the direction of pε p , αq in the σ ¡ α¡ space is tangential to the yield surface. We also know that the point pσ, pq moves along the yield surface, i.e. fF fσ pσ, αq : 9 σ fF fα pσ, αq 9 α 0.

(4.4)

In the following we denote f σ F : fF fσ pσ, αq and f α F : fF fα pσ, αq. Since 9 ε e satises Hooke's law, we can write

0 D : p9 ε ¡ 9 ε p q ¡ 9 σ. (4.5)
Multiplying this equation by f σ F , and using (4.3) and (4.4) we obtain

0 f σ F : D : 9 ε ¡ f σ F : D : 9 ε p ¡ f σ F : 9 σ f σ F : D : 9 ε ¡ pf σ F : D : f σ F pf α F q 2 qλ, leading to λ f σ F : D : 9 ε f σ F : D : f σ F pf α F q 2 .
Plugging this expression for λ into (4.3) and the so obtained ε p into (4.5) we get the stress-

strain relation 9 σ £ D ¡ D : f σ F f σ F : D f σ F : D : f σ F pf α F q 2 9 ε,
where for two second-order tensors a, b the tensor product yielding a fourth-order tensor is dened as pa bq ijkl : a ij b kl . The variation of the accumulated plastic strain as a function of 9 σ is directly obtained from (4.4).

Summing up, we can state the general elasto-plastic behavior law as follows:

• if pσ, αq are such that F pσ, αq 0 or F pσ, αq 0 and f σ F : • if pσ, αq are such that F pσ, αq 0 and f σ F :

9 σ ¤ 0, then 9 σ D : 9 ε, (4.6a)
9 σ ¡ 0, then 9 σ £ D ¡ D : f σ F f σ F : D f σ F : D : f σ F pf α F q 2 9 ε, (4.6c) 9 α ¡pf α F q ¡1 f σ F : 9 σ. (4.6d)

Poro-mechanical coupling

In this section we integrate the stress-strain relation (4.6) into the mechanical constitutive relation of the coupled hydro-mechanical problem introduced in Section 1.2. We express the eective stress σ I , which together with the stress σ p induced by the pressure yields the total stress σ. We recall that f : Ω ¢ p0, t F q Ñ R 3 and g : Ω ¢ p0, t F q Ñ R denote respectively the volumetric body force and a uid source, b ¡ 0 is the Biot-Willis coecient, c 0 the specic storage coecient, and κ : Ω Ñ rκ S , κ U s with 0 κ S κ U V the hydraulic conductivity. The goal is to determine the displacement eld u : Ω ¢ p0, t F q Ñ R 3 and the pressure p : Ω ¢ p0, t F q Ñ R verifying the equilibrium equations ¡∇ ¤ σ I pεpuq,

αq bpI ¨ f in Ω ¢ p0, t F q, (4.7a)
f t pc 0 p b∇ ¤ uq ∇ ¤ φppq g in Ω ¢ p0, t F q, ( 4.7b)
where the Darcy velocity is dened as φppq ¡κ∇p and σ I and α are derived from (4.6), with σ I instead of σ, taking into account the initial conditions up¤, 0q u 0 , αp¤, 0q 0 and pp¤, 0q p 0 in Ω,

with functions u 0 : Ω Ñ R 3 and p 0 : Ω Ñ R. For simplicity, we consider homogenous Dirichlet boundary conditions u 0 and p 0 on fΩ ¢ p0, t F q.

(4.7d)

As in Chapter 2, we next write the weak formulation. We dene the following sets and spaces

E tpσ, αq R 3¢3 sym ¢ R; F pσ, αq 0 or F pσ, αq 0 and f σ F : 9 σ ¤ 0u, P tpσ, αq R 3¢3 sym ¢ R; F pσ, αq 0 and f σ F : 9 σ ¡ 0u U H 1 p0, t F ; H 1 0 pΩqq, P H 1 p0, t F ; H 1 0 pΩqq.
where

D I D : f σ F f σ F : D f σ F : D : f σ F pf α F q 2 .
For hardening materials, existence and uniqueness of a solution for linear poroplasticity (i.e.

Numerical solution

In this section we explain how (4.9) is numerically solved in Code_Aster. We assume that the problem is well posed, using for example the regularization techniques of [58], see also [45,46]. We start by introducing some notation before we identify the space and time discretization, which is the same as in Chapter 2 in the two-dimensional case: an H 1 -conforming space discretization using the TaylorHood nite element spaces and an implicit Euler time stepping.

Notation

For the discretization of the time interval p0, t F q, we consider a sequence of discrete times pt n q 0¤n¤N with t 0 0 and t N t F , and such that t i t j if i j. For each 1 ¤ n ¤ N we dene the time step τ n : t n ¡ t n¡1 , the time interval I n : pt n¡1 , t n q, and the discrete backward dierencing operator f n t ϕ : τ ¡1 n pϕ n ¡ ϕ n¡1 q, where for a space-time function ϕ we denote ϕ n : ϕp¤, t n q.

At each time step n, let T n h be a conforming triangulation of Ω (i.e. a set of closed tetrahedra with union Ω and such that, for any distinct T 1 , T 2 T n h , their intersection T 1 T 2 is either a common face, a common edge, a vertex or the empty set) verifying the minimum angle condition (i.e. there exists α min ¡ 0 such that the minimum angle α T of a tetrahedron T T n h

satises α T ¥ α min ). V n h V n,int h V n,ext
h denotes the set of vertices, divided into interior vertices and exterior vertices lying on the boundary fΩ. For any subdomain ω Ω the set V n ω contains the vertices in ω, and for all vertices a V n h the set T a collects the elements sharing the vertex a. The corresponding open subset of Ω is ω a and is called a patch. For any element T T n h , h T denotes its diameter and n T its unit outward normal vector.

Discrete formulation

We next introduce the used discretization spaces. Let P k pTq be the space of trivariate polynomials in T T n h of total degree at most k ¥ 1, and let P k pT n h q tϕ L 2 pΩq; ϕ |T P k pTq dT

T n h u denote the corresponding broken space over T n h . The TaylorHood nite element spaces of degree k ¥ 1 are dened as U n h : P k 1 pT n h q H 1 0 pΩq and P n h : P k pT n h q H 1 0 pΩq. We indicate discrete space-time functions by the subscript hτ , whereas the discrete space function obtained when xing one time step t n is denoted as before by ϕ n h : ϕ hτ p¤, t n q.

Assumption 4.1. For simplicity, we assume that the source functions f and g are piecewise constant in time and in space, i.e. that the following holds:

• f and g are constant in each time interval. We denote f n : f |In and g n : g |In .

• For each 1 ¤ n ¤ N it holds f n P 0 pT n h q and g n P 0 pT n h q. We denote by p¤, ¤q h the discrete inner L 2 -product calculated using Gauss integration. We omit the subscript h, whenever the integral can be calculated exactly. The discrete formulation of the problem then reads: xing u 0 u 0 , p 0 p 0 and α 0 0, nd for each 1 ¤ n ¤ N with given u n¡1 h , p n¡1 h and α n¡1 the functions pu n h , p n h q U n h ¢ P n h and α n , such that

pσ I pεpu n h q, α n , εpv h qq h bpp n h , ∇ ¤ v h q pf n , v h q dv h U n h , (4.11a)
pρ ¡1 mpf n t u hτ , f n t p hτ q, q h q ¡ pφpp n h q, ∇q h q pg n , q h q dq h P n h ,

with the eective stress tensor and accumulated plastic strain increment at each Gauss point given by f n t σ I pεpu hτ q, α n q 5 D : f n t εpu hτ q if pσ I pεpu n h q, α n q, α n q E, pD ¡ D I q : f n t εpu hτ q if pσ I pεpu n h q, α n q, α n q P, (4.11c)

f n t α n 5 0 if pσ I pεpu n h q, α n q E, ¡pf α F q ¡1 f σ F : f n t σ I pεpu hτ q, α n q if pσ I pεpu n h q, α n q P, (4.11d)
where we denote by u hτ and p hτ the discrete space-time functions such that, at each time step 0 ¤ n ¤ N , it holds pu hτ , p hτ qp¤, t n q pu n h , p n h q and f t u n hτ :

f t u hτ |In τ ¡1 n pu n h ¡ u n¡1 h q and f t p n hτ : f t p hτ |In τ ¡1 n pp n h ¡ p n¡1 h q if n ¥ 1.

Linearization

In order to solve the nonlinear problem (4.11) numerically, we apply a linearization algorithm, leading at each 1 ¤ n ¤ N to a series of linear algebraic systems, whose solutions converge to the solution of (4.11). In Code_Aster, this is done using the NewtonRaphson method. It is applied only to the nonlinear term, expressing the dependency of the eective stress tensor on the primal variable u through the strain tensor εpuq. Note that, since ε depends linearly on u, it holds fσ I fu puq¤v fσ I fε pεpuqq : εpvq for any v U. Let pu n,0 h , p n,0 h q U n h ¢P n h be an initial guess and α n,0 α n¡1 . Then the linearized discrete problem at time step n and iteration i ¥

1 is to nd, for given u n¡1 h , p n¡1 h , α n¡1 , u n,i¡1 h , p n,i¡1 h and α n,i¡1 , the functions pu n,i h , p n,i h q U n h ¢ P n h such that, for all pv h , q h q U n h ¢ P n h , ¡ fσ I fε pεpu n,i¡1 h qq : εpu n,i h q, εpv h q © ¡ bpp n,i h , ∇ ¤ v h q £ fσ I fε pεpu n,i¡1 h qq : εpu n,i¡1 h q ¡ σ I pεpu n,i¡1 h q, α n,i¡1 q, εpv h q pf n , v h q, (4.12a)
pc 0 pf n t p n,i hτ , qq bp∇ ¤ f n t u n,i hτ , q h q ¡ pφpp n,i h q, ∇q h q pg n , q h q. (

This problem corresponds to a linear equation system A n,i¡1 U n,i G n,i¡1 . We stop the iterations, when

A n,i¡1 U n,i ¡ G n,i¡1 V ¤ Γ res G n,i¡1 V , (4.13)
where Γ res is typically in the order of 10 ¡6 . Let j denote the corresponding iteration of the Newton solver. We then set pu n h , p n h q pu n,j h , p n,j h q and move forward to the next time step.

Initial guess

It is well-known that the convergence of the Newton method depends on the choice of the initial guess pu n,0 h , p n,0 h q, especially if the yield function F is such that the convergence is not guaranteed. In any case, a good initial guess reduces in practice the number of iterations necessary to approximate the solution of the nonlinear problem (4.11) to a sucient accuracy.

In Code_Aster the initial guess at each time step n is obtained by linearizing equation (4.11) at pu n¡1 h , p n¡1 h q with respect to the time [44]. Since the only dependency on t is through pu, pq, we obtain, similarly to (4.12), ¡ fσ I fε pεpu n¡1 h qq : εpu n,0 h q, εpv h q © ¡ bpp n,0 h , ∇ ¤ v h q £ fσ I fε pεpu n¡1 h qq : εpu n¡1 h q ¡ σ I pεpu n¡1 h q, α n¡1 q, εpv h q pf n¡1 , v h q, pc 0 pf n t p n,0 hτ , qq bp∇ ¤ f n t u n,0 hτ , q h q ¡ pφpp n,0 h q, ∇q h q pg n¡1 , q h q, for all pv h , q h q U n h ¢ P n h .

Integration of the mechanical behavior law

To solve (4.12), we have to calculate fσ I fε pεpu n,i¡1 h qq and σ I pεpu n,i¡1 h q, α n,i¡1 q before each iteration. We start with σ I pεpu n,i¡1 h q, α n,i¡1 q, using (4.11c,4.11d). Therefore, we rst have to nd out if the displacement u n,i¡1 h causes elastic or plastic deformation in the material, taking into account the deformation history expressed by α n¡1 . This is done by testing if an elastic deformation leads to an admissible pair pσ I , αq, i.e. if F pσ I , αq ¤ 0. According to (4.11c,4.11d), the elastic deformation leads to σ I pεpu n,i¡1 h q, α n,i¡1 q D : εpu n,i¡1 h q and α n,i¡1 α n¡1 .

We then distinguish the following two cases:

(1) If this pair satises F pσ I pεpu n,i¡1 h q, α n,i¡1 q, α n,i¡1 q ¤ 0, then σ I pεpu n,i¡1 h q, α n,i¡1 q is given by the above expression and for any v h U n h it holds fσ I fε pεpu n,i¡1 h qq : εpv h q D : εpv h q.

(2) If F pσ I pεpu n,i¡1 h q, α n,i¡1 q, α n,i¡1 q ¡ 0, the deformation u n,i¡1 h cannot lead to an elastic deformation, we thus have to apply the plastic behavior law and obtain σ I pεpu n,i¡1 h q, α n,i¡1 q σ I pεpu n¡1 h q, α n¡1 q pD ¡ D I q : εpu n,i¡1 h ¡ u n¡1 h q α n,i¡1 α n¡1 ¡ pf α F q ¡1 f σ F : pσ I pεpu n,i¡1 h q, α n,i¡1 q ¡ σ I pεpu n¡1 h q, α n¡1 qq, and, for any v h U n h , it holds fσ I fε pεpu n,i¡1 h qq : εpv h q D : εpv h q ¡ fpD I : εq fε pεpu n,i¡1 h qq : εpv h q.

Equilibrated ux a posteriori error estimate

In this section we apply the results of Chapters 2 and 3 to the poro-plastic problem introduced in the previous section. We start by presenting the equilibrated uxes. The velocity reconstruction is the extension to three space dimensions of the reconstruction in Section 2.3.1. The stress tensor reconstructions are the same as presented in Section 3.3.2 for nonlinear elasticity problems. We then present the error measure, which, as in Chapter 2, is the dual norm of the residual of the weak formulation (4.9), and derive an upper bound on this error, calculated at each time step and each Newton iteration using the equilibrated ux reconstructions. In the estimate, we distinguish between the linearization error for the mechanical part and the space and time discretization errors for both the mechanical and the hydraulic parts. Finally, we explain the hybridization technique used in the implementation of the ux reconstruction. Comparing the properties of the continous and the discrete uxes we observe the same dierences as in the linear case: Recalling that σpu, p, αq σ I pεpuq, αq σ p ppq, where σ p ppq bpI, the uxes from (4.9) satisfy, for a. e. t p0, t F q, σpu, p, αqptq H s pdiv, Ωq, ¡∇ ¤ σpu, p, αqptq fptq, φppqptq Hpdiv, Ωq, ∇ ¤ φppqptq gptq ¡ f t p∇ ¤ u c 0 pqptq, wheras the discrete uxes from (4.12) at time step 1 ¤ n ¤ N and Newton iteration i ¥ 0 are in general such that σpu n,i hτ , p n,i hτ , α n,i q H s pdiv, Ωq, ¡∇ ¤ σpu, p, αq $f n,i , φpp n,i hτ q Hpdiv, Ωq, ∇ ¤ φpp n,i hτ q $g n,i ¡ f n t p∇ ¤ u n,i hτ c 0 p n,i hτ q.

For each vertex a V h , we dene the corresponding hat function ψ a P 1 pT h q as the piecewise linear function taking value one at the vertex a and zero on all other mesh vertices. In this section, we recall briey the two discrete reconstructions restoring these properties, starting with the Darcy velocity. For each vertex a V n h we then dene the spaces on the element patches around a. If a V n,int h we set W a h :tv h Hpdiv, ω a q; v h|T W h pTq dT T n a and v h ¤ n ωa 0 on fω a u, Q a h :tq h L 2 pω a q; q h|T P k¡1 pTq dT T n a and pq h , 1q ωa 0u, and if a V n,ext h , they are dened as W a h :tv h Hpdiv, ω a q; v h|T W h pTq dT T n a and v h ¤ n ωa 0 on fω a zfΩu, Q a h :tq h L 2 pω a q; q h|T P k¡1 pTq dT T n a u.

Equilibrated Darcy velocity reconstruction

q W a h ¢ Q a h such that for all pv h , q h q W a h ¢ Q a h φ a h , v h ¨ωa ¡ pr a h , ∇ ¤ v h q ωa ¡ ψ a φpp n,i h q, v h © ωa , (4.16a) ∇ ¤ φ a h , q h ¨ωa ¡ ψ a g n ¡ ψ a f n t p∇ ¤ u n,i hτ c 0 p n,i hτ q ∇ψ a ¤ φpp n h q, q h © ωa . (4.16b) Then set φ n,i h : °aV n h σ a h .
For interior vertices it holds ψ a g n ¡ ψ a f t p∇ ¤ u hτ c 0 p hτ q n ∇ψ a ¤ φpp n h q, 1 ¨ωa 0, guaranteeing that the pure Neumann problems (4.16) are well-posed for all a V n,int h . The resulting ux φ n h lies in Hpdiv, Ωq and satises the equilibrium ∇ ¤ φ n h g n ¡ f n t p∇ ¤ u hτ c 0 p hτ q.

(4.17)

The proof is the same as in Lemma 2.11, using also Assumption 4.1.

Equilibrated stress tensor reconstructions

For the three-dimensional computations we chose to rely on the stress reconstructions of introduced in Sections A. Σ T : P k 1 pTq, V T : P k pTq, Λ T : tµ P k pTq; µ ¡µ T u.

Like before, we dene for each vertex a V n h the mixed spaces on the patches around a, distinguishing between interior and boundary vertices. In the former case, we set Σ a h : tτ h Hpdiv, ω a q; τ h|T Σ T dT T n h and τ h n ωa 0 on fω a u, (

V a h : tv h L 2 pω a q; v h|T V T dT

V a h :tv h L 2 pω a q; v h|T V T dT T n h u, (4.19e)
Λ a h :tµ h L 2 pω a q; µ h|T Λ T dT T n h u.

a V n h nd pσ a h , r a h , λ a h q Σ a h ¢ V a h ¢ Λ a h such that for all pτ h , v h , µ h q Σ a h ¢ V a h ¢ Λ a h ,
pσ a h , τ h q ωa pr a h , ∇ ¤ τ h q ωa pλ a h , τ h q ωa pψ a σpu n,i h , p n,i h , α n,i q, τ h q ωa , p∇ ¤ σ a h , v h q ωa p¡ψ a f σpu n,i h , p n,i h , α n,i q∇ψ a ¡ y n,i disc , v h q ωa , pσ a h , µ h q ωa 0, where y n,i disc RM is the unique solution of py n,i disc , zq ωa ¡pf, ψ a zq ωa pσpu n,i h , p n,i h , α n,i q, εpψ a zqq ωa dz RM d .

Then set σ n,i h,disc : °aV h σ a h .

To express the linearization error, we observe that we can rewrite equation (4.12a) as pσ i¡1 pu n,i h , p n,i h , α n,i q, εpv h qq pf n , v h q, where the linear function σ i¡1 maps pu n,i h , p n,i h , α n,i q onto σ i¡1 pu n,i h , p n,i h , α n,i q : fσ I fε pεpu n,i¡1 h qq : εpu n,i h ¡ u n,i¡1 h q σ I pεpu n,i¡1 h q, α n,i¡1 q σ p pp n,i h q.

Construction 4.4 (Linearization error stress reconstruction). For each a V h nd pσ a h , r a h , λ a h q

Σ a h ¢ V a h ¢ Λ a h such that for all pτ h , v h , µ h q Σ a h ¢ V a h ¢ Λ a h ,
pσ a h , τ h q ωa pr a h , ∇ ¤ τ h q ωa pλ a h , τ h q ωa pψ a pσ i¡1 pu n,i h , p n,i h , α n,i q ¡ σpu n,i h , p n,i h , α n,i qq, τ h q ωa , p∇ ¤ σ a h , v h q ωa ppσ i¡1 pu n,i h , p n,i h , α n,i q ¡ σpu n,i h , p n,i h , α n,i qq∇ψ a y n,i disc , v h q ωa , pσ a h , µ h q ωa 0. On the space X n we dene the norm pv, qq 2

Xn :

» In pE εpvq q 2 pl ∇q q 2 dt.

Then the error measure e n,i for the error between the weak solution of (4.9) and the discrete solution pu n,i hτ , p n,i hτ , α n,i q of (4.12) at the time step 1 ¤ n ¤ N and iteration i ¥ 0 is dened by the dual norm of the residual e n,i : Rpu n,i hτ , p n,i hτ , α n,i q X I n sup pv,qqXn, pv,qq Xn 1 xRpu n,i hτ , p n,i hτ , α n,i q, pv, T T n h we dene the following local estimators of the space and time discretization error, as well as of the linearization error for the mechanical part,

η n,i sp,mec,T : E ¡1 σ n,i h,disc ¡ σpu n,i h , p n,i h , α n,i q T , (4.27a)
η n,i tm,mec,T ptq : E ¡1 σpu n,i h , p n,i h , α n,i q ¡ σpu n,i hτ , p n,i hτ , α n,i qptq T ,

η n,i lin,mec,T :

E ¡1 σ n,i h,lin T , (4.27c)
and the local spatial and temporal discretization error estimators for the hydraulic part:

η n,i sp,hyd,T :

t l φ n,i h ¡ φpp n,i h q T , (4.27d)
η n,i tm,hyd,T ptq : t l φpp n,i h q ¡ φpp n,i hτ qptq T .

(

For each of these error sources, the global estimator is given by η :

¤ ¥ 3 » In Ţ T n h η n,i ptq ¨2 dt 1 {2 (4.28)
We can then present the following local and global in time error estimates: Theorem 4.5 (Local in time a posteriori error estimate). Let pu, pq Y be the weak solution of (4.9), let Assumption 4.1 hold and let, at a time step 1 ¤ n ¤ N and a Newton iteration i ¥ 0, pu n,i h , p n,i h q U n h ¢ P n h be the discrete solution of the linearized system (4.12). Let e We present the approach for the stress reconstructions. The same procedure (ignoring the space of skew-symmetric matrices used in the stress reconstruction to impose the weak symmetry) can be applied to the velocity reconstruction and is presented for example in [22, Chapter 7.2] We start by replacing the space Σ a h in (4.19) by the broken BrezziDouglasMarini spaces, i.e.

for each a V n h we set Σ a h : tτ h L 2 pω a q; τ h|T Σ T dT T n h u.

The spaces V a h and Λ a h are still dened as in (4.19). The continuity of the normal components of the stress functions is thus no longer imposed by the choice of the discretization space.

Instead, we use Lagrange multipliers living in the space Ξ a h : P k 1 pF a q, where F a denotes the set of all faces in ω a if a V Γ a ψ a σpu n,i h , p n,i h , α n,i q, γ a ¡ψ a f σpu n,i h , p n,i h , α n,i q∇ψ a ¡ y n,i disc , whereas for the linearization error stress reconstruction they read Γ a pψ a pσ i¡1 pu n,i h , p n,i h , α n,i q ¡ σpu n,i h , p n,i h , α n,i qq, γ a pσ i¡1 pu n,i h , p n,i h , α n,i q ¡ σpu n,i h , p n,i h , α n,i qq∇ψ a y n,i disc .

Then, the equivalent problems consist in seeking pσ a h , r a h , λ a h , χ a

h q Σ a h ¢ V a h ¢ Λ a h ¢ Ξ a h such that for all pτ h , v h , µ h , ξ h q Σ a h ¢ V a h ¢ Λ a h ¢ Ξ a h , pσ a h , τ h q ωa Ţ Tω a pr a h , ∇ ¤ τ h q T pλ a h , τ h q ωa ¡ F Fa pχ a h , τ h n F q F pΓ a , τ h q ωa (4.33a) Ţ Tω a
p∇ ¤ σ a h , v h q T pγ a , v h q ωa , (4.33b) pσ a h , µ h q ωa 0, The problem (4.33a) corresponds to a matrix system of the form

¤ ¦ ¦ ¦ ¦ ¥ A B t 1 B t 2 C t B 1 B 2 C ¤ ¦ ¦ ¦ ¦ ¥ Σ R Λ X ¤ ¦ ¦ ¦ ¦ ¥ F G I 0 0 . (4.34)
Choosing basis functions equal to zero on all but exactly one element for the spaces Σ a h , V a h and Λ a h , the matrices A, B 1 and B 2 are block diagonal. Dening B pB 1 , B 2 q T , Y pR, Λq T and G pG I , 0q Y pBA ¡1 B t q ¡1 r¡BA ¡1 C t X BA ¡1 F ¡ Gs.

Dening now

D CA ¡1 B t pBA ¡1 B t q ¡1 BA ¡1 C t ¡ CA ¡1 C t , (4.36a)
H CA ¡1 B t pBA ¡1 B t q ¡1 rBA ¡1 F ¡ Gs ¡ CA ¡1 F, (4.36b)
the system becomes DX H.

The matrix D is symmetric and positive denite, owing to the fact that the sum of the two bilinear forms corresponding to the matrix B veries the discrete inf-sup inequality (c. f. [9,21]).

Adaptive Algorithm

We propose, similarly to Section 2.4.3 and 3.4.3, an adaptive algorithm balancing the error sources considered in the error estimate (4.29), namely the space and time discretization and the linearization error estimate. The algorithm is based on Algorithm 2.18 for the adaptation of the discretization, without considering an initial error estimate, since we assume that the initial conditions can be imposed exactly. Let therefore Γ tm ¡ 1 ¡ γ tm ¡ 0 be user-given weights and crit n , for all 1 ¤ n ¤ N , a chosen threshold that the error on the time interval I n should not exceed. We add an adaptive stopping criterion for the linearization loop. To this purpose, we introduce the user-given weight γ lin ¡ 0, which typically takes values around 0.1. In order to avoid an inuence of the reference parameters t and l on the adaptive stopping of the linearization iterations, we compare the linearization error estimate only to the discretization error estimate of the mechanical part and stop the iterations when it becomes negligible in this comparison. fε pεpu n,0 h qq and σ I pεpu n,0 h q, α n,0 q (d) Linearization loop i. Calculate pu n,i h , p n,i h q and the terms α n,i , fσ I fε pεpu n,i h qq and σ I pεpu n,i h q, α n,i q ii. Calculate the estimators η n,i sp,mec , η n,i tm,mec and η

Examples of elasto-plastic laws used in geomechanics

For isotropic materials, it is useful to describe the material independently of the used basis for the space variables. It is thus common not to express the yield as a function of σ directly, but to use a set of functions of σ I which are invariant under rotations. Some often used variables (cf. e.g. [17]) are the following:

• The principal stresses :

σ 1 , σ 2 , σ 3 R, such that σ 1 ¥ σ 2 ¥ σ 3 and £ σ 1 σ 2 σ 3
is the diagonalization of σ I .

This diagonalization always exists, since the image of σ is the space of real symmetric matrices.

• The invariants of the stress tensor :

I 1 pσ I q tr σ I , I 2 pσ I q 1 2 σ I : σ I , I 3 pσ I q detpσ I q.

• The equivalent von Mises stress depending on the deviator s : σ I ¡ 1 3 tr σ I I:

σ eq pσ I q 3 2 s : s 1 2 ppσ 1 ¡ σ 2 q 2 pσ 2 ¡ σ 3 q 2 pσ 3 ¡ σ 1 q 2 q.
Depending on the properties of the material, some of these variables can be more convenient than others for expressing the model. The rst invariant of the stress tensor expresses for example the volumetric part of the stress. Some materials, like metals, do not show any plastic deformation under hydrostatic stress (which corresponds to σ 1 σ 2 σ 3 , for example applying a uniform pressure on a sphere). Then, it is practical to use the invariant σ eq of the deviator, where the volumetric part of the stress is eliminated, since tr s 0. It is clear, that σ eq 0 corresponds to a hydrostatic stress. In the σ 1 ¡ σ 2 ¡ σ 3 ¡system these points correspond to the line tλp1, 1, 1q t ; λ Ru. The use of the invariant σ eq is convenient to describe materials like metals, which show an elastic behaviour for any hydrostatic stress but react by plastifying whenever the variation between σ 1 , σ 2 and σ 3 (called shear stress ) exceeds a certain level. The set tτ; σ eq pτq ku with k ¡ 0 describes an open cylinder with radius k around the line tτ; σ eq pτq 0u.

The following example is the yield function of metals, called von Mises criterion, which will help understand the slightly more complicated criteria used to modelize soils and rocks.

The von Mises criterion

Let σ y ¡ 0 be the tensile strength of the material. This is the limit from which on the material will undergo plastic deformations under uniaxial tension, and thus easy to determine in experiments. The von Mises yield function is dened as F pσ I q σ eq pσ I q ¡ σ y .

There are dierent ways of visualizing the corresponding yield surface. In the left picture in Figure 4.4 it is displayed in the space of the principal stress. The middle gure shows it in the deviatoric plane, which is the projection of space of principal stress onto the plane tpσ 1 , σ 2 , σ 3 q; σ 1 σ 2 σ 3 0u. In the right picture the criterion is visualized in the σ 1 ¡ σ 2 ¡plane. Unlike metals, geomaterials deform plastically under too high tension, whereas hydrostatic compression makes them more elastic under shear stress. Therefore, the DruckerPrager yield function is usually written in terms of σ eq and I 1 . In the I 1 ¡σ eq ¡ plane the set tτ; F pτq 0u is a half-line, as shown in the right picture of Figure 4.5. Introducing the parameter A ¡ 0 denoting its gradient, the yield function reads F pσ I q σ eq pσ I q AI 1 pσ I q ¡ σ y . We recall that softening corresponds to a decrease of the yield function F for constant σ I and increasing α. The way the yield surface moves depends once again on the material, and is expressed as a function R : R 0 Ñ R 0 of α.

σ y ¡ 1 ¡ ¡ 1 ¡ σ y,ult σy © α α ult © 2 if α ¤ α ult , σ y,ult if α ¡ α ult , (4.42)
where a, σ y,ult , α ult ¡ 0 with σ y,ult σ y are experimentally determined parameters. The yield criterion is then written as f pσ I , αq σ eq pσ I q AI 1 pσ I q ¡ Rpαq, which, in the case α 0, is equal to (4.41). At the cone tip, F is not dierentiable, implying that in the case of a plastic deformation, the direction of the latter is not dened. There are dierent possibilities of handling this singularity. One way consists in smoothening the cone tip, such that it becomes dierentiable.

The drawback of this approach is that it introduces a (non physical) regularization parameter the solution will depend on. In Code_Aster, a supplementary step is added to the steps dened in Section 4.3.5: after verifying if the pair pσ I pεpu n,i¡1 h q, α n,i¡1 q, α n,i¡1 q is elastic in step (1), we verify if it lies on the cone tip in step (1b):

(1b) Since at the cone tip it holds σ eq pσ I q 0 and I 1 pσ I q RpαqA ¡1 , the stress is given by σ I 1 3

Rpα n,i¡1 qA ¡1 I, where α n,i¡1 is approximated by solving

Rpα n,i¡1 ¡ α n¡1 q AI 1 pD : εpu n,i¡1 h ¡ u n¡1 h qq, i.e. by projecting the stress obtained using a linearly elastic stress increment onto the cone tip. We then verify if the corresponding plastic strain

ε p εpu n,i¡1 h q ¡ D ¡1 : σ I εpu n,i¡1 h q ¡ 1 3
Rpα n,i¡1 qA ¡1 D ¡1 : I lies inside the normal cone of the yield surface, i.e. if for any τ verifying F pτq ¤ 0 it holds τ : ε p ¤ 0. If this is the case, we already have the values for σ I pεpu n,i¡1 h q, α n,i¡1 q and α n,i¡1 , and, in order to apply Newton's method (4.12), we have to calculate fσ I fε pεpu n,i¡1 h qq : pεpu n,i h ¡ u n,i¡1 h q. In the case of the DruckerPrager law, it is possible to express this term only using the increment α n,i¡1 ¡ α n¡1 , for details we refer to [45, Section 2.2.3].

If ε p does not lie inside the normal cone the stress tensor violates the nomal condition and is not admissible. We then proceed as in step (2) in Section 4.3.5, assuming that the strain is plastic and that σ I does not lie on the cone tip, meaning that the plastic behavior law depending on f σ F is well-dened.

F pσ I q σ 3 pσ I q ¡ σ 1 pσ I q ¡ ¡mσ c σ 3 pσ I q sσ 2 c , (4.43)
where σ c is the uniaxial compressive strength of the intact rock, linked to the uniaxial tensile strength by σ c 2σ y pm ¡ pm 2 4sq 1 {2 q ¡1 , m ¡ 0 is a material parameter and 0 ¤ s ¤ 1 describes the damage and fracturation state, with s 1 corresponding to the intact rock.

The major dierence between the DruckerPrager and the HoekBrown criterion is that, in the latter, the curve tσ; F pσq 0u in the I 1 ¡ σ eq ¡space is not a half-line, but part of half a parabola, as can be seen in the right parts of Figures 4.5 and 4.6. For the HoekBrown criterion, the damage leading to softening is taken into account in the variables s and m.

Numerical results

In this section we present our numerical results. We start with a poro-elastic test with an analytical solution in order to compare the convergence of the space discretization error estimators to the convergence of the error. We then apply Algorithm 4.7 to an industrial test simulating the excavation of a tunnel, using dierent mechanical behavior laws.

10 0.3 10 0.4 10 0.5 10 0.6 10 0. τ 0.5 ¤ 10 ¡4 . We observe that the error estimators reect the convergence of these errors.

Tunnel excavation

In this test we apply the adaptive algorithm of Section 4.5 to a simulation of a tunnel excavation. The simulation is a three dimensional extension of the numerical test presented in Section 2.5.4 with only one gallery. We perform the same simulation using dierent mechanical behavior laws: one in linear elasticity, one using the DruckerPrager model and nally one using the viscoplastic L&K model [68] based on the HoekBrown criterion.

The domain Ω, illustrated in the left of Figure 4.8, is a 10m ¢ 40m ¢ 40m cutout of the rock, in which a tunnel is digged in the (horizontal) z-direction. The excavation time is t F 1.728 ¤ 10 6 s, corresponding to 20 days. Like in the two-dimensional case, we simulate the excavation by rst calculating the initial total equilibrium of the hole-free geometry, and then letting the boundary condition decrease linearly from the so obtained values to zero on the tunnel wall. At the top and the right of fΩ we set p p 0 4.5MPa and σn Ω σ 0 n Ω , where pσ 0,xx , σ 0,yy , σ 0,zz , σ 0,xy , σ 0,xz , σ 0,yz q p¡16.4MPa, ¡12.7MPa, ¡12.4MPa, 0, 0, 0q. On the rest of fΩ we apply symmetry conditions, i.e. u x 0 on the left, u y 0 at the bottom and u z 0 at the front and back of Ω. The initial conditions are pp¤, 0q p 0 and σp¤, 0q σ 0 , and the source terms f and g are equal to zero.

In all simulations we use the realistic parameter set pigure RFW ! gomprison of the glol error estimtes otined using stti nd n dptive lgorithm for the poroEelsti prolem For each of the three considered mechanical behavior laws we perform four static tests, meaning that we use a xed mesh and time step. Table 4.1 recapitulates for each test the number of tetrahedra and time steps. For test 2 we chose a ner mesh and a larger time step than for test 3, such that the number of space-time unknowns for both tests is similar. Our goal is to see if the spatial and temporal discretization error estimators will capture this dierence. For the two nonlinear mechanical behavior laws, the linearization is done using the Newton method of Section 4.3.3 with the convergence criterium Γ res 10 ¡6 in (4.13). If after ten iterations this convergence criterium is not veried, the current time step is divided by four.

Linear elastic mechanical behavior

We start with the linear elastic behavior law, where we compare the error estimates of the static algorithms to the estimates obtained applying Agorithm 4.7 which, in this case, only equilibrates the discretization errors by adapting the time step and the mesh, since we consider a linear problem. The adaptation parameters are chosen as γ tm 0.8 and Γ tm 1.3.

DruckerPrager behavior law

In the DruckerPrager yield function (4.41) we set A 0.33. We use parabolic softening with the parameters α ult 0.01, σ y 2.57MPa, and σ y,ult 0.57MPa in (4.42).

We can only present the results for the three static computations 2, 3 and 4, since in the rst test the algorithm did not converge. We start by comparing each of the static computations with a computation using the adaptive stopping criterion (4.37) with γ lin 0.1 for the Newton method, but without adapting the discretization. Table 4.2 shows the number of iterations, computation time and obtained error estimates for each of the three considered mesh and time step congurations for both the static and the adaptive case. We observe that the number of iterations is considerably reduced using the adaptive stopping criterion. pigure RFII ! gomprison of the glol error estimtes otined for the hruker!rger ehvior lw using stti nd n dptive lgorithm equilirting the spe nd time disretiztion error the error is quite expensive. This is, on the one hand, due to the size of the local problems for the equilibrated stress reconstructions. In the hybridized problems, there are 18 unknowns on each face of the patch, whereas the uid velocity reconstruction only requires one unknown per face. On the other hand, we integrated the computation into an industrial software, in which the resolution of the global problem is optimized, and whose architecture is not adapted for the use of mixed nite elements. We will detail this point in Section 4.7.3

We then perform a series of simulations, where we adapt the mesh and the time step such that the corresponding discretization error estimators satisfy (4.38) and (4.39). Again we compare the global error estimates obtained in the static computations to the ones from the adaptive algorithm. Figure 4.11 shows that the total number of space-time unknowns is reduced for comparable error estimate, and that the discretization error estimators are balanced using the adaptive algorithm. We observe that, without adaptation, the space discretization error estimators dominates the estimate, which was not the case for the linear elastic behavior law, and that the space discretization error estimate in test 2 is higher than in test 3, althouh the second mesh is ner than the third (cf. Table 4.1).

L&K behavior law

The L&K behavior law [46, 68, 91] is a viscoplastic law developped by EDF to modelize the behavior of rocks, and, in particular, the rock formation destined for the nuclear waste storage project Cigéo in the East of France (cf. Section 1.1). The plasticity criterion in this model is based on the HoekBrown yield function (4.43). We will not detail the model and its Total number of unknowns η sp , static η tm , static η sp , adaptive η tm , adaptive pigure RFIP ! gomprison of the glol error estimtes otined for the v8u ehvior lw using stti nd n dptive lgorithm equilirting the lineriztion nd the spe nd time disretiztion error iterations with the original convergence criterion (4.13) with Γ res 10 ¡6 and the adaptive criterion (4.37) with γ lin 0.1 in algorithms with xed time steps and meshes. In Table 4. 3 we see that we perform less iterations with the adaptive stopping criterion. Compared to the DruckerPrager test, we save less iterations but the error estimates are closer to the ones without the adaptive stopping of the Newton iterations. Again, we see that the computation time is much longer due to the calculation of the error estimate.

We then compare the error estimates obtained in the static computations to the error estimates obtained when using Algorithm 4.7 combining the discretization adaptation and the adaptive stopping criterion. In the left of Figure 4.12 we show the estimates as a function of the total number of unknowns. This number is calculated by summing up the number of unknowns of each performed iteration. Since the second static test requires a lot of iterations compared to the other tests (cf. Table 4.3), this number is higher than for test 4, which uses a ner mesh and more timsteps. Again, the estimated error for a comparable number of unknowns is reduced by the use of the adaptive algorithm. Figure 4.13 shows the distribution of the space discretization error estimators at t t F on the nest mesh of the static algorithms (in the right of Figure 4.8) and on an adaptively rened mesh with 2352 elements, which is displayed in the right of the gure. Again we observe that the distribution is more evenly on the adaptively rened mesh, and that, although the mesh le RFQ ! gomprison of the numer of itertionsD the omputtion time nd the glol estimte of the stti lgorithm nd n lgorithm using the dptive stopping riterion @RFQUA with γ lin 0.1 for the v8u model pigure RFIQ ! gomprison of sptil disretiztion error estimtors fter stti nd n dptive lgorithm using the v8u ehvior lw @leftA nd the dptively re(ned mesh @rightA has less elements, the maximum value is smaller than on the static mesh.

Comment on the implementation

The computation of the error estimators is implemented in the nite element software Code_Aster, but is not yet mature for industrial use. Code_Aster is a nite element code based on Lagrange nite elements, meaning that degrees of freedom can only exist on vertices. The computation of integrals on elements and their following assembling in the global matrix are optimized.

However, if degrees of freedom do not lie on vertices, like it is the case for the nite element spaces we use to reconstruct the equilibrated uxes, the assembly tools provided in the code can not be used. Owing to the hybridization of the local problems, the matrices A and B in (4.34) are block-diagonal and we can use the tools for element computations to calculate them. We only have to assembly the matrix C, whose element contributions are also computed using these tools. This choice is essential for the integration of the error estimators into the logic of the code, and also convenient since the elementary computation routines provide a fast calculation of the barycentric coordinates, which facilitates the computation of the RaviartThomas and BrezziDouglasMarini basis functions. The drawback of the use of these elementwise tools is that the resulting elementary matrices are anticipating their assembly not stored as a table one can access with one pointer. Values in an elementary matrix have to be accessed one by one by calling a routine providing their address. Especially for the stress tensor reconstruction, the recuperation of these matrices is extemely expensive. In order to save computation time, we construct and save for each vertex the matrix D in (4.36a), but also the necessary matrices to construct H and Σ in (4.36b) and (4.35) respectively, which requires a lot of memory. We then only have to compute the vectors F and G in (4.34), using again the elementary computation tools.

The introduction of degrees of freedom on mesh faces in Code_Aster is planned in the context of another PhD thesis at EDF R&D. The computation of the equilibrated ux reconstructions can then be reimplemented and industrialized in order to reduce the computation time. An-4.7. NUMERICAL RESULTS other possible improvement is to use parallelization, since the local problems are independent of each other.

Conclusion

The results we obtain balancing the dierent error sources in the excavation test are promising.

In our tests the adaptive stopping of the linearization iterations has an insignicant eect on the discretization error, while the number of total iterations is reduced by up to 80%. In each of the three test series (corresponding to the three considered mechanical behavior laws) the number of unknowns in the computation for a comparable error estimate was reduced by using adaptive algorithms.

A.1. INTRODUCTION

A.1 Introduction

We consider the linear elasticity problem on a simply connected polygon Ω R 2 : ¡∇ ¤ σpuq f in Ω, where u : Ω Ñ R 2 the displacement, and f : Ω Ñ R 2 the volumetric body force. The Cauchy stress tensor σ is given by Hooke's law σpuq λ trpεpuqqI 2 2µεpuq, where λ and µ are the Lamé parameters, and the symmetric gradient εpuq 1 2 pp∇uq T ∇uq describes the innitesimal strain.

In many applications, this problem is approximated using H 1 -conforming nite elements. It is well known that, in contrast to the analytical solution, the resulting discrete stress tensor does not have continous normal components across mesh interfaces, and its divergence is not locally in equilibrium with the source term f on mesh cells. In this paper we propose an a posteriori error estimate based on stress tensor functions which are reconstructed from the discrete stress tensor such that they verify both of the above properties.

A.2 Setting

We denote by L 2 pΩq the space of square-integrable functions taking values in R, and by p¤, ¤q and ¤ the corresponding inner product and norm. H 1 pΩq stands for the Sobolev space com-A.4. STRESS TENSOR RECONSTRUCTIONS Proof. From (A.4) and the symmetry of σpu ¡ u h q, we infer that u ¡ u h en ¢ σpu ¡ u h q, εpu ¡ u h q u ¡ u h en

¤ µ ¡ 1 {2 £ σpu ¡ u h q, εpu ¡ u h q ∇pu ¡ u h q ¤ µ ¡ 1 {2
sup vH 1 0 pΩq; ∇v 1 pσpu ¡ u h q, ∇vq.

(A.7)

Fix v H 1 0 pΩq, such that ∇v 1. Using the fact that u veries (A.2), and inserting p∇ ¤ σ h , vq pσ h , ∇vq 0 into the term inside the supremum yields pσpu ¡ u h q, ∇vq pf, vq ¡ pσpu h q, ∇vq pf ∇ ¤ σ h , vq pσ h ¡ σpu h q, ∇vq.

¤ σ h , vq § § ¤ § § § Ţ T h pf ∇ ¤ σ h , v ¡ Π 0 T vq T § § § ¤ Ţ T h h T π f ∇ ¤ σ h T ∇v T ,
whereas the CauchySchwarz inequality applied to the second term directly yields § § pσ h ¡ σpu h q, ∇vq § § ¤ Ţ T h σ h ¡ σpu h q T ∇v T .

Inserting these results in (A.7) and again applying the Cauchy-Schwarz inequality yields the result.

A.4 Stress Tensor Reconstructions

The set of vertices of the mesh T h is denoted by V h ; it is decomposed into interior vertices V int h and boundary vertices V ext h . For all a V h , T a is the patch of elements sharing the vertex a, ω a the corresponding open subdomain in Ω, n ωa its unit outward normal vector, and ψ a the piecewise ane hat function which takes the value 1 at the vertex a and zero at all the other vertices. For all T T h , V T denotes the set of vertices of T and h T its diameter.

From now on, u h denotes the solution of (A.3). The goal is to minimize the error estimate

S aw

h pω a q : tτ h Hpdiv, ω a q P k 2 sym pTq | τ h | T S aw T dT T a u,

V aw h pω a q : tv h rL 2 pω a qs 2 | v h | T V aw T dT T a u.

Let now k : p ¡ 1. We need to consider subspaces where a zero normal component is enforced on the stress tensor. Since the boundary condition in the exact problem prescribes the displacement and not the normal stress, we distinguish the case whether a is an interior vertex or a boundary vertex. For a V int h , we set Σ a h : tτ h S aw h pω a q | τ h n ωa 0 on fω a , τ h pbq 0 db V ext ωa u, (A.12a)

V a h : tv h V aw h pω a q | pv h , zq ωa 0 dz RM u, (A.12b) with V ext ωa V h fω a , and for a V ext h , we set Σ a h : tτ h S aw h pω a q | τ h n ωa 0 on fω a zfΩ, τ h pbq 0 db V ext ωa u, (A.13a)

V a h : V aw h pω a q, (A.13b) with V ext ωa V h pfω a zfΩq. As argued in [10], the nodal degrees of freedom on fω a are set to zero if the vertex separates two edges where the normal stress is zero.

Construction A.3 (AW stress reconstruction). Find σ a h Σ a h and r a h V a h such that for all pτ h , v h q Σ a h ¢ V a h , pσ a h , τ h q ωa pr a h , ∇ ¤ τ h q ωa pψ a σpu h q, τ h q ωa , (A.14a) p∇ ¤ σ a h , v h q ωa p¡ψ a f σpu h q∇ψ a , v h q ωa .

(A.14b)

Then, extending σ a h by zero outside ω a , set σ h : °aV h σ a h .

Using the denitions (A.12) and (A.13), the formulation (A.14) is equivalent to (A.9). For interior vertices, the source term in (A.14a) has to verify the Neumann compatibility condition On a patch ω a the global space S afw h pω a q is the subspace of Hpdiv, ω a q composed of functions belonging piecewise to S afw T . The spaces V afw h pω a q and Λ h pω a q consist of functions lying piecewise in V T and Λ T respectively, with no continuity conditions between two elements.

As for the previous construction, we dene subspaces with zero normal components enforced on the stress tensor, and distinguish between interior and boundary vertices that for all pτ h , v h , µ h q Σ a h ¢ V a h ¢ Λ a h , pσ a h , τ h q ωa pr a h , ∇ ¤ τ h q ωa pλ a h , τ h q ωa pψ a σpu h q, τ h q ωa , (A.18a) p∇ ¤ σ a h , v h q ωa p¡ψ a f σpu h q∇ψ a , v h q ωa , (A.18b) pσ a h , µ h q ωa 0.

(A.18c)

Then, extending σ a h by zero outside ω a , set σ h : °aV h σ a h .

pigure eFI ! ilement digrms for pS aw T , V aw T q with k 1 @leftA nd pS afw

2 1

 2 , contexte et structure du manuscrit . . . . . . . . . . . 1.2 Le problème de Biot . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Reconstruction de ux et estimateurs d'erreur pour le problème poro-élastique de Biot en 2D . . . . . . . . . . . . . . . . . . . . . . 1.4 Reconstruction de contraintes et estimateurs d'erreur pour l'élasticité non linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Algorithmes adaptatifs pour des problèmes poro-mécaniques en 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Le but de cette thèse est de développer des algorithmes adaptatifs pour la poro-mécanique et de les implémenter dans Code_Aster, le code éléments nis d'EDF R&D. Dans l'introduction, nous allons commencer par présenter les motivations et le contexte de la thèse, ainsi que le problème poro-mécanique. Ensuite, nous décrivons les idées et contributions principales de chaque chapitre de cette thèse. Pour conclure, nous allons donner quelques perspectives de ces travaux. .1. MOTIVATIONS, CONTEXTE ET STRUCTURE DU MANUSCRIT 1.1 Motivations, contexte et structure du manuscrit Nous commençons par présenter le contexte industriel de cette thèse. Ensuite nous expliquons ce que nous entendons par l'expression algorithmes adaptatifs et comment ceux-ci peuvent répondre aux problèmes rencontrés dans le contexte industriel. Nous résumons ensuite les résultats obtenus lors de la thèse et décrivons la structure du manuscrit.

Figure

  Figure IFI ! hém du prinipe des instlltions de gigéoD le entre industriel de stokge géologiqueF Source : ANDRA

4 1

 4 Figure IFP ! ixemple de géométries omplexe X tunnel et lvéolesF Source : ANDRA

2 i f spatial estimator scales vary considerably r e m e s h i fη

 2 Figure IFQ ! hém d9un lgorithme dpttif ve un ritère d9rrêt pour le solveur de xewton @ligne RA et l9dpttion de l disrétistion @lignes IH ! ISA en fontion des estimteurs η lin D ηsp et ηtm qui dénotent respetivement l9estimteur de l9erreur de linéristion et de l disrétistion sptile et temporelleF t et τn sont l9instnt et le ps de tempsF

1. 2

 2 Nous commençons au Chapitre 2 par développer une estimation d'erreur a posteriori pour des problèmes poro-élastiques linéaires en 2D. La reconstruction du ux hydraulique est celle présentée pour des problèmes de Darcy, que nous adaptons à la reconstruction du tenseur symétrique des contraintes en utilisant les espaces d'éléments nis mixtes d'ArnoldWinther qui imposent la symétrie du tenseur fortement. Dans le Chapitre 3, nous nous focalisons sur la partie mécanique et considérons l'élasticité non linéaire. Nous présentons également une nouvelle stratégie de reconstruction du tenseur de contraintes en utilisant les éléments d'ArnoldFalkWinther, qui orent des avantages d'implémentation en anticipant le cas 3D. Nous démontrons que, sous certaines hypothèses sur la loi de comportement, l'estimation multipliée par une constante est également une borne inférieure de l'erreur. Enn, dans le Chapitre 4, nous revenons au problème poro-mécanique et considérons d'un point de vue numérique des lois de comportement mécanique élasto-plastiques en 3D. Dans la partie restante de ce chapitre d'introduction, nous allons d'abord présenter brièvement les équations qui décrivent le comportement des sols. Dans les trois sections suivantes, nous résumons les Chapitres 2, 3, et 4. Nous nissons par donner quelques perspectives à nos travaux. Le problème de Biot Le problème de Biot modélise des milieux poreux constitués d'une part de grains solides formant un squelette déformable et un réseau de pores connectés entre eux, et d'autre part d'un uide occupant ces pores. Plusieurs matériaux présentent ce type de structure, par exemple le béton, certains matériaux organiques comme le bois ou ce qui nous intéresse dans cette thèse les géomatériaux comme les sols et les roches. Le modèle quasi-statique est basé sur deux principes physiques : l'équilibre mécanique et la conservation de la masse du uide. Soit

Figure

  Figure IFR ! ixemple d9un pth utour d9un noeud a du millge (gauche)D et l fontion hpeu d9un noeud (droite)F

Figure 2 .

 2 5. An important property a ψ a pigure PFS ! rt funtion for our constructions is the following partition of unity:

  pigure PFV ! veftX gomprison of the pressure in the qurter (veEspot prolem long the digonl etween the three lgorithms in le PFQF ightX en dpted mesh t time t 10 dys time. The last value of the analytical mechanical error under space renement is due to the error in time discretization which starts playing a role for the nest mesh. We observe that the orders of magnitude of the mechanical and hydraulic part are comparable in this test, and that the eectivity index is dominated by the hydraulic part under space renement, and by the mechanical part under time renement.

  quadrilateral, vertical cutout of the rock, in which two galleries are digged time-delayed in the z-direction, rst left, then right. Both excavations take 17.4 days (1.5¤10 6 s) and the second one starts 11.6 days (10 6 s) after the end of the rst one. For both excavations we rst calculate the initial total equilibrium of the hole-free geometry. Then the digging is simulated by linearly decreasing boundary conditions on the tunnel (convergence connement method). These are of Neumann type for the mechanical part and of Dirichlet type for the hydraulic part and start with the total stress measured at the equilibrium state and the pressure p 0 4.7MPa. Homo-pigure PFIH ! ptil disretiztion error estimtors t tF on (xed mesh @leftD PWDPUS dofsA nd on the lst mesh of n dptive lgorithm @rightD ISDHTR dofsA geneous Dirichlet boundary conditions for the y-component of the displacement and p p 0 are imposed on the bottom of Ω (except for the tunnel parts), while on the top, the left and the right sides of Ω, we set θn = θ ref n with pθ ref,xx , θ ref,yy , θ ref,xy q : p¡11MPa, ¡15.4MPa, 0q and p p 0 . These boundary conditions have to be taken into account for the stress reconstruction (cf. Remark 2.7) and the a posteriori error estimate (cf. Remark 2.15). The initial

Figure 2 .

 2 11, we plot the evolution of the error estimators in the two computations circled in the left graphic. Each mark stands for an iteration and shows the error estimate e n of the current time interval divided by τ n . For the plain algorithm, an iteration is gomprison etween stti lgorithm with (xed mesh nd time step nd the dptive lgoE rithm PFIV for the exvtion dmge test equal to a time step. The adaptive algorithm recalculates the solution at a time step whenever the error estimate lies over crit n (illustrated by the dashed line) by rening τ n or the mesh

  reconstruction and a posteriori error estimation for nonlinear elasticity This chapter consists of an article in preparation, written with Michele Botti. Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.2.1 Continuous setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Discrete setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Equilibrated stress reconstruction . . . . . . . . . . . . . . . . . . . 3.3.1 Patchwise construction in the ArnoldFalkWinther mixed nite element spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Discretization and linearization error stress reconstructions . . . . . 3.4 A posteriori error estimate and adaptive algorithm . . . . . . . . . 3.4.1 Guaranteed upper bound . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Distinguishing the dierent error components . . . . . . . . . . . . . 3.4.3 Adaptive algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.4 Local and global eciency . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 L-shaped domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Linear elasticity model . . . . . . . . . . . . . . . . . . . . . . . . . . HenckyMises model . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Notched specimen plate . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.1. INTRODUCTION

  The stress-strain law σ : Ω ¢ R d¢d sym Ñ R d¢d sym is assumed to satisfy regularity requirements inspired by [25, 86, 87]. Problem (3.1) describes the mechanical behavior of soft materials [103] and metal alloys [88]. Examples of stress-strain relations of common use in the engineering practice are given in Section 3.2.

  Winther mixed nite element spaces [9], leading to weakly symmetric tensors . In[94] this reconstruction is compared to a similar reconstruction introduced in [95] using the Arnold Winther nite element spaces [10], yielding a symmetric tensor, and very good agreement was observed while saving substantial computational eort. In Section 3.3 we apply this reconstruction to the nonlinear case by constructing two stress tensors: one playing the role of the discrete stress and one expressing the linearization error. They are obtained by summing up the solutions of constrained minimization problems on cell patches around each mesh vertex, so that they are Hpdivq-conforming and the sum of the two reconstructions veries locally the mechanical equilibrium (3.1a). The patch-wise equilibration technique was introduced in [28, 38] for the Poisson problem using the RaviartThomas nite element spaces. In[48] it is extended to linear elasticity without any symmetry constraint by using linewise Raviart Thomas reconstructions. Elementwise reconstructions from local Neumann problems requiring some pre-computations to determine the normal uxes to obtain an equilibrated stress tensor can be found in[4, 36, 73, 84], whereas in [82] the direct prescription of the degrees of freedom in the ArnoldWinther nite element space is considered.Based on the equilibrated stress reconstructions, we develop the a posteriori error estimate in Section 3.4 and prove that this error estimate is ecient, meaning that, up to a generic constant, it is also a local lower bound for the error. The idea goes back to [89] and was advanced amongst others by[START_REF]A posteriori error estimate[END_REF] 69, 71, 93] for the upper bound. Local lower error bounds are derived in[28, 38, 54, 56, 75]. Using equilibrated uxes for a posteriori error estimation oers several advantages. The rst one is, as mentioned above, the possible distinction and comparison of error components by expressing them in terms of uxes. Secondly, the error upper bound is obtained with fully computable constants. In our case these constants depend only on the parameters of the stress-strain relation. Thirdly, since the estimate is based on the discrete stress (and not the displacement), it does not depend on the mechanical behavior law

52 3 . 2 .

 32 SETTING as the Poisson problem in [27,56], linear elasticity in [48] and the related Stokes problem in [33] and recently in 3D in [57].

  satisfy the above assumptions. Hyperelasticity is a type of constitutive model for ideally elastic materials in which the stress is determined by the current state of deformation by deriving a stored energy density function Ψ : R d¢d sym Ñ R, namely σpτq : fΨpτq fτ .

Example 3 . 2 ( 2 trpτq 2 µ

 3222 Linear elasticity). The stored energy density function leading to the linear elasticity model is Ψ lin pτq : λ trpτ 2 q, (3.3) where µ ¡ 0 and λ ¥ 0 are the Lamé parameters. Deriving (3.3) yields the usual Cauchy stress tensor σpτq λ trpτqI d 2µτ.

(3. 4 )

 4 Being linear, the previous stress-strain relation clearly satises Assumption 3.1. Example 3.3 (HenckyMises model). The nonlinear HenckyMises model of [61, 80] corresponds to the stored energy density function Ψ hm pτq : α 2 trpτq 2 Φpdevpτqq, (3.5) where dev : R d¢d sym Ñ R dened by devpτq trpτ 2 q ¡ 1 d trpτq 2 is the deviatoric operator. Here, α p0, Vq and Φ : r0, Vq Ñ R is a function of class C 2 satisfying, for some positive constants C 1 , C 2 , and C 3 , C 1 ¤ Φ I pρq α, |ρΦ P pρq| ¤ C 2 and Φ I pρq 2ρΦ P pρq ¥ C 3 dρ r0, Vq. (3.6) We observe that taking α λ 2 d µ and Φpρq µρ in (3.5) leads to the linear case (3.3). Deriving the energy density function (3.5) yields σpτq λpdevpτqq trpτqI d 2μpdevpτqqτ,(3.7)with nonlinear Lamé functions μpρq : Φ I pρq and λpρq : α ¡ Φ I pρq. Under conditions(3.6) it can be proven that the previous stress-strain relation satises Assumption 3.1.In the previous example the nonlinearity of the model only depends on the deviatoric part of the strain. In the following model it depends on the term τ : C τ.

Example 3 . 4 (

 34 An isotropic reversible damage model). The isotropic reversible damage model of[34] can also be interpreted in the framework of hyperelasticity by setting up the energy density function asΨ dam pτq : p1 ¡ Dpτqq 2 τ : C τ ΦpDpτqq,(3.8)where D : R d¢d sym Ñ r0, 1s is the scalar damage function and C is a fourth-order symmetric and uniformly elliptic tensor, namely, for some positive constants C ¦ and C ¦ , it holdsC ¦ |τ| 2 d¢d ¤ C τ : τ ¤ C ¦ |τ| 2 d¢d , dτ R d¢d .(3.9)3.2. SETTINGThe function Φ : r0, 1s Ñ R denes the relation between τ and D by fφ fD 1 2 τ : C τ. The resulting stress-strain relation reads σpτq p1 ¡ DpτqqC τ.

( 3 .

 3 10)If there exists a continuous function f : r0, Vq Ñ ra, bs for some 0 a ¤ b ¤ 1, such that s r0, Vq Ñ sf psq is strictly increasing and, for all τ R d¢d sym , Dpτq 1 ¡ f pτ : C τq, the damage model constitutive relation satises Assumption 3.1.

( 3 .

 3 11)Owing to the growth assumption (3.2a), for all v, w H 1 pΩq, the form apv, wq : pσp∇ s vq, ∇ s wq (3.12a) is well dened and, from equation (3.11), we can derive the following weak formulation of (3.1):

  σ a h by zero outside ω a , set σ h : °aV h σ a h .For interior vertices, the source term in(3.19b) has to verify the Neumann compatibility condition p¡ψ a f σp∇ s u h q∇ψ a , zq ωa 0 dz RM d . (3.20) Taking ψ a z as a test function in (3.14), we see that (3.20) holds and we obtain the following result.

  .26) by applying (3.2b) and the Korn inequality for the left inequality, and the CauchySchwarz inequality and (3.2a) for the right one. In our case it holds C K c 2, owing to (3.1b). 60 3.4. A POSTERIORI ERROR ESTIMATE AND ADAPTIVE ALGORITHM Theorem 3.10 (Basic a posteriori error estimate). Let u be the analytical solution of (3.13)

  .27) Remark 3.11 (Constants C gro and C mon ). For the estimate to be computable, the constants C gro and C mon have to be specied. For the linear elasticity model (3.4) we set C gro : 2µ dλ and C mon : c 2µ, whereas for the HenckyMises model (3.7) we set C gro : 2μp0q d λp0q and C mon : 2μp0q. For the damage model (3.10) we take C gro : C ¦ and C mon : c C ¦ , where C ¦ and C ¦ are the constants appearing in (3.9). Following [94], we obtain a sharper bound in the case of linear elasticity, with µ ¡ 1 {2 instead of c 2C gro C ¡3 mon in (3.27).

  term of the right hand side of (3.29) we obtain, using (3.21) on each T T h to insert Π 0 T v, which denotes the L 2 -projection of v onto P 0 pTq, the CauchySchwarz inequality and the Poincaré inequality on simplices, § §

  theorem also holds for σ h : σ k h,disc σ k h,lin , where σ k h,disc and σ k h,lin are dened in Constructions 3.7 and 3.8 and we obtain the following result.

64 3 . 4 .

 34 A POSTERIORI ERROR ESTIMATE AND ADAPTIVE ALGORITHM Dene, for each T T h , pη k U,T q 2 :

  sym , |σpτ q ¡ σpηq| d¢d ¤ C Lip |τ ¡ η| d¢d . (Lipschitz continuity) (3.40) Notice that the three stress-strain relations presented in Examples 3.2, 3.3, and 3.4 satisfy the previous Lipschitz continuity assumptions. Owing to the denition of the functional R T T

  indicators, and second, to propose a stopping criterion for the Newton iterations based on the linearization error estimator. All the triangulations are conforming, since in the remeshing progress hanging nodes are removed by bisecting the neighboring element.

Figure 3 .

 3 Figure 3.2 compares the distribution of the error and the estimators on the initial and two adaptively rened meshes. The error estimators reect the distribution of the analytical error, which makes them a good indicator for adaptive remeshing. Figure 3.3 shows the global

µp3λ 2µq λ µ 3 ¤

 3 Left: xothed speimen plteF Right: nixil trtion urve experiments, yielding a relation between σ ii and ii for a space direction x i . Since we only consider isotropic materials, we can choose i 1. From this curve one can compute the nonlinear Lamé functions of (3.7) and the damage function in (3.10). Although the uniaxial relation is the same, the resulting stress-strain relations will be dierent. In our test, we use the σ 11 11 relation indicated in the right of Figure 3.5 with σ c 3 ¤ 10 4 Pa, E 10 8 Pa, E res 3 ¤ 10 7 Pa,

  stress σ and the elastic part ε e satisfy Hooke's law. There are dierent ways of classifying materials. In this work we only consider classical plasticity phenomena, meaning that the time and the speed of the deformation play an ancillary role; the deformation depends primarily on the current loading and possibly the dierent states the material went through. Counterexamples are phenomena of creep and fatigue or dynamic plasticity, involving additional strain parts in (4.1).

Figure 4 .

 4 Figure 4.2 shows, again on a one dimensional example, how softening changes the properties of a material. The left picture corresponds to the reaction of a material point under increasing

  if there exists a real forth-order tensor D ¦ such that 9 σ I D ¦ : 9 ε p D : ε e ) is proven for example in [37]. The well-posedness of the pure mechanical problem for a class of hardening elasto-plastic materials has been shown in [23]. As discussed in Section 4.1.2, for softening materials the problem does in general not have a unique solution. In this case, it is possible to add regularization parameters in the mechanical behavior law, as for example in [58].

4. 4 .

 4 EQUILIBRATED FLUX A POSTERIORI ERROR ESTIMATE pigure RFQ ! he virt!homs nd ernold!plk!inther (nite elements for k 1 4.4.1 Quasi-static ux reconstruction

Following

  the reconstruction in Section 2.3.1 we use the RaviartThomas nite elements of order k ¡ 1, illustrated in the left part of Figure 4.3 for k 1. On one tetrahedron T T n h they are given by W T P k¡1 pTq xP k¡1 pTq.

4 and 3 . 3 ,

 33 using the ArnoldFalkWinther nite elements illustrated in the right part of Figure 4.3. As in Section 3.3, we compute two stress tensor reconstructions: one corresponding to the discrete stress and one expressing the linearization error in terms of a stress tensor. This distinction will be useful in Section 4.4.3, where we estimate the dierent error sources separately. We recall the corresponding mixed nite element spaces and the reconstructions for the sake of completeness. The ArnoldFalkWinther spaces on one tetrahedron are the extension to tensors of the BrezziDouglasMarini space for the stress, and Lagrange elements for the Lagrange multipliers imposing the equilibrium and the weak symmetry of the stress tensor. They are dened as

( 4 .Construction 4 . 3 (

 443 19f )Then the two stress tensor reconstructions are dened as in Section 3.3.2: Discrete stress reconstruction). For each

4. 4 . 4

 44 Hybridization of the local problemsThe local problems (4.16), (4.20) and (4.21) result in saddle point problems, whose algebraic resolution is in general more costly than for symmetric positive denite systems of similar size. Following the idas of [59], the systems can be tranformed into symmetric positive denite systems by removing the inter-element continuities of the ux variables and imposing them using Lagrange multipliers. The system can then be reduced using static condensation, since ux and potential variables are discontinous from one element to the other.

  h n F F , ξ h q F 0.

( 4 .

 4 33d) Note that, in this formulation, we impose the homogenous Neumann condition of (4.19a) and (4.19d) weakly using the space Ξ a h and the denition of F a , which has the practical advantage that for any element T T n h the spaces Σ a h |T are the same, independently of the vertex a.

Algorithm 4 . 7 (

 47 Adaptive algorithm).1. Initialisation: choose an initial mesh T 0 h , an initial time step τ 0 , and set t 0 : 0 2. Time loop (a) Set n : n 1, T n h : T n¡1 h , and τ n : τ n¡1 (b) Calculate the initial guess pu n,0 h , p n,0 h q 94 4.6. EXAMPLES OF ELASTO-PLASTIC LAWS USED IN GEOMECHANICS (c) Calculate fσ I

( 4 .

 4 41)In the principal stress space, the yield surface is an open cone, as can be seen in the left picture of Figure4.5. As mentioned earlier, we only consider softening behavior in this work.

F

  pσq 0 pigure RFS ! he hruker!rger yield riterion

F

  pσq 0 pigure RFT ! he roek!frown yield riterion displyed from left to right in the prinipl stress speD on the plne tσ; σeqpσq 0u nd in the I1 ¡ σeq¡ plne 4.6.3 The HoekBrown criterion The HoekBrown criterion [6365] is one of the most used plasticity criteria in mining engineering. Its yield function, diplayed in Figure 4.6 is dened as

g 0 .

 0 ppt, x, yq ¡ cosp¡πtq sinpπxq sinpπyq sinpπzq with κ 1, c 0 0, and the stress-strain relation σ I pεq 2µε λ trpεqI, with the Lamé parameters µ λ 0.4. The resulting source terms f and g in (4.7) are given by f p3.6π 2 sinp¡πtq ¡ π cosp¡πtqq The exact solution is imposed as Dirichlet condition on the boundary fΩ. To evaluate the convergence of the space discretization error estimators we compare the mechanical estimators to the error e mec : he geometry of the rok with the hole for the tunnel nd the orsest nd (nest mesh for the stti omputtions

Figure 4 .

 4 Figure 4.7 shows the comparison of the spatial discretization error estimators to the mechanical and hydraulic errors obtained on a series of uniformly rened meshes with the xed time step

E 3 .

 3 62 ¤ 10 9 Pa, ν 0.3, b 0.6, c 0 9 ¤ 10 ¡11 Pa ¡1 , κ 10 ¡17 m 2 Pa ¡1 s ¡1 , where the Young modulus E and the Poisson ratio ν lead to the Lamé parameters µ 1.4 ¤ 10 9 Pa and λ 2.1 ¤ 10 9 Pa. The reference parameters are chosen as l 40m and t 86400s, corresponding to one day.

Figure 4 .

 4 Figure 4.9 shows that in our tests the use of the adaptive algorithm reduces signicantly the number of space-time unknowns and equilibrates the two discretization error sources. In the left part of Figure4.10 we compare the distribution of the spatial discretization error estimators at t t F on the nest static mesh, which is the right mesh in Figure4.8, to the ones of the last adaptive test (i.e. the one with the lowest error estimate). The corresponding

102 4 .

 4 7. NUMERICAL RESULTSpigure RFIH ! gomprison of sptil disretiztion error estimtors fter stti nd n dptive lgorithm using liner elti ehvior lw @leftA nd the dptively re(ned mesh @rightA mesh has 4356 elements and is displayed in the right of Figure4.10. The error estimators are distributed more evenly on the adaptively rened mesh, especially around the tunnel wall the values of the estimators are smaller than on the xed mesh.

  parameters here and refer to[46].As in the DruckerPrager test, we start by comparing the number of performed Newton 10

  Such equilibrated-ux a posteriori error estimates oer serveral advantages. First, error upper bounds are obtained with fully computable constants. Second, polynomial-degree robustness can be achieved for the Poisson problem in [27, 56], for linear elasticity in [48], and for the related Stokes problem in [33]. Third, they allow one to distinguish among various error components, e.g., discretization, linearization, and algebraic solver error components, and to equilibrate adaptively these components in the iterative solution of nonlinear problems [55]. An advantage for more general problems in solid mechanics is that the stress reconstruction is based on the discrete stress (not the displacement) and thus the estimate does not depend on the mechanical behaviour law. We present two stress reconstructions. Both use mixed nite elements on cell patches around mesh vertices, as proposed for the Poisson problem in [28, 38]. The rst one was introduced in [95] and uses the ArnoldWinther (AW) mixed nite element spaces [10] providing a symmetric stress tensor. The second one follows the same approach, but imposing the symmetry only weakly and using the ArnoldFalkWinther (AFW) mixed nite element spaces [9]. Elementwise reconstructions of equilibrated stress tensors from local Neumann problems can be found in [5, 67, 71], whereas direct prescription of the degrees of freedom in the AW nite element space is considered in [82].

(A. 8 )

 8 For the rst term in the right hand side of (A.8) we use (A.5) to insert the mean value Π 0 T v of v on T , the CauchySchwarz inequality, and the Poincaré inequality v ¡ Π 0 T v T ¤ h T π ∇v T on simplexes T T h , and obtain § § pf ∇

(A. 6 )

 6 avoiding global computations. As a result, both of the proposed reconstructions are based on local minimization problems on the patches ω a : σ a h : arg minτ h Σ a h ; ∇¤τ h ψaf τ h ¡ ψ a σpu h q ωa , (A.9)where we dene Σ a h separately for each construction and add a weak symmetry constraint in the second (AFW) construction. The global reconstructed stress tensor σ h is then obtained assembling the local solutions σ a h .A.4.1 ArnoldWinther Stress ReconstructionFor each element T T h , the local AW spaces of degree k ¥ 1 are dened by [10]S awT : tτ P k 2 sym pTq | ∇ ¤ τ rP k pTqu, V aw T : P k pTq, where P k sym pTq denotes the subspace of P k pTq composed of symmetric-valued tensors. Figure A.1 shows the corresponding 24 degrees of freedom for the symmetric stress tensor in the lowest-order case k 1: the values of the three components at each vertex of the triangle, the values of the moments of degree zero and 1 of the normal components each edge, and the value of the moment of degree zero of each component on the triangle. On a patch ω a , the AW mixed nite element spaces are dened as

112A. 4 .

 4 STRESS TENSOR RECONSTRUCTIONS p¡ψ a f σpu h q∇ψ a , zq ωa 0 dz RM. (A.15) Taking ψ a z as a test function in (A.3), we see that (A.15) holds. A.4.2 ArnoldFalkWinther Stress Reconstruction For each element T T h , the local AFW mixed nite element spaces [9] of degree k ¥ 1 hinge on the BrezziDouglasMarini mixed nite element spaces [29] for each line of the stress tensor and are dened by S afw T : P k pTq, V afw T : P k¡1 pTq, Λ T : tµ P k¡1 pTq | µ ¡µ T u.

  ¡∇ ¤ pκρF m q joue le rôle d'une source volumique de uide. Ωq tτ L 2 pΩq; ∇ ¤ τ L 2 pΩq et τ est symétriqueu, φ Hpdiv, Ωq tφ L 2 pΩq; ∇ ¤ φ L 2 pΩqu. Le tenseur σ étant symétrique, on peut remplacer εpvq par ∇v dans (1.7a), et on voit que ¡f L 2 pΩq et g ¡ f t pc 0 p ∇ ¤ uq L 2 pΩq vérient cette dénition.

	1.2. PROBLÈME PORO-ÉLASTIQUE EN 2D
	Plus généralement, on peut considérer des coecients de perméabilité qui sont des tenseurs
	d'ordre deux, ce qui permet de modéliser l'anisotropie du sol. En combinant (1.4) et (1.5), on
	peut réécrire (1.3) comme
	f t pc 0 p b∇ ¤ uq ∇ ¤ φppq g,	(1.6)
	où g Formulation faible

Pour la résolution des équations (1.1) et (1.6) dans un domaine Ω complétées par des conditions au bord et initiales nous supposons que les termes source f et g sont à presque tout temps t p0, t F q des fonctions de carré intégrable dans Ω, et que la norme correspondante dans Ω est également de carré intégrable dans l'intervalle de temps. Pour écrire la formulation variationelle nous supposons qu'à presque tout instant de temps les variables primales up¤, tq et pp¤, tq vivent dans les espaces de Sobolev H 1 pΩq et H 1 pΩq respectivement, et que les fonctions t Þ Ñ up¤, tq et t Þ Ñ pp¤, tq admettent une dérivée au sens faible. On cherche donc u H 1 p0, t F ; H 1 pΩqq et p H 1 p0, t F ; H 1 pΩqq vériant les conditions initiales, telles que pour presque tout t p0, t F q et pour toutes fonctions test v L 2 p0, t F ; H 1 pΩqq et q L 2 p0, t F ; H 1 pΩqq, » Ω σpεpuq, pq : εpvq dx

» Ω f ¤ v dx, (1.7a) » Ω f t pc 0 p b∇ ¤ uqq dx ¡ » Ω φppq ¤ ∇q dx » Ω gq dx.

(1.7b) En regardant la première intégrale dans (1.7b) on voit que, contrairement à u et p, il n'est pas nécessaire que les fonctions test soient dérivables en temps.

Le tenseur de contraintes σ σpεpuq, pq et la vitesse φ φppq résultants de (1.7) vérient deux propriétés importantes : premièrement, ils admettent pour presque tout t p0, t F q une divergence au sens faible. On a donc σ H s pdiv, Deuxièmement, ils satisfont les deux équations d'équilibre (1.1) et (1.6). Ces deux propriétés peuvent être facilement vériées, en rappelant qu'une fonction w L 2 pΩq est la divergence au sens faible de v L 2 pΩq si ³ Ω wψdx ¡ ³ Ω v ¤∇ψdx pour toute ψ DpΩq H 1 0 pΩq.

  IFS ! higrmmes des éléments (nis d9ernold!inther et de virt!homs de degré le plus s à travers les bords du patch. De telle manière, la somme des solutions locales est de nouveau Hpdivq-conforme. L'espace discret utilisé est l'espace de Raviart-Thomas de plus bas degré, dont les degrés de liberté sont représentés dans la partie droite de la Figure 1.5 avec ceux de .3. PROBLÈME PORO-ÉLASTIQUE EN 2D Cependant, la reconstruction des ux et notamment du tenseur de contraintes coûte cher en termes de temps de calcul. Ceci est d'une part dû au fait qu'on compare le temps de reconstruction avec celui de la résolution du problème initial par Code_Aster qui est, par le caractère industriel du logiciel, améliorée et optimisée depuis plusieurs années. Mais ces longs temps de calcul s'expliquent surtout par la taille et la nature des problèmes locaux à résoudre : les éléments d'ArnoldWinther sont très riches et, contrairement à d'autres familles d'éléments nis mixtes, l'hybridation des problèmes locaux n'est pas possible, à cause des degrés de liberté sur les sommets des éléments. Cette technique est utilisée pour transformer le problème de type point selle résultant de l'utilisation de méthodes mixtes en un problème déni positif, ce qui rend la résolution du problème plus facile et rapide en pratique. Une possibilité an d'améliorer les temps de calcul est donc de relâcher les contraintes de symétrie, ce qui sera discuté dans la section suivante.

	1.4 Reconstruction de contraintes et estimateurs d'erreur pour
	l'élasticité non linéaire
	Comme nous l'avons vu, l'originalité de cette thèse est le développement de reconstructions
	équilibrées du tenseur de contraintes mécaniques pour obtenir des estimateurs d'erreur a poste-
	riori pour le problème de Biot. Pour simplier l'analyse, nous nous focalisons dans le Chapitre
	3 sur des problèmes purement mécaniques. Plus précisement, nous nous plaçons dans un cadre
	[28, 38, 48, 56, 75].
	Il existe de nombreuses autres méthodes pour obtenir une estimation de l'erreur, une vue
	d'ensemble peut être trouvée dans [2, 3, 31, 32]. Quelques exemples hormis les estimateurs par
	résidus et équilibrage de ux sont la méthode de résidus équilibrés [1, 3, 71], les estimateurs
	fonctionels [81, 93], les techniques de lissage [116118], ou les méthodes hiérarchiques [3, 14]. Il

u h et p h ont en général des composantes normales discontinues à travers les faces du maillage utilisé pour la discrétisation, et ne vérient pas les équations (1.1) et (1.6). Cette représentation non-physique des contraintes et de la vitesse du uide sont la première raison pour laquelle nous souhaitons calculer des reconstructions de ux, le but étant d'obtenir des ux plus physiques (donc Hpdivq-conformes et équilibrés) et proches des ux discrets. La deuxième (et principale) motivation est l'utilisation de telles reconstructions de ux pour calculer des estimations d'erreur avec les avantages décrits en Section 1.1.

Pour la partie hydraulique, nous nous basons sur la reconstruction de ux introduite dans

[28, 38]

, en utilisant la formulation de

[56]

. L'idée intuitive constiste à chercher dans un espace discret Hpdivq-conforme approprié la fonction la plus proche du ux discret φpp h q qui vérie l'équilibre. Il s'agit donc de résoudre un problème de minimisation sous contrainte. Eectuer à chaque pas de temps ce calcul global sur tout Ω serait trop coûteux ; on se base donc sur des patchs autour des sommets du maillage (on rappelle qu'un patch est l'ensemble des mailles partageant un sommet donné) pour y résoudre le même type de problème de minimisation, et on additionne ensuite les solutions an d'obtenir un champ global. Pour que cette somme approche le ux discret, on multiplie φpp h q par la fonction chapeau du sommet, illustrée dans la Figure 1.4. Sur le bord des patchs on impose des conditions de Neumann homogènes pour pouvoir prolonger la solution par zéro sur le reste du domaine sans introduire des discontinuités 1.3. PROBLÈME PORO-ÉLASTIQUE EN 2D Figure l'espace utilisé pour imposer la divergence du ux. La contribution principale de ce chapitre est d'appliquer le même principe à la partie mécanique pour obtenir une reconstruction équilibrée, Hpdivq-conforme et symétrique du tenseur de contraintes totales. Les espaces discrets de tenseurs symétriques et Hpdivq-conformes sont les espaces d'Arnold-Winther introduits dans [10] pour les maillages triangulaires et dans [8] pour des maillages de tétraèdres en 3D. Dans l'élément triangulaire, il existe trois groupes de degrés de liberté (voir la partie gauche de la Figure 1.5). Le premier contient les moments jusqu'à un certain degré polynomial des composantes normales à travers les arêtes pour garantir leur continuité. Dans le deuxième groupe se trouvent les valeurs des composantes aux sommets du triangle pour garantir la symétrie du tenseur. Le troisième groupe de degrés de liberté contient les moments à l'intérieur du triangle de chaque composante. En utilisant les propriétés des reconstructions de ux ainsi obtenues, nous pouvons établir une borne supérieure garantie (donc sans constantes inconnues) sur la norme duale du résidu de la formulation faible du problème de Biot. Cette mesure d'erreur se prête naturellement au développement d'estimateurs d'erreur de type ux équilibrés. Pour pouvoir bien dénir le résidu et pour comparer les composantes de l'estimation provenant des parties hydraulique et mécanique du problème, nous avons introduit des paramètres d'adimensionnement. Ces paramètres n'interviennent que dans la phase de calcul d'estimateurs : le calcul initial de la solution discrète reste inchangé. Mise en ÷uvre dans Code_Aster Dans un premier temps, les reconstructions et le calcul des estimateurs d'erreur ont été implémentés de façon provisoire dans Code_Aster. Ce code développé par EDF R&D est principalement écrit en Fortran90 (initialement Fortran77), completé par des fonctions C pour réaliser quelques tâches impossibles en Fortran77 (comme l'allocation dynamique), et par des catalogues Python qui gèrent d'une part la communication entre l'utilisateur et le code et d'autre part les diérentes options de calculs élémentaires. Les paramètres d'entrée pour réa-Figure IFT ! gomprison de l distriution des estimteurs de disrétistion sptile sns et ve remillge dpttif @fF pigure PFIHA liser une étude avec Code_Aster sont un maillage et un chier de commandes s'appuyant sur le langage Python, qui enchaîne des commandes pour la lecture du maillage, la dénition du problème physique étudié, sa résolution numérique et le post-traitement. Pendant la résolution du problème exécutée par la commande STAT_NON_LINE, le ux hydraulique et le tenseur de contraintes sont sauvegardés à chaque étape du calcul dans une structure de données. Nous avons créé une commande de post-traitement réutilisant cette structure de données. A chaque pas de temps, elle calcule les estimateurs d'erreur de discrétisation spatiale et temporelle en fonction des ux discrets et les sauvegarde dans la même structure de données. Les termes sources f et g, ainsi que les conditions aux bords, doivent être renseignés dans les sources Fortran. Nous avons validé les estimateurs sur deux cas tests avec une solution analytique, dont un test purement mécanique. Ensuite, nous avons eectué deux tests d'application industrielle : le quart de ve-spot simulant l'injection et la production d'eau dans deux points distincts, et la simulation d'excavations de tunnels dans le contexte du stockage des déchets radioactifs (cf. Figure 1.6). Pour ces deux tests industriels nous avons écrit des chiers de commandes mettant en ÷uvre un algorithme adaptatif. Puisque le calcul des estimateurs est eectué dans une commande séparée de la commande STAT_NON_LINE, nous utilisons une boucle qui, pour chaque pas de temps, résout le problème en prenant le résultat de l'instant de calcul précédent comme donnée initiale, estime ensuite l'erreur, et utilise cette estimation pour adapter le pas de temps et le maillage. Dans la Figure 1.6, nous comparons la distribution des estimateurs de discrétisation spatiale de deux calculs l'un avec un maillage spatial et un pas de temps xé durant le calcul (à gauche), et l'autre en adaptant la discrétisation en fonction des estimateurs de l'erreur spatial et temporel (à droite). On voit que la distribution de l'erreur est signicativement plus homogène si on utilise l'algorithme adaptatif, et nous observons pour les deux tests industriels que, à erreur donné, le nombre de degrés de liberté à erreur donnée peut être réduit considérablement en utilisant un algorithme adaptatif. 1Bibliographie Nous terminons cette section avec un récapitulatif bibliographique. Le problème de Biot a été proposé par Biot [19] et von Terzaghi [110]. ení²ek [111] a démontré qu'il est bien posé (voir aussi [98]). La discrétisation H 1 -conforme utilisée dans Code_Aster utilise les éléments de Taylor-Hood [102], qui ont initialement été introduits pour le problème de Navier-Stokes. Son analyse d'erreur a priori est discutée dans [7779]. Une analyse d'erreur a posteriori pour le problème poro-élastique linéaire a été eectué dans [52, 76], en utilisant une approche d'estimation d'erreur par résidus [12, 13, 104], donc en mesurant les sauts des ux discrets. A notre connaissance, il n'existe pas d'estimation par équilibrage de ux pour le problème de Biot antérieure à celle proposée dans ce manuscrit. Ces méthodes ont été développées à partir de travaux de Ladevèze [69, 70], qui s'appuyent sur le théorème de Prager et Synge [89]. L'idée est de décomposer l'erreur en utilisant une reconstruction de ux qui est plus régulière que le ux discret et qui permet donc de remplacer la solution analytique inconnue dans l'erreur par des termes connus. Plusieurs méthodes pour obtenir une telle reconstruction de ux en résolvant des problèmes de Neumann sur chaque élément ont été développées dans [2, 4, 30, 71, 73, 84] et les références qui y sont contenues. Comme mentionné ci-dessus, nous allons utiliser des problèmes locaux de Neumann sur des patchs, ce qui a été proposé dans existe également des méthodes pour estimer l'erreur dans une certaine fonctionelle en utilisant des techniques de dualité, par exemple dans [16, 51, 83]. Dans le cas des problèmes non stationnaires, il est possible de distinguer les erreurs de dis-crétisation en espace et en temps, et d'adapter le pas de temps et le maillage (cf. Figure 1.3) pour équilibrer ces deux sources d'erreur, comme proposé dans [18, 54, 72, 106]. stationnaire de problèmes d'élasticité non linéaires, dans le but d'étendre les résultats à des problèmes élasto-plastiques. Du fait que l'état d'un matériau avec un comportement plastique à un instant ne dépend pas seulement des forces auxquelles il est sujet à cet instant, mais également de celles qui ont été appliquées avant, ces problèmes sont forcément non stationnaires. Nous les traitons dans le Chapitre 4. Dans le Chapitre 3 nous utilisons une nouvelle reconstruction du tenseur de contraintes, que nous avons développé dans l'Annexe A. Cette reconstruction est obtenue en imposant la contrainte de symétrie faiblement, et non plus fortement comme au chapitre précedent, ce qui facilite l'implémentation et nous permet d'améliorer les temps de calcul. Les deux contributions du Chapitre 3 sont les suivantes :

  and the constant C K,T is bounded by c 2psinpα T {4qq ¡1

CHAPTER 2. BIOT'S PORO-ELASTICITY PROBLEM 29 (cf [66]). Combining (2.10) and (2.12), and accounting for the bounds on the corresponding constants, we infer that, for all

  T,P ptq, t I n , by η n F,T,U ptq : E ¡1 θ n h ¡ θpu hτ , p hτ qptq T ,

	e n ¤	¤ ¥ »	h In Ţ T n	tpη n R,T,U	η n F,T,U ptqq 2 pη n R,T,P	η n F,T,P ptqq 2 udt	1 {2	.	(2.40)
	Proof. Let pv, qq X n . Recalling (2.37a), we have
	e n U pvq	»		θpu, pq ¡ θpu hτ , p hτ q, pvq ¨ptqdt
			» »	In In In	f, v ¨ptq ¡ θpu hτ , p hτ q, pvq ¨ptqdt ¡ f n ∇ ¤ θ n h , vptq loooooooooomoooooooooon θ n loooooooooooooooomoooooooooooooooon h ¡ θpu hτ , p hτ q, pvq ¨ptq © dt,	(2.41)
						T 1 ptq	T 2 ptq
										(2.39a)
				η n F,T,P ptq : t l φ n

h ¡ φpp hτ qptq T .

(2.39b) Theorem 2.13 (Local-in-time a posteriori error estimate). Let pu, pq Y be the weak solution of (2.8) and let pu hτ , p hτ q Y be the discrete solution of (2.15). Let 1 ¤ n ¤ N . Let e n be dened by (2.36) with estimators dened by (2.38) and (2.39). Then the following holds:

2.4. A POSTERIORI ERROR ANALYSIS AND SPACE-TIME ADAPTIVITY

where we have used (2.8a) to pass to the second line and we have inserted p∇ ¤ θ n hτ , vq pθ n hτ , pvqq 0 inside the integral to conclude. For the rst term we have, for a.e. t p0, t F

  When this is not the case, an additional data time-oscillation term appears in the right-hand side of the bound (2.40).

Remark 2.15 (Other types of boundary conditions). If we consider inhomogeneous Neumann boundary conditions θpu, pqn Ω θ N on fΩ N,U fΩ and φppq ¤ n Ω φ N on fΩ N,P fΩ, two more error estimators appear in (2.40). The details of how they are obtained for the hydraulic part are shown in [47] and can be directly applied to linear elasticity. For each T T n h , let E N,U T and E N,P T be the set of edges lying on fΩ N,U and fΩ N,P respectively, and let θ n h and φ n h be the ux reconstructions of Remarks 2.4 and 2.7. Then we set

  During the computation time of 30 days, we set p p 0 in the top right corner and p 4 ¤ 10 5 Pa in the bottom left corner, simulating the production and the injection respectively. The nondimensionalization parameters are l 140m and t 1h. The problem is dominated by hydraulic processes.

		# space-time unknowns	# iterations	η sp	η tm	η sp +η tm
	reference	13,754,520	120	0.204	0.315	0.519
	equivalent	973,620	45	1.14	1.54	2.68
	adaptive	846,174	71	0.462	0.507	0.969

  EQUILIBRATED STRESS RECONSTRUCTION pigure QFI ! ilement digrms for pΣT , VT , ΛT q in the se d q 2 a stress tensor σ h in a suitable (i.e. Hpdivq-conforming) nite element space by summing up these local solutions. The local problems are posed such that this global stress tensor is close to the discrete stress tensor σp∇ s u h q obtained from (3.14), and that it satises the mechanical equilibrium on each element.

	ArnoldFalkWinther mixed nite el-
	ement spaces
	Let us for now suppose that u h solves (3.14) exactly, before considering iterative linearization
	methods such as (3.15) in Section 3.3.2. For the stress reconstruction we will use mixed nite

element formulations on patches around mesh vertices in the spirit of

[94, 95]

. The mixed nite elements based on the dual formulation of (3.1a) will provide a stress tensor lying in Hpdiv, Ωq. A global computation is too expensive for this post-processing reconstruction, so we solve local problems on patches of elements around mesh vertices. The goal is to obtain 3.3.

In [95] the stress tensor is reconstructed in the ArnoldWinther nite element space [10], directly providing symmetric tensors, but requiring high computational eort. In this work, as in [94], we weaken the symmetry constraint and impose it weakly, as proposed in [9]: for each element T T h , the local ArnoldFalkWinther mixed nite element spaces of degree q ¥ 1 hinge on the BrezziDouglasMarini mixed nite element spaces

[29]

for each line of the stress tensor and are dened by

  R 3 u are CHAPTER 3. NONLINEAR ELASTICITY PROBLEM 57 the spaces of rigid-body motions respectively for d 2 and d 3.If a V ext : tτ h Σ h pω a q | τ h n ωa 0 on fω a zfΩu,

	h we set
	Σ a

h (3.18d)

  disc RM d is the unique solution of py k disc , zq ωa ¡pf, ψ a zq ωa pσp∇ s u k h q, ∇ s pψ a zqq ωa dz RM d .

	Construction 3.7 (Discrete stress reconstruction). For each a V h solve (3.19) with u k h in-stead of u h , σp∇ s u k h q instead of σp∇ s u k h q and the source term in (3.19b) replaced by
	¡ψ a f σp∇ s u k h q∇ψ a ¡ y k disc ,

where y k

(3.22)

  and applying the Green theorem, we see that ra h satises (3.44) for faces F F a zF ext h , since the normal components across F of a basis of ΣT span P p pFq. If a V ext h we can proceed in the same way for F F a F ext h to obtain(3.46). Finally, for a V int h it holds pr a h , zq ωa 0 for any z RM d by the denition (3.18b) of V a h , and by (3.54c) it follows that ra

	3.54) is well-posed. Plugging (3.54b) and (3.54c) into (3.53) yields Π ΣT p∇r a h q |T pσ a h ¡ ψ a σp∇ s u k h qq |T . Since the formulations (3.51) and (3.52) are equivalent, we can insert (3.55) and (3.54c) into (3.55) (3.51a) and obtain p∇r a We conclude that ra h lies in M a h . Proof of Theorem 3.16. We start by proving the local approximation property of the discrete h , τ h satises (3.45). stress reconstruction for any T T h

h q ωa pr a h , ∇ ¤ τ h q ωa 0 dτ h Σa h .

Choosing a basis function of

Σa h having zero normal trace across all edges except one edge F

  This solution is imposed as Dirichlet boundary condition on fΩ, together with f 0 in Ω. We perform this test for two dierent stress-strain relations. First on the linear elasticity model

			|T h | I eff
			34	1.00
	10 ¡0.5		84 130	1.01 1.02
			137	1.01
			214	1.05
	10 ¡1.5 10 ¡1	error, unif. estimate, unif. error, adap. estimate, adap.	239 293 429 524 601	1.01 1.00 1.01 1.01 1.01
			801	1.01
		10 2 |T h | number of mesh elements 10 3	1099 1142	1.02 1.02

(3.4)

, where we can compare the error estimate

(3.27)

to the analytical error u ¡ u h en . The second relation is the nonlinear HenckyMises model (3.7), for which we distinguish the discretization and linearization error components and use the adaptive algorithm from Section 3.4.3.

pigure QFQ ! vEshped domin with liner elstiity modelF Left: gomprison of the error estimte @QFPUA nd u ¡ u h en on two series of meshesD otined y uniform nd dptive remeshingF Right: i'etivity indies of the estimte for eh meshD with the meshes stemming from uniform re(nement highlighted in gryF

  QFR ! vEshped domin with renky!wises modelF Left: gomprison of the glol disretiztion nd lineriztion error estimtors on series of meshesD without nd with dptive stopping riterion for the xewton lgorithmF Middle: xumer of xewton itertions without nd with dptive stopping riterion for eh meshF Right: hisretiztion error estimte for uniform nd dptive remeshingF

					norm.	adap.
				0	4	2
	10 ¡2			1	4	2
	10 ¡5			4 2 3	5 5 5	3 3 3	10 ¡1
				5	6	4
	10 ¡11 10 ¡8	η disc , no η lin , no η disc , ad η lin , ad		6 7 8 9 10	6 6 7 7 8	4 4 4 5 5	10 ¡2	η disc , unif. η disc , adap.
		10 2 |T h | number of mesh elements 10 3	10 4	11	9	5		10 2 |T h | number of mesh elements 10 3	10 4
								3 2	μpρq,

corresponding to the Carreau law for elastoplastic materials (see, e.g. [60, 74, 97]), and we set a 1{20, b 1{2, and κ 17{3 so that the shear modulus reduces progressively to approximately 10% of its initial value. This model allows us to soften the singularity observed pigure in the linear case and to validate our results on more homogeneous error distributions. We apply Algorithm 3.14 with γ lin 0.1 and compare the obtained results to those without the adaptive stopping criterion for the Newton solver. In both cases, we use adaptive remeshing based on the spatial error estimators.

  The DruckerPrager criterion . . . . . . . . . . . . . . . . . . . . . . 96 4.6.3 The HoekBrown criterion . . . . . . . . . . . . . . . . . . . . . . . . 98 4.7 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.7.1 Analytical test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.7.2 Tunnel excavation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Linear elastic mechanical behavior . . . . . . . . . . . . . . . . . . . 101 DruckerPrager behavior law . . . . . . . . . . . . . . . . . . . . . . 102 L&K behavior law . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.7.3 Comment on the implementation . . . . . . . . . . . . . . . . . . . . 105 4.7.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 It is a well established concept that a solid body resists to applied forces by deforming itself.

							3.6. CONCLUSIONS
	10 ¡1			10 ¡1			HenckyMises nor. adap.	Damage nor. adap.
	10 ¡5 10 ¡3			10 ¡3 10 ¡5	η disc ,no η lin ,no η disc , ad	0 1 2 3	5 5 6 6	3 4 4 5	6 6 6 7	4 4 5 5
	10 ¡7			10 ¡7	η lin , ad	4 5	6 7	5 5	7 8	6 6
	10 ¡9			10 ¡9		6	6	5	8	6
	10 ¡11	10 3 |T h | number of mesh elements 10 4	10 5	10 ¡11	10 3 |T h | number of mesh elements 10 4
	pigure QFU ! xothed speimen plteF gomprison of the glol disretiztion nd lineriztion error estimtors
	without nd with dptive stopping riterion for the renky!wises model @leftA nd the dmge model @middleAD
	nd omprison of the numer of perfomed xewton itertions @rightAF

measure, the method can be applied to a wider range of behavior laws. Exploring both numerical tests, we have promising results for general plasticity and damage models. These results come at the price of solving local mixed nite element problems at each iteration of the linearization solver. In practice, the corresponding saddle point problems can be transformed into symmetric positive denite systems using the spaces of Section 3.4.4. Furthermore, these matrices (or their decomposition) can be computed once in a preprocessing stage, and only need to be recomputed if one or more elements in the patch have changed due to remeshing. 76 4.6.2

  Then the problem reads: nd pu, p, αq U ¢ P ¢ L 2 p0, t F ; L 2 pΩqq verifying (4.7c), such that 4.3. NUMERICAL SOLUTION for a. e. t p0, t F q , pσ I pεpuq, αq, εpvqq ¡ bpp, ∇ ¤ vq pf, vq if pσ I pεpuq, αq, αq E ¡pf α F q ¡1 f σ F : 9 σ I . if pσ I pεpuq, αq, αq P a.e. in Ω ¢ p0, t F q, (

			dv H 1 0 pΩq,	(4.9a)
	c 0 p 9 p, qq bp∇ ¤ 9 u, qq ¡ pφppq, ∇qq pg, qq σ I pεpuq, αq 9 5 D : 9 εpuq if pσ I pεpuq, αq, αq E pD ¡ D I q : 9 εpuq if pσ I pεpuq, αq, αq P 5	dq H 1 0 pΩq, a.e. in Ω ¢ p0, t F q, (4.9c) (4.9b)
	9 α	0

  T n h and pv h , zq ωa 0 dz RM u, : tµ h L 2 pω a q; µ h|T Λ T dT T n h u, RM : tb a ¢ x; a, b R 3 u is the space of rigid-body motions. If a V n,ext Hpdiv, ω a q; and τ h|T Σ T dT T n h and τ h n ωa 0 on fω a zfΩu,(4.19d)

	4.4. EQUILIBRATED FLUX A POSTERIORI ERROR ESTIMATE
	h	we set
	Σ a h :tτ h
		(4.19b)
	Λ a	(4.19c)

h where

  , which are not written in the same physical units. We use the Young modulus E, a reference time scale t and a reference length scale l . We denote X :L 2 p0, t F ; H 1 0 pΩqq ¢ L 2 p0, t F ; H 1 0 pΩqq, Y pU ¢ P ¢ L 2 p0, t F ; L 2 pΩqqq, and their restrictions to a time interval I n by X n and Y n . Then we can dene the map B : Y ¢X Ñ L 2 p0, t F ; Rq by Bppu, p, αq, pv, qqq : pσ I pεpuq, αq, εpvqq ¡ bpp, ∇ ¤ vq t pc 0 p 9 p, qq bp∇ ¤ 9 u, qq ¡ pφppq, ∇qqq, (4.23) and thus the weak formulation (4.9) is equivalent to nding pu, p, αq Y verifying the initial conditions (4.7c), such that for a. e. t p0, t F q , Bppu, p, αq, pv, qqq pf, vq t pg, qq dpv, qq X, We can then, similarly to Section 2.4.1, dene on every time interval I n the residual of (4.24) for any pu hτ , p hτ , αq Y n and pv, qq X n xRpu hτ , p hτ , αq, pv, qqy X I In Bppu hτ , p hτ , αq, pv, qqq ¡ pf n , vq ¡ t pg n , qqdt (4.25)

	4.4.2 Error measure
	We proceed as in Chapter 2 and introduce reference parameters allowing us to dene the
	residual of (4.9a,4.9b)(4.24)
	with (4.9c,4.9d).
	Then set σ n,i h,lin : °aV h σ a h .
	As shown in Lemma 3.9 together with Assumption 4.1, these two stress tensors lie in Hpdiv, Ωq
	and verify for all 1 ¤ n ¤ N
	¡∇ ¤ pσ k h,disc	σ k h,lin q f.	(4.22)

n ,Xn : »

  qqy X I Rpu hτ , p hτ , αq X I Bppu hτ , p hτ , αq, pv, qqq ¡ pf, vq ¡ pg, qqdt, where the equation in the rst line is shown in Corollary 2.16.

	90		4.4. EQUILIBRATED FLUX A POSTERIORI ERROR ESTIMATE
	4.4.3 A posteriori error estimate
	In this chapter, we directly present the a posteriori error estimate distinguishing dierent error
	sources. For each time step 1 ¤ n ¤ N and each iteration i ¥ 0, let φ n h , σ n,i h,disc and σ n,i h,lin be
	respectively the velocity and stress reconstructions from Constructions 4.2, 4.3 and 4.4. For
	all
				n ,Xn . (4.26)
	The global error is dened as
	e :	N ņ1	Rpu hτ , p hτ , αq X I
				» t F
		sup pv,qqX, pv,qq X 1	0

n

  Corollary 4.6 (Global in time a posteriori error estimate). Let pu, pq Y be the weak solution of (4.9), let Assumption 4.1 hold and let pu hτ , p hτ q be the discrete solution of(4.11). Let e ¤ η sp,mec η tm,mec η lin,mec η sp,hyd η tm,hyd .

	Applying the triangle inequality to distinguish the error sources, together with (4.28), yields
	e n,i ¤	¤ ¥	N ņ1	»	h In Ţ T n	¡	η n,i sp,mec,T	η n,i tm,mec,T ptq η n,i lin,mec,T	© 2 ¡	η n,i sp,hyd,T	η n,i tm,hyd,T ptq	© 2	1 {2
	¤ η n,i sp,mec		η n,i tm,mec	η n,i lin,mec		η n,i sp,hyd	η n,i tm,hyd .
	be dened by (4.27a) and the global error estimators given by (4.28) with (4.27). Then the
	following holds:
														(4.30)
	The proof of Corollary 4.6 is similar to the proof of Corollary 2.16.
														n,i be
	dened by (4.26) and the local error estimators given by (4.27). Then the following holds:
									e n,i ¤ η n,i sp,mec	η n,i tm,mec	η n,i lin,mec	η n,i sp,hyd	η n,i tm,hyd .	(4.29)
	Proof. Proceeding as in the proof of Theorem 2.13, using σ n,i h,disc	σ n,i h,lin instead of θ n h and
	taking into account (4.22) and (4.17), we obtain
	e n,i ¤	¤ ¥	N ņ1	»

e

  n,int h , and the set of all faces in ω a not lying 4.4. EQUILIBRATED FLUX A POSTERIORI ERROR ESTIMATE on the boundary fΩ if a V n,ext h . Using these spaces we can reformulate the local problems of Constructions (4.20) and(4.21). Note that the left hand sides of these two systems are equivalent. We introduce the functions Γ a and γ a , which for the discrete stress reconstructions are

  T , we can eliminate the variable Σ byΣ A ¡1 p¡B t Y ¡ C t X F q, ¡BA ¡1 C t ¡CA ¡1 B t ¡CA ¡1 C t

	CHAPTER 4. PORO-MECHANICAL PROBLEMS IN 3D	93
	leading to the system
	£ ¡BA ¡1 B t £ X Y	£ ¡BA ¡1 F G ¡CA ¡1 F	.
	Our next step is to eliminate Y , writing
			(4.35)

  α ¤ α ult , σ y aα ult if α ¡ α ult ,

				For the DruckerPrager yield criterion we can for
	example consider linear softening
		Rpαq	7 6 8	σ y aα
	or parabolic softening expressed by
		6 8
	Rpαq	7

if

  We consider the extension to three dimensions of the analytical solution of the linear elastic Biot problem of Section 2.5.2. Let Ω p0, 1q ¢ p0, 1q ¢ p0, 1q and t p0, 1 2 q and

	10 0
	10 ¡1	η sp,mec e mec η sp,hyd e hyd
		h ¡1	7 10 0.8
	pigure RFU ! gomprison of the spe disretiztion error estimte nd the error in the nlytil test se
	4.7.1 Analytical test

upt, x, y, zq sinp¡πtq

  The computation time of the static algorithms are measured without estimating the error. Since we performed each simulation only once, the given computation times should only be considered as an orientation, since they can vary for identic simulations. Nevertheless, we observe that the estimation of

		10 1.8	static				η sp , static
			adaptive				η tm , static
	Error estimation	10 1.2 10 1.4 10 1.6			10 1 10 1.5	η sp , adaptive η tm , adaptive
		10 1
			10 4.5	10 5	10 5.5	10 4		10 5
		Total number of space-time unknowns	Total number of space-time unknowns
		test	static iterations	static comp. time	static η	adaptive iterations	adaptive comp. time	adaptive η
		2	122	1min 9s	61.86	66		15min 38s	64.76
		3	170	1min 17s	55.30	31		4min 20s	60.03
		4	61	1min 6s	22.97	39		15min 38s	23.14

le RFP ! gomprison of the numer of itertionsD the omputtion time nd the glol estimte of the stti lgorithms nd lgorithms with the sme mesh nd time step using the dptive stopping riterion @RFQUA with γ lin 0.1 for the hruker!rger model

  q | τ h n ωa 0 on fω a if a V int h , τ h n ωa 0 on fω a zfΩ if a V ext h u,

	(A.17a)
	V a
	(A.17c)
	Construction A.4 (AFW stress reconstruction). Find σ a h Σ a h , r a h V a h and λ a h Λ a h such

. Let k : p and set Σ a h : tτ h S afw h pω a h : tv h V afw h pω a q | pv h , zq ωa 0 dz RM if a V int h u, (A.17b) Λ a h : Λ h pω a q.

  ΛT q with k 2 @rightA 116 BIBLIOGRAPHY IU tF fessonD qF gilletudD tF vF ghoheD nd F porestF Mécanique non linéaire des matériauxF riwi iene iurope vtdD PHHIF IV wF fietermn nd sF fu²kF he (nite element method for proli equtionsF Numer. Math.D RH@QAXQUQ!RHTD IWVPF IW wF eF fiotF qenerl theory of threeEdimensionl onsolidtionF J. Appl. Phys.D IPXISS!ITWD IWRIF PH hF fo0D wF fottiD nd hF eF hi ietroF e nononforming highEorder method for the fiot prolem on generl meshesD PHITF PI hF fo0D pF frezziD nd wF portinF edued symmetry elements in liner elstiityF Commun. Pur. Appl. Anal.D VXWS!IPID PHHWF PP hF fo0D pF frezziD nd wF portinF Mixed Finite Element Methods and ApplicationsD volume RR of Computational MathematicsF pringerD PHIQF PQ iF fonnetierF Mathematical Treatment of the Uncertainties Appearing in the Formulation of Some Models of PlasticityF hh thesisD he niversity of wrylnd t gollege rkD IWVVF PR fF foroomnd nd yF gF ienkiewizF eovery proedures in error estimtion nd dptivityF ssF edptivity in nonliner prolems of elstoEplstiity ehviourF Comput. Methods. Appl. Mech. Engrg.D IUTXIPU!IRTD IWWWF PS wF fottiD hF eF hi ietroD nd F ohlF e hyrid highEorder method for nonliner elstiityF rivXIUHUFHPISRD sumitted for pulitionD PHIUF PT F foussett nd vF pourmentF e posteriori error estimtion nd threeEdimensionl dptive remeshingX pplition to error ontrol of nonEstedy metl forming simultionsF sn International conference on numerical methods in industrial forming processesF yhio tteD xwspyw HRD IWWWF PU hF fressD F illweinD nd tF höerlF iquilirted residul error estimtes re pEroustF Comput. Methods Appl. Mech. Engrg.D IWVXIIVW!IIWUD PHHWF PV hF fress nd tF höerlF iquilirted residul error estimtor for edge elementsF Math. Comp.D UU@PTPAXTSI!TUPD PHHVF PW pF frezziD tF houglsD nd vF hF wriniF eent results on mixed (nite element methods for seond order ellipti prolemsF sn horodnitsyn flkrishnn nd vions idsFD editorsD Vistas in applied mathematics. Numerical analysis, atmospheric sciences, immunologyD pges PS!RQF yptimiztion oftwre snFD ulitions hivisionD xew orkD IWVTF QH F frink nd iF teinF e posteriori error estimion in lrgeEstrin elstiity using equilirted lol xeumnn prolemsF Comput. Methods Appl. Mech. Eng.D ITIXUU!IHID IWWVF QI gF grstensenF e unifying theory of posteriori (nite element error ontrolF Numer. Math.D IHH@SAXITQ! IUSD PHHSF QP gF grstensenD wF iigelD F rF F roppeD nd gF vöhrdF e review of uni(ed posteriori (nite element error ontrolF Numer. Math. Theory Methods Appl.D S@RAXSHW!SSVD PHIPF QQ wF germkD pF rehtD F ngD nd wF ohrlikF edptive inext itertive lgorithms sed on polynomilEdegreeEroust posteriori estimtes for the tokes prolemF rev reprint HIHWUTTPD suE mitted for pulitionD PHIUF QR wF gerverD wF ghiumentiD nd F godinF wixed stilized (nite element methods in nonliner solid mehnisX rt ssX trin loliztionF Comput. Methods in Appl. Mech. and Engrg.D IWW@QU!RHAXPSUI! PSVWD PHIHF QS eF vF ghillou nd wF uriF gomputle error estimtors for the pproximtion of nonliner prolems y linerized modelsF Comput. Meth. Appl. Mech. Engrg.D IWTXPIH!PPRD PHHTF QT vF ghmoinD F vdevèzeD nd pF ledF en enhned method with lol energy minimiztion for the roust posteriori onstrution of equilirted stress (eld in (nite element nlysisF Comput. Mech.D RWXQSU!QUVD PHIPF QU yF goussyF PoromechanicsF tohn iley nd onsD PHHRF QV F hestuynder nd fF wétivetF ixpliit error ounds in onforming (nite element methodF Math. Comput.D TV@PPVAXIQUW!IQWTD IWWWF QW hF eF hi ietro nd eF irnF e hyrid highEorder lokingEfree method for liner elstiity on generl meshesF Comput. Meth. Appl. Mech. EngrgD PVQXI!PID PHISF RH hF eF hi ietroD eF irnD nd tFEvF quermondF hisontinuous qlerkin methods for nisotropi semiE de(nite di'usion with dvetionF SIAM J. Numer. Anal.D RT@PAXVHS!VQID PHHVF RI hF eF hi ietroD iF plurudD wF ohrlíkD nd F ousefF e posteriori error estimtesD stopping riteriD nd dptivity for multiphse ompositionl hry )ows in porous mediF J. Comput. Phys.D PUTXITQ! IVUD PHIRF RP hF eF hi ietroD wF ohrlíkD nd F ousefF en posterioriEsedD fully dptive lgorithm for therml multiphse ompositionl )ows in porous medi with dptive mesh re(nementF Comput. and Math. with Appl.D TV@IPAXPQQI!PQRUD PHIRF RQ hF eF hi ietroD wF ohrlíkD nd F ousefF edptive regulriztionD lineriztionD nd disretiztion nd posteriori error ontrol for the twoEphse tefn prolemF Math. CompD VR@PWIAXISQ!IVTD PHISF RR godeester doument SFHQFHIF elgorithme non linéire qusiEsttiqueF www.code-aster.org/doc/ default/en/man_r/r5/r5.03.01.pdfF RS godeester doument UFHIFITF sntégrtion des omportements méniques élstoEplstiques de hruker!rgerD ssoié@drukprgerA et nonEssoié @drukprgnA et postEtritementsF www. code-aster.org/doc/default/en/man_r/r7/r7.01.16.pdfF RT godeester doument UFHIFPRF voi de omportement visoplstique viuF www.code-aster.org/ doc/default/en/man_r/r7/r7.01.24.pdfF RU F holej²íD eF irnD nd wF ohrlíkF hp!dption driven y polynomilEdegreeEroust posteriori error estimtes for ellipti prolemsF SIAM J. Sci. Comput.D QV@SAXeQPPH!eQPRTD PHITF RV F hörsek nd tF welenkF ymmetryEfreeD pEroust equilirted error indition for the hpEversion of the piw in nerly inompressile liner elstiityF Comput. Methods Appl. Math.D IQXPWI!QHRD PHIQF RW hFgF hruker nd F rgerF oil mehnis nd plsti nlysis or limit designF Quart. Appl. Math.D IH@PAXISU!ITSD IWSPF SH vF il elouiD eF irnD nd wF ohrlíkF qurnteed nd roust posteriori error estimtes nd lning disretiztion nd lineriztion errors for monotone nonliner prolemsF Comput. Methods Appl. Mech. Engrg.D PHHXPUVP!PUWSD PHIIF SI uF iriksson nd gF tohnsonF en dptive (nite element method for liner ellipti prolemsF Math. Comput.D SHXQTI!QVQD IWVVF SP eF irn nd F weunierF e posteriori error nlysis of iulerEqlerkin pproximtions to oupled elliptiE proli prolemsF ESAIM Math. Mod. Numer. Anal..D RQXQSQ!QUSD PHHWF SQ eF irnD sF mersD nd wF ohrlíkF iquilirted )ux posteriori error estimtes in L 2 pH 1 qEnorms for highEorder disretiztions of proli prolemsF preprint hlEHIRVWUPID PHIUF SR eF irn nd wF ohrlíkF e posteriori error estimtion sed on potentil nd )ux reonstrution for the het equtionF SIAM J. Numer. Anal.D RV@IAXIWV!PPQD PHIHF SS eF irn nd wF ohrlíkF edptive inext xewton methods with posteriori stopping riteri for nonliner di'usion hisF SIAM J. Sci. Comput.D QS@RAXeIUTI!eIUWID PHIQF ST eF irn nd wF ohrlíkF olynomilEdegreeEroust posteriori estimtes in uni(ed setting for onE formingD nononformingD disontinuous qlerkinD nd mixed disretiztionsF SIAM J. Numer. Anal.D SQ@PAXIHSV!IHVID PHISF SU eF irn nd wF ohrlíkF froken stle H 1 nd Hpdivq polynomil extensions for polynomilEdegreeE roust potentil nd )ux reonstrution in three spe dimensionsF preprint hlEHIRPPPHRD PHIUF SV F pernndesF Modélisation numérique objective des problèmes couplés hydromécaniques dans le cas des géomatériauxF hh thesisD niversité toseph pourier E qrenole sD PHHWF SW fF preijs de euekeF Displacement and equilibrium models in the nite element method. tohn iley nd onsD xew orkD IWTSF TH qF xF qtiD eF wrquezD nd F udolphF e priori nd posteriori error nlyses of ugmented twofold sddle point formultions for nonliner elstiity prolemsF Comput. Meth. Appl. Mech. Engrg.D PTR@IAXPQ!RVD PHIQF TI qF xF qti nd iF F tephnF e mixedEpiw formultion for nonliner inompressile elstiity in the plneF Numer. Methods Partial Dier. Equ.D IV@IAXIHS!IPVD PHHPF TP eF rnnukinenD F tenergD nd wF ohrlíkF e uni(ed frmework for posteriori error estimtion for the tokes prolemF Numer. Math.D IPPXUPS!UTWD PHIPF TQ iF roek nd iFF frownF impiril strength riterion for rok mssesF J. Geotech. Engng. Div.D IHTXIHIQ! IHQSD IWVHF TR iF roek nd iFF frownF Underground excavations in rockF vondonX snstnF winF wetllFD IWVHF TS iF roek nd iFF frownF he roek!frown filure riterion E IWVV updteF Proc. 15th Canadian Rock Mech. Symp.D pges QI!QVD IWVVF TT uFEF uimF qurnteed posteriori error estimtor for mixed (nite element methods of liner elstiity with wek stress symmetryF SIAM J. Numer. Anal.D RVXPQTR!PQVSD PHIIF TU uFEF uimF e posteriori error estimtor for liner elstiity sed on nonsymmetri stress tensor pE proximtionF J. Korean Soc. Ind. Appl. Math.D IT@IAXI!IQD PHIPF TV eF uleineF Modélisation numérique du comportement des ouvrages souterrains par une approche vis-coplastiqueF hh thesisD niversité de vorrineD PHHUF TW F vdevèzeF Comparaison de modèles de milieux continusF hh thesisD niversité ierre et wrie gurie @ris TAD IWUSF UH F vdevèzeF xouvelles proédures d9estimtion d9erreur reltive à l méthode des éléments (nis et pplitionsF Journées éléments nisD IWUUF UI F vdevèze nd hF veguillonF irror estimte proedure in the (nite element method nd pplitionsF SIAM J. Numer. Anal.D PHXRVS!SHWD IWVQF UP F vdevèze nd xF woësF e new posteriori error estimtion for nonliner timeEdependent (nite element nlysisF Comput. Meth. Appl. Mech. Engrg.D ISU@I!PAXRS!TVD IWWVF UQ F vdevèzeD tF F elleD nd F ougeotF irror estimtion nd mesh optimiztion for lssil (nite elementsF Engrg. Comp.D V@IAXTW!VHD IWWIF UR eF pF hF voul nd tF xF gF querreiroF pinite element nlysis of nonliner reeping )owsF Comput. Meth. Appl. Mech. Engrg.D UW@IAXVU!IHWD IWWHF US F vue nd fF sF ohlmuthF e lol posteriori error estimtor sed on equilirted )uxesF SIAM J. Numer. Anal.D RPXIQWR!IRIRD PHHRF UT F weunierF Analyse d'erreur a posteriori pour les couplages hydro-Mécaniques et mise en ÷uvre dans Code_AsterF hh thesisD Éole des onts risehD PHHUF WU hF ndriF ur l9pproximtion numérique des éoulements qusiExewtoniens dont l visosité suit l loi puissne ou l loi de grreuF Math. Modelling and Num. Anal.D PU@PAXIQI!ISSD IWWQF WV F iF howlterF hi'usion in poroEelsti mediF J. Math. Anal. Appl.D PSIXQIH!QRHD PHHHF WW iF tein nd F yhnimusF enisotropi disretiztionE nd modelEerror estimtion in solid mehnis y lol xeumnn prolemsF Comput. Meth. Appl. Mech. Engrg.D IUTXQTQ!QVSD IWWWF IHH iF teinD wF üterD nd F yhnimusF edptive (nite element nlysis nd modelling of solids nd struturesF pindingsD prolems nd trendsF Int. J. Numer. Methods Engrg.D THXIHQ!IQVD PHHRF IHI iF teinD wF üterD nd F yhnimusF irrorEontrolled dptive golEoriented modeling nd (nite element pproximtions in elstiityF Comput. Meth. Appl. Mech. Engrg.D IWTXQSWV!QTIQD PHHUF IHP gF ylor nd F roodF e numeril solution of the xvierEtokes equtions using the (nite element tehniqueF Comput. & FluidsD IXUQ!IHHD IWUQF IHQ vF F qF relorF The Physics of Rubber ElasticityF yxford niversity ressD eD IWUSF IHR F erfürthF A review of a posteriori error estimation and adaptive mesh-renement techniques. IWWTF IHS F erfürthF e review of posteriori error estimtion tehniques for elstiity prolemsF Comput. Meth. Appl. Mech. Engrg.D IUTXRIW!RRHD IWWWF IHT F erfürthF e posteriori error estimtes for (nite element disretiztions of the het equtionF CalcoloD RH@QAXIWS!PIPD PHHQF IHU wF ogeliusF en nlysis of the pEversion of the (nite element method for nerly inompressile mterilsF uniformly vlidD optiml error estimtesF Numer. Math.D RIXQW!SQD IWVQF IHV wF ohrlíkF yn the disrete oinré!priedrihs inequlities for nononforming pproximtions of the oolev spe H 1 F Numer. Funct. Anal. Opim.D PT@U!VAXWPS!WSPD PHHSF IHW wF ohrlíkF ni(ed priml formultionEsed priori nd posteriori error nlysis of mixed (nite element methodsF Math. Comp.D UWXPHHI!PHQPD PHIHF IIH uF von erzghiF Theoretical soil mechanicsF ileyD xew orkD IWRQF III eF ení²ekF he existene nd uniqueness theorem in fiot9s onsolidtion theoryF Aplikace MatematikyD PWXIWR!PIID IWVRF IIP wF illimsF yn the stress distriution t the se of sttionry rkF J. Appl. Mech.D PRXIHW!IIRD IWSUF IIQ F okooD uF mgtD nd rF xgokF pinite element method pplied to fiot9s onsolidtion theoryF Soils and FoundationsD IIXPW!RTD IWUIF IIR F ousefF A posteriori error estimates and adaptivity based on stopping criteria and adaptive mesh renement for multiphase and thermal ows. Application to steam-assisted gravity drainageF hh thesisD niversité ris TD PHIQF IIS yF gF ienkiewizD F gF viuD nd qF gF rungF irror estimtion nd dptivity in )ow formultion for forming prolemsF Internat. J. Numer. Methods Engrg.D PS@IAXPQ!RPD IWVVF IIT yF gF ienkiewiz nd tF F huF e simple error estimtor nd dptive proedure for prtil engineering nlysisF Internat. J. Numer. Methods Engrg.D QQ@PAXQQU!QSUD IWVUF IIU yF gF ienkiewiz nd tF F huF he superonvergent pth reovery nd posteriori error estimtesF prt IX he reovery tehniqueF Internat. J. Numer. Methods Engrg.D QQXIQQI!IQTRD IWWPF IIV yF gF ienkiewiz nd tF F huF he superonvergent pth reovery nd posteriori error estimtesF prt PX irror estimtes nd dptivityF Internat. J. Numer. Methods Engrg.D QQ@UAXIQTS!IQVPD IWWPF

T , V afw T ,

1.5. PROBLÈME PORO-MÉCANIQUE EN 3D

2.3. QUASI-STATIC FLUX RECONSTRUCTIONS

2.5. NUMERICAL RESULTS

http://web-code-aster.org

Remerciements

Abstract

We present an a posteriori error estimate for the linear elasticity problem. The estimate is based on an equilibrated reconstruction of the Cauchy stress tensor, which is obtained from mixed nite element solutions of local Neumann problems. We propose two dierent reconstructions: one using ArnoldWinther mixed nite element spaces providing a symmetric stress tensor, and one using ArnoldFalkWinther mixed nite element spaces with a weak symmetry constraint. The performance of the estimate is illustrated on a numerical test with analytical solution. 109 posed of L 2 pΩq functions with weak gradients in L 2 pΩq and H 1 0 pΩq for its zero-trace subspace. The weak formulation of problem (A.1) reads: nd u H 1 0 pΩq such that pσpuq, εpvqq pf, vq dv H 1 0 pΩq.

(A.

2)

The discretization of (A.2) is based on a conforming triangulation T h of Ω, verifying the minimum angle condition. We will use a conforming nite element method of order p ¥ 2. Let P p pT h q : tv L 2 pΩq | dT T h v| T P p pTqu, where P p pTq is the space of polynomials on T of degree less than or equal to p. For the sake of simplicity we assume that f lies in P p¡1 pT h q.

Then the discrete problem reads: nd u h H 1 0 pΩq P p pT h q such that pσpu h q, εpv h qq pf, v h q dv h H 1 0 pΩq P p pT h q.

(A.3)

A.3 A Posteriori Error Estimate

In this section, we derive an upper bound on the error between the analytical solution of (A.2)

and an arbitrary function u h H 1 0 pΩqP p pT h q. We will measure this error in the energy norm

where the last bound follows from λ ¥ 0 and Korn's inequality. Owing to (A.1b), we have C K 1 2 (this value would have been dierent if we had chosen mixed boundary conditions). We start by introducing reconstructed stress tensors that are more physical than σpu h q, which in general does not lie in Hpdiv, Ωq tτ L 2 pΩq | ∇ ¤ τ L 2 pΩqu and thus cannot verify the equilibrium equation (A.1a). Unlike σpu h q, however, these reconstructed tensors may not be symmetric.

Denition A.1 (Equilibrated stress reconstruction). We call equilibrated stress reconstruction any function σ h Hpdiv, Ωq constructed from σpu h q such that p¡∇ ¤ σ h , zq T pf, zq T dz RM dT T h ,

where RM :

Ru is the space of rigid body motions.

Theorem A.2 (A posteriori error estimate). Let u H 1 0 pΩq solve (A.2) and u h H 1 0 pΩq be arbitrary. Let σ h be a stress reconstruction verifying Denition A.1. Then

Using the denitions (A.17), the formulation (A.18) is equivalent to a modied version of (A.9), adding the weak symmetry constraint (A.18c). The condition (A.15) for all a V int h ensures that the constrained minimization problem (A.18) is well-posed.

A.4.3 Properties of the Stress Reconstructions

For both stress reconstructions we obtain the following result, recalling that we assume f to be piecewise polynomial of degree p ¡ 1.

Lemma A.5 (Properties of σ h ). Let σ h be prescribed by Construction A.3 or Construction A.4.

Then σ h Hpdiv, Ωq, and for all T T h , the following holds:

Proof. All the elds σ a h are in Hpdiv, ω a q and satisfy appropriate zero normal conditions so that their zero-extension to Ω is in Hpdiv, Ωq. Hence, σ h Hpdiv, Ωq. Let us prove (A.19). Since (A.15) holds for all a V int h , we infer that (A.14b) or (A.18b) is actually true for all v h V h pω a q. The same holds if a V ext h by denition of V a h . Hence, pψ a f ∇ ¤ σ a h , v h q ωa 0 for all v h V h pω a q and all a V h . Since V h pω a q is composed of piecewise polynomials that can be chosen independently in each cell T T a , and using σ h | T °aV T σ a h | T and the partition of unity °aV T ψ a 1, we infer that pf ∇ ¤ σ h , vq 0 for all v V T and all T T h . The fact that pf ∇ ¤ σ h q| T V T for any T T h , concludes the proof.

A.5 Numerical Results

We illustrate numerically our theoretical results on a test case with a known analytical solution.

We analyze the convergence rates of the error estimates and compare them to those of the analytical error. The computations were performed using the Code_Aster 1 software. The exact solution u pu x , u y q on the unit square Ω p0, 1q 2 is given by

with the Lamé parameters µ λ 1, and the corresponding body force f. The exact solution is imposed as Dirichlet condition on the whole boundary fΩ. The discretization is done on a series of unstructured grids with the polynomial degree p 2 in the conforming nite element method (A.3). For each computation, two error estimates are calculated, one for each stress reconstruction. The AFW reconstruction oers some advantages over the AW one:

it is cheaper (since by hybridization techniques we can avoid the resolution of saddle point problems), and the implementation for three-dimensional problems is easier (the lowest-order AW element in 3D has 162 degrees of freedom per element).

1 http://web-code-aster.org