Совместное разностное решение уравнений Даламбера и Максвелла. Двумерный случай Текст научной статьи по специальности «Математика»

СОВМЕСТНОЕ РАЗНОСТНОЕ РЕШЕНИЕ УРАВНЕНИИ ДАЛАМБЕРА И МАКСВЕЛЛА.

ДВУМЕРНЫЙ СЛУЧАИ

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

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

Аннотация

Предложена методика отыскания совместного разностного решения волнового уравнения и системы уравнений Максвелла, позволяющая использовать достоинства и избежать недостатков обоих упомянутых численных методов нанофотоники (сократить затраты памяти и применить известные методики задания падающей волны и наложения поглощающих слоёв). В двумерном случае на тестовых примерах продемонстрированы сходимость такого решения, возможность наложения РЫЬ-слоёв и задания падающей волны по технологии ТГ/БГ.

Ключевые слова: уравнение Даламбера, уравнения Максвелла, разностные схемы, РЫЬ-слой, методика ТГ/БГ.

Введение

Развитие нанофотоники [1] сопровождается как появлением новой элементной базы, так и разработкой численных методов расчёта дифракции света. Наибольшей популярностью среди последних пользуется FDTD (Finite - Difference Time - Domain) подход [2], основанный на разностном решении уравнений Максвелла и характеризующийся в силу этого универсальностью. Он позволяет моделировать распространение произвольных электромагнитных волн в любой среде, в том числе и с учётом их нелинейного взаимодействия. К сожалению, применение упомянутого метода сопровождают высокая вычислительная сложность и требования к объёму оперативной памяти, являясь платой за перечисленные достоинства.

Сокращение вычислительной сложности при реализации FDTD-метода традиционно связывают с использованием аппаратных средств, позволяющих организовывать параллельные и векторные вычисления по нему. Результаты многолетних исследований в этой области подробно представлены в монографиях [2-4] и использованы при разработке множества программных пакетов, среди которых хочется выделить свободно распространяемые с открытым кодом MEEP [5] (кластерные вычисления) и B-CALM [6] (вычисления на графических процессорах).

Проблеме сокращения объёма оперативной памяти при организации вычислений по рассматриваемому методу также посвящён ряд работ. Например, в [7] и [8] для этого используются подвижные и неравномерные сеточные области, в [9] предлагается декомпозиция области, в [10] - метод пирамид.

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

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

фракции ТЕ-волны [2] на цилиндрическом элементе (упрощение модели), приводящему к уравнению Да-ламбера. Разностное решение последнего (как будет показано далее) совпадает с результатами расчётов по FDTD-методу в указанном случае, но отличается хранением одной сеточной функции на двух временных слоях в отличие от FDTD (три сеточных функции на одном слое по времени), что позволяет на треть сократить требования на объём используемой оперативной памяти. Простота программной реализации также выгодно отличает выбранный подход.

Будучи в свое время в фокусе внимания исследователей, разностное решение уравнения Даламбера постепенно утратило популярность в связи с развитием FDTD-метода и применялось исключительно для формирования поглощающих граничных условий [11] для последнего. Возвращение интереса к нему в настоящее время [12] связано с отмеченными преимуществами, но происходит недостаточно быстро в силу отсутствия разработанного математического аппарата задания падающего поля, наложения поглощающих слоёв и учёта дисперсии материала. Последняя проблема недавно успешно решена в [13], а преодолению двух первых посвящена предлагаемая работа. Авторы находят разумным совместить разностные решения уравнения Даламбера (в области нахождения оптического элемента) и уравнений Максвелла (в поглощающих слоях), сочетая достоинства первого (сокращение требования к объёму оперативной памяти) и второго (разработанный математический аппарат ТГ/БГ-методики и РЫЬ-слоёв) методов.

1. Согласование разностных решений уравнений

Даламбера и Максвелла на границе сеточных подобластей

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

ЭЯ дЕ„ ЭН ЭЕ„

m

0 dt eoe"E"

