Применение ГРИД-технологий для решения нелинейного объемного сингулярного интегрального уравнения для определения эффективной диэлектрической проницаемости наноматериалов Текст научной статьи по специальности «Математика»

УДК 517.96+537.874.6

Ю. Г. Смирнов

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

Работа посвящена исследованию задачи определения эффективной диэлектрической проницаемости образцов наноматериалов произвольной геометрической формы, помещенных в прямоугольный волновод с идеально проводящими стенками. Задача сводится к решению нелинейного объемного сингулярного интегрального уравнения. Изучение интегрального уравнения опирается на результаты исследования соответствующей краевой задачи и теорему эквивалентности краевой задачи и интегрального уравнения. Доказаны теорема о существовании и единственности решений в Ь2 интегрального уравнения, сходимость численного метода Галеркина, получены результаты о гладкости решений. Предложен параллельный вычислительный алгоритм и процедура использования ГРИД-технологий для решения задачи.

Введение

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

Определение диэлектрических параметров нанокомпозитных материалов и сложных наноструктур с различной геометрией является актуальной задачей при использовании нанокомпозитных материалов и наноструктур на практике. Однако эти параметры не могут быть измерены экспериментально. Это приводит к необходимости применять методы математического моделирования и решать задачи численно с помощью компьютеров. При этом приходится решать трехмерные векторные задачи в полной электродинамической постановке. Решение таких задач является в настоящее время одной из самых актуальных проблем в электродинамике. Решение этих задач с приемлемой для практики точностью на электродинамическом уровне строгости математическими методами требует очень большого объема вычислений и часто невозможно даже на самых современных суперкомпьютерах. Особенно остро стоит проблема решения обратных электродинамических задач на сложной системе тел в резонансном диапазоне частот, возникающая при определении параметров нанокомпозитных материалов и наноструктур [1].

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

1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект 06-07-89063а).

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

От этих недостатков свободен метод объемных сингулярных интегральных уравнений [2, 3]. Здесь оператор получается эллиптическим, а интегральное уравнение решается только внутри тела (в области неоднородности). В отличие от [3], мы изучаем интегральное уравнение, опираясь в основном на результаты исследования соответствующей краевой задачи и теорему эквивалентности краевой задачи и интегрального уравнения. На этом пути удается доказать теорему о существовании и единственности решений в Ь2 интегрального уравнения, сходимость численного метода Галеркина, получить некоторые результаты о гладкости решений. После дискретизации получается конечномерная система уравнений с плотной (заполненной) матрицей. Таким образом, второй путь приводит к необходимости решать системы уравнений с плотными матрицами, но существенно меньших порядков (103-104).

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

Однако даже использование самых мощных отечественных кластеров (Т-60 в МГУ им. М. В. Ломоносова) оказывается недостаточным для ряда рассматриваемых задач. Поэтому необходимо применение ОМБ-технологий, позволяющих объединить мощности нескольких кластеров и выполнять на них распределенные вычисления [4].

1 Краевая задача дифракции для системы уравнений Максвелла

Рассмотрим следующую задачу дифракции. Пусть в декартовой системе координат Р = {х: 0 < Х1 < а, 0 < Х2 < Ь, -ж< Х3 < ^} - резонатор с идеально проводящей поверхностью дР . В резонаторе расположено объемное тело Q (Q с Р - область), характеризующееся постоянной магнитной проницаемостью ^0 и положительной 3 X 3-матрицей-функцией (тензором) диэлектрической проницаемости е(х). Компоненты е(х) являются ограниченными

функциями в области Q , ее ^), а также ё 1 е ^).

Граница дQ области Q кусочно-гладкая. Точнее, следуя [5], предположим, что для каждой точки границы Х0 е дQ существует окрестность 0

3 2 3

(в Я ) и С -диффеоморфизм этой окрестности на Я , при котором точка Х0

переходит в точку 0, а образом множества 0^Q является множество одного из следующих типов (ниже (Х1, Х2, Х3) - декартовы; (г, 0), г > 0, 0 е ^ - сфе-

3

рические координаты в R ). Либо xj > 0 ( Xq - точка гладкости границы); ли-

