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

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

Том 11, № 6, 2006

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

ОБЛАСТИ*

М.В. ПУЗАНОВ, Э.П. ШУРИНА Новосибирский государственный технический университет, Россия e-mail: misha.puzanov@gmail.com, shurina@online.sinor.ru

М. И. Эпов Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, Новосибирск, Россия e-mail: mepov@uiggm.nsc.ru

Modeling of time-dependent electromagnetic field in heterogenous media comprising insulating and conducting materials is considered. The proposed approach allows to interpret the source supporting conductor as an element of the media.

Введение

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

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

В [1, 2] представлены вариационные формулировки для магнитного поля в квазистационарной задаче о вихревом токе; в [3] дан обзор различных способов введения тока и напряжения в модель. В [4] предложено несколько способов учета тока и напряжения в модели токонесущего проводника, состоящего из множества тонких слоев.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-05-64528) и совместного международного проекта NWO и РФФИ (грант № 047.016.003).

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

Векторный метод конечных элементов [5-7], являясь современным методом расчета электромагнитных полей, представляет собой естественный инструментарий для сопряжения интегральных характеристик, таких как ток и напряжение, со значениями векторных характеристик электрического и магнитного полей.

Рассмотрим модель системы, состоящей из источника поля (генератора напряжения и кабеля, который образуют замкнутый контур, включающий генератор как элемент) и неоднородной окружающей среды, содержащей как непроводящие подобласти (воздух), так и подобласти, в которых a > 0 и е > 0 (земля).

Модель системы характеризуется следующими особенностями:

— источник поля требует (в силу размеров и неоднородности свойств) представления в виде элемента с конечной геометрией и, как следствие, включения источника в область моделирования;

— зависимость возбуждающего тока от времени произвольная;

— отсутствует априорная информация о пространственном распределении плотности тока в источнике.

Известные способы соединения модели электрической цепи с моделью электромагнитного поля используют специфичные вариационные формулировки для магнитного поля в сочетании с ограничениями на пространство, в котором ищется H [3]. Эти дополнительные ограничения (в частности, rot H = 0 для непроводящих подобластей в задаче о вихревом токе) требуют решения условной задачи минимизации с использованием множителей Лагранжа [8] либо выделения из расчетной сетки максимального дерева [9, 10]. Применение множителей Лагранжа предполагает наличие некоторой априорной информации о свойствах аппроксимирующего оператора, которая нужна для управления параметрами метода решения. Поиск же дерева в исходной сетке в трехмерном случае может оказаться достаточно трудоемким. Так или иначе оба эти метода существенно увеличивают затраты на решение задачи.

В представленной работе предложен подход, основанный на методе декомпозиции, использующий постановку Галёркина первого порядка [11] и не требующий условной минимизации функционала, что достигается за счет введения искусственных краевых условий на границе между источником поля и окружающей средой. Эти краевые условия (определение которых тривиально для относительно простой геометрии источника поля) строятся так, чтобы удовлетворять соотношениям, формулируемым для электрической цепи.

Предложенная модель описывает достаточно широкий класс систем, используемых в геоэлектрике.

1. Постановка задачи

Электромагнитное поле описывается системой уравнений Максвелла

rot E = -

dB ~dt

rot H

divB = divD :

dD

J + er-0,

P-

(1)

и материальными уравнениями В = ^И и Б = еЕ. В общем случае 3 = 3° + оЕ, где о — электрическая проводимость; е — диэлектрическая проницаемость; ^ — магнитная

Рис. 1. Область моделирования.

проницаемость; 3 будет интерпретироваться в зависимости от физического контекста.

Уравнения (1) рассматриваются в области П = П0 и П и П2 и П3, где Пг,г = 0, 3, — замкнутые множества; П0 — односвязная область (генератор напряжения); П1 — одно-связная область, такая, что П;п4 = П0 и П1 обладает свойством Вп(П0) = 0 Уп > 0, Вп(П1) = 0 Уп > 0, В1(П;п4) = 1, и Вп(П;п4) = 0 Уп > 1, где Вп(Д) — п-е число Бетти области Д Иными словами, П;п4 — петля с вставкой-генератором. Область Пех4 = П2 и П3 — окружающая среда, в которой П2 — воздух, П3 — земля. Область моделирования представлена на рис. 1.

