Proteção e estabilidade de sistemas elétricos_ Caio Graco

Prévia do material em texto

1 
PROTEÇÃO E ESTABILIDADE DE SISTEMAS ELÉTRICOS: ABORDAGEM POR 
FERRAMENTA MATEMÁTICA COMPUTACIONAL 1 
 
Caio Graco 2 
Eliel Roberto 3 
Judson Oliveira 4 
Mailton Ismar 5 
João Campos 6 
 
 
RESUMO 
 
O trabalho a seguir aborda de forma prática a aplicação das filosofias de 
proteção de sistemas elétricos de potência, com auxílio de ferramentas de simulação 
matemática foi possível comprovar efeitos e modificar parâmetros de projeto afim de 
resultados concisos. Dentro do contexto em estudo foi criado diagrama unifilar da 
estação trabalhada com contribuições de curtos-circuitos trifásicos, curto-circuito fase-
terra, além de abordar o dimensionamento e ajustes do transformador de corrente 
assim como dos relés de proteção 50/50N e 51/51N, demonstração do 
coordenograma de curvas de fase e neutro. É inerente ao próprio sistema elétrico de 
potência que ocorram falhas a rede, sendo assim impossível tornar o mesmo 
totalmente imune a essas ocorrências que por ventura serão minimizadas pelo estudo 
aqui em questão. Em um segundo momento deste trabalho é abordado o fluxo de 
potência do sistema elétrico, demonstrando a importância do mesmo no planejamento 
das condições de operação, controle e supervisão dos modelos de interesse de 
estudo. Esta abordagem matemática com ferramentas digitais é a proposta inicial do 
trabalho que segue. 
 
Palavras-chave: Relé de proteção. Fluxo de potência. Simulação. 
 
 
 
 
1Artigo apresentado à Universidade Potiguar - UnP, projeto integrante a disciplina de proteção e 
estabilidade de sistemas elétricos para o curso de Engenharia Elétrica. 
2Graduando em Engenharia Elétrica pela Universidade Potiguar – caio_gracolopes@hotmail.com 
3Graduando em Engenharia Elétrica pela Universidade Potiguar – eliel_sampa@yahoo.com.br 
4Graduando em Engenharia Elétrica pela Universidade Potiguar – judsonods@gmail.com 
5Graduando em Engenharia Elétrica pela Universidade Potiguar – ismareletrico@gmail.com 
6Prof. Dr. de engenharias e ciências exatas da Universidade Potiguar – joão.campos@unp.br 
2 
1 INTRODUÇÃO 
 
É fato que um dos parâmetros de desenvolvimento de um país é o seu sistema 
elétrico de potência que necessita ofertar aos seus usuários de forma múltipla e 
confiável um potencial energético que atenda aos requisitos básicos de economia e 
qualidade da energia, sendo intimamente ligados a esses requisitos a proteção e 
estabilidade do sistema como um todo. A proteção de sistemas elétricos atua com 
dois grandes objetivos; evitar falhas no sistema como curto circuitos que possam 
danificar equipamentos e materiais, promover o rápido reestabelecimento de energia, 
evitando danos aos consumidores e proporcionando uma qualidade no fornecimento 
da energia aos usuários (COTOSCK, 2007). 
Os estudos de proteção de sistemas elétricos devem assegurar aos usuários 
desses métodos a interrupção da eletricidade em casos de surtos na rede, visando 
proteger linhas, barramentos e equipamentos, mas além disso gerar dados com a 
intenção de posterior avaliação de causa das falhas ocorridas. Alguns aspectos são 
levados em consideração durante o projeto de proteção de sistemas, como custos de 
troca de relés eletromecânicos obsoletos por relés digitais, propagação do defeito 
gerado e tornar mínimo o tempo de inoperância do provimento de energia. 
A proteção é realizada por um conjunto de equipamentos, sendo os principais 
os transformadores de instrumentos, os relés e os disjuntores. Os transformadores de 
instrumentos, transformador de corrente e transformador de potencial, servem para 
traduzir para os níveis dos relés, a tensão e a corrente do sistema, respectivamente. 
Os relés servem para identificar o defeito no sistema, através da comparação do valor 
observado com o valor de seu ajuste, e acionar o disjuntor. O disjuntor serve para 
desconectar do sistema o elemento onde ocorreu o defeito, evitando dessa forma 
danos ao elemento e o funcionamento inadequado do sistema (SOUZA, 2010). 
O estudo adiante presente trata dos ajustes de projeto para relés de 
sobrecorrente 50/50N, 51/51N, dimensionamento do transformador de corrente, 
diagrama unifilar com curto circuitos, coordenação das curvas da fase e neutro. Em 
um segundo momento será apresentado o estudo de fluxo de carga de um sistema 
pré-definido pelo acadêmico Prof. Dr. João Campos. 
 
 
3 
2 DESENVOLVIMENTO 
 
O desenvolvimento do projeto se deu como um estudo de caso de uma indústria 
fictícia denominada “Industrial João Campos”, afim de modernização de seu sistema 
de proteção. Sendo assim os dados cedidos foram os seguintes: 
• Ponto de conexão COSERN 
o Z1 = 0,09 + j0,7 pu 
o Z0 = 0,13 + j1,6 pu 
o Sbase = 100 MVA (Potência de base) 
o Vbase = 13,8 kV (Tensão de base) 
o Transformador de corrente RTC 800/5 (Relé da COSERN) 
o Unidade 51: Tap 4,5A (secundário) curva 0,5 MI IEC 
o Unidade 50: Tap 14A (secundário) tempo 0,15s 
o Unidade 51N: Tap 0,45A (secundário) curva 0,17 NI IEC 
o Unidade 50N: Tap 1,5A (secundário) tempo 0,1s 
• Dados Indústria “João Campos”; 
o Trafo 2 MVA, 13.8 / 380v, Triângulo/Estrela aterrado, Z = 5% 
o Carga 1.5 MVA, FP=0.8 atrasado. 
o 2 linhas de transmissão com Z = 0.1 + j0.5 
o Gerador 0.2 MW, tensão de 1.0 pu 
Para um melhor entendimento do projeto, os cálculos e simulações 
matemáticas foram abordados passo a passo a seguir. 
 
2.1 CURTO-CIRCUITO 
 