3

бо xj > 0, Х2 > 0 ( Xq - точка на «выходящем» ребре); либо R \ {xj > 0, Х2 > 0} ( xq - точка на «входящем» ребре); либо r > 0,0 є Q', где Q' с S2 - односвязная область с кусочно-гладкой границей dQ' ( Xq - вершина «конуса с ребрами»). В частности, если dQ' - гладкая, то xq - коническая точка; если dQ' образована дугами больших окружностей, то Xq - вершина многогранного угла. Пусть Q -ограниченная область и каждая точка x є dQ принадлежит одному из этих типов. Тогда будем говорить, что Q - область с кусочно-гладкой границей. Будем также предполагать, что тело Q не касается стенок резонатора, dQ п ЭР = 0 . В Р \ Q среда изотропна и однородна с постоянными £q (> 0), (> 0).

Требуется определить электромагнитное поле Е, H є L2 ¡oc (Р), возбуж-

—iwt

даемое в резонаторе сторонним полем с временной зависимостью вида e . Источник стороннего поля - электрический ток j0 є Lnoc (Р). В области

Р с R стандартные дифференциальные операторы grad, div, rot понимаются в смысле обобщенных функций.

Будем искать «слабые» (обобщенные) решения системы уравнений Максвелла (ниже понятие решения будет уточнено):

rotH = —юєЕ + jE,

rOtE = 7Ю^0 H. (1)

Эти решения должны удовлетворять условиям на бесконечности [6]: поля E и H при |хз| > C для достаточно больших C > 0 имеют представление (+ соответствует , - соответствует —го ):

H ]=I

ч J p

Xp)n рЄз — /y«V2n p''

(±)Л(р2)| хз|

p e

p

-7ЮЄ0 (V2Пp )x ез ^ f ^0 (V2¥p )x e3 ^

^ pe3 — /Yp2)V2 T.

(2)

V Р Р 3 1Р 2 Р у

где ур'* ¿0 -АР) , 1ту(* > 0 или 1ту() - 0, ¿ур;) > 0 и , Пр( Х1,Х2);

А Р2) , Vр(х^,Х2) ( ¿2 =ю2ео^о) - полная система собственных значений и ортонормированных в ¿2( П) собственных функций двумерного оператора Лапласа -А в прямоугольнике П :={( х1; Х£ ):0 < х1 < а,0 < Х£ < Ь} с условиями Дирихле и Неймана, соответственно, и = е\ Э/ЭХ1 + д/дх2 . Для коэффициентов разложений (2) имеют место оценки

4^, е(±)= о( Рт), Р , (3)

для всех т е N.

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

Сформулируем обобщенные краевые условия на границе дР [5]. Если и - достаточно гладкое векторное поле в Р, то через yvu, yTu будем обозначать нормальную и касательную составляющие поля и на дР. В негладком случае дадим определение для равенств yvu = 0, yTu = 0. Пусть

и е L2ioc (Р;C ). Тогда, если div и е L2ioc (Р), то yvu = 0 означает, что

(u,grad v) = -(div и, v) Vv eH 1ттр(Р). (4)

Если rot и е L2 ioc (Р), то ути = 0 означает, что

(и, rot w) = (rot и, w) Vw e L2, comp (Р): rot w e L2, comp (Р), (5)

где H 1(Р) - пространство Соболева; обозначим

ит 1дР-=Ytu, uv IдР-=Yvu.

Для E , H должны выполняться краевые условия на стенках резонатора:

Et |дР = 0 Hv |дР = 0. (6)

Если выполняются уравнения Максвелла, то второе условие в (6) следует из первого, и его можно опустить. Но если рассматривать оператор Максвелла, порождаемый левой частью (1), то надо ставить оба условия.

Для и е H1 (Р) существуют граничные значения из пространства

1/2