Возбуждаемое генератором напряжение задается функцией п(£). Параметры среды имеют следующий вид:

' 0,х е П0 и П1, е0, х е П2 и П3,

^(х) = ^0,

е(х)

о(х) =

00, х е П0,

01, х е П1, 0,х е П2, 03,х е П3.

Будем обозначать через Г^ границу между областями П и Пj. Предполагаем, что индуктивностью системы можно пренебречь, поэтому она обладает только активным сопротивлением, следовательно, величина тока в проводнике П;п4 есть

I (*)

п(г)

Л1 + Д0'

где п(£) — заданная величина напряжения. Считаем, что ток во всей системе не терпит разрывов, а следовательно, распределение тока в кабеле можно найти исходя из:

— уравнений Максвелла;

— условия непротекания тока на границе области П;п4 (в воздух и в землю), т. е. условия 3п|ГокиГ1к = 0, к = 2, 3;

— условия сохранения полного тока в цепи J J ■ n = I, где S — некоторое сечение

S

проводника Qint.

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

Рассмотрим гильбертовы пространства

H(rot, П) = {u|u е (L2^))3, rot u e (L2^))3}

с нормой ||u|| = J lufdQ + J |rotu|2dn,

П П

H(div, П) = {u|u e (L2(Q))3, divu e L2(Q)} с нормой ||u|| = J |u12dQ + J(divu)2dn,

ПП

а также

H0(rot, П) = H(rot, П) П {u| u x т|dn = 0}, H0(div, П) = H(div, П) П {u| u ■ n|dn = 0}. Введем в H(rot, П) и H(div, П) скалярное произведение

(«,„) = /«■ „«m.

П

В области Пext магнитное поле определяется уравнением, которое получается из (1) исключением поля E:

dH д 2H т0 rot rot H + ^a——+ ^.е 2 = rot J , (2)

dt dt2

где е, a — константы.

В области П^ модель, связывающая магнитное поле и электрический ток, имеет вид

rot H = J,

, dH

rot (a 1J) = -^-дг,

/J.n = l, (3)

S

I J ■ nUnt = 0-

Здесь S — некоторая произвольная поверхность, являющаяся поперечным сечением области П;п^ Под поперечным сечением понимается поверхность Сейферта, двойственная единственному классу 1-циклов в П;п^

Выполнение условия J ■ n = I можно обеспечить специальными краевыми условиями

Sr

для поля H, для которых должно выполняться ограничение У H ■ d/ = I (т. е. ограниче-

dS

ние, налагаемое на циркуляцию магнитного поля по границе сечения S). Таким образом,

можно ввести краевые условия Н|ап = С(х,у,г), где С(х,у,г) подбирается исходя из соотношения / С(х,у,г) • = I. Эти условия могут быть заданы на части границы (на-

dS

пример, на дП0 — Г01), на остальной же части (в этом случае дП1 — Г01) задаются условия

3п|дП1-Г01 =

Из (3) получим

f rot H = J,

, dH

rot (a J) = —, v ; P dt

H |0По-Гох = C )>

J|

0П1-Г01

0.

Этим двум моделям соответствуют вариационные постановки. 1. Найти H G H(rot, Hext) такое, что W G H0(rot, Hext):

Г Г dH Г д2Н Г

rot rot H ■ V + ■ V + ■ V =

rot J0 ■ V.

(4)

2. Найти H 6 H(rot, Hint) и J G H0(div, Hint) такие, что VV G H0(rot, Hint) и VF G

H0(div, Hint):

Из соотношения

rot H ■ F

J ■ F,

dH

rot (a J) ■ V = — ■ V,

dt

H |дПо-Го1 = C (x,y,z) 0.

