Transducteur piézoélectrique pour émission-réception d'ultrasons

1. Introduction

Ce document montre comment mesurer l'impédance électrique d'un transducteur piézoélectrique de fréquence de résonance 40 kHz (émetteur-récepteur d'ultrasons). On verra comment déterminer les paramètres du modèle de Butterworth-Van Dyke, généralement utilisé pour représenter l'impédance au voisinage d'une résonance.

2. Mesure d'impédance

2.a. Montage expérimental

Soit Z le dipôle à étudier en fonction de la fréquence (Z désigne aussi son impédance complexe). L'enregistrement de la tension u(t) aux bornes de Z est fait avec la carte SysamSP5 en mode différentiel (entrées EA0 et EA4). Afin de mesurer le courant, on place en série avec Z une résistance R qui doit être de l'ordre de grandeur du module moyen de Z sur la plage de fréquence étudiée. L'enregistrement de la tension u_r(t) aux bornes de R est fait sur la voie EA1. La tension e(t) est sinusoïdale, de pulsation ω. Elle peut être délivrée par un générateur de fonctions (GBF) ou bien par la carte Sysam SP5 elle-même via la sortie SA1. Le courant maximal que la sortie d'un GBF peut fournir est 200 mA alors que celui de la carte Sysam est 50 mA. Dans le cas d'une très faible impédance, la résistance R pourra être nettement plus grande que |Z| afin d'éviter la saturation du courant de sortie du générateur. Par exemple pour la carte Sysam SP5, la tension maximale de la sortie SA1 étant +-10 V, on prendra R=200 Ω même si |Z| est de quelques ohms.

Figure pleine page

2.b. Traitement des signaux

L'intensité du courant dans le dipôle étudié est i(t)=u_r(t)/R. La moyenne de u(t) est calculée à partir du signal échantillonné par la moyenne arithmétique des N échantillons :

\bar{u} = \frac{1}{N} \sum_{n = 0}^{N - 1} u_{n} (1)

La fonction numpy.mean permet de calculer la valeur moyenne des échantillons stockés dans un tableau. Dans le cas cette présent, la tension u(t) a une moyenne nulle. L'amplitude efficace de u(t) est l'écart-type des N échantillons :

U_{e f f} = \sqrt{\frac{1}{N} \sum_{n = 0}^{N - 1} {(u_{n} - \bar{u})}^{2}} (2)

La fonction numpy.std() effectue le calcul de l'écart-type (standard deviation).

On obtient l'intensité du courant efficace I_eff de manière similaire et on en déduit l'impédance (module de l'impédance complexe) :

| Z | = \frac{U_{e f f}}{I_{e f f}} (3)

Considérons les expressions de l'intensité du courant et de la tension au borne de Z :

\begin{matrix} i (t) = I cos (ω t) & (4) \\ u (t) = | Z | I cos (ω t + φ) & (5) \end{matrix}

Afin de déterminer le déphasage $φ$ , considérons le produit suivant :

\begin{matrix} x (t) & = u (t) (i (t) - j i (t - \frac{T}{4})) & (6) \\ = | Z | I^{2} cos (ω t + φ) (cos (ω t) - j sin (ω t)) & (7) \\ = \frac{| Z | I^{2}}{2} (e^{j φ} + e^{- j (2 ω t + φ)}) & (8) \end{matrix}

Sa valeur moyenne est :

\bar{x} = \frac{| Z | I^{2}}{2} e^{j φ} (9)

Le déphasage $φ$ est donc l'argument de $\bar{x}$ .

Le calcul des échantillons x_n nécessite de disposer des échantillons de $i (t - \frac{T}{4})$ , c'est-à-dire i(t) en avance d'un quart de période. Il faut donc connaître le nombre d'échantillons correspondant à un quart de période, noté d, ce qui permet d'écrire :

x_{n} = u_{n} (i_{n} - j i_{n - d}) (10)

Si T_e est la période d'échantillonnage, on a en principe :

d = \frac{T}{4 T_{e}} (11)

Cette valeur n'est pas en général entière et il faut donc l'arrondir à l'entier le plus proche.

La détermination de d comporte deux difficultés :

Il faut connaître la période T. Si la tension e(t) est générée par la carte Sysam SP5 (comme expliqué plus loin), la période T est connue. Si la tension provient d'un générateur de fonctions numérique, la valeur précise de la fréquence est affichée sur l'appareil. Il est aussi possible de l'obtenir par analyse spectrale de u(t).
Lorsque le nombre d'échantillons par période est faible (de l'ordre d'une dizaine), l'échantillon d'indice n-d est très différent de l'échantillon décalé d'un quart de période, à cause de l'erreur d'arrondi sur la valeur de d.

La seconde difficulté peut être en grande partie résolue en augmentant la fréquence d'échantillonnage grace à une méthode d'interpolation. On procédera pour cela à une méthode d'interpolation par transformée de Fourier. Cette méthode consiste à calculer la transformée de Fourier discrète (TFD) des N échantillons du signal i_n, à séparer ses deux parties conjugées de taille N/2 puis à ajouter kN zéro entre ces deux parties. On obtient ainsi la TFD d'un signal qui comporte (k+1)N échantillons et qui est obtenu par la transformée de Fourier inverse. L'ajout de zéro dans le domaine fréquentiel est équivalent à une interpolation dans le domaine temporel (et réciproquement). La fréquence d'échantillonnage du nouveau signal est (k+1)f_e donc la période d'échantillonnage est T_e/(k+1). Il s'en suit que :

d = \frac{(k + 1) T}{4 T_{e}} (12)

Plus k est grand, plus la valeur de d est grande et moins elle est sujette aux erreurs d'arrondi.

2.c. Génération de la tension d'entrée

La tension est générée par la carte Sysam SP5, sur sa sortie SA1. Lorsqu'on utilise simultanément les entrées et une sortie, le convertisseur numérique-analogique de la carte fonctionne à la même fréquence d'échantillonnage que le convertisseur analogique-numérique. La fréquence d'échantillonnage maximale pour une utilisation simultanée des entrées et d'une sortie est 1  MHz. Nous pouvons donc en principe générer une sinusoïde jusqu'à une fréquence de 499 kHz. Pour des fréquences de l'ordre de 100 kHz, on pourrait s'inquiéter de la mauvaise qualité de la forme d'onde générée mais la nature statistique du traitement du signal détaillé ci-dessus le rend peu sensible aux inperfections de la forme d'onde, pourvu bien sûr que le critère de Nyquist-Shannon soit respecté.

Voyons comment sont générés les échantillons e_n de la tension e(t) pour une fréquence f souhaitée.

Le tableau contenant N échantillons (on prendra N=50000) doit correspondre à un nombre entier P de cycles du signal sinusoïdal. La fréquence est donc :

f = \frac{P}{N T_{e}} (13)

où T_e est la période d'échantillonnage. Le nombre d'échantillons par période est :

N_{p} = \frac{N}{P} (14)

Ce nombre est en général non entier, ce qui permet de faire varier finement la fréquence alors que la période d'échantillonnage ne peut prendre que des valeurs multiples d'une période de base, que nous fixons à 1 microseconde. Dans un premier temps, une valeur approximative de N_p est choisie, par exemple 100, mais cette valeur peut être abaissée lorsque la fréquence est trop grande.

Supposons que l'on ait fixé N et que l'on souhaite une valeur approximative Np et une fréquence freq. On commence par préciser la valeur de la période d'échantillonnage minimale :

teMin=1e-6

Lorsque le nombre Np d'échantillons par période est supérieur au rapport de la période par la période d'échantillonnage minimale, on doit corriger sa valeur :

Np = min(Np,1/(teMin*freq))

La période d'échantillonnage adoptée est un multiple de teMin qui permet d'obtenir environ la fréquence freq avec Np échantillons par période :

techant = int(1/(Np*freq*teMin))*teMin
if techant==0:
        techant = teMin

On calcule la valeur de P, le nombre de périodes pour les N échantillons :

P = int(freq*N*techant)

Pour finir, on calcule la fréquence effective, qui peut être légèrement différente de la fréquence demandée :

freq = P/(N*techant)

et la valeur effective du nombre de points par périodes :

Np=N/P

Cette valeur est différente de la valeur demandée car elle est en général non entière.

Voici comment se fait la génération des N échantillons de tension pour la sortie SA1 :

e1 = amp*np.cos(2*np.pi*P/N*np.arange(N))

Si l'on veut que la fréquence effective soit toujours très proche de la fréquence demandée, il faut que P soit grand et donc il faut que N soit très grand. La valeur N=50000 donne de très bons résultats. La quantité de mémoire interne de la carte Sysam SP5 qu'il faut utiliser est 3N mots de 12 bits, ce qui est largement inférieur à la valeur maximale (262144).

Pour l'obtention d'une courbe d'impédance en fonction de la fréquence, ce qui importe est la possibilité de faire varier finement la fréquence plus que la possibilité d'obtenir précisément une fréquence demandée.

2.d. Script Python

Ce script permet de faire une acquisition automatique de l'impédance pour une plage de fréquence choisie. La valeur de R doit être choisie afin qu'une amplitude de e(t) de l'ordre du volt conduise à une tension u_r(t) de l'ordre du volt également. La calibre des voies en entrée est fixé à 10 V. Compte tenu de la précision 12 bits du convertisseur A/N, l'échelon de tension est 4,9 mV et on peut donc sans problème accepter des amplitudes qui s'abaissent jusqu'à 100 mV.

analyseImpedance.py

import numpy as np
import matplotlib.pyplot as plt
import pycanum.main as pycan
import numpy.fft

La fonction interpol ci-dessous effectue l'interpolation par transformée de Fourier discrète. x est un tableau contenant le signal échantillonné. Le nombre d'échantillons est multiplié par ninter+1 (ninter contient le paramètre noté k ci-dessus).

def interpol(x,ninter):
    N = len(x)
    tfd = numpy.fft.fft(x)
    N1 = N//2
    tfd2 = np.concatenate((tfd[0:N1],np.zeros(N*ninter),tfd[N1:N]))
    y = np.real(numpy.fft.ifft(tfd2))*(ninter+1)
    return y

La fonction mesure présentée ci-dessous effectue l'acquisition et les calculs pour une fréquence donnée. Ses paramètres sont :

can : objet de la classe Sysam. Les deux entrées utilisées doivent être configurées au préalable.
R : valeur de la résistance.
freq : fréquence souhaitée, en hertz.
amp : amplitude de e(t), en volts.
delai : délai en secondes avant de commencer l'analyse des signaux (pour le régime transitoire).
N : nombre d'échantillons des signaux.
Np : nombre approximatif d'échantillons par période (maximal).
ninter : nombre d'échantillons ajoutés par interpolation.
plot : si True, les signaux e(t) et s(t) sont tracés.

Les valeurs renvoyées sont :

freq : fréquence effective du signal.
Z : impédance (module de l'impédance complexe).
phi : argument de l'impédance complexe, c.-à-d. déphasage entre i(t) et u(t).
techant : période d'échantillonnage, en secondes.
Np : nombre d'échantillons par période.
Ur : valeur efficace de u_r(t).

def mesure(can,R,freq,amp,delai,N=50000,Np=100,ninter=4,plot=False):
    teMin = 1e-6
    Np = min(Np,1/(teMin*freq))
    techant = int(1/(Np*freq*teMin))*teMin
    if techant==0:
        techant = teMin
    P = int(freq*N*techant)
    freq = P/(N*techant)
    Np = N/P
    e1 = amp*np.cos(2*np.pi*P/N*np.arange(N))
    can.config_echantillon(techant*1e6,N)
    can.acquerir_avec_sorties(e1,0)
    t=can.temps()[0]
    signaux = can.entrees()
    u=signaux[0]
    ur=signaux[1]
    u=u-u.mean()
    ur=ur-ur.mean()
    i=ur/R
    N = len(t)
    n1 = int(delai/techant)
    t=t[n1:N]
    t=t-t[0]
    u=u[n1:N]
    i=i[n1:N]
    N=len(u)
    r = ninter+1
    if ninter>0:
        u = interpol(u,ninter)
        i = interpol(i,ninter)
        t = np.arange(len(i))*t[len(t)-1]/len(i)
    d = int(Np*r/4)       
    Z = u.std()/i.std()
    print("d = %d"%d)
    phi = np.angle(np.mean(u[d:N]*(i[d:N]-1j*i[0:N-d])))
    if plot:
        plt.figure()
        plt.plot(t,u,'b')
        plt.plot(t,R*i,'r')
        plt.grid()
        plt.ylim(-10,10)
        plt.xlim(0,20/freq)
        plt.show()
    return(freq,Z,phi,techant,Np,ur.std())

Initialisation de l'interface, des listes et choix des fréquences :

can = pycan.Sysam("SP5")
can.config_entrees([0,1],[10.0,10.0],diff=[0])
liste_f = []
liste_Z = []
liste_phi = []
frequences = np.linspace(5e3,50e3,500)

choix des paramètres et de l'amplitude initiale :

R=1000 # résistance en série
delai = 2e-4 # délai du régime transitoire, à choisir en fonction du dipôle étudié
ninter=4
amp = 2.0

Nous avons fait le choix de fixer le gain des amplificateurs des convertisseurs A/N (calibre fixé à 10 V) mais d'ajuster l'amplitude de la tension e(t) afin que l'amplitude de la tension u_r(t) soit la plus grande possible et toujours inférieure à 9,5 V. Il est bon de contrôler la valeur efficace de la tension u_r(t) pendant l'acquisition. Si elle descend en dessous de 100 mV, il est préférable de reprendre le balayage avec une résistance R plus grande. Voici la boucle de traitement réalisant les mesures :

for f in frequences:    
    (f,Z,phi,te,Np,Ur) = mesure(can,R,f,amp,delai,ninter=ninter)
    while Ur*np.sqrt(2) > 9.5:
        amp *= 0.9
        (f,Z,phi,te,Np,Ur) = mesure(can,R,f,amp,delai,ninter=ninter)
    while Ur*np.sqrt(2) < 1.0 and amp < 9.5:
        amp *= 1.1
        (f,Z,phi,te,Np,Ur) = mesure(can,R,f,amp,delai,ninter=ninter)
    print("f = %f Hz, Z = %f, phi = %f, te = %f, Np = %f, Ur = %f"%(f,Z,phi,te,Np,Ur))
    liste_f.append(f)
    liste_Z.append(Z)
    liste_phi.append(phi)

can.fermer()
liste_phi = np.unwrap(liste_phi)
liste_Z = np.array(liste_Z)
liste_f = np.array(liste_f)

Les données sont sauvegardées :

np.savetxt('impedance-1.txt',np.array([liste_f,liste_Z,liste_phi]).T,header='f\t Z\t phi')

puis on trace le module et l'argument de l'impédance en fonction de la fréquence. Le tracé des parties réelle et imaginaire de l'impédance complexe est aussi intéressant.

plt.figure()
plt.plot(liste_f,liste_Z,'b-')
plt.grid()
plt.xlabel('f (Hz)')
plt.ylabel(r'$|Z|\ (\rm\Omega)$')
plt.figure()
plt.plot(liste_f,liste_phi/np.pi,'b-')
plt.xlabel('f (Hz)')
plt.ylabel(r'$\phi/\pi (\rm rad)$')
plt.grid()
plt.figure()
plt.plot(liste_f,liste_Z*np.cos(liste_phi),'b-')
plt.grid()
plt.xlabel('f (Hz)')
plt.ylabel(r'$Re(Z)\ (\rm\Omega)$')
plt.figure()
plt.plot(liste_f,liste_Z*np.sin(liste_phi),'b-')
plt.xlabel('f (Hz)')
plt.ylabel(r'$Im(Z)\ (\rm\Omega)$')
plt.grid()
plt.show()

2.e. Exemple

Le dipôle étudié est une bobine réalisée par enroulement de 13 spires sur un tore en ferrite. L'autoinductance totale est d'environ 1000 micro-henry. Nous souhaitons faire varier la fréquence entre 5 kHz et 50 kHz. L'impédance à 20 kHz devrait être d'environ 100 Ω. Nous choisissons donc cette valeur pour R.

On trace la partie réelle et la partie imaginaire de l'impédance complexe en fonction de la fréquence :

import numpy as np
from matplotlib.pyplot import *
from scipy.stats import linregress

[freq,Z_exp,phi_exp] = np.loadtxt('bobine-1.txt',unpack=True,skiprows=1) 
figure()
subplot(211)
plot(freq,Z_exp*np.cos(phi_exp),'k')
ylabel(r"$Re(Z)\ (\rm\Omega)$")
ylim(0,2)
grid()
subplot(212)
plot(freq,Z_exp*np.sin(phi_exp),'k')
xlabel("f (Hz)")
ylabel(r"$Im(Z)\ (\rm\Omega)$")
xlim(0,50e3)
ylim(0,300)
grid()

impedanceBobine.pdf

Si la bobine est modélisée par une impédance Z=r+jLω, la partie réelle est égale à r. Une régression linéaire sur la partie imaginaire permet d'obtenir le coefficient d'auto-inductance :

a, b, r_value, p_value, std_err = linregress(freq,Z_exp*np.sin(phi_exp))
L = a/(2*np.pi)

print((a,b))
--> (0.004687781673432709, 1.1386617514242516)

print(L)
--> 0.0007460836254624127

3. Étude du transducteur

3.a. Acquisition de l'impédance

La résistance choisie est R=1000 Ω. Un premier balayage est effectué avec 800 fréquences réparties linéairement de 20 kHz à 70 kHz. Un second balayage comporte 800 fréquences réparties de 70 kHz à 300 kHz. Voici les courbes obtenues :

[freq1,Z_exp1,phi_exp1] = np.loadtxt('piezo-1.txt',unpack=True,skiprows=1)
[freq2,Z_exp2,phi_exp2] = np.loadtxt('piezo-2.txt',unpack=True,skiprows=1)

figure(figsize=(16,8))
subplot(211)
plot(freq1,Z_exp1,'k')
plot(freq2,Z_exp2,'k')
ylabel(r"$|Z|\ (\rm\Omega)$")
grid()
subplot(212)
plot(freq1,phi_exp1/np.pi,'k')
plot(freq2,phi_exp2/np.pi,'k')
xlabel("f (Hz)")
ylabel(r"$\varphi/\pi\ (\rm rad)$")
grid()

impedancePiezoUS.pdf

Le premier minimum de |Z| (fréquence de résonance) se produit à f_r=40,2 kHz. C'est à cette fréquence qu'il faut exciter le transducteur pour obtenir une émission d'ultrasons car elle coïncide avec la fréquence de résonance mécanique. Il est suivi d'un maximum (anti-résonance) à f_ar=41,6 kHz. La deuxième résonance est à f_r'=50,0 kHz, suivie d'une anti-résonance à f_ar'=52,1 kHz. À beaucoup plus haute fréquence, on observe une résonance à 240 kHz suivie d'une anti-résonance à 253 kHz mais la forme de la courbe de phase montre qu'elle est de nature différente que les deux premières.

3.b. Modèle de Butterworth-Van Dyke

Ce modèle représente l'impédance électrique au voisinage d'une résonance et d'une antirésonance :

Figure pleine page

Le transducteur est constitué d'une lame d'un matériau piézoélectrique serrée entre deux électrodes métalliques. La capacité C_p est celle du condensateur constitué de ces deux électrodes. Elle est souvent donnée par le fabricant et se mesure facilement. Elle explique la tendance générale décroissante avec la fréquence de l'impédance. Pour notre transducteur, nous avons C_p=2,0 nF. La force appliquée à la lame est proportionnelle à la tension V. Dans l'effet piézoélectrique direct, l'application d'une contrainte engendre une tension. Dans l'effet piézoélectrique inverse, l'application d'une tension engendre une contrainte. Ces deux effets permettent au même transducteur d'être à la fois récepteur et émetteur d'ondes ultrasonores (ici à 40 kHz).

Les paramètres R_m,L_m,C_m proviennent du couplage électro-mécanique. L_m est proportionnelle à la masse équivalente mise en mouvement, C_m est proportionnelle au module de Young du matériau et R_m représente les pertes dissipatives et l'émission d'onde dans le milieu.

En l'absence de pertes (R_m=0) l'impédance s'écrit :

Z (R_{m} = 0) = \frac{- j}{C_{p} ω} \frac{\frac{1}{L_{m} C_{m}} - ω^{2}}{\frac{C_{p} + C_{m}}{L_{m} C_{m} C_{p}} - ω^{2}} (15)

La pulsation de résonance, qui annule l'impédance sans pertes, est aussi qualifiée de pulsation de résonance série car elle correspond à la condition de résonance de L_m et C_m en série. La pulsation de résonance pour R_m=0 est :

ω_{r} = \frac{1}{\sqrt{L_{m} C_{m}}} (16)

Cette fréquence est très proche de la fréquence de résonance mécanique. À cette fréquence, le piézoélectrique est équivalent à C_p et R_m en parallèle et la puissance transmise au transducteur lorsqu'il est utilisé en émetteur est :

P = \frac{V_{e f f}^{2}}{R_{m}} (17)

Cette puissance est essentiellement transmise sous forme d'onde acoustique mais une partie est dissipée dans l'émetteur.

La pulsation qui rend l'impédance sans pertes infinie est la pulsation d'anti-résonance, aussi appelée pulsation de résonance parallèle car elle fait intervenir C_p en parallèle avec C_m et L_m :

ω_{a r} = \sqrt{\frac{C_{p} + C_{m}}{L_{m} C_{m} C_{p}}} (18)

La simulation suivante montre l'influence de R_m sur la courbe représentant |Z| en fonction de la fréquence :

def impedance(Cp,Cm,Lm,Rm,f):
    w=2*np.pi*f
    Z1 = 1/(1j*Cp*w)
    Z2 = Rm+1j*Lm*w+1/(1j*Cm*w)
    return  1/(1/Z1+1/Z2)
    
Cp=2.1e-9
Cm=1.6e-10
Lm=0.10
f = np.linspace(38e3,43e3,1000)

figure()
Rm=200
plot(f,np.absolute(impedance(Cp,Cm,Lm,Rm,f)),label=r'$R_m=200\,\rm\Omega$')
Rm=400
plot(f,np.absolute(impedance(Cp,Cm,Lm,Rm,f)),label=r'$R_m=400\,\rm\Omega$')
Rm=600
plot(f,np.absolute(impedance(Cp,Cm,Lm,Rm,f)),label=r'$R_m=600\,\rm\Omega$')
Rm=800
plot(f,np.absolute(impedance(Cp,Cm,Lm,Rm,f)),label=r'$R_m=800\,\rm\Omega$')
grid() 
legend(loc='upper right')
xlabel('f (Hz)')
ylabel(r'$|Z|\ \Omega$')

courbesZ.pdf

La valeur de R_m a une grande influence sur la valeur maximale de l'impédance (à l'antirésonance) et dans une moindre mesure sur la valeur minimale (à la résonance). Une augmentation de R_m augmente légèrement la fréquence d'antirésonance et réduit légèrement la fréquence de résonance. La pulsation donnée par $(16)$ est donc légèrement plus grande que la pulsation de résonance réelle. La pulsation donnée par $(18)$ est légèrement plus petite que la pulsation d'antirésonance réelle.

3.c. Ajustement du modèle

La valeur de C_p étant connue, il s'agit de déterminer C_m,L_m,R_m. On commence par relever les fréquences de résonance et d'antirésonance sur la courbe expérimentale, f_r et f_ar. D'après les relations $(16)$ et $(18)$ , valables si R_m=0, on a :

\begin{matrix} C_{m} = ((\frac{ω_{a r}}{ω_{r}})^{2} - 1) C_{p} \\ L_{m} = \frac{1}{C_{m} ω_{r}^{2}} \end{matrix}

On se sert de ces deux relations pour déterminer C_m et L_m puis on ajuste R_m afin que la valeur maximale de l'impédance soit la même que la valeur expérimentale. Après ce premier réglage, il y un décalage de fréquence important entre la résonance du modèle et la résonance expérimentale (de même pour l'antirésonance), car la valeur de R_m a une influence notable sur ces fréquences. On est donc amené à augmenter légèrement la valeur de f_r par rapport à la valeur expérimentale et à réduire légèrement celle de f_ar. Après plusieurs essais et corrections, on parvient à un ajustement satisfaisant.

Pour le transducteur que nous avons étudié, la présence de deux couples résonance-antirésonance dans une plage de fréquence étroite nous oblige à considérer le modèle suivant pour représenter l'impédance entre 20 kHz et 70 kHz :

Figure pleine page

L'application de la méthode décrite ci-dessus conduit au résultat suivant :

Cp=2.0e-9
Rm=520
Rm2 = 520
fr=40.37e3
far=41.75e3
fr2=50.1e3
far2=51.75e3
Cm=((far/fr)**2-1)*Cp
Lm = 1/(Cm*4*np.pi**2*fr**2)
Cm2=((far2/fr2)**2-1)*Cp
Lm2 = 1/(Cm2*4*np.pi**2*fr2**2)
f = np.linspace(20e3,70e3,1000)
w=2*np.pi*f
Z1 = 1/(1j*Cp*w)
Z2 = Rm+1j*Lm*w+1/(1j*Cm*w)
Z3 = Rm2+1j*Lm2*w+1/(1j*Cm2*w)
Z = 1/(1/Z1+1/Z2+1/Z3)
figure(figsize=(16,6))
subplot(211)
plot(f,np.absolute(Z),'r',label='Modèle')
plot(freq1,Z_exp1,'k',label='Expérience')
xlabel("f (Hz)")
ylabel(r"$|Z|\ (\rm\Omega)$")
grid()
legend(loc='upper right')
subplot(212)
plot(f,np.angle(Z)/np.pi,'r')
plot(freq1,phi_exp1/np.pi,'k')
xlabel("f (Hz)")
ylabel(r"$\varphi/\pi\ (\rm rad)$")
grid()

ajustementModeleBVD2.pdf

print((Cm,Lm,Cm2,Lm2))
--> (1.3907226375005655e-10, 0.11175893139142654, 1.3390584101258618e-10, 0.0753642341250036)