Виведення методу FDTD

1. Вступ

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

Отже, ми починаємо з закону Ампера та закону Фарадея, представлених у диференціальній формі, розкладаємо їх на декартові компоненти, а потім адаптуємо до двовимірної оптичної задачі. Просторові похідні замінюються центральними скінченними різницями, часові похідні дискретизуються за допомогою стандартної схеми «стрибок-стрибок», а отримані вирази перегруповуються для отримання явних формул оновлення для компонентів електричного та магнітного полів.

Щоб зробити виклад більш зрозумілим, ми зводимо ці рівняння до двовимірного випадку, в якому \(\partial/\partial y = 0\) і магнітна провідність встановлюються в нуль, \(\sigma_m = 0\). Це значно спрощує алгебраїчні обчислення та наочно демонструє, як на практиці здійснюється зменшення розмірності. Однак у повній версії OghmaNano FDTD, електромагнітні рівняння розв’язуються у тривимірному просторі. Тому наведений тут двовимірний виклад має на меті бути чітким прикладом застосування методу, а не твердженням про те, що основний розв’язувач обмежений двовимірним простором.

На цій сторінці компоненти електричного поля обчислюються в цілі моменти часу \(n\Delta t\), а компоненти магнітного поля – у півцілі моменти часу \(\left(n+\tfrac{1}{2}\right)\Delta t\). Це стандартна схема «стрибкоподібної» сітки Йі, яка використовується в методі FDTD. На відміну від короткої оглядової сторінки, тут індексація вказана більш чітко, щоб виведення можна було сприймати як посібник з реалізації, а також як математичну довідку.

2. Початкові рівняння

Ми починаємо з рівнянь Максвелла для випадку лінійного середовища з електричною провідністю \(\sigma\), діелектричною проникністю \(\epsilon\), магнітною проникністю \(\mu\) і, для загальності, магнітною провідністю \(\sigma_m\):

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

У детермінованій формі завитки мають вигляд

\[ \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix}, \qquad \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix}. \]

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

3. Тензорна форма рівнянь для обчислення ротора.

Розкладання закону Ампера на декартові компоненти дає:

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = \frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}. \]

Узагальнена форма закону Фарадея має вигляд:

\[ -\sigma_m H_x - \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z}, \]

\[ -\sigma_m H_y - \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x}, \]

\[ -\sigma_m H_z - \mu \frac{\partial H_z}{\partial t} = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y}. \]

На цьому етапі ще не зроблено жодних припущень щодо розмірності. Наведені вище вирази є просто складовими частинами повних тривимірних рівнянь ротора.

4. Двовимірне зведення.

Щоб спростити процес отримання розв’язку, ми зараз зводимо задачу до двох вимірів, припускаючи, що немає жодних змін у напрямку \(y\), отже, \(\partial / \partial y = 0\). Цей крок робиться виключно для спрощення алгебраїчних обчислень і полегшення розуміння структури рівнянь оновлення. Цей самий метод можна безпосередньо розширити до трьох вимірів, і на практиці OghmaNano розв’язує повну тривимірну форму рівнянь Максвелла.

Виходячи з цього припущення, рівняння для електричного поля спрощуються до:

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}, \]

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

Якщо ми також встановимо магнітну провідність в нуль, \(\sigma_m = 0\), то рівняння магнітного поля набудуть вигляду:

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}, \]

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}, \]

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

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

5. Сітка Yee, індексування та позначення.

Для чисельного розв’язання наведених вище рівнянь як простір, так і час дискретизуються. Нехай крок обчислювальної сітки у напрямках \(x\) і \(z\) дорівнює відповідно \(\Delta x\) і \(\Delta z\), а крок часу – \(\Delta t\). Ми використовуємо цілочисельні індекси \(i\) і \(k\) для напрямків \(x\) і \(z\):

\[ x_i = i \Delta x, \qquad z_k = k \Delta z. \]

Час позначається цілим числом \(n\), щоб

\[ t^n = n \Delta t. \]

У схемі Йі електричні та магнітні поля не зберігаються в одних і тих же точках у просторі чи часі. Натомість вони розташовані в шаховому порядку. Компоненти електричного поля обчислюються в цілі моменти часу \(n\Delta t\), тоді як компоненти магнітного поля обчислюються в півцілих моментах часу \(\left(n+\tfrac{1}{2}\right)\Delta t\). Саме це є основою «стрибкоподібної» структури FDTD.

