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

УДК 621.3.011.013.

ИНТЕГРАЛЬНЫЕ УРАВНЕНИЯ ДЛЯ РАСЧЁТА ТРЁХМЕРНОГО МАГНИТНОГО ПОЛЯ В АНИЗОТРОПНОЙ СРЕДЕ

С.Н. Кадников д-р техн. наук, И.Е. Сергеева, асп.

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

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

Рассмотрим следующую модельную задачу. Магнитное поле в одно родной анизотропной среде создаётся некоторым распределением постоянных токов с плотностью 8, локализованных в области V, ограниченной поверхностью Б. Векторы индукции В и напряженности Н должны подчиняться уравнениям

гоЖ = 8; ^уВ = 0

(1) (2)

в области у и уравнениям

гоЖ = 0;

С^Б = 0 (4)

в неограниченной области Уе , внешней по отношению к V' . Выберем декартову систему координат, так чтобы оси х,у,2 были параллельны главным осям тензора абсолютной магнитной проницаемости Д а = д0 Ду, и пусть Д - диагональный тензор относительной магнитной проницаемости, а дх,ду,д2

- его диагональные компоненты (остальные равны нулю). Тогда

Б = 'До Дхнх + !доДуну + кДоД2Н2 ■ (5)

Введём векторный потенциал посредством соотношения Б = го1А . Учитывая, что Н = Да1Б , из (1) получаем уравнение относительно А

го^ a1rota = д0 8. (6)

Введя новое выражение для векторного потенциала А1 = ДА, положим С'уА1 = 0 . (То, что данное условие действительно может быть выполнено, установлено ниже). Тогда, после замены переменных х = д/ДДХХ, У = ^/ДуУ1 , 2 = •\/д721 , уравнение (6) можно записать в виде одного векторного уравнения

Пуассона

д 2А1 д2А 1 д 2А1 -

^Г + + = -Д 0ДхДУД28 (7)

Решение данного уравнения может быть записано в виде

% = Д0 Дх ДуДиТ- | 8

сч

4п ^ (8)

Переходя к исходным координатам х,у,2 и функции А, получаем для неё следующее выражение:

А ( ) Д-1 ^р

Ач = Д0 (хДуД2) 77.) 6р^, (9)

4п Я " Р

где

Ра =

(хд -хр)2 + (-Ур) , ( -2р)2^

Дх

(10)

Учитывая, что А 1= ДА, нетрудно заметить, что условие СмА1 ^ 0, с помощью которого из (6) получено (7) будет выполняться, если

СУр

" (11)

8р^- = 0

V

Если 8 = 0^а + Ь в объёме V, что следует из уравнения (1), то равенство (11) будет тождеством при любой дифференцируемой функции помещенной на место Р-1 [3]. Таким образом, формула

(9) действительно даёт решение уравнения (6), через которое могут быть выражены векторыБ и н в уравнениях (1)-(4). В частности, для вектора напряженности Н = Д-1 rot А получаем формулу

Н =—-ГадТ (12)

4п(Дх Ду Д2) V Ра

где г = !(хч -хр) + ](уч -ур) + к(ич-ир), а функция Ра определяется формулой (10). Для линейного тока силой I в замкнутом контуре I формула (12) принимает вид:

1 г гстр,г]

Нч =—-То5 У-рТ"^ (13)

4п(х Ду Д2) I Ра

Данная формула выражает закон Био-Савара-Лапласа для однородной анизотропной среды. При построении ИУ будет использоваться формула, аналогичная (9):

К / \ 0,5 Д-1 г сБр

Aq = Д0 (хДуДи) "4Л"У 'рр (14)

Б а

где ^ -плотность поверхностных токов.

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

д2 ф д2ф д2ф „

Дх + Ду ТТ + Д? ТГ = 0 дх2 ду2 д22

(15)

Фундаментальным решением данного уравнения является функция О,1, где Оа определена формулой (10). Если считать, что в некоторой области изотропной среды V расположены объемные магнитные заряды с плотностью р , то действуя также, как и при выводе (9), можно найти: - выражение для потенциала объемных зарядов в анизотропной среде:

_Гр ^

0,5 ^ Р о

Фа =

