Разработка комплекса программ решения электродинамических задач с использованием массивно-параллельных вычислительных систем тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Семенов, Алексей Николаевич

• Семенов, Алексей Николаевич
• кандидат наук
• 2013

Оглавление диссертации кандидат наук Семенов, Алексей Николаевич

Содержание

Введение

Глава 1. Создание комплекса программ численного моделирования уравнений Максвелла

1.1. Класс решаемых электромагнитных задач

1.2. Модели диэлектрической среды

1.3. Граиичные условия

1.4. Источники электромагнитных волн

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

1.6. Использование поглощающего слоя при моделировании открытых задач

1.7. Метод нолного/рассеяного поля

1.8. Численная реализация модели Друде для метаматериалов в FDTD схеме

1.9. Параллельный алгоритм численного решения уравнений Максвелла FDTD-методом

1.10. Гибридный подход к программной реализации параллельного алгоритма

1.11. 1D декомпозиция для трехмерных задач с протяженной областью

1.12. 3D декомпозиция для произвольных трехмерных задач

1.13. Особенности целевой архитектуры вычислительной системы

1.14. Эффективность распараллеливания разработанных алгоритмов численного решения

Глава 2. Применение комплекса программ для решения оптических задач с использованием метаматериалов

2.1. Апробация в базовых оптических задачах

2.2. Моделирование открытого микрорезонатора со слоем метама-териала

2.3. Моделирование идеальной линзы

Глава 3. Применение комплекса программ для решения задачи оценки коэффициента прохождения волноводной моды

в прямоугольном волноводе

3.1. Моделирование электромагнитного волновода с диэлектрической вставкой

3.2. Получение коэффициента прохождения из полноволнового решения

3.3. Расчёт коэффициента прохождения для однородной тонкой вставки

3.4. Исследование зависимости коэффициента прохождения от физических параметров неоднородности внутри вставки

Заключение

Литература

Введение диссертации (часть автореферата) на тему «Разработка комплекса программ решения электродинамических задач с использованием массивно-параллельных вычислительных систем»

Введение

Актуальность работы. При моделировании распространения электромагнитных волн в радиотехнике [22],[23],[8], микроэлектронике [28],[27], задачах нанооптики [29],[30], [10] и биодиагностики [25], [26], требуется решение полной трёхмерной системы уравнений Максвелла на сетках больших размеров. При реализации на вычислительной системе размер таких сеток составляет порядка терабайта и более. Большие сетки необходимы для описания исоднородностей, характерный размер которых значительно меньше длины волны, а общий размер исследуемой системы равен десяткам длин волн. Для решения таких задач необходимы вычислительные системы с большим объёмом памяти и высокой производительностью: суперкомпьютеры с массивно-параллельной архитектурой. Для их использования требуются специализированные комплексы программ, способные эффективно использовать терабайтные объёмы оперативной памяти и более чем тысячи процессоров вычислительной системы. При реализации комплекса программ необходимо также учитывать специфику архитектуры конкретной вычислительной системы. Комплекс должен быть многофункциональным и рассчитанным на решение задач как в ограниченной, так и в не ограниченной области.

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

Одним из наиболее распространённых способов численного решения пол-

ной системы уравнений Максвелла является явный метод конечных разностей во временной области, в иностранной литературе fmitc-difFcrencc timedomain (FDTD) [1]. Метод основан на использовании сетки Йи [4], в узлах которой располагаются компоненты электромагнитного поля. Основываясь на уравнениях Максвелла в интегральной форме, можно построить их конечно-разностную аппроксимацию. Поскольку FDTD-мстод решает задачу во временной области, он позволяет решать импульсные и негармонические задачи и получить результат для широкого спектра длин волн за один расчёт. Это необходимо при решении задач, в которых неизвестны резонансные частоты, или в случае моделирования широкополосных сигналов. Также FDTD-метод позволяет проследить временную эволюцию распространения волны в моделируемом объёме, что может быть необходимо для детального исследования процесса формирования устойчивой картины распределения электромагнитных волн. Метод позволяет непосредственно моделировать краевые эффекты и эффекты экранирования и удобен при задании анизотропных, дисперсных и нелинейных сред.

При численном решении FDTD-методом задач в неограниченном пространстве на границах расчётной области чаще всего используется искусственный поглощающий слой PML (Perfectly matched layer) [11], [9]. Рассеянная волна, уходя на бесконечность, поглощается анизотропным слоем на границе и, при правильном выборе параметров поглощающего слоя и разностной сетки, отражение от PML слоя минимально [5],[3].

Величина шага дискретизации по пространству должна быть значительно меньше исследуемых длин волн и размера исследуемой структуры. Это может потребовать сеток с маленьким шагом по пространству, что означает большие затраты памяти и возросшее время расчёта. При моделировании дисперсных и анизатропных сред также растут необходимые вычислительные ресурсы. Актуальной становится задача создания программного комплекса

на высокопроизводительных массивно-параллельных вычислительных системах с большим объёмом оперативной памяти, который позволит эффективно решать класс электродинамичских задач, описываемых полной системой уравнений Максвелла.

Существует ряд комплексов, реализующих параллельный FDTD-мстод на кластерных вычислительных системах([3б], [37]). Однако, использование таких комплексов, как правило, рассчитано на массивно-параллельные вычислительные системы общей архитектуры, что делает невозможным эффективно использовать их на специализированных суперкомпьютерах и полностью задействовать их вычислительные мощности. На факультете ВМК Московского государственного университета имени М.В. Ломоносова установлен суперкомпьютер IBM BlueGene/P, имеющий 8192 процессоров и 4 терабайта оперативной памяти. Разработанный в диссертации программный комплекс был эффективно реализован для этого суперкомпьютера и вычислительных систем схожей архитектуры.

