Анализ дифракции света на элементах цилиндрической микрооптики объединенным методом конечных элементов Галеркина и граничных элементов

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

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

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

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

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

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

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

Описание метода расчета

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

определить полное или рассеянное поле, характеризующие структуру.

Любое двумерное поле может быть разложено на E z -поляризованное и H z -поляризованное поля. В области дифракции поле описывается системой дифференциальных уравнений, различных для случаев ТЕ- и ТМ-поляризаций. Для TE-поляризации ( E ( x , у ) = ( 0,0, Ez ( x , у ) ) ) комплексная амплитуда u ( x , у ) обозначает полное электрическое поле E z ( x , у ), которое направлено вдоль оси z (вдоль образующей цилиндрического оптического элемента), координаты ( x , у ) лежат в плоскости нормального сечения. Для ТМ-поляризации ( H ( x , у ) = ( 0Д Hz ( x , у ) ) ) ком

В рассматриваемой задаче область расчета бесконечна. Однако, как известно, МКЭ применим только к конечной или ограниченной области. Таким образом, чтобы решить уравнение (1), бесконечная область Т, внешняя к рассеивателю, должна быть ограничена введением искусственной границы Г. Соответственно, для единственного решения задачи на данной искусственной границе должны быть введены граничные условия. Такие условия должны сделать границу как можно более прозрачной для рассеянного поля или, другими словами, должны минимизировать нефизические отражения от границы. Один из классов граничных условий, разработанных для этих целей, может быть получен из граничных интегральных уравнений, применяемых к внешней области. Эти граничные условия являются глобальными по своей природе, т.е. они соотносят поле в одном граничном узле с полем на всей границе. Данные граничные условия предотвращают отражение на границе для всех углов падения волн и приводят к точному решению.

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

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

И I ^- “ ^{A и} o Y ^- qk ² u о У" f i Y| ^{d o} = °.           ⁽³⁾

О ( P                   )

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

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

плексная амплитуда u ( x , у ) обозначает полное магнитное поле H z ( x , у ).

Полное поле u _О ( x , у ) в области О должно удовле

JJ P A Q d O = J P d Q d l = JJ V P V Q d O , о           г   d n      о

творять уравнению:

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

V-

1 p ^{( x} , у )

^V u о ^{( x} , у )

+ k 0 ² q ^{( x} , у ) и о ^{( x} , у ) = f о

d Q d n

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

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

где fо= jk0 Z0 J , p(x, у)=ц r , q(x, у)=Б r для TE-поляризации и д Г 1 ) д ( 1 )

f о =—^1 — J 1+^-1 — J x I , p ^{( x} , у )= ^е r , q ^{( x} , у )= ц r ^{д x} (^Б r ) Р у (^Б r )

ff I - ^V u о ( ^x , у ) ^VY ^- qk ² u о ( ^x , у ) Y I O ( p                           )

_JX du n^^ d _r = j_{f /} o Y d n .

г p   d ⁿ         n

для TM-поляризации. Константы ц _r и б _r представляют собой отношение магнитной и диэлектрической проницаемостей среды к аналогичным показателям свободного пространства, т.е., ц _r= ц / ц ₀ и б _r= б / б ₀, к ₀ представляет собой волновое число волны в свободном пространстве

/      X 1/2 to 2п к0 =ю(Ц0Б0) =- = —,                    (2)

c д

Z 0 =7ц 0 / Б 0 - импеданс свободного пространства, J - плотность электрического тока источника.

Систему базисных функций для О обозначим О ,л) ^{^} x ■ ^N y

