CC BY
20
Вычислительные технологии
Том 19, № 6, 2014
Анализ вычислительных схем для моделирования электромагнитного поля в средах с контрастными включениями в широком диапазоне частот
М.И. Эпов1, Э.П. Шурина1'2, Е.И. Михайлова1'2
Задача моделирования трёхмерного электромагнитного поля в областях с малыми включениями в частотном диапазоне возникает во многих приложениях физики, биологии, материаловедения и т.д. Однако из-за сложной геометрии, наличия внутренних границ, высокой контрастности используемых материалов решать такие задачи аналитическими методами невозможно. В статье рассматриваются современные модификации неконформного метода Галеркина, ориентированные на моделирование электромагнитных полей в областях с включениями, размеры которых значительно меньше габаритных размеров расчётной области, в диапазоне частот от 100 кГц до 1 ГГц.
Ключевые слова: векторный метод конечных элементов, разрывный метод Галеркина (DG-метод), композитный материал, микровключения.
Введение
Искусственные и естественные материалы, имеющие сложную геометрическую структуру и разрывные электрофизические свойства, в настоящее время находят широкое применение во многих областях промышленности и вызывают повышенный интерес исследователей. Электрофизические характеристики таких гетерогенных сред в целом могут значительно отличаться от характеристик материалов, из которых они состоят. По структуре гетерогенные среды можно разделить на двухкомпонентные (матрица/включения) и многокомпонентные [1]. Включения могут иметь различную геометрическую форму: сферические, пластинки-резонаторы, спирали, разорванные кольца и т. д. [2-5], и располагаться в матрице как регулярно, так и хаотично.
В силу сложной геометрии, высокой контрастности электрофизических характеристик матрицы и включений, рабочего диапазона частот от радиочастот до терагерце-вых, а также взаимодействий между включениями электромагнитные поля, наведённые в таких материалах, обладают сложной конфигурацией. Круг решаемых задач электромагнетизма для композитных структур достаточно широк: влияние формы включений и частоты на формирование волнового процесса [2, 3], нахождение эффективных электрофизических характеристик сложной среды [1-8] и др.
Классификация задач выполняется по различным признакам, например по частотной характеристике, рассматриваемой как отношение длины волны к геометрическим
1 Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, Новосибирск, Россия
2 Новосибирский государственный технический университет, Россия
Контактный e-mail: [email protected]
размерам образца. Различают высокие частоты (длина волны много меньше геометрических размеров области моделирования) [4], средние частоты (длина волны соизмерима с геометрическими размерами области моделирования) [9] и низкие частоты (длина волны много больше геометрических размеров области моделирования) [8]. Длину волны можно также соотносить с размерами включений [1].
Для моделирования электромагнитных полей в трёхмерных областях одним из наиболее эффективных является векторный метод конечных элементов (ВМКЭ) [10, 11]. В средах с разномасштабными включениями необходимо строить адаптивные разбиения расчётной области со значительным сгущением сеток во включениях, что приводит к резкому росту размерности получаемой системы линейных алгебраических уравнений (СЛАУ), особенно на базисах высоких порядков [12]. При повышении частоты, когда на одну длины волны должно приходиться минимум девять тетраэдров [13], этот рост ещё более значителен. Соответственно время решения СЛАУ и необходимый объём памяти увеличиваются.
Один из возможных путей решения данной проблемы — распараллеливание вычислительного процесса. Однако распараллеливание стандартных схем конформного ВМКЭ не позволяет эффективно преодолеть указанные выше проблемы. Возникает необходимость в построении вычислительных алгоритмов, обладающих естественным параллелизмом, что в свою очередь приводит к развитию неконформных методов: разрывный метод Галеркина (DG-метод) [14- 16], метод мортаров [17], метод декомпозиции области [18, 19].
Разрывный метод Галеркина, зарекомендовавший себя для моделирования течений [20, 21], движения теплового фронта [22, 23] и т.д., последние 10 лет применяется для широкого круга задач электромагнетизма. DG-метод позволяет достаточно просто комбинировать в одной сетке элементы разной геометрии и элементы разных порядков. Наиболее хорошо изученными являются схемы для системы уравнений Максвелла первого порядка во временной области [24, 25]. Для данных постановок наряду со сложностями, связанными с определением лифтинг-операторов (скачка и осреднения), возникают также проблемы с выбором оптимальной схемы аппроксимации по времени [25].
Вычислительных схем DG-метода для решения задач электромагнетизма в гармоническом режиме немного. Выделяют постановки для уравнений второго порядка относительно векторной комплекснозначной величины напряжённости электрического/магнитного поля (уравнение Гельмгольца) [25-28] и для системы уравнений первого порядка относительно векторных комплекснозначных величин напряжённости электрического и магнитного поля [29, 30].
При решении уравнения Гельмгольца отдельно рассматриваются высокие и низкие частоты, для которых используются разные вариационные формулировки: в [26, 28] для высоких частот (от 100 МГц), в [27] — для низких частот. В работе [27] предлагаются два подхода: в первом дивергентное ограничение учитывается с помощью регуляри-зационного члена, во втором применяется смешанный подход и вводятся множители Лагранжа.
Для данного класса задач более распространена IP-схема (внутреннего пенальти) [26-28], но могут быть использованы и другие, например постановка Бреззи [29].
Для гармонического режима в [27] приведены возможные вариационные формулировки и априорные оценки, но не показаны результаты численного решения. В работах [25, 31] численное моделирование выполняется для модельных задач, имеющих гладкое аналитическое решение. В [32] решается уравнение Гельмгольца для ^ = 1, е = 1.
Проведённые в [25, 31, 32] численные эксперименты позволяют верифицировать используемые схемы, однако не гарантируют переносимость этих схем на задачи с реальными физическими коэффициентами.
Целью настоящей работы является сравнительный анализ IP-схем разрывного метода Галеркина для моделирования электромагнитных полей в средах с сильно контрастными микровключениями в частотном диапазоне от 100 кГц до 1 ГГц. Сравнивается эффективность применения IIP-, NIP-, SIP-постановок и IP-схем без осреднений на низких, средних и высоких частотах.
1. Постановка задачи
Электрическое поле в гармоническом режиме описывается уравнением Гельмгольца
Vx/ VxE+k2E = 0 в ft,
(1)
где k2 = iua — u2e — волновое число, u = 2nf — циклическая частота, Гц, a — электропроводность среды, См/м, e = ere0 — диэлектрическая проницаемость, Ф/м, er — относительная диэлектрическая проницаемость, e0 = 8.85 x 10-12 Ф/м, / = /г/0 — магнитная проницаемость, Гн/м, /лг — относительная магнитная проницаемость, = 4п x 10-7, Гн/м, дft = Гт U Ге U Г0 — трёхмерная область с липшиц-непрерывной границей Гт = {у = 0} U {у = 0.07}, Го = {z = 0} U {z = 0.07} U {x = 0.07}, Ге = {x = 0} (рис. 1).
На границах расчётной области заданы краевые условия двух типов: • электрические
n x Е|Го = 0, (2)
n x EL = Ео, (3)
магнитные
Л 1 VxE x n|r = —iug.
(4)
Область моделирования ft неоднородна по электрофизическим характеристикам:
(г) (подобласти-вклю-
ftматрица (подобласть-матрица или вмещающая среда) и ft
включение (
а
Рис. 1. Области моделирования ft: а — сферические включения (объёмная доля включений меньше 1%), б — параллелепипеидальные включения (объёмная доля включений 17.5%)
чения) имеют различные значения диэлектрической проницаемости е, электропроводности а и магнитной проницаемости На межфрагментарных границах Г^ (матрица-включение) выполняются следующие условия [33]:
[n х E]r.j = 0, (5)
[n ■ (а + гше)Е]г.. = 0. (6)
Условие непрерывности (6) обеспечивает выполнение закона сохранения заряда в слабом смысле [33].
2. Вариационная постановка
2.1. Конформный векторный метод конечных элементов
Введём функциональные пространства [10, 11]
H(rot, П) = {u G L2(Q)| Vxu G L2(n)}, (7)
Ho(rot, П) = { u G H(rot, П)| u x n|dQ = 0}. (8)
В пространствах (7)-(8) определены норма и скалярное произведение.
Вариационная постановка принимает вид: найти Е G Ho (rot, П) + Е0 такое, что Vv G H0(rot, П) выполняется
У VxE ■ Vxvdn + j k2E ■ vdn = -iu J g ■ vdS. (9)
Q Q Гт
2.2. Неконформный разрывный метод Галеркина (DG-метод)
Рассматриваются IP-схемы разрывного метода Галеркина и их модификации для уравнения Гельмгольца (1) с краевыми условиями (2)-(4). Производится анализ вариационных формулировок на высоких частотах [26, 28].
Решение поставленной задачи разрывным методом Галеркина осуществляется в пространствах с частичной непрерывностью (7)-(8). Основные отличия от векторного метода конечных элементов заключаются в построении дискретной вариационной постановки. Законы сохранения в DG-методе выполняются локально на каждом конечном элементе, за пределами которого базисные функции нулем не доопределяются. Связи между конечными элементами устанавливаются путём введения специальных лифтинг-операторов — скачков и осреднений [28].
Рассмотрим далее построение дискретных вариационных постановок ВМКЭ и различных IP-схем DG-метода.
3. Дискретная вариационная постановка
Дискретизация области моделирования (см. рис. 1) выполняется на тетраэдральном адаптивном разбиении
Th = j П = N П, Vnk, П G Th,k = jnk П П=0 |. (10)
3.1. Конформный векторный метод конечных элементов
Для построения дискретной вариационной постановки конформного ВМКЭ вводятся конечномерное подпространство пространства Ho (rot, П)
H{}(rot, П) С Ho (rot, П) (11)
и базисные функции wk Е H^(rot, П). Дискретная вариационная постановка формулируется следующим образом: найти Eh Е H^(rot, П) + E0 такое, что Vvh Е H^(rot, П) выполняется
У ^-1VxEh ■ VxvhdП + J k2Eh ■ vhdП = -ш J g ■ vhdS. (12)
Q Q Гт
Отметим, что Eh удовлетворяет закону сохранения заряда в слабом смысле [33], так как для Hh(rot, П) выполняется условие вложения (теорема де Рама) [34]. Дискретный аналог напряжённости электрического поля E представим в виде
Eh (x) = £ qiwh(x), (13)
i
где qi — веса разложения Eh по базису пространства H|}(rot, П).
Для тетраэдральных конечных элементов в пространстве H(rot) предлагается достаточно большое количество различных как иерархических, так и интерполяционных векторных базисов [35-37]. Для векторных базисных функций носителями являются рёбра, грани и объём конечного элемента (в зависимости от порядка базиса). Функции делятся на соленоидальные (роторные) и градиентные: первые обеспечивают непрерывность тангенциальной компоненты поля, в то время как вторые отвечают за корректность скачка нормальной компоненты на межэлементных и межфрагментарных границах [33]. Следует отметить, что градиентные функции являются роторно-свободными и дают нулевые вклады в матрицу жёсткости, таким образом, матрица жёсткости для базисов полного порядка становится вырожденной.
В работе [36] выполнен сравнительный анализ векторных базисных функций для решения уравнения Гельмгольца в среде с контрастными электрофизическими свойствами подобластей. Далее в качестве точного решения будет рассматриваться численное решение, полученное на мелкой сетке на базисе первого полного порядка.
3.2. Неконформный разрывный метод Галеркина
Для построения дискретной вариационной постановки DG-метода на разбиении Th вводится пространство [28]
Vh = { v Е Ь2(П)3| v|k Е P1 (K)3VK Е Th} . (14)
Для заданного разбиения Th (10) вводятся обозначения: Fh — внутренняя граница между двумя соседними конечными элементами П+ и П-, n+ и n- — единичные внешние векторы-нормали к грани Fh относительно конечного элемента П+ и П- соответственно. Под соседними понимаются конечные элементы, имеющие общую грань. Тогда скачок есть
[v]T = n+ x v+ + n- x v- (15)
и среднее (осреднение)
{v} = (v+ + v-) /2. (16)
Соотношения (15)—(16) действуют на внутренних границах расчётной области. На внешней границе скачок и среднее определены следующим образом:
[v]T = n х v, (17)
{v} = v. (18)
В (15) и далее индекс T означает, что рассматривается тангенциальная компонента. Скачки и осреднения (15)-(18) вводятся по аналогии с DG-методами для скалярных эллиптических задач.
Дискретная вариационная постановка для DG-метода примет вид: найти Eh Е Vh такое, что для Wh Е Vh выполняется
aDG(Eh, vh) + (k2 Eh, vh)n = (g, vh)dn, (19)
где билинейная форма а^с(-, ■) для IP-постановок определяется как
aDG(u, v) = (j-1Vh х u, Vh х v)n - r J [u]T ■ {/-1Vh х v}ds-
Fh
- j [v]T ■ {/-1 Vh х u}ds + J a/-1 [u]T- [v]Tds, a > 0. (20)
Fh Fh
В зависимости от значения r получают различные варианты IP-постановки [28]:
• r = 0 — неполный IP (IIP);
• r = — 1 — несимметричный IP (NIP, стабилизированная версия метода Баумана — Одена);
• r =1 — классический (симметричный) IP (SIP).
Так как рассматриваются функциональные пространства, обладающие только частичной непрерывностью (тангенциальная непрерывность), то можно предположить, что в данной постановке осреднения на межфрагментарных границах излишни. Поэтому в работе рассматривается также вариационная постановка, построенная на основе IP-схемы, но не содержащая осреднений. Билинейная форма a^G для неё примет вид
aDG(u, v) = (j-1Vh х u, Vh х v)n + J a/-1[u]T-[v]Tds, a > 0. (21)
Fh
Стабилизирующий коэффициент a определяет скорость сходимости итерационного метода решения СЛАУ [27] и может значительно влиять на точность решения.
4. Результаты численного моделирования
Сравнительный анализ вычислительных схем разрывного метода Галеркина и векторного метода конечных элементов выполняется на трёх частотах: 100 кГц, 50 МГц и 1 ГГц.
Таблица 1. Электрофизические характеристики подобластей
Электрофизическая
характеристика Матрица Включения
е, Ф/м 10ео 80ео
а, См/м 10-2 10-6
Гн/м № №
Расчётные области приведены на рис. 1. Габаритные размеры расчётных областей ^матрица = [7см х 7см х 7см]. Диаметр каждого сферического включения 1 см, размеры
параллелепипеидальныхвключений: Пвключение1 = [1 смх5смх1 см], Пвключение2 = [2смх 5 см х 2 см], Пвключение3 = [0.5 см х 5 см х 4 см], Пвключение4 = [5 см х 5 см х 1 см]. В табл. 1 приведены электрофизические характеристики подобластей Пматрица и Пвключения.
Расчёты выполняются на достаточно грубой сетке. Для области с четырьмя сферическими включениями вычисления производятся на сетке из 2 076 тетраэдров. Общее количество степеней свободы для IP-постановок DG-метода N = 24912. Для области с четырьмя параллелепипеидальными включениями сетка состоит из 2 187 тетраэдров (N = 26 244).
В качестве точного решения принимаются результаты, полученные ВМКЭ на мелкой сетке: 171873 тетраэдра, N = 412 334 для области со сферами. Для расчётной области с параллелепипеидальными включениями используется сопоставимая по количеству тетраэдров сетка. Решение выводится по профилю Y = 0.035, Z = 0.02, x Е [0; 0.07], проходящему через центр включений для области со сферическими включениями, и по профилю Y = 0.035, Z = 0.0559, x Е [0; 0.07] для области с параллелепипеидальными включениями.
Оценка относительной погрешности решений, полученных различными модификациями DG-метода, выполняется в векторных нормах
5i
Ez- E'
ext I
E
ext I
M
M
E _Ei(ext)
M M
Ei(ext) z
(22)
S-2
Z
Ez E
ext I
E
ext I
M M
E _Ei(ext)
1/2
1/2
Ei(ext) z
(23)
Ez _ E
ext I
E
ext I
max
1iM
E _Ei(ext)
max
1iM
Ei(ext) z
(24)
где M = 100 — количество точек по профилю; Eext — действительная компонента Ez ("точное" решение); Ez — действительная компонента Ez, найденная численно.
На рис. 2 приведены действительные компоненты Ez для частоты 100 кГц. Габаритные размеры расчётной области значительно меньше длины волны во включении
1
z
1
2
2
2
z
2
z
5
оо
оо
z
оо
—1—'—'—'—1—'—'—'—Е_|_|_|_|_|_|_|_|_|_|_|_I ".
0.02 0.04 0.06 X о.02 0.04 0.06 X
Рис. 2. Действительные компоненты Ez, выведенные по профилю: а — сферические включения, б — параллелепипеидальные включения; сплошная линия — "точное" решение, штрих-пунктир — DG без осреднений, a =1. Частота 100 кГц
и в матрице. Волновой процесс на данной частоте не формируется. В табл. 2 показана погрешность решения каждой из рассмотренных IP-схем.
На графиках (см. рис. 2) введены следующие обозначения: "Точное" решение — действительная компонента поля Ez, принимаемая за "точное" решение, ВМКЭ на мелкой сетке (171873 тетраэдра); DG без осреднений — действительная компонента поля Ez, DG-метод в IP-постановке без осреднений. Если в DG-методе используется стабилизирующий коэффициент, то в легенде указывается его значение, например, DG без осреднений, a = 100.
IIP-постановка и постановка без осреднений дают практически одинаковые результаты (см. табл. 2), позволяют идентифицировать первое включение, но сглаживают второе. Для DG-метода в NIP-постановке введение стабилизирующего коэффициента даёт возможность повысить точность решения в 2 раза. SIP-постановка c a = 1 даёт нефизичное решение, на несколько порядков превышающее точное. Однако подбор стабилизирующего коэффициента, значительно превышающего 1, для SIP-схемы поз-
Таблица 2. Относительная погрешность в векторных нормах. Частота 100 кГц
Норма Параллелепипеидальные включения
NIP, a =1 NIP, a = 100 IP без осреднений IIP, a = 1 SIP, a = 1
S? 3.25E-01 1.36E-01 6.53E-02 5.46E-02 6.74E-01
52 3.28E-01 1.37E-01 6.61E-02 6.10E-02 7.06E-01
51 3.88E-01 1.50E-01 8.65E-02 1.09E-01 7.23E-01
Сферические включения
NIP, NIP, IP без IIP, a = 1 SIP,
a =1 a = 250 осреднений a = 1000
52 5.49E-01 2.98E-01 3.20E-01 2.97E-01 4.96E-01
52 5.54E-01 2.98E-01 3.25E-01 3.00E-01 5.09E-01
51 5.70E-01 3.46E-01 3.40E-01 3.19E-01 5.14E-01
-1-^-^-^-1-"-1---ьГИЗИ^кь. _,_,_,_|_,_,_,_|_,
0.02 0.04 0.06 х 0.02 0.04 0.06 X
Рис. 3. Действительные компоненты Ez, выведенные по профилю: а — сферические включения (a = 100), б — параллелепипеидальные включения (a = 250); сплошная линия — "точное" решение, штрихпунктир — DG без осреднений. Частота 50 МГц
воляет добиться решений, близких к IIP-постановке и постановке без осреднений. Для расчётной области с включениями-параллелепипедами наиболее эффективными оказываются DG-метод в постановке без осреднений и IIP-схема.
На рис. 3 и в табл. 3 приведены результаты моделирования на оптимальных для каждой постановки коэффициентах стабилизации для частоты 50 МГц. На данной частоте габаритные размеры расчётной области также меньше длины волны во включении и в матрице. Волновой процесс не формируется. Однако если на частоте 100 кГц доминировала матрица жёсткости, то при 50 МГц элементы матрицы массы на три порядка больше элементов матрицы жёсткости.
На частоте 50 МГц IIP-схема разрывного метода Галеркина и предложенная в статье постановка без осреднений при a = 1 так же, как SIP- и NIP-постановки, дают ложное решение и требуют подбора a. Для SIP-схемы по сравнению с остальными по-
Таблица 3. Относительная погрешность в векторных нормах. Частота 50 МГц
Норма Параллелепипеидальные включения
NIP, IP без IIP, a = 100 SIP,
a = 100 осреднений a = 250 a = 1500
S? 2.05E-01 2.23E-01 2.23E-01 2.55E-01
52 1.86E-01 2.17E-01 2.24E-01 8.43E-01
52 1.63E-01 2.02E-01 2.09E-01 6.40E-01
Сферические включения
NIP, IP без IIP, a = 500 SIP,
a = 500 осреднений a = 100 a = 1500
52 2.36E-01 2.89E-01 3.04E-01 4.48E-01
52 2.48E-01 2.85E-01 2.97E-01 4.42E-01
52 2.97E-01 2.63E-01 2.71E-01 3.51E-01
а б
0.5 0
-0.5
0.02 0.04 0.06 х 0.02 0.04 0.06 X
Рис. 4. Действительные компоненты Ez, выведенные по профилю: а — сферические включения (а = 1000), б — параллелепипеидальные включения (а = 2500); сплошная линия — "точное" решение, штрихпунктир — DG без осреднений. Частота 1 ГГц
становками необходимо вводить достаточно большой стабилизирующий коэффициент (см. табл. 3). Однако решение в этой постановке оказывается наиболее грубым. Наименьший коэффициент стабилизации требуется для IP-схемы без осреднений.
На рис. 4 приведены действительные компоненты Ez для частоты 1 ГГц. На данной частоте габаритные размеры расчётной области сопоставимы с длиной волны, начинает формироваться волновой процесс. В табл. 4 приведена относительная погрешность для DG-методов при оптимальном выборе стабилизирующего коэффициента а.
Даже при оптимальном выборе коэффициента а на частоте f = 1 ГГц NIP-, IIP-, и SIP-постановки для среды с небольшой объёмной долей включений (см. рис. 1, а) дают нефизичное решения. Для IP-постановки без осреднений и NIP-схемы не удаётся подобрать достаточно эффективные стабилизаторы в области с высоким содержанием включений (см. рис. 1, б).
Таблица 4. Относительная погрешность в векторных нормах. Частота 1 ГГц
Норма Параллелепипеидальные включения
NIP, IP без IIP, SIP,
а = 2500 осреднений а = 2500 а = 2250 а = 2500
sf 8.65E-01 1.15E+00 2.68E-01 3.52E-01
sf 8.37E-01 1.14E+00 2.66E-01 3.50E-01
sf 1.06E+00 1.51E+00 3.66E-01 4.04E-01
Сферические включения
NIP, IP без IIP, SIP,
а = 2000 осреднений а = 1000 а = 1500 а = 2000
sf 1.71E+00 4.09E-01 1.22E+00 1.66E+00
sf 1.61E+00 4.19E-01 1.12E+00 1.52E+00
S* 1.46E+00 4.76E-01 9.73E-01 1.31E+00
a 1200
1000
800
600
400
200
0
100 кГц 50 МГц 1 ГГц ш
Рис. 5. Зависимость стабилизирующего коэффициента a от частоты ш
При решении задачи модификациями DG-метода на приведенных в табл. 1 электрофизических параметрах в рассмотренном диапазоне частот решение СЛАУ итерационными методами затруднительно. Стандартные итерационные методы не обеспечивают сходимости СЛАУ. Двухуровневый решатель [33] не работает, так как для DG-методов дискретный комплекс де Рама не соблюдается: условие непрерывности тангенциальной компоненты поля выполняется глобально, но на уровне конечного элемента введённый неконформный базис не обеспечивает его выполнения. Для DG-методов в настоящей работе СЛАУ решается прямым методом (например методом Гаусса).
Заключение
В работе выполнен сравнительный анализ неконформного DG-метода и конформного векторного метода конечных элементов. Для IP-модификаций DG-метода определён диапазон частот, для которых эти модификации оказываются эффективными.
Изначально IP-схемы и рассмотренные лифтинг-операторы вводились для эллиптических задач. Предлагаемые для уравнения Гельмгольца постановки, построенные по аналогии, без стабилизации работают только для ограниченного диапазона частот, когда доминирует матрица жёсткости. Начиная с частоты 10 МГц преобладает матрица массы и данные постановки требуют введения стабилизирующего коэффициента.
Наиболее стабильно работают IIP-постановка и предложенная в статье постановка без осреднений. Эти схемы могут применяться при широком диапазоне частот для областей с малой и высокой долей включений в расчётной области и требуют наименьшей стабилизации. SIP-модификация работает наименее стабильно и требует создания специальных стабилизаторов. При повышении частоты для всех рассмотренных IP-постановок необходимо также увеличивать стабилизирующий коэффициент (рис. 5).
Для DG-методов, обладающих несимметричной плохообусловленной матрицей СЛАУ, применение стандартных итерационных методов решения СЛАУ нецелесообразно. Двухуровневые решатели также неприменимы. Здесь необходимы либо разработка специальных методов, ориентированных на неконформные схемы в частотном режиме, либо использование прямых методов (например метода Гаусса).
На основе выполненных исследований разработан модифицированный неконформный метод [38], сочетающий в себе основные преимущества ВМКЭ, DG-методов, метода суперэлементов [39] и многомасштабного метода конечных элементов [40].
Список литературы
[1] Емец Ю.П. Дисперсия диэлектрической проницаемости трёх- и четырёхкомпонентных матричных сред // Журнал техн. физики. 2003. Т. 73, № 3. С. 42-53.
[2] Butt H., Wilkinson T.D., Amaratunga G.A.J. FEM modeling of periodic arrays of multi-walled carbon nanotubes // Progress In Electromagnet. Res. 2012. Vol. 22. P. 1-12.
[3] Zhang Y., Cao L., Allegretto W., Lin Y. Multiscale numerical algorithm for 3D Maxwell's equations with memory effects in composite materials // Intern. J. Numer. Anal. Model. Ser. A. 2011. Vol. 1. P. 41-57.
[4] Михайлова Е.И. Моделирование трёхмерного электромагнитного поля в волноводах сложной геометрии // Сборник научных трудов НГТУ. 2012. № 4(70). С. 131-136.
[5] Вендик И.Б., Вендик О.Г. Метаматериалы и их применение в технике сверхвысоких частот (Обзор) // Журнал техн. физики. 2013. Т. 83, № 1. P. 3-28.
[6] Ouchetto O., Zouhdi S., Bossavit A. et al. Modeling of 3-D periodic multiphase composites by homogenization // Microwave Theory and Techniques. IEEE Transactions on. 2006. Vol. 54, No. 6. С. 2615-2619.
[7] Виноградов А.П., Дорофеенко А.В., Зухди С. К вопросу об эффективных параметрах метаматериалов // Успехи физ. наук. 2008. Т. 178, № 5. С. 511-518.
[8] Shurina E.P., Epov M.I., Shtabel N.V., Mikhaylova E.I. The calculation of the effective tensor coefficient of the medium for the objects with microinclusions // Engineering. 2014. Vol. 6, No. 3. P. 101-112. doi: 10.4236/eng.2014.63014.
[9] Kangarlu A., Robitaille P.M.L. Biological effects and health implications in magnetic resonance imaging // Concepts in Magnetic Resonance. 2000. No. 12. P. 321-359.
[10] Nedelec J.C. Mixed finite elements in R3 // Numerische Math. 1980. No. 35. P. 315-341.
[11] Nedelec J.C. A new family of mixed finite elements in R3 // Ibid. 1986. No. 50. P. 57-81.
[12] Schoberl J., Zaglmayr S. High order Nedelec elements with local complete sequence properties // COMPEL: The Intern. J. for Comput. and Math. in Electrical and Electronic Eng. 2005. Vol. 24, No. 2. P. 374-384.
[13] Monk P. Finite Element Methods for Maxwell's Equations. Oxford Univ. Press, 2003. 450 p.
[14] Antonietti P.F., Buffa A., Perugia I. Discontinuous Galerkin approximation of the Laplace eigenproblem // Comput. Methods in Appl. Mechanics and Eng. 2006. No. 195. P. 3483-3503.
[15] Hiptmair R., Moiola A., Perugia I. Error analysis of Trefftz-discontinuous Galerkin methods for the time-harmonic Maxwell equations // Math. of Comput. 2013. No. 82. P. 247-268.
[16] Hesthaven J., Warburton T. Discontinuous Galerkin methods for the time-domain Maxwell sequations // Appl. Comput. Electromagnet. Soc. Newsletter. 2004. No. 19. P. 12-30.
[17] Rapetti F., Bouillault F., Santandrea L. et al. Calculation of eddy currents with edge elements on non-matching grids in moving structures // Magnetics. IEEE Transactions on. 2000. Vol. 36(4). P. 1351-1355.
[18] Dolean V., Gander M.J., Lanteri S. et al. Effective transmission conditions for domain decomposition methods applied to the time-harmonic curl-curl Maxwell's equations // J. of Comput. Phys. (in press).
[19] Эпов М.И., Шурина Э.П., Архипов Д.А. Параллельные конечноэлементные вычислительные схемы в задачах геоэлектрики // Вычисл. технологии. 2013. Т. 18, № 2. С. 95-112.
[20] Cockburn B., Shi K. Conditions for superconvergence of HDG methods for Stokes flow // Math. of Comput. 2013. No. 82(282). P. 651-671.
[21] Hoteit H., Firoozabadi A. Multicomponent fluid flow by discontinuous Galerkin and mixed methods in unfractured and fractured media // Water Resources Res. 2005. No. 41(11). W11412.
[22] Li B.Q. Discontinuous Finite Elements in Fluid Dynamics and Heat Transfer. Springer, 2006. 578 p.
[23] Шокин Ю.И., Шурина Э.П., ИткинА Н.Б. Современные многосеточные методы. Часть I. Многомасштабные методы: Учеб. пособие. Новосибирск: НГТУ, 2009. 64 с.
[24] Canouet N., Fezoui L., Piperno S. Discontinuous Galerkin Time-Domain solution of Maxwell's equations on locally-refined nonconforming Cartesian grids // COMPEL. 2005. No. 24(4). P. 1381-1401.
[25] Lanteri S., Scheid C. Convergence of a discontinuous Galerkin scheme for the mixed time-domain Maxwell's equations in dispersive media // IMA J. Numer. Anal. 2011. RR-7634. P. 1-26.
[26] Houston P., Perugia I., Schotzau D. A review of discontinuous Galerkin methods for Maxwell's equations in frequency-domain // ECCOMAS. 2004. P. 1-16.
[27] Perugia I., Schotzau D. The hp-local discontinuous Galerkin method for low-frequency time-harmonic Maxwell equations // Math. Comput. 2003. Vol. 72. P. 1179-1214.
[28] Houston P., Perugia I., Schotzau D. Mixed discontinuous Galerkin approximation of the Maxwell operator // SIAM J. Numer. Anal. 2004. Vol. 42, No. 1. P. 434-459.
[29] Sarmany D. High-order Finite Element Approximations of the Maxwell Equations: Diss. Univ. of Twente, Netherlands, 2010. 142 p.
[30] Schneebeli A. Interior Penalty Discontinuous Galerkin Methods for Electromagnetic and Acoustic Wave Equations: Diss. Univ. Basel, 2006. 135 p.
[31] Houston P., Perugia I., Schotzau D. An a posteriori error indicator for discontinuous Galerkin discretizations of H (curl)-elliptic partial differential equations // IMA J. of Numer. Anal. 2007. Vol. 27(1). P. 122-150.
[32] Houston P., Perugia I., Schutzau D. Mixed discontinuous Galerkin approximation of the Maxwell operator: Non-stabilized formulation // J. of Sci. Comput. 2005. Vol. 22, No. 1-3. P. 315-346.
[33] Nechaev O.V., Shurina E.P., Botchev M.A. Multilevel iterative solvers for the edge finite element solution of the 3D Maxwell equation // Comput. and Math. with Appl. 2008. No. 55. P. 2346-2362.
[34] Hiptmair R. Finite elements in computational electromagnetism // Acta Numerica. 2002. Vol. 11. P. 237-339.
[35] Webb J.P. Hierarchal vector basis functions of arbitrary order for triangular and tetrahedral finite elements // Antennas and Propagat. IEEE Transactions. 1999. Vol. 47. P. 1244-1253.
[36] Михайлова Е.И., Шурина Э.П. Математическое моделирование высокочастотного электромагнитного поля в волноводных устройствах // Вестник НГУ. Математика. Механика. Информатика. 2013. Т. 13, № 4. C. 102-118.
[37] Yioultsis T.V., Tsiboukis T.D. Development and implementation of second and third order vector finite elements in various 3-D electromagnetic field problems // IEEE Trans. on Magnetics. 1997. Vol. 33, No. 2. P. 1812-1815.
[38] Шурина Э.П., Михайлова Е.И. Моделирование электромагнитных полей в гармоническом режиме неконформными численными методами // Теоретические основы и конструирование численных алгоритмов решения задач математической физики: Тез. докл. XX Всероссийской конф., посвящённой памяти К.И. Бабенко. М.: ИПМ РАН, 2014. С. 115.
[39] Федоренко Р.П. Введение в вычислительную физику. М.: Изд-во МФТИ, 1994. 528 с.
[40] Efendiev Y., Hou T.Y. Multiscale Finite Element Methods: Theory and Applications. Springer, 2009. 234 p.
Поступила в 'редакцию 13 мая 2014 г., с доработки —17 октября 2014 г.
Analysis of numerical schemes for modeling of electromagnetic field in media with contrast inclusions in a wide frequency range
Epov Mikhail I.1, Shurina Ella P.1'2, Mikhaylova Ekaterina I.1'2
Purpose. The purpose of this paper is a comparative analysis of variational formulations of the discontinuous Galerkin method (DG-method) for solution of the Helmholtz equation in a wide frequency range (from 100 kHz to 1 GHz) in a medium with contrast small inclusions of various shapes.
Methodology. The Helmholtz equation for electric field is solved on unstructured tetra-hedral mesh with vector basis functions of the first order in a medium with microinclu-sions. The following variational formulations are considered: the conforming vector Galerkin method, the nonconforming (discontinuous) Galerkin method based on various modifications of the IP approach.
Findings. On the basis of the developed computational schemes for domains with various types of microinclusions, we analyze the wave formation process depending on the elec-trophysical properties of the medium and the excitation frequency of the electromagnetic field. The estimation error for the solution was obtained by various modifications of the discontinuous Galerkin method.
Originality/value. Based on the results of numerical experiments we show that the IP DG-method without averaging and the IIP DG-method (the incomplete interior penalty method) are the most effective and stable approaches in the wide frequency range. We found the irrelevance of the iterative methods for solving the discrete analogues of variational formulations of discontinuous Galerkin method (in the space H(curl, Q)). Keywords: the vector finite element method, discontinuous Galerkin method, composite material, microinclusions.
Received 13 May 2014 Received in revised form 17 October 2014
1Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, 630090, Russia
2 Novosibirsk State Technical University, Novosibirsk, 630073, Russia Corresponding author: Shurina Ella P., e-mail: [email protected]
