Що таке FDTD?

1. Вступ

Метод скінченних різниць у часовій області (FDTD) безпосередньо розв’язує рівняння Максвелла, які описують поширення світла, у часовій області на дискретній просторовій сітці. На відміну від трасування променів, розв’язувачів мод або методів на основі матриць передачі, FDTD не робить жодних спрощуючих припущень щодо поширення або модальної структури. Натомість електричне та магнітне поля обчислюються явно в часі шляхом дискретизації рівнянь ротора – закону Ампера та закону Фарадея, щоб електромагнітні хвилі природним чином поширювалися, розсіювалися, інтерферували та дифрактували в межах області моделювання. Завдяки цій універсальності FDTD є надзвичайно точним і гнучким методом, але водночас потребує значних обчислювальних ресурсів з точки зору пам’яті та часу обчислень.

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

Оскільки метод FDTD практично не робить жодних припущень щодо напрямку поширення, структури пристрою чи поведінки поля, одна симуляція може охопити широкий спектр явищ, зокрема дифракцію та інтерференцію, ефекти ближнього поля, субхвильову геометрію та перехідні електромагнітні явища в межах єдиної моделі. Ця універсальність і робить метод FDTD таким потужним, але вона також має свою ціну: просторова сітка повинна враховувати найменший фізичний розмір, часовий крок суворо обмежений чисельною стабільністю, і тривимірні симуляції швидко стають ресурсомісткими як за обсягом пам’яті, так і за часом обчислень, навіть для геометрій, які здаються невеликими за розміром.

2. Чому не варто використовувати метод FDTD...

У більшості практичних задач оптики метод FDTD не повинен бути першим, до якого ви звертаєтесь. Він потребує значних обчислювальних ресурсів, і для багатьох систем існують набагато ефективніші підходи. Розробка оптичних систем, таких як камери або телескопи, моделювання масштабних систем візуалізації або відстеження шляху світла через макроскопічні оптичні пристрої набагато краще виконується за допомогою методів трасування променів. Так само, для плоских багатошарових структур, таких як тонкоплівкові сонячні батареї, OLED-стеки та оптичні покриття, набагато ефективніше використовувати методи матриці передачі. У цих випадках використання FDTD значно збільшує обчислювальні витрати, не надаючи при цьому додаткової фізичної інформації.

Більшість методів електромагнітного моделювання зрештою базуються на рівняннях Максвелла. Такі методи, як розв’язувачі мод і методи матриці передачі, розробляються шляхом введення спрощуючих припущень щодо геометрії, розмірності або структури поля з метою прискорення обчислень. Ці наближення є надзвичайно ефективними, але лише за умови, що їхні основні припущення є обґрунтованими.

Метод FDTD розташований на протилежному кінці цього спектру. Він надзвичайно потужний, особливо для складних тривимірних систем, нерегулярних геометричних форм і задач, у яких поле неможливо розкласти на прості складові або шари, наприклад, для фотонних кристалів, структур із розмірами, меншими за довжину хвилі, сильно розсіюючих середовищ і явищ поближнього поля. Ціною за таку гнучкість є обчислювальні витрати, як за обсягом пам’яті, так і за часом обчислень.

Головне – правильно підібрати інструмент. Використовуйте спрощені моделі, коли це дозволяють фізичні принципи, а метод FDTD – коли ці принципи перестають бути дійсними. Його перевага полягає у гнучкості та повноті фізичного опису; недолік – високі обчислювальні витрати.

3. Рівняння Максвелла, що лежать в основі методу FDTD.

У своїй основі метод FDTD є досить простим. Ви берете рівняння Максвелла, розкладаєте їх на складові для кожного просторового напрямку, а потім безпосередньо розв’язуєте ці рівняння на сітці. Тут немає нічого надприродного — метод концептуально простий і в багатьох аспектах простіший за більш складні чисельні методи, такі як розв’язувачі рівнянь переносу-дифузії. Для початку розкладемо закон Ампера. У лінійному середовищі з провідністю \(\sigma\), діелектричною проникністю \(\epsilon\) і магнітною проникністю \(\mu\) його можна записати так:

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

яке можна розширити до

Закон Фарадея формулюється наступним чином:

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

яке можна розширити, щоб отримати:

Рівняння, виділені синім і червоним кольорами, є визначенням методу FDTD і становлять його основу. Нічого складнішого за це не існує. Все інше – це реалізація: заміна похідних на скінченні різниці, вибір граничних умов, визначення джерел світла (збуджень), отримання результатів за допомогою детекторів і написання ефективного (часто паралельного) коду для моделювання зміни полів у часі. Ці додаткові компоненти роблять практичний розв’язувач FDTD складнішим, але вони не змінюють основну простоту базового методу.

4. Приклад іграшки, змодельованої за допомогою одномірного методу FDTD.

Щоб показати, наскільки простим є метод FDTD — як у плані концепції, так і в плані реалізації, — ми зараз зводимо повні рівняння до одномірної спрощеної моделі, яку можна реалізувати за допомогою кількох десятків рядків коду в MATLAB або Python. Мета тут полягає не в тому, щоб змінити фізичні принципи, а в тому, щоб усунути зайву складність, щоб структура методу була повністю зрозумілою.

