From 86631e4b1c66d469cb920f6ccabb38cb217f1085 Mon Sep 17 00:00:00 2001
From: fontchas5 <j.fontchastagner@gmail.com>
Date: Thu, 3 Mar 2022 23:33:16 +0100
Subject: [PATCH] v0.5 fin rajout PyFEMM

---
 TP-FEMM/src/accueil.md  |   7 +-
 TP-FEMM/src/chaufind.md |   4 +-
 TP-FEMM/src/msap.md     | 219 +++++++++++++++++++++++++++++++++++++++-
 TP-FEMM/src/ref.md      |   2 +
 4 files changed, 224 insertions(+), 8 deletions(-)

diff --git a/TP-FEMM/src/accueil.md b/TP-FEMM/src/accueil.md
index 9878015..ac6627b 100644
--- a/TP-FEMM/src/accueil.md
+++ b/TP-FEMM/src/accueil.md
@@ -1,8 +1,11 @@
 
 ## Modélisation & Simulation par EF des dispositifs magnétiques (BC7-EC2a ) : Travaux Pratiques sur le logiciel FEMM
 
-***Julien Fontchastagner - v 0.4 / mars 2022***
- 
+**Julien Fontchastagner - v 0.5 / mars 2022**  
+*Dernière mise à jour : Ajout de PyFEMM (interface avec python)*
+
+<br/>
+
 ---
 
 Bienvenue sur le mini-site dédié aux TP FEMM de l'EC S8-BC7-EC2a pour les 2A ENSEM Formation Énergie.     
diff --git a/TP-FEMM/src/chaufind.md b/TP-FEMM/src/chaufind.md
index bef320c..cbb026d 100644
--- a/TP-FEMM/src/chaufind.md
+++ b/TP-FEMM/src/chaufind.md
@@ -228,7 +228,7 @@ for n=1:npas
 end
 plot(t,Tmoy,'linewidth',2));
 xlabel("Temps (en s)");
-ylabel("Temperature moyenne (en K)");
+ylabel("Temperature moyenne du fond (en K)");
 title("Evolution de la temperature dans la casserole")
 grid on
 ```
@@ -321,7 +321,7 @@ for n in range(npas):
 plt.plot(t,Tmoy,'b',lw=2)
 plt.title('Évolution de la température dans la casserole')
 plt.xlabel('Temps (en s)')
-plt.ylabel('Température moyenne (en K)')
+plt.ylabel('Température moyenne du fond (en K)')
 plt.grid()
 plt.show()
 ```
diff --git a/TP-FEMM/src/msap.md b/TP-FEMM/src/msap.md
index 49e427b..14b8e24 100644
--- a/TP-FEMM/src/msap.md
+++ b/TP-FEMM/src/msap.md
@@ -44,8 +44,9 @@ Une figure représentant celle utilisée sur notre machine est donnée ci-dessou
 
 ### Inplantation (cadeau)
 
-Afin de vous évitez l'étape un peu fastidieuse de définition de la géométrie, je vous donne ci-dessous un exemple de programme définissant entièrement le problème et permettant de calculer la distribution de champ à vide (sans courant d'alimentation). 
+Afin de vous évitez l'étape un peu fastidieuse de définition de la géométrie, je vous donne ci-dessous un exemple de programme définissant entièrement le problème et permettant de calculer la distribution de champ à vide (sans courant d'alimentation).
 
+* *MATLAB :*
 ```Matlab
 close all 
 clear 
@@ -256,15 +257,225 @@ for n = 1:ntheta
 end
 ```
 
