# Modélisation des interactions champs - structures par une méthode de Galerkin discontinue pour des applications du génie électrique

Abelin Kameni

► To cite this version:

Abelin Kameni. Modélisation des interactions champs - structures par une méthode de Galerkin discontinue pour des applications du génie électrique. Energie électrique. Université de Paris-Saclay, 2019. tel-04534218

HAL Id: tel-04534218

https://theses.hal.science/tel-04534218

Submitted on 5 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
UNIVERSITÉ PARIS-SUD

École doctorale EOBE (ED 575)

Laboratoire Génie électrique et électronique de Paris (UMR 8507 CNRS)

Mémoire présenté pour l'obtention du

# Diplôme d'habilitation à diriger les recherches

Discipline : Génie électrique

par

Abelin KAMENI NTICHI

Modélisation des interactions champs - structures par une méthode de Galerkin discontinue pour des applications du génie électrique

Travaux soutenus le 28 juin 2019 devant le jury composé de :

<table><tr><td>CHRISTOPHE GEUZAINE</td><td>Professeur des universités</td><td>Examinateur</td></tr><tr><td>JEAN LEVEQUE</td><td>Professeur des universités</td><td>Examinateur</td></tr><tr><td>LIONEL PICHON</td><td>Directeur de recherche</td><td>Examinateur</td></tr><tr><td>FRANÇIS PIRIOU</td><td>Professeur des universités</td><td>Rapporteur</td></tr><tr><td>ALAIN REINEIX</td><td>Directeur de recherche</td><td>Rapporteur</td></tr><tr><td>HÉLÈNE ROUSSEL</td><td>Professeur des universités</td><td>Examinatrice</td></tr><tr><td>BENOIT VANDERHEYDEN</td><td>Professeur des universités</td><td>Rapporteur</td></tr></table>
# Table des matières

## 1 Informations Générales 3

<table><tr><td>1.1</td><td>Fiche synthétique</td><td>4</td></tr><tr><td>1.2</td><td>Curriculum Vitae</td><td>5</td></tr><tr><td>1.3</td><td>Thèse de doctorat</td><td>6</td></tr><tr><td>1.4</td><td>Résumé des Activités d'enseignement</td><td>7</td></tr><tr><td>1.4.1</td><td>Contenu des enseignements</td><td>7</td></tr><tr><td>1.4.2</td><td>Activités administratives</td><td>8</td></tr><tr><td>1.5</td><td>Résumé des activités de recherche</td><td>9</td></tr><tr><td>1.5.1</td><td>Contexte</td><td>9</td></tr><tr><td>1.5.2</td><td>Synthèse des contributions</td><td>10</td></tr><tr><td>1.5.3</td><td>Collaborations scientifiques</td><td>15</td></tr><tr><td>1.5.4</td><td>Activités d'intérêt collectif</td><td>17</td></tr><tr><td>1.5.5</td><td>Activités d'encadrement</td><td>17</td></tr><tr><td>1.5.6</td><td>Publications scientifiques</td><td>19</td></tr></table>

## 2 Calcul des champs rayonnés en compatibilité électromagnétique 28

<table><tr><td>2.1</td><td>Introduction</td><td>29</td></tr><tr><td>2.2</td><td>Formulation du système différentiel</td><td>31</td></tr><tr><td>2.2.1</td><td>Les équations de Maxwell</td><td>31</td></tr><tr><td>2.2.2</td><td>Formulation faible du problème</td><td>31</td></tr><tr><td>2.3</td><td>Flux caractéristiques aux interfaces</td><td>32</td></tr><tr><td>2.3.1</td><td>Invariants de Riemann</td><td>33</td></tr><tr><td>2.3.2</td><td>Flux numériques</td><td>35</td></tr><tr><td>2.3.3</td><td>Conditions aux bords</td><td>37</td></tr><tr><td>2.4</td><td>Condition d'absorption pour les problèmes de diffraction</td><td>38</td></tr><tr><td>2.4.1</td><td>Les parois absorbantes</td><td>39</td></tr><tr><td>2.4.2</td><td>Les couches parfaitement absorbantes</td><td>40</td></tr><tr><td>2.4.3</td><td>Optimisation des couches PML</td><td>41</td></tr><tr><td>2.5</td><td>Schéma numérique pour l'évolution temporelle</td><td>49</td></tr><tr><td>2.5.1</td><td>Les méthodes de Runge-Kutta implicites et explicites</td><td>49</td></tr><tr><td>2.5.2</td><td>Runge-Kutta explicite avec pas de temps local</td><td>51</td></tr><tr><td>2.6</td><td>Modélisation de boitiers de blindage</td><td>57</td></tr><tr><td>2.6.1</td><td>Interface simulant les plaques minces homogènes</td><td>57</td></tr><tr><td>2.6.2</td><td>Modélisation 3D de boitiers à parois composites</td><td>63</td></tr></table>
<table><tbody><tr><td>2.7</td><td>Diagnostic de réseau filaire</td><td>70</td></tr><tr><td>2.7.1</td><td>Modèle de réflectométrie pour la caractérisation des défauts non-francs</td><td>71</td></tr><tr><td>2.7.2</td><td>Evaluation des paramètres électriques équivalent</td><td>77</td></tr><tr><td><b>3</b></td><td><b>Calcul des champs induits dans les matériaux supraconducteurs</b></td><td><b>85</b></td></tr><tr><td>3.1</td><td>Introduction</td><td>86</td></tr><tr><td>3.2</td><td>Modélisation numérique des supraconducteurs</td><td>88</td></tr><tr><td>3.2.1</td><td>Loi de comportement électrique des supraconducteurs HTC</td><td>88</td></tr><tr><td>3.2.2</td><td>Formulations pour le calcul de champs électromagnétiques dans les supraconducteurs HTC</td><td>88</td></tr><tr><td>3.2.3</td><td>Linéarisation de la loi constitutive</td><td>90</td></tr><tr><td>3.3</td><td>Méthodes avec flux pour le calcul de champs dans les supraconducteurs</td><td>91</td></tr><tr><td>3.3.1</td><td>Méthode des volumes finis</td><td>92</td></tr><tr><td>3.3.2</td><td>Méthode de Galerkin discontinue</td><td>98</td></tr><tr><td>3.4</td><td>Modélisation des câbles supraconducteurs torsadés</td><td>114</td></tr><tr><td>3.4.1</td><td>Modèle équivalent des filaments torsadés</td><td>115</td></tr><tr><td>3.4.2</td><td>Réduction du modèle équivalent</td><td>120</td></tr><tr><td><b>4</b></td><td><b>Perspectives</b></td><td><b>127</b></td></tr><tr><td>4.1</td><td>Introduction</td><td>128</td></tr><tr><td>4.2</td><td>Modélisation de systèmes complexes</td><td>128</td></tr><tr><td>4.2.1</td><td>Méthode de décomposition de domaine</td><td>129</td></tr><tr><td>4.2.2</td><td>Méthode d'analyse d'incertitudes</td><td>129</td></tr><tr><td>4.3</td><td>Transmission d'énergie pour l'alimentation des systèmes</td><td>130</td></tr><tr><td>4.3.1</td><td>Transmission d'énergie sans fil</td><td>131</td></tr><tr><td>4.3.2</td><td>Transmission d'énergie sans pertes : câbles supraconducteurs</td><td>132</td></tr></tbody></table>
# Chapitre 1

## Informations Générales

### Sommaire

<table><tr><td>1.1 Fiche synthétique</td><td>4</td></tr><tr><td>1.2 Curriculum Vitae</td><td>5</td></tr><tr><td>1.3 Thèse de doctorat</td><td>6</td></tr><tr><td>1.4 Résumé des Activités d'enseignement</td><td>7</td></tr><tr><td>1.4.1 Contenu des enseignements</td><td>7</td></tr><tr><td>1.4.2 Activités administratives</td><td>8</td></tr><tr><td>1.5 Résumé des activités de recherche</td><td>9</td></tr><tr><td>1.5.1 Contexte</td><td>9</td></tr><tr><td>1.5.2 Synthèse des contributions</td><td>10</td></tr><tr><td>1.5.3 Collaborations scientifiques</td><td>15</td></tr><tr><td>1.5.4 Activités d'intérêt collectif</td><td>17</td></tr><tr><td>1.5.5 Activités d'encadrement</td><td>17</td></tr><tr><td>1.5.6 Publications scientifiques</td><td>19</td></tr></table>
1.1 Fiche synthétique

Abelin KAMENI NTICHI - 39 ans

Maître de conférences à l'Université Paris-Sud

Chercheur au Laboratoire de Génie électrique et électronique de Paris (GeePs)

Formations et diplômes

2005-2009 : Thèse de doctorat en Génie électrique, université Henri-Poincaré de Nancy

2004-2005 : Master Recherche de Mathématiques appliquées, université de Bordeaux I

2000-2004 : Ingénieur MATMECA, université de Bordeaux I

Activités d'enseignement

Depuis 2010 : Université de Paris-Sud 192h/an
UFR des Sciences et Polytech'Paris-Sud
Génie électrique, Electronique, Mécanique, Simulation numérique

2005-2009 : Université Henri-Poincaré de Nancy 300h
Génie électrique, Mécanique

Activités de recherche

Calcul de champs électromagnétiques pour le génie électrique

Mes travaux de recherches portent sur l'évaluation des champs électromagnétiques dans les dispositifs du génie électrique, plus particulièrement sur l'étude des problèmes de compatibilité électromagnétique rayonnée et sur la caractérisation des matériaux supra-conducteurs. Ils abordent à la fois, les aspects numériques à travers l'implémentation des outils numériques pour la simulation, et les applications par la mise en place des modèles permettant de répondre à certaines problématiques industrielles.

Publications

Articles dans revues internationales : 26

Actes de colloques internationaux avec comité de lecture : 37

Colloques nationaux avec actes : 7

Activités d'intérêt collectif

Responsable du recrutement des doctorants enseignants en Physique et Sciences de l'ingénieur.

Coordinateur de la 4ième année du cycle d'ingénieur de Polytech'Paris-Sud.

Rapporteur occasionnel pour les revues IEEE Tran. Magn, IEEE Appl. Supercond, IEEE EMC.
1.2 Curriculum Vitae

Abelin KAMENI NTICHI né le 04 novembre 1979 (39 ans)

Adresse professionnelle : Laboratoire Génie électrique et électronique de Paris (GeePs),
CNRS/ CentraleSupelec/ Univ. Paris-Sud/ Université Paris-Saclay/ Sorbonne Université,
3 & 11 rue Joliot-Curie, Plateau du Moulon, 91192 Gif-Sur-Yvette Cedex

Tel : 01.69.85.16.66 - Fax : 01.69.41.83.18

Maître de conférences à l'université de Paris-Sud

Chercheur au laboratoire Génie Electrique et Electronique de Paris (GeePs)

Formation

**2005 - 2009 Doctorat en génie électrique**
Groupe de Recherche en Electronique et Electrotechnique de Nancy (GREEN)

**2004 - 2005 Master Recherche de Mathématiques Appliquées**
Université de Bordeaux I

**2003 - 2004 Master Recherche de Mécanique Physique**
Université de Bordeaux I

**2000 - 2004 Cycle ingénieur MATMECA**
Ecole d'ingénieur en Modélisation Mathématique et Mécanique
Université de Bordeaux I

**1997 - 2000 Classes préparatoires scientifiques**
Lycée Kérichen de Brest

Expérience professionnelle

**Depuis 1/09/2010 : Maître de conférences de l'université de Paris-Sud**
Chercheur au laboratoire Génie électrique et électronique de Paris (GeePs)
Pôle physique et ingénierie de l'électromagnétisme (PIEM)

depuis 09/2013 Responsable du recrutement des doctorants enseignants en Physique/SI

depuis 02/2018 Responsable de la 4ième année du cycle ingénieur de Polytech 'Paris-Sud

**10/2009-09/2010 : Ingénieur de Recherche à l'université de Liège, Belgique**
Institut Montéfiore, équipe Applied Computational Electromagnetics (ACE)

**10/2008-08/2009 : Attaché Temporaire d'Education et de Recherche (ATER)**
Université Henri Poincaré de Nancy - GREEN

**10/2005-09/2008 : Allocataire de Recherche (MENRT) et Doctorant enseignant (Moniteur)**
Université Henri Poincaré de Nancy - GREEN
1.3 Thèse de doctorat

Période : 2005-2009

Spécialité : Génie électrique

Laboratoire : Groupe de Recherche en Electronique et Electrotechnique de Nancy (GREEN)

Titre : Etude par des méthodes analytiques et numériques de la répartition des champs induits dans les supraconducteurs à haute température critique

Soutenance : 24 juin 2009

Composition du jury :

Rapporteurs : Frédéric BOUILLAULT, Professeur, Université de Paris-Sud
Gérard MEUNIER, Directeur de recherche CNRS, ENSIEG-INPG Grenoble

Examinateurs : Christophe GEUZAINE, Professeur, Université de Liège
Bertrand DUTOIT, Docteur, Ecole Polytechnique Fédérale de Lausanne

Encadrants : Denis NETTER, Professeur, Université Henri-Poincaré de Nancy
Bruno DOUINE, Maître de conférences, Université Henri-Poincaré de Nancy

**Résumé** : Les travaux réalisés dans cette thèse ont pour but de décrire l'évolution des champs induits dans les matériaux supraconducteurs à haute température critique dont le comportement électrique est donné par une relation de type puissance E ~ Jⁿ. Ces travaux s'articulent autour de la résolution d'un problème de diffusion non linéaire issu du couplage de cette loi constitutive aux équations de Maxwell. Dans une première partie, on utilise une approche analytique basée sur le principe d'auto-similarité pour résoudre ce problème en 1D et caractériser la pénétration de la densité de courant dans une plaque supraconductrice. Cette solution nous permet de valider la méthode numérique mixte éléments finis-volumes finis (FEM-FVM) proposée pour les cas bi-dimensionnels. Dans cette approche numérique, une technique de linéarisation est introduite pour faire face aux difficultés dues à exposant n de la caractéristique E(J). Grâce à cette approche numérique on a pu mettre en évidence l'influence de la vitesse de variation d'une source imposée sur l'aimantation et la dissipation dans ces matériaux.

**Mots clés** : Equations de Maxwell, supraconducteurs à haute température critique, aimantation, pertes, volumes finis, éléments finis, auto-similarité.

**Valorisation** : [RI.27], [RI.26], [RI.25]
# 1.4 Résumé des Activités d'enseignement

J'ai découvert l'enseignement universitaire durant ma thèse de doctorat par le biais d'un contrat de moniteur portant sur un service de 64h/an. Durant ces 3 ans de monitorat, j'ai dispensé des Travaux Dirigés (TD) et Travaux Pratiques (TP) aux étudiants de Licence et Master du parcours Electricité Electronique et Automatique (EEA) de l'UFR des Sciences et Technique de l'université de Nancy. Cette expérience en tant que moniteur m'a permis de me faire une première idée du métier d'enseignant-chercheur, car il a fallu combiner travaux de recherche et mission d'enseignement. Par la suite, j'ai occupé pendant une année universitaire un demi-poste d'Attaché Temporaire d'Enseignement et de Recherche (ATER) portant sur un service de 96h/an.

Après un séjour post-doctoral en Belgique, j'ai été recruté comme maître de conférences à l'université de Paris-Sud avec affectation au département de physique de l'UFR des sciences pour la partie enseignement. Mon service de 192h/an se répartit entre Cours (CM), travaux Dirigés (TD), travaux Pratiques (TP) et Charges pédagogiques. Les étudiants concernés sont essentiellement ceux du cycle Licence.

## 1.4.1 Contenu des enseignements

- UFR des sciences et Techniques de l'université de Nancy 2005 – 2009 (288h/Total)

Parcours : Electricité Electronique et Automatique (L1, L2, L3, M1)

<table><thead><tr><td>Matières</td><td>Public</td><td>Nombre d'heures</td><td>Type de séances</td></tr></thead><tbody><tr><td>Circuits électriques</td><td>L1, L2</td><td>96h</td><td>TD, TP</td></tr><tr><td>Electromagnétisme</td><td>L2</td><td>48h</td><td>TD</td></tr><tr><td>Mécanique du point</td><td>L2</td><td>60h</td><td>TD</td></tr><tr><td>Electronique de puissance</td><td>L3</td><td>48h</td><td>TP</td></tr><tr><td>Outils informatique</td><td>M1</td><td>36h</td><td>TP</td></tr></tbody></table>

- UFR des sciences de l'université de Paris Sud, depuis 2010 (≥ 192h/an)

Parcours : 
Electricité Electronique et Automatique (L3 EEA)
Physique Chimie (L3PC)
Techniques de Production de l'Energie (L3pro TPE)
Physique Chimie Sciences de la Terre (L1 PCST)
Maths Physique Informatique (L1 MPI)

- Ecole d'ingénieur Polytech'Paris-Sud, depuis 2016 (32h/an)

Parcours : Cycle ingénieur 3ième année (Et3)
<table><thead><tr><th>Matières</th><th>Public</th><th>Nombre d'heures</th><th>Type de séances</th></tr></thead><tbody><tr><td>Electronique analogique</td><td>L3</td><td>96h</td><td>TD, TP</td></tr><tr><td>Mécanique du point</td><td>L1</td><td>52h</td><td>TD</td></tr><tr><td>Simulation numérique</td><td>L3, L3pro</td><td>64h</td><td>CM, TD, TP</td></tr></tbody></table>

<table><thead><tr><th>Matières</th><th>Public</th><th>Nombre d'heures</th><th>Type de séances</th></tr></thead><tbody><tr><td>Introduction aux ondes</td><td>Et3</td><td>32h</td><td>CM, TD</td></tr></tbody></table>

## 1.4.2 Activités administratives

1. Responsable du recrutement des Doctorants chargé d'une mission d'enseignement en Physique et Sciences de l'ingénieur (depuis 2013)

▷ Collecte des dossiers de candidature

▷ Organisation de la commission de sélection des candidats

▷ Gestion et suivi des affectations

▷ Interlocuteur Doctorants enseignants – Service administratif

2. Coordinateur 4ième année de l'école d'ingénieur Polytech' Paris-sud (depuis 2017)

▷ Mise en place et gestion de l'emploi du temps

▷ Gestion des réservations des salles d'enseignements

▷ Création et répartition des groupes TD et TP

▷ Organisation des jurys de semestres

▷ Interlocuteur Etudiants – Equipe de direction - Enseignants
**1.5 Résumé des activités de recherche**

J'effectue mes travaux de recherche au laboratoire Génie électrique électronique de Paris (GeePs) dans le pôle Physique et Ingénierie de l'Electromagnétisme (PIEM) sur le domaine de la modélisation des interactions entre les champs, les matériaux et les systèmes. Le fil conducteur de mes travaux de recherche est la mise en place de modèles numériques pour le calcul des champs électromagnétiques dans les dispositifs du génie électrique. Dans ce qui suit, je présente le contexte général de mes activités de recherche et les principales thématiques abordées.

**1.5.1 Contexte**

Avec l'augmentation de la population mondiale, les dispositifs électromagnétiques présents dans les domaines, du transport, de la communication, ou encore de la santé sont de plus en plus complexes. Leur élaboration nécessite la réalisation et la maitrise des étapes telles que, la caractérisation des matériaux utilisés, le dimensionnement et la fabrication des prototypes, ou encore les tests d'utilisation dans un environnement réel. Ce travail ne peut se faire uniquement sur la base de l'expérimentation car trop coûteuse en temps et en ressources. La mise en place des modèles de simulation pour aider à la caractérisation des matériaux, à la conception des systèmes, au dimensionnement des prototypes, est devenu indispensable pour le développement des applications du génie électrique. Concevoir des modèles de simulation, en s'appuyant d'une part sur une description mathématique des phénomènes physiques mis en jeu et d'autre part sur les outils et méthodes numériques permettant de résoudre les systèmes d'équations obtenus, est un enjeu important auquel doivent faire face les chercheurs, les ingénieurs et les industriels. Le développement et la diversification des outils de simulation numérique dans le domaine de l'électromagnétisme appliqué constitue aujourd'hui un axe important de la recherche en génie électrique. Le laboratoire de Génie électrique et électronique de Paris (GeePs) auquel j'appartiens est un acteur majeur de la recherche dans ce domaine en France et jouit d'une certaine renommée internationale. Les travaux qui y sont développés concernent les domaines de l'électronique, l'électromagnétisme et l'énergie. Depuis plusieurs années, ce laboratoire a acquis une expertise en modélisation dont on se sert pour traiter diverses problématiques industrielles. Celles liées à l'électromagnétisme sont traitées dans le pôle PIEM, dont les activités vont du matériau au système en couvrant les domaines de fréquences s'étalant du statique au TeraHertz [32, 118]. Dans ce pôle sont menés des travaux autour de la modélisation et la caractérisation des matériaux tels que, les matériaux magnétiques [33, 97], les matériaux composites[81, 85], ou encore des matériaux supraconducteurs [RI.16, RI.5]. Autour de ces travaux sur les matériaux, des études sont réalisées sur diverses problématiques, par exemple sur le contrôle électromagnétique des matériaux [31], la compatibilité électromagnétique des systèmes [7], ou encore l'exposition de
l'homme au rayonnement électromagnétique [34, 27], avec comme objectif final la réalisation de dispositifs permettant de répondre à des besoins industriels spécifiques [49]. Les travaux que je mène au sein du GeePs s'inscrivent parfaitement dans les problématiques traitées dans le pôle PIEM. La diversité des sujets traités dans ce pôle, des méthodes employées et des chercheurs avec qui je discute au quotidien constitue un apport considérable dont je bénéficie. Durant ces années passées au GeePs, j'ai ainsi pu acquérir une certaine expertise sur la modélisation des interactions entre les champs électromagnétiques et les structures.

### 1.5.2 Synthèse des contributions

Ma contribution porte particulièrement sur la modélisation des interactions entre les champs électromagnétiques avec les matériaux et les structures. Elle s'appuie essentiellement sur des modèles d'évaluation des champs électromagnétiques. Mon travail est axé sur le développement des outils numériques et la mise en place des modèles d'applications en relation avec des problématiques industrielles. Les développements réalisés concernent autant des applications en haute fréquence, comme par exemple la protection des systèmes d'électronique embarquée, que des applications en basse fréquence comme par exemple les câbles supraconducteurs de transport d'énergie électrique.

L'ensemble des travaux réalisés en haute fréquence nécessite l'évaluation des champs rayonnés, thème sur lequel j'ai travaillé et apporté des contributions notamment sur des problématiques liées à la compatibilité électromagnétique [RI.11, RI.21], le contrôle électromagnétique des systèmes [RI.7, RI.4, CI.9] ou encore les télécommunications [RI.8, RI.3]. Quant aux travaux réalisés en basse fréquence, ils se situent autour des matériaux supraconducteurs [RI.17, RI.5, RI.6], avec comme principal objectif, l'évaluation des grandeurs induites nécessaires au dimensionnement des dispositifs tels que les câbles de transport d'énergie, les moteurs supraconducteurs, les limiteurs de courant ou encore les aimants à forts champs utilisés dans certains dispositifs médicaux (appareils IRM).

#### Champs rayonnés : compatibilité et contrôle électromagnétique des systèmes.

Avec le développement des nouvelles technologies de plus en plus complexes, les problématiques liées à la maintenance et à la cohabitation des systèmes deviennent des points clés à maitriser pour s'assurer du bon fonctionnement des dispositifs. La maintenance d'un système passe par la détection, la localisation et la réparation des pannes qu'il peut subir. La maintenance peut s'avérer être cruciale, notamment lorsque la sécurité des personnes est en jeu. La maintenance hors fonctionnement des systèmes constitue le moyen classique pour éviter les pannes entrainant les dysfonctionnements. Malheureusement démonter l'ensemble du dispositif pour le contrôler s'avère très couteux voire impossible à réaliser dans de nombreux cas. Les techniques
de contrôle électromagnétique peuvent être mises en place afin de localiser la panne et ainsi
cibler la zone d'intervention, ce qui permet d'éviter le démontage total du système. Lorsque
les dispositifs cohabitent dans un même environnement, ils se comportent comme émetteurs et
récepteurs de signaux électromagnétiques pouvant créer des dysfonctionnements allant jusqu'à
l'arrêt complet du système.

Il est clair que l'étude des interactions entre les champs électromagnétiques avec les systèmes est au coeur des thématiques de contrôle et compatibilité électromagnétique. C'est pour répondre à des besoins industriels bien identifiés autour de ces thèmes que je développe des modèles numériques de calcul de champs électromagnétiques rayonnés.

La partie centrale de mes travaux dans sur ces thèmes porte sur le développement d'outils numériques pour l'évaluation des champs électromagnétiques en haute fréquence via la résolution des équations de Maxwell. Les méthodes classiques telles que la méthode des différences finies [53], la méthode des éléments finis [76], la méthode des volumes finis [83] ou encore la méthode des moments [120], sont celles généralement implémentées pour mener cette résolution. Avec l'augmentation de la taille des systèmes et la complexification des dispositifs à modéliser, ces approches numériques peuvent dans certain cas s'avérer inadaptées pour mercer à bien les simulations numériques souhaitées. En effet lorsque la résolution nécessite des inversions matricielles, comme c'est le cas avec la méthode des éléments finis et la méthode des moments, le coût de calcul est très impacté par l'augmentation du nombre d'inconnues du système. Pour réduire les coûts notamment des inversions matricielles, on a généralement recours à des développements supplémentaires, parmi lesquelles on peut citer le calcul parallèle sur architecture multi-coeurs de type CPU ou GPU [6, 13] et les techniques de compression matricielle [16]. La méthode des différences finies qui est la méthode la plus utilisée pour la simulation des phénomènes de propagation d'ondes en domaine temporel ne requiert pas d'inversion matricielle. Elle est généralement associée à des schémas numériques de type sautemouton et constitue une approche numérique rapide et simple à mettre en oeuvre [53]. Elle a été généralisée sur des maillages duaux et est connue sous le nom de "Finite Integration Technique" [28, 121]. Mais, l'efficacité des approches différences finies repose sur l'utilisation des maillages cartésiens réguliers qui malheureusement sont peu adaptés pour traiter les géométries complexes présentant notamment des courbures. Dans le cas de ces géométries complexes, pour préserver la convergence, la précision et la rapidité de ce type de méthode on a recourt à des développements supplémentaires comme par exemple, le raffinement local [30], l'élévation de l'ordre [122], la parallélisation sur architecture CPU ou GPU [71, 68].

D'une manière générale, on constate que l'utilisation des approches classiques que sont les
méthodes des éléments finis, des différences finies, des volumes finis et aussi celle des mo-
ments pour modéliser des systèmes complexes nécessitent souvent d'important développement
supplémentaires qui dans certains cas sont très difficiles à implémenter. Mettre en place de tels développements demande de monopoliser d'important moyens humains et matériels dans le domaine du calcul scientifique haute performance, ce qui est du ressort des laboratoires et équipes de recherche dont les activités et le personnel ont une forte orientation vers les mathématiques appliqués.

