Méthodologie d’utilisation de la MKME

Cet article à pour but d’expliquer comment poser une stratégie pour analyser, synthétiser puis résoudre les problèmes de compatibilité électromagnétique en utilisant les principes de la MKME. Une explication de la transformation en graphe d’un schéma classique puis deux exemples simples de lignes de transmission seront traités.

Analyse du problème

En premier lieu, nous devons analyser ce que nous voulons observer. Pour illustrer notre article, nous allons prendre l’exemple d’une ligne de transmission quelconque raccordée à un générateur d’impédance interne R_0 d’une part, et à une charge quelconque Z_L d’autre part. Le schéma représentatif d’une ligne est le suivant :

Figure 1 – Schéma électrique d’une ligne

Une ligne de transmission est susceptible d’être perturbée par des sources de rayonnement externes au circuit :

  • des couplages électrostatiques,
  • des couplages magnétostatiques,
  • des couplages en champ proche,
  • des couplages en champ lointain.

Il est nécessaire de bien comprendre le fonctionnement d’une ligne avant d’intégrer le fonctionnement de ces perturbations. C’est pourquoi nous allons détailler, dans cet article, l’utilisation de la MKME sans prendre en compte les rayonnements extérieurs. Cela sera traité dans un prochain article.

Le but de notre étude est de trouver une méthode permettant d’observer le fonctionnement d’une ligne. Pour cela, nous allons représenter notre schéma sous forme d’un graphe de Branin qui va nous permettre d’analyser différemment ce circuit :

Figure 2 – Graphe de Branin d’une ligne

Dans ce graphe, nous représentons le retard temporel induit par la ligne par une interaction (corde) entre deux mailles. Chaque composant est schématisé par les branches 1 à 4 et la ligne est divisée en deux branches (2 et 3) associées à deux sources e_g et e_d. Ces sources permettent d’induire les effets d’une des maille sur l’autre. Les équations qui régissent le circuit sont les suivantes :

