Warunki brzegowe w metodzie FDTD

1. Wprowadzenie

W symulacjach opartych na metodzie różnic skończonych w dziedzinie czasu (FDTD), warunki brzegowe określają, jak pola elektromagnetyczne zachowują się na zewnętrznych granicach obszaru obliczeniowego. W terminach fizycznych, warunki te określają, co się dzieje, gdy fala dociera do krawędzi komórki symulacyjnej. Jest to jeden z najważniejszych aspektów modelu FDTD, ponieważ obszar numeryczny jest zawsze skończony, podczas gdy problem fizyczny często odnosi się do propagacji w otwartej przestrzeni, powtarzalnych struktur lub zamkniętych przestrzeniach. Dlatego też, rolą warunku brzegowego jest zapewnienie poprawnego matematycznego przedłużenia pól poza ograniczoną siatkę.

Jeśli granica zostanie wybrana nieodpowiednio, fale będą się odbijać od niej i wracać do siatki, zakłócając obliczenia i powodując powstawanie wyraźnych artefaktów w postaci fal stojących. W najgorszym przypadku energia zostaje uwięziona wewnątrz obszaru symulacji, a obliczone rozłożenie pola nie odzwierciedla już zamierzonego działania urządzenia. Jest to szczególnie ważne w otwartych strukturach fotonicznych, problemach związanych z rozpraszaniem, systemach emitujących promieniowanie oraz w obliczeniach odpowiedzi impulsowej, gdzie wychodzące fale muszą swobodnie opuszczać obszar. Z kolei w strukturach periodycznych pożądanym zachowaniem nie jest absorpcja, lecz powtarzalność, a w niektórych problemach testowych warto wyraźnie zdefiniować wartości pól na granicy. Dlatego wybór warunku brzegowego musi być zgodny z fizyką modelu.

Z matematycznego punktu widzenia, równania Maxwella opisujące wirowość pola elektromagnetycznego rozwiązuje się na skończonej siatce Yee. Dla obszaru izotropowego, w którym nie występują źródła pola, można je zapisać w następującej postaci:

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

Równania aktualizacji wartości wewnątrz obszaru obliczeń wynikają bezpośrednio z tych zależności, ale w przypadku komórek znajdujących się na zewnętrznej granicy, obszar obliczeń wykracza poza zdefiniowaną siatkę. Warunek brzegowy dostarcza brakujących informacji niezbędnych do zakończenia aktualizacji. W praktyce oznacza to, że każdej zewnętrznej powierzchni obszaru symulacji należy przypisać regułę numeryczną, która określa, czy pola styczne są ustalane, absorbowane, okresowo powtarzane, czy też tłumione za pomocą odpowiedniej warstwy absorbującej.

W programie OghmaNano warunki brzegowe konfiguruje się za pomocą przycisku „Warunki brzegowe” na pasku narzędzi „Optyka” (patrz Rysunek ??). Odpowiedni edytor (patrz Rysunek ??) umożliwia przypisanie różnych warunków do każdej zewnętrznej powierzchni obszaru FDTD. Dostępne typy warunków brzegowych to: Dirichlet, Mur ABC, okresowe i PML.

Sześć możliwych do wyboru granic odpowiada sześciu ścianom symulacyjnej przestrzeni kartezjańskiej. W domyślnej konfiguracji wyświetlania, używanej w programie OghmaNano, \(y_0\) to górna ściana obszaru symulacji, \(y_1\) to dolna ściana, \(x_0\) to lewa ściana, \(x_1\) to prawa ściana, \(z_0\) to ściana znajdująca się najbliżej użytkownika, a \(z_1\) to ściana znajdująca się najdalej od użytkownika. Układ przyporządkowania kierunków zilustrowano na Rysunek ??.

Edytor warunków brzegowych przypisuje niezależnie jedną regułę do każdej z tych sześciu ścian. Dzięki temu możliwe jest tworzenie mieszanych obszarów, na przykład z warunkami brzegowymi okresowymi w kierunkach bocznych i warunkami brzegowymi absorbującymi w kierunku propagacji, lub z warstwą PML na wszystkich ścianach, aby symulować promieniowanie w przestrzeń. Po wybraniu warstwy PML, dla danej ściany dostępny jest dodatkowy parametr głębokości. Głębokość ta jest wyrażana w komórkach siatki i określa, ile zewnętrznych warstw FDTD jest zarezerwowanych dla medium absorbującego.

