HDR / ENS CACHAN - BRETAGNE présentée par

sous le sceau de l’Université européenne de Bretagne Grégory Vial

pour obtenir Préparée à l’Unité Mixte de Recherche 6625
L’HABILITATION À DIRIGER DES RECHERCHES Institut de recherche mathématique de Rennes
Mention : Mathématiques

Perturbations singulières Habilitation soutenue le 21 juin 2010

Composition du jury :

de problèmes elliptiques, Grégoire ALLAIRE,

Professeur à l’École Polytechnique (Palaiseau) / rapporteur

analyse asymptotique et Éric BONNETIER,

Professeur à l’Université Joseph Fourier (Grenoble) / rapporteur
Michael VOGELIUS,
calcul scientifique Professeur à Rutgers (Piscataway, USA) / rapporteur

François ALOUGES,
Professeur à l’École Polytechnique (Palaiseau) / examinateur
Monique DAUGE,
Directrice de Recherches à l’IRMAR (Rennes) / examinatrice
Arnaud DEBUSSCHE,
Professeur à l’ENS Cachan - Bretagne (Bruz) / examinateur
Michel PIERRE,
Professeur à l’ENS Cachan - Bretagne (Bruz) / examinateur
2
Table des matières

Introduction 5

1 Analyse asymptotique de problèmes elliptiques 7

1.1 Contexte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Superposition multi-échelle et développements raccordés . . . . . . . . . . 8
1.3 Problèmes de couches minces . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.1 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.2 Cas modèle : l’équation de Laplace . . . . . . . . . . . . . . . . . . . 12
1.3.3 Flexion d’une plaque mince avec raidisseur . . . . . . . . . . . . . . 14
1.4 Interaction entre perturbations . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.5 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2 Méthodes numériques pour la prise en compte de micro-défauts en mécanique 23

2.1 Contexte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2 Représentation asymptotique . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3 Couplage avec un modèle à discontinuité forte . . . . . . . . . . . . . . . . 28
2.4 Conditions aux limites approchées pour les profils . . . . . . . . . . . . . . 29
2.5 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3 Programmation pour le calcul scientifique 33

3.1 Contexte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.2 La bibliothèque éléments finis M ÉLINA . . . . . . . . . . . . . . . . . . . . 33
3.3 Calcul des profils pour la mécanique numérique . . . . . . . . . . . . . . . 34
3.4 Utilisation des éléments finis de haut degré . . . . . . . . . . . . . . . . . . 36
3.4.1 Problèmes de singularités . . . . . . . . . . . . . . . . . . . . . . . . 36
3.4.2 Problèmes d’oscillations . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.4.3 Problèmes d’interface . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.5 Exploitation de résultats numériques . . . . . . . . . . . . . . . . . . . . . . 40
3.5.1 Réseau d’électrodes et matriçage cellulaire . . . . . . . . . . . . . . 40
3.5.2 Partitions minimales et partitions nodales . . . . . . . . . . . . . . . 41
3.6 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

Références bibliographiques 49

3
4
Introduction

Ce rapport fait la synthèse de mon activité de recherche en analyse appliquée et calcul

scientifique. On retrouve une certaine homogénéité dans ces travaux, en particulier en
ce qui concerne les outils : analyse asymptotique des problèmes elliptiques, méthodes
d’éléments finis d’ordre élevé. Les applications, en revanche, sont le fait de rencontres
scientifiques et présentent naturellement une plus grande diversité.
La première partie de ce mémoire concerne l’analyse asymptotique de certains pro-
blèmes de perturbation singulière pour les EDP elliptiques. Ce travail s’inscrit dans la
continuité de mon travail de thèse sur les conditions aux limites approchées pour un
problème de couche mince en domaine non régulier. Les équations considérées sont clas-
siques (équation de Laplace, de Kirchhoff-Love ou de Navier), mais la difficulté réside
dans la géométrie sur laquelle elles sont posées : couches minces, inclusions, perturba-
tions du bord. La technique développée est l’analyse multi-échelle, rebaptisée superposi-
tion multi-échelle dans ce rapport pour ne pas laisser penser qu’elle est la seule à pouvoir
décrire un phénomène présentant plusieurs échelles. Une des contributions a d’ailleurs
consisté dans [DTV09] à comparer cette approche avec la méthode des développements
asymptotiques raccordés, très populaire dans les applications. De nombreux résultats
nouveaux ont pu être établis dans le cas des problèmes de remplacement d’une couche
mince par une condition aux limites approchée, voir [Via05, CCDV06, RV08, RV10].
D’autre part, des problèmes d’inclusions ou de perforation au bord ont été envisagés
avec ces techniques asymptotiques dans [DV05, BNDTV09]. Dans le cadre du projet ANR
Jeunes Chercheurs MACADAM, les résultats ont été appliqués aux équations de l’élas-
ticité linéaire. Il s’agit de prédire l’initiation d’une fissure dans un matériau présentant
des inhomogénéités de surface, voir [DV07, BDVV08, BNBD+ 10]. Un des points-clés de
ce travail concerne le calcul de profils définis dans des domaines non bornés. La recherche
de conditions aux limites approchées pour résoudre numériquement de tels problèmes
nous a amenés à considérer un problème de Ventcel sans condition de signe, nécessitant
une analyse pseudo-différentielle, voir [BNDHV10b].
De nombreux calculs numériques ont été menés dans le cas modèle des équations de
Laplace, et de l’élasticité linéaire (où la représentation asymptotique a été couplée avec
un modèle à discontinuité forte pour la description de l’endommagement). Ces simula-
tions, adossées à la bibliothèque de calculs éléments finis M ÉLINA [Mar10], ont nécessité
(et permis...) des développements nouveaux : éléments finis de haut degré [CDMV03],
intégrands, etc. Grâce à ces nouvelles fonctionnalités, le code a été enrichi de nouvelles

5
applications : problèmes de partitions minimales [BNHV10], analyse de la supraconduc-
tivité dans des sections polygonales [BNDMV07], résolution numérique de l’équation de
Cahn-Hilliard [GMV10] ou encore optimisation de la forme d’un réseau d’électrodes pour
le matriçage cellulaire [PV10].

6
Chapitre 1
Analyse asymptotique de problèmes
elliptiques

1.1 Contexte

Dans de nombreuses situations physiques, on est amené à considérer des problèmes

où intervient un petit paramètre ε, souvent lié à la géométrie du matériau considéré. C’est
le cas pour un objet conducteur recouvert d’une fine couche diélectrique, un assemblage
par une pellicule de colle, un canal dont le fond présente des petites rugosités, un mé-
tal avec des aspérités à sa surface ou un matériau contenant des inhomogénéités, etc.
Ces questions interviennent dans tous les champs d’application : propagation d’ondes,
mécanique des structures, acoustique, thermique, mécanique des fluides. . .
La complexité des problèmes a conduit à de riches développements mathématiques
sur ces sujets. De multiples techniques ont été introduites pour décrire le comportement
asymptotique de la solution uε d’une équation aux dérivées partielles dépendant d’un
petit paramètre ε. La finalité de l’étude consiste à obtenir un développement asympto-
tique de uε lorsque ε tend vers 0 (ou du moins les premiers termes). Parmi les techniques
propres aux problèmes aux limites elliptiques, on peut en distinguer deux : la méthode
des développements asymptotiques raccordées, et la superposition multi-échelle, détaillées au
paragraphe 1.2. Si la première, popularisée par Van Dyke [Van75], est très largement ré-
pandue dans les applications, la seconde est souvent utilisée dans les travaux mathéma-
tiques car plus facile à justifier.
On se concentre ici sur certains problèmes de perturbations singulières provenant
de la géométrie du domaine dans lequel est défini l’équation aux dérivées partielles. Le
livre [MNP00] constitue une référence générale sur ce type de problèmes pour des pertur-
bations locales (inclusions ou trous, congés de raccordement. . .) ou non locales, telles les
couches minces. Les perturbations locales peuvent être traitées dans le cadre de la théo-
rie du potentiel, voir [VV00, BHB06, ABC+ 08, AKLZ10] ou des développements multi-
échelle, en particulier en présence de singularités [NO93, CD96, BMNP01]. Quant aux
couches minces, elles ont fait l’objet d’une étude intensive en raison de leur importance
dans les applications en acoustique et électromagnétisme – voir [EN93, BL96, ALG98,
AHH00, HJ01, AKS06, CCDV06, JT06, BL08, Poi08, CH09] – ou en mécanique – voir [SP74,

7
LD91, LT92, GKL99, MS03, KMO04, AMM06, RV08, Rah08]. On peut citer aussi les pro-
blèmes de frontière rugueuse, qui relèvent de la même problématique cf. [APV98, JM01,
ACG06, MV07, BFNW08, Mil09, BM10]. Mentionnons enfin les problèmes d’homogénéi-
sation qui sont apparentés aux précédents même si des techniques de nature différente
sont généralement employées, voir par exemple [CD99, DM93, All02].
Le chapitre qui suit expose quelques contributions de l’auteur et ses collaborateurs à
l’étude asymptotique de perturbations géométriques pour des problèmes elliptiques. La
première partie concernant les techniques de développement asymptotique dans un cas
modèle est volontairement détaillée pour permettre de s’y appuyer dans la suite en vue
de traiter des problèmes plus complexes.
On fait ici la synthèse des articles [Via05, CCDV06, RV08, DTV09, RV10] avec G. Ca-
loz, M. Costabel, M. Dauge, L. Rahmani et S. Tordeux, qui font suite à mes travaux de
thèse, et [DV05, DV07, BNDTV07] effectués en collaboration avec V. Bonnaillie-Noël, M.
Dambrine et S. Tordeux dans le cadre du projet M ACADAM.

1.2 Superposition multi-échelle et développements raccordés

Considérons un problème aux limites, dont la solution est notée uε , mettant en jeu
deux échelles d’espace. L’échelle naturelle – en la variable x – décrit le comportement
régulier de la solution tandis que l’échelle rapide – en variable x/ε – correspond à un
comportement singulier de couche limite près de l’origine.
La technique des développements asymptotiques raccordés consiste à représenter uε de
deux manières différentes, selon qu’on est proche ou loin de l’origine :

k∑


 δk (ε)uk (x) loin de l’origine,
≥0
uε (x) ! (1.1)



 ∑ δk ( ε ) U k x
( ε ) près de l’origine,
k ≥0

où les fonctions uk et U k ne dépendent pas de ε, et la famille (δk (ε))k définit l’échelle

du développement asymptotique (par exemple δk (ε) = εk pour un développement en
puissances entières). Aucun des deux développements (1.1) n’est valide dans tout le do-
maine ; ils doivent coïncider dans une zone intermédiaire au moyen de conditions de rac-
cord, voir [Van75, Il’92].
La technique de superposition multi-échelle vise à obtenir un développement global
dans tout le domaine sous la forme

uε (x) ! ζ ( xε ) ∑k≥0 δk (ε)vk (x) + χ(x) ∑k≥0 δk (ε)V k ( xε ), (1.2)

où ζ et χ sont des fonctions de troncature régulières satisfaisant

% %
0 loin de l’origine, 1 loin de l’origine,
χ(x) = ζ (y) =
1 près de l’origine, 0 près de l’origine.

L’inconvénient majeur de la technique de superposition réside dans la dépendance

des termes vk et V k en les fonctions de troncature, alors que les termes uk et U k sont

8
intrinsèques au problème. En revanche, il est plus aisé de donner un sens précis à la
formule (1.2) en raison de son caractère global.
Dans [DTV09], avec M. Dauge et S. Tordeux, nous avons construit les développements
asymptotiques selon les deux techniques, et les avons comparés. Le problème modèle
étudié est volontairement simple : uε ∈ H10 (Ωε ) est la solution du problème de Laplace
−∆uε = f dans un ouvert bidimensionnel Ωε , polygone dont un des coins – supposé
situé en l’origine – a été “rabotté” ou “arrondi” à l’échelle ε, cf. F IG . 1.1. La solution u0 du
problème limite présente des singularités en l’origine, dépendant de l’angle α du domaine
Ω0 en ce sommet ; suivant la théorie de la régularité des problèmes à coins [Kon67, Gri85,
Dau88, NP94], elles sont données en coordonnées polaires par skλ (x) = r kλ sin(kλθ ), où
λ = π/α :
K
u0 ( x ) = ∑ bk0 skλ (x) + Ox→0 (|x|K ). (1.3)
k =1

Les scalaires bk0 sont appelés coefficients de singularité. L’exposant de singularité λ fournit
l’échelle du développement asymptotique, comme on va le voir ci-dessous.

Ω0 α Ωε

α H∞ α
r∗ r∗
• • •
0 0 R∗ 0 εR∗

F IGURE 1.1 – Exemple de domaine perturbé Ωε = Ω0 ∩ εH∞ .

Superposition multi-échelle
La méthode de superposition est constructive, et peut être assimilée à une prédiction-
correction : partant de la solution limite v0 = u0 , on forme le reste r0ε = uε − v0 qui satisfait
−∆r0ε = 0 dans Ωε mais présente une trace non-nulle au voisinage de l’origine. Précisé-
ment, en utilisant le développement (1.3) et l’homogénéité des fonctions singulières, on
peut écrire pour ε petit
r0ε (εX) = −u0 (εX) ! −ελ b10 sλ (X).
Pour corriger ce premier terme résiduel, on introduit un profil, défini dans le domaine
limite en variable rapide X = x/ε :
Ωε
H∞ = lim ,
ε →0 ε
solution du problème


 −∆V λ = 0 dans H∞ ,

