Расчет волноводов методом конечных элементов с использованием процедуры Банча-Кауфман Текст научной статьи по специальности «Математика»

ТЕОРЕТИЧЕСКАЯ И МАТЕМАТИЧЕСКАЯ ФИЗИКА

Расчет волноводов методом конечных элементов с использованием

процедуры Банча-Кауфман

Ю. В. Мухартова, Н.А. Боголюбова

Московский государственный университет имени М. В. Ломоносова, физический факультет, кафедра математики. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2. E-mail: a russell67@yandex.ru

Статья поступила 04.02.2014, подписана в печать 06.02.2014.

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

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

УДК: 519.63; 537.87. PACS: 02.60.Cb.

Введение

Метод конечных элементов (МКЭ) является одним из самых эффективный методов для численного решения задач математической физики и техники [1-7]. Первые попытки применения метода конечных элементов к расчету волноведущих систем относятся к середине шестидесятых годов прошлого века. К этому времени техника этого метода была хорошо разработана, и он успешно применялся для решения граничных задач механики. Были выписаны и исследованы вариационные функционалы для волноводов произвольного типа [8]. Первые работы были посвящены расчетам с помощью МКЭ металло-диэлектрических волноводов, а затем разработанная в этих работах техника была обобщена на открытые волноводные структуры [9-13].

Настоящая работа посвящена разработке и реализации эффективного численного алгоритма расчета основных характеристик волноводов с диэлектрическим заполнением на основе метода конечных элементов. Рассматривается полная векторная постановка задачи. В результате использования техники метода конечных элементов задача сводится к обобщенной алгебраической проблеме собственных значений [12]. Матрицы, фигурирующие в постановке этой задачи, являются незнакоопределенными симметричными ленточными матрицами высокого порядка. Их факторизация, необходимая для последующего использования итерационных методов, требует специальной стратегии, позволяющей учесть все вышеперечисленные особенности. В качестве таковой была разработана стратегия, основанная на одном из вариантов метода Банча - методе Банча-Кауфман, которая предоставляет возможность разрешить проблемы, связанные со спецификой данных матриц [13, 14]. Важное место в данном алгоритме занимает использование оптимальной схемы хранения элементов матриц, а также методики построения таких матриц, одним из вариантов которой является использование опорных матриц.

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

В регулярном волноводе прямоугольного поперечного сечения Э с идеально проводящими стенками [15]

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

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

Е(х, у, г, 0 = Е(х, у, г, )в-ш;

Н(х, у, г, Н(х, у, г, )в-^

и выражая вектор электрического поля Е через вектор магнитного поля Н, запишем уравнение для магнитной составляющей поля:

rot (е^1 rot И) _ k2^H = 0.

(1)

Граничное условие для вектора Е в случае идеальной проводящей стенки [Е х п]до = 0 переходит в два граничных условия для вектора Н:

1

rot И х n

Hx - dxHr

dD

= 0,

dD _ _

= ny (<%HZ _ Hy

(2)

dD

Граничные условия (2) являются для поставленной задачи естественными [3].

Рассматривая модовые решения вида Н(х, у, г) = = Н(х, у)в'вг, где в — постоянная распространения, и производя замену переменных Нх = ; Ну = ;

Нг = Нг, после несложных преобразований получим систему уравнений, записанную в матричном виде:

аН = в2вН. (3)

Здесь матрицы А и В имеют следующий вид:

A =

д_ дх

е I)

Н)

дУ v дУ/ м дУ V дх)

0

_дх Ю _ 0

0

0

Z

0

В

— -

— -

д_

дх

д_ дх

„-1

дх

ду

д_ дх д ду

е я— - - ю - ^

(4)

а Н ={Их,Ну,Нг}Т .

Щ = £ Фа Щ

а=1,т,п

