1、暂态稳定分析程序报告1 说明本程序是基于隐式梯形积分法的暂态稳定分析程序，采用 IEEE 3 机 9 节点模型，对其在 pt 时刻 7 号节点发生三相短路故障，并在ct 时刻将故障切除 (断开 5 号节点和 7 号节点之间的线路) 进行仿真，最后绘出 2 号和 1 号发电机的功角差( 21)与时间的关系曲线。发电机模型采用三阶模型，忽略阻尼，电枢电阻取为零；励磁系统仅有量测和放大环节，忽略反馈和电压调节器；负荷采用恒定阻抗负荷。程序是在 Matlab R2009b 上调试和运行的。2 模型推导及相关公式2.1发电机模型发电机采用三阶模型，其隐式梯形积分法推导如下：由发电机转子运动方程(1) d
2、T(P-)meJBwtA式中， =w-1，T J 和 t 单位为 s，w B=314rad/s，其他变量均为标幺wA值。根据梯形积分法则，得 tntn+1 时步的差分方程为(2)()(),1,1,2,()1hwPPnmnemneTJabwBab AAA其中(3),(),2,habPwwmneTJBaA由(2 )式，消去 ，其中 ，得到1nA,10mn(4),Pbwe式中 (5)()ba由发电机转子绕组暂态方程(6)()0TpEXIqdfd根据梯形积分法则差分化，得(7)(),1,1, ,1,20()2,hE IEnqqnf dnfdXIaEbqf 式中(8),()1210, , ,haXqqd
3、TdbEIEnqnf 由发电机定子电压方程(9)UXIrqaddI将其转化为 xy 坐标下为(10)(sincos)(cosin)0incosUXIIxyqxyEqd 又知 (11)PEIIeq将式(9) 代入式 (11)，并化为 xy 坐标得(12),111UIUIqxyentntnd将式(12 )代入式 (4)(13)aIbwt2.2励磁系统模型励磁系统传递函数框图如 图 1 所示：1pTR1AKTptUrefUf0E fTR RU图 1 励磁系统模型由 图 1 传递函数，可得励磁系统基本方程式为(14)=UT()0ptRRKTArefEf根据梯形积分法则，写出相应差分方程并整理得(15)
4、,1,1(),0UabTRUntnKAEfEfff式中(16),(2),2, ),habaUtnUTRTRTRTKKURnAAEffErefA 消去 UTR,n+1，U R,n+1，得(17),1,aUbtnffEf式中(18), 022,1,1,1,aKbKTRTRAAEffEfffUtnqnxynd 将式 (17)代入式(7)得到(19),(sincos)21EaUbqxyf EfIIqtn 联立式(13 )、 (10)和(19)，就得到系统的差分方程(略去下标 n+1)(20)0122,1,(sincos)i(sincos)03(i)4faIbwxyEUqf EfIIqf XIIxyxy
5、d 式中参数运用前面公式可求得。2.3负荷模型负荷模型看作是恒定阻抗，并入导纳矩阵。负荷的等值并联导纳为(21)*2SLYU其中 ， 是经过潮流计算的得到的负荷功率和负荷端节点电压。SLU2.4消去网络模型由于只有发电机端节点注入电流不为零，其他节点的注入电流全为零，因此可消去网络除发电机外的其他节点。消去过程如下：由于 (22)IYV其中 ，因此可以得到0In(23)0InrVY其中下标 n 代表发电机节点，r 代表其余节点。由公式(23 )展开并消去 Vr，得到(24)1()nYVIYn其中 就是消去网络节点后的导纳矩阵。Y消去网络节点后网络的方程如下：(25)56fIGUBxiijxij
6、yyji至此网络差分方程推导完毕，即式(20)和(25) 。对于 N 机系统有6N 个未知量，也有 6N 个方程，联立此方程运用牛顿迭代法即可求解。3 程序流程图本程序的主要流程图如下图所示：初始化导入数据和潮流计算结果初值计算引入故障并修改网络导纳和计算突变量形成雅可比矩阵并用求解差分方程根据修正量修改变量绘制发电机功角及功角差图T = T n ?T = c t ?切除故障 ， 修改网络导纳 、 计算突变量T = 0T = T + 1更新非定常参数YNNY结束T = p t ?是否失稳 ？NYNY图 2 暂态稳定分析程序流程图故障在 t=pt 时刻加入，并在 t=ct 时刻清除故障，如果系统
7、稳定令 t=t+1 继续进行牛顿法迭代，直到系统失稳，最后绘出发电机功角差与时间的关系曲线。4 变量说明LN：支路参数。格式如下：1 2 3 4 5 6首端节点 末端节点 电阻 r(pu) 电抗 x(pu) -B/2(或变比) 支路状态其中线路状态为 0 或 1。0 代表支路是线路，1 代表支路是变压器。GEN：发电机参数。格式如下：1 2 3 4发电机端节点 有功 P(pu) 无功 Q(pu)(或电压 U(pu) 发电机节点类型其中发电机节点类型为 0 或 1。0 代表节点是 PQ 节点，1 代表节点是 PV 节点。LOAD：负荷参数。格式如下：1 2 3负荷端节点 有功 P(pu) 无功
8、Q(pu)ROTOR：发电机暂态参数。格式如下：1 2 3 4 5 6 7节点名称 基准容量 漏抗 漏阻 D 轴电抗 D 轴暂态电抗D 轴次暂态电抗8 9 10 11 12 13 14D 轴开环时间常数D 轴次暂态时间常数Q 轴电抗Q 轴次暂态电抗Q 轴开环时间常数Q 轴次暂态时间常数惯性参数EXC：励磁系统参数。格式如下：1 2 3 4 5 6 7 8节点名称 TR KA T TB TC Efdmax EfdminAsp：定常参数。格式如下：1 2 3 4 5 6 7 8 9 10adita Pm0 aw aq1 aEf_1 aq2 aEf aUTR Uref Eq0a机械功率 wa1qEf
9、2qaEfUTRa励磁参考电压Q 初始轴电压nsp：非定常参数。格式如下：1 2 3 4 5 6 7IG0 Eq0_1 dita0 Pe0 UTR0 Ef0 Id0发电机电流Q 轴暂态电势发电机功角电磁功率励磁量测输出电压D 轴励磁电压D 轴电流8 9 10 11 12 13 14dw0 VR bdita_1 bEf_1 bq bUTR bEf旋转wA角速度励磁放大输出电压 bEfqbUTRbEf其中 ， ， ， ， ， ， ， ， ， ， , 是残aw1qEfa2qEfaUTREfqTREf差方程推导过程中产生的参数，可以用上述公式求出。其他变量定义如下：变量 含义 变量 含义pt 故障引入
10、时间 F 迭代残差ct 故障清除时间 deX 修正量slackbus 平衡节点 dita0，dita 存储发电机功角SG 发电机输出功率 ct 故障清除时间Base 系统基值 J 雅可比矩阵Un 节点电压 Yg 消去网络节点后的导纳矩阵Y 潮流计算的导纳矩阵fat 故障标志：fat=1 ，引入故障；fat=0，清除故障Yl 并入负荷后的导纳矩阵flag 失稳标志：flag=0 ，系统稳定；flag=1，系统失稳5 输出结果系统故障在 t=pt 时刻加入，故障是 7 号节点三相短路，在 ct 时刻将故障切除，断开 5 号和 7 号节点之间的支路。Matlab 仿真结果如下：0 0.5 1 1.5
11、 2 2.5-200-150-100-50050100150200time,sij degreesFig1 Plot of ij differences versus time(pt=0.05s and ct=0.167s) Fault cleared Fault introduced2131图 1 是当 pt=0.05s，ct=0.167s 时 2 号和 3 号发电机分别与 1 号发电机之间的功角差曲线图。由图可知，在 02.4s 的时间内，任何一个功角差均没有超过 180，系统是暂态稳定的；但 2.4s 以后就失去稳定了。0 0.2 0.4 0.6 0.8 1 1.2 1.4-200-15
12、0-100-50050100150200time,sij degreesFig2 Plot of ij differences versus time(pt=0.05s and ct=0.168s) Fault cleared Fault introduced2131图 2 是当 pt=0.05s，ct=0.168s 时 2 号和 3 号发电机分别与 1 号发电机之间的功角差曲线图。由图所示，系统在 2s 内就失去了稳定。进一步增加故障清除时间，得到图 3 所示图形，其中pt=0.05s，ct=0.4s。可以看出系统在清除故障前就已经失去了稳定。0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4-200-150-100-50050100150200time,sij degreesFig3 Plot of ij differences versus time(pt=0.05s and ct=0.4s) Fault cleared Fault introduced2131由此可以得出以下结论：系统发生故障时，随着故障切除时间
