Метод конечных разностей во временной области

FDTD относится к общему классу сеточных методов решения дифференциальных уравнений. Базовый алгоритм метода был впервые предложен Кейном Йи (Калифорнийский университет) в 1966 г. в статье «Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media» журнала «IEEE Transactions on Antennas and Propagation»[1]. Однако, название «Finite-difference time-domain» и аббревиатура FDTD были даны методу Алленом Тафловом (Северо-западный университет, штат Иллинойс).

В первоначальном узком смысле под FDTD подразумевалось использование базового алгоритма Йи для численного решения уравнений Максвелла. В современном более широком смысле FDTD включает в себя множество самых разнообразных возможностей: моделирование сред с дисперсными и нелинейными свойствами, применение различных типов сеток (помимо первично предложенной прямоугольной сетки Йи), использование методов постпроцессорной обработки результатов и т. д.

Примерно с 1990 г. метод конечных разностей стал основным для моделирования самых разных оптических приложений. Он может быть с успехом применен для решения широкого спектра задач: от моделирования сверхдлинных электромагнитных волн в геофизике (включая процессы в ионосфере) и микроволн (например для изучения сигнатурной радиолокации, расчёта характеристик антенн, разработки беспроводных устройств связи, в том числе цифровых) до решения задач в оптическом диапазоне (фотонные кристаллы, наноплазмоника, солитоны и биофотоника). К 2006 г. число публикаций, посвященных FDTD, достигло двух тысяч.

В настоящее время существует порядка 30 коммерческих программ FDTD, а также проекты с открытым исходным кодом (в числе которых несколько русских).

В уравнениях Максвелла изменение электрического поля E (частная производная) зависит от распределения в пространстве магнитного поля H (ротор). Аналогично, изменение поля H зависит от распределения в пространстве поля Е.

Поля в ячейке сетки FDTD. Из таких ячеек составляется пространственная трёхмерная сетка Йи

На этом наблюдении основан алгоритм Йи. Сетки для полей E и H смещены по отношению друг к другу на половину шага дискретизации времени и по каждой из пространственных переменных. Конечно-разностные уравнения позволяют определить поля E и H на данном временном шаге на основании известных значений полей на предыдущем.

При заданных начальных условиях алгоритм Йи дает эволюционное решение во времени от начала отсчета с заданным временным шагом.

Аналогичная (разделённая) сетка используется при решении задач гидродинамики (для давления и поля скорости).

Как и в любом другом разностном методе, в FDTD существует проблема неточного отображения границы тела на вычислительную сетку. Любая кривая поверхность, разделяющая соседние среды и геометрически не согласованная с сеткой, будет искажаться эффектом «лестничного приближения». Для решения данной проблемы можно использовать дополнительную сетку с большим разрешением в тех областях пространства, где расположены тела со сложной геометрической структурой[2]. Также можно видоизменять разностные уравнения в узлах сетки, находящихся вблизи границы между соседними телами[3]. Менее затратным методом является введение эффективной диэлектрической проницаемости вблизи границы между телами (subpixel smoothing) [4][5].

Численная схема FDTD не предполагает возможности табличного задания зависимости диэлектрической проницаемости от частоты. Однако, её можно представить в виде апроксимации (фитинга) членами Дебая, Друде, Лоренца или Лоренца с поглощением. Такая аппроксимация не обязательно имеет физический смысл, и может быть получена численно, например с помощью программы [6].

Генерация плоской волны в двумерном счетном пространстве, имеющем поглощающие граничные условия PML

Для того чтобы ограничить объем сетки, в FDTD нужны особые поглощающие граничные условия, которые моделируют уход электромагнитной волны на бесконечность. Для этого используются поглощающие граничные условия Мура или Ляо[7], или идеально согласованные слои (Perfect Matched Layers, PML). Условия Мура или Ляо намного проще, чем PML. Тем не менее, PML — строго говоря, являющиеся поглощающей приграничной областью, а не граничным условием как таковым — позволяют получить на порядки меньшие по величине коэффициенты отражения от границы.

