Chapitre 3 - Méthodes des éléments finis en électromagnétisme basse fréquence


Où on appliquera la MEF à nos problèmes d’électromagnétisme dans l’ARQS.

Sous-sections de MEF en Électromag. BF

Avant-propos

Maxwell House

Les éléments finis de Whitney développés à la fin du chapitre 2 nous fournissent toute une structure permettant d’accueillir parfaitement toutes les constituantes du champ électromagnétique en basse fréquence. Une figure valant mieux qu’un long discours, vous trouverez ci-dessous une version de la Maison de Maxwell (diagramme de Tonti 3D appliqué à l’électromagnétisme dans l’ARQS).

La Maison de Maxwell (Maxwell House) dans l'ARQS : les différents étages de chaque côté correspondent à un espace fonctionnel (en couleur), on passe d'un étage à l'autre par un opérateur différentiel (flèches noires), de l'avant à l'arrière par dérivation temporelle (flèches noires aussi), et d'un côté à l'autre par les lois de comportement (flèches grises).

Pour vous, outre l’aspect esthétique, le principal intérêt de ce diagramme est de rapidement trouver à quel espace appartient chaque champ scalaire ou vectoriel.

$\rightsquigarrow$ Dans la suite, nous allons voir comment formuler des problèmes simples du domaine avec la méthode développée précédemment.

Électrostatique

Dans le cas de l’électrostatique (en anglais electrostatics), nous nous intéresserons au calcul de la distribution du champ électrique ${\bf e}$ due à des charges statiques ou à différents niveaux de potentiel électrique.

Les relations intervenant en électrostatique sont :

  • 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 :

    • soit, d’après la condition sur la composante tangentielle de ${\bf e}$ : $$\left.{\bf n}\wedge{\bf e}\right|_{\Gamma_d} = {\bf 0}$$
    • soit, d’après la condition sur la composante normale de ${\bf d}$ : $$\left.{\bf n}\cdot{\bf d}\right|_{\Gamma_n} = 0$$

Formulation forte

Maxwell-Faraday en statique nous permet de définir le potentiel scalaire électrique ${\bf v}$ tel que : ${{\bf e} = -{\bf grad}\,v}$

Ainsi, dans le domaine d’étude $\Omega$, l’équation satisfaite par $v$ est obtenue par substitution dans l’équation de Maxwell-Gauss :

$$ \text{div}\left(\varepsilon\,{\bf grad}\,v \right) + \rho_q = 0$$
Important !

Attention : En électrostatique, les conducteurs sont supposés « parfaits », le champ électrique est nul à l’intérieur et donc $v$ y est constant.

En pratique, les conducteurs seront donc exclus du domaine d’étude et seules les conditions sur leurs bords $(\Gamma_k)$ interviendront :

  • soit en tant que condition de Dirichlet où : $$v|_{\Gamma_i} = v_i$$
  • soit en tant que condition de Neumann où : $$\left.\frac{\partial\,v}{\partial {\bf n}}\right|_{\Gamma_j} = {\bf grad}\,v\cdot {\bf n}\big|_{\Gamma_j} = -\frac{\sigma_{q_j}}{\varepsilon}$$

Aux frontières, nous aurons les mêmes conditions mais homogènes : $ v|_{\Gamma_d} = 0~$ ou $~{\bf grad}\,v\cdot {\bf n}\big|_{\Gamma_n} = 0$.

Finalement, la forme complète de la formulation forte à résoudre est donc :

$$\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.$$

Formulation faible

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$$

Exemples et applications

Condensateur cylindrique

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 :

Géométrie et maillage du condensateur cylindrique Géométrie et maillage du condensateur cylindrique

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 :

Constraint {	
  {Name TensionsImposees ;
    Case {      
      { Region Gamma1 ; Type Assign ; Value V1 ; }
      { Region Masse ; Type Assign ; Value 0 ; }
      { Region Dirichlet ; Type Assign ; Value 0 ; }
    }
  }
}

FunctionSpace {
  { Name Hgrad ; Type Form0 ; 
    BasisFunction { 
      { Name sn ; NameOfCoef vn ; Function BF_Node ; 
			Support Region[{Domaine}]  ; Entity NodesOf[All] ; }
    } 
    Constraint {
      { NameOfCoef vn ; EntityType NodesOf ; 
			NameOfConstraint TensionsImposees ; }
    }
  }
}

