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

На правах рукописи

СОТНИКОВА Наталья Юрьевна

АЛГОРИТМИЧЕСКОЕ И ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ

ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧ ЭЛЕКТРОМАГНИТНОГО РАССЕЯНИЯ НА ДИЭЛЕКТРИКЕ

Специальность 05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата технических наук

Москва 2013 005544873

005544873

Работа выполнена на кафедре «Прикладная математика» Федерального государственного бюджетного образовательного учреждения «Московский государственный технический университет радиотехники, электроники и автоматики» (МГТУ МИРЭА).

Научный руководитель: доктор физико-математических наук,

профессор, заведующий кафедрой прикладной математики МГТУ МИРЭА, заслуженный деятель науки РФ Самохин Александр Борисович Научный консультант: кандидат физико-математических наук,

доцент кафедры прикладной математики

МГТУ МИРЭА_

Куликов Сергей Павлович | Официальные оппоненты: доктор физико-математических наук,

профессор кафедры волновых процессов и систем управления Московского физико-технического института (государственного университета),

заслуженный деятель науки РФ Лукин Дмитрий Сергеевич; кандидат технических наук, доцент кафедры математического обеспечения и стандартизации информационных технологий МГТУ МИРЭА Григорьев Виктор Карлович Ведущая организация: Федеральное государственное бюджетное

учреждение науки Институт радиотехники и электроники им. В.А. Котельникова Российской академии наук

Защита состоится ЯО Ре ¡г о 2013 г. в ¡У Оо на заседании диссертационного совета Д.212.131.03 при МГТУ МИРЭА по адресу: 119454, г. Москва, пр-т Вернадского, д.78, ауд. Г-412.

С диссертацией можно ознакомиться в библиотеке МГТУ МИРЭА.

Автореферат разослан /Р НЭРсРп 2013 г. Ученый секретарь

Диссертационного совета Д.212.131.03, д.т.н., профессор

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы.

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

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

При использовании итерационных методов решения СЛАУ наблюдается значительный выигрыш по числу арифметических операций и объему оперативной памяти, требуемых для реализации алгоритмов. По сравнению с методом Гаусса время решения задачи будет в N / т раз меньше, где т - требуемое число итераций. Таким образом, решение СЛАУ размерностью N»1000 становится реально осуществимой задачей для персональных компьютеров.

Поэтому разработка и программная реализация быстросхо-дящихся итерационных алгоритмов, предназначенных для числен-

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

Цель работы.

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

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

• Исследование спектра рассматриваемых интегральных уравнений для различных классов задач.

• Алгоритмы нахождения оптимальных итерационных параметров для МОПИ.

• Разработка алгоритмов и подпрограмм, основанных на МОПИ, для решения СЛАУ, возникающих после дискретизации интегральных уравнений.

• Численные решения двумерных скалярных, а также двумерных и трехмерных векторных задач.

Научная новизна.

Научная новизна работы определяется следующими положениями:

• Численно исследованы спектральные картины интегрального оператора и их закономерности использованы для построения оптимального итерационного алгоритма.

• Определены явные формулы и алгоритмы для нахождения оптимального итерационного параметра для сходимости МОПИ.

• Исследована точность и область применения первого приближения МОПИ (оптимизированное борновское приближение), которое существенно использует информацию о спектре интегрального оператора. Показана эффективность этого метода, как в рэлеевском, так и в ближнем резонансном диапазонах длин волн.

Основные положения, выносимые на защиту.

• Комплекс программ, разработанный в среде МаШсас!, позволяющий получать решения задач электромагнитного рассеяния на диэлектрических структурах.

• Алгоритмы нахождения оптимальных итерационных параметров для МОПИ для различных практически важных случаев области локализации спектра матрицы перехода.

• Комплекс алгоритмов и программ на основе МОПИ для решения интегральных уравнений.

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

Методы исследования.

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

Практическая ценность работы.

• Исследования показали практическую возможность моделирования электромагнитного рассеяния на двумерных и трехмерных локально неоднородных объектах.

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

• Оптимизированное борновское приближение может служить предобуславливателем для нестационарных итерационных методов.

• Результаты исследований вошли в работы по проектам: «Исследование путей создания экспериментального образца аппаратно-программного комплекса для решения больших вычислительных задач (с приложениями в области электромагнетизма и гидроакустики)» (2013 г.), который проводился по Государственному контракту с Министерством образования и науки Российской

