Co to jest metoda FDTD?

1. Wprowadzenie

Metoda różnic skończonych w dziedzinie czasu (FDTD) rozwiązuje równania Maxwella – opisujące światło – bezpośrednio w dziedzinie czasu, wykorzystując dyskretną siatkę przestrzenną. W przeciwieństwie do metod śledzenia promieni, rozwiązywania równań modalnych lub metod macierzowych, FDTD nie wprowadza żadnych upraszczających założeń dotyczących propagacji lub struktury modalnej. Zamiast tego, pola elektryczne i magnetyczne są obliczane w sposób bezpośredni w czasie poprzez dyskretyzację równań rotacji – prawa Ampère’a i prawa Faradaya – co pozwala na naturalną propagację, rozpraszanie, interferencję i dyfrakcję fal elektromagnetycznych w całej dziedzinie symulacji. Ta ogólność sprawia, że FDTD jest bardzo dokładną i wszechstronną metodą, ale jednocześnie wymaga znacznych zasobów obliczeniowych pod względem pamięci i czasu przetwarzania.

\[ \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}. \]

Ponieważ metoda FDTD zasadniczo nie opiera się na żadnych założeniach dotyczących kierunku propagacji, struktury urządzenia ani zachowania pola, pojedyncza symulacja może uwzględnić szerokopasmową odpowiedź, dyfrakcję i interferencję, efekty bliskiego pola, geometrię o wymiarach mniejszych niż długość fali oraz przejściowe zjawiska elektromagnetyczne w ramach jednego modelu. Ta uniwersalność jest tym, co czyni metodę FDTD tak skuteczną, ale wiąże się również z pewnymi ograniczeniami: przestrzenną siatka musi uwzględniać najmniejszą fizyczną skalę, krok czasowy jest ściśle ograniczony ze względu na stabilność numeryczną, a trójwymiarowe symulacje szybko stają się wymagające pod względem pamięci i czasu obliczeń, nawet w przypadku geometrii, która wydaje się stosunkowo prosta.

2. Dlaczego nie należy stosować metody FDTD...

W większości praktycznych problemów optycznych metoda FDTD nie powinna być pierwszym wyborem. Jest ona zasobożerna obliczeniowo, a dla wielu systemów istnieją znacznie bardziej efektywne metody. Projektowanie układów optycznych, takich jak obiektywy do aparatów fotograficznych lub teleskopów, modelowanie rozbudowanych systemów obrazowania lub śledzenie promieni światła w makroskopowych układach optycznych, można znacznie lepiej realizować za pomocą metod śledzenia promieni. Podobnie, w przypadku płaskich struktur wielowarstwowych, takich jak cienkowarstwowe ogniwa słoneczne, układy OLED i powłoki optyczne, znacznie efektywniej stosuje się metody macierzowe. W takich przypadkach metoda FDTD generuje znaczne koszty obliczeniowe, nie przynosząc przy tym dodatkowych korzyści w postaci głębszego zrozumienia fizycznych zjawisk.

Większość metod symulacji elektromagnetycznych ostatecznie opiera się na równaniach Maxwella. Takie techniki, jak metody rozwiązywania problemów modalnych i metody macierzy przenoszenia, są opracowywane poprzez wprowadzenie upraszczających założeń dotyczących geometrii, wymiarowości lub struktury pola, aby przyspieszyć obliczenia. Te przybliżenia są niezwykle skuteczne, ale tylko wtedy, gdy spełnione są leżące u ich podstaw założenia.

Metoda FDTD znajduje się na przeciwnym krańcu tego spektrum. Jest niezwykle wydajna, szczególnie w przypadku złożonych trójwymiarowych układów, nieregularnych geometrii oraz problemów, w których pole nie można rozłożyć na proste mody lub warstwy – na przykład w przypadku fotonicznych kryształów, struktur o wymiarach mniejszych niż długość fali, silnie rozpraszających ośrodków oraz sprzężenia w bliskim polu. Ceną za tę elastyczność są wysokie wymagania obliczeniowe, zarówno pod względem pamięci, jak i czasu przetwarzania.

Kluczem jest wybór odpowiedniego narzędzia. W przypadkach, gdy pozwala na to model fizyczny, należy stosować uproszczone modele, a gdy te założenia przestają być aktualne, należy użyć metody FDTD. Jej zaletą jest elastyczność i pełna zgodność z zasadami fizyki, natomiast wadą jest wysoki koszt obliczeniowy.

