Estabilidad Transitoria (Ángulo)

Estabilidad de Sistemas de Potencia

Prof. Agustı́n Rafael Marulanda Guerra, PhD
Escuela Colombiana de Ingenierı́a Julio Garavito
Maestrı́a en Ingenierı́a Eléctrica

Octubre de 2018
Contenido

1 Introducción

2 Ecuación de Oscilación

3 Modelo Clásico de la Máquina Sı́ncrona

4 Criterio de Igualdad de Áreas

5 Solución por Integración Numérica

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 2 / 89

Introducción

Estabilidad Transitoria
El problema de estabilidad transitoria trata de la capacidad de los
generadores de seguir funcionando en sincronismo después de la ocurrencia
de una perturbación

Que los generadores se encuentren en sincronismo implica que las diferencias

de los ángulos eléctricos entre ellos permanece constante, es decir; sus
velocidades eléctricas son iguales (sincronismo)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 3 / 89

Introducción

Estabilidad Transitoria
Es un problema de estabilidad de gran perturbación; l magnitud de la
perturbación que tiene lugar es tal que las ecuaciones diferenciales que
describen el comportamiento dinámico del sistema no puede ser linealizado
para su análisis

Entra en el problema de estabilidad a corto plazo, ya que considera que las

dinámicas dominantes son las de los generadores sı́ncronos y, a lo sumo la de
sus sistemas de control primario (tensión y carga-velocidad)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 4 / 89

Introducción

Estabilidad Transitoria
En estudios de estabilidad transitoria el periodo de interés usualmente se
limita entre 3 y 5 segundos luego de la perturbación, aunque puede
extenderse alrededor de 10 segundos para sistemas muy grandes don modos
de oscilación dominante inter-areas

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 5 / 89

Introducción

Perturbación
En estudios de estabilidad interesa analizar el comportamiento del sistema
cuando éste ha sido sometio a una perturbación transitoria

La perturbación puede ser pequeña o grande:

Pequeña: Pequeñas variaciones en la carga o generación
Grandes: Cortocircuitos, perdida de carga, perdida de generación, o cambio brusco de la
potencia mecánica

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 6 / 89

Introducción

Pequeñas Perturbación
Las pequeñas perturbaciones de la carga ocurren continuamente en el
sistema, a pesar de esto; el sistema se ajusta por si mismo a las condiciones
cambiante

El sistema debe ser capaz de operar satisfactoriamente bajo esas condiciones

y satisfacer el máximo valor de la carga

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 7 / 89

Introducción

Gandes Perturbación
La estabilidad de un sistema viene determinada por la capacidad de sobrevivir
a numerosas perturbaciones de una naturaleza severa, tales como:
Cortocircuitos en una lı́nea de transmisión
Perdida de un generador de gran tamaño
Perdida de carga de gran tamaño
Perdida de una lı́nea de interconexión entre dos subestaciones
Cambio brusco de la potencia mecánica

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 8 / 89

Introducción

Comportamiento del ángulo del rotor para caso estable e inestable

2.5

2
Ángulo [rad.]

1.5

0.5

0
0 0.5 1 1.5 2 2.5 3 3.5
Tiempo [s.]

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 9 / 89

Ecuación de Oscilación
La ecuación de oscilación relaciona el movimiento rotatorio de la máquina
con el torque de inercia, resultante de la interacción de los torques eléctrico y
mecánicos en el rotor.

d2 θ
J = Ta [N.m]
dt2

J : Momento de inercia [Kg.m2 ]

θ : Ángulo mecánico del eje [rad]
Ta : Torque de aceleración [N.m]

Para un generador, el torque eléctrico, Te , es opuesto al torque mecánico Tm

de entrada;

d2 θ
J = Tm − Te [N.m]
dt2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 10 / 89

Ecuación de Oscilación

El ángulo mecánico se establece con una referencia sı́crona rotatoria,

moviendo a velocidad constante;

θ = ωt + δm [rad]
dθ dδm
=ω+
dt dt
2 2
d θ d δm
=
dt2 dt2

Entonces, la ecuación de oscilación se puede rescribir como;

d2 δ
J = Tm − Te [N.m]
dt2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 11 / 89

Ecuación de Oscilación
Si multiplicamos ambos lados por ωm , tenemos;

d2 δ
Jωm = ωm Tm − ωm Te [w]
dt2
2
d δ
Jωm 2 = Pm − Pe [w]
dt

Donde; M = Jωm se le denomina constante de inercia. Esta costante está

relacionada con la energı́a cinética de la masa rotacional del eje del
generador, Wk = 1/2Jωm . Por lo tanto;
2Wk
Jωm =
ωm
Entonces, la ecuación de oscilación se puede escribir como;

2Wk d2 δ
= Pm − Pe [w]
ωm dt2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 12 / 89

Ecuación de Oscilación
Se puede dividir por la potencia base para expresarlo en PU;

1 2Wk d2 δ Pm Pe
= −
Sbase ωm dt2 Sbase Sbase
2 2Wk d2 δ
= Pm (pu) − Pe (pu)
ωm Sbase dt2

Si definimos la constante de inercia de la máquina como;

Wk Energia almacenada a velocidad nominal
H= =
Sbase P otencia nomimal de la maquina

Entonces, la ecuación fundamental de la dinámica del sistema viene dada por;

2H d2 δ
= Pm (pu) − Pe (pu)
ωm dt2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 13 / 89

Ecuación de Oscilación

2H d2 δ
= Pm (pu) − Pe (pu)
ωm dt2

La ecuación de oscilación se puede resolver por medio de un sistema de

ecuaciones diferenciales de primer orden;

2H dωm
= Pm (pu) − Pe (pu)
ωm dt
dδ
= ωm − ωs
dt

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 14 / 89

Ecuación de Oscilación

ω y δ deben tener unidades consistentes; radianes (grados) eléctricos ó

P
mecánicos. δe = δm . Sı́ ω = 2πf , entonces;
2

H d2 δ
= Pm (pu) − Pe (pu) δ en radianes
πf dt2

H d2 δ
= Pm (pu) − Pe (pu) δ en grados
180f dt2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 15 / 89

Modelo Clásico de la Máquina Sı́ncrona
Consideraciones para el modelo:

1 Para estudios en el cual el perı́odo de análisis es pequeño en comparación con

0 0
Td0 se puede asumir que Eq se mantiene constante durante el perı́odo de
estudio.

2 Se ignora la saliencia asumiendo que las reactancias transitorias de eje directo