2. Edytor warunków brzegowych

Edytor warunków brzegowych umożliwia oddzielne sterowanie poszczególnymi powierzchniami \(y_{\min}\), \(y_{\max}\), \(x_{\min}\), \(x_{\max}\), \(z_{\min}\) i \(z_{\max}\), które odpowiadają odpowiednio powierzchniom \(y_0\), \(y_1\), \(x_0\), \(x_1\), \(z_0\) i \(z_1\) zgodnie z konwencją kierunkową przedstawioną na rysunku Rysunek ??. Dla każdej powierzchni wybrana opcja określa sposób aktualizacji stycznych pól elektrycznego i magnetycznego na zewnętrznej krawędzi siatki.

3. Warunek brzegowy Dirichleta

Warunek brzegowy Dirichleta ustala wartość pola na granicy. W najczęściej stosowanym przypadku metody FDTD wartość ta wynosi zero, co oznacza, że granica wymusza...

\[ \mathbf{E}_{\mathrm{tan}} = 0 \qquad \text{on the selected face,} \]

W związku z tym granica zachowuje się jak idealnie odbijająca ściana dla pola elektrycznego, a każda fala wychodząca z obszaru symulacji jest odbijana z powrotem do tego obszaru.

3. Warunek brzegowy Mur ABC

Warunki brzegowe typu „absorbująca ściana” mają na celu zmniejszenie odbić poprzez przybliżenie zachowania fali rozchodzącej się na zewnątrz na krawędzi obszaru obliczeniowego. Koncepcja polega na spełnieniu jednokierunkowego równania falowego, które jest normalne do granicy, tak aby energia opuszczająca obszar obliczeniowy nadal rozchodziła się na zewnątrz, zamiast być odbijana.

Dla fali jednowymiarowej rozchodzącej się w kierunku \(+x\), dokładne równanie opisujące rozchodzącą się falę ma postać:

\[ \frac{\partial u}{\partial x} + \frac{1}{c}\frac{\partial u}{\partial t}=0. \]

Pierwszorzędowy warunek brzegowy pochłaniający Mur’a uzyskuje się poprzez dyskretyzację tego równania na siatce FDTD. Na lewym brzegu stosuje się standardową formę.

\[ u_0^{n+1} = u_1^{n} + \frac{c\Delta t-\Delta x}{c\Delta t+\Delta x} \left( u_1^{n+1}-u_0^{n} \right), \]

na przeciwległej stronie występuje analogiczne wyrażenie. Tutaj \(u_i^n\) oznacza pole w punkcie o współrzędnych przestrzennych \(i\) i czasowych \(n\), \(\Delta x\) to odstęp między węzłami siatki, a \(\Delta t\) to krok czasowy w metodzie FDTD.

Warunki brzegowe typu „ściana absorbująca” są obliczeniowo tanie i wymagają niewielkiej dodatkowej pamięci, co czyni je atrakcyjnymi w przypadku prostych symulacji lub obliczeń wstępnych. Mają na celu przybliżenie wychodzącej fali, tak aby energia mogła opuścić obszar symulacji, zamiast być odbijana od zewnętrznej granicy.

Ponieważ warunek Mur’a jest jedynie przybliżeniem, zawsze występuje pewien stopień odbicia. Dokładność maleje w przypadku skośnego padania, szerokopasmowych impulsów oraz złożonych, wielowymiarowych pól rozpraszania. Dlatego warunki Mur’a najlepiej sprawdzają się w prostych problemach propagacji fal, gdzie fale padają na granicę mniej więcej prostopadle do powierzchni i nie jest wymagana bardzo wysoka dokładność.

