Renée El

Melhem Bat

Blaise Pascal

Insa : R Gourdon

M Stéphane

M Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

S Tordeux

Méthodes 2004

CHAPITRE 1. PRÉSENTATION DE LA THÈSE CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA

Keywords: Analyse multi-échelle Asymptotic analysis, multi-scale analysis, finite element method, inverted finite element method, linear elasticity, nonlinear elasticity

de ces années.

Table des matières

Table des matières vii

Liste des figures ix

Liste des tableaux xi

Liste des figures

1.1 Domaine initial Ω = Ω ε f ∪ Γ ε ∪ Ω ε m with Ω ε f =

Introduction

L'agencement de composants de géométries, propriétés et dimensions différentes est au coeur de tous les enjeux et requiert une importance capitale depuis que l'homme a commencé à modeler son environnement pour en faire une plateforme adaptée à ses besoins. L'évolution de nos besoins nous a alors conduit à devoir mener des études de plus en plus approfondies de l'intéraction entre ces composants surtout quand le rapport de leurs dimensions ou celui de leurs propriétés est très élevé. L'enjeu majeur qui motive ce genre d'études est l'apparition de contraintes souvent extrêmement élevées autour de l'interface de transition entre deux corps. Nous pouvons citer à titre d'exemple les cordons de soudure, les ailettes de refroidissement, les rivets dans les coques d'avions ou encore les renforts dans la gomme d'un pneumatique (notre cas d'intérêt). Dans le cas où les deux corps sont de dimensions caractéristiques très différentes D 1 D 2 → 0, l'approximation numérique des champs en question s'avère inefficace. Ainsi, à moins de mailler d'une façon excessivement fine autour des zones de petites dimensions, ce qui devient rapidement exhorbitant en termes de moyens et de temps de calcul, ou d'avoir recours à des techniques d'homogéneisation classiques ce qui requiert la présence d'une densité surfacique ou volumique de ces composants, l'approximation obtenue est de mauvaise qualité et n'est pas représentative de la réalité.

Toutes ces raisons nous amènent à vouloir proposer une modélisation mathématique la plus exacte possible et un shéma numérique robuste et efficace de l'intéraction entre deux milieux élastiques de rigidités différentes et dont le rapport des dimensions est très petit. L'objectif final est de proposer une stratégie numérique itérative de prise en compte de plusieurs inclusions dans un substrat élastique.

Contexte industriel

CHAPITRE 1. PRÉSENTATION DE LA THÈSE de gagner en précision, la présence de raffinements locaux peut être nécessaire, et l'efficacité du calcul pourra nécessiter une méthode itérative entre comportement global et comportement local au voisinage des inclusions.

Contexte scientifique

L'enjeu scientifique consiste en premier lieu à proposer une modélisation de l'influence d'une inclusion de dimension caractéristique ε petite par rapport à la dimension R du domaine dans lequel elle évolue. L'objectif est d'effectuer une analyse autour d'une seule inclusion c'est à dire sans passer par des techniques d'homogéneisation classiques qui requièrent une densité périodique d'inclusions, bien que ces dernières soient efficaces dans certains cas. Ainsi, donner une information précise sur l'influence d'une inclusion sur le substrat en termes de déplacements et surtout de contraintes donnera au concepteur un outil puissant lui permettant d'optimiser leur agencement, nombre et forme ce qui lui permettera de ne plus surdimensionner et impactera d'une façon significative le coût de production surtout quand il s'agit de specimens tels que les pneumatiques d'avions ou d'engins de génie civil. Cette modélisation permettera en second lieu de construire une stratégie numérique permettant d'approcher le problème d'une façon fiable et précise sans avoir à mailler d'une façon excessivement fine.

Etat de l'art

Certains travaux ont traité le cas des perturbations que ce soit pour les inclusions, pour des problèmes de transmission essentiellement pour la conductivité SÁNCHEZ-PALENCIA [1974] et l'équation de Helmholtz BENDALI et LEMRABET [1996]; HANSEN et collab. [2007] ou même pour les cavités et les domaines perforés. Le cas des domaines perforés a été étudié dans DAMBRINE et VIAL [2005,2007]; IL 'IN [1992]; LEWINSKI et SOKOLOWSKI [2000]; MAZ'YA et collab. [2000a,b] en utilisant des techniques s'appuyant sur la notion de profil, une solution normalisée de l'équation de Laplace dans le domaine extérieur est obtenue par explosion de la perturbation. En effet, l'influence d'une perforation sur le bord d'un domaine lisse est étudiée dans DAMBRINE et VIAL [2005] et un développement asymptotique de la fonction de forme en fontion d'un paramètre de taille en est dérivé.

Le cas des inclusions a été traité pour la première fois dans ESHELBY [1957] donnant une méthode pour trouver les contraintes obtenues dans un solide élastique lorsque l'inclusion subit un changement de forme (Pour une forme ellipsoïdale dont les constantes élastiques diffèrent de celles du reste du domaine). Le résultat le plus important est que le champ dans l'inclusion peut être déterminé sans avoir à chercher le champ en dehors de celle-ci. Le champ élastique résultant est obtenu à l'aide d'une séquence d'opérateurs imaginaires de coupe, de contrainte et de couplage et utilisant un potentiel biharmonique. La détermination du champ en dehors de l'inclusion a été brièvement traitée dans cet article mais puis a été approfondie dans ESHELBY [1959] où le champ élastique en dehors d'une inclusion ellipsoïdale s'exprime entièrement en termes du potentiel harmonique d'une ellipsoïde solide en remplaçant le potentiel biharmonique introduit dans ESHELBY [1957] par le potentiel harmonique d'une certaine distribution surfacique. Il est à noter que la solution proposée par Eshelby pour l'étude des inclusions est l'une des CHAPITRE 1. PRÉSENTATION DE LA THÈSE contributions les plus importantes à la micro-mécanique au vingtième siècle et repose sur l'hypothèse qu'une inclusion est intégrée dans un espace ambiant non borné. Schiavone a aussi abordé dans SCHIAVONE [2003] le cas d'une inhomogénéité elliptique incluse dans un domaine infini élastique et isotrope soumis à des déformations de cisaillement sous l'hypothèse d'un chargement non uniforme et a ensuite étendu l'étude avec Antipov dans ANTIPOV et SCHIAVONE [2003] au cas d'une inhomogénéité sous les mêmes hypothèses mais de forme arbitraire et d'interface imparfaite cette fois-ci. Une méthode y est présentée conduisant à la solution du problème en question avec l'identification de la forme de l'inhomogénéité et la forme de la fonction d'interface correspondante qui conduit à un champ de contraintes uniforme à l'intérieur. L'analyse est basée sur des méthodes à variables complexes. Li, Sauer et Wang ont approfondi dans LI et collab. [2005] la caractérisation précise des champs élastiques dûs à la présence d' inclusions dans des supports élastiques finis. Il est important de dire que, avant le travail effectué dans LI et collab. [2005], Kroner KRÖNER [1986,1990] et Mazilu MAZILU [1972] ont essayé d'étudier le problème d'inclusion en domaine fini en cherchant une fonction de Green de l'équation de Navier. La tentative s'est terminée par un échec en raison de la difficulté mathématique impliquée dans l'obtention d'une solution fermée de la fonction de Green en domaine fini. D'un autre côté, Ammari et Kang ont décrit dans AMMARI et KANG [2004] des techniques pour la reconstruction de petites inclusions en utilisant uniquement des mesures au bord. Ces techniques reposent sur des développements asymptotiques des perturbations sur le bord dues à la présence d'inclusions. L'approche utilisée pour les développements asymptotiques est basées sur des techniques de potentiel de couche permettant de traiter des bords rugueux. Par l'intermédiaire de la méthode de la variable complexe de Muskhelishvili et la technique de transformation conforme de Schwartz-Christoffel et motivés par les contraintes élevées qui se produisent à la jonction entre les fibres et la matrice dans des feuilles composites multi-fibres qui peuvent provoquer une rupture, CHAND et CONWAY [1968] a présenté une solution pour déterminer les contraintes dans une plaque élastique infinie contenant une inclusion rigide rectangulaire soumise à un champ de contraintes uniforme. En utilisant également la méthode de la fonction à variable complexe, le cas d'une inclusion de forme arbitraire intégrée dans un demi-espace est étudié dans SUN et PENG [2003]. Les solutions sont obtenues par le biais de techniques d'extension analytique, de continuation analytique et de cartographie conforme. Dans BERETTA et collab.

