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

УДК 535.42:537.86 П. А. Михеев1

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

В статье описаны два численных метода на основе метода Кирхгофа для решения задачи дифракции на системе отверстий в экране — решение по векторной модели и решение по скалярной модели. Проведено сравнение обоих методов и исследована область применимости скалярной модели.

Ключевые слова: дифракция на отверстиях, численное моделирование, голография, метод Кирхгофа, скалярная теория дифракции, векторная теория дифракции.

1. Введение. Задача дифракции электромагнитного поля на системе отверстий в экране возникает при расчете компьютерно синтезируемых голограмм. Технология компьютерно синтезируемой голографии привлекательна для использования в промышленности, так как одна голограмма способна заменить сложную оптическую систему. Однако практическое применение этой технологии осложнено большой вычислительной сложностью расчетов [1] (при расчете число отверстий в экране может достигать 1011). По этой причине для расчета используют приближенные модели дифракции.

Здесь описаны две приближенные модели, основанные на классическом методе Кирхгофа [2, §8.3], — более точная векторная модель и упрощенная скалярная модель. Скалярная модель не учитывает векторной природы электромагнитного поля и рассматривает его как скалярную величину. Это позволяет существенно упростить формулы и ускорить расчет. Однако из-за меньшей точности скалярная модель может применяться не всегда. Чтобы найти ее область применимости, было проведено сравнение решения по скалярной модели с решением по векторной модели.

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

2. Постановка задачи дифракции. Рассмотрим следующую задачу.

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

Выберем систему координат, как показано на рис. 1, так, чтобы экран находился в плоскости XY, границы отверстий были параллельны осям X и Y, а падающее поле приходило из полупространства Z < 0. Введем обозначения: Cj = (Xi,yi) — координаты центра г-го отверстия, щ — половина размера г-го отверстия по X, bi — половина размера г-го отверстия по Y, a Nr — нормаль к площадке наблюдения.

Обозначим через Tj область (—щ < X — Xi < щ) П (—bi < Y — < bi) П (Z = 0) (г-е отверстие), символом Т — совокупность отверстий IjTj, а символом L — плоскость XY без области Т.

Компоненты /•,' и II падающего поля обозначим через Е\у и //ц • В каждой точке экрана падающее поле аппроксимируем плоской волной:

Ew(pa + Ар) = VE(pQ) ■ А(ра) ■ + ^д^

Hw(pQ + Ар) = Vh(pq) ■ ■ А(ро) • +

V

где I = 1(ра) — вектор распространения волны в точке ро, Ve и V# задают поляризацию, а А — комплексная амплитуда. Потребуем, чтобы в пределах каждого отверстия I, Ve, и А менялись 1 Факультет ВМК МГУ, студ., e-mail: petrmikheevQya.ru

8 ВМУ, вычислительная математика и кибернетика, № 1

незначительно. Обозначим

Е0(р) = ЫР) ■ Жр),

Н0(р) = Ун(р)-,/-.А(р).

[1

Для электромагнитного поля в полупространстве Z > 0 будут выполняться уравнения Максвелла [3, гл. 1, § 1], условия излучения на бесконечности [3, гл. 1, § 5] и приближенные* граничные условия Кирхгофа [4, гл. 8, § 1]. Здесь и далее зависимость Е и Н от времени е~гиí будет опускаться. Получаем следующую систему уравнений для полупространства Z > 0 относительно Е и Н:

го! Ё = ъизуьЙ^ го! Н = —шеЁ,

\ёг х Ё] +

6<г X 6<г

Н]

= о I -г

оо,

Е\ь = 0, Е(р)\т = Е0(р)-е^.

Чтобы завершить постановку задачи, остается выразить искомую интенсивность 1(г) дифракционной картины через Е(г) и Н(г). Это будет действительная часть нормальной компоненты вектора Пойтинга [3, гл. 1, § 3]:

Цг) = Ые ( \[Ё

Н* 1 • Мя

3. Решение задачи дифракции

3.1. Векторная модель. Рассмотрим произвольную точку Я площадки наблюдения с радиус-вектором г. Пусть — поверхность пересечения полупространства Z > Ос шаром радиуса V с центром в Р^ — это часть поверхности, находящаяся в плоскости экрана Z = 0. Соответственно Р^ — сферическая часть поверхности в области Z > О, Ру = Р^ Ру.

Применим для Ру формулы Стреттона-Чу [3, гл. 1, § 3] (формулы Стреттона-Чу выводятся из уравнений Максвелла):

Ё{Е) = - (щЁ) grad^} йв,

Ру

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

Н(Е) = (Ц) | ^гаё ф, [п, Я]] - шеф[Е, п] - (п, Я) grad ^ |

Здесь п — внешняя нормаль к /',-• ф ^-«^(я.м)у¡^.(^(Д, М)); М — точка поверхности, по которой ведется интегрирование, с^(Д, М) — расстояние между двумя точками.

Перейдем к пределу при V ^ оо. Из условия излучения на бесконечности следует [5, гл. 1, § 2], что интегралы по равны нулю. Область /',.' в пределе переходит в плоскость экрана XV. После некоторых преобразований получим

Я(Д) = Я(Д) =