3. Równania Maxwella jako podstawa metody FDTD.

W swojej istocie metoda FDTD jest bardzo prosta. Wykorzystuje się równania Maxwella, rozkłada się je na składowe w każdym kierunku przestrzennym, a następnie rozwiązuje się te równania bezpośrednio na siatce. Nie ma w tym nic magicznego – metoda jest koncepcyjnie prosta i w wielu aspektach prostsza niż bardziej złożone metody numeryczne, takie jak metody rozwiązywania równań dyfuzji i dryfu. Na początek rozważmy prawo Ampère'a. W liniowym ośrodku o przewodności \(\sigma\), przenikalności elektrycznej \(\epsilon\) i przenikalności magnetycznej \(\mu\), można je zapisać w następujący sposób:

\[ \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} \]

które można rozwinąć w następujący sposób:

Prawo Faradaya można zapisać w następującej postaci:

\[ -\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} \]

co można rozwinąć, aby uzyskać:

Równania wyróżnione na niebiesko i czerwono stanowią definicję metody FDTD i są jej podstawą. Nie ma nic bardziej skomplikowanego. Wszystko inne to kwestie implementacji: zastępowanie pochodnych różnicami skończonymi, wybór warunków brzegowych, definiowanie źródeł światła (wzbudzeń), pobieranie wyników za pomocą detektorów oraz tworzenie wydajnego (często równoległego) kodu, który symuluje zmiany pól w czasie. Te dodatkowe elementy sprawiają, że praktyczne rozwiązanie FDTD jest bardziej złożone, ale nie zmieniają one podstawowej prostoty samej metody.

4. Przykładowa symulacja 1D FDTD z wykorzystaniem zabawki.

Aby pokazać, jak prosta jest metoda FDTD – zarówno pod względem koncepcji, jak i implementacji – uprościmy teraz pełne równania do jednowymiarowego modelu, który można zaimplementować w kilkudziesięciu liniach kodu w MATLAB-ie lub Pythonie. Celem nie jest zmiana zasad fizycznych, lecz usunięcie zbędnej złożoności, aby struktura metody była w pełni zrozumiała.

W tym celu rozpatrujemy zmiany jedynie wzdłuż kierunku \(x\) i wyraźnie ustawiamy pozostałe pochodne przestrzenne na zero: \[ \frac{\partial}{\partial y} = 0, \qquad \frac{\partial}{\partial z} = 0. \]. W rezultacie otrzymujemy układ elektromagnetyczny o polaryzacji poprzecznej, w którym występuje pole elektryczne \(E_z(x,t)\) i pole magnetyczne \(H_y(x,t)\).

Przy tych założeniach równania Maxwella dotyczące wirowości upraszczają się do:

\[ \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}. \]

Oto podstawowe równania ciągłe dla przypadku jednowymiarowego. Jednak w tej postaci nie można ich bezpośrednio rozwiązać na komputerze. Aby to zrobić, zastępujemy pochodne skończonymi różnicami na dyskretnej siatce. To jest sedno tej metody: ogólny przegląd można znaleźć w skończone różnice.

Wprowadzając siatkę przestrzenną \(x_i = i\Delta x\) i kroki czasowe \(t^n = n\Delta t\) oraz stosując standardową, przesuniętą w czasie (metodę „leapfrog”), równania aktualizujące przyjmują postać:

\[ 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]. \]

To już kompletna metoda FDTD w najprostszej formie. Pole elektryczne jest aktualizowane na podstawie różnic przestrzennych w polu magnetycznym, a pole magnetyczne jest aktualizowane na podstawie różnic przestrzennych w polu elektrycznym, przy czym oba pola zmieniają się naprzemiennie w czasie.

Dla ośrodka niemagnetycznego \(\mu = \mu_0\) przenikalność elektryczna jest powiązana z współczynnikiem załamania wzorem:

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

Oznacza to, że zmiana współczynnika załamania światła jest wprowadzana poprzez prostą modyfikację wartości \(\epsilon\) w całej siatce. W poniższym przykładzie rozważamy ośrodok, który przechodzi od stanu \(n=1\) do stanu \(n=3\), tak że impuls wprowadzony z lewej strony częściowo odbija się, a częściowo przechodzi przez granicę.

