Qu'est-ce que la méthode FDTD ?

1. Introduction

La méthode des différences finies dans le domaine temporel (FDTD) résout directement les équations de Maxwell – qui décrivent la lumière – dans le domaine temporel, sur une grille spatiale discrète. Contrairement au lancer de rayons, aux solveurs de modes ou aux approches par matrice de transfert, la FDTD ne fait aucune hypothèse simplificatrice concernant la propagation ou la structure modale. Au lieu de cela, les champs électriques et magnétiques sont calculés explicitement dans le temps en discrétisant les équations de rotation – la loi d’Ampère et la loi de Faraday – de sorte que les ondes électromagnétiques se propagent, se diffusent, interfèrent et se diffractent naturellement dans tout le domaine de simulation. Cette généralité rend la FDTD très précise et polyvalente, mais elle est également exigeante en termes de ressources de calcul, notamment en mémoire et en temps de traitement.

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

Étant donné que la méthode FDTD ne repose sur aucune hypothèse concernant la direction de propagation, la structure du dispositif ou le comportement du champ, une seule simulation peut prendre en compte la réponse sur une large bande de fréquences, la diffraction et l’interférence, les effets de champ proche, la géométrie à l’échelle sub-longueur d’onde et les phénomènes électromagnétiques transitoires, le tout dans un cadre unique. Cette généralité est ce qui rend la méthode FDTD si puissante, mais elle a aussi un inconvénient : la grille spatiale doit être suffisamment fine pour résoudre la plus petite échelle de longueur physique, le pas de temps est fortement limité par la stabilité numérique, et les simulations tridimensionnelles deviennent rapidement exigeantes en termes de mémoire et de temps de calcul, même pour les géométries qui semblent modestes en taille.

2. Pourquoi il ne faut pas utiliser la méthode FDTD...

Dans la plupart des problèmes d’optique pratiques, la méthode FDTD ne devrait pas être la première option envisagée. Elle est gourmande en ressources informatiques, et pour de nombreux systèmes, il existe des approches beaucoup plus efficaces. La conception de systèmes optiques tels que des appareils photo ou des télescopes, la modélisation de systèmes d’imagerie à grande échelle ou le suivi du trajet de la lumière à travers des assemblages optiques macroscopiques sont beaucoup mieux réalisés à l’aide de la méthode du tracé de rayons. De même, les structures planaires multicouches, telles que les cellules solaires à couches minces, les piles OLED et les revêtements optiques, sont traitées de manière beaucoup plus efficace à l’aide de méthodes de matrices de transfert. Dans ces cas, la méthode FDTD entraîne une augmentation significative des coûts de calcul sans apporter d’informations physiques supplémentaires.

La plupart des méthodes de simulation électromagnétique sont ultimement basées sur les équations de Maxwell. Des techniques telles que les méthodes de résolution des modes et les méthodes de matrice de transfert sont obtenues en introduisant des hypothèses simplificatrices concernant la géométrie, la dimensionnalité ou la structure du champ, afin d’accélérer les calculs. Ces approximations sont extrêmement efficaces, mais uniquement lorsque les hypothèses sous-jacentes sont valides.

La méthode FDTD se situe à l’extrémité opposée de cet éventail. Elle est extrêmement performante, en particulier pour les systèmes tridimensionnels complexes, les géométries irrégulières et les problèmes dans lesquels le champ ne peut pas être décomposé en modes ou en couches simples, comme les cristaux photoniques, les structures de taille inférieure à la longueur d’onde, les milieux fortement diffusants et le couplage en champ proche. Le prix de cette flexibilité est un coût de calcul élevé, tant en termes de mémoire que de temps de traitement.

L’élément essentiel est de choisir l’outil approprié. Utilisez des modèles simplifiés lorsque les lois de la physique le permettent, et utilisez la méthode FDTD lorsque ces hypothèses ne sont plus valables. Ses principaux avantages sont sa flexibilité et sa capacité à modéliser les phénomènes physiques de manière complète ; en contrepartie, il faut tenir compte du coût de calcul.

3. Les équations de Maxwell à la base de la méthode FDTD.

Au fond, la méthode FDTD est très simple. On prend les équations de Maxwell, on les décompose en leurs composantes pour chaque direction spatiale, puis on résout ces équations directement sur une grille. Il n’y a pas de magie : la méthode est conceptuellement simple et, à bien des égards, plus facile à mettre en œuvre que les approches numériques plus complexes, comme les solveurs de diffusion-dérive. Pour commencer, nous allons développer la loi d’Ampère. Dans un milieu linéaire doté d’une conductivité \(\sigma\), d’une permittivité \(\epsilon\) et d’une perméabilité \(\mu\), on peut l’écrire ainsi :

\[ \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t} = \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix} \]

ce qui peut être développé comme suit :

La loi de Faraday s’énonce comme suit :

\[ -\sigma_m \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t} = \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix} \]

ce qui peut être développé pour obtenir :

