SIMULAÇÃO NUMÉRICA DE UM ROTOR EM …monografias.poli.ufrj.br/monografias/monopoli10018281.pdf · these machines to evaluate their behavior, taking into account the presence of

SIMULAÇÃO NUMÉRICA DE UM ROTOR EM BANCADA DE TESTES

Rodrigo Martins de Oliveira

Projeto de Graduação apresentado ao Curso de

Engenharia Mecânica da Escola Politécnica,

Universidade Federal do Rio de Janeiro, como

parte dos requisitos necessários à obtenção do

título de Engenheiro Mecânico.

Orientador: Thiago Gamboa Ritto

Rio de Janeiro

Setembro de 2016

UNIVERSIDADE FEDERAL DO RIO DE JANEIRO

Departamento de Engenharia Mecânica

DEM/POLI/UFRJ



PROJETO FINAL SUBMETIDO AO CORPO DOCENTE DO DEPARTAMENTO DE

ENGENHARIA MECÂNICA DA ESCOLA POLITÉCNICA DA UNIVERSIDADE

FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS

NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE ENGENHEIRO MECÂNICO.

Aprovado por:

________________________________________________

Prof.º Thiago Gamboa Ritto, D.Sc.

________________________________________________

Prof.º Fernando Augusto de Noronha Castro Pinto, Dr.Ing.

________________________________________________

Prof.º Daniel Alves Castello, D.Sc.

_____________________________________________________

David Julian Gonzalez Maldonado, M.Sc.

RIO DE JANEIRO, RJ - BRASIL

SETEMBRO DE 2016

Oliveira, Rodrigo Martins de

Simulação numérica de um rotor em bancada de testes/ Rodrigo

Martins de Oliveira – Rio de Janeiro: UFRJ/ Escola Politécnica, 2016.

XII, 54 p.: il.; 29,7cm.


Projeto de Graduação – UFRJ/ Escola Politécnica/ Curso de

Engenharia Mecânica, 2016.

Referências Bibliográficas: p. 55.

1.Rotor. 2.Simulação numérica. 3.Método de elementos finitos.

4.Selo. 5.LAVI. I. Thiago Gamboa Ritto. II. Universidade Federal do Rio

de Janeiro, Escola Politécnica, Curso de Engenharia Mecânica. III.

Simulação numérica de um rotor em bancada de testes

iv

Resumo do Projeto de Graduação apresentado à Escola Politécnica / UFRJ como parte

dos requisitos necessários para a obtenção do grau de Engenheiro Mecânico.



Setembro/2016


Curso: Engenharia Mecânica

A aplicação de máquinas rotativas é bastante ampla no campo da engenharia. Desde

pequenos sistemas residenciais de bombeamento de água até grandes turbinas de aviões,

utilização em sistemas hidroelétricos ou na área de óleo e gás. Nesse contexto, é

necessário realizar estudo dinâmicos dessas máquinas para avaliar seu comportamento,

levando em consideração a presença de outros componentes envolvidos, tais como:

mancais, selos, discos, eixos, entre outros, para prever, evitar ou identificar falhas nos

equipamentos. Esse trabalho tem o objetivo principal de simular uma bancada de testes

desenvolvida pelo Laboratório de Acústica e Vibrações (LAVI) da COPPE em conjunto

da Petrobrás e verificar a influência dos parâmetros de rigidez e amortecimento de selos

mecânicos em um modelo de rotor. Para isso, é utilizado o programa LaviROT, baseado

em Matlab e em desenvolvimento no LAVI e o software ANSYS a fim de verificar os

resultados obtidos pelo primeiro.

Palavras-chave: rotor, simulação numérica, método de elementos finitos, selo, Lavi

v

Abstract of Undergraduate Project presented to POLI/UFRJ as a partial fulfillment of

the requirements for the degree of Engineer.

NUMERICAL SIMULATION OF A ROTOR IN TEST BENCH


September/2016

Advisor: Thiago Gamboa Ritto

Course: Mechanical Engineering

The application of rotating machines is very wide in the field of engineering. From small

residential water pumping systems to large aircraft turbines, hydroelectric systems in use

or in the oil and gas area. In this context, it is necessary to perform dynamic studies of

these machines to evaluate their behavior, taking into account the presence of other

components involved, such as bearings, seals, discs, shafts, among others, to anticipate,

prevent or identify equipment failures. This work aims, mainly, to simulate a test bench

developed by the Acoustics and Vibrations Laboratory (LAVI) COPPE together

Petrobras in order to verify the influence of the stiffness and damping parameters of

mechanical seals in a rotor model. For this, it is used the LaviROT, a program based on

Matlab under development in LAVI, and ANSYS software to verify the results of the

first.

Keywords: rotor, numerical simulation, finite-element analysis, seal, Lavi

vi

AGRADECIMENTOS

Agradeço a minha família por toda educação que me deram, carinho, apoio, e

sacrifícios que fizeram por mim. Todo meu crescimento pessoal e profissional e todas as

vitórias ao longo de meu caminho da vida devo a vocês.

À Daniela, que tem sido uma companheira maravilhosa na minha vida, me

apoiando, aconselhando e incentivando em todas as minhas escolhas. Obrigado por toda

a paciência, compreensão e carinho.

À Thaíssa, pela nossa grande amizade, e que mesmo distante, sempre me deu

apoio para seguir em frente em meu caminho.

Ao Professor Thiago Ritto por toda a passagem de conhecimento nos projetos que

realizei sob sua orientação.

À Equipe Minerva Baja UFRJ, onde passei três anos de minha formação e ganhei

uma segunda família, em especial Honda, Felipe Cristaldi, Adriano, Fred Fróes, Douglas,

Larissa, Fernando, César, Fredinho, Pedro Fusca, Casagrande e Vinicius. Devo a equipe

boa parte de meu aprendizado em engenharia.

À toda equipe do LAVI, que esteve comigo neste último ano, em especial David,

quem me ajudou bastante na realização deste projeto.

À toda equipe do Optimize, onde trabalhei por um ano e contribuiu para minha

formação, em especial Zélia e Elmer.

Ao médico oncologista Fladwmyr Barros e às fonoaudiólogas Nathalia Castellani

e Roberta Leite, que foram essenciais na minha recuperação, serei para sempre grato por

tudo.

E a todos os amigos e professores que fizeram parte da minha vida e contribuíram

para que eu alcançasse este momento. Obrigado a todos.

vii

SUMÁRIO

LISTA DE FIGURAS ..................................................................................................... ix

LISTA DE TABELAS ................................................................................................... xii

1 INTRODUÇÃO ............................................................................................................ 1

1.1 Contextualização .................................................................................................. 1

1.2 Objetivo ................................................................................................................ 1

1.3 Motivação ............................................................................................................. 1

1.4 Metodologia .......................................................................................................... 2

2 BASE TEÓRICA E ROTODINÂMICA ...................................................................... 3

2.1 Rotor de Jeffcott ................................................................................................... 3

2.2 Rotor com efeito giroscópico ............................................................................... 5

3 MODELAGEM POR ELEMENTOS FINITOS .......................................................... 9

3.1 Introdução ............................................................................................................. 9

3.2 Aplicação ao LaviROT ......................................................................................... 9

4 LaviROT ..................................................................................................................... 13

4.1 Introdução ........................................................................................................... 13

4.2 LaviROT versão 2015 ........................................................................................ 13

4.3 LaviROT versão 2016 ........................................................................................ 14

4.4 Novas opções de análises ................................................................................... 14

4.4.1 Análise de estabilidade ................................................................................. 14

4.4.2 Diagrama de Nyquist .................................................................................... 15

4.4.3 Análise do regime transiente e permanente do domínio do tempo .............. 15

5 COMPARAÇÃO ANSYS VS LaviROT ................................................................... 18

5.1 Modelo de teste inicial ........................................................................................ 18

5.2 Cálculo da rigidez do rolamento ......................................................................... 18

5.3 Simulação em ANSYS Workbench .................................................................... 22

5.3.1 Introdução ..................................................................................................... 22

viii

5.3.2 Geometria ..................................................................................................... 22

5.3.3 Geração de malha ......................................................................................... 23

5.3.4 Análise rotodinâmica .................................................................................... 24

5.3.5 Análise Modal .............................................................................................. 26

5.3.6 Análise estática ............................................................................................. 27

5.4 Simulação em Matlab – LaviROT ...................................................................... 28

5.4.1 Introdução ..................................................................................................... 28

5.4.2 Condições de contorno ................................................................................. 29

5.4.3 Resultados ..................................................................................................... 30

5.5 Comparação dos resultados ................................................................................ 32

6 SIMULAÇÃO DA BANCADA DE TESTES ........................................................... 33

6.1 Introdução ........................................................................................................... 33

6.2 Análise do rotor sem selos .................................................................................. 34

6.3 Análise do rotor com selos e sem força de pulso ............................................... 38

6.4 Análise do rotor com selos e com força de pulso ............................................... 42

6.5 Análise do rotor com selos de diferentes rigidezes diretas ................................. 44

6.6 Análise do rotor com selos de diferentes rigidezes acopladas............................ 46

6.7 Análise do rotor com selos de diferentes coeficientes de amortecimento diretos

51

7 CONCLUSÕES E TRABALHOS FUTUROS .......................................................... 54

8 REFERÊNCIAS BIBLIOGRÁFICAS ....................................................................... 55

APÊNDICE A – Desenho detalhado do eixo ................................................................. 56

ix

LISTA DE FIGURAS

Figura 1.1: Corte de um selo labirinto. ............................................................................. 1

Figura 2.1: A ilustração a direita mostra o movimento de precessão direta e a ilustração

da esquerda mostra o movimento de precessão retrógada. ............................................... 3

Figura 2.2: Modelo em 3D do rotor de Jeffcott [2]. ......................................................... 4

Figura 2.3: Vista lateral do modelo do rotor de Jeffcott [1]. ............................................ 4

Figura 2.4: Representação lateral de um rotor com o referencial inercial A e o referencial

móvel B. ........................................................................................................................... 6

Figura 2.5: Representação frontal de um rotor com os referenciais C e D. ...................... 6

Figura 3.1: Ilustração do elemento com 4 graus de liberdade. ....................................... 10

Figura 4.1: Interface do LaviROT 2015. ........................................................................ 13

Figura 4.2: Interface do LaviROT 2016. ........................................................................ 14

Figura 4.3: Detalhe do quadro de configuração da análise transiente. ........................... 17

Figura 5.1: Pontos de contato da esfera com as pistas interna e externa (adaptado de [10]).

........................................................................................................................................ 19

Figura 5.2: Representação do ângulo de contato (adaptado de [10]). ............................ 20

Figura 5.3: Geometria do eixo gerada no ANSYS. ........................................................ 23

Figura 5.4: Geometria do eixo com as seções visíveis. .................................................. 23

Figura 5.5: Detalhe das condições de contorno da análise rotodinâmica no ANSYS. ... 24

Figura 5.6: Detalhe dos dados do rolamento da direita. ................................................. 25