Mes travaux sur le développement des outils numériques s'appuient sur des collaborations avec des équipes de recherche plutôt très "numériciennes", notamment avec, l'équipe "Applied Computational Electromagnetics" du Pr Geuzaine de l'université de Liège et l'équipe "Applied mechanics and mathematics" du Pr Remacle de l'université de Louvain-la- neuve. Ces collaborations font suite à mon post-doc effectué à l'université de Liège, durant lequel j'ai participé à la mise en place d'outils de modélisation basés sur les méthodes éléments finis continus et discontinus. Ayant participé durant ce post-doc au développement d'un code de calcul Galerkin discontinue, j'ai adopté cette méthode pour les simulations des phénomènes de propagation en domaine temporel car elle présente de nombreux avantages, elle est notamment bien adaptée à l'utilisation des maillages non conformes et à la mise en place du calcul parallèle [51].

Cette technique qui a connu un essor important avec le développement des architectures multicoeurs des ordinateurs présente alors un intérêt certain pour la modélisation des grands systèmes. Elle a été introduite initialement pour la résolution de l'équation de transport en neutronique [82], puis a été implémentée en mécanique des fluides pour la résolution des équations de conservation [29, 87]. Son utilisation dans le domaine du génie électrique s'est intensifiée durant la dernière décennie. Cette approche combine les règles de discrétisation des méthodes éléments finis et volumes finis. Elle consiste à résoudre le problème localement sur chaque élément du maillage par une approche éléments finis, et à utiliser des termes de flux aux interfaces pour décrire les échanges entre éléments comme en volumes finis. L'aspect résolution locale sur chaque élément du maillage, simplifie l'utilisation des éléments d'ordre élevé ce qui permet d'éliminer les oscillations numériques et d'améliorer la précision des résultats.

Après l'implémentation de cette méthode, nous avons travaillé sur des techniques à lui associer afin de réduire le coût des simulations. Des couches parfaitement absorbantes que nous avons optimisées pour les simulations des problèmes de diffraction des ondes électromagnétiques sont introduites dans nos modèles pour tronquer le domaine de calcul afin de simuler uniquement les zones présentant un intérêt, ce qui réduit considérablement la taille des maillages et le temps de calcul [RI.22]. Dans de nombreux cas simulés, l'intégration temporelle se fait grâce à des schémas numériques explicites. Ces approches explicites possèdent une contrainte forte sur le pas de temps, car pour garantir la convergence du schéma numérique, le pas de temps utilisé dans les simulations est celui calculé à partir de la taille du plus petit élément du maillage. Cette contrainte conduit à des temps de calcul très important lorsque la géométrie à traiter présente des petits détails nécessitant l'utilisation d'éléments de maillage très petits. Pour réduire le coût
de calcul dû à cette contrainte, nous avons travaillé sur des schémas numériques dit à pas de temps local comme l'approche "multirate explicite" que nous avons implémentée afin de pouvoir faire des simulations explicites avec un pas de temps supérieur à celui calculé sur le plus petit élément du maillage [RI.12]. Ces outils numériques développés servent à la mise en place des modèles d'applications liés à des problématiques industrielles notamment sur les thématiques de contrôle et de compatibilité électromagnétique des systèmes.

Sur le thème de la compatibilité électromagnétique, les approches numériques mises en place ont servi d'outils de calcul pour les travaux de thèse de Mohamed Boubekeur sur la modélisation des boitiers dédiés à électronique embarquée dans l'automobile ou l'aéronautique [18]. Cette thèse a été réalisée pour répondre à la problématique du dimensionnement des nouveaux boitiers d'équipements dont les parois sont en matériaux composites. En effet, par soucis de gain de poids et de place, ces boitiers ont des parois de plus en plus fines et sont fabriqués avec des matériaux hétérogènes plus légers et ayant de bonnes propriétés de blindage. Pour réduire les coûts de développement de ce type de boitiers, la modélisation numérique est nécessaire. Afin de disposer de modèles numériques peut coûteux en temps de calcul, on doit faire face à la prise en compte de milieux de faible épaisseur qui nécessitent l'utilisation de maillage très fins conduisant à des temps de calcul très important. Dans cette thèse, nous avons proposé un modèle permettant d'étudier le coefficient de blindage de ce type de boitier en évitant de mailler ses parois. Le travail réalisé s'appuie sur la construction d'une condition d'interface qui remplace la paroi et reproduit son comportement en présence d'ondes électromagnétiques [RI.20, RI.11]. Ainsi dans les simulations proposées, on évite de mailler les parois du boitier ce qui permet de réduire considérablement la taille des systèmes et le temps de simulation.

Sur le thème du contrôle électromagnétique de structures, nous avons exploités le calcul des champs électromagnétiques en domaine temporel dans les travaux de thèse d'Antony Manet (Hossein Manesh) pour caractériser la signature des défauts non-francs pouvant apparaître dans les réseaux filaires [67]. Cette thèse a été réalisée dans le cadre d'une ANR sur le diagnostic des réseaux filaires en aéronautique (SODDA). L'une des problématiques liée était de pouvoir détecter et réparer les défauts naissant dans le réseau filaire de télécommunications avant que ceux-ci ne deviennent sévères et mettent en danger la sécurité de l'avion. L'analyse des signaux de réflectométrie est la technique communément utilisée pour la détection des défauts dans un réseau filaire [103]. Mais lorsqu'on est en face d'un réflectogramme, il faut pouvoir identifier la part du signal attribuée au défaut naissant. Compte tenu de la faible amplitude de la réponse d'un défaut naissant, celle-ci est très souvent noyée dans les bruits de mesures ou dans les signaux non souhaités comme ceux des connecteurs. Le travail réalisé a été d'une part de mettre en place une procédure permettant d'identifier la réponse d'un défaut naissant [RI.4], et d'autre
part d'établir et étudier les signatures des défauts non-francs en se basant sur une représen-
tation électrique équivalente dont les paramètres servent à alimenter les modèles de lignes de
transmission utilisés par les industriels pour le diagnostic de réseau filaire [RI.7].

Etude des matériaux supraconducteurs

Le développement de nouveaux dispositifs du génie électrique à conduit les acteurs de ce do-
maine à s'intéresser à de nouveaux matériaux pour remplacer les matériaux classiques devenus
souvent très couteux à produire ou ne permettant pas d'atteindre les niveaux de performances
exigés. C'est notamment ce qui se passe pour le transport d'énergie où l'on souhaite réduire
les pertes dues à l'effet joule, ou encore pour les appareils imagerie médicale où l'on souhaite
disposer des champs magnétiques intenses qui nécessiteraient des volumes très important de
matériaux conducteurs classiques comme le cuivre. Pour ces exemples, la solution est d'intro-
duire les matériaux supraconducteurs car étant fortement diamagnétiques et dissipant très peu
d'énergie [113]. Malheureusement l'introduction de ces matériaux dans les systèmes doit faire
face à de nombreuses difficultés, parmi lesquelles on peut citer, le système cryogénique résultant
du fait que c'est à très basse température que l'état supraconducteur se manifeste, et l'instabilité
de l'état supraconducteur qui traduit une grande sensibilité aux évolutions mécaniques, électro-
magnétiques et thermiques. Les avancées réalisées sur la caractérisation de ces matériaux ont
fait en sorte qu'il existe aujourd'hui de nombreux systèmes utilisant ces matériaux, par exemple
les moteurs [58], les limiteurs de courant [114], les aimants à fort champ [100] ou encore les
câbles de transport d'énergie [86]. Le développement de ces systèmes a été rendu possible grâce
aux nombreux travaux de modélisation réalisés pour proposer des solutions adéquates de dimen-
sionnement de ces systèmes [77, 115, 125]. Ces travaux de modélisation s'appuient sur des lois
de comportement du matériau établies par l'expérience, dont les plus utilisées sont le modèle de
Bean [5, 79] pour les supraconducteurs à basses température critique (LTS) et la loi puissance
pour les supraconducteurs à haute température critique (HTC) [15]. Aujourd'hui, toute l'atten-
tion des acteurs de ce domaine est tournée vers l'utilisation des supraconducteurs de type II à
haute température critique (HTC) car moins exigeant sur l'aspect cryogénique. L'expérimen-
tation sur ces matériaux étant couteuse et difficile à mettre en place, le travail de modélisation
numérique devient crucial. Mais celui-ci reste difficile à mettre en oeuvre à cause du compor-
tement non linéaire du matériau, ouvrant ainsi un axe de recherche sur le développement des
outils de modélisation des supraconducteurs HTC. Mes travaux sur cet axe de recherche ont dé-
buté lors de ma thèse de doctorat qui avait pour objectif de construire des modèles de calcul de
champs induits dans ces matériaux dans le but d'estimer l'impact de la vitesse de variation des
sources d'excitation utilisées lors de la caractérisation expérimentale de ceux-ci [56]. Après la
thèse, j'ai continué à travailler sur le développement d'outils numériques pour la modélisation
tri-dimensionnelle de ces matériaux [RI.6, RI.9, RI.17, RI.23]. Ces outils numériques servent
de base de travail pour aborder des problématiques du rencontrées dans l'industrie. C'est no-
tamment ce qui est fait dans les travaux de thèse de Makong Ludovic dont l'objectif a été de
proposer une approche pour l'évaluation des pertes dans les câbles supraconducteurs à brins
torsadés [66]. Ce type de câbles est utilisé dans les applications de moteurs et aussi pour le
transport d'énergie. Leur modélisation tri-dimensionnelle représente un challenge du point de
vue numérique car le rapport entre leur longueur et leur section est très élevé et les maillages
générés sont très volumineux, ce qui augmente la taille des systèmes à résoudre et conduit à
des temps de calcul très important pouvant atteindre plusieurs dizaines de jours [37]. Pour ré-
duire les coûts de calcul lors de l'évaluation des pertes sur ces câbles, nous avons proposé une
technique de réduction de modèle basée sur des transformations géométriques et les symétries
du système. Le modèle réduit proposé nous permet d'étudier une petite portion du câble afin
d'évaluer les pertes sur l'ensemble du câble, et ceci avec des temps de calcul de quelques heures
[RI.5] [RI.2].

1.5.3 Collaborations scientifiques

Thème compatibilité électromagnétique

— Institut Montéfiore de l'université de Liège en Belgique : équipe Applied Compu-
tational Electromagnetics (Pr Christophe Geuzaine), depuis 2010

La collaboration établie après mon post-doc dans cette équipe de recherche continue.
Elle s'articule autour du développement et l'utilisation des outils numériques pour les
simulations en l'électromagnétisme. Elle a donnée lieu à des publications communes
([RI.18], [RI.19], [RI.20], [RI.21], [RI.22]) sur le thème de la compatibilité électroma-
gnétique.

— Laboratoire ESYCOM de l'université de Marne la Vallée : Equipe Système RF et optiques (Pr Odile Picon), 01/2014 – 06/2017

Cette collaboration est née de la participation commune à l'ANR SODDA (2012–2016).
Elle a consisté au co-encadrement d'une thèse (Hossein Manesh) financée par l'ANR.
Elle s'est articulée autour de la comparaison des outils numériques pour le diagnostic de
réseau filaire en domaine temporel et a donné lieu à une publication en revue internatio-
nale ([RI.7]).

— Institute of Mechanics, Materials and Civil Engineering de l'université de Louvain-La-Neuve en Belgique : Equipe Applied mechanics and mathematics (Dr Jonathan Lambretchts), depuis 2010

Cette collaboration porte sur le maintien et la mise à jour du code de calcul Galerkin
discontinue en domaine temporel via une plateforme GIT. L'implémentation de ce code
de calcul a débuté durant mon post-doc à l'université de Liège en 2009.

— Laboratoire JLL de Sorbonne Université (Dr Bertrand Thierry), depuis 02/2018
C'est une collaboration qui vient de commencer. La première action engagée a été l'en-
cadrement conjoint d'un stage de Master Recherche réalisé par Xavier Adriaens sur
l'utilisation d'une une méthode de décomposition de domaines pour l'évaluation de
l'exposition des personnes aux champs électromagnétiques générés par les systèmes
de recharge inductive.

— CNES, Division des Lanceurs (Mr David Monchaux) , 02/2015-02/2017
Cette collaboration a duré 2ans et s'est faite dans le cadre d'un contrat post-doctoral
(Eva Dieudonné) sur l'évaluation de l'atténuation des ondes radio par les plasmas des
fusées. Elle a donnée lieu à deux publications en revues internationales ([RI.3], [RI.8]).

— Laboratoire LAGA de l'université de Paris 13 (Dr Mohamed Boubekeur), depuis 05/2018
C'est une collaboration naturelle qui se met en place avec Dr Mohamed Boubekeur (mon premier doctorant) qui occupe aujourd'hui un poste d'ingénieur de recherche CNRS au LAGA. Un travail commun sur la mise en place du calcul GPU pour les applications de la compatibilité électromagnétique est en cours d'élaboration.

Matériaux supraconducteurs

— Laboratoire GREEN de l'université de Lorraine (Pr. Jean Levêque, Pr. Bruno Douine, Dr. Smail Mezani, Dr. Kevin Berger), depuis 2009
Depuis ma thèse dans ce laboratoire, une collaboration s'est établie autour du développement et la validation des outils numériques dédiés à la modélisation des supraconducteurs. Elle porte actuellement sur la comparaison des résultats numériques issus de différents codes de calcul, sur le développement des techniques aidant à la caractérisation de matériaux supraconducteurs. Elle a donnée lieu de plusieurs publications dans des revues ([RI.13], [RI.16], [RI.23], [RI.24])

— AML Superconductivity and Magnets (Dr Philippe Masson), depuis 01/2014
L'origine de cette collaboration est la thèse de Makong Ludovic à l'université de Houston. Un encadrement commun de cette thèse (avec Dr Masson) a été réalisé et le doctorant a passé 50% du temps dans chaque entité. Dans le cadre de cette collaboration, j'ai séjourné 1 mois à l'université de Houston. Cette collaboration porte sur le développement d'outils numériques permettant de simplifier la modélisation des câbles supraconducteurs multi-filamentaires torsadés. De nombreuses publications en sont issues ([RI.5], [RI.6], [RI.9], [RI.10])
— Institut Montéfiore de l'université de Liège en Belgique : équipe Applied Computational Electromagnetics (Pr Christophe Geuzaine), depuis 2010

Une collaboration est en place également sur cette thématique. Elle porte sur le développement des outils numériques robustes pour le calcul de courant induit dans les supraconducteurs. Des publications communes en sont issues ([RI.9], [RI.17], [RI.23])

— Institute of Mechanics, Materials and Civil Engineering de l'université de Louvain-La-Neuve en Belgique : Equipe Applied mechanics and mathematics (Dr Jonathan Lambretchts), depuis 2010

La collaboration sur ce thème porte sur l'utilisation de la méthode Galerkin discontinue pour la modélisation des supraconducteurs. Des outils spécifiques (assemblage du système matriciel, solveur spécifique) dédiés aux supraconducteurs ont été développés. Des publications communes ont été réalisées ([RI.10], [RI.17], [RI.23])

### 1.5.4 Activités d'intérêt collectif

1. Expertise de projet ANR (depuis 2015) : 2 projets expertisés

2. Rapporteur d'articles de revues et de conférences

— Revues International (6 articles/an) : *IEEE Electromagnetic Compatibility, IEEE Applied Superconductivity, IEEE Transactions on Magnetics*

— Conférences internationales (8 articles/an) : *Conference on Electromagnetic Field Computation (CEFC), Symposium on Electric and Magnetic Fields (EMF), Conference on the Computation of Electromagnetic Fields (COMPUMAG), Conférence sur les Méthodes Numériques en Electromagnétisme (NUMELEC)*

3. Participation à un Groupe de Travail

— GDR SEEDS : membre du GT Supra autour de la modélisation et la caractérisation des matériaux supraconducteurs, depuis 2017

### 1.5.5 Activités d'encadrement

#### Stages Master Recherche

— 02/2011-06/2011 : Mohamed Boubekeur,

Formation : Master recherche en Mathématiques et Applications de l'université Pierre et Marie Curie

Sujet : *Calcul tridimensionnel de champs induits dans les supraconducteurs par une méthode Galerkin discontinue.*

Taux d'encadrement : 100%
Valorisation scientifique : 1 article de revue [RI.17]

— 03/2013-07/2013 : Smail Ziani,
Master recherche en Physique et Ingénierie de l'Energie de l'ENS CACHAN.

*Modélisation électromagnétique de défauts non-francs en vue du diagnostic de câbles.*

Taux d'encadrement : 100%

Valorisation : 1 article de conférence [CI.18]

— 03/2014-07/2014 : Xiang Liu,
Master recherche en Physique et Ingénierie de l'Energie de l'université de Paris-Saclay.

*Diagnostic filaire, impact des défauts non-francs sur les signaux de réflectometrie.*

Taux d'encadrement : 50%

— 03/2018-07/2014 : Xavier Adriaens,
3ième année cycle ingénieur CentraleSupelec.

*Méthode de décomposition de domaines appliquée à l'optimisation d'un capteur intra-corporel pour la transmission de données médicales.*

Taux d'encadrement : 50%

## Thèses

— **10/2011-12/2014 : Mohamed Boubekeur**

Sujet : *Evaluation de l'efficacité de blindage de structures avec plaques minces : modélisation par une méthode de Galerkin discontinue.*

Financement : Allocation doctorale MESR

Taux d'encadrement : 75%

Date de soutenance : 10 décembre 2014

Production scientifique : 4 articles de revues ([RI.11], [RI.18], [RI.19], [RI.20]), 4 communications en conférences internationales ([CI.16], [CI.26], [CI.31], [CI.34]), 2 communications en conférences nationales ([CN.6], [CN.7]).

Directeur de thèse : Lionel Pichon

Devenir du doctorant : Ingénieur de Recherche au LAGA, université de Paris 13.

— **01/2013-06/2016 : Hossein Manesh**

Sujet : *Etude des défauts non-francs sur des câbles en vue du diagnostic.*

Financement : ANR SODDA

Taux d'encadrement : 33%

Date de soutenance : 21 juin 2016

Production scientifique : 1 article de revue ([RI.7]), 1 communication en confé-
rence internationale ([CI.17]), 2 communications en conférences nationales ([CN.1], [CN.5]).

Directeurs de thèse : Lionel Pichon, Odile Picon

Devenir du doctorant : Ingénieur R&D, WireScan AS, Norvège.

— **01/2014-11/2017 : Ludovic Makong**

Sujet : 3D modeling of twisted multi-filamentary superconductors.

Financement : Bourse université de Houston

Taux d'encadrement : 50%

Date de soutenance : 28 novembre 2017

Production scientifique : 4 articles de revues ([RI.5], [RI.6], [RI.9], [RI.10]), 7 communications en conférences internationales ([CI.4], [CI.5], [CI.14], [CI.20], [CI.22], [CI.23]).

Directeur de thèse : Frédéric Bouillault

Devenir du doctorant : Ingénieur R&D, ORANO-Areva, France

— **10/2014-12/2017 : Xiang Liu**

Sujet : 3D Contrôle non destructif du sol et imagerie d'objets enfouis par des systèmes bi- et multi-statiques : de l'expérience à la modélisation.

Financement : Bourse fondation Supelec

Taux d'encadrement : 33%

Date de soutenance : 13 décembre 2017

Production scientifique : 6 communications en conférences internationales ([CI.1], [CI.2], [CI.9], [CI.10], [CI.11], [CI.12]), 2 communications en conférences nationales ([CN.3], [CN.4]).

Directeurs de thèse : Lionel Pichon, Marc Lambert

Devenir du doctorant : Ingénieur R&D, Valéo, France

### 1.5.6 Publications scientifiques

L'ensemble de mes travaux de recherche a fait l'objet de publications dont la liste détaillée est présentée ci-dessous

— 26 articles dans des revues internationales avec comité de lecture.

— 37 communications à des conférences internationales avec actes et comité de lecture.

— 7 communications à des conférences nationales avec actes et comité de lecture.
Articles dans revues internationales

[RI.1] Adel Razek, Lionel Pichon, Abelin Kameni, Ludovic Makong, and Sahand Rasm. Evaluation of human exposure owing to wireless power transfer systems in electric vehicles. *Athens Journal of Technology and Engineering*, accepted, 2019.

[RI.2] Abelin Kameni, Ludovic Makong, Frederic. Bouillault, and Philippe Masson. Reduced model to compute ac losses of twisted multi-filamentary superconductors. *IEEE Transactions On Applied Superconductivity*, 29 :5901306, 2019.

[RI.3] Eva Dieudonné, Abelin Kameni, David Monchaux, and Lionel Pichon. Radio frequency attenuation by a rocket plume using diffraction theory and finite element modeling. *Acta Astraunica*, 158 :334-341, 2019.

[RI.4] Abelin Kameni, Florent Loete, and Lionel Pichon. Time domain reflectometry model : analysis and characterization of a chafing defect in a coaxial cable. *European Physical Journal : Applied Physics*, 83 :30601, 2018.

[RI.5] Ludovic Makong, Abelin Kameni, Philippe Masson, and Frederic. Bouillault. 3d equivalent model to compute the electro-magnetic behavior of twisted multi-filamentary superconductors wires. *IEEE Transactions On Magnetics*, 54(3) :7202904, 2018.

[RI.6] Ludovic Makong, Abelin Kameni, Philippe Masson, and Frederic. Bouillault. h-formulation using the discontinuous galerkin method for the 3d modeling of superconductors. *IEEE Transactions On Magnetics*, 54(3) :7205204, 2018.

[RI.7] Antony Manet, Abelin Kameni, Florent Loete, Gerome. Genoulaz, Lionel Pichon, and Odile Picon. Equivalent circuit model of soft shield defects in coaxial cables using numerical modeling. *IEEE Transactions on Electromagnetic Compatibility*, 59(2) :533-536, 2017.

[RI.8] Eva Dieudonné, Abelin Kameni, David Monchaux, and Lionel Pichon. Radio frequency attenuation by rocket plume from ground study to in-flight prediction. *Applied Computational Electromagnetics Society*, 32(11) :1048-1055, 2017.

[RI.9] Ludovic Makong, Abelin Kameni, Philippe Masson, Frederic Bouillault, and Christophe Geuzaine. Nodal discontinuous galerkin method for high-temperature superconductors modeling based on the h-formulation. *International Journal for Numerical Methods in Engineering*, 31(2) :2298, 2017.

[RI.10] Ludovic Makong, Abelin Kameni, Jonathan Lambrechts, Philippe Masson, and Frederic Bouillault. 3d modeling of heterogeneous and anisotropic superconducting media. *IEEE Transactions On Magnetics*, 52(3) :7205404, 2016.
[RI.11] Mohamed Boubekeur, Abelin Kameni, and Lionel Pichon. Modeling of thin heterogeneous sheets in the discontinuous galerkin method for 3d transient scattering problems. *European Physical Journal : Applied Physics*, 73 :20901, 2016.

[RI.12] Abelin Kameni, Bruno Seny, and Lionel Pichon. Multirate technique for explicit discontinuous galerkin computations of time domain maxwell equations on complex geometries. *IEEE Transactions On Magnetics*, 52(3) :7204604, 2016.

[RI.13] Bruno Douine, Kevin Berger, Charles-Henri Bonnard, Frederic Sirois, Abelin Kameni, and Jean Levêque. Improved method for determining the n-value of hts bulks. *IEEE Transactions On Applied Superconductivity*, 26(3) :6800704, 2016.

[RI.14] Mohamed Lotfi Khene, Lotfi Alloui, Frederic Bouillault, Abelin Kameni, Souri Mohamed Mimoune, and Mouloud Feliach. 3-d numerical evaluation of trapped magnetic field and temperature of a rectangular gdbacuo bulk magnetized by mmpsc method. *IEEE Transactions On Magnetics*, 51(3) :9000204, 2015.

[RI.15] Lotfi Alloui, Mohamed Lotfi Khene, Abelin Kameni, Frederic Bouillault, Souri Mohamed Mimoune, and Mouloud Feliach. Analysis of trapped magnetic field during magnetization process by using spa method of bulk superconductor. *IEEE Transactions On Magnetics*, 51(3) :9000304, 2015.

[RI.16] Bruno Douine, Charles-Henri Bonnard, Kevin Berger, Frederic Sirois, Abelin Kameni, and Jean Levêque. Determination of $j_c$ and n-value of hts pellets by measurement and simulation of magnetic field penetration. *IEEE Transactions On Applied Superconductivity*, 25(4) :8001008, 2015.

[RI.17] Abelin Kameni, Mohamed Boubekeur, Lotfi Alloui, Frederic Bouillault, Jonathan Lambretchs, and Christophe Geuzaine. A 3-d semi-implicit method for computing the current density in bulk superconductors. *IEEE Transactions On Magnetics*, 50(2) :377-380, 2014.

[RI.18] Mohamed Boubekeur, Abelin Kameni, Laurent Bernard, Axel Modave, and Lionel Pichon. 3d modeling of thin resistive sheets in the discontinuous galerkin method for transient scattering analysis. *IEEE Transactions On Magnetics*, 50(2) :7012104, 2014.

[RI.19] Mohamed Boubekeur, Abelin Kameni, Lionel Pichon, Axel Modave, and Christophe Geuzaine. Analysis of transient scattering problems using a discontinuous galerkin method : application to the shielding effectiveness of enclosures with heterogeneous walls. *Journal of Numerical Modelling : Electronic Networks*, 27 :626-635, 2014.

[RI.20] Mohamed Boubekeur, Abelin Kameni, Axel Modave, Laurent Bernard, and Lionel Pichon. Modeling of weakly conducting thin sheets in the discontinuous galerkin method.
thod for shielding effectiveness evaluation. *Applied Computational Electromagnetics Society*, 28:931–937, 2013.

[RI.21] Abelin Kameni, Axel Modave, Mohamed Boubekeur, Valentin Préault, Lionel Pichon, and Christophe Geuzaine. Evaluation of shielding effectiveness of composite wall with a time domain discontinuous galerkin method. *European Physical Journal : Applied Physics*, 64(2):24508–24511, 2013.

[RI.22] Axel Modave, Abelin Kameni, Jonathan Lambrechts, Lionel Pichon, and Christophe Geuzaine. An optimum pml for scattering problems in the time domain. *European Physical Journal : Applied Physics*, 64(2):24502–24508, 2013.

[RI.23] Abelin Kameni, Frederic Bouillault, Christophe Geuzaine, Jean-François Remacle, Jonathan Lambrechts, and Smail Mezani. Discontinuous galerkin method for computing induced fields in superconductors materials. *IEEE Transactions On Magnetics*, 48(2):591–594, 2012.

[RI.24] Abelin Kameni, Smail Mezani, Jean Levêque, and Bruno Douine. Influence of speed variation of a transverse magnetic field on magnetization of an hts cylinder. *IEEE Transactions On Applied Superconductivity*, 21(4):3434–3441, 2011.

