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

Вычислительные технологии

Том 10, № 2, 2005

РЕШЕНИЕ НЕСТАЦИОНАРНЫХ УРАВНЕНИЙ МАКСВЕЛЛА ДЛЯ СРЕД С НЕОДНОРОДНЫМИ СВОЙСТВАМИ МЕТОДОМ КОНЕЧНЫХ

ОБЪЕМОВ*

А. С. Лебедев, М. П. Федорук, О. В. ШтыринА Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: sasa@ict.nsc.ru, mife@ict.nsc.ru

A numerical method to solve time-dependent Maxwell's equations is suggested. The method uses the finite volume approach on unstructured triangular grids. Numerical examples are given demonstrating the second order of accuracy for both homogeneous and inhomogeneous media.

Введение

Уравнения Максвелла являются фундаментальными уравнениями математической физики и имеют широкий спектр практических приложений. Конечно-разностные методы для численного решения нестационарных уравнений Максвелла начали интенсивно развиваться после выхода в свет пионерной работы [1], в которой был предложен метод второго порядка точности по временной и пространственным переменным, основанный на введении смещенных сеток.

К настоящему времени эти методы получили значительное развитие. Обзор конечно-разностных алгоритмов и примеры их применения для многочисленных практических приложений можно найти, например, в [2-5]. Однако данные алгоритмы применимы для численного моделирования в областях относительно простой формы, в которых можно ввести ортогональную систему координат: декартову, цилиндрическую или сферическую. Вместе с тем потребности фундаментальной науки и прикладных исследований привели к необходимости использовать уравнения Максвелла для моделирования электромагнитных полей в областях со сложной геометрией границ. В этих случаях границы расчетной области не совпадают с координатными линиями, что для конечно-разностных методов приводит к потере точности при постановке граничных условий и, как следствие, к потере точности получаемых численных решений. Использование преобразования координат и решение соответствующим образом преобразованных уравнений уже не в исходной физической области, а в расчетной логической области простой формы является лишь частичным выходом из положения. Кроме того, моделируемые объекты могут обладать

*Работа выполнена при финансовой поддержке Президента РФ (грант № НШ-2314.2003.1) и Российского фонда фундаментальных исследований (гранты № NWO-РФФИ 047.016.003, № 02-01-001029).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

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

Наиболее универсальные алгоритмы решения подобных задач получаются введением в моделируемой области неструктурированной сетки и применением метода конечных объемов или метода конечных элементов для численной дискретизации исходных уравнений Максвелла. Отметим, что к настоящему времени в литературе описаны алгоритмы для решения данной системы уравнений, основанные как на методе конечных элементов [6-8], так и методе конечных объемов [9, 10].

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

1. Численный алгоритм

Система уравнений Максвелла в отсутствие зарядов и токов имеет следующий вид:

д Б

Ж

- с гоШ = 0, Б = еЕ;

дВ

——+ с го1Е

0, В

&уБ = 0; &уВ = 0.

^Н;

(1)

(2)

(3)

(4)

Здесь Е и Н — напряженности электрического и магнитного полей соответственно; Б — вектор электрической индукции; В — вектор магнитной индукции; е — диэлектрическая проницаемость среды; ^ — магнитная проницаемость среды; с — скорость света.

Уравнения (1) и (2) покомпонентно в декартовой системе координат записываются следующим образом:

'дЯэ

дР1

~дГ дР2 д£ дРэ д£ дВ1

~дГ дВ2 д£ дВз

дЯ2

с

с

с

+ с

+ с

+ с

дХ2 дхз

д#1 дЯз

дхз дх1

дЯ2 дН1

дх1 дХ2

дЕз дЕ2

дХ2 дхз

'дЕ1 дЕз

дхз дх1

дЕ1

0,

(5)

дх дх2

0.

Далее везде полагаем ^ = 1, т. е. В = Н.

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

0

0

0

0

плоскости пространственных переменных (х\, х2). Запишем систему уравнений Максвелла (5) в консервативной форме

где и — вектор консервативных переменных; Е1 = Al.iV и Е2 = A2V — векторы потоков. В уравнениях (6) и (7) векторы И и V связаны матрицей перехода в = в(е): V = ©И. Конкретный вид матриц А1, А2, в и векторов И, V приводится далее при описании тестовых расчетов.

