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

С. М. Тиховод. Разработка схемных моделей метода численного решения уравнений

УДК 519.642.2: 621.3.011.713

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

С. М. Тиховод

Запорожский национальный технический университет, кафедра теоретической и общей электротехники E-mail: [email protected]

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

Ключевые слова: схемные модели, интегродифференциаль-ные уравнения, численные методы.

Scheme Models Development of Integro-Differential Equations Numeral Calculation of Processes Dynamics in Electric Circuits

S. M. Tykhovod

Zaporizhzhya National Technical University,

Chair of Theoretical and Over-All Electrical Engineering

E-mail: [email protected]

Scheme models of numeral calculation of integral-differential equations, describing transients in electric circuits are developed. It is shown that offered modeling has the best fast-acting concerning to the known calculations.

Key words: circuit diagram models, integro-differential equations, numerical methods.

ВВЕДЕНИЕ

Моделирование переходных процессов в электрических цепях является неотъемлемой частью работы по проектированию электротехнических устройств. Уравнения состояния, которые составляются по законам Кирхгофа в мгновенной форме, являются интегродифференциальными уравнениями. Обычно интегродифференциальные уравнения преобразуют в системы дифференциальных уравнений первого порядка, но при этом порядок системы уравнений возрастает. Реальные исследуемые цепи могут содержать несколько сотен элементов, что приводит к большим системам дифференциальных уравнений при моделировании. Для решения таких систем в настоящее время широко применяются многошаговые методы численного интегрирования дифференциальных уравнений с использованием полиномиальной аппроксимации решения. В результате получены разностные схемы [1, гл. 11], позволяющие вычислить значение искомой функции в некоторой временной точке по известным значениям функции (или ее производных) в нескольких предыдущих точках. Однако увеличение порядка систем дифференциальных уравнений приводит к увеличению времени и снижению точности моделирования. Поэтому модификация методов численного решения интегродифференциальных уравнений, приводящая к сокращению времени моделирования, является актуальной задачей.

Для инженеров-электриков, которым важен физический смысл математических действий, представляется более наглядным, если какая-либо математическая операция заменяется схемной моделью. Например, расчет переходных процессов операторным методом [2, гл. 10] сопровождается операторной схемой замещения или дифференциальный анализ электрических цепей [3] сопровождается соответствующими схемами замещения. При этом в схеме замещения соответствующий процесс должен быть такой, чтобы он полностью описывался алгебраическими уравнениями. Схемная модель позволяет от электрической цепи, в которой процессы описываются интегродифференциальными уравнениями, перейти к цепи с изображениями токов, для которых справедливы законы Кирхгофа, приводящие к алгебраическим уравнениям. Это открывает возможность использования всего многообразного аппарата теории цепей для работы с изображениями токов. Поэтому модификация численного метода, сопровождающаяся созданием адекватной схемной модели «дружественной» для инженеров-электриков, является ценной.

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

МАТЕРИАЛ И РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ

Рассмотрим одноконтурную цепь переменного тока, содержащую резистивный (Л), индуктивный (£) и емкостный (С) элементы, включенные последовательно. Пусть до коммутации конденсатор был заряжен до напряжения ис(0). При подключении при I = 0 источника переменной ЭДС е(£) в цепи происходит переходный процесс, который описывается линейным интегродифференциальным уравнением с постоянными коэффициентами:

+ Я1 + -1/ ¿(¿)еЙ + ис(0) = е(1). ш С Уд

(1)

Будем искать решение во временной области, состоящей из N одинаковых шагов Н. Узловые точки, являющиеся границами шагов, обозначим ¿0,¿1,^2, • • • ¿ы•

Решение для тока как функцию от времени в интервале времени [¿0, ¿ы] аппроксимируем полиномом N-й степени:

¿(¿) ~ р(Ъ) = а0 + а1 Ъ + а2¿2 +-----Ь аы¿ы• (2)

Для аппроксимирующего полинома (2) зададим условие, что в точках ¿к деления интервала изменения аргумента

¿(¿к ) = р(Ък) (3)

для к = 0,1, 2,..., N. Если условие (3) записать для каждой точки ¿к, то получим систему линейных алгебраических уравнений, если принять, что Ъ0 = 0 :

«о = ¿(¿о) = ¿0,

а0 + а1Н + а2 Н2 + ■ ■ ■ + аы Ны = ¿(¿1),

