课程设计报告

( 2011—20xx年度第二学期)

名 称：电力系统暂态上机 院 系：电气与电子工程学院 班 学 号： 1091420228 学生姓名： 雷 刚 指导教师： 麻 秀 范 老 师 设计周数：成 绩：

日期： 20xx年 7月3

日

课程 电力系统暂态上机课程设计报告 一、课程设计的目的与要求

巩固电力系统暂态分析理论知识，使学生掌握采用计算机对电力系统电磁暂态过程和机电暂态过程进行计算的基本方法，并进一步巩固计算机编程能力，为将来从事相关的技术工作打下必要的基础。

二、设计正文（详细内容见附录,用A4纸，页数不限）

1． 对称短路计算过程流程图和计算结果

2． 不对称短路计算过程流程图和计算结果

3． 静态稳定性计算过程流程图和计算结果

4． 暂态稳定性计算过程流程图和计算结果

5． 思考题

三、课程设计总结或结论

在这次电力系统暂态上机课程中，我们主要讨论了各种类型短路故障下（三相短路、单相接地短路、两相相间短路、两相接地短路、单相经电阻接地）系统网络中的电压、电流分布的计算以及电力系统的静态稳定、暂态稳定问题的分析。

通过本次课程设计，我对之前学习的知识有了一个巩固和提高，同时对电力系统故障分析有了更深刻的理解。电力系统的故障时，大部分电磁量将随时间变化，描述其特性的是微分方程，这给分析计算带来一定困难。在分析过程中通常尽量避免对微分方程直接求解，而是采用一定的工具和假设使问题得以简化，即把“微分方程代数化，暂态分析稳态化”。在分析不对称故障时，各相之间电磁量的耦合使问题的分析更为复杂，此时常用的分析方法是采用对称分量法将不对称问题转化为对称问题来求解。同时我对用来分析电力系统静态稳定的试探法，用来分析电力系统暂态稳定的改进欧拉法有了一些使用心得。这与手算系统短路电流时使用的网络化简方法大大不同。同时也明白了计算机编程方法对于电力系统稳态暂态计算的简便性、快速性和重要性。

最后，经过此次课程设计，对于我之前在数学建模中的学习到的Matlab知识起到了升华的作用，可以说这也算是一次比赛把，对较复杂的选择语句、循环语句的使用等有了更深的了解，为今后学习、工作中埋下一个伏笔，相信会受益匪浅。

四、参考文献

1.

2.

3.

4.

5.

6.

7.

1 《电力系统暂态分析》，李光琦，中国电力出版社，20xx年，第三版； 《电力系统分析》（上、下），何仰赞，华中科技大学出版社，19xx年，第二版。 《电力系统故障的计算机辅助分析》 重庆大学出版社 米麟书等 《电力系统潮流计算》 天津大学出版社 宋文南等 《电力系统故障分析》 清华大学出版社 周荣光 《短路电流实用计算方法》 电力工业出版社 西安交通大学等 《精通Matlab6.5》北京航空航天大学出版社，张志涌

课程 电力系统暂态上机课程设计报告 附录（设计流程图、计算结果、思考题答案）

1. 对称短路计算过程流程图和计算结果

流程图如下：

计算结果如下：

(1)导纳矩阵:

2

课程 电力系统暂态上机课程设计报告

Y1 =Y2=

0 -26.6667i 0 +10.0000i 0 +10.0000i

0 +10.0000i 0 -33.3333i 0 +10.0000i 0 +10.0000i 0 +10.0000i 0 -20.0000i Y0=

0 -29.0476i 0 + 5.0000i 0 + 5.0000i

0 + 5.0000i 0 -48.0228i 0 + 5.0000i 0 + 5.0000i 0 + 5.0000i 0 -10.0000i YY1=YY2=

0 -40.0000i 0 +10.0000i 0 +10.0000i 0 +20.0000i 0 0 +10.0000i 0 -60.0000i 0 +10.0000i 0 0 +40.0000i 0 +10.0000i 0 +10.0000i 0 -20.0000i 0 0 0 +20.0000i 0 0 0 -30.0000i 0 0 0 +40.0000i 0 0 0 -60.0000i

(2) 节点3发生三相短路故障 ?

故障点三序电流、三相电流、三相电流有效值：

I1= 0 - 9.8592i I2=0 I3=0 ;

