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

УДК 621.3.017.31

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

© 2014 г. И.Б. Подберезная, Ю.К. Ершов, А.В. Павленко

Подберезная Ирина Борисовна - канд. техн. наук, доцент, Podbereznaya Irina Borisovna - Candidate of Technical Sci-

кафедра «Электромеханика и электрические аппараты», ences, assistant professor, department«Electromecanics and

Южно-Российский государственный политехнический уни- Electric DevicesPlatov South-Russian State Polytechnic Uni-

верситет (НПИ) имени М.И. Платова. Тел. (86352)55-1-13. versity (NPI). Ph. (86352)55-1-13. E-mail: podbereznayaib@

E-mail: [email protected] mail.ru

Ершов Юрий Константинович - канд. техн. наук, доцент, Yershov Yury Konstantinovich - Candidate of the Technical

кафедра «Теоретическая электротехника и электрообору- Science, associate professor, department «Theoretical Electrical

дование», Южно-Российский государственный поли- Equipment and Electric Equipment», Platov South-Russian

технический университет (НПИ) имени М.И. Платова. State Polytechnical University (Novocherkassk Polytechnical

Тел. (86352)55-3-08. Institute). Ph. (86352)55-3-08.

Павленко Александр Валентинович - д-р техн. наук, про- Pavlenko Alexander Valentinovich - Doctor of Technical

фессор, зав. кафедрой «Электромеханика и электри- Sciences, professor, head of department «Electromecanics

ческие аппараты», Южно-Российский государственный and Electric Devices», Platov South-Russian State Polytech-

политехнический университет (НПИ) имени М.И. Платова. nic University (Novocherkassk Polytechnical Institute).

Тел. (86352)55-1-13. E-mail: [email protected] Ph. (86352)55-1-13. E-mail: [email protected]

Проведена оценка эффективности метода пространственных интегральных уравнений при расчете квазистационарных электромагнитных полей электротехнических устройств. Для этого в пространственно-временной области решалась задача расчета распределения напряженности магнитного поля в соленоиде с прямоугольным ферромагнитным сердечником. На основе проведенного ранее аналитического решения модельной задачи по расчету распределения электромагнитного поля в трехмерной постановке для системы, включающей в себя проводящие ферромагнитные участки, магнитные свойства которых определяются кривой намагничивания B(H), и обмотки, питаемой от источника тока, построен алгоритм расчета с использованием метода пространственных интегральных уравнений. Проанализированы возможные причины, влияющие на точность расчета. Предложенная модель расчета квазистационарного электромагнитного поля, основанная на методе пространственных интегральных уравнений, позволяет получать достаточно точные результаты. На её основе можно не только проводить анализ распределения электромагнитного поля внутри ферромагнитных областей, но и вычислять интегральные характеристики: потокосцепление, электромагнитную силу и др. Последние можно получать с приемлемой точностью и при укрупненной пространственной дискретизации. Это может быть полезным при проведении экспресс-анализа рабочих процессов в электромагнитных устройствах.

Ключевые слова: токовая катушка; ферромагнитный сердечник; распределение напряженности магнитного поля; трехмерная картина поля; метод пространственных интегральных уравнений.

The efficiency of the method of spatial integral equations in the calculation of quasi-stationary electromagnetic fields of electrical devices. To do this, in the space-time domain, the task of calculating the distribution of the magnetic field in a solenoid with rectangular ferromagnetic core. On the basis of earlier analytical solutions of the model problem for the calculation of electromagnetic field distribution in three-dimensional setting for the system, which includes portions of the conductive ferromagnetic , magnetic properties are determined by the magnetization curve B(H) and winding fed from a current source, an algorithm for the calculation using the spatial integral equations. The possible factors influencing the accuracy of the calculation. The proposed model for the calculation of quasi-stationary electromagnetic field, based on the method of spatial integral equations allows to obtain accurate results. Based on it, you can not only analyze the distribution of the electromagnetic field inside the ferromagnetic regions, but also to calculate the integral characteristics: the flux linkage, electromagnetic force, etc. The latter can obtain an acceptable accuracy and spatial sampling enlarged. This can be useful for rapid analysis of work processes in electromagnetic devices.

Keywords: current coil; the ferromagnetic core; distribution of intensity of a magnetic field; three-dimensional picture of a field; method of the spatial integrated equations.