[RI.25] Abelin Kameni, Denis Netter, Bruno Douine, Smail Mezani, and Jean Levêque. Scalling solution and n dependence of the eddy-current distribution in a flat superconductor. *IEEE Transactions On Applied Superconductivity*, 20(4):2248–2254, 2010.

[RI.26] Abelin Kameni, Denis Netter, Bruno Douine, Smail Mezani, and Jean Levêque. A 2d robust fe-fv mixed method to handle strong nonlinearities in superconductors. *IEEE Transactions on Magnetics*, 46(8):3445–3448, 2010.

[RI.27] Abelin Kameni, Bruno Douine, Jean Levêque, Frederic Sirois, and Denis Netter. New hybrid fe-fv method for computing current distribution in 2d superconductors : Application to an hts cylinder in transverse magnetic field. *IEEE Transactions On Applied Superconductivity*, 19(3):2423–2427, 2009.

Communications en conférences internationales

[CI.1] Xiang Liu, Mohamed Serhir, Abelin Kameni, Lionel Pichon, and Marc Lambert. Ground penetrating radar data imaging via kirchhoff migration method. In *Applied Computational Electromagnetics Society*, Florence, Italy, 2017.

[CI.2] Xiang Liu, Mohamed Serhir, Abelin Kameni, Lionel Pichon, and Marc Lambert. Buried targets detection from synthetic and measured b-scan ground penetrating radar data. In *11th European Conference on Antennas and Propagation*, Paris, France, 2017.
[CI.3] Eva Dieudonné, Abelin Kameni, Lionel Pichon, and David Monchaux. Radio frequency attenuation by a rocket plume using diffraction theory and finite element modeling. In *7th European Conference for Aeronautics and Space Sciences*, Milan, Italy, 2017.

[CI.4] Ludovic Makong, Abelin Kameni, Philippe Masson, and Frederic Bouillault. 3d equivalent model to compute the electro-magnetic behavior of twisted multi-filamentary superconductors wires. In *21th International Conference on Computation of Electromagnetic Fields*, Daejeon, South Korea, 2017.

[CI.5] Ludovic Makong, Abelin Kameni, Philippe Masson, and Frederic Bouillault. h-formulation using the discontinuous galerkin method for the 3-d modeling of superconductors. In *21th International Conference on Computation of Electromagnetic Fields*, Daejeon, South Korea, 2017.

[CI.6] Kevin Berger, Guillaume Escamez, Loïc Quéval, Abelin Kameni, Lotfi Alloui, and al. Benchmark on the 3d numerical modeling of a superconducting bulk. In *21th International Conference on Computation of Electromagnetic Fields*, Daejeon, South Korea, 2017.

[CI.7] Kevin Berger, Guillaume Escamez, Loïc Quéval, Abelin Kameni, Lotfi Alloui, and al. Modélisations 3d d'un cube supraconducteur, influence du maillage. In *9ème Conférence Européenne sur les Méthodes Numériques en Electromagnétisme*, Paris, France, 2017.

[CI.8] Abelin Kameni, Florent Loete, and Lionel Pichon. Modèle de réflectométrie en domaine temporel pour l'étude des défauts non-francs. In *9ème Conférence Européenne sur les Méthodes Numériques en Electromagnétisme*, Paris, France, 2017.

[CI.9] Xiang Liu, Mohamed Serhir, Abelin Kameni, Lionel Pichon, and Marc Lambert. Detectability of underground electrical cables junction with a ground penetrating radar : electromagnetic simulation and experimental measurements. In *European Geosciences Union General Assembly*, Vienne, Austria, 2016.

[CI.10] Xiang Liu, Mohamed Serhir, Abelin Kameni, Lionel Pichon, and Marc Lambert. Enhanced linear sampling method for gpr imaging. In *International Workshop on New Computational Methods for Inverse Problems*, Cachan, France, 2016.

[CI.11] Xiang Liu, Mohamed Serhir, Abelin Kameni, Lionel Pichon, and Marc Lambert. 3d dgtd method for antennas and gpr scenarios modeling. In *10th International Symposium on Electric and Magnetic Fields*, Lyon, France, 2016.
[CI.12] Xiang Liu, Mohamed Serhir, Abelin Kameni, Lionel Pichon, and Marc Lambert. Gpr imaging via multi-frequency linear sampling method. In *10th International Symposium on Electric and Magnetic Fields*, Lyon, France, 2016.

[CI.13] Eva Dieudonné, Abelin Kameni, Lionel Pichon, and David Monchaux. Comparison of 2d numerical methods for radio frequency attenuation by a rocket plume. In *10th International Symposium on Electric and Magnetic Fields*, Lyon, France, 2016.

[CI.14] Ludovic Makong, Abelin Kameni, Frederic Bouillault, Philippe Masson, and Christophe Geuzaine. Nodal discontinuous galerkin method for high-temperature superconductors modeling based on the h-formulation. In *10th International Symposium on Electric and Magnetic Fields*, Lyon, France, 2016.

[CI.15] Abelin Kameni and Frederic Bouillault. Influence de la direction du champ magnétique appliqué sur les pertes et l'aimantation des supraconducteurs. In *8ème Conférence Européenne sur les Méthodes Numériques en Electromagnétisme*, Saint-Nazaire, France, 2015.

[CI.16] Mohamed Boubekeur, Abelin Kameni, and Lionel Pichon. Modélisation d'une cavité en plaques hétérogènes avec une méthode de galerkin discontinue pour l'évaluation de son efficacité de blindage. In *8ème Conférence Européenne sur les Méthodes Numériques en Electromagnétisme*, Saint-Nazaire, France, 2015.

[CI.17] Hossein Manesh, Jerôme Genoulaz, Abelin Kameni, Florent Loete, and Lionel Pichon Odile Picon. Experimental analysis and modelling of coaxial transmission lines with soft shield defects. In *IEEE International Symposium on Electromagnetic Compatibility*, Dresden, Germany, 2015.

[CI.18] Abelin Kameni, Florent Loete, Smail Ziani, Khaled Kahalerras, and Lionel Pichon. Time domain modeling of soft faults in wiring system by a nodal discontinuous galerkin method with high-order hexahedral meshes. In *20th International Conference on Computation of Electromagnetic Fields*, Montreal, Canada, 2015.

[CI.19] Abelin Kameni, Bruno Seny, and Lionel Pichon. Multirate technique for explicit discontinuous galerkin computations of time domain maxwell equations on complex geometries. In *20th International Conference on Computation of Electromagnetic Fields*, Montreal, Canada, 2015.

[CI.20] Ludovic Makong, Abelin Kameni, Philippe Masson, Jonathan Lambrechts, and Frederic Bouillault. 3d modeling of heterogeneous and anisotropic superconducting media. In *20th International Conference on Computation of Electromagnetic Fields*, Montreal, Canada, 2015.
[CI.21] Bruno Douine, Kevin Berger, Jean Lévêque, Charles-Henri Bonnard, Frederic Sirois, and Abelin Kameni. Influence of n-value and rise rate on the magnetic field penetration in hts bulks. In *12th European Conference on Applied Superconductivity*, Lyon, France, 2015.

[CI.22] Ludovic Makong, Abelin Kameni, Philippe Masson, and Frederic Bouillault. 3d modeling of superconductors modeling based on e-h formulation. In *12th European Conference on Applied Superconductivity*, Lyon, France, 2015.

[CI.23] Ludovic Makong, Abelin Kameni, Philippe Masson, and Frederic Bouillault. 3d modeling of straight uncoupled multifilamentary superconductors subjected to elliptical field. In *12th European Conference on Applied Superconductivity*, Lyon, France, 2015.

[CI.24] Abelin Kameni, Ludovic Makong, Philippe Masson, and Frederic Bouillault. Influence of direction of the applied magnetic field on losses and magnetizations of superconductors. In *Applied Superconductivity Conference*, Charlotte, United States, 2014.

[CI.25] Abelin Kameni, Bruno Seny, and Lionel Pichon. Multirate technique for explicit discontinuous galerkin computations of time domain maxwell equations on complex geometries. In *6th International Conference on Advanced Computational Methods in Engineering*, Gent, Belgique, 2014.

[CI.26] Mohamed Boubekeur, Abelin Kameni, and Lionel Pichon. Modeling of thin heterogeneous sheets in the discontinuous galerkin method for transient problems. In *16th Bienal Conference on Electromagnetic Field Computation*, Annecy, France, 2014.

[CI.27] Mohamed Boubekeur, Abelin Kameni, Laurent Bernard, Axel Modave, and Lionel Pichon. 3d modeling of thin resistive sheets in the discontinuous galerkin method for transient scattering analysis. In *19th International Conference on Computation of Electromagnetic Fields*, Budapest, Hungary, 2013.

[CI.28] Abelin Kameni, Mohamed Boubekeur, Lotfi Alloui, and Frederic Bouillault. A 3d semi implicit method for computing the current density in bulk superconductors. In *19th International Conference on Computation of Electromagnetic Fields*, Budapest, Hungary, 2013.

[CI.29] Axel Modave, Christophe Geuzaine, Mohamed Boubekeur, Lionel Pichon, and Abelin Kameni. Evaluation of shielding effectiveness in the time domain using a dg method with an efficient pml. In *9th International Symposium on Electric and Magnetic Fields*, Bruges, Belgique, 2013.
[CI.30] Abelin Kameni, Frederic Bouillault, and Christophe Geuzaine. Semi-implicit 3d approach in a discontinuous galerkin method for computing current density in superconductors. In *9th International Symposium on Electric and Magnetic Fields*, Bruges, Belgique, 2013.

[CI.31] Mohamed Boubekeur, Abelin Kameni, Laurent Bernard, and Lionel Pichon. Analysis of transient scattering from nonperfectly conductive shields using a discontinuous galerkin method. In *9th International Symposium on Electric and Magnetic Fields*, Bruges, Belgique, 2013.

[CI.32] Axel Modave, Abelin Kameni, Eric Delhez, Lionel Pichon, and Christophe Geuzaine. An optimum pml for scattering problems in time-domain. In *7ème Conférence Européenne sur les Méthodes Numériques en Electromagnétisme*, Marseille, France, 2012.

[CI.33] Abelin Kameni, Axel Modave, Mohamed Boubekeur, Valentin Préault, Lionel Pichon, and Christophe Geuzaine. Evaluation de l'efficacité de blindage de parois hétérogènes par une méthode de galerkin discontinue en domaine temporel. In *7ème Conférence Européenne sur les Méthodes Numériques en Electromagnétisme*, Marseille, France, 2012.

[CI.34] Mohamed Boubekeur, Abelin Kameni, Axel Modave, Christophe Geuzaine, and Lionel Pichon. Modeling of thin resistive sheets in the discontinuous galerkin method for shielding evaluation. In *15th Biennial IEEE Conference on Electromagnetic Field Computation*, Oita, Japan, 2012.

[CI.35] Abelin Kameni, Jonathan Lambrechts, Jean-François Remacle, Smail Mezani, Frederic Bouillault, and Christophe Geuzaine. Discontinuous galerkin method for computing induced fields in superconducting materials. In *19th International Conference on Computation of Electromagnetic Fields*, Sidney, Australia, 2011.

[CI.36] Jean Levêque, Smail Mezani, and Abelin Kameni. Influence of the shape of the current on ac losses in superconducting tube and coil. In *Superconductivity Centennial Conference 2011 - EUCAS-ISEC-ICMC*, The Hague, Holland, 2011.

## Communications en conférences nationales

[CN.1] Hossein Manesh, Abelin Kameni, Jerôme Genoulaz, Florent Loete, and Lionel Pichon Odile Picon. Modélisation et analyse de défauts non francs dans le blindage des lignes coaxiales en vue du diagnostic. In *18ème Colloque International et Exposition sur la Compatibilité ÉlectroMagnétique*, Rennes, France, 2016.

[CN.2] Eva Dieudonné, Abelin Kameni, Lionel Pichon, and David Monchaux. Evaluation de l'atténuation des communications radiofréquences due aux jets d'échappement
des lanceurs spatiaux. In *18ème Colloque International et Exposition sur la Compatibilité ÉlectroMagnétique*, Rennes, France, 2016.

[CN.3] Xiang Liu, Abelin Kameni, Mohamed Serhir, Lionel Pichon, and Marc Lambert. 3d modeling with discontinuous galerkin time domain method for ground penetration radar application. In *Groupe de Recherche Ondes, Assemblée Générale "Interférences d'Ondes"*, Lyon, France, 2015.

[CN.4] Xiang Liu, Mohamed Serhir, Abelin Kameni, Lionel Pichon, and Marc Lambert. Gpr imaging via multi-frequency linear sampling method with limited aperture data. In *Groupe de Recherche Ondes, Assemblée Générale "Interférences d'Ondes"*, Lyon, France, 2015.

[CN.5] Hossein Manesh, Abelin Kameni, Jerôme Genoulaz, Florent Loete, and Lionel Pichon Odile Picon. Modélisation et analyse de défauts non-francs dans les lignes de transmission pour du diagnostic. In *Journées Nationales Microondes*, Bordeaux, France, 2015.

[CN.6] Mohamed Boubekeur, Abelin Kameni, and Lionel Pichon. Modélisation de parois minces conductrices dans une méthode galerkine discontinue pour l'évaluation de l'efficacité de blindages. In *17ème Colloque International et Exposition sur la Compatibilité ÉlectroMagnétique*, Clermont-Ferrand, France, 2014.

[CN.7] Mohamed Boubekeur, Abelin Kameni, and Lionel Pichon. Modélisation de parois minces conductrices dans une méthode galerkin discontinue pour l'Évaluation du blindage de boîtiers d'Équipements Électroniques. In *Symposium de Génie Electronique*, Cachan, France, 2014.
# Chapitre 2

## Calcul des champs rayonnés en compatibilité électromagnétique

### Sommaire

<table><tr><td>2.1</td><td>Introduction</td><td>29</td></tr><tr><td>2.2</td><td>Formulation du système différentiel</td><td>31</td></tr><tr><td>2.2.1</td><td>Les équations de Maxwell</td><td>31</td></tr><tr><td>2.2.2</td><td>Formulation faible du problème</td><td>31</td></tr><tr><td>2.3</td><td>Flux caractéristiques aux interfaces</td><td>32</td></tr><tr><td>2.3.1</td><td>Invariants de Riemann</td><td>33</td></tr><tr><td>2.3.2</td><td>Flux numériques</td><td>35</td></tr><tr><td>2.3.3</td><td>Conditions aux bords</td><td>37</td></tr><tr><td>2.4</td><td>Condition d'absorption pour les problèmes de diffraction</td><td>38</td></tr><tr><td>2.4.1</td><td>Les parois absorbantes</td><td>39</td></tr><tr><td>2.4.2</td><td>Les couches parfaitement absorbantes</td><td>40</td></tr><tr><td>2.4.3</td><td>Optimisation des couches PML</td><td>41</td></tr><tr><td>2.5</td><td>Schéma numérique pour l'évolution temporelle</td><td>49</td></tr><tr><td>2.5.1</td><td>Les méthodes de Runge-Kutta implicites et explicites</td><td>49</td></tr><tr><td>2.5.2</td><td>Runge-Kutta explicite avec pas de temps local</td><td>51</td></tr><tr><td>2.6</td><td>Modélisation de boitiers de blindage</td><td>57</td></tr><tr><td>2.6.1</td><td>Interface simulant les plaques minces homogènes</td><td>57</td></tr><tr><td>2.6.2</td><td>Modélisation 3D de boitiers à parois composites</td><td>63</td></tr><tr><td>2.7</td><td>Diagnostic de réseau filaire</td><td>70</td></tr><tr><td>2.7.1</td><td>Modèle de réflectométrie pour la caractérisation des défauts non-francs</td><td>71</td></tr><tr><td>2.7.2</td><td>Evaluation des paramètres électriques équivalent</td><td>77</td></tr></table>
### 2.1 Introduction

L'étude des interactions entre les champs électromagnétiques et les structures est un thème que l'on retrouve dans de nombreux domaines applicatifs. Les axes de recherches rencontrés sur cette thématique traduisent les enjeux liés à diverses problématiques notamment celui de la fiabilité des systèmes qui est un aspect à maitriser pour garantir la sécurité des personnes et des biens. La fiabilité d'un système repose sur différents facteurs parmi lesquels on trouve, la maitrise des procédés de fabrication, la qualité des éléments utilisés pour la fabrication notamment en terme de matériau, la protection du système contre toute agression extérieure pouvant nuire à son bon fonctionnement, et aussi la maintenance du système qui passe par la détection et réparation des pannes pouvant entrainer des dysfonctionnements. L'étude de ces deux derniers facteurs repose sur la connaissance des champs électromagnétiques à l'intérieur et à l'extérieur du système.

La question de la protection d'un système dans son environnement de fonctionnement est un sujet central dans le domaine de la compatibilité électromagnétique (CEM). Par exemple un système d'électronique embarquée, il est important de connaitre les niveaux de champs afin de pouvoir, caractériser les phénomènes de couplage aux équipements, respecter les normes environnementales notamment celles liées à la pollution électromagnétique, proposer des solutions adéquates de blindage pour la protection du dispositif. Sur la question de la maintenance des dispositifs, les signaux électromagnétiques émis et reçus sont d'importantes sources d'informations permettant de renseigner sur l'état du système. Leur analyse permet de mettre en évidence la présence de points de dysfonctionnements nécessitant une intervention pour réparation. L'évaluation des champs électromagnétiques rayonnés qui est un sujet central à la modélisation numérique dans ce domaine nécessite de disposer d'outils et méthodes de calculs robustes et efficaces.

De nombreuses méthodes et techniques pour le calcul de champs électromagnétiques en compatibilité et contrôle électromagnétique sont connues [3, 111]. Ces outils numériques permettent de réaliser des études autant dans le domaine fréquentiel que dans le domaine temporel. Dans le domaine fréquentiel, les résultats obtenus rendent compte de l'état du système en régime permanent mais ne fournissent aucune information sur l'évolution transitoire. Pour réaliser une étude sur une bande de fréquence importante, les calculs sont répétés pour des valeurs de fréquence judicieusement choisies et les résultats sont présentés sous forme discrète. Les simulations en domaine temporel s'avèrent interessantes car elles permettent d'observer l'évolution du système du régime transitoire vers le régime permanent. Pour réaliser une étude sur une bande de fréquence, une seule simulation avec des signaux à large bande comme le pulse gaussien suffit. Les résultats dans le domaine fréquentiel sont obtenus grâce à la transformée de Fourier, se présentent sous forme continue, et permettent d'avoir plus de détails sur la bande de fréquence
étudiée.

Parmi les approches numériques utilisées pour le calcul des champs électromagnétiques rayonnés en domaine temporel, c'est la méthode des différences finies qui est la technique la plus répandue. Cette méthode est simple à mettre en oeuvre, rapide et efficace. Mais cette efficacité repose sur l'utilisation des maillages cartésiens réguliers, qui malheureusement sont peu adaptés pour traiter les géométries présentant des courbures [53]. Les approches Galerkin discontinues en domaine temporel ont connu un essor important avec le développement des ordinateurs à architectures multi-coeurs et font partie aujourd'hui des méthodes de plus utilisées pour le calcul de champs électromagnétiques en haute fréquence.

Le principe d'une méthode de Galerkin discontinue est de traiter le problème sur chaque élément du maillage à l'aide d'une approche élément finis. Cette méthode peut alors se décrire comme une méthode élément finis standard dans lequel la grandeur calculée est discontinue. Pour apporter la continuité nécessaire au respect des phénomènes physiques mis en jeu, des termes de flux similaires à ceux utilisés dans les méthodes volumes finis sont appliqués aux bords de chaque élément du maillage. La construction de cette méthode repose sur la formulation faible établie sur chaque élément. Sur cet élément, on définit un espace d'approximation polynomiale telle sorte que les inconnues sont, les valeurs aux noeuds (approche nodale) [51], ou les coefficients caractérisant la projection de la solution sur la base polynomiale (approche modale) [57]. Dans ce type de méthode, les intégrales de frontières sur chaque élément sont remplacées par des expressions de "flux numériques" qui dépendent des valeurs de la solution sur les éléments voisins. Cette approche présente l'avantage d'être bien adaptée aux maillages non-structurés et non conformes, ce qui est un atout pour traiter les problèmes mettant jeu des géométries complexes [50]. Le caractère local de la discrétisation rend plus facile l'utilisation des espaces polynomiaux d'ordre élevé, ce qui permet de réduire les oscillations numériques et d'augmenter la précision des résultats [38, 39]. Le principal inconvénient de cette méthode est l'augmentation du nombre d'inconnues qui résulte du dédoublement des inconnues aux interfaces inter-éléments. L'impact de cet inconvénient est réduit par la mise en place du calcul parallèle dont l'implémentation est facilitée par la forme des matrices générées qui sont des matrices par blocs [12, 61].

Dans ce qui suit, on donne les principaux éléments de la discrétisation des équations de Maxwell en domaine temporel par une méthode de Galerkin discontinue à valeurs nodales.
## 2.2 Formulation du système différentiel

### 2.2.1 Les équations de Maxwell

Les équations de Maxwell sont des relations fondamentales de la physique qui régissent l'électromagnétisme. Ces équations au nombre de quatre, permettent de relier les grandeurs physiques que sont, le champ électrique **E**, le champ magnétique **H**, la densité de courant **J**, l'induction magnétique **B**, l'induction électrique **D** et la densité de charges électrique ρ. Elles s'écrivent

$$ \nabla \cdot \mathbf{D} = \rho \qquad (2.1a) $$

$$ \nabla \cdot \mathbf{B} = \mathbf{0} \qquad (2.1b) $$

$$ \frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = \mathbf{0} \qquad (2.1c) $$

$$ \frac{\partial \mathbf{D}}{\partial t} - \nabla \times \mathbf{H} = \mathbf{J} \qquad (2.1d) $$

La solution du système formé par ces 4 équations ne peut être trouvée que grâce à l'introduction de relations supplémentaires appelées lois de comportement qui permettent de relier l'induction électrique **D** au champ électrique **E**, la densité de courant **J** au champ électrique **E** et l'induction magnétique **B** au champ magnétique **H**. Dans nos travaux, les milieux étudiés sont caractérisés par les lois de comportement suivantes : $\mathbf{D} = \underline{\underline{\epsilon}}\mathbf{E}$ où $\underline{\underline{\epsilon}}$ est le tenseur de permittivité, $\mathbf{B} = \underline{\underline{\mu}}\mathbf{H}$ où $\underline{\underline{\mu}}$ est le tenseur de perméabilité, et $\mathbf{J} = \underline{\underline{\sigma}}\mathbf{E}$ où $\underline{\underline{\sigma}}$ est le tenseur de conductivité.

### 2.2.2 Formulation faible du problème

Pour simplifier cette présentation, on considère les équations de Maxwell dans un milieu non conducteur et isotrope, ce qui correspond à avoir $\mathbf{J} = \mathbf{0}$, $\underline{\underline{\epsilon}} = \underline{\underline{\epsilon}}\mathbf{I}$ et $\underline{\underline{\mu}} = \underline{\underline{\mu}}\mathbf{I}$, où $\mathbf{I}$ est le tenseur identité, $\underline{\underline{\epsilon}}$ est la permittivité électrique, et $\underline{\underline{\mu}}$ est la perméabilité magnétique.

Le système à résoudre pour déterminer les champs électromagnétiques rayonnés est formé par les equations de Maxwell-Faraday et Maxwell-Ampere (2.2). Les inconnues de ce système sont les composantes des champs électrique **E** et magnétique **H**

$$ \left\{ \begin{array}{l} \epsilon \frac{\partial \mathbf{E}}{\partial t} - \nabla \times \mathbf{H} = \mathbf{0} \\ \underline{\underline{\mu}} \frac{\partial \mathbf{H}}{\partial t} + \nabla \times \mathbf{E} = \mathbf{0} \end{array} \right. \qquad (2.2) $$

On rappelle que la méthode de Galerkin Discontinue consiste à résoudre le problème posé sur chaque élément d'un maillage du domaine de calcul. Les échanges entre éléments sont assurés par des termes de flux numériques qui sont construits de manière à respecter les propriétés
physiques, notamment la continuité des champs électromagnétiques dans ce cas.

Considérons un maillage $\mathcal{I} = \bigcup_{i=1}^{N_e} \Omega_e$ de notre domaine de calcul $\Omega$, où $N_e$ est le nombre d'élément du maillage. La discrétisation de type Galerkin discontinue s'opère sur la formulation faible établit sur chaque élément du maillage. Pour l'obtenir, on se sert de la formulation faible sur le domaine global $\Omega$. On décompose ensuite les intégrales sur $\Omega$ comme somme d'intégrales sur les éléments $\Omega_e$ et on applique le théorème du rotationnel. On établit alors que la formulation faible du système (2.2) sur chaque élément $\Omega_e$ s'écrit :

$$ \int_{\Omega_e} \frac{\partial \mathbf{E}}{\partial t} \cdot \phi d\Omega_e - \frac{1}{\epsilon} \int_{\Omega_e} \mathbf{H} \cdot (\nabla \times \phi) d\Omega_e - \frac{1}{\epsilon} \int_{\Gamma_e} (\nabla \times \mathbf{H}) \cdot \phi d\Gamma_e = \mathbf{0} \quad (2.3a) $$

$$ \int_{\Omega_e} \frac{\partial \mathbf{H}}{\partial t} \cdot \phi d\Omega_e + \frac{1}{\mu} \int_{\Omega_e} \mathbf{E} \cdot (\nabla \times \phi) d\Omega_e + \frac{1}{\mu} \int_{\Gamma_e} (\nabla \times \mathbf{E}) \cdot \phi d\Gamma_e = \mathbf{0} \quad (2.3b) $$

où $\phi$ est une fonction test, $\Gamma_e$ est le bord de l'élément $\Omega_e$ et $\nu$ sa normale extérieure. Cette écriture montre que le terme intégral sur le bord $\Gamma_e$ ne dépend que de la solution l'élément $\Omega_e$ considéré. Dans ce cas il n'y a pas de couplages entre les éléments. Pour que l'on puisse avoir interaction entre les différents éléments du maillage, ce terme intégrale est remplacé par un terme équivalent appelé flux numérique qui dépend de la solution sur les éléments voisins de $\Omega_e$. Les expressions de flux numériques généralement mises en place font apparaître les fonctions de saut et de moyenne à l'interface. Soit $\Omega_e$ et $\Omega'_e$ deux éléments voisins ayant une interface commune $\Gamma_e = \Omega_e \cap \Omega'_e$. Les fonctions saut et moyenne à l'interface $\Gamma_e$, notées respectivement $\llbracket \bullet \rrbracket$ et $\{\bullet\}$ sont données par :

$$ \llbracket u \rrbracket_{\Gamma_e}^{\Omega_e} = - \llbracket u \rrbracket_{\Gamma_e}^{\Omega'_e} = \frac{u_{\Omega'_e} - u_{\Omega_e}}{2} \quad (2.4a) $$

