电力系统下课程设计报告-短路电流计算

《电力系统分析》课程设计报告

题目： 3G9bus短路电流计算 

系    别     电气工程学院     

专业班级     10级电气四班     

学生姓名                

学    号               

指导教师                

提交日期    2012年12月10日  

一、设计目的

1.掌握电力系统短路计算的基本原理；

2.掌握并能熟练运用一门计算机语言（MATLAB语言或FORTRAN或C语言或C++语言）；

3.采用计算机语言对短路计算进行计算机编程计算。

二、短路电流计算的基本原理和方法

2.1电力系统节点方程的建立

利用节点方程作故障计算，需要形成系统的节点导纳（或阻抗）矩阵。一般短路电流计算以前要作电力系统的潮流计算，假定潮流计算的节点导纳矩阵已经形成，在此基础上通过追加支路的方式形成电力短路电流计算的节点导纳矩阵YN。

1）对发电机节点

在每一发电机节点增加接地有源支路 与串联

求短路稳态解:  

求短路起始次暂态电流解:  

一般情况下发电机定子绕组电阻忽略掉，并将与的有源支路转化成电流源与导纳并联的形式

2）负荷节点的处理

负荷节点在短路计一算中一般作为节点的接地支路，并用恒定阻抗表示，其数值由短路前瞬间的负荷功率和节点实际电压算出，即首先根据给定的电力系统运行方式制订系统的等值电路，并进行各元件标么值参数的计算，然后利用变压器和线路的参数形成不含发电机和负荷的节点导纳矩阵 YN 。

3）如果短路故障发生在输电线路内，在短路点要增加新节点，将一条输电线分成两段，并形成短路电流计算的节点导纳矩阵Y，最后形成包括所有发电机支路和负荷支路的节点方程如下：

一般，节点导纳矩阵 Y 阵与 YN 阶次相同。其差别只在于 YN 阵不含发电机和负荷的阻抗；节点注入电流向量 I中只有发电机端节点的电流不为零。有非零电流源注人的节点称为有源节点。系统中的同步调相机可按发电机处理。任进行起始次暂态电流计算时，大型同步电动机、感应电动机以及以电动机为主要成分的综合负荷，特别是在短路点近处的这些负荷，必要时也可以用有源支路表示，并仿照发电机进行处理。必须指出，在计算机已普遍应用的情况下，如果有必要的话，只要能提供短路计算所需的准确的原始数据，对短路进行更精确的计算并不存在什么障碍。

2.2利用节点阻抗矩阵计算短路电流

假定系统中的节点f经过渡阻抗zf发生短路。过渡阻抗渡阻抗zf，不参与形成网络的节点导纳（或阻抗）矩阵。图6-3中方框内的有源网络代表系统正常状态的单相等值网络。现在我们保持故障处的边界条件不变，把网络的原有部分同故障支路分开（见图6-3）。容易看出，对于正常状态的网络而言，发生短路相当于在故障节点f增加了一个注人电流一 If （短路电流以流出故障点为正，节点电流则以注入为正）。因此，网络中任一节点 i 的电压可表示为：    （6-3）

由式（6-3）可见，任一节点 i 的电压郁由两项叠加而成。第一项表示当 If=0时由网络内所有电源在节点i产生的电压，也就是短路前瞬间正常运行状态一 F 的节点电压，这是节点电压的正常分量，记为是电网的潮流解。第二项是当网络中所有电流源都断开，仅仅由短路电流If在节点i产生的电压，这就是节点电压的故障分童。上述两个分量的叠加，就等于发生短路后节点，的实际电压，即

     （6-4）

式中，是短路前故障点的正常电压；当i=f时，是故障节点 f 的自阻抗，也称输入阻抗。（边界条件方程）带入可得：

对于非变压器支路，令 k=1。从计算公式（6-7）和（6-8）可以看到，式中所用到的阻抗矩阵元素都带有列标f。这就是说，如果网络在正常状态下的节点电压为已知，为了进行短路计算，只须利用节点阻抗矩阵中与故障点f对应的一列元素。因此，尽管是采用了阻抗型的节点方程，但是并不需要作出全部阻抗矩阵。在短路的实际计一算中，一般只需形成网络的节点导纳矩阵．并根据具体要求，用第四章所讲的方法求出阻抗矩阵的某一列或某几列元素即可。在应用节点阻抗矩阵进行短路计算时，我们都将采用这种算法。