Figura 5.7: Detalhe dos dados do rolamento da esquerda. ............................................. 25

Figura 5.8: Diagrama de Campbell gerado pelo ANSYS. .............................................. 25

Figura 5.9: Gráfico Amplitude (mm) x Velocidade de rotação [RPM].......................... 26

Figura 5.10: Frequências naturais do sistema geradas no ANSYS. ............................... 26

Figura 5.11: Velocidades de rotação críticas e direção de precessão geradas no ANSYS.

........................................................................................................................................ 26

Figura 5.12: Frequências naturais do eixo livre geradas no ANSYS. ............................ 27

Figura 5.13: Detalhe do posicionamento das forças e dos seus respectivos módulos. ... 27

Figura 5.14: Resultado da análise estática gerada no ANSYS. ...................................... 28

Figura 5.15: Ilustração da separação das regiões do eixo. .............................................. 28

Figura 5.16: Interface do LaviROT com os dados para a simulação de teste. ............... 29

Figura 5.17: Diagrama de Campbell da simulação teste gerado no LaviROT. .............. 30

x

Figura 5.18: Gráfico da amplitude vs frequência para a simulação teste gerado no

LaviROT. ........................................................................................................................ 30

Figura 5.19: Resposta transiente da simulação teste gerado no LaviROT. .................... 31

Figura 5.20: Órbita da resposta transiente da simulação teste gerado no LaviROT. ..... 31

Figura 6.1: Diagrama de Campbell do primeiro modelo. ............................................... 34

Figura 6.2: Gráfico da amplitude vs frequência para disco 1 do primeiro modelo. ....... 35

Figura 6.3: Gráfico da amplitude vs frequência para disco 2 do primeiro modelo. ....... 35

Figura 6.4: Análise transiente do disco 1 do primeiro modelo. ...................................... 36

Figura 6.5: Órbita do disco 1 do primeiro modelo. ........................................................ 36

Figura 6.6: Análise transiente do disco 2 do primeiro modelo. ...................................... 37

Figura 6.7: Órbita do disco 2 do primeiro modelo. ........................................................ 37

Figura 6.8: Diagrama de Campbell do segundo modelo. ............................................... 38

Figura 6.9: Gráfico da amplitude em função da frequência para o disco 1 do segundo

modelo. ........................................................................................................................... 39

Figura 6.10: Gráfico da amplitude em função da frequência para o disco 2 do segundo

modelo. ........................................................................................................................... 39

Figura 6.11: Análise transiente para o disco 1 do segundo modelo. .............................. 40

Figura 6.12: Órbita do disco 1 do segundo modelo comparada com o do primeiro. ..... 40

Figura 6.13: Análise transiente para o disco 2 do segundo modelo. .............................. 41

Figura 6.14: Órbita do disco 2 do segundo modelo comparada com o do primeiro. ..... 41

Figura 6.15: Análise transiente para o disco 1 do terceiro modelo. ............................... 43

Figura 6.16: Análise transiente para o disco 2 do terceiro modelo. ............................... 43

Figura 6.17: Variação da amplitude máxima (em módulo) no regime transiente em função

da rigidez direta do selo. ................................................................................................. 44

Figura 6.18: Variação da amplitude máxima (em módulo) no regime permanente em

função da rigidez direta do selo. ..................................................................................... 45

Figura 6.19: Tempo de duração do regime transiente em função da rigidez direta do selo.

........................................................................................................................................ 45

Figura 6.20: Variação das frequências críticas de rotação em função da rigidez direta do

selo. ................................................................................................................................. 46

Figura 6.21: Variação da amplitude máxima no regime transiente em função da rigidez

acoplada do selo. ............................................................................................................ 47

Figura 6.22: Variação da amplitude máxima no regime permanente em função da rigidez

acoplada do selo. ............................................................................................................ 47

xi

Figura 6.23: Tempo de duração do regime transiente em função da rigidez acoplada do

selo. ................................................................................................................................. 48

Figura 6.24: Variação das frequências críticas de rotação em função da rigidez acoplada

do selo. ............................................................................................................................ 48

Figura 6.25: Gráfico amplitude em função da rigidez acoplada do selo para apontar

instabilidade. ................................................................................................................... 49

Figura 6.26: Diagramas de Nyquist para diferentes valores de rigidezes acopladas do selo.

........................................................................................................................................ 50

Figura 6.27: Variação da amplitude máxima no regime transiente em função do

amortecimento direto do selo. ........................................................................................ 52

Figura 6.28: Variação da amplitude máxima no regime permanente em função do

amortecimento direto do selo. ........................................................................................ 52

Figura 6.29: Tempo de duração do regime transiente em função do amortecimento direto

do selo. ............................................................................................................................ 52

xii

LISTA DE TABELAS

Tabela 5.1: Propriedades aplicadas ao eixo na simulação no LaviROT. ....................... 28

Tabela 6.1: Dados dos discos. ........................................................................................ 33

Tabela 6.2: Dados dos selos. .......................................................................................... 33

Tabela 6.3: Velocidades críticas e direções de precessão do primeiro modelo. ............. 38

Tabela 6.4: Velocidades críticas e direções de precessão do segundo modelo. ............. 42

1

1 INTRODUÇÃO

1.1 Contextualização

Máquinas rotativas têm uma utilização bastante ampla de maneira geral, seja em

aplicação de aparelhos domésticos como ventiladores e lavadoras, seja na aplicação da

indústria como turbinas e bombas.

A criação de novas máquinas acompanha o desenvolvimento de pesquisas nessa

área para o aumento da eficiência das mesmas, uma vez que envolvem geração e consumo

de energia, e por consequência um custo financeiro significativo.

Por esse motivo, a simulação de máquinas rotativas (rotores) é essencial para

avaliar o seu comportamento e verificar como cada parâmetro atua nesse aspecto.

Um desses parâmetros é o selo interno ou anular (Figura 1.1). Os selos anulares

têm sido objeto de estudo com a finalidade de caracterizá-los a partir de coeficientes

dinâmicos de rigidez e amortecimento. E assim, analisar como estes coeficientes atuam

na dinâmica do rotor.

Figura 1.1: Corte de um selo labirinto.

1.2 Objetivo

Esse trabalho tem o objetivo principal de simular uma bancada de testes

desenvolvida pelo Laboratório de Acústica e Vibrações (LAVI) da COPPE em conjunto

da Petrobrás e verificar a influência dos parâmetros de rigidez e amortecimento de selos

anulares em um modelo de rotor, e a verificação e ampliação do LaviROT, um software

baseado em Matlab, em desenvolvimento no LAVI para a simulação de rotores.

1.3 Motivação

A principal motivação deste trabalho é um projeto de pesquisa, em andamento,

com a Petrobrás para a caracterização de selos anulares do tipo labirinto, no caráter de

2

determinar coeficientes de amortecimento e rigidez experimentalmente através de uma

bancada de testes, para futuramente prever suas influências no comportamento do

sistema.

1.4 Metodologia

Para atingir esses objetivos, o trabalho foi dividido em três partes distintas.

A primeira foi implementar novos algoritmos ao LaviROT, ou seja, aumentar a

quantidade de resultados relevantes que o programa pode gerar, e também ampliar a

variedade de dados de entrada na simulação, podendo assim, criar novas situações reais

para simular.

A segunda fase foi comparar os resultados iniciais com um software comercial

difundido. Para isso, criou-se um modelo base de um eixo com um disco no centro

apoiado em mancais de rolamento em suas extremidades. Esse modelo foi simulado no

LaviROT e os resultados foram verificados com a utilização de outro software de

simulação por elementos finitos, o ANSYS.

A terceira fase foi a preparação e simulação da geometria do rotor da bancada de

testes com os mancais e selos posicionados utilizando o programa LaviROT.

3

2 BASE TEÓRICA E ROTODINÂMICA

Rotodinâmica é um termo usado para denominar a dinâmica de rotores. Para um

modelo simples de rotor, pode-se considerar um eixo rotativo, podendo ser rígido ou

flexível, com um disco acoplado e mancais de apoio, podendo ser rígidos ou flexíveis

também. O movimento descrito pelo centro geométrico do eixo segue a trajetória de uma

órbita (circular ou elíptica). Esta órbita pode seguir o mesmo sentido de giro do eixo

(precessão direta) ou pode seguir o sentido contrário desse giro (precessão retrógrada).

Esse comportamento é ilustrado na Figura 2.1.

Figura 2.1: A ilustração a direita mostra o movimento de precessão direta e a ilustração da

esquerda mostra o movimento de precessão retrógada.

2.1 Rotor de Jeffcott

O modelo de Rotor de Jeffcott é um modelo proposto por August Föppl (1895) e

por Henry Homan Jeffcott (1919) [1] para simplificar problemas complexos de

rotodinâmica. Basicamente, ele consiste em um disco rígido desbalanceado fixado em um

eixo flexível, uniforme e de massa desprezível, apoiado em mancais.

Neste modelo, podem ser feitas as seguintes considerações: o rotor gira com

velocidade constante 𝜔 junto com o disco desbalanceado de massa 𝑚 como mostra a

Figura 2.2. O centro geométrico 𝐶 do disco, dado pelas coordenadas 𝑈𝑥𝐶 , 𝑈𝑦𝐶 não

coincide com seu centro de massa 𝐺, dado pelas coordenadas 𝑈𝑥𝐺 , 𝑈𝑦𝐺 devido ao

desbalanceamento, mas ainda assim pertencem ao mesmo plano, como mostra a Figura

2.3.

4

Figura 2.2: Modelo em 3D do rotor de Jeffcott [2].

Figura 2.3: Vista lateral do modelo do rotor de Jeffcott [1].

Assumindo que a rigidez do disco não afeta a rigidez do eixo, pode-se determinar

a rigidez de deflexão no centro de uma viga uniforme bi apoiada pela seguinte equação

[3].

𝑘𝑠 =48 ∙ 𝐸 ∙ 𝐼

𝐿3 (2.1)

5

Onde, E é o módulo de elasticidade do material eixo, L o comprimento total do

eixo e I o momento de inércia de área do eixo, que para o caso de um eixo cilíndrico é

dado por:

𝐼 =𝜋 ∙ 𝐷4

64 (2.2)

Para desenvolver analiticamente o modelo, é necessário aplicar as leis do

movimento de Newton ao disco. Uma vez que a massa do eixo é considerada desprezível,

todas as forças atuantes no disco serão devido à inércia, à rigidez e ao amortecimento de

deflexão do eixo. As equações de movimento nas direções x e y são dadas por:

𝑚𝑥𝐺 = −𝑘𝑠𝑢𝑥𝐶 − 𝑐𝑠𝑥𝐶 (2.3)

𝑚𝑦𝐺 = −𝑘𝑠𝑢𝑦𝐶 − 𝑐𝑠𝑦𝐶 (2.4)

Também é possível reescrever as coordenadas do centro de massa do disco a partir

do centro geométrico e da posição angular do eixo 𝜔𝑡 no tempo t, gerando as seguintes

equações:

𝑢𝑥𝐺 = 𝑢𝑥𝐶 + 𝑒 ∙ cos (𝜔𝑡) (2.5)