Ia= 0 - 9.8592i Ib=-8.5383 + 4.9296i Ic=8.5383 + 4.9296i Ia_effective = Ib_effective = Ic_effective = 9.8592 ?

各节点三序电压、三相电压：

3

课程 电力系统暂态上机课程设计报告

2．不对称短路计算过程流程图和计算结果

不对称短路计算过程流程图如三相短路流程图。

(1) 节点3发生A相短路接地故障 ?

故障点三序电流、三相电流：

I1= 0 - 3.1239i I2= 0 - 3.1239i I3= 0 - 3.1239i ; Ia= 0 - 9.3718i Ib= 0 - 0.0000i Ic= 0 - 0.0000i ； ? 各节点三序电压、三相电压：

(2) 节点3发生A相经10Ω电阻接地故障 ? 故障点三序电流、三相电流：

I1= 0.0333 - 0.0004i I2= 0.0333 - 0.0004i I3= 0.0333 - 0.0004i ; Ia= 0.1000 - 0.0011i Ib= 0 Ic= 0 ；

4

课程 电力系统暂态上机课程设计报告

(3) 节点3发生b、c两相短路故障 ? 故障点三序电流、三相电流：

I1= 0 - 4.9697i I2= 0 + 4.9697i I3= 0 ;

Ia= 0 Ib= -8.6078 - 0.0000i Ic= 8.6078 + 0.0000i；

?

各支路三序电流、三相电流：

(4) 节点3发生b、c两相短路接地故障 ? 故障点三序电流、三相电流：

I1= 0 - 6.4473i I2= 0 + 3.4921i I3= 0 + 2.9552i ; Ia= 0 Ib= -8.6078 + 4.4327i Ic= 8.6078 + 4.4327i ；

3．静态稳定性计算过程流程图和计算结果

5

课程 电力系统暂态上机课程设计报告 静态稳定性计算过程流程图如下：

6

课程 电力系统暂态上机课程设计报告

计算结果：

Ke=0.100000,delta=91,P=1.345117,Eq=1.999966 Ke=0.200000,delta=92,P=1.363000,Eq=2.028274 Ke=0.300000,delta=93,P=1.380693,Eq=2.056939 Ke=0.400000,delta=95,P=1.398409,Eq=2.089197 Ke=0.500000,delta=96,P=1.415747,Eq=2.119390 Ke=0.600000,delta=97,P=1.432842,Eq=2.149961 Ke=0.700000,delta=98,P=1.449679,Eq=2.180919 Ke=0.800000,delta=99,P=1.466244,Eq=2.212270 Ke=0.900000,delta=101,P=1.482360,Eq=2.251089 Ke=1.000000,delta=102,P=1.498265,Eq=2.283967 Ke=1.100000,delta=103,P=1.513819,Eq=2.317233 Ke=1.200000,delta=104,P=1.529000,Eq=2.350886 Ke=1.300000,delta=105,P=1.543786,Eq=2.384924 Ke=1.400000,delta=106,P=1.558152,Eq=2.419341 Ke=1.500000,delta=107,P=1.572075,Eq=2.454133 Ke=1.600000,delta=108,P=1.585527,Eq=2.489291 Ke=1.700000,delta=109,P=1.598483,Eq=2.524806 Ke=1.800000,delta=108,P=1.613666,Eq=2.534328 Ke=1.900000,delta=105,P=1.626492,Eq=2.515421 Ke=2.000000,delta=103,P=1.635259,Eq=2.507344 Ke=2.100000,delta=100,P=1.637025,Eq=2.483645 Ke=2.200000,delta=98,P=1.637979,Eq=2.471608 Ke=2.300000,delta=96,P=1.635815,Eq=2.457976 Ke=2.400000,delta=94,P=1.630657,Eq=2.442916 Ke=2.500000,delta=92,P=1.622630,Eq=2.426585 Ke=2.600000,delta=90,P=1.611867,Eq=2.409132 Ke=2.700000,delta=88,P=1.598501,Eq=2.390698 Ke=2.800000,delta=86,P=1.582670,Eq=2.371416 Ke=2.900000,delta=84,P=1.564509,Eq=2.351412 Ke=3.000000,delta=82,P=1.544155,Eq=2.330805

