Быстрый метод расчета дифракции на цилиндрических диэлектрических микрообъектах Текст научной статьи по специальности «Физика»

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

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

Аннотация

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

Введение

Расчет поля дифракции электромагнитной волны на диэлектрических объектах можно проводить с помощью многих известных методов: разностных методов решения системы уравнений Максвелла [13], волнового уравнения или уравнения Гельмгольца [4], методов конечных и граничных элементов [5-8], прямых методов решения соответствующих интегральных уравнений Фредгольма второго рода [9-10]. Все перечисленные методы сводят двумерную задачу дифракции к решению алгебраической системы уравнений, размерность которой равна NxN, где N -число отсчетов рассчитываемого поля дифракции. Если выбрать поле дифракции размером 128x128 отсчетов, то размерность системы уравнений будет

равна 1282 х 1282. Даже если матрица такой системы уравнений имеет, как правило, 3-х диагональный вид, время решения такой задачи на компьютерах типа Pentium IV составит десятки и сотни минут.

В данной работе предлагается решать двумерную задачу анализа поля дифракции с помощью итеративного метода с использованием алгоритма быстрого преобразования Фурье. Это всегда возможно в тех случаях, когда функция амплитуды искомого светового поля выражена с помощью интегрального преобразования, ядром которого является функция Грина данной задачи. По самому определению функции Грина интегральное преобразование имеет вид свертки, для быстрого вычисления которой можно использовать алгоритм быстрого преобразования Фурье.

В работе описан итеративный алгоритм приближенного решения интегрального уравнения Фред-гольма 2-го рода, к которому сводятся задачи дифракции электромагнитной волны TE-поляризации на двумерном (цилиндрическом) диэлектрическом объекте.

Оптимальный параметр релаксации итеративного алгоритма выбирается подбором из достаточно широкого диапазона значений. Проведено сравнение полученного решения дифракции плоской волны на круглом цилиндре с известным аналитическим решением. Алгоритм позволяет рассчитывать поле дифракции размером 128x128 отсчетов за 120 итераций, которые на компьютере Pentium III Celeron 1000 MHz выполняются за 13 секунд.

1. Итеративный алгоритм На рис. 1 показана схема рассматриваемой зад

дачи. На цилиндрический объект ( — = 0 ) с диэлек-

dz

трической проницаемостью е ( ц = 1 - магнитная проницаемость), расположенный в вакууме ( е = 1), падает TE-поляризованная монохроматическая электромагнитная волна с длиной X.

Уп

8=1

Рис. 1. Схема двумерной дифракции монохроматической волны на цилиндрическом объекте с диэлектрической проницаемостью е

Тогда интегральное уравнение, связывающее поле дифракции внутри объекта площадью 8 и вне его, имеет следующий вид [9,11]:

ф(х,у) = ф0(х,у)-lk0(е 1 Цф (Ç,n)H02)ik<J(X-О2 + (y-n)2 Wn

4

(1)

где ф(х, у) = Ег (х, у) и фо(х, у) = Еог (х, у) - проекции на ось г (перпендикулярна к плоскости рис. 1) векторов напряженности электромагнитного поля рассеянной и падающей электромагнитных волн,

2п

k = — - волновое число света, (Ç, n) е S - декар-

товы координаты внутри объекта, H02)(x) - функ

S

ции Ханкеля 2-го рода нулевого порядка. Если рассматривать координаты (х,у) тоже принадлежащими только области объекта (х, у) е £, то выражение (1) можно рассматривать как интегральное уравнение Фредгольма 2-го рода относительно амплитуды поля дифракции ф(х,у), считая функцию ф0(х,у) заданной.

Двумерный интеграл в уравнении (1) может быть представлен в виде свертки:

Цф(х,y)H02)(х -5,y -Л)п =

где S (5, n) =

(2) (3)

[ п) -ф(5, п)]* н02)(5, п), [1,(5, п) е £, 1р,(5, п) г £,

- функции области объекта £ (рис. 1). Свертка (2) функций ф£ (х, у) = £(5, п)-ф(5, п) и Н02)(5, п) может быть вычислена с помощью трех двумерных преобразований Фурье:

ф£ *Н02 = Р(ф£)-Р(н02)],

где

