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

УДК 519.642

М. Ю. Медведик, Д. А. Миронов, Ю. Г. Смирнов

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

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

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

Abstract. Electromagnetic diffraction problem on dielectric body located in rectangular waveguide is considered. The problem is reduced to volume singular integral equation on the body. Numerical collocation method for solving the equation is considered. Taking info account complicated calculations supercomputer was used for solving the problem.

Keywords: boundary value problem, electromagnetic scattering, integral equations, numerical method.

Введение

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

1 Работа выполнена при поддержке гранта Минобрнауки РФ по ФЦП «Развитие потенциала высшей школы» № 2.1.1/1647.

Альтернативным подходом является применение метода объемных сингулярных интегральных уравнений [1]. Настоящая статья посвящена разработке численного метода для решения уравнения. Применяется метод коллокации с аналитическим суммированием медленно сходящихся рядов в функциях Грина.

1. Постановка задачи

Рассмотрим следующую задачу дифракции. Пусть в декартовой системе координат Р = {х : 0 < xi < a, 0 < x2 < b, - ^ < X3 < ^} - волновод с идеально проводящей поверхностью ЭР. В волноводе расположено объемное тело Q ( Q с Р - область), характеризующееся постоянной магнитной проницаемостью Цо и положительной 3 X 3 -матрицей-функцией (тензором) диэлектрической проницаемости ё(х). Компоненты ё(х) являются ограниченными функциями в области Q , êe L^ (Q), а также ê-1 е L(Q). Граница ЭQ области Q кусочно-гладкая. Требуется определить электромагнитное поле E, H е L2,

loe (P) , возбуждаемое в волноводе сторонним полем с временной зависимостью вида e lfflt. Источник стороннего поля - электрический ток j0 e L2,ioc(P) [2, 3].

В работах [2, 3] задача сводится к решению интегродифференциального уравнения для поля E :

E( x) = E0( x) + ko2 j Ge (r )

Q

ë( y)

e0

-1

E(y)dy +

- grad div J Ge (r )

Q

ê( y)

e0

-1

E(y )dy, x є Q,

(1)

, 2 2 где k0 =ю • e0 • |10 .

Кроме того,

E(x) = E0(x) + ko2 J Ge (r)

Q

ë( y)

e0

-1

E(y)dy +

- grad div J Ge (r )

Q

ë( y)

eo

-1

E(y)dy, x є P \ Q.

(2)

Формула (2) дает представление решения Е(х) в области Р \ Q , если Е(у), у е Q - решение уравнения (1). Поле Н выражается через решение (1) в виде

H( x) = H0( x) - г'ює0 rot J Ge (r)

Q

ë( y)

e0

-1

E(y)dy, x є P.

(З)

ЗЗ

Выражение для функции Грина имеет вид

оЕ = ^(оЕ, оЕ, оЕ),

где

^ ^ —V \Хо — Vo

e *пт\ 3 кп пт кп пт

---------cos—Xi sin-------x2 cos—yi sin-------y2 ;

abn^0 m^l Vnm(1 + 80n) a b a b

GE=abI I-

^2 ^ ^ ^ e VnmlX3 Уз\ . кп кт . кп кт

GE =—rI I-------------n , s , sin — X1 cos~rx2 sin — Vl cos^“У2 ;

ab “1 т=0 Vпт(1 + 80т) a b a b

_ ^ ^ —V X — V

3 2 ^ ^ e пт| 3 31 . кп . кт . кп . кт

Ge = — I I -sin—X1 sin — X2sin—y^m—У2 .

ab , , упт a b a b

п=1 т=1

В этих выражениях упт =

кп ^2 Г кт ^2 2

— I +1---------I - ко , при этом ветвь квад-

о ^ I Ь )

ратного корня выбирается так, чтобы 1т упт ^ 0 . После выделения особенности имеем

Go(r)=

Ge(r) = Go(r) + Goo(r) + g2(r),r =| X — y |;

’ЛоГ — 1 » . . 1 - -..12 3

. 1, Goo(r) = -—1, go(r) = d^gteo^o^ob

4кг 4кг

12 3

где £0, £0, g0 - гладкие функции.

Проинтегрируем функции Грина по параллелепипеду

Рц = |(хЬ x2, х3): »1 - - »1 + 1, »2 - ^ - г2 + 1 »3 - - »3 + ^

(4)

(5)

и обозначим их через О1, G2, Оз . Последовательно вычислим операции дивергенции и градиента:

grad div (G1J1, G2 J2, G3 J3) = grad

f 9Gi 9G;

9x1

9x19x2 9x19x3

x 1 +

9x29x1

9x22

9x29x3

x 2 +

--------U1 1------------.

9x39x1 9x39x2