$$ \{u\}_{\Gamma_e}^{\Omega_e} = \{u\}_{\Gamma_e}^{\Omega'_e} = \frac{u_{\Omega_e} + u_{\Omega'_e}}{2} \quad (2.4b) $$

## 2.3 Flux caractéristiques aux interfaces

Des expressions des flux numériques sont à construire pour remplacer les intégrales de bord de la formulation faible sur chaque élément (2.3). Pour les établir, on s'appuie sur la méthode des caractéristiques et notamment sur le calcul des invariants de Riemann. Lorsqu'on est en présence d'une équation aux dérivées partielles du premier ordre comme c'est le cas avec le système des équations de Maxwell (2.2), la méthode des caractéristiques consiste à déterminer des lignes le long desquelles le problème à résoudre se ramène à une équation différentielle ordinaire. Si le problème à résoudre est conservatif, alors les inconnues du problème différentiel ordinaire obtenu par la méthode des caractéristiques sont les invariants de Riemann [63, 64].
Ceux-ci ont la particularité d'être constants le long des lignes caractéristiques.

### 2.3.1 Invariants de Riemann

Soit une ligne caractéristique décrite par un vecteur $\alpha$ telle que $\|\alpha\| = 1$. Soit $\tau$ la variation de la solution le long de la ligne caractéristique. Les dérivées spatiales d'une grandeur scalaire $v$, peuvent s'écrire :

$$ \frac{\partial \mathbf{v}}{\partial x} = \frac{\partial \tau}{\partial x} \frac{\partial \mathbf{v}}{\partial \tau}; \quad \frac{\partial \mathbf{v}}{\partial y} = \frac{\partial \tau}{\partial y} \frac{\partial \mathbf{v}}{\partial \tau}; \quad \frac{\partial \mathbf{v}}{\partial z} = \frac{\partial \tau}{\partial z} \frac{\partial \mathbf{v}}{\partial \tau} \qquad (2.5) $$

où $\frac{\partial \mathbf{v}}{\partial \tau}$ représente la dérivée de $v$ le long de la ligne caractéristique de vecteur $\alpha = \nabla \tau$. En appliquant ce changement de variable au système (2.2), on obtient :

