Version 2.1 / Janvier 2026
Sur le webpoly de l’EC B8-11-EC1 : « Électromagnétisme BF analytique et numérique » du bloc de compétence « Modélisation & Simulation numérique en Génie Électrique » de 2ème année ENSEM (formation Énergie).
C’est un petit site internet à mi-chemin entre polycopié de cours et diaporama. L’idée est d’unifier tout le contenu pédagogique de l’EC avec quelque chose de plus synthétique qu’un polycopié de cours classique, mais plus détaillé que de simples transparents. Le but est donc d’avoir un support multi-usages pouvant s’adapter aux différents modes d’enseignement imposés selon les circonstances : présentiel, hybride ou distanciel (synchrone ou asynchrone).
Se (re)mettre à niveau en électromagnétisme en présence de milieux matériels (base fondamentale de l’électrotechnique / électromécanique).
Se perfectionner sur la méthode des éléments finis appliquée aux problèmes correspondants.
Comprendre et appréhender ce qui est caché derrière les interfaces des logiciels du domaine.
Connaître les formulations fortes de l’électromagnétisme dans l’approximation des régimes quasi stationnaires.
Savoir en déduire les formes faibles associées.
Maîtriser les particularités du domaine :
conditions aux limites correspondant au traitement de frontières lointaines,
types d’éléments,
conditions de jauge particulières,
formulations en potentiels.
Savoir interpréter, adapter et modifier un modèle développé avec l’interface logicielle ONELAB (mailleur libre Gmsh et solveur libre GetDP).
Où on reviendra sur les notions de base.
Rapide historique de la naissance de l'Électromagnétisme.
Quelles sont les grandeurs physiques correspondant au champ électromagnétique ?
Relations constitutives dans les milieux matériels considérés.
What else?
Formes particulières en basse fréquence.
Définition des potentiels, bien utiles pour calculer le champ.
Retrouvons les lois et théorèmes connus à partir des équations locales…
Détails des relations de passage entre deux milieux matériels distincts.
Comment calcule-t'on les grandeurs globales d'intérêt (Énergies, flux, forces, etc…) ?
Un peu de complexes dans ce monde réel…
Mettons en pratique toute cette partie théorique.
Bien que plus jeune que la mécanique, l’électromagnétisme est, avec cette dernière, une grande branche de la physique classique. Sur Wikipédia, on peut trouver un bon résumé historique. Plutôt que de le reformuler, je préfère le pomper odieusement :
Pendant longtemps les phénomènes électriques et les phénomènes magnétiques ont été considérés comme indépendants [1]. En 1600, William Gilbert explicite, dans son ouvrage De Magnete, la distinction entre corps électriques (il introduit ce terme) et magnétiques. Il assimile la Terre à un aimant, note les répulsions et attractions des aimants par leurs pôles et l’influence de la chaleur sur le magnétisme du fer. Il donne aussi les premières notions sur l’électricité, dont une liste des corps électrisables par frottement.
Les Grecs avaient seulement remarqué que des morceaux d’ambre frottés pouvaient attirer des corps légers, tels des copeaux ou de la poussière, et par ailleurs qu’il existait un minéral, la « pierre d’aimant » ou magnétite [2], capable d’attirer le fer et les métaux ferreux.
La découverte au XIXe siècle par Ørsted, Ampère et Faraday de l’existence d’effets magnétiques de l’électricité a conduit progressivement à envisager que les forces « électrique » et « magnétique » puissent être en fait unifiées, et Maxwell propose en 1860 une théorie générale de l’électromagnétisme classique, qui pose les fondements de la théorie moderne.
La perturbation des boussoles sous l’action de la décharge de la foudre était un phénomène bien connu au XVIIIe siècle. Cela créait un lien entre électricité et magnétisme, mais difficile à interpréter et impossible à reproduire. Par ailleurs les lois de l’électricité et du magnétisme énoncées par Charles Coulomb distinguaient bien l’électricité d’un côté et le magnétisme de l’autre, même si ces lois se présentaient sous la même forme mathématique.
En 1820, le Danois Hans Christian Ørsted fait une observation extraordinaire : un fil rectiligne parcouru par un courant continu dévie l’aiguille d’une boussole placée à proximité.
En 1820, André-Marie Ampère met en évidence les interactions entre courants électriques et assimile tout aimant, y compris le globe terrestre, à un ensemble de courants [3].
En 1831, Michael Faraday étudie le comportement d’un courant dans un champ magnétique, et s’aperçoit que celui-ci peut produire du travail. Ørsted avait découvert qu’un courant électrique produit un champ magnétique, Faraday découvre qu’un champ magnétique engendre un courant électrique. Il découvre ainsi le principe du moteur électrique, et donc la conversion du travail mécanique en énergie électrique, inventant ainsi la génératrice de courant. Dans un article de 1852 (« On the Physical Character of the Lines of Magnetic Force »), Faraday dévoile l’existence du champ magnétique en décrivant les « lignes de force » le long desquelles s’oriente la limaille de fer au voisinage de l’aimant.
En 1864, James Maxwell unifie les théories antérieures, comme l’électrostatique, l’électrocinétique ou la magnétostatique. Cette théorie unifiée explique, entre autres, le comportement des charges et courants électriques, des aimants, ou encore des ondes électromagnétiques telles la lumière ou les ondes radio qui apparaissent en fait comme la propagation de perturbations électromagnétiques. L’électromagnétisme est né.
⇝ Merci Wikipédia ! Mais pour le reste, on fera sans toi… 😉
De façon générale, à l’échelle macroscopique, le champ électromagnétique est constitué de cinq champs vectoriels et d’un champ scalaire :
${\bf h}$ : le champ magnétique en $\text{A}\cdot\text{m}^{-1}$
${\bf e}$ : le champ électrique en $\text{V}\cdot\text{m}^{-1}$
${\bf b}$ : l’induction magnétique en $\text{T}$ (= $\text{Wb}\cdot\text{m}^{-2}$)
${\bf d}$ : l’induction électrique en $\text{C}\cdot\text{m}^{-2}$
${\bf j}$ : la densité volumique de courant électrique en $\text{A}\cdot\text{m}^{-2}$
$\rho_q$ : la densité volumique de charge électrique en $\text{C}\cdot\text{m}^{-3}$
Dans ce cours, nous utiliserons les notations anglo-saxonnes pour les vecteurs : il seront donc affichés en gras.
Les six grandeurs ci-dessus sont liées entres-elles par deux types de relations :
Nous distinguerons toujours les inductions et les champs car nous serons en présence d’un ou plusieurs milieux matériels. En fait, nous nous plaçons dans le cadre de ce qu’on appelle l’Électrodynamique des milieux continus [Land84].
Dans un milieu matériel donné, les différents champs constituant le champ électromagnétique sont liés entre eux par les caractéristiques physiques du matériau. Ces relations sont modélisées par ce qu’on appelle les lois de comportement.
En ce qui nous concerne, on en distinguera trois types :
Dans un milieu non magnétique tel que le vide (ou l’air), mais aussi dans les isolants ou conducteurs électriques classiques (cuivre, aluminium, etc…), l’induction ${\bf b}$ est directement proportionnelle au champ magnétique ${\bf h}$ par :
$$\boxed{{\bf b} = \mu_0\,{\bf h}}$$Où $\mu_0$ est la perméabilité magnétique du vide :
$$ \mu_0 = 4 \pi~10^{-7}~~\text{H}\cdot\text{m}^{-1} $$En réalité, hormis le vide, un matériau est souvent légèrement diamagnétique (tendance à repousser le champ d’induction) ou paramagnétique (tendance à attirer le champ d’induction). On a une relation de type ${\bf b} = \mu_0\,\mu_r~{\bf h}$ avec $\mu_r$ très légèrement inférieur ou supérieur à 1. En pratique, ces effets sont suffisamment faibles pour être négligés.
Ce sont des matériaux dans lesquels apparaît une aimantation ${\bf m}$ (en $\text{A}\cdot\text{m}^{-1}$) lorsqu’ils sont soumis à un champ magnétique ${\bf h}$.
La loi correspondante est :
L’évolution de la valeur de l’induction, dans la direction privilégiée, en fonction de celle du champ décrit alors un cycle d’hystérésis comme sur la figure ci-dessous :
Exemple de cycle d'hystérésis
L’aimantation ${\bf m}$ s’exprime elle-même en fonction de ${\bf h}$ par :
$$ {\bf m} = {\bf m_0} + \chi_m\,{\bf h}$$avec :
On pourrait retrouver les expressions reliant ${\bf m_0}$ et $\chi_m$ à $B_r$ ou $H_c$.
Les matériaux doux sont caractérisés par un cycle très étroit, comme par exemple :
Exemple de cycle d'hystérésis pour un matériau doux
Les valeurs de champ coercitif ou d’induction rémanente sont suffisamment faibles pour que le cycle puisse être assimilé à une courbe passant par l’origine, dont l’allure peut se représenter comme ci-dessous :
Exemple de courbe B(H)
L’aimantation permanente du matériau est négligeable, et on en déduit :
$$\boxed{{\bf b} = \mu_0\,\underbrace{(1+\chi_m)}_{\displaystyle\mu_r}\,{\bf h} = \mu\,{\bf h}}$$$\mu$ est la perméabilité magnétique du matériau (en $\text{H}\cdot\text{m}^{-1}$), elle dépend généralement du champ : $\mu({\bf h})$.
Dans la suite, nous utiliserons deux types de modèles :
C’est ce type de matériaux qui est utilisé pour les circuits magnétiques des applications qui nous intéresseront : inductances, transformateurs, actionneurs, machines tournantes…
À l’inverse des précédents, les ferromagnétiques durs sont caractérisés par un cycle d’hystérésis très large, comme :
Exemple de cycle d'hystérésis pour un matériau dur
Dans la suite, on se situera toujours sur la partie linéaire correspondant à la partie haute du cycle. Par substitution, on a :
$$ {\bf b} = \mu_0\,({\bf h}+\chi_m\,{\bf h}+{\bf m_0})$$Qui s’écrit :
$$\boxed{{\bf b} = \mu_0\underbrace{(1+\chi_m)}_{\displaystyle\mu_{ra}}\,{\bf h}+\underbrace{\mu_0\,{\bf m_0}}_{\displaystyle{\bf b_r}}}$$Avec :
Elles sont valables pour des matériaux conducteurs électriques. On peut distinguer deux cas.
Cette catégorie correspondant à l’ensemble des conducteurs que nous verrons par la suite : principalement des métaux comme le cuivre, l’aluminium, le laiton, les aciers, etc…
La loi de comportement, y est simplement :
$$\boxed{{\bf e} = \rho\,{\bf j}}$$
où $\rho$ est la résistivité électrique du matériau (en $\Omega\cdot\text{m}$).
On pourra éventuellement lui préférer la relation inverse (strictement équivalente) :
avec $\sigma = \frac{1}{\rho}$, la conductivité électrique du milieu (en $\text{S}\cdot\text{m}^{-1}$).
Souvent, les deux relations ci-dessus ne sont pas linéaires car la résistivité et (donc) la conductivité dépendent fortement de la température. Par exemple, pour la conductivité du cuivre :
Évolution de la conductivité du cuivre en fonction de la température
Ce sont des matériaux qui, dans les bonnes conditions (de température et de champs) présentent une résistivité électrique nulle ou quasi-nulle. Il sont caractérisés par une loi de comportement du type :
$$ {\bf e} = E_c \left(\frac{ \lVert{\bf j} \lVert}{J_c}\right)^n \frac{{\bf j}}{ \lVert {\bf j} \rVert}$$
où $E_c$ et $J_c$ sont respectivement le champ et la densité de courant critiques du matériau.
En réalité, la loi ci-dessus fait intervenir d’autre phénomènes couplés car $J_c$ et $n$ dépendent des valeurs de l’induction et de la température.
Information
Mais, dans le cadre de ce cours, nous n’aborderons malheureusement pas l’étude des supraconducteurs faute de temps…
Un milieu diélectrique est un milieu isolant (non-conducteur du courant) qui peut se polariser électriquement lorsqu’il est soumis à un champ électrique (il est constitué de dipôles électriques). C’est le cas par exemple du verre, de céramiques (dont les céramiques piézoélectriques), du polypropylène, du mica, du téflon, etc… D’un point de vue « Loi de comportement », on peut faire une analogie avec les milieux magnétiques précédents.
Dans le cas du vide, ou de nombreux autres milieux non-diélectriques comme les conducteurs ou matériaux magnétiques. La relation entre l’induction et le champ électrique est simplement :
$$\boxed{{\bf d} = \varepsilon_0\,{\bf e}}$$où $\epsilon_0$ est la permittivité électrique du vide (en $\text{F}\cdot\text{m}^{-1}$), dont la valeur se déduit de la relation :
$$ \varepsilon_0\,\mu_0\,c^2 = 1$$En prenant $3\cdot 10^8~\text{m}\cdot\text{s}^{-1}$ pour la vitesse de la lumière, on approxime souvent $\varepsilon_0$ par :
$$ \varepsilon_0 = \frac{1}{36\,\pi} 10^{-9} ~~\text{F}\cdot\text{m}^{-1}$$Comme pour l’aimantation des matériaux magnétiques, il apparaît une polarisation électrique ${\bf p}$ dans les diélectriques lorsqu’ils sont soumis à un champ électrique ${\bf e}$.
La loi donnant l’induction ${\bf d}$ est alors :
La polarisation peut s’écrire sous la forme :
$$ {\bf p} = \varepsilon_0\,\chi_e\,{\bf e} + {\bf p_0}$$avec :
Ainsi, dans une direction privilégiée, la valeur de l’induction en fonction de celle du champ électrique décrit elle-aussi un cycle comme le montre la figure ci-dessous :
Exemple de cycle d'hystérésis d'un diélectrique
Dans le cadre de ce cours, nous ne intéresserons qu’aux cas linéaires où ${\bf p_0} = {\bf 0}$.
On aura donc :
où $\varepsilon_r$ est la permittivité relative du matériau considéré.
Nous avons volontairement éviter de trop détailler cette (déjà trop longue) partie car ce n’est pas l’objectif de l’EC. Nous aurions pu, par exemple, aborder des notions telles que l’anisotropie de certains matériaux ou des couplages multiphysiques supplémentaires (effet Hall, piézoélectricité, etc…), ou alors discuter des semi-conducteurs.
Pour des compléments, je vous renvoie à l’EC de S7 Matériaux pour l’Energie, partie sur les matériaux pour le Génie Électrique.
Les quatre équations de Maxwell sont les relations fondamentales régissant les lois de l’Électromagnétisme.
Leur forme universelle à l’échelle macroscopique est donnée ci-dessous :
S’il n’y avait qu’une seule chose à retenir en électromagnétisme et à connaître par cœur, c’est celle-là.
Remarque : En prenant la divergence de chaque terme de l’équation de Maxwell-Ampère, on obtient également l’équation de conservation de la charge :
$$ \text{div}~{\bf j} + \frac{\partial\,\rho_q}{\partial t} = 0$$Petit exercice pour faire une pause :
À partir de la forme générale des équations de Maxwell ci-dessus et des lois de comportement précédentes, retrouver comment sont modifiées les équations de Maxwell-Ampère et de Maxwell-Gauss que vous aviez l’habitude d’utiliser en premier cycle.
Rappelons que le but du présent module est de vous permettre de modéliser et simuler par éléments finis des dispositifs typiques de l’électrotechnique / électromécanique correspondant donc à de l’électromagnétisme basse-fréquence. La plage de fréquence concernée va donc du continu $(0~\text{Hz})$ à des valeurs de l’ordre de quelques centaines de $\text{kHz}$.
Considérons que le dispositif étudié est alimenté par une source sinusoïdale (de tension ou de courant) de pulsation $\omega$ et que les lois de comportement électrique et diélectrique des milieux matériels présents sont linéaires. En notant $\tau_c$ le rapport des normes des termes $\frac{\partial\,{\bf d}}{\partial t}$ et ${\bf j}$, on obtient après simplifications :
$$ \tau_c = \frac{\varepsilon_0\,\varepsilon_r\,\omega}{\sigma}$$La table ci-dessous donne quelques valeurs caractéristiques :
| Matériau | $\sigma~(\text{S}\cdot\text{m}^{-1})$ | $\varepsilon_r$ | $\tau_c$ @ $1\,\text{MHz}$ |
|---|---|---|---|
| Cuivre | $58~10^{6}$ | $1$ | $9,93~10^{-13}$ |
| Aluminium | $37,7~10^{6}$ | $1$ | $1,48~10^{-12}$ |
| Fer pur | $10,4~10^{6}$ | $1$ | $5,35~10^{-12}$ |
| Eau salée | $5$ | $80$ | $8,90~10^{-4}$ |
On constate que, dans tous ces cas, le terme $\displaystyle\frac{\partial\,{\bf d}}{\partial t}$ est clairement négligeable devant la densité de courant ${\bf j}$.
Dans la suite, nous le négligerons donc toujours, et nous nous placerons donc systématiquement dans ce qu’on appelle l’Approximation des Régimes Quasi-Statiques (ARQS).
Ainsi, la forme des équations de Maxwell que nous résoudrons dans ce cours est celle dans le cadre de l’ARQS, donnée ci-dessous :
Équations de Maxwell dans l’ARQS :
$$\left\{\begin{aligned} \text{div}\,{\bf b} &= 0 & \text{(MT)} \\ \text{div}\,{\bf d} &= \rho_q & \text{(MG)} \\ {\bf rot\,h} &= {\bf j} & \text{(MA)} \\ {\bf rot\,e} &= -\frac{\partial\,{\bf b}}{\partial t}~~~ & \text{(MF)} \end{aligned}\right.$$On remarque que désormais l’équation de Maxwell-Gauss est totalement découplée des autres.
La divergence appliquée à (MA), nous donne la loi des nœuds locale :
$$\boxed{\text{div}\,{\bf j} = 0}$$Plusieurs cas particuliers peuvent se présenter en fonction de la nature ou du type de source. Chacun définit un sous-domaine de l’électromagnétisme (en basse fréquence).
L’Électrostatique porte sur la distribution de champ électrique due à des charges statiques ou des différences de potentiel électrique.
Les équations qui y interviennent sont donc :
$$\left\{\begin{aligned} &\text{div}\,{\bf d} = \rho_q \\ &{\bf rot}\,{\bf e} = {\bf 0} \\ &{\bf d} = \varepsilon_0 \varepsilon_r\,{\bf e} \end{aligned}\right.$$Nous résoudrons ici certains problèmes correspondants.
L’Électrocinétique concerne la distribution de la densité de courant à l’intérieur des conducteurs (lorsque les effets magnétiques sont négligeables : en statique, ou à fréquence suffisamment faible).
Les équations qui y interviennent sont donc :
$$\left\{\begin{aligned} &{\bf rot}\,{\bf e} = {\bf 0} \\ &\text{div}\,{\bf j} = 0 \\ &{\bf j} = \sigma\,{\bf e} \end{aligned}\right.$$Nous résoudrons ici certains problèmes correspondants.
La Magnétostatique, quand à elle, traite de la distribution de champ magnétique créé par une densité de courant constante et/ou des aimants permanents.
Les équations qui y interviennent sont donc :
$$\left\{\begin{aligned} &{\bf rot}\,{\bf h} = {\bf j} \\ &\text{div}\,{\bf b} = 0 \\ &{\bf b} = \mu\,{\bf h} ~(+ {\bf b_r}) \end{aligned}\right.$$La densité de courant ${\bf j}$ sera soit imposée comme densité de courant source circulant dans des inducteurs (${\bf j_s}$), soit issue d’une résolution électrocinétique (cf ci-dessus).
Nous résoudrons ici certains problèmes correspondants.
Enfin, la Magnétodynamique s’occupe des distributions de champ magnétique et de courant induit dans les conducteurs issues de sources de courant variant dans le temps et/ou de sources statiques (courant continu ou aimants permanents) en mouvement.
Les équations qui y interviennent sont donc :
$$\left\{\begin{aligned} &{\bf rot}\,{\bf h} = {\bf j} \\ &\text{div}\,{\bf b} = 0 \\ &{\bf rot\,e} = -\frac{\partial\,{\bf b}}{\partial t} \\[0.5em] &{\bf b} = \mu\,{\bf h} ~(+ {\bf b_r}) \\ &{\bf j} = \sigma\,{\bf e} \end{aligned}\right.$$Nous verrons ici comment résoudre certains problèmes correspondants.
Pour traiter les problèmes avec des conducteurs en mouvement, il est parfois plus utile d’utiliser la loi d’Ohm locale généralisée :
$${\bf j} = \sigma\,({\bf e} + {\bf v}\wedge{\bf b})$$Celle-ci ne contredit pas la loi de comportement électrique car elle est exprimée dans le référentiel de l’observateur, il y a donc toujours $~{{\bf j} = \sigma\,{\bf e}}~$ dans le référentiel lié au conducteur (${{\bf v}={\bf 0}}$). L’intérêt est de pouvoir transformer, dans certains cas, le problème initial en un problème de magnétostatique plus simple à résoudre. Si cela vous intéresse, vous pouvez consulter l’exemple détaillé dans cet article.
Pour illustrer cette remarque, et faire une petite pause, je vous propose une vidéo expliquant le principe de fonctionnement d’un ralentisseur Telma, dispositif dont le fonctionnement est un parfait exemple de problème magnétodynamique qu’on peut résoudre avec une formulation statique avec la loi d’Ohm locale généralisée :
Il existe un cas particulier que nous n’aborderons pas dans le cadre de ce cours, mais dans lequel on ne peut pas négliger le terme $\frac{\partial\,{\bf d}}{\partial t}$ bien que correspondant à des applications en basse fréquence et donc entrant dans le cas quasi-statique.
En effet certains matériaux peuvent être à la fois conducteurs (avec des conductivités assez faibles mais non nulles) et posséder une perméabilité relative assez élevée. Les semi-conducteurs rentrent dans ce cadre par exemple, mais aussi beaucoup de tissus biologiques comme ceux des organes humains. On peut aussi s’intéresser à des problèmes comportant, à la fois, des milieux conducteurs et des isolants dans lesquels on cherche à calculer les distributions de champ électrique et de densité de courant en s’affranchissant complètement des aspects magnétique. On appelle parfois ce domaine d’étude Électrodynamique ou Modèle électrique quasi-statique.
Les équations qui y interviennent sont donc :
$$\left\{\begin{aligned} &{\bf rot}\,{\bf e} = {\bf 0} \\ &\text{div}\,\left({\bf j}+\frac{\partial\,{\bf d}}{\partial t}\right) = 0 \\ & {\bf j} = \sigma\,{\bf e} \\&{\bf d} = \varepsilon\,{\bf e} \end{aligned}\right.$$Mais, comme dit ci-dessus, les applications concernées n’entrent pas dans le cadre de ce module.
On se place dans un domaine $\Omega$ (ouvert étoilé quelconque de $\mathbb{R}^3$), et on notera $\Omega_c \subset \Omega$ le sous-domaine conducteur en présence de courants. Nous allons voir que les champs définis précédemment peuvent découler de différents potentiels.
Prenons comme point de départ l’équation de Maxwell-Thomson dans ce domaine et regardons ce qu’elle implique :
${\bf a}$ est appelé potentiel vecteur magnétique.
Attention le potentiel vecteur ${\bf a}$ ainsi défini n’est pas unique ! Il est défini à un champ de gradient près.
Le potentiel vecteur ${\bf a'}$ convient donc lui-aussi.
Pour le rendre unique, il est donc nécessaire de rajouter une condition de jauge. Dans l’ARQS, la plus connue est la jauge de Coulomb :
Important !
En effet, soit ${\bf a}$ tel que ${\bf b} = {\bf rot\,a}$, alors :
Mais il en existe d’autres parfois plus pratiques ou efficaces (d’un point de vue numérique), nous y reviendrons dans la troisième partie du module.
Même avec la condition de jauge, on remarque que le potentiel est encore défini à une constante près. Ceci ne sera pas un souci dès lors qu’on utilisera des conditions aux limites appropriées.
Astuce
Dans le cas de problèmes magnétiques, nous avons besoin de connaître les sources, donc ${\bf j}$, donc ${\bf e}$, dans les domaines conducteurs ($\Omega_c$).
Substituons alors l’expression précédente dans l’équation de Maxwell-Faraday :
$${\bf rot\,e} = -\frac{\partial}{\partial t}\,({\bf rot\,a}) \Leftrightarrow {\bf rot}\,\left({\bf e} + \frac{\partial\,{\bf a}}{\partial t}\right) = 0$$$$\Leftrightarrow \exists~v \in \Omega_c : {\bf e} + \frac{\partial\,{\bf a}}{\partial t} = -{\bf grad}\,v$$Soit :
$$\exists~v \in \Omega_c : {\bf e} = -{\bf grad}\,v - \frac{\partial\,{\bf a}}{\partial t}$$
$v$ est appelé potentiel scalaire électrique.
On remarquera que l’expression ci-dessus reste valable an remplaçant ${\bf a}$ par ${\bf a'} = {\bf a} + {\bf grad}\,u$, et $v$ par $\displaystyle v' = v + \frac{\partial\,u}{\partial t}$.
Pour des problèmes électrostatiques ou électrocinétiques, nous voulons déterminer ${\bf e}$ (et ${\bf d}$ ou ${\bf j}$) dans tout le domaine d’étude.
Alors, partant de l’équation de Maxwell-Faraday en statique, nous avons directement :
$${\bf rot\,e} = {\bf 0} \Leftrightarrow \exists~v \in \Omega : {\bf e} = -{\bf grad}\,v$$avec $v$, le potentiel scalaire électrique.
Nous avons déjà vu que l’équation Maxwell-Ampère impliquait la loi de nœuds locale. En partant de cette dernière, nous avons :
$$\text{div}\,{\bf j} = 0 \Leftrightarrow \exists~{\bf t} \in \Omega_c : {\bf j} = {\bf rot\,t}$$
${\bf t}$ est appelé potentiel vecteur électrique.
On veillera à ne pas confondre le potentiel ${\bf t}$ et la variable temps $t$…
Là encore, le potentiel vecteur ainsi défini n’est pas unique, et il sera nécessaire de rajouter des conditions pour pouvoir le déterminer.
Dans le troisième chapitre, Nous verrons aussi que, parfois, nous préférerons travailler avec un champ source (noté ${\bf h_s}$) défini lui-aussi par ${\bf j} = {\bf rot\,h_s}$ soit dans $\Omega$ entier, soit dans un sous-domaine contenant $\Omega_c$, soit dans $\Omega_c$ auquel on ajoutera d’éventuelles coupures.
En utilisant le potentiel ${\bf t}$, on peut développer Maxwell-Ampère comme :
$$ {\bf rot\,h} = \left\{\begin{array}{c l} {\bf rot\,t}~ &\text{dans}~ \Omega_{c} \\\\ {\bf 0}~ &\text{dans}~ \Omega\setminus\Omega_{c} \end{array}\right.$$Ainsi :
$$\left\{\begin{array}{r l} {\bf rot}\,({\bf h - t}) = {\bf 0}~ &\text{dans}~ \Omega_c\\\ {\bf rot\,h} = {\bf 0}~ &\text{dans}~ \Omega\setminus\Omega_c\end{array}\right.$$$$\Leftrightarrow\exists~\phi\in\Omega : \left\{\begin{array}{l l}{\bf h}={\bf t}-{\bf grad}\,\phi & \text{dans}~ \Omega_c\\\ {\bf h}= -{\bf grad}\,\phi & \text{dans}~ \Omega\setminus \Omega_c\end{array}\right.$$$\phi$ est appelé potentiel scalaire magnétique.
Dans le cas d’un problème sans courant, par exemple lorsque les seules sources de champ sont des aimants permanents, on a directement :
$$\exists~\phi\in\Omega : {\bf h} = -{\bf grad}\,\phi$$On considère un milieu matériel quelconque de volume $V$ de bord $\partial V = S$ représenté par la figure ci-dessous.
Représentation schématique du domaine considéré.
En appliquant le théorème de la divergence (théorème de Green-Ostrogradski) à l’induction magnétique ${\bf b}$ sur le volume $V$, on obtient :
$$\iiint_{V} \text{div}\,{\bf b}~\text{d} V = \oiint_{S} {\bf b}\cdot{\bf d S}$$D’après l’équation de Maxwell-Thomson, on a donc :
$$\boxed{\oiint_S {\bf b}\cdot{\bf d S} = 0}$$Ou encore :
$$\underbrace{\iint_{S_1} {\bf b}\cdot{\bf d S_1}}_{\displaystyle\Phi_1} = - \underbrace{\iint_{S_2} {\bf b}\cdot{\bf d S_2}}_{\displaystyle\Phi_2}$$Où $\Phi_1$ et $\Phi_2$ sont respectivement les flux de ${\bf b}$ à travers $S_1$ et $S_2$.
$\implies {\bf b}$ est à flux conservatif.
En appliquant le même théorème à l’induction électrique ${\bf d}$ sur le même volume $V$, et en y injectant l’équation de Maxwell-Gauss :
$$\oiint_S {\bf d}\cdot{\bf dS} = \iiint_V \text{div}\, {\bf d}~\text{d} V = \iiint_V \rho_q ~ \text{d} V $$Soit :
$$\boxed{\oiint_S {\bf d}\cdot{\bf dS} = Q_V}$$Où $Q_V$ est la charge électrique totale contenue dans le volume $V$.
Appliquons maintenant le théorème de Stokes au champ magnétique ${\bf h}$ sur le contour $C$ et reportons y l’équation de Maxwell-Ampère :
$$\oint_C {\bf h}\cdot{\bf d l} = \iint_{S_1} {\bf rot}\,{\bf h} \cdot {\bf d S_1} = \iint_{S_1} {\bf j} \cdot {\bf d S_1}$$Soit :
$$\boxed{\oint_C {\bf h}\cdot{\bf d l} = I_{S_1}}$$Où $I_{S_1}$ représente le courant total traversant la surface $S_1$ (courant enlacé par le contour $C$).
Appliquons le même théorème au champ électrique ${\bf e}$, toujours sur le contour $C$ :
$$\oint_C {\bf e}\cdot{\bf d l} = \iint_{S_1} {\bf rot}\,{\bf e} \cdot {\bf d S_1}$$En utilisant l’équation de Maxwell-Faraday, on obtient :
$$\oint_C {\bf e}\cdot{\bf d l} = - \iint_{S_1} \frac{\partial\,{\bf b}}{\partial t}\cdot{\bf d S_1}$$En considérant que $C$ reste immobile dans le référentiel d’étude, on en déduit :
$$\oint_C {\bf e}\cdot{\bf d l} = - \frac{\text{d}}{\text{d} t} \iint_{S_1} {\bf b}\cdot{\bf d S_1} $$Soit :
$$\boxed{\oint_C {\bf e}\cdot{\bf d l} = - \frac{\text{d}\,\Phi_1}{\text{d} t} }$$Où : $\oint_C {\bf e}\cdot{\bf d l}$ est la force électromotrice (f.é.m $\mathcal{E}$).
Une bonne généralisation aux cas de circuits en mouvement se trouve dans [Four85].
On a vu que la densité de courant ${\bf j}$ était elle-aussi à divergence nulle (par l’équation de Maxwell-Ampère). Par analogie avec l’induction magnétique, on obtient ainsi :
$$\boxed{\oiint_S {\bf j}\cdot{\bf d S} = 0}$$Soit :
$$\underbrace{\iint_{S_1} {\bf j}\cdot{\bf d S_1}}_{\displaystyle I_1} + \underbrace{\iint_{S_2} {\bf j}\cdot{\bf d S_2}}_{\displaystyle I_2} = 0$$$\Longrightarrow$ C’est la loi des nœuds !
Dans le cas particulier de milieux non-magnétiques (donc où ${\bf b} = \mu_0\,{\bf h}$), nous disposons d’une formule permettant de calculer une valeur locale de l’induction (au point ${\bf x_0}$) à partir de la distribution de courant source (supposée contenue dans le volume $V_s$) :
$$\boxed{{\bf b}({\bf x_0}) = \frac{\mu_0}{4\,\pi} \iiint_{V_s} \frac{{\bf j}({\bf x})\wedge({\bf x_0}-{\bf x})}{ \lVert {\bf x_0}-{\bf x} \rVert^3}~\text{d} V_{(\bf x)}}$$En réalité, la formule ci-dessus est une généralisation de la Loi de Biot et Savart initialement énoncée pour un circuit filiforme.
On dispose de formules analogues, pour :
En pratique, le principal intérêt de ce type de formule est de, parfois, permettre de calculer le champ source ou le potentiel vecteur électrique dans une formulation en potentiel scalaire magnétique.
Considérons deux milieux $\#1$ et $\#2$ séparés par une surface $S$, la normale à la surface étant dirigée de $\#1$ vers $\#2$. Pour un infiniment petit quelconque $u$, on définit un petit cylindre de volume $v_c$, de rayon $r = o(u)$ et de hauteur $l = o(u^2)$ centré sur $S$.
Une représentation schématique est donnée par la figure ci-dessous.
En notant $s_1$ et $s_2$ les bases du cylindre contenues respectivement dans $\#1$ et $\#2$, et $s_l$ sa surface latérale, on a :
$$\left\{\begin{aligned}s_1 = s_2 = s = o(u^2)\\ s_l = o(u^3)\end{aligned}\right.$$Appliquons la conservation du flux à notre cylindre :
$$\oiint_{\partial v_c} {\bf b}\cdot{\bf d S} = 0$$Soit, en développant et en multipliant par $\frac{1}{s}$ :
$$\frac{1}{s} \iint_{s_1} {\bf b}\cdot(-{\bf n_{12}})~ \text{d} S + \frac{1}{s} \iint_{s_2} {\bf b}\cdot{\bf n_{12}}~\text{d} S + \frac{1}{s} \iint_{s_l} {\bf b}\cdot{\bf d S} = 0$$Or :
$$\frac{1}{s} \iint_{s_l} {\bf b}\cdot{\bf d S} \leq \max_{s_l}( \lVert {\bf b} \rVert)\,\frac{s_l}{s} = o(u)$$Et :
$$\left\{\begin{aligned}\lim_{u \to 0}\left(\frac{1}{s} \iint_{s_1} {\bf b}\cdot(-{\bf n_{12}})~ \text{d} S\right) &= - {\bf b_1}\cdot{\bf n_{12}} \\ \lim_{u \to 0}\left(\frac{1}{s} \iint_{s_2} {\bf b}\cdot{\bf n_{12}}~\text{d} S\right) &= {\bf b_2}\cdot{\bf n_{12}}\end{aligned}\right.$$Où ${\bf b_1}$ et ${\bf b_2}$ sont les valeurs de l’induction au centre du cylindre approché respectivement depuis les milieux $\#1$ ou $\#2$.
Ainsi par passage à la limite dans l’expression générale, on obtient :
$$\boxed{{\bf n_{12}}\cdot({\bf b_2}-{\bf b_1}) = 0}$$On a donc continuité de la composante normale de l’induction magnétique ${\bf b}$ à travers la surface.
Appliquons le raisonnement précédent en partant du théorème de Gauss :
$$\oiint_{\partial v_c} {\bf d}\cdot{\bf dS} = \iiint_{v_c} \rho_q\,\text{d} V$$Le terme de gauche peut-être traité exactement de la même façon que pour ${\bf b}$. Pour le terme de droite, on peut écrire :
$$\frac{1}{s} \iiint_{v_c} \rho_q\,\text{d} V = \frac{1}{s} \iint_{s} \left(\int_{-\frac{l}{2}}^{\frac{l}{2}} \rho_q\,\text{d}{l}\right) \text{d} S$$En définissant localement la densité surfacique de charge $\sigma_q$ (densité de charge portée par la surface $s$, en $\text{C}\cdot\text{m}^{-2}$) par :
$$\sigma_q = \lim_{l \to 0} \int_{-\frac{l}{2}}^{\frac{l}{2}} \rho_q\,\text{d}{l}~,$$par passage à la limite ($u \to 0$), on obtient :
$$\boxed{{\bf n_{12}}\cdot({\bf d_2}-{\bf d_1}) = \sigma_q}$$On a donc un saut de la composante normale de ${\bf d}$ à travers $S$.
On trouve des densités surfaciques de charge dans tout conducteur à l’équilibre électrostatique.
En l’absence de densité surfacique de charge : on a évidemment continuité de la composante normale de l’induction électrique ${\bf d}$.
Dans l’ARQS, puisque $\text{div}\,{\bf j} = 0$, on peut appliquer le même raisonnement que pour ${\bf b}$, et ainsi :
$$\boxed{{\bf n_{12}}\cdot({\bf j_2}-{\bf j_1}) = 0}$$En particulier, sur le bord d’un conducteur isolé, ${\bf j_2}$ étant nulle à l’extérieur, on en déduit que la densité de courant est tangente à la surface.
En prenant une coupe de la configuration précédente, on obtient celle décrite par la figure ci-après.
Dans ce plan de coupe, on choisit arbitrairement un chemin rectangulaire $\mathscr{C}$, de sommets $A$, $B$, $C$ et $D$. Les segments $AB$ et $CD$ sont de longueurs $l_1 = l_2 = l = o(u)$ et portés par un vecteur unitaire quelconque orthogonal à la normale à la surface : ${\bf u_l}\perp{\bf n_{12}}$. Les segments $BC$ et $DA$ sont de longueur $\alpha = o(u^2)$ et portés par ${\bf n_{12}}$.
La loi de Lenz-Faraday appliquée au contour $\mathscr{C}$ et multipliée par un facteur $\displaystyle\frac{1}{l}$, nous donne :
$$\frac{1}{l} \oint_{\mathscr{C}} {\bf e}\cdot{\bf dl} = \frac{1}{l} \iint_{s_l} -\frac{\partial\,{\bf b}}{\partial t}\cdot{\bf d S}$$En développant le membre de gauche, on obtient :
$$\begin{aligned}\frac{1}{l} \oint_{\mathscr{C}} {\bf e}\cdot{\bf dl} = \frac{1}{l} \int_{l_1}{\bf e}\cdot{\bf u_l}\,\text{d} l ~+ &\frac{1}{l} \int_{l_2}{\bf e}\cdot(-{\bf u_l})\,\text{d} l \\ &+ \underbrace{\frac{1}{l}\left(\int_{B}^{C}{\bf e}\cdot{\bf n_{12}}\,\text{d} l + \int_{D}^{A}{\bf e}\cdot(-{\bf n_{12}})\,\text{d} l\right)}_{\displaystyle\lambda}\end{aligned}$$Or :
$$\lambda \leq \left(\max_{BC}( \lVert {\bf e} \rVert) + \max_{DA}( \lVert {\bf e} \rVert)\right)\frac{\alpha}{l} = o(u) $$De plus :
$$\begin{aligned} {\left| \frac{1}{l} \iint_{s_l} -\frac{\partial\,{\bf b}}{\partial t} \cdot {\bf d S} \right|} = \frac{1}{l} & \iint_{s_l} \frac{\partial\,{\bf b}}{\partial t}\cdot({\bf u_l}\wedge{\bf n_{l2}})\,\text{d} S \\ & \leq \max_{s_l}\left( \left\| \frac{\partial\,{\bf b}}{\partial t} \right\| \right)\,\alpha = o(u^2) \end{aligned}$$Finalement, on a :
$$\frac{1}{l} \int_{l_1}{\bf e}\cdot{\bf u_l}\,\text{d} l + \frac{1}{l} \int_{l_2}{\bf e}\cdot(-{\bf u_l})\,\text{d} l + o(u) = o(u^2)$$Par passage à la limite ($u \to 0$), on obtient :
$$({\bf e_1} - {\bf e_2})\cdot{\bf u_l} = 0$$Où ${\bf e_1}$ et ${\bf e_2}$ sont les valeurs du champ électrique au centre du rectangle approché respectivement depuis les milieux $\#1$ ou $\#2$.
La relation précédente étant valable quelque soit le vecteur ${\bf u_l}\perp{\bf n_{12}}$ choisi, on peut en déduire :
$$\boxed{{\bf n_{12}}\wedge({\bf e_2}-{\bf e_1}) = {\bf 0}}$$On a donc continuité de la composante tangentielle du champ électrique.
En appliquant le théorème d’Ampère sur $\mathscr{C}$ et en multipliant chaque terme par un facteur $\frac{1}{l}$ on obtient :
$$\frac{1}{l} \oint_{\mathscr{C}} {\bf h}\cdot{\bf dl} = \frac{1}{l} \iint_{s_l} {\bf j}\cdot{\bf d S}$$Le terme de gauche peut être traité exactement de la même manière que dans le paragraphe précédent.
Intéressons-nous au terme de droite :
$$\begin{aligned}\frac{1}{l} \iint_{s_l} {\bf j}\cdot{\bf d S} &= \frac{1}{l} \iint_{s_l} {\bf j}\cdot({\bf u_l}\wedge{\bf n_{12}})\,\text{d} S \\\ &= \frac{1}{l} \int_l \left(\int_{\frac{-\alpha}{2}}^{\frac{\alpha}{2}} {\bf j}\,\text{d} \alpha\right)\cdot({\bf u_l}\wedge{\bf n_{12}})\,\text{d} l\end{aligned}$$En définissant localement la densité surfacique de courant ${\bf j_s}$ (densité de courant portée par la surface $S$, en $\text{A}\cdot\text{m}^{-1}$) par :
$${\bf j_s} = \lim_{\alpha\to 0}\int_{-\frac{\alpha}{2}}^{\frac{\alpha}{2}} {\bf j}\,\text{d} \alpha ~~~,$$on obtient, par passage à la limite ($u\to 0$) :
$$({\bf h_1}-{\bf h_2})\cdot{\bf u_l} = {\bf j_s}\cdot({\bf u_l}\wedge{\bf n_{12}}) = {\bf u_l}\cdot({\bf n_{12}}\wedge{\bf j_s})$$Où ${\bf h_1}$ et ${\bf h_2}$ sont les valeurs du champ magnétique au centre du rectangle approché respectivement depuis les milieux $\#1$ ou $\#2$, et ${\bf j_s}$ y est également la valeur de la densité surfacique de courant.
Ainsi, pour tout ${\bf u_l}\perp{\bf n_{12}}$ :
$$({\bf h_1}-{\bf h_2})\cdot{\bf u_l} = ({\bf n_{12}}\wedge{\bf j_s})\cdot{\bf u_l}$$La densité superficielle de courant ${\bf j_s}$ étant portée par la surface $S$, et donc étant orthogonale à ${\bf n_{12}}$, on peut en déduire :
$$\boxed{{\bf n_{12}}\wedge({\bf h_2}-{\bf h_1}) = {\bf j_s}}$$On a donc un saut de la composante tangentielle de ${\bf h}$ à travers $S$.
En réalité, une densité surfacique de courant n’existe pas physiquement : on a toujours une densité de courant volumique. Cependant, dans certains cas, l’épaisseur de la zone où circule cette densité volumique est suffisamment faible pour l’assimiler à une densité surfacique. On peut citer deux cas fréquents : dans une plaque ou un tube très mince ; sur les bords d’un conducteur siège de courants induits à « haute fréquence ».
Sans densité superficielle de courant (donc souvent), il y a continuité de la composante tangentielle de ${\bf h}$.
De toutes les relations précédentes (Équations de Maxwell et lois de comportement), on peut normalement déduire les champs ${\bf h}$, ${\bf b}$, ${\bf e}$ et ${\bf d}$ à partir des sources ${\bf j}$ et $\rho_q$. Par contre, ils nous manque des informations si on désire déterminer les énergies et les forces dans le système étudié. C’est ce qui va nous intéresser dans cette partie.
Comme depuis le début, nous allons continuer de travailler à l’échelle locale. Pour ce faire, sortons quelque temps de l’ARQS et introduisons le vecteur de Poynting ${\bf \Pi_p}$, défini par :
$${\bf \Pi_p} = {\bf e}\wedge{\bf h}$$Pour un système électromagnétique contenu dans un volume $V$ de bord $\partial V = S$, la puissance énergétique $P_{\text{em}}$ apportée sous forme électromagnétique au système est donnée par le flux entrant de ${\bf \Pi_p}$, soit :
$$P_{\text{em}} = \oiint_S {\bf \Pi_p}\cdot {\bf n_e}\,\text{d} S$$où ${\bf n_e}$ est la normale entrante à la surface.
Ainsi, en réutilisant nos notations et en appliquant le théorème de la divergence :
$$P_{\text{em}} = - \oiint_S {\bf \Pi_p}\cdot {\bf dS} = \iiint_V \underbrace{-\text{div}\,\left({\bf e}\wedge{\bf h}\right)}_{\displaystyle p_{\text{em}}}\,\text{d} V$$Avec $p_{\text{em}}$ la densité volumique de puissance électromagnétique, qu’on peut développer avec l’identité vectorielle :
$$\text{div}\,\left({\bf e}\wedge{\bf h}\right) = ({\bf rot\,e})\cdot{\bf h}-{\bf e}\cdot({\bf rot\,h})$$Soit :
$$p_{\text{em}} = -{\bf h}\cdot({\bf rot\,e}) + {\bf e}\cdot({\bf rot\,h})$$En réinjectant les équations de Maxwell-Ampère et Maxwell-Faraday :
$$\begin{aligned}\displaystyle p_{\text{em}} &= {\bf h}\cdot\frac{\partial\,{\bf b}}{\partial t}+ {\bf e}\cdot\left({\bf j} + \frac{\partial\,{\bf d}}{\partial t}\right) \\ \displaystyle &= \underbrace{{\bf e}\cdot{\bf j}}_{\displaystyle p_j} + \underbrace{{\bf h}\cdot\frac{\partial\,{\bf b}}{\partial t}}_{\displaystyle p_{mag}} + \underbrace{{\bf e}\cdot\frac{\partial\,{\bf d}}{\partial t}}_{\displaystyle p_{el}}\end{aligned}$$Les différents termes étant :
Le premier terme du bilan de puissance précédent permet de calculer $P_J$ la puissance dissipée par effet Joule (production de chaleur). Il n’existe que dans les domaines conducteurs. Si nous regroupons ces derniers dans ce que nous appellerons $V_c$, on aura :
$$P_J = \iiint_{V_c} {\bf e}\cdot{\bf j}~\text{d} V$$Dans les conducteurs ohmiques, on sait que : ${\bf e} = \rho\,{\bf j}$, ou ${\bf j} = \sigma\,{\bf e}$.
Ainsi, la densité volumique de pertes Joule est :
$$p_j = \rho\, \lVert {\bf j} \rVert^2 = \frac{1}{\sigma}\, \lVert {\bf j} \rVert^2$$Et les pertes totales :
$$\boxed{P_J = \iiint_{V_c} \rho\, \lVert {\bf j} \rVert^2~\text{d} V}$$Dans le cas d’un conducteur portant un courant $I$, on peut directement identifier sa résistance $R$ grâce aux pertes Joule, car $P_J = R\,I^2$. Ainsi :
$$\boxed{R = \frac{1}{I^2} \iiint_{V_c} \rho\, \lVert{\bf j}\rVert^2~\text{d} V}$$En définissant $w_{mag}$ comme la densité d’énergie magnétique contenu dans le volume $V$, on a :
$$p_{mag} = \frac{\partial\,w_{mag}}{\partial t}$$La variation élémentaire de densité volumique d’énergie pendant $\text{d} t$ est ainsi :
$$\text{d}\,w_{mag} = p_{mag}\,\text{d} t = {\bf h}\cdot\frac{\partial\,{\bf b}}{\partial t}\,\text{d} t$$On peut alors en déduire la densité volumique d’énergie magnétique emmagasinée dans le système (à l’instant $t$ correspondant au point de fonctionnement magnétique $({\bf h},{\bf b})$) par :
$$w_{mag} = \int_{t_0}^t {\bf h}\cdot\frac{\partial\,{\bf b}}{\partial t}\,\text{d} t = \int_{0}^{b} {\bf h}\cdot{\bf d b}$$La borne d’intégration inférieure étant choisie arbitrairement pour correspondre à une « énergie stockée nulle » reflétant l’état initial du système.
Dans le cas particulier de milieux où la loi de comportement magnétique est linéaire (${\bf h}=\nu\,{\bf b}$), on a :
$$ w_{mag} = \int_{0}^{b} \nu\,{\bf b}\cdot{\bf d b} = \frac{\nu}{2} \lVert {\bf b} \rVert^2$$Soit :
$$\boxed{w_{mag} = \frac{\lVert {\bf b} \rVert^2}{2\,\mu}}$$Et, dans le cas général :
$$\boxed{w_{mag} = \int_{0}^{b} {\bf h}({\bf b})\cdot{\bf d b}}$$L’énergie magnétique contenue dans le volume $V$ est ainsi :
$$\boxed{W_{mag} = \left\{\begin{aligned} &\iiint_V \frac{\lVert{\bf b}\rVert^2}{2\,\mu} \,\text{d} V , &\text{en linéaire} \\ ~ & ~ \\ &\iiint_V \left(\int_{0}^{b} {\bf h}({\bf b})\cdot{\bf d b}\right) \,\text{d} V , &\text{en non-linéaire}\end{aligned}\right.}$$On peut également définir la densité volumique de coénergie magnétique $\widetilde{w}_{mag}$ par :
$$\widetilde{w}_{mag} + w_{mag} = {\bf h}\cdot{\bf b}$$Alors :
$$\begin{aligned}\text{d}\,\widetilde{w}_{mag} &= \text{d}\left({\bf h}\cdot{\bf b}\right) - \text{d}\,w_{mag} \\ &= {\bf d h}\cdot{\bf b} + \cancel{{\bf h}\cdot{\bf d b}} - \cancel{{\bf h}\cdot{\bf d b}} \\ &= {\bf b}\cdot{\bf d h}\end{aligned}$$On en déduit donc :
$$\widetilde{w}_{mag} = \int_{0}^{h} {\bf b}\cdot{\bf d h}$$Comme précédemment, la borne d’intégration inférieure est choisie arbitrairement pour correspondre à une « énergie stockée nulle » reflétant l’état initial du système.
Dans le cas particulier de milieux où la loi de comportement magnétique est linéaire (${\bf b}=\mu\,{\bf h}$), on a :
$$\widetilde{w}_{mag} = \int_{0}^{h} \mu\,{\bf h}\cdot{\bf d h} = \frac{\mu}{2} \lVert {\bf h} \rVert^2$$Soit :
$$\boxed{\widetilde{w}_{mag} = \frac{\mu\,\lVert{\bf h}\rVert^2}{2}}$$Et, dans le cas général :
$$\boxed{\widetilde{w}_{mag} = \int_{0}^{h} {\bf b}({\bf h})\cdot{\bf d h}}$$La coénergie magnétique contenue dans le volume $V$ est alors :
$$\boxed{\widetilde{W}_{mag} = \left\{\begin{aligned} &\iiint_V \frac{\mu\,\lVert{\bf h}\rVert^2}{2} \,\text{d} V , &\text{en linéaire} \\ ~ & ~ \\ &\iiint_V \left(\int_{0}^{h} {\bf b}({\bf h})\cdot{\bf d h}\right) \,\text{d} V , &\text{en non-linéaire}\end{aligned}\right.}$$Graphiquement, les deux densités volumique peuvent se représenter (dans le cas d’un matériau isotrope) par la figure ci-dessous.
En linéaire, les deux surfaces hachurées ci-dessus sont des triangles de même aire, et on a :
$$\boxed{w_{mag} = \frac{1}{2}\,{\bf h}\cdot{\bf b} = \widetilde{w}_{mag}}$$En linéaire dans l’ARQS, lorsque les sources de champs sont des courants situés dans un sous-volume conducteur $V_c$, nous disposons d’une autre formule de calcul de l’énergie magnétique créée par ces courants reposant sur le potentiel vecteur magnétique.
En effet, d’après ce qui précède, dans ce cas :
On peut développer avec une identité vectorielle :
Ainsi :
En réinjectant Maxwell-Ampère et en appliquant le théorème de la divergence, on obtient :
On notera que l’intégrale volumique se limite au domaine conducteur $V_c$ (seul siège de courants).
Ainsi, pour un système électromagnétique contenu dans un volume choisi suffisamment grand pour que l’influence des sources soient négligeable sur la frontière (${\bf a}$ et/ou ${\bf h}$ nul), on obtient finalement :
Cette formule pourra s’avérer très utile en pratique.
En présence de matériaux durs tels que les aimants, nous ne sommes plus sur une courbe $b(h)$ passant par l’origine, mais les formules générales précédentes restent valables.
Nous allons juste choisir des bornes inférieures différentes :
Graphiquement, on peut ainsi représenter les densités correspondantes sur la figure suivante :
Représentation des densités volumiques d'énergie et coénergie magnétiques dans un aimant
Alors :
$${w}_{mag} = \int_{b_r}^{b} {\bf h}\cdot{\bf d b},~ ~ \text{et} ~ ~\widetilde{w}_{mag} = \int_{-h_c}^{h} {\bf b}\cdot{\bf d h}$$Dans le cas classique de matériaux très durs où la courbe ci-dessus peut être assimilée à une droite, on a ainsi :
$$\begin{aligned}{w}_{mag} &= \int_{b_r}^{b} \frac{\bf (\bf b-b_r)}{\mu_a}\cdot{\bf d b} \\ & = \frac{1}{\mu_a} \left[\frac{({\bf b-b_r})^2}{2}\right]_{br}^b \\ &= \frac{1}{2\,\mu_a} ({\bf b-b_r})^2 = \frac{\mu_a}{2}\,{\bf h}^2 \end{aligned}$$On peut faire de même pour la densité volumique de coénergie :
$$\begin{aligned}\widetilde{w}_{mag} &= \int_{-h_c}^{h} \mu_a\,({\bf h+h_c})\cdot{\bf d h}\\\ &=\mu_a \left[\frac{({\bf h+h_c})^2}{2}\right]_{-hc}^h \\ &= \frac{\mu_a}{2} ({\bf h+h_c})^2 = \frac{{\bf b}^2}{2\,\mu_a}\,\end{aligned}$$Chose un peu surprenante, l’expression des densités d’énergie et coénergie magnétiques dans un aimant permanent est l’inverse de celles dans les autres matériaux. Soit :
$$\boxed{w_{mag} = \frac{\mu_a\, \lVert {\bf h} \rVert^2}{2}}~ ~ \text{et} ~ ~\boxed{\widetilde{w}_{mag}=\frac{ \lVert {\bf b} \rVert^2}{2\,\mu_a}}$$Tout ce que nous avons fait en magnétisme ci-avant, peut se développer également en électrostatique à partie de la densité volumique de puissance $p_{el}$. Plutôt que de tout redévelopper, nous allons procéder par analogie :
Ainsi, en ne considérant que le cas linéaire, nous obtenons, pour les densités volumiques d’énergie et coénergie électrostatiques :
$$\boxed{w_{el} = \frac{ \lVert {\bf d} \rVert^2}{2\,\varepsilon}=\frac{1}{2}\,{\bf e}\cdot{\bf d}=\frac{\varepsilon\, \lVert{\bf e} \rVert^2}{2}=\widetilde{w}_{el}}$$L’énergie électrostatique contenue dans le système est alors :
$$\boxed{W_{el} = \iiint_V \frac{1}{2}\,\varepsilon\,\lVert{\bf e}\rVert^2~\text{d} V = \widetilde{W}_{el}}$$Nous disposons également d’une formule pour l’énergie basée sur le potentiel scalaire électrique $v$. Considérant que la densité volumique de charge $\rho_q$ se situe dans un sous-volume $V_q$ de $V$, on a :
$$\boxed{W_{el} = \frac{1}{2} \iiint_{V_q} \rho_q\,v~\text{d} V}$$On aurait pu donner l’ensemble des relations possibles en considérant le cas non-linéaire, mais nous ne l’avons pas fait par souci de concision. On peut quand même fournir les densités volumiques d’énergie qui pourront, le cas échéant servir de point de départ :
$$w_{el} = \int_{0}^{d} {\bf e}\cdot{\bf d\,d},~ ~ ~ \text{et} ~ ~ ~ ~\widetilde{w}_{el} = \int_{0}^{e} {\bf d}\cdot{\bf d\,e}$$Pour calculer des forces (contre-)électromotrices (f.é.m), déterminer les grandeurs utiles lors d’un couplage avec les équations de circuits électrique, ou tout simplement identifier des inductances équivalentes, nous devrons calculer le flux magnétique traversant nos conducteurs (généralement des enroulements).
Par définition, le flux traversant une surface donnée $S$ est :
$$\varphi = \iint_S {\bf b}\cdot{\bf d S}$$On peut facilement le calculer en utilisant le potentiel vecteur. En effet, en utilisant le théorème de Stokes :
$$\varphi = \iint_S ({\bf rot\,a})\cdot{\bf d S} = \oint\limits_{C =\partial S} {\bf a}\cdot{\bf d l}$$Dans certains cas, cette formule est suffisante. Mais lorsque nous aurons des conducteurs avec des sections ne pouvant plus être considérées comme négligeables, il nous faut déterminer le « flux moyen » le traversant.
Prenons l’exemple simple d’une spire circulaire tel que schématisé ci-dessous :
Exemple simple d'une spire (en haut : circuit filiforme, en bas : géométrie réelle).
Si la section du conducteur $S_c$ est très faible, on peut considérer la spire comme un circuit filiforme (représenté en haut de la figure) et appliquer les formules ci-dessous pour calculer le flux à travers la surface en pointillés bleus.
Mais pour une géométrie plus réaliste telle que la vraie spire dessinée en rouge, on doit déterminer le flux moyen la traversant. Un moyen simple d’y parvenir est d’utiliser la valeur moyenne du potentiel vecteur ${\bf a}$ sur une section de conducteur. Soit :
$$\varphi = \int_C \left(\frac{1}{S_c} \iint_{S_c} {\bf a}~\text{d} S\right)\cdot{\bf dl} $$où $C$ est la longueur centrale de la spire.
On peut ainsi se ramener à l’intégrale volumique sur la spire complète (volume $V_s$) :
$$\varphi = \frac{1}{S_c} \iiint_{V_s} {\bf a}\cdot{\bf u}~\text{d} V$$${\bf u}$ étant un vecteur unitaire dirigé dans le sens du courant (orthogonal à $S_c$).
Pour un conducteur à section constante et portant un courant stationnaire $I$ :
Ainsi :
$$\boxed{\varphi = \frac{1}{I} \iiint_{V_s} {\bf a}\cdot{\bf j} ~\text{d} V}$$L’expression ci-dessus est particulièrement intéressante car elle nous permet de retrouver facilement la relation entre l’énergie magnétique et le flux ou l’inductance. En effet, on a directement :
$$\varphi = \frac{2}{I}\,W_{mag}$$Soit :
$$\boxed{W_{mag} = \frac{1}{2}\,\varphi\,I}$$En linéaire, puisque par définition, l’inductance $L$ du conducteur est :
$$ L = \frac{\varphi}{I}$$On a également l’expression bien connue (valable uniquement en linéaire):
$$\boxed{W_{mag} = \frac{1}{2}\,L\,I^2}$$Initialement énoncée pour une particule portant une charge élémentaire $q$, la force de Lorentz est donnée par :
$${\bf f} = q\,\left({\bf e} + {\bf v}\wedge{\bf b}\right)$$Cette expression peut se généraliser au cas d’une distribution volumique de sources (densités volumiques de charge $\rho_q$ et de courant ${\bf j}$) par :
$${\bf f_{em}} = \underbrace{\rho_q\,{\bf e}}_{\displaystyle{\bf f_C}} + \underbrace{{\bf j}\wedge{\bf b}}_{\displaystyle{\bf f_L}}$$Où la densité volumique de force ${\bf f_{em}}$, se décompose en deux termes :
La formule ci-dessus n’est valable qu’en l’absence de milieu diélectrique ou ferromagnétique. Pour l’utiliser dans le cas général (même dans l’ARQS), il faudrait rajouter les termes correspondant à la densité de charges liées pour la force de Coulomb, et ceux liés à la densité de courants d’aimantation pour celle de Laplace (les courants induits étant pris en compte dans ${\bf j}$). Soit quelque chose qui pourrait ressembler à :
$${\bf f_{em}} = (\rho_q-\text{div}\,{\bf p})\,{\bf e} + ({\bf j} + {\bf rot\,m})\wedge{\bf b}$$Ainsi, en dehors du cas particulier du vide (ou de milieux analogues comme l’air), cette formule est difficilement exploitable. C’est pourquoi nous préférerons utiliser une approche énergétique.
Le système considéré étant fermé (sans échange de matière avec l’extérieur), le premier principe de la thermodynamique nous permet d’écrire la variation de son énergie interne $\mathcal{U}$ :
$$\text{d}\,\mathcal{U} = \text{d}\,\mathcal{Q} + \text{d}\,\mathcal{W}$$avec :
Dans le cas de transformations réversibles, le deuxième principe nous donne :
$$\text{d}\,\mathcal{Q} = T~\text{d}\,\mathcal{S} $$$T$ et $\mathcal{S}$ étant respectivement la température et l’entropie du système.
L’énergie interne est alors une fonction dont une des variables naturelles est $\mathcal{S}$. Pour faciliter la résolution de certains problèmes, les thermodynamiciens ont l’habitude de définir d’autres grandeurs énergétiques en ajoutant des termes à l’énergie interne.
Ainsi, on peut définir l’énergie libre par :
$$\mathcal{F} = \mathcal{U} - T\,\mathcal{S}$$On a alors :
$$\begin{aligned}\text{d}\,\mathcal{F} &= \text{d}\,\mathcal{U} - \text{d}\,(T\,\mathcal{S})\\\ &= \cancel{T~\text{d}\,\mathcal{S}} + \text{d}\,\mathcal{W} - \cancel{T~\text{d}\,\mathcal{S}} - \mathcal{S}~\text{d}\,T\end{aligned}$$Soit :
$$\boxed{\text{d}\,\mathcal{F} = -\mathcal{S}~\text{d}\,T+\text{d}\,\mathcal{W}}$$Revenons à l’échelle locale en définissant les densités volumiques correspondant à nos grandeurs énergétiques :
$$\left\{\begin{array}{c} \displaystyle \mathcal{f} = \frac{\text{d}\,\mathcal{F}}{\text{d}\,V} \\[3ex] \displaystyle \mathcal{s} = \frac{\text{d}\,\mathcal{S}}{\text{d}\,V} \\[3ex] \displaystyle \mathcal{w} = \frac{\text{d}\,\mathcal{W}}{\text{d}\,V}\end{array}\right.$$Les variations temporelles correspondantes sont ainsi :
$$\begin{aligned}\frac{\partial\,\mathcal{f}}{\partial t} &= -\mathcal{s}\, \frac{\partial\,T}{\partial t} + \frac{\partial\,\mathcal{w}}{\partial t}\\[2ex] &= -\mathcal{s}\,\frac{\partial\,T}{\partial t} + \frac{\partial\,\mathcal{w_{em}}}{\partial t} + \frac{\partial\,\mathcal{w_{\ne em}}}{\partial t}\end{aligned}$$en séparant dans $\mathcal{w}$ la partie liée aux électromagnétiques $\mathcal{w_{em}}$ de celle liée aux autres phénomènes physiques $\mathcal{w_{\ne em}}$.
En réinjectant le bilan de puissance vu en début de section, on obtient finalement :
$$\boxed{\frac{\partial\,\mathcal{f}}{\partial t} = -\mathcal{s}\,\frac{\partial\,T}{\partial t} + {\bf e}\cdot{\bf j} + {\bf h}\cdot\frac{\partial\,{\bf b}}{\partial t} + {\bf e}\cdot\frac{\partial\,{\bf d}}{\partial t} + \frac{\partial\,\mathcal{w_{\ne em}}}{\partial t}}$$En considérant que notre système est le siège d’effets purement magnétiques, la variation de la densité d’énergie libre magnétique (qu’on notera dans ce cas $f_m$) ne contient plus que les termes suivants :
$$\text{d}\,\mathcal{f_m} = -\mathcal{s}~ \text{d} \,T + {\bf h}\cdot{\bf d b} = -\mathcal{s}~\text{d}\,T + \text{d}\,w_{mag}$$On remarque donc le lien entre énergie libre et énergie magnétique. À température fixée, les deux sont identiques. En fait, compte tenu de la différence d’ordre de grandeur des constantes de temps thermiques et électromagnétiques, on pourra toujours considérer que l’énergie magnétique telle que nous l’avons définie est l’énergie libre du système.
Dans le cadre d’une transformation à température et flux magnétique constants (et donc à induction constante), nous aurons :
$$\text{d} \,\mathcal{f_m} = 0 ~~ \Rightarrow ~~ \text{d}\,\mathcal{F_m} = 0$$Ainsi, si nous voulons calculer la force ${\bf F_m}$ d’origine magnétique s’exerçant sur un élément du système pouvant se déplacer selon la coordonnée généralisée ${\bf x_0}$ à flux constant, il suffit d’écrire la variation d’énergie correspondante (en utilisant le travail de cette force) :
$$\text{d}\,\mathcal{F_m} = {\bf F_m}\cdot{\bf d x_0} + {\bf grad_{x_0}}\,(\mathcal{F_m})\cdot{\bf d x_0}$$Puisque $\text{d}\,\mathcal{F_m} = 0$, et ${\bf grad_{x_0}}\,(\mathcal{F_m})= {\bf grad_{x_0}}\,(W_{mag})$, on en déduit l’expression de la force :
$${\bf F_m} = - {\bf grad_{x_0}}\,(W_{mag})$$Ou encore, en exprimant la composante dans la direction de ${\bf x_0}$ :
$$\boxed{F_m = - \left(\frac{\partial\,W_{mag}}{\partial\,x_0}\right)_{\varphi~\text{cst}}}$$Dans le cas d’une transformation où ce sont les courants qui sont maintenus constants, il convient de définir un nouveau potentiel thermodynamique $\mathcal{g_m}$ par :
$$\mathcal{g_m} = \mathcal{f_m} - {\bf h}\cdot{\bf b}$$Alors :
$$\text{d}\,\mathcal{g_m} = -\mathcal{s}~\text{d}\,T + \cancel{{\bf h}\cdot{\bf d b}} - \cancel{{\bf h}\cdot{\bf d b}} - {\bf b}\cdot{\bf d h}$$Soit :
$$\text{d}\,\mathcal{g_m} = -\mathcal{s}~\text{d}\,T - {\bf b}\cdot{\bf d h} = -\mathcal{s} ~ \text{d}\,T - \text{d}\,\widetilde{w}_{mag}$$Cette nouvelle grandeur énergétique $\mathcal{G_m}$, de densité volumique $\mathcal{g_m}$ est l’enthalpie libre magnétique (ou énergie libre de Gibbs) du système.
On notera le lien avec la coénergie telle que définie précédemment. En fait nous pourrons considérer que la coénergie magnétique est l’opposée de l’enthalpie libre magnétique (à température fixée).
Dans le cadre d’une transformation à température et courants constants (et donc à champ magnétique constant), nous aurons :
$$\text{d} \,\mathcal{g_m} = 0 ~~ \Rightarrow ~~ \text{d}\,\mathcal{G_m} = 0$$En procédant comme dans le paragraphe précédent pour calculer la force d’origine magnétique s’exerçant sur l’élément considéré, on obtient :
$$\text{d}\,\mathcal{G_m} = {\bf F_m}\cdot{\bf d x_0} + {\bf grad_{x_0}}\,(\mathcal{G_m})\cdot{\bf d x_0}$$Soit :
$${\bf F_m} = - {\bf grad_{x_0}}\,(\mathcal{G_m}) = + {\bf grad_{x_0}}\,(\widetilde{W}_{mag})$$Ou encore, en exprimant la composante dans la direction de ${\bf x_0}$ :
$$\boxed{F_m = + \left(\frac{\partial\,\widetilde{W}_{mag}}{\partial\,x_0}\right)_{I~\text{cst}}}$$Les formules précédentes peuvent aussi être utilisées pour calculer des couples. Si un couple $\Gamma$ d’origine magnétique s’exerce sur notre élément du système suivant la position angulaire généralisée $\alpha_0$. On aura :
$$\Gamma = - \left(\frac{\partial\,W_{mag}}{\partial\,\alpha_0}\right)_{\varphi~\text{cst}} = + \left(\frac{\partial\,\widetilde{W}_{mag}}{\partial\,\alpha_0}\right)_{I~\text{cst}}$$Dans le cas de phénomènes électrostatiques, la densité volumique d’énergie libre électrostatique $\mathcal{f_e}$ est, d’après ce qui précède :
$$\text{d}\,\mathcal{f_e} = -\mathcal{s}~ \text{d} \,T + {\bf e}\cdot{\bf d\,d} = -\mathcal{s}~\text{d}\,T + \text{d}\,w_{el}$$Dans le cadre d’une transformation à température et charges électriques constantes (et donc à induction constante), nous aurons :
$$\text{d} \,\mathcal{f_e} = 0 ~~ \Rightarrow ~~ \text{d}\,\mathcal{F_e} = 0$$Ainsi, si nous voulons calculer la force ${\bf F_m}$ d’origine magnétique s’exerçant sur un élément du système pouvant se déplacer selon la coordonnée généralisée ${\bf x_0}$ à flux constant, il suffit d’écrire la variation d’énergie correspondante (en utilisant le travail de cette force) :
$$\text{d}\,\mathcal{F_e} = {\bf F_e}\cdot{\bf d x_0} + {\bf grad_{x_0}}\,(\mathcal{F_e})\cdot{\bf d x_0}$$Puisque $\text{d}\,\mathcal{F_e} = 0$, et ${\bf grad_{x_0}}\,(\mathcal{F_e})= {\bf grad_{x_0}}\,(W_{el})$, on en déduit l’expression de la force :
$${\bf F_e} = - {\bf grad_{x_0}}\,(W_{el})$$Ou encore, en exprimant la composante dans la direction de ${\bf x_0}$ :
$$\boxed{F_e = - \left(\frac{\partial\,W_{el}}{\partial\,x_0}\right)_{Q~\text{cst}}}$$Dans le cas d’une transformation où ce sont les tensions qui sont maintenues constants, il convient de définir un nouveau potentiel thermodynamique $\mathcal{g_e}$ par :
$$\mathcal{g_e} = \mathcal{f_e} - {\bf e}\cdot{\bf d}$$Alors :
$$\text{d}\,\mathcal{g_e} = -\mathcal{s}~\text{d}\,T + \cancel{{\bf e}\cdot{\bf d\,d}} - \cancel{{\bf e}\cdot{\bf d\,d}} - {\bf d}\cdot{\bf d e}$$Soit :
$$\text{d}\,\mathcal{g_e} = -\mathcal{s} ~ \text{d}\,T - {\bf d}\cdot{\bf d e} = -\mathcal{s}~\text{d}\,T - \text{d}\,\widetilde{w}_{el}$$Cette nouvelle grandeur énergétique $\mathcal{G_e}$, de densité volumique $\mathcal{g_e}$ est l’enthalpie libre électrostatique (ou énergie libre de Gibbs) du système.
Dans le cadre d’une transformation à température et tensions constantes (et donc à champ électrique constant), nous aurons :
$$\text{d} \,\mathcal{g_e} = 0 ~~ \Rightarrow ~~ \text{d}\,\mathcal{G_e} = 0$$En procédant comme dans le paragraphe précédent pour calculer la force d’origine électrostatique s’exerçant sur l’élément considéré, on obtient :
$$\text{d}\,\mathcal{G_e} = {\bf F_e}\cdot{\bf d x_0} + {\bf grad_{x_0}}\,(\mathcal{G_e})\cdot{\bf d x_0}$$Soit :
$${\bf F_e} = -{\bf grad_{x_0}}\,(\mathcal{G_e}) = + {\bf grad_{x_0}}\,(\widetilde{W}_{el})$$Ou encore, en exprimant la composante dans la direction de ${\bf x_0}$ :
$$\boxed{F_e = + \left(\frac{\partial\,\widetilde{W}_{el}}{\partial\,x_0}\right)_{V~\text{cst}}}$$Comme très brièvement dit en début de page, l’expression de la densité volumique de force électromagnétique est difficilement exploitable en pratique.
En fait, il existe plusieurs façons de la formuler et aucune ne fait encore consensus. Quelle que soit l’approche utilisée, il faut considérer les opérateurs différentiels agissant sur nos champs au sens des distributions, et plusieurs termes surfaciques interviennent alors. Ces derniers correspondent aux éventuelles densités surfaciques de sources, mais aussi aux discontinuités de paramètres physiques (perméabilité, réluctivité, aimantation) dues aux changements de milieux.
Néanmoins, dans le cas d’un milieu homogène isotrope, on peut tout de même obtenir une formulation de cette densité volumique de force.
Repartons de la formule :
En y réinjectant les équations de Maxwell-Ampère et Maxwell-Gauss et en utilisant les lois de comportement ${\bf d} = \varepsilon_0\,{\bf e} + {\bf p}$ et ${\bf b} = \mu_0\,({\bf h} + {\bf m)}$, on obtient :
$$\boxed{{\bf f_{em}} = \varepsilon_0\,(\text{div}\,{\bf e})\,{\bf e} + \frac{1}{\mu_0}\,{\bf rot\,b}\wedge{\bf b}}$$Dans le cas de plusieurs milieux (et c’est généralement toujours le cas), il nous faut malheureusement rajouter les termes liés aux interfaces pour calculer la densité de force partout et pouvoir en déduire les forces ou couples s’exerçant globalement sur les pièces mobiles.
Nous disposons cependant d’une méthode plus pratique pour calculer ces forces globales : le tenseur de Maxwell.
Soit $\overline{\overline{{\bf T}}}$ champ de tenseur d’ordre 2, continûment différentiable sur le domaine d’étude $V$, tel que :
$${\bf f_{em}} = \overline{\text{div}}\,\overline{\overline{{\bf T}}}$$Ainsi, la force totale ${\bf F_{em}}$ s’exerçant sur un sous-domaine quelconque $\Omega \subset V$ peut se réduire à l’intégrale surfacique de $\overline{\overline{{\bf T}}}$ sur son bord $\partial \Omega$ via le théorème de la divergence :
$${\bf F_{em}} = \iiint_{\Omega} {\bf f_{em}} ~ \text{d} V = \iiint_{\Omega} \overline{\text{div}}\,\overline{\overline{{\bf T}}}~\text{d} V = \oiint_{\partial \Omega} \overline{\overline{{\bf T}}}\cdot{\bf d S}$$Le tenseur ainsi défini n’est pas unique et peut comprendre des termes difficiles à évaluer, en particulier aux interfaces entre les différents milieux. Mais on peut simplifier grandement les choses en :
Dans ce cas, la formule communément admise pour le champ de tenseur sur la frontière $\partial \Omega$ est :
$$\boxed{\overline{\overline{{\bf T}}} = \varepsilon_0\,{\bf e}\otimes{\bf e} + \frac{1}{\mu_0}\,{\bf b}\otimes{\bf b} - \frac{\delta}{2}\left(\varepsilon_0\,{\bf e}^2 + \frac{{\bf b}^2}{\mu_0}\right)}$$où $\otimes$ est le produit tensoriel et $\delta$ le tenseur unité d’ordre 2.
Les termes du tenseur s’expriment en $\text{N}\cdot\text{m}^{-2}$, ils sont donc homogènes à une pression (on l’appelle souvent « pression électromagnétique »). Les mécaniciens parlent plutôt de « contraintes », c’est pourquoi le tenseur est généralement appelé tenseur des contraintes électromagnétiques de Maxwell.
Dans la suite, nous traiterons séparément les parties magnétiques et électrostatiques.
Si on ne considère que les effets magnétiques, le terme général du tenseur est :
$$\left(T_{ij}\right) = \frac{1}{\mu_0}\,b_i\,b_j - \frac{\delta_{ij}}{2\,\mu_0} {\bf b}^2$$Et sa représentation matricielle dans le repère cartésien est donc :
$$\overline{\overline{{\bf T}}} = \frac{1}{\mu_0} \begin{pmatrix} b_x^2 - \frac{{\bf b}^2}{2} & b_x\,b_y & b_x\,b_z \\\\ b_x\,b_y & b_y^2 - \frac{{\bf b}^2}{2} & b_y\,b_z \\\\ b_x\,b_z & b_y\,b_z & b_z^2 - \frac{{\bf b}^2}{2}\end{pmatrix}$$$\longrightarrow$ Vérifions que cette définition est correcte.
La densité de force dans le milieu considéré est :
$${\bf f_m} = \begin{pmatrix}f_{m_x}\\\\f_{m_y}\\\\f_{m_z}\end{pmatrix}=\frac{1}{\mu_0}\,{\bf rot\,b}\wedge{\bf b}=\frac{1}{\mu_0}\,\begin{pmatrix}b_y\,(\partial_y b_x-\partial_x b_y)+b_z\,(\partial_z b_x - \partial_x b_z)\\\\ b_x\,(\partial_x b_y-\partial_y b_x)+b_z\,(\partial_z b_y - \partial_y b_z)\\\\b_x\,(\partial_y b_x - \partial_x b_y)+b_y\,(\partial_y b_z-\partial_z b_y)\end{pmatrix}$$Calculons la première ligne de $\overline{\text{div}}\,\overline{\overline{\bf T}}$ :
$$\begin{aligned}\big(\overline{\text{div}}\,\overline{\overline{\bf T}}\big)_x &= \frac{1}{\mu_0}\,\big(b_x\,\partial_x b_x - b_y\,\partial_x b_y - b_z\,\partial_x b_z + b_y\,\partial_y b_x + b_x\,\partial_y b_y + b_z\,\partial_z b_x + b_x\,\partial_z b_z\big)\\\\ &= \frac{1}{\mu_0}\,\big(b_x\,\cancel{(\partial_x b_x+\partial_y b_y+\partial_z b_z)}+ b_y\,(\partial_y b_x-\partial_x b_y)+b_z\,(\partial_z b_x - \partial_x b_z)\big)=f_{m_x}\end{aligned}$$Remarque : la simplification vient de $\text{div}\,{\bf b} = 0$. On pourrait faire de même pour la deuxième et troisième ligne.
Au final, on retrouve bien $\overline{\text{div}}\,\overline{\overline{{\bf T}}}={\bf f_m}$, le tenseur est bien défini dans le milieu où se situe la frontière du domaine d’intégration.
Du côté électrostatique, la densité volumique de forces est :
$${\bf f_{e}} = \begin{pmatrix}f_{e_x}\\\\f_{e_y}\\\\f_{e_z}\end{pmatrix}=\varepsilon_0\,\left(\text{div}\,{\bf e}\right)\,{\bf e} = \varepsilon_0\,\left(\partial_x e_x + \partial_y e_y + \partial_z e_z\right)\,\begin{pmatrix}{e_x}\\\\{e_y}\\\\{e_z}\end{pmatrix}$$Si on ne considère que les effets électrostatiques, le terme général du tenseur est :
$$\left(T_{ij}\right) = \varepsilon_0\,e_i\,e_j - \frac{\delta_{ij}}{2}\,\varepsilon_0\,{\bf e}^2$$Et sa représentation matricielle dans le repère cartésien est donc :
$$\boxed{\overline{\overline{{\bf T}}} = \varepsilon_0 \begin{pmatrix} e_x^2 - \frac{{\bf e}^2}{2} & e_x\,e_y & e_x\,e_z \\\\ e_x\,e_y & e_y^2 - \frac{{\bf e}^2}{2} & e_y\,e_z \\\\ e_x\,e_z & e_y\,e_z & e_z^2 - \frac{{\bf e}^2}{2}\end{pmatrix}}$$$\longrightarrow$ Vérifions que cette définition est correcte.
Calculons la première ligne de $\overline{\text{div}}\,\overline{\overline{\bf T}}$ :
$$\begin{aligned}\big(\overline{\text{div}}\,\overline{\overline{\bf T}}\big)_x &= \varepsilon_0\,\big(e_x\,\partial_x e_x - e_y\,\partial_x e_y - e_z\,\partial_x e_z + e_y\,\partial_y e_x + e_x\,\partial_y e_y + e_z\,\partial_z e_x + e_x\,\partial_z e_z\big)\\\\ &= \varepsilon_0\,\big(e_x\,(\partial_x e_x+\partial_y e_y+\partial_z e_z)+ e_y\,\cancel{(\partial_y e_x-\partial_x e_y)}+e_z\,\cancel{(\partial_z e_x - \partial_x e_z)}\big)=f_{e_x}\end{aligned}$$Remarque : les simplifications viennent directement du fait qu’en électrostatique : ${\bf rot\,e} = {\bf 0}$. On pourrait faire de même pour la deuxième et troisième ligne.
Au final, on retrouve bien $\overline{\text{div}}\,\overline{\overline{{\bf T}}}={\bf f_e}$, le tenseur est bien défini dans le milieu où se situe la frontière du domaine d’intégration.
Faisons une petite synthèse de ce qui précède sur un exemple purement théorique.
Considérons un système électromagnétique de volume $V$, contenant des domaines conducteurs $V_c$ parcourus par des courants, des domaines aimantés $V_a$, et des domaines ferromagnétiques $V_m$. Le reste du domaine est constitué d’air. On désire calculer la force magnétique ${\bf F}$ s’exerçant sur l’élément correspondant au volume $V_f$. Une représentation patatoïde est donnée ci-dessous.
Il nous suffit alors de définir un domaine $\Omega$ entourant $V_f$ dont la frontière $\partial\Omega$ passe uniquement dans l’air environnant $V_f$ et d’appliquer :
$${\bf F}=\begin{pmatrix}F_x\\\\ F_y\\\\ F_z\end{pmatrix} = \oiint\limits_{\partial\Omega} \overline{\overline{{\bf T}}}\cdot{\bf d S}$$avec :
$$\overline{\overline{{\bf T}}} = \frac{1}{\mu_0} \begin{pmatrix} b_x^2 - \frac{{\bf b}^2}{2} & b_x\,b_y & b_x\,b_z \\\\ b_x\,b_y & b_y^2 - \frac{{\bf b}^2}{2} & b_y\,b_z \\\\ b_x\,b_z & b_y\,b_z & b_z^2 - \frac{{\bf b}^2}{2}\end{pmatrix}$$Dans le cas particulier où les sources sont des grandeurs variant sinusoïdalement en fonction du temps $t$, nous pourrons utiliser une représentation complexe (dans le domaine fréquentiel) qui permet de s’affranchir de la dépendance à la variable $t$ : on l’appelle la représentation de Fresnel.
Soit $u(t)$ une fonction sinusoïdale de fréquence $f=\frac{1}{T}$ de la forme :
$$u(t) = u_{\text{eff}}\,\sqrt{2}\,\cos(\omega\,t+\varphi)$$où :
Nous pouvons alors définir le complexe $\underline{u}$ tel que : ${u(t) \mapsto \underline{u} ~ , ~ \text{tel que : } u(t) = \mathcal{Re}\left[ \sqrt{2} \, \underline{u} \, \text{e} ^{j \omega t}\right]}$, soit :
$$\boxed{\underline{u} = u_{\text{eff}} \, \text{e}^{j \varphi}}$$.
La valeur efficace et la phase suffisent à définir le nombre complexe représentatif du signal temporel, et sont respectivement son module et son argument.
Ainsi, la dérivation par rapport au temps dans le repère complexe devient une simple multiplication par le terme $j \omega$, c’est-à-dire :
$$\frac{\text{d}\,u}{\text{d} t}(t) \mapsto j \omega \,\underline{u}$$Nous pouvons le vérifier :
$$\begin{aligned}\frac{\text{d}\,u}{\text{d} t}(t) &= - \omega\,u_{\text{eff}}\,\sqrt{2}\,\sin(\omega\,t+\varphi) \\ &= \omega\,u_{\text{eff}}\,\sqrt{2}\,\cos(\omega\,t+\varphi+\frac{\pi}{2}) \\ &= \mathcal{Re}\left[\sqrt{2}\,\omega\,u_{\text{eff}}\, \text{e}^{j(\omega\,t\varphi+\frac{\pi}{2})}\right] \\ &= \mathcal{Re}\left[\sqrt{2}\, j \omega \,u_{\text{eff}} \,\text{e}^{j\varphi}\,\text{e}^{j \omega t}\right] \\ &= \mathcal{Re}\left[\sqrt{2}\, (j \omega \,\underline{u})\,\text{e}^{j \omega t}\right]\end{aligned}$$Il existe plusieurs façons de définir les grandeurs complexes à partir d’un même signal $u$. Nous avons choisi d’utiliser un « invariant de puissance » en choisissant la valeur efficace pour le module de $\underline{u}$ (conservation de la norme 2). Nous aurions pu choisir un « invariant d’amplitude » comme les anglo-saxons mais cela s’avère moins pratique à l’usage.
Les complexes ainsi définis ne doivent pas être une nouveauté pour vous : c’est la transformation utilisée depuis toujours en circuits électriques en régime sinusoïdal.
Si nous considérons un système électromagnétique alimenté par des tensions ou courants sinusoïdaux (de même pulsation $\omega$, pas trop élevée pour rester dans l’ARQS), les champs en résultant seront eux aussi sinusoïdaux en temps. Nous pouvons donc leur appliquer la transformation complexe précédente.
Par exemple, le champ magnétique sera :
$$\begin{aligned}{\bf h}(x,y,z,t) &= {\bf h}(x,y,z)\,\cos(\omega t + \varphi(x,y,z))\\ &= \sqrt{2} ~ {\bf h_{\text{eff}}}(x,y,z)\,\cos(\omega t + \varphi(x,y,z)) \\ &= \mathcal{Re}[\sqrt{2}~{\bf h_{\text{eff}}}(x,y,z)\,\text{e}^{j\varphi(x,y,z)}\,\text{e}^{j\omega\,t}]\end{aligned}$$Soit :
$$\underline{{\bf h}}(x,y,z) = {\bf h_{\text{eff}}}(x,y,z)\,\text{e}^{j\varphi(x,y,z)}$$avec
$${\bf h_{\text{eff}}}(x,y,z) = \begin{pmatrix} h_{x_\text{eff}} \\ h_{y_\text{eff}} \\ h_{z_\text{eff}}\end{pmatrix} = \begin{pmatrix} \sqrt{\frac{1}{T}\,\int_{0}^{T} h_{x}^2(x,y,z,t)\,\text{d} t}\\[0.5em]\sqrt{\frac{1}{T}\,\int_{0}^{T} h_{y}^2(x,y,z,t)\,\text{d} t}\\[0.5em]\sqrt{\frac{1}{T}\,\int_{0}^{T} h_{z}^2(x,y,z,t)\,\text{d} t}\end{pmatrix}$$Toutes les relations que nous avons vu jusqu’à maintenant pourront s’appliquer en complexe. Nous ne traiterons pas le cas de l’électrostatique (car statique comme son nom l’indique).
En régime harmonique, les équations de Maxwell sont :
$$\left\{\begin{aligned}\text{div}\,\underline{\bf b} &= 0 \\ {\bf rot\,\underline{h}} &= \underline{\bf j} \\ {\bf rot\,\underline{e}} &= -j \omega\,\underline{\bf b}\end{aligned}\right.$$En régime harmonique, nous aurons :
$$\begin{aligned} {\bf \underline{b}} &= \mu_0\mu_r\,{\bf \underline{h}},~ ~ ~ {\bf \underline{h}} = \nu\,{\bf \underline{b}} \\ {\bf \underline{j}} &= \sigma\,{\bf \underline{e}},~ ~ ~ {\bf \underline{e}} = \rho\,{\bf \underline{j}}\end{aligned}$$Les aimants étant des sources statiques, ils ne peuvent être considérés dans des problèmes harmoniques.
Nous nous intéresserons principalement qu’aux problèmes où les lois de comportement sont linéaires. En effet, le cas non-linéaire est mal défini. Par exemple si ${\bf h}$ est sinusoïdal, alors ${\bf b}$ ne l’est pas en non-linéaire, et inversement (cf figure ci-dessous). Dans ces cas, on ne considère que les valeurs du fondamental des grandeurs (hypothèse du premier harmonique).
Il existe des techniques permettant de traiter ce type de problèmes de façon approchée (en prenant également l’influence de la largeur du cycle d’hystérésis) en définissant par exemple des perméabilités complexes qui permettent de déduire la valeur du fondamental de l’induction associée (via son module) et le déphasage (via son argument), mais nous ne les traiterons pas dans le cadre de ce cours.
Nos potentiels peuvent se définir en complexe également :
$$\begin{aligned}{\bf \underline{b}} &= {\bf rot\,\underline{a}} \\ {\bf \underline{e}} &= -{\bf grad}\,\underline{v} - j\omega\,{\bf \underline{a}}\\ {\bf \underline{j}} &= {\bf rot\,\underline{t}} \\ {\bf \underline{h}} &= {\bf \underline{t}}-{\bf grad}\,\underline{\phi}\end{aligned}$$Dans l’expression ci-dessus, on a utilisé $^{\*}$ comme notation pour le complexe conjugué, et on a fait apparaître le produit hermitien. On pourrait facilement montrer que la densité de pertes ainsi définie correspond à :
$$p_j = \left< p_j(t) \right> = \frac{1}{T}\int_0^T \rho\,\lVert {\bf j}(x,y,z,t) \rVert^2\,\text{d}t$$Dans le cas d’une alimentation alternative mais non-sinusoïdale, le problème peut être résolu par décomposition en série de Fourier, puis en traitant harmonique par harmonique (résolution multi-harmoniques).
Calculer la résistance d’une barre de cuivre de section circulaire $S$ et de longueur $L$ alimentée par une tension $U$ constante.
Calculer les valeurs de capacités :
Les deux étant constitués de deux électrodes séparées par un milieu diélectrique de permittivité relative $\varepsilon_r$.
On considère une bobine rectangulaire à $N$ spires, de longueur $L$ très grande devant sa largeur, telle que représentée ci-dessous :
On désire calculer le champ d’induction créé par cette bobine, sans négliger le rayon des conducteur (on ne sera donc pas dans le cas d’un circuit filiforme).
Pour ce faire, nous procéderons par étapes :
Après avoir précisé les hypothèses, calculer le champ crée par un conducteur cylindrique de rayon $r_c$ parcouru par un courant I (ligne unifilaire).
Dans le langage de votre choix, tracer le champ et les lignes correspondantes :
En déduire celui crée par le conducteur précédent et son retour séparés d’une distance $d_{if}$.
En déduire finalement le champ total créé par la bobine et retrouver le tracé de la solution :
On considère une région de l’espace où règne un champ d’induction uniforme (imposé par une source extérieure) de la forme : ${\bf b} = B_{0}\,{\bf u_x}$.
On place à l’intérieur de cette zone un cylindre de rayon $R$ constitué d’un matériau de perméabilité relative $\mu_r$ constante.
Utiliser le langage de votre choix : MATLAB (ou GNU/Octave), Python, Julia, etc…
On considère un câble coaxial portant un courant $I$ tel que représenté sur la figure ci-dessous.
Le code MATLAB permettant de faire le calcul et le tracé est donné ci-dessous à titre indicatif :
On considère le transformateur élémentaire donné par la figure ci-dessous. Le primaire (en bleu) est constitué de $N_1$ spires, et le secondaire (en rouge) en comporte $N_2$.
Le code MATLAB permettant de faire le calcul et le tracé est donné ci-dessous :
On considère l’électroaimant représenté sur la figure ci-dessous :
Considérons tout d’abord que le circuit magnétique est infiniment perméable (hypothèse du $\mu_{\infty}$). Calculer la force $F$ s’exerçant sur la partie mobile (plaque du bas) par les dérivées de l’énergie et la coénergie, et mettre en évidence la différence de signe entre les deux expressions. Retrouver ensuite la valeur de la force par le tenseur de Maxwell.
Recalculer la force par le tenseur mais dans le cas où le circuit magnétique a une loi de comportement linéaire (perméabilité $\mu_r$). En déduire l’expression de la force de collage (à entrefer $e$ nul).
En utilisant les valeurs de paramètres et de courbes $B(H)$ données ci-dessous, superposer les tracés de la force en fonction du courant ($F(I_c)$) dans les 2 cas : linéaire et non-linéaire.
Cet exercice est adapté de la ventouse magnétique de fermeture de porte des TP FEMM de l’EC S8-B11-EC2.
Le code MATLAB pour faire le calcul :
Le but de l’exercice est de calculer l’inductance de fuite et la résistance d’une barre de machine asynchrone à cage en fonction de la fréquence du courant rotorique. Lorsque le rayon d’alésage de la machine est assez grand, le problème peut se ramener à un problème 2D cartésien tel que représenté par la figure :
Le courant I circulant dans la barre est supposé continu ($I = I_0$). Après avoir précisé vos hypothèses :
Le courant I est maintenant supposé sinusoïdal : $I(t) = I_{\text{eff}}\,\sqrt{2}\,\sin(\omega\,t)$.
Densité volumique de courant dans l'encoche à 50 Hz
Évolution de la valeur de l'inductance de fuite et de la résistance de la barre en fonction de la fréquence
La barre de cuivre de l’exercice Ex.1 est désormais alimentée en courant alternatif.
Calculer la densité de courant lorsque la barre est parcourue par un courant $I(t) = \sqrt{2}\,I_{\text{eff}}\,\cos(\omega\,t)$.
On utilise la transformation en complexes pour résoudre l’équation locale satisfaite par $\underline{\bf j}$ :
On considère une plaque d’acier (de perméabilité $\mu$ et de conductivité $\sigma$) de longueur $L_z$ et de largeur $l_y$ très grandes devant son épaisseur $2\,e$. Cette plaque est placée à l’intérieur d’un inducteur correspondant à un solénoïde de longueur $L$ à $N$ spires jointives de forme rectangulaire.
Celui-ci est parcouru par un courant $I(t)$ imposé. Une représentation schématique du dispositif et le repère associé sont donnés par la figure ci-dessous :
Dans un premier temps, on considère uniquement l’inducteur.
En précisant vos hypothèses, et en utilisant le théorème d’Ampère sur 3 contours judicieusement choisis, montrez que le champ
magnétique à l’intérieur du solénoïde est uniforme et s’écrit : $ {\bf h_s} = n\,I(t)\,{\bf u_z} $.
On explicitera la valeur de $n$ en fonction de $N$ et $L$.
Le courant d’alimentation est : $ I(t) = I_{\text{eff}} \sqrt{2}\,\cos(\omega\,t) $. Donner la forme complexe du courant $\underline{I}$, ainsi que du champ précédent $\underline{h_s}$.
Désormais, on insère la plaque dans l’inducteur. Donner, en le justifiant, l’équation générale de diffusion à laquelle obéit le champ magnétique à l’intérieur de celle-ci ( $x \in [−e; e]$).
À partir de l’équation précédente et de la question 2, donner la forme complexe du système permettant de calculer le champ complexe $\underline{\bf h}(x)$ dans la plaque. Soit : une équation différentielle ordinaire d’ordre 2 à coefficients constants complexes et 2 conditions aux limites.
Résoudre ce système (on fera intervenir l’épaisseur de peau $\delta$).
En déduire l’expression complexe de la densité de courant dans la plaque, puis celle de la densité volumique de puissance moyenne $p$ injectée dans la plaque par l’inducteur, et donc de la puissance moyenne totale $P$ injectée dans la plaque.
Où nous essaierons de faire le lien avec les notions de base vues dans l’EC de S7 Distributions et Équations aux Dérivées Partielles (Équations aux Dérivées Partielles), et définirons les notions théoriques dont nous aurons besoin.
Quelques remarques rapides.
Vers les formulations faibles de nos problèmes…
Approximation sur un partitionnement du domaine.
De toutes les méthodes numériques pouvant être utilisées pour résoudre les équations aux dérivées partielles relevant de l’électromagnétisme basse fréquence, la méthode des éléments finis (MEF) est sans nul doute la plus répandue.
Même si elle fut initialement développée pour la mécanique des structures, certains de ses avantages la rendent particulièrement bien adaptée à nos problèmes :
Il existe de nombreux codes de calcul par éléments finis (EF) dédiés à l’électromagnétisme BF. On peut citer par exemple (liste non exhaustive) :
Commerciaux / Licences propriétaires :
Freeware :
Libres (non spécifiquement dédiés à l’électromagnétisme) :
Dans les travaux pratiques de cet EC, nous utiliserons ONELAB.
Soit un domaine $\Omega$ de $\mathbb{R}^3$, de bord $\Gamma = \partial\Omega$ tel que $\Gamma = \Gamma_d \cup \Gamma_n$, comme présenté sur la figure suivante :
Nous considérons l’espace $\text{L}^2(\Omega)$ des champs scalaires $u$ de $\Omega$ défini par :
$$\text{L}^2(\Omega) = \left\{ u : \iiint_\Omega u^2({\bf x})\,\text{d} \Omega < \infty \right\}$$Il est naturellement muni d’un produit scalaire (produit interne), défini par :
$$(u,v)_{\Omega} = \iiint_\Omega u({\bf x})\,v({\bf x})\,\text{d} \Omega $$auquel est associée la norme :
$$\| u \|_{\text{L}^2(\Omega)} = (u,u)_{\Omega}^{\frac{1}{2}} = \left(\iiint_\Omega u^2({\bf x})\,\text{d} \Omega\right)^{\frac{1}{2}} $$Pour les champs vectoriels ${\bf u}$ de $\Omega$, nous pouvons faire de même avec l’espace $\textbf{L}^2(\Omega)$ :
$$\textbf{L}^2(\Omega) = \left\{ {\bf u} : \iiint_\Omega \|{\bf u}({\bf x})\|^2\,\text{d} \Omega < \infty \right\}$$Et le produit scalaire et la norme associée :
$$({\bf u},{\bf v})_{\Omega} = \iiint_\Omega {\bf u}({\bf x})\cdot{\bf v}({\bf x})\,\text{d} \Omega$$$$\| {\bf u} \|_{\textbf{L}^2(\Omega)} = ({\bf u},{\bf u})_{\Omega}^{\frac{1}{2}} = \left(\iiint_\Omega \| {\bf u}({\bf x})\|\,\text{d} \Omega\right)^{\frac{1}{2}} $$
Concrètement, ces espaces sont des espaces de Hilbert et correspondent aux champs de $\Omega$ à « énergie finie ». Ils conviennent donc parfaitement à nos champs électromagnétiques (et même à tout champ physique).
À partir des espaces précédents, on défini d’autres espaces d’intérêt par :
$$\text{H}^1(\Omega) = \left\{ u \in \text{L}^2(\Omega) : \partial_x u, \partial_y u,\partial_z u \in \text{L}^2(\Omega)\right\}$$$$\textbf{H}^1(\Omega) = \left\{ {\bf u} \in \textbf{L}^2(\Omega) : \partial_x {\bf u}, \partial_y {\bf u},\partial_z {\bf u} \in \textbf{L}^2(\Omega)\right\}$$Où $\partial_i$ désigne la dérivée partielle selon la coordonnée $i$ prise au sens des distributions.
Nos champs étant liés aux opérateurs ${\bf grad}$, $\text{div}$ et ${\bf rot}$, nous utiliserons plutôt les espaces définis ci-dessous :
$$\text{H}({\bf grad},\Omega) = \left\{ u \in \text{L}^2(\Omega) : {\bf grad}\,u \in \text{L}^2(\Omega)\right\} = \text{H}^1(\Omega)$$$$\textbf{H}({\bf rot},\Omega) = \left\{ {\bf u} \in \textbf{L}^2(\Omega) : {\bf rot\,u} \in \textbf{L}^2(\Omega)\right\}$$$$\textbf{H}(\text{div},\Omega) = \left\{ {\bf u} \in \textbf{L}^2(\Omega) : \text{div}\,{\bf u} \in \text{L}^2(\Omega)\right\}$$On peut remarquer que $\text{H}({\bf grad},\Omega)$ est $\text{H}^1(\Omega)$ (c’est juste la notation qui change), ainsi que :
Finalement, on voit que nos espaces sont liés par ce qu’on appelle le complexe de De Rham :
$$\text{H}^1(\Omega) = \text{H}({\bf grad},\Omega) \xrightarrow[]{{\bf grad}} \textbf{H}({\bf rot},\Omega) \xrightarrow[]{{\bf rot}} \textbf{H}(\text{div},\Omega) \xrightarrow[]{\text{div}} \text{L}^2(\Omega) $$Dans le cadre de ce cours, nous nous intéresserons principalement à deux types de problèmes détaillés ci-après. Nous verrons dans le chapitre 3 que nous pourrons toujours nous ramener à l’une ou l’autre de ces deux formulations.
Dans la suite, nous prendrons de gros raccourcis pour ne traiter que les aspects qui nous intéresseront pour résoudre nos problèmes. Vous trouverez un cadre et un développement mathématique plus rigoureux, en particulier sur les régularités des fonctions intervenant dans les formulations, dans le polycopié de votre cours de S7 (ou une version plus ancienne telle que [Anto21]).
Dans un premier temps, nous considérons le problème suivant. On cherche $u$ champ scalaire de $\Omega$ vérifiant la formulation forte :
$$\left\{\begin{aligned}\text{div}\left(\alpha\,{\bf grad}\,u\right) + \beta &= 0~, &\text{dans}~ \Omega \\ \left.u\right|_{\Gamma_d} &= 0~, &\text{sur}~ \Gamma_d \\ \left.\frac{\partial\,u}{\partial {\bf n}}\right|_{\Gamma_n} = \left.{\bf grad}\, u \cdot {\bf n} \right|_{\Gamma_n} &= \gamma~, &\text{sur}~ \Gamma_n \end{aligned}\right.$$où $\alpha, \beta \in \text{L}^2({\Omega})$ continûment différentiables sur $\Omega$ (voire $\overline{\Omega}$), et $\gamma \in \mathcal{C}^1(\Gamma_n)$.
Soit $v$ champ scalaire quelconque de $\text{H}({\bf grad},\Omega)$, qu’on appelle fonction test. En multipliant l’équation de la formulation forte et en intégrant sur tout le domaine, on obtient :
$$\iiint_{\Omega} \left[\text{div}\left(\alpha\,{\bf grad}\,u\right)\right]\,v~\text{d}\Omega + \iiint_{\Omega} \beta\,v ~\text{d}\Omega =0$$En utilisant la formule de Green en grad-div :
$$ \text{div}(a\,{\bf b}) = ({\bf grad}\,a)\cdot{\bf b} + (\text{div}\,{\bf b})\,a $$on peut développer le terme de gauche pour obtenir (intégration par parties) :
$$\iiint_{\Omega} -(\alpha\,{\bf grad}\,u)\cdot{\bf grad}\,v~\text{d}\Omega + \iiint_{\Omega} \text{div}(v\,\alpha\,{\bf grad}\,u)~\text{d}\Omega + \iiint_{\Omega} \beta\,v ~\text{d}\Omega =0$$On applique alors le théorème de la divergence :
$$\iiint_{\Omega} -(\alpha\,{\bf grad}\,u)\cdot{\bf grad}\,v~\text{d}\Omega + \oiint_{\Gamma} v\,\alpha\,{\bf grad}\,u \cdot {\bf n}~\text{d}\Gamma + \iiint_{\Omega} \beta\,v ~\text{d}\Omega =0$$On peut décomposer le terme surfacique en deux en utilisant la condition aux limites de Neumann pour obtenir :
On aboutit finalement à l’expression générale :
$$\iiint_{\Omega} (\alpha\,{\bf grad}\,u)\cdot{\bf grad}\,v~ \text{d}\Omega - \iint_{\Gamma_n} v\,\alpha\,\gamma~ \text{d}\Gamma - \iiint_{\Omega} \beta\,v ~ \text{d}\Omega = \iint_{\Gamma_d} v\,\alpha\,{\bf grad}\,u \cdot {\bf n}~ \text{d}\Gamma$$N’ayant aucune information sur le terme ${\bf grad}\,u$ sur $\Gamma_d$, ce terme peut être gênant pour déterminer $u$, mais nous pouvons nous en affranchir en choisissant la fonction test $v$ nulle sur $\Gamma_d$. C’est-à-dire en prenant $v \in \text{H}_0({\bf grad},\Omega)$ définit par :
$$\boxed{\text{H}_0({\bf grad},\Omega) = \{ u \in \text{H}({\bf grad},\Omega) : u\|_{\Gamma_d} = 0\} =\text{H}_0^1(\Omega)}$$Nous pouvons remarquer que $u \in \text{H}_0({\bf grad},\Omega)$ aussi.
Nous obtenons alors la formulation faible (ou variationnelle) de notre problème :
Cherchons $u \in \text{H}_0({\bf grad},\Omega)$, tel que :
$$\boxed{\forall v \in \text{H}_0({\bf grad},\Omega), \iiint_{\Omega} (\alpha\,{\bf grad}\,u)\cdot{\bf grad}\,v~ \text{d}\Omega - \iiint_{\Omega} \beta\,v ~ \text{d}\Omega - \iint_{\Gamma_n} \alpha\,\gamma\,v~ \text{d}\Gamma = 0}$$Pour simplifier la notation, on peut utiliser les produits scalaires définis en début de page, et le produit surfacique $\left<\cdot,\cdot\right>$ défini par :
$$\left< u,v \right>_{\Gamma} = \iint_{\Gamma} u\,v\,\text{d}\Gamma$$Notre formulation faible s’écrit alors :
$$\boxed{\forall v \in \text{H}_0({\bf grad},\Omega), \left(\alpha\,{\bf grad}\,u,{\bf grad}\,v\right)_{\Omega} - \left(\beta,v\right)_{\Omega} - \left<\alpha\,\gamma,v\right>_{\Gamma_n} = 0}$$On pourrait alors montrer (à l’aide du théorème de Lax-Milgram) que le problème correspondant à notre formulation faible est bien posé (existence et unicité de la solution), et qu’il est équivalent (à l’aide du théorème de trace) à la formulation forte, qui, elle aussi, correspond donc à un problème bien posé [Anto21], [Kern04], [Spit02].
Nous pouvons remarquer que, selon leur nature, les conditions aux limites sont traitées de façons totalement différentes.
Conditions de Dirichlet : Elles sont directement prises en compte via la définition de l’espace fonctionnel $\text{H}_0({\bf grad},\Omega)$ (espace des fonctions tests et de recherche de la solution).
Conditions de Neumann : Elles sont prises en compte via le terme surfacique $\left<\alpha\,\gamma,v\right>_{\Gamma_n}$ dans la formulation faible. Et en l’absence de ce terme, ce sont des conditions de Neumann homogènes qui s’appliquent sur $\Gamma_n$.
C’est pourquoi, en l’absence d’une condition particulière sur un bord du domaine, la condition aux limite par défaut est une condition de Neumann homogène $\frac{\partial\,u}{\partial {\bf n}} = 0$.
Considérons maintenant la formulation forte suivante. Cherchons un champ vectoriel ${\bf u}$ de $\Omega$ vérifiant :
$$\left\{\begin{aligned}{\bf rot}\left(\alpha\,{\bf rot}\,{\bf u}\right) + \boldsymbol{\beta} & = {\bf 0}~ , &\text{dans}~ \Omega\ \left.{\bf u}\wedge{\bf n}\right|_{\Gamma_d} &= {\bf 0}~ , &\text{sur}~ \Gamma_d \ \left.\frac{\partial\,{\bf u}}{\partial {\bf n}}\right|_{\Gamma_n} = \left.{\bf rot\,u}\wedge{\bf n}\right|_{\Gamma_n} &= \boldsymbol{\gamma}~ , &\text{sur}~ \Gamma_n \end{aligned}\right.$$où $\alpha \in \text{L}^2({\Omega})$, $\boldsymbol{\beta} \in \textbf{L}^2({\Omega})$ continûment différentiables sur $\Omega$ (voire $\overline{\Omega}$), et $\gamma \in \mathcal{C}^1(\Gamma_n)$.
Soit la fonction test ${\bf v}$ champ vectoriel quelconque de $\textbf{H}({\bf rot},\Omega)$. En prenant le produit scalaire de l’équation de la formulation forte par ${\bf v}$ et en intégrant sur tout le domaine, on obtient :
$$\iiint_{\Omega} \left({\bf rot}\left(\alpha\,{\bf rot}\,{\bf u}\right)\right)\cdot{\bf v}~ \text{d} \Omega + \iiint_{\Omega}\boldsymbol{\beta}\cdot{\bf v}~ \text{d} \Omega = {\bf 0}$$En utilisant la formule de Green en rot-rot :
$$\text{div}({\bf a}\wedge{\bf b}) = {\bf b}\cdot{\bf rot\,a} - {\bf a}\cdot{\bf rot\,b}$$on peut intégrer par parties le terme de gauche pour obtenir :
$$\iiint_{\Omega} \left(\alpha\,{\bf rot}\,{\bf u}\right)\cdot{\bf rot\, v}~ \text{d} \Omega - \iiint_{\Omega} \text{div}\,\left(\alpha\,{\bf rot\,u}\wedge{\bf v}\right)~ \text{d} \Omega + \iiint_{\Omega}\boldsymbol{\beta}\cdot{\bf v}~ \text{d} \Omega = {\bf 0}$$Après application du théorème de la divergence :
$$\iiint_{\Omega} \left(\alpha\,{\bf rot}\,{\bf u}\right)\cdot{\bf rot\, v}~ \text{d} \Omega - \oiint_{\Gamma} \left(\alpha\,{\bf rot\,u}\wedge{\bf v}\right)\cdot{\bf n}~ \text{d} \Gamma + \iiint_{\Omega}\boldsymbol{\beta}\cdot{\bf v}~ \text{d} \Omega = {\bf 0}$$En séparant les termes surfaciques et en remarquant que :
$$\begin{aligned}({\bf rot\,u}\wedge{\bf v})\cdot{\bf n} &= {\bf n}\cdot({\bf rot\,u}\wedge{\bf v}) = {\bf rot\,u}\cdot({\bf v}\wedge{\bf n}) \\ &= ({\bf n}\wedge{\bf rot\,u})\cdot{\bf v} = - ({\bf rot\,u}\wedge{\bf n})\cdot{\bf v}\end{aligned}$$on obtient :
$$\begin{aligned}\iiint_{\Omega} \left(\alpha\,{\bf rot}\,{\bf u}\right)&\cdot{\bf rot\, v}~ \text{d} \Omega + \iiint_{\Omega}\boldsymbol{\beta}\cdot{\bf v}~ \text{d} \Omega \\ &+ \iint_{\Gamma_n} \alpha\,\underbrace{\left({\bf rot\,u}\wedge{\bf n}\right)}_{\displaystyle \boldsymbol{\gamma}}\cdot{\bf v}~ \text{d} \Gamma = \iint_{\Gamma_d} \alpha\,{\bf rot\,u}\cdot\left({\bf v}\wedge{\bf n}\right)~ \text{d} \Gamma \end{aligned}$$Comme précédemment, on peut éliminer le second terme surfacique en prenant la fonction test ${\bf v}$ et la fonction inconnue ${\bf u}$ dans $\textbf{H}_0({\bf rot},\Omega)$ défini par :
$$\boxed{ \textbf{H}_0({\bf rot},\Omega) = \{{\bf v} \in \textbf{H}({\bf rot},\Omega) : \left.{\bf v}\wedge{\bf n}\right|_{\Gamma_d} = {\bf 0} \} }$$On aboutit à la formulation faible de notre problème :
Cherchons ${\bf u} \in \textbf{H}_0({\bf rot},\Omega)$ tel que :
$$\boxed{\begin{aligned}\forall\,{\bf v} \in \textbf{H}_0({\bf rot},\Omega),~ \iiint_{\Omega} &\left(\alpha\,{\bf rot}\,{\bf u}\right)\cdot{\bf rot\, v}~ \text{d} \Omega \ &+ \iiint_{\Omega} \boldsymbol{\beta}\cdot{\bf v}~ \text{d} \Omega + \iint_{\Gamma_n} \alpha\,\boldsymbol{\gamma}\cdot{\bf v}~ \text{d} \Gamma = 0\end{aligned}}$$Ou encore :
$$\boxed{\forall\,{\bf v} \in \textbf{H}_0({\bf rot},\Omega),~ \left(\alpha\,{\bf rot\,u},{\bf rot\,v}\right)_{\Omega} + \left(\boldsymbol{\beta},{\bf v}\right)_{\Omega} + \left<\alpha\,\boldsymbol{\gamma},{\bf v}\right>_{\Gamma_n} = 0}$$Les remarques précédentes sur les conditions aux limites restent valables :
Attention : la formulation rot-rot ainsi définie ne correspond pas obligatoirement à un problème bien posé. C’est le cas en 2D, mais en 3D nous devrons rajouter une condition de jauge (comme nous le verrons dans la fin du chapitre).
Dans les deux cas précédents, on peut écrire chaque formulation de la même façon, donnée par :
$$\displaystyle \forall v \in \text{H}_0({\bf grad},\Omega), ~ \underbrace{\left(\alpha\,{\bf grad}\,u,{\bf grad}\,v\right)_{\Omega}}_{\displaystyle a(u,v)} = \underbrace{\left(\beta,v\right)_{\Omega} + \left<\alpha\,\gamma,v\right>_{\Gamma_n}}_{\displaystyle \displaystyle L(v)}$$$$\forall\,{\bf v} \in \textbf{H}_0({\bf rot},\Omega),~ \underbrace{\left(\alpha\,{\bf rot\,u},{\bf rot\,v}\right)_{\Omega}}_{\displaystyle a({\bf u},{\bf v})} = \underbrace{- \left(\boldsymbol{\beta},{\bf v}\right)_{\Omega} - \left<\alpha\,\boldsymbol{\gamma},{\bf v}\right>_{\Gamma_n}}_{\displaystyle \displaystyle L({\bf v})} $$Où $a(\cdot,\cdot)$ est une forme bilinéaire continue et coercive, et $L(\cdot)$ une forme linéaire continue, définies sur $V \left(= \text{H}_0({\bf grad},\Omega) ~ \text{ou}~ \textbf{H}_0({\bf rot},\Omega)\right)$.
Pour simplifier, on ne conservera que la notation correspodant au premier cas ($u$ et $v$) mais ce que nous dirons sera aussi valable pour le deuxième (avec ${\bf u}$ et ${\bf v}$).
En fait, ce qu’on appelle formulation variationnelle du problème est cette forme, qu’on écrit :
$$\boxed{\left\{\begin{aligned}{\large \mathstrut} &\text{Trouver}~ u \in V ~ \text{tel que :} \\ &a(u,v) = L(v),~ ~\forall\,v \in V\end{aligned}\right.}$$La forme bilinéaire $a$ étant symétrique, résoudre la forme variationnelle est équivalent à trouver le minimum de la fonctionnelle énergétique (quadratique) définie par :
$$\boxed{J(v) = \frac{1}{2}\,a(v,v) - L(v)}$$Soit :
$$\boxed{u = \mathop{\text{argmin}}_{v \in V} J(v)}$$En effet, on a :
$$\begin{aligned}J(u+v)-J(u) &= \frac{1}{2}\,a(u+v,u+v) - L(u+v) - \frac{1}{2}\,a(u,u) + L(u) \\ &= \cancel{\frac{1}{2}\,a(u,u)} + \underbrace{\frac{1}{2}\,a(u,v) + \frac{1}{2}\,a(v,u)}_{\displaystyle a(u,v)} + \frac{1}{2}\,a(v,v) - \cancel{L(u)} - L(v) - \cancel{\frac{1}{2}\,a(u,u)} + \cancel{L(u)} \\ &= \underbrace{a(u,v) - L(v)}_{\displaystyle \text{linéaire en}~v} + \underbrace{\frac{1}{2}\,a(v,v)}_{\displaystyle o(v)} \end{aligned}$$On en déduit :
$$\text{d} J(u) v = a(u,v) - L(v)$$La différentielle de $J$ s’annule donc en $u$, solution de notre problème (ce qui correspond bien au minimum de $J$).
On peut donc résoudre notre formulation variationnelle en cherchant à minimiser la fonctionnelle énergétique (ou une approximation de cette dernière), mais en pratique, on préfère une autre méthode pour résoudre la formulation faible : la méthode de Galerkin (dont la MEF est un cas particulier).
L’idée de base derrière la méthode des éléments finis est l’approche classique consistant à décomposer un problème compliqué en une multitude de petits problèmes beaucoup plus simples à résoudre :
Pour illustrer les différents aspects abordés dans la suite, on utilisera l’exemple simple donné par la figure suivante :
Nous avons donc un problème de Poisson avec conditions de Dirichlet et Neumann, dont la formulation forte peut s’écrire :
$$\left\{\begin{aligned} \text{div} \left({\bf grad}\,u \right) + \beta &= 0, ~ \text{dans} ~ \Omega \\ u|_{\Gamma_d} &= 0 \\ \partial_n u|_{\Gamma_n} &= \gamma = 1 \end{aligned}\right.$$avec $\beta$ définie par morceaux par :
$$\left\{\begin{aligned} \beta &= 1 ~ ~ \text{ dans } ~ ~ \Omega_s \\ \beta &= 0 ~ ~ \text{ dans } ~ ~ \Omega/\Omega_s \end{aligned}\right.$$Et d’après la section précédente, la formulation faible est :
$$\begin{aligned}&\text{Trouver} ~ u \in \text{H}_0({\bf grad},\Omega) = \left\{ u \in \text{H}({\bf grad},\Omega) : u|_{\Gamma_d} = 0\right\} ~ \text{tel que :}\\\\ &\forall v \in \text{H}_0({\bf grad},\Omega), \iiint_{\Omega} {\bf grad}\,u \cdot{\bf grad}\,v ~\text{d}\Omega - \iiint_{\Omega_s} v ~ \text{d}\Omega - \iint_{\Gamma_n} v ~ \text{d}\Gamma = 0\end{aligned}$$Soit :
$$\boxed{\begin{aligned}&\text{Trouver} ~ u \in \text{H}_0({\bf grad},\Omega) ~\text{tel que :}\\\\ &\forall v \in \text{H}_0({\bf grad},\Omega), ~ ~ \underbrace{\left({\bf grad}\,u,{\bf grad}\,v\right)_{\Omega} }_{a(u,v)}+ \underbrace{\left(-1,v\right)_{\Omega_s} + \left<-1,v\right>_{\Gamma_n}}_{-L(v)} = 0\end{aligned}}$$$\rightsquigarrow$ C’est cette dernière que nous résoudrons par la suite.
Nous avons vu dans la section précédente que la forme faible de chaque problème considéré pouvait se mettre sous la forme :
$$\left\{\begin{aligned}~ &\text{Trouver}~ u \in V ~ \text{tel que :}\\ &a(u,v) = L(v),~ ~ \forall\,v \in V\end{aligned}\right.$$$V$ étant de dimension infinie, un tel problème est particulièrement difficile à résoudre. La méthode de Galerkin consiste à résoudre la formulation faible dans un sous-espace $V_h \subset V$ de dimension finie $n_h$ possédant les mêmes propriétés, soit :
$$\left\{\begin{aligned} ~ &\text{Trouver}~ u_h \in V_h ~ \text{tel que :} \\ &a(u_h,v_h) = L(v_h),~ ~ \forall\,v_h \in V_h\end{aligned}\right.$$où $u_h$ est une solution approchée de la solution exacte $u$.
L’idée est de résoudre cette formulation à partir d’une base $(\varphi_1,\varphi_2,\dots,\varphi_{n_h})$ de $V_h$. En décomposant $u_h$ dans cette base, on a :
$$u_h = \sum\limits_{j=1}^{n_h} u_j\,\varphi_j$$Le principe de la méthode de Galerkin est alors d’utiliser les fonctions de base $\varphi_i$ comme fonctions test dans la formulation faible. Ainsi :
$$\begin{aligned}\forall i \in [\![1,n_h]\!]~ ,~ &a\left(\sum\limits_{j=1}^{n_h} u_j\,\varphi_j,\varphi_i\right) = L(\varphi_i) \\ & \sum\limits_{j=1}^{n_h} u_j\,a(\varphi_j,\varphi_i) = L(\varphi_i)\end{aligned} $$Soit :
$$\boxed{\forall i \in [\![1,n_h]\!]~ ,~ \sum\limits_{j=1}^{n_h} a(\varphi_i,\varphi_j)\,u_j = L(\varphi_i)}$$En définissant les vecteurs ${\bf U_h} = \left( u_i \right)_{1\leq i \leq n_h}$ et ${\bf B_h} = \big( L(\varphi_i) \big)_{1\leq i \leq n_h}$, et la matrice ${{\bf A_h} = \big( a(\varphi_i,\varphi_j) \big)_{1\leq i,j\leq n_h}}$, le problème peut se mettre sous la forme matricielle :
$$\boxed{{\bf A_h} \, {\bf U_h} = {\bf B_h}}$$En résolvant ce système linéaire, on obtient alors les composantes $u_i$ de la solution approchée $u_h$ dans la base $(\varphi_i)_{1\leq i\leq n_h}$.
$\rightarrow$ La méthode des éléments finis permet d’appliquer ce principe de résolution.
On appelle élément fini le triplet $(K,P_K,\Sigma_K)$, où :
De plus, il faut qu’il existe une unique fonction de $P_K$ prenant des valeurs données aux degrés de liberté de $\Sigma_K$ : cette dernière condition définit l’unisolvance de l’élément fini $(K,P_K,\Sigma_K)$.
Nous allons maintenant voir comment mettre en place la méthode basée sur ces éléments.
Intéressons nous tout d’abords aux éléments géométriques $K$ obtenus par discrétisation du domaine d’étude $\Omega$, ce qu’on appelle le maillage. Le domaine (polyédrique) est alors partitionné avec le type d’entités élémentaires choisi (triangle, quadrangle, tétraèdre, hexaèdre, etc…).
Le maillage doit être conforme, c’est-à-dire :
Nous pouvons illustrer cela sur un exemple :
| Maillage conforme | Maillage non conforme |
|---|---|
|
|
|
Dans la suite, nous utiliserons le logiciel libre de maillage par éléments finis Gmsh développé par C. Geuzaine et J-F. Remacle.
Pour pouvoir discrétiser un problème donné, nous devons tout d’abord définir ses aspects purement géométriques. C’est-à-dire créer toutes les entités élémentaires :
Sur notre cas d’étude, on obtient par exemple :
Après avoir défini des tailles caractéristiques d’éléments et/ou discrétiser les lignes, on peut mailler notre géométrie. Ci-dessous, deux exemples :
| Triangulaire | Quadrangulaire |
|---|---|
|
|
Pour information, je vous donne en exemple les scripts permettant de générer ces maillages. Par habitude, j’utilise deux fichiers :
data.geo : servant à définir les paramètres qui pourront être partagés avec le solveur,poisson.geo : contenant la création de la géométrie.En plus des fonctions basiques de création des entités élémentaires, Gmsh possède des fonctions plus évoluées : transformations (symétries, translations, rotations, dilatations), extrusions. De plus, il intègre depuis quelques temps un noyau supplémentaire basé sur le moteur CAO OpenCascade permettant de créer directement des objets 2D ou 3D, d’effectuer des opérations booléennes, ou de manipuler des splines.
Pour plus d’informations sur le logiciel, vous pouvez :
$\rightarrow$ consulter le site officiel et son wiki ;
$\rightarrow$ suivre un très bon tutoriel fait par B. Thierry (ENSEM, promo 2007).
Ci-après, vous trouverez des exemples de maillages utilisés dans le cadre de travaux de mon équipe de recherche au laboratoire GREEN.
| Machine à courant continu | Ventouse magnétique |
|---|---|
|
|
| Accouplements magnétiques (maillages structurés) | Machines synchrones à griffes (mix structuré / non-structuré) |
|---|---|
|
|
|
|
$~$
Les exemples ci-dessus sont issus des thèses de :
Pour plus d’informations, ou si vous envisagez la possibilité de faire une thèse, n’hésitez pas à me contacter, ainsi que mes collègues : les Professeurs N. Takorabet ou D. Netter.
On peut également utiliser des maillages beaucoup plus complexes, comme par exemple les deux ci-dessous qui me permettent de calculer des courants induits dans des têtes humaines (par exemple pour des travailleurs opérant à proximité de sources basse fréquence comme des pinces à souder) :
| John Doe | Jane Doe |
|---|---|
|
|
$\boldsymbol{\rightarrow}$ Voir même un corps humain complet (oui, j’ai un peu galéré…) :
| John en entier (avec tous les organes) |
|---|
|
Maintenant que nous savons créer le maillage (noté $\mathcal{M}_h$), nous allons nous intéresser à la définition des fonctions qui permettront de constituer la base de l’approximation discrète $V_h$ de $V$ (fonctions de base).
Par souci de simplicité, nous ne nous intéresserons, dans un premier temps, qu’aux éléments nodaux permettant d’approximer $\text{H}^1_0(\Omega) \left(= \text{H}_0({\bf grad},\Omega)\right).$ On les appelle éléments de Lagrange. Pour clarifier un peu plus les choses, nous ne considérerons tout d’abord que des éléments s’appuyant sur des fonctions de base linéaires (polynômes d’ordre 1) permettant une approximation linéaire par morceaux de la fonction inconnue. On les appelle les éléments $\mathbf{P_1}$ (ou $\mathbf{Q_1}$ pour les quadrangles).
En fait, nous allons définir les fonctions de base sur ce qu’on appelle un élément de référence dont la forme géométrique est fixe et prédéfinie.
En dimension 1, les éléments sont des segments et l’élément de référence est le segment $[-1,1]$.
En dimension 2, nous pouvons donner le triangle et le quadrangle de référence :
Leur représentation graphique est donnée ci-dessous.
Nous nous contenterons de donner le tétraèdre et l’hexaèdre de référence :
Et leur représentation graphique est donnée ci-dessous.
Sur ces élément de références, nous pouvons définir des fonctions linéaires valant 1 sur un sommet (un nœud) et 0 sur les autres. Dans ce cas, on fera une approximation linéaire par morceaux de notre fonction inconnue. Ces fonctions seront notées $p_i$ pour $i \in [\![1,n_K]\!]$.
Nous aurons donc autant de fonctions de base que de sommets $n_K$ de l’élément de référence.
Dans le cas du triangle de référence, ces fonctions sont données par la table ci-dessous :
| Nœud $i$ | Fonction de base $p_i (u,v)$ |
|---|---|
| 1 | $p_1 (u,v) = 1 - u - v$ |
| 2 | $p_2 (u,v) = u$ |
| 3 | $p_3 (u,v) = v$ |
Une représentation graphique est donnée par :
On peut également les tracer en surélevé :
| $p_1$ | $p_2$ | $p_3$ |
|---|---|---|
|
|
|
Dans le cas du triangle de référence, ces fonctions sont données par la table ci-dessous :
| Nœud $i$ | Fonction de base $p_i (u,v)$ |
|---|---|
| 1 | $\displaystyle p_1 (u,v) = \frac{(1-u)(1-v)}{4}$ |
| 2 | $\displaystyle p_2 (u,v) = \frac{(1+u)(1-v)}{4}$ |
| 3 | $\displaystyle p_3 (u,v) = \frac{(1+u)(1+v)}{4}$ |
| 4 | $\displaystyle p_4 (u,v) = \frac{(1-u)(1+v)}{4}$ |
Une représentation graphique est donnée par :
On peut également les tracer en surélevé :
| $p_1$ | $p_2$ | $p_3$ | $p_4$ |
|---|---|---|---|
|
|
|
|
Dans le cas du tétraèdre de référence, ces fonctions sont données par la table ci-dessous :
| Nœud $i$ | Fonction de base $p_i (u,v,w)$ |
|---|---|
| 1 | $p_1 (u,v,w) = 1 - u - v - w$ |
| 2 | $p_2 (u,v,w) = u$ |
| 3 | $p_3 (u,v,w) = v$ |
| 4 | $p_4 (u,v,w) = w$ |
Et une représentation graphique de $p_1$ est donnée ci-dessous :
Dans le cas de l’hexaèdre de référence, ces fonctions sont données par la table ci-dessous :
| Nœud $i$ | Fonction de base $p_i (u,v,w)$ |
|---|---|
| 1 | $\displaystyle p_1 (u,v,w) = \frac{(1-u)(1-v)(1-w)}{8}$ |
| 2 | $\displaystyle p_2 (u,v,w) = \frac{(1+u)(1-v)(1-w)}{8}$ |
| 3 | $\displaystyle p_3 (u,v,w) = \frac{(1+u)(1+v)(1-w)}{8}$ |
| 4 | $\displaystyle p_4 (u,v,w) = \frac{(1-u)(1+v)(1-w)}{8}$ |
| 5 | $\displaystyle p_5 (u,v,w) = \frac{(1-u)(1-v)(1+w)}{8}$ |
| 6 | $\displaystyle p_6 (u,v,w) = \frac{(1+u)(1-v)(1+w)}{8}$ |
| 7 | $\displaystyle p_7 (u,v,w) = \frac{(1+u)(1+v)(1+w)}{8}$ |
| 8 | $\displaystyle p_8 (u,v,w) = \frac{(1-u)(1+v)(1+w)}{8}$ |
Et une représentation graphique de $p_1$ est donnée ci-dessous :
Considérons maintenant un élément réel $K$ de $\mathcal{M}_h$, et intéressons nous à la transformation permettant de passer de l’élément de référence associé $K_r$ à cet élément, ainsi que sa réciproque.
Soit ${\bf x}$ un point quelconque de $K$ de coordonnées $(x,y,z)$. Une paramétrisation possible de ${\bf x}$ est :
$${\bf x} = \begin{pmatrix}x\\y\\z\end{pmatrix} = \begin{pmatrix}x(u,v,w)\\y(u,v,w)\\z(u,v,w)\end{pmatrix} ={\bf x}({\bf u}) $$où ${\bf u} = (u,v,w) \in K_r$ permet de faire correspondre tout point $(x,y,z)$ de $K$ avec un point $(u,v,w)$ de $K_r$ de façon biunivoque (bijective). La transformation inverse est alors notée :
$${\bf u} = {\bf u}({\bf x})$$Une représentation schématique est donnée ci-dessous :
Les formules de différentiation des fonctions composées nous donnent alors :
$$\begin{pmatrix}\partial_u\\\partial_v\\\partial_w\end{pmatrix} = \underbrace{\begin{pmatrix}\partial_u\,x & \partial_u\,y & \partial_u\,z\\ \partial_v\,x & \partial_v\,y & \partial_v\,z\\ \partial_w\,x & \partial_w\,y & \partial_w\,z\end{pmatrix}}_{\displaystyle{\bf J}} \, \begin{pmatrix}\partial_x\\\partial_y\\\partial_z\end{pmatrix} $$$$\begin{pmatrix}\partial_x\\\partial_y\\\partial_z\end{pmatrix} = \underbrace{\begin{pmatrix}\partial_x\,u & \partial_x\,v & \partial_x\,w\\ \partial_y\,u & \partial_y\,v & \partial_y\,w\\ \partial_z\,u & \partial_z\,v & \partial_z\,w\end{pmatrix}}_{\displaystyle{\bf J^{-1}}} \, \begin{pmatrix}\partial_u\\\partial_v\\\partial_w\end{pmatrix} $$La matrice ${\bf J}$ est appelée matrice jacobienne de transformation. Les relations ci-dessus peuvent s’écrire sous forme plus compacte :
$$\partial_{\bf u} = {\bf J}\cdot\partial_{\bf x},~ ~ \text{avec :}~ ~ {\bf J} = \partial_{\bf u}\,{\bf x}$$$$\partial_{\bf x} = {\bf J^{-1}}\cdot\partial_{\bf u},~ ~ \text{avec :}~ ~ {\bf J^{-1}} = \partial_{\bf x}\,{\bf u}$$On rappelle au passage la formule de l’inverse :
$${\bf J^{-1}} = \frac{1}{\det {\bf J}}\,{}^{\operatorname {t} }\text{com}\,{\bf J}$$En pratique, on ne la calcule pas de cette façon. En effet, il est intéressant de remarquer que les colonnes de ${\bf J^{-1}}$ sont les gradients selon les coordonnées réelles ${\bf x}$ des fonctions $u({\bf x}),$ $v({\bf x}),$ et $w({\bf x})$. On les note : ${\bf grad_x}\,u({\bf x}),$ ${\bf grad_x}\,v({\bf x})$ et ${\bf grad_x}\,w({\bf x})$.
Nous pouvons calculer ces termes à partir des lignes de la matrice Jacobienne ${\bf J}$ qui sont ${}^{\operatorname t}\partial_u\,{\bf x},$ ${}^{\operatorname t}\partial_v\,{\bf x}$ et ${}^{\operatorname t}\partial_w\,{\bf x},$ et qui sont plus faciles à calculer.
Rappelons les deux propriétés de la comatrice suivantes :
$$\text{com}\,{}^{\operatorname t}{\bf M} = {}^{\operatorname t}\text{com}\,{\bf M}$$$${\bf M}\,{\bf a}\wedge{\bf M}\,{\bf b} = \text{com}\,{\bf M} \, ({\bf a}\wedge{\bf b})$$Ainsi, on a :
$$\begin{aligned}{\bf grad_x}\,u({\bf x}) &= {\bf J^{-1}}\,{\bf e_x} \\ &= \frac{1}{\det\,{\bf J}}\,{}^{\operatorname t}\text{com}\,{\bf J}\,{\bf e_x} \\ &= \frac{1}{\det\,{\bf J}}\,\text{com}\,{}^{\operatorname t}{\bf J}\,{\bf e_x}\\ &= \frac{1}{\det\,{\bf J}}\,\text{com}\,{}^{\operatorname t}{\bf J}\,({\bf e_y}\wedge{\bf e_z})\\ &= \frac{1}{\det\,{\bf J}}\, {}^{\operatorname t}{\bf J}\,{\bf e_y}\wedge{}^{\operatorname t}{\bf J}\,{\bf e_z}\\ &= \frac{1}{\det\,{\bf J}}\,\left(\partial_v\,{\bf x}\wedge\partial_w\,{\bf x} \right)\end{aligned}$$On peut faire de même pour les deux autres colonnes, au final, on obtient :
$$\boxed{\left\{\begin{aligned}{\bf grad_x}\,u({\bf x}) = \frac{1}{\det\,{\bf J}}\,\left(\partial_v\,{\bf x}\wedge\partial_w\,{\bf x} \right)\\[2ex] {\bf grad_x}\,v({\bf x}) = \frac{1}{\det\,{\bf J}}\,\left(\partial_w\,{\bf x}\wedge\partial_u\,{\bf x} \right) \\[2ex] {\bf grad_x}\,w({\bf x}) = \frac{1}{\det\,{\bf J}}\,\left(\partial_u\,{\bf x}\wedge\partial_v\,{\bf x} \right)\end{aligned}\right.}$$Le déterminant de la matrice jacobienne de transformation ($\det\,{\bf J},$ appelé Jacobien) peut se calculer par :
$$\det\,{\bf J} = \det\,{}^{\operatorname t}{\bf J} = \partial_u\,{\bf x}\cdot(\partial_v\,{\bf x}\wedge\partial_w\,{\bf x})$$Ce Jacobien est très important car il intervient directement dans les formules d’intégration sur l’élément $K$ qui sont (par la formule de changement de variable) :
$$\iiint_K f({\bf x})\,\text{d}{\bf x} = \iiint_{K_r} f\left({\bf x}({\bf u})\right)\,\det\,{\bf J}\,\text{d}{\bf u}$$Ci-dessous, vous trouverez une petite application permettant de calculer le Jacobien d’un triangle quelconque (avec l’aimable autorisation de son auteur B. Thierry) :
Déplacez les sommets du triangle pour modifier la valeur du Jacobien. Quand il est négatif cela signifie que le triangle est "retourné" par rapport au triangle de référence
Un élément isoparamétrique est un élément fini dont les fonctions de base nodales (qui permettent l’interpolation des champs scalaires) sont aussi utilisées pour paramétrer l’élément géométrique associé (fonctions d’interpolation géométrique) [Dula94].
Les fonctions de base vues ci-dessus sont définies par morceaux, sur chacun des éléments géométriques qui remplissent le domaine étudié, et permettent alors d’assurer leur continuité aux interfaces entre éléments. Ainsi, il n’y aura pas de discontinuité ni des champs scalaires interpolés, ni des coordonnées après transformation des éléments de référence vers les éléments réels. De telles fonctions de base sont dites conformes.
Considérons un élément fini nodal $(K,P_K,\Sigma_K)$. Si $n_K$ est l’ensemble des nœuds de $K$, de coordonnées ${\bf x_i},~i \in \mathbb{N},$ et si les $p_i({\bf u}),~i \in [\![1,n_K]\!],$ sont ses fonctions de base exprimées dans les coordonnées ${\bf u}$ de l’élément de référence $K_r$ associé à $K$, alors la paramétrisation de $K$ $\left({\bf x} = {\bf x}({\bf u})\right)$ est donnée par :
$$\boxed{{\bf x} = \sum\limits_{i = 1}^{n_K} {\bf x_i}\,p_i({\bf u})}$$et cet élément est isoparamétrique.
Dans ce cas, les fonctions de base permettent de passer de l’élément de référence à un élément réel et ainsi de déterminer la forme de ce dernier. C’est pour cela qu’on les appelle aussi fonctions de forme.
À l’aide de la transformation géométrique précédente et des fonctions de bases $p_i,~i \in [\![1, n_K]\!]$ définies sur l’élément de référence, nous pouvons définir les fonctions de base globales $\varphi_i,~i \in [\![1, n_h]\!]$ définies sur le maillage complet $\mathcal{M}_h$. Il y en aura autant que de nœuds et chacune sera liée à un nœud particulier. Elle seront non-nulles uniquement sur les éléments contenant le nœud considéré. Cet ensemble d’éléments est appelé le support de la fonction de base.
Finalement, d’après ce qui précède, la fonction de base $\varphi_k$ liée au nœud $k \in [\![1, n_h]\!]$ sera nulle sur tous les nœuds du maillage sauf le nœud $k$ où elle prendra la valeur $1$. Soit, en notant ${\bf x_j}$ les coordonnées du nœud $j$ :
$$\boxed{\forall i,j \in [\![1,n_h]\!], ~ \varphi_i({\bf x_j}) = \delta_{ij}}$$Une illustration interactive est donnée par l’application suivante (toujours avec l’aimable autorisation de son auteur) :
Application : Cliquez sur un sommet pour faire apparaitre la fonction de forme $\mathbf{P_1}$ associée. Les triangles où la fonction n’est pas nulle forment le support de la fonction de forme.
La famille $(\varphi_i)_{0\le i \le n_h}$ ainsi construite forme une base de $V_h$ qui est donc de dimension $n_h$ égale aux nombre de nœuds du maillage $\mathcal{M}_h$.
Soit $u_h$ une fonction de $V_h$. En la décomposant dans la base obtenue :
$$ u_h = \sum\limits_{i = 1}^{n_h} u_i\,\varphi_i$$Les coefficients $u_i$ sont en fait les valeurs de la fonctions au nœuds. En effet en appliquant en ${\bf x_j}$ :
$$u_h({\bf x_j}) = \sum\limits_{i = 1}^{n_h} u_i\,\varphi_i({\bf x_j}) = u_j$$Finalement, une fonction $u_h$ de $V_h$ se décompose donc dans la base $(\varphi_i)_{0\le i \le n_h}$ selon :
$$\boxed{u_h = \sum\limits_{i = 1}^{n_h} u_h({\bf x_i})\,\varphi_i}$$et elle est entièrement caractérisée par ses $n_h$ valeurs aux nœuds.
Dans le cas présenté, les fonctions de base de l’élément de référence servent aussi à paramétrer les éléments réels (fonctions d’interpolation géométrique). Cela correspond au cas particulier des éléments isoparamétriques, mais ces deux types de fonctions sont indépendants l’un de l’autre :
Les éléments présentés ci-dessus conviennent parfaitement pour représenter des champs scalaires, nous verrons plus tard lesquels utiliser pour les champs vectoriels.
Plus que les fonctions de bases globales, ce sont les fonctions de base sur les éléments de référence et la matrice Jacobienne qui sont utiles pour créer le système à résoudre, c’est ce que nous allons voir de suite.
On rappelle l’approximation de la formulation faible à résoudre obtenue par la méthode de Galerkin :
$$\forall i \in [\![1,n_h]\!]~ ,~ \sum\limits_{j=1}^{n_h} a(\varphi_i,\varphi_j)\,u_j = L(\varphi_i)$$Avec :
$$ a(\varphi_i,\varphi_j) = \iiint_{\Omega} \alpha\,{\bf grad}\,\varphi_i \cdot {\bf grad}\,\varphi_j~\text{d}\Omega$$Et :
$$\displaystyle L(\varphi_i) = \underbrace{\iiint_{\Omega} \beta\,\varphi_i~\text{d}\Omega}_{\displaystyle\left(\beta,\varphi_i\right)_{\Omega}} + \underbrace{\iint_{\Gamma_n} \alpha\,\gamma\,\varphi_i~\text{d}\Gamma}_{\displaystyle\left<\alpha\,\gamma,\varphi_i\right>_{\Gamma_n}}$$Cette formulation est équivalente à la forme matricielle :
$${\bf A_h}\,{\bf U_h} = {\bf B_h}$$avec :
$$\left\{\begin{aligned}&{\bf A_h} = \big(a(\varphi_i,\varphi_j)\big)_{1\le i,j\le n_h} = (A_{ij})_{1\le i,j\le n_h} \\[1ex] &{\bf B_h} = \big(L(\varphi_i)\big)_{1\le i\le n_h} = (B_i)_{1\le i\le n_h} \end{aligned}\right.$$Nous pourrions calculer les $n_h^2$ termes $A_{ij}$ séparément en effectuant une double boucle sur les nœuds. Cette façon de procéder n’est pas très pertinente car de nombreux termes sont nuls : typiquement, lorsque le nœud $i$ n’appartient pas à un élément du support de $\varphi_j$.
Une méthode plus efficace consiste à boucler sur tous les éléments de $\mathcal{M}_h$ et pour chaque élément $K_q$ calculer sa contribution $a_q(\varphi_k,\varphi_l),$ $(k,l) \in [\![1,n_{K_q}]\!]$ à la matrice globale ${\bf A_h}$ :
$$a_q(\varphi_k,\varphi_l) = \iiint_{K_q} \alpha({\bf x})\,{\bf grad}\,\varphi_k({\bf x}) \cdot {\bf grad}\,\varphi_l({\bf x})~\text{d}{\bf x} $$D’après notre formule de changement de base :
$$a_q(\varphi_k,\varphi_l) = \iiint_{K_r} \alpha\left({\bf x}({\bf u})\right)\,{\bf grad}\,\varphi_k\left({\bf x}({\bf u})\right) \cdot {\bf grad}\,\varphi_l\left({\bf x}({\bf u})\right)\,\text{det}\,{\bf J_q}~\text{d}{\bf u} $$En utilisant les relations entre les gradients et la matrice jacobienne de transformation géométrique :
$$a_q(\varphi_k,\varphi_l) = \iiint_{K_r} \alpha\left({\bf x}({\bf u})\right)~ \left({\bf J_q^{-1}}\,{\bf grad}\,p_k({\bf u})\right) \cdot \left({\bf J_q^{-1}}\,{\bf grad}\,p_l({\bf u})\right)\,\text{det}\,{\bf J_q}~\text{d}{\bf u}$$Finalement :
$$\boxed{a_q(\varphi_k,\varphi_l) = \text{det}\,{\bf J_q} \iiint_{K_r} \alpha\left({\bf x}({\bf u})\right)~ {}^{\operatorname t}\left({\bf grad}\,p_k({\bf u})\right)\,{}^{\operatorname t}({\bf J_q^{-1}})\,{\bf J_q^{-1}}\,{\bf grad}\,p_l({\bf u})~ \text{d}{\bf u}}$$Nous avons vu comment calculer ${\bf x}({\bf u})$, $\text{det}\,{\bf J_q}$ et ${\bf J_q^{-1}}$ à la page précédente, et les gradients des fonctions de base se déterminent facilement à partir de leurs expressions.
Ainsi, pour calculer les $n_{K_q}$ contributions de l’élément $K_q$, il nous reste plus qu’à savoir comment déterminer l’intégrale sur l’élément de référence $K_r$. Nous le verrons ci-après.
Intéressons nous maintenant au membre de droite en l’absence de condition de Neumann particulière.
Nous allons procéder de la même manière que pour ${\bf A_h}$ en calculant les contributions $L_q(\varphi_k)$ de chaque élément $K_q$ à ${\bf B_h}$ :
$$\begin{aligned}L_q(\varphi_k) &= \iiint_{K_q} \beta({\bf x})\,\varphi_k({\bf x})\,\text{d}{\bf x} \\ &= \iiint_{K_r} \beta\left({\bf x}({\bf u})\right)\,p_k({\bf u})\,\text{det}\,{\bf J_q}~\text{d}{\bf u}\end{aligned}$$Finalement :
$$\boxed{L_q(\varphi_k) = \text{det}\,{\bf J_q} \iiint_{K_r} \beta\left({\bf x}({\bf u})\right)\,p_k({\bf u})~\text{d}{\bf u}}$$Dans le cas d’une condition de Neumann non nulle, nous procédons comme ci-dessus mais en plus nous devons également calculer les contributions $L_q^n(\varphi_k)$ des éléments de la frontière $\Gamma_n$. Ceux-ci sont d’une dimension moindre que les éléments constituant le domaine $\Omega$. Ce seront donc des triangles ou des quadrangles en 3D, ou des segments en 2D.
Pour chaque éléments $K_q^n$ de $\Gamma_n$, nous calculons ainsi :
$$L_q^n(\varphi_k) = \iint_{K_q^n} \alpha({\bf x})\,\gamma({\bf x})\,\varphi_k({\bf x})\,\text{d}{\bf x}$$$$\boxed{L_q^n(\varphi_k) = \text{det}\,{\bf J_q^n} \iint_{K_r^n} \alpha\left({\bf x}({\bf u})\right)\,\gamma\left({\bf x}({\bf u})\right)\,p_k({\bf u})\,\text{d}{\bf u}}$$où $K_r^n$ est l’élément de référence associé aux éléments du bord du domaine.
Pour achever le calcul des coefficients des deux matrices, nous devons calculer les intégrales des différents termes sur le (ou les) élément(s) de référence. Pour cela nous allons utiliser une méthode de quadrature ; et généralement, c’est une formule de quadrature de Gauss.
Ainsi, on aura :
$$\iiint_{K_r} f({\bf u})\,\text{d}{\bf u} \approxeq \sum\limits_{m=1}^{n_g} \omega_m\,f({\bf x_{g_m}})$$où les ${\bf x_{g_m}}$ sont les coordonnées des $n_g$ points de Gauss choisis et $w_m$ les poids associés.
On notera que le calcul donne même la valeur exacte dans le cas de polynôme (à condition d’avoir le bon nombre de points de Gauss).
Dans la table ci-dessous, sont donnés quelques exemples :
| Élément | $n_g$ | Points de Gauss | Poids | Exacte pour |
|---|---|---|---|---|
| Segment | $1$ | $x_g = 0$ | $\omega =2$ | degré 1 |
| $2$ | $x_{g_1} = -\frac{1}{\sqrt{3}}$, $x_{g_2} = \frac{1}{\sqrt{3}}$ | $\omega_1 = \omega_2 = 1$ | degré 2 | |
| Triangle | $1$ | ${\bf x_g} = (\frac{1}{3},\frac{1}{3})$ | $\omega=\frac{1}{2}$ | degré 1 |
| $3$ | ${\bf x_{g_1}} = (\frac{1}{6},\frac{1}{6})$, ${\bf x_{g_2}} = (\frac{2}{3},\frac{1}{6})$, ${\bf x_{g_3}} = (\frac{1}{6},\frac{2}{3})$ | $\omega_1 = \omega_2 = \omega_3 = \frac{1}{6}$ | degré 2 | |
| Quadrangle | $1$ | ${\bf x_g} = (0,0)$ | $\omega=4$ | degré 1 |
| $4$ | ${\bf x_{g_1}} = (\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}})$, ${\bf x_{g_2}} = (-\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}})$, ${\bf x_{g_3}} = (\frac{1}{\sqrt{3}},-\frac{1}{\sqrt{3}})$, ${\bf x_{g_4}} = (-\frac{1}{\sqrt{3}},-\frac{1}{\sqrt{3}})$ | $\omega_1 = \omega_2$ $= \omega_3$ $= \omega_4$ $= 1$ | degré 2 | |
| Tétraèdre | $1$ | ${\bf x_g} = (\frac{1}{4},\frac{1}{4},\frac{1}{4})$ | $\omega=\frac{1}{6}$ | degré 1 |
| $4$ | ${\bf x_{g_1}} = (\alpha,\alpha,\alpha)$, ${\bf x_{g_2}} = (\alpha,\alpha,\beta)$, ${\bf x_{g_3}} = (\alpha,\beta,\alpha)$, ${{\bf x_{g_4}} = (\beta,\alpha,\alpha)}$, avec : ${\alpha = 0,138196601125}$, ${\beta = 0,585410196625}$ | $\omega_1 = \omega_2$ $= \omega_3 = \omega_4$ $= \frac{1}{24}$ | degré 2 |
On remarquera que dans tous les cas, la somme des poids est égale à la longueur, l’aire ou le volume de l’élément considéré.
Contrairement à ce qu’on pourrait penser, le choix du nombre de points de Gauss n’est pas primordial car l’erreur commise par la méthode d’intégration est faible devant les autres approximations de la méthode. Les seconds choix de la table précédente sont donc tout à fait pertinents (on n’utilisera que des fonctions de bases linéaires, voire quadratiques).
Après génération du système linéaire précédent, il ne reste plus qu’à résoudre. Il convient donc de choisir une méthode / un algorithme pouvant résoudre efficacement des systèmes de grande taille mais creux.
Pour les différentes méthodes de résolution (avec ou sans pré-conditionnement), je vous renvoie à l’EC d’« Analyse numérique », du tronc commun de S6.
Dans le cadre de ce cours, nous utiliserons GetDP pour définir et résoudre nos problèmes. Ce dernier s’appuie sur la suite PETSc configurée pour appeler le solveur linéaire direct multithreadé MUMPS (avec factorisation LU).
Après résolution, nous avons accès aux valeurs des degrés de liberté (valeurs du champ scalaire recherché) en chaque nœud du maillage. Nous pouvons ensuite calculer ce champ ou ses dérivés à partir de ces valeurs pour les tracer ou calculer tout autre grandeur locale ou globale. C’est ce qu’on appelle le post-traitement (post-processing).
Nous allons mettre en application ce qui précède sur notre exemple fil rouge, en écrivant le script GetDP permettant de le résoudre. Ce dernier doit porter le même nom que le fichier *.geo avec l’extension *.pro, ici : poisson.pro. Le déroulement suit les différentes étapes décrites dans les sections ci-avant, soit :
J’ai regroupé l’ensemble des fichiers permettant de résoudre notre exemple dans un petit modèle Onelab pouvant être utilisé directement dans Gmsh.
Vous pouvez même l’exécuter sur votre smartphone avec l’application Onelab (disponible sur les stores officiels) !
Edit (janv. 2026) : Pour ceux qui sont sous Android ou dérivés, l’app n’est malheureusement plus dans le Play Store. Vous pouvez néanmoins télécharger l’apk à partir de ce lien et l’installer directement (je procède comme ça sous GrapheneOS).
J’ai rajouté :
Exemple fil rouge à télécharger (compatible PC et smartphone) :
Nous venons de développer complètement le principe de la méthode des éléments finis pour des champs scalaires appartenant à $\text{H}^1(\Omega) \left(= \text{H}({\bf grad},\Omega)\right)$. Pour cela, nous avons utilisé des éléments de Lagrange bien adaptés pour ce type d’inconnues.
Cependant les différents champs constituant le champ électromagnétique sont également des grandeurs vectorielles pouvant appartenir à ${\bf H}({\bf rot},\Omega)$ ou ${\bf H}(\text{div},\Omega)$. Typiquement nous ne pouvons résoudre notre formulation rot-rot avec des éléments de Lagrange tels quels. Nous pourrions choisir de résoudre avec trois champs scalaires pour déterminer notre vecteur inconnu (un par composante), mais cette solution est loin d’être la plus pertinente en terme d’efficacité et d’aspects géométriques (respect de certaines conditions de continuité aux interfaces).
Parmi la grande variété d’éléments finis existants, une famille est très largement utilisée pour modéliser des problèmes électromagnétiques : les éléments finis de Whitney. En fait, ces derniers ont été spécifiquement créés dans ce but [Boss93]. Basés sur les formes différentielles du même nom, ils offrent un cadre idéal pour approximer nos champs en 3D car ils permettent de construire de manière séquentielle les espaces discrets auxquels ils appartiennent. Nous allons donc présenter succinctement ces éléments.
Contrairement aux éléments précédents, les éléments de Whitney sont initialement définis en 3D et s’appuient sur des degrés de libertés pouvant être reliés aux valeurs aux nœuds (éléments nodaux), mais aussi aux circulations le long des arêtes (éléments d’arête), aux flux au travers des facettes (éléments de facette), et aux intégrales sur les volumes (éléments de volume) des éléments du maillage. Compte tenu de ces différences de nature entre degrés de liberté, on les appelle parfois éléments finis mixtes.
Nous verrons dans la partie applicative qu’il est tout a fait possible de traiter des cas 2D. Selon les cas : soit directement, soit en utilisant certaines adaptations (arêtes ou facettes fictives orthogonales au plan d’étude).
Le principe de base est de construire successivement les espaces discrets $W^0$, $W^1$, $W^2$ et $W^3$ permettant d’approximer respectivement $\text{H}({\bf grad},\Omega)$, ${\bf H}({\bf rot},\Omega)$, ${\bf H}(\text{div},\Omega)$ et $L^2(\Omega)$, à partir d’une séquence de fonctions de base, liées respectivement aux nœuds $(\mathcal{N_h})$, aux arêtes $(\mathcal{A_h})$, aux facettes $(\mathcal{F_h})$ et aux volumes $(\mathcal{V_h})$ des éléments du maillage $\mathcal{M_h}$.
On peut synthétiser l’ensemble dans la table suivante :
| Espace | Approxime | Fonctions de bases | Propriété | Degré de liberté | Type d’élément |
|---|---|---|---|---|---|
| $W^0$ | $\text{H}({\bf grad},\Omega)$ | $\{s_n,~ n\in\mathcal{N_h}\}$ | $s_i({\bf x_j}) = \delta_{ij},~$ ${\forall i,j \in\mathcal{N}_h}$ | Valeur nodale | Élément nodal |
| $W^1$ | ${\bf H}({\bf rot},\Omega)$ | $\{ {\bf s_a},~ a\in \mathcal{A_h}\}$ | $\displaystyle \int_j {\bf s_i}\cdot{\bf dl} = \delta_{ij},~$ ${\forall i,j \in\mathcal{A_h}}$ | Circulation le long d’une arête | Élément d’arête |
| $W^2$ | ${\bf H}(\text{div},\Omega)$ | $\{ {\bf s_f},~ f\in \mathcal{F_h}\}$ | $\displaystyle \iint_j {\bf s_i}\cdot{\bf n}~ \text{d}S =\delta_{ij},~$ ${\forall i,j \in\mathcal{F}_h}$ | Flux à travers une face | Élément de facette |
| $W^3$ | $L^2(\Omega)$ | $\{s_v,~ v\in\mathcal{V_h}\}$ | $\displaystyle \iiint_j s_i~\text{d}V = \delta_{ij},~$ ${\forall i,j \in\mathcal{V}_h}$ | Intégrale volumique | Élément de volume |
En procédant ainsi, on respecte la séquence vue sur les espaces continus, à savoir :
$$W^0 \xrightarrow[]{{\bf grad}} W^1 \xrightarrow[]{{\bf rot}} W^2 \xrightarrow[]{\text{div}} W^3 $$Et l’ensemble nous permet également d’assurer les continuités aux interfaces :
Voyons maintenant comment définir ces fonctions de base (sur les éléments de références).
Puisque $W_0$ est une approximation de $\text{H}({\bf grad},\Omega) = \text{H}^1$, nous pouvons réutiliser les fonctions définies dans les parties précédentes. Nous aurons donc pour tout nœud $n_i = \{i\}$ de $\mathcal{N_h}$ :
$$s_{n_i} = p_i,~ ~ \forall\,i \in [\![1, n_K]\!]$$Tout ce qui précède reste donc valable. Les fonctions de base nodales des éléments de Whitney sont celles des éléments de Lagrange.
À une arête ${\bf a_{ij}} = \{i\rightarrow j\}$ de $\mathcal{A_h}$, on associe :
$${\bf s_{a_{ij}}} = p_j~{\bf grad}\left(\sum\limits_{r\in \operatorname{N}(j,\overline{i})} p_r \right) - p_i~{\bf grad}\left(\sum\limits_{r\in \operatorname{N}(i,\overline{j})} p_r \right) $$Où $\operatorname{N}(i,\overline{j})$ représente l’ensemble des nœuds de la facette $f$ contenant le nœud $i$ mais pas le nœud $j$.
Remarque : Dans le cas du tétraèdre, l’expression ci-dessus se ramène à :
$${\bf s_{a_{ij}}} = p_i~{\bf grad}\,p_j - p_j~{\bf grad}\,p_i$$À titre d’illustration, nous traçons un exemple pour les triangle, quadrangle, tétraèdre et hexaèdre de référence ci-dessous :
| ${\bf s_{a_{23}}}$ sur tétraèdre | ${\bf s_{a_{23}}}$ sur hexaèdre |
|---|---|
|
|
En considérant une facette fictive $f$ orthogonale au plan d’étude, on peut définir nos fonctions de base en 2D, et dans ce cas, on obtient :
| ${\bf s_{a_{23}}}$ sur triangle | ${\bf s_{a_{23}}}$ sur quadrangle |
|---|---|
|
|
À une facette orientée ${\bf f_{ijk}} = \{i,j,k\} = \operatorname{N}({\bf f_{ijk}})$ (si tétrahèdre) ou ${\bf f_{ijkl}} = \{i,j,k,l\} = \operatorname{N}({\bf f_{ijkl}})$ (si hexaèdre) de $\mathcal{F_h}$, on associe :
$${\bf s_{f_{ijk(l)}}} = a~{\large \sum\limits_{q \in \operatorname{N}({\bf f_{ijk(l)}})}} { p_q \left({\bf grad}\,\sum\limits_{r \in \operatorname{N}(q,\overline{q+1})} p_r \right)\wedge\left({\bf grad}\,\sum\limits_{r \in \operatorname{N}(q,\overline{q-1})} p_r \right) }$$où la liste des nœuds de la facette $\operatorname{N}({\bf f})$ est rendue circulaire, $\operatorname{N}(i,\overline{j})$ conserve sa définition précédente et $a$ est un coefficient tel que : $ a= 2$ si $\dim(\operatorname{N}({\bf f})) = 3$ et $a = 1$ si $\dim(\operatorname{N}({\bf f})) = 4$.
On représente graphiquement ces fonctions de facette en 2D (facette fictive perpendiculaire au plan d’étude) dans la table ci-dessous :
| ${\bf s_{f_{23\perp}}}$ sur triangle | ${\bf s_{f_{23\perp}}}$ sur quadrangle |
|---|---|
|
|
À l’élément $K$, on fait correspondre la fonction de base de volume :
$$s_{v_K} = \frac{1}{\operatorname{vol}(K)}$$qui est donc constante sur l’élément.
Avec les définitions précédentes on peut décomposer tout champ scalaire ou vectoriel selon son espace d’appartenance :
Soit $u \in \text{H}({\bf grad},\Omega)$, son approximation sur $\mathcal{M_h}$ est $u_h \in W^0$ tel que :
$$u_h = \sum\limits_{n\in\mathcal{N_h}} \sigma_n(u_h)\,s_n = \sum\limits_{n\in\mathcal{N_h}} u_n\,s_n$$avec :
$$u_n = u_h({\bf x_n})$$Soit ${\bf u} \in \textbf{H}({\bf rot},\Omega)$, son approximation sur $\mathcal{M_h}$ est ${\bf u_h} \in W^1$ tel que :
$${\bf u_h} = \sum\limits_{a\in\mathcal{A_h}} \sigma_a({\bf u_h})\,{\bf s_a} = \sum\limits_{a\in\mathcal{A_h}} {u_a}\,{\bf s_a}$$avec :
$$u_a = \int_a {\bf u_h}\cdot{\bf dl}$$Soit ${\bf u} \in \textbf{H}(\text{div},\Omega)$, son approximation sur $\mathcal{M_h}$ est ${\bf u_h} \in W^2$ tel que :
$${\bf u_h} = \sum\limits_{f\in\mathcal{F_h}} \sigma_f({\bf u_h})\,{\bf s_f} = \sum\limits_{f\in\mathcal{F_h}} {u_f}\,{\bf s_f}$$avec :
$$u_f = \iint_f {\bf u_h}\cdot{\bf n}~\text{d}S$$Soit $u \in L^2(\Omega)$, son approximation sur $\mathcal{M_h}$ est $u_h \in W^3$ tel que:
$$u_h = \sum\limits_{v\in\mathcal{V_h}} \sigma_v(u_h)\,s_v = \sum\limits_{v\in\mathcal{V_h}} u_v\,s_v$$avec :
$$u_v = \iiint_v u_h~\text{d}V$$Ainsi, en décomposant nos approximations dans les bons espaces, on peut appliquer la méthode de Galerkin en utilisant les fonctions de base ci-dessus comme fonctions tests. Par analogie avec les éléments nodaux de la section précédente, on se ramènera donc à la résolution d’un système linaire où le vecteur inconnu ${\bf U_h}$ correspondra aux différents degrés de liberté.
On rappelle la formulation faible continue à résoudre :
$$\left\{\begin{aligned}&\text{Trouver}~ {\bf u} \in\textbf{H}_0({\bf rot},\Omega) ~ \text{tel que :} \\ & \left(\alpha\,{\bf rot\,u},{\bf rot\,v}\right)_{\Omega} = - \left(\boldsymbol{\beta},{\bf v}\right)_{\Omega} - \left<\alpha\,\boldsymbol{\gamma},{\bf v}\right>_{\Gamma_n},~ \forall\,{\bf v} \in \textbf{H}_0({\bf rot},\Omega) \end{aligned}\right.$$Qui devient en discret :
$$\left\{\begin{aligned}&\text{Trouver}~ {\bf u_h} \in W_0^1 ~ \text{tel que :} \\ & \left(\alpha\,{\bf rot\,u_h},{\bf rot\,v_h}\right)_{\Omega} + \left(\boldsymbol{\beta},{\bf v_h}\right)_{\Omega} + \left<\alpha\,\boldsymbol{\gamma},{\bf v_h}\right>_{\Gamma_n} = 0,~ \forall\,{\bf v_h} \in W_0^1\end{aligned}\right.$$Où $W_0^1 \in W^1$ est obtenu en annulant les degrés de liberté sur les arêtes appartenant à la frontière $\Gamma_d$. On constate que cela implique bien notre condition de bord :
$$\left.{\bf u_h}\wedge{\bf n}\right|_{\Gamma_d} = {\bf 0}~ , ~ \text{sur}~ \Gamma_d$$La façon de procéder est identique à celle basée sur les éléments nodaux. Nous pouvons calculer les contributions de chaque arête élément par élément et les ajouter à la matrice de raideur globale ${\bf A_h}$, ou au vecteur second membre ${\bf B_h}$. On peut finalement résoudre le système obtenu, mais seulement après avoir implanté une condition de jauge comme expliqué ci-après. En effet, une telle formulation sous forme brute conduit à une matrice ${\bf A_h}$ singulière.
Nous savons malheureusement qu’un problème général (en 3D) comme celui-ci n’a pas de solution unique (vient de l’équivalence avec la formulation forte) et que le champ ${\bf u_h}$ solution est défini à un champ de gradient près. Il est alors impossible de résoudre le système obtenu sans avoir ajouté de condition de jauge particulière. Nous allons donc voir comment procéder.
Dans le domaine continu, pour rendre ${\bf u}$ unique, on peut ajouter une condition du type :
$${\bf u}\cdot{\bf w} = 0$$Où ${\bf w}$ est un champ de vecteurs dont les lignes de champ ne sont pas fermées et sont telles qu’elles peuvent relier toute paire de points quelconque du domaine d’étude $\Omega$ (par exemple un champ de gradient divergeant).
Dans le domaine discret, nous allons construire un arbre. C’est un ensemble d’arêtes qui relient tous les nœuds du maillage entre eux sans former de chemin fermé. Les arêtes du maillage qui ne sont pas dans l’arbre constituent le co-arbre associé. Un moyen de jauger ${\bf u_h}$ est ainsi d’annuler les composantes sur les arrêtes de l’arbre (condition équivalente à celle continue ci-dessus). Afin de ne pas dénaturer les conditions imposées sur $\Gamma_d$, il convient de commencer à construire l’arbre à partir des arrêtes de cette frontière. Une représentation graphique est donnée ci-dessous :
En rajoutant la condition de jauge d’arbre dans l’espace fonctionnel auquel appartient ${\bf u_h}$ (que nous noterons $W_{0,\text{JA}}^1 \subset W_0^1$), nous pouvons alors résoudre notre formulation faible rot-rot qui s’exprime :
$$\boxed{\left\{\begin{aligned}&\text{Trouver}~ {\bf u_h} \in W_{0,\text{JA}}^1 ~ \text{tel que :} \\ & \left(\alpha\,{\bf rot\,u_h},{\bf rot\,v_h}\right)_{\Omega} + \left(\boldsymbol{\beta},{\bf v_h}\right)_{\Omega} + \left<\alpha\,\boldsymbol{\gamma},{\bf v_h}\right>_{\Gamma_n} = 0,~ \forall\,{\bf v_h} \in W_{0,\text{JA}}^1\end{aligned}\right.}$$Cette façon de procéder est la plus efficace d’un point de vue numérique car elle permet de réduire la taille du système à résoudre (par annulation des degrés de liberté liés au arêtes de l’arbre). Cependant la solution obtenue perd son sens physique, seul son rotationnel pourra être interprété physiquement. Pour de plus amples compléments, je vous renvoie à [Dula94] .
Une autre Jauge possible est d’annuler la divergence de notre inconnue (Jauge de Coulomb) :
$$\text{div}\,{\bf u} = 0$$Pour obtenir la forme faible de cette condition, considérons un champ scalaire $\mu$ quelconque de $\text{H}({\bf grad},\Omega)$. On a :
$$\iiint_{\Omega}\text{div}({\bf u})\,\mu ~\text{d}V = 0$$En utilisant la formule de Green en grad-div :
$$\iiint_{\Omega} {\bf u}\cdot{\bf grad}\,\mu ~ \text{d}V - \iiint_{\Omega} \text{div}(\mu\,{\bf u})~ \text{d}V =0$$En appliquant le théorème de la divergence et en séparant les deux frontières :
$$\iiint_{\Omega} {\bf u}\cdot{\bf grad}\,\mu ~ \text{d}V - \iint_{\Gamma_d} \mu\,{\bf u}\cdot{\bf n}~ \text{d}S - \iint_{\Gamma_n} \mu\,{\bf u}\cdot{\bf n} ~ \text{d}V =0$$Nous ne considérerons que le cas d’une condition de Neumann homogène sur $\Gamma_n$ $(\boldsymbol{\gamma} = 0)$, soit : $\left.{\bf rot\,u}\wedge{\bf n}\right|_{\Gamma_n} = 0$. On peut donc en déduire : $\left.{\bf u}\cdot{\bf n}\right|_{\Gamma_n} = 0$. En choisissant $\mu = 0$ sur $\Gamma_d$, on annule le second terme surfacique.
Ainsi, notre condition de Jauge est équivalente à :
$$\boxed{\iiint_{\Omega} {\bf u}\cdot{\bf grad}\,\mu ~ \text{d}V = 0~ ,~ ~ \forall \, \mu \in \text{H}_0({\bf grad},\Omega)}$$Si nous ajoutons uniquement ce terme à la formulation faible initiale, nous nous trouverons face à une formulation avec une fonction inconnue et deux fonctions tests. Pour remédier à cela, nous devons introduire un multiplicateur de Lagrange $\lambda$. Ainsi la formulation variationnelle complète devient :
$$\left\{\begin{aligned}&\text{Trouver}~ {\bf u} \in\textbf{H}_0({\bf rot},\Omega),~ \lambda \in \text{H}_0({\bf grad},\Omega) ~ \text{tels que :} \\ &\left(\alpha\,{\bf rot\,u},{\bf rot\,v}\right)_{\Omega} + \left({\bf grad}\,\lambda,{\bf v}\right)_{\Omega} + \left(\boldsymbol{\beta},{\bf v}\right)_{\Omega} = 0,~ \forall\,{\bf v} \in \textbf{H}_0({\bf rot},\Omega)\\ &\left({\bf u},{\bf grad}\,\mu\right)_{\Omega} = 0, \forall \, \mu \in \text{H}_0({\bf grad},\Omega)\end{aligned}\right.$$La transcription dans le domaine discret est immédiate :
$$\boxed{\left\{\begin{aligned}&\text{Trouver}~ {\bf u_h} \in W_0^1,~ \lambda_h \in W_0^0~ \text{tels que :} \\ &\left(\alpha\,{\bf rot\,u_h},{\bf rot\,v_h}\right)_{\Omega} + \left({\bf grad}\,\lambda_h,{\bf v_h}\right)_{\Omega} + \left(\boldsymbol{\beta},{\bf v_h}\right)_{\Omega} = 0,~ \forall\,{\bf v_h} \in W_0^1 \\ &\left({\bf u_h},{\bf grad}\,\mu_h\right)_{\Omega} = 0, \forall \, \mu_h \in W_0^0 \end{aligned}\right.}$$Cette façon de procéder est plus lourde à résoudre : le système est beaucoup plus grand que le précédent car on a rajouté autant de degrés de liberté que de nœuds du maillage. Cependant, elle a le mérite de fournir une solution « plus physique » et plus facilement interprétable.
Où on appliquera la MEF à nos problèmes d’électromagnétisme dans l’ARQS.
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).
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.
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 :
Les équations de Maxwell en statique suivantes :
$$\left\{\begin{aligned}\text{div}\,{\bf d} &= \rho_q\\ {\bf rot}\,{\bf e} &= {\bf 0} \end{aligned}\right.$$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 :
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$$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 :
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 :
$$\left\{\begin{aligned}\text{div}\left(\varepsilon\,{\bf grad}\,v\right) + \rho_q &= 0~ , &\text{dans}~ \Omega\\ \left.v\right|_{\Gamma_d} &= 0~, &\text{sur}~ \Gamma_d \\ \left.v\right|_{\Gamma_i} &= v_i~, &\text{sur}~ \Gamma_i \\ \left.{\bf grad}\, v \cdot {\bf n}\right|_{\Gamma_n} &= 0~ , &\text{sur}~ \Gamma_n \\ \left.{\bf grad}\, v \cdot {\bf n}\right|_{\Gamma_j} &= -\frac{\sigma_{q_j}}{\varepsilon}~ , &\text{sur}~ \Gamma_j\end{aligned}\right.$$En appliquant le principe vu dans le chapitre précédent, nous en déduisons la formulation faible du problème :
Trouver $v \in \text{H}_{0i}({\bf grad},\Omega) = \{ u \in \text{H}({\bf grad},\Omega) : u|_{\Gamma_d} = 0,~ u|_{\Gamma_i} = v_i \}$, tel que :
$$\forall v' \in \text{H}_{0i}({\bf grad},\Omega),~~ \left(\varepsilon\,{\bf grad}\,v,\,{\bf grad}\,v'\right)_{\Omega} - \left(\rho_q,v'\right)_{\Omega} + \left<\sigma_{q_j},v'\right>_{\Gamma_j} = 0$$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 :
$$\forall v' \in \text{H}_{01}({\bf grad},\Omega),~~ \left(\varepsilon\,{\bf grad}\,v,\,{\bf grad}\,v'\right)_{\Omega} = 0$$Dans GetDP, cela se traduit par :
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.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 en .
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 :
.pro (oui, ça relève de la magie !).Comme je ne suis pas méchant, je vous donne la solution pour le condensateur cylindrique en .
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 :
Les équations de Maxwell en statique suivantes :
$$\left\{\begin{aligned}\text{div}\,{\bf j} &= 0 \\ {\bf rot}\,{\bf e} &= {\bf 0} \end{aligned}\right.$$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 :
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 :
$$ \text{div}\left(\sigma\,{\bf grad}\,v \right) = 0$$Sur les faces d’entrée et de sortie du courant (${\Gamma_d = \Gamma_{di}\cup\Gamma_{dj}}$), deux cas seront possibles :
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 :
$$\left\{\begin{aligned}\text{div}\left(\sigma\,{\bf grad}\,v\right) &= 0~ , &\text{dans}~ \Omega_c \\ \left.v\right|_{\Gamma_{di}} &= v_i~ , &\text{sur}~ \Gamma_{di} \\ \left.{\bf grad}\, v \cdot {\bf n}\right|_{\Gamma_{dj}} &= \pm \frac{j_{n_j}}{\sigma}~ , &\text{sur}~ \Gamma_{d_j} \\ \left.{\bf grad}\, v \cdot {\bf n}\right|_{\Gamma_n} &= 0~ , &\text{sur}~ \Gamma_n \end{aligned}\right.$$La formulation faible correspondant à la forme forte ci-dessus est obtenue rapidement par analogie avec l’électrostatique :
Trouver $v \in \text{H}_{0}({\bf grad},\Omega) = \{ u \in \text{H}({\bf grad},\Omega) : u|_{\Gamma_{di}} = v_i \}$, tel que :
$$\forall v' \in \text{H}_{0}({\bf grad},\Omega),~~ \left(\sigma\,{\bf grad}\,v,{\bf grad}\,v'\right)_{\Omega} + \left< j_{n_j},v' \right>_{\Gamma_{dj}} = 0$$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 en .
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 :
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}$
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; };
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) : .
Proposer des améliorations possibles du design d’un tel busbar (en terme de quantité de matière, masse, volume, pertes…).
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 :
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 :
${\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 :
$$\left\{\begin{aligned} {\bf rot}\left(\frac{1}{\mu_0}{\bf rot}\,{\bf a}\right) - {\bf j_s} &= 0 ~ ,\text{ dans }\Omega_c\\ {\bf rot}\left(\frac{1}{\mu_{ra}\mu_0}\left({\bf rot}\,{\bf a}-{\bf b_r}\right)\right) &= 0~ ,\text{ dans }\Omega_a\\ {\bf rot}\left(\frac{1}{\mu}{\bf rot}\,{\bf a}\right) &= 0 ~ ,\text{ ailleurs }\end{aligned}\right.$$Et nos conditions aux limites peuvent se traduire par :
$$\left\{\begin{aligned} \left.{\bf a}\wedge{\bf n}\right|_{\Gamma_d} &= {\bf 0}\\ \left.{\bf rot}\,{\bf a}\wedge{\bf n}\right|_{\Gamma_n} &= {\bf 0} \end{aligned}\right.$$En appliquant le principe détaillé dans le deuxième chapitre sur la formulation en rot-rot, on peut déduire la formulation faible de notre problème :
Trouver ${\bf a} \in \textbf{H}_{0}({\bf rot},\Omega) = \{ {\bf a} \in \textbf{H}({\bf rot},\Omega) : {\bf a}\wedge{\bf n}|_{\Gamma_{d}} = 0\}$, tel que :
$$\forall {\bf a'} \in \textbf{H}_{0}({\bf rot},\Omega),~~ \left(\mu^{-1}\,{\bf rot}\,{\bf a}\,,\,{\bf rot}\,{\bf a'}\right)_{\Omega} + \left(\mu^{-1}\,{\bf b_r}\,,\,{\bf rot}\,{\bf a'}\right)_{\Omega_a} + \left( -{\bf j_s} \,,\, {\bf a'}\right)_{\Omega_c}= 0$$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).
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 :
$$\left\{\begin{aligned}\text{div}\left(-\mu\,{\bf grad}\,\phi\right) &= 0~, &\text{dans}~ \Omega\setminus\Omega_a\\ \text{div}\left(-\mu\,{\bf grad}\,\phi + {\bf b_r}\right) &= 0~, &\text{dans}~ \Omega_a \end{aligned}\right.$$Et nos conditions aux limites peuvent se traduire par :
$$\left\{\begin{aligned} \left.\phi\right|_{\Gamma_n} &= {\bf 0} \\ \left.{\bf grad}\,\phi\cdot{\bf n}\right|_{\Gamma_d} &= {\bf 0} \end{aligned}\right.$$En procédant de la même façon qu’en électrostatique ou électrocinétique, on déduit la formulation faible associée :
Trouver $\phi \in \text{H}_{0}({\bf grad},\Omega) = \{ \phi \in \text{H}({\bf grad},\Omega) : \phi|_{\Gamma_{n}} = 0\}$, tel que :
$$\forall \phi' \in \text{H}_{0}({\bf grad},\Omega),~~ \left(-\mu\,{\bf grad}\,\phi\,,\,{\bf grad}\,\phi'\right)_{\Omega} + \left({\bf b_r}\,,\,{\bf grad}\,\phi'\right)_{\Omega_a}= 0$$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 :
Traiter numériquement l’exercice du cylindre plongé dans un champ uniforme :
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 :
Simuler le problème correspondant et donner une solution.
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) :
Un exemple de résultat possible pourrait être par exemple :
Apprès avoir détaillé son principe de fonctionnement, résoudre le problème correspondant à un coupleur à aimants tel que représenté ci-dessous :
Pour gagner du temps, vous pouvez télécharger la géométrie en .
Il faudra utiliser des conditions de périodicité puisque nous travaillerons sur un pôle.
Je vous donne un exemple d’implantation de ce type de conditions ci-dessous :
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 :
$$ \displaystyle \Gamma = \frac{2\,p~L_z}{\mu_0\,e} \iint_{S_e} \frac{(x^2-y^2)\,b_x b_y + x y\,(b_y^2-b_x^2)}{\sqrt{x^2+y^2}}\,dxdy$$Cette formule sera ensuite facile à implanter dans le post-processeur de GetDP.
Cet exemple est tiré de cet excellent article des Techniques de l’Ingénieur. Oui, je fais de l’auto-promo 😉
Si vous le souhaitez, vous pouvez le télécharger avec l’abonnement de l’UL en cliquant ici.
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.
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 :
$$\boxed{\left\{\begin{aligned}&\text{Trouver}~ {\bf a_h} \in W_0^1,~ \xi_h \in W_0^0~ \text{tels que :} \\ &\left(\mu^{-1}\,{\bf rot\,a_h},{\bf rot\,a_h'}\right)_{\Omega} + \left(\mu^{-1}\,{\bf b_r},{\bf rot\,a_h'}\right)_{\Omega_a} + \left(-{\bf j_s},{\bf a_h'}\right)_{\Omega} + \left({\bf grad}\,\xi_h,{\bf a_h'}\right)_{\Omega} = 0,~ \forall\,{\bf a_h'} \in W_0^1 \\ &\left({\bf a_h},{\bf grad}\,\xi_h'\right)_{\Omega} = 0, \forall \, \xi'_h \in W_0^0 \end{aligned}\right.}$$où $W_0^1$ et $W_0^0$ approximent respectivement $\textbf{H}_{0}({\bf rot},\Omega)$ et $\text{H}_{0}({\bf grad},\Omega)$.
Dans GetDP, il faudra d’abord définir ces deux espaces fonctionnels :
Puis la formulation correspondante (exemple pour un problème sans aimants avec lois de comportement magnétiques linéaires) :
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.
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 :
$$\boxed{\left\{\begin{aligned}&\text{Trouver}~ {\bf a_h} \in W_{0,\text{JA}}^1 ~ \text{tel que :} \\ & \left(\mu^{-1}\,{\bf rot\,a_h},{\bf rot\,a_h'}\right)_{\Omega} + \left(\mu^{-1}\,{\bf b_r},{\bf rot\,a_h'}\right)_{\Omega_a} + \left(-{\bf j_s},{\bf a_h'}\right)_{\Omega} = 0,~ \forall\,{\bf a_h'} \in W_{0,\text{JA}}^1\end{aligned}\right.}$$La taille du système sera ainsi fortement réduite par rapport au cas précédent puisqu’on s’affranchit du calcul de la grandeur nodale $\xi$.
En pratique, dans GetDP cela se définit de la façon suivante :
Et la formulation sera strictement la même que dans les cas 2D précédents.
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) à .
Avec, vous pourrez :
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 :
Pour gagner du temps, je vous fournis encore la géométie 3D en . Un exemple de géométrie possible est représenté ci-dessous :
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 équations qui nous concerneront seront :
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 :
$${\bf rot}\,\left({\bf e}+\frac{\partial\,{\bf a}}{\partial t}\right) = {\bf 0}$$.
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 :
$$ {\bf e} = -{\bf grad}\,v - \frac{\partial\,{\bf a}}{\partial t}$$Et finalement, l’équation de Maxwell-Ampère dans $\Omega_c$ donne :
$${\bf rot}\,\left( \mu^{-1}{\bf rot\,a}\right) = \left\{\begin{aligned}-\sigma\left({\bf grad}\,v + \frac{\partial\,{\bf a}}{\partial t}\right) ~ &~ \text{dans }\Omega_c \\ {\bf 0}~ &~ \text{ailleurs}\end{aligned}\right.$$La deuxième relation permettant de résoudre est la conservation de la densité de courant dans $\Omega_c$ :
$$\text{div}\left( \sigma\left({\bf grad}\,v + \frac{\partial\,{\bf a}}{\partial t}\right)\right) = 0$$Les conditions aux limites sont analogues à celles vues en magnétostatique ou électrocinétique.
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)$ :
Trouver ${\bf a} \in \textbf{H}_{0}({\bf rot},\Omega) = \{ {\bf a} \in \textbf{H}({\bf rot},\Omega) : {\bf a}\wedge{\bf n}|_{\Gamma_{d}} = 0\}$ et $v \in \text{H}_{0}({\bf grad},\Omega) = \{ u \in \text{H}({\bf grad},\Omega) : u|_{\Gamma_{di}} = v_i\}$, tels que :
$$\left\{\begin{aligned}\left(\mu^{-1}\,{\bf rot}\,{\bf a}\,,\,{\bf rot}\,{\bf a'}\right)_{\Omega} + \left(\sigma\,\partial_t\,{\bf a}\,,\,{\bf a'}\right)_{\Omega_c} + \left( \sigma\,{\bf grad}\,v \,,\, {\bf a'}\right)_{\Omega_c} = 0,~ ~ \forall\,{\bf a'} \in \textbf{H}_{0}({\bf rot},\Omega) \\ (\sigma\,{\bf grad}\,v \,,\, {\bf grad}\,v')_{\Omega_c} + (\sigma\,\partial_t\,{\bf a} \,,\, {\bf grad}\,v')_{\Omega_c} = 0,~ ~\forall\, v' \in \text{H}_{0}({\bf grad},\Omega)\end{aligned}\right.$$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 »).
La formulation faible à résoudre est alors :
Trouver $\underline{\bf a} \in \textbf{H}_{0}({\bf rot},\Omega) = \{ \underline{\bf a} \in \textbf{H}({\bf rot},\Omega) : \underline{\bf a}\wedge{\bf n}|_{\Gamma_{d}} = 0\}$ et $\underline{v} \in \text{H}_{0}({\bf grad},\Omega) = \{ \underline{u} \in \text{H}({\bf grad},\Omega) : \underline{u}|_{\Gamma_{di}} = \underline{v_i}\}$, tels que :
$$\left\{\begin{aligned}\left(\mu^{-1}\,{\bf rot}\,\underline{\bf a}\,,\,{\bf rot}\,\underline{\bf a'}\right)_{\Omega} + \left(\sigma\,j\omega\,\underline{\bf a}\,,\,\underline{\bf a'}\right)_{\Omega_c} + \left( \sigma\,{\bf grad}\,\underline{v} \,,\, \underline{\bf a'}\right)_{\Omega_c} = 0,~ ~ \forall\,\underline{\bf a'} \in \textbf{H}_{0}({\bf rot},\Omega) \\ (\sigma\,{\bf grad}\,\underline{v} \,,\, {\bf grad}\,\underline{v'})_{\Omega_c} + (\sigma\,j\omega\,\underline{\bf a} \,,\, {\bf grad}\,\underline{v'})_{\Omega_c} = 0,~ ~ \forall\, \underline{v'} \in \text{H}_{0}({\bf grad},\Omega)\end{aligned}\right.$$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 » :
Et la multiplication par $j\omega$ dans la formulation peut-être faite (au choix) :
Complex[0,1]*2*Pi*freq dans les expressions ;DtDof qui permet de définir une dérivée temporelle des degrés de liberté.À titre d’exemple, je vous propose de résoudre numériquement l’exercice sur la barre cylindre alimentée en alternatif.
Le modèle est .
Observer l’effet de peau ainsi que l’évolution de la résistance du conducteur en fonction de la fréquence :
Modifier les programmes précédents afin de modéliser une ligne triphasée et observer l’effet de proximité à 50 Hz :
Ça arrive bientôt…
Merci d’avoir choisi et suivi cet EC. J’espère que vous avez trouvé cela intéressant et qu’il vous a permis de mieux comprendre les aspects théoriques liés à la MEF appliquée à notre domaine, ainsi que ce qui est caché derrière les boîtes à clics interfaces des logiciels propriétaires.
[Anto21] X. Antoine, “Équations aux Dérivées Partielles, Théorie et résolution numérique”, Polycopié de cours de l’EC « EDP », ENSEM 2A Formation Énergie, 2021.
[Bast14] J. P. Bastos, N. Sadowski, “Magnetic Materials and 3D Finite Element Modeling”, CRC Press, Boca Raton, 2014. Disponible à la BU
[Boss93] A. Bossavit, “Électromagnétisme, en vue de la modélisation”, Mathématiques et Applications, Vol. 14, Springer-Verlag, Berlin, 1993. Disponible à la BU
[Dula94] P. Dular, “Modélisation du champ magnétique et des courants induits dans des systèmes tridimensionnels non linéaires”, Thèse de doctorat en Sciences appliquées, Université de Liège, 1994. Accessible en ligne via ORBi
[Four85] G. Fournet, “Électromagnétisme : à partir des équations locales”, Masson, Paris, 1985. Disponible à la BU
[Four93] G. Fournet, “Électromagnétisme”, Techniques de l’Ingénieur, traité Génie électrique, dossier D1020, 1993. Accessible en ligne via la BU
[Geuz01] C. Geuzaine, “High order hybrid finite element schemes for Maxwell’s equations taking thin structures and global quantities into account”, Thèse de doctorat en Sciences appliquées, Université de Liège, 2001. Accessible en ligne via ORBi
[Kern04] M. Kern, “Introduction à la méthode des éléments finis”, Polycopié de cours, École Nationale Supérieure des Mines de Paris, 2004. Accessible en ligne
[Land84] L. D. Landau, E. M. Lifshitz, L. P. Pitaevskii, “Electrodynamics of Continuous Media”, Course of Theoretical Physics, Vol. 8 (2nd ed.), Butterworth-Heinemann, Oxford, 1984. Disponible à la BU
[Nett18] D. Netter, “Le Blog de l’Énergie Électrique”, collection de notes de cours, ENSEM, version du 12 sept. 2018. Accessible en ligne via ARCHE
[Noga05] B. Nogarède, “Électrodynamique appliquée : fondements et principes physiques de l’électrotechnique”, Sciences sup. Dunod, Paris, 2005. Disponible à la BU
[Piri23] F. Piriou, and S. Clénet, “Modélisation Numérique en électromagnétisme Basse Fréquence Par la Méthode des éléments Finis”, 1st ed., ISTE Editions Ltd, London, 2023. Accessible en ligne via la BU
[Spit02] P. Spiteri, “Approche variationnelle pour la méthode des éléments finis”, Techniques de l’Ingénieur, traité Sciences fondamentales, Ref. AF503, 2002. Accessible en ligne via la BU
[Thie20] B. Thierry, “Main5 : Maillage et Éléments Finis”, Notes de cours, Polytech Sorbonne, 2020.