Задание падающей волны по технологии TF/SF при согласованном разностном решении уравнений Даламбера и Максвелла Текст научной статьи по специальности «Математика»

УДК 510.8

ЗАДАНИЕ ПАДАЮЩЕЙ ВОЛНЫ ПО ТЕХНОЛОГИИ TF/SF ПРИ СОГЛАСОВАННОМ РАЗНОСТНОМ РЕШЕНИИ УРАВНЕНИЙ ДАЛАМБЕРА И МАКСВЕЛЛА

© 2014 Л.В. Яблокова1, Е.Ю. Булдыгин1, Д.Л. Головашкин2

1 Самарский государственный аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет), 2 Институт систем обработки изображений РАН, г. Самара

Поступила в редакцию 17.12.2013

Предложена методика отыскания совместного разностного решения волнового уравнения и системы уравнений Максвелла, позволяющая совместить достоинства и избежать недостатков обоих упомянутых численных методов нанофотоники. В двумерном случае на тестовых примерах продемонстрированы сходимость такого решения, возможность задания падающей волны по технологии ТЕ/БЕ Ключевые слова: уравнение Даламбера, уравнения Максвелла, разностные схемы, методика ТЕ/БЕ

ВВЕДЕНИЕ

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

Указанный недостаток принято преодолевать учетом природы электромагнитного излучения (введением подвижных сеточных областей [1]), особенностей строения элементов нанофотоники (декомпозицией сеточной области [2]), разработкой новых приемов программирования (метод пирамид [3]) или обращением к разностному решению волнового уравнения [4]. Комбинированию последнего подхода с ЕОТЭ-мето-дом и посвящена предлагаемая работа.

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

Яблокова Людмила Вениаминовна, старший преподаватель кафедры прикладной математики. E-mail: [email protected]

Булдыгин Евгений Юрьевич, студент магистратуры факультета информатики. E-mail: [email protected] Головашкин Димитрий Львович, доктор физико-математических наук, доцент, ведущий научный сотрудник. E-mail: [email protected]

электрического поля вместо трех проекций электромагнитного) и на 10% снизить число арифметических операций (10 при вычислении поля в одном узле сеточной области вместо 11 при вычислении трех сеточных функций в ячейке Уее). Впрочем, при рассмотрении трехмерной задачи метод ЕЭТЭ становится более предпочтительным, что ограничивает область применения развиваемого подхода двумерным случаем.

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

Перенесение современных методик, разработанных для решения уравнений Максвелла, на разностное решение волнового уравнения не всегда приводит к упрощению задачи. Так, в работе [5] авторы обращаются к заданию более сложного "прозрачного" источника вместо хорошо зарекомендовавшей себя технологии разделения полей ТЕ/БЕ в силу неработоспособности последней при решении уравнения Даламбера.

Поэтому представляется целесообразным совместное разностное решение уравнений Далам-бера и Максвелла, где первое привлекается для снижения требований к системным ресурсам, а участие второго позволяет использовать без модификации современные эффективные методики. В настоящей статье это методика формирования падающей волны (ТЕ/БЕ), рассмотренная

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

1. СОГЛАСОВАНИЕ РАЗНОСТНЫХ РЕШЕНИЙ УРАВНЕНИЙ ДАЛАМБЕРА И МАКСВЕЛЛА НА ГРАНИЦЕ СЕТОЧНЫХ ПОДОБЛАСТЕЙ

Для уравнений Максвелла ограничимся рассмотрением случая Н - волны:

М0

£0£

дИ дЕ дИ дЕ

, М

'2 __X

дt дг ' 0 дt ду

дЕх _дИ дИу

(1)

дt ду дг

записав разностную схему Уее:

М

Ип+0,5 ИП _0'5 ЕП ЕП

ут ,к+0,5 ут ,к+0,5 хт ,к+1 хт ,к

И

И

М0

И п+0,5 Ип _0'5 Еп ЕП

гт+0,5,к г т+0,5,к хт+1,к хт ,к

к

к

£0£т ,к

ЕП +1 Еп Ип+0,5 Ип+0,5

хт ,к хт ,к г т+0,5,к г т-0,5 ,к

(2)

к

к

Ип+0,5 ип +0,5

Ут ,к+0,5 ут ,к-0,5

и

где сеточная проекция электрического поля на

ЕП „ ,, ц-п+0,5

„ „ х , магнитного на г и у - И;