dz

, mo

dt dy

dH dHr

(1)

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

ут,к+0,5 ут,к+0,5

Еп _ Еп

хт,к+1 хт,к

Нп+0,5 нп_0,5 еп еп

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

хт,к хт,к

Еп+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,5

ственно. Решая уравнение (1), введём Б (0<(<Т, 0<у<Ьу, 0<г<Ь), на которую наложим сеточную область Бь. В этом случае проекцию ЕП к определим в

узлах {(П, уш, г,): П=пЪь п=0, 1, .., Ы=Т/Ьь Ут=тЬу, т=0, .., М=Ьу/Ьу, гк=кЬг, к=0, .., К=Ь/Ьг}, проекцию - {(П+0'5' Ут" 2к+0'5): П+0,5=(п+0,5)ьь п=0, 1, .., N-1, ут=тЬу, т=1, .., М-1, гк+0,5=(к+0,5)Ьг, к=0, .., К-1}, проекцию Н^к - {(П+0,5, Ут+0,5, г)

П+О,5=(п+0,5)ьь п=0, 1, .., N-1, Ут+0,5=(т+0,5)Ьу, т=0, .., М-1, гк=кИг, к=1, .., К-1}. Выделим на Б подобласть № «<Т, Ь<у<Я, В<г<и), задав на ней Ц. Будем считать, что ЕП определена в узлах{(П ут, г)

хт,к

П=пЬь п=0, 1, .., Ы=Т/Ьь ут=тЬУ, т=Ь,...,Я, Ь=Ь/ЬУ, Я=Ь/Ъу, гк=Ш, к=В, .„и, и=Ь/Ьг, В=Ь/Ь}, Н^-

{((п+о,5, Ут, гк+о,5): П+о,5=(п+0,5)Ьь п=0, 1, .., N-1, Ут=тЬу, т=Ь+1, .., Я-1, гк+0,5=(к+0,5)Ьг, к=В, .., и-1},

НП^ - {(П+05, Ут+о,5, г) П+о,5=(п+0,5)Ьь п=0, 1, .., ^

1, Ут+о,5=(т+0,5)ЬУ, т=Ь, .., Я-1, гк=кЬг, к=В+1, .., К-1}. Тогда решение (1) по схеме (2) будет отыскиваться на Ц1 = Б\Ц, а сеточные уравнения схемы (2) решаться на БМ = БЬ\Ц.

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

э и

эи + э2и„

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

= 0

ит+1 _ 2и:л+ит_1 = _с_ и+и _ 2ит л+ит

Ь2 е_,.( Ы

(3)

ип _ 2ип + ип

^ т,к+1 т,к Т ^ т,к_1 )

+ ,2 ).

(4)

Решение (3) будем искать в Ц сеточной области Б^, наложенной на Ц, реализуя вычисления по (4)

на сеточной области Б^. Тогда проекция ит к определена в узлах {(П, ут, гк): П=пЬь п=0, 1, .., ^Т/Иь ут=тЬУ, т=Ь, .., Я, Ь=Ь/ЬУ, Я=Ь/ЬУ, гк=кЬг, к=В, ..,и, В= ЬЬ и=Ь/Ь2}.

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

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

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

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

правая граница Б^ т=Я, В<к<и

верхняя граница Б^ к=и, Ь<т<Я

нижняя граница Бь к=В, Ь<т<Я

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

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

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

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

т=0, 0<к<К левая граница БМ ЕП = 0 хт,к

к=0, 0<т<М нижняя граница БМ ЕП = 0 хт,к

к=К, 0<т<М верхняя граница Б^ ЕП = 0 хт,К

т=М, 0<к<К правая граница Б>М ЕП = 0 хМ ,к

0<т<Ь, Я<т<М, 0<к<К область БМ (2)

Ь<т<Я, 0<к<В, и<к<К область БМ (2)

Ь<т<Я, В<к<и область Б^ (4)

т=Ь,В<к<и левая граница Б™ (5)

т=Я,В<к<и правая граница Б^ (6)

к=В, Ь<т<Я нижняя граница Бь (7)

к=и, Ь<т<Я верхняя граница Б™ (8)

т=Ь, к=В нижний левый узел б: (9)

т=Ь, к=и верхний левый узел б: (10)

т=Я, к=и верхний правый узел б: (11)

т=Я, к=В нижний правый узел б: (12)

где формулы (5)-(12) имеют вид:

Ь

Ь

г

Ь

Ь

е0ет.к

Ь

Ь

+

U - 2Un + U

c Un - 2Un + En

c ,m+1,k ^^m,k ^ Em-1,k

^ m i-

Un - 2Un + Un

^ m,k+1 m,k m,k-1 )

