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

ИНФОРМАТИКА, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И УПРАВЛЕНИЕ

Филимонов С.И., ст. преп.

Белгородский государственный технологический университет им. В.Г. Шухова

РАЗРАБОТКА ИМИТАЦИОННОЙ МОДЕЛИ ПЕТЛИ ГИСТЕРЕЗИСА В ПРОГРАММНОМ КОМПЛЕКСЕ МЛТЬЛБ

sergey.filya@mail.ru

Разрабатывается математическая модель петли гистерезиса вида H=f(B). Предложена методика формирования частного несимметричного цикла перемагничивания. Реализована модель петли гистерезиса в модели силового трехфазного трансформатора в системе Matlab Simulink

Ключевые слова: петля гистерезиса, трехфазный трансформатор, симметричный и несимметричный цикл перемагничивания.

Моделирование петли гистерезиса, а так же частных петель для силовых трансформаторов и нелинейных индуктивностей в программном комплексе Matlab 8тиПпк возможно с помощью блока Powergui. Однако получающиеся при этом частные петли по форме значительно рас-20 15 10 5

$-5 -10 -15

-20

ходятся с фактическими. На рис. 1, а представлено семейство частных симметричных петель, которые при изменении напряженности магнитного поля (тока) практически не изменяют своей ширины, в отличие от фактических петель, показанных на рис. 1, б. [3]

8 4

гЧ

Ь 0

-4

-1.5

-0.5

/°А

0.5 1

1.5

-8 -11

^кА/м б

Рис. 1. Семейство частных петель гистерезиса (а) при моделировании, (б) фактические

Следовательно, для исследований некоторых специальных режимов работы трансформаторов необходимо использовать более корректную модель гистерезиса. Причем эта модель должна формировать как частные симметричные, так и несимметричные циклы перемагничивания. На сегодняшний день таким требованиям отвечают модели Джилса - Аттертона, Джона Чана и др. [2]. Наиболее простой для реализации является модель Джона Чана, так как она использует только основные параметры петли гистерезиса, такие как коэрцитивная сила (Hc), остаточная индукция и индукция

насыщения (Bs). Эти параметры для электротехнических сталей (ЭТС) легко найти, например, в справочнике [1].

Разработка гистерезисной модели сердечника для программного комплекса Matlab Sim-ulink, правильно реализующей частные петли гистерезиса, является актуальной задачей.

В модели гистерезиса Джона Чана верхняя и нижняя кривые петли предельного цикла пе-ремагничивания описываются соответственно уравнениями 1 и 2 [2].

в (т В*'(Н+ «( )" Н+нс\+Нс \в! /вг-1)

Вап (Н) =:

В, •( Н - нс)

+ Мо- Н,

- + Мо-Н,

(1)

(2)

Н -Нс\ + Нс-(В, /Вг-1)

где /л0= Ая -10 -7 - магнитная проницаемость

вакуума, Гн/м.

Вторым слагаемым в выражениях 1 и 2 можно пренебречь в силу малых значений.

Начальная кривая намагничивания (КН) при моделировании вычисляется как половина суммы верхней и нижней кривых по выражению 3:

В_

г^.Вцр (Н) + Вп (Н)

ЛН) ~ ^

(3)

2

а

Таким образом, предельный цикл перемаг-ничивания имеет вид, представленный на рис. 2.

вВ< Br; i

ш m^

Bup(H) Jhc H

Рис. 2. Петля гистерезиса и начальная кривая намагничивания

Однако в действительности предельный цикл - это один случай перемагничивания и чаще всего трансформатор работает по частным циклам, например, в режиме холостого хода получается симметричный частный цикл. В рассмотренной в статье [2] методике предлагается строить частные циклы путем переноса нижней и верхней ветвей предельного цикла соответственно вверх и вниз. При этом формируется траектория D-E-F-G-F, представленная на рис. 3, а. Участок ББ получается путем переноса нижней предельной ветки на величину Bd. Для моделирования КН были выбраны параметры ЭТС марки 3414 [1].

Частный цикл перемагничивания будет симметричным при 0 < Bd< Br. При Bd>Br участок Б-Б должен пройти выше точки А начала

1.5"

1.2

0.с

0.6

cq

-0.3 -0.6 -0.9 -1. -1.5

-120-100-80-60-40 -20 0 20 40 60 80 100 120 Я, А/м