п + 0 5 хтк _ ^ гт+0,5,к

ось х

^п+0, ут, к+0,5

введём

Ип

соответственно. Решая уравнение (1)

В : (0 < t < Т,0 < у < Ь ,0 < г < Ьг)

на

которую наложим сеточную область Вк. В этом

п

случае проекцию определим в узлах

{&, ут, г к): tn _ Пк, п _ 0, 1 N _ Т,

и

Ь

Ь

ут _пку,т_0, .„М_~г, гк _Щ,к _0, ..,К_

и

и

проекцию и;т°к+50,5 - Уm,гк+0,5) : ^+0,5 _(П+ 0,5)И,

П _ 0, 1, .., N-1, ут _ тИр _ 1, .., М-1,гк +0,5 _ (к + 0,5)иг к _ 0, ..,К-1} , проекцию - {(tn+0,5,Уm+0,5,гк):

С+0,5 _(П + 0,5)И,П_0, 1, ..,N-1,уи+0,5 _(т + 0,5)иу,

т _ 0, ..,М -1, гк _ кИг,к _ 1, ..,К-1}. Выделим на Б подобласть Вс (0</<Т,Ь<у<Я,В<г<Ц),

задав на ней ВИ. Будем считать, что Е; опре-

делена в узлах

{^ , у , г ): t _ ;к, ; _ 0, 1, ..,N _ Т /к,

4 п'«7 т ' к ' п ^ ? ? ?

у _ тИ , т _ Ь,...,Я, Ь _ Ь /И , Я _ Ь /И ,

ут у ? •• • ? I у? г уу

гк _ Щ ,к _ В, ..,и,и _ Ьи / Иг, В _ Ьъ / И },

И!0,; - №;+ 0,5 , у т+0^ гк ) : ^ +0,5 _(; + 0,5 ),

П _ 0, 1, .., N -1, у т+0,5 _(т + 0,5 ) Иу ,

т _ Ь, .., Я-1, гк _ Щ, к _ В +1, .., К-1).

Тогда решение (1) по схеме (2) будет отыскиваться на ВМ = В\Вс, а сеточные уравнения схемы (2) решаться на ВМ = В\В°ь .

Аналогично для волнового уравнения:

д2их д 2их д 2их

- +-2х - М£-2х _ 0

ду2 дг2 ' д/2 запишем разностную схему:

(3)

ип +1-2и\ + ип

т, к т ,к т, к

с ип +1к - 2ип, + ип ,,

/ т+1, к т,к т-1, к

И2 £

'Ч 1

ип -2ип + ип

^ т,к + 1 ^^ тк ^ т,к-1

И2

(4)

Решение (3) будем искать в Вс сеточной области ВИ , наложенной на Вс , реализуя вычисления по (4) на сеточной области В^ . Тогда проекция и;,к определена в узлах {(/п , ут, гк ):

t _ пИ , ; _ 0, 1 ..,N _ Т / И ,

ут _ тИу, т _ Ь, ..,Я,Ь _ Ь1 /Иу, Я _ Ьг / Иу,

г. _ кИ , к _ В, ..,и, В _ Ьъ / И ,и _ Ь / И }.

к г' ъ г' и г-*

Задавшись целью совместного отыскания разностного решения уравнений (1) и (3), согласуем вычисления на Б№ и Бм в табл. 1.

На границах объединенной области установим электрическую стенку. За начальное условие примем отсутствие поля в момент времени t _ 0.

Алгоритм перехода на следующий временной слой при реализации совместного решения состоит в определении сеточных функций в узлах В^

Таблица 1. Равные значения сеточных функций Е; к и и;к на границах сеточных областей

" ' вМ и ВУ

Границы области Узлы сеточных областей

левая граница ВУ т _ Ь, В < к < и

глУ правая граница Вк т _ Я, В < к < и

верхняя граница ВИ к _ и, Ь < т < Я

глУ нижняя граница ВИ к _ В, Ь < т < Я

Таблица 2. Расчёт совместного решения (1) и (3) на временном слое n+1

Узлы сеточных областей Подобласти вычислительных областей Решаемые уравнения

m _ 0,0 < k < K левая граница вкМ e;,, = 0

k _ 0,0 < m < M нижняя граница В^ E" = 0 Xm к

k _ K,0 < m <M верхняя граница ВМ E" = 0 xm ,K

