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

Литература

1. SaweyrB.A. Linear magnetic drive system. - U.S. Patent 3,735,231, May 22, 1973.

2. Quaid E., Hollis Ralph L. 3-DOF Closed-loop control for planar linear motors // Proceeding of the 1998 IEEE International Conference on Robotics & Automation. - Leuven, Belgium. - May, 1998. - P. 24882493.

3. Butler Zack J., Rizzi Alfred A., Hollis Ralph L. Integrated Precision 3-DOF Position Sensor for Planar Linear Motors // Proceeding of the 1998 IEEE International Conference on Robotics & Automation. - Leuven, Belgium. - May, 1998. - P. 2652-2658.

4. Fries Gregory A., Rizzi Alfred A., Hollis Ralph L. Fluorescent Dye Based Optical Position Sensing for Planar Linear Motors // Proceedings of the 1999 IEEE International Conference on Robotics & Automation. -Detroit, Michigan. - May, 1999. - P. 1614-1619.

5. Miller G.L. Capacitively incremental position measurement and motion control. - U.S. Patent 4,893,071, January 09, 1990.

6. Мухаметгалеев Т.Х. Разработка замкнутого по положению планарного дискретного электропривода. Кандидатская диссертация. - М.: МЭИ (ТУ), 1994. - 171 с.

Балковой Александр Петрович - Национальный исследовательский университет «Московский энерге-

тический институт», кандидат технических наук, ст. научный сотрудник, [email protected]

Тяпкин Михаил Геннадьевич - Национальный исследовательский университет «Московский энерге-

тический институт», аспирант, [email protected]

УДК 519.632.4

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

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

Введение

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

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

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

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

(1)

div B = 0,

где Н - вектор напряженности магнитного поля; J - вектор сторонних токов; В = цН + М - вектор индукции магнитного поля, ц - магнитная проницаемость, а М - намагниченность. Для ее решения используем подход с неполным скалярным магнитным потенциалом [1, 5]. Обозначим Н0 - решение задачи в однородном пространстве, т.е. удовлетворяющее системе уравнений

J rot H0 = J, [div H0 = 0.

(2)

Вычитая первое уравнение системы (2) из первого уравнения системы (1), получаем rot (H - H0 ) = 0, значит, H можно представить в виде

H = H0 - gradи . (3)

Подставляя (3) во второе уравнение системы (1), с учетом связи между H и B получаем

div(gradu) = div(|H0 + M) . (4)

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

J[u ]г. = а

(5)

= К + M)• n]г ,

ди

I— дп

где Г.. - граница между I-ой и.-ой подобластями.

Если внутри подобласти магнитную проницаемость и намагниченность можно считать константами, уравнение (4) вырождается в уравнение Лапласа, и для построения аппроксимации решения в этой подобласти можно использовать МГЭ. Такие подобласти почти всегда присутствуют в задаче, например, подобласть, содержащая воздушную среду. В остальных подобластях можно использовать МКЭ.

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

Рассмотрим применение метода граничных элементов к краевой задаче для уравнения Лапласа Дм = 0 . Пусть О - область трехмерного пространства с границей £ . Фундаментальное решение уравнения Лапласа определяется формулой

