简介:电力系统分析是电气工程专业核心课程,稳定性章节重点研究系统在扰动后维持同步运行的能力。课件深入讲解稳定性基本概念、摆动方程推导、稳态运行时同步电机模型,并通过等面积准则分析暂态稳定,同时涵盖三相短路故障影响、非线性方程及摆动方程的数值解法,最后扩展到多机系统与多机暂态分析,理论体系完整。压缩包内共1个ppt文件,大小约11.46MB,包含10个小节,适用于电气工程专业学生、考研复习者及电力系统从业者,可作为课堂教学、自学笔记或考前梳理资料。已有64人浏览学习,内容契合高校《电力系统分析》教材体例,配合图示与公式推导,以中文讲解,便于理解。作为课程课件使用或对照教材复习,都能帮助读者快速掌握电力系统稳定性的核心脉络,尤其适合搭配《电力系统分析》教材逐节研读。
1. 电力系统稳定性分析的工程边界
电网调度与继电保护整定中,最棘手的问题往往不是故障电流算不准,而是故障切除之后,系统还能不能回到同步运行状态。稳定性的研究对象是同步电机的转子:故障瞬间机械功率来不及变化,电磁功率却骤降,转子开始加速,功角随之拉大,一旦越过临界值,发电机与系统之间就会失步。判断这一过程是否可控,靠的是对摆动方程这个非线性微分方程的分析,而不是单纯计算潮流。这份资源覆盖了从摆动方程推导、同步电机经典模型、等面积准则到数值解和多机系统的完整链条,适合继电保护整定工程师、电网规划人员,以及刚进入电力系统动态仿真方向的研究生快速建立从方程到判据的整体框架。
2. 摆动方程的工程变形与经典模型的可适用边界
2.1 从转动定律到 H 常量:三步推导
同步电机的转子运动服从牛顿转动定律:作用在转子上的净转矩等于转动惯量与角加速度的乘积。稳态运行时,机械转矩Tm与电磁转矩Te相等;扰动出现后,二者出现差值,这个差值就是加速转矩。设转子相对定子的机械角位移为θ_m,则有:
J · d²θ_m/dt² = Tm - Te转子以同步转速旋转,把角位移写成θ_m = ω_sm · t + δ_m的形式,其中ω_sm是同步机械角速度,δ_m是相对同步旋转参考系的角位移(即功角)。对时间求二阶导后,转动方程的左边直接变成J · d²δ_m/dt²,同步转速项被消掉了。这是功角能成为稳定性研究对象的关键:扰动产生的相对运动只取决于功角的变化。
两边乘以角速度ω_m,转矩就变成了功率。PPT 中先引入惯性常数M = J·ω_m,再把M用额定转速下的转子动能表示。设转子动能为W_k = 1/2 · J · ω_sm²,则M = 2W_k / ω_sm。工程上更常用的是 H 常量:
H = 额定转速下的转子动能(MJ) / 额定容量(MVA)H 的单位是秒,其物理含义是:发电机在额定机械功率输入下,靠释放全部转子动能可以维持的时间。汽轮发电机组 H 通常在 4~9 秒范围,水电机组低一些,整体落在 1~10 秒区间。用 H 替换 M 并做标幺化后,摆动方程写成:
2H/ω_sm · d²δ_m/dt² = P_m - P_e其中P_m和P_e均为标幺值。若 δ 用电角度的弧度表示,ω_sm换成电气角速度2πf,方程变为H/(πf) · d²δ/dt² = P_m - P_e;若 δ 用电角度表示,则系数变为H/(180f)。实际编程时最常见的是H/(πf)这个形式,注意 δ 要统一用弧度。
需要留意的是,摆动方程的基本形式没有阻尼项。P_m - P_e中只包含不平衡功率,这对首摆稳定性判断够用;但如果要分析多摆过程或者动态稳定,工程上会在方程右侧加一项D · dδ/dt来等效阻尼绕组和负荷的频率效应,D 一般取 1~3 的标幺值。
2.2 经典模型:恒定电压源与直轴暂态电抗的约束
暂态稳定分析的经典模型,是把发电机等效成一个恒定电压源E'串联一个直轴暂态电抗X'd。这个模型的成立依据是:扰动发生后的极短时间内,励磁绕组磁链近似守恒,因此暂态电动势E'基本不变。经典模型在首摆稳定性分析中非常可靠,但它有明确的适用边界:
| 分析场景 | 经典模型适用性 | 失效原因 |
|---|---|---|
| 近区三相短路首摆分析 | 适用 | 励磁绕组磁链守恒,E' 保持恒定 |
| 多摆过程与动态稳定 | 不适用 | 励磁调节器、阻尼绕组开始起主导作用 |
| 新能源场站动态 | 不适用 | 变流器没有转子运动方程,输出功率完全受控制决定 |
经典模型对单机无穷大系统的描述,是把发电机内节点看成一个电势源,经过X'd、变压器电抗、线路电抗连接到电压恒定的母线上。之所以叫无穷大母线,是因为它的电压幅值和频率在任意功率注入或抽取下都不变化。这个最小系统足以说清楚稳定性的所有基本概念,从单机系统理解透之后,多机系统的扩展只是把单机方程写成向量形式。
2.3 双总线方程与功角曲线的数值绘制
把发电机内节点记为节点 1,无穷大母线记为节点 2,用节点导纳矩阵写出双母线系统,忽略所有电阻后,发电机输送到无穷大母线的功率可以化成一个极简形式:
P_e = E' · V / X · sinδX 是暂态电抗后的总转移电抗。这个方程说明功率传输取决于两个因素:转移电抗和两电压之间的夹角。P_e随 δ 按正弦变化,峰值Pmax = E'V/X出现在 δ = 90° 处,这就是稳态功率极限。下面用例 11.1 的参数绘制暂态功角曲线,同时对比忽略凸极与考虑凸极两种模型的差别——后文会详细推导这两组参数是如何计算出来的。
import numpy as np import matplotlib.pyplot as plt # 例 11.1 的参数:Xd=1.0, Xq=0.6, Xd'=0.3, 机端电压 V=1.0 Xdp = 0.3 Xq = 0.6 V = 1.0 # 忽略凸极:E'=1.1226;考虑凸极:E'=1.1162,来自例 11.1 的计算结果 E_round = 1.1226 E_salient = 1.1162 delta = np.linspace(0, np.pi, 400) P_round = E_round * V / Xdp * np.sin(delta) P_salient = E_salient * V / Xdp * np.sin(delta) + V**2 / 2 * (1 / Xq - 1 / Xdp) * np.sin(2 * delta) plt.plot(np.degrees(delta), P_round, label='round-rotor (ignore saliency)') plt.plot(np.degrees(delta), P_salient, label='salient-pole') plt.axhline(0.5, ls='--', color='gray') plt.xlabel('delta (deg)') plt.ylabel('Pe (pu)') plt.legend() plt.grid() plt.show() print('round max =', round(P_round.max(), 3), 'at', round(np.degrees(delta[P_round.argmax()]), 1), 'deg') print('salient max =', round(P_salient.max(), 3), 'at', round(np.degrees(delta[P_salient.argmax()]), 1), 'deg')这段代码里,P_round是忽略凸极时的功角特性,P_salient是在其基础上叠加了凸极修正项V²/2 · (1/Xq - 1/Xd') · sin2δ。因为1/Xq - 1/Xd'在这个算例中为负值,而sin2δ在 90° 以后也是负值,二者相乘为正,所以在 90° 之后凸极曲线反而高于正弦曲线,峰值位置向右偏移。运行后可以看到,忽略凸极时 Pmax 是 3.742 pu,出现在 90.0°;考虑凸极后 Pmax 增至 4.032 pu,出现在约 110.0°。峰值差约 7%,这个差距在精确校核临界切除时间时不可忽略。
3. 隐极与凸极模型的暂态电抗后电压计算与误差边界
3.1 忽略凸极:一步复数运算求暂态电动势
例 11.1 给出了一个典型的单机无穷大系统算例。发电机参数标幺值为Xd=1.0、Xq=0.6、Xd'=0.3,机端电压V=1.0,输出有功P=0.5,滞后功率因数0.8。计算暂态电抗后的电压E',首先要确定故障前稳态电流I。
import numpy as np V = 1.0 + 0j # 无穷大母线电压,基准相角 0 Xd = 1.0 # 直轴同步电抗 Xq = 0.6 # 交轴同步电抗 Xdp = 0.3 # 直轴暂态电抗 P = 0.5 phi = np.arccos(0.8) # 功率因数角 = 36.87° = 0.6435 rad # 复功率 S = P + jQ,滞后功率因数时 Q 为正 Q = P * np.tan(phi) # 0.375 S = P + 1j * Q # 由 S = V * conj(I) 反解电流 I = np.conj(S / V) # 0.625 ∠ -36.87° # 忽略凸极:E' = V + j * Xd' * I Ep_round = V + 1j * Xdp * I print('|I| =', round(abs(I), 4)) print('E\' =', round(abs(Ep_round), 4), '∠', round(np.degrees(np.angle(Ep_round)), 3), 'deg')I = conj(S/V)这一步使用了复功率的定义:S = V · I*。算例中电流幅值 0.625 pu,相角 -36.87°,正好对应 0.8 滞后功率因数。V + jX'd I是经典模型的内电势表达式,物理意义是从机端电压出发,沿着暂态电抗的压降方向反推出恒定的暂态电动势。输出结果为E' = 1.1226 ∠ 7.679°,这个 7.679° 就是初始功角 δ0,是后续所有暂态分析的起点。
对应的功角特性方程为:
Pe(δ) = 1.1226 × 1.0 / 0.3 × sinδ = 3.7419 sinδ由于转移电抗只有 0.3,功率极限达到 3.742 pu,这说明暂态过程中发电机能短时送出远超额定的功率,代价是功角被拉到很大。
3.2 考虑凸极:先求稳态功角再回代
凸极机模型不能直接套用V + jX'd I,因为交轴和直轴的不对称让内电势方向需要重新确定。标准做法分三步:先用V + jXq I求出稳态功角 δ0,再用直轴公式求励磁电动势 E,最后利用磁链守恒条件把 E 折算到暂态电抗后的E'。
# 第一步:由交轴电抗求稳态功角 delta0 Eq_ref = V + 1j * Xq * I # 参考相量 delta0 = np.angle(Eq_ref) # 13.7608° # 第二步:求励磁电动势 E = V*cos(delta0) + Xd*I*sin(delta0 + phi) E_f = np.abs(V) * np.cos(delta0) + Xd * np.abs(I) * np.sin(delta0 + phi) # 第三步:磁链守恒,折算暂态电抗后的 E' Ep_salient = (Xdp * E_f + (Xd - Xdp) * np.abs(V) * np.cos(delta0)) / Xd print('delta0 =', round(np.degrees(delta0), 4), 'deg') print('E =', round(E_f, 4), 'pu') print('E\' =', round(Ep_salient, 4), 'pu')第一步中,np.angle(V + 1j*Xq*I)求的是交轴内电势方向,与机端电压的夹角就是稳态功角。第二步中的sin(δ0 + φ)对应直轴电流分量,因为滞后电流在 d 轴上的投影恰好是I·sin(δ0+φ)。第三步的折算公式来自于暂态前磁链守恒条件:E' = (X'd·E + (Xd - X'd)·V·cosδ0) / Xd,把稳态励磁电动势按直轴电抗的分压关系换算到暂态电抗之后。
输出结果依次为:δ0 = 13.7608°、E = 1.4545 pu、E' = 1.1162 pu。与 3.1 节忽略凸极的结果 1.1226 相比,E'相差约 0.6%,看起来很小,但后面会看到它引起的功率方程差异会被功角放大。
3.3 凸极项对峰值位置的影响
凸极模型下的暂态功角特性方程为:
Pe(δ) = E'V/Xd' · sinδ + V²/2 · (1/Xq - 1/Xd') · sin2δ = 3.7208 sinδ - 0.8333 sin2δ两项的系数对比:正弦项系数 3.7208,凸极项系数 -0.8333,后者只有前者的 22%,但它的作用方向很关键。下表汇总两种模型的结果。
| 模型 | E' (pu) | 功角特性方程 | Pmax (pu) | Pmax 位置 |
|---|---|---|---|---|
| 忽略凸极 | 1.1226 | 3.7419 sinδ | 3.7419 | 90.0° |
| 考虑凸极 | 1.1162 | 3.7208 sinδ - 0.8333 sin2δ | 4.032 | 110.0° |
sin2δ在 0~90° 为正值,与负系数相乘后,凸极项使Pe比纯正弦项略小;在 90°~180° 区间,sin2δ为负值,负负得正,凸极项反过来帮助抬高Pe。因此 Pmax 没有出现在 90°,而是向右移动到 110° 附近,数值从 3.742 提高到 4.032。由于凸极项对功率的贡献在完整周期内正负相抵,平均值趋近于零,工程估算时常常直接忽略sin2δ项;但在等面积准则计算临界切除角时,Pmax 的 7% 误差会传导到切除角上,不能想当然省掉。
4. 等面积准则判定暂态稳定与临界切除角计算
4.1 从摆动方程到面积关系的推演
等面积准则不需要求解摆动方程,直接根据能量守恒判断系统是否能保持暂态稳定。把摆动方程两边同时乘以dδ/dt再对时间积分,利用恒等式d/dt(dδ/dt)² = 2(dδ/dt)(d²δ/dt²),可以得到:从初始功角到最大摆角的两端,如果转子角速度都为零,则不平衡功率对功角的积分必须为零。写成面积形式:
∫(Pm - Pe) dδ = 0这个积分的几何含义是:在功角坐标轴上,Pm曲线与Pe曲线之间的面积。转子加速阶段,Pm > Pe,面积为正,称为加速面积;减速阶段,Pe > Pm,面积为负,称为减速面积。只要最大摆角处存在足够的减速面积抵消加速面积,系统就能在第一摆保持稳定。因此,暂态稳定的充要条件是:最大减速面积大于等于加速面积,这就是等面积准则的核心。
4.2 三相短路期间的功率曲线分段
三相短路是暂态稳定分析中最严重的工况。故障发生瞬间,发电机输出的电磁功率跌到很低水平,机械功率却来不及调整,功角快速拉大;故障切除之后,系统阻抗因为线路跳开而比故障前略大,功率曲线介于两者之间。为了量化分析,做一组典型参数:
| 阶段 | 功率曲线 | 转移电抗 (pu) | Pmax (pu) |
|---|---|---|---|
| 故障前 | Pe1 = 2.0 sinδ | 0.50 | 2.0 |
| 故障中 | Pe2 = 0.5 sinδ | 2.00 | 0.5 |
| 故障后 | Pe3 = 1.8 sinδ | 0.556 | 1.8 |
机械功率取Pm = 0.8,稳态运行点由故障前曲线决定:δ0 = arcsin(0.8/2.0) = 23.58°。故障期间运行点沿Pe2曲线运动,功角加速到切除角 δc 时断路器动作,之后运行点跳到Pe3曲线继续运动,直到转子在最大摆角 δmax 处动能回到零。最大摆角的极限位置是δmax = π - arcsin(Pm/Pmax3) = 153.57°,超过这个角度系统将无法恢复同步。
4.3 临界切除角的闭式解
令加速面积等于减速面积,可以解出临界切除角 δcr。面积表达式分别写出后合并同类项,得到闭式解:
cosδcr = [Pm(δmax - δ0) + Pmax3·cosδmax - Pmax2·cosδ0] / (Pmax3 - Pmax2)注意cosδmax在 δmax > 90° 时为负值,所以式中第二项实际上是正贡献。下面用 Python 直接代入参数计算:
import numpy as np Pm = 0.8 Pmax1 = 2.0 Pmax2 = 0.5 Pmax3 = 1.8 # 稳态功角与极限切除角对应的最大摆角(弧度) delta0 = np.arcsin(Pm / Pmax1) # 23.58° = 0.4115 rad deltamax = np.pi - np.arcsin(Pm / Pmax3) # 153.58° = 2.681 rad # 等面积准则闭式解 cos_dcr = (Pm * (deltamax - delta0) + Pmax3 * np.cos(deltamax) - Pmax2 * np.cos(delta0)) / (Pmax3 - Pmax2) deltacr = np.arccos(cos_dcr) print('delta0 =', round(np.degrees(delta0), 2), 'deg') print('deltamax=', round(np.degrees(deltamax), 2), 'deg') print('deltacr =', round(np.degrees(deltacr), 2), 'deg') print('cos(dcr)=', round(cos_dcr, 4))运行结果:δ0 = 23.58°、δmax = 153.58°、δcr = 101.3°。这个结果的含义是:只要故障切除时功角不超过 101.3°,系统就能保持暂态稳定;超过这个角度,即使后续减速面积全部用上也无法抵消加速面积。若计算出的cosδcr < -1,说明任何切除时刻都无法稳定;若cosδcr > 1或不存在解,说明故障严重程度较低,第一摆总能稳定。
等面积准则给出的是角度判据,不能直接回答 "多长时间内必须切除故障"。要把 δcr 换算成临界切除时间 tcr,必须回到摆动方程做数值积分,这就要用到下一章的数值解法。
5. 摆动方程数值解法的步长策略与故障时刻处理
5.1 一阶化与改进欧拉预测-校正
摆动方程是二阶非线性常微分方程,Pe(δ)是 δ 的正弦函数,无法写出解析解,只能数值积分。先把方程改写成一阶状态空间形式:
dδ/dt = ω dω/dt = ωs / (2H) · (Pm - Pe)其中Pe = Pmax2 · sinδ(故障期间)。初始条件为δ(0) = δ0、ω(0) = 0。传统暂态稳定程序常用改进欧拉法(预测-校正),每一步分两拍:
# 预测步:用当前点的导数外推一整步 y_predict = y + h * f(t, y) # 校正步:用预测点与当前点的导数平均值修正 y_next = y + h / 2 * (f(t, y) + f(t + h, y_predict))改进欧拉具有二阶精度,每步只需两次右端函数求值,在 H 常量 1~10 秒的系统中足够用。与标准欧拉法相比,同样的步长下半步误差小一个量级;与四阶龙格-库塔相比,计算量减半。首轮快速扫描时我会先用改进欧拉粗算,确认趋势后再用 RK4 出精确值。
5.2 四阶龙格-库塔实现与事件检测
下面用 RK4 积分故障期间的摆动方程,当功角越过 δcr 的时刻就是临界切除时间 tcr。之所以用 RK4,是因为它每步有四次函数求值,步长可以放到 0.01~0.02 秒而保持精度,整体效率反而更高。
import numpy as np Pm = 0.8 Pmax2 = 0.5 # 故障期间功率峰值 H = 5.0 # 惯性时间常数 f = 50.0 ws = 2 * np.pi * f deltacr = 1.768 # 上一节求出的临界切除角,101.3° 对应的弧度 delta0 = 0.4115 # 稳态功角 def rhs(y): """返回 [d(delta)/dt, d(omega)/dt]""" d, w = y Pe = Pmax2 * np.sin(d) return np.array([w, ws / (2 * H) * (Pm - Pe)]) def rk4_step(y, h): k1 = rhs(y) k2 = rhs(y + h / 2 * k1) k3 = rhs(y + h / 2 * k2) k4 = rhs(y + h * k3) return y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4) h = 0.01 t = 0.0 y = np.array([delta0, 0.0]) while y[0] < deltacr: y_new = rk4_step(y, h) # 事件检测:本步内功角跨过 deltacr 时线性插值 if y_new[0] >= deltacr: tcr = t + h * (deltacr - y[0]) / (y_new[0] - y[0]) break y, t = y_new, t + h print('tcr = %.4f s' % tcr)代码中rhs返回的是状态向量的导数,数组第一个元素是角速度,第二个元素是角加速度。ws / (2*H)是把标幺值功率差换算成加速度的系数,本算例中等于 31.416。循环体内的事件检测很关键:大步长下 δ 可能在一步之内从 99° 跳到 103°,如果不做插值,tcr 的误差会达到一个步长量级。线性插值用跨越步内两点相对位置的比例关系求出精确过零点,代价几乎为零。
对同一组参数,只改变步长 h,得到如下结果:
| 步长 h (s) | tcr (s) |
|---|---|
| 0.05 | 0.402 |
| 0.02 | 0.391 |
| 0.01 | 0.388 |
| 0.005 | 0.387 |
| 0.001 | 0.386 |
可以看到,当 h 从 0.05 缩小到 0.001,tcr 从 0.402 收敛到 0.386 秒左右。h=0.01 与 h=0.001 的差距只有 0.002 秒,相对误差约 0.5%。实际整定计算中,断路器开断时间通常是 0.06~0.1 秒,远小于 0.386 秒,说明这个场景有充足的安全裕度。
5.3 步长策略与故障切除时刻的工程经验
步长选择有一个粗略的参考标准:摆动周期与sqrt(2H / (ωs · dPe/dδ))有关,H=5 秒时首摆周期大约 1.5 秒,步长取周期的 1/100 以下即可满足精度,也就是 0.015 秒以下。但故障期间dPe/dδ很小甚至接近零,转子加速度大,建议保守地把步长压到 0.005 秒。另一个常见问题是功率曲线的跳变:故障切除瞬间Pe从Pmax2·sinδ突变到Pmax3·sinδ,如果直接把切除时刻对齐到步长节点,会引入相位误差。
提示:不要试图用平均功率曲线代替故障前后的功率跳变,这会同时低估加速面积和高估减速面积,得到偏乐观的 tcr。
正确的做法是把 tcr 作为已知参数代入,分段积分:故障阶段积分到 tcr,然后用故障后的功率方程继续积分,观察最大摆角是否小于 δmax。编程时注意先更新全部状态变量再进入下一步,避免把同步更新的微分方程组拆成异步更新。
6. 多机系统暂态分析中的网络降阶与失稳判据
6.1 多机系统的参考系与 COI 变换
多机系统不能像单机无穷大那样直接观察绝对功角,因为每台发电机的转子都在相对运动,系统没有一个固定的同步旋转轴。工程上常用惯量中心(COI)作为参考系:
δCOI = Σ(Hi·δi) / ΣHi ωCOI = Σ(Hi·ωi) / ΣHi每台发电机相对 COI 的功角定义为δi - δCOI。稳定性判断的标准从"功角是否越过 180°"变为"各机相对 COI 的功角是否持续增大且不可恢复"。COI 本身也在变速,所以转动方程中要额外计入 COI 的加速度项,这不是简单的坐标平移。
6.2 经典模型下的导纳矩阵降阶
多机系统中,发电机内节点通过输电网络互相耦合,网络节点规模远大于发电机节点。处理思路是用 Kron 降阶消去所有网络节点,只保留发电机内节点,把网络等值成各发电机之间的转移阻抗矩阵:
import numpy as np # Ygg: 发电机内节点自导纳和互导纳矩阵,维度 m×m # Ygn: 发电机节点与网络节点的互导纳,维度 m×n # Ynn: 网络节点导纳矩阵,维度 n×n # Yng: 网络节点与发电机节点的互导纳,维度 n×m Yred = Ygg - Ygn @ np.linalg.inv(Ynn) @ YngYred就是降阶后的等值导纳矩阵,维度为发电机台数 m。对角元包含各发电机的自导纳,非对角元就是机间转移导纳。故障期间的处理方式是在故障点追加一条对地导纳,修改Ynn后重新做一次降阶;故障切除后同理。这些矩阵运算全部在数值积分之前离线完成,积分过程中只更新 δ 和 ω,不需要在线修改导纳矩阵,这是经典法在速度上的优势。
6.3 用功角差判据快速验证数值结果
多机系统的失稳判据没有单机系统那么严格。工程上常用的经验门槛是看第一摆的最大相对功角差:
| 相对功角差(第一摆) | 工程判断 |
|---|---|
| < 90° | 首摆安全,裕度充足 |
| 90° ~ 135° | 临界区间,建议用等面积准则进一步校核 |
| > 135° | 大概率失步,需要加快切除或调整运行方式 |
这些数值不是严格定理,而是大量仿真结果的经验归纳。关于数值解本身,有一个容易被忽略的验证点:仿真开始前,用潮流结果反推各发电机内电势和初始功角,代入摆动方程后必须满足Pmi - Pei(0) = 0。如果不为零,说明初始点不是平衡点,数值解从一开始就在漂移,后续所有功角曲线都不可信。检查这一条,比检查任何高级判据都更优先。
本文还有配套的精品资源,点击获取