• English
  • Русский (Russsian)

Это старая версия документа!


Здесь находится вводное описание метода Finite-Difference Time-Domain (FDTD).

FDTD

FDTD является одним из наиболее популярных методов для численного решения уравнений Максвелла. Базовый алгоритм этого метода был предложен Кейном Йи в 1966 г. Однако имеющуюся популярность FDTD приобрел только в 90х гг. прошлого столетия, когда он стал основным для моделирования самых разных оптических приложений. Вот некоторые причины его популярности:

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

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

Для реализации указанных возможностей требуется написание больших параллельных программ. Имеются коммерческие решения, однако они не могут удовлетворить исследователей, которые хотели бы понимать, как работает программа, а не иметь дело с черным ящиком. Это явилось причиной разработки открытой библиотеки EMTL (Electromagnetic Template Library).

Алгоритм Йи

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

В отсутствие свободных зарядов уравнения Максвелла в дифференциальной форме имеют вид:

tex:\frac{\partial \vec{B}}{\partial t} = -\nabla \times \vec{E} - \vec{M},

tex:\frac{\partial \vec{D}}{\partial t} = \nabla \times \vec{H} - \vec{J},

tex:\nabla \cdot \vec{D} = 0,

tex:\nabla \cdot \vec{B} = 0,

где tex:\vec E и tex:\vec H - напряженности электрического и магнитного полей, tex:\vec D - электрическая, tex:\vec B - магнитная индукция, tex:\vec J - плотность электрического тока, tex:\vec M - ее магнитный аналог, который мы включили для общности некоторых используемых далее выражений.

В целях простоты изложения мы ограничим рассмотрение только изотропными, недисперсными и линейными среды. Тогда вектора tex:\vec E и tex:\vec D, tex:\vec H и tex:\vec B связаны следующими соотношениями:

tex:\vec{D}=\varepsilon \vec{E}=\varepsilon_r \varepsilon_0 \vec{E}, \quad \vec{B}=\mu \vec{H}=\mu_r \mu_0 \vec{H},

где tex:\varepsilon, tex:\varepsilon_r, tex:\mu, tex:\mu_r - абсолютные и относительные диэлектрическая и магнитная проницаемости среды, tex:\varepsilon_0 и tex:\mu_0 - электрическая и магнитная постоянные.

Плотность электрического и магнитного токов представим в виде:

tex:\vec{J}=\vec{J}_{source}+\sigma \vec{E}, \quad \vec{M}=\vec{M}_{source}+\sigma^* \vec{H},

где tex:\vec{J}_{source} - плотность тока свободных источников, tex:\sigma - электропроводность, tex:\vec{M}_{source} и tex:\sigma^* - их магнитные аналоги.

Подставив эти соотношения в первую пару уравнений Максвелла, получим

tex:\frac{\partial \vec{H}}{\partial t} = -\frac{1}{\mu} \nabla \times \vec{E} - \frac{1}{\mu}\left(\vec{M}_{source} + \sigma^* \vec{H}\right),

tex:\frac{\partial \vec{E}}{\partial t} = \frac{1}{\varepsilon}\nabla \times \vec{H} - \frac{1}{\varepsilon}\left(\vec{J}_{source} + \sigma \vec{E}\right).

Согласно алгоритму Йи, уравнения Максвелла дискретизуются с использованием центрального разностного отношения для приближения пространственной и временной производных. Для этого сеточные узлы, в которых хранятся компоненты электрического и магнитного полей, смещены по отношению друг к другу на половину сеточного шага по каждой из пространственных переменных. В результате те узлы, которые соответствуют компонентам полей tex:\vec{E} расположены таким образом, что каждая компонента tex:\vec E окружена четырьмя компонентами tex:\vec H, и наоборот.


Положения сеточных узлов для компонент электрического и магнитного полей на сетке Йи.

Согласно алгоритму Йи соответствующие компонентам tex:\vec E и tex:\vec H узлы сдвинуты относительно друг друга по времени на половину временного шага (в качестве примера это продемонстрировано ниже для одномерной сетки Йи). Для расчета значений поля tex:\vec E на tex:n+1/2-ом временном шаге используются значения поля tex:\vec H на tex:n-ом. Аналогичным образом значения поля tex:\vec H на tex:n+1-ом шаге рассчитываются с использованием значений поля tex:\vec E на tex:n+1/2-ом шаге. Эта процедура продолжается до тех пор, пока расчет не будет закончен.


Развертка одномерной сетки Йи по времени: узлы, соответствующие электрическому и магнитному полям, сдвинуты друг относительно по пространству и по времени на половину шага.