Os curtos circuitos são alterações extremas da corrente que flui no sistema 
elétrico. Se não forem limitadas no seu modulo e no tempo, danificam os componentes 
elétricos por meio dos quais são conduzidos. Enquanto os tempos permitidos nos 
eventos de sobrecarga podem chegar a vários segundos, os tempos permitidos para 
a duração dos curtos-circuitos não devem superar o valor de 2 segundos. 
Normalmente, devem ser restritos entre 50 e 1000ms. Para tanto, os dispositivos de 
proteção devem ser extremamente velozes e os equipamentos de manobra, nos 
casos dos disjuntores e religadores, devem ter capacidade adequada para operar em 
condições extremas de corrente (MAMEDE, 2011). 
4 
2.1.1 Curto-circuito primário 
 
O código a seguir trata dos cálculos de curto-circuito aplicados a condição 
expressa no item 2.0 deste documento, as correntes abordadas são de curto-circuito 
no primário do transformador de 2,0 MVA. Este método calcula o curto-circuito 
simétrico no primário e o curto-circuito fase-terra máximo também no primário. 
 
% Cálculo das correntes de curto-circuito em instalações elétricas 
% Tensão de base primário do transformador 
Vbp = 13.8*10^3; 
% Tensão de base secundário do transformador 
Vbs = 380; 
% Potência de base estipulada pela COSERN 
Sb = 100*10^6; 
% Corrente de base do primário 
Ibp = Sb/(sqrt(3)*Vbp); 
% Corrente de base do secundário 
Ibs = Sb/(sqrt(3)*Vbs); 
% Impedância de base do primário 
Zbp = (Vbp^2)/Sb; 
% Impedância de base do secundário 
Zbs = (Vbs^2)/Sb; 
% Resistência de terra é apenas no secundário do transformador (R = 40m ohm) 
Z_terrapu = (40/Zbs); 
% Impedância do transformador em decimal 
Z_Transformador = 0.05; 
% Potência do transformador 
baseTransformador = 2*10^6; 
% Zpu do transformador 
Z_TransformadorPu = Z_Transformador*(Sb/baseTransformador)*i; 
% Sequência zero e positiva estipulada pela COSERN 
Z1_Distribuidora = 0.09 + 0.7*i; % pu em base 100 MVA 
Z0_Distribuidora = 0.13 + 1.6*i; % pu em base 100 MVA 
Zg1 = ((380)^2)/(0.2*10^6); % impedância do gerador 
Z1_geradorpu = Zg1/Zbs; % impedância do gerador em pu 
% A impedância já foi apresentada em pu 
ZL_pu = 0.1 + 0.5*i; % pu 
% Cálculo de curto-circuito 
fprintf('\n'); 
disp('Calculando curto-circuito primário do transformador' ); 
% z1_Equivalente = z2_equivalente 
fprintf('\n'); 
% Curto-circuito em retangular 
Icsp_retangular = (1/Z1_Distribuidora)*Ibp;% Curto-circuito simétrico no 
primário 
% Curto-circuito em polar 
Icsp = abs(Icsp_retangular); 
imp=imag(Icsp_retangular); 
rep= real(Icsp_retangular); 
Angulop= atand(imp/rep); 
fprintf('curto-circuito simétrico = '); 
fprintf('%d < %d',Icsp,Angulop); 
fprintf( ' A'); 
fprintf('\n'); 
fprintf('\n'); 
% curto-circuito F-N no primário 
5 
Iftp_retangular = (3/((2*Z1_Distribuidora)+Z0_Distribuidora))*Ibp; % 
Curto-circuito F-N no primário 
Iftp = abs(Iftp_retangular); 
imp2=imag(Iftp_retangular); 
rep2= real(Iftp_retangular); 
Angulop2= atand(imp2/rep2); 
fprintf('curto-circuito fase-terra = ' ); 
fprintf('%d < %d',Iftp,Angulop2); 
fprintf( ' A'); 
fprintf('\n'); 
fprintf('\n'); 
 
Resultado da execução do programa: 
 
Figura 1 - Curto-circuito no primário do transformador de potência 
Fonte: Autor 
 
 
2.1.2 Curto-circuito secundário 
 
Este código calcula o curto-circuito simétrico, fase-terra máximo e fase-terra 
mínimo no secundário do transformador. Com o relé posicionado a montante do 
transformador, ou seja, no primário do mesmo que possibilita a utilização dos valores 
calculados de curto-circuito do primário para a parametrização do relé de proteção 
com as funções 50/51 e 50/51N. sendo os demais dados utilizados para melhor 
entendimento do sistema elétrico de potência. 
 
% Curto-circuito simétrico, fase-terra para a pior situação, considerando Z1 = 
Z0 da linha de energização. 
disp('Calculando curto-circuito secundário do Transformador' ); % 
Barra da subestação 
fprintf('\n'); 
% ZTH = Z1_Distribuidora + Z_TransformadorPu. 
Icssc1 = (1/(Z1_Distribuidora+Z_TransformadorPu))*Ibs; 
Icssc2 = (1/(Z1_geradorpu+ZL_pu))*Ibs; % Pior situação possível do sistema 
% Curto-circuito em retangular 
Icssc_retangular = Icssc1+Icssc2; 
% Curto-circuito em polar 
Icssc = abs(Icssc_retangular); 
ims= imag(Icssc_retangular); 
res= real(Icssc_retangular); 
Angulos = atand(ims/res); 
fprintf('curto-circuito simétrico = '); 
fprintf('%d < %d',Icssc,Angulos); 
fprintf( ' A'); 
% Ia = 1/(z1_Equivalente+z2_Equivalente+z0_Equivalente+3zg) 
6 
fprintf('\n'); 
fprintf('\n'); 
% (Z1H,Z2H) Z0H 
% O pior caso possível 
Iftsc1= (3/((2*(Z1_Distribuidora+Z_TransformadorPu))+(Z_TransformadorPu)))*Ibs; 
% 
Fase-terra máximo 
Iftsc2 = (3/(3*(Z1_geradorpu)+ZL_pu))*Ibs; % fase-terra máximo 
% Curto-circuito em retangular 
Iftsc_retangular= Iftsc1 + Iftsc2; 
% Curto-circuito em polar 
Iftsc = abs(Iftsc_retangular); 
ims2= imag(Iftsc_retangular); 
res2= real(Iftsc_retangular); 
Angulos2= atand(ims2/res2); 
fprintf('curto-circuito fase-terra (máx.) = ' ); 
fprintf('%d < %d',Iftsc,Angulos2); 
fprintf( ' A'); 
fprintf('\n'); 
fprintf('\n'); 
% Curto-circuito fase terra mínimo (pior situação possível) 
% (Z1H,Z2H) Z0H 
ZF 
% Para o gerador a z1 = z0 = z2 
Iftsc1m = 
(3/((2*(Z1_Distribuidora+Z_TransformadorPu))+(Z_TransformadorPu+Z_terrapu) 
+3*(Z_terrapu)))*Ibs; % Fase-terra mínimo 
Iftsc2m = (3/(3*(Z1_geradorpu)+ZL_pu))*Ibs; % Fase-terra mínimo 
% Curto-circuito em retangular 
Iftsc_retangularm= Iftsc1m + Iftsc2m; 
% Curto-circuito em polar 
Iftscm = abs(Iftsc_retangularm); 
ims2m = imag(Iftsc_retangularm); 
res2m = real(Iftsc_retangularm); 
Angulos2m= atand(ims2m/res2m); 
fprintf('curto-circuto fase-terra (min) = ' ); 
fprintf('%d < %d',Iftscm,Angulos2m); 
fprintf( ' A'); 
fprintf('\n'); 
fprintf('\n'); 
 
