ISSN 2177-9139
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
ANÁLISE DE VIBRAÇÕES INDUZIDAS NA TORRE DE UM AEROGERADOR
Daniela Dalla Chiesa – [email protected]
Universidade Federal de Pelotas, Campus Capão do Leão, 96160-000 - Pelotas, RS, Brasil
Valdecir Bottega - [email protected]
Universidade Federal de Pelotas, Campus Capão do Leão, 96160-000 - Pelotas, RS, Brasil
Resumo. Este artigo tem por objetivo investigar vibrações que ocorrem na torre de
sustentação de um aerogerador, sujeito a uma força vertical harmônica, que pode ser
provocada por vários fatores, dentre eles pode-se citar o desbalanceamento das hélices,
peças defeituosas, desequilíbrio e a ressonância. A análise da vibração permite a
identificação desses problemas atuando como uma importante ferramenta no planejamento
da manutenção preventiva de um parque eólico. Dentre as suas principais vantagens tem-se a
diminuição das perdas de energia e custos operacionais provocados por interrupções
abruptas devido a quebras de componentes mecânicos dos aerogeradores. Neste trabalho,
uma torre de aerogerador foi discretizada pelo Método de Elementos Finitos e algoritmos
foram implementados com o auxílio do software Matlab, para determinar as frequências
naturais, os modos de vibração e a resposta do sistema estudado.
Palavras Chave: Aerogeradores, Elemento de barra, Elementos finitos, Vibrações.
1. INTRODUÇÃO
A notícia da instalação em Santa Vitória do Palmar de um dos maiores parques eólicos da
América Latina (Município de Santa Vitória do Palmar, 2014) e o crescente uso da energia
eólica na matriz energética, aumentou o interesse na análise de situações-problema que
possam interferir no funcionamento de um aerogerador.
As vibrações causadas na estrutura do aerogerador são prejudiciais ao seu funcionamento
e acarretam perda de energia e o aumento de custos, se não tratadas de forma preventiva e
com planejamento.
Uma das principais causas de interrupções do funcionamento de aerogeradores é
provocada por problemas na caixa multiplicadora (conjunto de engrenagens que aumentam a
velocidade de rotação entre o eixo da hélice e o eixo do gerador) (McNiff et al., 1990). Estes
problemas produzem vibrações que são transmitidas para a torre, sendo possível, associar as
vibrações da torre com as vibrações de componentes mecânicos defeituosos (Maalawi, 2007).
Desta forma, a análise de vibrações é uma ferramenta importante para detectar problemas
mecânicos num estado inicial, permitindo a realização de manutenções preventivas e evitando
35
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
paralisações e extensão de danos provocados por quebras abruptas de componentes mecânicos
de aerogeradores (Siqueira, 2012). Vibrações provocadas pelo efeito do vento são
predominantemente transversal ao eixo da torre, enquanto vibrações provocadas por falhas
mecânicas no sistema de transmissão de torque da hélice ao gerador são transmitidas ao eixo
axial da torre (Li, 2010).
Em virtude destes aspectos o presente trabalho visa estudar a vibração axial que ocorre na
torre de sustentação de um aerogerador sujeito a uma força harmônica, simulando vibrações
produzidas por falhas em componentes mecânicas do aerogerador. A importância de tal
análise permite identificar, prematuramente, problemas mecânicos.
Para fins de análise optou-se em utilizar o Método de Elementos Finitos (MEF) (Soriano,
2009), que atualmente vem sendo bastante empregado no estudo de estruturas. Segundo
Ferreira (2007) "o Método dos Elementos Finitos representa uma aproximação dum modelo
matemático que representa o mais fielmente possível o problema físico".
2. MODELO MATEMÁTICO
O objetivo desse artigo é realizar a modelagem da torre da estrutura de um aerogerador,
de forma que as outras massas sejam desconsideradas. A torre será modelada como uma
barra, conforme a Fig. 1, sujeita somente ao deslocamento axial, desconsiderando-se,
portanto, o deslocamento transversal.
x
F
L
E, A, ρ
Figura 1 - Idealização da estrutura de um aerogerador.
A barra possui densidade ρ, módulo de Young1 E, área da secção transversal A,
comprimento total L, eixo axial x e força F.
O primeiro passo é determinar as matriz de massa [M] e de rigidez [K] e o vetor de forças
{F} da equação do movimento [M ]{u}  [ K ]{u}  {F} , onde {u} representa o vetor de
acelerações e {u} representa o vetor de deslocamento. Para encontrar tais matrizes e vetores é
viável utilizar o Método de Elementos Finitos. Tal método consiste em discretizar a barra em
vários elementos e aproximar a solução por funções de interpolação. Para cada elemento
montam-se as matrizes individuais e posteriormente é feita a superposição dos resultados
encontrando as matrizes e vetores globais. De posse dessas informações é possível determinar
as frequências naturais e a realização da análise modal de todo o sistema.
Tomando um elemento de barra conforme mostrado na Fig. 2, onde Le representa o
comprimento do elemento, procura-se uma função aproximada u ( x) que terá como incógnitas
os deslocamentos nodais u1 e u2 do elemento.
1
Também conhecido por módulo de elasticidade.
36
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
Devido ao deslocamento ser axial há um nó em cada extremidade e, portanto, a função de
interpolação é linear.
nó 2
Le
Elemento de
barra
nó 1
Figura 2 - Elemento de barra.
A aproximação u ( x) em cada elemento é dada por:
2
u ( x)   Ni ( x)ui ,
(1)
i 1
onde N i representa as funções de interpolação do elemento dadas por:
x
,
Le
x
N 2 ( x)  e ,
L
N1 ( x)  1 
(2)
(3)
Desta maneira pode-se escrever u ( x) na forma matricial como:
u 
u ( x)   N1 ( x) N 2 ( x)  1   [ N ]{u e } ,
u2 
(4)
onde {u e } representa o vetor dos deslocamentos nodais do elemento e [N] representa a matriz
das funções de interpolação.
Será utilizada a energia cinética e de deformação elástica para a obtenção das matrizes de
massa e de rigidez, respectivamente.
A energia cinética do elemento é definida por (Prodonoff,1990):
2
1 Le
 du 
