CC BY
67
Вычислительные технологии
Том 18, № 4, 2013
Анализ векторных конечноэлементных аппроксимаций уравнений Максвелла в анизотропных средах*
Э.П. ШУРИНА1'2, Н. В. ШТАБЕЛЬ2, 1 Новосибирский государственный технический университет, Россия 2Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]
С использованием аппарата дифференциальных форм построена вычислительная схема на базе векторного метода конечных элементов, ориентированная на работу как в изотропных, так и в анизотропных средах. Выполнены расчёты электрического поля в указанных средах с применением неполного и полного векторных базисов первого порядка. Предложен способ интерполяции результатов расчётов между вложенными тетраэдральными сетками, основанный на элементах теории цепей.
Ключевые слова: уравнения Максвелла, дифференциальные формы, анизотропные среды.
Введение
Все увеличивающаяся сложность изучаемых и используемых в технических приложениях физических процессов приводит к необходимости разработки специализированных процедур математического моделирования. Вычислительные процедуры, построенные с использованием дискретных дифференциальных форм, — мощный современный инструмент для численного анализа аппроксимаций дифференциальных уравнений в частных производных.
Идея применения дифференциальных форм к теории электромагнетизма впервые была предложена ВезеЬашрз'ом в 1981 г. Такой подход позволял записать основные законы электромагнетизма (интегральные законы Ампера и Фарадея, закон Гаусса) в инвариантной относительно выбора системы координат форме. Кроме того, требования непрерывности компонент векторного поля (тангенциальных или нормальных) становятся естественными и обоснованными по сравнению с векторной формой записи.
Использование дифференциальных форм в задачах электромагнетизма получило широкое развитие в контексте метода конечных элементов. Основу этому заложил Х. Уитни в своей работе "Геометрическая теория интегрирования", опубликованной в 1957 г. и позднее изданной на русском языке [1]. Книга посвящена аппроксимации дифференциальных форм на сеточных симплициальных разбиениях. Развитие векторного метода конечных элементов связано с работами V. Girault и P.A. Raviart, F. Brezzi и M. Fortin, P.G. Giarlet и других авторов.
* Работа выполнена при поддержке Интеграционного проекта СО РАН № 98.
E. Tonti предложено приближение электромагнитного поля, основанное на использовании глобальных (интегральных) величин [2], позволяющих без привлечения дифференциальных уравнений получить так называемую конечную формулировку, являющуюся естественным расширением теории цепей применительно к электромагнитному полю. В 1988 г. A. Bossavit предложил расширение метода конечных элементов — векторный метод конечных элементов, который гармонично совмещал использование внешних дифференциальных форм и глобальные переменные поля [3].
Теория дискретных дифференциальных форм часто применяется в публикациях современных авторов (R. Hiptmair [4], F.L. Teixeira [5], D. Arnold [6]) для анализа конеч-ноэлементных аналогов уравнений Максвелла и получения оценок сходимости конеч-ноэлементных аппроксимаций, степеней свободы, свойств нуль-ядра оператора rot-rot. Исследования либо носят теоретический характер, либо расчёты и оценки представлены для задач не на тетраэдральных, а на регулярных параллелепипедальных разбиениях. При расчётах электромагнитных полей в анизотропных средах с тензорными коэффициентами в основном используют параллелепипедальные разбиения, поскольку базисные функции на таких сетках параллельны осям тензоров коэффициентов электропроводности, диэлектрической и магнитной проницаемости.
В настоящей работе представлена вычислительная схема для решения трёхмерного уравнения Гельмгольца, построенная с использованием теории дискретных дифференциальных форм. Предложенная схема позволяет решать уравнение Гельмгольца как в изотропных, так и в анизотропных средах на тетраэдральных неструктурированных разбиениях. Для получения оценок сходимости на тетраэдральных разбиениях в анизотропных средах, где точное решение неизвестно, предложен способ интерполяции поля между двумя сеточными разбиениями на основе теории цепей.
1. Система уравнений Максвелла в терминах дифференциальных форм
Интегральные законы Фарадея, Ампера и Гаусса, образующие систему уравнений Максвелла, в дифференциальных формах имеют вид
S
(1)
<
S
S
S
/
B = 0,
S
где S — поверхность, l — замкнутый контур, V — объём.
На языке дифференциальных форм в трёхмерном пространстве уравнения Максвелла в частотной области (при гармонической зависимости от времени eiwt) могут быть записаны в виде
dE = -iwB,
dH = iwD + J + Js,
(2)
dD = p, dB = 0,
где E и H — напряжённости электрического и магнитного полей, являющиеся 1-форма-ми (дифференциальные формы первого порядка); D и B —электрическая и магнитная индукция (2-формы); J, Js — плотности электрического тока (2-формы); p — 3-форма плотности заряда; d — оператор внешней производной, являющийся отображением k-формы в (k + 1)-форму [7]. Применительно к дифференциальным формам в метрике простраства R3 для 1-, 2-, 3-форм оператор d играет роль grad, rot и div операторов векторного исчисления.
Уравнения (2) не зависят от выбора R3, все метрические характеристики учитываются в дополнительных уравнениях, называемых материальными соотношениями:
D = E,
B = H,
J = *a E, (5)
где * — оператор Ходжа, представляющий собой отображение к-формы в (3 — к)-форму для трёхмерного пространства; е, ^ — диэлектрическая и магнитная проницаемость; а — тензор электропроводности. В общем случае е и ^ могут быть тензорами второго ранга с соответствующей мерой. Оператор Ходжа содержит в себе информацию о метрических характеристиках пространства, в котором определены уравнения Максвелла. Форма записи *еЕ означает воздействие тензора е на форму *Е [5]. Уравнения (2)-(5) представляют собой систему уравнений Максвелла.
Используя представления (3)-(5), можно перейти от системы уравнений Максвелла к уравнению Гельмгольца
6 dE — и2 *е Е + 1и*а Е = —{иЗ3. (6)
Интегральная форма уравнений Максвелла (1) позволяет естественным образом определить условия непрерывности поля на внутренних границах, разделяющих подобласти с разными свойствами:
[ Е = 0, [и = ! Зг, [ Б = ! рг, [в = 0, (7)
где I — замкнутый контур, проходящий перпендикулярно поверхности, разделяющей подобласти с разными свойствами, частично принадлежащий той или другой подобласти; Б — площадь поверхности, разделяющая подобласти с различными свойствами; Зг — наведённый ток на границе двух подобластей Г; рг — поверхностные заряды на границе подобластей.
2. Дифференциальные формы на сетках
Рассмотрим дискретные дифференциальные формы на симплициальных сеточных разбиениях. Симплициальную сетку можно представить в виде прямой суммы ячеек к-го порядка или к-ячеек вк, полностью заполняющих область в Я3 без наложения. Множество всех ячеек порядка к обозначим хк, тогда сеточное разбиение есть комплекс ячеек
п
х = 0 хк.
к=0
Для симплициальных тетраэдральных разбиений справедливо следующее соответствие: — точка, вершина тетраэдра; 51 — линия, ребро тетраэдра, в2 — плоскость, грань тетраэдра, в3 — объём тетраэдра. При заданной ориентации ячейки одного типа образуют к-цепь, определённую как линейная комбинация к-ячеек в \к:
Бк = ^ оцвк е хк, (8)
г
где аг — действительные целые числа.
Для дифференциальной формы П, интегрируемой по области 7, можно определить её сеточный аналог в терминах к-цепей Бк и сеточных форм вк
^П ^<Бк, 0к >. (9)
7
Тогда можно определить базис сеточных форм Ок, соответствующий базису к-цепей, так что < Бк, 0к >= , где — символ Кронекера. Таким образом, сеточная форма 0к определяется как
0к = Е вгОк, (10)
г
где веса вг е С, С — абелева группа [8].
Оператор внешней производной й определен на к-цепи следующим образом [9]:
< Бк ,й0к-1 >=< дБк, 0к-1 >, (11)
где дБк — граница цепи Бк — (к — 1)-цепь Бк-1.
3. Система уравнений Максвелла на сетках
Уравнения (2) в терминах дискретных дифференциальных форм принимают вид
< Б2,¿Е > = —г^<Б2,В>, (12) <Б2,йН> = 11и<Б2,В> + <£2,/> + <Б2^3 >, (13)
< Б3, ¿В > = 0, < Б3, ¿Б >=< Б3,р >,
где Б2, Б3 — соответственно 2-, 3-цепи тетраэдральной сетки, ассоциированные с формами Е и В, ячейки цепей Б2, Б3 называют первичной сеткой, бБ2, Б — цепи, образованные ячейками дуальной сетки, ассоциированной с формами Н и Б.
Комплексы ячеек первичной и дуальной сеток связаны соотношением: к-симплекс первичной сетки соответствует (3—к)-симплексу дуальной сетки. В результате получаем связи вида: точка — объём, ребро — грань, грань — ребро, объём — точка.
К построению дуального комплекса ячеек можно подходить по-разному. Некоторые авторы [10] строят дуальную сетку на основе барицентров симплексов, другие [11] предлагают использовать в качестве точек дуальной сетки центры описанных окружностей для симплексов как точки, равноудалённые от всех (к + 1)-точек к-симплекса. В любом случае ячейки дуальной сетки не являются симплексами, а есть выпуклые многогранники. Соотношения (3)-(5) в дискретном виде записываются как
<~Б2,Б>= <в),Е>,
3
< >= ^Ыгз <в) ,Н >,
<Б2,3>= £[*„]гз <в),Е>, (14)
3
где \ke\ij, [*и]з, ]з — элементы матриц Ходжа, [*£] : хк ^ Хк, Хк — дуальный комплекс ячеек.
Перейдем к уравнению Гельмгольца. Подействуем оператором *и-1 на обе части уравнения (12)
<Б\-к„-1 СЕ>=-гь < Б1,Н > . (15)
По определению оператора Ходжа *и-1 с1Е — 1-форма, определённая на дуальной сетке.
Возьмём внешний дифференциал от уравнения (15) и в соответствии с (11) получим
< Б2,С*и-1 СЕ >= -гь < Б2,СН > . Тогда, подставив в это выражение уравнение (13), имеем
< Б2,С*и-1 СЕ > -ь2 < Б2, Б > +гь < Б2, 7 >= -гь < Б2, 73 >,
или с учётом (14) —
< Б2,С*- СЕ > -ь2 < Б \*еЕ > +гь < Б Е >= -гь < Б2,,13 > . (16)
Рассмотрим выражение < Б2,С*и-1 СЕ > (**). По свойствам внешнего дифференциала
< Б2,С*и-1 СЕ >=< дБ2,*^СЕ >,
где дБ2 — граница 2-цепи дуальной сетки, т. е. 1-цепь дуальной сетки
дв1 = ^2 в3. 3
Отметим, что первый оператор С в (**) действует на форму, определённую на дуальной сетке. В работах [10, 12] показано, что выражение < Б2,С*и-1 СЕ > идентично матрице жёсткости, используемой в векторном методе конечных элементов для аппроксимации на рёбрах первичной сетки Б .
Дифференциальная форма З3 определена на рёбрах дуальной сетки как 2-форма с внешней ориентацией, т. е. по нормали к плоскости поперечного сечения источника
тока — петли. Как правило, в методе конечных элементов при решении задач электромагнетизма, в частности задач геоэлектрики, объём источника не учитывается. Источник тока индукционного типа представляют бесконечно тонким кольцом, по которому распределен ток = я*, где — 1-форма с внутренней ориентацией вдоль ребра первичной сетки. Более подробно ориентация дифференциальных форм на первичной и дуальной сетках рассмотрена в [13].
С учётом вышесказанного уравнение (16) принимает вид
< Б¿Е > —< Б\*£Е > < Б 1,*аЕ >= —г^ < Б1,3я* > . (17)
Чтобы найти конечноэлементную аппроксимацию Е, необходимо построить базис к-форм П,к — аналог базиса к-коцепей . Базис должен обладать следующими свойствами:
1) П,к = 0 вне вк и его соседей;
2) Пяк = ¿Па„£, согласуется с (11) и показывает отношения между Пяк для разных к.
^^ г г г
Построение базисов основано на этом пункте;
3) сеточная форма Огк и её непрерывный аналог П,к дают одинаковый результат интегрирования по соответствующим элементам х и X, т. е.
/ П4 =< вк,Ок >,
7
где 7 — к-мерная область в X, соответствующая к-ячейке в х.
Для симплициальных сеток базис, удовлетворяющий условиям 1-3, был предложен в работах [1, 14, 15], и базисные функции получили название уитни-форм. уитни-форма, ассоциированная с к-ячейкой вк, имеет вид
к
Пак = к! 1)-7Л^¿Лг0 Л ■ ■ ■ Л ¿Лг^-1 Л ¿Лу+1 Л ■ ■ ■ Л ¿Лгк,
г ¿=0
где Л^, 0 ^ ] ^ к, — барицентрические координаты симплекса вк, Л — внешнее произведение.
Для 1-формы Е базис на каждом ребре принимает вид привычной векторной базисной функции, ассоциированной с ребром тетраэдра:
Пя1 ^ ^ 1 = Лг°йЛг1 — Лг1 ¿Лг°, или = Л;°УЛл — Лг^ЛгО. (18)
] г г
Для задач электромагнетизма аппроксимация величины ¿Е, или го1Е, так же важна, как и аппроксимация самого поля. Известно, что поле можно представить в виде суммы двух частей: вихревой и градиентной. Базис (18) описывает вихревую часть поля, поскольку внешний дифференциал от уитни-функции не обращается в ноль. Для получения полного базиса, описывающего и роторную, и градиентную часть поля, базис (18) необходимо дополнить в соответствии со свойством 2 функциями, которые являлись бы градиентами сеточных форм меньшего порядка Пд,к. Так как уитни-формы определены на рёбрах в\, то двк есть узлы сетки и, значит, Пд, к есть полиномы от Лк-форм, определённых в узлах сетки. Тогда функции, дополняющие базис (18) до полного базиса первого порядка, имеют вид
С,1 = ¿(Лг°Лг1) = Аг0¿Агl + Лгl¿Лг0, или С,1 = Лг°УЛг1 + Лг1 УЛг°. (19)
Для полного базиса первого порядка на каждом ребре сетки определены две базисные функции: одна вида (18), вторая вида (19), соответствующие роторной и градиентной частям поля.
Используя уитни-формы, можно записать аппроксимацию поля Е на 1-цепи как сумму векторных базисных функций по всем симплексам б"1:
Е = ^ . (20)
аг -1 |т,-¿т,1 Л т,1 < вКе > —
Умножив (17) на т, 1, получим систему уравнений
^ а ^ ¿[*«-1 ]т, ¿т ,
г т
—т2 ^ а, ^ [*е]т, т,т Л т,1 < в1, е >+
г т
+гт ^ а, ^ Л < в1, е > = -гт ^ а, < в1, 7* « Л >. (21)
г т г
В терминах векторного исчисления слагаемое ^ а, ^ ¿[*«-1 ]т, Л т, 1 < в1, е > соот-
гт
ветствует матрице жёсткости, обычной для векторных конечноэлементных постановок [10, 16]. Выражение ^ а, ^[*е]т,Л т, 1 < в1, е > — конечноэлементная матрица
г т 3
массы.
Параметры среды в уравнениях Максвелла в общем случае — положительно определённые тензоры второго ранга. В настоящей работе рассматриваются изотропная среда по магнитной и диэлектрической проницаемости и анизотропная среда по электропроводности.
Электропроводность анизотропной среды может быть задана в виде билинейной формы [17]
а = аже! ё1 + ау 6262 + аг ёэёэ, (22)
где ах, ау, а^ — константы, ёг — единичные базисные векторы главных осей тензора. В этом случае матрица массы принимает вид
М, = У ■ а(к) ■ ^, (23)
где а(к) — значение тензора на элементе, — уитни-функция. Будем полагать, что тензор а не изменяется внутри элемента. С учётом представления тензора в виде билинейной формы уравнение (23) запишем как
Г м Г
М,- = ал(ё1 ■ Шг)(ё1 ■ + ^ / ау(б2 ■ Шг)(б2 ■ +
Пк Пк
„
+ ^ I аг(ёз ■ Шг)(ёз ■ . (24)
А:=Ь
Пк
Здесь выражение ■ — проекция базисной функции на соответствующую ось
Тензоры вида (22) обычно описывают трансверсально-изотропные среды с выделением одного из направлений. В общем случае произвольная анизотропная среда представляется билинейной формой
о = Рххёхёх + аХу е1е2 + ахх е^з + ОуХё2е,1 + Оуу (§2(52 + Оух §2ёз+ +Схх£ъ£\ + Оху ёзё2 + о гх езёз.
Для таких сред (23) выписывается в виде
(25)
м=УОхх(ё1 ■^)(ё1 ■)ЛПк+УОху(ё1 ■^ )(ё2 ■ ^№+
Пк Пк
^У Охх (§1 ■ Wi )(ёз ■ Wj )(10,к + J о ух (¿2 ■ Wi)(ël ■ Wj +
Пк Пк
+ / ОУУ(ё. ■ ■ ^№ / °ух(е2 ■ ^)(ё3 ■ № +
Пк Пк
+ /°хх(ёз ■ ^)(ё1 ■ ^№ + /°ху(ез ■ ^ ■ №+
Пк Пк
+ У Охх(ёз- Wi)(ëз• Wj.
Пк
(26)
4. Вычислительный эксперимент
Для получения оценок погрешности вычислительной схемы было проведено две серии вычислительных экспериментов. Первая серия — расчёты на структурированных вложенных сетках для двух типов базисных функций (полный и неполный базисы первого порядка) в области с изотропными свойствами. Вторая серия экспериментов проводилась на тех же сетках и базисах в анизотропном полупространстве.
Расчётная область — куб со стороной 2 м. Вложенные тетраэдральные сетки строились путём разбиения каждой стороны расчётной области на 10/20/40/80 частей с последующим делением каждого получившегося куба на шесть тетраэдров. Таким образом, тетраэдр грубой сетки состоит из восьми тетраэдров мелкой сетки. Размерности полученных сеточных разбиений приведены в табл. 1. Базисные функции, использованные в расчётах: 1 — уитни-функции первого порядка вида (18), 2 — полный базис первого порядка вида (18)-(19).
Таблица 1. Размерности сеточных разбиений
Сетка Число Число неизвестных Число неизвестных
тетраэдров для неполного базиса для полного базиса
10 6000 7930 15 860
20 48 000 59 660 119 320
40 384 000 462 520 925 040
80 3 072 000 3 641840 7283 680
Оценивание вычислительных схем на базе векторного метода конечных элементов для уравнений Гельмгольца представляет большой интерес для исследований. В работе [18] предложены теоретические оценки погрешности вычислительной схемы для уравнения Гельмгольца, аппроксимированного с помощью полного векторного базиса первого порядка. В норме пространства Ь2 оценка погрешности имеет вид
||Е - Ен\\0 ^ С 1П ||Е - УьЦо + к ||У х (Е - ^)||(
\VhGVh
где || ■ ||о — норма в пространстве Ь2; Е — точное решение; Е^ — численное решение; у^ — векторная базисная функция вида (18)-(19); Vh — пространство всех базисных функций; к, С — константы, зависящие от параметров среды и сеточного разбиения.
Для задач с известным аналитическим решением при оценивании погрешности в качестве Е, как правило, используют интерполяцию аналитической функции на сеточное разбиение. Сложность оценивания точности конечноэлементного решения для уравнения Гельмгольца с локальным источником заключается в том, что в отличие от задач с аналитическим решением точное решение в общем случае неизвестно, особенно если задача решается в областях с анизотропными свойствами. В настоящей работе предложены оценки погрешности, основанные на сравнении двух численных решений, полученных для различных параметров аппроксимации задачи (измельчение сетки и изменение базисных функций).
Для вычисления оценок погрешности для аппроксимации на разных базисах использовалась следующая формула:
Е ■= '> 'Т-Г р (27)
1-2
где || ■ || — евклидова норма; Е/ — значение поля Е в точке г, полученное на неполном
^2
"г1
базисе первого порядка; Ег2 — значение поля Е в точке г, полученное на полном базисе первого порядка для той же сетки; р — количество точек, используемых в расчёте. Для уитни-базиса применялась следующая форма оценки:
ЦЕС_—Е11| ||Е с
Е-/ ", (28)
где || ■ || — евклидова норма; Ес — вектор весов базисных функций на грубой сетке; ЕI — вектор весов базисных функций, интерполированных с мелкой сетки на грубую.
Аппроксимация области моделирования тетраэдральными разбиениями в изотропных и анизотропных средах при расчётах трёхмерных электромагнитных полей приводит к необходимости разработки специальных процедур интерполяции при вычислениях на вложенных регулярных или нерегулярных (неструктурированных) сетках.
Интерполяция с мелкой вложенной сетки на грубую для уитни-форм естественным образом вытекает из теории цепей (рис. 1). Действительно, аппроксимация на грубой сетке есть
<£С\ЕС >= ^ агрг <з1,вг >,
г
где аг — коэффициенты, отвечающие за ориентацию ребер в сетке; вг — веса базисных функций; 9г — сеточная форма — уитни-форма. По построению сетки ориентация всегда положительна, и коэффициенты аг можно опустить.
53
Рис. 1. Соотношения грубой и мелкой сеток
С другой стороны,
<£),Е' >= ^в <§1Л >, (29)
г
где 51 — ребра мелкой сетки. Заметим, что по построению сетки в1 = + Тогда, сгруппировав элементы (29), получим
< £},Е/ >= ^ < ^ ,вг0^0 + > .
г
Таким образом, в = вг0 + вг1 — формула для интерполяции весов уитни-форм с мелкой вложенной сетки на более грубую.
Для верификации предложенного способа интерполяции решения выполнены расчёты для задачи, имеющей аналитическое решение, на вложенных сетках. Было решено уравнение (3) с единичными коэффициентами среды для следующего аналитического решения и соответствующей правой части:
(х \ / (—и2 + ги)х \
у I , Я = I (-и2 + ш)у I . -2* / \ (2и2 - ¿2и)г + 2 у
В табл. 2 приведены оценки погрешности вычислений для задачи с аналитическим решением. Погрешность — оценка вида (28), где Ес — численное решение, Е ^ — интерполяция точного аналитического решения на ту же сетку. В последнем столбце табл. 2 приведена оценка погрешности £с_/ вида (28) на двух вложенных сетках
Таблица 2. Оценки погрешности для задачи с аналитическим решением относительно точного решения и при интерполяции с мелкой сетки на грубую £с_/
Сетка Сетка £с_/
10 3.9953Е-04 — —
20 4.7301Е-05 10-20 3.634Е-04
40 2.7431Е-05 20-40 3.6777Е-05
с применением сеточной интерполяции. При измельчении сетки 10 в два раза (сетка 20) погрешность Ес— уменьшается почти на порядок. Дальнейшее измельчение сетки позволяет повысить точность решения еще на 72 %. Уменьшение погрешности почти в 2 раза соответствует первому порядку сходимости базисных функций, как это было показано в [14, 15]. Из оценки £c-f с использованием сеточной интерполяции между парами вложенных сеток следует, что погрешность, вычисленная таким образом, не превышает соответствующую погрешность £с_4, полученную на более грубой сетке. Скорость сходимости, полученная по оценкам Ес— на сетках 10 и 20, имеет тот же порядок, что и вычисленная по оценкам интерполированной ошибки на парах сеток 10/20 и 20/40 (уменьшение в 8.4 и 9.8 раз соответственно). Размерности систем линейных алгебраических уравнений для использованных сеток приведены выше в табл. 1.
Для получения оценок аппроксимации поля Е была выбрана следующая модельная задача: расчётная область — изотропное или анизотропное полупространство, на границе воздух — среда находится токовый источник — квадратная петля, вся область — куб со стороной 2 м, сторона петли 0.2 м, частота, задаваемая в гармонической зависимости от времени, составляет 10 кГц, электропроводность изотропной среды 0.1 См/м. Электропроводность анизотропной среды задана тензором
/ 0.5118967549 -0.07529063325 0.02656734376 а = I -0.07529063325 0.1034841068 -0.005439385875 0.02656734376 -0.005439385875 0.1022944313
Представление результатов моделирования трёхмерных электромагнитных полей, полученных на тетраэдральных разбиениях, требует использования специальных процедур интерполяции, позволяющих перейти от векторного представления поля, ассоциированного с рёбрами сетки, к покомпонентному виду Ех, Еу, Ех в сечении или по профилю. На рис. 2 приведено распределение компоненты поля Ех в сечении г = 0,
У, м 2 Г
-I_1_I_I_I_1_и
J
X, М
X, м
Рис. 2. Компонента поля Ех в изотропной (а) и анизотропной (б) среде
а
полученное при расчётах в изотропной и анизотропной среде. Из рисунка видно влияние тензорного коэффициента электропроводности на электрическое поле. В анизотропной среде картина поля теряет свою симметричность, линии равного уровня деформируются по сравнению с картиной поля в изотропной среде.
В табл. 3 представлены оценки погрешности вычислений поля Е вида (28) на последовательности вложенных структурированных сеток в изотропной и анизотропной среде. Для задачи с изотропными коэффициентами среды измельчение сетки позволяет повысить точность решения на 48%, что соответствует теории. В анизотропной среде измельчение сетки в два раза позволило уменьшить погрешность на 31%, или в 1.416 раза. Решение уравнения Гельмгольца в анизотропной области более чувствительно к выбору сеточного разбиения и требует использования более мелких сеток для достижения большей точности по сравнению с задачей в изотропной среде.
Таблица 3. Оценки погрешности вычислений поля Е в изотропном и анизотропном полупространстве с использованием интерполяции
Параметр Изотропная среда Анизотропная среда
¿20—40 0.103 0.179
¿40—80 0.052 0.123
Таблица 4. Оценки погрешности вычислений поля Е в изотропном и анизотропном полупространстве с использованием интерполяции для базисов 1 и 2
Параметр Изотропная среда Анизотропная среда
Е1—2 на сетке 20 0.0121 0.138
¿4—2 на сетке 40 0.00676 0.0702
¿4—2 на сетке 80 0.00229 0.024
Рис. 3. Мнимые компоненты поля Ех и Еу в изотропной (а) и анизотропной (б) среде на структурированной (Ех стр, Еу стр) и неструктурированной (Ехнестр, Еунестр) сетках по профилю у = -0.6 м, г = 0
В табл. 4 показаны оценки погрешности вычисления электрического поля с использованием базисов двух типов: 1 — неполный базис уитни-функций вида (18), 2 — полный векторный базис первого порядка вида (18)-(19). Использование полного базиса для задачи в изотропной среде для одной и той же сетки позволяет улучшить решение во втором знаке для сетки 10. Дальнейшее измельчение сетки при увеличении базиса до полного даёт уменьшение оценки примерно в два раза с увеличением размерности сетки (см. табл. 1,3). Оценки для решения в анизотропной среде имеют тот же уровень сходимости (уменьшение примерно в два раза), однако оценка погрешности решения для самой грубой сетки в случае использования полного базиса меньше, чем при простом измельчении сетки (см. табл. 2, 3). При решении задач в анизотропных средах оптимальным выбором будет использование полного базиса совместно с дроблением сетки.
Приведённые выше результаты были получены на тетраэдральных структурированных сетках, однако при решении практических задач чаще всего используют тетраэдральные неструктурированные сетки с локальными сгущениями вблизи источника поля. На рис. 3 представлены графики мнимых компонент электрического поля Ex, Ey по профилю y = -0.6 м, z = 0 для изотропной и анизотропной среды на структурированной сетке 80 и неструктурированной тетраэдральной сетке. Неструктурированная сетка имеет 481 741 тетраэдр и 569 364 ребра со сгущением в области источника. Погрешность вычислений на структурированной и неструктурированной сетках не превышает 1.2 %, но вычислительные затраты для трёхмерного электрического поля на неструктурированных сетках значительно ниже.
Таким образом, с использованием аппарата дифференциальных форм построена вычислительная схема для моделирования трёхмерного электрического поля от индукционного источника, позволяющая работать как в изотропных, так и в анизотропных средах. Представление сеточного разбиения в виде комплекса ячеек и объединения k-цепей различного уровня позволяет интерполировать конечноэлементное решение с одной сетки на другую. Сеточная интерполяция на тетраэдральных сетках была применена для получения оценок погрешности для неполного базиса первого порядка (уитни-функ-ции). Использование полного базиса первого порядка для аппроксимации дифференциальных 1-форм по сравнению с неполным базисом первого порядка позволяет достичь более точного решения.
Список литературы
[1] Хасслер Уитни. Геометрическая теория интегрирования. М.: ИЛ, 1960.
[2] Tonti E. On the geometrical structure of electromagnetism // Gravitation. Electromagnetism and Geometrical Structures. Bologna: Pitagora Edit., 1995. P. 281-308.
[3] Bossavit A. Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism // IEEE Proc. 1988. 135. Pt. A. 8. P. 493-500.
[4] Hiptmair R. From E to edge elements // Berichte des IZWR. 2003. No. 1. P. 39-47.
[5] Teixeira F.L. Geometric aspects of the simplicial discretization of Maxwell's equations // Progress in Electromagnetics Res. 2001. Vol. 32. P. 171-188.
[6] Douglas N. Arnold, Richard S. Falk, Ragnar Winther. Finite element exterior calculus: from Hodge theory to numerical stability // Bulletin (New Series) of the Amer. Math. Soc. 2010. Vol. 47, No. 2. P. 281-354.
[7] Bossavit A. Differential geometry for the students of numerical methods in electromag-netism // Electricite de France. Etudes et Recherchers. August, 1991.
[8] Винберг Э.Б. Курс алгебры. 3-е изд. М.: Факториал Пресс, 2002. 544 с.
[9] Mathieu Desbrun, Eva Kanso, Yiying Tong. Discrete differential forms for computational modeling // Eds. A.I. Bobenko, P. Schroder, J.M. Sullivan and G.M. Ziegler. Oberwolfach Seminars. Discrete Differential Geometry. Birkhaiiser Verlag Basel/Switzerland, 2008. Vol. 38. P. 287-323.
[10] Bo He, Teixeira F.L. Geometric finite element discretization of Maxwell equations in primal and dual spaces // Phys. Lett. A. 2006. Vol. 349. P. 1-14.
[11] Mathieu Desbrun, Anil N. Hirani, Melvin Leok, and Jerrold E. Marsden. Discrete Exterior Calculus // http://arxiv.org/pdf/math/0508341.pdf
[12] Bo He, Teixeira F.L. Differential forms, galerkin duality, and sparse inverse approximations in finite element solutions of maxwell equations // IEEE Trans Antennas Propag. 2007. Vol. 55, No. 5.
[13] Teixeira F.L., Chew W.C. Lattice electromagnetic theory from a topological viewpoint // J. Math. Phys. 1999. Vol. 40, No. 1.
[14] Nedelec J.C. Mixed finite elements in R3 // Numer. Math. 1980. Vol. 35, No. 3. P. 315-341.
[15] Nedelec J.C. A new family of mixed finite elements in R3 // Ibid. 1986. Vol. 50, No. 1. P. 57-81.
[16] Burkay Donderici, Teixeira F.L. Mixed finite-element time-domain method for transient maxwell equations in doubly dispersive media // IEEE Trans. Microw. Theory Tech. 2008. Vol. 56, No. 1. P. 113-120.
[17] Абрамов А.А. Введение в тензорный анализ и риманову геометрию. М.: Физматлит, 2004.
[18] Zhong L.Q., Shu S., Wittum G., Xu J.C. Optimal error estimates for Nedelec edge elements for time-harmonic Maxwell's equation //J. Comp. Math. 2009. Vol. 27, No. 5. P. 563-572.
Поступила в редакцию 27 февраля 2013 г., с доработки —11 июня 2013 г.