^ДхДуДи) , V

(16)

- и аналогичное выражение для потенциала простого слоя зарядов, распределенных на поверхности Б :

1 А ЬБр

Фа = 4"7-) Нт"

4п(Д хД уД 2) б а Отметим, что потенциал (16) удовлетворят уравнению (15) всюду вне Б .

Рассмотри теперь методику построения ИУ, используя следующую модельную задачу. Постоянные токи с заданной плотностью ст 0 расположены в неограниченной области Уе анизотропной среды с диагональным тензором магнитной проницаемости Д ае = д0Д е (рис. 1).

К Л,

Рис. 1.

Внутренняя ограниченная область у также заполнена анизотропной средой с диагональным тензором проницаемости Да = д0Д¡ (Д¡ ф Де). Постоянные токи локализованы в области У0 е Уе . Вектора Н,Не вторичного поля, вызванного вторичными источниками на Б , должны подчиняться уравнениям: вектор Не -уравнению (1) в области У0 и уравнению (3) в области Уе - У0 ; вектор Ц- уравнению (3) в области У|.

Вектор В во всем пространстве, исключая Б, должен удовлетворять уравнению (2).

На границе раздела сред должны выполняться условия непрерывности касательных составляющих, векторов напряженности полного поля: Н| = Н + Н0, Н'е = Не + Н0е . Вектор Н0 поля заданных токов определен только в области у при условии, что все пространство заполнено однородной анизотропной средой с тензором Да, аналогично вектор Н0е определен только в области Уе при условии, что все пространство заполнено анизотропной средой с тензором ДХае . Согласно данному представлению внешнего поля, касательные составляющие векторов вторичного поля должны подчиняться граничному условию

[П,Не-Н ] = [п,Ци - Ное ] (17)

Кроме того, на поверхности Б должны быть непрерывны нормальные составляющие индукции полного поля В', Ве . Это даст еще одно граничное условие

(п, ДеНе - Д¡Н|) = (п, Д¡Н0| - ДеНое ) (18)

Для построения ИУ в данной задаче необходимо применять метод разделения областей, в соответствии с которым в каждой из областей У',Уе следует использовать специальные представления

для векторов поля Ц,Не . В области Уе искомый вектор представим в виде Не = rot Ае, где Ае определим формулой (14). В результате

^] nme t Rae

где me = (XeHye Uze) ^ ; Rae - определено формулой (10) при цх = |xe , |y = |ye , lz = ize , X = Xe , У = y

Heq = 1-Ф1 P3e J dSp (19)

eq 4nme J Rae P

В области у вектор H определим как градиент потенциала (16):

H iq = inm bR3;dSp (20)

где mj = (xi lyiIzi) ' , Rai - определено формулой (10) при ix = ixi, ly = lyi, lz = izi, x = Xj, у = yj, z = zi •