0 0
y eje en cuadratura son iguales, Xd = Xq . Esto implica un flujo concatenado
constante.
0 0
3 Xd = Xq también implica que la tensión detrás de la reactancia transitoria
0
ra + jXd tiene una magnitud constante.

4 Devanados de amortiguación despreciables y potencia mecánica constante.

5 Rotor del generador y eje del motor primario forman un sólido rı́gido

6 Resistencias del devanado del estator y rotor despreciables

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 16 / 89

Modelo Clásico de la Máquina Sı́ncrona

Ll L fd

Id Ifd
El flujo concatenado en pu en el directo
viene dado por; d ad Lad fd

ψad = −Lad id + Lad if d

ψd = ψad − Ll id Figura: Equivalente del eje directo
ψf d = ψad + Lf d if d

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 17 / 89

Modelo Clásico de la Máquina Sı́ncrona

Ll L fd
De la ecuación de flujo de campo
despejamos la corriente de campo; Id Ifd

d ad Lad
ψf d − ψad fd
if d =
Lf d

Sustituyendo en ψad , tenemos; Figura: Equivalente del eje directo

Lad
ψad = −Lad id + (ψf d − ψad )
Lf d

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 18 / 89

Modelo Clásico de la Máquina Sı́ncrona

Despejando ψad , tenemos;

Ll L fd
 
0 ψf d Id Ifd
ψad = Lad −id +
Lf d
d ad Lad fd
Donde:
0 1 0
Lad = = Ld − Ll
1 1 Figura: Equivalente de eje directo
−
Lad Lf d

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 19 / 89

Modelo Clásico de la Máquina Sı́ncrona

Ll L1q
Similarmente para el eje en cuadratura;
  Iq I1q
0 ψ1q
ψaq = Laq −iq +
L1q q aq Laq 1q

Donde:
0 0
Laq = Ld − Ll
Figura: Equivalente para el eje en q

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 20 / 89

Modelo Clásico de la Máquina Sı́ncrona
La tensión de eje directo viene dada por;

ed = −ra id − ωψq
= −ra id + ω (Ll iq − ψaq )

Sustituyendo ψaq , tenemos;

 
0 ψ1q
ed = −ra id + ωLl iq − ωLaq −iq +
L1q
 
 0
 0 ψ1q
= −ra id + ω Ll + Laq iq − ωLaq
L1q
0 0
= −ra id + Xq iq + Ed
donde;  
0 0 ψ1q
Ed = −ωLaq
L1q

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 21 / 89

Modelo Clásico de la Máquina Sı́ncrona
Similarmente para el eje en cuadratura;
0 0
ed = −ra iq − Xd id + Eq
donde;  
0 0 ψf d
Eq = −ωLad
Lf d
0 0
Despreciando la saliencia, Xd = Xq , la tensión en los terminales del estator
es:
 0 0
 0
ed + jeq = Ed + jEq − ra (id + jiq ) + Xd (iq − jid )
 0 0
 0
= Ed + jEq − ra (id + jiq ) − jXd (id + jiq )
0
 0

Et = E − ra + jXd It

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 22 / 89

Modelo Clásico de la Máquina Sı́ncrona

Entonces, el modelo clásico de la máquina sı́ncrona viene dado por;

0
 0

Et = E − ra + jXd It

Despreciando la resistencia, tenemos;

0 0
Et = E − jXd It

donde,
0 0 0
E = Ed + jEd
 
ψ1q ψf d
= L0 ad − +j
L1q Lf d

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 23 / 89

Modelo Clásico de la Máquina Sı́ncrona

I
q

E'

E'q E'd

R
E'R

Figura: Modelo clásico de la máquina

sı́ncrona Figura: Representación fasorial

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 24 / 89

Modelo Clásico de la Máquina Sı́ncrona

La tensión detrás
0
de la reactancia transitoria
está definida E

El ángulo del rotor δ representa el ángulo en

0 0
el cual E adelanta a Et I
q

0 E'
La magnitud de E se puede determinar d

calculado su valor antes de la perturbación

E'q E'd

R
E'R

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 25 / 89

Modelo Clásico de la Máquina Sı́ncrona

Cuando el sistema tiene una perturbación, la

0
magnitud de E permanece constante al valor
previo de la perturbación, y δ cambia en
proporción a la desviación de la velocidad del q
I

rotor con respecto a la velocidad de sincronismo E'

E'q E'd

R
E'R

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 26 / 89

Modelo Clásico de la Máquina Sı́ncrona

Si se desprecia la resistencia, la potencia

eléctrica entregada por la máquina viene
dada por;
0
E Et
Pe = 0 Senδ = Pmax senδ
Xd

donde;
0
E Et
Pmax = 0
Xd

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 27 / 89

Modelo Clásico de la Máquina Sı́ncrona

Dado que se desprecia la resistencia del estator, la Pe representa tanto la

potencia del entrehierro como la potencia en terminales de la máquina

Caracteristica Potencia-Ángulo
2.5

2
Potencia normalizada (Pr/Po)

1.5

0.5

0
0 0.5 1 1.5 2 2.5 3 3.5
Ángulo de la carga

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 28 / 89

Modelo Clásico de la Máquina Sı́ncrona

Durante una perturbación, la oscilación de δ está superpuesta a la velocidad

de sincronismo ω0 , pero la desviación de la velocidad (∆ωr = dδ/dt) es
mucho más pequeña que ω0 .

La velocidad del generador es prácticamente igual a ω0 , y el torque en pu del

entrehierrro se puede considerar igual a la potencia en pu del entrehierro.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 29 / 89

Modelo Clásico de la Máquina Sı́ncrona

Entonces, la ecuación de oscilación puede escribirse como;

2H d2 δ
= Pm − Pmax senδ
ω0 dt2
d2 δ ω0
= [Pm − Pmax senδ]
dt2 2H

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 30 / 89

Variación en la Potencia Mecánica

Entonces, la ecuación de oscilación puede escribirse como;

2H d2 δ
= Pm − Pmax senδ
ω0 dt2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 31 / 89

Variación en la Potencia Mecánica
P

P e= P sen
max

c
Área A 1 Área A 2
b
Pm1

Pm0 a

0 1 m

t(s)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 32 / 89

Criterio de igualdad de áreas

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

d2 δ ω0
= (Pm − Pe )
dt2 2H
Pe es una función no lı́neal de δ, y por lo tanto no es posible resolverla de
forma directa. Si multiplicamos ambos lados por dδ/dt, tenemos;

dδ d2 δ ω0 dδ
2
= (Pm − Pe )
dt dt 2H dt
El lado izquierdo de la ecuación se puede trabajar como la derivada del
cuadrado de la frecuencia del sistema,
 2