{ to kl ( x , у ) } kl 0 и систему базисных функций для Г обозначим { to { , ( x , у ) }" 1 , где N x , N y - число узлов сеточного покрытия прямоугольной области О по оси x и у, соответственно, М - число узлов сеточного покрытия границы Г.

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

A u + B v = C f ,

где u = ( u 1 , ..., и

-.Ny У

- вектор, составленный из

N x ’ N y коэффициентов { u n ( k )+ _; = u k1 }      разложения:

где Q _k , j - область разбиения области Q, включающая узлы сети k и j .

В качестве кусочно-линейного базиса была определена функция вида:

^toQ , 1 ( ^x . y ) =

Nx, Ny uQ( x, y ) = E uk, i ®“/(x, y )■                     (6)

k , 1 = 0

Вектор f = ( f ,..., f NxNy ) T - вектор, составленный из коэффициентов разложения:

Nx, Ny f Q( x, y )= E fy ®Q1(x, y).                   (7)

k , 1 = 0

Хотя равенства (6) и (7) действительны для всех точек ( x , y ) в области Q , необходима отдельная обработка величин поля и его частных производных на границе Г от значений во внутренней области. Разложение, аналогичное (6) и (7), для поля и его част-

1 -	xk - h	x	y 1 - y h ’		если x, y е Q h1 1
1 —	x k - h	x			если x, y	eQ h,, k,1,
1 +	y 1 - h	y			если x, y	eQ h.1.
1 +	x k - h	x	+ yl-h	y ,	если x, y	eQ h.1. 4
1 +	x k - h	x	,		если x, y	eQ h ,5
1 -	y 1 - h	y	,		если x, y	eQ h,1, 6

ных производных на границе имеет вид:

M u Г( x ’ y ) = E u m ® m ( x ’ y )’ m = 1	(8)
M v Г ( x ’ y ) = E v m ® m ( x ’ y )’ m = 1	(9)
M f Г ( x ’ y ) = E f m ® m ( x ’ y )’ m = 1	(10)

где h - длина сегмента сети покрытия.

Элементы ( a k ’ j ) матрицы А , элементы ( b_m , s ) матрицы В и элементы ( c k ^{, j} ) матрицы С вычисляются из уравнений (11), (12) и (13), соответственно. Тогда систему уравнений (10) можно представить в виде:

[ ^A Q , q ] [ ^А Г,П ]   ⁰

.[ A Q,r ]   [ A г , г ] [ В ]

где ( x , y ) e Г, v= ( v 1 , ..., v_M ) T - вектор, составленный из коэффициентов разложения v_k = д u_k /д п .

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

a Ny ( k )+ 1 , Ny ( i )+ j

	f	dtoQ ( x ’ y ) dto Q j ( x ’ y ) ,
= JJ Q k , 1	1	д x       д x
	p ( x ’ y )	+ дю Q > 1 ( x , y ) дю Q j -( x ’ y )	(11)
	(	['     д y         д y    ]

^- к 2 q ⁽ x , y ) m “⁽ x , y ) m _kj ⁽ x , y ) ) d ^Q ,

k, i =[1, Nx ], 1, j = [1, Ny ], где Qk, j - область разбиения области Q, включающая узлы сети k и j.

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

b ms = - f ® m ®^{Г d 1} ’

Г m,                                           (12)

m , s = [ 1, M ] ,

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

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

c Ny ( k ) + 1 , Ny ( i ) + j = JJ ^to k , 1 ⁽ x ’ y Ж; j ⁽ x , ^y ) ^{d Q} ,

^Q k,1                                           (13)

k , i = [ 1, N x ] , 1 , j = [ 1, N y ] ,

Система уравнений (15) не имеет единственного решения, так как она состоит из N равенств с N+M неизвестными: N = N_xN_y - общее число узловых величин поля u_k , ₁ ( x , y ) в области Q и M производных по нормали на граничных узлах v_k , ₁ ( x , y ).

Определим поля в свободном пространстве Т (вне области Q с ее границей Г). Так как это пространство однородное, то поле может быть сформулировано в терминах граничных интегралов с соответствующей функцией Грина. Полное поле u _Т( x , y ) в области Т должно удовлетворять уравнению:

V-

- V u _¥ (У ^p

+ k 0 ² qu _v (У = ^f г ’

^еУ,

где f у = jk 0 Z 0 J _z⁴' , p ( x , y ) = и r , q ( x , y ) = e r для TE-поляризации и

= —Ef j у ^{д x} L^e r y

Wf J^

) д y (e r )

p (x ’ y ) = e r, q (x ’ y ) = и r для TM-поляризации. J Т - плотность электрического тока источника в свободном пространстве.

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

дем функцию Грина и * , удовлетворяющую условиям излучения Зоммерфельда и являющуюся фундаментальным решением уравнения Гельмгольца:

V ² и * ( L , п ) + к ² и * ( L , п ) = -5 ( L п ), пеУ . (17)

Фундаментальное решение для уравнения Гельмгольца в двумерном однородном пространстве известно и равно и * = (i/4) H 01)( кг),                                  (18)

H 0¹) ( kr ) = J ₀ ( kr + iY₀ ( kr ) ) - функция Ханкеля первого рода и нулевого порядка, где J ₀ - функция Бесселя первого рода нулевого порядка, Y 0 - функция Неймана нулевого порядка.

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

J [V ² и * ( L , п ) + к ² и * ( L , п ) ] и ( п ) d V = у

= ^- J q ⁽ п ) и * &  пЖ +                      (19)

Г

⁺ J и ⁽ п ) q * ( i пЖ⁺ J f v ⁽ п ) и * ( i n )^{d V} Г                   V

Функции в обоих интегралах в правой части уравнения (19) - это q(п) = dи (п)/dn' - нормальные при ^еГ (8) и (9), получим соотношение (23), которое может быть представлено в матричном виде (24) с элементами матриц [D] и [G] вида (25) и (26).

