Моделирование анизотропных метаматериалов с помощью параллельной реализации метода конечных объемов для решения нестационарных уравнений Mаксвелла

Вычислительные технологии

Том 14, № 3, 2009

Моделирование анизотропных метаматериалов с помощью параллельной реализации метода конечных объемов для решения нестационарных уравнений Ыаксвелла*

Л. Ю. ПРОКОПЬЕВА Учреждение Российской академии наук Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: ludmila.prokopeva@ict.nsc.ru

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

Ключевые слова: уравнения Максвелла, метод конечных объемов, параллельные вычисления, анизотропная среда, метаматериалы, гиперлинза.

Введение

Одно из актуальных направлений для приложения численных методов решения уравнений Максвелла обусловлено развитием экспериментальных и теоретических исследований в области нанофотоники. Сейчас разработки наноструктурированных оптических материалов с нетривиальными свойствами для создания новых современных оптических приборов ведутся во многих научно-исследовательских лабораториях. В последние годы исследователям удается получать материалы с кардинально новыми электромагнитными свойствами, не наблюдаемыми ранее в природе и даже не известными прежде классической электродинамике, а потому названными метаматериалами (от греч. “мета” — за пределами, сверх). Интересно, что родоначальником одного из таких направлений стал российский ученый В.Г. Веселаго, впервые обобщивший классические законы физики для материалов, имеющих отрицательный показатель преломления [1]. Эти его работы получили широкое признание только после экспериментального подтверждения существования таких материалов с применением металлических нанорезонаторов, периодически расположенных в диэлектрике [2-9]. Другой, не менее революционной, развивающейся технологией с использованием метаматериалов, имеющей непосредственное отношение к данной работе, является попытка преодолеть дифракционный предел

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 09-01-00352-а) и Интеграционного проекта СО РАН № 113.

© ИВТ СО РАН, 2009.

оптической линзы [10- 14]. И, наконец, самым любопытным, пожалуй, экспериментом в области метаматериалов оказался проект создания “плаща-невидимки” (cloaking device). Основная идея последнего заключается в создании такой наноструктурированной оптической среды вокруг маскируемого объекта, которую оптический импульс проходил бы без поглощения и отражения, фактически огибая объект, что сделало бы его невидимым для стороннего наблюдателя [15-20].

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

На данный момент в лабораториях, занимающихся исследованиями метаматериалов, для проведения численных экспериментов чаще всего прибегают к использованию коммерческого программного обеспечения, такого как, например, COMSOL MULTIPHYSICS или ANSYS, основанных на методе конечных элементов. Однако современные требования к численным расчетам включают в себя возможность сеток эффективно покрывать сложные объекты с нетривиальной геометрией, адекватное описание физических процессов в композитных средах диэлектрик—металл, возможность проведения вычислительного эксперимента с использованием многопроцессорных вычислительных систем. В связи с этим повышенное внимание в последнее время стало уделяться методам численного решения нестационарных уравнений Максвелла на неструктурированных сетках, таким как метод конечных объемов [21-25].

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

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

Во-первых, в конечно-объемной формулировке, следующей из формулы Гаусса— Остроградского, потребуется вычислять значение потока поля на ребрах (гранях — для трехмерного случая). Поскольку значения, полученные аппроксимацией дискретного поля из центров правого и левого конечных объемов, могут сильно различаться, особенно на границах раздела сред, то для нахождения потока точно решается одномерная задача Римана вдоль нормали к ребру (грани).

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

В-третьих, для моделирования реальных физических задач предсказать даже дальнее поле1 удается не всегда, и чтобы проводить расчеты в ограниченной области, требуется аккуратная постановка поглощающих граничных условий. В конечно-объемные методы техника введения дополнительного поглощающего слоя (PML) пришла из конечноразностных методов [27], и эффективность ее различных модификаций активно обсуждается, например, в [28]. Однако в данной работе успешно использовались более простые локальные граничные условия (Silver-Mueller condition), в общем случае имеющие первый порядок точности.

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

В-пятых, проблема, связанная с возможным недостатком вычислительных ресурсов при необходимости расчетов сколько-нибудь существенного образца наноструктуриро-ванного метаматериала, фактически решена с помощью разработанной параллельной версии метода конечных объемов [29]. Параллельный алгоритм основан на декомпозиции расчетной области на подобласти, каждая из которых предназначена для проведения расчетов на отдельном процессоре, а обмен необходимыми значениями между процессорами реализован с помощью библиотеки MPI, что позволяет выполнять вычисления на современных многопроцессорных системах различной архитектуры. Таким образом, вопрос о вычислительной сложности упирается лишь в возможности современных суперкомпьютеров, производительность которых, надо сказать, вот уже 40 лет четко следуя “закону Мура”, удваивается каждые 18 месяцев, и уже превосходит 1015 операций с плавающей запятой в 1 с (“RoadRunner” [30]).