m _ M,0 < k < K правая граница ВМ E" = 0 XM ,k

0 < m < L, R < m <M, 0 < k < K область ВМ (2)

L < m < R,0 < k < B,U < k < K область БМ (2)

L < m < R, B < k < U область В^ (4)

m _ L, B < k < U левая граница Вк (5)

m _ R, B < k < U правая граница Вк (6)

k _ B, L < m < R нижняя граница Вк (7)

k _ U, L < m < R верхняя граница Вк (8)

m _ L, k _ B нижний левый узел Вк (9)

m _ L, k _ U верхний левый узел Вк (10)

m _ R, k _ U верхний правый узел Вк (11)

m _ R, k _ B нижний правый узел Вк (12)

и Dh в соответствии с табл. 2, где формулы (5) (12) имеют вид:

U"+1 — 2U" + U"—1 U" — 2U" + E"

Um,к m,к U m,к _ C ,U m + 1,k ¿Um,k m—1,k +

Tb. _ V Tb.

T jn+1 _jn Jjn—1 Jjn —IT!" 4- T?"

m ,k m,k m, k C ^ m + 1,k m Em—1,k +

m, k __(

U", +l — 2U", + U\ ,

+ m ,k +1 m, k m ,k—1 )

i 2

(5)

Tjn+1 _r)T jn \Tjn—1 1ч" —IT jn A-Tjn

m ,k LKJ m ,k m, k C /JZm+1,k LKJ m, k U m—1,k

U". +1 —2U". + U". 1

+ m ,k +1 m, k m ,k—1 )

i 2

(6)

UTk —2Um,k +Um, k1 c Um+1,k—2Um,k +Um—1,k

—+

U", +1 — 2U", + E", 1

+ m ,k +1 m, k m ,k—1 )

2

(7)

Um,k—2Um,k +Um,к _ c u+1,k—2Um,k +Um—u +

i 2 _ V i 2

E" k+1 —2U"k + U"k 1

+ m, k +1 m, k m ,k—1 )

i 2

(8)

TT "+1 — it t" 4-ti" tj" —itj" 4- Г"

m ,k LKJmkk^Um, k C / m+1, k m,k ^m—1,k

U"k +1 —2U"k + E\ 1

+ m, k +1 m, k m ,k—1 )

(9)

E" k+1 —2U"k + U"k 1

+ m, k+1 m,k m,k—1 )

i 2

(10)

U"+1 —2U" + U"—1 E" —2U" + U"

m,k LKJ m,k m,k C /JZm+1,k m,k^ U m—1,k

h 2

h2

E" k+1 —2U"k + U"k 1

+ m,k+1 m,k m, k—1 )

i 2

(11)

и"+1_9/"t" a_tj"~1 f" - it!" 4- tj"

m,k LKJ m,k m,k C ^m+1,k m,^ U m—1,k

h 2

h2

U"k+1 —2U", + E", 1

+ m,k +1 m,k m,k—1 )

(12)

h

Для подтверждения правомерности предложенного подхода к совместному решению на одной вычислительной области уравнений Максвелла и волнового авторы провели серии вычислительных экспериментов. При этом использовались операционная система Windows 7 Professional SP1, вычислительный комплекс Matlab 7.0.1 и процессор Pentium (R) Dual - Core CPU T 4400 @ 2.20 GHz, RAM 3.00 CB.

На рис. 1 уравнения Максвелла решаются относительно сеточных функций Hz (обозначены треугольниками), H (крестиками) и Ex (закрашенными квадратиками); уравнения Даламбе-

Рис. 1. Объединение сеточных областей без учета дискретизации по времени.

ТУ С" и" + °>5 и*

Квадратам соответствуют проекции , треугольникам - н г , крестикам - н у

ра - относительно сеточных функций Ех (не закрашенными квадратиками); на общей границе сеточных областей расположены значения фун- кого источника. Их совпадение свидетельствует

распространении плоской однородной электромагнитной волны в вакууме при задании жёст-

кции напряжённости электрического поля (наполовину закрашенные квадратики).

