ii

WEGeuro Indústria Eléctrica S. A.

Departamento de Projecto Mecânico

Eduardo Afonso Ribeiro

ANÁLISE DINÂMICA DE MOTORES

ELÉCTRICOS: MODELO E CONTROLO

Projecto apresentado como conclusão do período de estágio curricular na empresa WEGeuro, desenvolvido na disciplina de Projecto do curso de Mestrado Integrado em Engenharia Mecânica, do Departamento de Engenharia Mecânica e Gestão Industrial da Faculdade de Engenharia da Universidade do Porto.

Orientadores:

André F. C. Rodrigues

Pedro M. L. Ribeiro

PORTO

04/07/2009

iii

À minha família.

iv

RESUMO

A ciência tem avançado com extrema velocidade e a tecnologia ultrapassa os limites antigos a cada minuto. As necessidades do mercado exigem cada vez mais máquinas rotativas operando em elevadas rotações, com o objectivo de aumentar o ganho de potência através de motores mais leves e mais eficientes. A análise dinâmica destas estruturas é muito importante, assim como a precisão dos resultados gerados pelo modelo numérico é da mesma forma necessária. Um estudo teórico e uma modelagem dinâmica de rotores utilizando o método numérico dos elementos finitos, através de elementos de viga de Timoshenko com elementos de 3 nós e funções de interpolação quadráticas é realizado. A metodologia utilizada na obtenção das matrizes de massa, rigidez, amortecimento e giroscópica é apresentada. Técnicas usuais de controlo, tais como mancais de filme de óleo e alteração na geometria do veio são apresentadas. Outras alternativas de controlo como absorsores de vibração e utilização de técnicas de optimização para determinação da equilibragem óptima também estão presentes neste trabalho. Finalmente casos de estudo são realizados e devidamente discutidos.

v

ABSTRACT

The science has advanced quickly and each minute technology surpasses the old limits. Industrial market needs require more and more rotating machines running in high speeds, with the objective of increasing the power gain through more efficient and lighter motors. The dynamic analysis of these machines is very important, as well as the accuracy of the results generated by the numerical models employed in the anlyses. A theoretical study and a dynamic model for rotors is developed in this dissertation. The numerical approach uses a finite element model, based on a Timoshenko beam quadratic finite element, with 3 nodes per element, four degrees of freedom per node and quadratic shape functions of the C0 class. The methodology used to derive the mass, stiffness, damping and gyroscopic matrices is presented. Usual techniques of vibration control, such as oil film bearings and mass test balancing are also presented in this work. Additionaly, other alternatives of vibration control, like vibration absorbers and optimization techniques for determination of the balance are presented. Finally, case studies are carried out and properly discussed.

vi

AGRADECIMENTOS

Agradecimentos a toda a minha família pelo apoio incondicional em todos os momentos.

Agradecimentos a todos os professores e engenheiros com os quais adquiri grande parte de meus

conhecimentos, em especial aos professores Carlos Alberto Bavastri, Dr Eng. Marco Antônio

Luersen, Dr. Eng. e Jucélio Tomás Pereira Dr. Eng.

Agradecimentos ao meu colega Aleksander Kokot, autor de [1], pela ajuda na realização deste

trabalho.

Agradecimentos também à WEG que nos ajudou nas questões relativas ao estágio na WEGeuro,

em especial à Hilton Penha Silva, e Silvana Tecila. Também aos colaboradores Francisco José

Doubrawa Filho e Hideraldo Luís Vasconcelos dos Santos.

À WEGeuro por propiciar a realização do estágio em suas instalações, à toda equipe do Projecto

Mecânico, em especial à aos engenheiros André Rodrigues e Luís Araújo.

Agradecimentos ao laboratório pela colaboração na parte experimental, em especial a Domingos

Manuel Silva Oliveira Duarte e Rui Manuel Borges Moura Guedes.

Agradecimentos a todas as pessoas que tornaram este projecto possível, em especial prof. Pedro

Leal Ribeiro da FEUP.

vii

ÍNDICE DE CONTEÚDOS

1. INTRODUÇÃO ........................................................................................................ 9 1.1 Descrição dos capítulos ................................................................................................................ 12 1.2 Motores Eléctricos ........................................................................................................................ 13

2. MODELO DE ELEMENTOS FINITOS ................................................................ 14 2.1 Obtenção da Matriz de Rigidez .................................................................................................... 14

2.1.1 Deformação no veio ............................................................................................................. 16 2.1.2 Energia Potencial de Deformação do Veio .......................................................................... 18 2.1.3 Efeito do Carregamento Axial ............................................................................................. 19 2.1.4 Formulação do Elemento Finito de Viga com Interpolação Quadrática (3 nós) .................. 19

2.2 Obtenção da Matriz de Massa e Giroscópica ............................................................................... 22

2.2.1 Composição de Rotações ..................................................................................................... 22 2.2.2 Energia Cinética do Disco ................................................................................................... 23 2.2.3 Formulação do Elemento Finito para o Disco ..................................................................... 24 2.2.4 Energia Cinética do Veio ..................................................................................................... 24 2.2.5 Formulação do Elemento Finito com Interpolação Quadrática (3 nós) ............................... 24

2.3 Modelo para Mancais e Fundação ................................................................................................ 26 2.4 Modelo de Interferência entre Pacote de Chapas e Veio .............................................................. 29

3. METODOLOGIA DE CÁLCULO ......................................................................... 30

4. ALTERNATIVAS DE CONTROLO DE VIBRAÇÃO ......................................... 31 4.1 Absorvedores de Vibração ....................................................................................................... 31 4.2 Equilíbrio Óptimo......................................................................................................................... 37

4.3 Apoios ......................................................................................................................................... 39

5. RESULTADOS E COMPARAÇÕES .................................................................... 44 5.1 Interface Gráfica ........................................................................................................................ 44

5.2 Comparação dos Resultados Numéricos com os Dados da RENK .................................... 45

5.3 Modelo de Elementos Finitos Aplicado ao Cálculo dos Parâmetros Modais de Rotores: Validação Experimental ............................................................................................................................. 52

5.4 Equilibragem Teórica ................................................................................................................ 55

5.5 Projecto de Absorsor de Vibrações Viscoelástico ................................................................ 62

5.6 Controle de Vibração Através de Selecção de Chumaceiras ............................................... 65

6. CONCLUSÕES ....................................................................................................... 76

7. TRABALHOS FUTUROS ...................................................................................... 77

8. REFERÊNCIAS ...................................................................................................... 78

APÊNCICE A .................................................................................................................. 81

APÊNDICE B .................................................................................................................. 85

viii

APÊNCICE C ................................................................................................................... 86

APÊNDICE D .................................................................................................................. 88

9

1. INTRODUÇÃO

Em equipamentos mecânicos em geral observa-se a presença de vibração. Este fenómeno

pode ser agravado se a estrutura atingir a ressonância, já que as amplitudes do movimento

oscilatório aumentam consideravelmente. A ressonância ocorre quando a frequência de

excitação coincide com uma das frequências naturais, para sistemas com baixo amortecimento,

o que deve ser evitado, do mesmo modo que em sistemas amortecidos nos quais a ressonância

ocorre em frequências ligeiramente abaixo das frequências naturais. Em máquinas rotativas, faz-

se extremamente necessário controlar os níveis de vibração, pois, actualmente cada vez mais se

têm máquinas maiores operando em rotações mais elevadas, e quanto maior a velocidade de

rotação maiores as influências de forças actuantes no sistema e maiores os níveis de vibração,

logo alternativas de controlo devem ser incorporadas ao sistema.

Baseado no modelo geral apresentado por [1] é possível predizer o comportamento

dinâmico de um rotor. A partir deste modelo numérico é apresentada uma metodologia de

cálculo para determinação dos parâmetros modais de um sistema rotativo. O cálculo dos

parâmetros modais de um rotor pode tornar-se muito complexo a medida que se considera não

linearidades do sistema tais como, mancais de filme de óleo, materiais viscoelásticos,

enrijecimento devido ao esforço axial, entre outros. As alternativas de controlo de vibração de

sistemas girantes apresentadas neste trabalho dependem de uma modelagem adequada do

sistema, logo, a precisão na predição do comportamento dinâmico do sistema afecta

directamente a eficiência dos mecanismos de controlo propostos.

Um dos métodos de controlo de vibrações mais difundido e utilizado para rotores é a

equilibragem. Na década de 60 já se utilizava para rotores rígidos em baixa rotação a técnica do

balanceamento modal. Este método consiste na correcção da componente do desequilíbrio

através da adição ou subtracção de massas de acordo com os sucessivos modos de vibrar, sendo

equilibrado um modo de cada vez sem que o equilíbrio de um determinado modo influencie no

modo posterior. Para cada alteração necessita-se de um teste de eficiência da massa

acrescentada ou removida, segundo [2].

Na década de 70 a correcção do desequilíbrio através da técnica dos coeficientes de

influência tornou-se uma alternativa de controlo das vibrações. Diferentemente da equilibragem

modal, a técnica baseada nos coeficientes de influência da mecânica consiste na determinação

das massas de correcção através de várias medições experimentais realizadas com diferentes

massas posicionadas nos diversos planos de correcção pré definidos. Nota-se, que esta técnica

exige mais testes relativamente á técnica do balanceamento modal (referência [2]).

10

Em [3] pode-se observar uma comparação entre a equilibragem de rotores flexíveis

utilizando os métodos da equilibragem modal e da equilibragem através de coeficientes de

influência. A proposta inicial dos autores é obter uma única expressão unificando os dois

métodos de equilibragem, além da obtenção de vectores independentes da rotação de

equilibragem, sendo possível então equilibrar rotores flexíveis em baixa rotação em múltiplas

velocidades, porém, existe ainda a necessidade de aplicar a equilibragem em altas rotações em

alguns rotores específicos, como turbinas de aviões.

A equilibragem em múltiplas velocidades, ou modo a modo, também é estudada em [4].

Com planos de equilibragem bem definidos e um modelo através de elementos finitos, baseado

no conceito de que os melhores resultados para a equilibragem são obtidos, adoptando planos no

local onde se encontra o desequilíbrio, vários testes com bons resultados são apresentados

verificando a técnica.

Métodos mais recentes, como a equilibragem automática derivada de outras técnicas de

balanceamento dinâmico, são propostos para reduzir vibrações em rotores. O balanceamento

automático consiste em posicionar esferas livres que podem girar em torno de discos para

compensar o desequilíbrio do rotor. Em [5] os autores aplicam o método no modelo de rotor

conhecido como Stodola - Green [6, 7] que possibilita a consideração de um rotor flexível

composto por um disco (corpo rígido).

Para a equilibragem automática pode-se derivar a equação do movimento e obter uma nova

equação autónoma para cada sistema. A referência [5] utiliza-se de um modelo que proporciona

bons resultados em rotores flexíveis equilibrados entre a primeira e a segunda rotação crítica,

quando se têm parâmetros modais estáveis. Um modelo geral para rotor rígido e uma

formulação lagrangeana é utilizada em [8].

Outra alternativa inovadora que consiste na equilibragem sem massas de testes é

apresentada por [9]. O método se baseia na determinação das massas e fases do desequilíbrio,

através de uma técnica de optimização. Para isso é necessário um método numérico capaz de

predizer o comportamento dinâmico do sistema com precisão. O método numérico utilizado

pelo autor é o método da matriz de transferência, aliada a uma técnica de optimização como os

algoritmos genéticos.

Neste trabalho a técnica de equilibragem apresentada baseia-se no trabalho proposto por [9],

diferenciando-se pelo método numérico utilizado, neste caso o método dos elementos finitos e

pela técnica de optimização de algoritmos genéticos elitistas e de nicho compartilhado.

11

Outra alternativa no controlo de vibrações é através de absorsores, que são dispositivos

utilizados para dissipar energia do sistema vibratório. A referência [10] apresenta mecanismos

optimizados e muito eficientes no controle de vibração flexional em rotores, chamados

absorsores de vibração viscoelásticos. Com elevado amortecimento estes sistemas massa – mola

- amortecedor oferecem a estruturas girantes a capacidade de, mesmo que originalmente com

baixo amortecimento, trabalharem acima de suas rotações críticas. Através de um modelo de

elementos finitos preciso do sistema a ser controlado, o absorsor pode ser projectado para

reduzir as amplitudes de vibração do sistema na ressonância.

Seguindo a mesma metodologia apresentada em [10], [11] apresenta absorsores de vibração

torsional. Os parâmetros do rotor são obtidos também através do método numérico dos