Les équations mises en évidence en bleu et en rouge représentent la définition de la méthode FDTD et constituent le cœur de cette méthode. Rien de plus complexe. Tout le reste consiste en des étapes de mise en œuvre : remplacer les dérivées par des différences finies, choisir… conditions aux limites, définissant sources de lumière (excitations), en extrayant les résultats à l’aide de détecteurset en rédigeant un code efficace (souvent parallèle) pour faire évoluer les champs dans le temps. Ces éléments supplémentaires rendent un solveur FDTD pratique plus complexe, mais ils ne modifient pas la simplicité fondamentale de la méthode de base.

4. Un exemple de simulation 1D FDTD à des fins ludiques.

Pour montrer à quel point la méthode FDTD est simple, tant sur le plan conceptuel que sur le plan de sa mise en œuvre, nous allons maintenant réduire les équations complètes à un modèle simplifié unidimensionnel qui peut être implémenté en quelques dizaines de lignes de code MATLAB ou Python. L’objectif ici n’est pas de modifier les principes physiques, mais de supprimer la complexité inutile afin que la structure de la méthode soit parfaitement claire.

Pour ce faire, nous ne prenons en compte que la variation dans la direction \(x\) et nous fixons explicitement les dérivées spatiales restantes à zéro : \[ \frac{\partial}{\partial y} = 0, \qquad \frac{\partial}{\partial z} = 0. \]. Il en résulte un système électromagnétique transversal doté d’un champ électrique \(E_z(x,t)\) et d’un champ magnétique \(H_y(x,t)\).

En partant de ces hypothèses, les équations de Maxwell se simplifient et se réduisent à :

\[ \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} \]

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x}. \]

Voici les équations continues fondamentales pour le cas unidimensionnel. Cependant, sous cette forme, elles ne peuvent pas être résolues directement à l’aide d’un ordinateur. Pour ce faire, nous remplaçons les dérivées par des différences finies sur une grille discrète. C’est là l’essence de la méthode : voir différences finies pour une vue d’ensemble.

En introduisant une grille spatiale \(x_i = i\Delta x\) et des pas de temps \(t^n = n\Delta t\), et en utilisant la disposition standard décalée (méthode de Verlet), les équations de mise à jour deviennent :

\[ E_z^{n+1}(i) = E_z^n(i) + \frac{\Delta t}{\epsilon(i)\,\Delta x} \left[ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2}\right) \right] \]

\[ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2}\right) = H_y^{n-\frac{1}{2}}\!\left(i+\frac{1}{2}\right) + \frac{\Delta t}{\mu\,\Delta x} \left[ E_z^n(i+1) - E_z^n(i) \right]. \]

Il s’agit déjà de la méthode FDTD dans sa forme la plus simple. Le champ électrique est mis à jour en fonction des différences spatiales du champ magnétique, et le champ magnétique est mis à jour en fonction des différences spatiales du champ électrique, les deux champs se succédant dans le temps.

Pour un milieu non magnétique, \(\mu = \mu_0\), la permittivité est liée à l’indice de réfraction par la relation suivante :

\[ \epsilon(x) = \epsilon_0 n(x)^2. \]

Cela signifie qu’une variation de l’indice de réfraction est introduite en modifiant simplement la valeur de \(\epsilon\) sur la grille. Dans l’exemple ci-dessous, nous considérons un milieu qui passe de \(n=1\) à \(n=3\), de sorte qu’une impulsion lancée depuis la gauche est partiellement réfléchie et partiellement transmise à l’interface.

Dans la pratique, les équations ci-dessus sont simplement intégrées dans une boucle temporelle : le champ magnétique est mis à jour à partir du champ électrique, les conditions aux limites sont appliquées, puis le champ électrique est mis à jour à partir du champ magnétique, et le processus se répète. Les champs sont décalés dans le temps de manière à ce qu’ils évoluent de manière simple et explicite, chaque mise à jour ne dépendant que des valeurs locales des points voisins. Cela est illustré schématiquement dans ??, et c’est essentiellement tout ce que l’ordinateur effectue à chaque pas de temps.

5. Un exemple de simulation 1D FDTD pour jouet, réalisé en Matlab/Python.

Le même schéma FDTD unidimensionnel peut être mis en œuvre en quelques dizaines de lignes de code seulement. Ci-dessous, vous trouverez des implémentations équivalentes en MATLAB et en Python. Les deux programmes génèrent une impulsion gaussienne et montrent la réflexion à l’endroit d’un saut d’indice de réfraction.

clear;
clc;

% Physical constants
c0   = 299792458;
mu0  = 4*pi*1e-7;
eps0 = 1/(mu0*c0^2);

% Grid
Nx = 400;
dx = 20e-9;
dt = dx/(2*c0);

% Material profile
n = ones(1, Nx);
n(220:end) = 3.0;
eps_r = n.^2;

% Fields
Ez = zeros(1, Nx);
Hy = zeros(1, Nx-1);

% Source
source_pos = 80;
t0 = 40;
spread = 12;

Nt = 600;

figure;