Разрабатываемый численный метод конечных объемов для решения нестационарных уравнений Максвелла тестировался на задачах дифракции плоской электромагнитной волны на металлическом цилиндре и диэлектрическом градиентном цилиндре [31]. В данной работе метод обобщен для анизотропных сред и с его помощью проведено моделирование цилиндрических гиперлинз, разработкой которых в последнее время занимается одна из ведущих в исследовании метаматериалов лаборатория в университете Purdue (США, штат Индиана). Результаты расчетов, о которых пойдет речь в этой работе, были представлены на конференции [32], здесь же подробно излагаются применяемый алгоритм и численный метод.

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

1. Гиперлинза

Одной из разработок в области метаматериалов является изобретение под названием гиперлинза — прибор, основанный на использовании анизотропных метаматериалов, который способен различать предметы менее половины длины волны видимого света (дифракционный предел) и проецировать изображение в дальнее поле, где с помощью стандартной оптики изображение может быть увеличено далее до нужных размеров. Оригинальная идея прибора описана в [10, 11], последовавшие экспериментальные подтверждения работы прибора можно найти, например, в [12, 13]. Однако в приведенных работах недостатком предлагаемого дизайна являлось то, что на границах гиперлинзы импеданс п = \/№*/£ претерпевал существенный скачок. Это приводило к частичному отражению оптического сигнала и значительно уменьшало его интенсивность на выходе. Поэтому позже в [14] были разработаны модели “внутренней” и “внешней” гиперлинз, построенные на условии равенства импедансов на внешней и внутренней границах прибора соответственно и при этом являющиеся немагнитными средами (изготовление метаматериалов, имеющих магнитный отклик на оптических частотах, представляет собой отдельную, сложную в эксперименте, задачу). В данной работе будет проведено моделирование этих приборов и численно продемонстрирована их работа.

В модели цилиндрическая гиперлинза находится в вакууме (ег = 1) и представляет собой бесконечный в направлении оси г цилиндр (в сечении изображен на рис. 1) с выколотой сердцевиной:

{

Аналитическое выражение для диэлектрической проницаемости гиперлинз выводится в [14]. Для внутренней гиперлинзы условие равенства импедансов при р = а приводит к соотношениям

єР

Аналогично, для внешней гиперлинзы из условия равенства импедансов при р = I соотношения примут вид

Єр

Здесь, как ив (2), єр,Єф — ненулевые компоненты тензора относительной диэлектрической проницаемости в цилиндрических координатах єг = ^ Є Є ^ , (Пр,Пф)Т =

Рис. 1. Гиперлинза: внутренний радиус а = 600 нм, внешний радиус I = 3 мкм; на окружность, отмеченную штриховой линией, помещены источники

£0£г (Ер,Еф)Т; г(р) = т 1(р — /) + Ь — выбранное линейное преобразование координат, которое взаимно однозначно отображает кольцо (а < р < /} на кольцо (а < г < Ь},

причем г(а) = а, г(/) = Ь; коэффициент т = (/ — а)/(Ь — а); параметры линзы а = 600 нм, Ь = 610 нм, I = 3 мкм (рис. 1).

Далее, на дугу окружности

помещаются пять жестких когерентных монохроматических источников света с длиной волны Л = 632 нм и амплитудой равной:

Теоретически такая конструкция должна точно проецировать пять источников с интерфейса р = а на р = /, при этом, поскольку импеданс не терпит скачка только на одной границе (внутренней либо внешней границе гиперлинзы), следует ожидать, что часть энергии светового излучения отразится от другой границы (внешней либо внутренней соответственно). Для того чтобы качественно оценить работу приборов и уровень потерь вследствие отражения, смоделируем распространение оптического импульса (4) во внутренней и внешней гиперлинзах (2), (3) численно, используя метод конечных объемов во временной области.

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

В общем случае распространение электромагнитных волн описывается полными уравнениями Максвелла, кроме того, для оптической немагнитной среды сразу будем предполагать ^г = 1 и отсутствие свободных зарядов Ю = 0:

| (р, ф) | р = а — 30 нм,

k=0 4 18 ’ 4 18 , 3п/4 < ф < 5п/4.

—— rot H = 0, где D(w) = e0£rE,

< dB

——+ rot E = 0, где B = u0H,

div D = 0, div B = 0.

\

В данном случае задача имеет решение в виде ТМ-волны E = (Ex, Ey, 0), H = (0 , 0, Hz), для которой уравнения (5) в матричной форме примут вид

где

Рис. 2. Треугольная сетка, адаптированная к геометрии линзы

Здесь Є**1 — компоненты обратного тензора относительной диэлектрической проницаемости в декартовых координатах, которые могут быть выражены через известные Єр, Єф

(2), (3):

_-, cos2 ф sin2 ф cos2 ф sin2 ф /1 1 \

єхх =—------ + —----, ЄУУ=—------- + ~----, Єху =sin ф COs ф

Єр Єф Єф Єр \ЄР ЄФ /

Далее в расчетной области r < 4 мкм строится неструктурированная треугольная сетка, учитывающая геометрию гиперлинзы и форму источников (рис. 2)2.

Для построения численной схемы уравнения (6) интегрируются по каждой треугольной ячейке Ai и преобразуются к контурному интегралу по формуле Гаусса—Остроградского:

dtj Udxidx2 + AUdr = 0,

A¿ k=1rk

где Гк — ребра ячейки, А = п1А1 + п2А2, п = (п1; п2) — внешняя нормаль к ребру Г Наконец, аппроксимируя уравнение (9), получаем окончательную схему:

ттп+1 ттга 3

Б ---1 + Ё и,к = 0,

Т *—' ’

к=1

где —П — значения поля в барицентре хь треугольного элемента Аі в момент времени Ьп = пт; т — шаг по времени; Еі^ — значение потока Аи в центре к-го ребра г-го треугольника; іі^ — длина ребра; Бі — площадь г-го треугольника.

Неструктурированная сетка была построена А.С. Лебедевым, которому автор выражает свою благодарность.

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

д д

- и + а—и = 0,

дЬ дп

и, следовательно, вычисляется по формуле

?г,к = А+Иь + А-Ия,

где

А

±

Здесь Б — матрица, столбцы которой являются правыми собственными векторами матрицы А; \к — собственные числа матрицы А; Иь,я — значения поля И, аппроксимированные на ребро из левого и правого примыкающих треугольников по формуле Тейлора, при этом потребуются значения градиента поля (градиент поля на неструктурированной сетке аппроксимируется специальным образом, подробно описанным в

[26]).

В данном случае (7) матрицы А± примут вид:

/ асп2

А+ = -

5

асщ

~

а

е0

Ьсп2

Ьсп\

~

Ь

е0

п2 \

^о

Щ

^0

с5

А-

/ асп2

~Т~

асщ

5

а

е0

Ьсп2

~

Ьсщ

Т~

Ь

е0

п2

^0

щ

^0

— — с5

где а = п,ех!/‘ — Ь = п1буу1 — ще,-', 5 = ^ п1е— — 2п1п2е-У + ф-.

Однако из численного эксперимента следует, что для разрывных сред схема будет работать точнее, если аппроксимацию значений выполнять в непрерывных переменных. В данном случае непрерывными являются компоненты V = (Ор ,Еф ,ИХ) (нормальная к линии разрыва компонента индукции и касательная компонента поля). И затем только Иь,я = ТV ь,я использовать для вычисления потоков в (12). Здесь Т — матрица перехода, вычисленная в середине ребра:

сое ф —е0еф вт ф 0

Т

е0еф сое ф 0

0

^0

Тогда удобно ввести новое обозначение С± ваться формулой ¥г,к = С+Уь + С-Уя.

А±Т и вместо формулы (12) пользо-

1

2

3. Численные результаты

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

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

Прежде всего проведем моделирование источников (4) в вакууме (рис. 3, а). Видно, что сигнал в результате интерференции когерентных источников фокусируется, затем рассеивается. Теперь добавим внутреннюю линзу (2) и внешнюю (3) (рис. 3, б и 4, а). Как и ожидалась, помимо того, что расчет верифицирует работу цилиндриче-

а б

Внутренняя гиперлинза

х, мкм х, мкм

Рис. 3. Источники в вакууме, амплитуда поля И

х, мкм

Рис. 4. Внешняя линза, амплитуда поля И% (а); сравнение изображений от пяти источников на внутреннем и внешнем интерфейсах внешней линзы (б). Часть энергии уходит в виде отраженных от внутреннего интерфейса волн, чем объясняется шум на графике

Рис. 5. Фрагменты распространения излучения во внутренней линзе

ской гиперлинзы (изображение пяти источников из ближней зоны точно проецируется с увеличением в дальнюю зону, преодолевая тем самым дифракционный предел для Л = 632 нм), появляются также и отраженные от границ цилиндрические моды.

В результате при сравнении изображений на внутреннем и внешнем интерфейсах будут наблюдаться некоторое ослабление сигнала и возникновение шума (рис. 4, б). В целом качество работы такого прибора сравнимо с “идеальной линзой” [14], а его изготовление будет значительно упрощено за счет использования немагнитного материала.

Кроме того, анимация расчета во временной области дает визуальное представление работы прибора. Например, на рис. 5 представлена магнитная компонента поля (амплитуда которой была ранее изображена на рис. 3, б) в два выбранных момента времени. Основной сигнал здесь выглядит в виде пяти радиальных лучей, которые возвышаются над остальным шумом. Присутствие шума здесь объясняется цилиндрическими модами, отраженными от внешней границы. При просмотре анимации целиком видно, что шум представляет собой стоячие волны.

4. Заключение

Проведенные в данной работе расчеты по большей части являются демонстрационными и повторяют оригинальные расчеты с помощью коммерческого программного обеспечения из [14]. Тем самым показана эффективность метода конечных объемов решения нестационарных уравнений Максвелла при моделировании частной задачи по исследованию метаматериалов, а также продемонстрирована модификация метода, разработанного ранее в [26-29, 31-35], для анизотропных сред. Кроме того, лишь незначительные модификации двумерного метода конечных объемов требуются для численного решения полных трехмерных уравнений Максвелла на тетраэдральных сетках, что будет показано в дальнейших публикациях. Отметим также, что в трехмерных задачах будет использоваться параллельная версия алгоритма [29]. В дальнейшем разработанные параллельные методы планируется применять для моделирования электромагнитных явлений в новых, актуальных для исследования метаматериалов, конструкциях в су-

щественно трехмерных постановках.

Список литературы

[1] VESELAGO V.G. The electrodynamics of substances with simultaneously negative values of e and ^ // Sov. Phys. Usp. 1968. Vol. 10, N 4. P. 509-514.

[2] Pendry J.B. Negative refraction makes a perfect lens // Phys. Rev. Lett. 2000. Vol. 85, N 18.

P. 3966-3969.

[3] Smith D.R., Padilla W.J., Vier D.C. et al. Composite medium with simultaneously negative permeability and permittivity // Phys. Rev. Lett. 2000. Vol. 84, N 18. P. 4184-4187.

[4] Shelby R.A., Smith D.R., Schultz S. Experimental verification of a negative index of refraction // Science. 2001. Vol. 292. P. 77-79.

[5] Parazzoli C.G., Greegor R.B., Nielsen J.A. Performance of a negative index of

refraction lens // Appl. Phys. Lett. 2004. Vol. 84, N 17. P. 3232-3234.

[6] Shalaev V.M., Cai W., Chettiar U.K. et al. Negative index of refraction in optical

metamaterials // Opt. Lett. 2005. Vol. 30. P. 3356-3358.

[7] Zhang S., Fan W., Panoiu N.C. et al. Experimental demonstration of near-infrared

negative-index metamaterials // Phys. Rev. Lett. 2005. Vol. 95. Id. 137404.

[8] Kildishev A.V., Cai W., Chettiar U.K. et al. Negative refractive index in optics of metal-dielectric composites // J. Opt. Soc. Am. B. 2006. Vol. 23. P. 423-433.

[9] Zhang S., Fan W., Malloy K.J., Brueck S.R.J. et al. Demonstration of metal-dielectric negative-index metamaterials with improved performance at optical frequencies // J. Opt. Soc. Am. B. 2006. Vol. 23, P. 434-438.

[10] Jacob Z., Alekseyev L.V., Narimanov E. Optical Hyperlens: Far-field imaging beyond the diffraction limit // Opt. Express. 2006. Vol. 14. P. 8247-8256.

[11] Salandrino A., Engheta N. Far-field subdiffraction optical microscopy using metamaterial crystals: Theory and simulations // Phys. Rev. B. 2006. Vol. 74. Id. 075103

[12] Liu Z., Lee H., Xiong Y. et al. Optical hyperlens magnifying sub-diffraction-limited objects // Science. 2007. Vol. 315. P. 1686.

[13] Smolyaninov I., Hung Y., Davis C. Magnifying superlens in the visible frequency range // Science. 2007. Vol. 315. P. 5819.

[14] Kildishev A.K., Narimanov E.E. Impedance-matched hyperlens // Opt. Lett. 2007. Vol. 32, N 23. P. 3432-3434.

[15] Greenleaf A., Lassas M., Uhlmann G. Anisotropic conductivities that cannot be detected by EIT // Physiol. Meas. 2003. Vol. 24. P. 413-419.

[16] Pendry J.B., Schurig D., Smith D.R. Controlling electromagnetic fields // Science. 2006. Vol. 312. P. 1780.

[17] Schurig D., Mock J.J., Justice B.J. et al. Metamaterial electromagnetic cloak at

microwave frequencies // Science. 2006. Vol. 314. P. 977-980.

[18] Zhang B., Chen H., Wu B.—I. et al. Response of a cylindrical invisibility cloak to

electromagnetic waves // Phys. Rev. B. 2007. Vol. 76. Id. 121101.

[19] Cai W., Chettiar U.K., Kildishev A.V., Shalaev V.M. Optical cloaking with metamaterials // Nat. Photonics. 2007. Vol. 1. P. 224-227.

[20] Cai W., Chettiar U.K., Kildishev A.V., Shalaev V.M. Designs for optical cloaking with high-order transformations // Optics Express. 2008. Vol. 16. P. 5444-5452.

[21] Shankar V., Mohammadian A.H., Hall W.F. A time-domain, finite-volume treatment for the Maxwell equations // Electromagnetics. 1990. Vol. 10. P. 127-145.

[22] Bonnet P., Ferrieres X., Paladian F. et al. Electromagnetic wave diffraction using a finite volume method // Electron. Lett. 1997. Vol. 33, N 1. P. 31-32.

[23] Fumeaux C., Baumann D., Leuchtmann P., Vahldieck R. A generalized local timestep scheme for efficient FVTD simulations in strongly inhomogeneous meshes // IEEE Trans. Microwave Theory Tech. 2004. Vol. 52, N 3. P. 1067-1076.

[24] Fumeaux C., Baumann D., Bonnet P., Vahldieck R. Developments of finite-volume techniques for electromagnetic modeling in unstructured meshes // 17th International Zurich Symposium on Electromagnetic Compatibility, 2006.

[25] Firsov D. et al. High-Order FVTD on unstructured grids using an object-oriented computational engine // ACES J. 2007. Vol. 22, N 1.

[26] Лебедев А.С., Федорук М.П., Штыринл О.В. Конечно-объемный алгоритм решения нестационарных уравнений Максвелла на неструктурированной сетке // Журн. вычисл. математики и мат. физики. 2006. Т. 47, № 7. С. 1286-1301.

[27] Berenger J.-P. A perfectly matched layer for the absorption of electromagnetic waves // J. of Comput. Phys. 1994. Vol. 114, N 2. P. 185-200.

[28] Kaufmann T., Sankaran K., Fumeaux C., Vahldieck R. A review of perfectly matched absorbers for the Finite-Volume Time-Domain method // ACES J. 2008. Vol. 23, N 3. P. 184-192.

[29] Прокопьева Л.Ю., Шокин Ю.И., Лебедев А.С., Федорук М.П. Параллельная реализация метода конечных объемов для решения нестационарных уравнений Максвелла на неструктурированной сетке // Вычисл. технологии. 2007. Т. 12. Спецвыпуск 4. С. 59-69.

[30] www.top500.org

[31] Prokopyeva L.Yu., Shokin Yu.I., Lebedev A.S. et al. Parallel numerical modeling of modern optics devices // Proc. 3nd Russian-German Workshop, Novosibirsk, 2007. P. 93-107.

[32] Prokopeva L.J., Lebedev A.S., Fedoruk M.P., Kildishev A.V.

FVTD simulations of nano-structured plasmonic metamaterials / / 24th

Annual Review of Progress in Applied Computational Electromagnetics. March 30 — April 4, 2008. Canada, Niagara Falls. P. 562-566.

[33] Лебедев А.С., Федорук М.П., Штырина О.В. Решение нестационарных уравнений Максвелла для сред с неоднородными свойствами методом конечных объемов // Вычисл. технологии. 2005. Т. 10, № 2. С. 60-73.

[34] Shokin Yu.I., Lebedev A.S., Shtyrina O.V., Fedoruk M.P. Solution of Maxwell’s equation on partially unstructured meshes // Notes on Numerical Fluid Mechanics and Multidisciplinary Design. 2006. Vol. 93. P. 1-13.

[35] Fedoruk M.P., Lebedev A.S., Shokin Yu.I. Finite volume algorithm for nonstationary Maxwell equations on an unstructured grid // RJNAMM. 2007. Vol. 22, N 1. P. 1-18.

Поступила в редакцию 9 октября 2008 г., в переработанном виде — 30 декабря 2008 г.