О выборе численных методов интегрирования уравнений переходных процессов электроэнергетических систем

U_d = –ψ_q·ω₀
U_q = ψ_d·ω₀ (2)
Основанные на (2) уравнения электромагнитных переходных процессов, которые применяются при обычных проектных и эксплуатационных расчетах устойчивости, получается в следующем виде:
при r = 0
U_q = E_q’ – x_d’·i_d (3)

где		–	переходным реактивным СГ;
		–	ЭДС за переходным реактивным.

Из векторной диаграммы СГ следуют следующие связи:
U_q = E_q - x_d·i_d
U_d = x_q·i_q
тогда получается
E_q = E_q’ + (x_d – x_d’)·i_d (4)
тоже на основании (2) можно получить
(5)

где	Eq = ψf = xad·if	–	ЭДС холостого хода;
	Eqe = xad·if0		характеризует стационарную ЭДС холостого хода;
		–	постоянная времени обмоток и возбуждения при разомкнутой обмотке статора.

от (4) и (5) получается
(6)
При учете переходных процессов в демпферных обмотках, также расположенных на роторе, вводят в рассмотрение сверхпереходную ЭДС .
Сверхпереходные реактивные сопротивления по продольной и поперечной осям (x_d” и x_q”), тогда сверхпереходные процессы в СГ по продольной оси по аналогии с (2) записывается в виде
(7)
Поскольку эффект демпферных обмоток сказывается и по поперечной оси, то сверхпереходные процессы в СГ имеют место и поперечной оси
(8)
Окончательно, система упрощенная уравнениями электромагнитных переходных процессов СГ:
T_d0 dE_q’ /dt= E_qe – E_q’ – (x_d –x_d’)·i_d
T_q0 dE_d”/dt = –E_d” + (x_q – x_q”)·i_q (9)
T_d0”dE_q”/dt = E_q’ – E_q” – (x_d’ – x_d”)·i_d
Уравнения движения ротора [5]:
(10)

Уравнения АРВ и АРС учтены двухзвенной моделью [5].
Система дифференциальных уравнений описывающих переходные процессы в возбудителе и АРВ, получено согласно [5].
T_e dE_qe /dt = U_p + E_qe(0) – E_qe; E_qmin ≤ E_qe ≤ E_qmax
T_p dU_p /dt = U_Вх – U_p; U_pmin ≤ U_p ≤ U_pmax (11)

где
	Eqe	–	напряжение возбудителя, это синхронная ЭДС статора в установившемся режиме при напряжении возбудителя Uf;
	Eqe(0)	–	уставка возбудителя;
	Up	–	напряжение возбуждения возбудителя (напряжение на выходе регулятора);
	UВх	–	напряжение на входе регулятора;
	U(0)	–	уставка регулятора напряжения;
	Te, Tp	–	постоянные времени возбудителя и регулятора;
	Ku,K’u	–	коэффициенты усиления по каналам отклонения, напряжения и первой производной напряжения.
	Ks,K’s	–	коэффициенты усиления по каналам отклонения, скольжения и первой производной скольжения.

Уравнение АРС:

0
где
(12)


где ;
В различных областях науки и техники, широко применяются для решения дифференциальных уравнений метод Рунге-Кутта [1, 2]. Эти методы относятся к одношаговым методам, что в них интегрирование не на каждом шаге начинается заново, а скорости изменения интегрируемы переменных, целиком используются только внутри шага.
Среди многочисленных методов Рунге-Кутта наиболее широко используется вариант одношагового метода, где по четырем изменениям скорости каждой интегрируемой переменной вычисляется средневзвешенная скорость z_ср для определения приращения переменной Δy_n+1 в пределах одного шага интегрирования h т.е.
y_n+1 = y_n + hz_ср (13)
Наиболее распространена следующая формула Рунге-Кутта
(14)
где K₁ = hz(t,y_n); ; ;
.
С помощью этих формул значения y_n+1 получают исходя из значения y_n после четырех подстановок переменных в каждое дифференциальное уравнение системы. Для сложных уравнений это может потребовать значительного машинного времени на каждом шаге. При использовании в расчетах переходных процессов ЭЭС явного метода Рунге-Кутта приведенные выше (9 12) уравнения должны быть сведены к форме Коши. Отметим, что данная процедура связана с определенными трудностями и для сложных ЭЭС не всегда выполнима.
Особенности уравнений переходных процессов в ЭЭС позволяют предложить другой путь интегрирования, основанный на предварительном получении разностных уравнений отдельных элементов ЭЭС с применением неявных методов и последующим объединением этих уравнений в единую систему алгебраических нелинейных уравнений [3].
В соответствии с предложенным в [3] подходом запишем систему дифференциальных уравнений
(15)
Для m+1 -го расчетного шага применением формул неявных методов можно представить разностное соотношение:
, (16)