dδ dδ ω0 dδ
= (Pm − Pe )
dt dt 2H dt

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 33 / 89

Criterio de igualdad de áreas

Integrando a ambos lados, tenemos;

P
 2 Z
dδ ω0
= (Pm − Pe ) dδ P e= P
dt 2H c
max
sen

Área A 1 Área A 2
b
Pm1

Pm0 a

Inicialmente la velocidad de la desviación

dδ
ω= es igual a cero, ésta puede varı́a 0 1 m

dt
como resultado de una perturbación

Para una operación estable, la desviación del

ángulo δ debe estar entre unos lı́mites,
alcanzando un valor máximo y luego t(s)

cambiando de dirección

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 34 / 89

Criterio de igualdad de áreas

Para que δ cambie de dirección es necesario P

dδ
que la velocidad de desviación ω = sea P e= P
max
sen

dt c
igual a cero en algún momento después de Área A 1
b
Área A 2

Pm1
la perturbación. Por lo tanto;
Pm0 a

Z δm
ω0
(Pm − Pe ) dδ = 0
δ0 2H 0 1 m

Esto indica que el área entre Pm y Pe debe

ser igual a cero para que el sistema sea
estable t(s)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 35 / 89

Criterio de igualdad de áreas
Para que ocurra la estabilidad después de la perturbación, la aceleración debe
ser cero;
 2 Z δm
dδ ω0
=0⇒ (Pm − Pe ) dδ = 0
dt 2H δ0

P e= P sen
max

c
Área A 1 Área A 2
b
Pm1

Pm0 a

0 1 m

t(s)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 36 / 89

Criterio de igualdad de áreas

Paraque la acelación sea igual a cero, el área

A1 debe ser igual al área A2
P

El área A1 representa la anergı́a cinética

P e= P sen
max
ganada por el rotor durante la aceleración c

cuando δ cambia de δ0 a δ1 ; Pm1

Área A 1
b
Área A 2

Z δ1 Pm0 a

A1 = (Pm − Pe ) dδ
δ0
0 1 m

El área A2 representa la energı́a cinética

perdida por el rotor durante la
desaceleración cuando δ cambia de δ1 a δm ;
t(s)
Z δm
A2 = (Pe − Pm ) dδ
δ1

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 37 / 89

Criterio de igualdad de áreas

El criterio de igual de áreas permite estudiar las oscilaciones de una máquina

con respecto a una barra infinita, sin necesidad de recurrir a la solución
numérica de la ecuación de oscilación.

Si el sistema es inestable δ continua en aumento de manera indefinida con

respecto al tiempo, y la máquina pierde sincronismo.

Si el sistema es estable δ tiene oscilaciones (no senoidales) cuyas amplitudes

disminuyen, en la práctica a causa de los términos de amortiguamiento que
no se incluyen en la ecuación de oscilación.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 38 / 89

Criterio de igualdad de áreas

Los casos tı́picos que se estudian con este criterio son:

Cambio repentino en la entrada mecánica.

Efecto del tiempo de desconexión luego de una falla.

Pérdida repentina de una lı́nea en paralelo.

Cortocircuito repentino en una lı́nea.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 39 / 89

Ángulo Crı́tico y Tiempo Crı́tico

Se denomina ángulo crı́tico, δc , de despeje de falla al valor lı́mite del ángulo

para el cual deben actuar las protecciones y garantizar la permanencia en
sincronismo.
Su valor se obtiene de aplicar el criterio de igualdad de áreas. Estro es;
 2 Z δm
dδ ω0
=0⇒ (Pm − Pe ) dδ = 0
dt 2H δ0

Y el instante en que se alcanza el ángulo crı́tico se denomina tiempo crı́tico y

se obtiene al integral la ecuación,

2H d2 δ
= Pm − Pe = Pa
ω0 dt2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 40 / 89

Ángulo Crı́tico y Tiempo Crı́tico

Dado que durante la falla la potencia eléctrica, Pe, es nula entonces Pa=Pm;
el torque de aceleración y la aceleración son estrictamente constante. Ası́,
partiendo de la ecuación de oscilación:

2H d2 δ
= Pm − Pe = Pa
ω0 dt2
Dado que la potencia eléctrica es igual a cero, entonces;

2H d2 δ
= Pm
ω0 dt2
Integrado con respecto al tiempo
Z t
dδ ω0 ω0
= Pm dt = Pm t
dt 2H 0 2H

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 41 / 89

Ángulo Crı́tico y Tiempo Crı́tico

Integrando nuevamente, tenemos

ω0
δ= P m t2
2H
El comportamiento del ángulo del rotor se define por:
ω0
δ = δ0 + Pm t 2
2H
De esta forma, el águlo crı́tco viene dado por:

ω0
δc = δ0 + Pm t2c
2H

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 42 / 89

Ángulo Crı́tico y Tiempo Crı́tico

Mientras que el tiempo crı́tico se obtiene conociendo el valor de δc ,

ω0
δc = δ0 + Pm t2c
2H
Despejando tc , tenemos;
s
2H (δc − δ0 ) ω0
tc = δ0 + Pm t2c
ω 0 Pm 2H

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 43 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

En el sistema mostrado en la figura, el generador está entregando una

potencia activa de 1 pu a través de transformador elevador y dos lı́neas en
paralelo a una barra infinita de tensión 1 pu. La tensión en el lado de alta del
transformador se regula a un valor de 1.05 pu.
Se produce una situación de falla trifásica de una de las lı́neas que parte del
transformador, inmediatamente después de la barra, despejado por la
protección mediante la apertura de la lı́nea.
Se desea conocer el tiempo crı́tico de despeje de la falla utilizando el modelo
clásico de la máquina

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 44 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

En el sistema mostrado en la figura, el generdor está entregando una

potencia activa de 1 pu a través de transformador elevador y dos lı́neas en
paralelo a una barra infinita de tensión 1 pu. La tensión en el lado de alta del
transformador se regula a un valor de 1.05 pu.
Se produce una situación de falla trifásica de una de las lı́neas que parte del
transformador, inmediatamente después de la barra, despejado por la
protección mediante la apertura de la lı́nea.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 45 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Se desea conocer el tiempo crı́tico de despeje de la falla utilizando el modelo

clásico de la máquina

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 46 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Los parámetros del sistema se indican en la siguiente tabla

Tabla: Parámetros del sistema