У просторовому вимірі компоненти поля також зміщуються таким чином, щоб кожна центрована скінченно-різницева похідна була природно визначена між сусідніми точками сітки. У повному програмному коді кожна компонента поля має свою пару просторових індексів, наприклад \(E_x^n(i,k)\), \(H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2})\) і так далі. Щоб зберегти зрозумілість обчислень, деякі повторювані «спостережні» індекси не завжди явно вказуються в наступних кроках. Однак інтерпретація завжди однакова:

\(i\) визначає \(x\)-напрямок.
\(k\) визначає \(z\)-напрямок.
\(n\) визначає дискретні часові рівні для цілочисельних значень напруженості електричного поля, і
\(n+\tfrac{1}{2}\) показує часові рівні, що відповідають зміні магнітного поля на півтону.

Наприклад, централізована скінченно-різницева апроксимація для \(\partial H_y / \partial z\) має вигляд:

\[ \frac{\partial H_y}{\partial z} \approx \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}, \]

і аналогічно, центрована апроксимація для \(\partial H_z / \partial x\) є

\[ \frac{\partial H_z}{\partial x} \approx \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Завдяки чіткому визначенню цієї познаки, наведені нижче обчислення можна однозначно інтерпретувати як належний алгоритм оновлення сітки Йі.

6. Виведення рівнянь для обчислення зміни електричного поля

Електричне поле обчислюється в цілі моменти часу \(n\Delta t\), тоді як магнітне поле обчислюється в моментах часу, що є половиною цілого числа \(\left(n+\tfrac{1}{2}\right)\Delta t\). Для отримання рівнянь оновлення, похідна за часом апроксимується за допомогою центральної скінченної різниці, а провідний член обчислюється на основі середнього значення старих і нових електричних полів. Це забезпечує чисельно стабільне трапецієподібне наближення провідності.

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

6.1 Оновіть рівняння для \(E_x\).

Почніть з...

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}. \]

Оцініть значення коефіцієнта провідності, використовуючи середнє значення \(E_x^{n+1}\) і \(E_x^n\), і оцініть похідну за часом, застосовуючи метод скінченних різниць із центральною апроксимацією:

\[ \sigma \frac{E_x^{n+1}(i,k) + E_x^n(i,k)}{2} + \epsilon \frac{E_x^{n+1}(i,k) - E_x^n(i,k)}{\Delta t} = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Зберіть усі члени, що містять \(E_x^{n+1}\), зліва:

\[ \left( \frac{\sigma}{2} + \frac{\epsilon}{\Delta t} \right) E_x^{n+1}(i,k) = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} + \left( \frac{\epsilon}{\Delta t} - \frac{\sigma}{2} \right) E_x^n(i,k). \]

Отже.

\[ E_x^{n+1}(i,k) = \frac{ \left( \dfrac{\epsilon}{\Delta t} - \dfrac{\sigma}{2} \right) E_x^n(i,k) - \dfrac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} }{ \dfrac{\epsilon}{\Delta t} + \dfrac{\sigma}{2} }. \]

Якщо помножити чисельник і знаменник на \(2\Delta t\), отримаємо знайому форму коефіцієнта:

\[ E_x^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_x^n(i,k) - \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

6.2 Оновіть рівняння для \(E_y\).

Почніть з...

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}. \]

Дискретизуйте обидва члени похідної на сітці Йі:

\[ \sigma \frac{E_y^{n+1}(i,k) + E_y^n(i,k)}{2} + \epsilon \frac{E_y^{n+1}(i,k) - E_y^n(i,k)}{\Delta t} = - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Перестановка дає змогу

\[ E_y^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_y^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \left[ - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} \right]. \]

Це єдиний компонент електричного поля в двовимірному зведенні, який містить внески від двох окремих похідних магнітного поля, і тому найлегше допустити в ньому помилку в коді. Зокрема, відносний знак між членами похідної \(x\) і \(z\) безпосередньо випливає з розкладання ротора і не повинен змінюватися.

6.3 Оновіть рівняння для \(E_z\).

Почніть з...

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

Дискретизація дає.

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

Отже, звідси випливає, таким чином.

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

На цьому етапі отримано дискретні рівняння для обчислення зміни електричного поля у їхній стандартній формі з урахуванням загасання. Наступний крок – вивести відповідні рівняння для обчислення зміни магнітного поля, використовуючи проміжні значення часу.

