Estabilidade Transitória em Sistemas Elétricos de Potência

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

UNIVERSIDADE ESTADUAL DO OESTE DO PARANÁ 
CENTRO DE ENGENHARIAS E CIÊNCIAS EXATAS 
CAMPUS DE FOZ DO IGUAÇU 
CURSO DE ENGENHARIA ELÉTRICA 
 
 
 
TRABALHO DE CONCLUSÃO DE CURSO 
 
 
 
 
ESTUDO DE ESTABILIDADE TRANSITÓRIA EM SISTEMAS 
ELÉTRICOS DE POTÊNCIA 
 
 
 
GUILHERME AUGUSTO MORELLO CAGNINI 
 
 
 
 
FOZ DO IGUAÇU - PR 
2018 
 
 
 
 
GUILHERME AUGUSTO MORELLO CAGNINI 
 
 
 
 
 
 
 
 
 
ESTUDO DE ESTABILIDADE TRANSITÓRIA EM SISTEMAS 
ELÉTRICOS DE POTÊNCIA 
 
 
 
 
Trabalho de Conclusão de Curso apresentado 
ao Curso de Engenharia Elétrica da 
Universidade Estadual do Oeste do Paraná, 
como parte dos requisitos para a obtenção do 
título de Engenheiro Eletricista. 
 
Orientador: Prof. M. Sc. Robson Almir de 
Oliveira 
 
 
 
 
 
 
 
FOZ DO IGUAÇU - PR 
2018 
 
 
 
 
 
 
 
 
RESUMO 
CAGNINI, G. A. M. (2018). ESTUDO DE ESTABILIDADE TRANSITÓRIA EM 
SISTEMAS ELÉTRICOS DE POTÊNCIA. Monografia de Trabalho de Conclusão de Curso 
(Graduação) – Curso de Engenharia Elétrica, Universidade Estadual do Oeste do Paraná – 
UNIOESTE, Foz do Iguaçu, 2018. 
 
Neste trabalho serão discutidos inicialmente, o problema de estabilidade e suas classificações. 
Como o tema é amplo e complexo, o problema será limitado ao estudo de estabilidade 
transitória. Desta forma serão apresentados os modelos e métodos matemáticos para análise de 
estabilidade transitória em sistemas elétricos de potência. Uma breve comparação entre esses 
métodos será realizada a fim de determinar qual o método mais adequado na resolução do 
problema de estabilidade. 
Com a consolidação do conhecimento do problema, dos modelos e dos métodos numéricos será 
formulado um fluxograma expondo o funcionamento do algoritmo a ser implementado. Este 
algoritmo será desenvolvido a partir da plataforma MATLAB e tem por objetivo determinar a 
estabilidade transitória em sistemas elétricos de potência (SEP). 
A partir do programa elaborado, será estudada a estabilidade transitória com base em 
simulações de eventos para diferentes sistemas e cenários. Essas análises permitirão determinar 
se o sistema é estável ou instável e também determinar o tempo crítico de abertura, caso exista. 
 
Palavras-chave: Estabilidade transitória; Métodos numéricos; Tempo crítico de abertura. 
 
 
 
ABSTRACT 
CAGNINI, G. A. M. (2018). STUDY OF TRANSIT STABILITY IN ELECTRICAL 
POWER SYSTEMS. End of Course Paper (Graduation) – Course of Electrical Engineering, 
State University of Western Parana – UNIOESTE, Foz do Iguaçu, 2018. 
 
In this work the stability problem and its classifications will be discussed initially. As the theme 
is broad and complex, the problem will be limited to the study of transient stability. In this way 
will be presented the mathematical models and methods for analysis of transient stability in 
electrical power systems. A brief comparison between these methods will be performed to 
determine the most appropriate method for solving the stability problem. 
With the consolidation of the knowledge of the problem, the models and the numerical methods 
will be formulated a flow chart exposing the operation of the algorithm to be implemented. This 
algorithm will be developed from the MATLAB platform and aims to determine the transient 
stability in electrical power systems (EPS). 
From the developed program, the transient stability will be studied based on event simulations 
for different systems and scenarios. These analyzes will determine whether the system is stable 
or unstable and also determine the critical opening time, if any. 
 
Keywords: Transient stability; Numerical methods; Critical opening time. 
 
 
 
Lista de Abreviaturas e Siglas 
CAI Critério das Áreas Iguais 
E’ Força eletromotriz da máquina 
EDO Equação diferencial ordinária 
ESP Electrical Power System 
H Relação entre a energia cinética armazenada em megajoule na velocidade 
síncrona pela potência nominal da máquina em MVA 
I Corrente elétrica fornecida pela máquina 
J Momento de inércia total das massas do rotor, em kg.m2 
M Constante de inércia 
p. u Por unidade 
Pa Potência de aceleração 
Pe Potência elétrica no entreferro 
Pm Potência de entrada do eixo da máquina menos perdas rotacionais 
SEP Sistema elétrico de potência 
T Tempo em segundos 
Ta Torque de aceleração resultante, em N.m 
TCA Tempo crítico de abertura 
Te Torque elétrico ou eletromagnético, em N.m 
Tm Torque do eixo ou torque mecânico, em N.m 
 𝑉𝑡 
 
Tensão nos terminais da máquina 
𝑋𝑑
′ Reatância transitória de eixo direto 
Z Impedância característica da linha 
γ Constante de propagação 
δm Deslocamento angular do rotor, em radianos mecânicos. 
 
 
θm Medida absoluta do ângulo do rotor que cresce com o tempo, em radianos 
mecânicos 
ωsm Velocidade síncrona da máquina 
 
 
 
Lista de Figuras 
Figura 2.1 – Classificação do problema de estabilidade. ................................................ 17 
Figura 2.2 – Sistema máquina-barra infinita. .................................................................. 22 
Figura 2.3 – Características da curva ângulo-potência ................................................... 23 
Figura 2.4 – Representação do modelo clássico do gerador síncrono ............................ 24 
Figura 2.5 – Modelo da carga para o fluxo de potência. ................................................. 26 
Figura 2.6 – Modelo da carga de impedância constante. ................................................ 26 
Figura 2.7 – Modelo de linha curta. ................................................................................ 27 
Figura 2.8 – Modelo de linha média. .............................................................................. 28 
Figura 2.9 – Modelo de linha longa. ............................................................................... 28 
Figura 2.10 – Solução do sistema particionada. .............................................................. 30 
Figura 2.11 – Representação gráfica do método de Euler .............................................. 33 
Figura 2.12 – Curva das equações preditora e corretora ................................................. 35 
Figura 2.13 - Inclinações usadas pelo método de Runge-Kutta de quarta ordem ........... 38 
Figura 2.14 – Exemplo de matriz jacobiana.................................................................... 41 
Figura 3.1 – Representação do SEP para análise de estabilidade transitória. ................. 43 
Figura 3.2 – Sistema em condições normais de operação ............................................... 45 
Figura 3.3 – Sistema de potência com falta no meio da linha......................................... 45 
Figura 3.4 – Abertura dos disjuntores e remoção da linha .............................................. 46 
Figura 3.5 – Fluxograma do algoritmo de estabilidade transitória. ................................ 49 
Figura 4.1 – Sistema teste seis barras, duas máquinas e barra infinita ........................... 51 
Figura 4.2 – Ângulos dos rotores em valores absolutos para falta na barra 1 sem o 
desligamento de linha. .............................................................................................................. 53 
Figura 4.3 – Ângulos dos rotores em valores absolutos para falta na barra 2 com o 
desligamento da linha 2 – 4. ..................................................................................................... 54 
Figura 4.4 – Ângulos dos rotores em valores absolutos e diferença entre os ângulos para 
falta no meio da linha 1 – 6 com o desligamento da linha 1 – 6. Tempo de abertura 1.1s. ...... 55 
Figura 4.5 – Ângulos dos rotores em valores absolutos e diferença entre os ângulos para 
falta no meio da linha 1 – 6 com o desligamento da linha 1 – 6. Tempo de abertura 1.2s. ...... 56 
Figura 4.6 – Ângulos dos rotores em valores absolutos e diferença entre os ângulos para 
falta
no meio da linha 1 – 6 sem o desligamento da linha 1 – 6. Tempo de abertura 1.1s. ...... 57 
Figura 4.7 – Ângulos dos rotores em valores absolutos e diferença entre os ângulos para 
falta no meio da linha 1 – 6 sem o desligamento da linha 1 – 6. Tempo de abertura 1.2s. ...... 57 
Figura 4.8 – Sistema teste sete barras, três máquinas e barra infinita ............................. 59 
 
 
Lista de Tabelas 
 
Tabela 4.1 – Dados de carga ........................................................................................... 51 
Tabela 4.2 – Dados de linha incluindo os transformadores ............................................ 51 
Tabela 4.3 – Dados de gerador ........................................................................................ 52 
Tabela 4.4 – Resultados do fluxo de carga ..................................................................... 52 
Tabela 4.5 – Tempos de abertura para diferentes cenários ............................................. 58 
Tabela 4.6 – Dados de carga ........................................................................................... 60 
Tabela 4.7 – Dados de linha incluindo os transformadores ............................................ 60 
Tabela 4.8 – Dados de gerador ........................................................................................ 60 
Tabela 4.9 – Resultados do fluxo de carga ..................................................................... 61 
Tabela 4.10 – Tempos de abertura para diferentes cenários ........................................... 61 
 
 
 
SUMÁRIO 
Capitulo 1 – Introdução ............................................................................................................ 12 
1.1 Objetivos ............................................................................................................... 13 
1.1.1 Objetivo geral ................................................................................................. 13 
1.1.2 Objetivos específicos ..................................................................................... 13 
1.2 Justificativa da necessidade de estudos de estabilidade transitória ....................... 14 
1.3 Estrutura do trabalho ............................................................................................. 15 
Capítulo 2 – Fundamentação Teórica ....................................................................................... 16 
2.1 Introdução ............................................................................................................. 16 
2.2 O problema de estabilidade em sistemas elétricos de potência ............................. 16 
2.2.1 Classificação da Estabilidade ......................................................................... 17 
2.2.2 Estabilidade Angular ...................................................................................... 19 
2.3 Modelagem matemática dos componentes necessários ao estudo de estabilidade 
transitória e fluxo de carga ................................................................................................... 24 
2.3.1 Máquina síncrona ........................................................................................... 24 
2.3.2 Carga .............................................................................................................. 25 
2.3.3 Linhas de Transmissão ................................................................................... 27 
2.3.4 Transformadores ............................................................................................ 29 
2.4 Modelagem da rede reduzida ................................................................................ 29 
2.5 Métodos numéricos de soluções de equações diferenciais ................................... 31 
2.5.1 Método de Euler ............................................................................................. 32 
2.5.2 Método de Heun ............................................................................................. 33 
2.5.3 Método de Runge-Kutta de quarta ordem ...................................................... 35 
2.5.4 Comparação dos métodos .............................................................................. 39 
2.6 Fluxo de Potência .................................................................................................. 39 
Capítulo 3 – Estruturação do problema e sistemática do algoritmo ......................................... 42 
3.1 Resolução do Fluxo de potência ........................................................................... 42 
3.2 Determinação da Estabilidade do SEP .................................................................. 46 
Capítulo 4 – Estudo de caso ..................................................................................................... 50 
4.1 Introdução ............................................................................................................. 50 
4.2 Configuração do sistema teste ............................................................................... 50 
4.3 – Eventos simulados .............................................................................................. 52 
 
 
4.3.1 Curto-circuito na barra 1 sem abertura de linha ................................................. 53 
4.3.2 Curto-circuito na barra 2 com o desligamento da linha 2 – 4 ............................ 54 
4.3.3 Curto-circuito no meio da linha 1 – 6 com desligamento da linha 1 – 6 ........... 55 
4.3.4 Curto-circuito no meio da linha 1 – 6 sem o desligamento da linha 1 – 6 ......... 56 
4.3.5 Resumo das simulações realizadas ..................................................................... 58 
4.4 Consequências da adição de um gerador ao sistema ........................................... 59 
Capítulo 5 – Conclusões ........................................................................................................... 63 
5.1 Trabalhos futuros .................................................................................................. 64 
Referências ............................................................................................................................... 65 
Anexos ...................................................................................................................................... 69 
 
 
 