Федерации; «Математическое моделирование взаимодействия электромагнитного поля со сложными трехмерными структурами на основе сингулярных объемных интегральных уравнений» ведомственной научной программы «Развитие научного потенциала высшей школы» (2009 г.). Практическая ценность выполненных разработок подтверждена актами о внедрении.

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

Апробация работы.

Основные положения и результаты работы докладывались на:

• Международном симпозиуме «Progress In Electromagnetics Research Symposium» (Москва 2009 г.) - одном из крупнейших ежегодно проводящихся симпозиумов по применению и прогрессивному развитию теоретической и прикладной электродинамики и разнообразных ее приложений во всех частотных диапазонах; симпозиумы PIERS проводятся при поддержке Академии электромагнетизма с 1989 г.;

• Международной научно-практической конференции «Актуальные проблемы и перспективы развития радиотехнических и инфокоммуникационных систем» (Москва 2013 г.);

• Международной научно-практической конференции «Современные информационные технологии и ИТ-образование» (Москва 2011 г.);

• Конкурсе прикладных разработок и исследований в области компьютерных технологий «Компьютерный континуум: от идеи до воплощения» (Москва 2011 г.) - работа отмечена дипломом;

• Международных конференциях «Информационно-вычислительные технологии в науке» (Москва 2007, 2008, 2009 гг.);

• Научно-практических конференциях «Современные информационные технологии в управлении и образовании» (Москва 2010,2011 гг.);

• Научно-технических конференциях МГТУ МИРЭА (Москва 2007, 2008, 2009, 2011 гг.);

• Конкурсе «Лучшая научная работа студентов и молодых ученых» (Москва 2008 г.) — работа отмечена дипломом 3-й степени.

Публикации.

По материалам диссертационной работы опубликовано 10 работ, из них 1 публикация в ведущем рецензируемом научном издании, рекомендованном ВАК РФ, 1 публикация в трудах международного симпозиума «Progress in Electromagnetics Research Symposium», проводящегося при поддержке Академии электромагнетизма с 1989 г. (включен в систему цитирования Web of Science - Conference Proceedings Citation Index) - может быть причислена к публикациям в изданиях, рекомендованных ВАК РФ, 2 статьи в материалах международных конференций, 6 статей в сборниках трудов конференций.

Личный вклад соискателя.

Определены оптимальные итерационные параметры для МОПИ для различных конфигураций спектра на комплексной плоскости. Разработаны комплексы алгоритмов и программ для решения интегральных уравнений для двумерного и трехмерного случаев. Проведены численные исследования.

Структура и объём работы.

Диссертационная работа состоит из введения, четырех глав, выводов, библиографического списка из 34 наименований. Работа содержит 127 страниц, включая 60 рисунков, 1 таблицу и 2 приложения.

СОДЕРЖАНИЕ РАБОТЫ

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

В первой главе описана задача рассеяния электромагнитных волн на неоднородном диэлектрическом теле.

Пусть в свободном пространстве существует некая область Q, внутри которой магнитная проницаемость /л равна магнитной проницаемости свободного пространства ц0. Диэлектрическая прони-

цаемость е внутри области задана достаточно произвольным распределением. Вне области диэлектрическая проницаемость равняется Ео-

Задача рассеяния электромагнитных волн на неоднородном диэлектрическом теле описывается интегральными уравнениями в области <2 относительно полного электрического поля Е.

Интегральные уравнения имеют вид:

• в двумерном скалярном случае -

Ё(х) = Ё° (х) + к201 [£(у) - \]Ё(у)С(г) с!у,

2

• в двумерном и трехмерном векторных случаях —

Ё(х) = Ё\х) + ф(х) - 1]ОД + к1\[е{у) -\Щу)в(г)с1у +

е

+ у.р. | §£(у) - 1}Ё(у), ётс1)ёгас1 С(г) с1у, (1)

е

где Ё(х) - неизвестное поле;

Ё°(х) — заданное первичное поле от сторонних источников в отсутствии диэлектрика;

к0 - волновое число свободного пространства;

е - комплекснозначное распределение относительной диэлектрической проницаемости по локально неоднородному телу б;

О(г) = —\Н™(к0г) - функция Грина в двумерном пространстве, или С(г) =--функция Грина в трехмерном простран-

4-ж-г

стве, где г=|х-у| - расстояние между точкой наблюдения х и точкой истока (интегрирования) у;

и = —1/2 - для двумерного случая, или

и = -1/3 - для трехмерного случая;