Ke=3.100000,delta=81,P=1.535174,Eq=2.323407 Ke=3.200000,delta=79,P=1.511558,Eq=2.301767 Ke=3.300000,delta=77,P=1.486103,Eq=2.279814 Ke=3.400000,delta=76,P=1.474144,Eq=2.271027 Ke=3.500000,delta=74,P=1.445989,Eq=2.248507 Ke=3.600000,delta=73,P=1.432473,Eq=2.239073 Ke=3.700000,delta=71,P=1.401945,Eq=2.216241 Ke=3.800000,delta=70,P=1.387071,Eq=2.206335 Ke=3.900000,delta=69,P=1.371732,Eq=2.196226 Ke=4.000000,delta=67,P=1.338431,Eq=2.173194 Ke=4.100000,delta=66,P=1.321987,Eq=2.162830 Ke=4.200000,delta=65,P=1.305171,Eq=2.152344 Ke=4.300000,delta=63,P=1.269626,Eq=2.129492 Ke=4.400000,delta=62,P=1.251928,Eq=2.118929 Ke=4.500000,delta=61,P=1.233937,Eq=2.108312 Ke=4.600000,delta=60,P=1.215672,Eq=2.097660 Ke=4.700000,delta=59,P=1.197153,Eq=2.086987 Ke=4.800000,delta=58,P=1.178398,Eq=2.076310 Ke=4.900000,delta=56,P=1.139925,Eq=2.054329 Ke=5.000000,delta=55,P=1.120623,Eq=2.043828 Ke=5.100000,delta=54,P=1.101144,Eq=2.033369 Ke=5.200000,delta=53,P=1.081504,Eq=2.022963 Ke=5.300000,delta=52,P=1.061715,Eq=2.012622 Ke=5.400000,delta=51,P=1.041792,Eq=2.002355 Ke=5.500000,delta=50,P=1.021746,Eq=1.992173 Ke=5.600000,delta=49,P=1.001589,Eq=1.982085 Ke=5.600000,delta=49,P=1.001589,Eq=1.972100

最终选择放大倍数Ke=2.2

7

课程 电力系统暂态上机课程设计报告

4．暂态稳定性计算过程流程图和计算结果

暂态稳定性计算过程流程图如下：

1

课程 电力系统暂态上机课程设计报告 计算结果如下：

0.15s时切除故障的摇摆曲线

Delta

t 00.20.40.60.811.21.41.61.82

0.15s时切除故障的摇摆曲线 0.25s时切除故障的摇摆曲线

根据摇摆曲线判断：0.15s时切除故障系统暂态稳定，0.25s时切除故障系统的功角无限增大，系统失去稳定性。

另据试探，0.2s时系统临界稳定，为保证系统暂态稳定，切除时间最大值约为0.2s。

2

课程 电力系统暂态上机课程设计报告

程序：

Matrix导纳阵计算程序

clc;

clear;

%---------------输入已知条件----------%

bus_Num1=3; %节点数

bus_Num2=5; %包括发电机节点的节点数

branch_Num1=3; %线路数

branch_Num2=5; %包括发电机支路的支路数

branch1_No1=[1,1,2];%支路首节点

branch1_No2=[2,3,3];%支路末节点

branch2_No1=[1,1,2,1,2];

branch2_No2=[2,3,3,4,5];

%输入支路各序阻抗，z1_branch表示支路正序阻抗，z2_branch表示支路负序阻抗，z0_branch表示支路零序阻抗

z1_branch(1)=j*0.1;

z1_branch(2)=j*0.1;

z1_branch(3)=j*0.1;

z2_branch=z1_branch;

z0_branch(1)=j*0.2;

z0_branch(2)=j*0.2;

z0_branch(3)=j*0.2;

%---------------第一步：不考虑发电机节点计算节点导纳矩阵----------%

3

课程 电力系统暂态上机课程设计报告 %节点导纳矩阵，Y1表示不计发电机节点的正序网络节点导纳阵，Y2表示不计发电机节点的负序网络节点导纳阵，Y0表示不计发电机节点的零序网络节点导纳阵，

Y1=zeros(bus_Num1);

Y1(1,1)=1/(j*0.15);

Y1(2,2)=1/(j*0.075);

%请同学们求正序节点导纳矩阵

for i=1:branch_Num1

y=1/z1_branch(i);

ii=branch1_No1(i);