$$ \left\{ \begin{array}{ll} \frac{\partial \mathbf{E}}{\partial t} - \epsilon^{-1} \left( \alpha \times \frac{\partial \mathbf{H}}{\partial \tau} \right) = 0 & \\ \\ \frac{\partial \mathbf{H}}{\partial t} + \mu^{-1} \left( \alpha \times \frac{\partial \mathbf{E}}{\partial \tau} \right) = 0 & \end{array} \right. \qquad (2.6) $$

que l'on peut mettre sous la forme d'une équation de type :

$$ \frac{\partial \mathbf{u}}{\partial t} + A \frac{\partial \mathbf{u}}{\partial \tau} = 0 \qquad (2.7) $$

$$ \text{avec, } A = \begin{pmatrix} 0 & 0 & 0 & 0 & \epsilon^{-1}\alpha_z & -\epsilon^{-1}\alpha_y \\ 0 & 0 & 0 & -\epsilon^{-1}\alpha_z & 0 & \epsilon^{-1}\alpha_x \\ 0 & 0 & 0 & \epsilon^{-1}\alpha_y & -\epsilon^{-1}\alpha_x & 0 \\ 0 & -\mu^{-1}\alpha_z & \mu^{-1}\alpha_y & 0 & 0 & 0 \\ \mu^{-1}\alpha_z & 0 & -\mu^{-1}\alpha_x & 0 & 0 & 0 \\ -\mu^{-1}\alpha_y & \mu^{-1}\alpha_x & 0 & 0 & 0 & 0 \end{pmatrix} \quad \text{et } \alpha = \begin{pmatrix} \alpha_x \\ \alpha_y \\ \alpha_z \end{pmatrix} $$

et $\mathbf{u} = (E_x, E_y, E_z, H_x, H_y, H_z)^T$.

La matrice $A$ admet trois valeurs propres doubles, $\lambda_1 = 0$, $\lambda_1 = \frac{1}{\sqrt{\mu\epsilon}}$ et $\lambda_2 = -\frac{1}{\sqrt{\mu\epsilon}}$.

On rappelle que, si $\mathbf{u}_{\lambda_i}$ est un vecteur propre associé à la valeur propre $\lambda_i$ ($i = 1, 2, 3$), alors ce vecteur propre doit satisfaire à la relation suivante :

$$ \frac{\partial \mathbf{u}_{\lambda_i}}{\partial t} + \lambda_i \frac{\partial \mathbf{u}_{\lambda_i}}{\partial \tau} = 0 \qquad (2.8) $$

et les invariants de Riemann que l'on cherche, sont ces vecteurs propres associés aux valeurs propres de la matrice $A$.
Le vecteur propre $\mathbf{u}_{\lambda_1}$ s'obtient en faisant le produit scalaire des équations du système (2.6) par $\alpha$. Cela conduit à :

$$ \left\{ \begin{array}{l} \frac{\partial (\alpha \cdot \mathbf{E})}{\partial t} = 0 \\ \frac{\partial (\alpha \cdot \mathbf{H})}{\partial t} = 0 \end{array} \right. \implies \mathbf{u}_{\lambda_1} = \begin{pmatrix} \alpha \cdot \mathbf{E} \\ \alpha \cdot \mathbf{H} \end{pmatrix} \qquad (2.9) $$

Afin de construire les vecteurs propres $\mathbf{u}_{\lambda_2}$ et $\mathbf{u}_{\lambda_3}$, on suppose que les composantes $\alpha_x$, $\alpha_y$ et $\alpha_z$ sont des constantes, et on fait le produit vectoriel des équations du système (2.6) par $\alpha$ :

$$ \frac{\partial (\alpha \times \mathbf{E})}{\partial t} - \epsilon^{-1} \frac{\partial [\alpha \times (\alpha \times \mathbf{H})]}{\partial \tau} = 0 \qquad (2.10a) $$

$$ \frac{\partial (\alpha \times \mathbf{H})}{\partial t} + \mu^{-1} \frac{\partial [\alpha \times (\alpha \times \mathbf{E})]}{\partial \tau} = 0 \qquad (2.10b) $$

Puisque $\alpha$ est un vecteur unitaire, le produit vectoriel des equations (2.10a) et (2.10b) respectivement par $\sqrt{\frac{\epsilon}{\mu}}\alpha$ et $\sqrt{\frac{\mu}{\epsilon}}\alpha$ conduit à :

$$ \frac{\sqrt{\frac{\epsilon}{\mu}} \partial [\alpha \times (\alpha \times \mathbf{E})]}{\partial t} + \frac{1}{\sqrt{\mu\epsilon}} \frac{\partial (\alpha \times \mathbf{H})}{\partial \tau} = 0 \qquad (2.11a) $$

$$ \frac{\sqrt{\frac{\mu}{\epsilon}} \partial [\alpha \times (\alpha \times \mathbf{H})]}{\partial t} - \frac{1}{\sqrt{\mu\epsilon}} \frac{\partial (\alpha \times \mathbf{E})}{\partial \tau} = 0 \qquad (2.11b) $$

Lorsqu'on additionne (2.10a) et (2.10b) respectivement à (2.11b) et (2.11a), on obtient les deux relations suivantes :

$$ \frac{\partial \left[ \alpha \times \mathbf{E} + \sqrt{\frac{\mu}{\epsilon}} \alpha \times (\alpha \times \mathbf{H}) \right]}{\partial t} - \frac{1}{\sqrt{\mu\epsilon}} \frac{\partial \left[ \alpha \times \mathbf{E} + \sqrt{\frac{\mu}{\epsilon}} \alpha \times (\alpha \times \mathbf{H}) \right]}{\partial \tau} = 0 \qquad (2.12a) $$

$$ \frac{\partial \left[ \alpha \times \mathbf{H} + \sqrt{\frac{\epsilon}{\mu}} \alpha \times (\alpha \times \mathbf{E}) \right]}{\partial t} + \frac{1}{\sqrt{\mu\epsilon}} \frac{\partial \left[ \alpha \times \mathbf{H} + \sqrt{\frac{\epsilon}{\mu}} \alpha \times (\alpha \times \mathbf{E}) \right]}{\partial \tau} = 0 \qquad (2.12b) $$

Et lorsqu'on fait la différence de (2.10a) et (2.10b) respectivement à (2.11b) et (2.11a), on
obtient deux relations supplémentaires données par :

$$ \frac{\partial \left[ \alpha \times \mathbf{E} - \sqrt{\frac{\mu}{\epsilon}} \alpha \times (\alpha \times \mathbf{H}) \right]}{\partial t} + \frac{1}{\sqrt{\mu\epsilon}} \frac{\partial \left[ \alpha \times \mathbf{E} - \sqrt{\frac{\mu}{\epsilon}} \alpha \times (\alpha \times \mathbf{H}) \right]}{\partial \tau} = 0 \quad (2.13a) $$

$$ \frac{\partial \left[ \alpha \times \mathbf{H} - \sqrt{\frac{\epsilon}{\mu}} \alpha \times (\alpha \times \mathbf{E}) \right]}{\partial t} - \frac{1}{\sqrt{\mu\epsilon}} \frac{\partial \left[ \alpha \times \mathbf{H} - \sqrt{\frac{\epsilon}{\mu}} \alpha \times (\alpha \times \mathbf{E}) \right]}{\partial \tau} = 0 \quad (2.13b) $$

On remarque que les relations (2.12a) et (2.13b) donnent le vecteur propre $\mathbf{u}_{\lambda_2}$ associé à la valeur propre $\lambda_2 = -\frac{1}{\sqrt{\mu\epsilon}}$ :

$$ \mathbf{u}_{\lambda_2} = \begin{pmatrix} \alpha \times \mathbf{E} + \sqrt{\frac{\mu}{\epsilon}} \alpha \times (\alpha \times \mathbf{H}) \\ \alpha \times \mathbf{H} - \sqrt{\frac{\epsilon}{\mu}} \alpha \times (\alpha \times \mathbf{E}) \end{pmatrix} \quad (2.14) $$

De même, les relations (2.13a) et (2.12b) donnent le vecteur propre $\mathbf{u}_{\lambda_3}$ associé à la valeur propre $\lambda_3 = \frac{1}{\sqrt{\mu\epsilon}}$ :

$$ \mathbf{u}_{\lambda_3} = \begin{pmatrix} \alpha \times \mathbf{E} - \sqrt{\frac{\mu}{\epsilon}} \alpha \times (\alpha \times \mathbf{H}) \\ \alpha \times \mathbf{H} + \sqrt{\frac{\epsilon}{\mu}} \alpha \times (\alpha \times \mathbf{E}) \end{pmatrix} \quad (2.15) $$

### 2.3.2 Flux numériques

Soit deux éléments du maillage que l'on note par $\Omega_e^-$ pour désigner l'élément courant sur lequel on cherche la solution, et $\Omega_e^+$ l'élément adjacent. Ces deux éléments on une interface commune $\Gamma_e = \Omega_e^- \cap \Omega_e^+$, dont normale extérieure est notée $\nu$. On suppose que les lignes caractéristiques sont dirigées suivant la normale $\nu$. Puisque les vecteurs propres $\mathbf{u}_{\lambda_2}$ et $\mathbf{u}_{\lambda_3}$ sont invariants le long des caractéristiques, on peut supposer que :

1. Les invariants qui traversent l'interface $\Gamma_e$ en se propageant de l'élément courant $\Omega_e^-$ vers l'élément adjacent $\Omega_e^+$ (sens des vitesses positives) prennent leurs valeurs dans l'élément courant, ce qui se traduit par $\mathbf{u}_{\lambda_3}^{\text{num}} = \mathbf{u}_{\lambda_3}^{-}$ :

$$ (\nu \times \mathbf{E})^{\text{num}} - \mathbf{Z}^{-} \nu \times (\nu \times \mathbf{H})^{\text{num}} = \nu \times \mathbf{E}^{-} - \mathbf{Z}^{-} \nu \times (\nu \times \mathbf{H}^{-}) \quad (2.16a) $$

$$ (\nu \times \mathbf{H})^{\text{num}} + \mathbf{Y}^{-} \nu \times (\nu \times \mathbf{E})^{\text{num}} = \nu \times \mathbf{H}^{-} + \mathbf{Y}^{-} \nu \times (\nu \times \mathbf{E}^{-}) \quad (2.16b) $$
ce qui s'écrit encore

$$Y^{-} (\mathbf{v} \times \mathbf{E})^{\text{num}} - \mathbf{v} \times (\mathbf{v} \times \mathbf{H})^{\text{num}} = Y^{-} \mathbf{v} \times \mathbf{E}^{-} - \mathbf{v} \times (\mathbf{v} \times \mathbf{H}^{-}) \quad (2.17a)$$

$$Z^{-} (\mathbf{v} \times \mathbf{H})^{\text{num}} + \mathbf{v} \times (\mathbf{v} \times \mathbf{E})^{\text{num}} = Z^{-} \mathbf{v} \times \mathbf{H}^{-} + \mathbf{v} \times (\mathbf{v} \times \mathbf{E}^{-}) \quad (2.17b)$$

$$\text{avec } Z^{-} = \sqrt{\frac{\mu^{-}}{\epsilon^{-}}} \text{ et } Y^{-} = \frac{1}{Z^{-}}$$

2. Les invariants qui traversent l'interface $\Gamma_e$ en se propageant de l'élément adjacent $\Omega_e^+$ vers l'élément courant $\Omega_e^-$ (sens des vitesses négatives) prennent leurs valeurs dans l'élément adjacent, ce qui se traduit par $\mathbf{u}_{\lambda_2}^{\text{num}} = \mathbf{u}_{\lambda_2}^+$ :

$$(\mathbf{v} \times \mathbf{E})^{\text{num}} + Z^{+} \mathbf{v} \times (\mathbf{v} \times \mathbf{H})^{\text{num}} = \mathbf{v} \times \mathbf{E}^{+} + Z^{+} \mathbf{v} \times (\mathbf{v} \times \mathbf{H}^{+}) \quad (2.18a)$$

$$(\mathbf{v} \times \mathbf{H})^{\text{num}} - Y^{+} \mathbf{v} \times (\mathbf{v} \times \mathbf{E})^{\text{num}} = \mathbf{v} \times \mathbf{H}^{+} - Y^{+} \mathbf{v} \times (\mathbf{v} \times \mathbf{E}^{+}) \quad (2.18b)$$

ce qui s'écrit encore

$$Y^{+} (\mathbf{v} \times \mathbf{E})^{\text{num}} + \mathbf{v} \times (\mathbf{v} \times \mathbf{H})^{\text{num}} = Y^{+} \mathbf{v} \times \mathbf{E}^{+} + \mathbf{v} \times (\mathbf{v} \times \mathbf{H}^{+}) \quad (2.19a)$$

$$Z^{+} (\mathbf{v} \times \mathbf{H})^{\text{num}} - \mathbf{v} \times (\mathbf{v} \times \mathbf{E})^{\text{num}} = Z^{+} \mathbf{v} \times \mathbf{H}^{+} - \mathbf{v} \times (\mathbf{v} \times \mathbf{E}^{+}) \quad (2.19b)$$

$$\text{avec en notant } Z^{+} = \sqrt{\frac{\mu^{+}}{\epsilon^{+}}} \text{ et } Y^{+} = \frac{1}{Z^{+}}.$$

Lorsqu'on additionne (2.17a) à (2.19a) et (2.17b) à (2.19b) on établit les expressions suivantes

$$(\mathbf{v} \times \mathbf{E})^{\text{num}} = \mathbf{v} \times \frac{\{\mathbf{Y}\mathbf{E}\}}{\{\mathbf{Y}\}} + \mathbf{v} \times \frac{(\mathbf{v} \times \llbracket \mathbf{H} \rrbracket)}{\{\mathbf{Y}\}} \quad (2.20a)$$

$$(\mathbf{v} \times \mathbf{H})^{\text{num}} = \mathbf{v} \times \frac{\{\mathbf{Z}\mathbf{H}\}}{\{\mathbf{Z}\}} - \mathbf{v} \times \frac{(\mathbf{v} \times \llbracket \mathbf{E} \rrbracket)}{\{\mathbf{Z}\}} \quad (2.20b)$$

Ces expressions (2.20a) et (2.20b) constituent des flux numériques de type décentré avant, en-
core appelé "flux upwind". Les schemas numériques construits à partir de ce type de flux sont
non conservatifs et dissipatifs. Les flux centrés sont obtenus lorsque la fonction saut à l'interface
est nulle, ce qui conduit aux expressions ci-dessous (2.21a) et (2.21b). Avec ces flux centrés, les
schemas numériques sont conservatifs et dispersifs.

$$(\mathbf{v} \times \mathbf{E})^{\text{num}} = \mathbf{v} \times \frac{\{\mathbf{Y}\mathbf{E}\}}{\{\mathbf{Y}\}} \quad (2.21a)$$

$$(\mathbf{v} \times \mathbf{H})^{\text{num}} = \mathbf{v} \times \frac{\{\mathbf{Z}\mathbf{H}\}}{\{\mathbf{Z}\}} \quad (2.21b)$$
### 2.3.3 Conditions aux bords

Dans l'approche Galerkin discontinue présentée, des termes de flux sont introduits pour décrire les échanges aux interfaces entre éléments. Pour des éléments du maillage possédant une interface sur le bord du domaine, des termes de flux spécifiques doivent être construits pour décrire l'interaction avec l'extérieur. Ceux-ci doivent permettent d'imposer les conditions aux limites physiques souhaitées. Dans nos travaux, pour obtenir les termes de flux correspondant aux conditions aux bords, on adapte ceux mis en place précédemment (2.20) et (2.21), en redéfinissant les fonctions saut et moyenne. Ci-dessous, on présente les expressions de flux utilisées pour imposer un champ électrique ou magnétique dans nos simulations.

#### Dirichlet électrique

C'est une condition qui permet d'imposer un champ électrique $E_{\text{ext}}$ sur toute ou partie du bord. Elle correspond à :

- Imposer le champ électrique extérieur : $\nabla \times E_D = \nabla \times E_{\text{ext}}$.

- Imposer la continuité de la composante tangentielle du champ magnétique : $\nabla \times H_D = \nabla \times H_{\text{ext}}$.

- Utiliser les termes de flux spécifiques suivant :

▷ Les flux décentrés "upwind"

$$ \{\mathbf{u}\} = \begin{pmatrix} \frac{\mathbf{E}_D + \mathbf{E}^-}{2} \\ \mathbf{H}^- \end{pmatrix}, \quad \llbracket \mathbf{u} \rrbracket = \begin{pmatrix} \frac{\mathbf{E}_D - \mathbf{E}^-}{2} \\ 0 \end{pmatrix} \quad \text{et} \quad \{Z\} = Z^+ = Z^- \quad \text{ce qui conduit à :} $$

$$ \left\{ \begin{array}{l} (\nabla \times \mathbf{E})^{\text{num}} = \nabla \times \frac{\mathbf{E}_D + \mathbf{E}^-}{2} \\ (\nabla \times \mathbf{H})^{\text{num}} = \nabla \times \mathbf{H}^- - \nabla \times \left( \nabla \times \frac{\mathbf{E}_D - \mathbf{E}^-}{2Z^-} \right) \end{array} \right. \qquad (2.22) $$

▷ Les flux centrés

$$ \{\mathbf{u}\} = \begin{pmatrix} \mathbf{E}_D \\ \mathbf{H}^- \end{pmatrix} \quad \text{et} \quad \{Z\} = Z^+ = Z^-, \text{ ce qui donne} $$

$$ \left\{ \begin{array}{l} (\nabla \times \mathbf{E})^{\text{num}} = \nabla \times \mathbf{E}_D \\ (\nabla \times \mathbf{H})^{\text{num}} = \nabla \times \mathbf{H}^- \end{array} \right. \qquad (2.23) $$

On note que si $E_{\text{ext}} = 0$, on obtient la condition de Dirichlet homogène électrique qui est connue sous le nom de mur électrique. Elle est couramment utilisée pour traiter les conducteurs parfaits.
### Dirichlet magnétique

C'est une condition qui permet d'imposer un champ électrique $H_{ext}$ sur toute ou partie du bord. Elle correspond à :

- Imposer le champ magnétique extérieur : $\nabla \times H_D = \nabla \times H_{ext}$.

- Imposer la continuité de la composante tangentielle du champ électrique : $\nabla \times E_D = \nabla \times E_{ext}$.

- Utiliser les termes de flux spécifiques suivant :

▷ Les flux décentrés "upwind"

$$ \{\mathbf{u}\} = \begin{pmatrix} \mathbf{E}^{-} \\ \frac{\mathbf{H}_{D} + \mathbf{H}^{-}}{2} \end{pmatrix}, \quad \llbracket \mathbf{u} \rrbracket = \begin{pmatrix} 0 \\ \frac{\mathbf{H}_{D} - \mathbf{H}^{-}}{2} \end{pmatrix} \quad \text{et} \quad \{Y\} = Y^{\text{ext}} = Y^{-}, \text{ ce qui conduit} $$

à :

$$ \left\{ \begin{aligned} (\nabla \times \mathbf{E})^{\text{num}} &= \nabla \times \mathbf{E}^{-} - \nabla \times \left( \nabla \times \frac{\mathbf{H}_D - \mathbf{H}^{-}}{2Y^{-}} \right) \\ (\nabla \times \mathbf{H})^{\text{num}} &= \nabla \times \frac{\mathbf{H}_D + \mathbf{H}^{-}}{2} \end{aligned} \right. \qquad (2.24) $$

▷ Les flux centrés

$$ \{\mathbf{u}\} = \begin{pmatrix} \mathbf{E}^{-} \\ \mathbf{H}_{D} \end{pmatrix}, \text{ ce qui donne} $$

$$ \left\{ \begin{aligned} (\nabla \times \mathbf{E})^{\text{num}} &= \nabla \times \mathbf{E}^{-} \\ (\nabla \times \mathbf{H})^{\text{num}} &= \nabla \times \mathbf{H}_{D} \end{aligned} \right. \qquad (2.25) $$

On note que si $H_{ext} = 0$, on obtient la condition de Dirichlet homogène magnétique qui est connue sous le nom de mur magnétique.

## 2.4 Condition d'absorption pour les problèmes de diffraction

Une onde électromagnétique peut se propager dans un domaine homogène infini sans ren- contrer d'obstacle. Mais, si elle en rencontre un, alors il apparaitra un phénomène de diffraction qui est caractérisé par la création d'une onde diffractée. L'étude et la caractérisation de ces ondes diffractées peut permettre d'avoir des informations sur l'object diffractant. C'est notamment ce principe qui est exploité lors du contrôle électromagnétique des systèmes notamment sur les problématiques de détection d'objets enfouis, ou encore de diagnostic de réseaux filaires par réflectométrie.

La modélisation du phénomène de diffraction des ondes électromagnétiques nécessite une bonne description des ondes diffractées. Lorsque la situation réelle à simuler impose un domaine de
calcul borné, par exemple par des parois de type murs électrique ou murs magnétique, les ondes diffractées interfèrent avec celles réfléchies par les parois, ce qui rend leur étude d'autant plus difficile. Mais lorsque la situation réelle impose de considérer un domaine non borné, par exemple la propagation dans l'espace libre, on doit face à l'impossibilité de construire un tel milieu. Il est alors nécessaire d'établir un problème équivalent défini sur un domaine borné. La question que cela soulève est de savoir comment simuler la continuité de la propagation hors des limites d'un milieu borné, de telle sorte que cela soit équivalent à la propagation dans un milieu infini c'est à dire sans création d'ondes réfléchies aux bords du domaine. Pour faire cela, deux variétés de techniques ont été développées, les parois absorbantes et les couches absorbantes.

### 2.4.1 Les parois absorbantes

On appelle paroi parfaitement absorbante, une frontière fictive sur laquelle il n'y a pas de phénomène de reflection des ondes. Sur cette frontière fictive, on définit une condition absorption qui ne doit pas modifier les données initiales du problème. Ce type de condition aux frontières est défini à l'aide d'un opérateur pseudo-différentiel global qui doit permettre d'absorber les ondes au niveau de la paroi [47]. Dans le cas des équations de Maxwell elle se définit comme une équation du type :

$$ \mathcal{B}(\mathbf{v} \times \mathbf{E}, \mathbf{v} \times \mathbf{H}) = 0 \qquad (2.26) $$

où $\mathcal{B}$ est un opérateur pseudo différentiel s'exprimant en fonction des composantes tangentielles des champs $\mathbf{E}$ et $\mathbf{H}$. En pratique, il est très difficile d'introduire ce genre opérateur dans les schémas numériques et on a alors recours à des approximations. L'approximation au premier ordre de l'équation (2.26), encore connu sous le nom de condition de Silver-Muller, est très répandue dans les travaux de modélisation en électromagnétisme. Celle-ci s'appuie sur l'hypothèse selon laquelle les champs électromagnétiques à l'extérieur du domaine de calcul se propagent comme des ondes planes à incidence normale, ce qui se traduit par : $\mathbf{v} \times \mathbf{E} = \mathbf{Z}\mathbf{H}$ et $\mathbf{v} \times \mathbf{H} = -\mathbf{Y}\mathbf{E}$, avec $\mathbf{Z} = \sqrt{\frac{\mu}{\epsilon}}$ et $\mathbf{Y} = \frac{1}{\mathbf{Z}}$. En effectuant un produit vectoriel ces deux expressions par la normale $\mathbf{v}$, on obtient deux relations qui représentent chacune une condition d'absorption :

$$ \mathbf{v} \times \mathbf{E} - \mathbf{Y}\mathbf{v} \times (\mathbf{v} \times \mathbf{H}) = \mathbf{0} \qquad (2.27a) $$

$$ \mathbf{v} \times \mathbf{H} + \mathbf{Z}\mathbf{v} \times (\mathbf{v} \times \mathbf{E}) = \mathbf{0} \qquad (2.27b) $$

On remarque que ces deux relations (2.27) correspondent à annuler l'invariant de Riemann $u_{\lambda_3}$ (2.15) associé à la valeur propre $\lambda_3 = \frac{1}{\sqrt{\mu\epsilon}}$. Cette condition veut aussi dire qu'aucune information n'est transportée par les caractéristiques entrantes. Du côté pratique, cette condition permet de très bien absorber les ondes à incidence normale. Cette absorption se dégrade lorsque les
ondes sont à incidence oblique et on assiste à l'apparition des faibles réflexions pouvant modifier les résultats escomptés. Néanmoins, cette condition est satisfaisante pour de nombreux modèles. Dans nos travaux, nous l'avons adaptée de manière à pouvoir l'utiliser dans des simulations nécessitant de modéliser un champ incident venant de l'extérieur du domaine. La condition implémentée permet d'imposer un champ entrant dans le domaine de calcul et d'absorber les champs sortant du domaine de calcul.

L'implémentation réalisée correspond à :

- Imposer la continuité de la composante tangentielle des champs électrique et magnétique : $v \times \mathbf{E}_D = v \times \mathbf{E}_{\text{ext}}$ et $v \times \mathbf{H}_D = v \times \mathbf{H}_{\text{ext}}$

- Utiliser les termes de flux spécifiques suivant :

▷ Les flux décentrés "upwind"

$$ \{\mathbf{u}\} = \begin{pmatrix} \frac{\mathbf{E}_D + \mathbf{E}^-}{2} \\ \frac{\mathbf{H}_D + \mathbf{H}^-}{2} \end{pmatrix}, \quad \llbracket \mathbf{u} \rrbracket = \begin{pmatrix} \frac{\mathbf{E}_D - \mathbf{E}^-}{2} \\ \frac{\mathbf{H}_D - \mathbf{H}^-}{2} \end{pmatrix} \quad \text{et} \quad \{Y\} = Y^{\text{ext}} = Y^-, \text{ ce qui conduit} $$

à :

$$ \left\{ \begin{aligned} (v \times \mathbf{E})^{\text{num}} &= v \times \frac{\mathbf{E}_D + \mathbf{E}^-}{2} - v \times \left( v \times \frac{\mathbf{H}_D - \mathbf{H}^-}{2Y^-} \right) \\ (v \times \mathbf{H})^{\text{num}} &= v \times \frac{\mathbf{H}_D + \mathbf{H}^-}{2} + v \times \left( v \times \frac{\mathbf{E}_D - \mathbf{E}^-}{2Z^-} \right) \end{aligned} \right. \quad (2.28) $$

▷ Les flux centrés

$$ \{\mathbf{u}\} = \begin{pmatrix} \mathbf{E}_D \\ \mathbf{H}_D \end{pmatrix}, \text{ ce qui donne} $$

$$ \left\{ \begin{aligned} (v \times \mathbf{E})^{\text{num}} &= v \times \mathbf{E}_D + Zv \times [v \times (\mathbf{H}_D - \mathbf{H}^-)] \\ (v \times \mathbf{H})^{\text{num}} &= v \times \mathbf{H}_D - Yv \times [v \times (\mathbf{E}_D - \mathbf{E}^-)] \end{aligned} \right. \quad (2.29) $$

### 2.4.2 Les couches parfaitement absorbantes

Pour simuler la propagation en espace libre, une autre technique consiste à utiliser un domaine de calcul fini que l'on entoure par des couches dans lesquelles les ondes sortantes seront très fortement atténuées par l'ajout de termes dissipatifs. Ce concept a servi de base au développement des couches absorbantes parfaitement adaptées (Perfectly Matched Layers ou PML) qui ont été introduites par J-P Bérenger en 1994 pour traiter des problèmes de diffraction des ondes électromagnétiques [8]. Ces couches PML, désignent des zones qui absorbent les ondes sortantes sans aucune reflexion à l'interface avec le domaine de calcul. On dit alors qu'une couche est parfaitement adaptée lorsqu'une onde incidente est parfaitement transmise du domaine de calcul à la couche, c'est à dire lorsque le coefficient de réflexion à l'interface entre la PML et le
domaine de calcul est nul :

$$R_{PML/Domaine} = \left| \frac{u_{ref}}{u_{inc}} \right| = 0 \qquad (2.30)$$

Dans les PML, les équations de propagation valides dans le domaine de calcul sont modifiées dans le but d'atténuer les ondes sortantes. Dans les développements de J-P Bérenger portant sur le cas cartésien, le problème à résoudre est obtenu en décomposant les inconnues initiales comme la somme de contributions dans chaque direction de l'espace. Dans les couches PML, on rajoute des termes de dissipation définis à partir des fonctions d'absorption $\sigma_x(x), \sigma_y(y), \sigma_z(z)$ assimilables à des conductivités [8, 9]. En utilisant les solutions en ondes planes, on établit que la condition de nullité du coefficient de reflexion se traduit sur les fonctions d'absorption $\sigma(\sigma_x, \sigma_y, \sigma_z)$ par :

$$\int_{0}^{\delta} \sigma(s) ds = +\infty \qquad (2.31)$$

avec $\delta$ l'épaisseur de la couche dans la direction considérée. On dira qu'une couche est parfaitement absorbante lorsque l'onde qui a été transmisse à la couche PML devient nulle en atteignant sa frontière extérieure. Ceci est garanti lorsque la fonction $\sigma$ est singulière sur cette frontière. L'utilisation de couches PML autre que cartésiennes dans le domaine temporel, a été rendu possible grâce aux développement des couches PML convexes [60, 110]. La mise en place du système d'équations modifié à résoudre suit un principe légèrement différent de celui proposé par J-P Bérenger. Dans ce cas, les équations de Maxwell, avec des termes de dissipation pour la zone PML sont d'abord transposées au régime frequentiel. Dans ce régime, le passage des équations valides dans le domaine de calcul à celles valides dans les couches PML peut s'interpréter comme un changement de métrique de l'espace [109]. Après un passage en coordonnées curviligne des équations en domaine fréquentiel, et grâce à une transformée de Fourier inverse, on établit le nouveau système à résoudre dans le domaine temporel. Celui-ci contient des équations et des variables supplémentaires qui sont introduites pour s'assurer de l'équivalence entre les équations en domaine fréquentiel et les équations en domaine temporel [74, 75].

### 2.4.3 Optimisation des couches PML

Dans la pratique, il est difficile d'implémenter des couches PML permettant d'avoir une absorption parfaite de l'onde sans aucune réflexion. L'efficacité d'une couche PML va alors dépendre de sa capacité à pouvoir absorber efficacement l'onde tout en minimisant la réflexion à l'interface PML/Domaine. Pour optimiser une couche PML, il faut jouer sur les grandeurs que sont, la fonction d'absorption $\sigma$, son épaisseur $\delta$, et sa discrétisation spatiale, dans le but de minimiser l'erreur par rapport à une solution de référence évaluée sur un domaine très grand
non tronqué. Généralement, l'optimisation de PML est réalisée pour chaque cas de simulation, ce qui peut rendre fastidieux le travail de modélisation. L'objectif du travail réalisé a été de déterminer les paramètres optimum de PML à utiliser pour l'ensemble de nos simulations. Pour cela, l'efficacité des PML construites avec des fonctions d'absorption polynomiales et hyperboliques ont été comparées et leurs paramètres ont été optimisés pour des problèmes de diffraction rencontrés en CEM. Ces travaux ont mis en évidence l'existence d'un paramètre indépendant du problème traité et commun aux fonctions d'absorption hyperboliques couramment utilisées dans les PML. Il a été prouvé dans ces travaux que l'une des deux fonctions d'absorption suivante peut être utilisée pour l'ensemble de nos simulations portant sur le calcul de l'efficacité de blindage : $\sigma_h(d) = \frac{c_0}{\delta - d}$, $\sigma_{sh}(d) = \frac{c_0}{\delta - d} - \frac{c_0}{\delta}$, avec $c_0$ la vitesse de propagation et $\delta$ l'épaisseur de la PML. Ces travaux sont présentés dans l'article ci-après **An optimum PML for scattering problems in the time domain** publié dans le journal *European Physical Journal Applied Physics* [RI.22].
# An optimum PML for scattering problems in the time domain*

Axel Modave<fnref n="1" />,ª, Abelin Kameni<fnref n="2" />, Jonathan Lambrechts<fnref n="3" />, Eric Delhez<fnref n="4" />, Lionel Pichon<fnref n="2" />, and Christophe Geuzaine<fnref n="1" />

Received: 5 October 2012 / Received in final form: 30 November 2012 / Accepted: 2 January 2013

Published online: 6 November 2013 – © EDP Sciences 2013

**Abstract.** In electromagnetic compatibility, scattering problems are defined in an infinite spatial domain, while numerical techniques such as finite element methods require a computational domain that is bounded. The perfectly matched layer (PML) is widely used to simulate the truncation of the computational domain. However, its performance depends critically on an absorption function. This function is generally tuned by using case-dependent optimization procedures. In this paper, we will present some efficient functions that overcome any tuning. They will be compared using a realistic scattering benchmark solved with the Discontinuous Galerkin method.

## 1 Introduction

With the increasing use of electrical, electronic and electromagnetic systems, the study of their undesirable interactions and side effects is an important aspect to be considered by the design engineer. In this context, numerical methods are intensively used to solve problems of electromagnetic compatibility (EMC).

Many EMC problems dealing with radiated perturbations are defined on large or infinite spatial domains. When finite element-type numerical methods are used, a major challenge is then to truncate the computational domain without altering too much the accuracy of the solution.

Different techniques have been developed to handle the truncation of the computational domain. The perfectly matched layer (PML) technique, introduced by Bérenger [1], is one of them and is widely used in both frequency and time domains. The idea is to surround the truncated computational domain with a specific layer. Inside it, the physical medium is adapted in such a way that outgoing waves are progressively damped and leave the domain without reflection. To achieve this, the differential equations of the problem are transformed and an absorption function $\sigma(d)$ governing the damping is introduced.

This paper deals with the choice of the function $\sigma(d)$ that is still tricky when numerical methods are used. We present functions that provide good results without any tuning, unlike the commonly used absorption functions that require case-dependent optimization procedures.

A set of PML equations and a nodal Discontinuous Galerkin (DG) method are introduced in Section 2 for a two-dimensional scattering problem in the time domain. Functional forms of $\sigma$ are presented and compared with case-dependent-optimized polynomial profiles in Section 3. Two examples of scattering applications are finally presented in Section 4.

## 2 Mathematical and numerical framework

### 2.1 Scattering problem in unbounded domain

We consider the two-dimensional scattering of waves by an obstacle in an unbounded domain. The transverse magnetic case is treated, where the electric field has only one component in the z-direction, i.e., $e = [0, 0, e_z]$, while the magnetic field is perpendicular to the z-direction, i.e., $h = [h_x, h_y, 0]$.

In order to perform the numerical simulations with a bounded computational domain, the unbounded domain is truncated, and the resulting domain $\Omega$ is surrounded by a PML $\Omega_{pml}$. The domain $\Omega$ is assumed to be convex. Its boundary $\Gamma$ must be sufficiently regular. The PML thickness $\delta$ is constant.

<footnote n="1">¹ Department of Electrical Engineering and Computer Science (Institut Montéfiore), Université de Liège, Grande Traverse 10, 4000 Liège, Belgium</footnote>

<footnote n="2">² Laboratoire de Génie Électrique de Paris, UMR 8507 CNRS, Supelec, Universities Paris VI and Paris XI, 11 rue Joliot Curie, 91192 Gif-sur-Yvette, France</footnote>

<footnote n="3">³ Institute of Mechanics, Materials and Civil Engineering (iMMC), Université Catholique de Louvain-la-Neuve, avenue George Lemaître 4-6, 1348 Louvain-la-Neuve, Belgium</footnote>

<footnote n="4">⁴ Department of Aerospace and Mechanical Engineering, Université de Liège, Grande Traverse 12, 4000 Liège, Belgium</footnote>

<footnote>* Contribution to the Topical Issue “Numelec 2012”, Edited by Adel Razek.</footnote>

<footnote>ª e-mail: a.modave@ulg.ac.be</footnote>
The non-zero components of the fields e and h are governed by Maxwell's equations in $\Omega$. An incident field is prescribed using a total/scattered field approach with transition at $\Gamma$ (see, e.g., [2], Sect. 5.6). The scattered field formulation is then used in $\Omega_{pml}$ while a total field formulation is considered in $\Omega$.

## 2.2 A general PML for convex truncated domain

In the PML $\Omega_{pml}$, the scattered fields are governed by the modified equations

$$ \varepsilon \frac{\partial e}{\partial t} - \nabla \times h = -\varepsilon (\sigma - \sigma_{\kappa}) q - \varepsilon \sigma_{\kappa} e, \quad (1) $$

$$ \mu \frac{\partial h}{\partial t} + \nabla \times e = -\mu (\sigma - \sigma_{\kappa}) t(\mathbf{t} \cdot \mathbf{h}) - \mu \sigma_{\kappa} \mathbf{h}. \quad (2) $$

The additional vector field $q$ introduced in equation (1) has only one component along $e_z$, i.e., $q = [0, 0, q_z]$, and is governed by

$$ \varepsilon \frac{\partial q}{\partial t} - \nabla_n \times h = -\varepsilon \sigma q, \quad (3) $$

with $\nabla_n = [n(n \cdot \nabla)]$. In these equations, $\sigma(d)$ is the absorption function, assumed to be positive in $\Omega_{pml}$, and

$$ \sigma_{\kappa}(\mathbf{x}) = \frac{1}{\kappa^{-1} + d} \int_{0}^{d} \sigma(r) dr, \quad (4) $$

where $d(\mathbf{x}) = \text{dist}(\mathbf{x}, \Gamma)$ is the distance to $\Gamma$ and $\kappa(\mathbf{x})$ is the curvature of $\Gamma$ at the point closest to $\mathbf{x}$. The unit vector fields $n(\mathbf{x})$ and $t(\mathbf{x})$ give the normal and the tangent directions to $\Gamma$, at the same point. The vector field $n$ then gives the stretching direction of the PML.

These equations are built using a procedure similar to the one described by Collino and Monk [3]. Maxwell's equations are firstly written in the frequency domain with curvilinear coordinates. A complex substitution of $d$ is then made (see [3], Sect. 3). Finally, the time-dependent equations (1) and (2) are obtained using an inverse Fourier transform in space and defining the additional equation (3). These equations are equivalent to those given in [3] for the transverse electric mode and a circular truncation boundary.

If $\sigma(d)$ is singular at the outer boundary of the PML, i.e., at $d = \delta$, the solution must satisfy a homogeneous Dirichlet condition on one of the fields for the problem to be well posed. Alternatively, a Neumann condition or even a radiation condition could be implemented to improve the efficiency of the layer. In this paper, a homogeneous Dirichlet condition on $h$ is always prescribed.

## 2.3 Numerical solution

The equations above are solved using the nodal DG finite element method as described below.

The domain $\Omega$ and the PML $\Omega_{pml}$ are decomposed into triangular elements. Over each element $\Omega_e$, each non-zero Cartesian component of each field is approached by a sum of Lagrange basis functions, i.e.,

$$ u(t, \mathbf{x}) \approx u_e^{\text{num}}(t, \mathbf{x}) = \sum_{n} u_{e,n}^{\text{num}}(t) \phi_{e,n}(\mathbf{x}), \quad \forall \mathbf{x} \in \Omega_e, $$

where $u$ is the considered component, $u_e^{\text{num}}$ is its discrete approximation in the element $\Omega_e$, $u_{e,n}^{\text{num}}$ is the value of $u_e^{\text{num}}$ at the node $n$ and $\phi_{e,i}$ is the corresponding Lagrange function that is equal to one at the node $n$ and equal to zero at the others. The global numerical solution is then discontinuous at the cell interfaces [4–6].

The numerical scheme is obtained by considering a weak form of the equations. Multiplying the equations (1)–(3) by test functions $\tilde{e}$, $\tilde{h}$ and $\tilde{q}$, integrating the resulting equations over each cell $\Omega_e$ and using integration by part, one obtains

$$ \int_{\Omega_e} \left[ \varepsilon \frac{\partial e}{\partial t} \cdot \tilde{e} - \mathbf{h} \cdot (\nabla \times \tilde{e}) - S_e \cdot \tilde{e} \right] d\Omega_e - \int_{\Gamma_e} [n_e \times \mathbf{h}]^{\text{num}} \cdot \tilde{e} d\Gamma_e = 0, \quad (5) $$

$$ \int_{\Omega_e} \left[ \mu \frac{\partial \mathbf{h}}{\partial t} \cdot \tilde{\mathbf{h}} + \mathbf{e} \cdot (\nabla \times \tilde{\mathbf{h}}) - S_{\mathbf{h}} \cdot \tilde{\mathbf{h}} \right] d\Omega_e + \int_{\Gamma_e} [n_e \times \mathbf{e}]^{\text{num}} \cdot \tilde{\mathbf{h}} d\Gamma_e = 0, \quad (6) $$

$$ \int_{\Omega_e} \left[ \frac{\partial q}{\partial t} \cdot \tilde{q} - \mathbf{h} \cdot (\nabla_n \times \tilde{q}) - S_q \cdot \tilde{q} \right] d\Omega_e + \int_{\Omega_e} \left[ \frac{1}{\kappa^{-1} + d} (\mathbf{n} \times \mathbf{h}) \cdot \tilde{q} \right] d\Omega_e + \int_{\Gamma_e} [(\mathbf{n}_e \cdot \mathbf{n})(\mathbf{n} \times \mathbf{h})]^{\text{num}} \cdot \tilde{q} d\Gamma_e = 0, \quad (7) $$

where $S_e$, $S_h$ and $S_q$ are the right-hand sides of the equations (1)–(3), $\Gamma_e$ is the boundary of the cell $\Omega_e$ and $n_e$ is the exterior normal to $\Gamma_e$. The fields are discontinuous at the interfaces between adjacent cells. The values at the two sides of a given interface are both considered to compute the boundary terms. For equations (5) and (6), the two following options are used:

1. A centered discretization, i.e.,

$$ [n_e \times \mathbf{h}]^{\text{num}} = \frac{n_e \times \{Zh\}}{\{Z\}}, $$

$$ [n_e \times \mathbf{e}]^{\text{num}} = \frac{n_e \times \{Ye\}}{\{Y\}}. $$

2. An upwind discretization, i.e.,

$$ [n_e \times \mathbf{h}]^{\text{num}} = \frac{n_e \times \{Zh\}}{\{Z\}} - \frac{1}{\{Z\}} n_e \times (n_e \times [e]), $$

$$ [n_e \times \mathbf{e}]^{\text{num}} = \frac{n_e \times \{Ye\}}{\{Y\}} + \frac{1}{\{Y\}} n_e \times (n_e \times [h]), $$
while a Lax-Friedrichs discretization is considered for the last equation (7), i.e.,

$$[(n_e \cdot n)(n \times h)]^{\text{num}} = (n_e \cdot n) \frac{n \times \{Zh\}}{\{Z\}} - c[\mathbf{q}],$$

where c, Y and Z are the phase velocity, the admittance and the impedance of the medium. {$u$} = $(u^+ + u^-)/2$ and $[\mathbf{h}] = (u^+ - u^-)/2$ are the mean and the jump of the quantity $u$ at the interface. The superscript “_” corresponds to the interior value, while “+” is for the exterior value [6].

The spatial integration scheme is obtained by injecting the discrete fields in the equations and using a numerical integration. The time stepping is made with the Crank-Nicholson method.

## 3 Absorption functions

The success of the PML technique is mainly due to the perfectly matching property at the interface between the PML medium and the physical domain in the continuous context. It ensures the perfect transmission of waves from one medium to the other without reflection. Moreover, for a PML with a finite thickness $\delta$, reflections are avoided if $\sigma$ is singular at the outer boundary of the PML (see [7] or [8]). However, this condition does not apply to discrete contexts (e.g., after discretization by finite elements) and the perfect matching is not ensured.

When numerical methods are used, the properties of the PML depend critically on $\sigma$, on the thickness $\delta$ and on the discretization. The continuous properties of the PML can be nearly recovered by taking a large thickness $\delta$ or by using a very accurate discretization with a fine mesh or high-order elements. This approach is however not feasible because it induces a substantial increase of the computational cost.

The function $\sigma(d)$ is the most crucial parameter. Indeed, the performance of the PML decreases dramatically when the damping of the solution induced by large values of $\sigma$ cannot be captured by the discretization grid [9]. Therefore, $\sigma$ must be chosen in such a way as to introduce enough damping of outgoing waves without inducing a too sharp decrease of the fields in the PML. Polynomial profiles of $\sigma$ allow such a progressive damping:

$$\sigma_p(d) = \alpha \left(\frac{d}{\delta}\right)^n, \qquad (8)$$

with $n > 0$ and where $\alpha$ is the value of the absorption function at the outer side of the layer. This profile is widely used (see, e.g., [7, 8, 10, 11]). However, no general rule exists to choose the parameters $\alpha$ and $n$. Therefore, they are usually tuned using expensive and case-dependent optimization procedures. As an alternative, hyperbolic profiles have been proposed by Bermúdez et al. [7]:

$$\sigma_h(d) = \frac{\alpha}{(\delta - d)^n}, \qquad (9)$$

$$\sigma_{sh}(d) = \frac{\alpha}{(\delta - d)^n} - \frac{\alpha}{\delta^n}, \qquad (10)$$

with again two free dimensionless parameters: $\alpha$ and $n$. Fortunately, these two profiles do not require any tuning with $n = 1$. Bermúdez et al. [7] show in a continuous finite element context with a straight Cartesian PML and time-harmonic equations that the solution is indeed close to the optimum when $\alpha$ is equal to the propagation velocity of the medium. We extended this result to the finite difference and the DG contexts with the time-dependent equations in [9, 12].

## 4 Numerical simulations

We show in this section that the results of the previous works [7, 9, 12] are still valid for more realistic scattering benchmarks.

### 4.1 Scattering by a cavity with an aperture

Consider the scattering by the perfectly conducting cavity studied by Ojeda and Pichon [13] in the time-harmonic context. This cavity is rectangular (600 mm × 400 mm) of

<figure><img src="image_1.png" /><figcaption>Fig. 1. Geometry and mesh of the first benchmark.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>Fig. 2. Relative mean error as a function of $\alpha$ with the different absorption functions for the first benchmark. Elements of order 1 are used with centered (a) or upwind numerical fluxes (b). Vertical lines indicate $\alpha = c_0$.</figcaption></figure>

thickness 20 mm with an aperture of 160 mm in the middle of the left side (Fig. 1). The truncated computational domain is rectangular with rounded corners, surrounded by a PML of thickness $\delta = 0.2$ m. A Gaussian pulse is prescribed with the incident field, i.e.,

$$e^i = E^i \sqrt{\frac{\mu_0}{\varepsilon_0}} \exp \left[ - \left( \frac{(x-x_0) - c_0 t}{\ell} \right)^2 \right] \hat{z}, \quad (11)$$

with $x_0 = -1200$ mm, $\ell = 50$ mm, $c_0 = 1/\sqrt{\mu_0\varepsilon_0}$ and $E^i = 0.5\sqrt{\mu_0/\varepsilon_0}$. It describes a pulse coming from the

<figure><img src="image_3.png" /><figcaption>Fig. 3. Relative mean error as a function of $\alpha$ with the different absorption functions for the first benchmark. Elements of order 2 are used with centered (a) or upwind numerical fluxes (b). Vertical lines indicate $\alpha = c_0$.</figcaption></figure>

left side of the domain and propagating to the right. The final time of the simulation is $T = 250\Delta t$, with the time step $\Delta t = 510^{-11}$ s.

The efficiency of the PML is quantified by computing the relative mean error of the simulation, i.e.,

$$\frac{\int_0^T \int_0^T \left( \frac{\varepsilon_0}{2} \|e - e_{\text{ref}}\|^2 + \frac{\mu_0}{2} \|h - h_{\text{ref}}\|^2 \right) d\Omega dt}{\int_0^T \int_0^T \left( \frac{\varepsilon_0}{2} \|e_{\text{ref}}\|^2 + \frac{\mu_0}{2} \|h_{\text{ref}}\|^2 \right) d\Omega dt}, \quad (12)$$
where $e_{\text{ref}}$ and $h_{\text{ref}}$ are reference fields computed with a larger computational domain. Figures 2 and 3 show this relative mean error as a function of $\alpha$ for different absorption functions: parabolic $\sigma_2$, cubic $\sigma_3$, hyperbolic $\sigma_h$ (with $n = 1$) and shifted-hyperbolic $\sigma_{sh}$ (with $n = 1$).

<figure><img src="image_4.png" /><figcaption>Fig. 4. Geometry and mesh of the second benchmark.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5. Relative mean error as a function of $\alpha$ with the different absorption functions for the second benchmark. Elements of order 1 are used with centered (a) or upwind numerical fluxes (b). Vertical lines indicate $\alpha = c_0$.</figcaption></figure>

For each case, the different absorption functions considered in this numerical experiment give equivalent results when they are optimized, i.e., at the minimum of the curves of Figures 2 and 3, except the hyperbolic function $\sigma_h$ that does not perform so well at order 1 with upwind numerical fluxes.

The smallest error with $\sigma_h$ and $\sigma_{sh}$ is always obtained with a value of $\alpha$ close to the propagation velocity $c_0$.

## 4.2 Scattering by an electromagnetic shield

As a second application, consider the scattering of a plane wave by the wall represented in Figure 4. The wall is made of a non-conducting medium and is filled with circular conducting inclusions stored in squared form. The thickness of the wall is 1 mm and the surfacic inclusion rate is $\tau = 20\%$. The permittivity is $5\varepsilon_0$ in the non-conducting medium and $\varepsilon_0$ in the inclusions. The permeability is $\mu_0$ throughout the wall. The conductivity of the inclusions is $10^3$ S/m.

In this problem, we are not interested in the electromagnetic field in the air. Therefore, two PMLs are directly put against the walls. In this simulation, the incident field (11) is prescribed at the interface between the left PML and the wall using a total/scattered field transition, with $x_0 = -5$ mm, $\ell = 1$ mm and $E^i = 0.01\sqrt{\mu_0/\varepsilon_0}$. Reflected waves are damped in the left PML, while the right PML absorbs transmitted waves.

The relative mean error (12), where the integration is carried out over the shield and the air layers, is plotted in Figure 5. The duration of the simulation and the time step are $T = 3000 \Delta t$ and $\Delta t = 4 \times 10^{-14}$ s, respectively.

The curves confirm the conclusions obtained with the previous benchmark. The optimum $\alpha$ is close to $c_0$ for hyperbolic functions just like for the first test case. While an optimum value can also be identified for polynomial profiles, this value differs significantly from the optimum value found for the cavity problem. It differs also from the optimum values identified in previous works.

## 5 Conclusion

A set of PML equations have been presented for two-dimensional scattering EMC problem in the time domain, to simulate a two-dimensional convex truncation boundary. The computational domain could then be reduced to decrease the computational cost.

Concerning the absorption function $\sigma(d)$, the results of this paper are coherent with those of previous works [7,9,12]: hyperbolic and shifted-hyperbolic functions (9) and (10) can be used with $\alpha = c$ and $n = 1$. They are efficient, without any case-dependent costly tuning. By contrast, while polynomial profiles perform equally well, the optimum values of the corresponding coefficients are strongly case dependent and cannot be set on physical grounds only. The hyperbolic and shifted-hyperbolic profiles are therefore potentially much better
suited to large-scale industrial EMC problems, e.g., three-
dimensional simulations of the shielding effectiveness of
complex structures.

This work was supported in part by the Belgian Science Policy (IAP P6/21 and P7), Belgian French Community (ARC 09/14-02) and Walloon Region (WIST3 No. 1017086 "ONELAB"). Jonathan Lambrechts is a Postdoctoral Researcher with the Belgian National Fund for Scientific Research (FNRS).

References

1. J.P. Bérenger, J. Comput. Phys. **114**, 185 (1994)

2. A. Taflove, S.C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 3rd edn. (Artech House, Boston, 2005)

3. F. Collino, P.B. Monk, SIAM J. Sci. Comput. **19**, 2061 (1998)

4. G. Cohen, X. Ferrieres, S. Pernet, J. Comput. Phys. **217**, 340 (2006)

5. L. Fezoui, S. Lanteri, S. Lohreng, S. Piperno, Math. Model. Numer. Anal. **39**, 1149 (2005)

6. J.S. Hesthaven, T. Warburton, J. Comput. Phys. **181**, 186 (2002)

7. A. Bermúdez, L. Hervella-Nieto, A. Prieto, R. Rodríguez, SIAM J. Sci. Comput. **30**, 312 (2007)

8. F. Collino, P.B. Monk, Comput. Methods Appl. Mech. Eng. **164**, 157 (1998)

9. A. Modave, E. Deleersnijder, E.J.M. Delhez, Ocean Dyn. **60**, 65 (2010)

10. J.P. Bérenger, *Perfectly Matched Layer (PML) for Computational Electromagnetics* (Morgan & Claypool, 2007)

11. S. Gedney, in *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, edited by A. Taflove (Artech House, 2005), pp. 273–328

12. A. Modave, E.J.M. Delhez, C. Geuzaine, in *Proceedings of the 10th International Conference on Mathematical and Numerical Aspects of Waves Propagation, Vancouver, Canada, 2011*, pp. 591–594

13. X. Ojeda, L. Pichon, J. Electromagn. Waves Appl. **19**, 1375 (2005)
## 2.5 Schéma numérique pour l'évolution temporelle

La discrétisation temporelle d'une équation s'opère sur un intervalle fini $[0, T]$. Elle s'applie sur une suite croissante de N points équi-répartis $(t_n)_{0 \le n \le N}$, telle que $0 = t_0 < t_1 < \dots < t_{N-1} < t_N = T$. A partir de cette suite, on définit le pas de temps par $\delta t = t_{n+1} - t_n$ et une solution discrete par $U^n = U(t_n) = U(n\delta t)$. La discrétisation temporelle consiste à approcher l'opérateur de dérivation par une formule de dérivée à droite : $\frac{\partial U}{\partial t} = \frac{U^{n+1} - U^n}{\delta t}$.

Le système d'équations issu de la discrétisation spatiale des équations de Maxwell peut se mettre sous la forme d'un système d'équations différentielles ordinaires :

$$ \frac{\partial U}{\partial t} = F(U, t) \qquad (2.32) $$

avec $U = \begin{pmatrix} \mathbf{E} \\ \mathbf{H} \end{pmatrix}$ et $F(U, t) = \begin{pmatrix} F_E(\mathbf{H}) \\ F_H(\mathbf{E}) \end{pmatrix}$ dans le cas d'un flux centré et $F(U, t) = \begin{pmatrix} F_E(\mathbf{E}, \mathbf{H}) \\ F_H(\mathbf{E}, \mathbf{H}) \end{pmatrix}$ dans le cas d'un flux décentré upwind. Dans mes travaux, pour résoudre ce système d'équations ordinaires, j'ai principalement travaillé avec les méthodes de Runge-Kutta. Dans les sections suivantes, je propose une brève présentation des méthodes Runge-Kutta implicites et explicites, et un résumé de ma contribution à la réduction du temps de simulation grâce à l'hybridation d'une approche Runge-Kutta explicite.

### 2.5.1 Les méthodes de Runge-Kutta implicites et explicites

La méthode de Runge-Kutta sert à calculer itérativement la solution discrète $U^{n+1}$ à partir de la connaissance de la solution $U^n$ en utilisant un ensemble de points intermédiaires définis sur l'intervalle $[t_n, t_{n+1}]$. Les schémas numériques associés à cette approche peuvent être explicites ou implicites.

- Dans les cas explicites, la stabilité et la convergence des schémas numériques dépendent de la condition (Courant - Friedrichs - Levy) imposée sur le pas de temps. On rappelle que cette condition contraint le pas de temps $\delta t$ afin de garantir la convergence. Celui-ci doit satisfaire à l'inégalité suivante :

$$ \max_{\Omega_e \in \mathcal{I}} \left(\frac{c}{h}\right) \delta t \le A \qquad (2.33) $$

où h désigne le rayon de courbure des éléments du maillage, A une constante et c la vitesse de propagation de l'onde.

- Dans les cas implicites, les schémas numériques sont inconditionnellement stables sans
aucune restriction sur le pas de temps, mais nécessitent la résolution de systèmes linéaires ce qui est un inconvénient en terme de temps de calcul à cause des opérations d'inversions matricielles.

Une approche de Runge-Kutta à s étapes intermédiaires repose sur la connaissance des coefficients d'une matrice $A \in M_{s \times s}$ et de deux vecteurs b et c de $\mathbb{R}^s$. Ceux-ci sont généralement présentés sous la forme d'un tableau appelé tableau du Butcher.

$$ \begin{array}{c|c} c & A \\ \hline & b^1 \end{array} $$

<caption>TABLE 2.1 – Tableau de butcher</caption>

### 1. Runge-Kutta explicite

L'algorithme explicite consiste à calculer $U^{n+1}$ à partir de $U^n$ en effectuant les opérations suivantes :

(a) Solutions intermédiaires :

$$ U^k = U^n + \delta t \sum_{j=1}^{k-1} A_{kj} K_j \text{ avec } 1 \le k \le s $$

$$ K_k = F(U^k, t_n + c_k \delta t) $$

(b) Solution à $t_{n+1}$ :

$$ U^{n+1} = U^n + \delta t \sum_{k=1}^{s} b_k K_k $$

Les schémas Runge-Kutta explicites couramment rencontrés sont ceux à 1, 2, 3, 4 étapes intermédiaires dont la matrice $A$ et les vecteurs $b$ et $c$ sont respectivement donnés par :

$$ — s = 1 \text{ (Euler explicite) : } A = (0), b = (0) \text{ et } c = (1) $$

$$ — s = 2 : A = \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix}, b = \begin{pmatrix} 1/2 \\ 1/2 \end{pmatrix} \text{ et } c = \begin{pmatrix} 0 \\ 1 \end{pmatrix} $$

$$ — s = 3 : A = \begin{pmatrix} 0 & 0 & 0 \\ 1/2 & 0 & 0 \\ -1 & 2 & 0 \end{pmatrix}, b = \begin{pmatrix} 1/6 \\ 2/3 \\ 1/6 \end{pmatrix} \text{ et } c = \begin{pmatrix} 0 \\ 1/2 \\ 1 \end{pmatrix} $$

$$ — s = 4 : A = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 1/2 & 0 & 0 & 0 \\ 0 & 1/2 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix}, b = \begin{pmatrix} 1/6 \\ 1/3 \\ 1/3 \\ 1/6 \end{pmatrix} \text{ et } c = \begin{pmatrix} 0 \\ 1/2 \\ 1/2 \\ 1 \end{pmatrix} $$

### 2. Runge-Kutta implicite

L'Algorithme implicite à s étapes est quasiment identique à celui du cas explicite :

(a) Solutions intermédiaires :
$$U^k = U^n + \delta t \sum_{j=1}^{k} A_{kj} K_j \text{ avec } 1 \le k \le s$$

$$K_k = F(U^k, t_n + c_k \delta t)$$

(b) Solution à $t_{n+1}$ :

$$U^{n+1} = U^n + \delta t \sum_{k=1}^{s} b_k K_k$$

La différence se trouve dans le calcul des solutions intermédiaires $U^k$. Puisque la grandeur $K_k$ depend de $U^k$, il est alors nécessaire de résoudre un système matriciel à chaque étapes intermédiaires. Les schémas Runge-Kutta implicites couramment rencontrés sont ceux à 1, 2 étapes intermédiaires dont la matrice $A$ et les vecteurs $b$ et $c$ sont donnés par :

— $s = 1$ étape intermédiaire (Euler implicite) : $A = (1), b = (1)$ et $c = (1)$

— $s = 2$ étapes intermédiaires : $A = \begin{pmatrix} 1 - 1/\sqrt{2} & 0 \\ 1/\sqrt{2} & 1 - 1/\sqrt{2} \end{pmatrix}$, $b = \begin{pmatrix} 1/\sqrt{2} \\ 1 - 1/\sqrt{2} \end{pmatrix}$ et

$$c = \begin{pmatrix} 1 - 1/\sqrt{2} \\ 1 \end{pmatrix}$$

En pratique, pour de grands systèmes, les schémas explicites sont préférés aux schémas implicites car ne nécessitent pas la résolution de problèmes matriciels aux étapes intermédiaires. Malheureusement la convergence de ces approches explicites est liée au pas de temps donné par la condition CFL. Généralement celui utilisé est calculé sur le plus petit élément du maillage, ce qui conduit à des coûts de calcul très importants notamment lorsque le système modélisé contient de petits détails nécessitant des maillages fins.

### 2.5.2 Runge-Kutta explicite avec pas de temps local

Dans de nombreuses applications, il est nécessaire de raffiner localement le maillage pour pouvoir capter certain phénomènes. Dans ces cas, l'utilisation du pas de temps calculé sur le plus petit élément du maillage conduit à des coûts de calcul très importants. Afin de réduire ces coûts, il faut pouvoir utiliser un pas de temps supérieur que celui calculé par le plus petit élément.

Des approches dites à pas de temps local font parties des techniques développées pour modifier localement le pas de temps dans un schéma numérique [35, 59]. Dans nos travaux nous nous sommes intéressés à une approche dite "multirate" Runge-Kutta explicite qui a été initialement mise en place pour résoudre les équations de la dynamique des fluides discrétisées par une approche de Galerkin discontinue [96, 99]. Notre travail a consisté à transposer cette approche à la résolution des équations de Maxwell pour le calcul des champs électromagnétiques rayonnés. L'approche "multirate" Runge-Kutta repose sur le choix d'une méthode Runge-Kutta de base,
et sur le tri des éléments du maillage en fonction de leur taille avec constitution de plusieurs groupes d'éléments. Pour chaque groupe, la méthode Runge-Kutta explicite de base est stable pour le pas de temps calculé sur son plus petit élément. D'une manière générale, cette méthode consiste à résoudre le problème algébrique issu de la discrétisation spatiale des équations, sur différents groupes d'éléments, en utilisant des pas de temps différents, et en s'assurant que les échanges d'informations entre les différents groupes se fassent aux mêmes moments. Pour mettre cela en place, il faut modifier pour chaque groupe les étapes intermédiaires de l'approche Runge-Kutta explicite de base, et introduire des zones dite tampon qui serviront à l'échange d'informations entre les groupes. Dans l'article ci-après **Multirate Technique for Explicit Discontinuous Galerkin Computations of Time-Domain Maxwell Equations on Complex Geometries** publié dans *IEEE Transactions On Magnetics* nous présentons une approche "Multirate" basée sur la méthode Runge-Kutta explicite à deux étapes intermédiaires [RI.12]. Dans cet article, on montre notamment que pour des simulations sur des géométries complexes, on peut réduire considérablement le temps de calcul par rapport à l'approche Runge-Kutta de base utilisant le pas de temps calculé sur le plus petit élément du maillage.
# Multirate Technique for Explicit Discontinuous Galerkin Computations of Time-Domain Maxwell Equations on Complex Geometries

A. Kameni,¹ B. Seny,² and L. Pichon¹

¹Group of Electrical Engineering-Paris, CentraleSupelec, Université Paris Sud, Université Pierre et Marie Curie, Gif-sur-Yvette 91192, France

²Institute of Mechanics, Materials and Civil Engineering, Université Catholique de Louvain, Louvain-la-Neuve B-1348, Belgium

This paper presents a multirate technique to improve and optimize the time step of a second-order two-stage explicit Runge-Kutta (ERK) scheme. In this technique, the mesh elements are stored in different groups according to their stable time steps. These groups are sorted into two classes: 1) the bulk groups, where the two-stage ERK method is applied once or repeatedly; 2) the buffer groups that are used to accommodate the transition between two bulk groups. This technique is proposed for accelerating the explicit discontinuous Galerkin computations of the time-domain Maxwell equations. An application example on a human skull is proposed to show the efficiency of this technique and to simulate wave propagation on complex geometries.

**Index Terms**—Discontinuous Galerkin method, explicit Runge-Kutta (ERK) schemes, Maxwell equations, multirate technique.

## I. INTRODUCTION

EXPLICIT methods, such as the explicit Runge-Kutta (ERK) schemes, are commonly used to perform the electromagnetic modeling of complex systems. The maximum allowable time step to ensure stability and convergence depends on the Courant-Friedrichs-Lewy (CFL) condition. The time step must be proportional to the element size and inversely proportional to the maximum eigenvalue of the Jacobian of the system. Unstructured meshes are often required to capture a wide spectrum of physical scales. Local mesh refinements may be responsible for the significant gaps between the smallest and average element sizes. In such cases, the efficiency of the traditional single-rate explicit time-stepping methods may be drastically low with respect to the problem size.

The discontinuous Galerkin methods are well suited to be used in combination with a local time-stepping strategy to reduce the expensive computations by adapting the time step under local stability conditions [1], [2]. The multirate methods are a subset of the local time-stepping schemes family. They allow the use of different time steps that are integer ratios of each other to solve a discrete system [3]. This technique has been used with a multistep Adams-Bashforth scheme to improve the discontinuous Galerkin computations of the Maxwell equations when the mesh is sorted in two groups [4]. The first group is constituted by fine elements, where the Adams-Bashforth scheme is applied repeatedly, and the second group is constituted by coarse elements, where the Adams-Bashforth scheme is applied once.

The development and the implementation of the ERK multirate schemes have been recently proposed to accelerate geophysical flow computations when the mesh is sorted in few groups [5]. An extension of these strategies to the parallel framework has also been proposed in [6].

In this paper, a multirate technique based on a second-order accurate two-stage ERK method (Heun's method) is presented and applied on the Maxwell equations. This approach consists of gathering the mesh elements in different groups that satisfy the CFL condition for a certain range of time steps. These groups are sorted into two classes: 1) the bulk groups, where a two-stage ERK method is applied once or repeatedly; 2) the buffer groups that are used to accommodate the transition between two bulk groups, where a two-stage ERK scheme is adapted to coincide with the ERK schemes used in the bulk groups. The use of these buffer groups allows an increase in the reduction in the computational cost. The efficiency of this approach is shown through the computations of an electric field on a human skull.

## II. MAXWELL EQUATIONS AND DISCRETE SYSTEM

Let $E$, $H$, and $J$ be the electric field, the magnetic field, and the current density, respectively. They satisfy the Maxwell equations given by

$$ \begin{cases} \epsilon \partial_t E - \nabla \times H = -J \\ \mu \partial_t H + \nabla \times E = 0 \end{cases} \quad (1) $$

where $\epsilon$ and $\mu$ are the permittivity and the permeability of the medium, respectively. In a conductive medium, $J = \sigma E$, with $\sigma$ being the conductivity.

The discontinuous Galerkin methods are introduced for solving the conservative form of partial differential equations. This method consists of discretizing the variational formula-

Manuscript received July 5, 2015; revised September 4, 2015; accepted September 9, 2015. Date of publication September 22, 2015; date of current version February 17, 2016. Corresponding author: A. Kameni (e-mail: abelin.kameni@lgep.supelec.fr).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TMAG.2015.2480537
<table><caption>TABLE I
(LEFT) BUTCHER TABLEAU OF RK22. (CENTER) BUTCHER TABLEAU OF RK22<sup>4a</sup> ISSUED FROM APPLYING TWICE RK22. (RIGHT) BUTCHER TABLEAU OF THE EXTENDED METHOD RK22<sup>4b</sup></caption><thead><tr><th>0</th><th>1</th><th>1/2</th><th>1/2</th><th>0</th><th>1/2</th><th>1/4</th><th>1/4</th><th>1/2</th><th>1/4</th><th>1/4</th><th>1/4</th><th>1/4</th><th>1/4</th><th>1/4</th><th>1/4</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1/2</td><td>1/2</td><td>1</td><td>1/4</td><td>1/4</td><td>1/4</td><td>1/2</td><td>1/4</td><td>1/4</td><td>1/4</td><td>1/4</td><td>1/4</td><td>1/4</td><td>1/4</td></tr></tbody></table>

tion of (2) on each mesh element T of the domain $\Omega = \cup T$

$$ \left\{ \begin{aligned} & \int_T \epsilon \partial_t E \phi - \int_T H \times \nabla \phi - \int_{\partial T} (n \times H) \phi = - \int_T \sigma E \phi \\ & \int_T \mu \partial_t H \psi + \int_T E \times \nabla \psi + \int_{\partial T} (n \times E) \psi = 0 \end{aligned} \right. \quad (2) $$

where $\phi$ and $\psi$ are the test functions.

In each T, a finite-element method is applied, and the mapping technique is used to facilitate the use of high-order mesh elements. The interface terms are replaced by numerical flux expressions as in a finite-volume method. Different formulations of the flux expressions exist [7]. Expression (3), shown below, resulting in different numerical schemes, is implemented. For $\alpha = 0$, the centered fluxes are obtained, and the numerical schemes are dispersive. For $\alpha = 1$, the upwind fluxes are obtained, and the numerical schemes are dissipative

$$ \left\{ \begin{aligned} (n \times H)^{\text{num}} &= n \times \frac{\left\{ \sqrt{\frac{\mu}{\epsilon}} H \right\}}{\left\{ \sqrt{\frac{\mu}{\epsilon}} \right\}} - \alpha \left( n \times \frac{(n \times [E])}{\left\{ \sqrt{\frac{\mu}{\epsilon}} \right\}} \right) \\ (n \times E)^{\text{num}} &= n \times \frac{\left\{ \sqrt{\frac{\epsilon}{\mu}} E \right\}}{\left\{ \sqrt{\frac{\epsilon}{\mu}} \right\}} + \alpha \left( n \times \frac{(n \times [H])}{\left\{ \sqrt{\frac{\epsilon}{\mu}} \right\}} \right) \end{aligned} \right. \quad (3) $$

where $[u] = (u^+ - u^-/2)$ and $\{u\} = (u^+ + u^-/2)$, with $u = (E, H)$. The superscript “-” denotes the values for the fields in the current element, while the superscript “+” is for the adjacent element.

The spatial discretization of Maxwell's equations leads to an ordinary differential equation on each mesh element T that can be presented as a Cauchy problem

$$ \left\{ \begin{aligned} \frac{du_i}{dt} &= f(u_i(t), t) \\ u_i(0) &= u_i^0. \end{aligned} \right. \quad (4) $$

The explicit Runge-Kutta methods are set to integrate the solution in time.

### III. MULTIRATE STRATEGY

#### A. Multirate Time Integration

In this paper, the multirate time integration is based on Heun's method, which is a second-order accurate two-stage ERK method, that we note RK22, whose Butcher tableau is given on the left side of Table I. It allows the computation of

<figure><img src="image_6.png" /><figcaption>Fig. 1. Example group partitioning for multirate. (Left) Initial bulk groups P with time step $2\Delta t_m$ and R with time step $\Delta t_m$. (Right) Buffer group G with time step $2\Delta t_m$ is inserted between the two bulk groups.</figcaption></figure>

the next step solution $u^{n+1}$ from the previous one $u^n$

$$ (s1): u^{(1)} = u^n \quad \text{and} \quad K^1 = f(u^{(1)}, t^n) $$

$$ (s2): u^{(2)} = u^n + \Delta t K^1 \quad \text{and} \quad K^2 = f(u^{(2)}, t^n + \Delta t) $$

$$ u^{n+1} = u^n + \frac{1}{2} \Delta t (K^1 + K^2). \quad (5) $$

The aim of the multirate time integration is to compute the discrete solution in groups whose stable time steps are different. Let us consider a mesh example shown on the left side of Fig. 1. It has two bulk groups of stable time steps: 1) $\Delta t_m$ and 2) $2\Delta t_m$, noted R and P. With a traditional single-rate method, the time step $\Delta t_m$ is used for all the mesh elements. Our goal is to run the elements of group R with the time step $\Delta t_m$ and those of P with $2\Delta t_m$ to reduce the computational cost. The idea is to apply RK22 on P once, and successively on R twice.

The use of RK22 in group R with the time step $\Delta t_m$ has an influence on the integration scheme in group P. This influence is limited at the two connected neighboring elements, because the base method has two stages. If RK22 with $2\Delta t_m$ is used on these elements, the conservation of the fluxes at the interfaces is not satisfied. That is why a buffer region is introduced between the groups of different stable time steps, as shown on the right side of Fig. 1, where a buffer group, noted G, is inserted between the bulk groups P and R.

Let us denote $u_0$ as the solution to an element of R; $u_1$ the solution on an element of G and, neighbor of R; $u_2$ the solution to an element of G and neighbor of P; and $u_3$ the solution to an element of P.

Applying RK22 with the time step $\Delta t_m$ leads to $u_0^{n+(1/2)}$

$$ (s1): u_0^{(1)} = u_0^n; \quad K_0^1 = f(u_0^{(1)}, t^n) $$

$$ \begin{aligned} (s2): u_0^{(2)} &= u_0^n + \frac{1}{2} \Delta t K_0^1; \quad K_0^2 = f \left( u_0^{(2)}, t^n + \frac{1}{2} \Delta t \right) \\ u_0^{n+\frac{1}{2}} &= u_0^n + \frac{1}{4} \Delta t (K_0^1 + K_0^2). \end{aligned} \quad (6) $$

Applying again RK22 allows the computation of $u_0^{n+1}$

$$ (s3): u_0^{(3)} = u_0^{n+\frac{1}{2}}; \quad K_0^3 = f(u_0^{(3)}, t^{n+\frac{1}{2}}) $$

$$ (s4): u_0^{(4)} = u_0^{n+\frac{1}{2}} + \frac{1}{2} \Delta t K_0^3; \quad K_0^4 = f \left( u_0^{(4)}, t^{n+\frac{1}{2}} + \frac{1}{2} \Delta t \right) $$

$$ u_0^{n+1} = u_0^n + \Delta t \left( \frac{1}{4} K_0^1 + \frac{1}{4} K_0^2 + \frac{1}{4} K_0^3 + \frac{1}{4} K_0^4 \right). \quad (7) $$
These four steps constitute a four-stage ERK method, noted RK22<sup>4a</sup>, which is stable for the time step $\Delta t$, whose Butcher tableau is given in Table I (center).

Since RK22<sup>4a</sup> has four stages, the base method RK22 is extended to a four-stage method, noted RK22<sup>4b</sup>, which is stable for the time step $\Delta t$, whose Butcher tableau is given on the right side of Table I. For the elements of $P$ whose integration schemes are not influenced by the use of RK22<sup>4a</sup>, the solution $u_3^{n+1}$ is computed by RK22<sup>4b</sup> in four stages

$$ (s1): u_3^{(1)} = u_3^n; \quad K_3^1 = f(u_3^{(1)}, t^n) $$

$$ (s2): u_3^{(2)} = u_3^n + \Delta t K_3^1; \quad K_3^2 = f(u_3^{(2)}, t^n + \Delta t) $$

$$ (s3): u_3^{(3)} = u_3^{(1)}; \quad K_3^3 = K_3^1 $$

$$ (s4): u_3^{(4)} = u_3^{(2)}; \quad K_3^4 = K_3^2 $$

$$ u_3^{n+1} = u_3^n + \Delta t \left( \frac{1}{4} K_3^1 + \frac{1}{4} K_3^2 + \frac{1}{4} K_3^3 + \frac{1}{4} K_3^4 \right). \quad (8) $$

The solutions in the buffer group $G$ are obtained using the coefficients of the Butcher tableau [see the right side of Table I].

The solution $u_1^{n+1}$ is obtained by

$$ (s1): u_1^{(1)} = u_1^n; \quad K_1^1 = f(u_1^{(1)}, t^n) $$

$$ (s2): u_1^{(2)} = u_1^n + \Delta t K_1^1; \quad K_1^2 = f(u_1^{(2)}, t^n + \Delta t) $$

$$ (s3): u_1^{(3)} = u_1^{(1)}; \quad K_1^3 \neq K_1^1 $$

$$ (s4): u_1^{(4)} = u_1^{(3)} + \Delta t K_1^3; \quad K_1^4 = f(u_1^{(4)}, t^n + \Delta t) $$

$$ u_1^{n+1} = u_1^n + \Delta t \left( \frac{1}{4} K_1^1 + \frac{1}{4} K_1^2 + \frac{1}{4} K_1^3 + \frac{1}{4} K_1^4 \right). \quad (9) $$

Note that $K_1^3 \neq K_1^1$, because it depends on the interface term with an element of group $R$, whose stage (s3) in (7) is computed at $t^{n+(1/2)}$.

The solution $u_2^{n+1}$ is obtained by

$$ (s1): u_2^{(1)} = u_2^n; \quad K_2^1 = f(u_2^{(1)}, t^n) $$

$$ (s2): u_2^{(2)} = u_2^n + \Delta t K_2^1; \quad K_2^2 = f(u_2^{(2)}, t^n + \Delta t) $$

$$ (s3): u_2^{(3)} = u_2^{(1)}; \quad K_2^3 = K_2^1 $$

$$ (s4): u_2^{(4)} = u_2^{(2)}; \quad K_2^4 \neq K_2^2 $$

$$ u_2^{n+1} = u_2^n + \Delta t \left( \frac{1}{4} K_2^1 + \frac{1}{4} K_2^2 + \frac{1}{4} K_2^3 + \frac{1}{4} K_2^4 \right). \quad (10) $$

In this case, $K_2^4 \neq K_2^2$, because it depends on the interface term with an element of group $G$, whose stage (s4) in (9) has been influenced by stage (s3) in (7) computed at $t^{n+(1/2)}$.

RK22<sup>4a</sup> and RK22<sup>4b</sup> have the same number of stages and the same weighting coefficients, as shown in Table I. This allows preservation of flux conservation at the interfaces between these two groups. The resulting multirate scheme is a second-order accurate method that is, stable for the time step $\Delta t = 2\Delta t_m$.

## B. Construction of the Multirate Groups

Let us consider a mesh of the domain defined by the tessellation $\Omega = \cup T$. In this paper, the multirate groups are built for having their stable time steps of ratio $\kappa = 2$. Let us denote $\delta t_m$ and $\delta t_M$ as the time steps according to the CFL of the smallest and biggest elements, respectively.

<figure><img src="image_7.png" /><figcaption>Fig. 2. Example of the sorting of mesh elements into different groups. (Left) Time steps of each element according to their CFL condition. (Right) Sorting of elements in four groups of tag 0, 3, 6, and 9.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 3. Introduction of buffer groups. (Left) First buffer elements are introduced around the group of tag 0. (Right) Four bulk groups of tag 0, 3, 6, and 9, and three buffer groups of tag 2, 5, and 8 are obtained.</figcaption></figure>

The mesh elements are gathered in different groups that satisfy the CFL condition for a certain range of time steps.

If the time step of the multirate strategy is $\Delta t \in [\Delta t_m, \Delta t_M]$, the number of multirate groups $N_g = z_g + 1$, where $z_g$ is an integer defined by $z_g = \log_2(\Delta t/\Delta t_m)$. The multirate group number 0 is constituted by the elements whose time steps are in $[\Delta t, \Delta t_M]$. The multirate group number $z$ is constituted by the elements whose time steps are in $[(\Delta t/2^{z+1}), (\Delta t/2^z)]$, with $z \le z_g$. The time step of the multirate strategy $\Delta t$ that guarantees the stability of the method has to be chosen, such that $\Delta t \in [\Delta t_m, 2^z \Delta t_m]$, where $2^z \Delta t_m \le \Delta t_M$.

For the construction of multirate groups, we choose in the first step to assign a tag $\theta = 3(z_g - z)$ at each group. The example mesh presented in Fig. 2 shows the organization of the elements in four multirate groups. The tag 9 is assigned to group number 0, whose elements have stable time steps in the interval $[8\Delta t_m, 10\Delta t_m]$. The tag 6 is assigned to group number 1, whose elements have stable time steps in the interval $[4\Delta t_m, 8\Delta t_m]$. The tag 3 is assigned to group number 2, whose elements have stable time steps in the interval $[2\Delta t_m, 4\Delta t_m]$. The tag 0 is assigned to group number 3, whose elements have stable time steps in the interval $[\Delta t_m, 2\Delta t_m]$.

The next step is to introduce the buffer groups that serve to accommodate the transition between bulk groups. These buffers have two connected elements, because the base method ERK22 has two stages. To distinguish groups, the tag $\theta = 3(z_g - z) + 2$ is attributed to buffer groups, and the tag $\theta = 3(z_g - z)$ is conserved for bulk groups. The group of tag 0 remains the same, and the buffer elements are successively introduced. During this process, it is ensured that two neighboring elements have neighboring tags.

The illustration presented in Fig. 3 shows the introduction of buffer elements. The first buffer elements are introduced around the bulk group of tag 0. Four bulk groups of tag 0, 3, 6, and 9 and three buffer groups of tag 2, 5, and 8 are constituted.
<figure><img src="image_9.png" /><figcaption>Fig. 4. (Left) Location of the unnecessary buffer elements. (Right) Final multirate partitioning.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 5. (Left) Simulation domain. (Right) Illustration of the electric field that propagates.</figcaption></figure>

<table><caption>TABLE II
CPU TIME OF MULTIRATE SIMULATIONS AND THE GAINS VERSUS
LeapFrog SCHEME (GAIN<sub>LF</sub>) AND THE SINGLE-RATE
SCHEME (GAIN<sub>RK22</sub>)</caption><thead><tr><th>N<sub>g</sub></th><th>Δt(s)</th><th>buffer elts</th><th>CPU<sub>tot</sub>(s)</th><th>Gain<sub>LF</sub></th><th>Gain<sub>RK22</sub></th></tr></thead><tbody><tr><td>2</td><td>2Δt<sub>m</sub></td><td>1%</td><td>137100</td><td>1.02</td><td>1.85</td></tr><tr><td>3</td><td>2<sup>2</sup>Δt<sub>m</sub></td><td>14%</td><td>71100</td><td>1.97</td><td>3.57</td></tr><tr><td>4</td><td>2<sup>3</sup>Δt<sub>m</sub></td><td>34%</td><td>49200</td><td>2.84</td><td>5.16</td></tr><tr><td>5</td><td>2<sup>4</sup>Δt<sub>m</sub></td><td>38%</td><td>45200</td><td>3.09</td><td>5.62</td></tr><tr><td>6</td><td>2<sup>5</sup>Δt<sub>m</sub></td><td>39%</td><td>42040</td><td>3.33</td><td>6.04</td></tr></tbody></table>

The last step consists of transferring the unnecessary buffer elements to the bulk group with an inferior tag. These elements appear at the boundary and are not located between the two bulk groups. Fig. 4 shows the unnecessary buffer elements corresponding to the groups of tag 2 and 8 (dotted elements), which are transferred to groups 0 and 6.

## IV. NUMERICAL EXAMPLE

The simulation domain is a box $B$ ($\epsilon_B = \epsilon_0$, and $\sigma_B = 0$ S/m), containing a skull $S_K$ ($\epsilon_{S_K} = 32\epsilon_0$, and $\sigma_{S_K} = 0.8$ S/m) that is inside a sphere $S_1$ ($\epsilon_{S_1} = 14\epsilon_0$, and $\sigma_{S_1} = 0.1$ S/m), as shown in Fig. 5. Inside the skull, there is another sphere $S_2$ ($\epsilon_{S_2} = 49\epsilon_0$, and $\sigma_{S_2} = 1.15$ S/m). The permeability of vacuum $\mu_0$ is assumed in all the domains.

The faces of the box oriented along the y-axis are magnetic walls, and those oriented along the z-axis are electric walls. On the left-hand side, an absorbing Silver-Müller boundary condition is set. An incident field $\tilde{E}_{\text{inc}} = (0, 0, E_z^{\text{inc}})$ propagates from the right-hand side. It is a Gaussian-modulated pulse of center frequency $f_0 = 1.2$ GHz, such that $E_z^{\text{inc}}(t) = \sin(2\pi f_0 t)e^{-a(t-t_0)^2}$, $t_0 = 2.5$ ns and $a = 15 \cdot 10^8$. The simulations are carried out for a duration $T_f = 25$ ns. They are performed on an 8core–15Go/ RAM–2.3 GHz computer.

<figure><img src="image_11.png" /><figcaption>Fig. 6. Comparison of $E_z$ recorded at the center of the sphere $S_2$.</figcaption></figure>

A tetrahedral mesh of 279 802 third-order spatial elements is used. The minimum and maximum elementwise stable time steps according to the CFL condition are $\Delta t_m = 24 \cdot 10^{-14}$ s and $\Delta t_M = 17 \cdot 10^{-12}$ s, respectively. The single-rate RK22 and the LeapFrog scheme are first used with $\Delta t_m$. Their CPU times are $\text{CPU}_{\text{RK22}} = 254 100$ s and $\text{CPU}_{\text{LF}} = 140 100$ s. They are compared with those obtained by the multirate strategy in Table II to show the computational gain. The setup of the groups takes a few seconds, which are neglected in the CPU times. When $N_g$ increases, the percentage of the buffer elements also increases, and the gain is improved. The higher stable time step is $2^5 \Delta t_m$, because $2^6 \Delta t_m > \Delta t_M$. The electric field recorded at the center of the sphere $S_2$ is plotted in Fig. 6 to show that the results of the three methods are very close.

## V. CONCLUSION

In this paper, a multirate technique to reduce the computational cost of the explicit computations of the time-domain Maxwell equations on complex geometries is presented. A numerical example shows the gain obtained, compared with the simulations performed by the classical LeapFrog and Heun methods.

## REFERENCES

[1] S. Schomann, N. Gödel, T. Warburton, and M. Clemens, “Local timestepping techniques using Taylor expansion for modeling electromagnetic wave propagation with discontinuous Galerkin-FEM,” *IEEE Trans. Magn.*, vol. 46, no. 8, pp. 3504–3507, Aug. 2010.

[2] S. Descombes, S. Lanteri, and L. Moya, “Locally implicit time integration strategies in a discontinuous Galerkin method for Maxwell’s equations,” *J. Sci. Comput.*, vol. 56, no. 1, pp. 190–218, Jul. 2013.

[3] E. M. Constantinescu and A. Sandu, “Multirate timestepping methods for hyperbolic conservation laws,” *J. Sci. Comput.*, vol. 33, no. 3, pp. 239–278, Dec. 2007.

[4] N. Gödel, S. Schomann, T. Warburton, and M. Clemens, “GPU accelerated Adams–Bashforth multirate discontinuous Galerkin FEM simulation of high-frequency electromagnetic fields,” *IEEE Trans. Magn.*, vol. 46, no. 8, pp. 2735–2738, Aug. 2010.

[5] B. Seny, J. Lambrechts, R. Comblen, V. Legat, and J.-F. Remacle, “Multirate time stepping for accelerating explicit discontinuous Galerkin computations with application to geophysical flows,” *Int. J. Numer. Methods Fluids*, vol. 71, no. 1, pp. 41–64, Jan. 2013.

[6] B. Seny, J. Lambrechts, T. Toulorge, V. Legat, and J.-F. Remacle, “An efficient parallel implementation of explicit multirate Runge–Kutta schemes for discontinuous Galerkin computations,” *J. Comput. Phys.*, vol. 256, pp. 135–160, Jan. 2014.

[7] L. Fezoui, S. Lanteri, S. Lohrengel, and S. Piperno, “Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes,” *ESAIM, Math. Model. Numer. Anal.*, vol. 39, no. 6, pp. 1149–1176, 2005.
## 2.6 Modélisation de boitiers de blindage

En compatibilité électromagnétique, la protection des systèmes contre les aggressions ex- térieures est une question importante pour en garantir la fiabilité des systèmes. Dans le cas des dispositifs électromagnétiques du génie électrique, cette protection consiste en des solu- tions pour contrôler les niveaux d’émission et de reception d’ondes électromagnétiques. Pour plusieurs systèmes comme par exemple l’électronique embarquée, cette protection est réalisée par des boitiers fabriqués à partir de matériaux homogènes ou hétérogènes. Les parois de ces boitiers sont généralement fines et présentent des ouvertures, servant au passage des câbles de connections ou encore à la ventilation. Le dimensionnement de ces boitiers s’appuie sur l’éva- luation des champs électromagnétiques à l’intérieur et à l’extérieur ce qui permet de déterminer l’efficacité de blindage de ceux-ci. Mais lorsque les fréquences de fonctionnement sont éle- vées, la présence de plaques fines conductrices introduit une difficulté de modélisation à cause des faibles valeurs de l’épaisseur de peau. En effet l’évaluation rigoureuse des champs élec- tromagnétiques dans ces cas passe par la prise en compte de la décroissance des champs dans l’épaisseur de peau. Dans les modèles numériques d’étude de boitiers, la solution la plus simple consiste à mailler finement ces plaques, ce qui malheureusement conduit à des systèmes très couteux à résoudre. Des approches ont été proposées pour éviter cette difficulté lors de la mo- délisation de ces boitiers. L’une des plus répandues consiste en la représentation du boitier sous forme de circuit électrique équivalent et à l’utilisation de la méthode des lignes de transmission. Elle a donné lieu à des développements analytiques [88, 119] et numériques [4, 48]. Malheureu- sement, dans cette approche, les parois conductrices sont considérées comme des conducteurs parfaits et seul le rayonnement à travers les ouvertures est modélisé. Une autre pour simplifier la modélisation de ces boitiers est basée sur l’utilisation de conditions d’impédances de surfaces. Ces conditions sont couramment utilisées pour traiter les problèmes de diffractions d’ondes électromagnétiques par des objets partiellement ou totalement pénétrables [36, 98, 123]. C’est autour de ce type de conditions que nous avons réalisé nos développements.

### 2.6.1 Interface simulant les plaques minces homogènes

L’utilisation des conditions d’impédances dans les travaux de modélisation en électroma- gnétisme a fait l’objet de nombreux développements numériques [54, 69, 92]. Ce type de condi- tion qui permet d’éviter d’avoir à mailler les parties conductrices de faibles épaisseurs, a été implémentée en éléments finis et différences finies, pour faciliter la modélisation de nombreux dispositifs en domaine fréquentiel et temporel [40, 41, 93].

Notre contribution a été d’introduire ce type de condition dans une méthode Galerkin discon- tinue en domaine temporel en vue du calcul de l’efficacité de blindage des boitiers [RI.20].
La construction de cette condition repose sur l'hypothèse selon laquelle la propagation des grandeurs électromagnétiques dans une plaque mince conductrice est mono-dimensionnelle et peut être décrite à partir d'expressions analytiques. Notre travail d'implémentation a consisté à introduire, dans la méthode de Galerkin discontinue, les expressions analytiques permettant d'évaluer les champs réfléchis et transmis par une plaque. Pour construire ces expressions, on utilise les résultats du modèle analytique 1D d'une plaque conductrice en domaine fréquentiel que l'on transporte dans le domaine temporel grâce à une transformée de Fourier inverse. Comme l'évolution des champs électromagnétiques dans la plaque est décrite analytiquement, celle-ci est enlevée et remplacée par une interface dotée d'une condition de type impédance de surface. L'utilisation d'une telle interface permet de réduire la taille du maillage (Fig.2.1) et Fig.2.2), ce qui a un impact sur le temps de simulation.

En pratique, introduire cette condition dans la méthode Galerkin discontinue consiste à définir des termes de flux spécifiques pour l'interface qui remplace la plaque conductrice. Pour cela, établit des relations entre les champs de part et d'autre de l'interface ($\mathbf{E}^{-}$, $\mathbf{H}^{-}$) et ($\mathbf{E}^{+}$, $\mathbf{H}^{+}$), en utilisant les expressions analytiques reliant les champs aux bords intérieurs de la plaque ($\mathbf{E}^{0}$, $\mathbf{H}^{0}$) et ($\mathbf{E}^{d}$, $\mathbf{H}^{d}$) (Fig.2.3).

Ces travaux ont fait l'objet de publications dont celle présentée ci-après **3-D Modeling of Thin Sheets in the Discontinuous Galerkin Method for Transient Scattering Analysis** parue dans *IEEE Transactions On Magnetics* qui présente la construction, l'implémentation et la validation de cette condition d'interface.

<figure><img src="image_12.png" /><figcaption>FIGURE 2.1 – Paroi et interface</figcaption></figure>
# 3-D Modeling of Thin Sheets in the Discontinuous Galerkin Method for Transient Scattering Analysis

Mohamed Boubekeur¹, Abelin Kameni¹, Laurent Bernard¹, Axel Modave², and Lionel Pichon¹

¹Laboratoire de Genie Electrique de Paris, SUPELEC, Université Paris Sud et Université Pierre et Marie Curie, Gif sur Yvette cedex 91192, France

²Department of Electrical Engineering and Computer Science, Montefiore institute, Université de Liège, Liège 4000, Belgique

This paper presents a modeling of thin sheets. An interface condition based on analytical solution is used to avoid a fine mesh. This condition is integrated in a time-domain discontinuous Galerkin method to evaluate the shielding effectiveness. This approach is validated by a comparison with analytical solution. 2-D and 3-D cavities are simulated to illustrate the efficiency of the condition.

**Index Terms**—Discontinuous Galerkin (DG) methods, electromagnetic compatibility, interface condition, thin sheet.

## I. INTRODUCTION

MANY problems in electromagnetic compatibility require adequate numerical approaches to evaluate shielding effectiveness. The ability to model features that are small relative to the cell size is often important in electromagnetic simulations. This may lead to subsequent increase in memory and execution time due to a refined mesh around small details of the geometry. In addition, the quality of the mesh can be strongly affected. The impact of thin sheets is even more so consequent for the time-domain computations.

To avoid the spatial discretization of thin sheets, different interface conditions have been proposed. In the frequency domain, interface condition is developed using analytical solutions [1]. This condition can be included in a 3-D model [2]. In the time domain, an inverse Fourier or Laplace transform is combined to a convolution product. Many papers have been devoted to implement this approach in the FDTD method [3]–[5]. Nevertheless, the stair casing error present in the FDTD method may affect significantly the numerical results.

The discontinuous Galerkin (DG) method is a powerful approach for solving time-dependent problems [6]. It is based on the local solution of the equations in each cell and uses flux terms to connect adjacent elements. It has the advantage of the unstructured mesh and high spatial order scheme unlike the conventional FDTD. Such a high spatial scheme can reduce the dispersive error induced by the low level of the spatial approximation in the FDTD.

In this paper, a sheet interface condition is built to replace a thin sheet. A general time-domain formulation is proposed to consider the physical characteristics of the sheet. The approach is based on rational approximation and a recursive algorithm to compute the convolution products. For low frequencies (when the skin depth is larger than the thickness), a simplified expression is derived.

This interface condition is implemented in a DG module of the software GMSH [7].

## II. PROBLEM FORMULATION

Maxwell's equation is solved in the time domain

$$ egin{cases} \epsilon \partial_t \vec{E} - \nabla \times \vec{H} = -\vec{J} \\ \mu \partial_t \vec{H} + \nabla \times \vec{E} = 0 \end{cases} \qquad (1) $$

where $\epsilon$ is the permittivity of the medium and $\mu$ is its permeability. The current density of the conductive medium is such as $\vec{J} = \sigma \vec{E}$, with $\sigma$ the conductivity.

The DG variational formulation of (1) in each element T is given by [6]

$$ egin{cases} \int_T \epsilon \partial_t \vec{E} \cdot \vec{\phi} - \int_T \vec{H} \cdot \nabla \times \vec{\phi} - \int_{\partial T} (\vec{n} \times \vec{H})^{\text{num}} \cdot \vec{\phi} = - \int_T \sigma \vec{E} \cdot \vec{\phi} \\ \int_T \mu \partial_t \vec{H} \cdot \vec{\psi} + \int_T \vec{E} \cdot \nabla \times \vec{\psi} + \int_{\partial T} (\vec{n} \times \vec{E})^{\text{num}} \cdot \vec{\psi} = 0 \end{cases} \qquad (2) $$

where $\vec{\phi}$ and $\vec{\psi}$ are the test functions.

The interface terms between adjacent elements are evaluated using numerical fluxes $(\vec{n} \times \vec{E})^{\text{num}}$ and $(\vec{n} \times \vec{H})^{\text{num}}$ given by [6]

$$ egin{cases} (\vec{n} \times \vec{H})^{\text{num}} = \vec{n} \times \frac{\{Z\vec{H}\}}{\{Z\}} - \alpha(\vec{n} \times \frac{\{\vec{n} \times \{\vec{E}\}\}}{\{Z\}}) \\ (\vec{n} \times \vec{E})^{\text{num}} = \vec{n} \times \frac{\{Y\vec{E}\}}{\{Y\}} + \alpha(\vec{n} \times \frac{\{\vec{n} \times \{\vec{H}\}\}}{\{Y\}}) \end{cases} \qquad (3) $$

where $Z = 1/Y = \sqrt{\mu/\epsilon}$, $[u] = u^+ - u^-/2$ and $\{u\} = u^+ + u^-/2$.

The – denotes the values in the current element and + for the adjacent element. For $\alpha = 0$, centered fluxes are obtained and the numerical scheme is only dispersive [8]. For $\alpha = 1$, upwind fluxes are obtained and the numerical scheme is additionally dissipative [6].

For simulations, an implicit second-order Runge-Kutta scheme that permits a large time step is considered for time discretization, and a first spatial order with upwind fluxes is used for the spatial scheme.

Manuscript received June 28, 2013; revised September 19, 2013; accepted September 25, 2013. Date of current version February 21, 2014. Corresponding author: M. Boubekeur (e-mail: mohamed.boubekeur@lgep.supelec.fr). Digital Object Identifier 10.1109/TMAG.2013.2284505
<figure><img src="image_13.png" /><figcaption>Fig. 1. Geometry of the thin sheet.</figcaption></figure>

### III. CONSTRUCTION OF THE INTERFACE CONDITION

The electromagnetic field equations in the sheet are considered in 1-D and in the frequency domain [1]

$$ \begin{cases} \partial_x E_y = -j\omega\mu H_z \\ \partial_x H_z = -(\sigma + j\omega\epsilon)E_y \end{cases} \qquad (4) $$

where $\epsilon$, $\mu$, and $\sigma$ are, respectively, the permittivity, the permeability, and the conductivity of the sheet (Fig. 1), its thickness is $d$.

Using the analytical 1-D solution of (4), the electromagnetic fields on the both sides of the sheet are connected by

$$ \begin{pmatrix} E_y(\omega, d) \\ H_z(\omega, d) \end{pmatrix} = \begin{pmatrix} \cosh(\gamma d) & \eta \sinh(\gamma d) \\ \eta^{-1} \sinh(\gamma d) & \cosh(\gamma d) \end{pmatrix} \begin{pmatrix} E_y(\omega, 0) \\ H_z(\omega, 0) \end{pmatrix} \quad (5) $$

where $\gamma$ is the planar propagation constant

$$ \gamma = \sqrt{\mu\omega(j\sigma - \epsilon\omega)} \qquad (6) $$

and $\eta$ is the intrinsic impedance of the sheet

$$ \eta = \sqrt{\frac{\mu\omega}{\epsilon\omega - j\sigma}} \qquad (7) $$

This relation can be rewritten as an impedance relation for a layer

$$ \begin{pmatrix} E_y(\omega, 0) \\ E_y(\omega, d) \end{pmatrix} = \begin{pmatrix} \hat{z}_1(\omega) & -\hat{z}_2(\omega) \\ \hat{z}_2(\omega) & -\hat{z}_1(\omega) \end{pmatrix} \begin{pmatrix} H_z(\omega, 0) \\ H_z(\omega, d) \end{pmatrix} \quad (8) $$

with

$$ \begin{cases} \hat{z}_1(\omega) = -\frac{\eta}{\tanh(\gamma d)} \\ \hat{z}_2(\omega) = -\frac{\eta}{\sinh(\gamma d)}. \end{cases} \qquad (9) $$

The time-domain relation of (8) is a convolution product

$$ \begin{pmatrix} E_y(t, 0) \\ E_y(t, d) \end{pmatrix} = \begin{pmatrix} z_1(t) - z_2(t) \\ z_2(t) - z_1(t) \end{pmatrix} * \begin{pmatrix} H_z(t, 0) \\ H_z(t, d) \end{pmatrix} \quad (10) $$

where $z_1$ and $z_2$ are, respectively, the inverse Fourier transform (IFT) of $\hat{z}_1$ and $\hat{z}_2$.

Let note $(\vec{E}^0, \vec{H}^0)$ and $(\vec{E}^d, \vec{H}^d)$ the fields on the sides of the sheet and $\vec{n}$ the outward unit normal. The relation (10) can be reformulated with tangential components of the fields in the general case and becomes

$$ \begin{cases} \vec{n} \times \vec{E}^0(t) = \vec{n} \times [\vec{n} \times (z_1 * \vec{H}^0(t) - z_2 * \vec{H}^d(t))] \\ \vec{n} \times \vec{E}^d(t) = \vec{n} \times [\vec{n} \times (z_2 * \vec{H}^0(t) - z_1 * \vec{H}^d(t))] \end{cases} \quad (11) $$

### IV. IMPLEMENTATION OF THE SHEET CONDITION

#### A. General Sheet Condition

The implementation of the relation (11) is performed in two steps, the first step consists to compute the IFT of the impedance functions $\hat{z}_1$ and $\hat{z}_2$. The second step consists in computing the convolution production between the impedance functions and the magnetic fields.

The impedance functions $\hat{z}_1$ and $\hat{z}_2$ are characterized by a wideband spectrum. Evaluation of their discrete IFT has a significant time cost and sometimes the results are not very accurate. The equivalent rational approximations are considered in the frequency domain using the vector fitting procedure [9].

The rational approximations of the impedance functions defined in (9) are

$$ \left\{ \begin{aligned} \hat{z}_1(\omega) &= z_{1,0} + \sum_{k=1}^{N_1} \frac{R_{1,k}}{j\omega - p_{1,k}} \\ \hat{z}_2(\omega) &= z_{2,0} + \sum_{k=1}^{N_2} \frac{R_{2,k}}{j\omega - p_{2,k}} \end{aligned} \right. \qquad (12) $$

where $z_{i,0}$ are the asymptotic high-frequency value of $\hat{z}_i(\omega)$, $R_{i,k}$ is the $k$th residue, $p_{i,k}$ is the $k$th pole, and $N_i$ the number of poles of the rational approximation.

The IFT is calculated analytically from the rational approximation, the transient impedance functions are

$$ \left\{ \begin{aligned} z_1(t) &= z_{1,0}\delta(t) + \sum_{k=1}^{N_1} R_{1,k} \exp(p_{1,k}t) \\ z_2(t) &= z_{2,0}\delta(t) + \sum_{k=1}^{N_2} R_{2,k} \exp(p_{2,k}t). \end{aligned} \right. \qquad (13) $$

The convolution integrals are determined numerically using the piecewise linear recursive formulation proposed in [10]. The numerical formulation of the relation (11) for the $m$th time iteration is written as

$$ \left\{ \begin{aligned} \vec{n} \times \vec{E}^{0,m} &= \vec{n} \times \left[ \vec{n} \times \left( \vec{E}_1^{0,m} - \vec{E}_2^{0,m} \right) \right] \\ \vec{n} \times \vec{E}^{d,m} &= \vec{n} \times \left[ \vec{n} \times \left( \vec{E}_2^{0,m} - \vec{E}_1^{d,m} \right) \right] \end{aligned} \right. \qquad (14) $$

with

$$ \left\{ \begin{aligned} \vec{E}_1^{j,m} &= z_{i,0} \vec{H}_1^{j,m} + \sum_{k=1}^{N_i} \vec{\xi}_{i,k}^{j,m} \\ \vec{\xi}_{i,k}^{j,m} &= q_{1,k} \vec{\xi}_{i,k}^{j,m-1} + q_{2,k} \vec{H}_1^{j,m-1} + q_{3,k} \vec{H}_1^{j,m} \\ q_{1,k} &= \exp(\beta_k) \\ q_{2,k} &= \frac{\alpha_k}{\beta_k} [1 + (\beta_k - 1) \exp(\beta_k)] \\ q_{3,k} &= \frac{\alpha_k}{\beta_k} [\exp(\beta_k) - \beta_k - 1] \\ \alpha_k &= \frac{R_{i,k}}{p_{i,k}} \\ \beta_k &= p_{i,k} \Delta t \end{aligned} \right. \qquad (15) $$

for $j \in \{-, +\}$, $i \in \{1, 2\}$, $k \in \{1, 2, ..., N_i\}$ and the time step $\Delta t$.

For an interface representing the sheet, let consider $(\vec{E}^-, \vec{H}^-) = (\vec{E}^0, \vec{H}^0)$ the fields for the current element $-$, and $(\vec{E}^+, \vec{H}^+) = (\vec{E}^d, \vec{H}^d)$ the fields for the adjacent
element +. The resulting flux terms for the current element for this interface are obtained by taking only the magnetic field for the current element, and the electric field is taken from the second expression of the relation (14)

$$ \left\{ \begin{array}{l} (\vec{n} \times \vec{H})^{\text{num},n} = \vec{n} \times \vec{H}^{-,n} \\ (\vec{n} \times \vec{E})^{\text{num},n} = \vec{n} \times \vec{E}^{+,n}. \end{array} \right. \qquad (16) $$

## B. Simplified Sheet Condition

For the low frequencies, the electric fields are supposed constant inside the sheet, and $\epsilon\omega \ll \sigma$.

When $\omega \to 0$, $\omega\epsilon$ can be neglected. The impedance and the propagation constant can be rewritten as

$$ \left\{ \begin{array}{l} \eta = \sqrt{\frac{\jmath\omega\mu}{\sigma}} \\ \gamma = \sqrt{\jmath\omega\mu\sigma}. \end{array} \right. \qquad (17) $$

The asymptotic behavior of $(\hat{z}_1 + \hat{z}_2)$ is

$$ (\hat{z}_1 + \hat{z}_2) = -\frac{2}{\sigma d}. \qquad (18) $$

The relation (10) is rewritten using the average of the electric field

$$ \left\{ \begin{array}{l} \frac{\vec{E}^d(\omega)}{2} = \vec{E}^0(\omega) \\ \frac{\vec{E}^d(\omega) + \vec{E}^0(\omega)}{2} = \frac{1}{2}(\hat{z}_1 + \hat{z}_2) * (\vec{H}^0 - \vec{H}^d). \end{array} \right. \qquad (19) $$

A simple form of (11) is obtained by supposing the average of the electric fields on the sides of the sheet

$$ \left\{ \begin{array}{l} \vec{n} \times \vec{E}^0(t) = \frac{1}{\sigma d} \vec{n} \times \left[ \vec{n} \times (\vec{H}^d(t) - \vec{H}^0(t)) \right] \\ \vec{n} \times \vec{E}^d(t) = \frac{1}{\sigma d} \vec{n} \times \left[ \vec{n} \times (\vec{H}^d(t) - \vec{H}^0(t)) \right]. \end{array} \right. \qquad (20) $$

The resulting flux terms for this simplified sheet condition are

$$ \left\{ \begin{array}{l} (\vec{n} \times \vec{H})^{\text{num}} = \vec{n} \times \vec{H}^{-} \\ (\vec{n} \times \vec{E})^{\text{num}} = \frac{1}{\sigma d} \vec{n} \times \vec{n} \times (\vec{H}^{+} - \vec{H}^{-}). \end{array} \right. \qquad (21) $$

These flux terms are independent of the recursive procedure, and are the same complexity as an ordinary flux terms. This condition is valid when the thickness of the sheet is smaller than the skin depth.

# V. NUMERICAL RESULTS

## A. 1-D Validation Test

To validate the general sheet condition, $\vec{E} = (0, 0, E_z)$, $\vec{H} = (0, H_y, 0)$, and two vacuum mediums separated by a conducting sheet are considered. The permittivity, the permeability, the conductivity, and the thickness of this sheet are, respectively, $\epsilon = 5\epsilon_0$, $\mu = 5\mu_0$, $\sigma = 2000$ S/m, and $d = 1$ mm, where $\epsilon_0$ and $\mu_0$ are the permittivity and the permeability of the vacuum medium (Fig. 1). This sheet is illuminated by an harmonic plane wave for different frequencies. The skin depth varies between 1 and 0.22 mm. Let note $\vec{E}_i, \vec{E}_t$ the incident and the transmitted electric fields. The shielding effectiveness for the electric field ($SE_E$) is defined by the following expression:

$$ SE_E = 20 \log_{10} \left( \frac{|\vec{E}_i|}{|\vec{E}_t|} \right). \qquad (22) $$

<figure><img src="image_14.png" /><figcaption>Fig. 2. Shielding effectiveness.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 3. Transient transmitted field.</figcaption></figure>

Fig. 2 shows the comparison of the $SE_E$ between the solution computed with the general sheet condition and the analytical solution calculated in the frequency domain. The transmitted electric field is taken in the harmonic regime to calculate $SE_E$. The difference is negligible.

The same sheet is illuminated by an incident Gaussian pulse. The minimum of the skin depth is 0.3 mm. The computed solution using the general sheet interface is compared with the solution obtained when the sheet is meshed in Fig. 3. The difference is also negligible and demonstrates that the transient regime is also considered.

This condition is also valid for low conductivities, even if the wavelength is smaller than the skin depth. For $\sigma = 10^{-4}$ S/m, the error of the transmitted electric field between analytical and computed solutions remains less than 1%.

These examples validate the sheet condition for the general case.

## B. 2-D Cavity

Let us consider 50 cm × 40 cm cavity. A 10-cm aperture is located on the front side. The characteristics of the walls are $\epsilon = \epsilon_0$, $\mu = \mu_0$, $\sigma = 50$ S/m, and $d = 1$ mm (Fig. 4). This cavity is illuminated by the Gaussian pulse used in the 1-D example. The incident field is normal to the aperture. This pulse excites the frequencies up to 1 GHz and the minimum of the skin depth is 4.5 mm.

The solution is computed for two cases. The first case using the simplified sheet condition and the second when the sheet is meshed. For the first case, 3000 elements are used. For the second case, 13000 elements are used with 10000 elements only for the walls of 1-mm thickness. A personal computer
<figure><img src="image_16.png" /><figcaption>Fig. 4. 2-D cavity.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 5. Electric field at the center of the cavity in the time domain.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 6. 3-D cavity.</figcaption></figure>

3.05-GHz Pentium is used. The time cost of simulation is 4 h for the first case and 20 h for the second.

In Fig. 5, the electric field in the center of the cavity is shown. This example illustrates the efficiency of the condition for problems with thin conductive walls.

### C. 3-D Cavity

A transient scattering problem is studied. Let consider a 3-D cavity (Fig. 6) [11], whose dimensions are $a = 300$ mm, $b = 120$ mm, $d = 300$ mm, $l = 100$ mm, $w = 5$ mm, and $t = 1$ mm.

This cavity is illuminated by the previous incident Gaussian pulse. Different cases of conductivity of the walls are compared: $\sigma = 200$ S/m, $\sigma = 50$ S/m and perfect electric conductor (PEC). For finite conductivities, the simplified sheet condition is used. A fast Fourier transform is applied to the electric field in the center of the cavity, and the shielding effectiveness is plotted in Fig. 7. The result obtained for the PEC case is similar to the analytical formulation [11]. For the

<figure><img src="image_19.png" /><figcaption>Fig. 7. Shielding effectiveness of the electric field in the center of the cavity in the frequency domain.</figcaption></figure>

finite conductivities, the shielding effectiveness is lower and the resonance peaks are attenuated.

## VI. CONCLUSION

An interface condition to avoid the mesh of thin sheets is presented. It allows a reduction of the computational time cost. The interface condition is implemented in a 3-D time-domain DG method and is validated. Future developments will extend the interface condition for the case of multilayer and heterogeneous material.

## REFERENCES

[1] K. M. Mitzner, "Effective boundary conditions for reflection and transmission by an absorbing shell of arbitrary shape," *IEEE Trans. Antennas Propag.*, vol. 16, no. 6, pp. 706-712, Nov. 1968.

[2] F. Bocquet, L. Pichon, and A. Razek, "3D FEM analysis of electromagnetic wave scattering from a dielectric sheet in EMC problems," *IEEE Trans. Magn.*, vol. 34, no. 5, pp. 2791-2794, Sep. 1998.

[3] M. S. Sarto, "A new model for the FDTD analysis of the shielding performances of thin composite structures," *IEEE Trans. Electromagn. Compat.*, vol. 41, no. 4, pp. 298-306, Nov. 1999.

[4] M. Feliziani, F. Maradei, and G. Tribellini, "Field analysis of penetrable conductive shields by the finite-difference time-domain method with impedance network boundary conditions (INBCs)," *IEEE Trans. Electromagn. Compat.*, vol. 41, no. 4, pp. 307-319, Nov. 1999.

[5] M. Feliziani, "Subcell FDTD modeling of field penetration through lossy shields," *IEEE Trans. Electromagn. Compat.*, vol. 54, no. 2, pp. 299-307, Apr. 2012.

[6] J. S. Hesthaven and T. Warburton, "Nodal high-order methods on unstructured grids: I. Time-domain solution of Maxwell's equations," *J. Comput. Phys.*, vol. 181, no. 1, pp. 186-221, Sep. 2002.

[7] C. Geuzaine and J.-F. Remacle, "Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities," *Int. J. Numer. Methods Eng.*, vol. 77, no. 11, pp. 1309-1331, Sep. 2009.

[8] L. Fezoui, S. Lanteri, S. Lohrengel, and S. Piperno, "Convergence and stability of a discontinuous galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes," *Math. Model. Numer. Anal.*, vol. 39, no. 6, pp. 1149-1176, Nov. 2005.

[9] B. Gustavsen and A. Semlyen, "Rational approximation of frequency domain response by vector fitting," in *Proc. IEEE Power Eng. Soc. Winter Meeting*, Feb. 1998, pp. 1-9.

[10] K. S. Oh and J. E. Shutt-Aine, "An efficient implementation of surface impedance boundary condition for the finite-difference time domain method," *IEEE Trans. Antennas Propag.*, vol. 43, no. 7, pp. 660-666, Jul. 1995.

[11] M. P. Robinson, T. M. Benson, C. Christopoulos, J. F. Dawson, M. D. Ganley, A. C. Marvin, et al., "Analytical formulation for the shielding effectiveness of enclosures with apertures," *IEEE Trans. Electromagn. Compat.*, vol. 40, no. 3, pp. 240-248, Aug. 1998.
<figure><img src="image_20.png" /><figcaption>FIGURE 2.2 – Maillage paroi et interface</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 2.3 – éléments paroi</figcaption></figure>

### 2.6.2 Modélisation 3D de boitiers à parois composites

Le développement et l'utilisation des matériaux composites dans l'industrie offrent des perspectives intéressantes en terme de gain de poids et de place pour de nombreux dispositifs. Ils sont progressivement introduits dans l'industrie aéronautique et automobile notamment pour la fabrication des boitiers servant à la protection des équipements. Les matériaux composites utilisés dans ces boitiers sont généralement constitués de fibres (carbone, verre,...), noyés dans une matrice epoxy [55]. Si ces matériaux offrent l'immense avantage d'être légers tout en étant
résistants, il n’en demeure pas moins que leur modélisation est plus compliquée que celle des matériaux homogènes. Si l’utilisation des matériaux composites, pour des applications mécaniques a été rendue possible par l’expertise acquise dans ce domaine depuis plusieurs décennies [52], l’emploi des composites pour des applications de blindage électromagnétique est en revenue freinée par le manque de modèles. Des modèles analytiques d’homogénéisation pour la modélisation électromagnétique ont été développés au sein de GeePs à travers plusieurs travaux de thèses [80, 84]. Grâce à ces travaux, le milieu hétérogène peut être remplacé par un milieu homogène équivalent dont les paramètres effectifs (conductivité, perméabilité, permittivité) qui dépendent de la fréquence tiennent compte de la structure et des propriétés locales du milieu d’origine.

Pour modéliser en régime fréquentiel les boitiers dont les parois sont en matériaux composites, les propriétés effectives du milieu homogénéisé équivalent sont simplement introduites dans un modèle numérique [80]. Mais lorsque la modélisation est réalisée en domaine temporel comme dans nos travaux, l’utilisation de ces propriétés effectives nécessite des traitements supplémentaires. Notre contribution pour la modélisation de ce type de boitiers a été d’introduire les résultats d’homogénéisation de matériaux composites dans l’approche Galerkin discontinue pour des simulations en domaine temporel [RI.11, RI.19]. Le travail réalisé est une extension de la condition d’interface pour paroi homogène décrite dans la section précédente. Cette implémentation nous à permis de proposer un modèle de calcul de l’efficacité de blindage d’un boitier dont les parois sont faites de deux plaques composites ayant leurs fibres orientées perpendiculairement. Ces résultats sont présentés dans l’article ci-après **Modeling of thin heterogeneous sheets in the discontinuous Galerkin method for 3D transient scattering problems** paru dans *European Physical Journal Applied Physic*.
# Modeling of thin heterogeneous sheets in the discontinuous Galerkin method for 3D transient scattering problems*

Mohamed Boubekeurª, Abelin Kameni Ntichi, and Lionel Pichon

Group of Electrical Engineering – Paris (GeePs), UMR CNRS 8507, CentraleSupelec, Univ. Paris-Sud, Sorbonne Université, UPMC Univ. Paris 06, 11 rue Joliot Curie, Plateau de Moulon, 91192 Gif-sur-Yvette Cedex, France

Received: 30 July 2015 / Received in final form: 6 November 2015 / Accepted: 12 January 2016
Published online: 22 February 2016 – © EDP Sciences 2016

**Abstract.** This paper presents a modeling of heterogeneous sheets in the time domain discontinuous Galerkin method. An homogenization model combined to a sheet interface condition is used to avoid the mesh of these sheets in order to study the transient response of heterogeneous enclosures. The validation of this approach is based on a comparison with the case when the sheet is meshed. To illustrate the efficiency of the interface condition, the simulation of a 3D cavity is performed.

## 1 Introduction

Thanks to the gain of the mass, composite materials are increasingly used in automotive and aircraft industry. If their mechanic behaviour is well controlled, it is not the same for the electromagnetic behaviour, especially for electromagnetic compatibility (EMC) applications. Electronic equipments enclosures are generally made with metal sheets. The presence of metal is easily taken into account in numerical methods by considering the material as a perfect electric conductor (PEC). Natural boundary conditions avoid the mesh of metallic sheets. The replacement of metal by composite materials, such as a resin reinforced by a conductive fibers must be accompanied by appropriate modeling tools to predict their shielding effectiveness (SE). Unfortunately different difficulties arise in the modeling of such materials. For example, take into account features that are small relative to the cell size is difficult. This leads to subsequent increase in memory and execution time due to a refined mesh around small details of geometries. Also the quality of the mesh can be strongly affected. The impact of thin heterogeneous sheets is even more so consequent for the time domain computations. To avoid the mesh of thin sheets, an interface condition based on 1D analytical solution is proposed [1] for homogeneous sheets. This condition has been included in a three dimensional model based on finite element [2]. Thanks to an inverse Fourier transform this condition can be introduced in the time domain computations. It is implemented in FDTD for multilayer sheet [3], and recently included in

the time domain discontinuous Galerkin (DG) method [4]. Another difficulty concerns the dependence of the sheet electromagnetic properties to the microstructures. The effective parameter of composite materials have to be established. Many papers have been devoted for determining the effective parameter of heterogeneous materials in the transverse electric case [5] or in the transverse magnetic case [6,7]. In the transverse magnetic case, the Maxwell-Garnett model is the most popular model to compute the effective parameter in electromagnetic excitation. It provides satisfying results as long as the wavelength is larger than the size of heterogeneities [8]. In references [9,10], a new model is proposed for composite materials illuminated by an electromagnetic wave to provide the effective parameter in high frequency. These works take into account the influence of conductive fibers diluted into a dielectric host matrix under dynamic conditions.

In this paper, the interface condition implemented in the time domain DG method [4] is extended for heterogeneous sheet composed by a resin and reinforced by a conductive fibers. The effective parameters are computed with the dynamic homogenization model presented in reference [10] or the Wiener model [5], and used in the time domain thanks to a Fourier transform combined with a convolution product.

## 2 Problem formulation

Maxwell's equations are solved in the time domain:

$$ \left\{ \begin{array}{l} \epsilon \partial_t E - \nabla \times H = -J \\ \mu \partial_t H + \nabla \times E = 0 \end{array} \right. \qquad (1) $$

<footnote>ª e-mail: mohamed.boubekeur@centralesupelec.fr</footnote>

<footnote>* Contribution to the topical issue “Numelec 2015 – Elected submissions”, edited by Adel Razek</footnote>
where $\epsilon$ is the permittivity of the medium and $\mu$ is the permeability. The current density of the conductive medium is such as $J = \sigma E$, with $\sigma$ the conductivity.

The DG variational formulation of (1) in each element $T$ is given by [11]:

$$ \left\{ \begin{aligned} & \int_T \epsilon \partial_t \mathbf{E} \cdot \mathbf{\phi} - \int_T \mathbf{H} \cdot \nabla \times \mathbf{\phi} \\ & - \int_{\partial T} (\mathbf{n} \times \mathbf{H})^{\text{num}} \times \mathbf{\phi} = - \int_T \sigma \mathbf{E} \cdot \mathbf{\phi} \\ & \int_T \mu \partial_t \mathbf{H} \cdot \mathbf{\psi} + \int_T \mathbf{E} \times \nabla \times \mathbf{\psi} \\ & + \int_{\partial T} (\mathbf{n} \times \mathbf{E})^{\text{num}} \cdot \mathbf{\psi} = 0 \end{aligned} \right. \qquad (2) $$

where $\mathbf{\phi}$ and $\mathbf{\psi}$ are test functions.

The interface terms between adjacent elements are evaluated using numerical fluxes $(\mathbf{n} \times \mathbf{E})^{\text{num}}$ and $(\mathbf{n} \times \mathbf{H})^{\text{num}}$ given by [11]:

$$ \left\{ \begin{aligned} (\mathbf{n} \times \mathbf{H})^{\text{num}} &= \mathbf{n} \times \frac{\{ZH\}}{\{Z\}} - \mathbf{n} \times \frac{(\mathbf{n} \times [\mathbf{E}])}{\{Z\}} \\ (\mathbf{n} \times \mathbf{E})^{\text{num}} &= \mathbf{n} \times \frac{\{Y\mathbf{E}\}}{\{Y\}} + \mathbf{n} \times \frac{(\mathbf{n} \times [\mathbf{H}])}{\{Y\}} \end{aligned} \right. \qquad (3) $$

where $Z = \frac{1}{Y} = \sqrt{\frac{\mu}{\epsilon}}$, $[u] = \frac{u^+ - u^-}{2}$ and $\{u\} = \frac{u^+ + u^-}{2}$.

The “-” denotes the values in the current element, and “+” for the adjacent element.

### 3 The interface condition

Let's consider the one dimensional electromagnetic field equations in the frequency domain. The electromagnetic fields on the both sides of the sheet are connected by:

$$ \begin{pmatrix} E_y(\omega, 0) \\ E_y(\omega, d) \end{pmatrix} = \begin{pmatrix} \hat{z}_1(\omega) - \hat{z}_2(\omega) \\ \hat{z}_2(\omega) - \hat{z}_1(\omega) \end{pmatrix} \begin{pmatrix} H_z(\omega, 0) \\ H_z(\omega, d) \end{pmatrix} \qquad (4) $$

with:

$$ \left\{ \begin{aligned} \hat{z}_1(\omega) &= -\frac{\eta}{\tanh(\gamma d)} \\ \hat{z}_2(\omega) &= -\frac{\eta}{\sinh(\gamma d)} \end{aligned} \right. \qquad (5) $$

and

$$ \left\{ \begin{aligned} \gamma &= \sqrt{\mu\omega(\gamma\sigma - \epsilon\omega)} \\ \eta &= \sqrt{\frac{\mu\omega}{\epsilon\omega - \gamma\sigma}} \end{aligned} \right. \qquad (6) $$

with $d$ the thickness of the sheet, and $\epsilon, \mu, \sigma$ are respectively the permittivity, the permeability and the conductivity of the sheet (Fig. 1).

The time domain equivalent relation of (4) is a convolution product:

$$ \begin{pmatrix} E_y(t, 0) \\ E_y(t, d) \end{pmatrix} = \begin{pmatrix} z_1(t) - z_2(t) \\ z_2(t) - z_1(t) \end{pmatrix} * \begin{pmatrix} H_z(t, 0) \\ H_z(t, d) \end{pmatrix}, \qquad (7) $$

where $z_1$ and $z_2$ are respectively the inverse Fourier transform (IFT) of $\hat{z}_1$ and $\hat{z}_2$.

<figure><img src="image_22.png" /><figcaption>Fig. 1. Geometry of the thin sheet.</figcaption></figure>

Let's note $(\mathbf{E}^0, \mathbf{H}^0)$ and $(\mathbf{E}^d, \mathbf{H}^d)$ the fields on the sides of the sheet and $\mathbf{n}$ the outward unit normal. The relation (7) can be reformulated with tangential components of the fields and becomes:

$$ \left\{ \begin{aligned} \mathbf{n} \times \mathbf{E}^0(t) &= \mathbf{n} \times [\mathbf{n} \times (z_1 * \mathbf{H}^0(t) - z_2 * \mathbf{H}^d(t))] \\ \mathbf{n} \times \mathbf{E}^d(t) &= \mathbf{n} \times [\mathbf{n} \times (z_2 * \mathbf{H}^0(t) - z_1 * \mathbf{H}^d(t))] \end{aligned} \right. \qquad (8) $$

In the DG formulation, the interface condition representing a thin sheet is constructed as follows: let's consider $(\mathbf{E}^-, \mathbf{H}^-) = (\mathbf{E}^0, \mathbf{H}^0)$ the fields for the current element “-”, and $(\mathbf{E}^+, \mathbf{H}^+) = (\mathbf{E}^d, \mathbf{H}^d)$ the fields for the adjacent element “+”. The magnetic field is taken from the adjacent element, and the electric field is taken from the second expression of the relation (8).

$$ \left\{ \begin{aligned} (\mathbf{n} \times \mathbf{H})^{\text{num},n} &= \mathbf{n} \times \mathbf{H}^{-,n} \\ (\mathbf{n} \times \mathbf{E})^{\text{num},n} &= \mathbf{n} \times \mathbf{E}^{+,n}. \end{aligned} \right. \qquad (9) $$

### 4 The dynamic homogenization model

The homogenization consists in defining a fictive homogeneous medium with the same behaviour than the heterogeneous medium. The effective parameter are determined in the frequency domain. Let's note $\epsilon^*$ the complex permittivity:

$$ \epsilon^* = \epsilon - j\frac{\sigma}{\omega}. \qquad (10) $$

For a medium constituted by a resin reinforced by a conductive fibers (Fig. 2), the effective permittivity in direction perpendicular to the fibers is given by [9]:

$$ \epsilon_{\perp}^{*} = \frac{\epsilon_{1}^{*} \frac{f_{1}}{\epsilon_{\infty}^{*} + \epsilon_{1}^{*}} + \epsilon_{2}^{*} \frac{f_{2}}{\epsilon_{\infty}^{*} + \epsilon_{2}^{*}}}{\frac{f_{1}}{\epsilon_{\infty}^{*} + \epsilon_{1}^{*}} + \frac{f_{2}}{\epsilon_{\infty}^{*} + \epsilon_{2}^{*}}}, \qquad (11) $$

where $\epsilon_1^*$ et $f_1$ are respectively the complex permittivity and volume fraction of the resin, and $\epsilon_2^*$ and $f_2$ are respectively the complex permittivity and volume fraction of the fibers, $\epsilon_\infty^*$ is the complex permittivity of the infinite medium, suggested in reference [12].
<figure><img src="image_23.png" /><figcaption>Fig. 2. Geometry of the composite medium.</figcaption></figure>

In this study, the expression of $\epsilon_{\infty}^*$ proposed in reference [10], is given by:

$$ \epsilon_{\infty}^* = \epsilon_1^* - j \frac{\sigma_{\infty}}{\omega}, \quad (12) $$

with

$$ \frac{\sigma_{\infty}}{\omega} = \frac{2}{R^2}(1 + f_2) \int_0^R r \left| \frac{\mathcal{J}_1(r\gamma_2)}{\mathcal{J}_0(R\gamma_2)} \right| dr, \quad (13) $$

where R is the radius of the fiber, $\gamma_2$ the propagation constant of the fiber defined in equation (6), $\mathcal{J}_\alpha$ is the first kind Bessel function of order $\alpha$.

The effective permittivity of the homogenized medium $\epsilon_\perp$ and his conductivity $\sigma_\perp$ are obtained from equation (11) and are used to compute the impedance functions $\hat{z}_1$ and $\hat{z}_2$ given in equation (5). These parameters depend on the singular frequency $\omega$ and the properties of the heterogeneous medium ($f_1, f_2, R, \epsilon_1, \epsilon_2, \sigma_1, \sigma_2$).

## 5 Validation example

Let's consider $E = (0, E_y, 0)$, $H = (0, 0, H_z)$, a composite sheet constituted by a resin ($\epsilon_1 = 2\epsilon_0, \mu_0, \sigma_1 = 0.01$ S/m) reinforced by 30 conductive fibers ($\epsilon_2 = \epsilon_0, \mu_0, \sigma_2$) located in the vacuum ($\epsilon_0, \mu_0, \sigma_0 = 0$ S/m). The volume fraction of the resin is $f_1 = 80.37\%$, and the volume fraction of the conductive fibers is $f_2 = 19.63\%$. The radius of the fiber is $R = 0.05$ mm and the thickness of the sheet is $d = 6$ mm (Fig. 3). This sheet is illuminated by a Gaussian pulse including frequencies up to 100 GHz. Two values of $\sigma_2$ are considered, $\sigma_2 = 1000$ S/m and $\sigma_2 = 10000$ S/m. The fibers are oriented in the direction $e_z$, and the incident wave is oriented in the direction $e_x$.

The solutions are computed for two cases. The first case uses the interface condition with the effective permittivity and the effective conductivity computed with the dynamic homogenization model. The size of the mesh is 1346 elements, The effective parameters are computed for the composite sheet described previously (Fig. 3). The effective permittivity and the effective conductivity are computed for the same configuration presented in Figure 3. They are presented in Figures 4 and 5. The effective permittivity can be considered constant, but the effective conductivity is not constant.

In the second case, the sheet is meshed. The size of the mesh is 4350 elements. A computer with 32 CPU

<figure><img src="image_24.png" /><figcaption>Fig. 3. Computational domain for the validation example.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 4. Effective permittivity of the sheet for $\sigma_2 = 10000$ S/m.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 5. Effective conductivity of the sheet for $\sigma_2 = 10000$ S/m.</figcaption></figure>

(2 GHz Pentium Xeon, Intel, Santa Clara, California, USA) is used. The CPU time is 4 min for the first case and 25 min for the second case.

For $\sigma_2 = 1000$ S/m, the transmitted electric field in the time domain is the same for the two cases (Fig. 6). The shielding effectiveness of the electric field ($\text{SE}_E$) is defined by the following expression:

$$ \text{SE}_E = 20 \log_{10} \left( \frac{|\mathbf{E}_i|}{|\mathbf{E}_t|} \right). \quad (14) $$

The $\text{SE}_E$ is presented in Figure 7. Three cases are compared. The solution computed when the sheet is meshed, the solution computed with the interface condition, and the analytical solution computed using the effective parameters. The difference is negligible. This example validates the interface condition for heterogeneous sheet in the time domain.

For $\sigma_2 = 10000$ S/m, the transmitted fields for the case when the interface condition and the case when the sheet is meshed are plotted in Figure 8. They have the same behaviour.

The comparison of the shielding effectiveness shows (Fig. 9) that the interface condition is similar to the analytical solution using the effective parameters.