October 21, 2010¡cticas... · CHAPTER 1. CRITERIO DE LAS ÁREAS. 2 1.1.1 Cálculos eléctricos....

Estabilidad transitoria de gran perturbación. Criterio de las áreas

October 21, 2010

Chapter 1

Criterio de las áreas.

Permite analizar la estabilidad de un generador conectado a un nudo de potencia in�nita cuando el factor de amortiguación esnulo.

En el anexo teneís las fórmulas correspondientes a la resolución de este problema.

Se utilizará el programa matlab.

Se creará un �chero matlab (�chero.m), para ir introduciendo las operaciones.

% signi�ca un comentario, el cual no hay falta introducir en el �chero.

1.1 Caso 1.

Consideramos un generador conectado a un nudo de potencia in�nita a través de un transformador elevador y una línea:

Vamos a considerar que el generador, el transformador y la línea se representan por impedancias inductivas:

Tendremos como datos:

� Generador: potencia nominal Sg (MVA), tensión nominal Ug, H (s), reactancia síncrona (por unidad) referida a la potenciay tensión del generador.

� Transformador: tensiones entrada-salida; reactancia transformador Xt (por unidad).

� Líneas: reactancia línea Xl (por unidad)

� Nudo potencia in�nita: U∞

Las operaciones se realizarán en el entorno Matlab. Se introduciran los datos en un �chero, para tener posibilidad de cambiarcosas.

1

CHAPTER 1. CRITERIO DE LAS ÁREAS. 2

1.1.1 Cálculos eléctricos.

� Vamos a calcular los datos del circuito eléctrico (tensiones, ángulos, ...)

Se parte como dato del factor de potencia en bornes del generador como dato cosf y del valor de intensidad, que suponemos quees a plena carga (1 por unidad)

% Datos

H=3; xg=0.3; xt=0.15; xl=0.1;cosf=0.8;

%desfase intensidad. La tensión en bornes del generador se toma a 0º grados.

angulo=acos(cosf); angulogrados=angulo*180/pi;

%intensidad compleja

intensidad=1*(cos(-angulo)+j*sin(-angulo));

%cálculo fuerza electromotriz e = u+ jxg.i

e=1+j*xg*intensidad; emodulo=abs(e); eangulo=angle(e); eangulogrados=eangulo*180/pi;

%cálculo tensión barra potencia in�nita e = u∞ + jxtotal.i

xtotal=xg+xt+xl;

uin�nita=e-1j*xtotal*intensidad;

uin�nitamodulo=abs(uin�nita);

uin�nitaangulo=angle(uin�nita); uin�nitaangulogrados=uin�nitaangulo*180/pi;

%ángulo de potencia. Desfase entre e, u∞

delta0=eangulo-uin�nitaangulo; delta0grados=delta0*180/pi;

1.1.2 Representación de las curvas.

%función potencia eléctrica (se representa por medio de vectores)

pmax=(emodulo*uin�nitamodulo)/xtotal; % valor máximo de la senoide

delta=linspace(0,pi,1000); % puntos a representar en la grá�ca (eje de abcisas)

pe=pmax*sin(delta); plot(delta,pe); hold on %función potencia eléctrica

%función potencia mecánica (se representa por medio de vectores y es una línea recta)

pmo=pmax*sin(delta0);

pm=pmo+3*delta-3*delta; plot(delta,pm);

1.1.3 Cálculo del ángulo crítico.

Se van a calcular áreas mediante aproximaciones (función trapz en matlab). Deberemos saber la situación de cada punto en losvectores que representan las curvas. Para ello utilizaremos una función llamada posicion_vector (anexo) para tal �n.

%dimension vectores pe (potencia eléctrica) y pm (potencia mecánica)

dimension=length(pe);


%posición que corresponde a δo (delta0) en el vector pe (potencia eléctrica) cuando la potencia mecánica es pmo

posdelta0=posicion_vector(dimension,pe,pmo);

%cálculo ángulo crítico mediante formula exacta

deltacriticaexacta=acos(sin(delta0)*(pi-2*delta0)-cos(delta0));