координат. В этом случае перенос нижней ветви предельного участка вверх приведет к ломаной несимметричной частной петле, так как точка Б окажется на верхней ветви. Для получения ровного несимметричного цикла, когда, например, происходит изменение H от -40 до +120 А/м (рис. 3, б), была разработана методика, которая отличается от предложенной в статье [2] более удобной реализацией.

При увеличении напряженности до величины H = 0 (участок Б-Б) нижняя ветвь частного цикла формируется простым переносом нижней ветви предельного цикла как в предыдущем случае. Затем в точке H = 0 необходимо зафиксировать величину текущей индукции B. При дальнейшем увеличении H до +120 А/м необходимо сохранять пропорцию между расстояниями до предельной петли и до средней (начальной) кривой. Согласно рис. 3, б получается пропорция 4 из соответствующих отрезков, а кривые при изменении H от 0 до + Hmax и до - Hmax описываются соответственно уравнениями 5 и 6:

AE AB.,

A1E1 A' B,'

(4)

Вргт(Н) = (Бир(Н)-Бте,(Н})• Ба /4 + БяЫ(Н), (5) Б^(Н) = Бтс1 (Н) - (Бте,(Н) - Бп(Н)) • Ба / Бг, (6)

где Ба = АЕ - значение индукции при H = 0.

Таким образом, образуется частный цикл по траектории Б-Б-К-Вг-Б.

1.5 1.2

0.9

0.6 * °

-0.3 -0.6 -0.9 -1.2 -1.5

Br'