Прежде чем получить конкретные разностные уравнения, используемые в алгоритме Йи, введем следующие обозначения. Будем ставить в соответствие каждому сеточному узлу три целых числа tex:i, tex:j и tex:k, определяющих положение этого узла в пространстве как

tex:(i,j,k)=(i\Delta x,j\Delta y,k\Delta z), где tex:\Delta x, tex:\Delta y и tex:\Delta z суть сеточные шаги по соответствующим направлениям. Будем обозначать произвольную сеточную функцию tex:u как

tex:u(i \Delta x, j \Delta y, k \Delta z, n \Delta t)=u_{i,j,k}^n,

где tex:\Delta t - шаг по времени, а tex:n - текущий временной шаг.

Как уже упоминалось, в алгоритме Йи для аппроксимации присутствующих в уравнениях Максвелла производных используется центральное разностное отношение. Так, для производной по координате tex:x в точке tex:(i \Delta x, j \Delta y, k \Delta z) в момент времени tex:n \Delta t имеем

tex:\frac{\partial u}{\partial x}(i \Delta x, j \Delta y, k \Delta z, n \Delta t)=\frac{u_{i+1/2,j,k}^n-u_{i-1/2,j,k}^n}{\Delta x}+O\left[(\Delta x)^2\right].

Аналогично для производной по времени tex:t

tex:\frac{\partial u}{\partial t}(i \Delta x, j \Delta y, k \Delta z, n \Delta t)=\frac{u_{i,j,k}^{n+1/2}-u_{i,j,k}^{n-1/2}}{\Delta t}+O\left[(\Delta t)^2\right].

Применяя эти приближения для производных, участвующих в проекции закона Ампера на ось tex:x

tex:\frac{\partial E_x}{\partial t}=\frac{1}{\varepsilon} \left[ \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z} - \left(J_{source,x}+\sigma E_x\right)\right],

получаем

tex:\varepsilon_{i,j+1/2, k+1/2}\frac{E_x|_{i,j+1/2, k+1/2}^{n+1/2}-E_x|_{i,j+1/2, k+1/2}^{n-1/2}}{\Delta t}=

tex:\frac{H_z|_{i,j+1, k+1/2}^n-H_z|_{i,j, k+1/2}^n}{\Delta y} -\frac{H_y|_{i,j+1/2, k+1}^n-H_y|_{i,j+1/2, k}^n}{\Delta z} -J_x|_{i,j+1/2, k+1/2}^n - \sigma_{i,j+1/2, k+1/2}E_x|_{i,j+1/2, k+1/2}^n.

Стоящие в правой части переменные берутся на временном шаге tex:n, включая поле tex:E_x. Поскольку для значения tex:E_x в момент времени tex:n данных на сетке нет, нужно для него пользоваться каким-либо приближением. Хорошим выбором является усреднение по соседним временным слоям:

tex:E_x|_{i,j+1/2, k+1/2}^n=\frac{E_x|_{i,j+1/2, k+1/2}^{n+1/2}+E_x|_{i,j+1/2, k+1/2}^{n-1/2}}{2}.

Положим tex:\Delta=\Delta x=\Delta y=\Delta z. Тогда мы можем явно выразить tex:E_x на tex:n+1/2 шагу:

tex:E_x|_{i,j+1/2, k+1/2}^{n+1/2}= C_{a,E_x}|_{i,j+1/2, k+1/2} E_x|_{i,j+1/2, k+1/2}^{n-1/2} +

tex:C_{a,E_x}|_{i,j+1/2, k+1/2}
(
H_z|_{i,j+1, k+1/2}^n-H_z|_{i,j, k+1/2}^n + H_y|_{i,j+1/2, k}^n-H_y|_{i,j+1/2, k+1}^n - J_x|_{i,j+1/2, k+1/2}^n \Delta
),

где tex:C_{a}|_{i,j,k}=\left( 1-\frac{\sigma_{i,j,k}\Delta t}{2 \varepsilon_{i,j,k}} \right)\left( 1+\frac{\sigma_{i,j,k}\Delta t}{2 \varepsilon_{i,j,k}} \right)^{-1}, tex:C_{b}|_{i,j,k}=\left( \frac{\Delta t}{\varepsilon_{i,j,k} \Delta} \right)\left( 1+\frac{\sigma_{i,j,k}\Delta t}{2 \varepsilon_{i,j,k}} \right)^{-1}.

Мы получили разностное уравнение, соответствующее проекции закона Ампера на ось tex:x, которое вместе с пятью оставшимися аналогичными разностными уравнениями и формируют алгоритм Йи.

