Eacfault.m
上传用户:eighthdate
上传日期:2014-05-24
资源大小:270k
文件大小:3k
源码类别:

其他行业

开发平台:

Matlab

  1. function eacfault(Pm, E, V, X1, X2, X3)
  2. % This program obtains the power angle curves for a one-machine system
  3. % before fault, during fault and after the fault clearance.
  4. % The equal area criterion is applied to find the critical clearing angle
  5. % for the machine to stay synchronized to the infinite bus bar
  6. %
  7. % Copyright (c) 1998 by H. Saadat
  8. if exist('Pm')~=1
  9. Pm = input('Generator output power in p.u. Pm = '); else, end
  10. if exist('E')~=1
  11. E = input('Generator e.m.f. in p.u. E = '); else, end
  12. if exist('V')~=1
  13. V = input('Infinite bus-bar voltage in p.u. V = '); else, end
  14. if exist('X1')~=1
  15. X1 = input('Reactance before Fault in p.u. X1 = '); else, end
  16. if exist('X2')~=1
  17. X2 = input('Reactance during Fault in p.u. X2 = '); else, end
  18. if exist('X3')~=1
  19. X3 = input('Reactance aftere Fault in p.u. X3 = '); else, end
  20. Pe1max = E*V/X1; Pe2max=E*V/X2; Pe3max=E*V/X3;
  21. delta = 0:.01:pi;
  22. Pe1 = Pe1max*sin(delta); Pe2 = Pe2max*sin(delta); Pe3 = Pe3max*sin(delta);
  23. d0 =asin(Pm/Pe1max); dmax = pi-asin(Pm/Pe3max);
  24. cosdc = (Pm*(dmax-d0)+Pe3max*cos(dmax)-Pe2max*cos(d0))/(Pe3max-Pe2max);
  25.   if abs(cosdc) > 1
  26.   fprintf('No critical clearing angle could be found.n')
  27.   fprintf('system can remain stable during this disturbance.nn')
  28.   return
  29.   else, end
  30. dc=acos(cosdc);
  31. if dc > dmax
  32. fprintf('No critical clearing angle could be found.n')
  33.   fprintf('System can remain stable during this disturbance.nn')
  34.   return
  35.   else, end
  36. Pmx=[0  pi-d0]*180/pi; Pmy=[Pm Pm];
  37. x0=[d0 d0]*180/pi; y0=[0 Pm]; xc=[dc dc]*180/pi; yc=[0 Pe3max*sin(dc)];
  38. xm=[dmax dmax]*180/pi; ym=[0 Pe3max*sin(dmax)];
  39. d0=d0*180/pi; dmax=dmax*180/pi; dc=dc*180/pi;
  40. x=(d0:.1:dc);
  41. y=Pe2max*sin(x*pi/180);
  42. y1=Pe2max*sin(d0*pi/180);
  43. y2=Pe2max*sin(dc*pi/180);
  44. x=[d0 x dc];
  45. y=[Pm y Pm];
  46. xx=dc:.1:dmax;
  47. h=Pe3max*sin(xx*pi/180);
  48. xx=[dc xx  dmax];
  49. hh=[Pm h  Pm];
  50. delta=delta*180/pi;
  51. if X2 == inf
  52. fprintf('nFor this case tc can be found from analytical formula. n')
  53. H=input('To find tc enter Inertia Constant H, (or 0 to skip) H = ');
  54.    if H ~= 0
  55.    d0r=d0*pi/180; dcr=dc*pi/180;
  56.    tc = sqrt(2*H*(dcr-d0r)/(pi*60*Pm));
  57.    else, end
  58. else, end
  59. %clc
  60. fprintf('nInitial power angle     = %7.3f n', d0)
  61. fprintf('Maximum angle swing     = %7.3f n', dmax)
  62. fprintf('Critical clearing angle = %7.3f nn', dc)
  63. if X2==inf & H~=0
  64. fprintf('Critical clearing time  = %7.3f sec. nn', tc)
  65. else, end
  66. h = figure; figure(h);
  67. fill(x,y,'m')
  68. hold;
  69. fill(xx,hh,'c')
  70. plot(delta, Pe1,'-', delta, Pe2,'r-',  delta, Pe3,'g-', Pmx, Pmy,'b-', x0,y0, xc,yc, xm,ym), grid
  71. Title('Application of equal area criterion to a critically cleared system')
  72. xlabel('Power angle, degree'), ylabel(' Power, per unit')
  73. text(5, 1.07*Pm, 'Pm')
  74. text(50, 1.05*Pe1max,['Critical clearing angle = ',num2str(dc)])
  75. axis([0 180  0 1.1*Pe1max])
  76. hold off;