Совместное разностное решение уравнений Даламбера и Максвелла. Одномерный случай

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

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

По сравнению с разностным решением уравнений Максвелла, аналогичное решение уравнения Даламбера позволяет в двумерном случае на треть сократить объём выделяемой оперативной памяти (сохраняются два временных сечения электрического поля вместо трёх проекций электромагнитного) и на 10% снизить число арифметических операций (10 при вычислении поля в одном узле сеточной области вместо 11 при вычислении трёх сеточных функций в ячейке Yee). Впрочем, при рассмотрении трёхмерной задачи метод FDTD становится более предпочтительным, что ограничивает область применения развиваемого подхода двумерным случаем.

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

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

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

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

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

Для одномерных уравнений Максвелла

дHy = dEx ;   SEx д t     Sz ’° St

^д H y дz

запишем разностную схему Yee [6]

Ц о

£ 0

H n + 0,5 y k + 0,5

E nn ^{+ 1} x_k

H

ht

-

ht

n - 0,5

y k + 0,5

En xk+1

- Enn xk

hz

Eⁿ x_k

-

_H n + 0,5 y k + 0,5

-

_H n + 0,5 yk - 0,5

hz

где сеточная проекция электрического поля на ось X – E_xⁿ _k определена в узлах

{( t n , z k ): t n = nh t , n =0, 1, .., N = T / h t , z k = L^W + kh z , k =0,.., K = L^M / h z },

а сеточная проекция магнитного поля на ось Y –

H +0,5 в узлах yk+0,5

