无调节器的隐式发电机¶
模型描述¶
本示例考虑了汽轮发电机-变压器机组通过双回路架空线路 (OHL) 到无限功率母线 (IPB) 的运行情况,不带励磁和速度控制器。在模型中 5 秒的时间瞬间,其中一条架空线路的起始点发生了三相短路 (SC)。采用面积法对短路跳闸时限进行了分析计算,并与数学建模结果[1-3]进行了比较。
演示了使用命令控制从脚本开发环境启动和设置模型的过程、仿真结果的处理、仿真结果的可视化以及建议的模型独立工作脚本。对信号进行记录,并显示其时间图。模型外部视图:
本示例中使用的主要模块及其参数: 1.Synchronous Machine Round Rotor - 隐极同步发电机 ($P_{ном} = 500 МВт$,$U_{ном} = 24 кВ$). 2.电压源(三相)* - 用于模拟 SWBM 的三相电压源,通过指定有效线电压和相移来设置稳态 ($E = 500\angle0° кВ$,$S_{кз} = 10 ГВт$). 3.3. 故障(三相) - 短路。 4. 来自 Library of Passive Elements 的三相输电线路 Three-Phase PI Section Line ($L = 200км$, AC 500/64, 2 circuits) 和双绕组三相变压器 Two-Winding Transformer (Three-Phase) (TD-630000/500) 块。
用面积法计算跳闸角限值¶
导入处理图形所需的模块
using Plots;
gr();
GMS 参数:
# номинальное напряжение, В
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 + Xdd * 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 * (Xdd * 1im * Us ^ 2 / Sg)
Z12d = (Ua + I1 * (Xd * 1im * Us ^ 2 / Sg)) / I2
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);
使用暂态 EMF$E`_q$ 计算角度特性、 在短路时间的初始时刻和跳闸后保持不变:
# вырабатываемая генератором мощность в нормальном режиме
P0 = 400e6;
Q0 = 0;
# синхронная ЭДС по оси q
Eq = sqrt((Us + Q0 * abs(Z12d) / Us) ^ 2 + (P0 * abs(Z12d) / Us) ^ 2)
# угол ротора в исходном режиме
d0 = atan(P0 * abs(Z12d) / (Us ^ 2 + Q0 * abs(Z12d)))
# переходная ЭДС
Ep = sqrt((Us + Q0 * abs(Z12) / Us) ^ 2 + (P0 * abs(Z12) / Us) ^ 2)
# угол переходной ЭДС
d0p = atan(P0 * abs(Z12) / (Us ^ 2 + Q0 * abs(Z12)))
# переходная ЭДС по оси q в нормальном режиме
Eqq = Ep * cos(d0p - d0)
# переходная ЭДС по оси q в аварийном и послеаварийном режимах
Eqq_pav = Eqq
Eqq_av = Eqq
# максимальная активная мощность в нормальном режиме
Pmax = (Eqq ^ 2 / abs(Z11) * sin(alpha11) + Eqq * Us / abs(Z12)) / 1e6
# угол переходной ЭДС по оси q
d0p = asin(P0 / Pmax / 1e6)
# максимальная мощность в аварийном режиме
Pmax_av = 0
# максимальная мощность в послеаварийном режиме
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 * (pi - d0_pav - d0p) + Pmax_pav * cos(pi - d0_pav) - Pmax_av * cos(d0p)) / (Pmax_pav - Pmax_av));
区域法的可视化:
# функции угловых характеристик
delta = 0:0.1:180
Pmax = x -> (Eqq ^ 2 / abs(Z11) * sin(alpha11) .+ (Eqq * Us / abs(Z12)) * sin.(x .- alpha12)) / 1e6
Pmax_av = x -> fill(0, size(x))
Pmax_pav = x -> (Eqq_pav ^ 2 / abs(Z11_pav) * sin(alpha11_pav) .+ Eqq_pav * Us / abs(Z12_pav) * sin.(x .- alpha12_pav)) / 1e6
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 = [L"Норм.реж." L"Авар.реж." L"П/АВ.реж." L"Мощн.турбины"], ylabel = L"P, МВт", xlabel = L"\delta\prime, град", legend=false,
xaxis = (0:30:180), size = (900,440), left_margin=5Plots.mm, bottom_margin=5Plots.mm)
plot!((d0p:0.01:d_pr) .* 180 / pi, Pmax_av((d0p:0.01:d_pr)), fillrange = Pt((d0p: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!([d0p,d_pr,dkr] .* 180 / pi, c = 5)
annotate!([(90, Pmax(pi / 2) - 50, (L"Норм.реж.", 10, :black)),
(50, Pmax_av(pi / 2) .+ 30, (L"Авар.реж.", 10, :black)),
(90, Pmax_pav(pi / 2) - 50, (L"ПА/В.реж.", 10, :black)),
(d0p * 180 / pi + 20, P0 * 1e-6 - 200, (L"F_{уск}", 15, :black, :left)),
(d_pr * 180 / pi + 15, P0 * 1e-6 + 200, (L"F_{торм}", 15, :black, :left)),
(0, P0 * 1e-6 + 50, (L"P_{турб}", 10, :black, :left)),
(d0p * 180 / pi, 30, (L"\delta\prime_0 = %$(round(d0p * 180 / pi, digits=2))", 12, :black, :right)),
(d_pr * 180 / pi, 30, (L"\delta\prime_{пред} = %$(round(d_pr * 180 / pi, digits=2))", 12, :black, :left)),
(dkr * 180 / pi, 30, (L"\delta\prime_{крит} = %$(round(dkr * 180 / pi, digits=2))", 12, :black, :right))])
短路跳闸的极限时间:
Tj = 7.05
tpr = sqrt((d_pr - d0p) * 2 * Tj / (0.8 * 2 * pi * 50))
print("t_пред = $(round(tpr, digits = 3))")
建模¶
加载模型
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.24 c,以保持动态稳定性:
step = "Отключение КЗ";
# отключение цепи ВЛ
engee.set_param!(model_name * "/" * step, "Time" => 5.24);
运行加载的模型:
results = engee.run(model_name);
为导入仿真结果,已提前启用了所需信号的日志记录并设置了其名称。让我们从 results 变量中转换模拟结果:
# вектор времени симуляции
sim_time = results["w"].time;
# вектор скорости вращения ротора
w1 = results["w"].value;
# вектор тока генератора
P1 = results["P"].value;
短路持续时间增加到 0.2402 秒时的模拟结果:
# отключение цепи ВЛ
engee.set_param!(model_name*"/"*step, "Time" => 5.2402);
results = engee.run(model_name);
# вектор скорости вращения ротора
w2 = results["w"].value;
# вектор тока генератора
P2 = results["P"].value;
发电机有功功率图:
plot(sim_time, P1, title = L"Мощность\; генератора", ylabel = L"P, о.е.", xlabel="Время, c", label = L"t_{КЗ} = 0.24 с", legendfontsize = 10)
plot!(sim_time, P2, label = L"t_{КЗ} = 0.2402 с", ylims = (-1, 2),size = (900,440), left_margin=5Plots.mm, bottom_margin=5Plots.mm)
转子转速图
plot(sim_time, w1, title = L"Скорость\; ротора", ylabel = L"\omega, о.е.", xlabel=L"Время, c", label = L"t_{КЗ} = 0.24 с", legendfontsize = 10)
plot!(sim_time, w2, label = L"t_{КЗ} = 0.2402 с", ylims = (0.95, 1.1),size = (900,440), left_margin=5Plots.mm, bottom_margin=5Plots.mm)
在第二次实验中,发电机的动态稳定性受到破坏,因此模型中的短路跳闸极限时间是根据经验得出的。结果与理论计算趋同,误差很小。
println("Предельное время отключения аналитически = "*string(round(tpr, digits=3)) * "с")
println("Предельное время отключения практически = "*string(0.24) * "с")
附录¶
尝试自行更改以下模型参数,并研究其对动态稳定性的影响:
- 将架空线路长度增加 200 km;
- 故障(三相) 块中的短路类型;
- 将初始模式下发电机的有功功率降低 200 兆瓦。
结论¶
在本示例中使用了 Engee 模型的指令控制和仿真结果上传工具,并展示了 Plots 模块的工作。考虑了隐极发电机通过变压器和双回路架空线路在 SWBM 上的运行,不带励磁和速度控制器。对动态稳定性进行了分析计算,并将计算结果与模拟数据进行了比较。
参考文献¶
1.Venikov V. A. Transient electromechanical processes in electrical systems: Textbook for electric power engineering specialities of universities.- 4 ed., revision and supplement. 2.举例说明电气系统的瞬态:大学教材》(V. V. Ezhkov, N. I. Zelenohat, I. V. Litkens 等;V. A. Stroev 编辑)。- 莫斯科:Znak,1996 年。- 224 pp. 3.Dolgov, A.P. Stability of electrical systems : a textbook / Dolgov A.P..- Novosibirsk : Novosibirsk State Technical University, 2010.- 177 c.- ISBN 978-5-7782-1320-3.