Elemento Parámetro
0
Generador Xd = 0,15 pu, H = 3 s
Transformador Xcc = 0,1 pu
Lı́neas XL = 0,6 pu (cada una)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 47 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

A partir de las ecuaciones de flujo se puede

obtener la corriente inyectada en la barra:
V2 V1
P = sen(θ2 − θ1 )
0,3

1,05
P = senθ2 −→ θ2 = 16,602 grados
0,3

V2 V1
Q= cos(θ2 − θ1 )
0,3

1,05
Q= cos16,602 −→ Q = 0,021 pu
0,3

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 48 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

En el sistema mostrado en la figura, el generdor está entregando una

potencia activa de 1 pu a través de transformador elevador y dos lı́neas en
paralelo a una barra infinita de tensión 1 pu. La tensión en el lado de alta del
transformador se regula a un valor de 1.05 pu.
Se produce una situación de falla trifásica de una de las lı́neas que parte del
transformador, inmediatamente después de la barra, despejado por la
protección mediante la apertura de la lı́nea.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 49 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Se desea conocer el tiempo crı́tico de despeje de la falla utilizando el modelo

clásico de la máquina

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 50 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Los parámetros del sistema se indican en la siguiente tabla

Tabla: Parámetros del sistema

Elemento Parámetro
0
Generador Xd = 0,15 pu, H = 3 s
Transformador Xcc = 0,1 pu
Lı́neas XL = 0,6 pu (cada una)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 51 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

A partir de las ecuaciones de flujo se puede

obtener la corriente inyectada en la barra:
V2 V1
P = sen(θ2 − θ1 )
0,3

1,05
P = senθ2 −→ θ2 = 16,602 grados
0,3

V2 V1
Q= cos(θ2 − θ1 )
0,3

1,05
Q= cos16,602 −→ Q = 0,021 pu
0,3

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 52 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

La corriente del circuito viene dada por la potencia;

P − jQ
It =
V1∗
1 − j0,021
=
1 − j0
= 1∠ − 1,2
La tensión en los terminales del generador, viene dada por;
Vt = It (Xcc + XL /2) + V1
= 1∠ − 1,2 (j0,1 + j0,3) + 1∠0
= 1,085∠21,63

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 53 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

La tensión interna del generador se puede

escribir como;

0 0
E = Vt + It Xd
= 1,085∠21,63 + 1,0∠ − 1,2(j0,15)
= 1,151∠28,53 pu

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 54 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

La Curva de Potencia Eléctrica antes de la falla, es;

0
E V1
Pe = sen(δ − θ1 )
(0,15 + 0,1 + 0,3)
= 2,093senδ

Donde antes de la falla δo = 28,53 y Pa = Pm − Pe = 0

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 55 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Cuando se produce la falla, la potencia eléctrica se hace cero manteniéndose

la potencia mecánica en el eje del generador, por lo tanto; la potencia de
aceleración es igual a la potencia mecánica Pa = Pm − Pe = Pm

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 56 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Al despejar la falla, se abre una de las lı́neas y la potencia eléctrica posfalla

viene dada por la nueva configuración del sistema, es decir;

0
E V1
Pepf = sen(δ − θ1 )
(0,15 + 0,1 + 0,6)
= 1,354senδ

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 57 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

P 2.093 sen

1.354 sen

Área A2

Pm

Área A1

0 1 c m

Figura: Curvas de potencia-ángulo pre y posfalla

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 58 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

El ángulo crı́tico se determina como el valor

de δ que iguale el área 1 con el área 2 P 2.093 sen

El valor del área 1 viene dada por la 1.354 sen

diferencia entre el ángulo crı́tico y el ángulo Área A2

Pm
inicial por la potencia mecánica:
Área A1
A1 = Pm (δc − δo )
= δc − 0,4984 0 1 c m

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 59 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

El área 2 viene dada por la curva de

potencia-ángulo después de la falla y la
potencia mecánica, esto es;

P 2.093 sen
Z δm
A2 = (1,354senδ − Pm ) dδ 1.354 sen
δc
Área A2
Z δm Z δm
Pm
= 1,354 senδdδ − Pm dδ
δc δc
Área A1

= 1,354 (Cosδc − Cosδm ) − (δm − δc )

0 1 c m

A2 = 1,354Cosδc − 1,354Cosδm − δm + δc

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 60 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

δ1 , viene dado por la ángulo que iguala

la potencia eléctrica con la mecánica
después que la falla ha sido despejada.
Esto es;
P 2.093 sen
Pm = 1,354senδ1
 
1
δ1 = 1,354Sen−1
1.354 sen

1,354 Área A2

Pm
δ1 = 0,831 rad
Área A1
Cómo δm = π − δ1 , entonces;
1 c m
δm = π − δ1 0

= π − 0,831
= 2,311 rad

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 61 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Sustituyendo δ1 en la ecuación del área 2,

tenemos;

A2 = 1,354Cosδc − 1,354Cosδm − δm + δc P 2.093 sen

= 1,354Cosδc + δc − 1,354Cos2,311 − 2,311

1.354 sen
= 1,354Cosδc + δc − 1,398 Área A2

Pm

Área A1
El ángulo crı́tico corresponde al valor de δ
cuando el área de aceleración es igual al 0 1 c m

área de desacelaración. Esto es;

A1 = A2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 62 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Entonces, el ángulo crı́tico viene dado por;

P 2.093 sen

A1 = A2
1.354 sen
δc − 0,4984 = 1,354Cosδc + δc − 1,398
Área A2

−0,4984 = 1,354Cosδc − 1,398 Pm

0,6646 = Cosδc Área A1

δc = 0,844 rad
0 1 c m

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 63 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

El comportamiento de δ viene dado por la

ecuación dinámica del rotor;

2H d2 δ
= Pm − P e
ωo dt2
2.093 sen
d2 δ ωo P

= Pm − Pe
dt2 2H Z 1.354 sen
dδ ωo
= Pm dt Área A2

dt 2H Pm

dδ 377
= t Área A1
dt 2(3)
Z
0 1 c m

δ = 62,8 tdt
1
δ= 62,8t2
2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 64 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Como δo = 0,498 rad, el comportamiento de

δ se puede expresar como; P 2.093 sen

1.354 sen
1 ωo 2
δ = δo + t Área A2

2 2H Pm

1
= 0,498 + 62,8t2 Área A1
2
0 1 c m

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 65 / 89

 Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

El tiempo crı́tico se puede conocer a partir

de la ecuación de δ, conociendo δc . Esto es;
1
δc = 0,498 + 62,8t2c 2.093 sen
2 P

1 2
δc − 0,498 = 62,8tc 1.354 sen
2
s Área A2

