Применение метода конечных разностей для решения задачи дифракции Н-волны на двумерных диэлектрических решетках Текст научной статьи по специальности «Физика»

ПРИМЕНЕНИЕ МЕТОДА КОНЕЧНЫХ РАЗНОСТЕЙ ДЛЯ РЕШЕНИЯ ЗАДАЧИ ДИФРАКЦИИ Н-ВОЛНЫ НА ДВУМЕРНЫХ ДИЭЛЕКТРИЧЕСКИХ РЕШЕТКАХ

Д.Л. Головашкин, Н.Л. Казанский, В.Н. Сафина * Институт систем обработки изображений РАН, * Самарский государственный аэрокосмический университет

Аннотация

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

Введение

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

На сегодняшний день развитие технологии позволяет изготавливать субволновые диэлектрические решетки, период которых сравним с длиной волны, а глубина канавок - с периодом. Традиционно, для анализа дифракции света на таких решетках используется «дифференциальный метод» («differential method») [1], суть которого состоит в определении поля в зоне модуляции посредством решения системы обыкновенных дифференциальных уравнений второго порядка (полученных из уравнения Гельмгольца) с последующей «сшивкой» решения с разложением Релея вне зоны модуляции. Обладая рядом достоинств данный подход не всегда удобен. При формировании системы обыкновенных дифференциальных уравнений приходится ограничиваться определенным количеством элементов ряда разложения диэлектрической проницаемости на периоде решетки, иначе система получится бесконечной. Вместе с тем, число обусловленности матрицы такой системы будет зависеть от периода решетки, а не от формы рельефа. Следовательно, если на периоде рельеф представим в виде быстро осциллирующей или разрывной функции, то, выбирая малое число элементов ряда разложения диэлектрической проницаемости, получим решение задачи для сглаженного профиля, лишь отдаленно напоминающего профиль решетки. А при учете всех необходимых элементов ряда будем иметь дело с плохо обусловленной матрицей, которую невозможно корректно обратить даже известными методами факторизации, применяемыми для плохо обусловленных матриц.

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

для уравнений Максвелла [3], характеризуется многократным (т раз, расстояние в одну длину волны фронт плоской монохроматической волны преодолевает за т шагов по времени) применением быстрого преобразования Фурье (БПФ) для перехода в частотную область. Переформулировав данный алгоритм, авторы избавились от необходимости использования БПФ, существенно сократив длительность вычислений.

1. Применение разностной схемы Уев для моделирования работы дифракционной решетки Ограничимся рассмотрением двумерной диэлектрической решетки, на которую нормально (по направлению Z) падает Н-волна, характеризующаяся напряженностями ЕХ,НУ,Н (рис. 1).

Г2

D /

Г, г3

у

у У — z Г4

Рис. 1. Область вычислительного эксперимента Б и граничные условия Г для задачи дифракции Н-волны на двумерной диэлектрической решетке

Поместим вычислительную область Б между параллельными отрезками Г2 и Г4, отстоящими друг от друга на период ё решетки, и потребуем равенства компонент поля на этих границах ЕХ/Г2=ЕХ/Г4, НУ/Г2=НУ/Г4, Н/Г2= Н/Г4, что соответствует периодичности поля по направлению У. Волна, покинувшая область через границу Г2, появится со стороны границы Г4, как пришедшая от соседнего фрагмента решетки. Такая постановка соответствует решетке, бесконечной по напрвлению У. На Г! зададим условия ввода излучения в область

ЕХ=8ш(юО (1)

- нормальное падение на решетку в момент времени t (te(0,TJ, где T-длительность эксперимента) плоской равномерной Я-волны с циклической частотой ю. Г3 положим электрической стенкой, на которой

Переходя к дискретной модели, наложим на D сеточную область Dh (рис. 2.).

0 12 К-1 К

-Я—>-

ит

- и:

.211—А—II—^

ЛI—А—11—11—А

II—А-

Рис. 2. Сеточная область Квадратами обозначены узлы, в которых определено значение

Ет , треугольниками - Нт ,

х] ,к У],к+1/2

кругами - Нт .

+1/2,к

Тогда сеточная функция Ет задана на об-

х],к

ласти {(^,уА): 1т=шЬ4, т=0, 1, .., М=ШЬ у^у, ]=0,

.., 1=Ьу^у, Zk=khz, к=0, .., K=Lz/hz+1}, Нт - на

у],к+1/2

области {(иуА+ш): tm=mht, т=0, 1, .., М, у^у,

j=0, .., I, zk=khZ, к=0, .., К-1} и Нт - на

+1/2,к

{(иу+1/2А): tm=mht, п=0, 1, .., М, у^у, j=0, .., I, Zk=khz, к=1, .., К-1}.

На заданной сетке запишем разностную схему Уее [3].

