Application de la méthode des éléments finis à l’électromagnétisme

Raphaël du Roy de Chaumaray

Florian Lavrilleux
23 août 2023

1 Introduction
La méthode des éléments finis est une méthode de résolution numérique d’équations aux dérivées par-
tielles. Grâce à cette méthode, on peut décrire approximativement le comportement de n’importe quel système
physique, peu importe sa complexité, tant qu’il est continu et qu’il obéit à une EDP linéaire. On souhaite
résoudre l’équation de Helmholtz et ainsi donner l’expression du champ électrique associé à une onde plane qui
se propage dans le vide et qui rencontre un objet diffractant quelconque. On se place en régime harmonique
afin de se libérer de la dépendance temporelle du champ électrique. On peut ainsi tracer des cartes de champ
en fonction de la position uniquement.

2 Problème de diffraction 1D
2.1 Cadre mathématique du problème
On souhaite modéliser le champ électrique associé à une onde plane incidente rencontrant un objet
diffractant d’épaisseur e. Le problème étant symétrique par translation suivant ŷ et ẑ, on se ramène au cas d’un
problème 1D. Le champ électrique incident a pour expression : ui (x) = exp(ik0 x). Notre but est de donner la
solution ut (x) sur le segment [0, L], où ut est le champ dit "total". On considère dans ce problème un milieu
amagnétique, c’est-à-dire de perméabilité magnétique relative µr = 1. Les champs ui et ut vérifient tous deux
l’équation de propagation de Helmholtz 1D :

u′′ (x) + k 2 u(x) = 0 (1)

On note εr (x) la fonction

