Références

Classical electrodynamics

John David Jackson

Wiley Computational Electromagnetics

Anders Bondeson, Thomas Rylander, Pär Ingelström

Springer

E-book (accessible sur le campus) :

http://dx.doi.org/10.1007/b136922 Ce cours sur le web : http://bruno.lepetit.pagesperso-orange.fr/teaching.html

Une question ? [email protected]

1

modélisation multiphysique à l’échelle méso-macroscopique Plan

1. Introduction : que veut-on calculer ?

2. Rappels d’électromagnétisme 3. Une méthode de résolution des équations de Maxwell : Différences finies domaine temporel 4. Equation de la chaleur : traitement des non-linéarités 5. Applications industrielles Mécanique : Stéphane Guinard, EADS Electromagnétisme : Ivan Revel, EADS 2

Les matériaux : métalliques et composites… 787 Dreamliner Composite Profile

Les composites : fibres et résines Délaminages

Composite ; un empilement complexe : plis, protection de surface, peinture

Un matériau doit assurer une multitude de fonctions : tenue mécanique (efforts, chocs), thermique, protection électromagnétique (champs forts, foudre effets directs et indirects, électrostatique, furtivité…), tenue à la corrosion, au feu, esthétique…

La modélisation aide à dimensionner le matériau en amont dans le développement d’un programme, soit en amont d’essais, soit pour économiser des essais.

Plusieurs types de contraintes se prêtent bien à une modélisation de type physique, Passant par la résolution d’équations aux dérivées partielles : Mécanique, thermique, électromagnétique…

Modélisation électromagnétique

Que cherche-t-on ?

Effet d’une agression extérieure (champ fort, foudre) sur les systèmes embarqués (câblages, équipements…).

Une modélisation à 2 niveaux :

1. Modélisation locale

Ex : matériau et efficacité de blindage , câblages, fentes…

2. Modélisation globale

Des outils utilisant les équation de Maxwell permettent de réaliser ces modèles.

Modélisation thermique

Tenue en température de matériaux Tenue à la foudre (effets directs) de matériaux : électro-thermo-mécanique

EQUATIONS DE MAXWELL

E, H : champs D, B : inductions

Conservation de la charge : Milieux linéaires : permittivité électrique perméabilité magnétique Loi d’Ohm : 11

Liens entre les différentes équations de Maxwell

Je prends la divergence de l’équation de Faraday :  .

  

E

   

t

 .

B

 0 Par ailleurs je prends la divergence de l’équation d’Ampère :  .

  

H

  .

J

  

t

 .

D

 0 Par comparaison avec l’équation de conservation de la charge :  

t

 .

D

   

t

Conclusion : si à l’instant initial les champs satisfont :  .

D

  Et  .

B

 0 Alors à tout instant suivant, il suffit que les champs satisfassent les équations d’Ampère et de Faraday (et de conservation de la charge) pour satisfaire les autres équations de Maxwell.



On ne propage numériquement que les équations de Faraday et Ampère.

12

Propagation dans un milieu conducteur

Expulsion de la charge :

 .

E

     

t

    0    0

e

 

t

 0       10  9 36  1 5 .

10 7  0 .

2

as

Equations de propagation :

   

E

   

H



t

  

H

  

E

   

E



t

Solution recherchée sous la forme : 

E



H

  

E



H

0

e

0

i e i

( 

t

( 

t

  

k r

) 

k

 

r

) 13



i k

  

E



i k

  

H

 

i

 

H



i

 ( 1  

i

 ) 

E k

  (

k

  

E

)    2 ( 1  

i

 ) 

E k

  

k n k

 (  ) 2 1  1  

i

 1 2 (relation de dispersion) Vitesse de phase : Indice du milieu : Impédance d’onde :

v

 

k

 1 (  ) 2 1 1  

i

 1 2  1 (  0  0 ) 1 2

n Z

  ( 

r

 

E



H r

 ) 2 1 1 

n



i

 

E

 

H

  1 2

k

  1  1 

i

 1 ( 

r



r

) 2 1 1  

i

 1 2    1 2 1 

c n

1  

i

 1 2 14

Cas particulier : mauvais conducteur :    1

v

 

k Z

  1   1 2   1 2 Vide (SI):  0  10  9 36   0