Цель диссертационной работы. Разработка параллельных алгоритмов и комплекса программ для численного решения полной системы трёхмерных уравнений Максвелла на современных многопроцессорных вычислительных системах и суперкомпьютерах для моделирования электродинамических задач радиотехники, микроэлектроники, нанооптики и биодиагностики, требующих для их описания и расчёта терабайтных данных.

Научная новизна. Разработан программный алгоритм для решения крупномасштабных электродинамических задач, который с помощью гибридных технологий MPI и ОрепМР может эффективно выполнятся на современных массивно-параллельиых вычислительных системах. На основе предложенного алгоритма создан программный комплекс на многоядерных вычислительных системах и суперкомпьютерах серии IBM Blue Gene/P. Показано, что комплекс имеет хорошую масштабируемость и может эффективно

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

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

На защиту выносятся следующие основные результаты и положения:

• Разработан программный алгоритм для решения электродинамических задач, использующий гибридные технологии MPI и ОрепМР. Создан многофункциональный комплекс программ на многоядерных вычислительных системах и суперкомпьютере IBM BlueGene/P, позволяющий решать широкий круг электродинамических задач. Показано, что созданный комплекс имеет хорошую масштабируемость и может эффективно выполняться на массивно-параллельных системах, оперируя с терабайтами вычислительных данных.

• С помощью созданного комплекса, исследована зависимость субволновой разрешающей способности идеальной линзы из метаматериала от её физических параметров. Найдены значения параметров, при которых возможно получить субволновую разрешающую способность.

• С помощью комплекса впервые решена задача вычисления добротности

открытого резонатора Фабри-Перо со вставкой из метаматериала, ранее проводившегося только для не субволновых размеров резонатора.

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

Апробация работы. Основные результаты диссертации докладывались на:

• Международной конференции "Progress in Electromagnetics Research Symposium"(Stockholm, Sweden, 2013 год)

• Международной конференции "Progress in Electromagnetics Research Symposium" (Moscow, Russia, 2012 год)

• Международном научном семинаре "Inverse Problems and Large-Scale Modeling with Applications in Electromagnetics" (Karlstad, Sweden, 2012 год)

• VI Всероссийской конференция "Радиолокация и радиосвязь" Российская академия наук Институт радиотехники и электроники им. В. А. Котелышкова (Москва, 2012 год)

• Объединенном Фельдовском семинаре по электродинамике и антеннам. (Москва, 2012 год) Российская академия наук Институт радиотехники и электроники им. В. А. Котелышкова

• XIII Всероссийской школы-семинара "Физика и применение микроволн" (Звенигород, МГУ им. М.В. Ломоносова, 2011 год)

• Международной суперкомпьютерной конференции "Научный сервис в сети Интернет: экзафлопеное будущее" (Новороссийск, 2011 год)

• Научной конференции "Ломоносовские чтения"(Москва, МГУ им. М.В. Ломоносова, 2011 год)

Публикации. Материалы диссертации опубликованы в 9 печатных работах, из них б статей в журналах списка ВАК ([40]-[45]) и 3 тезиса докладов ([46]-[48]).

Структура и объем диссертации Диссертация состоит из введения, 3 глав, заключения и библиографии. Общий объем диссертации 100 страниц, из них 89 страниц текста, включая 23 рисунка. Библиография включает 60 наименований на 7 страницах.

Глава 1

Создание комплекса программ численного моделирования уравнений Максвелла

В диссертационной работе автором был создан комплекс программ ЕМЛУБокегЗБ, предназначенный для численного решения полной системы уравнений Максвелла.

1.1. Класс решаемых электромагнитных задач

При разработке комплекса программ ЕМ\¥8о1уегЗО исходной постановкой задачи было найти решение системы уравнений Максвелла в оптической задаче в отсутствии внешних зарядов:

гоШ

дВ

Зм — 0яВ

ал)

дг

