Объединенный метод конечных элементов Галеркина и граничных элементов для анализа дифракции ТМ-поляризованной плоской волны на цилиндрических оптических элементах Текст научной статьи по специальности «Физика»

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

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

Аннотация

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

Введение

Для задач моделирования дифракции света на оптических элементах с размерами порядка длины волны свет должен рассматриваться как электромагнитное излучение. Основная часть численных методов моделирования дифракции могут классифицироваться как дифференциальные [1], разностные [2,3], интегральные [4-8], вариационные [9-12], дискретных источников [13], рассмотрение хода лучей [14].

Для определения электромагнитных полей в точке пространства интегральные методы объединяют вклады в эту точку от поля источников по объему или на поверхности. Популярность интегральных методов основывается на их способности решать неограниченные полевые задачи, т. к. условие излучения Зоммерфельда безусловно удовлетворяется в формулировке задачи. Более того, интегральные методы требуют знания поля только на поверхности дифракционного элемента, а не полного поля в пространстве, что минимизирует число неизвестных. В работах [15, 16] представлен объединенный метод на основе метода граничных элементов. В [17] разработан метод расчета дифракции на плоскопараллельной пластине с неоднородностью на основе интегральных уравнений, связанный с численным решением тензора Грина. Недостаток таких методов в том, что они приводят к полностью заполненным матрицам и, следовательно, требуют большего объема компьютерной памяти и длительного времени вычисления. Также границу применения методов необходимо брать по физической границе объекта. Если неоднородность имеет нетривиальную форму, то это приводит к увеличению числа неизвестных.

Разностное решение дифференциальных уравнений Максвелла рассматривалось в работах [18-21, 3]. В работе [22] описано разностное решение волнового уравнения. Недостатками такого подхода являются невозможность использования условий излучения, ограничения на шаги сетки. Для моделирования стационарных задач прохождения излучения разностными схемами используется конечное количество длин волн падающего импульса, что искажа-

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

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

Вариационные методы в задачах с ограниченной областью задачи определяют решения уравнения Гельмгольца путем минимизации функционального соотношения. В работе [10] уравнение Гельм-гольца решалось методом конечных элементов Галеркина с использованием граничных условий сложного вида, зависящих от неизвестного параметра, что потребовало применение границы, определенной формы. Кроме того, данный метод также не включает в себя условия излучения Зоммерфельда.

В работе [24] представлен гибридный метод на основе метода конечных элементов, сформулированного через метод Ритца, и метода граничных элементов. В данном гибридном методе метод конечных элементов применяется для решения уравнения Гельмгольца во внутренней части неоднородного элемента микрооптики, и применяется интегральная методика, метод граничных элементов - к области, внешней к элементу микрооптики, где должны удовлетворяться условия излучения. Оба метода соединяются на границе между внешней и внутренней частями, с удовлетворением условий непрерывности поля. Использование метода конечных элементов для определения поля внутри объекта приводит к матрице трехдиагонального вида, что требует меньше компьютерной памяти и времени вычисления, чем методы объемных интегралов [25]. Результатом использования метода граничных элементов для определения поля на границе является более точное решение, чем применение метода конечных элементов с условиями поглощающей границы. Но применение метода Ритца к решению уравнения Гельм-

гольца некорректно, т.к. он накладывает требование положительности оператора решаемого уравнения. Вывод о знакоопределенности оператора уравнения Гельмгольца сделать нельзя.

В предыдущих работах авторов [26-28] разработаны методы анализа и синтеза элементов двухмерной (цилиндрической) микрооптики, основанные на объединенном методе, который объединяет метод конечных элементов Галеркина (МКЭГ) и метод граничных элементов (МГЭ) для случая ТЕ-поляризации. В данной работе приводится разновидность объединенного метода МКГЭ-МГЭ для случая ТМ-поляризации.

Описание метода расчета для случая ТМ-поляризации Для ТМ-поляризации комплексная амплитуда u(x, у) обозначает полное магнитное поле Н(х, у), которое направлено вдоль оси г (вдоль образующей цилиндрического оптического элемента), координаты (х,у) лежат в плоскости нормального сечения. Полное поле и(х, у) представляет собой сумму рассеянного поля и!!с(х,у) и падающего поля и'"с(х, у).

