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

ЭНЕРГЕТИКА И ЭЛЕКТРОТЕХНИКА

УДК 621.3.017.31

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

ВХОДНОГО СИГНАЛА

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

Подберезная Ирина Борисовна - канд. техн. наук, доцент, кафедра «Электромеханика и электрические аппараты», Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова. Тел. (86352)55-1-13. E-mail: [email protected]

Ершов Юрий Константинович - канд. техн. наук, доцент, кафедра «Теоретическая электротехника и электрооборудование», Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова. Тел. (86352)55-3-08.

Павленко Александр Валентинович - д-р техн. наук, профессор, зав. кафедрой «Электромеханика и электрические аппараты», Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова. Тел. (86352)55-1-13. E-mail: [email protected]

Podbereznaya Irina Borisovna - Candidate of Technical Sciences, assistant professor, department «Electromecanics and Electric Devices», Platov South-Russian State Polytechnic University. (NPI). Ph. (86352)55-1-13. E-mail: podbereznayaib @mail.ru

Yershov Yury Konstantinovich - Candidate of the Technical Science, associate professor department «Theoretical Electrical Equipment and Electric Equipment», Platov South-Russian State Polytechnic University (NPI). Ph. (86352)55-3-08.

Pavlenko Alexander Valentinovich - Doctor of Technical Sciences, professor, head of the department «Electromecanics and Electric Devices» Platov South-Russian State Polytechnic University (NPI). Ph. (86352)55-1-13. E-mail: [email protected]

Целью данной работы стало исследование эффективности метода пространственных интегральных уравнений применительно к расчету несинусоидальных периодических полей в пространственно-временной области. Исследования проводились на примере расчета квазистационарного электромагнитного поля в призме прямоугольного сечения при различных формах гармонически изменяющегося сигнала в виде нарастающих синусоиды, треугольного сигнала и прямоугольных импульсов. Проблемы численного дифференцирования, связанные с ошибками, вносимыми в вычисления округлением и усечением, удалось преодолеть посредством применения достаточно точного метода Риддерса, который дает примерно 7 - 8 верных значащих цифр первой производной. Его применение позволило решить задачу для частот до 50 - 300 Гц и с приемлемой точностью отслеживать процессы, происходящие в сердечнике с обмоткой, питаемой от источника периодического несинусоидального тока.

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

The aim of this work was the study of the effectiveness of the method of spatial integral equations for computation of non-sinusoidal periodic fields in the space- time domain. Studies were conducted on the example of the calculation of quasi-stationary electromagnetic field in a rectangular prism with various forms of harmonically varying signal in the form of increasing sinusoidal signal of triangular and rectangular pulses. Numerical differentiation problems related to errors introduced in the calculation rounding and truncation, was overcome by applying a sufficiently accurate method Ridders which gives about 7-8 correct significant digits of the first derivative. Its use has allowed to solve the problem for frequencies up to 50 - 300 Hz, and with reasonable accuracy to track the processes occurring in the core with winding

Keywords: quasistationary electromagnetic field; current coil; not sinusoidal harmonious signal; distribution of intensity of a magnetic field.

В работах [1, 2] были рассмотрены аналитический расчет и реализация метода пространственных интегральных уравнений (ПИУ) на примере расчета квазистационарного электромагнитного поля в призме прямоугольного сечения.

Продолжим исследования и рассмотрим в качестве источника тока ) периодическую несинусоидальную функции времени.

Решение предполагает разложение периодической несинусоидальной функции времени в ряд Фурье. Этот подход эффективен в случае линейных характеристик материала. В электромагнитных системах с нелинейными характеристиками предполагается последовательное решение во времени с шагом дискретизации At.

Целью данной работы является исследование эффективности метода ПИУ применительно к расчету несинусоидальных периодических полей в пространственно-временной области.

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

Электромагнитное поле рассматривается в квазистационарном приближении. Расчетная область может состоять из совокупности отдельных подобластей, включающих ферромагнитную (объём V^) и проводящую (объём Vy) среды, а также токонесущие элементы (объём Vj). Для рассматриваемого примера (рис. 1) V совпадает с Vy.

Выберем длину сердечника такой, чтобы она была значительно больше размера его поперечного сечения S0 (примерно на порядок). Полученную область обозначим V, поверхность, ограничивающую V, через S .

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

Если обмотка питается от источника синусоидального тока i (t) = Imaxsinrot, то i (t)= Imaxrocos rot, аналитическое решение имеет вид [1]:

Hz (X ^t) = 11 maxro c0s roT

1 -ЕЕ

16

-a(i-т)