𝑢𝑦𝐺 = 𝑢𝑦𝐶 + 𝑒 ∙ sen(𝜔𝑡) (2.6)

Substituindo as coordenadas do centro de massa do disco obtidas nas equações

2.5 e 2.6 nas equações 2.3 e 2.4, obtém-se as equações de movimento do rotor de Jeffcott

em relação ao seu centro geométrico [4]:

𝑚𝑥𝐶 + 𝑘𝑠𝑢𝑥𝐶 + 𝑐𝑠𝑥𝐶 = 𝑚 ∙ 𝑒 ∙ 𝜔2 ∙ cos(𝜔𝑡) (2.7)

𝑚𝑦𝐶 + 𝑘𝑠𝑢𝑦𝐶 + 𝑐𝑠𝑦𝐶 = 𝑚 ∙ 𝑒 ∙ 𝜔2 ∙ sen(𝜔𝑡) (2.8)

2.2 Rotor com efeito giroscópico

Nas Figuras 2.4 e 2,5, pode-se observar os referencias adotados, sendo A (base

𝑎1, 𝑎2, 𝑎3 ) o referencial inercial e B (𝑏1, 𝑏2, 𝑏3), que admite um giro em torno de 𝑎2.

Além dessas bases, criou-se uma base C (𝑐1, 𝑐2, 𝑐3) e uma base D (𝑑1, 𝑑2, 𝑑3). A base

C realiza um giro 𝛽 em torno de 𝑏3 e a base D realiza um giro 𝜃 em torno de 𝑐1, que seria

a rotação do disco.

𝝎𝐷𝐴 = 𝜔𝐵𝐴 + 𝜔𝐶𝐵 + 𝜔𝐷𝐶 (2.9)

6

Figura 2.4: Representação lateral de um rotor com o referencial inercial A e o referencial móvel B.

Figura 2.5: Representação frontal de um rotor com os referenciais C e D.

Para manter o modelo simplificado, algumas hipóteses devem ser tomadas, são

elas: velocidade de rotação () do disco constante, 𝛼 e 𝛽 são considerados pequenos, logo

seus respectivos cossenos são aproximadamente 1 e seus respectivos senos são

aproximadamente 𝛼 e 𝛽. Os termos multiplicativos entre eles serão aproximadamente 0.

Considerando essas bases, pode-se escrever a velocidade angular 𝜔 da seguinte

maneira:

𝜔𝐷𝐴 = 𝒃𝟐 + 𝒄𝟑 + 𝒄𝟏 (2.11)

Escrevendo o vetor velocidade angular na base C, tem-se:

7

𝝎𝐷𝐵𝐴 = [ 𝑇𝑎𝑏 ] 𝝎𝐷

𝐴𝐴 (2.11)

𝝎𝐷𝐶𝐴 = [ 𝑇𝑏𝑐 ] 𝝎𝐷

𝐵𝐴 (2.12)

Onde:

[ 𝑇𝑎𝑏 ] = [𝑐𝑜𝑠𝛼 0 𝑠𝑒𝑛𝛼

0 1 0−𝑠𝑒𝑛𝛼 0 𝑐𝑜𝑠𝛼

] (2.13)

[ 𝑇𝑐𝑏 ] = [𝑐𝑜𝑠𝛽 −𝑠𝑒𝑛𝛽 0𝑠𝑒𝑛𝛽 𝑐𝑜𝑠𝛽 0

0 0 1

] (2.14)

Logo, a velocidade angular pode ser escrita da seguinte maneira:

𝝎 = (𝑠𝑒𝑛𝛽 + )𝒄𝟏 + 𝑐𝑜𝑠𝛽𝒄𝟐 + 𝒄𝟑 (2.15)

A derivada de 𝝎 no referencial 𝑪 é dada por:

= (𝑠𝑒𝑛𝛽 + 𝑐𝑜𝑠𝛽)𝒄𝟏 + (𝑐𝑜𝑠𝛽 − 𝑐𝑜𝑠𝛽)𝒄𝟐 + 𝒄𝟑 (2.16)

Para obter os momentos que atuam no disco, utiliza-se as equações de Euler do

movimento. A velocidade angular Ω do sistema é diferente da velocidade angular 𝜔 do

corpo axissimétrico. As equações são dadas por:

𝜔𝐷𝐴 = Ω = (𝑠𝑒𝑛𝛽)𝒄𝟏 + (𝑐𝑜𝑠𝛽)𝒄𝟐 + 𝒄𝟑 (2.17)

𝑀𝑥 = 𝐼𝑥𝑥𝑥 + 𝐼𝑧𝑧𝜔𝑧Ω𝑦 + 𝐼𝑦𝑦𝜔𝑦Ω𝑧 (2.18)

𝑀𝑦 = 𝐼𝑦𝑦𝑦 + 𝐼𝑥𝑥𝜔𝑥Ω𝑧 + 𝐼𝑧𝑧𝜔𝑧Ω𝑥 (2.19)

𝑀𝑧 = 𝐼𝑧𝑧𝑧 + 𝐼𝑦𝑦𝜔𝑦Ω𝑧 + 𝐼𝑥𝑥𝜔𝑥Ω𝑦 (2.20)

Aplicando as equações 2.15 da velocidade e 2.16 da aceleração às equações 2.18,

2.19 e 2.20, tem-se:

𝑀𝑥 =1

2𝑚𝑟²(𝑠𝑒𝑛𝛽 + 𝑐𝑜𝑠𝛽) (2.21)

𝑀𝑦 =1

4𝑚𝑟2(𝑐𝑜𝑠𝛽 − 𝑐𝑜𝑠𝛽) +

1

2𝑚𝑟2(𝑠𝑒𝑛𝛽 + ) −

1

4𝑚𝑟² (2.22)

𝑀𝑧 =1

4𝑚𝑟2 +

1

4𝑚𝑟22𝑠𝑒𝑛𝛽𝑐𝑜𝑠𝛽 −

1

2𝑚𝑟2(𝑠𝑒𝑛𝛽 + )𝑐𝑜𝑠𝛽 (2.23)

O momento resistivo é proporcional aos ângulos 𝛼 e 𝛽 pela constante de rigidez

torcional 𝐾𝑡 e proporcional às velocidades angulares e pela constante de

amortecimento C. Sendo assim, a equação do momento aplicado é:

𝑀 = −(𝐶 + 𝐾𝑡𝛼)𝒃𝟐− (𝐶 + 𝐾𝑡𝛽)𝒄

𝟑 (2.24)

Fazendo 𝑀𝑐 = [ 𝑇𝑏𝑐 ]𝑀𝑏, tem-se:

𝑀 = −(𝐶𝑠𝑒𝑛𝛽 + 𝐾𝑡𝛼𝑠𝑒𝑛𝛽)𝒄𝟏 − (𝐶𝑐𝑜𝑠𝛽 + 𝐾𝑡𝛼𝑐𝑜𝑠𝛽)𝒄𝟐 − 𝐾𝑡𝛽𝒄𝟑 (2.25)

8

Igualando a equação 2.25 às equações 2.21, 2.22 e 2.23, obtém-se as equações

diferenciais que representam o movimento do disco.

𝐼 = 𝐶 + 𝐾𝛼 + 2𝐼𝑝𝜔 (2.26)

𝐼 = 𝐶 + 𝐾𝛽 − 2𝐼𝑝𝜔 (2.27)

Onde:

𝐼 =1

2𝑚𝑟2 (2.28)

𝐼𝑝 =1

4𝑚𝑟4 (2.29)

9

3 MODELAGEM POR ELEMENTOS FINITOS

3.1 Introdução

A análise por elementos finitos é um método para aproximar soluções de equações

diferencias, amplamente utilizado em softwares de simulação numérica, que consiste em

discretizar o sistema sob análise de vários elementos. O método subdivide um problema

maior em partes menores e mais simples, chamadas de elementos finitos. As equações

simples que modelam estes elementos são então montadas em um sistema maior de

equações que descrevem todo o problema. Elementos finitos usa métodos de cálculo

variacional para aproximar uma solução, minimizando uma função de erro associada.

3.2 Aplicação ao LaviROT

O LaviROT é um programa que também faz uso desse método como estratégia

para resolver o problema rotor-mancal seguindo o modelo de viga de Euler-Bernoulli [5].

Neste modelo, a área da seção transversal do eixo permanece ortogonal à linha neutra do

eixo. A equação deste modelo para uma viga sem amortecimento, em um plano de duas

dimensões, é escrita como [4]:

ρA∂²𝑣(x, t)

∂t2+ EI

∂4v(x, t)

∂x4= f(x, t) (3.1)

onde v é o deslocamento transversal do eixo, E o módulo de elasticidade do material, 𝜌 é

a massa específica do material do eixo, A é a área da seção transversal, I é o momento de

inércia desta área, f é uma força externa por unidade de comprimento, x e t são parâmetros

de posição e tempo definidos pelos intervalos [0,L] e [0,T], onde L é o comprimento do

eixo.

De maneira análoga, é possível calcular o deslocamento 𝑤(𝑥, 𝑡) na direção

ortogonal a 𝑣(𝑥, 𝑡). Devido a rotação do eixo, um acoplamento giroscópico existe entre

os deslocamentos v e w [4]. As equações resultantes são discretizadas pelo método de

elementos finitos [4, 6]. Cada elemento tem 2 nós com 4 graus de liberdade cada um, duas

translações e duas rotações, como mostra a Figura 3.1.

10

Figura 3.1: Ilustração do elemento com 4 graus de liberdade.

O deslocamento do elemento é dado por:

𝑢(𝑒) = [𝑁] 𝑢(𝑒) (3.2)

onde:

𝑢(𝑒) = (𝑣1𝜃𝑧1 𝑤1𝜃𝑦1 𝑣2𝜃𝑧2 𝑤2𝜃𝑦2)𝑇 (3.3)

são deslocamentos generalizados dos nós e [N] é composto pelas funções de forma, que

são as funções cúbicas Hermitianas neste caso. Supõe-se que o equipamento mantém uma

velocidade constante de rotação Ω, em torno do eixo X.

A coordenada do elemento é escrita como 𝜉 = 𝑥/𝑙𝑒, onde 𝑙𝑒 é o comprimento do

elemento. Supondo pequenos deslocamentos para os graus de liberdade de translação e

rotação. As matrizes elementares de massa, giroscópica e de rigidez do eixo são dadas

por:

[𝑀𝑠](𝑒) = 𝜌𝐴∫(𝑁𝑣

𝑇𝑁𝑣

1

0

+ 𝑁𝑤𝑇𝑁𝑤) ∙ 𝑙𝑒𝑑𝜉 (3.4)

[𝐺𝑠](𝑒) =

𝜌𝑙𝑝𝑙𝑒2

Ω∫(𝑁𝑣′𝑇𝑁′𝑣 −

1

0

𝑁𝑤′𝑇𝑁′𝑤) ∙ 𝑙𝑒𝑑𝜉 (3.5)

[𝐾𝑠](𝑒) =