H (дР) в смысле теории следов. Почти везде на дР определен вектор нормали. Поэтому можно говорить о равенствах следов yvu = 0, yTu = 0, что будет равносильно этим равенствам в смысле данного выше определения.

Пусть также E0 и H0 - решения рассматриваемой краевой задачи в отсутствие неоднородного тела Q , e(x) = е01, x е Р (I - единичный тензор):

rotH0 =-7Ю£0 E0 + j'E, rotE0 = 7ЮЦ0 H0 (7)

с краевыми условиями

EtV = 0, HvV = 0. (8)

Эти решения могут быть выражены аналитически через jE с помощью введенного в п. 2 тензора Грина. Решения не обязаны удовлетворять условиям на бесконечности. Например, E0 и H 0 могут быть ТМ- или ТЕ-модой этого волновода.

Имеют место результаты о гладкости решений задач (1)-(6) и (7)-(8) при более гладких данных [2]. Сформулируем один из таких результатов.

Утверждение 1. Пусть ]'Е є Н^Іос(Р). Тогда Е0,Н0 є Н^Іос(Р). Пусть,

2 1 — ^ ^ 1 кроме того, дОє С , єєС (О), тогда сужения Е |д,Н |дє Н (0 и

Е \р\д,Н |р\дє Н 1Іос(Р \ О). Кроме того, справедливы условия сопряжения на дО:

[ Ет ]\дЄ = 0 [Нт ]\дЄ = ° =

где [ • ] означает разность следов с разных сторон д<2.

В предположениях утверждения 1 краевые условия на дР и условия сопряжения на д<2 понимаются в смысле равенства следов элементов из 1/2 1/2

Н Іос(дР) и Н (д<2). Ясно, что при первоначальных общих предположениях о тензоре £ такие условия сопряжения не имеют смысла.

2 Тензорная функция Грина прямоугольного резонатора

Построим диагональный тензор Грина Ое , компоненты которого являются фундаментальными решениями уравнения Гельмгольца в Р с ко-2 2

эффициентом ко =ю £0^0 и удовлетворяют краевым условиям первого или второго рода на дР, обеспечивающим обращение в нуль тангенциальных составляющих напряженности электрического поля на стенках волновода. Его компоненты имеют вид (см. [6])

п=0 т=1

2(1 + 8ри) , аЬу пт^(У птс)

(пп Л . (пт Л (пп Л . (пт Л -у |Хз-Уз|

х 0081—Х1 181ИI ~Ь~Х2 IС081—У1 I81П I ~^У2 |е

оЕ=22

п=1 т=0

2(1 + 80т ) .

аЬУ пт^ птс)

. (пп

х81ПI Х1 IС081 ~Ь~Х2 |8т

пт

пп

пт

-у 1008 I—у2 1.

У пт х3

-Уз\ •

ое=22

п=1 т=1

аЬУ пт^( У птс)

пп

пт

пп

пт

-Уз\

(9)

В этих выражениях упт = ратного корня выбирается так, чтобы 1т упт > 0 .

1

пп Л2 ( пт Л2 , 2

— | + \--------| - к0 , при этом ветвь квад-

а I І Ь I

Запишем Gm с выделенной особенностью при x = у :

1 JK\x-У1

Gm = Т-1-Т + Sm (x,У), x, У е Р , (10)

4- | x - у |

где функция gm е C“(QXР) [6, стр. 132]. Отсюда и в силу симметрии функций Грина Gm (x,у) = Gm (у,x) (m = 1,2,3) имеем

Утверждение 2. Тензор Грина Ge допускает представление „ 1 e^ko|x-У| „

Ge =~—:-------т I + g (x, У), x, у е Р, (11)

4- | x - у |

где матрица-функция (тензор) g е C“ (Q X Р) и g е C“ (Р X Q).

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

Отметим, что функции Грина имеют единственную особенность вида 1 eik0\x-у|

-----------и не имеют других особенностей в силу сделанного нами пред-

4- | x - у |

положения о том, что тело не касается поверхности волновода.

3 Объемное сингулярное интегральное уравнение

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

Пусть решения краевых задач (1)-(6) и (7), (8) существуют и единственны. Перепишем (1) в эквивалентной форме:

rotH = -7Ю£0E + jm, rotE = 7Ю|10H , (12)

где

Je = Jp + jE . (13)

В последнем равенстве jp = -7ю( e(x) -EqI')Ё - электрический ток поляризации.

Нетрудно проверить, что решение краевой задачи (6), (12) имеет вид

E = 7'ю^0Ae - ———grad div .Jp, H = rotAp, (14)

7Ю£

0

где

Ае = |&Е (г)'Е (у)Лу - (15)

Р

векторный потенциал электрического тока. Потенциал Ае удовлетворяет уравнению

аАе + ¿о Ае = - 1е . (16)

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

Однако формулы (14) не дают явного решения задачи (6), (12), т.к. ток

]е зависит от Е . Из соотношений (13)-(15) для поля Е следует интегро-дифференциальное уравнение

Е (х) = Е 0( х) + к21 Ое (г) <2

м

е0

-1

Ё ( у )4У +

+graddiv j Ое (г)

2

ё( У)

е0

-1

е(у)<яу, х є 2.

Кроме того,

Ё(х) = Ё0(х) + ко21 Ое (г)

+graddiv j Ое (г)

2

2

~ё( у) е0

ё( У)

е0

-1

Ё (у )<^у +

-1

Ё(у)йу, х є Р \ 2 .

(17)

(18)

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

Гё( У)

Н(х) = Н0 (х) - 7ЮЄ0ГОҐ| Ое (г)

2

-1

е0

Е(у)йу, х є Р .

(19)

Сведем полученное выше интегро-дифференциальное уравнение к объемному векторному сингулярному интегральному уравнению.

Представим функцию Грина в виде

Оё (г) = О0 (г) + О (г) + 02 (г), г =| х - у |

(20)

00(г)=

?ік0г -1

I, 01(г) =±-1, в2(г) = ^{£\£2,£3}. (21)

4кг 4кг

Применяя теорему [7] о дифференцировании интеграла с ядром, имеющим слабую особенность, придем к известному [3] представлению:

-----Ї-------------ип (у )сІу = V. р. Г

дх,1 дхп 4кг пПУ У 12 п 2

4~ип {у)йу - і81пип (х). (22)

ох^охп 4кг 3

Используя полученные соотношения, переходим от интегро-дифференциального уравнения (16) к векторному сингулярному интегральному уравнению:

Е(х) +1 — -1 Е(х) - V. р.Г Гт(х, у) ё(у) -1 Е(у)ёу -

■Г Г(х, у) ^ -1 ЕШу - Г Г2(х, у) ^ -1 ЁШу = Е 0( х). (23)

Q I ё0 ] Q V -о \

Здесь тензоры Г, Г1, Г 2 имеют вид:

Г( х, у) = кО&Е (г) + ( • ,егаа)егааао(г), Г1( х, у) = ( • ^га^гааа^ г),

(24)

(25)

(26)

Вопрос о разрешимости уравнения (23) и об эквивалентности краевой задачи дифракции и сингулярного интегрального уравнения устанавливается в следующей теореме.

Теорема 1. Пусть тело Q с кусочно-гладкой границей дQ характеризуется положительным тензором диэлектрической проницаемости ёеЬх^) иё-1 е Ьх^). Пусть Е,Н и Е0,Н0 - единственные решения краевых задач (1)-(6) и (7), (8), соответственно. Тогда существует и единственно решение Е е ¿^(Я) уравнения (23). Обратно, если Е е ¿2^) - решение интегрального уравнения (23), то формулы (13)-(15) (или (18), (19)) дают решение краевой задачи для системы уравнений Максвелла (1), удовлетворяющее условию (6).

Одним из наиболее эффективных методов численного решения интегральных уравнений является метод Галеркина.

Для уравнения Аф = /, (ср, / е X) в пространстве X метод формулируется следующим образом. Пусть конечномерные подпространства Хп с X являются линейными оболочками базисных функций: Хп = ърап{у1, •••, vn}, Рп : X — Хп - ортопроекторы. Потребуем, чтобы для Vk выполнялось условие аппроксимации

где ( • , • )x - скалярное произведение в X. Представив приближенное ре-

4 Метод Галеркина

Ух е X Нш т£ | х - х |= 0.

п——^ xеXn

(27)

Метод Галеркина записывается следующим образом: (АФп, V) X =(/, V) X, 1 = ^..^ n,

(28)

п

шение в виде фп = 2 , получим систему линейных алгебраических

к=1

уравнений для отыскания неизвестных коэффициентов Ск :

п

2 Ск (^к, V ) X = (/, V ) X, 1 = ^..^ п. (29)

(29)

к=1

Определение 1. Метод Галеркина будем называть сходящимся для оператора А, если существует число N такое, что для каждого /є 1тА

приближенное уравнение (28) имеет единственное решение фп є Хп для всех п > N, и если эти решения сходятся фп ^ ф при п к единственному решению ф уравнения Аф = /.

В этом случае имеет место квазиоптимальная оценка скорости сходимости [8]:

||фп -ф|И С іпґ 1к-ф|| •

¥єХп

(30)

Рассмотрим вопрос о сходимости метода Галеркина для уравнения (23). Сформулируем лемму.

Лемма 1 [8]. Предположим, что А : X ^ X есть ограниченный оператор, имеющий ограниченный обратный, и что проекционный метод сходится для А. Пусть В - линейный ограниченный оператор, А + В инъективен. Оператор В удовлетворяет любому из двух условий:

а) 8ир

nєN

А-1Р

^п * п

И = ч < 1

или

б) В компактен.