v.p.\ обозначает сингулярные интегралы, то есть интегралы по площади (для двумерных задач) или объему (для трехмерных задач) Q за вычетом бесконечно малого круга (для двумерных задач) или шара (для трехмерных задач) в окрестности точки у=х.

Для решения задачи интегральное уравнение дискретизирует-ся с помощью метода конечных сумм и сводится к СЛАУ относи-

тельно значений неизвестной подынтегральной функции. На область рассеяния наносится расчетная вычислительная сетка, и далее задача уже состоит в нахождении значений неизвестной функции в узлах расчетной сетки.

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

Описаны алгоритмы и подпрограммы формирования сетки для цилиндрической и декартовой расчетных моделей.

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

Запишем исходную СЛАУ в следующем виде

Ах =(2) Здесь А - невырожденная матрица системы линейных уравнений, х - неизвестный вектор, ? - правая часть. Запишем (2) в виде, более удобном для вычисления итераций

\-Тх + ф. (3)

В (3) Т— матрица перехода, используя которую вычисляются итерации по следующим формулам

5(* +1) = к = одд (4)

Для случая Т = Е - А выражения (4) представляют собой запись метода простой итерации. Для случая Т = —И'1 (Ь + £/) (4) является записью метода Якоби. Здесь исходная матрица А представляется в виде суммы нижней треугольной матрицы Ь, диагональной матрицы £> и верхней треугольной матрицы II.

Метод простой итерации для (4) сходится к решению (2) при условии, что радиус сходимости матрицы перехода (максимум модулей собственных чисел) меньше единицы, то есть

р = р{Т)< 1, (5)

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

введения в матрицу перехода комплексного параметра к дает возможность существенно расширить область сходимости. При этом модифицированная матрица перехода и новая правая часть МОПИ имеют вид

г г г \ Т-Ы ~ Ф (ел

Т = Т{к) =-, (р = ——. (6)

\-к \-к

Нашей целью является нахождение такого итерационного параметра к0 для новой матрицы перехода Т, что

р(Т(к0)) = птр{Т{к))<\. (7)

1с

Известен путь нахождения оптимального параметра сходимости МОПИ в (7). Если выпуклая оболочка комплексного спектра исходной матрицы перехода Т не содержит точку 1, то можно построить сходящийся ряд вида (4) с матрицей перехода вида (6). Оптимальным параметром сходимости лг0 является центр круга о)0, который содержит выпуклую оболочку спектра матрицы перехода Т и «виден» из точки 1 на комплексной плоскости под наименьшим углом 2 • а. При этом скорость сходимости ряда оптимальной простой итерации определяется как р - эт( а).

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

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

Если спектр матрицы перехода располагается внутри треугольника, не содержащего точку 1, алгоритм нахождения опти-

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

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

В четвертой главе даны численные результаты на основе разработанных алгоритмов и программ.

Сходимость метода простой итерации зависит от спектральных свойств матрицы перехода. Проводятся численные исследования свойств матрицы перехода. При больших значений N спектр матрицы стремится к спектру интегрального оператора.

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

В векторных случаях задачи рассеяния электромагнитных волн, а именно в двумерном случае Н-поляризации и общем трехмерном случае, спектр интегрального оператора перехода Т имеет 2 составляющие: 1) непрерывную - зависит от диэлектрической проницаемости среды (для однородной среды с диэлектрической проницаемостью в непрерывный спектр есть отрезок от точки 0 до

точки 1-е); 2) дискретную - с точками накопления на непрерывном спектре.

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

к = (\-е)/2, (8)

где £ - относительная диэлектическая проницаемость среды.

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