Интегралы в равенствах (25) и (26) могут быть оценены численно. Объединяя системы уравнений (15) и (24), получим замкнутую систему линейных алгебраических уравнений для решения задачи дифракции плоской волны на цилиндрическом микрообъекте (27).

cu mm

M

= Z h s=1

■ 1

J ^ I P s ^{+ [} P s + 1 ^- P s ft) ^U *^X

. 0

^{X (} P m , P s ^{+ [} P s + 1 ^- P s ^] Y ^)d Y ⁺ J ^' ( P s - 1 ⁺

⁺ ( ^[ P s ^- P s - 1 ] y) и ‘ ⁽ P m , P s - 1 ^{+ [} P s ^- P s -Ж ]"

^{- U} s [ J 0 ™ s '( P s ^{+ [} P s + 1 ^- P s Ю^Х

_X ^{d U} * ⁽ P m , P s ^{+ [} P s + 1 ^- P s ] Y ) _d + d n '              ^Y

⁺ J ₀ ¹ ®^Г ( P s - 1 ^{+ [} P s - P s - 1 ^] Y ) ^X

_X ^d u ‘ ⁽ P m , P s - 1 ^{+ [} P s ^- P s - 1 1 Y ) d d n '                ^Y

m = [ 1, M ] .

[ D ] ^+[ G ] V _Г = й Г ” .

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

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

d m = - ^{h [} J 0 ® s "( P s ^{+ [} P s + 1 - P s ^{] 7} ) ^Х

c ( L ) и ( L ) = - J f у ( п ) и * ( L , р )d V +

V

+ J q ( п ) и * ( L , n )d Г - J и ( п ) q * ( L , п ) d r .

Г                      Г

_X ^{d U} ' ⁽ P m , P s ^{+ [} P s + 1 - P s ^] Y ) d + d n '              ^Y

+ J 0 ® s "(p s - 1 ^{+ [} P s - P s - 1 ^] y )^x

_X ^{d u} ' ⁽P m = P s - 1 ^{+ [} P s - P s - 1 ^{] L )} _d d n '               ^Y

C m 5

ms ,

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

g m , s

= ^h J ^ ( p s ^{+ [} P s +1

^- P s ^] y )^x

c (L) = 1 -    Ф ,

2п

где ф - внутренней угол кусочно-линейной границы в точке L Первое слагаемое в правой части является полем, создаваемым источником f ^ в свободном пространстве, и может быть обозначено как падающее поле и in .

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

^{X u} * ⁽ P m , P s ^{+ [} P s +1 ^- P s ^] Y ^)d Y ⁺ + J ® Г ( P s -1 + [ P s - P s -1 ] y ) x 0

^{X u} ' ⁽ P m , P s -1 ^{+ [} P s ^- P s -1^] Y ^)d Y] ,

c ( L ) и ( L ) = J q ( п ) и * ( L , n )d Г -

Г

- J и ( п ) q • ( L , n )d Г + и X).

Г

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

m , s = [ 1, M ] .

,

где подматрица А Ω,Ω размерностью ( N – M) ×( N – M) включает в себя коэффициенты соотношений поля во внутренних узлах сети разбиения. Подматрицы А _Ω,Г и А _Г,Ω размерностью ( N – M) × M и M ×( N – M), соответственно, включают в себя коэффициенты связи поля в граничных и внутренних узлах. Подматрица А Г, Г размером M × M включает в себя коэффициенты связи поля в граничных узлах. Подматрица В размером M × M включает в себя коэффициенты соотношений между производными поля на границе и полем во внутренних узлах сети разбиения. Подматрица D размером M × M включает в себя коэффициенты связи поля свободного пространства в граничных узлах. Подматрица G размером M × M включает в себя коэффициенты соотношений между производными поля на границе и полем свободного пространства во внутренних узлах сети разбиения. Подматрица С Ω,Ω размерностью ( N – M) ×( N – M) включает в себя коэффициенты соотношений источников поля во внутренних узлах сети разбиения. Подматрицы С Ω,Г и С Г,Ω размерностью ( N – M) × M и M ×( N – M), соответственно, включают в себя коэффициенты соотношений источников поля в граничных и внутренних узлах. Подматрица С Г, Г размером M × M включает в себя коэффициенты соотношений источников поля в граничных узлах. Подматрица E – единичная матрица размером M × M. Векторы u Ω и u Г – вектора напряженности поля во внутренних и граничных узлах сети, v _Г – векторы нормальных производных поля в граничных узлах сети. Векторы f _Ω и f _Г – векторы источников поля во внутренних и граничных узлах сети, u ⁱ _Γ ⁿ – вектор напряженности внешнего падающего поля в граничных узлах сети. Таким образом, размерность системы уравнений (27) – ( N + M )×( N + M ).

После определения значений поля и его производных на границе Г поле в любой точке области Ψ определяется соотношением (22), где с ( ξ )=1.

