CC BY
44
УДК 621.3.012.8::621.3.011.719:681.3
С. М. Тиховод канд. техн. наук Запорожский национальный технический университет
Моделирование динамических процессов в нелинейных магнитоэлектрических цепях с использованием диакоптики и параллельных вычислений
Разработана методика моделирования динамических процессов в магнитоэлектрических цепях с использованием диакоптического подхода и параллельных вычислений. Приведен пример расчета переходного электромагнитного процесса в катушке с ферромагнитным сердечником, используя параллельные вычисления.
моделирование, электромагнитный процесс, ферромагнитный сердечник, параллельные вычисления
Введение
Разработка сложных электромагнитных устройств, как правило, выполняется при активном использовании компьютерного имитационного моделирования динамических электромагнитных процессов. Для эффективного моделирования указанных процессов используются различные магнитоэлектрические схемы замещения, содержащие электрические, а также магнитные ветви [1]. Схемы замещения таких устройств, как мощные трансформаторы или реакторы, могут быть весьма сложны, поэтому известные системы моделирования PSpice или Simulink становятся малоэффективными. Для успешного моделирования динамических электромагнитных процессов с помощью магнитоэлектрических схем замещения разработан программный ком -плекс Colo [1]. Комплекс разработан в предположении, что исследуемые цепи состоят только из двухполюсников. Однако усложнение конструкций трансформаторов и схем, в которые они включаются, а также большая детализация конструктивных особенностей в моделях данных изделий привели к тому, что размер главной матрицы, используемой в расчете на каждом шаге интегрирования, может достигать нескольких тысяч. Все матрицы имеют тенденцию: чем больше размер матрицы, тем больше число обусловленности матрицы и тем меньше точность расчета системы уравнений. Для устранения данного недостатка в работе [2] предложено использовать диакоптический подход создания подсхем и выполнена модификация программного комплекса Colo для поддержки подсхем. Это позволяет снизить размер используемых матриц, а также упрощает подготовку исходных данных. Однако время выполнения моделирования больших задач может быть весьма значительным. В последних версиях системы Matlab, с помощью которой разработан программный комплекс Colo, добавлена возможность параллельных вычислений. Для эффективного использования таких возможностей не© С. М. Тиховод 2010 р.
обходимы кластеры или компьютеры с многоядерными процессорами. В последнее время такие компьютеры появились в массовом пользовании, и их цена соизмерима с ценой обычных компьютеров. Поэтому использование диакоптических и параллельных вычислений для моделирования динамических электромагнитных процессов в сложных магнитоэлектрических схемах замещения является актуальной задачей.
Целью статьи является разработка методики использования диакоптических подходов и параллельных вычислений в системе Matlab для моделирования динамических электромагнитных процессов в сложных магнитоэлектрических схемах замещения.
Постановка задачи
Поддержку подсхем можно организовать различными способами. Представим себе, что часть схемы мы хотим представить в виде подсхемы (рис. 1).
Рис. 1. Основная схема исследуемой цепи с выделенной подсхемой
В работе [2] показан способ замещения подсхемы эквивалентными источниками тока, число которых равно числу входов подсхемы. В данной статье показано использование другого способа организации подсхем:
1. Подсхему, показанную на рис. 1, заменим источниками ЭДС (рис. 2). Количество источников равно ко-
Рис. 2. Замена подсхемы в основной схеме эквивалентными источниками ЭДС
Рис. 3. Замена основной схемы в подсхеме эквивалентными источниками тока
личеству входов подсхемы.
2. Зададим начальные приближения токов входов подсхемы.
3. Для цепи, показанной на рис. 2, выполним расчет токов входов подсхемы /1, /2, /3.
4. Для цепи, показанной на рис. 3, выполним расчет напряжений входов подсхемы м1, м2, м3 при условии, что значения токов источников тока приняты вычис-
ленным значениям /1, /2, /3.
5. Для цепи, показанной на рис. 2, зададим значения напряжений источников ЭДС, равные вычисленным значениям м 1, м2, м3, и выполним расчет токов входов подсхемы /1, /2, /3.
6. Переход на п. 4.
7. Итерационный цикл завершается при достижении сходимости вычислительного процесса, то есть если значения токов /1, /2, /3 на двух смежных итерациях соответственно отличаются не более чем на заданную малую величину е.
Если моделируемое устройство содержит нелинейные элементы, то пункты 4 и 5 должны содержать внутренние итерационные циклы.
Для иллюстрации предложенной методики и описания практического применения параллельных вычислений рассмотрим практический пример. К электрической цепи, содержащей источник синусоидальной ЭДС, подключена катушка с ферромагнитным сердечником (рис. 4). Требуется промоделировать переходный процесс токов во всех ветвях, а также магнитного потока после подключения источника с помощью ключа.
Катушка с ферромагнитным сердечником может быть представлена в виде модели с помощью магнитоэлектрических схем замещения, составленных различными методами, использующими:
- явный магнитный поток;
- производную магнитного потока и дифференциальные магнитные сопротивления;
- производную магнитного потока и «магнитные» конденсаторы [2].
Остановимся на последнем методе.
В работе [3] введено понятие магнитного тока смещения /ц = ЙФ / Л по аналогии с электрическим током смещения плотностью сСБ / Й/. Подобно электрическому току смещения, протекающему через емкостной элемент, магнитный ток смещения должен протекать через магнитный емкостной элемент С . Рассмотрим
Рис. 5. Магнитоэлектрическая схема замещения катушки с сердечником с помощью магнитного конденсатора
взаимодействующие электрическую и магнитную ветви магнитоэлектрической схемы замещения. Пусть в отрезке магнитной ветви, на которой расположена катушка с числом витков Ы, проходит сигнал, значение которого равно производной магнитного потока. За-
Рис. 6. Магнитоэлектрическая схема замещения, разделенная на две подсхемы
меним магнитную ветвь виртуальным магнитным конденсатором Ст (рис. 5).
Тогда напряжение на магнитном конденсаторе, ум -ноженное на емкость, равно:
umCm =
• dФ dt
dt = Ф.
(1)
Отсюда и = Б/С . Следовательно, емкость магнитного конденсатора равна обратной величине магнитного сопротивления Ст = 1/ Ят.
В магнитную ветвь можно также включать резистивный элемент, который характеризует потери в стали [3]. Таким образом получается компактная схема замещения. С учетом сказанного катушку с ферромагнитным сердечником, показанную на рис. 4, можно представить магнитоэлектрической схемой замещения. Ди-акоптический подход предполагает разделение большой цепи на подсхемы. Разделим моделируемое устройство, показанное на рис. 4, по линии А1-А2. В результате получим две схемы замещения, изображенные на рис. 6. В магнитоэлектрических схемах для единообразия все электрические токи и магнитные токи смещения йФ/ й будем обозначать символом «і».
Составим по законам Кирхгофа систему уравнений для цепи, показанной на рис. 6, а:
— /’1 + i0 — /4 = 0,
uc4 + R0i0 + R4i4 = e(t), di1
L— + R1i1 — uc4 — R4i4 = —ES . dt
C
du
c4
dt
— i4 = 0
В матричной форме система (2) имеет вид:
0 0 -1 0 1 -1
0 0 0 1 * о *4
1 0 -1 0 0? -
0 С4 0 0 0 -1
*1
ис4 1 0 1
є(ґ)
ис4 - Е
/0 0
4
(3)
или
МІ.-XI = Р1.. ,
кіг кіг 7
(37)
Где символом «'» обозначена производная по времени соответствующей переменной.
Чтобы матричное уравнение (3) имело единственное решение, необходимо, чтобы количество строк в матрицах М1к11, Р1к.г равнялось количеству переменных. Для этого систему уравнений (2) необходимо дополнить для каждой переменной состояния уравнением численного метода интегрирования обыкновенных дифференциальных уравнений. В простейшем случае это метод Эйлера, но более эффективными являются многошаговые неявные методы. Для неявного метода Эйлера дополнительные уравнения имеют вид:
с4 Иис4 — ис4(к-1)1
(4)
В матричной форме система (4) имеет вид:
И 1
- И 1
• XI — /1(к-1)
ис4(к-1)
(5)
где к - шаг интегрирования численного метода, а индекс (к-1) означает, что значение переменной соответствует предыдущему шагу интегрирования (номер текущего шага равен к).
Или
М1 -XI = F1
шеі ШЄІ
Объединим уравнения (3) и (5) в одну систему:
" Мк/г ' • XI — ^1к/г
_М теі _ F1me^ _
(6)
Система (6) имеет единственное решение, если на каждом шаге интегрирования известны значения, входящие в правую часть, то есть ЭДС источника в данный момент времени, значения переменных состояния на предыдущем шаге интегрирования и значение ЭДС виртуального источника Е.
Составим теперь по законам Кирхгофа систему уравнений для цепи, показанной на рис. 6, б:
?2 - /3 + /4 — 0 ,
Я2/2 + ист + ^3 + ис3 — 0 ' — ис3 — Я3/3 — N/7 — 0,
/7 + /'о — 3,
- Я6/6 + Я1/1 + N/2 — 0, Яб1б - и — 0
(7)
К уравнениям Кирхгофа добавим уравнения, связывающие напряжение и ток конденсатора С3, а также уравнения, связывающие напряжение ист, магнитный ток /2 и магнитный поток Б магнитного конденсатора С :
- (Стист ) + /2 — 0
- Стист + Ф — 0
- С3ис3 + /3 — 0
(8)
В первом уравнении системы (8) емкость магнитного конденсатора С - величина переменная и поэтому ее нельзя выносить за знак производной.
Объединим системы (7) и (8) в одну систему и запишем ее в матричной форме:
(Стист )
0 0 0 0 0 1 -1 1 0 0 0 " ис3 0
0 0 1 1 0 *2 Я3 0 0 0 0 ист 0
0 0 0 -1 - N 0 Я3 0 0 0 0 ис3 0
-1 0 0 0 0 1 0 0 0 0 0 /7 0
0 - С3 0 0 0 0 1 0 0 0 0 0 /2 — 3
0 0 -С т 0 0 0 0 0 1 0 0 Ч 0
0 0 0 0 1 0 0 0 0 1 0 Ч 0
0 0 0 0 Я7 N 0 0 0 - Я6 0 Ф Ф(к-1)
0 0 0 0 0 0 0 0 0 Яб -1 ‘в и ис3(к-1)
(9)
Или
М2-Х2 = Р2. ,
кіг кіг7
(90
В системах (7), (8) переменными состояния являются напряжение на конденсаторе и и магнитный поток Ф. Для переменных состояния составим аналогично (5) уравнения численного метода:
М2 Х2 = F2 ,.
шеі ШЄІ
(10)
Объединим уравнения (9) и (10) в одну систему:
" М2кіГ' її сч * ^2кгг
_М2~" _ _^теі _
(11)
Система (11) имеет единственное решение, если на каждом шаге интегрирования известны значения, вхо-
дящие в правую часть, то есть значения переменных состояния на предыдущем шаге интегрирования и значение тока виртуального источника 3. Емкость магнитного конденсатора С - величина переменная, зависит от магнитного потока. Поэтому на каждом шаге интегрирования должен быть организован итерационный процесс, в котором значение С вычисляется как предельная величина.
Исходя из изложенного, разработан алгоритм решения данной задачи. Блок-схема алгоритма представлена на рис. 7. В соответствии с численным методом интегрирования обыкновенных дифференциальных урав -нений весь процесс объединен циклом при изменении номера шага интегрирования т. На каждом шаге интегрирования организован итерационный цикл согласо-
Рис. 7. Блок-схема алгоритма решения задачи диакоптическим способом с использованием параллельных вычислений
вания процессов с параметром цикла п. Далее идет разветвление на два процесса. 1аЫМех - системная переменная, задающая номер параллельного процесса. В каждом параллельном процессе вычисления выполняются аналогично. В первом процессе для получения текущего решения на данном шаге интегрирования в правой части системы уравнений используется значение напряжения согласования Екоторое принимается от второго процесса с помощью оператора Е х = 1аЪЯесе1уе(2).
В начальный момент времени второй процесс еще не пошел и значение напряжения согласования Е первому процессу не передано. Поэтому используется начальное значение напряжения согласования Е0. Далее вычисляется вектор текущего решения Х14. Если электрическая цепь нелинейная, то вектор текущего решения вычисляется многократно в итерационном цикле, до тех пор, пока на двух смежных итерациях разность токов на нелинейных элементах не станет меньше заданной малой величины. После завершения итерационного процесса из вектора решения Х1( извлекается значение тока согласования 3 и с помощью оператора 1аЪ8еМ(./,2) его значение пересылается второму процессу.
Во втором процессе оператором 3 = 1аЪЯесе1уе(1) производится прием значения тока согласования, которое вычислено в первом процессе и передано второму процессу. Затем вычисляется правая часть системы уравнений (9) и вектор решения Х2(. Вектор решения многократно вычисляется и корректируется емкость магнитного конденсатора Ст в итерационном процессе. После завершения итерационного процесса из вектора решения Х2( извлекается значение напряжения
согласования E s и с помощью оператора labSend(Es,1) его значение пересылается первому процессу.
Такие вычисления проводятся на каждом шаге интегрирования в цикле согласования, который охватывает оба процесса. Как показывают вычисления, для достижения практической точности достаточно 3-6 циклов согласования.
В системе Matlab параллельные процессы можно запускать различными способами. Рассмотрим технологию pmode [4]. Чтобы запустить два параллельных процесса, в командном окне, введем команду: pmode start 2. Появится окно параллельных процессов (рис. 8).
В командном окне курсор приглашения изменится на P>>. В командном окне введем вызов нашей программы, составленной в соответствии с изложенным алгоритмом. В результате расчетов в каждом окне появятся результаты вычислений. Для дальнейшего использования результатов их необходимо передать из рабочих процессов в клиентскую станцию. В нашем случае это главное рабочее пространство workspace системы Matlab. Например, передача массива решений из второго процесса X2 и вектора времени tm в клиентскую станцию производится командами: pmode lab2client X2 2; pmode lab2client tm 2.
Далее полученные данные можно обрабатывать всеми средствами системы Matlab. Например, график изменения тока 17, протекающего через катушку, можно построить командой:
plot(tm, X2(5,:),’k-’); grid; xlabel(‘t, c’);ylabel(‘i7,A’).
В результате получим график изменения тока, показанный на рис. 9.
Рис. 8. Окно параллельных процессов системы Matlab
О 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
с
Рис. 9. График изменения во времени тока катушки после подключения источника
Разработана аналогичная программа, позволяющая выполнять данный расчет с применением диакоптики, но без аппарата использования параллельных вычислений. Измерение процессорного времени операторами йсЛос при заданной одинаковой погрешности показало, что при использовании параллельных вычислений время расчета уменьшается на 20 %. Такое незначительное уменьшение времени расчета объясняется тем, что матрицы в данной задаче небольшие и время расчета соизмеримо со временем передачи данных от одного рабочего процесса другому. Существенного сокращения времени расчета при использовании параллельных вычислений следует ожидать при решении больших задач на компьютерах с большим количеством ядер.
Выводы
Разработана методика моделирования динамических процессов в магнитоэлектрических цепях с использованием диакоптического подхода и параллельных вычислений. Разработанная методика предназначена для моделирования динамических процессов в цепях, имеющих большое число нелинейных элементов. Она может использоваться на компьютерах с многоядерными процессорами.
Список литературы
1. Тиховод С. М. Система компьютерного моделирования динамических процессов в нелинейных магнитоэлектрических цепях / Тиховод С. М. // Технічна електродинаміка. - 2008. - № 3. - С. 16-23.
2. Тиховод С. М. Разработка аппарата поддержки подсхем в программном комплексе Colo / Тиховод С. М. // Електротехніка та електроенергетика. - 2005.
- № 2. - С. 24-25.
3. Шакиров М. А. Магнитоэлектрические схемы замещения катушек индуктивности и трансформаторов / Шакиров М. А. // Электричество. - 2003. -№ 11. - С. 34-45.
4. Оленев Н. Н. Параллельное программирование в Matlab и его приложения // Оленев Н. Н., Печен-кин Р. В., Черенцов А. М. - М. : Вычислительный центр им. А. А. Доронина, 2007. - 120 с.
Поступила в редакцию 17.12.09 г.
Тиховод С. М. Моделювання динамічних процесів в нелінійних магнітоелектричних колах з використанням діакоптіки і паралельних обчислень.
Розроблена методика моделювання динамічних процесів в магнітоелектричних колах з використанням діакоптичного підходу і паралельних обчислень. Приведений приклад розрахунку перехідного електромагнітного процесу в котушці з феромагнітним осердям з використанням паралельних обчислень.
моделювання, електромагнітний процес, феромагнітне осердя, паралельні обчислення
Tykhovod S. M. Modeling of dynamic processes in nonlinear magnetoelectric circuits by the use diakoptic and parallel culculation.
Modeling of dynamic processes in the magnetoelectric circuits by the use diakoptic procedure and parallel calculation is developed. Instant of parallel calculation of electromagnetic transient in the coil with the ferromagnetic-core is adduced.
modeling, electromagnetic process, ferromagnetic core, parallel calculations