jj=branch1_No2(i);

Y1(ii,ii)=Y1(ii,ii)+y;

Y1(ii,jj)=Y1(ii,jj)-y;

Y1(jj,ii)=Y1(jj,ii)-y;

Y1(jj,jj)=Y1(jj,jj)+y;

end

Y2=Y1;%负序等于正序

Y0=zeros(bus_Num1);

Y0(1,1)=1/(j*0.0525);

Y0(2,2)=1/(j*0.0263);

%请同学们求零序节点导纳矩阵

for i=1:branch_Num1

y=1/z0_branch(i);

ii=branch1_No1(i);

jj=branch1_No2(i);

Y0(ii,ii)=Y0(ii,ii)+y;

Y0(ii,jj)=Y0(ii,jj)-y;

Y0(jj,ii)=Y0(jj,ii)-y;

Y0(jj,jj)=Y0(jj,jj)+y;

end

4

课程 电力系统暂态上机课程设计报告 %---------------第二步：考虑发电机节点计算节点导纳矩阵----------%

%节点导纳矩阵，YY1表示计及发电机节点的正序网络节点导纳阵，YY2表示计及发电机节点的负序网络节点导纳阵

z1_branch(4)=j*0.05;

z1_branch(5)=j*0.025;

z2_branch=z1_branch;

YY1=zeros(bus_Num1);

YY1(4,4)=1/(j*0.1);

YY1(5,5)=1/(j*0.05);

%请同学们求正序节点导纳矩阵

for i=1:branch_Num2

y=1/z1_branch(i);

ii=branch2_No1(i);

jj=branch2_No2(i);

YY1(ii,ii)=YY1(ii,ii)+y;

YY1(ii,jj)=YY1(ii,jj)-y;

YY1(jj,ii)=YY1(jj,ii)-y;

YY1(jj,jj)=YY1(jj,jj)+y;

end

YY2=YY1;%负序等于正序

%---------------第三步：计算节点阻抗矩阵----------%

Z1=inv(Y1);

Z2=inv(Y2);

Z0=inv(Y0);

ZZ1=inv(YY1);

ZZ2=inv(YY2);

5

课程 电力系统暂态上机课程设计报告

短路电流计算程序shortcircuit

clc;

clear;

%数据来源于教材《电力系统暂态分析》P77例（3-4），P143例（5-7）

Matrix %计算节点导纳矩阵、节点阻抗矩阵，形成全局变量

Fault_Node=input('输入短路点编号 ;\n Fault_Node=');

Fault_Type=input('输入短路类型 ;\n(1)Fault_Type=0为三相短路;\n(2)Fault_Type=1为a相接地短

路;\n(3)Fault_Type=2为a相经10欧姆电阻接地短路\n(4)Fault_Type=3为bc两相相间短路\n(5)Fault_Type=4为bc两相短路接地\nFault_Type=')

a=-0.5+j*sqrt(3)/2;

T=[1 1 1 % T为对称分量法的变换矩阵，见P87公式（4-4）

a^2 a 1

a a^2 1];

%---------------第一步：计算短路点的序电流，相电流---------------%

%计根据故障类型选择不同的计算公式，计算故障点各序电流

if Fault_Type==0 %三相短路

I_Fault1=1/Z1(Fault_Node,Fault_Node);

I_Fault2=0;

I_Fault0=0;

elseif Fault_Type==1 %a相接地短路

I_Fault1=1/(Z1(Fault_Node,Fault_Node)+Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node));

6

课程 电力系统暂态上机课程设计报告 I_Fault2=I_Fault1;

I_Fault0=I_Fault1;

elseif Fault_Type==2 %a相经10欧姆电阻接地短路

Zf=10

I_Fault1=1/(Z1(Fault_Node,Fault_Node)+Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node)+3*Zf); I_Fault2=I_Fault1;

I_Fault0=I_Fault1;

elseif Fault_Type==3 %bc两相相间短路

I_Fault1=1/(Z1(Fault_Node,Fault_Node)+Z2(Fault_Node,Fault_Node));

I_Fault2=-1*I_Fault1;

I_Fault0=0;

elseif Fault_Type==4 %bc两相短路接地

I_Fault1=1/(Z1(Fault_Node,Fault_Node)+Z2(Fault_Node,Fault_Node)*Z0(Fault_Node,Fault_Node)/(Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node)));

