r,n+1/2 - bruno lepetit

Report
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  
 .B  0
t


Par ailleurs je prends la divergence de l’équation d’Ampère :  .  H   . J   . D  0
t
Par comparaison avec l’équation de conservation de la charge :
Conclusion : si à l’instant initial les champs satisfont :

t
 .D  
 .D 
Et

t
 .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 :

t
 .E 



0

   0e

 
t

0


 
10
9
1
36  5 . 10
7
 0 . 2 as

Equations de propagation :


H
  E  
t



E
  H  E  
t
Solution recherchée
sous la forme :

 i (  t  kr )
E  E0e

 i (  t  kr )
H  H 0e
13

 ik 

 ik 


E   i  H

 
 
2
k  ( k  E )    (1 
)E
i 

 
H  i  (1 
)E
i 


k  kn
k  (  )  1 
1

1
2
2
i 
(relation de dispersion)
Vitesse de phase :
v

k
1

1
(  ) 2 1 

Indice du milieu :
Impédance d’onde :

n  ( r r ) 1 

1
1
(  0 0 )
2
i 
1
2
1
1
2
1
( r r ) 2 1 

1
c

n
2
i 
1
2
i 


E
E
k
Z      

H
nH
1
1
 1

i 
  2
 
 
1
1

i 
1
2
14

Cas particulier : mauvais conducteur :
v


k

 1
1
0 
1
   2
Vide (SI):
1
  2
Z  
 
k (
1
) (1  i ) 
2
1 i
2
Z (

2
1
) (1  i ) 
2
9
36 
 0  4  10
7
Z  120   377 
Cas particulier : bon conducteur :

10

1 i



 1
1
 2 2
δ : épaisseur de peau   


(décroit quand ω croit) :
  
 R  iL 
15
Propagation avec atténuation :

 i (  t  kz )
 i ( t  z )  z

E  E0e
 E0e
e 

 i (  t  kz )
 i ( t  z )  z

H  H 0e
 H 0e
e 
Exemples pour f=100 kHz (foudre…) :
- Métal : σ=5 107 S/m  δ≈0.2 mm
- Composite Carbone (CFC, CFRP…) σ= 104 S/m  δ≈15 mm
16
Interfaces
Discontinuités aux interfaces
1
2
n
ρs (C/m2) , Js (A/m): densités surfaciques
17
Discontinuités à la surface d’un conducteur
Parfait
- 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
Réel
- 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

1


( 2 ) 2
A(k )e
i (  ( k ) t  kx )
dk
ω(k) : relation de dispersion
A(k ) 
1
1
( 2 )



u ( x , t  0 ) e dx
ikx
2
Δx.Δk≈π
Onde plane infinie : Δk=0 , Pulse infiniment étroit : Δk=+∞
19
Dispersion
Supposons une loi de dispersion linéaire
(vide, diélectrique…) :
i ( 0 
u ( x, t ) 
e
d
dk
k0 )t
0
1
( 2 ) 2

 ik ( x 


A(k )e
d
dk
t)
0
 (k )   (k0 ) 
dk  e
i ( 0 
d
dk
k0 )t
0
u(x 
d
dk
d
dk
(k  k0 )
0
t, t  0)
0
Le paquet d’onde se propage sans déformation
à la vitesse de groupe :
vg 
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


E
avec   H  
t
j’obtiens (équation de Helmholtz)


H
  E  
t

 E
2
t
2

 E
2
c
2
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  E0e
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

 x 
y
Arc sin  R sin   
R
 2 

R=1.
2
R 
c t
z
y
R=1.2
R=0.9
R=0.01
z
c
x  k z
R=1 : on retrouve la loi de dispersion physique : ω=ck
R>1 : pour certains intervalles en k, ω devient complexe
eiωt 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
1
( 2 ) 2



u ( x, t  0)e
ikx
dx
1
A(k ) 
1
1
( 2 )
2

a
a
 2  2 sin( ka )
ikx
e dx    a
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
t
y


1 H

y
z
1 E x
 z
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
E x
t
H
Différences finies
Au point (r,n) :
Ex
n 1
 Ex
r

r
2t
n 1
H
y
r
H
2t
n
1 H
y
r 1

r

y
z
1 E x
 z
H
n
y
r 1
2z
n 1
y


y
t
n 1
1 H

1 Ex

n
r 1
 Ex
n
r 1
2z
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 Ex : (r,n)
Grille Hy : (r+1/2,n+1/2)
E x
t
H
t
y


1 H

y
z
1 E x
 z
Au point (r,n+1/2) :
Au point (r+1/2,n) :
- 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
Ex : points (p+1/2, q, r, n)
Ey : points (p, q+1/2, r, n)
Ez : points (p, q, r+1/2, n)
Composantes du champ H : milieu des faces
Temps demi-entiers
Hx : points (p, q+1/2, r+1/2, n+1/2)
Hy : points (p+1/2, q, r+1/2, n+1/2)
Hz : 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

 i (  t  kr )
E  E0e

 i (  t  kr )
H  H 0e
Quelle relation ω(k) pour avoir une solution du type
Expressions des opérateurs différentiels
f
Si f(x)=e i(ωt-kx) alors
t
f
x

 
F  dtF
t
 