Z

  4  10  7 120   377  Cas particulier : bon conducteur :

k

 (  2 ) 2 1 ( 1 

i

)  1  

i

   1

δ

: épaisseur de peau (décroit quand ω croit) :     2    1 2

Z

 (  2  ) 2 1 ( 1 

i

)  1  

i



R



iL

 15

Propagation avec atténuation : 

E



H

 

E

0

e i

( 

t



kz

)  

H

0

e i

( 

t



kz

)  

E

0

e i

( 

t

 

z

)

e

 

z

 

H

0

e i

( 

t

 

z

)

e

 

z

Exemples pour f=100 kHz (foudre…) : - Métal : σ=5 10 7 S/m 

δ≈0.2 mm

- Composite Carbone (CFC, CFRP…) σ= 10 4 S/m 

δ≈15 mm

16

Interfaces

Discontinuités aux interfaces

1 2

n

ρ s

(C/m 2 )

, J s (A/m): densités surfaciques

17

Discontinuités à la surface d’un conducteur Parfait Réel

- La composante normale de E peut connaitre une discontinuité pour un conducteur chargé.

-La composante tangentielle de H peut être discontinue s’il y a des courants de surface - Les autres composantes sont continues - Seule la composante normale de E peut être discontinue - Les autres décroissent avec l’effet de peau en pénétrant dans le conducteur.

18

Propagation d’un paquet d’ondes

u(x,t) : une composante de

E

ou

H

à une dimension Transformée de Fourier :

u

(

x

,

t

)  1 ( 2  ) 2 1    

A

(

k

)

e i

(  (

k

)

t



kx

)

dk ω(k)

: relation de dispersion

A

(

k

)  1 ( 2  ) 2 1    

u

(

x

,

t

 0 )

e ikx dx

Δx.Δk≈π

Onde plane infinie :

Δk=0

, Pulse infiniment étroit :

Δk=+∞

19

Dispersion

Supposons une loi de dispersion linéaire (vide, diélectrique…) :  (

k

)   (

k

0 ) 

d



dk

0 (

k



k

0 )

u

(

x

,

t

) 

e i

(  0 

d



dk

0

k

0 )

t

( 2  ) 2 1     

A

(

k

)

e



ik

(

x



d



dk

0

t

)

dk



e i

(  0 

d



dk

0

k

0 )

t u

(

x



d



dk

0

t

,

t

 0 ) Le paquet d’onde se propage sans déformation à la vitesse de groupe :

v g



d



dk

0 Dans le cas général (conducteur...) la relation de dispersion n’est pas linéaire. Il y a donc déformation du paquet d’onde durant sa propagation. 20

Dispersion numérique

Dans un milieu diélectrique (dispersion linéaire) combinant avec   

H

   

E



t

j’obtiens (équation de Helmholtz)   2 

t

 

E

2 

E

 

c

2    2  

H

 

E t



z

2 Question : Quelle relation de dispersion ?

Grille spatiale et temporelle utilisée pour la résolution numérique 21

L’équation de Helmholtz, Pour une composante de champ, en différences finies, devient : Le champ à l’instant n+1 est donc donné par la récurrence : Je cherche la relation de dispersion que doit satisfaire une onde du type

E



E

0

e j

( 

t



kz

) Pour être solution de l ’équation discrétisée

Nouvelle relation de dispersion :

QUESTION : Limite pour des intervalles tendant vers 0 ?

22

y

 2

R Arc

sin  

R

sin

x

2  

R



c



t



z y

  

z c

R=1.2

R=1.

R=0.9

R=0.01

x



k



z R=1 :

on retrouve la loi de dispersion physique :

ω=ck R>1

: pour certains intervalles en

k, ω



e iωt

devient complexe a une composante exponentielle réelle  amplification non physique de cette contribution. Instabilité

R<1

: dispersion non linéaire sur les

k

élevés  déformation du signal 23

CRENEAU

A

(

k

)  1 ( 2  ) 2 1   

u

 (

x

,

t

 0 )

e ikx dx A

(

k

)  1 ( 2  ) 2 1   

a a e ikx dx

 2  1 2

a

sin(

ka

)

ka

24

PROPAGATION NUMERIQUE DES EQUATIONS DE MAXWELL CAS 1D