12 
 
Capitulo 1 – Introdução 
Em Sistemas Elétricos de Potência (SEP) de grande porte, como o sistema elétrico 
brasileiro, a operação é efetuada para que o atendimento às cargas ocorra de forma econômica 
e segura, respeitando as diversas restrições elétricas que possam existir. As restrições referem-
se a limites físicos em equipamentos e limites de segurança dinâmicos para intercâmbio entre 
regiões (SILVA, 2009). 
Para o atendimento a essas restrições, é necessário que o sistema elétrico seja 
constantemente avaliado por meio de resultados de simulações, uma vez que está 
constantemente sujeito às ocorrências que causam perturbações no seu estado normal. Estas 
perturbações alteram as grandezas elétricas (corrente, tensão e frequência), provocando 
violações nas restrições operativas. As perturbações mais comuns e também as mais severas 
são os curtos-circuitos (geralmente chamados de defeitos ou faltas), que ocorrem em 
decorrência da ruptura da isolação entre as fases ou entre a fase e terra (FREITAS, DIAS, et 
al., 2012). 
A propriedade do sistema que permite as máquinas síncronas do sistema responder a um 
distúrbio (falta), a partir de uma condição normal de operação, de tal maneira a retomarem uma 
condição de operação novamente normal é caracterizada como estabilidade de sistema de 
potência. Estudos de estabilidade são usualmente classificados em três tipos, dependendo da 
natureza e ordem de grandeza do distúrbio. Estes estudos são chamados estudo de estabilidade 
transitória, dinâmica e em regime permanente (COSTA, 2010).
Conforme será abordado no capítulo 2, a análise de estabilidade do sistema elétrico de 
potência constitui-se de um estudo muito amplo e complexo, portanto este trabalho procurou se 
limitar ao estudo da estabilidade do sistema elétrico em relação à variação do ângulo do rotor 
dos geradores síncronos a grandes perturbações, também denominada análise de estabilidade 
transitória. 
O estudo de estabilidade transitória consiste nos casos envolvendo variações grandes e 
bruscas (impactos) de gerações e/ou cargas, as quais podem provocar perdas de sincronismo 
entre as máquinas síncronas ligadas no sistema (LOPES e SILVA, 2012). 
 
 
13 
 
 
1.1 Objetivos 
Dentro do contexto apresentado na subseção anterior, são apresentados nesta subseção, 
os objetivos gerais e específicos deste trabalho. 
 
1.1.1 Objetivo geral 
Este trabalho tem por objetivo verificar o comportamento do ângulo do rotor da máquina 
síncrona, de maneira que após uma contingência do sistema de potência possa se determinar a 
estabilidade ou instabilidade dele. Para isso, será estudado o problema de estabilidade 
transitória de ângulo do rotor e desenvolvido um algoritmo na plataforma MATLAB apto a 
realizar esta análise. Através dos resultados será determinado o tempo crítico de abertura, tempo 
máximo para que a proteção (disjuntores das linhas) atue sem levar as máquinas conectadas ao 
SEP a instabilidade. 
 
1.1.2 Objetivos específicos 
Para a análise de estabilidade de ângulo do rotor será necessária a resolução da equação 
swing das máquinas síncronas, a qual define o comportamento do ângulo do rotor (SILVA, 
2009). Para tal, apresenta-se neste trabalho três métodos numéricos. Baseando-se na literatura 
será efetuada uma breve comparação entre esses métodos e determinado qual o método mais 
adequado. Sendo os métodos estudados: Euler, Heun e Runge-Kutta de quarta ordem. 
Além disso, independentemente do método adotado para a resolução da equação da barra 
swing, algumas informações para seu desenvolvimento serão necessárias, como por exemplo, 
o conhecimento da força eletromotriz (f.e.m) do gerador e a potência ativa fornecida por ele. 
Para determinar estes parâmetros, a partir da plataforma MATLAB será desenvolvido um 
algoritmo para a obtenção dos dados de fluxo de carga. Sendo que o objetivo do fluxo de 
potência é encontrar as magnitudes e fases das tensões de cada barra em um sistema a partir dos 
dados de cargas ativas e reativas, da tensão dos geradores e das potências ativas dos geradores 
(MATOS, 2012), através de equacionamentos matemáticos, como será demonstrado no capítulo 
14 
 
3, se torna possível determinar as f.e.ms dos geradores. Já a potência mecânica é diretamente 
determinada pelo fluxo de carga, pois é igual a potência ativa fornecida pelos geradores. 
Por fim, será escolhido um sistema no qual serão realizadas simulações e análises para 
validar o método designado para a solução da equação da barra swing e do algoritmo de fluxo 
de carga, assim como, proposta uma solução para o problema. 
 
1.2 Justificativa da necessidade de estudos de estabilidade 
transitória 
Conforme DIAS E PILONO (2010), segundo BRETAS e ALBERTO (2000), um sistema 
elétrico de potência é formado basicamente por geradores síncronos e cargas interligados por 
linhas de transmissão. O fluxo de potência ativa nas linhas está relacionado às diferenças entre 
os ângulos de fase dos geradores. Em condições normais de operação os alternadores estão na 
velocidade síncrona e existe um balanceamento entre a potência elétrica e a potência mecânica 
proveniente das fontes primárias. Ou seja, o sistema opera em regime permanente, mantendo 
assim as diferenças entre os ângulos de fases constantes para que o fluxo de potência também 
seja constante. 
 Entretanto, ainda conforme DIAS E PILONO (2010), em situações decorrentes de uma 
grande perturbação, como: curto-circuito, abertura de uma linha de transmissão, perda de 
geração, ou uma combinação desses eventos, o equilíbrio entre a potência elétrica gerada e a 
potência mecânica de entrada é drasticamente rompido, principalmente nos geradores mais 
próximos ao defeito. Em consequência os rotores das máquinas sofrem diferentes acelerações, 
levando eventualmente algumas delas a perda do sincronismo com o resto do sistema. 
De acordo com PINHEIRO, SILVA, et al (2017), durante uma condição de perda de 
sincronismo há uma grande variação na corrente e na tensão com relação à frequência da 
máquina afetada. A amplitude de corrente e a frequência de operação acima do valor nominal 
podem resultar em stress dos enrolamentos e torques pulsantes, consequentemente provocando 
vibrações mecânicas que podem ser prejudiciais ao gerador. 
Já as consequências ao SEP, segundo PINHEIRO, SILVA, et al (2017), estão ligadas com 
a sua perda de estabilidade pela perda de sincronismo do gerador com o sistema. Quando o SEP 
15 
 
não volta ao seu ponto de equilíbrio é dito que ele está instável e esse estado do sistema precisa 
ser mitigado rapidamente para evitar possíveis faltas de energia e danos aos seus equipamentos. 
 
Logo, o estudo de estabilidade transitória é de grande importância prática, pois aborda 
fenômenos de grande impacto para o sistema elétrico. Possibilitando a previsão do seu 
comportamento, após a condição atípica ao qual ele foi submetido. Podendo dessa forma, 
antecipar as devidas mudanças que devem ser efetuadas na rede para que ocorrida uma falta o 
sistema seja capaz de se recuperar e alcançar uma condição estável de equilíbrio (SILVA, 
2017). 
1.3 Estrutura do trabalho 
No capítulo 2 é apresentada a classificação do problema de estabilidade, a modelagem 
matemática dos componentes do sistema elétrico de potência, o cálculo do fluxo de carga, teoria 
fundamental para as condições iniciais do SEP e o comportamento dinâmico da máquina 
síncrona e sua representação nos estudos de estabilidade transitória, com atenção para o ângulo 
do rotor. Por fim, os métodos numéricos para a resolução das equações diferenciais no domínio 
do tempo provenientes do problema de estabilidade transitória. 
No capítulo 3 são expostas a sistemática e modelagem envolvidas no desenvolvimento 
do algoritmo de análise de estabilidade transitória. 
No capítulo 4 são apresentados os resultados e simulações utilizando os modelos e 
metodologias apresentados nos capítulos 2 e 3. 
Por último, no capítulo 5 são apresentadas as conclusões obtidas ao longo deste trabalho, 
assim como sugestões para trabalhos futuros. 
 
 
16 
 
Capítulo 2 – Fundamentação Teórica 
2.1 Introdução 
De acordo com STOTT (1979), a análise de estabilidade é uma das três mais comuns 
efetuadas em sistemas elétricos de potência (SEP), sendo as outras duas, a do fluxo de potência 
e de curto-circuito. Todavia, a análise de estabilidade, em termos de modelagem e de processo 
de solução, é a mais complexa. Por conseguinte, neste capítulo, foram estudados os tipos de 
estabilidade existentes nos SEPs, dando destaque à estabilidade transitória. 
Para a obtenção de variáveis de interesse na solução do problema de estabilidade, como 
o módulo e ângulo da força eletromotriz dos geradores, assim como a potência mecânica 
fornecida por eles, também é realizada uma breve revisão sobre o fluxo de potência e o método 
para resolve-lo. 
Além disso, devido à grande variedade de elementos existentes nos sistemas elétricos de 
potência, também foi realizado um estudo mais detalhado nas formas de modelagem e 
princípios de funcionamento dos principais elementos constituintes do sistema, dados por: 
geradores, transformadores de potência, linhas de transmissão e cargas (LUZ, 2015). 
 
2.2 O problema de estabilidade em sistemas elétricos de potência 
Estabilidade de um sistema elétrico de potência pode ser definida como uma propriedade 
que o mesmo apresenta em se manter operacionalmente em equilíbrio sob condições normais e 
após perturbações.
A estabilidade depende das condições iniciais do sistema e da severidade da 
perturbação. De acordo com IEEE/CIGRE (2004): a estabilidade de um sistema elétrico de 
potência é a habilidade do sistema, para uma determinada condição de operação, de atingir um 
estado de equilíbrio após estar sujeito a alguma perturbação física com a maioria de suas 
variáveis limitadas, isto é, dentro de limites operacionais, de tal forma que praticamente todo o 
sistema permaneça intacto (SOHN, 2014). 
Quando o sistema não retorna ao seu ponto de equilíbrio na ocorrência de um distúrbio, 
ele é dito instável. Instabilidades devem ser mitigadas durante a operação, evitando, assim, 
períodos de falta de energia e, também, danos aos equipamentos do Sistema Elétrico 
17 
 
 
(BORDEIRA, 2011). Exemplos críticos de falta de energia ocasionados por problemas de 
instabilidade são os blecautes dos EUA e Europa (Inglaterra, Países Nórdicos e Itália) em 2003, 
citados em (BARBOSA, 2013). 
 
