《MATLAB实验电力系统暂态稳定分析.pdf》由会员分享，可在线阅读，更多相关《MATLAB实验电力系统暂态稳定分析.pdf（12页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、.实验三 电力系统暂态稳定分析 电力系统暂态稳定计算实际上就是求解发电机转子运动方程的初值问题，从而得出-t和-t 的关系曲线。每台发电机的转子运动方程是两个一阶非线性的常微分方程。因此，首先介绍常微分方程的初值问题的数值解法。一、常微分方程的初值问题（一）问题及求解公式的构造方法 我们讨论形如式（3-1）的一阶微分方程的初值问题 00)(),()(yxybxayxfxy （3-1）设初值问题（3-1）的解为)(xy，为了求其数值解而采取离散化方法，在求解区间ba,上取一组节点 bxxxxxanii110 称iiixxh1（1,1,0ni）为步长。在等步长的情况下，步长为 nabh 用iy表示
2、在节点ix处解的准确值)(ixy的近似值。设法构造序列 iy所满足的一个方程（称为差分方程）),(1hyxhyyiiii （3-2）作为求解公式，这是一个递推公式，从（0 x，0y）出发，采用步进方式，自左相右逐步算出)(xy在所有节点ix上的近似值iy（ni,2,1）。在公式（3-2）中，为求1iy只用到前面一步的值iy，这种方法称为单步法。在公式（3-2）中的1iy由iy明显表示出，称为显式公式。而形如（3-3）),(11hyyxhyyiiiii （3-3）的公式称为隐式公式，因为其右端中还包括1iy。如果由公式求1iy时，不止用到前一个节点的值，则称为多步法。由式（3-1）可得 dy=d
3、xyxf),(（3-4）两边在ix，1ix上积分，得.1)(,()()(1iixxiidxxyxfxyxy （3-5）由此可以看出，如果想构造求解公式，就要对右端的积分项作某种数值处理。这种求解公式的构造方法叫做数值积分法。（二）一般的初值问题的解法 1 欧拉法和改进欧拉法 对于初值问题（3-1），采用数值积分法，从而得到（3-5）。对于（3-5）右端的积分用矩形公式(取左端点)，则得到 1)(,()(,(iixxiixyxfhdxxyxf 进而得到（3-1）的求解公式（3-2）),(1iiiiyxfhyy （i=0，1，2，n-1）（3-6）此公式称为欧拉（Euler）格式。如果对式（3-5
4、）右端的积分用梯形公式)(,()(,(2)(,(111iixxiixyxfxyxfhdxxyxfii 则可以得到初值问题（3-1）的梯形求解公式如式（3-7）),(),(2111iiiiiiyxfyxfhyy （i=0，1，2，n-1）（3-7）式（3-7）是个隐式公式。可以采取先用欧拉格式求一个)(1ixy的初步近似值，记作1iy，称之为预报值，然后用预报值1iy替代式（3-7）右端的1iy，再计算得到1iy，称之为校正值，这样建立起来的预报校正方法称为改进欧拉格式),(),(2),(1111iiiiiiiiiiyxfyxfhyyyxfhyy （3-8）2 龙格库塔方法 在单步法中，应用最广
5、泛的是龙格库塔（Runge-kutta）法，简称 RK 法。下面直接给出一种四阶的龙格库塔法的计算公式（3-9）.),()21,2()21,2(),()22(61342312143211KyhxfhKKyhxfhKKyhxfhKyxfhKKKKKyyiiiiiiiiii （3-9）它也称为标准（古典）龙格库塔法。例 3-1 研究下列微分方程的初值问题 0)0(21122yyxy 解：这是一个特殊的微分方程，其解的解析式可以给出，为 21xxy 应用龙格库塔法，取h=0.25，根据式（3-9）编写一段程序，由零开始自左相右逐步算出)(xy在所有节点ix上的近似值iy。计算结果见表 3-1。计算结
6、果表明，四阶龙格库塔方法的精度是较高的。表 3-1 nx ny nnyxy)(2.0 0.39995699 4.3e-5 4.0 0.23529159 2.5e-6 6.0 0.16216179 3.7e-7 8.0 0.12307683 9.2e-8 实际上，MATLAB 为常微分方程提供了很好的解题指令，使得求解常微分方程变得很容易，并且能将问题及解答表现在图形上。因此，我们可以不用根据式（3-9）编写较复杂的程序，而只需应用 MATLAB 提供的常微分方程解题器来解决问题。下面给出用 MATLAB编写的解题程序。首先编写描述常微分方程的 ODE 文件，文件名为myfun，便于解题器调用它
7、。function dy=myfun(x,y)dy=zeros(1,1);dy=1/(1+x2)-2*y2;再编写利用解题器指令求解y的程序。clear x0=0;for i=1:4 xm=2*i;.y0=0;x,y=ode45(myfun,x0 xm,y0);format long y(length(y)end plot(x,y,-)运行上述程序，在得到几个点的函数值的同时，也得到函数 y 的曲线，如图 3-1 所示。0246800.10.20.30.40.5xy 图 3-1 根据运算结果画出 y 的曲线 二、简单电力系统的暂态稳定性（一）物理过程分析 某简单电力系统如图 3-2(a)所示，
8、正常运行时发电机经过变压器和双回线路向无限大系统供电。发电机用电势E作为其等值电势，则电势E与无限大系统间的电抗为 212TLTdxxxxx （3-10）这时发电机发出的电磁功率可表示为 sinsinMPxUEP （3-11）如果突然在一回输电线路始端发生不对称短路，如图 3-2(b)所示。故障期间发电机电势E与无限大系统之间的联系电抗为 xxxxxxxxxxTLTdTLTd)2)()2()(2121 （3-12）在故障情况下发电机输出的电磁功率为 sinsinMPxUEP （3-13）在短路故障发生之后，线路继电保护装置将迅速断开故障线路两端的断路器，如图 3-2(c)所示。此时发电机电势E
9、与无限大系统间的联系电抗为.21TLTdxxxxx （3-14）发电机输出的功率为 sinsinMPxUEP （3-15）U=cE(a)GT1T2LLjxLjx1Tjx2Tjxdx j UE(b)LjxLjx1Tjx2Tjxdx j UjxE(c)Ljx1Tjx2Tjxdx j U 图 3-2 简单电力系统及其等值电路（a）正常运行方式及其等值电路；（b）故障情况及其等值电路；（c）故障切除后及其等值电路 如果正常时发电机向无限大系统输送的有功功率为0P，则原动机输出的机械功率TP等于0P。假定不计故障后几秒种之内调速器的作用，即认为机械功率始终保持0P。因此，可以得到此简单电力系统正常运行、
10、故障期间及故障切除后的功率特性曲线如图 3-3 所示。0kcmh0PPTPPPPakh 图 3-3 简单系统正常运行、故障期间及故障切除后的功率特性曲线 对于上述简单电力系统，我们可以根据等面积定则求得极限切除角。但是，实际工作需要知道在多少时间之内切除故障线路，也就是要知道与极限切除角对应的极限切除时间。要解决这个问题，必须求解发电机的转子运动方程。（二）求解发电机的转子运动方程.求解发电机转子运动方程可以得出-t 和-t 的关系曲线。其中-t 曲线一般称为摇摆曲线。在上述简单电力系统中故障期间的转子运动方程为)sin(1)1(1MTJPPTdtddtd （3-16）式中，功率角，其单位为弧
11、度；转子角速度，标幺值；1转子的同步角速度，即1=f2=314.16，其单位为弧度/秒；JT发电机的惯性时间常数，其单位为秒；TP、MP分别为机械和电磁功率，标幺值。这是两个一阶的非线性常微分方程，它的起始条件是已知的，即 t=0t=0；=0=1.0；=0=MTPP1sin 故障切除后，由于系统参数改变，以致发电机功率特性发生变化，必须开始求解另一组微分方程：)sin(1)1(1MTJPPTdtddtd （3-17）式中变量含义同前述，其中MP也为标幺值。这组方程的起始条件为 t=ct；=c；=c 其中ct为给定的切除时间；c、c为与ct时刻对应的和，它们可由故障期间的-t 和-t 的关系曲线
12、求得（和都是不突变的）。一般来说，在计算故障发生后几秒种的过程中，如果始终不超过 180，而且振荡幅值越来越小，则系统是暂态稳定的。当发电机与无限大系统之间发生振荡或失去同步时，在发电机的转子回路中，特别是阻尼绕组中将有感应电流而形成阻尼转矩（也称为异步转矩）。当作微小振荡时，阻尼功率可表达为：DP=D=)1(D （3-18）式中，D称为阻尼功率系数；为转子角速度的偏移量，标幺值；为转子角速度，标幺值。阻尼功率系数D除了与发电机的参数有关外，还和原始功角、的振荡频率有关。在一般情况下它是正数。在原始功角较小，或者定子回路中有串联电容使定子回路总电阻相对于总电抗较大时，D 可能为负数。如果考虑阻
13、尼功率的影响，则故障后的转子运动方程又可表达为.sin)1(1)1(1MTJPDPTdtddtd （3-19）电力系统暂态稳定计算包括两类问题，一类是应用数值计算法得出故障期间的曲线后，根据曲线找到与极限切除角对应的极限切除时间，此时只需要求解微分方程（3-16）；另一类是已知故障切除时间，需要求出摇摆曲线来判断系统的稳定性，此时需要分段分别求解微分方程（3-16）和（3-17）。如果考虑阻尼转矩的影响，则此时需要分段分别求解微分方程（3-16）和（3-19）。三、例题 例 3-2 某简单电力系统如图 3-4 所示，取基准值BS=220MVA，BU=209KV。换算后的参数已经标在图中，其中一
14、回线的电抗Lx=0.486，JT=8.18 秒。设电力线路某一回的始端发生两相接地短路。假定E=常数。（1）计算保持暂态稳定而要求的极限切除角。（2）计算极限切除时间，并且作出在 0.15 秒切除故障时的-t 曲线。GT1T2Lkdx1TxLx2Tx0P0QU=0.122=1.0=0.486=0.138=0.432=0.295LLxx40=1.0=0.22x 图 3-4 某简单电力系统的接线图 解：计算系统正常运行方式，决定E和0。由 3-3(a)的正序网络可得，此时系统的总电抗为 x=0.295+0.138+0.243+0.122=0.798 发电机的暂态电势为：E=22)0.1798.00
15、.1()0.1798.02.00.1(=1.41 0=1tg798.02.00.1798.0=34.53（2）故障后的功率特性 又由 3-3(b)的负序、零序网络可得故障点的负序、零序等值电抗为 2x=)122.0243.0()138.0432.0()122.0243.0()138.0432.0(=0.222.0 x=)122.0972.0(138.0)122.0972.0(138.0=0.123 所以在正序网络故障点上的附加电抗为：079.0123.0222.0123.0222.0 x 于是故障时等值电路如图 3-3(c)所示，则 80.2079.0365.0433.0365.0433.0
16、x 因此，故障期间发电机的最大功率为：504.08.20.141.1xUEPM（3）故障切除后的功率特性 故障切除后的等值电路如图 3-3(d)所示 041.1122.0486.0138.0295.0 x 此时最大功率为 35.1041.10.141.1xUEPM 0102.13235.10.1sin180h j 0.295j 0.138j 0.243j 0.122U=1.00.10P2.00Q41.1E0053.34j 0.138j 0.243j 0.122j 0.432j 0.138j 0.1224243.0jj 0.122j 0.243j 0.138j 0.295j 0.295j 0.1
17、38j 0.486j 0.122)1(f(a)2(f(b)0(f(c)(d41.1E41.1Ej 0.079U=1.0U=1.0 图 3-5 例題 7-12 的等值电路（a）正常运行等值电路；（b）负序和零序等值电路；（c）故障时等值电路；（d）故障切除后等值电路（4）计算极限切除角.MMMhMhTcmPPPPP00coscos)(cos=504.035.153.34cos504.02.132cos35.1)53.342.132(1800.100=0.458 cm074.62（5）找出极限切除时间cmt 根据（3-16），首先计算初值 6027.018053.340，0.10 令 y(1)=，
18、y(2)=。编写描述故障期间转子运动方程的 ODE 文件，文件名为myequ。function dy=myequ(t,y)dy=zeros(2,1);f=50;w1=2*pi*f;dy(1)=(y(2)-1)*w1;dy(2)=(1/8.18)*(1.0-0.504*sin(y(1);再编写利用解题器指令求解 y 的程序。clear t0=0;tm=0.25;d0=(34.53/180)*pi;w0=1;T,Y=ode45(myequ,t0 tm,d0 w0);plot(T,(Y(:,1)/pi)*180,-,0.194,62.76,*)text(0.194,60,delta_cmax=62.
19、76circ,FontSize,10)text(0.194,56,t_cmax=0.194s,FontSize,10)00.050.10.150.20.2530405060708090?图 3-6 例題 7-12 的-t 曲线 图 3-6 给出短路发生后 0 秒到 0.25 秒期间的-t 计算曲线，根据最大切除角cm.（074.62）找到极限切除时间cmt为 0.194 秒。由图 3-6 可见，如果故障切除时间大于 0.194秒，则发电机的功角将不断地增大，最终失去暂态稳定。在极限切除时间之前切除故障，发电机的摇摆曲线的状况将在下面作计算、分析。（6）不考虑阻尼转矩的影响，当故障切除时间为 0
20、.15 秒时通过计算得出-t 曲线 首先编写描述故障期间转子运动方程的 ODE 文件，文件名为”myfun01”。function dy=myfun01(t,y)f=50;w1=2*pi*f;TJ=8.18;Pt=1.0;P2m=0.504;dy=zeros(2,1);dy(1)=(y(2)-1)*w1;dy(2)=(1/TJ)*(Pt-P2m*sin(y(1);再编写描述故障切除后转子运动方程的 ODE 文件，文件名为”myfun02”。function dy=myfun02(t,y)f=50;w1=2*pi*f;TJ=8.18;Pt=1.0;P3m=1.35;dy=zeros(2,1);d
21、y(1)=(y(2)-1)*w1;dy(2)=(1/TJ)*(Pt-P3m*sin(y(1);编写利用解题器指令求解y的小程序。clear t0=0;tc=0.15;tm=2.0;d0=(34.53/180)*pi;w0=1.0;T1,Y1=ode45(myfun01,t0 tc,d0 w0);dc=Y1(length(Y1),1);wc=Y1(length(Y1),2);T2,Y2=ode45(myfun02,tc tm,dc wc);plot(T1,(Y1(:,1)/pi)*180,-,T2,(Y2(:,1)/pi)*180,-,tc,(dc/pi)*180,*)text(0.28,50,
22、itt_c=0.15s,FontSize,8)text(0.28,43,itdelta_c=51.71circ,FontSize,8)xlabel(itt)ylabel(itdelta)计算结果表明，功角沿着故障切除后的功角特性曲线根据等面积定则作等幅振荡，如图 3-7 所示。实际上，由于阻尼转矩的影响，振荡的幅度是逐渐衰减的，功角最终运行在k=47.8。因此，发电机能够保持暂态稳定。.00.511.52020406080100?图 3-7 不考虑阻尼转矩影响，当 0.15 秒切除故障时发电机的t曲线（7）考虑阻尼转矩的影响，当故障切除时间为 0.15 秒时通过计算得出-t 曲线 描述故障期间
23、转子运动方程的 ODE 文件与（6）相同，文件名也为”myfun01”。重新编写描述故障切除后转子运动方程的 ODE 文件，文件名为”myfun03”，阻尼功率系数 D 取为 15。function dy=myfun03(t,y)f=50;w1=2*pi*f;TJ=8.18;Pt=1.0;P3m=1.35;D=15;dy=zeros(2,1);dy(1)=(y(2)-1)*w1;dy(2)=(1/TJ)*(Pt-D*(y(2)-1)-P3m*sin(y(1);再编写利用解题器指令求解y的小程序。clear t0=0;tc=0.15;tm=4;d0=34.53*3.14/180;w0=1;T1,
24、Y1=ode45(myfun01,t0 tc,d0 w0);dc=Y1(length(Y1),1);wc=Y1(length(Y1),2);T2,Y2=ode45(myfun03,tc tm,dc wc);plot(T1,(Y1(:,1)/pi)*180,-,T2,(Y2(:,1)/pi)*180,-,tc,(dc/pi)*180,*)text(0.3,50,itt_c=0.15s,FontSize,8)text(0.28,43,itdelta_c=51.71circ,FontSize,8)xlabel(itt)ylabel(itdelta).0123430405060708090?图 3-8 不考虑阻尼转矩影响，当 0.15 秒切除故障时发电机的t曲线 计算结果表明，功角沿着故障切除后的功角特性曲线作减幅振荡，如图 3-8 所示。功角最终运行在k=47.8。因此，发电机能够保持暂态稳定。