F  dF
f (t 

t
2
t
, x)
2
x
)  f (t , x 
2
x
2
x
dx
 
d  dy
d
 z





t
2 i sin(
f ( x, t )  d t f ( x, t )
t
k x
 2 i sin(
)
)
2

t
f (t , x 

, x )  f (t 
)
2

f ( x, t )  d x f ( x, t )
x
2 i sin(
dt 
t
2
t
 2 i sin(
dx 
...
)
k xx
)
2
x
31
Si

 i (  t  kr )
E  E0e

 i (  t  kr )
H  H 0e

Alors les équations
de Yee donnent :
Résultat
Relation de dispersion 3D
t

E

t
 


   H  d t E  d  H

 

  d E  d  ( d  E )

 

donc
d  (d  E )  d

 
Donc
d  (d  E ) 

 d  d
 


   E   d t H   d  E
2
t
Combinant ces 2 équations :
2
t

H
2


Par ailleurs
E  d
  
2 
d dE  d E

 
d  (d  E )
 dx  dy  dz
2
2
2
2
2
2

k yy 

k xx 
k z  z  


2
 sin(
)
)
)
 sin(
 sin(
t 

2
2
2
2
 
 
 
)   (c t ) 
 sin(



2 
x
y
z



 


 

 



 




32
Condition de stabilité :
2
2
2

k

y


k

x
k z  z  


y
  sin( x
 sin(
)
)
)
 sin(
2
2
2
2
 
 
 1
 ( k x , k y , k z ), ( c  t ) 



x
y
z


 

 

 



 





  1 2  1 2  1 2 
2
  
  
 1
Qui est satisfaite si : ( c  t )  
  x 
  z  
 y 

Condition CFL (Courant-Friedrichs-Levy) :
t 
1
1
2 2
  1 2  1 

1


  
c 
  
 
  x 
  z  
 y 

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 ≈ 10003=109 (en fait on a moins besoin pour la hauteur)
soit 12*8=96 Go d’occupation mémoire
Pas de temps max : Δx/31/2c≈0.2 ns, soit Nt=500 000
Nombre d’opérations : de l’ordre de 6*500 000*109=3 1015
Temps de calcul sur un processeur à 30 Gflops : 105 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



*
H
*

  E  
 H

1 1
Impédance
t
E
i 
  2
d’onde :
Z    



(cf page 14)
H
 
E

  H  E  
1
t
i 
 
Pour une onde normale à la surface le coefficient de réflexion est :
1
2
1
2
Z  Z0
Z0  Z
Si Z=Z0 : 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
MILIEU Z0
Pas de discontinuité à l’interface
(Courant réparti dans l’épaisseur)
Ei  Er  Et
-ki
kt
Ht
t
Hr
Z 0 H i  Z 0 H r  ZH
Hi  Hr  Ht
 
Et
Er
Hi  Hr  Ht
Donc :
ki
Hi
Hi  Hr  Ht
Z 0 H i  Z 0 H r  ZH
Ei
MILIEU Z
Hi 
t
Hr 
Z  Z0
2Z 0
Z  Z0
2Z 0
Ht
Ht
Z  Z0
Z  Z0
37
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
 E xz
t
 E yz
t
 E yx
t
E z
t

H z
y
H


y
z
H x
  E xz
  E yz
z



H z
H
x



x
y

H x
y

 H xy
t
 H xz
t
H
yz
t
H
yx
t
H z
t

E y  E yx  E yz
H
y
E y
z
E x
z


H x  H xy  H xz
y
 H
yx
H
yz
E z


E x  E xy  E xz
  H xz
*
 H
*
E z
x
E y
x

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.
E x
y
38
Dans la pratique, il faut considérer des murs perpendiculaires à x, y ,z.
…et il faut augmenter progressivement les conductivités, par exemple :
 z  z0 
 (z)   0

L


2
 ( 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  modèle basse fréquence
(pas d’effet de peau = pas d’induction)
Pas de
dépendance
en z
Question :
S1 , S 2 ?
40
Loi d’Ohm sur la plaque résistive :
L
Plaque carrée :
R 
e
1 L
 Le

1
e
= Impédance de surface
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
-
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


H
  E  
t
E r
z
E z (r )  E z (a )  

r

E z
r
 

t
r
H  d 
E


t
z
a
H 
r
d
a
On a Ez(a)=0 (pourquoi ?) 43
Hypothèses sur les champs
(basse fréquence, et localement)
H 
(Théorème d’Ampère en magnétostatique)
I
2 r
(Théorème de Gauss en électrostatique)
Er 
E z (r )  

r

t
H  d 
a
2  r
r
E

z
r
d
a
ln(
E z (r )  

Q
r
 I
1 Q 

E z  L 

  t   z 
)
a   I  1  Q 


2    t   z 
À une distance moyenne R
Equation de continuité
Q
dQ ( z ) dz  I ( z ) dt  I ( z  dz ) dt
Q
t

I
I(z)
Q(z)dz
I(z+dz)
t

I
z
0
z
44
Modification de l’algorithme de Yee avec fil (vertical)
Intensité
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)
- 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 / 2
pqr  1 / 2
t 
E

z
L 
n
pqr  1 / 2
n
n
Q pqr  1  Q pqr 

L

  z

Yee sans changement
2. Calculer H au temps n+1/2, à partir de H au temps n-1/2
et de E au temps n
3. Calculer E au temps n+1, à partir de E au temps n
et de H et de I au temps n+1/2
 I
n 1 / 2
pqr  1 / 2
Yee avec rajout du courant du fil
n 1 / 2

I pqr  1 / 2
 X . Y
46
4. Calculer Q au temps n+1, à partir de Q au temps n
et de I au temps n+1/2
n 1
Q pqr  Q pqr 
n
t
z
I
n 1 / 2
pqr  1 / 2
n 1 / 2
 I pqr 1 / 2


similar documents