% Obs: é de suma importância salientar que Z0 > Z1, então o curto-circuito fase- 
terras será maior que o simétrico 
 
Resultado da execução do programa: 
 
Figura 2 - Curto-circuito no secundário do transformador de potência 
Fonte: Autor 
 
 
 
7 
2.2 RELÉ DE PROTEÇÃO (UNIDADE TEMPORIZADA DE FASE) 
 
São relés que operam quando o valor da corrente do circuito ultrapassa um 
valor pré-fixado ou ajustado. Os relés de sobrecorrente podem ser instantâneos (ANSI 
50) ou temporizados (ANSI 51), a característica dos relés de sobrecorrente é 
representada pelas suas curvas tempo versus corrente. Estas curvas variam em 
função do tipo de relé (MARDEGAM, 2010). 
Os tipos de curva abordadas neste estudo de proteção e seletividade das 
unidades temporizadas do relé foram as curvas muito inversa tanto para a unidade 
consumidora quanto para a unidade geradora, ambos apresentaram o mesmo TMS. 
O critério utilizado para a seletividade foi o cálculo do múltiplo consumidor e múltiplo 
da unidade geradora, sendo que esses múltiplos são responsáveis pela seletividade 
dos dois relés, tanto da concessionária quanto da indústria João Campos. 
Em relação do dimensionamento do transformador de corrente, é necessário 
que o mesmo atenda a condição de no mínimo 20 vezes o valor da corrente nominal 
da barra. 
 
=========================================================================== 
% ==== Estabilidade de Sistemas de Potência================================ 
% ==== Nome: Caio,Eliel,Judson e Mailton ================================== 
% ==== Data: 20/11/2018 =================================================== 
% ==== Orientador: Dr.João Campos ========================================= 
% ==== Parametrização do Relé 50/51 P ===================================== 
% ==== A Cosern forneceu os seguintes dados do ponto de conexão =========== 
% ==== • Z1= 0,09 +j 0,7 pu; ============================================== 
% ==== • Z0 = 0,13 +j 1,6 pu; ============================================= 
% ==== • Sbase: 100 MVA; Vbase: 13,8 kV =================================== 
% ==== Transformador de corrente associado ao relé da Cosern ============== 
% ==== RTC 800/5 ========================================================== 
% ==== Fase: ============================================================== 
% ==== Unidade 51: Tap 4,5 A (secundário) curva 0,5 MI IEC================= 
% ==== Unidade 50: Tap 14 A (secundário) Tempo: 0,15 s ==================== 
% ========================================================================= 
% ========================================================================= 
 
fprintf('\n\n\n'); 
fprintf('\n Corrente nominal no primário do transformador de potência:\n'); 
baseTransformador = 2*10^6; % 2,0 MVA 
Vbp = 13800; % Tensão de base do primário 
InP = (baseTransformador/(sqrt(3)*Vbp)); % Corrente nominal do trafo 
fprintf(' In = %.2f A \n',InP); % Corrente nominal 
fprintf('\n Cálculo do transformador de corrente:\n'); 
Icssc = 5927.92; % Corrente de curto circuito no barramento da COSERN 
In_Tc = Icssc/20; % Estimativa do transformador de corrente 
fprintf(' In_Tc = %.2f A \n',In_Tc); % Corrente nominal do transformador de 
corrente 
% É aconselhável utilizar o tc com corrente nominal de 300A, e RTC = 60 
fprintf('\n Unidade temporizada de fase (F51P) do consumidor:\n'); 
8 
% Sabendo que a corrente nominal da subestação do consumidor é de 83.67 A, 
em 13,8 kV, e considerando um fator de crescimento de carga igual a 1,2 (1,2 
- 1,5) tem-se: 
Kf_consumidor = 1.2; 
RTC = 60; 
Tap_consumidor51 = (Kf_consumidor*InP)/RTC; % tap do consumidor 
fprintf('\n Secundário:\n'); 
fprintf(' Tap_consumidor = %.2f A \n',Tap_consumidor51 );% Corrente nominal 
do transformador de corrente 
fprintf('\n Primário:\n'); 
fprintf(' Tap_consumidor = %.2f A \n',(Tap_consumidor51*RTC) ); % Corrente 
nominal do transformador de corrente para calcular os tempos adequados de 
atuação do relé do consumidor, é necessário realizar os cálculos dos múltiplos 
(m): 
% Múltiplo da concessionária 
fprintf('\n Múltiplo da COSERN - fase (51P):\n'); 
RTC_cosern = 800; % RTC abordado no problema 
Tap_cosern51 = 4.5; % Tap abordado no problema 
M_cosern = (Icssc/(RTC_cosern*(Tap_cosern51))); 
fprintf(' M_cosern = %.2f \n',M_cosern); % Múltiplo da COSERN 
% Múltiplo do consumidor 
fprintf('\n Múltiplo do consumidor - fase (51P):\n'); 
M_consumidor = (Icssc/(RTC*(Tap_consumidor51)));fprintf(' M_consumidor = %.2f \n',M_consumidor); % Múltiplo do consumidor 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Cálculos dos tempos %%%%%%%%%%%%%%%%%%%%%%%%%% 
 
fprintf('\n Tempo do Relé da COSERN (s):\n'); 
% Curva adotada é a curva muito inversa -> IEC 
Tms_cosern = 0.5; 
T_cosern = (13.5/((M_cosern^1)-1))*Tms_cosern; % Tempo em segundos 
fprintf(' Tempo = %.2f s \n',((T_cosern))); % Tempo em ms 
% Curva adotada é a curva muito inversa -> IEC 
fprintf('\n Tempo do Relé do consumidor (s):\n'); 
 
