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

Научный вестник НГТУ. - 2010. - № 2(39)

УДК 517.946:518.12:538.3:538.5

Об использовании граничных элементов при моделировании электромагнитных процессов с существенным влиянием вихревых токов*

М.М. КОРСУН, И.М. СТУПАКОВ, М.Э. РОЯК

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

Ключевые слова: метод конечных и граничных элементов, электромагнитные процессы ВВЕДЕНИЕ

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

Для моделирования нестационарных задач электромагнетизма в некоторых случаях можно использовать метод граничных элементов [5]. К преимуществам этого метода можно отнести значительное упрощение алгоритмов построения сеток, поскольку нет необходимости построения сетки внутри области. Кроме того, метод граничных элементов естественным образом учитывает неограниченные подобласти.

Однако применение метода граничных элементов для моделирования нестационарных электромагнитных процессов в чистом виде затрудняет необходимость использования векторного потенциала в проводящих областях и возможная зависимость коэффициентов уравнения от искомого поля. Кроме того, метод граничных элементов дает плотную СЛАУ.

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

* Статья получена 24 февраля 2010 г.

с помощью векторного метода конечных элементов, а в непроводящих - с помощью метода граничных элементов.

Эффективность численного моделирования на основе вычислительной схемы с совместным использованием векторного и скалярного методов конечных элементов [1] и вычислительной схемы с совместным использованием векторного метода конечных элементов и метода граничных элементов будет сопоставляться на примере решения задачи в расчётной области, которая представляет собой модель дипольного магнита с расположенным внутри магнитным экраном. Этот магнит является частью схемы экстракции пучка заряженных частиц из циклического ускорителя, разрабатываемой в ИЯФ СО РАН под руководством Н.А. Винокурова [4]. На рис. 1 изображена схема конструкции магнита и экрана в проекциях на плоскость ZY (слева) и на плоскость XY (справа), размеры указаны в миллиметрах.

Материал экрана Ш Электротехническая сталь Токовая обмотка

Рис. 1. Схема конструкции магнита: проекции в плоскостях 11 (слева) и XI (справа)

0,01

0,005

0

0,005

0,01

0,015 X

Рис. 2. Половина магнитного экрана, сечение в плоскости 1 = 0

Y

0

Конструкция симметрична относительно плоскостей Z = 0 и X = 0, поэтому в расчётную область можно включить только часть пространства Y > 0 , Z > 0, а симметрию обеспечить соответствующими краевыми условиями на плоскостях Z = 0 и X = 0.

Модель магнитного экрана показана на рис. 2. Экран содержит по два кольца меди и стали; ближайшее к центру кольцо экрана стальное. Толщина кольца меди составляет 0.1 мм, стали 0.08 мм, толщина воздушного зазора 0.11 мм. Относительная магнитная проницаемость стали 547, коэффициент проводимости меди 57 • 106 См/м, стали 3.3 • 106 См/м. Относительная магнитная проницаемость электротехнической стали 1000. Ток в обмотках изменяется по следующему закону:

I = 5580 • sin_4 j А, t е[0,500] мкс.

1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Система уравнений Максвелла, используемая для описания нестационарных электромагнитных процессов, при несущественности токов смещения принимает вид

rot Й = Jст +сЁ, (1)

Л 8B

rot E =--, (2)

8t

div B = 0, (3)

где H - напряженность магнитного поля; JCT - вектор плотностей сторонних токов; E -напряженность электрического поля; с - удельная проводимость среды; B - индукция магнитного поля (B = цН , где ц - коэффициент магнитной проницаемости, который в общем

случае может зависеть от поля H).

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

(возможно, неодносвязных) подобластей Qc и Q0, причём в подобласти Q0 удельная проводимость с = 0, отсутствуют сторонние токи и нельзя построить контур, по которому можно было бы обойти ненулевой ток [1].

Тогда в подобласти Qo напряженность магнитного поля H0 можно представить в виде

H0 =- gradU, (4)

где U - полный скалярный магнитный потенциал. Соответственно индукция магнитного поля B0 в подобласти Q0 представляется в виде

B0 = -ц gradU. (5)

Следовательно, для выполнения уравнения (3) необходимо, чтобы скалярный потенциал U в области Q0 удовлетворял уравнению

- div (ц gradU ) = 0. (6)

В подобласти Qc индукцию магнитного поля Bc представим через вектор-потенциал A с помощью соотношения

