Calcul des champs produits par des courants variant dans le temps ou des sources statiques en mouvement.
Sous-sections de MEF en Électromag. BF
Avant-propos
Maxwell House
Les éléments finis de Whitney développés à la fin du chapitre 2 nous fournissent toute une structure permettant d’accueillir parfaitement toutes les constituantes du champ électromagnétique en basse fréquence.
Une figure valant mieux qu’un long discours, vous trouverez ci-dessous une version de la Maison de Maxwell (diagramme de Tonti 3D appliqué à l’électromagnétisme dans l’ARQS).
La Maison de Maxwell (Maxwell House) dans l'ARQS : les différents étages de chaque côté correspondent à un espace fonctionnel (en couleur), on passe d'un étage à l'autre par un opérateur différentiel (flèches noires), de l'avant à l'arrière par dérivation temporelle (flèches noires aussi), et d'un côté à l'autre par les lois de comportement (flèches grises).
Pour vous, outre l’aspect esthétique, le principal intérêt de ce diagramme est de rapidement trouver à quel espace appartient chaque champ scalaire ou vectoriel.
$\rightsquigarrow$ Dans la suite, nous allons voir comment formuler des problèmes simples du domaine avec la méthode développée précédemment.
Électrostatique
Dans le cas de l’électrostatique (en anglais electrostatics), nous nous intéresserons au calcul de la distribution du champ électrique ${\bf e}$ due à des charges statiques ou à différents niveaux de potentiel électrique.
Les relations intervenant en électrostatique sont :
La loi de comportement dans les milieux diélectriques :
$${\bf d} = \varepsilon\,{\bf e}$$
Pour les conditions aux limites, c’est-à-dire aux frontières du domaine d’étude (${\partial\Omega = \Gamma = \Gamma_d\cup \Gamma_n}$), nous aurons :
soit, d’après la condition sur la composante tangentielle de ${\bf e}$ : $$\left.{\bf n}\wedge{\bf e}\right|_{\Gamma_d} = {\bf 0}$$
soit, d’après la condition sur la composante normale de ${\bf d}$ : $$\left.{\bf n}\cdot{\bf d}\right|_{\Gamma_n} = 0$$
Formulation forte
Maxwell-Faraday en statique nous permet de définir le potentiel scalaire électrique ${\bf v}$ tel que : ${{\bf e} = -{\bf grad}\,v}$
Ainsi, dans le domaine d’étude $\Omega$, l’équation satisfaite par $v$ est obtenue par substitution dans l’équation de Maxwell-Gauss :
$$ \text{div}\left(\varepsilon\,{\bf grad}\,v \right) + \rho_q = 0$$
Important !
Attention : En électrostatique, les conducteurs sont supposés « parfaits », le champ électrique est nul à l’intérieur et donc $v$ y est constant.
En pratique, les conducteurs seront donc exclus du domaine d’étude et seules les conditions sur leurs bords $(\Gamma_k)$ interviendront :
soit en tant que condition de Dirichlet où : $$v|_{\Gamma_i} = v_i$$
soit en tant que condition de Neumann où : $$\left.\frac{\partial\,v}{\partial {\bf n}}\right|_{\Gamma_j} = {\bf grad}\,v\cdot {\bf n}\big|_{\Gamma_j} = -\frac{\sigma_{q_j}}{\varepsilon}$$
Aux frontières, nous aurons les mêmes conditions mais homogènes : $ v|_{\Gamma_d} = 0~$ ou $~{\bf grad}\,v\cdot {\bf n}\big|_{\Gamma_n} = 0$.
Finalement, la forme complète de la formulation forte à résoudre est donc :
Nous disposons de tous les outils pour résoudre numériquement l’exemple classique du condensateur cylindrique du premier chapitre (cf ici).
Compte tenu des symétries, le problème pourra se résoudre en 2D axisymétrique sur une demie-géométrie, telle que représentée sur la figure ci-dessous montrant également les conditions aux limites associées ainsi qu’un exemple de maillage :
Dans ce cas, la formulation faible à résoudre se ramène à :
Trouver $v \in \text{H}_{01}({\bf grad},\Omega) = \{ u \in \text{H}({\bf grad},\Omega) : u|_{\Gamma_d} = 0,~ u|_{\Gamma_1} = V_1 \}$, tel que :
On définit tout d’abord les contraintes sur le potentiel $v$.
Celles-ci permettent de définir ensuite l’espace fonctionnel dans lequel se trouve $v$. On l’appelle ici Hgrad (correspond au $W^0$ du chapitre précédent) qui approxime $\text{H}_{01}({\bf grad},\Omega)$. On utilise des éléments nodaux avec les fonctions de base $s_n$ vues précédemment.
Enfin, nous définissons la forme faible du problème.
Après résolution, on obtient la distribution du potentiel scalaire électrique :
Et on peut calculer la valeur de la capacité à partir de l’énergie électrostatique.
Vous pouvez télécharger les fichiers complets encliquant ici.
Condensateurs plans
Adapter l’exemple précédent pour résoudre le problème similaire du condensateur plan à armatures rectangulaires :
modifier la géométrie ;
changer de jacobien ;
adapter les formules de calcul de l’énergie et de capacité.
Étudier également le cas d’un condensateur plan mais, cette fois, à armatures circulaires.
Potentiel Flottant
Dans des problèmes tels que définis ci-dessus, nous n’avons pas directement accès à la charge portée par les armatures (il faut nécessairement passer par des calculs post-processeur). Nous ne pouvons pas non plus proprement considérer cette charge comme une entrée du problème pour calculer la différence de potentiels associée.
Nous pourrions penser à utiliser une condition de Neumann non homogène sur l’armature portant le potentiel $V_1$ mais la valeur du potentiel résultant du calcul ne serait plus constante le long de l’armature.
Pour pallier cet inconvénient, les auteurs de GetDP ont implanté une fonctionalité très pratique : les Global quantities (grandeurs globales) qui permettent de définir ce qu’on appelle un potentiel flottant (floating potential) via un terme global (global term) dans la formulation.
Plutôt que de faire une explication un peu trop appoximative, je préfère vous renvoyer à l’exemple du wiki de ONELAB détaillant son principe de fonctionnement : Electrostatics with floating potentials :
Télécharger les fichiers et lire attentivement les commentaires du .pro (oui, ça relève de la magie !).
Adapter les exemples précédents pour les résoudre en potentiel flottant.
Comme je ne suis pas méchant, je vous donne la solution pour le condensateur cylindrique en cliquant ici.
Électrocinétique
Dans le cas de l’électrocinétique (en anglais electrokinetics ou current flow), nous nous intéresserons au calcul de la distribution de la densité volumique statique de courrant électrique ${\bf j}$ dans des conducteurs.
Les relations intervenant en électrocinétique sont :
La loi de comportement dans les milieux conducteurs électriques :
$${\bf j} = \sigma\,{\bf e}$$
Pour les conditions au limites, c’est-à-dire aux frontières du domaine d’étude $\Omega_c$, nous aurons :
soit, d’après la condition sur la composante tangentielle de ${\bf e}$ sur les faces d’entrée ou sortie du courant $(\Gamma_d)$ : $$\left.{\bf n}\wedge{\bf e}\right|_{\Gamma_d} = {\bf 0}$$
soit, d’après la condition sur la composante normale de ${\bf j}$ sur les autres bords du conducteurs $(\Gamma_n)$ : $$\left.{\bf n}\cdot{\bf j}\right|_{\Gamma_n} = 0$$
Formulation forte
Comme en électrostatique, nous pouvons définir le potentiel scalaire électrique $v$ tel que : ${{\bf e} = -{\bf grad}\,v}$.
En réinjectant cette relation couplée à la loi de comportement dans la loi des nœuds locale, on obtient :
Sur les faces d’entrée et de sortie du courant (${\Gamma_d = \Gamma_{di}\cup\Gamma_{dj}}$), deux cas seront possibles :
Soit une condition de potentiel imposé par une condition de Dirichlet sur $\Gamma_{di}$ de type : $$v|_{\Gamma_{di}} = v_i$$
Soit une condition de Neumann non-homogène imposant la densité de courant normale à la surface $\Gamma_{dj}$ : $$\left.\frac{\partial\,v}{\partial {\bf n}}\right|_{\Gamma_{dj}} = {\bf grad}\,v\cdot {\bf n}\big|_{\Gamma_{dj}} = \pm \frac{j_{n_j}}{\sigma}$$
Sur les autres bords du domaine $(\Gamma_n)$, nous aurons des conditions de Neumann homogènes : ${{\bf grad}\,v\cdot {\bf n}\big|_{\Gamma_{n}} = 0}$
Finalement, la forme complète de la formulation forte à résoudre est donc :
On considère une résistance de sole de four électrique telle que représentée ci-dessous :
On désire calculer la distribution de la densité de courant à l’intérieur lorsqu’on applique une tension de 230 V à ses bornes, afin d’en déduire sa résistance via les pertes Joule associées.
L’implantation dans GetDP ressemble beaucoup au cas précédent :
Après calcul, on obtient pour la répartition du potentiel scalaire électrique :
Vous pouvez télécharger les fichiers complets encliquant ici.
Pour la petite histoire, voici une photo de celle de mon four perso après démontage pour changement suite à un défaut d’isolement :
Résistance de charge
Traiter le problème analogue sur une forme correspondant aux résistances de charge utilisées sur nos bancs expérimentaux :
En dimensionner une permettant de dissiper 5 kW sous 48 V continus. On pourra utiliser comme conducteur du Constantan, par exemple du 45Ni-55Cu de résistivité électrique : $\rho = 14,9\cdot 10^{-7}~\Omega\cdot\text{m}$
Astuce
Pour la géométrie, on pensera à utiliser une extrusion combinant à la fois une translation selon le vecteur $(v_{t_x},v_{t_y},v_{t_z})$ et une rotation d’un angle $\alpha$ autour de l’axe de vecteur directeur $(v_{a_x},v_{a_y},v_{a_z})$ et passant par le point $(x_a,y_a,z_a)$. Dans Gmsh, cela se définit de façon la plus simple : Extrude{ {vtx,vty,vtz}, {vax,vay,vaz}, {xa,ya,za}, alpha } { Surface{numero}; }
ou de façon un peu plus élaborée (si on veut structurer le maillage) : sortie[] = Extrude{ {vtx,vty,vtz}, {vax,vay,vaz}, {xa,ya,za}, alpha } { Surface{numero}; Layers{Ncouches}; Recombine; };
Busbar avec potentiel flottant
Tout comme dans le cas électrostatique de la section précédente, il est possible de définir, dans l’espace fonctionnel d’approximation $W^0$ , des grandeurs globales associées à notre inconnue $v$. Dans le cas présent, elles correspondront aux potentiels globaux $V_i$ imposés sur certaines faces du bords du domaine, ou aux courants $I_j$ traversant certaines autres.
En utilisant (ou pas) ces grandeurs globales, calculer la distribution de courant à l’intérieur du busbar en cuivre représenté ci-dessous lorsque qu’un courant de 2000 A arrive par le haut et se sépare en 2 fois 1000 A qui sortent par le bas.
Pour éviter de perdre trop de temps, je vous fournis la géométrie (libre à vous de l’utiliser ou non) : Cliquer ici pour la télécharger.
Proposer des améliorations possibles du design d’un tel busbar (en terme de quantité de matière, masse, volume, pertes…).
Magnétostatique
La magnétostatique (en anglais magnetostatics) correspond au calcul de la distribution du champ ou de l’induction magnétiques produit par des courants continus ou des aimants.
Considérons le domaine $\Omega$ de $\mathbb{R}^3$ constitué de plusieurs sous-domaines schématisé par la figure :
$\Omega_m$ correspondant à des matériaux magnétiques ;
$\Omega_c$ correspondant à des conducteurs électriques parcourus par une densité de courant imposée par une source extérieure ${\bf j_s}$ ;
$\Omega_a$ correspondant à des aimants permanents d’induction rémanente ${\bf b_r}$ ;
une zone d’air entourant ses différents sous-domaines.
Et sa frontière $\Gamma (=\partial\Omega)$ est divisée en deux morceaux $\Gamma_d$ et $\Gamma_n$.
Les relations régissant la magnétostatique sont :
Les équations de Maxwell en statique suivantes : $$\left\{\begin{aligned}\text{div}\,{\bf b} &= 0 \\ {\bf rot}\,{\bf h} &= {\bf j_s}~\text{dans } \Omega_c, \text{ et } {\bf 0} \text{ ailleurs}\end{aligned}\right.$$
Les lois de comportements magnétiques : $$\left\{\begin{aligned}{\bf b} &= \mu_0\,{\bf h}&\text{dans l'air et }\Omega_c \\ {\bf b} &= \mu\,{\bf h} &\text{dans~}\Omega_m \\ {\bf b} &= \mu_0\mu_{ra}\, {\bf h}+{\bf b_r} &\text{dans~}\Omega_a\end{aligned}\right.$$
Les contions aux limites, par exemple :
continuité de la composante normale de l’induction : ${\bf n}\cdot{\bf b} = 0$ sur $\Gamma_d$ ;
continuité de la composante tangentielle du champ : ${\bf n}\wedge{\bf h} = {\bf 0}$ sur $\Gamma_n$.
Formulation en potentiel vecteur magnétique
Formulation forte
${\bf b}$ étant à divergence nulle, nous avons vu dans le premier chapitre, qu’il était possible de définir un potentiel vecteur magnétique ${\bf a}$ tel que :
$${\bf b} = {\bf rot}\,{\bf a}$$
En lui associant les lois de comportement et en reportant dans l’équation de Maxwell-Ampère, on obtient :
Rappelons ici que la formulation ci-dessus ne conduit pas à une solution unique comme vu dans les chapitres précédents, il nous faudra rajouter une condition de jauge sur ${\bf a}$ pour pouvoir le calculer en 3D (jauge de Coulomb ou jauge d’arbre).
Formulation en potentiel scalaire magnétique
Formulation forte
Dans le cas particulier d’un problème sans densité de courant source $({\bf j_s} = {\bf 0})$, on a : ${\bf rot}\,{\bf h} = {\bf 0}$ dans tout le domaine.
On peut alors définir un potentiel scalaire magnétique $\phi$ tel que :
$${\bf h} = -{\bf grad}\,\phi$$
En lui associant les lois de comportement et en reportant dans l’équation de Maxwell-Thomson, on obtient :
Dans le cas de problèmes 2D plans ou axisymétriques, la formulation en potentiel vecteur magnétique se ramène à calculer uniquement ses composantes $a_z(x,y)$ ou $a_\theta(r,z)$.
Dans ces cas, la jauge de Coulomb $(\text{div}\,{\bf a} = 0)$ est automatiquement vérifiée et n’a pas besoin d’être spécifiquement imposée.
Pour définir $W^1$ approximant $\textbf{H}({\bf rot},\Omega)$, il existe dans GetDP des fonctions de base spécialement dédiées. Elles sont associées à des arrêtes fictives s’appuyant sur les nœuds du domaine et perpendiculaires au plan d’étude. En pratique, nous pourrons définir notre espace d’aproximation comme suit :
Trouver les conditions aux limites à imposer aux frontières du domaine (cubique) permettant d’obtenir un champ horizontal uniforme de $1~\text{T}$ à l’intérieur.
Ajouter un cylindre au centre du domaine pour traiter les 3 cas possibles en fonction du type de matériau :
Dimensionnement d’une inductance pour électronique de puissance
Pour des besoins d’électronique de puissance (ne m’en demandez pas plus), on désire concevoir une inductance de 200 µH sous 15 A avec :
des noyaux de ferrites de type E 42/21/20 dont les caractéristiques sont fournies par cette datasheet ;
du fil de cuivre émaillé de 2,5 mm² (diamètre de 1,8 mm)
Simuler le problème correspondant et donner une solution.
Information
Après avoir traité le problème avec une loi de comportement magnétique linéaire dans la ferrite ($\mu_r$ constant), il faudra sûrement résoudre le problème en prenant en compte la non-linéarité de cette loi. Pour ce faire, je vous renvoie à cette explication très bien faite sur le wiki de GetDP. Vous pouvez également consulter cette page.
Je vous fournis également la définition complète du matériau (N27) de la ferrite (en utilisant les mêmes notations que dans les exemples du wiki) :
Pour le calcul du couple, une bonne approche consiste à moyenner, sur le volume d’entrefer sous un pôle $(V_e = S_e\,L_z)$, l’expression obtenue par le tenseur de Maxwell. Montrer préalablement (sur papier) que cette moyenne peut s’écrire :
Conditions de Jauge sur le potentiel vecteur magnétique
Nous avons déjà vu plusieurs fois précédemment que le potentiel vecteur
magnétique ${\bf a}$ était défini à un champ de gradient près et qu’il était
donc nécessaire de le jauger pour assurer l’unicité de la solution.
Dans le chapitre précédent, deux méthodes ont été présentées
et nous allons maintenant voir leur application à nos problèmes ainsi
que la manière de les définir dans GetDP.
Jauge de Coulomb
Par application directe de la méthode du chapitre 2, nous allons définir un champ scalaire $\xi$ dont le gradient permettra d’assurer $\text{div}\,{\bf a} = 0$ en l’introduisant comme multiplicateur de Lagrange dans la formulation.
Ainsi la formulation ccomplète dans le domaine discret sera :
Le principal inconvénient de cette méthode est d’aboutir à une taille de système
conséquente conduisant à des temps de calcul élévés, c’est pourquoi on lui préfèrera
généralement la technique suivante.
Jauge d’arbre
Comme vu dans cette partie,
nous allons construire un arbre sur l’ensemble des arrêtes du domaine en partant
des surfaces sur lesquelles sont imposées des conditions particulières
(dans le cas présent, les conditions de Dirichlet) et imposer une circulation nulle
de ${\bf a}$ le long de cet arbre directement dans l’espace fonctionnel associé.
Ce dernier sera alors noté $W_{0,\text{JA}}^1$, et notre formulation devient ainsi :
Et la formulation sera strictement la même que dans les cas 2D précédents.
Applications
Retour sur notre inductance pour électronique de puissance
Dans le problème de dimensionnement d’inductance de la page précédente, nous avons
utilisé un modèle 2D alors que l’épaisseur de la ferrite utilisée est relativement faible.
Il serait donc judicieux de valider (ou non) cette hypothèse par un calcul 3D.
Pour ce faire, je vous fournis le modèle complet (réduit à un huitième de géométrie gràce aux symétries) à Télécharger ici.
Avec, vous pourrez :
tester les deux types de jauge en linéaire et non-linéaire ;
observer les valeurs des champs tels que l’induction dans le circuit magnétique par exemple :
tester votre design et vérifier si vous respectez le cahier des charges (200 µH sous 15 A) ;
tracer la valeur de l’inductance en fonction du courant :
Retour sur notre accouplement à aimants permanents
L’accouplement étudié page précédente reposait lui aussi sur l’hypothèse que les
effets de bords sont négligeables et il serait opportun de le vérifier.
Résoudre le problème en 3D et comparer les deux formulations :
en potentiel vecteur magnétique (avec la jauge que vous voulez) ;
en potentiel scalaire magnétique.
Pour gagner du temps, je vous fournis encore la géométie 3D en cliquant ici. Un exemple de géométrie possible est représenté ci-dessous :
Représentation schématique du coupleur considéré
Magnétodynamique
La magnétodynamique (en anglais magnetodynamics) correspond au calcul de la distribution du champ (ou de l’induction) magnétique et des courants induits produits par des courants variables dans le temps et/ou par des sources en mouvement (aimants ou courants).
Considérons un domaine $\Omega$ de $\mathbb{R}^3$ comportant un sous-domaine $\Omega_c$ conducteur électrique, dont la frontière $\Gamma (=\partial\Omega)$ est comme précédemment divisée en deux morceaux $\Gamma_d$ et $\Gamma_n$.
Les lois de comportements : $$\left\{\begin{aligned}{\bf b} &= \mu\,{\bf h}&\text{dans }\Omega \\ {\bf j} &= \sigma\,{\bf e} &\text{dans~}\Omega_c\end{aligned}\right.$$
Les contions aux limites, par exemple :
continuité de la composante normale de l’induction : ${\bf n}\cdot{\bf b} = 0$ sur $\Gamma_d$ ;
continuité de la composante tangentielle du champ : ${\bf n}\wedge{\bf h} = {\bf 0}$ sur $\Gamma_n$.
Formulation forte
Le potentiel vecteur magnétique ${\bf a}$ est défini comme précédemment et donc ${\bf b} = {\bf rot\,a}$. En le repportant dans l’équation de Maxwell-Faraday, on obtient :
On peut donc définir dans $\Omega_c$ un champ scalaire $v$, potentiel scalaire électrique tel que ${{\bf e}-\frac{\partial\,{\bf a}}{\partial t} = - {\bf grad}\,v}$, soit :
Les conditions aux limites sont analogues à celles vues en magnétostatique ou électrocinétique.
Formulation faible
En combinant les différentes approches vues dans les sections précédentes, on en déduit la formulation faible complète du problème, dite $({\bf a},v)$ :
Dans le cas particulier où les sources sont en régime sinusoïdal forcé de pulsation $\omega$, nous pourrons utiliser la transformation complexe vue dans le chapitre 1 et résoudre directement en complexe : c’est ce qu’on appelle la magnétoharmonique (certains l’appellent « magnétostatique complexe »).
L’implantation dans GetDP n’est pas plus compliquée que ce que nous avons vu jusqu’à présent. Le passage en complexe et la fréquence associée sont précisés dans la partie « Résolution » :