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

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

Том 5, № 6, 2000

ОБ АЛГОРИТМАХ РЕШЕНИЯ УРАВНЕНИЙ МАКСВЕЛЛА НА НЕСТРУКТУРИРОВАННЫХ СЕТКАХ*

Э.П. Шурина, М.Ю. Великая, М. П. Федорук Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: mife@net.ict.nsc.ru

The methods for solving the time-domain two-dimensional Maxwell equations are presented. These methods are based on the finite volume and the finite element schemes with Lagrange multipliers. Some test results are described.

Введение

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

д E 2 „ 1 .

— - c rotB =-----------j,

at s0

(1)

дВ

~di

+ rotE = 0,

Р

divE = —

So

(2)

(3)

divB = 0,

(4)

где Е — вектор напряженности электрического поля; В — вектор магнитной индукции; j — плотность токов проводимости; р — плотность зарядов в области; с — скорость света в вакууме; е0 — электрическая постоянная; — магнитная постоянная, причем константы удовлетворяют соотношению £0^0с2 = 1. Плотности тока и заряда в этих уравнениях считаются заданными функциями координат и времени.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №98-02-17115а.

© Э.П. Шурина, М.Ю. Великая, М.П. Федорук, 2000.

При реализации таких алгоритмов возникает принципиальная трудность, связанная с нарушением дискретного аналога соотношения (3). Это обусловлено тем, что плотности заряда р и тока j вычисляются по алгоритмам, не зависящим от уравнений Максвелла и друг от друга. Значит, закон сохранения заряда

(& ) + ^ = 0 (5)

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

В литературе описан ряд приемов, обеспечивающих выполнение дискретного аналога закона сохранения заряда в основном применительно к алгоритмам метода частиц на декартовой сетке [7-10]. Так, в работе [7] для решения уравнения (5) предложено на каждом временном шаге вносить поправку в продольную часть электрического поля. Это приводит к необходимости на каждом шаге по времени решать уравнение для поправки к электростатическому потенциалу:

А8^ = ^у^Е — —.

£о

В [8] описан очень простой способ коррекции разностного закона Гаусса, основанный на введении в уравнение (1) члена, пропорционального величине УГ, который получил название псевдотока (^ = ^уЕ — р/е0). В [9, 10] разработан алгоритм, в котором для выполнения разностного аналога уравнения непрерывности предложено находить вклад в плотность тока от каждой макрочастицы, непосредственно подсчитывая потоки через стороны (грани) прямоугольных ячеек. Поскольку этот закон автоматически выполняется для отдельной макрочастицы, он справедлив и для всей системы в целом. Обобщение данного подхода на случай двумерной неструктурированной треугольной сетки дано в работе [11].

В [1, 2] развит еще один подход к решению уравнений Максвелла. В его основе лежит использование неопределенных множителей Лагранжа, благодаря которым уравнения (3),

(4), накладывающие ограничения на дивергенцию полей, включаются в решаемую систему совместно с уравнениями (1), (2). Для этой цели в уравнения (3), (4) вводятся “фиктивные” переменные — множители Лагранжа <^(х, £) и р(х, £), которые осуществляют корректировку электромагнитных полей, обеспечивая тем самым выполнение уравнений (3), (4). При этом множители Лагранжа удовлетворяют однородным условиям Дирихле на границе области.

В настоящей работе описан алгоритм для численного решения двумерных уравнений Максвелла на неструктурированной сетке, основанный на идеях, предложенных в работах [1, 2]. Рассмотрены два метода записи исходных дифференциальных уравнений в дискретной форме: метод конечных элементов и метод модифицированного конечного объема [12] и приведены результаты некоторых тестовых расчетов.

1. Вариационная формулировка задачи для волновых уравнений

Перепишем систему уравнений (1) - (4) в форме, включающей неопределенные множители Лагранжа [1]:

— У^ — с2го1В = —(6)

ОЬ £о

д В

—— Ур + го1Е = 0, (7)

а1уЕ = (8)

£0

а1уВ = о. (9)

Очевидно, что при выполнении закона сохранения заряда с учетом краевых условий

Дирихле для множителей Лагранжа, функции <^(х, Ь) и р(х,Ь) тождественно равны нулю во всей области и полученная система уравнений (6) - (9) эквивалентна исходной системе уравнений Максвелла (1)-(4). Систему уравнений (6)-(9) запишем в форме волновых уравнений второго порядка:

я2е 1 д 1

—— + с2го1го1Е — УФ =-----—, (10)

дЬ2 £0 дЬ

д2В 1

— + с2го1го1В — УР =— го1 (11)

оЬ2 £о