Полное поле иО(х, у) в О (внутренняя область без граничных точек) должно удовлетворять уравнению Гельмгольца:

1

Р( x, y)

Vun (x, y)

+ ko2q(x, y)uQ (X, y) = 0 :

(1)

где р(х, у)=Цг и q(x, у)=ег для ТЕ поляризации и р(х, у)=ег и q(x, у)=цг для ТМ поляризации. Константы /иг и ег представляют собой отношение магнитной и диэлектрической проницаемостей среды к аналогичным показателям свободного пространства, т.е. /иг=^/^00 и ег=е/е0, к0 представляет собой волновое число волны в свободном пространстве:

ко

'(ho

ч 1/2 = ю = 2п с Х0

(2)

Таким образом, обе поляризации удовлетворяют уравнению (1), но с использованием замены переменных:

Ez ^ Hz, Ц ^ 6Г , 6Г .

(3)

Так как падающее поле и1пс(х, у) известно, то уравнение (1) должно быть решено для рассеянного поля и!!с(х, у).

Рассеянное поле и*с в свободном пространстве ведет себя в соответствии с условием излучения Зоммерфельда

dn

■ — ikoU

(4)

на бесконечности и удовлетворяет уравнению Гельмгольца (1) во внешней области с граничными условиями вида:

ди!С _

dn ^ на поверхности Г,

(5)

1 диО (х, у) ди'пс (х, у) где g (х, у)=--^^--при (х, у) е Г

е дп дп

для ТМ поляризации и при условии прохождения границы Г по линии раздела сред (объект и внешнее пространство).

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

Метод Галеркина вместо решения уравнения Гельмгольца (1) основан на решении соотношений вида:

[[— д[иОс + иОпс ]у - qk2 и + иОпс ] у^ йО = 0, (6)

где у - произвольная функция из области определения уравнения (1).

Используя первую формулу Грина:

Ц рд< ао = | р^П. ш - Ц" vрvQ ао,

о г йп о

для функций Р и <, где О - область плоскости х, у; Г - ее граница, обходимая против часовой стрелки;

dn

производная в направлении внешней нормали

к кривой Г, получим:

Я ( - V [< (X, y) + «г (X, y)] Vy - qk2 [uQc (x, y) -

Q V Р

+uQc (x, y)]y) dQ-

"i

y_ d[«Qc (x, y) + u'Qc (x,y)] df=0 p dn

(7)

Разделяя переменные, можно записать: JJI - V«Qc (x, y )Vy- qk 2»J (x, y)y| dQ -

"J

y d«Qc (x, y)dr =

dn

(8)

1

+

JJ(P V«^ (x,y)Vy - qk«; (x,y)y | dQ+ У duQ (x, y)

dn

-df .

Систему базисных функций для Q обозначим (®Qk,/(x, y)} —0y и систему базисных функций для Г

{afm(x,y)}f=x обозначим {«fy^x,y)} Mf{kJ>l, гдеfk, l),

однозначно сопоставляющая точке (xk, y) точку на Г с порядковым номером m, может иметь вид:

I, если k — 1

N + k, если I — N и k Ф 1

3N + 2 - /, если k — N и I Ф N

4N + 2 -I, если k — N и I Ф 1 и I Ф N.

f (k, I) —

(9)

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

Q

метода Галеркина можно записать систему линейных уравнений:

Аих + Ву*с = - (Аитс + Вутс), (10)

где и=(щ, ..., и + 2 )Т - вектор, составленный из коэффициентов {ищк)+1 = иу) разложения:

иЬ(х, у) = X ик,1 ю?,1(х У).

(11)

Хотя равенство (11) действительно для всех точек (х, у) в 2, необходима отдельная обработка величин поля и его частных производных на Г от значений во внутренней области. Разложение, аналогичное (11), для поля и его частных производных на границе имеет вид:

Л(Х, у) =

= X ик,1 юГи(х у) = X и,г

к, 1=0 т=1

М

Ах, У) = XУтЮ1, (х, У)

, Ю

Г' (т) Г1

(х, у),

(12) (13)

где (х, у) е Г, у=(уь ..., УМ) - вектор, составленный из коэффициентов разложения {ут = Удад) , дик

ук = к

дп

Элементы матрицы А вычисляются по формулам:

аЫу(к)+1, Ыу(1)+]_

г

1

я

Р( х, у )

дю",(х, у) дю",(x,у)

дх

дх

дю" (х, у) дю22 (х, у)

(14)

ду ду

"к? <?(х, у )Юу (х, у)®22,. (х, у)) й 2,

к, / = 1, ..., Ых; 1,. = 1, ..., Му, где 2к,. - треугольная область, включающая узлы сети к и ,.

Элементы матрицы В вычисляются по формулам:

(15)

ьт,р = - Ф ®ттр юр"р й1,

т, р = 1, ... , М,

где Гт, р - линейная область границы Г, включающая узлы границы т и р.

Введем обозначение а= ащк-)+1, Ыу(,)+,-. Матрица А является блочной трехдиагональной .

Каждая из матриц Ак.к является трехдиагональ-ной матрицей порядка N. Более точный анализ показывает, что матрицы {Ак,к-1) Ы=2 и {Ак,к+1) к^1 - двух-диагональные.

Вычислим элементы (а \\) матрицы А и элементы (Ьт^) матрицы В для частного случая кусочно-линейного базиса вида:

'(x, у ) =

1 -1 -1 + 1 + 1+ 1-

хк - х у1 - у

Т

хк- х

Т

у/ - у

Т

хк - х

Т

хк х

Т

у! - у

если х, у е 2

если х, у е 2

если х, у е 2

если х, у

если х, у е Цу_5

если х, у е 2

В силу симметрии А, трехдиагональности матриц {Акк }}= и двухдиагональности матриц

{Ак,к-1}} и {Ак,к+1 }кЫ=11 достаточно указать формул^1

для вычисления элементов:

а к ,1 а к ,1 -1 а к-1,1 а к-1,1+1 а к ,1 , а к ,1 , а к ,1 , а к ,1 ,

1<к, 1<Ы.

