Например, Бобцов

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

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

УДК 519.632.4
ИСПОЛЬЗОВАНИЕ БЫСТРОГО МЕТОДА ГРАНИЧНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ ЗАДАЧ МАГНИТОСТАТИКИ
И.М. Ступаков, М.Э. Рояк

Рассматривается схема решения задач магнитостатики методом граничных элементов и его ускорение, основанное на быстром методе мультиполей. Такой подход позволяет значительно снизить вычислительные затраты и решать задачи большей размерности. Благодаря использованию симметричной галеркинской постановки метод можно использовать совместно с методом конечных элементов для решения нелинейных задач магнитостатики. Приводятся результаты вычислительных экспериментов, демонстрирующие эффективность рассматриваемого подхода. Ключевые слова: метод граничных элементов, метод конечных элементов, магнитостатика, неполный скалярный магнитный потенциал, сферические гармоники.

Введение

Задачи магнитостатики часто возникают при моделировании стационарных магнитных процессов. К ним относится расчет магнитных полей, возбуждаемых постоянным током или постоянными магнитами. Основным подходом для численного решения таких задач в настоящее время является метод конечных элементов (МКЭ) [1]. Для его применения область, в которой решается задача, должна быть разбита на элементы, в роли которых обычно выступают многогранники. Результат решения задачи сильно зависит от качества построения этого разбиения. По этой причине в ситуациях, когда необходимо решать задачи с большим количеством трехмерных объектов (или с объектами сложной формы), построение достаточно качественной конечно-элементной сетки может стать нетривиальной проблемой. В таких случаях может быть целесообразным использовать метод граничных элементов (МГЭ), в котором требуется сетка только на границе между объектами, что значительно упрощает ее построение.
Основными недостатками МГЭ по сравнению с МКЭ является невозможность эффективного учета нелинейных свойств объектов и плотная структура системы линейных алгебраических уравнений (СЛАУ), получаемой в результате применения метода. Первый недостаток можно обойти путем совместного использования МКЭ и МГЭ. Для устранения второго недостатка было предложено несколько различных схем быстрого МГЭ [2]: крестовая аппроксимация матриц, использование вейвлетов в качестве базисных функций, использование быстрого метода мультиполей. Далее рассматривается именно последний вариант, основанный на разложении фундаментального решения в ряд по сферическим функциям и иерархическом разбиении пространства. Впервые такой подход был предложен для быстрого моделирования систем частиц в [3], дальнейшее развитие метода изложено в [4]. Благодаря использованию этого метода можно значительно снизить вычислительные затраты в МГЭ.

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

Задача магнитостатики может быть описана системой уравнений

rot H  J, div B  0,

(1)

70 Научно-технический вестник информационных технологий, механики и оптики,
2012, № 5 (81)

И.М. Ступаков, М.Э. Рояк

где H ‒ вектор напряженности магнитного поля; J ‒ вектор сторонних токов; B  H  M ‒ вектор ин-

дукции магнитного поля,  ‒ магнитная проницаемость, а M ‒ намагниченность. Для ее решения ис-

пользуем подход с неполным скалярным магнитным потенциалом [1, 5]. Обозначим H0 ‒ решение задачи в однородном пространстве, т.е. удовлетворяющее системе уравнений

droivt HH00

 

J, 0.

Вычитая первое уравнение системы (2) из первого уравнения системы (1), получаем

(2)

rot H  H0   0 ,
значит, H можно представить в виде H  H0  grad u . Подставляя (3) во второе уравнение системы (1), с учетом связи между H и B получаем

(3)

div  grad u   div H0  M .

(4)

Для решения задачи удобно разбить область на несколько непересекающихся подобластей так,

чтобы всем границам между материалами соответствовали границы разбиения. Из уравнения (4) можно

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

u  0,

   



ij
u n

 

ij

 H0

 M  nij

,

(5)

где ij ‒ граница между i-ой и j-ой подобластями.

Если внутри подобласти магнитную проницаемость и намагниченность можно считать константа-

ми, уравнение (4) вырождается в уравнение Лапласа, и для построения аппроксимации решения в этой

подобласти можно использовать МГЭ. Такие подобласти почти всегда присутствуют в задаче, например,

подобласть, содержащая воздушную среду. В остальных подобластях можно использовать МКЭ.

Метод граничных элементов

Рассмотрим применение метода граничных элементов к краевой задаче для уравнения Лапласа u  0 . Пусть  ‒ область трехмерного пространства с границей S . Фундаментальное решение урав-
нения Лапласа определяется формулой

U

r,r



1 4

1 r r

.

Если подставить фундаментальное решение в формулу Грина

 vu


 uvdV




S

v

u n

u

v n

dS

,