где	α	–	коэффициенты разностных уравнений, зависящие от коэффициентов исходных уравнений и шага интегрирования;
		–	вектор характеристика предысторий, элементы которого определяются теми же параметрами, что и и значений переменных на предшествующим шаге.

Для указанных выше неявных методов выражения коэффициентов (16) имеют следующий вид:
- метод Эйлера
, ,
- формула трапеций


- метод Эйлера с подгоночными коэффициентами
;
; ,

h -

–

шаг интегрирования. малый параметр при производной.

В соответствии с изложенным способом уравнения синхронного генератора (9 для m+1 -го расчетного шага преобразуем к виду (индекс шага m+1 здесь опущен):

(17)


(18)


(19)

μ_t= α_ψ ψ + μ_{t δ} (20)
μ _ρ = α _τ μ_τ + μ _{ρ δ}


(21)

Уравнения состояния электрической сети формируется на основе I- закона Кирхгофа в виде узловых уравнений (УУ) в форме баланса токов [6], и в синхронно вращающихся осях имеют вид:

, (22)
где Y - квадратная матрица собственных и взаимных узловых
проводимостей порядка N*N ( N – число узлов в схеме) ;
U – матрица – столбец узловых напряжений:
I - матрица – столбцы узловых токов.
Применение разностного подхода для решения жестких дифференциальных уравнений сводит задачу расчета ПП во времени к формированию обобщенных уравнений исследуемой системы в виде системы линейных и нелинейных алгебраических уравнений и совместному их решению на каждом временном шаге численного интегрирования.
Для совместного решения на временном шаге токи от генераторных узлов могут быть представлены в виде:
. (23)
Здесь коэффициенты и вектор предыстории получаются из системы (17 последовательным исключением внутренних переменных синхронного генератора (СГ) и имеют вид:


(24)

Таким образом, для коэффициентов уравнений (24) можно записать

; ; ; .
Введем преобразование перехода между координатами i -го СГ и синхронно вращающимися координатами П:


Cos (δ _i ) Sin (δ _i )
П _i = Sin (δ _i ) - Cos (δ _i ) (25)

где δ _i - угол между осью ротора и синхронно вращающимися осями, узловые токи от генераторов в координатах сети можно представить в виде

I . (26)
где - квазидиагональная квадратная матрица, элементы которой получаются из (23) с учетом П:

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

Для исследования численных методов были разработаны программы расчета ПП в ЭЭС при различных возмущениях (наброс мощности на валу турбины, К.З., отключения и включения ЛЭП).
В предложенных моделях расчетов электромеханических переходных процессов также реализован алгоритм выбора величины шага численного интегрирования. Его значение корректируется на основе анализа скорости изменения медленных переменных – скольжения СГ. В качестве критерия изменения шага принято отклонение в значениях скорости изменения углов СГ, полученных при итерационном решении (22), (26) и численным дифференцированием по предшествующим значениям.
Как видно из результатов, приведенных на рис. 2. а, б, метод Эйлера, являющийся простейшим неявным методом, имеет большую погрешность и не обеспечивает необходимую точность даже при уменьшении шага интегрирования. Повышение точности этого метода с помощью подгоночных коэффициентов позволяет производить расчеты с значительно большим шагом. Однако как показывает анализ, величина шага интегрирования в методе Эйлера с подгонкой все же ограничена точностью вычисления экспоненты в . Вместе с тем, при реализации предложенной стратегии выбора шага