% cálculo mediante aproximaciones de las integrales. En el anexo está la función calculodeltacritica.

deltacritica=calculodeltacritica(pmax,pmo);

1.1.4 Dibujar el área de aceleración y de desaceleración

% valor de potencia eléctrica para el ángulo crítico δcritica

pedeltacritica=pmax*sin(deltacritica);

% posición que corresponde en el vector de potencia eléctrica al valor del ángulo crítico δcritica

posdeltacritica=posicion_vector(dimension,pe,pedeltacritica);

if deltacritica>(pi/2) %comprobar si el ángulo crítico es superior a 90º

posdeltacritica=dimension-posdeltacritica;

end

%dibujar área aceleradora

cero=0.0001+3*delta-3*delta; % función 0 (eje de abcisas)

puntos=[cero(posdelta0:posdeltacritica),�iplr(pm(posdelta0:posdeltacritica))];

t1=[delta(posdelta0:posdeltacritica),�iplr(delta(posdelta0:posdeltacritica))];

f1=�ll(t1,puntos,1);

%dibujar área desaceleradora

posdeltapi0=dimension-posdelta0;

puntos_2=[pm(posdeltacritica:posdeltapi0),�iplr(pe(posdeltacritica:posdeltapi0))];

t1_2=[delta(posdeltacritica:posdeltapi0),�iplr(delta(posdeltacritica:posdeltapi0))];

f1=�ll(t1_2,puntos_2,3);

1.1.5 Cálculo del tiempo crítico

%tiempo duración fallo alcanza deltacritica (fórmula exacta)

tdeltacritica=sqrt((4*H)/(2*pi*50*pmo)*(deltacritica-delta0));

1.2 Caso 2.

Sea un generador conectado a un nudo de potencia in�nita a través de un transformador y dos líneas en parelelo. Se va a estudiarla estabilidad cuando se produce un fallo al comienzo de una de las líneas y la falta se despeja por apertura de la línea afectada.



Vamos a calcular los datos del circuito eléctrico (tensiones, ángulos, ...)


%Datos

H=3; xg=0.3; xt=0.15; xl1=0.2; xl2=0.2; cosf=0.8;

%paralelo líneas

xl=(xl1*xl2)/(xl1+xl2);


intensidad angulo=acos(cosf); angulogrados=angulo*180/pi;



%cálculo fuerza electromotriz. e = ug + jxg.i

e=1+j*xg*intensidad;

emodulo=abs(e);

eangulo=angle(e); eangulogrados=eangulo*180/pi;

%cálculo tensión barra potencia in�nita. e = u∞ + jxtotal.i

xtotal=xg+xt+xl;




%ángulo de potencia. Desfase e, u∞


1.2.2 Representación de las curvas.

%función potencia eléctrica con las dos líneas funcionando (se representa por medio de vectores)



pe=pmax*sin(delta); plot(delta,pe); hold on

%potencia mecánica (se representa por medio de vectores)

pmo=pmax*sin(delta0); pm=pmo+3*delta-3*delta; plot(delta,pm);

%dimension vectores



%posición vector pe a deltao

posdelta0=posicion_vector(dimension,pe,pmo);

%curva potencia eléctrica con una sola línea funcionando

xtotal2=xg+xt+xl1;

pmax2=(emodulo*uin�nitamodulo)/xtotal2; pe2=pmax2*sin(delta); plot(delta,pe2);

1.2.3 Cálculo del ángulo crítico

Se van a calcular áreas mediante aproximaciones (función trapz en matlab). Deberemos saber las situación de cada punto en losvectores que representan las curvas. Para ello utilizaremos una función llamada posicion_vector para tal �n.

% posición vector potencia electrica funcionando solo una línea para el ángulo inicial δo

delta02=asin(pmo/pmax2);

%cálculo ángulo crítico

deltacritica2=calculodeltacritica(pmax2,pmo);

%potencia eléctrica cuando tenemos una línea en el punto δcritica

pedeltacritica2=pmax2*sin(deltacritica2);

%posición en el vector del ángulo crítico

posdeltacritica2=posicion_vector(dimension,pe2,pedeltacritica2);

