Bonjour à tous,

Cela fait un petit moment que je n’ai pas rédigé d’article mais il est parfois difficile d’allier le travail, la vie de famille et le blog. Et vu que le mois de mai est synonyme de ponts et reliquats de congés, j’ai enfin un peu de temps…

Couplage de la foudre sur une ligne

Une fois que l’on a compris le fonctionnement de la MKME (voir le précédent article sur la méthodologie d’utilisation de la MKME), on s’aperçoit que les possibilités sont immenses pour se représenter n’importe quel problème, qu’il soit d’origine CEM, électrique ou même mécanique.
Mais là où la MKME prend toute son ampleur, c’est pour se représenter rapidement et simplement des interactions électromagnétiques tels que les effets de la foudre sur un câble par exemple. Nous allons observer les effets de la foudre sur une ligne seule puis sur une ligne enfermée dans une cavité métallique. Pour commencer, nous regardons les effets de la foudre sur une ligne nue de 1 m de hauteur et de longueur 5 m. La distance entre la foudre et la ligne est de 100 m. Ce système est représenté par le schéma suivant :

La maille de gauche représente le circuit de la foudre. Nous modéliserons l’arc foudre par une FEM E_f ayant une forme double-exponentielle d’amplitude maximum 900 MV. La capacité C_{ns} représentant la capacité entre les nuages et le sol sera fixée à 87,5 nF. La résistance R_{ns} sert simplement à limiter le courant que nous fixons à 4 k\Omega, ce qui nous fera un courant maximum d’environ 200 kA. La maille de droite représente une petite ligne adaptée (R_0 = R_L = 50 \,\Omega) avec E_i la tension induite par le champ magnétique de la foudre B_0. L_L représente l’inductance du fil. \beta représente l’interaction, en l’occurrence l’induction magnétique crée par la première maille sur la deuxième.\\

E_f est de la forme :

 E_f= A(e^{-\frac{t}{\tau_1}} - e^{-\frac{t}{\tau_2}})

Avec \tau_1 : la largeur à mi-hauteur de l’impulsion.

et \tau_2 : le temps de montée (10-90%) de l’impulsion.

Nous partons du principe que le temps de propagation est très inférieur au temps de montée de la foudre. La place la matrice impédance et le vecteur source sont les suivant :

 Z=\begin{bmatrix}  R_{ns}+1/C_{ns}p & 0 \\  0 & R_0+R_L+L_Lp  \end{bmatrix}

 E=\begin{bmatrix}  e_f & e_i  \end{bmatrix}

Afin de calculer les courants, il nous faut encore déterminer la tension e_i. Nous partons des équations suivantes avec Q, le retard dû à la propagation du champ entre les deux mailles, S la surface de couplage et R la distance entre la foudre et la ligne :

    \[ B_0(t) = \frac{\mu_0i_1(t)}{2\pi R} \]

    \[ e_i(t) = -S\frac{\mathrm{d}B_0(t-Q)}{\mathrm{d}t} \]

Et nous développons afin de déterminer i_2 en mettant à gauche les éléments connus :

    \[e_i(t) = (R_0+R_L)i_2(t)+L_L\frac{\mathrm{d}i_2(t)}{\mathrm{d}t} = (R_0+R_L)i_2(t)+L_L\frac{\big(i_2(t)-i_2(t-1)\big)}{\mathrm{d} t} \]

    \[ -S\frac{\big( B_0(t-Q)-B_0(t-1-Q)\big)}{\mathrm{d} t} + \frac{L_L}{\mathrm{d} t}i_2(t-1) = (R_0+R_L+\frac{L_L}{\mathrm{d} t})i_2(t) \]

Nous calculons B_0(t-Q) en fonction des courants passés à l’aide d’une fonction « retard ». Le simulation (programme donnée en Annexe 1) donne les résultats suivants :

Blindage de la ligne

Après avoir calculé le courant sur une ligne nue, nous allons insérer cette même ligne, avec les même charges, dans une cavité métallique parallélépipédique de longueur L_b = 10 m, de hauteur et de profondeur h_b = l_b = 3 m et d’épaisseur de parois W_b = 2 mm. Nous la considérons en acier de conductivité \sigma=5,9.10^{-6}S.m^{-1}. Une cavité peut être comparée à un filtre passe-bas sur des champs EM. Une cavité possède une fréquence de coupure qui lui est propre en fonction de sa taille et en-deçà de laquelle il n’y a aucune atténuation :

    \[ f_0 = \frac{1}{2\pi}\frac{1}{\mu_0.\sigma.W_b.\frac{S_1}{L_b}} = \frac{1}{2\pi}\frac{R_h}{L_h} = 30 Hz \]