Можно показать, что в случае непроводящей среды (tex:\sigma=0) и при отсутствии источников тока tex:\vec{J} два оставшиеся незадействованными нами уравнения Максвелла выполняются автоматически, а именно, дивергенция полей tex:\vec{E} и tex:\vec{H} всегда равна нулю. При наличии источников тока справедливость двух оставшихся уравнений не вытекает непосредственно из разностной схемы, но с хорошей точностью подтверждается в численном эксперименте.

Выбор значений tex:\Delta x, tex:\Delta y, tex:\Delta z обуславливается геометрией задачи и спектральным составом излучения. Можно дать следующую рекомендацию: на характерный размер объекта (радиус шариков, толщину экрана и т. п.) должно приходится не менее нескольких сеточных шагов, а на характерную длину волны - от десяти и больше. Величина значения tex:\Delta t ограничена сверху условием Куранта

tex:\Delta t<\frac{1}{c\sqrt{\displaystyle{\frac{1}{(\Delta x)}^2+\frac{1}{(\Delta y)^2}+\frac{1}{(\Delta z)^2}}}},

выполнение которого необходимо для устойчивости разностной схемы.

В этом параграфе мы рассмотрели реализацию алгоритма Йи для случая, когда среда характеризуется не зависящими от частоты скалярными значениями tex:\varepsilon, tex:\mu, tex:\sigma, tex:\sigma^*. Численное моделирование анизотропных, дисперсных и нелинейных сред требует модификации этого алгоритма и применения вспомогательных разностных уравнений.

Отметим также, что численная схема FDTD не предполагает возможности табличного задания зависимости диэлектрической проницаемости от частоты. Как моделировать в FDTD дисперсные материалы, можно прочитать в разделе Фитинг диэлектрической проницаемости.

Численный эксперимент FDTD

Перечислим «основных участников» численного эксперимента FDTD.

  • Материальные тела, оптические свойства которых мы исследуем.
  • Источник электромагнитной волны. Самый простой способ задания источника заключается в задании временной зависимости величины tex:J_{source} в уравнениях Максвелла. Такой тип источника обычно используется при моделировании диполей. Для генерации плоской волны более удобен другой тип источника, реализуемый с помощью метода полного и рассеянного поля (total-field / scattered-field method).
  • Детекторы, снимающие поля на сетке в течение всего численного эксперимента. Детекторам не отвечают какие-то реальные тела в пространстве, под ними подразумевается лишь то, что мы записываем значения полей в каких-то точках внутри вычислительного объема в файл. По окончании расчета на основании этого файла можно будет восстановить ход численного эксперимента. Детекторы могут размещаться в произвольном месте вычислительного объема и не привязаны к сетке. Значения полей на них получаются путем интерполяции по соседним сеточным узлам.

Обычный сценарий численного эксперимента FDTD выглядит так:

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

Поясним некоторые использованные нами понятия.

Применение FDTD

FDTD применяется в самых разнообразных приложениях: для расчета рассеяния падающей волны на различных структурах, для моделирования излучения источника в заданном электромагнитном окружении, для изучения оптических свойcтв резонаторов и волноводов. Ниже мы изложим, какие величины при этом измеряются.

Предварительные замечания

Напомним, что любое решение уравнений Максвелла tex:{\bf F}({\bf r},t) (tex:\bf F это либо tex:\bf E, либо tex:\bf H) в случае линейных сред и при отсутствии свободных зарядов и источников тока может быть представлено в виде суперпозиции гармонических полей:

tex:{\bf F}({\bf r},t)={\bf A}({\bf r},\omega)\cos{\omega t}+{\bf B}({\bf r},\omega)\sin{\omega t}.

Удобно рассматривать вектор tex:{\bf F} как действительную часть от комплексного вектора tex:{\bf F_c}={\bf C}\exp(-i\omega t), где tex:{\bf C}={\bf A}+i{\bf B}:

tex:{\bf F}={\rm Re}({\bf F_c}).

Комплексная зависимость tex:\exp(-i\omega t) не несет никакого физического смысла и вводится лишь для удобства.

Вектор плотности потока электромагнитной энергии (вектор Пойнтинга) равен tex:{\bf P} = {\bf E} \times {\bf H}. В монохроматической волне его величина быстро осциллирует, чего большинство приборов отследить не в силах. Вместо этого измеряют среднее вектора Пойнтинга по времени tex:{\bf S}:

tex:{\bf S}=<{\bf P}>=\frac{1}{\tau}\int_t^{t+T}{\bf P}(\tau) d\tau,

где временной интервал tex:T существенно больше частоты tex:1/\omega.

Можно показать, что для монохроматической волны

