简介电力系统分析是电气工程专业核心课程稳定性章节重点研究系统在扰动后维持同步运行的能力。课件深入讲解稳定性基本概念、摆动方程推导、稳态运行时同步电机模型并通过等面积准则分析暂态稳定同时涵盖三相短路故障影响、非线性方程及摆动方程的数值解法最后扩展到多机系统与多机暂态分析理论体系完整。压缩包内共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串联一个直轴暂态电抗Xd。这个模型的成立依据是扰动发生后的极短时间内励磁绕组磁链近似守恒因此暂态电动势E基本不变。经典模型在首摆稳定性分析中非常可靠但它有明确的适用边界分析场景经典模型适用性失效原因近区三相短路首摆分析适用励磁绕组磁链守恒E 保持恒定多摆过程与动态稳定不适用励磁调节器、阻尼绕组开始起主导作用新能源场站动态不适用变流器没有转子运动方程输出功率完全受控制决定经典模型对单机无穷大系统的描述是把发电机内节点看成一个电势源经过Xd、变压器电抗、线路电抗连接到电压恒定的母线上。之所以叫无穷大母线是因为它的电压幅值和频率在任意功率注入或抽取下都不变化。这个最小系统足以说清楚稳定性的所有基本概念从单机系统理解透之后多机系统的扩展只是把单机方程写成向量形式。2.3 双总线方程与功角曲线的数值绘制把发电机内节点记为节点 1无穷大母线记为节点 2用节点导纳矩阵写出双母线系统忽略所有电阻后发电机输送到无穷大母线的功率可以化成一个极简形式P_e E · V / X · sinδX 是暂态电抗后的总转移电抗。这个方程说明功率传输取决于两个因素转移电抗和两电压之间的夹角。P_e随 δ 按正弦变化峰值Pmax EV/X出现在 δ 90° 处这就是稳态功率极限。下面用例 11.1 的参数绘制暂态功角曲线同时对比忽略凸极与考虑凸极两种模型的差别——后文会详细推导这两组参数是如何计算出来的。import numpy as np import matplotlib.pyplot as plt # 例 11.1 的参数Xd1.0, Xq0.6, Xd0.3, 机端电压 V1.0 Xdp 0.3 Xq 0.6 V 1.0 # 忽略凸极E1.1226考虑凸极E1.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, labelround-rotor (ignore saliency)) plt.plot(np.degrees(delta), P_salient, labelsalient-pole) plt.axhline(0.5, ls--, colorgray) 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 给出了一个典型的单机无穷大系统算例。发电机参数标幺值为Xd1.0、Xq0.6、Xd0.3机端电压V1.0输出有功P0.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 jXd 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 jXd 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 (Xd·E (Xd - Xd)·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(δ) EV/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.12263.7419 sinδ3.741990.0°考虑凸极1.11623.7208 sinδ - 0.8333 sin2δ4.032110.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.502.0故障中Pe2 0.5 sinδ2.000.5故障后Pe3 1.8 sinδ0.5561.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.050.4020.020.3910.010.3880.0050.3870.0010.386可以看到当 h 从 0.05 缩小到 0.001tcr 从 0.402 收敛到 0.386 秒左右。h0.01 与 h0.001 的差距只有 0.002 秒相对误差约 0.5%。实际整定计算中断路器开断时间通常是 0.06~0.1 秒远小于 0.386 秒说明这个场景有充足的安全裕度。5.3 步长策略与故障切除时刻的工程经验步长选择有一个粗略的参考标准摆动周期与sqrt(2H / (ωs · dPe/dδ))有关H5 秒时首摆周期大约 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。如果不为零说明初始点不是平衡点数值解从一开始就在漂移后续所有功角曲线都不可信。检查这一条比检查任何高级判据都更优先。本文还有配套的精品资源点击获取