Te    A   dx,
2 o
 dt 
(5)
du
representa a primeira derivada, em relação ao tempo t, da função u ( x) .
dt
Substituindo Eq. (4) em Eq. (5) obtém-se:
onde
Le
1
Te  {u e }T  A  [ N ]T [ N ]dx{u e } ,
2
0
(6)
37
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
onde {u e } representa a primeira derivada em relação ao tempo, do vetor dos deslocamentos
nodais do elemento, {u e }T é vetor transposto, da primeira derivada em relação ao tempo, dos
deslocamentos nodais do elemento e [ N ]T é a matriz transposta das funções de interpolação.
Substituindo as funções de interpolação Eq. (2) e Eq. (3) em Eq. (6) obtém-se a energia
cinética:
1
Te  {u e }T [ M e ]{u e } ,
2
(7)
onde [ M e ] é a matriz de massa do elemento e pode ser encontrada igualando-se Eq. (7) com
Eq. (6):
Le
[ M ]   A  [ N ]T [ N ]dx .
e
(8)
0
Resolvendo a integral tem-se a matriz de massa do elemento que é dada por:
[M e ] 
 ALe  2 1 
6
1 2  .


(9)
A energia de deformação do elemento é (Prodonoff,1990):
2
1 Le
 du 
Ve   AE   dx,
o
2
 dx 