, J Uni-ro! J rot u ■ v = J rot v ■ u + J div(u X v),

n int n int n int

а также учитывая, что V G H0(rot, Hint) и, следовательно,

div((a J) X V)= / ((a J) X V) ■ n = 0

окончательно получаем систему

rot H ■ F = / J ■ F,

a-1 J ■ rot V = — I ■ V

dH

dt

n int n int

H |dno-roi = C (x,y,z)

, J Uni-roi = 0.

(5)

n

n

n

n

n

n

n

n

n

n

n

Задача "склеивания" по непрерывности поля Н постановок (5) и (4) решается с использованием метода декомпозиции области с наложением [12, 13]. Далее для удобства будут использоваться термины "внутренняя" и "внешняя задачи", обозначающие соответственно задачи (3) и (2).

3. Дискретизация

Определим дискретные подпространства W(П) = spanj^ G H(rot, П), i = 1, N}, а также = W(П) П H0(rot, П), F(П) = spanj^- G H(div, H),i = !,M}, F°(H) = F(П) П H0(div, П).

Для внутренней задачи будем использовать следующую схему пространственной дис-

nint

кретизации. Представляя магнитное поле в виде H(t) = вг(^) фг и ток в виде

г=1

m

J(t) = E a,(t) у,, где фг G W(П^), у, G F(П^), nint = dim(W(П^)), m = dim(F(П^)), j=i

приходим к системе уравнений

r nint f m Г

rot фг ■ Ур = Z) a, у, ■ Щр,

г=1

j=1

Oint nint

Ё a, j a 1 у, ■ rot фк = - ^ ^dt^ / ^ Фг ' Фк,

j О. i=1 Q.

O int О int

H |дО0-Г01 = C (x,y,z ),

J|

dOi-roi

0,

где Щр G F0(Пlnt) и фк G W0(Пlnt). В матричном виде (6) примет вид

(6)

Св = Fa,

A a = -Ef, dt

(7)

Срг = J rot фг ■ Щр, Eki = у ^ фг ■ фк,

О

О -1

Fp, = J у, ■ Щр, Ak, = J a у, ■ rot фк.

ОО

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

Свn = Fan,

A (Aan + (1 - A)an-1) = -или в векторно-матричном виде

tn - tn-1

— E (вп - вп-1)

(8)

SnE AA " " вп " " SnEen-1 - (1 - A)Aan-1 "

С -F an 0

1

где sn =-1--, Л е [0,1].

г"! n ¿n — ¿n—1 ' L ) J

Дискретизация внешней задачи выполняется с помощью дискретного аналога слабой

next

формы Галёркина. Пусть H(t) = где ^ е W(Hext), next = dim(W(Hext)), тогда

i=1

дискретная форма постановки (4) имеет вид е W0(Hext):

n

Yi I rot ^ ■ rot + ' + "dt^ I ' = I rot J0 '

next next next

или в векторно-матричной форме

+ M' Ж + M * ^ = f

где

Rij = J rot ^i ■ rot ^j; Mjj = J ^a^i ■ ^j; Mjj = J ^j ■ ^j.

next next next

Вводя обозначения

k 1 = 1 k = - 1

r^n— 1 - /....-....-IN/...,-. .41 r^n -

(tn+1 - tn—1)(tn+1 - tn) ' (tn - tn—1)(tn+1 - tn)

1

к

n+1

(tn+1 - tn—1)(tn - tn—1)' для неявной схемы аппроксимации по времени получаем

[R + Sn+1 + kn+1M*] 7n+1 = fn + Sn+1MCTYn - M* [kn—1 Yn—1 + knYn] • (10)

n

4. Декомпозиция области