ф \ ( г,к - [К, [п, Е]] - шц[Н, п] - (п, Е)Н ^гк - -

ф\ ( г,к - [К, [п, //]] + ше[Е, п] - (п, Н)Н ^гк - -

йв, йв,

где р = Д ^ М, р = |р|, Н = (Нх,ку,кг) = р/р.

Теперь используем граничные условия Кирхгофа. В приближении Кирхгофа поле на поверхности экрана со стороны Z > 0 вне отверстий равно нулю, так что требуется интегрировать только по отверстиям. В соответствии с граничными условиями в отверстиях будет выполняться равенство Е\т = Ец/. Используя связь между I, Е и Я, получим также, что //1 -/• = Нц/.

Найдем решение Е', Н' для единственного отверстия с размерами 2а х 2Ь, находящегося в начале координат. Общее решение в дальнейшем будет получено суммированием по отверстиям. Формулы для одного отверстия имеют вид

а Ь

Е'(Я) = - 11 ф

— а —Ь а Ь

Я'(Д) = -11'ф

—а —Ь

1\ - ^ - - ^ ^ - - / 1 гк--1 [/г, [п, Ец/}} — ш[г[Нц/,Щ — (п, Ец/)к I гк--

1\ - ^ - - ^ ^ - 1 гк--1 [/г, [п, //„•]] + ше[Ец/,Щ — (п, Нц/)к I гк--

йАхйА у, йАхйАу.

Используем приближенные равенства* гк — р

-1

гк, К ~ ш, р

-1

г 1 и подставим явные значения

для /ц| . П\\ . Будем считать, что в пределах отверстия волновой фронт меняется мало, и вынесем Еа и Я0 за интеграл:

Е'(Я) = Я'(Д) =

1

Ажг 1

Ажг

1\ ^ ^ ^ ^ ^ - ^ / 1 гк--I [т, [п, Еа}} — ш/л[Но,п] — (п, Еа)пг Ык--

гк

1\ ^ ^ - - ^ - _ / 1

[т, [п, Я0]] + п] — (п, Я0)т I гк--

г

а Ь

© = у у е^е^Д^.+Дгг,) ¿АхёАу.

—а —Ь

После дальнейших преобразований получим

#(д) = _ А

4-7ГГ

Я'(Д) =

Ажг

тгЕ0у ^тхЕ0у - туЕ0х_

тгНйу ^тхН0у - туН0х

-Ноу Нах О

~Е0у

Еах О

Е,

Ог

Я

Ог

*Эти приближенные равенства могут применяться при условии г тах(2а, 26). 9 ВМУ, вычислительная математика и кибернетика, № 1

е,

О = е}кг [ е^Н1х-тх)Ах [ егк{1у-ту)Ау ^д _ ^кг8Шаа • sin/ЗЬ

.1 .1 а/3 '

— а —Ь

Здесь а = к(1х — тх), /3 = к(1у — шу).

Теперь проведем суммирование по отверстиям и выпишем окончательные расчетные формулы векторной модели:

N

Ё(г) = Сг,Ё0(Сг),Н0(Сг),\,аг,Ьг^ ,

г=1

N

Н(г) = Яуейог (г - Сг, Ё0(В¿), Я0(с¿), А, Щ, Ь, г=1

/vector (г) = Re ( ^ Е(г) х Н*(г)

NR

где /vector (г) — рассчитанная интенсивность дифракционной картины в точке f, a EvectOT и HvectOT определяются следующим образом:

Evector(г, Eq,Hq, A, a, b) =

iketkr

7гг

HvectOT (г, Eq,Hq, А, а, Ь) :

ikeikr

ИТ

mzE0x / LI mx

mzEQy / HQx + EQZ- my

^mxEQy - myEQx_ V £ 0 mz_

mzHQx J £ ~E0y mx

mzHQy E ox + HQz ■ my

—mxHQy - myHQx_ V M 0 mz

sinaa • sin/ЗЬ

a/3 '

sinaa • sin/ЗЬ

a/3

Здесь используются обозначения г = |f|, т = r/r, a = k(lx — mx), /3 = k(ly — my).

3.2. Скалярная модель. Для сведения к скалярной задаче используем приближение. Будем считать, что во всем полупространстве Z > 0 вектор Е параллелен некоторому единичному вектору Ые, а вектор Н параллелен единичному вектору После этого /•. и II выражаются через одну скалярную функцию II:

/•;^ Л;/, • I// ^ Л;/, • ./-Г.

V

В полупространстве И > 0 для ¿7 выполняются уравнение Гельмгольца [3, гл. 1, § 1], выводимое непосредственно из уравнений Максвелла, условие излучения Зоммерфельда [6, гл. 3, § 3] и граничные условия Кирхгофа:

А[/ + к2[/ = 0,

ди ,1ТГ _(\

—--гки = о - ) , г ^ оо,

дг \г

U\l = 0, U(p)\T = UQ-eik^.

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

наблюдения: Щ = Ау 11 ■ Жд^

Искомая интенсивность дифракционной картины выражается через II по формуле

/зса1аг(г) = ^-|^)|2-

Система уравнений для II линейна, поэтому, как и в случае векторной модели, будем искать решение для единственного отверстия с размерами 2а х 2Ь, находящегося в начале координат. Общее решение в дальнейшем будет получено суммированием по отверстиям.

Рассмотрим произвольную точку Д площадки наблюдения с радиус-вектором г. Пусть /'■ — поверхность пересечения полупространства Я > 0 с шаром радиуса V с центром в Д. /',.' — это часть поверхности, находящаяся в плоскости экрана Z = 0. Соответственно Р; — сферическая часть поверхности в области Я > 0, Ру = / ',.' и

Применим для / ',. интегральную формулу Кирхгофа-Гельмгольца [4, гл. 8, § 1] (это можно сделать, так как выполнено уравнение Гельмгольца):

Pv

Здесь п — внешняя нормаль к /',.; ф ^-«^(я.м)у(^.(^(Д, М)); М — точка поверхности, по которой ведется интегрирование; с^(Д, М) — расстояние между двумя точками.

Перейдем к пределу при V ^ сю. Из условия излучения Зоммерфельда [6, гл. 3, § 3] следует, что интеграл по равен нулю. Область /'.' в пределе переходит в плоскость экрана XV.

Вычислим д'ф/дп. Для этого введем обозначения: р = Ё— М. р = |р|, к = р/р. Имеем

д'ф др дф др 1\ дп дп др дп \ р)

Аналогично выводу для векторной модели используем приближенные равенства* %к — р~1 и %к, к и т, р~1 и г-1.

В приближении Кирхгофа поле и на экране вне отверстия равно нулю, так что теперь требуется интегрировать только по отверстию. Получим

а Ь

и (К) = J ! ф^т^ -1ктхЦ^ йАхйАу.

—а —Ь

Теперь продифференцируем dU/дп:

ЯТГ ятг

(g) = ^,zikUQeik^ll

ятт

U\T(q) = UQ-eik^l\ ^

dU

г

Таким образом, будем иметь

а Ь

и(д) = (lz + mz)Uo ■ [ dAxdAy = -+ mz)UQ ■ (^4eikrS[naa 'f11^) .

47гг J J 47гг \ ap J

—a —b

Проведем суммирование по отверстиям и выпишем окончательные расчетные формулы скалярной модели:

Ui = A(ci)^(l(ci)-NRy),

N

u(r) = ¿/scalar (f- cj, Uq, А, Щ, b^j , i= 1

где /scalar — рассчитанная интенсивность дифракционной картины в точке f, a i/scaiar определяется следующим образом:

тт ikelkr s'maa-s'mßb

гЛса1аг(г, (70,l,A,a,5) =--(lz + mz)Uo--ъ-.

тхт ар

Используются обозначения г = |r|, m = r/r, а = k(lx — mx), ß = k(ly — тпу).

*Эти приближенные равенства могут применяться при условии г тах(2а, 2Ь).

4. Сравнение моделей дифракции

4.1. Случай одного отверстия. Обозначим через Т отношение интенсивности дифракционной картины, рассчитанной по скалярной модели, к интенсивности, рассчитанной по векторной модели. Для случая дифракции на одном отверстии Т удалось представить в виде

т ^scalar & + т^ • (Г- "*) ;

/vector х Рн] •

где Ре и Рн выражаются через Ve, Vh и т.

С использованием этой формулы было доказано, что:

1) при падении волны на экран вдоль нормали при расположенной параллельно экрану площадки наблюдения скалярная и векторная модели вообще не различаются;

2) в случае одного отверстия при переходе от векторной модели к скалярной не смещаются дифракционные максимумы;

3) функция Т различия между моделями не зависит от длины волны, размеров отверстия и расстояния до точки наблюдения.

Было проведено численное исследование зависимости Т от трех параметров — угла между нормалью к экрану и направлением распространения падающей волны, угла дифракции и наклона площадки наблюдения. Результаты показаны в табл. 1 и 2.

Таблица 1

Значения Т при отклонении падающей волны от нормали

Угол дифракции, град Отклонение падающей волны от нормали, град

0 < 10 < 20 < 30 < 40 < 50

0 1.0 (0.99, 1.0) (0.95, 1.0) (0.9, 1.0) (0.83, 1.0) (0.75, 1.0)

< 10 1.0 (0.97, 1.01) (0.93, 1.04) (0.86, 1.08) (0.77, 1.1) (0.69, 1.15)

< 20 1.0 (0.95, 1.02) (0.9, 1.05) (0.82, 1.09) (0.73, 1.2) (0.62, 1.32)

< 30 1.0 (0.95, 1.04) (0.87, 1.07) (0.79, 1.15) (0.7, 1.32) (0.59, 1.56)

< 40 1.0 (0.93, 1.06) (0.87, 1.1) (0.77, 1.21) (0.68, 1.5) (0.55, 4.06)

Таблица 2

Значения Т при наклоне площадки наблюдения

Угол дифракции, град Наклон площадки наблюдения, град

0 < 10 < 20 < 30 < 40 < 50

0 1.0 1.0 1.0 1.0 1.0 1.0

< 10 1.0 (0.99, 1.02) (0.97, 1.03) (0.94, 1.06) (0.92, 1.09) (0.88, 1.12)

< 20 1.0 (0.97, 1.04) (0.94, 1.07) (0.9, 1.11) (0.86, 1.17) (0.81, 1.24)

< 30 1.0 (0.95, 1.05) (0.91, 1.11) (0.86, 1.19) (0.81, 1.28) (0.74, 1.47)

< 40 1.0 (0.93, 1.07) (0.88, 1.16) (0.82, 1.27) (0.76, 1.43) (0.69, 1.77)

4.2. Случай системы отверстий. В случае системы отверстий при переходе от векторной модели к скалярной могут сдвинуться дифракционные максимумы. По этой причине отношение /scalar//vector не годится для сравнения моделей. Вместо него введем функцию различия 8:

jj _ /scalar — /vector

^(/vector)

где S — оператор гауссова размытия с а = А/л/2. Значение S в какой-либо точке площадки наблюдения — это отношение разности рассчитанной по обеим моделям интенсивности к усредненной интенсивности в окрестности этой точки.

Исследовалась зависимость максимальных значений \5\ от показанных на рис. 2 параметров. Для каждого набора параметров была проведена серия численных экспериментов с различными случайными системами отверстий.

Падающая волна Рис. 2. Варьируемые параметры

Рис. 3. Пример распределения 6 по площадке наблюдения для случайной системы отверстий

В качестве примера на рис. 3 показано одно из рассчитанных распределений 6 по площадке наблюдения для апертуры 90°. В центре угол дифракции равен нулю, по углам достигает 22°. Были получены следующие результаты:

1) функция 6 является быстро меняющейся функцией от координат на площадке наблюдения с характерным размером неоднородностей порядка А;

2) максимумы \5\ возрастают с увеличением угла дифракции и апертуры.

Оценки для максимумов \5\ приведены в табл. 3.

Таблица 3

Зависимость максимумов \5\ от угла дифракции и апертуры

Угол дифракции, град

Апертура, град

< 30

< 50

< 70

< 90

< 110

< 1

0.03

0.08

0.12

0.18

0.22

< Ю

0.05

0.1

0.23

0.36

0.47

< 20

0.08

0.2

0.26

0.4

0.59

< 30

0.1

0.23

0.4

0.57

0.68

< 40

0.14

0.27

0.5

0.59

0.82

5. Результаты. В широком диапазоне параметров решение для векторной модели очень близко к точному решению. Было проведено сравнение решения для векторной модели с решением для скалярной. Обычно скалярная модель используется только при небольших апертурах, но представленные в статье результаты свидетельствуют о том, что в некоторых практических случаях ее допустимо применять и при высоких апертурах тоже. Эти результаты имеют большое значение для технологии компьютерно синтезируемой голографии, поскольку вычислительная сложность расчетов очень велика, а скалярная модель позволяет использовать значительно более простые формулы.

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

1. Князьков Д. Ю. Эффективные методы расчета электромагнитных полей // Вычиел. методы и программирование. 2012. 13. № 1. С. 181-188.

2. Born М., Wolf Е. Principles of Optics. Oxford: Pergamon Press, 1968 (Борн M., Вольф Э. Основы оптики. M.: Наука, 1973).

3.Ильинский A.C., Кравцов В. В., Свешников А. Г. Математические модели электродинамики. М.: Высшая школа, 1991.

4. Виноградова М.Б., Руденко О.В., Сухоруков А. П. Теория волн. М.: Наука, 1979.

5. Дмитриев В. И., Березина Н. И. Численные методы решения задач синтеза излучающих систем. М.: Изд-во Моск. ун-та, 1986.

6. Goodman J.W. Introduction to Fourier Optics. N.Y.: McGraw-Hill, 1968 (Гуд мен Дж. Введение в фурье-оптику. М.: Мир, 1970).

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

NUMERICAL SOLUTION OF ELECTROMAGNETIC FIELD DIFFRACTION ON HOLES GRATING

Mikheev P. A.

The article describes scalar and vector numerical methods based on Kirchhoff approximation for modelling electromagnetic field diffraction on holes grating. These methods were compared and range of the scalar method applicability was investigated.

Keywords: diffraction on holes, numerical modelling, holography, Kirchhoff approximation, scalar diffraction theory, vector diffraction theory.