Метод FDTD — розв'язання рівнянь Максвелла на сітці

1. Навіщо розв'язувати рівняння Максвелла напряму?

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

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

Це робить FDTD надзвичайно універсальним: будь-яку геометрію, яку можна представити на сітці (довільно складні провідники, діелектрики, анізотропні чи дисперсійні матеріали), можна симулювати тим самим базовим алгоритмом. Ціна — обчислювальна вартість: сітки, достатньо дрібні для розв'язання деталей масштабу довжини хвилі у 3D, можуть вимагати мільярдів вузлів сітки та мільйонів кроків часу для електрично великих задач.

Відправна точка — два ротор-рівняння Максвелла в безджерельному, лінійному, ізотропному середовищі:

Це пов'язана пара диференціальних рівнянь у частинних похідних першого порядку: похідна за часом від H залежить від просторового ротора E, а похідна за часом від E залежить від просторового ротора H (плюс член втрат на провідність). FDTD наближає кожну похідну — за часом і за всіма трьома просторовими напрямками — центрованою скінченною різницею, а потім обчислює поле на наступному кроці часу.

2. Сітка Йі: розсування E та H у просторі

Ключова ідея Кейна Йі, опублікована в статті IEEE 1966 року, полягала в тому, щоб не зберігати E і H в одних і тих самих точках сітки. Натомість компоненти електричного та магнітного поля зсунуті на пів комірки сітки в кожному напрямку — схема, яку тепер називають коміркою Йі або решіткою Йі.

У 3D-комірці Йі кожна компонента поля E розташована в центрі ребра куба, а кожна компонента поля H — у центрі грані куба, спрямована по нормалі до цієї грані:

Це розсування — не косметичний вибір, а саме те, що змушує центровані просторові різниці E точно збігатися з розташуванням H, і навпаки. Центрована різниця E_y та E_z навколо точки H_x дає наближення другого порядку точності до (∇ × E)_x саме там, де розташовано H_x. У результаті закон Фарадея та закон Ампера виконуються з точністю другого порядку без будь-якої інтерполяції.

Розсування має й глибокий фізичний зміст: воно автоматично забезпечує виконання умов дивергенції ∇·B = 0 та ∇·D = ρ у кожній точці сітки в будь-який момент часу, оскільки дискретний оператор ротора на сітці Йі має ту саму властивість нульового простору, що й неперервний оператор ротора. Саме тому FDTD не потребує окремого кроку «очищення дивергенції», на відміну від багатьох інших розв'язувачів рівнянь у частинних похідних.

У двох вимірах (поширений випадок для навчання й для багатьох плоских задач) сітка Йі зводиться до двох незалежних поляризацій. Мода TM_z пов'язує E_z, H_x, H_y; мода TE_z пов'язує H_z, E_x, E_y. 2D-симуляція потребує оновлення лише трьох компонент поля замість шести — саме це зазвичай реалізує інтерактивна FDTD-демонстрація в браузері.

3. Рівняння оновлення: крокування «жабкою» в часі

Час дискретизується аналогічно: значення поля E обчислюються на цілих кроках часу (t = 0, Δt, 2Δt, ...), а значення поля H — на напівцілих кроках (t = Δt/2, 3Δt/2, ...). Це схема «жабки» (leapfrog): E та H чергуються, кожне оновлюється, використовуючи найсвіжіші значення іншого, які вже відомі з попереднього напівкроку.

Для найпростішого випадку — 1D-симуляції плоскої хвилі, що поширюється вздовж z, лише з компонентами E_x та H_y — рівняння оновлення зводяться до:

Кожне оновлення повністю явне: нове значення поля обчислюється простою зваженою сумою вже відомих значень без розв'язання системи рівнянь. Саме це робить FDTD швидким на кожному кроці часу й тривіально паралельним — оновлення кожної комірки сітки залежить лише від найближчих сусідів, що робить метод ідеальним для прискорення на GPU та декомпозиції області на кластерах.

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

4. Межа стійкості Куранта

Явна схема маршового просування в часі корисна лише тоді, коли вона чисельно стійка — якщо малі похибки не зростають експоненційно з кожним кроком. Для FDTD стійкість вимагає, щоб крок часу Δt був досить малим, аби інформація не могла поширюватися більш ніж на одну комірку сітки за крок, оскільки алгоритм обмінюється інформацією лише між сусідніми вузлами сітки Йі.

Це умова Куранта–Фрідріхса–Леві (CFL), і для FDTD вона набуває конкретної форми, відомої як межа Куранта:

де c — швидкість світла в середовищі. На практиці інженери застосовують запас безпеки — зазвичай обираючи Δt на рівні 90–99% від максимально допустимого значення (число Куранта S = c·Δt/Δx тримають трохи нижче 1 в 1D, або нижче 1/√розмірність у вищих вимірах).