в качестве функции v можно получить интегральное представление гармонической (т.е. удовлетворяю-

щей уравнению Лапласа) функции

ur



U
S

r,r

u n

r dSr




S

U nr

r,ru r dSr

,

(6)

где первый интеграл называется потенциалом простого слоя, а второй – потенциалом двойного слоя. При

получении уравнения (6) используется тот факт, что u  0 и U r,r  r  r .

Функция, удовлетворяющая уравнению (6), автоматически удовлетворяет уравнению Лапласа, значит, для решения задачи остается найти значения функции u и ее потоки на границе области так,
чтобы выполнялись уравнения (5). Для этого в МГЭ используют граничные интегральные уравнения, которые связывают значение и потоки на границе.
Для получения первого граничного уравнения точку, в которой вычисляется значение функции u
в формуле (6), устремим к границе S области  и возьмем предел. При этом необходимо учесть, что
потенциал двойного слоя терпит разрыв на границе. Запишем результат в операторной форме [2]:

u



V

u n



 

1 2

I



K

 

u

,

(7)

где потенциалы простого и двойного слоя обозначены соответственно символами V и K , тождествен-

ный оператор ‒ символом I .

Научно-технический вестник информационных технологий, механики и оптики, 2012, № 5 (81)

71

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

Для получения второго уравнения перед вычислением предела возьмем от выражения (6) произ-

водную по направлению нормали к границе S. Учитывая разрыв производной потенциала простого слоя

на границе, можно получить второе граничное уравнение

u n



 

1 2

I



K



 

u n



Du

,

(8)

где символами K  и D обозначены сопряженный потенциал двойного слоя и гиперсингулярный инте-

гральный граничный оператор соответственно [2].

Если выразить из уравнения (7) потоки и подставить в уравнение (8), можно получить оператор

Стеклова‒Пуанкаре в симметричной форме

u n



 

D



 

1 2

I



K V

1

 

1 2

I



K

 

 

u

,

на основе которого наиболее удобно производить решение задачи, разбитой на несколько подобластей.

На практике вместо обращения оператора V решают систему уравнений (7) и (8), предполагая, что по-

токи в левой части уравнения (8) даны. Для ее решения будем использовать метод Галеркина.

Запишем систему уравнений (7) и (8) в слабой форме:

Vw,


v





  

 

1 2

I



K

 

u,

v

  



0,

v



H

1/

2



S



,

 

w,

 

1 2

I



K

 

v

  





Du,

v





 

u n

,

v

 

,

v



H

1/

2



S



,

(9)

где скалярное произведение u, v  uvdS , а  H 1/ 2 S и H 1/ 2  S  ‒ пространства Соболева– S
Слободецкого [2]. Для построения дискретного аналога системы (9) представим функции w и u в виде

разложения по финитным полиномиальным базисным функциям

wr  pii r , u r  qii r . ii
Подставляя это разложение в систему (9) и выбирая в качестве пробных функций v , функции

i для первого уравнения и i для второго, получаем СЛАУ

V





1 2

M



K

T





1 2

M D



K



 

  

p q

  



  

0 F

  

,

(10)

где

   Vij  V  j , i  i rU r,r  j r dSrdSr , SS

(11)

   Kij 

K  j , i



S

S

i

r

U nr

r, r 

j

r dSrdSr

,

(12)

   Dij  D j , i  U r,r curl i r  curl  j r dSrdSr , SS

(13)

  Mij   j , i  i r  j r dSr , S
а F – вектор правой части, определяемый условиями сопряжения на границе подобласти (или краевыми

условиями). Система (10) удобна тем, что для стыковки метода на границе с МКЭ (или другой подобла-

стью граничных элементов) достаточно выбрать одинаковые базисные функции i и подставить по-

лученные дискретизации в уравнение неразрывности потоков.

Структура СЛАУ (10) является плотной, а значит, вычислительные затраты на ее сборку растут

квадратично с ростом числа элементов. Чтобы этого избежать, используем быстрый мультипольный ме-

тод. Он основан на разложении фундаментального решения уравнения Лапласа в ряд

     1
r r



 n0

rn r n1

n
Ynm
mn

r

Yn m

r

,

(14)

где Ynm ‒ сферические гармоники порядка n [3]. Подставляя это разложение в формулу (11), можно избавиться от вычисления двойного интеграла и считать интегралы по r и по r независимо:

   
       Vij 

n

 

r n Ynm

n0 mn S

r

j

r

dSr



 



S

Ynm r

r
n 1

i

r

 dSr  .

72 Научно-технический вестник информационных технологий, механики и оптики,
2012, № 5 (81)

