Semi Final Project.docx

Rapport de Projet de Méthode Numérique Master Spécialisé Telecommunications Systems Engineering Sujet :

Méthode FDTD ; Stabilité et Convergence.

SOMMAIRE FDTD 1. 2. 3. 4. 5. 6.

Introduction, Technique de FDTD Algorithme FDTD Convergence de la méthode FDTD Stabilité de la méthode FDTD 1D FDTD sur Matlab

1. FDTD : 1.1 Introduction,

FDTD est l'acronyme de l'expression anglaise Finite Difference Time Domain. C'est une méthode de calcul de différences finies dans le domaine temporel, qui permet de résoudre des équations différentielles dépendantes du temps. Cette méthode est couramment utilisée en électromagnétisme pour résoudre les équations de Maxwell. Elle était appliquée avec succès dans des différents problèmes , comme la diffusion à partir des objets métalliques et des diélectriques , des antennes , des microstriplines et aussi l’absorption ElectroMagnétique dans le corps humain exposé à des radiations . La raison principale pour le succès de cette méthode est que la méthode elle même est la simplicité d’implémentation . Elle a été premièrement proposé par K.Yee et après améliorée par des autres dans les années 70. On a trois étages pour résolution avec FDTD La région de solution est divise aux grilles de nœuds. Approximer les équations différentielles à l’équivalent de différence finie, qui relie les variable dépendante (temporelles ou spatiale) a un point dans la solution à un autre point au voisinage. Appliquer les conditions aux limites et les values initiales pour résoudre les différentes équations.

1.2 Les équations constitutives On considère un milieu linéaire, isotrope, et non-dispersif, les équations de maxwell peuvent s’écrites : 𝜌 ∇. 𝐸⃗ = ⁄𝜀

---------------(1)

⃗ =0 ∇. 𝐻

---------------(2) ⃗

𝜕𝐻 𝛻 ∧ 𝐸⃗ = −𝜇 𝜕𝑡

⃗ =𝜀 𝛻∧𝐻

⃗ 𝜕𝐷 𝜕𝑡

---------------(3) ---------------(4)

2.Methode de différance finie : Approximation à droite, à gauche et centrale Soit une fonction f(x), on peut trouver la dérivée de trois points sur l’arc sur ce droit ; l’un à droite, l’autres à gauche, et le troisième au centre :

𝑓 ′ (𝑥 ) ≈

𝑓(𝑥+𝛥𝑥)−𝑓(𝑥)

𝑓 ′ (𝑥 ) ≈

𝑓(𝑥)−𝑓(𝑥−∆𝑥)

𝑓 ′ (𝑥 ) ≈

𝑓(𝑥+𝛥𝑥)−𝑓(𝑥−∆𝑥)

𝛥𝑥

𝛥𝑥

2𝛥𝑥

=== > à droite === > à gauche === > au centre

3. Applications aux équations de maxwell : 3.1 Discrétisation temporelle : On va utiliser la même technique de différence finie central sur les équations de Maxwell. 𝛻 ∧ 𝐸(𝑡) = −𝜇

𝜕𝐻(𝑡) 𝜕𝑡

𝛻 ∧ 𝐻(𝑡) = 𝜀

𝐻(𝑡 + 𝛥𝑡) − 𝐻(𝑡) 𝛥𝑡

(3.1)

𝐸(𝑡 + 𝛥𝑡) − 𝐸(𝑡) 𝛥𝑡

(3.2)

𝛻 ∧ 𝐸(𝑡) = −𝜇

𝜕𝐸(𝑡)

𝛻 ∧ 𝐻(𝑡) = 𝜀

𝜕𝑡

Les parties droites des équations (4) existent au temps t et les parties gauches existent au temps 𝑡 + ∆𝑡. Cette formulation est instable car on n’a pas une simultanéité entre les deux champs E et H. C’est-a-dire pas de synchronisation. Pour rendre les équations stables, on doit ajuster le temps pour avoir les parties droites et les parties gauches à exister au même point de temps. C’est-à-dire, on va échantillonner le champ E et le champs H dans le temps. E existe aux temps entiers de (0, ∆𝑡, 2∆t, ….) et H existe à demi pas de temps (

∆𝑡 2

∆

∆

2

2

, 𝑡 + , 2𝑡 + , … ..)

Donc, on a les équations 3.1 et 3.2 modifiées comme suit : 𝛻 ∧ 𝐸|𝑡 = −𝜇

𝐻|𝑡+𝛥𝑡/2 − 𝐻|𝑡−∆𝑡/2 𝛥𝑡

𝛻 ∧ 𝐻|𝑡+𝛥𝑡∕2 = 𝜀

𝐸|𝑡+𝛥𝑡 − 𝐸|𝑡 𝛥𝑡

(4.1)

(4.2)

Maintenant on a les équations de mis à jour pour les champs E et H dans le temps futur :

𝐻 |𝑡+𝛥𝑡∕2 = 𝐻 |𝑡−𝛥𝑡∕2 − 𝐸 |𝑡+𝛥𝑡 = 𝐸 |𝑡 +