On ne propage que les eq. d’Ampère et Faraday 1D Elles ont les 2 la même structure : Dérivée spatiale=dérivée temporelle 

E x



t



H y



t

    1  1  

H y



z



E



z x

Premières idées : 1. Différences finies centrées 2. Utiliser la même grille pour E et H (celle déjà utilisée). 25

Différences finies, grille unique

Différences finies Au point

(r,n)

: 

E



t



H x y



t

    1  1  

H y



z



E



z x E x H y n

 1 

E x r r

2 

t n

 1  2 

t H y n

 1

r n

 1

r

    1  1 

H y E x r n

  1 2 

z H r y n

 1  2 

z E x n r

 1

n r

 1 Algorithme du type « saut de grenouille » (leapfrog) : je passe directement du temps n-1 au temps n+1 en passant par-dessus n.

Inconvénient de l’algorithme : précision réduite du fait que, en temps comme en distance, On fait des sauts de 2 pas.

En utilisant des grilles décalées (staggered) pour E et H, on peut réduire le saut à un seul pas !

26

Différences finies, grilles décalées

Grille

E x

Grille

H y

: (r,n) : (r+1/2,n+1/2) Au point

(r,n+1/2)

: Au point

(r+1/2,n)

: 

E x



t



H y



t

    1  1  

H y



z



E



z x

- Les sauts sont bien d’un seul pas - Les points où sont utilisées les équations n’appartiennent ni à la grille E, ni à la grille H - QUESTION : quelle relation de dispersion ?

27

PROPAGATION NUMERIQUE DES EQUATIONS DE MAXWELL CAS 3D : ALGORITHME DE YEE

(K.S. Yee, IEEE Trans. Ant. Propag., AP-14, 302-307, 1966) Composantes du champ

E : milieu des arêtes Temps entiers

E x E y E z

: points (p+1/2, q, r, n) : points (p, q+1/2, r, n) : points (p, q, r+1/2, n) Composantes du champ

H : milieu des faces Temps demi-entiers

H x H y H z

: points (p, q+1/2, r+1/2, n+1/2) : points (p+1/2, q, r+1/2, n+1/2) : points (p+1/2, q+1/2, r, n+1/2) 28

CALCUL DES CHAMPS H AU TEMPS n+1/2

  

H



t

   

E

Au point (p, q+1/2, r+1/2, n) Au point (p+1/2, q, r+1/2, n) Au point (p+1/2, q+1/2, r, n) 29

CALCUL DES CHAMPS E AU TEMPS n+1

  

E



t

   

H

Au point (p+1/2, q, r, n+1/2) Au point (p+1/2, q+1/2, r, n+1/2) Au point (p, q, r+1/2, n+1/2) 30

ALGORITHME DE YEE : RELATION DE DISPERSION 3D

Quelle relation

ω(k)

pour avoir une solution du type 

E



H

  

E



H

0

e

0

i e i

( 

t

( 

t

  

k r

)  

k r

 )

Expressions des opérateurs différentiels

Si

f(x)=e i(ωt-kx)

alors 

f



t

 

f



x



f f

(

t

(

t

, 

x



t

2  ,

x

)  

t



x

)  2 

x f f

(

t

  2 (

t

,

x



t

,

x

) 

x

) 2    

t

2

i

sin(  ) 2 

t

2

i

sin( 

x k



x

) 2

f

(

x

,

t

) 

f

(

x

,

t

)

d t f



d x

(

x

,

t

)

f

(

x

,

t

)  

t



F

  

F

 

d t



F d

  

F d

      

d d d x y z

    

d t d x

   

t

2

i

sin(  ) 2 

t

2

i

sin( 

x k x



x

) 2 ...

31

Si 

E



H

  

E

0 

H

0

e i e i

( 

t

( 

t



k

 

r

  

k r

) ) Alors les équations de Yee donnent :   

H



t

    

E

  

E



t

   

H



d

t



H

 

d

t



E

 

d

  

E



d

  

H

Combinant ces 2 équations :  

d

t

2 

E



d

  (

d

  

E

)

d

  (

d

  

E

) Donc 

d



d

  (

d

 donc  

E

)  

E d

 

d

 

d

 

E

 Par ailleurs

d

 2 

E