V λ (X) = −b10 sλ (X) pour X ∈ ∂H∞ , (1.4)


 λ
V (X) → 0 lorsque |X| → ∞.
Ainsi, on peut former un début de développement asymptotique sous la forme

u0 (x) + χ(x)εV λ ( xε ),

9
où la fonction χ permet de localiser le profil au voisinage de l’origine afin de conserver
la condition de Dirichlet homogène sur le reste du bord du domaine. Posant r1ε = uε −
u0 (x) − χ(x)εV λ ( xε ), on vérifie aisément que

r1ε = Ox→0 (|x|2λ ),

si bien que la condition de Dirichlet est satisfaite à un ordre supérieur par rapport à r0ε .
Bien sûr, la troncature induit une erreur dans l’équation de Laplace, qu’il faut corriger à
son tour. On peut développer le profil V λ lorsque X tend vers l’infini :
K
V λ (X) = ∑ Bk1 s−kλ (X) + O|X|→∞ |X|−K , (1.5)
k =1

avec Bk1 ∈ R, voir [Néd01, Kon67] et en déduire

−∆r1ε ! −ε2λ B11 [∆, χ]s−λ .

Ici encore, un cran est gagné dans les puissances de ε. En itérant le processus, on construit
le développement complet sous la forme
K K
uε (x) = ∑ εkλ vkλ (x) + χ(x) ∑ εkλ V kλ ( xε ) + rKε ,
k =0 k =1