ао + а1(2Н) + «2(2Н)2 + ■ ■ ■ + аы (2Н)ы = ¿(¿2),

(4)

^0 + а^Н) + а2^Н)2 +-----Ь аы^Н)ы = ¿(¿ы)•

Вычтем из уравнений системы (4) первое уравнение и получим сокращенную систему:

а1Н + а2 Н2 + ■ ■ ■ + аы Ны = ¿1 (¿1) — ¿0, «1(2Н) + а2(2Н)2 + ■ ■ ■ + аы (2Н)ы = ¿1^2) — ¿0,

(5)

ка1(Ж) + а2(NН)2 + ■ ■ ■ + аы (NН)ы = ¿ы (¿ы — ¿0 )• В матричной форме система (5) имеет вид

Н Н2 ■ ■ Ны

2Н (2Н)2 ■ ■ (2Н)ы

(NН)2 ■ ■ (NН)ы

а1 ¿(Н) — ¿0

а2 = ¿(2Н) — ¿0

аы

или

V ■ А = I

= I — ¿0,

(6)

(7)

где V — матрица Вандермонда без первой строки и первого столбца; А = [а1? а2, • • •, аы]т — вектор коэффициентов аппроксимирующего полинома; I = [¿(Н), ¿(2Н), • • •, ¿^Н)]т — вектор значений тока в опорных точках 1, 2, •••^. Будем считать, что номер k отрезка, на которые разделен интервал изменения аргумента, совпадает с номером точки деления ¿к, расположенной справа отрезка.

Продифференцируем выражение (2):

Ш