W przypadku wymagających symulacji optycznych, struktur emitujących promieniowanie lub sytuacji, w których fale padają na granicę pod różnymi kątami, należy zastosować warstwy idealnie dopasowane (PML). Granice warstw PML wprowadzają obszar pochłaniający, który znacznie skuteczniej i przy znacznie niższym współczynniku odbicia tłumi fale wychodzące, niż w przypadku warunków Mur.

4. Warunek brzegowy okresowy

Warunek brzegowy okresowy identyfikuje jedną ścianę obszaru symulacji z przeciwległą ścianą, tak aby pole wychodzące z jednej strony wchodziło ponownie przez drugą. Jest to odpowiednie, gdy fizyczna struktura powtarza się w nieskończoność w przestrzeni, a komórka obliczeniowa reprezentuje pojedynczą komórkę elementarną tej powtarzającej się geometrii.

W najprostszej formie okresowość narzuca...

\[ \mathbf{E}(x+L_x,y,z,t)=\mathbf{E}(x,y,z,t), \qquad \mathbf{H}(x+L_x,y,z,t)=\mathbf{H}(x,y,z,t), \]

Podobnie, w kierunkach \(y\) lub \(z\), gdy te granice są oznaczone jako okresowe. Numerycznie, wartości pól wymagane bezpośrednio za jedną ze stron obszaru pobierane są z przeciwległej strony. Oznacza to, że przeciwległe powierzchnie należy traktować jako sparowaną parę. Fizycznie oznacza to, że rozwiązanie odpowiada nie pojedynczemu obiektowi, lecz nieskończonej siatce kopii symulowanej komórki.

Warunkowe granice są powszechnie stosowane w przypadku struktur, które powtarzają się w przestrzeni, takich jak siatki dyfrakcyjne, kryształy fotoniczne, elementy metamateriałów oraz periodyczne geometrie światłowodów. Modelując tylko jeden element, można znacznie zmniejszyć nakład obliczeniowy, zachowując jednocześnie możliwość odtworzenia zachowania całej powtarzającej się struktury, pod warunkiem, że źródło i geometria są zgodne z założoną periodycznością.

W przypadku zastosowania warunków brzegowych okresowych, pola, które opuszczają jedną stronę obszaru symulacji, wracają przez przeciwległą stronę. W związku z tym energia nie jest absorbowana, lecz przepływa w obie strony wewnątrz komórki obliczeniowej. Warunków brzegowych okresowych nie należy zatem stosować do modelowania przestrzeni otwartej. Jeśli warunki brzegowe okresowe zostaną zastosowane na wszystkich powierzchniach obszaru symulacji, energia nie będzie mogła opuścić tego obszaru, a całkowita energia elektromagnetyczna może się z czasem kumulować. W wielu praktycznych symulacjach FDTD warunki brzegowe okresowe są zatem stosowane w kierunkach poprzecznych, natomiast warunki brzegowe absorbujące, takie jak PML, są stosowane w kierunku propagacji.

W najbardziej ogólnym przypadku, periodyczne pola mogą obejmować przesunięcie fazy między sąsiednimi komórkami. Zjawisko to opisuje warunek Blocha.

\[ \mathbf{E}(x+L_x,y,z,t)=\mathbf{E}(x,y,z,t)e^{ik_xL_x}, \]

Podobna zależność występuje również w przypadku pola magnetycznego. Standardowe warunki brzegowe okresowe, zastosowane w tym modelu, odpowiadają szczególnemu przypadkowi \(k_x=0\), co oznacza, że pola powtarzają się identycznie na przeciwległych ścianach obszaru symulacji.

5. Warunek brzegowy PML

Idealnie dopasowana warstwa (PML) to standardowa, wysokowydajna warstwa absorbująca, stosowana we współczesnych symulacjach metodą FDTD. Jej zadaniem jest tłumienie fal elektromagnetycznych wychodzących z obszaru symulacji przy minimalnym odbiciu, w szerokim zakresie częstotliwości i dla szerokiego zakresu kątów. Zamiast stosować pojedynczą formułę brzegową na zewnętrznej krawędzi, PML wprowadza sztuczne medium absorbujące, otaczające fizyczny obszar symulacji. Fale wchodzące do tej warstwy ulegają tłumieniu, bez występowania nagłych zmian impedancji na granicy.