où le reste rKε satisfait (rKε (H1 (Ωε ) = O (εKλ ), cette dernière étant obtenue à l’aide d’esti-
mations a priori 1 . Notons que seule une fonction de troncature est ici nécessaire grâce à
l’inclusion Ωε ⊂ Ω0 . Dans le cas général, la partie du développement en variable lente
est multipliée par ζ ( xε ) pour en gommer le comportement local en l’origine.

Développements raccordés
L’échelle (εkλ )k ayant été mise en évidence par la méthode précédente, on postule un
ansatz de la forme
& +∞
&
& uε (x) !
&
&
∑ εkλ ukλ (x) loin de l’origine,
k =1
& (1.6)
& +∞
&
& uε (x) ! ∑ ε U ( ε ) près de l’origine.
kλ kλ x
&
k =0

En injectant ces expressions dans le problème satisfait par uε , et en utilisant des dévelop-
pements du type (1.3) et (1.5), on obtient
% %
∆ukλ = f δk0 dans Ω0 , ∆U kλ = 0 dans H∞ ,
ukλ = 0 sur ∂Ω0 . U kλ = 0 sur ∂H∞ .

On a bien sûr u0 = u0 , solution du problème limite. Quant aux termes suivants, ils ne
sauraient être variationnels sans être nuls... Il faut donc autoriser une explosion en l’ori-
gine pour ukλ et une croissance à l’infini pour V kλ . En développant selon les fonctions
1. Insistons sur le fait que les développements asymptotiques présentés sont généralement des séries
non-convergentes. La notion de convergence au sens des développements asymptotiques signifie que le
reste d’ordre K converge vers 0 lorsque ε tend vers 0, ce d’autant plus rapidement que K est grand.

10
singulières (croissantes et décroissantes), on peut rechercher les termes sous la forme
& ' (
& kλ
& u = ∑ akp s− pλ + ckp s pλ ,
& p ≥1
&
& ' (
& U kλ =
&
&
∑ p A k pλ
s + C k − pλ
p s .
p ≥1

Si l’on souhaite que les deux développements (1.6) coïncident dans la zone intermédiaire
|x| * 1 et | xε | + 1, on obtient les conditions de raccord (en utilisant l’homogénéité des
fonctions singulières) &
& k (k− p)
& a p = Cp (0 si k < p),
&
& k
& A p = c(pk− p) (0 si k < p).

Le développement (1.6) peut-être justifié par une technique proche de la superposition

multi-échelle. En effet, on peut construire une approximation globale dans Ωε via la for-
mule
K ) * K
uε = ζ ( ηxε ) ∑ εkλ ukλ (x) + 1 − ζ ( ηxε ) ∑ εkλ U kλ ( xε ) + ρεK (x),
k =0 k =0

où ηε définit la région intermédiaire : ηε → 0 et ηε /ε → +∞. L’estimation du reste obtenue

√
dans le cas le plus favorable (ηε = ε) s’écrit

(ρεK (H1 (Ωε ) = O (εKλ/2 ).

Comparaison des deux techniques

On a déjà souligné la dépendance des termes vkλ , V kλ du développement par super-
position multi-échelle vis-à-vis des fonctions de troncature, comme le montre l’appari-
tion de commutateurs [∆, χ] dans la construction. On a aussi vu que seule une technique
de superposition permet de définir une approximation globale, et par là d’obtenir des
estimations d’erreur. L’analyse rigoureuse dans le cas des développements raccordés ap-
paraît pour la première fois dans [Tor04, JT06].
Chaque technique présente des spécificités, mais il est possible de relier les deux dé-
veloppements obtenus. Il est facile de vérifier que les termes coïncident localement (en
raison de l’estimation optimale de l’erreur sur le reste rKε ) :

ukλ = vkλ loin de l’origine et U kλ = V kλ près de l’origine.

Plus précisément, ils diffèrent de la partie non-variationnelle des termes des développe-
ments raccordés :

 k −1

 vkλ (x) = ukλ (x) − (1 − χ(x)) ∑ akp s− pλ (x),


p =1
 k
 kλ kλ
 V (y) = U (y) − χ(y)

 ∑ Akp s pλ (y).
p =1

Il existe toute une famille de développements situés entre la situation variationnelle des
vkλ , V kλ et le cas totalement non-variationnel donné par ukλ , U kλ .

11
1.3 Problèmes de couches minces
1.3.1 Motivations
Contrairement au problème modèle du “coin arrondi”, les problèmes de couche mince
ne sont pas des perturbations locales de la géométrie. Elles ne sont pas pour autant néces-
sairement plus difficiles à traiter, notamment dans le cas d’un domaine régulier. Il s’agit
d’effectuer un changement d’échelle pour se ramener à une couche d’épaisseur fixe, le
paramètre ε apparaissant alors dans les équations (dont la complexité peut s’accroître
notablement en raison des termes dus à la courbure).
La motivation principale pour obtenir un développement asymptotique dans le cas
des couches minces est la recherche d’un modèle approché ne faisant plus intervenir la
couche. Cette dernière est remplacée par une condition aux limites approchée (parfois
appelée condition d’impédance) dont l’effet est le même à un ordre donné. Le problème
approché est plus simple à résoudre numériquement car ne nécessite pas de maillage à
l’échelle de la couche mince, voir par exemple [EN93, BL96, HJ01, CCDV06].

1.3.2 Cas modèle : l’équation de Laplace

De manière habituelle, le problème a d’abord été étudié dans le cadre de l’équation
de Laplace. Il s’agit de trouver uε ∈ H10 (Ωε ) satisfaisant
− div ( a(x)∇uε (x)) = f ( x ) dans Ωε ,
où le domaine est défini par Ωε = Ωint ∪ Γ ∪ Ωext
ε , la conductivité a ( x ) est constante sur

chaque sous-domaine :
a(x) = aint dans Ωint , a(x) = aext dans Ωext
ε
,
et la donnée f est nulle dans Ωextε . La couche mince Ωε est obtenue par dilatation de
ext
la frontière Γ dans la direction normale sortante, voir Fig. 1.2. Il s’agit d’un problème de

Ωext
ε Ωext
ε

• ε •

Γ Ωint •0 α

Ωint
ε Γext
ε
Γ

Γext
ε

F IGURE 1.2 – Domaines avec couche mince (cas régulier et singulier).

transmission fortement coercif relevant de la théorie variationnelle, satisfaisant l’estima-

tion a priori suivante, avec une constante C indépendante de ε :
(uε (H1 (Ωε ) ≤ C ( f (L2 (Ωε ) .

Cas d’un domaine régulier

À l’aide des coordonnées (s, y) du repère de Frénet (s est l’abscisse curviligne sur Γ), la
ε peut être représentée par le produit cartésien Γ × (0, ε ). En posant Y =
couche mince Ωext

12
y/ε, le laplacien se développe en puissances de ε : ∆x = ∑k≥0 εk−2 Lks,Y , où les opérateurs
différentiels Lks,Y ne dépendent pas de ε. En injectant l’ansatz

 ∑ ε u (x)
 k k
 dans Ωint ,
k ≥0
uε (x) = (1.7)

 ∑ ε U (s, y/ε)
 k k
 dans Ωext
ε ,
k ≥0

dans les équations obtenues, on obtient une famille de problèmes définissant uk et U k via
une résolution alternée entre les domaines intérieur et extérieur (les opérateurs associés
sont ∆x et L0s,Y = ∂Y2 , respectivement). À l’aide d’estimations a priori (indépendantes de
ε), on montre des estimations optimales du reste du développement.
L’étude des premiers termes permet d’identifier (et de justifier) des conditions aux
limites approchées rendant compte de l’effet de la couche mince à un ordre plus ou moins
+ ,
élevé. Les trois premiers problèmes approchés s’écrivent − div a(x)∇viε (x) = f (x) dans
Ω (i = 0, 1, *2), avec la condition suivante sur Γ : v0ε = 0 (ordre 0), v1ε + γε∂n v1ε = 0 (ordre 1),
)
κ (x)
1 + 2 ε2 v2ε + γε∂n v2ε = 0 (ordre 2). Le coefficient γ = aint /aext est le ratio des valeurs
de la fonction a(x) entre les deux sous-domaines, et κ (x) la courbure au point x ∈ Γ.
Par comparaison du développement asymptotique de uε avec ceux des viε , on peut
montrer les estimations suivantes en norme d’énergie
) *
(uε − viε (H1 (Ωint ) = O εi+1 .

Cas d’un domaine à coin

Le cas où la frontière Γ n’est pas régulière est beaucoup plus délicat. Dans [CCDV06],
on a considéré un domaine présentant une singularité conique d’ouverture α, cf. F IG. 1.2.
Dans la suite, on suppose que le quotient π/α n’est pas rationnel, pour des raisons tech-
niques (des termes supplémentaires en log ε peuvent apparaître dans le cas général). La
dérivation d’un développement asymptotique comme dans le cas régulier n’est plus pos-
sible en raison des singularités de coin apparaissant dans les termes successifs. Une ana-
lyse de régularité permet de comprendre ce phénomène : le problème intérieur présente
les singularités du problème de Dirichlet, correspondant aux exposants π/αZ, alors que
le problème de transmission a ses singularités propres (pour lesquelles les exposants ne
sont pas connus analytiquement). De ce fait, la résolution alternée entre les domaines in-
térieur et extérieur ne permet pas de capturer les singularités réelles du problème résolu
par uε : le coin doit être vu de l’intérieur et de l’extérieur. Des situations similaires avaient
été rencontrées pour d’autres problèmes, cf. [CD96, BMNP01].
On a montré qu’il était possible de construire un développement (en puissances non-
entières de ε) à l’aide d’une correction au coin par des profils adaptés. Ces profils sont
définis dans le domaine non-borné H∞ obtenu par blow up au voisinage du coin 0 :

Ωε
H∞ = lim .
ε →0 ε

Il s’agit d’un secteur angulaire d’ouverture α, avec couche d’épaisseur 1, voir F IG. 1.3.

13
 Gext

H∞,ext
0 α
H∞,int

F IGURE 1.3 – Le domaine infini H∞ .

Les premiers termes du développement s’écrivent dans Ωint :

uε (x) = u0 (x) + ελ χ(x)W λ ( xε ) + · · ·

où u0 est la solution du problème limite (comme dans le cas régulier), λ = π/α le premier
exposant de singularité, χ une fonction de troncature de localisation au coin, et W λ est
défini par ) *
− div A(X)∇W λ (X) = 0 dans H∞ ,

avec la condition à l’infini W λ (X) ∼ Rλ sin(λθ ), en coordonnées polaires. La fonction A

vaut aint dans H∞,int , aext dans H∞,ext .
Insistons sur le fait que la construction des profils est nettement plus délicate ici que
dans le cas du problème du “coin arrondi”. Le domaine ne possède en effet pas la même
homogénéité car il mêle un domaine conique avec des “bandes” cartésiennes. Cette dif-
férence – apparemment technique – soulève des questions non triviales pour l’obtention
du développement des profils en l’infini car l’outil de la transformée de Mellin doit être
adapté.
Le développement obtenu permet d’analyser la performance de la condition d’impé-
dance d’ordre 1 dans le cas non-régulier. Dans [Via05], on a montré l’estimation suivante
) *
(uε − v1ε (H1 (Ωint ) = O εmin(3,2λ) ,

qui fait apparaître une perte d’ordre par rapport au cas régulier pour des angles ren-
trants (α > π). Ce résultat peut toutefois être considéré positivement, puisqu’il valide
la convergence des conditions aux limites approchées pour des domaines non-réguliers
(certes avec un taux moins élevé), situation dans laquelle elles sont utilisées fréquemment
dans les applications.
La question de modifier la condition d’impédance pour retrouver une convergence
en ε2 est toujours ouverte.

1.3.3 Flexion d’une plaque mince avec raidisseur

Avec L. Rahmani, on a généralisé dans [RV08, RV10] ces résultats au cas d’une plaque
en flexion munie de raidisseurs minces sur son bord. Le contexte géométrique est celui

14
 Ωext
ε Γext
ε

O1 ε Ωext
ε
O
• • • 2
Γ
Γ Ωint
Γ0
Ωint
Γ0
ε Γext
ε

F IGURE 1.4 – The plate Ωint stiffened by Ωext

ε .

de la figure 1.4 (cas régulier, à gauche) : le domaine Ωint représente la plaque, et Ωext
ε

le raidisseur d’épaisseur ε. On se place dans le cadre du modèle de Kirchhoff-Love qui

décrit le déplacement en flexion de la plaque uε ∈ H2 (Ωε ) :


 Dint ∆2 uε = f dans Ωint ,





 Dext ∆2 uε = 0 dans Ωext ε ,

[ M(uε )] = [ T (uε )] = 0 sur Γ, (1.8)





 M(uε ) = T (uε ) = 0 sur Γext ε ,



uε = ∂n uε = 0 sur Γ0 ,

la condition d’encastrement sur Γ0 est ajoutée pour éviter la prise en compte – technique
– des déplacements rigides. Les opérateurs M et T sont dits moment de flexion et force de
cisaillement, respectivement :
- + ,.
M = D ∆ + (1 − ν) 2n1 n2 ∂12 − n21 ∂22 − n22 ∂21 ,
- + ,.
T = D ∂n ∆ + (1 − ν)∂τ (n21 − n22 )∂12 + n1 n2 (∂22 − ∂21 ) .

Le coefficient D, constant sur chaque sous-domaine, s’exprime en fonction du module

d’Young E et du coefficient de Poisson ν (caractéristiques du matériau composant chaque
partie) sous la forme
2E
D= .
3(1 − ν2 )
Pour prendre en compte le fort différentiel entre les matériaux composant la plaque et le
raidisseur, on considère le cas où

Eint = O(1) et Eext = O(ε−1 ),

le raidisseur est ainsi à la fois très mince et très raide. Le problème (1.8) est associé à une
forme bilinéaire coercive sur H2 (Ωε ). Ici encore, on a une estimation a priori uniforme
en ε.

Cas d’un raidisseur entourant un domaine régulier

Dans [RV08], on considère la situation de gauche sur la figure 1.4 : le raidisseur est pré-
sent sur la totalité du bord Γ de Ωint , supposé régulier. On a généralisé les techniques utili-
sées pour le laplacien au problème de flexion (1.8). La technicité est bien plus grande, mais

15
on peut encore construire un développement asymptotique complet de la forme (1.7). Il
permet d’identifier et de valider, de la même manière que dans [Via05], des conditions
aux limites approchées. Le première d’entre-elles s’écrit

M ( v ε ) + Q0 ( v ε ) = 0 et T (vε ) + P0 (vε ) = 0 sur Γ,

où les opérateurs P0 et Q0 sont donnés par

- 2
+ ,.
Q0 = − Dext 2(1 − νext )∂s (∂s ∂n − c(s)∂s ) − (1 − νext )c(s) ∂2s + c(s)∂n ,
- 2
+ , .
P0 = − Dext (1 − νext )∂2s ∂2s + c(s)∂n + 2(1 − νext )∂s [c(s) (∂s ∂n − c(s)∂s )] ,

la variable s désignant l’abscisse curviligne sur Γ et ∂n la dérivée normale extérieure. Le

problème ainsi obtenu pour vε est associé à une forme bilinéaire coercive sur H2 (Ωint ). Il
est donc bien posé, et on a l’estimation suivante :

(uε − vε (H2 (Ωint ) = O(ε).

Cas d’un raidisseur sur une portion du bord

Dans [RV10], on s’intéresse au cas où le raidisseur n’est présent que sur une partie
du bord Γ, situation qui s’apparente au cas irrégulier de la figure 1.2 (l’angle à prendre
en compte est alors α = π). Le comportement singulier du problème limite ε = 0 est
dû au changement de conditions au bord aux points de raccord O1 et O2 , cf. F IG . 1.4.
On se restreint au cas où le module d’Young est constant de l’ordre O(1) dans chaque
sous-domaine (voir § 1.5 pour le cas – plus réaliste – où Eext = O(ε−1 )). Par ailleurs, pour
simplifier les calculs, on suppose que le raidisseur est situé sur une portion plane du bord.
En raison des singularités du problème limite, une correction de type profil doit ici
encore être apportée pour construire le développement asymptotique. Le domaine mo-
dèle dans lequel sont définis les profils est un demi-plan muni d’une demi-bande infinie
d’épaisseur 1, voir F IG . 1.5. La construction de ces profils est plus délicate pour cet opé-

Gext
G1 1 H∞,ext
G0 • G
O

H∞,int

F IGURE 1.5 – Le domaine H∞ dans le cas d’un raidisseur partiel.

rateur que pour le laplacien. L’ensemble des exposants de singularité du problème limite
est donné par – voir [Kon67, BR80, Gri85, Dau88, Nic93]
/ 0
1
S= + k ± iη0 ; k ∈ Z ,
2
où η0 est défini par
 34 
52
1 4 4
η0 = log  + −1 .
π (1 − ν ) (3 + ν ) (1 − ν ) (3 + ν )

16
Pour un exposant de singularité λ ∈ S donné, on note sλ la fonction singulière associée.
Le profil V λ résout alors le problème


 Dint ∆2 V λ = f dans H∞,int ,





 Dext ∆2 V λ = 0 dans H∞,ext ,



 [ M(V λ )] = [ T (V λ )] = 0 sur G,
(1.9)

 M ( V λ ) = T (V λ ) = 0 sur G ext ,





 V λ = ∂n V λ = 0 sur G0 ∪ G1 ,




V λ ∼ sλ à l’infini.

On peut montrer que ce problème est bien posé dans l’espace à poids
%
V ∇V
V= V; 2
∈ L2 (H∞ ), ∈ L2 (H∞ ), ∂α V ∈ L2 (H∞ ) pour |α| = 2,
2x3 2x3
8
et V = 0 ; ∂ν V = 0 sur G+ ∪ G− .

À l’aide d’une transformation de Mellin (polaire dans H∞,int , cartésienne dans H∞,ext ), on
peut obtenir le développement suivant pour V λ (écrit ici dans le domaine intérieur) :

V λ = sλ + ∑ V λ,µ + O( R− P ), lorsque R → ∞,
µ ∈S
− P<Re µ<Re λ

où V λ,µ sont donnés en coordonnées polaires par V λ,µ (r, θ ) = r µ ∑!, finie φ! (θ ) log! r avec
des fonctions φ! régulières (tous les exposants de S n’apparaissent en fait pas dans la
somme).
Le développement asymptotique de la solution uε du problème (1.8) prend alors la
forme suivante : pour tout entier K (fixant la précision-cible),
9 :
K 2 K
k,N −k λ ri
uε,int (x) = ∑ εk uint (x) + ∑ χ (ri ) ∑ ε k
∑ ck,λ [log ε] ε λ
Vint ( ε , θi )
k =0 i =1 k =0 λ ∈S( N − k )

K
+ rε,int ( x ),

9 :
K 2 K
uε,ext (x) = ∑ε k k,N −k
Uext ( x, yε ) + ∑ χ (ri ) ∑ ε k
∑ ck,λ [log ε] ε λ λ ri
Vext ( ε , θi )
k =0 i =1 k =0 λ ∈S( N − k )

K
+ rε,ext ( x ),
k,N −k
où x = ( x, y) et (ri , θi ) sont les coordonnées polaires centrées en Oi . Les termes uint et
k,N −k N −k
Uext sont réguliers jusqu’aux bord, plats au voisinage des points Oi comme ri . Les
scalaires ck,λ [log ε] présentent une dépendance polynomiale en ε. Enfin, on a les estima-
tions suivantes (optimales) pour les restes d’ordre K :
K
(rε,int (H2 (Ωε ) = O(εK logK ε).

17
1.4 Interaction entre perturbations
Considérons le problème simple suivant : soient Ω0 et ω deux domaines bornés conte-
nant l’origine. Pour ε > 0 petit, on pose

Ωε = Ω \ εω.

Le domaine Ωε obtenu est une perturbation du domaine initial Ω0 par une petite inclu-
sion de taille ε et de forme ω, cf. F IG . 1.6 (gauche).

d
x−
ε•
2ηε
•0 •
Ωε εω Ωε xε +εω
− −
0 •x+
ε
x+
ε + εω +

Γ Γ

F IGURE 1.6 – Petites inclusions.

On considère le problème modèle : trouver uε ∈ H1 (Ωε ) satisfaisant



 −∆uε = f dans Ωε ,

uε = 0 sur Γ = ∂Ω0 , (1.10)



∂n uε = 0 sur ε∂ω.

Le problème est très similaire à celui présenté au paragraphe 1.2, le domaine “explosé”
étant ici le domaine extérieur H∞ = R2 \ ω.
Dans [DV05], on a décrit l’influence d’une telle perturbation au bord, i.e. l’origine 0
est ici supposée située au bord du domaine Ω0 , sur l’énergie de Dirichlet
;
1
j(ε) = − |∇uε |2 dx.
2 Ωε

Lorsque le bord du domaine est plan au voisinage de l’inclusion, un développement

asymptotique complet peut être construit. Les profils mis en jeux sont définis dans le
demi plan perturbé à l’échelle 1 : H∞ = R2+ \ ω. On montre alors la formule asympto-
tique
j(ε) = j(0) + ε2 AH∞ |∇u0 (0)|2 + O (ε2 ), (1.11)

où le scalaire AH∞ est une quantité géométrique ne dépendant que de la forme du do-
maine modèle ω.
Le cas d’un bord courbe est plus délicat puisque le domaine Ω0 n’est alors plus lo-
calement stable par l’homothétie X = x/ε. Dans [DV07], on a validé un développement
asymptotique à l’ordre 2 qui permet de justifier encore la formule (1.11).
Ces développements sont fréquemment utilisés dans les calculs numériques pour
l’optimisation de forme, via la méthode dite optimisation topologique, voir [SŻ99, LS00].
Il s’agit d’évaluer la pertinence d’effectuer une inclusion en un point du domaine dans
le but de faire décroître une certaine quantité. Une formule asymptotique du type (1.11)

18
permet de répondre à cette question 2 . Dans la pratique, le calcul est effectué pour tous
les points de Ω0 situés sur une grille, et un trou est percé en chacun de ceux où la formule
asymptotique prédit que c’est favorable. Toutefois, on peut s’interroger sur le fait que
les contributions des différentes perforations ont été traitées individuellement. Lorsque
deux inclusions sont proches, un couplage peut avoir lieu, et avoir une influence sur
l’asymptotique. Cette remarque, ainsi que d’autres motivations issues des applications en
mécanique cf. § 2, nous ont conduits à étudier les interactions entre inclusions proches.
Il est aisé d’obtenir quelques premiers résultats en exploitant l’analyse précédente
du cas d’une seule inclusion. Considérons deux inclusions εω + , εω − distantes de 2ηε
(centrées en x− +
ε et xε , respectivement), cf. F IG . 1.6 (droite).

Cas d’inclusions éloignées

Si ηε = O(1), on peut justifier qu’un développement asymptotique à l’ordre 2 est
fourni par la superposition des contributions individuelles des deux inclusions :
' +
(
− x−x− + )V + ( x−xε ) + O( ε2 ),
u ε ( x ) = u0 ( x ) + ε χ ( x − x −
ε ) V ( ε
ε
) + χ ( x − x ε ε (1.12)

où la fonction χ est une troncature qui vaut 1 près de 0 et les profils V ± sont définis par


 −∆V ± = 0 dans R2 \ ω ± ,

∂n V ± = −∇u0 (0) · n sur ∂ω ± , (1.13)



V± → 0 à l’infini.

Cas d’inclusions aussi proches que petites

Si ηε = O(ε), on peut considérer les deux inclusions comme formant un seul motif
et raisonner comme s’il n’y avait qu’une seule inclusion. Le développement (1.11) n’est
plus valable, mais on doit écrire

uε (x) = u0 (x) + εχ(x)W ( xε ) + O(ε2 ), (1.14)

où W résout un problème similaire à (1.13) dans le plan privé des deux inclusions à
l’échelle 1.

Cas d’inclusions modérément proches

Dans [BNDTV09], on s’est intéressé au cas intermédiaire où les deux inclusions sont
distantes de 2ηε = εα , avec 0 < α < 1.
On a montré qu’un développement asymptotique complet pouvait être construit en
puissances entières de ε et εα , de la forme (le développement est ici écrit sans fonction de

2. Dans le cas de l’énergie de Dirichlet, la formule (1.11) laisse penser que la réponse est la même quelle
que soit la position du trou, i.e. de l’origine 0. En fait, il faut garder à l’esprit que les problèmes d’optimisation
de forme comportent généralement des contraintes, par exemple de volume et qu’il faut donc prendre en
compte la perte de volume entre Ω0 et Ωε .

19
troncature, les traces sur le bord ∂Ω0 étant relevées par les termes v p+αq ) :
' + −, + +,
(
uε (x) = u0 (x) + ε V0− x−εxε + V0+ x−εxε
' ' + −, + +,
( ' + −, + +,
((
+ ε εα Vα− x−εxε + Vα+ x−εxε + ε V1− x−εxε + V1+ x−εxε
' + −, + +,
(
+ ∑ ε p+αq v p+αq (x) + ε ∑ ε p+αq Vp−+αq x−εxε + Vp++αq x−εxε rεn (x),
( p,q)∈Kn ( p,q)∈Kn

avec
< =
Kn = ( p, q) ∈ Z2 | p ≥ 0, q ≥ − 32 p + 2, q ≥ − p et p + αq ≤ n .
On a l’estimation suivante du reste :

(rεn (H1 (Ωε ) = O (εn ).

Les termes du développement correspondent aux contributions suivantes :

4 profils en variable rapide Vp±+αq : ils prennent en compte le développement de Tay-
lor de u0 en 0, ainsi que leur interaction mutuelle. V ± relève la trace de V ∓ sur le
bord de l’autre inclusion xε± + ε∂ω ± ,
4 correcteurs en variable lente v p+αq : ils relèvent la trace des profils sur le bord exté-
rieur ∂Ω0 .
En particulier, on remarque que le début du développement coïncide avec (1.12), mais le
reste est détérioré. En effet, le premier reste rε1 est d’ordre εmin(1+α,3−2α) + ε2 . On peut
distinguer deux régimes :
4 pour α < 2/3, les inclusions sont peu proches l’une de l’autre. Le terme principal
du reste rε1 provient du développement de Taylor de u0 en l’origine 0,
4 pour 2/3 < α < 1, les inclusions sont plus proches. Le reste est principalement dû
à l’interaction des profils V0− et V0+ .
On a aussi traité le cas où l’une des deux inclusions est proche du bord du domaine,
problème similaire à ceux étudiés dans [BV00, CV03, BHB06, AKLZ10]. Une extension
aux équations de l’élasticité linéaire est présentée dans la partie suivante, donnant lieu à
des calculs numériques en mécanique de la rupture.

1.5 Perspectives
De nombreuses pistes restent à explorer. Dans le cadre de l’interaction entre inclu-
sions, le cas de deux inclusions très proches, i.e. de taille ε et distantes de εα avec α > 1
est en cours. Trois domaines limites sont à considérer :
4 le domaine sans inclusion :
Ω0 = lim Ωε ,
ε →0

4 le domaine “explosé” à l’échelle ε, consistant en les deux inclusions à l’échelle 1,

tangentes l’une à l’autre :
Ωε
H∞ = lim ,
ε →0 ε

20
 4 le domaine “explosé” à l’échelle εα , qui est une bande d’épaisseur 1 (les droites la
définissant correspondant aux plans tangents des inclusions) :

Ωε
B∞ = lim ,
ε →0 εα
Par ailleurs, on peut poursuivre l’étude de l’asymptotique du problème (1.8) dans
le cas non-régulier avec un module d’Young en ε−1 dans le raidisseur. Les problèmes
définissant les profils dépendent alors eux-mêmes du paramètre, et l’on doit gérer un
double développement (ε → 0 et X → +∞). La question de l’efficacité des conditions
aux limites approchées doit aussi être quantifiée dans ce cas, comme il a été fait pour le
laplacien.
En outre, même dans le cas de l’équation de Laplace, il serait particulièrement utile
pour les applications de déterminer une modification des conditions aux limites appro-
chées au voisinage d’un point singulier, afin de retrouver un ordre de convergence plus
élevé. Cette question est difficile et pourrait être envisagée d’un point de vue numérique,
en lien avec les développements à l’infini des profils.

21
22
Chapitre 2
Méthodes numériques pour la prise en
compte de micro-défauts en
mécanique

2.1 Contexte
Ce travail s’inscrit dans le cadre du projet ANR Jeunes Chercheurs M ACADAM et plus
particulièrement dans une collaboration avec D. Brancherie, mécanicienne à l’UTC. Il
concerne la mécanique de la rupture, domaine particulièrement actif du point de vue nu-
mérique. De nombreuses méthodes pour décrire l’initiation et la propagation de fissures
ont été développées, voir par exemple [WI90, Oli95, DBDI05, CLP+ 07, KMZ08, BFM08,
CFM09]. La très utilisée méthode XFEM consiste en une méthode d’éléments finis où l’es-
pace d’approximation est enrichi de fonctions indicatrices pour traiter les discontinuités
dues à la fissuration, et de fonctions singulières en front de fissure.
Notre objectif est d’étudier le comportement à rupture des structures comportant des
petits défauts de surface. Ces perturbations ont un effet important sur l’apparition des
zones de fissuration. Afin de ne pas mailler la géométrie réelle – qui en nécessiterait une
description fine – on souhaite utiliser un modèle asymptotique (voir § 2.2) dont l’approxi-
mation peut être calculée sur un maillage assez grossier. L’information obtenue est alors
couplée à un modèle d’endommagement (cf. § 2.3). Au paragraphe 2.4, on s’intéresse
plus particulièrement à l’approximation des profils apparaissant dans la représentation
asymptotique.
Ce chapitre fait la synthèse des articles [DV07, BNDTV09, BNDHV10b] et [BDVV08,
BNBD+ 10], en collaboration avec V. Bonnaillie-Noël, D. Brancherie, M. Dambrine, F. Hé-
rau, S. Tordeux et P. Villon.

2.2 Représentation asymptotique

On considère un matériau soumis à des contraintes mécaniques, et présentant des
aspérités, cavités ou inhomogénéités de surface. Lorsque les sollicitations augmentent,
on souhaite déterminer où vont avoir lieu les fissurations, ainsi que leur évolution. La

23
géométrie comporte deux échelles et, comme on l’a déjà souligné, un calcul numérique
direct requiert un raffinement de maillage important au voisinage des défauts, qui peut
alourdir fortement la durée de simulation. On préfère recourir à un modèle simplifié
faisant intervenir le matériau sans défaut.
Dans [DV07, BDVV08], on a étendu les résultats asymptotiques présentés dans le pa-
ragraphe 1.4 aux équations de l’élasticité linéaire. Précisément, on considère les équations
de Navier dans un domaine borné Ωε présentant un défaut en son bord :


 −µ∆uε − (λ + µ)grad div uε = f sur Ωε ,

uε = u D sur Γ D , (2.1)



σ (u) · n = g sur ∂Ωε \ Γ D .

Ici encore, le domaine Ωε est obtenu à partir du domaine limite Ω0 à l’aide d’un motif
auto-similaire centré en 0, donné par le domaine semi-infini H∞ : Ωε = Ω0 \ εH∞ , cf.
Fig. 2.1. Pour des raisons techniques, on suppose que le bord de Ω0 coïncide avec un

ΓD ΓD
Ω0 H∞ Ωε

• •
0 •0 0

F IGURE 2.1 – Les domaines Ω0 , H∞ et Ωε .

segment de droite au voisinage du point 0, et que les chargements f et g s’annulent au

voisinage de ce point. Au premier ordre, uε s’écrit
' )x* ) x *(
u ε ( x ) ! u0 ( x ) − ε α1 V1 + α2 V2 , (2.2)
ε ε
où u0 est la solution sur le domaine non perturbé Ω0 , et les réels α1 , α2 sont donnés par

α1 = σ11 (u0 )(0) et α2 = σ12 (u0 )(0). (2.3)

Les profils V1 et V2 sont solutions de problèmes de Navier sans chargement distribué

dans H∞ avec des conditions de Neumann (dans le problème suivant, ! = 1 ou 2) :


 −µ∆V! − (λ + µ)grad div V! = 0 sur H∞ ,

σ (V! ) · n = G! sur ∂H∞ , (2.4)



V! → 0 en l’infini.

avec G1 = (n1 , 0) et G2 = (0, n1 ), où n1 est la première composante de la normale exté-

rieure à H∞ .
La méthode numérique proposée repose sur la représentation (2.2). Le terme limite u0
est calculé sur le domaine Ω0 (qui peut être maillé relativement grossièrement). Bien sûr
on doit encore résoudre les problèmes (2.4) pour déterminer les profils V! . Il faut veiller
à ce que ce calcul ne soit pas trop coûteux, sous peine de perdre tout le bénéfice de la
simplification du modèle. Si, dans le cas de l’équation de Laplace en dimension 2, des

24
méthodes ad hoc peuvent être employées cf. [BNDTV09], ce n’est pas le cas pour les équa-
tions de Navier. Notons aussi que les solutions de (2.4) ne sont pas connues explicitement
même dans le cas où l’inclusion est semi-circulaire.
Pour résoudre (2.4), on tronque le domaine H∞ au rayon R > 0 en posant
R
H∞ = { x ∈ H ∞ ; | x | < R }.

Nous avons pour l’instant utilisé la méthode la plus naïve, consistant à imposer une
condition de Dirichlet homogène sur la frontière artificielle Γ R = ∂H∞ R \ ∂H . Cette
∞
condition peut être considérée comme une condition aux limites approchée d’ordre 0,
voir § 2.4 pour une discussion autour de conditions d’ordre plus élevé. On pourrait en-
visager d’autres approches, vu la décroissance des profils en l’infini, comme le couplage
avec une méthode de représentation intégrale, voir par exemple [LVLH92].
La formule (2.2) fournit une approximation à l’ordre ε2 du champ de déplacement
à partir de celui dans le domaine non perturbé et des profils, qui sont intrinsèques au
motif de perturbation. Dans le cas de plusieurs inclusions, on peut étendre les résultats
obtenus dans [BNDTV09] pour l’équation de Laplace au cas des équations de Navier :
si les inclusions sont relativement proches les unes des autres – i.e. à distance εα avec
α < 1 – il suffit de superposer leurs contributions respectives (l’erreur commise est certes
détiorée, mais reste négligeable vis-à-vis de ε).
Afin de valider la stratégie numérique, on considère le cas-test de la figure 2.2, où le
matériau est sollicité en traction selon l’axe horizontal. Les deux défauts en 01 , 02 sont
de forme semi-circulaire, de rayons respectifs de 2 mm et 1 mm ; la distance 01 02 vaut
30 mm. Notons que les profils sont calculés une seule fois pour les deux inclusions, car

200 mm

u
100 mm

Ωε
E = 38 × 104 MPa, ν = 0.18

01 02

FFigure 1: Définition du problème modèle.

IGURE 2.2 – Définition du problème test.

elles partagent le même motif H∞ . Il suffit donc de superposer les contributions selon la
formule (2.2).
Le calcul des profils a été effectué avec la bibliothèque de calculs éléments finis M É -
LINA [Mar10] (voir chapitre suivant pour plus de détails). L’approximation du déplace-
ment uε a utilisé un code de calcul développé par D. Brancherie. Si la formule de super-
position (2.2) est commode du point de vue mathématique, il est apparu qu’il était plus
simple techniquement de reconstruire l’approximation de uε un peu différemment. On
a procédé à un enrichissement cinématique de la base éléments finis par les profils, la
perturbation étant prise en compte par des formules de quadrature adaptées localement.

25
Afin d’assurer la stabilité du schéma, on doit relier les degrés de liberté additionnels par
une stratégie “maître-esclave”. Notons que ce calcul ne nécessite pas l’évaluation des
contraintes limites α! , voir formules (2.3).
La figure 2.3 montre deux maillages : le premier, fin, discrétise le domaine perturbé Ωε
et est utilisé pour calculer la solution de référence uε à des fins de comparaison. Le second
maille Ω0 assez grossièrement et permet de mettre en œuvre la méthode numérique basée
sur la représentation asymptotique (2.2). Les solutions sont comparées sur la figure 2.4 :
l’erreur relative (calculée à l’aide de la solution de référence uε ) est de l’ordre de 0.25%.

100

90

80

70

60

50

40

30

20

10

0
0 20 40 60 80 100 120 140 160 180 200

(a) maillage fin

100

80

60

40

20

0
0 50 100 150 200
(b) maillage grossier

F IGURE 2.3 – Maillages utilisés

26
 (a) Calcul de référence sur maillage fin
(a) Calcul de référence (b) Calcul avec enrichiss

(b) Calcul avec enrichissement

Calcul de référence (a) Calcul de(b) Calcul avec
référence enrichissement
(b) Calcul (c) enrichissement
Erreur relative cinématique
cinématique
avec

Figure 3 – Déplacements ux obtenus par un calcul standard et avec enrich

relative

stresses !xx (ref. computation) stresses ! (enr

xx
100 500 100

80 400 80

60 300 60

40 200 40
(c) Carte de l’erreur relative
(c) Erreur relative (c) Erreur relative
20 100 20

F IGURE 2.4 – Déplacements

0
0 20 40
u x obtenus
60 80 100
: (a)
120
par
140
calcul
160
standard
180 200
sur
0 un maillage fin,
0 et (b)
0 20 40 60 80 1
avec enrichissement. Le graphe (c) présente la carte de l’erreur relative.
nts
gureux3obtenus par un calcul
– Déplacements standard
ux obtenus paretun
avec enrichissement,
calcul carteenrichissement,
standard et avec de l’erreur carte de l
ative (a) Calcul de référence (b) Calcul avec enric
27
(ref. computation) stresses !xx (ref. computation) stresses ! (enriched computation) log10 of stresses
energy relative error computation)
! (enriched
xx xx
xx
100 500 100 500 100 100 500 !1

!2
400 80 400
 Two-scale appro

opening and discontinuity line orientation. The discontinu

tiated at the element containing the perturbation and propa
2.3 Couplage avec un modèle à discontinuité forte
to the top edge of the specimen. The mesh objectivity of the
La représentation simplifiée (2.2) permet d’obtenir une approximation à moindre coût
on Figure
de la solution 10 :(2.1).
u du problème
ε theCe responses and
calcul est utilisé pour total
décrire dissipation
les zones de loca- evaluat
discretisations are similar.
lisation de contraintes à l’aide d’un modèle à discontinuité forte (Strong Discontinuity Ap-
duites.
proach, Lavoir
SDA), Figure 5(b)
[SOA93, Oli95, donne,
DBDI05]. en finne de
Cette méthode chargement,
vise pas à représenter la les orient
géométrie réelle des fissures, mais indique, pour chaque élément du maillage que ren-
ainsi
contre laque
fissure,l’ouverture de(rien
une direction de fracture fissure. La perturbation
n’assure la continuité de la ligne globale). de plus gra
σ1 (maximum principal stress)
La description est moins fine qu’avec d’autres types de méthodes, mais présente l’avan- σ1 (max
l’origine 100 de la fissure conduisant à la rupture. Le développeme
tage d’être peu sensible au maillage.
212

210
100

d’un déchargement
Comme dans les méthodes XFEM,
80 élastique de
la base nodale desl’ensemble du
éléments finis est enrichie
fonctions de Heavyside pour prendre en compte les discontinuités de fissures. De ce fait,
domaine,
de la p 208 80
206

pas une
l’espace deuxième
d’approximation
60 fissuration.
ne coïncide pas avec celui utilisé pour le calcul de u , où l’en- 204ε 60

richissement utilisait les profils pré-calculés. On a donc adopté une stratégie 202
de transfert
40
de champ entre les deux discrétisations, basée sur un critère énergétique. 40
200

198
20 !1 (maximum principal stress) 20
100 196
100
0 194
210 0
0 20 40 60 80 100 120 140 160 180 200 0 20 40 60 80
80 192 80

60 (a) Structured mesh 205

60 (b) Un
40 200
40

20

Figure 8. Maximum principal stress associated with the r

195 20

field u
0 20 40 60 80 100 120 140 160 180 200 0
0 20

F IGURE 2.5 – Contraintes principales maximales reconstruites.

(a) Contraintes principales maximales reconstruites (b) Disconti
crack opening 0.00E+00 crack o
4.35E−01
8.70E−01
1.30E+00
Figure 5 – Différence entre le champ de déplacement reconstru 1.74E+00
2.17E+00
non perturbé, contraintes principales maximales recontruites 2.61E+00
3.04E+00
3.48E+00

replacements PSfrag replacements 3.91E+00

4.35E+00
4.78E+00
5.22E+00

5 Conclusion
(a) Structured mesh
F IGURE 2.6 – Discontinuités introduites et ouverture de fissure.
(b) U

Nous avons présenté28une stratégie permettant de prend

Figure
défauts 9. comportement
sur le Discontinuity line and crack
jusqu’à opening
rupture d’une structure
 0 160 180 200 0 20 40 60 80 100 120 140 160 180 200
192 192

mesh (b) Unstructured mesh

cipal stress associated with

Les figures 2.5, 2.6 thelesreconstructed
présentent résultats obtenus dans ledisplacement
cas d’un maillage cartésien :
la ligne de fissuration s’ouvre à partir de l’inclusion la plus grosse (la seconde inclusion
ne donne pas lieu à fissuration). Sur la figure 2.7, on présente un résultat dans un cas
non structuré où l’on voit que la ligne de discontinuité reste quasi verticale malgré la
mauvaise qualité du maillage.

0.00E+00 crack opening 0.00E+00

4.35E−01 4.61E−01
8.70E−01 9.23E−01
1.30E+00 1.38E+00
1.74E+00 1.85E+00
2.17E+00 2.31E+00
2.61E+00 2.77E+00
3.04E+00 3.23E+00
3.48E+00 3.69E+00

frag replacements 3.91E+00

4.35E+00
4.15E+00
4.61E+00
4.78E+00 5.07E+00
5.22E+00 5.54E+00

mesh (b) Unstructured mesh

F IGURE 2.7 – Ouverture de fissure dans le cas d’un maillage non structuré.

ne and crack2.4opening
Conditions aux limites approchées pour les profils
On a vu plus haut que le calcul des profils – voir équations (1.13) et (2.4) – est un
point essentiel des méthodes numériques proposées. Il faut tout à la fois une approxi-
mation précise et peu coûteuse. Le fait d’utiliser une conditions aux limites artificielle de
type “Dirichlet homogène” impose de prendre de grands rayons de troncature, et donc
un temps de calcul plus important. On résume ici un travail effectué dans [BNDHV10b]
qui représente un premier pas vers l’obtention et la validation d’une condition aux li-
mites améliorée pour les équations de Navier. La recherche de telles conditions est liée
aux problèmes de couches minces, voir [EN93, BL96, HJ01], ou de parois rugueuses,
voir [APV98, JM01, ACG06], ou encore de conditions absorbantes pour les phénomènes
o methods and their coupling allowing to describe the be-
propagatifs, voir par exemple [EM77, Giv91, HR95].

ttle structures. The

Cas de main
l’équation feature of this work is that com-
de Laplace
at a macroscopic levelle problème
On considère while modèle
the suivant
models dans letake
domaine into account
extérieur H = R \ω: ∞
2

characteristics. This is achieved 




−∆Vthanks
= 0 todansanH enrichment
, ∞ of
t space with the profiles computed 


from
∂ V = g · n the
nsur ∂ω,asymptotic ana- (2.5)
t rid of a fine mesh of the real geometry. V → 0 à l’infini,

où ω est un domaine borné de R2 , n désigne le vecteur normal extérieur, et g est une

o consider more general loadings, geometries of the struc- R = H ∩ B le domaine tronqué (B
fonction vectorielle régulière donnée. On note H∞ ∞ R R

e defects. We expect no major difficulties

29 in reaching such
ioned above are designed in that spirit.
désigne la boule centrée en 0 de rayon R), et V 0 la solution du problème avec conditions
de Dirichlet : 

 −∆V 0 = 0 dans H∞ R,

∂n V 0 = g · n sur ∂ω, (2.6)



V0 = 0 sur ∂BR .

Au vu de la décroissance V (X) = O(|X|−1 ) en l’infini, on peut montrer l’estimation

(V − V 0 (H1 (H∞R ) = O( R−1 ),

d’où un choix de R très grand si l’on souhaite obtenir une bonne précision. Une méthode
pour déterminer des conditions aux limites approchées d’ordre plus élevé consiste à uti-
liser un développement précisé en l’infini du type (1.5), qui peut être obtenu dans le cas
particulier de l’équation de Laplace à l’aide de séries de Fourier :

a1 b1 a2 b2
V (X) = cos(θ ) + sin(θ ) + cos(2θ ) + sin(2θ ) + · · ·
|X| |X| | X |2 | X |2

On cherche alors une relation différentielle qui annule les premiers termes de ce déve-
loppement. Le domaine étant borné par une boule, la dérivée normale est radiale si bien
qu’un calcul simple conduit aux conditions suivantes :

4 Ordre 1 : V 1 + R∂n V 1 = 0,
3R R2
4 Ordre 2 : V 2 + 2 ∂n V
2 − 2 ∆τ V
2 = 0.

La condition d’ordre 1 est de type Robin-Fourier, celle d’ordre 2 de type Ventcel. Toutes-
deux conduisent à des problèmes variationnels bien posés, et on a les estimations

(V − V i (H1 (H∞R ) = O( R−i−1 ) (i = 0, 1, 2).

Dans [BNDTV09], des calculs numériques ont été menés pour comparer ces conditions
(une méthode basée sur l’inversion conforme z 6→ z−1 est aussi testée).

Vers une condition d’ordre 1 en élasticité

Dans le cas des équations de Navier – cf. problème (2.4) – on peut encore écrire un
développement en l’infini et chercher à améliorer la condition à imposer sur la frontière
artificielle. La condition obtenue s’écrit alors – voir [BNDHV10a]
  > ?
1 1 1
1 E  0 1 E ( 1 − ν ) 0 0
σ(V )n + V + ∆τ V1 = 0. (2.7)
R 1+ν 1−ν R 2(1 + ν)(1 − 2ν) 0 1
0 1

La forme bilinéaire associée n’est plus coercive et le problème ne relève plus de la théo-
rie variationnelle. Dans [BNDHV10b], on a étudié un cas modèle qui présente la même
difficulté : %
−∆u = f dans Ω,
(2.8)
∂n u + αu + β∆τ u = 0 sur ∂Ω,

30
avec α, β > 0. Il s’agit encore d’un problème de type Ventcel, mais le coefficient β a le
“mauvais” signe. Si Ω est le disque unité de R2 , un calcul direct à l’aide de séries de Fou-
rier montre que le problème (2.8) admet une unique solution sous la condition (nécessaire
et suffisante)
α∈/ { βn2 − n ; n ∈ N }.
Pour étudier le cas général, on note Λ l’opérateur de Dirichlet à Neumann sur ∂Ω : si v
désigne le relèvement harmonique de g dans Ω, alors Λg = ∂n v|∂Ω . Après relèvement du
second membre f , le problème (2.8) se récrit comme une équation non-locale sur le bord
∂Ω :
β∆τ u + Λu + αu = ϕ. (2.9)
De cette manière, le problème apparaît comme une perturbation pseudo-différentielle de
l’opérateur de Laplace Beltrami. La théorie de Fredholm s’applique et le problème (2.9)
admet une unique solution en dehors d’une condition sur les coefficients, du type α 7= αn ,
où (αn ) est une suite croissant vers l’infini (qui dépend de β, supposé fixé).
Dans le contexte des conditions aux limites artificielles, on considère le problème

 −∆u = 0 dans H∞ R,

∂n u + αu + β∆τ u = 0 sur ∂BR , (2.10)

 u = g sur ∂ω,

où la condition de Ventcel précédente a été utilisée comme condition “absorbante” sur

la frontière ∂BR . On a montré que le problème (2.10) est bien posé si R est suffisamment
grand. Des simulations numériques laissent penser qu’il n’est pas nécessaire de choisir R
très grand pour avoir existence et unicité d’une solution.

2.5 Perspectives
Il reste à étudier la condition (2.7) dans le cadre de l’élasticité linéaire. Les techniques
pseudo-différentielles utilisées pour le problème modèle Laplace-Ventcel doivent pou-
voir s’étendre, et on attend des résultats similaires. Une étude numérique est aussi en
cours pour valider l’ordre de la condition aux limites approchée. Le calcul des profils
ainsi amélioré sera alors intégré à la méthode numérique pour l’étude complète de fissu-
ration.
Une extension au cas tri-dimensionnel est aussi envisagée, permettant de considérer
des problèmes plus proches de la réalité. Si les aspects théoriques semblent se prolonger
dans difficulté majeure, il n’en va pas de même des codes de calcul pour lesquels un
investissement important est nécessaire.

31
32
Chapitre 3
Programmation pour le calcul
scientifique

3.1 Contexte
Les travaux de nature théorique présentés dans les deux parties précédentes trouvent
leur motivation dans la nécessité de simplifier le modèle mathématique en vue du calcul
numérique. Ainsi une grande partie de mon activité de recherche consiste en la mise
en place de codes de calcul pour tester les méthodes proposées ou valider des résultats
asymptotiques.
La plupart des calculs numériques effectués s’appuie sur la bibliothèque de calcul élé-
ments finis M ÉLINA [Mar10], développée par D. Martin. Si une partie des résultats pré-
sentés ici concernent naturellement les thèmes développés précédemment (mécanique de
la rupture, modèles asymptotiques – voir § 3.3 et § 3.4.1), d’autres calculs ont été menés,
issus de rencontres scientifiques et n’ayant d’autre lien entre-eux que la nécessité d’un
calcul robuste et précis. Ainsi, ont été abordés des problèmes aussi différents que l’ana-
lyse de la supraconductivité en domaine polygonal (§ 3.4.2), la séparation de phase dans
un mélange binaire (§ 3.4.3), l’optimisation de forme pour le matriçage cellulaire (§ 3.5.1)
ou l’utilisation de simulations numériques pour tester des conjectures dans un problème
de partitions minimales (§ 3.5.2).
On fait ici la synthèse des articles [CDMV03, Via05, BNDMV07, BDVV08, BNDTV09,
BNHV10, PV10, GMV10], en collaboration avec V. Bonnaillie-Noël, D. Brancherie, M. Cos-
tabel, M. Dambrine, M. Dauge, L. Goudenège, B. Helffer, D. Martin, M. Pierre, S. Tordeux,
P. Villon.

3.2 La bibliothèque éléments finis M ÉLINA

Développée au sein de l’IRMAR et de l’ENSTA, M ÉLINA s’inscrit dans le cadre des bi-
bliothèques universitaires pour le calcul par éléments finis comme FreeFem, GetFEM++,
OFELI, etc. C’est la version Fortran 77 (1986-2010) qui a été utilisée ici, mais signalons
qu’une version C++ totalement refondue est déjà en partie opérationnelle, et sera pro-
chainement disponible.

33
 La caractéristique principale de M ÉLINA est sa modularité, qui permet à l’utilisateur
de personnaliser cet outil en un logiciel dédié à une tâche particulière. Les problèmes
sont définis directement à l’aide de leur formulation faible, et ne sont pas choisis parmi
une liste pré-définie. De nouveaux intégrands peuvent être introduits, de même que de
nouveaux éléments finis ou d’autres méthodes de résolution. Ainsi, il est possible de
traiter une gamme très large de problèmes.
Parmi les bibliothèques éléments finis, M ÉLINA est une des rares à proposer des élé-
ments de haut (voire très haut) degré (jusqu’à 64 pour les éléments tensoriels). Cette fonc-
tionnalité a été introduite lors de mes travaux de thèse, et a été utilisée sur de nombreux
problèmes. Ce chapitre expose les résultats obtenus avec cet outil dans différents do-
maines d’application.
M ÉLINA a été utilisé dans bien d’autres applications (en électromagnétisme en parti-
culier). Pour plus de détails, on renvoie au site web [Mar10], section “Pressbook”.

3.3 Calcul des profils pour la mécanique numérique

On a décrit dans le chapitre 2 une méthode numérique pour le calcul à rupture d’une
structure fragile. On détaille ici quelques aspects liés au calcul numérique des profils.
Le calcul des profils en domaine non borné (voir problèmes (2.4) ou (2.5) pages 24
et 29, respectivement) est un ingrédient essentiel des méthodes asymptotiques propo-
sées. Les profils sont intrinsèques à la géométrie des inclusions, ils peuvent être consi-
dérés comme le résultat du couplage entre le motif local de perturbation et la géométrie
sectorielle à l’infini (dans le cas d’une inclusion, l’angle du secteur vaut 2π). De fait, il
faut prendre en compte de manière précise à la fois la géométrie de la perturbation, et le
comportement en l’infini.
Comme il a été signalé plus haut, le comportement à l’infini est généralement pris
en compte par une condition sur une frontière circulaire fictive au rayon R. Quant à la
géométrie locale, elle nécessite une approximation fine qui requiert un grand nombre
d’éléments si l’on se restreint à une approximation de degré 1 ou 2.

(a) Maillage du domaine tronqué à (b) Maillage pour la méthode d’in-

R = 10. version.

F IGURE 3.1 – Maillages utilisés pour le calcul des profils.

34
 Sur la figure 3.1 (gauche), on représente un maillage utilisé pour une méthode de
conditions aux limites approchées avec R = 10 (dans le cas d’une inclusion circulaire).
La représentation géométrique est réalisée avec des polynômes Q8 . L’utilisation d’un
maillage de haut degré permet la réduction du nombre d’éléments (ici 128) tout en as-
surant une excellente précision géométrique.
Dans le cas de l’équation de Laplace (2.5), on a testé une méthode particulière basée
sur l’inversion ϕ : z 6→ z−1 . L’équation −∆V = 0 reste inchangée, la donnée de Neu-
mann est transportée, et la condition à l’infini devient une condition de Dirichlet ponc-
tuelle en l’origine. L’intérêt est évidemment d’être ramené à un domaine borné, dont un
maillage Q8 est représenté sur la figure 3.1 (droite).
Sur la figure 3.2, on représente l’erreur commise (par rapport à la solution exacte,
connue explicitement dans le cas d’une inclusion circulaire) dans le calcul du profil en
fonction du degré d’interpolation de la méthode de Galerkin (le degré d’approximation
géométrique, quant à lui, reste égal 8). Les conditions aux limites approchées sont celles
données page 30.

0
10

−1
10

−2
Relative maximum error

10

−3
10

−4
10

−5
10

Artificial boundary method (order 0)

−6 Artificial boundary method (order 1)
10
Artificial boundary method (order 2)
Inversion method
−7
10
1 2 3 4 5 6 7 8
Degree of interpolation [prop. to sqrt(number of DOF) ]

F IGURE 3.2 – Comparaison inversion / conditions aux limites approchées.

On observe un phénomène de saturation pour la méthode utilisant des conditions

aux limites approchées, à un niveau d’autant plus bas que la condition est précise. La
méthode d’inversion est plus précise, mais est difficilement généralisable à d’autres cas,
par exemple en élasticité.
Pour le couplage représentation asymptotique / modèle à discontinuité forte (voir 2.3
page 28), les valeurs des profils sont utilisées pour l’enrichissement cinématique de l’es-
pace d’approximation. Bien sûr, les maillages ne coïncident pas, et un transfert est néces-
saire sur la grille utilisée par la méthode à discontinuité forte. Si les profils ont été calculés
sur une grille régulière (comme le maillage polaire à gauche de la figure 3.1), la projection
n’est pas difficile. Dans le cas général, une méthode utilisant un bucket sort a été mise en
œuvre.

35
3.4 Utilisation des éléments finis de haut degré
Dans cette partie, on présente trois situations où l’utilisation d’éléments finis de haut
degré a permis un calcul précis en temps court, là où l’emploi d’éléments P1 aurait né-
cessité un temps de calcul prohibitif.

3.4.1 Problèmes de singularités

Dans [CDMV03], nous avons utilisé pour la première fois les possibilités offertes par
les éléments de haut degré pour le calcul de modes propres pour les équations de Max-
well harmoniques. Si l’on considère la formulation variationnelle
; ;
rot u · rot v dx = λ2 u · v dx,
Ω Ω

la discrétisation directe nécessite l’emploi des éléments d’arête introduits dans [Néd80].
Une méthode classique pour permettre l’utilisation d’éléments finis nodaux consiste en la
régularisation par un terme en div u div v, cf. [HL96]. Toutefois, dans le cas d’un domaine
Ω non-régulier, il a été montré dans [Cos91] que l’espace d’approximation discret n’est
pas dense dans l’espace continu. Une méthode introduisant un poids lié à la distance aux
points de singularité dans le terme de régularisation est proposée dans [CD02].
La figure 3.3 montre les résultats d’approximation obtenus pour le calcul de la pre-
mière valeur propre dans un domaine bi-dimensionnel avec un coin rentrant pour diffé-
rents degrés d’approximation géométrique et fonctionnel. Weighted Regularization of Maxwell Equations 7

0 1
10

2
2
3
-2 3
10 4 5
6 7 8 9 10
4

Geo. deg. 1
5
-4
10 Geo. deg. 2
6

Geo. deg. 3
7
Geo. deg. 4 8
-6
10 8 9 10
Geo. deg. 5
9
Geo. deg. 6 9
10
-8 10
10 2 3 4
10 10 10

(a) Maillage du domaine. (b) Erreur relative sur la première valeur propre.
Fig. 6. Relative errors for the 1st eigenvalue

F IGURE 3.3 – Calculs de régularisation10à0 poids pour les équations de Maxwell.

1

-2 2
3 4
10 2
5 6 7 8 9 10

Ici encore, le fait de travailler avec des maillages de haut degré (la frontière du do-
3
maine discret est polynomiale par morceaux) 10
-4
permet de retrouver une approximation
4
aussi bonne que prévue par la théorie. Comme souvent, la mauvaise représentation
Geo. deg. 1
5
5
6 de la
7 8 9 10
-6
géométrie au niveau discret est un facteur limitant. Notons aussi que les points
10 Geo. deg. 2
6 d’interpo-
lation sont choisis aux abscisses de Gauss-Lobatto
-8
pour éviter
Geo. deg. 3
un phénomène 7de Runge
8 9 10
10 7 8
(qui apparaît à partir des degrés 6-8) détériorant significativement les résultats. 9
Geo. deg. 4

-10 Geo. deg. 5

10 10
10

36 Geo. deg. 6

-12
10 2 3 4
10 10 10

Fig. 7. Relative errors for the 2nd eigenvalue

 Mentionnons enfin que le raffinement géométrique au voisinage du coin, lié à une
interpolation de haut degré, permet de retrouver des résultats de convergence similaires
à ceux existant pour des domaines réguliers, voir [Mel02]. Cette stratégie a été utilisée
dans [Via05] pour valider l’ordre de convergence des conditions d’impédance dans les
domaines à coins.

3.4.2 Problèmes d’oscillations

Dans [BNDMV07], il s’agit de calculer les modes propres de l’opérateur de Schrödin-
ger avec champ magnétique en limite semi-classique dans un domaine borné Ω :

Ph φh = λh φh avec Ph = −(h∇ − i A)2 (+ conditions naturelles).

√
Lorsque h tend vers 0, la première fonction propre se concentre à l’échelle h près du
point de courbure maximale (cas où Ω est régulier) ou du sommet d’angle minimum
(cas polygonal). Des oscillations à l’échelle h apparaissent aussi, rendant les calculs nu-
mériques extrêmement délicats. Une approche standard avec éléments finis P1 et raffi-
nement échoue. En effet, la phase des fonctions propres oscille fortement dans tout le
domaine (pas seulement là où le module se concentre) et il est indispensable de capturer
toutes ces oscillations pour obtenir une approximation raisonnable du module, pourtant
non-oscillant... On a montré qu’une stratégie en p-version (augmentation du degré des
éléments finis, plutôt que réduction du pas de maillage) fournit d’excellents résultats.
Les figures 3.4 et 3.5 montrent la première fonction propre obtenue avec environ 4000
élements Q1 et 4 éléments Q20 (1764 degrés de liberté).

F IGURE 3.4 – Première fonction propre avec 3969 éléments Q1 (module et phase).

F IGURE 3.5 – Première fonction propre avec 4 éléments Q20 (module et phase).

37
 On voit ainsi que la méthode numérique Q1 fournit un résultat qui n’est pas celui
attendu : la fonction approchée correspond à un mode propre faiblement oscillant, mais
d’un niveau plus élevé que la fonction oscillante obtenue par calcul Q20 .

3.4.3 Problèmes d’interface

Dans [GMV10], on présente des résultats de simulations numériques pour la sépara-
tion de phase d’un mélange binaire. Le phénomène est modélisé par une fonction d’état u
(valant 1 et −1 pour les deux phases pures) satisfaisant l’équation de Cahn-Hilliard dans
un domaine borné Ω : 
 ∂ u = ∆w dans Ω,
 t

w = ψ(u) − ε2 ∆u dans Ω,



∂n u = ∂n w = 0 sur ∂Ω,
où le terme non linéaire ψ(u) est donné par
1+u
ψ(u) = αu + β log (α, β ∈ R ).
1−u
Dans beaucoup de travaux, on remplace cette expression par un développement limité
du logarithme, évitant ainsi les difficultés liées aux singularités en ±1, dont le modèle à Logarithmic energies under polynomial non!linearity.

l’ordre 3 s’écrit !1

ψ(u) = u3 − u.
!2

!3

On a utilisé une méthode d’éléments finis mixte pour éviter la programmation d’éléments
!4
Log10 (Energy2n ! Energylog)

conformes-H2 (par exemple nodaux C 1 , présents dans très peu de codes numériques).
!5

!6
Parmi les aspects abordés, on en distingue ici deux. f !7 4
f6
!8 f8
Influence du remplacement du logarithme par un polynôme dans le terme non linéaire !9
f10
f12
hal-00461231, version 1 - 4 Mar 2010

f14

On considère le cas mono-dimensionnel Ω = (0, 1). La figure 3.6 représente l’erreur

!10
f16

commise par une approximation polynomiale du logarithme (l’abscisse est le temps, l’or-
!11
0 100 200 300 400 500 600
Time iterations
700 800 900 1000

donnée le logarithme décimal Figure

de l’erreur en norme L2 ) pour différents degrés. Il apparaît
10: Polynomial energies versus logarithmic energy.

en particulier que l’approximation par un polynôme de degré 3 est peu précise.

0

!1

!2

!3

f4
f6
!4 f8
f10
f12
f14
!5 f16

!6
0 100 200 300 400 500 600 700 800 900 1000

Figure 11: L2 errors between the polynomial solutions and the logarithmic solution.

F IGURE 3.6 – Erreur L2 entre la solution avec logarithme et la solution avec approximation
polynômiale f k de degré k − 1. 14

38
 Décomposition spinodale
Sur la figure 3.7, on présente l’état du mélange en 6 instants différents. Partant d’une
répartition aléatoire des deux composants, on assiste à une séparation de phase, nommée
décomposition spinodale. Par la suite, la longueur de l’interface va tendre à être minimisée,
jusqu’à devenir une ligne droite. Notons que les calculs doivent être précis (on a utilisé ici
des éléments finis quadrangulaires de degré jusqu’à 10), ce d’autant que le paramètre ε
est petit (ε mesure l’épaisseur de l’interface).

(a) (a)
t=0t=0 (b) (b)
t=0.01
t=0.01
4 Mar 2010
4 Mar 2010

(a) (a)
t=0t=0 (b) (b)
t=0.01
t=0.01

(a) t=0 (b) t=0.01

1 -2010
20101 - 2010

(a) (a)
t=0t=0 (b) (b)
t=0.01
t=0.01
4 Mar
Mar
1 -2010
version
version
- 4
1 - 4 Mar
version 1 - 4 Mar
1
versionversion
version
hal-00461231,
hal-00461231,hal-00461231,
hal-00461231,
hal-00461231,

(c) (c)
t=0.05
t=0.05 (d) (d)
t=0.2
t=0.2
hal-00461231,

(c) t=0.05 (d) t=0.2

(c) (c)
t=0.05
t=0.05 (d) (d)
t=0.2
t=0.2

(c) (c)
t=0.05
t=0.05 (d) (d)
t=0.2
t=0.2

(e) t=0.6 (f) t=1

(e) (e)
t=0.6
t=0.6 (f) (f)
t=1t=1
F IGURE 3.7 – Décomposition spinodale
Figure 3: Spinodal
Figure decomposition
3: Spinodal under
decomposition a logarithmic
under potential.
a logarithmic potential.
(e) (e)
t=0.6
t=0.6 (f) (f)
t=1t=1

Figure
Figure
3: Spinodal
3: Spinodal
decomposition
decomposition
under
under
a logarithmic
a logarithmic
potential.
potential.
(e) (e)
t=0.6
t=0.6 (f) (f)
t=1t=1
8398
Figure
Figure
3: Spinodal
3: Spinodal
decomposition
decomposition
under
under
a logarithmic
a logarithmic
potential.
potential.

8 8
3.5 Exploitation de résultats numériques
3.5.1 Réseau d’électrodes et matriçage cellulaire
Dans [PV10], on s’est intéressé à un microsystème permettant le positionnement de
cellules selon un réseau en vue d’un traitement individuel (transfert de gène, par exem-
ple), voir [HP91, MHG99, FFP+ 03]. La méthode s’appuie sur la force de diélectropho-
rèse que subit une cellule lorsqu’elle est soumise à un champ électrique non-uniforme.
Précisément, pour une cellule sphérique de rayon r dans un milieu de permittivité ε m , la
force de diélectrophorèse s’exprime en fonction du carré du champ appliqué E comme

F = 2πε m r3 Re [K ]∇( E2 ),

où le signe de la partie réelle de K (facteur de Clausius-Mosotti, dépendant des permittivités

de la cellule et du milieu) détermine le sens de la force. Dans le cas de la diélectrophorèse
négative (Re [K ] < 0), les cellules ont tendance à se positionner dans les minima de champ,
cf. F IG . 3.8.

F IGURE 3.8 – Un réseau d’électrodes rondes interdigitées.

Dès lors, il s’agit de déterminer la géométrie du réseau d’électrodes pour minimiser

le temps d’accès des cellules aux minima de champ. Du point de vue mathématique,
il est possible d’exprimer le temps de positionnement en fonction du carré du champ
électrique appliqué. On montre que les conditions d’optimalité peuvent se réduire à un
problème mono-dimensionnel et conduisent à l’équation intégrale d’Abel
; 1
u(t)
@ dt = cste,
0 | x − t|

d’où on peut extraire u, et donc le champ optimum à appliquer. On a ensuite mis en

œuvre une méthode de variation de frontière pour déterminer une géométrie d’élec-
trodes permettant d’approcher ce champ optimum (au sens des moindres carrés). Un
résultat avec une spline comportant 5 points de contrôle est visible sur la figure 3.9. Cette
approche pourtant relativement élémentaire permet un gain d’environ 20% sur le temps
de positionnement.

40
 •
•
• • • •

• • • •

F IGURE 3.9 – Électrodes rondes de référence (gauche) et électrodes optimisées (droite).

3.5.2 Partitions minimales et partitions nodales

Dans [BNHV10], on s’intéresse à un problème de partitionnement optimal. Soit Ω un
domaine borné de R2 fixé, et soit D = ( Di )i=1,...,k une partition de Ω en k sous-domaines.
Les valeurs propres de l’opérateur de Laplace-Dirichlet sur un domaine ω sont notées :

λ1 ( ω ) < λ2 ( ω ) ≤ λ3 ( ω ) ≤ · · ·

et Λ(D) est défini comme la plus grande première valeur propre sur la partition D :

Λ(D) = max λ1 ( Di ).
i =1,...,k

On considère alors le problème d’optimisation

Lk := min {Λ(D) ; D est une k − partition de Ω} . (3.1)

Ce problème a été largement étudié d’un point de vue théorique, voir [BBH98, CTV03,
HHOT09] par exemple. On peut montrer qu’il existe une partition minimale selon la dé-
finition précédente, ainsi que des propriétés qualitatives : les frontières de la partition
minimale sont des courbes régulières, sauf en un nombre fini de points, où elles sont des
unions de telles courbes se coupant à angles égaux.
Une des questions intéressantes concerne les liens entre partitions minimales et parti-
tions nodales (i.e. partitions fournies par les changements de signe des vecteurs propres
sur Ω). S’il est facile de vérifier que ces notions coïncident lorsque k = 2 :

L2 = λ 2 ( Ω ) ,

la question reste ouverte dès que k ≥ 3, même pour des géométries simples comme le
carré ou le disque. Dans [BNHV10], on a effectué des calculs numériques visant à étudier
cette question.
On décrit ici un des résultats obtenus : considérons un domaine Ω présentant un
axe de symétrie. On se restreint aux partitions symétriques (il n’est pas évident toutefois
que la partition minimale le soit...). En coupant selon l’axe de symétrie, on considère le
problème mixte Dirichlet-Neumann posé sur le demi-domaine Ω+ , voir F IG . 3.10 :


 −∆ϕ = λϕ dans Ω+ ,

∂n ϕ = 0 sur [ x0 , b], (3.2)



ϕ = 0 ailleurs.

41
 Ω+

a• • •b
x0

F IGURE 3.10 – Le problème mixte Dirichlet-Neumann.

Si on parvient à trouver x0 tel que la ligne nodale du deuxième vecteur propre pour
le problème (3.2) parte de x0 pour atteindre le bord de Ω, alors on aura un candidat pour
la 3-partition minimale de Ω après symétrisation.
Dans le cas du carré unité Ω = (0, 1)2 , les partitions données sur la figure 3.10 sont
proposées comme partitions minimales du carré unité – elles ont même valeur Λ(D).

F IGURE 3.11 – Partitions candidates pour L3 .

3.6 Perspectives
Une version C++ de la bibliothèque M ÉLINA sera bientôt rendue publique par D. Mar-
tin. Un travail de test sur des applications variées sera nécessaire, ainsi que des contribu-
tions dans le domaine du maillage : création et visualisation de maillages de haut degré,
projection entre maillages, etc. Ce travail de soutien d’un code de calcul pérenne est es-
sentiel, il permet la mise en œuvre de calculs numériques conséquents sans en imposer
la programmation complète.
Concernant les applications décrites, de nombreux prolongements sont envisagés ;
certains sont déjà en cours de développement. Citons, entre autres, le cas de l’équation
de Cahn-Hilliard avec bruit blanc, la validation d’une condition aux limites approchée
d’ordre 1 pour le calcul de profils en élasticité, l’étude de la supraconductivité en dimen-
sion 3, etc.

42
Bibliographie

[ABC+ 08] H. Ammari, E. Bonnetier, Y. Capdeboscq, M. Tanter, and M. Fink. Electri-

cal impedance tomography by elastic deformation. SIAM J. Appl. Math.,
68(6) :1557–1573, 2008.
[ACG06] Y. Amirat, G. A. Chechkin, and R. R. Gadyl8 shin. Asymptotics of simple
eigenvalues and eigenfunctions for the Laplace operator in a domain with
oscillating boundary. Zh. Vychisl. Mat. Mat. Fiz., 46(1) :102–115, 2006.
[AHH00] H. Ammari, L. Halpern, and K. Hamdache. Asymptotic behaviour of thin
ferromagnetic films. Asymptot. Anal., 24(3-4) :277–294, 2000.
[AKLZ10] H. Ammari, H. Kang, M. Lim, and H. Zribi. Layer potential techniques in
spectral analysis. part i : complete asymptotic expansions for eigenvalues
of the laplacian in domains with small inclusions. Trans. Amer. Math. Soc.,
362 :2901–2922, 2010.
[AKS06] H. Ammari, H. Kang, and F. Santosa. Scattering of electromagnetic waves
by thin dielectric planar structures. SIAM J. Math. Anal., 38(4) :1329–1342
(electronic), 2006.
[ALG98] H. Ammari and C. Latiri-Grouz. Approximate boundary conditions for
thin periodic coatings. In Mathematical and numerical aspects of wave propa-
gation (Golden, CO, 1998), pages 297–301. SIAM, Philadelphia, PA, 1998.
[All02] G. Allaire. Shape optimization by the homogenization method, volume 146 of
Applied Mathematical Sciences. Springer-Verlag, New York, 2002.
[AMM06] A. Ait Moussa and J. Messaho. Limit behavior of an oscillating thin layer.
In Proceedings of the 2005 Oujda International Conference on Nonlinear Analy-
sis, volume 14 of Electron. J. Differ. Equ. Conf., pages 21–33 (electronic), San
[Link], TX, 2006. Southwest Texas State Univ.
[APV98] Y. Achdou, O. Pironneau, and F. Valentin. Effective boundary condi-
tions for laminar flows over periodic rough boundaries. J. Comput. Phys.,
147(1) :187–218, 1998.
[BBH98] D. Bucur, G. Buttazzo, and A. Henrot. Existence results for some optimal
partition problems. Adv. Math. Sci. Appl., 8(2) :571–579, 1998.

43
[BDVV08] D. Brancherie, M. Dambrine, P. Villon, and G. Vial. Effect of surface de-
fects on structure failure : a two-scale approach. Eur. J. Comp. Mech., 17(5–
7) :613–624, 2008.
[BFM08] B. Bourdin, G. A. Francfort, and J.-J. Marigo. The variational approach to
fracture. J. Elasticity, 91(1-3) :5–148, 2008.
[BFNW08] D. Bucur, E. Feireisl, Š. Nečasová, and J. Wolf. On the asymptotic limit of
the Navier-Stokes system on domains with rough boundaries. J. Differen-
tial Equations, 244(11) :2890–2908, 2008.
[BHB06] M. F. Ben Hassen and E. Bonnetier. An asymptotic formula for the voltage
potential in a perturbed 6-periodic composite medium containing mispla-
ced inclusions of size 6. Proc. Roy. Soc. Edinburgh Sect. A, 136(4) :669–700,
2006.
[BL96] A. Bendali and K. Lemrabet. The effect of a thin coating on the scattering
of a time-harmonic wave for the Helmholtz equation. SIAM J. Appl. Math.,
56(6) :1664–1693, 1996.
[BL08] A. Bendali and K. Lemrabet. Asymptotic analysis of the scattering of a
time-harmonic electromagnetic wave by a perfectly conducting metal coa-
ted with a thin dielectric shell. Asymptot. Anal., 57(3-4) :199–227, 2008.
[BM10] D. Bresch and V. Milišić. High order multi-scale wall-laws, part i : The
periodic case. Quart. Appl. Math., posted on M.h 10, 2010 PII S0033-569X-
10-01135-0 (to appear in print), 2010.
[BMNP01] M. Bourlard, A. Maghnouji, S. Nicaise, and L. Paquet. Asymptotic ex-
pansion of the solution of a mixed Dirichlet-Ventcel problem with a small
parameter. Asymptot. Anal., 28(3-4) :241–278, 2001.
[BNBD+ 10] V. Bonnaillie-Noël, D. Brancherie, M. Dambrine, S. Tordeux, and G. Vial.
Effect of micro-defects on structure failure : coupling asymptotic analysis
and strong discontinuity approach. Eur. J. Comp. Mech., à paraître, 2010.
[BNDHV10a] V. Bonnaillie-Noël, M. Dambrine, F. Hérau, and G. Vial. Absorbing condi-
tions for the linear elasticity equations. en préparation, 2010.
[BNDHV10b] V. Bonnaillie-Noël, M. Dambrine, F. Hérau, and G. Vial. On generalized
Ventcel’s type boundary conditions for laplace operator in a bounded do-
main. SIMA, à paraître, 2010.
[BNDMV07] V. Bonnaillie-Noël, M. Dauge, D. Martin, and G. Vial. Computations of the
first eigenpairs for the Schrödinger operator with magnetic field. Comput.
Methods Appl. Mech. Engrg., 196(37-40) :3841–3858, 2007.
[BNDTV07] V. Bonnaillie-Noël, M. Dambrine, S. Tordeux, and G. Vial. On modera-
tely close inclusions for the Laplace equation. C. R. Math. Acad. Sci. Paris,
345(11) :609–614, 2007.
[BNDTV09] V. Bonnaillie-Noël, M. Dambrine, S. Tordeux, and G. Vial. Interactions
between moderately close inclusions for the Laplace equation. M3AS,
19(10) :1853–1882, 2009.

44
[BNHV10] V. Bonnaillie-Noël, B. Helffer, and G. Vial. Numerical simulations for nodal
domains and spectral minimal partitions. COCV, 16(1) :221–246, 2010.
[BR80] H. Blum and R. Rannacher. On the boundary value problem of the bihar-
monic operator on domains with angular corners. Math. Methods Appl. Sci.,
2(4) :556–581, 1980.
[BV00] E. Bonnetier and M. Vogelius. An elliptic regularity result for a composite
medium with “touching” fibers of circular cross-section. SIAM J. Math.
Anal., 31(3) :651–677 (electronic), 2000.
[CCDV06] G. Caloz, M. Costabel, M. Dauge, and G. Vial. Asymptotic expansion of
the solution of an interface problem in a polygonal domain with thin layer.
Asymptotic Analysis, 50(1, 2) :121–173, 2006.
[CD96] M. Costabel and M. Dauge. A singularly perturbed mixed boundary value
problem. Comm. Partial Differential Equations, 21 :1919–1949, 1996.
[CD99] D. Cioranescu and P. Donato. An introduction to homogenization, volume 17
of Oxford Lecture Series in Mathematics and its Applications. The Clarendon
Press Oxford University Press, New York, 1999.
[CD02] M. Costabel and M. Dauge. Weighted regularization of Maxwell equations
in polyhedral domains. A rehabilitation of nodal finite elements. Numer.
Math., 93(2) :239–277, 2002.
[CDMV03] M. Costabel, M. Dauge, D. Martin, and G. Vial. Weighted regularization of
Maxwell equations : computations in curvilinear polygons. In Numerical
mathematics and advanced applications, pages 273–280. Springer Italia, Milan,
2003.
[CFM09] A. Chambolle, G. A. Francfort, and J.-J. Marigo. When and how do cracks
propagate ? J. Mech. Phys. Solids, 57(9) :1614–1622, 2009.
[CH09] S. Chun and J. S. Hesthaven. High-order accurate thin layer approxima-
tions for time-domain electromagnetics. I. General metal backed coatings.
J. Comput. Appl. Math., 231(2) :598–611, 2009.
[CLP+ 07] E. Chahine, P. Laborde, J. Pommier, Y. Renard, and M. Salaün. Study of
some optimal XFEM type methods. In Advances in meshfree techniques, vo-
lume 5 of Comput. Methods Appl. Sci., pages 27–38. Springer, Dordrecht,
2007.
[Cos91] M. Costabel. A coercive bilinear form for Maxwell’s equations. J. Math.
Anal. Appl., 157(2) :527–541, 1991.
[CTV03] M. Conti, S. Terracini, and G. Verzini. An optimal partition problem related
to nonlinear eigenvalues. Journal of Funct. Anal., 198 :160–196, 2003.
[CV03] Y. Capdeboscq and M. Vogelius. A general representation formula for
boundary voltage perturbations caused by internal conductivity inhomo-
geneities of low volume fraction. M2AN, 37(1) :159–173, jan 2003.
[Dau88] M. Dauge. Elliptic Boundary Value Problems in Corner Domains – Smooth-
ness and Asymptotics of Solutions. Lecture Notes in Mathematics, Vol. 1341.
Springer-Verlag, Berlin, 1988.

45
[DBDI05] N. Dominguez, D. Brancherie, L. Davenne, and A. Ibrahimbegović. Pre-
diction of crack pattern distribution in reinforced concrete by coupling a
strong discontinuity model of concrete cracking and a bond slip of reinfor-
cement model. Engin. Comput., 22(5–6) :558–581, 2005.
[DM93] G. Dal Maso. An introduction to Γ-convergence. Progress in Nonlinear Diffe-
rential Equations and their Applications, 8. Birkhäuser Boston Inc., Boston,
MA, 1993.
[DTV09] M. Dauge, S. Tordeux, and G. Vial. Self-similar perturbation near a corner :
matching versus multiscale expansions for a model problem. To appear,
2009.
[DV05] M. Dambrine and G. Vial. On the influence of a boundary perforation on
the dirichlet energy. Control and Cybernetics, 34(1) :117–136, 2005.
[DV07] M. Dambrine and G. Vial. A multiscale correction method for local sin-
gular perturbations of the boundary. M2AN Math. Model. Numer. Anal.,
41(1) :111–127, 2007.
[EM77] B. Engquist and A. Majda. Absorbing boundary conditions for the nume-
rical simulation of waves. Math. Comp., 31(139) :629–651, 1977.
[EN93] B. Engquist and J. C. Nedelec. Effective boundary conditions for electro-
magnetic scattering in thin layers. Rapport interne 278, CMAP, 1993.
[FFP+ 03] M. Frénéa, S. Faure, B. Pioufle, P. Coquet, and H. Fujita. Positioning li-
ving cells on a high-density electrode array by negative dielectrophoresis.
Materials Science and Engineering : C, 23(5) :597–603, 2003.
[Giv91] D. Givoli. Nonreflecting boundary conditions. J. Comput. Phys., 94(1) :1–29,
1991.
[GKL99] G. Geymonat, F. Krasucki, and S. Lenci. Mathematical analysis of a bonded
joint with a soft thin adhesive. Math. Mech. Solids, 4(2) :201–225, 1999.
[GMV10] L. Goudenège, D. Martin, and G. Vial. High order finite element calcula-
tions for the deterministic Cahn-Hilliard equation. en révision, 2010.
[Gri85] P. Grisvard. Boundary value problems in non-smooth domains. Pitman, Lon-
don, 1985.
[HHOT09] B. Helffer, T. Hoffmann-Ostenhof, and S. Terracini. Nodal domains and
spectral minimal partitions. Ann. Inst. H. Poincaré Anal. Non Linéaire,
26 :101–138, 2009.
[HJ01] H. Haddar and P. Joly. Effective boundary conditions for thin ferroma-
gnetic coatings. Asymptotic analysis of the 1D model. Asymptot. Anal.,
27(2) :127–160, 2001.
[HL96] C. Hazard and M. Lenoir. On the solution of time-harmonic scattering
problems for Maxwell’s equations. SIAM J. Math. Anal., 27(6) :1597–1630,
1996.
[HP91] Y. Huang and R. Pethig. Electrode design for negative dielectrophoresis.
Measurement Science and Technology, 2(12) :1142–1146, 1991.

46
[HR95] L. Halpern and J. Rauch. Absorbing boundary conditions for diffusion
equations. Numer. Math., 71(2) :185–224, 1995.
[Il’92] A. M. Il’lin. Matching of asymptotic expansions of solutions of boundary
value problems. Translations of Mathematical Monographs, 1992.
[JM01] W. Jäger and A. Mikelić. On the roughness-induced effective boundary
conditions for an incompressible viscous flow. J. Differential Equations,
170(1) :96–122, 2001.
[JT06] P. Joly and S. Tordeux. Matching of asymptotic expansions for wave pro-
pagation in media with thin slots i : The asymptotic expansion. Multiscale
Modeling and Simulation : A SIAM Interdisciplinary Journal, 5(1) :304–336,
2006.
[KMO04] F. Krasucki, A. Münch, and Y. Ousset. Mathematical analysis of nonlinear
bonded joint models. Math. Models Methods Appl. Sci., 14(4) :535–556, 2004.
[KMZ08] D. Knees, A. Mielke, and C. Zanini. On the inviscid limit of a model for
crack propagation. Math. Models Methods Appl. Sci., 18(9) :1529–1569, 2008.
[Kon67] V. A. Kondrat8 ev. Boundary value problems for elliptic equations in do-
mains with conical or angular points. Trans. Moscow Math. Soc., 16 :227–
313, 1967.
[LD91] H. Le Dret. Problèmes variationnels dans les multi-domaines, volume 19 of
Recherches en Mathématiques Appliquées [Research in Applied Mathematics].
Masson, Paris, 1991. Modélisation des jonctions et applications. [Modeling
of junctions and applications].
[LS00] T. Lewiński and J. Sokołowski. Topological derivative for nucleation of
non-circular voids. The Neumann problem. In Differential geometric me-
thods in the control of partial differential equations (Boulder, CO, 1999), volume
268 of Contemp. Math., pages 341–361. Amer. Math. Soc., Providence, RI,
2000.
[LT92] K. Lemrabet and D. Teniou. Vibrations d’une plaque mince avec raidisseur
sur le bord. Maghreb Math. Rev., 2(1) :27–41, 1992.
[LVLH92] M. Lenoir, M. Vullierme-Ledard, and C. Hazard. Variational formulations
for the determination of resonant states in scattering problems. SIAM J.
Math. Anal., 23(3) :579–608, 1992.
[Mar10] D. Martin. M ÉLINA, bibliothèque de calculs éléments finis.
[Link] 1986–2010.
[Mel02] J. M. Melenk. hp-finite element methods for singular perturbations, volume
1796 of Lecture Notes in Mathematics. Springer-Verlag, Berlin, 2002.
[MHG99] H. Morgan, M. Hughes, and N. Green. Separation of submicron biopar-
ticles by dielectrophoresis. Biophysical journal, 77(1) :516–525, 1999.
[Mil09] V. Milišić. Very weak estimates for a rough Poisson-Dirichlet problem with
natural vertical boundary conditions. Methods Appl. Anal., 16(2) :157–185,
2009.

47
[MNP00] Vladimir G. Maz’ya, Sergey A. Nazarov, and B. A. Plamenevskij. Asymp-
totic theory of elliptic boundary value problems in singularly perturbed domains.
Birkhäuser, Berlin, 2000.
[MS03] M. Mabrouk and H. Samadi. Linear and semi-linear reinforcement pro-
blems by thin layers. Z. Angew. Math. Phys., 54(2) :349–375, 2003.
[MV07] A. Madureira and F. Valentin. Asymptotics of the Poisson problem in do-
mains with curved rough boundaries. SIAM J. Math. Anal., 38(5) :1450–
1473 (electronic), 2006/07.
[Néd80] J.-C. Nédélec. Mixed finite elements in R3 . Numer. Math., 35(3) :315–341,
1980.
[Néd01] J.-C. Nédélec. Acoustic and electromagnetic equations, volume 144 of Applied
Mathematical Sciences. Springer-Verlag, New York, 2001. Integral represen-
tations for harmonic problems.
[Nic93] S. Nicaise. Polygonal interface problems, volume 39 of Methoden und Verfahren
der Mathematischen Physik [Methods and Procedures in Mathematical Physics].
Verlag Peter D. Lang, Frankfurt am Main, 1993.
[NO93] S. A. Nazarov and M. V. Olyushin. Perturbation of the eigenvalues of the
Neumann problem due to the variation of the domain boundary. Algebra i
Analiz, 5(2) :169–188, 1993.
[NP94] S. A. Nazarov and B. A. Plamenevsky. Elliptic problems in domains with pie-
cewise smooth boundaries, volume 13 of de Gruyter Expositions in Mathematics.
Walter de Gruyter & Co., Berlin, 1994.
[Oli95] J. Olivier. Continuum modelling of strong discontinuities in solid mecha-
nics using damage models. Comput. Mech., 17(1–2) :49–61, 1995.
[Poi08] C. Poignard. Asymptotics for steady-state voltage potentials in a bidimen-
sional highly contrasted medium with thin layer. Math. Methods Appl. Sci.,
31(4) :443–479, 2008.
[PV10] M. Pierre and G. Vial. Best design for a fastest cells selecting process.
DCDS-S, à paraître, 2010.
[Rah08] L. Rahmani. The effect of a thin layer on a nonlinear thermoelastic plate.
Appl. Math. Sci. (Ruse), 2(49-52) :2489–2499, 2008.
[RV08] L. Rahmani and G. Vial. Reinforcement of a thin plate by a thin layer. Math.
Methods Appl. Sci., 31(3) :315–338, 2008.
[RV10] L. Rahmani and G. Vial. Multi-scale asymptotic expansion for a singular
problem of a free plate with thin stiffener. en préparation, 2010.
[SOA93] J.C. Simo, J. Oliver, and F. Armero. An analysis of strong discontinuity in-
duced by strain softening solutions in rate-independent solids. J. Comput.
Mech., 12 :277–296, 1993.
[SP74] E. Sánchez-Palencia. Problèmes de perturbations liés aux phénomènes de
conduction à travers des couches minces de grande résistivité. J. Math.
Pures Appl. (9), 53 :251–269, 1974.

48
[SŻ99] J. Sokołowski and A. Żochowski. On the topological derivative in shape
optimization. SIAM J. Control Optim., 37(4) :1251–1272 (electronic), 1999.
[Tor04] S. Tordeux. Méthodes asymptotiques pour la propagation des ondes dans
les milieux comportant des fentes. Thèse de doctorat., 2004.
[Van75] M. VanDyke. Perturbation methods in fluid mechanics. The Parabolic Press.,
1975.
[Via05] G. Vial. Efficiency of approximate boundary conditions for corner domains
coated with thin layers. CR. Acad. Sci., 340 :215–220, January 2005.
[VV00] M. S. Vogelius and D. Volkov. Asymptotic formulas for perturbations in
the electromagnetic fields due to the presence of inhomogeneities of small
diameter. M2AN Math. Model. Numer. Anal., 34(4) :723–748, 2000.
[WI90] E.L. Wilson and A. Ibrahimbegović. Use of incompatible displacement
modes for the calculation of element stiffnesses or stresses. Finite Elem.
Anal. Des., 7 :229–241, 1990.

49
École normale supérieure de Cachan - Antenne de Bretagne
Campus de Ker Lann - Avenue Robert Schuman - 35170 BRUZ

Tél : +33(0)2 99 05 93 00 - Fax : +33(0)2 99 05 93 29