Table des matières Python

Aimant et bobine : flux et force - Étude expérimentale

1. Introduction

Ce document montre une étude expérimentale de l'interaction entre un aimant permanent et une bobine. Lorsque la bobine n'est pas alimentée, on mesure la force électromotrice (f.é.m.) dans la bobine causée par le mouvement de l'aimant, dont on déduit le flux magnétique généré par l'aimant dans la bobine pour différentes positions relatives. Lorsqu'un courant électrique connu circule dans la bobine, on mesure la force sur l'aimant pour différentes positions relatives.

Un modèle bidimensionnel correspondant à ce problème est traité dans Aimant et bobine : flux et force.

Dans l'hypothèse où l'aimantation de l'aimant est constante (aimant idéal), celui-ci se comporte comme une distribution de courants surfaciques, c'est-à-dire en fin de compte comme une bobine parcourue par un courant d'intensité constante. En l'absence de pièces en ferromagnétique doux, l'interaction entre l'aimant et la bobine est donc similaire à l'interaction entre deux bobines (système linéaire). Si l'on note Φ le flux magnétique dans la bobine produit par l'aimant (ce flux ne contient pas le flux propre) et I l'intensité du courant dans la bobine, L l'auto-inductance de la bobine et La celle de l'aimant, l'énergie du champ magnétique s'écrit :

Ec=IΦ+12LI2+12LaIa2(1)

Comme démontré dans Flux magnétique généré par un aimant permanent, la résultante des forces magnétiques sur l'aimant est égale au gradient de cette énergie par rapport à la position de l'aimant (résultat valable pour un déplacement des bobines à courant constant) :

F=gradEc(2)

En l'absence de pièces en ferromagnétique doux, les auto-inductances ne varient pas et on a donc :

F=grad(IΦ)(3)

L'objectif de l'étude expérimentale est de vérifier la pertinence de ce résultat théorique, sachant que l'aimant réel n'est pas idéal.

En présence de pièces en ferromagnétique doux qui se déplacent avec l'aimant, ce résultat n'est plus valable car le coefficient L varie lorsque la pièce en fer se déplace par rapport à la bobine. L'auto-inductance de la bobine est cependant mesurable pour différentes positions relatives. Les deux termes de l'énergie qui varient sont donc accessibles par des mesures électriques (de forces électromotrices). En revanche, lorsque une pièce en fer doux se déplace par rapport à l'aimant (par exemple le noyau de la bobine), le terme faisant intervenir La varie à cause de la variation de cette auto-inductance. Or ce terme de l'énergie n'est évidemment pas accessible par des mesures électriques.

2. Flux magnétique dans la bobine

La tension aux bornes de la bobine est enregistrée avec l'Analog Discovery 2 et l'interface Python présentée dans Interface Python pour Analog Discovery.

Le mouvement de l'aimant (fait à la main dans cette partie) est filmé avec une caméra Basler déclenchée par un signal émis par un Arduino (Déclenchement d'une caméra Basler). L'enregistrement vidéo et celui de la tension sont déclenchés simultanément par pression sur un bouton poussoir.

Voici le script qui effectue l'enregistrement de la tension :

acquisitionTension.py
import numpy as np
import matplotlib.pyplot as plt
from analog_0_1 import Device,AnalogInput,AnalogOutput


device = Device(-1)
device.open()

analog = AnalogInput(device)
analog.channels([0],[0.5])
fSamp = 400
size = 8000
t = analog.sampling(fSamp,size)
analog.externalTrigger(1,timeout=0,position=0.0)

voltage = analog.record()
u0 = voltage[0,:]
plt.figure()
plt.plot(t,u0)
plt.grid()
np.savetxt("mouvement-1.txt",np.array([t,u0]).T)
plt.show()
device.close()
            

L'enregistrement contient 8000 échantillons et la fréquence d'échantillonnage est 400 Hz, soit une durée de 20 secondes. L'acquisition est déclenchée par un front montant sur l'entrée T1 de l'AD, reliée au bouton poussoir qui commande aussi le démarrage de l'enregistrement vidéo. La position du déclenchement est en fait au milieu de la fenêtre de 20 secondes (position=0.0), ce qui fait que l'enregistrement vidéo démarre à t=10s par rapport au début de l'enregistrement de la tension. On laissera donc une durée d'au moins 10 secondes entre le lancement du script et l'appuie sur le bouton poussoir qui provoque le déclenchement. Les 10 premières secondes de l'enregistrement se font sans déplacement de l'aimant, celui-ci étant assez loin pour que le flux magnétique soit négligeable. Cette partie du signal permettra de faire une compensation d'offset précise, indispensable pour l'intégration de la tension. Le mouvement de l'aimant se fait pendant la seconde partie, d'une durée de 10 secondes également.

Le signal de déclenchement image par image de la caméra Basler a une fréquence 40 images par secondes ou 20 images par secondes. La fréquence d'échantillonnage de la vidéo est donc 10 fois ou 20 fois plus petite que la fréquence d'échantillonnage de la tension.

La tension aux bornes de la bobine est égale à la force-électromotrice e(t). D'après la loi de Faraday, l'intégration de celle-ci à partir d'une position où le flux est nul donne le flux :

Φ(t)=-0te(t')dt'(4)

L'intégration numérique est faite avec le filtre récursif défini par :

yn=yn-1+Te2(xn+xn-1)(5)

Te est la période d'échantillonnage. L'intégration se fait en fait à partir de l'instant t=10s (milieu de l'enregistrement du signal). Pour éviter une dérive du flux due à une tension numérique non nulle lorsque la f.é.m. est nulle, on retranche au signal numérique sa moyenne calculée sur les 10 premières secondes.

Le traitement de la vidéo et du signal numérique consiste à calculer le flux par intégration et à tracer une courbe de flux synchronisée avec le mouvement de l'aimant visible sur la vidéo. Le traitement de la vidéo est fait avec OpenCV. Le script suivant comporte une classe Python qui permet de tracer une courbe sur une image, qui peut être incluse dans une vidéo :