2.2.1 Classificação da Estabilidade 
O estudo de estabilidade de um SEP se divide em três classes: Estabilidade de ângulo de 
rotor, Estabilidade de frequência e Estabilidade de tensão. Segundo BARBOSA (2013), esta 
classificação tem por base os seguintes fatores: o fenômeno que caracteriza o tipo de 
instabilidade e as causas físicas que conduzem à sua ocorrência. 
A Figura 2.1 indica que além dessas três classes o problema de estabilidade de sistemas 
elétricos de potência pode ser dividido em diferentes subproblemas. 
 
Fonte: (SOHN, 2014). 
Figura 2.1 – Classificação do problema de estabilidade. 
 
As diferentes divisões de estabilidade mostradas na Figura 2.1, são classificadas segundo 
três vertentes: 
18 
 
 
1 – Tempo {
Curto prazo
Médio prazo
Longo prazo
 
2 – Variáveis de interesse {
Ângulo
Tensão
Frequência
 
3 – Tipos de perturbação {
Grandes Perturbações
Pequenas Perturbações
Variações lentas e previsíveis de carga
 
A estabilidade de frequência refere-se à capacidade de um SEP manter a frequência após 
a ocorrência de um incidente severo, resultando no desequilíbrio entre a produção e a carga e 
dependerá da capacidade do sistema para manter/restaurar o equilíbrio entre a produção e a 
carga, com o mínimo de desligamento de cargas (BARBOSA, 2013). 
Já a estabilidade de tensão consiste da capacidade do SEP de manter as tensões nas barras 
em torno de seus valores nominais durante operação em regime permanente e também durante 
ocorrência de distúrbios. Para a caracterização de uma condição de instabilidade de tensão, o 
distúrbio deve ser capaz de reduzir progressivamente a magnitude da tensão em pelo menos 
uma das barras do sistema, de forma que não seja possível controlar esse parâmetro 
(BORDEIRA, 2011). 
Enfim, a estabilidade angular, também conhecida como estabilidade rotórica, é a 
capacidade da máquina síncrona de um SEP interligado se manter em sincronismo após a 
ocorrência de uma perturbação. Conforme BORDEIRA (2011), para a caracterização de uma 
condição de instabilidade angular, o ângulo do rotor de pelo menos uma das máquinas do 
sistema aumenta indefinidamente em relação aos demais sem que seja possível controlá-lo. 
Como exposto, o problema de análise de estabilidade de um SEP é bem abrangente, pois 
engloba a análise de características distintas em diferentes amplitudes e tempos. 
Por conseguinte, como o objetivo deste trabalho consiste em verificar o comportamento 
do ângulo do rotor da máquina síncrona, esta monografia se limita a análise de estabilidade 
transitória referente a grandes perturbações. 
19 
 
 
2.2.2 Estabilidade Angular 
Há dois tipos de estabilidade angular, sendo eles: “Estabilidade angular para Pequenas 
perturbações” e “Estabilidade angular para Grandes perturbações”, ou “Transitória”. 
O primeiro, de acordo com BORDEIRA (2011), se refere ao estudo de pequenas 
alterações (como a entrada de uma pequena carga) nos parâmetros do sistema durante a 
operação em regime permanente. Para ser enquadrada nesse conceito, a análise pode ser feita 
através da linearização das equações. 
O segundo tipo, ainda conforme BORDEIRA (2011), refere-se à capacidade dos 
geradores síncronos de permanecerem em sincronismo, após a ocorrência de grandes 
perturbações. Ao analisar a estabilidade a grandes perturbações deseja-se investigar a 
capacidade do sistema elétrico em absorver grandes impactos causados por modificações 
estruturais sensíveis, como a ocorrência de curto-circuito, chaveamento de linhas ou entrada de 
um grande bloco de carga (WALANTUS, 2014). Nesses casos, no entanto, a linearização das 
equações não permite realizar a análise da operação, já que, as não-linearidades inerentes ao 
SEP não podem ser desprezadas. A resposta resultante do sistema envolve grandes excursões 
dos ângulos dos geradores e é influenciada pela relação não linear entre potência e ângulo 
(KUNDUR, 1994). 
Para determinar se um sistema é estável transitoriamente, é necessário que este esteja 
convenientemente descrito através de equações matemáticas. As equações diferenciais que 
descrevem o comportamento dinâmico do sistema podem ser obtidas através de um balanço de 
potência em cada máquina do sistema. A máquina é movimentada por um elemento primário 
que lhe fornece potência mecânica. Parte dessa energia mecânica é convertida em energia 
elétrica, que por sua vez é entregue a rede. A parte que não é convertida em potência elétrica 
transforma-se em potência de aceleração do rotor da máquina (WALANTUS, 2014). 
 A equação que descreve o movimento do rotor de uma máquina síncrona está baseada 
no princípio elementar da dinâmica que diz ser o torque da aceleração igual ao produto do 
momento de inércia do rotor pela sua aceleração angular (STEVENSON, 1982), ou seja: 
 𝐽
𝑑2𝜃𝑚
𝑑𝑡2
= 𝑇𝑎 = 𝑇𝑚 − 𝑇𝑒 𝑁.𝑚 
Eq.(2.1) 
O torque mecânico tem origem no agente motor (água em hidrelétricas, vapor em 
termelétricas, por exemplo), e a potência elétrica exigida pelas cargas gera torques elétricos 
20 
 
 
através dos campos magnéticos (BRETAS e ALBERTO, 2000). O torque mecânico Tm é 
considerado positivo para o gerador síncrono, ou seja, atua no sentido de acelerar o rotor do 
gerador. Já o torque elétrico Te atua no sentido contrário. Isto signigica que Ta é o torque 
resultante do eixo do gerador, que tende a acelerar o rotor no sentido positivo da rotação de θm, 
quando o torque mecânico for maior que o elétrico e tende a desacelerar o rotor quando o torque 
mecânico for menor que o elétrico. 
Como θm é medido com respeito a um eixo de referência estacionária sobre o rotor, 
constitui-se em uma medida absoluta do ângulo do rotor. Consequentemente, cresce 
continuamente com o tempo, embora ainda em velocidade síncrona constante. Simplificando, 
o ângulo 𝜃𝑚 é variável no tempo e representa desvios do deslocamento angular do rotor em 
relação a posição angular síncrona. Como a velocidade do rotor relativa à síncrona é de 
interesse, é mais conveniente medir a posição angular do rotor com respeito a um eixo de 
referência que gira em velocidade síncrona (STEVENSON, 1982). Portanto, definimos a 
posição angular do rotor pela seguinte equação: 
𝜃𝑚 = 𝜔𝑠𝑚 + 𝛿𝑚 Eq.(2.2) Eq.(2.3) 𝜃𝑚 = (𝜔𝑟 − 𝜔𝑠)𝑡 + 𝛿𝑚 Eq.(2.3) 
Na qual ωsm é velocidade síncrona da máquina e δm é o deslocamento angular do rotor, 
em radianos mecânicos, a partir do eixo de referência da rotação síncrona. As derivadas da 
Equação (2.2) com respeito ao tempo são: 
𝑑𝜃𝑚
𝑑𝑡
= 𝜔𝑠𝑚 +
𝑑𝛿𝑚
𝑑𝑡
 
Eq.(2.3) 
𝑑2𝜃𝑚
𝑑𝑡2
=
𝑑2𝛿𝑚
𝑑𝑡2
 
Eq.(2.4) 
Realizando uma substituição de equações e lembrando da dinâmica elementar que 
potência é igual ao torque vezes a velocidade angular, temos: 
𝐽𝜔𝑠𝑚
𝑑2𝛿𝑚
𝑑𝑡2
= 𝑃𝑎 = 𝑃𝑚 − 𝑃𝑒 
Eq.(2.5) 
Em que, Pm é a potência de entrada do eixo da máquina menos perdas rotacionais, Pe é a 
potência elétrica no entreferro e Pa é a potência de aceleração
que leva em conta qualquer 
desequilibrio entre essas duas quantidades. 
Sabendo que 𝐽𝜔𝑠𝑚 é o momento angular do rotor e que na velocidade síncrona é 
representado por M, chamado de constante de inércia, sendo M: 
21 
 
 
𝑀 =
2𝐻
𝜔𝑠𝑚
𝑆𝑚𝑎𝑞 
Eq.(2.6) 
Em que: Smaq é a potência mecânica nominal da máquina em MVA e H é a relação entre 
a energia cinética armazenada em megajoule na velocidade síncrona pela potência nominal da 
máquina em MVA. Logo temos: 
2H
ωsm
d2δm
dt2
= Pa = Pm − Pe 
Eq.(2.7) 
Para um sistema com uma frequência elétrica de f Hz, sendo ωsm = 2𝜋𝑓 , a equação 
anterior torna-se: 
𝐻
𝜋𝑓
𝑑2𝛿𝑚
𝑑𝑡2
= 𝑃𝑎 = 𝑃𝑚 − 𝑃𝑒 
Eq.(2.8) 
Conforme STEVENSON (1982), a Equação (2.8), chamada de equação de oscilação da 
máquina, é a equação fundamental que governa as dinâmicas rotacionais das máquinas 
síncronas em estudos de estabilidade. Notamos que esta equação é diferencial de segunda ordem 
e que pode ser escrita como duas equações diferenciais de primeira ordem: 
2𝐻
𝜔𝑠𝑚
𝑑𝜔
𝑑𝑡
= 𝑃𝑚 − 𝑃𝑒 
Eq.(2.9) 
𝑑𝛿
𝑑𝑡
= 𝜔 − 𝜔𝑠𝑚 
Eq.(2.10) 
Segundo MASIERO, GURSKI e CASTRO (2016), as Equações (2.9) e (2.10) são 
utilizadas para descrever o comportamento dinâmico do sistema em estudo e, por apresentarem 
a equação de oscilação através de duas equações de primeira ordem, trazem maior facilidade 
quanto à implementação de métodos numéricos de solução. Observa-se que a primeira equação 
do sistema, Eq. (2.9), apresenta a primeira derivada do termo ω, ou a primeira derivada do 
desvio de velocidade do rotor, que representa a aceleração angular do rotor. Esse 
equacionamento então possibilita uma avaliação rápida do comportamento da máquina em 
decorrência das potências elétrica e mecânica, da seguinte forma: 
i) Para 𝑃𝑚 > 𝑃𝑒 , o gerador acelera; 
ii) Para 𝑃𝑚 < 𝑃𝑒 , o gerador desacelera. 
22 
 
 
2.1.2.1 Critério das Áreas Iguais 
De acordo com FANG e JI-LAI (2009), o critério das áreas iguais (CAI) pode ser usado 
para classificar rapidamente a estabilidade transitória, calcular rapidamente o tempo crítico de 
eliminação de falhas e o grau de estabilidade. A vantagem deste método é que ele evita a divisão 
de grupos de máquinas, a equivalência e a redução de rede, e precisa de menos informações, 
menos tempo de computação e tem a mesma precisão que a simulação no domínio do tempo na 
avaliação da estabilidade transitória. 
Combinando a curva característica do ângulo de potência relativo do par de geradores, o 
CAI é aplicado diretamente ao par de geradores para avaliar a estabilidade rapidamente. 
O sistema máquina-barra infinita e suas características de ângulo de potência são 
mostrados nas Figuras (2.2) e (2.3). Nele, 𝛿0 é o ângulo de equilíbrio inicial do rotor, 𝛿1 é o 
ângulo de equilíbrio do rotor após a falta, 𝛿𝑚á𝑥 é o ângulo do rotor pós-falta no ponto de 
equilíbrio instável e 𝛿𝑐𝑟𝑖𝑡 é o ângulo máximo para se liberar a falta sem que haja a perda de 
sincronismo. As curvas de potência elétrica durante os períodos pré-falta, pós-falta e durante a 
falta são representadas de acordo com a legenda multicor na Figura 2.3. 
 
Fonte: (FANG e JI-LAI, 2009). 
Figura 2.2 – Sistema máquina-barra infinita. 
 
23 
 
 
 
Fonte: Adaptado de (LOTERO, 2012). 
Figura 2.3 – Características da curva ângulo-potência 
 
As áreas A e B mostradas na Figura 2.3 são computadas da seguinte forma: 
𝐴 = ∫ (𝑃𝑚 − 𝑃𝑒
𝑑𝑢𝑟𝑎𝑛𝑡𝑒 𝑎 𝑓𝑎𝑙𝑡𝑎
)𝑑𝛿
𝛿𝑐𝑟𝑖𝑡
𝛿0
 
Eq.(2.11) 
𝐵 = ∫ (𝑃𝑒
𝑝ó𝑠−𝑓𝑎𝑙𝑡𝑎
 − 𝑃𝑚)𝑑𝛿
𝛿𝑚á𝑥
𝛿𝑐𝑟𝑖𝑡
 
Eq.(2.12) 
O CAI expressa que: se 𝐴 < 𝐵, o sistema máquina barra infinita é estável, se 𝐴 > 𝐵 o 
sistema é instável e se 𝐴 = 𝐵 o sistema é criticamente estável. 
As áreas A e B também podem ser caracterizadas por áreas de aceleração e frenagem, 
respectivamente. No momento em que a potência mecânica é maior que a potência elétrica a 
máquina está acelerando e quando é menor a máquina está desacelerando. 
24 
 
 
2.3 Modelagem matemática dos componentes necessários ao 
estudo de estabilidade transitória e fluxo de carga 
Neste item se abordará como os elementos da rede de um SEP são descritos 
matematicamente para a solução do problema de estabilidade. 
 
2.3.1 Máquina síncrona 
Para a representação da máquina síncrona no estudo de estabilidade transitória será 
adotado, neste trabalho, o denominado “modelo clássico”, no qual a máquina é modelada 
através de uma fonte de tensão atrás da reatância transitória. Essa tensão é constante em módulo, 
porém varia sua posição angular (SILVA, 2009). O modelo de gerador pode ser simplificado 
através da representação deste por uma fonte de tensão interna de amplitude constante (𝐸𝑔) 
conectado em série com a reatância transitória de eixo direto (𝑋𝑑′), conforme a Figura 2.4 
(GLOVER, SARMA e OVERBYE, 2011). 
Para um sistema de potência realístico, esse modelo é válido quando o tempo de 
simulação é pequeno e o regulador de tensão lento. Ainda assim, este modelo é bastante 
utilizado, devido sua facilidade de implementação e economia de tempo de processamento. 
Além disto, segundo (SILVA, 2009), os resultados obtidos podem ser considerados coerentes, 
de forma qualitativa, com os obtidos com modelos mais completos. 
 
Fonte: Adaptado de (GLOVER, SARMA e OVERBYE, 2011) 
Figura 2.4 – Representação do modelo clássico do gerador síncrono 
 
25 
 
 
𝐸′ = 𝑉𝑡 + 𝑗𝑋𝑑
′ 𝐼 Eq.(2.13) 
Em que: 
𝐸′ – Força eletromotriz da máquina 
𝑉𝑡 – Tensão nos terminais da máquina 
𝑋𝑑
′ – Reatância transitória de eixo direto do enrolamento do estator 
𝐼 – Corrente elétrica fornecida pela máquina 
δg – É o ângulo de carga do gerador, este é a defasagem entre a tensão interna (𝐸′) e a 
tensão terminal (Vt) do gerador 
Esta representação é obtida através das hipóteses simplificadoras listadas a seguir: 
i) A máquina opera sob condições de carga trifásica balanceada em sequência 
positiva; 
ii) Não há variação na excitação da máquina, ou seja, a amplitude de 𝐸′é constante 
durante o período transitório; 
iii) Parâmetros como perdas, saliência dos polos e saturação na máquina são 
desprezados. 
 