e v ' x

1 n=1 l nmn

. nnx . mny x sin-sin-

. nnx . mny x sin-sin-

d T = 1 max ro

1 ОТ ОТ

— sin rot - Е Е

m=1 n=1 l nmn

16

2

b а2 +ro

x™2 ('

а cos rot + ro sin rot -ae

'))],

(1)

где а =

А п,т » (пп] (

-; А пт =1 — I +1 — I - характеристика У ' V а ) V Ь )

ческие числа оператора Лапласа.

Используем в качестве исходных данных размеры сердечника а, м, Ь, м, частоту f ,Гц, период

Т = 1f, с, круговую частоту го = , амплитудное

значение тока в обмотке /тах, А, число витков на единицу длины , электрическую проводимость материала у , относительную магнитную проницаемость ц и следующие их значения: а = Ь = 0,02 м;

w0 = 1; f = 50 Гц; T = 0,02c;y = 1,26•Ю7^

(Ом • м)

цa = = 6,283•Ю-4; ц = 500 ; m = п = 21.

Принимаем, что ток изменяется по закону

i (t ) =

/2 Л

1 - e Т0

Imax sin rot (рис. 2), тогда

21 te т0 i (t )= 21 max te-sin rot +

,2 Л

1 - e T°

I max ro cos rot,

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

формула (1) примет вид

X

a

ro

a

Г

2

0

Hz (x, y, t) =

( x2

=J

21 max X в T° •

-max-Sin ЮХ +

2

1 - в T°

1 maxЮ C0S rox

1 -ЕЕ

16

1 n=i l nm%

-a(t-x) . nnx . mny

в v 'Sin-Sin-—

dx, (2)

где Imax = 0,3А ;То = 0,16с2.

Аналитическое вычисление интеграла в (2) достаточно трудоёмко, но его можно взять, используя численные методы. В данном случае он вычислялся по квадратурным формулам Гаусса - Кронрода.

i(t), А 0,267

0,133

0

-0,133 -0,276 -0,400

0 0,01 0,02 0,03 0,04 0,05 t, с Рис. 2. Закон изменения тока в катушке

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

На рис. 4 более детально представлены процессы, протекающие в центре сердечника (точка 4 рис. 3), рассчитанные с использованием формулы (2). Здесь график 1 - без учета влияния вихревых токов; 2 - с учетом влияния вихревых токов; 3 - напряженность магнитного поля, наведенная вихревым током.

И, А/м 0,12 0,09 0,06 0,03 0

-0,03 -0,06 -0,09 -0,12 -0,15

s

t / \

t \

1 /_ / \

2 '__/ _ _ _

3

0 0,005 0,01 0,015 г, с

Рис. 4. Изменение во времени напряженности магнитного поля в центре сердечника

Для расчета распределения магнитного поля на основе ПИУ рассматривалась в работе [2] система уравнений (3), представленная в преобразованном для момента времени (k +1) в виде (3а) - (3д).

/

Л

1

2 3 *---*----- 4 /

< y k. 7

0,24 0,18 0,12 0,06 0

-0,06 -0,12 -0,18 -0,24 -0,30

Л / \

2 3 Л\ / Ml f — J \ \

... / / \ -Л-........ 1 i_J

1 \ ч —1—У ■ • \ Ч/У . / \ \ 1 / i x i V ч 1

4 \ / I ^ 1 \

0,01

0,02

0,03

0,04

0,05

t, с

Рис. 3. Распределение напряженности магнитного поля с учетом влияния вихревых токов по ширине сердечника: 1 - И2 (0;0,01; 0,05; г); 2 - И2 (0,001; 0,01; 0,05; г); 3 - И2 (0,005; 0,01; 0,05; г) ; 4 - И2 (0,01; 0,01; 0,05; г)- графики изменения напряженности магнитного поля с учетом влияния вихревых токов, рассчитанные по формуле (2)

X

X

0

X

a

z

0

x

A = *> 4л

JJJ £вих dv + JJJ ^^ dV -

.v r V r

-§ ^^ dS +JJJ ^ dVJ

H

r 1

4л

JJJ

3(M • r)r M

dV -

-JJJ ^^ dV-JJJ ^ dVj

V r

dA

J вих = УЕ = У| " grad Фе | =

= У

—--^ #§f dS^ dt 4ле0 S r3

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

-+-L- #§f

dt 4ле0 S r3

(3)

где М - намагниченность в точке, принадлежащей объёму V ; Е - напряженность электрического поля;

Фе = —1— <й — dS - скалярный электрический потен-

4Л60 S г

циал простого слоя электрических зарядов с плотностью — на границе S области V ; п - внешняя нормаль к S; ц - относительная магнитная проницаемость ферромагнетика; г - радиус-вектор, проведенный из элемента dV в «точку наблюдения» Q, в которой определяются значения векторов поля.

A, =^0 1 4л

Н = —

1 4л

3N 3N 3^ат

Е g1 J i +Е g 2 j,Mt + Е g 3 J i =1 i=1 i =1

3N 3N 3Nrax

Е g5 ßMt -Е g 4 J г - Е g6 J i =1 i =1 i =1

; (3а)

;(3б)

J . =-y

вих 1

7 ji ?,

л

At

4ле

0 i=1

; (3в)

§m =-2sС

Ax - Ax Ay - Ay

Am(k+1) Am(k) m(k+1) лш( k)

nx + ny +

At

Am (k+1) Am (k) 1 ^

+nz ( ^-— + Пх

At

4ле,

At

е g 7(m § i (k)+

0 1 =1

1 ns , 1 nS ,,

+n„-ЕЕ g 7(ym) § 1 (k) + n^— ЕЕ g 7jm §

•У 4ле0 1^1(k^ z 4ле

0 1=1

M1 =(Ц- 1) H1;

m (k)

; (3г) (3д)

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

Дискретные уравнения в векторном виде понимаются как проекции векторов в декартовой системе

координат X y, г , а g1 }1, g2}1, g3 }1 , g4}1, g5 }1,

g 6 ц , g 7 ц есть элементы матриц G1, G2, G3, G4,

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

Для повышения точности расчетов производная по времени векторного магнитного потенциала в (3в), (3г) определялась по методу Риддерса [3], который в данном случае был модифицирован. Основные идеи метода и результаты исследования приведены в прил. 1.

Результаты расчетов представлены на рисунках ниже. На рис. 5 график 1 - напряженность магнитного поля в центре сердечника (см. точку 4 на рис. 3) с учетом влияния вихревых токов, рассчитанная по программе на основе ПИУ, график 2 - то же, но полученная с использованием формулы (2) и принятая за эталонную. Н, А/м 0,004 0,003 0,002 0,001

0

-0,001 -0,002 -0,003 -0,004 -0,005

i

\\

/ Л / / \ \__ 2 1 \ \

у > \ ! \ \ 1 \ \ \\

^ / ' \ / \ \ / /

h 1 \ ' \ \ / / 1 t \\ w

//

ч1

р

0 0,012 0,024 0,036 0,048 t, с

Рис. 5. Изменение напряженности магнитного поля в центре сердечника с учетом влияния вихревых токов

Погрешность в вычислении амплитуды напряженности магнитного поля в центре не превышает 3 %, а сдвиг фаз 0,0024 с для случая разбивки на N = Nx х Ny х Nz = 7 х 7 х 33 = 1617 элементарных объемов [2].

На рис. 6 представлен треугольный сигнал (4), а на рис. 7 - результаты расчетов на основе ПИУ для этого случая

i(t) =

(

1 - e Т0

i — Е

0 2 Е

л п=1

k-1 (-1) ^

k2

sin k гot

(для k = 1,3,5.).

(4)

j

5

3

r

r

i(t), А

0,267 0,133 0

-0,133 -0,267 -0,004

0 0,01 0,02 0,03 0,04 0,05 г, с Рис. 6. Треугольный сигнал

H, А/м 0,267 0,133 0

-0,133 -0,267 -0,400

Л

1 з Л /\ / \ / \ / \ / /\

/ /х / N 2 \ // \

"--_ __. ■ \ / \ \ !

V \ / \ / \/

0,01 0,02 0,03 0,04 0,05 t, с

Рис. 7. Графики составляющих И при треугольном сигнале

На рис. 8 - сигнал в виде прямоугольных импульсов

(

2 \

i(t) =

i(t), А 0,25 0,20 0,15 0,10 0,05 0

1 - е т»

2 (я/2 ® 1 . k -я/ 2

10- 1Ч-+1Т sin—TT я l 4 k=1 k 2

(для k = 1,2,3...).

cos k rot

(5)

/1

0,01 0,02 0,03 0,04 0,05 г, с

Рис. 8. Сигнал в виде прямоугольных импульсов

На рис. 9 - результаты расчетов на основе ПИУ для (5).

H, А/м 0,267

0,133 0

-0,133 -0,267 -0,400

0

0,01 0,02 0,03

0,04

0,05 t, с

Рис. 9. Графики составляющих И2 при сигнале в виде прямоугольных импульсов

На рис. 7 и 9 график 1 - составляющая Иг напряженности магнитного поля с учетом влияния вихревых токов у края сердечника (точка 1 на рис. 3), график 2 - то же, но в центре сердечника (точка 4 на рис. 3), график 3 - составляющая Иг напряженности магнитного поля в центре сердечника, обусловленная вихревыми токами.

Исследования процессов, происходящих в сердечнике при питании обмотки от источника периодического несинусоидального тока, проводились при частотах f = 50^300Гц, расчет осуществлялся по методу ПИУ в пространственно-временной области. Рассмотренная модель позволила с приемлемой точностью отследить все процессы. При расширении границ частот периодических сигналов (свыше f = 50 Гц) численное дифференцирование (3в), (3г), основанное на методе Риддерса [3], обеспечило получение более точных результатов.

В следующих работах будут приведены результаты исследований, призванные оценить границы применимости метода ПИУ при вариациях объемной дискретизации расчетной области, проводящих и магнитных свойств материальной среды и частотного диапазона источников поля.

Приложение 1

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

Будем искать производные функции F (х) в точке х численным методом по формуле дифференцирования с центральной разностью

^х) = F(х + к) - F(х - к) + Ег , 2Н

1

0

где Ег = -1 (h)2 F (х) - погрешность дифференцирования.

Будем уменьшать шаг дифференцирования при каждой итерации в р раз, тогда

F (х) = А1 + Ег = А2 + ЕГ2 = А3 + Ег3 = А4 + Ег4 = ....

При этом

A =

F (x + h) - F (x - h) 2h

A2 =

F (x + h) - F (x - h) P P .

2

h

A3 =

h h F (x + 4) - F (x -h)

P P

2

h

Er Er2 Er3 2

Er2 Er3 Er4

Используем метод Ромберга (1955), применив полиномиальную интерполяцию Ричардсона (1911):

F (t) :

P2 A2 - Al p2 -1

=Bi.

Построим таблицу оценок по методу Ромберга для порядка оценки производной т = 4 (табл. П1, П2).

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

F (x) =-

--, той же, что была предложена в [3].

sin х - х

Вычислим F (1) с начальным шагом численного дифференцирования h = 0,01. Истинное значение -140,737735571297. В табл. П3 представлены результаты работы программы при р = V2, в табл. П4 -

при p = 2 , в табл. П5 - при p =

1 -^5

2

. Абсолютная

погрешность вычисления в первом случае составляет 1 • 10-12, а в двух других - и того меньше.

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

Таблица П1

Таблица оценок по методу Ромберга

x

2

Шаг, h h h h h

порядок —

оценки P P P P

m = 0 Ai A2 A3 A4 A5

m = 1 Bi B2 B3 B4

m = 2 Q C2 C3

m = 3 Di d2

m = 4 Ei

Таблица П2

Формулы оценок при различных порядках оценки

m = i m = 2 m = 3 m = 4

ll P m - i JA C = P2mB,+i - Bt J P2m -1 n2mC - C D =P +1 J P2m -1 JE = P - i JD

Таблица П3

Пример работы программы при p

Шаг, порядок оценки h h V? h h h

&)' (Л)' {л r

m = 0 141,678097131387 141,20636406831 140,971663667478 140,854603319728 140,796145400314

m = 1 140,734631005234 140,736963266645 140,737542971978 140,737687480901

m = 2 140,737740687116 140,737736207089 140,737735650542

m = 3 140,737735567085 140,737735571035

m = 4 140,737735571298

Пример работы программы при p = 2

Таблица П4

Шаг, порядок оценки h h 2 h (2)2 h (2)3 h (2)4

m = 0 141,678097131387 140,971663667478 140,796145400314 140,752333523322 140,741384777835

m = 1 140,736185846175 140,73763931126 140,737729564324 140,737735196006

m = 2 140,737736208932 140,737735581195 140,737735571451

m = 3 140,737735571231 140,737735571296

m = 4 140,737735571297

Пример работы программы при p = 1-V5 2 Таблица П5

Шаг, 2h 4h 8h 16h

порядок оценки h 1 -V? (1 -Я)' (1 -Я f (1-Г5)'

m = 0 141,678097131387 141,095457089646 140,874160547471 140,789814340667 140,757623381506

m = 1 140,735365740644 140,737391762813 140,73768551804 140,737728274614

m = 2 140,737728274614 140,737735697424 140,737735578309

m = 3 140,737735570461 140,737735571279

m = 4 140,737735571297

Литература

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

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

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

3. Ridders C.J.F. Advances in Engineering Software. 1982. 4(2) 75 - 76.

30 мая 2014 г.