𝐸𝐼

𝑙𝑒4

∫(𝑁𝑣′′𝑇𝑁′′𝑣

1

0

+ 𝑁𝑤′′𝑇𝑁′′𝑤) ∙ 𝑙𝑒𝑑𝜉 (3.6)

onde I é o momento de inércia de área e 𝑙𝑝 é o momento polar de inércia. As funções de

interpolação são definidas para um elemento padrão com coordenada 𝜉 ∈ [0,1].

𝑁𝑣 = [𝑁𝑣1 𝑁𝑣2 0 0 𝑁𝑣3 𝑁𝑣4 0 0] (3.7)

𝑁𝑣′ = [𝑁𝑣1

′ 𝑁𝑣2′ 0 0 𝑁𝑣3

′ 𝑁𝑣4′ 0 0] (3.8)

11

𝑁𝑤 = [0 0 𝑁𝑤1 − 𝑁𝑤2 0 0 𝑁𝑤3 − 𝑁𝑤4] (3.9)

𝑁𝑤′ = [0 0 − 𝑁𝑤1

′ 𝑁𝑤2′ 0 0 − 𝑁𝑤3

′ 𝑁𝑤4′ ] (3.10)

e as funções individuais são dadas por:

𝑁𝑣1 = 1 − 3𝜉2 + 2𝜉³ (3.11)

𝑁𝑣2 = 𝑙𝑒(𝜉 − 2𝜉2 + 𝜉3) (3.12)

𝑁𝑣3 = 3𝜉2 − 2𝜉³ (3.13)

𝑁𝑣4 = 𝑙𝑒(−𝜉2 + 𝜉3) (3.14)

Após a integração, chega-se às matrizes elementares. A matriz de massa do

elemento do eixo é dada por:

[𝑀𝑠](𝑒) = 𝜌𝐴𝑙𝑒

[ 𝑀01 𝑀03 0 0 𝑀04 𝑀06 0 0𝑀03 𝑀02 0 0 −𝑀06 𝑀05 0 00 0 𝑀01 −𝑀03 0 0 𝑀04 −𝑀06

0 0 −𝑀03 𝑀02 0 0 𝑀06 𝑀05

𝑀04 −𝑀06 0 0 𝑀01 −𝑀03 0 0𝑀06 𝑀05 0 0 −𝑀03 𝑀02 0 00 0 𝑀04 𝑀06 0 0 𝑀01 𝑀03

0 0 𝑀06 𝑀05 0 0 𝑀03 𝑀01 ]

(3.15)

onde 𝑀01 = 78/210, 𝑀02 = (8/840)𝑙𝑒2, 𝑀03 = (44/840)𝑙𝑒, 𝑀04 = 27/210, 𝑀05 =

(−6/840)𝑙𝑒2, 𝑀06 = (−26/840)𝑙𝑒. A matriz giroscópica é dada por:

[𝐺𝑠](𝑒) = 𝜌𝑙𝑒𝑙𝑝Ω

[

0 0 −𝐺01 𝐺02 0 0 𝐺01 𝐺02

0 0 −𝐺02 𝐺03 0 0 𝐺02 −𝐺04

𝐺01 𝐺02 0 0 −𝐺01 𝐺02 0 0−𝐺02 −𝐺03 0 0 𝐺02 𝐺04 0 0

0 0 𝐺01 −𝐺02 0 0 −𝐺01 −𝐺02

0 0 −𝐺02 −𝐺04 0 0 𝐺02 𝐺03

−𝐺01 −𝐺02 0 0 𝐺01 −𝐺02 0 0−𝐺02 𝐺04 0 0 𝐺02 −𝐺03 0 0 ]

(3.16)

onde 𝐺01 = 6/(5𝑙𝑒2), 𝐺02 = 1/(10𝑙𝑒), 𝐺03 = 4/30, 𝐺04 = 1/30. E a matriz de rigidez

do eixo é dada por:

[𝐾𝑠](𝑒) =

𝐸𝐼

𝑙𝑒

[ 𝐾01 𝐾03 0 0 𝐾04 𝐾03 0 0𝐾03 𝐾02 0 0 −𝐾03 𝐾05 0 00 0 𝐾01 −𝐾03 0 0 𝐾04 −𝐾03

0 0 −𝐾03 𝐾02 0 0 𝐾03 𝐾05

𝐾04 −𝐾03 0 0 𝐾01 −𝐾03 0 0𝐾03 𝐾05 0 0 −𝐾03 𝐾02 0 00 0 𝐾04 𝐾03 0 0 𝐾01 𝐾03

0 0 𝐾03 𝐾05 0 0 𝐾03 𝐾01 ]

(3.17)

onde 𝐾01 = 12/𝑙𝑒2, 𝐾02 = 4, 𝐾03 = 6/𝑙𝑒, 𝐾04 = −12/𝑙𝑒

2, 𝐾05 = 2.

Mancais, selos e discos são representados apenas em um nó, logo a matriz de cada

um deles tem a dimensão 4 x 4. As matrizes de rigidez e amortecimento dos mancais e

selos são dadas por:

12

[𝐾𝑏](𝑛) = [

𝑘𝑦𝑦 0 𝑘𝑦𝑧 0

0 0 0 0𝑘𝑧𝑦 0 𝑘𝑧𝑧 0

0 0 0 0

] (3.18)

[𝐶𝑏](𝑛) = [

𝑐𝑦𝑦 0 𝑐𝑦𝑧 0

0 0 0 0𝑐𝑧𝑦 0 𝑐𝑧𝑧 0

0 0 0 0

] (3.19)

onde 𝑘𝑦𝑦, 𝑘𝑦𝑧, 𝑘𝑧𝑦 e 𝑘𝑧𝑧 são coeficientes de rigidez de mancais ou selos e 𝑐𝑦𝑦, 𝑐𝑦𝑧, 𝑐𝑧𝑦 e

𝑐𝑧𝑧 são coeficientes de amortecimento. As matrizes de massa e efeito giroscópico dos

discos são dadas por:

[𝑀𝑑](𝑛) = [

𝑀𝑑 0 0 00 𝐼𝑑 0 00 0 𝑀𝑑 00 0 0 𝐼𝑑

] (3.20)

[𝐺𝑑](𝑛) = [

0 0 0 00 0 0 −2𝐼𝑑Ω0 0 0 00 2𝐼𝑑Ω 0 0

] (3.21)

onde 𝑀𝑑 é a massa do disco e 𝐼𝑑 é o momento de inércia do disco relacionado aos eixos

Y e Z.

No laviROT, as matrizes de massa, amortecimento e rigidez tem suas dimensões

definidas a partir do número de graus de liberdade que o problema de elementos finitos

apresenta, e este está diretamente relacionado ao número de nós.

O número de nós do modelo é definido pela equação abaixo:

𝑛𝑛𝑜 = 𝑛𝑒𝑙 + 1 (3.22)

onde 𝑛𝑛𝑜 é o número de nós e 𝑛𝑛𝑒𝑙 é o número de elementos

Por sua vez, o número de graus de liberdade total (𝑛𝐺𝐷𝐿) é definido pelo número

de nós vezes o número de graus de liberdade por nó (𝑛𝑛𝑜𝑒𝑙):

𝑛𝐺𝐷𝐿 = 𝑛𝑛𝑜 ∙ 𝑛𝑛𝑜𝑒𝑙 (3.23)

13

4 LaviROT

4.1 Introdução

A ferramenta LaviROT é um programa em desenvolvimento, na plataforma do

Matlab, baseado em elementos finitos para a simulação de rotores. Seu objetivo é ser um

programa capaz de realizar cálculos e diagramas para avaliação de modelos de rotores,

oferecendo ao LAVI uma ferramenta para desenvolvimento de programa na área de

rotodinâmica. A ferramenta permite a simulação de eixos com mancais, discos e selos

acoplados e suas respectivas rigidezes e amortecimentos. São levados em conta a inércia,

rigidez, amortecimento viscoso, movimentos laterais e também a rotação, acoplamento

giroscópico, inércia de rotação, mancais e selos na direção radial.

Nas seções 4.2 e 4.3 serão apresentadas as funcionalidades antigas (versão 2015)

e novas (versão 2016) do programa respectivamente para mostrar a evolução do

programa.

4.2 LaviROT versão 2015

A primeira versão do programa oferecia ao usuário a possibilidade de criar uma

geometria de eixo para o rotor, selecionar dois tipos de materiais diferentes, aplicar

condições de rigidez e amortecimento para mancais e selos, posicionar discos e denotar

propriedades a eles. Além disso, o programa gerava duas respostas: diagrama de

Campbell e a resposta no domínio da frequência. A Figura 4.1 abaixo ilustra a interface

original do LaviROT.

Figura 4.1: Interface do LaviROT 2015.

14

4.3 LaviROT versão 2016

A versão 2016 do LaviROT visa não apenas ampliar a gama de análises que o

usuário pode realizar, mas também a variedade de dados de entrada e condições de

contorno para possibilitar a simulação de diferentes situações. A Figura 4.2 mostra a

interface do LaviROT versão 2016.

Figura 4.2: Interface do LaviROT 2016.

4.4 Novas opções de análises

Foram adicionadas ao programa as seguintes opções de análises: Análise de

Estabilidade, Diagrama de Nyquist e do regime transiente e permanente do domínio do

tempo.

4.4.1 Análise de estabilidade

Esta análise se baseia no problema de autovalor de um sistema de n graus de

liberdade. Ela visa verificar se o sistema é estável ou não. Isso feito a partir da verificação

da parte real dos autovalores da matriz abaixo:

[𝐴(Ω)] = |[0] [𝐼]

[𝑀]−1 ∙ [𝐾] [𝑀]−1 ∙ ([𝐶] + [𝐺])| (4.1)

onde I é a matriz identidade, M a matriz massa, K a matriz de rigidez, C a de

amortecimento e G a matriz de efeito giroscópico. Se a parte real de algum autovalor for

positiva, o sistema é instável. Se todos os autovalores tiverem parte real negativa, o

sistema será estável.

15

4.4.2 Diagrama de Nyquist

Um gráfico de Nyquist é outro nome dado a um diagrama polar de uma função. O

diagrama polar de uma função de transferência senoidal 𝐺(𝑗𝜔) é um gráfico do módulo

de 𝐺(𝑗𝜔) versus o ângulo de fase de 𝐺(𝑗𝜔) em coordenadas polares com 𝜔 variando de

0 a infinito [7].

No programa, é gerado um diagrama de Nyquist para o usuário. O programa

calcula a partir da representação de espaço de estados de um sistema linear (equação 4.2

e equação 4.3).

(𝑡) = 𝐴(𝑡)𝑥(𝑡) + 𝐵(𝑡)𝑢(𝑡) (4.2)

𝑦(𝑡) = 𝐶(𝑡)𝑥(𝑡) + 𝐷(𝑡)𝑢(𝑡) (4.3)

onde 𝑥(𝑡) é o vetor de estado, 𝑦(𝑡) é o vetor de saída (output), 𝑢(𝑡) é o vetor de entrada