𝛥𝑡 (𝛻 ∧ 𝐸 | 𝑡 ) 𝜇

𝛥𝑡 (𝛻 ∧ 𝐻 |𝑡+∆𝑡/2 ) 𝜀

(5.1)

(5.2)

3.2 La grille de Yee pour FDTD La grille de Yee est un schéma explicite utilisant la méthode des différence finie . Les champs E et H sont continus et donc on ne peut pas stocker cette information dans la mémoire. Alors, proposé par Yee, la région de solution est divisée aux cellules discrètes appelé cellule de Yee.

Figure 2 : La grille ou maillage de Yee montrant une cellule unité de cotes ∆𝑥, ∆𝑦 𝑒𝑡 ∆𝑧.

On va échelonner la position de chaque composant des champs E et H dans une cellule de la grille comme indique sur la figure 4 ci-dessous.

Figure 3 : La grille de Yee 1D, 2D, 3D.

--- Approximations de la grille de Yee pour les équations Maxwell : -Les équations de divergences sont déjà satisfaites avec l’implémentations de la grille de Yee. 𝛻𝐵 = 0, 𝛻𝐸 = 0

Normalisation des champs On normalise le champs H pour avoir le même ordre magnitude pour H et E. ⃗

𝜕𝐻 𝛻 ∧ 𝐸⃗ = −[𝜇] 𝜕𝑡

⃗ = [𝜀] 𝛻∧𝐻

̃ = √𝜇0 𝐻 𝐻 𝜀0

𝛻∧𝐸 =−

̃ [𝜇𝑟 ]𝜕𝐻 𝑐0 𝜕𝑡

̃ = [𝜀𝑟 ]𝜕𝐸 𝛻∧𝐻

𝜕𝐸⃗ 𝜕𝑡

(6.1) (6.2)

(7.1) (7.2)

𝑐0 𝜕𝑡

Développement spatial des équations Donc, si on développe les équations (7.1) et (7.2) pour (x, y, z ; t) : (8.1)

𝜕𝐸𝑧 ⅆ𝑦

𝛻∧𝐸 =−

̃ [𝜇𝑟 ]𝜕𝐻 𝑐0 𝜕𝑡

𝜕𝐸𝑥 ⅆ𝑧 𝜕𝐸𝑦 ⅆ𝑥

𝜕𝐸𝑦

−

𝜕𝑧 𝜕𝐸𝑧

−

𝜕𝑥 𝜕𝐸𝑥

−

𝜕𝑦

̃𝑧 𝜕𝐻 ⅆ𝑦

̃= 𝛻∧𝐻

[𝜀𝑟 ]𝜕𝐸 𝑐0 𝜕𝑡

̃𝑥 𝜕𝐻 ⅆ𝑧

−

̃𝑦 𝜕𝐻 ⅆ𝑥

−

=−

1 𝐶0

̃𝑦 𝜕𝐻 𝜕𝑧

𝜕𝑥 ̃𝑥 𝜕𝐻 𝜕𝑦

= =

