Различия

Здесь показаны различия между выбранной ревизией и текущей версией данной страницы.

ru:fdtd [2012/07/03 20:40]
deinega [Приложения]
ru:fdtd [2013/08/08 21:27] (текущий)
deinega [Приложения]
Строка 1: Строка 1:
  
 +Здесь находится вводное описание метода Finite-Difference Time-Domain (FDTD).
 +
 +====== FDTD (Finite-Difference Time-Domain) ======
 +
 +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>
 +
 +<tex>
 +\frac{\partial \vec{D}}{\partial t} = \nabla \times \vec{H} - \vec{J},
 +</tex>
 +
 +<tex>
 +\nabla \cdot \vec{D} = 0,
 +</tex>
 +
 +<tex>
 +\nabla \cdot \vec{B} = 0,
 +</tex>
 +
 +где <tex>\vec E</tex> и <tex>\vec H</tex> - напряженности электрического и магнитного полей, <tex>\vec D</tex> - электрическая, <tex>\vec B</tex> - магнитная индукция, <tex>\vec J</tex> - плотность электрического тока, <tex>\vec M</tex> - ее магнитный аналог, который мы включили для общности некоторых используемых далее выражений.
 +
 +В целях простоты изложения мы ограничим рассмотрение только изотропными, недисперсными и линейными среды. Тогда вектора <tex>\vec E</tex> и <tex>\vec D</tex>, <tex>\vec H</tex> и <tex>\vec B</tex> связаны следующими соотношениями:
 +
 +<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>
 +
 +где <tex>\varepsilon</tex>, <tex>\varepsilon_r</tex>, <tex>\mu</tex>, <tex>\mu_r</tex> - абсолютные и относительные диэлектрическая и магнитная проницаемости среды, <tex>\varepsilon_0</tex> и <tex>\mu_0</tex> - электрическая и магнитная постоянные.
 +
 +Плотность электрического и магнитного токов представим в виде:
 +
 +<tex>
 +\vec{J}=\vec{J}_{source}+\sigma \vec{E}, \quad \vec{M}=\vec{M}_{source}+\sigma^* \vec{H},
 +</tex>
 +
 +где <tex>\vec{J}_{source}</tex> - плотность тока свободных источников, <tex>\sigma</tex> - электропроводность, <tex>\vec{M}_{source}</tex> и <tex>\sigma^*</tex> - их магнитные аналоги.
 +
 +Подставив эти соотношения в первую пару уравнений Максвелла, получим
 +
 +<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>
 +
 +<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>
 +
 +Согласно алгоритму Йи, уравнения Максвелла дискретизуются с использованием центрального разностного отношения для приближения пространственной и временной производных. Для этого сеточные узлы, в которых хранятся компоненты электрического и магнитного полей, смещены по отношению друг к другу на половину сеточного шага по каждой из пространственных переменных. В результате те узлы, которые соответствуют компонентам полей <tex>\vec{E}</tex> расположены таким образом, что каждая компонента <tex>\vec E</tex> окружена четырьмя компонентами <tex>\vec H</tex>, и наоборот. 
 +
 +<html><center></html>
 +{{:fdtd:fdtd_yee.png?250}}\\
 +//Положения сеточных узлов для компонент электрического и магнитного полей на сетке Йи.//\\
 +<html></center></html>
 +
 +Согласно алгоритму Йи соответствующие компонентам <tex>\vec E</tex> и <tex>\vec H</tex> узлы сдвинуты относительно друг друга по времени на половину временного шага (в качестве примера это продемонстрировано ниже для одномерной сетки Йи). Для расчета значений поля <tex>\vec E</tex> на <tex>n+1/2</tex>-ом временном шаге используются значения поля <tex>\vec H</tex> на <tex>n</tex>-ом. Аналогичным образом значения поля <tex>\vec H</tex> на <tex>n+1</tex>-ом шаге рассчитываются с использованием значений поля <tex>\vec E</tex> на <tex>n+1/2</tex>-ом шаге. Эта процедура продолжается до тех пор, пока расчет не будет закончен.
 +
 +<html><center></html>
 +{{:fdtd:fdtd_time.png?250}}\\
 +//Развертка одномерной сетки Йи по времени: узлы, соответствующие электрическому и магнитному полям, сдвинуты друг относительно по пространству и по времени на половину шага.//\\
 +<html></center></html>
 +
 +Прежде чем получить конкретные разностные уравнения, используемые в алгоритме Йи, введем следующие обозначения. Будем ставить в соответствие каждому сеточному узлу три целых числа <tex>i</tex>, <tex>j</tex> и <tex>k</tex>, определяющих положение этого узла в пространстве как
 +
 +<tex>
 +(i,j,k)=(i\Delta x,j\Delta y,k\Delta z),
 +</tex>
 +где <tex>\Delta x</tex>, <tex>\Delta y</tex> и <tex>\Delta z</tex> суть сеточные шаги по соответствующим направлениям. Будем обозначать произвольную сеточную функцию <tex>u</tex> как
 +
 +<tex>
 +u(i \Delta x, j \Delta y, k \Delta z, n \Delta t)=u_{i,j,k}^n,
 +</tex>
 +
 +где <tex>\Delta t</tex> - шаг по времени, а <tex>n</tex> - текущий временной шаг.
 +
 +Как уже упоминалось, в алгоритме Йи для аппроксимации присутствующих в уравнениях Максвелла производных используется центральное разностное отношение. Так, для производной по координате <tex>x</tex> в точке <tex>(i \Delta x, j \Delta y, k \Delta z)</tex> в момент времени <tex>n \Delta t</tex> имеем
 +
 +<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>
 +
 +Аналогично для производной по времени <tex>t</tex>
 +
 +<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>
 +
 +Применяя эти приближения для производных, участвующих в проекции закона Ампера на ось <tex>x</tex>
 +
 +<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>
 +
 +получаем
 +
 +<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>
 +
 +<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>
 +
 +Стоящие в правой части переменные берутся на временном шаге <tex>n</tex>, включая поле <tex>E_x</tex>. Поскольку для значения <tex>E_x</tex> в момент времени <tex>n</tex> данных на сетке нет, нужно для него пользоваться каким-либо приближением. Хорошим выбором является усреднение по соседним временным слоям:
 +
 +<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>
 +
 +Положим <tex>\Delta=\Delta x=\Delta y=\Delta z</tex>. Тогда мы можем явно выразить <tex>E_x</tex> на <tex>n+1/2</tex> шагу:
 +
 +<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>
 +
 +<tex>
 +C_{b,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>
 +
 +где 
 +<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>
 +<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>
 +
 +Мы получили разностное уравнение, соответствующее проекции закона Ампера на ось <tex>x</tex>, которое вместе с пятью оставшимися аналогичными разностными уравнениями и формируют алгоритм Йи.
 +
 +Можно показать, что в случае непроводящей среды (<tex>\sigma=0</tex>) и при отсутствии источников тока <tex>\vec{J}</tex> два оставшиеся незадействованными нами уравнения Максвелла выполняются автоматически, а именно, дивергенция полей <tex>\vec{E}</tex> и <tex>\vec{H}</tex> всегда равна нулю. При наличии источников тока справедливость двух оставшихся уравнений не вытекает непосредственно из разностной схемы, но с хорошей точностью подтверждается в численном эксперименте.
 +
 +Выбор значений <tex>\Delta x</tex>, <tex>\Delta y</tex>, <tex>\Delta z</tex> обуславливается геометрией задачи и спектральным составом излучения. Можно дать следующую рекомендацию: на характерный размер объекта (радиус шариков, толщину экрана и т. п.) должно приходится не менее нескольких сеточных шагов, а на характерную длину волны - от десяти и больше. Величина значения <tex>\Delta t</tex> ограничена сверху условием Куранта
 +
 +<tex>
 +\Delta t<\frac{1}{c\sqrt{\displaystyle{\frac{1}{(\Delta x)}^2+\frac{1}{(\Delta y)^2}+\frac{1}{(\Delta z)^2}}}},
 +</tex>
 +
 +выполнение которого необходимо для устойчивости разностной схемы.
 +
 +В этом параграфе мы рассмотрели реализацию алгоритма Йи для случая, когда среда характеризуется не зависящими от частоты скалярными значениями <tex>\varepsilon</tex>, <tex>\mu</tex>, <tex>\sigma</tex>, <tex>\sigma^*</tex>. Численное моделирование анизотропных, дисперсных и нелинейных сред требует модификации этого алгоритма и применения вспомогательных разностных уравнений.
 +
 +Отметим также, что численная схема FDTD не предполагает возможности табличного задания зависимости диэлектрической проницаемости от частоты. Как моделировать в FDTD дисперсные материалы, можно прочитать в разделе [[fitting]].
 +
 +===== Численный эксперимент FDTD =====
 +
 +Перечислим "основных участников" численного эксперимента FDTD.
 +
 +  *Материальные тела, оптические свойства которых мы исследуем.
 +  *Источник электромагнитной волны. Самый простой способ задания источника заключается в задании временной зависимости величины <tex>J_{source}</tex> в уравнениях Максвелла. Такой тип источника обычно используется при моделировании диполей. Для генерации плоской волны более удобен другой тип источника, реализуемый с помощью метода полного и рассеянного поля (total field / scattered field method).
 +  *Детекторы, снимающие поля на сетке в течение всего численного эксперимента. Детекторам не отвечают какие-то реальные тела в пространстве, под ними подразумевается лишь то, что мы записываем значения полей в каких-то точках внутри вычислительного объема в файл. По окончании расчета на основании этого файла можно будет восстановить ход численного эксперимента. Детекторы могут размещаться в произвольном месте вычислительного объема и не привязаны к сетке. Значения полей на них получаются путем интерполяции по соседним сеточным узлам.
 +
 +Обычный сценарий численного эксперимента FDTD выглядит так: 
 +  * Внутри вычислительного объема, задаваемого размером используемой сетки, помещаются материальные тела, источник и детекторы (детекторам не отвечают какие-то реальные тела в пространстве, под ними подразумевается лишь то, что мы записываем значения полей в каких-то точках в файл). На границе вычислительного объема должны быть заданы граничные условия, обычно, периодические или поглощающие.
 +  * Источник генерирует конечную во времени электромагнитную волну, спектральный состав которой должен покрывать интересующий нас диапазон частот. Далее, волна падает на тела, перерассеивается на них, и, при наличии поглощающих граничных условий, через какое-то время уходит из вычислительного объема. История распространения волны фиксируется детекторами.
 +  * С помощью преобразования Фурье записанные на детекторах значения полей переводятся в частотное представление. Далее, обрабатывая их (например, интегрируя поток энергии поля через какую-либо поверхность), мы получаем интересующие нас оптические характеристики рассматриваемой структуры тел.
 +
 +Поясним некоторые использованные нами понятия. 
 +
 +===== Метод полного и рассеянного поля =====
 +
 +Этот метод используется для моделирования бесконечно удаленного источника плоской волны. Этот метод основан на линейности уравнений Максвелла и вытекающего из них принципа суперпозиции. А именно, мы предполагаем, что измеряемое (полное) поле <tex>\vec E_{total}</tex> и <tex>\vec H_{total}</tex> может быть представлено в виде суммы
 +
 +<tex>
 +\vec E_{total} = \vec E_{inc} + \vec E_{scat}, \quad \vec H_{total} = \vec H_{inc} + \vec H_{scat}.
 +</tex>
 +
 +<tex>\vec E_{inc}</tex>, <tex>\vec H_{inc}</tex> соотвествует полю падающей волны, которая предполагается известной во всех точках пространства в любой момент времени. Это та волна, которая распространялась бы в пространстве, если бы в нем не существовало никаких тел. <tex>\vec E_{scat}</tex>, <tex>\vec H_{scat}</tex> соответствует рассеянной волне, представляющей собой результат взаимодействия падающей волны с телами. Значение рассеянной волны заранее неизвестно.
 +
 +Разностные уравнения FDTD могут независимо применяться как для полного поля, так и для падающего или рассеянного полей. Это позволяет разбить вычислительный объем на две области: область полного поля, в которой рассчитывается полное поле, и область рассеянного поля, в которой рассчитывается рассеянное поле. Эти две области разделяет виртуальная граница, которая служит для генерации плоской волны в область полного поля. Разностные уравнения, используемые для расчета компонент поля в прилегающих к этой границе сеточных узлах, отличаются от исходных наличием дополнительных слагаемых, в которых учитывается значение поля падающей волны.
 +
 +<html><center></html>
 +{{:fdtd:fdtd_tfsf.png?250}}\\
 +//Схема численного эксперимента FDTD с использованием метода полного и рассеянного поля.//\\
 +<html></center></html>
 +
 +===== Поглощающие граничные условия ===== 
 +
 +Для устранения нефизичного переотражения электромагнитной волны от границы вычислительного объема и моделирования таким образом ухода волны на бесконечность в FDTD должны использоваться особые поглощающие граничные условия. В настоящее время наиболее успешной реализацией этих условий является помещение вдоль границы вычислительного объема тонкого слоя специального материала, называемого идеально согласованным слоем (Perfectly Matched Layer - PML). Этот материал в идеале полностью поглощает все падающие на него волны без какого-либо отражения независимо от угла падения и длины волны.
 +
 +Понятие идеально согласованного слоя (PML) было введено Дж.-П. Беренгером в 1994 году. Работа такого слоя основывалась на разбиении исходных полей <tex>\vec E</tex> и <tex>\vec H</tex> на две компоненты, для каждой из которых должны решаться свои уравнения. Впоследствии были предложены усовершенствованные формулировки PML эквивалентные первоначальной формулировке Беренгера. Так, в одноосном PML (Uniaxial PML) используется анизотропный поглощающий материал, что позволяет не вводить дополнительные переменные и остаться в рамках исходных уравнений Максвелла. Однако одноосный PML, как и PML в формулировке Беренгера, не удобны тем, что в них отсутствует поглощение затухающих волн, что не позволяет помещать PML близко к рассеивающим телам. Этого недостатка лишен оборотный PML (Convolutional PML), основанный на аналитическом продолжении уравнений Максвелла в комплексную плоскость таким образом, что их решение экспоненциально затухает. CPML также удобнее в ограничении бесконечных проводящих и дисперсных сред. Помимо этого математическая формулировка CPML обладает большей наглядностью и доступностью для понимания.
 +
 +В некоторых случаях использование PML приводит к расходимости FDTD.
 +Эту проблему можно устранить путем использования Additional back absorbing layers technique [[http://www.sciencedirect.com/science/article/pii/S0010465510001839|http]]{{:deinega_-_long-time_behavior_of_pml_absorbing_boundaries_for_layered_periodic_structures.pdf|PDF}}
 +
 +===== Периодические граничные условия ===== 
 +
 +FDTD может применяться не только для конечных, но и для бесконечных периодических структур. Для их моделирования используются периодические граничные условия по одному или по нескольким направлениям. В частности, фотонно-кристаллические пластинки и антиотражающие покрытия являются планарными периодическими структурами: они периодичны по двум направлениям и имеют ограниченную протяженность по оставшемуся направлению. Реальные экспериментальные образцы, конечно, имеют ограниченный размер и по первым двум направлениям, но поскольку этот размер существенно больше их протяженности по непериодическому направлению, влиянием конечности образца с хорошей степенью точности можно пренебречь.
 +
 +В практических приложениях, как правило, наиболее интересны спектры прохождения и отражения от рассматриваемых структур. Для их получения в случае нормального падения в FDTD используется следующий численный эксперимент. Граница между областями полного и рассеянного поля генерирует плоскую волну в форме ограниченного по времени импульса, который падает на структуру (на рисунке ниже этот импульс генерируется на границе 1 и движется слева направо). Часть падающего импульса отражается от структуры, пересекает плоскость, в которой находятся детекторы, меряющие отраженную волну (плоскость 2 на рисунке), и поглощается задним слоем PML. Другая часть проходит сквозь структуру, пересекает плоскость, в которой расположены детекторы, меряющие прошедшую волну (плоскость 2' на рисунке), и поглощается передним слоем PML. Численный эксперимент продолжается до тех пор, пока сигнал не выйдет из вычислительного объема. Время выхода сигнала зависит от конкретной геометрии эксперимента. При увеличении продольной протяженности структуры это время увеличивается по причине увеличения расстояния, которое должен пройти сигнал, чтобы из нее выйти. Время выхода сигнала зависит также и от оптических свойств исследуемой структуры, например, при наличии поглощения оно меньше.
 +
 +<html><center></html>
 +{{:fdtd:fdtd_period.png?300}}\\
 +//Схема численного эксперимента FDTD для расчета спектров прохождения и отражения планарной периодической структуры в случае нормально падения. 1 -- генерирующая электромагнитный импульс граница между областями полного (справа от 1) и рассеянного (слева от 1) поля; 2, 2' - плоскости, в которых расположены детекторы для измерения отраженной и прошедшей волны соответственно.//\\
 +<html></center></html>
 +
 +По окончании эксперимента записанные детекторами значения полей <tex>\vec E(t)</tex>, <tex>\vec H(t)$</tex> переводятся в частотное представление с помощью дискретного преобразования Фурье. Далее, рассчитывая среднее значение вектора Пойнтинга по времени <tex>\vec S(\omega)</tex> на каждом из детекторов по формуле
 +
 +<tex>
 +\vec S(\omega ) = \frac{1}{2} {\rm Re} \left((\vec E(\omega ) \times \vec H(\omega )^* \right),
 +</tex>
 +
 +мы можем посчитать средний поток энергии <tex>W(\omega)</tex>, проходящий через поверхность <tex>A</tex>, на которой расположены детекторы:
 +
 +<tex>
 +W(\omega ) = \int\limits_A {\vec S(\omega)\vec {dA}}.
 +</tex>
 +
 +Нормируя поток энергии прошедшей и отраженной волны на падающий импульс, мы получаем интересующие нас спектры прохождения <tex>T(\omega)</tex> и отражения <tex>R(\omega)</tex>. Поглощение вычисляется как <tex>A(\omega)=1-T(\omega)-R(\omega)</tex>.
 +
 +Схема FDTD для моделирования наклонного падения плоской волны на периодическую структуру может быть найдена здесь [[http://www.opticsinfobase.org/abstract.cfm?uri=ol-33-13-1491|http]]{{:valuev_-_iterative_technique_for_analysis_of_periodic_structures_at_oblique_incidence_in_the_finite-difference_time-domain_method.pdf|PDF}}
 +
 +===== Метод подсеточного сглаживания ===== 
 +
 +Как и в любом другом разностном методе, в FDTD существует проблема неточного отображения границы тела на вычислительную сетку. В непосредственной близости от границы двух сред уравнения Максвелла должны решаться с учетом граничных условий для векторов <tex>\vec E</tex> и <tex>\vec H</tex>. Любая кривая поверхность, разделяющая соседние среды и геометрически не согласованная с сеткой, будет искажаться эффектом "лестничного приближения". В результате порядок точности FDTD понижается с изначального второго до первого.
 +
 +Для решения данной проблемы были предложены различные методы. Первая группа методов основана на изменении способа дискретизации уравнений Максвелла с той целью, чтобы сетка максимально соответствовала рассматриваемой геометрии. Например один из таких методов заключается в увеличении разрешения сетки в тех областях пространства, где расположены тела со сложной геометрической структурой. Также можно видоизменять разностные уравнения в узлах сетки, находящихся вблизи границы между соседними телами. Более радикальным шагом является введение нерегулярной неортогональной сетки, которая полностью согласовывалась бы с рассматриваемой геометрией. Такая сетка автоматически генерируется программой, которая следит за расстановкой всех тел в пространстве. Общим недостатком упомянутых методов является увеличение объема требуемой памяти и порой существенное снижение производительности.
 +
 +Другая группа методов основывается на введении эффективной диэлектрической проницаемости <tex>\varepsilon</tex> вблизи границы между телами. Рассмотрим некоторый контрольный объем <tex>{\delta}x\times{\delta}y\times{\delta}z</tex>, окружающий выбранный узел сетки, и предположим, что он содержит границу раздела между двумя средами со значениями диэлектрической проницаемости <tex>\varepsilon_1</tex> и <tex>\varepsilon_2</tex>. Тогда граничные условия для векторов <tex>\vec{E}$ и $\vec{D}</tex> могут быть реализованы на сетке путем ввода тензора обратной диэлектрической проницаемости в форме:
 +
 +<tex>
 +\hat{\bf \varepsilon}^{-1}={\bf P}<\varepsilon^{-1}>+({\bf 1}-{\bf P})<\varepsilon>^{-1},
 +</tex>
 +
 +где <tex>{\bf P}</tex> - это матрица <tex>P_{ij}=n_in_j</tex> соответствующая проектору на вектор нормали <tex>\vec{n}</tex> к границе между двумя средами, а <tex><></tex> есть усреднение по объему.
 +
 +Более подробно о реализации этого метода в FDTD можно почитать здесь
 +[[http://www.opticsinfobase.org/abstract.cfm?URI=ol-32-23-3429|http]]{{:deinega_-_subpixel_smoothing_for_conductive_and_dispersive_media.pdf|PDF}}
 +
 +====== Применение FDTD ======
 +
 +FDTD применяется в самых разнообразных приложениях: для расчета рассеяния падающей волны на различных структурах, для моделирования излучения источника в заданном электромагнитном окружении, для изучения оптических свойcтв резонаторов и волноводов. Ниже мы изложим, какие величины при этом измеряются.
 +
 +===== Предварительные замечания =====
 +
 +Напомним, что любое решение уравнений Максвелла <tex>{\bf F}({\bf r},t)</tex> (<tex>\bf F</tex> это либо <tex>\bf E</tex>, либо <tex>\bf H</tex>) в случае линейных сред и при отсутствии свободных зарядов и источников тока может быть представлено в виде суперпозиции гармонических полей:
 +
 +<tex>
 +{\bf F}({\bf r},t)={\bf A}({\bf r},\omega)\cos{\omega t}+{\bf B}({\bf r},\omega)\sin{\omega t}
 +</tex>.
 +
 +Удобно рассматривать вектор <tex>{\bf F}</tex> как действительную часть от комплексного вектора 
 +<tex>
 +{\bf F_c}={\bf C}\exp(-i\omega t)
 +</tex>, 
 +где
 +<tex>
 +{\bf C}={\bf A}+i{\bf B}
 +</tex>:
 +
 +<tex>
 +{\bf F}={\rm Re}({\bf F_c})
 +</tex>.
 +
 +Комплексная зависимость <tex>\exp(-i\omega t)</tex> не несет никакого физического смысла и вводится лишь для удобства.
 +
 +Вектор плотности потока электромагнитной энергии (вектор Пойнтинга) равен <tex>{\bf P} = {\bf E} \times {\bf H}</tex>. В монохроматической волне его величина быстро осциллирует, чего большинство приборов отследить не в силах. Вместо этого измеряют среднее вектора Пойнтинга по времени <tex>{\bf S}</tex>:
 +
 +<tex>
 +{\bf S}=<{\bf P}>=\frac{1}{\tau}\int_t^{t+T}{\bf P}(\tau) d\tau
 +</tex>,
 +
 +где временной интервал <tex>T</tex> существенно больше частоты <tex>1/\omega</tex>.
 +
 +Можно показать, что для монохроматической волны
 +
 +<tex>
 +{\bf S}=\frac{1}{2}{\rm Re}({\bf E_c}\times{\bf H^*_c})
 +</tex>.
 +
 +Абсолютная величина вектора <tex>{\bf S}</tex> называется интенсивностью света <tex>I</tex>.
 +
 +===== Рассеяние на теле =====
 +
 +Рассмотрим задачу о рассеянии падающей волны на заданном теле. Поле вне тела можно представить в виде суммы падающей и отраженной волн:
 +
 +<tex>
 +{\bf E}={\bf E_i}+{\bf E_{s}}, \qquad
 +{\bf H}={\bf H_i}+{\bf H_{s}}
 +</tex>.
 +
 +{{fdtd:flux.png?300}}
 +
 +В задаче о рассеянии интересуются количеством энергии рассеянной волны, ее угловой зависимостью, а также количеством поглощаемой телом энергии.
 +
 +Окружим тело воображаемой замкнутой поверхностью <tex>A</tex>. Средняя скорость, с которой энергия притекает в ограниченный ею объем, равна
 +
 +<tex>
 +W_{abs}=-\int{\bf S}\cdot{\bf n} dA
 +</tex>,
 +
 +где <tex>{\bf n}</tex> - нормаль к поверхности.
 +
 +Если <tex>W_{abs}>0</tex>, энергия поглощается внутри объема. Если внешняя среда является непоглощающей, то <tex>W_{abs}</tex> есть энергия, поглощаемая телом.
 +
 +В нашем случае среднее по времени от вектора Пойнтинга <tex>{\bf S}</tex> может быть записано в виде суммы трех слагаемых (нижний индекс <tex>c</tex> для комплексных векторов <tex>{\bf E}</tex> и <tex>{\bf H}</tex> мы опускаем)
 +
 +<tex>
 +{\bf S}=\frac{1}{2}{\rm Re}({\bf E}\times{\bf H^*})={\bf S_i}+{\bf S_{sca}}+{\bf S_{ext}}
 +</tex>,
 +
 +где каждое слагаемое по отдельности равно
 +
 +<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>.
 +
 +Наличие <tex>{\bf S_{ext}}</tex> есть следствие интерференции падающей и отраженной волн.
 +
 +Учитывая только что сказанное, для <tex>W_{abs}</tex> получаем
 +
 +<tex>
 +W_{abs}=W_i-W_{sca}+W_{ext}
 +</tex>,
 +
 +<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>.
 +
 +В случае непоглощающей внешней среды <tex>W_i=0</tex>, и
 +
 +<tex>
 +W_{ext}=W_{abs}+W_{sca}.
 +</tex>
 +
 +Поскольку рассматриваемые величины линейно зависит от интенсивности падающей волны <tex>I</tex>, имеет смысл рассматривать отношения
 +
 +<tex>
 +C_{ext}=\frac{W_{ext}}{I}, \qquad
 +C_{abs}=\frac{W_{abs}}{I}, \qquad
 +C_{sca}=\frac{W_{sca}}{I}
 +</tex>,
 +
 +имеющее размерность площади, и называемые сечениями экстинции, поглощения и рассеяния
 +
 +Можно определить эффективности экстинции, рассеяния и поглощения
 +
 +<tex>
 +Q_{ext}=\frac{C_{ext}}{G}, \qquad
 +Q_{abs}=\frac{C_{abs}}{G}, \qquad
 +Q_{sca}=\frac{C_{sca}}{G}
 +</tex>,
 +
 +где <tex>G</tex> - площадь проекции тела на плоскость перпендикулярную направлению падающей волны (для сферы радиуса <tex>r</tex> имеем <tex>G=\pi r^2</tex>). На длинах волн больших размеров тела, тело может рассеивать и поглощать больше света, чем падает на его геометрическое сечение. Соответствующие эффективности при этом будут больше единцы.
 +
 +===== Матрица амплитуд рассеяния =====
 +
 +Помимо энергии рассеянной волны часто интересуются ее угловым распределением, которое характеризуется матрицей амплитуд рассеяния. Вводится она следующим образом.
 +
 +Пусть на тело падает монохроматическая линейно поляризованная волна. 
 +
 +{{fdtd:far.png?300}}
 +
 +Направим ось <tex>z</tex> вдоль направления падения. Выберем какую-либо точку внутри тела в качестве начала координат <tex>O</tex>, а оси <tex>x</tex>, <tex>y</tex> направим произвольно. Нашей системе координат соответствуют ортонормированный базис векторов <tex>{\bf e_x}</tex>, <tex>{\bf e_y}</tex>, <tex>{\bf e_z}</tex>. 
 +
 +Направление рассеяния <tex>{\bf e_r}</tex> вместе с направлением падения <tex>{\bf e_z}</tex> определяют плоскость рассеяния, которая однозначно задается азимутальным углом <tex>\phi</tex>, если только <tex>{\bf e_r}</tex> не параллельно оси <tex>z</tex>. В этих двух случаях (<tex>{\bf e_r}=\pm{\bf e_z}</tex>) любая плоскость параллельная оси <tex>z</tex> может считаться плоскостью рассеяния. 
 +
 +Вектор <tex>{\bf E_i}</tex> падающего электрического поля, лежащий в плоскости <tex>xy</tex>, удобно разложить на две компоненты, одна из которых параллельна плоскости рассеяния, а другая перпендикулярна ей
 +
 +<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>
 +
 +Вектора ортонормированного базиса 
 +
 +<tex>
 +{\bf e_{\perp{}i}}=\sin\phi{\bf e_x}-\cos\phi{\bf e_y}
 +</tex>,
 +
 +<tex>
 +{\bf e_{\parallel{}i}}=\cos\phi{\bf e_x}+\sin\phi{\bf e_y}
 +</tex> 
 +
 +образуют правую тройку с вектором <tex>{\bf e_z}</tex>:
 +
 +<tex>
 +{\bf e_{\perp{}i}}\times{\bf e_{\parallel{}i}}={\bf e_z}
 +</tex>
 +
 +Кроме того, имеем
 +
 +<tex>
 +{\bf e_{\perp{}i}}=-{\bf e_{\phi}}, \qquad {\bf e_{\parallel{}i}}=\sin\theta{\bf e_r}+\cos\theta{\bf e_{\theta}}
 +</tex>,
 +
 +где <tex>{\bf e_r}</tex>, <tex>{\bf e_{\theta}}</tex>, <tex>{\bf e_{\phi}}</tex> образуют базис, связанный с полярной системой координат <tex>(r, \theta, \phi)</tex>.
 +
 +Если компоненты падающего поля в направлениях x и y обозначить через <tex>E_{xi}</tex> и <tex>E_{yi}</tex>, то
 +
 +<tex>
 +E_{\parallel{}i}=\cos\phi{}E_{xi}+\sin\phi{}E_{yi}
 +</tex>,
 +
 +<tex>
 +E_{\perp{}i}=\sin\phi{}E_{xi}-\cos\phi{}E_{yi}
 +</tex>.
 +
 +Известно, что на больших расстояниях от тела (<tex>{\bf kr}\gg{}1</tex>) (в так называемой дальней зоне) рассеянное поле приближенно является поперечным (<tex>{\bf e_r} \cdot {\bf E_s} \sim 0 </tex>), и, может быть записано как
 +
 +<tex>
 +{\bf E_s}=E_{\parallel{}s}{\bf e_{\parallel{}s}}+E_{\perp{}s}{\bf e_{\perp{}s}}
 +</tex>,
 +
 +где 
 +
 +<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>.
 +
 +Базисный вектор <tex>{\bf e_{\parallel{}s}}</tex> параллелен, а <tex>{\bf e_{\perp{}s}}</tex> перпендикулярен плоскости рассеяния. Заметим, однако, что <tex>{\bf E_s}</tex> и <tex>{\bf E_i}</tex> определены в разных системах базисных векторов. Тем не менее, вследствие линейности уравнений Максвелла, связь между ними можно представить в матричном виде
 +
 +{{fdtd:sc_matrix.png?250}},
 +
 +где комплексные числа <tex>S_j</tex> (<tex>j=1,2,3,4</tex>) составляют амплитудную матрицу рассеяния, зависящую от угла рассеяния <tex>\theta</tex> и азимута <tex>\phi</tex>.
 +
 +Экспериментальное измерение матричных элементов <tex>S_j</tex> затруднительно, поэтому обычно меряют элементы матрицы рассеяния, связывающей так называемые параметры Стокса для падающей и рассеянной волн. Это действительная матрица <tex>4\times{}4</tex>, из 16 элементов которой независимы только семь, которые в свою очередь могут быть выписаны через абсолютные значения <tex>\left|S_j\right| (j=1,2,3,4)</tex> и разницы в фазе между <tex>S_j</tex>. Желающие узнать об этом более подробно, могут ознакомиться с главой 3.3 книги ((C. F. Bohren and D. R. Huffman: Absorption and Scattering of Light by Small Particles, Wiley-Interscience, New York (1983) )).
 +
 +===== Коэффициенты прохождения и отражения для плоскопериодических структур =====
 +
 +===== Радиационное затухание источника =====
 +
 +====== Приложения ======
 +
 +Программа EMTL использовалась для решения различных оптических задач на многих кластерных вычислительных системах.
 +Ниже представлен лист избранных публикаций:
 +
 +  * A. Deinega, I. Valuev, B. Potapkin and Yu. Lozovik, "Minimizing light reflection from dielectric textured surfaces," JOSA A 28, 770 (2011) [[http://www.opticsinfobase.org/abstract.cfm?URI=josaa-28-5-770|http]]{{:deinega_-_minimizing_light_reflection_from_dielectric_textured_surfaces.pdf|PDF}}
 +  * S. Zalyubovskiy et. al., "Theoretical limit of localized surface plasmon resonance sensitivity to local refractive index change and its comparison to conventional surface plasmon resonance sensor", JOSA A 29, 994 (2012) [[http://www.opticsinfobase.org/josaa/abstract.cfm?uri=josaa-29-6-994|http]]{{:zalyubovskiy_-_theoretical_limit_of_localized_surface_plasmon_resonance_sensitivity.pdf|PDF}}
 +  * A. Deinega, S. John, "Solar power conversion efficiency in modulated silicon nanowire photonic crystals", J. Appl. Phys. 112, 074327 (2012) [[http://jap.aip.org/resource/1/japiau/v112/i7/p074327_s1|http]]{{:deinega_-_solar_power_conversion_efficiency_in_modulated_silicon_nanowire_photonic_crystals.pdf|PDF}}
 +  * S. Belousov et. al., "Using metallic photonic crystals as visible light sources", Phys. Rev. B 86, 174201 (2012) [[http://prb.aps.org/abstract/PRB/v86/i17/e174201|http]]{{:belousov_-_using_photonic_crystals_as_visible_light_sources.pdf|PDF}}
 +  * A. Deinega, S. Eyderman, S. John, "Coupled optical and electrical modeling of solar cell based on conical pore silicon photonic crystals", J. Appl. Phys. 113, 224501 (2013) [[http://jap.aip.org/resource/1/japiau/v113/i22/p224501_s1|http]]{{:deinega_-_coupled_optical_and_electrical_modeling_of_solar_cell_based_on_conical_pore_silicon_photonic_crystals.pdf|PDF}}
 +Смотрите также {{:deinega_thesis.pdf|диссертацию Алексея Дейнеги}}.
 
/var/www/fdtd.kintechlab.com/docs/data/pages/ru/fdtd.txt · Последние изменения: 2013/08/08 21:27 — deinega     Наверх