2(δc − 0,498) Pm

tc =
62,8 Área A1
s
2(0,844 − 0,498)
tc = 0 1 c m

62,8
tc = 0,105 s.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 66 / 89

Solución por Integración Numérica

La ecuación dinámica del rotor viene dada por;

d2 δ(t) ωo
2
= (Pm − Pe )
dt 2H
Para describir el comportamiento de ángulo del roto δ(t), se debe resolver la
ecuación diferencial que describe la dinámica del rotor

d2 δ(t) ωo
= (Pm − Pmax Sen δ)
dt2 2H
δ(0) = δ0

Se utiliza la integración numérica para resolver la ecuación diferencial

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 67 / 89

Solución por Integración Numérica

Generalmente, en el análisis de sistemas de potencia, las ecuaciones

deferenciales a resolver son ecuaciones diferenciales no lineales de donde se
conoce sus valores iniciales,
dx
= f (x, t)
dt
x: Es el vector de estados de n variables dependientes
t: Variable independiente (tiempo)

EL objetivo es calcular a x como una función de t, conociendo sus respectivos

valores iniciales de x = x0 y t = t0 .

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 68 / 89

Solución por Integración Numérica
Método de Euler

Considere la ecuación diferencial de primer orden

dx
= f (x, t)
dt
con x = x0 en t = t0 , tal como se muestra en la figura

x
Solución real

Tangente

x1
x
x0

t
t
t0 t1

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 69 / 89

Solución por Integración Numérica
Método de Euler

En el punto p(x0 , t0 ) se puede aproximar la

curva representando la solución por su
tangente y su pendiente; x
Solución real

dx
|x=x0 = f (x0 , t0 )
dt
Tangente
Ası́, el incremento en x se obtiene como;
x1
x
dx x0
∆x =|x=x0 ∆t
dt
Por lo tanto, el valor de x en t = t1 , (donde t
t1 = t0 + ∆t) está dado por; t
t0 t1

dx
x1 = x0 + ∆x = x0 + |x=x0 ∆t
dt

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 70 / 89

Solución por Integración Numérica
Método de Euler

El método de Euler es equivalente a usar la

expansión de Taylor de primer grado para x
en el punto pto(x0 , t0 )
x
Solución real

∆t2 ∆t3 ...

x1 = x0 +∆t (ẋ0 )+ (ẍ0 )+ ( x 0 )+· · ·
2! 3!
Tangente
Luego de calcular x1 para t = t1 , podemos x1
x
tomar otro valor de t usando el incremento x0
del tiempo: t2 = t1 + ∆t. Esto es;
dx t
x2 = x1 +|x=x1 ∆t
dt t0 t1
t

Aplicando esta técnica sucesivamente se

puede determinar el valor correspondiente de
x para diferentes valores de t.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 71 / 89

Solución por Integración Numérica
Método de Euler

La forma generaliza dell método de Euler se

puede escribir como;
x
Solución real

dx
xn+1 = xn + |x=xn ∆t
dt
Tangente

x1
El método sólo considera la primera deriva x
x0
de x (método de primer orden)

Para reducir el error entre el valor real y el t

t
calculado, el paso (∆t) debe ser lo t0 t1

suficientemente pequeño para que el ajuste

se aproxime a la curva real

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 72 / 89

Solución por Integración Numérica
La ecuación dinámica del rotor viene dada por;
d2 δ(t) ωo
2
= (Pm − Pmax Sen δ)
dt 2H
Como el método de Euler se aplica a ecuaciones diferenciales de primer orden;
dω(t) ωo
= (Pm − Pmax Sen δ)
dt 2H
dδ(t)
= ω(t)
dt

donde;
δ(0) = δ0
ω(0) = 0

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 73 / 89

Solución por Integración Numérica
Ejemplo

En el sistema mostrado en la figura, el generdor está entregando una

potencia activa de 1 pu a través de transformador elevador y dos lı́neas en
paralelo a una barra infinita de tensión 1 pu. La tensión en el lado de alta del
transformador se regula a un valor de 1.05 pu.
Se produce una situación de falla trifásica de una de las lı́neas que parte del
transformador, inmediatamente después de la barra, despejado por la
protección mediante la apertura de la lı́nea.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 74 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Los parámetros del sistema se indican en la siguiente tabla

Tabla: Parámetros del sistema

Elemento Parámetro
0
Generador Xd = 0,15 pu, H = 3 s
Transformador Xcc = 0,1 pu
Lı́neas XL = 0,6 pu (cada una)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 75 / 89

Solución por Integración Numérica
Ejemplo

Se desea conocer el tiempo crı́tico de despeje de la falla utlizando la

integración numérica (Método de Euler)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 76 / 89

Solución por Integración Numérica
Ejemplo

Según lo calculado en el ejercicio previo;

Pepref = 2,093Sen δ P 2.093 sen

Pef alla =0
1.354 sen

Peposf = 1,354Sen δ Área A2

Pm

Área A1

En este caso, el sistema se presenta por tres

ecuaciones de oscilación. Con H = 3 y 0 1 c m

Pm = 1:

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 77 / 89

Solución por Integración Numérica
Ejemplo

Antes de la falla (t < 0);

dω(t) ωo
= (Pm − Pmax Sen δ)
dt 2H
dω(t)
= 62,832 − 131,51Sen δ
dt
P 2.093 sen

En el momento de la falla (o < t < tc ) 1.354 sen

Área A2

dω(t) Pm

= 62,832
dt Área A1

0 1 c m

Luego de la Falla (t > tc )

dω(t)
= 62,832 − 85,07Sen δ
dt
©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 78 / 89
Solución por Integración Numérica
Ejemplo

Deseamos simular el comportamiento de δ

por un intervalo de 1 s con 10 muestras. Por
lo tanto;
P 2.093 sen

tf − t0
∆t = = 0,1
n 1.354 sen

Área A2
Antes de la falla PA = 0 ⇒ Pm = Pe . Pm

Entoces;
Área A1
 
Pm
δ0 = arcsen 0 1 c m

Pmec
= 0,4981 rad

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 79 / 89

Solución por Integración Numérica
Ejemplo

Para t0 = 0; el ángulo inicial viene dado por;

δ(t = 0) = δ0 = 0,4981 rad 2.093 sen

P

= 0,4981 rad
1.354 sen

Área A2

Mientras que la velocidad inicial Pm

corresponde a;
Área A1

dδ0 0 1 c m

ω(t = 0) =
dt
= 0 rad/s.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 80 / 89