Порушення межі Куранта не пробачається: навіть один крок часу вище порогу спричиняє необмежене зростання найвищочастотних мод сітки, і за кілька кроків уся симуляція розходиться у видимий чисельний шум або переповнення. Цей зв'язок між просторовою роздільною здатністю та максимальним стійким кроком часу означає, що подрібнення сітки для розв'язання дрібних геометричних деталей автоматично вимагає менших (численніших) кроків часу — головний чинник вартості в 3D-моделях FDTD із тонкими деталями.

5. Чисельна дисперсія та роздільна здатність сітки

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

Для 1D-сітки точне дисперсійне співвідношення дискретних рівнянь оновлення має вигляд:

що зводиться до ідеального ω = ck лише в границі Δx, Δt → 0. За скінченного кроку сітки хвилі з коротшою довжиною хвилі (більшим k відносно сітки) поширюються трохи повільніше за справжню швидкість світла, а хвилі, що йдуть діагонально через сітку, поширюються з дещо іншою швидкістю, ніж хвилі, вирівняні з осями. На великих відстанях поширення це накопичується у видиме розтягування імпульсу та фазову похибку.

Стандартне інженерне правило для утримання похибки дисперсії на прийнятно малому рівні (зазвичай кілька відсотків за фазою) таке:

де λ_min — найкоротша довжина хвилі в симуляції (найвища частота інтересу, або довжина хвилі в матеріалі з найвищим показником заломлення). Оскільки пам'ять і час обчислення масштабуються з кількістю комірок сітки — а в 3D це означає куб лінійної роздільної здатності — вибір кількості комірок на довжину хвилі зазвичай є найбільшим важелем впливу на вартість симуляції.

6. Введення джерел: жорсткі, м'які та TFSF