Formulation {
  { Name Electrostat ; Type FemEquation ;
    Quantity {
      { Name v  ; Type Local  ; NameOfSpace Hgrad ; }
    }
    Equation {
      Galerkin { [ eps[]*Dof{Grad v}  , {Grad v} ] ;
		In Domaine ; Jacobian Jvol ; Integration Integ ; }
    }
  }
}

Explications :

  1. On définit tout d’abord les contraintes sur le potentiel $v$.
  2. Celles-ci permettent de définir ensuite l’espace fonctionnel dans lequel se trouve $v$. On l’appelle ici Hgrad (correspond au $W^0$ du chapitre précédent) qui approxime $\text{H}_{01}({\bf grad},\Omega)$. On utilise des éléments nodaux avec les fonctions de base $s_n$ vues précédemment.
  3. Enfin, nous définissons la forme faible du problème.

Après résolution, on obtient la distribution du potentiel scalaire électrique : v dans le domaine v dans le domaine

Et on peut calculer la valeur de la capacité à partir de l’énergie électrostatique.

Vous pouvez télécharger les fichiers complets en cliquant ici .


Condensateurs plans

  1. Adapter l’exemple précédent pour résoudre le problème similaire du condensateur plan à armatures rectangulaires :
    • modifier la géométrie ;
    • changer de jacobien ;
    • adapter les formules de calcul de l’énergie et de capacité.
  2. Étudier également le cas d’un condensateur plan mais, cette fois, à armatures circulaires.

Potentiel Flottant

Dans des problèmes tels que définis ci-dessus, nous n’avons pas directement accès à la charge portée par les armatures (il faut nécessairement passer par des calculs post-processeur). Nous ne pouvons pas non plus proprement considérer cette charge comme une entrée du problème pour calculer la différence de potentiels associée.
Nous pourrions penser à utiliser une condition de Neumann non homogène sur l’armature portant le potentiel $V_1$ mais la valeur du potentiel résultant du calcul ne serait plus constante le long de l’armature.

Pour pallier cet inconvénient, les auteurs de GetDP ont implanté une fonctionalité très pratique : les Global quantities (grandeurs globales) qui permettent de définir ce qu’on appelle un potentiel flottant (floating potential) via un terme global (global term) dans la formulation.

Plutôt que de faire une explication un peu trop appoximative, je préfère vous renvoyer à l’exemple du wiki de ONELAB détaillant son principe de fonctionnement : Electrostatics with floating potentials :

  1. Télécharger les fichiers et lire attentivement les commentaires du .pro (oui, ça relève de la magie !).
  2. Adapter les exemples précédents pour les résoudre en potentiel flottant.

Comme je ne suis pas méchant, je vous donne la solution pour le condensateur cylindrique en cliquant ici .

Électrocinétique

Dans le cas de l’électrocinétique (en anglais electrokinetics ou current flow), nous nous intéresserons au calcul de la distribution de la densité volumique statique de courrant électrique ${\bf j}$ dans des conducteurs.

Les relations intervenant en électrocinétique sont :

  • 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 :

    • soit, d’après la condition sur la composante tangentielle de ${\bf e}$ sur les faces d’entrée ou sortie du courant $(\Gamma_d)$ : $$\left.{\bf n}\wedge{\bf e}\right|_{\Gamma_d} = {\bf 0}$$
    • soit, d’après la condition sur la composante normale de ${\bf j}$ sur les autres bords du conducteurs $(\Gamma_n)$ : $$\left.{\bf n}\cdot{\bf j}\right|_{\Gamma_n} = 0$$

Formulation forte

Comme en électrostatique, nous pouvons définir le potentiel scalaire électrique $v$ tel que : ${{\bf e} = -{\bf grad}\,v}$.

En réinjectant cette relation couplée à la loi de comportement dans la loi des nœuds locale, on obtient :

$$ \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 :

  1. Soit une condition de potentiel imposé par une condition de Dirichlet sur $\Gamma_{di}$ de type : $$v|_{\Gamma_{di}} = v_i$$
  2. Soit une condition de Neumann non-homogène imposant la densité de courant normale à la surface $\Gamma_{dj}$ : $$\left.\frac{\partial\,v}{\partial {\bf n}}\right|_{\Gamma_{dj}} = {\bf grad}\,v\cdot {\bf n}\big|_{\Gamma_{dj}} = \pm \frac{j_{n_j}}{\sigma}$$