I_Fault2=-1*I_Fault1*Z0(Fault_Node,Fault_Node)/(Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node));

I_Fault0=-1*I_Fault1*Z2(Fault_Node,Fault_Node)/(Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node));

end;

str='短路电流';

str

Iabc=T*[I_Fault1 I_Fault2 I_Fault0].' %相量

Iabc_effective=abs(Iabc) %有效值

7

课程 电力系统暂态上机课程设计报告

%---------------第二步：计算各个节点的序电压，相电压---------------% %故障分量

Ifault_node1=zeros(bus_Num1,1);

Ifault_node2=zeros(bus_Num1,1);

Ifault_node0=zeros(bus_Num1,1);

for m=1:bus_Num1

if m==Fault_Node

Ifault_node1(m,1)=-1*I_Fault1; %故障电流分量

Ifault_node2(m,1)=-1*I_Fault2;

Ifault_node0(m,1)=-1*I_Fault0;

else

Ifault_node1(m,1)=0;

Ifault_node2(m,1)=0;

Ifault_node0(m,1)=0;

end

end

Ufault_node1=zeros(bus_Num1,1);

Ufault_node2=zeros(bus_Num1,1);

Ufault_node0=zeros(bus_Num1,1);

Ufault_node1=Y1\Ifault_node1; %故障电压分量

Ufault_node2=Y2\Ifault_node2;

Ufault_node0=Y0\Ifault_node0;

%正常分量

Unormal_node1=ones(bus_Num1,1);

Unormal_node2=zeros(bus_Num1,1);

Unormal_node0=zeros(bus_Num1,1);

%各个节点的序电压

8

课程 电力系统暂态上机课程设计报告 Uall_node1=zeros(bus_Num1,1);

Uall_node2=zeros(bus_Num1,1);

Uall_node0=zeros(bus_Num1,1);

Uall_node1=Unormal_node1+Ufault_node1; %叠加原理

Uall_node2=Unormal_node2+Ufault_node2;

Uall_node0=Unormal_node0+Ufault_node0;

str='各节点序电压';

str

U120=[Uall_node1 Uall_node2 Uall_node0].' % 序分量

str='各节点相电压';

str

Uabc=T*[Uall_node1 Uall_node2 Uall_node0]' %相量('表示求转置矩阵）

Uabc_effective=abs(Uabc) %有效值

%---------------第三步：计算各个支路的序电流，相电流---------------%

Ibranch_1=zeros(branch_Num1,1);

Ibranch_2=zeros(branch_Num1,1);

Ibranch_0=zeros(branch_Num1,1);

for i=1:branch_Num1

Ibranch_1(i,1)=(Uall_node1(branch1_No1(i),1)-Uall_node1(branch1_No2(i),1))/z1_branch(i); %支路电压/支路阻抗

Ibranch_2(i,1)=(Uall_node2(branch1_No1(i),1)-Uall_node2(branch1_No2(i),1))/z2_branch(i);

Ibranch_0(i,1)=(Uall_node0(branch1_No1(i),1)-Uall_node0(branch1_No2(i),1))/z0_branch(i);

end

9

课程 电力系统暂态上机课程设计报告 str='各支路序电流';

str

Ibranch_120=[Ibranch_1 Ibranch_2 Ibranch_0].' % 序分量

str='各支路相电流';

str

Ibranch_abc=T*[Ibranch_1 Ibranch_2 Ibranch_0]' %相量('表示求转置矩阵） Ibranch_abc_effective=abs(Ibranch_abc) %有效值

%---------------第四步：计算发电机节点的序电压，相电压---------------% %故障分量

I2fault_node1=zeros(bus_Num2,1);

I2fault_node2=zeros(bus_Num2,1);

I2fault_node0=zeros(bus_Num2,1);

for m=1:bus_Num2

if m==Fault_Node

I2fault_node1(m,1)=-1*I_Fault1;

I2fault_node2(m,1)=-1*I_Fault2;

I2fault_node0(m,1)=-1*I_Fault0;

else

Ifault_node1(m,1)=0;

Ifault_node2(m,1)=0;

Ifault_node0(m,1)=0;

end

end

10

课程 电力系统暂态上机课程设计报告

U2fault_node1=zeros(bus_Num2,1);