Для цього ми розглядаємо зміни лише вздовж напрямку \(x\) і явно встановлюємо значення решти просторових похідних рівним нулю: \[ \frac{\partial}{\partial y} = 0, \qquad \frac{\partial}{\partial z} = 0. \]. Таким чином, у нас залишається поперечна електромагнітна система з електричним полем \(E_z(x,t)\) і магнітним полем \(H_y(x,t)\).

За цих умов рівняння Максвелла, що описують ротор, спрощуються до:

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

Це основні диференціальні рівняння для одномірного випадку. Однак, у такому вигляді їх неможливо безпосередньо розв’язати за допомогою комп’ютера. Для цього ми замінюємо похідні на скінченні різниці на дискретній сітці. В цьому полягає суть методу: загальний огляд можна знайти в скінченні різниці.

Вводячи просторову сітку \(x_i = i\Delta x\) та часові інтервали \(t^n = n\Delta t\) і використовуючи стандартну «шахову» схему, рівняння оновлення набувають такого вигляду:

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

Це вже повна версія методу FDTD у найпростішому вигляді. Електричне поле оновлюється на основі просторових різниць у магнітному полі, а магнітне поле оновлюється на основі просторових різниць в електричному полі, причому ці два поля чергуються в часі.

Для немагнітного середовища \(\mu = \mu_0\) діелектрична проникність пов’язана з показником заломлення наступним чином:

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

Це означає, що зміна показника заломлення досягається шляхом простої зміни значення \(\epsilon\) у межах заданої області. У наведеному нижче прикладі ми розглядаємо середовище, яке переходить від \(n=1\) до \(n=3\), так що імпульс, що поширюється зліва, частково відбивається і частково проходить крізь межу поділу.

На практиці наведені вище рівняння просто розміщують у циклі, що моделює рух у часі: магнітне поле оновлюється на основі електричного поля, застосовуються граничні умови, потім електричне поле оновлюється на основі магнітного поля, і процес повторюється. Поля змінюються в часі таким чином, щоб вони послідовно змінювали один одного в простій, явній схемі, причому кожне оновлення залежить лише від локальних сусідніх значень. Це схематично показано в ??, і по суті, це все, що робить комп’ютер на кожному кроці часу.

5. Простий приклад моделювання за допомогою методу FDTD в Matlab/Python.

Ту саму одномірну схему FDTD можна реалізувати лише в декількох десятках рядків коду. Нижче наведено еквівалентні реалізації в MATLAB і Python. Обидві схеми генерують гаусів імпульс і демонструють відбиття від різкої зміни показника заломлення.

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. Запуск моделювання.

Поведінка моделювання показана на ??–??. Спочатку гаусів імпульс поширюється безперешкодно через область \(n=1\). Коли він досягає межі зміни показника заломлення, поле розщеплюється на дві складові: відбитий імпульс, який рухається назад до області з нижчим показником заломлення, і пропущений імпульс, який проникає в область з вищим показником заломлення. Оскільки показник заломлення збільшується до \(n=3\), пропущена хвиля сповільнюється та стискається за довжиною хвилі, тоді як відбита складова зберігає початкові характеристики поширення. Це саме та фізична поведінка, яку ми очікували, і вона демонструє, що навіть ця мінімальна реалізація відображає основні фізичні принципи поширення, відбиття та проходження хвиль.

11. Куди піти далі?

Наведені вище приклади для MATLAB і Python демонструють, наскільки простим є метод FDTD в теорії, але це лише спрощені реалізації. На практиці для реалістичного моделювання оптичних процесів потрібні великі обчислювальні сітки, тривалий час обчислень і ретельне опрацювання граничних умов і джерел. Щоб застосовувати метод FDTD до практичних, реальних задач, потрібен вирішувач, який був би надзвичайно оптимізований, паралелізований і здатний працювати як на центральних процесорах (CPU), так і на графічних процесорах (GPU).

Саме це й пропонує OghmaNano: повноцінний тривимірний розв’язувач FDTD, написаний мовою C з використанням технології OpenCL для прискорення обчислень, що дозволяє безперешкодно запускати одну й ту саму модель на багатоядерних процесорах або графічних процесорах. Замість того, щоб створювати власний розв’язувач з нуля, ви можете використати представлені тут ідеї та негайно застосувати їх до реалістичних фотонних структур, пристроїв і задач хвильової оптики.

Щоб розпочати, ви можете переглянути сторінку FDTD в OghmaNano, щоб отримати загальне уявлення про можливості, повторно ознайомитися з інформацією про те, що таке FDTD, або з повним математичним обґрунтуванням, або перейти безпосередньо до практичної реалізації, використовуючи джерела світла, граничні умови та детектори.

Також доступний широкий спектр готових прикладів, зокрема: поширення в просторі, інтерференція на подвійній щілині, фотонні кристалічні хвилеводи, дифракційні ґратки Брегга, резонатори Фабрі–Перо, нахилені грані хвилеводів, конічні структури, пристрої Маха–Цендера, кільцеві резонатори та антиблікові покриття.

Коротко кажучи: основний метод простий, як показано вище, але OghmaNano розроблено таким чином, щоб використати цю простоту та розширити її, перетворивши на практичний, високоефективний інструмент для реального оптичного моделювання.