9x

3

x 3,

где J = (J1, J2, J3). 34

д 2g д2G

Значения функций Gi, G2, G3 и вторых производных 1 1

д^2 ’ дx2 дx]

1

д 2G1 д 2G2 д 2G2 д2G2 д 2G3 д 2G3 д 2G

-----—, -------—, ---------------—,---------—,-----------—, --представлены в [4].

дxздxl дxlдx2 д^ дx3дx2 дxlдxз дx2дxз дх3

Нетрудно видеть, что

Gk (i1 , г2, Jl, J2) = Gk ( Jl, J2, г1, г2) ; (6)

2 2

(i1, г2, Jl, /2) = ( Jl, J2, г1, i2) ; (7)

дx^2 дx2

д2^ ч д2а

— (г1,Ь Jb J2) = Л -,2 (Jl, J2,г1,i2) , (8)

дx2дxl дxlдx2

k = 1,..,3 ;

д G- (i1, г2, Jl, J2) = ? Gl (i1, J2, jl, i2) ; (9)

дxздxl ’ ’ ’ дxздxl

дд G3 Оь^ jl,J2)=дг"G^L(il,j2, jl,i2) ; (10)

дxlдxз дxlдxз

д2g д2G

-(гЪг2, jl, J2) = - ^ (jl,^гЬ J2) ; (11)

дxздxl ’ ’ ’ дxlдxз

дд G (ib^ jl, J2) = -f V* (jl,г2,J2) ; (12)

дXздX2 дXздX2

д2G3 ґ ч д2G3 ,

- (ib г 2, J l, J 2) = Л -, ( J l, ^ J 2) ; (13)

дx2дxз дx2дxз

д2G2 . _ д2G3

Д Д (гЪ^ Jb J2) = - а (ib J2, Jbi2) . (14)

дxздx2 дx2дxз

Формулы симметрии используются далее в реализации метода колло-кации.

2. Метод коллокации

Рассмотрим теперь вопрос о построении схемы для метода коллокации рассматриваемой задачи дифракции. Будем формулировать метод для интег-

родифференциального уравнения. Предположим, что матрица

V £0

об-

£( x)

ратима в Q ,-----------1 є Lœ (Q), I - единичная матрица.

є0

Введя обозначения

e(x)

£0

х-1

-I

, J :=

ё( х)

£0

-I

E.

перейдем к следующему уравнению:

AJ = £(x)J(x) -k0 J G(x, y)J(y)dy -graddiv J G(x, y)J(y)dy = E0(x), x e ß .

ß ß

Представим это уравнение в виде системы трех скалярных уравнений:

Y^-j (х) - kl j G(x, y)Jl (y)dy-A div х j G (x, y)J (y)dy=E01 (x),l=1,2,3. (15)

3

Y

i=1

Q

Q

12 3

Определим компоненты приближенного решения Jп = (Jn, Jn , Jn) следующим образом:

Jn = £ahfl(x), J2n = £bkfk2(x), J3 = £Ckf3(x), k=1 k =1 k=l

где fk - базисные функции-«ступеньки».

Ниже проводится построение функций fk . Будем считать, что Q - параллелепипед: Q = {x: a < xi < ö2, bi < Х2 < b2, c < X3 < C2}. Разобьем Q на элементарные параллелепипеды:

Пk/от = {x : x1,k < x1 < x1,k+1, x2,l < xl < x2,l+1, x3,m < x3 < x3,m+1};

a2 - a1 7 7 b2 - b1 c2 - c1

x1,k = a1 +-k, x2, l = b1 +-l, x3,m = C1 +------------------m,

n n n

где k, l, m = 0,..., n-1.