Ет+1 _ Ет

S0Sj,k '

j,k_

К

и™+1 ититит

Zj+1/2,k Zj-1/2,k yj,k+1/2 yj ,k-1/2

(2)

j k

К

Em - Em

x1,k Xj

Записанная схема не предусматривает определения проекций Ну Н2 на Г1, Г2 и Г3. На Г4 сетка не налагается вовсе. Следовательно, для корректного задания граничных условий достаточно наложить условия Дирихле только на проекцию Ех при х=0 и х=Ьх.

Как известно, схема Уее характеризуется первым порядком аппроксимации (по пространству и времени) исходной дифференциальной задачи и условием устойчивости [3]:

1

ht <-

1

И2у + И2,

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

Для волны Н-типа проекция электрического поля Ех полностью характеризует всю электромагнитную волну (проекции магнитного поля вьгражаются через Ех). Следовательно, изучая результирующее поле (при 1=Т), можно ограничиться рассмотрением Ех. Исследуя прохождение монохроматической Я-волны через различные среды, поле традиционно представляют в виде комплексной амплитуды

X Хуе X

Тогда

полагая Ex = Re{?x exp(- irat)}.

Ех= Ex Cosrat+ Ex Sinrot.

-Will

Определить комплексную амплитуду поля можно посредством двух замеров Ех в разные моменты времени, решая уравнения:

Ех1= Ex Cosrat1+ Ex, Sinrot1

(4)

ех2= ЕЕх СО8ю/2+ ЕЕх 8шю/2 (5)

х2 хге 2 хт 2 4 '

Условимся брать /1 и /2 =Т такими, что

ю/2=ю/1+п/2. (6)

Учитывая, что Со8ю/1=Б1пю/2 и Б1пю/1=-Со8ю/2, вместо (4), (5) получим:

Ех1= Ех Бтю/2 - Е х СОБю/2,

2 " ^x

д....

-Цо^

ит+1 - ит Em - Em

yj ,k+1 / 2 yj ,k+1 /2 _ xj ,k+1 xj ,k

■j, k

и»+1 - ит Em - Em

Zj+1/2,k Zj+1/2,k _ xj ,k+1 xj ,k .

j ,k

К

h„

причем в узлах j=0 (2) примет вид

Em+1 - Em

К

И^т+1 Ит+1 Ит+1 ит+1

z1/2,k ZJ-1/2, k y0,k+1/2 z0,k-1/2

h„

К

а (3) в узлах j=J запишется как

(3)

Ех2= Ех СО8ю/2+ ЕЕх 8шю/2.

хГе 2 х^т 2

Возведя в квадрат оба уравнения, сложим их:

(Ех1 )2 +(Ех2 )2 = (хГе I2 + (х^т 2=1

где I - искомая величина, пропорциональная интенсивности света.

Величина /2 в (6) не несет в себе смысла Т (или любого другого момента 0<4<Т) из сеточной области. Строго говоря, в вычислительном эксперименте мы не порождаем монохроматическую волну, так как до начала эксперимента излучение в области отсутствовало, но к моменту Т волну, в интересующей нас части области, можно с высокой точностью счи-

тать монохроматической. Для такой волны начало отсчета, от которого откладывается можно выбирать произвольно. В дальнейшем удобно брать ю^=л/2, чтобы Собю^=0, Бтю^=1. Тогда

Exre Ех1 Exim Ех2

(10)

Вышеизложенное описывает переход в частотную область, от Ех к Ex, при котором, как видно, можно ограничиться выбором двух временных слоев, не используя БПФ на последних слоях по времени с дальнейшим усреднением результата, как это предложено в [2]. Такой подход существенно упрощает реализацию алгоритма и сокращает длительность вычислений.

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

Наиболее простой тестовый вычислительный эксперимент заключается в моделировании распространения плоской равномерной электромагнитной волны в вакууме. Положим длину волны Х=1 мкм, размер вычислительной области ограничим величинами Ly=2 мкм, L2=5 мкм. Отведем 100 узлов сетки на один микрон (по направлениям У и Z), тогда ./=200, К=500. Дискретизацию по времени определим из следующих соображений: пусть за время Т волна дойдет от линии 2=0 до 2= L2 (распространяясь по направлению Z), при этом расстояние в один микрон преодолевается за т=100 шагов по времени. Учитывая, что скорость плоской электромагнитной волны в вакууме равна с, зададим Т^2/с, Ы= L2*т. На границе Г3 обратим напряженности поля в нуль, так как за все время эксперимента оно там не появится. Границу Г свяжем с условиями излучения плоской равномерной Т - волны, которые сформулируем как дискретизацию по времени условия (1):

ЕХ=ът(2%*И*с*т/Х). (7)

За погрешность эксперимента в первых трех тестовых примерах примем величину

е = тах

I о - тах I

I о

I о - тт I

I о

х 100%, (8)

где ^ - теоретически ожидаемое значение I в области регистрации. В первом и втором тестовых экспериментах ^ 1 (В/м)2, областью регистрации является вся сеточная область. При выбранных параметрах дискретизации погрешность вычислительного эксперимента составила 0,01% .

Во втором тестовом примере такая же волна распространялась в среде с показателем преломления п=1,5 (оптическое стекло). При той же дискретизации по пространству и Т=n*Lz/c, погрешность возросла до е=0,05%, что вызвано уменьшением скорости волны в среде по сравнению с вакуумом (длина волны описы-

вается меньшим количеством точек по направлению

В следующем вычислительном эксперименте моделируется прохождение нормально падающей волны через плоскую границу раздела среда (п=1,5)/вакуум. Известно, что при нормальном падении доля энергии прошедшей волны в падающей

составляет 1 -

(п -1)2 = (п +1)2

=96%. Далее положим в (8)

!0=0,96*п (В/м) при краевом условии (7).

Выберем L2=3,33 мкм (К=667, 200 узлов на микрон) с таким расчетом, чтобы при t=T и границе сред на отметке 2=1,33 мкм, отраженная волна не успела дойти до 2=0, а прошедшая до 2= Lz. Для т=100 наблюдаем неустойчивость разностной схемы. Увеличивая т до 200, получаем е=0,3%. При т=500 е=0,1%.

На рис. 3 выделим три области: от 2=0 до 2=1,33 мкм, от 2=1,33 до 2=1,83 мкм (область регистрации) и от 1,83 до 3,5 мкм. В третьей области волна еще не устоялась, ее нельзя считать монохроматической. В средней области наблюдаем монохроматическую волну с постоянным значением I, как и полагается прошедшей через границу раздела монохроматической волне.

3 г, мкм

Рис. 3. Значение I при y=const в эксперименте с прохождением нормально падающей волны через плоскую границу раздела среда/воздух Здесь и вычислялись приведенные погрешности. Для первой области характерна интерференционная картина падающей и отраженной от границы раздела волн.

Четвертый вычислительный эксперимент связан с исследованием дифракционной решетки. Была выбрана бинарная решетка с периодом ё=2,5Х, шириной ступеньки а=ё/2 и высотой ступеньки X (рис. 4).

Известно, что непосредственно за периодическим оптическим элементом комплексную амплитуду, как периодическую функцию, всегда можно разложить в ряд Фурье (разложение по плоским волнам):

Е

- (у ) = i

g=-<

С

2п

где

С8 = 7 0 Ех (У2 ехР(_ '8 У'

среда

вакуум

1 2 3 4 5 6 7 8 2,мкм

Рис. 4. Вычислительная область эксперимента с бинарной решеткой. Цифрами обозначены стадии распространения фронта излучения

При этом интенсивности прошедших распространяющихся порядков вычисляются как

|2

18 = Се1 Со^ф .

8

(9)

Угол ф для порядка 8 удовлетворяет формуле решетки

Бтф = 8

(10)

Распространяющимися считаются только те

порядки, для которых 8

< 1, иначе выражение

(10) теряет смысл.

Получив в результате эксперимента Ех (у), вычислим Сг, затем по формуле (9) находим Нормируя интенсивности прошедших распространяющихся порядков, запишем результат: 10=0,055; 11=1-1=0,3729; 12=1-2=0,0727. Найденные значения характеризуют волну, прошедшую через дифракционную решетку.

Исследование той же решетки в рамках скалярной теории света приводит к значениям: 10=0; 11=1-1=0,3711; 12=1-2=0. Если первый и минус первый порядки отличаются от ранее найденых всего на 0,48%, то остальные порядки разняться существенно. Скалярная теория света не позволяет оценить их с приемлемой точностью.

Интересно сравнить данные численного моделирования с результатом, полученным с помощью «дифференциального метода» [1]. Положив в этом методе количество слагаемых, которые участвуют в разложении профиля диэлектрической проницаемости в ряд Фурье, N=5, получим: 10=0,0405; 11=1-1=0,3781; 12=1-2=0,086. Очевидно, для разложения ступенчатой функции пяти слагаемых недоста-

точно. При N=7 (10=0,0549; 11=1-1=0,3558; 12=1_ 2=0,0736) для данного метода перестает выполняться закон сохранения энергии, а при N=17 (10=0,0551; 11=1-1=0,3499; 12=1-2=0,0756) появляется численная неустойчивость. Авторы затрудняется в выборе какого-либо из приведенных результатов для сравнения, представляя читателю возможность самостоятельно определиться с наиболее достоверными данными и сравнить их с числами, полученными посредством разностного подхода.

2. Разностное решение волнового уравнения для моделирования работы дифракционной решетки

Стремясь повысить точность экспериментов, авторы обратились к моделированию распространения Я-волны посредством разностного решения волнового уравнения.

Для неоднородной среды волновое уравнение относительно вектора напряженности электрического поля записывается как [4]:

72 2, еЦ д2Е

V2Е-- 2 2 .

с2 дл2

+ [ 8?ай (1п ц)) х голЕ + +8гаё [ Е • 8гай (1п е) = 0.

(11)

Учитывая, что для Н-волны в выбранной системе координат Е = (Ех ,0,0), среда характеризуется постоянным значением магнитной проницаемости по всем направлениям и неизменным значением диэлектрической проницаемости по направлению X, (11) можно переписать в виде:

V2 ЕхдЕ = 0.

х с2 дЛ2

(12)

Составим простейшую явную разностную схему для (12) на сеточной области Бь реализуя те же краевые условия, что и для схемы Уее:

Ет+1 _ 2 Ет + Ет-1

и2

-1, к>-

-2 е:

(13)

-2 Е"

■2 е:

■2 е:

(14)

■2 е:

У

-е:

с

2

И

У

+

И

z

2

И

2

с

2

+

И

- 2Em + em

h2

■ 2 Em + Em

(15)

- 2Em + Em

к

Схема (13)-(15) характеризуется вторым (в отличии от схемы Уее) порядком аппроксимации (по времени и пространству) и условием устойчивости

[5]

h

<-

1

(при hx=hy=h), аналогичным условию

h сл/2

устойчивости схемы Yee. Сравнивая представленные разностные схемы по вычислительной сложности, следует отдать предпочтение схеме (13)-(15). Полагая длительность операций умножения, деления, вычитания и сложения одинаковой (все они производятся за один такт процессора Pentium 4), получаем, что использование схемы Yee влечет 19%-е увеличение количества арифметических операций по сравнению со схемой (19)-(21).

Успешно проведя тестирование схемы (13)-(15) авторы провели численное моделирование работы бинарной решетки (рис. 4), исследуя зависимости интенсивностей прошедших порядков от ширины ступеньки (табл. 1).

Табл.1. Зависимость интенсивностей прошедших порядков от ширины ступеньки бинарной решетки

Ширина ступеньки дифракционной решетки Нормированные интенсивности дифракционных порядков

Io I1=I-1 I2=I-2

d/2 0,055 0,3729 0,0727

5/8d 0,2 0,2591 0,0884

3/4d 0,5402 0,0861 0,0535

7/8d 0,7917 0,0307 0,0274

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

Выводы

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

новываясь на прямом решении системы уравнении Максвелла или уравнения д'Аламбера. В силу теоремы о сходимости разностного решения к точному можно гарантировать сходимость предложенных разностных схем при сгущении сетки, что выгодно отличает данныи подход от распространенного «диференциального метода». Единственным недостатком предложенных разностных схем являются значительные затраты времени при производстве вычислений (каждый из экспериментов в табл.1. длился в течении 7 часов на компьютере Р-1У-2400 с использованием языка программирования Matlab 5.2). В дальнейшем, для снижения затрат времени авторы намерены использовать условия Беренгера [6] и Мура [7] в качестве поглощающих граничных условий на Г1 и Г3.

Благодарности Работа выполнена при поддержке Министерства образования РФ, Администрации Самарской области и Американского фонда гражданских исследований и развития (CRDF) в рамках российско-американской программы «Фундаментальные исследования и высшее образование» (BRHE), а также гранта Президента РФ № НШ-1007.2003.1 и грантов РФФИ № 01-01-00097 и № 03-01-00109.

Литература

1. Electromagnetic Theory on Gratings: Topics in current physics (22, Ed. By R. Petit, N.Y.: SpringerVerlag, 1980).

2. Hiroyuki Ichikawa, Electromagnetic analysis of diffraction gratings by the finite-difference timedomain method // J. Opt. Soc. Am.-1998. Vol. 15, N. 1. P. 152-157.

3. Yee K.S. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media // IEEE Trans. Antennas Propag. AP-14, 1966. Р. 302-307

4. Неганов В.А., Раевский С.Б., Яровой Г.П. Линейная макроскопическая электродинамика // Под ред. Неганова В. А. Т.1, М.: Радио и связь, 2000. 509 с.

5. Берс Л., Джон Ф., Шехтер М. Уравнения с частными производными // Пер. с англ. Ю.В. Егорова под ред. О. А. Олейник, М.:Мир, 1966. 352 с.

6. Berenger Jean-Pierre A perfectly matched layer for the absorption of electromagnetic waves // Journal of computational physics.1994. N. 114. P. 185-200.

7. Mur G. Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations // IEEE Trans. Elec-tromagn. Compat. EMC-23, 1981. Р. 377-382.

У