+ ,2

U - 2Un + U

m,k m,k m,k

c Em+i,k- 2Um,k+Um~i,i

em„(

Un - 2Un + Un

^ m,k+1 m,k m,k-1 )

+ 72

U - 2Un + U

m,k m,k m,k

c Mm+i,k- 2Um,k+Um-i,k

^ m i-

J Tn _ rs T Tn . Tjn

+ Um,k+1 2Um,k Em,k-1 ) + ,2

Un+1 - 2Un + Un-1

m,k m,k m,k

c Um+1,k- 2Um,k+Um-u

e„t(

En _ 2Un + Tin + Em,k+1 2Um,k U m,k-1 )

Un+1 - 2Un + Un-1

m,k m,k m,k

c ,Um+1,k-2Um,k+Em-1,1

JTn - 2 U jn + En

+ U m,k+1 m,k Em,k-1 )

i 2

Un+1 - 2Un + Un-1

m,k m,k m,k

c Un - 2Un + En

c , m+1,k ^^m,k ^ Em-1,k

em„(

pn - 2 U fn + rrn

+ Em,k+1 2Um,k U m,k-1 )

U - 2Un + U

m,k m,k m,k

c ,K+1,k - 2Um,k + Um-1,l

em„(

pn _ llfn+TTn

+ Em,k+1 2Um,k Um,k-1 )

T2

Un+1 - 2Un + U

m,k m,k m,k

c En - 2Un + Un

^ m i-

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

I jn —Ol jn ~pn

+ UmA+1 2Um,k Em,k-1 )

hZ *

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

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

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

Lz Я Я и 1 хххз ■ Д ■ А ■ А 1 X X X J 1 ■ ■ ■ - - - ■ ■ 1 е х х х • • • х х : IA НА ■ А ■■■■■А ■ А 1 е х х х ■ • ■ х х : 1 в в : х х 1 А В А ■ : х х

х х х : ■ А ■ А ■ А 1 X X X 3 : х х х ■ • ■ х х : iahahah -iahai е х х х • ■ ■ х х : : X X 1 А ■ А В X X