Эти формулы, согласно (14) и (16), имеют вид (для простоты используются обозначения

^^ ь^к ,1 ,/•-'•

=я

(

дюк ,1(х у)

дх

(

дюкЛх1У)_ ду

-к0 Ъ ю",12( х, у)) йхйу =

Ц —— йхйу + Ц —— йхйу +

Т

2, и 2^1,2

24И 2,-^4,5

+ Ц —— йхйу + Ц —Р—йхйу

р О1И 26 р1,6

-к,

23 И24 Р3,4

Л ||' -

хк + у1 Г + 2 [1 - хк + у1

<( х + у) + -1 (х + у )2 | йхйу +

+1Х «5 (('+Т I-211+Т)х+£ Ь'+ +л «6 (Í1 - ТТ+2 (■ - Т) у+Т2 Ь

дю",1(х у) дю",1 -1(х у)

ак-1-1 =

я

дх

дх

+ дю",1(х у) дю",1 -1(х у)

ду ду -ко ^ ю",1(х, у)юк>1 -1(х, у)) йхйу =

т=1

а

2

2

2

Г

т Ц ~~ dxdy - k2

Q,U Q6 pi,6

íí qi ÍÍ1

h+lL If 1 + У-! |+1 (i + У-! | x + h Jl h J h( h

+i ( x + У' + У'-1

h V h

►Я «6 Ifi

■j У - ht ( + У2 Í} dxdy +

+xt+yvi |(i - £ l-1Í1 - У. | x+ h Jl h J h V h

hVJy - F((x+y2 Z}dxdy

-,l — íí

droQ,i(x, y) droQ-i,i(x y)

ôx

ôx

droQ,i(x, y) 5roQ-i,i(л; y)

dy

dy

- kó2 qi roQ,i( x У)>

-i,i(x,У)Zdxdy — -íí -p—dxdy -k(

Q,U Q2 pi,2

kó x

íí «'ÍÍ1 -

^JÍ'+f J+i('+Y i-y+

+h VJ x -¿(>x+x: Z}dxdy + +f «2 |(' +

xk-1 + (i lh x

1 -ZL |-1 (i + iL |y +

h J h V h

i ( J x - F ( + x! Z} dxdy

(17)

Элемент а к,1+1 находится по формуле для вычисления а к''1 -1, но интегрирование производится по областям О и О \ 14, элемент а \+1'1 находится по формуле для вычисления а \-1'1, но интегрирование производится по областям О \ 14 и О \ 15, элемент

' k ,i , k-i,i+i

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

тям « i,5 и « i,6-

Отсюда сразу следует, что матрицы {Лк,к-1} f=2 и ^^¿.t+i} fji' являются диагональными. Элементы матрицы В вычисляются по формулам:

bm.m= Ф raLdl= á>\ 1 -^^ 1 dx+

гh

Г'

h J

+ ^(l + xL_x J dx = |h,

bm,m-1= Ф rom rom-' dl=

= <j>jl -

x x jí1+xk-^ |dx =

h JV h

(x - xk-i)(xk - x)

h

i

dx = —h, 6

Г'

bm,m+1 bm,m-1.

Тогда вместо системы уравнений (ió) получим:

'[ A.J [ A^. J ó

[ A. [ A^J [B]

[ A.J [ A^. J ó

[ A. ,J [ A^J [B]

(19)

В равенстве (19) коэффициенты связи между значениями поля только во внутренних узлах представлены подматрицей [Ах,8] размерностью (Ы-М)х(Ы-Ы), где N - количество узлов разбиения области О, М - количество узлов на границе Г. Связь между величинами поля только на граничных узлах представлена подматрицей [АГ, Г] размерностью МхМ, и связь между величинами поля на внутренних и на граничных узлах представлена подматрицей [АГ,Х] размерностью (^М)хМ и транспонированной [АХ,Г] размерностью (^М)хМ. Взаимодействие между производными величин поля на граничных узлах представлено подматрицей [В] размерностью МхМ.

К сожалению, система уравнений (19) не имеет единственного решения, так как она состоит из N равенств с N+M неизвестными: N узловых величин поля и*ск,1(х, у) и М производных по нормали на граничных узлах уХк,(х, у).

Для построения граничного интегрального уравнения для поля и его нормальной производной воспользуемся теоремой Грина в следующем виде:

í [V2u* (I, n) + k2u* (I, n) J u(n)dQ(x) —

Q

— -í q(n)u* (I, n)d Г( x) + í u(n)q* (I, n)df( x).

(2ó)

Функции в обоих интегралах в правой части

уравнения (2ó) - это q(n) —

du(n) dn

- нормальные про-

изводные амплитуды поля.

Пусть функция u является фундаментальным решением уравнения Гельмгольца:

V2u * (I, n) + к 2u * (|, n) = -8(|, n). (21)

Фундаментальное решение для уравнения Гельмгольца в двумерном однородном пространстве известно и равно

u* — (i /4)(kr)

(22)

x

Г

Q

Q

a

2ó

где г = 7[(П) - х1(|)]2 +[(п) - х^)] 2 , Н 0« (г) =

^(г) + /Т0(г) - функция Ханкеля первого рода и нулевого порядка, где ^ - функция Бесселя первого рода нулевого порядка, У0 - функция Неймана нулевого порядка.

Подставляя равенство (21) в уравнение (20) и осуществляя предельный переход точки наблюдения % из внутренней в граничную, найдем

c(5)u(5) +|u(n)q (5, Г(x) =

г

= J q(n)M* (5, n)d Г( x).

(23)

Это уравнение обеспечивает функциональную связь между функцией и и ее нормальной производной « на границе Г. Функция с в уравнении (23) равна:

c(5) = 1 -±Ф, 2п

(24)

где ф - внутренней угол кусочно-линеинои границы в точке 5.

Подставляя вместо функции комплексной амплитуды на границе ее аппроксимацию базисными кусочно- линейными функциями, получим вместо (23):

м (Г 1

cm< = X h { < J0 ®Г ( + [Ps+1 - Ps ]y) X

i=i vL

JU" (Pm , Ps + [C s+1 - C s ]Y)

дп'

-d у +

+< J0 ®г (сs-1 + [с s - сs-1]у)

Xdu (с m , с s-1 + [с s - с s-1]Y )

дп'

1

J®r (C s + [C s+1 - с s ]Y)

d y

(25)

X и* (с т , С , + [С ,+! - С , Ю й| +

1

+{< (с,-1 +[с* -с*-1]у)

0

Х и (С т , С *-1 + [С , - С *-1]У) й У]} ,

т = [1, М ].

Равенство (25) может быть представлено в матричном виде:

[Б] < +[Е] уГс = 0, (26)

где элементы матриц [Б] и [Е]:

Йтд =- Ю"(С, + [С*+1 - С* ]У)Х

Хди* (с т , С * + [С *+1 - С * ]У)

дп'

-d y-

-hj0 ®Г (С s-1 + [С s - с s-1]Y)

(с m , с s-1 + [с s - с s-1]5)

(27)

дп'

d y + c 5

I m m

em,s = +h

1

j®,r(c s +[с s+1 - с s ]y)

X U (C m , C s + [C s +1 - C s ]Y)dY +

1

+Jror(C s-1 + [C s - C s-1]Y)X 0

X U (Cm , C s-1 + [C s - C s-1]Y)dY] .

Интегралы в равенствах (27) и (28) могут быть оценены численно.

Принимая во внимание скачок нормальной производной поля на реальной границе рассеивающего объекта для падающей ТМ поляризованной

волны

+ v^

и полагая для определенности относительную диэлектрическую проницаемость внешней к объекту

11пс тс

, и учитывая равенство у1 = у2 , можно

записать

S. V 2 S9 V,

v'c — ZL\r'C 12 V1 _ 2

vsc =S1 < +(S1 -1) v"

или

vsc - (S1 - 1)v

(29)

Тогда, объединяя системы уравнений (19) и (26), получим замкнутую систему линейных алгебраических уравнений для решения задачи дифракции плоской ТМ поляризованной волны на двухмерном (цилиндрическом) диэлектрическом элементе микрооптики:

[ Ах,, ] [ Лг,х ] 0 " [А х ,г] [Лг>г] [В]

0 [О] [Е]

L A*, * ] L A г,* ] 0

La*,г] LAr,r] [B]

[E]

(30)

где - это область 2 без своей границы Г, подматрица размерностью ЫхЫ включает в себя коэффициенты вкладов внутренних узлов, подматрицы А5Г и Аг,5 размерностью ЫхМ и МхЫ соответственно включают в себя коэффициенты вкладов внутренних и граничных узлов, подматрица АГ, Г размером МхМ включает в себя коэффициенты вкладов граничных узлов, и5 и иГ - вектора напряженности поля во внутренних и граничных узлах сети.

e

m,s

v " =

2

а

Система линейных уравнений (30) может быть представлена в матричной записи:

тис=иии, (31)

и-=Т"1ииш, (32)

где и"с=

"с 1"

"с 1"

иг , и = иг

"с 1"

V V

для случая ТМ поляризации

[ Ах, ] [ А^ ] 0 [ А, ,г] [Лг,г] [В]

т=

[В] 8-[Е]

и = -

[ А ,,, ]

[ А , ,г]

[ Аг.х ] [ Аг,г]

0

[В]

0

0

1

- 1|[Е]

Число коэффициентов в матрице Т в системе уравнений (31) с учетом трехдиагональности матрицы А составляет Ы1Ы2(12М+ 10Ы2) + 16(М+ Ы2)2, где М, Ы2 — число узлов разбиения прямоугольной области О по каждой стороне прямоугольника. Количество операций, требуемых для решения системы (31), можно оценить как О(Ы12Ы2(12Ы1+ 10Ы2)).

С помощью объединенного метода МКЭГ-МГЭ, представленного системой линейных уравнений (30), (31), производится расчет рассеянного поля как вдоль границы, так и во всех внутренних узлах. Для определения рассеянного поля во внешних по отношению к объекту точках, величины поля в граничных узлах линейно интерполируются, используя равенства (12), частные производные поля в граничных узлах линейно интерполируются, используя равенство (13), и интерполированные данные подставляются в уравнение (23), но £ уже является внешней точкой, и функция с(£)=1.

Полное поле представляет собой сумму рассеянного и падающего полей [29]:

и (х, у) = и'" (х, у) - ф

ди"

ди* дп

й г,

(33)

дп

(х, у) ¿П .

Для ТМ поляризации уравнение (33) в матрич-

ной записи имеет вид:

и=и"^и"с+

' 1 -8, ^

Иуш=и"^Т-1

хЦи'" +

' 1 -8 ^

/

/

ИУ"

(34)

где матрица W имеет следующий вид:

W =

0

[г] 1 [н]

81

где элементы матриц Б и И находятся из уравнений

/р (x, у) = юр (с р +[с р+1 - с р ]у)х ди*((х y), с р +[с р+1 - ср ]у)

дп

-й у +

(35)

й у,

у)

(36)

+40 юр (с р-1+[с р- с р-1]^)х

ди*((x, y), с р-1 +[ср - с р-1]У)

х-------

дп

кр(x, у) = - к{0' юр (с р +[с р+1 - с р

хи * ((х, у), с р + [с р+1 - с р ]у) й у-

-к|<!юр (ср-1 +[ср -ср-1]^)х хи * ((х, у), с р-! + [с р - с р-1 ]у) йу, р = 1, ..., М.

Численные примеры Рассмотрим дифракцию плоской ТЕ-поляризованной волны на диэлектрическом однородном цилиндре с круглым сечением. Пусть плоская волна падает на цилиндр слева направо вдоль оси х, длина волны 10=1 мкм. Относительная диэлектрическая проницаемость цилиндра е=2, радиус равен 0,5 мкм. Окружающее цилиндр однородное пространство имеет параметры 8=^=1. Расчет производился объединенным методом МКЭГ-МГЭ [26-28]. В качестве области П выбиралась квадратная область размером 1х1 мкм, а контура , - периметр области П. После нахождения проекции на ось г электрического вектора Ег внутри области П и на ее периметре ,, рассчитывалось поле вне области П в области Т с помощью теоремы Грина по формуле (33).

На рис. 1 представлена зависимость от числа точек разбиения П средних квадратичных отклонений сформированного поля от поля, рассчитанного аналитически [28], по интенсивности и амплитуде. В суммах из функции Бесселя и Ханкеля, аппроксимирующих комплексную амплитуду, учитывалось 30 слагаемых, что давало погрешность порядка 0,001 %. Средние квадратичные отклонения по интенсивности и амплитуде вычисляются соответственно по формулам:

81 =

1

£ ¡1 (х,, у,) -Гс (х,, у,)|2

£[ 1'"с (х,, у,)]2

(37)

8А =

£1 |и( хк, у1 ^ - |и'"с (хк, у,

£[| и'"с (х,, у, )| ]2

(38)

Ошибка при уменьшении расстояния между узлами к, равного отношению размера области расчета I к числу узлов дискретизации Ы, до 0,03Х быстро спадает и далее медленно уменьшается.

X

5 амплитуды ........8 интенсивности

40 50 60 70 80 Число узлов разбиения

Рис. 1. Зависимость среднего квадратичного отклонения амплитуды и интенсивности рассчитанного поля (от аналитически рассчитанного поля) от дискретизации сетки У противоположного падению волны края цилиндра наблюдается фокусировка поля, максимум которого лежит на границе раздела сред. На рис.2 приведена зависимость значений максимума интенсивности в зависимости от числа узлов дискретизации. Из графика видно, что значения максимума с увеличением числа узлов дискретизации стремятся к значению максимума, рассчитанному аналитически.

— численное решение ^ ........аналитическое решение

30 40 50 60 70 80 Число узлов дискретизации Рис. 2. Зависимость максимума интенсивности рассчитанного поля от дискретизации сетки Распределение амплитуды поля дифракции на рассчитываемом диэлектрическом цилиндре показано на рис. 3. Число точек разбиения области О было равно Ж=50х50. Распределение амплитуды по оси х, проходящей через центр цилиндра, показано на рис. 4. Цилиндр располагается на координатах от 0 мкм до 1 мкм.

Ниже для сравнения приводятся результаты расчетов для ТЕ- и ТМ-поляризаций.

Рис. 3. Распределение амплитуды поля дифракции плоской ТЕ-поляризованной волны на круговой цилиндр, полученное объединенным методом

О

--- 1 1 1 1 I /1 /| N 1 1 . I ---

/\! 1 'V/ |_ 1

--- V ¡ч/47 1 V 1 1 1 1 ---

-1 0 1 2 3

X, мкм

Рис. 4. Распределение амплитуды по оси х для рис. 3 Расчет поля дифракции ТЕ-поляризованной волны на диэлектрическом цилиндре с квадратным сечением 1х1 мкм и относительной диэлектрической проницаемостью е =2 представлен на рис. 5. Число точек разбиения области О равно Ж=50х50. Распределение амплитуды по оси х, проходящей через центр, показано на рис. 6. Цилиндр располагается от 0 мкм до 1 мкм.

Рис. 5. Распределение амплитуды поля дифракции плоской ТЕ-поляризованной волны на квадратном цилиндре, полученное объединенным методом

О

1 _ 1 1 1 ---

/| ---

Л1- 4 \ 1

1 у 1 \/ | V 1 1 1 ---

0 1 2 3

-1

X, мкм

Рис. 6. Распределение амплитуды по оси х для рис. 5

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

Рассмотрим дифракцию плоской ТМ-поляризо-ванной волны на диэлектрическом однородном цилиндре с круглым сечением. Пусть плоская волна падает на цилиндр слева направо вдоль оси х, длина волны Л0 = 1 мкм. Относительная диэлектрическая проницаемость цилиндра е=2, радиус равен 0,5 мкм. Окружающее цилиндр однородное пространство имеет параметры е=^=1. Расчет производился объединенным методом МКЭГ-МГЭ по формулам (32). В качестве области О выбирался круг радиусом 1 мкм, а контура - периметр области О. Число точек разбиения на 1 мкм бралось 50. После нахождения проекции на ось г магнитного вектора И внутри области О и на ее периметре 5, рассчитывалось поле вне области О в

области Т с помощью теоремы Грина по формулам (33). Распределение амплитуды показано на рис. 7. На рис. 8 показано распределение амплитуды вдоль оси х. Начало координат находится на левом крае цилиндра. Рассчитанное распределение амплитуды совпадает с аналитическим решением с ошибкой по амплитуде 5 =2%.

Рис. 7. Распределение амплитуды светового поля при дифракции плоской волны ТМ-поляризации на круглом цилиндре

Рис. 8. Распределение амплитуды по оси х для рис. 7

Рис. 9. Распределение амплитуды поля дифракции плоской ТМ-поляризованной волны на квадратном цилиндре, полученное объединенным методом

Рис. 10. Распределение амплитуды по оси х для рис. 9

Расчет поля дифракции ТМ поляризованной волны на диэлектрическом цилиндре с квадратным сечением 1х1 мкм и относительной диэлектрической проницаемостью е=2 представлен на рис. 9. Число точек разбиения области П равно N=50x50. Распределение амплитуды по оси х, проходящей через центр, показано на рис. 10. Цилиндр располагается от 0 мкм до 1 мкм.

Картины распределения амплитуды поля для цилиндров с круглым и квадратным сечением отличаются по характеру отраженной волны. В случае квадратного сечения рассеяние назад меньше, что можно объяснить влиянием кривизны границы среды. Максимальное значение амплитуды поля также наблюдается в случае круглого сечения. Можно заметить, при ТМ-поляризации максимум интенсивности находится уже не на границе раздела сред, как в случае ТЕ-поляризации, а внутри цилиндра.

Заключение

В работе полученные следующие результаты.

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

Работоспособность метода продемонстрирована хорошим согласием (1-2%) полей дифракции плоской волны на круглом цилиндре, рассчитанных данным методом и по известным аналитическим формулам.

Литература

[1] Montiel F., Neviere M. Differential theory of gratings: extension to deep gratings of arbitrary profile and permittivity through the R-matrix propagation algorithm // Journal of Optical Society of America. - 1994.-Vol. 11.-P. 3241-3250.

[2] Taflove A. Computational electromagnetics: the finite-difference time domain method // Artech House, Boston. - 1995.

[3] Головашкин Д.Л., Дегтярев А.А., Сойфер В.А. Моделирование волноводного распространения света оптического излучения в рамках электромагнитной теории // Компьютерная оптика, 1997, №17. - с.5-9.

[4] Brebbia C.A. The boundary Element Method for Engineers // Pentech Press, London; Halstead Press, New York, 1978 (Second edition, 1980).

[5] Choi M. K. Numerical calculation of light scattering from a layered sphere by the boundary-element method // Journal of Optical Society of America. -2001.- Vol. 18.-N. 3.-P. 577-583.

[6] Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов // Пер. с англ. -М.: Мир, 1987. - 524 с.

[7] Colton D., Kress R. Integral equation methods in scattering theory, J.Willey&Sons, New York, 1983, 311p.

[8] Бреббия К., Уокер С. Применение метода граничных элементов в технике: Пер. с англ. - М.: Мир, 1982. - 248 с.

[9] Davies J. B. Finite element analysis of waveguides and cavities - a review // IEEE Trans. Magn. - 1993. - Vol. 29. - P. 1578-1583.

[10] Lichtenberg B., Gallagher N. Numerical modeling of diffractive devices using the finite element method // Optical Engineering. -1994. - Vol. 33. -No. 11. P. 3518-3526.

[11] Михлин С.Г. Вариационные методы в математической физике. М.: Техн.-теорет.лит., 1957, с.476.

[12] Сьярле Ф. Метод конечных элементов для эллиптических задач. - М.: Мир, 1980. -512 с.

[13] Blaike R. J., McNab S. J. Evanescent interfer-ometric lithography // Applied Optics. - 2001. - Vol. 40.

- N. 4. - P. 1692-1698.

[14] Voznesensky N. Simulation model for light propagation through nanometer-sized structures // Optical Memory and Neural Networks. -2000. -Vol. 9. -N. 3. -P. 175-183.

[15] Prather D. W., Mirotznik M. S., Mait J. N. Boundary integral methods applied to the analysis of diffractive optical elements // Journal of Optical Society of America. - 1997. - Vol. 14. - P. 34-43.

[16] Tanaka M., Tanaka K. Computer simululation for two-dimensional near-field optics with use of a metal-coated dielectric probe // Journal of Optical Society of Amereca. - 2001.- Vol. 18.-N. 4.-P. 919-925.

[17] Paulus M., Martin O. J. F. Light propagation and scattering in stratified media: a Green's tensor approach // Journal of Optical Society of Amereca. - 2001.

- Vol. 18.-N. 4.-P. 854-861.

[18] Dou W. B., Yung E. K. N. Diffraction of an electromagnetic beam by an aperture in a conducting screen // Journal of Optical Society of America. - 2001.-Vol. 18. -N. 4.-P. 801-806.

[19] Lee J.-F., Palandech R., and Mittra R. Modeling three-dimensional discontinuities in waveguides using nonorthogonal FDTD algorithm // IEEE Trans. Microwave Theory Tech. - 1992. -Vol. 40. -P. 346-352.

[20] Prather D. W., Shi S. Formulation and application of the finite-difference time-domain method for the analysis of axially symmetric diffractive optical elements // Journal of Optical Society of America. - 1999, Vol. 16.-N. 5.- P. 1131-1142.

[21] Shi S., Tao X., Yang L., Prather D. W. Analysis of diffractive optical elements using a nonuniform finite-difference time-domain method // Opt. Eng.-2001.-Vol. 40.-N. 4.-P. 503-510.

[22] Gruzdev V., Gruzdeva A. Finite-difference time-domain modeling of laser beam propagation and scattering in dielectric materials // Proceedings of SPIE. -2001. -Vol. 4436. -P. 27- 38.

[23] Berenger G.P. A perfectly matched layer for the absorption of electromagnetic waves // Journal of Computational Physics. - 1994. - vol. 114. - P. 185-200.

[24] Mirotznik M., Prather D., Mait J. A hybrid finite element-boundary element method for the analysis of diffractive elements // Journal of Modern Optics. -1996. - Vol. 43. - N. 7. - P. 1309-1321.

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

[26] Kotlyar V.V., Nesterenko D.V. A finite element method in the problem of light diffraction by micro-optics // Optical Memory and Neural Networks. -2000. -Vol. 9. -No. 3., p. 209-219.

[27] Kotlyar V.V., Nesterenko D.V. Analysis of light diffraction by binary micro-optics using a combination of boundary element method and finite element method // Proceedings of SPIE. -2001. -Vol. 4242. -P. 125 - 132.

[28] Котляр В.В., Нестеренко Д.В. Дифракция электромагнитной волны на круговом диэлектрическом цилиндре: расчет по аналитическим формулам и методом конечных элементов - граничных элементов // Физика волновых процессов и радиотехнические системы, Самара, ПГАТИ, том 3, №3-4, 2000. -С. 25-28.

[29] Colton D., Kress R. Integral equation methods in scattering theory // John Wiley&Sons, New York, 1983.