elementos finitos e materiais viscoelásticos são os elementos dissipadores de energia.

Neste trabalho o projecto de um absorsor dinâmico viscoelástico para um rotor flexível é

apresentado, e o comportamento dinâmico destes tipos de material é brevemente discutido. O

projecto é realizado através da atribuição de uma massa considerada eficiente para a massa total

do rotor controlado e seus parâmetros obtidos de acordo com limitações geométricas impostas

pelos demais componentes do motor eléctrico. A proposta consiste em tentar controlar as duas

primeiras velocidades rotações críticas.

Mancais de filme de óleo são uma outra forma de controlo das amplitudes de vibração de

um motor eléctrico, porém é difícil calcular os coeficientes de rigidez e amortecimento que são

muito influenciados pela temperatura e velocidades de trabalho, além da influência da carga a

que o rotor está submetido.

Na referência [12] apresentam-se alguns fenómenos típicos que ocorrem em rotores

apoiados sobre mancais de filme de óleo, e ainda uma análise de estabilidade. O autor

demonstra a utilização dos números de Sommerfeld e Ocvirk com o objectivo de determinar a

posição do centro do eixo em relação ao centro do mancal, pois esta é predominantemente

determinada pelo tipo de escoamento dentro do mancal. Na sequência [13] faz uma

diferenciação simplificada entre mancais curtos e longos.

Em [14] Farias e Ribas apresentam uma comparação entre os resultados obtidos através da

teoria descrita por [12] com os valores calculados pela rotina comercial de cálculo da fabricante

de mancais de filme de óleo RENK.

Os dados referentes a dimensões dos mancais, excentricidade, tipo de lubrificante, folga,

reacções nos mancais e temperatura, influenciam nos coeficientes de rigidez e amortecimento,

assim como a frequência de rotação. Com foco na redução das amplitudes de vibração do rotor,

12

estes parâmetros são modificados de forma conveniente a produzir coeficientes de

amortecimento e rigidez que resultem na redução das amplitudes de vibração.

Neste trabalho são apresentadas metodologias de controlo de vibração, ainda baseada na

sequência de equações para modelagem numérica de um rotor, [1], apresenta-se a interface

gráfica implementada. O controlo de vibração é dependente de uma modelagem precisa, logo

este trabalho pode ser entendido como uma sequencia à referência [1]. Simulações numéricas

são realizadas assim como seus respectivos resultados são discutidos.

Em resumo este trabalho tem por objectivos propor um modelo de elementos finitos para

predição do comportamento dinâmico de rotores, baseado na metodologia de desenvolvimento

apresentada em [1]. Porém o modelo de elementos finitos, neste trabalho, considera elementos

com três nós consequentemente os elementos possuem 12 graus de liberdade, e funções de

forma quadráticas. A proposta de desenvolvimento e implementação dos elementos quadráticos,

pode ser considerada basicamente como uma verificação da precisão dos elementos lineares,

considerando-se que os elementos quadráticos representam melhor o fenômeno físico.

Três propostas de controlo de vibração fazem parte dos objectivos, dentre elas,

balanceamento óptimo, e controlo através de mancais de filme de óleo, recebem maior ênfase na

parte dos resultados. Da proposta de absorsores de vibração é apresentado apenas o seu modelo

e projecto, pois, devido à escassez de tempo não foi possível avaliar a sua influência no sistema,

tanto em simulações numéricas quanto experimentais.

1.1 Descrição dos Capítulos

Na secção 3 os objectivos do trabalho são descritos, assim como a necessidade do estudo de

métodos de controlo é explicada. Na secção 4 o modelo de elementos finitos para um rotor é

apresentado. Baseado na teoria de viga de Timoshenko um elemento com funções de

interpolação quadráticas e 3 nós por elemento finito é mostrado.

No capítulo 5 a metodologia de cálculo utilizada para a obtenção dos parâmetros modais do

rotor é apresentada. Na sexta secção um estudo mais detalhado das alternativas de controlo

apresentadas é realizado.

No capítulo 7 os resultados são apresentados a fim de validar o modelo teórico descrito no

trabalho. As conclusões obtidas após a aplicação dos métodos de controlo são mostradas e o

projecto de um absorsor é apresentado.

13

1.2 Motores Eléctricos

A gama de motores produzidos pela WEGeuro é composta em sua maioria por rotores

rígidos, porém para a produção de motores com maiores rendimentos, rotores mais leves para

maiores potências produzidas faz-se necessário o projecto de rotores flexíveis. Para rotores em

que se tornam inviáveis modificações geométricas, ou a utilização de mancais especiais, pode-se

fazer uso de mecanismos de controlo de vibração.

Além do projecto mais preciso que deve ser realizado quando se tem rotores flexíveis, são

produzidos na WEGeuro motores anti-deflagrantes, capazes de trabalhar em atmosferas

explosivas Figura 1. Estes tipos de motores, devido ao ambiente no qual estão inseridos

possuem diversas particularidades, das quais, pode-se citar o espaço entre ferro reduzido. Este

espaço reduzido, distância entre o rotor e o estator, justifica o facto de as amplitudes de vibração

do veio deverem ser controladas a fim de evitar um toque entre estator e rotor, o que poderia

provocar uma explosão.

Logo, adopta-se como foco para o controlo das vibrações nos motores eléctricos o veio.

Componentes como carcaça, tampas, massa rotórica, entre outros citados em [1] estão presentes

na maioria dos motores eléctricos, e apesar de influenciarem no comportamento vibratório do

sistema, controlo de vibração aplicado a estes dispositivos são geralmente ineficientes devido ao

facto de que a principal e predominante fonte de vibração em motores eléctricos é o veio no

qual devem ser feitas as intervenções.

14

Figura 1 - Motor anti-deflagrante W22XC 355MLJ IMB3T [15].

2. MODELO DE ELEMENTOS FINITOS

Na modelagem de um rotor se faz uso de diversos métodos matemáticos e numéricos, tais

como, dentre outros, dinâmica de meios contínuos, modos assumidos, elementos finitos.

Um dos métodos numéricos mais utilizados é o método dos elementos finitos. Ruhl e

Booker [16] já utilizavam este método para modelar um rotor com base na teoria de Euler.

Ainda considerando a teoria de viga de Euler Bernoulli, em [17] é inserido um factor de

correcção directamente na matriz de rigidez para considerar o efeito do corte transverso. No

modelo os efeitos giroscópico e de inércia de rotação são considerados. Seguindo a mesma

metodologia, em [18] é utilizado para os elementos finitos, uma formulação com teoria de viga

de Timoshenko, e funções de forma da classe C0 e ainda com a consideração dos efeitos do

amortecimento viscoso e histerético.

Neste trabalho, sequência de [1], é utilizado o método dos elementos finitos. A partir da

teoria de Timoshenko com funções de forma da classe C0, considerando os efeitos giroscópico e

inércias de rotação, obtêm-se as equações do movimento. Os mancais são modelados de forma

geral, assim pode-se considerar um número qualquer de mancais localizados em qualquer

posição aleatória ao longo do comprimento do veio, e ainda propriedades variáveis com a

rotação. De acordo com resultados obtidos em [1], o modelo mais adequado para considerar o

enrigecimento no veio devido a inserção do pacote de chapas com interferência é apresentado.

A variação da rigidez do sistema devido ao carregamento axial e ao campo magnético são

15

considerados e seus modelos brevemente apresentados. As forças de excitação consideradas são:

o desequilíbrio e excitações externas harmónicas.

2.1 Obtenção da Matriz de Rigidez

A matriz de rigidez para um elemento finito de viga, como apresentado por [1], pode ser

obtida à partir do estudo da energia potencial de deformação. Para o entendimento da

deformação sofrida pela viga pode-se considerar os estados não deformado R e deformado

*R de um segmento de recta PQ de um corpo apresentados na Figura 2.

y

X

Z

Q(x+dx, y +dy, z+dz)

P(x, y, z)

Q* (x*+dx, y* +dy, z*+dz)

P* (x*, y*, z*)

ds

ds*

z

Y

x O

R*

R

Figura 2: Segmento linear PQ no corpo não deformado (-) e deformado (--) [1].

Assim pode-se obter:

222

2

1

x

w

x

v

x

u

dx

duxx , 2.1

222

2

1

y

w

y

v

y

u

dy

dvyy , 2.2

222

2

1

z

w

z

v

z

u

dx

dwzz , 2.3

16

y

w

x

w

y

v

x

v

y

u

x

u

y

u

x

vyxxy 2

1 , 2.4

z

w

x

w

z

v

x

v

z

u

x

u

z

u

x

wzxxz 2

1 2.5

e

z

w

y

w

z

v

y

v

z

u

y

u

z

v

y

wzyyz 2

1 . 2.6

Devido à hipótese de consideração de pequenos deslocamentos, os termos quadráticos

podem ser desconsiderados. xx , yy e zz são as deformações longitudinais do elemento

diferencial dS paralelas à x , y e z respectivamente.

Estas deformações estão deduzidas para um estado tridimensional de deformações, na

secção seguinte as mesmas serão escritas para o caso de um veio. Baseado nas equações das

deformações para um veio, pode-se obter a energia potencial de deformação do mesmo.

2.1.1 Deformação no veio

A Figura 3 apresenta os deslocamentos da secção transversal de um veio aos sistemas

coordenados global XYZ e móvel xyz , sendo tΩ o ângulo de rotação, e B um ponto

qualquer da secção do veio, com de coordenadas 11, zx . u e w são os deslocamentos

relativos às coordenadas globais X e Z enquanto *u e *w são os deslocamentos relativos

às coordenadas locais x e z .

C

x1

z1

w

w*

u*

u

Z

X

x

z

O

Ωt B

Figura 3 – Coordenadas do centro geométrico C e um ponto arbitrário B no veio, numa vista transversal

ao seu comprimento [1].

17

Para relacionar os deslocamentos entre os dois sistemas de referência apresentados na

Figura 3, utiliza-se uma matriz de transformação de coordenadas e de acordo com [1] pode-se

chegar a:

twtuu sin*cos* , 2.7

e

twtuw cos*sin* . 2.8

Dos deslocamentos mostrados nas equações 2.7 e 2.8 obtem-se as deformações do veio.

Duas hipóteses distintas para determinar essas deformações são mostradas em [1], e dependem

directamente da teoria de viga utilizada.

A Figura 4 mostra uma hipótese para deformação de uma viga, conhecida como teoria das

vigas finas, ou teoria de Euler Bernoulli na qual a rotação é igual à derivada do deslocamento

transversal. A Figura 5 apresenta a hipótese de Timoshenko que considera a influência do corte

transversal, assim a rotação da secção transversal do veio já não é igual a derivada do

deslocamento transversal, devido a existência das deformações de corte .

Figura 4 – Deslocamento em uma viga de Euler - Bernoulli, para as coordenadas y, z e y, x [1].

18

Figura 5 – Deslocamento em uma viga de Timoshenko, para as coordenadas y, z e y, x [1].