if deltacritica2>(pi/2)

posdeltacritica2=dimension-posdeltacritica2;

end



cero=0.0001+3*delta-3*delta;

puntos=[cero(posdelta0:posdeltacritica2),�iplr(pm(posdelta0:posdeltacritica2))];

t1=[delta(posdelta0:posdeltacritica2),�iplr(delta(posdelta0:posdeltacritica2))];

f1=�ll(t1,puntos,1)


posdelta02=posicion_vector(dimension,pe2,pmo);

posdeltapi0=dimension-posdelta02;

puntos_2=[pm(posdeltacritica2:posdeltapi0),�iplr(pe2(posdeltacritica2:posdeltapi0))];

t1_2=[delta(posdeltacritica2:posdeltapi0),�iplr(delta(posdeltacritica2:posdeltapi0))];



1.2.5 Cálculo del tiempo crítico

%tiempo duración fallo alcanza deltacritica

tdeltacritica=sqrt((4*H)/(2*pi*50*pmo)*(deltacritica2-delta02));

1.3 Caso 3.

Sea un generador conectado a un nudo de potencia in�nita a través de un transformador y dos líneas en parelelo. Se va a estudiarla estabilidad cuando se produce un fallo en un punto intermedio de una de las líneas y la falta se despeja por apertura de lalínea afectada.


Vamos a calcular los datos del circuito eléctrico (tensiones, ángulos, ...)


% Datos

H=3; xg=0.3; xt=0.15; xl1=0.2; xl2=0.2; xl2mitad=0.1; cosf=0.8;


angulo=acos(cosf); angulogrados=angulo*180/pi;



%cálculo fuerza electromotriz e = ug + jxg.i

e=1+j*xg*intensidad; emodulo=abs(e); eangulo=angle(e); eangulogrados=eangulo*180/pi;

%cálculo tensión barra potencia in�nita. e = u∞ + jxtotal.i

xtotal=xg+xt+xl;




%ángulo de potencia. Desfase e, u∞


1.3.2 Representación de las curvas

%función potencia eléctrica




pe=pmax*sin(delta); plot(delta,pe); hold on %función potencia eléctrica

%potencia mecánica

pmo=pmax*sin(delta0); pm=pmo+3*delta-3*delta; plot(delta,pm);

%dimension vectores pe (potencia eléctrica) y pm (potencia mecánica)


%posición que corresponde a δo (delta0) en el vector pe (potencia eléctrica)

posdelta0=posicion_vector(dimension,pe,pmo); % posicion_vector en una función (ver anexo)

%curva potencia eléctrica durante el fallo (en medio de una de las líneas)

porcentajefallo=0.5;

xtotal2=xg+xt+xl1+((xg+xt)*xl1)/(xl2*porcentajefallo);

pmax2=(emodulo*uin�nitamodulo)/xtotal2; pe2=pmax2*sin(delta); plot(delta,pe2)

%curva potencia eléctrica con una sola línea

xtotal3=xg+xt+xl1;

pmax3=(emodulo*uin�nitamodulo)/xtotal3;

pe3=pmax3*sin(delta); plot(delta,pe3);

1.3.3 Cálculo del ángulo crítico

Se van a calcular áreas mediante aproximaciones (función trapz en matlab). Deberemos saber las situación de cada punto en losvectores que representan las curvas. Para ello utilizaremos una función llamada posicion_vector para tal �n.

% cálculo mediante aproximaciones de las integrales. En el anexo está la función calculodeltacritica2lineas.

deltacritica=calculodeltacritica2lineas(delta0,pmax2,pmax3,pmo);


%valor de potencia eléctrica durante el fallo para el ángulo inicial δo

pedelta0curva2=pmax2*sin(delta0);

%posición vector de potencia eléctrica durante el fallo para δo

posdelta0curva2=posicion_vector(dimension,pe2,pedelta0curva2);

%valor de potencia eléctrica durante el fallo para el valor de δcritica

pedeltacriticacurva2=pmax2*sin(deltacritica);

%posición vector de potencia eléctrica durante el fallo para δcritica