В случае однородной среды с потерями отрезок непрерывного спектра комплексный - [0,1 - е\ и оптимальный параметр определяется с помощью формул

с=1?: (9)

Здесь с - середина комплексного отрезка (так как рассматривается диэлектрик с потерями Ьг^б") < 0, то Яе(с) < 0, 1т(с)>0 и область расположения непрерывного спектра - вторая четверть комплексной плоскости), К - радиус окружности, проведенной через концы отрезка и точку сингулярности - точку 1, к- оптимальный параметр МОПИ.

В общем случае, для определения спектральной картины требуются значительные машинные затраты. Однако, проведенные исследования показали, что спектральная область и оптимальный параметр с хорошей точностью могут быть найдены на редкой сетке при еще недостаточной аппроксимации исходного интегрального уравнения для получения решения. В этом случае вычислительные затраты не столь велики. Другими словами, итерационный параметр находится на редкой сетке (2-3 узловых точки на длину волны), а решение интегральных уравнений с использованием найденного параметра ищется на более густой сетке, которая необходима для аппроксимации интегральных уравнений.

В примерах, представленных в работе, определение спектральных картин и оптимальных параметров, осуществлялось для значений размерности матрицы СЛАУ N порядка 100...200, так как при расчетах с большими значениями N изменений практически не наблюдалось.

На рис.] представлено решение скалярной задачи Е-поляризации частотой 0.9 гГц на длинном цилиндре радиуса 0.15м из мышечных тканей с диэлектрической проницаемостью £■ = 58-61/ и костным включением с £■ = 14 — 3/. Следует отметить, что представленная на рис.1,а спектральная картина является характерной для случая резонансного рассеяния, когда область неоднородности соизмерима или больше длины волны в среде, при этом линейные размеры области локализации спектра гораздо больше 1, а дискретный спектр располагается очень близко к точке 1, что ведет к ухудшению сходимости процесса итераций.

а) б) в)

Рис.1. Численное решение двумерной задачи рассеяния на неоднородном цилиндре.

Оптимальный параметр к — —11- 27/ был рассчитан на редкой сетке с N = 108, что заняло примерно 5 секунд. Для сравнения на более густой сетке с ;V = 1200 время расчета оптимального параметра составило порядка 10 минут, причем расчет этого варианта находится на грани исчерпания резервов оперативной памяти компьютера. Однако, как было сказано выше, в действительности уточнения значения оптимального параметра не требуется.

Найденный на редкой сетке оптимальный параметр использовался для расчета задачи поглощения на более густой сетке (10 узловых точек на длину волны). Решение было получено за т = 464 итерации с точностью £ = 10~4 по электрическому полю. На рис. 1,6 представлена сеточная картина двумерной области рассеяния, на рис. 1 ,в - диаграмма рассеяния.

На рис.2 представлено решение трехмерной задачи рассеяния плоской электромагнитной волны на диске радиуса к0с1 = 0.75Я и диэлектрической проницаемостью е = 4. Спектральная картина и оптимальный параметр к = —1-1.7/, представленные на рис.2,а, рассчитывались на редкой сетке в 2 узла на длину волны в среде с числом узловых точек N - 324.

№

Яс^),^!).!*««)

а) б)

Рис.2. Рассеяние электромагнитной волны на диске.

Диаграмма рассеяния, представленная на рис.2,б, была получена при найденном на редкой сетке значении оптимального параметра. Число узловых точек при решении интегрального уравнения составило N = 1296, число итераций т = 53.

Для оценки приближенного решения интегрального уравнения рассеяния (1) исследовалась точность и область применения первого приближения МОПИ, которое существенно использует информацию о спектре интегрального оператора Т:

т - к-1

Ет = +--. (10)

1 -к \-к

Здесь Еа) - значение поля на 1-ой итерации, Ет - начальное приближение, к - комплексный параметр, оптимальное значение для которого зависит от области локализации спектра интегрального оператора перехода Т в (10). В случае к = 0 приближение (10) известно как борновское приближение (БП), которое с хорошей точностью дает решение для малых отклонений диэлектрической проницаемости от проницаемости вакуума «1.

Проведенные численные исследования применимости оптимизированного борновского приближения (ОБП) для оценки поля в двумерных и трехмерных векторных задачах рассеяния электромагнитных волн на диэлектрическом теле показали эффективность этого метода как в рэлеевском, так и в ближне-резонансном диапазонах длин волн, а также в смешанном рэлеевско-резонансном по разным координатам диапазоне.

0.6

й 0.4

0.2

\ 1

Л

1к#

а) б)

Рис.3. Численное решение трехмерной задачи рассеяния на тонкой однородной резонансной оболочке с помощью борновского и оптимизированного борновского приближений.

Так, на рис.3,а представлен спектр тонкой и однородной резонансной оболочки с £ = 2. Радиус оболочки 2Х, толщина ЛУ60. Оболочка занимает сектор углов -л/3 < ф < л/3 в направлении распространения плоской волны. На рис.3,6 представлена симметричная диаграмма направленности. Ошибка БП достаточно значительная - 35%, ошибка же ОБП всего 1% и диаграмма направлен-

ности ОБП сливается с диаграммой направленности решения, полученного на основе интегрального уравнения и МОПИ.

В заключении сформулированы основные результаты работы и выводы.