tex:{\bf S}=\frac{1}{2}{\rm Re}({\bf E_c}\times{\bf H^*_c}).

Абсолютная величина вектора tex:{\bf S} называется интенсивностью света tex:I.

Рассеяние на теле

Рассмотрим задачу о рассеянии падающей волны на заданном теле. Поле вне тела можно представить в виде суммы падающей и отраженной волн:

tex:{\bf E}={\bf E_i}+{\bf E_{s}}, \qquad
{\bf H}={\bf H_i}+{\bf H_{s}}.

В задаче о рассеянии интересуются количеством энергии рассеянной волны, ее угловой зависимостью, а также количеством поглощаемой телом энергии.

Окружим тело воображаемой замкнутой поверхностью tex:A. Средняя скорость, с которой энергия притекает в ограниченный ею объем, равна

tex:W_{abs}=-\int{\bf S}\cdot{\bf n} dA,

где tex:{\bf n} - нормаль к поверхности.

Если tex:W_{abs}>0, энергия поглощается внутри объема. Если внешняя среда является непоглощающей, то tex:W_{abs} есть энергия, поглощаемая телом.

В нашем случае среднее по времени от вектора Пойнтинга tex:{\bf S} может быть записано в виде суммы трех слагаемых (нижний индекс tex:c для комплексных векторов tex:{\bf E} и tex:{\bf H} мы опускаем)

tex:{\bf S}=\frac{1}{2}{\rm Re}({\bf E}\times{\bf H^*})={\bf S_i}+{\bf S_{sca}}+{\bf S_{ext}},

где каждое слагаемое по отдельности равно

tex:{\bf S_i}=\frac{1}{2}{\rm Re}({\bf E_i}\times{\bf H_i^*}), \qquad
{\bf S_{sca}}=\frac{1}{2}{\rm Re}({\bf E_{s}}\times{\bf H_{s}^*}), \qquad
{\bf S_{ext}}=\frac{1}{2}{\rm Re}({\bf E_i}\times{\bf H_{s}^*}+{\bf E_{s}}\times{\bf H_i^*}).

Наличие tex:{\bf S_{ext}} есть следствие интерференции падающей и отраженной волн.

Учитывая только что сказанное, для tex:W_{abs} получаем

tex:W_{abs}=W_i-W_{sca}+W_{ext},

tex:W_i=-\int{\bf S_i}\cdot{\bf n} dA, \qquad
W_{sca}=\int{\bf S_{sca}}\cdot{\bf n} dA, \qquad
W_{ext}=-\int{\bf S_{ext}}\cdot{\bf n} dA.

В случае непоглощающей внешней среды tex:W_i=0, и

tex:W_{ext}=W_{abs}+W_{sca}.

Поскольку рассматриваемые величины линейно зависит от интенсивности падающей волны tex:I, имеет смысл рассматривать отношения

tex:C_{ext}=\frac{W_{ext}}{I}, \qquad
C_{abs}=\frac{W_{abs}}{I}, \qquad
C_{sca}=\frac{W_{sca}}{I},

имеющее размерность площади, и называемые сечениями экстинции, поглощения и рассеяния

Можно определить эффективности экстинции, рассеяния и поглощения

tex:Q_{ext}=\frac{C_{ext}}{G}, \qquad
Q_{abs}=\frac{C_{abs}}{G}, \qquad
Q_{sca}=\frac{C_{sca}}{G},

где tex:G - площадь проекции тела на плоскость перпендикулярную направлению падающей волны (для сферы радиуса tex:r имеем tex:G=\pi r^2). На длинах волн больших размеров тела, тело может рассеивать и поглощать больше света, чем падает на его геометрическое сечение. Соответствующие эффективности при этом будут больше единцы.

Матрица амплитуд рассеяния

Помимо энергии рассеянной волны часто интересуются ее угловым распределением, которое характеризуется матрицей амплитуд рассеяния. Вводится она следующим образом.

Пусть на тело падает монохроматическая линейно поляризованная волна.

Направим ось tex:z вдоль направления падения. Выберем какую-либо точку внутри тела в качестве начала координат tex:O, а оси tex:x, tex:y направим произвольно. Нашей системе координат соответствуют ортонормированный базис векторов tex:{\bf e_x}, tex:{\bf e_y}, tex:{\bf e_z}.

Направление рассеяния tex:{\bf e_r} вместе с направлением падения tex:{\bf e_z} определяют плоскость рассеяния, которая однозначно задается азимутальным углом tex:\phi, если только tex:{\bf e_r} не параллельно оси tex:z. В этих двух случаях (tex:{\bf e_r}=\pm{\bf e_z}) любая плоскость параллельная оси tex:z может считаться плоскостью рассеяния.