Тогда проекционный метод также сходится для оператора А + В. Перепишем интегральное уравнение (23) для электрического поля в виде

(I + 5 - К)Е = Е0, где операторы 5 и К определяются в соответствии с (23):

(31)

(5Е)( х) =

£ (х)

£0

-1

Е(х) - V.р| Г\(х, у)

2

£( у)

£0

-1

Е( у)оУ;

(КЕ)(х) = | Г( х, у) <2

£( у)

£0

-1

Е ( у^у + | Г 2( х, у)

2

£( у)

£0

-1

Е(у)йу • (32)

Применяя Лемму 1, получаем следующий результат.

Теорема 2. Пусть однородное уравнение (31) имеет только тривиальное решение и тензор диэлектрической проницаемости таков, что

о-|1/2

ess 8ир

хє2

3

2

I ,п=1

е1п( х)

-5-

0

1п

(

1 +

-1

(33)

и выполнено условие аппроксимации (27). Тогда уравнение (31) однозначно разрешимо для любой правой части Е0 є ¿2 (2) и метод Галеркина сходится для уравнения (31).

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

уравнения (17). Этот подход оказывается эффективным в силу более удобно-

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

е( х)

-1

обратима в Q,

£( х)

е0

л-1

-1

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

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

£ =

£( х)

е0

-I

J :=

£( х)

е0

-I

перейдем от (17) к следующему уравнению:

_ ТГ div x i (x, У) J (y)dy = E 01 (x), l = 1,2,3.

(34)

(35)

Q

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

3

2 V (x) - k0 i G1e (x, y) Jl (y)dy -

1= Q

- AJ = J (x) - £q i GeJ (y)dy - grad div i GeJ (y)dy = E0 (x). (36)

Q Q

Определим компоненты приближенного решения J :

_ N _ N _ N

J1 = 2aj}(x), J2 = 2hfhx), J3 = 2fx), (37)