и (г, г')=— у ; 4л |г - г'|

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

[ уДм - мДуйУ = [ V — - м — Ж, сП сП

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

м (г) = | и (г, г')^ (г')^-.[|и (г, г> (г')£, , (6)

х сП X С®Г'

где первый интеграл называется потенциалом простого слоя, а второй - потенциалом двойного слоя. При получении уравнения (6) используется тот факт, что Дм = 0 и Д и (г, г') = 8 (г - г').

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

Для получения первого граничного уравнения точку, в которой вычисляется значение функции м в формуле (6), устремим к границе X области О и возьмем предел. При этом необходимо учесть, что потенциал двойного слоя терпит разрыв на границе. Запишем результат в операторной форме [2]:

м = У ^ + (11 - К1 м , (7)

ап ^ 2 )

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

Для получения второго уравнения перед вычислением предела возьмем от выражения (6) производную по направлению нормали к границе & Учитывая разрыв производной потенциала простого слоя на границе, можно получить второе граничное уравнение

ды (1 Л ды — = 1 — I + K' I— + Du dn I 2 ) dn

(8)

где символами К' и Б обозначены сопряженный потенциал двойного слоя и гиперсингулярный интегральный граничный оператор соответственно [2].

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

| = ( Б + (2 I + К'], -(1 I + К

на основе которого наиболее удобно производить решение задачи, разбитой на несколько подобластей. На практике вместо обращения оператора V решают систему уравнений (7) и (8), предполагая, что потоки в левой части уравнения (8) даны. Для ее решения будем использовать метод Галеркина. Запишем систему уравнений (7) и (8) в слабой форме:

(Ум>,V) -21 + К]и,V] = 0,Уу е Н-1/2 (5),

'II + К] V^ + (Би, V) = ^, V], Уу е Н1/2 (5)

где скалярное произведение (и,V) = |тй5, а Н 1/2 (5) и Н1/2 (5) - пространства Соболева-

Слободецкого [2]. Для построения дискретного аналога системы (9) представим функции V и и в виде разложения по финитным полиномиальным базисным функциям

™ (г ) = X p - Ф, (г ), и (г ) = X q■ ^, (г ).

/ -

Подставляя это разложение в систему (9) и выбирая в качестве пробных функций V, функции {ф,.} для первого уравнения и {у,.} для второго, получаем СЛАУ

w,

(9)

V -((M + K) D

где

(10)

(11)

(12) (13)

(( M + K)

Vj = (Фj,Ф,) = Цф,- (r')U(r,r>, (r)dSrdSr, ,

S S

Kj = ((,Ф,) = Цф, (r')dU(r,r')vyj (r)dSrdSr' ,

S S dnr

Dj = (Dvj, y, ) = JJ U (r, r ') curl y, (r ')• curl yj (r) dSrdS, ,

S S

M =(v j, Ф,)^Ф,( r )v j (r )dSI,

S

а F - вектор правой части, определяемый условиями сопряжения на границе подобласти (или краевыми условиями). Система (10) удобна тем, что для стыковки метода на границе с МКЭ (или другой подобластью граничных элементов) достаточно выбрать одинаковые базисные функции {vy,} и подставить полученные дискретизации в уравнение неразрывности потоков.

Структура СЛАУ (10) является плотной, а значит, вычислительные затраты на ее сборку растут квадратично с ростом числа элементов. Чтобы этого избежать, используем быстрый мультипольный метод. Он основан на разложении фундаментального решения уравнения Лапласа в ряд

1

^X^ X Y*" (rKm (r').

(14)

|г - Г | п=о |Г ' | .....

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

=х XX (щт- (г)ф, (-к ^ ф,- (г ' №

п=0 щ=-п ^ 5 у! 5 г

n

XI

Матрицы (12) и (13) можно вычислять аналогично. Разумеется, необходимо учитывать область сходимости ряда (14), который сходится только при |г| < |г'|, что приводит к необходимости древовидного разбиения пространства [3], в котором вклады от пар элементов, расположенных в одном листе дерева, считаются напрямую по формулам (11)—(13), а в остальных случаях используется вычисление вкладов через суммирование ряда (14). Это позволяет вычислять интегралы для каждой базисной функции только один раз и, следовательно, избежать квадратичного роста затрат.

Результаты вычислительных экспериментов

Для оценки точности разрабатываемого подхода сравним результаты вычисления искажений магнитного поля вокруг металлической балки, изображенной на рис. 1, при постоянном внешнем поле с результатами, полученными с помощью МКЭ. На рис. 2 приводятся графики напряженности магнитного поля, нормированные на модуль внешнего поля. Как видно из графиков, решение, полученное МГЭ на грубой сетке, точнее, чем решение, полученное МКЭ на удвоенной сетке.

Отметим, что если брать сферические гармоники до 10-го порядка, то решения, полученные быстрым и обычным МГЭ, практически не отличаются.

Для тестирования эффективности предложенной схемы решим задачу о вычислении искажений магнитного поля земли металлическими объектами в стенах помещения. Расположение металлических объектов представлено на рис. 3.

Рис. 1. Металлическая балка (с габаритами 0,2*0,2*1,0 м) и гранично-элементная сетка на ее поверхности

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

При решении этой задачи даже на очень грубой сетке обычным МГЭ затраты времени и памяти становятся значительно больше, чем при использовании ускоренного подхода. Так, прямой метод требует только на сборку матриц более часа времени и более 4 ГБ памяти. Использование МКЭ для решения этой задачи также затруднительно, поскольку необходимо построение сетки в воздухе между объектами. Применение же быстрого МГЭ позволяет решить эту задачу даже на относительно подробной сетке за приемлемое время. Как видно из таблицы, для быстрого МГЭ рост затрат времен и памяти при дроблении сетки близок к линейному, тогда как у обычного МГЭ он является квадратичным.

%Н

0,88 —— —— 0,87 0,86 0,85 0,84 0,83 0,82 0,81 0,8 0,79

О 0.1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 I МКЭ, начальная сегка — — МКЭ, удвоенная сетка МГЭ, начальная сетка

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

\ ■

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

Количество степеней свободы Время сборки матриц, с Полное время решения, с Память, МБ

22364 561 1338 2150

41156 845 2007 3014

77768 1435 3023 5433

151676 2730 5944 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.

Ступаков Илья Михайлшович - Новосибирский государственный технический университет, аспирант,

[email protected]

Рояк Михаил Эммануилович - Новосибирский государственный технический университет, доктор

технических наук, доцент, профессор, [email protected]