(input), 𝐴(𝑡) é a matriz de estados, 𝐵(𝑡) é a matriz de entrada (input), 𝐶(𝑡) é a matriz de

saída (output), 𝐷(𝑡) é a matriz de transmissão direta e S é um vetor que localiza os nós

onde há forçamento.

𝐴 = [[0]𝑛𝑥𝑛 [𝐼]𝑛𝑥𝑛

[𝑀]𝑛𝑥𝑛−1 ∙ [𝐾]𝑛𝑥𝑛 [𝑀]𝑛𝑥𝑛

−1 ∙ ([𝐶] + [𝐺])𝑛𝑥𝑛] (4.4)

𝐵 = [

[𝑍]𝑛𝑥𝑛 ∙ [𝑆]𝑛𝑥1

[𝑀]𝑛𝑥𝑛−1 ∙ [𝑆]𝑛𝑥1

] (4.5)

𝐶 = [

00⋮

1000

]

2𝑛𝑥1

(4.6)

𝐷 = 0 (4.7)

4.4.3 Análise do regime transiente e permanente do domínio do tempo

Esta análise gera gráficos de amplitude de vibração no tempo e permite observar

os regimes transiente e permanente. O programa realiza os cálculos a partir da equação

4.8 [4], utilizando o método de iteração trapezoidal ODE23t do Matlab.

= [[0] [𝐼]

−[𝑀]−1 ∙ [𝐾] −[𝑀]−1 ∙ ([𝐶] + [𝐺])] ∙ 𝑄 + [

[0]

[𝑀]−1 ∙ [𝐹]] (4.8)

onde 𝐺 é a matriz giroscópica e tem as mesmas dimensões da matriz de massa (𝑀𝑛𝑥𝑛), e

𝐹 é o vetor de força (𝐹𝑛𝑥1), que tem o mesmo número de linhas da matriz de massa. Essa

análise exige uma configuração prévia, que será explicada a seguir

No âmbito das entradas de dados, foram adicionadas duas tabelas distintas, uma

para a inserção de forçamentos harmônicos e outra de forças de pulso. Na primeira,

entram valores dos módulos das forças harmônicas (𝑓𝑦 e 𝑓𝑧) nas direções Y e Z do sistema

16

de coordenadas e suas respectivas frequências de excitação, w, e ângulos de fase, 𝜑 (já

existe uma defasagem padrão de 90° entre as duas direções). Esse forçamento estará

presente desde o tempo inicial de análise até o tempo final. As forças harmônicas são

definidas pelas equações abaixo:

𝐹𝑦 = 𝑓𝑦 ∙ cos (𝑤 ∙ 𝑡 + 𝜑) (4.9)

𝐹𝑧 = 𝑓𝑧 ∙ 𝑠𝑒𝑛(𝑤 ∙ 𝑡 + 𝜑) (4.10)

Na segunda, entram valores de forças na forma de pulso nas direções Y e Z. Esta

função não utiliza delta de Dirac para calcular a resposta ao pulso. Essa função aplica

uma força durante um período de tempo muito curto (na ordem de 10−4s) definido pelo

passo da análise transiente, e o usuário poderá escolher em que instante de tempo o pulso

ocorre.

Vale ressaltar que é possível adicionar todos os tipos de forçamentos numa mesma

simulação. Os forçamentos serão adicionados à matriz de F de forças na equação 4.8.

Todos os vetores de força citados são criados automaticamente a partir da entrada

de dados de forças imposta pelo usuário. No código do programa, foi inserida uma rotina

que interpreta os nós com carregamento e aplica esses carregamentos nos graus de

liberdade correspondentes a cada nó.

Em conjunto a opção de análise transiente foi inserido no programa um painel de

configurações desta opção (Figura 4.3). Nela, é possível escolher a velocidade de rotação

do eixo, a qual será realizada a simulação, ativar ou não o efeito da gravidade, escolher o

nó de interesse da geometria para gerar os gráficos resposta (seja nos nós dos discos,

mancais, selos ou onde há forças atuando), opção de o programa calcular as órbitas para

cada um desses nós citados e o intervalo de tempo da análise e seu respectivo passo (esses

dois últimos itens são apenas parâmetros para geração de gráficos, o passo inserido não é

o mesmo passo utilizado pelo comando ode23t do Matlab, uma vez que esta função utiliza

um passo adaptativo).

17

Figura 4.3: Detalhe do quadro de configuração da análise transiente.

O efeito da gravidade foi inserido de tal forma que o peso de cada elemento (e não

de cada região) é dividido pelos seus respectivos nós. Esse peso é calculado através do

comprimento, diâmetro e massa específica associados a cada elemento discretizado e da

aceleração da gravidade. Cada nó da malha terá metade do peso de seus elementos

associados, por exemplo, o nó 2, que é compartilhado pelos elementos 1 e 2, terá uma

força aplicada nele que equivale à metade do peso do elemento 1 somado a metade do

peso do elemento 2.

18

5 COMPARAÇÃO ANSYS VS LaviROT

O programa LaviROT é uma ferramenta ainda em desenvolvimento. Por este

motivo, ele ainda deve passar por testes para verificar e validar seus resultados, se

provando confiável, antes que seja difundido em qualquer campo, seja no laboratório ou

para uso comercial. Este capítulo tem o propósito de criar um modelo de teste no

LaviROT e em um outro software comercial difundido, no caso o ANSYS, para comparar

os resultados obtidos em ambos os programas, e assim, validar os resultados do LaviROT.

5.1 Modelo de teste inicial

Como modelo de teste, utilizou-se a geometria do eixo do rotor que está em uma

bancada de testes no LAVI. No modelo, o eixo foi apoiado por dois mancais de rolamento

em suas extremidades. Além disso, um disco foi adicionado na parte central da geometria.

A fim de aproximar mais da realidade, tanto este teste quanto as simulações

posteriores, a rigidez direta dos mancais de rolamento foi calculada de acordo com a seção

5.2.

5.2 Cálculo da rigidez do rolamento

O deslocamento de um eixo, a partir de sua posição original concêntrica, depende

da deformação elástica das pistas e das esferas dos rolamentos que sustentam o eixo. Essa

deformação elástica depende da geometria e do material do rolamento e da folga interna

ou pré-carga (interferência). A folga é um fator importante porque determina a área de

tensão nas pistas. Quanto menor a folga (ou maior a interferência), mais elementos irão

compartilhar a carga externa aplicada.

Para uma dada força, o tamanho da área de contato determina a magnitude da

tensão nos componentes do rolamento. Pela teoria de contato de Hertz, para dois corpos

em contato, submetidos à uma força de compressão, o ponto de contato inicial entre eles

torna-se uma área (uma elipse para o caso de duas esferas). Essa força é dada pela equação

abaixo:

𝐹 = 𝐾𝑝 ∙ 𝛿𝑝3/2

(5.1)

Onde o subscrito p representa o ponto de contato da esfera com a pista (seja interna, no

ponto I ou externa, ponto O), 𝛿𝑝 é a deformação no ponto de contato e 𝐾𝑝 é a constante

força-deformação para um único ponto de contato, que depende da geometria e do

19

material da superfície [8, 9]. Nesse caso, há deformação tanto na pista externa quanto na

interna (Figura 5.1).

Figura 5.1: Pontos de contato da esfera com as pistas interna e externa (adaptado de [10]).

A deformação total é dada por:

𝛿 = 𝛿𝑝𝑖 ∙ 𝛿𝑝𝑜 (5.2)

𝛿 = (1

𝐾𝑝𝑖+

1

𝐾𝑝𝑜) ∙ 𝐹0,67 (5.3)

onde "i" representa a pista interna e "o" a pista externa. Assim, a equação 3 pode ser

rescrita como:

𝐹 = 𝐾𝑝𝑖𝑜 ∙ 𝛿1,5 (5.4)

onde 𝐾𝑝𝑖𝑜 é a constante força-deformação para o contato da esfera com as pistas dos

rolamentos em 2 pontos.

𝐾𝑝𝑖𝑜 =(𝐾𝑝𝑖 + 𝐾𝑝𝑜)

(𝐾𝑝𝑖0,67 + 𝐾𝑝𝑜

0,67)1,5 (5.5)

𝐾𝑝𝑖𝑜 depende principalmente da razão de curvatura f (Figura 5.2), do diâmetro da

esfera, do diâmetro do rolamento e do ângulo formado entre os pontos de contato da esfera

do rolamento e a vertical. A constante é definida pela equação 5.6.

𝑓 =2 ∙ 𝑟 + 𝐷𝑏

𝐷𝑏 (5.6)

onde r é o raio de curvatura e 𝐷𝑏 é o diâmetro da esfera do rolamento.

20

Figura 5.2: Representação do ângulo de contato (adaptado de [10]).

Segundo [11], ao desprezar a influência do diâmetro do rolamento e do ângulo de

contato alpha e para valores de f menores que 0,1 (maioria dos casos para rolamento de

esferas rígidas) tem-se a seguinte equação:

𝐾𝑝𝑖𝑜 =34300

𝑓0,35∙ 𝐷𝑏

0,5 (5.7)

A equação 5.4 indica a relação da deformação elástica com uma única força de

compressão aplicada ao elemento rolante. De fato, existe mais de um elemento rolante

suportando essa força de compressão, e neste caso, o efeito de cada elemento precisa ser

levado em consideração.

Quando a pista interna de um rolamento é deslocada de sua posição concêntrica

original à uma distância 𝑥𝑚 referente ao centro da pista externa, parte dessa distância

consiste na folga radial interna 𝑐𝑟. Assim, a deformação elástica na direção da carga radial

aplicada é:

𝛿0 = 𝑥𝑚 − 𝑐𝑟 (5.8)

E a deformação elástica em qualquer elemento rolante a partir de um ângulo 𝜑

com a vertical é dado por:

𝛿(𝜑) = 𝑥𝑚 ∙ cos (𝜑) ∙ 𝑐𝑟 (5.9)

Fazendo o equilíbrio de forças nas pistas do rolamento, tem-se:

21

𝐹𝑟 = ∑𝐹(𝜑𝑗)

𝑍

𝑗=1

(5.10)

Onde Z é o número de esferas do rolamento.

Para um dado rolamento, com as posições angulares das esferas conhecidas, um

deslocamento na pista interna 𝑥𝑚 pode ser assumido, resultando numa deformação

elástica 𝛿(𝜑), em cada esfera, calculada pela equação 8.

Se a fração de carga radial que é transmitida através da esfera diretamente em

linha com a carga aplicada é conhecida, então, o deslocamento resultante 𝑥𝑚 da pista

interna pode ser calculado diretamente. Da equação 5.4 tem-se:

𝑥𝑚 = 𝛿0 + 𝑐𝑟 (5.11)

Com:

𝛿0 = (𝐹𝑚

𝐹𝑝𝑖𝑜)

0,67

(5.12)

Onde 𝐹𝑚 é a força na esfera diretamente em linha com a carga radial aplicada.

Palmgren [11] faz aproximações negligenciando efeitos de folga e geometria do

rolamento:

𝑥𝑚 = 4,36 ∙ 10−8 ∙ (𝐹𝑚

2

𝐷𝑏)

0,33

(5.13)

Antes da equação 5.11 ser aplicada, a máxima força de compressão, 𝐹𝑚, em um

único elemento rolante deve ser determinada. A relação entre 𝐹𝑚 e 𝐹𝑟 é dada por:

𝐹𝑚 =4,37

𝑍∙ 𝐹𝑟 (5.14)

Substituindo a equação 5.14 na equação 5.11 tem-se:

𝑥𝑚 − 𝑐𝑟 = (4,37 ∙ 𝐹𝑟

𝑍 ∙ 𝐾𝑝𝑖𝑜)

11,08

(5.15)

Por fim, a relação entre a força radial e deformação elástica do rolamento foi

estabelecida e a rigidez pode ser calculada. Nota-se que a rigidez é não-linear com o

deslocamento da pista interna em relação à externa.

Substituindo 𝐹 = 𝐹𝑚 e 𝛿 = 𝑥𝑚 − 𝑐𝑟 na equação 5.4, tem-se:

𝐹𝑚 = 𝐾𝑝𝑖𝑜 ∙ (𝑥𝑚 − 𝑐𝑟)1,5 (5.16)

A rigidez do rolamento 𝐾(𝑥𝑚) pode ser obtida diferenciando a equação 5.16 em

relação ao deslocamento 𝑥𝑚:

22

𝐾(𝑥𝑚) =𝑑

𝑑𝑥𝑚𝐹𝑟 (5.17)

𝐾(𝑥𝑚) =𝑍

4,37∙

𝑑

𝑑𝑥𝑚𝐹𝑚 (5.18)

𝐾(𝑥𝑚) = 1,5 ∙𝑍

4,37∙ 𝐾𝑝𝑖𝑜 ∙ (𝑥𝑚 − 𝑐𝑟)

0,5 (5.19)

Para o caso da bancada de teste, o rolamento utilizado é o rolamento de esferas

SKF 6304-2Z. Pelo catálogo da SKF [12] tem-se os valores 𝐷 = 9,55mm; Z = 7 esferas;

𝐷𝑏 = 9,52mm. A folga radial da pista interna é de 0,01mm e a deformação elástica 𝑥𝑚

foi estimada em 0,008mm. Para esses valores, o resultado da equação 19. Desta forma, o

valor da rigidez do rolamento é de 1,771 ∙ 108 N/m.

5.3 Simulação em ANSYS Workbench

5.3.1 Introdução

O software de simulação por método de elementos finitos, ANSYS, apresenta uma

ferramenta para estudo de rotodinâmica. Ela é um complemento da análise modal e

permite não só a observação gráfica dos modos naturais de vibração do modelo, mas

também dos valores das frequências e a criação de um diagrama de Campbell.

Através de uma análise transiente, é possível gerar um gráfico de resposta em

frequência para a amplitude do deslocamento. E em uma análise estática, calcula-se o

deslocamento gerado pelo pulso de dois atuadores.

5.3.2 Geometria

O eixo foi desenhado na própria plataforma de desenho do ANSYS (Design

Modeler). A fim de gerar uma geometria com características semelhantes ao modelo feito

no LaviROT e que o programa gerasse elementos de viga automaticamente, a geometria

foi feita a partir de pontos, linhas de corpo (line bodies) entre estes pontos e seções

circulares. Todas as dimensões do desenho estão de acordo com o desenho técnico no

Apêndice A. A geometria final é apresentada abaixo nas Figuras 5.3 e 5.4.

23

Figura 5.3: Geometria do eixo gerada no ANSYS.

Figura 5.4: Geometria do eixo com as seções visíveis.

5.3.3 Geração de malha

Para analisar o eixo, foi escolhido o elemento de viga que se assemelhasse ao

usado no LaviROT. Sendo assim o tipo de elemento adotado foi o BEAM188. Segundo

o guia de usuário do ANSYS [13], esse elemento é adequado para análises lineares, desde

estruturas delgadas às vigas grossas. O elemento BEAM189 não foi utilizado porque é

um elemento adequado para análises não lineares. O BEAM188 é baseado na teoria de

viga de Timoshenko [14], e possui 6 graus de liberdade: rotações em 3 eixos e translações

em 3 eixos. Já o elemento EB-Ritto [5] é baseado na teoria de viga de Euler-Bernoulli.

Além do BEAM188, outros elementos estão presentes na análise, tais como MPC184,

para realizar a união entre os corpos de linha; COMBI214, que é um elemento

representando o sistema mola-amortecedor dos mancais; e MASS21, para representar a

massa do disco desbalanceado.

24

O tamanho médio de cada elemento foi selecionado de forma que o número total

de elementos BEAM188 fosse o mesmo da análise no LaviROT, 24 elementos. Por isso,

o tamanho escolhido foi de 46mm. Nota-se que não necessariamente todos os elementos

têm o mesmo comprimento, o próprio programa regula isto de forma aos elementos se

adequarem à geometria.

5.3.4 Análise rotodinâmica

Condições de contorno

Como condição de contorno, foram utilizados juntas do tipo rolamento (bearing)