(𝜇𝑦𝑥

1 𝐶0

=

̃𝑥 𝜕𝐻

(𝜇𝑥𝑥

𝐶0

=−

̃𝑍 𝜕𝐻

−

1

=−

(𝜇𝑧𝑥

1 𝐶0 1 𝐶0

𝐶0

̃𝑥 𝜕𝐻 𝜕𝑡 ̃𝑥 𝜕𝐻

(𝜀𝑥𝑥

(𝜀𝑦𝑥

1

𝜕𝑡

(𝜀𝑧𝑥

𝜕𝑡

𝜕𝐸𝑥 𝜕𝑡 𝜕𝐸𝑥 𝜕𝑡 𝜕𝐸𝑥 𝜕𝑡

+ 𝜇𝑥𝑦 + 𝜇𝑦𝑦

+ 𝜀𝑧𝑦

̃𝑦 𝜕𝐻 𝜕𝑡

𝜕𝑡

+ 𝜀𝑥𝑦

+ 𝜇𝑥𝑧

𝜕𝑡

̃𝑦 𝜕𝐻

+ 𝜇𝑧𝑦

+ 𝜀𝑦𝑦

̃𝑦 𝜕𝐻

𝜕𝐸𝑦 𝜕𝑡 𝜕𝐸𝑦 𝜕𝑡 𝜕𝐸𝑦 𝜕𝑡

+ 𝜇𝑦𝑧 + 𝜇𝑧𝑧

+ 𝜀𝑥𝑧 + 𝜀𝑦𝑧 + 𝜀𝑧𝑧

Si on considère un tenseur diagonal, on a les matrices suivantes : 𝜀𝑥𝑥 [𝜀]=| 0 0

0 𝜀𝑦𝑦 0

0 0| 𝜀𝑧𝑧

𝜇𝑥𝑥 [𝜇]=| 0 0

0 𝜇𝑦𝑦 0

0 0 | 𝜇𝑧𝑧

Les 6 équations se réduisent à :

̃𝑥 𝜕𝐸𝑧 𝜕𝐸𝑦 1 𝜕𝐻 − = − 𝜇𝑥𝑥 𝜕𝑦 𝜕𝑧 𝐶0 𝜕𝑡

(9.1)

̃𝑦 𝜕𝐻 𝜕𝐸𝑥 𝜕𝐸𝑧 1 − = − 𝜇𝑦𝑦 𝜕𝑧 𝜕𝑥 𝐶0 𝜕𝑡

(9.2)

̃𝑧 𝜕𝐸𝑦 𝜕𝐸𝑥 1 𝜕𝐻 − = − 𝜇𝑧𝑧 𝜕𝑥 𝜕𝑦 𝐶0 𝜕𝑡

(9.3)

̃𝑦 ̃𝑧 𝜕𝐻 𝜕𝐻 1 𝜕𝐸𝑥 − = 𝜀𝑥𝑥 𝜕𝑦 𝜕𝑧 𝐶0 𝜕𝑡

(9.4)

̃𝑥 𝜕𝐻 ̃𝑍 𝜕𝐸𝑦 𝜕𝐻 1 − = 𝜀𝑦𝑦 𝜕𝑧 𝜕𝑥 𝐶0 𝜕𝑡

(9.4)

̃𝑦 𝜕𝐻 ̃𝑥 𝜕𝐻 1 𝜕𝐸𝑧 − = 𝜀𝑧𝑧 𝜕𝑥 𝜕𝑦 𝐶0 𝜕𝑡

(9.6)

̃𝑧 𝜕𝐻 𝜕𝑡 ̃𝑧 𝜕𝐻 𝜕𝑡 ̃𝑧 𝜕𝐻

(8.3)

)

(8.4)

)

(8.5)

)

(8.6)

𝜕𝐸𝑧

𝜕𝐸𝑧 𝜕𝑡 𝜕𝐸𝑧 𝜕𝑡

) (8.2)

)

𝜕𝑡

𝜕𝑡

)

L’implémentation de ces 6 équations dans la grille de Yee.  Pour 𝐻𝑥 :

𝑖,𝑗+1,𝑘

𝐸𝑧

𝑖,𝑗,𝑘

| −𝐸𝑧 𝑡

|

𝑡

𝛥𝑦

𝑖,𝑗,𝑘+1

−

𝐸𝑦

𝑖,𝑗,𝑘

| −𝐸𝑦 𝑡

|

𝑡

𝛥𝑧

=−

̃ 𝑖,𝑗,𝑘 𝑖,𝑗,𝑘 𝐻𝑥 |𝑡+𝛥

𝑢𝑥𝑥 𝑐0

𝑡/2

̃𝑥𝑖,𝑗.𝑘 | −𝐻

𝑡−𝛥𝑡/2

(10.1)

𝛥𝑡

Pour Hy :

𝑖,𝑗,𝑘+1

𝐸𝑥

𝑖,𝑗,𝑘

| −𝐸𝑥 𝑡

𝛥𝑧

|

𝑡

𝑖+1,𝑗,𝑘

−

𝐸𝑧

𝑖,𝑗,𝑘

| −𝐸𝑧 𝑡

𝛥𝑥

|

𝑡

=−

̃𝑦𝑖,𝑗,𝑘 | 𝑖,𝑗,𝑘 𝐻 𝑡+𝛥

𝑢𝑦𝑦 𝑐0

𝑡/2

̃𝑦𝑖,𝑗,𝑘 | −𝐻

𝛥𝑡

𝑡−𝛥𝑡/2

(10.2)

Pour Hz :

𝑖+1,𝑗,𝑘

𝐸𝑦

𝑖,𝑗,𝑘

| −𝐸𝑦 𝑡

𝛥𝑥

|

𝑡

𝑖,𝑗+1,𝑘

−

𝐸𝑥

𝑖,𝑗,𝑘

| −𝐸𝑥 𝑡

|

𝑡

=−

𝛥𝑦

̃𝑧𝑖,𝑗,𝑘 | 𝑖,𝑗,𝑘 𝐻 𝑡+𝛥

𝑢𝑦𝑦

𝑡/2

𝑐0

̃𝑧𝑖,𝑗.𝑘 | −𝐻

(10.3)

𝑡−𝛥𝑡/2

𝛥𝑡

Pour Ex :

̃𝑧𝑖,𝑗,𝑘 | 𝐻

𝑡+𝛥𝑡/2

̃𝑧𝑖,𝑗−1,𝑘 | −𝐻

𝑡+𝛥𝑡/2

𝛥𝑦

̃𝑦𝑖,𝑗,𝑘 | 𝐻

−

𝑡+𝛥𝑡/2

̃𝑦𝑖,𝑗,𝑘−1 | −𝐻

𝑡+𝛥𝑡/2

𝛥𝑧

𝑖,𝑗,𝑘

=