Эксперименты проводились при различных значениях дискретизации сеточной области О, й иQT, где первый параметр характеризовал число узлов сеточной области по пространству (приходящееся на одну длину волны); второй - количество узлов по времени (приходяще-

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

2. ЗАДАНИЕ ПАДАЮЩЕЙ ВОЛНЫ ПО ТЕХНОЛОГИИ TF/SF ПРИ СОВМЕСТНОМ РЕШЕНИИ

Моделирование распространения излучения

еся на временной интервал, за который плоский через оптический элемент кроме наложения се-

волновой фронт в вакууме пройдет расстояние в одну длину волны); третий - "длительность" запускаемого цуга в длинах волн. При этом они менялись от (10,20,5) - минимальных значений, удовлетворительно описывающих распростране-

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

ние плоской однородной волны в свободном про- торого он изготовлен, но и от вида падающей

странстве, до (100,200,15) - соответствующих электромагнитной волны - распределения ком-

весьма низким величинам погрешности. В каче- плексных амплитуд пр°ещий её вект°р°в в пр°-

стве вычислительной области О брался квадрат странстве и времени.

с длиной стороны 20Х , при этом шаги по пространству полагались равными.

В первой серии экспериментов были повторены результаты из работы [5], полученные при

Непосредственное использование модели жесткого источника (вполне приемлемой в одномерном случае) в двумерном варианте оказывается неудачным, при необходимости корректной ра-

боты с отраженной от оптического элемента волны. Наиболее популярным приемом, в этом случае, является задание падающей волны по технологии ТР/БР, связанной с искусственным разделением поля на результирующее (в оптическом элементе и его непосредственной окрестности) и рассеянное (в остальной области). Выражения, описывающие указанное разделение полей, содержат слагаемые с аналитически (или численно) заданным падающим полем. Таким образом, происходит учёт последнего.

Реализация методики ТБ^Б связана со следующей модификацией разностной схемы Уее:

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

~п+0,5 ~п-0,5

НУк+0,5 - НЛ+0,5 ЕХщ - Е.

п ~п

'к+1 Е Хк

к

к

к = 0, К -1, п = 0, N -1,

~п+1 ~п -—п+0, 5 -—■ п+0,5

Е Хк — Е Хк Н у к+0,5 — Н у к-0,5

к

к

(14)

ЕТ1 - ЕХп нп+0,5 - (нп+0,5 + Н^ ;;"к)

Х1-1,к Х1-1, к 21-0,5,к 4 21-1,5,к 1 ;, ,к

п+0,5 гттп+0,5

к

к

к = 0, К, к = 0, N -1, Е в-1 = Бт а>пк.

Н п+0,5 - н п+0,5

у1-1,к+0,5 у1-1 ,к-0,5 .

к ' Е-;1 - Епх (щ +0,5 + Н У++1'5к) - н;+0,5

ХК+1,к ХК+1,к 4 гК+1,5,к -■к+1,5,к/ гк+0,5,к

к

Я п+0,5 ттп+0,5

у - н у

уК+1,к+0,5 уК+1 ,к-0,5

к

к

т-п+1 т?п

Еп+1 Еп Н п +0,5 Н п+0,5

Хт ,и+1 Хт ,и+1 2т+0,5,и+1 2т-0,5,и+1

к

к

у