Avec R_h = \frac{L}{\sigma .W_b.h_b} : la résistance rencontrée par les courants de surface.

et L_h = \frac{\mu_0.S}{h} : l’inductance de la cavité.

et S : la surface sur laquelle se couple le champ EM.

Maintenant que nous avons définis notre boite, nous la représentons dans le schéma suivant :

Il s’agit du même schéma que précédemment mais avec la cavité insérée entre les deux autres mailles.

Nous commençons par poser les matrices suivantes :

Z_b=\begin{bmatrix}  R_{ns} & 0 & 0 & 0 & 0 & 0 & 0\\  0 & C_{ns} & 0 & 0 & 0 & 0 & 0\\  0 & 0 & R_b & 0 & 0 & 0 & 0 \\  0 & 0 & 0 & L_b & 0 & 0 & 0 \\  0 & 0 & 0 & 0 & R_0 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & L_L & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & R_L \\  \end{bmatrix}

C=\begin{bmatrix}  1 & 0 & 0 \\  1 & 0 & 0 \\  0 & 1 & 0 \\  0 & 1 & 0 \\  0 & 0 & 1 \\  0 & 0 & 1 \\  0 & 0 & 1 \\  \end{bmatrix}

Z=\begin{bmatrix}  R_{ns}+1/C_{ns}p & 0 & 0 \\  0 & R_b+L_b & 0 \\  0 & 0 & R_0+R_L+L_b \\  \end{bmatrix}

Nous décrivons les interactions \alpha et \beta du système :

    \[ \alpha \to e_b(t) = -S_1\frac{\mathrm{d}B_0(t-Q)}{\mathrm{d}t} = -h_b.L_b\frac{\mathrm{d}B_0(t-Q)}{\mathrm{d}t}\]

    \[\beta \to e_i = -S_2\frac{\mathrm{d}B_t(t)}{\mathrm{d}t} = -h_l.L_L\frac{\mathrm{d}B_t(t)}{\mathrm{d}t} \]

    \[B_t(t)=B_0(t)-B_r(y_0,z_0,t)\]

Avec B_t : le champ sous la ligne.

et B_r : le champ créé par le courant induit sur la boite.

Nous devons donc calculer B_r en tout point à l’intérieur de la cavité à un instant donné. L’intégration dans le temps se fera dans un second temps. Nous schématisons la cavité selon la figure suivante :

D’après la loi de Biot et Savart, nous pouvons dire que le champ élémentaire créé par l’élément \overrightarrow{\mathrm{d}l} (centré en P) au point M est le suivant :

    \[\mathrm{d}Br(z_0,y_0) = \frac{\mu_0}{4\pi} \frac{i\overrightarrow{\mathrm{d}l} \wedge \overrightarrow{u_{MP}}}{MP^2} \]

    \[\overrightarrow{u_{MP}} = \cos\phi.\overrightarrow{u_y} - \sin\phi.\overrightarrow{u_z} \]