а1уЕ = —, (12)

£0

а1уВ = о, (13)

8^ др где = дь; = дь.

Начальные и краевые условия для системы уравнений (10)-(13) формулируем следующим образом:

Е(Ь = 0) = Е0, В(Ь = 0) = В0,

(Ь = 0) = У р(Ь = 0) + с2гсШо — —1(Ь = 0),

дЬ £0

д В

— (Ь = 0) = Ур(Ь = 0) — го1Ео, а условия на границе области Г:

Ф = 0, Р = 0, Е х п = 0, го1В х п = ^0} х п.

Вариационную постановку для волновых уравнений определим в слабой форме Галер-

кина, для чего вводим функциональные пространства

Н(го1, П) = {V € Ь2(П)3, го^ € Ь2(П)3},

н(а1у, п) = {V € ь2(п)3, € ь2(п)3},

Н*(П) = {р € Ь2(П), Ур € Ь2(П)}.

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

Y = H(rot, П) ^ H(div, П),

Yo = {F Є Y, F x n = 0 на S},

Z = H1 (П)3,

Z0 = {F Є Z, F x n = О на S}.

Вариационная постановка для электромагнитного поля в слабой форме с учетом введенных обозначений имеет вид [l]: для магнитного поля

найти B Є Z, P Є L2(П), являющиеся решением уравнений

-2B VdQ + c2

I VB : VVdn + I PdivVdn + V I (VBa x n) ■ (u« x V)d7

Jo Jo 1 Jr

./п дЬ2 ./п ./п — ^