Получим формулы для f{dm , i = 1, 2, 3 :

[1, x ^ nklm,

0, x £ Пкіт.

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

3

мому условию аппроксимации в ¿2 = ¿2 X ¿2 X ¿2.

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

' 41 A 2 A 3 B1Л

A21 A22 A23 B2

ч A31 A32 A33 B3 J

Элементы столбцов Бу и матриц Ац определяются из соотношений

4 = ¿0 (х);

4 =1и/1 ( ) - §£/£о | Ок (ху, у) /I (у)ёу -А | -^-О1 (ху, у)/} (у)ёу, (16)

д Ч д Х1

где координаты точки коллокации

Х = {х71,х72,х73 ), х/1 = ( + 1/2) х72 = ( + 1/2))ь х73 =(3 +

к, I = 1,2,3; 71,72, ¿3, У1, У2, У3 = 0, •••, п -1.

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

3. Реализация метода. Результаты расчетов

С учетом симметрии (6)-(14) общее количество коэффициентов Ау , которое необходимо вычислить, определяется по формуле

N = 3 ■п4 •(п2 + 1) + п6 + 2 ■п5 •(п +1) . (17)

2 2

Без учета симметрии это число равно N = 9п6, где п - количество интервалов разбиения всей области по одной координате.

Для получения результатов были установлены значения п = 9,10 .

Зависимость количества коэффициентов АЦ и необходимого объема

оперативной памяти для их хранения прямо пропорциональная. Результаты расчета необходимого объема оперативной памяти при установленных значениях п с учетом симметрии (6)-(14) представлены в табл. 1. Также в табл. 1

показан коэффициент уменьшения количества значений АЦ при учете симметрии по отношению к количеству вычисляемых коэффициентов без учета симметрии [4].

Таблица 1

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

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

9 2735937 41,747 1,748

10 5130000 78,278 1,754

Для вычисления необходимых коэффициентов Aj в суммах рядов

с количеством членов, равным 500, использованы ресурсы многопроцессорных вычислительных комплексов и составлена программа с использованием программного интерфейса MPI. MPI [5] - удобный стандартный API для использования в прикладных задачах ресурсов многопроцессорных комплексов. На каждом вычислительном многопроцессорном комплексе используется одна или несколько реализаций (компиляторов) MPI.

Для упрощения вычислений коэффициентов Aj и передачи их между

процессами [5] память была выделена в виде следующих шести одномерных массивов:

G

1) массив коэффициентов Gj и-----г1 ;

Эл^

Э 2g

2) массив коэффициентов G^ и-----г2 ;

Эх2

Э 2g

3) массив коэффициентов G3 и------3 ;

Эх3

4) массив коэффициентов

5) массив коэффициентов

6) массив коэффициентов

Э G

dx^dx2

э 2G3

Эх^Эхз

э2Сз

Эх2Эх3

Использована схема [6] для распределения вычислений коэффициентов каждого массива между процессами.

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

C =

Ni

Ni

1 J Ni I

+ 1, если номер процесса меньше < — f ;

I Р \

, если номер процессабольше или равен ■

где N - общее количество коэффициентов в массиве с номером і (і = 1,..., 6); { } - остаток от целочисленного деления; [ ] - целая часть деления; р - количество выделенных процессов на задачу.

Программа для вычисления коэффициентов АЦ, реализованная с использованием МР1-функций, была запущена на суперкомпьютерам комплексе СКИФ МГУ «Чебышев». Основные характеристики комплекса представлены в [6]. Количество выделенных процессоров на задачу - 100.

На примере работы процесса с номером 0, где производится максимальное количество вычислений, оценим уменьшение времени работы программы для вычисления коэффициентов матрицы при учете симметрии и без учета времени записи на диск. Исходные данные и оценка для п = 9,10 при количестве процессов равном 100 представлены в табл. 2.

Таблица 2

Количество коэффициентов для процесса с номером 0 и оценка коэффициента уменьшения времени вычисления коэффициентов при учете симметрии матрицы и разбиения на несколько массивов при количестве процессов равном 100

п 9 10

Количество коэффициентов в соответствии с распределением, использованным в [6] 47830 90000

Количество коэффициентов с учетом симметрии и с разбиением на несколько массивов 19294 36156

Оценка коэффициента уменьшения времени работы программы 2,48 2,49

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

Таблица 3

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

п 9 10

Время выполнения (с) 4673 8812

Время выполнения (мин) 77,88 146,87

Время выполнения (ч) 1,30 2,44

Коэффициент уменьшения времени выполнения программы 2,29 2,26

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

На рис. 1 представлено абсолютное значение третьей компоненты, полученное решением СЛАУ методом сопряженных градиентов [7, 8] при а = 2, Ь = 1, с = 2, к0 = 2,5, п = 10, х3 = 1,5 (сечение по координате).

Получены результаты по решению объемного сингулярного интегрального уравнения задачи дифракции на диэлектрическом теле, размещенном в слое при п = 9,10, где п - количество интервалов разбиения всей области по одной координате. По результатам табл. 3 видно, что при использовании симметрии матрицы и разделении вычислений коэффициентов по типу получен выигрыш по времени вычисления коэффициентов матрицы в 2,29 и 2,26 раза при п, равном 9 и 10 соответственно.

Рис. 1

4. Субиерархический алгоритм