2.3.2 Carga 
De acordo com MATOS (2012) a curva de carga de um dado barramento de distribuição 
decorre de hábitos de consumo, temperatura, nível de renda, forma de tarifação, etc. Em vista 
disso, dentre os vários parâmetros de um sistema elétrico de potência, a carga dos consumidores 
é o de determinação mais difícil. 
A representação da carga depende muito do tipo de estudo realizado. A carga pode ser 
representada por potência constante, corrente constante ou impedância constante. É importante 
que se conheça a variação das potências ativas e reativas com a variação da tensão, pois embora 
seja exato considerar as características PV e QV de cada tipo de carga para simulação de fluxo 
de carga e estabilidade, o tratamento analítico é muito complicado (BORGES, 2005). Para os 
cálculos envolvidos, existem maneiras tradicionais para representação de carga, as quais serão 
apresentadas nesta seção. 
26 
 
 
2.3.2.1 Representação da carga para o fluxo de potência 
Para o algoritmo de fluxo de potência a carga foi definida como uma carga de potência 
constante, pois este modelo é invariante com a tensão e permite representar cargas que 
apresentam variação desprezível de potência. 
A Figura 2.5 mostra a representação da carga como potência ativa e reativa constantes. 
 
Fonte: (MATOS, 2012). 
Figura 2.5 – Modelo da carga para o fluxo de potência. 
 
𝑆𝑐𝑎𝑟𝑔𝑎 = 𝑃𝐿 + 𝑗𝑄𝐿 Eq.(2.14) 
 
2.3.2.2 Representação da carga para o estudo de estabilidade 
O problema mais crítico na análise de estabilidade transitória é quando a falta está 
próxima a barra do gerador, em razão da impedância vista do gerador diminuir com a 
proximidade da falta em relação a ele. Consequentemente, se torna interessante representar a 
carga como uma parte do sistema,
ou seja, como uma impedância constante, uma vez que o 
problema de estabilidade está relacionado com a impedância “vista” dos geradores. 
A Figura 2.6 indica a representação do tipo impedância constante. 
 
Fonte: (MATOS, 2012). 
Figura 2.6 – Modelo da carga de impedância constante. 
27 
 
 
2.3.3 Linhas de Transmissão 
Simplificadamente, as linhas de transmissão são os elementos que tem por função 
interligar os centros geradores aos centros consumidores (LUZ, 2015). 
Podemos classificar as linhas de transmissão como linhas curtas, médias ou longas. As 
linhas de transmissão curtas são aquelas de até 80 quilômetros, as linhas médias são aquelas 
que possuem de 80 a 240 quilômetros e as linhas longas são as maiores de 240 quilômetros 
(JÚNIOR e ROCHA, 2012). 
 
2.3.3.1 Linha Curta 
As linhas de transmissão são classificadas curtas quando a capacitância em derivação é 
tão pequena que pode ser inteiramente desprezada sem perda apreciável de precisão e é 
suficiente considerar apenas a resistência e a indutância em série (JÚNIOR e ROCHA, 2012). 
A Figura 2.7 representa seu modelo. 
 
Fonte: (MATOS, 2012). 
Figura 2.7 – Modelo de linha curta. 
2.3.3.2 Linha Média 
Já o modelo de linha média é composto diretamente por uma impedância série e uma 
admitância shunt total da linha, que é dividida em duas partes iguais colocadas junto às barras 
emissoras e receptoras, conforme Figura 2.8. Este circuito recebe o nome de π-nominal. 
28 
 
 
 
Fonte: Adaptado de (MATOS, 2012). 
Figura 2.8 – Modelo de linha média. 
 
2.3.3.3 Linha Longa 
Por fim, as linhas longas possuem uma modelagem considerando a indutância e a 
capacitância distribuídas uniformemente ao longo da linha de transmissão e podem ser 
representadas pela Figura 2.9, que apresenta os parâmetros de linha longa adaptados para o 
modelo π-nominal. 
 
Fonte: (MATOS, 2012). 
Figura 2.9 – Modelo de linha longa. 
 
Para a obtenção dos parâmetros de linha longa adaptados, as seguintes expressões podem 
ser utilizadas (MATOS, 2012): 
𝑍𝑒𝑞 = 𝑍 ×
𝑠𝑒𝑛ℎ(𝛾 × 𝑙)
𝛾 × 𝑙
 
 
Eq.(2.15) 
29 
 
 
𝑦𝑒𝑞 = 𝑌 ×
𝑡𝑎𝑛ℎ(𝛾 × 𝑙 2⁄ )
(𝛾 × 𝑙 2⁄ )
 
Eq.(2.16) 
Em que: 
𝛾 = √𝑧 × 𝑦 Eq.(2.17) 
𝑍 = 𝑧 × 𝑙 Eq.(2.18) 
𝑌 = 𝑦 × 𝑙 Eq.(2.19) 
Em que: Z é a impedância característica da linha e γ é a constante de propagação. Ambas 
γ e Z são quantidades complexas (WALANTUS, 2014). 
Como exposto, é feita uma correção do modelo π-nominal, de forma a se obter um modelo 
com parâmetros concentrados, mas que pode ser utilizado para linhas longas quando se tem o 
interesse somente nos valores de tensão e correntes nas extremidades da linha. 
Este modelo pode ser utilizado para representar linhas curtas, médias e longas. Devido a 
consideração da distribuição uniforme de seus parâmetros ele é mais exato para representa-las, 
portanto, foi o modelo adotado tanto para a o desenvolvimento do algoritmo do fluxo de carga 
como para a análise de estabilidade. 
 
2.3.4 Transformadores 
Cada gerador possui um transformador elevador associado, cujo modelo simplificado é 
apenas uma reatância série indutiva (𝑋𝑡). 
 
2.4 Modelagem da rede reduzida 
A matriz das impedâncias de barra e admitância de barra de um sistema elétrico real são 
muito grandes e tem dimensão da ordem de milhares de linhas. Nos estudos não é necessário 
se conhecer a tensão em todas as barras do sistema. Logo, seguem técnicas para reduzir a 
dimensão da rede, eliminando-se trechos não prioritários da rede para o estudo em questão 
(BORGES, 2005). 
30 
 
 
Particiona-se a matriz e ordenam-se as equações de tal forma que todas as barras sem 
fonte fiquem juntas e na parte inferior da matriz, conforme Figura 2.10 a seguir. 
 
Fonte: (BORGES, 2005). 
Figura 2.10 – Solução do sistema particionada. 
 
Supondo-se que 𝐼𝑏 = 0 vem: 
[
𝐼𝐴
𝐼𝐵
] = [
𝑌𝐴𝐴 𝑌𝐴𝐵
𝑌𝐴𝐵
𝑡 𝑌𝐵𝐵
] × [
𝑉𝐴
𝑉𝐵
] 
 
Eq.(2.20) 
𝐼𝐴 = 𝑌𝐴𝐴 × 𝑉𝐴 + 𝑌𝐴𝐵 × 𝑉𝐵, 
 
Eq.(2.21) 
𝐼𝐵 = 𝑌𝐴𝐵
𝑡 × 𝑉𝐴 + 𝑌𝐵𝐵 × 𝑉𝐵 = 0 → 𝑉𝐵 = 𝑌𝐵𝐵
−1 × 𝑌𝐴𝐵
𝑡 × 𝑉𝐴 
 
Eq.(2.22) 
Substituindo-se o valor de 𝑉𝐵 na equação de 𝐼𝐴 vem: 
𝐼𝐴 = 𝑌𝐴𝐴 × 𝑉𝐴 − 𝑌𝐴𝐵 × 𝑌𝐵𝐵
−1 × 𝑌𝐴𝐵
𝑡 × 𝑉𝐴 
 