Tms_consumidor = 0.5; 
T_consumidor = (13.5/((M_consumidor^1)-1))*Tms_consumidor; % Tempo em 
segundos 
fprintf(' Tempo = %.2f s \n',((T_consumidor)));% Tempo em ms 
% Dessa forma, a diferença de tempo entre a atuação do relé da concessionária 
e do consumidor é: 
fprintf('\n Tempo de seletividade:\n'); 
delta_tempo = T_cosern-T_consumidor; 
fprintf(' Tempo = %.2f s \n',((delta_tempo))); 
 
%%%%%%%%%%%%%%%%%%%%%%% Função 50 – fase %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% Como a unidade instantânea de fase (50) não deve atuar para as correntes 
de energização dos transformadores (inrush), tem-se que: 
fprintf('\n Função 50 (secundário do TC) - fase:\n'); 
Kf_consumidor50 = 8; % Considerar curto-circuito = 8*In 
RTC50 = 60; 
Tap_consumidor50 = (Kf_consumidor50*InP)/RTC50; 
fprintf(' Tap_consumidor50 = %.2f A \n',Tap_consumidor50); 
% Sabendo que o tap do consumidor no secundário do TC é 5A, então no primário 
do transformador de corrente a corrente máxima de curto-circuito é de 
aproximadamente 300A 
fprintf('\n Tempo do consumidor função 50 - fase:\n'); 
Tms_consumidor = 0.5; 
M_consumidor50 = Icssc/(Tap_consumidor50*RTC50); 
9 
T_consumidor50 = (13.5/((M_consumidor50^1)-1))*Tms_consumidor; % Tempo em 
segundos 
fprintf(' Tempo = %.2f s \n',((T_consumidor50)));% Tempo em ms 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Consumidor %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
TMS1 = 0.5; 
Iaj1 = 1.67; % Corrente de ajuste do relé da função 51P 
I1 = (2.07:0.01:11.16) ; % Vetor intervalo no secundário 
for ii = 1: length(I1) % Percorre o vetor de acionamento das funções 
t1(ii) = (13.5*TMS1) / (((I1(ii)/Iaj1).^1) -1); % Vetor tempo de atuação 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Cosern %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
TMS2 = 0.5; 
Iaj2 = 4.5; % Corrente de ajuste do relé da função 51P 
I2 = (4.9:0.01:14); % Vetor intervalo no secundário 
for i2 = 1: length(I2) % Percorre o vetor de acionamento das funções 
t2(i2) = (13.5*TMS2) / (((I2(i2)/Iaj2).^1) -1); % Vetor tempo de atuação 
end 
% Plote a curva de corrente vs tempo 
figure(1); 
plot(I1,t1,'Color','b','LineWidth',5.0); 
hold on; 
plot(I2,t2,'Color','r','LineWidth',5.0); 
xlabel('Corrente elétrica(A)'); 
ylabel('Tempo de atuação(s)'); 
title ('Curva T X I Relé 51F'); 
legend ('Indústria = 0,5','Cosern = 0,5'); 
grid on; 
hold off; 
figure(2); 
plot(I1,t1,'Color','b','LineWidth',5.0); 
xlabel('Corrente elétrica(A)'); 
ylabel('Tempo de atuação(s)'); 
title ('Curva T X I Relé 51F'); 
legend ('Indústria = 0,5'); 
grid on; 
figure(3); 
plot(I2,t2,'Color','r','LineWidth',5.0); 
xlabel('Corrente elétrica(A)'); 
ylabel('Tempo de atuação(s)'); 
title ('Curva T X I Relé 51F'); 
legend ('Cosern = 0,5'); 
grid on; 
 
% Obs: foi utilizado nas duas situações a mesma configuração de curva, pois 
ao calcular a seletividade dos relés utilizando os múltiplos, é possível 
perceber que mesmo com a mesma curva e Tms o relé da indústria irá respeitar 
a seletividade do sistema elétrico de potência (tempo > 400ms). É importante 
salientar que esses gráficos estão sendo trabalhados com níveis de corrente 
adequados ao secundário do transformador de corrente. 
 
 
 
10 
Resultado da execução do programa: 
 
 
 
 
Figura 3 - Parametrização para função 50/51P (relé de proteção) 
Fonte: Autor 
 
 
2.3 RELÉ DE PROTEÇÃO (UNIDADE TEMPORIZADA DE NEUTRO) 
 
O código abaixo aplica-se a limitação de curto-circuito monofásico na barra de 
entrada situada a jusante do transformador principal, sendo que algumas situações o 
curto-circuito fase-terra será maior que o curto-circuito simétrico, porém, neste caso o 
curto-circuito fase-terra apresenta menor intensidade. 
Em relação do dimensionamento do transformador de corrente, é necessário 
que o mesmo atenda a condição de no mínimo 20 vezes o valor da corrente nominal 
da barra. 
11 
 