+* *Python :*
+```Python
+import numpy as np
+import femm
+
+# Donnees probleme (toutes les longueurs sont en cm)
+p = 3             # nombre de paires de poles
+Lz = 20           # longueur
+Ra = 10           # rayon d'alesage
+alphae = 0.5      # larg. ang. encoche
+alphai = alphae/4 # larg. ang. isthme
+he = 3            # hauteur encoche
+hi = 0.25         # hauteur isthme
+ecs = 1           # epaisseur culasse stator
+e = 0.2           # entrefer
+ea = 1.5          # epaisseur aimant
+ecr = 2           # epaisseur culasse rotor
+alphaa = 0.7      # coef. arc polaire
+
+taup = np.pi/p   # pas polaire
+taud = taup/6 # pas dentaire
+
+N = 20        # nombre de snp.pires
+Ieff = 0      # valeur eff courant d'alim.
+w = 2*np.pi*50   # pulsation
+t = 0
+I1 = np.sqrt(2)*Ieff*np.cos(w*t)
+I2 = np.sqrt(2)*Ieff*np.cos(w*t-2*np.pi/3)
+I3 = np.sqrt(2)*Ieff*np.cos(w*t+2*np.pi/3)
+
+Omega = w/p  # vitesse
+angRot = 0   # position rotor
+angStat = 0  # on laisse 0 : centre sur phase a
+
+femm.openfemm()
+femm.newdocument(0)
+femm.mi_probdef(0,'centimeters','planar',1e-8,Lz) 
+
+# Materiaux
+femm.mi_getmaterial('Air')
+femm.mi_getmaterial('M-19 Steel')
+femm.mi_modifymaterial('M-19 Steel',0,'Toles Stator')
+femm.mi_getmaterial('1010 Steel')
+femm.mi_modifymaterial('1010 Steel',0,'XC10')
+femm.mi_getmaterial('Copper')
+femm.mi_modifymaterial('Copper',0,'Bobine')
+femm.mi_getmaterial('N30')
+femm.mi_modifymaterial('N30',0,'NdFeB')
+
+# Sources
+femm.mi_addcircprop('Phase a',I1,1)
+femm.mi_addcircprop('Phase b',I2,1)
+femm.mi_addcircprop('Phase c',I3,1)
+
+# CL
+femm.mi_addboundprop('Fond culasse',0,0,0,0)
+
+femm.mi_addboundprop('antiper1',0,0,0,0,0,0,0,0,5,0,0)
+femm.mi_addboundprop('antiper2',0,0,0,0,0,0,0,0,5,0,0)
+femm.mi_addboundprop('antiper3',0,0,0,0,0,0,0,0,5,0,0)
+femm.mi_addboundprop('antiper4',0,0,0,0,0,0,0,0,5,0,0)
+
+femm.mi_addboundprop('Bande Roulement',0,0,0,0,0,0,0,0,7,angRot,angStat)
+ # Attention : bug dans femm.mi_addboundprop.m, on doit rajouter :   
+femm.mi_modifyboundprop('Bande Roulement',10,angRot) 
+
+# stator
+R1 = Ra+hi
+R2 = R1+he
+R3 = R2+ecs
+theta1 = (1-alphai)*taud/2
+theta2 = (1-alphae)*taud/2
+theta3 = (1+alphae)*taud/2
+theta4 = (1+alphai)*taud/2
+femm.mi_addnode(Ra,0)
+femm.mi_drawarc(Ra,0,Ra*np.cos(theta1),Ra*np.sin(theta1),theta1*180/np.pi,1)
+femm.mi_drawline(Ra*np.cos(theta1),Ra*np.sin(theta1),R1*np.cos(theta1),R1*np.sin(theta1))
+femm.mi_drawline(R1*np.cos(theta1),R1*np.sin(theta1),R1*np.cos(theta2),R1*np.sin(theta2))
+femm.mi_drawline(R1*np.cos(theta2),R1*np.sin(theta2),R2*np.cos(theta2),R2*np.sin(theta2))
+femm.mi_drawarc(R2*np.cos(theta2),R2*np.sin(theta2),R2*np.cos(theta3),R2*np.sin(theta3),
+           alphae*taud*180/np.pi,5)
+femm.mi_drawline(R2*np.cos(theta3),R2*np.sin(theta3),R1*np.cos(theta3),R1*np.sin(theta3))
+femm.mi_drawline(R1*np.cos(theta3),R1*np.sin(theta3),R1*np.cos(theta4),R1*np.sin(theta4))
+femm.mi_drawline(R1*np.cos(theta4),R1*np.sin(theta4),Ra*np.cos(theta4),Ra*np.sin(theta4))
+
+femm.mi_drawarc(Ra*np.cos(theta4),Ra*np.sin(theta4),Ra*np.cos(taud),Ra*np.sin(taud),
+           theta1*180/np.pi,1)
+r1 = 0.1
+femm.mi_createradius(R1*np.cos(theta2),R1*np.sin(theta2),r1)
+femm.mi_createradius(R1*np.cos(theta3),R1*np.sin(theta3),r1)
+r2 = 0.4
+femm.mi_createradius(R2*np.cos(theta2),R2*np.sin(theta2),r2)
+femm.mi_createradius(R2*np.cos(theta3),R2*np.sin(theta3),r2)
+femm.mi_addsegment(R1*np.cos(theta1),R1*np.sin(theta1),R1*np.cos(theta4),R1*np.sin(theta4))
+femm.mi_selectgroup(0)
+femm.mi_copyrotate(0, 0, taud*180/np.pi, 5 )
+femm.mi_drawline(Ra,0,R3,0)
+femm.mi_drawline(Ra*np.cos(taup),Ra*np.sin(taup), R3*np.cos(taup),R3*np.sin(taup))
+femm.mi_drawarc(R3,0,R3*np.cos(taup),R3*np.sin(taup),taup*180/np.pi,5)
+femm.mi_drawline(Ra,0,Ra-e/3,0)
+femm.mi_drawline(Ra*np.cos(taup),Ra*np.sin(taup),(Ra-e/3)*np.cos(taup),(Ra-e/3)*np.sin(taup))
+femm.mi_drawarc((Ra-e/3),0,(Ra-e/3)*np.cos(taup),(Ra-e/3)*np.sin(taup),taup*180/np.pi,1)
+
+femm.mi_addblocklabel((Ra+R1)/2*np.cos(taud/2),(Ra+R1)/2*np.sin(taud/2))
+femm.mi_selectlabel((Ra+R1)/2*np.cos(taud/2),(Ra+R1)/2*np.sin(taud/2))
+femm.mi_setblockprop('Air',1,0,0,0,0,0)
+femm.mi_clearselected()
+femm.mi_addblocklabel((R2+R1)/2*np.cos(taud/2),(R2+R1)/2*np.sin(taud/2))
+femm.mi_selectlabel((R2+R1)/2*np.cos(taud/2),(R2+R1)/2*np.sin(taud/2))
+femm.mi_setblockprop('Bobine',1,0,'Phase a',0,2,N)
+femm.mi_clearselected()
+femm.mi_addblocklabel((R2+R1)/2*np.cos(taud*1.5),(R2+R1)/2*np.sin(taud*1.5))
+femm.mi_addblocklabel((R2+R1)/2*np.cos(taud*2.5),(R2+R1)/2*np.sin(taud*2.5))
+femm.mi_selectlabel((R2+R1)/2*np.cos(taud*1.5),(R2+R1)/2*np.sin(taud*1.5))
+femm.mi_selectlabel((R2+R1)/2*np.cos(taud*2.5),(R2+R1)/2*np.sin(taud*2.5))
+femm.mi_setblockprop('Bobine',1,0,'Phase b',0,3,-N)
+femm.mi_clearselected()
+femm.mi_addblocklabel((R2+R1)/2*np.cos(taud*3.5),(R2+R1)/2*np.sin(taud*3.5))
+femm.mi_addblocklabel((R2+R1)/2*np.cos(taud*4.5),(R2+R1)/2*np.sin(taud*4.5))
+femm.mi_selectlabel((R2+R1)/2*np.cos(taud*3.5),(R2+R1)/2*np.sin(taud*3.5))
+femm.mi_selectlabel((R2+R1)/2*np.cos(taud*4.5),(R2+R1)/2*np.sin(taud*4.5))
+femm.mi_setblockprop('Bobine',1,0,'Phase c',0,4,N)
+femm.mi_clearselected()
+femm.mi_addblocklabel((R2+R1)/2*np.cos(taud*5.5),(R2+R1)/2*np.sin(taud*5.5))
+femm.mi_selectlabel((R2+R1)/2*np.cos(taud*5.5),(R2+R1)/2*np.sin(taud*5.5))
+femm.mi_setblockprop('Bobine',1,0,'Phase a',0,2,-N)
+femm.mi_clearselected()
+femm.mi_addblocklabel((R2+R3)/2*np.cos(taup/2),(R2+R3)/2*np.sin(taup/2))
+femm.mi_selectlabel((R2+R3)/2*np.cos(taup/2),(R2+R3)/2*np.sin(taup/2))
+femm.mi_setblockprop('Toles Stator',1,0,0,0,1,0)
+femm.mi_clearselected()
+femm.mi_selectsegment(R2,0)
+femm.mi_selectsegment(R2*np.cos(taup),R2*np.sin(taup))
+femm.mi_setsegmentprop('antiper1',0,1,0,1)
+femm.mi_clearselected()
+femm.mi_selectsegment((Ra-e/6),0)
+femm.mi_selectsegment((Ra-e/6)*np.cos(taup),(Ra-e/6)*np.sin(taup))
+femm.mi_setsegmentprop('antiper2',0,1,0,1)
+femm.mi_clearselected()
+femm.mi_selectarcsegment(R3*np.cos(taup/2),R3*np.sin(taup/2))
+femm.mi_setarcsegmentprop(5,'Fond culasse',0,1) 
+femm.mi_clearselected()
+
+# rotor
+R1 = Ra-2*e/3
+R2 = Ra-e
+R3 = R2-ea
+R4 = R3-ecr
+theta1 = (1-alphaa)*taup/2
+theta2 = (1+alphaa)*taup/2
+femm.mi_drawline(R4,0,R3,0)
+femm.mi_drawline(R3,0,R1,0)
+femm.mi_drawline(R4*np.cos(taup),R4*np.sin(taup),R3*np.cos(taup),R3*np.sin(taup))
+femm.mi_drawline(R3*np.cos(taup),R3*np.sin(taup),R1*np.cos(taup),R1*np.sin(taup))
+femm.mi_drawline(R3*np.cos(theta1),R3*np.sin(theta1),R2*np.cos(theta1),R2*np.sin(theta1))
+femm.mi_drawline(R3*np.cos(theta2),R3*np.sin(theta2),R2*np.cos(theta2),R2*np.sin(theta2))
+femm.mi_addarc(R4,0,R4*np.cos(taup),R4*np.sin(taup),taup*180/np.pi,5)
+femm.mi_addarc(R3,0,R3*np.cos(theta1),R3*np.sin(theta1),theta1*180/np.pi,5)
+femm.mi_addarc(R3*np.cos(theta2),R3*np.sin(theta2),R3*np.cos(taup),R3*np.sin(taup),
+               theta1*180/np.pi,5)
+femm.mi_addarc(R3*np.cos(theta1),R3*np.sin(theta1),R3*np.cos(theta2),R3*np.sin(theta2),
+               alphaa*taup*180/np.pi,5)
+femm.mi_addarc(R2*np.cos(theta1),R2*np.sin(theta1),R2*np.cos(theta2),R2*np.sin(theta2),
+               alphaa*taup*180/np.pi,1)
+femm.mi_addarc(R1,0,R1*np.cos(taup),R1*np.sin(taup),taup*180/np.pi,1)
+
+femm.mi_addblocklabel((R2+R3)/2*np.cos(theta1/2),(R2+R3)/2*np.sin(theta1/2))
+femm.mi_selectlabel((R2+R3)/2*np.cos(theta1/2),(R2+R3)/2*np.sin(theta1/2))
+femm.mi_setblockprop('Air',1,0,0,0,5,0)
+femm.mi_clearselected()
+femm.mi_addblocklabel((R4+R3)/2*np.cos(taup/2),(R4+R3)/2*np.sin(taup/2))
+femm.mi_selectlabel((R4+R3)/2*np.cos(taup/2),(R4+R3)/2*np.sin(taup/2))
+femm.mi_setblockprop('XC10',1,0,0,0,5,0)
+femm.mi_clearselected()
+femm.mi_addblocklabel((R2+R3)/2*np.cos(taup/2),(R2+R3)/2*np.sin(taup/2))
+femm.mi_selectlabel((R2+R3)/2*np.cos(taup/2),(R2+R3)/2*np.sin(taup/2))
+femm.mi_setblockprop('NdFeB',1,0,0,taup/2*180/np.pi,5,0)
+femm.mi_clearselected()
+femm.mi_selectsegment(R2,0)
+femm.mi_selectsegment(R2*np.cos(taup),R2*np.sin(taup))
+femm.mi_setsegmentprop('antiper3',0,1,0,5)
+femm.mi_clearselected()
+femm.mi_selectsegment((R4+R3)/2,0)
+femm.mi_selectsegment((R4+R3)/2*np.cos(taup),(R4+R3)/2*np.sin(taup))
+femm.mi_setsegmentprop('antiper4',0,1,0,5)
+femm.mi_clearselected()
+femm.mi_selectarcsegment(R4*np.cos(taup/2),R4*np.sin(taup/2))
+femm.mi_setarcsegmentprop(5,'Fond culasse',0,1) 
+femm.mi_clearselected()
+
+femm.mi_selectarcsegment((Ra-2*e/3)*np.cos(taup/2),(Ra-2*e/3)*np.sin(taup/2))
+femm.mi_selectarcsegment((Ra-e/3)*np.cos(taup/2),(Ra-e/3)*np.sin(taup/2))
+femm.mi_setarcsegmentprop(1,'Bande Roulement',0,1)
+femm.mi_clearselected()
+
+femm.mi_zoomnatural()
+femm.mi_saveas('masap.FEM')
+
+ntheta = 30
+dtheta = 2 
+for n in range(ntheta):
+    femm.mi_analyze()
+    femm.mi_loadsolution()
+    femm.mo_hidepoints()
+    femm.mo_showdensityplot(0,0,2,0,'mag')
+    femm.mo_savebitmap("imagesMag{}.bmp".format(n))
+    angRot = angRot+dtheta
+    femm.mi_modifyboundprop('Bande Roulement',10,angRot)
+```
+
+
 Avec les résultats obtenus, on peut réaliser la petite animation suivante : 
 
 ![Flux à vide](./figures/masap.gif "Flux à vide dans la machine")
 