posdeltacriticacurva2=posicion_vector(dimension,pe2,pedeltacriticacurva2);

if deltacritica>(pi/2)


posdeltacriticacurva2=dimension-posdeltacriticacurva2;

end


puntos=[pe2(posdelta0curva2:posdeltacriticacurva2),�iplr(pm(posdelta0curva2:posdeltacriticacurva2))];

t1=[delta(posdelta0curva2:posdeltacriticacurva2),�iplr(delta(posdelta0curva2:posdeltacriticacurva2))];

f1=�ll(t1,puntos,1)


pedeltacriticacurva3=pmax3*sin(deltacritica);

posdeltacriticacurva3=posicion_vector(dimension,pe3,pedeltacriticacurva3);

if deltacritica>(pi/2)

. . . posdeltacriticacurva3=dimension-posdeltacriticacurva3;

end

posdelta0curva3=posicion_vector(dimension,pe3,pmo);

posdeltapi0curva3=dimension-posdelta0curva3;

puntos_2=[pm(posdeltacriticacurva3:posdeltapi0curva3),�iplr(pe3(posdeltacriticacurva3:posdeltapi0curva3))];

t1_2=[delta(posdeltacriticacurva3:posdeltapi0curva3),�iplr(delta(posdeltacriticacurva3:posdeltapi0curva3))];


1.3.5 Cálculo del tiempo crítico.

La evolución del ángulo de potencia corresponde a la siguiente ecuación

d2δ

dt2+wo2H

.pmax2.sinδ =wo2H

pmo

Ecuación diferencial, que para resolverla numericamente, tendremos que resolver el sistema siguiente:

y′

1 = y2 = dδdt

d2δdt2 = y

′

2 = wo

2H pmo −wo

2H .pmax2.sinδ

Se escribe in �chero fallo.m:

function dy=fallo(t,y)

dy = zeros(2,1);

wo=2*pi*50; H=6; pmo=0.8; pmax2=0.6784;

% a column vector

dy(1) = y(2);

dy(2) = ((wo)/(2*H))*(pmo-pmax2*sin(y(1)));

Una vez escrito, resolvemos el sistema de ecuaciones diferenciales donde valores iniciales a las variables y1 e y2 de cero y simulandodesde 0 a 1 segundo

[t,y]=ode45(@fallo,[0,1],[0 0]);

% representamos el resultado

plot(t,y)

Para saber el instante en que se produce la δcritica

for k1=2:length(y)

if ((y(k1,1,:))>deltacritica) && ((y(k1-1,1,:))<deltacritica)

salida=k1; break;

end

end

tiempo=salida*(1/length(y))

Chapter 2

Anexos

2.1 Formulación empleada.

La potencia eléctrica suministrada por el generador es:

Pe =E.U∞Xe

senδ

Potencia aceleradora y desaceleradora:

Pace = Pm − Pe Pdec = Pe − Pm

La ecuación de oscilación se puede escribir como:

dδ

dtdw =

wo2H

Pacedδ

El generador será estable si la variación de velocidad en un cierto instante de tiempo del proceso transitorio es nula. Se puedeexpresar como:

ˆ δdec

δ0

Pacedδ =

ˆ δmax

δdes

Pdecdδ Aace = Adec

Ángulo crítico para cuando la función potencia eléctrica no cambia:

ˆ δcri

δ0

(Pm − Pe)dδ =

ˆ Π−δ0

δcri

(Pe − Pm)dδ

La evolución de la velocidad angular es:

ˆ w

w0

dw = w − wo =

ˆ t

0

w0

2H(Pm − Pe)dt

El ángulo del rotor es:

ˆ δ

δo

dδ = δ − δ0 =

ˆ t

o

(w − wo)

9

CHAPTER 2. ANEXOS 10

2.2 Función posicion_vector.

function posicion = posicion_vector(dimension,y,valor)

%dimension: dimensión del vector

%y: vector (será la potencia eléctrica)

%valor: el valor que toma la potencia eléctrica en un punto.

% Se calcula la posición para un valor de ángulo menor de π/2

for k=2:dimension/2;

if (y(k)>valor) && (y(k-1)<valor)

