Неявнополюсный генератор без регуляторов

Описание модели

     В этом примере рассматривается работа блока турбогенератор-трансформатор через двухцепную воздушную линию (ВЛ) на шину бесконечной мощности (ШБМ) без регуляторов возбуждения и скорости вращения. В момент времени 5 с. в модели происходит трёхфазное короткое замыкание (КЗ) в середине одной из цепей ВЛ. Проводится аналитический расчёт динамической устойчивости системы с помощью метода площадей и сравнение с результатами математического моделирования [1,2].

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

    Основные блоки используемые в данном примере:

1. Synchronous Machine Round Rotor - неявнополюсный синхронный генератор.
2. Voltage Source (Three-Phase) - трёхфазный источник напряжения, используется для моделирования ШБМ, установившийся режим выставлен с помощью задания действующего линейного напряжения и фазового сдвига.
3. Fault (Three-Phase) - короткое замыкание.
4. Из библиотеки пассивных элементов блоки трёхфазной ЛЭП Three-Phase PI Section Line и двухобмоточного трехфазного трансформатора Two-Winding Transformer (Three-Phase).

    Параметры системы:

Элемент	Параметр
ШБМ	Эквивалентная ЭДС $E = 500\angle0° кВ$ $S_{кз} = 10 ГВт$
Генератор	$P_{ном} = 500 МВт$ $U_{ном} = 24 кВ$ $P_{турбины} = 400 МВт$
ВЛ	$U_{ном} = 500 кВ$ $L = 200км$ АС 500/64 2 цепи
ТР	ТД-630000/500/24

Расчёт предельного угла отключения методом площадей

    Импортируем необходимые модули для работы с графиками:

using Plots;
plotlyjs();

    Параметры ШБМ:

# номинальное напряжение, В
Us = 500e3
# мощность КЗ, ВА
Skz = 10e9
# полное сопротивление системы, Ом
Zkz = Us ^ 2 / Skz
# соотношение активного и индуктивного сопротивления
XR = 15
# индуктивное сопротивление системы, Ом
Xkz = Zkz * sin(atan(XR))
# активное сопротивление системы, Ом
Rkz = Zkz * cos(atan(XR))
# комплексное сопротивление системы, Ом
Zs = Rkz + 1im * Xkz;

    Параметры ЛЭП:

# удельное сопротивление прямой последовательности Ом/км
Z0 = 0.01967 + 0.2899im
# количество цепей
n = 2
# длина ВЛ, км
L = 200
# удельная емкостная проводимость, См/км
b0 = 3.908e-6
bc = b0 * L * n
# полное сопротивление ВЛ, Ом
Zl = Z0 * L / n
# емкостное сопротивление половины ВЛ, Ом
Zc = -2im / bc;

    Параметры трансформатора:

# активное сопротивление ТР, Ом
Rt = 0.4514
# индуктивное сопротивление ТР, Ом
Xt = 30.62
Zt = Rt + Xt * 1im;

    Параметры генератора:

# номинальное напряжение, В
Ug = 24e3
# номинальная мощность, ВА
Sg = 500e6
# активное сопротивление, о.е.
Ra = 0.003
# синхронное индуктивное сопротивление по оси d, о.е.
Xd = 1.81
# переходное индуктивное сопротивление по оси d, о.е.
Xdd= 0.3
# комплексное сопротивление генератора, Ом
Zg = (Ra + Xd * 1im) * Us ^ 2 / Sg;

    Параметры нагрузки:

# активная мощность нагрузки, Вт
Pn = 20e6
# активное сопротивление нагрузки, Ом
Zn = Us ^ 2 / Pn;

    Расчёт собственных и взаимных сопротивлений $Z_{11}$ и $Z_{12}$ в нормальном режиме методом единичного тока:

I2 = 1
Uc = I2 * Zs
Ic0 = Uc / Zc
Ibc = I2 + Ic0
Ub = Uc + Zl * Ibc
Ib0 = Ub / Zc
Iab = Ib0 + Ibc
Ua = Ub + Iab * Zt
Ia0 = Ua / Zn
I1 = Ia0 + Iab
U1 = Ua + I1 * Zg
Z12 = U1 / I2
Z11 = U1 / I1
alpha11 = pi / 2 - angle(Z11)
alpha12 = pi / 2 - angle(Z12);

    Расчёт собственных и взаимных сопротивлений $Z_{11}$ и $Z_{12}$ в послеаварийном режиме методом единичного тока:

I2 = 1
Uc = I2 * Zs
Ic0 = Uc / (2 * Zc)
Ibc = I2 + Ic0
Ub = Uc + 2 * Zl * Ibc
Ib0 = Ub / (2 * Zc)
Iab = Ib0 + Ibc
Ua = Ub + Iab * Zt
Ia0 = Ua / Zn
I1 = Ia0 + Iab
U1 = Ua + I1 * (Xdd * 1im * Us ^ 2 / Sg)
Z12_pav = U1 / I2
Z11_pav = U1 / I1
alpha11_pav = pi / 2 - angle(Z11_pav)
alpha12_pav = pi / 2 - angle(Z12_pav);

    Расчёт собственных и взаимных сопротивлений $Z_{11}$ и $Z_{12}$ в аварийном режиме методом единичного тока:

I2 = 1
Uc = I2 * Zs
Ic0 = Uc / (2 * Zc * Zl / (2 * Zc + Zl))
Ibc = I2 + Ic0
Ub = Uc + Zl * 2 * Ibc
Ib0 = Ub / (2 * Zc * Zl / (2 * Zc + Zl))
Iab = Ib0 + Ibc
Ua = Ub + Iab * Zt
Ia0 = Ua / Zn
I1 = Ia0 + Iab
U1 = Ua + I1 * (Xdd * 1im * Us ^ 2 / Sg)
Z12_av = U1 / I2
Z11_av = U1 / I1
alpha11_av = pi / 2 - angle(Z11_av)
alpha12_av = pi / 2 - angle(Z12_av);

    Расчёт угловых характеристик производится со следующими допущениями:

1. Переходная ЭДС $E`_q$ остаётся неизменной в начальный момент времени и после отключения КЗ.
2. Угловые характеристики для аварийного и послеаварийного режмов построены при $E`_q=E`_{q0}=const$.
3. Учитываются активные сопротивления всех элементов, зарядные емкости ВЛ.

# исходный режим
# вырабатываемая генератором мощность в нормальном режиме
P0 = 400e6;
Q0 = 0;
# напряжение на шинах генератора в исходном режиме
Ug = Us * 25.5 / 24
# синхронная ЭДС по оси q
Eq = sqrt((Us + Q0 * Xd * Us ^ 2 / Sg / Us) ^ 2 + (P0 * Xd * Us ^ 2 / Sg / Us) ^ 2)
# переходная ЭДС по оси q в аварийном и послеаварийном режимах
Eqq_pav = sqrt((Us + Q0 * Xdd * Us ^ 2 / Sg / Us) ^ 2 + (P0 * Xdd * Us ^ 2 / Sg / Us) ^ 2)
# максимальная активная мощность в нормальном режиме
Pmax = (Eq ^ 2 / abs(Z11) * sin(alpha11) + Eq * Us / abs(Z12)) / 1e6
# угол ротора в исходном режиме
d0 = asin(P0 / Pmax / 1e6)
# максимальная мощность в аварийном режиме
Pmax_av = (Eqq_pav ^ 2 / abs(Z11_av) * sin(alpha11_av) + Eqq_pav * Us / abs(Z12_av)) / 1e6
# максимальная мощность в послеаварийном режиме
Pmax_pav = (Eqq_pav ^ 2 / abs(Z11_pav) * sin(alpha11_pav) + Eqq_pav * Us / abs(Z12_pav)) / 1e6
# угол ротора по работе на послеаварийной характеристике
d0_pav = asin(P0 / Pmax_pav / 1e6)
# критический угол
dkr = pi - d0_pav
# предельный угол отключения
d_pr = acos((P0 * 1e-6 * (dkr - d0) + Pmax_pav * cos(dkr) - Pmax_av * cos(d0)) / (Pmax_pav - Pmax_av));

    Визуализация метода площадей:

# функции угловых характеристик
delta = 0:0.1:180
Pmax = x -> (Eq ^ 2 / abs(Z11) * sin(alpha11) + Eq * Us / abs(Z12)) / 1e6 * sin.(x .- alpha12)
Pmax_av = x -> (Eqq_pav ^ 2 / abs(Z11_av) * sin(alpha11_av) + Eqq_pav * Us / abs(Z12_av)) / 1e6 * sin.(x .- alpha12_av)
Pmax_pav = x -> (Eqq_pav ^ 2 / abs(Z11_pav) * sin(alpha11_pav) + Eqq_pav * Us / abs(Z12_pav)) / 1e6 * sin.(x .- alpha12_pav)
Pt = x -> fill(P0 * 1e-6, size(x))
# визуализация метода площадей
plot(delta, [Pmax(delta .* pi ./ 180) Pmax_av(delta .* pi ./ 180) Pmax_pav(delta .* pi ./ 180) Pt(delta)],
label = ["Норм.реж." "Авар.реж." "П/АВ.реж." "Мощн.турбины"], ylabel = "P, МВт", xlabel = "delta, град", legend=false,
xaxis = (0:30:180))
plot!((d0:0.01:d_pr) .* 180 / pi, Pt((d0:0.01:d_pr)), fillrange = Pmax_av((d0:0.01:d_pr)), fillalpha = 0.2, c = 2, label = "Fуск")
plot!((d_pr:0.01:dkr) .* 180 / pi, Pmax_pav((d_pr:0.01:dkr)), fillrange = Pt((d_pr:0.01:dkr)), fillalpha = 0.2, c = 1, label = "Fторм")
vline!([d0,d_pr,dkr] .* 180 / pi, c = 5)
annotate!([(90, Pmax(pi / 2) + 50, ("Норм.реж.", 10, :black, :left)),
        (90, Pmax_av(pi / 2) - 50, ("Авар.реж.", 10, :black, :left)),
        (90, Pmax_pav(pi / 2) - 50, ("ПА/В.реж.", 10, :black, :left)),
        (d0 * 180 / pi + 15, P0 * 1e-6 - 100, ("Fуск.", 10, :black, :left)),
        (d_pr * 180 / pi + 15, P0 * 1e-6 + 100, ("Fторм.", 10, :black, :left)),
        (0, P0 * 1e-6 + 50, ("Pтурб.", 10, :black, :left)),
        (d0 * 180 / pi, 0, ("d0=" * string(round(d0 * 180 / pi, digits=2)), 10, :black, :right)),
        (d_pr * 180 / pi, 0, ("dпред=" * string(round(d_pr * 180 / pi, digits=2)), 10, :black, :right)),
        (dkr * 180 / pi, 0, ("dкрит=" * string(round(dkr * 180 / pi, digits=2)), 10, :black, :right))])

    На графике площадей отображены значения предельного и критического углов ротора:

println("Предельный угол по методу площадей = "*string(round(d_pr * 180 / pi,digits=2)) * "°")
println("Критический угол по методу площадей = "*string(round(dkr * 180 / pi,digits=2))*"°")

Предельный угол по методу площадей = 116.86°
Критический угол по методу площадей = 157.25°

Моделирование

    Загрузка модели:

model_name = "synchronous_generator_round_rotor";
model_name in [m.name for m in engee.get_all_models()] ? engee.open(model_name) : engee.load( "$(@__DIR__)/$(model_name).engee");

    Настройка модели на время КЗ 0,3062 c при котором ротор почти достигнет критического угла с помощью командного управления:

step = "Отключение КЗ";
# отключение цепи ВЛ
engee.set_param!(model_name * "/" * step, "Time" => 5.3062);

    Запуск загруженной модели:

results = engee.run(model_name);

    Для импорта результатов моделирования было заранее включено логгирование необходимых сигналов и заданы их имена. Преобразуем результаты моделирования из переменной results:

# вектор времени симуляции
sim_time = results["Rotor angle"].time;
# вектор угла ротора
rotor_angle = results["Rotor angle"].value;
d_pr_model = rotor_angle[sim_time .== 5.3062][1];
# вектор тока генератора
I = results["I"].value;

    График тока генератора:

plot(sim_time, I, title = "Ток генератора", ylabel = "I, о.е.", xlabel="Время, c", legend = false)

    График угла ротора относительно угла напряжения ШБМ:

plot(sim_time, rotor_angle, title = "Угол ротора", ylabel = "delta, град", xlabel="Время, c", legend = false)
plot!([4.5, 5.3062], [130, rotor_angle[sim_time .== 5.3062][1]], c=:black)
plot!([5.3062], [rotor_angle[sim_time .== 5.3062][1]], seriestype=:scatter, c=:orange, markerstrokewidth = 0)
annotate!(4.5, 130, ("Отключение 1-й цепи ВЛ", 10, :black, :right))
plot!([4.5, 5], [30, rotor_angle[sim_time .== 5][1]], c=:black)
plot!([5], [rotor_angle[sim_time .== 5][1]], seriestype=:scatter, c=:orange, markerstrokewidth = 0)
annotate!(4.5, 30, ("Возникновение КЗ", 10, :black, :right))

     При длительности КЗ 0,3062 с динамическая устойчивость генератора сохраняется. Отключение КЗ происходит при угле ротора 115 градусов. В следующем опыте увеличим длительность КЗ до 0,3064 секунды.

     Увеличение временени КЗ до 0,3064 с и запуск модели:

# отключение цепи ВЛ
engee.set_param!(model_name*"/"*step, "Time" => 5.3064);
results = engee.run(model_name);
# вектор времени симуляции
sim_time = results["Rotor angle"].time;
# вектор угла ротора
rotor_angle = results["Rotor angle"].value;
# вектор тока генератора
I = results["I"].value;

    График тока генератора:

plot(sim_time, I, title = "Ток генератора", ylabel = "I, о.е.", xlabel="Время, c", legend = false)

    График угла ротора относительно угла напряжения ШБМ:

plot(sim_time, rotor_angle, title = "Угол ротора", ylabel = "delta, град", xlabel="Время, c", legend = false)
plot!([4.5, 5.3064], [130, rotor_angle[sim_time .== 5.3064][1]], c=:black)
plot!([5.3064], [rotor_angle[sim_time .== 5.3064][1]], seriestype=:scatter, c=:orange, markerstrokewidth = 0)
annotate!(4.5, 130, ("Отключение 1-й цепи ВЛ", 10, :black, :right))
plot!([4.5, 5], [30, rotor_angle[sim_time .== 5][1]], c=:black)
plot!([5], [rotor_angle[sim_time .== 5][1]], seriestype=:scatter, c=:orange, markerstrokewidth = 0)
annotate!(4.5, 30, ("Возникновение КЗ", 10, :black, :right))

     Динамическая устойчивость генератора нарушилась, таким образом были найдены предельный угол отключения КЗ и предельное время отключения. Результат с небольшой погрешностью сходится с теоретическим расчётом по методу площадей.

println("Предельный угол по методу площадей = "*string(round(d_pr * 180 / pi, digits=1))*"°")
println("Предельный угол из модели = "*string(round(d_pr_model, digits=1))*"°")

Предельный угол по методу площадей = 116.9°
Предельный угол из модели = 115.0°

Дополнение

     Попробуйте самостоятельно изменить следующие параметры модели и исследовать, как это влияет на динамическую устойчивость:

1. увеличить длину ВЛ на 200 км;
2. вид КЗ в блоке Fault (Three-Phase);
3. снизить активную мощность генератора в исходном режиме на 200 МВт.

Выводы

    В данном примере были использованы инструменты для командного управления моделью Engee и выгрузки результатов моделирования, показана работа с модулем Plots. Была рассмотрена работа неявнополюсного генератора через трансформатор и двухцепную ВЛ на ШБМ без регуляторов возбуждения и скорости вращения. Проведен аналитический расчёт динамической устойчивости с помощью метода площадей и выполнено сравнение результатов с данными моделирования.

Ссылки

1. Веников В. А. Переходные электромеханические процессы в электрических системах: Учеб. для электроэнергет. спец. вузов. — 4-е изд., перераб. и доп. — М.: Высш. шк., 1985
2. Переходные процессы электрических систем в примерах и иллюстрациях: Учебное пособие для вузов (В. В. Ежков, Н. И. Зеленохат, И. В. Литкенс и др.; Под ред. В. А. Строева). - М.: Знак, 1996. - 224с., ил.