Метод объёмных интегральных уравнений в задачах магнитостатики

Математические методы и приложения

для моделирования сложных систем

УДК 519.64

Метод объёмных интегральных уравнений в задачах

магнитостатики

П. Г. Акишин, А. А. Сапожников

Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, Дубна, Московская область, Россия, 141980

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

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

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

Ключевые слова: магнитостатика, метод объёмных интегральных уравнений, метод конечных элементов, дискретизация, нелинейная проблема, итерационный процесс.

1. Интегральные постановки задачи магнитостатики

Пусть Н(а), В (а), ММ (а) есть напряжённость, индукция и намагниченность магнитного поля в точке а . Имеет место следующее интегральное уравнение:

Н(а) = Нь (а) + — Уау (М(х), Уа ^а—щ) ^ ^

а

где Н^(а) — поле от токовых обмоток, О — область, заполненная ферромагнетиком. Величины Н, В, ММ связаны следующими нелинейными соотношениями:

Н

где ^о — абсолютная магнитная проницаемость вакуума, ^(х) — магнитная проницаемость.

Статья поступила в редакцию 30 декабря 2013 г.

2. СЕ^ метод

N

Пусть область С разбита на подобласти С = У С^. Мера пересечения ^

г=1

/ xAva:

0,г = ]. а,г = | . Положим В(х) в каждом постоянным и равным В^. Тогда [1]:

N

V«

Нг=н з ^

3 = 1

Ка> • ^)

с,

, г = 1,М,

где В^ и М^ удовлетворяют (2). Пусть