Для оценки эффективности метода пространст- технических устройств и определения влияния степе-венных интегральных уравнений (ПИУ) при расчете ни дискретизации численных моделей в пространст-квазистационарных электромагнитных полей электро- венно-временной области на точность получаемых

результатов решалась задача расчета распределения напряженности магнитного поля в соленоиде с прямоугольным ферромагнитным сердечником. Длина сердечника - I, сечение - а х Ь , проводимость материала - у , режим работы - ненасыщенный, с постоянной магнитной проницаемостью - , число витков обмотки на единицу длины сердечника - w0. Обмотка питается от источника тока /(^) (рис. 1).

dB

rot H = J ; rot E =--,

dt '

следствия из них

div B = 0 ; div J = 0, а также уравнения связи

B = ц0 (H + M) ; J = уE

и краевые условия

Rot E = 0 ; Rot H = 0; Div J = 0 ; Div B = 0;

(2)

(3)

|E (M, t )|-

->0; H (M, t )■

->0. (4)

Здесь Н, В - соответственно вектора напряженности и индукции магнитного поля; М - вектор намагниченности вещества; J = Jвих + Jст - в общем случае сумма плотностей вихревого и стороннего тока; Е -вектор напряженности электрического поля; ц0, е0 -соответственно магнитная и диэлектрическая проницаемости вакуума; у - удельная электрическая проводимость; п - единичный вектор нормали к поверхности, направленный из среды 1 в среду 2.

Введением скалярного фе и векторного А электродинамических потенциалов краевая задача (1) - (4) с помощью тождеств векторного анализа [1 - 4] сводится к следующей системе уравнений:

Рис. 1. Эскиз призмы прямоугольного сечения

Электромагнитное поле рассматривается в квазистационарном приближении. Расчетная область состоит из совокупности отдельных подобластей: ферромагнитной (объём У^) и одновременно проводящей

(объём Уу) среды, а также токонесущих элементов

(объём Уу). Для рассматриваемого примера (рис. 1)

Vц совпадает с Уу.

Ограничим длину сердечника таким образом, чтобы она была значительно больше размера поперечного сечения сердечника So (примерно на порядок). Обозначим объем полученной области V, поверхность, ограничивающую V, - как £ .

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

Используем основные уравнения Максвелла [1, 2]

A = Ь. 4л

Ш J вих

„z rrr rot M

dV + JJJ-dV -

- § ^^ dS + JJJ dVj

s r vj r

H

1

4л

3(M • r)r M

dV -

-JJJ ^^ dV-JJJ ^ dVj

V r V, r

(5)

= У

—--^ ds'

dt 4ле0 S r3

§ = 2eo M = (ц-1) H

dA §Ц- dS dt 4ле0 S r3

(1)

где Jвих - вектор плотности вихревого тока; Jст-вектор плотности тока катушки; М - намагниченность в точке, принадлежащей объёму dV ; V - проводящая и одновременно ферромагнитная область; £ - поверхность, ограничивающая область V; VJ -объем, занимаемый катушкой; £ - поверхностная плотность свободных зарядов; п - внешняя нормаль к поверхности £; ц - относительная магнитная проницаемость ферромагнетика; г - радиус-вектор, проведенный из элемента dV в «точку наблюдения» Q , в которой определяются значения векторов поля.

Г

5

3

r

r

J

Система (5) позволяет определить неизвестные вторичные источники: вихревые токи, намагниченность ферромагнитных деталей, наведенные на поверхности электрические заряды.

Алгоритм расчета электромагнитного поля

Решение системы (5) проводится численным методом последовательных интервалов времени А/. Производные по времени представляются в конечно-разностном виде. На каждом интервале времени решение системы (5) также проводится численно, для чего осуществляется пространственная дискретизация расчетных областей £, V , VJ в виде элементарных объемов AV (прямоугольных призм) и площадок А£ (прямоугольников), ориентированных по координатным осям х, у, z (рис. 1). Объем V разбивается на N призм, объем VJ на Жкат призм, поверхность £ на Ns площадок.

Система уравнений (5) для момента времени (k +1) в этом случае представляется в виде:

Л: = ^ 1 4л

Н =—

1 4л

3N 3N 3NKaT

Z glji JßHx i + Z g2j,Mt + Z g3J i=1 i=1 i=1

3NK

Z g5 jM -z g4J i - Z g6J

