CC BY
24
2018, Т. 160, кн. 3 С. 477-494
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)
УДК 519.633.2
О ЧИСЛЕННЫХ МЕТОДАХ РЕШЕНИЯ КВАЗИСТАЦИОНАРНЫХ УРАВНЕНИЙ МАКСВЕЛЛА В НЕОДНОРОДНЫХ СРЕДАХ
А.А. Арбузов1, Р.З. Даутов2, Е.М. Карчевский2, М.М. Карчевский2, Д.В. Чистяков1
1 Компания TGT Oilfield Services, г. Казань, 420108, Россия 2Казанский (Приволжский) федеральный университет, г. Казань, 420008, Россия
Аннотация
Описаны способы сведения квазистационарной системы уравнений Максвелла к системе уравнений относительно либо поля электрической напряженности, либо поля магнитной напряженности. Сформулировано понятие обобщенного решения краевых задач в ограниченных областях для полученных систем уравнений. Устанавливаются условия существования и единственности обобщенных решений. Специально рассмотрен случай осевой симметрии исходных задач, на этой основе предлагается класс тестовых примеров. Эти примеры имеют точные решения, передающие основные особенности решений исходных задач. Сконструированы приближенные методы решения, основанные на конеч-ноэлементных аппроксимациях пространственных операторов. Особое внимание уделено методам на тетраэдральных сетках. Использованы как узловые (лагранжевы) элементы, так и элементы Неделека нулевого и первого порядков. При помощи сконструированных тестовых примеров проведено сравнение вычислительной эффективности предложенных методов конечных элементов. В случае небольших перепадов коэффициентов уравнений и использования регулярных сеток конечных элементов метод, использующий элементы Неделека первого порядка демонстрирует определенные преимущества по точности и трудоемкости.
Ключевые слова: уравнения Максвелла, квазистационарное приближение, метод конечных элементов, тестовые примеры
Введение
Многие прикладные задачи электродинамики, в частности задачи электромагнитной диагностики, приводят к необходимости решения начально-краевых задач для системы уравнений Максвелла в так называемом квазистационарном приближении. Понятно, что, вообще говоря, это может быть выполнено только при помощи численных методов. Для областей сложной формы наиболее эффективными здесь оказываются методы, основанные на конечноэлементной аппроксимации пространственного оператора задачи. Цель настоящей работы - выполнить сравнение двух классов методов конечных элементов, основанных на лагранжевом кусочно-линейном приближении искомого векторного поля и на элементах Неделека. Последние имеют более сложную конструкцию, но позволяют поэлементно учесть необходимое условие на дивергенцию векторного поля. Предварительно мы указываем способы сведения квазистационарной системы уравнений Максвелла к задаче относительно векторного поля электрической напряженности Е либо векторного поля магнитной напряженности Н. Описывается также способ сведения
задачи к скалярному вырождающемуся параболическому уравнению второго порядка в случае, когда исходная задача обладает осевой симметрией. На этой основе строится класс тестовых примеров, передающих основные особенности решений в средах с кусочно-постоянными характеристиками. Далее выполняется сравнение различных вариантов метода конечных элементов. При этом мы моделируем поведение приближенных методов при решении нестационарных задач на основе неявных аппроксимаций по времени.
1. Постановка задачи. Математическая модель
Распространение электромагнитных полей в квазистационарном приближении (см., например, [1-7]) описывается математической моделью, основанной на системе уравнений Максвелла (в гауссовой системе единиц),
4П
rot H =— j (закон Ампера), (1)
c
1 дв
rot E =---— (закон Фарадея), (2)
div D = 4п р (закон Кулона), (3)
div B = 0 (закон Гаусса), (4)
законе Ома
j = a(x) E + js(x,t), (5)
материальных соотношениях
D = e(x) E, (6)
B = f(x) H, (7)
где e, f, a, js - заданные функции - диэлектрическая проницаемость, магнитная проницаемость, проводимость и сторонний ток соответственно, c = const -скорость света.
Считаем, что уравнения определены при 0 < t < ж в ограниченной односвяз-ной области И С R3 с кусочно-гладкой границей дИ.
1.1. Вывод уравнения для определения поля E. Система уравнений (1)-(7) содержит три неизвестных функции: E, H, р. Сокращение числа уравнений обычно достигается за счет процедуры продолжения системы, иначе говоря, за счет повышения порядка уравнений (см., например, [4, 5]). Опишем эту процедуру. Считая, что в начальный момент времени H = 0, из (2) и (7) получаем
t
c
H(t) = -Сrot [ E(r) dr.
f J
0
Подставляя это выражение в (1) и учитывая (5), (7), приходим к уравнению относительно векторного поля E
aE + к rot I fro t J E(r) dr I = -js, (8)
где / = 1//, к = c2/(4n) - постоянная. Полагая теперь, что
t
A(t) = -c j Е(т) dr, (9)
получаем
Ясно, что
Е = - (Ю)
H = /rot A. (11)
A(0)=0. (12)
Используя (10) и (8), приходим к уравнению для векторного поля A
dA
+ к rot (/trot A) = cjs. (13)
Если A - решение уравнения (13), то для E, H, определяемых формулами (10), (11), очевидно, выполняются уравнения (1), (2), (4); функцию р можно определить затем из уравнений (3), (6). В дальнейшем будем полагать, что
div js =0. (14)
Отсюда вследствие (13) получаем, что div aA не зависит от t, и тогда, учитывая (9), будем иметь
div aA = 0. (15)
Сформируем определение обобщенного решения уравнения (13) и соответствующие начальные и граничные условия. При этом, как обычно, мы будем опираться на формулу интегрирования по частям1
/rot* • ф,, = /* • rot$,, + /(. х Ф) • ф(16)
Q Q dQ
которая дает возможность поставить в соответствие уравнению (13) так называемое интегральное тождество
j adA ■ Ф dx + к j /rot A ■ rot Ф dx = c j js ■ Ф dx. (17)
Q Q Q
Это тождество естественно положить в основу определения обобщенного решения (см., например. [8-10]). Будем говорить, что функция A - обобщенное решение второй краевой задачи для уравнений (13), (15) с начальным условием (12), если она удовлетворяет начальному условию (12), интегральному тождеству (17) при любой функции Ф.
Будем говорить, что вектор функция A - решение первой граничной задачи для уравнений (13), (15) с начальным условием (12), если она удовлетворяет граничному условию
v х A = 0, x е dQ, (18)
начальному условию (12) и интегральному тождеству (17) для любой функции Ф, удовлетворяющей тому же граничному условию.
1 Здесь и далее символы «•», « X » обозначают скалярное и векторное произведения трехмерных векторов соответственно.
1.2. Вывод уравнения для определения поля Н. Используя (5), запишем уравнение (1) в виде
c
— <г rot H = E + j (19)
4п
где а = 1/а. Вычислим ротор от обеих частей уравнения (19) и затем используем (2). В результате получим уравнение
—Н
+ к rot(arot H) = rot(ajs). (20)
Если Н - решение уравнения (20), а вектор функция E определена затем при помощи соотношения (19), то для E, Н выполняются уравнения (1), (2). Далее, если положить
Н (0) = 0, (21)
то, как нетрудно убедиться, уравнение (4) есть следствие уравнения (20), а функция р может быть найдена из уравнений (3), (6).
Уравнение (20) по форме совпадает с уравнением (13). Поэтому обобщенные формулировки краевых задач даются для него аналогично предыдущему случаю. Будем говорить, что вектор функция Н - обобщенное решение первой краевой задачи для уравнения (20), если для нее выполнено начальное условие (21), граничное условие
^ х Н = 0, ж е дQ, (22)
а также выполнено интегральное тождество
j l1—^ • Ф dx + к j arot Н • гс^Ф dx = j ajs • rot Ф dx (23)
n n n
при любой функции Ф , удовлетворяющей граничному условию (22)
Будем говорить, что вектор функция Н - обобщенное решение второй краевой задачи для уравнения (20), если для нее выполнено начальное условие (21), а также интегральное тождество (23) при любой функции Ф .
1.3. О корректности краевых задач для квазистационарной системы уравнений Максвелла. Сформулированные выше краевые задачи для уравнений (13), (20) естественно рассматривать в общем виде как краевые задачи для уравнения
—bu
—t + rot (a rot u) = rot f = f, x е Q, 0 <t<T, (24)
с дополнительным условием
div(bu) =0, x е Q, 0 <t <T. (25)
При соответствующем выборе функций a, b мы получаем либо уравнение (13), либо уравнение (20), так что излагаемые ниже результаты будут относится к постановке уравнений Максвелла в терминах как поля E, так и поля Н. Исследование существования и единственности обобщенных решений этих задач в соответствующих функциональных пространствах может выполнено методами [10-13].
2. Осесимметричные задачи для квазистационарной системы
уравнений Максвелла
В этом разделе описывается методика решения начально-краевой задачи —bu
+rot(arotu) = f, x е Q, 0 <t<T, (26)
div(bu) =0, ж е О, 0 <t<T, (27)
u х v = 0, ж е д О, 0 <t<T, (28)
u(x, 0) = 0, ж е О, (29)
при специальных предположениях об исходных данных, а именно предполагается, что область О - круговой цилиндр высоты 2L и радиуса R. Ось жз декартовой системы координат xi, Ж2, жз направлена по оси цилиндра. Начало координат лежит в центре цилиндра. Введем в рассмотрение также цилиндрическую систему координат r, ф, z, связанную с декартовой системой координат обычным образом:
ж1 = r cos ф, ж2 = r sin ф, жз = z. (30)
Будем считать в дальнейшем функции a, b кусочно постоянными, поверхности разрыва этих функций могут быть расположены только на координатных поверхностях вида r = const или z = const. Будем полагать также, что правая часть уравнения (24) в представлении декартовыми координатами имеет вид
f (x,t) = (-x2g(r, жз, t), xig(r, жз, t), 0), r = (ж1 + ж2)1/2. (31)
Функция g предполагается кусочно непрерывной и кусочно непрерывно дифференцируемой по ж, причем поверхности разрыва функции g лежат на координатных поверхностях вида r = const или z = const. Ясно, что div f = 0 в областях гладкости функции g. Нетрудно видеть, что в каждой точке ж е О вектор f (ж, t) лежит в касательной плоскости к указанным координатным поверхностям. Поэтому элементарные вычисления, основанные на формуле Грина
Jf ■Vndx = -J п div Fdx + J nF ■ vdx, (32)
w w y
где F, п дифференцируемы в области ш, y - граница области ш, v - единичная внешняя нормаль к y , приводят к тождеству
J f(x,t) ■Уф(ж) dx = 0 V ф е C°° (О). (33)
п
Пример. Для задачи о возбуждении электромагнитного поля генерирующей катушкой с обмоткой, занимающей полый цилиндр
Оо = {ro < r < ri, 0 < ф < 2п, zo < z < zi}
и описываемой уравнением (13) в терминах потенциала A, правая часть уравнения (26) определяется заданием функции g при помощи соотношений
Г1 /r, (r,z) е Оо, 0 <t<T,
g(r, z, t) =
[0, (r,z) / О0, 0 <t<T. Отметим, что при этом \f \ = 1, если (r,z) е Оо, 0 <t <T. Решение задачи (26)-(29) будем разыскивать в виде
u(x,t) = (—x2v(r,xз,t),Xlv(r,xз,t), 0), (34)
где v(r,xз,t) - подлежащая определению функция. Нетрудно проверить, что соотношение (34) определяет осесимметричное векторное поле, ортогональное оси жз .
Функция rv есть физическая ф-компонента этого поля. Предполагается, что при любом t € (0, T) функция v как функция декартовых координат принадлежит пространству H¿(П). Аналогично предыдущему, используя формулу (16), нетрудно убедиться, что тогда div(bu) =0, х € П, при любом t € (0, T). Как следствие получаем, что вне поверхностей разрыва функции b выполнено равенство div u = = 0, поэтому, используя хорошо известную формулу rot(rot u) = — Ди + V(div u), получаем, что в областях, где функции a, b сохраняют постоянные значения, для вектор-функции u выполняется уравнение
ди
b— = аДи + f. (35)
dt v '
Отметим прежде всего, что для любой вектор функции и = (u\,u2,u^), представленной декартовыми компонентами, выполнено равенство
Дu = (Дul, Дu2, Дuз). (36)
В правой части (36) фигурирует скалярный оператор Лапласа. Подставляя (31) и (34) в (35), получим уравнения
dv
xkb^t = аДхиv + xkg, k = 1, 2. (37)
Отсюда вследствие (30) будем иметь
,dv idv
cos ф b — = аД cos фv + cos ф g, sin ф b — = аД sin фv + sin ф g, (38)
dt dt
где v = rv, V = rg.
Записывая скалярный оператор Лапласа в цилиндрических координатах:
д 1 d f du\ +1 d2u + d2u (39)
r dr \ dr J r2 дф2 dz2 '
и выполняя элементарные выкладки, каждое из уравнений (38) приведем к виду
dv í 1 d í dv\ d2v V\ v (4)
dt \r dr \ dr) dz2 r2) g
Отметим, что дифференциальный оператор в правой части (40) есть ф-компонента векторного оператора Лапласа в цилиндрических координатах в осесимметричном случае. Уравнению (40) можно придать более компактную форму. Вводя новую неизвестную функцию w = rV и принимая во внимание тождественные преобразования
d í d / w\\ w di dw w\ w d2w 1 dw d /1 (41)
dr V dr \r J I r2 dr V dr r ) r2 dr2 r dr dr \r dr I '
получим, что
dw d 1 dw d2w
bdt = rdr(ardr + adZ2 +p, (42)
где р = т2д. Решая уравнение (42) с однородными граничными условиями Дирихле на границе области {0 < т < К, — Ь < г < Ь}
w(0,z,t) = •ш(Я,г,г)=0, — Ь < г < Ь, 0 < К ж, (43)
т(г, -Ь,г) = т(т,Ь,г)=0, 0 < г < К, 0 < г< ж, с начальным условием
т(г,г, 0) = 0, 0 < г < К, -Ь < г < Ь,
(44)
(45)
и выполняя обычные условия сопряжения (непрерывности функции т и «потоков»),
0,
1 дт г дг
0,
дт дг
0,
(46)
где $ - поверхность разрыва коэффициента а, [ • обозначает скачок при переходе через поверхность $, и используя затем формулу (34), получим вектор-функцию и, удовлетворяющую всем условиям задачи (26)-(29). В силу известной теоремы единственности (см. п. 1.3) других решений эта задача иметь не может.
5
5
3. Тестовый пример
Решается задача
ди
— + гсЛ(аrotи) = х е О, 0 <г <Т,
и = 0, х е о, 0 <г < ж, и(х,г) = 0, х е до, 0 <г < ж, 0,
(47)
(48)
(49)
(50)
и(х, 0) = 0, х е О.
Область О - прямой круговой цилиндр, см. разд. 2. Функция а = а(г) предполагается положительной, кусочно-постоянной и определяется следующими соотношениями:
{аь 0 < г < К1,
а2, К1 < г < К2, (51)
а3, К2 < г < К.
Здесь К\, К2 , 0 < К\ < К2 < К - заданные постоянные. Для описания функций ], и введем в рассмотрение функцию X = X(г), определяемую соотношениями
схг
2
0 < г < К1,
X(г) = { -г3/(3а2) + г2С2 + сз, К1 < г < К22,
(52)
С4(К2 - г2),
К2 < г < К,
где вектор постоянных с = (с1, С2, сз, С4) - решение системы линейных алгебраических уравнений Ас = Ь,
К21
А
-К12
1
0
2а 1К1 -2а2 К1 0
\
0
К22 2а2К2
1 К22 - К2 0 2а3К2
(-К\/(3а2 Д -К21 К3/(3а2)
К22
\
/
Нетрудно убедиться, что функция X = X(г) - решение задачи
(53)
г(аг X')' + х(г) = 0, г е (0,К), г = К1,г = К2,
(54)
0
0
Ь
I r, r е [Д1,Д2], (55)
*RiR](r) = [o, rе [Ri,R2], (55)
[X]\r=R1 =[X]\r=R2 =0, [ar-1X']\r=Ri = [ar-1X']|r=R2 = 0, (56)
X (0)= X (R)=0. (57)
Пусть, далее, Z(z), z е [—L, L], - дважды непрерывно дифференцируемая функция, ©(t), t е [0, то] - непрерывно дифференцируемая функция, причем ©(0) = 0, Z(—L) = Z(L) = 0. Тогда
u(x,t) = ©(t)Z(z) (—x2X(r)/r2, x1X(r)/r2, 0) (58)
есть решение задачи (47)-(50) при
f (x,t) = (— x2g(r,z,t),xig(r,z,t), 0), (59)
где
g(r,z,t)= p/r2, (60)
p = (©'(t)Z(z)X(r) + ©(t)Z(z)x[Ri,R2](r) — a(r)©(t)Z"(z)X(r)) . (61)
Можно положить, например, ©(t) = 1 — e-at, а > 0, Z(z) = 1 — z2/L2 . Элементарными вычислениями проверяется, что функция
w(r,z,t) = X (r)Z (z)©(t)
и функция p, определяемая соотношением (61), удовлетворяют всем условиям задачи (42)-(46) при b = 1, поэтому они могут быть использованы при тестировании приближенных методов ее решения. Рис. 1-3 иллюстрируют результаты приближенного решения тестового примера (параметры определены в разд. 7, в котором описываются численные эксперименты).
4. Метод конечных элементов
Опишем метод конечных элементов (МКЭ) (см., например, [14]) для приближенного решения квазистационарной задачи. Исходная задача в терминах поля E имеет следующий вид:
дА
а — +rot(xrotА) = cjs, x е Q, t е (0,Т], (62)
A\t=0 =0, А х = ° (63)
Здесь
2
с" 1 дА ^ X =4^' -~сдъ = К (64)
Будем исходить из обобщенной формулировки этой задачи: найти вектор функцию А, удовлетворяющую условиям (63) и тождеству
j адА • А' dx + j xrot А • rot А' dx = j cjs • А'
dx, (65)
для любого А х V = 0 при всех Ь € (0, Т].
Метод конечных элементов полностью определяется конструкцией конечномерного подпространства основного пространства решаемой задачи. Обозначим это подпространство через V^ .
Рис. 1. Первая компонента ух приближенного решения тестового примера при использовании линейных узловых элементов и сетки с N = 6358 узлами на плоскости х3 = 0 (слева) и погрешность \их — ух\ ее вычисления (справа)
|и1-у1| Гог х3=0
0.05 0
-0.05
1
-1 1
-0.5 0 х.
Рис. 2. Первая компонента ух приближенного решения тестового примера при использовании элементов Неделека нулевого порядка и сетки с N = 6358 узлами на плоскости х3 = 0 (слева) и погрешность \их — ух\ ее вычисления (справа)
*3=°
2
Рис. 3. Первая компонента ух приближенного решения тестового примера при использовании элементов Неделека первого порядка и сетки с N = 6358 узлами на плоскости х3 = 0 (слева) и погрешность \их — ух\ ее вычисления (справа)
Приближенное решение задачи (65) по методу конечных элементов определяется следующим образом: найти вектор функцию Ah G Vh, Ah\t=o = 0, такую,
что QA
j a A • A' dx + j xr°t Ah • rot A' dx = j cjs • A' dx (66)
Qh Qh Qh
для любого A' G Vh при всех t G (0, T]. Здесь Qh - область определения вектор функций из Vh - аппроксимация области Q.
Пространство Vh определяется обычно как кусочно-полиномиальное пространство: область определения решения Q С R3 разбивается на совокупность элементарных подообластей Ki, i = 1, 2,...,Nk , а элементы из Vh определяются как
Nk
полиномы на каждом Ki. Полагается, что Qh = U Ki. При этом подобласти Ki,
i=i
называемые также конечными элементами, а также полиномы на них определяются так, чтобы обеспечить вложение Vh С V в случае, когда Q = Qh. В нашем случае это равносильно требованию непрерывности касательных компонент вектор-функций из Vh на межэлементных границах.
Конечные элементы Ki, как правило, являются многогранниками в R3 с плоскими гранями. Наиболее распространены на практике тетраэдры (четырехгранники с четырьмя треугольными гранями), призмы (пятигранники с двумя треугольными и тремя четырехугольными гранями) и гексаэдры (шестигранники с шестью четырехугольными гранями). Наиболее универсальными являются тетраэдры, поскольку области достаточно сложной формы легко разбиваются на такие элементы. Ясно, что Qh = Q только в том случае, когда Q сама является многогранником с плоскими гранями.
Разбиение области на конечные элементы, как и множество
Th = {Ki,K2,...,Knk },
принято называть триангуляцией области. В МКЭ обычно предполагается, что любые два элемента Th либо могут не иметь общих точек, либо могут иметь общую грань, или общее ребро, или общую вершину. Грани элементов будем обозначать буквой F (face), ребра буквой E (edge), вершины буквой V (vertex). Произвольный конечный элемент будем обозначать буквой K .
Таким образом, допустимое для задачи (65) пространство Vh, называемое пространством конечных элементов, определяется следующим образом:
Vh = {vh G Ho(rot;Qh) : vh\x G P(K) VK G Th}. (67)
Здесь H0(rot;Q) = {u G H(rot;Q) : v x u = 0, x G dQ}, H(rot; Q) = {u G [L2(Q)]3 : rot u G [L2(Q)]3}, P(K) = P(K) xP(K) x P(K), P(K) есть некоторое пространство полиномов, определенных на элементе K.
5. Система уравнений МКЭ
Опишем алгоритм построения системы уравнений для определения приближенного решения задачи, то есть вектор-функции Ah. Пусть вектор-функции образуют базис в Vh . Это означает, что произвольная функция из Vh может быть однозначно представлена в виде линейной комбинации этих функций. Тогда
Nd
Ah(x,t) = ^ aj (t)Vj (x), x G Qh, t G [0,T]. (68)
j=i
Поскольку тождество (66) справедливо для всех А' € Ун, то оно должно быть справедливо и при А = фг для всех г = I,..., Ыл. Подставляя в тождество (66) разложение (68) и А = фг, придем к системе обыкновенных дифференциальных уравнений (ОДУ) для определения неизвестных а^ (г), 3 = 1,..., Ыл,
^ . Nd
( / а ф: • фг ¿х^ а^(г) + ( / хго^з • rotфi ¿х^ а^(г) = • фг ¿х, (69)
где г = 1, 2, ... , , а: (г) = ¿а^ (г)/¿г.
Запишем (69) в векторно-матричном виде. Для этого определим вектор-столбцы неизвестных, их производных, а также вектор правой части:
а(г) = (а1(г),а2(г),.. .,ащ(г))т, а'(г) = (а1(г),а2 (г),...,аЩ (г))т,
I = (У 3 • фг X .
Пн
Определим также матрицы М и А размера М х М:
(Г . N (Г . N
М = = J а ф^ • фг ¿х ^ г ^ , А = |Аг:- = J xrotфj • rotфг ^ ^ .
^н ^н
Тогда (69) примет вид автономной системы ОДУ
Ма'(г) = ^(а(г)), г е (0,т], (70)
где ^(а) = I — Аа. В силу начальных условий Ан|4=0 = 0 получаем
а(0) = 0. (71)
Решив задачу Коши (70), (71), применяя тот или иной численный метод, приближенное решение исходной задачи найдем по формуле (68). Например, неявный метод Эйлера (неявная схема первого порядка аппроксимации) для решения задачи Коши (70), (71) имеет вид
(М + т: А)аР+1 = рт, з =0, 1,...,т — 1,
где = Ма: + т: (I — Аа:). Неявная схема является является безусловно устойчивой, и т:, шаги сетки по времени, могут выбираться лишь из соображений точности.
6. Пространства тетраэдральных конечных элементов
Опишем используемые далее пространства конечных элементов, которые основаны на разбиении области на тетраэдры.
6.1. Пространство линейных узловых элементов. Пространство линейных узловых элементов Р1 определяется следующим образом:
Ун = {ун е С(Пн) : Ук\к = ак + Акх VК еТн, х у\дПн = 0},
где ак € R3, Ак € R3x3, x = (x,y,z). Таким образом, компоненты вектор-функций из Vh являются полиномами первого порядка на каждом конечном элементе. Поскольку касательные компоненты функций из Vh непрерывны на межэлементных границах, то Vh С H(rot; Qh).
Определим базис в Vh. Для этого вершины тетраэдров назовем узлами сетки и пронумеруем их от 1 до N в каком-либо порядке. Координаты узла с номером i обозначим через Xi. С каждым узлом сетки Xi свяжем кусочно-линейную непрерывную в Qh скалярную функцию ^ = ф^х) такую, что фi(xi) = 1 и фi(xj) = 0 при i = j. Тогда Nd = 3N базисных вектор-функций vj (х), j = 1, 2,..., 3N, определяются следующим образом:
Vi = (фг, 0,0)т, vn+i = (0, фг, 0)T, V2N+i = (0, 0, фг)т,
где i = 1, 2,..., N. Разложение произвольной функции Ah из Vh по этому базису, имеющее вид
Nd
Ah (х) = Y1 aj Vj ^
j = l
можно записать также в виде
Ah(x)
(A i,h(x)\ A2,h(x) \A3,h(x)J
Nd
У/aiфi(x) i=i Nd
+iфi (x)
i=i Nd
^a2N+iфi(x)
Из определения функций ф^ следует, что коэффициенты разложения aj имеют следующий смысл: для г = 1, 2, . . . , N получаем а^ = А\, аN= А2,н(хг),
а2И+г = Аз , ъ,(Хг).
6.2. Пространства элементов Неделека. В качестве второго примера рассмотрим пространство конечных элементов Неделека (см. [15, 16]) нулевого порядка N0, которое определяется следующим образом:
Vh = {vh € H(rot; Qh) : VhK = aK + ск x x VK € %, Vh x v|
snh
0},
где а к, ск € К3. Таким образом, компоненты вектор-функций из Vh являются неполными линейными полиномами на каждом конечном элементе. Вектор-функции из Vh, вообще говоря, не являются непрерывными во всей расчетной области, а терпят разрыв на межэлементных границах.
Базисные функции этого пространства связаны не с вершинами, а с ребрами тетраэдров. Ребра всех тетраэдров пронумеруем от 1 до NE в каком-либо порядке. Каждому ребру ej ставится в соответствие базисная функция так, что в разложении произвольной функции Ah из Vh по этому базису, имеющему вид
NE
^(х) = ^ ^ (x), j=1
коэффициенты разложения aj имеют следующий смысл:
У Ah • (xk - xi) ds,
a
j
где (х^, xi) - граничные точки ребра ej . Такой выбор коэффициентов обеспечивает включение Vh С H(rot;Qh).
Пусть K - некоторый тетраэдр, e = {i,j} - его ребро с вершинами xi и Xj. Тогда сужение базисной функции фе на элемент K и его ротор имеют вид
= Ai VAj - Xj VAi, rot = 2VAi xVAj. (72)
Здесь Ak, 1 < fc < 4 - барицентрические координаты элемента K.
В качестве третьего примера укажем пространство конечных элементов Неделе-ка первого порядка N1, [15, 16]. Пространство полиномов P(K) на элементе K (см. (67)) определяется как множество полиномов первого порядка. Размерность этого пространства равна 2Ne и с каждым ребром связываются две базисные функции, а именно, наряду с функцией (72) с ним связывается также функция
фе = AiVAj + Xj VAi, rot фе = 0.
7. Анализ точности МКЭ для уравнений с разрывными коэффициентами
Приведем результаты решения стационарной задачи для уравнения Максвелла описанными выше тремя методами. Эта задача иллюстрирует поведение приближенных методов при решении нестационарной задачи на основе неявных аппроксимаций по времени. Мы ограничимся рассмотрением регулярных сеток конечных элементов и небольших перепадов коэффициентов уравнения. Сравнению схем при сильной нерегулярности сеток и больших перепадах коэффициентов уравнения при решении нестационарной задачи предполагается посвятить отдельную работу.
7.1. Тестовый пример. Изменим вектор функции u и f, заданные формулами (58) и (59), так, чтобы они не зависели от времени, то есть положим ©(t) = 1. Тогда эти функции удовлетворяют условиям
bu + rot(a rot u) = g, x G Q, (73)
u = 0, x G дQ, (74)
где g = f + bu, область Q - цилиндр радиуса R, высоты 2L,
' a1, 0 < r < R1, a(r) = <( a2, Ri < r < R2, (75)
^3, R2 < r < R.
Пусть b = 1, ai = аз = 1, a2 = 0.5, L =1, R =1, Ri = 0.2, R2 = 0.5. Задача (73)-(75) решалась в предположении, что вектор-функция u в левой части (73) неизвестна. Численное решение обозначаем в дальнейшем через y.
7.2. Результаты вычислений. Расчеты, основанные на применении линейных узловых элементов и элементов Неделека, проводились для описанного в п. 7.1 тестового примера с разрывными коэффициентами последовательно на четырех сетках. Основные характеристики конечно-элементных сеток и объем требуемой памяти для элементов P1, N0 и N1 приведены в табл. 1, 2. При использовании элементов N0 возникают матрицы M и A существенно больших порядков, чем при применении элементов P1, но места в памяти компьютера для хранения матрицы A + M требуется меньше.
Табл. 1
Характеристики конечно-элементных сеток для элементов Р1, N0 и N1
Число узлов, N 291 485 918 6358
Число степеней свободы для Р1 Число степеней свободы для N0 Число степеней свободы для N1 873 1534 3068 1455 2800 5600 2754 5491 10982 19074 41129 82258
Табл. 2
Требуемая память для сетки с 6358 узлами, элементы Р1, N0 и N1
Тип элементов Р1 N0 N1
Число степеней свободы, Nd 19074 41129 82258
Размер матрицы М (в КЬ) 298 7759 30553
Размер матрицы А (в КЬ) 9028 6733 6894
Размер матрицы А + М (в КЬ) 9028 7759 30553
Табл. 3
Скорость сходимости
Число узлов, N 291 485 918 6358 Скорость сходимости
Р1, \и - у\\ь2 0.0209 0.0105 0.0077 0.0036 Ъ1АА
N0, \\и - у\\ь2 0.0266 0.0227 0.0117 0.0056 к1-ьа
N1, \\и - у\\ь2 0.0426 0.0424 0.0038 0.0010 к2-0'
Табл. 4
Сетка с 6358 узлами, элементы Р1, N0; сетка с 918 узлами, элементы N1
Тип элементов Р1 N0 N1
Размер матрицы М (в КЬ) 298 7759 3920
Размер матрицы А (в КЬ) 9028 6733 888
Размер матрицы А + М (в КЬ) 9028 7759 3920
В табл. 3 сравнивается скорость сходимости численных схем, основанных на применении элементов Р1, N0 и N1, при решении тестового примера. Скорость сходимости для элементов Р1 несколько ниже, чем для элементов N0 и существенно ниже, чем для N1. Примерно одинаковый уровень точности достигается для элементов Р1 и N0 для сетки с 6358 узлами, а для элементов N1 - всего с 918. В табл. 4 сравниваются требуемые вычислительные ресурсы. Видно, что при использовании элементов N1 их требуется существенно меньше.
Приведенные результаты расчетов позволяют сделать вывод о том, что в случае небольших перепадов коэффициентов уравнения и использования регулярных сеток конечных элементов все три метода работоспособны. Для достижения точности одного порядка при использовании элементов N1 требуется меньше вычислительных ресурсов, чем при использовании элементов Р1 и N0.
Благодарности. Работа выполнена за счет средств субсидии, выделенной Казанскому федеральному университету для выполнения государственного задания в сфере научной деятельности, проект № 1.12878.2018/12.1.
Литература
1. Ландау Л.Д., Лифшиц Е.М. Электродинамика сплошных сред. - М.: Наука, 1982. -621 с.
2. Дюво Г., Лионе Ж.-Л. Неравенства в механике и физике - М.: Наука, 1980. - 382 с.
3. Тамм И.Е. Основы теории электричества. - М.: Наука, 1976. - 616 с.
4. Галанин М.П., Попов Ю.П. Квазистационарные электромагнитные поля в неоднородных средах. - М.: Наука, Физматлит, 1995. - 315 с.
5. Acevedo R., Meddahi S., Rodriguez R. An E-based mixed formulation for a time-dependent eddy current problem // Math. Comput. - 2009. - V. 78, No 268. - P. 1929-1949.
6. Калинин А.В., Морозов С.Ф. Система уравнений Максвелла в квазистационарном магнитном приближении // Вестн. ННГУ. Сер. Матем. моделирование и оптимальное управление. - 2001. - Вып. 1. - С. 97-106.
7. Rapettia R., Rousseaux B. On quasi-static models hidden in Maxwell's equations // Appl. Numer. Math. - 2014. - V. 79. - P. 92-106. - doi: 10.1016/j.apnum.2012.11.007.
8. Ладыженская О.А. Краевые задачи математической физики. - М.: Наука, 1973. -407 с.
9. Михлин С.Г. Линейные уравнения в частных производных. - М.: Высш. шк., 1977. -423 с.
10. Карчевекий М.М., Павлова М.Ф. Уравнения математической физики. Дополнительные главы. - СПб.: Лань, 2016. - 276 с.
11. Калинин А.В. Математические задачи физической диагностики. Корректность задач электромагнитной теории в стационарном и квазистационарном приближении. -Н. Новгород: Изд-во Нижегор. ун-та, 2007. - 121 с.
12. Гаевекий X., Грегер К., Захариае К. Нелинейные операторные уравнения и операторные дифференциальные уравнения. - М.: Мир, 1978. - 336 с.
13. Zeidler E. Nonlinear Functional Analysis and its Applications. II/A: Linear monotone operators. - N. Y.: Springer-Verlag, 1990. - XVIII, 467 p. - doi: 10.1007/978-1-46120985-0.
14. Сьярле Ф. Метод конечных элементов для эллиптических задач. - М.: Мир, 1980. -512 с.
15. Nedelec J.C. A new family of mixed finite elements in R3 // Numer. Math.- 1986. -V. 50, No 1. - P. 57-81. - doi: 10.1007/BF01389668.
16. Rodriguez A.A., Valli A. Eddy Current Approximation of Maxwell Equations. Theory, Algorithms and Applications. - Springer-Verlag, 2010. - XIII, 347 p. - doi: 10.1007/97888-470-1506-7.
Поступила в редакцию 15.03.18
Арбузов Андрей Александрович, кандидат физико-математических наук Компания TGT Oilfield Services
ул. Магистральная, д. 59, г. Казань, 420108, Россия E-mail: info@tgtoil.com
Даутов Рафаил Замилович, доктор физико-математических наук, профессор кафедры вычислительной математики
Казанский (Приволжский) федеральный университет
ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: Rafail.Dautov@kpfu.ru
Карчевский Евгений Михайлович, доктор физико-математических наук, профессор кафедры прикладной математики
Казанский (Приволжский) федеральный университет
ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: Evgenii.Karchevskii@kpfu.ru
Карчевский Михаил Миронович, доктор физико-математических наук, профессор кафедры вычислительной математики
Казанский (Приволжский) федеральный университет
ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: Mikhail.Karchevsky@kpfu.ru
Чистяков Дмитрий Владимирович, кандидат физико-математических наук Компания TGT Oilfield Services
ул. Магистральная, д. 59, г. Казань, 420108, Россия E-mail: info@tgtoil.com
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online) UCHENYE ZAPISKI KAZANSKOGO UNIVERSITETA. SERIYA FIZIKO-MATEMATICHESKIE NAUKI (Proceedings of Kazan University. Physics and Mathematics Series)
2018, vol. 160, no. 3, pp. 477-494
On Numerical Methods for Time-Dependent Eddy Current Problems for the Maxwell Equations in Inhomogeneous Media
A.A. Arbuzova* , R.Z. Dautovb**, E.M. Karchevskiib***, M.M. Karchevskiibm, D.V. Chistyakova*
aTGT Oilfield Services Company, Kazan, 420108 Russia bKazan Federal University, Kazan, 420008 Russia E-mail: *info@tgtoil.com, **Rafail.Dautov@kpfu.ru, ***Evgenii.Karchevskii@kpfu.ru, **** Mikhail.Karchevsky@kpfu.ru
Received March 15, 2018 Abstract
E-based and H-based formulations for time-dependent eddy current problems for the Maxwell equations in inhomogeneous media have been derived. The concept of generalized solutions for boundary-value problems in bounded regions for obtained systems of equations has been formulated. Conditions for the existence and uniqueness of the generalized solutions have been established. Axisymmetric problems have been thoroughly considered, and a class of test problems has been proposed. Their exact solutions have the same key features as the solutions of the original problems. Numerical methods based on finite element approximations of the three-dimensional operators have been constructed. Particular attention has been paid to the methods on tetrahedral elements. Linear Lagrange elements, zero-order and first-order Nedelec elements have been used. The computational efficiency of the proposed finite element approximations has been analyzed by solving the constructed test problems. For small gaps in the coefficients of equations and regular finite element meshes, the first-order Nedelec elements have certain advantages in terms of accuracy and computational costs.
Keywords: Maxwell equations, eddy current approximation, finite element method, test problems
Acknowledgments. The research was funded by the subsidy allocated to Kazan Federal University for the state assignment in the sphere of scientific activities, project no. 1.12878.2018/12.1.
Figure Captions
Fig. 1. The first component y1 of the approximate solution of the test problem (the left panel) and the error \u1 — y1\ in its calculation (the right panel) on the plane x3 = 0. Here, the linear Lagrange elements and the mesh with N = 6358 nodes were used.
Fig. 2. The first component y1 of the approximate solution of the test problem (the left panel) and the error \u1 — y1\ in its calculation (the right panel) on the plane x3 = 0. Here, the zero-order Nedelec elements and the mesh with N = 6358 nodes were used.
Fig. 3. The first component y1 of the approximate solution of the test problem (the left panel) and the error \u1 — y1\ in its calculation (the right panel) on the plane x3 = 0. Here, the first-order Nedelec elements and the mesh with N = 6358 nodes were used.
References
1. Landau L.D., Pitaevskii L.P., Lifshitz E.M. Electrodynamics of Continuous Media. Butterworth-Heinemann, 1984. 460 p.
2. Duvant G., Lions J.L. Inequalities in Mechanics and Physics. Berlin, Heidelberg, SpringerVerlag, 1976. xvi, 400 p. doi: 10.1007/978-3-642-66165-5.
3. Tamm I.E. Osnovy teroii elektrichestva [Fundamentals of The Theory of Electricity]. Moscow, Mir, 1979. 684 p. (In Russian)
4. Galanin M.P., Popov Yu.P. Kvazistatsionarnye elektromagnitnye polya v n,eod,n,orod,n,ykh sredakh [Quasi-Stationary Electromagnetic Fields in Inhomogeneous Media]. Moscow, Nauka, Fizmatlit, 1995. 315 p. (In Russian)
5. Acevedo R., Meddahi S., Rodriguez R. An E-based mixed formulation for a time-dependent eddy current problem. Math. Comput., 2009, vol. 78, no. 268, pp. 1929-1949. doi: 10.1090/S0025-5718-09-02254-6.
6. Kalinin A.V., Morozov S.F. The Maxwell equations system in the quasi-stationary magnetic approximation. Vestn. Nizhegorod. Gos. Univ. Ser. Mat. Model. Optim. Upr., 2001, no. 1, pp. 97-106. (In Russian)
7. Rapettia R., Rousseaux B. On quasi-static models hidden in Maxwell's equations. Appl. Numer. Math., 2014, vol. 79, pp. 92-106. doi: 10.1016/j.apnum.2012.11.007.
8. Ladyzhenskaya O.A. The Boundary Value Problems of Mathematical Physics. New York, Springer-Verlag, 1985. xxx, 322 p. doi: 10.1007/978-1-4757-4317-3.
9. Mikhlin S.G. Lineinye uravneniya v chastnykh proizvodnykh [Linear Partial Differential Equations]. Moscow, Vyssh. Shk., 1977. 423 p. (In Russian)
10. Karchevskii M.M., Pavlova M.F. Uravneniya matematicheskoi fiziki. Dopolnitel'nye glavy [Equations of Mathematical Physics. Additional Chapters]. St. Petersburg, Lan', 2016. 276 p. (In Russian)
11. Kalinin A.V. Matematicheskie zadachi fizicheskoi diagnostiki. Korrektnost' zadach elektro-magnitnoi teorii v statsionarnom i kvazistatsionarnom priblizhenii [Mathematical Problems of Physical Diagnostics. The Correctness of the Problems of the Electromagnetic Theory in the Stationary and Quasi-Stationary Approximations]. Nizhny Novgorod, Izd. Nizhegorod. Univ., 2007. 121 p. (In Russian)
12. Gaevsky H., Greger K., Zakharias K. Netinei,nye operatornye uravneniya i operatornye dif-ferentsial'nye uravneniya [Nonlinear Operator Equations and Operator Differential Equations]. Moscow, Mir, 1978. 336 p. (In Russian)
13. Zeidler E. Nonlinear Functional Analysis and Its Applications. II/A: Linear Monotone Operators. New York, Springer-Verlag, 1990. xviii, 467 p. doi: 10.1007/978-1-4612-0985-0.
14. Ciarlet P.G. The Finite Element Method for Elliptic Problems. SIAM, 2002. xxiii, 529 p. doi: 10.1137/1.9780898719208.
15. Nedelec J.C. A new family of mixed finite elements in R3. Numer. Math., 1986, vol. 50, no. 1, pp. 57-81. doi: 10.1007/BF01389668.
16. Rodriguez A.A., Valli A. Eddy Current Approximation of Maxwell Equations. Theory, Algorithms and Applications. Milan, Springer-Verlag, 2010. xiii, 347 p. doi: 10.1007/97888-470-1506-7.
Для цитирования: Арбузов А.А., Даутов Р.З., Карчевский Е.М., Карчев-/ ский М.М., Чистяков Д.В. О численных методах решения квазистационарных урав-\ нений Максвелла в неоднородных средах // Учен. зап. Казан. ун-та. Сер. Физ.-матем.
науки. - 2018. - Т. 160, кн. 3. - С. 477-494.
For citation: Arbuzov A.A., Dautov R.Z., Karchevskii E.M., Karchevskii M.M., / Chistyakov D.V. On numerical methods for time-dependent eddy current problems ( for the Maxwell equations in inhomogeneous media. Uchenye Zapiski Kazanskogo \ Universiteta. Seriya Fiziko-Matematicheskie Nauki, 2018, vol. 160, no. 3, pp. 477-494.
(In Russian)