7. Виведення рівнянь для обчислення зміни магнітного поля

Магнітне поле оновлюється з інтервалами, що відповідають півцілим значенням часу. У випадку \(\sigma_m = 0\) рівняння для магнітного поля не містять члена, що описує загасання, тому формули оновлення набувають простішої форми прямого обчислення. Ця частина обчислення є алгебраїчно коротшою, ніж у випадку з електричним полем, оскільки тут немає трапецієподібного члена провідності, який потрібно було б враховувати.

7.1 Оновіть рівняння для \(H_x\).

Почніть з...

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}. \]

Дискретизація за часом між проміжними кроками \(n-\tfrac{1}{2}\) і \(n+\tfrac{1}{2}\) та в просторі з використанням центральної різниці в напрямку \(z\) дає:

\[ \mu \frac{ H_x^{n+\frac{1}{2}}(i,k+\frac{1}{2}) - H_x^{n-\frac{1}{2}}(i,k+\frac{1}{2}) }{\Delta t} = \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}. \]

У компактній формі запису з використанням половинного індексу, яка часто використовується в математичних викладках, той самий результат можна записати так:

\[ \mu \frac{ H_x^{n+\frac{1}{2}} - H_x^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Отже.

\[ H_x^{n+\frac{1}{2}} = H_x^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

7.2 Оновіть рівняння для \(H_y\).

Почніть з...

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

Дискретизація дає.

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

Знову ж таки, у скороченій формі індексу, що містить лише половину необхідних значень,

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

Отже, звідси випливає, таким чином.

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

7.3 Оновіть рівняння для \(H_z\).

Почніть з...

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

Дискретизація дає змогу

\[ \mu \frac{ H_z^{n+\frac{1}{2}}(i+\frac{1}{2},k) - H_z^{n-\frac{1}{2}}(i+\frac{1}{2},k) }{\Delta t} = - \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

Іншими словами, у стислому вигляді:

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

Отже.

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

Після виконання всіх цих замін, стандартні рівняння оновлення для методу FDTD, що використовує схему «стрибок-стрибок», застосовуються безпосередньо.

8. Остаточні рівняння для оновлення.

Узагальнюючи отримані результати, можна сказати, що двовимірні рівняння оновлення для методу FDTD мають такий вигляд:

\[ H_x^{n+\frac{1}{2}}(i,k+\tfrac{1}{2}) = H_x^{n-\frac{1}{2}}(i,k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}, \]

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

\[ H_z^{n+\frac{1}{2}}(i+\tfrac{1}{2},k) = H_z^{n-\frac{1}{2}}(i+\tfrac{1}{2},k) - \frac{\Delta t}{\mu} \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

У такому вигляді викладена суть схеми «стриб-скок», яка стає цілком зрозумілою. Електричне поле в момент часу \(n+1\) визначається на основі магнітного поля в момент часу \(n+\tfrac{1}{2}\), тоді як магнітне поле в момент часу \(n+\tfrac{1}{2}\) визначається на основі електричного поля в момент часу \(n\). Таким чином, поля змінюються в часі в певній послідовності.

9. Зауваження щодо впровадження.

У практичному коді коефіцієнти, що множать попередні значення поля та члени, пов’язані з ротором, часто обчислюються заздалегідь і зберігаються, особливо коли значення \(\epsilon\), \(\mu\) або \(\sigma\) змінюються від комірки до комірки. Наприклад, оновлення електричного поля зазвичай відбувається за такою схемою:

\[ E^{n+1} = C_a E^n + C_b (\nabla \times H), \]

де локальні коефіцієнти \(C_a\) і \(C_b\) кодують параметри матеріалу та крок часу. Це спрощує реалізацію та дозволяє уникнути повторного обчислення одних і тих самих факторів на кожному кроці часу.

Варто також наголосити, що наведена вище формула описує лише основну схему оновлення FDTD. Для проведення практичного моделювання додатково потрібні джерела, граничні умови, поглинальні шари, такі як ідеально узгоджені шари (PML), і крок часу, який задовольняє умові стабільності Куранта. Жодне з цих доповнень не змінює фундаментальну структуру, отриману тут, але кожне з них модифікує спосіб застосування рівнянь оновлення до конкретних точок сітки.

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