Dans la page consacrée aux équations de Maxwell, nous avons énoncé, parmi ces quatre équations, les deux équations qui définissent l'évolution du champ électrique et du champ magnétique. Dans le vide et en l'absence de charges et de courants électriques, elles sont :
\( \vec{rot} \: \vec{E} = -\frac{\partial \vec{B}}{\partial t} \)
\( \vec{rot} \: \vec{B} = \mu_0 \epsilon_0 \frac{\partial \vec{E}}{\partial t} \)
Réarrangeons-les pour faire apparaître les variations temporelles de chaque champ :
\( \frac{\partial \vec{E}}{\partial t} = \frac{1}{\mu_0 \epsilon_0} \vec{rot} \: \vec{B} \)
\( \frac{\partial \vec{B}}{\partial t} = - \vec{rot} \: \vec{E} \)
Pour facilier l'exercice, nous allons choisir de travailler en dimension 1, c'est à dire avec une propagation des champs sur un seul axe, que par tradition je vais nommer l'axe Oz. Je vais également polariser mon champ électrique sur l'axe Ox et mon champ magnétique sur l'axe Oy.
Sous ces hypothèses, reprenons l'expression de la variation temporelle du champ électrique, en utilisant les coordonnées cartésiennes du rotationel (voir si besoin la page sur les opérateurs différentiels) :
\( \frac{\partial \vec{E}}{\partial t} = \frac{1}{\mu_0 \epsilon_0} \left| \begin{array}{cc} \vec{e_x}&\vec{e_y}&\vec{e_z}\\ \frac{\partial}{\partial x}&\frac{\partial}{\partial y}&\frac{\partial}{\partial z}\\ B_x&B_y&B_z \end{array} \right| \)
J'ai supposé que le champ \( \vec{E} \) n'avait qu'une seule composante sur Ox, ce qui me donne donc :
\( \frac{\partial E_x}{\partial t} \vec{e_x} = \frac{1}{\mu_0 \epsilon_0} (\frac{\partial B_z}{\partial y} - \frac{\partial B_y}{\partial z}) \vec{e_x} \)
soit en projetant sur Ox, et en remarquant que \( \frac{\partial B_z}{\partial y} = 0 \) :
\( \frac{\partial E_x}{\partial t} = -\frac{1}{\mu_0 \epsilon_0} \frac{\partial B_y}{\partial z} \)
Avec les mêmes calculs et en se rappelant que le champ \( \vec{B} \) n'a qu'une seule composante sur Oy, j'obtiens l'autre équation :
\( \frac{\partial B_y}{\partial t} = -\frac{\partial E_x}{\partial z} \)
J'obtiens ainsi deux équations qui décrivent l'évolution des champs électrique et magnétique.
Nous allons utiliser la méthode FDTD, pour Finite Difference Time Domain, et un algorithme particulier dû à Kane Yee en 1966. La discrétisation des équations se fait de manière habituelle en utilisant l'approximation des différences centrées de second ordre. La particularité de l'algorithme de Yee est de permettre de résoudre l'évolution des deux champs sur la même grille et avec le même schéma. La subtilité consiste à utiliser des demi-pas pour discrétiser le temps et l'espace, et d'imbriquer les calculs sur le domaine spatial et temporel. Vous aurez d'ailleurs noté que dans notre système d'équations, le temps, l'espace et les deux champs étaient étroitement imbriqués ...
Plus précisément, le champ électrique est calculé pour des demi-pas de temps et des pas entiers d'espace, alors que le champ magnétique est calculé pour des pas entiers de temps et des demi-pas d'espace.
Avec j l'indice spatial et n l'indice temporel, avec \( \Delta z \) et \( \Delta t \), respectivement le pas spatial et le pas temporel, je noterai chaque pas de champ électrique \( E_x(z,t) = E_x(j\Delta z,n\Delta t) = E_x^{j,n}\), avec bien sur la même notation pour le champ magnétique.
En appliquant l'algorithme de Yee et en utilisant les notations définies ci-dessus, j'obtiens pour le champ électrique l'équation discrète :
\( \frac{E_x^{j,n+\frac{1}{2}} - E_x^{j,n-\frac{1}{2}} }{\Delta t} = -\frac{1}{\mu_0 \epsilon_0} \frac{B_y^{j+\frac{1}{2},n} - B_y^{j-\frac{1}{2},n}}{\Delta z} \)
et pour le champ magnétique :
\( \frac{B_y^{j+\frac{1}{2},n+1} - B_y^{j+\frac{1}{2},n} }{\Delta t} = -\frac{E_x^{j+1,n+\frac{1}{2}} - E_x^{j,n+\frac{1}{2}}}{\Delta z} \)
En arrangeant ces deux équations pour obtenir un schéma explicite et en posant \( c^2 = \frac{1}{\mu_0 \epsilon_0}\) j'obtiens :
\( E_x^{j,n+\frac{1}{2}} = E_x^{j,n-\frac{1}{2}} - \frac{c^2 \Delta t}{\Delta z }(B_y^{j+\frac{1}{2},n} - B_y^{j-\frac{1}{2},n}) \)
\( B_y^{j+\frac{1}{2},n+1} = B_y^{j+\frac{1}{2},n} - \frac{\Delta t}{\Delta z } (E_x^{j+1,n+\frac{1}{2}} - E_x^{j,n+\frac{1}{2}}) \)
Vous le savez maintenant, les schémas explicites peuvent être instables numériquement. Leur stabilité dépend du choix des pas spatiaux et temporels, qui doivent correspondre à une réalité physique. Considérons le cas du pas spatial : pour que notre calcul est un sens, il faut que ce pas soit petit devant la longueur d'onde de l'ondulation du champ électromagnétique.
Nous avons déjà vu que la stabilité du schéma était conditionnée au respect du critère CFL, qui dans notre cas veut que \( \frac{c^2 \Delta t}{\Delta z } \le \frac{1}{2} \). Il faudra en tenir compte lors de l'implémentation du modèle. Dans le script, je poserai le critère CFL \( \alpha = \frac{c^2 \Delta t}{\Delta z } \)
L'algorithme de Yee utilise des indices fractionnaires. Fort bien, mais comment implémenter ces indices dans un programme informatique ? A priori, je ne connais pas de langage qui admette des indices de tableau fractionnaires ... En fait, il faut se souvenir que ces indices fractionnaires sont un artifice de calcul pour symboliser l'imbrication des champs électrique et magnétique. L'évolution des deux champs doit être calculée simultanément, c'est à dire informatiquement dans la même itération de temps. Mais évidemment, dans nos programmes, cette itération sera entière.
Les conditions initiales de notre système d'équations sont fixées par la nature de la source des champs électrique et magnétique. J'ai choisi une source sinusoïdale qui sera définie dans la fonction Source(). Vous pouvez modifier les paramètres de cette source sinusoïdale et même définir un autre type de source, par exemple une gaussienne.
La propagation d'un champ électromagnétique peut s'effectuer dans un espace quasi-infini ou bien au contraire dans une cavité fermée. Ce sont les conditions aux limites qui déterminent le comportement du champ EM dans ce cas. Mais même dans le cas d'une propagation dans un espace infini, il nous faudra bien borner notre domaine, car la mémoire de nos ordinateurs n'a rien d'infini, malheureusement ! Il va donc falloir simuler un espace (et donc une grille de calcul) infini.
La solution la plus simple, c'est de poser \(\frac{c^2 \Delta t}{\Delta z } = \frac{1}{2} \), soit \( \Delta t = \frac{\Delta z }{2c^2} \). Dans ce cas, le champ est complétement absorbé aux limites. Il existe une manière plus élégante de définir les conditions aux limites, qui a été proposée par Gerrit Mur et que je n'aborderai pas ici (voir par exemple ici).
Visualiser la propagation du champ électromagnétique nécessite de disposer de moyens d'animation. MatPlotLib fournit un module de ce genre (animate), mais il ne me satisfait pas pour l'application que je compte faire ici. Nous allons donc utiliser le module visual de Python, qui est bien sous tous rapports mais qui nécessite de disposer de l'IDE VPython.
Je vous encourage donc à charger cet IDE depuis SourceForge pour votre système d'exploitation (Linux, MacOS ou windows) et votre version de Python (la 2.7.9 pour moi). Son utilisation est en tout point semblable à celle de l'IDE Python "normal".
Il est semble-t-il possible d'appeler le module visual depuis anaconda, mais je n'y suis pas arrivé. Je n'ai pas non plus passé beaucoup de temps à comprendre et résoudre le problème ...
Je ne vais pas décrire tout le code, que vous trouverez dans le script OEMPropagation_V.py, mais seulement les points marquants. Presque la moitié du code consiste d'ailleurs à mettre en place les définitions de l'animation et l'affichage des deux courbes de champ.
Arrêtez-nous sur la fonction qui définit la source des champs
def Source(x):
A0 = 0.1
omega = 2*pi/100.0
s = A0*sin(omega*x)
return s
C'est une sinusoïde, dont vous pourrez modifier l'amplitude et la pulsation. Attention toutefois aux effets de bord, que vous pouvez deviner en lisant attentivement les paragraphes précédants...
La définition des conditions initiales se fait par :
X = arange(Xmax)
Ex[:Xmax,0] = Source(X)
By[:Xmax,0] = Source(X)
Après avoir tracé le champ EM initial par :
PlotChamp(0)
la boucle de calcul et d'affichage du champ est lancée. On commence par appliquer le schéma défini ci-dessus :
Ex[1:Xmax-1,1] = Ex[1:Xmax-1,0] + Alpha*(By[0:Xmax-2,0] - By[2:Xmax,0])
By[1:Xmax-1,1] = By[1:Xmax-1,0] + Alpha*(Ex[0:Xmax-2,0] - Ex[2:Xmax,0])
puis on applique les conditions aux limites (les conditions de bord) sur les deux champs E et B :
Ex[0,1] = Ex[0,0] + Alpha*(By[Xmax-2,0] - By[1,0])
Ex[Xmax-1,1] = Ex[Xmax-1,0] + Alpha*(By[Xmax-2,0] - By[1,0])
By[0,1] = By[0,0] + Alpha*(Ex[Xmax-2,0] - Ex[1,0])
By[Xmax-1,1] = By[Xmax-1,0] + Alpha*(Ex[Xmax-2,0] - Ex[1,0])
la valeur calculée des deux champs est affichée :
PlotChamp(1)
puis on swappe les valeurs de champ qui viennent d'être calculées pour qu'elles deviennent le nouveau support de calcul :
Ex[:Xmax,0] = Ex[:Xmax,1]
By[:Xmax,0] = Ex[:Xmax,1]
La boucle tourne indéfiniment, jusqu'à ce que vous arrêtiez le script en fermant la fenêtre VPython. Vous pouvez régler la vitesse de défilement en modifiant le paramètre FreqAff.
Après lancement de votre script sous VPython, vous obtiendrez l'animation suivante, qui montre la propagation liée du champ électrique et magnétique :
Le script Python OEMPropagation_V.py est disponible dans le package OEMPropagation.zip.
Contenu et design par Dominique Lefebvre - tangenteX.com décembre 2016 -- Pour me joindre par ou PhysiqueX ou 
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.