(1)    \begin{equation*} \left\{ \begin{array}{cl} e_g=(V_C+Z_Ci_2)e^{-\tau/p} \\ e_d=(V_A-Z_Ci_1)e^{-\tau/p} \end{array} \right. \end{equation*}

Mise en équation

Maintenant que nous avons notre schéma et nos équations, il ne nous reste plus qu’à développer afin de discrétiser notre système pour la programmation. Nous cherchons à observer les courants des deux mailles. Pour commencer, il nous faut déterminer les sources de chaque maille et les intégrer dans le covecteur E :

(2)    \begin{equation*} E= \begin{bmatrix}e_0 & e_0e^{-\tau/p} \end{bmatrix} \end{equation*}

Ensuite, il nous faut déterminer une métrique (une matrice d’impédances), Z_{nm} qui puisse transformer chaque fém en courant de la manière suivante :

(3)    \begin{equation*} e_n=Z_{nm}i^m \end{equation*}

Pour déterminer Z_{nm}, nous devons déterminer toutes les impédances du circuit. Nous effectuons une somme directe de ces éléments que nous intégrons dans une matrice que nous appellerons « métrique du réseau » :

(4)    \begin{equation*} Z_b=\begin{bmatrix} R_0 & 0 & 0 & 0 \\ 0 & Z_C & 0 & 0 \\ 0 & 0 & Z_C & 0 \\ 0 & 0 & 0 & Z_L \end{bmatrix} \end{equation*}

Cette métrique est définie dans l’espace des branches et pour calculer les courants de mailles, nous devons la définir dans l’espace des mailles. Pour cela nous utilisons une matrice de connectivité C qui va nous permettre d’associer ces impédances dans l’espace des mailles :

(5)    \begin{equation*} C=\begin{bmatrix} 1 & 0 \\ 1 & 0 \\ 0 & 1 \\ 0 & 1 \end{bmatrix} \end{equation*}

Chaque ligne représente une branche et chaque colonne une maille. Le chiffre « 1 » indique une connexion entre la branche et la maille donnée. En effectuant les produits matriciels entre la connectivité transposée, la matrice Z_b et la connectivité, nous obtenons la métrique dans l’espace des mailles :

(6)    \begin{equation*} Z=C^T*Z_b*C \end{equation*}

(7)    \begin{equation*} Z=\begin{bmatrix} R_0+Z_C & 0 \\ 0 & Z_C+Z_L \\ \end{bmatrix} \end{equation*}

Sachant que nous connaissons les fém, il faut inverser ce tenseur d’impédances par un tenseur d’admittances :

(8)    \begin{equation*} \begin{array}{cl} i^m=Z_{nm}^{-1}e_n\\ i^m=Y^{mn}e_n \end{array} \end{equation*}

Cas d’une charge résistive

Il nous faut, pour finir, déterminer les sources des mailles en fonction d’éléments connus. Dans un premier exemple, nous donnerons à Z_L une impédance résistive que nous renommons R_L :

(9)    \begin{equation*} \begin{array}{cl} e_0(t)-e_g(t)=(R_0+Z_C)i_1(t)\\ e_d(t)=(Z_C+R_L)i_2(t) \end{array} \end{equation*}

Nous arrangeons les équations du Branin (1) :

(10)    \begin{equation*} \begin{array}{cl}e_g(t)=(R_L-Z_C)i_2(t-\tau)\\ e_d(t)=e_0(t-\tau)+(Z_C-R_0)i_1(t-\tau) \end{array} \end{equation*}

Puis nous remplaçons et développons en mettant les termes connus à gauche. En effet, nous admettons pour la suite fonctionner en mode implicite, nous connaissons donc les termes en t-\tau :

(11)    \begin{equation*} \begin{array}{cl} e_0(t)-(Z_L-Z_C)i_2(t-\tau)=(R_0+Z_C)i_1(t)\\ e_0(t-\tau)+(Z_C-R_0)i_1(t-\tau)=(Z_C+R_L)i_2(t) \end{array} \end{equation*}

Nous avons donc deux équations à deux inconnues à résoudre à l’aide du logiciel Spyder.

 Programmation du cas résistif

Nous avons, à présent, toutes les données en main pour commencer la programmation Python. Le programme complet et commenté est donné en Annexe 1. Comme dans tout bon programme classique, nous déterminons en premier les fonctions appelées par la suite dans le programme ainsi que les composants et paramètres d’entrée. Dans notre cas, nous créons une fonction permettant de créer un retard ret :

" Fonction ret : retard"
" ----------------------------------------"
def ret(n,m):
    if(n<m):
        return 0
    if(n>=m):
        return n-m
" ----------------------------------------"
C=3E8               # Célérité
d=10.               # Longueur de la ligne
tau=d/C             # Vitesse de propagation
Zc=200.             # Impédance de la ligne
Ro=10E-3            # Résistance interne
Rl=1E6              # Résistance de charge
dt=2E-10            # Pas temporel
N=2000              # Nombre d’échantillons
tt=np.zeros((N,1),dtype=float)
J=np.zeros((N,2),dtype=float)
eo=np.zeros((N,1),dtype=float)

Puis viens ensuite la boucle de calcul qui va nous permettre, à partir d’une fonction de Heaviside ou gaussienne en e_0(t), de calculer les courants dans les différentes mailles. Nous reprenons les résultats des équations trouvées précédemment :

for t in range(1,N):
    s=t*dt
    Z=[[Ro+Zc,0],[0,Zc+Rl]]
    #eo[t]=np.exp(-((s-4e-9)/2e-9)**2.)
    eo[t]=1-np.exp((-s/2e-10)) 
    q=int(tau/dt) 
    eg=(Rl-Zc)*J[ret(t,q),1] 
    ed=eo[ret(t,q)]+(Zc-Ro)*J[ret(t,q),0] 
    E=[eo[t][0]-eg,ed] 
    tt[t]=s y=np.linalg.inv(Z) 
    J[t,:]=np.transpose(np.dot(y,np.transpose(E)))

Dans ce programme, nous écrivons les deux sources trouvées théoriquement e_g et e_d et, à la suite, nous indiquons le covecteur E contenant la somme de chacune des sources des mailles (selon leur sens, bien entendu). Tout cela donne les résultats suivants selon le type d’impulsion en e_0 :

Figure 3 – Résultats pour une impulsion Heaviside
Figure 4 – Résultats pour une impulsion gaussienne

Cas d’une charge RLC

Dans ce cas, nous suivons la même logique que pour le cas précédent mais le raisonnement est un peu plus complexe car la présence d’une capacité et d’une self-inductance viens ajouter deux mailles supplémentaires, ce qui change notre connectivité ainsi que notre tenseur d’impédances. Comme nous travaillons en temporel, ces éléments vont également ajouter des opérations en dérivées et en intégrales. Nous travaillerons donc avec le schéma suivant :

Figure 5 – Schéma d’une ligne avec une charge bouchon

De la même manière que précédemment, nous transformons en graphe de Branin :

Figure 6 – Schéma du Branin avec une charge bouchon

Nous commençons par calculer les matrices suivantes :

(12)    \begin{equation*} Z_b=\begin{bmatrix} R_0 & 0 & 0 & 0 & 0 & 0 \\ 0 & Z_C & 0 & 0 & 0 & 0 \\ 0 & 0 & Z_C & 0 & 0 & 0 \\ 0 & 0 & 0 & R_L & 0 & 0 \\ 0 & 0 & 0 & 0 & Lp & 0 \\ 0 & 0 & 0 & 0 & 0 & 1/Cp \\ \end{bmatrix} \end{equation*}

(13)    \begin{equation*} C=\begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \end{equation*}

(14)    \begin{equation*} Z=\begin{bmatrix} R_0+Z_C & 0 & 0 & 0 \\ 0 & Z_C+R_L & Z_C & Z_C \\ 0 & Z_C & Z_C+Lp & Z_C \\ 0 & Z_C & Z_C & Z_C+\frac{1}{Cp} \\ \end{bmatrix} \end{equation*}

 

La lettre p est l’opérateur de Laplace remplaçant l’expression j\omega. Reste donc à trouver les sources des mailles du circuit en fonction des éléments connus :

(15)    \begin{equation*} \left\{ \begin{array}{cl} e_0(p)-e_g(p)=(R_0+Z_C)i_1(p)\\ e_d(p)=(Z_C+R_L)i_2(p)+(Z_C+Lp)i_3(p)+(Z_C+\frac{1}{Cp})i_4(p) \end{array}\right. \end{equation*}

Nous arrangeons les équations du Branin (1) :

(16)    \begin{equation*} \left\{ \begin{array}{cl} e_g(t)=(R_L-Z_C)i_2(t-\tau)-Z_C[i_3(t-\tau)+i_4(t-\tau)]\\ e_d(t)=e_0(t-\tau)+(Z_C-R_0)i_1(t-\tau) \end{array}\right. \end{equation*}

Puis nous remplaçons et développons en mettant les termes connus à gauche. Nous devons également remplacer C_p et L_p par leurs équivalents temporels :

(17)    \begin{equation*} \left\{ \begin{array}{cl} e_0(t)-(R_L-Z_C)i_2(t-\tau)-Z_C[i_3(t-\tau)+i_4(t-\tau)]&=(R_0+Z_C)i_1(t)\\ e_0(t-\tau)+(Z_C-R_0)i_1(t-\tau)&=(Z_C+R_L)i_2(p)+(Z_C+Lp)i_3(p)+(Z_C+\frac{1}{Cp})i_4(p) \end{array}\right. \end{equation*}

(18)    \begin{equation*} \left\{ \begin{array}{cl} e_0(t)-(R_L-Z_C)i_2(t-\tau)-Z_C[i_3(t-\tau)+i_4(t-\tau)]&=(R_0+Z_C)i_1(t)\\ e_0(t-\tau)+(Z_C-R_0)i_1(t-\tau)&=R_Li_2(t)+Z_C[i_2(t)+i_3(t)+i_4(t)]+L\frac{\mathrm{d}i_3(t)}{\mathrm{d}t}+\frac{1}{C}\int_t i_4(t) \mathrm{d}t \end{array}\right. \end{equation*}

Nous pouvons discrétiser ces équations tout en mettant les termes connus à gauche. Nous utilisons la méthode de l’intégration discrète pour échantillonner l’intégrale. Cette méthode induit une erreur mais il s’agira de choisir un pas de temps \mathrm{d}t adéquat. Nous introduisons donc Q_4(t) =\sum_{0}^{t}i_4(t\,\mathrm{d}t), que nous utilisons de la manière suivante : Q(t-\tau) =\sum_{0}^{t-\tau}i_4(t\,\mathrm{d}t). Les équations deviennent donc :

(19)    \begin{equation*} \left\{ \begin{array}{cl} e_0(t)-(R_L-Z_C)i_2(t-\tau)-Z_C[i_3(t-\tau)+i_4(t-\tau)]=(R_0+Z_C)i_1(t)\\ e_0(t-\tau)+(Z_C-R_0)i_1(t-\tau)-\frac{L}{\delta}i_3(t-\tau)-\frac{\delta}{C}Q_4(t-\tau)\\ =(Z_C+R_L)i_2(t)+Z_C[i_3(t)+i_4(t)]+\frac{L}{\delta}i_3(t)+\frac{\delta}{C}i_4(t) \end{array}\right. \end{equation*}

De manière préalable à la partie programmation et à un instant t = N et un retard \tau = R, nous pouvons écrire ceci, en faisant apparaître les sources e_g et e_d retardées :

Programmation du cas RLC

Nous voyons que les expressions sont un peu plus compliquées mais il n’est beaucoup plus compliqué de programmer ce cas. Le programme complet et commenté est donné en Annexe 2. Nous énumérons les paramètres d’entrée comme précédemment :

" Fonction ret : retard"
" ----------------------------------------"
def ret(n,m):
    if(n<m):
        return 0
    if(n>=m):
        return n-m
" ----------------------------------------"
import numpy as np
import pylab as plt
R0=200.      # Résistance interne du générateur
Zc=200.      # Impédance de la ligne
Rl=1000.     # Résistance de charge
L=1E-6       # Self de charge
C=1E-9       # Capacité de charge
d=10.        # Longueur de la ligne
dt=0.6E-10   # Pas temporel
N=5000       # Nombre d’échantillons
tau=d/3E8    # Vitesse de propagation
q=int(tau/dt)# Calcul du retard
J=np.zeros((N,4),dtype=float)  # Tableau des courants
tt=np.zeros((N,1),dtype=float) # Tableau du temps
Q=np.zeros((N,1),dtype=float)  # Tableau somme des courants dans C
W=np.zeros((1,4),dtype=float)  # Tableau des sources internes
e0=np.zeros((N,1),dtype=float) # Tableau de la source e0
Jt=np.zeros((N,1),dtype=float) # Somme des courants des i2, i3, i4

Et nous mettons en place la boucle de calcul avec les fonctions Heaviside ou gaussienne, au choix. Il nous faut également ajouter les calculs des courants dans C que nous appelons Q[t] comme décrit dans la partie théorique (Q_4). Q permet d’additionner l’ensemble des courants déjà calculés précédemment et nous indique donc les courants dans C. Nous utilisons une méthode un peu différente que pour le cas précédent pour écrire les sources. Nous n’écrivons pas directement les sources e_g et e_d . En effet, pour simplifier la compréhension, nous avons séparé les sources externes, avec le covecteur E, des sources internes, avec le covecteur W. Les sources e_g et e_d sont intégrées dans W. Il n’y a qu’une seule source externe e_0 qui apparaît tel quel dans la maille 1, et qui est retardée dans les mailles 2, 3 et 4 comme illustré ci-dessous :

Les sources internes, correspondent aux éléments restant dans les équations théoriques. En dernier lieu, nous effectuons la somme de ces deux éléments dans le covecteur source complet T qui nous servira à calculer les courants J. Nous avons également ajouté l’élément J_T , la somme des courants i_2, i_3 et i_4, qui est finalement l’information utile que nous afficherons par la suite :

for t in range(N):
    e0[t]=np.exp(-((t*dt-4e-9)/2e-9)**2.)  # Impulsion gaussienne
    #e0[t]=1.-np.exp(-t*dt/50E-12)    # Impulsion échelon
for t in range(1,N):
    tt[t]=t*dt                          # Calcul du temps
    Z=[[R0+Zc,0   ,0      ,0      ],    # Matrice impédances
    [0    ,Zc+Rl,Zc    ,Zc     ],
    [0    ,Zc  ,Zc+L/dt,Zc     ],
    [0    ,Zc  ,Zc     ,Zc+dt/C]]
    E=[e0[t,0],e0[ret(t,q),0],e0[ret(t,q),0],e0[ret(t,q),0]]      
    W[0][0]=-Rl*J[ret(t,q),1]+Zc*(J[ret(t,q),1]+J[ret(t,q),2]+J[ret(t,q),3])
    W[0][1]=(Zc-R0)*J[ret(t,q),0]
    W[0][2]=(Zc-R0)*J[ret(t,q),0]+L/dt*J[(t-1),2]
    W[0][3]=(Zc-R0)*J[ret(t,q),0]-dt/C*Q[(t-1),0]

    T=E+W[0,:]                     # Covecteur source complet
    y=np.linalg.inv(Z)             # Inversion de la matrice Z    
    J[t,:]=np.dot(T,y)             # Calcul des courants de mailles
    Q[t]=Q[t-1]+J[t,3]             # Calcul du courant dans C
    Jt[t]=J[t,1]+J[t,2]+J[t,3]     # Somme des courants i2, i3, i4

Ce programme donne les résultats suivants selon le type d’impulsion en e_0 :

Figure 7 – Résultats pour une impulsion Heaviside
Figure 8 – Résultats pour une impulsion gaussienne

 Organigramme générique

Ces deux exemples relativement simples nous ont permis d’appréhender la résolution de problèmes CEM en utilisant la MKME. D’une manière générale, nous pouvons schématiser sous la forme d’un l’organigramme (voir figure 9) les étapes à suivre afin de résoudre n’importe quel problème. Bien entendu, ces étapes sont génériques et il faut adapter ses études aux différents problèmes rencontrés, mais cela permet d’avoir un plan adaptable.

Organigramme générique

Annexe 1

Modèle de Branin (équivalent ligne de transmission)
    avec charge résistive     

@author: Thomas Raynaud

import numpy as np
import pylab as plt

" Fonction ret : retard"
" ----------------------------------------"
def ret(n,m):
    if(n<m):
        return 0
    if(n>=m):
        return n-m       
" ----------------------------------------"
C=3E8               # Célérité
d=10.               # Longueur de la ligne
tau=d/C             # Vitesse de propagation 
Zc=200.             # Impédance de la ligne
Ro=10E-3            # Résistance interne du générateur
Rl=1E6              # Résistance de charge

dt=2E-10            # Pas temporel
N=2000              # Nombre d'échantillons

tt=np.zeros((N,1),dtype=float) # Tableau du temps
J=np.zeros((N,2),dtype=float)  # Tableau des courants
eo=np.zeros((N,1),dtype=float) # Tableau de la source e0

for t in range(1,N):
    s=t*dt                           # Calcul du temps
    Z=[[Ro+Zc,0],[0,Zc+Rl]]           # Matrice impédance (mailles)
    #eo[t]=np.exp(-((s-4e-9)/2e-9)**2.)     # Impulsion gaussienne
    eo[t]=1-np.exp((-s/2e-10))         # Ou impulsion échelon
    q=int(tau/dt)                     # Calcul du retard
    eg=(Rl-Zc)*J[ret(t,q),1]                # Mise en forme de eg
    ed=eo[ret(t,q)]+(Zc-Ro)*J[ret(t,q),0]   # Mise en forme de ed
    E=[eo[t][0]-eg,ed]                 # Covecteur source
    tt[t]=s                             # Vecteur temps
    y=np.linalg.inv(Z)               # Inversion de la matrice Z
    J[t,:]=np.transpose(np.dot(y,np.transpose(E)))  # Calcul des courants de mailles

Annexe 2

Modèle de Branin (équivalent ligne de transmission)
    avec charge RLC parallèle (Bouchon)
    
@author: TR

" Fonction ret : retard"
" ----------------------------------------"
def ret(n,m):
    if(n<m):
        return 0
    if(n>=m):
        return n-m       
" ---------------------------------------"
import numpy as np
import pylab as plt

R0=200.      # Résistance interne du générateur
Zc=200.      # Impédance de la ligne
Rl=1000.     # Résistance de charge
L=1E-6       # Self de charge
C=1E-9       # Capacité de charge
d=10.        # Longueur de la ligne
dt=0.6E-10   # Pas temporel
N=5000       # Nombre d'échantillons
tau=d/3E8    # Vitesse de propagation 
q=int(tau/dt)# Calcul du retard

J=np.zeros((N,4),dtype=float)  # Tableau des courants
tt=np.zeros((N,1),dtype=float) # Tableau du temps
Q=np.zeros((N,1),dtype=float)  # Tableau somme des courants dans C
W=np.zeros((1,4),dtype=float)  # Tableau des sources internes
e0=np.zeros((N,1),dtype=float) # Tableau de la source e0
Jt=np.zeros((N,1),dtype=float) # Somme des courants des i2, i3, i4

for t in range(N):
    e0[t]=np.exp(-((t*dt-4e-9)/2e-9)**2.)  # Impulsion gaussienne
    #e0[t]=1.-np.exp(-t*dt/50E-12)  # Impulsion échelon
    
for t in range(1,N):
    tt[t]=t*dt                          # Calcul du temps
    Z=[[R0+Zc,0   ,0      ,0      ],    # Matrice impédances
       [0    ,Zc+Rl,Zc     ,Zc     ],
       [0    ,Zc  ,Zc+L/dt,Zc     ],
       [0    ,Zc  ,Zc     ,Zc+dt/C]]
      
    E=[e0[t,0],e0[ret(t,q),0],e0[ret(t,q),0],e0[ret(t,q),0]]
    W[0][0]=-Rl*J[ret(t,q),1]+Zc*(J[ret(t,q),1]+J[ret(t,q),2]+J[ret(t,q),3])
    W[0][1]=(Zc-R0)*J[ret(t,q),0]
    W[0][2]=(Zc-R0)*J[ret(t,q),0]+L/dt*J[(t-1),2]
    W[0][3]=(Zc-R0)*J[ret(t,q),0]-dt/C*Q[(t-1),0]
    
    T=E+W[0,:]                     # Covecteur source complet
    y=np.linalg.inv(Z)             # Inversion de la matrice Z
    
    J[t,:]=np.dot(T,y)             # Calcul des courants de mailles
    Q[t]=Q[t-1]+J[t,3]             # Calcul du courant dans C
    Jt[t]=J[t,1]+J[t,2]+J[t,3]     # Somme des courants i2, i3, i4

plt.plot(tt,J[:,0])
plt.plot(tt,Jt)
plt.grid(True)
plt.title('Courants de mailles')
plt.xlabel('Temps (s)')
plt.ylabel('Courant (A)')

Laisser un commentaire