% =========== Estabilidade de Sistemas de Potência ======================== 
% =========== Nome: Caio,Eliel,Judson e Mailton =========================== 
% =========== Data: 20/11/2018 ============================================ 
% =========== Orientador: Dr.João Campos ================================== 
% =========== Parametrização do Relé 50/51 P ============================== 
% =========== A Cosern forneceu os seguintes dados do ponto de conexão ==== 
% • Z 1= 0,09 +j 0,7 pu;=================================================== 
% • Z0 = 0,13 +j 1,6 pu;=================================================== 
% • Sbase: 100 MVA; Vbase: 13,8 kV ======================================== 
% Transformador de corrente associado ao relé da Cosern =================== 
% RTC 800/5 =============================================================== 
% Fase: =================================================================== 
% Unidade 51N: Tap 0,45 A (secundário) curva 0,17 NI IEC ================== 
% Unidade 50N: Tap 1,5 A (secundário) Tempo: 0,1 s ======================== 
% Sabendo que a corrente nominal da subestação do consumidor é de 83.67 A, 
em 13,8 kV e considerando um fator de desequilíbrio de 20%, tem-se que as 
seguintes condições devem ser satisfeitas: 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Função 50/51 N %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
fprintf('\n\n\n'); 
fprintf('\n Corrente nominal no primário do transformador de Potência:\n'); 
baseTransformador = 2*10^6; % 2,0 MVA 
Vbp = 13800; % Tensão de base do primário 
InP = (baseTransformador/(sqrt(3)*Vbp)); % Corrente nominal do trafo 
fprintf(' In = %.2f A \n',InP); % Corrente nominal 
fprintf('\n Cálculo do Transformador de corrente:\n'); 
Iftscmax = 4161.54; % Corrente de curto circuito no barramento da cosern 
In_Tc = Iftscmax/20; % Estimativa do transformador de corrente 
fprintf(' In_Tc = %.2f A \n',In_Tc); % Corrente nominal do TC 
% É aconselhavel utilizar o tc com corrente nominal de 250A, e RTC = 50 
fprintf('\n Unidade temporizada de fase (F51N) do consumidor:\n'); 
% Sabendo que a corrente nominal da subestação do consumidor é de 83.67 A, 
em 13,8 kV, e considerando um fator de crescimento de carga igual a 0,2 (0,1 
- 0,45) tem-se: 
Kf_consumidorN = 0.2; % considera o fator de serviço de 20% 
RTCN = 50; 
Tap_consumidor51N = (Kf_consumidorN*InP)/RTCN; % tap do consumidor 
fprintf('\n Segundario:\n'); 
fprintf(' Tap_consumidor = %.2f A \n',Tap_consumidor51N );% Corrente nominal 
do transformador de corrente 
fprintf('\n Primario:\n'); 
fprintf(' Tap_consumidor = %.2f A \n',(Tap_consumidor51N*RTCN) ); % Corrente 
nominal do transformador de corrente 
% Para calcular os tempos adequados de atuação do relé do consumidor, é 
necessário realizar os cálculos dos múltiplos (m): 
% Múltiplo da concessionária 
fprintf('\n Múltiplo da cosern - Neutro (51N):\n'); 
RTC_cosernN = 800; % RTC abordado no problema 
Tap_cosern51N = 0.45; % tap abordado no problema 
M_cosernN = (Iftscmax/(RTC_cosernN*(Tap_cosern51N))); 
fprintf(' M_cosern = %.2f \n',M_cosernN); % Múltiplo da cosern 
% Múltiplo do consumidor 
fprintf('\n Múltiplo do consumidor - Neutro (51N):\n'); 
M_consumidorN= (Iftscmax/(RTCN*(Tap_consumidor51N))); 
fprintf(' M_consumidor = %.2f \n',M_consumidorN); % Múltiplo do consumidor 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Cálculos dos tempos %%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
12 
fprintf('\n Tempo do Relé da cosern 51N (s):\n'); 
% curva adotada é a curva normalmente inversa -> IEC 
Tms_cosernN = 0.17; 
T_cosernN = (0.14/((M_cosernN^0.02)-1))*Tms_cosernN; % tempo em segundos 
fprintf(' Tempo = %.2f s \n',((T_cosernN))); % tempo em ms 
% curva adotada é a curva muito inversa -> IEC 
fprintf('\n Tempo do Relé do consumidor 51N (s):\n'); 
Tms_consumidorN = 0.17; 
T_consumidorN = (13.5/((M_consumidorN^1)-1))*Tms_consumidorN; % tempo em 
segundos 
fprintf(' Tempo = %.2f s \n',((T_consumidorN)));% tempo em ms 
% Dessa forma, a diferença de tempo entre a atuação do relé da concessionária 
e do consumidor é: 
fprintf('\n Tempo de seletividade 51N:\n'); 
delta_tempoN = T_cosernN-T_consumidorN; 
fprintf(' Tempo = %.2f s \n',((delta_tempoN))); 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%% Função 50 - Neutro %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Como a unidade instantânea de neutro (50) não deve atuar para as correntes 
de energização dos transformadores (inrush), tem-se que: 
 
fprintf('\n Função 50 (secundário do TC) - Neutro:\n'); 
Kf_consumidor50N = 8; % considerar curto-circuito = 8*In 
RTC50N = 50; 
Tap_consumidor50N = (Kf_consumidor50N*InP)/RTC50N; 
fprintf(' Tap_consumidor50 = %.2f A \n',Tap_consumidor50N); 
% sabendo que o tap do consumidor no secundário do TC é 13.39 A, então no 
primário do transformador de corrente a corrente máxima de curto-circuito é 
de aproximadamente 107,12 A 
fprintf('\n Tempo do consumidor função 50 - Neutro:\n'); 
% Curva muito inversa -> IEC 
Tms_consumidorN = 0.17; 
M_consumidor50N = Iftscmax/(Tap_consumidor50N*RTC50N); 
T_consumidor50N = (13.5/((M_consumidor50N^1)-1))*Tms_consumidorN; % tempo em 
segundos 
fprintf(' Tempo = %.2f s \n',((T_consumidor50N)));% tempo em ms 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%% Consumidor %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% Curva muito inversa -> IEC 
TMS1N = 0.17; 
Iaj1N = 0.33; % corrente de ajuste do relé da função 51P 
I1N = (0.43:0.0001:13.39) ; % vetor intervalo no secundário 
for ii = 1: length(I1N) % percorre o vetor de acionamento das funções 
t1N(ii) = (13.5*TMS1N) / (((I1N(ii)/Iaj1N).^1)-1); % vetor tempo de atuação 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Cosern %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% Normalmente inversa -> IEC 
TMS2N = 0.17; 
Iaj2N = 0.45; % corrente de ajuste do relé da função 51P 
I2N = (0.49:0.0001:1.5) ; % vetor intervalo no secundário 
for i2 = 1: length(I2N) % percorre o vetor de acionamento das funções 
t2N(i2) = (0.14*TMS2N) / (((I2N(i2)/Iaj2N).^0.02)-1); % vetor tempo de 
atuação 
end 
13 
% Plote a curva de corrente vs tempo 
figure(1); 
plot(I1N,t1N,'Color','b','LineWidth',5.0); 
hold on; 
plot(I2N,t2N,'Color','r','LineWidth',5.0); 
xlabel('Corrente elétrica(A)'); 
ylabel('tempo de atuação(s)'); 
title ('Curva T X I Relé 51N'); 
legend ('Indústria = 0,17','Cosern = 0,17'); 
grid on; 
hold off; 
figure(2); 
plot(I1N,t1N,'Color','b','LineWidth',5.0); 
xlabel('Corrente elétrica(A)'); 
ylabel('tempo de atuação(s)'); 
title ('Curva T X I Relé 51N'); 
legend ('Indústria = 0,17'); 
grid on; 
figure(3); 
plot(I2N,t2N,'Color','r','LineWidth',5.0); 
xlabel('Corrente elétrica(A)'); 
ylabel('tempo de atuação(s)'); 
title ('Curva T X I Relé 51N'); 
legend ('Cosern = 0,17'); 
grid on; 
 
% obs.: foi utilizado nas duas situações a mesma configuração de curva, pois 
ao calcular a seletividade dos relés utilizando os múltiplos, é possível 
perceber que com mesmo Tms e curvas diferente o relé da indústria irá 
respeitar a seletividade do sistema elétrico de potência (tempo> 400ms). É 
importante salientar que esses gráficos estão sendo trabalhados com níveis 
de corrente adequados ao secundário do transformador de corrente. 
 