[2012] un développement asymptotique de la différence entre le déplacement à la frontière avec et sans inclusion est présenté (le corps et l'inclusion peuvent être anisotropes). Pour dériver le développement asymptotique, les auteurs suivent l'approche introduite par Capdeboscq et Vogelius dans CAPDEBOSCQ et VOGELIUS [2003] (voir aussi CAPDE-BOSCQ et VOGELIUS [2004]) pour le problème de conductivité. Notons que la formule obtenue généralise ceux déjà obtenues dans le cas de corps homogènes isotropes à faible inertie de diamètres faibles.

Pour notre problème d'inclusion, nous allons nous inspirer des techniques utilisées pour des domaines perforés introduite précedemment où un développement asymtotique en termes d'un paramètre de taille et de correcteurs à différents ordres est dérivé. Cette méthode des correcteurs a une alternative qui est celle des développements asymptotiques raccordés et qui, bien que plus difficile à mettre en oeuvre, offre certains avantages. Développée et présentée dans BENDALI et collab. [2012,2016,2013,2011] , elle consiste à proposer un développement asymptotique dit intérieur au voisinage de l'in-CHAPITRE 1. PRÉSENTATION DE LA THÈSE homogénéité et un deuxième dit extérieur assez loin du domaine de celle-ci et de raccorder les deux d'une façon adéquate. La construction des deux développements se base sur l'introduction d'un champ global, dit extérieur, décrivant le comportement global et tronqué au voisnage immédiat de l'inhomogénéité et un champ défini à l'intérieur d'une couche limite entourant celle-ci. Les conditions de raccord lient les développements intérieur et extérieur et permettent leur détermination complète. Ces conditions expriment ce qui est connu sous le nom du principe de Van Dyke imposant que ces deux développements correspondent à la même solution. Les expressions spécifiques des conditions de raccord dépendent, bien évidemment, des problèmes considérés. Un des avantages de la méthode raccordée est qu'elle met en évidence les singularités au voisinage du centre des hétérogénéités en variable lente (singularité par exemple en 1 r ) alors que la méthode des correcteurs les laisse implicites dans le comportement en variable rapide (singularité par exemple en 1 R ). Cependant la problématique industrielle qui consiste à déterminer l'influence de l'inclusion sur le domaine dans lequel elle évolue et ne porte pas grand intérêt aux champs à l'intérieur de l'inclusion oriente notre choix vers la méthode des correcteurs qui est plus facile à mettre en oeuvre.

Formulations mathématiques

Toutes les fonctions et distributions sont définies sur l'espace Euclidien réel bidimensionnel R 2 .

Soit

r = |x| = (x 2 1 + x 2 2 )
1 2 la distance entre le point x = (x 1 , x 2 ) et l'origine. Rappelons que D(R 2 ) désigne l'espace des fonctions indéfiniment différentiables à support compact et D (R 2 ) son dual associé appellé espace des distributions. Avec λ = (λ 1 , . . . ,

λ n ) ∈ N n un multi-indice D λ = D λ 1 1 • • • D λ n n est l'opératuer différentiel d'ordre |λ| = λ 1 + • • • + λ n . L 2 (R 2 ) est l'espace classique des fonctions mesurables tel que ( R 2 |u| 2 dx) < ∞. C'est un espace de Banach pour la norme u = ( R 2 |u| 2 dx) 1 2 . Avec ρ = (1 + r 2 )
1 2 l'espace de Sobolev à poids, adapté à notre cas d'étude, est introduit par :

W 1,2 0,0 (R 2 ) = u ∈ D (R 2 ), ρ -1 (lg r ) -1 u ∈ L 2 (R 2 ), ∇u ∈ (L 2 (R 2
)) 2 , qui est un espace de Banach reflexif équipé de sa norme naturelle :

u W 1,2 0,0 (R 2 ) = ( ρ -1 (lg r ) -1 u 2 L 2 (R 2 ) + ∇u 2 L 2 (R 2 ) ) 1 2 .
Nous définissons aussi la semi norme :

|u| W 1,2 0,0 (R 2 ) = ∇u L 2 (R 2 ) .
Nous posons lg(r ) = ln(2 + r 2 ), B R = B(0, R) la boule ouverte de centre 0 et de rayon R et B R = (B R ) c l'extérieur B R . Pour finir, nous définissons P 0 comme étant l'espace des fonctions constantes. Nous utiliserons aussi, certaines fois, la notation f g à la place de : il existe une constante C tel que f ≤ Cg .

Nous introduisons maintenant les deux problèmes de transmission auquels nous allons nous intéresser : une approximation scalaire anti-plane et une approximation linéaire en déformation plane. Nous commençons par définir le domaine d'intérêt. Soit f est aussi considérée régulière. Le problème anti-plan réduit le problème mécanique à un problème de Poisson et est caractérisé par un déplacement réduit à sa composante normale au plan (O, x 1 , x 2 ) sous la forme u 3 = u(x 1 , x 2 ). Les conditions de transmission qui sont des conditions de continuité sont imposées sur l'interface Γ ε et une condition de Dirichlet non homogène est posée sur le bord Γ.

La formulation forte du problème de transmission anti-plan s'écrit : Le second problème de transmission traité est un problème d'élasticité linéaire en déformation plane avec deux tenseurs de rigidité d'ordre quatre E f et E m dans les domaines de l'inclusion et de la matrice. La loi de comportement reliant le tenseur des contraintes à celui des déformations est :

α f ∆u ε f = 0 dans Ω ε f , (1.1a) α m ∆u ε m = 0 dans Ω ε m , (1.1b) u ε f = u ε m sur Γ ε , (1.1c) α f ∂ n u ε f = α m ∂ n u ε m sur Γ ε , (1.1d) u ε m = U sur Γ.
σ k = E k : k , k ∈ { f ,
m} (où ":" est le produit doublement contracté). Les conditions de transmission qui sont des conditions de continuité des déplacements et des contraintes sont imposées sur l'interface Γ ε et une condition de Dirichlet non homogène est imposée sur le bord Γ. La formulation forte du problème en élasticité linéaire s'écrit : .

div σ f (u ε f ) = 0 dans Ω ε f , (1.2a) div σ m (u ε m ) = 0 dans Ω ε m , (1.2b) u ε f = u ε m sur Γ ε , (1.2c) σ f (u ε f ) • n f = -σ m (u ε m ) • n m sur Γ ε , (1.2d) u ε m = U sur Γ. (1.2e) où u ε f : Ω ε f → R 2 est la solution du problème dans le domaine de l'inclusion Ω ε f , u ε m : Ω ε m → R 2 est
Γ ε 0 2ε Γ Ω ε m R Ω ε f . FIGURE 1.1 -Domaine initial Ω = Ω ε f ∪ Γ ε ∪ Ω ε m with Ω ε f = εΩ 1
f . Noter que Γ = ∂Ω n'est pas nécessairement de forme circulaire.

Résumé des travaux

Cette thèse s'insère dans le cadre d'une bonne approximation numérique de problèmes de transferts entre une inclusion et un substrat de dimension très supérieure à celle-ci. Les deux milieux sont considérés comme étant élastiques et la contribution de l'inclusion sera présentée sous forme de correcteurs à différents ordres indépendants de la taille caractéristique de l'inclusion ε. La thèse est organisée en trois chapitres. Le premier, celui en cours, est introductif présentant le cadre industriel et scientifique et les résultats importants obtenus. Le second est consacré à une modélisation mathématique de l'influence d'une inclusion fine sur un substrat élstique de rigidité différente pour un problème de Poisson et un problème en élasticité linéaire homogène et isotrope moyennant une analyse asymptotique multi-échelle. Cette première phase nous sert à préparer les outils nécessaires à une approximation numérique dans les normes. Le troisième chapitre s'intéresse à l'approximation numérique de cette influence par le biais de la méthode des éléments finis (FEM) pour certains correcteurs et la méthode des éléments finis inversés pour d'autres (IFEM). Ce chapitre comprend surtout une stratégie itérative de prise en compte de plusieurs inclusions. Le quatrième et dernier chapitre traitte d'une adaptation de la méthode à des lois de comportement non linéaires en commençant par des potentiels simples pour finir avec un potentiel quelconque. Toutes les simulations présentées ont été effectuées en utilisant la bibliothèque d'éléments finis GetFEM++ RENARD et POMMIER [2014]. Nous nous intéressons dans ce chapitre à la construction d'un développement asymptotique pour les deux problèmes présentés dans la Section 1.5 à un ordre arbitraire K sous la forme suivante (k ∈ { f , m}) :

u [K] ε,k (x) = u (0) (x) + K i =2 ε i v (i ) (x) + K i =1 ε i V (i -1) k ( x ε ), (1.3)
où u (0) est la solution du problème (1.1) quand ε = 0 (i.e sans inclusion), les fonctions v (i ) sont définies sur Ω et les V (i ) sont définies sur tout R 2 .

Un aspect important pour les applications est que les correcteurs v (i ) et V (i ) obtenus ne dépendent pas de ε et dépendent linéairement de u (0) par le développement de son gradient en x = 0 ce qui permet d'en calculer une base, une fois pour toutes, par forme d'inclusion.

Nous construisons le développement asymptotique étape par étape en commençant par introduire la différence

d ε 0,k = u ε k -u (0) pour k ∈ { f , m} qui résout le problème suivant : α f ∆d ε 0, f (x) = 0 dans Ω ε f , (1.4a) α m ∆d ε 0,m (x) = 0 dans Ω ε m , (1.4b) d ε 0, f (x) = d ε 0,m (x) sur Γ ε , (1.4c) α f ∂ n d ε 0, f (x) = α m ∂ n d ε 0,m (x) + (α m -α f )∂ n u (0) (x) sur Γ ε , (1.4d) d ε 0,m (x) = 0 sur Γ.
(1.4e)

Le développement est construit en compensant le terme qui apparait sur Γ ε à chaque coup et d'autres termes qu'on verra apparaitre sur le bord extérieur Γ pour les ordres suivants.

Pour cela nous transformons le problème initial en un problème posé sur une géométrie indépendante du petit paramètre ε. Nous figeons en premier lieu l'interface moyennant une remise à l'échelle 'scaling' y = x ε (cette technique nous a été inspirée par DAMBRINE et VIAL [2005]; VIAL [2003]

) et introduisons pour k ∈ { f , m} D ε 0,k (y) = d ε 0,k (x)
. Nous faisons ensuite tendre ε vers 0 afin d'éliminer toute dépendance en ε et obtenir le problème suivant :

α f ∆D ε 0, f (y) = 0 dans Ω 1 f , (1.5a) α m ∆D ε 0,m (y) = 0 dans Ω ∞ , (1.5b) D ε 0, f (y) = D ε 0,m (y) sur Γ 1 , (1.5c) α f ∂ n D ε 0, f (y) = α m ∂ n D ε 0,m (y) + ε(α m -α f )∂ n u (0) (εy) sur Γ 1 , (1.5d) D ε 0,m (y) → 0 quand |y| → ∞. (1.5e)
La figure 1.2 illustre la géométrie dilatée. Ceci nous permet d'introduire un problème posé sur la géométrie dilatée dont les solutions sont les fonctions V (p) intervenant dans le développement (1.3). Les régularités du domaine nous permettent d'écrire un développement de Taylor de ∇ ∇ ∇u (0) (εy) pour y ∈ Γ 1 :

∇ ∇ ∇u (0) (εy) = K i =0 ε i ω (i ) 0 (y) + O(ε K+1 ). (1.6) CHAPITRE 1. PRÉSENTATION DE LA THÈSE . 0 2 Ω 1 f Γ 1 Ω ∞ . FIGURE 1.2 -Géometrie dilatée.
Nous introduisons une série de problèmes posés sur la géométrie dilatée dont les solutions vont se charger de compenser le terme qui apparait sur le bord Γ 1 :

α f ∆V (p) f (y) = 0 dans Ω 1 f , (1.7a) α m ∆V (p) m (y) = 0 dans Ω ∞ , (1.7b) V (p) f (y) = V (p) m (y) sur Γ 1 , (1.7c)
α f ∂ n V (p) f (y) = α m ∂ n V (p) m (y) + ϕ (p) (y) sur Γ 1 , (1.7d) V (p) m (y) → 0 quand |y| → ∞.
(1.7e) L'approximation au premier ordre est construite en compensant le premier terme du développement de (α m -α f )∂ n u (0) (εy). Ainsi avec V (0) solution du Problème (1.7) pour ϕ (0) (y) = (α m -α f )ω (0) 0 (y) • n nous actualisons le résidu :

d ε 1,k (x) = u ε k (x) -u (0) (x) -εV (0) k ( x ε ), k ∈ { f , m}. (1.8)
Et nous en tirons l'approximation au premier ordre :

u [1] ε,k (x) = u (0) (x) + εV (0) k ( x ε ), k ∈ { f , m}.
(1.9)

Le problème résolu par d ε 1,k est le suivant :

α f ∆d ε 1, f (x) = 0 dans Ω ε f , (1.10a) α m ∆d ε 1,m (x) = 0 dans Ω ε m , (1.10b) d ε 1, f (x) = d ε 1,m (x) sur Γ ε , (1.10c) α f ∂ n d ε 1, f (x) = α m ∂ n d ε 1,m (x) + (α m -α f )(∂ n u (0) (x) -ω (0) 0 (y)) sur Γ ε , (1.10d) d ε 1,m (x) = -εV (0) m ( x ε ) sur Γ.
(1.10e)

CHAPITRE 1. PRÉSENTATION DE LA THÈSE
Nous introduisons un problème équivalent au problème résolu par V (p) dont le rôle est de compenser le terme qui est apparu sur Γ :

α m ∆ x v (q) (x) = 0 dans Ω, (1.11a) v (q) (x) = ψ (q) (x)
sur Γ.

(1.11b)

Nous prouvons que V (0) m se développe à l'infini sous la forme suivante :

V (0) m (y) = K i =1 a (i ) 0 (θ) |y| i + O 1 |y| K+1 = K i =1 ε i a (i ) 0 (θ) |x| i + O(ε K+1 ),
x ∈ Γ.

(1.12) où θ est la variable angulaire en coordonnées polaires. Maintenant, avec V (1) solution du problème pour

ϕ (1) (y) = (α m -α f )ω (1) 0 (y) • n = (α m - α f )ω (1) 0 (y) et v (2) solution du problème pour ψ (2) (x) = - a (1) 0 (θ)
|x| nous actualisons le résidu encore une fois pour en tirer le développement au second ordre :

u [2] ε,k (x) = u (0) (x) + εV (0) k ( x ε ) + ε 2 v (2) (x) + ε 2 V (1) k ( x ε ), k ∈ { f , m}. (1.13)
Nous allons jusqu'au quatrième ordre afin de déterminer les formes de ϕ (p) et ψ (p) et pouvoir construire le développement à un ordre arbitraire K annoncé dans (1.3). Le résidu à l'ordre K résout alors le problème suivant :

α f ∆d ε K, f (x) = 0 dans Ω ε f , (1.14a) α m ∆d ε K,m (x) = 0 dans Ω ε m , (1.14b) d ε K, f (x) = d ε K,m (x) sur Γ ε , (1.14c) α f ∂ n d ε K, f (x) = α m ∂ n d ε K,m (x) + O(ε K ) sur Γ ε , (1.14d) d ε K,m (x) = O H 1 2 (Γ) (ε K+1 )
sur Γ.

(1.14e)

Et nous établissons l'estimation d'erreur suivante :

Theorème 1.6.1. Il existe une constante C > 0 independente de ε tel que :

d ε K H 1 (Ω) ≤ Cε K , pour tout K ∈ N, où d ε K est le résidu à l'ordre K.
Afin de valider le modèle construit nous résolvons en premier lieu le Problème (1.1) par la méthode de séparation des variables puis nous résolvons le Problème résolu par le premier correcteur V (0) (1.7) par la même méthode et ceci pour une géométrie circulaire de l'inclusion. Nous constatons que la première contribution de l'inclusion est en ε 2 dans les deux cas ce qui représente une première validation. L'adaptation à l'élasticité linéaire est effectuée en introduisant des correcteurs analogues à ceux introduits pour le problème de Poisson à part le fait qu'il y a un terme supplémentaire à compenser sur le domaine Ω ε f -div x (E f : (u (0) ))(x) dû à la différence des rigidités. Pour ce terme là nous introduirons un troisième problème de compensation

,k = u ε k -u (0) est le suivant : div x σ f (d ε 0, f )(x) = -div x (E f : (u (0) ))(x) dans Ω ε f , (1.15a) div x σ m (d ε 0,m (x)) = 0 dans Ω ε m , (1.15b) d ε 0, f (x) = d ε 0,m (x) sur Γ ε , (1.15c) σ f (d ε 0, f (x)) • n f = -σ m (d ε 0,m (x)) • n m +[(E m -E f ) : (u (0) (x))] • n f sur Γ ε , (1.15d) d ε 0,m (x) = 0 sur Γ.
(1.15e)

Ceci nous amène au développement à l'ordre K suivant :

u [K] ε,k (x) = u (0) (x) + K i =2 ε i v (i ) (x) + K i =1 ε i V (i -1) k ( x ε ) + K i =1 ε i +1 Z (i -1) k ( x ε ).
(1.16)

Le résidu à l'ordre K résout le problème suivant :

div σ f (d ε K, f )(x) = O H -1 (Ω;R 2 ) (ε K+1 ) dans Ω ε f , (1.17a) div σ m (d ε K,m )(x) = 0 dans Ω ε m , (1.17b) d ε K, f (x) = d ε K,m (x) sur Γ ε , (1.17c) σ f (d ε K, f )(x) • n f = -σ m (d ε K,m )(x) • n m + O(ε K ) sur Γ ε , (1.17d) d ε K,m (x) = O H 1 2 (Γ,R 2 ) (ε K+1 )
sur Γ.

(1.17e)

Et nous établissons l'estimation d'erreur suivante :

Theorème 1.6.2. Dans le cas isotrope E m : = λ trace( )I+2µ , il existe une constante C > 0 indépendante de ε tel que :

d ε K H 1 (Ω;R 2 ) ≤ Cε K , pour tout K ∈ N.
Nous prouvons aussi que les fonctions

V (p) m et Z (p)
m se développent à l'infini comme tel :

K i =1 a (i ) 0 (θ) |y| i + O 1 |y| K+1 = K i =1 ε i a (i ) 0 (θ) |x| i + O(ε K+1 ),
x ∈ Γ.

(1.18)

Pour finir, nous établissons l'existence et l'unicité respectivement des fonctions V (p)

m et V (p) m , Z (p)
m posées sur des domaines infinis intervenant dans les développements (1.3) et(1.16) dans les espaces adéquats suivants :

V = v ∈ W 1,2 0,0 ; Γ R 0 vdσ = 0 et V = v ∈ (W 1,2 0,0 ) 2 ; Γ R 0 v dσ = 0 où 0 < R 0 < R est un réel
arbitraire et W 1,2 0,0 est l'espace de Sobolev à poids introduit précédemment.

Chapitre 3 : Une stratégie numérique de prise en compte de l'influence d'une inclusion fine sur un substrat élastique

Notre sujet d'intérêt dans ce chapitre est d'approximer numériquement l'influence d'une fine inclusion sur un domaine élastique pour les deux problèmes introduits précédemment.

CHAPITRE 1. PRÉSENTATION DE LA THÈSE

A moins de mailler d'une façon extrêmement fine en vue de sa petite dimension, ce que nous voulons éviter, la méthode des éléments finis classique demeure incapable de capter l'inclusion. Nous faisons alors appel au formalisme introduit dans la Section 1.6.1 et nous nous contentons ici d'approcher numériquement le développement asymptotique au premier ordre (1.9) et son analogue pour l'élasticité plane (la première contribution étant en ε 2 les autres ordres sont jugés négligeables pour une inclusion assez fine). Ce développement comporte deux termes : u (0) qui sera approché par la méthode des éléments finis classique (FEM) et V (0) qui sera approché par la méthode des éléments finis inversés (IFEM) que nous adapterons aux deux problèmes de transmission bidimensionnels résolus par les correcteurs posés sur la géométrie dilatée (tout R 2 ) pour le problème de Poisson et le problème en élasticité linéaire. Nous commençons par faire l'inventaire des méthodes d'approximation sur des domaines infinis et motivons le choix de la méthode des éléments finis inversés par le fait qu'elle soit générique et applicable à une multitude de problèmes, exactement conforme et n'utilisant aucune condition artificielle ou condition de bord. Ensuite le cadre de la méthode est présenté en commençant par introduire la notion de simplexe infini et son simplexe fini associé dans le cas bidimensionnel :

T ∞ (a 0 | a 1 , a 2 ) = {λ 0 a 0 + λ 1 a 1 + λ 2 a 2 | λ 0 ≤ 0, λ 1 ≥ 0, λ 2 ≥ 0, λ 0 + λ 1 + λ 2 = 1} et S(a 0 , a 1 , a 2 ) = {λ 0 a 0 + λ 1 a 1 + λ 2 a 2 | λ 0 ≥ 0, λ 1 ≥ 0, λ 2 ≥ 0, λ 0 + λ 1 + λ 2 = 1} illustrés par la figure 3.2 : T ∞ h a 0 a 1 a 2 S FIGURE 1.
3 -Une illustration du simplexe infini T ∞ et de son simplexe fini associé S dans le cas bidimensionnel.

Le domaine dilaté est ensuite partitionné en une partie bornée et une partie infinie :

R 2 = Ω b ∪ Ω ∞ avec Ω b ∩ Ω ∞ = , (1.19)
Le domaine infini Ω ∞ sera partitionné en un nombre fini de simplexes infinis :

Ω ∞ = ∪ M i =1 T ∞ i . (1.20) CHAPITRE 1. PRÉSENTATION DE LA THÈSE Nous introduisons ensuite le domaine ficitif Ω * = i nt (∪ M i =1 S i )\{a 0 } et l'inversion polygo- nale Φ : Ω ∞ → Ω * qui transforme le domaine infini Ω ∞ en le domaine fini Ω * définie par : ∀y ∈ Ω ∞ , Φ(y) = |h i | 4 (y.h i ) 2 y.
(1.21) où h i est le vecteur altitude associé au simplexe infini T ∞ i . Nous considérons maintenant un ensemble de triangulations :

J (R 2 ) = (T , T * )|T ∈ J 1 (Ω b ), T * ∈ J 2 (Ω * ) .
tel que :

-T est une triangulation conforme de Ω b .

-T * est une triangulation conforme de Ω * .

-Pour tout K ∈ T * , il existe i ≤ M tel que K ⊂ S i (en d'autres termes, T * est une union conforme des triangulations des sous-domaines S 1 , ..., S M ). Avec h K le diamètre d'un triangle arbitraire K de T ∪ T * nous posons h = h(T , T * ) = max K∈T ∪T * h K et nous établissons, sous certaines conditions, l'estimation d'erreur suivante :

u -u h W 1,2 0,0 (R 2 ) ≤ C * h k , (1.22)
Pour une constante C * ne dépendant pas de h.

Moyennant la transformation Φ nous transformons les intégrales posées sur Ω ∞ en des intégrales posées sur le domaine Ω * pour les problèmes en formulations faibles du Laplacien et de l'élasticité linéaire comme suit : Pour le problème de Poisson (1.7) et dans le cas où l'inclusion est circulaire, l'expression obtenue par séparation de variables de V (0) m nous permet en premier lieu d'effectuer une comparaison par rapport à la solution obtenue numériquement en traçant celle-ci dans certaines directions en comparant les deux solutions point par point (le correcteur tend bien vers 0 à l'infini et la coïncidence avec la solution analytique est parfaite) puis d'effectuer, en second lieu, une analyse de convergence dans l'espace à poids W 1,2 0,0 . Les ordres de convergence numériques sont proches des ordres théoriques conformément à l'estimation (1.22). Pour finir, nous présentons une stratégie numérique itérative de prise en compte d'une CHAPITRE 1. PRÉSENTATION DE LA THÈSE ou de plusieurs inclusions en ajoutant à chaque itération une condition de Neuman sur un bord ficitif représentant la zône d'influence de chaque inclusion en terme de correcteur au premier ordre et ceci jusqu'à atteindre l'équilibre. A chaque itération, ∇ ∇ ∇u (0) (0) est calculé pour le problème en Laplacien et (u (0) )(0) pour le problème en élasticité linéaire puis injectés comme termes sources dans les problèmes associés résolus par V (0) ensuite nous interpolons ε∂ n V (0) ou ε (V (0) ) sur le bord de correction fictif. En pratique, une base du correcteur est précalculée par forme d'inclusion (deux éléments pour le problème en Laplacien et quatre pour celui en élasticité lineaire) ce qui rend les itérations très peu couteuses. Nous évaluons les normes L 2 et H 1 de la différence entre la solution corrigée à chaque itération et une solution de référence (solution obtenue avec un maillage raffiné contenant une ou plusieurs inclusions selon les cas). Dans tous les cas traités la correction est significative dès la première itération pour se stabiliser au bout de quelques itérations.

Ω∞ α m ∇ ∇ ∇ψ 1 (y) •∇ ∇ ∇ψ 2 (y) dy = Ω * α m (∇ ∇ ∇ * ψ1 (y * )) T B -1 B -T ∇ ∇ ∇ * ψ2 (y * )|d et (B)| dy * , et Ω ∞ (ψ 1 (y)) : E m : (ψ 2 (y)) dy = Ω * (B -T ∇ ∇ ∇ * ψ1 (y * ) +∇ ∇ ∇ * ψ1 (y * )B -1 ) : E m : (B -T ∇ ∇ ∇ * ψ2 (y * ) +∇ ∇ ∇ ψ2 (y * )B -1 )
W(I 1 ) = µ k 2 (I 1 -3) sur R ε k (k ∈ { f , m}). ( 1
τ = ∂W(F) ∂F -qF -T = µ k F -qF -T dans R ε k (k ∈ { f , m}).
(1.26)

Le système fort à résoudre est alors le suivant :

Di v (τ) = 0 dans R ε k (k ∈ { f , m}), (1.27a) x f (X 1 , X 2 , X 3 ) = x m (X 1 , X 2 , X 3 ) ∀(X 1 , X 2 ) ∈ Γ ε et ∀X 3 ,
(1.27b)

τ f (X 1 , X 2 , X 3 )N = τ m (X 1 , X 2 , X 3 )N ∀(X 1 , X 2 ) ∈ Γ ε et ∀X 3 ,
(1.27c)

u ε m = U ∀(X 1 , X 2 ) ∈ Γ et ∀X 3 sur le bord extérieur R ε .
(1.27d)

La résolution de ce système nous amène au système suivant :

∆u ε f = d µ f dans Ω ε f , (1.28a) ∆u ε m = d µ m dans Ω ε m , (1.28b
)

u ε f = u ε m sur Γ ε , (1.28c) µ f ∂ n u ε f = µ m ∂ n u ε m sur Γ ε , (1.28d) u ε m = U sur Γ.
(

∆d ε 0, f (x) = d µ f - d µ m dans Ω ε f , (1.29a) ∆d ε 0,m (x) = 0 dans Ω ε m , (1.29b) d ε 0, f (x) = d ε 0,m (x) sur Γ ε , (1.29c) µ f ∂ n d ε 0, f (x) = µ m ∂ n d ε 0,m (x) + (µ m -µ f )∂ n u (0) (x) sur Γ ε , (1.29d) d ε 0,m (x) = 0 sur Γ.
(1.29e)

Nous obtenons ainsi le développement asymptotique à l'ordre K suivant :

u [K] ε,k (x) = u (0) (x) + K i =2 ε i v (i ) (x) + K i =1 ε i V (i -1) k ( x ε ) + ε 2 W (0) ( x ε ).
(1.30) La seconde partie concerne une généralisation à des potentiels quelconques et nous y présentons une manière d'approximer numériquement un premier ordre. Nous commençons par présenter le problème fort résolu par le résidu à l'ordre zéro d ε 0,k :

Di v (τ(u (0) + d ε 0,k ) -τ m (u (0) )) = 0 dans Ω ε k (k = f , m), (1.31a) τ f N f = τ m N m sur Γ ε , (1.31b) d ε 0, f = d ε 0,m sur Γ ε , (1.31c) d ε 0, f = 0 sur Γ.
(1.31d) Après scaling, certains développements et linéarisations puis en explicitant le résidu à l'ordre 0 par D ε 0,k = εV (0) k + O(ε 2 ) nous obtenons la formulation faible résolue par le premier correcteur V (0) k (y) suivante :

Trouver V (0) k : R 2 -→ R 3 tel que : Ω 1 f T f (u (0) )(0)∇ ∇ ∇ y V (0) f (y) : ∇ ∇ ∇ y H f (y) dy + Ω ∞ T m (u (0) )(0)∇ ∇ ∇ y V (0) m (y) : ∇ ∇ ∇ y H m (y) dy = Ω 1 f (τ f (u (0) )(0) -τ m (u (0) )(0)) : ∇ ∇ ∇ y H f (y) dy, (1.32a)
H f = H m sur Γ 1 , (1.32b)
V (0) f (y) = V (0) m (y) sur Γ 1 ,
(1.32c)

V (0) m → 0 quand |y| → ∞, (1.32d)
H m → 0 quand |y| → ∞.
(1.32e)

Pour toute fonction test H satisfaisant les conditions au dessus.

Pour obtenir le premier ordre V (0) k (y) il nous faudra l'estimer alors numériquement par les méthodes des éléments finis et des éléments finis inversés comme ceci a été fait précédemment.

Remarque 1.6.1. Nous revenons ici à un développement analogue au cas linéaire pour de l'élasticité anisotrope. La non linéarité se limite juste au calcul du champ sans inclusion u (0) .

Conclusion et perspectives

Ce travail a fait l'objet de l'étude de l'intéraction entre une inclusion fine et le domaine dans lequel celle-ci évolue. L'étude a été motivée par les applications dans les pneumatiques où l'agencement inclusion-caoutchouc représente une pièce maitresse de la phase de conception mais le champ d'applications des résultats obtenus reste très vaste. L'étude effectuée est subdivisée en deux grands volets. Le premier concerne une partie modélisation où il a été question de construire en premier lieu un développement CHAPITRE 1. PRÉSENTATION DE LA THÈSE asymptotique multi-échelle modélisant l'influence d'une inclusion en partant du champ sans inclusion et en corrigeant à chaque ordre par des correcteurs posés sur des domaines indépendants de la dimension très petite de l'inclusion. L'étude a été effectuée d'une façon rigoureuse pour un problème de Poisson et un problème en élasticité linéaire puis des généralisations à des lois de comportement non linéaires ont été proposées (l'hyperélasticité) à titre prospectif. Tout au long du travail la géométrie a été réduite au cas bidimensionnel.

Le deuxième volet a concerné l'approximation du premier ordre obtenu dans le cas des comportements linéaires par les méthodes des éléments finis classique et celle des éléments finis inversés. Une analyse de convergence a été effectuée et une stratégie itérative de prise en compte de plusieurs inclusions présentée. Les résultats numériques sont satisfaisants et valident le modèle construit.

Les perspectives de ce travail sont de continuer l'étude entamée pour l'hyperélasticité et de présenter des résultats numériques associés à différents ordres. L'extension à une géométrie tridimensionnelle est aussi envisagée. L'objectif final convoité est d'étendre l'étude à un cas tridimensionnel non linéaire sous sollicitations dynamiques.

Références

Introduction

The aim of this Chapter is to present an asymptotic expansion of the influence of a small inclusion of different stiffness in an elastic media. The applicative interest of this study is to provide tools to take into account this influence from the deformation without inclusion thanks to additive terms that can be precalculated and which depend only on the shape of the inclusion. We treat two problems : an anti-plane linearized elasticity problem and a plane strain one. On every expansion order we provide corrective terms modeling the influence of the inclusion using techniques of scaling and multi-scale asymptotic expansions. The built expansion is validated by comparing it to a test case obtained by solving the Poisson transmission problem in the case of a circular shape of the inclusion using variable separation method. Proofs of existence and uniqueness on unbounded domains of our fields are also adapted to the bidimensional Poisson problem and the linear elasticity one. This Chapter has been accepted for publication in journal of elasticity ARFAOUI et collab.

[2016] and its outline is as follows : Asymptotic analyses for two transmission problems of interest are presented in Sections 2.2 and 2.3, respectively. For each problem, existence, uniqueness and error estimates are established. Since we are led to deal with functions defined on infinite domains we reconsidered proof of existence and uniqueness using the generalized Hardy's inequalities and adequate weighted Sobolev spaces. Finally, some conclusions and perspectives are presented in Section 2.4.

Asymptotic expansion for the antiplane problem

The aim of this section is to build an asymptotic expansion u [K] ε,k of a given order K for the solution to the antiplane problem (1.1) having the following form for k ∈ { f , m} :

u [K] ε,k (x) = u (0) (x) + K i =2 ε i v (i ) (x) + K i =1 ε i V (i -1) k ( x ε ), (2.1)
where u (0) is the solution to Problem (1.1) when ε = 0 (i.e. without any inclusion), the functions v (i ) are defined on Ω and the functions V (i ) on the whole R 2 .

Remark 2.2.1. An important aspect for the applications is that the functions v (i ) and V (i ) , solutions to some problems which will be detailled further on, do not depend on ε and depend linearly on u (0) only by the expansion of its gradient at x = 0. For a given shape of the inclusion, this make it possible to pre-compute once for all a basis for each v (i ) and V (i ) . This pre-computation can be done either analitically in the circular case (see Section

2.2.5.2), by an integral formulation using a Green kernel (as for instance in BONNETIER et TRIKI [2009]) or by a numerical approximation with a boundary element method or a Galerkin method.

The rest of this section is subdivised into the following sub-sections : The asymptotic expansion construction is detailed in Section 2.2.1. Then, existence and uniqueness of the V (i ) is stated in Section 2.2.2 and their behavior at infinity is specified in Section 2.2.3. This allows, in Section 2.2.4, to deduce an optimal error estimate for asymptotic expansion 2.1. Finally in Section 2.2.5 an explicit computation is performed using the separation of variable method in the case of a circular body and a circular inclusion. In addition to giving the explicit form of the expansion in this case, it allows us to compare the solution to its expansion and conclude that it is in accordance with the theoretical results.

Building the asymptotic expansion

We build the asymptotic expansion step by step. We need first to introduce the difference

d ε 0,k = u ε k -u (0) for k ∈ { f , m} which is solution to : α f ∆d ε 0, f (x) = 0 in Ω ε f , (2.2a) α m ∆d ε 0,m (x) = 0 in Ω ε m , (2.2b) d ε 0, f (x) = d ε 0,m (x) on Γ ε , (2.2c) α f ∂ n d ε 0, f (x) = α m ∂ n d ε 0,m (x) + (α m -α f )∂ n u (0) (x) on Γ ε , (2.2d) d ε 0,m (x) = 0 on Γ.
(2.2e)

Note the emergence of the term (α m -α f )∂ n u (0) in the transmission equation (4.37d).

The original problem is posed on domains depending on ε. In order to work with a fixed domain and a fixed interface, we introduce as in DAMBRINE et VIAL [2005]; VIAL

[2003] the scaling y = x ε . Let Ω m,ε = Ω ε m ε . Then, for k ∈ { f , m}, D ε 0,k (y) = d ε 0,k (x)
is the solution to the following scaled problem :

α f ∆D ε 0, f (y) = 0 in Ω 1 f , (2.3a)
α m ∆D ε 0,m (y) = 0 in Ω m,ε , (2.3b)
D ε 0, f (y) = D ε 0,m (y) on Γ 1 , (2.3c)
α f ∂ n D ε 0, f (y) = α m ∂ n D ε 0,m (y) + ε(α m -α f )∂ n u (0) (εy) on Γ 1 , (2.3d) D ε 0,m (y) = 0 on Γ ε . (2.3e)
Obtaining equation (2.3d) deserves to be detailed. Indeed, using the scaling y = x ε we obtain the relation between the derivation with respect to the fast variable y and that with respect to the slow variable x is : ∂ ∂y = ε ∂ ∂x

. From there and using that the transformation preserves the normal we get that :

α f ∂ n D ε 0, f (y) = εα f ∂ n d ε 0, f (x).
Finally, using equation (2.2d) we obtain the desired result. Now, as ε is considered to be small compared to the size of the domain and to eliminate completely the dependence on ε, we approximate Problem (2.3) by the following problem posed on the dilated domain Ω ∞ = lim ε→0 Ω m,ε as shown in Fig. 1.2 :

α f ∆D ε 0, f (y) = 0 in Ω 1 f , (2.4a) α m ∆D ε 0,m (y) = 0 in Ω ∞ , (2.4b) D ε 0, f (y) = D ε 0,m (y) on Γ 1 ,
(2.4c)

α f ∂ n D ε 0, f (y) = α m ∂ n D ε 0,m (y) + ε(α m -α f )∂ n u (0) (εy) on Γ 1 , (2.4d) D ε 0,m (y) → 0 when |y| → ∞. (2.4e)
This allows us now to introduce the problem on the dilated geometry whose solutions

CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA
will be the functions V (p) of the expansion (2.1) :

α f ∆V (p) f (y) = 0 in Ω 1 f , (2.5a) α m ∆V (p) m (y) = 0 in Ω ∞ , (2.5b) V (p) f (y) = V (p) m (y) on Γ 1 ,
(2.5c)

α f ∂ n V (p) f (y) = α m ∂ n V (p) m (y) + ϕ (p) (y) on Γ 1 , (2.5d) V (p) m (y) → 0 when |y| → ∞.
(2.5e)

The functions V (p) are introduced to offset terms appearing on the boundary of Γ ε where the data ϕ (p) (y) will be adapted for each order p. In fact, the smaller these terms, the smaller the residual and the expansion is driven to a higher order. The existence and uniqueness of the functions V (p) will be established later in Section 2.2.2. From the regularity assumptions on the domain, we can write a Taylor expansion of ∇u (0) (εy) for y ∈ Γ 1 :

∇ ∇ ∇u (0) (εy) = K i =0 ε i ω (i ) 0 (y) + O(ε K+1 ).
(2.6)

First-order approximation

The first-order approximation is constructed by compensating the first term in the expansion of (α

m -α f )∂ n u (0) (εy). With V (0) solving Problem (2.5) for ϕ (0) (y) = (α m -α f )ω (0) 0 (y)• n, where ω (0) 0 (y) = ∇ ∇ ∇u (0) (0)
, we update the residual :

d ε 1,k (x) = u ε k (x) -u (0) (x) -εV (0) k ( x ε ), k ∈ { f , m}.
(2.7)

The update of the difference allows us to obtain the first-order approximation :

u [1] ε,k (x) = u (0) (x) + εV (0) k ( x ε ), k ∈ { f , m}. (2.8) Then, for k ∈ { f , m}, d ε 1,k is solution to : α f ∆d ε 1, f (x) = 0 in Ω ε f , (2.9a) α m ∆d ε 1,m (x) = 0 in Ω ε m , (2.9b) d ε 1, f (x) = d ε 1,m (x) on Γ ε , (2.9c) α f ∂ n d ε 1, f (x) = α m ∂ n d ε 1,m (x) + (α m -α f )(∂ n u (0) (x) -ω (0) 0 (y)) on Γ ε , (2.9d) d ε 1,m (x) = -εV (0) m ( x ε ) on Γ.
(2.9e)

Second-order approximation

We have built a first-order approximation of the solution with u (0) and V (0) independent of ε. The next step is to build the second-order approximation. At this stage, we must do the same thing we did for the term that appeared on the interface Γ ε with the term that appears this time on Γ. In other words, we have two terms to compensate :

(α m -α f )(∂ n u (0) - ω (0)
0 ) which has already been compensated in the first-order and that has yet to be, and the new term -εV (0) m ( x ε ) on Γ. In fact, by analyzing things a little deeper, we can see that

CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA the term on Γ ε is of order ε then it is in O(ε 2 ) according to the H -1 2 norm ( see (2.31)).
The term on Γ is of order ε 2 (expansion of V at infinity see Section 2.2.3) and there is no scaling here so it is in O(ε 2 ) according to the H 1 2 norm. By doing so we got the optimal error estimate in the previous first order approximation. So, to continue having this optimality the two terms have to be compensated in the same time from now on.

The functions introduced by Problem (2.5) are used to correct the term on Γ ε and we will introduce equivalent functions to play the same role but in this case on the boundary Γ.

To this aim, we consider the following problem :

α m ∆ x v (q) (x) = 0 in Ω, (2.10a) v (q) (x) = ψ (q) (x)
on Γ.

(2.10b)

First, let us expand V (0) m at infinity as follows :

V (0) m (y) = K i =1 a (i ) 0 (θ) |y| i + O 1 |y| K+1 = K i =1 ε i a (i ) 0 (θ) |x| i + O(ε K+1 ), x ∈ Γ, (2.11)
where θ is the angular polar coordinate. This expansion at infinity is justified in Section 2.2.3. Now, with V (1) solution of Problem (2.5) for ϕ

(1) (y) = (α m -α f )ω (1) 0 (y) • n = (α m -α f )ω (1) 0 (y) and v (2) solution of Problem (2.10) for ψ (2) (x) = - a (1) 0 (θ) |x|
we write the second-order term :

d ε 2,k (x) = d ε 1,k (x) -ε 2 v (2) (x) -ε 2 V (1) k ( x ε ), k ∈ { f , m}.
(2.12)

The update of the difference allows us to obtain the second-order approximation :

u [2] ε,k (x) = u (0) (x) + εV (0) k ( x ε ) + ε 2 v (2) (x) + ε 2 V (1) k ( x ε ), k ∈ { f , m}. (2.13)
Then, for k ∈ { f , m}, d ε 2,k solves the following problem :

α f ∆d ε 2, f (x) = 0 in Ω ε f , (2.14a) α m ∆d ε 2,m (x) = 0 in Ω ε m , (2.14b) d ε 2, f (x) = d ε 2,m (x) on Γ ε , (2.14c) α f ∂ n d ε 2, f (x) = α m ∂ n d ε 2,m (x) + (α m -α f )(∂ n u (0) (εy) -ω (0) 0 (y) -εω (1) 0 (y)) +ε 2 (α m -α f )(∂ n v (2) (εy)) on Γ ε , (2.14d) d ε 2,m (x) = ε(ε a (1) 0 (θ) |x| -V (0) m ( x ε )) -ε 2 V (1) m ( x ε ) on Γ.
(2.14e)

Third-order approximation

We expand as done before ∇v (2) (εy) as follows :

∇ ∇ ∇v (2) (εy) = K-2 i =0 ε i ω (i ) 2 (y) + O(ε K+1 ). (2.15)
We expand too V (1) m (see Section 2.2.3) :

V (1) m (y) = K i =1 a (i ) 1 (θ) |y| i + O 1 |y| K+1 = K i =1 ε i a (i ) 1 (θ) |x| i + O(ε K+1 ), x ∈ Γ.
(2.16)

CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA With V (2) solution of Problem (2.5) for ϕ (2) (y) = (α m -α f )(ω (2) 0 (y) + ω (0) 2 (y)) • n = (α m - α f )(ω (2)
0 (y) + ω (0) 2 (y)) and v (3) solution of Problem (2.10) for ψ (3)

(x) = - a (1) 1 (θ) |x| - a (2) 0 (θ)
|x| 2 we deduce the third-order residual :

d ε 3,k (x) = d ε 2,k (x) -ε 3 v (3) (x) -ε 3 V (2) k ( x ε ), k ∈ { f , m}.
(2.17)

The update of the difference allows us to obtain the third-order approximation :

u [3] ε,k (x) = u (0) (x) + εV (0) k ( x ε ) + ε 2 V (1) k ( x ε ) + ε 2 v (2) (x) + ε 3 v (3) (x) + ε 3 V (2) k ( x ε ), k ∈ { f , m}.
(2.18) Then, for k ∈ { f , m}, d ε 3,k solves the following problem :

α f ∆d ε 3, f (x) = 0 in Ω ε f , (2.19a) α m ∆d ε 3,m (x) = 0 in Ω ε m , (2.19b) d ε 3, f (x) = d ε 3,m (x) on Γ ε , (2.19c) α f ∂ n d ε 3, f (x) = α m ∂ n d ε 3,m (x) +(α m -α f )(∂ n u (0) (εy) -ω (0) 0 (y) -εω (1) 0 (y) -ε 2 ω (2) 0 (y)) +ε 2 (α m -α f )(∂ n v (2) (εy) -ω (0) 2 (y)) +ε 3 (α m -α f )(∂ n v (3) (εy)) on Γ ε , (2.19d) d ε 3,m (x) = ε(ε 2 a (2) 0 (θ) |x| 2 -ε a (1) 0 (θ) |x| -V (0) m ( x ε )) +ε 2 (ε a (1) 1 (θ) |x| -V (1) m ( x ε )) -ε 3 V (2) m ( x ε ) on Γ.
(2.19e)

Fourth-order approximation

In order to conclude about the general form of functions ϕ (n) involved in problems solved by V (n) and functions ψ (n) involved in problems solved by v (n) we push the expansion to the fourth order. We expand as done before ∇ ∇ ∇v (3) (εy) as follows :

∇ ∇ ∇v (3) (εy) = K-3 i =0 ε i ω (i ) 3 (y) + O(ε K+1 ).
(2.20) and V (2) m as follows (see Section 2.2.3) :

V (2) m (y) = K i =1 a (i ) 2 (θ) |y| i + O 1 |y| K+1 = K i =1 ε i a (i ) 2 (θ) |x| i + O(ε K+1 ), x ∈ Γ.
(2.21) With V (3) solution of Problem (2.5) for ϕ (3) (3) 0 (y)+ω (1) 2 (y)+ω (0) 3 (y)) and v (4) solution of Problem (2.10) for ψ (4)

(y) = (α m -α f )(ω (3) 0 (y) + ω (1) 2 (y) + ω (0) 3 (y)) • n = (α m -α f )(ω
(x) = - a (1) 2 (θ) |x| - a (2) 1 (θ) |x| 2 - a (3) 0 (θ)
|x| 3 we deduce the fourth-order residual :

d ε 4,k (x) = d ε 3,k (x) -ε 4 v (4) (x) -ε 4 V (3) k ( x ε ), k ∈ { f , m}. (2.22) CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA
The update of the difference allows us to obtain the fourth-order approximation :

u [4] ε,k (x) = u (0) (x) + εV (0) k ( x ε ) + ε 2 V (1) k ( x ε ) + ε 2 v (2) (x) + ε 3 V (2) k ( x ε ) + ε 3 v (3) (x) +ε 4 v (4) (x) + ε 4 V (3) k ( x ε ), k ∈ { f , m}. (2.23)
Then, for k ∈ { f , m}, d ε 4,k solves the following problem :

α f ∆d ε 4, f (x) = 0 in Ω ε f , (2.24a) α m ∆d ε 4,m (x) = 0 in Ω ε m , (2.24b) d ε 4, f (x) = d ε 4,m (x) on Γ ε , (2.24c) α f ∂ n d ε 4, f (x) = α m ∂ n d ε 4,m (x) +(α m -α f )(∂ n u (0) (εy) -ω (0) 0 (y) -εω (1) 0 (y) -ε 2 ω (2) 0 -ε 3 ω (3) 0 (y)) +ε 2 (α m -α f )(∂ n v (2) (εy) -ω (0) 2 (y) -εω (1) 2 (y)) +ε 3 (α m -α f )(∂ n v (3) (εy) -ω (0) 3 (y)) + ε 4 (α m -α f )(∂ n v (4) (εy)) on Γ ε , (2.24d) d ε 4,m (x) = ε(ε 3 a (3) 0 (θ) |x| 3 + ε 2 a (2) 0 (θ) |x| 2 + ε a (1) 0 (θ) |x| -V (0) m ( x ε )) +ε 2 (ε 2 a (2) 1 (θ) |x| 2 + ε a (1) 1 (θ) |x| -V (1) m ( x ε )) +ε 3 (ε a (1) 1 (θ) |x| -V (2) m ( x ε )) -ε 4 (V (3) m ( x ε )) on Γ. (2.24e)

The K-th order approximation

From the previous section we are able to deduce the general forms of ϕ (n) and ψ (n) so we can update our residual to the order n ≥ 2 as follows :

d ε n,k (x) = d ε n-1,k (x) -ε n v (n) (x) -ε n V (n-1) k ( x ε ), k ∈ { f , m}, (2.25)
where V (n-1) solves Problem (2.5) for ϕ

(n-1) (y) = (α m -α f ) p+q=n-1 p≥0,q =1 ω (p) q (y) • n and v (n) solves Problem (2.10) for ψ (n) (x) = -p+q=n-1 p≥1,q≥0 a (p) q (θ)
|x| p . We are now able to build an approximation u [K] ε (x) to the order K ≥ 2 :

u [K] ε,k (x) = u (0) (x) + K i =2 ε i v (i ) (x) + K i =1 ε i V (i -1) k ( x ε ), k ∈ { f , m}. (2.26)
Remark 2.2.2. According to the power of ε this is not the expansion to the order K in the true sense. Indeed, it is the expansion to the order K + 1 truncated of a term. To clearly see this, let us take the example where K is equal to 1. By analyzing the expression of u [1] ε (x), we notice that initially the expansion has a term in ε 2 (εV (0) ( x ε )). However, an ε 2 order term will appear at the second order (ε 2 v (2) ). So, in order to obtain all the terms in ε K , we have to push the expansion at the order K and truncate the term

ε K V K . CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA Now, replacing u [K] ε (x) in Problem (1.1) we obtain the following problem ( f ε (x) = g ε (y)) : α f ∆u [K] ε (x) = 0 in Ω ε f , (2.27a) α m ∆u [K] ε (x) = 0 in Ω ε m , (2.27b) u [K] ε (x) = 0 on Γ ε , (2.27c) α∂ n u [K] ε (x) = f ε (x) on Γ ε , (2.27d) u [K] ε (x) = h ε (x) + U on Γ, (2.27e)
By setting u (0) = v (0) and using the expansion (2.26) we obtain :

g ε (y) = K k=0 k =1 ε k (α m -α f )∂ n v (k) (εy) + K k=1 ε k-1 α∂ n V (k-1) (y) Γ 1 .
(2.28)

Using the expressions of ϕ (i -1) (y) we obtain the expression of the jump

α∂ n V (k-1) (y) Γ 1 = (α f -α m ) p+q=k-1 p≥0,q =1 ω (p) q (y)
• n and by a simple recurrence we deduce the expressions of

∂ n v (k) (εy) = K-k l =0 ε l ω (l ) k (y) • n + O(ε K+1-k
). Now replacing these expressions in (2.28) we obtain :

g ε (y) = K k=0 k =1 ε k [(α m -α f ) K-k l =0 ε l ω (l ) k (y) • n + O(ε K+1-k )] + K k=1 ε k-1 [(α f -α m ) p+q=k-1 p≥0,q =1 ω (p) q (y) • n]. (2.29)
By operating two changes of idex n = k + l for the first part of (2.29) and n = k -1 for the second part of it we obtain :

g ε (y) = K n=0 ε n (α m -α f ) k+l =n l ≥0,k =1 ω (l ) k (y) • n + O(ε K+1 ) + K-1 n=0 ε n [(α f -α m ) p+q=n p≥0,q =1 ω (p) q (y) • n] = ε K f (y). (2.30) with f (y) = ε K (α m -α f ) k+l =K p≥0,k =1 ω (l ) k (y) • n.
Then, we obtain :

g ε (.) L 2 (Γ 1 ) = O(ε K ). (2.31)
Using one more time the expansion (2.26) we obtain :

h ε (x) = K k=2 ε k (- p+q=k-1 p≥1,q≥0 a (p) q (θ) |x| p ) + K k=1 ε k ( K-k l =1 ε l a (l ) k (θ) |x| l + O(ε K+1-k )). (2.32)
By applying exactly the same changes of indices we obtain :

h ε (x) = ε K+1 ψ(x).
(2.33) and finally :

h ε (x) H s (Γ) = O(ε K+1 ). (2.34) Consequently, d ε K,k = u ε k -u [K] ε,k (x), k ∈ { f , m} satisfy : α f ∆d ε K, f (x) = 0 in Ω ε f , (2.35a) α m ∆d ε K,m (x) = 0 in Ω ε m , (2.35b) d ε K, f (x) = d ε K,m (x) on Γ ε , (2.35c) α f ∂ n d ε K, f (x) = α m ∂ n d ε K,m (x) + O(ε K ) on Γ ε , (2.35d) d ε K,m (x) = O H 1 2 (Γ) (ε K+1 )
on Γ.

(2.35e)

Existence and uniqueness in unbounded domains for the Poisson problem

A weak formulation of Problem (2.5) writes :

α f Ω 1 f ∇ ∇ ∇V (p) f •∇ ∇ ∇w dx + α m Ω ∞ ∇ ∇ ∇V (p) m •∇ ∇ ∇w dx = Γ 1 ϕw dσ, ∀w ∈ V , (2.36) where V = v ∈ W 1,2 0,0 ; Γ R 0 v = 0 with 0 < R 0 < R an arbitrary real number and W 1,2 0,0
is introduced in the appendix. In fact, the space W 1,2 0,0 contains constant functions and the Poincaré-type inequality is established on W 1,2 0,0 /P 0 , the quotient space of the adequate weighted Sobolev space with constant functions. We can see that the bilinear form

a(z, w) = α f Ω 1 f ∇ ∇ ∇z •∇ ∇ ∇w dx + α m Ω ∞ ∇ ∇ ∇z •∇ ∇
∇w dx is coercive on V . In addition, the bilinear form a(•, •) and the linear one l (w) = Γ 1 ϕw dσ are continuous on V . Then, the existence and uniqueness of the functions V (p) in V are guaranteed by the Lax-Milgram theorem. is smooth and can be written as the sum of its Fourier series :

Expansion of the functions V

V (p) m (R, θ) = α (p) 0 + n≥1 α (p) n cos(nθ) + β (p) n sin(nθ).
(2.37) Using Poisson's Kernel, we then get that :

V (p) m (r, θ) = α (p) 0 + n≥1 R n r n α (p) n cos(nθ) + β (p)
n sin(nθ) .

(2.38)

The behavior of the function V

Error estimate for the antiplane problem

In this section we establish an error estimate for the studied problem by increasing the H 1 -norm by powers of ε (characteristic length of the inclusion). First, we begin by

CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA considering the following problem for a given f ε ∈ H -1 2 (Γ ε ) and h ε ∈ H 1 2 (Γ) ( f ε (x) = g ε (y)) : α f ∆ x z f (x) = 0 in Ω ε f , (2.39a) α m ∆ x z m (x) = 0 in Ω ε m , (2.39b) z f (x) = z m (x) on Γ ε ,
(2.39c)

α f ∂ n z f (x) = α m ∂ n z m (x) + f ε (x) on Γ ε , (2.39d) z m (x) = h ε (x)
on Γ.

(2.39e)

By passing to a weak form we obtain the following equivalent problem : Find z ∈ H 1 (Ω)

with z = h on Γ such that

α f Ω ε f ∇ ∇ ∇z •∇ ∇ ∇w dx + α m Ω ε m ∇ ∇ ∇z •∇ ∇ ∇w dx = Γ ε f ε w dσ, ∀ω ∈ H 1 0 (Ω). ( 2
satisfying z H 1 (Ω) ≤ C h ε H 1 2 (Γ) + g ε L 2 (Γ 1 ) ,
where C > 0 is a constant independent of ε.

Proof. The surjectivity of the trace operator guarantees the existence of a continuous lif-

ting operator H ∈ H 1 (Ω) of h ∈ H 1 2 (Γ) (see e.g. ERN et GUERMOND [2004]; FORTIN et GARON [2011]) : Then z = z -H ∈ H 1 0 (Ω)
solves the following variational equation :

α f Ω ε f ∇ ∇ ∇z •∇ ∇ ∇w dx + α m Ω ε m ∇ ∇ ∇z •∇ ∇ ∇w dx - Γ ε f ε w dσ + α f Ω ε f ∇ ∇ ∇H •∇ ∇ ∇w dx + α m Ω ε m ∇ ∇ ∇H •∇ ∇ ∇w dx = 0, ∀w ∈ H 1 0 (Ω). (2.41)
The bilinear form associated to this weak formulation [START_REF] Fortin | Weighted Sobolev spaces for Laplace's equation in R n[END_REF]. Then, from the Lax-Milgram theorem we deduce the existence and uniqueness of z solution of (2.41). The coercivity of this bilinear form yields :

a(u, v) = α f Ω ε f ∇ ∇ ∇u • ∇ ∇ ∇v dx + α m Ω ε m ∇ ∇ ∇u •∇ ∇ ∇v dx is coercive, i.e. α v 2 H 1 (Ω) ≤ a(v, v), see ERN et GUERMOND [2004]; FOR- TIN et
α z 2 H 1 (Ω) ≤ Γ ε f ε z dσ -α f Ω ε f ∇ ∇ ∇H •∇ ∇ ∇z dx -α m Ω ε m ∇ ∇ ∇H •∇ ∇ ∇z dx, (2.42)
and hence using the estimation (A.17) demonstrated in the appendix A.2 we obtain

α z 2 H 1 (Ω) ≤ C g ε L 2 (Γ 1 ) z H 1 (Ω) + ∇H L 2 (Ω) ∇z L 2 (Ω) .
(2.43)

Using the facts that ∇z L 2 (Ω) ≤ C z H 1 (Ω) and ∇H L 2 (Ω) ≤ C H H 1 (Ω) we obtain :

α z 2 H 1 (Ω) ≤ C g ε L 2 (Γ 1 ) + H H 1 (Ω) z H 1 (Ω) .
(2.44)

The continuity of the lifting operator gives us that

H H 1 (Ω) ≤ C h H 1 2 (Γ)
, so we can write :

z H 1 (Ω) ≤ C g ε L 2 (Γ 1 ) + h ε H 1 2 (Γ)
.

(2.45)

CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA

We are now ready to state the following result :

Theorem 2.2.4. There exists a constant C > 0 independent of ε such that :

d ε K H 1 (Ω) ≤ Cε K , for all K ∈ N, where d ε K is the residual introduced in Section 2.2. Proof. With z = d ε K we have g ε ∈ O L 2 (Γ 1 ) (ε K ) and h ε ∈ O H 1 2 (Γ) (ε K+1 ) (see Section 2.2.1.5).
Remark 2.2.5. The estimation obtained is not optimal but the existence of the expansion at any order and the techniques developed in TORDEUX [2004] may allow to establish an optimal estimate outside a neighborhood r ≥ r 0 > 0 of the inclusion.

Solution for a circular inclusion using the separation of variable method

In this section, we solve Problem (1.1) by the separation of variable method then we solve Problem (1.7) by the same method and compare the two results for the first-order of the asymptotic expansion. We will use the same geometrical configuration but with a circular shape for the interface Γ ε and the boundary Γ 1 (see Fig. 2.1).

.

0 2ε Γ Ω ε m R Ω ε f Γ ε . FIGURE 2.
1 -Circular inclusion and domain.

Antiplane problem solve using the separation of variable method

Considering the polar coordinates (r, θ) and the following particular form for the

u ε k , k ∈ { f , m} : u ε k (r, θ) = f k (r )g k (θ), (2.46)
we obtain from equations (4.35a) and (4.35b) the following equations :

r 2 f k (r ) f k (r ) + r f k (r ) f k (r ) = - g k (θ) g (θ) , k ∈ { f , m}. (2.47)
The left hand sides of these equations are independent of θ and the right ones are independent of r . Then they are both equal to the same constant, called separation constant,

CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA
that can be proven to be equal to n 2 k with n k ∈ N. This leads to the system of two decoupled equations :

r 2 f k (r ) + r f k (r ) -n 2 k f k (r ) = 0, (2.48a) g k (θ) + n 2 k g k (θ) = 0. (2.48b)
Using the equation (4.35c) we can see that g 1 (θ) = g 2 (θ) and then n k = n. We observe now that we can search solutions of the forms of r α and thus we obtain :

f f (r ) = C 1 r n + cr -n , (2.49a) f m (r ) = C 2 r n + C 3 r -n . (2.49b)
The solution must be bounded at the origin, hence c = 0. The next step is to use the continuity of the fields and their normal derivatives throughout the interface Γ ε to determine the constants C 1 , C 2 and C 3 . In fact, from equations (4.35c) and (4.35d) we obtain with

β = α m -α f α m +α f and h = ε R the following expressions : C 1 (n) = R -n (1 + β) 1 + βh 2n , (2.50a) C 2 (n) = R -n 1 + βh 2n , (2.50b) C 3 (n) = R n βh 2n 1 + βh 2n .
(2.50c)

By expanding C 1 , C 2 and C 3 with respect to the small parameter h and taking n = 1 we can bring up the first contribution of the inclusion in the expression of matrix field :

u ε m = u (0) m + h 2 u (2) m .
(2.51)

Since u (0) m , u (2) m are independent of ε so we can conclude that the first contribution of the inclusion is of order ε 2 .

Solving V (0) by the separation of variable method

Thanks to the separation of variable method and using polar coordinates (r s , θ) (r s = r ε ) we decompose the functions V (0) k as follows :

V (0) k (r s , θ) = n≥1 f (k) n (r s ) cos(nθ) + g (k) n (r s ) sin(nθ), k ∈ { f , m}.
(2.52) Injecting (2.52) in (4.38a) and (4.38b) and solving the differential equations we obtain :

f k n = A k r n s + a k r -n s ,
(2.53a)

g k n = B k r n s + b k r -n s .
(2.53b)

Using the interface conditions prescribed by (4.38c) and (4.38d) we obtain :

V (0) m (r s , θ) = α m -α f r s a α f + α m cos θ + b α f + α m sin θ . (2.54)
with ∇ ∇ ∇u (0) (0) = (a b) t . Thus, Solving directly Problem (1.1) by the separation of variable method we obtain that the first contribution of the inclusion is in ε 2 . Now injecting (2.54) in (2.8) we obtain :

u [1] ε,m (r, θ) = u (0) + ε 2 r (α m -α f )a α f + α m cos θ + (α m -α f )b α f + α m sin θ . (2.55)
The first contribution is indeed in ε 2 which is in accordance with the result presented in Section 2.2.5.1.

Asymptotic expansion for the plane strain linear elasticity problem

We now propose an adaptation of the construction of the asymptotic expansion to the vector valued case of the plain strain elasticity problem (1.2). One of the main differences will be the appearance of an additional term in the expansion related to the difference of the constitutive laws in the inclusion and in the matrix. The form of expansion is as follows :

u [K] ε,k (x) = u (0) (x) + K i =2 ε i v (i ) (x) + K i =1 ε i V (i -1) k ( x ε ) + K i =1 ε i +1 Z (i -1) k ( x ε ). (2.56)
where u (0) denotes the solution to Problem (1.2) still for ε = 0, the functions v (i ) are defined on Ω and the V (i ) k and Z (i ) k functions on the whole R 2 .

The structure of this section globally follows the one for the Poisson problem in Section 2.2 apart from the fact that we do not propose a resolution by separation of variables, mainly because it gives no easily exploitable result in that case. Of course, an approximation by boundary elements or a Galerkin method can still be considered. The result on the behavior at infinity of functions V (i ) and Z (i ) is given in section 2.3.3 but only in the isotropic case. This also limits to the isotropic case the error estimate given in Section 2.3.4.

Building the asymptotic expansion

The difference

d ε 0,k = u ε k -u (0) (k ∈ { f , m}) is now solution to : div x σ f (d ε 0, f )(x) = -div x (E f : (u (0) ))(x) in Ω ε f , (2.57a) div x σ m (d ε 0,m (x)) = 0 in Ω ε m , (2.57b) d ε 0, f (x) = d ε 0,m (x) on Γ ε , (2.57c) σ f (d ε 0, f (x)) • n f = -σ m (d ε 0,m (x)) • n m + [(E m -E f ) : (u (0) (x))] • n f on Γ ε , (2.57d) d ε 0,m (x) = 0 on Γ. (2.57e)
This time, we note the emergence of two terms

[(E m -E f ) : (u 0 )(x)] • n f on Γ ε , which is similar to the term (α m -α f )∂ n u (0) (x)
m,ε = Ω ε m ε , D 0,k (y) = d 0,k (x) (k ∈ { f , m}) is solution to the following scaled problem : div y σ f (D ε 0, f )(y) = -ε 2 div x (E f : (u (0) ))(εy)) in Ω 1 f , (2.58a) div y σ m (D ε 0,m )(y) = 0 in Ω m,ε ,
(2.58b)

D ε 0, f (y) = D ε 0,m (y) on Γ 1 ,
(2.58c)

σ f (D ε 0, f )(y) • n f = -σ m (D ε 0,m )(y) • n m +ε[(E m -E f ) : x (u (0) (εy))] • n f on Γ 1 , (2.58d) D ε 0,m (y) = 0 on Γ ε . (2.58e)
Now, as ε is considered to be small compared to the size of the domain, we approximate Problem (2.58) by the following problem defined on the dilated domain Ω ∞ = lim ε→0 Ω m,ε :

div y σ 1 (D ε 0, f )(y) = -ε 2 div x (E f : (u (0) ))(εy)) in Ω 1 f , (2.59a) div y σ m (D ε 0,m )(y) = 0 in Ω ∞ ,
(2.59b)

D ε 0, f (y) = D ε 0,m (y) on Γ 1 , (2.59c) σ f (D ε 0, f )(y) • n f = -σ m (D ε 0,m )(y) • n m +ε[(E m -E f ) : x (u (0) (εy))] • n f on Γ 1 , (2.59d) D ε 0,m (y) → 0 when |y| → ∞. (2.59e)
Let us then introduce two problems which will allow us to construct the different orders of the asymptotic expansion for the contribution of the inclusion. Both of these problems are defined on the dilated domain Ω ∞ . For each order p, we consider the functions V (p) and Z (p) solutions to :

div y σ f (V (p) f )(y) = 0 in Ω 1 f , (2.60a) div y σ m (V (p) m )(y) = 0 in Ω ∞ ,
(2.60b)

V (p) f (y) = V (p) m (y) on Γ 1 , (2.60c) σ f (V (p) f )(y) • n f = -σ m (V (p) m )(y) • n m + ϕ (p) (y) on Γ 1 , (2.60d) V (p) m (y) → 0 when |y| → ∞, (2.60e)
and

div y σ f (Z (p) f )(y) = κ (p) (y) in Ω 1 f , (2.61a) div y σ m (Z (p) m )(y) = 0 in Ω ∞ ,
(2.61b)

Z (p) f (y) = Z (p) m (y) on Γ 1 , (2.61c) σ f (Z (p) f )(y) • n f = -σ m (Z (p) m )(y) • n m on Γ 1 ,
(2.61d)

Z (p) m (y) → 0 when |y| → ∞.
(2.61e)

The functions V (p) and Z (p) are introduced to offset terms appearing on the boundary Γ ε and in the domain Ω ε f where the data φ (p) (y) and κ (p) (y) will be adapted to each order p.

The existence and uniqueness of the functions V (p) and Z (p) are established later in Section 2.3.2. From the regularity assumptions on the domain we can write Taylor expansions of x (u (0) )(εy) for y ∈ Γ 1 and -div x (E f : (u (0) ))(εy) for y ∈ Ω 1 f :

x (u (0) (εy)) = K i =0 ε i ω (i ) 0 (y) + O(ε K+1 ), (2.62)
and

-div x (E f : (u (0) ))(εy) = K i =0 ε i γ (i ) (y) + O(ε K+1 ).
(2.63)

First-order approximation

The first-order approximation is constructed by compensating the first term in the expansion of [(E m -E f ) : (u (0) )(εy)] • n f and the first one in the expansion of -div x (E f : (u (0) ))(εy). With V (0) solving Problem (2.60) for ϕ (0)

(y) = [(E m -E f ) : ω (0) 0 (y)]
• n f and Z (0) solving Problem (2.61) for κ (0) (y) = γ (0) (y) we update the residual :

d ε 1,k (x) = u ε k (x) -u (0) (x) -εV (0) k ( x ε ) -ε 2 Z (0) k ( x ε ), k ∈ { f , m}.
(2.64)

The update of the difference allows us to obtain the first-order approximation :

u [1] ε,k (x) = u (0) (x) + εV (0) k ( x ε ) + ε 2 Z (0) k ( x ε ), k ∈ { f , m}. (2.65) Then, for k ∈ { f , m}, d ε 1,k is solution to : div x σ f (d ε 1, f )(x) = -div x (E f : (u (0) ))(εy) -γ (0) 0 (y) in Ω ε f , (2.66a) div x σ m (d ε 1,m )(x) = 0 in Ω ε m , (2.66b) d ε 1, f (x) = d ε 1,m (x) on Γ ε , (2.66c) σ f (d ε 1, f )(x) • n f = -σ m (d ε 1,m )(x) • n m +(E m -E f ) : [ x (u (0) (εy)) -ω (0) 0 (y)] • n f on Γ ε , (2.66d) d ε 1,m (x) = -εV (0) m ( x ε ) -ε 2 Z (0) m ( x ε ) on Γ.
(2.66e)

Second-order approximation

We have built a first-order approximation of the solution with u (0) , V (0) and Z (0) independent of ε. The next step is to build the second-order approximation by compen-

sating (E m -E f ) : [ (u (0) (εy)) -ω (0) 0 (y)] • n f on Γ ε , -div x (E f : (u (0) ))(εy) -γ (0) 0 (y) in Ω ε f and -εV (0) m ( x ε ) on Γ.
For that we introduce equivalent functions that play the same role as those introduced before but this time on Γ :

div x σ m (v (q) )(x) = 0 in Ω, (2.67a) v (q) (x) = ψ (q) (x) on Γ.
(2.67b)

We need to expand V (0) m at infinity :

V (0) m (y) = K i =1 a (i ) 0 (θ) |y| i + O 1 |y| K+1 = K i =1 ε i a (i ) 0 (θ) |x| i + O(ε K+1 ),
x ∈ Γ.

(2.68)

The previous expansion is proven in Section 2.3.3 in the isotropic case. So, with V (1) solution of Problem (2.60) for ϕ (1)

(y) = [(E m -E f ) : ω (1) 0 (y)] • n f , Z (1)
solution Problem (2.61) for κ (1) (y) = γ (1) (y) and v (2) solution of Problem (2.67) for ψ (2) (x) = -a (1) 0 (θ) |x| we update our residual as follows :

d ε 2,k (x) = d ε 1,k (x) -ε 2 V (1) k ( x ε ) -ε 2 v (2) (x) -ε 3 Z (1) k ( x ε
).

(2.69)

The update of the difference allows us to obtain the second-order approximation :

u [2] ε,k (x) = u (0) (x) + εV (0) k ( x ε ) + ε 2 Z (0) k ( x ε ) + ε 2 V (1) k ( x ε ) + ε 2 v (2) (x) + ε 3 Z (1) k ( x ε ), k ∈ { f , m}.
(2.70) Then, for k ∈ { f , m}, d ε 2,k solves the following problem :

div x σ f (d ε 2, f )(x) = -div x (E f : (u (0) ))(εy) -γ (0) 0 (y) -εγ (1) 0 (y) in Ω ε f , (2.71a) div x σ m (d ε 2,m )(x) = 0 in Ω ε m , (2.71b) d ε 2, f (x) = d ε 2,m (x) on Γ ε , (2.71c) σ f (d ε 2, f )(x) • n f = -σ m (d ε 2,m )(x) • n m +(E m -E f ) : [ (u (0) (εy))) -ω (0) 0 (y) -εω (1) 0 (y)] • n f +ε 2 (E m -E f ) : [ (v (2) (εy))] • n f on Γ ε , (2.71d) d ε 2,m (x) = ε(ε a (1) 0 (θ) |x| -V (0)( x ε )) m -ε 2 (Z (0) m ( x ε ) + V (1) m ( x ε )) -ε 3 Z (1) m ( x ε ) on Γ.
(2.71e)

Third-order approximation

Following the same approach, we expand as done before (v (2) )(εy) as follows :

(v (2) )(εy) = K i =0 ε (i ) ω (i ) 2 (y) + O(ε K+1 ).
(2.72) Now, we need to expand V (1) m + Z (0) m at infinity since both functions have a factor of ε 2 :

V (1) m (y) + Z (0) m (y) = K i =1 a (i ) 1 (θ) |y| (i ) + O 1 |y| K+1 = K i =1 ε i a (i ) 1 (θ) |x| i + O(ε K+1 ), x ∈ Γ.
(2.73)

The previous expansion is proven in Section 2.3.

3 in the isotropic case. With V (2) solution of Problem (2.60) for ϕ (2) (y) = (E m -E f ) : (ω (2) 0 (y) + ω (0) 2 (y)) • n f , Z (2)
solving Problem (2.61) for κ (2) (y) = γ (2) (y) and v (3) solution of Problem (2.67) for ψ (3) = -a (1) 1 (θ)

|x| -

a (2) 0 (θ)
|x| 2 we write the third-order update :

d ε 3,k (x) = d ε 2,k (x) -ε 3 v (3) (x) -ε 3 V (2) k ( x ε ) -ε 4 Z (2) k ( x ε ), k ∈ { f , m}.
(2.74)

The update of the difference allows us to obtain the third-order approximation :

u [3] ε,k (x) = u (0) (x) + εV (0) k ( x ε ) + ε 2 Z (0) k ( x ε ) + ε 2 V (1) k ( x ε ) + ε 2 v (2) (x) + ε 3 Z (1) k ( x ε ) + ε 3 v (3) (x) + ε 3 V (2) k ( x ε ) + ε 4 Z (2) k ( x ε ), k ∈ { f , m}. (2.75) CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA
Then, for k ∈ { f , m}, d ε 3,k solves the following problem :

div x σ f (d ε 3, f )(x) = -div x (E f : (u (0) ))(εy) -γ (0) 0 (y) -εγ (1) 0 (y) -ε 2 γ (2) 0 (y) in Ω ε f , (2.76a) div x σ 2 (d ε 3,m )(x) = 0 in Ω ε m , (2.76b) d ε 3, f (x) = d ε 3,m (x) on Γ ε , (2.76c) σ f (d ε 3, f )(x) • n f = -σ m (d ε 3,m )(x) • n m +(E m -E f ) : [ε(u (0) )(εy) -ω (0) 0 (y) -εω (1) 0 (y) -ε 2 ω (2) 0 (y)] • n f +ε 2 (E m -E f ) : [ (v (2) (εy)) -ω (0) 2 (y)] • n f +ε 3 (E m -E f ) : [ (v (3) (εy))] • n f on Γ ε , (2.76d) d ε 3,m (x) = ε(ε 2 a (2) 0 (θ) |x| 2 + ε a (1) 0 (θ) |x| -V (0) m ( x ε )) +ε 2 (ε a (1) 1 (θ) |x| -V (1) m ( x ε ) -Z (0) m ( x ε )) -ε 3 (V (2) m ( x ε ) + Z (1) m ( x ε )) -ε 4 Z (2) m ( x ε )
on Γ.

(2.76e)

Fourth-order approximation

As done for the Poisson problem, we go there to the fourth-order to deduce the general form of the function ψ. Now we expand (v (3) )(εy) :

(v (3) )(εy) = K i =0 ε i ω (i ) 3 (y) + O(ε K+1 ).
(2.77)

We need to expand too V (2) m + Z (1) m at infinity since both functions have a factor of ε 3 :

V (2) m (y) + Z (1) m (y) = K i =1 a (i ) 2 (θ) |y| i + O 1 |y| K+1 = K i =1 ε i a (i ) 2 (θ) |x| i + O(ε K+1 ), x ∈ Γ.
(2.78)

The previous expansion is proven in Section 2.3.3 in the isotropic case.

Finally, with V (3) solution of Problem (2.60) for ϕ (3) 3) solving Problem (2.61) for κ (3) (y) = γ (3) 0 (y) and v (3) solution of Problem (2.67) for ψ (3) = -a (1) 2 (θ)

(y) = (E m -E f ) : (ω (3) 0 (y) + ω (1) 2 (y) + ω (0) 3 (y)) • n f , Z (
|x| -

a (2) 1 (θ)
|x| 2 -a (3) 0 (θ) |x| 3 we write the fourth-order update :

d ε 4,k (x) = d ε 3,k (x) -ε 4 v (4) (x) -ε 4 V (3) k ( x ε ) -ε 5 Z (3) k ( x ε ), k ∈ { f , m}. (2.79) CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA
We obtain the following equations for

d ε 4,k (x), k ∈ { f , m} : div σ f (d ε 4, f )(x) = -div x (E f : (u (0) ))(εy) -γ (0) 0 (y) -εγ (1) 0 (y) -ε 2 γ (2) 0 (y) -ε 3 γ (3) 0 (y) in Ω ε f , (2.80a) div σ m (d ε 4,m )(x) = 0 in Ω ε m ,(2.80b) d ε 4, f (x) = d ε 4,m (x) on Γ ε , (2.80c) σ f (d ε 4, f )(x) • n f = -σ m (d ε 4,m )(x) • n m +(E m -E f ):[( (u (0) (εy)) -ω (0) 0 (y) -εω (1) 0 (y) -ε 2 ω (2) 0 (y) -ε 3 ω (3) 0 (y)] • n f +ε 2 (E m -E f ) : [( (v (2) (εy)) -ω (0) 2 (y)) -εω (1) 2 (y))] • n f +ε 3 (E m -E f ) : [ (v (3) (εy)) -ω (0) 3 (y)] • n f +ε 4 (E m -E f ) : [ (v (4) (εy))] • n f on Γ ε , (2.80d) d ε 4,m (x) = ε(ε 3 a (3) 0 (θ) |x| 3 + ε 2 a (2) 0 (θ) |x| 2 + ε a (1) 0 (θ) |x| -V (0) m ( x ε )) +ε 2 (ε 2 a (2) 1 (θ) |x| 2 + ε a (1) 1 (θ) |x| -V (1) m ( x ε ) -Z (0) m ( x ε )) +ε 3 (ε a (1) 2 (θ) |x| -V (2) m ( x ε ) -Z (1) m ( x ε )) -ε 3 (V (3) m ( x ε ) + Z (2) m ( x ε )) -ε 5 Z (3) m ( x ε ) on Γ. (2.80e)

The K-th order approximation

We can update our residual at the order n ≥ 3 as follows :

d ε n,k (x) = d ε n-1,k (x) -ε n v (n) k (x) -ε n V (n-1) k ( x ε ) -ε n+1 Z (n-1) k ( x ε ), k ∈ { f , m}. (2.81) With v (n) solves Problem (2.67) for ψ (n) = -p+q=n-1 p≥1,q≥0 a (p) q (θ) |x| p , V (n-1) the solution of Problem (2.60) with ϕ (n-1) (y) = (E m -E f ) : p+q=n-1 p≥1,q≥0 ω (p)
q (y) • n f , and Z (n-1) solving Problem (2.61) for κ (n-1) (y) = γ (n-1) (y). We are now able to build an approximation u [K] ε (x) to the order K ≥ 3 (k ∈ { f , m}) :

u [K] ε,k (x) = u (0) (x) + K i =2 ε i v (i ) (x) + K i =1 ε i V (i -1) k ( x ε ) + K i =1 ε i +1 Z (i -1) k ( x ε ).
(2.82)

We obtain the following equations for

d ε K,k (k ∈ { f , m}) : div σ f (d ε K, f )(x) = O H -1 (Ω;R 2 ) (ε K+1 ) in Ω ε f , (2.83a) div σ m (d ε K,m )(x) = 0 in Ω ε m , (2.83b) d ε K, f (x) = d ε K,m (x) on Γ ε , (2.83c) σ f (d ε K, f )(x) • n f = -σ m (d ε K,m )(x) • n m + O(ε K ) on Γ ε , (2.83d) d ε K,m (x) = O H 1 2 (Γ,R 2 ) (ε K+1 ) on Γ.
(2.83e)

Existence and uniqueness in unbounded domains for the linear elasticity problem

A weak formulation of Problem (2.60) writes :

Ω 1 f (V (p) f ) : E f : (ω) dx + Ω ∞ (V (p) m ) : E m : (ω) dx = Γ 1 ϕ • w dσ,
for all w ∈ V.

(2.84) Introducing the space V = v ∈ (W 1,2 0,0 ) 2 ; Γ R 0 v dσ = 0 we can see that the bilinear form

Ω 1 f (V (p) f ) : E f : (ω) dx + Ω ∞ (V (p)
m ) : E m : (ω) dx is coercive on V. In addition, this bilinear form and the linear form Γ 1 ϕ • w dσ are continuous on V. Then, the existence and uniqueness of the function V (p) in V solution to (2.84) follow from the Lax-Milgram theorem. The existence and uniqueness of Z (p) are proven in the same way.

Expansion of the functions V (p) m and Z (p) m at infinity

We restrict here to the isotropic case where E m : = λ trace( )I + 2µ , for λ, µ the Lamé coefficients. In that case, according to BREZIS [1996] the functions m which takes the following form :

ψ (s) 1 (θ) = sin((s + 1)θ) cos((s + 1)θ) , (2.85)
ψ (s) 2 (θ) = cos((s + 1)θ) -sin((s + 1)θ) , (2.86)
ψ (s) 3 (θ) = (1 -s + 2ρ) sin((s -1)θ) -(1 + s + 2ρ) cos((s -1)θ) , ( 2
V (p) m = Z (p) m = s∈N * r -s a s ψ (-s) 1 + b s ψ (-s) 2 + c s ψ (-s) 3 + d s ψ (-s) 4 .
(2.89)

Error estimate for the linear elasticity problem

In this section we establish an error estimate for the studied problem by increasing the H 1 (Ω; R 2 ) norm by powers of ε (characteristic length of the inclusion). First, we begin by considering the following problem for a given

f ε ∈ H -1 2 (Γ ε ; R 2 ), h ε ∈ H 1 2 (Γ; R 2 ) and s ε ∈ H -1 (Ω ε f ; R 2 ) ( f ε (x) = g ε (y)) : div σ f (z f )(x) = -s ε (x) in Ω ε f , (2.90a) div σ m (z m )(x) = 0 in Ω ε m , (2.90b) z f (x) = z m (x) on Γ ε , (2.90c) σ f (z f )(x) • n f = -σ m (z m )(x) • n m + f ε (x) on Γ ε , (2.90d) z m (x) = h ε (x)
on Γ.

(2.90e)

CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA

Its weak form can be written as : Find z ∈ H 1 (Ω; R 2 ) with z = h on Γ such that :

Ω ε f (z) : E f : (ω) dx + Ω ε m (z) : E m : (ω) dx = Γ ε f ε • ω dσ + Ω ε f s ε • ω dx, (2.91)
for all ω ∈ H 1 (Ω; R 2 ) such that ω = 0 on Γ.

Lemma 2.3.1. There exists a unique solution z to the weak form (2.91) of Problem (2.90)

satisfying z H 1 (Ω) ≤ C h ε H 1 2 (Γ;R 2 ) + g ε L 2 (Γ 1 ;R 2 ) + s ε H -1 (Ω ε f ;R 2 ) , where C > 0 is a constant independent of ε.
Proof. The surjectivity of the trace operator guarantees the existence of a continuous lifting operator H ∈ H 1 (Ω; R 2 ) of h ε (see [START_REF] Fortin | Weighted Sobolev spaces for Laplace's equation in R n[END_REF]). Then, z = z -H solves the following equation for all ω ∈ H 1 (Ω; R 2 ) such that ω = 0 on Γ :

Ω ε f (z) : E f : (ω) dx + Ω ε m (z) : E m : (ω) dx + Ω ε f (H) : E f : (ω) dx + Ω ε m (H) : E m : (ω) dx - Γ ε f ε • ω dσ = Ω ε f s ε • ω dx. (2.92) The bilinear form a(u, v ) = Ω ε f (u) : E f : (v ) dx + Ω ε m (u) : E m : (v ) dx associated to the weak formulation (2.92) is coercive, i.e. α v 2 H 1 (Ω;R 2 )
≤ a(v , v ), see e.g. ERN et GUERMOND [2004]; FORTIN et GARON [2011]. Then, the Lax-Milgram theorem guarantees the existence and uniqueness of z solution of (2.92) and hence existence and uniqueness of solution to (2.92) follow. Using the coercivity of the bilinear form we obtain :

α z 2 H 1 (Ω;R 2 ) ≤ Γ ε f ε .z dσ - Ω σ(H) : (z) dx + Ω ε f s ε .z dx, (2.93)
hence, using the estimation (A.17) demonstrated in the appendix A.2 we obtain

α z 2 H 1 (Ω;R 2 ) ≤ C g ε L 2 (Γ 1 ;R 2 ) z H 1 2 (Γ ε ;R 2 ) + (H) L 2 (Ω;R 2 ) (z) L 2 (Ω;R 2 ) + s ε H -1 (Ω;R 2 ) z H 1 (Ω;R 2 ) ,
(2.94) and therefore,

α z 2 H 1 (Ω;R 2 ) ≤ C g ε L 2 (Γ 1 ;R 2 ) + H H 1 (Ω;R 2 ) + s ε H -1 (Ω;R 2 ) z H 1 (Ω;R 2 ) .
(2.95)

The continuity of the lifting operator gives us that

H H 1 (Ω;R 2 ) ≤ C h ε H 1 2 (Γ;R 2 )
, so we can write :

α z H 1 (Ω;R 2 ) ≤ C g ε L 2 (Γ 1 ;R 2 ) + h ε H 1 2 (Γ;R 2 ) + s ε H -1 (Ω;R 2 ) .
(2.96)

We can now state the following result :

Theorem 2.3.2. In the isotropic case E m : = λ trace( )I + 2µ , there exists a constant C > 0 independent of ε such that :

d ε K H 1 (Ω;R 2 ) ≤ Cε K , for all K ∈ N. CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIA Proof. With z = d ε K we have s ε ∈ O H -1 (Ω;R 2 ) (ε K+1 ), g ε ∈ O L 2 (Γ 1 ,R 2 ) (ε K ) and h ε ∈ O H 1 2 (Γ,R 2 ) (ε K+1 )
using Lemma 2.3.1 we obtain the announced error estimate.

Remark 2.3.3. The limitation of the error estimate to the isotropic case is only due to the limitation to this case of the expansion at infinity given in Section 2.3.3. An extension to some non-isotropic cases may require an adapted result of expansion at infinity of functions V

Conclusion and perspectives

The purpose of this work was to study the interaction between an inclusion of a small size compared to the size of the elastic body in which it operates. We were able to build an asymptotic expansion starting from the field without inclusion and correcting by additional functions that model the presence of the inclusion and its influence on the mechanical fields (asymptotic expansion at an arbitrary order K given by u [K] ε (x) in (2.26) for the Poisson problem and by the displacement u [K] ε (x) in (2.82) for linear elasticity). The model had been validated by comparing the obtained result given by (2.55) with the one obtained by solving the problem using a separation of variable method (2.51) (for the Poisson problem and with a circular shape of the inclusion). Error estimates had been established in Sections 2.2.4 and 2.3.4. The main result for the applications is that the functions v (i ) and V (i ) , solutions to some problems do not depend on ε and depend linearly on u (0) only by the expansion of its gradient at x = 0. For a given shape of the inclusion, this make it possible to pre-compute once for all a basis for each v (i ) and V (i ) . Natural perspectives of this work is the design of a numerical strategy first to compute a few terms of the expansion obtained and then to compute the influence of an inclusion in an arbitrary shaped body and the influence of several inclusions. Extension to large strain elasticity and to three-dimensional problems are also of great interest. Finally, the extension to inhomogeneous and anisotropic cases should be done without major difficulties.

Introduction

In the previous chapter we have studied the influence of the presence of an inclusion of very small size compared to the dimensions of the elastic medium. In this work, we did not use any homogenization or averaging technique which requires a periodic distribution of inclusions density and the used techniques, techniques of scaling and multi-scale asymptotic expansions, led us to obtain corrective terms modeling the influence of an inclusion on the elastic domain (asymptotic expansion of the displacement field around the characteristic parameter of the inclusion ε). Some of these correctors solve problems posed on unbounded domains which leads us to propose an adequate numerical strategy to approximate these correctors and thus validate the modeling built beforehand. The second type of correctors solve problems on bounded domains and hence does not require a particular study. Thus, the major difficulty lies in the approximation of the functions posed on infinite domains (We shall content ourselves in this work with the first-order corrector).

The outline of this chapter is as follows : We begin by presenting in Section 3.2 the goals and the process to follow, we then proceed in Section 3.4 to the choice of the appropriate numerical method. The validation of the method is done in Section 3.5 by comparing the obtained results with those given analytically for test cases (using the separation of variables method) and we propose an iterative numerical strategy taking into account the presence of several inclusions into the same medium. Finally, we propose in Section 3.5.3 an adaptation for the linear elasticity problem.

Main motivation and process to follow

As mentioned previously in the introduction, our subject of interest is to numerically approximate a mechanical problem whose configuration is in the form of an arrangement of an elastic domain containing an inclusion of different stiffness and of size extremely smaller than the dimensions thereof. The major problem we want to circumvent is the inability to take into account an inclusion by the classical finite-element method, due to the smallness of its size compared to the size of the whole domain. Thus, unless meshing in an extremely fine manner around the zone containing the inclusion, which can quickly become prohibitive in terms of computation means and time, the obtained solutions are not representative of the problem that they are desired to approximate. To avoid this we call upon the previously introduced formalism and we are content here to seek to approximate the first-order corrector that appears in the first-order asymptotic expansion (1.9) (the higher-order correctors can be approximated by the same technique). Let us first mention the first gain which is that our correctors are posed on domains independent of ε which saves us from meshing in an excessively fine way in order to be able to capture the inclusion. The second point to mention is that we have two types of correctors : Some are placed on bounded domains and we use the classical finite-element method to obtain the approximation, the second kind is posed on unbounded domains and we need an appropriate method to approximate them. The ultimate contribution lies in the fact that our first-order corrector does not depend on ε and depends only on u (0) (the solution without inclusion) by the expansion of its gradient at x = 0 of the blow they will be computed only once (thanks to the linearity of our equations) which gives us a gain in the management of several inclusions interacting in the same domain by pre-computing a basis once for all our correctors for a given shape of inclusion.

As already mentioned, the final goal is to present an iterative strategy of taking into account several inclusions starting from an arbitrarily shaped matrix and correcting by Neumann conditions on an artificial interface representing the limit of the zones of influence of the inclusions (Γ N i ). These Neumann conditions are obtained from a pre-computed base of our first-order corrector (see Fig. 3.1).

Γ D Γ N 2 Γ N 4 Γ N 1 Γ N 3 Ω 1 int Ω 2 int Ω 3 int Ω 4 int Ω ext

A numerical method for the problem in infinite domain

The approximation of solution to partial differential equations in infinite domains can be done with several methods having different characteristics concerning accuracy, sophistication and flexibility. Most of these methods rely either on an integral representation of the solution and the use of boundary elements (see, e.g., BENDALI et collab. [1999]; BENDALI et SOUILAH [1994]; BOUBENDIR et collab. [2008]; COLTON et KRESS [1983]; GI-ROIRE [1987]; GIROIRE et NÉDÉLEC [1978]; LEROUX [1977]; NÉDÉLEC [1976,1998]; ROKH-LIN [1983]) ; -On replacing the unbounded domain by a sufficiently large bounded domain enclosed by a Perfectly Matched Layer (PML) (see BÉRENGER [1994,1996]), or on the boundary of which an artificial boundary condition is prescribed ; -Or on a polar expansion of the solution like in spectral methods (see, e.g., CANUTO et collab. [1985]) or in infinite-elements methods.

Because of its flexibility we choose here to adapt a method proposed by BOULME-ZAOUD [2005] which consists in mapping the unbounded domain into a bounded one called Inverted Finite-Element Method (IFEM). The classical finite-element method is then used on the two domains (the bounded domain and the mapped unbounded domain into a bounded one by means of a polygonal inversion). In the version we propose, the two domain solutions are linked by a Lagrange multiplier. The method is of an arbitrary degree and is exactly conforming. An adequate weighted Sobolev space is used to describe the decay of our functions at infinity. Finally, this method does not use any artificial boundary or conditions.

CHAPITRE 3. A NUMERICAL STRATEGY FOR THE ACCOUNT OF SMALL INCLUSIONS IN ELASTIC MEDIA

Our goal is to use an inverted finite-element method to approximate solutions of (2.5) completed with a condition prescribed on a fictitious arbitrary contour Γ R 0 V = 0. V m is also required to satisfy the following condition :

V m (y) = o(log |y|) 1 2
when |y| → ∞.

(3.1)

In the next section, we prepare the appropriate framework to formulate the IFEM method introduced by Boulmezaoud in BOULMEZAOUD [2005] and adapted in BHOWMIK et collab.

[2016] for the case of a bidimensional perforated unbounded domain with a homogeneous Dirichlet condition.

The adequate logarithmic weighted Sobolev space and some properties

Let us remember that our aim is to approximate Problem (2.5) which is posed on the dilated geometry. For that, we recall the spaces

W 1,2 0,0 = u ∈ D (R 2 ), ρ -1 (lg r ) -1 u ∈ L 2 (R 2 ),∇ ∇ ∇u ∈ (L 2 (R 2 )) 2
, and

V = v ∈ W 1,2 0,0 ; Γ R 0 vdσ = 0 , introduced previously. It is well know that D(R 2 ) is dense in W 1,2 0,0 (see AMROUCHE et collab. [1994]). Lemma 3.3.1. Let V k ∈ W 1,2 0,0 (R 2 ) (k ∈ { f , m}). Then, V m (r, θ) L 2 (S 1 ) log r = 0, when r → ∞, (3.2)
where S 1 is the unit circle of R 2 .

Proof. Let r > 1. Then,

R 2 \B r |ϕ(y)| 2 (|y| 2 )(log(|y|) 2 dy = S 1 +∞ r |ϕ(r, θ)| 2 r (log r ) 2 dr dθ = S 1 |ϕ(r, θ)| 2 log r dθ + 2 S 1 +∞ r ϕ(r, θ) ∂ϕ ∂r (r, θ) log r dr dθ Cauchy-Shwartz inequality gives ϕ(r, .) 2 L 2 (S 1 ) log r ≤ 2 ϕ |y| log(|y|) 2 L 2 (R 2 \B r ) + ∇ϕ 2 L 2 (R 2 \B r ) ≤ 2 ϕ 2 W 1,2 0,0 (R 2 \B r )
.

By density, the latter inequality remains valid for ϕ ∈ W 1,2 0,0 (R 2 ). Taking the limit when r → ∞ ends the proof.

It is also demonstrated in Section 2.2.3 of Chapter 2 that the functions V m can be expanded at ∞ in the following form :

V (0) m (y) = K i =1 a (i ) 0 (θ) |y| i + O 1 |y| K+1 = K i =1 ε i a (i ) 0 (θ) |x| i + O(ε K+1 ),
x ∈ Γ.

Inverted finite-element method : Details and convergence

This section is devoted to the discretization of the two bidimensional elliptic Problems (2.5) and (2.60) by means of the inverted finite-element method.

T ∞ h a 0 a 1 a 2 S FIGURE 3.
2 -An illustration of the altitude vector h, the finite simplex and the supporting hyperplane associated to a 2D infinite simplex.

First, let us recall the definition of infinite simplices or infinite (triangles) according to [START_REF] Boulmezaoud | Inverted finite elements : A new method for solving elliptic problems in unbounded domains[END_REF] adapted to the bidimensional case. Given three points a 0 , a 1 and a 2 of R 2 , the infinite triangle T ∞ is defined as follows :

T ∞ (a 0 | a 1 , a 2 ) = {λ 0 a 0 + λ 1 a 1 + λ 2 a 2 | λ 0 ≤ 0, λ 1 ≥ 0, λ 2 ≥ 0, λ 0 + λ 1 + λ 2 = 1}.
a 0 is called the fictitious vertex of T ∞ (a 0 | a 1 , a 2 ), while a 1 and a 2 called its real vertices (see Fig. 3.2). We associate to this infinite triangle T ∞ (a 0 | a 1 , a 2 ) the triangle S whose vertices are a 0 , a 1 and a 2 :

S(a 0 , a 1 , a 2 ) = {λ 0 a 0 + λ 1 a 1 + λ 2 a 2 | λ 0 ≥ 0, λ 1 ≥ 0, λ 2 ≥ 0, λ 0 + λ 1 + λ 2 = 1}.
The altitude vector of T ∞ (a 0 | a 1 , a 2 ) is defined as h = Πa 0a 0 where Πa 0 is the orthogonal projection of the fictitious vertex a 0 on the line (a 1 a 2 ). It follows that |h| is the distance between a 0 and the line (a 1 a 2 ).

The first step in deploying inverted finite elements consists of partitioning the whole R 2 (see Fig. 1.2) as :

R 2 = Ω b ∪ Ω ∞ with Ω b ∩ Ω ∞ = , (3.3) with Ω b a bounded open subset of R 2 , Ω ∞ an unbounded open subset of R 2 .
We assume in addition that the unbounded component Ω ∞ can be partitioned as the union of a finite number of infinite simplices as follows :

Ω ∞ = ∪ M i =1 T ∞ i , (3.4)
where T 1 , . . . , T M are infinite triangles having a common fictitious vertex a 0 and satisfying :

∀i = j , int(T ∞ i ) ∩ int(T ∞ j ) = . (3.5)
Without loss of generality, we suppose that

a 0 = 0. (3.6)
We also assume that the intersection of T ∞ i and T ∞ j , i = j , is either a common infinite edge or empty. We denote by h 1 , . . . , h M the altitude vectors of T ∞ 1 , . . . , T ∞ M , and by S 1 , . . . , S M the associated triangles. We set

Ω * = int(∪ M i =1 S i )\{a 0 } (3.7)
We may observe that Ω * is a bounded polygonal open subset of R 2 . It is a fictitious domain which must be distinguished from Ω b . As we shall see, Ω * will play a prominent role in the development of the inverted finite-element method. We define the polygonal radius r (•) :

y ∈ R 2 → R + as follows : for 1 ≤ i ≤ M and y ∈ T ∞ i ∪ S i r (y) = y • h i |h i | 2 . (3.8) Denoting (a (i ) 0 = 0, a (i ) 1 , a (i )
2 ) the vertices of S i and by λ (i ) 0 , λ (i ) 1 , λ (i ) 2 their corresponding barycentric coordinates. Then,

r (y) = λ (i ) 1 (y) + λ (i ) 2 (y) = 1 -λ (i ) 0 (y) for y ∈ T ∞ i ∪ S i . (3.9) Indeed, for y = λ (i ) 0 (y)a (i ) 0 + λ (i ) 1 (y)a (i ) 1 + λ (i ) 2 (y)a (i ) 2 ∈ T ∞ i ∪ S i , we have y • h i = λ (i ) 1 (y) • h i + λ (i ) 2 (y) • h i = (λ (i ) 1 (y) + λ (i ) 2 (y))|h i | 2 .
From (3.9), it becomes clear that r is continuous between two neighboring triangles T ∞ i and T ∞ j (or S i and S j ). In other words, r is continuous over R 2 . Moreover, there exist two constants c 1 > 0 and c 2 > 0 such that

∀y ∈ R 2 , c 1 |y| ≤ r (y) ≤ c 2 |y|. (3.10)
Consider also the polygonal inversion Φ : Ω ∞ → Ω * defined as follows

∀y ∈ Ω ∞ , Φ(y) = y r (y) 2 . (3.11)
It is quite clear that Φ is a bijective transformation from Ω ∞ into Ω * . One can easily prove the following :

-Φ is continuous ; (3.12) -for all y ∈ Ω ∞ : r (Φ(y)) • r (y) = 1; -Φ(y) = y if Ω ∞ ∩ Ω * = Ω ∞ ∩ Ω b ; -Φ(T ∞ i ) = S i \{0}, i = 1, . . . , M;
-Φ is a bijection from Ω ∞ to Ω * and

∀y * ∈ Ω * Φ -1 (y * ) = y * r (y * ) 2 .

CHAPITRE 3. A NUMERICAL STRATEGY FOR THE ACCOUNT OF SMALL INCLUSIONS IN ELASTIC MEDIA

For convenience, we also define the piecewise constant altitude vector :

h(x) = h i for x ∈ T ∞ i ∪ S i , i = 1, . . . , M. (3.13)
In order to clarify these notions, let us consider a particular choice of the decompositions (3.3) and (3.4). Define the box

Q = (y 1 , y 2 ) ∈ R 2 | |y 1 | ≤ R 1 and |y 2 | ≤ R 2 = [-R 1 , R 1 ] × [-R 2 , R 2 ],
where R 1 > 0 and R 2 > 0 are two real numbers chosen such that

Ω 1 f ⊂ Q. Let Ω b = R 2 ∩ int(Q), Ω ∞ = R 2 \Q = {(x 1 , x 2 ) ∈ R 2 | |y 1 | > R 1 or |y 2 | > R 2 }. (3.14) We can write Ω ∞ = T ∞ 1 ∪ T ∞ 2 ∪ T ∞ 3 ∪ T ∞ 4 , (3.15) with T ∞ 1 = (y 1 , y 2 ) ∈ R 2 | y 1 ≥ R 1 and y 1 ≥ R 1 R 2 |y 2 | , T ∞ 2 = (y 1 , y 2 ) ∈ R 2 | y 2 ≥ R 2 and y 2 ≥ R 2 R 1 |y 1 | , T ∞ 3 = (y 1 , y 2 ) ∈ R 2 | y 1 ≤ -R 1 and y 1 ≤ - R 1 R 2 |y 2 | , T ∞ 4 = (y 1 , y 2 ) ∈ R 2 | y 2 ≤ -R 2 and y 2 ≤ - R 2 R 1 |y 1 | .
It follows that

Ω * = (y 1 , y 2 ) ∈ R 2 | |y 1 | < R 1 and |y 2 | < R 2 }\{(0, 0) .
Triangles T ∞ 1 , . . . , T ∞ 4 have the origin as a common fictitious vertex. The corresponding altitude vectors are

h 1 = R 1 e 1 , h 2 = R 2 e 2 , h 3 = -R 1 e 1 , h 4 = -R 2 e 2
with e 1 = (1, 0) and e 2 = (0, 1) are the unit vectors of the canonical basis of R 2 . We find

r (x) = max |y 1 | R 1 , |y 2 | R 2 .
The polygonal inversion is 2

Φ(y 1 , y 2 ) =      R 2 1 1 y 1 , y 2 y 2 1 if | y 1 R 1 | ≥ | x 2 R 2 |, R 2 2 y 1 y 2 2 , 1 y 2 if | y 1 R 1 | ≤ | y 2 R 2 |.
Ω 1 f Γ 1 Ω ∞ T 1 T 2 T 3 T 4 S 2 S 3 S 4 S 1 Φ Φ Φ Φ Ω b Ω * FIGURE 3.

-Inverted finite element method decomposition and correspondence

Now, we consider a set of pairs of triangulations

J (R 2 ) = (T , T * )|T ∈ J 1 (Ω b ), T * ∈ J 2 (Ω * ) such that -T is a conformal triangulation of Ω b . -T * is a conformal triangulation of Ω * .
-For each K ∈ T * , there exists i ≤ M such that K ⊂ S i (in other words, T * is a conforming union of triangulations of the subdomains S 1 , . . . , S M ).

For each pair of triangulations (T , T * ) ∈ J (Ω d ), we set

T * * = K ∈ T * |0 ∉ K ,
(in other words, triangles of T * * do not touch the common fictitious vertex a 0 = 0), and

h = h(T , T * ) = max K∈T ∪T * h K ,
where h K is the diameter of K where K is an arbitrary triangle of T ∪ T * . We also set for each

K ∈ T * d k = inf x∈K |y|,
(that is, d K denotes the distance between 0 and K ).

From the geometric point of view, we also assume that this set of triangulations satisfies the following assumption : there exists a real µ ∈ [0, 1] and four constants c 0 > 0, c * 1 > 0, c * 2 > 0 and c * 3 > 0 such that for each (T , T * ) ∈ J (Ω) we have max

K∈T ∪T * h K ρ K ≤ c 0 , (3.16) max K∈T * * h K d 1-µ K ≤ c * 1 h, (3.17) max K∈T * \T * * h K ≤ c * 2 h 1 µ , ( 3
u ∈ V such that ∀w ∈ V , a(u, w) = 〈ϕ, w〉, (3.20)
where Let k ≥ 1 be an integer. Consider the discrete space

a(u, w) = α f Ω 1 f ∇ ∇ ∇V (p) f •∇ ∇ ∇w dy + α m Ω ∞ ∇ ∇ ∇V (p) m •∇ ∇ ∇w dy. ( 3
V h = {v h ∈ C 0 (R 2 ) | v h|K ∈ P k (K), ∀K ∈ T , vh|K * ∈ P k (K * ), ∀K * ∈ T * , vh (0) = 0}.
Lemma 3.4.1. If θ > 0, then the following inclusion holds

V h → V . (3.26)
Proof. Let us observe first that for all v h ∈ V h one has

v h ∈ L 2 l oc (R 2 ), ∇v h ∈ L 2 l oc (R 2 ) 2 .
Moreover, since v h is continuous, piecewise affine and vanishing at 0, we have

∀y * ∈ Ω * , | vh (y * )| ≤ C h r (y * ).
It follows that where C and C are two constants not depending on x. Thus, using polar coordinates, and since θ > 0, we get

∀y ∈ Ω ∞ , |v h (y)| ≤ Cr (y * ) θ = Cr (y) -θ ≤ C|y| -θ (y * = Φ(y)),
R 2 |v| 2 |y| 2 (log(|y| 2 + 2) dy < ∞.
Similarly, since |∇ * vh | ≤ C 2 (∇ * vh is the gradient of vh with respect to y * ), we get

|∇v h (y)| ≤ |θ -1|r (y) -θ | v(y * )| + r (y) θ-1 |∇Φ(y)∇ * v(y * )|.
Since

∇Φ(y) = 1 r (y) 2 δ i , j -2 y j h(y) • e i r (y)|h| 2 1≤i , j ≤1
, with y = Φ(y * ).

(3.27)

In view of (3.10), we get

∇Φ(y) ≤ C|y| -2 ,
and

|∇v h (y)| ≤ Cr (y) -θ-1 ≤ C |y| -θ-1 . Hence, if θ > 0, then ∇u h ∈ L 2 (R 2 ).
In the sequel, we assume that θ > 0.

(3.28)

The discrete problem can be written as : find u h ∈ V h such that Now, before stating the following theorem we introduice the space W m α (R 2 ) as follows :

∀v h ∈ V h , a(u h , v h ) = 〈ϕ, v h 〉, ( 3
W m α (R 2 ) = u ∈ D (R 2 ), ρ |µ|-m+α D µ u ∈ L 2 (R 2 ), ∀|µ| ≤ m .
Theorem 3.4.2. The discrete problem (3.29) has one and only one solution u h ∈ V h . Moreover, suppose that u ∈ H k+1 l oc (R 2 ) and

∀|λ| ≤ k + 1, (∂ λ u)(|y|, .) L 2 (S 1 ) ≤ C |y| θ+|λ| f or |y| ≥ R, (3.30)
for two constants C > 0 and R > 0. Then, for each real number ε such that 0 < ε < mi n(1, θ), there exists a constant C ε > 0, not depending on u and h such that

u -u h W 1,2 0,0 (R 2 ) ≤ C ε h k u H k+1 (Ω b ) + h k min( µ 0 µ ,1) u W k+1 k+θ-ε (Ω ∞ ) , (3.31)
where

µ 0 = θ -ε. If in addition µ < θ k , then u -u h W 1,2 0,0 (R 2 ) ≤ C * h k , (3.32)
for some constant C * not depending on h.

CHAPITRE 3. A NUMERICAL STRATEGY FOR THE ACCOUNT OF SMALL INCLUSIONS IN ELASTIC MEDIA

Proof. Since u ∈ H k+1 l oc (R 2 ) and in view of condition (3.30),we have

u ∈ W k+1 k+θ-ε (R 2 ),
for each ε > 0.

Existence and uniqueness are standard. From Céa's Lemma, we have

u -u h W 1,2 0,0 (R 2 ) inf v h ∈V h u -v h W 1,2 0,0 (R 2 ) . (3.33) Since V h ⊂ W 1 0 (R 2 ) → W 1,2 0,0 (R 2 ), we have u -u h W 1,2 0,0 (R 2 ) inf v h ∈V h u -v h W 1 0 (R 2 ) . ( 3

.34)

Consider now the following interpolation mapping :

Π h : C 0 (R 2 ) → V h which associates to each v ∈ C 0 (R 2 ) the unique element z ∈ W h satisfying z |K = π k v for all K ∈ T , ẑ|K * = π k v for all K * ∈ T * , where -for K ∈ T or K * ∈ T * , π K denotes the local P k interpolation operator in K, -for K * ∈ T * \T * * (that is, a 0 ∈ K * ), π K * denotes the local • P k interpolation operator in K * , with • P k = {p ∈ P k | p(a 0 ) = 0}.
It is proven in BOULMEZAOUD [2005] (see Theorem 1) that

u -Π h u W 1,2 0,0 (R 2 ) h k u H k+1 (Ω b ) + h k min( µ 0 µ ,1) u W k+1 k+θ-ε (Ω ∞ )
. This ends the proof of (3.31). For getting (3.32), it suffices to choose ε > 0 such that ε < mi n(1, θkµ).

Implementation of the method

Our objective in this section is to give some details concerning the linear system arising from the discrete system (3.29). We first exhibit a basis of the space V h . Let (x i ) i ∈J (resp. (x * ,i ) i ∈J * ) be the nodes of the mesh T (resp. T * ). We define ω i , i ∈ J the function of V h satisfying ω i (x j ) = δ i , j for all j ∈ J.

We define also

ω * i , i ∈ J * the function of V h satisfying ω * i (x j ) = δ i , j for all j ∈ J * .
We propose here to detail the calculation of all these entries. Let ψ 1 and ψ 2 be two elements of V h . We have

R 2 α(y)∇ψ 1 (y) • ∇ψ 2 (y) dy = Ω b α(y)∇ψ 1 (y) • ∇ψ 2 (y) dy + Ω∞ α m ∇ψ 1 (y) • ∇ψ 2 (y) dy
The first term can be written as In what follows, we shall detail only the integral in the unbounded region Ω ∞ . We consider ∇ψ 1 and ∇ψ 2 as column vectors ; the scalar product of ∇ψ 1 and ∇ψ 2 will be written as (∇ψ 1 ) t ∇ψ 2 . We have We present in Fig. 3.4 the domain Ω b containing a circular inclusion of radius 1 (it is the scaled problem). Note that we construct a conformal mesh around the contour of the inclusion to be able to apply the rigidity α f = 100 exactly on the domain of the inclusion (red color). On the rest of the domain is applied the rigidity α m =10 (blue color). We present

Ω b α(y)∇ψ 1 (y) • ∇ψ 2 (y) dy = K∈T K α(y)∇ψ 1 (y) • ∇ψ 2 (y) dy,
ψ i (y) = ψi (y * ) = ψ i • Φ(y * ),

Computational results and discussions

The aim of this section is to present numerical results for the approximated displacement field u [1] ε,k . We begin by a convergence analysis of the first-order corrector V (0) k according to the norm of the weighted Sobolev space W 1,2 0,0 . Then, we assemble the two terms u (0) (approximated by the classical finite-element method) and εV (0) k (approximated by the inverted finite-element method). Finally, we present an iterative strategy taking into account several inclusions by adding Neumann conditions at each iteration. The Built strategy is validated by comparing it to a reference solution obtained when meshing in an extremely fine way.

Convergence analysis

For a circular shape of the contour Γ 1 of the inclusion, we solve in Chapter 2 Section 2.2.5.2 Problem (2.5) using the variable separation method. We thus obtain an exact solution which will enable us to study the convergence. Recall that the exact solution in this case is given by :

V (0) m (r, θ) = α m -α f r a α f + α m cos θ + b α f + α m sin θ .
where ϕ = (α m -α f )(a b) t • n = -90(10 100) t • n is an arbitrary source term imposed on the contour Γ 1 .

We start by plotting the approximated solution for θ = 0 The curve shows a decay towards zero at infinity with a perfect coincidence of its values with the exact solution.

We now present two convergence curves for P 1 and P 2 elements in W 1,2 0,0 . The theoretical convergence rate given by the error estimate (3.32) for P 1 elements is 1 and we obtain numerically a convergence rate of 1.0347 see Fig. 3.7 which represents a light superconvergence. For the P 2 elements the convergence rate is 1.874 see Fig. 3.8 which is in agreement with the theoretical rate 2.

Dirichlet-Neumann iterative strategy

In this Section we present an iterative method of taking into account the influence of several inclusions by correcting by Neumann conditions on a fictitious boundary representing roughly the zone of influence of the first-order corrector. The strategy is presented for a first case which is the case of a single circular inclusion of radius 0.1 placed at the center of a square matrix of dimension 20 then extended to several inclusions of same type. Finally we present some results for a non-circular inclusion.

We first numerically precalculate a basis V 10 , V 01 of our first-order corrector. V 10 is precalculated for a source term ϕ 10 = (α m -α f )(1 0) t • n and V 01 for a source term ϕ 01 = (α m -α f )(0 1) t • n. The obtained fields are then interpolated on a circular geometry of radius 10 and then stored in two files. The choice of this value is due to the abrupt decay of these correctors observed on the curves in different directions. That is to say that for a radius of 10 we recover almost all the information. Using the scaling x = εy, the zone of influence of the first-order corrector interpolated on the matrix is a circle of radius 1. Outside this influence zone, the inclusion has no longer any significant influence and the field is close to u (0) . Inside it the considered field is u (0) + εV (0) .

At each iteration, the gradient of u (0) that we denote by (a b) t is computed in x = 0 (the position where the inclusion is centered) and then V (0) is updated (V (0) = aV 10 +bV 01 ). Then, at the next iteration, a source term is added at the interface of the influence zone which corresponds to a Neumann condition in order to ensure the continuity of the normal derivative. In practice the term ε∂ n V (0) is added on the square inscribed in the circle of radius 1. We reiterate until we reach equilibrium. Indeed, this approach comes from a method of domain decomposition Dirichlet-Neumann without overlap of Schwarz type. The domain Ω fictitiously decomposed into two domains where Ω int representing the area of influence of the inclusion and Ω ext the rest of the domain. We note u int and u ext respectively the fields in Ω int and Ω ext . We begin by computing the field without inclusion u (0) on the matrix Ω to obtain u 1 ext = u (0) on Γ ext . At the iteration l , we solve a Dirichlet problem on Ω int with the condition u l ext imposed on Γ int to obtain the field u l +1 int + εV l +1 . Finally, we solve on Ω ext a Neumann problem with the condition α m ∂ n u l +1 int + ε∂ n V l +1 on Γ ext . Fig. 3.1 explicates the decomposition. Now, since Ω int ⊂ Ω, we have Γ int = Γ ext = Γ N and u l ext = u l ext on Γ N which reduces the problem solved at each iteration to the next weak problem :

Find u = U on Γ d such that : α m Ω ∇ ∇ ∇u •∇ ∇ ∇w dx = α m ε Γ N ∂ n Vw dσ dσ, (3.35)
with the particular form for the Dirichlet condition U = sin(θ) on all ∂Ω (Γ d ).

Table 3.1 summarizes and details the different steps of an iteration of the method for the Poisson problem. Fig. 3.9 illustrates the α m = 10 stiffness matrix in blue and the interface of Neumann in red with, of course, with the Dirichlet condition U on the outer boundary of the matrix as illustrated in Fig. 3.1.

We then evaluate the L 2 -and H 1 -norms of the difference between the reference solution (a refined solution containing an inclusion) and the first-order corrected solution. For the L 2 -norm this difference for a non corrected solution is 2.27 * 10 -3 then decreases Step Operation carried out 0 interpolate the precalculated V 10 and V 01 on their influence area using scaling around inclusion position p i to obtain V i 01 and V i 10 1 compute u without inclusion For the H 1 -norm this difference for a non corrected solution is 1.84 * 10 -2 then decreases to 8.22 * 10 -3 after the first iteration to stabilize also around the value of 8.16 * 10 -3 after 4 iterations as shown in Table 3.2.

2 compute ∇ ∇ ∇(u)(p i ) = (a i b i ) t at the positions of the inclusions p i 3 actualize u = u + i a i V i 10 + b i V i 01 4 prescribe εα m (a i ∂ n V i 10 + b i ∂ n V i 01 )
In the next iterations the norms oscillate weakly around these two values. This may be due to the value of our corrector on the Neumann interface which is theoretically zero but never will be in practice since it would necessitate the addition of another corrective term.

We then present the case of two circular inclusions placed at positions (0, 0) and (5, 0). Corrections are made by two Neumann interfaces as shown in Fig. 3.10. We evaluate also the L 2 -and H 1 -norms of the difference between the reference solution (a refined solution containing two inclusion) and the first-order corrected solution. For the L 2 -norm this difference for a non corrected solution is 3.92 * 10 -3 then decreases to 2.76 * 10 -3 after the first iteration to stabilize around the value of 2.06 * 10 -3 after 4 iterations. For the H 1 -norm this difference for a non corrected solution is 2.60 * 10 -2 then decreases to 1.18 * 10 -2 after the first iteration to stabilize also around the value of 1.17 * 10 -2 after 4 iterations as

Adaptation for the plain strain elasticity problem

The adaptation to the linear elasticity is done without major difficulty. We begin by transforming as we have done for the Poisson problem the integral posed on the infinite domain Ω ∞ in an integral posed on the bounded domain Ω * as following :

Ω ∞ (ψ 1 (y)) : E m : (ψ 2 (y)) dy = Ω * (B -T ∇ ∇ ∇ * ψ1 (y * ) +∇ ∇ ∇ * ψ1 (y * )B -1 ) : E m : (B -T ∇ ∇ ∇ * ψ2 (y * ) +∇ ∇ ∇ ψ2 (y * )B -1 )| det(B)| dy * .
We then summarize the iterative strategy presented for the Poisson problem by presenting the new basis to be precalculated and which is composed this time of four elements.

In fact, we precalculate :

CHAPITRE 3. A NUMERICAL STRATEGY FOR THE ACCOUNT OF SMALL INCLUSIONS IN ELASTIC MEDIA
to 1, 9 * 10 -2 after the first iteration to stabilize also around the value of 1.9 * 10 -2 after 4 iterations as shown in Table 3.5.

Iterations 0 1 4 L 2 -norm 0.0041788229 0.0031518115 0.0030774998 H 1 -norm 0.029179903239 0.019063487882 0.019055390244 TABLEAU 3.5 -The evolution of the L 2 -and H 1 -error norms as a function of the number of iterations for the case of one circular inclusion for the linear elasticity problem.

We then present the case of two circular inclusions placed at positions (0, 0) and (5, 0). Corrections are made by two edges of Neumann as shown in Fig. 3.10.

We evaluate also the L 2 and H 1 norms of the difference between the reference solution (a refined solution containing two inclusion) and the first-order corrected solution. For the L 2 -norm this difference for a non corrected solution is 5.97 * 10 -3 then decreases to 4.57 * 10 -3 after the first iteration to stabilize around the value of 4.45 * 10 -3 after 4 iterations. For the H 1 -norm this difference for a non corrected solution is 4.08 * 10 -2 then decreases to 2.7 * 10 -2 after the first iteration to stabilize also around the value of 2.7 * 10 -2 after 4 iterations as shown in Table 3 Finally, we present the case of a beam embedded on its right side, subjected to a constant loading on its upper face and reinforced by two square inclusions of side 2ε placed near the connection of embedment as shown in Fig. 3.11. The obtained results are presented in Table 3.7.

We notice that the corrections are more important since the inclusions are placed close enough to the embedded connection which stiffens the behavior of the structure. Iterations 0 1 4 L 2 -norm 0.0088638967 0.0044846349 0.0045123176 H 1 -norm 0.071276543794 0.028118468676 0.02800230894 TABLEAU 3.7 -The evolution of the L 2 -and H 1 -error norms as a function of the number of iterations for the case of two square inclusions for the linear elasticity problem.

Conclusions and perspectives

The purpose of this chapter was to approximate the first-order asymptotic expansion which decomposes in the form of the field without inclusion corrected by the first-order corrector. The field without inclusion is approximated by the classical finite-element method without major difficulty.

We begin in Section 3.4 by presenting the numerical method (Inverted Finite-Element Method) that allows us to approximate the first-order corrector, an error estimate is also established (3.32). Then, according to this estimate, we perform in Section 3.5 the convergence analysis of our first-order corrector in the appropriate weighted Sobolev space. The convergence rates are close enough to the theoretical rates given by the estimate. Finally, we conclude by presenting an iterative strategy taking into account several inclusions. We find that the iterative strategy is efficient and provides a satisfactory correction. Nevertheless, to absorb all the error we will have to go to higher orders.

Natural perspectives of this work are to extend the approximation to higher order correctors then to non linear behavior laws.

Introduction

Nowadays, inclusion reinforcing rubber and rubber-like materials and structures are in increasing use in different industrial products and engineering applications. These materials and structures have multiple scales and are heterogeneous. Therefore, we have to deal with problems related to composite materials and composite structures, since they are composed by different phases, constituents or components at certain spatial scales of observation. At the sample length scale, fillers (inclusions) such as black carbon and silica have to be added to rubber (matrix) materials to enhance the mechanical properties. This reinforcement induces change in modulus and stiffness, stress softening effect, wearing and tearing resistance, etc HAMED [2000]; JEAN [2009b]; JEAN et collab. [2011]; MARK et [START_REF] Mark | The Science and Technology of Rubber[END_REF]; TRELOAR [1975]. At a larger length scale, rubber structures are reinforced by different types of inclusions (long fibers, laminates, etc) to improve stiffness and strength properties among others. As an example, we cite rubber tires reinforced with steel GENT et WALTER [1985]. Generally, the inclusion reinforcement length scale is small enough compared to the reinforced matrix length scale. When the stiffness of reinforced rubber and rubber-like materials and structures is reasonably well understood, their strength is less understood and a catastrophic failure can occur caused by a local stress concentration at the matrix/inclusion interface. Thus, the availability of detailed information on stress and deformation fields around inclusions is of a paramount importance in the design of performing reinforced materials and structures. Nevertheless, in contrast to metallic materials and structures, the deformation and stress fields around inclusions in rubber-like materials and structures are more complicated due to their geometrical and material nonlinear properties, which makes investigating elastostatic fields relatively difficult and multiply the analysis methods : experimental MARTIN-BORRET [1998]; THIBO-DEAU et WOOD [1938]; TRELOAR [1975], theoretical ESHELBY [1957,1959]; [START_REF] Lignon | Modélisation multi-échelles de nappes fibrées en compression[END_REF]; SEITZ et TURNBULL [1956], and numerical JEAN [2009a]; [START_REF] Lignon | Modélisation multi-échelles de nappes fibrées en compression[END_REF] among others. Following these analyses, physical experiments as well as numerical analysis by the finite-element method (FEM) computations show that very high stress concentrations can occur in the vicinity of inclusions and near interfaces where the material parameters are discontinuous causing crack development. In fracture mechanics, the different fracture criteria are based on stress field in the body and therefore a good knowledge of its analytical/numerical form is important. However, experimental and theoretical approaches are generally available for academic simplified problems. Furthermore, a direct numerical approach, based on the finite-element method, can be applied by modeling each constituent or component of the reinforced material or structure. Hence, both properties and geometrical distributions should be taken into account by refining the mesh near inclusions. Consequently, this modeling approach induces prohibitive computational costs. Moreover, it is well known that the convergence rates of standard FEM decrease due to the presence of stress singularities. Here, the probable knowledge of the stress field can be used to develop improved algorithms where e.g. special singular functions are included in the FE-spaces. To conclude, the direct numerical method will encounter difficulties when used to solve these multi-scale problems.

In Chapters 2 and 3 of this thesis, we have developed mathematical and numerical approaches to treat a multi-scale inclusion/matrix problem for linear elastic behavior. The objective of this chapter is to extend this work to nonlinear elastic behavior and to give a first approach to analyze a rubber matrix reinforced by a small inclusion problem.

A literature overview provides several methods for the evaluation and the prediction CHAPITRE 4. MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALS of the elastostatic fields for the inclusion/matrix problem depending in the manner in which the fine heterogeneous scale is treated. Namely, two methods shall be drawn : (1) solving the boundary-value inclusion/matrix problem at a certain length scale of observation when the problem can be considered as continuum, (2) the multi-scale method when two length scales are defined at least.

In the linearized setting, the pioneering elegant work of Eshelby ESHELBY [1957, 1959]; SEITZ et TURNBULL [1956] on ellipsoidal inclusion is considered as masterpieces of the classical theory of elasticity. Following Eshelby's seminal work, an intense research on two-dimensional inclusion problems was done by numerous scientists combining the Eshelby idea with the powerful complex-variable formulation MUSKHELISHVILI [1953] adding complexities such as arbitrary shape geometry, imperfect interface, nonhomogenous boundary conditions at infinity, etc MURA [1987]; PARNELL [2016]; SHAHZAD et collab. [2017] (see these references for an extensive bibliography). Within the framework of finite deformation and nonlinear elasticity which is the scope of this chapter OGDEN [1997], during the past five decades, only few works have been focused on the analysis of the strain and stress fields around a defect, inclusion, crack, notch, etc. This is due to the formidable complexity of the mathematical problem [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF] which makes the boundary-value problem equations highly nonlinear and very difficult to solve analytically or even numerically, as they depend essentially on the form of the hyperelastic potential. For harmonic materials introduced by JOHN [1960], the elastostatic fields can be expressed by analytic functions of a complex variable making possible the use of [START_REF] Muskhelishvili | Some Basic Problems of the Mathematical Theory of Elasticity[END_REF] complex formalism to solve some nonlinear boundary-value problems VARLEY et CUMBERBATCH [1980]. The extension of complex formalism to nonlinear plane problem was elegantly formulated by RU [2002] and new analytical solutions have been obtained for a variety of inclusions with different shapes (inclusions and nanoinclusions) and infinite boundary conditions (uniform and nonuniform) KIM et SCHIA-VONE [2009]; WANG et collab. [2004]. For others hyperelastic potentials, some scientists studied the problem of a finite/semi-infinite two-phase composite reinforced by circular/elliptic cylinder and determined the elastostatic fields in each phase (matrix and fiber) for neoHookean materials and some particular boundary conditions (anti-plane, torsion, etc) DEBOTTON [2005]; DEBOTTON et HARITON [2005]. Another way to model the inclusion/matrix problem is to consider the medium as continuum and to use the anisotropic hyperelastic formalism for solving the boundary-value problem SPENCER [1984].

The different methodologies reviewed above for the inclusion/matrix problems are addressed for particular geometrical shapes, model behavior and boundary conditions at a single scale of observation. The link between different scales is the objective of the multi-scale method. There are different approaches treating multi-scale problems depending on their classifications. In this work, the multi-scale problems are classified on two categories. The first one, following the idea of equivalent inclusion method proposed by Eshelby ESHELBY [1957, 1959], has the objective to elaborate a macroscopic constitutive modeling based on the microscopic models of each phases, constituents or components and the prediction of the effective properties KANOUTÉ et collab. [2009]. Examples of applications of this methodology concern the determination of effective mechanical properties of the hyperelastic effective potential [START_REF] Ogden | Extremum principles in non-linear elasticity and their application to composites[END_REF] and some of them are exact results for particular microstructures (i.e., constitutive law, size, shape, and orientation of the fibers) DEBOTTON [2005]; DEBOTTON et HARITON [2005]. Even if this approach has succeeded in determining the average properties of heterogeneous materials and structures, it could not predict the stress concentration which preclude the failure. The second approach of the multi-scale method treats the class of multi-scale problems dealing with defects near which the macroscopic effective models deduced from the first approach are invalid (inclusions, cracks, dislocations, perturbations, etc). It is based on the asymptotic multi-scale method AMMARI et KANG [2004]; AMMARI et collab. [2002]; DAUGE [1988]. The main idea of this method is to consider the solution of matrix/inclusion problem as a perturbed one from the solution of the unperturbed matrix geometry. Taking into account the smallness of the characteristic length scale of the inclusion compared to the matrix one, multi-scale asymptotic expansion can then be performed. Numerous contributions in the linearized setting are done BENDALI et collab. [2009]; BONNETIER et TRIKI [2013]; DAMBRINE et VIAL [2007]; VIAL [2003].

The aim of this chapter is to analyze a paradigm problem of infinitely long cylindrical perfectly bonded fiber/matrix hyperelastic solid in anti-plane transformation. In fact, it is well known that the anti-plane transformation is the simplest possible setting within the finite deformation elastic theory. Such deformations have been the subject of much attention and some necessary and sufficient mathematical conditions, that have restricted the hyperelastic potential form, are given by KNOWLES [1976,1977] for incompressible and compressible materials to admit non-trivial states of anti-plane shear. This paradigm problem is much simpler to study analytically/numerically and serve as a pilot problem which can introduce the more complex planar or three-dimensional deformations. In Section 4.2, the boundary value problem is formulated in a fully nonlinear Lagrangian framework and in unified formalism for compressible and incompressible hyperelastic behavior. Section 4.3 is concerned with particular hyperelastic potentials : neo-Hookean and Hadamard materials for incompressible and compressible behaviors, respectively. We give the explicit boundary-value problem equations and we show that is nearly similar to the Laplace boundary-value problem treated in Chapter 2. The boundary-value problem associated to a general hyperelastic potential is presented in Section 4.4 which is highly nonlinear to determine the elastostatic fields analitically. Then, its weak formulation coupled with the multi-scale asymptotic expansion of the displacement field is presented.

Formulation of the general anti-plane shear problem

Consider a long cylindrical composite body made up by two homogeneous isotropic incompressible hyperelastic materials : a long cylindrical fiber perfectly bonded to a long cylindrical matrix. Classically, we say that the matrix is reinforced by a fiber, which is considered as an inclusion. In its undeformed free configuration, this body occupies an infinite cylindrical region R ε ⊂ R Ω ε m , connected by a common interface Γ ε , have characteristic lengths respectively 2ε and 2R, respectively, in such a way that ε is small enough compared to R.

The rectangular Cartesian coordinate system (X 1 , X 2 , X 3 ) is chosen in a such manner CHAPITRE 4. MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALS that X 3 is parallel to the R ε generator. These regions are defined as :

R ε = {X | (X 1 , X 2 ) ∈ Ω,-∞ < X 3 < +∞} = R ε f ∪ R ε m , (4.1a) R ε f = X | (X 1 , X 2 ) ∈ Ω ε f , -∞ < X 3 < +∞ , (4.1b) R ε m = X | (X 1 , X 2 ) ∈ Ω ε m , -∞ < X 3 < +∞, Γ = ∂Ω ε m , (4.1c
)

Ω = Ω ε f ∪ Ω ε m ∪ Γ ε , (4.1d)
where X = (X 1 , X 2 , X 3 ) is the position vector of the particle in the undeformed configuration

and Ω denotes a connected cross-section of R ε with a smooth exterior boundary ∂Ω ε m = Γ as already presented in Fig. 1.1. The cylinder R ε is subjected to prescribed boundary conditions on its lateral boundary such that the position vector of the particle X in the undeformed configuration R ε is mapped to x in the deformed configuration R ε * with a transformation parallel to the generators of the cylinder and independent of its axial position. Hence, the deformable solid body occupying the region R ε in its undeformed configuration is said to undergo an anti-plane deformation such that :

x 1 = X 1 , (4.2a)
x 2 = X 2 , (4.2b)

x 3 = X 3 + u ε (X 1 , X 2 ). (4.2c)
The transformation (4.2) is admissible if the out-of-plane displacement u ε (X 1 , X 2 ) is oneto-one and twice continuously differentiable function at least in the open cross-section Ω ε k of the cylindrical region R ε k . To describe the geometrical deformation, the secondorder gradient tensor F, quantity that measures the deformation in the neighborhood of X, is introduced :

F = ∇ ∇ ∇(x), F i , j = ∂x i ∂X j = x i , j (i , j = 1, 2, 3) in Ω ε k (k ∈ { f , m}), ( 4.3)
where ∇ ∇ ∇(.) is the gradient operator with respect to material coordinates. Furthermore, the physical requirement of the global material impenetrability imposes to the transformation x to be one-to-one on R ε k . Thus,

x(X) = x(X ) ∀X, X ∈ R ε k if only if X = X (k ∈ { f , m}). (4.4)
The local form of (4.4) imposes that the deformation, defined by relation (4.3), must satisfy the following condition :

J = det(F) > 0 in Ω ε k (k ∈ { f , m}), ( 4.5)
in which J = det(F) is the Jacobian of the deformation. The geometrical inequality constraint (4.5) ensures that the transformation (4.2) is invertible and admissible for compressible material behavior. If the material behavior is considered incompressible, condition (4.5) is reduced to the so-called constraint of incompressibility :

J = det(F) = 1 in Ω ε k (k ∈ { f , m}), ( 4
F i j =    1 0 0 0 1 0 u ε ,1 u ε ,2 1    (i , j = 1, 2, 3) in Ω ε k (k ∈ { f , m}). (4.7)
Since the transformation gradient (4.7) induced by the antiplane transformation expression (4.2) has J = det(F) = 1 regardless of material behavior, the transformation (4.2) is an isochoric transformation.

The left Cauchy-Green second-order deformation tensor B and its components matrices associated to the transformation (4.2) are :

B = FF T , B i j =    1 0 u ε ,1 0 1 u ε ,2 u ε ,1 u ε ,2 1 + |∇ ∇ ∇u ε | 2    in Ω ε k (k ∈ { f , m}). ( 4
I 1 = tr(B) = 3 + |∇ ∇ ∇u| 2 , ( 4.9a)
I 2 = 1 2 [(tr(B)) 2 -tr(B 2 )] = 3 + |∇ ∇ ∇u| 2 , ( 4.9b)
I 3 = det(B) = 1. (4.9c)
Where tr(• • •) denotes the trace operator. In an undeformed state I 1 = I 2 = 3 and I 1 ≥ 3, I 2 ≥ 3 for others states of deformations. For a hyperelastic material, the existence of an elastic potential function per unit undeformed volume is assumed depending only on the three invariants of B (4.9)

W = W(I 1 , I 2 , I 3 ) in R ε k (k ∈ { f , m}). ( 4

.10)

The mechanical response of the material is then characterized by a strain energy (4.10). The first Piola-Kirchoff stress tensor τ (force per unit undeformed area) is related to the deformation by [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF] :

τ = ∂W(F) ∂F = 2 ∂W(B) ∂B F, (4.11)
for compressible materials and :

τ = ∂W(F) ∂F -qF -T = 2 ∂W(B) ∂B F -qF -T , (4.12)
for incompressible materials while q(X 1 , X 2 , X 3 ) is an unknown Lagrange multiplier field resulting from the incompressibility constraint. The Cauchy stress tensor, which expresses the force per unit deformed area, can be deduced by the use of : (I 1 I -B)B -qI, (4.17)

σ = JτF T . ( 4
for the Cauchy stress. These expressions are analogous to those associated with incompressible constitutive law where q(X 1 , X 2 , X 3 ) is an unknown Lagrange multiplier field resulting from the constraint of incompressibility and depending on (X 1 , X 2 , X 3 ). If the material behavior is compressible, the field q takes the following form :

q = -2 ∂W(I 1 , I 2 , I 3 ) ∂I 3 I 3 =1 , (4.18)
and is linked to the out-of-plane displacement u ε (X 1 , X 2 ) though (4.18) and the invariant relations (4.9) and depending only on (X 1 , X 2 ). Equilibrium equations in the absence of a body force can be expressed as :

Di v (τ) = 0 in R ε k (k ∈ { f , m}), ( 4.19)
where Di v (• • •) is the divergence operator with respect to material coordinates. At this stage, the bond and boundary conditions will be specified. The transformation (4.2) and its out-plane displacement are required to satisfy the bond conditions at the interface Γ ε and the boundary conditions at the external boundary Γ.

The bond conditions at the interface impose the continuity of the transformation and the first Piola-Kirchoff stress across the interface Γ ε :

x f (X 1 , X 2 , X 3 ) = x m (X 1 , X 2 , X 3 ) ∀(X 1 , X 2 ) ∈ Γ ε and ∀X 3 , (4.20a) τ f (X 1 , X 2 , X 3 )N = τ m (X 1 , X 2 , X 3 )N ∀(X 1 , X 2 ) ∈ Γ ε and ∀X 3 , (4.20b)
where N is the unit outward normal vector to the interface Γ ε defined in the initial configuration.

On the external boundary Γ, the boundary conditions imposed on the transformation x should be compatible with the kinematic loading conditions (anti-plane loading) as follows :

u ε = U ∀(X 1 , X 2 ) ∈ Γ and ∀X 3 on lateral boundary of R ε . (4.21)
Generally, the anti-plane shear boundary conditions are applied uniformly for most studies in particular for linear behavior. Nevertheless, the nonuniformity in the applied boundary conditions can help us to predict and understand the highly varying nonlinear elastostatic fields near the matrix-inclusion interface. To do this, boundary conditions (4.21) allow nonuniform loading.

µ k > 0, (4.23a) h(1) = 0, (4.23b) ∂h(J) ∂J J=1 = - µ 2 , (4.23c) ∀J > 0, ∂ 2 h(J) ∂ 2 J ≥ 0. (4.23d)
The constitutive law can be expressed. Then, the first Piola-Kirchhoff stress tensor can be deduced from (4.16) and the hyperelastic potential (4.22) :

τ = ∂W(F) ∂F -qF -T = µ k F -qF -T in R ε k (k ∈ { f , m}). (4.24)
In which the components are :

τ 11 = τ 22 = τ 33 = µ -q, ( 4.25a)
τ 12 = τ 21 = 0, (4.25b)

τ 13 = qu ε ,1 , (4.25c
)

τ 23 = qu ε ,2 , (4.25d)
τ 31 = µu ε ,1 , (4.25e)
τ 32 = µu ε ,2 , (4.25f)
and q = µ, (4.26)

for compressible material behavior.

In the deformed configuration, the Cauchy stress tensor σ is expressed in terms of potential energy by

σ = ∂W(F) ∂F F T -qI = µ k B -qI in R ε * k (k ∈ { f , m}). ( 4

.27)

Its components are :

σ 11 = σ 22 = µ -q, ( 4.28a)
σ 12 = σ 21 = 0, (4.28b)

σ 13 = σ 31 = µu ε ,1 , (4.28c
)

σ 23 = σ 32 = µu ε ,2 , (4.28d
)

σ 33 = µ(1 + |∇ ∇ ∇u ε | 2 ) -q. (4.28e)
Bearing in mind the equilibrium equations (4.19) together with the constitutive law (4.25) and the transformation gradient expression (4.7), one obtain for α = 1, 2 :

(µ k -q) ,α + q ,3 u ,α = 0 in Ω ε k (k ∈ { f , m}), (4.29a) (µ k u ,α) ,α = q ,3 in Ω ε k (k ∈ { f , m}). (4.29b)
The resulting equations (4.29) constitute a system of three differential equations for two unknown functions u ε and q. [START_REF] Knowles | On finite anti-plane shear for incompressible elastic materials[END_REF][START_REF] Knowles | A note on anti-plane shear for compressible materials in finite elastostatics[END_REF] has shown that the Neo-Hookean For compressible behavior, the field q is determined by relation (4.26). For incompressible behavior and since u ε does not depend on X 3 , then (4.29b) requires that q ,3 is a linear function of X 3 : q

(X 1 , X 2 , X 3 ) = d k X 3 + p(X 1 , X 2 ) in R ε k (k ∈ { f , m}). (4.30)
Then, we substitute this expression in (4.29b), which gives :

p(X 1 , X 2 ) = µ k + d k u ε (X 1 , X 2 ) + b k in Ω ε k (k ∈ { f , m}), ( 4.31)
where d k and b k are constants depending on the domain Ω ε k . Hence, the equilibrium equations (4.29) reduce to :

q(X 1 , X 2 ) = µ k + d k (X 3 + u ε (X 1 , X 2 )) + b k , (4.32a) µ k ∆u ε (X 1 , X 2 ) = d k in R ε k (k ∈ { f , m}), ( 4.32b)
where d k = b k = 0 for compressible behavior.

To satisfy bonds conditions on interface Γ, the following conditions are deduced :

u ε f (X 1 , X 2 ) = u ε m (X 1 , X 2 ) on Γ ε , (4.33a) µ f ∂ n u ε f (X 1 , X 2 ) = µ m ∂ n u ε m (X 1 , X 2 ) on Γ ε , (4.33b) µ f -p f = µ m -p m , (4.33c
)

p k = c st (k ∈ { f , m}). ( 4.33d)
So we deduce from (4.33) and the expression of q given in (4.32a) that : Finally, the boundary-value problem can be stated :

d f = d m = d , (4.34a
∆u ε f = d µ f in Ω ε f , (4.35a) ∆u ε m = d µ m in Ω ε m , (4.35b
)

u ε f = u ε m on Γ ε , (4.35c) µ f ∂ n u ε f = µ m ∂ n u ε m on Γ ε , (4.35d) u ε m = U on Γ. (4.35e)
The Cauchy stress tensor σ components are :

σ 11 = σ 22 = -d X 3 + d u ε (X 1 , X 2 ) + b, ( 4.36a)
σ 12 = σ 21 = 0, (4.36b)

σ 13 = σ 31 = µ k u ε ,1 , (4.36c
)

σ 23 = σ 32 = µ k u ε ,2 , (4.36d
)

σ 33 = µ k |∇ ∇ ∇u ε | 2 -d X 3 -d u ε (X 1 , X 2 ) -b. (4.36e)
The boundary-value problem announced below is identical to the one of isotropic linear elasticity except of the existence of an axial normal stress σ 33 due to the nonlinear effect contrary to the linear case and a non homogeneous term that appears in the inclusion domain. It belongs to the family of elliptic transmission problems GRISVARD [1992].

General multi-scale asymptotic expansion

Problem (4.35) is quite similar to Problem (1.1), the only difference is that there is one additional term to be offset for incompressible behavior (for compressible behavior this term vanishes since d = 0).

In fact, we build the asymptotic expansion as done before by introducing the difference

d ε 0,k = u ε k -u (0) for k ∈ { f , m} which is solution to : ∆d ε 0, f (x) = d µ f - d µ m in Ω ε f , (4.37a) ∆d ε 0,m (x) = 0 in Ω ε m , (4.37b) d ε 0, f (x) = d ε 0,m (x) on Γ ε , (4.37c) µ f ∂ n d ε 0, f (x) = µ m ∂ n d ε 0,m (x) + (µ m -µ f )∂ n u (0) (x) on Γ ε , (4.37d) d ε 0,m (x) = 0 on Γ. (4.37e)
To compensate the term that appears on Ω ε f we introduce a similar problem to Problem (2.5) :

∆W (0) f (y) = d µ f - d µ m in Ω 1 f , (4.38a) µ m ∆W (0) m (y) = 0 in Ω ∞ , (4.38b)
W (0) f (y) = W (0) m (y) on Γ 1 , (4.38c)
µ f ∂ n W (0) f (y) = µ m ∂ n W (0) m (y) on Γ 1 , (4.38d)
W (0) m (y) → 0 when |y| → ∞, (4.38e)
and obtain the following K-th order expansion : f is constant and does not need to be expanded as it was the case in all the compensations made beforehand.

u [K] ε,k (x) = u (0) (x) + K i =2 ε i v (i ) (x) + K i =1 ε i V (i -1) k ( x ε ) + ε 2 W (0) ( x ε ). ( 4

Multi-Scale asymptotic analysis for generalized Neo-Hookean materials 4.4.1 Associated boundary value problem

For hyperelastic material, the existence of an elastic potential function W(I 1 , I 2 , I 3 ) per unit undeformed volume is assumed. Since not all incompressible isotropic hyperelastic potentials can sustain nontrivial an anti-plane transformation KNOWLES [1976,1977], we introduce a particular class of elastic materials KNOWLES [1976,1977] and JIANG et KNOWLES [1991] depending only on the first and the third invariants I 1 = tr(B) and I 3 = det(B), respectively : where the functions f and g are twice differentiable and satisfy the following requirements :

W(I 1 , I 3 ) = W(I 1 , I 3 = 1) + ∂W(I 1 , I 3 ) ∂I 1 I 3 =1 f (I 3 ) + g (I 3 ) in R ε k (k ∈ { f , m}), ( 4
f (I 3 = 1) = 0 in R ε k (k ∈ { f , m}), (4.41a) ∂f (I 3 ) ∂I 3 I 3 =1 = -1 in R ε k (k ∈ { f , m}), ( 4.41b
)

g (I 3 = 1) = 0 in R ε k (k ∈ { f , m}), (4.41c) ∂g (I 3 ) ∂I 3 I 3 =1 = 0 in R ε k (k ∈ { f , m}), ( 4.41d)
to recover the infinitesimal linear elastic behavior and to satisfy the stress free initial configuration. In this case of compressible behavior, the field q (4.18) is reduced to :

q = -2 ∂W(I 1 , I 2 , I 3 ) ∂I 3 = 2 ∂W(I 1 , I 2 ) ∂I 1 I 3 =1 . ( 4

.42)

For incompressible elastic behavior, the strain energy (4.40) is reduced as well to the generalized Neo-Hookean hyperelastic potential by setting JIANG et KNOWLES [1991] :

f (I 3 = 1) = 0 in R ε k (k ∈ { f , m}), (4.43a
)

g (I 3 = 1) = 0 in R ε k (k ∈ { f , m}). ( 4
τ = 2 ∂W(I 1 , I 2 , I 3 ) ∂I 1 I 3 =1 F -qF -T in R ε k (k ∈ f , m). ( 4

.44)

In which the components are on each Ω ε k (k = f , m) :

τ i i = 2 ∂W ∂I 1 -q, i ∈ {1, 2, 3}, (4.45a
)

τ α3 = qu ε ,α , α ∈ {1, 2}, (4.45b
)

τ 3α = 2 ∂W(F) ∂I 1 u ε ,α , α ∈ {1, 2}, (4.45c
)

τ 12 = τ 21 = 0, (4.45d)
where the field q(X1, X2, X3) is an unknown Lagrange multiplier field for incompressible material behavior and is linked to the displacement u by the relation (4.42) for compressible behavior. The Cauchy stress tensor can be deduced :

σ = ∂W(F) ∂F F T -qI = 2 ∂W(I 1 , I 2 , I 3 ) ∂I 1 I 3 =1 B -qI in R ε * k (k ∈ { f , m}). ( 4
σ i i = 2 ∂W ∂I 1 -q, i ∈ {1, 2}, (4.47a
)

σ α3 = σ 3α = 2 ∂W ∂I 1 u ε ,α , α ∈ {1, 2}, (4.47b
)

σ 33 = 2 ∂W ∂I 1 (1 + |∇ ∇ ∇u ε | 2 ) -q, ( 4
∂W ∂I 1 -q ,α + q ,3 u ε ,α = 0 in R ε k (k ∈ { f , m}), (4.48a) 2 ∂W ∂I 1 u ε ,α ,α = q ,3 in R ε k (k ∈ { f , m}). (4.48b)
For compressible behavior the field q is done by (4.42). For incompressible behavior and since only the field q(X 1 , X 2 , X 2 ) depends on X 3 in equations (4.48) it takes the form :

q(X 1 , X 2 , X 3 ) = d k X 3 + p(X 1 , X 2 ) in R ε k (k ∈ { f , m}). ( 4.49)
Using this expression in the first of (4.48) :

p(X 1 , X 2 ) = 2 ∂W ∂I 1 + d k u ε (X 1 , X 2 ) + b k in Ω ε k (k ∈ { f , m}), ( 4
(X 1 , X 2 , X 3 ) = 2 ∂W ∂I 1 + d k (X 3 + u ε (X 1 , X 2 )) + b k in R ε k (k ∈ { f , m}), (4.51a) 2 ∂W ∂I 1 u ε ,α ,α = d k in R ε k (k ∈ { f , m}), ( 4.51b)
where d k = b k = 0 for compressible behavior.

To satisfy bond conditions at the interface Γ ε (4.20), the following conditions are deduced :

u ε f (X 1 , X 2 ) = u ε m (X 1 , X 2 ), (4.52a) 2 ∂W f ∂I 1 -q f = 2 ∂W m ∂I 1 -q m , (4.52b) 2 ∂W f ∂I 1 ∂ N u ε f (X 1 , X 2 ) = 2 ∂W m ∂I 1 ∂ N u ε m (X 1 , X 2 ) ∀X 3 on Γ ε . (4.52c)
So we deduce from (4.52) and the expression of q given in (4.51a) that :

d f = d m = d , (4.53a) b f = b m = b. ( 4
∂W ∂I 1 u ε ,α ,α = d k in Ω ε k , (4.54a) u ε f = u ε m on Γ ε , (4.54b) ∂W f ∂I 1 ∂ n u ε f = ∂W m ∂I 1 ∂ n u ε m on Γ ε , (4.54c) u ε m = U on Γ. (4.54d)
The components of the Cauchy stress tensor σ are :

σ 11 = σ 22 = -d X 3 -d u ε (X 1 , X 2 ) -b, ( 4.55a)
σ 12 = σ 21 = 0, (4.55b)

σ 13 = σ 31 = ∂W ∂I 1 u ε ,1 , (4.55c
)

σ 23 = σ 32 = ∂W ∂I 1 u ε ,2 , (4.55d
)

σ 33 = 2 ∂W ∂I 1 |∇ ∇ ∇u ε | 2 -d X 3 -d u ε (X 1 , X 2 ) -b. (4.55e)
The boundary value-problem announced above is identical to the one of isotropic linear elasticity except of the existence of eventually a non homogeneous term in the right of equation (4.54a) for incompressible behavior. It belongs to the family of elliptic transmission problems GRISVARD [1992]. We remark the existence of an axial normal stress σ 33 due to the nonlinear effect contrary to the linear case.

First-order multi-scale asymptotic expansion 4.4.2.1 From strong to weak formulations

We begin by presenting the equations satisfied by the solution without inclusion in the strong form : Di v τ(u (0) ) = 0 in Ω (0) , (4.56a) u (0) = U on Γ, (4.56b) and its associated weak formulation :

Find u (0) : Ω (0) -→ R 3 such that : Ω (0) τ(u (0) ) : ∇ ∇ ∇h dx = 0 in Ω (0) , (4.57a) u (0) = U on Γ,

u ε m = U on Γ, (4.59d
)

h m = 0 on Γ, (4.59e)
for all test functions h m and h f satisfying the conditions above.

First Piola-Kirchhoff tensor expansion

For this we expand, first of all, the first Piola-Kirchoff tensor :

The scaled weak formulation

We present here a scaled weak formulation using the same scaling introduced previously y = x ε and introducing the scaled difference D 0,k (y) = d 0,k (x) (k ∈ { f , m})). Recall that the goal is to work on a fixed geometry.

We obtain the following weak scaled formulation : Let us now also expand T k (u (0) )(εy), k ∈ { f , m} :

T f (u (0) )(εy) = T f (u (0) )(0) + O(ε), (4.66a)
T m (u (0) )(εy) = T m (u (0) )(0) + O(ε). The idea now is to obtain an approximation of the first-order corrector using the finiteelement method and the inverted finite elements one as already done for the isotropic linear elasticity case. Only the determination of the field without inclusion u (0) remains non linear. ) solved by the first order corrector has a variational formulation that may not be coercive. This is due to the loss of coercivity of T f and T m and also to the loss of isotropy which implies that the decay of the first-order corrector to 0 at ∞ is no longuer assured. This loss of coercivity implies that the first-order expansion may not exist, at least in this form where the exponents are integers. One way to remedy this is to say thet there may exist an expansion but not in integer exponents (see for example CALOZ et collab. [2006] for which an analogous phenomenon exists when there is a singularity of geometry).

Conclusion and perspectives

This Chapter has the particularity to be exploratory in the sense that we have investigated the possibility of extending the work done in linear elasticity to hyperelasticity. The boundary-value problem associated with isotropic hyperelastic inclusion/matrix materials under the anti-plane deformation kinematic is formulated. The first application, with a Neo-Hookean hyperelastic potential, allowed us to deduce the explicit analytical multiscale asymptotic expansion for the anti-plane displacement. It was shown that there exists an axial Cauchy stress due to the problem geometrical nonlinearity. The boundary-value problem of the more complicated general isotropic hyperelastic potential is also formulated. The complexity of the associated partial differential equations does not enable us to determine the analytical elastostatic fields. To overcome this problem, a weak formulation of the strong problem with the use of the anti-plane displacement multi-scale asymptotic expansion is deduced. A numerical strategy similar to the one formulated in Chapter 3 can be adopted and developed. ANNEXE A.

A.2 An estimation for small domain

For f ε (x) of the form f ε (x) = g ε ( x ε ) with g ε (y) bounded uniformly in ε ∈ (0, 1) and y ∈ Γ 1 we state the following lemma (established under the assumption that Ω 1

f = Ω 1 f ε is a Lipschtiz domain) : Lemma A.2.1. Let f ε ∈ L 2 (Γ ε ). Then Γ ε f ε v dσ ≤ Cε -1 2 f ε L 2 (Γ ε ) v H 1 (Ω ε f ) (A.16)
where C is a constant independant of ε, f ε and of v ∈ H 1 (Ω ε f ).

Proof. We write that

Γ ε f ε v dσ ≤ f ε L 2 (Γ ε ) v L 2 (Γ ε )
and then that

v L 2 (Γ ε ) = ε -1 2 ε 1 2 v L 2 (Γ ε ) + [v] 2,Γ ε with [v] 2,Γ ε = Γ ε XΓ ε | f (x) -f (x )| 2 |x -x | 2 dσ x dσ x 1 2
.

It is established in the page 212 of that there exists a constant C independant of ε and of v ∈ H 1 (Ω ε f ) such that

ε 1 2 v L 2 (Γ ε ) + [v] 2,Γ ε ≤ C inf w∈H 1 0 (Ω ε f ) v + w H 1 (Ω ε f ) .
Hence the estimation (A.16) by taking w = 0.

A useful estimation can then been deduced from (A.16). Indeed, by simply noticing that by changing the variable y = x ε , we have

ε -1 2 f ε L 2 (Γ ε ) = g ε (.) L 2 (Γ 1 )
and using the fact that for v ∈ H 1 (Ω), we have

v H 1 (Ω ε f ) ≤ v H 1 (Ω)
we reach the following estimation Γ ε f ε v dσ ≤ C g ε (.) L 2 (Γ 1 ) v H 1 (Ω) .

(A.17) Ce travail de thèse a concerné la modélisation mathématique et l'approximation numérique de l'influence d'une inclusion très fine sur un substrat élastique de différente rigidité. L'étude est motivée par les applications dans les pneumatiques et ne se base pas sur des techniques d'homogénéisation classiques. En effet, l'objectif a été de traiter l'intéraction entre une seule inclusion et son milieu élastique et non une densité d'inclusions. L'étude a comporté trois volets, le premier concernant une modélisation mathématique pour des lois de comportement linéaires aboutissant à une expression de la contribution de l'inclusion sous la forme du champ sans inclusion corrigé par des correcteurs à différents ordres. Ces correcteurs sont indépendants de la taille caractéristique de l'inclusion. Le second a concerné l'approximation numérique de cette influence moyennant la méthode des éléments finis et celle des éléments finis inversés. Une stratégie numérique de prise en compte de l'influence de plusieurs inclusions y est aussi présentée. Le dernier volet est prospectif et discute de la possibilité de l'extension de l'approche pour des lois de comportement non linéaires.

A.3 Références

Mots clés :

Analyse asymptotique, analyse multi-échelle, méthode des éléments finis, méthode des éléments finis inversés, élasticité linéaire, élasticité non linéaire

Asymptotic and numerical analysis of thin inclusions in elastic media

Abstract :

This work focused on mathematical modeling and numerical approximation of the influence of a very thin inclusion on an elastic substrate of different stiffness. The study is motivated by applications in tires and is not based on conventional homogenization techniques. Indeed, the objective was to treat the interaction between a single inclusion and its elastic medium and not a density of inclusions. The study consisted of three parts, the first concerning mathematical modeling for linear behavior laws leading to an expression of the contribution of the inclusion in the form of the inclusion-free field corrected by correctors at different orders. These correctors are independent of the characteristic size of the inclusion. The second relates to the numerical approximation of this influence by means of the finite element method and that of the inverted finite elements. A numerical strategy for taking into account the influence of several inclusions is also presented. The last part is prospective and discusses the possibility of extending the approach for nonlinear behavioral laws. Ce travail de thèse a concerné la modélisation mathématique et l'approximation numérique de l'influence d'une inclusion très fine sur un substrat élastique de différente rigidité. L'étude est motivée par les applications dans les pneumatiques et ne se base pas sur des techniques d'homogénéisation classiques. En effet, l'objectif a été de traiter l'interaction entre une seule inclusion et son milieu élastique et non une densité d'inclusions. L'étude a comporté trois volets, le premier concernant une modélisation mathématique pour des lois de comportement linéaires aboutissant à une expression de la contribution de l'inclusion sous la forme du champ sans inclusion corrigé par des correcteurs à différents ordres. Ces correcteurs sont indépendants de la taille caractéristique de l'inclusion, Le second a concerné l'approximation numérique de cette influence moyennant la méthode des éléments finis et celle des éléments finis inversés. Une stratégie numérique de prise en compte de l'influence de plusieurs inclusions y est aussi présentée. Le dernier volet est prospectif et discute de la possibilité de l'extension de l'approche pour des lois de comportement non linéaires.

MOTS-CLÉS

εΩ 1 f« 3 3 1. 4

 1334 . Noter que Γ = ∂Ω n'est pas nécessairement de forme circulaire. . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Géometrie dilatée. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Une illustration du simplexe infini T ∞ et de son simplexe fini associé S dans le cas bidimensionnel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Circular inclusion and domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Taking into account several inclusions using artificial Neumann conditions. 3.2 An illustration of the altitude vector h, the finite simplex and the supporting hyper-plane associated to a 2D infinite simplex. . . . . . . . . . . . . . . . . . 3.3 Inverted finite element method decomposition and correspondence . . . . 3.4 Conforming mesh of the contour of the inclusion. . . . . . . . . . . . . . . . . 3.5 The fictitious domain Ω * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Behavior of V m at ∞ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Convergence curve l og (wei g ht ed -nor m) = f (l og (h)) for P 1 finite elements. 3.8 Convergence curve l og (wei g ht ed -nor m) = f (l og (h)) for P 2 finite elements. 3.9 The correction edge of Neumann . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 The case of two inclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 The case of two square inclusions placed close to the mounting connection ix Liste des tableaux 3.1 The operations performed during the iterations for the Poisson problem. . 3.2 The evolution of the L 2 -and H 1 -error norms as a function of the number of iterations for the case of one inclusion for the Poisson problem. . . . . . . . 3.3 The evolution of the L 2 -and H 1 -error norms as a function of the number of iterations for the case of two inclusions for the Poisson problem. . . . . . . . 3.4 The operations performed during the iterations for the linear elasticity problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 The evolution of the L 2 -and H 1 -error norms as a function of the number of iterations for the case of one circular inclusion for the linear elasticity problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 The evolution of the L 2 and H 1 error norms as a function of the number of iterations for the case of two circular inclusions for the linear elasticity problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 The evolution of the L 2 -and H 1 -error norms as a function of the number of iterations for the case of two square inclusions for the linear elasticity problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Rien de grand ne s'est accompli dans le monde sans passion » Contexte scientifique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Etat de l'art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Résumé des travaux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.6.1 Chapitre 2 : Développement asymptotique multi-échelle de l'influence d'une inclusion fine sur un substrat élastique. . . . . . . . . 8 1.6.2 Chapitre 3 : Une stratégie numérique de prise en compte de l'influence d'une inclusion fine sur un substrat élastique . . . . . . . . 11 1.6.3 Chapitre 4 : Développement asymptotique multi-échelle de l'influence d'une inclusion fine sur un substrat élastique pour l'hyperélasticité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.7 Conclusion et perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.8 Références . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

CHAPITRE 1 .

 1 PRÉSENTATION DE LA THÈSE Ω ε f ⊂ R 2 un domaine borné et connexe de dimension charactéristique ε représentant l'inclusion (aussi appellée fibre) de taille très petite par rapport à la dimension charactéristique de la matrice représentée par Ω ε m ⊂ R 2 de rayon R et de bord éxtérieur régulier Γ. Le développement asymptotique étant effectué par rapport au paramètre ε, nous désignons par Ω 1 f le domaine tel que Ω ε f = εΩ 1 f et nous considérons Ω ε f comme un ensemble paramétré de domaines. L'interface Γ ε = ∂Ω ε

  (1.1e) où α f et α m sont les modules de cisaillement des deux matériaux, u ε f est la solution du problème dans le domaine de l'inclusion Ω ε f , u ε m est la solution du problème dans le domaine de la matrice Ω ε m , U est un déplacement préscrit sur Γ et ∂ n représente la dérivée normale (extérieure pour Ω ε f , intérieure pour Ω ε m ). Pour U ∈ H 1 2 (Γ), l'existence et l'unicité d'une solution faible du problème précédent est garantie par le théorème de Lax-Milgram.

  la solution du problème dans le domaine de la matrice Ω ε m et U un déplacement préscrit. Les vecteur n f et n m dénotent les normales unitaires sortantes à Ω ε f et Ω ε m respectivement. L'existence et l'unicité du problème précédent est donnée par le théorème de Lax-Milgram pour une donnée U ∈ H 1 2 (Γ, R 2 ). Le tenseur des déformations linéaire est défini par = 1 2 (∇ ∇ ∇u+∇ ∇ ∇ t u) et le tenseur des contraintes CHAPITRE 1. PRÉSENTATION DE LA THÈSE peut s'écrire σ(u) = λt r ( )I + 2µ où tr et I sont respectivement la trace et le tenseur identité et (•) t est la transposée de (•). Fig. 1.1 rend explicite la géométrie.

CHAPITRE 1 .

 1 PRÉSENTATION DE LA THÈSE 1.6.1 Chapitre 2 : Développement asymptotique multi-échelle de l'influence d'une inclusion fine sur un substrat élastique.

CHAPITRE 1 .

 1 PRÉSENTATION DE LA THÈSE analogue à celui résolu par les V (p) k dont les solutions Z (p) k sont posées sur des domaines infinis. Le prolème résolu par le résidu à l'ordre zéro d ε 0

2. 4

 4 Conclusion and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Références . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  (p) m at infinity Recall that the functions V (p) m solve Problem (2.5). Let B(O, R) be the ball of center O and radius R > 0 such that Ω 1 f ⊂ B(O, R). Then, the trace of V (p) m on the curve ∂B(O,R)∩Ω ∞

  of the antiplane problem, and a new term -div (E f :(u (0) ))(x) in Ω ε f . CHAPITRE 2. MULTI-SCALE ASYMPTOTIC EXPANSION FOR A SMALL INCLUSION IN ELASTIC MEDIAStill using the scaling y = x ε and denoting Ω

  s + 2ρ) cos((s -1)θ) -(1 + s + 2ρ) sin((s -1)θ) , (2.88) whith ρ = µ (λ+µ) , generate the expansion at infinity of the functions V

Chapitre 3 A 46 3. 4

 3464 numerical strategy for the account of small inclusions in elastic media « On résout certains problèmes par approximation, en négligeant de petites quantités. » Jean le Rond d'Alembert Sommaire 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.2 Main motivation and process to follow . . . . . . . . . . . . . . . . . . . . 44 3.3 A numerical method for the problem in infinite domain . . . . . . . . . . 45 3.3.1 The adequate logarithmic weighted Sobolev space and some properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Inverted finite-element method : Details and convergence . . . . . . . . 47 3.4.1 Implementation of the method . . . . . . . . . . . . . . . . . . . . . . 53 3.5 Computational results and discussions . . . . . . . . . . . . . . . . . . . . 55 3.5.1 Convergence analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.5.2 Dirichlet-Neumann iterative strategy . . . . . . . . . . . . . . . . . . 57 3.5.3 Adaptation for the plain strain elasticity problem . . . . . . . . . . . 59 3.6 Conclusions and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.7 Références . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

FIGURE 3 . 1 -

 31 FIGURE 3.1 -Taking into account several inclusions using artificial Neumann conditions.

Fig. 3 .

 3 Fig. 3.3 illustrates the decomposition.

  0

  .21) and 〈ϕ, w〉 = Γ 1 ϕw dσ.

(3. 22 )

 22 We now fix a parameter θ > 0. Given a real-valued function v defined on Ω ∞ , we define the function v over Ω * as follows v(y * ) = r (y * ) -θ+1 v(Φ -1 (y * )) for y * ∈ Ω * . (3.23) We may observe that v(y) = r (y) -θ+1 v(Φ(y)) for y ∈ Ω ∞ , (3.24) and v(y * ) = v(y * ) for y * ∈ Ω * ∩ Ω ∞ . (3.25)

CHAPITRE 3 .

 3 A NUMERICAL STRATEGY FOR THE ACCOUNT OF SMALL INCLUSIONS IN ELASTIC MEDIA

  .29) with a(•, •) and ϕ(•) are given by (3.21) and (3.22).

CHAPITRE 3 .

 3 A NUMERICAL STRATEGY FOR THE ACCOUNT OF SMALL INCLUSIONS IN ELASTIC MEDIAand can be computed as in the finite-element method.

  with y = Φ(y * ). Thus, Ω∞ α m ∇ψ 1 (y) • ∇ψ 2 (y) dy = Ω * α m (∇ * ψ1 (y * )) T B -1 B -T ∇ * ψ2 (y * )| det(B)| dy * where B is the Jacobian matrix of the transformation Φ defined by B = ∂Φ i ∂y * j 1≤i , j ≤2 . The edges of the two domains Ω * and Ω b are coupled via a Lagrange multiplier. The implantation and the calculation of B and det(B) had been developed and taken into account using the open source finite-element library GetFEM++ RENARD et POMMIER [2014].

FIGURE 3 . 4 -FIGURE 3 . 5 -

 3435 FIGURE 3.4 -Conforming mesh of the contour of the inclusion.

FIGURE 3 . 6 -FIGURE 3

 363 FIGURE 3.6 -Behavior of V m at ∞

FIGURE 3

 3 FIGURE 3.8 -Convergence curve l og (wei g ht ednor m) = f (l og (h)) for P 2 finite elements.

CHAPITRE 3 .

 3 A NUMERICAL STRATEGY FOR THE ACCOUNT OF SMALL INCLUSIONS IN ELASTIC MEDIA

FIGURE 3 . 9 -

 39 FIGURE 3.9 -The correction edge of Neumann

TABLEAU 3 . 3 -

 33 FIGURE 3.10 -The case of two inclusions

FIGURE 3 .

 3 FIGURE 3.11 -The case of two square inclusions placed close to the mounting connection

CHAPITRE 4 .

 4 MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALS

  3 composed by an open bounded cylindrical regionR ε f ⊂ R 3 with a connected cross section Ω ε f ⊂ R 2 surrounded by an open bounded cylindrical region R ε m ⊂ R 3 with a connected cross section Ω ε m ⊂ R 3 . The cross sections Ω ε f ,

  .8) Here | • • • | denotes the Cartesian norm of a vector and the notation (• • •) T is used to represent the transpose operator. The scalar fundamental invariants of the left Cauchy-Green secondorder tensor deformation B combined with the transformation (4.2) are :

  .13) Combining the relation (4.13) and the τ expression (4.12) gives .Since the transformation (4.2) is isochoric, the constitutive law relating stress to deformation can be written in the unifying way :

CHAPITRE 4 .

 4 MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALS et KNOWLES [1991] :

CHAPITRE 4 .

 4 MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALSisotropic hyperelastic material admits non-trivial states of anti-plane shear and then the equations (4.29) have non-trivial solutions.

  ) b f = b m = b.(4.34b)

CHAPITRE 4 .

 4 MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALS

.39) Remark 4 . 3 . 1 .

 431 The compensation is made only once since the term appearing in Ω ε

  .40) CHAPITRE 4. MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALS

  .50) where d k and b k are constants depending on the domain Ω ε k (k ∈ { f , m}). Hence the equilibrium equations (4.48) reduce to, q

( 4 .

 4 57b) for all test functions h such that h = 0 on Γ.We present then the strong formulation associated to the differenced ε 0,k = u ε k -u (0) for k ∈ { f , m} : Di v (τ(u (0) + d ε 0,k ) -τ m (u (0) )) = 0 in Ω ε k (k ∈ { f , m}), (4.58a)τ f N f = τ m N m on Γ ε , (4.58b) d ε 0, f = d ε 0,m on Γ ε ,(4.58c)d ε 0, f = 0 on Γ,

( 4 .

 4 58d) CHAPITRE 4. MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALS and its associated weak formulation :Ω ε f (τ(u (0) + d ε 0, f ) -τ m (u (0) )) : ∇ ∇ ∇h f dx + Ω ε m (τ(u (0) + d ε 0,m ) -τ m (u (0) )) : ∇ ∇ ∇h m dx = 0,(4.59a)h f = h m on Γ ε ,

  τ(u (0) + d ε 0,k ) = τ(u (0) ) + T k (u (0) )∇ ∇ ∇d ε 0,k + O( d ε with T k (u (0) ) = ∂ 2 W ∂F 2 (u (0) ).Using d ε 0,k the first-order approximation of d ε 0,k and (4.60), the weak formulation reads :Find d ε 0,k : Ω -→ R 3 such that : Ω ε f T f (u (0) )∇ ∇ ∇d ε 0, f : ∇ ∇ ∇h f dx + Ω ε m T m (u (0) )∇ ∇ ∇d ε 0,m : ∇ ∇ ∇h m dx = Ω ε f (τ m (u (0) ) -τ f (u (0) )) : ∇ ∇ ∇h m dx,(4.61a)h f = h m on Γ ε , (4.61b) d ε 0, f = d ε 0,m on Γ ε ,(4.61c)d ε 0,m = 0 on Γ,(4.61d)h m = 0 on Γ.

TT

  f (u (0) )(εy)∇ ∇ ∇D ε 0, f (y) : ∇ ∇ ∇ y H f (y) dy+ Ω ε m ε T m (u (0) )(εy)∇ ∇ ∇ y D ε 0,m (y) : ∇ ∇ ∇ y H m (y) dy] = Ω 1 f (τ m (u (0) )(εy) -τ f (u (0) )(εy)) : ∇ ∇ ∇ y H f (y) dy. (4.62) CHAPITRE 4.MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALSTo lighten the writings we set : (u (0) )(εy) -τ f (u (0) )(εy)) : ∇ ∇ ∇ y H f (y) dy.The weak scaled formulation has the form : f (u (0) )(εy)∇ ∇ ∇D ε 0, f (y) : ∇ ∇ ∇ y H f (y) dy+ Ω ε m ε T m (u (0) )(εy)∇ ∇ ∇ y D ε 0,m (y) : ∇ ∇ ∇ y H m (y) dy] = A. (4.63)The next step is to expand τ k (u (0) )(εy), k ∈ { f , m} :τ f (u (0) )(εy) = τ f (u (0) )(0) + O(ε),(4.64a)τ m (u (0) )(εy) = τ m (u (0) )(0) + O(ε). (4.64b) We thus obtain the expression of A : A = Ω 1 f (τ m (u (0) )(0) -τ f (u (0) )(0)) : ∇ ∇ ∇ y H f (y) dy + O(ε). (4.65)

TT

  f (u (0) )(0)∇ ∇ ∇D ε 0, f (y) : ∇ ∇ ∇ y H f (y) dy + Ω ε m ε T m (u (0) )(0)∇ ∇ ∇ y D ε 0,m (y) : ∇ ∇ ∇ y H m (y) dy] = Ω 1 f (τ m (u (0) )(0) -τ f (u (0) )(0)) : ∇ ∇ ∇ y H f (y) dy. (4.67)We finally arrive to this new formulation :f (u (0) )(0)∇ ∇ ∇D ε 0, f (y) : ∇ ∇ ∇ y H f (y) dy + Ω ε m ε T m (u (0) )(0)∇ ∇ ∇ y D ε 0,m (y) : ∇ ∇ ∇ y H m (y) dy] = Ω 1 f (τ m (u (0) )(0) -τ f (u (0) )(0)) : ∇ ∇ ∇ y H f (y) dy,(4.68a)H f = H m on Γ 1 , (4.68b) D ε 0, f = D ε 0,m on Γ 1 ,(4.68c)D ε 0,m = 0 on Γ ε ,(4.68d)H m = 0 on Γ ε . (4.68e) CHAPITRE 4. MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALSfor all test function H satisfying the conditions above. Now, as ε is considered to be small compared to the size of the domain and to eliminate completely the dependence on ε, we approximate Problem (4.68) by the following problem posed on the dilated domain :Find D ε 0,k : R 2 -→ R 3 such that : 1 ε [ Ω 1 f T f (u (0) )(0)∇ ∇ ∇D ε 0, f (y) : ∇ ∇ ∇ y H f (y) dy + Ω ∞ T m (u (0) )(0)∇ ∇ ∇ y D ε 0,m (y) : ∇ ∇ ∇ y H m (y) dy] = Ω 1 f (τ m (u (0) )(0) -τ f (u (0) )(0)) : ∇ ∇ ∇ y H f (y) dy,(4.69a)H f = H m on Γ 1 ,(4.69b)D ε 0, f = D ε 0,m on Γ 1 ,(4.69c)D ε 0,m → 0 when |y| → ∞,(4.69d)H m → 0 when |y| → ∞. (4.69e)for all test function H satisfying the conditions above.The residue D ε 0,k solves a problem of linear elasticity but this time anisotropic. We then propose, by analogy with what has been done for the isotropic linear elasticity, to write the residue D ε 0,k in the form of a first corrector weighted by epsilon and a rest :D ε 0,k = εV (0) k + O(ε 2 ). (4.70)By replacing the first-order expansion of D ε 0,k (4.70) in Problem (4.69), we obtain the weak formulation satisfied by the first ordet corrector V(0) k :Find V (0) k : R 2 -→ R 3 such that :Ω 1 f T f (u (0) )(0)∇ ∇ ∇ y V (0) f (y) : ∇ ∇ ∇ y H f (y) dy + Ω ∞ T m (u (0) )(0)∇ ∇ ∇ y V (0) m (y) : ∇ ∇ ∇ y H m (y) dy = Ω 1 f (τ f (u (0) )(0) -τ m (u (0) )(0)) : ∇ ∇ ∇ y H f (y) dy,(4.71a)H f = H m sur Γ 1 , m → 0 quand |y| → ∞,(4.71d) H m → 0 quand |y| → ∞. (4.71e) for all test function H satisfying the conditions above.

CHAPITRE 4 .

 4 MULTI-SCALE ASYMPTOTIC ANALYSIS FOR ISOTROPIC HYPERELASTIC MATRIX-INCLUSION MATERIALS Remark 4.4.1. Problem (4.71

  FOLIO ADMINISTRATIF THESE DE L'UNIVERSITE DE LYON OPEREE AU SEIN DE L'INSA LYON NOM : BEN HASSINE DATE de SOUTENANCE : 26/09/2017 Prénoms : Mohamed Rafik TITRE : Étude asymptotique et numérique d'inclusions fines dans des domaines élastiques NATURE : Doctorat Numéro d'ordre : 2017LYSEI086 École doctorale : EDA512 InfoMaths Spécialité : Mathématiques et applications RÉSUME :

1.6.3 Chapitre 4 : Développement asymptotique multi-échelle de l'in- fluence d'une inclusion fine sur un substrat élastique pour l'hy- perélasticité

	Ce chapitre fait l'objet d'une généralisation du travail précedemment effectué pour
	des lois de comportement non linéaires. En effet, le comportement d'un pneumatique,
	certes élastique, est loin d'être linéaire et ses composants composés de matériaux caou-
	tchoutiques sont régis par des lois dites hyperélastiques dérivant d'un potentiel comme
	suit :	σ =	∂W(F) ∂F	F T -qI.	(1.24)
	La recherche du potentiel approchant au mieux le comportement réel de ces matériaux
	est d'actualité et le restera toujours faisant l'objet d'un travail fastidieux alliant essais mé-
	caniques et méthodes d'identification pour donner naissance à des potentiels représen-
	tant de la manière la plus fiable le comportement réel des matériaux en question. Les
	expressions des potentiels obtenues sont de plus en plus compliquées comportant un
	nombre de paramètres qui ne cesse d'augmenter et sont dans la plupart des cas gardées
	secrètes par les industriels. L'enjeu devient alors de construire un développement asymp-
	totique multi-échelle représentant l'influence d'une inclusion fine sur un domaine hy-
	perélastique indépendant de l'expression du potentiel W en question. Ceci donnera à la
	modélisation effectuée un degrés de liberté supplémentaire permettant surtout au simu-
	lateur d'être indépendant. La perspective d'intégrer la méthode dans un code de calcul est
	envisagée en créant un objet inclusion par exemple dont les entrées seraient taille, forme,
	position et potentiel. Il ne restera plus alors à l'opérateur que d'introduire l'expression du
	potentiel.
	Ce chapitre est subdivisé en deux parties, la première concernant une classe particulière
	de matériaux gouvernée par un potentiel Neo-Hookéen simple :

3 Asymptotic expansion for the plane strain linear elasticity problem . . .

  Lyon. 7 SÁNCHEZ-PALENCIA, E. 1974, «Problèmes de perturbations liés aux phénomènes de conduction à travers des couches minces de grande résistivité», J. Math. Pures Appl., vol. 53, n o 9, p. 251-269. 3 .3.4 Error estimate for the linear elasticity problem . . . . . . . . . . . .

	Chapitre 2
	Multi-scale asymptotic expansion for a
	small inclusion in elastic media
	(p) m and Z	(p)

AMMARI, H. et H. KANG. 2004, Reconstruction of small inhomogeneities from boundary measurements, Springer-Verlag, Berlin Heidelberg. 4 ANTIPOV, Y. A. et P. SCHIAVONE. 2003, «On the uniformity of stresses inside an inhomogeneity of arbitrary shape», Int. J. Eng. Sci., vol. 68, p. 299-311. 4 BENDALI, A., P. H. COCQUET et S. TORDEUX. 2012, «Scattering of a scalar time-harmonic wave by n small spheres by the method of matched asymptotic expansions», Numerical Analysis and Applications, vol. 5, n o 2, p. 116-123. 4 BENDALI, A., P.-H. COCQUET et S. TORDEUX. 2016, «Approximation by multipoles of the multiple acoustic scattering by small obstacles in three dimensions and application to the foldy theory of isotropic scattering», Archive for Rational Mechanics and Analysis, vol. 219, n o 3, p. 1017-1059. 4 BENDALI, A., M. FARES, E. PIOT et S. TORDEUX. 2013, «Mathematical justification of the rayleigh cavity model with the method of matched asymptotic expansions», SIAM J. Appl. Math., vol. 71, n o 1, p. 438-459. 4 BENDALI, A. et K. LEMRABET. 1996, «The effect of a thin coating on the scattering of a time-harmonic wave for the Helmholtz equation», SIAM J. Appl. Math., vol. 56, n o 6, p. 1664-1693. 3 BENDALI, A., A. MAKHLOUF et S. TORDEUX. 2011, «Field behavior near the edge of a microstrip antenna by the method of matched asymptotic expansions», Q. Appl. Math., vol. 69, p. 691-721. 4 CHAPITRE 1. PRÉSENTATION DE LA THÈSE MAZILU, P. 1972, «On the theory of linear elasticity in statically homogeneous media», Rev. Roum. Math. Pur. Appl., vol. 17, p. 261-273. 4 MAZ'YA, V., S. NAZAROV et B. PLAMENEVSKIJ. 2000a, Asymptotic Theory of Elliptic Boundary Value Problems in Singularly Perturbed Domains, Birkhauser Verlag, Basel. 3 MAZ'YA, V., S. NAZAROV et B. PLAMENEVSKIJ. 2000b, Asymptotic Theory of Elliptic Boundary Value Problems in Singularly Perturbed Domains., Birkhauser Verlag, Basel. 3 RENARD, Y. et J. POMMIER. 2014, Getfem++ Short User Documentation, Release 4.3, SCHIAVONE, P. 2003, «Neutrality of the elliptic inhomogeneity in the case of non-uniform loading», IMA journal of applied mathematics, vol. 8, p. 161-169. 4 SUN, Y. et Y. PENG. 2003, «Analytic solutions for the problems of an inclusion of arbitrary shape embedded in a half space», Applied Mathematics and Computation, vol. 140, p. 105-113. 4 VIAL, G. 2003, Analyse multi-échelle et conditions aux limites approchées pour un probleme avec couche mince dans un domaine à coin, thèse de doctorat, Université de Rennes 1. 8 « Dans la nature, tout a toujours une raison. Si tu comprends cette raison, tu n'as plus besoin de l'expérience. » Leonard de Vinci Sommaire 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Asymptotic expansion for the antiplane problem . . . . . . . . . . . . . . 2.2.1 Building the asymptotic expansion . . . . . . . . . . . . . . . . . . . 2.2.2 Existence and uniqueness in unbounded domains for the Poisson problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Expansion of the functions V (p) m at infinity . . . . . . . . . . . . . . . 2.2.4 Error estimate for the antiplane problem . . . . . . . . . . . . . . . . 2.2.5 Solution for a circular inclusion using the separation of variable method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3.1 Building the asymptotic expansion . . . . . . . . . . . . . . . . . . . 2.3.2 Existence and uniqueness in unbounded domains for the linear elasticity problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Expansion of the functions V m at infinity . . . . . . . . . . 2

  In other words, condition(3.16) means that the triangulation T of Ω b must be shape regular in the classical sense. Conditions (3.17)-(3.19) mean that the triangulations T * are graded : the triangles touching the fictitious vertex a 0 have a size of order h

	CHAPITRE 3. A NUMERICAL STRATEGY FOR THE ACCOUNT OF SMALL
	INCLUSIONS IN ELASTIC MEDIA
	min K∈T * * h K ≥ c * 3 h	1 µ .	(3.19)
		1
		µ . Conditions
	(3.18)-(3.19) express a refinement around a 0 in order to capture more or less information
	at infinity.
	We are now in position to approximate Problem (2.5).
	The continuous problem is : find
			.18)

  The evolution of the L 2 -and H 1 -error norms as a function of the number of iterations for the case of one inclusion for the Poisson problem.

	Iteration	0	1	4
	L 2 -norm	0.0022772836	0.0017580537	0.0012025580
	H 1 -norm	0.018486736759	0.008221498939	0.008162436664
	TABLEAU 3.2 -

  .6. The evolution of the L 2 and H 1 error norms as a function of the number of iterations for the case of two circular inclusions for the linear elasticity problem.

	Iterations	0	1	4
	L 2 -norm	0.0059758968	0.0045773017	0.0044582372
	H 1 -norm	0.040875404871	0.027011468644	0.02701120983
	TABLEAU 3.6 -

asymptotique et numérique d'inclusions fines dans des domaines élastiques Résumé :

AMROUCHE, C., V. GIRAULT et J. GIROIRE. 1994, «Weighted Sobolev spaces for Laplace's equation in R n », J. Math. Pures Appl., vol. 73, p. 579-606. 87, 88 Analyse

  :Analyse asymptotique, analyse multi-échelle, méthode des éléments finis, méthode des éléments finis inversés, élasticité linéaire, élasticité non linéaire

	Laboratoire (s) de recherche : Institut Camille Jordan
	Directeur de thèse: Yves Renard
	Président de jury :

Composition du jury : Abderrahmane BENDALI (rapporteur) ,Faouzi TRIKI (rapporteur), Monique DAUGE (examinatrice), Eric

Bonnetier (examinateur)

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI086/these.pdf © [M.R. Ben Hassine], [2017], INSA Lyon, tous droits réservés

Les matériaux composés de couches de fibres sont très présents dans l'industrie. Ils sont notamment utilisés dans les sports automobiles et dans l'aéronautique, secteurs de pointe dans lesquels le compromis rigidité/poids est très important. Ces composites sont constitués de deux parties : la matrice et des renforts fibrés. Les fibres assurent la solidité et la rigidité du matériau. La matrice, quant à elle, sert uniquement à maintenir les fibres solidaires les unes aux autres. Elle est donc souvent choisie peu résistante et assez légère. Un matériau composite est généralement stratifié, c'est à dire que plusieurs couches de fibres (appelées "plis") sont empilées. Ces plis peuvent être constitués soit de fibres placées toutes dans le même sens (plis unidirectionnels), soit d'un tissage de fibres (comme pour un tissu classique). Ce type d'agencement (fibres matrice) est présent dans le pneumatique et joue un rôle primordial dans la conception et la fabrication de ce dernier. L'inclusion de renforts de toutes dimensions et formes au sein de la structure du pneu (ou autre) est un élément de conception clé permettant de tirer le meilleur parti d'agencements composites. Bien sûr, de tels renforts ou inclusions, pour être intéressants, auront forcément des propriétés de raideur très différentes de celles du substrat. En découlent des singularités locales de contrainte, susceptibles de générer des endommagements locaux. Aussi, il importe d'être capable de simuler de tels agencements de manière efficace, en décorrelant les discrétisations du substrat et des inclusions afin de gagner en flexibilité sur l'évaluation des solutions de conception. Cela impose la définition de schémas de couplage empêchant tout verrouillage, avec un niveau de précision au rendez-vous. Afin

Remerciements

V 1000 for ϕ 1000 = (E m -E f ) :

1 0 0 0

• n,

where λ and µ are the Lamé coefficients. The obtained fields are then interpolated on a circular geometry of radius 10 and then stored in four files. The zone of influence is the same. It is the square inscribed in the circle of radius 1.

At each iteration, (u (0) )(0) noted a b c d is computed at x = 0 (the position where the inclusion is centered) and then V (0) is updated (V (0) = aV 1000 + bV 0100 + cV 0010 + d V 0001 ) before imposing its epsilon-weighted gradient on the Neumann interface. We reiterate until we reach equilibrium.

The weak problem solved at each iteration is the following :

Find u (0) = U on Γ d such that :

with the particular form for the Dirichlet condition U = (si n(θ) si n(θ)) t on all ∂Ω (Γ d ).

Table 3.4 summarizes and details the different steps of an iteration of the method for the linear elasticity problem.

Step

Operation carried out 0 interpolate the precalculated V 1000 , V 0100 , V 0010 and V 0001 on their influence area using scaling around inclusion position p i to obtain V i 1000 , V i 0100 , V i 0010 and

go to step 1 until the norms reach equilibrium TABLEAU 3.4 -The operations performed during the iterations for the linear elasticity problem.

We now test, as done for the Poisson problem, the case of one circular inclusion placed at x = 0 for the linear elasticity problem with λ f = 1000, µ f = 2000, λ m = 10 and µ m = 20.

We evaluate the L 2 and H 1 norms of the difference between the reference solution (a refined solution containing an inclusion) and the first-order corrected solution. For the L 2 -norm this difference for a non corrected solution is 4.17 * 10 -3 then decreases to 3.15 * 10 -3 after the first iteration to stabilize around the value of 3.07 * 10 -3 after 4 iterations. For the H 1 -norm this difference for a non corrected solution is 2.91 * 10 -2 then decreases The boundary value problem may now be defined as follows. Given an hyperelastic potential W(I 1 , I 2 , I 3 ) (4.10) with the constitutive law (4.16)-(4.17) we seek the unknown anti-plane displacement u ε on Ω ε k and eventually the unknown Lagrange multiplier field q on each R ε k for an incompressible material behavior such that the equilibrium equations (4.19) hold and satisfy the bond conditions (4.20) and the boundary conditions (4.21). This boundary-value problem is an over determined system of three partial differential equations with one or two unknown fields. In order to solve this over determined boundary-value problem and to have a nontrivial solution, KNOWLES [1976,1977] gave some restrictions on the form of the strain energy W(I 1 , I 2 , I 3 ) which are necessary and sufficient conditions.

The solution of the boundary value problem announced below depend essentially on the form of the strain energy W(I 1 , I 2 , I 3 ) (4.10). Since no restriction has been imposed on the nature of this functional, it may be questionable if the predicted response done by the constitutive behavior (4.12) or (4.14) is physically reasonable and if the solution of the finite elasticity problem exist and is unique. This is the so-called Truesdell's Haupt problem and it is concerned with stability, thermodynamics and some mathematical conditions. The most known requirements on the strain energy are those of Baker-Ericksen inequality, the Colleman-Noll inequality, ellipticity condition, and polyconvexity condition.

In the following, in order to alleviate mathematical expressions the index (k), unless it is necessary, will be omitted.

Multi-Scale asymptotic analysis for particular hyperelastic matrix-inclusion materials 4.3.1 Formulation of the associated boundary value problem

We introduce a particular class of materials governed by the hyperelastic potential OGDEN [1997] :

where µ f and µ m are the shear modulus of the fiber/inclusion and matrix hyperelastic materials for infinitesimal deformations, respectively. If the functional h(J) = 0, the strain energy defined by (4.22) governs a compressible material behavior and corresponds to the special cases of the Hadamard materials called reduced-order Hadamard hyperelastic potential (compressible Neo-Hookean potential) and studied extensively by [John 1966].

If the material behavior is incompressible, i.e. J = 1, the functional h(J) = 0, the strain energy defined by (4.22) corresponds to the Neo-Hookean potential. Strain energy (4.22) is the simplest phenomenological and molecular constitutive model function of rubber like materials. However, its capability to predict experimental data is poor especially at high values of deformation. Nevertheless, the strain energy model (4.22) is the most used model in finite elasticity to deduce analytic solutions of the associated boundary-value problems [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF]. To satisfy the strong ellipticity condition and to recover the h(J) infinitesimal behavior law, the functional should satisfy the following requirements JIANG

Annexe A A.1 Poincaré-type inequality

We establish the existence and uniqueness of the profile functions V (k) introduced previously in the asymptotic expansion at order K given by (2.26). The functions V (k) are defined on R 2 (two-dimensional geometry) and are solutions to the family of Problems (2.5). For that, we will introduce weighted Sobolev spaces W m,p α,β (W 1,2 0,0 in our case) similar to classical ones W m,p (H 1 in our case) but with weights that describe the growth or the decay of functions at infinity. The idea of using weighted spaces arises naturally from Hardy's inequalities and will allow us to establish a Poincaré inequality relating to norms of functions and to that of their derivatives.

The use of these weights is necessary to obtain a Poincaré-type inequality and eliminates the drawbacks of spaces defined by the closure of D(R 2 ) for the Dirichlet norm and which are not always spaces of distributions. The used weights arises naturally from Hardy's inequality or from a generalized Hardy's inequality (see AMROUCHE et collab. [1994]) and the classical ones are of the form ρ = (1 + |x| 2 ) 1 2 but there is appearance of a logarithmic factor in our case.

A.1.1 An intermediate result

The following result is an intermediate result to prove the equivalence of the norm and the semi-norm.

Lemma A.1.1. For any large enough real number R, there exists a constant C R such that :

Démonstration. Let ϕ belongs to D R . First, observe that, owing to the support of ϕ, all integrals in the norm and semi-norm are taken on B R instead of R 2 . Hence, since the origin is in the interior of B R , we can use r and ln r instead of ρ(r ) and lg(r ) in the expression of the norm and seminorm. Then using ∂ϕ ∂r

= ∇ϕ • r r we can write :

Let θ be the angular variable, then we have :

Now, assuming that R is large enough, we apply the generalized Hardy's inequality (see AMROUCHE et collab. [1994]) (with γ = -2) to the function r → ϕ(r, θ). Integrating with respect to θ and applying (A.1) we obtain :

The needed result can now be proven.

A.1.2 A Poincaré-type inequality

) is a norm on W 1,2 0,0 (R 2 )/P 0 , and that :

Thus, we only have to prove that there exits c > 0 such that :

The proof proceeds in two steps. The first step consists in eliminating the quotient norm by choosing an adequate representative of the class of u. To this end, we fix a bounded open domain of R 2 , with positive measure, say O, and we choose the representative U of u in W 1,2 0,0 (R 2 ) that satisfies the system of equations :

It is easy to see that (A.4) determines U uniquely and that : u W 1,2 0,0 (R 2 )/P 0 ≤ U W 1,2 0,0 (R 2 ) .

(A.5) Therefore, the second step consists in proving that there exists a constant C such that the following bound holds for all U in W 1,2 0,0 (R 2 ) satisfying (A.4) :

We shall prove it by contradiction. If (A.6) is not true, there exists a sequence (U ν ) of elements of W 1,2 0,0 (R 2 ) satisfying (A.5) and such that :

Hence the sequence (U ν ) is bounded in W 1,2 0,0 (R 2 ) and since this is a reflexive Banach space, we can extract a subsequence, still denoted by (U ν ), that converges weakly to an element U * of W 1,2 0,0 (R 2 ) and it is easy to check from this weak convergence that U * also satisfies (A.5). But since |u| W 1,2 0,0 (R 2 ) tends to 0, the lower semi-continuity of the norm implies that |u| W 1,2 0,0 (R 2 ) = 0. Thus, U * is a constant (polynomial of P 0 ) and the fact that U * satisfies (A.5) implies that U * = 0. ANNEXE A. Now, we need a strong convergence to conclude by contradiction, but we cannot use a standard compactness argument on an unbounded domain. Instead, we shall derive a strong convergence via an adequate partition of unity that will enable us to consider separately a bounded domain where the topologies of W 1,2 0,0 and H 1 coincide and the exterior of a ball, where Lemma A.1.1 can be applied.

Let R denote a real number, large enough to apply the generalized Hardy's inequality. Let ϕ and ψ be two functions of C ∞ (R 2 ) such that :

Since for fixed R, W 1,2 0,0 (B R+1 ) is isomorphic to R, H 1 (B R+1 ), we have that U µ converges weakly to 0 in H 1 (B R+1 ). Since H 1 (B R+1 ) is compactly embedded into L 2 (B R+1 ), it follows that U ν → 0 strongly in L 2 (B R+1 ), (A.9)

In addition, as |U ν | W 1,2 0,0 (R 2 ) tends to 0, it follows that U ν → 0 strongly in W 1,2 (B R+1 ), (A.10) so that ϕU ν → 0 strongly in W 1,2 0,0 (B R+1 ). (A.11)

Now, let us examine the behavior of ψU ν . For fixed ν, let (θ j ) be a sequence of functions of D(R 2 ) that tends to U ν in W 1,2 0,0 (R 2 ). Then, ψθ j belongs to D(B R ) and we can apply to it Lemma (A.1.1) ψθ j W 1,2 0,0 (R 2 ) ≤ C R |ψθ j W 1,2 0,0 (B R ) .

(A.12)

Then, letting j tend to infinity and using the fact that ψ is identically one outside B R+1 , we obtain : Next, let ν tend to infinity and observe that ψU ν tends to zero strongly in W 1,2 0,0 (B R ∩ B R+1 ) because B R ∩ B R+1 is bounded and W 1,2 0,0 (B R ∩ B R+1 ) is isomorphic to W 1,2 (B R ∩ B R+1 ) ; we derive that : ψU ν → 0 strongly in W 1,2 0,0 (B R ). (A.14)

Since U ν = ϕU µ + ψU µ , we obtain :

U µ → 0 strongly in W 1,2 0,0 (R 2 ), (A.15) which contradicts the assumption (A.7) that :

U µ W 1,2 0,0 (R 2 ) = 1.