\overrightarrow{\mathrm{d}l} peut être orienté en \pm\overrightarrow{u_y} et en \pm\overrightarrow{u_z}, on peut donc écrire :

    \begin{align*} \left\{ \begin{array}{cl} \mathrm{d}Br(z_0,y_0) = \frac{\mu_0}{4\pi} \frac{i \cos\phi}{R^2} \mathrm{d}z $\hspace{1 cm} pour \hspace{0.1 cm}$ z \, \mathcal{2} \, [0,h_b] \\ \mathrm{d}Br(z_0,y_0) = \frac{\mu_0}{4\pi} \frac{i \sin\phi}{R^2} \mathrm{d}y $\hspace{1 cm} pour \hspace{0.1 cm}$ y \, \mathcal{2} \, [0,L_b] \end{array}\right. \end{align*}

Nous intégrons sur \mathrm{d}z et \mathrm{d}y :

    \begin{align*} Br(z_0,y_0) = \frac {\mu_0} {4\pi} \left[ \int_{z(0)=0}^{h_b} \frac{\cos\phi} {R^2} \mathrm{dz} + \int_{z(L_b)=0}^{h_b} \frac{\cos\phi} {R^2}\mathrm{dz} + \int_{y(0)=0}^{L_b} \frac{\sin\phi} {R^2} \mathrm{dy} + \int_{z(h_b)=0}^{L_b} \frac{\sin\phi} {R^2}\mathrm{dy} \right] \end{align*}

    \begin{align*} Br(z_0,y_0) = \frac {\mu_0}{4\pi} \int_{z(0)=0}^{h_b} \frac{y_0} {(y_0^2+(z-z_0)^2)^{1/2} (y_0^2+(z-z_0)^2)^2} \mathrm{dz} \\ + \frac {\mu_0}{4\pi} \int_{z(L_b)=0}^{h_b} \frac{L_b-y_0} {(y_0^2+(z-z_0)^2)^{1/2} (y_0^2+(z-z_0)^2)^2} \mathrm{dz} \\ + \frac {\mu_0}{4\pi} \int_{y(0)=0}^{L_b} \frac{z_0} {(z_0^2+(y-y_0)^2)^{1/2} (z_0^2+(y-y_0)^2)^2} \mathrm{dy} \\ + \frac {\mu_0}{4\pi} \int_{z(h_b)=0}^{L_b} \frac{h_b-z_0} {(z_0^2+(y-y_0)^2)^{1/2} (z_0^2+(y-y_0)^2)^2} \mathrm{dy} \end{align*}

    \begin{align*} Br(z_0,y_0) = \frac {\mu_0}{4\pi} \int_{z(0)=0}^{h_b} (y_0^2+(z-z_0)^2)^{-3/2}) y_0. \mathrm{dz} + \frac {\mu_0}{4\pi} \int_{z(L_b)=0}^{h_b} (y_0^2+(z-z_0)^2)^{-3/2}) (L_b-y_0) \mathrm{dz} \\ + \frac {\mu_0}{4\pi} \int_{y(0)=0}^{L_b} (z_0^2+(y-y_0)^2)^{-3/2}) z_0. \mathrm{dy} + \frac {\mu_0}{4\pi} \int_{z(h_b)=0}^{L_b} (z_0^2+(y-y_0)^2)^{-3/2}) (h_b-z_0) \mathrm{dy} \end{align*}

Et nous établissons les équations afin de déterminer i_1, i_2 et i_3 en mettant à gauche les éléments connus :

    \[e_3(t) = -S_2 \frac{\mathrm{d}}{\mathrm{dt}} (B_0(t) - B_r(z_0,y_0,i_2(t),\mathrm{dy},\mathrm{dz})) \]

sur y \, \mathcal{2} \, [A,B]

