modélisation numérique du chauffage par induction: approche éléments finis et calcul parallèle, Thèse de Mécanique

Télécharge modélisation numérique du chauffage par induction: approche éléments finis et calcul parallèle et plus Thèse au format PDF de Mécanique sur Docsity uniquement! HAL Id: tel-00443740 https://pastel.archives-ouvertes.fr/tel-00443740 Submitted on 4 Jan 2010 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. Modélisation numérique du chauffage par induction : approche éléments finis et calcul parallèle Valérie Labbé To cite this version: Valérie Labbé. Modélisation numérique du chauffage par induction : approche éléments finis et calcul parallèle. Mécanique [physics.med-ph]. École Nationale Supérieure des Mines de Paris, 2002. Français. ￿NNT : 2002ENMP1085￿. ￿tel-00443740￿ THESE présentée à L’ECOLE NATIONALE SUPERIEURE DES MINES DE PARIS par Valérie Labbé Pour obtenir le grade de Docteur de l’Ecole des Mines de Paris Spécialité Mécanique Numérique Modélisation numérique du chauffage par induction Approche éléments finis et calcul parallèle soutenue le 22 avril 2002, devant le jury composé de : Pr Michel BERNADOU ..............................Président Dr François-Xavier ROUX ...........................Rapporteur Pr Rachid TOUZANI...................................Rapporteur Dr François BAY..........................................Examinateur Dr Alain BOSSAVIT....................................Examinateur Pr Jean-Loup CHENOT...............................Examinateur Dr Serge PIPERNO.......................................Examinateur A mes parents… Ein Strudel von nie geahnter Seligkeit hat mich ergriffen… (Kleist) Index des figures 3 Index des figures Figure 1: Principe du chauffage par induction........................................................................ 12 Figure 2: Représentation de la profondeur de peau ................................................................ 12 Figure 3: Exemples d’inducteurs à une spire........................................................................... 13 Figure 4: Evolution de la distribution de puissance électromagnétique avec la température au passage de la température de Curie pour une billette d’acier de diamètre 120mm et une fréquence d’alimentation de 300Hz. ......... 19 Figure 5: Représentations des modes TE (a) et TM (b) pour une symétrie axiale. ................ 21 Figure 6: Partitionnement du domaine initial en un sous-domaine linéaire et un non- linéaire .............................................................................................................. 38 Figure 7: Conditions aux limites standards ............................................................................. 42 Figure 8: Conditions aux limites de notre modèle ................................................................... 43 Figure 9: Induction heating setup ............................................................................................ 47 Figure 10: Typical magnetization curve B(H) obtained with the Frohlich-Kenelly model. B is expressed in Tesla and H in A/m. Coefficients α and β are respectively 1.5 and 4500. ................................................................................ 57 Figure 11: Example of mesh for a simulation with a coil displacement: the area of displacement is delimited and meshed.............................................................. 71 Figure 12: Organisation of the coupling procedure ................................................................ 82 Figure 13: Static long inductor case: geometry and mesh (12000 nodes)............................... 84 Figure 15: Electrical conductivity, specific heat and thermal conductivity versus Temperature...................................................................................................... 84 Figure 16: Long inductor case: effective electric field isovalues ............................................ 86 Figure 18 : a) Evolution with respect to time of the experimental intensity in the coil ; b) Evolution with respect to time of the numerically computed electrical field in the part.................................................................................................. 87 Figure 20: Long inductor case: comparison between experimental and computed temperature evolutions ..................................................................................... 88 Figure 22: Long inductor case: Temperature field at two given time steps ............................ 89 Figure 23: Long inductor case: nodal velocities in the mesh at a given time step.................. 89 Figure 24: Static short inductor case: geometry and mesh ..................................................... 90 Figure 25 : Evolution of the electric field versus time for a non magnetic material at three given locations in the coil, the part and in the space in between. The electrical field is maximum in the coil....................................................... 91 Figure 26: Static short inductor case: effective electric field isovalues for two different permeabilities of the part at a frequency of 60 Hz. a) µrelative =1 b) µrelative =90.................................................................................................... 92 Figure 27: Effective electrical field profiles on the radial axis for 3 different frequencies f=700Hz, 1000Hz and 5000Hz. a)The part is a non magnetic material (relative magnetic permeability of the part equal 1), b). The part is a ferromagnetic material (relative magnetic permeability of the part equal 90). .................................................................................................................... 92 Figure 28: Static short inductor case at 60Hz: magnetic lines for rµ =1 and rµ =90........... 93 Figure 29: Temperature at two time steps at a frequency of 60 Hz for a non magnetic part (µrelative =1); a) t=5s b)t=25s................................................................... 93 4 Index des figures Figure 30: Temperature at two time steps for a magnetic part (µrelative =90); a) t=1s b)t=5s................................................................................................................ 94 Figure 31: nodal velocities in the workpiece due to dilatation effects for a frequency of 60 Hz; a) non magnetic material (µrelative =1), b) non magnetic material (µrelative =90)...................................................................................................... 94 Figure 32: Moving inductor case: geometry and initial mesh ................................................. 95 Figure 33: Moving inductor case: temperature isovalues at a given displacement step of the inductor....................................................................................................... 95 Figure 34: Géométrie et maillage ............................................................................................ 98 Figure 35: paramètres procédés .............................................................................................. 98 Figure 36: Données physiques des différents matériaux : σ, la conductivité électrique, µr la perméabilité magnétique, ρc la chaleur spécifique et k la conductivité thermique...................................................................................... 98 Figure 37 : Evolution temporelle du champ électrique en un point de la surface de la pièce. ................................................................................................................. 99 Figure 38: Evolution temporelle du champ électrique calculé à partir de deux initialisations différentes du schéma d’intégration de Lees. .......................... 100 Figure 39: Evolution temporelle du champ électrique calculé à partir de deux initialisations différentes du schéma d’intégration de Dupont....................... 101 Figure 40 : Evolution temporelle de la température en surface de la pièce pour différents schémas d’intégration en temps ..................................................... 102 Figure 41: Evolution temporelle du champ électrique en un point de la surface de la pièce pour différents pas de temps électromagnétique, T étant la période du signal électrique......................................................................................... 103 Figure 42 : Evolution temporelle du champ électrique en différents points du domaine pour une fréquence de 500Hz, un pas de temps dt=T/128, avec la période T=0.002s......................................................................................................... 104 Figure 43 : Modification du profil électrique avec le pas de temps à l’intérieur de la pièce ................................................................................................................ 105 Figure 44 : Profils radiaux de champ électrique calculés à partir de modèles considérant le terme volumique supplémentaire ou non ................................ 113 Figure 45 : Evolutions de la température en surface de la pièce calculées à partir de modèles considérant le terme volumique supplémentaire ou non .................. 114 Figure 46 : Variation du profil du champ électrique calculé avec le modèle complet pour différentes perméabilités magnétiques relatives .................................... 115 Figure 47 : Comparaison des profils du champ électrique effectif (code du Cemef) et de la norme du champ électrique complexe (code de LNMS). Les profils sont tracés dans la pièce de rayon 20mm............................................................... 116 Figure 48 : Evolutions temporelles de la température en surface de la pièce, sur l’axe de symétrie et entre les deux. a) solution analytique tirée de l’article de Wang [55] ; b) solution numérique ................................................................ 118 Figure 49 : cette photo montre la pièce avant traitement thermique (pièce de droite) et après (pièce de gauche) .................................................................................. 119 Figure 50 : Géométrie du levier de boîte de vitesse. Les zones en rouge représentent les zones subissant un traitement de surface. Les zones de déplacement de l’inducteur sont également représentées. ....................................................... 120 Figure 51: Iso valeurs du champ électrique effectif pour deux positions de l’inducteur....... 123 Figure 52: Iso-valeurs des températures au cours du temps. ................................................ 124 Figure 53: Architecture à mémoire partagée sur un calculateur à quatre processeurs........ 130 Index des figures 5 Figure 54: Architecture à mémoire distribuée sur un calculateur à quatre processeurs...... 130 Figure 55: Architecture à mémoire hiérarchique sur un calculateur à quatre processeurs ..................................................................................................... 131 Figure 56: Décomposition du domaine Ω en deux sous-domaines avec recouvrement ........ 140 Figure 57: Décomposition du domaine Ω en deux sous-domaines { }21 ,ΩΩ sans recouvrement .................................................................................................. 143 Figure 58: Presentation of the domain of study Ω and its boundaries. ................................. 155 Figure 59: Flow chart of the induction heating code............................................................. 161 Figure 60 : Partitioning method strategy............................................................................... 166 Figure 61: Exemple of distance value between two elements................................................ 167 Figure 62: Splitting of an initial multi-material mesh of 714 elements in three sub- meshes of 240, 236 and 238 elements............................................................. 168 Figure 63: The domain partitioning algorithm...................................................................... 170 Figure 64: General form of the parallel matrix-vector product ............................................ 171 Figure 65: geometry and process parameters........................................................................ 173 Figure 66: Effective electrical field (V.m-1) computations on 4 processors for a mesh of 2804 elements. ................................................................................................ 174 Figure 67: Efficiency of the parallel code with the parallel diagonal preconditioner with a mesh of 2804 nodes on the 32 bi-processors cluster. .......................... 174 Figure 68: Efficiency of the parallel code for 20 thermal iterations on a short inductor case on a mesh of 19966 nodes on the 32 bi-processors cluster. ................... 175 Figure 69: Efficiency of the parallel code for 20 thermal iterations on a short inductor case on a mesh of 29966 nodes on the 32 bi-processors cluster. ................... 175 Figure 70: Computational time versus the number of processors for three different mesh sizes ................................................................................................................. 176 Figure 71: Efficiency versus the number of processors for different mesh sizes ................... 176 Figure 72: Handling of the electro-thermal coupling ............................................................ 187 Figure 73: Algorithm of the direct model............................................................................... 187 Figure 74: General algorithm for the optimization procedure coupled to the induction heating model.................................................................................................. 191 Figure 75: The domain partitioning algorithm...................................................................... 194 Figure 76: General organization of the parallel optimization algorithm.............................. 196 Figure 77: Geometry of the gas bottle and of the ten inductors ............................................ 197 Figure 78: Evolutions of the current densities and of the cost function versus the number of iterations........................................................................................ 198 Figure 79: ISO curves of the effective electrical field on a partition into 8 sub-meshes from the initial mesh of size 70000 nodes....................................................... 199 Figure 80: ISO curves of temperature field on a partition into 8 sub-meshes from the initial mesh of size 70000 nodes. The colored regions are parts of the gas container. ........................................................................................................ 200 Figure 81: Efficiency of the parallel optimization code with a mesh of 19834 nodes on the 32 bi-processors cluster............................................................................ 201 Figure 82: Efficiency of the parallel optimization code with a mesh of 27784 nodes on the 32 bi-processors cluster............................................................................ 201 Figure 83 : Evolution of the CPU time with respect to the number of processors for a mesh of about 20000 and 28000 nodes........................................................... 202 8 Introduction electro-thermomécanique. Ce projet a impliqué le Cemef, le laboratoire Lnms de l’université de Ljubljana, le centre de recherche anglais EA-Technology, Transvalor ainsi qu’un forgeron anglais:UEF Chesterfield Cylinders, un groupe industriel italien SIAP-TQT et une PME danoise BL-Maskinfabrik. Ce projet a comporté quatre grands axes de recherche : − Modélisation directe des procédés de chauffage par induction afin de reproduire le plus fiablement possible les champs électriques et thermo- mécaniques à partir de paramètres procédés donnés, − Identification des paramètres physiques par analyse inverse : une modélisation réaliste nécessite en entrée des paramètres physiques de bonnes qualités. Des mesures fines de caractérisation physique ont été réalisées en Angleterre chez EA-Technology. Néanmoins, certains paramètres comme la perméabilité magnétique sont difficiles à mesurer, notamment pour des températures élevées. Ainsi une procédure d’identification des paramètres physiques par analyse inverse à été développée conjointement au Cemef et au Lnms, − Optimisation automatique des procédés : un algorithme d’optimisation automatique des procédés de chauffage par induction a été étudié et implémenté au Cemef. Il utilise une méthodologie de type contrôle optimal couplé à un algorithme de type gradient conjugué, − Calcul parallèle : le modèle direct et l’algorithme d’optimisation ont été parallélisés afin de réduire les temps de simulation ou permettre des calculs sur des maillages de grandes tailles. Les développements ont été effectués au Cemef. Les partenaires industriels ont apporté leurs expériences dans le domaine afin de mettre au point un modèle numérique adapté, qui réponde efficacement aux besoins des utilisateurs. Nous effectuons, en premier lieu, une analyse bibliographique des différents modèles mathématiques existants pour modéliser les phénomènes électromagnétiques avec les différentes approximations, plus ou moins fortes, couramment utilisées. Le domaine d’étude de ces phénomènes électromagnétiques étant par nature non borné, nous avons étudié les méthodes numériques employées et posé les choix stratégiques en vue d’obtenir un modèle performant. Nous décrivons dans le deuxième chapitre les modèles mathématiques électromagnétiques, thermiques et mécaniques que nous utilisons ainsi que leurs couplages. Des éléments de validation du modèle sont présentés. Introduction 9 Le troisième chapitre décrit la stratégie employée pour paralléliser le modèle direct. Il débute par une étude bibliographique des différentes méthodes existantes. La méthode de partitionnement de maillage ainsi que la méthode parallèle et son implémentation sont décrites. Des mesures de performance en termes de temps de temps de calcul sont présentées. Enfin la quatrième et dernière partie traite de la méthode de parallélisation de l’algorithme d’optimisation, présenté brièvement. Un cas d’optimisation des densités de courants pour des inducteurs fixes est présenté afin de servir de base à des mesures de performances en termes d’efficacité et d’accélérations par rapport à la version séquentielle. Au cours de cette thèse, trois articles ont été rédigés. Ils sont actuellement en cours d’évaluation par les comités de lecture. Nous avons inséré ces communications dans ce manuscrit dans l’état où elles étaient lors de l’envoi. Ainsi certaines parties sont rédigées en langue anglaise. Nous nous sommes efforcés d’insérer harmonieusement les publications dans le manuscrit. Chapitre 1: Etude bibliographique 13 Afin de transmettre la plus grande partie de l'énergie à la pièce à traiter, plusieurs paramètres sont à prendre en considération: - la disposition relative des inducteurs et de la pièce (couplage, longueurs respectives), - la fréquence d'alimentation et l'effet de peau qui caractérisent la répartition des courants induits dans la pièce : plus la fréquence augmente, plus les courants induits se concentrent en surface Cette notion fondamentale est déterminée par la profondeur de pénétration encore appelée épaisseur de peau. Typiquement, les inducteurs sont alimentés par des courants alternatifs de fréquence variant de quelques dizaines de Hertz à plusieurs centaines de milliers de Hertz, - les propriétés magnétiques (perméabilité relative), électriques (résistivité) et thermiques (conductibilité) des pièces à chauffer, variant pour la plupart avec la température, - le type d'inducteur (géométrie, nature du conducteur, technologie). Les géométries d’inducteurs peuvent être très variées, allant de la simple spire à des inducteurs multi-spires de formes complexes, Figure 3. Figure 3 : Exemples d’inducteurs à une spire 14 Chapitre 1: Etude bibliographique 1.2 Les applications industrielles Le procédé de chauffage par induction est de plus en plus utilisé et ceci de manière croissante dans les milieux industriels pour la préchauffe de pièces avant mise en forme à chaud (forgeage, matriçage, laminage, brasage), pour les traitements thermiques (trempe) ou encore pour des opérations de soudure entre pièces métalliques. Les traitements de surface recouvrent des opérations très diverses : - dégraissage, décapage, séchage, - galvanisation et étamage, - cuisson de vernis et peintures, plastification. Dans le cas de l'utilisation du chauffage par induction, le transfert thermique du revêtement s'opère du support vers l'extérieur, ce qui est favorable aux opérations de séchage et de cuisson (évacuation des solvants et vapeurs). Ce mode de chauffage permet donc d'obtenir : - une meilleure adhérence, - un meilleur aspect de surface, - une bonne reproductibilité, critère important pour le séchage des peintures colorées, - une grande souplesse d'utilisation par le choix des températures de traitement, - enfin, une ligne de production plus compacte et susceptible de fonctionner de façon discontinue, en l'absence de toute inertie thermique. Les applications dans ce domaine sont très vastes. Par exemple, on peut citer : - la polymérisation de vernis intérieur sur tubes aérosols, - la cuisson de joints d'étanchéité, - la polymérisation de vernis sur fils et méplats de cuivre, - le revêtement, - la ligne de galvanisation, - le recuit. Un autre type d’application qui tend à se développer récemment au sein des industries verrières, chimiques, céramiques, environnementales et chez les réfractoristes est la fusion de verre et d’oxydes par induction en autocreuset. En effet, la résistivité électrique des oxydes (1 à 10 Ω.cm à 1500°C) due à une conduction ionique est compatible avec la fusion par Chapitre 1: Etude bibliographique 15 induction. Leur faible conduction thermique aux basses températures et une résistivité décroissante avec la température permet d’utiliser la technique de l’induction directe en autocreuset avec une profondeur de peau égale au rayon de la charge. Cet autocreuset constitué du même matériau solide que l’on cherche à fondre, se forme grâce au refroidissement optimal de l’inducteur (mono-spire) et permet d’atteindre des températures supérieures à 2500°C sans contact du bain avec l’inducteur (pas de pollution du produit). Les applications sont les suivantes : - Fusion de cristal, - Fusion de verres spéciaux ou techniques, - Fusion d’oxydes réfractaires, - Elaboration de phosphates, - Vitrification de déchets. Quelle que soit la nature des applications industrielles, le chauffage par induction présente un certain nombre d'avantages intrinsèques qui expliquent son développement croissant : - rapidité de chauffage liée à la possibilité d'obtenir des densités de puissance très élevées, - localisation précise de l'effet thermique grâce à une conception d'inducteur et une fréquence de fonctionnement adaptée à la pièce à chauffer, - possibilité de chauffer à des températures très élevées avec un rendement pratiquement indépendant de la température. Ce procédé répondant parfaitement aux exigences industrielles de la moyenne et grande série : - facilité d'automatisation des équipements, - absence d'inertie thermique (démarrage rapide), - bonne reproductibilité des opérations effectuées, - rendement de chauffage souvent très élevé, - absence de pollution par la source de chaleur (source froide), - bonnes conditions de travail. 18 Chapitre 1: Etude bibliographique 2.1 Equations générales Les équations de Maxwell permettent de décrire tous phénomènes électromagnétiques. Elles sont au nombre de quatre et sont applicables sans aucune restriction à tous les milieux matériels : Équation du flux magnétique: 0=⋅∇ B , (2) Équation de Maxwell-Gauss: ( )libreliétotalE ρρρε +==⋅∇ )( , (3) Équation de Maxwell-Faraday: t B E ∂ ∂ −=×∇ , (4) Équation de Maxwell-Ampere: t E JH ∂ ∂ +=×∇ )(ε , (5) où B est l’induction magnétique, E est le champ électrique, H est le champ magnétique, J est la densité de courant électrique associée aux charges libres et totalρ est la densité de charge totale regroupant les charges liées et les charges libres. Les paramètres physiques sont µ , la perméabilité magnétique, et ε , la permittivité du milieu au point considéré. En particulier pour les métaux, les charges libres sont les électrons de conduction et les charges liées sont représentées par les cations du réseau cristallin. En régime permanent ou dans tout le domaine des fréquences hertziennes, on estime qu’il n’y a pas d’excédent local de charge et donc la densité totale de charge est considérée comme étant nulle. On peut alors réécrire l’équation de Maxwell-Gauss (3) : 0)( =⋅∇ Eε . (6) Pour des milieux isotropes, l’excitation magnétique H est relié à l’induction magnétique B par la relation constitutive: HB rr µ= . (7) Il est courant de décomposer la perméabilité magnétique par : Chapitre 1: Etude bibliographique 19 0µµµ r= , (8) où rµ est sans dimension et représente la perméabilité magnétique relative. Pour des milieux paramagnétiques et diamagnétiques, rµ est une constante très proche 1. En revanche, pour des milieux ferromagnétiques, la relation liant les champs B et H n’est plus linéaire : la perméabilité magnétique relative rµ est fonction de la norme de H et de la température T : HTHB r 0),( µµ= . (9) La dépendance de la perméabilité magnétique pour un ferromagnétique par rapport à la température est importante, notamment lors de la transition de Curie : le matériau devient amagnétique avec une perméabilité relative constante et proche de un. Les profils électromagnétiques dans la pièce vont être considérablement modifiés, voir la Figure 4. Figure 4 : Evolution de la distribution de puissance électromagnétique avec la température au passage de la température de Curie pour une billette d’acier de diamètre 120mm et une fréquence d’alimentation de 300Hz. Si on applique un champ électromagnétique de fréquence fixée aux frontières d’un matériau paramagnétique ou diamagnétique, donc linéaire, la réponse du milieu sera linéaire et les champs électromagnétiques internes au matériau oscilleront à la même fréquence bien que pouvant être déphasés. En revanche, pour des matériaux ferromagnétiques, des harmoniques secondaires de fréquence nouvelles sont générés déformant la forme de l’onde 20 Chapitre 1: Etude bibliographique électromagnétique : ces matériaux sont non linéaires. Cette non- linéarité se traduit mathématiquement par une dépendance de la perméabilité magnétique par rapport à H . La dernière équation nécessaire est la loi d’Ohm : EJ σ= , (10) où σ est la conductivité électrique dépendante de la température. Finalement le système d’équations initiales s’écrit: équation du flux magnétique 0=⋅∇ B , (11) équation de Maxwell-Gauss 0)( =⋅∇ Eε , (12) équation de Maxwell-Faraday t B E ∂ ∂ −=×∇ , (13) équation de Maxwell-Ampere t E JH ∂ ∂ +=×∇ ε , (14) relation intrinsèque au matériau HB rr µ= , (15) loi d’Ohm EJ σ= . (16) 2.2 Formulation mathématique tridimensionnelle La procédure standard consiste à remplacer le système d’équations différentielles du premier ordre (11)-(16) par une équation différentielle du second ordre de type équation de propagation des ondes. Par élimination de l’induction magnétique B dans l’équation de Maxwell-Faraday (13) à l’aide des équations (14) à (16) et en supposant le matériau isotrope, on obtient l’équation vectorielle suivante pour le champ électrique : Chapitre 1: Etude bibliographique 23 t E E r H z H zr ∂ ∂ += ∂ ∂ − ∂ ∂ θ θ εω , (27) t E E r rH r z z ∂ ∂ += ∂ ∂ εσθ )(1 . (28) Les équations (23), (25) et (27) sont découplées des équations (24), (26) et (28). D’autre part, il n’est pas nécessaire d’introduire les conditions de divergences nulles car elles sont déjà prises en compte indépendamment. Ainsi les champs [ zr HHE ,,θ ] sont indépendants des champs [ zr EEH ,,θ ] dans une configuration bidimensionnelle. Les systèmes d’équations peuvent donc être résolus indépendamment et leurs solutions ajoutées. L’ensemble de solutions [ ]zr HHE ,,θ constitue le mode TM, représentatif et adapté à un problème de chauffage par induction avec une symétrie axiale. Les champs [ ]zr EEH ,,θ représentent le mode TE et sont adaptés à la description d’un problème de conduction électrique ou chauffage par effet Joule, [13]. 2.3.2 Les inconnues du problème Pour une configuration bidimensionnelle à symétrie axiale, il est na turel de décrire le procédé de chauffage par induction avec le champ électrique θE en éliminant les champs zr HH , dans les équations (23), (25) et (27). Finalement, l’équation décrivant le champ électromagnétique peut être réduite à une équation scalaire pour le champ )0),,(,0(),( zrEzrE θ= r : 0 1 2 2 =×∇×∇+ ∂ ∂ + ∂ ∂ )E µ ( t E s t E e rr . (29) Néanmoins une formulation en champ électrique se rencontre rarement dans la littérature. La raison principale, mis à part des raisons historiques, provient certainement du fait que le terme source, par exemple un champ imposé ou une densité de courant source, n’apparaît pas explicitement dans l’équation (29). Le terme source doit être introduit manuellement en décomposant le champ électrique en une contribution par courant induit et une contribution par courant source imposé, ce qui permet d’introduire un terme supplémentaire. Un autre moyen de faire consiste simplement à imposer des conditions aux limites spécifiques sur la surface de la pièce ou de l’inducteur. 24 Chapitre 1: Etude bibliographique Plus communément dans la littérature, on rencontre plutôt le potentiel magnétique )0,,0( θAA = car dans ce cas le terme source est introduit de manière “naturelle” comme nous allons le voir au paragraphe 2.5.1 : sJVAt A t A =∇−=×∇×∇+ ∂ ∂ + ∂ ∂ σ µ σε ) 1 ( 2 2 r . (30) Dans une configuration plane avec une symétrie de translation, le mode TE est mieux adapté pour décrire les phénomènes électromagnétiques et dans ce cas l’inconnue du problème est réduite à la composante orthogonale du champ magnétique )),(,0,0( yxHH z= r . D’autre part, de par sa complexité, l’équation (29) ou (30) n’est jamais gardée telle quelle dans la littérature. Les approximations utilisées sont décrites dans la section suivante. 2.3.3 Les approximations standards a/ L’approximation des régimes quasi permanents (ARQP) Une approximation standard et couramment utilisée est de négliger les courants de déplacement (31) dans l’équation de Maxwell-Ampère [9], [48], [55]. Cette approximation des lois générales de l’électromagnétisme est valide pour des distributions ne variant pas trop rapidement dans le temps. t E J D ∂ ∂ = ε (31) Domaine de validité de l’ARQP Nous allons introduire les potentiels électriques et magnétiques obtenus à partir des équations de Maxwell (11) et (13) : - un champ vectoriel ),( trA rr appelé potentiel vecteur magnétique AB rrr ×∇= , (32) Chapitre 1: Etude bibliographique 25 - un potentiel scalaire électrique V tel que V t A E ∇− ∂ ∂ −= . (33) Pour plus de détail, voir les manuels de base sur l’électromagnétisme (par exemple [23]). En remplaçant dans les équations (3) et (5) les champs E r et B r par leurs expressions (32) et (33) et en utilisant la jauge de Lorenz : 0 1 2 vrr = ∂ ∂ +∇ t V c A. , (34) où c est la vitesse des ondes dans le vide, on arrive aux équations de Poisson : 0 1 02 2 2 vr r r =+ ∂ ∂ −∆ j t A c A µ , (35) 0 1 0 2 2 2 =+ ∂ ∂ −∆ ε ρ t V c V , (36) dont une solution physiquement acceptable pour une distribution de dimension finie est la solution des potentiels retardés : v 4 1 0 d r ) c r ?(t pe V(M,t) ∫ − = , (37) v 4 0 d r ) c r (tj p µ (M,t)A ∫ − = r r . (38) De même, dans le cadre de l’ARQP, on obtient de la même manière avec les équations réduites les solutions : ds r t tMV ∫= )( 4 1 ),( 0 ρ πε , (39) 28 Chapitre 1: Etude bibliographique ( )tierBRetrB ω)(~),( rr = , (49) ( )tierHRetrH ω)(~),( rr = , (50) ( )tierERetrE ω)(~),( rr = , (51) ( )tierARetrA ω)(~),( rr = , (52) où AEHB ~,~,~,~ représentent les champs électromagnétiques complexes et Re(.) la partie réelle. De cette manière, le problème est réduit à une équation stationnaire. Les modules des champs complexes calculés représentent les amplitudes efficaces des champs sinusoïdaux réelles et la phase des champs complexes donne la différence de phase avec le signal périodique source imposé aux bornes de l’inducteur. Cette approximation devient inadaptée si: - la source de courant appliquée aux bornes de l’inducteur n’est plus de forme sinusoïdale car dans ce cas les champs électromagnétiques eux-mêmes sont de nature différente, - le matériau employé a un comportement magnétique non linéaire car dans ce cas, une source de courant, même sinusoïdale, donne naissance à des champs électromagnétiques non sinusoïdaux. Néanmoins, malgré ces limitations importantes, la grande majorité des auteurs ont utilisé cette approche et ont favorisé, par une approche stationnaire, la réduction des temps de calcul. c/ Conclusion sur les approximations standards La quasi-totalité des modèles mathématiques utilisés pour modéliser le chauffage par induction négligent les courants de déplacement, approximation raisonnable dans le domaine de fréquence employé typiquement sur les installations industrielles. Une large majorité des auteurs a également choisi d’appliquer en plus l’approximation harmonique malgré le fait qu’elle soit mal adaptée aux matériaux ferromagnétiques. En effet pour ces matériaux, les harmoniques secondaires des champs électromagnétiques, de fréquences différentes, ne sont Chapitre 1: Etude bibliographique 29 pas calculées et donc n’apparaîtront pas au niveau du calcul de la puissance Joule injectée dans le calcul thermique. Bien sûr, pour des matériaux non magnétiques et une source sinusoïdale, cette approche est de loin la meilleure car elle amène des réductions significatives en terme de temps de calcul. 2.4 Les modèles électromagnétiques standards en symétrie axiale Les différentes équations que l’on peut trouver dans la littérature pour décrire le champ électromagnétique au cours d’un procédé de chauffage par induction axisymétrique ont été répertoriées. Tous les modèles sont fondés sur le calcul du potentiel vecteur magnétique réduit dans ce cas à une composante scalaire perpendiculaire au domaine d’étude : )0),,(,0( zrAA θ= . En remplaçant B par son expression ArotB r = dans l’équation (44) et en utilisant la loi d’Ohm (46) on obtient : EJA σ µ ==×∇×∇ ) 1 ( . (53) La densité de courant J représente la contribution des courants sources, aussi bien que la contribution des courants induits. Il est important de noter que des courants sont induits aussi bien dans la pièce que dans l’inducteur. Il suffit alors de remplacer E r par son expression en terme de potentiel (33) pour obtenir l’équation vectorielle : sJVAt A =∇−=×∇×∇+ ∂ ∂ σ µ σ ) 1 ( , (54) où VJ S ∇−= σ représente la densité de courant imposée dans l’inducteur. En coordonnées cylindriques, l’équation (54) devient : ( ) sJz A z rA rrrt A =      ∂ ∂ ∂ ∂ −      ∂ ∂ ∂ ∂ − ∂ ∂ θ θ θ µµ σ 11 . (55) Au lieu de développer l’équation (54) en coordonnées cylindriques, de nombreux auteurs [17], [34], [41], [48], [55] utilisent la formule vectorielle : ) 1 () 1 () 1 ( AAA ∇⋅∇−⋅∇∇=×∇×∇ µµµ . (56) 30 Chapitre 1: Etude bibliographique Associée à la condition de jauge de Coulomb 0=⋅∇ A , (54) devient une équation de type diffusion : sJAt A =∇⋅∇− ∂ ∂ ) 1 ( µ σ . (57) Si on réécrit l’équation (57) en coordonnées cylindriques et que l’on compare avec l’équation (55), on trouve qu’il manque deux termes dans l’équation (57) : 2r A µ θ , (58) r A r θ µ        ∂ ∂ − 1 . (59) Le terme (59) peut être significatif pour des matériaux ferromagnétiques sous la température de Curie ou en présence de gradients de température locaux. Cependant ce terme est toujours négligé dans les modèles mathématiques pouvant être trouvé dans la littérature. De même, peu d’auteurs prennent en compte le terme supplémentaire (58), [36]. Il serait intéressant d’évaluer l’importance de ce terme dans le calcul du potentiel magnétique. L’approximation harmonique L’inducteur est parcourue par un courant sinusoïdal de pulsation ω . La formulation harmonique consiste à remplacer le potentiel magnétique par son expression complexe (52) dans l’équation électromagnétique (57). La dérivée en temps du potentiel magnétique est remplacée par le produit du potentiel par ωj . On obtient les équations stationnaires suivantes où le terme supplémentaire (58) est présent ou non suivant les auteurs [36] : S J)A µ (Ajs rrrr =∇⋅∇− 1 ω , (60) sJr A AAj =+∇⋅∇− 2 ) 1 ( µµ σω . (61) Une équation similaire peut être développée pour le champ électrique réduit à sa composante orthoradiale scalaire. Chapitre 1: Etude bibliographique 33 3.1.1 Formulation variationnelle Deux procédures distinctes permettent d’obtenir une formulation faible pour le problème électromagnétique. La première utilise la méthode des résidus pondérés, la deuxième est basée sur la minimisation d’une fonction variationnelle. Nous allons montrer l’application de la formulation faible sur l’équation dépendante du temps (61) qui comprend le terme supplémentaire (58). Les autres approximations, harmoniques ou sans le terme (58) en sont uniquement des cas particuliers. Soit l’équation générale: S J r A A t A =+           ∇⋅∇− ∂ ∂ 2 1 µµ σ θ θ θ , (70) où le Laplacien s’écrit en coordonnés cylindriques: ) 1 ( θµ A∇⋅∇ = ) 1 ()( 1 θθ µµ A zz A r r rr ∂ ∂ ⋅ ∂ ∂ + ∂ ∂ ⋅ ∂ ∂ . On définit l’espace V par :    = ∂ ∂ Ω∈Ω∈    = 0 v ),(L r v ),(HvV 21 θ . La formulation faible de l’équation (70) s’écrit, V∈∀w : v v v 1 v 2 dwJdwr A dAwdw t A S∫∫∫∫ ΩΩΩΩ =+      ∇⋅∇− ∂ ∂ µµ σ θ θ θ , (71) où l’élément de volume dv s’écrit dv=2πrdrdz. En intégrant par partie le second terme de l’équation (71) sur le domaine Ω , nous obtenons la formulation faible du problème électromagnétique [34], [55]: ∫∫∫∫ ∂ ∂ +=∇∇++ ∂ ∂ G0O 2 ds µ 1 dv dv 1 dv w n A wJAw. µ ) w µr A t A (s ? S O ? O ?? , (72) 34 Chapitre 1: Etude bibliographique où Γ est la frontière du domaine fermé Ω, 0µ est la perméabilité magnétique sur la frontière extérieure Γ , soit celle de l’air. Le dernier terme dans le second membre disparaît si des conditions aux limites de Neumann nulles ou de Dirichlet sont appliquées. Il est intéressant de montrer que dans le cadre de l’approximation harmonique, la fonctionnelle I à minimiser est complexe : ∫∫∫ ΓΩΩ ∂ ∂ + ∂ ∂ −+− ∇ +∇= ds n A A n A AdsAJAJds r A AI ss )( 1 )()( 1 * * 0 ** 2 2 2 θ θ θ θθθ θ θ µµ . (73) Cette fonctionnelle doit être minimisée relativement au complexe conjugué *θA pour obtenir l’équation régissant le potentiel vecteur magnétique, [38]. 3.1.2 Formulation des équations intégrales Une autre alternative à l’approche variationnelle consiste à formuler le problème électromagnétique en termes d’intégrales sur des sous-domaines fermés du domaine d’étude ou bien sur les frontières de ces sous-domaines. Il existe principalement deux méthodes intégrales pour notre problème électromagnétique. La plus répandue dans la littérature s’appuie sur l’application de la seconde relation de Green et conduit à une formulation intégrales sur les contours frontières, [3], [4]. La deuxième méthode consiste à calculer le champ magnétique à l’aide de la loi intégrale de Biot et Savart. a/ Formulation intégrale aux frontières – Relation de Green Un problème électromagnétique où les phénomènes physiques peuvent être décrits par un potentiel scalaire φ sont généralement gouvernés par une équation de type Laplace (74) ou par une équation de type Helmholtz (75) : Ts=∇ φ2 , (74) Tsfßf 22 =+∇ , (75) où Ts représente le terme source imposé. Pour la modélisation des procédés de chauffage par induction, l’équation de Laplace gouverne généralement le comportement du potentiel φ dans l’air alors que l’équation de Helmholtz décrit le comportement du potentiel φ avec une Chapitre 1: Etude bibliographique 35 dépendance en temps harmonique dans des matériaux linéaires avec une conductivité électrique )(Tσ et une perméabilité magnétique )(Tµ constantes ou dépendantes de la température. Une formulation intégrale aux frontières est obtenue en introduisant le noyau de Green et en utilisant le théorème de Green. L’équation de Helmholtz dans le cadre de l’approximation harmonique pour le potentiel magnétique et pour un matériau linéaire s’écrit : s?2 2 ? 2 Jµ)A r 1 (aA −=+−∇ , (76) où le nombre complexe ωσµα j=2 , µ et σ dépendent uniquement de la température (matériau linéaire) et SJ est la densité de courant source dans l’inducteur ( 0=SJ hors de l’inducteur). Nous introduisons la fonction de Green z)(r,G ax qui représentent le potentiel vecteur généré par une ligne de courant localisée dans le plan d’étude au point ),( 00 zr . Cette fonction est définie par l’équation suivante : ),() 1 ( 002 22 zzrrG r G axax −−−=+−∇ δα . (77) La solution analytique pour la fonction de Green en configuration axisymétrique s’écrit :       −−= )()() 2 1( 1 ),,,( 2 0 00 mEmK k r r k zzrrGax π , (78) avec 2 0 2 0 02 )()( 4 zzrr rr k −+− = , (79) 21 km −= , (80) où )(mK et )(mE représentent les intégrales elliptiques du premier et second degré. Les valeurs de la fonction de Green et de son gradient dépendent uniquement des paramètres géométriques et leurs calculs sont rapides. En appliquant la seconde identité de Green (81) à l’équation (76), on arrive à l’équation (82) : ∫∫ ΩΩ ∂ ∂ − ∂ ∂ =∆−∆ δ ds n A G n G AdvAGGA )()( , (81) 38 Chapitre 1: Etude bibliographique ∫ Ω = rdsNJS jj 0 . (88) Dans le cadre d’une approximation harmonique, le système discrétisé devient : [ ] [ ] [ ]SAK = , (89) ∫ Ω ∇∇++= rdsNNNN r jK jijiij ) 1 ) 1 (( 2 µ σ µ ω . (90) Une discrétisation spatiale intéressante a été proposée par Marchand et Foggia [35], basée sur une partition puis un couplage entre les domaines linéaires et non- linéaires. Le domaine global Ω est divisé en un domaine linéaire comprenant l’air et le ou les inducteur(s) et un domaine constitué de la pièce. Ce dernier sera non linéaire pour des métaux ferromagnétiques. Ces deux sous-domaines sont résolus séparément de manière indépendante. Le couplage est réalisé par une condition de continuité de la dérivée (91) sur l’interface commune aux sous- domaines. n rA ∂ ∂ )( θ (91) Figure 6: Partitionnement du domaine initial en un sous-domaine linéaire et un non-linéaire La méthode des éléments finis permet de modéliser le problème électromagnétique de manière précise notamment par une bonne prise en compte du comportement non linéaire du champ électromagnétique dans la pièce si le matériau est ferromagnétique. Une limitation de la méthode provient du fait que le domaine naturel d’étude est infini. Fermer le domaine par une frontière arbitraire peut conduire à l’apparition d’ondes réfléchies sur les frontières ce qui va modifier le résultat du calcul. Pour que cette méthode reste fiable, Non linéaire Pièce Inducteur Domaine linéaire Non-linéaire Interface commune Chapitre 1: Etude bibliographique 39 il faut que le domaine soit pris “assez grand”. D’autre part, étant donné que tout le domaine -y compris l’air- est maillé, le nombre d’éléments nécessaires peut facilement croître, ce qui va entraîner des besoins en stockage et en temps de calcul importants. De plus, pour modéliser le déplacement de l’inducteur, de nombreux remaillages vont être nécessaires. 3.2.2 Méthode des éléments frontières La discrétisation des équations intégrales de frontière utilise les techniques classiques de discrétisation. Les domaines pièce et inducteur uniquement sont discrétisés. Le potentiel PAθ au nœud ),( ii zrP = du contour c est la somme des intégrales de chaque élément frontière : dlr n G AdlrG n A rdsGJAr fC ax NelC f k fnbnoC kfC ax kNelC f fCnbno k Nel e nbnoe m e ax m s i i ∫∑ ∑∫∑ ∑∑ ∑ ∫ ∂ ∂ − ∂ ∂ += == Ω = Ω = Ω 111 12 1 θ θ θ µ , (92) où ΩNel est le nombre d’éléments du domaine Ω , CNel le nombre d’éléments de la frontière c du domaine Ω , Ωnbnoe le nombre de nœuds du domaine Ω et CeNbno le nombre de nœuds du contour c. La difficulté avec cette méthode provient des intégrales singulières lorsque le point p où le potentiel est calculé appartient à l’élément frontière où l’on cherche à évaluer l’intégrale. En effet la fonction de Green n’est alors pas définie. Cette difficulté est contournée en calculant numériquement ces intégrales à l’aide de fonctions logarithmiques pondérées [3]. La présence de singularités géométriques, tels des angles, complique également substantiellement l’évaluation des intégrales, [29]. Finalement, on arrive au système suivant : [ ] { } [ ] { }T n A NAM =       ∂ ∂ + , (93) avec ∫ ∂ ∂ += jC ax ijiij rdln G rM δ 2 1 , (94) dlrGN jC axij ∫−= , (95) 40 Chapitre 1: Etude bibliographique       ∈ = ∫ sinon. 0 inducteur,O si rdsG Jµ T j jO axs ij (96) La méthode des éléments de frontière présente de nombreux avantages : − prise en compte naturelle du domaine ouvert infini, − seuls les contours du domaine sont discrétisés. Les inconnues du problème sont uniquement calculées sur les frontières. Ainsi on passe d’un problème à 2n inconnues à un problème à n inconnues, − le potentiel en chaque point interne du domaine est ensuite calculé rapidement par simple intégration numérique, − prise en compte aisée du déplacement des inducteurs. Néanmoins, les désavantages sont loin d’être négligeable: − les matériaux magnétiques ne sont pas pris en compte car seul les milieux linéaires peuvent être traités, − la matrice générée par la discrétisation spatiale est pleine, − il faut calculer des intégrales singulières. On voit que l’inconvénient principal provient de l’incapacité de la méthode à traiter les matériaux non linéaires. C’est ce point essentiel qui a motivé l’intérêt pour les méthodes mixtes. 3.2.3 Les méthodes mixtes a/ Eléments finis / éléments frontières Cette formulation regroupe une formulation variationnelle sur les domaines fermés, représentant les matériaux métalliques, alliée à une formulation intégrale du champ électromagnétique dans l’air. Les principes généraux sont décrits notamment par Salon [46], [47], Zienkiewicz [57] ou Brebbia [3]. La méthode des éléments de frontière permet de Chapitre 1: Etude bibliographique 43 Figure 8 : Conditions aux limites de notre modèle 3.4 Discussion Nous avons vu qu’il existe différents choix de variables d’état et d’approximations possibles pour la résolution du modèle électromagnétique. Il est inutile, dans la gamme de fréquence utilisée par le procédé, de résoudre les équations de Maxwell complètes, trop complexes. Nous avons décidé de prendre en compte l’approximation des régimes quasi-permanents qui permet de passer d’une équation de type propagation des ondes à une équation de type diffusion. Cette approximation est quasiment toujours présente dans tous les modèles rencontrés dans la littérature. En revanche, l’approximation harmonique se base sur des hypothèses fortes étant donné qu’elle n’est réellement valide que pour des matériaux linéaires et donc non ferromagnétiques. Elle est cependant largement utilisée par de nombreux auteurs car elle apporte un gain de temps de calcul très important, du fait que l’on obtient au final une équation stationnaire complexe modélisant le comportement de l’amplitude du champ électromagnétique. Nous avons choisi de considérer une équation dépendante en temps, beaucoup plus riche, qui prend en compte tous les phénomènes non linéaires et notamment la génération d’harmoniques secondaires pour les ferromagnétiques. Cependant, ce modèle est plus complexe à résoudre. Au niveau du choix de la variable d’état, nous avons vu que le champ électrique et le potentiel magnétique sont réduits à une unique composante orthoradiale. Les formulations avec l’une ou l’autre variable sont quasiment équivalentes. Néanmoins, nous n’avons pas vu de modèle résolu en champ électrique dans la littérature, ceci certainement pour des raisons historiques. Les principales différences entre les deux formulations se résument premièrement à la prise en compte de la non- linéarité dans l’équation, due à la dépendance de la perméabilité magnétique avec le champ magnétique, lui même calculé (de manière différente) à partir du champ 0= ∂ ∂ z rAθ Conditions de symétrie 0=θA 0=+ ∂ ∂ θ θ A r rA Condition aux limites de type Robin : 44 Chapitre 1: Etude bibliographique électrique ou du potentiel magnétique. Une deuxième différence réside dans la prise en compte de la densité de courant source. Cette dernière est soit introduite directement dans l’équation en potentiel magnétique, soit prise en compte à travers sa dérivée pour l’équation en champ électrique. Enfin, le terme source pour l’équation de la chaleur (issu de la puissance joule) est soit obtenu directement dans le cas d’une résolution en champ électrique (car proportionnel au carré du champ électrique), soit dérivé numériquement par rapport au temps dans le cas d’une résolution en potentiel magnétique. Dans l’un ou l’autre cas, nous avons des imprécisions dues aux dérivées temporelles numériques. Cependant, étant donné que la densité de courant source est souvent introduite de manière analytique (par exemple sous la forme d’une sinusoïde), sa dérivée va être exacte, ce qui favorise une résolution en champ électrique. D’autre part, nous trouvons qu’elle apporte une représentation physique plus intuitive. Concernant le choix de la méthode numérique, nous avons vu que la méthode éléments finis était nécessaire pour traiter des matériaux non linéaires. D’autre part, le choix du modèle électromagnétique tenant compte des non- linéarités est rédhibitoire pour une méthode basée uniquement sur les intégrales de frontière. Les méthodes mixtes possèdent des qualités certaines avec une bonne prise en compte du domaine ouvert, des matériaux non linéaires et un déplacement aisé de l’inducteur. Néanmoins, la résolution de la matrice du système, pleine, peut s’avérer extrêmement coûteuse suivant la taille des maillages traités. D’autre part, dans une optique de calcul parallèle, indispensable pour la simulation de cas réels conséquents, la parallélisation d’une méthode mixte va être difficile à mettre en place. Pour ces différentes raisons, nous nous sommes orientés vers une méthode numérique “tout éléments finis”. Chapitre 2 : Modélisation du chauffage par induction 45 Chapitre 2 Modélisation mathématique et numérique des procédés de chauffage par induction en configuration axisymétrique Nous présentons dans ce chapitre la stratégie de modélisation des procédés de chauffage par induction développée au centre de Mise en Forme des Matériaux. Les modèles physiques, électromagnétique, thermique et mécanique sont décrits. Puis nous verrons comment les systèmes d’équations aux dérivées partielles issus de ces modèles sont discrétisés d’abord spatialement puis temporellement. Les phénomènes électromagnétiques et thermiques présentent des échelles de temps fondamentalement différentes. Le couplage multi-physique va devoir tenir compte de ces caractéristiques propres. Ces différents points font l’objet de la publication insérée ci après qui va également présenter les premiers éléments de validation du logiciel face à des publications ou à des mesures expérimentales effectuées sur site industriel. Suite à la publication insérée, des études complémentaires sur le choix du schéma d’intégration optimal sont décrites. Enfin le mode de stockage utilisé ainsi que le solveur utilisé pour résoudre les systèmes matriciels sont détaillés. 1 Le modèle mathématique et sa résolution numérique Le modèle mathématique du procédé de chauffage par induction et les méthodes de résolution numérique sont présentés dans l’article suivant, soumis à l’International Journal for Numerical Methods in Engineering (IJNME) 48 Chapitre 2 : Modélisation du chauffage par induction 1.2 Potential applications of a numerical model Industrial use of induction heating processes aims at achieving various objectives. We can name, among others: - reaching a uniform prescribed temperature, - achieving control of the grain size, - getting a certain hardness level at the surface, - …. Industrials try to reach most of these goals by trial-and-error procedures. But these procedures cannot most of the time achieve fine and accurate control of these goals and are furthermore time consuming. It is clear that, since induction heating processes involves and couples various physical phenomena, the goals will be achieved only through an in-depth understanding of the industrial process, and therefore through numerical modeling. The numerical model presented here has been developed having in mind the fact that a fine and efficient numerical model coupled with optimization techniques [21] could be a major step in determining optimal process parameters (frequency, current density, coil geometry, …) suitable to reach the industrial objectives previously mentioned. The frequency is an important control parameter to achieve either surface heating (high frequencies) or a rather uniform heating in the workpiece with a skin depth of about the same order of magnitude or slightly larger than the geometric depth (low frequencies). 1.3 The various numerical models Several numerical models have been developed up to this day in order to model the induction heating process. These models carry out a simulation of the induction heating process but with various restrictions. Some models do not carry out full coupling between electromagnetism, heat transfer and solid mechanics, but focus only on the electromagnetism/heat transfer coupling. Most numerical models use a harmonic approximation to evaluate the electromagnetic field in the workpiece. This approximation is valid for linear magnetic materials but can become inaccurate when dealing with non- linear magnetic materials [10], [44]. Coupling between inductor and workpiece is carried either through the use of a boundary element method [47], [10] or through finite elements [43]. In this last method, the air domain is also meshed and the electromagnetic equation is solved on the global domain made out of the air, the inductor and the workpiece. Chapitre 2 : Modélisation du chauffage par induction 49 2. The mathematical model As we have stated previously, induction heating processes involve three main different physical phenomena, which are related to: - electromagnetism, - heat transfer, - solid mechanics. A complete mathematical model for induction heating processes therefore needs to consider and tightly couple the equations modeling these phenomena. We present here a general mathematical model which couples: - the Maxwell equations - in order to access the electric field in the workpiece, - the heat transfer equation - the heat dissipated in the part by the eddy currents being the main source of temperature evolution in the workpiece, - the mechanical equilibrium equations - which will enable us to access the stress and strain rates in the workpiece. 2.1. The electromagnetic model 2.1.1. The Maxwell equations The global system of equations modeling electromagnetic wave propagation and any kind of electromagnetic phenomena was established by Maxwell more than a century ago. This system is made up of the four following equations: Magnetic flux equation: 0=∇ B. rr , (99) Maxwell-Gauss equation: 50 Chapitre 2 : Modélisation du chauffage par induction 0=∇ E. rr , (100) Maxwell-Faraday equation: t B E ∂ ∂ −=×∇ r rr , (101) Maxwell-Ampere equation: t D jH ∂ ∂ +=×∇ r rrr , (102) where H is the magnetic field, B the magnetic induction, E the electric field, D the electric flux density, j r the electric current density associated with free charges, and x denotes the vector product - thus E rr ×∇ is the curl vector of E . This system of equations needs to be completed by relations which take into account material properties (magnetic permeability µ, electrical conductivity σ, dielectric constant ε). These relations are the three following ones: ED rr ε= , (103) HHTB rrr ),(µ= , (104)      = air. in the 0 ,conductors in the r r r Es(T) j , (105) relation (105) being commonly known as the Ohm law. It should be pointed out here that the magnetic permeability µ, and the electrical conductivity σ do strongly depend on temperature. Furthermore, for non- linear ferromagnetic materials, the magnetic permeability µ may also depend on the magnetic field strength H r . Equations , (99) - (105) are valid for any range of frequencies and for any kind of geometric configuration. Chapitre 2 : Modélisation du chauffage par induction 53 t S J EvE t E ∂ ∂ −=×∇×−      ×∇×∇+ ∂ ∂ r rrrr r 1 σ µ σ , (116) or in cylindrical coordinates: t S J ? E z . z sv r ? E µrr ? E µ? E µ . t ? E s 1 2 11 ∂ ∂ −= ∂ ∂ −      ∂ ∂ −+      ∇∇− ∂ ∂ . (117) As the inductor is moving in parallel with the axis of symmetry, its velocity writes )),(,0,0(),( trvtrv Z rrr = . Indeed, a component in the axial direction of the velocity would imply a transformation of the inductor size and geometry, which is not meaningful in industrial induction heating applications. If tc is the characteristic time of the electromagnetic phenomena sTtC 001.032/ <<≈ for f = 50Hz, E the characteristic electrical field value, CL is the characteristic length, 0µµµ r= where rµ is the non-dimensional relative permeability and 0µ is the magnetic permeability of the air. By replacing in relation (117), CC LrrtttEEE ~,~, ~ ===θ where rtE ~,~,~ are non dimensional, we get the non-dimensional equation: t J ct J z E cL zvE r E rrcL E r E rcL EE rcL E t E ct E ~ ~ ~ ~ . ~ ~1 ~2 0 2~ ~1 2 0 ~~1.~ 2 0 ~ ~ ∂ ∂−= ∂ ∂ −         ∂ ∂−+        ∇∇− ∂ ∂ σ µµµµµµ σ (118) which becomes t J E J z E cL z v c t r E rrcL ct r E rcL ctE rcL ct t E ~ ~ ~ ~ . ~ ~1 ~2 0 2~ ~1 2 0 ~~1.~ 2 0 ~ ~ ∂ ∂ −= ∂ ∂ −         ∂ ∂ −+         ∇∇− ∂ ∂ σ µσµµσµµσµ (119) 54 Chapitre 2 : Modélisation du chauffage par induction If we compare the two non dimensional numbers 2 0 cL ct σµ and cL zvct , (120) we have: zv cL >> σµ0 1 (121) as 1.05.0 −≤ smv z and mLC 1≈ typically. Hence the term due to the coil velocity may be neglected in a first approximation. 2.1.5. The source term The translation of experimental measurements of the power or current intensity into numerical input data is one of the most complex tasks in the simulation of induction heating processes. There are basically two ways to enter the input for the electromagnetic equation – either by entering a current density or a tension. a) The current density source term The electromagnetic equation with )0),,,(,0(),,( tzrJtzrJ SS θ= r writes: t S? J r ? E µrr ? E µ? E µ . t ? E s ∂ ∂ −=      ∂ ∂ −+      ∇∇− ∂ ∂ 1 2 11 . (122) The variation with time of the current intensity may be of any shape (sinusoidal or more complex shapes). For a sinusoidal source current, we get: Chapitre 2 : Modélisation du chauffage par induction 55 )cos(),(),,( ϕ+= wtzrJtzrJ OS , (123) with the pulsation fπω 2= . b) The tension source term At a local point, the total current density writes V∇− ∂ ∂ −= r r r s(r,z) t A s(r,z)(r,z,t)j , (124) where tAjinduced ∂∂−= rr σ denotes the induced current density, V∇−= rr σsJ the imposed current density and V the electric scalar potential. As we do not know the local potential value of the gradient V∇ r , we rather use a uniform gradient over the coil section. This uniform gradient equals the total voltage applied at both ends of the coil divided by the total length of the coil: turns_of_number*r2 V L V J coil tot coil tot S π σ −=σ−= . (125) 2.1.6. The boundary conditions The choice of boundary conditions when carrying out a global finite element simulation needs to be carried out carefully. Finite element computation of electromagnetic field in the air can be altered by “artificial reflections” on the boundary of the air domain. Hence the outer boundaries need to be far enough from the induction heating set-up to reduce the approximations. A fairly good approximation can be obtained if a Robin boundary condition is applied on an enclosing “box” for the surround ing air. The dimensions of the “box” being “large enough” in order not to perturb the magnetic field lines. This Robin like condition writes in our case: 58 Chapitre 2 : Modélisation du chauffage par induction This relation works well when one knows the magnetic field H r and needs simply to compute the magnetic permeability ( )TH ,rµ or the magnetic induction field Br . Conversely, to compute the H r field one needs to have an equivalent relation: ( ),TBµ B H r r r = . (129) This relation is obtained from the first order Frohlich-Kenelly relation (130) by seeking the physically realistic solution of the second order equation (131) in H r : ( ) ( )HTH r r + += β α µµ 0, , (130) ( ) HHB r r r         + += β α µ 0 . (131) Hence we get: ( )       +−++−−= BBBH rrrr βµβµαβµα µ 0 2 00 0 4 2 1 . (132) 2.1.7.2 Formulation in H and T in the time domain decomposition The Frohlich-Kenelly formulation is still a good approximation when temperature dependency is added. In this case, coefficients α and β are no longer constant, they are temperature-dependent, [5]. In order to use an easily implementable formulation, one would rather separate the temperature dependency from the magnetic field dependency. This dependency separation has been used by various authors. Listing [10], [49] and [27] giving the following general form: ( ) ( ) ( )TfH H H TH ., 20           + − + += r r r r β α β α µµ . (133) Chapitre 2 : Modélisation du chauffage par induction 59 Several formulations f(T) can be found in literature. We have chosen to use the following formulation which has been used for instance by [2]: ( ) ? c T T Tf         −= 1 , (134) where CT is the Curie temperature and γ is the temperature sensibility parameter. Relation (132) becomes with the temperature dependency: ( )       +−++−−= BBffBH rrrr βµβµαβµα µ 0 2 00 0 4 2 1 . (135) 2.1.7.3 Formulation in the heat transfer cycle As will be explained in the coupling strategy paragraph, at the end of each heat transfer time- step solving, differences between the electromagnetic parameters used to compute the electric field and the ones calculated using the newly calculated temperature maps are evaluated. If these differences exceed a given threshold, a new cycle of electromagnetic computations is performed; if not, the heat transfer solving goes on. Indeed, the magnetic field is no more evaluated at each electromagnetic time step. The magnetic permeability has therefore to be defined as in harmonic models. A few authors dealing with harmonic models take into account the non- linearity of the magnetization curve. [24] gives a wide bibliographical research on the best way to handle this non- linearity. [12] and [52] use a mean magnetic permeability that can be defined as the harmonic mean value over one period. The main problem with the time-averaged value is the lack of physical meaning. To overcome this lack, [6] and [54] calculate an equivalent magnetization curve using some energy equivalency methods. Using energy equivalency methods, assuming that the magnetic field is sinusoidal, and searching a linear equivalent magnetic permeability that induces, over one period, the same magnetic energy yields to: ( ) ∫ ∫= pf tB m eq dtHdB Hfp µ 2 0 0 2 2 4 , (136) where f is the electromagnetic frequency and Hm the maximum value of the sinusoidal magnetic field. Since we use a time domain decomposition, the electromagnetic fields may be non-sinusoidal. Using Frohlich-Kenelly relation and the same approach than in [32] - except 60 Chapitre 2 : Modélisation du chauffage par induction that time is no longer taken into account, we have found that the equivalent magnetic permeability becomes, with respect to temperature and the maximal magnetic field : ( )               −            +β βαβ + α +µ=µ γ C 2 maxmax 0maxlin T T 1 H ln H 2 H 2 T,H . (137) 2.1.8 Computation of the magnetic field H r The magnetic field values need to be determined when dealing with the heating of non- linear magnetic materials, since the magnetic permeability depends in that case on the strength of the magnetic field. The magnetic induction field is computed from the electric field values by using the relation t B E ∂ ∂ −=×∇ r rr . (138) Using an explicit time integration scheme, the magnetic induction field at time t can be derived from the one at time t-δt writes:                    ∂ ∂ −+−= ∂ ∂ +−= ∂ ∂ −−= r E r E t?t)?(tB)E(r. rr ?t t)?(tB(t)B E z .?tt)?(tB(t)B t)z,(r,B ?? z?zz ?rr r (139) The magnetic field H r is then computed from relation (132) or (135). 2.2 The thermal model Modeling of induction heating processes involves modeling of heat transfer - which is mainly governed by the heat dissipated by eddy currents in the workpiece - and heat fluxes at the interface between the workpiece and the air. The model will have therefore to include the heat transfer equation, the heat source term due to eddy currents, as well as appropriate boundary conditions. Chapitre 2 : Modélisation du chauffage par induction 63 The classical equilibrium equations for a solid undergoing mechanical loads can be expressed in a local form by: 0 rr =+ )f?div(s , (145) where σ denotes the stress field in the solid and f r the body forces per mass unit. It needs to be noted that the material behavior is not taken into account at this stage. 2.3.2 The constitutive law The material behavior is assumed to obey a thermo-elastic-viscoplastic constitutive law. The material strain rate ε& is decomposed as the sum of three fractions – namely, the elastic fraction elε& , the plastic fraction plε& , and the thermal one thε& : htplel ε+ε+ε=ε &&&& . (146) The elastic fraction elε& is related to the stress field tensor σ& through the Hooke law : ij 1 ds E ? s E ? e kkij el ij &&& − + = , (147) where E and ν respectively denote the Young modulus and the Poisson coefficient. It should be noted here that we do not need to use a Jaumann derivation for σ& since the rotations involved are quite small. The deformation thε& due to thermal expansion can be expressed as: ij th ij dTae && = , (148) where α denotes the volume thermal expansion coefficient. The equivalent stress is defined through the Von Mises yield criterion: . 2223 222 2 12 ) xz s yz s xy (s zz s xx s zz s yy s yy s xx s)eq(s +++                   −+      −+      −= (149) The equivalent stress σeq needs to comply with the material yield stress Re: 64 Chapitre 2 : Modélisation du chauffage par induction e eq Rs ≤ . (150) In case the equivalent stress is equal to the material yield stress Re, plastic deformation will occur. The relation between the plastic strain rate and the stress field is governed by the plastic flow law: 0), ≥= ?(s s f ?e pl ∂ ∂& , (151) where f denotes the yield criterion and λ the plastic multiplier. The model can include strain- hardening, by using a power law: ne s Kss eq == 0 , (152) where n denotes the strain-hardening exponent and ε the equivalent deformation defined by: 2 1 2 3 2 )e(e i,j ij∑= . (153) In order to have a general description which can fit material behavior at all temperature ranges, the model has been generalized into an elastic-viscoplastic model, by including material strain rate sensitivity through a Norton-Hoff power law. We therefore introduce this new term in equation (152) which then becomes: mene s Ks . = 0 , (154) where ε& denotes the strain rate and m the strain rate sensitivity. 2.3.3. The physical parameters The material parameters involved in the mechanical model include Young’s modulus and Poisson coefficient, the yield stress and the strain-hardening coefficient, as well as the material consistency and strain rate sensitivity. All these parameters depend on temperature and need to be well determined, especially when dealing with heat treatment applications where the highly localized heat source can induce stresses and distortion that can interact with metallurgy. The evolution of these parameters can here again be provided in two ways: either through the use of constant values or tabulated values of the parameters at some given Chapitre 2 : Modélisation du chauffage par induction 65 temperatures, or by using analytical expressions defining the evolution of these parameters with respect to temperature. 3. The numerical methods Up to this day, most numerical methods for induction heating process modeling have mainly focused on the electromagnetic and thermal computations. Most of these methods are based on one of the three following approaches: - finite elements method, - boundary element method, - mixed finite elements/ boundary elements method. We wish here to carry out in a coupled way electromagnetic, thermal and mechanical computations. More specifically, we need to carry out the following computations: - in the part, electromagnetic, thermal and mechanical, - in the inductor and in the air, electromagnetic and thermal. We have chosen to carry out a complete finite element approach – using six-nodes quadratic triangles - for these three kind of computations. This means that the coupling between the part and the inductor for the electromagnetic computations is carried out through the computation of the electromagnetic field in the air. We detail in the following subsections the specificities and implementation details in each case of our approach. For each of the three problems, we shall build the weak formulation – which we need to have in order to use finite elements, and then space and time discretization. 3.1. The electromagnetic problem 3.1.1. The weak formulation We consider the domain Ω as being: 68 Chapitre 2 : Modélisation du chauffage par induction A basis of the discretised space is provided by the shape functions Ni associated to each node i of the mesh:     = = otherwise.0 ,1 ji j)(nodeN i (162) This means that equation (161) needs only to be verified for all shape functions Ni , and that the unknown Eh can be expressed as: )(.)(.)()( . 1 . 1 xNExNjnodeExE j nodesnb j h jj nodesnb j hh ∑∑ == == . (163) Equation (161) thus becomes: nb.nodes1,i 11 =∀=+ ∂ ∂ ∑∑ == )l(N),N(t).NEb(),N(t).N t E a( iij h j nb.nodes j ij nb.nodes j h j , (164) which can be rewritten: nb.nodes1,i 11 =∀=+ ∂ ∂ ∑∑ == )l(N),N(t).b(NE),N(t).a(N t E iij h j nb.nodes j ij nb.nodes j h j , (165) which is equivalent to the following linear system: [ ] [ ]{ } { }ememem BE(t)K(t) t E C =+       ∂ ∂ , (166) With [ ] [ ] { } ∑ ∫ ∑ ∫∫∫ ∑ ∫ = = =       ∂ ∂ −==       ∂ ∂++∇∇== == nb.elts 1elt elt i s i i em nb.elts 1elt elt ij elt ij2 elt ijij ij em nb.elts 1elt elt ijij ij em dv.N t J )l(NB )dvN(N rµr 1dvNN µr 1dvN.N µ 1)N,b(NK dvNs N)N,a(NC (167) We wish to draw attention here to the fact that using finite elements for electromagnetic computations in the air requires to mesh a large enough domain to avoid any artificial Chapitre 2 : Modélisation du chauffage par induction 69 reflection problems. It will be shown in the numerical result section that a good quality computation must display magnetic field lines that do not intersect the boundary of the global domain. Another point which helps avoid artificial reflection problems is to have well-suited boundary conditions for the air domain. 3.1.3. Time discretization Since we have chosen to solve the time-dependent model, we now need to integrate numerically in time the space-discretized electromagnetic equation (166). For accuracy purposes, we have selected and implemented a second-order two time step finite difference scheme detailed here. In a first stage, the system is solved at time t* such that t < t* < t+δ t2 : 0 with 32123211 =+++++−= aaa)dt(tata)dt(tat* , (168) where δt1 denotes the previous time step and δt2 the current one. The electric field E* at time t* and its time derivative write: 21 321 dtttdtt* EaEaEaE +− ++= , (169) 12 12 1 dt EE )(? dt EE ? t E tdtttdtt * − −+ − = ∂ ∂ −+ . (170) System (166) is written at time t*. E* and its derivative are replaced by expressions (169) and (170). The system is solved for the unknown variable E* : [ ] [ ] { } [ ] [ ] 123 1 2 23 2 12 1 01 23 1 1 1 dt ? dta ?a c dta ?a dt ? dt ? c ECcECcBEKC dta ? dtt*emt*em*em**em*em − −= + − += ++=        + − (171) 70 Chapitre 2 : Modélisation du chauffage par induction The two-time step scheme we have used needs to solve a non- linear equation, as the matrix [Cem] is dependant on the magnetic field. In order to avoid an iterative procedure of resolution of the system, the matrix are linearized [58] in order to depend only on their values at time t and t-δt1 : [ ] [ ] [ ]td tt* C))dt dt (a(aC) dt dt a(aC 1 2 32 1 1 2 31 1 +++−= − . (172) The second and last stage consists in the computation of {E} at time t+δt2: { } { } { } { } )EaEaE( a E tdtt*dtt 2 1 1 3 1 2 −−= −+ . (173) We have chosen here to have δt1 = δt2 = δtem for the electromagnetic time step. The choice of a good value for δtem is not obvious; however, a value of around (T/32) or (T/64) - where T stands for the period of the power current supply - seems to be a good compromise between computation time costs and results precision. We now have a complete numerical scheme which enables us to compute the response in terms of electric field for a given current power supply. As the matrix of the linear system (171) is symmetric and positive definite, we can solve the linear system by using an iterative pre-conditioned conjugate gradient solver; we will use here a diagonal preconditioner. 3.1.4. The moving inductor case In heat treatment applications, we are often faced with the case of inductors moving continuously along the z axis. In the matrix system (171), the [Kem] matrix is modified with the supplementary term Ev rrr ×∇× σ due to the displacement of the inductor as follows: [ ] ∑ ∫ ∫ ∫∫ =               ∂ ∂ − ∂ ∂ ++∇∇ == nb.elts elt elt elt i j zij elt ij elt ij ijij em dvN z N v)dvN(N rµr dvNN µr dvN.N µ ),Nb(NK 1 2 1 11 (174) Chapitre 2 : Modélisation du chauffage par induction 73 3.2.1. The integral formulation The basic approach will be the same as for the electromagnetic problem. We consider the domain Ω ΤΗ as being: inductorpartTH Ω∪Ω=Ω . (177) We define the following boundaries: - ΓΤΗ as being the outer boundary of the ΩΤΗ domain, - Γ0 ΤΗ part of the outer boundary where the temperature is prescribed – Dirichlet boundary condition, - Γ1 ΤΗ part of the outer boundary where the heat flux is prescribed – Neumann boundary condition, - Γ2 ΤΗ part of the outer boundary where there is a convection-radiation boundary condition. We need to establish a weak formulation of equation (140) along with the boundary conditions specified in equations (142)-(144). The functional space V is defined here as:    Γ== ∂ ∂ Ω∈    = TH 0 1 on 0,0 ),(HV ψ θ ψ ψ . (178) After multiplying equation (140) by a test function ψ belonging to the functional space V, integrating on the whole domain, and using the Green theorem, we get: V,? ds, .?Thds .?Fdv ?Qds T.?hdv ?T.kdv ? t T ?C ext G2THG1TH prescribed O em G2THOO ∈∀++=+∇∇+ ∂ ∂ ∫∫∫∫∫∫ & rr (179) which writes: 74 Chapitre 2 : Modélisation du chauffage par induction ds. .?Thds .?Fdv ?Ql(? ( ds, T.?hdv ?T.k? )b(T, dv, ? t T ?C? ), t T a( ,V?),l(? )b(T,? ), t T a( ext G2THG1TH prescribed O em G2THO O ∫∫∫ ∫∫ ∫ ++= +∇∇= ∂ ∂ = ∂ ∂ ∈∀=+ ∂ ∂ & ψ (180) Here again, if we assume ρc and k to be independent on temperature, this problem belongs to the class of parabolic problems, since: - the bilinear form a(.,.) is continuous and positive definite, - the linear form l(.) is continuous, - the bilinear form b(.,.) is continuous and V-elliptic. This problem is also a Cauchy-Dirichle t-Neumann problem and the existence and uniqueness of a solution T(t,x) belonging to L2([0,t];V) can be proven, [43]. 3.2.2. Finite element space discretization We use the same mesh for the part and inductor as the one used for electromagnetic computations. The functional space V is approximated by the discretized space Vh , the test functions ψ by ψh and the unknown T by Th. The discretized version of equation (180) is then: hhhhh nb.nodes j h nb.nodes j h V?)l(?)(t),?Tb()(t),? t T a( ∈∀=+ ∂ ∂ ∑∑ == 11 . (181) Equation (181) needing only to be verified for all shape functions Ni , and the unknown Th being expressed as: Chapitre 2 : Modélisation du chauffage par induction 75 )(.)(.)()( . 1 . 1 xNTxNjnodeTxT j nodesnb j h jj nodesnb j hh ∑∑ == == . (182) Equation (181) thus becomes: nb.nodes1,i 11 =∀=+ ∂ ∂ ∑∑ == )l(N),N(t).NTb(),N(t).N t T a( iij h j nb.nodes j ij nb.nodes j h j , (183) which can be rewritten: nb.nodes1,i 11 =∀=+ ∂ ∂ ∑∑ == )l(N),N(t).b(NT),N(t).a(N t T iij h j nb.nodes j ij nb.nodes j h j . (184) We then get the following equations discretized in space: [ ] [ ]{ } { }ththth BT(t)K(t) t T C =+       ∂ ∂ , (185) With [ ] [ ] { } .NhTNF.NQ)l(NB ,NhNN.Nk)N,b(NK dv, N?CN)N,a(NC nb.elts 1elt Gelt jext Gelt jprescribed elt iemi i th nb.elts 1elt Gelt ij elt ijij ij th nb.elts 1elt elt ijij ij th 21 2 ∑ ∫∫∫ ∑ ∫∫ ∑ ∫ = ∩∂∩∂ = ∩∂ =         ++==         +∇∇== == dsdsdv dsdv & (186) 3.2.3. Time discretisation We now need to integrate numerically in time the space-discretized thermal equation (185). We use the same second-order two-time step finite difference scheme as for the time integration of the electromagnetic problem. This leads us again to a two-stage solving procedure: Stage 1: We define time t* as: 78 Chapitre 2 : Modélisation du chauffage par induction ∫ ∫ = = ∈∀= O O *f.*)l( *)s:e(*)a(u, ?*)l(*)a(u, dv vv dv vv with Vvv (196) The existence and uniqueness of a solution u(x) is proven classically when the material behavior is linear elastic; when the material behavior becomes non- linear, this raises more questions which we will not detail here. 3.3.2. Finite element space discretisation We use the same mesh for the part as the one used for electromagnetic and thermal computations. The functional space V is approximated by the discretized space Vh , the test functions v* by v*h and the unknown displacement field v by vh. The discretized version of equation (196) is then: hhhhh nodes.nb 1j V*v)*v(l)*v),t(v(a ∈∀=∑ = . (197) Equation (197) needing only to be verified for all shape functions Ni , and the unknown vh being expressed as: (x).Nv(x).Nj) (nodev(x)v j nb.nodes 1j h jj nb.nodes 1j hh ∑∑ == == , (198) equation (197) thus becomes: nodes.nb,1i)N(l)N,N).t(v(a iij h j nodes.nb 1j =∀=∑ = , (199) which can be rewritten: nb.nodes1,i)l(N)N,(t).a(Nv iij h j nb.nodes 1j =∀=∑ = . (200) We then get the following equations discretized in space: Chapitre 2 : Modélisation du chauffage par induction 79 [ ]{ } { }meme B)t(vK = , (201) with [ ] { } )N(lB )N,N(aK ii me ijij me = = (202) In case of non- linear mechanical behavior, the matrix KME in system (201) depends on the solution. A Newton-Raphson algorithm is used in that case to converge towards the solution. 3.3.3. Time integration At each time step, system (201) provides us with the velocity field. The geometry of the part is then updated using an exp licit time integration procedure: { } { } { } t.VXX tdttdtt δ+= ++ . (203) 3.4. The coupling procedure As we have said previously, numerical modeling of the induction heating procedure requires the development of a coupling procedure between: - electromagnetic computations - thermal computations - mechanical computations We detail here the strategy which we have selected in order to carry out this coupling. The computation takes place in an incremental way. Two time steps have been defined: - one time step δtem for electromagnetic computations - one time step δtth for thermal and mechanical computations At each time step, the electromagnetic field distribution and thus the induced currents depend explicitly on the thermal field of the body, ruled by the heat transfer equation. Inversely, the 80 Chapitre 2 : Modélisation du chauffage par induction eddy currents calculated from the electromagnetic solution are used as the thermal sources for the thermal finite element analysis. The handling of the coupling is carried out in two different stages: ♦ The criterion from electromagnetic to thermal calculations relying on the stabilization of the mean heat source term over the electromagnetic periods. Once the electromagnetic field has been calculated, the rate of heat generation emQ& for the heat equation needs to be evaluated at every integration points. As the electromagnetic time step is far smaller than the thermal one, we do not consider the instantaneous Joule power calculated at a given time at every integration points. We rather consider a mean Joule power averaged over one or several periods of the electromagnetic field: dttEt T nTQ nT Tn em ∫ − = )1( 2)(int,)(int, 1 int),( θσ , (204) where int is the considered integration point, T is the period of the power supply currents, n is number of periods considered and )(int, tEθ is the value at time t of the electric field interpolated at the integration point int. At the end of each electromagnetic period, the newly calculated mean power is compared to the one calculated at the previous period until it stabilizes. The following convergence test is conducted at every integration points ε< −+ )( )())1(( nTQ nTQTnQ em emem . (205) If (205)is verified then thermal computations are started with the stabilized thermal source power calculated at (n+1)T, where ε is the magnetic convergence parameter. ♦ The criterion from thermal calculations back to electromagnetic calculations is based on the variations of the magnetic and electric parameters with temperature. Thermal computations can use the same source term derived from the electromagnetic computations as long as the physical magnetic parameters such as the magnetic Chapitre 2 : Modélisation du chauffage par induction 83 Billet height = 120mm Billet diameter = 44mm Billet centre hole diameter = 3mm Thickness of Kaowool insulation = 3mm Height of ceramic former = 120mm Internal diameter of ceramic former = 50mm External diameter of ceramic former = 60mm Ceramic material = Fused Alumina Number of layers for copper windings = 3 Number of turns on inner layer = 55 Number of turns on middle layer = 54 Number of turns on outer layer = 55 Diameter of copper wire = 2mm Length of coil on former = 117mm Two thermocouples have been placed: the first one on the surface on the middle of the part and the second one in a hole bored in the center of the part so to measure temperature inside the part. 84 Chapitre 2 : Modélisation du chauffage par induction Figure 13: Static long inductor case: geometry and mesh (12000 nodes) 4.1.1.1. Magnetic and thermal data of the EN3 steel The physical data are given by a set of physical values for different temperatures and computed at a given temperature by linear interpolation. 0.0E+00 2.0E+06 4.0E+06 6.0E+06 8.0E+06 1.0E+07 1.2E+07 0 200 400 600 800 1000 1200 1400 Temperature (°C) Electrical conductivity Specific heat 0 10 20 30 40 50 60 0 200 400 600 800 1000 1200 1400 Temperature (°C) Thermal conductivity Figure 14: Electrical conductivity, specific heat and thermal conductivity versus Temperature Chapitre 2 : Modélisation du chauffage par induction 85 The magnetic permeability of the EN3 steel is dependent on the magnitude of the magnetic field and follows the law: 20 )( ),( H H H HT r r r r + − + += β α β αµµ , where H r is the magnitude of the magnetic field, α and β are Frohlich coefficients (α = 2.31695 and β = 5769.55), and 0µ is the magnetic permeability of the air. The data we get on the EN3 steel for the magnetic permeability did not take into account the temperature dependency of the magnetic permeability. 4.1.1.2. Process parameters Frequency (Hz) 500 Current density (Amps/m2) 8.108 Electromagnetic time step (s) T/64 Thermal time step (s) 0.5 4.1.2. Results and validation In order to enable a visualization of the joule power, source term for the thermal computations, an effective electrical field was defined at a given node as the integration over an electromagnetic period of the square of the real instantaneous electrical field: ∫ + = T1)(n Tn 2 eff dt*(t)E T 1 E 88 Chapitre 2 : Modélisation du chauffage par induction 0.0 50.0 100.0 150.0 200.0 250.0 0.00 5.00 10.00 15.00 20.00 25.00 30.00 35.00 Time (s) T em pe ra tu re ( °C ) Experimental surface Temperatures Experimental heart Temperatures Numerical surface temperature numerical heart temperature Figure 17: Long inductor case: comparison between experimental and computed temperature evolutions Figure 18 shows the temperature results at two given time steps. Heat propagates here rather uniformly from the surface towards the center. This is typical of a global homogeneous heating application. time=16 s time=40 s Chapitre 2 : Modélisation du chauffage par induction 89 Figure 18: Long inductor case: Temperature field at two given time steps Figure 19 shows the results obtained in terms of velocity field for the coupled mechanical computation. One can observe the dilatation effects taking place in the workpiece while it is being heated. Figure 19: Long inductor case: nodal velocities in the mesh at a given time step 4.2. A static short inductor case 4.2.1. Description of the case The geometry of the case studied here has been presented in a paper by Wang [55], which details the analytical solutions for two different constant magnetic permeability of the part: 0µµ = and 0*90 µµ = . Only the analytical expression for a relative magnetic permeability of one is valid. For a relative magnetic permeability greater than 1, a supplementary due to the 90 Chapitre 2 : Modélisation du chauffage par induction radial derivative of the magnetic permeability term and which is missing in the previously cited article, needs to be added, especially on the interface between the part and the air, and considerably modifies the electromagnetic field. The electromagnetic and temperature fields are shown as well as the mechanical deformations. Figure 20 shows the mesh which has been refined along the workpiece surface in order to have several elements within the electromagnetic skin-depth. Figure 20: Static short inductor case: geometry and mesh The physical and numerical parameters are summarized in the following tables: Process parameters Current density (A.m-2) 7 109 Electromagnetic time step (s) T/32 (T= period) Thermal time step (s) 1 Magnetic and thermal data The physical properties are given for steel and copper (respectively part and inductor). In this case, the properties are constant in time. STEEL (PART) COPPER (COIL) AIR σ (Ω-1 .m-1) 3 106 60 106 0 Chapitre 2 : Modélisation du chauffage par induction 93 Figure 24 provides us with a display of magnetic field lines for a magnetic and non-magnetic material. Figure 24: Static short inductor case at 60Hz: magnetic lines for rµ =1 and rµ =90 The temperature computations are carried over the whole domain. Heat propagates from the surface centre towards the rest of the workpiece and is also conducted in the air. Figure 25 and Figure 26 show that for a given frequency, temperature distributions can be totally different depending on the nature of the material. Figure 25: Temperature at two time steps at a frequency of 60 Hz for a non magnetic part (µrelative =1); a) t=5s b)t=25s 94 Chapitre 2 : Modélisation du chauffage par induction Figure 26: Temperature at two time steps for a magnetic part (µrelative =90); a) t=1s b)t=5s A thermo-mechanical coupled analysis has been carried out. As it can be expected, the plot of the nodal velocity field shows dilatation of the part in the heated regions and contraction on both edges. a) Figure 27: Nodal velocities in the workpiece due to dilatation effects for a frequency of 60Hz; a) non magnetic material (µrelative =1), b) non magnetic material (µrelative =90) Chapitre 2 : Modélisation du chauffage par induction 95 4.3. A moving short inductor case This case deals with a moving short inductor - a case representative of what takes place in an induction heat treatment process. As explained previously, the displacement of the induc tor is modeled through a continuous variation of physical properties of the mesh. Figure 28 shows the geometry and the initial location of the inductor in the mesh. Figure 28: Moving inductor case: geometry and initial mesh Figure 29 shows the temperature distribution at a given time step; one can notice the influence of the inductor displacement on the temperature field in the workpiece. Figure 29: Moving inductor case: temperature isovalues at a given displacement step of the inductor 98 Chapitre 2 : Modélisation du chauffage par induction Les calculs ont été effectués sur la géométrie du cas présenté dans l’article de Wang [55] : Figure 30: Géométrie et maillage Le maillage est composé d’éléments quadratiques à 6 nœuds et est composé de 3556 nœuds. Les paramètres procédé sont les suivants : Fréquence 60 Hz Densité de courant source 7 109 A.m-2 Pas de temps électromagnétique T/32 (T : période) Pas de temps thermique 0.5 s Figure 31: les paramètres procédés Les données physiques, magnétiques et thermiques, sont données pour un acier type et pour le cuivre qui composent respectivement la pièce et l’inducteur. Dans ce cas idéalisé, les propriétés physiques ne dépendent pas de la température. Acier Cuivre Air σ (Ω-1 .m-1) 3 106 60 106 0 µr 1 1 1 ρc (J.Kg-1.K-1) 4.87 106 3.43 106 1.3 103 k (W.m-1.K-1) 35 389 2.63 10-2 Figure 32: Données physiques des différents matériaux : σ est la conductivité électrique, µr la perméabilité magnétique, ρc la chaleur spécifique et k la conductivité thermique. Chapitre 2 : Modélisation du chauffage par induction 99 Un courant sinusoïdal est injecté dans l’inducteur. Etant donné que la pièce est non ferromagnétique, le champ électromagnétique obtenu par un calcul linéaire, devrait avoir la même configuration sinusoïdale avec toutefois un déphasage possible. Ce déphasage est de l’ordre de 4/π par rapport à la source de courant. Nous avons testé les différents types de schéma d’intégration en temps. Le champ électrique montre des oscillations parasites de grande amplitude pour des schémas plutôt explicites (Crank-Nicolson ou Lees) et ont tendance à diminuer, voir à complètement disparaître avec un schéma de Dupont ou complètement implicite, Figure 33. Ainsi un calcul électromagnétique, demandé avec une précision de 2% entre deux périodes électromagnétiques consécutives, va converger en 16 périodes avec un schéma de Crank Nicolson et va nécessiter plus de 60 périodes avec un schéma de Lees! -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 0 0.01 0.02 0.03 0.04 0.05 0.06 time (s) E le ct ri ca l f ie ld ( V /m ) CN LEES DUPONT implicite Figure 33 : Evolution temporelle du champ électrique en un point de la surface de la pièce. Ces oscillations parasites sont dues en partie à une estimation exagérée de l’amplitude du champ électrique aux pas de temps de la première période. Le schéma de Dupont et le schéma implicite se montrent stables et le calcul va converger en 3 périodes. 100 Chapitre 2 : Modélisation du chauffage par induction 2.2 Initialisation des schémas d’intégration en temps Un autre point important concerne l’initialisation du schéma d’intégration en temps. En effet, le schéma d’intégration à deux pas de temps nécessite de connaître les valeurs du champ électrique imposé 0=tEθ au temps initial t=0 mais aussi au temps ttt ∆−= . Une première solution consiste à initialiser le champ électrique ttE ∆−θ au temps ttt ∆−= aux mêmes valeurs que le champ tEθ . Une autre solution est d’utiliser un schéma de Crank-Nicolson pour le premier pas de temps. Les Figure 34 et Figure 35 montrent l’influence du type d’initialisation sur les schémas d’intégration en temps de Lees et de Dupont. Pour un schéma de Lees, une initialisation par un calcul de Crank-Nicolson donne un meilleur résultat en réduisant l’amplitude des oscillations. En revanche pour le schéma de Dupont, une initialisation directe des champs donne un meilleur résultat en ne surestimant pas l’amplitude du champ électrique au premier pas de temps. -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 0 0,005 0,01 0,015 0,02 0,025 0,03 0,035 0,04 0,045 0,05 Temps (s) C h am p é le ct ri q u e (V /m ) Initialisation par Crank-Nicholson Initialisations par valeurs imposées Figure 34: Evolutions temporelles du champ électrique calculées pour deux initialisations différentes du schéma d’intégration de Lees. Chapitre 2 : Modélisation du chauffage par induction 103 3.1 Cas linéaire Différents pas de temps, définis comme des fractions de la période ont été testés. Un pas de temps de l’ordre de T/20 est le minimum suffisant pour commencer à avoir des résultats fiables. Un pas de temps trop grand va avoir tendance à exagérer l’amplitude. Le signal est de plus en plus déphasé au cours du temps par rapport aux solutions réelles, ce qui montre bien que la fréquence elle-même est mal reproduite. -25 -20 -15 -10 -5 0 5 10 15 20 25 0 0.01 0.02 0.03 0.04 0.05 0.06 Temps (s) C ha m p él ec tr iq ue ( V /m ) dt=T/10 dt=T/20 dt=T/50 dt=T/100 Figure 37: Evolution temporelle du champ électrique en un point de la surface de la pièce pour différents pas de temps électromagnétique, T étant la période du signal électrique. 3.2 Cas non linéaire Le choix du pas de temps est primordial. Un pas de temps trop sous-estimé va complètement modifier la solution du calcul, voire même entraîner un problème de convergence. En effet, l’induction magnétique B r est calculée de manière explicite à partir du champ électrique E r . Une bonne estimation du champ B r est nécessaire car sa norme va permettre le calcul de la norme du champ magnétique H r dont va ensuite dépendre le calcul de la perméabilité magnétique ),( TH r µ pour le calcul du nouveau champ électrique. Un pas de temps assez 104 Chapitre 2 : Modélisation du chauffage par induction petit est indispensable. En revanche un pas de temps excessif va considérablement alourdir les temps de calculs en retardant la convergence des calculs et peut conduire à une perte de précision. Nous avons testé différents pas de temps et leur influence sur le calcul électromagnétique et thermique. Il est difficile d’établir une règle précise pour le choix du pas de temps, ce dernier dépendra de la fréquence, de l’intensité du courant appliqué à l’inducteur mais également de la valeur de perméabilité magnétique du matériau de la pièce. Une valeur trop grande, en général supérieure à dt=T/64 va générer une instabilité numérique après quelques pas de temps électromagnétiques. Ceci se traduit par une augmentation très importante du nombre d’itérations du solveur, voire même une divergence. Le cas test utilisé est celui décrit dans la publication 1, § 4.1.1, se référer à la Figure 13. -8 -6 -4 -2 0 2 4 6 8 0.04 0.041 0.042 0.043 0.044 0.045 0.046 Temps (seconde) C h am p é le ct ri q u e (V /m ) surface inducteur air surface de la pièce champ interne à la pièce Figure 38 : Evolution temporelle du champ électrique en différents points du domaine pour une fréquence de 500Hz, un pas de temps dt=T/128, avec la période T=0.002s. Le champ électromagnétique se déforme et ne s’apparente plus à une sinusoïde. L’amplitude maximum du champ électrique se situe au niveau de la surface de la pièce et non plus dans Chapitre 2 : Modélisation du chauffage par induction 105 l’inducteur comme pour un matériau linéaire. Les déformations restent néanmoins de faible amplitude dans l’air et dans l’inducteur. ch am p é le ct ri q u e (V /m ) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0.04 0.041 0.042 0.043 0.044 0.045 0.046 temps (seconde) dt=T/128 dt=T/256 Figure 39 : Modification du profil électrique avec le pas de temps à l’intérieur de la pièce 4 Stockage des données et méthode de résolution du système matriciel La résolution numérique des systèmes matriciels résultants des calculs électromagnétique et thermique est effectuée par un solveur itératif. Le choix d’une méthode itérative plutôt que directe est justifié d’abord par un gain en temps de calcul mais également par un gain important en termes de stockage de données, les solveurs itératifs ne faisant appel qu’à des produits matrice-vecteur. Enfin, ce choix nous a permis de nous orienter vers une stratégie de type partitionnement de domaine pour le développement de la version parallèle. La matrice du système linéaire résultant du modèle électromagnétique ou thermique étant symétrique et définie positive, nous nous sommes naturellement orientés vers un solveur de type gradient conjugué préconditionné. Concernant le mode de stockage des données, étant donné que le solveur itératif utilise des produits matrice-vecteur, nous avons pu utiliser un stockage compact de type Morse pour la matrice du système. 108 Chapitre 2 : Modélisation du chauffage par induction où AxxAAxbxr −=−=)( est le résidu. La fonctionnelle E est minimisée par une méthode de descente qui revient à déterminer, à la kième itération, la direction de descente 0≠kp et un scalaire kα telle que : )()( 1 kk xExE <+ , avec kkkk pxx *1 α+=+ , 1* −+= kkkk prp β . Nous présentons ici l’algorithme de gradient conjugué préconditionné. Il converge théoriquement en au plus n itérations, où n est le nombre d’itérations. Chapitre 2 : Modélisation du chauffage par induction 109 Algorithme du gradient conjugué préconditionné Initialisations X0 donné r0=b-A x0 C p0= r0 z0= p0 Minimisation de J(xk+1) kkk1k kkk1k kk kk k pAarr paxx )p,p(A )z,(r a −= += = + + k=k+1 Calcul de la nouvelle direction k1k1k1k kk 1k1k 1k k1k pßzp )z,r( )z,(r ß rz C +++ ++ + ++ += = = 1 Critère d’arrêt ε≤ +1k r ?? NON OUI FIN 110 Chapitre 2 : Modélisation du chauffage par induction Les performances de l’algorithme de gradient conjugué dépendent du conditionnement de la matrice du système. Plus )(Acond est proche de 1, plus vite convergera l’algorithme. Pour un système mal conditionné ; il est nécessaire d’utiliser un préconditionneur C, matrice non- singulière, symétrique et définie positive. On considère alors le système suivant : bCAxC 11 −− = . 4.3 Le préconditionnement Un bon préconditionneur possède les propriétés suivantes : − )()( 1 ACondACCond <− : en théorie, le choix optimal pour la matrice C serait C=A de telle manière que 1)( 1 =− ACCond . Il faut déterminer C-1 aussi proche de A-1 que possible pour que )( 1ACCond − soit aussi proche que possible de 1. − C doit être aussi creuse que A pour ne pas trop alourdir les tailles de stockage. − C doit être facilement calculable et inversible pour ne pas trop alourdir les temps de calculs. En pratique, il s’agit de trouver le bon compromis entre 1−C le plus proche possible de 1−A et des calculs de 1−C pas trop coûteux. Il existe plusieurs types de préconditionneurs, que l’on va décrire ci-dessous. 4.3.1 Le préconditionneur diagonal Ce preconditionneur consiste simplement à reprendre les termes diagonaux de la matrice A : )(AdiagC = . Il est particulièrement intéressant en termes de stockage, en temps et coût de calcul. Chapitre 2 : Modélisation du chauffage par induction 113 4.4 Discussion Un solveur itératif a été implémenté et testé avec succès sur les systèmes matriciels électromagnétiques et thermiques. Il a permis de gagner en temps de calcul CPU mais également en place mémoire, étant donné que lui est associé un stockage compact de type Morse. Par ailleurs, au vu des performances des différents préconditionneurs, nous avons décidé d’utiliser un préconditionneur diagonal. 5 Influence du terme de dérivée de la perméabilité magnétique Le terme de dérivée radiale (59) de la perméabilité magnétique est rarement pris en compte dans les modèles électromagnétiques vus dans la littérature. Néanmoins il peut changer considérablement la forme de la solution. Reprenons le cas test présenté Figure 30 avec simplement la perméabilité magnétique relative de la pièce qui est maintenant égale à 50. On lance deux calculs s’appuyant sur deux modèles différents : le premier prend en compte le terme supplémentaire et le deuxième le néglige. On trace le profil radial du champ électrique dans le domaine -20 -10 0 10 20 30 40 50 60 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 coordonnée radiale (mètre) C ha m p él ec tr iq ue (V /m ) avec le terme volumique sans le terme volumique Figure 40 : Profils radiaux de champ électrique calculés à partir de modèles considérant le terme volumique supplémentaire ou non Inducteur Pièce 114 Chapitre 2 : Modélisation du chauffage par induction La Figure 41 montre la différence de montée en température que cela implique en surface de la pièce. C’est en effet en surface de la pièce que la différence de température, due à la prise en compte de ce terme, est la plus exacerbée (par rapport aux zones internes). 0 500 1000 1500 2000 2500 3000 0 1 2 3 4 5 6 7 8 Temps (seconde) T em p er at u re ( ° C ) sans le terme volumique avec le terme volumique Figure 41 : Evolutions de la température en surface de la pièce calculées à partir de modèles considérant le terme volumique supplémentaire ou non Enfin la Figure 42 nous montre l’évolution du profil radial du champ électrique calculé avec le modèle complet pour différentes perméabilités. On remarque que dès que notre matériau devient magnétique, l’amplitude maximum du champ électrique dans tout le domaine se retrouve sur la surface de la pièce. Chapitre 2 : Modélisation du chauffage par induction 115 -20 -10 0 10 20 30 40 50 60 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 coordonée radiale (mètre) ch am p é le ct ri q u e (V /m ) mu=1 mu=50 mu=90 Figure 42 : Variation du profil du champ électrique calculé avec le modèle complet pour différentes perméabilités magnétiques relatives 6 Eléments de validation supplémentaires La publication précédente donne des exemples de validation de notre modèle par rapport à des mesures expérimentales. Des études de validation supplémentaires ont été menées d’une part par comparaison avec le code de nos partenaires slovènes et d’autre part par comparaison avec une solution analytique, tirée de l’article de Wang [55]. 6.1 Comparaison entre les codes du CEMEF et de LNMS Nous avons comparé des résultats de simulation entre les codes de CEMEF et de LNMS (Laboratory for Numerical Modeling and Simulation) de l’université de Ljubljana. Le code développé par LNMS résout le problème électromagnétique en s’appuyant sur une formulation harmonique. La méthode de résolution est une méthode mixte éléments finis – éléments frontières. Les codes sont comparés sur un cas linéaire avec une pièce non magnétique. En effet dans ce cas, la formulation harmonique est tout à fait valable et les deux codes sont censés donner exactement les mêmes résultats. Inducteur Pièce 118 Chapitre 2 : Modélisation du chauffage par induction 0 100 200 300 400 500 600 700 800 900 1000 0 5 10 15 20 25 Temps (s) Te m pe ra tu re ( °C ) point en surface point au milieu point sur axe de symetrie a) b) Figure 44 : Evolution temporelle de la température en surface de la pièce, sur l’axe de symétrie et entre les deux. a) solution analytique tirée de l’article de Wang [55] ; b) solution numérique 7 Applications à un cas industriel de traitement thermique d’un levier de boîte de vitesse Le cas industriel que nous présentons concerne le traitement thermique d’un levier de boîte de vitesse. Ce test nous a été proposé par le groupe italien TQT-SIAP, partenaire du projet européen Heatmaster. Le chauffage par induction est ici utilisé pour durcir par traitement thermique certaines zones de la pièce (Figure 45) afin d’augmenter la dureté en surface, la résistance à la fatigue et à l’usure. Les zones traitées doivent subir des évolutions temporelles précises de températures, ces dernières étant dictées par les courbes de transformation isotherme ou courbes TTT. Les fréquences couramment utilisées par TQT sont de l’ordre de 10kHz. Les zones traitées, initialement de structures austénitiques, sont chauffées rapidement à une température comprise entre 900°C – 950°C puis refroidies brusquement afin d’obtenir une struc ture martensitique. Chapitre 2 : Modélisation du chauffage par induction 119 Figure 45 : cette photo montre la pièce avant traitement thermique (pièce de droite) et après (pièce de gauche) Les caractéristiques géométriques, ainsi que les zones à traiter et les zones de déplacement des inducteurs sont représentées sur la Figure 46. Les deux zones en bas de pièce sont traitées l’une après l’autre successivement, puis l’inducteur est retiré, la pièce est retournée automatiquement et le traitement de la dernière zone s’effectue. 120 Chapitre 2 : Modélisation du chauffage par induction Figure 46 : Géométrie du levier de boîte de vitesse. Les zones en rouge représentent les zones subissant un traitement de surface. Les zones de déplacement de l’inducteur sont également représentées. Chapitre 2 : Modélisation du chauffage par induction 123 Figure 47 : Iso valeurs du champ électrique effectif pour deux positions de l’inducteur. La Figure 48 montre l’évolution des iso de température au cours du déplacement de l’inducteur. 124 Chapitre 2 : Modélisation du chauffage par induction Figure 48 : Iso-valeurs des températures au cours du temps. La simulation permet alors de modifier les différents paramètres procédés, afin de se rapprocher des objectifs en termes d’évolutions temporelles de température. L’épaisseur de la couche matensitique peut ainsi être calculée à partir d’un diagramme TTT. Chapitre 2 : Modélisation du chauffage par induction 125 8 Conclusion Dans ce chapitre, nous avons présenté les modèles électromagnétiques, thermo-mécaniques, ainsi que leurs discrétisations temporelles et spatiales. La stratégie de couplage a été décrite ainsi que le choix du solveur, en l’occurrence un gradient conjugué avec un préconditionneur diagonal. Des études complémentaires sur les pas de temps, les schémas en temps ainsi que des tests de validation ont été nécessaire afin d’obtenir un modèle direct fiable et robuste. Ces qualités sont essentielles pour l’optimisation du procédé par analyse inverse, [22]. Ce modèle s’avère être très souple, pouvant traiter n’importe quel matériau et n’importe quelle intensité d’entrée dans l’inducteur. Néanmoins, il peut se révéler coûteux en temps de calcul, notamment au niveau de la résolution magnétique qui peut nécessiter un nombre important de périodes avant de converger, ou bien dès lors que l’on prend en compte un déplacement de l’inducteur. La prochaine étape nécessite la parallélisation du code afin de réduire les temps de calcul et de permettre l’étude d’applications avec des tailles de maillage importantes. 128 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction Les choix du niveau de parallélisation d’un code séquentiel et d’une méthode de parallélisation vont être dictés par - le type d’architecture de la machine parallèle (parallélisme massif ou à grains moyens ou grossiers) - la méthode numérique employée (éléments finis, méthode des intégrales de frontière) ainsi que la nature du système matriciel issu des discrétisations spatiale et temporelle du modèle physique Le but est d’obtenir les meilleures performances possibles. En effet un même code parallèle peut se montrer efficace sur un calculateur parallèle donné mais être complètement inadapté à un autre configuration de machine. D’autre part, une qualité essentielle d’un code commercial doit être sa portabilité, c’est à dire sa capacité à rester efficace d’une architecture à l’autre ou du moins sur des types ou famille de calculateurs donnés. Nous allons commencer par passer en revue dans un premier paragraphe les notions de base en calcul parallèle, les types d’architecture des calculateurs parallèles ainsi que les caractéristiques primordiales d’un logiciel parallèle performant en termes de gain en temps de calcul. Comme nous le verrons, nous avons décidé de nous orienter vers une stratégie de programmation de type « Message Passing » réellement adaptée aux méthodes éléments finis implicites, qui sont celles qui nous concernent. Nous détaillons, dans le second paragraphe, les différentes méthodes existantes basée sur le paradigme de partitionnement de maillage et de programmation par échange de messages. La méthode utilisée est ensuite exposée dans la publication que nous avons inséré et les résultats en termes de mesure de performances sont détaillés. 1 Généralités sur le calcul parallèle 1.1 Mesure des performances d’un code parallèle Pour estimer les performances et donc la qualité d’un code parallèle, il faut quantifier les gains en temps de calcul du code parallèle par rapport au code séquentiel. Ces mesures de performances peuvent varier suivant le nombre de processeurs et la taille du domaine étudié en termes de nombre de nœuds. En effet un algorithme parallèle donné peut s’avérer moins efficace si le nombre de processeurs utilisés devient trop important, ceci notamment à cause des coûts de communications. Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 129 Deux grandeurs sont communément utilisées pour quantifier les performances d’une implémentation parallèle : − l’accélération ou speed-up C’est le rapport entre le temps d’exécution séqt du code séquentiel et le temps Npt d’exécution du code parallèle sur pN processeurs Np seq t t S = (208) Le but est d’obtenir un code parallèle offrant un facteur d’accélération maximal. En particulier, pour un calcul lancé sur pN processeurs, l’idéal serait d’obtenir un algorithme s’exécutant pN fois plus vite que le meilleur algorithme séquentiel. Quand de telles performances sont atteintes, S= pN ou du moins S de l’ordre de pN , on dit que l’algorithme parallèle atteint un facteur d’accélération linéaire. - l’efficacité C’est le rapport de l’accélération sur le nombre pN de processeurs utilisés 100(%) ×= pN S E (209) L’efficacité d’un code parallèle peut chuter assez rapidement avec le nombre de processeurs, on dit dans ce cas que l’algorithme n’est pas extensible ( ou encore scalable). 1.2 Architecture des calculateurs parallèles Il existe une grande variété de calculateurs parallèles, regroupant dans la même catégorie aussi bien des supercalculateurs vectoriels que des réseaux de stations de travail. Chaque type de calculateur possède une architecture différente qui va être caractérisée par l’organisation de sa mémoire, sa granularité, la topologie des connections entre les processeurs et par le mode d’organisation des tâches des processeurs. 130 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction a/ Granularité La granularité d’une machine est le rapport entre la puissance totale de calcul et le nombre de processeurs. Un calculateur de granularité fine (élevée) va posséder un grand nombre de processeurs à faible puissance de calcul. Ce type de machine va permettre de faire du parallélisme massif. Les calculateurs à grain moyen ou grossier possèdent inversement un petit nombre de processeurs à forte puissance de calculs. b/ L’organisation de la mémoire On dit que la mémoire est partagée si tous les processeurs accèdent à une mémoire commune (Figure 49) et ne disposent pas de mémoire locale. Figure 49 : Architecture à mémoire partagée sur un calculateur à quatre processeurs Inversement, la mémoire est dite distribuée si chaque processeur dispose d’une mémoire locale (Figure 50) Figure 50 : Architecture à mémoire distribuée sur un calculateur à quatre processeurs Des organisations intermédiaires, dites à mémoire hiérarchique (Figure 51), existent. Les processeurs ont accès à une mémoire partagée et possèdent chacun une mémoire locale qui leur est propre. Cette dernière catégorie regroupe la plupart des calculateurs parallèles, la mémoire locale étant au minimum la mémoire cache. MEMOIRE P1 P2 P3 P4 P1 P2 P3 P4 M1 M2 M3 M4 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 133 Il existe principalement deux niveaux de langage de programmation parallèle : - Langage à parallélisme de données « data parallel » : ce langage est adapté à des calculateurs SIMD à grains fins et à mémoire partagée. Le parallélisme est géré automatiquement dès la compilation, la mémoire centrale est partagée virtuellement par le compilateur. Le programmeur doit organiser la séquence d’exécutions de manière à ce que les données appartenant à une même zone mémoire ne soient pas modifiées simultanément. On peut citer le langage HPF (High Performance Fortran) ou OpenMP. Le lecteur peut se référer à [56] pour une présentation complète de ce langage. - « Message passing » : les communications et échanges de données entre processeurs sont gérés par le développeur qui va utiliser des fonctions de base, pour initialiser le réseau, envoyer et réceptionner des messages de manière synchrone ou asynchrone. Ce mode de programmation est très répandu sur les machines MIMD à mémoire distribuée, mais peut aussi bien être utilisé sur les machines à mémoire partagée. Dans ce cas, les messages ne sont pas explicitement envoyés, la librairie gérant alors simplement l’accès à la mémoire. Deux librairies sont majoritairement utilisées : PVM (Parallel Virtual Machine) et MPI (Message Passing Interface). Le programmeur a le choix entre plusieurs stratégies de programmation : - MIMD (Multiple Instruction Multiple Data) : les séquences de code sont développés différemment pour chaque processeur qui effectue des tâches différentes - Maître-Esclave : un processeur maître distribue les tâches aux processeurs esclaves - SPMD (Single Program Multiple Data) : tous les processeurs exécutent le même programme. Une programmation basée sur les librairies standards explicites PVM ou MPI va générer un code facilement portable sur différent types d’architectures. En revanche un code généré spécifiquement par un langage à parallélisme de données va se révéler plus efficace sur une machine à mémoire partagée mais va être difficilement portable. 134 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 1.4 Optimisation du code parallèle Les règles d’or pour pouvoir limiter au maximum les pertes de performances consistent à : a/ paralléliser au maximum le code : En effet, la loi d’Amdhal (211) montre que l’accélération maximum pouvant être atteinte par un code parallèle est limitée par la fraction de temps de calcul effectuée en séquentiel, [56]. Pour démontrer cette loi, on peut décomposer le temps d’exécution d’un code séquentiel en une partie séquentielle et en une partie parallélisable : // seq seq seqseq ttt += , où seqseqt est le temps d’exécution séquentiel de la fraction de code non parallélisable et // seqt est le temps d’exécution séquentiel de la fraction de code parallélisable. De même, on peut décomposer le temps d’exécution de code parallèle sur Np processeurs en une partie séquentielle et une partie parallèle : // Np seq NpNp ttt += , où seqNpt est le temps d’exécution parallèle de la fraction de code non parallélisable et // Npt est le temps d’exécution parallèle de la fraction de code parallélisable. On a égalité entre les temps d’exécution séquentiels seqseqt et seq Npt . D’autre part, pour un code parallèle idéal, on aurait une accélération linéaire : Np t t seqNp // // = . On introduit λ, la proportion de calcul parallélisable : seq seq t t // =λ . (210) L’accélération (208) devient : PN S λ λ +− = )1( 1 . (211) Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 135 Si λ=1, le code est complètement parallélisable et NpS = . Si λ=0, le code n’est pas du tout parallélisable et S=1. Ainsi une parallélisation partielle d’un code séquentiel n’est pas une bonne approche. Cette fraction de temps séquentielle se traduit soit par un calcul redondant sur tous les processeurs (la tâche séquentielle est effectuée par tous les processeurs sans communication), soit le calcul est effectué sur un processeur, les autres étant en attente jusqu’à ce que le résultat de la tâche séquentielle leur soit communiqué. On s’aperçoit que dans l’une ou l’autre méthode, la perte d’efficacité va dépendre directement du nombre de tâches à effectuer en séquentiel. En conséquence, il s’agit de minimiser au maximum la fraction non parallélisable d’un code séquentiel. b/ Minimiser les temps de communications : le temps de transfert des données d’un processeur vers un autre se décompose classiquement de la manière suivante : octettamponinitioncommunicat tNttt ++= , (212) Avec ioncommunicatt : temps de transfert du message de taille N octets, initt : temps d’initialisation de la communication, tampont : temps d’affectation des données dans une mémoire tampon, octett : temps de transmission d’un octet. Ces temps d’initialisation et de débit vont dépendre des caractéristiques intrinsèques du calculateur. Si le temps d’initialisation est important, on aura intérêt à condenser les communications pour envoyer un minimum de messages longs plutôt que beaucoup de messages courts. Le débit des données va dépendre de la topologie du réseau. Un réseau en cross-bar va transférer simultanément les données vers tous les processeurs alors qu’un réseau en anneau va transférer les données de manières séquentielles avec un temps de latence non négligeable. c/ répartir uniformément les charges de travail entre les différents processeurs : Les processeurs moins sollicités, une fois leur tâche immédiate terminée vont devoir attendre passivement le dernier processeur actif. Un déséquilibre des charges se traduit clairement par 138 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction Nous allons spécifiquement nous intéresser à l’application de ces méthodes à un problème de type parabolique. Pour des problèmes de type elliptiques, nous renvoyons le lecteur aux thèses de E. Perchat [40], S Marie [33] et à l’article de Chan and al [11]. Nous allons décrire les algorithmes de décomposition de domaine appliqués à la résolution du système d’équations obtenu par une discrétisation implicite du problème parabolique initial suivant pour ],0[),( Ttr ×Ω∈ r :         ×Ω∂= Ω= ×Ω=− ∂ ∂ ],0[0),( )()0,( ],0[ 0 Tsurtru surruru TsurguL t u r rr (213) où ),( rtgg r = , )(0 ru r sont donnés et L est un opérateur elliptique. Bien que les méthodes présentées s’appliquent à un opérateur elliptique L quelconque, nous allons nous intéresser plus spécifiquement à la résolution du système parabolique suivant (214) avec L de la forme .)( ∇⋅∇= rr aL , le paramètre a représentant soit le coefficient de transfert thermique k dans le cadre de la thermique, soit a=1/µ avec µ la perméabilité magnétique dans le cadre de la résolution électromagnétique. On précise que pour le modèle électromagnétique, des termes supplémentaires, de type b.u viennent s’ajouter. Leur présence ne modifiant pas la méthode décrite, nous ne les mentionnons pas ici dans un souci de clarté.            ×∂= = ×=∇⋅∇− ∂ ∂ T][0,Osur0t),ru( Osur)r(u,0)r(u T][0,Osurgu)(a t u 0 r rr rr (214) On considère une discrétisation éléments finis en espace et en différences finies en temps. Par exemple avec un schéma d’Euler implicite, on obtient le système linéaire : { } { } { } { } { } { }00 111 /)( uu guMdtuu n nnnn = +−=− = +++ , (215) Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 139 où M est la matrice creuse, symétrique et définie positive issue de la discrétisation éléments finis de )( ua ∇⋅∇ rr , et dt est le pas de temps. Des équations similaires sont obtenues avec un schéma de Crank-Nicolson ou à deux pas temps. A chaque itération, le système linéaire suivant doit être résolu : { } { } { } 11)( ++ +=+ nnn gdtuuMdtI , (216) où I est la matrice identité. Pour des raisons de simplicité d’écriture, nous allons considérer le système linéaire : { } { } 11 ++ = nn fuA , (217) avec )( MdtIA += et { } { } { } 11 ++ += nnn gdtuf , A est une matrice creuse, symétrique. Il existe principalement trois types de méthodes de décomposition de domaine que nous allons passer en revue : deux méthodes primales et une méthode duale − Méthode alternative de Schwarz − Méthode du complément de Schur primale − Méthode du complément de Schur duale La description des ces méthodes est illustrée sur une partition du domaine initial en deux sous-domaines pour facilité la compréhension. Une extension de ces méthodes à un nombre plus important de sous-domaines est immédiate. Nous verrons les problèmes qu’elle engendre. 2.1.1 Les méthodes de décomposition de domaine primales a/ Approche avec recouvrement des sous-domaines : méthode alternative de Schwarz Le domaine Ω est décomposé en plusieurs sous-domaines se recouvrant partiellement. Par exemple, soit deux sous-domaines { }21 ,ΩΩ recouvrant Ω, Figure 52. 140 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction Figure 52 : Décomposition du domaine Ω en deux sous-domaines avec recouvrement Avec 2,1, =Ω iiδ représentant la frontière du sous-domaine iΩ et 2,1, =Γ ii la frontière du sous-domaine iΩ intérieure au domaine Ω. La méthode de Schwarz consiste à découpler les problèmes locaux sur chaque sous-domaine et à effectuer un point fixe sur le recouvrement. Les conditions aux limites sur les frontières extérieures iΩδ \ iΓ sont les mêmes que celles du problème global sur ces mêmes frontières. En revanche, sur les frontières internes 2,1, =Γ ii , on impose des conditions de Dirichlet où u prend les valeurs en frontière du calcul précédemment effectué sur le sous-domaine adjacent. La résolution du problème global s’effectue par une méthode itérative sur les sous-domaines de type Gauss- Seidel par blocs (méthode multiplicative de Schwarz), ou de type Jacobi par blocs (méthode additive de Schwarz). L’écriture de la forme discrétisée des algorithmes nécessite l’introduction pour le domaine iΩ de la matrice TiR , constituée de 1 et de 0, qui va transformer un vecteur { }ix de taille in ( in étant le nombre d’inconnus du domaine iΩ ) en un vecteur de taille n, où n est le nombre d’inconnues du domaine initial global Ω, en complétant par des zéros : { } { } nkavec denoeudunpasestnksi denoeudunestksix xR i iki ki T i ,1'0 =    Ω Ω = . Inversement, la matrice iR va réduire un vecteur { }x de taille n du domaine initial Ω, en un vecteur de taille in , gardant uniquement les indices des nœuds appartenant au domaine iΩ . 2Ωδ Γ 1Γ 2Ω 1Ω Ω 1Ωδ 2Γ Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 143 ∑ = −− = Np i ii T iadd RARP 1 11 (222) b/ Sans recouvrement : Méthode du complément de Schur primale Le domaine initial Ω est divisé en deux sous-domaines { }21 ,ΩΩ sans recouvrement, Figure 53, tel que 21 Ω∪Ω=Ω et 021 /≠Ω∩Ω . soit Ω∩Ω=Γ δδ 1 , l’interface entre les deux sous- domaines { }21 ,ΩΩ . Figure 53 : Décomposition du domaine Ω en deux sous-domaines { }21 ,ΩΩ sans recouvrement La solution u du problème continu global se décompose en ( )ΓΩΩ= uuuu ,, 21 avec 1Ωu (respectivement 2Ωu puis Γu ) représentant la restriction de u à 1Ω (respectivement à 2Ω et à Γ). Ainsi le problème initial continu peut se décomposer en deux sous problèmes locaux sur { }21 ,ΩΩ Γ 2Ωδ Γ 2Ω 1Ω 1Ωδ 144 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction               ×= ×∂= = ×=∇⋅∇− ∂ ∂ = T][0,Gsuruu T][0,G)\O(sur0u Osuruu T][0,Osurg)u(a t u 1,2ipour GOi iOi i 0 OiOi iOii Oi rr (223) auquel il faut ajouter une condition de continuité au passage de l’interface Γ. Dans le cas spécifique de l’équation de diffusion de la chaleur ou bien du modèle électromagnétique, cette condition de continuité s’écrit : )u(an)u(an O222O111 ∇⋅=∇⋅ rrrr sur Γ, (224) avec in r la normal extérieure sur Γ au domaine. Si on suppose Γu connu sur l’interface Γ, alors les systèmes locaux (223) à chaque sous-domaines iΩ peuvent être résolus complètement et ce, de manière indépendante. Le point important à résoudre reste la détermination des valeurs de u sur l’interface Γ en s’appuyant sur la condition de continuité (224). La formulation discrète de la méthode du complément de Schur s’obtient en décomposant notre système matriciel (217) suivant les contributions des deux sous-domaines indépendamment et de l’interface : { } { }fuA = devient           =                     Γ Ω Ω Γ Ω Ω f f f u u u AAA AA AA TT 2 1 2 1 332313 2322 1311 0 0 (225) où l’on a supposé que les nœuds interface étaient numérotés en dernier, iiA étant les restrictions de la matrice système aux sous-domaines iΩ et à l’interface Γ. Les contributions 12A et TA12 s’annulent étant donné qu’il n’existe pas de couplage direct entre les nœuds des domaines 1Ω et 2Ω avec une méthode éléments finis. Les restrictions du vecteur { }u aux domaines 1Ω et 2Ω sont obtenus très facilement, à partir du moment où { }Γu est fixé, grâce aux deux premières lignes du système : Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 145 { } { } { })( 1311111 ΓΩ−Ω −= uAfAu et { } { } { })( 2321222 ΓΩ−Ω −= uAfAu . (226) A partir de la troisième ligne, et en substituant { }1Ωu et { }2Ωu par les relations (226), il vient : { } { }3~fuS =Γ , (227) avec la matrice S définie par : 23 1 222313 1 111333 AAAAAAAS TT −− −−= , (228) et le second membre { }3~f s’écrit : { } { } { } { }21222311111333~ fAAfAAff TT −− −−= . (229) La matrice S est généralement appelée complément de Schur de la matrice 33A dans A. La résolution du système (227) peut s’effectuer soit à partir d’une méthode directe mais qui peut se révéler très coûteuse de part la nature dense de la matrice S , soit par une méthode itérative. Ainsi, généralement, une méthode de type gradient conjugué préconditionné est utilisée. Le problème condensé, mieux conditionné que le problème global, converge plus rapidement. L’une ou l’autre méthode vont nécessiter la résolution des deux sous problèmes locaux, résolution pouvant être effectuée en parallèle. Cette étape est obligatoire pour permettre soit l’assemblage partiel de S (méthode directe), soit l’accès aux produits matrice S - vecteur présents dans un méthode itérative. c/ Discussion sur les méthodes primales appliquées à un problème parabolique Les algorithmes des méthodes de décomposition de domaine primales présentés ici pour une partition du domaine initial en deux sous-domaines, peuvent être utilisés identiquement pour un grand nombre de partitions, seulement la gestion des recouvrements s’avère plus délicate. Notre problème parabolique présente un conditionnement cond(A) borné par )( 2−Ο hdt . L’application des méthodes de décomposition de domaine à un problème de nature parabolique s’en trouve simplifiée sur deux points : • Pour un pas de temps de discrétisation dt assez petit, il n’est pas nécessaire d’utiliser une méthode de sous structuration pour le transit des informations, quasi obligatoire pour les problèmes elliptiques. En effet, dès que le pas de temps de discrétisation temporel devient 148 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction Il s’agit de trouver )(),,( 22121 Ω××∈ ΩΩΩΩ LVVuu λ solutions de : );v,(vsupinf 21) O2 v O1 (v γ γ ΩΩ L . (234) La formulation discrète du Lagrangien par une méthode de discrétisation de Galerkine s’écrit : )T 2O2 T 1O1 T iO T iOiO 2 1i ii T iO2O1,O vR-v(R?)bvvA(v 2 1 );v(v +−= ∑ = γL , (235) où iR est la matrice de restriction à l’interface du domaine iΩ : { } { } n,kavec denoeudunpasn'estksi denoeudunestksix xR ki ki T i 1 0 =     Γ Γ = . De manière équivalente à (234) en formulation discrète, { } { } { } ),,( 21 TTT uu λΩΩ est solution du système d’équations : { } { } { } { } { } { } { } { } 0 2211 2222 1111 =− −= −= ΩΩ Ω Ω uRuR RbuA RbuA T T λ λ (236) En remplaçant dans la troisième ligne les { }iuΩ par leurs expressions définies dans les 2 premières lignes, on arrive finalement au système { } { } { }( ) { } { } { }( ) { } { } { }11112122 22 1 22 11 1 11 bARbARD RbAu RbAu T T −− − Ω − Ω −= −= −= λ λ λ (237) Avec la matrice TT RARRARD 2 1 221 1 11 −− += . Le système est résolu de manière itérative en appliquant par exemple une méthode de type gradient conjugué pour obtenir la solution discrète { }λ et nécessite à chaque itération la résolution des sous problèmes locaux. La méthode FETI présente des bonnes propriétés de convergence et peut être utilisée avec des préconditionneurs moins coûteux que ceux nécessaires aux méthodes primales (C. Farhat et F. X. Roux, [20]). Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 149 2.2 Méthode de partitionnement de domaine Cette méthode consiste à construire localement sur chaque sous-domaines les sous systèmes matriciels locaux. Puis la solution globale est obtenue par une résolution itérative du système matriciel global à l’aide d’un solveur itératif parallèle. La matrice du sys tème global n’est pas explicitement construite. Le système matriciel global doit être suffisamment bien conditionné pour pouvoir être résolu à l’aide d’un solveur itératif. Cette méthode modifie peu la version séquentielle du code d’origine, simplement au lieu de travailler sur le domaine global, elle va être restreinte aux partitions pour l’assemblage des matrices et la construction du second membre du système. L’effort parallèle va se concentrer au niveau du solveur itératif. Les opérations de bases, (produit scalaire, produit matrice- vecteur,..) vont ainsi être modifiées pour permettre une résolution globale. La même séquence de code tourne sur les différents processeurs simplement sur des données différentes : c’est une programmation de type SPMD (Single Program Multiple Data). La méthode de partitionnement de domaine présente deux variantes : • Méthode SPMD par éléments : Un élément appartient à un sous-domaine et à un seul. La méthode par éléments étant celle que nous avons utilisée pour la parallélisation du code éléments finis de simulation des procédés de chauffage, nous renvoyons le lecteur à la publication qui suit pour sa description. • Méthode SPMD par nœuds : En théorie, un nœud appartient à un sous-domaine et à un seul. En pratique cependant, la maillage initial est divisé par élément, comme pour la méthode SPMD par élément, avec un recouvrement d’une bande d’elements (se référer à A. Issman et G. Degrez [26] pour plus de détails). 2.3 Discussion Nous avons vu que les méthodes primales et duales constituent une méthode de résolution efficace pour les problèmes mal conditionnés qui supporte mal un solveur itératif. Néanmoins, elles imposent une condition sur le pas de temps discrétisé pour rester performantes et surtout scalables sans passer par des méthodes de sous structurations, difficiles à mettre en place pour des maillages non structurés. D’autre part, elles nécessitent des préconditionneurs de type 150 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction Schwarz par exemple, souvent coûteux en temps de calcul et en communications entre les processeurs. Concernant les méthodes de partitionnement de domaine, elles sont uniquement envisageables pour des systèmes assez bien conditionnés supportant une résolution itérative. Elles donnent en général de bons résultats. Leur scalabilité va dépendre du choix du préconditionneur. Si un préconditionneur diagonal accélère suffisamment la vitesse de convergence de l’algorithme itératif de résolution, alors la scalabilité sera théoriquement excellente puisque alors le coût d’assemblage est minime et surtout est indépendant du nombre de partitions. Nous avons testé avec succès et donc appliqué un solveur de type gradient conjugué préconditionné à la résolution de nos systèmes paraboliques issus des discrétisations spatiales et temporelles des équations électromagnétiques et thermiques. Une approche de type partitionnement de domaine était alors envisageable pour la parallélisation du modèle. 3 Stratégie de parallélisation SPMD Concernant le type de calculateurs visés, nous nous sommes placés dans le cadre le plus général : machines de type MIMD à grains moyens ou grossiers et à mémoire distribuée. Une architecture à mémoire distribuée est plus difficile à prendre en compte qu’une architecture à mémoire partagée. En effet elle nécessite de répartir les données et les tâches entre les différents processeurs et ce de manière équilibrée afin d’avoir la même charge de travail sur tous les processeurs. Au niveau du partitionneur de maillage, un partitionneur existant, développé par T. Coupez au Centre de Mise en Forme de Matériaux de l’Ecole des Mines de Paris, a été adapté à notre code multimatériaux, nécessitant des développements conséquents. Un maillage initial multimatériau en triangle 6 nœuds obtenu par une méthode de Delaunay est partitionné en Np domaines. Ce maillage possédant une spécificité multimatériau, est plus riche qu’un maillage typique P1 de Forge2. Outre les données élémentaires comme : − le nombre de nœuds − le nombre d’éléments − le nombre de nœuds frontières Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 153 µsfp d 1 = , (238) where f is the frequency, σ the electrical conductivity and µ the magnetic permeability. For instance for a typical steel at 20°C, the skin depth may range from 3mm at a frequency of 50 Hz down to 0.04mm at 105 Hz. It is clear from equation (1) to see why high frequencies are used to achieve surface heating, while low frequencies are favored for initial uniform pre- heating. The accuracy of the simulation of heat treatment or of homogeneous induction heating will depend on the quality of the electromagnetic computations in the skin depth, which can be very small for high frequencies and for a magnetic material. The mesh needs to be extremely refined locally on the surface of the part in order to have enough elements in the width of the skin depth, so as to describe properly the electromagnetic field. This condition imposes the use of large meshes, which can be highly consuming in terms of memory storage and computational time. Parallel computation, by dividing the memory needs between the processors and the reducing the CPU time brings a solution to enable precise modeling of realistic industrial cases. 1. The sequential model Induction heating processes couple electromagnetic induction and thermal diffusion phenomena, [13]. We first present the electromagnetic model and its discretization. The thermal problem and the coupling procedure are then described. 1.1 The electromagnetic problem The electromagnetic model is based on simplified Maxwell’s equations, where the displacement current term tD ∂∂ r has been neglected in the Maxwell-Ampere equation (242). This is known as the magneto-quasi-static approximation. The Maxwell equation system writes Magnetic flux equation 0. =∇ B rr (239) Maxwell-Gauss equation 0. =∇ E rr (240) Maxwell-Faraday equation t B E ∂ ∂ −=×∇ r rr (241) Maxwell-Ampere equation t D JH ∂ ∂ +=×∇ r rrr (242) 154 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction where H is the magnetic field, B the magnetic induction, E the electric field, D the electric flux density, J the electric current density associated with free charges, and x denotes the vector product - thus E rr ×∇ is the curl vector of E . This system of equations needs to be completed by relations which take into account material properties (magnetic permeability µ, electrical conductivity σ, dielectric constant ε). These relations are the three following ones: EeD rr = , (243) HHTB rrr ),(µ= , (244)     = airthein conductorstheinET J 0 )( r r r σ . (245) By combining equations (241) (242) (244) and (245), we get t J E ∂ ∂ −=      ×∇×∇ r r µ 1 . (246) The total current density is decomposed into the sum of the induced current density and of the prescribed one: S J induced JJ rrr += . (247) Finally we get in cylindrical coordinates the electromagnetic model : t SJ r E rr E E t E ∂ ∂ −=      ∂ ∂ −+      ∇∇− ∂ ∂ θ µ θ µθµ θσ 1 2 11 . . (248) The choice of boundary conditions when carrying out a global finite element simulation needs to be carried out carefully. Finite element computation of electromagnetic field in the air can be altered by “artificial reflections” on the boundary of the air domain. We define Ω as the enclosed domain of study, including the part to be heated, the inductor and the air gap. 0Γ is the axis of symmetry such that 0=r and Γ the outer boundary, see Figure 54. A fairly good approximation can be obtained if a Robin boundary condition (249) is prescribed on the enclosing box for the surrounding air. nE ∂∂ θ denotes the derivative of θE along the unit outward normal direction n r and re r is the unit radial direction. A null Dirichlet boundary condition for the θE variable is prescribed on the axis of symmetry. Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 155 Gn.e r E n E r ?? on 0=+ ∂ ∂ rr (249) 0 on 0 GE = θ (250) Figure 54: Presentation of the domain of study Ω and its boundaries. The functional space V in which we are searching the solution is defined by (251) ,where H1(Ω) is a Sobolev space (252) :    = ∂ ∂ Ω∈Ω∈    = 0),(L r ),(HV 21 θ ψψ ψ , (251) { })(L),(L)(H 221 Ω∈ψ∇Ω∈ψ=Ω . (252) We then multiply equation (16) by a test function Ψ belonging to the functional space V and we integrate on the whole domain. After using the Green theorem, we get the weak formulation : V 0 1 1 2 1 ∈∀         + ∂ ∂ ∫ − + ∫ ∂ ∂ −=∫ ∂ ∂ ∫ ++∫ ∇∇+∫ ∂ ∂ ψ?dGn r e r ? E n ? E GG µ ?dv O t S J O ?)dv ? (E rµrO ?dv µr ? E dv O ? E? . µ dv O ? t ? E s rr (253) Chapitre 1 Part Inductor Air Ω Γ Γ Γ 158 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction T being the period of the electromagnetic source term. Taking this average value is justified when one compares the scale of the electromagnetic time step (less than 0.001 s.) with the average induction heating time scale (of the order of seconds). Boundary conditions may be convection and radiation at the interface between the workpiece and the air, (264) where h denotes the convection coefficient, extT the room temperature (in Kelvin), emiε the material emissivity, and Steσ the Stefan constant, or prescribed heat flux or prescribed temperature. Radiation and convection )TT()TT(hn.Tk 4ext 4 Steemiext −σε+−=∇− rr (264) Prescribed heat flux prescribedn.Tk Φ=∇− rr (265) Prescribed temperature prescribedTT = (266) The numerical approach will be the same as for the electromagnetic problem. We consider the domain Ω ΤΗ as being: inductorpartTH Ω∪Ω=Ω . (267) We define the following boundaries: - ΓΤΗ as being the outer boundary of the ΩΤΗ domain, - Γ0 ΤΗ part of the outer boundary where the temperature is prescribed – Dirichlet boundary condition, - Γ1 ΤΗ part of the outer boundary where the heat flux is prescribed – Neumann boundary condition, - Γ2 ΤΗ part of the outer boundary where there is a convection-radiation boundary condition. We need to establish a weak formulation of equation (262) along with the boundary conditions (264)-(266). The functional space V is defined here as:    Γ== ∂ ∂ Ω∈    = TH 0 1 on 0,0 ),(HV ψ θ ψ ψ . (268) After multiplying equation (262) by a test function ψ belonging to the functional space V, integrating on the whole domain, and using the Green theorem, we get the weak formulation : Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 159 V,? ds, .?Thds .?Fdv ?Qds T.?hdv ?T.kdv ? t T ?C ext G2THG1TH prescribed O em G2THOO ∈∀++=+∇∇+ ∂ ∂ ∫∫∫∫∫∫ & rr . (269) We use the same mesh for the part and inductor as the one used for electromagnetic computations. The functional space V is approximated by the discretized space Vh , the test functions ψ by ψh and the unknown T by Th. The discretized version of equation (269) writes : [ ] [ ]{ } { }ththth BtTKt t T C =+       ∂ ∂ )()( , (270) [ ] [ ] { } ∑ ∫∫∫ ∑ ∫∫ ∑ ∫ = ∩∂∩∂ = ∩∂ =         ++=         +∇∇= = nb.elts elt Gelt jext Gelt jprescribed elt iemi th nb.elts elt Gelt ij elt ij ij th nb.elts elt elt ijij th dsNhTdsNFdv.NQB dsNhNdvN.NkK dvN?CNC 1 1 1 21 2 & (271) The time discretization is identical as for the electromagnetic problem. We obtain the following procedure : - Stage 1: system (272) is solved at time t* for T* [ ] [ ] { } { } [ ] { } [ ] 123 1 2 23 2 12 1 tt *th 0 t*th 1 *th**th*th 23 t 1 t c tt 1 t c TCcTCcBTKC t 1 δ −γ − δα γα = δα γα + δ −γ + δ γ =       ++=      + δα γ δ− (272) Matrices are linearized and do only depend on their values at time t and t-δt1 : 160 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction [ ] [ ] [ ] 1 1 2 32 1 1 2 31 1 d ttd tt * C)) dt dt (a(aC) dt dt a(aC −− +++−= , [ ] [ ] [ ] 1 1 2 32 1 1 2 31 1 d ttd tt * K)) dt dt (a(aK) dt dt a(aK −− +++−= . (273) - Stage 2: The temperature field at time t+δt2 is computed: { } { } { } { } )TaTaT( a T tdtt*dtt 2 1 1 3 1 2 −−= −+ . (274) 1.3 Electromagnetic and thermal coupling Due to the difference of time scale between the electromagnetic and thermal phenomena, a strong coupling between both models is not appropriate. At each time step, the electromagnetic field distribution and thus the induced currents depend explicitly on the thermal field in the workpiece, ruled by the heat transfer equation. Inversely, the eddy currents calculated from the electromagnetic solution are used as the heat sources for the thermal finite element analysis. The handling of the coupling is carried out in two different stages: - The criterion from electromagnetic to thermal calculations relying on the stabilization of the mean heat source term over the electromagnetic periods. Once the electromagnetic field has been calculated, the rate of heat generation emQ& for the heat equation needs to be evaluated at every integration point. As the electromagnetic time step is far smaller than the thermal one, we do not consider the instantaneous Joule power calculated at a given time at every integration points. We rather consider a mean Joule power averaged over one or seve ral periods of the electromagnetic field: dttEt T nTQ nT Tn em ∫ − = )1( 2)(int,)(int, 1 int),( θσ , (275) where int is the considered integration point, T is the period of the power supply currents, n is number of periods considered and )(int, tEθ is the value at time t of Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 163 2 ),( 1 bbA− − (279) It is equivalent to minimize J or to minimize E defined by: ))(),(()( 1 xrAxrxEx n −=ℜ∈∀ , (280) where AxxAAxbxr −=−=)( is the residual vector. In order to minimize the functional E, the descent method are constructed by choosing at the kth iteration a descent direction 0≠kp and a scalar kα such that )()( 1 kk xExE <+ with kkkk pxx α+=+1 and 1−+= kkkk prp β . The performance of the iterative method may strongly benefit from a low condition number of the system matrix. In order to lower the condition number of the system matrix, we may use a preconditioner, i.e., a non-singular, symmetric and positive definite matrix C and consider the equivalent system: bCAxC 11 −− = . 164 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction The preconditioned conjugate algorithm writes : Initializations Initial guess x0 Residual r0=b-A x0 Descent direction 0 1 0 rCp −= Auxiliary vector 00 pz = Minimization of J(xk+1) ( )kk kk k ppA zr , ),( =α Upgrade of kkkk pxx α+=+1 Upgrade of kkkk pArr α−=+1 Upgrade of 1 1 1 + − + = kk rCz k=k+1 Evaluation of the new direction ( )kk kk k zr zr , ),( 11 1 ++ + =β Upgrade of kkkk pzp 111 +++ += β Stopping criterion ε≤kr ?? where (..,..) means the scalar product. A compact storage is used for the symmetric matrix. The diagonal preconditioned conjugate gradient shows a relative good behavior with a fast convergence rate for both the electromagnetic and thermal systems. In theory, this saving of time should become even more important with the increase of number of nodes in the mesh. The detail of the conjugate gradient solver shows that only scalar and Matrix-vector products are needed to converge to the solution. These mathematical operations may be easily transformed to local operations by adding of the different contributions at specific locations. This is the foundation of the parallel strategy. Furthermore a diagonal preconditioner is well adapted to a parallel method as it is evaluated locally. NO YES EXIT Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 165 2. SPMD parallel method A large panel of parallel computers are available today, ranging from super computers to station networks. These different types of parallel computers are characterized by different architectures depending on the memory location (distributed or shared or hierarchical), on the number of processor and their computational power, and on the topology of the connections between the processors. An important notion for the choice of an adapted parallel strategy is the notion of granularity defined by the ratio of the processor numbers over their power. Granularity differentiates massive parallel computers (very large number of low power processors) from coarse grain machines (small number of high power processors). The level of parallelization and hence the parallel method strategy will be different depending on the type of computer considered. We have decided to focus here on a “divide and conquer” strategy adapted to coarse grain computers with distributed memory and based on the message-passing paradigm. In that field, several domain decomposition methods are available, among them the alternative Schwarz method, the Schur complement method [11] or the FETI method [20]. Furthermore for well- conditioned global problem that can be solved with an iterative solver, a domain partitioning method can be applied which offers the advantage of being easily portable which is a non- negligible aspect for a commercial code. We have considered this last method to parallelise the finite element code described in the previous section. It consists in separating an initial large mesh in separated non-overlapping sub-meshes. The local sub-system contributions are built and are used within the parallel iterative solver to obtain the global solution fie ld, see Figure 56. We start by describing the mesh partitioning algorithm and its adaptation with our multi-material properties. The parallel method and the preconditioning strategy are then presented. Finally numerical results and performance measurements are shown. 168 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction ii. The greedy algorithm A different color is allocated to each germ. At each iteration, all the non colored elements at a distance equal to 1 to one color, do take this color, and so on until all elements are colored. iii. Optimization of the sub-meshes This step consists of minimizing the interfaces size by changing locally the colors of one element so to minimize a cost function that takes into account the number of elements in each color (load par processor) and the number of edges shared by elements of different colors which represent the interfaces. An element is hence visited, its color is changed or not, in order to optimize the interface’s size. For more details, see [16]. This partitioner acts on a linear mesh. Once the topology of the non-structured sub-meshes are obtained, the linear elements are converted into quadratic triangular elements and the interfaces node local numbering are recomposed along with the multi-material specifications. The main difference of the sub-meshes with a non partitioned mesh concerns the external boundaries which are considered as global only if they are boundaries of the initial mesh. Otherwise they are considered as domain interfaces. For each domain, the number of interfaces is given along with their nodes global numbering from the initial mesh. Figure 58 illustrates an optimized partitioning of a mesh. Figure 58: Splitting of an initial multi-material mesh of 714 elements in three sub-meshes of 240, 236 and 238 elements. Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 169 For each sub-domain, a file describing the interfaces and their nodes is also created. The major aspect to be careful with is that the nodes allocated on a common interface of two sub- domains need to be described in the two corresponding files in exactly the same order. 2.2 SPMD parallel strategy The induction heating code has been parallelized with a SPMD domain partitioning method in connection with the MPI (Message Passing Interface) library. The iterative diagonal preconditioned conjugate gradient solver has been applied with success on the global resolution of both parabolic equations describing the electromagnetic and thermal phenomena. Iterative solvers like the conjugate gradient only use matrix-vector products and scalar products, if we neglect in a first place the preconditioning procedure. Those mathematical operations can be run locally on each partitioned mesh and their scalar or vector values simply need to be actualized or added at the interface nodes. The SPMD (Single Program Multiple Data) approach consists in running the same sequential code on all the different processors, but with different data. Data actualizations are needed at different steps of the code and especially in the matrix system resolution procedure. The domain partitioning method basically remains to parallelized the iterative solver. With respect to the sequential code, two data structures specific to the parallel version were added. The first data structure concerns the number of processors, their name and the correspondence between their number and the sub- domain allocated to each one and inversely. The second data structure handles for each sub- domain - The number of neighbor sub-domains which correspond to the number of interfaces - The number of each neighbor sub-domain - The number of nodes per interface - The numbering of the interface nodes In a detailed manner, a local matrix system with the boundary conditions of the global problem is built on each sub-domains by the allocated processor. At the end of this step, the local matrices and local load vectors are available. The parallel iterative solver algorithm is described with ku being either the electric field or the temperature field at iteration k depending on whether the electromagnetic or the thermal system is being solved. 170 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction Figure 59: The domain partitioning algorithm Developments specific to the parallel code (with respect to the sequential code) are put in italic. We will know describe more particularly the three parallel operations i. Vector field actualization We consider the vector V that needs to be actualized. On the interior and on the global boundaries of each local subdomain, the values of V are identical to the global sequential vector V. However, they are different on the interfaces between two or more subdomains, where the different contributions need to be added. We will illustrate this procedure with the MPI library. - Opening of the reception order on all processors with a non blocking MPI- IRECV procedure. Each processor waits for incoming input vectors from all neighbor processor with the size of the common interface. - For all processors, sending of the local interface values to all the processors which share the same interface. The MPI_ISEND procedure is used. - Synchronization of all processors with the procedure MPI_WAITALL - On all processors, summation on all interface nodes of the local vector values with the received values Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 173 Process parameters Frequency (Hz) 50 Current density 7.0 10-9 Electromagnetic time step (s) T/32 Thermal time step (s) 0.5 Magnetic parameters Steel (part) Copper (coil) air σ (Ω-1 .m-1) 3 106 60 106 0 µrelative 1 1 1 Figure 61: geometry and process parameters The physical properties are given for steel and copper (respectively part and inductor) and the air. In this case, the properties are not temperature dependent. Different meshes have been considered with a number of nodes ranging from 2804 to 30000 nodes. The limitation came from the memory available to run the sequential code which has limited us to meshes lower than approximately 30000 nodes. The meshes have been partitioned into two, four and up to twenty submeshes. For example, Figure 62 shows the ISO values of the effective electric field over a partition into 4 submeshes of the 2804-nodes mesh. The effective electrical field gives a good representation of the joule power dissipated in the part and is defined at a given node as the integration over an electromagnetic period of the square of the real instantaneous electrical field: ∫ + = Tn Tn eff dttnodeET nodeE )1( 2 *),( 1 )( 174 Chapitre 3: Calcul parallèle en modélisation du chauffage par induction Figure 62: Effective electrical field (V.m-1) computations on 4 processors for a mesh of 2804 elements. The code runs a complete electromagnetic computation and 20 thermal time step computations. The corresponding CPU time has been measured for each case. Sequential 2 Processors 4 Processors 10 Processors 20 Processors Number of nodes 2804 nodes 1440 nodes in average 735 nodes in average 290 nodes in average 150 nodes in average CPU Speed Up Efficiency 55.75 s 1 100 % 53.66 s 1.04 52 % 41.1 s 1.35 34% 47.49s 1.17 12% 42.5s 1.32 6.6% Figure 63: Efficiency of the parallel code with the parallel diagonal preconditioner with a mesh of 2804 nodes on the 32 bi-processors cluster. Sequential 2 processors 4 processors 10 processors 20 processors Chapitre 3: Calcul parallèle en modélisation du chauffage par induction 175 Number of nodes 19966 nodes 9990 nodes in average 5000 nodes in average 2000 nodes in average 1000 nodes in average CPU Speed Up Efficiency 1180.75 s 1 100 % 794.16 s 1.5 75 % 391.04 s 3.0 75 % 147.52s 7.54 75 % 97s 12.2 61 % Figure 64 : Efficiency of the parallel code for 20 thermal iterations on a short inductor case on a mesh of 19966 nodes on the 32 bi-processors cluster. Wang case Number of nodes Sequential 29966 nodes 2 processors 15000 nodes in average 4 processors 7900 nodes in average 10 processors 3070 nodes in average 20 processors 1500 nodes in average CPU Speed Up Efficiency 2077.1s 1 100% 1377.26s 1.52 76% 708.2s 2.95 74% 266.2 7.86 79% 95s 21.8 109% Figure 65: Efficiency of the parallel code for 20 thermal iterations on a short inductor case on a mesh of 29966 nodes on the 32 bi-processors cluster. The Speed up and efficiency are quite limited for the small mesh of 2804 nodes, see Figure 63. For a larger mesh, the speed up and efficiency become far more interesting, as showed in Figure 64 and Figure 65 The parallel code will become really useful for large mesh with a number of degree of freedom greater than 20 000. As no remeshing procedures are needed during the computation even with a non static inductor, the sequential part of the code is reduced to the its minimum and consists simply of one processor getting the data file names and sending them to the other processors. Time measurements with a non preconditioned conjugate gradient algorithm was not done as the thermal problem has real difficulties to converge without the diagonal preconditioner due to its bad conditioning number. In Figure 65, the efficiency for 20 processors is extremely good, over 100%. This phenomenon can be explained with the fact that the sequential model with 30000 nodes reaches the limit of the hardware of in term of virtual and local memories, which is quite a disadvantage in term of computational time. For the same initial mesh sizes and for the same number of partition, the parallel performance can be quite different, as it depends in a larger extent on the quality of the mesh and of the submeshes. Chapitre 4: Calcul parallèle en optimisation 179 Chapitre 4 Calcul parallèle pour l’optimisation du procédé de chauffage par induction Ce chapitre concerne la parallélisation d’une procédure d’optimisation automatique des paramètres contrôle des procédés de chauffage par induction mise au point par Yann Favennec [22] sur la base du modèle direct présenté dans le chapitre 2. L’objectif d’un code de simulation numérique est de décrire de manière aussi précise et fiable que possible les phénomènes physiques entrant en compte dans le procédé simulé dans une optique de conception ou bien d’amélioration du procédé. La méthode traditionnelle qui consiste à tester de nombreux jeux de paramètres et à quantifier leurs effets sur les champs finaux peut se révéler fastidieuse, consommatrice de temps, de moyens humains et informatiques. Des méthodes plus récentes sont basées sur une optimisation automatique des paramètres de contrôle. Différents objectifs d’optimisation sont possibles, comme par exemple réduire les temps de chauffe ou bien l’énergie consommée durant le procédé. Dans le cadre du projet européen, nous avons plutôt cherché à améliorer la qualité du chauffage: obtenir un jeu de paramètre de contrôle afin que les évolutions de température soient aussi proches que possibles d’évolutions de température optimales prescrites par l’utilisateur en vue entre autre de contrôler la métallurgie finale du produit. L’algorithme utilisé est un algorithme de gradient conjugué pour lequel les sensibilité sont calculées à partir d’une méthode adjointe discrète. L’algorithme complet d’optimisation nécessite de lancer de nombreuses fois le modèle direct et de résoudre plusieurs fois les systèmes électromagnétiques et thermiques adjoints. Ces opérations sont coûteuses en temps de calculs. La parallélisation de la procédure d’optimisation permet de diminuer considérablement les temps de calcul, notamment pour des domaines d’étude de grandes tailles, qui d’ailleurs seraient inenvisageables en séquentiel pour des raisons de capacité mémoire. 180 Chapitre 4: Calcul parallèle en optimisation Ce chapitre est présenté sous la forme d’une publication, soumise au journal of parallel and distributed computing. Le modèle direct electro-thermique est brièvement décrit, puis la méthode et l’algorithme d’optimisation sont détaillés (se référer également à l’annexe 1). Puis la stratégie de parallélisation, basée sur celle développée pour le modèle direct, est expliquée. Enfin des mesures de speed-up et d’efficacité permettent de mesurer les performances du code parallèle par rapport au code séquentiel en termes de temps de calcul. Chapitre 4: Parallélisation de la procédure d’optimisation 183 t S J r E µ 1 r2r E µ 1 E µ 1 . t E s ∂ ∂ −=      ∂ ∂ −+           ∇∇− ∂ ∂ r on ]t,[t f0 ×Ω , (284) 0t)E(0, = on the axis of symmetry, ]t,[ttandrfor0t),rE( f0 ∈∞→= rr , Ω= onE,0)rE( 0 r , where E is the orthoradial component of the electrical field, σ is the electrical conductivity, µ is the magnetic permeability and SJ r , the input current density (may be of any shape). The weak formulation of equation , (284) is discretized spatially with a Galerkin finite element method [43]. [ ] [ ]{ } { }EEE BE(t)K(t) t E C =+       ∂ ∂ , (285) whith {E} the vector of unknowns {E}={E1, …, En } and: [ ] [ ] { } ds.N t J B ds)N(N rrµ 1 dsNN rµ 1 dsN.N µ 1 K dsNNsC elt i s nb.elts 1elt i E elt ij elt ij2 elt ij nb.elts 1elt ij E nb.elts 1elt elt ijij E ∫∑ ∫∫∫∑ ∑ ∫ ∂ ∂ −= ∂ ∂ ++∇∇= = = = = The time discretization uses a second-order two time step finite difference scheme. Finally we end with a matrix system to solve: ( ) ( ) ]t,[tt 0TE,u,R such that xt,EE Find f0 E ∈∀== , (286) where u are the process parameters (for instance the current density JS). b/ The thermal model We consider the domain Ω with the following boundaries: Γ as being the outer boundary of the Ω domain, Γ0 part of the outer boundary where the temperature is prescribed (Dirichlet boundary condition), Γ1 part of the outer boundary where the heat flux is prescribed 184 Chapitre 4: Parallélisation de la procédure d’optimisation (Neumann boundary condition) and Γ2 part of the outer boundary where there is a convection-radiation boundary condition. Temperature evolution in the workpiece is governed by the classical heat transfer equation: em QT)div(k t T ?C & r =∇− ∂ ∂ on ]t,[t f0 ×Ω , (287) prescribedTT = on Γ0 , prescribedn.Tk Φ=∇− rr on Γ1 , )T(Tse)T(ThnT.k 4 ext 4 Steemiextcv −+−=∇− rr on Γ2 , where ρ is the workpiece density, C the specific heat, k the thermal conductivity, (both C and k being temperature-dependent). hcv denotes the convection coefficient, extT the room temperature (in Kelvin), emiε the material emissivity, and Steσ the Stefan constant. emQ& is the local heat density rate generated by the Joule losses due to eddy currents: 2EQ em σ=& with dt(t)Es(t) P 1 Es 1)P(n nP 22 ∫ + = , with P being the period of the electromagnetic source term. Going through the weak formulation and the spatial discretisation, one ends with the following linear matrix system to solve [ ] [ ]{ } { }TTT BT(t)K(t) t T C =+       ∂ ∂ , (288) with [ ] [ ] { } ∫∫∫∑ ∫∫∑ ∑ ∫ ∩∂∩∂= ∩∂= = ++= +∇∇= = 21 2 Gelt jext Gelt jprescribed elt iem nb.elts 1elt i T Gelt ij elt ij nb.elts 1elt ij T nb.elts 1elt elt ijij T NhTNF.NQB NhNN.NkK N?CNC & Chapitre 4: Parallélisation de la procédure d’optimisation 185 The time discretization uses a second-order two time step finite difference scheme. We finally end with a matrix system to solve: ( ) ( )    ∈∀== f0 T t,tt0TE,R:thatsuchxt,TTFind . (289) c/ The electro-thermal coupling procedure The electromagnetic and thermal model are tightly coupled. The handling of the coupling is carried out by checking two different stages: ♦ Stabilization of the mean heat source term over the electromagnetic periods. The mean Joule heating rate is computed over one period P of the electromagnetic field at every integration points: dtE(t)s(t) P 1 (nP)Q nP 1)P(n 2 em ∫ − = . Electromagnetic computations are run until the heat source term converges to a stabilized value: e (nP)Q (nP)Q1)P)((nQ em emem < −+ . ♦ Stabilizations of the magnetic and electric parameters values with respect to temperature : Thermal computations can use the same source term derived from the electromagnetic computations as long as the physical magnetic parameters such as the magnetic permeability and the electric conductivity do not change “too much”. Their variations with temperature are tested after each new thermal computation. The following criteria are used for each mesh element: andn nn %5 )( )()( max max 1 max < Τ Τ−Τ + σ σσ %5 )( )()( max max 1 max < Τ Τ−Τ + n nn µ µµ , where 1max +Τn is the maximum value of the temperature field in the mesh element at time therdtt + and n maxΤ is the maximum value of the temperature field in the same element at the current time t . 188 Chapitre 4: Parallélisation de la procédure d’optimisation ( ) ( ) ( ) 2 O optTtx,TThTg ∫ −== , (293) where the operator . is a norm. Hence, the optimization problem writes: ( ) ( ) ( ) ( )ujminujsuch that Uu Find 0TE,R 0TE,u,R ,Uuad T E ad = = ∈ =∈ . (294) Electromagnetic control parameters u are the prescribed exciting frequency and the current input circulating within the coils. The velocity of inductors with respect to the billet to be heated is also a parameter that can be controlled. The control space adU is then to be defined. This depends on the kind of the optimization problem dealt with. For the real optimization case dealt with in the last section of the paper, the proper control space adU is given. Several minimization methods can be used to solve the optimization problem (294). They are detailed and compared in [22]. In the general case, when the functional j is differentiable, a necessary condition for u to be solution of (294) is given by: ( ) ( ) adUuuuuj ∈∀≥− 0.' . (295) The determination of J∇ r , gradient of J with respect to all controls kuk ∀, enables to use gradient type minimization methods. The optimal control approach is used here. The first step is to introduce the Lagrangian L given in (296). ( ) ( ) ( ) [ ] ( ) [ ]ff tt TT tt EETE TERTEuRTJTEuL ,, 00 ,,,,,,,,, ×Ω×Ω ++= λλλλ , (296) where Eλ and Tλ are adjoint variables related to electromagnetic and thermal problems, and where: [ ] ∫ ∫ Ω ×Ω = f f t t tt dtdxuvvu 0 0 , , . (297) The second step is to derive the Lagrangian with respect to controls ku and integrate by parts all scalar products. It can be shown that if Tλ and Eλ are solution of the reverse adjoint problem: Chapitre 4: Parallélisation de la procédure d’optimisation 189             =        ∈∀ ∂ ∂ −=+ ∂ ∂ −       −=        ∈∀      −=+ ∂ ∂ − ∫∑ ∫∑ = = 0t?C t,tt? E B ?K t ? C NTTt?C t,ttNTT?K t ? C f EE f0 T T EE E E i elt opt nbelts 1eltf TT f0i elt opt nbelts 1elt TT T T (298) the objective function gradients k du dj are given by: ( ) [ ]ftt E k E k u TEuR du dj ,0 , ,, ×Ω ∂ ∂ = λ (299) When dealing with frequency and current amplitude as control parameters, the first term of the scalar product involved in (299) is calculated deriving analytically the direct electromagnetic problem. The choice for the perturbation is not trivial. Experience has shown that the perturbation depends, among other things, on the element size of the finite element mesh. As for the direct induction heating model, a time integration of the adjoint problem (298) within the whole time range [ ]ftt ,0 is not practically possible. Therefore a weak coupling strategy has been developed (see [22] for more details). Once the objective function gradient has been calculated, the determination of u is done using a gradient type method. A first choice Ou being given, one builds the series defined as: nnnn Duu α+= −1 , (300) where nD is the descent direction calculated using, in the present case, the conjugate gradient method, and nα is the descent step defined by: ( )nnn Duj αα += −1minarg . (301) The global coupling procedure, encompassing the direct and the reverse model is schematically shown in Figure 70. 190 Chapitre 4: Parallélisation de la procédure d’optimisation Input physical data, geometry… Initialization of E,H, B, T E.M. computation Period completed? Convergence? Thermal computation Finished? End direct model Initialization ofλE, λT Derivation of costfunc. wrt temp Thermal adjoint computation E.M. adjoint computation Period completed? Convergence? Finished? End adjoint model Parabolic interpolation research algorithm Control stabilized? End Initial guessed controls NO NO NO NONO NO Figure 70: General algorithm for the optimization procedure coupled to the induction heating model. 3. The parallel strategy The method employed to parallelize the global optimization procedure is based on a domain partitioning method which has been already successfully applied for the parallelization of the direct induction heating model, [31]. The global SPMD (Single Program Multiple Data) approach based on a domain partitioning method tends to solve directly the global linear systems with the contribution of all processors. This strategy is only usable if the different matrix systems constructed during the Chapitre 4: Parallélisation de la procédure d’optimisation 193 Figure 71: The domain partitioning algorithm The parallelization of the operations, such as scalar or matrix-vector products, are detailed in [31]: 1. Vector field actualization We consider the vector V that needs to be actualized. On the interior and on the global boundaries of each local subdomain, the values of V are identical to the global sequential vector V. However, they are different on the interfaces between two or more subdomains, where the different contributions need to be added. We will illustrate this procedure with the MPI library. - Opening of the reception order on all processors with a non blocking MPI- IRECV procedure. Each processor waits for incoming input vectors from all neighbor processor with the size of the common interface. - For all processors, sending of the local interface values to all the processors which share the same interface. The MPI_ISEND procedure is used. 194 Chapitre 4: Parallélisation de la procédure d’optimisation - Synchronization of all processors with the procedure MPI_WAITALL - On all processors, summation on all interface nodes of the local vector values with the received values 2. Scalar product (x, y)Ω The scalar product of two global vectors are computed locally on each sub- domain. The different contributions are then added and the global scalar product value is sent back to all processors. However, this way, the contribution of one interface node is added twice or more when the considered node is shared by two sub-domains or more. To avoid this, a weighting factor is introduced - for an internal node, it is equal to one - for an interface node, it is equal to the inverse of the number of domains owning the considered node Hence the scalar products between two global vector x and y writes: ∑ ∑ = = ΩΩΩ         = ΩPROCESSOR iNODEN i N node ii nodeweightnodeynodexyx 1 1 )(.)(.)(),( , (302) with     = )(__ 1 1 )( nodedomainofnumber nodeweight For an internal node for an interface node (303) 3. Matrix-vector product Ax The global matrix A is the sum of the local matrices iAΩ (extended to the size of the global problem). On the interface nodes, the value of the global matrix A corresponds to the sum of the contributions of the local values on all domains which share the given interface node. Hence the local matrix-vector product is done on all sub-domains, followed by an actualization of the resulting vector at the interface nodes. Chapitre 4: Parallélisation de la procédure d’optimisation 195 The general organization of the parallel optimization algorithm is described in Figure 72. The parallel tasks are in blue. Figure 72 : General organization of the parallel optimization algorithm. Direct model Direct model Direct model Iterative parallel solver Computation of the cost function Adjoint models Adjoint models Adjoint models Iterative parallel solver Computation of the cost function gradient Descent direction Descent direction Descent direction Linear research algorithm Direct model Direct model Direct model Iterative parallel solver 198 Chapitre 4: Parallélisation de la procédure d’optimisation 35 resolutions of the direct model and 5 resolutions of the adjoint model have been necessary to determine the optimal current densities. The cost function has decreased by a factor of 60 between the first and the last iterations. We see that the current densities of the eight center coils are quite identical. However the two extreme coils do have quite different current densities, which have been only slightly modified during the iterations. The explanation is that they are further away and thus their contributions in term of heating rate are less significant. Hence the derivative of the cost function with respect to their current densities are small. Figure 75 : ISO curves of the effective electrical field on a partition into 8 sub-meshes from the initial mesh of size 70000 nodes. Figure 75 shows the iso values of the effective electrical field defined by: ∫ + = T1)(n Tn 2 eff dt*(t)E T 1E As the surface of the part and the inductors were meshed with a very refined mesh, the partitions are really inhomogeneous in size thus having the number of finite elements. Chapitre 4: Parallélisation de la procédure d’optimisation 199 Figure 76 : ISO curves of temperature field on a partition into 8 sub-meshes from the initial mesh of size 70000 nodes. The colored regions are parts of the gas container. Figure 76 shows the temperature field. This time the relative size of the different submeshes are not real. We can see that the temperatures are quite homogenous in the part. The test described above has been computed on a mesh of 70000 nodes for precision purposes; this mesh has been partitioned into eight meshes of size slightly less than 10000 nodes, see Figure 75. The total computation has lasted more than 24 hours. Performance measurements were run on meshes of about 20000 and 28000 nodes to enable comparisons in term of CPU time with the sequential code, as the sequential code could not run with meshes larger than 28000 nodes. We present results in term of speed up obtained on parallel computers with a distributed architecture. The parallel platform available at our laboratory and which we have used for CPU measurements is a cluster of 32 bi-processors Pentium III at 1GHz and 512 Mo of RAM. The data transfer velocity is 2 Go/s. 200 Chapitre 4: Parallélisation de la procédure d’optimisation The measurement have been made only over one complete optimization iteration, which has runs both the electromagnetic and thermal adjoint computations and nine time the direct electromagnetic and thermal models. Sequential 2 Processors 4 Processors 10 Processors 20 Processors Number of nodes 19834 9980 in average 4990 in average 2000 in average 1000 in average CPU(s) Speed Up Efficiency 22374.82 1 100% 16282.55 1.4 70% 10829.1 2.1 52% 5802.9 3.9 39% 3152.2 7.1 35% Figure 77: Efficiency of the parallel optimization code with a mesh of 19834 nodes on the 32 bi-processors cluster. 27784 nodes Sequential 2 Processors 4 Processors 10 Processors 20 Processors CPU(s) Speed Up Efficiency 38131.54 1 100% 26929.09 1.42 71% 15757.0 2.42 60% 6920.3 5.51 55% 3963.8 9.62 48% Figure 78: Efficiency of the parallel optimization code with a mesh of 27784 nodes on the 32 bi-processors cluster. The performances in terms of computational time restriction are good. Yet, the efficiency is rather disappointing when the number of processors increases. This is due to the fact that a computation of the electromagnetic or thermal adjoint model requires to store all the nodal values at all the direct model iterations of the corresponding field. Hence those field values are written in data files after each iteration of the direct model and then read once before the adjoint model resolution. The length of these files is proportional to the number of nodes of the considered mesh. Those files become smaller when the number of submeshes increases. Although the time needed for these operations gets largely reduced when the number of partitions increases, the time consumption due to the opening and closing of the files added to the time needed to initialize writing and reading procedures remains the same and is Conclusion générale et perspectives 203 Conclusion générale et perspectives La première étape de ce travail a consisté à établir, développer et valider un modèle performant pour modéliser les procédés de chauffage par induction, que ce soit en préchauffe ou pour des traitements thermiques. Ce procédé est complexe de par sa nature multi-physique et nécessite le couplage entre des modèles : - électromagnétique, - thermique, - éventuellement thermo-mécanique. Le choix du modèle électromagnétique est primordial. De nombreuses approximations basées sur des hypothèses plus ou moins fortes existent. Nous avons seulement utilisé l’approximation des régimes quasi-permanents. Nous avons vu que cette première approximation, qui revient à négliger le phénomène de propagation des ondes, est valable dans la gamme de fréquences utilisée lors des procédés de chauffage par induction, les plus hautes fréquences étant largement inférieures au mégahertz. La propagation des ondes est alors considérée comme instantanée, ce qui au vu de la taille caractéristique des installations (quelques mètres) par rapport à la célérité de la lumière (3.105 m/s) est tout à fait raisonnable. En revanche, nous avons choisi d’écarter l’approximation harmonique des champs électromagnétiques. Cette approximation découple les évolutions spatiales et temporelles du champ et revient à calculer une amplitude complexe pour le champ électromagnétique à partir d’une équation stationnaire. L’avantage d’une telle approximation est le gain souvent important en temps de calcul. Seulement, on perd une précision importante sur l’évolution temporelle et sur la déformation des champs électromagnétiques lorsqu’il s’agit d’un matériau ferromagnétique. En effet, les harmoniques secondaires ne sont pas prises en compte. Afin de pouvoir représenter les phénomènes physiques le plus réellement possible, le modèle électromagnétique utilisé est dépendant du temps. Néanmoins, afin de n’être pas trop pénalisant en temps de calcul, des compromis entre la précision des calculs et le temps de calcul nécessaire ont été étudiés. Ils se situent au niveau : - du nombre de calculs électromagnétiques nécessaires pour bien décrire l’évolution temporelle d’une période électromagnétique, 204 Conclusion générale et perspectives - du nombre de périodes électromagnétiques nécessaires pour arriver à une solution stable, - du nombre de calculs électromagnétiques complets nécessaires au cours de l’évolution du champ de température. Ces points importants, ainsi que des échelles de temps caractéristiques électromagnétiques et thermiques présentant un rapport allant de 10-2 à 10-6 ont nécessité la mise en place d’un couplage faible, basé sur la stabilisation du terme de puissance Joule moyennée sur une période électromagnétique ainsi que sur la stabilisation des paramètres électromagnétiques au cours de la montée en température. La méthode numérique employée, de type éléments finis, est fiable et robuste. Néanmoins, elle nécessite une bonne compréhension des phénomènes physiques électromagnétiques inhérents au procédé. En effet, modéliser un espace ouvert par une méthode éléments finis nécessite la fermeture du domaine et l’imposition de conditions aux limites artificielles. L’utilisateur doit estimer la taille du domaine étudié qui doit être assez grand pour ne pas venir tronquer les lignes du champ électromagnétique et ainsi les modifier. Son avantage par rapport à une méthode mixte est que la matrice du système est creuse et symétrique. La résolution du problème est facilitée et se prête mieux à des développements en calcul parallèle. Enfin, une nouvelle stratégie a été développée pour simuler le déplacement de l’inducteur: ses propriétés se déplacent virtuellement dans l’air. Cette méthode a donné de très bons résultats et ne nécessite aucun remaillage. Les perspectives de recherche sont multiples. Au niveau des données, le modèle accepte actuellement une tension ou une dens ité de courant source uniforme dans l’inducteur. Suite à un calcul électromagnétique complet, la répartition de courants est connue dans l’inducteur et permet une évaluation de l’intensité réelle circulant dans les spires. Il serait intéressant de mettre au point un outil de transfert des données électrotechniques vers nos paramètres d’entrées. Un autre point, plus académique, serait d’effectuer des comparaisons pour des matériaux ferromagnétiques entre un modèle harmonique et le nôtre, dépendant en temps. En effet nous avons vu que ces deux modèles donnent des solutions identiques pour des matériaux amagnétiques. Tout l’intérêt de notre modèle dépendant en temps apparaît par son analyse beaucoup plus riche des matériaux non linéaires. Nous avons vu que le signal périodique peut être grandement déformé et ne ressemble alors plus du tout à une sinusoïde. Néanmoins, il Conclusion générale et perspectives 205 n’est pas forcément évident que la puissance Joule, issue du calcul électromagnétique et obtenue par intégration sur une période électromagnétique, soit très différente de celle obtenue par une analyse harmonique. Cette différence serait très intéressante à quantifier. Enfin des comparaisons entre les méthodes numériques ‘tout’ éléments finis et mixtes permettraient de quantifier la précision des méthodes suivant les tailles des éléments finis, les tailles du domaine de fermeture, ainsi que les différences en temps de calculs. Un autre axe de ce travail a consisté à étudier et à implémenter une stratégie de parallélisation du modèle direct et de la procédure d’optimisation. Nous avons commencé par tester des solveurs itératifs préconditionnés sur nos différents modèles de type parabolique. Ceux ci donnant des résultats satisfaisants par rapport notamment à un solveur direct, nous avons pu nous orienter vers une méthode de parallélisation SPMD de type partitionnement de domaine. Cette méthode, simple et efficace, donne de très bons résultats au niveau du modèle direct, avec une bonne efficacité et une bonne scalabilité. La parallélisation de l’optimisation montre une efficacité convenable sur deux et quatre processeurs mais qui tend à chuter rapidement avec le nombre de processeurs: la scalabilité est relativement moyenne. Ce problème fait apparaître une thématique de recherche intéressante en calcul parallèle appliqué aux méthodes adjointes: améliorer la scalabilité de l’optimisation parallèle en développant une meilleure stratégie d’accès aux données, en rééquilibrant les données stockées et les données à recalculer. Enfin les perspectives à plus long terme consisteraient à développer un modèle analogue tridimensionnel. 208 Références bibliographiques [10] Chaboudez C., Clain S., Glardon R., Mari D., Rappaz J., Swierkosz M., Numerical modeling in induction heating for axisymmetric geometries, IEEE Trans. Magn., Vol. 33, No. 1, pp. 739-745, Jan. 1997. [11] Chan T. F., Mathew T. P., Domain decomposition algorithms, Acta Numerica, pp.61- 143, 1994. [12] Clain S., Rappaz J., Swierkosz M., Touzani R., Numerical modelling of induction heating for two-dimensional geometries, Mathematical models and methods in applied sciences, Vol. 3, No. 6, pp. 805-822, 1993. [13] Davies E. J., Conduction and Induction heating, London:P. Peregrinus Ltd, 1990. [14] Digonnet H., Repartitionnement dynamique, mailleur parallèle et leurs applications à la mise en forme des matériaux, Thèse de l’Ecole des Mines de Paris, 2002. [15] Dawson C. N. and Dudd Q., A domain decomposition method for parabolic equations based on finite elements, in 4th Int Symp. On Domain Decomposition Methods for partial differential Equations, SIAM, Philadelphia, PA, 1991. [16] DRAMA Consortium, Final DRAMA Cost Model. Deliverable D1.1b, DRAMA Project, 1999. [17] Egan L. R., Furlani E. P., A computer simulation of an induction heating system, IEEE Trans. Magn., Vol. 27, No. 5, pp. 4343-4354, 1991. [18] Evans D. J. , The use of preconditioning in iterative methods for solving linear equations with symmetric positive definite matrices, J. Inst. Maths Applics, Vol 4, pp 295-314. [19] Farhat C.et Roux F. X., A method of finite element tearing and interconnecting and its parallel solution algorithm, J. Comput. Systems Eng., 2 :149-156, 1991. Références bibliographiques 209 [20] C. Farhat et F. X. Roux, Implicit parallel processing in structural mechanics, Comp. Mechan. Advances, 2 :1-124, 1994. [21] Favennec Y., Labbé V., Bay F., Induction heating processes optimisation – A General Optimal Control Approach, submitted to Journal of comput. Phys.l, 2001. [22] Favennec Y., Modelisation numerique en chauffage par induction. Analyse inverse et optimisation, Thèse de Doctorat, Ecole Nationale Supérieure des Mines de Paris, 2002. [23] Gie H., J. P. Sarmant, Electromagnetisme, Volumes 1 et 2, collection des sciences physiques, Technique et documentation, Lavoisier, 1985. [24] Guerin C., Détermination des pertes par courants de Foucault dans les cuves de transformateur, Thèse de Docteur-Ingénieur, INP Grenoble, 9/1994. [25] Hogge M. A., A comparison of two and three level integration schemes for non-linear heat conduction, Numerical Methods in heat transfer, Ed. R. W. Lewis, K. Morgan, O. C. Zienkiewicz, John Wiley and sons, pp. 75-90, 1981. [26] Issman E. and Degrez G., Non overlapping preconditioners for a parallel implicit Navier Stokes solver, Future Generation computer Systems, 13: 303-313, 1998. [27] Jufer M., Radulescu M., Inclusion of magnetic hysteresis losses and eddy current dynamic losses in the modelling of induction heating for axisymmetric geometries, International report LEME, Ecole Polytechnique Fédérale de Lausanne, 1994. [28] El-Kaddah N., Craen R., Loue W.,A mathematical model of induction heating of thixoformable aluminium billets, Light metals 1998: Proc. of the technical sessions of the 127th TMS annual meeting, San Antonio, Feb 15-19,1998. 210 Références bibliographiques [29] Krähenbühl N., Nicolas A., Méthode des équations intégrales de frontières, développement des techniques et formulation axisymétriques, RGE March, pp. 222-226, 1985. [30] Kuznetsov Y. A., New algorithms for approximate realization of implicit difference schemes, Sov. J. Numer. Anal.Math. Modell. 3, 99-114, 1988. [31] Labbé V., Favennec Y., Bay F. , Numerical modeling of heat treatment and induction heating processes using a SPMD parallel computational method, Soumis à Engineering Computations, Int. J. Computer-Aided Engineering and Software, Février 2002. [32] Labridis D., Dokopoulos P., Calculation of eddy current losses in nonlinear ferromagnetic materials., IEEE Trans. Magn., Vol. 25, No. 3, pp. 2665-2669, 1989. [33] Marie S., Un modèle de parallélisation SPMD pour la simulation numérique de procédés de mise en forme des matériaux, Thèse de Doctorat, Ecole Nationale Supérieure des Mines de Paris, 1997. [34] Ter Maten E.J.W., Melissen J.B.M., Simulation of inductive heating., IEEE Trans. Magn., vol. 28, No.2, pp. 1277-1290, 1992. [35] Marchand C., Foggia A., 2-D finite element program for magnetic induction heating., IEEE Trans. Magn., Vol. 19, No. 6, pp. 2647-2649, 1983. [36] Masse Ph., Morel B., Breville T., A finite element prediction correction scheme for magnetothermal coupled problem during Curie transition, IEEE Trans. Magn., Vol. 21, No. 5, pp. 181-183, 1985. [37] Meijerik J.A., Van Der Vost H., An iterative solution method for linear systems of which the matrix is a symmetric matrix, Math. Comp., Vol 31, n° 137, pp.148-162.