В приложениях представлены копии актов о внедрении научных результатов диссертационной работы в научно-исследовательских работах, проводимых в МГТУ МИРЭА и ЗАО НВК «ВИСТ».

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ И ВЫВОДЫ

1. Разработан пакет алгоритмов и программ для численного решения интегрального уравнения рассеяния в двумерном и трехмерном случаях.

2. Проведено тестирование разработанного пакета алгоритмов и программ.

3. Численно исследованы спектральные картины интегрального оператора и их закономерности использованы для построения оптимального итерационного алгоритма.

4. Определены явные формулы и алгоритмы для оптимального итерационного параметра МОПИ, показана высокая эффективность этого метода.

5. Показана практическая возможность численного моделирования электромагнитного рассеяния на двумерных и трехмерных локально неоднородных объектах.

СПИСОК ПУБЛИКАЦИЙ ПО ТЕМЕ ДИССЕРТАЦИИ

В изданиях, рекомендованных ВАК РФ

1. Куликов С.П., Сотникова Н.Ю. Оптимизированное борновское приближение для векторного электромагнитного рассеяния на диэлектрике. — Электромагнитные волны и электронные системы. 2012. Т. 17. №3. С. 11-17.

В трудах международного симпозиума «Progress in Electromagnetics Research Symposium», проводящегося при поддержке

Академии электромагнетизма с 1989 г. (включен в систему

цитирования Web of Science - Conference Proceedings Citation

Index)

2. S.P. Kulikov, N.Y. Voronina. Numerical Analysis of Scattering and Absorption Problems of Electromagnetic Waves of a Mobile Communication Range on Non-uniform Biological Structures. -Progress In Electromagnetics Research Symposium Proceedings, 2009, pp. 1905-1909.

В других изданиях

3. Сотникова Н.Ю. Численное решение задач электромагнитного рассеяния методом оптимальной простой итерации. - Актуальные проблемы и перспективы развития радиотехнических и инфоком-муникационных систем. Сборник научных трудов международной научно-практической конференции. - М.: МГТУ МИРЭА, 2013. 4.2. С.285-288.

4. Куликов С.П., Сотникова Н.Ю. Оптимизированное борновское приближение для оценки электромагнитного поля дифракции на диэлектрике. - Современные информационные технологии и ИТ-образование. Сборник избранных докладов научно-практической конференции. Под ред. проф. В.А. Сухомлина. - М.: ИНТУИТ.РУ, 2011. С.811-823.

5. Куликов С.П., Сотникова Н.Ю. Численное решение краевой задачи для уравнения Пуассона с использованием метода релаксации и других итерационных методов. - Современные информационные технологии в управлении и образовании: Сборник научных трудов. -М: ФГУП НИИ «Восход», 2011. 4.2. С. 103-109.

6. Воронина Н.Ю., Куликов С.П. Численное исследование задачи поглощения электромагнитных волн диапазона мобильной связи на неоднородных биологических тканях. — Современные информационные технологии в управлении и образовании. Сборник научных трудов. М: ФГУП НИИ «Восход», 2010. 4.2. С.47-55.

7. Куликов С.П., Сотникова Н.Ю. Численное решение краевой задачи для уравнения Пуассона с использованием метода сопряженных градиентов и других итерационных методов. - 60 Научно-техническая конференция. Сборник трудов. 4.2. Физико-математические науки./ МИРЭА. - М: 2011. С.54-58.

8. Воронина Н.Ю., Куликов С.П. Алгоритмическое и программное обеспечение для численного решения интегрального уравнения рассеяния в трехмерном случае. - 58 Научно-техническая конференция. Сборник трудов. 4.2. Физико-математические науки./ МИРЭА. - М: 2009. С.35-41.

9. Волошанюк A.B., Воронина Н.Ю., Куликов С.П., Сафиули-на A.A. Использование алгоритмов поиска оптимального значения комплексного параметра модифицированного метода простой итерации для численного решения интегрального уравнения рассеяния. - 57 Научно-техническая конференция. Сборник трудов. 4.2. Физико-математические науки. / МИРЭА. - М: 2008. С.60-66.

10. Волошанюк A.B., Воронина Н.Ю., Куликов С.П., Сафиули-на A.A. Алгоритмы поиска оптимального значения комплексного параметра в модифицированном методе простой итерации. — 56 Научно-техническая конференция, посвященная 60-летию МИРЭА. Сборник трудов. 4.2. Физико-математические науки. / МИРЭА.-М: 2007. С. 11-17.