Понятие идеально согласованных слоев (PML) было введено Жаном Пьером Беренже в статье журнала «The Journal of Computational Physics» в 1994 г.[8] Идея PML Беренже основывалась на разбиении исходных полей E и H на две компоненты, для каждой из которых должны решаться свои уравнения. Впоследствии были предложены усовершенствованные формулировки PML эквивалентные первоначальной формулировке Беренже. Так, в одноосном PML (Uniaxial PML) используется анизотропный поглощающий материал, что позволяет не вводить дополнительные переменные и остаться в рамках исходных уравнений Максвелла[9]. Однако одноосный PML, как и PML в формулировке Беренже, не удобны тем, что в них отсутствует поглощение затухающих волн, что не позволяет помещать PML близко к рассеивающим телам. Этого недостатка лишен оборотный PML (Convolutional PML), основанный на аналитическом продолжении уравнений Максвелла в комплексную плоскость таким образом, что их решение экспоненциально затухает[10]. CPML также удобнее в ограничении бесконечных проводящих и дисперсных сред. Помимо этого математическая формулировка CPML обладает большей наглядностью и доступностью для понимания.

В некоторых случаях использование PML приводит к расходимости расчета FDTD. Эту проблему можно устранить путём помещения за PML дополнительной поглощающей стенки[11].

Ход расчёта FDTD выглядит следующим образом:

• Задается счетная область, разрешение сетки и граничные условия. Граничные условия могут быть поглощающими или периодическими. Последние применяются для моделирования нормального падения плоской волны на периодическую структуру. Схема FDTD для моделирования наклонного падения требует периодических условий со сдвигом по времени, которые могут быть реализованы с помощью разных методов [12][13][14].
• Внутри счетной области помещаются материальные тела с заданными оптическими свойствами (диэлектрическая проницаемость и магнитная проводимость).
• Задается источник. Самый простой способ задания источника заключается в задании временной зависимости плотности тока J в уравнении Ампера. Такой тип источника обычно используется при моделировании диполей. Для генерации плоской волны более удобен другой тип источника, реализуемый с помощью метода полного и рассеянного поля (Total Field / Scattered Field).
• Источник генерирует конечную во времени электромагнитную волну, спектральный состав которой должен покрывать интересующий диапазон частот. Далее, волна падает на тела, перерассеивается на них, и, при наличии поглощающих граничных условий, через какое-то время уходит из счетной области. История распространения волны сохраняется.
• С помощью преобразования Фурье записанные значения полей переводятся в частотное представление. Далее, обрабатывая их (например, интегрируя поток энергии поля через какую-либо поверхность), можно получить оптические характеристики рассматриваемой структуры тел. Используя метод преобразования ближнего поля в дальнее (Near to Far Transformation), можно получить значения полей за пределами счетной области на основании эволюции поля внутри счетной области[15].

Как и любой другой численный метод, FDTD имеет свои достоинства и недостатки.

Достоинства:

• FDTD — это простой и интуитивно понятный метод.
• Поскольку FDTD работает во временной области, он позволяет получить результат для широкого спектра длин волн за один расчет. Это может быть полезно при решении задач, в которых не известны резонансные частоты или в случае моделирования широкополосных сигналов.
• FDTD позволяет создавать анимированные изображения распространения волны в моделируемом объеме.
• FDTD удобен при задании анизотропных, дисперсных и нелинейных сред.
• Метод позволяет непосредственно моделировать краевые эффекты и эффекты экранирования, причем поля внутри и вне экрана могут быть рассчитаны как напрямую, так и нет.

Недостатки:

• Величина шага дискретизации по пространству должна быть значительно меньше исследуемых длин волн и типичных размеров исследуемой структуры. В некоторых случаях (инверсные опалы с маленькими перегородками между шариками) это может потребовать сеток с маленьким шагом, что означает большие затраты памяти и большое время расчета.
• FDTD рассчитывает поля внутри счетной области. Если требуется найти поле на большом расстоянии от источника, то необходимо увеличение счетной области и времени расчета. Существуют модификации метода для нахождения поля на удалении, но они требуют постобработки.

• Уравнения Максвелла
• Электродинамика
• Разностная схема
• Численные методы