U2fault_node2=zeros(bus_Num2,1);

U2fault_node0=zeros(bus_Num1,1);

U2fault_node1=YY1\I2fault_node1;

U2fault_node2=YY2\I2fault_node2;

%正常分量

U2normal_node1=ones(bus_Num2,1);

U2normal_node2=zeros(bus_Num2,1);

U2normal_node0=zeros(bus_Num2,1);

%各个节点的序电压

U2all_node1=zeros(bus_Num2,1);

U2all_node2=zeros(bus_Num2,1);

U2all_node0=zeros(bus_Num2,1);

U2all_node1=U2normal_node1+U2fault_node1;

U2all_node2=U2normal_node2+U2fault_node2;

U2120=[U2all_node1 U2all_node2 U2all_node0]'; % 序分量

U2abc=T*[U2all_node1 U2all_node2 U2all_node0]' ; %相量('表示求转置矩阵） U2abc_effective=abs(U2abc); %有效值

str='发电机节点序电压';

str

U4120=U2120(:,4)

U5120=U2120(:,5)

11

课程 电力系统暂态上机课程设计报告

str='发电机节点相电压相量及有效值';

str

U4abc=U2abc(:,4)

U5abc=U2abc(:,5)

U4abc=U2abc_effective(:,4)

U5abc=U2abc_effective(:,5)

%---------------第五步：计算发电机支路的序电流，相电流---------------%

Ibranch4_1=Ibranch_1(1,1)+Ibranch_1(2,1);

Ibranch4_2=Ibranch_2(1,1)+Ibranch_2(2,1);

Ibranch4_0=Ibranch_0(1,1)+Ibranch_0(2,1);

Ibranch5_1=-Ibranch_1(1,1)+Ibranch_1(3,1);

Ibranch5_2=-Ibranch_2(1,1)+Ibranch_2(3,1);

Ibranch5_0=-Ibranch_0(1,1)+Ibranch_0(3,1);

str='发电机支路序电流';

str

Ibranch4_120=[Ibranch4_1 Ibranch4_2 Ibranch4_0]' % 序分量

Ibranch5_120=2\[Ibranch5_1 Ibranch5_2 Ibranch5_0]' % 序分量

str='发电机支路相电流';

str

Ibranch4_abc=T*[Ibranch4_1 Ibranch4_2 Ibranch4_0]' %相量('表示求转置矩阵） Ibranch4_abc_effective=abs(Ibranch4_abc) %有效值

Ibranch5_abc=2\T*[Ibranch5_1 Ibranch5_2 Ibranch5_0]' %相量('表示求转置矩阵）

12

课程 电力系统暂态上机课程设计报告 Ibranch5_abc_effective=abs(Ibranch5_abc) %有效值

静态稳定计算stability_smallsignal

clear;

clc;

U=1; %系统电压

Eq0=1.972; %空载电动势

UG0=1.21; %机端电压

Xe=0.504; %线路电抗值

Xd=0.982; %同步电抗

Xdd=0.344; %暂态电抗

Xd_all=1.486; %系统电抗(包含同步电抗)

Xdd_all=0.848; %系统电抗(包含暂态电抗)

Tj=10; %惯性时间常数

Td0=10; %励磁绕组时间常数

Ke=0.1;

while(Ke<5.7)

for delta=49:110

%请同学们求各个参数

Eq1=Xdd_all*Eq0/Xd_all+U*(1-Xdd_all/Xd_all)*cos(delta*pi/180); UGd0=U*sin(delta*pi/180)*Xd/Xd_all;

UGq0=Eq0*Xe/Xd_all+Xd*U*cos(delta*pi/180)/Xd_all;

UG=sqrt(UGd0^2+UGq0^2);

a=Ke^2*Xe^2/(Xd_all^2)-1;

b=2*(Ke^2*Xd*Xe*U*cos(delta*pi/180)/(Xd_all^2)+(Eq0+Ke*UG0));

13

课程 电力系统暂态上机课程设计报告

c=Ke^2*Xd^2*U^2*cos(delta)*cos(delta)/(Xd_all^2)+Ke^2*Xd^2*U^2*sin(delta)*sin(delta)/(Xd_all^2)-(Eq0+Ke*UG0)^2;

Eq=(-b+sqrt(b*b-4*a*c))/(2*a);%求解方程