Resultado da execução do programa: 
 
 
14 
 
 
Figura 4 - Parametrização para função 50/51N (relé de proteção) 
Fonte: Autor 
 
 
 
3 FLUXO DE CARGA 
 
O problema de fluxo de carga é calcular a magnitude de tensão e ângulo de 
fase em cada barramento da rede elétrica de potência em condições de estado. Os 
estudos de fluxo de carga normalmente chamados estudos de fluxo de potência, são 
extremamente importantes para o desenho, planejamento e controle de sistemas 
elétricos de potência, partindo do diagrama unifilar e dos dados de entrada dos 
barramentos, linhas de transmissão e transformadores (GALVEZ, 2009). 
 
3.1 ADMITÂNCIA 
 
 A determinação da matriz admitância nodal (Y) tem grande importância para os 
cálculos de rede elétrica em sistemas de potência, a matriz admitância relaciona as 
tensões elétricas nodais com as correntes elétricas injetadas ao sistema através de 
geradores (BENEDITO, 2010). 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
============= Estabilidade de Sistemas de Potência ======================== 
============= Nome: Caio,Eliel,Judson e Mailton =========================== 
============= Data: 15/11/2018 ============================================ 
============= Orientador: Dr.João Campos ================================== 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% A impedância da linha em Pu foi proposta na questão 
% ZLpu = 0.1000 + 0.5000*i 
% Digite as linhas de transmissão propostas no sistema 
% | Da | Para | R | X | Bshunt | 
% | Barra | Barra | | | | 
% Matriz dados linhas 
d_l = [ 1 2 0.0000 2.5000 0.0000; 
 2 3 0.1000 0.5000 0.0000; 
 2 4 0.1000 0.5000 0.0000; ]; 
% Monta Matriz admitância de barras Y 
15 
fb = d_l(:,1); % Da barra numero 
tb = d_l(:,2); % Para a barra numero 
r = d_l(:,3); % Resistencia, R 
x = d_l(:,4); % Reatância, X 
b = d_l(:,5); % Admitância B Shunt 
z = r + i*x; % Matrix Z 
y = 1./z; % Inverte cada elemento de Z 
b = i*b; % Muda Bshunt para imaginário 
nbarra = max(max(fb),max(tb)); % identifica número de barras 
nramos = length(fb); % número de ramos 
Ybarra = zeros(nbarra,nbarra); % Inicializa YBarra 
% Completa Elementos fora da diagonal 
for k=1:nramos 
 Ybarra(fb(k),tb(k)) = -y(k); 
 Ybarra(tb(k),fb(k)) = Ybarra(fb(k),tb(k)); 
end 
% Completa os elementos da Diagonal 
for m=1:nbarra 
 for n=1:nramos 
 if fb(n) == m | tb(n) == m 
 Ybarra(m,m) = Ybarra(m,m) + y(n) + b(n); 
 end 
 end 
end 
fprintf('\n\n'); 
fprintf('\n MATRIZ ADMITÂNCIA DE BARRA: \n\n'); 
disp(Ybarra); 
fprintf('\n\n'); 
 
Resultado da execução do programa: 
 
Figura 5 - Matriz admitância de barra (sistema de potência proposto) 
Fonte: Autor 
 
3.2 FLUXO DE CARGA 
 
O código abaixo calcula o fluxo de carga do sistema proposto, em anexo estão 
os resultados por simulação matemática em Octave e simulação em software de 
sistema elétrico de potência Power World. 
=========================================================================== 
=========== Estabilidade de Sistemas de Potência ========================== 
=========== Nome: Caio,Eliel,Judson e Mailton ============================= 
=========== Data: 15/11/2018 ============================================== 
=========== Orientador: Dr.João Campos ====================================%% digite o percurso das barras “DE(Barra) PARA(Barra)” 
d_l = [ 1 2 ; 
16 
 2 3 ; 
 2 4 ; ]; 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%% matriz admitância de barra %%%%%%%%%%%%%%%%%%%%%% 
Ybarra= [ 
 
 0.00000 - 0.40000*i -0.00000 + 0.40000*i 0.00000 + 0.00000*i 
0.00000 + 0.00000*i; 
 -0.00000 + 0.40000*i 0.76923 - 4.24615*i -0.38462 + 1.92308*i 
-0.38462 + 1.92308*i; 
 0.00000 + 0.00000*i -0.38462 + 1.92308*i 0.38462 - 1.92308*i 
0.00000 + 0.00000*i; 
 0.00000 + 0.00000*i -0.38462 + 1.92308*i 0.00000 + 0.00000*i 
0.38462 - 1.92308*i; ]; 
fb = d_l(:,1); 
tb = d_l(:,2); 
nbarra = max(max(fb),max(tb)); % identifica número de barras 
nramos = length(fb); % número de ramos 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Dados fornecidos das barras 
% Tipos: 1- Referencia; 2-PV(tensão controlada); 3-PQ(barra de carga) 
 
% | Barra | Tipo | V | delta | PGi | QGi | PDi | 
QDi | 
 
dados_barra = [ 
 1 1 1.0000 0.0000 0.0000 0.0000 0.0000 
0.0000 ; 
 2 3 1.0000 0.0000 0.0000 0.0000 0.0000 
0.0000 ; 
 3 3 1.0000 0.0000 0.0000 0.0000 0.0120 
0.0090 ; 
 4 2 1.0000 0.0000 0.0020 0.0000 0.0000 
0.0000 ; 
 ]; 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% Assume-se que a Barra 1 é a barra de referência ou "Slack bus". 
% Leitura e interpretação dos dados fornecidos, separando em vetores 
 
barra = dados_barra(:,1); % Número da barra. 
tipo = dados_barra(:,2); % Tipo da barra 1-Slack, 2-PV, 3-PQ. 
 V = dados_barra(:,3); % Tensões iniciais nas barras 
 del = dados_barra(:,4); % Ângulos iniciais nas barras 
 Pg = dados_barra(:,5); % PGi, Potência Ativa injetada na barra 
 Qg = dados_barra(:,6); % QGi, Potência Reativa injetada na barra 
 Pl = dados_barra(:,7); % PDi, Potência Ativa demandada na barra 
 Ql = dados_barra(:,8); % QDi, Potência Reativa demanda na barra 
 P = Pg - Pl; % Pi = PGi - PDi, Potencia Ativa liquida na 
barra 
 Q = Qg - Ql; % Qi = QGi - QDi, Potencia Reativa liquida na 
barra 
 Psp = P; % P especificadas 
 Qsp = Q; % Q especificadas 
 G = real(Ybarra); % Matriz Condutância 