И.М. Ступаков, М.Э. Рояк
Матрицы (12) и (13) можно вычислять аналогично. Разумеется, необходимо учитывать область сходимости ряда (14), который сходится только при r  r , что приводит к необходимости древовидного разбиения пространства [3], в котором вклады от пар элементов, расположенных в одном листе дерева, считаются напрямую по формулам (11)–(13), а в остальных случаях используется вычисление вкладов через суммирование ряда (14). Это позволяет вычислять интегралы для каждой базисной функции только один раз и, следовательно, избежать квадратичного роста затрат.
Результаты вычислительных экспериментов Для оценки точности разрабатываемого подхода сравним результаты вычисления искажений магнитного поля вокруг металлической балки, изображенной на рис. 1, при постоянном внешнем поле с результатами, полученными с помощью МКЭ. На рис. 2 приводятся графики напряженности магнитного поля, нормированные на модуль внешнего поля. Как видно из графиков, решение, полученное МГЭ на грубой сетке, точнее, чем решение, полученное МКЭ на удвоенной сетке. Отметим, что если брать сферические гармоники до 10-го порядка, то решения, полученные быстрым и обычным МГЭ, практически не отличаются. Для тестирования эффективности предложенной схемы решим задачу о вычислении искажений магнитного поля земли металлическими объектами в стенах помещения. Расположение металлических объектов представлено на рис. 3.
Рис. 1. Металлическая балка (c габаритами 0,2×0,2×1,0 м) и гранично-элементная сетка на ее поверхности
В модели задано 319 металлических прутьев с квадратным сечением 2×2 см, магнитная проницаемость металла принята равной 1000.
При решении этой задачи даже на очень грубой сетке обычным МГЭ затраты времени и памяти становятся значительно больше, чем при использовании ускоренного подхода. Так, прямой метод требует только на сборку матриц более часа времени и более 4 ГБ памяти. Использование МКЭ для решения этой задачи также затруднительно, поскольку необходимо построение сетки в воздухе между объектами. Применение же быстрого МГЭ позволяет решить эту задачу даже на относительно подробной сетке за приемлемое время. Как видно из таблицы, для быстрого МГЭ рост затрат времен и памяти при дроблении сетки близок к линейному, тогда как у обычного МГЭ он является квадратичным.

Рис. 2. Нормированные модули напряженности магнитного поля, полученные из расчетов МКЭ и МГЭ

Научно-технический вестник информационных технологий, механики и оптики, 2012, № 5 (81)

73

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

Рис. 3. Металлические объекты в стенах помещения, общая длина – 15 м; ширина – 9 м; высота – 6 м

Количество степеней свободы
22364 41156
77768
151676

Время сборки матриц, с
561 845 1435 2730

Полное время решения, с 1338
2007
3023
5944

Память, МБ
2150 3014 5433 9699

Таблица. Вычислительные затраты быстрого метода граничных элементов
Заключение
Полученные результаты подтверждают, что метод граничных элементов с использованием для ускорения быстрого мультипольного метода эффективен для решения линейных задач магнитостатики. Он позволяет даже на достаточно грубых сетках получать решение с высокой точностью и при этом существенно упрощает процедуры построения сетки по сравнению с методом конечных элементов. При этом он остается достаточно перспективным для решения и нелинейных задач совместно с методом конечных элементов, поскольку в этом случае процедура генерации сетки усложнится несущественно – сетки во всех ферромагнитных объектах могут строиться независимо друг от друга.
Литература
1. Соловейчик Ю.Г., Рояк М.Э., Персова М.Г. Метод конечных элементов для решения скалярных и векторных задач: Учеб. пособие. – Новосибирск: Изд-во НГТУ, 2007. – 896 с.
2. Steinbach O. Numerical approximation methods for elliptic boundary value problems. – New York: Springer Science, 2008. – 386 p.
3. Greengard L., Rokhlin V. A fast algorithm for particle simulations // J. Comput. Phys. – 1987. – V. 73. – P. 325–348.
4. Cheng H., Greengard L., Rokhlin V. A Fast Adaptive Multipole Algorithm in Three Dimensions // J. Comput. Phys. – 1999. – V. 155. – P. 468–498.
5. Ступаков И.М., Корсун М.М., Рояк М.Э. Об учете источников электромагнитного поля в совместном методе конечных и граничных элементов // Научно-технический вестник СПбГУ ИТМО. – 2010. – № 5 (69). – С. 67–71.

Ступаков Илья Михайлшович Рояк Михаил Эммануилович

– Новосибирский государственный технический университет, аспирант, istupakov@gmail.com
– Новосибирский государственный технический университет, доктор технических наук, доцент, профессор, royak@fpm.ami.nstu.ru

74 Научно-технический вестник информационных технологий, механики и оптики,
2012, № 5 (81)