for t = 1:Nt

    % Update H
    for i = 1:Nx-1
        Hy(i) = Hy(i) + (dt/(mu0*dx))*(Ez(i+1) - Ez(i));
    end

    % Update E
    for i = 2:Nx-1
        Ez(i) = Ez(i) + (dt/(eps0*eps_r(i)*dx))*(Hy(i) - Hy(i-1));
    end

    % Source
    Ez(source_pos) = Ez(source_pos) + exp(-((t - t0)/spread)^2);

    % Simple absorbing boundaries
    Ez(1)   = Ez(2);
    Ez(end) = Ez(end-1);

    % Plot every few steps
    if mod(t,5) == 0
        plot(1:Nx, Ez, 'b', 'LineWidth', 1.5);
        hold on;
        plot([220 220], [-1.5 1.5], '--r', 'LineWidth', 1.0);
        hold off;
        xlabel('Grid index');
        ylabel('E_z');
        title(['1D FDTD pulse reflecting from an index step, t = ', num2str(t)]);
        axis([1 Nx -1.5 1.5]);
        grid on;
        drawnow;
    end
end

import numpy as np
import matplotlib.pyplot as plt

c0 = 299792458.0
mu0 = 4.0e-7 * np.pi
eps0 = 1.0 / (mu0 * c0**2)

Nx = 400
dx = 20e-9
dt = dx / (2 * c0)

n = np.ones(Nx)
n[220:] = 3.0
eps_r = n**2

Ez = np.zeros(Nx)
Hy = np.zeros(Nx - 1)

source_pos = 80
t0 = 40
spread = 12

Nt = 600

plt.ion()
fig, ax = plt.subplots()

for t in range(Nt):

    Hy += (dt / (mu0 * dx)) * (Ez[1:] - Ez[:-1])
    Ez[1:-1] += (dt / (eps0 * eps_r[1:-1] * dx)) * (Hy[1:] - Hy[:-1])

    Ez[source_pos] += np.exp(-((t - t0) / spread) ** 2)

    Ez[0] = Ez[1]
    Ez[-1] = Ez[-2]

    if t % 5 == 0:
        ax.clear()
        ax.plot(Ez)
        ax.axvline(220, linestyle='--')
        ax.set_xlim(0, Nx)
        ax.set_ylim(-1.5, 1.5)
        plt.pause(0.01)

plt.ioff()
plt.show()

6. Lancer la simulation.

Le comportement de la simulation est illustré dans les figures ?? et ??. Au départ, l’impulsion gaussienne se propage sans problème à travers la région \(n=1\). Lorsqu’elle atteint la zone de changement d’indice, le champ se divise en deux composantes : une impulsion réfléchie qui se propage en sens inverse dans la région à faible indice, et une impulsion transmise qui pénètre dans la région à indice plus élevé. Comme l’indice de réfraction augmente jusqu’à atteindre la valeur \(n=3\), l’onde transmise ralentit et sa longueur d’onde diminue, tandis que la composante réfléchie conserve les caractéristiques de propagation d’origine. Il s’agit précisément du comportement physique attendu, ce qui démontre que même cette implémentation minimale permet de reproduire les principes fondamentaux de la propagation, de la réflexion et de la transmission des ondes.

11. Quelle sera la prochaine étape ?

Les exemples MATLAB et Python ci-dessus illustrent la simplicité du principe de la méthode FDTD, mais il ne s’agit que d’implémentations simplifiées. En pratique, les simulations optiques réalistes nécessitent des grilles de grande taille, de longs temps de calcul et une gestion rigoureuse des limites et des sources. Pour appliquer la méthode FDTD à des problèmes concrets et utiles, il faut un solveur hautement optimisé, capable de fonctionner en parallèle et sur des processeurs (CPU) et des cartes graphiques (GPU).

C’est précisément ce qu’offre OghmaNano : un solveur 3D FDTD complet, écrit en C et doté d’une accélération OpenCL, qui permet d’exécuter le même modèle de manière fluide sur des processeurs multicœurs ou des cartes graphiques. Au lieu de créer votre propre solveur à partir de zéro, vous pouvez vous inspirer des idées présentées ici et les appliquer immédiatement à des structures photoniques, des dispositifs et des problèmes d’optique ondulatoire réalistes.

Pour commencer, vous pouvez explorer les Méthode FDTD dans OghmaNano. page présentant un aperçu des fonctionnalités, à consulter de nouveau qu’est-ce que la méthode FDTD ou le dérivation mathématique complèteou passez directement à la configuration pratique en utilisant sources de lumière, conditions aux limiteset détecteurs.

Une vaste gamme de exemples prêts à l’emploi sont également disponibles, notamment : propagation dans l’espace libre, interférence à double fente, guides d’ondes à cristaux photoniques, réseaux de Bragg, Cavités de Fabry-Pérot, facettes d’un guide d’ondes inclinées, structures fuselées, Dispositifs de Mach-Zehnder, résonateurs annulaireset revêtements antireflets.

En résumé : la méthode de base est simple, comme illustré ci-dessus, mais OghmaNano est conçu pour exploiter cette simplicité et la développer en un outil pratique et performant pour la simulation optique réelle.