Solución por Integración Numérica
Ejemplo

Para t1 ;
dδ(t)
δ1 = δ0 + ∆t |δ=δ0
dt
2.093 sen
= δ0 + ∆tω0 P

= 0,4981 + 0,1 [0] 1.354 sen

= 0,4981 rad Área A2

Pm

Área A1
dω(t)
ω1 = ω0 + ∆t |ω=ω0
dt 0 1 c m

= ω0 + ∆t [62,833 − 85,07senδ0 ]
= 2,2188 rad/s.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 81 / 89

Solución por Integración Numérica
Ejemplo

Para t2 ;
dδ(t)
δ2 = δ1 + ∆t |δ=δ1
dt
2.093 sen
= δ1 + ∆tω1 P

= 0,4981 + 0,1 [2,2188] 1.354 sen

= 0,72 rad Área A2

Pm

Área A1
dω(t)
ω2 = ω1 + ∆t |ω=ω1
dt 0 1 c m

= ω0 + ∆t [62,833 − 85,07senδ1 ]
= 4,4378 rad/s.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 82 / 89

Solución por Integración Numérica
Ejemplo

Para t3 ;
dδ(t)
δ3 = δ2 + ∆t |δ=δ2
dt
2.093 sen
= δ2 + ∆tω2 P

= 0,72 + 0,1 [4,4378] 1.354 sen

= 1,1638 rad Área A2

Pm

Área A1
dω(t)
ω3 = ω2 + ∆t |ω=ω2
dt 0 1 c m

= ω2 + ∆t [62,833 − 85,07senδ2 ]
= 5,1117 rad/s.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 83 / 89

Solución por Integración Numérica
Ejemplo

Se continua el proceso de integración

numérica, hasta completar el tiempo de
simulación, con lo cual se obtienen los
2.093 sen
siguientes valores; P

1 t delta ( t ) omega ( t )
1.354 sen
2
3 0.0000 0.4981 0.0000 Área A2
4 0.1000 0.4981 2.2185
5 0.2000 0.7200 4.4370 Pm
6 0.3000 1.1637 5.1106
7 0.4000 1.6747 3.5817
8 0.5000 2.0329 1.4034 Área A1
9 0.6000 2.1732 0.0715
10 0.7000 2.1804 −0.6550
11 0.8000 2.1149 −1.3469 0 1 c m

12 0.9000 1.9802 −2.3426

13 1.0000 1.7459 −3.8638

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 84 / 89

Solución por Integración Numérica
Ejemplo

Comportamiento de δ y ω con ∆t = 0,1

Curva de Ángulo y Velocidad para tiempo critico de 0.01 s.

2.5

2
Ángulo [rad.]

1.5

0.5

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Tiempo [s.]

4
Velocidad [rad/s.]

-2

-4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Tiempo [s.]

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 85 / 89

Solución por Integración Numérica
Ejemplo

Comportamiento de δ y ω con ∆t = 0,001

Curva de Ángulo y Velocidad para tiempo critico de 0.01 s.

1.4

1.2
Ángulo [rad.]

0.8

0.6

0.4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Tiempo [s.]

2
Velocidad [rad/s.]

-1

-2

-3
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Tiempo [s.]

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 86 / 89

Solución por Integración Numérica
Ejemplo

Código Matlab
1 clear all ;
2 t o =0; % Tiempo i n i c i a l
3 t ( 1 )=t o ;
4 tn =1.0; % Tiempo f i n a l de s i m u l a c i o n
5 Tc = 0 . 0 1 ; % Tiempo c r i t i c o
6 n =10; % S u b i n t e r v a l o s de i n t e g r a c i o n
7 D t =(tn−t o ) /n ; % Paso de i n t e g r a c i o n
8 Pmaxo = 2 . 0 9 3 ; % P o t e n c i a maxima p r e f a l l a
9 Pmec=1; % Potencia macanica
10 d e l t a 0=a s i n ( Pmec/Pmaxo ) ; % Condicion i n i c i a l del angulo
11 v e l 0= 0 ; % V e l o c i d a d i n i c i a l de l a maquina
12 H=3; % C o n s t a n t e de i n e r c i a
13 f =60; % F r e c u e n c i a Hz
14 Pmaxpref = 2 . 0 9 3 ; % P o t e n c i a máxima P r e f a l l a
15 Pmaxf =0; % P o t e n c i a maxima d u r a n t e l a f a l l a
16 Pmaxposf = 1 . 3 5 4 ; % P o t e n c i a maxima p o s t f a l l a
17 t ( 1 )=t o ; omega ( 1 )=v e l 0 ;
18 d e l t a ( 1 )=d e l t a 0 ;
19 f o r k =1: n
20 t ( k+1)=t ( k )+D t ;
21 d e l t a ( k+1)=d e l t a ( k )+ D t * omega ( k ) ;
22 i f t ( k+1)<Tc
23 Pmax=Pmaxf ;
24 end
25 i f ( t ( k+1)>=Tc )
26 Pmax=Pmaxposf ;
27 end
28 omega ( k+1)=omega ( k ) +( p i * f /H* (Pmec−Pmax* s i n ( d e l t a ( k ) ) ) ) * D t ;
29 end

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 87 / 89

Solución por Integración Numérica
Ejemplo

Comportamiento de δ y ω con ∆t = 0,001, para distintos valor de tiempo

crı́tico.

Curva de Ángulo y Velocidad para diferentes tiempos de despeje

10

t =0.135
8 c
Ángulo [rad.]

0
0 0.2 0.4 0.6 0.8 1 1.2
Tiempo [s.]

30
Velocidad [rad/s.]

t =0.135
20 c

10

-10
0 0.2 0.4 0.6 0.8 1 1.2
Tiempo [s.]

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 88 / 89

Solución por Integración Numérica
Ejemplo

Comportamiento de δ y ω con ∆t = 0,001, para distintos valor de tiempo

crı́tico.

Curva de Ángulo y Velocidad para diferentes tiempos de despeje

70

60

50
Ángulo [rad.]

40

30

20

10

0
0 0.5 1 1.5 2 2.5
Tiempo [s.]

100

80
Velocidad [rad/s.]

60

40

20

-20
0 0.5 1 1.5 2 2.5
Tiempo [s.]

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 89 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

En el sistema de la figura, la barra 3 es de potencia infinita y tensión 1 pu,

mientras que la barra 2 tiene una tensión de 1,01 + j0,2. La potencia
demanda por la carga conectada en la barra 2 es de 1 + j0,4

Determine si el sistema es estable al desconectar bruscamente toda la carga