17 
 B = imag(Ybarra); % Matriz Susceptância 
 
 pv = find(tipo==2 | tipo==1); % Barras PV 
 pq = find(tipo == 3); % Barras PQ 
 npv = length(pv); % Numero de barras PV 
 npq = length(pq); % Numero de barras PQ 
k_max = 3; 
 
%% Para atingir o número de interações desejado foi arbitrado a tolerância 
de erro de 0,01 utilizando o método Newton Raphson 
 
for iteracao = 1:k_max 
 
 P = zeros(nbarra,1); 
 Q = zeros(nbarra,1); 
 
 % Calcula P e Q 
 
 for i = 1:nbarra 
 for k = 1:nbarra 
 P(i) = P(i) + V(i)* V(k)*(G(i,k)*cos(del(i)-del(k)) + B(i,k)*sin(del(i)-
del(k))); 
 Q(i) = Q(i) + V(i)* V(k)*(G(i,k)*sin(del(i)-del(k)) - B(i,k)*cos(del(i)-
del(k))); 
 end 
 end 
 % Calcula a mudança para os valores especificados 
 dPa = Psp-P; 
 dQa = Qsp-Q; 
 k = 1; 
 dQ = zeros(npq,1); 
 for i = 1:nbarra 
 if tipo(i) == 3 
 dQ(k,1) = dQa(i); 
 k = k+1; 
 end 
 end 
 dP = dPa(2:nbarra); 
 M = [dP; dQ]; % Vetor de diferenças/variações de potencias 
 % Montagem da Matriz Jacobiana 
 % Monta a Matriz J1 - Derivada parcial da Potência Ativa pelos Ângulos 
 J1 = zeros(nbarra-1,nbarra-1); 
 for i = 1:(nbarra-1) 
 m = i+1; 
 for k = 1:(nbarra-1) 
 n = k+1; 
 if n == m 
 for n = 1:nbarra 
 J1(i,k) = J1(i,k) + V(m)* V(n)*(-G(m,n)*sin(del(m)-del(n)) + 
B(m,n)*cos(del(m)-del(n))); 
 end 
 J1(i,k) = J1(i,k) - V(m)^2*B(m,m); 
 else 
 J1(i,k) = V(m)* V(n)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-
del(n))); 
 end 
 end 
end 
 
% Monta Matriz J2 - Derivada parcial da Potência Ativa pela tensão 
18 
 
 J2 = zeros(nbarra-1,npq); 
 for i = 1:(nbarra-1) 
 m = i+1; 
 for k = 1:npq 
 n = pq(k); 
 if n == m 
 for n = 1:nbarra 
 J2(i,k) = J2(i,k) + V(n)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-
del(n))); 
 end 
 J2(i,k) = J2(i,k) + V(m)*G(m,m); 
 else 
 J2(i,k) = V(m)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n))); 
 end 
 end 
end 
 
 
% Monta Matriz J3 - Derivada parcial da Potência Reativa pelos Ângulos 
 
 J3 = zeros(npq,nbarra-1); 
 for i = 1:npq 
 m = pq(i); 
 for k = 1:(nbarra-1) 
 n = k+1; 
 if n == m 
 for n = 1:nbarra 
 J3(i,k) = J3(i,k) + V(m)* V(n)*(G(m,n)*cos(del(m)-del(n)) + 
B(m,n)*sin(del(m)-del(n))); 
 end 
 J3(i,k) = J3(i,k) - V(m)^2*G(m,m); 
 else 
 J3(i,k) = V(m)* V(n)*(-G(m,n)*cos(del(m)-del(n)) - B(m,n)*sin(del(m)-
del(n))); 
 end 
 end 
end 
 
 
% Monta matriz J4 - Derivada parcial da Potência Reativa pela Tensão 
 
 J4 = zeros(npq,npq); 
 for i = 1:npq 
 m = pq(i); 
 for k = 1:npq 
 n = pq(k); 
 if n == m 
 for n = 1:nbarra 
 J4(i,k) = J4(i,k) + V(n)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-
del(n))); 
 end 
 J4(i,k) = J4(i,k) - V(m)*B(m,m); 
 else 
 J4(i,k) = V(m)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)- del(n))); 
 end 
 end 
end 
 
 J = [J1 J2; J3 J4]; % Matriz Jacobiana 
19 
 X = inv(J)*M; % Vetor de correção dos valores 
 
 dTh = X(1:nbarra-1); % Modifica Angulo da Tensão 
 dV = X(nbarra:end); % Modifica Modulo da Tensão 
 
 del(2:nbarra) = dTh + del(2:nbarra); % Angulo da tensão 
 
 k = 1; 
 for i = 2:nbarra 
 if tipo(i) == 3 
 V(i) = dV(k) + V(i); % Modulo da tensão na barra da carga 
 k = k+1; 
 end 
 end 
 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
fprintf('\n Método: Newton-Raphson \n'); 
fprintf('\n NUMERO DE ITERACÕES REALIZADAS: %d \n',iteracao); 
fprintf('\n RESULTADO DOS CALCULOS DAS TENSOES, ANGULOS NAS BARRAS\n'); 
fprintf('\n'); 
 
for n=1:nbarra 
 
 fprintf(' V%d = %.4f < %.3f [deg]\n',n,V(n),((180/pi)*del(n))); 
 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fluxos %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
fprintf('\n\n FLUXOS DE POTÊNCIA NAS LINHAS:'); 
Vm = V.*cos(del) + j*V.*sin(del); 
ri=1; 
for rf=1:nramos 
% fluxo da barra d_l(ri,1) para a barra d_l(rf,2) 
 S(d_l(ri,1),d_l(rf,2))=Vm(d_l(ri,1))*((conj(Vm(d_l(ri,1)))-
conj(Vm(d_l(rf,2))))/Ybarra(d_l(ri,1),d_l(rf,2))^-1); 
 if S(d_l(ri,1),d_l(rf,2))~=0 
 fprintf('\n\n P%d%d=%.4f[pu] Q%d%d=%.4f[pu] 
',d_l(ri,1),d_l(rf,2),real(S(d_l(ri,1),d_l(rf,2))),d_l(ri,1),d_l(rf,2),imag
(S(d_l(ri,1),d_l(rf,2)))); 
 fprintf(' P%d%d=%.4f[pu] Q%d%d=%.4f[pu]\n',d_l(rf,2),d_l(ri,1),-
real(S(d_l(ri,1),d_l(rf,2))),d_l(rf,2),d_l(ri,1),imag(S(d_l(ri,1),d_l(rf,2)
))); 
 end 
 ri=ri+1; 
 if ri>nramos 
 break; 
 end 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