Eq.(2.23) 
Agrupando-se em termos vem: 
𝐼𝐴 = ⌈𝑌𝐴𝐴 × 𝑉𝐴 − 𝑌𝐴𝐵 × 𝑌𝐵𝐵
−1 × 𝑌𝐴𝐵
𝑡 ⌉ × 𝑉𝐴 = [𝑌𝐴 × 𝑉𝐴] 
 
Eq.(2.24) 
A ordem da matriz 𝑌𝐴 neste exemplo é a do número de barras com fonte de corrente. 
Colocando-se de forma escalar, tem-se que a eliminação da barra n é: 
𝑌𝑖𝑗
′ = 𝑌𝑖𝑗 −
𝑌𝑖𝑗 × 𝑌𝑛𝑗
𝑌𝑛𝑛
 
 
Eq.(2.25) 
Que é chamada de eliminação de Kron. 
Para maior eficiência computacional deve-se evitar a inversão da matriz 𝑌𝐵𝐵. O 
procedimento é então o de eliminar uma barra por vez, aplicando-se a eliminação de Kron tantas 
31 
 
 
vezes quanto o número de barras a serem eliminadas (BORGES, 2005). Para o problema de 
estabilidade essa metodologia é aplicada k vezes, na qual k representa o número de barras de 
carga. Isso acontece porque na análise de estabilidade a dinâmica do sistema está relacionada 
aos geradores síncronos. 
 
2.5 Métodos numéricos de soluções de equações diferenciais 
Para determinar se um sistema é estável transitoriamente, primeiramente é necessária a 
modelagem matemática do mesmo em todos os momentos da ocorrência da falta (LUZ, 2015), 
uma vez que a dinâmica deste se modifica de acordo com a alteração da configuração da rede, 
durante a falta e após a atuação da proteção (MASIERO, GURSKI e CASTRO, 2016). 
Como a configuração da rede muda pela ação do curto-circuito, a potência fornecida, que 
está relacionada com a configuração da rede, também muda. Durante esse tempo, o 
comportamento do sistema é descrito pelas suas equações diferencias ordinárias (EDOs). 
Chama-se equação diferencial uma equação em que a incógnita é uma função e apresenta 
uma relação com as derivadas desta função. As equações diferenciais são utilizadas na 
resolução de problemas de modelagem matemática e quando a função desconhecida depende 
de uma única variável independente, são chamadas de equações diferenciais ordinárias (SILVA, 
2017). 
Há várias maneiras de se resolver analiticamente uma EDO, entretanto, cabe ressaltar que 
nem toda EDO pode ser resolvida analiticamente, pois os métodos analíticos são aplicáveis 
apenas a certas formas especiais de funções, mas toda EDO pode ser resolvida através da 
utilização de métodos numéricos (BRITO e AMARAL, 2012). 
Resumidamente, métodos numéricos são conjuntos detalhados e operações sequenciadas 
que nos levam a uma estimativa do problema. Quando estes são utilizados com o auxílio de 
computadores propiciam uma melhor agilidade no estudo das equações diferenciais, pois 
produzem tabelas e gráficos que possibilitam fazer um estudo mais detalhado dos dados obtidos 
na resolução. 
Na seção 2.1.2 foram desenvolvidas as equações de oscilação da barra swing, 
caracterizadas como EDOs. A solução formal de tais equações não pode ser explicitamente 
encontrada. Mesmo no caso de uma simples máquina oscilando com respeito a um barramento 
32 
 
 
infinito é bastante difícil obter soluções na forma literal e, portanto, métodos usando 
computação digital são normalmente usados. 
Para resolver a equação diferencial do problema proposto neste trabalho (estabilidade 
transitória), determinados métodos numéricos serão estudados. Além disso, será realizada uma 
comparação, baseada na literatura, entre esses métodos numéricos a fim de identificar o melhor 
artifício para a solução do problema. 
Além dos métodos numéricos para a solução do problema de estabilidade também será 
apresentado o método de Newton-Raphson para a determinação do fluxo de carga. 
 
2.5.1 Método de Euler 
Segundo LEITE (2006), o método de Euler é um procedimento numérico para aproximar 
a solução da equação diferencial: 
𝑥′ = 𝑓(𝑡, 𝑥) Eq.(2.26) 
que satisfaz a condição inicial 𝑥(𝑡0) = 𝑦0. (Admite-se que a equação
é suficientemente bem-
comportada, de forma a garantir a existência de uma única solução num intervalo que contenha 
o ponto 𝑡0 e os pontos 𝑡𝑖 usados na construção seguinte.) 
Sabe-se que o gráfico da solução passa pelo ponto (𝑡0 , 𝑥0) com inclinação igual a 𝑥
′(𝑡0) 
(ou seja, com inclinação igual a 𝑓(𝑡0 , 𝑥0)). Isto serve de ponto de partida para achar uma 
aproximação da solução. 
Começando pelo ponto (𝑡0 , 𝑥0), seguindo na direção dada pela inclinação. Usando um 
pequeno passo h, segue-se ao longo da reta tangente até chegar ao ponto (𝑡1 , 𝑥1), em que: 
𝑡1 = 𝑡0 + ℎ 𝑒 𝑥1 = 𝑥0 + ℎ𝑓(𝑡0 , 𝑥0) Eq.(2.27) 
O método de Euler consiste na repetição deste processo e gera a sucessão de pontos: 
𝑡𝑖+1 = 𝑡𝑖 + ℎ 𝑒 𝑥𝑖+1 = 𝑥𝑖 + ℎ𝑓(𝑡𝑖 , 𝑥𝑖) 𝑖 = 0,1,2… Eq.(2.28) 
Como este método considera apenas a primeira derivada de 𝑥, este é referido como 
método de primeira ordem. 
Uma representação do método de Euler é dada na Figura 2.11. 
Do desenvolvimento da série de Taylor tem-se um erro local como sendo: 
33 
 
 
𝜀ℎ ≅
(ℎ)2
2
𝑥′′ Eq.(2.29) 
Já o erro global vai se acumulando até que 𝜀 ≅ ℎ. Esta relação mostra que o erro tende 
a zero conforme h → 0 e, portanto, o método é convergente. De forma geral, um dado método 
é convergente quando: 
𝜀𝑖+1 = (ℎ)
𝑝 𝑝 > 0 
Eq.(2.30) 
Neste caso, o método possui ordem p, portanto, o método de Euler é um método de 
primeira ordem. Métodos de ordem superior convergem mais rapidamente para a solução exata 
(não exigem passos tão pequenos), sendo por isso mais indicados que os métodos de primeira 
ordem. 
A utilização de um passo muito pequeno (tendendo a zero) é impraticável, pois exigiria 
um tempo computacional demasiadamente alto, além de fazer com que os erros de 
arredondamento aumentassem rapidamente (FONTANA, 2018). 
 
Figura 2.11 – Representação gráfica do método de Euler 
 
2.5.2 Método Euler aprimorado (Fórmula de Heun) 
Este método faz parte da categoria de passos múltiplos, onde os valores da derivada em 
mais de um ponto são utilizados para determinar o próximo ponto (FONTANA, 2018). 
34 
 
 
O método de Euler possui uma velocidade de convergência baixa porque utiliza a 
derivada do início do intervalo como média para o intervalo todo. Para melhorar a estimativa 
da inclinação, o método de Heun determina duas derivadas para o intervalo. A primeira derivada 
é tomada no início do intervalo e equivale àquela usada no método de Euler. A segunda é 
tomada no ponto final do intervalo. A partir dessas duas derivadas, utiliza-se a média delas para 
obter uma estimativa melhorada da inclinação no intervalo (GARCIA, 2012). 
Do método de Euler, tínhamos que: 
𝑥𝑖
′ =
𝑑𝑥
𝑑𝑡
= 𝑓(𝑡𝑖, 𝑥𝑖) 
 Eq.(2.31) 
era a inclinação do início do intervalo. A fórmula de Heun é uma modificação simples do 
método de Euler e consiste em uma abordagem preditiva-corretiva. O passo inicial consiste 
determinar um valor preditor com base no método de Euler (FONTANA, 2018): 
𝑥𝑖+1
0 = 𝑥𝑖 + 𝑓(𝑡𝑖, 𝑥𝑖)ℎ Eq.(2.32) 
No método de Euler, essa equação é a usada para solucionar o problema das EDOs. 
Contudo, no método de Heun, esta estimativa é usada como aproximação inicial do valor, 
chamado de “predita”. Assim, essa equação é chamada de equação preditora, pois ela fornece 
uma estimativa de 𝑥𝑖+1 que permite o cálculo da inclinação na extremidade final do intervalo. 
A partir do valor obtido da Equação (2.32), é calculada a etapa corretora com base no 
valor médio da função avaliada em 𝑥𝑖 e no valor preditor 𝑥𝑖+1
0 : 
�̅�′ =
𝑥𝑖
′ + 𝑥𝑖+1
′
2
=
𝑓(𝑡𝑖, 𝑥𝑖) + 𝑓(𝑡𝑖+1, 𝑥𝑖+1
0 )
2
ℎ 
 Eq.(2.33) 
Essa inclinação média é então usada para extrapolar a fórmula de Euler: 
𝑥𝑖+1 = 𝑥𝑖 +
ℎ
2
(𝑓(𝑡𝑖, 𝑥𝑖) + 𝑓(𝑡𝑖+1, 𝑥𝑖+1
0 )) 
 Eq.(2.34) 
Esta fórmula representa uma correção do valor estimado anteriormente 𝑥𝑖+1
0 . 
Esta abordagem de Heun é chamada do tipo preditor-corretor e é apresentada na Figura 
2.12. 
35 
 
 
 
Fonte: Adaptado de (GARCIA, 2012). 
Figura 2.12 – .Curva das equações preditora e corretora 
 
 
2.5.3 Método de Runge-Kutta de quarta ordem 
O método de Runge-Kutta é provavelmente um dos métodos mais populares e o de quarta 
ordem é um dos mais utilizados para obter soluções aproximadas de valor inicial. Cada método 
de Runge-Kutta consiste em comparar um polinômio de Taylor apropriado para eliminar o 
cálculo das derivadas, fazendo-se várias avaliações da função a cada passo. O método de 
Runge-Kutta pode ser entendido como um aperfeiçoamento do método de Euler, com uma 
melhor estimativa da derivada da função. No método de Euler a estimativa do valor de 𝑡𝑖+1 é 
realizado com o valor de 𝑡𝑖 e com a derivada no ponto 𝑥𝑖. No método de Runge-Kutta, busca-
se uma melhor estimativa da derivada com a avaliação da função em mais pontos no intervalo 
[𝑥𝑖 , 𝑥𝑖+1] (STERZA e BRANDI, 2016). 
O método de Runge-Kutta, do mesmo modo que o método de Euler, pode ser obtido a 
partir da Série de Taylor. Entretanto, no caso do Método de Euler consideram-se apenas os dois 
primeiros termos da série. Já, nos Métodos de Runge-Kutta, considera-se uma quantidade maior 
36 
 
 
de termos da série. Por esse motivo, em geral, quando aplicados, estes últimos levam a 
resultados mais precisos (RISO e NOTARE, 2011). 
Polinômio de Taylor: 
𝑝(𝑥) = 𝑓(𝑎) + 𝑓′(𝑎)
(𝑥 − 𝑎)1
1!
+ 𝑓′′ (𝑎)
(𝑥 − 𝑎)2
2!
+ ⋯+ 𝑓𝑛(𝑎)
(𝑥 − 𝑎)𝑛
𝑛!
 Eq.(2.35) 