De acordo com [1] a expressão do deslocamento axial pode ser escrita como:

)()(* yzyxvv , 2.9

na qual as rotações podem ser escritas como:

)(*

yy

wuz

2.10

e

)(*

yy

uyz

. 2.11

yx e yz são as deformações de corte nos planos yx e yz respectivamente. Assim

Substituindo os deslocamentos, dados pelas equações 2.7, 2.8 e 2.9 na equação da deformação

longitudinal do veio (equação 2.2), chega-se a:

22

*

2

1*

2

1

y

w

y

u

yz

yxyy

. 2.12

Sendo, o primeiro e o segundo termos da equação 2.12 parcelas lineares da deformação ( l ),

enquanto o terceiro e quarto termos representam as quantidades não lineares ( nl ).

2.1.2 Energia potencial de deformação do veio

19

A expressão da energia potencial de deformação de um corpo elástico pode ser representada

através da equação 2.13 [19], assim considerando-se as tensões predominantemente presentes

em um veio sob flexão obtém-se a equação 2.14.

dVUV

Tσε2

1, 2.13

V

zyzyyxyxyyyyU 2

1. 2.14

Ao adoptar a hipótese de material isotrópico e linear, para pequenos deslocamentos tem-se

ao aplicar a lei de Hooke a energia potencial de deformação do veio em termos matriciais dada

por [1]:

dyy

wyGS

y

wy

dyy

wyψGS

y

wyψ

dyy

yEI

y

y

y

yψEI

y

yψU

L T

L T

L

xz

0

**

0

**

0

)()(2

1

)()(2

1

)()()()(

2

1

. 2.15

A partir da equação 2.15 é possível obter a matriz de rigidez do sistema.

2.1.3 Efeito do carregamento axial

Uma viga sob carregamento axial tem sua rigidez à flexão alterada. Esta variação de rigidez

se obtém através dos termos não lineares presentes na equação 2.12, termos estes referentes aos

graus de liberdade de deslocamentos transversais [1]. Este efeito ainda se mostra proporcional

ao carregamento aplicado e pode até se tornar negativo se o veio estiver submetido à

compressão. Este fenómeno é muito comum em motores verticais.

2.1.4 Formulação do elemento finito de viga com interpolação quadrática (3 nós)

Nesta secção é realizada uma aproximação por elementos finitos do tipo quadráticos, com o

intuito de se obter a matriz de rigidez, a metodologia para a obtenção da matriz de rigidez é a

proposta em [1]. O elemento finito correspondente possui três nós e quatro graus de liberdade

por nó, sendo 1u , 2u e 3u os graus de liberdade de deslocamento na direcção x nos nós 1, 2 e

3 respectivamente. Analogamente tem-se 1w , 2w e 3w para a direcção z . Os graus de

20

liberdade de rotação para os nós 1, 2 e 3 em relação aos eixos x e z respectivamente são 1 ,

2 , 3 , 1 , 2 e 3 conforme mostra a Figura 6.

w1 1 Z 2

u1 1 u2

2 X

w2 3

3 u3

w3

Y

L

Figura 6 – Representação dos graus de liberdade no elemento finito de viga de 3 nós.

Aproximando as funções deslocamento através de deslocamentos nodais e das funções de

interpolação iN , tem-se:

θNTyNyNyNy 332211 )()()()( , 2.16

ψNTyNyNyNy 332211 )()()()( , 2.17

uNTuyNuyNuyNyu 332211 )()()()(* 2.18

e

wNTwyNwyNwyNyw 332211 )()()()(* . 2.19

sendo θ , ψ , u , w os vectores contendo as rotações e os deslocamentos nodais e N o vector

das funções de interpolação para elemento quadrático, dado por:

L

y

L

y

L

y

L

y

L

yT 21

21

21

2

N . 2.20

Introduzindo as equações 2.16, 2.17, 2.18 e 2.19 na equação 2.15, considerando B a

derivada de N , e manipulando os termos de maneira análoga com a [1] tem-se:

21

wBBwθBNw

wNBθθNNθ

uBBuψBNu

uNBψψNNψ

ψBBψθBBθ

109

87

65

43

21

KK

KK

KK

KK

KK

dyGSdyGS

dyGSdyGS

dyGSdyGS

dyGSdyGS

dyEIdyEIU

LTT

LTT

LTT

LTT

LTT

LTT

LTT

LTT

L

zTT

L

xTT

00

00

00

00

00

2

1

2

1

2

1

2

1

2

1

2

1

2

1

2

1

2

1

2

1

. 2.21

na qual iK para 1,...,10i são as matriz de rigidez relativas aos graus de liberdade a que

estão associadas. Segundo [20], a parcela de rigidez da matriz elementar relativa a flexão é

obtida à partir da integração completa de 1K e 2K e os termos de 3K a 10K relativos ao corte

são sub - integrados à partir de uma integração gaussiana selectiva dando origem à contribuição

do corte na matriz de rigidez elementar [1].

A matriz de rigidez elementar é obtida seguindo a metodologia utilizada por [1], aplica-se as

equações de Lagrange (equação 2.22) em todas as parcelas iK .

iq

iiii

Fqq

U

q

T

q

T

dt

d

, 2.22

sendo T a energia cinética, U a energia potencial de deformação, a função de dissipação de

Rayleigh, iqF os esforços generalizados em relação as coordenadas generalizadas iq e iq as

velocidades generalizadas do sistema. Através da expansão do termo relativo a 1K , que

relaciona apenas o grau de liberdade , e aplicando as equações de Lagrange, com a

consideração de que 121

112 KK ,

131

113 KK e

132

123 KK , tem-se:

1133

1122

1111

1

KKKU

, 2.23

22

1233

1222

1211

2

KKKU

2.24

e

1333

1322

1311

3

KKKU

. 2.25

As equações 2.23, 2.24 e 2.25 representam as relações entre os graus de liberdade 1 , 2 e

3 e os coeficientes de rigidez da matriz 1K . Considerando que os graus de liberdades são

dispostas na matriz elementar como mostrado na equação 2.26, tem-se a equação 2.27.

33331111 wuwuT δ 2.26

0000000

0000

0000000

0000

0000

0000000

0000000

133

132

131

123

122

121

113

112

111

3

3

2

2

1

1

1

3322111

KKK

KKK

KKK

w

uwu

. 2.27

Este procedimento é realizado de maneira análoga para todas as parcelas iK , e assim é

possível obter a matriz elementar do sistema (APÊNDICE A).

Contribuição dos carregamentos externos na rigidez do sistema

Os modelos de elementos finitos para os efeitos de alteração da rigidez do sistema através

de carregamentos axial e magnético podem analogamente ser deduzidos como em [1] para o

elemento de viga de Timoshenko com 3 nós. As matrizes que representam estes efeitos podem

ser visualizadas no (APENDICE B).

2.2 Obtenção da Matriz de Massa e Giroscópica

Nesta secção apresenta-se resumidamente a obtenção da matriz de massa e giroscópica, o

procedimento para obtenção destas matrizes é o mesmo apresentado em [1].

23

2.2.1 Composição de rotações

A Figura 7 apresenta a composição de rotações de um disco, o mesmo gira em torno dos

eixos coordenados do sistema inercial e móvel. A sequencia de rotações apresentada na Figura 7

começa por uma rotação em torno do eixo Z , na sequencia uma em torno do eixo 'x e

finalmente em torno de ''y .

Figura 7 – Eixos de referência para um disco que sofre rotação num veio flexível [1].

Assim baseado na Figura 7 e de acordo com [1] obtém-se:

T

T

z

y

x

sentt

sentt

tsen

t

coscos

coscos

10ω . 2.28

na qual 10ω é o vector de rotações com as componentes das rotações nos eixos coordenados

inerciais.

2.2.2 Energia cinética do disco

24

A energia cinética total do disco em termos da sua massa DM e de seus momentos de

inércia DxI , DyI e DzI pode ser escrita como:

z

y

x

Dz

Dy

Dx

zyxD

D

I

I

IvM

T

00

00

00

2

1

2

2

, 2.29

sendo que a velocidade de translação é dada por

22

2

t

w

t

uv . 2.30

Assim assumindo pequenos deslocamentos e simetria do veio, segundo [1] pode-se chegar

a:

tI

ttI

t

w

t

uMT

Dy

Dx

DD

22

1

2

1

2

1

2

22

22

. 2.31

Na equação 2.31 nota-se a presença do efeito giroscópico representado pelo termo

proporcional à velocidade e à coordenada generalizada , o qual é dado por 22

1DyI .

Esta parcela da energia cinética é responsável pela geração da matriz giroscópica. Note-se que

neste trabalho são utilizadas as duas notações para velocidade, ou seja, derivada do

deslocamento em relação ao tempo, faz-se isso com o objectivo de facilitar a escrita das

equações e para colocá-las neste trabalho da maneira mais conveniente.

2.2.3 Formulação de elementos finitos para o disco

As matrizes de massa e giroscópica do disco, que serão posteriormente adicionadas às

matrizes elementares do sistema, são obtidas através da aplicação das equações de Lagrange

equação 2.22.

Assim tem-se:

25

w

u

I

I

w

u

I

I

M

M

Dy

Dy

Dx

Dx

D

D

agiroscópic Matrizmassa de Matriz

000

000

0000

0000

000

000

000

000

. 2.32

2.2.4 Energia cinética do veio

A energia cinética para o veio é obtida através da integração dos infinitos elementos

diferencias nos quais se divide o veio [1], assim integrando ao longo do comprimento pode-se

obter:

dyt

IIL

ytt

Iyt

w

t

uST

V

VV

L

V

LL

S

0

2

0

22

0

22

2

2

1

2

1

. 2.33

sendo a massa específica do veio, VL o comprimento total do veio, S a área da secção

transversal, I o momento de segunda ordem de área .

2.2.5 Formulação do elemento finito com interpolação quadrática (3 nós)

Para obter as matrizes de massa é giroscópica uma aproximação através de elementos finitos

com funções de forma quadráticas é apresentada, seguindo a mesma metodologia apresentada

por [1]. O elemento finito da Figura 6 é utilizado também para a dedução nesta secção. As

funções de velocidade são aproximadas através de velocidades nodais e funções de interpolação

iN , semelhantes as funções utilizadas para os deslocamentos, mostrados nas equações 2.16,

2.17, 2.18 e 2.19.

Se θ ,ψ , u e w são os vectores das velocidades nodais a equação 2.33 pode ser escrita na

forma,

26

)(

0)(

2

)(

0

)(

0

)(

0

)(

0

2

2

1

2

1

2

1

2

1

vi

LTTT

v

V

iv

LTTT

iii

LTTT

ii

LTTT

i

LTTT

S

V

VV

VV

dyIIL

dySdyS

dySdyST

θNψN

ψNψNθNθN

wNwNuNuN

, 2.34

e analisando os termos )(i , )(ii , )(iii e )(iv da equação 2.34, correspondentes à matriz de

massa tem-se:

ψNNψθNNθ

wNNwuNNu

43

21

00

00

2

1

2

1

2

1

2

1

M

LTT

M

LTT

M

LTT

M

LTT

VV

VV

dyIdyI

dySdyS

, 2.35

A aplicando as equações de Lagrange obtém-se a matriz de massa para o veio (APÊNDICE

C). Os termos da matriz de massa são organizados na matriz seguindo a ordem do vector de

acelerações nodais

][ 333322221111 wuwuwu , 2.36

Para a análise do efeito giroscópico, o termo )(vi da equação 2.34 é desenvolvido.

Considerando a aproximação da função velocidade para o grau de liberdade e a função

deslocamento para o grau de liberdade , tem-se:

θNNψ

θNψNθψ

G1

LTT

LTTT

LT

dyI

dyIdyyyI

0

00

2

22

. 2.37

Com G1 obtém-se a matriz anti-simétrica proporcional a e a , chamada de matriz

giroscópica, proporcional a velocidade de rotação. Através da manipulação da equação 4.37 e

27

aplicação das equações de Lagrange ao problema, da mesma maneira como apresentado para

um elemento de dois nós em [1], para o elemento de três a matriz ordenada segundo a mesma

disposição dos graus de acelerações nodais, porém para velocidades nodais fica da forma

apresentada no APÊNDICE D.

2.3 Modelo para Mancais e Fundação

Neste trabalho utiliza-se o modelo para mancais e fundação proposto por [1] que considera o

mancal apenas como uma massa concentrada, com rigidez e amortecimento associados e que

são inseridas no respectivo nó com graus de liberdade de deslocamento u e w sem restrições

aos graus de liberdade e , a matriz de rigidez, amortecimento e massa relativas aos

mancais e que são adicionadas às matrizes globais do sistema são mostradas respectivamente

pelas equações 2.38, 2.39 e 2.40.

w

u

kk

kkwu

K

zzzx

xzxxmi

. 2.38

w

u

cc

ccwu

C

zzzx

xzxxmi

. 2.39

w

u

m

mwu

M

mi

mimi

0

0 . 2.40

Os termos da matriz de rigidez xxk , zzk , xzk e zxk representam a rigidez directa nas

direcções x e z e os termos cruzados para as mesmas direcções. Da mesma forma são

formadas as matrizes de amortecimento, com os termos directos e cruzados xxc , zzc , xzc e zxc ,

na respectivas direcções x e z . Por fim, mim representa a massa do mancal.

Considerando mais um nó acoplado em cada mancal e apenas com graus de liberdade de

deslocamento a fundação é modelada, seguindo a proposta de [21] e também mostrada em [1].

Na Figura 8, 1mk representa a rigidez do primeiro mancal e 2mk a rigidez do segundo

mancal. 1fk e 2fk são os coeficientes de rigidez das respectivas fundações. Embora a Figura 8

28

apresente apenas dois mancais e duas fundações, o modelo proposto pode introduzir n mancais

e fundações [1].

Rotor

Mancais

Fundação

Nó 1

Nó n+1 Nó n+2

Nó n

Nó n+4

x

y z Nó n+3

1mk

1fm

1fk

1mc

1fc

2mk

2fm

2fk

2mc

2fc

Figura 8 - Representação do modelo simplificado de fundação utilizado, [1].

A rigidez dos mancais é compartilhada por dois nós, 11 nenó . Enquanto a rigidez da

fundação está sujeita ao deslocamento relacionado com o nó 1n pois o nó 3n está

engastado. Então, tomando-se apenas os graus de liberdade de deslocamento u , a matriz de

rigidez associada aos mancais e fundação é dada por:

1

1

111

11

11

nfumumu

mumu

n

ad

u

u

kkk

kkuu

K , 2.41

sendo 1muk e 1fuk a rigidez do mancal 1 e da fundação 1 na direcção do deslocamento u ,

respectivamente. Para w tem-se que:

1

1

111

11

11

nfwmwmw

mwmw

n

ad

w

w

kkk

kkww

K . 2.42

A rigidez do mancal 1 e da fundação 1 na direcção do deslocamento w são,

respectivamente, 1mwk e 1fwk . Para os deslocamentos u e w , a ordem da matriz global de

rigidez será incrementada em duas vezes o número de fundações, conforme mostrado pela

matriz da equação 2.43, para um sistema com dois mancais e duas fundações [1].

29

44222

222

111

111

2

22

2

11

11

2

2

1

1

1

1

221111

000000

000000

000000

000000

0000

000

000

000

nnfwmwmw

fumumu

fwmwmw

fumumu

mw

mwmu

mu

mwmw

mumu

n

n

n

n

n

n

nnnnnn

kkk

kkk

kkk

kkk

k

kk

k

kk

kk

w

u

w

u

w

u

w

uwuwuwuwu

K

. 2.43

Para amortecimento intrínseco do veio nulo, admite-se amortecimento apenas nos graus de

liberdade referentes aos mancais e à fundação. Na equação 2.44 é apresentada a matriz de

amortecimento relacionada com o veio, mancais e fundação [1].

44222

222

111

111

2

22

2

11

11

2

2

1

1

1

1

221111

000000

000000

000000

000000

0000

000

000

000

nnfwmwmw

fumumu

fwmwmw

fumumu

mw

mwmu

mu

mwmw

mumu

n

n

n

n

n

n

nnnnnn

ccc

ccc

ccc

ccc

c

cc

c

cc

cc

w

u

w

u

w

u

w

uwuwuwuwu

C

, 2.44

Os coeficientes 1muc e 1mwc são relativos ao amortecimento dos mancais, 1fuc e 1fwc são os

coeficientes de amortecimento da fundação, nas direcções dos deslocamento u e

w respectivamente.

Para a matriz de massa, é inserida a massa da fundação de acordo com os nós do sistema e

consequentemente aumenta-se a dimensão da matriz [1]. 1fm relativa à massa da fundação sob

o primeiro mancal e 2fm para o segundo mancal. Para os graus de liberdade relativos ao eixo

de coordenadas x , a matriz de massa relativa à fundação do primeiro mancal e que será

adicionada na matriz global é dada por:

1

1

1

11

0

00

nf

n

ad

u

u

m

uu

M . 2.45

Então a matriz de massa global para a estrutura em estudo fica na forma [1]

30

442

2

1

1

2

2

1

1

1

1

221111

0000000

0000000

0000000

0000000

0

0000A

0000S

S

0000A

0000M

nnf

f

f

f

n

n

n

n

n

n

nnnn

m

m

m

m

w

u

w

u

w

u

w

uwuwuwu

M

, 2.46

2.4 Modelos de interferência entre pacote de chapas e veio

O modelo de utilizado para considerar as alterações das propriedades do rotor após a

inserção do pacote de chapas no veio com ajuste de interferência é proposto por [1]. Este

modelo considera propriedades equivalentes no escalonamento que contém o pacote de chapas,

essas propriedades de massa e elasticidade são então quantificadas pelas equações 2.47, 2.48,

2.49 e 2.50. Na qual o módulo de Young é representado por E , a densidade de massa por , o

momento de inércia de área por I , a área de secção transversal por S e módulo de elasticidade

transversal por G .

discodiscoveioveioeq SSS . 2.47

discodiscoveioveioeq III . 2.48

discodiscoveioveioeq IEIEEI . 2.49

discodiscoveioveioeq SGSGGS . 2.50

31

3. METODOLOGIA DE CÁLCULO

Os parâmetros modais do rotor são obtidos através de uma análise que consiste basicamente

no cálculo de valores e vectores próprios. Através destes dados pode-se determinar o

comportamento dinâmico do rotor ao longo da frequência de rotação.

Após a determinação das matrizes elementares pela aplicação das equações de Lagrange e

método dos elementos finitos, secções 2.1 e 2.2, pode-se através de uma assemblagem obter as

matrizes globais do sistema, as quais são função da discretização e das propriedades do sistema,

tais como geometria e propriedades dos materiais. Assim, através das matrizes obtidas tem-se o

sistema de equações global:

tttt fqKqGCqM . 3.1

Nesta equação q é o vector das coordenadas generalizadas, com o mínimo de coordenadas

necessárias para descrever a posição de um corpo no espaço. Estas coordenadas correspondem a

todos os graus de liberdade do sistema. E as forças generalizados são representados por )(tf .

O amortecimento intrínseco do material metálico que compõem o rotor é desprezado, porém

o amortecimento dos mancais na maioria dos modelos é considerado, e um modelo através de

amortecimento proporcional não é adequado para mancais de filme de óleo que possui

coeficientes muito elevados, então por ser considerado um modelo de amortecimento do tipo

viscoso, uma mudança de variáveis do espaço de configurações para o espaço de estado é

realizada, para diminuir a ordem das equações diferenciais e possibilitar a solução de um

problema de valores e vectores próprios.

A solução do problema de valores e vectores próprios é mais detalhada em [1]. Os

resultados do problema de valores e vectores próprios, são as frequências naturais do problema

e as respectivas formas de vibrar. Neste tipo de problema existe a matriz giroscópica, que é

dependente da velocidade de rotação, então para que seja possível entender o comportamento da

máquina em toda faixa de operação é construído o chamado diagrama de Campbell, e ainda

considerando os vectores próprios, pode-se ainda obter os modos de flexão do rotor [1].

Passando o sistema do espaço de estado para um espaço modal do espaço de estado pode-se

obter a resposta do sistema, outra ferramenta importante para análise de vibrações em sistemas

rotativos [1].

32

4. SOLUÇÕES PARA CONTROLO DE VIBRAÇÃO

Poucos são os métodos óptimos que podem ser utilizados no controle de vibração de

sistemas girantes devido a sua complexidade de projecto. Os autores das referências, [10] e [11]

mostram bons resultados de pesquisas desenvolvidas utilizando absorsores óptimos de vibração,

porém neste trabalho, como citado anteriormente o projecto do absorsor não é realizado com

base em uma rotina de optimização. O autor da referência [10] assim como seus orientadores,

pioneiros na aplicação de absorsores dinâmicos viscoelásticos de vibração para rotores,

colocaram esta proposta no registro de patentes. Portanto para qualquer aplicação deste

componente em rotores os autores devem ser devidamente consultados. O balanceamento

teórico óptimo mostrado por [9] e mancais de filme de óleo (referencia [14]) também são

alternativas aqui estudadas.

4.1 Absorsores de Vibração

Absorsores de vibração são dispositivos adicionados ao sistema, com o objectivo de

dissiparem energia se possuírem elemento absorsor, e absorverem vibração do mesmo. O

princípio básico adoptado no projecto de absorsores é que este sistema secundário acoplado

possua a mesma frequência natural, que se deseja neutralizar no sistema primário (Figura 9).

Assim uma simplificação corriqueira pode ser adoptada, como mostrada na Figura 9. O sistema

primário é representado como sendo o conjunto massa mola amortecedor, 1, e o absorsor é

representado pelo conjunto massa mola amortecedor 2.

Figura 9 – Representação do absorsor, sistema secundário, acoplado ao sistema primário.

Para acrescentar ao modelo a influência dos absorsores pode-se utilizar um modelo de

parâmetros equivalentes generalizados, como é mostrado em [10]. Este modelo também é

aplicado para modelar sistemas apoiados em fundação flexível [21].

m2

m1

k2 c2

k1 c1

x2

x1

33

Outra metodologia que pode ser utilizada é o acréscimo de mais graus de liberdade ao

sistema. Este modelo é mais simples, pois possui equações que não dependem da frequência,

excepto para absorsores viscoelásticos. Aplicando-se a segunda lei de Newton ao sistema

representado na Figura 9, pode-se chegar as matrizes que se deve adicionar ao sistema de

equações original,

122

221

1

n

i

uu

uuuni

u

ad

u

u

kk

kkkuu

K , 4.1

tal que uadK é a matriz elementar de rigidez que deve ser adicionada ao sistema devido à

presença do absorsor, 1k é a rigidez do sistema primário e 2k é a rigidez do absorsor. Nota-se

que a matriz deve ser adicionada na posição do i-ésimo nó, respectivo ao grau de liberdade de

deslocamento, neste caso u , grau de liberdade relativo ao deslocamento horizontal. A mesma

análise para o grau de liberdade de deslocamento na direcção vertical pode ser feita, logo tem-se

a seguinte relação:

122

221

1

n

i

uu

uuuni

abs

ad

w

w

kk

kkkww

K . 4.2

As matrizes de amortecimento (equações 4.3 e 4.4) são semelhantes às matrizes de rigidez,

porém posicionadas de acordo com as velocidades nodais.

122

221

1

n

i

uu

uuuni

u

ad

u

u

cc

cccuu

C 4.3

e

122

221

1

n

i

ww

wwwni

w

ad

w

w

cc

cccww

C 4.4

A matriz de massa adicionada ao sistema, devido ao absorsor, é mostrada na equação 4.5, a

mesma é obtida a partir dos termos de massa proporcionais às acelerações nodais.

34

12

1

1

0

0

n

i

u

uni

u

ad

u

u

m

muu

M . 4.5

Analogamente para o grau de liberdade w tem-se:

12

1

1

0

0

n

i

w

wni

w

ad

w

w

m

mww

M . 4.6

Assim as matrizes globais do sistema com o aumento de graus de liberdade para o

acréscimo de um absorsor no sistema são representadas nas equações 4.7, 4.8 e 4.9. Nota-se que

o número de linhas e colunas acrescentadas é proporcional ao número de absorsores acoplados

ao sistema.

2222

22

221

221

1

1

11

000

000

00

0

0

nnn

n

n

n

i

i

nnnnii

kk

kk

kkk

kkk

w

u

w

u

w

u

wuwuwu

K. 4.7

35

2222

22

221

221

1

1

11

000

000

00

0

0

nnn

n

n

n

i

i

nnnnii

cc

cc

ccc

ccc

w

u

w

u

w

u

wuwuwu

C. 4.8

222

2

1

1

1

1

11

0000

0000

00

00

00

nnn

n

n

n

i

i

nnnnii

m

m

m

m

w

u

w

u

w

u

wuwuwu

M. 4.9

Quando perfeitamente sintonizados, os absorsores do tipo massa mola absorvem totalmente

a amplitude de vibração para a frequência natural do sistema primário a qual se deseja

neutralizar, no entanto duas novas frequências naturais são geradas, antes e depois da frequência

neutralizada, este comportamento é mostrado para o sistema da Figura 10, sem o absorsor, e

com o absorsor massa mola na Figura 11, na qual se pode observar que na ressonância o sistema

apresentará grande amplitude de vibração. Para absorsores que possuem o elemento dissipador

de energia os picos serão amortecidos, o que é uma grande vantagem na utilização deste tipo de

componente, evitando falhas por fadiga e atenuando fortemente a vibração [11].

36

Figura 10 – Sistema primário sem absorsor.

Figura 11 – Sistema primário com absorsor massa mola.

Os materiais viscoelásticos são muito empregados na confecção de absorsores de vibração.

Isso se deve ao comportamento elástico - viscoso destes materiais dissipar muita energia, assim

esses tipos de absorsores quando adicionados ao sistema geram uma curva semelhante a Figura

12.

1 n 2

0

1x

n

0

1x

37

Figura 12 – Comparação entre os sistemas com absorsor sem amortecimento - e com amortecimento - -.

Os materiais viscoelásticos têm propriedades muito dependentes da frequência e da

temperatura, como apresentado na Figura 13, sendo relativamente complexo modelar seu

comportamento matematicamente.

Figura 13 – Nomograma de material viscoelástico 242. (Ref. 3M Ultra-Pure Viscoelastic

Damping Polymer)

No caso de absorsores viscoelásticos as matrizes representadas pelas equações 4.1, 4.2, 4.3,

4.4, 4.7 e 4.8, têm coeficientes dependentes da frequência. O cálculo dos parâmetros modais e

do diagrama de Campbell pode então ser obtido através da metodologia apresentada por [22]. O

método consiste na solução de um problema de valores e vectores próprios interno (Figura 14 -

(b)). Esse procedimento de cálculo deve ser executado, pois no problema se tem a matriz

giroscópica que é função da frequência de rotação rpm e para cada rotação 1 se tem novas

matrizes de rigidez e amortecimento função da frequência . Assim o diagrama de Campbell

1

n 2 0

1x

38

interno consiste na determinação das frequências naturais do sistema em ao longo da frequência

através de um loop interno, ou seja dentro do loop para frequência de rotação rpm . (Figura

14).

Figura 14: Diagramas de Campbell e Campbell interno, [22] .

Baseado nestes procedimentos de cálculo e num modelo matemático do material

viscoelástico desenvolvido em [10] uma metodologia para projecto óptimo de absorsores

viscoelásticos para rotores. A viscoelasticidade, que determina a lei de variação dos coeficientes

de amortecimento e rigidez do material em relação à frequência, é modelada matematicamente

através de derivadas fraccionárias. Os absorsores são adicionados ao sistema através de

parâmetros equivalentes generalizados a fim de construir uma função objectivo a ser

optimizada.

O projecto de um absorsor é apresentado neste trabalho usando uma metodologia simples

sem optimização de parâmetros como rigidez e massa do absorsor. Valores de massa são

adoptados para um valor adequado para que o absorsor não perca muita eficiência, e a posição

do componente no veio é escolhida através da análise dos modos de vibrar gerados pelo código

numérico e de acordo com as limitações impostas pelos outros componentes do motor. Com o

valor de massa adoptado a espessura da manta de borracha é calculada, mantendo o valor da

frequência natural do absorsor igual à frequência para qual que se deseja neutralizar a vibração.

Na Figura 15 pode ser observado um esquema de um absorsor de vibração composto por uma

39

massa e um material viscoelástico com espessura h e quando sujeito a uma força F sofre um

deslocamento x devido a um deslocamento angular .

Figura 15: Absorsor com material viscoelástico.

Como materiais viscoelásticos dissipam mais energia quando sujeita ao corte, a espessura do

material é calculada baseado no equilíbrio de forças aplicado à Figura 15. Logo:

hkxkF sen11 . 4.10

Considerando G o módulo de elasticidade transversal, 1k a rigidez do absorsor e a A

área de corte. Então para pequenas deformações pode-se obter a seguinte relação:

Ak

Gh

1

. 4.11

A rigidez do componente é determinada considerando o absorsor como um sistema de um

grau de liberdade, então:

12

1 mk n . 4.12

4.2 Equilíbrio óptimo

Diz-se que uma estrutura está desequilibrada quando seu centro de massa (linha) não

coincide com o eixo de rotação do sistema. Este fenómeno é resultante da não homogeneidade

do material (porosidade, inclusões e outros), imperfeições do fabrico e montagem.

As massas desequilibradas geram forças centrífugas, perpendiculares ao eixo do rotor. Esta

força harmónica excita o sistema, e o seu comportamento pode ser observado numa análise da

chamada Resposta Fundamental do Rotor ao Desequilíbrio.

h

x

γ

F

40

O objectivo da equilibragem é reduzir a magnitude de vibração para valores não maiores do

que os aceitáveis.

Normas internacionais classificam rotores de acordo com os seus requisitos de equilibragem

e estabelecem métodos para avaliação do desequilíbrio residual. De acordo com a norma ISO

11342:1998, em algumas situações nas quais o amortecimento é muito pequeno, as amplitudes

de vibração são muito elevadas se a máquina opera próxima ou na condição de ressonância.

Nestes casos pode ser mais vantajoso aumentar o amortecimento ou modificar a frequência de

ressonância, do que efectuar um processo de equilibragem.

A norma ainda sugere dois tipos principais de equilibragem: a equilibragem em multi-

planos, que é qualquer procedimento de equilibragem, para rotores rígidos ou flexíveis, que

envolva mais do que dois planos de correcção, ou também a equilibragem modal, descrita pelo

processo de equilíbrio de rotores flexíveis, no qual as correcções são feitas para modos de vibrar

distintos.

A distribuição do desequilíbrio é função do processo de fabrico e montagem do rotor

podendo variar consideravelmente mesmo para rotores idênticos. A distribuição do

desequilíbrio possui maior influência para rotores flexíveis do que para rotores rígidos porque

esta distribuição determina o grau de intensidade com que um determinado modo será excitado.

A equilibragem em rotores flexíveis é geralmente mais complexa em relação a equilibragem de

rotores rígidos, devido às elevadas amplitudes de vibração e altas rotações nominais.

Uma metodologia para a equilibragem de rotores flexíveis é apresentada. A técnica é

semelhante ao procedimento apresentado por [9] diferenciando-se pelo método numérico

utilizado e pelo tipo de algoritmo genético. A abordagem nesta secção é muito diferente quando

comparada com a secção 4.1, pois nesta secção não serão adicionados dispositivos ao sistema a

fim de controlar a vibração, mas sim a intervenção é realizada directamente na força que actua

no sistema, neste caso o desequilíbrio (Figura 16).

Figura 16: Desequilíbrio do rotor.

e

F(Ω)

41

Na equação 4.13 observa-se, considerando um desequilíbrio com massa dm , excentricidade

e e fase , a força de excitação resulta em:

irpmd em exp2F . 4.13

Pode-se assim calcular a resposta em frequência em termos de amplitudes complexas. Com

estas amplitudes pode-se construir uma função objectivo a ser optimizada, equação 6.14.

12

12

12

11

11

12

12

12

11

11

nmx

nm

nm

nmx

nm

nm

p

HEw

HEu

HEw

HEu

HEw

HEu

Hw

Hu

Hw

Hu

Hw

Hu

V

. 4.14

As respostas complexas calculadas para as direcções horizontal ( Hu ) e vertical ( Hv ), para

m rotações e n pontos no veio podem então ser comparadas com medições experimentais,

HEu e HEw . A minimização desta adição resulta na correcção para o desequilíbrio, logo

através de um algoritmo de optimização, por exemplo algoritmos genéticos, pode ser encontrado

o desequilíbrio óptimo e através destes valores determinar a correcção para o rotor.

4.2.1 Algorítmos Genéticos

A técnica de optimização não linear chamada algoritmos genéticos segue o mesmo princípio

da evolução de Darwin, em que os indivíduos mais aptos gerados têm maior probabilidade de

prosseguir ao longo das gerações. Diversos são os parâmetros utilizados nos códigos de

optimização através de algoritmos genéticos, e segundo [23] não existem valores específicos ou

óptimos para todos os problemas, em cada caso deve-se escolhe-los de acordo com a

experiência.

Os operadores geralmente introduzidos num código de algoritmos genéticos são: operadores

de selecção, que direccionam a escolha para melhores regiões do espaço de busca dando

preferência aos indivíduos mais aptos; operadores de cruzamento, que são responsáveis por

manter as características dos indivíduos passados e operadores de mutação que restabelecem a

diversidade genética eventualmente perdida durante os processos de cruzamento. Para melhorar

a eficiência dos algoritmos pode-se também escolher adequadamente o tamanho da população,

42

percentual de mutação, percentual de cruzamento, número de casais, número de gerações,

pressão, dentre outros. Alguns destes parâmetros de entrada são descritos em [23].

A Figura 17 representa um processo de optimização com algoritmos genéticos. Segundo [9]

o primeiro passo numa optimização é a etapa de codificação, na qual se descreve o problema de

uma forma matemática e estabelece os domínios para a execução do código. Após a

inicialização das colónias de forma aleatória, começa o mecanismo da selecção. Este

procedimento se baseia nos conceitos da biologia de que os indivíduos mais aptos têm maior

probabilidade de sobrevivência. Posteriormente são realizadas as operações genéticas nos

indivíduos seleccionados previamente, estes indivíduos são então cruzados entre si e neste

cruzamento é adicionada uma taxa de mutação. Após a geração dos indivíduos a função é

calculada, e de acordo com o número de gerações, que também pode ser entendido como o

número de vezes ou iterações que se deseja executar o algoritmo, a população óptima é

mostrada ou o processo se repete com a população inicial sendo a população resultante do

procedimento anterior.

Figura 17: Sequência lógica de execução de um código de algoritmos genéticos.

Diferentemente de [9], utiliza-se neste trabalho uma classe de algoritmos genéticos

chamados elitistas. Tais algoritmos são chamados desta maneira por trabalharem apenas com as

melhores populações assim as convergências ao mínimo global são mais prováveis quando

comparadas aos algoritmos genéticos comuns.

Pode-se observar ainda na equação 4.14 que esta técnica é capaz de considerar quantos

planos de correcção forem necessários, pode-se também considerar nos cálculos várias rotações

e pontos de medição, assim quanto mais completo é o fornecimento dos dados, melhores são os

resultados da equilibragem teórica.

43

4.3 Apoios

Os apoios além de obviamente sustentar o sistema rotativo, podem ser trabalhados com o

intuito de reduzir as vibrações do sistema. Considerando os rotores usualmente fabricados,

pode-se concluir que o amortecimento do sistema é quase todo inserido pelos mancais,

desconsiderando-se os efeitos de amortecimento provenientes da “almofada de ar” no espaço

entre ferro e o atrito entre as chapas da massa rotórica que podem estar presentes. Rotores que

possuem mancais de rolamento, são passíveis de altas amplitudes de vibração ao passar por uma

velocidade de rotação crítica, devido ao baixo amortecimento do conjunto. No trabalho

apresentado por [22] são acrescentadas sob os apoios mantas de material viscoelástico

caracterizadas através de um modelo de derivadas fraccionárias. O diagrama de Campbell é

obtido este modelo poder utilizado como ferramenta de controlo de vibração, pois reduz as

amplitudes de resposta do rotor devido ao elevado amortecimento introduzido no sistema.

Outro modelo de mancais existente, e que é capaz de introduzir um elevado amortecimento

no sistema é o chamado mancal de filme de óleo. Este tipo de apoio sustenta o veio através de

um filme fino de óleo, com espessura geralmente da ordem de um para mil do diâmetro [14]. Os

mancais hidrodinâmicos, assim chamados, tem características muito distintas em relação aos

outros tipos de mancais existentes, suas propriedades de rigidez e amortecimento são altamente

anisotrópicas e dependentes da velocidade de rotação da máquina, além de que seus coeficientes

de amortecimento são extremamente elevados.

O comportamento deste tipo de apoio pode ser estimado utilizando um estudo aprofundado

sobre o escoamento do fluido que está contido no interior do mancal. A Figura 18 mostra a

secção de um mancal de filme de óleo, uma força F representa o carregamento radial no veio na

secção que está apoiada no mancal, representa a rotação do veio, e é a excentricidade entre

o veio e o centro do mancal, a linha tracejada roxa representa a distribuição de pressão no

mancal, minh é a distância mínima entre o veio e o mancal, é o ângulo entre o eixo vertical e

o ponto de espessura mínima de fluido, maxP é o ângulo entre o ponto de máxima pressão e o

eixo vertical, maxP é a pressão máxima e 0P é o ângulo entre o eixo vertical e o ponto onde a

pressão é nula.

44

Figura 18: Mancal de filme de óleo.

Geralmente o estudo se baseia na determinação do perfil de pressão interno no mancal

através da solução da equação de Reynolds, que considerada para regime permanente e

espessura de filme constante na direcção axial do veio, segundo [12] é dada por:

hURz

ph

zx

ph

R

1

61 33

2. 4.15

Duas soluções usuais são conhecidas para esta equação. Para mancais curtos, ou seja,

relações 1/ DL , [24] propôs a solução dada pela equação 6.16:

3

22

2 )(1

)(3

4

cos

seny

L

rC

Up

r

. 4.16

Na qual é a viscosidade dinâmica do óleo, U é a velocidade tangencial do veio, L é o

comprimento do mancal e y é a coordenada axial.

A solução proposta por Sommerfeld e descrita em [25] considera uma relação 1/ DL ,

assim uma boa aproximação proposta é a consideração de um perfil de pressão constante ao

longo de y , logo, o segundo termo da equação 6.16 pode ser desprezado. Assim tem-se a

solução, conhecida como solução para mancais longos:

222 )(12

)()(26

cos

sencos

C

Urp

r

. 4.17

45

Esta solução para mancais longos é a utilizada pela rotina de cálculos SBCALC 3.0, fornecida

pela RENK. Os coeficientes de rigidez e amortecimento são estimados na rotina através de uma

relação entre o carregamento radial, a velocidade de rotação, a folga diametral e o número de

Sommerfeld que contém a influência da pressão.

Assim como apresentado por [1] para mancais simples, em para mancais de filme de óleo os

graus de liberdade associados aos mancais são apenas deslocamentos u e w . Então, as

matrizes de rigidez ( miK ) e amortecimento ( miC ) adicional devidas ao mancal posicionado no i-

ésimo nó podem ser escritas como apresentado pelas equações 4.18 e 4.19.

i

i

rpmzzrpmzx

rpmxzrpmxx

ii

mi

w

u

kk

kkwu

K

)()(

)()( . 4.18

w

u

cc

ccwu

C

rpmzzrpmzx

rpmxzrpmxxmi

)()(

)()( . 4.19

Os termos da matriz de rigidez )( rpmxxk , )( rpmzzk , )( rpmxzk e )( rpmzxk

representam a rigidez directa nas direcções x e z e os termos cruzados para as mesmas

direcções. Da mesma forma são formadas as matrizes de amortecimento, com os termos directos

e cruzados )( rpmxxc , )( rpmzzc , )( rpmxzc e )( rpmzxc , na respectivas direcções x e z .

A massa dos mancais pode ser incluída no modelo através da consideração de fundação, e então

aos nós em que estão posicionados os mancais são adicionados dois novos graus de liberdade

para cada novo nó.

Como as propriedades dos mancais são dependentes da frequência de rotação, o problema de

valores e vectores próprios deve ser calculado a cada actualização das matrizes de rigidez e

amortecimento.

46

5. RESULTADOS E COMPARAÇÕES

Nesta secção serão apresentados os resultados obtidos com a formulação de elementos

finitos proposta. A validação do modelo consiste na comparação com resultados de um

fabricante de mancais de filme de óleo, assim como uma comparação com resultados

experimentais na forma de curvas de resposta em frequência.

As propostas de controlo são também testadas e seus resultados apresentados. O

balanceamento teórico é comparado com medições experimentais, e sua atenuação na vibração

do sistema é medida de forma numérica. Simulações numéricas com diferentes parâmetros de

mancais hidrodinâmicos são realizadas com o objectivo de avaliar se é possível utilizar mancais

de filme de óleo para controlo de vibração do sistema. É apresentado ainda o absorsor

projectado segundo os parâmetros modais identificados do sistema (motor BFN6 500 H2 SP14),

porém não foi possível obter resultados devido ao curto tempo disponível para fabricação e

implementação numérica do modelo numérico do dispositivo.

O cálculo do comportamento do sistema e seus resultados são obtidos através de uma rotina

de cálculos foi desenvolvida a partir do código computacional existente ROTORDYN. O código

original foi desenvolvido através de projectos e termos de cooperação estabelecidos entre a

empresa WEG Indústrias Eléctricas e o Laboratório de Vibrações da Universidade Tecnológica

Federal do Paraná. Este código foi totalmente transcrito para que fosse possível sua execução

em um software livre e também foi adaptado para a realidade da WEGeuro.

5.1 Interface Gráfica

O trabalho realizado para alterar a rotina de cálculo ROTODYN, consistiu em duas etapas.

Situações corriqueiras, como mancais de filme de óleo, esforço axial, vibração torsional,

enrijecimento devido ao pacote de chapas e excitação externa harmónica, foram implementadas

através de funções auxiliares e assim colocada junto às funções pré existentes. A segunda etapa

consistiu no desenvolvimento de uma interface gráfica para entrada de dados. As

implementações foram todas feitas no software Scilab, pois todo o código teve de ser transcrito

para uma linguagem freeware.

Através de janelas interactivas, o usuário pode de maneira mais didáctica inserir os dados e

executar as simulações. A interface inicial é mostrada na Figura 19. Pode-se observar as funções

para entrada de dados do veio, discos, costelas, mancais com propriedades constantes ou

47

variáveis, fundação e tipo excitação do tipo desequilíbrio ou externa harmónica. Tem-se ainda a

ferramenta para a discretização iterativa e cálculo estático que gera como resultados as reacções

nos apoios e a linha elástica ou deformada do veio. As opções de cálculo são: diagrama de

Campbell, rotações críticas à flexão e torção, além dos modos de vibrar de flexão e torção. A

resposta no tempo e as respostas em frequência para as excitações do tipo desequilíbrio ou

externa harmónica também podem ser calculadas. Nota-se a presença da função que

corresponde a local de acesso aos dados necessários para o cálculo dos parâmetros para o

balanceamento teórico.

Figura 19: Interface ROTORDYN.

5.2 Comparação dos Resultados Numéricos com Dados da RENK

Resultados existentes no estudo do motor BFN6 500 H2 SP14 fornecidos pela fabricante

RENK foram comparados a fim de validar os resultados obtidos com o código computacional

elaborado. O modelo discretizado do veio foi o mesmo presente no estudo do motor realizado

pela RENK, Figura 20. Porém não se tem conhecimento do tipo de elemento utilizado pela

RENK, o elemento utilizado neste trabalho é o proposto nas secções anteriores e seus resultados

são comparados tanto com os resultados da RENK, quanto com os resultados de [1].

48

Figura 20: Veio do motor BFN6 discretizado com 29 elementos finitos, mancais representados por

triângulos azuis, veio por rectângulos brancos, discos por rectângulos verdes e desequilíbrios em vermelho.

Nesta simulação o fabricante de chumaceiras RENK não utiliza qualquer modelo de

aumento de rigidez devido ao pacote de chapas, as propriedades dos mancais são consideradas

para uma frequência de rotação de ][3000 rpm , portanto constantes. Na Figura 21 pode-se

observar a interface da rotina de cálculos para as propriedades dos mancais de filme de óleo.

49

Figura 21: Interface da rotina de cálculos RENK para obter as propriedades dos mancais de filme de óleo.

No modelo numérico do rotor o pacote de chapas é representado por quatro discos

consecutivos para que a massa rotórica seja distribuída de maneira mais uniforme nos nós dos

elementos, e é considerada como fundação elástica do sistema a parte estrutural do mancal de

filme, dados estes que são da mesma forma fornecidos pela RENK. A Tabela 1 apresenta as

dimensões e posicionamento dos quatro discos que representam o pacote de chapas do rotor. Os

discos são considerados como sendo de material homogéneo, com densidade de massa

3/7850 mkg . O veio maciço também possui a mesma densidade dos discos e as

dimensões são apresentada na Tabela 2. As demais propriedades do veio são: módulo de

elasticidade GPaE 210 e coeficiente de Poisson 3.0 .

Tabela 1: Dimensões dos discos que simula o pacote de chapas.

Posição dos discos [mm]

Diâmetro externo [mm]

Diâmetro interno [mm]

Espessura[mm]

982,5 482,6 190 245

1227,5 482,6 190 245

1472,5 482,6 190 245

1717,5 482,6 190 245

Tabela 2: Dimensões do veio do rotor que compõem o motor BFN6.

Posição final do escalonamento no veio [mm] Diâmetro do escalonamento [mm]

170 80

480 100

595 160

1900 190

2114 200

2253 160

2618 100

2678 85

Para efeito de comparação dos resultados, as curvas de resposta em frequência devido ao

desequilíbrio (Tabela 3) são apresentadas. As curvas de respostas representam a resposta em

50

frequência na direcção vertical, calculadas a ][860 mm e ][1840 mm do bordo de ataque do

veio.

Tabela 3: Dados para o desequilíbrio hipotético do rotor.

Local do desequilíbrio 1º Ponto – 860 [mm] 2º Ponto – 1840[mm]

em [g mm] 313 182

][º 0 0

Os modelos de mancais de filme de óleo são os fornecidos pela RENK, EFNLK/Q 9-90,

com lubrificação através de fluido ISO-VG 32, para carregamento de ][26,8 kN . A Tabela 4

apresenta os dados do mancal de filme de óleo.

Tabela 4: Propriedades dos mancais de filme de óleo a 3000 [rpm].

Coeficientes de Rigidez [N/m] Posição [m]

xxk xzk zzk zxk

142525000 34919000 756050000 420066000 0,36

134940000 39250000 782649000 413385000 2,38

Coeficientes de Amortecimento [Ns/m]

xxc xzc zzc zxc

318000 589000 2663000 598000 0,36

288000 565000 2636000 574000 2,38

Na Figura 22 observa-se o diagrama de Campbell obtido pela RENK para o motor BFN6

500 H2 SP14. Para comparar com o diagrama fornecido pela RENK na Figura 23 é mostrado

um diagrama de Campbell calculado através do código desenvolvido neste trabalho. O objectivo

é validar os resultados calculados pelo programa através de dados para o mesmo motor, com a

geometria conforme apresentado na Tabela 1 e na Tabela 2.

51

Figura 22: Diagrama de Campbell, dados do estudo da RENK sobre o motor BNF6 500.

52

Figura 23: Diagrama de Campbell para o rotor em estudo do motor BNF6, com mancais de filme de óleo, teoria de viga de Timoshenko com 3 nós por elemento, funções de interpolação quadráticas.

A curva de resposta em frequência retirada do estudo do motor realizado pela fabricante de

mancais de filme de óleo RENK também é comparada com a curva gerada pela rotina de

cálculos implementada. Na Figura 24 tem-se a curva retirada do estudo do rotor fornecido pela

RENK, a curva de resposta apresenta valores calculados para vários pontos ao longo do veio.

Figura 24: Curva de resposta em frequência para o rotor em estudo, com chumaceira EFNLK, utilizando

a teoria de viga de Timoshenko com 3 nós por elemento, funções de interpolação quadráticas, obtida do estudo da RENK.

Por sua vez, a Figura 25 apresenta a curva de resposta do rotor obtida com o modelo de

elementos finitos proposto por este trabalho, e calculada através do um código numérico

ROTORDYN. Neste gráfico estão presentes as respostas na direcção vertical, nos pontos 2 HA

e 3HA (Figura 24), considerando excitação do tipo desequilíbrio com propriedades mostradas

na Tabela 3.

53

Figura 25: Curva de resposta em frequência para o rotor em estudo, com chumaceira EFNLK, utilizando

a teoria de viga de Timoshenko com 3 nós por elemento, funções de interpolação quadráticas, obtida numericamente.

Com os parâmetros de rigidez e amortecimento dos mancais apresentados na Tabela 4,

calculados pelo programa RENK’s SBCALC 3.0, os resultados fornecidos pela RENK e os

obtidos com o auxílio do código implementado se mostraram muito próximos conforme mostra

mostrados a Tabela 5, as diferenças relativas entre os resultados obtidos com o código numérico

ROTORDYN e dos dados da RENK são calculados por:

100*1

21..

valor

valorvalorreldif

. 5.1

Tabela 5: Comparação entre as rotações críticas.

Resultado RENK ][rpm ROTORDYN ][rpm dif.rel. [%]

1ª Rotação Crítica 1820 1862 2,3

2ª Rotação Crítica 2220 2209 0,5

Resultado RENK ][ m ROTORDYN ][ m

Amplitude – 1ª Rotação Crítica 2,6 3,7

Amplitude – 2ª Rotação Crítica 7 9,7

54

A Tabela 6 apresenta uma comparação entre os resultados obtidos em [1] com os resultados

obtidos neste trabalho. Observa-se uma boa proximidade entre os resultados.

Tabela 6: Comparação elementos com 2 e 3 nós.

Resultado ROTORDYN -

elemento 3 nós ][rpm

ROTORDYN – elemento 2 nós [1]

][rpm dif.rel. [%]

1ª Rotação Crítica 1862 1867 0,3

2ª Rotação Crítica 2209 2217 0,4

Resultado ROTORDYN - elemento 3

nós ][ m ROTORDYN -

elemento 2 nós [1] ][ m

Amplitude – 1ª Rotação Crítica 3,7 3,6

Amplitude – 2ª Rotação Crítica 9,7 9,7

5.3 Modelo de Elementos Finitos Aplicado ao Cálculo dos Parâmetros Modais de Rotores: Validação Experimental

O modelo de elementos finitos proposto e implementado no programa ROTORDYN é

utilizado para o cálculo dos parâmetros modais do rotor que compõem o motor BFN6 500 H2

SP14. Os resultados são obtidos utilizando a discretização apresentada na Figura 26. Deve-se

notar que a discretização mostrada possui mais nós do que estão apresentados na Figura 26,

devido a formulação constituída de elementos finitos quadráticos, logo deve-se considerar que

existe um nó intermediário a cada dois nós (elemento).

55

Figura 26: Discretização do veio.

No diagrama de Campbell apresentado na Figura 27 obtido com a discretização proposta,

foram considerados também parâmetros do mancal de filme de óleo EFNLK variáveis com a

frequência de rotação, e um modelo de interferência para a zona afectada pelo pacote de chapas.

Esta configuração é a que será utilizada para calcular os resultados das alternativas de

controlo propostas.

56

Figura 27: Diagrama de Campbell para o rotor em estudo do motor BNF6, com mancais de filme de óleo,

utilizando a teoria de viga de Timoshenko com 3 nós por elemento, funções de interpolação quadráticas.

Considerando as hipóteses estabelecidas para o cálculo do diagrama da Figura 27, a resposta

do rotor também foi calculada e é mostrada na Figura 28. O desequilíbrio considerado é o

mostrado na Tabela 3 os pontos 1 e 2 de cálculo da resposta são mm860 e mm1840 do

bordo de ataque.

Figura 28: FRF para o rotor em estudo, com chumaceira EFNLK, utilizando a teoria de viga de Timoshenko

com 3 nós por elemento, funções de interpolação quadráticas, obtida numericamente.

57

As linhas de cor azul escuro e vermelho são as respostas nos pontos 1 e 2 respectivamente,

na direcção vertical. As linhas de cor verde e azul clara são as respostas nos pontos 1 e 2

respectivamente, na direcção horizontal.

Na Tabela 7 é possível observar a diferença entre o cálculo numérico e as medições

experimentais, que são apresentadas na secção 5.4.

Tabela 7: Comparação entre as rotações críticas, obtidas numericamente e experimentalmente.

Resultado Experimental

][rpm Numérico

][rpm dif.rel. [%]

1ª Rotação Crítica 2304 2208 4,2

2ª Rotação Crítica 2940 3086 4,9

5.4 Equilibragem Teórica

Utilizando a função de equilibragem teórica através de algoritmos de optimização do código

computacional ROTORDYN, recorrendo ao processo inverso, ou seja, ao invés de obter o

equilíbrio, foi obtido o desequilíbrio e a forma de comparação dos resultados é através de curvas

de resposta em frequência. Os planos de balanceamento são considerados a mm860 e

mm1840 do bordo de ataque do veio. As curvas de resposta em frequência obtidas

experimentalmente para comparação interna no código e posterior optimização foram geradas

com medições em mm250 e mm2479 em relação ao bordo de ataque do veio, para as

direcções horizontais e verticais.

A partir das curvas mostradas nas Figura 29, Figura 32, Figura 35 e Figura 38 são calculadas

amplitudes complexas nas respectivas rotações críticas, conforme apresentadas na Tabela 8.

Essas amplitudes são comparadas internamente pelo código e assim é possível calcular o

equilíbrio (correcção) ou o desequilíbrio do veio, em termos de massa e fase.

58

Tabela 8: Amplitudes complexas.

Ponto (1 - 260 mm, 2 – 2480 mm)

Rotação ][rpm Valor * 10-6

1- horizontal 2304 -1,260+0,337i

1- vertical 2304 -0,600+0,350i

2 – horizontal 2304 +0,630-0,440i

2 – vertical 2304 +0,180-0,180i

1- horizontal 2940 -0,450-0,315i

1- vertical 2940 -0,130+0,350i

2 – horizontal 2940 -0,071+0,071i

2 – vertical 2940 +0,205-0,205i

A Tabela 9 mostra o resultado da optimização a partir das amplitudes complexas calculadas

anteriormente, na sequência, de posse deste resultado, são calculadas curvas de resposta em

frequência para comparação com as curvas experimentais (Figura 30 a Figura 40).

Tabela 9: Resultado da optimização, desequilíbrio.

Plano Massa ][ mkg * 10-6 Fase ][º

1º 170 20,6

2º 147 168

Figura 29: Medição experimental de resposta do rotor ao desequilíbrio, direcção vertical, a 260 [mm] do

bordo de ataque.

59

Figura 30: Resposta numérica ao desequilíbrio, direcção vertical, a 260 [mm] do bordo de ataque,

amplitude.

Figura 31: Resposta numérica ao desequilíbrio, direcção vertical, a 260 [mm] do bordo de ataque, fase.

60

Nas curvas Figura 29, Figura 30 e Figura 31, nota-se uma grande semelhança entre as

curvas numérica e experimental, porém a curva numérica possui amplitudes de resposta

inferiores.

Da mesma maneira que para a direcção vertical, a curva de resposta em frequência

horizontal do mancal que está localizado no bordo de ataque também foi medida (Figura 32).

Quando esta é comparada com a resposta numérica (Figura 33 e Figura 34) nota-se uma boa

semelhança.

Figura 32: Medição experimental de resposta do rotor ao desequilíbrio, direcção horizontal, a 260 [mm]

do bordo de ataque.

61

Figura 33: Resposta numérica ao desequilíbrio, direcção horizontal, a 260 [mm] do bordo de ataque,

amplitude. .

Figura 34: Resposta numérica ao desequilíbrio, direcção horizontal, a 260 [mm] do bordo de ataque, fase.

62

As curvas mostradas na Figura 35, Figura 36 e Figura 37 são respectivamente a medição

experimental no bordo oposto e a amplitude de resposta e fase calculadas numericamente no

mesmo ponto. Nota-se que nesta simulação as curvas apresentam a mesma ordem de grandeza

em relação a amplitude e comportamento ao longo da frequência semelhante, com uma

diferença de comportamento em relação à primeira rotação crítica.

Figura 35: Medição experimental de resposta do rotor ao desequilíbrio, direcção vertical, a 200 [mm] do

bordo oposto ao ataque.

63

Figura 36: Resposta numérica ao desequilíbrio, direcção vertical, a 200 [mm] do bordo oposto ao ataque, amplitude.

Figura 37: Resposta numérica ao desequilíbrio, direcção vertical, a 200 [mm] do bordo oposto ao ataque,

fase.

Finalmente na direcção horizontal tem-se as curvas, experimental Figura 38 e numéricas

Figura 39 e Figura 40, nesta comparação nota-se que não existe uma boa semelhança entre as

curvas como mostrado anteriormente.

64

Figura 38: Medição experimental de resposta do rotor ao desequilíbrio, direcção horizontal, a 200 [mm]

do bordo oposto ao ataque.

Figura 39: Resposta numérica ao desequilíbrio, direcção horizontal, a 200 [mm] do bordo oposto ao

ataque, amplitude.

65

Figura 40: Resposta numérica ao desequilíbrio, direcção horizontal, a 200 [mm] do bordo oposto ao

ataque, fase.

5.5 Projecto de Absorsor de Vibração Viscoelástico

O projecto de um absorsor de vibração baseia-se no sistema que será corrigido. Assim os

parâmetros modais do sistema são de extrema importância para dimensionamento, determinação

da geometria, e selecção dos materiais dissipadores de energia que compõe o absorsor. Os

parâmetros modais do sistema apresentado de forma discretizada pela Figura 26, são calculados

tornando possível o projecto do absorsor.

A partir do diagrama de Campbell e dos modos de vibrar, Figura 27, Figura 41 e Figura 42,

é possível determinar tanto a posição quanto a frequência natural do absorsor.

A melhor eficiência no controle de vibração se obtém com o absorsor totalmente

sintonizado e localizado na posição de maior deslocamento. Em motores eléctricos estas regiões

não estão normalmente disponíveis, logo deve-se obter o melhor material viscoelástico, com

maior factor de perda possível para a temperatura e frequência em questão, a fim de compensar

o facto de não ser possível colocar o absorsor no ponto de maior deslocamento, ponto no qual o

absorsor teria eficiência máxima.

66

Figura 41: 1º Modo de vibrar.

Figura 42: 2º Modo de Vibrar.

Baseado nas curvas expostas anteriormente a ainda no dimensionamento mostrado na

secção 4.1 foi possível dimensionar os componentes do absorsor. A posição escolhida de acordo

com os modos de vibrar das duas primeiras rotações críticas foi 2500 mm do bordo de ataque. A

67

região ideal seria na metade do veio, porém a massa rotórica impede que seja posicionado ali

um absorsor.

A escolha do material viscoelástico é fundamental no projecto do absorvedor, pois a rigidez

e o amortecimento do mesmo depende muito do tipo de material, assim como da temperatura e

frequência, como exemplificado anteriormente. Ainda o material deve possuir uma boa

integridade de suas propriedades, ou seja, estas devem estar sempre de acordo com o

comportamento apresentado pelo seu respectivo nomograma, para não dessintonizarem o

absorsor. As massas tem uma influência característica, delas depende muito a eficiência do

absorsor. Absorvedores com massas muito pequenas são geralmente ineficientes para sistemas

de médio a grande porte.

A geometria do absorvedor é determinada para que o sistema secundário trabalhe em regime

de um grau de liberdade, na frequência a qual se deseja absorver a resposta. No caso de rotores

o absorsor foi projectado para que cada corpo tivesse dois graus de liberdade, porém com as

mesmas frequências naturais de vibração, assim é possível reduzir as amplitudes de vibração nas

duas direcções, vertical e horizontal. No caso do motor BFN6 têm-se duas frequências naturais

na banda de frequência de trabalho, logo o sistema projectado tem quatro corpos com dois graus

de liberdade cada um e duas frequências de trabalho diferentes. O dispositivo foi projectado

para ser posicionado no bordo oposto ao ataque, entre a tampa e o ventilador. Após a montagem

do absorvedor no flange de sustentação tem-se o dispositivo apresentado pela Figura 43.

Figura 43: Absorvedor de vibração viscoelástico acoplado ao flange de sustentação.

Para evitar que o conjunto sofra qualquer rotação de acordo com o veio, o interior do flange

deve abrigar um rolamento de esferas, e ainda em cada bordo do flange um sistema de fixação é

68

posicionado para manter o absorvedor em sua posição estática, ou seja, a flange não deve

rotacionar com o veio.

5.6 Controle de vibração através da selecção da chumaceira

Como os coeficientes de rigidez e amortecimento das chumaceiras influenciam muito nas

características do sistema, os parâmetros modais foram obtidos com a implementação da

metodologia de cálculo que permite considerar as propriedades variáveis dos mancais.

Com dados obtidos através da rotina de cálculo do fabricante de chumaceiras RENK, é

possível, através de uma interpolação cúbica do tipo spline, gerar curvas (Tabela 11 e Tabela

12) para cada coeficiente variante com a frequência de rotação. Esses dados são inseridos no

ciclo de cálculo de valores e vectores próprios, e são actualizados a cada frequência de rotação.

Pode-se observar que as curvas de amortecimento (Tabela 12) tem um padrão de

variação, todas apresentam um decaimento dos valores de amortecimento com o aumento da

frequência de rotação, estas curvas são obtidas através dos dados da Tabela 10, assim como as

curvas para os coeficientes de rigidez. Os coeficientes de rigidez por sua vez não apresentam um

padrão devido à curva do coeficiente xxK que varia de maneira contrária aos outros coeficientes

de rigidez.

Tabela 10: Coeficientes de rigidez e amortecimento para o mancal de filme de óleo EFNLK, 9-90, ISO VG-32, Força radial 8,26[kN].

Rot. [rpm]

Kxx [kN/m]

Kxz [kN/m]

Kzz [kN/m]

Kzx [kN/m]

Cxx [kN/m]

Cxz [kN/m]

Czz [kN/m]

Czx [kN/m]

500 159957 55813 1,03E+06 513548 1927 4090 19784 4159

1000 160775 35337 809083 463587 1101 1977 8786 2006

2000 161192 24907 697549 438138 585 971 4111 985

3000 161125 26574 715374 442205 386 649 2771 658

4000 161065 28074 731420 445867 287 488 2098 495

5000 160973 30384 756113 451501 227 392 1704 398

6000 160895 32346 777096 456289 187 328 1437 333

7000 160806 34565 800831 461704 158 282 1249 286

8000 160730 36455 821043 466316 137 248 1106 251

9000 160663 38146 839122 470441 120 221 993 224

10000 160614 39363 852132 473410 107 199 900 202

69

Tabela 11: Gráficos da variação dos parâmetros de rigidez com a frequência de rotação.

[rpm]

[rpm]

70

Tabela 12: Gráficos da variação dos parâmetros de amortecimento com a frequência de rotação.

71

Baseado nos gráficos da Tabela 11 e Tabela 12 é calculada uma curva de resposta em

frequência (Figura 44), que servirá como referência para comparar a redução das amplitudes de

vibração, com a variação dos parâmetros geométricos e de filme de óleo do mancal. O

desequilíbrio calculado na secção 5.4 é o utilizado como excitação para cálculo das respostas

em frequência. Os pontos de medição estão posicionados a 860 [mm] e 1840 [mm] do bordo de

ataque. O rotor é o mesmo apresentado na Figura 26.

Figura 44: Curva de resposta em frequência para chumaceira sem modificações (referência) Chumaceira -

EFNLK- 9 90 40ºC _iso vg 32, autom, 10610rpm _fr=8,26.

Com a redução da folga, Figura 45, pode-se observar que a amplitude sofre uma grande

atenuação na segunda rotação crítica, porém na primeira rotação crítica a atenuação é menor.

Pode-se observar ainda que a primeira rotação crítica apesar te ter uma leve atenuação em sua

amplitude, a razão de amortecimento é menor, característica que pode ser observada nas curvas.

72

Figura 45: Curva de resposta em frequência para chumaceira - EFNLK- 9 90 40ºC _iso vg 32, autom -1,

10610rpm _fr=8,26.

O efeito na amplitude de vibração ao considerar-se um fluido mais viscoso pode ser

observado na Figura 46, devido ao amortecimento mais elevado pode-se observar uma leve

redução de amplitude, para ambas as rotações críticas.

73

Figura 46: Curva de resposta em frequência para chumaceira - EFNLK- 9 90 40ºC _iso vg 220, autom,

10610rpm _fr=8,26.

Na Figura 47 tem-se o efeito da diminuição do diâmetro na amplitude de vibração, observa-

se que a atenuação é muito pequena comparadas às outras alterações.

74

Figura 47: Curva de resposta em frequência para chumaceira - EFNLK- 9 80 40ºC _iso vg 32, autom,

10610rpm _fr=8,26.

Na Figura 48 o tipo de chumaceira é alterado, nota-se uma forte redução de amplitude para a

segunda rotação crítica, porém para a primeira observa-se um aumento, tanto de amplitude

quanto de frequência.

75

Figura 48: Curva de resposta em frequência para chumaceira - EFZLK- 9 80 40ºC _iso vg 220, autom -1,

10610rpm _fr=8,26.

76

6. CONCLUSÕES

Neste trabalho foi proposto um modelo numérico para análise dinâmica de rotores baseado

no método dos elementos finitos e na teoria de viga de Timoshenko, sendo que o elemento

proposto possui três nós e assim funções de forma quadráticas. Os efeitos giroscópicos e de

inércia de rotação estão presentes no modelo, assim como os efeitos de fundação, força

magnética, força axial e mancais com propriedades variáveis. Simulações numéricas foram

realizadas, e o modelo se mostrou bem preciso. Tanto as comparações com outros modelos

numéricos como apresentado anteriormente, quanto comparações com medições experimentais,

não geraram diferenças relativas maiores do que 4.9% para os valores das primeiras duas

rotações críticas. A comparação com a referência [1] mostrou que os dois modelos estão

convergindo, porém o modelo deste trabalho exige um maior tempo computacional.

Propostas de controlo de vibração foram apresentadas. No que condiz ao equilíbrio óptimo

as simulações geraram resultados coerentes porém uma das curvas não mostrou um

comportamento semelhante ao experimental como comentado anteriormente. Porém três das

quatro curvas numéricas se mostraram bem coerentes em relação às curvas experimentais. Em

relação ao método de optimização, não foi possível observar uma repetibilidade de resultados,

porém com várias simulações foram encontrados bons resultados, isso pode ser explicado pelo

tipo de algoritmo utilizado. A proposta de alteração da geometria das chumaceiras gerou bons

resultados, chegando a reduzir algumas amplitudes de vibração pela metade. Uma dificuldade

encontrada foram as condições de fronteira para algumas chumaceiras, ou seja, condições estas

nas quais não existe perfil de pressão para sustentação do veio, ou o escoamento se torna

turbulento. Resultados melhores poderiam ter sido obtidos com chumaceiras que possuem

realimentação de óleo, ou resfriamento por convecção através de água, pois assim os

coeficientes de amortecimento são maiores devido ao controlo de temperatura.

A proposta de absorsor de vibração infelizmente não foi possível obter resultados devido a

escassez de tempo, tanto para implementação numérica quanto para construção do dispositivo

para realização de experimentos. Porém mesmo sem resultados foi decidido apresentar nesta

tese um modelo para acrescentar o efeito de um absorsor em um sistema rotativo, assim como o

seu respectivo projecto.

77

7. TRABALHOS FUTUROS

Em relação à determinação do desequilíbrio óptimo, bons resultados podem ser gerados com

outras técnicas de optimização, tais como técnicas de optimização não lineares (método quasi

Newton). As posições dos planos de equilíbrio poderiam ser consideradas variáveis a serem

optimizadas, e ainda a consideração de mais planos de equilíbrio poderiam gerar bons

resultados.

Um modelo preciso para determinação dos coeficientes de rigidez e amortecimento das

chumaceiras, com o objectivo de optimizar a geometria das mesmas ou a viscosidade do fluido,

para obter a menor amplitude de resposta possível, pode produzir boas condições de serviço

para máquinas rotativas.

78

8. REFERÊNCIAS

[1] A. Kokot, Modelo de Elementos Finitos para Caracterização Dinâmica de Rotores

Aplicado a Motores Eléctricos, Projecto, (2009).

[2] A. Muszynska, Rotordynamics, Taylor & Francis, Nevada, 2005.

[3] S. -G. Tan, X. -X. Wang, A theoretical introduction to low speed balancing of

flexible rotors: unification and development of the modal balancing and influence

coefficient techniques, Journal of Sound and Vibration, 168, (1993), pp 385-394.

[4] F. Sève, M. A. Andrianoely, A. Berlioz, R. Dufour, M. Charreyron, Balancing of

machinery with a flexible variable-speed rotor, Journal of Sound and Vibration, 264,

(2003), pp 287-302.

[5] J. Chung, I. Jang, Dynamic response and stability analysis of an automatic ball

balancer for flexible rotor, Journal of Sound and Vibration, 259, (2003), pp 31-43.

[6] A. Stodola, Steam and Gas Turbines, New York, McGraw-Hill, (1927).

[7] R. Green, Gyroscopic effects of critical speeds of a flexible rotor, American Society

of Mechanical Engineers Journal of Applied Mechanics, 15, (1948), pp 369-376.

[8] D. J. Rodrigues, A. R. Champneys, M. I. Friswell, R. E. Wilson, Automatic two-

plane balancing for rigid rotors, International Journal of Non-Linear Mechanics,

Elsevier, (2008).

[9] B. Xu, L. Qu, R. Sun, The optimization technique-based balancing of flexible rotors

without test runs, Journal of Sound and Vibration, 238, (2000), pp 877-892.

79

[10] F. J. Doubrawa Filho, Controle de vibração flexional em sistemas girantes

utilizando neutralizadores dinâmicos viscoelásticos, Curitiba, UTFPR, Monografia,

(2008).

[11] E. R. Ceccon, Controle de vibração torcional usando neutralizadores dinâmicos

viscoelásticos, Curitiba, UTFPR, Monografia, (2008).

[12] D. Childs, Turbomachinery Rotordynamics: Phenomena, Modeling, and Analysis,

New York, John Wiley & Sons, Inc., (1993).

[13] R. L. Norton, Projeto de Máquinas: Uma Abordagem Integrada, Porto Alegre,

Bookman, (2004).

[14] F. M. N. Farias, R. A. Ribas, Caracterização dinâmica de rotores trabalhando a

elevadas rotações: comportamento de mancais hidrodinâmicos, Curitiba, UTFPR,

Monografia, (2005).

[15] WEG, Hazardous Areas Motors: Low, Medium and High Voltage, (2009).

[16] K. E. Ruhl, J. L. Booker, A Finite Element Model for Distributed Parameter

Turborotor Systems, ASME Journal of Engineering for Industry, 94, (1972), pp 126-

132.

[17] M. Lalanne, G. Ferraris, Rotordynamics Prediction in Engineering, New York,

John Wiley & Sons Ltd, (1990).

[18] L. -W. Chen, D. -M. Ku, Analysis of whirl speeds of rotor-bearings systems with internal damping by C0 finite elements, Finite Elements in Analysis and Design, 9, (1991) pp. 169-176.

[19] R. C. Hibbler, Resistência dos Materiais, Rio de Janeiro, LTC, (2000).

[20] T. J. R. Hughes, The Finite Element Method, Mineola, Dover Publications, (2000).

[21] E. A. Ribeiro, A. Kokot, H. L. V. Santos, M. A. Luersen, C. A. Bavastri, Modelo

de fundação em máquinas rotativas utilizando parâmetros equivalentes generalizados,

Salvador, CONEM, CON08-1817, (2008).

80

[22] C. A. Bavastri, E. M. S. Ferreira, E. M. O. Lopes, Modeling of dynamic rotors with

flexible bearing using viscoelastic materials, XI DINAME, Ouro Preto, (2005).

[23] A. F. Pinho, J. A. B. Montevechi, F. A. S. Marins, Análise de aplicação de projeto

de experimentos nos parâmetros dos algoritmos genéticos, S&G Revista Eletrônica,

Sistemas & Gestão, 2, (2007), pp 319-331.

[24] F. W. Ocvirk, H. H. Mabie, Mechanisms and Dynamics of Machines, John Wiley

& Sons, New York, (1975).

[25] B. J. Hamrock, S. R. Schmid, B. O. Jacobson, Fundamentals of Film Fluid

Lubrication, CRC, New York, (2004).

81

APÊNDICE A – DETERMINAÇÃO DA MATRIZ DE RIGIDEZ

Os termos 1K e 2K da equação 4.21 são integtrados de forma completa como mostrado no

Quadro A 1.

Quadro A 1 – Termos de rigidez relativos à flexão.

Os termos relativos ao corte devem ser sub - integrados com o objectivo de evitar o shear

locking, ou bloqueio de solução. Para proceder a integração gaussiana selectiva deve-se

modificar os limites de integração para adequar o problema ao modelo:

dL

dy

yLL

2

11

2/2/

, A/1

e aplicando a quadratura gaussiana com peso w igual a 1, conforme a Tabela A 1, tem-se:

82

23

11

23

11

23

1

23

121

0

1

1

Lf

Lf

Lfw

LfwLdfdyyf

L

. A/2

Então pode-se calcular os termos relativos ao corte (Quadro A 2).

Tabela A 1 – Posição dos pontos e pesos para a quadratura gaussiana no intervalo de -1 a 1.

N (nº de pontos) Localização dos pontos Pesos

1 0.0 2.0

2 1

3 1.0

83

Quadro A 2 – Integração dos termos relativos ao corte.

84

Posicionando todos os termos na matriz global de rigidez, pode-se chegar na matriz de

rigidez, mostrado no Quadro A 3 a parte correspondente a flexão e mostrado no Quadro A 4 a

parte correspondente ao corte.

Quadro A 3 – Matriz global de rigidez devido a flexão.

Quadro A 3 – Matriz global de rigidez devido ao corte.

85

APÊNDICE B – FORÇA AXIAL

Ao realizar a integração, pode mostrar que a matriz de rigidez devido aos esforços axiais

pode ser obtida como mostrada no Quadro B 1.

Quadro B 1 – Matriz de rigidez devido aos esforços axiais.

Para o elemento de Timoshenko utilizando funções de interpolação quadráticas a matriz de

rigidez devido às forças magnéticas é mostrada no Quadro B 2.

Quadro B 2 – Matriz de rigidez devido às forças magnéticas.

86

APÊNDICE C – DETERMINAÇÃO DA MATRIZ DE MASSA

O Quadro C 1 apresenta o desenvolvimento das integrais das parcelas da equação 4.35.

Quadro C 1 – Parcelas correspondentes a matriz de massa elementar.

Baseado na ordenação dos graus de liberdade apresentada na equação 4.36, expandindo os

termos e aplicando as equações de Lagrange obtém-se a matriz elementar de massa (Quadro C

2).

87

Quadro C 2 – Matriz elementar de massa

88

APÊNDICE D – DETERMINAÇÃO DA MATRIZ GIROSCÓPICA

A partir do resultado obtido com o cálculo da integral da parcela correspondente aos efeitos

giroscópicos (equação 4.37) pode-se apresentar o Quadro D 1.

Quadro D 1 – Integração do termo correspondente ao efeito giroscópico.

Executando o procedimento de expansão dos termos e aplicando as equações de Lagrange

obtém-se a matriz giroscópica elementar (Quadro D 2).

Quadro D 2 – Matriz giroscópica.