𝑖,𝑗,𝑘 𝜀𝑥𝑥 𝐸𝑥

𝑐0

|

𝑖,𝑗,𝑘

𝑡+∆𝑡

−𝐸𝑥

𝛥𝑡

|

𝑡

(10.4)

Pour Ey :

̃𝑥𝑖,𝑗,𝑘 | 𝐻

𝑡+𝛥𝑡/2

̃𝑥𝑖,𝑗,𝑘−1 | −𝐻

𝑡+𝛥𝑡/2

𝛥𝑧

̃𝑧𝑖,𝑗,𝑘 | 𝐻

−

𝑡+𝛥𝑡/2

̃𝑧𝑖−1,𝑗,𝑘 | −𝐻

𝑡+𝛥𝑡/2

𝛥𝑥

𝑖,𝑗,𝑘

=

𝑖,𝑗,𝑘 𝜀𝑦𝑦 𝐸𝑦

𝑖,𝑗,𝑘

|

𝑡+∆𝑡

𝑐0

−𝐸𝑦

|

𝑡

𝛥𝑡

(10.5)

Pour Ez :

̃𝑦𝑖,𝑗,𝑘 | 𝐻

𝑡+𝛥𝑡/2

̃𝑦𝑖−1,𝑗,𝑘 | −𝐻

𝑡+𝛥𝑡/2

𝛥𝑥

̃𝑥𝑖,𝑗,𝑘 | 𝐻

−

𝑡+𝛥𝑡/2

̃𝑥𝑖,𝑗−1,𝑘 | −𝐻

𝑡+𝛥𝑡/2

𝛥𝑦

𝑖,𝑗,𝑘

=

𝑖,𝑗,𝑘 𝜀𝑧𝑧 𝐸𝑧

𝑐0

|

𝑖,𝑗,𝑘

𝑡+∆𝑡

−𝐸𝑧

𝛥𝑡

|

𝑡

(10.6)

POUR 1D FDTD ∆x = ∆y = 0

Alors les équations 10.1-10.6 se réduisent à 6 équations suivantes : 𝑖,𝑗,𝑘+1

−

𝐸𝑦

𝑖,𝑗,𝑘

| −𝐸𝑦 𝑡

|

𝑡

𝛥𝑧 𝑖,𝑗,𝑘+1

𝐸𝑥

=−

𝑖,𝑗,𝑘

|𝑡 − 𝐸𝑥

|𝑡

𝛥𝑧

̃ 𝑖,𝑗,𝑘 ̃𝑥𝑖,𝑗.𝑘 | −𝐻 𝑖,𝑗,𝑘 𝐻𝑥 |𝑡+𝛥 𝑡−𝛥𝑡/2 𝑡/2

(11.1)

𝑢𝑥𝑥 𝑐0

𝛥𝑡 𝑖,𝑗,𝑘 ̃𝑖,𝑗,𝑘 𝑖,𝑗,𝑘 ̃ 𝑢𝑦𝑦 𝐻𝑦 |𝑡+𝛥𝑡/2 − 𝐻𝑦 |𝑡−𝛥𝑡/2

=−

𝑐0

𝛥𝑡

̃𝑧𝑖,𝑗,𝑘 = 0 𝐻 ̃𝑦𝑖,𝑗,𝑘 | 𝐻

𝑡+𝛥 𝑡

̃𝑦𝑖,𝑗,𝑘−1 | −𝐻

𝑡+

2

−

𝛥𝑡 2

𝛥𝑧

̃𝑥𝑖,𝑗,𝑘 | 𝐻 𝑡+𝛥

𝑖,𝑗,𝑘

=

𝛥𝑧 𝑖,𝑗,𝑘

𝐸𝑧

(11.3)

𝑖,𝑗,𝑘 𝜀𝑥𝑥 𝐸𝑥 = 𝑐0

̃𝑥𝑖,𝑗,𝑘−1 | −𝐻 𝑡+𝛥𝑡/2

𝑡/2

(11.2)

𝑖,𝑗,𝑘 𝜀𝑦𝑦

|𝑡+∆𝑡 − 𝐸𝑥𝑖,𝑗,𝑘 |𝑡 𝛥𝑡

𝑖,𝑗,𝑘 𝐸𝑦 | 𝑡+∆𝑡

𝑐0

(11.4)

−

𝑖,𝑗,𝑘 𝐸𝑦 | 𝑡

(11.5)

𝛥𝑡

=0

(11.6)

On remarque les équations de Maxwell sont couplés en deux équations Les équations (11.1) et (11.5) représente le mode Hx/Ey. Même les équations 11.2 et 11.4 représente le mode Ex/Hy. Les composants longitudinal Ez et Hz sont nulles. Puisqu’on est dans une seule dimension les indices i et j peuvent être éliminés à partir d’ici. On arrive alors, aux équations suivantes :

Mode Ey/Hx : −

𝐸𝑦𝑘+1 | −𝐸𝑦𝑘 | 𝑡

