UNIVERSITE LIBRE DE BRUXELLES Année académique 2009- 2010 Faculté des Sciences appliquées Service Aéro-Thermo-Mécanique Méthode de couplage éléments finis/éléments de frontière avec accélération pour l'équation d'induction. Directeur de Mémoire : Prof. Gérard Degrez Personne ressource : Ir. Thomas Cordaro Mémoire de fin d'étude présenté par Xavier Dechamps en vue de l'obtention du diplôme de Master en Sciences de l’Ingénieur Civil Mécanicien à finalité Aéronautique. Résumé du travail Dans ce présent travail, il sera traité de l’utilisation de méthodes d’éléments de frontière pour la résolution de l’équation d’induction du champ électrique en électromagnétisme. Ce travail s’insère dans le cadre du sujet de doctorat de Thomas Cordaro visant à la modélisation numérique de l’écoulement au sein d’une torche à plasma. Cette torche à plasma, élaborée au von Kármán Institute, fournit les résultats expérimentaux indispensables à la validation du code numérique. L’intérêt de l’utilisation de la méthode des éléments de frontière dans le domaine de l’électromagnétisme est évident. En effet, dans le cas de la méthode des éléments finis, il est indispensable de tenir compte du domaine extérieur en plus du domaine intérieur à la torche afin de respecter la condition de rayonnement infiniment loin de la torche. Or seul l’intérieur de la torche constitue le domaine d’intérêt, ce qui induit une taille du système d’équations beaucoup plus grande qu’elle ne devrait l’être. Dans le cas de la méthode des éléments de frontière, seul le domaine intérieur doit intervenir. Il serait tentant d’en déduire que cette réduction du nombre d’inconnues va mener à un temps de calcul réduit, mais il sera montré dans la suite qu’il n’en est pas toujours le cas. En effet, la méthode des éléments de frontière nécessite énormément de temps pour assembler les matrices correspondant à l’équation intégrale dirigeant le problème sur la frontière de la torche. Dans un but d’accélérer cet assemblage, la méthode multipôles semble toute indiquée. La méthodologie suivie dans ce travail est la suivante : Encoder et valider la méthode des éléments finis (Finite Element Method - • FEM) dont la solution servira de référence pour la suite (point de vue temps de calcul). Le chapitre II se concentrera sur l’étude des résultats issus de cette méthode. Encoder et valider une méthode de couplage entre les méthodes des éléments • finis et des éléments de frontière (Boundary Element Method - BEM). La solution issue de la méthode BEM sur le contour de la torche servira de condition aux limites pour la méthode FEM, utilisée pour déterminer les inconnues au sein de la torche. Ceci constituera le sujet du chapitre III. Le principal désavantage de la méthode BEM réside en la production de • matrices pleines et donc un temps d’assemblage et d’inversion assez conséquent. Une manière d’accélérer le processus est de se tourner vers la méthode multipôles, méthode émergente depuis quelques années. Le chapitre IV sera dédié à l’étude de cette méthode. Les différents codes sont écrits sous format Matlab®. Ce choix a été grandement facilité par l’aisance de l’utilisation des matrices par ce logiciel. Par contre, comme l’affichage d’un grand nombre de surfaces élémentaires est relativement malaisé avec Matlab®, la représentation des différents résultats sera faite à l’aide de Gmsh ([GEU09]), un outil de maillage tridimensionnel et de post-traitement. La génération des différents maillages se fera également par l’intermédiaire de ce dernier logiciel. Remerciements Je voudrais commencer en remerciant tous ceux qui ont contribué à leur manière à l’élaboration de ce travail. Tout d’abord, je tiens à remercier mon promoteur, le professeur Gérard Degrez, source intarissable de connaissances, de conseils et d’anecdotes. Je suis constamment étonné par son dynamisme et sa capacité à fournir des explications aux situations les plus complexes. Que serait mon mémoire sans Thomas Cordaro ? Merci beaucoup, Thomas ! Je te serai éternellement reconnaissant pour ta patience toujours présente. J’espère seulement ne pas en avoir abusé. Petit mot personnel : mon antédiluvien PC a quand même tenu le coup jusque la fin, même s'il m'a causé quelques frayeurs. Vient maintenant le tour de ma famille. Je me rends compte maintenant des efforts réalisés par mes parents pour supporter ma mauvaise humeur malheureusement assez souvent présente. Par chance, mon grand frangin a passé le cap du mémoire avant moi, mes parents ont ainsi pu exercer leur talent de remonteurs de moral. Merci. Merci beaucoup, 'pa et 'man. Merci également à toi, Yves, pour ton soutien et tes conseils. Jessica, Jessica, Jessica. Tu t'es souvent retrouvée aux premières lignes de mes explosions de sentiments lors de cette dernière année d’études. Je ne saurais dire dans quel état je serais si tu n’étais pas là à mes côtés. Je te remercie de tout mon cœur. Je tiens également à remercier les autres membres de l’équipe ATM. Axel, pour tes blagues et piques toujours présentes, Matthew, pour tes chemises toujours aussi originales, ainsi que tous les autres pour leur magnifique accueil. Table des matières HISTOIRE ET PHYSIQUE DU PROBLEME................................................................................................1 1. UN PEU D'HISTOIRE...............................................................................................................................1 2. PHYSIQUE DU PROBLEME......................................................................................................................1 METHODE DES ELEMENTS FINIS..............................................................................................................5 1. INTRODUCTION.....................................................................................................................................5 2. DISCRETISATION PAR LA METHODE DES ELEMENTS FINIS......................................................................6 2.1. Discrétisation de l’équation du champ induit.............................................................................7 2.2. Discrétisation de l’équation du champ total...............................................................................8 3. ASSEMBLAGE DES MATRICES ELEMENTAIRES.......................................................................................9 4. DEVELOPPEMENTS DES MATRICES ELEMENTAIRES..............................................................................10 5. PARTICULARITES DU MAILLAGE..........................................................................................................12 6. CONDUCTIVITE ELECTRIQUE NULLE DANS TOUT LE DOMAINE.............................................................14 6.1. Détermination de la taille du domaine extérieur.......................................................................14 6.2. Validation du code implémenté.................................................................................................16 6.3. Performances temporelles du code FEM sur le problème du champ créé par les spires..........19 7. CONDUCTIVITE ELECTRIQUE NON NULLE A L’INTERIEUR DE LA TORCHE.............................................21 7.1. Conséquences de la présence d’une conductivité électrique.....................................................21 7.2. Influence de la conductivité électrique - 27.6 MHz...................................................................22 7.3. Influence de la fréquence d’excitation - σ = 5000 S/m..............................................................25 7.4. Performances temporelles du code FEM sur le problème du champ induit..............................26 8. CONCLUSION......................................................................................................................................27 METHODE DES ELEMENTS DE FRONTIERE........................................................................................28 1. INTRODUCTION...................................................................................................................................28 2. UN PEU D'HISTOIRE.............................................................................................................................29 3. COMPARAISON DES METHODES DES ELEMENTS FINIS ET DES ELEMENTS DE FRONTIERE......................29 4. FONDEMENTS MATHEMATIQUES DE LA METHODE DES ELEMENTS DE FRONTIERE...............................31 4.1. Fonction de Green.....................................................................................................................31 4.2. Equation intégrale pour le champ électrique............................................................................32 4.2.1. Cas général tridimensionnel.................................................................................................................32 4.2.2. Caractère axisymétrique.......................................................................................................................34 4.2.3. Domaine de dimensions infinies............................................................................................................35 5. DEVELOPPEMENTS NUMERIQUES DE LA METHODE DES ELEMENTS DE FRONTIERE...............................36 6. COUPLAGE FAIBLE SUR LE FLUX.........................................................................................................39 7. MODIFICATIONS A APPORTER A LA METHODE DES ELEMENTS FINIS....................................................40 8. CONDUCTIVITE ELECTRIQUE NULLE A L'INTERIEUR DE LA TORCHE - VALIDATION DU CODE...............41 8.1. Validation de la méthode BEM sans couplage..........................................................................42 8.2. Validation de la méthode de couplage FEM / BEM..................................................................43 9. CONDUCTIVITE ELECTRIQUE NON NULLE A L'INTERIEUR DE LA TORCHE.............................................44 9.1. Influence de la conductivité électrique - 27.6 MHz...................................................................45 9.2. Influence de la fréquence d'excitation - σ = 5000 S/m..............................................................48 9.3. Conductivité électrique variable dans le plasma.......................................................................50 10. PERFORMANCES TEMPORELLES DE LA METHODE DE COUPLAGE FEM / BEM.....................................53 11. CONCLUSION......................................................................................................................................56 METHODE MULTIPOLES...........................................................................................................................58 1. INTRODUCTION...................................................................................................................................58 2. UN PEU D’HISTOIRE.............................................................................................................................58 3. DEVELOPPEMENTS MATHEMATIQUES.................................................................................................59 3.1. Référentiel centré sur la spire de courant.................................................................................61 3.2. Référentiel centré sur le centre multipôles................................................................................62 3.3. Moments multipôles...................................................................................................................63 4. DEVELOPPEMENTS NUMERIQUES DE LA METHODE MULTIPOLES.........................................................66 4.1. Discrétisation des moments multipôles.....................................................................................66 4.2. Regroupement des éléments de frontière...................................................................................66 4.3. Méthodologie de la méthode MM-BEM....................................................................................67 5. RESULTATS NUMERIQUES...................................................................................................................69 5.1. Validation de la méthode de couplage MM-BEM/FEM............................................................69 5.2. Influence de la position des centres multipôles.........................................................................70 5.3. Influence du nombre de groupes d'éléments de surface............................................................72 5.4. Performances temporelles.........................................................................................................73 6. CONCLUSION......................................................................................................................................75 CONCLUSION - PERSPECTIVES...............................................................................................................76 BIBLIOGRAPHIE...........................................................................................................................................77 ANNEXE A - DETERMINATION DE L’EQUATION D’INDUCTION DU CHAMP ELECTRIQUE.79 ANNEXE B - FORMULATION THEORIQUE DES CHAMPS ELECTRIQUE ET MAGNETIQUE CREES PAR UNE SPIRE.......................................................................................................81 ANNEXE C - INTEGRATION NUMERIQUE.............................................................................................85 ANNEXE D - FONCTION DE GREEN DE L'EQUATION SCALAIRE DE LAPLACE EN COORDONNEES CYLINDRIQUES.....................................................................................87 ANNEXE E - EQUATION BEM DE L’EQUATION SCALAIRE DE LAPLACE A PARTIR DU CAS TRIDIMENSIONNEL.............................................................................................................91 ANNEXE F - POST-TRAITEMENT - CHAMP MAGNETIQUE..............................................................97 Table des figures Figure 1 : Mini-torche de la VKI lors de son fonctionnement. Sa fréquence d’excitation est de 27.6MHz ([ABE00]).................................................................................................................................2 Figure 2 : Géométrie et dimensions de la mini-torche à plasma ([ABE00])........................................................2 Figure 3 : Méthode des éléments finis. Topologie du maillage et des propriétés des milieux......................6 Figure 4 : Transformation géométrique - définition de l’élément parent...................................................11 Figure 5 : Modèle de maillage adapté aux propriétés physiques du problème...........................................13 Figure 6 : Influence de la taille du domaine extérieur sur la qualité de la solution - Evolution de la solution selon la coordonnée radiale pour une coupe au niveau du conducteur central...........15 Figure 7 : Influence de la taille du domaine extérieur sur la qualité de la solution - Evolution de l’erreur relative selon la coordonnée radiale pour une coupe au niveau du conducteur central...............................................................................................................................................15 Figure 8 : Influence de la taille du domaine extérieur sur la norme L de l’erreur pour les nœuds situés 2 à l’intérieur de la torche..................................................................................................................16 Figure 9 : Evolution de la norme du champ électrique théorique (I.11) [V/m] dans le domaine extérieur à la torche.........................................................................................................................17 Figure 10 : Emission du champ par les spires - Norme du champ électrique total [V/m]. Comparaison entre la solution théorique (au-dessus) et la solution numérique FEM (en dessous).................18 Figure 11 : Variation de la norme L de l’erreur par rapport au nombre de degrés de liberté à 2 l’intérieur de la torche.....................................................................................................................19 Figure 12 : Temps caractéristiques obtenus sur le problème du champ créé par les spires en fonction du nombre de degrés de liberté......................................................................................................20 Figure 13 : Norme du champ électrique total [V/m]. 27.6 MHz - σ = 1000 S/m. Mise en évidence du phénomène d’effet pelliculaire........................................................................................................23 Figure 14 : Phase (°) du champ électrique total. 27.6 MHz - σ = 1000 S/m.................................................23 Figure 15 : Influence de la conductivité électrique. 27.6 MHz. Norme du champ électrique total en fonction de la coordonnée radiale [V/m]. Coupe transversale au niveau du conducteur central...............................................................................................................................................24 Figure 16 : Influence de la fréquence d’excitation. σ = 5000 S/m. Norme du champ électrique total [V/m] en fonction de la coordonnée radiale. Coupe transversale au niveau du conducteur central...............................................................................................................................................25 Figure 17 : Temps caractéristiques obtenus sur le problème du champ induit en fonction du nombre de degrés de liberté..........................................................................................................................26 Figure 18 : Couplage entre les méthodes FEM et BEM. Domaines de résolution pour les méthodes BEM et FEM....................................................................................................................................28 Figure 19 : Définition des domaines intervenant dans le développement de l’équation intégrale au point x ..............................................................................................................................................33 0 Figure 20 : Définition des notations intervenant dans le développement de l’équation intégrale au point x pour un domaine de dimensions infinies...................................................................................35 0 Figure 21 : Fonctions d’interpolation du champ électrique et du flux sur la frontière Γ..........................37 Figure 22 : Conventions pour la description de la condition faible sur le flux...........................................39 Figure 23 : Validation de la méthode des éléments de frontière sans couplage. Représentation de la solution sur la frontière...................................................................................................................42 Figure 24 : Validation du code BEM sans couplage - Norme L de l’erreur en fonction du nombre de 2 nœuds situés sur la frontière...........................................................................................................43 Figure 25 : Validation de la méthode de couplage FEM/BEM - Norme du champ électrique total [V/m]. Solution théorique (au-dessus) et la solution numérique FEM/BEM (en dessous)........43 Figure 26 : Variation de la norme L de l’erreur en fonction du nombre de degrés de libertés à 2 l’intérieur de la torche. Résultats FEM, BEM (6 points d’intégration - vert) et BEM (nombre variables de points d’intégration - rouge)......................................................................................44 Figure 27 : Norme du champ électrique total [V/m]. σ = 1000 S/m - 27.6 MHz. Comparaison des résultats FEM (au dessus) et FEM/BEM (en dessous)..................................................................45 Figure 28 : Phase (°) du champ électrique total. σ = 1000 S/m - 27.6 MHz. Comparaison des résultats FEM (au dessus) et FEM/BEM (en dessous).................................................................................45 Figure 29 : Influence de la conductivité électrique. 27.6 MHz. Norme du champ électrique total [V/m] en fonction de la coordonnée radiale. Solutions FEM et FEM/BEM..........................................46 Figure 30 : Influence de la conductivité électrique. 27.6 MHz. Norme du champ électrique total [V/m] en fonction de la coordonnée radiale. Solutions FEM et FEM/BEM. Echelle logarithmique...47 Figure 31 : Comparaison des résultats FEM avec imposition de la condition de rayonnement à des distances variables et du résultat FEM/BEM. 27.6 MHz et σ = 1000 S/m..................................48 Figure 32 : Influence de la fréquence d’excitation. σ = 5000 S/m. Norme du champ électrique total [V/m] en fonction de la coordonnée radiale. Solutions FEM et BEM.........................................49 Figure 33 : Influence de la fréquence d’excitation. σ = 5000 S/m. Norme du champ électrique total [V/m] en fonction de la coordonnée radiale. Solutions FEM et BEM. Echelle logarithmique..49 Figure 34 : Conductivité électrique non uniforme au sein de la torche. Topologie du problème..............50 Figure 35 : Conductivité électrique non uniforme au sein de la torche. Norme du champ électrique total [V/m]. Solution FEM (au-dessus) et FEM/BEM (en dessous).............................................51 Figure 36 : Conductivité électrique non uniforme au sein de la torche. Phase (°) du champ électrique total. Solution FEM (au-dessus) et FEM/BEM (en dessous)........................................................51 Figure 37 : Conductivité électrique non uniforme au sein de la torche. Norme du champ électrique total [V/m] selon une coupe transversale au niveau du conducteur central (z=0.053m)...........52 Figure 38 : Conductivité électrique non uniforme au sein de la torche. Norme du champ électrique total [V/m] selon une coupe transversale au niveau du conducteur central. Echelle logarithmique...................................................................................................................................52 Figure 39 : Performances temporelles. σ = 1000 S/m - 27.6 MHz. Méthode FEM (II.8) et (II.9). Méthode FEM/BEM avec un nombre de points d'intégration fixe et variable..........................53 Figure 40 : Temps caractéristiques obtenus par la méthode FEM / FEM avec un nombre fixe de points d’intégration en fonction du nombre de degrés de liberté...........................................................54 Figure 41 : Taux d’occupation de la matrice globale [K] [%] en fonction du nombre de degrés de liberté intérieurs à la torche............................................................................................................55 Figure 42 : Schéma explicatif du fonctionnement des méthodes MM-BEM et FMM ([YOS01])....................60 Figure 43 : Notations utilisées pour une spire de courant centrée en l'origine du repère en coordonnées sphériques ([HIR10]).............................................................................................................................61 Figure 44 : Notations utilisées pour le centre du repère sphérique centré sur le centre multipôles.........62 Figure 45 : Regroupement des éléments de surface dans le plan méridien [SIN08].......................................67 Figure 46 : Assemblage de la matrice [A’] par les méthodes BEM traditionnelle et multipôles...............68 Figure 47 : Validation de la méthode de couplage MM-BEM/FEM - Norme du champ électrique total [V/m]. Solution théorique (au-dessus) et la solution numérique MM-BEM/FEM (en dessous) avec 50 groupes d'éléments de frontière........................................................................................70 Figure 48 : Proportion [%] d’opérations effectuées par la méthode multipôles lors de l’assemblage des matrices de l’équation intégrale en fonction de la position des centres multipôles sur l’axe longitudinal......................................................................................................................................71 Figure 49 : Proportion [%] d’opérations effectuées par la méthode multipôles lors de l’assemblage des matrices de l’équation intégrale en fonction du nombre de groupes d’éléments de surface.....72 Figure 50 : Temps caractéristiques obtenus par la méthode MM-BEM / FEM avec un nombre fixe de points d’intégration en fonction du nombre de degrés de liberté................................................73 Figure 51 : Performances temporelles. σ = 1000 S/m - 27.6 MHz. Méthodes FEM (II.8) et (II.9), FEM/BEM avec un nombre de points d'intégration fixe et variable, MM-BEM / FEM...........74 Figure 52 : Topologie du problème du champ créé par une spire...............................................................81 Figure 53 : Définition des coordonnées aréolaires et positions des points d’intégration dans le repère local...................................................................................................................................................85 Figure 54 : Définition du repère cylindrique par rapport au repère cartésien...........................................88 Figure 55 : Norme du champ magnétique théorique [T] créé par les spires dans le vide..........................98 Tableaux Tableau 1 : Profondeurs de peau [mm] théoriques (II.29) pour les différentes configurations envisagées dans les études correspondantes................................................................................22 Tableau 2 : Influence de la conductivité électrique - Profondeurs de peau [mm] théorique (II.29) et numérique......................................................................................................................................24 Tableau 3 : Influence de la fréquence d’excitation - Profondeurs de peau [mm] théorique (II.29) et numérique......................................................................................................................................25 Tableau 4 : Validation du code BEM sans couplage - Norme L de l’erreur de la solution numérique...43 2 Tableau 5 : Influence de la conductivité électrique - Profondeurs de peau [mm] théorique (II.29) et numériques pour les méthodes FEM et FEM/BEM...................................................................46 Tableau 6 : Influence de la fréquence d’excitation - Profondeurs de peau [mm] théorique (II.29) et numériques pour les méthodes FEM et FEM/BEM...................................................................50 Tableau 7 : Temps d’inversion [s] des matrices pour les méthodes FEM et FEM/BEM...........................55 Tableau 8 : Pourcentage du nombre de points d'intégration en fonction du nombre de nœuds situés sur la frontière. σ = 1000 S/m - 27.6 MHz..........................................................................................56 Tableau 9 : Proportions de temps [%] demandées pour l’assemblage des matrices correspondant à l’équation intégrale par rapport à la durée totale. Méthodes FEM/BEM et MM- BEM/FEM......................................................................................................................................74 Tableau 10 : Points d’intégrations sur un élément triangulaire - Table de Hammer................................86 Tableau 11 : Points d’intégrations sur un élément unidimensionnel - Table de Gauss-Legendre([WAR08])86 Chapitre I Histoire et physique du problème 1. Un peu d'histoire L’interaction entre l’écoulement d’un fluide conducteur et un champ électromagnétique, également appelé magnétohydrodynamique (MHD), constitue un sujet d’étude assez ancien. Les premières traces d’expérience à ce sujet remontent à 1832, année à laquelle Michael Faraday tenta de mesurer le débit de la Tamise ([FAR32]). En effet, l’expérience était rendue possible grâce à la présence d’une grande quantité de sel dans l’eau qui interagissait avec le champ électromagnétique terrestre et produisait une différence de potentiel entre les deux rives du fleuve. L’intensité du courant à mesurer étant trop faible avec les moyens de l’époque, l’expérience se solda donc par un échec. En 1942, le terme « magnétohydrodynamique » fut pour la première fois introduit par Hannes Alfvén ([ALF42]). Sa contribution au domaine lui valut le prix Nobel de physique en 1970. Une application technologique de cette théorie est la torche à plasma (plasma wind tunnel). Historiquement, deux types de technologies différentes se sont développés à partir des années 1960 ([ABE00]). Les Européens et Américains se basèrent sur la génération d'un arc électrique entre deux électrodes pour augmenter la température du gaz jusqu'environ 10.000K par effet Joule ([AUW99]). Les Russes développèrent de leur côté le plasmatron, forme du procédé actuellement utilisée par les Européens ([GOR99], [REE61]) et également connue sous le nom de Inductively Coupled Plasma (ICP) dont le fonctionnement est expliqué dans le point suivant. A partir de 1970, des modèles de simulations des ICP hautes pressions apparurent peu à peu sous l'influence qu'eut le travail de Miller et Ayen ([MIL96]). Des contributions importantes furent celles de Boulos (1976), Mostaghimi (1985) et McKelliget (1986). 2. Physique du problème Le cadre du travail se repose entièrement sur la description suivante du problème. La von Kármán Institute (VKI) a mis en oeuvre une mini-torche à plasma ICP (Inductively Coupled Plasma, figure 1) dont les dimensions sont données sur la figure 2 ([ABE00]). Le principe de fonctionnement d’une telle torche est assez simple en soi : un gaz est injecté dans un tube en quartz autour duquel plusieurs spires sont placées. Ces spires sont excitées par un courant alternatif à très grande fréquence (de l’ordre du MHz), ce qui induit un courant secondaire au sein du gaz. L’apparition de ce courant secondaire va élever la température du gaz (~10.000K), par effet Joule, jusqu’à sa transformation partielle en un plasma faiblement ionisé. En raison de la très grande pureté du plasma produit, les Chapitre I - Histoire et physique du problème 1