p qui, à tout x, associe la permittivité diélectrique du milieu. On associe alors à k
la fonction k(x) = k0 εr (x), où k0 = 2π λ0 , avec λ0 la longueur d’onde de la lumière dans le vide. Dans notre
problème, une lame de permittivité εmr est placée dans le vide, entre les points l et l + e. On a donc :
(
ε0r , ∀x ∈ [0, l[∪]l + e, L]
εr (x) =
εm
r , ∀x ∈ [l, l + e]

En soustrayant l’équation (1) appliquée au champ ut à celle appliquée au champ ui , on obtient :

   
u′′t (x) − u′′i (x) + k02 εr (x)ut (x) − ε0r ui (x) = 0 (2)

On définit le champ diffracté comme ud := ut − ui . En transformant l’équation (2), on obtient l’équation de

propagation vérifiée par le champ diffracté :

u′′d (x) + k02 εr (x)ud (x) = k02 εm

 
r − εr (x) ui (x) (3)

Le terme de droite est équivalent à un terme source. On le notera S(x). Tout se passe comme si une source
(équivalente à un courant ⃗j) était localisée dans la lame diffractante. Le champ diffracté satisfait donc la
condition d’onde sortante de Sommerfeld, qui, en 1D, est la suivante :
 
lim u′d ∓ ik0 ud = 0
x→±∞

1
Cela signifie que loin de la source (donc à une distance D ≫ e de la lame), le champ ud vérifie la relation
u′d = ±ik0 ud (selon qu’on soit à droite ou à gauche de la lame). Cela nous sera utile par la suite. Le but est
de résoudre numériquement l’équation (3), qui est la formulation dite faible du problème. Cette formulation est
équivalente à la formulation dite forte (on va l’admettre), qui est la suivante :
Z L Z L Z L
2 ′′ 2
∀w ∈ L (R), ud (x)w(x)dx + k0 εr (x)ud (x)w(x)dx = S(x)w(x)dx
0 0 0

En intégrant par parties le premier terme, cela nous donne :

h iL Z L Z L Z L
u′d (x)w(x) − u′d (x)w′ (x)dx + k02 εr (x)ud (x)w(x)dx = S(x)w(x)dx (4)
0 0 0 0

Pour la programmation, un problème intervient : on a dim(L2 ) = +∞, donc on ne peut pas tester l’égalité
sur l’ensemble des fonctions-test de L2 . Pour remédier à ce problème, on va diviser le segment [0, L] en N − 1
sous-ensembles Ii . Il y a donc N noeuds xi , i = 0, ..., N . On va introduire l’espace des fonctions-triangle, définies
comme suit : x−x
i−1 x − xi−1
 = si x ∈ Ii
 xi − xi−1 hi


φi (x) = xi+1 − x = xi+1 − x si x ∈ Ii+1
 x − xi hi+1
 i+1


0 sinon
On va aussi faire en sorte que les points x = l et x = l + e soient des noeuds. Pour cela, soient N1 , N2 et N3 le
nombre de noeuds dans les ensembles [0, l], [l, l + e] et [l + e, L], on a :
     
L−e p 
0
e p 
m
N1 = N3 = Re × n εr N2 = Re × n εr
2λ0 λ0
n est un entier arbitraire. Plus n est grand, plus l’approximation de la fonction ud sera précise. En effet, le but
du programme est d’approximer ud par une fonction continue par morceaux, soit une combinaison linéaire de
fonctions triangles. En supposant qu’elle soit dérivable terme à terme, et en injectant la nouvelle expression de
ud dans (4), on obtient :
n
! n
!
h iL Z L X Z L X Z L
′ ′ ′ 2
∀j ∈ [0, N ], ud φj − αi φi φj dx + k0 εr (x) αi φi φj dx = S(x)φj dx (5)
0 0 0 0
i=0 i=0

Si on considère que l’épaisseur de la lame est "négligeable" devant la longueur de l’intervalle [0, L], on peut
appliquer la condition de Sommerfeld au terme entre crochets, ce qui donne :
h iL    
u′d φj = ik0 αN φN (L)φj (L) + α0 φ0 (0)φj (0) = ik0 αN φj (L) + α0 φj (0) , ∀j
0

 qui est non nulle en L est φN et la seule

La seule fonction triangle


qui est non nulle est 0 est φ0 , donc le terme
entre crochets vaut ik0 αN φN (L) + α0 φ0 (0) = ik0 αN + α0 . Cela peut encore se réécrire comme :
   
ik0 0 ... 0 α0
 ..   .. 
0 0 ... .    . 
 
..  ·  ..  = BV

 . .. ..
 .. . . .   . 
0 ... . . . ik0 αN
L’équation (5) se ramène donc à une équation matricielle d’inconnue V qui est la suivante :
(B − S + k02 M )V = C
En effet, comme ud est approximée par une fonction continue par morceaux, la seule connaissance des αi
est nécessaire pour déterminer numériquement ud . La matrice S est appelée matrice de raideur ("stiffness"
Z L
en anglais), de coefficients Sij = φ′i φ′j dx. La matrice M est appelée matrice de masse, de coefficients
Z L 0 Z L
Mij = εr (x)φi φj dx. Enfin, le vecteur C a pour coefficients Cj = S(x)φj dx.
0 0

2
2.2 Résolution numérique
2.2.1 Calcul numérique d’intégrales par la méthode de Simpson
On considère une fonction f dont on veut évaluer l’intégrale entre deux points a et b. Prenons n + 1
points équidistants x0 = a < x1 < . . . < xn = b et considérons un polynôme interpolant f en ces points. Le but
est d’approximer l’intégrale de f par une intégrale sur ce polynôme entre a et b (car les polynômes sont faciles
à intégrer). Supposons que l’on peut définir n + 1 polynômes Pi d’ordre n, tels que :

Pi (xj ) = δij

Le polynôme P approximant f peut s’exprimer sur la base des Pi :

n
X
P (x) = f (xi )Pi (x)
i=0

P est un polynôme d’ordre au plus n puisque c’est une combinaison linéaire de polynômes d’ordre n et on a bien
P (xj ) = f (xj )Pj (xj ) = f (xj ), ce qui signifie que P passe par tous les points d’interpolation. On peut montrer
que les Pi ont pour expression :
n
Y (x − xj )
Pi (x) =
j=0
(x i − xj )
j̸=i

On a donc :
Z b Z b Z n
bX
I= f (x)dx ≈ P (x)dx = f (xi )Pi (x)dx
a a a i=0
n Z b
X
= f (xi )Pi (x)dx
i=0 a

Xn Z b
= f (xi ) Pi (x)dx
i=0 |a {z }
ωi

Les coefficients ωi se calculent en effectuant le changement de variable x −→ y = (x − a)/h, où h = (b − a)/n

est la longueur d’un intervalle [xi , xi+1 ] :
n
bY n
nY n
nY
(x − xj ) b−a (yh + a − xj ) b−a
Z Z Z
A
ωi = dx = dy = dy
a j=0
(xi − xj ) n 0 j=0
(xi − xj ) n 0 j=0
B
j̸=i j̸=i j̸=i

On a A = (yh + a − a)(yh + a − (a + h)) . . . (yh + a − b) = (yh + a − a)(yh + a − (a + h)) . . . (yh + a − (a + hn)) =

n
Y
hn+1 (y − k) et B = (xi − a)(xi − a − h) . . . (xi − a − nh) = (a + ih − a)(a + ih − a − h) . . . (a + ih − a − nh) =
k=0
k̸=i
Yn
hn+1 (i − k).
k=0
k̸=i
On voit que le terme B est négatif si n − i est impair et positif si n − i est pair. On a donc B = (−1)n−i i!(n − i)!.
Le coefficient ωi a donc pour expression :
n
n Y
b − a (−1)n−i
Z
ωi = (y − k)dy
n i!(n − i)! 0 k=0
k̸=i

La méthode de Simpson est l’application pour n = 2. L’intégrale s’exprime comme suit :

     
b−a 4 b − a −8 a+b b−a 4 b−a a+b
I= × f (a) − × f + × f (b) = f (a) + 4f + f (b)
4 6 2 6 2 4 6 6 2

3
2.2.2 Application à la matrice de masse
La matrice M fait intervenir le produit φi φj . Le support de ces deux fonctions est commun pour |i−j| ≤ 1.
Cela veut dire que la matrice M est tridiagonale. On a :
Z L Z xi Z xi+1
Mii = εr (x)φ2i (x)dx = εr (x)φ2i (x)dx + εr (x)φ2i (x)dx
0 xi−1 xi

Comme εr (x) est constante sur tous les intervalles Ii , on a :

Z xi Z xi+1
2
Mii = εr (Ii−1 ) φi (x)dx + εr (Ii ) φ2i (x)dx
xi−1 xi
hi 
2  hi+1  1 2

= εr (Ii−1 ) × 0 + 4 21 + 1 + εr (Ii )


× 1+4 2 +0
6 6
hi hi+1
= εr (Ii−1 ) + εr (Ii ) , i = 1, ..., N − 1
3 3
h1 hN
Pour i = 0 et i = N , on a M00 = et MN N = (puisque les fonctions triangles φ0 et φN sont "tronquées"
3 3
et que ε vaut 1 aux bords du domaine). On voit que la matrice est symétrique. En appliquant un raisonnement
similaire, on obtient :
Z L Z xi+1
hi+1
Mi,i+1 = Mi+1,i = εr (x)φi (x)φi+1 (x)dx = εr (x)φi (x)φi+1 (x)dx = εr (Ii )
0 xi 6

La matrice M prend donc la forme suivante :

 
h1 εr (I0 )h1
 3 6 
 εr (I0 )h1 εr (I0 )h1 +εr (I1 )h2 εr (I1 )h2 
 6 3 6 
 
 εr (I1 )h2 εr (I1 )h2 +εr (I2 )h3 εr (I2 )h3 
 6 3 6 
M =
 .. .. ..


 . . . 
 

 εr (IN −2 )hN −1 εr (IN −2 )hN −1 +εr (IN −1 )hN εr (IN −1 )hN 

 6 3 6 
εr (IN −1 )hN hN
6 3
     
2 1
1 2   2 1   
     
εr (I0 )h1   εr (I1 )h2  1 2  + . . . + εr (IN −1 )hN
  
=  +   
6 


 6 


 6 



     2 1
1 2

2.2.3 Application à la matrice de raideur

On suit le même raisonnement. La dérivée de la fonction triangle a pour expression :

1 si x ∈ Ii

′
φi (x) = −1 si x ∈ Ii+1

0 sinon


On en déduit donc que :

Z xi Z xi+1
′2 hi hi+1 1 1
Sii = φi (x)dx + φ′2
i (x)dx = + ; S00 = ; SN N =
xi−1 xi 6 6 h1 hN

1
Si,i+1 = Si+1,i = −
hi+1

4
La matrice S prend donc la forme suivante :
 
1
− h11
 h11
   
− h 1 1
− h12
 1 −1
 1 h1 + h2

 −1 1   
 − h12 1
h2 + h3
1     
= 1  1
     
S= .. .. ..  h1 
 + ... +  

 . . .  

 hN 




1 1 1

1 
   1 −1

 − hN −1 hN −1 + hN − hN  −1 1
− h1N 1
hN

2.2.4 Application au vecteur "source" C

Ici, on va utiliser la méthode des trapèzes. C’est la méthode de résolution numérique des intégrales pour
n = 1. Cela nous donne :
−1 1 b − a 
I = ω0 f (a) + ω1 f (b) = −(b − a) × f (a) + (b − a) × f (b) = f (a) + f (b)
2 2 2
La composante i du vecteur C s’exprime donc comme :
Z xi+1
Ci = S(x)φi dx
xi−1
Z xi Z xi+1
= S(x)φi dx + S(x)φi dx
xi−1 xi
hi   h 
i+1

= S(xi−1 ) φi (xi−1 ) +S(xi ) φi (xi ) + S(xi ) φi (xi ) +S(xi+1 ) φi (xi+1 )
2 | {z } | {z } 2 | {z } | {z }
=0 =1 =1 =0
S(xi )
= [hi + hi+1 ]
2
Cela implique que :
 
S(x0 )h1 /2      
 S(x1 )(h1 + h2 )/2
 S(x0 )
 h 
 + h2 S(x1 ) + . . . + hN
1 S(x1 )

C= ..    
=
 
.
 
  2   2 S(x )
2  2 S(xN −1 )
S(xN −1 )(hN −1 + hN )/2 S(xN )
S(xN )hN /2

2.3 Convergence
 2
 At 
Nous avons vérifié la convergence de notre programme en calculant le coefficient de réflexion R =   ,
Ai
avec Ai et At respectivement les amplitudes du champ incident et du champ transmis. Nous l’avons comparé à
√ n2 − n1
sa valeur théorique grâce au lien entre indices de réfraction et permittivités n = ε. On a donc R = =
√ √ n2 + n1
√ε2 − √ε1 en incidence normale. Cela nous permet de tracer le graphe 1. On voit que la convergence est rapide
ε2 + ε1
avec un écart relatif inférieur à 5% dès n = 80.

2.4 Résultats
Dans le cas 1D, on obtient le résultat présenté figure 2. On observe un champ sinusoïdal atténué, qui est
conforme à nos attentes.

5
(a) Indice de réfraction en fonction du facteur de précision (b) Écart relatif à la valeur théorique en fonction du
facteur de précision

Figure 1 – Validation du modèle à l’aide du calcul du coefficient de réflexion

Figure 2 – Champ total pour les paramètres suivants :

λ0 = 700 nm; e1 = e3 = λ0 ; e2 = 3λ0 ; ε = 9 + 2i
Les barres noires délimitent l’obstacle.

6
3 Problème de diffraction 2D
3.1 Cadre mathématique du problème
3.1.1 Génération du maillage
On considère un domaine borné Ω ⊂ R2 . Un maillage M est une partition de Ω en éléments polygonaux
K qui vérifie les propriétés suivantes :
[
1. Ω = K : les K recouvrent Ω,
K∈M
◦ ◦
2. ∀(Ki , Kj ) ∈ {K} × {K}, i ̸= j, Ki ∩ Kj = ∅ : l’intersection de l’intérieur de deux éléments distincts est
vide.
Le maillage est dit conforme lorsque l’intersection de deux éléments distincts est soit vide soit un élément
commun aux deux (sommet ou arête). Soient K1 , K2 , . . . , Kn les n arêtes d’un polygone K ∈ M, de longueurs
h1 , . . . , hn , on définit localement la taille caractéristique du maillage comme :

hK = max hi
i=1,...,n

La longueur caractéristique globale du maillage se définit alors naturellement comme suit :

h = max hK
K∈M

Par la suite, les éléments K du maillage seront des triangles quelconques. Nous allons maintenant voir comment
est généré un maillage triangulaire. On a pour cela besoin de deux données : les nombres np de points et
nt d’éléments. Les points sont numérotés de 0 à np − 1 et les triangles de 0 à nt − 1. Comme nous sommes
en 2D, l’ordre de la numérotation est arbitraire. On se place dans la base canonique B = (ex , ey ) de R2 .
Un point Pi , i ∈ {0, . . . , np − 1} est repéré par les coordonnées xi et yi dans cette base. Soit un triangle
(0) (1) (2)
Ki , i ∈ {0, . . . , nt − 1}, on note ses sommets Pi , Pi et Pi . On créé alors les matrices P et T , nommées
respectivement matrices de points et de connectivité, et définies comme suit :
 
(0) (0) (0)
P P . . . P nt −1 
   0 1
x x1 . . . xnp −1
P = 0 et T = P0(1) P1(1) . . . Pn(1)
 
y0 y1 . . . ynp −1 t −1 


(2) (2) (2)
P0 P1 . . . Pnt −1

Ces matrices vont être nécessaires dans notre programme, et nous permettront d’alléger les notations. Le domaine
total sera noté Ω. On choisit un cercle de rayon L. L’objet diffractant sera associé au domaine noté Ωd , et on
a Ωd ⊂ Ω. On n’impose aucune condition sur la nature du bord ∂Ωd . Il peut être autant lisse que polygonal.
Le milieu ambiant est donc Ω\Ωd . On a besoin de séparer ces deux milieux pour deux raisons. D’une part, on
aura en général h(Ωd ) ̸= h(Ω\Ωd ). D’autre part, on a besoin que le maillage soit "cohérent", c’est à dire qu’il
recouvre Ωd et Ω\Ωd indépendamment pour éviter les superpositions de maillages.

3.1.2 Mise en forme des équations

Cette fois-ci, le champ incident a pour expression :

ui (r) = exp(ik · r) (6)

Dans la base B, on identifie les vecteurs position r = [x, y] et d’onde k = [k0 cos θ, k0 sin θ], où θ est l’angle
d’incidence. On se place dans le cas d’une polarisation TE (transverse électrique), ce qui signifie que l’on a
µr = 1 dans tout l’espace, où µr est la perméabilité magnétique relative. Soit u(r) le champ total. Il satisfait
l’équation de Helmholtz, qui est la suivante :

∆u(r) + k02 εr (r)u(r) = 0

En posant, comme dans la première partie, ud = u − ui , on aboutit à l’équation :

∆ud (r) + k02 εr (r)ud (r) = ∆ui (r) + k02 εm

 
r − εr (r) ui (r)

7
  
On définit le terme source S(r) = ∆ui (r) + k02 εm r − εr (r) ui (r). On est donc dans la même situation que le
cas 1D. L’onde rayonne depuis l’objet diffractant, et satisfait la condition de rayonnement de Sommerfeld, qui
en 2D s’exprime comme suit :
√
 
∂
lim r − ik0 ud (r) = 0 (7)
r→∞ ∂r
p
r = x2 + y 2 est la distance à l’interface de l’objet diffractant. On note dΩ = dxdy l’élément infinitésimal
d’aire. La formulation forte du problème est alors la suivante :
Z Z Z
∀w ∈ L2 (R2 ), ∆ud w dΩ + k02 εr (r)ud w dΩ = Sw dΩ (8)
Ω Ω Ω
 
Le laplacien de ud peut s’écrire ∆ud = ∇ · ψ, où ψ = ∇ud . On définit aussi n(x, y) = nx (x, y), ny (x, y) le
vecteur unitaire normal au bord ∂Ω, dirigé vers l’extérieur. Le premier terme devient alors :
Z Z
∆ud w dΩ = (∇ · ψ) w dΩ
Ω
ZΩ Z
∂ψ ∂ψ
= w dΩ + wd Ω
∂x Ω ∂y
ZΩ Z Z Z
∂w ∂w
= ψwnx dΩ − ψ dΩ + ψwny dΩ − ψ dΩ
∂Ω Ω ∂x ∂Ω Ω ∂y
Z Z
= (n · ψ)w dΩ − ψ · (∇w) dΩ
Z∂Ω  
Ω
Z
= n · (∇ud ) w dΩ − (∇ud ) · (∇w) dΩ
∂Ω Ω

Maintenant, comme dans la première partie, on introduit un espace discret de fonctions {φi } (que l’on appelera
"fonctions chapiteau") semblable à L2 (R2 ), de dimension np , le nombre de points du maillage. Ces fonctions
satisfont la relation suivante :
∀i, j ∈ {0, nj − 1}, φi (Pj ) = δij
np −1
X
On pose que ud = αi φi . L’équation (8) devient alors :
i=0

np −1 np −1
Z Z ! Z ! Z
  X X
∀j ∈ {0, np −1}, n·(∇ud ) φj dΩ− ∇ α i φi ·∇φj dΩ+k02 εr (r) αi φi φj dΩ = Sφj dΩ
∂Ω Ω i=0 Ω i=0 Ω
(9)
Transformons le premier terme. En utilisant la règle de dérivation à la chaine, on obtient :
∂ud ∂ud ∂ud ∂r ∂ud ∂r
n · ∇ud = nx + ny = nx + ny
∂x ∂y ∂r ∂x ∂r ∂y
p
Comme r = x2 + y 2 , on obtient :
( )
∂ud x y ∂ud  n · r 
n · ∇ud = nx p + ny p =
∂r x2 + y 2 x2 + y 2 ∂r r

On a défini le domaine Ω comme un cercle de rayon L. Si on suppose que L ≫ λ, où λ est la taillecaractéristique

n · r
de l’objet diffractant, on peut alors appliquer la condition (7) afin d’obtenir n · ∇ud = ik0 ud . Enfin,
r
comme on étudie la fonction sur le bord du domaine, on a n · r = ∥n∥ ∥r∥ = r (puisqu’on a posé que n était de
norme unité). Finalement, le premier terme de l’équation (9) devient :
np −1
Z  Z Z !
 X
n · ∇ud φj dΩ = ik0 ud φj dΩ = ik0 αi φi φj dΩ
∂Ω ∂Ω ∂Ω i=0
 T
On obtient alors une équation matricielle d’inconnue V = α0 ... αnp −1 identique à celle du cas 1D :

(B − S + k02 M ) · V = C

8
3.2 Résolution numérique
3.2.1 Calcul numérique d’intégrales doubles sur un triangle
On considère une fonction f (x, y) et un triangle quelconque K formé par les points Pi (xi , yi ), i = 1, 2, 3.
Calculons l’aire de ce triangle, notée AK . Pour cela, on considère un quatrième point P4 qui est tel que le
quadrilatère P1 P2 P4 P3 soit un parallélogramme. L’aire de ce parallélogramme, notée A, a pour expression :
 
x2 − x1 x3 − x1
A = ∥P1 P2 ∧ P1 P3 ∥ = |(x2 − x1 )(y3 − y1 ) − (y2 − y1 )(x3 − x1 )| = det
y2 − y1 y3 − y1
On obtient l’aire du triangle en utilisant le fait que A = 2AK .
On considère le triangle K ′ formé par les points p1 (0, 0), p2 (1, 0) et p3 (0, 1). On définit un nouveau système de
coordonnées (ξ, η) tel que pi (ξ, η) = Pi (x, y). Pour cela on construit deux applications linéaires, dépendantes du
(K) (K)
triangle K, de R2 dans R : ϕ1 et ϕ2 , qui sont telles que :
(
(K) (K) (K)
ϕ1 (0, 0) = x1 , ϕ1 (1, 0) = x2 , ϕ1 (0, 1) = x3
(K) (K) (K)
ϕ2 (0, 0) = y1 , ϕ2 (1, 0) = y2 , ϕ2 (0, 1) = y3
On montre que : (
(K)
ϕ1 (ξ, η) = x1 + (x2 − x1 )ξ + (x3 − x1 )η
(K)
ϕ2 (ξ, η) = y1 + (y2 − y1 )ξ + (y3 − y1 )η
Le jacobien du changement de coordonnées (ξ, η) −→ (x, y) s’exprime donc comme suit :
 
(K) (K)
∂ϕ1 ∂ϕ1
  
x2 − x1 x3 − x1

J =  ∂ξ ∂η  =

 ∂ϕ(K) ∂ϕ(K)  y2 − y1 y3 − y1

2 2
∂ξ ∂η
On a alors : |det J| = A = 2AK
On souhaite maintenant calculer l’intégrale de f sur le triangle K. On utilise la formule du changement de
variables pour obtenir :
Z Z   Z  
(K) (K) (K) (K)
I= f (x, y)dxdy = f ϕ1 (ξ, η), ϕ2 (ξ, η) |det J|dξdη = 2AK f ϕ1 (ξ, η), ϕ2 (ξ, η) dξdη
K K′ K′
′
Avec le triangle K tel qu’on l’a défini, on peut définir les bornes d’intégration. On a :
Z 1 Z 1−ξ  
(K) (K)
I = 2AK f ϕ1 (ξ, η), ϕ2 (ξ, η) dξdη (10)
0 0
Pour la suite, on utilise les méthodes d’intégration 1D présentées précédemment.

3.2.2 Pour la matrice B

Z
La matrice B est la matrice de coefficients Bij = ik0 φi φj dΩ. On a considéré que le bord du domaine
∂Ω
Ω était un cercle, mais en réalité c’est un ensemble de N points Pi , i = 0, . . . , N − 1, avec N ≪ np , reliés un
à un par des arêtes de longueur L. Considérons deux points Pi et Pj consécutifs, de coordonnées respectives
(xi , yi ) et (xj , yj ). On a donc : q
Lij = (xi − xj )2 + (yi − yj )2
Considérons ces deux mêmes points Pi et Pj . Sur le bord de Ω, deux fonctions φi et φj ont un support commun
si |i − j| ≤ 1. Donc les termes non nuls de B seront Bii , Bij et Bij ′ où φj et φ′j sont les deux fonctions ayant
un support commun avec φi . On utilise les résultats de la partie 1D pour la matrice de masse et on obtient :
 Z Pi Z Pj !
L ′ Lij 
 2 2 j i
B = ik φ dΩ + φ dΩ = ik +

ii 0 0
 i i
3 3

Pj′


 Pi
 Z Pj
ik0 Lij

Bij = Bji = ik0 φi φj dΩ =

 Pi 6
Pi
 Z
ik 0 Lj ′ i


 Bij ′ = Bj ′ i = ik0 φi φj dΩ =


 ′
Pj 6

9
L’ensemble P2 (N −1) est l’ensemble des combinaisons sans répétition de deux éléments de l’ensemble {0, . . . , N −
1}. La matrice B aura alors la forme suivante :

ik0 X
B= Lij B(i, j)σ(i, j)
6
i,j∈P2 (N −1)

On a introduit l’ensemble P2 (N − 1) afin de ne pas prendre en compte les répétitions d’indices. La matrice
B(i, j) et la fonction σ(i, j) sont telles que :

i j
.
.. .. .. 
. .
i... 2 ... 1 ...
 

 .. . . ..  1 si Pi si Pj sont consécutifs
B(i, j) =  . . .  et σ(i, j) =
  0 sinon
j
... 1 ... 2 ... 
.. .. ..
. . .

3.2.3 Pour la matrice S

Z
La matrice S a pour coefficients Sij = ∇φi · ∇φj dΩ. Les fonctions φi et φj ont un support commun
Ω
si Pi et Pj sont consécutifs, donc s’ils appartiennent au même triangle. D’où :
 Z
∇φi · ∇φj dΩ si Pi si Pj sont consécutifs, appartiennent au même triangle K


Sij = K
 0 sinon


(0) (1)
Considérons un triangle d’indice T , où i ∈ {0, . . . , nt − 1}, et dont les trois sommets sont notés PT , PT et
(2)
PT . On a donc trois fonctions φ0 , φ1 et φ2 , qui sont telles que :
(j)

φi PT = δij , i, j ∈ {0, 1, 2}

Ces fonctions étant linéaires, on peut les exprimer comme suit :

φi (x, y) = ai + bi x + ci y

On en déduit donc que :

Z
∂φi ∂φj ∂φi ∂φj
∇φi · ∇φj = + = bi bj + ci cj =⇒ ∇φi · ∇φj dΩ = (bi bj + ci cj )AK
∂x ∂x ∂y ∂j K

où AK est l’aire du triangle K. Il reste à déterminer l’expression des coefficients bi et ci . Si l’on note Pi (xi , yi ), Pj (xj , yj )
et Pk (xk , yk ) les trois points du triangle K, on a :

ai + bi xj + ci yj = 0
φi (Pj ) = φi (Pk ) = 0 =⇒ =⇒ bi (xj − xk ) = ci (yk − yj )
ai + bi xk + ci yk = 0

=⇒ bi = λ(yk − yj ) et ci = λ(xj − xk ), λ ∈ R
On montre que le facteur de normalisation est λ = A2K . Les expressions de bj , bk , cj et ck sont obtenues par
permutation circulaire sur les indices i, j, k. La matrice S a donc pour expression :
X
S= S(i, j, k)AK
K

10
i, j et k sont les indices des points composant le triangle K. La matrice S(i, j, k) est telle que :
i j k
 .. .. .. 
. . .
i ... b2i + c2i ... bi bj + ci cj ... b i bk + c i c k ...
 
.. .. ..
 
 
 . . . 
b2j + c2j
 
j
S(i, j, k) =  ... bi bj + ci cj ... ... bj bk + cj ck ...

 .. .. .. 

 . . . 

k... bi bk + ci ck ... bj bk + cj ck ... b2k + c2k ...

.. .. ..
. . .

3.2.4 Pour la matrice M

Z
La matrice M a pour coefficients Mij = εr (r)φi φj dΩ. D’abord, on a :
Ω

ε si (x, y) ∈ Ωd
εr (x, y) =
1 si (x, y) ∈ Ω\Ωd

Étant donné que le maillage M recouvre aussi Ωd , un triangle K ∈ M ne peut appartenir à la fois à Ωd et
à Ω\Ωd . Cela signifie que εr (r) est constante sur tous les triangles du maillage. En raisonnant comme pour la
matrice S, on obtient donc :
 Z
 εr (K)

φi φj dΩ si Pi si Pj sont consécutifs, appartenant au même triangle K
Mij = K
 0 sinon


On utilise l’équation (10) pour en déduire que :

Z Z 1 Z 1−ξ  
Iij = φi φj dΩ = 2AK φi φj xi + (xj − xi )ξ + (xk − xi )η, yi + (yj − yi )ξ + (yk − yi )η dξdη
K 0 0
Z 1 "Z 1−ξ  
#
= 2AK φi φj xi + (xj − xi )ξ + (xk − xi )η, yi + (yj − yi )ξ + (yk − yi )η dη dξ
0 0
Z 1
= 2AK I(ξ)dξ
0   
AK 1
= I(0) + 4I + I(1)
3 2
On a utilisé la méthode de Simpson pour la dernière égalité. Déterminons I(ξ) à l’aide de la méthode du point
du milieu. On a :  
1−ξ
I(ξ) = (1 − ξ)φi φj ξ,
2
Cela implique que :
    
AK 1 1 1
Iij = φi φj 0, + 4φi φj ,
3 2 2 4
   
AK xi + xk yi + yk x
i xj xk yi yj yk 
= φi φj , + 4φi φj + + , + +
3 2 2 4 2 4 2 4 4
       
AK P i + Pk Pi + Pj Pj + Pk
= φi φj + 4 φi φj + φi φj
3 2 4 4
Le passage de la deuxième à la troisième ligne (pour le second terme) est justifié par le fait que les fonctions φ
sont linéaires. On a donc aussi φ(PA /n) = φ(PA )/n, ∀n ∈ R. D’où :
       
Pi + Pj Pi + Pj Pi + Pj 1 Pi + Pj
φi φj = φi φj = φi φj
4 4 4 4 2

11
On obtient donc :
        
AK Pi + P k Pi + Pj P j + Pk AK δij 1 AK
Iij = φi φj + φi φj + φi φj = + +0 = (δij + 1)
3 2 2 2 3 4 4 12

La matrice M s’exprime finalement comme suit :

X AK X AK
M =ε M (i, j, k) + M (i, j, k)
12 12
K∈Ωd K∈Ω\Ωd

Où la matrice M (i, j, k) est telle que :

i j k
 .. .. .. 
 . . . 
i... 2 ... 1 ... 1 ...
.. .. .. 
. . .
 
 
M (i, j, k) = j  . . . 1 ... 2 ... 1 ...
 
 .. .. .. 
. . .
 
 
k . . . 1 ... 1 ... 2 ...
 
.. .. ..
. . .

3.2.5 Pour la matrice C

Z
On a Ci = S(x, y)φi (x, y) dxdy. Soit Pi le point tel que φi (Pi ) = 1. Ce point est commun à plusieurs
Ω
triangles : 4 s’il est dans Ω et 2 s’il est au bord. On a donc :
X Z X (K)
Ci = Sφi dxdy = Ci
K,Ni ∈K K K,Ni ∈K

(K) 1
En utilisant les mêmes méthodes d’intégration que pour la matrice M , on montre que Ci = S(Pi )AK . Ainsi,
3
soient Pi , Pj et Pk les trois points formant le triangle K, on obtient :

X AK X AK  i j k
(K)

C= C = ... S(Pi ) . . . S(Pj ) . . . S(Pk ) . . .
3 3
K K

3.2.6 Validation du modèle

Pour le cas 2D, nous avons vérifié nos résultats sur un exemple calculé exactement par nos camarades du
groupe "Diffraction sur un cylindre". Nous avons effectué le calcul avec les mêmes paramètres et avons obtenus
les résultats présentés figure 3. On voit que les deux champs sont compatibles. Le fait que notre champ soit
moins marqué est une simple conséquence du fait qu’il s’agisse d’une approximation. Alors que dans le modèle
testé par nos collègues, il s’agit de l’équation exacte. Il nous faudrait donc augmenter la qualité du maillage
pour arriver à un résultat comparable. Cependant, cela a un coût élevé en performance, et nous voyons ainsi
que la limite de ce programme correspond aux limites de l’appareil utilisé.

12
 (a) Partie réelle du champ modèle (b) Partie imaginaire du champ modèle

Real part of the scattering field Complex part of the scattering field
1.5
1.2
750 750 1.2
0.9
0.9
500 0.6 500
0.6
250 0.3 250
0.3
0.0
y(nm)

0 0 0.0
0.3
250 250 0.3
0.6
500 500 0.6
0.9
0.9
750 1.2 750
1.2
1.5
1000 750 500 250 0 250 500 750 1000 1000 750 500 250 0 250 500 750 1000
x(nm)
(c) Partie réelle du champ testé (d) Partie imaginaire du champ testé

Figure 3 – Test de la validité du programme sur un exemple de paramètres : λ = 350 nm; r = 175 nm; ε = 9+2i
On voit que les deux champs sont identiques bien que le champ testé soit moins précis.

3.3 Convergence (qualitative)

Contrairement au cas 1D, il est difficile de tester la convergence de ce programme. De plus celle-ci dépend
de chaque objet et il nous est donc impossible de réaliser cette étude formellement. Nous allons donc étudier
qualitativement la convergence de ce programme sur un exemple. On représente sur la figure 4 les maillages
associés au champ obtenus après exécution du programme. On voit que la convergence s’obtient pour une
longueur caractéristique assez faible, avant laquelle le champ est peu représentatif de la réalité. Cela montre
une nouvelle fois la limite de ce programme. Heureusement, cette limite est en général atteignable pour un
ordinateur classique.

3.4 Résultats
L’avantage principal de ce programme, comparé à celui de nos camarades du projet "Diffraction par
un cylindre", est le fait de pouvoir faire diffracter n’importe quelle forme. On représente sur la figure 5 en
annexe, des exemples de résultats obtenus à partir de ce programme. Pour les obtenir, nous avons superposer
sur Inkscape, les cartes de champ obtenues sur Python, et l’objet dessiné sur Gmsh.

13
 Real part of the scattering field
3.0
2.4
200
1.8
100 1.2
0.6
0

y(nm)
0.0
100 0.6
1.2
200
1.8

300 2.4
200 100 0 100 200
x(nm)
(a) Paramètres : F = 2000, D = 500

Real part of the scattering field

3.6

200 2.7
1.8
100 0.9
0.0
0
y(nm)

0.9
100 1.8
2.7
200 3.6
4.5
300 200 100 0 100 200
x(nm)
(b) Paramètres : F = 200, D = 50

Real part of the scattering field

3.6
200
2.4

100 1.2

0 0.0
y(nm)

1.2
100
2.4
200
3.6
300 4.8
300 200 100 0 100 200
x(nm)
(c) Paramètres : F = 20, D = 5

Figure 4 – Comparaison des champs (à droite) obtenus à partir de maillages (à gauche) de différentes longueurs
caractéristiques
On a posé ici F = h(Ω\Ωd ) et D = h(Ωd ).

14
4 Conclusion
Notre but était de modéliser le champ électrique provenant d’une onde plane, après diffraction par un
objet de forme quelconque. A l’aide de ce programme, nous avons atteint cet objectif et avons présenté quelques
possibilités qu’il offre en Annexe. Nous avons vu qu’il a cependant ces limites puisqu’il s’agit d’une approxima-
tion, et que l’obtention d’une carte précise nécessite un ordinateur performant. Les application de cette méthode
sont diverses. Dans le domaine de l’électromagnétisme, on peut simuler une "cape d’invisibilité". Pour cela, on
cherche à donner l’impression que le champ électrique n’est pas affecté par la présence d’un objet diffractant
(celui que l’on cherche à rendre invisible). Outre l’électromagnétisme, cette méthode trouve de nombreux do-
maines d’application. Par exemple, lorsque l’on veut étudier l’effet de diverses contraintes appliquées à un objet
complexe. Nous pourrions étudier une de ces applications dans un prochain projet.

15
5 Annexe

Real part of the scattering field Complex part of the scattering field
3.2
3.00

1000 1000 2.4

2.25

1.50 1.6
500 500
0.75 0.8
y(nm)

0 0.00 0 0.0

0.75
0.8
500 500
1.50
1.6

1000 2.25 1000

2.4
3.00

1000 500 0 500 1000 1000 500 0 500 1000

x(nm)

(a) Paramètres : λ = 150 nm, θ = π/2, ε = 9 + 2j

Real part of the scattering field Complex part of the scattering field
3.0
1000 1000 1.80
2.4
1.35
800 800
1.8
0.90
600 1.2 600
0.45
0.6
y(nm)

400 400 0.00

0.0
200 200 0.45
0.6
0.90
0 1.2 0
1.35
1.8
200 200
1.80
2.4

0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200
x(nm)

(b) Paramètres : λ = 90 nm, θ = π, ε = 20 + 2j

Real part of the scattering field Complex part of the scattering field
4.8
3.6

200 200 3.6

2.4

2.4
100 1.2 100

1.2
0.0
y(nm)

0 0
0.0
1.2
100 100 1.2
2.4
2.4
200 200
3.6
3.6
300 4.8 300
300 200 100 0 100 200 300 200 100 0 100 200
x(nm)

(c) Paramètres : λ = 50 nm, θ = −π/4, ε = 9 + 2j

Real part of the scattering field Complex part of the scattering field
2.0
600 600 1.6
1.5
1.2
400 400
1.0
0.8

200 0.5 200

0.4
y(nm)

0 0.0 0 0.0

0.5 0.4
200 200

1.0 0.8
400 400
1.2
1.5
600 600 1.6
2.0
600 400 200 0 200 400 600 600 400 200 0 200 400 600
x(nm)

(d) Paramètres : λ = 200 nm, θ = −0.404, ε = 50 + 2j

Figure 5 – Cartes du champ électrique total obtenues pour différents objets, matériaux, angles d’incidence, et
longueurs d’onde

16