𝛥𝑧

̃𝑥𝑘 | 𝐻 𝑡+𝛥

𝑡/2

𝑡

=−

̃𝑥𝑘 | ̃𝑥𝑘 | −𝐻 𝑘 𝐻 𝑡+𝛥𝑡/2 𝑡−𝛥𝑡/2 𝑢𝑥𝑥 𝑐0

̃𝑥𝑘−1 | −𝐻 𝑡+𝛥𝑡/2 𝛥𝑧

(12.1)

𝛥𝑡 𝑘 𝑘 𝐸𝑘 | 𝜀𝑦𝑦 𝑦 𝑡+∆𝑡 − 𝐸𝑦 |𝑡 = 𝑐0 𝛥𝑡

(12.2)

Mode Ex/Hy : ̃𝑦𝑘 | 𝐻

𝑡+𝛥 𝑡 2

−

̃𝑦𝑘−1 | 𝛥𝑡 −𝐻 𝑡+ 2

𝛥𝑧

𝑘

=

(12.3)

𝑘

𝑘 𝐸𝑥 | −𝐸𝑥 |𝑡 𝜀𝑥𝑥 𝑡+∆𝑡

𝑐0

𝛥𝑡

̃𝑦𝑘 | ̃𝑦𝑘 | −𝐻 𝑘 𝐻 𝑢𝑦𝑦 𝐸𝑥𝑘+1 |𝑡 − 𝐸𝑥𝑘 |𝑡 𝑡+𝛥𝑡/2 𝑡−𝛥𝑡/2 =− 𝛥𝑧 𝑐0 𝛥𝑡

(12.4)

En résolvant les équations (12.1) , (12.2) , (12.3) et (12.4) pour trouver les nouvelles ̃𝑥𝑘 | ̃𝑦𝑘 | équations de mis à jour 𝐸𝑦𝑘 |𝑡+∆𝑡 et 𝐻 ; 𝐸𝑥𝑘 |𝑡+∆𝑡 et 𝐻 : 𝑡+𝛥 𝑡+𝛥 𝑡/2

𝑡/2

̃𝑥𝑘 | ̃𝑥𝑘−1 | 𝐻 −𝐻 𝑐 𝛥𝑡 𝑡+∆𝑡/2 𝑡+∆𝑡/2 0 = 𝐸𝑦𝑘 |𝑡 + ( 𝑘 ) ( ) 𝛥𝑧 𝜀𝑦𝑦

𝐸𝑦𝑘 |𝑡+𝛥

𝑡

̃𝑥𝑘 | 𝐻 𝑡+𝛥 /2 𝑡

𝐸𝑥𝑘 |𝑡+𝛥

𝑡

𝐸𝑦𝑘+1 |𝑡 − 𝐸𝑦𝑘 |𝑡 𝑐 𝛥𝑡 0 ̃𝑥𝑘 | =𝐻 + ( 𝑘 )( ) 𝑡−∆𝑡/2 𝛥𝑧 𝜇𝑥𝑥

(13.1)

(13.2)

̃𝑦𝑘 | ̃𝑦𝑘−1 | 𝐻 −𝐻 𝑐 𝛥𝑡 𝑡+∆𝑡/2 𝑡+∆𝑡/2 0 = 𝐸𝑥𝑘 |𝑡 − ( 𝑘 ) ( ) (13.3) 𝛥𝑧 𝜀𝑥𝑥

̃𝑦𝑘 | 𝐻 𝑡+𝛥𝑡 /2

=

̃𝑦𝑘 | 𝐻 𝑡−∆𝑡/2

𝑘+1 𝑘 𝑐0 𝛥𝑡 𝐸𝑥 |𝑡 − 𝐸𝑥 |𝑡 − ( 𝑘 )( ) 𝛥𝑧 𝜇𝑦𝑦

(13.4)

C0, ∆𝑡 𝑒𝑡 ∆𝑧 sont des constantes. Donc, on peut définir les deux coefficients suivants : 𝑚𝐸𝑘𝑦 =

𝑐0 ∆𝑡 𝑘 𝜀𝑦𝑦

𝑘 𝑘 𝑚𝐻 =𝑐0 ∆𝑡/𝜇𝑥𝑥 𝑥

𝑚𝐸𝑘𝑥 =

𝑐0 ∆𝑡 𝑘 𝜀𝑥𝑥

𝑘 𝑘 𝑚𝐻 =𝑐0 ∆𝑡/𝜇𝑦𝑦 𝑦

Condition aux limites parfaite (Perfect Boundary Conditions) : Soit le cas de 1D pour le mode Ey/Hx ci-dessus, et on suppose qu’aux extrémités du matériau, l’espace libre a des caractéristiques identiques à ce dernier. On prend l’extension de la grille de z=0 à z= 𝑁𝑧 ∆𝑧 . Donc, 𝑁𝑧 est le nombre d’itérations.