Вектор tex:{\bf E_i} падающего электрического поля, лежащий в плоскости tex:xy, удобно разложить на две компоненты, одна из которых параллельна плоскости рассеяния, а другая перпендикулярна ей

tex:{\bf E_i}=\left( E_{0\parallel}{\bf e_{\parallel{}i}}+E_{0\perp}{\bf e_{\perp{}i}} \right)
\exp(i\vec{k}\vec{z}-i\omega{}t)=
E_{\parallel{}i}{\bf e_{\parallel{}i}}+E_{\perp{}i}{\bf e_{\perp{}i}}

Вектора ортонормированного базиса

tex:{\bf e_{\perp{}i}}=\sin\phi{\bf e_x}-\cos\phi{\bf e_y},

tex:{\bf e_{\parallel{}i}}=\cos\phi{\bf e_x}+\sin\phi{\bf e_y}

образуют правую тройку с вектором tex:{\bf e_z}:

tex:{\bf e_{\perp{}i}}\times{\bf e_{\parallel{}i}}={\bf e_z}

Кроме того, имеем

tex:{\bf e_{\perp{}i}}=-{\bf e_{\phi}}, \qquad {\bf e_{\parallel{}i}}=\sin\theta{\bf e_r}+\cos\theta{\bf e_{\theta}},

где tex:{\bf e_r}, tex:{\bf e_{\theta}}, tex:{\bf e_{\phi}} образуют базис, связанный с полярной системой координат tex:(r, \theta, \phi).

Если компоненты падающего поля в направлениях x и y обозначить через tex:E_{xi} и tex:E_{yi}, то

tex:E_{\parallel{}i}=\cos\phi{}E_{xi}+\sin\phi{}E_{yi},

tex:E_{\perp{}i}=\sin\phi{}E_{xi}-\cos\phi{}E_{yi}.

Известно, что на больших расстояниях от тела (tex:{\bf kr}\gg{}1) (в так называемой дальней зоне) рассеянное поле приближенно является поперечным (tex:{\bf e_r} \cdot {\bf E_s} \sim 0), и, может быть записано как

tex:{\bf E_s}=E_{\parallel{}s}{\bf e_{\parallel{}s}}+E_{\perp{}s}{\bf e_{\perp{}s}},

где

tex:{\bf e_{\parallel{}s}}={\bf e_{\theta}}, \qquad
{\bf e_{\perp{}s}}=-{\bf e_{\phi}}, \qquad
{\bf e_{\perp{}s}} \times {\bf e_{\parallel{}s}} = {\bf e_r}.

Базисный вектор tex:{\bf e_{\parallel{}s}} параллелен, а tex:{\bf e_{\perp{}s}} перпендикулярен плоскости рассеяния. Заметим, однако, что tex:{\bf E_s} и tex:{\bf E_i} определены в разных системах базисных векторов. Тем не менее, вследствие линейности уравнений Максвелла, связь между ними можно представить в матричном виде

,

где комплексные числа tex:S_j (tex:j=1,2,3,4) составляют амплитудную матрицу рассеяния, зависящую от угла рассеяния tex:\theta и азимута tex:\phi.

Экспериментальное измерение матричных элементов tex:S_j затруднительно, поэтому обычно меряют элементы матрицы рассеяния, связывающей так называемые параметры Стокса для падающей и рассеянной волн. Это действительная матрица tex:4\times{}4, из 16 элементов которой независимы только семь, которые в свою очередь могут быть выписаны через абсолютные значения tex:\left|S_j\right| (j=1,2,3,4) и разницы в фазе между tex:S_j. Желающие узнать об этом более подробно, могут ознакомиться с главой 3.3 книги 1).

Коэффициенты прохождения и отражения для плоскопериодических структур

Радиационное затухание источника

Приложения

Программа EMTL использовалась для решения различных оптических задач на многих кластерных вычислительных системах. Ниже представлен лист избранных публикаций:

  • A. Deinega, I. Valuev, B. Potapkin and Yu. Lozovik, «Minimizing light reflection from dielectric textured surfaces,» JOSA A 28, 770 (2011).

Смотрите также диссертацию Алексея Дейнеги.

1) C. F. Bohren and D. R. Huffman: Absorption and Scattering of Light by Small Particles, Wiley-Interscience, New York (1983)
 
/home/kintechlab/fdtd.kintechlab.com/docs/data/attic/ru/fdtd.1329089722.txt.gz · Последние изменения: 2012/02/13 03:35 — deinega     Наверх