В выражениях (19), (20) вектора Hi,He определены вне S . Чтобы использовать граничные условия (17), (18), необходимо найти предельные значения выражений [n,He ] , ^n,Hi ] , (( |eHe) , (( |iHi) на поверхности S • Будем считать, что точка q находится на нормали к поверхности S в области Ve (вне S). Считая S поверхностью Ляпунова [4], рассмотрим выражение

[n q'Hq ] = i^fe^ Sp <21)

e S ae

Используя формулу [ a | b,c]] = b(a'C)-c(a,b), где аДе - произвольные векторы, приведем (21) к

виду

|n H 1 1 X ■ (nq'Fe)dS 1 Г F (nq, 1)

|nq eq]=4nmef pi^Sp-4nmef^

dSp (22)

Второй член в данном выражении при q е Б является несобственным интегралом, что можно пока-

х

зать с помощью теории потенциала [4]. В первом члене произведем замену переменных: х1 = - е

/ix

ye ze

y1 = , z1 = ,— , в результате чего он примет вид

д/M-ye V^ze

1 Г г (n1q'R1 )JO Si 1

где индекс «1» указывает на использование координат х1,у1,21; а

( 2 2 2 0,5

Р1 = ((х1с1 - х1р) + (у1с) - у1р) + (1с) - 21р) I . Данное выражение является производной по нормали потенциала простого слоя с плотностью 1р . Его предельное значение на Б известно [4]:

1 г х (М1), Б1 1

Возвращаясь к переменным хе,уе,2е и подставляя полученные выражения в (22), получим, что

¡¡q 1s -г (n1q'R1)

= 2 + 4П $ ^HR? ^ (23)

q^S-, S1 '

на S

1

Гп H 1 q + 1 nq [ pre]]t Lnq'Heq ] = 2 + (f-R3-1

О „ T 3 dSp (24)

2 4nme S R3ae p

Предел выражения |п,й| ^ на поверхности Б , где вектор Ц представлен формулой (20), является

сингулярным интегралом, существующим в смысле главного значения [5]. В итоге, подставляя полученные выражения для [п,й], |п,йе ] в граничные условия (17), получаем первое из искомых интегральных уравнений:

* + ^^ ^ - ¥ ^Бр =2 [«Л - йое, ] (Ь)

е б ае I б а1

Для построения второго интегрального уравнения необходимо вычислить предельное значение на Б выражения (п, ДЦ), где й| определено формулой (20). Действуя так же, как и при выводе (23),

получим, что на Б

((, Дй I ) = + ^^ (26)

Предел на Б выражения (п, Дейе), где йе выражено формулой (19), является сингулярным интегралом. В результате, подставляя (п, Д|й|) и (24) в граничное условие (18), получаем второе интегральное уравнение:

ap- 2nm У ap

1 S

(ñ q'^T)

- k1

1 i- (ñqA¡ \ b'^e 1) / - - \

dSp + -¡-¿ГT R3 dSp = 2 (H "AeH)e) (27)

2nme S Rae

Уравнения (25), (27) образуют искомую систему ИУ для решения данной модельной задачи. К ядру первого интеграла добавлена константа к1, что равносильно выполнению условия

стрЬБР = 0 (28)

Б

которое необходимо для единственности решения системы (25), (27).

Чтобы доказать единственность решения (25), (27), что необходимо для ее численного решения, положим правые части равными нулю и покажем, что полученная однородная система будет иметь только нулевое решение. Отметим вначале, что краевое условие (18), из которого было получено уравнение (27), после добавления константы к1 к ядру первого интеграла может быть переписано в виде

к1Д0

(( ДaeHe " ДaiHi)" ¡ПО ^pdSp = (( Дa¡H0¡ -ДaeH0e )

¡ S

2nm|

Проинтегрируем данное равенство по S . Учитывая, что ДatHe = Be = rot Ae , Ц = -Уф|,

div fir+д°д* IF )=0

согласно (15), после применения теоремы Гаусса получаем условие (28), которому будет удовлетворять любое решение (27). Иными словами, при анализе системы (25), (27) можно рассматривать уравнения без учета константы k1, но при дополнительном условии (28).

Перепишем однородное условие (17) в следующем виде: [п, Д-1rotAe] + [п, V^] = 0 . Умножим его

скалярно на Ae и после циклической перестановки в смешанном векторном произведении и интегрирования по S получим

ф(a>tAe,Ae] п) dS + ф([v9j,Ae], П) dS = 0 (29)

S S

Применяя теорему Гаусса для первого члена, найдем

ф (|[дд-e rot A e,AA e ] n )dS = - j" (rot Дrot A e,A e ) dV + J (Дrot A e, rot A e) dV =

Ve

= До j (xe^x + AyeHy + AzeH-z) dV = j (iBeHe)dV

(30)

Здесь было использовано векторное тождество d¡v [a,bj = (brota - arotb) и уравнение rot Д-1 rot Ae = 0 . Второй член в (29) может быть преобразован следующим образом:

ф((Ae] n) dS = ф(rot (9iAe),n) dS - фф (rot Ae,n) dS = -<£ф (BjdS) = - J div (Bj) dV =

S S S S Ve

e (31)

= Д0J( + HiyHjy +^iZH2z) dV = {(Bi,Hi) dV

V V

В результате, складывая (30), (31), согласно (29), получим

J(BeHe) dV + J(BHi) dV = 0

Ve V

Подынтегральные выражения в каждом члене данного равенства неотрицательные (они выражают удвоенную плотность энергии магнитного поля) и поэтому Hi = 0,He = 0 . В области V вектор Hi определен равенством Hi = -Vci; и поэтому ci = c1 = const. При этом следует учесть, что ci является потенциалом простого слоя, согласно формуле (16). Определим его в области Ve, считая, ее заполненной магнитной средой с тензором проницаемости Дai. Переходя к координатам x1,y1,z1, можно показать, [5] что полученный потенциал (обозначим его ce ) при условии (28) будет равен нулю. Поскольку плотность простого слоя на S определяется разностью нормальных производных потенциалов ci и ce на S, приходим к выводу, что ст = 0 . В результате, вместо уравнения (25), правая часть которого равна нулю, получаем уравнение

i+dtf^-0 ,32,

Необходимо проверить, имеет ли данное уравнение ненулевые решения. Ранее было установлено, что He = rot Ae = 0 . Будем считать, что все пространство заполнено однородной средой с тензором проницаемости Дae и определим вектор Ai в области V согласно формуле (14). При условии divS i = 0 , которое следует из уравнения (32), rot Hj = 0 [6]. Кроме того, div ДeH' = 0 , и поэтому вектор H' можно определить формулой H' = -Vc'. Потенциал c' будет удовлетворять в V уравнению (15). Поскольку He = 0 , то на внутренней стороне S согласно краевому условию (18) (при H0i = Hae = 0)

Bin = 0 или (п, Д eVc') = 0 (33)

Далее воспользуемся формулой

div ф'Д eVcl = ф' div Д eVcl + (Vtf, Д eVc') (34)

Первый член справа в этой формуле, согласно (15), равен нулю. Считая функцию ф' дифференцируемой везде в V, интегрируя (34) по объему и применяя теорему Гаусса, можно показать, что ф' = ci = const и Hj = 0 . Разность касательных составляющих Hj и He на S равна i , отсюда вытекает,

что i = 0 и уравнение (32) имеет только нулевое, а система (25), (27) - единственное решение. Однако этот вывод справедлив только для односвязных областей, в которых любой замкнутый контур можно стянуть в точку. Для многосвязных областей данный вывод не имеет места.

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

любому замкнутому контуру внутри обмотки равна числу ампер-витков. Применительно к данной задаче это равносильно следующей цепочке равенств:

<jj H[di = -<jj V9¡dí = ф+ - ф[_ = фindl = I ф 0 (35)

1 ii l¡ где ^ - любой замкнутый контур, охватывающий S изнутри; li - контур такого же типа, лежащий на внутренней стороне S ; ф[+, ф'_ - значения потенциалов на противоположных сторонах вспомогательной поверхности S1, пересекающей тор (рис.2); in - нормальная к li составляющая тока i . Из соотношений (35) следует, что для того, чтобы поле внутри S было равно нулю, необходимо условие

<jj( T,ni )dl = 0 (36)

i¡

где n¡ - нормаль к li, лежащая в плоскости, касательной к S . Покажем, что данное условие достаточно для того, чтобы поле внутри S было равно нулю.

Проинтегрируем (34) по области у с исключенной поверхностью S1 и, применяя теорему Гаусса, получим

|ф;(д^ф;,П) dS + |ф+(д^ф,п) dS-[ф-(^ф,п) dS = JC^v^) dv

s s1 s1 v¡

Согласно уравнению непрерывности div Д eV^ = 0 на S1 (Д eVф¡+ ,П) = (Д ^ф^,П). Кроме того, на поверхности S нормальная составляющая вектора ДeVф¡ равна нулю. Поэтому вместо (37) получим

iJ(^¡,n) ds = J(ДeVф;) dv.

S1 v¡

Согласно условию (36), I = 0 и поэтому H¡ = 0 . Поскольку разность касательных составляющих Hi и He на S равна i , то i = 0 , т.е. уравнение (32) на двусвязной поверхности S не будет иметь ненулевых решений, что обеспечивает единственность решения системы (25), (27). Данный вывод справедлив, очевидно, при любой форме двусвязной поверхности S .

При численном решении (25), (27) условие (36) целесообразно внести в уравнение (25). Отметим, что в условии (36) контур l, охватывающий S, может быть выбран произвольно в силу условия

divS i = 0 [6], вытекающего из уравнения (25). Зафиксировав контур li, умножим (36) на k1ni, где k1 -произвольный численный коэффициент, и сложим с уравнением (25):

^ + + Miq^i) ^¿Ь iR^P = 2 [Пq,H0i _H0e] (38)

e s Rae li i s Rai

Данное уравнение и уравнение (27) образуют систему, решающую поставленную краевую задачу. Чтобы в этом убедиться, перепишем (38) с учетом (17) в виде

[а Д_¡rotÁe] + k2Lniq<j> ((ni) dl _ [ VФi] = [nq,Hüi _ H0e ]

Умножим данное равенство векторно на П , затем скалярно на орт-вектор t¡ , касательный к контуру li, и проинтегрируем по li. В результате все члены, кроме содержащего интеграл, обратятся в нуль и данное соотношение сведется к условию (36). Таким образом, сумма уравнения (25) и условия (36) равносильна уравнениям (25) и (36), т.е. система (25), (27), (36) равносильна системе (38), (27). Если численное решение (43),(32) производится путём редукции к системе линейных алгебраических уравнений (СЛАУ), то процесс вычислений можно организовать так, что добавление условия (36) к уравнения (25) будет равносильно добавлению одной или двух строк к матрице СЛАУ.

В другом варианте данной краевой задачи, когда источники внешнего поля расположены в области V, искомая система уравнений может быть построена по той же методике. В области у напряжённость вторичного поля определяется как Hi = Д-j1 rot Á¡ , где

A ¡ = _h¡L. ф A ¡ 4nm¡ T p Ra¡

В области Ve - He = ^ф,, где

1 х dSp

Фе 4nme |apRae Формулы для предельных значений [п,И|] и (n, ДaeHe) на S имеют вид

[ n^ ]-!+iL f Ч^ (39)

S

(nq,|ДaeHe) = ^ + i ар (nq,re^dSp (40)

v ^H-ae e) 2 4nme у p R3 p v у

S ae

После подстановки выражений для H0,He в граничные условия (17), (18) и использования соотношений (39),(40) получится следующая система уравнений, аналогичная по структуре системе (25), (27):

и -2nmSdSp + 2^S°p^dSp =2 [П,,*.-Hoe] (41)

i S ai es ae

а + _L_ & а iL & (П q [ M)dS =

+ 2nme PIRT^ 2nm,f r3 dSp = (42)

= 2 ( nq,Hoi - Hoe)

Данная система имеет единственное решение, если область V| - односвязная. Интегральный оператор относительно а во втором уравнении в корректировке не нуждается, поскольку интегральное уравнение внешней задачи Неймана имеет единственное решение. Если область V двухсвязная и ограничена поверхностью тороидального типа, то нужно использовать дополнительное условие вида

$(T,ne ) dl = 0 ,

le

l п

где e - замкнутый контур, лежащий на внешней стороне S; e - орт

l S

вектор нормали к e, лежащий в плоскости, касательной к S (рис. 2). Данное условие обеспечивает

H

равенство нулю циркуляции вектора e по любому контуру, охватывающему поверхность S. После K n

умножения на 1e оно должно быть добавлено к уравнению (41),что обеспечит единственность решения системы (41),(42).

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

Литература

1. Тозони О.В., Маергойз И.Д. Расчёт трёхмерных электромагнитных полей. - Киев:Техника,1974. - 352 с.

2. Маергойз И.Д. Расчёт статических полей в кусочно-однородных анизотропных средах // Изв. АНСССР. Энергетика и транспорт. - 1972. - № 2. - С. 117-125.

3. Тамм И.Е. Основы теории электричества. - М.: Наука,1976.- 616 с.

4. Гюнтер Н.М. Теория потенциала и приложение к основным задачам математической физики. -М.:Гостехиздат,1953. - 360 с.

5. Михлин С.Г. Линейные уравнения в частных производных. - М.: Высш. шк., 1977. - 431 с.

6. Кадников С.Н. Метод интегральных уравнений для расчета электромагнитного поля / Иван. гос. энерг. ун-т. - Иваново, 2003. - 340 с.