Личный вклад автора. В работах [1], [2], [4], [8] имеет место неразделимое соавторство. В [5] и [7] автору принадлежат численные результаты решения краевой задачи для уравнения Пуассона. В [6] автору принадлежит проведение численных расчетов для задачи поглощения электромагнитных волн диапазона мобильной связи на неоднородных биологических тканях. В [9] и [10] автору принадлежит разработка алгоритмов поиска оптимального параметра для метода оптимальной простой итерации.

Подписано в печать: 9.11.13 Тираж: 120 экз. Заказ № 1642

Типография ООСГАй-клуб"(Печатныйсалон МДМ) 119146,г.Москва,Комсомольский пр-т, д.28 Тел.8(495)782-88-39

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ «МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ РАДИОТЕХНИКИ, ЭЛЕКТРОНИКИ И АВТОМАТИКИ»

04201451795 На правах рукописи

Сотникова Наталья Юрьевна

АЛГОРИТМИЧЕСКОЕ И ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧ ЭЛЕКТРОМАГНИТНОГО

РАССЕЯНИЯ НА ДИЭЛЕКТРИКЕ

Специальность 05.13.18 -Математическое моделирование, численные методы и комплексы программ

Диссертация на соискание ученой степени кандидата технических наук

Научный руководитель: д.ф.-м.н., профессор, заслуженный деятель науки РФ Самохин Александр Борисович

Научный консультант: к.ф.-м.н., доцент

Куликов Сергей Павлович

Москва 2013

ОГЛАВЛЕНИЕ

Введение...............................................................................................................................3

1. Задача электромагнитного рассеяния на диэлектрике............................................9

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

1.2. Расчетная вычислительная сетка.............................................................................13

2. Итерационные методы для решения задач электромагнитного рассеяния.........21

2.1. Решение задач методом простой итерации............................................................23

2.2. Обобщение метода простой итерации.....................................................................30

2.3. Оптимальный итерационный параметр для сходимости МОПИ.........................38

2.4. Метод релаксации.....................................................................................................44

3. Разработка алгоритмов и программ для решения интегральных уравнений задач рассеяния волн с помощью МОПИ..................................................................................47

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

3.2. Описание подпрограмм определения оптимального параметра.........................50

3.3. Описание подпрограммы МОПИ для решения систем линейных уравнений.... 53

4. Численные исследования..........................................................................................55

4.1. Двумерные скалярные задачи рассеяния................................................................56

4.2. Двумерные векторные задачи рассеяния................................................................87

4.3. Трехмерные векторные задачи рассеяния...............................................................95

4.4. Исследование задач рассеяния на диэлектрике с помощью борновского приближения....................................................................................................................101

4.5. Двумерные краевые задачи электростатики.........................................................114

Заключение.......................................................................................................................119

Литература........................................................................................................................120

Приложение 1...................................................................................................................125

Приложение 2...................................................................................................................127

Введение

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

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

Указанные задачи изучались в работах отечественных ученых - А.Б. Самохина, С.П. Куликова, Ю.Г. Смирнова, В.И. Дмитриева, A.C. Ильинского, Е.В. Захарова [1-7] и других отечественных и зарубежных специалистов. В данной работе мы основываемся на полученных ими результатах.

При использовании итерационных методов решения СЛАУ наблюдается значительный выигрыш по числу арифметических операций и объему оперативной памяти, требуемых для реализации алгоритмов. По сравнению с методом Гаусса время решения задачи будет в N / т раз меньше, где т — требуемое число итераций. Таким образом, решение СЛАУ размерностью И» 1000 становится реально осуществимой задачей для персональных компьютеров.

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

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

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

• Исследование спектра рассматриваемых интегральных уравнений для различных классов задач.

• Алгоритмы нахождения оптимальных итерационных параметров для МОПИ.

• Разработка алгоритмов и подпрограмм, основанных на МОПИ, для решения СЛАУ, возникающих после дискретизации интегральных уравнений.

• Численные решения двумерных скалярных, а также двумерных и трехмерных векторных задач.

Научная новизна работы определяется следующими положениями:

• Численно исследованы спектральные картины интегрального оператора и их закономерности использованы для построения оптимального итерационного алгоритма.

• Определены явные формулы и алгоритмы для нахождения оптимального итерационного параметра для сходимости МОПИ.

• Исследована точность и область применения первого приближения МОПИ (оптимизированное борновское приближение), которое существенно использует информацию о спектре интегрального оператора. Показана эффективность этого метода, как в рэлеевском, так и в ближнем резонансном диапазонах длин волн.

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

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

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

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

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