Задача определения электромагнитного поля внутри прямоугольного волновода решается для заполненной секции волновода методом коллокации. Рассмотрим алгоритм расчета электромагнитного поля в прямоугольном волноводе для тела произвольной формы. Будем предполагать, что решение задачи для заполненной секции получено, и в нашем распоряжении находится матрица, составленная методом коллокации. Для решения задачи дифракции на теле сложной формы необходимо, чтобы тело целиком вмещалось в используемую нами секцию прямоугольного волновода [9, 10]. Предлагаемый метод позволяет получить матрицу для новой фигуры, пользуясь матрицей, составленной для заполненной секции волновода. Скорость построения новой матрицы будет напрямую зависеть от размера фигуры и размера сетки. На рис. 2 это рассматривается на примере одной фигуры, описываемой двумя разными сетками. Слева находится фигура, в центре приведен пример неоптимальной сетки, справа построена оптимальная сетка. С точки зрения скорости счета правая сетка предпочтительнее.

I П П I ЕЕЕ

Рис. 2

В построенной фигуре введем новую нумерацию элементарных параллелепипедов. Произведя полный перебор элементарных параллелепипедов, получаем новую сетку. Эту сетку будем использовать для расчета поля на теле сложной формы. Далее, решая СЛАУ для матрицы, составленной с использованием новой сетки, находим значения поля внутри фигуры сложной формы.

Приведем результаты расчета электрического тока на фигуре сложной формы на сетке 9x9x9. Волна падает вдоль оси 02. Правая часть равна 1. Конфигурацию фигуры (рис. 3) представим по слоям, перпендикулярным оси

02. Форма первых трех слоев представлена в левой части (рис. 3), три центральных слоя представлены в центре и форма последних трех слоев фигуры представлена в правой части. Значения поля на втором, пятом и седьмом слое соответственно (слева на право), представлены на рис. 4.

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

реализованный с использованием суперкомпьютерных вычислений.

Список литературы

1. Самохин, A. Б. Интегральные уравнения и итерационные методы в электромагнитном рассеянии / A. Б. Самохин - М. : Радио и Связь, 1998.

2. Медведик, М. Ю. Применение ГРИД-технологий для решения объемного сингулярного интегрального уравнения для задачи дифракции на диэлектрическом теле субиерархическим методом / М. Ю. Медведик, Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2008. - № 2. - С. 2-14.

3. Васюиии, Д. И. Метод коллокации решения объемного сингулярного интегрального уравнения в задаче определения диэлектрической проницаемости материала / Д. И. Васюнин, М. Ю. Медведик, Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2009. - № 3. -С. 68-78.

4. Медведик, М. Ю. Численное решение объемного сингулярного интегрального уравнения методом коллокации / М. Ю. Медведик, Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. -2009. - № 4. - С. 55-71.

5. MPI: A Message - Passing Interface Standart. Version 1.0. - University of Tennessee. -1994. - May, 5.

6. Миронов, Д. А. Применение суперкомпьютерных вычислительных сред для решения объемного сингулярного интегрального уравнения задачи дифракции на диэлектрическом теле / Д. А. Миронов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2009. - № 2. - С. 14-24.

7. Ор тега, Дж. Введение в параллельные и векторные методы решения линейных систем / Дж. Ортега. - М. : Мир, 1991.

8. Миронов, Д. А. Применение суперкомпьютерных вычислительных комплексов для решения объемного сингулярного интегрального уравнения задачи дифракции на диэлектрическом теле / Д. А. Миронов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2008. - № 3. -С. 55-62.

9. Медведик, М. Ю. Субиерархический метод решения интегрального уравнения на плоских экранах произвольной формы / М. Ю. Медведик // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. -2009. - № 4. - С. 49-55.

10. Медведик, М. Ю. Субиерархический параллельный вычислительный алгоритм для решения задач дифракции электромагнитных волн на плоских экранах / М. Ю. Медведик, Ю. Г. Смирнов // Радиотехника и электроника. - 2007. - Т. 6. -№ 4. - С. 1-6.

Медведик Михаил Юрьевич

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

Medvedik Mikhail Yuryevich Candidate of physical and mathematical sciences, associate professor, sub-department of mathematics and supercomputer modeling,

Penza State University

E-mail: [email protected]

Миронов Денис Алексеевич аспирант, Пензенский государственный университет

E-mail: [email protected]

Смирнов Юрий Геннадьевич

доктор физико-математических наук, профессор, заведующий кафедрой математики и суперкомпьютерного моделирования, Пензенский государственный университет

E-mail: [email protected]

УДК 519.б42 Медведик, М. Ю.

Субиерархический подход для решения объемного сингулярного интегрального уравнения задачи дифракции на диэлектрическом теле в волноводе методом коллокации I М. Ю. Медведик, Д. А. Миронов, Ю. Г. Смирнов II Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2Q1Q. - № 2 (14). - С. З2-4З.

Mironov Denis Alekseevich

Postgraduate student,

Penza state university

Smirnov Yury Gennadyevich Doctor of physical and mathematical sciences, professor, head of sub-department of mathematics and supercomputer modeling, Penza State University