De modo geral, o método de Runge-Kutta é dado por: 
𝑡𝑖+1 = 𝑡𝑖 + ℎ(𝑖𝑛𝑐𝑙𝑖𝑛𝑎çã𝑜) Eq.(2.36) 
em que a i𝑛𝑐𝑙𝑖𝑛𝑎cã𝑜 é uma constante que é obtida através do cálculo da inclinação em vários 
pontos no interior do subintervalo. A ordem do método indica o número de pontos usado em 
um subintervalo para determinar o valor da i𝑛𝑐𝑙𝑖𝑛𝑎cã𝑜, isto é, o método de Runge-Kutta de 
segunda ordem utiliza a inclinação em dois pontos, o método de terceira ordem utiliza três 
pontos, e assim por diante (STERZA e BRANDI, 2016). 
O procedimento de Runge-Kutta de quarta ordem consiste em encontrar constantes 
apropriadas de tal forma que: 
𝑡𝑖+1 = 𝑡𝑖 + ℎ(𝑐1𝐾1 + 𝑐2𝐾2 + 𝑐3𝐾3 + 𝑐4𝐾4) Eq.(2.37) 
sendo 𝑐1, 𝑐2, 𝑐3 e 𝑐4 constantes a serem determinadas e 
𝑘1 = 𝑓(𝑥𝑖, 𝑡𝑖) Eq.(2.38) 
𝑘2 = 𝑓(𝑥𝑖 + 𝛼1ℎ, 𝑡𝑖 + 𝛽1ℎ𝑘1) Eq.(2.39) 
𝑘3 = 𝑓(𝑥𝑖 + 𝛼2ℎ, 𝑡𝑖 + 𝛽2ℎ𝑘1 + 𝛽3ℎ𝑘2) Eq.(2.40) 
𝑘4 = 𝑓(𝑥𝑖 + 𝛼3ℎ, 𝑡𝑖 + 𝛽4ℎ𝑘1 + 𝛽5ℎ𝑘2 + 𝛽6ℎ𝑘3) Eq.(2.41) 
Igualando a Equação (2.34) com o polinômio de Taylor de grau 4, resultará um sistema 
não-linear de 11 equações e 13 incógnitas, ou seja, haverá infinitas soluções, onde as 13 
incógnitas são 𝑐1, 𝑐2, 𝑐3, 𝑐4, 𝛼1, 𝛼2, 𝛼3, 𝛽1, 𝛽2, 𝛽3, 𝛽4, 𝛽5 e 𝛽6 (STERZA e BRANDI, 2016). 
No entanto há uma solução, em particular, que é muito utilizada, fazendo: 
𝑐1 = 𝑐4 =
1
6⁄ ; 𝛼3 = 𝛽6 = 1; 
 𝑐2 = 𝑐3 =
1
3⁄ ; 𝛼1 = 𝛼2 = 𝛽1 = 𝛽3 =
1
2⁄ ; 𝑐2 = 𝑐3 =
1
3⁄ ; 𝛼1 = 𝛼2 = 𝛽1 = 𝛽3 =
1
2⁄ ; 𝑐2 = 𝑐3 =
1
3⁄ ; 𝛼1 = 𝛼2 = 𝛽1 = 𝛽3 =
1
2⁄ ; 
 
 𝛽2 = 𝛽4 = 𝛽5 = 0; 
 
37 
 
 
Dessa forma a solução de x para o método de Runge-Kutta de quarta ordem é calculada 
através da seguinte equação: 
𝑥𝑖+1 = 𝑥𝑖 +
1
6
(𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4) 
 
Eq.(2.42) 
Em que 𝑥𝑖+1 é a aproximação de 𝑥(𝑡𝑖+1) e 
𝑘1 = 𝑓(𝑥𝑖, 𝑡𝑖) Eq.(2.43) 
𝑘2 = 𝑓 (𝑥𝑖 +
1
2
ℎ, 𝑡𝑖 +
1
2
ℎ𝑘1) 
Eq.(2.44) 
𝑘3 = 𝑓 (𝑥𝑖 + +
1
2
ℎ, 𝑡𝑖 + +
1
2
ℎ𝑘2) 
Eq.(2.45) 
𝑘4 = 𝑓(𝑥𝑖 + ℎ, 𝑡𝑖 + ℎ𝑘3) Eq.(2.46) 
Em conformidade com KUNDUR (1994), a interpretação física dessas soluções é: 
 𝑘1 – (Inclinação no início da etapa) 
 𝑘2 – (Primeira aproximação da inclinação na etapa intermediária,
usando a inclinação 
𝑘1 para determinar o valor de x no ponto 𝑡𝑖 + ℎ/2) 
 𝑘3 – (Segunda aproximação da inclinação na etapa intermediária, usando a inclinação 
𝑘2 para determinar o valor de x no ponto 𝑡𝑖 + ℎ) 
 𝑘4 – (Inclinação no final da etapa) 
Na qual ℎ é o valor incremental de x, dado pela média ponderada das estimativas baseadas 
nas inclinações, no início, meio e fim, do intervalo de tempo (DIAS e PILONI, 2010), conforme 
Figura 2.13. 
38 
 
 
 
 Inclinações usadas pelo método de Runge-Kutta de quarta ordem 
 
No método de Runge-Kutta os incrementos nos valores das variáveis dependentes são 
calculados a partir de um conjunto de fórmulas expressas, em termos das derivadas calculadas 
num conjunto de pontos pré-fixados. Uma vez que cada valor de 𝑥 é determinado pelas fórmulas 
de uma maneira unívoca, este método não requer repetidas aproximações. O seu principal 
inconveniente é o tempo que demora a calcular várias vezes a lª derivada, sobretudo se a sua 
expressão for complexa. Outro inconveniente é não se poder saber, ao longo do cálculo, os 
valores dos erros. A principal vantagem reside em que é suficiente conhecer o valor da função 
num único ponto (o ponto inicial), para se poder determinar os seus valores nos pontos seguintes 
(BARBOSA, 2013). 
39 
 
 
2.5.4 Comparação dos métodos 
O método de Euler tem uma aplicação bastante restrita devido à sua pequena precisão, 
isto é, normalmente precisamos escolher um passo h muito pequeno para obter soluções de boa 
qualidade, o que implica um número elevado de passos e, consequentemente, alto custo 
computacional (JUSTO, SAUTER, et al., 2018). 
O método de Heun é também simples de aplicar e tem a vantagem de se poder controlar 
o grau de aproximação das sucessivas aproximações obtidas para 𝑥. Tem, contudo, também 
uma precisão limitada e exige pequenos intervalos h para a variável independente (BARBOSA, 
2013). 
 O método de Runge-Kutta (4ª ordem) exige um elevado número de operações 
aritméticas, mas os resultados em geral são mais precisos, pois são utilizados mais elementos 
da Série de Taylor na sua implementação. Portanto, é o método adotado para este trabalho. 
2.6 Fluxo de Potência 
O cálculo de fluxo de potência em uma rede de energia elétrica consiste essencialmente 
na determinação do estado da rede, da distribuição dos fluxos e de algumas outras grandezas de 
interesse (MONTICELLI, 1983). 
Os estudos de fluxo de carga possuem grande importância no planejamento da expansão 
futura de sistemas de potência como também na determinação da melhor operação de sistemas 
existentes. Este trabalho utiliza o fluxo de carga para encontrar as condições pré-falta do 
sistema. A principal informação obtida do estudo de fluxo de carga é o módulo e o ângulo de 
fase da tensão em cada barra e as potências ativa e reativa que circulam em cada linha 
(STEVENSON, 1986). 
O objetivo do fluxo de potência é encontrar as magnitudes e fases das tensões de cada 
barra em um sistema a partir dos dados de cargas ativas e reativas, da tensão dos geradores e 
das potências ativas dos geradores (MATOS, 2012). 
A solução para o problema de fluxo de potência começa com a identificação das variáveis 
conhecidas e desconhecidas no sistema, isto depende da definição do tipo de barramento. Uma 
barra sem geradores conectados a ela é chamada de barra de carga (barra PQ). Barras com pelo 
menos um gerador conectado são chamadas de barras de geração (barra PV). Para que o 
conjunto de equações que resolvem esse sistema tenha solução, é necessário arbitrar a um 
40 
 
 
barramento módulo e ângulo da tensão para ser a referência do sistema (barra Vθ), este 
barramento deve possuir um gerador e é chamado de barra de referência (MATOS, 2012). 
Para cada barra de carga, tanto a magnitude da tensão e ângulo são desconhecidas e devem 
ser encontrados. Para cada barra geradora, o ângulo de tensão deve ser encontrado. 
Há diversos métodos computacionais para a determinação do fluxo de potência, 
entretanto neste trabalho se optou pelo método de Newton-Raphson, Anexo D. 
Como descrito, o método de Newton-Raphson utiliza como base a expansão da série de 
Taylor. Para solucionar o problema de fluxo de carga, as equações de potência ativa e reativa a 
serem resolvidas possuem a seguinte estrutura: 
𝑃𝑘
𝑐𝑎𝑙𝑐 = ∑|𝑉𝑘 × 𝑉𝑚 × 𝑌𝑘𝑚| cos(𝜃𝑘𝑚 + 𝛿𝑘 − 𝛿𝑚)
𝑁
𝑛=1
 
 
Eq.(2.47) 
𝑄𝑘
𝑐𝑎𝑙𝑐 = − ∑|𝑉𝑘 × 𝑉𝑚 × 𝑌𝑘𝑚| sen(𝜃𝑘𝑚 + 𝛿𝑘 − 𝛿𝑚)
𝑁
𝑛=1
 
 
Eq.(2.48) 
A barra de referência (barra 𝑉𝜃) é omitida da solução iterativa para determinar tensões, 
pois tanto o módulo quanto o ângulo da tensão são especificados. São especificados P e Q em 
todas as barras, exceto a barra 𝑉𝜃, e será estimado o módulo e ângulo em todas as barras exceto 
a barra 𝑉𝜃. Utilizam-se estes valores estimados para calcular os valores de 𝑃𝑘
𝑐𝑎𝑙𝑐 e 𝑄𝑘
𝑐𝑎𝑙𝑐, e 
definem - se: 
∆𝑃𝑘 = 𝑃𝑘
𝑒𝑠𝑝 − 𝑃𝑘
𝑐𝑎𝑙𝑐 Eq.(2.49) 
∆𝑄𝑘 = 𝑄𝑘
𝑒𝑠𝑝 − 𝑄𝑘
𝑐𝑎𝑙𝑐
 Eq.(2.50) 
em que, 𝑃𝑘
𝑒𝑠𝑝
e 𝑄𝑘
𝑒𝑠𝑝
 representam as potências especificas nas barras de carga e nos 
geradores. 
Para a solução do problema de fluxo de carga é necessária a montagem da matriz 
jacobiana, que advém das derivadas parciais de P e Q em relação a cada uma das variáveis da 
equação 2.47 e 2.47. A Figura 2.15 demonstra a matriz jacobiana de um sistema composto por 
3 barras onde a barra 1 é a barra de oscilação, dai, os elementos da matriz representados pelos 
índices 2 e 3. 
41 
 
 
Assim, são definidos os subsistemas 1 e 2, dos quais, o primeiro trata do cálculo da tensão 
e do ângulo das barras e o segundo calcula as potências ativa e reativa na barra de referência, e 
potências reativas nas barras PV. A dimensão do subsistema 1 consiste do número de barras PV 
mais duas vezes o número de barras PQ. A dimensão do subsistema 2 é número de barras PV 
mais 2. 
 