Kluczową ideą jest dopasowanie warstwy PML do właściwości ośrodka w płaszczyźnie wejścia, dzięki czemu fala wchodząca nie napotyka nieciągłości, która sama w sobie generowałaby odbicie. Po wejściu do warstwy, pola są tłumione w sposób wykładniczy. W formie ciągłej można to interpretować za pomocą złożonego rozciągania współrzędnych, na przykład.

\[ x \;\rightarrow\; \int_0^x s_x(\xi)\,d\xi, \qquad s_x(\xi)=\kappa_x(\xi)+\frac{\sigma_x(\xi)}{j\omega\varepsilon_0}, \]

wraz z odpowiadającymi im współczynnikami rozciągania w pozostałych kierunkach. Tutaj \(\sigma_x\) to sztucznie wygenerowany profil przewodności, a \(\kappa_x\) to współczynnik skalowania stosowany w formułach PML z rozciągniętymi współrzędnymi. W efekcie, rozprzestrzeniające się pola są tłumione podczas przechodzenia przez warstwę, przy jednoczesnym zachowaniu w przybliżeniu braku odbić na granicy z rzeczywistą dziedziną.

W programie OghmaNano, wybranie PML dla danej powierzchni aktywuje tę warstwę pochłaniającą po tej stronie obszaru symulacji. Powiązany parametr głębokości PML określa liczbę komórek siatki przypisanych do tej warstwy. Grubsza warstwa PML zazwyczaj zapewnia silniejsze tłumienie, ponieważ pole ma większą odległość, na której może zanikać. Jeśli warstwa PML jest zbyt cienka, część energii może dotrzeć do zewnętrznej krawędzi i się odbić. Jeśli jest wystarczająco gruba, wychodzące pole jest tłumione do wartości pomijalnej, zanim dotrze do zewnętrznej krawędzi.

Wygodnym sposobem na zrozumienie zależności grubości od danego parametru jest zauważenie, że składowa pola wewnątrz materiału absorbującego maleje w przybliżeniu zgodnie z następującym wzorem:

\[ E(d)\sim E(0)\exp(-\alpha d), \]

gdzie \(d\) to odległość, jaką fala pokonuje wewnątrz warstwy PML, a \(\alpha\) to efektywny współczynnik tłumienia, który zależy od zmiany przewodności. Zwiększenie liczby komórek PML powoduje zatem zmniejszenie amplitudy sygnału resztkowego docierającego do zewnętrznej granicy obszaru obliczeniowego. W praktyce głębokość warstwy PML powinna być wystarczająco duża, aby odbicia od tylnej krawędzi warstwy były pomijalne w porównaniu z fizycznym sygnałem, który nas interesuje.

Metoda PML jest zazwyczaj preferowanym warunkiem brzegowym w otwartych symulacjach fotonicznych i elektromagnetycznych, obejmujących promieniowanie w światłowodach, rozpraszanie na pojedynczych obiektach, emisję przypominającą promieniowanie antenowe oraz propagację impulsów. Jest znacznie bardziej niezawodna niż proste warunki brzegowe absorbujące, gdy pola padają na granicę pod kątem lub zawierają szerokie spektrum częstotliwości. Z tego powodu metoda PML jest domyślnym wyborem w wielu modelach FDTD, w których oczekuje się, że fale będą opuszczać obszar obliczeniowy.

Nawet przy stosowaniu PML, przestrzeganie dobrych zasad modelowania pozostaje ważne. Warstwę pochłaniającą nie należy umieszczać zbyt blisko obszaru o silnie zanikającym lub wysoko rezonansowym polu, ponieważ pola bliskie, które się nie rozprzestrzeniają, mogą wchodzić w interakcję z materiałem pochłaniającym i zmieniać wynik symulacji. Ogólnie rzecz biorąc, fizyczną strukturę należy oddzielić od PML wystarczającą ilością wolnej przestrzeni, aby wychodząca fala mogła się swobodnie formować. W ten sposób stosowana, PML zapewnia najlepsze przybliżenie w metodzie FDTD otwartej, nieodbijającej domeny zewnętrznej.