Результаты численного моделирования

Рассмотрим дифракцию плоской волны (TE- и ТМ-поляризаций) на диэлектрическом и проводящем однородных цилиндрах с круглым сечением для экспериментального исследования сходимости объединенного метода. Сходимость алгоритма зависит от длины сегмента h , длины волны источника λ , µ _r и ε _r среды. Поскольку магнитная и диэлектрическая проницаемости среды являются переменными в задаче, рассмотрим зависимость решения от параметра λ/ h , определяющего количество сегментов сети на одной длине волны. Пусть плоская волна падает на цилиндр, длина волны λ 0 =0,5 мкм. Радиус цилиндра равен 0,25 мкм. Относительная диэлектрическая проницаемость проводящего цилиндра из алюминия ε= – 41,4+ i 11,9. Относительная диэлектрическая проницаемость диэлектрического цилиндра ε=2,25. Окружающее цилиндр однородное пространство имеет параметры ε = µ =1. Расчет производился объединенным методом МКЭГ-МГЭ. В качестве области Ω выбиралась квадратная область, контур Г – периметр области Ω . Источники внутри

Ω отсутствуют. Область Ω была покрыта квадратной сетью 105×105 узлов. Процесс расчета занял 10 минут на персональном компьютере Pentium4.

На рис. 1 представлены результаты численного моделирования дифракции TE- и ТМ-волны на диэлектрическом цилиндре. На рис. 2 – результаты моделирования дифракции на проводящем цилиндре.

Рис. 1. Распределение интенсивности поля дифракции на диэлектрическом цилиндре (инвертировано): ТЕ-поляризация (а), ТМ-поляризация (б)

Рис. 2. Распределение интенсивности поля дифракции на проводящем цилиндре (инвертировано): ТЕ-поляризация (а), ТМ-поляризация (б)

Для оценки дифракционных процессов используем диаграмму направленности рассеяния, зависящей от угловой координаты ϕ , определенную в бесконечно удаленных точках ( ρ→∞ ) как

σ ( ϕ ) = lim2 πρ ρ→∞

I _usc I ¹

I uⁱⁿ I 2.

В частности, диаграмма направленности определялась в прямом направлении ( ϕ =0), в обратном направлении ( ϕ = π ) и в поперечном направлении ( ϕ = π /2) при ρ =10000λ. Ниже приведены зависимости σ / λ от параметра λ/ h . Результаты представлены в дБ ( σ дБ =10 log 10 σ ). На рис. 3 представлены результаты оценки дифракции TE- и ТМ-волны на диэлектрическом цилиндре.

На рис. 4 представлены зависимости значений диаграммы направленности дифракции TE- и ТМ-волны на проводящем цилиндре.

На рис. 5 приведена зависимость от параметра λ/ h относительного отклонения значений диаграммы направленности рассеяния TE- и ТМ-волны на диэлектрическом цилиндре от значений при λ h = 100 для углов ϕ равных 0, π /2, π . На рис. 6 приведена зависимость от параметра λ/ h относительного отклонения значений диаграммы направленности рассеяния TE- и ТМ-волны на проводящем цилиндре.

Рис. 4. Зависимость диаграммы направленности рассеяния на проводящем цилиндре от параметра Ж для ТЕ-поляризации (а) и ТМ-поляризации (б)

Дст/а,% 0,30

0,25

0,20

0,15

0,10

0,05

а)

20     40     60

Рис. 6. Зависимость относительного отклонения значений диаграммы направленности рассеяния на проводящем цилиндре от параметра Ж для ТЕ-поляризации (а) и ТМ-поляризации (б)

Рис. 5. Зависимость относительного отклонения значений диаграммы направленности рассеяния на диэлектрическом цилиндре от параметра Ж для ТЕ-поляризации (а) и ТМ-поляризации (б)

Серия экспериментов с диэлектрическим цилиндром показала, что относительное отклонение значений диаграммы направленности рассеяния становится менее 5 % при λ/h>40 и менее 1% при λ/h>80 для обеих поляризаций. Эксперименты с проводящим цилиндром по- казали, что относительное отклонение значений диаграммы направленности рассеяния становится менее 5% при λ/h>30 для ТЕ-поляризации и при λ/h>50 для ТM-поляризации и менее 1% при λ/h>50 для ТЕ-поляризации и при X/h > 80 для ТМ-поляризации.

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

Заключение

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

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

Проведено исследование зависимости относительной погрешности результатов моделирования от отношения длины волны падающего света к длине сегмента сети покрытия.

Работа выполнена при поддержке российско-американской программы «Фундаментальные исследования и высшее образование» (грант CRDF RUX0-014-SA-06), а также грантов РФФИ 05-0850298, 07-07-97601 и 07-07-97600.