(10)
du
representa a primeira derivada, em relação a x, da função u ( x) .
dx
Então, para um dado elemento de barra a energia de deformação é dada substituindo Eq.
(4) em Eq. (10) de modo que:
onde
Le
1
Ve  {u e }T AE  [ N ' ]T [ N ' ]dx{u e },
2
0
(11)
onde [ N ']T representa a matriz transposta, da primeira derivada em relação a x, das funções de
interpolação e [ N ' ] representa a matriz da primeira derivada, em relação a x, das funções de
interpolação.
Substituindo as funções de interpolação Eq. (2) e Eq. (3) em Eq. (11) encontra-se a da
energia de deformação:
1
Ve  {u e }T [ K e ]{u e } ,
2
(12)
onde [ K e ] é a matriz de rigidez do elemento e pode ser obtida igualando-se Eq. (12) com Eq.
(11):
38
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
Le
[ K ]  EA  [ N ' ]T [ N ' ]dx .
e
(13)
0
Resolvendo a integral encontra-se a matriz de rigidez do elemento:
[K e ] 
EA  1 1
.
Le  1 1 
(14)
Tendo encontrado as matrizes de massa e rigidez do elemento, resta encontrar o vetor de
forças.
O trabalho realizado pelas forças distributivas do elemento é expresso como (Ferreira,
2007):
Le
 W e   p udx ,
o
(15)
onde p é o carregamento aplicado por unidade de comprimento e  u o deslocamento virtual.
Substituindo Eq. (4) em Eq. (15), obtém-se:
Le
 W e   {u e }T  p[ N ]T dx .
o
(16)
ou W e   {u e }T {F e } (Ferreira, 2007) então pode-se escrever:
Le
{F e }   p[ N ]T dx .
o
(17)
onde {F e } representa a o vetor de forças externas aplicada em cada elemento.
Substituindo a matriz das funções de interpolação em Eq. (17) e supondo uma força
constante p, chega-se ao vetor:
{F e } 
pLe
2
1
1 .

(18)
Sendo assim, a equação do movimento do elemento é dada por:
 ALe  2 1   u1  EA  1
6
1 2 u   Le  1

 2

1  u1  pLe 1
.
 
1  u2 
2 1
(19)
3. MATRIZES E VETORES GLOBAIS DA BARRA
Será montada a equação [M ]{u}  [ K ]{u}  {F} para todo o sistema e isto é feito
realizando a superposição de todos os elementos da barra, com o intuito de exemplificar o
procedimento será feita uma discretização em quatro elementos e posteriormente os cálculos
serão estendidos a toda a estrutura, portanto, toma-se uma barra conforme mostrado na Fig. 3.
39
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
x
F
5
Elemento 4
4
Elemento 3
3
Elemento 2
2
1
Elemento 1
Figura 3 - Barra dividida em quatro elementos.
Observando a Fig. 3 pode-se ver que a barra está dividida em quatro elementos. O
elemento 1 é composto pelos nós 1 e 2, o elemento 2 é composto pelos nós 2 e 3, o elemento
3 é composto pelos nós 3 e 4 e o elemento 4 é composto pelos nós 4 e 5. Verifica-se que cada
elemento possui um nó que o conecta com o outro.
Montam-se as matrizes de massa e rigidez e o vetor de forças para cada um dos
elementos, neste momento, pode-se considerar que cada um dos elementos tenha módulo de
elasticidade, densidade, área da seção transversal e comprimento diferentes.
A matriz de massa do sistema é obtida a partir da superposição das matrizes de massa de
cada elemento, como segue abaixo:
 1 A1 Le1
 3

 1 A1 Le1
 6

[M ]   0



 0

 0

1 A1 Le1
0
0
2 A2 Le 2
 2 A2 Le 2
3
6
0
6
 AL
e
1 1 1
3

 2 A2 Le 2
 2 A2 Le 2
6
3

3 A3 Le3
0
3 A3 Le3
3
6
3 A3 L 3
3 A3 L 3
6
3
e
e

 4 A4 Le 4
3
 4 A4 L 4
e
0
0
6




0 

 (20)
0 

 4 A4 Le 4 
6 

 4 A4 Le 4 
3 
0
A matriz de rigidez é obtida a partir da superposição das matrizes de rigidez de cada
elemento, como segue abaixo:
40
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
 E1 A1
 Le
1

 E1 A1
  Le
1


[K ]   0


 0


 0

E1 A1
Le1

0
E1 A1 E2 A2
 e
Le1
L2

E2 A2
Le 2

E2 A2
Le 2
E2 A2 E3 A3
 e
Le 2
L3

0
0
0
E3 A3
Le3
0
0

E3 A3
Le3
E3 A3 E4 A4
 e
Le3
L4

E4 A4
Le 4




0 


0 

E4 A4 
 e 
L4 
E4 A4 

Le 4 
0
(21)
Os vetores de aceleração e de deslocamento são dados, respectivamente por:
{u}  {u1 u2
u3 u4
u5 }T
{u}  {u1 u2
u3 u4
u5 }T
(22)
E o vetor de forças é dado por:
 p Le
{F }   1 1
 2
p1Le1 p2 Le 2

2
2
p2 Le 2 p3 Le3

2
2
p3 Le3 p4 Le 4

2
2
T
p4 Le 4 

2 
(23)
4. MÉTODO DE ELEMENTOS FINITOS
Os cálculos feitos acima se referem a um problema genérico de uma barra, contudo
precisam-se adaptar as matrizes as condições de contorno do problema em questão, portanto
há a necessidade da reformulação das matrizes globais de massa, de rigidez e do vetor de
forças do sistema.
A barra da Fig. 1 tem área de secção transversal constante A, densidade ρ, módulo de
elasticidade E e comprimento total L e será discretizada em 9 elementos de igual
L
comprimento, portanto o comprimento Le de cada elemento é .
9
Deve-se, também, observar que a barra está engastada em x=0 e livre em x=L, portanto
nas matrizes [M] e [K] podemos zerar alguns componentes das primeiras linhas e colunas de
forma que o deslocamento nodal u1  0 . Assim a matriz de massa e rigidez são reescritas,
respectivamente, por:
41
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
1
0

0

0
 AL  0
[M ] 

54  0
0

0
0

 0
1
0

0

0
9 EA 0
[K ] 

L 0
0

0
0

0
0
0
0
0
0
0
0
0
4
1
0
0
0
0
0
0
1
4
1
0
0
0
0
0
0
1
4
1
0
0
0
0
0
0
1
4
1
0
0
0
0
0
0
0
4
1
0
0
0
0
0
0
1
4
1
0
0
0
0
0
0
1
4
1
0
0
0
0
0
0
1
4
0
0
0
0
0
0
0
1
0
0 
0

0
0

0
0

0
1

2 
(24)
0
0 
1 2 1 0
0
0
0
0
0

0 1 2 1 0
0
0
0
0
0
0 1 2 1 0
0
0
0

0
0
0 1 2 1 0
0
0
0
0
0
0 1 2 1 0
0

0
0
0
0
0 1 2 1 0 
0
0
0
0
0
0 1 2 1

0
0
0
0
0
0
0 1 1 
0
0
0
0
0
0
0
0
2
1
0
0
0
0
0
0
(25)
Há uma força harmônica que atua no último elemento da barra tomada por
F  F0 cos(wt ) , onde F0 é a amplitude e w é a frequência da excitação, dessa maneira a
matriz de forças é reescrita com a excitação aplicada na última linha, conforme segue abaixo:
T


L
{F }  0 0 0 0 0 0 0 0 0
F0 cos( wt ) 
18


(26)
5. ANÁLISE MODAL
No Método de análise modal é usado o teorema de expansão e os deslocamentos das
massas são expressos como uma combinação linear dos modos normais do sistema. Essa
transformação linear desacopla as equações de movimento, de modo que obtemos um
conjunto de n equações diferenciais de segunda ordem não acopladas. (Rao, 2008).
O primeiro passo é encontrar a solução do problema de autovalor e determinar as
frequências naturais e os modos de vibração. Em vibrações livres desconsideram-se forças
externas, ou seja, tem-se a seguinte equação:
[M ]{u}  [ K ]{u}  0,
(27)
Substituindo u(t )  Xsen(wt   ) na Eq. (27) obtém-se o problema de autovalor
generalizado:
[ K ]  w2 [ M ] X  0,
(28)
42
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
onde w2 são os autovalores, X os autovetores e  o ângulo de fase.
A Eq. (28) é conhecida como problema de autovalor ou do valor característico. A única
maneira de ser ter uma solução não nula para a Eq. (28) é que [ K ]  w2 [ M ] seja uma
matriz singular, dessa forma:
det [ K ]  w2 [M ]  0.
(29)
Tal condição permite encontrar a equação característica para a determinação das
frequências naturais e modos de vibração do sistema. As frequências naturais são obtidas
pelos autovalores da equação característica e as formas modais pelos autovetores.
Para encontrar tais valores, considerou-se uma torre de 18m de altura, com área de secção
transversal A= 12m2, densidade ρ = 7.830 kg/m3 e módulo de Young E= 210 Gpa. Com o
auxílio do software Matlab chegam-se aos resultados das frequências naturais e dos modos de
vibração, conforme Tabela 1 abaixo e Fig. 4, respectivamente.
Tabela 1- Frequências naturais e modos de vibração
Frequência
Valor em
Hz
w1
72,02
w2
218,26
w3
371,13
w4
535,06
w5
713,80
w6
908, 20
w7
1.110,61
w8
w9
1.294,95 1.411, 49
Figura 4 – Primeiros 4 modos de vibração
A equação modal do sistema é dada por (Rao, 2008):
q(t )  [2 ]q(t )  Q(t ) .
(30)
onde q(t ) é o vetor de aceleração em coordenadas generalizadas e q(t ) é o vetor de
deslocamento em coordenadas generalizadas e:
[2 ]  [X]T [K][X] e Q(t )  [ X ]T {F}
(31)
43
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
onde [ X ]T é a matriz transposta modal normalizada pela massa e {F } é o vetor de forças
tomado por:
{F}  [0 0 0 0 0 0 0 0 0 10cos(20t )]
(32)
Tendo encontrado todas as matrizes e vetores da Eq. (30), basta calcular a solução do
sistema de equações. Pode-se perceber que esta equação é uma equação diferencial de
segunda ordem não-acoplada, que descreve o movimento de um sistema não-amortecido com
um grau de liberdade. Sua solução é determinada com a utilização do software Matlab,
contudo a equação deve ser expressa como um sistema de equações diferenciais de primeira
ordem acopladas. Para tanto, realiza-se a transformação da referida equação em espaço de
estado.
Considerou-se o problema com condições iniciais de deslocamento nulas e um intervalo
de tempo de 0,5s. Com a implementação do algoritmo no software Matlab, empregou-se o
Método de Runge Kutta de quarta ordem, encontrando a resposta da equação modal,
conforme mostrado na Fig. 5, representando as vibrações induzidas na torre por uma força
harmônica, simulando a força provocada por um componente mecânico defeituoso no
aerogerador.
Figura 5 - Resposta da vibração da torre sujeita a uma força harmônica.
6. CONSIDERAÇÕES FINAIS
Neste trabalho, aplicou-se o Método dos Elementos Finitos na análise de estruturas,
possibilitando a determinação das matrizes e vetores para realização da análise modal da torre
de um aerogerador. Consequentemente, o MEF em conjunto com a análise modal, mostrou ser
bastante viável para o estudo realizado.
A utilização do software Matlab possibilitou a implementação de algoritmos que
realizassem os cálculos requeridos e simulações, mostrando as vibrações induzidas na torre de
sustentação de um aerogerador por uma força harmônica no sentido vertical. Este estudo
complementa resultados experimentais encontrados na literatura, onde sensores são colocados
na torre de aerogeradores para detectar vibrações provocadas especificamente por problemas
44
XX EREMAT - Encontro Regional de Estudantes de Matemática da Região Sul
Fundação Universidade Federal do Pampa (UNIPAMPA), Bagé/RS, Brasil. 13-16 nov. 2014.
mecânicos (Siqueira, 2012) e mostra que a modelagem matemática destes sistemas pode ser
usada para projetos de manutenção preventiva, reduzindo custos e aumentando a
confiabilidade de sistema de energia eólica.
REFERÊNCIAS
FERREIRA, A. J. M. Problemas de Elementos Finitos em Matlab. Fundação Calouste
Gulbenkian, 2007.
LI, J. C. J.; CHEN. X. Aerodynamic analysis of wind turbines. Journal of Mechanical
Science and Technology, vl. 25, n. 1 p. 89-95, 2010.
MAALAWI, K. Y. A model for yawing dynamic optimization of a wind turbine structure.
International Journal of Mechanical Sciences, v. 49, n. 10, p. 1130–1138, 2007.
MCNIFF, B.; MUSIAL, W. D.; ERRICHELLO, R. Variations in Gear Fatigue Life for
Different Wind Turbine Braking Stategies, Golden, Colorado USA. Solar Energy Research
Institute, 1990.
MUNICÍPIO DE SANTA VITÓRIA DO PALMAR. Energia Eólica Desenvolvimento que
vem com os ventos. Acessado em 10 jul. 2014. Online. Disponível em:
http://www.santavitoria.rs.gov.br/portal1/municipio/noticia.asp?iIdMun=100143344&iIdNoti
cia=216801.
PRODONOFF,V. Vibrações Mecânicas – Simulação e Análise. Rio de Janeiro: Maity
Comunicação e Editora Ltda, 1990.
RAO, S. Vibrações Mecânicas. São Paulo: Pearson Prentice Hall, 2008.
SIQUEIRA, C. D. A análise de vibrações como ferramenta para a melhoria da
manutenção em Aerogeradores. 2012. 181f. Dissertação para obtenção do Grau de Doutor
em Engenharia Mecânica - Faculdade de Ciências e Tecnologia, Universidade Nova de
Lisboa.
SORIANO, H. L. Elementos Finitos: Formulação e Aplicação na Estática e Dinâmica das
Estruturas. Rio de Janeiro: Ciência Moderna, 2009.
45