W praktyce powyższe równania umieszcza się w pętli iteracyjnej, w której obliczenia wykonywane są dla kolejnych momentów czasu: pole magnetyczne jest aktualizowane na podstawie pola elektrycznego, stosowane są warunki brzegowe, a następnie pole elektryczne jest aktualizowane na podstawie pola magnetycznego, po czym proces się powtarza. Pola są przesunięte w czasie, dzięki czemu w prosty i bezpośredni sposób następują po sobie, a każda aktualizacja zależy tylko od lokalnych wartości w sąsiedztwie. Sposób działania tego procesu przedstawiono schematycznie w ??, i zasadniczo to wszystko, co komputer robi w każdym kroku czasowym.

5. Przykładowa symulacja jednowymiarowego modelu FDTD w Matlabie/Pythonie.

Tę samą, jednowymiarową metodę FDTD można zaimplementować w zaledwie kilkudziesięciu liniach kodu. Poniżej znajdują się przykładowe implementacje w językach MATLAB i Python. Obie symulują propagację impulsu Gaussa i pokazują efekt odbicia od granicy o różnym współczynniku załamania.

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. Uruchamianie symulacji.

Zachowanie symulacji przedstawiono na rysunkach ??–??. Początkowo impuls Gaussa rozchodzi się swobodnie przez obszar \(n=1\). Gdy dociera do obszaru o zmieniającym się współczynniku załamania, pole rozdziela się na dwie składowe: impuls odbity, który przemieszcza się z powrotem do obszaru o niższym współczynniku załamania, oraz impuls przechodzący, który wchodzi do obszaru o wyższym współczynniku załamania. Ponieważ współczynnik załamania wzrasta do wartości \(n=3\), prędkość fali przechodzącej maleje, a jej długość fali ulega skróceniu, natomiast składowa odbita zachowuje pierwotne właściwości rozchodzenia się. Jest to dokładnie takie zachowanie, jakiego można się spodziewać, i pokazuje, że nawet ta minimalna implementacja oddaje podstawowe zasady fizyczne dotyczące rozchodzenia się fali, odbicia i transmisji.

11. Co zobaczyć dalej? / Gdzie udać się następnym razem?

Przedstawione powyżej przykłady w MATLAB i Pythonie ilustrują, jak prosta jest w zasadzie metoda FDTD, ale są to jedynie uproszczone implementacje. W praktyce realistyczne symulacje optyczne wymagają dużych siatek obliczeniowych, długiego czasu obliczeń oraz starannego zarządzania warunkami brzegowymi i źródłami. Aby móc wykorzystać metodę FDTD do rozwiązywania praktycznych, realnych problemów, potrzebny jest solver, który jest wysoce zoptymalizowany, zrównoleglony i może działać zarówno na procesorach (CPU), jak i na procesorach graficznych (GPU).

Dokładnie to oferuje OghmaNano: kompletny trójwymiarowy solver oparty na metodzie FDTD, napisany w języku C i wykorzystujący akcelerację OpenCL, co pozwala na płynne uruchamianie tego samego modelu na wielordzeniowych procesorach lub kartach graficznych. Zamiast tworzyć własny solver od podstaw, można wykorzystać przedstawione tutaj pomysły i natychmiast zastosować je do realistycznych struktur fotonicznych, urządzeń i problemów związanych z optyką falową.

Aby rozpocząć, możesz zapoznać się ze stroną FDTD w OghmaNano, która zawiera ogólny opis możliwości, wrócić do sekcji co to jest FDTD lub pełne wyprowadzenie matematyczne, albo przejść bezpośrednio do praktycznych ustawień, wykorzystując źródła światła, warunki brzegowe i detektory.

Dostępna jest również szeroka gama gotowych przykładów, obejmujących m.in.: propagację w przestrzeni swobodnej, interferencję w układzie dwóch szczelin, fotoniczne przewodniki krystaliczne, siatki Bragg’a, rezonatory Fabry’ego-Pérota, pochyłe powierzchnie przewodników, struktury zwężające się, układy Macha-Zehndera, rezonatory pierścieniowe oraz powłoki antyrefleksyjne.

Krótko mówiąc: podstawowa metoda jest prosta, jak pokazano powyżej, ale OghmaNano zostało zaprojektowane tak, aby wykorzystać tę prostotę i rozwinąć ją w praktyczne, wydajne narzędzie do rzeczywistych symulacji optycznych.