i=1 i=1 i=1

где fk - базисные функции-«крышки», существенно зависящие лишь от переменной xk .

Ниже проводится построение функций f1. Будем считать, что Q - параллелепипед: Q = {x: a < x1 < 02, b < x2 < ¿2, C < x3 < ^2}, Q с P . Разобьем Q параллелепипедами

Пijm = {x : x1,i—1 < x1 < x1,i+1, x2, j < x2 < x2, j+1, x3,m < x3 < x3,m+1};

a2 - a , b2 - b - c

x1 i = a1 + —-Li, x2 j = ¿1 + 2—----- j, x3m = C1 + 2—---L m,

n n n

где i = 1,..., n -1; j, m = 1, ..., n/2-1.

Обозначив h1 :=| x1 i - x1 i- |, получим формулы для fjm :

Z’1 =

Jijm

1 1 I x1 x1,i I, xє Пу'т ;

x гП1

ijm •

0

2 3

Функции ^т, ^т, зависящие от переменных Х2 и Х3 соответственно,

определяются аналогичными соотношениями. Построенное множество базисных функций удовлетворяет требуемому условию аппроксимации в ¿2 [9].

12 3

Перенумеруем базисные функции fi ,, г = 1,..., N. Расширенную матрицу для нахождения неизвестных коэффициентов аг, Ъ, сг удобно представить в блочной форме:

A11 12 A13 B1Л

21 A A22 A23 B2

A A32 A33 B3 J

элементы колонок Bfr и матриц определяются из соотношений:

Bk = ( eQ , fk);

j GkE (x, y) f/(y)dy, fk

Q

Q

j G (x, y) -^—flj(y)dy, -^—fi

J VXk nxi

— fk dxk

\

k,l = 1,2,3; i, j = 1,..., Ж Здесь функция G имеет вид

4

G=-n

n=1 m=1

aby „mSh(Y nmc)

. f nn ^ . f Km ^ . f nn ^ . f кт Л -y X-VA x sin I—xi I sin 1-^X2 I sin I—yi I sin I —ь—y2 I e 3 .

Предложенный метод Галеркина реализован для решения ряда задач дифракции.

Часто интерес представляют задачи рассеяния в среде, характеризующейся постоянной во всем объеме резонатора диэлектрической проницаемостью (£ = £qI) и тензорной магнитной проницаемостью Д в Q (вне

Q Д = ДоI). В этом случае краевая задача сводится к объемному сингулярному интегральному уравнению (такого же типа) для магнитного поля и выражению для электрического поля через решение этого уравнения

H (X) = H Q( x) + kQ j GH (r)

Q

Д (y)

Д( y) До

- i

H (y)dy +

+graddiv j Gh (r)

Q

-1

H(y)dy, x є Q;

E(х) = E0 (х) + іощогоґ J Gh (r)

Q

Д( y) _ j

. До .

H(y)dy, x є P .

В последних формулах Gh (x, y) - тензорная функция Грина прямоугольного волновода, отвечающая произвольному распределению источников магнитного поля. Как и для рассматривавшейся функции Грина Ge (х, y), имеет место представление в виде суммы сингулярного слагаемого того же вида и гладкой функции. Следовательно, для задачи о возбуждении резонатора магнитным током верны все теоремы, сформулированные выше.

5 Обратная краевая задача

Мы будем рассматривать обратную краевую задачу для определения эффективной диэлектрической проницаемости образца наноматериала, расположенного в волноводе. Рассмотрим изотропный случай и будем считать, что є(х) = є, где є - неизвестная константа (эффективная диэлектрическая проницаемость) образца [1]. Предположим, что я/a < ko < я/b . В этом случае в волноводе может распространяться только одна мода, потому что

Im= 0, y1^ =y¡ko _—2/a2 > 0 и Im у(р') > 0 для всех p, j за исключением p = 1 и j = 2 . Мы также предполагаем, что

і?0 ¡ \ /((+)■ Я • Ях1 _іу12^х3

E (х) = ^2А 'iro^-sin— e 1 3.

a a

Здесь A(+) - (известная) амплитуда распространяющейся волны,

^1 = cos ЯХ1/a . Следовательно, gE ^ 0 и G— ^ 0 равномерно по y є Q при |х3| . Мы также получаем, что

G2 1 • Ях1 • яу1 —іу12)х3_y3\ , 0

Ge------sin—Lsin—-e 1 1 3 3 ^ 0

abY10 a a

равномерно по y є Q при |х^ ^го. Затем, мы имеем: divGe ^ 0 равномерно

dG2

по y є Q при 1х3 (потому что —— ^ 0 равномерно по y є Q при

дх2

|х3 ). Вычислив предел при |х3 в (18), получим уравнение

E(х) = E0 (х) + k02 (-^_ 1 IJG— (х,y)—2 (y)dy, хє Q, (38)

є" ) Q

є0

и, принимая во внимание условие на бесконечности (2), при |х3| -

e2Q1(+)e-iY' х3iro^0 — sin—1 = a a

-iY(2),

3 ІГОД0

a a

.(+) _іу(2)х^ Я . яхі

= e2 A 'e 1 3 ІЮД^—sin—- +

1\к0в^ Гчш - Уз)

V е0 )

