Dinamica e estabilidade de sistemas eletricos de potencia

Prévia do material em texto

Considere o sistema radial representado a seguir, composto por um turbogerador. 
 
Figura 1 - sistema estudado 
Os parâmetros do sistema são mostrados na Tabela 1: 
 
 
Para a condição de normal de operação (antes da aplicação da falta): 
 
O circuito elétrico equivalente: 
 
 
Figura 2 - Circuito elétrico equivalente na condição pré-defeito. 
A potência elétrica fornecida pela maquina e a mesma que chega ate o barramento 
infinito, pois na análise clássica desconsidera-se os elementos resistivos das linhas de 
transmissão 
 Assim para calcula a tensão interna da maquina 𝐸∠𝛿 
 
𝑃 = 
𝑉1 ∙ 𝑉2
𝑋𝑡 +
𝑋𝐿𝑇
2
sin(𝛿1) 
 
0,9 = 
1,05 ∙ 1
0,15 +
0,4
2
∙ sin(𝛿1) 
𝛿1 = sin
−1 (
0,9
3
) = {
17,46°
0,304𝑟𝑎𝑑
 
Desta forma 
𝑉1∠𝛿1 = 1,05∠17,46° 𝑝𝑢 
Como a tensão da barra (2) é conhecida e está na referência, determina-se a corrente 𝐼 
do circuito elétrico da seguinte forma: 
 
𝐼 = 
𝑉1 − 𝑉2
𝑗 ∙ (𝑋𝑡 +
𝑋𝐿𝑇
2 )
 
 
𝐼∠𝜃 = 
1,05∠17,46° − 1∠0°
𝑗 ∙ (0,15 +
0,4
2 )
= 0,9∠ − 0,3° 𝑝𝑢 
 
Finalmente pode-se escrever que: 
 
𝐸∠𝛿 = 𝑋′𝑑 ∙ 𝐼∠𝜃 + 𝑉1∠𝛿1 
 
𝐸∠𝛿 = (0,2∠90° ) ∙ (0,9∠ − 0,3° ) + (1,05∠17,46°) 
𝐸∠𝛿 = (0,18∠89,7° ) + (1,05∠17,46°) 
𝐸∠𝛿0 = 1,118∠26,28° 𝑝𝑢 0.458668𝑟𝑎𝑑 
 
 
Diagrama fasorial do sistema estudado 
 
A equação da potência Elétrica: 
 
𝑃𝑒 =
𝐸 ∙ 𝑉2
𝑋𝑒𝑞
sin(𝛿) 
Onde 
𝑋𝑒𝑞 = 𝑋′𝑑 + 𝑋𝑡 +
𝑋𝐿𝑇
2
= 0,55 𝑝𝑢 
Então 
𝑃𝑒 =
1,118 ∙ 1
0,55
sin(𝛿) = 2,033 sin(𝛿) 𝑝𝑢 
 
A potência máxima transferida: 
 
𝑃𝑚𝑎𝑥 = 2,033 𝑝𝑢 
 
 
 
O coeficiente de potência sincronizante 
 
𝑃𝑠 = 𝑃𝑚𝑎𝑥 cos(𝛿0) = 1.82263pu 
 
Característica da potencia em função do ângulo 
 
 
 
A equação de oscilação da máquina síncrona 
 
𝑑2𝛿
𝑑𝑡2
= 37,7(0,9 − 2,033 sin(𝛿)) 
 
E por fim temos as equações de estado do sistema 
 
𝑑𝜔
𝑑𝑡
= 37.7(0,9 − 2,033 sin(𝛿) 
 
𝑑𝛿
𝑑𝑡
= 𝜔 − 𝜔𝑠 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Parte 2 
 
• Equação de estado pré falta 
𝑃max _𝑝𝑟𝑒 = 2,033 sin(𝛿) 𝑝𝑢 
 
𝑑𝜔
𝑑𝑡
= 37.7(𝑃𝑚 − 2.033 sin(𝛿)) 
 
𝑑𝛿
𝑑𝑡
= 𝜔 − 𝜔𝑠 → 
𝑑𝛿
𝑑𝑡
= 𝜔 − 377 
 
 
 
• Equação de estado durante a falta 
 
𝑑𝜔
𝑑𝑡
= 37.7(𝑃𝑚 − 0) 
 
𝑑𝛿
𝑑𝑡
= 𝜔 − 𝜔𝑠 → 
𝑑𝛿
𝑑𝑡
= 𝜔 − 377 
 
• Equações de estado pós falta 
 
𝐸′ = 1.1181∠26.2766° 𝑝𝑢 
𝑋𝑒𝑞 = 0,15 + 0,4 + 0,2 = 0.75𝑗 
 
𝑃𝑚á𝑥_𝑝𝑜𝑠 =
|𝐸′||𝑉∞|
𝑋𝑒𝑞
=
|1.1181||1|
|0.75|
= 1.4908 𝑝𝑢 
 
𝑃𝑒 = 1.4908 sin 𝛿 
 
𝑑𝜔
𝑑𝑡
= 37.7(𝑃𝑚 − 1.4908 sin(𝛿)) 
 
𝑑𝛿
𝑑𝑡
= 𝜔 − 𝜔𝑠 → 
𝑑𝛿
𝑑𝑡
= 𝜔 − 377 
 
 
 
 
• Giagrama de simulação 
 
 
Determinar o tempo crítico de abertura dos disjuntores da LT2 de maneira 
a manter o sistema estável. Para isto, utilizar o critério das áreas iguais e, 
posteriormente, comprovar os resultados através da simulação no tempo. 
𝐸∠𝛿0𝑝𝑟𝑒 = 1,118∠26,28° 𝑝𝑢 0.458668𝑟𝑎𝑑 
 
 
𝛿𝑐 = cos
−1 [
𝑃𝑚(𝛿𝑚á𝑥 − 𝛿0𝑝𝑟𝑒) + 𝑃𝑝ó𝑠 cos 𝛿𝑚á𝑥
𝑃𝑝ó𝑠 − 0
] 
 
Onde 
𝛿0𝑝𝑜𝑠𝑓𝑎𝑙𝑡𝑎 = sin
−1 (
𝑃𝑚
𝑃𝑚𝑎𝑥𝑝𝑜𝑠
) = 0.648137𝑟𝑎𝑑 
 
𝛿𝑚á𝑥 = 𝜋 − 𝛿0𝑝𝑜𝑠𝑓𝑎𝑙𝑡𝑎 
 
 
𝛿𝑚á𝑥 = 𝜋 − 0,6481 = 2.493455𝑟𝑎𝑑 
 
 
𝛿𝑐 = cos
−1 [
0.9(2.493455 − 0.4586) + 1.4908 cos(2.493455)
 1.4908
] 
= {
1.1249𝑟𝑎𝑑
64.45479°
 
 
 
Tempo critico 
𝑡𝑐 = √
4𝐻(𝛿𝑐 − 𝛿0𝑝𝑟𝑒)
2𝜋𝑓0 ∗ 𝑃𝑚
= √
4𝐻(1.1249 − 0.4586)
2𝜋60 ∗ 0.9
= 0.198185 
 
 
 
 
0 1 2 3 4 5 6 7 8 9 10
0
0.5
1
1.5
TEMPO DE FALTA < Tc ; Dm = 0 ; H = 5s
P
 [
p
u
]
0 1 2 3 4 5 6 7 8 9 10
20
40
60

[r
a
d
]
0 1 2 3 4 5 6 7 8 9 10
374
376
378
380
TEMPO [S]

[r
a
d
/S
]
0 1 2 3 4 5 6 7 8 9 10
-2
0
2
TEMPO DE FALTA > Tc ; Dm = 0 ; H = 5s
P
 [
p
u
]
0 1 2 3 4 5 6 7 8 9 10
0
5
10
x 10
4

[r
a
d
]
0 1 2 3 4 5 6 7 8 9 10
0
500
1000
TEMPO [S]

[r
a
d
/S
]
• Considere agora o efeito dos amortecimentos mecânicos no diagrama de 
simulação. Para D=2, como ficam as equações de estado? O que acontece com 
o tempo crítico de eliminação da falta? 
 
Nota-se que com o amortecimento o tempo critico aumenta. 
 
• Equação de estado pré falta 
𝑃max _𝑝𝑟𝑒 = 2,033 sin(𝛿) 𝑝𝑢 
 
𝑑𝜔
𝑑𝑡
= 37.7(𝑃𝑚 − 2.033 sin(𝛿) − 2(𝜔 − 377)) 
 
𝑑𝛿
𝑑𝑡
= 𝜔 − 𝜔𝑠 → 
𝑑𝛿
𝑑𝑡
= 𝜔 − 377 
 
 
• Equação de estado durante a falta 
 
𝑑𝜔
𝑑𝑡
= 37.7(𝑃𝑚 − 0) 
 
𝑑𝛿
𝑑𝑡
= 𝜔 − 𝜔𝑠 → 
𝑑𝛿
𝑑𝑡
= 𝜔 − 377 
 
• Equações de estado pós falta 
0 1 2 3 4 5 6 7 8 9 10
0
0.5
1
TEMPO DE FALTA > Tc ; Dm = 2 ; H = 5s
P
 [
p
u
]
0 1 2 3 4 5 6 7 8 9 10
25
30
35
40

[r
a
d
]
0 1 2 3 4 5 6 7 8 9 10
377
377.5
378
TEMPO [S]

[r
a
d
/S
]
 
𝐸′ = 1.1181∠26.2766° 𝑝𝑢 
𝑋𝑒𝑞 = 0,15 + 0,4 + 0,2 = 0.75𝑗 
 
𝑃𝑚á𝑥_𝑝𝑜𝑠 =
|𝐸′||𝑉∞|
𝑋𝑒𝑞
=
|1.1181||1|
|0.75|
= 1.4908 𝑝𝑢 
𝑃𝑒 = 1.4908 sin 𝛿 
𝑑𝜔
𝑑𝑡
= 37.7(𝑃𝑚 − 1.4908 sin(𝛿) − 2(𝜔 − 377)) 
𝑑𝛿
𝑑𝑡
= 𝜔 − 𝜔𝑠 → 
𝑑𝛿
𝑑𝑡
= 𝜔 − 377 
 
De acordo com a simulação, é possível avaliar que com o amortecimento o tempo 
critico aumenta. 
 
 
 
 
• Faça variações nos parâmetros/condições do sistema (ao menos algumas). Por 
exemplo: 
➢ O que acontece com a resposta do sistema se o valor de H for alterado? 
H tc 
1 0.089 
5 0.198 
10 0.280 
 
De acordo com os calculo o tempo critico aumenta se o H aumentar 
 0 1 2 3 4 5 6 7 8 9 10
0
0.5
1
1.5
 Tc=0.089 ; Dm = 0 ; H = 1s
P
 [
p
u
]
 
 
Quando H aumenta o sistema a resposta diminui a frequência oscilatória 
 
➢ E se for considerada a construção de uma terceira LT em paralelo com 
mesmo valor de impedância? 
𝑥𝐿𝑇𝑝𝑟𝑒 = 0,133𝑝𝑢 
 
0,9 = 
1,05 ∙ 1
0,15 + 0,133
∙ sin(𝛿1) 
 
𝛿1 = sin
−1 (
0,9
3,71
) = {
14,04°
0,245𝑟𝑎𝑑
 
 
𝐼∠𝜃 = 
1,05∠14,04° − 1∠0°
𝑗 ∙ (0,283)
= 0,9∠ − 4,18° 𝑝𝑢 
 
𝐸∠𝛿 = (0,2∠90° ) ∙ (0,9∠ − 4,18 ) + (1,05∠14,04°) 
𝐸∠𝛿0 = 1,1194∠22,83° 𝑝𝑢 0.398𝑟𝑎𝑑 
Considerando que ocorra alguma falta no início de qualquer uma das linhas 
𝑥𝐿𝑇𝑑𝑢𝑟𝑎𝑛𝑡𝑒 = 𝑖𝑛𝑓 
𝑥𝐿𝑇𝑝𝑜𝑠 = 0,2𝑝𝑢 
 
𝑃𝑒 =
𝐸 ∙ 𝑉2
𝑋𝑒𝑞
sin(𝛿) 
0 1 2 3 4 5 6 7 8 9 10
0
0.5
1
1.5
 Tc = 0.198 ; Dm = 0 ; H = 5s
P
 [
p
u
]
0 1 2 3 4 5 6 7 8 9 10
0
0.5
1
1.5
 Tc = 0.198 ; Dm = 0 ; H = 10s
P
 [
p
u
]
TEMPO [S]
𝑋𝑒𝑞 = 𝑋′𝑑 + 𝑋𝑡 +
𝑋𝐿𝑇
2
= 0,4833 𝑝𝑢 
𝑃𝑒 =
1,1194 ∙ 1
0,4833
sin(𝛿) = 2,316 sin(𝛿) 𝑝𝑢 
A construção de uma terceira linha aumento 12% a potência máxima da linha de 
transmissão 
 
➢ E se for considerada a compensação da LT1? 
O caso de compensação por capacitores irá reduzir a impedância da LT, causando um 
efeito análogo ao da construção de uma nova linha d transmissão 
 
➢ E se fosse permitido variar o valor de E´ através da atuação do regulador de 
tensão? 
O valor de potência transmitida seria maior pois o E’ está diretamente ligado 
com a formula de transferência de potência abaixo, além do que teríamos que 
recalcular nosso tempo e ângulo crítico. Caso eles não sejam recalculados o sistema se 
tornaria instável. 
𝑃𝑚á𝑥 =
|𝐸||𝑉2|
𝑋𝑒𝑞
𝑠𝑒𝑛(𝛿) 
➢ Avalie se o tempo crítico para eliminação da falta trifásica franca forneceria um 
bom indicador quando da ocorrência de uma falta não franca. Esta abordagem 
seria muito conservadora, precitada ou adequada? Pode-se também avaliar o 
valor da resistência de falta nesta afirmação, isto é, verificar se o valor 
considerado influência na conclusão obtida. 
Essa abordagem seria precipitada, pois para uma falta não franca o sistema pode 
manter a transferência de potencia sem perder a estabilidade, logo não existirá um 
tempo critico a ser calculado. 
 
➢ E se o valor da potência mecânica for alterado? 
O valor da potência mecânica interfere diretamente no ângulo da tensão do barramento 
V1, influenciando no ângulo do E’ (delta 0), desta maneira variando delta 0 temos que 
recalcular o tempo critico, ângulo máximo e delta crítico que vão influenciar na 
estabilidade do sistema. Caso eles não sejam recalculados o sistema se tornaria instável. 
 
➢ Mostreque as expressões/conceitos de estabilidade de regime permanente 
podem ser validadas através da simulação no domínio do tempo. 
De acordo com o Saadat 2ªed 
ESTABILIDADE EM ESTADO ESTACIONÁRIO - PEQUENOS DISTÚRBIOS 
O movimento do sistema é livre e a estabilidade é garantida se o sistema retornar ao 
seu estado original. Tal comportamento pode ser determinado em um sistema linear 
examinando a equação característica do sistema. Presume-se que os controles 
automáticos, como regulador de tensão e governador, não estejam ativos. As ações do 
governador e do sistema de saída e dos dispositivos de controle são discutidas no 
capítulo 12 ao lidar com a estabilidade dinâmica. 
Para ilustrar o problema de estabilidade em regime permanente, consideramos o 
comportamento dinâmico de um sistema de uma máquina conectado a um barramento 
infinito, conforme mostrado na figura 11.1. Substituir a energia elétrica de 11.29 na 
equação de balanço dada em 11.21 resulta em 
𝑃𝑒 = 𝑃𝑚𝑎𝑥 sin(𝛿) 11.29 
 
𝐻
𝜋𝑓0
𝑑2𝛿
𝑑𝑡2
= 𝑃𝑚 − 𝑃𝑒 11.21 
 
𝐻
𝜋𝑓0
𝑑2𝛿
𝑑𝑡2
= 𝑃𝑚 − 𝑃𝑚𝑎𝑥 sin(𝛿) 11.36 
A equação de balanço é uma função não linear do ângulo de potência. No entanto, para 
pequenos distúrbios, a equação de oscilação pode ser linearizada com pouca perda de 
precisão como segue. Considere um pequeno desvio no ângulo ∆𝛿 de potência do ponto 
de operação inicial 𝛿0, isto é, 
𝛿 = 𝛿0 + ∆𝛿 
Substituindo na equação 11.36 
𝐻
𝜋𝑓0
𝑑2(𝛿0 + ∆𝛿)
𝑑𝑡2
= 𝑃𝑚 − 𝑃𝑚𝑎𝑥 sin(𝛿0 + ∆𝛿) 
 
𝐻
𝜋𝑓0
𝑑2𝛿0
𝑑𝑡2
+
𝐻
𝜋𝑓0
𝑑2∆𝛿
𝑑𝑡2
= 𝑃𝑚 − 𝑃𝑚𝑎𝑥(sin 𝛿0 cos ∆𝛿 + sin ∆𝛿 cos 𝛿0) 
 
Uma vez que o ângulo ∆𝛿 é pequeno, cos ∆𝛿 ≅ 1 e sin ∆𝛿 ≅ ∆𝛿 
Então temos 
𝐻
𝜋𝑓0
𝑑2𝛿0
𝑑𝑡2
+
𝐻
𝜋𝑓0
𝑑2∆𝛿
𝑑𝑡2
= 𝑃𝑚 − 𝑃𝑚𝑎𝑥 sin 𝛿0 − 𝑃𝑚𝑎𝑥∆𝛿 cos 𝛿0 
Desde o estado operacional inicial 
𝐻
𝜋𝑓0
𝑑2𝛿
𝑑𝑡2
= 𝑃𝑚 − 𝑃𝑚𝑎𝑥 sin(𝛿) 
A equação acima se reduz a equação linearizada em termos de mudanças incrementais 
no ângulo de potência 
𝐻
𝜋𝑓0
𝑑2∆𝛿
𝑑𝑡2
+ 𝑃𝑚𝑎𝑥∆𝛿 cos 𝛿0 = 0 11.38 
 
𝑃𝑚𝑎𝑥 cos 𝛿0 é a inclinação da curva 𝑃 − 𝛿 no ponto 𝛿0 . 
 
É conhecido como coeficiente de sincronização, denotado por 𝑃𝑠. Este coeficiente 
desempenha um papel importante na determinação da estabilidade do sistema e é dado 
por 
𝑃𝑠 = 𝑃𝑚𝑎𝑥 cos 𝛿0 11.39 
 
Substituindo na equação 11.38 
𝐻
𝜋𝑓0
𝑑2∆𝛿
𝑑𝑡2
+ 𝑃𝑠∆𝛿 = 0 
a solução da equação diferencial de segunda ordem acima depende das raízes da 
equação de características dada por 
𝑠2 = −
𝜋𝑓0
𝐻
𝑃𝑠 
quando Ps é negativo, temos uma raiz na metade direita do plano s, a resposta está 
aumentando exponencialmente e a estabilidade é perdida. Quando Ps é positivo, temos 
duas raízes no eixo j-w, e o movimento é oscilatório e sem amortecimento. O sistema é 
marginalmente estável com uma frequência natural de oscilação dada por 
𝜔𝑛 = √
𝜋𝑓0
𝐻
𝑃𝑠 
pode ser visto na figura 11.3 que a faixa onde Ps (ou seja, a inclinação dp / dd) é positiva 
encontra-se entre 0 e 90 com um valor máximo sem carga (ângulo inicial = 0). 
 
Enquanto houver uma diferença na velocidade angular entre o rotor e o campo de 
entreferro giratório resultante, a ação do motor de indução ocorrerá entre eles, e um 
torque será configurado no rotor tendendo a minimizar a diferença entre as duas 
velocidades angulares. Isso é chamado de torque de amortecimento 𝑃𝑑. a potência de 
amortecimento é aproximadamente proporcional ao desvio da velocidade. 
𝑃𝑑 = 𝐷
𝑑𝛿
𝑑𝑡
 
o coeficiente de amortecimento D pode ser determinado a partir de dados de projeto 
ou por teste. Torques de amortecimento adicionais são causados pela característica de 
velocidade / torque do motor principal e a dinâmica da carga, que não são considerados 
aqui. Quando o coeficiente de potência de sincronização 𝑃𝑠 for positivo, por causa da 
potência de amortecimento, as oscilações eventualmente irão amortecer, e a operação 
no ângulo de equilíbrio será restaurada. Nenhuma perda de sincronismo ocorre e o 
sistema é estável. 
Se o amortecimento for levado em consideração, a equação de oscilação linearizada 
torna-se 
𝐻
𝜋𝑓0
𝑑2∆𝛿
𝑑𝑡2
+ 𝐷
𝑑∆𝛿
𝑑𝑡
+ 𝑃𝑠∆𝛿 = 0 
Ou 
𝑑2∆𝛿
𝑑𝑡2
+
𝜋𝑓0
𝐻
𝐷
𝑑∆𝛿
𝑑𝑡
+
𝜋𝑓0
𝐻
𝑃𝑠∆𝛿 = 0 
ou em termos da equação diferencial de segunda ordem padrão, temos: 
𝑑2∆𝛿
𝑑𝑡2
+ 2𝜉𝜔𝑛
𝑑∆𝛿
𝑑𝑡
+ 𝜔𝑛
2∆𝛿 = 0 
𝜔𝑛 é a frequência natural de oscilação e 𝜉 e definido como taxa de amortecimento 
adimensional 
𝜔𝑛 = √
𝜋𝑓0
𝐻
𝑃𝑠 
 
𝜉 = 
𝜋𝑓0𝐷
2𝐻𝜔𝑛
=
𝜋𝑓0𝐷
2𝐻
√
𝐻
𝜋𝑓0𝑃𝑠
=
𝐷
2
√
𝜋𝑓0
𝐻𝑃𝑠
 
A equação característica é 
𝑠2 + 2𝜉𝜔𝑛𝑠+𝜔𝑛
2 = 0 
para condições normais de operação 𝜉 < 1, e as raízes da equação característica são 
complexas 
𝑠1, 𝑠2 = −𝜉𝜔𝑛 ± 𝑗 𝜔𝑛√1 − 𝜉2 
𝑠1, 𝑠2 = −𝜉𝜔𝑛 ± 𝑗 𝜔𝑑 
Onde 𝜔𝑑 é definido como frequência amortecida de oscilação 
𝜔𝑑 = 𝜔𝑛√1 − 𝜉2 
É claro que para amortecimento positivo, as raízes da equação característica têm parte 
real negativa se o coeficiente de potência de sincronização Ps for positivo. A resposta é 
limitada e o sistema é estável. 
 
Agora escrevemos na forma de variável de estado. isso torna possível estender a análise 
para sistemas multimáquinas. 
𝑥1 = ∆𝛿 e 𝑥2 = ∆𝜔 = ∆�̇� então 
𝑥1̇ = 𝑥2 e 𝑥2̇ = −𝜔𝑛²𝑥1 − 2𝜉𝜔𝑛𝑥2 
 
[
𝑥1̇
𝑥2̇
] = [
0 1
−𝜔𝑛² −2𝜉𝜔𝑛
] [
𝑥1
𝑥2
] 
 
𝑥(𝑡)̇ = 𝑨𝑥(𝑡) 
𝑨 = [
0 1
−𝜔𝑛² −2𝜉𝜔𝑛
] 
esta é a equação da variável de estado não forçada ou a equação de estado homogênea. 
se as variáveis de estado x1 e x2 são a resposta desejada, definimos o vetor de saída y 
(t) 
𝑦(𝑡) = [
1 0
0 1
] [
𝑥1
𝑥2
] 
Or 
𝑦(𝑡) = 𝑪𝑥(𝑡) 
Aplicando a transformada de Laplace 
𝑠𝑋(𝑠) − 𝑥(0) = 𝑨𝑋(𝑠) 
𝑋(𝑠) = (𝑠𝑰 − 𝑨)−1𝑥(0) 
Fazendo os cálculos 
𝑋(𝑠) = 
[
𝑠 + 2𝜉𝜔𝑛 1
−𝜔𝑛
2 𝑠
] 𝑥(0)
𝑠2 + 2𝜉𝜔𝑛𝑠+𝜔𝑛2
 
quando o rotor é repentinamente perturbado por um pequeno ângulo ∆𝛿0 
𝑥1(0) = ∆𝛿0 e 𝑥2(0) = ∆𝜔0 = 0 
 
∆𝛿(𝑠) = 
(𝑠 + 2𝜉𝜔𝑛)∆𝛿0
𝑠2 + 2𝜉𝜔𝑛𝑠+𝜔𝑛2
 
E 
∆𝜔(𝑠) = 
𝜔𝑛
2∆𝛿0
𝑠2 + 2𝜉𝜔𝑛𝑠+𝜔𝑛2
 
Fazendo a transformada inversa de Laplace 
∆𝛿 = 
∆𝛿0
√1 − 𝜉2
𝑒−𝜉𝜔𝑛𝑡 ∙ sin(𝜔𝑑𝑡 + 𝜃) 
E 
∆𝜔 = −
𝜔𝑛∆𝛿0
√1 − 𝜉2
𝑒−𝜉𝜔𝑛𝑡 ∙ sin(𝜔𝑑𝑡) 
E 𝜃 e dado por 
𝜃 = cos−1(𝜉) 
o movimento do rotor em relação ao campo giratório síncrono é 
𝛿 = 𝛿0 + 
∆𝛿0
√1 − 𝜉2
𝑒−𝜉𝜔𝑛𝑡 ∙ sin(𝜔𝑑𝑡 + 𝜃) 
e a frequência angular do rotor é 
𝜔 = 𝜔0 −
𝜔𝑛∆𝛿0
√1 − 𝜉2
𝑒−𝜉𝜔𝑛𝑡 ∙ sin(𝜔𝑑𝑡) 
a constante de tempo de resposta é 
𝜏 =
1
𝜉𝜔𝑛
=
2𝐻
𝜋𝑓0𝐷
 
e a resposta se estabelece em aproximadamente quatro constantes de tempo, e a 
resolução 
𝑡𝑠 = 4𝜏 
de 11.42 e 11.47, notamos que conforme a constante de inércia H aumenta, a frequência 
natural e a taxa de amortecimento diminuem, resultando em um longo tempo de 
acomodação. Um aumento no coeficiente de potência de sincronização Ps resulta em 
um aumento na frequência natural e uma diminuição na taxa de amortecimento 
➢ Faça simulações e verifique o aconteceria com tempo crítico de eliminação da 
falta se este mesmo curto-circuito fosse eliminado sem a retirada da linha de 
transmissão (em vez da abertura da LT como ocorre no enunciado). 
Utilizando o programa eacfault.m (anexo no fim do relatório) podemos analisar qual 
será o tempo critico caso a linha de transmissão não seja retirada do SEP 
>> Eacfault 
Generator output power in p.u. Pm = 0.9 
Generator e.m.f. in p.u. E = 1.1181 
Infinite bus-bar voltage in p.u. V = 1 
Reactance before Fault in p.u. X1 = 0.55 
Reactance during Fault in p.u. X2 = inf 
Reactance aftere Fault in p.u. X3 = 0.55 
For this case tc can be found from analytical formula. 
To find tc enter Inertia Constant H, (or 0 to skip) H = 5 
Initial power angle = 26.277 
Maximum angle swing = 153.723 
Critical clearing angle = 84.946 
Critical clearing time= 0.246 sec. 
Logo, o tempo critico seria maior para uma falta onde a LT afetada voltasse após o 
período de falta 
 
ANEXO 
 
function eacfault(Pm, E, V, X1, X2, X3) 
% This program obtains the power angle curves for a one-
machine system 
% before fault, during fault and after the fault clearance. 
% The equal area criterion is applied to find the critical 
clearing angle 
% for the machine to stay synchronized to the infinite bus 
bar 
% 
% Copyright (c) 1998 by H. Saadat 
 
if exist('Pm')~=1 
Pm = input('Generator output power in p.u. Pm = '); else, 
end 
if exist('E')~=1 
E = input('Generator e.m.f. in p.u. E = '); else, end 
if exist('V')~=1 
V = input('Infinite bus-bar voltage in p.u. V = '); else, 
end 
if exist('X1')~=1 
X1 = input('Reactance before Fault in p.u. X1 = '); else, 
end 
if exist('X2')~=1 
X2 = input('Reactance during Fault in p.u. X2 = '); else, 
end 
if exist('X3')~=1 
X3 = input('Reactance aftere Fault in p.u. X3 = '); else, 
end 
 
Pe1max = E*V/X1; Pe2max=E*V/X2; Pe3max=E*V/X3; 
delta = 0:.01:pi; 
Pe1 = Pe1max*sin(delta); Pe2 = Pe2max*sin(delta); Pe3 = 
Pe3max*sin(delta); 
d0 =asin(Pm/Pe1max); dmax = pi-asin(Pm/Pe3max); 
cosdc = (Pm*(dmax-d0)+Pe3max*cos(dmax)-
Pe2max*cos(d0))/(Pe3max-Pe2max); 
 if abs(cosdc) > 1 
 fprintf('No critical clearing angle could be found.\n') 
 fprintf('system can remain stable during this 
disturbance.\n\n') 
 return 
 else, end 
dc=acos(cosdc); 
if dc > dmax 
fprintf('No critical clearing angle could be found.\n') 
 fprintf('System can remain stable during this 
disturbance.\n\n') 
 return 
 else, end 
Pmx=[0 pi-d0]*180/pi; Pmy=[Pm Pm]; 
x0=[d0 d0]*180/pi; y0=[0 Pm]; xc=[dc dc]*180/pi; yc=[0 
Pe3max*sin(dc)]; 
xm=[dmax dmax]*180/pi; ym=[0 Pe3max*sin(dmax)]; 
d0=d0*180/pi; dmax=dmax*180/pi; dc=dc*180/pi; 
x=(d0:.1:dc); 
y=Pe2max*sin(x*pi/180); 
y1=Pe2max*sin(d0*pi/180); 
y2=Pe2max*sin(dc*pi/180); 
x=[d0 x dc]; 
y=[Pm y Pm]; 
xx=dc:.1:dmax; 
h=Pe3max*sin(xx*pi/180); 
xx=[dc xx dmax]; 
hh=[Pm h Pm]; 
delta=delta*180/pi; 
if X2 == inf 
fprintf('\nFor this case tc can be found from analytical 
formula. \n') 
H=input('To find tc enter Inertia Constant H, (or 0 to 
skip) H = '); 
 if H ~= 0 
 d0r=d0*pi/180; dcr=dc*pi/180; 
 tc = sqrt(2*H*(dcr-d0r)/(pi*60*Pm)); 
 else, end 
else, end 
%clc 
fprintf('\nInitial power angle = %7.3f \n', d0) 
fprintf('Maximum angle swing = %7.3f \n', dmax) 
fprintf('Critical clearing angle = %7.3f \n\n', dc) 
if X2==inf & H~=0 
fprintf('Critical clearing time = %7.3f sec. \n\n', tc) 
else, end 
 
h = figure; figure(h); 
fill(x,y,'m') 
hold; 
fill(xx,hh,'c') 
plot(delta, Pe1,'-', delta, Pe2,'r-', delta, Pe3,'g-', 
Pmx, Pmy,'b-', x0,y0, xc,yc, xm,ym), grid 
 
title('Application of equal area criterion to a critically 
cleared system') 
xlabel('Power angle, degree'), ylabel(' Power, per unit') 
text(5, 1.07*Pm, 'Pm') 
text(50, 1.05*Pe1max,['Critical clearing angle = 
',num2str(dc)]) 
axis([0 180 0 1.1*Pe1max]) 
hold off;