= — [ го^ ■ У^П +-------[ (1 х п) ■ УdY, УУ € 2,

£о ,/п £о ./г

[ а1уВд^П = 0, Уд € Ь2(П), п

дВ

В(Ь = 0) = Во, (Ь = 0) = Ур(Ь = 0) — го1Ео;

для электрического поля

найти Е € 20, Ф € Ь2(П), являющиеся решением уравнений

д 2E

■ VdQ + с2

дt2

з

[ VE : VVdQ + [ ФdivVdП + V [(VEa x n) ■ (u« x V)d7

,/n Jo a=1 Г

П ,/п

— I I1 ■ У^П + — I ра1уУ^П, УУ € 2о,

£о ип дЬ £о Уп

[ а1уЕд^П = — [ рqdП, Уд € Ь2(П),

п £0 п

Е(Ь = 0) = Е0, (Ь = 0) = У<р(Ь = 0) + с2го1В0-------1(Ь = 0).

ОЬ £о

Известно, что такого рода задачи корректные, если пространства, в которых ищутся решения, удовлетворяют условию Ладыженской — Бабушки — Бреззи (ЛББ) [11]:

Зр> 0 : эир *Уу^П > в||д||ь2(п), Уд € Ь2(П).

VeY ||У||у

В нашем случае два пространства (У, Ь2(П)) удовлетворяют этому условию [1] и вариационная задача является корректно поставленной.

Рассмотрим двумерный случай, когда векторы полей имеют только две ненулевые компоненты. Не нарушая общности, положим ненулевыми компоненты х, у векторных полей Е = (Ех,ЕУ), В = (ВХ,ВУ), а в качестве пробной возьмем функцию вида У = (УХ,УУ),

функции Ух, Уу выберем из одного пространства и обозначим их одним символом “ У ”. Вариационные постановки для компонент х, у имеют вид: для магнитного поля

[ ^ 0Х УdП + с2 [ УBxУУdП+ [ Pд—dП =— [ гotxjУdП +--------------[ (1 х п)хУdY,

Уп дЬ2 ,/п ]п дх £о и п £о Л х

[ _ 0У УdП + с2 [ УВУУУdП+ [ P-—dП = — [ гotyjУdП^--------------[ (1 х п)уУd/y,

Уп дЬ2 ,/п Jn дУ £о Уп £о Л У

Г /ЗВ + сШЛ qdп = 0,

./п \ дх ду )

для электрического поля

УExУУdП+ [ Ф-—dП = — — [ УdП + — [ рд—dП,

п х £0 п Ь £0 п х

УEyУУdП+ [ Фд—dП = — — [ VdП + — [ рд—dП,

( п ) у £0 п Ь £0 п у

Г /^ qdП = 1 [ рqdП.

п х у £0 п

Шаг по времени АЬ считаем постоянным. Электрические поля вычисляем в целые моменты времени Ьп = пАЬ, а магнитные поля — в полуцелые моменты времени Ьп+1/2 = (п + 1/2)АЬ. Множители Лагранжа находим в текущий момент времени Ьп+1/2 для Р и Ьп для Ф. Считаем также, что триангуляция расчетной области задана в виде множества узлов хк = (Хк ,Ук) и множества треугольников еп.

Выбор дискретных конечномерных подпространств необходимо осуществить так, чтобы соблюдалось условие Ладыженской — Бабушки — Бреззи [11]. Для этого достаточно, чтобы степень функций из конечномерных подпространств 2), 2) и V), различалась на единицу. Введем дискретные конечномерные подпространства

2) = {У) : УПп €0н,Ук € Р1},

2° = {У) € ^У) х п = 0Ухк € Г)},

Ь2 = {qh : УПп € ЧН € Р0}

которые определены таким образом, что функции из подпространств 2), 2°) принадлежат классу кусочно-линейных функций, а функции из подпространства V) являются кусочнопостоянными. Такой выбор подпространств обеспечивает выполнение условия ЛББ.

Помимо этого, при использовании базисных функций первого и нулевого порядков полученные в результате пространственной дискретизации матрицы систем линейных алгебраических уравнений имеют разреженную структуру. Значения физических переменных Е, В аппроксимируем линейными функциями, тогда как множители Лагранжа, являющиеся вспомогательными, определим в классе кусочно-постоянных функций. Любую функцию У) € 2) можно представить в виде линейной комбинации финитных базисных функций Ф; : У) = ^ У^ = УТФ, где У; — значения функции в узлах дискретизации, Ф; — кусочно-линейные базисные функции пространства, характеризующиеся свойством Ф;(хк) = 5гк ($гк — символ Кронекера). Любую функцию Ч) € V) можно представить аналогично, в виде линейной комбинации базисных функций ф;: ди = д;ф; = qTф, где —

сеточные значения функции; фг — кусочно-постоянные базисные функции пространства

Ь^, обладающие свойством

фг(хк) =

1, хк € эирр(фг) 0, хк € эирр(фг)

(под эирр(фг) понимается множество всех треугольников, являющихся носителем фг). Дискретный аналог вариационной постановки для магнитного поля имеет вид

МБ ( п=1 I,

1 фтФdП вп+1/2 +

_*/ Оп

МБ

рп+1/2 I = ^

п=1

{Аг2

фт

_*/ °п

Вп-1/2 _

с2 С 43 [> т 1> Бп—1/2 1 [ фтФdП

_0 °п х Аг2 _0 °п

вп—3/2+

+1 [ фт^п+1/2^П + - I фт(jn+1/2 X п^7 £0 лоп £0^п

(14)

_±_

^ | Аг2

п=1 ^

фт фdП

Вп+1/2 + У '

МБ

р"+1/2 И £{ Аг*

п=1

фт фdП

Вп-1/2 _ У

с2 С 43 1> т [> Бп—1/2 1 [ фтФdП

_0 °п у Аг2 _0 °п

В п-3/2+ У

+ 1 [ фтМуГ+1/2^П + - / фт(jn+1/2 X п)у<1ч\ , £0 ,/пп

МБ Г Г дфт

£/.

<9ж

фйПБ,

■п+1/2

+

дфт

дУ

ф(т^п+1/2

0

(15)

(16)

(ЖЕ — число треугольников в триангуляции, ф — локальные базисные функции). В векторно-матричной форме уравнения (14)-(16) принимают вид

М

А2

М

Вп+1/2 + QT р п+1/2 = ррп+1/2 Бп+1/2 + QT р п+1/2 = рт+1/2

где М

ь =

Qx

ф^- ф^П

Оп

М

г,^'=1

Аг2 у ^у

QxБn+1/2 + Qy Буп+1/2 = 0

М

Х^Х

|Пп|

3

матрица массы;

<9ф^- <9фг + <9ф^- дфг

М

г,^'=1

<9ж <9ж ду ду

М

5ф,

<9ж

фг

Ь I dП } — матрица жесткости;

г,^'=1

ф

М

Q

водных;

г,^=1

п„ \ дУ

грп+1/2 = ( 2М _ с2 ь | вп-1/2 _ М Бп—3/2 + рп+1/2

^ х I Аг2 с ь I Бх Аг2 Бх + Р х ;

фг ) dП > — матрицы первых произ-

г,^'=1

П

П

п

п

п

п

П

п

р п+1/2

У

Т ±.

ф (ІП + —

Ф (ІП + —

(5п+1/2 х п)Тф,

-[ ^Г+1/2 х п)ТФгі7|

£0іїп )

векторы правых ча-

стей. Для электрического поля дискретная вариационная постановка имеет аналогичный вид.

На заданной триангуляции матрицы Ь, Qx, Qy имеют одинаковый портрет и разреженную структуру, что позволяет применять экономичную схему хранения — разреженный строчный формат. Использование приближенной, а не точной аналитической формулы для вычисления вкладов в матрицу массы позволяет привести ее к диагональному виду, что дает возможность применить эффективные алгоритмы для нахождения решений полученных матричных систем уравнений без операции обращения матрицы.

2. Конечно-объемная аппроксимация

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

область с некоторой погрешностью аппроксимации ( и 0у ). Будем исходить из того,

что определена триангуляция области . Тогда контрольным объемом 0у, окружающим ]-й узел, является область, ограниченная участками медиан [13]. Дискретизация области на такие контрольные объемы (дуальная сетка) представлена на рис. 1.

Рис. 1. Дискретизация фрагмента расчетной области (а), граничный узел (б) и конечный элемент еп (в). Штриховой линией показаны контрольные объемы, принадлежащие внутренним (Э,Р, Я, т, *) и граничным (к, I) узлам; 1-3 номера узлов элемента еп.

Введем следующие обозначения: 0у — контрольный объем узла ]; еп — треугольный конечный элемент п; Б у — граница контрольного объема; БЩ — часть границы Б у, образованная медианой, выходящей из узла г элемента еп. Зафиксируем некоторую внутреннюю

а

6

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

Г Б'П+1/2 Г дРп+1/2 Г вП-1/2

к~Чп, =2' — *+

I I Бп—3/2

+с2 divVBn-1/2dn - х 2 ^ + — I rot:rjn+1/2dn, (17)

./Пі ./Пі д

Пі д

1 f

єо Пі

Rn вУ

/* B п+1/2 /* dP n+1/2

L, -ЪтdQ -L, —dn=2k^~dQ+

r r Bn-3/2 І Г

+c2 I divVByp^dfi - y dQ + — rotyjn+1/2dQ, (IS)

д єо Jn,

П і

Іг V Пі w t/ Пі

дх ду

Аналогичные выражения для электрического поля имеют вид

г еп+1 Г дФп+1 Г Еп

к л я -к—ш = 2к дхж+

г г Еп-1 1 Г д1 п+1 г2 дрп+1

+с2 ^уУ-М - -VdП-------1-dП----------------^dQ,

Уп4 х к, Дг2 £0]йг дг £0 дх ’

Г Еп+1 Г дфп+1 Г Еп

/й, -Дртd п -/й, —лп = 2.к ДУ2 ^+

+с2 [ - [ —-2Гdtt - — [ %+- dn - —dtt,

4 у 4 Дг2 £о]йг дг £о ду ’

Найдем вклады от треугольного конечного элемента в глобальные матрицы конечнообъемной системы линейных алгебраических уравнений, получаемой из уравнений (17)-(19).

Матрица конвективно-диффузионного члена L

Используя формулу Остроградского — Гаусса, получим

[ divVBxdQ = [ (VBx)ndj, (20)

JQi J Гi

dBx

где (VBx)n = —----производная по направлению внешней нормали. Если т — касатель-

ная к некоторой границе Г и а — угол между касательной и положительным направлениями оси x, то cos(n,x) = — sin а, cos(n,y) = cos а и интеграл, входящий в уравнение (20), запишется так:

[ dBx , dBx ,

J (VBx)„dY = J -gxrdy — -gydx

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

[ ауУВ^П, « [ ——хdy-----------—хdx. (21)

{ ^ А" дх дУ

{щвп содержит узел г} пЬ пЬ

Здесь внешняя сумма означает учет вкладов всех треугольников, примыкающих к узлу г, а пЬ — локальные номера узлов треугольника с номером п, смежных с узлом г. Таким образом, каждая внутренняя сумма состоит из двух слагаемых — потоков через грани Б'Пь, пЬ £ {1, 2, 3}\{г}. В методе конечного объема пространственные интегралы заменяются интегралами по границе и задача состоит в построении аппроксимации потоков через грани. Тогда в предположении некоторого профиля искомых величин на гранях контрольного объема выражение (21) дает дискретный аналог для узла г, который связывает значения искомой функции в узле г и смежных ему узлах триангуляции ],т, ..., д. Получим этот дискретный аналог для внутреннего узла г. Рассмотрим некоторый треугольник (г,],т), содержащий узел г (рис. 2).

Узлы триангуляции (г,],т) имеют локальные номера 1, 2, 3 (после перехода к сборке матрицы и вектора правой части по конечным элементам это ограничение станет несущественным). Вклад треугольника (г,],т) в выражение для потока через границу контрольного объема Бг содержит потоки через два отрезка БП, БЩ. Используем следующие обозначения:

(Хпь, Упь) — основание медианы, выходящей из узла пЬ, пЬ £ {2, 3};

(хс,ус) — координаты центроида треугольника;

1пЬх = ХпЬ - хс, 1пЬу = У)пЬ - ус — длины проекций отрезка медианы БПЬ на оси х, у;

У = кпьх + Ьпь — уравнение медианы, выходящей из узла пЬ.

С дВ

Пусть 3’ппЬ = - х — поток через отрезок медианы Бпь, пЬ £ {2, 3} (в индекса-

18п дп

пЬ

ции потока 3’ппЬ 1 означает, что рассматривается контрольный объем узла с локальным номером 1). Построим аппроксимацию потока 3{Лъ. Согласно формуле (21), получим

Г [ХпЬ дВх(х, кпьх + Ьпь) , ('УпЬ д—х((у - Ь,пъ)/кпъ, у) , 1 п

=\ I ------------Гу-------* ~1, -------------Щ)----------М01*

где коэффициент 01пъ определяется как

01, пЬ = sgn

1 Х1 У1

1 хпЬ упЬ

1 х{2,3}\{пЬ} У{2,3}\{пЬ}

(22)

В каждом из контрольных объемов выберем положительное направление обхода. Для удобства вычислений примем, что интегрирование производится от центра масс (хс, ус) до основания медианы (хпъ,упъ), но соответствующие интегралы входят в балансовые соотношения со знаками “+” (01пъ = 1) или “ -” (01,пъ = -1) в зависимости от того, совпадает ли направление интегрирования с положительным направлением обхода. На практике не придется вычислять определитель (22), поскольку, как показано в [12], согласованные ориентации конечных элементов и контрольных объемов связаны между собой и, таким образом, определенный способ хранения номеров узлов конечно-элементной сетки позволяет избежать дополнительных вычислений для учета положительного направления обхода. Для получения окончательного вида аппроксимации 3ппъ необходимо предположение

Рис. 2. Локальная нумерация вершин треугольника еп и отрезков медиан.

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

Вх(х, у) = < + а\х + а^у, (х, у) £ еп,

значения коэффициентов могут быть определены обращением матрицы Б71 через значения в узлах треугольника:

"ап" ' В1' К 3 II "ь ) 1—1“ ь 3 II ■ 1 х1 у1"

а1 = Кп вх 1 х2 у 2

а2 Вх3 1 хз Уз _

(23)

где Вх = (вх,вх,—х)т — локальный вектор значений функции в узлах треугольника.

Аппроксимация частных производных на элементе еп имеет вид что поток 31ппъ может быть вычислен приближенно:

дВх

дх

а2,

дВ

ду

х п

- — а3, так

(*хпЬ

а^/х -

а^у\ 0\,пъ = {ап1пъх - а21пъу} °1,пъ = 3™

1пъ■

(24)

Согласно (23), входящие в уравнение (24) коэффициенты ап,аЛ выразим через компо^

пп 2 , а3

ненты матрицы Кп и локального вектора неизвестных Вх = (Вх, Вх2, Вх)

т

3

ап =

Екп Вт ап = кп Вт

к2тВх , а3 =2.^ к3тВх

т=1

т=1

Тогда

Тп = 31пъ

1пъ,^2 кптВ? - и £ кптВ'т\ = (!п,В„ пЬ £ {2, 3},

(25)

т=1

т=1

где вектор-строка (1пъ состоит из коэффициентов, с которыми входят в выражение (24) компоненты локального вектора неизвестных Вх = (Вх, Вх2, Вх3)т. Компоненты (1пъ имеют вид

^Хпъ)? = {к3т 1пъх - к2т1пъу}01)Ггъ, т = 1, 2, 3- (2б)

Выражения (25), (26) представляют собой аппроксимацию потоков 312, 3^3 через часть границы Б2Л; У Бп контрольного объема Пг, принадлежащую элементу еп. Любой треугольник е содержит части трех контрольных объемов. Потоки через часть границ контрольных объемов, принадлежащую еп, войдут в балансовые соотношения контрольных объемов Тг,

с

с

У

3

Пу, Пт, а коэффициенты их аппроксимаций дадут вклады в строках и столбцах г,], т. Поэтому вклад треугольного конечного элемента (г,],т), соответствующий аппроксимации интеграла в (21), в матрицу Ь можно представить в виде локальной матрицы 1еп

($1

$П

$П

где каждая из вектор-строк (и = 1, 2, 3) получается суммированием двух вектор-строк,

соответствующих отрезкам медиан, окружающим узел V:

$п

пЬє{1,2,3}\{^}

V = 1, 2, 3.

Согласно формуле (26), при генерации локальных вектор-строк $Лпъ, где V, пЬ = 1, 2, 3, необходим учет положительного направления обхода, выбранного в каждом из контрольных объемов ПДля получения согласованной ориентации всех контрольных объемов достаточно согласовать ориентации конечных элементов [13]. После выбора порядка обхода треугольников, например левостороннего, аппроксимация потока через часть границы бп Р| Бп контрольного объема Пг примет вид

Г дВх

п S■§ дх

$8 & —

г-У 2

аП $х —

аП $у\ +

г у 3

аП $х —

аП $У

( — к3тІ2х + к2тІ2у + к3тІ3х — к2тІ3у) = $п В:

т=1

Здесь вектор-строка определяется как

к31( — І2х + І3х) + к21( — І3у + І2у) к32( — І2х + І3х) + к22( І3у + І2у) к33( І2х + І3х) + к23( — І3у + І2у)

Т

Сделаем некоторые упрощения с учетом того, что

хп хк ^ (хк хп),

Іпх — Ікх = хп — хк = 2(Хк — Хп)’

Таким образом,

$П

1 пу — Іку = Уп — ук = 2(У к — Уп )'

0.5к31 (Х2 — Х3) + 0.5к21 (У3 — У2)

0-5к32(Х2 — Х3) + 0.5к22(У3 — У2)

0.5к33 (Х2 — Х3) + 0.5к23(У3 — У2)

Т

Вектор-строки ,$п можно получить из (27) циклической перестановкой индексов.

І

п

с

с

с

с

у

у

Матрица массы М

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

(Пг Р| еп — часть контрольного объема, принадлежащая конечному элементу еп). Нетрудно

Матрицы первых производных Qx, Qy

Операторы первого порядка аппроксимируем, используя предположение о линейности подынтегральной функции. Тогда

где К'п — вторая строка из матрицы коэффициентов (23); Вх — сеточные значения функции. Элементы локальных матриц имеют вид

(Ь1} Ь2, Ь3 — Ь-координаты на конечном элементе еП). Как показано в [11], элементы локального вектора правой части можно вычислить по формуле

заметить, вид т.= ■

^ У І 0 еп

[ кп<т = |П гП епКПВх =1 ЫкпВх

У.П 3

Вектор правой части

Рассмотрим аппроксимацию вектора правой части в общем виде

На каждом элементе представим функцию f (х) линейным интерполянтом

3

f (х) = X! f (Хк)Ьк

к=1

Дискретизация уравнений (17)-(19) методом конечных объемов с использованием линейных и кусочно-постоянных базисных функций на треугольниках дает три матричных уравнения:

М Вх

_^п+1/2 _ Q рп+1/2 _ рп+1/2

м

м2

рп+1/2 _ Q рп+1/2 _ рп+1/2

QxBn+1/2 + QyВп+1/2 _ 0

1,

где м _ < ~ІЄпІ5із

I 3 ) і,і=1

матрицы первых производных;

N

— матрица массы; Qx _ ^ —1^1^

N

і,3 = 1

Qь

рп+1/2 _ I 2М _ г.2т\ Вп-1/2 _ М Вп-3/2 + рп+1/2.

гх \ Д^2 т І Вх Д^2 Вх р X .

N

і,3=1

рп+1/2 _ / 2М - г2т\ Вп-1/2 - М вп-3/2 + рп+1/2.

РУ V Д^2 с т I Ву Д£2 ВУ + ру .

Рх+1/2, Рп+1/2 — векторы правых частей [12, 13]. Для электрического поля система уравнений имеет аналогичный вид.

3. Алгоритмы решения матричных систем линейных уравнений

Алгоритм Удзавы

Рассмотрим алгоритм Удзавы [14], ориентированный на решение задач о нахождении седловой точки. Пусть Н1, Н2 — конечномерные гильбертовы пространства со скалярным произведением (•, •), которое определяется в зависимости от класса рассматриваемых функций. Сформулируем абстрактную задачу о седловой точке: найти X Є Н1, У Є Н2, являющиеся решением матричной системы

' м Qт 1 X Е

_ Q 0 У О

Здесь Е Є Н1, О Є Н2, М : Н1 ^ Н1 — линейный, симметричный, положительно определенный оператор, QT : Н2 ^ Н1 — линейный оператор, сопряженный с Q : Н1 ^ Н2. Применяя блочное исключение, из (28) получим

QM-^тУ _ QM-1Е _ О.

Очевидно, что матрица QM-1QT симметрична и положительно определена, так как матрица М — симметричный, положительно определенный оператор. В [14] показано, что

(Ям-^т V. V) _ ^ (мци!

Следовательно, необходимым и достаточным условием единственности решения (28) является выполнение условия Ладыженской — Бабушки — Бреззи. Пусть заданы начальные приближения Х0, У0. Определим последовательность {(Хг, Уг)}, г = 1, 2, ... , п [14]:

Хг+1 = Хг + М-1(Е - (МХг + ЯтУг)), уг+1 = т(№+1 - о)

(т — некоторое действительное число). Рассмотрим применение алгоритма Удзавы на примере уравнений для магнитного поля. Волновые уравнения Максвелла эквивалентны трем матричным уравнениям на каждом временном шаге:

д^2 Вх + Ях Р = Рх, (29)

д£г Ву + Яу Р = Ру, (30)

ЯхВх + Яу Ву = °. (31)

Систему (29)-(31) представим в виде

М 0 ЯТ' Вх Рх

0 М ЯТ Ву = Ру

1 Я Яу 0 Р 0

(32)

где М = М/Аі2. Введем обозначения

М

М 0

0 М

Я = [Ях,Яу], В

Вх , Б = Рх

Ву .Ру.

и систему (32) запишем в виде

' М Я ' В ' ' Р '

Я 0 Р 0

(33)

(матрица М — диагональная). Как отмечено выше, при использовании функций первого порядка для аппроксимации физических переменных и функций нулевого порядка для множителей Лагранжа условие ЛББ выполняется. Следовательно, для решения системы (33) можно применить алгоритм Удзавы. Пусть заданы начальные приближения ВХ,В0,Р0. Определим последовательность

ВХ+1 = ВХ + М-1(Р* - (МВХ + ЯТРг)), ву+1 = Ву + М-1(Ру - (Мву + ЯТРг)), Рг+1 = Рг + т (ЯхВх + Яу ву),

где т — некоторое действительное число. Так как матрица массы диагонализирована, вычисление обратной к ней матрицы не требует больших затрат.

Регуляризация по Тихонову

Рассмотрим при фиксированном моменте времени систему матричных уравнений (29) -

(31), полученных конечно-элементной аппроксимацией. Из первых двух уравнений выра-

зим Вх, Ву:

Вх = (Д|) (-°ХР + Рх) • (34)

= (^1) ‘(-ОХР + Р») • (35)

Полученные выражения подставим в уравнение (31):

°х( (-ОХ Р + Рх) + °у(" (-°Х Р + Ру) =° (36)

и разрешим (36) относительно множителя Лагранжа:

°х( Д» )-1°т Р+°у (Д» )-1°т Р = °х ()"1рх + я (Д» )"Ч

Так как матрица М диагонализирована, она перестановочна с матрицами Ох,Оу. Тогда уравнение относительно Р примет вид

АР = Р, (37)

где А = (ОхОХ + ЯуОХ), Р = ОхРх + ОуРу. Решив (37), из (34), (35) получим Вх,Ву. Однако система (37) плохо обусловлена, и задача нахождения множителя Лагранжа требует регуляризованного алгоритма [15]. Найдем решение Р£, доставляющее минимум функционалу

Ма[Р,Р] = ||АР - Р||2 + е||Р||2. дМа [Р, Р]

Из условия экстремума ---------- = 0 для определения регуляризованного решения Р£

дР

получим

(АтА + е/) Р£ = АтР. (38)

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

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

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

В первом тесте [4] в начальный момент времени задавалось распределение ^-компоненты магнитного поля

Вх(х, 0) = В0 сов(аж) сов(бу), (39)

тП А ПП -Р О 1 1 і п

где а = ----, Ь = —. В расчетах полагалось В0 = 1,п = т = 1,Х0 = у0 = 1- В расчет-

Хо Хо

ной области вводилась равномерная квадратная сетка с шагом Н = 0.01, каждая ячейка которой делилась диагональю на два треугольника для получения треугольной сетки. Результаты вычислений, а также соответствующее аналитическое решение для компоненты электрического поля Ех приведены на рис. 3.

Рис. 3. Зависимость электрического поля Ех от времени (а) и распределение его в сечении у = 0 (б) в точке (0.3, 0.4) при Аі = 4 • 10-10 е, Н = 0.1 (слева) и Аі = 4 • 10-11 е, Н = 0.01 (справа).

(—) — аналитическое решение; (□) — метод конечного объема;

(о) — метод конечных элементов + регуляризация; (+) — то же + алгоритм Удзавы.

Как видно из графиков, соответствие между результатами численного моделирования и аналитическим решением достаточно хорошее.

Во втором тесте [3] в начальный момент времени задавалось распределение ^-компоненты электрического поля

Ег(х, 0) = Е0 сов(аж) сов(Ьу),

тп пп

где а =------, Ь = —. В расчетах полагалось Е0 = 1, п = т =1, х0 = у0 = 1. Магнитное поле

Хо Хо

имеет ненулевыми Х-, у-компоненты. На рис. 4 представлено аналитическое и численное решения компонент магнитного поля Вх.

Рис. 4. Зависимость магнитного поля Вх от времени (а) и распределение его в сечении х = 0 (б) в точке (0.3, 0.4) при Дг = 4 • 10-10 с, Н = 0.1 (слева) и Дг = 4 • 10-11 с, Н = 0.001 (справа). Условные обозначения см. в подрисуночной подписи к рис. 3.

Заключение

Представленные в статье методы решения нестационарных уравнений Максвелла в физических переменных реализованы и исследованы в серии расчетов высокочастотных электромагнитных полей с высокой и низкой амплитудами при различных параметрах дискретизации области и временного интервала. Сравнение методов показало, что предпочтение следует отдать алгоритму Удзавы, поскольку этот метод обладает лучшей сходимостью и точностью, требует меньшее число операций умножения матрицы на вектор. Полученные результаты хорошо согласуются с известными аналитическими решениями тестовых задач [3, 4], что позволяет сделать вывод о возможности применения данного подхода для расчета нестационарных электромагнитных полей.

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

[1] Assous F., Degond P., Heitze E., Raviart P. A., Segre J. On a finite-element method for solving the three-dimensional Maxwell equations // J. Comp. Phys. 1993. V. 109. P. 222-237.

[2] Assous F., Degond P., Segre J. Numerical approximation of the Maxwell equations in inhomogeneous media by P1 conforming finite element method // J. Comp. Phys. 1996. V. 128. P. 363-380.

[3] Munz C.-D., Schneider R., Vos U. A finite-volume method for the maxwell equations in the time domain. Forschungszentrum Karlsruhe “Technik und Umwelt”. Karlsruhe, 1996.

[4] Hermeline F. Two coupled particle-finite volume methods using delaunay-voronoi meshes for the approximation of vlasov-poisson and vlasov-maxwell equations // J. Comp. Phys. 1993. V. 106. P. 1-18.

[5] Trowbridge C.W., Bryant C.F., Emson C.R. I. Some developments in electromagnetic field computation // Proc. 7th MAFELAP Conf. on the Math. of Finite Elements and Appl., 1990.

[6] Березин Ю. А., Вшивков В. А. Метод частиц в динамике разреженной плазмы. Новосибирск: Наука, 1980.

[7] Boris J. P. Relativistic plasma simulation — optimization of a hybrid code coordinates // Proc. 4th Intern. Conf. on the Numerical Simulation of Plasmas, Washington, Sept. 1970. P. 3-6.

[8] Marder B. A method for incorporating Gauss’ law into electromagnetic PIC codes // J. Comp. Phys. 1987. V. 68. P. 48-51.

[9] Villasenor J., Buneman O. Rigorous charge conservation for local electromagnetic field solvers // Comp. Phys. Commun. 1992. V. 69. P. 306-361.

[10] Березин Ю. А., Брейзман Б.Н., Вшивков В. А. Численное моделирование ин-жекции мощного электронного пучка в вакуумную камеру с сильным магнитным полем // ПМТФ. 1981. №1. С. 3-9.

[11] Coullietteand D.L., Koch М. On the difficulties and remedies in enforcing the div = 0 condition in the finite element analysis of thermal plumes with strongly temperature-dependent viscosity // Intern. J. Num. Methods in Fluids. 1994. V. 18. P. 189-214.

[12] Шурина Э.П., Войтович Т. В. Анализ алгоритмов методов конечных элементов и конечного объема на неортогональных сетках при решении уравнений Навье — Стокса// Вычисл. технологии. 1997. Т. 2, №4. С. 84-104.

[13] Рояк М.Э., СоловЕйчик Ю.Г., Шурина Э.П. Сеточные методы решения краевых задач математической физики: Учеб. пособие. Новосибирск: Изд-во НГТУ, 1998.

[14] Bramble J. Н., Pasciak J. E., Vassillv A. Т. Analysis of the inexact Uzawa algorithm for saddle point problems // SIAM J. Num. An. 1997. V. 34, N 3. P. 1072-1092.

[15] Тихонов А. Н., Арсенин В.Я. Методы решения некорректных задач. M.: Наука, 1974.

Поступила в редакцию 14 февраля 2000 г.