no segundo e no penúltimo nó da geometria (a rigidez dos rolamentos foi calculada como

mostra a seção 5.2; velocidade de rotação numa faixa variando de 1000RPM até

10000RPM; e um ponto de massa no centro do eixo representando o disco desbalanceado.

A localização e as propriedades de rigidez dos mancais e localização e as propriedades

de inércia e desbalanceamento do disco são mostradas nas figuras 5.5, 5.6 e 5.7.

Posteriormente, a velocidade de rotação é removida para calcular as frequências naturais

do eixo parado.

Figura 5.5: Detalhe das condições de contorno da análise rotodinâmica no ANSYS.

25

Figura 5.6: Detalhe dos dados do rolamento da direita.

Figura 5.7: Detalhe dos dados do rolamento da esquerda.

Resultados

Como resultados, foram gerados o diagrama de Campbell (Figura 5.8), o diagrama

de amplitude em função da velocidade de rotação (Figura 5.9), e as frequências naturais

do eixo sem a rotação (Figura 5.10).

Figura 5.8: Diagrama de Campbell gerado pelo ANSYS.

26

Figura 5.9: Gráfico Amplitude (mm) x Velocidade de rotação [RPM].

Figura 5.10: Frequências naturais do sistema geradas no ANSYS.

A Figura 5.11, gerada pelo ANSYS, abaixo mostra os valores das velocidades

críticas de cada modo, bem como a direção de precessão e se o modo é estável ou não.

Figura 5.11: Velocidades de rotação críticas e direção de precessão geradas no ANSYS.

5.3.5 Análise Modal


Para calcular a frequência natural do corpo livre de apoios ou forçamentos

externos, foi feita uma análise modal com o eixo sem as condições de contorno aplicadas

anteriormente. O objetivo disso é comparar o resultado com o valor experimental, da

primeira frequência natural, obtido no laboratório com o auxílio de um osciloscópio e um

martelo de impacto.

27

Resultados

A Figura 5.12 abaixo mostra os valores da primeira frequência natural para as duas

direções transversais do eixo.

Figura 5.12: Frequências naturais do eixo livre geradas no ANSYS.

5.3.6 Análise estática


Nesta análise, o objetivo foi calcular o deslocamento gerado pelo pulso de uma

força em 4 passos de 5N, 10N, 15N e 20N para simular o pulso dos atuadores que estarão

na bancada de teste. Essa simulação não levou em conta a velocidade de rotação do eixo.

Para isso, usou-se as mesmas condições de fixação por rolamento descritas no sub-seção

5.3.4 e aplicou-se duas forças com aumento gradativo do seu módulo do sentido negativo

de Y, nos locais demonstrados pela Figura 5.13 abaixo.

Figura 5.13: Detalhe do posicionamento das forças e dos seus respectivos módulos.

Resultados

A Figura 5.14 abaixo mostra o resultado da simulação. Nota-se que o

deslocamento máximo gerado pelo atuador é de aproximadamente 0,09mm

aproximadamente, na região central do eixo, para uma força de 20N, bem abaixo do limite

máximo de 0,3mm permitido pelo batente da carcaça.

28

Figura 5.14: Resultado da análise estática gerada no ANSYS.

5.4 Simulação em Matlab – LaviROT

5.4.1 Introdução

No ambiente do LaviROT, a geometria foi dividida em 13 regiões, como mostra

a Figura 5.15. Cada região contendo seus dados de massa específica, coeficiente de

Poisson, módulo de elasticidade, comprimento e diâmetro. Os valores dessas

propriedades são apresentados na Tabela 5.1.

Figura 5.15: Ilustração da separação das regiões do eixo.

Tabela 5.1: Propriedades aplicadas ao eixo na simulação no LaviROT.

Região Mod. Elást.

(Pa)

Massa

esp.

(kg/m³)

Coef.

Poisson

Comprimento

(m)

Diâmetro

(m)

1 2,000e+11 7850 0,29 0,015 0,006

2 2,000e+11 7850 0,29 0,0139 0,015

3 2,000e+11 7850 0,29 0,23345 0,015

4 2,000e+11 7850 0,29 0,0485 0,018

5 2,000e+11 7850 0,29 0,05275 0,020

6 2,000e+11 7850 0,29 0,065 0,040

29

7 2,000e+11 7850 0,29 0,065 0,040

8 2,000e+11 7850 0,29 0,05275 0,020

9 2,000e+11 7850 0,29 0,0485 0,018

10 2,000e+11 7850 0,29 0,23345 0,015

11 2,000e+11 7850 0,29 0,01 0,015

12 2,000e+11 7850 0,29 0,0411 0,010

5.4.2 Condições de contorno

Como condições de contorno, nos nós 2 e 12 foram adicionadas rigidezes diretas

(𝐾𝑥𝑥 e 𝐾𝑦𝑦) na tabela de mancais e selos com os mesmos valores calculados na seção 5.2.

Na tabela de discos foi inserido um disco no nó 7 com 3kg de massa, 80mm de diâmetro

e 0,03mm de excentricidade. A Figura 5.16 mostra a interface do programa com os dados.

Figura 5.16: Interface do LaviROT com os dados para a simulação de teste.

As análises de interesse para o comparativo são: Campbell, Bode e Transiente. A

velocidade máxima de rotação (opção RPM MAX no LaviROT) configurada para realizar

o diagrama de Campbell e Bode foi 10.000 RPM. Para a análise transiente, a velocidade

de rotação escolhida foi 1458 RPM. O motivo da escolha desta velocidade será explicado

na seção de resultados a seguir. E por fim, optou-se por observar os resultados no nó onde

o disco está localizado, devido à hipótese de este nó apresentar a maior amplitude de

deslocamento.

30

5.4.3 Resultados

A partir desta simulação, obteve-se os diagramas de Campbell (Figura 5.17), o

gráfico da amplitude em função da frequência (Figura 5.18) o gráfico da amplitude, do

nó onde está o disco, ao longo do tempo (Figura 5.19) e sua respectiva órbita (Figura

5.20).

Figura 5.17: Diagrama de Campbell da simulação teste gerado no LaviROT.

Figura 5.18: Gráfico da amplitude vs frequência para a simulação teste gerado no LaviROT.

31

Figura 5.19: Resposta transiente da simulação teste gerado no LaviROT.

Figura 5.20: Órbita da resposta transiente da simulação teste gerado no LaviROT.

O gráfico da análise transiente foi gerado após o resultado das frequências críticas

pois era desejado utilizar a frequência natural para gerar a resposta transiente com ela.

Então, pelo diagrama de Campbell foi observado o valor de 1458 RPM e então colocado

no programa para calcular a resposta transiente.

32

Pelos resultados, percebe-se que os gráficos apresentam uma coerência entre

velocidade crítica e amplitude máxima, ou seja, a velocidade crítica no diagrama de

Campbell coincide com a frequência correspondente ao pico de amplitude no gráfico da

amplitude em função da frequência; e a análise transiente mostra um comportamento de

um sistema ressonância.

5.5 Comparação dos resultados

Pelos resultados gerados no ANSYS e no LaviROT, pode-se notar que ambos

apresentaram o mesmo valor para a primeira velocidade crítica de rotação no diagrama

de Campbell, mas havendo diferença no valor da segunda frequência crítica. Ambos os

programas mostraram o pico de amplitude em relação a frequência crítica para o mesmo

valor de velocidade de rotação (1458 RPM). A principal hipótese para essa diferença é o

uso de diferentes modelos de viga (Euler-Bernoulli e Timoshenko), uma vez que no

modelo de Euler-Bernoulli, a deformação devida ao esforço cortante é desprezada,

enquanto no modelo do Timoshenko não. Por este motivo, o primeiro modelo é mais

rígido e as frequências naturais são maiores.

Em relação aos valores de amplitude os programas apresentam resultados

semelhantes comparando o gráfico da Figura 5.9 e os gráfico das Figuras 5.18 e 5.19. O

primeiro apresenta o valor de 0,1183mm para a amplitude máxima, o segundo 0,110mm

e o terceiro 0,125mm aproximadamente. O erro relativo entre esses resultados é de 5,4%,

que é um erro pequeno e corresponde a erros numéricos na simulação.

Por esta validação fica claro que o LaviROT é um programa confiável para a

simulação de rotores e é de fato uma alternativa viável a outros softwares do mercado.

33

6 SIMULAÇÃO DA BANCADA DE TESTES

6.1 Introdução

Neste capítulo, será simulado o rotor em bancada de testes do LAVI. A estratégia

geral adotada para tal foi utilizar o LaviROT, a geometria do eixo da Figura 5.15 com os

dados da Tabela 5.1, mancais posicionados nos nós 2 e 12, com valor de rigidez direta

calculada na seção 5.2, taxa de amortecimento de 0,05 e rotação máxima de 15000 RPM.

As análises transientes foram realizadas para velocidades 8000 RPM, com um passo de

0,0002s e o efeito da gravidade foi considerado também.

Como primeiro caso, calcularam-se as amplitudes para o sistema sem selos,

levando em consideração, além dos mancais, apenas dois discos nos nós 6 e 8 da

geometria. O objetivo desta primeira análise é servir de referência para comparação com

as próximas simulações.

Nos casos seguintes, selos foram posicionados nos nós vizinhos aos discos (nós 5

e 9). Como ainda não há informações concretas sobre valores de coeficientes de rigidez

e amortecimento dos selos que serão utilizados na bancada, esses dados foram estimados

inicialmente por [15] e depois esses coeficientes foram alterados a fim de observar a

influência desses parâmetros dos selos no comportamento do rotor. Também aplicou-se

forças de pulso para simular os atuadores que estarão presentes na bancada e observar se

os deslocamentos respeitam os limites da bancada de testes. Os dados dos discos são

apresentados na Tabela 6.1 e os dados de referência dos selos são apresentados na Tabela

6.2.

Tabela 6.1: Dados dos discos.

Disco Localização Massa Diâmetro Excentricidade

1 Nó 6 5 Kg 100 mm 0,03 mm

2 Nó 8 5 Kg 100 mm 0,03 mm

Tabela 6.2: Dados dos selos.

Selo Localização Rigidez direta Rigidez

acoplada

Amortecimento

direto

1 Nó 5 100000 N/m 80000 N/m 250 N.s/m

2 Nó 9 100000 N/m 80000 N/m 250 N.s/m

34

A carcaça do rotor tem um batente interno (touchdown) que apenas permite um

deslocamento de 0,3mm. Logo, qualquer resultado que mostre uma amplitude de vibração

superior a este valor, significa que haverá um choque do rotor com a carcaça.

6.2 Análise do rotor sem selos

Neste primeiro caso, o modelo simulado levou em consideração o eixo, com os

dois mancais, dois discos com as propriedades da Tabela 6.1, nenhum selo posicionado e

as demais condições citadas na seção 6.1.

Como resultados iniciais, obteve-se o diagrama de Campbell (Figura 6.1), o

gráfico da amplitude em função da frequência para os nós dos discos 1 e 2 (Figuras 6.2 e

6.3 respectivamente) e os gráficos da amplitude em função do tempo para os nós dos

discos 1 e 2 (Figuras 6.4 e 6.6) e suas órbitas (Figuras 6.5 e 6.7).

Figura 6.1: Diagrama de Campbell do primeiro modelo.

35

Figura 6.2: Gráfico da amplitude vs frequência para disco 1 do primeiro modelo.

Figura 6.3: Gráfico da amplitude vs frequência para disco 2 do primeiro modelo.

36

Figura 6.4: Análise transiente do disco 1 do primeiro modelo.

Figura 6.5: Órbita do disco 1 do primeiro modelo.

37

Figura 6.6: Análise transiente do disco 2 do primeiro modelo.

Figura 6.7: Órbita do disco 2 do primeiro modelo.

Pelo gráfico da Figura 6.1, tem-se as seguintes velocidades de rotação críticas e

direções de precessão (Tabela 6.3).

38

Tabela 6.3: Velocidades críticas e direções de precessão do primeiro modelo.

Modos críticos Velocidade crítica (RPM) Direção de precessão

1 949 Retrógrada

2 950 Direta

3 4900 Retrógrada

4 5330 Direta

Os gráficos das Figuras 6.2 e 6.3 mostram um pico de 0,00025m ou 0,25mm de

amplitude para o primeiro modo crítico. Esse valor ainda respeita o limite de 0,3mm

disponível na carcaça.

Já pelos gráficos da análise transiente mostram que o primeiro pico do regime

transiente atinge um valor de 0,43mm. Vale lembrar que, neste modelo, o rotor acelera

de 0 RPM até 8000 RPM instantaneamente, por isso as amplitudes tendem a ser mais

elevadas. Outro aspecto a ser observado é o tempo que o sistema leva até atingir o regime

permanente. Neste caso, percebe-se que o sistema leva aproximadamente 1s.

6.3 Análise do rotor com selos e sem força de pulso

Neste segundo caso, o modelo simulado levou em consideração o eixo, com os

dois mancais, dois discos e dois selos nas posições mencionadas na seção 6.1 e com as

propriedades das Tabelas 6.1 e 6.2.

Para observar o comportamento desse modelo, foram gerados os gráficos das

Figuras 6.8 (Diagrama de Campbell), 6.9 e 6.10 (amplitude em função da frequência para

cada disco).

Figura 6.8: Diagrama de Campbell do segundo modelo.

39

Figura 6.9: Gráfico da amplitude em função da frequência para o disco 1 do segundo modelo.

Figura 6.10: Gráfico da amplitude em função da frequência para o disco 2 do segundo modelo.

Além desses resultados, foram gerados também gráficos que mostram as

amplitudes de vibração nas direções Y e Z no regime transiente para os nós que

apresentam discos e suas respectivas órbitas, como mostram as Figuras 6.11, 6.12, 6.13 e

6.14.

40

Figura 6.11: Análise transiente para o disco 1 do segundo modelo.

Figura 6.12: Órbita do disco 1 do segundo modelo comparada com o do primeiro.

41

Figura 6.13: Análise transiente para o disco 2 do segundo modelo.

Figura 6.14: Órbita do disco 2 do segundo modelo comparada com o do primeiro.

Pelo gráfico da Figura 6.8, tem-se as seguintes velocidades de rotação

críticas e direções de precessão (Tabela 6.4).

42

Tabela 6.4: Velocidades críticas e direções de precessão do segundo modelo.

Modos críticos Velocidade crítica (RPM) Direção de precessão

1 1086 Retrógrada

2 1872 Direta

3 5087 Retrógrada

4 5743 Direta

Pelos resultados obtidos nas Figuras 6.9 e 6.10, nota-se que este limite de batente

interno é sempre respeitado (no regime permanente) ao longo de toda a faixa de

velocidades de rotação e que o selo reduz quase pela metade a amplitude de deslocamento.

Observando os resultados da análise transiente, percebe-se que a inclusão dos

selos faz com que o sistema se comporte de maneira muito diferente do sistema da seção

6.2. Pelos gráficos das órbitas no regime permanente (Figuras 6.12 e 6.14) nota-se que,

com a adição dos selos, o centro de giro dos nós dos discos translada 0,04mm no sentido

positivo de Y, e a amplitude máxima na direção Z foi reduzida de 0,2mm para 0,114

(redução de 44%). E pelos gráficos das Figuras 6.11 e 6.13 tem-se que o maior pico de

amplitude no regime transiente é de 0,23mm, e o tempo até atingir o regime permanente

é aproximadamente 0,9s, uma redução de 47% na amplitude e 10% no tempo em relação

aos resultados da seção 6.2.

6.4 Análise do rotor com selos e com força de pulso

Nessa análise foi considerada a situação de uma força de pulso de 20N atuando

nas direções Y e Z. Como já foi dito anteriormente, o objetivo desta análise é simular o

atuador que estará presente na bancada. O pulso foi dado no nó central (nó 7) da

geometria, e em um instante de tempo em que o regime do sistema já estivesse

permanente (foi adotado o instante 1s). As Figuras 6.15 e 6.16 mostram os resultados das

amplitudes nos nós dos discos.

43

Figura 6.15: Análise transiente para o disco 1 do terceiro modelo.

Figura 6.16: Análise transiente para o disco 2 do terceiro modelo.

44

Os resultados mostram que o pulso de 20N de magnitude não representou um

aumento significativo da amplitude do deslocamento (apenas 0,04mm), somente criando

um novo regime transiente para o sistema.

6.5 Análise do rotor com selos de diferentes rigidezes diretas

Nesta seção, é avaliada, por uma análise discreta, a influência do coeficiente de

rigidez direta 𝐾𝑦 e 𝐾𝑧 dos selos no comportamento do sistema, sendo que os demais

coeficientes se mantêm constante. Para isso variou-se o coeficiente de rigidez direta dos

selos de 100.000N/m à 500.000N/m, e assim gerou-se um gráfico (Figura 6.17) que

mostra a evolução da amplitude máxima no regime transiente em função desse

coeficiente, outro para o regime permanente (Figura 6.18), um terceiro que mostra a

duração do regime transiente (Figura 6.19) e por último a influência sobre as velocidades

de rotação críticas (Figura 6.20).

Figura 6.17: Variação da amplitude máxima (em módulo) no regime transiente em função da

rigidez direta do selo.

45

Figura 6.18: Variação da amplitude máxima (em módulo) no regime permanente em função da

rigidez direta do selo.

Figura 6.19: Tempo de duração do regime transiente em função da rigidez direta do selo.

46

Figura 6.20: Variação das frequências críticas de rotação em função da rigidez direta do selo.

Os resultados mostram que o aumento da rigidez direta dos mancais contribui para

a redução das amplitudes de vibração de modo não-linear (quase exponencial).

Observando o gráfico da Figura 6.20 pode-se afirmar que as quatro primeiras

frequências críticas de rotação do rotor variam de forma diretamente proporcional em

relação à 𝐾𝑦𝑧 e 𝐾𝑧𝑦.

6.6 Análise do rotor com selos de diferentes rigidezes acopladas


rigidez acoplado 𝐾𝑦𝑧 e 𝐾𝑧𝑦 dos selos no comportamento do sistema, sendo que os demais

coeficientes se mantêm constante. Para isso variou-se o coeficiente de rigidez acoplado

dos selos de 10.000N/m à 100.000N/m, e assim gerou-se um gráfico (Figura 6.21) que

mostra a evolução da amplitude máxima no regime transiente em função desse

coeficiente, um para o regime permanente (Figura 6.22) para a duração do regime

transiente (Figura 6.23) e para a influência sobre as velocidades de rotação críticas (Figura

6.24). Por fim, para a análise de estabilidade foram gerados dois gráficos, um apontando

a partir de qual valor de 𝐾𝑦𝑧 e 𝐾𝑧𝑦 o sistema se torna instável, através das amplitudes de

47

vibração (Figura 6.25), e outro, mostrando os diagramas de Nyquist para diferentes

valores desses parâmetros (Figura 6.26).

Figura 6.21: Variação da amplitude máxima no regime transiente em função da rigidez acoplada

do selo.

Figura 6.22: Variação da amplitude máxima no regime permanente em função da rigidez acoplada

do selo.

48

Figura 6.23: Tempo de duração do regime transiente em função da rigidez acoplada do selo.

Figura 6.24: Variação das frequências críticas de rotação em função da rigidez acoplada do selo.

49

Figura 6.25: Gráfico amplitude em função da rigidez acoplada do selo para apontar instabilidade.

50

Figura 6.26: Diagramas de Nyquist para diferentes valores de rigidezes acopladas do selo.

Os resultados mostram que a rigidez acoplada não influencia o sistema de maneira

semelhante à rigidez direta. Comparando as Figuras 6.21 e 6.22, percebe-se que o

coeficiente de rigidez acoplada dos selos exerce papéis diferentes sobre o regime

transiente e o permanente. Enquanto que a amplitude máxima no regime transiente

permanece praticamente constante, no regime permanente, a amplitude máxima de

vibração tende a aumentar conforme os valores de 𝐾𝑦𝑧 e 𝐾𝑧𝑦 aumentam. Já em relação a

duração do regime transiente, percebe-se que o tempo de duração não apresenta variação

significativa em função desses parâmetros.

51

Para as velocidades críticas é possível observar comportamentos linear em relação

aos coeficientes de rigidez acoplados. A primeira e a segunda frequências críticas são

mais sensíveis a variação de 𝐾𝑦𝑧 e 𝐾𝑧𝑦, sendo que a primeira varia negativamente e a

segunda positivamente com o aumento de 𝐾𝑦𝑧 e 𝐾𝑧𝑦. Já as demais, não apresentam

variação significativa.

Os gráficos da Figura 6.25 mostram diagramas de Nyquist para 5 situações

diferentes. O primeiro (para 𝐾𝑦𝑧 = 80000N/m) é o diagrama para o segundo modelo (rotor

com selo com as propriedades da Tabela 6.2), o último (parte inferior da figura) é o

diagrama para o caso o qual a rigidez acoplada atinge seu valor crítico e o sistema torna-

se instável e os demais gráficos mostram diagramas para valores de 𝐾𝑦𝑧 próximos ao

valor crítico, mas ainda dado com estável.

A Figura 6.26 mostra, através das amplitudes de deslocamento, para qual valor de

𝐾𝑦𝑧 o sistema se torna instável. Pode-se perceber que quando a rigidez atinge,

aproximadamente, o valor de 160500N/m, a amplitude tende a aumentar muito,

caracterizando instabilidade.

6.7 Análise do rotor com selos de diferentes coeficientes de amortecimento

diretos


amortecimento direto 𝐶𝑦𝑧 e 𝐶𝑧𝑦 dos selos no comportamento do sistema, sendo que os

demais coeficientes se mantêm constante. Para isso variou-se o coeficiente de

amortecimento direto dos selos de 100N.s/m à 400N.s/m, e assim gerou-se um gráfico

(Figura 6.27) que mostra a evolução da amplitude máxima no regime transiente em

função desse coeficiente, outro para o regime permanente (Figura 6.28) um terceiro que

mostra a duração do regime transiente (Figura 6.29).

52

Figura 6.27: Variação da amplitude máxima no regime transiente em função do amortecimento

direto do selo.

Figura 6.28: Variação da amplitude máxima no regime permanente em função do amortecimento

direto do selo.

Figura 6.29: Tempo de duração do regime transiente em função do amortecimento direto do selo.

53

Pelos gráficos das Figuras 6.27 e 6.28, pode-se observar que a variação do

coeficiente de amortecimento do selo não exerce influência significativa nas amplitudes

de deslocamento do eixo, pois elas permanecem praticamente constantes.

Já o gráfico 6.29 mostra que o tempo de duração do regime transiente do rotor

aumenta diretamente proporcional ao aumento do coeficiente de amortecimento direto do

selo. Essa tendência mostra um comportamento superamortecido. Porém, para confirmar

essa hipótese, seriam necessárias outras análises para o sistema.

54

7 CONCLUSÕES E TRABALHOS FUTUROS

Primeiramente, a comparação entre as simulações feitas no ANSYS e no

LaviROT, mostrou que o LaviROT é um programa confiável para análises de

rotodinâmica. A diferença numérica entre os resultados é muito pequena, correspondendo

apenas aos erros numéricos. Embora o programa ainda esteja em desenvolvimento, ele é

uma alternativa aos softwares comerciais que existem no mercado.

De maneira geral, a simulação da bancada de teste mostrou que em nenhum

momento haveria choque do rotor com a carcaça para os casos com selos presentes.

Os resultados da análise de variação da rigidez direta dos selos mostraram que as

amplitudes de vibração reduzem e as frequências críticas de rotação aumentam conforme

esse fator aumenta.

Os resultados da análise da variação da rigidez acoplada mostram que esse fator

atua de maneira diversa no sistema: a amplitude máxima no regime transiente permanece

praticamente constante e no regime permanente, ela tende a aumentar conforme os valores

de 𝐾𝑦𝑧 e 𝐾𝑧𝑦 aumentam. Em relação as velocidades críticas, foi possível observar

comportamentos linear em relação aos coeficientes de rigidez acoplados. Sendo a

primeira e a segunda frequências críticas mais sensíveis; A primeira variando

negativamente e a segunda positivamente com o aumento de 𝐾𝑦𝑧 e 𝐾𝑧𝑦. Já as demais, não

apresentam variação significativa. E em relação a análise de estabilidade, conclui-se que

o valor crítico para a rigidez acoplada do selo é de 160500N/m. Valores superiores a esse

tornam o sistema instável.

Já com relação ao amortecimento direto do selo, conclui-se que este fator não

exerceu influência significativa nas amplitudes máximas. Sua influência se deu no tempo

de duração do regime transiente, onde este aumenta conforme o amortecimento aumenta.

Para melhorar o programa, seria interessante incluir o modelo de viga de

Timoshenko, aprimorar a análise de estabilidade utilizando o critério de Nyquist ou até

aplicando outros métodos (Lyapunov, por exemplo) um modo torcional de vibração, uma

ferramenta de análise estática de deflexão do eixo. Outra sugestão é migrar a plataforma

base do LaviROT para uma plataforma gratuita (Python, por exemplo), pois o Matlab é

um software pago.

Por último, seria importante a comparação dos resultados gerados pelo LaviROT

com testes experimentais para aumentar ainda mais a confiabilidade do programa.

55

8 REFERÊNCIAS BIBLIOGRÁFICAS

[1] YOON, Se Young; LIN, Zongli; ALLAIRE, Paul E., Control of Surge in

Centrifugal Compressors by Active Magnetic Bearings: Theory and

Implementation., Springer Science Business Media, 2012

[2] NPTEL Mechanical Engineering (on-line), Disponível em: http://nptel.ac.in/courses/

[3] CRANDALL, S. H.; LARDNER, T. J.; ARCHER, R. R.; COOK, N. H.; DAHL, N.

C. An introduction to the mechanics of solids. McGraw-Hill, 1978

[4] CHILDS, Dara W, Turbomachinery rotordynamics: phenomena, modeling, and

analysis., John Wiley Sons, 1993.

[5] Ritto, Thiago G., Sampaio, Rubens., Rotor Dynamics with MATLAB - Finite

Element Analysis, PUC-Rio, 2007.

[6] INMAN, Daniel J.; SINGH, Ramesh Chandra., Engineering vibration., Upper

Saddle River: Prentice Hall, 2001.

[7] OGATA, Katsuhiko, Engenharia de Controle Moderno, 5ed. São Paulo, Pearson

Prentice Hall, 2010.

[8] CHANGSEN, W., Analysis of Rolling Element Bearings, Mechanical Engineering

Publications Ltd., London, 1991.

[9] HARRIS, T.A., Rolling Bearing Analysis, Wiley, New York, 2001.

[10] TIWARI, R., Analysis and Identification in Rotor-Bearing Systems, Indian

Institute of Technology Guwahati, 2010.

[11] PALMGREN, A., Ball and Roller Bearing Engineering, 3rd ed., Burbank, 1959.

[12] Catálogo SKF (on-line), Disponível em: <http://www.skf.com/>

[13] ANSYS Mechanical APDL Element Reference, ANSYS Inc, 2012

[14] TIMOSHENKO, Stephen P.; GERE, James M.; Mecânica dos Sólidos, Vol. 1,

LTC, 1993

[15] CHILDS, D. W.; RAMSEY, Christopher., Seal-rotordynamic-coefficient test

results for a model SSME ATD-HPFTP turbine interstage seal with and without

a swirl brake., ASME, Transactions, Journal of Tribology, v. 113, 1991

56

APÊNDICE A – Desenho detalhado do eixo

Documents

SIMULAÇÃO NUMÉRICA DE UM ROTOR EM …monografias.poli.ufrj.br/monografias/monopoli10018281.pdf · these machines to evaluate their behavior, taking into account the presence of