Опишем подход к моделированию ПП в ТТ со стальным сердечником. Применение таких моделей ТТ необходимо для изучения поведения средств РЗиА при экстремальных переходных процессах, когда измерительные ТТ входят в насыщение и начинают работать со значительной погрешностью, что может обусловить неправильную работу устройств РЗиА. Ключевой особенностью моделирования является способучёта нелинейности основной кривой намагничивания стали сердечника ТТ.
Для начала на наиболее простом примере – на примере одиночно работающего ТТ – рассмотрим порядок составления и решения системы дифференциальных уравнений, описывающих ПП во вторичных цепях ТТ. В дальнейшем укажем подход к математическому описанию ТТ, работающих в группе. Математические модели силовых трансформаторов (например, для изучения броска тока намагничивания) или трансформаторов напряжения (например, для изучения явления феррорезонансных перенапряжений)могут быть составленыпо аналогии.
Рассмотрим схему ТТ, приведённую на рис. 1. В простейшем исполнении ТТ представляет собой две обмотки, намотанные на магнитный шихтованных сердечник, выполненный из электротехнической стали с площадью поперечного сечения S и длиной средней магнитной линии l ср. Предполагается, что материал, из которого выполнен сердечник, известен, то есть известна его основная характеристика намагничивания, представляющая собой зависимость магнитной индукции в стали сердечник от напряжённости магнитного поля B = f (H). Использование основной кривой намагничивания без учёта динамической петли гистерезиса является один из ключевых допущений при моделировании переходных процессов в ТТ, так как он позволяет значительно упростить математическое описание моделируемого объекта.Учёт влияния явления гистерезиса при моделировании процессов, при которых имеет место глубокое насыщение сердечников, не выполняется, а используется основная кривая намагничивания: при рассматриваемых глубоких насыщениях сердечника динамическая петля гистерезиса столь узка, что ею можно пренебречь. При необходимости гистерезис может быть учтён, однако этот случай здесь не рассматривается.
Сам сердечник и обмотки изолированы друг от друга и от земли и находятся внутри заполненного маслом бака или могут быть залиты в эпоксидном компаунде (пунктирная линия на рис.1), что исключает непосредственный доступ к ним.
Выводы первичной обмотки, имеющей число витков W 1и включаемой последовательно в цепь, где предполагается измерение тока, маркируются как «Л1» и «Л2» (название этих клемм образовано от слова «линия»). Выводы вторичной обмотки, имеющей число витков W 2и к которой подключаются измерительные приборы, средства учёта, защиты и автоматики и т.п., маркируются как «И1» и «И2» (название этих клемм образовано от слова «измерение»). Правило маркировки выводов обмоток следующее: если ток в первичной обмотке i 1в определённый момент времени t втекает в клемму «Л1» и вытекает из клеммы «Л2», то в этот же самый момент ток во вторично обмотке i 2должен вытекать из клеммы «И1» и втекать в клемму «И2». Клеммы на рис.1 обозначены в соответствии с данным правилом. Нетрудно видеть, что зажимы «Л1» («Л2») и «И1» («И2») являются одноимёнными.
Число витков вторичной обмотки W 2, как правило, велико, поэтому необходим учёт её активного сопротивления R обм и индуктивности рассеяния L обм. Ёмкостные токи между элементами обмоток (витки и катушки) и между обмотками и сердечником в обычных условиях работы (при частотах ниже 1000÷5000 Гц) весьма малы, и ими можно пренебречь.
Рисунок 1 – Схема трансформатора тока
Нагрузка вторичных цепей ТТ (в расчётную нагрузку включается также сопротивление контрольных кабелей), как правило, активно-индуктивная и может быть представлена параметрами R нагр и L нагр.Так как активные и индуктивные сопротивления обмоток и нагрузки соединены последовательно, их можно объединить эквивалентными суммарными значениями R Σ и L Σ:
(1)
(2)
Обычно при рассмотрении высоковольтных ТТ считают, что ток i 1не зависит от нагрузки, подключенной ко вторичной обмотке ТТ, так как её сопротивление мало и ТТ работает в режиме, близком к режиму КЗ, в результате чего падение напряжения между зажимами «Л1» и «Л2» незначительно. Это позволяет не рассматривать процессы в первичной обмотке, полагая закон изменения тока i 1 известным и заданным в виде уравнения i 1 = f (t).Тогда неизвестными величинами, изменяющимися во времени и подлежащими определению, являются ток во вторичной обмотке ТТ i 2и магнитная индукция в стали сердечника B. Составим и решим дифференциальные уравнения, описывающие изменения во времени двух указанных параметров – ток i 2 и магнитную индукцию B.
Обойдём контур, образованный вторичной обмоткой ТТ, по 2-му закону Кирхгофа, направление обхода контура выберем совпадающим с выбранным направлением тока i 2(рис.1) и учтём, что при выбранном направлении тока i 2и магнитного потока в сердечнике Фэлектродвижущая сила индукции в сердечнике направлена так, что способствует протеканию тока в выбранном направлении, а потому слагаемое, учитывающее магнитную связь, входит в баланс падений напряжения с отрицательным знаком:
(3)
Выразим суммарный магнитный поток через магнитную индукцию в сердечнике B и перепишем вышеприведённое уравнение:
(4)
Полученное уравнение содержит две неизвестные переменные (i 2и B), стоящиепод дифференциалом. Чтобы найти закон изменения двух величин во времени, необходимо составить ещё одно уравнение, связывающее параметры i 2и B между собой. Его можно получить, рассматривая магнитную цепь ТТ. Активными потерями в стали сердечника допустимо пренебречь ввиду непродолжительного характера исследуемых ПП.
При выбранном направлении токов i 1и i 2 магнитодвижущие силы направлены встречно друг другу, причем магнитодвижущая сила, созданная катушкой с током i 1, имеет то же направление, что и магнитный поток Ф. Применим закон полного тока к контуру, совпадающему со средней магнитной линией длиной l ср, обходя контур по положительному направлению потока Ф:
(5)
Продифференцируем по времени полученное значение H:
(6)
Производная di 1 /dt предполагается известной величиной, поскольку закон изменения тока i 1во времени также известен. Свяжем между собой величины B и Н, которые определяют характеристику намагничивания стали сердечника, для чего перепишем величину dB/dt из правой части уравнения (4) следующим образом:
(7)
гдеμ д = dB/dH –дифференциальная абсолютная магнитная проницаемость, определяемая как производная магнитной индукции B по напряжённости магнитного поля H для текущего значения B в сердечнике. Для нахождения μ д требуется информация о виде характеристики B=f(H).
Уравнения (7) и (4) полностью описывают процессы в электрических и магнитных цепях рассматриваемого ТТ. Теперь необходимо определиться со способом решения этих уравнений.Поиск решения в значительной степени осложняется резкой нелинейностью характеристики намагничивания сердечника ТТ, так как даже при пренебрежении явлением гистерезиса строгое и, главное, точное математическое описание характеристики одновременно в насыщенной и ненасыщенной её части с помощью некоторой универсальной формулы затруднительно. Как видно, из (7), нас интересует даже не сама характеристика намагничивания, а её производная в тот или иной момент времени. Это означает, что даже при относительно удачном способе описания характеристики B=f(H) на всех её участках в виде некоторой общей формулы, погрешность неминуемо возрастёт, когда потребуется вычислить производную этой функции. В этой связи следует отказаться от описания характеристики B=f(H) в виде одной формулы и вместо этого использовать аппроксимацию на небольших участках характеристики B=f(H), так как в этом случае как сама функция B=f(H), так её производная могут быть вычислены с требуемой точностью. Наилучшие результаты получаются при использовании кубической, сплайновой или Эрмитовой интерполяции. Однако отказ от использования общей формулы лишает нас возможности получить строгое математическое описание процессов в цепях ТТ и приводит к необходимости решать систему дифференциальных уравнений численными методами. Чтобы численно решить систему дифференциальных уравнений, все уравнения должны быть приведены к форме Коши.
Приведём к форме Коши уравнение, описывающее изменение тока i 2.В соответствии с преобразованиями, уравнение (7) равно сумме падений напряжения в обмотках ТТ и на нагрузке:
(8)
Перенесём в левую часть слагаемые, имеющие множитель di2/dt, а не имеющие его перенесём вправо:
(9)
Выразим отсюда di2/dt:
(10)
Выражение (10), как видно, получилось достаточно громоздим, а потому судить о правильности его вывода на первый взгляд затруднительно. Однако поиск способа доказать корректность модели является неотъемлемым этапом моделирования, поэтому внимательно проанализируем выражение (10), чтобы установить его непротиворечивость реальным физическим явлениям. Как известно, трансформатор передаёт электрическую энергию только на переменной токе. В самом деле, если ток в первичной обмотке постоянный, то его производная di 1 /dt равна нулю и выражение (10) приобретает форму, которой соответствует уравнение, описывающее переходные процессы с RL -контуре, которым является вторичная обмотка ТТ:
(11)
В знаменателе правой части выражения (11) стоит сумма индуктивности внешней цепи и индуктивности катушки со сталью. Также укажем на то, что множитель, стоящий за di 1 /dt в выражении (10), представляет собой коэффициент взаимной индукции между двумя обмотками, который тем меньше, чем меньше витков на первичной и/или вторичной обмотке, чем меньше сечение магнитопровода и больше его длина. Размерность у этого коэффициента та же, что и у индуктивности, а при умножении на di 1 /dt он приобретает ту же размерность, что и слагаемое i 2 R Σ. Наконец, убедимся в правильности размерности физических величин, стоящих по обе стороны знака равенства. В самом деле,[А/с] = [В/Гн]. Убедившись таким образом в корректности составленного первого выражения, выразим в форме Кошиуравнение для магнитной индукции B. Как видно из (7):
(12)
Подставим полученное нами выражение (10) в выражение (12):
(13)
Выражение для магнитной индукции B получилось достаточно громоздким даже для случая простейшей схемы ТТ.В случае, когда приходится иметь дело с трёхфазной группой ТТ, число уравнений и их громоздкость возрастает настолько, что приводить конечный результат преобразования просто бессмысленно и займёт чрезмерно много места. Получить запись большого числа уравнений в форме Коши и не ошибиться в процессе преобразования можно с использованием компьютерных технологий, о чём сказано далее.
Решение задачи Коши предполагает наличие начальных значений переменных i 2и B. Начальное значение тока i 2можно принять равным нулю или равным приведённому к виткам вторичной обмотки значению первичного тока в начальный момент времени, чтобы имитировать работы ТТ в режиме источника тока. Выбор начального значения магнитной индукции в сердечнике B является более сложной задачей. Задание ненулевого значения B в начальный момент времени является способом учёта остаточной индукции сердечника в момент возникновения моделируемой ситуации. В основном именно этот параметр определяет скорость насыщения стали сердечника ТТ и срыв вторичного тока при ПП, что способно привести к нежелательному или замедленному действию средств РЗиА. Однако задание ненулевого значения B при пренебрежении явлением гистерезиса является, строго говоря, некорректным действием, но тем не менее при моделировании такой шаг допустим и позволяет полноценно исследовать ПП во вторичных цепях ТТ без внесения значимых погрешностей по сравнению с реальной ситуацией вследствие узкости гистерезисной петли в сердечниках реальныхТТ.
Любая модель всегда является некоторым упрощением реального объекта и не описывает всех явлений, присущих исследуемому объекту.Модель лишь воспроизводит только наиболее значимые свойства объекта с точки зрения задачи, решаемой при помощи моделирования, а потому при попытке использовать модель в области задач, для которой она не предназначена, можно получить результаты, расходящиеся с поведением реального объекта. Укажем на ограничения, свойственные разработанной модели ТТ. При составлении модели мы пренебрегли явлением гистерезиса, что допустимо только при рассмотрении таких процесс в ТТ, когда магнитная индукция при перемагничивании достигает значений, близких к индукции насыщения, что свойственно токам аварийного режима. Однако в условиях нормального режима ТТ работает на начальном участке характеристики намагничивания для обеспечения требуемой точности измерений. Данному режиму свойственны малые значения B, когда пренебрегать гистерезисом нельзя. Таким образом, модель не предназначена для исследования поведения ТТ в нормальном режиме работы энергосистемы, когда по обмоткам ТТ протекают токи, сопоставимые с номинальными. Однако в таком случае не происходит насыщенияТТ, и задача исследования ПП во вторичных цепях отпадает, так как поведение РЗиА в этих условиях предсказуемо.
Пример 1. Рассмотрим порядок решения системы дифференциальных уравнений, описывающих ПП в ТТ, на примере ТТ типаТПЛ-10, используя возможности программного пакета MATLAB. Рассматриваемый ТТ имеет конструкцию, схожую с таковой, представленной на рис.1. Параметры ТТ разных конструкций приводятся в [1]. Номинальная нагрузка вторичной обмотки класса точности «Р» Z ном = 0,6 Ом; номинальный первичный ток I ном1 = 400 А; номинальный вторичный ток I ном2 = 5 А;номинальная предельная кратность первичного тока (наибольшее значениекратностипервичноготокапо отношения к номинальному его значению, при котором полная погрешность при нагрузке, равной Z ном, не превышает 10%) K 10ном = 13; площадь поперечного сечения сердечника S = 11,2·10-4 м2;длина средней магнитной линии l ср = 0,48 м; число витков вторичной обмотки W 2 = 159; активное сопротивление вторичной обмотки R обм = 0,28 Ом; индуктивное сопротивление вторичной обмотки X обм = 0,1 Ом. Число витков первичной обмотки равно W 1 = 2 (отношение W 2/ W 1несколько меньше отношения I ном1/ I ном2 вследствие витковой коррекции, выполненной с целью уменьшения погрешности измерения). Для решения уравнений ПП в ТТ (а именно – для определения μ д) необходима информация о значениях B для любого значения Н. В таблице 1 приведены значения функции B=f(H) в нескольких точках для электротехнической стали марки Э330А, на рис.2 приведена кривая B=f(H), для её построения использовалась интерполяция кубическими полиномами (табличные значения точек B(H) указаны на графике в виде кружков).Используя большее количество экспериментальных точек B(H), можно при необходимости повысить точность интерполяции.
Таблица 1 – Характеристики основной кривой намагничивания электротехнической стали Э330А
| Магнитная индукция В, Тл | 0,0 | 0,4 | 0,8 | 1,1 | 1,4 | 1,5 | 1,6 |
| Напряжённость магнитного поля Н, А/м |
Продолжение таблицы 1
| Магнитная индукция В, Тл | 1,7 | 1,85 | 1,9 | 1,95 | 2,0 |
| Напряжённость магнитного поля Н, А/м |
Рисунок 2 – Основная кривая намагничивания стали Э330А, полученная в результате интерполяции: а) на интервале изменения напряжённости магнитного поля от 0 до 1000 А/м; б) на интервале изменения напряжённости магнитного поля от 0 до 30000 А/м
Ниже приведён программный код решения с комментариями.
Тело функции example.m:
% Команда очистки командной строки, удаления текущих переменных % и закрытия созданных программных окон соответственно:
clc,clear, closeall
% Объявление глобальных переменных:
global S lcp W1 W2 R_sumL_sum di1 t2 BB HH
%% ПОДГОТОВКА ИСХОДНЫХ ДАННЫХ
fs=4800; % Шаг дискретизации, Гц
t=0:1/fs:0.08; % Массив времени моделирования, с
alpha=-90; %Начальный угол синусоидального сигнала первичного тока, градусы
% Массив мгновенных значений первичного тока, А (имитируется ток % аварийного режима заданной кратности с максимальной % апериодической слагающей):
i1=10*400*sqrt(2)*(sin(2*pi*50*t+alpha*pi/180)+1*exp(-t/0.1));
% Параметры моделируемого ТТ и его нагрузки:
S=11.2e-4; % Площадь поперечного сечения сердечника, кв.м.
lcp=0.48; % Длина средней магнитной линии, м
W1=2; % Число витков первичной обмотки
W2=159; % Число витков вторичной обмотки
Robm=0.28; % Активное сопротивление вторичной обмотки, Ом
Lobm=0.1/314; % Индуктивность вторичной обмотки, Гн
Rload=0.55; % Активное сопротивление нагрузки, Ом
Lload=0.05/314; % Индуктивность нагрузки, Гн
%Определение характеристики намагничивания стали сердечника:
B=[0 0.4 0.8 1.1 1.4 1.5 1.6 1.7 1.85 1.9 1.95 2 3.5];
H=[0 28 56 75 100 130 200 800 2500 5000 10000 30000 1161831];
% Задание параметров точности решения дифференциальных % уравнений. При рассмотрениидлительных ПП может быть% целесообразно задать более высокие параметры точности% (например, установить параметр 'AbsTol' равным 1e-9, а % параметр 'RelTol'установить равным 1e-6), чтобы избежать % накопления ошибок причисленном интегрировании, однако % следует учесть, что время ожидания решения уравнений в этом % случае возрастает.
options=odeset('AbsTol',1e-6,'RelTol', 1e-3);
% Начальные значения переменных 'i2' (A) и 'B' (Тл) % соответственно:
y0=[i1(1)*W1/W2 0];
%% РЕШЕНИЕ УРАВНЕНИЙ ПЕРЕХОДНЫХ ПРОЦЕССОВ
% Взятие производной первичного тока, А/с:
di1=diff(i1)*fs;
% Создание массива времени той же размерности, что и массив % значений 'di1',и в том же временном интервале, что имеет % массив 't':
t2=linspace(t(1),t(end),length(di1));
% Суммарное активное сопротивление и суммарная индуктивность % вторичной обмотки:
R_sum=Robm+Rload;
L_sum=Lobm+Lload;
% Поскольку при перемагничивании значения 'В' и 'Н' могут % приобретать как положительные, так и отрицательные значения, % следует определить характеристику намагничивания при % отрицательных значениях 'В' и 'Н', что легко сделать в силу её % симметричности, пользуясь функцией 'fliplr', меняющей порядок % следования элементов массива на противоположный:
BB=[-fliplr(B(2:end)), B];
% Учитывать первый элемент массива 'В' при его расширении не % будем, так как он равен нулю, аналогично сделаем с массивом % 'Н'.
HH=[-fliplr(H(2:end)), H];
% Построим характеристику намагничивания, пользуясь % интерполяцией кубическими полиномами.
% Значения 'В', в которых мы хотим доопределить функцию 'Н':
B_=-B(end):0.05:B(end);
% Воспользуемсяфункцией'interp1'.На вход функции подаются 3 % массива данных: значения Х и Y исходной функции, также массив % Х1, содержащий точки, в которых надо определить промежуточные % значения функции Y(X).Навыходе функции мы получаем массив% значений функции Y1 в точках Х1. Такжедополнительно % указывается способ интерполяции, здесь – используетсяпараметр % 'cubic'(интерполяция кубическими полиномами).
H_=interp1(BB,HH,B_,'cubic');
% Построим график функции 'В(Н)':
figure, plot(H_,B_),grid, title ('Основная кривая намагничивания стали Э330А'), xlabel ('Напряжённость магнитного поля Н, А/м'), ylabel ('Магнитная индукция В, Тл')
% Очень важно проанализировать вид кривой намагничивания, % полученный после интерполяции, прежде чем приступать к % следующему этапу. В результате интерполяции не должно % получиться никаких резких изломов, выпуклостей или иных % очевидных искажений исходной кривой, так как в противном % случае результат вычисления ПП также будет некорректным. Если % же кривая намагничивания после интерполяции искажена и по % форме не похожа на настоящую, то следует или применить иной % способ интерполяции,или задать большее количество исходных % точек, или и то и другое сразу. В дальнейшем при вычислении % дифференциальной абсолютной магнитнойпроницаемости следует % использовать именно тот способ интерполяции,результаты точной % работы которого не вызывают сомнений.
% Начальная и конечная точки интегрирования, с:
ts=[t(1) t(end)];
% Используем стандартную функцию решения дифференциальных % уравнений'ode23t'. Результатом работы данной функции является % массив данных, содержащий значения времени 'T' и значения % интегрируемых функций 'y'. Входными параметрами функции % 'ode23t' являются: ссылка на функцию, содержащую запись % дифференциальных уравнений в форме Коши (в нашем случае это % функция 'trans_eq', которая описана далее),пределы % интегрирования 'ts', начальные значения переменных 'y0' и % прочие параметры (переменная 'options'),определяющие работу % функции:
[T,y]=ode23t('trans_eq',ts,y0,options);
%% ПОСТРОЕНИЕ ГРАФИКОВ
% Первичный ток, приведённый к числу витков вторичной обмотки, и % вторичныйток:
figure, plot (t,i1*W1/W2, T, y(:,1)), grid, title ('Осциллограммы токов'), xlabel ('Время, с'), ylabel ('Токи, А')
% Магнитная индукция в сердечнике:
figure, plot (T,y(:,2)), grid, title ('Зависимость B(t) в сердечнике'), xlabel ('Время, с'), ylabel ('Магнитная индукция, Тл')
%% ЗАПИСЬ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ В ФАЙЛ
% Создание данной программы является подготовительным этапом для % исследования поведения устройств РЗА, которые будут некоторым % образом обрабатывать сигналы вторичного тока ТТ. Однако % следует отметить, что результаты моделирования, полученные на % выходе функции ode23t() непосредственно не пригодны для % исследования поведения цифровых устройств РЗА, и связано это % со следующим. После решения системы дифференциальных уравнений % значения вторичного тока распределены по времени % неравномерно, что связано с особенностью применённой функции-% решателя, которая изменяет шаг численного интегрирования в % зависимости от точности решения, так какэто позволяет % значительно уменьшить время решения. Однако при исследовании % устройств РЗА следует помнить, что такие устройства сначала % производят дискретизацию входного сигнала с некоторой частотой % дискретизации fs, которая может отличаться в зависимости от % типа устройства. Только после данной процедуры следует% цифровая обработка электрических сигналов.В настоящее время в % цифровых устройствах РЗА частота дискретизации принята равной % 1200 Гц, в дальнейшем планируется её увеличение до 4800 Гц. % Чтобы имитировать работу цифровой защиты, необходимо% исследуемый сигнал привести к равномерной сетки времени с % заданным параметром частоты дискретизации. Выполним указанное % преобразование и сохраним результаты моделирования в отдельном % файле, который можно открыть как с помощью встроенных функций % MATLAB, так и в обычном текстовом редакторе. В дальнейшем при % имитации работы устройств РЗА будем загружать результаты % моделирования в новую программу из этого файла, чтобы % разделить длительные по времени процедуры моделирования ТТ и % имитации работы РЗА.
% Создание равномерной сетки времени c новым шагом дискретизации % fs:
fs=1200; % Гц
X=linspace(0,t(end),t(end)*fs);
% Для определения значений электрических сигналов в точках, % заданных переменной Х, воспользуемся сплайновой % интерполяцией:
Y1=interp1(t,i1*W1/W2,X,'spline'); % Приведение первичного тока
Y2=interp1(T,y(:,1),X,'spline'); % Приведение вторичного тока
% Запись результатов приведения в файл с указанными именем. Файл % будет сохранён в текущй папке MATLAB%
dlmwrite('record.mat',[Y1;Y2]);
% Чтобы в дальнейшем считать этот файл, следует воспользоватьсяфункцией dlmread().
Функция «example» обращается к функции «trans_eq», описание которой приводится ниже. Эта функция должна помещаться в ту же директорию, что и исходная функция «example» или же в иной директории, находящейся в пути поиска среды MATLAB.Тело функции trans_eq.m:
functiondy=trans_eq(T,y)
% Приём переменных из функции example.m:
global S lcp W1 W2 R_sumL_sum di1 t2 BB HH
% Кратко поясним работу данной функции. Эта пользовательская % функция располагается в отдельном файле с именем 'trans_eq', и % к ней обращается функция-решатель 'ode23t'в процессе % численного решения дифференциальных уравнений. Функция-% решатель посылает в функцию 'trans_eq' текущие значения % времени интегрирования 'Т', а также значения переменных 'y' (в % нашем случае -- это значения тока 'i2' и индукции 'В') в % данный момент времени 'Т'. Решатель будет последовательно % обращаться к этой функции до тех пор, пока текущее значение % времени 'Т' не достигнет выставленного ранее конечного времени % интегрирования. Шаг, с которым меняется текущее значение % времени 'Т', заранее неизвестен, так как функция-решатель сама % определяет оптимальное значение шага, чтобы, с одной стороны, % решить уравнения с требуемой точностью, с другой -- чтобы % ускорить время решения. Выходным параметром функции'trans_eq' % являются уравнения ПП в форме Коши, записанные в переменную-% массив 'dy'. При записи уравнений чрезвычайно важно правильно % соблюсти последовательность обозначений переменных. Ранее в % переменной в 'у0' были записаны начальные значения тока 'i2' и % индукции 'В', причём первым элементом мы условились обозначить % ток 'i2', а вторым элементом -- индукцию 'В'. Следовательно, в % этом же порядке следует в данной функции написать уравнения % для приращений тока и индукции: первым элементом переменной-% массива 'dy' должно бытьуравнение приращения тока 'i2', % вторым -- уравнение приращения индукции'B'. Аналогично, % первым элементом переменной-массива 'у' будет являться текущее % значение тока 'i2', вторым -- текущее значение индукции 'B'. % Ясно, что с течением времени будут меняться значения % производной первичного тока и дифференциальная абсолютная% магнитная проницаемость. Обе эти величины входятв % дифференциальные уравнения, подлежащие решению. Текущие% значения этих величин для момента времени 'T' следует % рассчитывать в теле функции 'trans_eq', прежде чем подставить % их в выходную переменную 'dy'.
% Текущее значение производной первичного тока найдём с помощью
% интерполяции кубическими полиномами:
di1_T=interp1(t2,di1,T,'cubic');
% Определим значение 'dB/dH' для текущего значения переменной % 'В', для чего численно возьмём производную функции 'В(Н)' c% фиксированным шагом 'deltaB':
deltaB=0.0002; % Шаг взятия производной, Тл
dH=interp1(BB,HH,(y(2)+deltaB),'cubic')-interp1(BB,HH,(y(2)-deltaB),'cubic');
mu_dif=(y(2)+deltaB)/dH-(y(2)-deltaB)/dH;
% Записываем правые части диф. уравнений в форме Коши:
di2=((di1_T*mu_dif*S*W1*W2/lcp)-y(1)*R_sum)/(L_sum+(mu_dif*S*W2^2)/lcp);
dB=mu_dif*(di1_T*W1-di2*W2)/lcp;
%Вносим в правильном порядке в выходную переменную 'dy' правые % части диф.уравнений:
dy=[di2; dB];
end
Финальный этап – проверка корректности воспроизведения ПП с помощью составленной программы. Лучший способ проверки – сравнение осциллограмм ПП, полученных в результате моделирования, с осциллограммами ПП, полученными на реальном объекте в тех же условиях. Осциллограммы токов в обмотках и зависимость изменения магнитной индукции во времени в насыщенном ТТ приводятся в литературе[1, 2], и эти зависимости совпадают с результатами работы программы, представленными на рис.3. На рис.3, а штриховой линией показан первичной ток КЗ, приведённый к вторичной обмотке ТТ, сплошной линией – ток во вторичной обмотке ТТ. Как видно из опыта, насыщение ТТ происходит в первом полупериоде ПП при нулевом значении остаточной индукции до КЗ. Если бы значение остаточной индукции было больше нуля, насыщение сердечника при прочих равных условиях произошло бы быстрее; если бы остаточная индукция имела отрицательное значение – насыщение произошло бы позже, а в некоторых случаях насыщение могло вовсе не произойти, и ТТ работал бы без существенных погрешностей в течение всего времени ПП, то есть до окончательного затухания апериодической слагающей. Большую роль в скорости насыщения сердечника играет также нагрузка, подключенная к вторичным цепям ТТ. Провести оценку влияния указанных факторов на скорость насыщения сердечника ТТ в ПП полезно самостоятельно.
Отметим одну важную особенность. Выше в таблице 1 по опытным данным мы определили значения кривой B(H) до 2 Тл включительно. Однако как можно видеть из текста программы «example.m», мы также определили дополнительно точку со значением 3,5 Тл, которая продлевает кривую B(H) по линейному закону в области насыщения до указанного значения индукции. Дело в том, что для интерполяции значения B(H) был выбран параметр «cubic» в функции «interp1». Интерполяция кубическими полиномами дала наилучшую форму зависимости B(H), однако используя этот параметр, нельзя проводить экстраполирование зависимости B(H), если расчётное значение B оказывается больше 2 Тл. Попытка в таком случае решить уравнения ПП приведёт к сбою программы-расчётчика и расхождению решения. В случае использования интерполяции сплайнами (то есть если применить параметр «spline» в функции «interp1») эта проблема исчезает, однако сплайновая интерполяция для решения этой задачи оказывается не слишком точной.
Рассмотрим важный инструмент программной среды MATLAB, позволяющий легко приводить систему дифференциальных уравнений к форме Коши, что необходимо для их дальнейшего решению любыми численными методами. Как можно было видеть, приведение к нужному виду уравнений даже для случая одиночно работающего ТТ оказалось непростой задачей, которая становится на порядке сложнее в случае рассмотрения более сложных объектов. Однако символьный процессор программной среды MATLAB способен свести к минимуму сложности, связанные с преобразованием исходных уравнений.
Рисунок 3 – Изменение токов обмотках ТТ и магнитной индукции в его сердечнике в режиме КЗ в условиях наличия в токе повреждения максимальной апериодической слагающей
Отметим, что возможности символьного процессора для производства расчётов в общем виде среды MATLABвесьма обширны (например, возможно проводить интегрирование и дифференцирование функций, записанных в буквенной форме, получая ответ также в виде формулы), однако в данной работе мы заострим внимание только на тех его возможностях, которые необходимы при решении поставленных нами задач.
Пример 2. Представим уравнения, описывающие ПП в цепях, в форме Коши, пользуясь символьным процессором MATLAB.Для начала необходимо составить дифференциальные уравнения, полностью описывающие процессы в электрических и магнитных цепях ТТ. В нашем случае это уравнения (4) и (7). Запишем эти уравнения в матричной форме следующим образом: оставим в левой части значения di2/dt, dB/dt икоэффициенты перед ними, остальное перенесём в правую часть. Тогда получим:
Из этой записи видно, что уравнения будут приведены к форме Коши, если записанная матрица будет приведена к диагональному виду. Таким образом, можно использовать стандартную имеющуюся в пакете MATLAB функцию решения систем линейных алгебраических уравнений, сводящий поиск решения к приведению матрицы к диагональному виду (этот метод решения известен как метод Жордана-Гаусса). Ниже приведён код программы, осуществляющей преобразование исходных уравнений:
% Объявляем переменные, содержащиеся в исходных уравнениях, в % буквенном (символьном) виде. Последовательно введём % переменные, отображающие суммарную индуктивность вторичной % обмотки 'L_sum', суммарное активное сопротивление 'R_sum', % дифференциальную абсолютную магнитную проницаемость'mu_dif', % число витков первичной'W1'и вторичной обмотки 'W2', сечение % магнитопровода'S', длину средней магнитной линии'lcp', % значение производной первичного тока'di1dt', ток во вторичной % обмотке'i2':
symsL_sumR_summu_difW1W2Slcpdi1dti2
% Запишем коэффициенты перед di2/dt и dB/dt в матричной форме:
A=[L_sum, -S*W2;
mu_dif*W2/lcp, 1 ];
%Аналогичным образом записываем правые части:
Z=[ -i2*R_sum;
di1dt*mu_dif*W1/lcp ];
% Вывод решенияна командную строку:
Q=A\Z
В результате работы программы получаем следующее:
Q =
-(R_sum*i2*lcp-S*W1*W2*di1dt*mu_dif)/(S*mu_dif*W2^2+L_sum*lcp)
(mu_dif*(L_sum*W1*di1dt+R_sum*W2*i2))/(S*mu_dif*W2^2+L_sum*lcp)
Первое уравнение является правой частью для di2/dt, второе уравнение является правой частью для dB/dt. Полученная запись является корректной, несмотря на её явную несхожесть с выведенными ранее уравнениями. Чтобы в этом убедиться, подставить полученные правые части дифференциальных уравнений в функцию «trans_eq», не забыв при этом сделать подстановку: вместо переменной i2 следует подставитьпеременную y(2),так как именно в этой переменной хранится текущее значение тока во вторичной обмотке ТТ; аналогично, вместо переменной di1dt следует подставить текущее значение производной первичного тока, которое вычисляется заранее в теле функции. В итоге функция «trans_eq» будет иметь вид (первый поясняющий комментарий для краткости опускаем):
functiondy=trans_eq(T,y)
global S lcp W1 W2 R_sumL_sum di1 t2 BB HH
% Текущее значение производной первичного тока найдём с помощью
% интерполяции:
di1_T=interp1(t2,di1,T,'cubic');
% Определим значение dB/dH для текущего значения переменной В, % для чегочисленно возьмём производную функции В(Н) c % фиксированным шагом deltaB:
deltaB=0.0002;
dH=interp1(BB,HH,(y(2)+deltaB),'cubic')-interp1(BB,HH,(y(2)-deltaB),'cubic');
mu_dif=(y(2)+deltaB)/dH-(y(2)-deltaB)/dH;
% Записываем правые части диф. уравнений в форме Коши:
i2=y(1);% Подстановка
di1dt= di1_T;% Подстановка
% Правая часть для di2/dt:
di2=-(R_sum*i2*lcp - S*W1*W2*di1dt*mu_dif)/(S*mu_dif*W2^2 + L_sum*lcp);
% Правая часть для dB/dt:
dB=(mu_dif*(L_sum*W1*di1dt + R_sum*W2*i2))/(S*mu_dif*W2^2 + L_sum*lcp);
%Вносим в правильном порядке в выходную переменную 'dy' правые части диф.уравнений:
dy=[di2; dB];
end
После этого запускаем функцию «example» и убеждаемся, что работа программы не изменилась.
ЛИТЕРАТУРА
1. Королев, Е. П. Расчеты допустимых нагрузок в токовых цепях релейной защиты / Е. П. Королев, Э. М. Либерзон. – Москва: Энергия, 1980. – 208 с.
2. Электрические цепи с ферромагнитными элементами в релейной защите / А. Д. Дроздов, А. С. Засыпкин, С. Л. Кужеков [и др.]; под ред. В. В. Плато