ну = £ ф„ (Ж

а=1,т,п

Щ = £ Фа (Н

А1Н1 = в2В1Н1,

(6)

где

Ну

; ; (Ч ; (^)т ; (^)т ; (^т ;

Т

(Ч; (Ч,; ЧЬ (7)

2. Построение матриц

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

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

Построим элементные матрицы. Аппроксимируем компоненты вектора магнитного поля Н на элементе Т1 следующим образом:

(5)

где — значение Нц в а-м узле, а = I, т, п —

вершины треугольника Т1, ц = х, у, г, Фг = х/Нх; Фт = у/Ьу; Фп = 1 — х/Нх — у/Ну — линейные базисные функции.

Рассмотрим уравнение (3), матрицы А и В которого имеют вид (4), на элементе Т1 и заменим вектор Н

аппроксимирующим его вектором Н] = |Щ, Щ, Н^ .

Применяя стандартную технику метода конечных элементов и учитывая вид матриц (4) и формулы (5), получим уравнение (3) на элементе Т1 в виде

а матрицы А1 и В1 размером 9 х 9 имеют блочную структуру. Блоками являются квадратные матрицы: А'а1, а, 7 = г, т, п, размером 3 х 3. Матрицы А1 и В1, связанные с элементом Т1 , называются элементными [16].

Для элемента Т2 проведем аналогичную процедуру построения элементных матриц АН2 и ВН2 .

Рассмотрим теперь процесс сборки матриц. Обозначим протяженность сечения О вдоль оси х через ¿х = Нх х Ых, а вдоль оси у — через Ьу = Ну х Ыу. Тогда сечение О будет разбито на N элементов, где N = 2Ых х Ыу. Занумеруем элементы снизу вверх, переходя от столбца к столбцу слева направо. Аналогичным образом пронумеруем и узлы нашей сетки (число узлов Ыи = (Ых + 1) х (Ыу + 1)). Каждому узлу

ставится в соответствие три величины Нх, Ну, Нг — компоненты вектора магнитного поля.

Переформулируем задачу (6), (7) таким образом, чтобы полученная задача включала в себя все элементы вида Т1 или Т2.

Произвольному элементу Т1 с номером £ соответствуют узлы с номерами N1, М + 1, Ы\+Ыу+1. Уравнение (6) будет иметь на элементе Т1 вид

А! Нк = в2 В1 Нк,

где

Нк = [ Щ

Ы1+Ыу + 1

; [Ну

; К

М+^ + 1

Нх) ; [НА ; [Нг

М + 1 V /М+1

г, ; (Нх) ; (нЛ ; (Нг) }

N+1 V /М \ /М \ / N1)

Теперь уравнение (6) можно записать следующим образом:

А £ Ни = в2Вк Ни. (8)

Здесь вектор Ни имеет вид

Ни н {Щх) ■; 1; (Н) 1; (Н) 2;...;(Н),

а блочная матрица А£

(в ••• в в

А* =

вв

в в

в в

АН1 АН1

пп п

АН1 АН1

т п гг

АН1 АН1

вв

АН1 Аш

АН1 А

АН1

в

• в • в

•в

•в

где блоки А'а1, а, 7 = г, т, п, размером 3 х 3 расположены на пересечении строк и столбцов с номерами N1, М + 1 и N1+Ny+1, а в — нулевая матрица

0

£

0

£

д

д

а

а

а

а=1.т.п

размером 3 х 3. Матрица Вк имеет вид, аналогичный матрице Ак. Элементные матрицы оказываются, таким образом, «вмороженными» в нулевые опорные матрицы.

Проводя аналогичные рассуждения для всех элементов типа Т1, приходим к Лх х Лу уравнениям типа (8) для элементов с номерами к = 1,3, ...,2ЛХ х Ny — \. Аналогично, проводя для элементов типа Т2 подобную процедуру, снова получаем Лх х Лу уравнений для элементов с номерами к = 2,4,..., 2Лх х Лу.

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

AHU = в2ВНи,

(9)

где

'2NXN„

'2NXN„

A = IE AJ ; H = |£H

k=1

k=1

Сборка матрицы жесткости A и матрицы Н сводится к простому суммированию опорных матриц.

Фигурирующие в уравнении (9) матрицы A и Н

обладают следующими свойствами. Матрицы A и Н являются симметричными, однако не являются знако-

определенными. Матрица A — вырожденная, поэтому прибавим к обеим частям равенства (9) выражение -к2^НННи, где Н = max е. В результате получим задачу на собственные значения, матрица которой уже не является вырожденной. Матрицы A и Н являются ленточными матрицами, причем для приведенного разбиения прямоугольного сечения D принятая нумерация узлов является оптимальной с точки зрения минимизации ширины ленты, которая оказывается равной Sh = 3(Ny + 2).

3. Факторизация матрицы жесткости с использованием стратегии Банча-Кауфман

Исходная задача (3) с помощью техники метода конечных элементов сводится к обобщенной алгебраической проблеме собственных значений (9). Для ее решения применяется метод обратных итераций со сдвигом [17], что приводит к решению итерационной системы уравнений вида

A - \kH) yk+i = Hxk, k = 1,2,...

(10)

При факторизации матрицы ^А — ХкВ^ возникает

ряд проблем. Поскольку матрицы А и В, а следовательно, и матрица Лк = А — ХкВ — незнакоопределенные, то невозможно применение метода Холецкого [12]. Применение методов с пивотированием (выбором главного элемента) по строкам либо по столбцам, также невозможно, поскольку матрица Лк уже перестает быть ленточной. Нельзя использовать и метод Гаусса без пивотирования, поскольку на главной диагонали матрицы в процессе факторизации могут появиться малые элементы.

Наиболее целесообразным в данном случае оказывается применение метода, основанного на предложенной в диссертации Джемса Банча (Калифорнийский университет в Беркли) в 1969 г. стратегии устойчивого

алгоритма для обработки симметричных незнакоопре-деленных матриц и ее модификации, известной как стратегия Банча-Кауфман [13].

Пусть симметричная невырожденная матрица А = = (а^), ¡, Ц = 1, п, представлена в блочном виде:

[б сЛ

А = , где Б — невырожденная подматри-

\С в/

ца размера к х к, к < п. Тогда матрица А может быть представлена в виде произведения трех матриц А = Ь1Б1ЬТ, причем Ь1 — нижнетреугольная матрица,

а В1 = | блочно-диагональная матрица.

V в

Производя аналогичную процедуру с матрицей В1, т.е. представляя ее в виде В1 = Ь2Б2ЬТ2, записываем матрицу А как произведение матриц Ь1Ь2Б2Ь12 Ь\. Здесь В2, В2 — матрицы размера (п — к) х (п — к);

В В (Б о\

Ь2 — нижнетреугольная матрица, В2 = _ —

\0 Д/

блочно-диагональная, Б1 — матрица размера к1 х к1, а матрицы Ь2, В2 размера п х п имеют блочную структуру и содержат в качестве блоков матрицы Б, Б1,

¿2 и ¿2.

Повторяя затем эту процедуру, производим разложение матрицы А и получаем ее представление в факто-ризованном виде: А = Ь1Ь2 • • • ¿¡-^¡¿¡¿ТЬ^ 1 • • • ¿Т¿Т = = ¿¿¡¿Т, где Ь — нижнетреугольная матрица.

В разработанной Банчем и Кауфман стратегии предполагается в качестве Бк на к -м шаге выбирать матрицу 1 х 1 (скалярный шаг исключения) или 2 х 2 (блочный шаг) в зависимости от некоторого условия.

Применяя теперь описанный алгоритм к матрицам Д, В2 и так далее, производим разложение матрицы А до тех пор, пока получаемая матрица В не будет иметь порядок 1 х 1 или 2 х 2.

Матрицы А и задачи (6) являются разреженными матрицами ленточной структуры. Но поскольку при факторизации матриц происходит заполнение внутри ленты, компактное хранение их в строчном формате одним из стандартных методов оказывается неэффективным. Возможно также при обработке матриц увеличение ширины ленты на единицу. Наиболее соответствующей применяемому методу факторизации матриц, основанному на процедуре Банча-Кауфман, является диагональная схема хранения симметричных матриц [16]. В виде прямоугольного массива хранится нижняя половина ленты. Если матрица А размером N х N такова, что полуширина ее ленты равна ЛБ, то с учетом возможного увеличения ЛБ на единицу при факторизации для ее хранения создается прямоугольный массив ЬБ размером N х (ЛБ + 1). При записи матрицы А в массив ЬБ пересчет индексов осуществляется следующим образом: пусть элемент ац е А записывается как элемент Ь&Н е ЬБ. Тогда индексы (¡, Ц) и (¿, Н) связаны соотношениями g = I, Н = /+ЛБ+1 — I, а обратный переход от (¿, Н) к (¡, Ц) производится по формулам I = g, Ц = g + Н — (ЛБ + 1).

Тестирование реализующей предложенный алгоритм программы происходило на волноводе без заполнения

и на волноводе с кусочно-постоянным заполнением. На рис. 1 показаны варианты заполнения сечения волновода, а на рис. 2-4 показаны дисперсионные кривые основной моды для пустого волновода и волновода с различным кусочно-постоянным заполнением. На всех рисунках точному решению соответствует сплошная кривая. На рис. 2 приведена дисперсионная кривая основной моды пустого волновода. Значения, рассчитанные с помощью предлагаемого алгоритма с использованием процедуры Банча-Кауфман (БК-ал-горитма), отмечены кружками. На рис. 3 приведена

а

б

Рис. 1. а — Симметричное трехкомпонентное сечения волновода: диэлектрическая проницаемость центральной части равна £1, диэлектрическая проницаемость боковых частей равна £. б — Двухкомпонентное сечения волновода: диэлектрические проницаемости частей равны £1 и £

0.80

0.60

0.40

0.20

3.0 3.5 4.0 4.5 5.0 5.5 к!2л

Рис. 2. Дисперсионная кривая основной моды пустого волновода. Сплошная кривая соответствует точному решению, рассчитанные значения отмечены кружками

1.20

1.00

0.80

0.9 к!2к

Рис. 3. Дисперсионная кривая основной моды волновода с симметричным трехкомпонентным сечением при следующих значениях геометрических и материальных параметров: d/a = 0.5; £1 = 2; £ = 3. Сплошная кривая соответствует точному решению, рассчитанные значения отмечены кружками

0.7 к/ 271

Рис. 4. Дисперсионная кривая основной моды волновода с двухкомпонентным сечением при следующих значениях геометрических и материальных параметров: d/a = 0.5; £1 = 2; £ = 1. Сплошная кривая соответствует точному решению, рассчитанные значения отмечены кружками

дисперсионная кривая основной моды волновода, сечение которого имеет вид, изображенный на рис. 1, а). Кружками отмечены значения, рассчитанные на основе БК-алгоритма. На рис. 4 приведена дисперсионная кривая основной моды волновода, сечение которого имеет вид, изображенный на рис. 1, б). Кружками отмечены значения, полученные с использованием БК-алгоритма. Результаты тестирования показали, что применение для факторизации возникающих в методе конечных элементов матриц процедуры, основанной на стратегии Бан-ча-Кауфман, позволяет построить весьма эффективный алгоритм расчета дисперсионных кривых металло-ди-электрических волноводов со сложным заполнением. При этом важное место в данном алгоритме занимает использование оптимальной схемы хранения элементов матриц, а также методики сборки таких матриц, основанной на применении опорных матриц.

Заключение

Использование методики, основанной на применении при факторизации матриц метода Банча — Кауфман в сочетании с диагональной формой хранения элементов матриц и применения опорных матриц при матричной сборке, позволяет построить достаточно простой и надежный алгоритм численного решения спектральных задач математической теории волноводов в полной векторной постановке. Он может быть распространен на случай анизотропного, кирального и биизотропного заполнения, а также использован для построения блока расчета прямой задачи в алгоритме решения задачи синтеза волноводов. Задачи синтеза ставятся как задачи математического проектирования для определения основных параметров волновода, при которых он обладает требуемыми техническими и эксплуатационными характеристиками. Характерной особенностью таких задач является то, что, в отличие от обратных задач распознавания, в задачах синтеза, как правило, отсутствует требование единственности решения. Наиболее полный и универсальный подход к решению задач синтеза волноводов заключается в рассмотрении таких задач как математически некорректных, с применением для их решения метода регуляризации А.Н.Тихонова. Такой подход был предложен в работах А. Г. Свешникова и А. С. Ильинского [18-19].

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

Работа выполнена при поддержке РФФИ (проект № 12-01-00479).

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

1. Зенкевич О. Метод конечных элементов в технике. М.., 1975. Zienkiewicz O.C. The finite method in engineering science. McGraw-Hill, 1973.

2. Норри Д., де Фриз Ж.. Введение в метод конечных элементов. М., 1981. Norrie D.H., de Vries G. An introduction to finite element analysis. Academic Press, 1978.

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

4. Боголюбов А.Н., Боголюбов Н.А., Свешников А.Г. // Физические основы приборостроения. 2013. 2, № 1. С. 10.

5. Боголюбов А.Н., Делицын А.Л., Красильникова А.В. и др. // Зарубежная радиоэлектроника. Успехи современной радиоэлектроники. 1998. № 5. С. 39.

6. Боголюбов А.Н., Буткарев И.А., Минаев Д.В. // Радиотехника и электроника. 2005. 50, № 2. С. 140.

7. Боголюбов А.Н., Делицын А.Л., Лавренова А.В. // Электромагнитные волны. 2004. 9, № 8. С. 22.

8. Никольский В.В. Вариационные методы для внутренних задач электродинамики. М., 1967.

9. Bossavit A. // IEEE Proc. 1988. 135, N 8. P. 493.

10. Svedin J. // IEEE Trans. Microwave Theory Tech. 1989. 37. P. 1708.

11. Bermudes A, Pedreira D.G. // Numer. Math. 1992. 61, N 2. P. 39.

12. Уилкинсон Дж. Х. Алгебраическая проблема собственных значении. М., 1970. Wilkinson J.H. The Algebraic Eigenvalue Problem. Clarendon Press, 1965.

13. Bunch J.R., Kaufman L. // Mathematics of Computation. 1977. 31 . P. 163.

14. Bunch J.R., Rose D.J. Sparse matrix computations. Academic Press, 1976.

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

16. Писсанецки С. Технология разреженных матриц. М., 1988. Pissanetzky S. Sparse Matrix Technology. Academic Press, 1984.

17. Калиткин Н.Н., Альшина Е.А. Численные методы. Кн. 1.Численный анализ. М., 2013.

18. Боголюбов А.Н., Красильникова А.В., Минаев Д.В., Свешников А.Г. // Радиотехника. 1997. № 1. С. 81.

19. Боголюбов А.Н., Красильникова А.В., Минаев Д.В., Свешников А.Г. // Математическое моделирование. 2000. 12, № 1. С. 13.

Calculation of waveguides by the finite element method using the Bunch-Kaufman procedure Yu.V. Mukhartova, N.A. Bogolyubova

Department of Mathematics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia. E-mail: a russell67@yandex.ru.

Abstract—An algorithm for the numerical solution of the problem for eigenvalues in a waveguide in the complete vector statement using the finite-element method and the Bunch—Kaufman procedure for the factorization of the matrix of the derived set of linear algebraic equations was built and implemented.

Keywords: metal-dielectric waveguide, finite element method, Bunch-Kaufman method, factorization of a matrix. PACS: 02.60.Cb. Received 2 February 2014.

English version: Moscow University Physics Bulletin 3(2014). Сведения об авторах

1. Мухартова Юлия Вячеславовна — канд. физ.-мат. наук, ст. науч. сотрудник; тел.: (495) 939-10-33; e-mail: muhartova@yandex.ru.

2. Боголюбов Николай Александрович — аспирант; тел.: (495) 939-10-33, e-mail: russell67@yandex.ru.