гоШ = —- + 3 - (7ЕВ дЬ

(НуВ = О (НуВ = О В = цН Б = еЕ

где £ = е(г, а;), /2 = /г(г,а;) — диэлектрическая и магнитная проницаемости, зависящие от частоты с электрическими потерями, а J и Ям ~ плотность электрического тока и плотность тока магнитного заряда соответственно. Так же будем считать что, материалы могут иметь магнитные и электрические потери энергии, определяемые не зависящими от частоты значениями электрической и магнитной проводимостей ое и Ое-

Граничные условия Начальные условия Модели сред

Условия первого рода. Любая компонента электромагнитного поля равняется заданной функции Ё{х,у,г,1)-Дх,у,2Л геС Й(х,у,2,1)-гдеС-гратцаобластц Задание начального профиля электрического и магнитного ПОЛЯ. Е(х,у, /,Г = 0) = /(х,>%2); Н(х;у,гЛ = 0) = §(х:у}г); е=е(г), М=Р(Г)

£ = = Ц(г,ь») Дисперсная модель Друде. Позволяет моделировать метаматериалы.

Периодические условия. Метод полного/рассеянного поля. Решение задачи может искаться в виде суммы Е полное =Епадэющее+Ерассеянное, что позволяет моделировать падение из бесконечности плоскопараллельной волны под произвольным заданным углом.

Поглощающий слой иРМ1. поглощающий слой, позволяющий решать задачи е неограниченной области.

Таблица 1. Функциональность комплекса программ ЕМ\Уво1иегЗВ

Комплекс позволяет рассматривать постановки задач, где граничные условия и параметры среды могут иметь вид, представленный в Таблице 1. Реализуемая функциональность комплекса программ позволяет моделировать широкий круг задач, таких как рассеивание света на нанообъектах в открытой области, распространение электромагнитных волн в волноводе с диэлектрическими вставками, рассеивание света па фотонных кристаллах и задач, связанных с материалами, имеющими отрицательный коэффициент преломления. Например, задачу моделирования рассеивания волн в открытых микрорезонаторах на основе слоистых метаматериалов [40], а также задачу прохождения электромагнитной волны через идеальную линзу.

1.2. Модели диэлектрической среды

При описании свойств материалов с отрицательным коэффициентом преломления нужно иметь в виду, что эти материалы должны обладать частотной дисперсией. Действительно, если е и /х оба отрицательны, то при отсутствии дисперсии плотность энергии электрического и магнитного поля,

равная [15]:

IV =1-{еЕ2 + /х#2), (1.2)

будет отрицательной. При наличии частотной дисперсии плотность энергии (1.2) имеет вид [15]:

дЩи1\Е2 + дЫимн2

2 дш ди

В ряде работ [13] [18] предложено описывать метаматериалы дисперсной моделью Друде [13]. Для модели Друде е{и>) и ¡1{и}) принимают вид:

е(и) = £0 I 1 - ре

О)2 — гШ'Уе

и? \ {1А)

= /¿о ( 1--—- '

•'рт

,2 ^ГТГ

- гш'Ут у

где МрсШрт - электрическая и магнитная плазменные частоты, а 7е,7т - частоты столкновений частиц, характеризующие электромагнитные потери в материале. Заметим, что эта модель удовлетворяет условиям

аКыМ>0а|р(цМ>0- (15)

дсо ди

1.3. Граничные условия

При решении задач с наличием границ используются заданные граничные условия. Например, для описания идеального проводника (РЕС)[1], используется условия равенства нулю тангенциальных компонент электромагнитного поля. Также, в рамках комплекса реализованы периодические граничные условия.

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

используются поглощающие граничные условия. Существует ряд аналитических условий на бесконечности для уравнений Максвелла [31], [32]. Большинство из них имеет узкую область применения и ряд существенных ограничений [34] [33]. При численном решении задач на бесконечности также существуют различные методы вычисления значений поля на границе, зачастую привносящие значительные искажения в полученное решение [1]. Хорошие результаты при численном моделировании различных задач электродинамики показали поглощающие граничные условия Perfect matched layer(PML)[1]. PML требует аккуратного выбора параметров поглощающего слоя и разностной сетки для получения приемлемого результата.

Есть несколько эквивалентных формулировок PML. Первоначальная формулировка J. P. Berenger [2] в зарубежной литературе именуется split-field PML. В ней решение искусственным образом разбивается на сумму двух опять же искусственных компонент электромагнитного поля. На данный момент, более общей формулировкой считается uniaxial PML или UPML [5], в котором в области PML решается тоже волновое уравнение, но в комбинации с искусственной анизатропной поглощающей средой[3]. Обе эти формулировки были выведены напрямую из решения задачи, для электромагнитной волны, падающей на границу поглощающей среды под произвольным углом, и с последующим введением дополнительных условий отсутсвия отражения падающей волны от границы раздела сред. Так же UPML поглощает экспоненциально затухающие (эвапесцентные) волны, чего не наблюдается в split-field формулировке, и что является важным свойством при моделировании идеальной линзы, поскольку именно эвапесцентные волны обеспечивают субволновую разрешающую способность линзы. Тем не менее, UPML подход является достаточно трудоемким при обобщении его на другие уравнения и другие координатные системы. Также из рассмотрения теряется тот важный факт, что PML работает (не производит отражений) для неоднородных сред,

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

1.4. Источники электромагнитных волн

В современных задачах электродинамики источники электромагнитных волн могут быть устроены крайне сложно, что затрудняет их моделирование, если учитывать каждую деталь такого источника. Необходимый шаг дискретизации может быть настолько мал, что для моделирования потребуется не только очень большой объем оперативной памяти, но и очень долгое время работы самой программы. В таких случаях необходимо либо упрощение рассматриваемой модели, либо использование массивно-параллельных вычислительных систем. Корректное задание источника является важным вопросом при моделировании с помощью РБТБ-метода. Существует множество видов возможных реализаций источников. Один из них - так называемый "жесткий" источник, который создаст в своей окрестности силовые поля, задаваемые некоторой временной функции. Преимуществом такого подхода является то, что падающая волна имеет ту же форму, что и возбуждающий сигнал. Поля в окрестности источника, известные на каждом временном шаге, не учитываются в явной РБТБ схеме при вычислении компоненты поля. Поскольку при таком подходе возникают нежелательные отраженные волны из области источника, такой тип возбуждения для РОТО-мстода в настоящее время используется редко. Другой тип возбуждения, реализуемый с помощью выражения для источника тока и заряда в уравнении Максвелла, широко используется на практике поскольку он не вызывает отраженных волн в окрестности источника. Два типа источников, описанные выше, могут быть отнесены к "локальным" источникам, которые лежат полностью в расчетной

области.

Для задач рассеяния волна падает из бесконечности, и моделировать такую задачу с помощью "локальных" источников затруднительно. Можно выделить два распространенных подхода к решению этой проблемы. В одном из них используется метод рассеянного поля, при котором неизвестной величиной является рассеянное поле, а поле, порождаемое падающей волной, представляется в виде выражения для источника в уравнении, вычисляемом па каждом временном шаге FDTD-метода. Второй метод - это метод полного/рассеянного поля, при котором в качестве неизвестного внутри некоторой подобласти рассматривается полное поле, а вне ее - рассеянное поле. При таком подходе поле падающей волны рассчитывается только на границе двух областей. Стоит отметить, что эти два метода могут применять не только к источникам в виде плоской волны, но и к любому другому типу источников, включающему в себя локальный источник.

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

FDTD-метод встречается в русскоязычных статьях как метод КР-ВО(конечных разностей во временной области). Метод является явной разностной схемой, позволяющей выразить значения компонент электромагнитного поля на следующем временном слое через значения на текущем слое. FDTD-метод основан на использовании ячеек Йи (Yee lattice) [4], составляющих разностную сетку. Обозначим узел пространственной сетки как:

(г, j, к) = (iAx, jAy, kAz) (1.6)

Здесь Ах, А у, Аг шаги однородной сетки по направлениям х, у, г, а г — О ..М, 7 - О..Аг, Л = 0..Р.

Значение сеточной функции и в узле в момент времени £п обо-

значим как:

и(гАх,зАу,кАх,пА1) = ип(ъ,^к) (1.7)

где At шаг по времени. Все компоненты электромагнитного поля Ех, Еу, Ег, Нх, Ну, Н2 берутся в разных точках ячейки, компоненты магнитного поля Н находятся в центрах граней, а компоненты электрического поля Е - посередине ребер. Ячейка Йи изображена на Рис. 1.1(а).

Рис. 1.1а. Расположение компонент поля в ячейке потоковой сетки.

¿Hz(i+1/2,j+1/2,k)

-Hy(i+1/2,j,k-1/2)

1 / / / / т р

« * * ч\\ЁХ \ » » 1 \ » » » \ ч \ \ ч* > л\\\ LXxJ ' У / ' »' / / / . ' / / ' / / / / » t » /

Ну (i+1/2,j,k+1/2)

i

<6)

Рас. 1.16. Расположение контура и поверхности интегрирования и si

для компоненты Ех.

M(U+U) -Ez(i,j,k-1/2)

' / ' . * / * / • / / / 'л * ф / * f / ф t % » WV WV » % \ » ,\\V

\ Ч \ ' ) \\ \ \ / » » % % ' % » » % \ % % % » . % % % » ч \ \ ' ///. / / # » в л * * / / * • / • s / / ' / * // / /

p(i+1,j+1,k)

Ez (i+1,j,k-1/2) p(i+1,j+1,k-1)

i

W

Рис. l.le. Расположение контура и поверхности интегрирования /2 и S2 для компоненты Ну и расположение е и ц.

Будем считать, что диэлектрическая проницаемость непрерывна внутри ячейки Йи, а магнитная проницаемость постоянна в ячейке сдвинутой на —1/2 по х, у и г. Значение е{1,], к) = ё{г + + \,к + \) известно в центре ячеек, а = в узлах. При этом е может изменять-

ся лишь на гранях ячеек сетки Йи, а магнитная проницаемость ^ лишь на гранях сдвинутых ячеек сетки Йи, т.е. на координатных плоскостях, проходящих через центры ячеек сетки. Тогда е может быть разрывна на плоскостях х = Х{,у — = гк, где Х{ = гAx,yj = уАу,Хк — г Ах, а ¡л разрывна на х = Х(+г,у = = Аналогично, электрическая проводимость непре-

рывна внутри ячейки Йи, а магнитная проводимость непрерывна в ячейке сдвинутой на —1/2 по х, у и г. Значение оя(г, 3, к) = ое{ъ 4-3 4-к +известно в центре ячеек, а <тя(г,:/, к) = он{г,3, к) в узлах. Тогда оЕ может быть разрывна на плоскостях х = гсг-, у = уу, г = где Х{ = {Ах, yj = 3Ау, Хк = гАг, а он разрывна на х = х^1,у = у^к, % =

Рассмотрим уравнения Максвелла в интегральной форме:

д_ дь

д

Шз = фШ/-7; -—

дЬ

Вйв = о Ес11 + 1М;

(1.8)

где I = (7 - (тЕ), 1м = Ым ~ °н)- В точке {г + ^З, к) - е и оЕ разрывны но направлениям у иг. Поскольку тангенциальная к границе раздела компонента Ех непрерывна, то из уравнений (1.8), взятых по контуру 1\ и области 51 (см. Рис. 1.1(6)), для Их в соответствии с расположением компонент поля в ячейке Йи (см. Рис. 1.1(а)), перепишем левую часть уравнения (1.8) как:

дг(£+0.5) Ау{з+0.5)

Бх (у, г)(1у(1г ~ { при Бх = е(г)Ех имеем } «

д_ дЬ

дг(а;-0.5) д1/(7-0.5)

9 „ д А

где< £{1,3, к) >= (ё(г+1/2^-1/2,а;-1/2)+£(г+1/2^-1/2^+1/2)+£(ш/2^+1/2,а;-1/2) +

^(г+1/2Л-1/2,А;+1/2))/4. Такое усреднение следует из условий непрерывности для компоненты Ех на поверхности ¿ь В правой части уравнения (1.8), в соответствии с расположением компонент поля в ячейке Йи (см. Рис. 1.1(6)), интеграл по контуру ¿1 примет вид:

()Нdl = (Ну ¿+0.5J,Jfc+0.5 — Ну j+o.5,i,A;-0.5)A^— h

№ г+0.5 j+0.5,fc г+0.5 j—0.5,

(1.10)

Для тока / аналогично (1.9) имеем:

(J - (ТЕ) « J{i,j, к)- < (TEliJik) >, (1.11)

ГДе < CTE(iJ, к) >= {(fE(i+l/2J-l/2,k-l/2)+^E(i+l/2J-l/2Ml/2)^E{i+l/2,j+l/2,k-l/2) +

<7E(m/2j+i/2,fc+i/2))/4. Подставляя (1.9), (1.10), (1.11) в (1.8), получаем: $

< £(i,j,k) > AyAz = (Ну i+0.5,i,A;+0.5 ~ Ну ¿+0.5,j',fc-0.5)Az~

(Hz ¿+0.5J+0.5,fc - Hz j+0.5j-0.5,i;)A7/ + J{l, j, k)~ < <JE(iJtk) > ■

(1.12)

Выражение для других компонент поля Е имеет аналогичный вид.

Для магнитного поля Ну в точке (г + к + магнитная проницае-

¿л ¿л

мость \х разрывна по направлениям хи z, тогда аналогично (1.9) из уравнений (1.8) в соответствии с Рис 1.1(b). для компоненты Ну по контуру ¡2, имеем: Дх(г+1) Az(k+1)

Ву j(x, z)dzdx « { при By = цНу имеем } «

Лх(г) Лг(к)

~< V(i,j,k) > §1НУ j&zAx,

д_

dt

д *

— Byí¿S2 = (Ez ¿+1 0.5 — ^z i,i,fc-0.5)A^— ¿2

(Ex i+0.5,j,k—l — Ex i+0.5J,k)Az,

(1.14)

где < 1Ц1,з,к) >= + ++/^(»-»-1 Анало-

гично < <7^,*) >= (<тя( ~+1,к) + ^я«^^, + ^я^Т,^,*) + ^я(г+1~+1,,+1))/4. Нспре-рывность тангенциальных составляющих Е и Н сохраняется вдоль границы раздела сред. Используемая аппроксимация не требует специальных условий сшивки на границах раздела, поскольку выполнение этих условий следует из уравнений Максвелла в интегральной форме. Разностная дивергенция роторов полей Е и Н сохраняется точно.

Перепишем полученные вихревые уравнения в покомпонентную форму для декартовой системы координат. Получаем следующую систему из 6 скалярных уравнений, эквивалентных уравнениям, полученны[ из балансных соотношений (1.8):

дНх 1 гдЕу dEz .

= ^—d я --Jm*+ <°н> Нх) ,

dt < ¡i> az ду

дНц 1 rdEz дЕт . _ тт

«* = ^- ~df - ^+ < > Н»)]■ дН,_ 1 ОЕ, _ ОЕ, _ < ан >

dt < ¡i> ду дх

дЕх 1 rdHz дНу

1.15)

1.16)

1.17)

1.18) 1.19)

я, ^ ^я я V ■ — ' — ч1-20)

dt < е > дх ду

Эта система из 6 уравнений в частных производных представляет из себя основу для численного FDTD алгоритма, который позволяет описать общее

RP - ^ - V+ < °Е > ЕХ)Ъ

dt < £ > ду dz dEv 1 гдНх dHz , т

eEí_ i ав>ЕЛ

взаимодействие электромагнитных волн с трехмерными объектами в ограниченной области [6]. Для дискретизации уравнений (1.15)-(1.20) как по времени так и по пространству будем использовать центральные разности, имеющие второй порядок аппроксимации по пространству и времени [12]. Выражение для первой производной по времени будет иметь вид:

п+1/2 _ п-1/2

= +0{А12). (1.21)

Применим теперь конечные разности к уравнению (1.15). Производные заменим конечными разностями в соответствии с расположением полей в решетке Йи. Тогда получим:

тр 1П+1/2 _ р |П—1/2 гг 1П _ тт \п

хН+\/2,з,к хп+\/2,з,к __1___2п+\/2,з+\/2,к *н+\/2,з—\/2,к

Л* < £,-+1/2,з,к >

Ч \п — и \п

ПУН+1/2^,к+1/2 ПУЧ+1/2,],к-\/2 ,.„ „ |П ч

-д^--1г+1/2,.7,& < &Е{+\/2,з,к > ^х\{+\/2,з,к) •

(1.22)

Заметим, что все значения в правой части берутся по времени в точке п, но на предыдущих шагах Ех вычислялось не в точке п, чтобы не хранить лишние шаги по времени сделаем замену:

г, 171+1/2 р | тг—1/2

г, \п _ ^1г+1/2^А: ^ ^¿+1/2,^ (л 0<>ч

2чг+1/2,./,А; ~ 2 '

Подставив выражение (1.23) в (1.22), приведя подобные слагаемые и выразив значение на следующем временно шаге, в итоге получаем:

д*

г, 171+1/2 _ ( 2 < £<+1/2,з,к > \ Г ,71-1/2 (_ < £{+1/2,1к>

^хи+л/<1п к ~ / ^ Л+ I х¡2-1-1 /2.?.к

хН+1/2,з,к - ^ _ < > Ы) ^ ^ _ < > '

2 < £{+\/2,з,к > 2 < £{+1/2л,к >

' И \п — И \п И \п — И \п \

п*\г+\/2,з+\/2,к пЛг+\/2,з-\/2,к , 11УН+1/2,^+1/2 пУЧ+\/2,з,к-\/2 т[п \

К Ы) + Аг

Таким же образом получим выражение для оставшихся компонент электрического поля:

< °^+1/2<к > дг Д«

,1 -

Р |»+1/2 _ ( 2 < £г,7+1/2,А; > \ „ .гг-1/2 / £»¿+1/2,к_\

у4^+1/2,к - ^ _ < > Аг) *укт/2,к + ^ _ < >А1)'

2 < £¿¿+1/2,к > 2 < £{¿+1/2,к >

/ И \п —И \п И \п — И \п \

( 1г,7+1/2,к+1/2 Пх\^+1/2,к-1/2 дг1»+1/2^+1/2,к П *Ч-1/2,з+1/2,к ,|п Л

^ Дг До; ^ 1^+1/2,^'

(1.25)

_ < <УЕы+1/2 > А*

Г 1П+1/2 _ / 2 < к+1/2 > \ р 1П—1/2 / < ^¿¿+1/2 > \

^к^+1/2 - ^ _ < аЕи к+1/2 > Аг ) + ^ _ < аЕ^к+1/2 >А1)'

2 < Ег,з,к+ 1/2 > 2 < £¿¿,¿+1/2 >

I Д^ ^ к+1/2

(1.26)

Аналогично получим явную схему для нахождения компонент магнитного поля:

^ _ < Р"Яг,¿+1/2,к+1/2 > Аг

тт |п+1 _ ( 2 < /¿¿¿+1/2,Л;+1/2> \ гт \п ,

х 'г¿+1/2,А;+1/2 " ^ _ < ^+1/2^+1/2 > А*/

2 < ^»¿+1/2,Л+1/2 >

Аг

{ < Мг,7+1/2,£+1/2 > \

< 0ям+1/2,„+1/2 > Ау

2 < /-4.7+1/2,£+1/2 > р .п+1/2 _ р .п+1/2 ^ 1П+1/2 _ р .п+1/2

( ^УН,з+1/2,к+1 ^УЧ,5+1/2,к _ гН,з+\,к+\/2 к+1/2 _ т .п+1/2 А

Д АХ АУ ^ Н,з+\/2,к+1/2^ '

тт т+1

11У Н+1/2^,к+1/2

1

<

_— у I •*•/ — ц у» I - / —_

2 < №+1/2,з,к+1/2 > ^ „ |п ^ и п¿-ц/2,^+1/2 ^ ^^ / 2 < Мг+1/2^',А;+1/2 >

лг

< №+1/2,з,к+1/2 >

2 < №+1/2,з,к+1/2 > р 1П+1/2 _ 1П+1/2 г! 1П+1/2 _ р .п+1/2

гН+1,3,к+1/2 гН,з,к+1/2 _ хп+1/2^,к+1 х\{+1/2,з,к Ах Аг

(1.28)

.1 -

ТТ 171+1

1 -

2 < №+1/2^+1/2,к > \ гт 171 ,

^ ^1+1/2,2+1/2,к /

2 < ^{+1/2,3+1/2,к > Аг

< №+1/2,3+1/2,к >

р 171+1/2 _ р -р \ть+-х./4 _ р 171+1/^

^1г+1/2,^+1^ ¿+1/2,7,/г ^1г+1,^+1/2,/с 1г,7+1/2,А;

171+1/2

2 < №+1/2,3+1/2,к > 71+1/2

171+1/2

Ду

Да;

- /

171+1/2

МН+\/2,з+\/2,ку ' (1.29)

На основе представленной схемы реализован комплекс программ численного решения трехмерной системы уравнений Максвелла ЕМ\¥8о1уегЗО.

В работе [1] получено следующее уравнение численной дисперсии РОТБ-мстода в трехмерном случае для £ = ц — 1, &е — = 0:

Список литературы диссертационного исследования кандидат наук Семенов, Алексей Николаевич, 2013 год

Литература

1. A.Taflove and S.C.Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method.Norwood, MA:Artech, 2000.

2. J.P.Berenger, "A perfectly matched layer for the absorption of electromagnetic waves, "J.Comput.Phys., vol.114, no.l, pp.185-200, 1994.

3. Z.S.Sacks, D.M.Kingsland, R.Lee, and J.F.Lee, "A perfectly matched anisotropic absorber for use as an absorbing boundary condition, "IEEE Trans.Antennas and Propagation, vol.43 ,no.l2, pp.1460-1463, 1995.

4. Kane Yee. "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media". Antennas and Propagation, IEEE Transactions on 14: 302-307

5. S.D. Gedney (1996). "An anisotropic perfectly matched layer absorbing media for the truncation of FDTD latices". Antennas and Propagation, IEEE Transactions on 44: 1630-1639.

6. K. R. Umashankar and A. Taflove, "A novel method to analyze electromagnetic scattering of complex objects,"IEEE Trans. Electromagn. Compat., vol. 24, pp. 397-405, Nov. 1982.

7. A. Taflove and M. E. Brodwin, "Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell's equations,"IEEE Trans. Microwave Theory Tech., vol. 23, pp. 623-630, Aug. 1975.

8. A. Taflove, J. Dabkowski, and M. Genge, "Mitigation of buried pipeline voltages due to 60 Hz ac inductive coupling. Part I: Design of joint rights

of way; Part II: Pipeline grounding methods,"IEEE Trans. Power Apparatus Systems, vol. 98, pp. 1806-1823, Sept./Oct. 1979.

9. А. Зеленин. "Решение уравнений Максвелла методом FDTD."// По материалам сайта http://zfdtd.narod.ru

10. А.Н. Боголюбов, И.А. Буткарев, Ю.С. Дементьева. "Численное моделирование двумерных фотонных кристаллов."// "ЖУРНАЛ РАДИОЭЛЕКТРОНИКИ'^ 11,2006

11. David М. Hockanson // Perfectly Matched Layers Used as Absorbing Boundaries in a Three-dimensional FDTD Code

12. Самарский A.A. Теория разностных схем. M.: Наука, 1983.

13. Yang Нао, Raj Mittra - FDTD Modeling of Metamaterials Theory and Applications Artech House Publishers, 2009

14. Cui, Tie Jun; Smith, David; Liu, Ruopeng (Eds.), "Metamaterials: Theory, Design, and Applications 2010, XXIV, 368 p.

15. В. Г. Веселаго, "Электродинамика веществ с одновременно отрицательными значениями е и ц"// УФН. - 1967. - Т.92 - N7 - С.517

16. В. Г. Веселаго О формулировке принципа Ферма для света, распространяющегося в веществах с отрицательным преломлением. // УФН. — 2002. - Т. 172. - № 10. - С. 1215.

17. D. Schurig, J.B. Pendry, D.R. Smith, Opt. Express 14, 9794 (2006)

18. R. W. Ziolkowski and E. Heyman, "Wave propagation in media having negative permittivity and permeability Phys. Rev. E, 64, 056625, 2001.

19. Semenov A. N., Smirnov A. P., Ignatyeva D. O., Sukhorukov A. P., "Mathematical Modeling of an Open Microcavity with a Layer of Metamaterial"// Izvestiya Rossiiskoi Akademii Nauk. Seriya Fizicheskaya, 2011, Vol. 75, No. 12, pp. 1741-1744.

20. Dolling G., Wegener M., Soukoulis C.M., Linden S. "Negative-index metamaterial at 780 nm wavelength."// Opt Lett, 2007, 32:53-55

21. Ландау Л.Д., Лившиц E.M., "Электродинамика сплошных сред"// Наука, 1992

22. Simpson. I I . and A Taflove. "Three-dimensional FDTD modeling of impulsive ELF propagation about the Eanh-sphere."IEEE Trans. Antennas Propagai, Vol. 52, 2004. pp 443-451

23. Simpson. J J. and A. Taflove. "Efficient modeling of impulsive ELF antipodal propagation about the Earth-sphere using an optimized two-dimensional geodesic FDTD grid."IEEE Antennas Wireless Propagai Leu.. Vol. 3. 2004. pp. 215-218.

24. Huang. Y . Simulation of Semiconductor Materials Using FDTD Method. M S thesis. Northwestern Univ.. Evanston. IL, 2002.

25. Hagness. S C. A. Taflove. and J. E Bridges. 'Three-dimensional FDTD analysis of a pulsed microwave cocfocal system for breast cancer detection design of an antenna-array element "IEEE Trans Antennas Propagai.. Vol 47. 1999. pp. 783-791.

26. Bond. E J. X Li. S. C. Hagness. and В D Van Veen. "Microwave imaging via ipace time beamforming for early detection of breast cancer.** IEEE Trans Antennas Propagat. Vol 51. 2003. pp 1690-1705

27. Simpson. J. J.. A Tafiove. i A Mu. and H. Heck. "Computational and experimental study of a microwave electromagnetic bandgap structure with waveguiding defect for potential use as a bandpass wireless interconnect.** IEEE Microwave Wireless Components Lett. Vol 14. 2001. pp 345-345

28. Simpson. J. J.. A Tafiove. J. A. Mix. and H Heck. "Advances in hypcrspecd digital interconnects using electromagnetic bandgap technology: Measured low-loss 43-GHz passband centered at 50 GHz." Proc IEEE Antennas Propagai. Soc Inti Svmp . Washington. D.C . 2005

29. Park. H.-C., S.-H. Kim. S.-H Kwon. Y.-O. Ju. J.-K. Yang. J.-H. Back. S.-B Kim. and Y -II Lee, "Electrically driven tingle-cell photonic crystal laser,"Science, Vd 305. 2004. pp 1444-1447

30. Yanik. M F . S Fan. M. Soljactc. and J. D. Joannopoulos. "All-optical transistor action with bistable switching in a photonic crystal cross-waveguide geometry."Optics Lett., Vol. 28. 2003. pp 2506-2508.

31. Ильинский А.С. Распространение электромагнитных волн в нерегулярных волноводах переменного сечения. М.: Изд-во Моск. университета, 1970.

32. Ильинский А.С., Галишникова Т.Н. Математическое моделирование процесса отражения плоской электромагнитной волны от волнистой поверхности // Радиотехника и электроника. 1999. Т. 44, № 7. С. 773-786.

33. Захаров Е.В., Дмитриева И.В., Орлик С.И. Уравнения математической физики // Академия Москва, 2010, 320 с.

34. А.Н.Тихонов, А.А.Самарский, Уравнения математической физики // Издательство: МГУ, 2004.

35. Е.А. Глазкова, H.H. Попова, Анализ эффективности гибридного параллельного программирования на примере системы Blue Gene/P // Научный сервис в сети Интернет: масштабируемость, параллельность, эффективность: Труды Всероссийской суперкомпьютерной конференции (21-26 сентября 2009 г., г. Новороссийск). - М.: Изд-во МГУ, 2009. - 524 с.

36. Ardavan F. Oskooi, David Roundy, Mihai Ibanescu, Peter Bermel, J. D. Joannopoulos, and Steven G. Johnson, "MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method," Computer Physics Communications 181, 687-702 (2010).

37. A. V. Zakirov, V. D. Levchenko CFmaxwell - The Program Code for FDTD Modelling of Very Large Size Problems //in Progress in Electromagnetics Research Symposium, Stockholm, Sweden, August 2013

38. Wenhua Yu, Raj Mittra, Tao Su, Yongjun Liu, Xiaoling Yang "Parallel Finite-Difference Time-Domain Method"// Artech House; 1 edition (2006)

39. Yu. G. Smirnov, Yu. V. Shestopalov and E. D. Derevyanchuk, "Permittivity Reconstruction of Layered Dielectrics in a Rectangular Waveguide From the Transmission Coefficients at Different Frequencies,"in Inverse Problems Workshop, Karlstad, Sweden, May 2012.

40. Семенов A. H., Смирнов А.П., Игнатьева Д. О., Сухоруков А.П. Математическое моделирование открытого микрорезонатора со слоем метамате-риала // Известия РАН. Серия физическая 75(12), 2011, с. 163 7-1640

41. Семенов А.Н., Смирнов А.П. Численное моделирование уравнений Максвелла с дисперсными материалами //Математическое моделирование, 2013, том 25, номер 12, с. 19-32

42. Smirnov A.P., Semenov A.N., Ignatyeva, D.O., Sukhorukov. A.P. Full-wave modeling of open subw avelength resonator with metamaterial // Progress in Electromagnetics Research Symposium , Moscow, 2012, pp. 1254-1258

43. Smirnov, A.P., Semenov, A.N. Full wave Maxwell's equations solver EMWSolver3D // Progress in Electromagnetics Research Symposium, Moscow, 2012, pp. 252-255

44. Smirnov A. P., Semenov A.N., Shestopalov Y. Modeling of Electromagnetic Wave Propagation in Guides with Inhomogeneous Dielectric Inclusions // Inverse problems and large-scale computations/ Eds. Larisa Beilina, Yury V. Shestopalov, Springer

45. Smirnov A.P., Semenov A.N., Shestopalov Y. FDTD Simulation of Waveguide with Non-uniform Dielectric Slab // Progress in Electromagnetics Research Symposium , Stockholm, 2013, pp. 76-78

46. Семенов A. H. Параллельная реализация численного решения уравнений Максвелла Fdtd методом для больших задач // Научный сервис в сети Интернет: экзафлопсное будущее: Труды Международной суперкомпьютерной конференции М.: Изд-во МГУ, 2011. с.530-533

47. Smirnov А.P., Semenov A.N., Shestopalov Y. Modeling of electromagnetic wave propagation in guides with inhomogeneous dielectric inclusions // The Workshop proceedings (tentative title "Inverse Problems and Large-Scale Modeling with Applications in Electromagnetics. Proceedings of conferences supported by the Visby Program"), Karlstad, Sweden, 2012, pp. 13-14

48. Семенов A.H., Смирнов А.П. Комплекс программ численного решения уравнений Максвелла EMWSOLVER3D //VI Всероссийская конферен-

ция "Радиолокация и радиосвязь" , JRE (Журнал Радиоэлектроники) -ИРЭ им. В.А.Котельникова РАН Москва, 19-22 ноября 2012 г, с. 146-150

49. Schuster, A. An introduction to the Theory of Optics. London, E. Arnold. 1904.

50. Мандельштам JI.И., Лекции по оптике, теории относительости и квантовой механике. М: Наука. 1972.

51. Smith D.R., Padilla W.J., Vier D.C. Nemat-Nasser S.C., Schultz S. // Physical Review Letters. Vol. 84. No. 18. pp. 4184(3). 2000.

52. Shalaev V.M.// Nature Photonics, Vol. 1, pp. 41-48. 2007.

53. Xiao S., Chettiar U.K., Kildishev A.V., Drachev V.P., Shalaev V.M. // Optics Letters. Vol. 34. No. 22. pp. 3478-3480. 2009.

54. Pendry J.B. // Phys. Rev. Lett. 2000. V. 85. P. 3966.

55. Панфилова H.O., Д.О. Сапарина, А.П. Сухоруков // Изв. РАН. Сер. физ. 2006. Т. 70. № 12. С. 1722.

56. Engheta N. // IEEE Antennas and Wireless Propagation Letters. 2002. V. 1. № 1. P. 10.

57. Д. О. Сапарина, А. П. Сухоруков // Известия вузов. Прикладная нелинейная динамика. - 2009. - Т. 17, N 3. - С. 3-16

58. Saparina, D. О.; Sukhorukov, А. Р. // Laser Physics. 2009. Vol. 19. No. 5. pp.1125-1130

59. Zhao Y., Belov P., Hao Y. // Phys. Rev. E. 2006. V. 75. P.037602.

60. P. Chatelain et al. «Billion vortex particle direct numerical simulations of aircraft wakes» // Comput. Methods Appl. Mech. Engrg. — 2008. — Vol. 197. - Pp. 1296-1304. - doi: 10.1016/j.cma.2007.11.016.