K1=Eq1*U/Xdd_all*cos(delta*pi/180)+U*U*(Xdd_all-Xd_all)/(Xdd_all*Xd_all)*cos(2*delta*pi/180); K2=U*sin(delta*pi/180)/Xdd_all;

K3=Xdd_all/Xd_all;

K4=U*sin(delta*pi/180)*(Xd_all-Xdd_all)/Xdd_all;

K5=UGd0*U*Xd*cos(delta*pi/180)/(UG*Xd_all)-UGq0*U*Xdd*sin(delta*pi/180)/(UG*Xdd_all); K6=UGq0*(Xd_all-Xd)/(UG*Xdd_all);

Kemin=-(K1-K2*K3*K4)/(K3*(K1*K6-K2*K5));

Kemax=-K4/K5;

A=zeros(3);

%请同学们求P186公式7-14状态方程各个参数

A(1,2)=2*pi*50;

A(2,1)=-K1/Tj;

A(2,3)=-K2/Tj;

A(3,1)=-(K4+Ke*K5)/Td0;

A(3,3)=-(1/K3+Ke*K6)/Td0;

V=eig(A);%求特征值

V1=real(V(1));%取出特征根的实部

V2=real(V(2));

V3=real(V(3));

if (V1<0&V2<0&V3<0)%判断稳定性

Ke_result=Ke;

delta_result=delta;

P_result=Eq*U*sin(delta*pi/180)/Xd_all;

14

课程 电力系统暂态上机课程设计报告 else

s=sprintf('Ke=%f,delta=%d,P=%f,Eq=%f\n',Ke_result,delta_result,P_result,Eq); disp(s);

break;

end;

end;

Ke=Ke+0.1;

end;

暂态稳定计算stability_transient

clear;

clc;

f=50; %系统额定频率

Tj=8.47; %归算后的发电机惯性时间常数

PT=1; %正常运行时发电机向无穷大系统传输的有功功率

P2M=0.48; %故障存在时发电机的最大功率

P3M=1.38; %故障切除后发电机的最大功率

%下面是利用改进欧拉法进行逐段计算

%需要注意故障切除前后电磁功率有跃变

h=0.001; %设置步长0.001s

Duration=2; %设置计算时段长度2s

CutTime=input('输入故障切除时间\n');

%故障发生时的功角变化过程

Delta(1)=33.92; %初始功角

Omega(1)=1; %初始转速

t(1)=0;

for i=1:round(CutTime/h)

dD(1)=(Omega(i)-1)*360*f;%时段始的变化率

dO(1)=(PT-P2M*sin(Delta(i)*pi/180))/Tj;

15

课程 电力系统暂态上机课程设计报告

Delta(i+1)=Delta(i)+dD(1)*h;%时段莫的估计值

Omega(i+1)=Omega(i)+dO(1)*h;

dD(2)=(Omega(i+1)-1)*360*f;%时段末对应于估计值的变化率 dO(2)=(PT-P2M*sin(Delta(i+1)*pi/180))/Tj;

dD(2)=(dD(1)+dD(2))/2;%时段中的平均变化率

dO(2)=(dO(1)+dO(2))/2;

Delta(i+1)=Delta(i)+dD(2)*h;%时段末的值

Omega(i+1)=Omega(i)+dO(2)*h;

t(i+1)=i*h;

end;

%故障切除后的功角变化过程

for i=round(CutTime/h)+1:round(Duration/h)

dD(1)=(Omega(i)-1)*360*f;%时段始的变化率

dO(1)=(PT-P3M*sin(Delta(i)*pi/180))/Tj;

Delta(i+1)=Delta(i)+dD(1)*h;%时段莫的估计值

Omega(i+1)=Omega(i)+dO(1)*h;

dD(2)=(Omega(i+1)-1)*360*f;%时段末对应于估计值的变化率 dO(2)=(PT-P3M*sin(Delta(i+1)*pi/180))/Tj;

dD(2)=(dD(1)+dD(2))/2;%时段中的平均变化率

dO(2)=(dO(1)+dO(2))/2;

Delta(i+1)=Delta(i)+dD(2)*h;%时段末的值

Omega(i+1)=Omega(i)+dO(2)*h;

16

课程 电力系统暂态上机课程设计报告

t(i+1)=i*h;

end;

plot(t,Delta);

17