= а1 + 2а2£ +-----Ь Naы¿(ы

(8)

С. М. Тиховод. Разработка схемных моделей метода численного решения уравнений

Если для точек ¿1 , £2,...,£м подставить в (8) значение времени выраженное через шаг и номер точки, то получим следующую систему линейных уравнений:

ai + +-----h NonhN-1 = ¿i (ti),

ai + 2a2(2h) + ■ ■ ■ + NaN (2h)N-1 = ¿i(t2 ),

(9)

а + + ■ ■ ■ + (Ж)м-1 = ¿М (¿м).

В матричной форме система (9) имеет вид

" 1 2h ■ ■ NhN-i " ai

¿2 = 1 2 ■ 2h ■ ■ N(2h)N-i a2

JN- _ 1 2Nh ■ ■ N(Nh)N-i_ aN

(10)

или

I' = Т1 ■ А, (11)

где I' = [¿'(Л,), ¿'(2Д),..., ¿'— вектор производных тока для точек с номером к = 1, 2,..., N. Проинтегрируем выражение (2) от нуля до к-й точки при изменении номера к от 1 до N:

rtk

p(t)dt =

rtk

ao + ait + a2t* + ■ ■ ■ + aN ^ dt = aotk + y tk2 + y tk3 + ■ ■ ■ + N+1 tkN+1.

2N

ao + ait + a2

oo

Подстановка в (12) значений к от 1 до N дает:

' rti

(12)

o

rto

p(t)dt = aoh + ai h2 + | h3 + ■ ■ ■ + N+y hN+ = Ji,

p(t)dt = 2aoh + f- (2h)2 + f (2h)3 + ■ ■ ■ + (2h)N+

= J2,

(13)

/o p(t)dt = Naoh + -2(Nh)2 + ^(Nh)3 + ■ ■ ■ + N^W

= J

N.

Если учесть, что а0 = ¿(¿о) = ¿о, то из (13) получим следующую систему в матричной форме:

1 Лм+1

N+1 Л

+1 (2Л.)М+1

Ji г 2 h2 55" со

J2 = i (2h)2 i (2h)3

JN j (Nh)2 i (Nh)3

N +i

+Г (Nh)

N+i

N+i

ai h

a2 2h

= ¿o ■

aN .Nh.

или

(14)

(15)

3 = Т2 ■ А + ¿оН,

где I — вектор интегралов (14) для значений к = 1, 2,..., N.

Распишем уравнение (1) с учетом аппроксимации (2) для точек к = 1, 2, ...^. Получим в матричной форме выражение:

Ы' + Ш + Ш + исо = е, (16)

где В = 1/С, ис0 — значение напряжения на конденсаторе в начале интервала, е — вектор значений ЭДС источника в точках 1, 2,..., N текущего временного интервала.

Если подставить в выражение (16) матрицы I (7) , I' (11), 3 (15) , то получим:

(LTx + RV + BT2)A = e - uco - Rio - BHio.

(17)

Уравнение (17) можно интерпретировать следующим образом. Пусть в контуре исходной цепи R-L-C-e протекает ток ¿(¿). Тогда согласно уравнению (17) исходной цепи соответствует цепь замещения, в которой протекает сигнал А, изображающий ток ¿(¿). При этом в цепи замещения резистивный

o

t

N

элемент имеет операторное сопротивление КУ и последовательно с ним навстречу току включается источник ЭДС Ля0 (рис. 1).

Индуктивный элемент имеет операторное сопротивление ¿Т, а емкостный элемент имеет операторное сопротивление ВТ2 и последовательно с ним навстречу току включается источник ЭДС ВНЯ0 + ис0. Источнику ЭДС е(£) в оригинальной цепи соответствует изображающий вектор е в цепи замещения.

Докажем, что в узлах схемы замещения для изображений А соблюдается закон токов Кирхгофа. Для этого воспользуемся уравнением (7). В любом узле электрической цепи для токов ветвей, принадлежащих узлу, в любой момент времени, а следовательно, токов в начале интервала Я0 и для векторов токов I, выполняется закон токов Кирхгофа:

Рис. 1. Схема замещения контура R-L-C-e

¿(Ik - iok) = ^ VAk = 0,

(18)

k = 1

k=1

где Ь — количество ветвей, сходящихся к узлу, к — текущий номер ветви, сходящейся к узлу. Если уравнение (18) умножить на матрицу, обратную матрице У, то получим:

¿(Ak) = 0.

(19)

k=1

Что и требовалось доказать.

Из изложенного материала сделаем выводы.

Реальному току i(t) соответствует векторное изображение А в схеме замещения, показанной на рис. 1. Все изображения тока А в схеме замещения удовлетворяют законам Кирхгофа, если схема замещения составляется по следующему правилу:

- источник ЭДС заменяется векторным источником e, содержащим значения ЭДС в N опорных точках;

- резистивный элемент имеет операторное сопротивление RV и последовательно с ним навстречу току включается дополнительный источник ЭДС Ri0;

- индуктивный элемент имеет операторное сопротивление LT1;

- емкостный элемент имеет операторное сопротивление BT2 и последовательно с ним навстречу току включается дополнительный источник ЭДС BHi0 + uc0.

Таким образом, в схеме замещения электрической цепи изображения Ak оригиналов токов ik(t) удовлетворяют законам Кирхгофа. Следовательно, при известных значениях токов ветвей i0k и напряжений на конденсаторах uc0k в начале интервала [t0,tN] система уравнений, составленная по законам Кирхгофа, для всех узлов без одного и для всех главных контуров имеет единственное решение. В результате решения системы линейных алгебраических уравнений получаем векторы полиномиальных коэффициентов Ak для всех ветвей. Зная для любой ветви коэффициенты полинома и значение i0 в начальной точке t0, мы можем получить значение тока и напряжения на конденсаторе во всех произвольных точках любого из N отрезков в интервале времени [t0, ¿n].

Рассмотрим, как ведет себя полиномиальная аппроксимация, если изменять число опорных точек N. Как показано в работе [4, chap. 3] максимальная погрешность R аппроксимации полиномом p(x) степени N +1 некоторой функции f (x), имеющей ограниченные непрерывные производные до степени N +1 на участке [a, b] , при условии, что в N +1 различной точке полином совпадает с функцией f (x), определяется выражением

(b - a)

N+1

R = maX |f (x) - P(x)| -

a<x<b

max |f(N+1)(x)|

a<x<b

(N + 1)!

(20)

С. Л1. Тиховод. Разработка схемных моделей метода численного решения уравнений_

Чтобы практически воспользоваться условием (20) представим, что при разложении в ряд Фурье функция /(х) реально содержит Б гармоник. Тогда согласно теореме С. Берштейна [5, гл. 3] значение производной /(*+1) (х) по модулю имеет ограничение:

|/(*+1)(х)| < Б*+1 М, (21)

где М = тах |/(х)|.

а<х<Ъ

С учетом выражения (21) погрешность (20) при равномерном расположении опорных точек преобразуем к относительной форме:

Л < N *+1 Б N+1 (22)

5 = М < N + 1)! Б . (22)