de barra 2, quedando el transformador y la lı́nea conectados
Si en la situación anterior, se desconecta permanentemente la lı́nea, es
estable el sistema?

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 90 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Antes de la perturbación, las corrientes por

el sistema vienen dadas por:
V2 − V3
IL = = 0,501∠ − 2, 862
XL

 ∗
S2
I2 = Figura: Esquema del sistema antes
V2 de la falla
= 1,046∠ − 10, 601

It = IL + I2 = 1, 544∠ − 8, 098

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 91 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

La interna del geneador se obtiene

por:
0 0
E = V2 + It (xd + xcc )
= 1, 01 + j0, 2 + 1, 544∠ − 8, 098(j0, 5)
Figura: Esquema del sistema antes de la
= 1, 477∠40, 759 falla

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 92 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

La potencia activa antes de la falla,

viene dada por;
V2 V3
Pe,1 = Sen(θ2 − θ3 )
XL
1, 096
= Sen(11, 2◦ )
0, 4
= 0, 5
Figura: Esquema del sistema antes de la
Pe = Pe,1 + Pe,2 = 1, 5 falla
Como en el estado antes de la falla,
la potencia de aceleración es igual a
cero; la Pm = P e = 1, 5

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 93 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Si se desconecta toda la carga;

E 0 V3
Pedesc = 0 Sen(δ)
xd + xcc + xL
= 1, 641Sen δ

Pe = Pe,1 + Pe,2 = 1, 5

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 94 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Para que el sistema sea estable se debe cumplir que el área de la zona de
aceleración (A1 ) es menor o igual al área de la zona de desaceleración. Es
decir;
A1 = A2
Z δ1 Z δ2
(Pm − 1, 641Sen δ) dδ = (Pm − 1, 641Sen δ) dδ
δ0 δ1

Figura: Curva P − δ luego de la pérdida de carga

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 95 / 89
Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Donde, los lı́mites de integración

viene dados por:

δ0 = 40, 759◦ (0, 7113 rad.)

 
1, 5
δ1 = arcsen
1, 461
= 66, 075◦ (1, 153 rad)
δ2 = π − δ1
Figura: Curva P − δ luego de la pérdida de
= 113, 925◦ (1, 988 rad) carga

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 96 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

El área de aceleración es;

Z δ1
A1 = (Pm − 1, 641Sen δ)
δ0
Z δ1 Z δ1
= Pm dδ − 1, 641 Sen δ dδ
δ0 δ0

= δPm |δδ10 +1, 641Cos δ |δδ10

= 1, 5(1, 153 − 0, 7113) + 1, 641(Cos66, 075◦ − Cos40, 754◦ )
= 0, 085

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 97 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

El área de desaceleración es;

Z δ2
A1 = (1, 641Sen δ − Pm )
δ1
Z δ2 Z δ2
= 1, 641 Sen δ dδ − Pm dδ
δ1 δ1

= −1, 641Cos δ |δδ21 +δPm |δδ21

= 1, 641(−Cos1, 988 + Cos1, 153) − 1, 5(1, 988 − 1, 153)
= 0, 078

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 98 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

A1 = 0,85 > A2 = 0, 078

Como A1 > A2 , el área de

aceleración es mayor al área de
dessacelaración, por lo tanto; el
sistema no es estable. Bajo esta
condición, las protecciones del
Figura: Curva P − δ luego de la pérdida de
generador desconectan la máquina
carga
para evitar la aceleración
incontrolada del mismo

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 99 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Caso 2: Apertura de la lı́nea de transmisión

Figura: Esquema de la red durante la falla

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 100 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Se asume durante la falla, un

modelo de admitancia constante
para la carga. Esto es;

| V 2 |2
Z2 =
S∗
| 1, 03 |2
=
1 − j0, 4
Figura: Esquema de la red durante la falla
= 0, 915 + j0, 366

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 101 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

La corriente por el sistema viene dada por;

| E0 |
| I2 | = 0
| xd + xcc + Z2 |
1, 477
=
| j(0, 35 + 0, 15) + 0, 915 + j0, 366 |
1, 477
=
1, 259
= 1, 173
Figura: Esquema de la red durante
La potencia eléctrica es; la falla

Pe =| I2 |2 < (Z2 )
= 1, 258

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 102 / 89

Ángulo Crı́tico y Tiempo Crı́tico
Ejemplo

Asumiendo que la potencia

mecánica permanece constante.
Entonces;

PA = Pm − Pe
= 1, 5 − 1, 258
= 0, 242

Por lo tanto, el sistema es inestable Figura: Curva de igualdad de áreas

ya que el área de aceleración crece
infinitamente

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 103 / 89

Factores que Afectan la Estabilidad Transitoria

Podemos identificar los siguientes factores que afectan a la estabilidad

transitoria de un generador sńcrono:

La carga del generador: Mientras más carga tenga el generador, menos estable
es el sistema

La distancia a la falla: Cuanto mayor es la distancia, mayor es la potencia

eléctrica inyectada por el generador durante la falla y más estable es el sistema

El tiempo de despeje de la falla: A menor tiempo de despeje mayor estabilidad

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 104 / 89

Factores que Afectan la Estabilidad Transitoria

Podemos identificar los siguientes factores que afectan a la estabilidad

transitoria de un generador sı́ncrono:

La selectividad de las protecciones: Cuanto mayor es, mayor es la sinusoide que

representa la potencia inyectada después de la falla y por lo tanto; más estable
es el sistema

La inercia del generador: Mientras mayor sea su inercia, menos se acelera la

máquina y más estable es el sistema

La tensión de la red de transporte y la tensión interna del generador: Mientras

mayores sean éstas, será mayor la sinusoide que representa la potencia
inyectada en el sistema después de la falla, y por lo tanto; más estable el
sistema

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 105 / 89

Estabilidad Multimáquina

El comportamiento de una máquina sincrona conectada a una barra de potencia

infinita ante una perturbación, es sin duda un mecanismo académico para estudiar
el comportamiento del ángulo de potencia en el tiempo, δ(t).

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 106 / 89

Estabilidad Multimáquina

En este caso el ángulo de potencia, δ(t), es la diferencia entre el ńgulo interno de

la máquina y el ángulo de referencia rotŕica; el cual se supone que es referencia a
la barra de potencia infinita.
dδ(t)
= ω − ω0
dt

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 107 / 89

Estabilidad Multimáquina

Para un sistema de tan solo dos máquinas, se requieren dos ecuaciones de

oscilación; una para cada generador. Donde una de ellas es la referencia.