{(tn+0,5, zk+0,5): tn+0,5 = (n +0,5)ht, n =0, 1, .., N–1, zk+0,5= LW +(k +0,5)hz, k =0,.., K–1} сеточной области DhM , наложенной (рис. 1) на область вычислительного эксперимента DM (0 < t < T, LW

Аналогично для одномерного волнового уравнения д2Ех   2 д2Е

---xL - c‘_x д t2дx

= 0

запишем разностную схему [7]

Enn+1 - 2 Enn + Enn-1       R"x   - 2 Enn + Enn xmxm_______xm __   2   xm+1_________xm ht2                                 hz2

где функция E_xⁿ_m определена в узлах

{(tn, zm): tn = nht, n =0, 1, .., N = T/ht, zm = mhz, m =0, .., M = LW/hz} сеточной области DhM , наложенной (рис. 1) на DW

(0<t<T, 0<z<LW).

тлм n^       f"Uh         ^

----------uh--------xO 0,5 1 1,5 2 ... K-0,5 K

I--------------U^--------- | >

^{IV L} aD^{M L}

Рис. 1. Объединение сеточных областей.

Квадратами изображены узлы для проекции электрического поля, окружностями – для магнитного

Задавшись целью совместного отыскания разностного решения уравнений (1) и (3), согласуем вычисления на DW и DM, полагая значения En , опре-xM делённые на DhW при m = M (крайний правый узел), и Exn0 , определённые на DhM при k =0 (крайний ле- вый узел), равными (рис. 1). На границах объединённой области D = DWuDM установим электрическую стенку. За начальное условие примем отсутствие поля в момент времени t =0.

Алгоритм перехода на следующий временной слой при реализации совместного решения состоит в определении сеточных функций в узлах D_h^W и D_h^M в соответствии с табл. 1, где (5) имеет вид:

En+1 - 2 En + En-1    JEn xM        xM      xM   c 2   x1

^ht²

- 2 En + En xm xM-1

^hz²

.

В (5) Е^₁ определено на области DM, Е"_{м 1}- на D^W, остальные значения относятся к обеим областям.

Для подтверждения правомерности предложенного подхода к совместному решению на одной вычислительной области уравнений Максвелла и вол- нового авторы провели четыре серии вычислительных экспериментов. При этом использовались операционная система Windows XP Professional SP3, вычислительный комплекс Octave 2.9.15 и процессор AMD Opteron 246.

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

Узлы сеточных областей	Подобласти вычислительных областей	Решаемые уравнения
m = 0	z =0, левая граница D и DW	Е"+1= 0 x0
0 < m< M–1	0 < z< LW, DW без границ	(4)
m = M и k = 0	z = LW, общая граница DW и DM	(5)
0< k< K	LW< z< LW+LM, DM без границ	(2)
k = K	z = LW+ LM, правая граница D и DM	Е"+1= 0 xK

Эксперименты проводились при различных значениях дискретизации сеточной области Q, Qt и QT, где первый параметр характеризовал число узлов сеточной области по пространству (приходящееся на одну длину волны); второй – количество узлов по времени (приходящееся на временной интервал, за который плоский волновой фронт в вакууме пройдёт расстояние в одну длину волны); третий – «длительность» запускаемого цуга в длинах волн. При этом они менялись от (10, 20, 5) – минимальных значений, удовлетворительно описывающих распространение плоской однородной волны в свободном пространстве, до (100, 200, 15) – соответствующих весьма низким величинам погрешности. Протяжённость области D составляла 20λ, что позволяло не учитывать влияние отражённой от границ волны в области регистрации.

Первая серия вычислительных экспериментов связана с разностным решением уравнений Максвелла (L^W =0, D = D^M) и имела своей целью получить значения погрешностей, относящихся исключительно к замене производных разностными отношениями. Далее указанные погрешности использовались для сравнения с результатами остальных экспериментов как величины, значительное отклонение от которых нежелательно и свидетельствует о несостоятельности выбранной в настоящей работе стратегии.

Излучение вводилось в область посредством «жёсткого» источника [7], расположенного в центре. Равномерная погрешность б = max 11 -1Ak 11 х 100% характеризовала максимальное отклонение от аналитического решения (модуль комплексной амплитуды распространяющейся в свободном пространст- ве плоской однородной монохроматической волны равен 1) на отрезке [10,5λ, 11,5λ].

Таблица 2. Погрешности ε первой и второй серий экспериментов

Q, Qt	QT
	5	10	15
10, 20	0,7268	0,2434	0,1275
20, 40	0,2217	0,1080	0,0371
50, 100	0,0716	0,0234	0,0113
100, 200	0,0234	0,0083	0,0046

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

Во второй серии решалось волновое уравнение при L^M = 0 и D = D^W. Дискретизация сеточной области для электрической компоненты поля, расположение и тип источника, размеры вычислительной области и области регистрации погрешности совпадали с предыдущим случаем.

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

Третья серия относится к совместному решению уравнений (1) и (3). Совокупная протяжённость объединённых вычислительных областей D^W и D^M составила 20λ (по 10λ на каждую). Источник располагался в области решения волнового уравнения на расстоянии половины длины волны от границы раздела (рис. 2а), область регистрации [10λ, 11λ] помещалась в D^M. Таким образом, расстояние от источника до этой области не менялось по сравнению с предыдущими экспериментами.

Dw (---— ^   J DM

Рис. 2. Вычислительные области третьей (a) и четвёртой (б) серий экспериментов. Звёздочкой обозначен «жёсткий» источник, в обведённой пунктиром области регистрируется погрешность

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

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

Таблица 3. Погрешности ε третьей серии экспериментов

Q,Qt	QT
	5	10	15
10, 20	0,5533	0,1907	0,0926
20, 40	0,3531	0,0562	0,0313
50, 100	0,0687	0,0217	0,0112
100, 200	0,0235	0,0081	0,0047

Для определения значения напряжённости электрического поля EXM' на стыке DW и D^¹ по (5) используется центральная по пространству разность, в которую входят: общее значение EnM , E"M 1 из DW и EX из D’M . При этом Eⁿn_M 1 определяется из (4) по тому же шаблону, а E_xⁿ₁– уже через сеточные значения магнитного поля (2), взятые на другом расстоянии от узла k = 1 (хотя шаблон при этом тоже центральный).

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

При расположении источника в области D^M на расстоянии половины длины волны от границы раздела (четвёртая серия экспериментов, рис. 2б), а области регистрации [9λ, 10λ] в D^W результаты отличаются от представленных в таблицах 2 и 3 весьма незначительно – тоже на одну десятитысячную процента для густой сеточной области.

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

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

Представляя уравнения Максвелла в поглощающем слое вместо (1) записывают:

дH  .

Цо —r + _cHy д t

дEx   дE    г

-    , s0    + GE дz   0 д t       x

⁵Hy дz ’

процента от модуля амплитуды (теоретического) падающей волны (табл. 4).

где второе слагаемое в левой части определяет поверхностную плотность электрического и магнитного (воображаемого) токов. Следующее соотношение электрической и магнитной проводимостей g/e₀ = g /ц₀

и плавное возрастание значения g от начала слоя к концу (при правильно подобранных параметрах) обеспечивает поглощение излучения без отражения [9].

Простейшая разностная аппроксимация уравнений Максвелла в поглощающем слое имеет вид:

Цо

Н^т+0,5 yk+0,5

Hm-0,5 yk +0,5

ht

* rjm-0,5

^{+ G}k+0,5 Hy +0,₅

^^^^^^B

Em xk+1

^^^^^^B

hz

E^m x_k

Таблица 4. Погрешности δ*100 пятой серии вычислительных экспериментов

Q, Qt	QT
	5	10	15
10, 20	0,3546	0,1032	0,1463
20, 40	0,1713	0,1238	0,1261
50, 100	0,1127	0,1068	0,1068
100, 200	0,0129	0,0115	0,0115

^E0

Em⁺¹

xk

^^^^^^B

E^m x_k

ht

+g, Em

k xk

^^^^^^B

_H т+0,5 yk+0,5

^^^^^^B

_H т+0,5 yk-0,5

hz

,

в котором сеточные функции проводимостей (электрическая и магнитная) определены в тех же узлах D_h^M , что и соответствующие напряжённости.

При постановке вычислительных экспериментов расположим поглощающий слой длины λ у правого края D.

Тогда σ задаётся на отрезке [L^W + L^M – λ, L^W + L^M] по формуле:

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

При отыскании совместного разностного решения (1) и (3) (шестая серия экспериментов) примем L^W = 10λ, L^M =2λ, расположенный в том же месте вычислительной области «прозрачный» источник теперь находится в D^W (рис. 3б).

g = g

W M^q

z + X — L — L |

max

X

где G_max и q зависят от дискретизации сеточной области согласно табл. 4 из [11].

Верифицируя предположение о возможности совместного решения волнового уравнения и системы уравнений Максвелла с поглощающим слоем, сначала (для получения эталонного решения) рассмотрим случай D = D^M (пятая серия экспериментов, рис. 3а); L^W = 0, L^M = 12λ. Для исследования отражённой волны зададим «прозрачный» источник [2] (в точке z = 9,5λ) и будем регистрировать в подобласти [10Z, 11Z] величину 5 = max| АГ”j (Аотр'- комплексная амплитуда отражённой волны), являющуюся абсолютной погрешностью вычислительного эксперимента (отражённой от PML-слоя волны в идеальном случае быть не должно).

Таблица 5. Погрешности δ*100 шестой серии вычислительных экспериментов

Q, Qt	QT
	5	10	15
10, 20	0,1287	0,1312	0,1483
20, 40	0,1645	0,1232	0,1322
50, 100	0,1136	0,1068	0,1068
100, 200	0,0127	0,0115	0,0115

•                                                                                                                                  ।_______________II

Х/2 XX б>-----@-------( II

D1* DM

Рис. 3. Вычислительные области пятой (a) и шестой (б) серий экспериментов. Символом «@» обозначен «прозрачный» источник, в обведённой пунктиром области регистрируется погрешность, в закрашенной расположен поглощающий слой

Отметим удовлетворительное функционирование поглощающего слоя; даже на самой грубой сеточной области (Q = 10, Qt =20) при работе с не вполне монохроматической волной (QT = 5) модуль комплексной амплитуды отражённого поля составил треть

Сравнивая результаты вычислительных экспериментов (табл. 5) с данными пятой серии (табл. 4), отметим их совпадение в диапазоне Q = 50, Qt = 100, QT = 10– Q = 100, Qt = 200, QT = 15 и незначительность отличий при остальных параметрах, что свидетельствует о возможности совместного разностного решения уравнений Даламбера и Максвелла в случае, когда последние используются для моделирования поглощающего слоя.

3. Задание падающей волны по технологии TF/SF при совместном решении

Формирование падающей волны через задание «прозрачного» источника (как это было сделано в предыдущем параграфе) не всегда удобно. Излучение от такого источника распространяется не только в желаемом исследователю направлении, но и в противоположном (этим, например, характеризуется пакет MEEP [12], реализующий FDTD-метод). Кроме того, поглощению у границ области в этом случае подлежит результирующее поле, состоящее из падающего и рассеянного. Всё это зачастую приводит к увеличению погрешности вычислительного эксперимента.

Задание падающей волны по классической технологии TF/SF (total field/scattered field [13]) в рамках FDTD-метода позволяет избежать распространения в противоположном направлении и появления падаю-

щего поля в области расположения PML-слоёв (там окажется только рассеянная волна).

Основой технологии является разделение полей при решении (1) на результирующее (в подобласти kL ≤ k ≤ kR области D_h^M ) и рассеянное (k< kL и k > kR), при котором в соответствующих узлах вместо (2) записываются разностные уравнения

Цо

_H n+0,5

ykL-0,5

^^^^^^B

_H n-0,5 ykL-0,5

te

n

^ykL

^^^^^^B

Ek)

^^^^^^B

-

Таблица 6. Погрешности ε седьмой серии экспериментов при аналитическом задании падающей волны

Q, Qt	QT
	5	10	15
10, 20	4,0678	3,1662	2,2889
20, 40	1,2341	0,733	0,7927
50, 100	0,4902	0,2658	0,1752
100, 200	0,126	0,0899	0,0793

ht

hx

;

^0

_En+1 x_kL

^^^^^^B

Eⁿ x_kL

H

• n+0,5

ykL+0,5

^^^^^^B

(H

• n+0,5

ykL-0,5

+ Hn"И ykL-0,5

ht

hx

;

Ц

Hⁿ+0,5 ykR +0,5

^^^^^^B

Hn-0,5

ykR +0,5

Enn  - (E xkR +1

n

^xkR

^^^^^^B

E^{ɶn x}kR

ht

hx

;

^0

_En+1 ^xkR

^^^^^^B

En+1   ( h xkR  _\

• n+0,5

ykR+0,5

+Hin+0,5) ykR+0,5

^^^^^^B

H

• n+0,5

ykR-0,5

ht

hx

,

где сеточные функции под тильдой соответствуют аналитически [13] или численно [14] заданному падающему полю.

В работе [10] высказывается предположение о возможности использования приведённой методики при решении волнового уравнения. Однако авторам настоящей статьи не известно ни одной успешной реализации (включая собственную попытку [5]) этого замысла. Выходом из сложившейся ситуации может стать совместное решение уравнений Даламбера и Максвелла, при котором область D^W содержит исключительно общее поле, а к D^M относится граница раздела общего и рассеянного полей.

Предваряя экспериментальную проверку этого утверждения, проведём серию вычислительных экспериментов (седьмую) для определения погрешностей FDTD-метода при задании падающей волны по методике TF/SF. Примем для этого L^W =0, L^M =20λ, D = D^M, задав результирующее поле на отрезке [9λ,11λ] (рис. 4а).

DM

• a)------#-------1                  !-------#-----

• рассеянное результирующее рассеянное поле              ’ поле                поле

• б)------#------1               I------#-----

• D^ Dw D^

Рис. 4. Вычислительные области седьмой (a)

и восьмой (б) серий экспериментов. Символом «#» обозначена граница раздела результирующего и рассеянного полей, в обведённой пунктиром области регистрируется погрешность

Регистрация погрешности ε производится на отрезке [9,5λ, 10,5λ] в центре D. Поглощающие слои у границ области не используются. При численном задании падающей волны в (6) результаты совпадают с представленными в табл. 2, для аналитического способа задания – приведены в табл. 6.

Существенное отличие (на порядок) результатов из обеих таблиц традиционно объясняется [14, 15] взаимной компенсацией погрешностей двух разностных решений (основного и вспомогательного) задачи (1), используемых в рамках численной методики TF/SF. При аналитическом задании волны в (6) такой компенсации не происходит.

Совместное решение (1) и (3) (восьмая серия экспериментов) сопроводим разбиением D^M на две подобласти DM _ DLM и DRM , где DLM определена на

M   M       M  W  M W     W отрезке [0,L ], DR – на [L +L , 2L +L ], а D длиной LW поместим между ними (рис. 4б). Примем при этом LM = 9,5λ, LW =λ и будем формировать значение ε на том же отрезке, что и ранее. Теперь он совпадает с DW, а левая и правая границы результирующего поля по-прежнему принадлежат DM: левая находится в DLM , правая – в DRM . Таким образом, в рамках одного подхода реализуется технология TF/SF через (2), (6) и возможность исследования дифракции на оптическом элементе (его можно поместить в DW) посредством решения (4).

Отметим отличие результатов, связанных с аналитическим и численным заданием падающей волны (табл. 7), присущее также предыдущей серии экспериментов и вызванное, очевидно, той же причиной.

Таблица 7. Погрешности ε восьмой серии экспериментов

Q, Qt	Задание падающей волны в (6)
	Аналитическое
	QT
	5	10	15
10, 20	3,2787	2,1234	2,0421
20, 40	1,7547	1,0815	0,5867
50, 100	0,4257	0,2717	0,1958
100, 200	0,1527	0,0949	0,0709
Q, Qt	Численное
	QT
	5	10	15
10, 20	0,6168	0,3302	0,1740
20, 40	0,3476	0,0662	0,0537
50, 100	0,0668	0,0189	0,0112
100, 200	0,024	0,0081	0,0046

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

(табл. 6, 7), при численном – несколько десятитысячных (табл. 2, 7). Очевидно, большая разница в первом случае обусловлена неудачной (аналитической) методикой задания падающей волны, которая сама по себе характеризуется значительной погрешностью. Наличие упомянутой разницы (впрочем, весьма несущественной, что свидетельствует о правомерности совместного решения) объясняется причиной, выявленной при анализе третьей серии экспериментов.

Заключение

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

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

Исследования выполнены при поддержке грантов Президента РФ МД-6809.2012.9, НШ-4128.2012.9 и грантов РФФИ №11-07-13164-офи-м-2011-РЖД, № 10-07-00553.