Расчеты по (22) показывают, что при фиксированных значениях Б относительная погрешность уменьшается с ростом количества опорных точек и степени полинома. Однако трудно заранее предсказать, сколько реальных гармоник содержит искомое решение. В работе [6, гл. 4] авторы предостерегают, что для некоторых функций значения производных порядка N растут быстрее, чем N!. Поэтому рекомендуют практически пользоваться интерполянтами степени не больше пятой. Однако это касается некоторых искусственных функций, которые на практике встречаются редко. На участке [а, Ь] погрешность ведет себя неравномерно. Зависимость погрешности от номера точки k на [а, Ь] при равномерном расположении опорных точек определяется следующим выражением [7, гл. 2]:

к!^ - 1 - к!) ■ тах |/(*+1)(х)|

Л <--+1 • (23)

Если в выражении (23) учесть ограничение (21), то получим относительную погрешность в зависимости от номера шага к на [а, Ь] :

Л < И^ - 1 - к!) ■ Б*+1 Л„+1 (24)

Л < (N +1)! ' • ( )

Оценка (24) показывает, что на краю участка [а, Ь] погрешность интерполяции максимальна. Поэтому, к примеру, для снижения погрешности можно применить полином степени 10, но использовать интерполяцию в девяти точках.

Однако увеличивать степень полинома бесконечно нельзя. С ростом числа N матрицы становятся плохо обусловленными. Матрица Т2 в выражении (14) содержит шаг интегрирования в степени N + 1. Чтобы при малом шаге Н матрицы У, Т1 и особенно Т2 не стали плохо обусловленными необходимо шаг интегрирования умножить на нормирующий коэффициент такой, чтобы нормированный шаг был близок к единице. Тогда на нормирующий коэффициент нужно умножить значения всех ин-дуктивностей и емкостей, входящих в схему, а значение частоты нужно разделить на нормирующий коэффициент. Это улучшит обусловленность матриц, но увеличение степени полинома больше десяти может в некоторых случаях привести к неадекватному решению.

На больших интервалах изменения независимой переменной Ь >> NН уравнение (17) можно решать методом циклической прогонки, увеличивая каждый раз текущее время на N шагов. В результате определим коэффициенты аппроксимирующего полинома для всего интересующего интервала времени. В каждом цикле необходимо вычислять новые значения ¿0 и ис0, которые будут использованы в следующем цикле. Можно эти значения вычислить в конце текущего участка, но для улучшения точности интерполяции можно часть последних шагов не использовать:

¿0 = ¿0пред + а1 (^ - ^)Н) + а2(^ - ^)Н)2 + ■ ■ ■ + а*(^ - ^)Н)*. (25)

ис0 = ис0пред + В[(^ - адн) + у (^ - N)Н)2 + ■ ■ ■ + ^а+у (^ - ^)Н)*+1], (26) где ¿0пред и ис0пред равны значениям ¿0 и ис0, вычисленным в предыдущем цикле.

Рис. 2. Схема модельной цепи

Для испытания разработанной методики в системе МаНаЬ составлено несколько компьютерных программ для расчета переходных процессов в электрических цепях. Рассмотрим расчет переходного процесса некоторой модельной цепи (рис. 2). Требуется рассчитать переходный процесс изменения токов после замыкания ключа.

Согласно правилу преобразования, показанному на рис. 1, составим схему замещения для изображений токов, показанных на исходной схеме (рис. 2). Схема замещения для изображений показана на рис. 3.

На схеме замещения дополнительные источники ЭДС имеют следующие обозначения индексами: первый индекс «0» означает, что берется значение тока в точке к = 0 текущего временного интервала; второй индекс обозначает номер ветви в цепи. Для схемы замещения, показанной на рис. 3, система уравнений, составленная по законам Кирхгофа для изображений токов, имеет следующий вид:

Рис. 3. Схема замещения модельной цепи А0 - А1 - А4 = 0, А1 - А2 - Аз = 0,

ШоУАо + Ш4УА4 + В4Т2А4 = е - Шо^оо - £4Н«о4 - иС04,