posicion=k

end

end

2.3 Función calculodeltacritica.

function [ solucioncritica ] = calculodeltacritica(pemax,pmo)

delta0=asin(pmo/pemax); % cálculo valor de δ0

intervalo=0.0001; %incremento de δ que se utiliza durante el cálculo

delta_1=delta0; %valor inicial de δ

%cálculo de δcritica

for k=1:1000000

delta_1=delta0+k*intervalo; % incremento de δ

x=linspace(delta_1,pi-delta0,10000); %puntos que se utilizan para representar la curva

f=pemax*sin(x)-pmo; %curva

error=trapz(x,f)-pmo*(delta_1-delta0); %área desaceleración-área aceleradora

if abs(error)<0.001 %Aaceleracion = Adesaceleracion

solucioncritica=delta_1; break;

end

end

end

2.4 Función calculodeltacritica2lineas.

function [ solucioncritica ] = calculodeltacritica2lineas(delta0,pemaxcurva2,pemaxcurva3,pmo)

% cálculo valor de δ0 correspondiente a la curva de potencia eléctrica durante el fallo

delta03=asin(pmo/pemax3);

intervalo=0.0001; %incremento de δ que se utiliza durante el cálculo

delta_1=delta0; %valor inicial de δ

%cálculo de δcritica

for k=1:1000000

delta_1=delta0+k*intervalo; % incremento de δ

xace=linspace(delta0,delta1,10000); %puntos que se utilizan para representar la curva de aceleración

xdes=linspace(delta_1,pi-delta03,10000); %puntos que se utilizan para representar la curva de desaceleración


face=pemax2*sin(xac); %curva aceleración

fdes=pemax3*sin(xdes); %curva desaceleración

areaace=pmo*(delta_1-delta0)-trapz(xace,face); % Área aceleradora

areades=trapz(xdes,fdes)-pmo*(delta_1-delta0); % Área desaceleradora

error=areaace-areades;

if abs(error)<0.001 %Aaceleracion = Adesaceleracion

solucioncritica=delta_1; break;

end

end

end


Datos (valores por unidad):

Grupo H Xg Xt XL1 XL2 cosφ intensidad

1 3 0,3 0,15 0,1 0,2 0,8 12 3,2 0,31 0,14 0,11 0,21 0,75 0,93 3,4 0,32 0,14 0,12 0,2 0,7 0,924 3,6 0,33 0,13 0,13 0,22 0,81 0,995 3,8 0,34 0,12 0,14 0,23 0,83 0,986 4 0,35 0,13 0,11 0,24 0,78 17 2,6 0,36 0,13 0,1 0,22 0,83 0,998 2,8 0,31 0,15 0,12 0,23 0,79 0,959 3,1 0,32 0,12 0,13 0,2 0,84 0,9810 3,3 0,33 0,13 0,14 0,24 0,85 111 3,5 0,34 0,14 0,12 0,25 0,83 0,9312 2,7 0,35 0,16 0,12 0,22 0,82 0,9813 2,9 0,36 0,17 0,11 0,21 0,81 114 2,5 0,37 0,15 0,1 0,2 0,8 1

Obtener y representar:

Caso 1 (cortocircuito-eliminación del cortocircuito):

Fuerza electromotriz: e

Tensión de la red de potencia in�nita: u∞

Ángulo de potencia inicial: δo

Ángulo crítico: δcritica

Ángulo �nal área desaceleradora: δfinal

Tiempo critico: tcritica

Dibujar curvas con las áreas de aceleración y desaceleración.

Caso 2 (cortocircuito al comienzo de la línea 2 - eliminación del cortocircuito mediante apertura de la línea 2):








Caso 3 (cortocircuito en la mitad de la línea 2 - eliminación del cortocircuito mediante apertura de la línea 2):








October 21, 2010¡cticas... · CHAPTER 1. CRITERIO DE LAS ÁREAS. 2 1.1.1 Cálculos eléctricos....

Documents

Transcript of October 21, 2010¡cticas... · CHAPTER 1. CRITERIO DE LAS ÁREAS. 2 1.1.1 Cálculos eléctricos....