Fonte: (MATOS, 2012) 
Figura 2.13 – Exemplo de matriz jacobiana. 
 
Equações que envolvem mais barras são resolvidas invertendo a jacobiana. A matriz 
jacobiana representa a derivada primeira (inclinação) da Equação (2) do Anexo D. Os valores 
encontrados para Δk e ∆|𝑉𝑘| são adicionados aos valores anteriores do módulo e ângulo da 
tensão, de modo a se obter os novos valores para 𝑃𝑘,𝑐𝑎𝑙𝑐 
(1)
e 𝑄𝑘,𝑐𝑎𝑙𝑐
(1)
, e começar a próxima iteração. 
O proceso é repetido até que a precisão escolhida seja satisfeita quando aplicada a todos os 
valores de cada coluna da matriz. Entretanto, para melhorar a convergência, as estimativas 
iniciais da tensão devem ser razoáveis.
 
 
42 
 
Capítulo 3 – Estruturação do problema e sistemática do 
algoritmo 
Equações de um sistema multimáquinas podem ser escritas similarmente a de um sistema 
conectado a barra infinita. A fim de reduzir a complexidade da análise da estabilidade 
transitória, suposições simplificadoras similares são feitas da seguinte maneira: 
– Presume-se que as potências de entrada permaneçam constantes durante todo o período 
de simulação; 
– Usando as tensões de barramento pré-faltas, todas as cargas são convertidas em 
admitância equivalentes à terra e são assumidas constantes. 
 
3.1 Resolução do Fluxo de potência 
Dadas as suposições, o primeiro passo na determinação de análise de estabilidade 
transitória é resolver o fluxo de potência, para determinar as magnitudes e ângulos de fase das 
tensões nas barras. Após, são determinadas as correntes dos geradores antes da perturbação, 
que são calculadas a partir de: 
𝐼𝑖 =
𝑆𝑖
∗
𝑉𝑖
∗ =
𝑃𝑖 − 𝑗𝑄𝑖
𝑉𝑖
∗ , 𝑖 = 1,2, … ,𝑚 
 
Eq.(3.1) 
Em que: 
– m é o número de geradores; 
– V𝑖 é a tensão no terminal do i-ésimo gerador; 
– 𝑃𝑖 e 𝑄𝑖 são as potências ativa e reativa do gerador. 
Determinado o valor da corrente dos geradores e desprezando as resistências de armadura, 
a força eletromotriz de cada gerador é obtida como consequência da equação (2.13). 
Em
seguida, é montada a matriz 𝑌𝑏𝑎𝑟𝑟𝑎 do sistema para a análise do problema de 
estabilidade transitória, que é diferente da matriz utilizada no fluxo de potência, pois no tema 
de estabilidade transitória os estudos são realizados baseados na f.e.m do gerador e não na 
tensão do seu barramento terminal. 
43 
 
 
 
Os elementos diagonais da matriz 𝑌𝑏𝑎𝑟𝑟𝑎 são a soma das admitâncias conectadas a ele, e 
os elementos não pertencentes a diagonal representam as admitâncias entre os nós. Para incluir 
a força eletromotriz interna das máquinas na equação, nos nós adicionais são inclusas as 
reatâncias transitórias das máquinas. Segundo SAADAT (1999), para inclui-las à solução do 
problema, m barramentos são adicionadas a rede do sistema de potência de n barramentos. A 
rede equivalente, já com as cargas convertidas em admitâncias é representada na Figura 3.1. 
Todas as cargas são convertidas para admitância equivalentes e adicionadas a matriz 
𝑌𝑏𝑎𝑟𝑟𝑎, usando a relação: 
𝑦𝑖0 = 
𝑆𝑖
∗
⌊𝑉𝑖⌋2
= 
𝑃𝑖 − 𝑗𝑄𝑖
⌊𝑉𝑖⌋2
 
 
Eq.(3.2) 
Na qual o índice i indica o número da barra. 
 
Fonte: Adaptado de (SAADAT, 1999) 
Figura 3.1 – Representação do SEP para análise de estabilidade transitória. 
 
Os nós n+1, n+2, ..., n+m representam os barramentos internos da máquina, ou seja, os 
barramentos por trás das reatâncias transitórias. A equação de tensão dos nós para este sistema, 
considerando o nó zero como referência é: 
44 
 
 
 
 
[
 
 
 
 
 
 
 
 
 
 
 
𝐼1
𝐼2
.
.
.
𝐼𝑛
_____
𝐼𝑛+1
.
.
.
𝐼𝑛+𝑚]
 
 
 
 
 
 
 
 
 
 
 
=
[
 
 
 
 
 
 
 
 
 
 
 
 
𝑌11 … 𝑌1𝑛 𝑌1(𝑛+1) … 𝑌1(𝑛+𝑚)
𝑌21 … 𝑌2𝑛 𝑌2(𝑛+1) … 𝑌2(𝑛+𝑚)
. … . . … .
. … . . … .
. … . . … .
𝑌21 … 𝑌2𝑛 𝑌2(𝑛+1) … 𝑌2(𝑛+𝑚)
__________________________________________________________
𝑌(𝑛+1)1 … 𝑌(𝑛+1)𝑛 𝑌(𝑛+1)(𝑛+1) … 𝑌(𝑛+1)(𝑛+𝑚)
. … . . … .
. … . . … .
. … . . … .
𝑌(𝑛+𝑚)1 … 𝑌(𝑛+𝑚)𝑛 𝑌(𝑛+𝑚)(𝑛+1) … 𝑌(𝑛+𝑚)(𝑛+𝑚)]
 
 
 
 
 
 
 
 
 
 
 
 
[
 
 
 
 
 
 
 
 
 
 
 
𝑉1
𝑉2
.
.
.
𝑉𝑛
_____
𝐸𝑛+1
′
.
.
.
𝐸𝑛+𝑚
′ ]
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Eq.(3.3) 
ou, 
𝐼𝑏𝑎𝑟𝑟𝑎 = 𝑌𝑏𝑎𝑟𝑟𝑎𝑉𝑏𝑎𝑟𝑟𝑎 Eq.(3.4) 
Em que 𝐼𝑏𝑎𝑟𝑟𝑎 é o vetor de correntes injetadas e 𝑉𝑏𝑎𝑟𝑟𝑎 é o vetor de tensões medidas a 
partir do nó de referência. 
Como citado por MASIERO, GURSKI e CASTRO (2016), para determinar se um sistema 
é estável transitoriamente, primeiramente é necessária a modelagem matemática do mesmo em 
todos os momentos da ocorrência da falta, uma vez que a dinâmica deste se modifica de acordo 
com a alteração da configuração da rede durante a falta e após a atuação da proteção. 
Desta forma, é comum dividir a análise em três períodos distintos: período pré-falta, 
período em falta e período pós-falta, conforme a seguir. 
 
3.1.1 Pré-falta 
 O período pré-falta é definido como o instante de operação do sistema imediatamente 
antes da ocorrência de uma grande perturbação. Nesse momento o sistema opera em sua 
condição nominal. 
A Figura 3.2 representa um sistema em condições nominais de operação. 
 
45 
 
 
 
 
Fonte: Adaptado de (KUNDUR, 1994) 
Figura 3.2 – Sistema em condições normais de operação 
 
3.1.2 Durante a falta 
O período em falta compreende desde a ocorrência de uma falta até a atuação dos 
equipamentos de proteção que a isolam. Neste caso, as características do sistema se alteram de 
suas condições nominais. 
Esse período pode ser caracterizado por dois tipos de falta, sendo elas: falta na barra ou 
falta na linha de transmissão. 
De acordo com KUNDUR (1994), o estudo clássico de estabilidade transitória é baseado 
na aplicação de uma falta trifásica. Uma falta trifásica sólida no barramento k na rede resulta 
em 𝑉𝑘 = 0. Isso é simulado removendo-se a linha e a coluna k da matriz 𝑌𝑏𝑎𝑟𝑟𝑎 de pré-falta. 
Diferentemente, no caso de uma falta trifásica na linha de transmissão não se resulta em 
𝑉𝑘 = 0, assim, neste caso não ocorre a eliminação da linha e coluna do barramento k. Para 
representar a nova matriz 𝑌𝑏𝑎𝑟𝑟𝑎 no momento da falta em uma linha de transmissão, é retirada 
a impedância existente entre o nós das faltas da matriz 𝑌𝑏𝑎𝑟𝑟𝑎 pré-falta e, conforme o ponto da 
falta, a impedância é distribuída entre as barras, representada como uma admitância shunt. 
A Figura 3.3 representa um sistema para uma falta no meio da linha. 
 
Fonte: Adaptado de (KUNDUR, 1994) 
Figura 3.3 – Sistema de potência com falta no meio da linha 
46 
 
 
 
3.1.3 Pós-falta 
 No período pós-falta, a falta é eliminada com a atuação dos dispositivos de proteção e o 
sistema passa a operar com novas características. 
Quando a falha é eliminada, pode envolver a remoção da linha defeituosa e a matriz 
𝑌𝑏𝑎𝑟𝑟𝑎 é recalculada para refletir a alteração nas redes. A eliminação da linha consiste em zerar 
a impedância 𝑍𝑖𝑗 existente entre as barras de falta e posteriormente montar a matriz 𝑌𝑏𝑎𝑟𝑟𝑎 após 
a falta, conforme a Figura 3.4. 
 
Fonte: Adaptado de (KUNDUR, 1994) 
Figura 3.4 – Abertura dos disjuntores e remoção da linha 
 
3.2 Determinação da Estabilidade do SEP 
Determinadas as matrizes pré-falta, durante a falta e pós-falta é aplicada a redução de 
Kron as barras de carga e as novas matrizes reduzidas são calculadas, eliminando todos os nós, 
exceto os nós acoplados aos geradores. Com isso, é determinada a potência elétrica de saída de 
cada máquina, para cada período (pré-falta, durante a falta e após a falta). Essa potência é 
expressa em termos de suas tensões internas: 
𝑆𝑒𝑖 = 𝐸𝑖
′∗𝐼𝑖 Eq.(3.5) 
ou, 
𝑃𝑒𝑖 = 𝑅𝑒(𝐸𝑖
′∗𝐼𝑖) 
 
Eq.(3.6) 
Em que, 
47 
 
 
 
𝐼𝑖 = ∑𝐸𝑗
′𝑌𝑖𝑗
𝑚
𝑗=1
 
Eq.(3.7) 
Enfim, expressando as tensões e admitâncias na forma polar, resulta-se em: 
𝑃𝑒𝑖 = ∑|𝐸𝑖
′||𝐸𝑗
′||𝑌𝑖𝑗|cos (𝜃𝑖𝑗 − 𝛿𝑖 − 𝛿𝑗)
𝑚
𝑗=1
 