Основные положения, выносимые на защиту:

• Комплекс программ, разработанный в среде МаШсас!, позволяющий получать решения задач электромагнитного рассеяния на диэлектрических структурах.

• Алгоритмы нахождения оптимальных итерационных параметров для МОПИ для различных практически важных случаев области локализации спектра матрицы перехода.

• Комплекс алгоритмов и программ на основе МОПИ для решения интегральных уравнений.

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

Практическая ценность работы заключается в следующих положениях:

• Исследования показали практическую возможность моделирования электромагнитного рассеяния на двумерных и трехмерных локально неоднородных объектах.

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

• Оптимизированное борновское приближение может служить предобу-славливателем для нестационарных итерационных методов.

• Результаты исследований вошли в работы по проектам: «Исследование путей создания экспериментального образца аппаратно-программного комплекса для решения больших вычислительных задач (с приложениями в области электромагнетизма и гидроакустики)» (2013 г.), который проводился по Государственному контракту с Министерством образования и науки Российской Федерации; «Математическое моделирование взаимодействия электромагнитного поля со сложными трехмерными структурами на основе сингулярных объемных интегральных уравнений» ведомственной научной программы «Развитие научного потенциала высшей школы» (2009 г.). Практическая ценность выполненных разработок подтверждена актами о внедрении, представленными в приложениях 1 и 2.

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

Апробация работы. Основные положения и результаты работы докладывались на:

• Международном симпозиуме «Progress In Electromagnetics Research Symposium» (Москва 2009 г.) - одном из крупнейших ежегодно проводящихся симпозиумов по применению и прогрессивному развитию теоретической и прикладной электродинамики и разнообразных ее приложений во всех частотных диапазонах; симпозиумы PIERS проводятся при поддержке Академии электромагнетизма с 1989 г.;

• Международной научно-практической конференции «Актуальные проблемы и перспективы развития радиотехнических и инфокоммуникацион-ных систем» (Москва 2013 г.);

• Международной научно-практической конференции «Современные информационные технологии и ИТ-образование» (Москва 2011 г.);

• Конкурсе прикладных разработок и исследований в области компьютерных технологий «Компьютерный континуум: от идеи до воплощения» (Москва 2011 г.) - работа отмечена дипломом;

• Международных конференциях «Информационно-вычислительные технологии в науке» (Москва 2007, 2008, 2009 гг.);

• Научно-практических конференциях «Современные информационные технологии в управлении и образовании» (Москва 2010, 2011 гг.);

• Научно-технических конференциях МГТУ МИРЭА (Москва 2007, 2008, 2009, 2011 гг.);

• Конкурсе «Лучшая научная работа студентов и молодых ученых» (Москва 2008 г.) - работа отмечена дипломом 3-й степени.

1. Задача электромагнитного рассеяния на диэлектрике

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

Пусть в свободном пространстве существует некая область Q, внутри которой магнитная проницаемость ¡л равна магнитной проницаемости свободного пространства /л0. Диэлектрическая проницаемость е внутри области задана достаточно произвольным распределением. Вне области диэлектрическая проницаемость постоянна и равняется е0.

Задача рассеяния электромагнитных волн на неоднородном диэлектрическом теле описывается интегральными уравнениями в области относительно полного электрического поля Е.

Интегральные уравнения электромагнитного рассеяния имеют вид:

• в двумерном скалярном случае -

Ё(х) = Ё\х) + кЩ [е(у) -1 }Ё{у)С{г) аУ,

е

• в двумерном и трехмерном векторных случаях -