Bc= rot A . (7)

При таком представлении B в Qc автоматически выполняются уравнения (2) и (3) из системы уравнений Максвелла. Уравнение (1) с учетом соотношения (7) преобразуется к виду

rot

(1

—rot A

ц

+ с— = JCT (8)

8t

VI

Таким образом, в подобласти с нулевой проводимостью индукция магнитного поля Б0 определяется через скалярный потенциал и и магнитное поле описывается уравнением (6), а в подобласти индукция Б° определяется через вектор-потенциал А и электромагнитное поле описывается уравнением (8). Чтобы индукция Б, определяемая в подобласти соотношением (5), а в подобласти соотношением (7), удовлетворяла системе (1)-(3) во всей расчетной области О = и , на границе S между подобластями и должны выполняться условия непрерывности нормальных составляющих Б и тангенциальных составляющих Н :

Б0 • п = п , (9)

Н0 х п = Н°х п, (10)

где п - любая нормаль к рассматриваемой границе S.

2. СОВМЕСТНОЕ ИСПОЛЬЗОВАНИЕ ВЕКТОРНОГО МЕТОДА

КОНЕЧНЫХ ЭЛЕМЕНТОВ И МЕТОДА ГРАНИЧНЫХ ЭЛЕМЕНТОВ

Аппроксимация по времени уравнения (8) приводит к векторному уравнению

(

rot

1

—rot A Ц

Л

+ yA = F :

(11)

где коэффициент у и вектор-функция F определяются разностной схемой аппроксимации по времени.

Пусть П - внешняя нормаль по отношению к , тогда с учетом соотношений (5) и (7) условия непрерывности (9)-(10) принимают вид

Г 1 X ^

—rot A х n

VЦ

-Ц

и

= -( gradU х n )|(

:(rot A • n )| ^ .

(12)

(13)

Будем считать, что на границах расчетной области О = О0 и в качестве краевых условий могут быть заданы краевые условия равенства нулю касательных либо нормальных составляющих магнитной индукции. Обозначим через ^х участок границы, на котором задано условие равенства нулю касательных составляющих В, а через Sn - участок границы расчетной области О, на котором задано условие равенства нулю нормальных составляющих В .

Для скалярных произведений введем следующие обозначения, соответствующие интегралам по объему и границе: (V, w) = | V • wdV, (V, w) =|V • wdГ . С учетом введенных

V Г

обозначений эквивалентная вариационная постановка для векторного уравнения (11) и условия сопряжения (12) принимает вид [1]

f

—rot A,rot Ф Ц

Г rot ,

л

+ (gradU х n, Ф)^+(уА, Ф) = (F, Ф) , УФ е Hr0ot (Q0). (14)

Под H0ot (Qa) понимается пространство вектор-функций Ф, определённых на Qa, для которых функция rot Ф является суммируемой с квадратом, при этом касательные всех функций Ф должны быть равны нулю на границе Sn n Sа (Sа - граница области Qa).

Если коэффициент магнитной проницаемости в области Q0 постоянный, уравнение (6) эквивалентно уравнению Лапласа

- div (grad U ) = 0. (15)

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

U (N )=i J^lr(M) dSM - i iU (M У

(

1

Л

V rnm J

dS,

M

(16)

где RNM - расстояние между точками N и М , Пм - нормаль в точке М . Предполагается, что точка N находится внутри области О0, точка М , по которой ведется интегрирование -

s

s

Q

а

на границе S. Представление (16) справедливо и для неограниченной области. Для простоты изложения будем полагать, что граница области О0 является бесконечно удаленной (за исключением участка границы с областью ).

Первый интеграл в формуле (16) называется потенциалом простого слоя. Введем связанный с ним граничный интегральный оператор V :

Vv = Нт — Г—V(М).

N—Nе 4п1 ^ ^ ' М

(17)

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

К =1V (N') + Нт — Г-

1 У ' И V АГ'^С /I тт 1

(

1

Л

V ^м )

V (м) ё8м .

(18)

Введем также сопряженный с оператором К оператор К , связанный с нормальной производной потенциала простого слоя

( Л \

К\ = -1V(N')+ lim пn,• grad — Г—V(М)dSl

V 4л 8 ^М

М

(19)

Оператор D, связанный с нормальной производной потенциала двойного слоя, введем следующим образом:

(л Л \ \

V (М)dSм

Dv = - lim п ы, • grad

N—N 'е8

- V

V4л 8ж м

1

V ^М )

(20)

Отметим, что при вычислении предела интегралы в формулах (17)-(19) становятся несобственными, а интеграл в формуле (20) понимается в смысле конечной части по Адамару [6].

Пусть I - тождественный оператор. Тогда, используя операторы (17)-(20), можно вы, ди

числить значения и\„ и потоки - как пределы функции и, представленной в виде (16),

8 дп 8

при N — N' е 8 и записать систему граничных интегральных уравнений:

8

ди

ли ( '1т

V- + 1 -I -

ей 8 2

(1 ^ ди

-I + К

V 2

+ DU

8 '

(21) (22)

Отметим, что уравнения (21) и (22) выполняются в случае гладкой границы 8 . В случае кусочно-гладкой границы их следует понимать в обобщенном (слабом) смысле.

При решении задач с использованием декомпозиции области удобно использовать эти уравнения не напрямую, а предварительно определив оператор Стеклова-Пуанкаре, который в явном виде связывает значения и потоки на границе. Для этого выразим из уравнения (21) слагаемое с потоком

V

ди

= 21 + К |и15-

(23)

8

8

Оператор V является обратимым [6], поэтому можно выразить из этого уравнения поток и подставить его в правую часть уравнения (22):

dU

аТ

11 + К*Л V"1 i11 + K | + D

U

S '

(24)

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

Используя оператор Стеклова-Пуанкаре, можно сформулировать эквивалентное условию (13) вариационное уравнение

f

11 + K*l V"1 ( 11 + K | + D

Л

Us ,v

:( rot A • n, v) , Vv e H1 ( S ),

(25)

где Н1 (S) - гильбертово пространство скалярных функций V , определенных на границе S и

имеющих интегрируемые с квадратом первые производные.

Однако использовать уравнение (25) для построения аппроксимации мешает обратный

лт-1

оператор V , поэтому введем вспомогательное уравнение, являющееся вариационным аналогом уравнения (23):

(

V

8U

m

Л

, v

S Js

ff1 Л i ^

-I + K ,V

S

Vv e H 0 (S ),

(26)

S

где Н0 (S) - гильбертово пространство скалярных функций V , определенных на границе S

и интегрируемых с квадратом.

Подставив уравнение (22) в условие (13), получим следующее вариационное уравнение:

(

V i+K-18U

1 2 ) 8n

Л

-(-MDU|s , v) =(rot A • n, v) , Vv e H1 (S). (27)

s Js s s

Система уравнений (26) и (27) эквивалентна уравнению (25), но не содержит обратного

оператора V .

Используя сопряженность операторов К и К* , уравнение (27) можно записать в виде

(

Ц „Л 8U

21+K)v, я

s Js

+ (-^DU|S, v) = (rot A • n, v) , Vv e H1 (S).

(28)

Уравнения (26) и (28) позволяют построить матричные аппроксимации граничных интегральных операторов. В ряде случаев гораздо эффективнее собирать СЛАУ в виде, аналогичном уравнению (25), с обращением матрицы оператора V, что позволяет уменьшить размерность СЛАУ, исключив из нее неизвестные, отвечающие за потоки.

Если сначала пронумеровать неизвестные, отвечающие за аппроксимацию скалярного потенциала и, а затем неизвестные, отвечающие за аппроксимацию векторного потенциала

А, то получится блочная матрица вида

(PT2 Л

vT1 F J

(29)

где Р - матрица аппроксимации оператора Стеклова-Пуанкаре; F - конечно-элементная матрица, соответствующая объёмным интегралам в левой части уравнения (14); 7 - матрица,

S

S

соответствующая поверхностному интегралу в уравнении (14); - матрица, соответствующая поверхностному интегралу в правой части уравнения (28). Заметим, что в (29) только матрица F зависит от коэффициентов уравнения и шага по времени, что позволяет не только не пересчитывать сами матрицы Р, Т и Т2 при выполнении итераций по нелинейности и шагов по времени, но и однократно выполнять разложение (как полное, так и неполное) соответствующей части матрицы (29).

3. РЕЗУЛЬТАТЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

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

Анализ эффективности методов проведем путем сравнения вычисления компоненты Вх магнитной индукции вдоль двух линий съема, ориентированных по оси 2 от 0 до 0.11 м; начальная точка первой линии съема (0.001, -0.0145, 0) м, начальная точка второй линии (0.001, -0.0075, 0) м. Измерения по заданным линиям проведем в зазоре дипольного магнита, где возмущение магнитным экраном однородного поля наиболее существенно. Высокое возмущение поля может привести к смещению пучка на магнитный экран и вследствие этого - к потере части частиц или всего пучка при экстракции. Поэтому к точности моделирования поля в этой области предъявляются особые требования.

На рис. 3, а и б показана относительная разница вычисления компоненты Вх магнитной индукции в процентах относительно решения методом конечных элементов на подробной сетке на одном из времен моделирования (t = 50 мкс).

1

0,5 0

-0,5 -1 -1,5 -2 -2,5

Вх, % -0,5

-1

-2,5

0,02

0,08

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

г - решение, полученное совместным использованием метода конечных и граничных

элементов на подробной сетке X - решение, полученное методом конечных элементов на грубой сетке

Рис. 3. Относительная разница вычисления компоненты Вх : а — измерения вдоль первой и б — вдоль второй линий съема

В,, %

1,5

1,5

2

3

0

0,04

0,06

0,1

0,08

ь

0,06

0

0,02

0,04

0,1

ь

б

а

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

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

Характеристика вычислительных затрат Смешанный метод Векторно-скалярный МКЭ

грубая сетка подробная сетка грубая сетка подробная сетка

Время решения задачи на одном временном слое, с 2 21 15 360

Время для предварительных вычислений, мин 1.3 33 - -

Число неизвестных в СЛАУ 4522 29958 110517 833583

Размер оперативной памяти, требуемой для решения, Мб 35 483 726 12256

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

Следует отметить, что для решения СЛАУ во всех случаях использовался параллельный прямой метод решения PARDISO из библиотеки Intel MKL, оказавшийся быстрее итерационных методов. Расчеты проводились на компьютере со следующими характеристиками: процессор - Intel Xeon 8-ядерный с тактовой частотой 3.0 ГГц; размер оперативной памяти 32 Гб. Очевидно, что в условиях ограниченной памяти итерационные методы более предпочтительны.

ЗАКЛЮЧЕНИЕ

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

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

СПИСОК ЛИТЕРАТУРЫ

[1] Соловейчик Ю.Г., Рояк М.Э. Совместное использование узловых и векторных конечных элементов для расчета трехмерных нестационарных электромагнитных полей // Сиб. журн. индустр. математики. - 2004. - Т. 7. -№ 3(19). - С. 132-147.

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

[3] Тихонов А.Н., Самарский А.А. Уравнения математической физики: учеб. пособие. - 6-е изд., испр. и доп. -М.: Изд-во МГУ, 1999.

[4] Bondarenko A.V., Vinokurov N.A. Beam extraction from a synchrotron through a magnetic shield // Digest reports of the 17 Int. Synchr. Radiation Conf. - Novosibirsk, 2008. - P. 4-4.

[5] Hiptmair R. Boundary element methods for eddy current computation // Lect. Notes in Appl. and Comp. Mech. -2006. - N. 29. - P. 213-249.

[6] Langer U., Steinbach O. Coupled Finite and Boundary Element Domain Decomposition Methods // Lect. Notes in Appl. and Comp. Mech. - 2006. - N. 29. - P. 61-96.

Корсун Мария Михайловна, аспирант кафедры прикладной математики Новосибирского государственного технического университета. Основное направление научных интересов - математическое моделирование электромагнитных процессов. Имеет 11 публикаций.

E-mail: maria.korsun@gmail.com

Ступаков Илья Михайлович, аспирант кафедры прикладной математики Новосибирского государственного технического университета. Основное направление научных интересов - математическое моделирование электромагнитных процессов. Имеет 2 публикации.

E-mail: istupakov@gmail.com

Рояк Михаил Эммануилович, доктор технических наук, профессор кафедры прикладной математики Новосибирского государственного технического университета. Основное направление научных интересов -математическое моделирование электромагнитных процессов. Имеет более 50 публикаций.

E-mail: royak@fpm.ami.nstu.ru

M.M. Korsun, I.M. Stupakov, M.E. Royak

On using of boundary element for electromagnetic field with significant influence of eddy currents simulation

The calculation scheme based on coupled vector finite element and boundary element method are considered. This scheme allows simulating of electromagnetic fields with eddy current in conductor object. The efficiency of considering method is shown by solving of model problem.

Key words: boundary and finite element method, electromagnetic process.