Схема декомпозиции области с наложением приведена на рис. 2, где изображено сечение прямого участка петли. Петля разбита на две непересекающихся подобласти - П^ и П^, где П^ представляет собой слой петли толщиной Для решения на г-м временном слое предлагается итерационная процедура, ориентированная на вычисление наведенного в петле тока. При этом внутренняя задача решается в подобласти П^ и П^, а внешняя — в области Пех^ В качестве краевых условий для внешней задачи используется значение поля, полученное при решении внутренней задачи, и наоборот. Подробный алгоритм приведен ниже (через вк обозначим решение задачи (7) на к-м временном шаге, через 7к — решение задачи (9) на к-м временном шаге). Будем полагать начальные значения полей нулевыми.

1. в0

2. y0

3. к = 2.

0.

Y1

0.

вк+1 - вk

4. while к = N V

> £t

5. n = 0.

k

Q.J,

Рис. 2. Декомпозиция области (сечение прямого участка петли.)

6. (ßj

[п], a[n

7. while

= CInt (в

e[fcn+l] — ß[n]

k-1

8. 9.

Yn]

ßN -l

> do.

„k-2 ok >ß[.

Ext (rv ^ [n]

efn+^afn+i]) = Int (ßk-1

7[fcn]

10. n = n +1.

11. done.

12. (ak,ßk,7k) = (afnpßinpYfn]}

13. k = k + 1.

14. done. Здесь:

yM

Ext ^y k 1,y k 2,ß[n^ означает решение краевой задачи (4) с использованием

схемы (10) в области next U при условии, что ток в области (в тонком слое петли) рассматривается как сторонний. Решение задачи включает в себя магнитное поле, рассчитанное в области П^;

— (ß[n+1]' = Int ^ßk-1,Yjn^ — это решение внутренней задачи (5) в области

Пз U n^d с использованием схемы (7) при условии, что в качестве краевых условий для J на каждой итерации берутся условия Jn|r = 0, а в качестве краевых условий для поля Н

берется Н™+1] задачи(;

ext H[n+1]

, т. е. значение H на границе, полученное при решении внешней

— = CInt (yßk ^ означает решение внутренней задачи (5) без наведенного

тока по схеме (7), при этом граничные значения для H задаются из соотношения

I(t) = y J(t)ndS = J H(t) dL,

S dS

где H|r = const. Иными словами, предполагаем, что поле на границе постоянно вдоль каждого контура, охватывающего поперечное сечение петли, а циркуляция его равна полному току, который легко определить из закона Ома. Такие краевые условия налагают

г

более сильные требования, чем ограничение I(t) = J J(t)n dS, однако они верны, если

S

пренебречь влиянием противоположных частей петли друг на друга.

5. Численный эксперимент

На данном этапе наибольший интерес представляет решение внутренней задачи (5), поэтому численный эксперимент проводился для дискретной формулировки (6).

Взяты две параллелепипеидальные сетки, используемые для дискретизации области Oint, первая содержит 9 600 ребер и 8 500 граней (количество ячеек сетки вдоль внутреннего ребра контура 20, число ячеек вдоль стороны сечения 5), вторая — 40 920 ребер и 38400 граней (количество ячеек сетки вдоль внутреннего ребра контура 20, число ячеек вдоль стороны сечения 10). Общая размерность задачи при использовании векторного базиса первого порядка для первой сетки составляет 18 100, для второй — 79 320. Источник имеет геометрические размеры, близкие к размерам реальной физической системы. При моделировании использовались следующие параметры: a0 = о\ = 107, u(t) = e-at sin(2nvt), где v = 500, а = 200. Шаг в схеме интегрирования по времени в схеме (8) выбирался равным

Экспериментально установлено, что подходящим с точки зрения устойчивости процесса моделирования во времени для схемы (8) является значение параметра Л = 0.5. Для

5 • 10-5.

Рис. 3. Краевые условия для поля H.

Рис. 4. J в одном из углов проводника.

Z

0.02

О

-0.02

0.50 0.52 0.54 У

Рис. 5. Сечение плоскостью x = 0.

решения системы линейных алгебраических уравнений на каждом временном шаге использовался метод GMRES [14], размерность подпространства Крылова равна 100, критерием выхода было уменьшение нормы невязки в 107 раз.

На рис. 3 схематично показан вид краевых условий для поля H. На рис. 4 приведен вид поля J в одном из углов проводника для второй сетки. Аналогичная качественная картина тока наблюдается для любого из четырех углов.

Для исследования поведения решения при измельчении сетки использовались значения компонент тока Jx(P1,t), Jx(P2,t), Jx(P3,t) и магнитного поля Hz(P1,t), Hy(P2,t), Hz(P3,t) в точках P1 = (0, 0.52, 0), P2 = (0,0.52, 0.015) и P3 = (0, 0.535,0) (рис. 5).

В качестве критерия точности решения применялся следующий:

u(t)/R «У Jn(t)n,

S

т. е. интеграл от плотности тока по сечению проводника должен с достаточной точностью совпадать с исходной силой тока, используемой для создания краевых условий. В качестве сечения S выбрано одно из П0 П П1, т.е. связная часть границы между "вставкой" и остальным проводником.

В таблице приведены исследуемые значения компонент поля H и плотности тока J в указанных точках, а также значения вычисленной силы тока In(t) на двух сетках. I(t) — исходное значение силы тока. Значения t выбирались случайно, чтобы проиллюстрировать величину ошибки при различных значениях силы тока.

Из таблицы видно, что Hz (P1, t), которая в теории должна быть нулевой, уменьшается с измельчением сетки и точность вычисления H предположительно выше, чем J (меньше разница между соответствующими значениями для сеток M1 и M2) (рис. 6, 7).

Поскольку для аппроксимации уравнений первого порядка в пространствах H (rot) и H(div) используется векторный базис первого порядка, теоретическая скорость сходимости имеет первый порядок. Порядок аппроксимации также первый [6, 7]. Так как в работе приведено решение реальной задачи, для оценки точности использовалась величина полного тока в замкнутой петле In(t) на двух сетках. Точность его вычисления позволяет говорить, что практическая скорость сходимости не ниже теоретической.

Из рис. 8 видно, что согласно проведенным расчетам можно предполагать отсутствие скин-эффекта в проводнике, поскольку значения Jx(P1), Jx(P2) и Jx(P3) можно считать практически равными. Как видно из рис. 9, Hz(P2) и Hy(P3) с большой точностью совпадают, что просто объясняется симметрией сечения проводника. Наличие небольшой ненулевой составляющей Hz (P1) обусловлено погрешностью вычислений. На рис. 10 показаны зависимости исходной силы тока от времени, а также вычисленные на сетках M1 и M2.

Значения компонент H(t), J(t) и I(t)

M1 M2

t I(t) Jx(P1,t) Hz (P1,t) Hz (P2, t) I n(t) Jx(P1,t) Hz (P1,t) Hz (P2, t) I n(t)

0.00015 345.9 -181206 13.1 —1639.8 289.9 -221491 0.61 —1604.2 354.4

0.00275 320.3 -204377 1.57 —1510.1 327.0 -196248 -3.16 —1500.7 313.9

0.0038 -215.8 121518 -9.08 997.1 -194.4 134282 0.07 1024.5 -214.8

0.0045 319.2 -135662 -24.5 -1386 217.0 -197995 2.65 -1488.8 316.7

0.0061 71.62 -54540.6 25.3 -351.1 87.26 -44348.7 2.16 -320.72 70.95

0.0073 -147.5 91213.4 -3.81 684.3 -145.9 98219.6 2.78 684.01 -157.1

Рис. 6. Зависимость Н от времени в точке Р2 на сетках М1 и М2; Н2(Ь),М1 соответствует Нг(í) в точке Р2 на сетке М1; Н2(Ь), М2 — Нг (í) в точке Р2 на сетке М2.

Рис. 7. Зависимость И от времени в точке Р3 на сетках М1 и М2; H2(í),Ml соответствует Ну(¿) в точке Р3 на сетке М1; H2(í), М2 — Ну (¿) в точке Р3 на сетке М2.

500000 400000 300000 200000 100000 0

-100000 -200000 -300000 -400000 -500000

0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008

Рис. 8. Зависимость J от времени в точках Р\, Р2 и Р3 на сетке М2; Jl(t) соответствует ,1Х(Ь) в точке Р1; J2(t) — Jx(t) в точке Р2; J3(t) — Jx (Ь) в точке Р3.

3000 2000 1000 0

-1000 -2000 -3000 -4000

0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008

Рис. 9. Зависимость Н от времени в точках Р1, Р2 и Р3 на сетке М2; Н1 (Ь) соответствует Нх(Ь) в точке Р1; Щ(Ь) — Нх(Ь) в точке Р2; Нз(Ь) — Ну(Ь) в точке Р3.

Рис. 10. Зависимость от времени I(Ь) и 1п(Ь) на сетках М1 и М2; 1п(Ь),М1 соответствует 1п(Ь) на сетке М1; Iп(Ь),М2 — 1п(Ь) на сетке М2; I(Ь) — истинное значение силы тока.

Рис. 11. Зависимость J от времени в точке Р1 на сетках М1 и М2; ,11(Ь),М1 соответствует Jx(t) в точке Р1 на сетке М1; Jl(t), М2 — Jx(t) в точке Р1 на сетке М2.

Видно, что совпадение на сетке M2 лучше. На рис. 6 и 7 показаны зависимости поля H от времени в точках P2 и P3 соответственно, вычисленные на сетках M1 и M2. На рис. 11 представлены зависимости J от времени в точке Pi, вычисленные на сетках M1 и M2.

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

[1] Dular P., Geuzaine C., Legros W. A natural method for coupling magnetodynamic h-formulations and circuit equations // IEEE Trans. Magn. 1999. Vol. 35, N 3. P. 1626-1629.

[2] Buffa A. Some numerical and theoretical problems in computational electromagnetism // Ph.D. Thesis. Univ. of Milano, 2000.

[3] Hiptmair R., Sterz O. Current and voltage excitation for the Eddy Current model // Intern. J. Numer. Modelling. 2005. Vol. 18, N 1. P. 1-21.

[4] Geuzaine C. High order hybrid finite element schemes for maxwell's equations taking thin structures and global quantities into account // Ph. D. Thesis. Univ. of Li'ege, Belgium, Faculty of Appl. Sci., 2001.

[5] Hiptmair R. Finite elements in computational electromagnetism // Acta Numerica. 2002. Vol. 11. P. 237-339.

[6] Nedelec J.C. Mixed finite elements in R3 // Numerische Mathematik. 1980. Vol. 35. P. 315-341.

[7] Nedelec J.C. A new family of mixed finite elements in R3 // Numerische Mathematik. 1986. Vol. 50. P. 57-81.

[8] Alonso-Rodriguez A., Hiptmair R., Valli A. Mixed finite element approximation of eddy current problems // IMA J. Numer. Anal. 2004. Vol. 24. P. 255-271.

[9] Bossavit A. Computational Electromagnetism Variational Formulation, Complementary, Edge Elements. San Diego, London, Boston, New York, Sydney, Tokio, Toronto: Acad. Press, 1997.

[10] Gross P., Kotiuga P. Finite element-based algorithms to make cuts for magnetic scalar potentials: Topological constraints and computational complexity // Geometric Methods for Comp. Electromagnetics / F. Teixeira (Ed). 2001. Vol. 32. Cambridge: EMW Publ., 2001. P. 207-245.

[11] White D.A. Discrete time vector finite element methods for solving Maxwell's equations on 3d unstructured grids // Ph.D. Thesis. Univ. of California, Livermore, 1997.

[12] Stefanica D. Domain decomposition methods for mortar finite elements // Ph.D. Thesis. Courant Institute of Math. Sci. N.Y. Univ., 1999.

[13] Domain Decomposition Methods for Vector Field Problems // Ph.D. Thesis. Courant Institute of Math. Sci. N.Y. Univ., 1999.

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

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