І 81И—е

аЬую ' а а

Е2 (у)ф. (39)

Из этого следует

е((+) = А+) + ко2 Г — _ 1І --------------------------181И М. е"^ Уз Е2 (у )<у. (40)

V Є0 ) 1 л/'пгпи.л Л а

•( 2)У

Мы предполагаем, что коэффициент е| известен из эксперимента. Таким образом, мы имеем

С

_ 1 = -

где

С = -

Єо ( І, Е )

/люцо*Уіо ( б{+)_ А (+)

(41)

/2

ко

г ■ Яу1 г'у(2)у3

І = Є2 81П —е П ^3 .

а скобки обозначают скалярное произведение в пространстве ¿2 ( Q):

( {,Е) = |Г(.у)Е(^)у .

е

(42)

(43)

(44)

6 Итерационный метод для решения обратной краевой задачи

Подставляя (34) и (41) в формулу (35), мы получаем нелинейное объемное интегральное уравнение:

(, Е )

С

)Е(х)_Ео ( х)) = ко21Ое(х,у)Е(у)у-

Є

+ егаааіуІОе ( х,у)Е(у)у, хє Q.

е

(45)

Введем линейный интегральный оператор:

4)Е := к)1 СЕ( х, у )Е( у )с!у + graddiv | Ое( х, у )Е( у )с!у . (46)

е

е

Рассмотрим итерационный процесс для решения нелинейного интегрального уравнения (45):

С

( Еи( х)_ Ео( х )) = ( А)Еи+1 )(х).

(47)

При п = 0, 1, ... на каждом шаге приходится решать линейное объемное интегральное уравнение. Метод решения уравнения описан в п. 4. После решения уравнения (45) с заданной точностью с помощью итерационной процедуры (47) по формулам (41)-(43) находим неизвестную диэлектрическую проницаемость е.

В предложенной процедуре определения е наиболее сложным этапом является решение уравнения (47) на каждом шаге итераций. Решению объемных сингулярных интегральных уравнений посвящены работы [2-4]. Ниже описывается параллельный вычислительный алгоритм и процедура применения ГРИД-технологий для решения этих уравнений.

7 Субиерархический параллельный алгоритм

При численном решении уравнений (47) можно использовать параллельный алгоритм для многопроцессорных кластеров и распределенных вычислительных систем [1)]. Неизбежность использования подобных алгоритмов вызвана большим объемом вычислительной работы. Последнее обусловлено, прежде всего, трехмерным и векторным характером задачи, численное решение которой приводит к системам линейных алгебраических уравнений (СЛАУ) большой размерности (порядка ста тысяч и более). При решении этих систем (за приемлемое время) эффективно применение параллельных версий решения СЛАУ методом сопряженных градиентов [11]. Однако при решении задачи наибольшую трудность представляет не решение системы, а ее заполнение: элементы матрицы представляются через шестимерные интегралы от тензора Грина, компоненты которого представляются в виде рядов. Решение данной проблемы описано в [4]. Заметим, что процесс составления матрицы СЛАУ и ее решения упрощается за счет использования теплицевой структуры матрицы.