Будем считать, что в расчетной области построена сетка с треугольными ячейками Дг таким образом, что если две различные ячейки пересекаются, то они имеют общую вершину или общее ребро (рис. 1).

Интегрируя уравнение (7) по треугольному элементу Д, получим уравнение

где Г — граница ячейки; п = (п1, п2) — внешняя нормаль.

Определим вектор Ип в барицентре Хв треугольного элемента, т. е. в точке пересечения медиан в момент времени Ьп = пт, где т — шаг по времени. Аппроксимируя в уравнении (8) производную по времени конечной разностью, для ¿-й ячейки получаем уравнение

(6)

или

(7)

(8)

где — площадь ячейки.

Рис. 1. Две соседние ячейки треугольной сетки Д^ и Ад с барицентрами Х^ и ХД соответственно и общей стороной, середина которой находится в точке Xе, в плоской геометрии (X = [хг,{ = 1, 2}).

Чтобы из (9) найти и™+1, надо вычислить интегралы по ребрам Гд треугольника Д^. Эти интегралы представляют собой потоки искомых величин через ребра. Предположим, что в треугольном элементе искомые функции изменяются линейно. Тогда по значениям этих функций в барицентре треугольника и по известным (предварительно найденным) градиентам этих функций можно найти значения функций в центре ребра со стороны ¿-го треугольника. Аналогичным образом находим значения функций в центре ребра со стороны треугольника, находящегося по другую сторону ребра. Далее, используя значения функций по разные стороны ребра, находим потоки в центре ребра. Тогда интегралы в (9) приближенно вычисляем по формуле средних прямоугольников.

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

1.1. Нахождение градиентов функций в ячейке

Рассмотрим задачу нахождения градиентов компонент электромагнитных полей. Пусть / (Хв) = /в — значение функции в барицентре ячейки сетки, тогда / (Х^) = / (ж^7, ж^7^ = /В — значение / в барицентре ]-го треугольника.

Уравнение плоскости, которая проходит через точку (Хв, /В и имеет наклон, задаваемый значениями функции /В, ] = 1,2,3 в трех не лежащих на одной прямой точках

ХВ, ] = 1, 2, 3, записывается следующим образом:

/ (X) = / (Х1, Х2) = /оВ + а (Ж1 - жВ0) + в (Ж2 - жВ0) , (10)

где

а = -Дз1 (жВ2 - жВ1) + Д21 (жВ3 - жВО : в = Дз1 (жВ2 - жВ1) - Д21 (жВ3 - жВ1) ,

Д = /В /В Д = /зВ /В Д = / (жВ2 - жВ1) (жВ?2 - жВ?1)

Д21 = ' Дз1 = ' Д = ае^ (жВз - жВ1) (жВз - жВ1)

Предварительное вычисление коэффициентов а, в для треугольной ячейки с координатами барицентра ХВ = {жВо, г = 1, 2} зависит от того, сколько соседей имеет эта ячейка:

1) ячейка имеет трех соседей (не является граничной), тогда а,в находятся по значениям функции в барицентрах этих соседей;

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

Д <8,

р (ХВ - хв) р (ХВ - ХВ)

где р(Х - У) = у/(ж1 - У1)2 + (ж2 - У2)2;

3) имеет одного соседа, тогда в качестве а, в берутся значения а, в для этого соседа. Для каждого ]-го треугольника по найденным а, в вычисляются значения / (ХР) = /Р в узлах сетки ХР = {жР}. В одном узле ХР получаются столько значений / (ХР), сколько треугольников имеют данный узел в качестве одной из своих вершин. Затем величина / (ХР) находится как среднее арифметическое всех этих значений.

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