Br —""" к

D 1 / E EV7 / A' /

""A A f

с

Рис. 3. (а) симметричный, (б)

Перейдем к реализации рассмотренной модели петли гистерезиса в модели силового трехфазного трансформатора в Matlab Simulink. Из множества библиотечных моделей наиболее подходящей для внедрения, а так же для исследований, связанных с магнитной системой трансформатора, является «Three-phase Two-Winding Saturable Transformer - Core Type». В отличие от других моделей здесь обмотки

-120-100-80-60-40 -20 0,20 40 60 80 100 120 Я А/м

б

несимметричный частный цикл

трансформатора разных фаз магнитосвязаны. Особенностью структуры этой модели является реализация КН в виде функции H=f(B). То есть, для применения модели петли гистерезиса, необходимо получить из имеющихся уравнений B(H) уравнения вида H(B).

Воспользовавшись математическим вычислительным комплексом Maple, после простых преобразований уравнений 1 и 2 имеем:

а

Н (В) = 2Н -Н-((В -ВгУ В~вв + Вг- 1В1), (7)

и с с В -(IВ - В,)

((В - В )-В - В В + В -|В|) , ч Нпп(В) = -Нс-г) „ ' В + г ' 1) . (8)

Решением задачи по преобразованию начальной КН, описанной уравнением 3, являются восемь корней, три из которых 9, 10, 11 определяют ее характер:

(В) =

Вг - (IВ - В,)

Нс - (В2 + 2В - (Вг - В,)- ВВ•( 4В2 - 8ВВ, + 5 В2)+ ВА - 2ВгВ3)

Нте, 2 (В ) =-

-( в

2Вг •( в - В,)

НсВ, - (Вг - В, + В, - Вг )2 + 4В2)

2ВВ.

Нтва 3( В ) =

Нс -(2В - (Вг - В,)-В2 + ВВ +^В1 -(4В2 + 8ВВ, + 5В2) + В4 - 2В,Д3)

-2ВГ -(В + В,)

(9)

(10) (11)

Причем каждое из полученных уравнений действует на своем интервале B, то есть:

НтеЛ( В)

Нтгй (В) = Нтгй 2 (В) \Г

Нтгй 3(В)

>Втга (Нс );

^ (Нс ) < В < Втеа (Нс );

<-Втга (Нс У;

(12)

Здесь Bmed(Hc) величина магнитной индукции, определяемая по формуле 3 при H = Нс:

Втей (Нс ) ="

В.-Н,

(

2Н + Н -

(13)

В-1

V В /

Общий вид обратной петли гистерезиса и начальной КН представлен на рис. 4, а. Отдельный вид начальной КН, образованной тремя графиками, представлен на рис. 4, б.

1 Ншеё1(Б) Ншеё2(Б) \

-Бшеё(Нс) \

1 1 — ~—" _ 1 1 -► Б

Бшеё(Нс)

\ Нте<В(Б)

Рис. 4. Общий вид (а) обратной

Теперь, чтобы получить частный симметричный цикл перемагничивания необходимо при увеличении индукции В сместить влево кривую Н^В) на величину Ва при уменьшении В сместить вправо кривую Нир(В) на ту же величину. В результате формируется частный симметричный цикл перемагничивания по траектории Б-Б-Р-О-Р (рис. 5, а).

Для получения частных несимметричных петель применяется аналогичная методика, как и

б

петли гистерезиса, (б) начальной кривой

с прямой петлей. Однако теперь необходимо вычислять Н. Поэтому когда индукция В увеличивается и Н меньше 0 кривая Нир(В) просто смещается вправо (участок ББ на рис. 5, б), а когда Н начинает входить в интервал от 0 до +Нтах, то необходимо соблюдать соотношение 4. В итоге формируется цикл по траектории Б-Б-К-Бг-Б (рис. 5, б).

а

120 100 80 60 40 120

-20 -40 -60 -80 -100 -120

ц

/ // А' . / в, Ч ... 1-

Г °

¡/с

120 100 80 60 40 ^20

-20 -40 -60 -80 -100 -120

к!

АЛ^........ •••' " * J

7/Е' А /Е У Вг' Вг

Б

ус

-1.5 -1.2 -0.9 -0.6 -0.3 0 0.3 0.6 0.9 1.2 1.5

В. Тл

-1.5 -1.2 -0.9 -0.6 -0.3 0 0.3 0.6 0.9 1.2 1.5

В. Тл

а б

Рис. 5. (а) Симметричный, (б) несимметричный частный цикл обратной петли

Получим уравнения для кривой при несим- Однако только два определяют характер кривой.

метричном цикле, когда величина Н изменяется Причем уравнения 14 и 15 описывают поведе-

от 0 до +Нтах. Для этого выразим переменную Н ние кривой в первом квадранте, а 16 и 17 - в

из уравнений 5 и 6. При преобразовании каждо- третьем квадранте: го из уравнений появляется по восемь корней.

Ищы 1( В ) = ■

ИЛ -(в,-Вг + Вг - В, )2 + (4В,Ва1р - 4 В )•( Ва1рВг - В)

ИгШ 2(В) =

2В •( В В , - В)

г \ 5 а1.р !

И А •( В,ВГ + 2В •( В, - Вг)-В2 + 4В2В2г + ВВ .(5ВД - 2В, - 8ВгВ)+4В,ВгВа (В,В + Вг (В, - В)-В2))

И1ф( В ) = ■

2Вг ■ (В, -В)

ИВ АВ -В + .1( В -В )2 + (4В В , + 4В )•( В ,В + В)

с , [г г / V 5 а1р ! \ а1р г !

2В • (В В , + В )

^ \ 5 а1р !

ИШ(В) =

ИсВ, -(в,2 + 2В • (В, - Вг)-ВВг+ 4В2Вг2 + В,Вг-(5В,.Вг - 2В, + 8ВгВ)+4ВАВа(Вг (В, + В)-В,В - В,)

2Вг- (В, + В)

(14)

(15)

(16) (17)

Общий вид части обратной петли гистере- образованная двумя графиками по уравнениям зиса и начальной КН в первом квадранте, а так 14 и 15, при изменении Н от 0 до +Нтах, пред-же кривая частного несимметричного цикла, ставлены на рис. 6.

В, В

Рис. 6. Петля гистерезиса и частный несимметричный цикл в первом квадранте

0

Как видно из рис. 6 графики Нгщ1й1(Б) и Нгщш(Б) пересекаются в точке А с координатами Бгщм01 и НгщЬЮ1. Следовательно, общий график Нгщы(Б) изменяется сначала по кривой НгщЬН(Б),

затем в точке А переходит на кривую Нгщш(Б). Координаты точки А в первом и третьем квадранте описываются соответственно уравнениями 18, 19, 20, 21:

B.

BsBr-(1 + B„)

right 01

В + B„

H •(Bs + Br)•

(

H

Br - Bs +1

|B.2 - 2BBB + B2h

s r a s

B + B

r s

right 01

2Br •(Br -BsBa) B B •(l + B )

R = s r У a'

01 " B + B

s r

H •(B + B )

c \ s r t

(

H

B -B +J

rs

B2 - 2B B B + B

2\\

s r a s

B + B

ght 01

2Br •( BsBa- Br)

(18)

(19)

(20)

(21)

Полученную модель обратной петли гистерезиса можно представить в виде S - функции, описанной в Matlab. Для ее применения выбранный трансформатор необходимо преобразовать, убрав блоки с КН, а вместо них поставить блоки с S - функцией.

Рассмотрим подробнее работу модели трехфазного трансформатора с разработанной петлей гистерезиса. Помимо B-H характеристики стали в данной модели учитывается геометрия магнитной цепи. Структура модели блока «Three-phase Two-Winding Saturable Transformer - Core Type» представлена на рис. 7. Выводы первичных обмоток обозначены индексом 1, а

вторичных - 2. Магнитная связь между обмотками, которые располагаются на едином сердечнике, моделируется в подсистеме «Core». На вход подсистемы подаются значения магнитных потоков «Fluxes», а на выходе получаются магнитодвижущие силы МДС. Магнитные потоки в такой модели трансформатора вычисляются из напряжений «Voltages», которые снимаются через датчики напряжений «V», подключенных к сопротивлениям Rm. Для управления источниками тока «Imag», которые моделируют ток намагничивания в обмотках, полученные МДС преобразуются в токи «Currents».

Рис. 7. Структура модели трехфазного трансформатора

В подсистеме «Core», которая представлена на рис. 8, применяются подсистемы (блоки магнитного сопротивления) нелинейной магнитной цепи, состоящей из трех стержней («L1 1», «L12», «L13») и четырех частей ярем («L21», «L2 2», «L2 3», «L2 4»), а так же блоки воздушных элементов «Air», моделирующие потоки рассея-

ния и поток нулевой последовательности. Установленные датчики напряжения «FMM» измеряют величину МДС каждой из шести катушек, а результаты поступают на блок «State - Space», который решает дифференциальные уравнения обобщенной электрической машины.

Рис. 8. Структура подсистемы «Core»

В подсистемах нелинейной магнитной цепи стандартной модели «Three-phase Two-Winding Saturable Transformer - Core Type» используется нелинейный блок «Look - Up Table», в котором

задается В - Н характеристика КН трансформатора, рис. 9, а. Разработанная модель петли гистерезиса в виде S-функции представлена на рис.9,б блоком «ЕипеНоп1».

XYGrapM

а б

Рис. 9. Структура подсистемы нелинейной магнитной цепи (а) с КН, (б) с петлей гистерезиса

Произведем количественную оценку работы моделей трансформаторов с гистерезисом и без него. Примем к рассмотрению силовой трансформатор типа ТМЭ-400/10 со следующими параметрами: Sн = 400 кВА, ин1 = 6300 В, 1н1 = 38.49 А, 1хх = 1.52 %, АРх = 820 Вт, АРКЗ = 6062 Вт. Схема и группа соединения обмоток Б/Ун-11. Магнитопровод трансформатора выполнен из стали ЭТ-А-3408 ГОСТ 21427.1-83. Длины средних линий магнитопровода - междуосевое расстояние по соседним стержням Ь2 = 325 мм, по ярмам Ь1 = 635 мм. Площадь сечения магни-топровода S = 225.264 см2. Количество витков обмотки НН w2 = 27, ВН w1 = 425.

Рассмотрим пуск трансформатора на холостом ходу. Структурная схема модели представлена на рис. 10. Здесь модель трансформатора представлена подсистемой «О/У», которая подключается к блоку «Трехфазный источник напряжения» через блок «Выключатель». Для пуска трансформатора подается сигнал с блока «Задание» на включение выключателя. Величина и форма фазных токов измеряется в блоке «Трехфазный измеритель и и 1» и выводятся на блок «Осциллограф». Результаты измерений для модели трансформатора без петли гистерезиса и с петлей представлены на рис. 11 а и б.

П

Завам«

Солп1 С ОТО? СоопЗ уаьс А 1аьс В л >

Л» ■-- 1 А1 И А2-

V III ш С Трея&амый источим Сопи^0""4 Ьшжхател ИЗ С с Грехфаэиый меритель и ь 82 В2- С1 С? с

От/

Рис. 10. Структурная модель пуска трансформатора

Рис. 11. Осциллограммы токов в обмотках ВН в модели (а) без петли, (б) с петлей

Как видно из рис. 11, б бросок токов намагничивания при включении оказывается гораздо больше, чем в модели без петли гистерезиса, рис. 11, а. Для фазы А их величина составляет соответственно 375 А и 143 А. Помимо этого видно, что затухание токов в модели с гистерезисом происходит явно быстрее. За время 0.6 с ток на рис. 11, б изменяется от 375 А до 50 А (в 7,5 раз), в то время как ток на рис. 11, а за то же время изменяется от 143 до 50 А (в 3 раза). Это объясняется тем, что наличие петли гистерезиса увеличивает потери мощности на перемагничи-вание.

400

Об адекватности разработанной модели силового трансформатора свидетельствуют данные, полученные экспериментальным путем. В ходе эксперимента комплектное распределительное устройство КРУ подключало сразу два трансформатора ТМЭ-400/10 со схемой соединения обмоток Б/Ун и У/Ун. Осциллограммы суммарных фазных токов первичных обмоток снимались с трансформаторов тока (ТТ), установленных в КРУ. Результаты измерения для фазы А представлены на рис. 12.

320 240 160 80 ^ 0 -80 -160

Г

0 0.5 1 1.5 2 2.5 3 /. с-

Рис. 12. Осциллограмма тока в фазе А при эксперименте

Аналогично эксперименту структурная модель пуска (рис. 10) была дополнена вторым трансформатором. Для сравнения на рис. 13, б представлена осциллограмма такой же фазы А, как и в эксперименте. Как известно ТТ не пропускают постоянную составляющую измеряемого тока, поэтому на осциллограмме рис. 12 наблюдаются отрицательные полуволны, в от-

личие от смоделированных осциллограмм на рис. 13, б, где использовались датчики тока ДТ. Для оценки адекватности было проведено восстановление кривой первичного тока с помощью интегрирующего - форсирующего звена с равными постоянными времени, результат приведен на рис. 13, а.

400 320 240 160 80 0 80

О 0.5

1.5

/, с

2.5

400 320 240 160

80 , О -80

О 0.5

1.5

/, с

2.5

б

Рис. 13. Осциллограммы тока в фазе А (а) восстановленного, (б) смоделированного

Как видно на рис. 13 в обоих случаях за время 0,5 с ток затухает до величины 70 А. Отличие состоит только в первом броске тока. Так на рис. 13, а бросок тока соответствует 400 А, когда на рис. 13, б он достигает величины 320 А. Это объясняется тем, что применяемый ТТ в первый момент времени насыщается. Для количественной оценки всего переходного процесса были посчитаны амплитуды токов за 16 периодов и построены затухающие экспоненциальные зависимости без первого пульса. В результате оценки средних значений было выявлено, что экспериментальный ток затухает на 5.5 % быстрее, чем смоделированный, что говорит об адекватности имитационной модели.

Выводы:

1. Разработана модель обратной петли гистерезиса по преобразованным уравнениям Джона Чана и др. для сердечников трансформаторов и нелинейных индуктивностей;

2. Разработана и предложена методика формирования частных симметричным и

несимметричных циклов перемагничивания стали сердечника;

3. Получены рекомендации по применению модели петли гистерезиса в существующих моделях трансформаторов программного комплекса Matlab^;

4. Доказана адекватность математической модели.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Молотилов Б.В., Миронов Л.В., Петренко А.Г. и др. Холоднокатанные электротехнические стали: справочник. М.: Изд. Металлургия, 1989. 168 с.

2. Володин В. Гистерезисная модель нелинейной индуктивности симулятора LTspice // Силовая электроника. 2010. №1. С. 56-60.

3. Михайлова М. Ю., Яковлев А. О., Нестеров А. М. Теоретические основы электротехники: методические указания к выполнению лаб. работ для студентов специальности 140211-Электроснабжение. Белгород: из-во. БГТУ им. В. Г. Шухова, 2012. 54 с._

Filimonov S.I.

HYSTERESIS LOOP SIMULATION MODEL DEVELOPMENT IN THE MATLAB PROGRAMM COMPLEX

Is developed mathematical model of hysteresis loop H=f(B) kind. Method of formatting partly asymmetric cycle of magnetization reversal is proposed. There is realized model of three-phase transformer hysteresis loop in the Matlab Simulink system.

Key words: hysteresis loop, three-phase transformer, symmetrical and asymmetrical cycle of magnetization reversal.

1

2

1

2

3

3

а

Филимонов Сергей Игоревич, старший преподаватель кафедры электроэнергетики Белгородский государственный технологический университет им. В.Г. Шухова. Адрес: Россия, 308012, Белгород, ул. Костюкова, д. 46. E-mail: sergey.filya@mail.ru