Ш1УА1 + ЫТ1А1 + Ш2УА2 - В4Т2А4 - Ш4УА4 = £4Шо4 - Ш1«01 - Ш¿02 + иС04 + Ш^04, -ШУА2 + ШзУАз + £зТ2Аз = ^¿02 - Шз^оз - £зШоз,

(27)

где ¿00, ¿01, ¿02, ¿оз, ¿04, ис0з, ис04 — обозначения токов ветвей 0-4 и напряжений на конденсаторах 3, 4 в начале текущего цикла, Н — вектор шагов (см. выражение (6)).

Система уравнений (27) для изображений токов и напряжений является системой алгебраических уравнений с постоянными коэффициентами, и имеет единственное решение — векторы Аъ, где номер ветви Ъ = 0•• • 4.

По программе, составленной согласно предложенному алгоритму, выполнен расчет переходного процесса при следующих значениях исходных данных: е(£) = 100вт(27г/г + тг/4), / = 50 Гц, Яг = 1 Ом, Д2 = 200 Ом, Л3 = 1 Ом, = 3 Ом, С3 = 10 мкФ, С4 = 50 мкФ. График тока -¡г^.), полученный в результате расчета, представлен на рис. 4.

Для оценки точности вычислений по предложенному методу выполнен также точный аналитический расчет переходного процесса при тех же значениях исходных данных. Полученное аналитическое выражение для тока ¿1 (¿) имеет вид

300

1, гпэ

Рис. 4. Зависимость от времени тока, полученная в результате моделирования (звездочками показаны точки, полученные точным расчетом)

С. M. Тиховод Разработка схемных моделей метода численного решения уравнений_

¿1 (t) = DiePl* + D2ep2t + D3eP3* + /im sin(wt + ф), (28)

где Di = -0.355 - 1.01j, D2 = -0.355 + 1.013j, D3 = 0.137, pi = (-348.51 + 3143j) ■ 103, p2 = (-348.51 - 3143j) ■ 103, рз = -5.548 ■ 103, /im = 0.592A, ф = 1.315.

Шаг интегрирования при расчете предложенным методом выбран таким, чтобы относительное отклонение значения тока от соответствующего точного значения в точках локальных максимумов не превышало ±10% (при уменьшении шага интегрирования погрешность уменьшается). Выполнен также расчет модельной задачи при использовании метода Гира с максимальным постоянным шагом интегрирования таким, чтобы относительное отклонение значения тока от соответствующего точного значения в точках локальных максимумов также не превышало ±10%. С помощью операторов tic/toc оценивалось время расчета. Сравнение процессорного времени расчета модельной задачи по предложенному методу и по методу Гира показало, что рассматриваемый метод имеет быстродействие более чем в четыре раза лучшее, чем многошаговый метод Гира. Это можно объяснить следующими соображениями:

- в предложенном методе уже известна форма решения (полиномоидальная), поэтому достаточно найти решение только в редких опорных точках, а в промежуточных точках оно определяется по коэффициентам полинома;

- за каждый цикл выполнялся расчет не одного, а N = 8 шагов интегрирования;

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

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

Выводы

1. Предложенная методика построения схем замещения электрической цепи позволяет непосредственно вычислять коэффициенты полиномиальной аппроксимации искомых токов в переходных режимах.

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

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

4. Для модельной задачи отмечено лучшее быстродействие расчета по сравнению с известными численными методами.

Библиографический список

1. Современные численные методы решения обыкновенных дифференциальных уравнений / под ред. Дж. Холл, Дж. Уатт. М., 1979. 312 с.

2. Демирчян К. С., Нейман Л. Р., Коровкин Н. В., Че-чурин В. Л. Теоретические основы электротехники : в 2 т. СПб., 2003. Т. 2. 567 с.

3. Пухов Г.Е. Дифференциальный анализ электрических цепей. Киев, 1982. 496 с.

4. Shampine L.F., Allen R.C., Pruess S. Fundamentals

of numerical computing. N.Y.; Chichester; Brisbane; Toronto; Singapore, 1997. 268 c.

5. ЗиГмунд А. Тригонометрические ряды : в 2 т. М., 1965. Т. 2. 538 с.

6. Каханер Д., Моулер К., Нэш С. Численные методы и программное обеспечение : пер. с англ. М., 1998. 575 с.

7. Буслов В. А., Яковлев С. Л. Численные методы. Исследование функций : курс лекций. СПб., 2001. 59 с.