Résultat 

d

t

2 

d

 2 

d

x

2 

d

y

2 

d

z

2

d

  (

d

  

E

)

Relation de dispersion 3D

sin(  

t

) 2 2  (

c



t

) 2     sin(

k x



x

2 

x

)   2     sin(

k y



y



y

2 )    2    sin(

k z

 

z

2

z

)   2   32

Condition de stabilité :  (

k x

,

k y

,

k z

), (

c



t

) 2     sin(

k x



x

2 

x

)   2    sin( 

k y



y

) 

y

2  2       sin(

k z



z

2 

z

)   2    1 Qui est satisfaite si : (

c



t

) 2   1 

x

2    1 

y

  2  1 

z

2    1

Condition CFL (Courant-Friedrichs-Levy) :



t



c

   1 1 

x

2    1 

y

  2  1 

z

2   1 2 33

Choix de la taille de maille : - Description suffisamment précise des détails géométriques - Critère de longueur d’onde ≈ 10 points par longueur d’onde Choix du pas de temps : critère CFL Occupation mémoire : (2 pas de temps stockés)*(6 composantes)*Nx.Ny.Nz=12 Nx.Ny.Nz

Temps de calcul : proportionnel à 6 Nx.Ny.Nz.Nt

 Proportionnel à (fréquence)**4  Impossibilité pratique du calcul à très hautes fréquences : autres méthodes (asymptotiques : optique, théorie de rayons…) 34

Exemple :

on veut modéliser un avion. Taille : 50 m.

On considère un coup de foudre qui dure 100 μs.

On veut représenter de façon suffisamment fine les détails  taille de maille : 0.1 m Nombre de mailles ≈ 1000 3 =10 9 (en fait on a moins besoin pour la hauteur) soit 12*8=96 Go d’occupation mémoire Pas de temps max : Δx/3 1/2 c≈0.2 ns, soit Nt=500 000 Nombre d’opérations : de l’ordre de 6*500 000*10 9 =3 10 15 Temps de calcul sur un processeur à 30 Gflops : 10 5 s ≈ la journée Fréquence max traitée correctement dans ce calcul : longueur d’onde = (10 points par longueur d’onde)*0.1m=1m, soit 300 MHz.

35

Limitations du volume maillé

Il faut absorber l’onde en bords de volume. Si on mettait un conducteur (effet de peau…) il y aurait des réflexions très importantes.

On ajoute une couche absorbante de matériaux fictifs (conductivité magnétique !)  en bords de volume  

E

     

H t

    

H

  

E

  * 

H

 

E



t

Impédance d’onde : (cf page 14)

Z

 

E



H

   1 2 1  1 

i



i

  *  1 2 1 2 Pour une onde normale à la surface le coefficient de réflexion est :  

Z Z

0  

Z

0

Z

Si Z=Z 0 : pas de réflexion ! Adaptation d’impédance On obtient cette adaptation d’impédance si :       0 0   *    36

CALCUL DU COEFFICIENT DE REFLEXION A L’INTERFACE

Pas de discontinuité à l’interface (Courant réparti dans l’épaisseur)

Z

E H

i i

 

E

r

H

r

 

E

t

H

t H i

0

H



i Z



H r Z

0 

H

0

i H i



H H r t

 

H r Z

0 

ZH t H r H t



ZH t

MILIEU Z 0

E

i k i

H

i

E

r -k i

H

r

H H r i

 

Z



Z Z

2

Z

0 

Z

0 2

Z

0 0

H t H t

Donc :  

Z Z

 

Z

0

Z

0 MILIEU Z

E

t

H

t

37

k t

Cas général : onde incidente non normale

J.P. Bérenger, J. Comput. Phys. 114, 185 (1994) Séparer les composantes du champ parallèle à la couche absorbante en 2 morceaux (ici, la couche absorbante en perpendiculaire à z)      

E xy



t

 

H



y z



E xz



t



E yz



t



E yx



t

 

E t z

   

H



z



H



z x y

   

H



x



H



x z y

  

E xz

 

E yz



H



y x

          

H xy



t

 

E



y z E E x y



H



t



H xz yz



t



H yx

  

t H t z

   

E



z



E



z x

   

E y



x



E



x z



y

    *

H yz



E



y x