cvplot_1.py
class CvPlot:
    def __init__(self,width,height,left,right,top,bottom,extent):
        self.frame = np.zeros((height,width,3),np.uint8)
        self.width = width
        self.height = height
        self.left = left
        self.right = right
        self.top = top
        self.bottom = bottom
        self.plot_width = width-left-right
        self.plot_height = height-top-bottom
        self.xmin = extent[0]
        self.xmax = extent[1]
        self.ymin = extent[2]
        self.ymax = extent[3]
        self.ax = self.plot_width / (self.xmax-self.xmin)
        self.ay = self.plot_height / (self.ymax-self.ymin)
        self.font = cv.FONT_HERSHEY_SIMPLEX
        self.ticksize = 5
        self.background = (0,0,0)
        self.framecolor = (255,255,255)
        

    def xy2pixel(self,x,y):
        px = (x-self.xmin)*self.ax+self.left
        py = self.height-self.bottom-(y-self.ymin)*self.ay
        return int(px),int(py)

    def clear(self):
        self.frame = np.zeros((self.height,self.width,3),np.uint8)
        for c in range(3):
            self.frame[:,:,c] = self.background[c]
        

    def drawFrame(self,linewidth=1):
        yy = [self.height-self.bottom,self.top]
        for y in yy:
            cv.line(self.frame,(self.left,y),(self.width-self.right,y),self.framecolor,linewidth)
        xx = [self.left,self.width-self.right]
        for x in xx:
            cv.line(self.frame,(x,self.top),(x,self.height-self.bottom),self.framecolor,linewidth)

    def drawCurve(self,xx,yy,color=(255,255,255),linewidth=1):
        N = len(xx)
        if N >= 2:
            px0,py0 = self.xy2pixel(xx[0],yy[0])
            px1,py1 = self.xy2pixel(xx[1],yy[1])
            cv.line(self.frame,(px0,py0),(px1,py1),color,linewidth)
            for k in range(2,N):
                px0,py0 = px1,py1
                px1,py1 = self.xy2pixel(xx[k],yy[k])
                cv.line(self.frame,(px0,py0),(px1,py1),color,linewidth)

    def labels(self,xlabel,ylabel,size=0.5):
        cv.putText(self.frame,xlabel,(self.width//2,self.height-self.bottom//4),self.font,size,self.framecolor)
        cv.putText(self.frame,ylabel,(5,self.height//2),self.font,size,self.framecolor)

    def xticks(self,N,fmt="%f",linewidth=1,size=0.5,full=False):
        xx = np.linspace(self.xmin,self.xmax,N+1)
        for x in xx:
            px,py = self.xy2pixel(x,self.ymin)
            if full:
                cv.line(self.frame,(px,py),(px,py-self.plot_height),self.framecolor,linewidth)
            else:
                cv.line(self.frame,(px,py),(px,py-self.ticksize),self.framecolor,linewidth)
            s = fmt%(x)
            cv.putText(self.frame,s,(px-int(10*len(s)*size),py+int(30*size)),self.font,size,self.framecolor)

    def yticks(self,N,fmt="%f",linewidth=1,size=0.5,full=False):
        yy = np.linspace(self.ymin,self.ymax,N+1)
        for y in yy:
            px,py = self.xy2pixel(self.xmin,y)
            cv.line(self.frame,(px,py),(px+self.ticksize,py),self.framecolor,linewidth)
            if full:
                cv.line(self.frame,(px,py),(px+self.plot_width,py),self.framecolor,linewidth)
            else:
                cv.line(self.frame,(px,py),(px+self.ticksize,py),self.framecolor,linewidth)
            s = fmt%(y)
            cv.putText(self.frame,s,(px-int(len(s)*20*size),py),self.font,size,self.framecolor)
            
    def getFrame(self):
        return self.frame               
                

Voici le script de traitement du signal échantillonné et de la vidéo. On ajoute le tracé du flux magnétique dans la bobine en fonction du temps et de la f.é.m. en fonction du temps. Le fichier vidéo généré est encodé en format MP4 avec le codec XVID. Le nom du fichier vidéo source doit être identique au nom du fichier texte qui contient le signal, le premier ayant une extension TXT, le second une extension AVI.

analyseMouvementAimant.py
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
from cvplot_1 import CvPlot

nom = "mouvement-6-20fps"
fechant = 20 # période d'échantillonnage de la vidéo
techant = 1/fechant
video = cv.VideoCapture("%s.avi"%nom)
n= 0
success,frame = video.read()
cv.namedWindow("frame")
(height,width,c)=frame.shape
[t,e] = np.loadtxt("%s.txt"%nom,unpack=True)
full_width = width
full_height = 3*height

i = 0
t_trigger = 10.0
while t[i] < t_trigger: i+=1
te = t[i+1]-t[i] # période d'échantillonnage de la tension
e -= e[0:i].mean()
print("i = %d"%i)
plt.figure()
plt.plot(t,e)
flux = [0]
fem = [0]
for j in range(1,len(e)-i):
    flux.append(flux[j-1]-te*0.5*(e[i+j]+e[i+j-1]))
    fem.append(e[i+j])
flux = np.array(flux)
fem = np.array(fem)
plt.figure()
temps = np.arange(len(flux))*te
plt.plot(temps,flux)
plt.plot(temps,fem)
plt.show()
step = int(techant/te)
print("step = %d"%step)
t = 0
color = (255,0,0)

flux_max = np.max(np.absolute(flux))*1.5
fem_max = np.max(np.absolute(fem))*1.5
flux_max = 0.06
fem_max = 0.25
plot_flux = CvPlot(width,height,100,50,50,60,[0,10,-flux_max,flux_max])
plot_fem = CvPlot(width,height,100,50,50,60,[0,10,-fem_max,fem_max])
fourcc = cv.VideoWriter_fourcc(*'XVID')
#fourcc = cv.VideoWriter_fourcc(*'MP4V') #même codec
out = cv.VideoWriter("%s.mp4"%nom, fourcc, 30.0, (full_width,full_height))
i=0
while success:
    full_frame = np.zeros((full_height,full_width,3),np.uint8)
    plot_flux.clear()
    plot_flux.drawFrame()
    plot_fem.clear()
    plot_fem.drawFrame()
    plot_flux.labels("t (s)","Phi (Wb)",size=0.5)
    plot_flux.xticks(10,"%0.1f")
    plot_flux.yticks(2,"%0.2g",full=True)
    plot_fem.labels("t (s)","fem (V)",size=0.5)
    plot_fem.xticks(10,"%0.1f")
    plot_fem.yticks(2,"%0.2g",full=True)
    if i>=1:
        if i>=len(fem): i = len(fem)-1
        temps = np.arange(i)*te
        plot_flux.drawCurve(temps,flux[0:i],(0,0,255),linewidth=2)
        plot_fem.drawCurve(temps,fem[0:i],(0,0,255),linewidth=2)
    
    full_frame[0:height,0:width,:] = frame
    full_frame[height:2*height,0:width,:] = plot_flux.getFrame()
    full_frame[2*height:3*height,0:width,:] = plot_fem.getFrame()
    i += step
    cv.imshow('frame',full_frame)
    out.write(full_frame)
    key = cv.waitKey()
    if key == ord('q'):
        break
    success,frame = video.read()
    n += 1
    t += techant

out.release()
                 

Le fichier MP4 obtenu est converti en fichier WEBM (codec VP9) pour la diffusion sur le web. Voici la commande FFMPEG permettant de faire cette conversion :

ffmpeg -i mouvement-6-20fps.mp4 -b:v 500k -r 20 mouvement-6-20fps.webm

L'option -b:v indique le taux de bits (bitrate) et l'option -r indique la fréquence d'image finale, ici 20 Hz. L'attribution d'une fréquence d'image égale à la fréquence de prise de vue n'est pas obligatoire mais elle permet d'obtenir une lecture de la vidéo à la vitesse réelle. Bien entendu, si la fréquence de prise de vue est élevée (par exemple 100 fps), on peut prendre -r 30 de manière à obtenir une lecture au ralenti.

Voici un exemple : un aimant droit (AlNiCo) est déplacé sur l'axe d'une bobine de 250 spires :

Voici le même déplacement effectué plus rapidement :

Voici le déplacement d'un aimant dans la direction transversale par rapport à l'axe de la bobine. L'aimant est fait d'un assemblage de deux aimants plats en ferrite :

Voici le déplacement d'un aimant constitué de trois aimants plats :

3. Distance aimant-bobine

Afin de déterminer la distance entre l'aimant et la bobine, le chariot sur lequel la jauge de contrainte et l'aimant sont fixés est entrainé par une vis sans fin reliée à un moteur pas à pas. Celui-ci est piloté par un driver dont la commande consiste à choisir le sens de rotation (DIR) et à envoyer une impulsion sur STEP pour chaque pas. Le déplacement est de 10 mm pour 1000 pas. Avec une fréquence d'impulsions de 100 Hz, on parvient donc à un déplacement de 1 mm par secondes, soit 100 s pour réaliser un déplacement de 100 mm, suffisant pour observer une grande variation de la force.

L'Arduino (fixé sur le chariot) pilote le moteur et effectue l'enregistrement de la force. Le démarrage du moteur et de l'enregistrement sont déclenchés simultanément par appuie sur un bouton poussoir, qui commande également un autre arduino chargé de générer les impulsions pour le pilotage de la caméra.

Les impulsions pour la commande du moteur pas à pas sont générées par le Timer 1. La lecture les données fournies par l'amplificateur HX711 sur lequel la jauge de contrainte est reliée est faite tous les 100 ms, au moyen d'interruptions générés par le Timer 3. Dans le gestionnaire de cette interruption, la donnée 24 bits est lue puis transmise immédiatement au PC sous forme binaire (non convertie en newton). L'appuie sur le bouton poussoir (relié à l'entrée START) provoque par interruption l'exécution de la fonction start, qui déclenche le Timer 1 et l'acquisition de la force, avec transmission par port série au PC.

railMoteurPAP-HX711.ino
#define STEP 11
#define DIR 8
#define EN 9
#define START 3
#define DOUT 4
#define SCK 5
#define GAIN 1 // 1 = 128 channel A, 3 = 64 channel A, 2 = 32 channel B
#define PERIOD 10000 // période en microsecondes
#define NSTEPS 10000 // 1000 pas par cm, 0.1 s par secondes
#define COMPTEUR true

#define START_TRANSMISSION 100
#define STOP_TRANSMISSION 101

uint16_t diviseur[6] = {0,1,8,64,256,1024};
volatile uint32_t compteur;
bool actif;
int dir;
volatile bool transmission;
bool triggered;
volatile int32_t transbuffer[0];
volatile uint8_t d[3];
volatile unsigned long valeur;
volatile uint16_t indice;
volatile uint8_t filler = 0x00;

void init_pwm_timer1(uint32_t period) {
    char clockBits;
    uint32_t icr;
    cli();
    TCCR1A = 0;
    TCCR1A |= (1 << COM1A1); //Clear OC1A on compare match when upcounting, set OC1A on compare match when downcounting
    TCCR1A |= (1 << COM1B1);
#if defined(__AVR_ATmega2560__) || defined(__AVR_ATmega32U4__)
    TCCR1A |= (1 << COM1C1);
#endif
    TCCR1B = 1 << WGM13; // phase and frequency correct pwm mode, top = ICR1
    int d = 1;
    icr = (F_CPU/1000000*period/2);
    while ((icr>0xFFFF)&&(d<5)) { // choix du diviseur d'horloge
        d++;
        icr = ((F_CPU/1000000)*period/2/diviseur[d]);
   }
   clockBits = d;
   ICR1 = icr; // valeur maximale du compteur
   TIMSK1 = 1 << TOIE1; // overflow interrupt enable
   TCNT1 = 0; // mise à zéro du compteur
   OCR1A = icr*0.5;
   TCCR1B |= clockBits; // déclenchement du compteur
   sei(); // activation des interruptions
}  

void init_pwm_timer3(uint32_t period) {
    char clockBits;
    uint32_t icr;
    cli();
    TCCR3A = 0;
    TCCR3A |= (1 << COM1A1); //Clear OC1A on compare match when upcounting, set OC1A on compare match when downcounting
    TCCR3A |= (1 << COM1B1);
#if defined(__AVR_ATmega2560__) || defined(__AVR_ATmega32U4__)
    TCCR3A |= (1 << COM1C1);
#endif
    TCCR3B = 1 << WGM13; // phase and frequency correct pwm mode, top = ICR1
    int d = 1;
    icr = (F_CPU/1000000*period/2);
    while ((icr>0xFFFF)&&(d<5)) { // choix du diviseur d'horloge
        d++;
        icr = ((F_CPU/1000000)*period/2/diviseur[d]);
   }
   clockBits = d;
   ICR3 = icr; // valeur maximale du compteur
   TIMSK3 = 1 << TOIE1; // overflow interrupt enable
   TCNT3 = 0; // mise à zéro du compteur
   OCR3A = icr*0.5;
   TCCR3B |= clockBits; // déclenchement du compteur
   sei(); // activation des interruptions
}    
    		  
ISR(TIMER1_OVF_vect) { // impulsion de commande du moteur pas à pas
    if (COMPTEUR) {
      compteur +=1;
      if (compteur == NSTEPS) {
          OCR1A = 0;
          digitalWrite(EN,HIGH);
          actif = false;
          
      }
    }
}

ISR(TIMER3_OVF_vect) { // lecture de la force
  
  
  
   //while (digitalRead(DOUT)==HIGH) {;}
   d[2] = shiftIn(DOUT,SCK,MSBFIRST);
   d[1] = shiftIn(DOUT,SCK,MSBFIRST);
   d[0] = shiftIn(DOUT,SCK,MSBFIRST);
   for (unsigned int i=0; i<GAIN; i++) {
      digitalWrite(SCK,HIGH);
      delayMicroseconds(1);
      digitalWrite(SCK,LOW);
      delayMicroseconds(1);
   }
  if (d[2] & 0x80) {
    filler = 0xFF;
  } else {
    filler = 0x00;
  }
   valeur = ( static_cast<unsigned long>(filler) << 24
     | static_cast<unsigned long>(d[2]) << 16
      | static_cast<unsigned long>(d[1]) << 8
      | static_cast<unsigned long>(d[0]) );
   long v = static_cast<long>(valeur);
   if (transmission&&triggered) {
        transbuffer[0] = v;
        Serial.write((uint8_t *)transbuffer,4);
     }
   
}      		  



void start() {
  if (actif) return;
  if (dir) {
    digitalWrite(DIR,HIGH);
    dir = 0;
  }
  else {
    digitalWrite(DIR,LOW);
    dir = 1;
  }
  dir = 0;
  digitalWrite(EN,LOW);
  compteur = 0;
  init_pwm_timer1(PERIOD);
  actif = true;
  triggered = true;
}

void start_transmission() {
  char c;
  c = 0;
  Serial.write(c);
  c = 255;
  Serial.write(c);
  c = 0;
  Serial.write(c);
  transmission = true;
  triggered = false;
}

void stop_transmission() {
  transmission = false;

} 

void setup() {
  Serial.begin(115200);
  Serial.setTimeout(0);
  char c = 0;
  Serial.write(c);
  c = 255;
  Serial.write(c);
  c = 0;
  Serial.write(c);
  pinMode(STEP,OUTPUT);
  digitalWrite(STEP,LOW);
  pinMode(DIR,OUTPUT);
  digitalWrite(DIR,LOW);
  pinMode(EN,OUTPUT);
  digitalWrite(EN,HIGH);
  pinMode(START,INPUT);
  pinMode(DOUT,INPUT);
  pinMode(SCK,OUTPUT);
  digitalWrite(SCK,LOW);
  triggered = false;
  attachInterrupt(digitalPinToInterrupt(START),start,RISING);
  actif = false;
  transmission = false;
  dir = 1;
  init_pwm_timer3(100000);
}

void loop() {
   char com;
   if (Serial.available()>0) {
      com = Serial.read();
      if (com==START_TRANSMISSION) start_transmission();
      else if (com==STOP_TRANSMISSION) stop_transmission();
   }
}    		  
       

Voici le fichier contenant la classe permettant de récupérer la force au cours du temps. La constante ECHELLE a été déterminée par un étalonnage préalable au moyen d'un poids de 50 g (voir Jauge de déformation). Elle permet d'obtenir la force en mg mais sans correction du décalage.

arduinoHX711.py

import numpy
import serial
import time

class Arduino():
    def __init__(self,port):
        self.ser = serial.Serial(port,baudrate=115200)
        c_recu = self.ser.read(1)
        while ord(c_recu)!=0:
            c_recu = self.ser.read(1)
        c_recu = self.ser.read(1)
        while ord(c_recu)!=255:
            c_recu = self.ser.read(1)
        c_recu = self.ser.read(1)
        while ord(c_recu)!=0:
            c_recu = self.ser.read(1)
        self.START_TRANSMISSION = 100
        self.STOP_TRANSMISSION = 101
        self.FECHANT = 10.0
        self.ECHELLE = 0.07520
        
    def close(self):
        self.ser.close()
    def write_int8(self,v):
    	char = int(v&0xFF) # nécessaire pour les nombres négatifs
    	self.ser.write((char).to_bytes(1,byteorder='big'))
    def write_int16(self,v):
        v = numpy.int16(v)
        char1 = int((v & 0xFF00) >> 8)
        char2 = int((v & 0x00FF))
        self.ser.write((char1).to_bytes(1,byteorder='big'))
        self.ser.write((char2).to_bytes(1,byteorder='big'))
    def write_int32(self,v):
        v = numpy.int32(v)
        char1 = int((v & 0xFF000000) >> 24)
        char2 = int((v & 0x00FF0000) >> 16)
        char3 = int((v & 0x0000FF00) >> 8)
        char4 = int((v & 0x000000FF))
        self.ser.write((char1).to_bytes(1,byteorder='big'))
        self.ser.write((char2).to_bytes(1,byteorder='big'))
        self.ser.write((char3).to_bytes(1,byteorder='big'))
        self.ser.write((char4).to_bytes(1,byteorder='big'))
    def start_transmission(self):
        self.write_int8(self.START_TRANSMISSION)
        c_recu = self.ser.read(1)
        while ord(c_recu)!=0:
            c_recu = self.ser.read(1)
        c_recu = self.ser.read(1)
        while ord(c_recu)!=255:
            c_recu = self.ser.read(1)
        c_recu = self.ser.read(1)
        while ord(c_recu)!=0:
            c_recu = self.ser.read(1)
    def stop_transmission(self):
        self.write_int8(self.STOP_TRANSMISSION)
    
    def lecture(self):
        buf = self.ser.read(4)
        x = buf[3]*0x1000000 + buf[2]*0x10000 + buf[1]*0x100 + buf[0]
        if buf[3]&0x40: #bit de signe = 1
            x = x-0x100000000
        x = x*self.ECHELLE
        return x
       

Le script suivant permet de faire une acquisition de la force. La durée est 100 secondes, égale à la durée de fonctionnement du moteur. L'acquisition ne démarre que lorsqu'on appuie sur le bouton poussoir

arduinoHX711-animate.py

import numpy
from matplotlib.pyplot import *
import matplotlib.animation as animation
from arduinoHX711 import *

ard = Arduino("COM3")


fe = 10.0
duree = 100
te = 1/fe
N = int(duree/te)

t = numpy.linspace(0,duree,N)
F = numpy.zeros(N)

fig,ax = subplots()
line0, = ax.plot(t,F)
ax.grid()
ax.set_xlabel("t (s)")
ax.set_ylabel("F (mg)")


ard.start_transmission()
n=0
temps = numpy.array([])
force = numpy.array([])
t = 0.0
def animate(i):
    global ax,line0,ard,n,N,te,temps,force,t
    if n<N:
        f = ard.lecture()
        temps = numpy.append(temps,t)
        force = numpy.append(force,f)
        t += te
        n += 1
        line0.set_xdata(temps)
        line0.set_ydata(force)
        ax.axis([0,temps[n-1],force.min(),force.max()])

ani = animation.FuncAnimation(fig,animate,N,interval=te*1000*0.8)
show()
ard.stop_transmission()
ard.close()
numpy.savetxt("force-1.txt",numpy.array([temps,force]).T,delimiter="\t",fmt="%.6e",header="t (s)\t F (mg)")   		 
    		 
        

La force (en mg) et le temps sont enregistrés dans un fichier texte.

L'expérience est conduite de la manière suivante. Initialement, la face de l'aimant se trouve à 100 mm de la face d'entrée de la bobine (le déplacement est de 100 mm); le courant dans la bobine est nul mais l'alimentation est réglée en régulation de courant à 2 A, avec une tension nulle. Le script arduinoHX711-animate.py est lancé (l'enregistrement vidéo aussi) puis on appuie sur le bouton poussoir, ce qui déclenche le déplacement motorisé de l'aimant, l'acquisition de la force (et l'enregistrement vidéo à une fréquence de 4 images par secondes). Après quelques secondes, on augmente la tension de l'alimentation de la bobine pour faire passer le courant de 2 A. La première partie de l'enregistrement de la force (moins de 10 secondes) servira pour le calcul du décalage à appliquer à la force (puisque la force est nulle en l'absence de courant dans la bobine). Bien sûr, cette procédure ne sera pas valable en présence d'un noyau en fer doux dans la bobine : il faudra initialement placer la bobine et le noyau assez loin de l'aimant pour que la force soit négligeable. Le réglage du zéro, c'est-à-dire la détermination de la valeur donnée par le convertisseur A/N qui correspond à une force nulle, pourrait être faite une fois pour toute au préalable.

Le script suivant effectue le calcul du décalage (à partir des 3 premières secondes) et trace la force (en mN) en fonction de la distance (en mm) dans une vidéo :

analyseMouvementAimant-2.py
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
from cvplot_1 import CvPlot

fechant = 4
techant = 1/fechant
nom = "mouvement-4"
force = "force-4"
video = cv.VideoCapture("%s-4fps.avi"%nom)
n= 0
success,frame = video.read()
cv.namedWindow("frame")
(height,width,c)=frame.shape
full_width = width
plot_height = int(height*0.7)
full_height = height + plot_height
[temps,F] = np.loadtxt("%s.txt"%force,unpack=True,skiprows=1)
F *= 1e-3*9.81 # conversion en mN
F -= np.mean(F[0:30])
# déplacement de 100 mm à 100 pas par secondes (1000 pas pour 1 cm), soit 1 mm par seconde
x = temps
plt.figure()
plt.plot(x,F)
plt.xlabel('x (mm)')
plt.ylabel('F (mN)')
plt.grid()
plt.show()
te = temps[1]-temps[0]


plot_force = CvPlot(width,plot_height,100,50,50,60,[0,100,0,500])
fourcc = cv.VideoWriter_fourcc(*'MP4V')
out = cv.VideoWriter("%s.mp4"%nom, fourcc, 30.0, (full_width,full_height))

t = 0
while success:
    full_frame = np.zeros((full_height,full_width,3),np.uint8)
    plot_force.clear()
    plot_force.drawFrame()
    plot_force.labels("x (mm)","F (mN)",size=0.5)
    plot_force.xticks(10,"%0.1f")
    plot_force.yticks(1,"%0.2g",full=True)
    n = min(int(t/te),len(F)-1)
    plot_force.drawCurve(x[0:n],F[0:n],(0,0,255),linewidth=2)
    full_frame[0:height,0:width,:] = frame
    full_frame[height:height+plot_height,0:width,:] = plot_force.getFrame()
    cv.imshow('frame',full_frame)
    out.write(full_frame)
    key = cv.waitKey()
    if key == ord('q'):
        break
    success,frame = video.read()
    t += techant

out.release()

           

Voici la vidéo obtenue. La vitesse de déroulement est de 30 images par seconde alors que la prise de vue est à 4 images par secondes. Le mouvement est donc accéléré d'un facteur 30/4.

Voici la force en fonction de la distance entre l'aimant et la bobine :

import numpy as np
from matplotlib.pyplot import *
[temps,F] = np.loadtxt("force-4.txt",unpack=True,skiprows=1)
F *= 1e-3*9.81 # conversion en mN
F -= np.mean(F[0:30])
# déplacement de 100 mm à 100 pas par secondes (1000 pas pour 1 cm), soit 1 mm par seconde
d = 100-temps
figure(1)
plot(d,F)
xlabel('d (mm)')
ylabel('F (mN)')
grid()
           
fig1fig1.pdf

Le déplacement motorisé de l'aimant permet d'obtenir la force en fonction de la position de l'aimant mais il ne permet pas d'obtenir le flux magnétique car le déplacement avec le moteur pas-à-pas est beaucoup trop lent pour donner une f.é.m. mesurable.

4. Flux et force en fonction de la distance

L'objectif est d'obtenir la force et le flux en fonction de la distance. Il est aisé d'obtenir la force et le flux en fonction de temps au cours du déplacement. Pour obtenir la position de l'aimant, nous optons pour une analyse de l'enregistrement vidéo du mouvement. Nous utilisons deux aimants plats en ferrite superposés, maintenus par la force magnétique sur une plaque en fer, elle-même vissée sur la jauge de contrainte qui mesure la composante horizontale de la force. L'autre extrémité de la jauge est fixée sur un support mobile guidé en translation par deux rails. L'électronique est aussi fixée sur ce support. La position de l'aimant est repérée par une pastille autocollante rouge. Afin de pouvoir faire un étalonnage précis de la distance sur les images de la vidéo (surtout reproductible d'une vidéo à l'autre), nous avons placé deux réglettes en polyéthylène (d'épaisseur environ 2 mm), l'une juste devant la face d'entrée de la bobine, l'autre à 100 mm. Voici une image extraite d'une des vidéos, dans la position innitiale de l'aimant :

aimant-bobine

Sur cette photo, la jauge de contrainte n'est pas visible car elle est masquée par les aimants. La caméra est une Basler acA1300-200uc (capteur 1280x1024 pixels de 6.14x4.92 mm) avec un objectif de focale 16 mm. La région effectivement filmée fait 1024x600 pixels. L'axe optique est placé au milieu du trajet de l'aimant. L'incertitude sur la position de l'aimant peut être estimée à 0,5 mm (incertitude type), ce qui est largement plus grand que la taille correspondant à un pixel (environ 0,14 mm).

Dans la première série d'expériences, l'enregistrement de la force électromotrice aux bornes de la bobine est fait avec l'enregistrement vidéo, comme expliqué dans la partie 1 ci-dessus. Le signal de f.é.m. comporte 800 échantillons à 400 Hz, soit une durée de 20 s. L'enregistrement vidéo est démarré exactement à t=10 s lorsqu'on appuie sur un bouton poussoir, ce qui génère un front montant permettant de déclencher l'aquisition de la tension. La caméra est pilotée par 400 impulsions de période 25 ms, soit une durée de 10 s, correspondant aux instants t=10 à t=20 s. On doit bien évidemment effectuer le déplacement manuel de l'aimant pendant la prise de vue vidéo.

L'extraction du mouvement de la pastille rouge sur la vidéo est faite avec l'application ColorTracking (page de description à venir), qui permet de suivre un objet repéré par sa couleur. Cette application permet aussi de déterminer l'échelle en millimètre sur la vidéo. Elle fournit un fichier contenant le temps et la position du centre de la pastille. L'origine du repère est placée sur la face gauche de la réglette placée contre la face d'entrée de la bobine. L'extraction de la position est stoppée à un moment où la pastille rouge est encore entièrement visible; en conséquence, la dernière distance obtenue est positive (environ 5 mm). Ce qui importe n'est pas la distance absolue (qui est définie avec une origine arbitraire) mais la reproductibilité de la mesure de distance d'une expérience à l'autre.

L'étape suivante consiste à calculer le flux magnétique et à générer une vidéo contenant la courbe de flux en fonction de la distance et de la f.é.m. en fonction de la distance. Le script suivant fait ce traitement et enregistre un fichier texte contenant la position de l'aimant et le flux.

analyse-mouvement-flux.py
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
from cvplot import CvPlot

nom = "mouvement-87"
fechant = 40 
techant = 1/fechant # période d'échantillonnage de la vidéo
video = cv.VideoCapture("%s-40fps.mp4"%nom)
t,x,y = np.loadtxt("%s-xy.txt"%nom,unpack=True,skiprows=1)
n= 0
success,frame = video.read()
n += 1
cv.namedWindow("frame")
(height,width,c)=frame.shape
[t,e] = np.loadtxt("%s.txt"%nom,unpack=True)
full_width = width*2
full_height = 3*height
#e = -e #changement de signe de la fem
i = 0
t_trigger = 10.0
while t[i] < t_trigger: i+=1
te = t[i+1]-t[i] # période d'échantillonnage de la fem
print("te = ",te)
e -= e[0:i].mean()
print("i = %d"%i)
plt.figure()
plt.plot(t,e)
flux = [0]
fem = [0]
for j in range(1,len(e)-i):
    flux.append(flux[j-1]-te*0.5*(e[i+j]+e[i+j-1]))
    fem.append(e[i+j])
flux = np.array(flux)
fem = np.array(fem)
plt.figure()
temps = np.arange(len(flux))*te
plt.plot(temps,flux)
plt.plot(temps,fem)
plt.show()
step = int(techant/te)
print("step = %d"%step)
t = 0
color = (255,0,0)

flux_max = np.max(np.absolute(flux))*1.2
fem_max = np.max(np.absolute(fem))*1.2

plot_flux = CvPlot(full_width,height,150,50,50,100,[0,100,-flux_max,flux_max])
plot_fem = CvPlot(full_width,height,150,50,50,100,[0,100,-fem_max,fem_max])

fourcc = cv.VideoWriter_fourcc(*'MP4V')
out = cv.VideoWriter("%s.mp4"%nom, fourcc, 30.0, (full_width,full_height))

i=0 # indice pour le tableau du flux
position = []
while success and n<len(x)-1:
    full_frame = np.zeros((full_height,full_width,3),np.uint8)
    plot_flux.clear()
    plot_flux.drawFrame()
    plot_fem.clear()
    plot_fem.drawFrame()
    plot_flux.labels("x (mm)","Phi (Wb)",size=1)
    plot_flux.xticks(10,"%0.1f",size=1)
    plot_flux.yticks(2,"%0.2g",full=True,size=1)
    plot_fem.labels("x (mm)","fem (V)",size=1)
    plot_fem.xticks(10,"%0.1f",size=1)
    plot_fem.yticks(2,"%0.2g",full=True,size=1)
    if i>=1:
        if i>=len(fem): i = len(fem)-1
        dx = (x[n+1]-x[n])/step
        for p in range(step):
            position.append(x[n]+p*dx)
        plot_flux.drawCurve(position,flux[0:i],(0,0,255),linewidth=2)
        plot_fem.drawCurve(position,fem[0:i],(0,0,255),linewidth=2)
    
    full_frame[0:height,0:width,:] = frame
    full_frame[height:2*height,0:full_width,:] = plot_flux.getFrame()
    full_frame[2*height:3*height,0:full_width,:] = plot_fem.getFrame()
    i += step
    cv.imshow('frame',full_frame)
    out.write(full_frame)
    key = cv.waitKey()
    if key == ord('q'):
        break
    success,frame = video.read()
    n += 1
    t += techant
out.release()
position = np.array(position)
np.savetxt("%s-flux(x).txt"%nom,np.array([position,flux[0:len(position)]]).T,header="x (mm)\ Flux(Wb)")
[x,phi] = np.loadtxt("%s-flux(x).txt"%nom,unpack=True,skiprows=1)
plt.figure()
plt.plot(x,phi)
plt.xlabel("x (mm)")
plt.ylabel("flux (Wb)")
plt.grid()
plt.show()
                
                       

La vidéo ci-dessous montre le résultat :

Voici le flux en fonction de la distance pour 5 expériences réalisées dans les mêmes conditions :

 
fig2= figure(2)
for n in [85,86,87,89,90]:
    [x,flux] = np.loadtxt("mouvement-%d-flux(x).txt"%n,unpack=True,skiprows=1)
    plot(x,flux,label="Flux "+"%d"%n)
xlabel("x (mm)")
ylabel("flux (Wb)")
grid()
legend(loc='upper right')    
                         
fig2fig2.pdf

La seconde expérience consiste à mesurer la force sur l'aimant lorsqu'un courant traverse la bobine. L'enregistrement vidéo et celui de la force sont déclenchés simultanément. Pendant les trois premières secondes, le courant dans la bobine est nul, ce qui permettra un ajustement du zéro de la force. La durée de l'enregistrement est de 30 secondes. La fréquence d'échantillonnage de la force est de 10 Hz, celle de la vidéo de 40 Hz.

L'électronique associée à la jauge de déformation et les programmes Arduino et Python permettant d'enregistrer la force en fonction du temps sont décrits dans Jauge de déformation. La fréquence d'échantillonnage est 10 Hz. La détermination de l'échelle de conversion entre la valeur numérique (24 bits) et la valeur de la force en grammes a été déterminée au préalable en plaçant le dispositif verticalement de manière à placer un poids de 50 g sur la plaque fixée à la jauge. Cette valeur d'étalonnage est codée en dur dans le programme Arduino.

Le script suivant génère une vidéo contenant la force en fonction de la position et un fichier texte contenant la position et la force.

analyse-mouvement-force.py
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
from cvplot import CvPlot

nom = "mouvement-77"
force = "force-77.txt"
fechant = 40 
techant = 1/fechant # période d'échantillonnage de la vidéo
video = cv.VideoCapture("%s-40fps.mp4"%nom)
t,x,y = np.loadtxt("%s-xy.txt"%nom,unpack=True,skiprows=1) # position obtenue par analyse de la vidéo
n= 0 # indice d'image
success,frame = video.read()
n += 1
cv.namedWindow("frame")
(height,width,c)=frame.shape
full_width = width
full_height = 2*height
[temps,F] = np.loadtxt(force,unpack=True,skiprows=1)

plt.figure()
plt.plot(temps,F)
g = 9.81
F *= g*1e-6
F -= np.mean(F[0:10]) # réglage du zéro (début sans courant dans la bobine)
plt.xlabel("t (s)")
plt.ylabel("F (N)")
plt.grid()
plt.show()
te = temps[1]-temps[0] # période d'échantillonnage de la force (1/10 s)
step = techant/te # = 0.25
print("step = %f"%step)



Fmax = np.max(np.absolute(F))*1.2
plot_force = CvPlot(width,height,100,50,50,100,[0,100,0,Fmax])
fourcc = cv.VideoWriter_fourcc(*'MP4V')
out = cv.VideoWriter("%s.mp4"%nom, fourcc, 30.0, (full_width,full_height))


i = 0 # indice pour le tableau de force
force = []
t = 0

while success and n < len(x):
    full_frame = np.zeros((full_height,full_width,3),np.uint8)
    plot_force.clear()
    plot_force.drawFrame()
    plot_force.labels("x (mm)","F (N)",size=1)
    plot_force.xticks(10,"%0.1f",size=1)
    plot_force.yticks(1,"%0.2g",full=True,size=1)
    k = int(i)
    force.append(F[k]+(F[k+1]-F[k])*(i-k)) # interpolation linéaire
    plot_force.drawCurve(x[0:n],force,(0,0,255),linewidth=2)
    i = i+step
    full_frame[0:height,0:width,:] = frame
    full_frame[height:2*height,0:width,:] = plot_force.getFrame()
    cv.imshow('frame',full_frame)
    out.write(full_frame)
    key = cv.waitKey()
    if key == ord('q'):
        break
    success,frame = video.read()
    n += 1
    t += techant
out.release()
force = np.array(force)

np.savetxt("%s-F(x).txt"%nom,np.array([x[0:len(force)],force],dtype=object).T,header="x (mm)\t F(N)")
[x,force] = np.loadtxt("%s-F(x).txt"%nom,unpack=True,skiprows=1)
plt.figure()
plt.plot(x,force,'.')
plt.xlabel("x (mm)")
plt.ylabel("F (N)")
plt.grid()
plt.show()                
                          

Voici le résultat pour une expérience (le courant dans la bobine est de 1 A) :

Voici la force en fonction de la distance pour 4 expériences réalisées dans les mêmes conditions (courant de 1 A) :

figure(3)
for n in [77,78,79,80]:
    [x,force] = np.loadtxt("mouvement-%d-F(x).txt"%n,unpack=True,skiprows=1)
    Fx = -force #composante de la force sur l'axe x
    plot(x,Fx,label="Force "+"%d"%n)
xlabel("x (mm)")
ylabel("Fx (N)")
grid()
legend(loc='upper right')                 
                         
fig3fig3.pdf

La force est échantillonnée à 10 Hz alors que la position (extraite de la vidéo) est échantillonnée à 40 Hz. Nous avons augmenté la fréquence d'échantillonnage de la force au moyen d'une interpolation linéaire.

Nous voulons vérifier la relation (3), qui dans la cas présent s'écrit :

Fx=IdΦdx(6)

Pour la position la plus lointaine (x=100mm), la force et le flux sont négligeables. On peut donc effectuer une intégration depuis cette position :

xmaxxFx(x)dx=IΦ(x)(7)

Voici l'intégration numérique de la force :

 
#intégration de la force depuis x=100 mm
figure(2)
for n in [77,78,79,80]:
    [x,force] = np.loadtxt("mouvement-%d-F(x).txt"%n,unpack=True,skiprows=1)
    force = -force #la force est en fait attractive
    # attention : x n'est pas échantillonné régulièrement
    # attention : x est décroissant dans le tableau
    N = len(x)
    integrale = np.zeros(N)
    s = 0
    for i in range(1,N):
        s += force[i]*(x[i]-x[i-1])*1e-3
        integrale[i] = s
    I = 1.0
    plot(x,integrale/I,'k--',label="Force "+"%d"%n)
title("I = 1,0 A")
grid()
xlabel("x (mm)")
ylabel("integrale de Fx (N.m) / I (A)")
legend(loc='upper right')
grid()                        
                          
fig4fig4.pdf

Cette figure montre les intégrales de la force divisées par l'intensité du courant dans la bobine (en trait noir interrompu) et les flux expérimentaux (en couleur). Pour un courant de 1 A dans la bobine, l'accord entre les deux courbes est bon, avec un écart qui ne dépasse pas 10% .

Voici les résultats pour un courant de 2 A dans la bobine :

 
figure(4)
for n in [85,86,87,89,90]:
    [x,flux] = np.loadtxt("mouvement-%d-flux(x).txt"%n,unpack=True,skiprows=1)
    plot(x,flux,label="Flux "+"%d"%n)
legend(loc='upper right')   
for n in [81,82,83,84]:
    [x,force] = np.loadtxt("mouvement-%d-F(x).txt"%n,unpack=True,skiprows=1)
    force = -force #la force est en fait attractive
    # attention : x n'est pas échantillonné régulièrement
    # attention : x est décroissant dans le tableau
    N = len(x)
    integrale = np.zeros(N)
    s = 0
    for i in range(1,N):
        s += force[i]*(x[i]-x[i-1])*1e-3
        integrale[i] = s
    I = 2.0
    plot(x,integrale/I,'k--',label="Force "+"%d"%n)
title("I = 2,0 A")
xlabel("x (mm)")
ylabel("integrale de Fx (N.m) / I (A)")
legend(loc='upper right')
grid()                        
                          
fig5fig5.pdf

Ces résultats expérimentaux valident la relation (3), au moins pour un courant dans la bobine jusqu'à 2 A.

Voici les résultats d'expériences réalisées avec un noyau en fer inséré dans la bobine sur environ la moitié de sa longueur. Le courant dans la bobine est de 0,5 A.

 
figure()
for n in [91,92,93,94,95,96]:
    [x,flux] = np.loadtxt("mouvement-%d-flux(x).txt"%n,unpack=True,skiprows=1)
    plot(x,flux,label="Flux "+"%d"%n)
legend(loc='upper right')   
for n in [97,99,100,101]:
    [x,force] = np.loadtxt("mouvement-%d-F(x).txt"%n,unpack=True,skiprows=1)
    force = -force #la force est en fait attractive
    # attention : x n'est pas échantillonné régulièrement
    # attention : x est décroissant dans le tableau
    N = len(x)
    integrale = np.zeros(N)
    s = 0
    for i in range(1,N):
        s += force[i]*(x[i]-x[i-1])*1e-3
        integrale[i] = s
    I = 0.5
    plot(x,integrale/I,'k--',label="Force "+"%d"%n)
title("Avec noyau, I=0,50 A")
xlabel("x (mm)")
ylabel("integrale de Fx (N.m) / I (A)")
legend(loc='upper right')
grid()                        
                          
fig6fig6.pdf

Ces résultats confirment que la relation (3) n'est pas valable en présence d'un pièce en ferromagnétique doux dont la distance à l'aimant varie. La force calculée au moyen de cette relation sous-estime la force réelle. Ce résultat peut se comprendre comme la conséquence du fait que l'aimantation du noyau en fer doux n'est pas constante : elle augmente lorsque l'aimant s'approche. Autrement dit, si le noyau est modélisé par une bobine, le courant dans celle-ci n'est pas constant.

Voici les résultats avec le même noyau et un courant de 1 A :

 
figure()
for n in [91,92,93,94,95,96]:
    [x,flux] = np.loadtxt("mouvement-%d-flux(x).txt"%n,unpack=True,skiprows=1)
    plot(x,flux,label="Flux "+"%d"%n)
legend(loc='upper right')   
for n in [102,103,104,105]:
    [x,force] = np.loadtxt("mouvement-%d-F(x).txt"%n,unpack=True,skiprows=1)
    force = -force #la force est en fait attractive
    # attention : x n'est pas échantillonné régulièrement
    # attention : x est décroissant dans le tableau
    N = len(x)
    integrale = np.zeros(N)
    s = 0
    for i in range(1,N): 
        s += force[i]*(x[i]-x[i-1])*1e-3
        integrale[i] = s
    I = 1.0
    plot(x,integrale/I,'k--',label="Force "+"%d"%n)
title("Avec noyau, I=1,0 A")
xlabel("x (mm)")
ylabel("integrale de Fx (N.m) / I (A)")
legend(loc='upper right')
grid()                         
                          
fig7fig7.pdf

Voici les résultats avec le même noyau et un courant de 2 A :

 
figure()
for n in [91,92,93,94,95,96]:
    [x,flux] = np.loadtxt("mouvement-%d-flux(x).txt"%n,unpack=True,skiprows=1)
    plot(x,flux,label="Flux "+"%d"%n)
legend(loc='upper right')   
for n in [106,107,108,109,110]:
    [x,force] = np.loadtxt("mouvement-%d-F(x).txt"%n,unpack=True,skiprows=1)
    force = -force #la force est en fait attractive
    # attention : x n'est pas échantillonné régulièrement
    # attention : x est décroissant dans le tableau
    N = len(x)
    integrale = np.zeros(N)
    s = 0
    for i in range(1,N): 
        s += force[i]*(x[i]-x[i-1])*1e-3
        integrale[i] = s
    I = 2.0
    plot(x,integrale/I,'k--',label="Force "+"%d"%n)
title("Avec noyau, I=2,0 A")
xlabel("x (mm)")
ylabel("integrale de Fx (N.m) / I (A)")
legend(loc='upper right')
grid()                         
                          
fig8fig8.pdf

L'écart entre la courbe de force et celle du flux se réduit lorsque l'intensité du courant augmente. On peut comprendre ce résultat : lorsque le courant dans la bobine augmente, l'effet de l'approche de l'aimant sur l'aimantation du noyau devient relativement moins important.

Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.