d2 δ1 (t) ω0
2
= [Pm1 − Pe1 ]
dt 2H1
d2 δ2 (t)2 ω0
2
= [Pm2 − Pe2 ]
dt 2H2

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 108 / 89

Estabilidad Multimáquinas

En cada caso, el ángulo de potencia entre las dos máquinas depende de los
ángulos internos de cada máquina y el águlo de referencia rotatorio

Máquina 1 Máquina 2

dω1 (t) ω0 dω2 (t) ω0

= [Pm1 − Pe1 ] = [Pm2 − Pe2 ]
dt 2H1 dt 2H2
dδ1 (t) dδ2 (t)
= ω1 − ω0 = ω2 − ω0
dt dt

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 109 / 89

Estabilidad Multimáquina

Suponga que se tiene un cierto sistema de potencia con n generadores sı́ncronos,

los cuales inyectan potencias alimentando r número de cargas estáticas a través
de un sistema de transmisión de energı́a.

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 110 / 89

Estabilidad Multimáquina

1 Se toma como referencia la barra 0, el neutro del istema

2 Las barras 1, 2, 3, ..., n; son los nodos internos de las máquinas, es decir; los
buses a los que se les aplica las tensiones detrás de las reactancias transitoria

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 111 / 89

Estabilidad Multimáquina

1 El primer paso para realizar el estudio de estabilidad multimáquina, es

conocer las condiciones del sistemas antes de la perturbación

2 Se realiza un estudio de flujo de carga, para conocer los valores de potencia,

Pk y Qk , y la tensión, magnitud y ángulo, en cada nodo del sistema

3 El resultado del flujo de carga determina las condiciones del sistema (estado
estable) antes que ocurra una perturbción

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 112 / 89

Estabilidad Multimáquina

Modelo de las cargas

1 Las admitancias que representan las cargas se determinan de las condiciones
previas a la perturbación (flujo de carga)

2 Si una determinada barra de carga tiene una tensión VL , y potencia, SL ,

Entonces;  ∗
SL
YL =
VL2

3 Las admitancias en paralelo equivalentes para cada carga del sistema, se

añaden a la correspondiente diagonal de la matriz admitancia del sistema

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 113 / 89

Estabilidad Multimáquina
Modelo de los Generadores
1 El resultado del flujo de carga permite calcular la mgnitud y ángulo de las
tensiones internas de las máquinas previo a la ocurrencia de una perturbación.
Esta magnitud se considera constante a lo largo del estudio de estabilidad.
Pn − Qn
Itn =
V t∗n
0 0
En = V tn − jxdn

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 114 / 89

Estabilidad Multimáquina

Modelo de la Red
1 Las tensiones y corrientes del sistemas están relacionadas por la matriz de
admitancia del sistemas;
[I] = [Y ][V ]

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 115 / 89

Estabilidad Multimáquina

Para incluir la tensión detras de la reactancia transitoria, agregamos m nodos

de generación interna a la matriz de admitancia del sistema de tamaño mxm,
para formar una matriz de tamaño n + m.

La idea es reducir la matriz, de forma tal que solo queden los nodos
correspondientes a los generadores (reducción de Kron)

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 116 / 89

Estabilidad Multimáquina

    " [V ] #
[0] [Ynn ] [Ynm ] h n0 i
= t
[Im ] [Ynm ] [Ymm ] Em
h 0 i
[0] = [Ynn ] [Vn ] + [Ynm ] Em
h 0 i
t
[Im ] = [Ynm ] [Vn ] + [Ymm ] Em

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 117 / 89

Estabilidad Multimáquina

Despejamos la tensiónes en los nodos sin generación,

h 0 i
[0] = [Ynn ] [Vn ] + [Ynm ] Em
h 0 i
−1
[Vn ] = − [Ynn ] [Ynm ] Em

Sustituimos [Vn ] en la ecuación de corrientes;

h 0 i
t
[Im ] = [Ynm ] [Vn ] + [Ymm ] Em
h 0 i h 0 i
−1 t
= − [Ynn ] [Ynm ] Em [Ynm ] + [Ymm ] Em
h ih 0 i
t −1
= [Ymm ] − [Ynm ] [Ynn ] [Ynm ] Em
h 0 i
= [Ygen ] Em

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 118 / 89

Estabilidad Multimáquina

De esta forma, la corriente que entra a la barra de generación k, viene dada

por;

h 0 i
[Im ] = [Ygen ] Em
0 0 0
Ik = yk1 E1 + yk2 E2 + · · · + ykn En
n
X 0
= yk,r Er
r=1

También se sabe que la potencia compleja que inyecta el generador k, viene

dado por;
0
Sk = Pk + Qk = Ek Ik∗

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 119 / 89

Estabilidad Multimáquina

Sustituyendo la corriente en la ecuación de potencia compleja, tenemos;

0
Sk = Ek Ik∗
n
!∗
0 X 0
= Ek yk,r Er
r=1
n
0 X 0
∗
= Ek yk,r Er∗
r=1

Por lo tanto, la potencia activa entregada por el generador k viene dada por;

Pk = Re{Sk }
( n
)
0 X 0
∗
= Re Ek yk,r Er∗
r=1

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 120 / 89

Estabilidad Multimáquina

Si la admitancia y la tensión interna de la máquina se representan en forma

polar, tenemos;

yk,r = | yk,r | ∠θk,r

0 0
Er = | Er | ∠δr

Sustituyen obtenemos;
( n
)
0 X 0
∗
Pk = Re Ek yk,r Er∗
r=1
( n
)
0 X 0
= Re | Ek | ∠δk | yk,r | ∠θk,r | Er | ∠δr
r=1
( n )
X 0 0
= Re | Ek | | yk,r | | Er | ∠(θk,r − δk + δr )
r=1

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 121 / 89

Estabilidad Multimáquina

De esta forma, la potencia inyectada por el generador k viena dada por;

n
X 0 0
Pk = | Ek | | yk,r | | Er | ∠Cos (θk,r − δk + δr )
r=1

Y las respectivas ecuaciones de oscilación;


dωk (t) ω0 h Pn 0 0
i

 = Pmec,k − r=1 | Ek | | yk,r | | Er | ∠Cos (θk,r − δk + δr )
dt 2Hk
ωk (t)
= ωk − ω0


dt

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 122 / 89

Estabilidad Multimáquina

Sistema IEEE-9 Bus

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 123 / 89

Estabilidad Multimáquina

Sistema IEEE-9 bus

©Agustı́n Marulanda, PhD Estabilidad Octubre de 2018 124 / 89