F [f ] = JJ f (х, У)е-г2п( х5+yn)dxdy

(4)

(5)

- прямое преобразование Фурье, Р- - обратное преобразование Фурье.

Итеративный алгоритм решения уравнения (1) имеет следующий вид:

ф п+1(х у) = Уф 0(х у) +(1 - У)ф п(х у) -

-гРЦфп (5, п)Н02ГЦ(х-5)2 + (у-п)2 )с!^п, (6)

где (х, у) е £; фп+1 и фп - амплитуда поля дифракции внутри области £ на (п+1) -ом и п-ом шагах итераций; у - постоянная релаксации алгоритма, которая регулирует скорость его сходимости,

ik o2(s-1)

4

Р =

. В качестве начального приближения

можно брать падающее поле:

ф:( х, у) = ф 0(х, у). (7)

Чтобы алгоритм (6) обладал свойством релаксации, оператор Ё в уравнении (6) должен быть «сжимающим» [12]:

||/п+1 - Яп+11| = |/ - 4п | 4 /п - Я 4 , (8)

где /и я - произвольные функции, а норма функции ||/|| определяется выражением:

1 = 1 Цf (х, У)2 dxdy .

(9)

С учетом уравнения (6) левую часть уравнения (8) можно записать в виде:

||fn+1 - ' Яп+1

= 1(1 - Y)(fn - Sn) - rP( fsn - Ssn) * H02)||

(10)

где /£ = £/, £ - функция описывающая область задания объекта (3). Используя неравенство треугольника и равенство Парсеваля можно получить оценку для нормы (10):

||/п+1 - Яп+М 4-У|+ИИт3Х }/п - Яп|| , (11)

где | р| тах - максимальное значение модуля функции Р(5, п), которая согласно (10) равна Фурье-образу от функции Ханкеля:

Р[/] = | I Н02)ГФ2 + у2 У2П(х5+^ёхёу .(12)

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

центре при г = д/х2 + у2 = 0 . Фурье образ функции Ханкеля будет функция, имеющая особенность на радиусе р = Х-1. Действительно, интеграл в (12) в полярных координатах может быть вычислен [13]:

Р (5, п) = 2п| Н 02 (кг) 0 (2прг )г^г =

i

П 2 [р2-/

In 2 -р !1

рХ < 1, рХ > 1,

(13)

где 3 0 (х) - функция Бесселя первого рода нулевого , 2п

порядка, к = —; г и р - полярные координаты:

Х

■-2 +n2

Поэтому, строго говоря, неравенство (8) может выполняться только при Y = 0, так как F = ж .

А ' II max

Но на практике интеграл в (12) вычисляется как двойная сумма на конечной сетке отсчетов, и вместо бесконечного значения функции Ханкеля в нуле выбирается ее конечное значение в близкой к нулю точке. Поэтому на практике |f|max < ж и неравенство (11) имеет смысл.

Чтобы выполнялось неравенство (8) требуется, чтобы в неравенстве (11) выражение в фигурных скобках было меньше единицы:

1 -У|+ИИmax < 1. (14)

Если выбрать параметр у в интервале (0,1): 0< у <1, то неравенство (14) будет выполняться при любом у, если выполняется условие:

MlF max < 1. (15)

S

со

ж

S

Ж

Это условие накладывает ограничение на диапазон параметров задачи, при которых она может быть решена интегральным методом (6). Например, если зафиксировать число отсчетов сетки и длину волны X, то неравенство (15) дает ограничение сверху на диэлектрическую проницаемость, которую может иметь объект:

4

е < 1 +-

* о И

(16)

Если X выбрать из интервала (1,2): 1<у <2 (при у >2 и у <0 неравенство (14) не выполняется ни при каких условиях), то неравенство (14) будет выполняться при условии более жестком, чем (15):

IN Fl

< - -1 < i.

Y

(17)

Таким образом, из анализа неравенства (14) следует, что итеративный алгоритм (6) будет обладать свойством релаксации (8) при выборе у из некоторого диапазона на интервале (0,1) при условии, что параметры задачи удовлетворяют неравенству (16). Если параметры задачи не удовлетворяют условию (16), то ни при каких у алгоритм (6) не будет работоспособен.

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

Для проверки работоспособности алгоритма (6) проводилось численное сравнение полей дифракции плоской ТЕ-поляризованной волны на круглом диэлектрическом цилиндре, рассчитанных методом (6) и с помощью известных аналитических формул [14].

На рис. 2 показана зависимость среднеквадратичного отклонения о интенсивности полей дифракции, рассчитанных методом (6) и с помощью ряда из цилиндрических функций [14], от числа итераций:

£(ф * (П, m) -ф Ideai (n, m)f n,m=1

(18)

Еф l(n m)

0,4 Л 1 1 ~~ 1 ~~ 1

0,3 \ Л г 1 1 ~ 1 ~~ 1

0,2 \ V j 1 _| 1— L

0,1 1 1 1 k- -1 L V Г*—----- L __1_ 1 L 1

0 20 40 60 80 100 120

Рис. 2. График зависимости среднеквадратичного отклонения о от числа итераций п

Время расчета 120 итераций поля дифракции 128x128 отсчетов на компьютере Celeron 1000 MHz составляет около 13 секунд. Значение невязки ст за 120 итераций достигло значения 0,0034 (рис. 2). Параметры расчета объекта следующие: диаметр круглого цилиндра D=1 мкм, длина волны X =1 мкм, е =2 - диэлектрическая проницаемость цилиндра. Рассчитывалось поле дифракции 2,5x2,5 мкм (или 128x128 отсчетов). Постоянная релаксации была выбрана оптимальной и равнялась у =0,35.

На рис. 3а приведены графики зависимости невязки ст для численного эксперимента с теми же параметрами, что и на рис. 2, но в зависимости от у при одинаковом значении числа итераций N=100. Видно, что в широком диапазоне 0,15< у <0,5 невязка не превосходит 1%.

0,20 0,16 0,12 0,08 0,04 0

—I.

__1_

—1-

-_1_

0,7 0,8 0,9 1,0 1,1 1,2

Рис. 3. График зависимости невязки о от параметра релаксации у для численного эксперимента с теми же параметрами, что и на

рис. 2 (а); и график зависимости числа итераций(о <1%) от диаметра цилиндра в мкм, при е =2 (б)

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

ст =

n,m=1

Z(k-\(п,т) -Фк(n,m))2

N

(19)

Еф2-1(n'm)

где к - номер итерации.

Сходимость алгоритма при дифракции волны на круглом диэлектрическом цилиндре зависит от

X

параметра —, где Б - диаметр цилиндра. На рис. 3б

показан график зависимости количества итераций (для достижения ошибки ст меньшей 1%) от диаметра цилиндра (в мкм) при длине волны X =1 мкм.

X

Алгоритм сходится успешно, если >0,9, что

справедливо при диэлектрической проницаемости е < 2. При е =2,5 график, аналогичный показанному

X

на рис. 3б, можно получить при условие: >1,28, а

при е =3 при условии — >1,8.

г>

а) I

I б)

ч . 1 1 1 1 1 1 Л! -1- ----1---- L 1 J -1-1- ---_l----1.--- J_ _L I I I I I _L _l J_ _L

VI V \ ---1----\\т~ I V Ч I | | | ----I---- U—-1---- ---1----f--- ---_l----1.---

10 20 30 40

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

На рис. 4 показаны результаты расчета поля дифракции с помощью алгоритма (6): расчет амплитуды дифракции плоской TE-поляризованной волны на цилиндрическом квадрате (а) и полукруге (в), и зависимости ошибки ст2 от числа итераций для

квадрата (б) и полукруга (г). Для квадрата решение стабилизировалось после 174 итераций ( ст 2 =0,000196, время расчета на компьютере Celeron 1000 MHz составило 17 секунд), а для полукруга - после 42 итераций ( ст 2 =0,000196, время расчета 4 секунды). Параметры эксперимента те же, что и для рис. 2: сторона квадрата (а) и диаметр полукруга (в) были равны длине волны 1 мкм.

На рис. 5 показана амплитуда поля дифракции непараксиального гауссова пучка на диэлектрическом цилиндре с круговым сечением ( е =2). Все поле было размером 10x10 мкм (или 512x512 отсчетов), диаметр круга равен длине волны 1 мкм, радиус перетяжки гауссова пучка (в центре картины на рис. 5) равен 1 мкм. На рис. 5 показаны: распределение амплитуды вблизи перетяжки (а); распределение амплитуды дифракции гауссова пучка на цилиндре с круглым сечением, расположенным в центре перетяжки (б), смещенным вниз от центра на 0,25 мкм (в), и на 1 мкм (г).

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

Амплитуда падающего непараксиального пучка рассчитывалась по формуле [15]:

EZ 0( ^ У) =■

1

0>/2л

j expJ - -Ц- + ik[ xi; + yVï-Ç2

I 2ст0 V

(20)

n,m=1

СТо =

2

n,m=1

ст

где ст0 - параметр гауссова пучка, связанный с ра-

1

диусом его перетяжки ю соотношением: ст0 = —.

кю

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

• Разработан итеративный алгоритм для быстрого расчета поля дифракции TE-поляризованной электромагнитной волны на цилиндрических диэлектрических микрообъектах.

• С помощью быстрого преобразования Фурье за 100 итераций, которые выполняются на компьютере Celeron 1000 MHz за 11 секунд для массива отсчетов 128x128, было рассчитано поле дифракции плоской волны на круглом цилиндре, отличающееся от точного решения в среднем на 0,34%.

• Численно показано, что параметр релаксации итеративного алгоритма можно варьировать в широком диапазоне 0,15< у <0,5 для достижения практически приемлемой точности решения задачи дифракции.

• Рассчитано поле дифракции непараксиального гауссова пучка на диэлектрическом круглом цилиндре, диаметр которого сравним с диаметром перетяжки и расположен вокруг этой перетяжки.

Благодарности

Работа поддержана грантом Президента РФ № НШ-1007.2003.1, грантом РФФИ № 02-01-08055 и российско-американской программой «Фундаментальные исследования и высшее образование» («BRHE»).

Литература

1. Taflove A. Computation electromagnetics: The finite-difference time-domain method // Artech House, Boston, 1995.

2. Prather D.W., Shi S. Formulation and application of the finite-difference time-domain method for the analysis of axially symmetric diffractive optical elements // J. Opt. Soe. Am. A., 1999. V.16. No. 5. P. 1131-1142.

3. Головашкин Д.Л., Сойфер В.А. Анализ прохождения электромагнитного излучения через дифракционную линзу // Автометрия, 1999. В. 6. С. 119-121.

4. Gruzdev V., Gruzdeva A. Finite-difference timedomain modeling of laser beam propagation and scattering in dielectric materials // Proceedings A SPIE, 2001. V. 4436, P. 27-38.

5. Davies J.B. Finite elements analysis of waveguides and cavities - a review // IEEE Trans. Magn., 1993. V. 29, P.1508-1583.

6. Mirotznik M., Prather D., Mait J. A hybrid finite element-boundary element method for the analysis of diffractive elements // J. Mod. Opt., 1996. V. 43. No. 7. P. 1309-1321.

7. Kotlyar V.V., Nesterenko D.V. A finite elements method on the problem of light diffraction by micro-optics // Opt. Mem. and Neur. Net., 2000. V. 9. No 3. P. 209-219.

8. Dong B, Lin J., Gu B., Yang G. Rigorous electromagnetic analysis of a microcylindrical axilens with long focal depth and high transverse resolution // J. Opt. Soc. Am. A, 2001. V. 8. No 7. P. 1465-1470.

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

10. Kotlyar V.V., Lichmanov M.A. Analysis of light diffraction by micro-optics using finite elements method // Opt. Mem. and Neur. Net., 2001. V. 10. No. 4, P. 257-265.

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

12. Василенко Г.И., Тараторин А.М. Восстановление изображений // М.: Радио и связь, 1986.

13. Прудников А.П., Бричков Ю.А., Маричев О.И. Интегралы и ряды. Специальные функции // М.: Наука, 1983.

14. Ваганов Р.Б., Каценеленбаум Б.З. Основы теории дифракции // М.: Наука, 1982. С. 272.

15. Petersson L.E., Smith G.S. Three-dimensional electromagnetic diffraction of a Gaussian beam by a perfectly conducting half-plane // J. Opt. Soc. Am. A, 2002. V.19. No.11. P. 2265-2280.