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é