Eq.(3.8) 
Tratando as tensões de excitação dos geradores durante a falta e pós-falta como 
constantes, ou seja, iguais as tensões de pré-falta, a energia elétrica do i-ésimo gerador em 
termos das novas matrizes admitâncias reduzidas podem ser obtidas da equação (3.8). 
Para facilitar o entendimento, consideremos o sistema da figura 3.2, com uma falta no 
meio da linha, uma potência fornecida pela máquina igual a 1 p.u., a f.e.m do gerador 1,05 p.u., 
a tensão da barra infinita 1 p.u e a constante do gerador H = 5 MJ/MVA. Com esses dados têm-
se as seguintes equações de potência elétrica: 
𝑃𝑒𝑙é𝑡𝑟𝑖𝑐𝑜
𝑝𝑟é−𝑓𝑎𝑙𝑡𝑎
= 2,1 sin 𝛿 Eq.(3.9) 
𝑃𝑒𝑙é𝑡𝑟𝑖𝑐𝑜
𝑑𝑢𝑟𝑎𝑛𝑡𝑒−𝑓𝑎𝑙𝑡𝑎
= 0,808 sin 𝛿 Eq.(3.10) 
𝑃𝑒𝑙é𝑡𝑟𝑖𝑐𝑜
𝑝𝑟é−𝑓𝑎𝑙𝑡𝑎
= 1,5 sin 𝛿 Eq.(3.11) 
Para determinar a estabilidade do sistema é aplicada a equação de oscilação da máquina 
síncrona, Equação (2.13), equação fundamental que governa as dinâmicas rotacionais das 
máquinas síncronas. Entretanto, para a solução do problema de análise transitória pelo método 
de Runge-Kutta, é necessário representar essa equação por duas equações de estado para cada 
gerador, de maneira a se obter a função do deslocamento angular e da velocidade angular em 
função do tempo: 
2𝐻
𝜔𝑠
𝑑𝜔𝑖
𝑑𝑡
= 𝑃𝑚𝑖 − 𝑃𝑒𝑖 𝑝. 𝑢. 
 Eq.(3.12) 
𝑑𝛿𝑖
𝑑𝑡
= 𝜔𝑖 − 𝜔𝑠 
Eq.(3.13) 
Com as condições iniciais, para o caso da Figura 3.3. 
48 
 
 
 
{
𝛿(0) = 0,496 𝑟𝑎𝑑
𝜔(0) = 0 𝑟𝑎𝑑/𝑠
 
 Eq.(3.14) 
Tem-se as seguintes equações: 
60𝜋
5
𝑑𝜔𝑖
𝑑𝑡
= 1 − 2,1 sin 𝛿 Eq.(3.15) 
60𝜋
5
𝑑𝜔𝑖
𝑑𝑡
= 1 − 0,818 sin 𝛿 Eq.(3.16) 
60𝜋
5
𝑑𝜔𝑖
𝑑𝑡
= 1 − 1,5 sin 𝛿 Eq.(3.17) 
Agora é aplicado o método de Runge-Kutta, o qual basicamente consiste em calcular cada 
ponto da função com o valor atual mais um conjunto de inclinações representadas por 𝑘1, 𝑘2, 𝑘3 
e 𝑘4. 
No algoritmo (Anexo A) as inclinações 𝑘1, 𝑘2, 𝑘3 e 𝑘4 tratam-se como vetores coluna de 
dimensão 2x1, sendo que os elementos das posições (i,1) representam as inclinações para a 
equação de 𝛿 e os elementos da posição (i+1,1) representam as inclinações para as equações de 
𝜔. 
Em que, 𝑖 = 1,3,5, … 
Os tempos para cada equação são discriminados, sendo que o período em que a equação 
3.15 é resolvida é para quando o tempo percorrido no algoritmo é menor que o tempo de atuação 
da falta. Já a equação 3.16 é solucionada quando o tempo percorrido é maior ou igual ao tempo 
de falta e menor que o tempo de abertura (eliminação) da falta. Por fim, para o desenvolvimento 
da equação 3.17 o tempo percorrido no algoritmo deve ser maior que o tempo de abertura. 
Com isso, a simulação opera para determinar a estabilidade do sistema, até que as parcelas 
revelem uma tendência definida quanto à estabilidade ou instabilidade. Normalmente, a solução 
é realizada para duas oscilações para mostrar que a segunda oscilação não é maior que a 
primeira. Se as diferenças de ângulo não aumentam, o sistema é estável. Se alguma das 
diferenças de ângulo aumentar indefinidamente, o sistema fica instável (DUY, MINH e LOC, 
2000). 
O fluxograma apresentado na Figura 3.5 retrata a sistemática apresentada e assessora em 
seu entendimento e define a modelagem do algoritmo. 
 
49 
 
 
 
 
Figura 3.5 – Fluxograma do algoritmo de estabilidade transitória.
 
 
50 
 
Capítulo 4 – Estudo de caso 
Neste capítulo serão apresentadas análises de estabilidade transitória para diferentes 
cenários de simulações. 
4.1 Introdução 
A estabilidade é avaliada pela variação angular de cada gerador com todos os demais 
geradores do sistema, de acordo com SILVA (2009). Para que o caso seja aprovado (sistema 
estável) todos os geradores deverão permanecer em sincronismo e as variações angulares entre 
as máquinas deverão ficar restritas a um certo limite. Segundo KUNDUR (2009), acima de um 
certo limite, um aumento na separação angular entre os geradores é acompanhado de uma 
diminuição na transferência de potência. Isso aumenta ainda mais a separação angular e leva a 
instabilidade. 
O tempo máximo para eliminar a falta de forma que o sistema continue estável é definido 
como tempo crítico de abertura. Dessa forma, a eliminação do defeito antes do tempo crítico 
leva a um sistema estável, e após o tempo crítico a um sistema instável (BRETAS e ALBERTO, 
2000). Logo, o tempo que os equipamentos de proteção levam para extinguir uma falta é 
essencial para conservar uma condição estável de operação de um sistema. Neste item serão 
determinados esses tempos para diferentes cenários. 
Com a abertura da proteção, o período de ocorrência da falta é na faixa de milissegundos, 
dessa forma não há tempo suficiente para a atuação do regulador de frequência da máquina 
primária (LUZ, 2015), portanto o regulador foi desconsiderado neste trabalho. 
 
4.2 Configuração do sistema teste 
Para realizar as simulações foi escolhido um sistema de 6 barras e 2 geradores. O 
diagrama unifilar do sistema é apresentado na figura 4.1. 
51 
 
 
 
 
Fonte: Adaptado de (MISHRA e UMREDKAR, 2013) 
Figura 4.1 – Sistema teste seis barras, duas máquinas e barra infinita 
 
Para obter os dados do fluxo de potência e realizar a análise de estabilidade é necessário 
se ter os dados das interligações e das barras. Tais dados são apresentados nas tabelas 4.1, 4.2 
e 4.3. 
 
Tabela 4.1 – Dados de carga 
Barra Carga (MW) Carga (Mvar) 
1 0 0 
2 0 0 
3 0 0 
4 100 70 
5 90 30 
6 160 110 
Fonte: (MISHRA e UMREDKAR, 2013) 
Tabela 4.2 – Dados de linha incluindo os transformadores 
De Barra Para Barra R (p.u) X (p.u) B/2 (p.u) 
1 4 0,035 0,225 0,0065 
1 5 0,025 0,105 0,0045 
1 6 0,400 0,215 0,0055 
2 4 0,000 0,035 0,0000 
3 5 0,000 0,042 0,0000 
4 6 0,028 0,125 0,0035 
5 6 0,026 0,175 0,0300 
Fonte: (MISHRA e UMREDKAR, 2013) 
52 
 
 
 
Tabela 4.3 – Dados de gerador 
Barra Magnitude da 
tensão (p.u) 
Geração (MW) Impedância 
transitória (p.u) 
H 
(MW.s/MVA) 
1 1,06 --- --- --- 
2 1,04 150 0,15 10 
3 1,03 100 0,10 5 
Fonte: Adaptado de (MISHRA e UMREDKAR, 2013) 
A partir daí o programa NR6BARRAS desenvolvido em MATLAB e apresentado no 
ANEXO B se realiza o cálculo do fluxo de carga, gerando os dados apresentados na tabela 4.4. 
Tabela 4.4 – Resultados do fluxo de carga 
Barra 
Tipo 
Tensão Geração Carga 
Módulo 
(p.u) 
Ângulo 
(graus) 
Ativo 
(p.u) 
Reativo 
(p.u) 
Ativo 
(p.u) 
Reativo 
(p.u) 
1 𝑉𝜃 1,0600 0,0000 1,1202 0,7526 0,000 0,000 
2 𝑃𝑉 1,0400 1,2012 1,5000 1,2734 0,000 0,000 
3 𝑃𝑉 1,0300 0,3857 1,0000 0,5565 0,000 0,000 
4 𝑃𝑄 0,9984 -1,8068 0,000 0,000 -1,0000 -0,7000 
5 𝑃𝑄 1,0081 -2,0207 0,000 0,000 -0,9000 -0,3000 
6 𝑃𝑄 0,8903 -6,3938 0,000 0,000 -1,6000 -1,1000 
 
4.3 – Eventos simulados 
Com base nos resultados do fluxo de carga se torna possível realizar a análise de 
estabilidade representada pelo algoritmo no Anexo A. Para tanto foram consideradas 
simulações para quatro tipos de falta: 
– Curto-circuito trifásico em uma barra, havendo o desligamento permanente de uma 
linha de transmissão. 
– Curto-circuito trifásico em uma barra, com eliminação do curto sem a retirada de uma 
linha. 
– Curto-circuito na metade da linha de transmissão, havendo o desligamento permanente 
de uma linha. 
– Curto-circuito na metade da linha de transmissão, com eliminação do curto sem a 
retirada de uma linha. 
53 
 
 
 
Os quatro tipo de simulações foram realizados para um passo de integração de 0,0025s, 
um tempo de simulação de 2,5s e uma margem de teste de 100ms entre cada intervalo. 
 
4.3.1 Curto-circuito na barra 1 sem abertura de linha 
 Nesta simulação, é considerado um curto-circuito trifásico próximo a barra 1 do sistema. 
Essa falta ocorre em t = 1s e cessa em t = 1,1s. Não há retirada da linha de falta, logo a matriz 
admitância após o periodo de falta permanece igual à do periodo atencessor a falta. A resposta 
dos ângulos do rotor é mostrada na figura 4.2. 
 
Figura 4.2 – Ângulos dos rotores em valores absolutos para falta na barra 1 sem o 
desligamento de linha com tempo de abertura em 2,4s. 
 
Percebe-se que os ângulos 𝛿 dos geradores 1 e 2 variam senoidalmente, contudo 
amortecidos devido as caracteristicas resistivas existentes na linha. Tal variação não caracteriza 
uma perda de sincronismo, já que se mantém dentro dos limites de estabilidade. 
Foram simulados tempos de abertura da falta até 2,5s (tempo máximo de simulação), quer 
dizer, a falta se manteve por mais tempo, e respostas similares a da Figura 4.2 foram obtidas. 
Pode-se dessa forma perceber que a situação em estudo não apresenta um tempo crítico de 
abertura. Ou seja, do ponto de vista da estabilidade transitória, qualquer tempo de abertura da 
54 
 
 
 
proteção garante o sincronismo entre as máquinas do sistema, uma vez que não há uma 
divergência entre os ângulos das máquinas. 
 
4.3.2 Curto-circuito na barra 2 com o desligamento da linha 2 – 4 
 Para esta simulação, é considerado um curto-circuito trifásico próximo a barra 2 do 
sistema. Essa falta ocorre em t = 1s e cessa em t = 1,1s. O defeito é eliminado com o 
desligamento da linha de transmissão que interliga as barras 2 e 4. Logo, para este caso a matriz 
de admitância no período pós-falta é diferente da matriz pré-falta. A resposta dos ângulos do 
rotor é mostrada na figura 4.3. 
 
Figura 4.3 – Ângulos dos rotores em valores absolutos para falta na barra 2 com o 
desligamento da linha 2 – 4. 
 
O comportamento da resposta