Sur les autres bords du domaine $(\Gamma_n)$, nous aurons des conditions de Neumann homogènes : ${{\bf grad}\,v\cdot {\bf n}\big|_{\Gamma_{n}} = 0}$

Finalement, la forme complète de la formulation forte à résoudre est donc :

$$\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.$$

Formulation faible

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$$

Exemples et applications

Résistance de sole de four électrique

On considère une résistance de sole de four électrique telle que représentée ci-dessous : Résistance de sole Résistance de sole

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 :

Constraint {	
  {Name Dirichlet ;
    Case {      
      {Region Sentree ; Type Assign ; Value Veff ;}
      {Region Ssortie ; Type Assign ; Value 0. ;}
    }
  }
}

FunctionSpace {
  { Name Hgrad_v ; Type Form0 ; 
    BasisFunction { 
      { Name sn ; NameOfCoef vn ; Function BF_Node ; 
		Support Region[{Domaine}]  ; Entity NodesOf[All] ; }
    } 
    Constraint {
      { NameOfCoef vn  ; EntityType NodesOf ; NameOfConstraint Dirichlet ; }
    }
  }
}


Formulation {
  { Name Electrocinetique ; Type FemEquation ;
    Quantity {
      { Name v  ; Type Local  ; NameOfSpace Hgrad_v ; }
    } 
    Equation {
      Galerkin { [ Sigma[]*Dof{Grad v}  , {Grad v} ] ;
		In Domaine ; Jacobian Jvol ; Integration Integ ; }
     }
  }
}

Après calcul, on obtient pour la répartition du potentiel scalaire électrique : v dans résistance de sole v dans résistance de sole

Vous pouvez télécharger les fichiers complets en cliquant ici .

Pour la petite histoire, voici une photo de celle de mon four perso après démontage pour changement suite à un défaut d’isolement : Ma resistance de sole perso Ma resistance de sole perso


Résistance de charge

  1. Traiter le problème analogue sur une forme correspondant aux résistances de charge utilisées sur nos bancs expérimentaux : “Résistances de charges pour manip” “Résistances de charges pour manip”

  2. En dimensionner une permettant de dissiper 5 kW sous 48 V continus.
    On pourra utiliser comme conducteur du Constantan, par exemple du 45Ni-55Cu de résistivité électrique : $\rho = 14,9\cdot 10^{-7}~\Omega\cdot\text{m}$

Astuce

Pour la géométrie, on pensera à utiliser une extrusion combinant à la fois une translation selon le vecteur $(v_{t_x},v_{t_y},v_{t_z})$ et une rotation d’un angle $\alpha$ autour de l’axe de vecteur directeur $(v_{a_x},v_{a_y},v_{a_z})$ et passant par le point $(x_a,y_a,z_a)$. Dans Gmsh, cela se définit de façon la plus simple :
Extrude{ {vtx,vty,vtz}, {vax,vay,vaz}, {xa,ya,za}, alpha } { Surface{numero}; }
ou de façon un peu plus élaborée (si on veut structurer le maillage) :
sortie[] = Extrude{ {vtx,vty,vtz}, {vax,vay,vaz}, {xa,ya,za}, alpha } { Surface{numero}; Layers{Ncouches}; Recombine; };


Busbar avec potentiel flottant

Tout comme dans le cas électrostatique de la section précédente, il est possible de définir, dans l’espace fonctionnel d’approximation $W^0$ , des grandeurs globales associées à notre inconnue $v$. Dans le cas présent, elles correspondront aux potentiels globaux $V_i$ imposés sur certaines faces du bords du domaine, ou aux courants $I_j$ traversant certaines autres.

  1. 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. “Busbar” “Busbar” Pour éviter de perdre trop de temps, je vous fournis la géométrie (libre à vous de l’utiliser ou non) : Cliquer ici pour la télécharger .

  2. Proposer des améliorations possibles du design d’un tel busbar (en terme de quantité de matière, masse, volume, pertes…).

Magnétostatique

La magnétostatique (en anglais magnetostatics) correspond au calcul de la distribution du champ ou de l’induction magnétiques produit par des courants continus ou des aimants.