(хрк, /Р^, к = 1, 2, 3, где Хрк = |Хр ^ — вершины ^-го треугольника, к — номер верши-

1.2. Нахождение компонент электромагнитных полей в центре ребра

Введем в рассмотрение матрицы А1 = вЛ1 и Л2 = ©Л2. Тогда система уравнений (6) переписывается в виде

д д д - V + Л^-— V + Л2—V = 0. дъ дх1 дх2

Обозначим через Хв барицентр треугольного элемента, а через Xе — центр его ребра. Тогда величину V(Xе) можно вычислить по следующей формуле со вторым порядком точности по времени и пространству:

V(Xе) = V(XB) + 5 ■ (Xе - Хв) - Т ( Л^V(XB) + V(XB)

дд

2 V 1 дх! 2 дх2

матрица градиентов компонент электромагнитных полей

х=хв '

где 5 = {в-} = { ъх,

в барицентре ячейки.

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

VL(Xе) = V(XB) + Бь ■ (Xе - хв) - 2 (л!дх;V(XB) + Л2V(XB)

Т / д д

) = V(XR) + Бл ■ (Xе - хB) - - (А1 — V(XR) + Л 2—V(XB)

1.3. Вычисление потоков через границу ячейки

Запишем уравнение (6) в ортогональной системе координат (п, в)

д д д

дх1 = П1 — дп дв

д д д

дх2 = П2Т-дп + П1 —, дв

где п = (п;,п2) — вектор единичной нормали к общему ребру треугольников Д^ и Дд,

д

направленный от треугольника Д^ к треугольнику Дд; —— — производная по нормали,

дп

д

—--производная по касательной.

дв

Тогда получим уравнение

ди + (Ащ! + А2П2) -пV + (-а4Ш2 + А42П^ дV = 0.

Рассмотрим одномерную задачу Римана для уравнения

д д

т и+Адпу = (11)

где А = А1п1 + А2п2, и через ее решение вычислим поток сквозь границу ячейки. Пусть 4 = \Je-l4д. Уравнение (11) запишем в виде

дд

а;и + ВдПи = 0 В = Ав<4).

Пусть Я — матрица, столбцы которой являются собственными векторами матрицы В. Тогда О = Я-1ВЯ, О = О- + О+, где О — диагональная матрица, составленная из собственных значений матрицы В, — диагональные матрицы, полученные из О заменой всех отрицательных (положительных) собственных чисел нулями. Имеем

В = ЯОЯ-1 = ДР+Я-1, + ДР-Я-1у = В+ + в-.

в+ в-

Для нахождения потоков будем решать систему уравнений

д ( ) д в-1 (4) ^У + (В+ + В-) в-1 (4 ) дПV = 0

или, вводя обозначения С = В±в 1 (4),

в-1 <4> 1у + (С + + С-) дПУ = а

Таким образом, поток вектора и через общее ребро треугольников Д^ и Дд в направлении нормали п (см. рис. 1) равен Е* = С+У* (Xе) + С-УД (Xе). Из уравнения (9) получаем

з

И+1 = И - ьЕ*,

з

ИД+1 = ид + д ЕД

или

з

П+1 — ЛТп__1к

УП+1 = V - ^в(еь^/АЬЕ*

Д=1

з

УД+1 = уд + — в(ед^/А л ЕД

где 5АЬ и 5Ал — площади треугольников Д* и Дд соответственно; т — шаг по времени; /АЬ и /Ад — длины ребер с номером к (локальная нумерация по к независимая); Е* и ЕД — потоки через ребра с номером к.

2. Результаты тестовых расчетов

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

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

( ЗЕл дНз

дЬ дх

2

дЕ2 дН3 < — + 3

0,

0,

дЬ дх

дНз дЕх + дЕ2

0.

дЬ дх2 дхх Система (12) принимает консервативный вид для вектора

(12)

и = V

Е1 Е2 Н3

Матрица перехода

100 ей о 1 о |.

0 0 1

Перепишем систему (12) в матричном виде:

д д д - V + Ах-— V + А2т— V = 0, дь дхх дх2

где

000 Ах ^(001 010

А2

0 0 -1 000 1 0 0

Одномерную задачу Римана для нахождения потоков будем решать для уравнения

дЬ V + А дп V 0'

где

/ 0 0 -п2

А= 0 0 п1

\ -п2 П1 0

Определим матрицы

с +=2

п22 -п1п2 -п2 с-=2 ( -п22 п1п2 -п2

-п1п2 п21 п1 п1 п2 - п 21 п1

-п2 п1 1 -п2 п1 -1

Тогда поток из L-треугольника в R-треугольник равен

Fl = C+VL (Xc) + C-VR (Xc

Система уравнений (12) имеет аналитическое решение:

а

E1 = — — sin (ах2) cos (вж1 — —í); Р

E2 = cos (ах2) sin (вх1 — —í);

H3 = — cos (ах2) sin (вх1 — —í) ,

Р

(13)

(14)

(15)

где а2 + в2 = , а = в = п.

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

0.75'

0.5

0.25

0.75-

0.5-

0.25-

00 025 0:1 0.75

00 0.25 0.5 0.75

Рис. 2. Расчетная область покрыта регулярной (а) и нерегулярной (б) треугольными сетками с 10 ячейками вдоль одного ребра единичного квадрата.

а

x

x

2

2

1

1

Рис. 3. Сетка в многосвязной области с 5 ребрами вдоль каждой стороны единичного квадрата (934 ячейки во всей области).

Рис. 4. Сетка в многосвязной области с 10 ребрами вдоль каждой стороны единичного квадрата (3784 ячейки во всей области).

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

В общем случае для областей сложной формы применялся разработанный авторами алгоритм триангуляции, итерационным образом строящий сетку внутри области по заданным на ее границе узлам. Примеры сеток с плавно изменяющимся размером ячеек для квадрата с вырезом показаны на рис. 3, 4.

Расчеты в области с вырезом проводились на последовательности из четырех сеток, содержащих 934, 3784, 14 915 и 57072 ячеек. На рис. 5 показаны изолинии первой компоненты электрического поля Е\ после трех периодов, полученные по схеме первого порядка на трех различных сетках, содержащих 934, 3784 и 14 915 ячеек (а — в соответственно). Аналогичные результаты, полученные по схеме второго порядка, показаны на рис. 6. Из этих рисунков видно качественное превосходство схемы второго порядка.

Как следует из рис. 5, 6, данный алгоритм позволяет выполнять расчеты в областях произвольной формы, в том числе многосвязных.

Для оценки отклонения численного решения от аналитического введем величину ин-

0.75-

0.5-

0.25-

0.75

0.5

0.25

Ь

"о \ ю с^ /

-Ахч У\) о оэ ¿о / / со / о' /

\ \ г

^ / ¡р I (V о / Д1Л-

О" <р

-0.125—

"025 0.5 075"

0.75

0.5

0.25

___-0.12;

"025 05 аП

Рис. 5. Изолинии первой компоненты напряженности электрического поля Е\ после трех периодов распространения, полученные по схеме первого порядка на трех различных сетках: а — 934 ячейки, б — 3784 ячейки, в — 14 915 ячеек.

0.75-

0.5-

0.25-

0.75

0.5

0.25

0.75

0.5

0.25

а

в

X

X

X

2

2

2

0

0

X

1

1

а

в

X

X

X

2

2

2

1

1

Рис. 6. Изолинии первой компоненты напряженности электрического поля Е\ после трех периодов распространения, полученные по схеме второго порядка на трех различных сетках: а — 934 ячейки, б — 3784 ячейки, в — 14 915 ячеек.

тегральной погрешности 8 в пространстве Ь2 в момент времени ¿га = пт:

"В\ лгп

8

V™ — V™ IX1

' у гамт, V ) V ^11*2

|УП„„ (Хв )ИЬ2

\

т

г=1

Е К (ХВ) - Е? (ХВ))2 + к (ХВ) - Е? (ХВ))2 + (и? (ХВ) - я? (хв))

А,-

т

г=1

Е (Е? (ХВ))2 + (Е? (ХВ))2 + (я? (ХВ)) 02

А,

где Т — общее количество треугольных ячеек в вычислительной области; (Хв)

и У?гаа (Хв) — вычисленные и точные (см. уравнения (13)-(15)) значения электромагнитных полей в барицентре ячейки соответственно. В табл. 1 показаны максимальные Ь2-погрешности 8 для различных вычислительных схем и расчетных сеток.

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

Тест 2. В качестве второго теста рассмотрим распространение ТМ-волны в среде с разрывным значением диэлектрической проницаемости 4:

4 = 4 (х1)

4 при |х1| < 0.5,

1 при |х1| > 0.5. Область, в которой проводились расчеты, показана на рис. 9.

2

Т аб л и ц а 1

Максимальные вычислительные ¿2-погрешности 5 в расчете ТЕ-волны в однородной среде в декартовых координатах

Тип сетки Тип области Число 5 (нулевые 5 (линейная

треугольников градиенты) аппроксимация

по внешней градиентов)

границе

Нерегулярная Односвязная 10 х 4 0.1339 0.008995

20 х 4 0.07153 0.002270

40 х 4 0.03694 0.000571

80 х 4 0.018855 0.000144

160 х 4 0.00953 0.00002038

Многосвязная 5 х 4 0.179808 0.0217624

10 х 4 0.095329 0.0052643

20 х 4 0.0498402 0.0014573

40 х 4 0.025615 0.0003555

Регулярная Односвязная 16 х 4 0.072144 0.002047

32 х 4 0.037155 0.00051

64 х 4 0.018876 0.0001274

128 х 4 0.00952 0.00003184

256 х 4 0.00478 0.00000796

0.02

N=20, т=0.007

10"3

; Л/\/\/\, I/ N=40,1=0.0035 / N=40,1=0.0035 7

: / :/ ? Д,,. ! : , I ' N=80,1=0.0018 "( "' , '"' Г '", '"'' Г '", '"' 1 * 10"4 - N=80,1=0.0018 "............

N=10,1=0.014

N=20,1=0.007

8 10

6 8

10

Рис. 7. Эволюция ¿2-погрешности в расчете ТЕ-волны на последовательности нерегулярных треугольных сеток: а — нулевые градиенты, квадрат; б — линейная аппроксимация градиентов, квадрат.

у'УУУУУУ

N=20,1=0.0005 N=40,1=0.00025

8 102

10

ю

10'

N=5, 1=0.002

N=10,1=0-001 N=20,1=0.0005

N=40,1=0.00025

t

10

а

а

Рис. 8. Эволюция ¿2-погрешности в расчете ТЕ-волны на последовательности нерегулярных треугольных сеток: а — нулевые градиенты, квадрат с вырезом; б — линейная аппроксимация градиентов, квадрат с вырезом.

8=1

8 = 4

8 = 4

8 = 1

-1 -1/2 0 1/2 1 Х1

Рис. 9. Расчетная область распространения ТМ-волны в среде с плоским диэлектрическим слоем.

В случае ТМ-волны Е = = Н3 = 0, исходная система уравнений (6) в безразмерных переменных имеет следующий вид:

0, 0.

Векторы и, V и матрица перехода в (е) равны:

3 5Я1

+

5£ 5^2

5Я1 5Ез

+

Ы 5^2

5Я_ 5Ез

1 5£ 5ж1

0,

(16)

и

Н1

V

Ез

Н1

в(е)

е 0 0 0 1 0 0 0 1

Запишем систему (16) в матричном виде относительно вектора V:

0,

5 5 5

- V + Л— V + V 5г 5ж1 5Ж2

где

А1

( 0 0 - ^

е

0 0 0 -1 0 0

/0 1 0^

е

100 000

Одномерную задачу Римана для нахождения потоков будем решать для уравнения

5£ И + 5п V 0

с матрицей

А

0

П2 0

-П1 0

—п1 0 0

Определим матрицы

( Vе

1 2

С + = 1

V

П1

п2

—п1 \

П1П2

'е

П1П2

п

2

е

С-

( — ^ П2

V7! V7! /

V

— П1

е

П1П2

—п1 \

П1П2

Поток из ¿-треугольника в Д-треугольник равен = С+VL (Xе) + С VR (Xе Система уравнений (16) имеет следующее аналитическое решение:

Ез

2п

ехр

3

'-у/Эп

1

2 сое — хН (Кух2) сое , где |х1| < -, х2 > 0,

1

3 (1 — 2|х1 |) 1 вт (КуХ_) сое , где | > ^, X > 0,

1

Таблица2

Максимальные вычислительные ¿2-погрешности 5 в расчете ТМ-волны в среде с плоским диэлектрическим слоем в декартовых координатах

Тип сетки Число треугольников по внешней границе 5 (нулевые градиенты) 5 (линейная аппроксимация градиентов)

Нерегулярная 16 х 6 0.21503 0.0145

32 х 6 0.1171 0.00374

64 х 6 0.06084 0.000943

128 х 6 0.0313 0.000256

Регулярная 16 х 6 0.185 0.00791

32 х 6 0.0989 0.00200

64 х 6 0.0512439 0.000502

128 х 6 0.0262 0.000126

Н1

2Ку /2п \ . . . , , 1

-сое I — жН сое (Куж2) 81п (и£) , где |х1| < -, ж2 > 0,

и V 3 / 2

Ку

и

ехр

5п

3

(1 — 2|ж1|) 1 сов (Куж2) 81п (и£) , где |ж1| > -, ж2 > 0

2

4п /2п \ . ^ . . . . . 1 — — 81п —- Ж1 й1п(Куж_)81п(и£) , где |Ж11 < -, Ж2 > 0, 3и V 3 / 2

2у/3п ^|ж1| 3и ^ж1

п

1

ехр | о (1 — 2|ж1|) ) 81п (Куж2) 81п (и£) , где |ж11 > ^, ж2 > 0.

2п 13

Здесь Ку = —У у , и

4п

3^ '

0.2 0.18 0.16 0.14 0.12 0.1 0.08

0.06 0.04

0.02

- ;' \ . N=10,1=0.009 ^

¥ А \ , N=20,1=0.0044 V /\ 1\ 1\

/IV / 1 " А \ \1 \1 \/ \ N=40,1=0.0022

! \ \ \ ■ [/ 1 \/ V ; : Г 1 N=80,1=0.001 1 , 1 , 1

5 10"2

10"3

10

10-

N=20,1=0.0044

N=40,1=0.0022 N=80,т=0.001

а

Рис. 10. Эволюция ¿2-погрешности в расчете ТМ-волны на последовательности нерегулярных треугольных сеток: а — нулевые градиенты; б — линейная аппроксимация градиентов.

В табл. 2 показаны максимальные Ь2-погрешности 8, полученные для алгоритмов первого и второго порядков и различных расчетных сеток.

Эволюция ошибки 8 по норме Ь2 для различных сеток показана на рис. 10. Видно, что, как и в предыдущем случае, линейная аппроксимация функций из барицентров треугольных элементов дает алгоритм второго порядка точности.

Заключение

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

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

[1] Yee K.S. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media // IEEE Trans. Antennas and Propagat. 1966. Vol. 17. P. 585-589.

[2] Taflove A. Computational Electrodynamics. The Finite-Difference Time-Domain Method. Boston: Artech House, 1995.

[3] Advances in Computational Electrodynamics. The Finite-Difference Time-Domain Method / A. Taflove (Ed.). Boston: Artech House, 1998.

[4] Computational Electrodynamics. The Finite-Difference Time-Domain Method / A. Taflove, S.C. Hagness (Eds.). Boston: Artech House, 2000.

[5] Sullivan D.M. Electromagnetic Simulation Using the FDTD Method. N. Y.: IEEE Press, 2000.

[6] Assous F., Degond P., Segre J. Numerical approximation of the Maxwell equations in inhomogeneous media by p1 conforming finite element method //J. Comput. Phys. 1996. Vol. 128. P. 363-380.

[7] Hiptmair R. Finite elements in computational electromagnetism // Acta Numerica. 2002. P. 237-330.

[8] Ziming Chen, Qiang Du, Jun Zou. Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients // SIAM J. Numer. Anal. 2000. Vol. 37. P. 1542-1570.

[9] Hermeline F. Two coupled particle-finite volume methods using Delaunay — Voronoi meshes for the approximation of Vlasov — Poisson and Vlasov — Maxwell equations //J. Comput. Phys. 1993. Vol. 106. P. 1-18.

[10] Cioni J.-P., Fezoui L., Issautier D. High order upwind schemes for solving time-domain Maxwell equations //La Recherche Aerospatiale. 1994. N 5. P. 319-328.

Поступила в редакцию 26 августа 2004 г., в переработанном виде —16 ноября 2004 г.