20 
4 CONCLUSÃO 
 
Este projeto teve oobjetivo de apresentar a implantação via método 
matemático de Newton Raphson em um sistema elétrico de potência composto por 4 
barras, utilizando apenas as 3 primeiras iterações do sistema, sendo esse método 
amplamente usado na prática por ser mais eficiente e convergir mais rapidamente aos 
resultados. Além disso, foi implementado nesse memorial de cálculo algoritmos de 
parametrização para dispositivo digital de proteção (Relé), tendo como ênfase as 
funções 50/51P e 50/51N. 
A interface do código foi comentada passo a passo para evitar duvidas nas 
simulações, os resultados apresentados foram satisfatórios quando comparados aos 
resultados obtidos via software Power World (fluxo de carga). Já para a 
parametrização do relé de proteção foi utilizado apenas o Octave para execução dos 
cálculos e gerar coordenogramas. Segue em anexo os resultados obtidos pela 
simulação no Power World, resolução do fluxo de carga com Octave e manuscrito, 
tabela geral de resultados e os coordenogramas proposto pelo o contratante. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
21 
PROTECTION AND STABILITY OF ELECTRICAL SYSTEMS: approach by 
computational mathematical tool 
 
ABSTRACT 
 
The following work deals in a practical way with the application of power protection 
systems, with the aid of mathematical simulation tools, it was possible to verify effects 
and modify design parameters to obtain concise results. Within the universe under 
study was created a single-line diagram of the station worked with contributions of 
three-phase short circuits, phase-to-ground short circuit, besides addressing the 
current transformer design, dimensioning and adjustments of 50 / 50N and 51 / 51N 
protection relays, coordinate of phase and neutral curves. It is inherent to the electrical 
power system itself that faults occur in the network, so it is impossible to make the 
same totally immune to those occurrences that will be minimized by the study in 
question. In a second moment of this work it is approached the power flow of the 
electrical system, demonstrating the importance of the same in the planning of the 
conditions of operation, control and supervision of the models of interest of study. This 
mathematical approach with digital tools is the initial proposal of the paper that follows. 
Keywords: Protection relay. Power flow. Simulation. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
22 
REFERÊNCIAS 
 
LIMA, Rogério Lúcio. Proteção de Sistemas Elétricos. 2016. Universidade do 
Estado de Mato Grosso, Faculdade de Ciências Exatas e Tecnologias. Disponível 
em:<http://sinop.unemat.br/site_antigo/prof/foto_p_downloads/fot_13024pse_aula_0
1_pdf_PSE_Aula_01.pdf>. Acesso em: 18 nov. 2018. 
 
COTOSCK, Kelly Regina. Proteção de Sistemas Elétricos: Uma Abordagem 
Técnico Pedagógica. 2007. 109 f. Dissertação (Mestrado) - Curso de Engenharia 
Elétrica, Centro de Pesquisa e Desenvolvimento em Engenharia Elétrica, 
Universidade Federal de Minas Gerais, Belo Horizonte, 2007. Disponível em: 
<https://www.ppgee.ufmg.br/defesas/381M.PDF>. Acesso em: 18 nov. 2018. 
 
SOUZA, Marcos Paulo de Melo Gomes de. Coordenação de proteção de um 
sistema elétrico de potência interligado utilizando relés de sobrecorrente e de 
distância. 2010. 177 f. TCC (Graduação) - Curso de Engenharia Elétrica, 
Departamento de Engenharia Elétrica da Escola Politécnica, Universidade Federal do 
Rio de Janeiro, Rio de Janeiro, 2010. Disponível em: 
<http://monografias.poli.ufrj.br/monografias/monopoli10000578.pdf>. Acesso em: 18 
nov. 2018. 
 
MAMEDE FILHO, João; MAMEDE, Daniel Ribeiro. Proteção de sistemas elétricos 
de potência. 8. ed. Rio de Janeiro: Grupo Editorial Nacional, 2011. 805 p. 
 
MARDEGAM, Claudio. Dispositivos de proteção. 2010. Proteção e seletividade. 
Disponível em: <http://www.osetoreletrico.com.br/wp-
content/uploads/2010/05/Ed50_marco_protecao_seletividade_capIII.pdf>. Acesso 
em: 20 nov. 2018. 
 
GALVEZ, Juan Daniel Carrillo. Modelagem, simulação e analise de fluxo de carga 
da rede elétrica de transporte de guatemala, utilizando software de livre 
acesso. 2009. 193 f. TCC (Graduação) - Curso de Engenharia Elétrica, Faculdade de 
Engenharia, Universidade de São Carlos de Guatemala, Guatemala, 2009. 
 
 
 
 
 
 
 
 
 
 
 
 
23 
ANEXOS 
[1] Configuração número de interações Power World 
 
 
[2] Diagrama Power World 
 
 
24 
[3] Transformador Power World 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
25 
[4] Linhas de transmissão Power World 
 
 
 
 
 
 
 
 
 
 
26 
[5] Barra de referência 01, Power World 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
27 
[6] Barra 02, Power World 
 
 
 
 
 
 
 
 
 
 
 
28 
[7] Barra 03, Power World 
 
 
29 
[8] Barra 04, Power World 
 
30 
 
[9] Resultados no Octave do fluxo de carga 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
31 
[10] Diagrama unifilar do sistema 
 
 
32 
[11] Coordenograma de fase P-1 
 
[12] Coordenograma de fase P-2 
 
 
33 
[13] Coordenograma de fase P-3 
 
[14] Coordenograma de neutro N-1 
 
 
34 
[15] Coordenograma de neutro N-2 
 
[16] Coordenograma de neutro N-3 
 
 
 
35 
[17] Tabela de ajustes calculados para o relé - fase 
 
FUNÇÃO PARÂMETROS VALOR 
CTR fase (Ia, Ib, Ic) RTC 60 
50P1P Fase Instantânea 11,16 
51P1P Fase temporizada 1,67 
51P1C Tipo de curva C2 (muito inversa-IEC) 
51P1TD Tms (multiplicador de 
tempo) 0,5 
51P1RS Reset “N” (not) 
 
[18] Tabela de ajustes calculados para o relé - neutro 
FUNÇÃO PARÂMETROS VALOR 
CTR neutro (In) RTC 50 
50N1P Fase Instantânea de neutro 13,39 
51N1P Fase temporizada de 
neutro 
0,33 
51N1C Tipo de curva C2 (muito inversa-IEC) 
51N1TD Tms (multiplicador de 
tempo) 0,17 
51N1RS Reset “N” (not)