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

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

Том 9, № 5, 2004

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

О. В. Нечаев, Э.П. Шурина Новосибирский государственный технический университет, Россия

e-mail: [email protected]

М. П. ФЕДОРУК

Институт вычислительных технологий СО РАН, Новосибирск Россия

e-mail: [email protected]

An algorithm based on the vector finite element method for numerical solution of quasi-steady Maxwell's equations is considered. This algorithm allows taking into account the continuous condition of the tangential components of electric field on the boundary of media with different properties. The results of test calculations are presented.

Введение

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

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

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

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

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

Максвелла. Отметим, что к настоящему времени в литературе описан ряд алгоритмов, посвященных решению данной системы уравнений на основе как метода конечных элементов [6-11], так и метода конечных объемов [12, 13].

В работе развит вычислительный алгоритм решения квазистационарных уравнений Максвелла на основе векторного метода конечных элементов, использующего так называемые пространства Неделека (Nedelec) [8, 9]. Применение данного алгоритма для расчета электрических полей в материалах с разрывными свойствами позволяет естественным образом учесть условие непрерывности тангенциальных компонент электрического поля на границе раздела двух сред.

1. Математическая модель

Поведение гармонического электромагнитного поля описывается системой уравнений Максвелла

rot Ё + iwB = 0; (1)

rot H - iwD — J = Jo; (2)

div ff = p; (3)

div B = 0; (4)

div J + iwp = 0, (5)

где ё — напряженность электрического поля; ff — электрическая индукция; H — напряженность магнитного поля; B — магнитная индукция; J — плотность тока; J0 — плотность тока источника; p — плотность электрических зарядов; ш — циклическая частота; е — диэлектрическая проницаемость; ц — магнитная проницаемость; a — электрическая проводимость; i — мнимая единица. В линейном случае

ff = еЁ; (6)

B = цЙ; (7)

J = aJ. (8)

Пусть выполняется условие разрешимости системы уравнений Максвелла

div J0 = 0. (9)

Учитывая (3), (6) и (8), закон сохранения (5) можно преобразовать к следующему виду:

div ((a + iwe)J) = 0. (10)

Используя (6)-(8), систему (1), (2) можно записать в виде уравнения второго порядка относительно поля EJ

rot —rotЁ + k2 Ё = -iwJo, (11)

где k2 = iwa — ш2е.

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

[n х J]|r = 0; (12)

[П • (а + iwe)E]|r = 0. (13)

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

n х EU = 0. (14)

2. Векторная вариационная постановка

Для расчетной области П введем пространства Неделека [8, 9] H(rot;Q) = {V G [¿2(П)]3 : rot V G [¿2(П)]3}, Ho (rot; П) = {v G H (rot; П) : v х v|r = 0}

с нормой

I n12ot Q = n • ПМП + / rot n • rot ПМП.

n q

Определим скалярное произведение

( n, V) = I n • vMa

Для задачи (11), (14) сформулируем следующую вариационную постановку [10]: Найти EE G H0(rot; П) такое, что для W G H0(rot; П)

rot —rot E, v ) + (k2E, v) = -i(wJ0, v). (15)

V J

В соответствии с векторной формулой Грина

У [p (rot n • rot v) — n • rot (p rot v)]dn = Jp [( v х rot V) • n]dS.

Q dQ

Уравнение (15) можно представить в следующем виде:

I rot1 rot E • vdn = I 1 rot E • rot vdn — / 1 [(V х rot E) • v]dS. (16)

J V J V J V

Q Q dQ

Из свойств введенных пространств следует, что второе слагаемое в правой части (16) равно нулю.

В результате получаем векторную вариационную постановку Найти EE G H0(rot; П) такое, что для W G H0(rot; П)

1 rotE,rot V ) + (k2E, v) = —¿(wj0, v). (17)

V

Q

Пусть имеет место следующее свойство вложения

gradф € Н(ш^П), Уф € Ь2(П). (18)

Вариационная постановка (17) выполняется для всех V € Н(гс^; П), согласно (18) возьмем V = gradф, ф € Ь2(П), и тогда (17) имеет вид

^-гс^ Ё,rotgradф^ + (к2Ё,gradф) = -г(^50,gradф) Уф € Ь2(П).

Учитывая (9) и свойство rotgrad ф = 0, получим

((^2е + г^)Ё, grad ф) = 0 Уф € Ь2(П). (19)

Уравнение (19) является вариационным аналогом закона сохранения (10), т.е. решение вариационной задачи (17) удовлетворяет закону сохранения заряда (5) в слабом смысле.