Si  les champs propagent de z=0 vers z= 𝑁𝑧 ∆𝑧, c-à-d vers l’extérieur  les matériaux aux deux extrémités doivent avoir les mêmes caractéristiques (linéaires, homogène, non-dispersive et même indice de réfraction 𝜂 ) 𝜂𝛥𝑧  2𝛥𝑡 = , l’onde propage 2∆𝑡 par cellule. 𝑐0

A la limite de z initiale : Seulement le champ E qui doit être modifier ici, c’est-à-dire équation (13.1). ℎ3 = ℎ2

̃𝑥1 ℎ1 = 𝐻

ℎ2 = ℎ1

1

̃ 𝑥| 𝐻 𝐸𝑦1 |𝑡+𝛥𝑡

=

𝐸𝑦1 |𝑡

+

𝑚𝐸𝑘𝑦

(

𝑡−∆𝑡/2

− ℎ3 )

𝛥𝑧

A la limite de z fin : Seulement le champ H qui doit être modifier ici, c’est-à-dire équation (13.2). 𝑒3 = 𝑒2

𝑒2 = 𝑒1 𝑁𝑧

𝑁

𝑒1 = 𝐸𝑦 𝑧 𝑁𝑧

𝑧 𝑒3 − 𝐸𝑁 𝑦 |

̃𝑥 | ̃𝑥 | 𝐻 =𝐻 + 𝑚𝑘𝐻𝑥 ( 𝑡+∆𝑡/2 𝑡−∆𝑡/2

𝛥𝑧

𝑡

)

2. Convergence Pour assurer la convergence du résultat

 On doit trouver la longueur d’onde la plus petite 𝜆𝑚𝑖𝑛 dans la grille. 𝜆𝑚𝑖𝑛 =

𝑐0 𝑓𝑚𝑎𝑥 𝑛𝑚𝑎𝑥

𝑛𝑚𝑎𝑥 = 𝑙′ 𝑖𝑛𝑑𝑖𝑐𝑒𝑑𝑒 𝑟𝑒𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛 𝑚𝑎𝑥𝑖𝑚𝑎𝑙 𝑡𝑟𝑜𝑢𝑣𝑒𝑟 𝑑𝑎𝑛𝑠 𝑙𝑎 𝑔𝑟𝑖𝑙𝑙𝑒 On va résoudre l’onde au moins 10 cellules 𝛥𝜆 ≤

𝜆𝑚𝑖𝑛⁄ 10

C’est-à-dire, on a 10 cellules par longueur d’onde.  La résolution doit être suffisante à résoudre la plus petite dimension sur la grille. Soit 𝑑𝑚𝑖𝑛 la plus petite structure dans une cellule dans la grille.

Donc, on peut

𝛥ⅆ ≈

ⅆ𝑚𝑖𝑛 𝑁𝑑

𝑜𝑛 𝑝𝑟𝑒𝑛𝑑 1 ≤ 𝑁ⅆ ≤ 4 𝑒𝑠𝑡 𝑙𝑎 𝑝𝑙𝑢𝑠 𝑒𝑓𝑓𝑖𝑐𝑎𝑐𝑒𝑠 𝑣𝑎𝑙𝑒𝑢𝑟𝑠

prendre : 𝛥𝑥 = 𝛥𝑦 = 𝛥𝑧 = min[𝛥𝜆 , 𝛥ⅆ ]

3. Stabilité Une OEM propageant dans l’espace libre ne doit pas dépasser la vitesse de la lumière 𝑐0 .Donc, le premier problème à régler pour éviter l’instabilité

numérique est le fait qu’en une itération temporelle , un point quelconque de l’onde ne doit pas propager dans plus une cellule . C’est-à-dire : Pour que l’onde propage dans la distance ∆𝑧 d’une cellule, on a besoin d’une durée de ∆𝑡 =

𝜂∆𝑧 𝑐0

.

ou’ 𝜂 =l’indice de réfraction la plus petite dans la grille. Dans l’espace libre on aura ∆𝑡 =

Pour 2D ∆𝑡 =

𝜂∆𝑧 √2𝑐0

𝜂∆𝑧 𝑐0

et pour 3D ∆𝑡 =

=

∆𝑧 𝑐0

𝜂∆𝑧 √3𝑐0

avec 𝜂 = 1

.

Pour assurer la stabilité de FDTD, l’incrémentation de temps ∆𝑡 doit satisfaire le critère de stabilité de Courant-Friedrichs-Lewy dans 3D exprimé par : 1

𝛥𝑡 ≤ 𝑣√

1 1 1 + 2+ 2 2 𝛥𝑥 𝛥𝑦 𝛥𝑧

On a 𝑣 =

𝑐0 𝜂

𝛥𝑡 ≤

𝜂 1 1 1 𝑐0 √ 2 + 2 + 2 𝛥𝑥 𝛥𝑦 𝛥𝑧

Si 𝛥𝑥 = 𝛥𝑦 = 𝛥𝑧=𝛿 : 𝛥𝑡 ≤

𝜂𝛿 𝑐0 √𝑛

