## MODÉLISATION

# PAR LA MÉTHODE DES ÉLÉMENTS FINIS AVEC ONELAB

**Guillaume DEMESY<fnref n="1" />*, André NICOLET<fnref n="1" />, Frédéric ZOLLA<fnref n="1" />, Christophe GEUZAINE<fnref n="2" />**

guillaume.demesy@fresnel.fr

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

Nous présentons ici le logiciel open source ONELAB de modélisation numérique par la méthode des éléments finis pour les applications photoniques. Nous illustrons à l'aide de quelques exemples une bibliothèque évolutive de modèles paramétrables couvrant une large gamme de dispositifs rencontrés en nanophotonique. Celle-ci permet d'aborder facilement la simulation d'applications réalistes tout en permettant au spécialiste de développer ses propres modèles avancés.

https://doi.org/10.1051/photon/202010040

Article publié en accès libre sous les conditions définies par la licence Creative Commons Attribution License CC-BY (http://creativecommons.org/licenses/by/4.0), qui autorise sans restrictions l'utilisation, la diffusion, et la reproduction sur quelque support que ce soit, sous réserve de citation correcte de la publication originale.

En ce début de 21ᵉ siècle, la technologie permet de façonner la matière avec une précision de quelques nanomètres, ce qui est bien inférieur aux longueurs d'onde de la lumière (en fait on parle ici aussi bien de la lumière visible que de l'infrarouge et de l'ultraviolet) qui peut ainsi être manipulée avec une subtilité inégalée jusqu'à présent. Tout comme la miniaturisation de l'électronique a donné lieu à la microélectronique, se dessine aujourd'hui le paysage très novateur de la nanophotonique. Pour y

concevoir de nouveaux systèmes, il faut prendre en compte la nature fondamentale de la lumière : c'est une onde électromagnétique interagissant avec la matière. L'étude de dispositifs tels que des réseaux de diffraction, des métasurfaces, des fibres optiques (qui peuvent être microstructurées) implique la prise en compte de géométries complexes où seul un outil numérique permet des calculs précis. En considérant des propriétés matérielles simples, il s'agit alors de résoudre une équation d'onde vectorielle pour le champ électrique (ou le champ magnétique).

Il n'existe pas de méthode numérique universellement reconnue comme étant meilleure que les autres mais ce que nous présentons ici est un outil très général et facile à utiliser, et ce, dans une large gamme de problèmes électromagnétiques tout en restant performant en temps de calcul. Cet outil logiciel est fondé sur la méthode des éléments finis et en particulier la méthode *finite element frequency domain* (FEFD). Un des principaux avantages de cette méthode est sa souplesse : les éléments finis permettent de représenter des géométries quelconques avec une grande

<footnote n="1">¹ Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France</footnote>

<footnote n="2">² Université de Liège, Institut Montefiore, Belgique</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.**
Modèle grating3D.pro : vecteur de Poynting dans une cellule solaire nanostructurée (à gauche) composée d'un réseau 2D de nano-tiges de sections circulaires. Une seule période est représentée. Le modèle permet de calculer l'efficacité de la cellule solaire. Le même modèle permet de simuler la réponse d'un réseau 2D de plots d'or à maille hexagonale (en bas).</figcaption></figure>

O X Y Z C 1:1 S K D D Done

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

précision (dans certaines struc-
tures, le diable se cache parfois
dans les détails...) et permet d'intro-
duire facilement des matériaux aux
propriétés complexes, anisotropes
et dispersifs en temps par exemple.
L'approche fréquentielle plutôt que
temporelle est, quant à elle, très
naturelle dans le contexte phy-
sique de la photonique où les pro-
priétés du système vont dépendre
de façon critique de la fréquence
de travail à travers les propriétés
des matériaux et les résonances du
système. La méthode des éléments
finis nécessite le maillage de la géo-
métrie en éléments simples (des
téträèdres, par exemple, mais qui
peuvent avoir des arêtes et des faces
courbes pour mieux s'adapter à la
géométrie réelle) et conduit à des
systèmes d'équations algébriques
de très grandes tailles mais égale-
ment très creux (la proportion de
coefficients non nuls est très faible).
Au cours des dernières décennies,
la puissance de calcul des calcula-
teurs a fortement augmenté mais
surtout différents obstacles s'op-
posant à l'utilisation pratique de
la FEFD pour la propagation des
ondes ont été surmontés de façon

satisfaisante : les outils de maillage,
les méthodes numériques pour les
très grands systèmes creux, le fait
de pouvoir ramener les sources à
distance finie dans la zone à mailler,
et les méthodes permettant de simu-
ler les conditions d'onde sortante,
ainsi que des interpolations spé-
cifiquement adaptées aux champs
vectoriels de l'électromagnétisme
(éléments d'arêtes) rendent cette
méthode actuellement très efficace.

LA MÉTHODE DES ÉLÉMENTS
FINIS EN RÉGIME FRÉQUENTIEL
(FEFD) POUR L'OPTIQUE

Techniquement, la méthode des élé-
ments finis repose sur deux piliers :
la formulation faible et une approxi-
mation de la solution par morceaux
à l'aide d'une décomposition de la
géométrie en parties finies (les élé-
ments finis) dont la taille est petite
par rapport à l'échelle de variation
spatiale des fonctions à approximer.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Modèle scattering.pro : onde plane illuminant un anneau torique fendu. La norme du champ électrique total est représentée en échelle de couleurs. Une coupe de la partie réelle du vecteur de Poynting est représentée sur un plan passant par l'entrefer. Une PML sphérique est implémentée dans le modèle.</figcaption></figure>

Commençons par la formulation faible : soit notre équation à résoudre qui est une équation aux dérivées partielles du second ordre, ici l'équation des ondes vectorielles dans un milieu inhomogène. Dans le cas de matériaux linéaires, en régime harmonique et loin des sources, cette équation des ondes vectorielles prend la forme suivante en termes du phaseur de champ électrique $E(x,y,z)$ :

$$ - \operatorname{rot} \mu^{-1} \operatorname{rot} \mathbf{E} + \omega^2 \varepsilon \mathbf{E} = 0, $$

où $\omega$ est la pulsation et où la perméabilité $\mu$ et la permittivité $\varepsilon$ dépendent de la position et peuvent être définies par morceaux pour modéliser l'assemblage de plusieurs matériaux. Pour obtenir la formulation faible, on multiplie les deux membres de l'équation par une fonction de pondération (une « fonction test ») et on intègre le résultat sur le domaine de définition du problème, ce qui fournit une égalité scalaire (un résidu pondéré). L'étape cruciale consiste alors à intégrer par parties l'expression trouvée de manière à reporter une partie de la dérivation de la fonction inconnue sur la fonction test, ce qui fait apparaître un terme sur le bord du domaine qui sera utile pour spécifier les conditions

aux limites, et permet de traiter de manière naturelle des discontinuités matérielles. On obtient ainsi une expression où les fonctions inconnues sont dérivées une seule fois, au lieu de deux, et l'exigence de régularité sur ces fonctions est ainsi affaiblie (d'où la notion de formulation faible).

Appliquons maintenant cette technique pour obtenir la méthode des éléments finis. On procède au maillage de la géométrie de notre problème, c'est-à-dire que l'on divise l'espace en tétraèdres (c'est un exemple mais c'est aussi la façon la plus courante de procéder) d'une manière telle que deux tétraèdres soient, soit disjoints, soit possèdent un sommet, une arête ou une face triangulaire en commun. Le maillage est alors dit conforme. On s'assurera également que les surfaces de discontinuité dans la géométrie du problème (une interface entre deux matériaux, par exemple) soient approximées correctement par des faces triangulaires du maillage. Il ne reste plus qu'à choisir, sur chaque tétraèdre, des fonctions très simples des coordonnées interpolant les grandeurs physiques inconnues, et à les injeter dans la formulation faible pour obtenir un système algébrique dont la

solution numérique donne l'approximation attendue.

Il reste à imposer les conditions aux limites correspondant au problème physique (ici aussi, les éléments d'arêtes ont des propriétés adéquates). Parmi les cas pratiques intéressants, on a le cas des structures périodiques que l'on modélise en imposant des conditions quasi-périodiques de Floquet-Bloch (cf. figure 3). Ces conditions permettent de réduire l'étude de la structure (infinie puisque périodique) à une seule cellule. Un autre problème est la prise en compte de l'espace libre. Une partie des ondes électromagnétiques sort de la zone d'intérêt et va se perdre à l'infini (considéré comme un milieu homogène) ce qui se traduit mathématiquement par une condition dite d'onde sortante. Le numéricien se trouve ici dans la même situation que l'expérimentateur dans sa cage de Faraday : pour travailler en espace confiné en ayant l'illusion d'un espace sans limite, il faut tapisser les bords de cet espace avec un matériau parfaitement absorbant. C'est ce qui est réalisé numériquement à l'aide des couches parfaitement adaptées (PML, *perfectly matched layers*, cf. figure 2) obtenues
en faisant un changement de coordonnées à valeurs complexes (pour rappel, nous sommes ici en régime harmonique et toutes nos grandeurs physiques sont complexes... c'est maintenant le tour des coordonnées !) qui fournissent des solutions souvent satisfaisantes. Il existe d'autres approches telles que les conditions aux limites absorbantes (ABC, absorbing boundary conditions).

Terminons en disant que l'on peut étudier la réponse d'une structure à l'illumination par une onde électromagnétique (un problème direct associé à la résolution d'un système linéaire) mais on peut aussi étudier des champs électromagnétiques propres à la structure, correspondant alors à des résonances ou à des propagations guidées en l'absence de toute source de champ. Le système algébrique obtenu par la méthode des éléments finis a pour cela une forme très agréable : il se met naturellement sous la forme d'un polynôme du second degré de la pulsation $\omega$ (on peut aussi prendre le nombre d'onde $k$ dans le cas des modes guidés) avec des coefficients matriciels, forme qui est très bien adaptée pour le calcul numérique des valeurs propres (et des vecteurs propres associés).

Pour une introduction courte mais détaillée avec des formules expliques, allant du cas élémentaire 1D à l'électromagnétisme, le lecteur pourra se reporter au chapitre correspondant dans la référence [3]. Plusieurs logiciels commerciaux sont fondés sur cette méthode FEFD. Parmi les plus populaires pour les applications photoniques, on peut citer COMSOL et Ansys HFSS.

## ONELAB

Nous présentons ici une alternative utilisant uniquement des logiciels open source, libres et gratuits, reposant sur la plateforme ONELAB (http://onelab.info).

ONELAB permet de piloter de manière interactive tous les éléments de la chaîne de calcul FEFD : définition

paramétrique des géométries, maillages, calculs et post-traitements. Pour les applications photoniques, deux logiciels principaux sont pilotés par ONELAB :

* Gmsh (http://gmsh.info [1]), pour la description paramétrique des géométries, le maillage et le post-traitement des résultats ;

* GetDP (http://getdp.info [2]), pour le calcul par éléments finis proprement dit.

Des scripts Python peuvent également être pilotés par ONELAB, par exemple pour certaines opérations répétitives, de post-traitement avancé, ou d'optimisation automatique. ONELAB est disponible sur toutes les plateformes de calcul (Windows, Linux et macOS), et une version mobile peut également être téléchargée pour tablette Android et iOS.

Une bibliothèque gratuite d'outils de modélisation spécifiques à la photonique a été développée par les auteurs de cet article pour cette plateforme, pouvant être utile aussi bien au chercheur académique qu'à l'ingénieur qui développe une application industrielle. L'intégralité des modèles, depuis la description géométrique jusqu'aux opérations de post-traitement, en passant par les équations aux dérivées partielles et les propriétés matériaux, sont décrites dans les fichiers d'entrée en format texte. La grande force de la bibliothèque est qu'elle peut être utilisée à plusieurs niveaux, du mode le plus élémentaire ne nécessitant (presque) aucune connaissance mathématique jusqu'au mode expert, où des connaissances à la fois des mathématiques sous-jacentes et de la syntaxe de certains métalangages sont nécessaires :

* on peut se contenter de faire valier des paramètres prédéfinis dans l'interface graphique, et procéder ainsi à une étude exploatoire interactive d'un type d'application (réseau de diffraction 2D ou 3D, fibre microstructurée...);

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

# L'OPTIQUE EST NOTRE AVENIR™

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

## Séparateurs de faisceaux par Edmund Optics®

* Large gamme de longueurs d'ondes allant de 250 nm à 8 µm

* Grande variété de formes, de tailles et de rapports de division

* Options de polarisation avec des rapports d'extinction allant jusqu'à 10 000:1

* Vaste inventaire de pièces standard disponibles pour une livraison rapide.

Edmund Optics® fabrique et fournit des séparateurs de faisceau standard et personnalisés pour une large gamme d'applications.

Pour en savoir plus :

www.edmundoptics.fr/beamsplitters

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

Rendez-vous à Strasbourg :

SPIE.

PHOTONICS
EUROPE

29.03 - 02.04.2020 | Stand 208

+33 (0) 820 207 555
sales@edmundoptics.fr

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

Edmund
optics | worldwide
<figure><img src="image_9.png" /><figcaption>Figure 3.
Mode de Bloch dans un cristal photonique et relation de dispersion (modèle rhombus.pro).</figcaption></figure>

* on peut également spécifier les paramètres en « batch » et faire tourner les simulations sur des supercalculateurs, pour des calculs de très grande taille ou pour de l'optimisation automatique ;

* enfin, en se plongeant dans le métalangage de Gmsh et GetDP, on peut modifier les exemples et définir soi-même ses propres modèles mathématiques, sans passer par un langage de bas niveau.

Tous les modèles décrits dans la suite sont disponibles sur la page http://onelab.info/photonics.

## DIFFRACTION

Un des problèmes les plus classiquement rencontrés en nanophotonique concerne la diffraction d'un champ électromagnétique par une structure dont la géométrie présente des dimensions de l'ordre de grandeur de la longueur d'onde incidente. On souhaite alors connaître le champ électromagnétique proche du motif diffractant, par exemple dans des applications d'exaltation de champ, et les grandeurs énergétiques associées au champ lointain, par exemple dans des applications de type filtrage spectral. Dans ce dernier cas l'utilisateur

souhaite connaître, dans le cas périodique les efficacités de diffraction et les pertes Joule, et dans le cas de résonateurs isolés le diagramme de rayonnement et les sections efficaces de diffusion, absorption, etc.

Par exemple, le modèle grating3D. pro est fourni avec une batterie de géométries pré-programmées, de sorte qu'il est possible de passer en quelques clics de l'étude de la réponse spectrale d'une cellule solaire nano-structurée (nanofils de silicium sur le volet de gauche de la figure 1) à celle d'un réseau à pavage hexagonal de plots d'or sur un substrat de verre (structure plasmonique sur volet de droite de la figure 1). Les détails plus théoriques concernant ces modèles sont disponibles dans les références [4,5]. Grâce aux éléments d'arête de haut ordre implémentés dans GetDP, il est possible d'atteindre de très bonnes précisions sur les grandeurs énergétiques d'intérêt.

Un autre exemple (scattering. pro) concerne les résonateurs isolés, comme l'anneau fendu présenté en figure 2. Ce modèle permet de déterminer la réponse en champ proche et en champ lointain d'un résonateur quelconque éclairé par une onde plane ou par un point source (il

s'agit alors de la fonction de Green), les quasi-modes du résonateur ou encore sa matrice T.

## ANALYSES MODALES

Nous présentons deux autres exemples afin d'illustrer les capacités de ONELAB pour des problèmes modaux. L'exemple rhombus.pro concerne le calcul des modes propagatifs dans un cristal photonique constitué de tiges diélectriques infinies parallèles de section circulaire. On se ramène à l'étude de la section droite d'une seule cellule avec conditions de Floquet-Bloch. Bien que de géométrie bi-dimensionnelle, les modes sont des champs de vecteurs à trois composantes (cas conique, cf. figure 3 où neuf cellules élémentaires sont représentées pour illustrer la nature quasi-périodique du champ). La relation de dispersion entre le vecteur d'onde de Bloch et la fréquence est tracée en temps réel dans ONELAB et fait apparaître une bande interdite photonique. Le mode représenté peut être sélectionné interactivement en double-cliquant sur le point correspondant des courbes de dispersion.

L'exemple NonLinearEVP.pro illustre la capacité de ONELAB à résoudre des problèmes modaux en
<figure><img src="image_10.png" /><figcaption>Figure 4. Mode fondamental d'un réseau de tiges circulaires de permittivité dispersive régie par un modèle de Drude (modèle **NonLinearEVP.pro**). Le problème aux valeurs propres est alors non-linéaire.</figcaption></figure>

présence de matériaux dispersifs, c'est-à-dire dont la permittivité dépend de la valeur propre recherchée. Le problème aux valeurs propres résultant est alors non-linéaire. La figure 4 montre le mode fondamental d'un réseau de diffraction constitué de tiges faites d'un métal de Drude. Ces calculs de modes quasi-normaux avec des matériaux dispersifs (DQNM [7]) ont fait l'objet d'importants développements ces dernières années et sont aujourd'hui à la pointe de la modélisation numérique en photonique.

## CONCLUSION

La présentation succincte de quelques modèles illustre les capacités de ONELAB pour la modélisation en photonique. Pour aller plus loin, le lecteur est invité à consulter les références ci-dessous et surtout à télécharger et à tester le code sur sa machine favorite.

## REMERCIEMENTS :

Nous tenons à remercier tous les développeurs et la communauté ONELAB, GetDP et Gmsh, avec une pensée particulière à Patrick Dular, co-créateur de GetDP qui nous a quittés bien trop tôt. Les auteurs remercient l'ANR dans le cadre du financement ANR-16-CE24-0013, de même que la Région Wallonne et la Fédération Wallonie-Bruxelles de Belgique pour les financements WIST3 No 1017086 ONELAB t ARC WAVES 15/19-03.

## RÉFÉRENCES

[1] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. *International journal for numerical methods in engineering* 79 (2009): 1309-1331.

[2] P. Dular, C. Geuzaine, F. Henrotte, W. Legros, *A general environment for the treatment of discrete problems and its application to the finite element method*. *IEEE Transactions on Magnetics* 34 (1998): 3395-3398.

[3] F. Zolla, G. Renversez, A. Nicolet, *Foundations of photonic crystal fibres*. *World Scientific* (2005).

[4] G. Demésy, F. Zolla, A. Nicolet, *A ONELAB model for the parametric study of mono-dimensional diffraction gratings*, arXiv:1710.11451

[5] G. Demésy, S. John, *Solar energy trapping with modulated silicon nanowire photonic crystals*. *Journal of Applied Physics* 112 (2012): 074326.

[6] N. Marsic, H. De Gersem, G. Demésy, A. Nicolet, C. Geuzaine, *Modal analysis of the ultrahigh finesse Haroche QED cavity*. *New Journal of Physics* 20 (2018): 043058.

[7] P. Lalanne et al. *Quasinormal mode solvers for resonators with dispersive materials*. *JOSA A* 36 (2019): 686-704.

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