3. Дискретная вариационная постановка

Построим в расчетной области гексаэдральную сетку, на ячейках которой определим базисные edge-функции, ассоциированные с ребрами сетки конечномерного подпространства Hh(rot; П) С H(rot; П).

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

N1 = 1 - y)(1 - z)l, N

N3 = 1 - y)z1, N

N5 = 1 - z)(l - N

N7 = 1 - z)xj, N

N9 = 1 - x)(1 - y)k, N

N11 = (1 - x)yk N

2 = y(1 - z)1,

e •

4 = yz:b

6 =z(1- x)j

zxj,

x(l - y)k

e _ 8 = e 10 e 12

xyk.

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

Также построим конечномерное подпространство Ь2(П)^ пространства Ь2(П), в качестве базисных функций возьмем трилинейные узловые функции

Ф1 = (1 - x)(1 - y)(1 - z),

ф3 = (1 - x)y(1 - z), ф5 = (1 - x)(1 - yK

ф7 = (1 - x)yz

x(1 - y)(1 - z), xy(1 - z) ■6- x(1 - y^

ф8 = xyz-

2

Ф4 ФФ

e

Покажем выполнение свойства вложения (18) для конечномерных подпространств H(rot; n)h и L2(n)h, т.е. покажем, что градиент от базисной функции конечномерного подпространства L2(H)h (трилинейной узловой функции) есть линейная комбинация базисных функций подпространства H(rot; H)h

gradФ1 = -((1 - y)(1 - z), (1 - x)(1 - z), (1 - x)(1 - y))

T

-N1 - N5 - N9.

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

grad = N1 - 147 - 1410, grad = — NN2 + NN5 - NN1е, grad ф5 = — NN 3 — N 6 + NN 9,

grad фф grad ф[

N 2 + N 7 — N е N 3 — N 8 + N е

grad ф7 = —N4 + N6 + grad ф8 = N4 + N8 + N

те

е2.

Сформулируем дискретную вариационную постановку: найти Ен € Н^(го1 ; П) такое,

что для V€ ; П)

-го1 Его1 "М + (к2ЕV*) = — г(^Ло, V*).

(20)

В результате получаем следующую систему линейных алгебраических уравнений (СЛАУ):

(21)

А + В — С Ег " /г "

С А + В Ег . /г .

где Ег и Ег — веса в разложении по базису действительной и мнимой компонент поля Е соответственно, а элементы матриц В и С определяются соотношениями

[А]

г,3

-ГО1 NN ГО1 NN ь

ц г' 3

(22)

[В}г,3 = — («^Ч* (23)

[С]г,з = (<7сЛЧ* N*).

Система линейных алгебраических уравнений (21) является несимметричной и плохо обусловленной. Применение стандартных решателей и предобусловливателей, основанных на неполном ЬИ-разложении матрицы, не дает желаемого эффекта. Например, GMNR.ES очень медленно сходится, а BiCGSTAB начинает стагнировать. Поэтому для решения системы уравнений была разработана специальная итерационная процедура, в основе которой лежит метод BiCGSTAB [14].

4. Результаты численного моделирования

Скорость сходимости предложенной векторной конечно-элементной аппроксимации была исследована на тестовой задаче, имеющей аналитическое решение

Область П = [—0.5, 0.5]3

го1 го1 Е — к2Е = ", к = сопэ^ в П,

div Е = 0,

Е |

дП

(24)

Е

9"

Е"

— 2 * еов(ж) * вт(у) * sin(z) sin(x) * оов(у) * sin(z) sin(x) * sin(y) * cos(y)

Таблица1

Зависимость относительной погрешности приближенного решения в точке (0.1, 0.1, 0.4) от шага сетки при к2 = 1

) 0.1(3630) 0.05(26460) 0.025(201720)

|Е*-Е2| Ет 5.1 ■ 10-5 1.4 ■ 10-5 3.7 ■ 10-6

|Еу-Е0| Е„ 5.1 ■ 10-5 1.5 ■ 10-5 3.5 ■ 10-6

\Ег-ЕЛ Е2 6.1 ■ 10-5 1.7 ■ 10-5 4.8 ■ 10-6

В области П построена последовательность вложенных параллепипидальных равномерных сеток. Полученная после дискретизации задачи (24) СЛАУ решалась с фиксированной точностью 10-10.

Введем следующие обозначения: к — шаг сетки, N — размерность СЛАУ, Е — точное решение, Е^ — приближенное решение.

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

В области П определим следующую задачу:

го^го! Е - к2Е = 0 в П, V _

ее0Е = 0,

Е|Г = Ед,

(25)

где к2 = ^2ее0, е

относительная диэлектрическая проницаемость;

Ед

е-1 * (-2) * еов(ж) * вт(у) * в1п(г) вт(ж) * оов(у) * в1п(г) б1П(х) * вт(у) * еов(г)

(26)

¿(ж, у, г)

1, ж < 0, 3, ж > 0.

На решении задачи (25) было исследовано выполнение условия непрерывности которое в данном случае примет вид

(27) (13),

[п ■ ¿еоЕ]|х=о = 0 или [¿Ех]|х=о = 0.

Таблица2

Зависимость в от шага сетки Н и циклической частоты и в точке (0, -0.11, 0.11), относительная точность решения СЛАУ 10-14

Н(Ж) 0.1(3630) 0.05(26460) 0.025(201720)

и = 103 1.5738-10-2 9.7811-10-1 1.9512

и = 105 1.071210-2 5.317610-3 2.771310-3

и = 107 8.7611 10-3 4.3480-10-3 2.2171-10-3

ТаблицаЗ

Зависимость в от шага сетки Н и циклической частоты и в точке (0, -0.11, 0.11), относительная точность решения СЛАУ 10-19

Н(Ж) 0.1(3630) 0.05(26460) 0.025(201720)

и = 103 1.0699 10-2 5.3854-10-3 2.7273-10-3

и = 105 1.0698-10-2 5.391210-3 2.7642-10-3

и = 107 8.7610-10-3 4.3480-10-3 2.2171-10-3

В табл. 2 и 3 приведены значения относительной погрешности вычисления нормальной компоненты электрической индукции ¿¿0Е^ на границе разрыва относительной диэлектрической проницаемости ¿:

в

- ¿Е^ |ж=о

где Е^+, Е^- — значения электричекого поля (25) слева и справа от границы разрыва соответственно.

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

0.5

0.25

-0.25

-0.5

-0.5

0.5

0.5

0.25

-0.25

-0.5

X

-0.5

0.5

X

0

0

0

0

Рис. 2. ЕХт, Еут в плоскости г = 0.4, частота и = 1 кГц.

Рис. 3. (ЕХт, Е^т) в плоскости г = 0.4, частота и = 1 МГц.

X X

Рис. 4. Егт на прямой г = 0.4, у = —0.5, Рис. 5. Егт на прямой г = 0.4, у = —0.5, ча-

частота ш = 1 кГц: сплошная линия — Е^™, стота ш = 1 МГц: сплошная линия — Е^ штриховая — Еут. штриховая — Еут.

На рис. 2-5 приведены результаты моделирования гармонического электрического поля в области [—1.0,1.0]3, состоящей из двух подобластей (см. рис. 1), со значениями удельной проводимости ае = 1(Ом-м)-1, а2 = 10(0м-м)-1 и относительной диэлектрической проницаемости ве = в2 = 1. В подобласти 1 расположена генераторная катушка, ток в которой равен 1. Границей между подобластями является плоскость х = 0.

Как следует из рис. 4 и 5, нормальная к поверхности разрыва а компонента электрического поля является разрывной, а тангенциальная — непрерывной.

Заключение

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

На основании вычислительных экспериментов можно сделать следующий вывод. Предложенная вариационная постановка позволяет правильно учесть разрывные свойства материалов. Для задач в непроводящих областях с частотами менее 100 кГц использование стандартных решателей не позволяет достаточно точно учесть скачок нормальной компоненты электрического поля Е" или требует для этого трудно достижимой на практике точности решения СЛАУ. Для решения таких задач необходима разработка специальных предобусловливателей, позволяющих учитывать разрывные свойства материалов при низких частотах. Альтернативой является использование вариационной постановки с множителями Лагранжа [11].

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

[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] Trowbridge C.W., Bryant C.F., Emson C.R.I. Some Developments in Electromagnetic Field Computation // Proc. 7th MAFELAP Conf. on The Mathematics of Finite Elements and Appl. 1990.

[8] Nedelec J.C. Mixed Finite Elements in R3 // Numer. Math. 1980. Vol. 35, N 3. P. 315-341.

[9] Nedelec J.C. A New Family of Mixed Finite Elements in R3 // Numer. Math. 1986. Vol. 50, N 1. P. 57-81.

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

[11] Chen Z., Du Q., Zou J. Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients // SIAM J. Numer. Anal. 2000. Vol. 137, N 5. P. 1542-1570.

[12] 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.

[13] 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.

[14] Saad Y. Iterative Methods for Sparse Linear Systems. L.: PWS Publ. Comp., 1996.

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