𝑛 = 𝑙𝑎 𝑑𝑖𝑚𝑒𝑛𝑠𝑖𝑜𝑛 𝑑𝑒 𝑙𝑎 𝑔𝑟𝑖𝑙𝑙𝑒 ; 𝜂 = 𝐿′ 𝑖𝑛𝑑𝑖𝑐𝑒 𝑑𝑒 𝑟𝑒𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛 𝑑𝑒 𝑙𝑎 𝑔𝑟𝑖𝑙𝑙𝑒

Dans notre cas de 1D, la propagation se fait dans chaque 2∆𝑡 par cellule : 2𝛥𝑡 ≤

𝜂𝛥𝑧 𝑐0

Pendant une itération, un changement dans un champ électrique peut être remarqué seulement par un champ magnétique adjacent. On doit prendre 2∆𝑡 avant qu’on remarque le changement du champ électrique suivant, car les équations sont mises à jour l’une par l’autre.

Source d’excitation Des impulsions courtes sont utilisées dans l’excitation de FDTD. L’impulsion Gaussien est la plus utilisée dans l’excitation de FDTD. Elles nous permettent à simuler plusieurs fréquences en même temps. 𝑡 − 𝑡0 2 ) ) 𝑔(𝑡 ) = 𝑒𝑥𝑝 (( 𝜏 𝜏 𝑒𝑠𝑡 𝑙𝑎 𝑙𝑎𝑟𝑔𝑒𝑢𝑟 𝑑𝑒 𝑙′𝑖𝑚𝑝𝑢𝑙𝑠𝑖𝑜𝑛 Trouvant la transformer de Fourier de 𝑔(𝑡 ) : 𝑡 2 1 −𝑓 2 𝑔(𝑡 ) = 𝑒𝑥𝑝 (( ) ) ⇔ 𝐺 (𝑓) = exp ( 2 ) 𝜏 𝐵 √𝛱𝐵 𝐺 (𝑓) 𝑒𝑠𝑡 𝑎𝑢𝑠𝑠𝑖 𝑢𝑛𝑒 𝑎𝑢𝑠𝑠𝑖 𝑢𝑛𝑒 𝑓𝑜𝑛𝑐𝑡𝑖𝑜𝑛 𝐺𝑎𝑢𝑠𝑠𝑖𝑒𝑛

B est la fréquence maximale de 𝐺 (𝑓). Donc, 𝐵 =

1 𝜋𝜏

et on prend aussi 𝑡0 ≈ 6𝜏

 Déterminer 𝑓𝑚𝑎𝑥  Trouver la largeur de l’impulsion en utilisant le 𝑓𝑚𝑎𝑥 1 1 𝐵 = 𝑓𝑚𝑎𝑥 = → 𝜏 ≤ 𝜋𝜏

𝜋𝑓𝑚𝑎𝑥

Pour ajouter la source aux champs, tout simplement on ajoute la source au champ a un point donne sur la grille après la mise en jeux des champs. Cette méthode est appelée Simple Soft Source. Elle est plus efficace que la méthode Simple Hard Source ou les champs sont remplacé par la source a un point sur la grille. La SSS permet l’onde à se propager dans les deux directions.

̃𝑥𝑘 | ̃𝑥𝑘 | 𝐻 =𝐻 + 𝑔𝐻 | 𝑡+∆𝑡/2 𝑡+∆𝑡/2

Avec 𝑔𝐻 |𝑘

𝑘

𝐸𝑦𝑘 |𝑡+𝛥𝑡 = 𝐸𝑦𝑘 |𝑡+𝛥𝑡 + 𝑔𝐻 |𝑘

est la fonction Gaussienne.

Calcul du nombre de pas temporel : Le temps de propagation des champs dans la grille entière est : 𝑛𝑚𝑎𝑥 𝑁𝑧 𝛥𝑧 𝑡𝑝𝑟𝑜𝑝 = 𝑐0 T doit inclut le temps propagation complet dans l’impulsion : 𝑇 ≥ 12𝜏 On prend 5 temps de propagation d’une simulation sur la grille : 𝑇 ≥ 5𝑡𝑝𝑟𝑜𝑝 Donc le temps de propagation complet maintenant est : 𝑇 = 12𝜏 + 5𝑡𝑝𝑟𝑜𝑝 Alors, si ∆𝑡 est le pas de temps de propagation, on peut facilement trouver le nombre d’itérations : 𝑇 𝑁 = 𝑟𝑜𝑢𝑛𝑑 ↑ [ ] 𝛥𝑡

2.2 Algorithme pour FDTD :

ℎ3 = ℎ2 ̃𝑥1 ℎ3 = ℎ2 𝐻 ̃𝑥1 ℎ1 = 𝐻

Computer la résolution de la grille : 𝑧 ′ = 𝑚𝑖𝑛 [ 𝛥𝑧 =

𝜆𝑚𝑖𝑛 ⅆ𝑚𝑖𝑛 , ], 𝑁𝜆 𝑁𝑑

𝑁 = 𝑟𝑜𝑢𝑛𝑑 ↑ [

FIN?

ⅆ𝑧 ] 𝛥𝑧 ′

ℎ2 = ℎ1

Terminé.

ℎ2 = ℎ1 ℎ3 = ℎ2 ̃𝑥1 ℎ1 = 𝐻

Mise en jeux H à partir de E :

ⅆ𝑧 𝑁

ℎ3 = ℎ2 Le pas de temps :

ℎ2 = ℎ1

Initialise les H et E a zero :

𝑛𝑚𝑎𝑛 ∆𝑧 𝛥𝑡 = 2𝑐0

E=H=0

Initialise les H et E a zero : Computer la source : 𝜏=

ℎ1 = ℎ2 = ℎ1

0⋅5 ,𝑡 𝑓𝑚𝑎𝑥 0

E=H=0

= 6𝜏, 𝑔(𝑡) =

𝑡−𝑡 𝑒𝑥𝑝 [− ( 𝜏 0 )]

Initialise les H et E a zero : E=H=0 Calcul des coefficients : 𝑚𝐸𝑘𝑦 =

𝑐0 𝛥𝑡 𝑘 ∆𝑧 𝜀𝑦𝑦

𝑐 𝛥𝑡

𝑘 𝑚𝐻 = 𝜇 0 ∆𝑧 𝑥 𝑥𝑥

Initialise les H et E a zero : E=H=0

Initialise les H et E a zéro : 𝐸𝑦𝑘 =𝐻𝑥𝑘 =0

E=H=0

̃𝑥1 ℎ1 = 𝐻

Figure 1 : Algorithme pour FDTD. 𝐻|𝑡+𝛥𝑡/2 𝐸|𝑡 𝐸|𝑡+𝛥𝑡 𝐻|𝑡+∆𝑡/2 𝐻|𝑡+3𝛥𝑡/2 𝐸|𝑡+∆𝑡 𝐸|𝑡+2𝛥𝑡 𝐻|𝑡+3∆𝑡/2 𝐻|𝑡+5𝛥𝑡/2 𝐸|𝑡+2∆𝑡 …...

Implémentation de 1D FDTD Matlab On peut modifie notre algorithme dans la figure 1 pour 1D.

Computer les coefficients de mise a jeux :

mEy=

𝑐0 ∆𝑡 𝑘 𝜀𝑦𝑦

𝑘 mHy=𝑐0 ∆𝑡/𝜇𝑥𝑥

Initialise les H et E a zéro : E=H=0

𝛻

FIN?

Terminé.

Calcule H à partir de E en utilisant l’équation (6.1)

𝐻 |𝑡+𝛥𝑡∕2 = 𝐻 |𝑡−𝛥𝑡∕2 −

𝛥𝑡 𝜇

(𝛻 ∧ 𝐸 | 𝑡 )

Calcule E à partir de H en utilisant l’équation (6.2)

𝐸 |𝑡+𝛥𝑡 = 𝐸 |𝑡 +

𝛥𝑡 𝜀

(𝛻 ∧ 𝐻 |𝑡+∆𝑡/2 )

Figure 1 : Algorithme pour FDTD.

Les etapes de codages : 1. 2. 3. 4. 5. 6.

Implémenter l’algorithme FDTD Ajouter la source Ajouter la condition aux limites Modifier la source pour être dans une seule direction Calculer les coefficients de transmission et de réflexion. Ajouter un composant

Conditions aux limites numérique (Conditions aux limites de Dirichlet): Si Nz est le nombre d’itération

On prendre et 𝑐0 𝛥𝑡 ̃𝑥𝑘 | ̃𝑥𝑘 | 𝐻 = 𝐻 + ( 𝑘 )( 𝑡+𝛥 /2 𝑡−∆𝑡/2

̃𝑥𝑘 | ̃𝑥𝑘 | 𝐻 =𝐻 + mHy ( 𝑡+𝛥 /2 𝑡−∆𝑡/2 𝑡

𝐸𝑦𝑘 |𝑡+𝛥 = 𝐸𝑦𝑘 |𝑡 + ( 𝑡

𝑡

𝑡

𝜇𝑥𝑥

𝑡

𝐸𝑦𝑘 |𝑡+𝛥

𝐸𝑦𝑘+1 | −𝐸𝑦𝑘 |

=

𝐸𝑦𝑘 |𝑡

𝑐0 𝛥𝑡 𝑘 𝜀𝑦𝑦

𝛥𝑧

0−𝐸𝑦𝑘 |

𝑡

𝛥𝑧

𝑡

) pour k
) pour k=Nz

̃𝑥𝑘 | ̃𝑥𝑘−1 | 𝐻 −𝐻 𝑡+∆𝑡/2 𝑡+∆𝑡/2

)(

𝛥𝑧

̃𝑥𝑘 | 𝐻 −0 𝑡+∆𝑡/2

+ mEy (

𝛥𝑧

) Pour k>1

) Pour k=1