1）解潮流计算，

2）修正解潮流的YN形成Y

3）指定短路点f

4）计算节点阻抗矩阵第k列

5）计算短路电流：

6）计算节点电压： 

7）计算支路电流：

8）输出计算结果

三、3G9bus短路电流在计算机的编程

3.1、三机九节点系统

图1  三机九节点系统

表1  九节点系统支路参数

3.2 程序设计

主函数

Sbase_MVA=100.

fid=fopen('Nodedata.txt');

N=textscan(fid, '%s %u %d %f %f %f %f %f %f')

fclose(fid);

busnumber=size(N{1},1)

for i=1:busnumber

    Bus(i).name=N{1}(i);

    Bus(i).type=N{2}(i);

    Bus(i).no=i;

    Bus(i).Base_KV=N{3}(i);

    Bus(i).PG=N{4}(i);

    Bus(i).QG=N{5}(i);

    Bus(i).PL=N{6}(i);

    Bus(i).QL=N{7}(i);

    Bus(i).pb=N{8}(i);

    Bus(i).V=1.0;

    Bus(i).angle=0;

end

fid=fopen('Aclinedata.txt');

A=textscan(fid, '%s %s %f %f %f %f')

fclose(fid);

aclinenumber=size(A{1},1)

for i=1:aclinenumber

    Acline(i).fbname=A{1}(i);

    Acline(i).tbname=A{2}(i);

    Acline(i).Base_KV=A{3}(i);

    Acline(i).R=A{4}(i);

    Acline(i).X=A{5}(i);

    Acline(i).hB=A{6}(i);

    for k=1:busnumber

       if strcmp(Acline(i).fbname, Bus(k).name)

       Acline(i).fbno=Bus(k).no;

       end

       if strcmp(Acline(i).tbname, Bus(k).name)

       Acline(i).tbno=Bus(k).no;

       end

    end    

end

fid=fopen('Transdata.txt');

T=textscan(fid, '%s %f %f %s %f %f %f %f')

fclose(fid);

tansnumber=size(T{1},1)

for i=1:tansnumber

    Trans(i).fbname=T{1}(i);

    Trans(i).fbBase_KV=T{2}(i);

    Trans(i).fbrated_KV=T{3}(i);

    Trans(i).tbname=T{4}(i);

    Trans(i).tbBase_KV=T{5}(i);

    Trans(i).tbrated_KV=T{6}(i);

    Trans(i).R=T{7}(i);

    Trans(i).X=T{8}(i);

    for k=1:busnumber

       if strcmp(Trans(i).fbname, Bus(k).name)

       Trans(i).fbno=Bus(k).no;

       end

       if strcmp(Trans(i).tbname, Bus(k).name)

       Trans(i).tbno=Bus(k).no;

       end

    end

    Trans(i).k=Trans(i).tbrated_KV*Trans(i).fbBase_KV/Trans(i).fbrated_KV/Trans(i).tbBase_KV;

    tempx=Trans(i).fbrated_KV^2/Trans(i).fbBase_KV^2;

    Trans(i).X=tempx*Trans(i).X;

    Trans(i).R=tempx*Trans(i).R;

end

%N=0

%Trans(1)

%Trans(2)

% for Y=G+jB matrix

[G,B,B2]=FormYmatrix(Bus,busnumber,Acline,aclinenumber,Trans,tansnumber);   

%B:=B';B2:=B"