Considérons le domaine $\Omega$ de $\mathbb{R}^3$ constitué de plusieurs sous-domaines schématisé par la figure : Domaine magnétostatique Domaine magnétostatique

  • $\Omega_m$ correspondant à des matériaux magnétiques ;
  • $\Omega_c$ correspondant à des conducteurs électriques parcourus par une densité de courant imposée par une source extérieure ${\bf j_s}$ ;
  • $\Omega_a$ correspondant à des aimants permanents d’induction rémanente ${\bf b_r}$ ;
  • une zone d’air entourant ses différents sous-domaines.

Et sa frontière $\Gamma (=\partial\Omega)$ est divisée en deux morceaux $\Gamma_d$ et $\Gamma_n$.

Les relations régissant la magnétostatique sont :

  • Les équations de Maxwell en statique suivantes : $$\left\{\begin{aligned}\text{div}\,{\bf b} &= 0 \\ {\bf rot}\,{\bf h} &= {\bf j_s}~\text{dans } \Omega_c, \text{ et } {\bf 0} \text{ ailleurs}\end{aligned}\right.$$
  • Les lois de comportements magnétiques : $$\left\{\begin{aligned}{\bf b} &= \mu_0\,{\bf h}&\text{dans l'air et }\Omega_c \\ {\bf b} &= \mu\,{\bf h} &\text{dans~}\Omega_m \\ {\bf b} &= \mu_0\mu_{ra}\, {\bf h}+{\bf b_r} &\text{dans~}\Omega_a\end{aligned}\right.$$
  • Les contions aux limites, par exemple :
    • continuité de la composante normale de l’induction : ${\bf n}\cdot{\bf b} = 0$ sur $\Gamma_d$ ;
    • continuité de la composante tangentielle du champ : ${\bf n}\wedge{\bf h} = {\bf 0}$ sur $\Gamma_n$.

Formulation en potentiel vecteur magnétique

Formulation forte

${\bf b}$ étant à divergence nulle, nous avons vu dans le premier chapitre, qu’il était possible de définir un potentiel vecteur magnétique ${\bf a}$ tel que :

$${\bf b} = {\bf rot}\,{\bf a}$$

En lui associant les lois de comportement et en reportant dans l’équation de Maxwell-Ampère, on obtient :

$$\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.$$

Formulation faible

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$$
Important !

Rappelons ici que la formulation ci-dessus ne conduit pas à une solution unique comme vu dans les chapitres précédents, il nous faudra rajouter une condition de jauge sur ${\bf a}$ pour pouvoir le calculer en 3D (jauge de Coulomb ou jauge d’arbre).


Formulation en potentiel scalaire magnétique

Formulation forte

Dans le cas particulier d’un problème sans densité de courant source $({\bf j_s} = {\bf 0})$, on a : ${\bf rot}\,{\bf h} = {\bf 0}$ dans tout le domaine.

On peut alors définir un potentiel scalaire magnétique $\phi$ tel que :

$${\bf h} = -{\bf grad}\,\phi$$

En lui associant les lois de comportement et en reportant dans l’équation de Maxwell-Thomson, on obtient :

$$\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.$$

Formulation faible

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$$

Sous-sections de Magnétostatique

Problèmes 2D

Simplifications

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 :

FunctionSpace {
  { Name Hrot ; Type Form1P ; 
    BasisFunction { 
      { Name se ; NameOfCoef ae ; Function BF_PerpendicularEdge ; 
	 	  Support Region[{Domaine}]  ; Entity NodesOf[All] ; }  } 
    Constraint {
      { NameOfCoef ae  ; EntityType NodesOf ; NameOfConstraint Dirichlet ; }
    }
  }  
}

Applications

Un premier problème simple sans source : un cylindre plongé dans un champ uniforme

Traiter numériquement l’exercice du cylindre plongé dans un champ uniforme :

  1. Trouver les conditions aux limites à imposer aux frontières du domaine (cubique) permettant d’obtenir un champ horizontal uniforme de $1~\text{T}$ à l’intérieur.
  2. Ajouter un cylindre au centre du domaine pour traiter les 3 cas possibles en fonction du type de matériau : cylindre_dans_champ_uniforme cylindre_dans_champ_uniforme

Dimensionnement d’une inductance pour électronique de puissance

Pour des besoins d’électronique de puissance (ne m’en demandez pas plus), on désire concevoir une inductance de 200 µH sous 15 A avec :

  • des noyaux de ferrites de type E 42/21/20 dont les caractéristiques sont fournies par cette datasheet ;
  • du fil de cuivre émaillé de 2,5 mm² (diamètre de 1,8 mm)

Simuler le problème correspondant et donner une solution.

Information

Après avoir traité le problème avec une loi de comportement magnétique linéaire dans la ferrite ($\mu_r$ constant), il faudra sûrement résoudre le problème en prenant en compte la non-linéarité de cette loi. Pour ce faire, je vous renvoie à cette explication très bien faite sur le wiki de GetDP. Vous pouvez également consulter cette page.

Je vous fournis également la définition complète du matériau (N27) de la ferrite (en utilisant les mêmes notations que dans les exemples du wiki) :

Function{		
  If(NonLineaire)
    Mat_core_b = {0.000000, 0.068851, 0.079231, 0.096630, 0.104285, 0.117920, 0.123374, 
      0.150912, 0.191749, 0.202635, 0.239024, 0.243936, 0.245556, 0.267244, 
      0.290565, 0.296190, 0.312896, 0.313356, 0.335942, 0.364839, 0.368482, 
      0.376766, 0.477117, 0.494881, 0.666990, 0.828539, 0.832720};

    Mat_core_h = {0.000000, 44.635198, 47.996938, 54.360481, 57.254136, 63.513926, 66.205407, 
      83.184012, 128.638642, 150.321852, 280.295679, 307.761975, 317.501481, 
      486.037901, 775.608272, 869.477778, 1359.072840, 1373.269136, 2431.451852, 
      5699.176543, 6389.154321, 8458.007407, 61532.617284, 73419.641975, 
      197003.703704, 321027.901235, 324350.000000};

    Mat_core_b2 = Mat_core_b()^2;
    Mat_core_h2 = Mat_core_h()^2;
    Mat_core_nu = Mat_core_h() / Mat_core_b();
    Mat_core_nu(0) = Mat_core_nu(1);
    Mat_core_nu_b2  = ListAlt[Mat_core_b2(), Mat_core_nu()] ;
    nu_core[] = InterpolationLinear[ SquNorm[$1] ]{ Mat_core_nu_b2() } ;
    dnudb2_core[] = dInterpolationLinear[SquNorm[$1]]{ Mat_core_nu_b2() } ;
    h_core[] = nu_core[$1] * $1 ;
    dhdb_core[] = TensorDiag[1,1,1] * nu_core[$1#1] + 2*dnudb2_core[#1] * SquDyadicProduct[#1]  ;
    nu[ Ferrite ]  = nu_core[$1];
    dhdb[ Ferrite ] = dhdb_core[$1];
  Else
    nu[ Ferrite ]  = 1./(mur*mu0);
  EndIf
}

Un exemple de résultat possible pourrait être par exemple :

Inductance Inductance


Accouplement magnétique à aimants permanent

Apprès avoir détaillé son principe de fonctionnement, résoudre le problème correspondant à un coupleur à aimants tel que représenté ci-dessous :

Représentation schématique d'un pôle de coupleur considéré

Pour gagner du temps, vous pouvez télécharger la géométrie en cliquant sur ce lien .

Remarque

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 :

Constraint {	
  {Name Anticyclique;
    Case {
       { Region Gauche; Type Link; RegionRef Droite;
		  Coefficient (-1) ; Function Rotate[ Vector[$X,$Y,$Z], 0, 0, -Pi/p ] ; } 
    }
  }
}
Astuce

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.

Information

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.

Problèmes 3D

Conditions de Jauge sur le potentiel vecteur magnétique

Nous avons déjà vu plusieurs fois précédemment que le potentiel vecteur magnétique ${\bf a}$ était défini à un champ de gradient près et qu’il était donc nécessaire de le jauger pour assurer l’unicité de la solution.

Dans le chapitre précédent, deux méthodes ont été présentées et nous allons maintenant voir leur application à nos problèmes ainsi que la manière de les définir dans GetDP.

Jauge de Coulomb

Par application directe de la méthode du chapitre 2, nous allons définir un champ scalaire $\xi$ dont le gradient permettra d’assurer $\text{div}\,{\bf a} = 0$ en l’introduisant comme multiplicateur de Lagrange dans la formulation.

Ainsi la formulation ccomplète dans le domaine discret sera :

$$\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 :

Constraint {	
  {Name Dirichlet ;
    Case {      
      { Region Dirichlet ; Type Assign ; Value 0. ; }
    }
  }
}

FunctionSpace {
  { Name Hrot ; Type Form1 ; // W_0^1
    BasisFunction { 
      { Name se ; NameOfCoef ae ; Function BF_Edge ; 
	     Support Region[{Domaine}]  ; Entity EdgesOf[All] ; }  } 
    Constraint {
      { NameOfCoef ae  ; EntityType EdgesOf ; NameOfConstraint Dirichlet ; }
    }
  }  
  { Name Hgrad ; Type Form0 ; // W_0^0
    BasisFunction { 
      { Name sn ; NameOfCoef xin ; Function BF_Node ; 
	       Support Region[{Domaine}]  ; Entity NodesOf[All] ; }  } 
    Constraint {
      { NameOfCoef xin  ; EntityType NodesOf ; NameOfConstraint Dirichlet ; }
    }
  } 
}

Puis la formulation correspondante (exemple pour un problème sans aimants avec lois de comportement magnétiques linéaires) :

Formulation {
  { Name Magnetostat ; Type FemEquation ;
    Quantity {
      { Name a  ; Type Local  ; NameOfSpace Hrot ; }
      { Name xi ; Type Local  ; NameOfSpace Hgrad ; }
    }
    Equation {
      Integral { [ nu[] * Dof{d a} , {d a} ];
	     In Domaine; Jacobian Jvol; Integration Integ; }
      
      Integral { [ -js[] , {a} ];
	     In Bobines; Jacobian Jvol; Integration Integ; }
      
      Integral{ [ Dof{Grad xi} , {a}  ];
	     In Domaine; Jacobian Jvol; Integration Integ; }
      
      Integral { [ Dof{a} , {Grad xi} ];
	     In Domaine; Jacobian Jvol; Integration Integ; }
    }
  } 
}

Le principal inconvénient de cette méthode est d’aboutir à une taille de système conséquente conduisant à des temps de calcul élévés, c’est pourquoi on lui préfèrera généralement la technique suivante.

Jauge d’arbre

Comme vu dans cette partie, nous allons construire un arbre sur l’ensemble des arrêtes du domaine en partant des surfaces sur lesquelles sont imposées des conditions particulières (dans le cas présent, les conditions de Dirichlet) et imposer une circulation nulle de ${\bf a}$ le long de cet arbre directement dans l’espace fonctionnel associé. Ce dernier sera alors noté $W_{0,\text{JA}}^1$, et notre formulation devient ainsi :

$$\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 :

Constraint {	
  {Name Dirichlet ;
    Case {      
      { Region Dirichlet ; Type Assign ; Value 0. ; }
    }
  }
  { Name Jauge_arbre ; Type Assign ;
    Case {
      {Region Domaine ; SubRegion Dirichlet ; Value 0. ; }
    }
  }
}

FunctionSpace {
  { Name Hrot ; Type Form1 ;  // W_(0,JA)^1
    BasisFunction { 
      { Name se ; NameOfCoef ae ; Function BF_Edge ; 
	      Support Region[{Domaine}]  ; Entity EdgesOf[All] ; }  } 
    Constraint {
      { NameOfCoef ae  ; EntityType EdgesOf ; NameOfConstraint Dirichlet ; }
      { NameOfCoef ae  ; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
	       NameOfConstraint Jauge_arbre ; }
    }
  }  
}

Et la formulation sera strictement la même que dans les cas 2D précédents.


Applications

Retour sur notre inductance pour électronique de puissance

Dans le problème de dimensionnement d’inductance de la page précédente, nous avons utilisé un modèle 2D alors que l’épaisseur de la ferrite utilisée est relativement faible. Il serait donc judicieux de valider (ou non) cette hypothèse par un calcul 3D.

Pour ce faire, je vous fournis le modèle complet (réduit à un huitième de géométrie gràce aux symétries) à Télécharger ici .

Avec, vous pourrez :

  • tester les deux types de jauge en linéaire et non-linéaire ;
  • observer les valeurs des champs tels que l’induction dans le circuit magnétique par exemple : Champ dans l’inductance en 3D Champ dans l’inductance en 3D
  • tester votre design et vérifier si vous respectez le cahier des charges (200 µH sous 15 A) ;
  • tracer la valeur de l’inductance en fonction du courant : Inductance en fonction du courant Inductance en fonction du courant

Retour sur notre accouplement à aimants permanents

L’accouplement étudié page précédente reposait lui aussi sur l’hypothèse que les effets de bords sont négligeables et il serait opportun de le vérifier.

Résoudre le problème en 3D et comparer les deux formulations :

  1. en potentiel vecteur magnétique (avec la jauge que vous voulez) ;
  2. en potentiel scalaire magnétique.

Pour gagner du temps, je vous fournis encore la géométie 3D en cliquant ici . Un exemple de géométrie possible est représenté ci-dessous :

Représentation schématique du coupleur considéré

Magnétodynamique

La magnétodynamique (en anglais magnetodynamics) correspond au calcul de la distribution du champ (ou de l’induction) magnétique et des courants induits produits par des courants variables dans le temps et/ou par des sources en mouvement (aimants ou courants).

Considérons un domaine $\Omega$ de $\mathbb{R}^3$ comportant un sous-domaine $\Omega_c$ conducteur électrique, dont la frontière $\Gamma (=\partial\Omega)$ est comme précédemment divisée en deux morceaux $\Gamma_d$ et $\Gamma_n$.

Les équations qui nous concerneront seront :

  • Les équations de Maxwell : $$\left\{\begin{aligned}\text{div}\,{\bf b} &= 0 \\ {\bf rot}\,{\bf h} &= {\bf j}~\text{dans } \Omega_{c}, \text{ et } {\bf 0} \text{ ailleurs} \\ {\bf rot\,e} &= -\frac{\partial\,{\bf b}}{\partial t} ~\text{dans }\Omega_c\end{aligned}\right.$$
  • Les lois de comportements : $$\left\{\begin{aligned}{\bf b} &= \mu\,{\bf h}&\text{dans }\Omega \\ {\bf j} &= \sigma\,{\bf e} &\text{dans~}\Omega_c\end{aligned}\right.$$
  • Les contions aux limites, par exemple :
    • continuité de la composante normale de l’induction : ${\bf n}\cdot{\bf b} = 0$ sur $\Gamma_d$ ;
    • continuité de la composante tangentielle du champ : ${\bf n}\wedge{\bf h} = {\bf 0}$ sur $\Gamma_n$.

Formulation forte

Le potentiel vecteur magnétique ${\bf a}$ est défini comme précédemment et donc ${\bf b} = {\bf rot\,a}$. En le repportant dans l’équation de Maxwell-Faraday, on obtient :

$${\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.

Formulation faible

En combinant les différentes approches vues dans les sections précédentes, on en déduit la formulation faible complète du problème, dite $({\bf a},v)$ :

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.$$

Sous-sections de Magnétodynamique

Magnétoharmonique

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 » :

Resolution {
  { Name MagnetoHarmonique_2D;
    System {
      { Name A; NameOfFormulation MagnetoHarmonique_2D;
			Type ComplexValue; Frequency freq;
      }
    }
    Operation {
      Generate[A];
      Solve[A]; 
      SaveSolution[A];
    }
  }
}

Et la multiplication par $j\omega$ dans la formulation peut-être faite (au choix) :

  • directement via un terme Complex[0,1]*2*Pi*freq dans les expressions ;
  • ou en utilisant DtDof qui permet de définir une dérivée temporelle des degrés de liberté.

Applications

Barre cylindrique

À 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 téléchargeable ici .

Observer l’effet de peau ainsi que l’évolution de la résistance du conducteur en fonction de la fréquence :

“J efficace dans barre” “J efficace dans barre”

“Résistance en fonction de la fréquence” “Résistance en fonction de la fréquence”

Câble triphasé

Modifier les programmes précédents afin de modéliser une ligne triphasée et observer l’effet de proximité à 50 Hz :

“J efficace dans cable triphasé” “J efficace dans cable triphasé”

Prise en compte du mouvement


Ça arrive bientôt…