Симуляції потрібен спосіб вводити енергію в сітку. Найпростіший підхід, жорстке джерело (hard source), просто перезаписує значення поля в одній точці сітки на кожному кроці часу — наприклад, примусово задає E_x[джерело] = sin(ωt). Це легко реалізувати, але воно діє як ідеальний розсіювач: будь-яка хвиля, що пізніше повертається до цієї точки (відбиття від межі чи об'єкта), відбивається назад у сітку замість того, щоб пройти крізь неї, що забруднює результати паразитними повторними відбиттями.

М'яке джерело (soft source) натомість додає член джерела до наявного значення поля замість перезапису, що дозволяє вихідним і повернутим хвилям природно накладатися в цій точці без штучного розсіювання, властивого жорсткому джерелу. Це переважний метод для більшості неперервних та імпульсних збуджень.

Для задач розсіювання — обчислення того, як хвиля розсіюється на об'єкті, з чітким відокремленням розсіяного поля від величезного падаючого поля, — стандартна техніка це формулювання повне поле/розсіяне поле (total-field/scattered-field, TFSF). Сітку поділяють на внутрішню область «повного поля» (падаюча + розсіяна хвиля) і зовнішню область «розсіяного поля» (лише розсіяна хвиля), з'єднані тонкою віртуальною межею, де корекційні члени точно вводять падаюче поле. Це дозволяє плоскій хвилі чисто освітлювати ціль з одного напрямку, тоді як зовнішня область фіксує лише розсіяну відповідь об'єкта — необхідно для розрахунків ефективної площі розсіювання радара (RCS) та діаграм спрямованості у дальній зоні.

Часова форма джерела також важлива. Поширений вибір — гаусів модульований імпульс, який має гладкий, обмежений за смугою частотний спектр і уникає збудження паразитного високочастотного шуму сітки, який спричинило б розривне джерело у вигляді сходинки:

7. Поглинальні межі: ідеально узгоджений шар (PML)

Оскільки комп'ютерна сітка обов'язково скінченна, а більшість електромагнітних задач моделюють випромінювання в необмежений вільний простір, зовнішні краї сітки мають поводитися так, ніби простір триває нескінченно — поглинаючи вихідні хвилі без відбиття. Наївне обрізання (просто зупинити сітку) відбиває хвилі назад у область, як ідеальне дзеркало, спотворюючи всю симуляцію.

Сучасне рішення, запропоноване Жаном-П'єром Беранже 1994 року, — ідеально узгоджений шар (Perfectly Matched Layer, PML): штучна поглинальна область навколо розрахункової області, у якій хвильовий імпеданс математично узгоджений з вільним простором під будь-яким кутом падіння та на будь-якій частоті, тож у неперервній (нерозбитій на сітку) границі коефіцієнт відбиття точно дорівнює нулю незалежно від кута падіння. Усередині PML амплітуда поля експоненційно спадає завдяки штучному профілю провідності, що плавно наростає від нуля (на внутрішньому краю PML) до максимального значення (на зовнішньому, ідеально провідному завершенні):

Поліноміальне наростання (а не різкий стрибок до повної провідності) — це принципово: різкий розрив провідності сам по собі спричиняє відбиття, тому σ наростає достатньо плавно, щоб дискретизована сітка бачила майже неперервне узгодження імпедансу. Типовий PML має товщину 8–20 комірок і, за правильного налаштування, досягає коефіцієнтів відбиття нижче -60...-80 дБ — далеко за межею, потрібною для точних симуляцій дальнього поля й антен. Оскільки реальний PML реалізовано на дискретній сітці, а не в справді неперервному середовищі, невелике залишкове відбиття завжди зберігається, але воно нехтовно мале для практичних інженерних цілей.

8. Мінімальна 1D-реалізація FDTD

Ядро FDTD-розв'язувача напрочуд коротке. Ось повна 1D-симуляція гаусового імпульсу, що поширюється у вільному просторі, з простою втратною поглинальною межею, написана мовою JavaScript:

const N = 400;          // кількість комірок сітки
const dx = 1e-3;       // просторовий крок (м)
const c0 = 3e8;        // швидкість світла (м/с)
const dt = dx / (2 * c0); // межа Куранта, S = 0.5

let Ex = new Float64Array(N);
let Hy = new Float64Array(N);
const eps0 = 8.854e-12, mu0 = 1.2566e-6;

function step(t) {
  // оновити H з E (напівкрок "жабки")
  for (let k = 0; k < N - 1; k++) {
    Hy[k] += (dt / (mu0 * dx)) * (Ex[k + 1] - Ex[k]);
  }
  // оновити E з H (напівкрок "жабки")
  for (let k = 1; k < N; k++) {
    Ex[k] += (dt / (eps0 * dx)) * (Hy[k] - Hy[k - 1]);
  }
  // м'яке гаусове джерело в комірці 50
  const t0 = 40, tau = 12;
  Ex[50] += Math.exp(-((t - t0) ** 2) / (2 * tau * tau));
}

Виконання цього циклу протягом кількох сотень кроків часу з побудовою графіка Ex відносно позиції на сітці для кожного кадру показує, що гаусів імпульс розщеплюється на дві зустрічно-поширювані фронти хвиль, які рухаються назовні (дуже близько) зі швидкістю світла — саме те, що передбачають рівняння Максвелла, отримане без ручного запису розв'язку хвильового рівняння. Розширення того самого скелету до 2D і 3D, додавання матеріалів з різними ε та μ для кожної комірки та додавання PML на краях — саме так побудовано промислові FDTD-розв'язувачі, такі як MEEP, Lumerical FDTD та openEMS, лише з набагато більшою оптимізацією та значно більшою кількістю компонент поля.

9. Застосування: антени, фотоніка та радар

Оскільки FDTD обробляє довільну геометрію та довільні властивості матеріалів (включно з частотно-залежними, анізотропними та нелінійними середовищами, через розширення базових рівнянь оновлення), він став одним із двох домінантних методів в обчислювальній електродинаміці поряд із частотним методом скінченних елементів (FEM).

• Проєктування антен: прогнозування діаграм випромінювання, вхідного імпедансу та смуги пропускання антен телефонів, патч-антен і антенних решіток, вбудованих у реалістичний корпус пристрою.
• Фотоніка й нанофотоніка: симуляція поширення світла через хвилеводи, кільцеві резонатори, фотонні кристали та плазмонні наноструктури, де розмір елементів порівнянний з оптичною довжиною хвилі або менший за неї.
• Ефективна площа розсіювання радара (RCS): обчислення того, скільки енергії радара літак, корабель чи транспортний засіб розсіює назад до радара, за допомогою описаної вище техніки TFSF.
• Біоелектромагнетизм: моделювання того, як радіочастотна енергія мобільних телефонів чи котушок МРТ поглинається тканинами, виражене через питому швидкість поглинання (SAR).
• Мікрохвильові й радіочастотні схеми: вилучення S-параметрів фільтрів, відгалужувачів і роз'ємів безпосередньо з даних поля в часовій області через перетворення Фур'є записаних напруг і струмів портів.

Головна перевага FDTD над методами в частотній області полягає в тому, що одна широкосмугова симуляція в часовій області, за якою слідує перетворення Фур'є записаних полів, дає частотну характеристику для всієї смуги інтересу за один прогін — тоді як частотний розв'язувач необхідно перезапускати окремо для кожної частотної точки. Головний недолік — пам'ять і час обчислення для електрично великих 3D-задач, оскільки роздільна здатність сітки має розв'язувати найкоротшу довжину хвилі, присутню в усій області, навіть у регіонах далеко від будь-яких дрібних геометричних деталей — обмеження, що спонукало появу нерівномірних і адаптивних варіантів сітки.