Основным вычислительным узлом в методе сопряженных градиентов, как в любом итерационном методе, является процедура умножения матрицы на вектор, к которой и применяется параллельный алгоритм. Совместно с применением параллельного алгоритма умножения матрицы на вектор, для сокращения времени счета задачи используется субиерархический подход [4]. Этот подход позволяет использовать предварительно вычисленные элементы матрицы для канонической геометрической фигуры - куба, из которого пользователем «вырезается» рассчитываемая фигура образца. В начале расчета пользователем посредством web-интерфейса производится выбор параметров счета и геометрии задачи. При каждом умножении матрицы на вектор производится заполнение нулями соответствующих векторов в зависимости от выбранной геометрии. Таким образом, задача решается для образца произвольной геометрической формы, «вырезанного» из исходного куба.

8 Применение ГРИД-технологий

Как отмечалось раннее, решение поставленной задачи может оказаться очень емким с точки зрения вычислительного процесса. Использование суб-иерархического подхода позволяет решать задачи на телах с произвольной геометрией разной вычислительной сложности: от очень простых (когда выбрано несколько десятков носителей) до очень сложных (с числом носителей порядка ста тысяч). Для быстрого решения простых задач достаточно ресурсов небольшого кластера, в то время как решение сложных задач требует ис-

пользования ресурсов самых современных кластеров. В связи с этим естественно использовать ГРИД-технологию для решения столь сложной вычислительной задачи [4]. Использование web-интерфейса позволяет пользователю задать его параметры счета и геометрию задачи, перед решением задачи на кластере производить анализ сложности запускаемой задачи. Если для решения задачи требуются небольшие вычислительные ресурсы, она решается на мини-кластере, в случае использования больших вычислительных ресурсов решение производится на более мощных кластерах. В качестве кластера для небольших вычислительных задач предполагается использовать миникластер ПГУ, для задач со средней сложностью - кластер ЮУрГУ, для самых больших задач предполагается использование кластера НИВЦ МГУ. Задачи решаются в режиме метакомпьютинга. Функциональная схема использования ресурсов при решении задачи изображена на рис. 1.

Получение результатов

Рис. 1 Использование данных кластеров и ресурсов СКИФ-ГРИД полигона производится в соответствии с проектом Союзного государства СКИФ-ГРИД, государственный контракт от 16 июля 2007 г. № СГ-2/07

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

1. Shestopalov, Yu. V. Volume singular integral equation method for determination of effective permittivity of meta-and nano-materials / Yu. V. Shestopalov, Yu. G. Smirnov, V. V. Yakovlev // Proc. Progress in Electromagnetics Research Symposium (PIERS'2008) (Jule 2-6, 2008). - Cambridge, MA, 2008. - P. 291-292.

2. Смирнов, Ю. Г. Исследование электромагнитной задачи дифракции на диэлектрическом теле методом объемного сингулярного интегрального уравнения / Ю. Г. Смирнов, А. А. Цупак // Журнал вычислительной математики и математической физики. - 2004. - Т. 44. - № 12. - С. 2264-2279.

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

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

5. Бирман М. Ш., Соломяк М. З. // Успехи математических наук. - 1987. -Т. 42. - Вып. 6. - С. 61-75.

6. Ильинский, А. С. Дифракция электромагнитных волн на проводящих экранах / А. С. Ильинский, Ю. Г. Смирнов. - М. : ИПРЖР, 1996.

7. Михлин, С. Г. Многомерные сингулярные интегралы и интегральные уравнения / С. Г. Михлин. - М. : Физматгиз, 1962.

8. Kress, R. Linear Integral Equations / R. Kress // Applied Mathematical Sciences. -Vol. 82. - Springer-Verlag. New-York Inc., 1989.

9. Марчук, Г. И. Введение в проекционно-сеточные методы / Г. И. Марчук, В. И. Агошков. - М. : Наука, 1981.

10. Воеводин, В. В. Параллельные вычисления / В. Воеводин, Вл. В. Воеводин. -СПб. : БХВ-Петербург, 2002.

11. Ортега. Параллельное программирование для многопроцессорных вычислительных систем / Ортега. - СПб. : БХВ-Петербург, 2002.