+<br/>
 
 **➤ Analysez le programme fourni et commentez les différentes parties.**
 
-**➤ Ayant volontairement laissé la discrétisation automatique, regardez l'influence sur les résultats (induction, flux, énergie, coénergie) ; et choisissez un bon compromis finesse / pertinence des résultats.**
-
 **➤ Calculez et tracez les flux à vide, puis déduisez-en les fem à vide.**
 
 
@@ -298,7 +509,7 @@ Le calcul du couple exercé sur le rotor se fait de la même façon que celui de
 
 <br/>
 
->**Félicitations !**
+> **Félicitations !**
 >
 >**Si vous avez tout assimilé et tout compris vous êtes devenu un bon utilisateur de FEMM. Avec de la pratique, vous pourrez même devenir un *power user* !**
 
diff --git a/TP-FEMM/src/ref.md b/TP-FEMM/src/ref.md
index f68c0f4..244ee97 100644
--- a/TP-FEMM/src/ref.md
+++ b/TP-FEMM/src/ref.md
@@ -6,6 +6,8 @@
 
 * D. Meeker, [Finite Element Method Magnetics: OctaveFEMM, User’s Manual](https://www.femm.info/wiki/Files/files.xml?action=download&file=octavefemm.pdf), 17 mars 2018.
 
+* D. Meeker, [Finite Element Method Magnetics: pyFEMM 0.1.3, User’s Manual](https://www.femm.info/wiki/pyFEMM/manual.pdf), 20 juillet 2021. 
+
 <br/>
 
 # Crédits
-- 
GitLab