et z \, \mathcal{2} \, [0,C]

    \begin{align*} \left\{ \begin{array}{cl} e_f(t) = (R_{ns}+ \frac{1}{C_{ns}} \int \mathrm{dt})i_1(t) \\ e_b(t) = (R_b+L_b\frac{\mathrm{d}}{\mathrm{dt}})i_2(t) \\ e_i(t) = (R_0+R_L+L_L\frac{\mathrm{d}}{\mathrm{dt}})i_3(t) + \beta' \frac{\mathrm{d}}{\mathrm{dt}}i_2(t) \end{array}\right. \end{align*}

    \begin{align*} \left\{ \begin{array}{cl} e_f(t) - \frac{\delta t}{C_{ns}}Q[t-1] = (R_{ns}+ \frac{1\delta t}{C_{ns}})i_1(t) \\ -S_1(B_0(t-Q)-B_0(t-1-Q)) + \frac{\mathrm{L_b}}{\mathrm{dt}}i_2(t-1) = (R_B + \frac{\mathrm{L_b}}{\mathrm{dt}})i_2(t) \\ -S_2 \frac{\beta'}{\delta t}i_2(t-1) + \frac{\mathrm{L_L}}{\mathrm{dt}}i_3(t-1) = (R_0 + R_L + \frac{\mathrm{L_L}}{\mathrm{dt}})i_3(t) - S_2 \frac{\beta'}{\delta t}i_2(t) \end{array}\right. \end{align*}

Nous programmons ces équations dans le script de l’annexe 2 qui nous permet d’observer les courbes suivantes :

Dans le même raisonnement que précédemment, nous observons une densité de courant de l’ordre de 1 A/m2 sur la cavité et un courant d’environ 20 μA sur la ligne. Nous observons donc très bien le phénomène de blindage autour de la ligne en comparant les courants avec et sans cavité. Nous calculons une atténuation d’environ 115 dB ce qui nous laisse penser que le résultat est toujours cohérent. Même si les valeurs des résultats sont discutable, on pourra toujours ajuster telle ou telle variable ou améliorer le calcul du champ sous la ligne par exemple. Mais dans l’ensemble cette simulation nous montre bien le phénomène d’atténuation et d’induction magnétique. Ici, nous traitons d’une boite parfaite mais on pourrait également dans une étape supplémentaire ajouter le couplage d’une ou plusieurs fente représentant une porte CEM par exemple. Cette étape sera peut-être traitée dans un futur article…

Voilà pour aujourd’hui, en espérant que vous avez apprécié cette petite démonstration et que la partie théorique ne vous a pas trop barbé.

A bientôt.

Annexe 1

Couplage d'un champ magnétique foudre sur une ligne de petite dimension 
 
@author: Thomas Raynaud
"""
" 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

N = 10000 # Nombre de points de calculs
dt = 10E-8 # Pas temporel
A = 900E6 # Amplitude
tm = 1.54E-6 # Temps de montée
td = 88E-6 # Temps de descente
R = 1E2 # Distance entre la ligne et la foudre
L = 5. # Longueur de la ligne
H = 1. # Hauteur de la ligne
S = H*L # Surface de couplage de la ligne
Ll = L*1E-6 # Indutance de la ligne
R0 = 50. # Résistance du générateur
Rl = 50. # Résistance de charge
Rns= 4000. # Résistance nuage-sol
Cns= 87.5E-9 # Capacité nuage-sol
mu0 = np.pi*4E-7 # Perméabilité du vide
q = int(R/3E8/dt) # Retard

tt = np.zeros((N,1),dtype=float) # Temps
I = np.zeros((N,2),dtype=float) # Courants
W = np.zeros((N,2),dtype=float) # Sources internes
E = np.zeros((N,2),dtype=float) # Sources foudre
Q = np.zeros((N,1),dtype=float) # Somme des courants dans C
B = np.zeros((N,1),dtype=float)
dB = np.zeros((N,1),dtype=float)

for t in range(1,N):
  tt[t] = t*dt
 
  Z = [[Rns+dt/Cns,0],
      [0,R0+Rl+Ll/dt]] 
 
  E[t,0] = A*(np.exp(-t*dt/td)-np.exp(-t*dt/tm))
  E[t,1] = S*dB[ret(t,q)]
 
 
  W[t,0] = -dt/Cns*Q[t-1,0]
  W[t,1] = Ll/dt*I[t-1,1]
  T = E[t,:]+W[t,:] # Vecteur sources complet
 
  y = np.linalg.inv(Z)
  I[t,:] = np.dot(T,y)

  Q[t,0] = Q[(t-1),0]+I[t,0]
  B[t] = mu0*I[t,0]/(2*R*np.pi)
  dB[t] = (B[t,0]-B[t-1,0])/dt

plt.subplot(2,1,1)
plt.plot(tt*1E6,I[:,0]/1E3)
plt.grid()
plt.title('Courant de foudre')
plt.xlabel('Temps (us)')
plt.ylabel('Courant (kA)')
plt.subplot(2,1,2)
plt.plot(tt*1E6,I[:,1])
plt.grid()
plt.title('Courant induit sur la ligne')
plt.xlabel('Temps (us)')
plt.ylabel('Courant (A)')
plt.show()

Annexe 2

Couplage d'un champ magnétique foudre sur une ligne de petite dimension 
 dans une boite faradisée
 
@author: Thomas Raynaud
"""
" Fonction ret : retard"
" ----------------------------------------"
def ret(n,m):
 if(n<m):
 return 0
 if(n>=m):
 return n-m 
" ----------------------------------------"

from scipy.integrate import quad 
import numpy as np
import pylab as plt

N = 10000              # Nombre de points de calculs
dt = 1E-8              # Pas temporel
u0 = np.pi*4E-7        # Permeabilité du vide
sigma= 5.9E6           # Conductivité acier carbone
A = 900E6              # Amplitude
tm = 1.54E-6           # Temps de montée
td = 88E-6             # Temps de descente
R1 = 1E2               # Distance entre la maille 1 et 2
R2 = 1.                # Distance entre la maille 2 et 3
Lb = 7.                # Longueur de la boite
lb = 3.                # Profondeur de la boite
hb = 3.                # Hauteur de la boite
Wb = 0.002             # Epaisseur du matériau
hl = 1.                # hauteur de la ligne
ll = 5.                # longueur de la ligne
S1 = 1.                # Surface de couplage de la boite
S2 = ll*hl             # Surface de couplage de la ligne
Ll = ll*1E-6           # Indutance de la ligne
R0 = 50.               # Résistance du générateur
Rl = 50.               # Résistance de charge
Rns= 4000.             # Résistance nuage-sol
Cns= 87.5E-9           # Capacité nuage-sol
Rh = 2*(Lb+lb)/(sigma*Wb) # Résistance de la paroi de la boite
Lh = u0*Lb*hb/lb       # Inductance de la boite
q = int(R1/3E8/dt)     # Retard entre maille 1 et 2
y0 = 1.                # Coordonnées point de départ
z0 = 1.                # Coordonnées point de départ
dy = 0.1               # Pas graphique suivant y
dz = 0.1               # Pas graphique suivant z


tt=np.zeros((N,1),dtype=float) # Temps
I=np.zeros((N,3),dtype=float)  # Courants
W=np.zeros((N,3),dtype=float)  # Sources internes
E=np.zeros((N,3),dtype=float)  # Sources externes
Q=np.zeros((N,1),dtype=float)  # Somme des courants dans C
B0=np.zeros((N,1),dtype=float)
Bt=np.zeros((N,1),dtype=float)
dB0=np.zeros((N,1),dtype=float)
Br=np.zeros((N,1),dtype=float)
dBt=np.zeros((N,1),dtype=float)
Br=np.zeros((N,1),dtype=float)

" Calcul de Br0 : champ normalisé crée par le courant circulant sur les parois de la boite"
" dans l'axe de la ligne à l'intérieur de la boite "
" ----------------------------------------"
def Bz0(z): 
 return (dz*y0*(y0**2+(z-z0)**2)**-1.5)
def BzLb(z): 
 return (dz*(Lb-y0)*((Lb-y0)**2+(z-z0)**2)**-1.5)
def By0(y): 
 return (dy*z0*(z0**2+(y-y0)**2)**-1.5)
def Byhb(y): 
 return (dy*(hb-z0)*((hb-z0)**2+(y-y0)**2)**-1.5)

Br1,err1 = quad(Bz0,0,hb)
Br2,err2 = quad(BzLb,0,hb)
Br3,err3 = quad(By0,0,Lb)
Br4,err4 = quad(Byhb,0,Lb)
Br0 = u0*(Br1+Br2+Br3+Br4)/(4*np.pi)
" ----------------------------------------"

for t in range(1,N):
  tt[t] = t*dt
  B0[t] = u0*I[ret(t,q),0]/(2*R1*np.pi) # Calcul du champ crée par la foudre
  dB0[t]= (B0[t]-B0[t-1])/dt            # Calcul de dérivée de B0
 
  Z1 = [[Rns+dt/Cns,0,0],               # Matrice impédance
        [0,Rh+Lh/dt,0],
        [0,0,R0+Rl+Ll/dt]] 
 
  E[t,0] = A*(np.exp(-t*dt/td)-np.exp(-t*dt/tm))
  E[t,1] = S1*dB0[t] 
  E[t,2] = S2*Bt[t-1]
 
  W[t,0] = -dt/Cns*Q[t-1]
  W[t,1] = Lh*I[t-1,1]/dt
  W[t,2] = Ll*I[t-1,2]/dt 
  T = E[t,:]+W[t,:]            # Vecteur sources complet
 
  y = np.linalg.inv(Z1)
  I[t,:] = np.dot(T,y)

  Q[t] = Q[t-1]+I[t,0]    # Calcul du courant dans Cns
  Br[t]= I[t,1]*Br0       # Calcul du champ crée par le courant I1
  Bt[t]= B0[t]-Br[t]      # Calcul du champ résultant dans la boite

plt.figure(1)
plt.subplot(3,1,1)
plt.plot(tt*1E3,I[:,0]/1E3)
plt.grid()
plt.title('Courant de foudre')
plt.ylabel('Courant (kA)')
plt.subplot(3,1,2)
plt.plot(tt*1E3,I[:,1]/(2*Lb*hb))
plt.grid()
plt.title('Densite de courant induit sur les parois de la boite')
plt.ylabel('Courant (A/m2)')
plt.subplot(3,1,3)
plt.plot(tt*1E3,I[:,2]*1E6)
plt.grid()
plt.title('Courant induit sur la ligne dans la boite')
plt.xlabel('Temps (ms)')
plt.ylabel('Courant (uA)')
plt.show()

Laisser un commentaire