dlmwrite('Gmatrix.txt', G, 'delimiter', '\','precision', 6);

dlmwrite('Bmatrix.txt', B, 'delimiter', '\','precision', 6);

G

B

B2

pause

[JP,JQ]=FormJPQmatrix(Bus,B,B2,busnumber);

JP

iJP=-inv(JP)

JQ

iJQ=-inv(JQ)

pause

%maxiteration=0

for i=1:busnumber

   NodeV(i)=Bus(i).V;

   Nodea(i)=Bus(i).angle;

   VX(i)=Bus(i).V*cos(Bus(i).angle);

   VY(i)=Bus(i).V*sin(Bus(i).angle);

   dQGQL(i)=Bus(i).QG-Bus(i).QL;

   dPGPL(i)=Bus(i).PG-Bus(i).PL;      

end

NodeV=NodeV'

Nodea=Nodea'

%VX=VX'

%VY=VY'

dQGQL=dQGQL'

dPGPL=dPGPL'

pause

%for nointer=1:10

maxdP=1.;

maxdQ=1.;

epsilon=0.000001;

noiteration=0;

while (maxdP>epsilon)&(maxdP>epsilon)

[deltaP,deltaQ,maxdP,maxdQ]=FormdPQvector(Bus,NodeV,Nodea,dQGQL,dPGPL,B,G,busnumber);

deltaP;

deltaQ;

maxdP;

maxdQ;

da=iJP*deltaP;

dV=iJQ*deltaQ;

Nodea=Nodea+da;

NodeV=NodeV+dV;

noiteration=noiteration+1;

if noiteration>20

    break

    end

end

for i=1:busnumber

   Bus(i).V=NodeV(i);

   NodeV(i)=NodeV(i)*Bus(i).Base_KV;

   Bus(i).angle=Nodea(i);

   Nodea(i)=Nodea(i)*180/pi;

end

noiteration

Nodea=Nodea'

NodeV=NodeV'

Clear

子函数

%生成G、B矩阵

function [G,B,X]=FormYmatrix(Bus,busnumber,Acline,aclinenumber,Trans,tansnumber)

Y=zeros(busnumber);

X=zeros(busnumber);

for i=1:busnumber

Y(i,i)=Y(i,i)+Bus(i).pb*j;

end

for i=1:aclinenumber

    f=Acline(i).fbno;

    t=Acline(i).tbno;

    Y(f,f)=Y(f,f)+Acline(i).hB*j+1/(Acline(i).R+Acline(i).X*j);

    Y(t,t)=Y(t,t)+Acline(i).hB*j+1/(Acline(i).R+Acline(i).X*j);

    Y(f,t)=Y(f,t)-1/(Acline(i).R+Acline(i).X*j);

    Y(t,f)=Y(t,f)-1/(Acline(i).R+Acline(i).X*j);

    X(f,f)=X(f,f)-1/Acline(i).X;

    X(t,t)=X(t,t)-1/Acline(i).X;

    X(f,t)=1/Acline(i).X;

    X(t,f)=1/Acline(i).X;    

end

for i=1:tansnumber

    f=Trans(i).fbno;

    t=Trans(i).tbno;

    Y(f,f)=Y(f,f)+1/(Trans(i).R+Trans(i).X*j);

    Y(t,t)=Y(t,t)+1/(Trans(i).R+Trans(i).X*j)/Trans(i).k^2;

    Y(f,t)=Y(f,t)-1/(Trans(i).R+Trans(i).X*j)/Trans(i).k;

    Y(t,f)=Y(t,f)-1/(Trans(i).R+Trans(i).X*j)/Trans(i).k;

    X(f,f)=X(f,f)-1/Trans(i).X;

    X(t,t)=X(t,t)-1/Trans(i).X;

    X(f,t)=1/Trans(i).X;

    X(t,f)=1/Trans(i).X;  

end

G=real(Y);

B=imag(Y);

end

%生成JP、JQ矩阵

function [JP,JQ]=FormJPQmatrix(Bus,B,B2,busnumber)

JP=B;

JQ=B2;

for i=1:busnumber

    if Bus(i).type==1

        for k=1:busnumber

        JQ(i,k)=0.;

        JQ(k,i)=0.;

        JP(i,k)=0.;

        JP(k,i)=0.;            

        end

        JQ(i,i)=1.;

        JP(i,i)=1.;

    end

     if Bus(i).type==3

        for k=1:busnumber

        JQ(i,k)=0.;

        JQ(k,i)=0.;     

        end

        JQ(i,i)=1.;

    end 

end

end

%计算偏节点PQ差量

function [deltaP,deltaQ,maxdP,maxdQ]=FormdPQvector(Bus,NodeV,Nodea,dQGQL,dPGPL,B,G,busnumber)

       deltaQ=dQGQL;

       deltaP=dPGPL;

       maxdP=0.;

       maxdQ=0.;    

for i=1:busnumber

       if Bus(i).type==1

       deltaQ(i)=0.;

       deltaP(i)=0.;

       end

       if Bus(i).type==3

        deltaQ(i)=0.;

        %y1=0;

        %y2=0;

        y3=0;

        for k=1:busnumber

            if (B(i,k)~=0|G(i,k)~=0)

            %y1=y1+(G(i,k)*VX(k)-B(i,k)*VY(k));

            %y2=y2+(G(i,k)*VY(k)+B(i,k)*VX(k));

            y3=y3+NodeV(k)*(G(i,k)*cos(Nodea(i)-Nodea(k))+B(i,k)*sin(Nodea(i)-Nodea(k)));

            end

        end

        deltaP(i)=deltaP(i)-y3*NodeV(i);

        %deltaP2(i)=(deltaP2(i)-(y1*VX(i)+y2*VY(i)))/Bus(i).V;

       end

       if Bus(i).type==2

        %y1=0;

        %y2=0;

        y3=0;

        y4=0;

        for k=1:busnumber

            if (B(i,k)~=0|G(i,k)~=0)

            %y1=y1+(G(i,k)*VX(k)-B(i,k)*VY(k));

            %y2=y2+(G(i,k)*VY(k)+B(i,k)*VX(k));

            y3=y3+NodeV(k)*(G(i,k)*cos(Nodea(i)-Nodea(k))+B(i,k)*sin(Nodea(i)-Nodea(k)));

            y4=y4+NodeV(k)*(G(i,k)*sin(Nodea(i)-Nodea(k))-B(i,k)*cos(Nodea(i)-Nodea(k)));

            end

        end

        deltaP(i)=deltaP(i)-y3*NodeV(i);

        %deltaP2(i)=(deltaP2(i)-(y1*VX(i)+y2*VY(i)))/Bus(i).V;

        deltaQ(i)=deltaQ(i)-y4*NodeV(i);

        %deltaQ2(i)=(deltaQ2(i)-(y1*VY(i)-y2*VX(i)))/Bus(i).V;

       end

if maxdP           maxdP=abs(deltaP(i)); 

       end

if maxdQ           maxdQ=abs(deltaQ(i)); 

       end

     deltaP(i)=deltaP(i)/NodeV(i);

     deltaQ(i)=deltaQ(i)/NodeV(i);

end

end

3.3输出并计算结果

由以上程序可以得出：

（1）进行系统正常运行状态的潮流计算，求得：

己知公式＝1-（Zif/(Zff+Zf)）因为金属性短路时Zf＝0公式为

＝1-（Zif/Zff）再根据公式可得，

V1&(0) =0.2774 ； V2&(0) ＝0.6770； V3&(0)＝0.6770； V4&(0) ＝0；

V5&(0)＝0.1870； V6&(0)＝0.1934；  V7&(0) ＝0.54197；V8&(0)＝0.5456；

V9&(0)＝0.5511

（2）不含发电机和负荷的节点导纳矩阵YN：

(3)形成包括所有发电机支路和负荷支路的节点导纳矩阵Y，即在YN中的发电机节点和负荷节点的自导纳上分别增加发电机导纳和负荷导纳：

（4），计算节点阻抗矩阵，从而得到阻抗矩阵中的第f列；

（5）计算短路电流

因为Zf＝0，所以短路电流公式为：If&(0)=1/ Zff

母线对地短路电流标幺值：If = - 9.1239i

母线对地短路电流有名值：Ifymzh＝2.2903

（6）计算各支路的短路电流

己知公式 Ipq&(0)=(k*Vp&(0)- Vq&(0))/Zpq 所以

I14＝4.85i；  I27＝2.16i； I39＝2.15i；

I45＝-0.7399+6.207i；  I46＝-0.37562+2.03277i；

I57＝-0.4215+2.12099i； I69＝-0.458575+1.998i；

I78＝-0.00587+0.04972i； I＝-0.00635+0.05381i；

四．总结

MATLAB是一个数学建模和计算软件，通过编程能够将复杂的矩阵计算转移到计算机上进行演算，极大的降低了人手计算的强度，给计算工作带来了极大的便利。但是对于初学的我来说也存在着一定的难度，不过通过网上丰富的资源和同学的帮助我慢慢的掌握了这个软件的基本操作，对我以后工作打了一定的基础。它不仅有强大的元素按功能，他还有强大的绘图功能，我对它的了解也仅仅就是一点点，或许说还没入门。

通过这次的课程设计，锻炼到电力系统潮流分析和短路电流计算的能力。3机9节点的计算，通过计算机和软件演算的辅助，可以大大节省时间，计算的结果也非常精确。通过短路电流计算，可以判断回路电气设备能不能承受这个数值的短路电流，也分析出电气设备中的短路保护装置能不能分断这个数值的电流。