Lu XXX ■ А ■ А ■ А [ XXX ■ А ■ Д ■ А [ XXX 1 □ □ □ ■ ■ ■ □ □ 1 ■ □ □ □ - - - □ □ 1 X X 1 А ■ А ■ X X 1 А ■ А В X X

XXX ■ А ■ А ■ А [ XXX ■ Д ■ А ■ А I XXX 1 □ □ □ ■ ■ ■ □ □ [ Dw 1 □ □ □ ■ ■ ■ □ □ 1 X X 1 А ■ А В X X ■ А ■ А В X X

Ьв * * * ^ ■ А ■ А ■ А 1 X X X 3 : х х х ■ ■ ■ х х з ia U Л Я А Я ■■■ Я А Я ai е х х х ■ • ■ х х : X X 1 А В А ■ : х х

X X X 3 ■ А ■ Д ■ А 1 X X X 3 ■ ■ ■ 1 Ll\ : х х х ■ • • х х : 1АВДВДВ---ВДВД1 ; х х х ■ ■ • х х з 1 В В ■ ■ ■ ■ В ■ 1 LR\ х * DM 1 А В А В X X 1 в в Ly

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

гп+0,5

Яи +0,5 ттп+0,\

z t , крестикам - H у

ym,i+0,5

Эксперименты проводились при различных значениях дискретизации сеточной области ((,(>1 и (Т, где первый параметр характеризовал число узлов сеточной области по пространству (приходящееся на одну длину волны); второй - количество узлов по времени (приходящееся на временной интервал, за который плоский волновой фронт в вакууме пройдёт расстояние в одну длину волны); третий - «длительность» запускаемого цуга в длинах волн. При этом они менялись от (10, 20, 5) минимальных значений, удовлетворительно описывающих распространение плоской однородной волны в свободном пространстве, до (100, 200, 15), соответствующих весьма низким величинам погрешности. В качестве вычислительной области Б брался квадрат с длиной стороны 201, при этом шаги по пространству полагались равными.

В первой серии экспериментов были повторены результаты из работы [14], полученные при распространении плоской однородной электромагнитной волны в вакууме при задании жёсткого источника. Их совпадение свидетельствует о правомерности излагаемого выше подхода в случае свободного пространства.

2. Задание падающей волны по технологии ТЕ/БЕ при совместном решении

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

+

2

2

h

h

+

2

2

h

h

+

2

2

h

h

+

2

2

h

h

+

2

2

h

h

+

2

2

h

h

+

2

2

h

h

z

+

2

2

h

h

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

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

E - Б"

xL—1,k xL—1, k

Hn+0,5 — (Нл+0'5

zL—0,5,k zL— i,5,k

H n+0,5

+ H ZL—1,5,k )

h

Hn+0,5 — Hn+0,5

yL— i,k+0,5 yL—1,k—0,5

En+1 — En

xR+1,k xR+1,k

(ял+0'5 + HyR+15k) — ял+0'5

V ZR+i,5,k yR+1,5,^ ZR+0,5,k

H"n+0,5 — H"n+0,5

yR+1,k+0,5 yR+1,k—0,5

E — E

Hn+°'5 — H'

zm+0,5,U+1

n+ 0,5

m—0,5, U+1

(Hn

hh n+0,5

+ HУт,и+1,5 ) — H-

n+0,5

ym,U+0,5

En+1 — En

xm,B—1 xm,B—1

Hn+0,5 — Hn+0,5

zm+0,5,B—1 zm—0,5,B—1

Hn+0,5 _(Hn+0,5

ym,B—0,5 ym,B—0,5

H n+0,5

+ Hy»,B—1,5 )

m

m

m

m

Hn+0,5 — Hn—0,5

zL—1,5,k zL—1,5,k

Hn+0,5 — Hn—0,5

zR+1,5,k zR+1,5,k

Hn+0,5 — Hn—0,5

ym,B—1,5 ym,B—1,5

Hn+0,5 — Hn—0,5

ym,U+1,5 ym,U+1,5

(En xL-1,k - EHxL-1,k ) - Exn xL— 2,k .

Exn xR+2,k (En V xR+1,k H xR+1,k )

EHnx ) En 1 xm,B 2

Exm,u+2 (exlHx»,U+1)

(13)

А 4

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

Hn+0,5 — ТГ0,5 E" — En

H yk+05 Hyk+0,5 Exk+1 Ext

h

k = 0, K —1, n = 0, N — 1,

e-

Hn+1 нп Hn+0'5 _ Hn+0,5

Üxk — Üxk _ Hyk+0,5 H yk—0,5

h

(14)

k = 0, K, k = 0, N —1, E b—1 = sin wnA,.

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

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

|л — nj

(eR = maxJ-1, Nk- значения, полученные разно-

k Л

стными методами, Ak - соответствующие аналитическому решению), отыскивая её на оптической оси элемента, варьируя k = B,U и m = 10Q + 1.

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

Число шагов QT = 5 QT = 10 QT = 15

Q Qt eR

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

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

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

h

h

h

h

h

h

h

e

h

h

h

h

h

h

h

h

h

h

3. Наложение поглощающего слоя при совместном решении

Важной задачей при моделировании дифракции электромагнитного излучения разностным методом является помещение поглощающих слоёв (PML - perfectly matched layer [16]) у границ вычислительной области для исключения отражения от последних. По замыслу авторов настоящей работы, рассмотрение задачи дифракции на оптическом элементе уместно связывать с решением волнового уравнения в области непосредственного нахождения этого элемента, а организацию поглощения у границ вычислительной области (вне элемента) - с решением уравнений Максвелла, моделирующим PML-слой. При этом сокращаются требования к системным ресурсам ЭВМ по сравнению с FDTD-методом (в двумерном случае) и используется современная методика наложения PML-слоя вместо устаревшего подхода Мура [11], традиционно сопровождающего решение волнового уравнения.

Один из способов записи уравнений Максвелла в поглощающем слое (предложен Беренже в 1994 году [16]) характеризуется введением, наряду с электрической проводимостью о, магнитной о*, связанных соотношением о = (e0 / m0)s*. Тогда уравнения (1) принимают вид:

ЭЯ, _ dE„

m-

m

e

dt ЭЯ dt

E dt

7 + о* Hy = -

+ о* H = -

dz

Э7

(15)

+ оE„ =

H H7

Эу Эг

Представляя распространение и затухание поля вдоль разных направлений отдельными уравнениями, производят расщепление электрической компоненты и проводимостей [16], записывая

m

m

ЭН, + о*Я =-Э( E" - ) ;

dt z 7 dz

ЭН d(Exz + E )

e-

t

E dt E

+оН =

dy

dH

+ 0 E =--7;

z ^ dz

H

(16)

t

где Ex = E}

+ 0 E =

7 xy

dy

+ Ey , а отличные от нуля проекции р(рУ

= s0о* и fi(pz = s0о* обеспечат затухание волны в соответствующих направлениях.

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

(

о = о

y max

о = о

y max

Ly - У

L

в областях 1, 2, 3;

У J

К - К+у

, в областях 5, 6, 7;

f L - z

V h J

в областях 1, 8, 7;

f L-1 + z V

z z

v z

в областях 3, 4, 5, где qe R.

h Iy-LY

Ly 0

7 6 5

8 4

1 2 3

О Ьг 1г

Рис. 2. Схема расположения РМЬ

Соответствующая конечно-разностная аппроксимация (16) в поглощающих слоях [16]:

гл+0,5 ггл-0,5

mo ну

- Яn

ym,k+0,5 ym,k+0,5 . _* ГГЛ-0,5

--+о7 НУ

и zm,k+0,5 ym,k+0,5

ht

E" - E" + E" - E"

xym,k+1 xy.m,k xzm,k.+l xzm,k .

m° Hn;°'5 - Hn-0-5

r 0 z

m+0,5,k . _* T rn-0,5

--+о y H

ym+0,5,k zm+0,5,k

E" - E" + E" - E"

xym+1,k xym,k xzm+1,k xzm,k

(17)

En+1 - En

zm,k

h

- + о E" = -

zm,k xzm,k

Нл+0,5 - нл+0,5

ym,k+0,5 ym,k-0,5

e

E"+1 - E"

xym,k xym,k . TJU _

m,k h о 7m,k Xym,k =

Нл+0,5 - Нл+0,5

zm+0,5,k zm-0,5,k

h„

Упрощая программную реализацию вычислений по схеме (17) и оптимизируя её векторизацию (увеличивая длину векторов в РЫЬ-области [17]), авторы переходят от краевых условий первого рода (электрической стенки) к периодическим условиям [18], что позволяет объединить поглощающие слои. В частности, прежние слои с нечётными номерами объединятся в один, также попарно объединятся слои 2, 6 и 4, 8 (рис. 3).

Таким образом, на рис. 3 будем различать границы вычислительной области (чёрная линия), поглощающих слоёв (тонкая чёрная линия), рассеянного и результирующего полей (пунктирная линия) и сеточных областей ОМ и Щ (линия из точек).

В ходе третьей серии вычислительных экспериментов определим толщину объединённых поглощающих слоёв в две длины волны, длину квадратной

о = о

z max

о = о

z max

h

h

h

e

m,k

h

L

вычислительной области (ЬУ=Ь) положим равной четырём длинам волн (сократив по сравнению с предыдущими экспериментами на порядок за счёт использования РМЬ).

О 12-2Ь2 12-Ь2 1г

О

lr-2Lr

Iy-Ly lr

Волновое уравнение

Уравнения максвелла

Результирут .....ющее поле

Рассеянное поле

Рис. 3. Взаимное расположение ТР/БР, РМЬ и области решения волнового уравнения при циклических граничных условиях

Параметры поглощающих слоёв ашах и ц будем варьировать в зависимости от дискретизации сеточной области, следуя рекомендациям из [19]. Границу раздела полей расположим непосредственно за граничными условиями (левая и верхняя части границы раздела) и РМЬ-слоями (нижняя и правая части) на расстоянии одного шага сеточной области по пространству от них. На таком же расстоянии от границы полей расположим границу раздела сеточных областей. Центр диэлектрического цилиндра разместим в точке

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

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

Число шагов QT =10 QT = 15

Q Qt £r

10 20 0,1629 0,1576

20 40 0,0628 0,0592

50 100 0,0206 0,0221

100 200 0,0087 0,0088

Отметим, что в отличие от предыдущего случая (табл. 3) увеличение длительности временного интервала с (Т=10 до (Т=15 не привело к заметному снижению погрешности (табл. 4). Если при (Т=10 результаты из табл. 3 и 4 почти совпадают (применение РМЬ не внесло значимых искажений в результат), то для (Т=15 они отличаются в два-три раза. Следовательно, величины из последней колонки табл. 4 характеризуют в значительной степени погрешность, связанную с использованием РМЬ-слоёв.

Заключение

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

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

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

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

Литература

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

Taflove, A. Computational Electrodynamics: The Finite-Difference Time-Domain Method / A. Taflove, S. Hagness. ed. 3-d. - Boston: Arthech House Publishers, 2005. - 852 p. Yu, W. Parallel finite-difference time-domain method / Wenhua Yu, Raj Mittra, Tao Su, Yongjun Liu, Xiaoling Yang - Artech house, Inc - Boston, 2006. - 262 p. Elsherbeni, A. The Finite Difference Time Domain Method for Electromagnetics: With MATLAB Simulations / Atef El-sherbeni, Demir Veysel. - SciTech Publishing, 2009. - 450 р. http: // ab-initio.mit.edu / wiki / index.php / Meep. http: // b-calm.sourceforge.net .

Fidel, B. Hybrid ray - FDTD moving window approach to pulse propagation / B. Fidel, E. Heyman, R. Kastner, R.W. Zioklowski // Journal of Computational Physics. - 1997. -Vol. 138, Issue 2. - P. 480 - 500.

Shi, S. Analysis of diffractive optical elements using a nonuniform finite - difference time - domain method / S. Shi, X. Tao, L. Yang, D.W. Prather // Optical Engineering. -2001. - Vol. 40, № 4. - P. 503 - 510. Головашкин, Д.Л. Декомпозиция сеточной области при разностном решении уравнений Максвелла / Д.Л. Головашкин, Н.Л. Казанский // Математическое моделирование. - 2007. - Т. 19, № 2. - C. 48 - 58.

10. Головашкин, Д.Л. Решение сеточных уравнений на графических вычислительных устройствах. Метод пирамид / Д. Л. Головашкин, А. В. Кочуров // Вычислительные технологии. - 2012. - Т. 17, N 3. - C. 55 - 69.

11. Mur, G. Absorbing boundary conditions for the finite -difference approximation of the time - domain electromagnetic field equations / G. Mur / IEEE Trans. Electromagnetic Compability. - 1981. - Vol. 23. - P. 377 - 382.

2.

3

4

8

9

12. Козлова, Е.С. Моделирование распространения короткого двумерного импульса света / Е.С. Козлова, В.В. Котляр // Компьютерная оптика. - 2012. - Т. 36, № 2. - С. 158-164.

13. Козлова, Е.С. Моделирование предвестников Зоммер-фельда и Бриллюэна в среде с частотной дисперсией на основе разностного решения волнового уравнения / Е.С. Козлова, В.В. Котляр // Компьютерная оптика. -2013. - Т. 37, №2. - C. 146-154.

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

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

16. Berenger, J.-P. A perfectly matched layer for the absorption of electromagnetic waves / Jean-Pierre Berenger // Journal of Computational Physics. - 1994, № 114. - P.185-200.

17. Ортега, Д.М. Введение в параллельные и векторные методы решения линейных систем / Д.М. Ортега; пер. с англ. Х.Д. Икрамова, И.Е. Капорина; под ред. Х.Д. Ик-рамова. - М.: Мир, 1991 - 364 с.

18. Головашкин, Д.Л. Постановка излучающего условия при моделировании работы цилиндрических дифракционных оптических элементов методом разностного решения уравнений Максвелла / Д.Л. Головашкин // Математическое моделирование. - 2007. - Т. 19, № 3. - C. 3-14.

19. Головашкин, Д.Л. Методика формирования падающей волны при разностном решении уравнений Максвелла (двумерный случай) / Д.Л. Головашкин, Н.Л. Казанский // Автометрия. - 2007. - Т. 43, № 6. - С. 78-88.

References

1. Diffraction nanophotonics / edited by V.A. Soifer - Moscow: "Fizmatlit" Publisher, 2011. - 680 p. - (In Russian).

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

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

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

5. http: // ab-initio.mit.edu / wiki / index.php / Meep.

6. http: // b-calm.sourceforge.net .

7. Fidel, B. Hybrid ray - FDTD moving window approach to pulse propagation / B. Fidel, E. Heyman, R. Kastner, R.W. Zioklowski // Journal of Computational Physics. - 1997. -Vol. 138, Issue 2. - P. 480-500.

8. Shi, S. Analysis of diffractive optical elements using a nonuniform finite - difference time - domain method / S. Shi, X. Tao, L. Yang, D.W. Prather // Optical Engineering. -2001. - Vol. 40, № 4. - P. 503-510.

9. Golovashkin, D.L. Mesh Domain Decomposition in the Finite-Difference Solution of Maxwell's Equations / D.L. Golovashkin, N. L. Kazanskiy // Optical Memory & Neural Networks (Information Optics). - 2009. - Vol. 18, № 3P. P. 203-211.

10. Golovashkin, D.L. Solution of difference equations in graphical computing devices. Method of pyramids. / D.L. Golovashkin, A.V. Kochurov // Computing Technologies. - 2012. - Vol. 17, N 3. - P.55-69. - (In Russian).

11. Mur, G. Absorbing boundary conditions for the finite -difference approximation of the time - domain electromagnetic field equations / G. Mur / IEEE Trans. Electromagnetic Compability. - 1981. - Vol. 23. - P. 377-382.

12. Kozlova, E.S. Model operation of distribution of a short twodimensional impulse of light / E.S. Kozlova, V.V. Kotlyar // Computer Optics. - 2012. - Vol. 36, N 2. - P. 158-64. - (In Russian).

13. Kozlova, E.S. Simmulations of Sommerfeld and Brillouin pre-cursorsin the medium with frequency dispersion using numerical method of solving wave equations / E.S. Kozlova, V.V. Kotlyar // Computer Optics. - 2013. - Vol. 37, N 2. - P. 146-154. - (In Russian).

14. Golovashkin, D.L. Joint finite-difference solution of the Dalamber and Maxwell's equations. One - dimensional case / D.L. Golovashkin, L.V. Yablokova // Computer Optics. - 2012. - Vol. 36, N 4. - P. 527-533. - (In Russian).

15. Vaganov, R.B. Fundamentals of the theory of diffraction / R.B. Vaganov, B.Z. Kacenelenbaum - Moskow: "Nauka" Publisher, 1982 - 272 p. - (in Russian).

16. Berenger, J.-P. A perfectly matched layer for the absorption of electromagnetic waves / Jean-Pierre Berenger // Journal of Computational Physics. - 1994. - № 114. - P. 185-200.

17. Ortega, M. Introduction to Parallel and Vector Solution of Linear Systems / M. Ortega. - Plenum Press, New York, USA, 1988. - 364 p.

18. Golovashkin D.L. Statement of radiating conditions at work modeling of cylindrical diffractive optical elements using differential solutions of Maxwell's equations / D.L. Golovashkin // Mathematic Modeling. - 2007. -Vol 19, N 3. - P. 3-14. - (In Russian).

19. Golovashkin, D.L. Incident wave source conditions for the finite-difference time-domain method: Two-dimensional formulation / D.L. Golovashkin, N.L. Kazanskiy // Optoelectronics, Instrumentation and Data Processing. - 2007. - V. 43(6). - P. 547-555.

JOINT FINITE - DIFFERENCE SOLUTION OF THE D'ALEMBERT AND MAXWELL'S EQUATIONS.

TWO - DIMENSIONAL CASE

E. Yu. Buldygin 2, D.L. Golovashkin 1, L. V. Yablokova 2

1 Image Processing Systems Institute of the Russian Academy of Sciences,

2 S.P. Korolyov Samara State Aerospace University (National Research University)

Abstract

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 (to reduce expenses of memory and to apply well-known TF/SF and PML techniques). In a two - dimensional case on test examples convergence of such decision, possibility of imposing of PML-ayers and a task of an incident wave on the TF/SF technology are shown.

Key words: D'Alembert equation, Maxwell's equations, difference scheme, PML-layer, TF/SF technique.

Сведения об авторах

Головашкин Димитрий Львович, доктор физико-математических наук, доцент, ведущий научный сотрудник Института систем обработки изображений РАН. Область научных интересов: разностное решение уравнений Максвелла (FDTD-метод), дифракционная оптика, векторные и параллельные матричные вычисления.

E-mail: [email protected] .

Dimitry Lvovich Golovashkin, Doctor of Physical and Mathematical Sciences, Associate Professor, Leader Researcher of Image Processing Systems Institute of the Russian Academy of Sciences. Scientific interests: FDTD-method, sub wave optics, vector and parallel algorithms of matrix computation.

Булдыгин Евгений Юрьевич, 1990 года рождения. В 2012 году получил степень бакалавра по направлению «Прикладная математика и информатика». В настоящее время является студентом магистратуры факультета информатики Самарского государственного аэрокосмического университета. Область научных интересов: разностное решение уравнений Максвелла.

E-mail: mathcuprum@gmail. com .

Evgenii Yurievich Buldygin, (b. 1990). Graduated (2012) with a bachelor's degree in Applied Mathematics and Informatics. At present he is a student of MA course at S.P. Korolyov Samara State Aerospace University. The main area of interests is finite difference methods.

Яблокова Людмила Вениаминовна, 1978 года рождения. В 2000 году окончила Самарский государственный педагогический университет. Работает в СГАУ (Самарский государственный аэрокосмический университет имени академика С.П. Королёва) старшим преподавателем кафедры прикладной математики. Область научных интересов: разностные схемы.

E-mail: lyablokova @yandex.ru .

Lyudmila Veniaminovna Yablokova, (b. 1978) graduated from Samara State Pedagogical University in 2000. She is senior teacher of Applied Mathematics sub-department of S.P. Korolyov Samara State Aerospace University. Research interests: finite-difference approximations.

Поступила в редакцию 29 ноября 2013 г.