*

H xz



E xy



E yx



E xz



E yz H H x y



H xy



H yx



H xz



H yz

Quand

σ=σ*=0,

on retrouve les équations dans le milieu physique.

Seule les morceaux en z (cad correspondant à une propagation en z) sont modifiés. La propagation en x,y est inchangée.

38

Dans la pratique, il faut considérer des murs perpendiculaires à x, y ,z.

…et il faut augmenter progressivement les conductivités, par exemple : 2  (

z

)   0

L z

0  (

z



z

0 )  0 ,  (

z



z

0 

L

)   0 39

Plaques minces : cas 2D (pour simplifier)

L.K. Wu et L.T. Han, IEEE Transactions on Antennas and Propagation, vol. 40, 628, 1992 On ne veut pas mailler l’épaisseur de la plaque, cela conduirait à des maillages énormes !

On remplace la plaque par une résistance équivalente  (pas d’effet de peau = pas d’induction) modèle basse fréquence Pas de dépendance en z

Question :

S 1 , S 2 ?

40

Loi d’Ohm sur la plaque résistive : L e L Continuité de la composante tangentielle du champ et loi d’Ohm dans la plaque Discontinuité du champ H tangentiel aux interfaces en présence de courant de surface Plaque carrée :

R

 1 

L Le

 1 

e

= Impédance de surface Ce qui rend caduque eq. 5c + 41

Algorithme initial Algorithme modifié

sur la plaque

2R

En supposant une dépendance linéaire en temps autour de

n+1/2

42

Prise en compte de fils minces

R. Holland, L. Simpson, IEEE Transactions on Electromagnetic Compatibility, vol. 23, 88, 1981 - Une perturbation EM est susceptible d’induire des courants parasites sur les systèmes - Il faut donc coupler les équations de propagation des 6 composantes de champ à des équations sur les courants et charges induits sur les fils

Etablissement des équations sur I et Q

  

E

    

H



t



E r



z

 

E



r z

   

H

 

t E z

(

r

) 

E z

(

a

)    

t a



r H



d

   

z a



r E r d

 On a

E z (a)=0

(pourquoi ?) 43

Hypothèses sur les champs (basse fréquence, et localement)

H

 

I

2 

r

(Théorème d’Ampère en magnétostatique)

E r



Q

2 

r

(Théorème de Gauss en électrostatique)

E z

(

r

)    

t a



r H



d

   

z a



r E r d



E z

(

r

)  

r

ln( 2 

a

)   

I



t

 1  

Q



z

 

E z



L

  

I



t

 1  

Q



z

  À une distance moyenne R

Equation de continuité

dQ

(

z

)

dz

 

Q



t

  

I



z I

(

z

)

dt



I

(

z



dz

)

dt I(z) Q(z)dz I(z+dz)



Q



t

 

I



z

 0 44

Modification de l’algorithme de Yee avec fil (vertical)

Intensités

I :

milieu des arêtes Temps demi-entiers

I

: points (p, q, r+1/2, n+1/2)

Charge

Charges Q

:

coins du cube, temps entiers Q: points (p, q, r, n)

Intensité

- Par cellule : 8 quantités (6 composantes de champ+Q+I) à propager - Il faut donc coupler les équations de propagation des 6 composantes de champ à des équations sur les courants et charges induits sur les fils 45

Modification de l’algorithme de Yee avec fils

1. Calculer I au temps n+1/2, à partir de I au temps n-1/2 et de E et de Q au temps n

I n

 1 /

pqr

 2 1 / 2 

I n

 1 /

pqr

 2 1 / 2  

t L E z n pqr

 1 / 2 

L Q n pqr

 1   

z Q n pqr

2. Calculer H au temps n+1/2, à partir de H au temps n-1/2 et de E au temps n Yee sans changement 3. Calculer E au temps n+1, à partir de E au temps n et de H et de I au temps n+1/2 Yee avec rajout du courant du fil 

I n

 1 /

pqr

 2 1 / 2 

X

.



Y

4. Calculer Q au temps n+1, à partir de Q au temps n et de I au temps n+1/2 46

n

 1

Q pqr



Q n pqr

 

t



z



I n

 1

pqr

/  2 1 / 2 

I n

 1

pqr

/  2 1 / 2 