(тгп+0,5 +Ттп+0,5 \_ тт п+0,5

(Нут,и+1,5 + Нут,и+1,5) Н у^

к

Еп+1 Еп Н п+0,5 Н п+0,5

Хт ,В-1 Хт ,В-1 2т+0,5 ,В-1 2т-0,5,В-1

(13)

к

к

тт п+0,5 _ / тт п+0,5 +ттп+0,5 \

Н утВ-05 (Нут,В-0,5 + Н утВ-15 ) _

к '

щ+0 5 - Нп-0,5 (Еп - ЕХ11 к) - Е"х

21-1,5,к 21-1,5,к 4 Х1-1,к 1 1,к ' Х.

к

к

г п+0,5 и п-0,5 т?п

Н п+0,5 НГ

2Я+1,5 ,к 2Я+1,5,к ХЯ+2,к 4 ХЯ+1,к

Еп (Еп „„,, Е Х«+1,к )

к

к

ттп+0,5 _ тт п-0,5 / трп _ ТГ ^ — рп

Пут,В-Х,Ь Пут,В-Х,Ь (ЕХт,В-1 Е Хт,В-1 ) ЕХт,В-2

к

Н п+0,5 Н п-0,5 ут ,и+1,5 ут,и+1,5

к

ЕХ - (Е" - ЕХт и+1 )

Хт,и+2 4 Хт ,и+1 т ,и +1 '

к

к

Замысел второй серии экспериментов состоял в верификации возможности задания падающей волны по методике ТБ^Б при совместном решении уравнений Максвелла и Даламбера.

Усложняя задачу (по сравнению с первой серией экспериментов) исследуем дифракцию на бесконечном однородном диэлектрическом цилиндре кругового сечения, совместив его центр с центром области , определив радиус равный половине длины волны, показатель преломления п = 1,5. Остальные параметры вычислительного эксперимента остались аналогичными первой серии. Сравнивая полученные результаты с известным аналитическим решением [6], для оценки погрешности воспользуем-

ся величиной 8

(£„ = тах

4 К к

Л - N

, N

зна-

чения, полученные разностными методами, Лк -соответствующие аналитическому решению), отыскивая её на оптической оси элемента варьируя к = В, и и т = 100 +1.

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

Таблица 3. Результаты второй серии вычислительных экспериментов

Число шагов 5 = бт = 10 бт =

0 б,

10 20 0,1338 0,1424 0,0809

20 40 0,0608 0,0682 0,0183

50 100 0,0307 0,0215 0,0075

100 200 0,0274 0,0096 0,0037

8

к

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

Очевидно, для получения удовлетворительного результата (сочтём за него решение с не более чем 2% погрешностью) достаточно использовать цуг в 10 длин волн. При QT = 15 удовлетворительный результат достигается при дискретизации Q, Qt равной (20,40). При меньшей длине цуга поле в исследуемой области не успевает устояться, волна не является монохроматической в силу чего, результаты характеризуются высокой погрешностью. Таким образом, можно заключить, что разработанный метод является достойной альтернативой "чистому" БЭТЭ.

ЗАКЛЮЧЕНИЕ

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

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

При совместном численном решении была использована методика TF/SF, применение которой при использовании только волнового уравнения является малоэффективной. Результаты использования методики при совместном решении уравнений Максвелла и волнового уравнения практически неотличимы от аналогичного решения уравнений Максвелла, что позволяет говорить об эквивалентности сеточных задач и, соответственно, получаемых решениях.

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

СПИСОК ЛИТЕРАТУРЫ

1. Дифракционная нанофотоника [под ред. В.А. Сой-фера]. М.: Физматлит, 2011, 680 с.

2. Taflove A., Hagness S. Computational Electrodynamics: The Finite-Difference Time-Domain Method: 3ed. ed. Boston: Arthech House Publishers, 2005. P. 852.

3. Parallel finite-difference time-domain method / Wenhua Yu, Raj Mittra, Tao Su, Yongjun Liu, Xiaoling Yang // ARTECH HOUSE, INC - Boston, 2006. P. 262.

4. Elsherbeni A. The Finite Difference Time Domain Method for Electromagnetics: With MATLAB Simulations. Atef Elsherbeni, Demir Veysel - SciTech Publishing, 2009. 450 р.

5. Головашкин Д.Л., Яблокова Л.В. Совместное разностное решение уравнений Даламбера и Максвелла. Одномерный случай // Компьютерная оптика. 2012. Т.36. №4. С. 527 - 533.

6. Ваганов Р.Б. Каценеленбаум Б.З. Основы теории дифракции. М. : Наука, 1982. 272 c.

TASK OF FALLING WAVE ON TECHNOLOGY OF TF/SF AT THE FINITE-DIFFERENCE SOLUTION OF THE D'ALEMBERT AND MAXWELL'S EQUATIONS

© 2014 L.V. Yablokova1, E.Y. Buldygin1, D.L. Golovashkin2

1 Samara State Aerospace University named after Academician S.P. Korolyov (National Research University) 2 Image Processing Systems Institute of the RAS, Samara

The technique of searching of the collateral finite - difference solution of a wave equation and set of equations of Maxwell is offered, allowing to combine advantage and to avoid shortcomings of both made mention numerical methods of nanophotonics. In a two - dimensional case on test examples convergence of such decision, possibility a task of an incident wave on the TF/SF technology are shown. Key words: D'Alembert equation, Maxwell's equations, difference scheme, TF/SF technique.

Lyudmila Yablokova, Senior Lecturer at the Applied Mathematics Department. E-mail: [email protected] Evgenii Buldygin, Student of MA course. E-mail: [email protected]

Dimitry Golovashkin, Doctor of Physical and Mathematical Sciences, Associate Professor, Leader Researcher. E-mail: [email protected]