в = {В1,В2,...,вну , М(В) = {^М(В1),^оМ(В2),...,^М(Вн))т,

Н 3 = (/юН 3 (а1),^оН8 (а2)1...1^оН3 (ам ))Т. -[¿•и] ••• [Аш У

Пусть матрица [Л] есть [Л]

[Ам 1] • • • ].

. Для любого постоянного

вектора М справедливо соотношение:

V«

[Агз]М = ^

I (^М, Уа ^ - ^

с,

Учитывая (2), систему (3) можно записать в виде:

В = Н3 + ([Л] + [Е]) М(В), где [Е] — единичная матрица размерности [3К х 3Щ.

3. Метод усреднения по элементу

Вариацией СЕЙМ метода будет замена точки наблюдения йг в С^ интегрированием по [2]:

Нг 1= IН3(а№а + ^ У ^ I (м,, V а ^

Сг Сг 3 ^Ог С,

г = 1,К.

Пусть

В= (В1,В2,...,ВМ )Т, М(В) = (^оМ(В1),^М(В2),...^оМ(Вм))'

Н3 = \iioJ Н3(а)6ьа,^о у Н3(о)ёиаН3(а)ёг>а

Ог 02

[4

[¿11]

[ Аш ]

_[ An 1]    [ Ann ]_ Для любого постоянного вектора Я справедливо соотношение:

[А,ЛМ = 1<Ъ. V ¡(я, V. ^ iv,

Gi Gj

Пусть [ С] — диагональная матрица размера [3 N х 3 Щ] с элементами по диагонали (а, с1, с\, с2, с2, с2,..., с^, с^, с^), где ^ = / ё V. Учитывая (2), систему (5) можно

Gi

записать в виде:

[ С]В = Н 3 +

4. Кусочно-линейная аппроксимация намагниченности

Интегрирование по элементу позволяет использовать аппроксимацию намагниченности, отличную от постоянной, например, линейную. Рассмотрим линейную конечно-элементную аппроксимацию намагниченности в пределах объёмного элемента. Для этого будем предполагать, что область можно приблизить объеди-

N

нением тетраэдров: С = и

¿=1

Потребуем, чтобы тетраэдры {5^} удовлетворяли следующим условиям:

1) мера пересечения S^ с Sj равна нулю при г = ]]

2) вершины одного тетраэдра не могут быть внутренними точками грани или ребра другого тетраэдра, т.е. если два тетраэдра касаются, то они касаются или только по целому ребру, или только по целой грани, или только по одной вершине.

Пусть Р^, к = 1,.., Ь — набор всех вершин тетраэдров {5^}. Обозначим Н(Рк) = Нк, М(Рк) = Мк, В(Рк) = Вк. Пусть /к(х) — функция формы, ассоциированная с вершиной Рк:

_ Г 1, если к = I, Л( Рг) = \0, если к = I.

Функция /к(х) на каждом тетраэдре есть линейная функция. Используя эти обозначения, имеем:

N

Y, I fi(й)Нdv, = j («)dv, +

J=1 G G

N

G

v„

+£/Л(a)V Jи,V^

J=1G Lg

dv,, i = 1,N.

Пусть матрица [С] есть матрица размерности [3L х 3 L] :

'[Си] ••• [СШ ]

[ С] =

L[ См 1]

[Cnn ].

где — диагональная матрица размерности [3 х 3] вида

'1 0 0^

0 1 0

с \0 0 1

Пусть [ А] есть матрица размерности [3 Ь х 3Ь]:

Л А11] ••• [ Аш ]-

[Сц] = / П(а)и ( 0 1 0 ) .

\ П П 1 /

[А] =

1[Ам 1]

где [А^] — матрица размерности [3 х 3], такая, что для любого постоянного вектора М справедливо соотношение:

[ Агз]М = I а)ёу-аУ У ¡д (х)(М, УааV*

с [в

Пусть В = (В1,В2,...,Вь)Т, М(В)=(1юМ(В1),1юМ(В2),...,»0М(ВЬ))Т ,

Н5 = ( ^^ к(а)Нь(а)^-а,..., ^о ! ¡ь(а)Нь(а)ЛУа

с с

Учитывая (2), систему (8) можно записать в виде:

[ С ]В = Н5 + ([А] + [С ])М(В).

5. Итерационные методы решения нелинейных задач

Основная трудность реализации кусочно-постоянной (5) и кусочно-линейной

(8) дискретизаций связана с вычислением большого числа шестикратных, в общем случае сингулярных интегралов. Для упрощения объёмные интегралы (6),

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

Вк+1 = н3 + ([А] + [С])М(Вк), Во = 0, к = 1, 2,...

Итерационный процесс заканчивался, когда невязка уравнений (3) становилась меньше наперёд заданного числа . Для системы уравнений (5) использовался итерационный процесс:

[ С]Вк+1 = Н3 + ([А] + [С])М(Вк), Во = 0, к = 1,2,.

Обращение матрицы [ С] не вызывало проблем, поскольку [ С] — диагональная матрица. Итерационный процесс (12) использовался также для решения системы

уравнений (8). Однако в этом случае матрица [С] не являлась диагональной матрицей. Для решения линейной системы уравнений [С]х = у использовался метод неполного разложения Холецкого в сочетании с методом сопряжённых градиентов [4].

6. Результаты моделирования

Рассмотренные подходы дискретизации объёмных интегральных уравнений использовались для моделирования дипольного магнита эксперимента СВМ (С81, Дармштадт).

На рис. 1 приведено разбиение магнита на тетраэдры. В расчётах учитывалась дипольная симметрия поля, которая позволила сократить в 8 раз число неизвестных. Восьмая часть магнита разбивалась на 3959 тетраэдров, общее число вершин которых равнялось 1112.

Рис. 1. Разбиение магнита на тетраэдры

На рис. 2 приведено распределение компоненты поля Ву вдоль оси магнита. Результаты расчётов, использующие метод коллокаций (3), обозначены кривой 2. Результаты, базирующиеся на дискретизациях (5) и (8), соответствуют кривым 1 и 3.

Ж о 1

*// *// "У У е'УУ °УУ ""УУ 7/ -'УУ .V .V V V <ь° 5У> ¡ь> С> <£> ъ- й- ^ 4 4

Рис. 2. Распределение Ву компоненты поля вдоль оси магнита

Литература

1. Armstrong A. G. et al., 1976. — GFUN3D User Guide. — RL-76-029/A.

2. Акишин П. Г. Метод интегральных уравнений в задачах магнитостатики: Автореф. дис. канд. физ.-мат. наук. — Дубна: ОИЯИ, 1983. — 11-83-558. [Akishin P. G. The Integral Equation Method in Magnetostatic Problems: Abstract of a PhD Thesis. — Dubna: JINR, 1983. — (in russian). ]

3. Pasciak J. E. An Iterative Algorithm for the Volume Integral Method for Mag-netostatics Problems // Comp. Maths. with Appls. — 1982. — Vol. 8, No 8. — Pp. 283-290.

4. Meijerink J. A., van der Vorst H. A. An Iterative Solution Method for Linear Systems, of Which the Coefficient Matrix is a Symmetric M-matrix // Math. Com-put. — 1977. — Vol. 31, No 137. — Pp. 148-162.

UDC 519.64

The Volume Integral Equations Method in Magnetostatics

Problems

P. G. Akishin, A. A. Sapozhnikov

Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow region, Russia, 141980

In this article, use of volume integral equations method for calculations of magnetic systems is considered. GFUN program package based on integral approach to magnetostatics applies the method of collocations and piecewise constant approximations of unknown variables within the elements to discretization of equations. The limitation of this approach is related to singularity of the integral equations kernel.

An alternative to collocation method, replacing observation point by integration over discretization elements, is considered in this article. This approach enables one to use higher order approximations for unknown variables. In the context of the finite element method, the piecewise constant and piecewise linear approximations of unknown variables are considered. The problem of computing matrix elements for discretized systems of equations can be reduced to evaluation of sixth-order integrals, singular ones in the general case, over two different elements of the computational region. Possible methods are proposed for calculating this kind of integrals. Iterative processes for solving the arising nonlinear systems of dis-cretized equation are discussed. The proposed approach enables one to build discretizations with higher precision of approximation for the initial volume integral equations of magneto-statics. The proposed method was used for 3D modeling of a dipole magnet.

Comparison of results obtained for simulation of the dipole magnet using different versions of integral magnetostatics problem discretization are given.

Key words and phrases: magnetostatics, volume integral equations method, the finite element method, discretization, nonlinear problem, iterative process.