Ё(х) = Ё\х) + и[е(х) -1]£(» + кЦ[Б{у) -1 ]Ё(у)С(г)с1у +

е

+ ър. I (ИдО -1 Щу), ра^роЛ 0{г) йу, (1.1.1)

е

где Ё{х) - неизвестное поле;

Ё° (х) — заданное первичное поле от сторонних источников в отсутствии диэлектрика;

б - комплекснозначное распределение относительной диэлектрической проницаемости по локально неоднородному телу <2;

(/(г) = Н^2\к0г) - функция Грина в двумерном пространстве, или

(7(г) =--функция Грина в трехмерном пространстве, где г=|х-у| -

4 • я • г

расстояние между точкой наблюдения х и точкой истока (интегрирования) у; и = -1/2 - для двумерного векторного случая, или V — -1/3 - для трехмерного векторного случая;

обозначает сингулярные интегралы, то есть интегралы по площади

(для двумерных задач) или объему (для трехмерных задач) за вычетом бесконечно малого круга (для двумерных задач) или шара (для трехмерных задач) в окрестности точки у=х.

Приведем уравнение (1.1.1) в трехмерном случае к виду, удобному для дальнейших вычислений, с помощью введения тензорной функции Грина:

т = -\{е{р)-№{рУ\-о{к5))+ / к\е{Ч)-\)ТШ¥1] , (1.1.2)

где Г = (1 + Г2УУ)С = в((\-1(кЯу1 ~(кК)-2)1 + (3(кЯу2 +31(кКу] -1)М). (1.1.3) Здесь С-С{К) - функция Грина в трехмерном пространстве, I -

единичныи тензор,

R =

rp-rq

, R = R/R, RR - тензор второго ранга, элементы

которого в декартовой системе координат имеют вид:

RR = -V R

if \2 > (Лг) АхАу Ах Az

AyAx (Ay) AyAz AzAx AzAy (Az)

(1.1.4)

у

Здесь Ах = xp-xq, Ay = yp-yq, Az = zp-zq\ для цилиндрической сетки

Ах = xm = rm cos ((pPiq), Ау = ум = fp q Sin ((pp<q).

Интеграл в смысле главного значения заменяется интегралом по области VjQ.p§, где Q.p S — шар исчезающе малого радиуса 8 вокруг точки наблюдения р. Известно, что

lim f FEdva = 0 ,

*-*°а /а

где Q.pr - шар любого радиуса г >5 вокруг точки р. Поэтому в (1.1.2) при

интегрировании по ячейке сетки при q = р остается взять численно интеграл по области между объемом ячейки и объемом шара, вписанного в нее.

Аналогичное (1.1.2) выражение получается и в двумерном случае Н-поляризации. Отличие касается первого слагаемого, коэффициент при котором и = -1/2. Тензорная функция Грина в этом случае имеет вид

где G = G(R) - функция Грина в двумерном пространстве; Hq2\JcR), H[2\kR)

А А

- функции Ханкеля 2-го рода нулевого и первого порядка; RR - тензор второго ранга с четырьмя элементами, аналогичными (1.1.4).

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

элементов ац матрицы Т является -^-(£„-1), и это значение не изменяется с

измельчением сетки.

Рассеянное поле в дальней зоне определяется на основе (1.1.2), причем в (1.1.2) первое слагаемое отсутствует ввиду отсутствия особенности в

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

E^=\k2(s(q)-l)TEdVq (1.1.5)

v

Тензор Г = (I + r2VV)G представлен в (1.1.3), (1.1.4).

Поэтому, окончательно в трехмерном случае для поля в дальней зоне получим выражение

Е

scat

1

ATTR

рх\k2(e(q)-W(q)exp(ik(v,rq)dVq

Здесь Е - поле уже в основной системе координат, связанной с телом, и имеет, вообще говоря, все три отличные от нуля компоненты, а поле Е5Са' имеет угловые компоненты в сферической системе координат.

Аналогично, в двумерном случае Н-поляризации в дальней зоне:

Е.

scat

= ^ yj-^jj PxJ*2 (e(q) ~ 1)Е (q) exp (ikrq cos (q> - tpq )dSq

Диаграмма направленности (рассеяния) в двумерном и трехмерном случаях определяется соответственно как

ДН(^) = lim IkR

Я-> 00

Е

scat

2 »

ДЩв,<р)=\т AkR2

R-* оо

Е

. scat

Е"

Если источник поля - падающая плоская волна, то в векторном двумерном случае получаем

1

ДВД =

4 к

I х J к2 (e(q) - l)E(^) exp(ikrq cos(<p - (pq )dSq

(1.1.6)

а в трехмерном случае

ДН (в,<р) =

Алк

р х J кг (s(q) - l)E(g) ехр01(р, г, )dVq

(1.1.7)

1.2. Расчетная вычислительная сетка

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

Для расчетов распределения поглощения и рассеяния электромагнитного поля были разработаны цилиндрическая и декартовая расчетные модели. Сетка адаптирована к геометрической форме тела: для округлых и цилиндрических тел используется цилиндрическая система координат, для прямоугольных и кубических тел - декартовая (рис. 1.2.1).

Рис. 1.2.1. Цилиндрическая и декартовая системы координат.

Остановимся более подробно на цилиндрической расчетной модели. В этом случае для расчетов была принята равномерная цилиндрическая расчетная сетка, со