i =1 i =1 i =1

J =-Y

вих 1 l

Im =-2s0 Л

Л

Л(k+1)- Л1(к) + 1 7 , -+-— Z g 7 p

4лбо i=1

At

; (5a)

; (56)

; (5в)

Лх - Лх Лу - Лу

Лт(к+1) Лт(к) лт(к +1) лт( k)

nx + ny +

At

At

+n.

- Лг

m(k+1) m(k)

At

+n

1 Ns

+n

1 Ns

4ле

z g 7{m г, (k)+nz^~ z g 7(;m i

Z g 1 ^ p (k)+ j=1

1 Ns

x . L-, b ' jm^j(k)

4лбо j=1

0 j=1

4ле

о j=1

( z ) г jm j(k)

Mj =(^- 1 Hi

; (5r)

(5д)

ределяемых геометрией магнитной системы при её пространственной дискретизации.

Примеры вычисления матриц и их коэффициентов представлены в [5].

Поскольку геометрия системы неизменна, а магнитная проницаемость ц принята в данной модельной задаче постоянной, то уравнения (5б) и (5д) сводятся к одному уравнению вида

3N 1 1

Z g5 M--- Mj = —

ц-1

i=1

4л

3Nks

Z g4 jiJBUX i + Z g6 jiJCT i

i =1 i =1

или в матричном виде

G5 - D

H--1

M =

-1G4

4л

J вих +

-1G6

4л

(6)

где матрица D является диагональной с коэффициен-

том

1

в диагонали.

Ц-1

Расчет основан на решении системы уравнений (5а) - (5д) для момента времени (tk+1), где k = 0,1,2,... . Начальные условия имеют вид: /0 = 0 ,

Jвих (*0) = 0, M(/о) = 0, ) = 0, Л(^) = 0.

На каждом временном шаге в момент tk+1 задается плотность тока, например

(

Jст (tk+1 ) =

(tn+1)

1 - e т

2 Л

('k+1 ) .

где j = 1,2,3,... 3N - номер проекции векторов в точке наблюдения; i = 1,2,3,... 3N - номер проекции векторов в элементарном объёме; т = 1,2,3,. Ns - номер внешней элементарной площадки; k, k +1 - соответствующие номера моментов времени.

Полученные дискретные уравнения в векторном виде следует понимать как обобщенную форму записи для вычислений проекций векторов в декартовой системе координат х, у, z . Кроме этого g1 ^, g2^, g3 ^ ,

g 4 ^ , g5 ^ , g 6 ^ , g 7 ^ в полученной системе есть

элементы матриц G1, G2, G3, G4, G5, G6, G7 , оп-

Решение системы (5а) - (5д) осуществляется последовательной численной реализацией каждого из входящих в систему уравнений:

- уравнения (5б) и (5д) решаются совместно в виде (6), при этом заданными считаются Jст и Jвих .

Значения J ст задаются для момента времени (tk+1), а

Jвих берутся с предыдущего момента времени (/¡^);

- по результатам решения (6) (M (/¡,+1)) определяется Л (/¡,+1) - формула (5а);

- используя найденные Л(tk+1), Л(tk) и соотношение (5г), находим распределение % (tk+1). При этом

производные векторного потенциала определяются по формуле

Щг - -ъ (Л (tk+1)-Л ("-■));

- соотношение (5в) уточняет распределение

Jвих (tk+1 ) ;

- распределение напряженности магнитного поля Н +1) получается из (5б).

Для интегральной оценки погрешности численного расчета проверяется выполнение закона сохранения энергии

1

J

ст

J

1 Ш{Н(tk+i)B(tk+1) + Aвих (tk+1) Jвих (tk+1)}dV

2 V

=1 ju{a ct (tk+i) J CT (tk+i)} dv =

2 v

= 1Ш{Н CT (tk +1) B CT (tk+1)} dV,

где A C

H,

V

Br

векторы, рассчитанные от дей-

ствия стороннего тока, плотностью J с

H, B

W

1ЕНет (tk+1)Bct (tk+1 )-f1 е{н(tk+1)b(tk+1)-

2 N V 2 N

+ Авих (tk+1 ) Jвих (tk+1 )}

< 0,0001.

i (t )=

,2 \

1 - e

ад ад

ЕЕ

m=1 n=1

Xn

(-1)

16

2

-e

Ца У

nmn~ a

n, m = 1, 3, 5,... .

. nn . mn . , sin—x sin-y > +1

b '

(7)

В качестве исходных данных использовались следующие значения: а = Ь = 0,02 м; 1тах = 0,38 А ;

Т = 0,02с; х = 0,16 с2; у = 1,26 107

Ц a = = 6,283 -10

-4 .

/Омхм; Ц = 500; m = n = 21; толщина обмот-

ки 0,001м (для ПИУ); длина ферромагнитного сердечника I = 0,1м (для ПИУ).

Результаты, полученные по формуле (7), были приняты за эталонные при решении этой же задачи, но с использованием метода пространственных интегральных уравнений, описанного выше.

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

векторы, в которых учтено действие как стороннего, так и вихревого тока; А вих - вектор, обусловленный

действием только вихревого тока, плотностью J вих .

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

(

тока. Закон изменения тока / (/) =

2

1-e

1 max И ег0

производная i (t) представлены на рис. 2.

/'(i). А;

Рассмотрим решение задачи на примере линейной модели, предложенной в [6], поскольку в такой постановке алгоритм может быть проверен с наибольшей точностью. В работе было представлено аналитическое решение модельной задачи нахождения распределения напряженности магнитного поля с учетом влияния вихревых токов в бесконечно длинном соленоиде с ферромагнитным сердечником сечением а х Ь и проводимостью материала сердечника у , работающем в ненасыщенном режиме с постоянной магнитной проницаемостью ца, с обмоткой, питаемой от источника тока ), с числом витков обмотки на единицу длины сердечника w0 (рис. 1). Используем

полученное решение для исследования поведения электромагнитного поля в ферромагнитном сердечнике прямоугольного сечения при питании обмотки от источника тока ), изменяющегося по закону:

i (О. А/с 03

0.6

0.4

02

/ \ 2

/ ! \

1 ! \ 1 */

/ X %

02

0.4

М

0 8

t,c

Рис. 2. Закон изменения тока (1) и его производной (2)

Для удобства сравнения полученных результатов внутри ферромагнитного сердечника выбраны три контрольные точки с координатами: 1-я -(0,0043; 0,01; 0,05) - ближе к краю серединного сечения; 2-я - (0,0071; 0,01; 0,05)- промежуточная точка и 3-я - (0,01; 0,01; 0,05) - геометрический центр.

Результаты расчета при заданном токе представлены на рис. 3 и в табл. 1.

Я, А/м

V у

В этом случае распределение напряженности магнитного поля принимает вид:

Н (х, у, 0 = Щ) х

-0.1

/¿¿--4

2 v ■/ у/' f

i \ /// 3

--- 5

0.72

t, с

Рис. 3. Зависимости напряженности магнитного поля при изменении тока

График 1 - напряженность магнитного поля в центре катушки без учета влияния вихревых токов, 2 - с учетом влияния вихревых токов, 3 - напряжен-

т

2

т

X

ность магнитного поля, наведенная вихревым током, вычисленные по формуле (7); 4 - напряженность магнитного поля в центре катушки с учетом влияния вихревых токов и 5 - напряженность магнитного поля, наведенная вихревым током, рассчитанные по программе на основе ПИУ.

Погрешность для случая разбивки на элементы 7 х 7 х 33 = 16177, при изменении тока по закону

(

i (t ) =

л \

1 - e

К

контрольной

точке

V /

И,, (0,01;0,01; 0,05; /) с учетом влияния вихревых токов представлена на рис. 4.

Д, А/м

0,007

О 0.24 0.« 0.72 0.96 Рис. 4. Абсолютная погрешность

t, с

Потоки через поперечное сечение в центре сердечника, рассчитанные двумя разными методами, практически совпадают (рис. 5), отклонение не превышает 0,2 %.

Ф,х10

0 24

043

0.72

0.96

t, с

Рис. 5. Потоки через поперечное сечение в центре сердечника

На рис. 6 представлено распределение векторов плотности тока в катушке, а также распределение векторов плотности вихревого тока в поперечном сечении сердечника для случая разбивки на элементы 7 х 7 х 33 = 1617 .

Таблица 1

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

т

B

Шаг по времени / = 0,21 с при значении вихревого тока, близкого к максимальному значению (см. рис. 3, графики 3 и 5)

Координаты контрольной точки, заданной в программе на основе ПИУ с учетом толщины обмотки катушки (х, у, 2, /), А/м Модельная задача, формула (7) Программа на основе ПИУ, количество элементов 7 х 7 х 33 = 1617 Программа на основе ПИУ, количество элементов 9 х 9 х 33 = 2673

Значение A, % Значение A, %

И,(0,011; 0,011; 0,05; /) - центр, без учета вихревых токов 0,0915 0,0915 0,0 0,0915 0,0

И2(0,0053; 0,011; 0,05; /) - ближе к краю, с учетом вихревых токов 0,0730 0,0727 0,41 0,0729 0,14

И2 (0,0081; 0,011; 0,05; /) - ближе к центру, с учетом вихревых токов 0,0650 0,0647 0,46 0,0648 0,27

И2(0,011; 0,011; 0,05; /) - центр с учетом вихревых токов 0,0630 0,0626 0,63 0,0628 0,35

Шаг по времени / = 0,48 с , при значении абсолютной погрешности, близкой к максимальному значению, рис. 4

Координаты контрольной точки, заданной в программе на основе ПИУ с учетом толщины обмотки катушки (х, у, 2, /), А/м Модельная задача, формула (7) Программа на основе ПИУ, количество элементов 7 х 7 х 33 = 1617 Программа на основе ПИУ, количество элементов 9 х 9 х 33 = 2673

Значение A, % Значение A, %

И2(0,024; 0,011; 0,05; /) - центр, без учета вихревых токов 0,2890 0,2890 0,0 0,2890 0,0

И2(0,0053; 0,011; 0,05; /) - ближе к краю, с учетом вихревых токов 0,2810 0,2767 1,53 0,2778 1,14

И2(0,0081; 0,011; 0,05; /) - ближе к центру, с учетом вихревых токов 0,2780 0,2724 2,0 0,2743 1,32

И2(0,011; 0,011; 0,05; /) - центр с учетом вихревых токов 0,2770 0,2703 2,4 0,2725 1,62

Таблица 2

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

N dt Контрольная точка, Hz(x, y, z, t), А/м Элементарный объем, с контрольной точкой в программе на основе ПИУ Я2(х, у, г, А/м Относительная погрешность, % Полное время счета, с/ч

Количество элементов 7 х 7 х 33 = 1617, количество решаемых уравнений 4851

1 0,0024 0,2770 0,2700 2,53 2942/0,81

2 0,0012 0,2770 0,2703 2,4 356^0,99

3 0,0006 0,2770 0,2711 2,1 5447/1,51

4 0,0003 0,2770 0,2714 2,02 8686/2,41

Количество элементов 9 х 9 х 33 = 2673 , количество решаемых уравнений 8019

5 0,0024 0,2770 0,2722 1,71 12135/3,37

6 0,0012 0,2770 0,2725 1,62 13695/3,80

7 0,0006 0,2770 0,2726 1,58 15255/4,24

8 0,0003 0,2770 0,2732 1,37 16815/4,67

Рис. 6. Скриншот распределения векторов плотности тока в катушке и вихревого тока в поперечном сечении сердечника

Проверка зависимости расчетных параметров от шага по времени dt

В качестве контрольной точки выбрана точка с координатами х = 0,011м; у = 0,011м; г = 0,05м; / = 0,48 с (точка в центре дросселя), закон изменения

тока /

(t ) =

2

1-e

Imax . Результаты расчета про-

граммы при различных временных шагах представлены в табл. 2.

Вывод

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

Литература

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

2. Тамм И.Е. Основы теории электричества. М; Л., 1946. 660 с.

3. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М., 1970. 720 с.

4. Курбатов П.А., Аринчин С.А. Численный расчёт электромагнитных полей. М., 1984. 168с.

5. Подберезная И.Б. Применение пространственных интегральных уравнений для расчета квазистационарных электромагнитных полей в электромеханических устройствах // Изв. ЮФУ. Техн. науки. 2014. № 3(152). С. 250 - 264.

6. Подберезная И.Б., Ершов Ю.К., Павленко А.В. Аналитический расчет распределения магнитного поля в бесконечно длинной призме прямоугольного сечения с токовой обмоткой // Изв. вузов. Электромеханика. 2012. № 1. С. 12 - 15.

Поступила в редакцию

20 декабря 2013 г.

т