Thiago G. Ritto
Análise de Vibrações deSistemas Lineares e
Não-Lineares no Contextoda Formulação Fraca,
Análise Modal eDecomposição deKarhunen-Loève
DISSERTAÇÃO DE MESTRADO
DEPARTAMENTO DE ENGENHARIAMECÂNICA
Programa de Pósgraduação em
Engenharia Mecânica
Rio de JaneiroSetembro de 2005
Thiago G. Ritto
Análise de Vibrações de SistemasLineares e Não-Lineares no Contexto da
Formulação Fraca, Análise Modal eDecomposição de Karhunen-Loève
Dissertação de Mestrado
Dissertação apresentada como requisito parcial paraobtenção do grau de Mestre pelo Programa de Pósgraduação em Engenharia Mecânica do Departamento deEngenharia Mecânica da PUCRio
Orientador: Prof. Rubens Sampaio
Rio de JaneiroSetembro de 2005
Thiago G. Ritto
Análise de Vibrações de SistemasLineares e Não-Lineares no Contexto da
Formulação Fraca, Análise Modal eDecomposição de Karhunen-Loève
Dissertação apresentada como requisito parcial paraobtenção do grau de Mestre pelo Programa de Pósgraduação em Engenharia Mecânica do Departamento deEngenharia Mecânica do Centro Técnico Cientíco daPUCRio.Aprovada pela Comissão Examinadora abaixoassinada.
Prof. Rubens SampaioOrientador
Departamento de Engenharia Mecânica PUCRio
Prof. Fernando Alves RochinhaDepartamento de Engenharia Mecânica UFRJ
Prof. Edson Luiz Cataldo FerreiraDepartamento de Matemática Aplicada UFF
Prof. Marcelo Túlio PiovanUniversidad Tecnológica Nacional, Facultad Regional
Bahía Blanca
Prof. José Eugênio LealCoordenador Setorial do Centro Técnico Cientíco
PUCRio
Rio de Janeiro, 9 de Setembro de 2005
Todos os direitos reservados. É proibida a reproduçãototal ou parcial do trabalho sem autorização da univer-sidade, do autor e do orientador.
Thiago G. RittoGraduouse em Engenharia Mecânica e ProduçãoMecânica pela Pontifícia Universidade Católica do Riode Janeiro (Rio de Janeiro, RJ). Trabalhou duranteum ano e meio na Companhia Siderúrgica de Tubarão(Vitória, ES) onde exerceu atividades relacionadas àmanutenção preditiva e análise de vibrações. Este tra-balho de mestrado gerou um artigo publicado no18th International Congress of Mechanical Engineering(COBEM 2005), [32]. Na presente data trabalha na em-presa FURNAS CENTRAIS ELÉTRICAS SA (Rio deJaneiro, RJ) onde faz parte da Coordenação de Quali-dade Total.
Ficha CatalográcaRitto, Thiago G.
Análise de Vibrações de Sistemas Lineares e Não-Lineares no Contexto da Formulação Fraca, AnáliseModal e Decomposição de Karhunen-Loève/ Thiago G.Ritto; orientador: Rubens Sampaio. Rio de Janeiro: PUCRio, Departamento de Engenharia Mecânica,2005.
214 f. ; 30 cm
Dissertação (mestrado) - Pontifícia UniversidadeCatólica do Rio de Janeiro, Departamento de EngenhariaMecânica.
Inclui referências bibliográcas.
1. Engenharia Mecânica - Teses. 2. Análise de Vi-brações. 3. Formulação fraca. 4. Método de Galerkin.5. Análise modal. 6. Base de Karhunen-Loève. I. Sam-paio, Rubens. II. Pontifícia Universidade Católica do Riode Janeiro. Departamento de Engenharia Mecânica. III.Título.
CDD: 621
Agradecimentos
Inicialmente e principalmente agradeço ao meu pai Ritto, à minha
mãe Nazareth e ao meu irmão Fabio, maiores incentivadores e fontes pri-
mordiais das minhas energias.
Agradeço ao meu orientador, Rubens, pelas inúmeras aulas sobre
a vida, aprendizados muito maiores do que os que estão contidos neste
trabalho.
Agradeço à Companhia Siderúrgica de Tubarão (CST) através do
gerente Rúben Pinasco, por conceder todo o suporte que eu precisava, do
engenheiro Álvaro Pio, pelo apoio, e dos engenheiros Jorge Pires e Weber
Batista, por me ensinarem muito sobre análise de vibrações aplicada na
indústria.
Agradeço ao Conselho Nacional de Desenvolvimento Cientíco e Tec-
nológico (CNPQ) pelo auxílio concedido no primeiro ano do mestrado.
Finalmente agradeço a todos os amigos de Vitória e do Rio e aos
colegas e professores do Departamento de Engenharia Mecânica da PUC-
Rio que me deram a força e o incentivo que eu precisava.
Resumo
Ritto, Thiago G.; Sampaio, Rubens. Análise de Vibrações deSistemas Lineares e Não-Lineares no Contexto da Formu-lação Fraca, Análise Modal e Decomposição de Karhunen-Loève. Rio de Janeiro, 2005. 214p. Dissertação de Mestrado Departamento de Engenharia Mecânica, Pontifícia UniversidadeCatólica do Rio de Janeiro.
Neste trabalho a Análise de Vibrações é tratada no contexto da formulação
fraca. Um sistema contínuo é formulado abstratamente em um espaço de
Hilbert e uma base de projeção é escolhida para a dinâmica. Um esquema de
convergência para a aproximação é garantido à medida em que se aumenta o
número de funções da base usada para representar a resposta do problema.
Esta é a idéia por traz de métodos como o Método dos Elementos Finitos e
o Método dos Modos Supostos, que derivam do Método de Galerkin. Esta
estratégia é diferente do que comumente é ensinado nos cursos de vibrações,
onde um sistema massamola é analisado, e sistemas discretos formados por
massas, molas e amortecedores são discutidos. Nestes casos não se sabe qual
é o erro cometido na análise numérica. A Análise de Vibrações é muito usada
na manutenção preditiva de máquinas rotativas. Alguns fenômenos obser-
vados nesses equipamentos motivaram o desenvolvimento de um modelo
numérico que pudesse reproduzir tais fenômenos para melhor entendê-los.
Um sistema rotormancal é modelado e sua resposta dinâmica comparada
qualitativamente com a resposta dinâmica captada através de acelerômetros
xados nos mancais de um exaustor da Companhia Siderúrgica de Tubarão
(CST). Durante o trabalho diversos programas foram desenvolvidos através
da plataforma MATLAB.
PalavraschaveAnálise de Vibrações, Formulação Fraca, Método de Galerkin, Análise
Modal, Base de Karhunen-Loève.
Abstract
Ritto, Thiago G.; Sampaio, Rubens. V. Rio de Janeiro, 2005.214p. MSc. Dissertation Departamento de Engenharia Mecânica,Pontifícia Universidade Católica do Rio de Janeiro.
Vibration Analysis is treated in the context of weak formulation. A con-
tinuous system is formulated in the Hilbert space and one base is selected
to project the dynamics. An approximation scheme is guaranteed by in-
creasing the number of functions in the base used to represent the response.
This is the idea behind methods like the Finite Element Method and As-
sumed Modes Method, which derive from Galerkin Method. This strategy
is dierent from what is commonly taught in vibration courses, where a
massspring system is analyzed and discrete systems composed by masses,
springs and dashpots are discussed. In those cases the error of the numerical
analysis is not known. Vibration Analysis is very used in predictive mainte-
nance of rotating machines. Some phenomenons observed in those machines
motivated the development of a numerical model that could reproduce such
phenomenons to better understand them. A rotorbearing system is mod-
elled and its dynamic response is qualitative compared to the dynamic re-
sponse captured by accelerometers xed on the bearings of a blower of the
steel company Companhia Siderúrgica de Tubarão (CST). During this work
several programs were developed using MATLAB software.
KeywordsVibration Analysis, Weak Formulation, Galerkin Method, Modal Ana-
lysis, Karhunen-Loève Basis.
Conteúdo
1 Introdução 151.1 Manutenção Mecânica 151.2 Motivação 171.3 Estruturação 18
2 Formulação Fraca 202.1 Introdução 202.2 Espaços de Hilbert 202.3 Equacionamento 252.4 Problema Modelo 282.5 Formulação Fraca 392.6 Operadores Diferenciais 432.7 Convergência das Aproximações 482.8 Conclusões 50
3 Escolha da Base de Projeção 513.1 Introdução 513.2 Método dos Modos Supostos 513.3 Método dos Elementos Finitos 593.4 Carregamento 723.5 Condições de Contorno 753.6 Conclusões 85
4 Análise Modal 864.1 Introdução 864.2 Modos Normais 874.3 Modos Complexos 934.4 Conclusões 97
5 Transformada Rápida de Fourier e Função Resposta em Freqüência 985.1 Introdução 985.2 Transformada Rápida de Fourier 985.3 Função Resposta em Freqüência 1065.4 Conclusões 113
6 Decomposição de KarhunenLoève 1146.1 Introdução 1146.2 Base de KL e Modos de Vibração 1186.3 Barra chocando-se contra um ostáculo 1256.4 Conclusões 131
7 Inuência do desbalanceamento e folga nos mancais de um rotor embalanço 132
7.1 Introdução 132
7.2 Motivação 1357.3 Equacionamento 1387.4 Simulação numérica 1427.5 Conclusões 148
8 Conclusões e Trabalhos Futuros 149
Referências Bibliográcas 152
A Análise de vibrações na internet 157
B Mais Exemplos 159
C Lista dos Programas Desenvolvidos 174
D Relação entre os Programas 178
E Manual dos Programas Desenvolvidos 184
F Arquivos do CALFEM 204
Lista de Figuras
1.1 Ponte de Tacoma, julho de 1940 151.2 Equipamentos relevantes em uma indústria: (a) Turbogerador;
(b) Compressor; (c) exaustor; (d) peneira vibratória. 16
2.1 Representação do Teorema da Projeção 252.2 Barra xa em um extremidade e livre na outra 282.3 Elemento diferencial de uma barra 282.4 Modo de corpo rígido de uma barra com as duas extremidades
livres 352.5 Dez primeiros modos de vibração de uma barra xa em uma
extremidade e livre na outra 372.6 Dez primeiros modos de vibração de uma barra xa nas duas
extremidades 382.7 Espaço de soluções 412.8 Erro ortogonal ao espaço formado pelas funções teste usadas na
aproximação 432.9 Algoritmo para convergência 50
3.1 Viga de EulerBernoulli 513.2 Viga engastadalivre 533.3 Viga engastada, excitada na extremidade livre 563.4 Resposta dinâmica na extremidade livre 573.5 Resposta dinâmica na extremidade livre 583.6 Aproximação linear 613.7 Polinômios de interpolação diferentes de zero apenas em uma
pequena parte do domínio 633.8 Elemento com aproximação (a) Linear, 2 nós, (b) Quadrática, 3
nós e (c) Cúbica, 4 nós 643.9 Aproximação quadrática 643.10 Elemento de viga 653.11 Funções Hermitianas 663.12 Viga engastada, excitada na extremidade livre 683.13 Resposta dinâmica em x=L 683.14 Resposta dinâmica em x=L 693.15 Eixos local e global para um elemento de barra 703.16 Elemento com três graus de liberdade em cada nó 713.17 Força concentrada 723.18 Força distribuída 733.19 Força concentrada na extremidade de uma viga 733.20 Força distribuída uniformemente ao longo de uma viga 743.21 Condições de contorno no extremo L. (a) ligação elástica; (b)
massa pontual 753.22 Cinco primeiros modos de uma barra xada em uma extremidade
e livre na outra 78
3.23 Convergência da aproximação 783.24 Cinco primeiros modos de uma barra xa em uma extremidade
e com apoio elástico na outra 803.25 Convergência da aproximação 803.26 Primeiro modo de vibração 813.27 Segundo modo de vibração 813.28 Cinco primeiros modos barra xamassa 833.29 Convergência da aproximação 833.30 Primeiro modo de vibração 843.31 Segundo modo de vibração 84
4.1 Separação dos parâmetros modais 864.2 Movimento de corpo rígido 884.3 Modo normal 914.4 Tentativa de visualização de um modo complexo 96
5.1 Freqüências naturais e modos de vibrações 995.2 Sinal real (esquerda) e sinal digitalizado (direita) 1005.3 (a) Sinal teórico (verde) e coletado (azul) (b) FFT dos sinais 1015.4 (a) Sinal teórico (verde) e coletado (azul) (b) FFT dos sinais 1025.5 (a) Sinal teórico (verde) e coletado (azul) (b) FFT dos sinais 1025.6 (a) Sinal teórico (verde) e coletado (azul) (b) FFT dos sinais 1035.7 Vazamento 1035.8 Janela Hanning 1045.9 Coleta de sinais (a) Acelerômetro; (b) Sensor de proximidade 1045.10 Programa usado para auxiliar a análise dos dados de vibrações 1055.11 Desbalanceamento 1055.12 Problema uidodinâmico 1065.13 FRF Teste com martelo de impacto 1085.14 FRF perto de uma freqüência natural (a) amplitude (b) Fase 1085.15 FRF Modos de vibração 1095.16 Medição da FRF em uma estrutura 1105.17 FRF do sistema 1115.18 Peneira vibratória industrial 1125.19 (a) FRF (b) Ângulo de fase 113
6.1 Sistema massamolaamortecedor com dois graus de liberdade 1166.2 Coecientes temporais a1(azul) e a2(verde) 1176.3 Reconstrução com 1 POM (a) Resposta x1 (b) Resposta x2 1186.4 Reconstrução com 2 POMs (a) Resposta x1 (b) Resposta x2 1186.5 Rotor com dois graus de liberdade, x1 e x2, e velocidade de
rotação constante, Ω 1216.6 Resposta dinâmica (a) x1 (b) x2 1226.7 Resposta dinâmica (a) x1 (b) x2 1236.8 Reconstrução (em azul) com apenas 1 POM 1236.9 Diagrama de Campbell para um disco 1246.10 Diagrama de Campbell para um cilindro 1256.11 Barra chocandose contra um obstáculo 125
6.12 (a) Convergência da aproximação, variando ∆t; (b) Convergên-cia da aproximação, variando N 127
6.13 Resposta do ponto da extremidade livre da barra 1276.14 (a) Força média x dist; (b) Força média x ka 1286.15 Modos de vibração x base de KL 1286.16 Aproximação da resposta dinâmica em x = L 1296.17 Aproximação da resposta dinâmica perto da região de choque 1306.18 Primeiro modo normal x Primeiro modo empírico 1306.19 Derivada segunda do primeiro modo normal x Derivada segunda
do primeiro modo empírico 131
7.1 Corte transversal do mancal 1357.2 Exaustor Principal da Máquina de Lingotamento Contínuo 1367.3 Espectro de freqüências. (a) 30 de janeiro de 2005, (b) 9 de
fevereiro de 2005 1367.4 Exaustor da dessulfuração 1377.5 Espectro de freqüências. (a) 2 de julho de 2001, (b) 23 de maio
de 2002 1377.6 Sistema rotodinâmico considerado 1397.7 Viga engastada em uma extremidade e com uma massa na outra,
discretizada pelo MEF 1407.8 (a) Convergência da aproximação, variando dt; (b) Convergência
da aproximação, variando N 1437.9 Conguração do eixo quando em operação 1447.10 Resposta dinâmica no mancal número 3. (a) Resposta transiente
u1 (azul) e u2 (verde); (b) Órbita após transiente 1447.11 Força média no mancal número 3 (azul) e 4 (verde) 1457.12 Força média no mancal número 4 1467.13 (a) Força média no mancal número 3 versus G; (b) Força média
no mancal número 4 versus G 1467.14 FFT do deslocamento no mancal número 3: G = 2 (azul) e
G = 20 (magenta) 1477.15 (a) Força média no mancal número 3 versus folga; (b) Força
média no mancal número 3 (azul) e 4 (verde) versus folga 148
8.1 Etapas da análise 1498.2 Algoritmo para convergência 1508.3 Análise Modal e Decomposição de KarhunenLoève 150
B.1 Viga engastada em uma extremidade 159B.2 Cinco primeiros modos de uma viga engastada em uma extrem-
idade 160B.3 Viga biapoiada 161B.4 Cinco primeiros modos de uma viga biapoiada 161B.5 Estrutura 162B.6 Estrutura discretizada 163B.7 Primeiro modo de vibração 165B.8 Cinco primeiros modos de uma Viga Engastadalivre 166B.9 Convergência da aproximação 167
B.10 Cinco primeiros modos de uma viga Biapoiada 168B.11 Convergência da aproximação 169B.12 Viga conguração L 169B.13 Cinco primeiros modos da conguração L 170B.14 Convergência da aproximação 171B.15 Viga conguração triângulo 171B.16 Cinco primeiros modos da conguração triângulo 172B.17 Convergência da aproximação 173
C.1 Diagrama de uxo 174
D.1 Primeiro grupo de programas 179D.2 Segundo grupo de programas 180D.3 Terceiro grupo de programas 181D.4 Quarto grupo de programas 182D.5 Simulação de um rotor 183
Lista de Tabelas
1.1 Comparação de custos entre manutenção corretiva, preventiva epreditiva 17
3.1 Comparação entre valores obtidos para o deslocamento 75
6.1 Dados usados no programa barra_choque 126
7.1 Parâmetros do programa rotor_choque 1407.2 Dados usados no programa rotor_choque 143
"A melhor maneira de prever o futuro é criá-lo"
Peter Drucker, Looking Ahead: Implications of the present.
1Introdução
1.1Manutenção Mecânica
Em todo o mundo, a todo tempo, máquinas estão funcionando e
estruturas estão sendo projetadas. Em um projeto de uma estrutura, por
exemplo, as características dinâmicas devem ser levadas em consideração
para que não ocorram surpresas inesperadas, como foi o caso da Ponte
de Tacoma (1940), no Estado de Washington (EUA), que entrou em
ressonância com rajadas de vento, Figura 1.1.
Figura 1.1: Ponte de Tacoma, julho de 1940
Vibração é uma oscilação em torno de um movimento de referência.
Na área de manutenção, a Análise de Vibrações é uma ferramenta ecaz e
necessária. Máquinas rotativas importantes para um processo, como uma
turbomáquina, por exemplo, são monitoradas através de sensores que
registram os níveis de vibração de pontos estratégicos (os mancais). Uma
máquina é desligada se os valores de vibração ultrapassarem um limite
estipulado, evitando, assim, falhas catastrócas e prejuízos astronômicos.
A Figura 1.2 mostra alguns equipamentos relevantes em uma indústria
como é o caso de um turbogerador, um compressor, um exaustor e uma
peneira vibratória.
16
(a) (b)
(c ) (d)
Figura 1.2: Equipamentos relevantes em uma indústria: (a) Turbogerador;
(b) Compressor; (c) exaustor; (d) peneira vibratória.
Normalmente separase a manutenção em três losoas: manutenção
corretiva, manutenção preventiva e manutenção preditiva. Na manutenção
corretiva uma ação só é tomada depois que uma falha acontece. Na
manutenção preventiva algumas peças ou equipamentos, previamente se-
lecionados, são substituídos preventivamente em intervalos de tempos de-
terminados.
Na manutenção preditiva existe um acompanhamento da tendência
de parâmetros, como vibração, desgaste ou temperatura, para que só haja
intervenção em uma máquina caso seja necessário. Através desta losoa,
podese agir próativamente, visto que as tendências do funcionamento do
equipamento estão sendo acompanhadas. Mais do que isso, com ferramentas
de Análise de Sinais e Análise de Vibrações, é possível identicar a fonte do
problema de uma vibração alta.
A Tabela 1.1 mostra uma comparação de custos entre essas diversas
losoas.
17
ABRAMAN∗ NMW Chicago∗∗
Manutenção Corretiva 35 a 40 % 17 a 18 US$ / CV.ano
Manutenção Preventiva 25 a 30 % 11 a 13US$ / CV.ano
Manutenção Preditiva 11 a 17 % 7 a 9 US$ / CV.ano∗ABRAMAN Associação Brasileira de Manutenção, pesquisa realizada entre 1985 a 2003∗∗ NMW National Manufacturing Week, dados de 1998
Tabela 1.1: Comparação de custos entre manutenção corretiva, preventiva
e preditiva
Na Tabela 1.1, o custo percentual (%) está relacionado com o custo
total de manutenção e CV é CavaloVapor. A losoa mais adequada deve
ser analisada para cada situação. O custo para a estruturação de uma
manutenção preditiva é bastante alto.
Sosticando um pouco mais a análise, o próximo passo é fazer um
modelo quantitativo. Com um modelo, e eventuais ajustes, se pode tentar
prever o comportamento de um equipamento para fazer as modicações
necessárias na hora certa.
O presente texto aborda a Análise de Vibrações no contexto da
formulação fraca, que em Mecânica é conhecida como a abordagem em
termos das potências virtuais.
1.2Motivação
Podemse citar duas grandes motivações que permearam o desenrolar
deste trabalho:
1 • Fazer uma apostila didática, usando MATLAB, que seja usada por
alunos de graduação.
2 • Investigar alguns fenômenos observados em máquinas rotativas.
A primeira motivação se deve ao fato de se estar abordando a Análise
de Vibrações no contexto da formulação fraca. Normalmente nos cursos de
vibrações um sistema massamola é analisado e sistemas discretos formados
por massas, molas e amortecedores são discutidos. Nestes casos não se
sabe qual é o erro cometido na análise numérica. A estratégia usada é
formular um problema de vibrações em um espaço de Hilbert e, escolhendo
criteriosamente uma base de projeção, fazer um esquema de aproximação
em dimensão nita de modo a atingir uma precisão pré-estabelecida. O
esquema de convergência para a aproximação é garantido à medida em
18
que se aumenta o número de funções da base usada para representar a
resposta do problema. Esta é a idéia por traz de métodos como o Método
dos Elementos Finitos e o Métodos dos Modos Supostos, que derivam do
Método de Galerkin.
A segunda motivação se deve ao fato da Análise de Vibrações ser muito
usada na manutenção preditiva de máquinas rotativas. Alguns fenômenos
observados nesses equipamentos motivaram o desenvolvimento de um mo-
delo numérico que pudesse reproduzir tais fenômenos para melhor entendê-
los. No capítulo 7 um sistema rotormancal é modelado e sua resposta
dinâmica comparada qualitativamente com a resposta dinâmica captada
através de acelerômetros xados nos mancais de um exaustor da Compan-
hia Siderúrgica de Tubarão (CST).
Em cada etapa do trabalho, pelo menos um exemplo é feito e, para
a maioria deles, um programa foi desenvolvido através da plataforma
MATLAB. A lista dos programas desenvolvidos se encontra no anexo
C. A relação entre os programas se encontra no anexo D. E o manual
dos programas se encontra no anexo E. Adaptando uma frase de Albert
Schweitzer1 digo: Dar um exemplo não é a melhor forma de ensinar, é a
única forma de ensinar, e também de aprender.
1.3Estruturação
No capítulo 2, a partir de sistema mecânico contínuo, os conceitos de
modos de vibrações e freqüências naturais são apresentados. A formulação
fraca surge como uma expressão da dinâmica deste sistema. A teoria é
formulada abstratamente em espaços de Hilbert. Essa forma de trabalhar o
problema facilita o estudo da convergência das aproximações. A equação
originada de sistemas dinâmicos, em geral, não tem como solução uma
expressão analítica, daí a necessidade de discretizar o sistema para buscar
uma aproximação. O método usado para a discretização é o Método de
Galerkin. As matrizes de massa, de amortecimento e de rigidez do sistema
surgem naturalmente de integrações.
No capítulo 3, duas escolhas de funções teste (funções que formam uma
base para a projeção da dinâmica) para o Método de Galerkin: o Método
dos Modos Supostos e o Método dos Elementos Finitos, são discutidas. Os
1Albert Schweitzer (18751965). Example is not the main thing in inuencing others.It is the only thing.
19
modos são as funções teste que geram a melhor base para a projeção da
dinâmica, como será visto.
No capítulo 4, a Análise Modal e suas limitações são examinadas. Além
da Análise Modal só se aplicar à sistemas lineares, dependendo das matrizes
do sistema, ocorrem coisas curiosas como é o caso dos modos complexos.
A Transformada Rápida de Fourier e a Função Resposta em Freqüên-
cia, ferramentas essenciais para a análise de estruturas e equipamentos, são
assuntos do capítulo 5.
As diculdades se tornam maiores em sistemas nãolineares. No
capítulo 6 a Decomposição de Karhunen-Loève é introduzida como uma
extensão à Análise Modal para sistemas dinâmicos não-lineares. A base de
Karhunen-Loève é a melhor base para a projeção da dinâmica, como será
visto.
No capítulo 7 um sistema rotormancal é modelado e sua resposta
dinâmica comparada qualitativamente com a resposta dinâmica captada
através de acelerômetros xados nos mancais de um exaustor da Companhia
Siderúrgica de Tubarão (CST).
2Formulação Fraca
2.1Introdução
Este capítulo apresenta um panorama geral de como um problema de
vibrações será tratado neste trabalho. Os espaços de Hilbert são apresenta-
dos e um problema modelo de uma barra em movimento longitudinal, xa
em uma extremidade e livre na outra, é resolvido. Neste problemamodelo
os conceitos de modos de vibração, freqüências naturais e superposição são
introduzidos. A seguir é feito uma formulação fraca do problema que in-
corpora todos os vínculos do problema. A formulação fraca é o equaciona-
mento do problema de vibrações em um espaço de Hilbert. Depois se usa o
Método de Galerkin para discretizar o problema e reduzílo ao cálculo de
um número nito de parâmetros. Os operadores de um sistema mecânico
são obtidos e uma discussão é feita com relação à dissipação. No nal do
capítulo é mostrado como a convergência das aproximações será medida.
2.2Espaços de Hilbert
Nesta seção serão apresentados os conceitos básicos dos espaços de
Hilbert, porém antes serão apresentados os espaços vetoriais e os espaços
métricos. Esta seção tem como referências básicas Kreyszig [23] e Strang
[42, 41].
Um espaço vetorial é um conjunto X formado por vetores x1,x2, ..
sobre os quais são denidas as operações algébricas de soma vetorial e
multiplicação de vetores por escalares. X é denido sobre um corpo K que
pode ser real, K = R, ou complexo, K = C.
21
A soma é comutativa e associativa:
x1 + x2 = x2 + x1
x1 + (x2 + x3) = (x1 + x2) + x3
(2-1)
Para quaisquer x1,x2 ∈ X e α, β ∈ K:
1. α(βx1) = (αβ)x1
2. 1x = x
3. α(x1 + x2) = αx1 + αx2
4. (α+ β)x1 = αx1 + βx1
(2-2)
Além dessas operações, existe um vetor 0, chamado vetor nulo, e para
qualquer vetor x, existe um vetor (−x), de maneira que:
x+ 0 = x
x+ (−x) = 0(2-3)
Um espaço vetorial familiar é o espaço Euclidiano Rn, cujas operações
são denidas da seguinte forma:
x1 + x2 = (x11 + x21, ..., x1n + x2n)
αx1 = (αx11, ..., αx1n)(2-4)
Outro exemplo de espaço vetorial é o espaço funcional de funções
contínuas em um intervalo [t1, t2] denotado por C[t1, t2]. Este espaço é
formado por funções reais e contínuas denidas no intervalo [t1, t2]. As
operações são denidas por:
(x1 + x2)(t) = x1(t) + x2(t)
(αx1)(t) = αx1(t)(2-5)
Um espaço métrico é um par (X, d), sendo X um conjunto e d uma
métrica, isto é, uma função denida em X x X tal que para quaisquer
x1,x2,x3 ∈ X têmse os seguintes axiomas:
1. d(x1,x2) é real e não negativa
2. d(x1,x2) = 0 ⇐⇒ x1 = x2
3. d(x1,x2) = d(x2,x1)
4. d(x1,x2) ≤ d(x1,x3) + d(x2,x3)
(2-6)
O quarto axioma é conhecido como desigualdade triangular. No espaço
22
Euclidiano tridimensional, R3, uma métrica é dada por:
d(x1,x2) =√
(x11 − x21)2 + (x12 − x22)2 + (x13 − x23)2 (2-7)
O espaço funcional C[t1, t2] também é um espaço métrico, com métrica
denida como:
d(x1,x2) = maxt ∈ [t1,t2]
|x1(t)− x2(t)| (2-8)
Uma εvizinhança de um elemento x0 de um espaço métrico X é um
conjunto aberto:
B(x0, ε) = x ∈ X | d(x,x0) < ε (2-9)
Vizinhança é qualquer subconjunto deX (aberto ou fechado) contendo
uma εvizinhança de x0.
Outro conceito importante é o de espaço métrico completo. Uma
seqüencia xn em um espaço métrico X é dita convergente se existe um
x ∈ X tal que:
limn→∞
d(xn,x) = 0 (2-10)
Sendo x denominado o limite da seqüencia xn. Uma seqüencia xné dita de Cauchy se, para qualquer ε > 0, existe um N(ε) > 0 tal que:
d(xn,x) < ε , para ∀ n > N(ε) (2-11)
Podese provar, Kreyszig [23], que toda seqüência convergente em um
espaço métrico é de Cauchy. Um espaço X é dito completo se toda seqüencia
de Cauchy em X converge.
Um espaço produto interno é um espaço vetorial X com um produto
interno denido em X. Um produto interno é uma aplicação com domínio
X x X e contradomínio no campo escalar K determinado por X, isto é,
para cada par de vetores x1 e x2 há um escalar associado:
〈x1,x2〉 = β , β ∈ K (2-12)
O produto interno satisfaz as seguintes operações:
1. 〈(x1 + x2),x3〉 = 〈x1,x3〉+ 〈x2,x3〉2. 〈αx1,x2〉 = α 〈x1,x2〉3. 〈x1,x2〉 = 〈x2,x1〉4. 〈x1,x1〉 ≥ 0
〈x1,x1〉 = 0 ⇐⇒ x1 = 0
(2-13)
23
A barra denota complexo conjugado. O produto interno dene a
seguinte norma associada:
‖x1‖ =√〈x1,x1〉 (2-14)
E a seguinte métrica:
d(x1,x2) = ‖x1 − x2‖ =√〈(x1 − x2), (x1 − x2)〉 (2-15)
Um espaço de Hilbert H é um espaço produto interno completo. O
espaço Euclidiano Rn é um espaço de Hilbert de dimensão nita.
Um dos conceitos mais importantes que emerge da denição de pro-
duto interno é o de ortogonalidade. Um elemento x1 ∈ X é dito ortogonal
a x2 ∈ X se:
〈x1,x2〉 = 0 (2-16)
A ortogonalidade entre x1 e x2 é representada por x1 ⊥ x2. Um
conjunto ortogonal Q em um espaço produto interno X é um subconjunto
Q ⊂ X cujos elemento são ortogonais aos pares. Os elementos nãonulos
de um conjunto ortogonal são linearmente independentes, ou seja, dados
os elementos x1,x2, ..,xn e os escalares α1, α2, .., αn, a combinação linear
α1x1 + α2x2 + ..αnxn, só será igual a zero se α1 = α2 = .. = αn = 0.
Um conjunto ortonormal P é um conjunto ortogonal cuja norma dos
elementos é unitária:
〈x1,x2〉 =
0 se x1 6= x2
1 se x1 = x2(2-17)
Os conjuntos ortonormais de interesse em espaços de Hilbert são
aqueles que têm um número suciente de elementos, de maneira que todo
elemento do espaço possa ser representado, ou aproximado, através do
uso deste conjunto. Denese um conjunto ortonormal total P o conjunto
capaz de gerar todas as combinações lineares dos elementos do espaço.
Podese provar, Kreyszig [23], que sempre existe um conjunto ortonormal
total em um espaço de Hilbert H 6= 0. Vamos tratar aqui apenas dos
espaços separáveis, ou seja, aqueles que têm uma base enumerável de vetores
ortonormais.
Qualquer elemento x ∈ H pode ser escrito da forma:
x =∞∑
n=1
αnen (2-18)
24
Sendo αn uma seqüência de escalares e en uma base ortonormal
en ⊂ H. en é um subconjunto enumerável de H. A equação (2-18)
também pode ser escrita como:
limn→∞
‖x− (α1e1 + α2e2 + ..+ αnen)‖ = 0 (2-19)
en forma uma base de representação de qualquer elemento de H. A
convergência da expansão é garantida pela relação de Parseval:
∞∑n=1
| 〈x, en〉 |2 = ‖x‖2 (2-20)
A vantagem de se ter uma base é que os coecientes de expansão da
equação (2-18) podem ser determinados fazendose o produto interno do
vetor por cada um dos elementos da base.
〈x, ei〉 =
⟨(∞∑
n=1
αnen
), ei
⟩=
∞∑n=1
αn 〈en, ei〉 = αi (2-21)
Desta forma a equação (2-18) pode ser reescrita como:
x =∞∑
n=1
〈x, en〉 en (2-22)
〈x, en〉 é chamado de coeciente de Fourier de x.
Dadas duas funções f e g contínuas por partes denidas no intervalo
[0, L], um produto escalar pode ser denido por:
〈f, g〉 =
∫ L
0
f(x)g(x)dx (2-23)
A barra denota complexo conjugado. A partir do produto interno
podese obter a seguinte norma:
‖f‖ =
(∫ L
0
|f(x)|2dx)1/2
(2-24)
Sendo f(x)f(x) = |f(x)|2.Um espaço de Hilbert H pode ser representado como uma soma direta
de um subespaço fechado Y e seu complemento ortogonal Z: Y ⊥ = z ∈H |Z ⊥ Y :
X = Y ⊕ Z (2-25)
25
Este resultado também é referido como Teorema da Projeção, Kreyszig
[23]. Explicando de outra forma: estando a solução em um espaço de Hilbert:
x =∞∑
n=1
αnen, ao tentar aproximála com um número determinado de
elementos, N , do espaço:
x =N∑
n=1
αnen︸ ︷︷ ︸xN
+erroN (2-26)
O erro da aproximação é ortogonal ao espaço gerado pelos termos da
aproximação, Figura 2.1.
Solução
Aproximação - XN
erroN
Figura 2.1: Representação do Teorema da Projeção
O Método de Galerkin, que será descrito na seção 2.5, nada mais é do
que o Teorema da Projeção.
2.3Equacionamento
Esta seção tem como referências básicas Rosenberg [33] e Meirovitch
[26]. Das três leis enunciadas por Isaac Newton (1642 1727), a segunda
é a mais utilizada e a mais importante. Antes de Newton, Galileo (1564
1642) já havia desenvolvido os conceitos de aceleração, inércia e referencial
inercial.
Segunda lei de Newton:
Esta lei foi formulada para uma única partícula, mas pode ser esten-
dida para um sistema de partículas e corpos rígidos. O somatório das forças
que atuam em uma partícula é igual à massa da partícula (mi) vezes a
26
aceleração da partícula (ai) em relação a um referencial inercial (xyz). A
equação de movimento de uma partícula mi é dada por:
Fi + fi +N∑
j=1j 6=i
fij = miai (2-27)
Sendo:
N o número de partículas;
fi as forças de vínculo;
fij a força entre as partículas i e j
Fi a força externa.
A equação (2-27) representa a relação entre as forças resultantes e
o movimento de cada partícula de um sistema. A Mecânica Newtoniana
é também conhecida como Mecânica Vetorial por lidar com grandezas
vetoriais como deslocamento e força.
Quando um sistema tiver restrições de movimento (vínculos) apare-
cerão forças de reação (terceira lei de Newton), fi. Dependendo da com-
plexidade do sistema estas forças podem ser complicadas de calcular. O
Princípio de D'Alembert resolve este problema.
Princípio de D'Alembert:
D'Alembert (17171783), de uma maneira engenhosa, estendeu o
Princípio dos Trabalhos Virtuais para sistemas dinâmicos, produzindo assim
o primeiro Princípio Variacional da Dinâmica. Rearrumando a equação (2-
27), Segunda Lei de Newton:
Fi + fi +N∑
j=1j 6=i
fij −miai = 0 (2-28)
Em um sistema em equilíbrio estático, o Princípio dos Trabalhos
Virtuais, originado por Bernoulli, diz que o trabalho realizado pelas forças
de vínculos em um deslocamento virtual, compatível com os vínculos
(restrições) do sistema, é zero, ou seja:
N∑i=1
fi.δri = 0 (2-29)
Sendo δri o deslocamento virtual ou variação admissível. Unindo o
Princípio dos Trabalhos Virtuais da Estática com o Princípio de D'Alembert
27
da Dinâmica temse:
N∑i=1
(Fi −miai).δri =N∑
i=1
fi.δri (2-30)
Segundo o Princípio de D'Alembert o trabalho realizado pelas forças
de vínculos pode ser descartado do problema dinâmico de um sistema de
partículas, ou seja:N∑
i=1
(Fi −miai).δri = 0 (2-31)
O princípio de D'Alembert, assim como a segunda lei de Newton, tem
uma abordagem vetorial e usa coordenadas cartesianas para descrever as
dinâmicas.
Equações de Movimento de Lagrange:
Lagrange (17361813) desenvolveu um tratamento geral de sistemas
dinâmicos formulado por meio de quantidades escalares: Energia cinética
(T ) e trabalho das forças externas (W ). Dependendo da complexidade do
sistema a aplicação direta da segunda lei de Newton se torna impraticável.
Na Mecânica Lagrangeana, ao contrário da Mecânica Newtoniana, o
conceito de coordenadas inclui também coordenadas mais abstratas, não
necessariamente cartesianas. Para descrever um sistema dinâmico com n
graus de liberdade são necessárias pelo menos n coordenadas indepen-
dentes. Tais coordenadas independentes são chamadas de coordenadas gen-
eralizadas.
Em muitos problemas a descrição do movimento por coordenadas
cartesianas xi, yi, zi (i = 1, 2, .., n) não gera um sistema com todas coor-
denadas independentes. Nestes casos é aconselhável descrever o movimento
através de coordenadas generalizadas.
As Equações de Lagrange são de extrema importância na Mecânica
Analítica. Elas representam as equações de movimento de um sistema
dinâmico em termos das coordenadas generalizadas. O Lagrangeano de um
sistema, L, pode ser denido como:
L = T − V (2-32)
No caso de sistemas conservativos a equação de Lagrange é dada por:
d
dt
∂L
∂qi− ∂L
∂qi= 0 (2-33)
28
Sendo qi a iésima coordenada generalizada. Para sistemas não con-
servativos, a equação de Lagrange é dada por:
d
dt
∂L
∂qi− ∂L
∂qi= Qi (2-34)
SendoQi a força não conservativa associada à coordenada generalizada
qi.
2.4Problema Modelo
Um problema modelo é analisado nesta seção. O objetivo é introduzir
os conceitos de modos de vibração, freqüências naturais e superposição.
Considere a barra representada na Figura 2.2:
A
x=0 x=L
f(x,t)
Figura 2.2: Barra xa em um extremidade e livre na outra
A equação diferencial que descreve o movimento longitudinal de uma
barra pode ser obtida aplicandose a segunda lei de Newton em um elemento
diferencial, Figura 2.3.
i i i
Figura 2.3: Elemento diferencial de uma barra
ρ(x)A(x)dx∂2u(x, t)
∂t2=
(Pi(x, t) +
∂Pi(x, t)
∂xdx
)−Pi(x, t)+f(x, t) (2-35)
Sendo
u(x, t) o deslocamento longitudinal;
Pi(x, t) a força longitudinal interna;
29
f(x, t) a força longitudinal externa;
A(x) a área da seção transversal;
ρ(x) a massa por unidade de comprimento.
O modelo de barra leva em consideração apenas o deslocamento na
direção longitudinal da barra.
Pela lei de Hooke:Pi(x, t)
A(x)= E(x)ε(x, t)
Sendo E(x) o módulo de elasticidade do material e ε(x, t) =∂u
∂xa
deformação, logo:
ρ(x)A(x)∂2u(x, t)
∂t2=
∂
∂x
(A(x)E(x)
∂u(x, t)
∂x
)+ f(x, t) (2-36)
Considerando a seção transversal da barra (A) e as propriedades do
material (E e ρ) constantes em x, e considerando f(x, t) = 0, podese
escrever:∂2u(x, t)
∂t2= c2
∂2u(x, t)
∂x2(2-37)
Sendo c2 = E/ρ. A equação (2-37) é a equação de movimento da barra.
As condições de contorno são dadas por:
u(0, t) = 0 EA∂u
∂x(L, t) = 0 (2-38)
As condições iniciais são dadas por:
u(x, 0) = u0(x)∂u
∂t(x, 0) = v0(x) (2-39)
A primeira condição de contorno, u(0, t) = 0, é uma condição de
contorno essencial (ou de Dirichlet, ou cinemática, ou geométrica), [41].
Esta condição deve ser imposta. A barra está xa na extremidade x = 0,
ou seja, este ponto não se move.
A segunda condição de contorno,∂u
∂x(L, t) = 0, é uma condição de
contorno natural (ou de Neumann, ou dinâmica), [41]. A barra está livre na
extremidade x = L, ou seja, a força de reação neste ponto é zero.
As condições de contorno serão tratadas de forma distinta, como será
visto na seção 2.5. Uma delas, u(0, t) = 0, será incorporada na denição do
espaço.
O problema (2-37, 2-38 e 2-39) será resolvido em três passos:
30
1. O Método de Separação de Variáveis é usado para obter, em princípio,
duas EDOs.
2. As soluções para as duas Equações Diferenciais Ordinárias (EDOs)
que satisfazem as condições de contorno são determinadas.
3. A solução do problema é composta através da superposição de
soluções. Esta etapa é necessária porque, individualmente, as pro-
postas não satisfazem as condições iniciais.
1. Primeiro passo Método de Separação de Variáveis.
Solução proposta da forma:
u(x, t) = X(x)T (t) (2-40)
A variável u que depende de x e t é separada no produto de duas
funções: X que só depende de x e T que só depende de t. Substituindo a
equação (2-40) na equação (2-37) obtémse:
XT = c2X ′′T ∴X ′′
X=
T
c2T= −λ2 (2-41)
Sendo T =dT
dte X ′ =
dX
dx
A equação inicial, formulada como uma Equação Diferencial Parcial
(EDP), foi separada em duas EDOs:
T + λ2c2T = 0 (2-42)
X ′′ + λ2X = 0 (2-43)
2. Segundo passo Determinação das soluções das equações (2-42) e
(2-43) que satisfaçam as condições de contorno do problema:
u(0, t) = X(0)T (t) = 0 ∴ X(0) = 0
EA∂u(L, t)
∂x= EA
dX(L)
dxT (t) = 0 ∴
dX
dx
∣∣∣∣x=L
= 0
(2-44)
Esse problema: Achar os pares (λ,X), para X 6= 0, tais que X ′′ +
λ2X = 0, X(0) = 0 e X ′(0) = 0; é um problema de autovalor. (X(x) pode
31
se anular em alguns pontos do intervalo [0, L]). A solução geral da equação
(2-43) é dada por:
X(x) = C1cos(λx) + C2sen(λx) (2-45)
Sendo C1 e C2 coecientes a serem determinados a partir das condições
de contorno. Substituindo a primeira condição de contorno, X(0) = 0, na
equação (2-45) obtémse C1 = 0, logo:
X(x) = C2sen(λx) (2-46)
Pela segunda condição de contorno:
EAdX
dx
∣∣∣∣x=L
= (EA)λC2cos(λL) = 0 (2-47)
Como buscase uma solução não trivial, necessariamente, C2 6= 0 e
λ 6= 0. Neste ponto surge uma restrição que não fôra imposta anteriormente:
λ deve satisfazer a seguinte equação,
cos(λL) = 0 ∴ λnL =(2n− 1)π
2, sendo n = 1, 2, 3, .. (2-48)
Procuravase uma solução, mas neste ponto notase que existem
innitas soluções para o problema, pois n = 1, 2, 3, ...
λn =(2n− 1)π
2L(2-49)
As soluções da equação (2-43) que satisfazem as condições de contorno
(2-38) são:
Xn(x) = C2nsen (λnx) n = 1, 2, 3, ... (2-50)
O resultado contempla os λn possíveis, tal que X 6= 0. A função
nula não interessa pois procurase uma solução que permita satisfazer as
condições iniciais.
Reescrevendo as equações (2-42) e (2-43):
Tn + λ2nc
2Tn = 0 (2-51)
X ′′n + λ2
nXn = 0 (2-52)
32
A solução geral da equação (2-51) é dada por:
Tn(t) = C3ncos(λnct) + C4nsen(λnct) (2-53)
As soluções são dadas por: un(x, t) = Xn(x)Tn(t). Multiplicando as
equações (2-50) e (2-53):
un(x, t) = sen(λnx)(C5ncos(cλnt) + C6nsen(cλnt)) (2-54)
Sendo C5n = C2nC3n e C6n = C2nC4n coecientes a serem determina-
dos a partir das condições iniciais.
Essas equações, 2-54, não satisfazem, necessariamente, as condições
iniciais. Para contornar este fato, surge a id'eia de somar todas as soluções,
ou seja, fazer uma superposição das soluções.
3. Terceiro Passo Superposição.
A solução geral da equação de movimento da barra (2-37) que satisfaz
as condições de contorno (2-38) para qualquer condição inicial é dada pela
superposição das soluções:
un(x, t) =∞∑
n=1
sen(λnx) (C5ncos(cλnt) + C6nsen(cλnt)) (2-55)
Sendo os coecientes C5n e C6n determinados a partir das condições
iniciais. Considere a primeira condição inicial u(x, 0) = u0(x). Substituindo
a na equação (2-55) chegase a:
u0(x) =∞∑
n=1
C5nsen(λnx) (2-56)
Multiplicando a equação (2-56) por sen(λmx) e integrando no intervalo
[0, L], obtémse:
∫ L
0
u0(x)sen(λmx)dx =
∫ L
0
∞∑n=1
C5nsen(λnx)sen(λmx)dx (2-57)
Pela propriedade de ortogonalidade da função seno, podese escrever:
33
∫ L
0
sen
((2n− 1)π
2Lx
)sen
((2m− 1)π
2Lx
)dx =
L/2 para n = m
0 para n 6= m
(2-58)
Todos os termos dentro dos somatórios da equação (2-57) serão zero,
exceto o termo n = m. Logo o coeciente C5n pode ser escrito como:
C5n =2
L
∫ L
0
u0(x)sen(λnx)dx (2-59)
Agora considere a condição inicial u(x, 0) = v0(x). Substituindoa na
equação (2-55) chegase a:
v0(x) =∞∑
n=1
C6ncλnsen(λnx) (2-60)
Multiplicando a equação (2-60) por sen(λmx) e integrando no intervalo
[0, L], obtémse:
∫ L
0
v0(x)sen(λmx)dx =
∫ L
0
∞∑n=1
C6ncλnsen(λnx)sen(λmx)dx (2-61)
Usando novamente as propriedades de ortogonalidade chegase ao
coeciente C6n:
C6n =2
Lcλn
∫ L
0
v0(x)sen(λnx)dx (2-62)
Portando a solução geral da equação de movimento de uma barra
xa em uma extremidade e livre na outra é dada pela equação (2-55) com
coecientes estabelecidos pelas expressões (2-59) e (2-62).
Reescrevendo a equação (2-55) de uma forma compacta:
u(x, t) =∞∑
n=1
an(t)φn(x) (2-63)
Sendo an(t) = (C5ncos(cλnt) + C6nsen(cλnt)) coecientes temporais
que dependem das condições iniciais, e φn(x) = sen(λnx) os modos de
vibração do sistema.
Modos de vibração e freqüências naturais:
Modo de vibração é uma conguração na qual todos os pontos de um
sistema vibram em uma freqüência única (freqüência natural). No problema
34
analisado os modos de vibração são dados por:
φn(x) = sen
((2n− 1)πx
2L
)= sen (λnx) (2-64)
Freqüências naturais são chamadas de freqüências de ressonância.
Se o sistema é excitado em uma destas freqüências ele tem sua resposta
amplicada. As freqüências naturais são dadas pela expressão:
ωn =(2n− 1)πc
2L= cλn (2-65)
Os modos de vibração e as freqüências naturais carregam informações
fundamentais de sistemas dinâmicos. No capítulo 4 os modos e as freqüên-
cias naturais são obtidos para um sistema dinâmico discreto.
Modo de corpo rígido:
Sistemas nãoxados, apresentam movimento de corpo rígido, isto é,
um deslocamento do corpo sem que haja deformação elástica. Por exemplo,
uma barra livre nas duas extremidades, cujas condições de contorno são
dadas por:
EA∂u
∂x(0, t) = 0 EA
∂u
∂x(L, t) = 0 ∴ EA
dX
dx
∣∣∣∣x=0,L
= 0 (2-66)
A derivada da equação (2-45) é dada por:
X ′(x) = −C1λsen(x) + C2λcos(x) (2-67)
Substituindo a primeira condição de contorno, X ′(0) = 0, na equação
(2-67):
X ′(0) = −C1λsen(0)︸ ︷︷ ︸=0
+C2λ cos(0)︸ ︷︷ ︸=1
= 0 ∴ C2 = 0 (2-68)
Substituindo a segunda condição de contorno, X ′(L) = 0, na equação
(2-67):
X ′(L) = −C1λsen(λL) = 0 (2-69)
O valor de λ pode ser zero, pois nesse caso X|λ=0 6= 0. O sistema é
positivo semidenido e devido a este fato é possível haver autovalor igual
a zero, ou seja, freqüência natural igual a zero. Os modos relacionados à
freqüência zero são denominados modos de corpo rígido. A equação (2-43):
35
X ′′ + λ2X = 0, para λ = λ0 = 0 e X = X0, pode ser escrita como:
d2X0(x)
dx2= 0 (2-70)
Integrando duas vezes a equação (2-70) obtémse:
X0(x) = A0 (2-71)
Sendo A0 uma constante de integração. Normalizando o modo de
forma que∫ L
0ρX2
0dx = 1:
X0(x) =1√ρL
(2-72)
X0 representa um modo de corpo rígido correspondente à freqüência
natural, ω0 = 0, Figura 2.4. Estes modos ocorrem em sistemas sem restrição
onde não há forças exercidos por suportes.
Sistema Discreto
Resposta
PrecisãoSatisfeita?
Aumentar número de elementosda base usados na aproximação
ok!
Não
Sim
x=0
x=L
Lr
1
Figura 2.4: Modo de corpo rígido de uma barra com as duas extremidades
livres
Aproximação:
Os innitos modos de vibrações do sistema, φn, são funções do espaço
C[0, L]. Os modos de vibração geram um espaço de Hilbert, que é o
completamento de C[0, L].
A solução do problema de vibração longitudinal da barra está neste
espaço e ela é aproximada quando se escolhe um número nito de modos
para representar a solução:
uN(·, t) =N∑
n=1
an(t)φn(x) (2-73)
36
Sendo uN(·, t) a aproximação de u(·, t) com N modos de vibração.
limN→∞
‖uN(·, t)− u(·, t)‖ = 0 (2-74)
A convergência é garantida pelas propriedades dos espaços de Hilbert,
introduzidas na seção 2.2. Dependendo da dinâmica em análise serão
necessários mais ou menos modos para satisfazer uma determinada precisão.
O erro da aproximação será a soma dos termos que não forem considerados
no cálculo.
erroN(·, t) =∞∑
n=N+1
un(·, t) (2-75)
Sendo erroN o erro da resposta dinâmica para n = N . Estabelecida
uma precisão, é possível calcular o número de termos necessários para se
obter uma aproximação dentro da precisão requerida.
Exemplo 2.1: Considere uma barra de aço (E=200 MPa e
ρ=7850kg/m3) com um metro de comprimento, xa em uma extremidade e
livre na outra. Os modos e freqüências naturais podem ser calculados pelas
seguintes expressões, encontradas em [18] e [3]:
Os modos de vibração:
φn(x) = sen(λnx) (2-76)
Sendo n = 1, 2, 3, .. e:
λn =(2n− 1)π
2L(2-77)
As freqüências naturais:
ωn = cλn (2-78)
A Figura 2.5 mostra os dez primeiros modos de vibração calculados.
37
Figura 2.5: Dez primeiros modos de vibração de uma barra xa em uma
extremidade e livre na outra
As dez primeiras freqüências naturais calculadas (rad/s) são:
7.930
23.790
39.640
55.500
71.360
87.220
103.070
118.930
134.790
150.640
Essas aproximações foram obtidas utilizandose o programa
modal_barra (ver Anexo E). Este programas também calcula os mo-
dos de vibração e as freqüências naturais de uma barra xa nas duas
extremidade pelas seguintes expressões encontradas em [18] e [3]:
38
Modos de vibração:
φn(x) = sen(λnx) (2-79)
Sendo:
λn =nπ
L(2-80)
Freqüências naturais:
ωn = cλn (2-81)
A Figura 2.6 mostra os dez primeiros modos de vibração calculados.
Figura 2.6: Dez primeiros modos de vibração de uma barra xa nas duas
extremidades
As dez primeiras freqüências naturais calculadas (rad/s) são:
39
15.860
31.710
47.570
63.430
79.290
95.140
111.000
126.860
142.720
158.570
2.5Formulação Fraca
As referências básicas usadas nesta seção foram Strang [41] e Hughes
[17]. Para introduzir o conceito de formulação fraca será usado um exemplo
unidimensional de um problema de valor de contorno, equação (2-82).
− d
dx
(p(x)
du
dx
)+ q(x)u = f(x) , 0 ≤ x ≤ 1 (2-82)
Sendo as funções u, p, q e f denidas no intervalo [0, 1]. As condições
de contorno do problema são dadas por:
u(0) = 0du
dx(1) = 0 (2-83)
Este é um problema pertencente à classe de problemas de Sturm-Liouville.
A forma mais simples deste problema é, para p(x) = 1 e q(x) = 0:
d2u
dx2+ f(x) = 0 (2-84)
Será suposto que f tenha energia nita, ou seja:∫ 1
0
(f(x))2 dx <∞ (2-85)
A expressão (2-85) diz que f pode ser qualquer função contínua por
partes.
O espaço de funções que satisfaz a equação (2-85) é notado por H0
(=L2(0, L)), sendo o índice sobrescrito um indicativo do número requerido
de derivadas com energia nita para a função f . Neste caso, f não precisa
ter derivada com energia nita.
40
Uma função f que pertence ao espaço de funções denotado por H1C ,
signica que a primeira derivada de f deve ter energia nita, ou seja:∫ 1
0
(df
dx
)2
dx <∞ (2-86)
O índice subscrito, C, se refere às condições de contorno u(0) = 0 e
u′(L) = 0.
O problema da equação (2-84) com condições de contorno (2-83)
como está denido até agora é conhecido como formulação forte. Alguns
métodos de aproximação usam a formulação forte do problema. O mais
conhecido e notável é o Método das Diferenças Finitas. O MEF requer uma
formulação diferente, a formulação fraca.
Partindo de uma formulação forte, para chegar à formulação fraca
são necessários 2 passos:
1. Produto interno da equação (2-84) com w.
⟨d2u
dx2, w
⟩+ 〈f, w〉 = 0 (2-87)
Sendo w uma função teste. A equação (2-87) também pode ser escrita
como: ∫ 1
0
d2u(x)
dx2w(x)dx+
∫ 1
0
f(x)w(x)dx = 0 (2-88)
2. Integração por partes da equação (2-88). Lembrando:∫ 1
0udv =
uv|10 −∫ 1
0vdu:
∫ 1
0
du
dx
dw
dxdx− w
du
dx
∣∣∣∣10
=
∫ 1
0
fwdx (2-89)
Observe o termo em destaque da equação (2-89). As condições de
contorno são incorporadas na equação através deste termo. Anteriormente
foram denidas duas condições de contorno, a natural e a essencial, agora
será explicada a importâncias dessas denições.
A função teste, w, deve satisfazer a condição de contorno essencial,
u(0) = 0, mas não precisa satisfazer a condição de contorno natural,
u′(1) = 0. Isto dene o espaço de funções admissíveis.
Outro fato vantajoso da formulação fraca é que a exigência quanto às
derivadas da função u é de derivada de primeira ordem, enquanto que na
formulação forte u precisava ter derivada de segunda ordem.
41
Portanto, na formulação fraca, as funções u e w pertencem ao espaço
denotado por H1E. O índice subscrito E se refere à condição de contorno
essencial.
De uma forma geral, se na formulação fraca a exigência é de derivada
de ordem N , as condições de contorno essenciais aparecem com derivadas,
no máximo, de ordem (N − 1). Comparativamente, na formulação forte
a exigência seria de derivadas de ordem 2N . Resumindo, o problema na
formulação fraca é denido como:
Determinar u ∈ H1E, tal que:∫ 1
0
du
dx
dw
dxdx− w(1)
du
dx(1)︸ ︷︷ ︸
=0
+w(0)︸︷︷︸=0
du
dx(0) =
∫ 1
0
fwdx , ∀w ∈ H1E
(2-90)
Ou, ∫ 1
0
du
dx
dw
dxdx =
∫ 1
0
fwdx , ∀w ∈ H1E(0, L)
w ∈ H1(0, L) : w(0) = 0
(2-91)
Pelas condições de contorno naturais e essenciais, os termos que
carregam as condições de contorno da equação são zero.
Através da formulação fraca, é possível aproximar a solução de um
problema através de funções mais simples, com exigência menor quanto ao
número de derivadas. H1E ⊂ H2
C , observe a Figura 2.7.
Espaço de funções para formulação fraca
Espaço de funções para formulação forte
Figura 2.7: Espaço de soluções
Dado um problema na formulação fraca, caso não haja como obter
uma solução em forma analítica, devese aproximar u. Será considerada
42
uma aproximação da forma:
uN(x) =N∑
i=1
aiφi(x) (2-92)
Sendo ai coecientes a serem determinados a partir de condições su-
plementares e φi funções linearmente independentes, previamente denidas.
N é a dimensão do espaço gerado por φ1, φ2, .., φN . Este espaço de funções
está contido no espaço H1E.
Esta aproximação gera um erro:
u(x) = uN(x) +∞∑
i=N+1
aiφi(x)︸ ︷︷ ︸erroN
(2-93)
Reescrevendo a equação (2-90) com a aproximação:∫ 1
0
duN(x)
dx
dw(x)
dxdx−
∫ 1
0
f(x)w(x)dx =
∫ 1
0
d(erroN(x))
dx
dw(x)
dxdx (2-94)
A função teste, w, será aproximada por ψ1, ψ1, .., ψN . Assim como
φi, estas funções são linearmente independentes e geram um espaço de
funções que está contido no espaço H1E. Substituindo a aproximação (2-92)
na equação (2-94) chega-se a:
ai
∫ 1
0
dφi(x)
dx
dψi(x)
dxdx−
∫ 1
0
fψi(x)dx =
∫ 1
0
d(erroN(x))
dx
dψi(x)
dxdx (2-95)
No Método dos Resíduos Ponderados, [26], o erro deve ser ortogonal
ao espaço formado pelas funções teste usadas para representar o problema.
Esta funções teste formam uma base para a projeção da dinâmica. Resíduos
Ponderados é um nome genérico dado a uma família de métodos que utiliza
o Teorema da Projeção, descrito na seção 2.2.
43
Espaço das funções teste
Erro
Figura 2.8: Erro ortogonal ao espaço formado pelas funções teste usadas na
aproximação
Isto signica que:∫ 1
0
d(erroN(x))
dx
dψi(x)
dxdx = 0 (2-96)
A equação (2-95) ca:
ai
∫ 1
0
dφi(x)
dx
dψi(x)
dxdx−
∫ 1
0
fψi(x)dx = 0 (2-97)
Os diversos métodos se distinguem na escolha das funções teste, ψi e
das funções aproximantes φi.
O método de BubnovGalerkin é conhecido como o método simples de
Galerkin. Neste método o espaço das funções teste e funções aproximantes
é idêntico:
ψi = φi (2-98)
A equação (2-99) mostra a formulação resultante após a aplicação do
Método de Galerkin no problema da equação (2-84):
ai
∫ 1
0
dφi(x)
dx
dφi(x)
dxdx−
∫ 1
0
fφi(x)dx = 0 (2-99)
2.6Operadores Diferenciais
44
2.6.1Operador Autoadjunto
Um operador L é dito autoadjunto se:
〈Lu,w〉 = 〈u, Lw〉 (2-100)
Sendo u e w funções teste que satisfazem todas as condições de
contorno. Se L é simétrico L = LT .
2.6.2Sistema Conservativo
Considere o problema da barra xa em uma extremidade e livre na
outra discutido na seção 2.4. Aplicando o Método de Galerkin na equação
de movimento de uma barra:
ρA∂2u(x, t)
∂t2− EA
∂2u(x, t)
∂x2= f(x, t) (2-101)
Primeiro a equação de movimento da barra é projetada na base
formada pelas funções teste:∫ L
0
[ρA
∂2u(x, t)
∂t2− EA
∂2u(x, t)
∂x2− f(x, t)
]φj(x)dx = 0 (2-102)
As funções teste, φi, geram um espaço de Hilbert H1[0, L].
Integrando por partes a equação (2-102):
ρA
∫ L
0
∂2u(x, t)
∂t2φj(x)dx+ EA
∫ L
0
∂u(x, t)
∂x
dφj(x)
dxdx+
−EA ∂u(x, t)
∂xφj(x)
∣∣∣∣L0
=
∫ L
0
f(x, t)φj(x)dx
(2-103)
Desenvolvendo o termo em destaque na equação (), chegase a:
−EA ∂u(x, t)
∂xφj(x)
∣∣∣∣L0
= −EA ∂u
∂x(x = L)︸ ︷︷ ︸
=0
φj(L) + EA∂u
∂x(x = 0)φj(0)︸ ︷︷ ︸
=0
= 0
(2-104)
Devido à condição de contorno natural,∂u
∂x(L) = 0. E devido à
condição de contorno essencial, φj(0) = 0. Fazendo uma aproximação na
45
forma:
uN(x, t) =N∑
i=1
ai(t)φi(x) (2-105)
Substituindo a aproximação (2-105) na equação (2.6.2), já incorpo-
rando as condições de contorno e considerando o erro ortogonal aos termos
da aproximação, chegase a:
ai(t)ρA
∫ L
0
φi(x)φj(x)dx+ ai(t)EA
∫ L
0
dφi(x)
dx
dφj(x)
dxdx =
=
∫ L
0
fi(x, t)φj(x)dx
(2-106)
Ou:
M(φi, φj)ai(t) +K(φi, φj)ai(t) = F (φj) (2-107)
O operador bilinear de massa ou inercial, M(φi, φj), o operador
bilinear de deformação ou elástico, K(φi, φj) e o operador linear das forças
externas, F (φj), são denidos da seguinte forma:
M(φi, φj) = ρA
∫ L
0
φi(x)φj(x)dx K(φi, φj) = EA
∫ L
0
dφi(x)
dx
dφj(x)
dxdx
F (φj) =
∫ L
0
f(x, t)φj(x)dx
(2-108)
Com a discretização, os operadores, M(φi, φj) e K(φi, φj) se tornam
matrizes e F (φj) um vetor:
Mjiai(t) +Kjiai(t) = Fi(t) ∴
M a(t) +Ka(t) = F(t)
(2-109)
Sendo i = 1, 2, .., N e j = 1, 2, .., N . Fazendo o paralelo com sistemas
discretos, sistemas auto-adjuntos são aqueles que possuem matrizes de
massa (M) e de rigidez (K) simétricas, [26, 34]. Isto implica no sistema
gerar autovalores e autovetores reais, propriedade muito importante que
será explorada mais adiante.
2.6.3Sistema com Dissipação
Modelar dissipação é um grande desao e não existe, até hoje, uma
teoria satisfatória. Geralmente, em uma abordagem inicial, fazse um mo-
delo supondo que o sistema é conservativo e a dissipação é introduzida
46
a posteriori. Essa é a abordagem que será usada. Será considerado que a
dissipação é viscosa, o caso mais simples, pois leva a um modelo linear.
Existe uma diculdade grande para conseguir identicar o amortecimento
de uma estrutura. Em muitos casos aproximações são feitas considerando o
amortecimento proporcional ou de Rayleigh (1877):
Cp = αM + βK (2-110)
No capítulo 4 serão comentadas outras vantagens dessa aproximação.
Trabalhos como o de Adhikari [1] e Phani [29] lidam com a mode-
lagem de dissipações e maneiras de identicar o amortecimento tanto ex-
perimentalmente quanto numericamente. Nesses trabalhos modelos mais
complicados de amortecimento são propostos, como amortecimentos não
proporcionais e amortecimentos locais.
Nesta dissertação só será visto o amortecimento viscoso. No capítulo
4 será considerado um caso de amortecimento nãoproporcional e suas
implicações para o sistema dinâmico.
Inserindo uma dissipação na equação de uma barra:
ρA∂2u(x, t)
∂t2+ c
∂u(x, t)
∂t− EA
∂2u(x, t)
∂x2= f(x, t) (2-111)
Sendo c o coeciente de amortecimento.
Aplicando o Método de Galerkin, primeiro a equação de movimento
da barra é projetada na base formada pelas funções teste:∫ L
0
[ρA
∂2u(x, t)
∂t2+ c
∂u(x, t)
∂t− EA
∂2u(x, t)
∂x2− f(x, t)
]φj(x)dx = 0
(2-112)
As funções teste, φi, geram um espaço de Hilbert H1[0, L].
Integrando por partes a equação (2-112):
ρA
∫ L
0
∂2u(x, t)
∂t2φj(x)dx+ EA
∫ L
0
∂u(x, t)
∂x
dφj(x)
dxdx
+c
∫ L
0
∂u(x, t)
∂tφj(x)dx −EA
∂u(x, t)
∂xφj(x)
∣∣∣∣L0
=
∫ L
0
f(x, t)φj(x)dx
(2-113)
47
Desenvolvendo o termo em destaque na equação (2-113), chegase a:
−EA ∂u(x, t)
∂xφj(x)
∣∣∣∣L0
= −EA ∂u
∂x(x = L)︸ ︷︷ ︸
=0
φj(L) + EA∂u
∂x(x = 0)φj(0)︸ ︷︷ ︸
=0
= 0
(2-114)
Devido à condição de contorno natural,∂u
∂x(L) = 0. E devido à
condição de contorno essencial, φj(0) = 0.
Fazendo uma aproximação na forma:
uN(x, t) =N∑
i=1
ai(t)φi(x) (2-115)
Substituindo a aproximação (2-115) na equação (2-113), já incorpo-
rando as condições de contorno e considerando o erro ortogonal aos termos
da aproximação, chegase a:
M(φi, φj)ai(t) + C(φi, φj)ai(t) +K(φi, φj)ai(t) = F (φj) (2-116)
O operador de massa ou inercial,M(φi, φj), o operador de deformação
ou elástico, K(φi, φj), o operador de dissipação, C(φj), e o operador das
forças externas, F (φj), são denidos da seguinte forma:
M(φi, φj) = ρA
∫ L
0
φi(x)φj(x)dx , K(φi, φj) = EA
∫ L
0
dφi(x)
dx
dφj(x)
dxdx
F (φj) =
∫ L
0
f(x, t)φj(x)dx , C(φi, φj) = c
∫ L
0
φi(x)φj(x)dx
(2-117)
Repare que em (2-117) temse que: C(φi, φj) = αM(φi, φj), ou seja, o
modelo de dissipação considerado é o de amortecimento proporcional. Com
a discretização o sistema matricial ca:
Mjiai(t) + Cjiai(t) +Kjiai(t) = Fi(t) ∴
M a(t) + Ca(t) +Ka(t) = F(t)
(2-118)
Sendo i = 1, 2, .., N e j = 1, 2, .., N .
48
2.7Convergência das Aproximações
Usando a Teorema da Projeção e fazendo as aproximações como
foram demonstradas nas seções anteriores, garantese que a aproximação da
solução é representada por elementos de um espaço de Hilbert. Quanto mais
elementos forem incorporados na aproximação, mais próxima a aproximação
está da solução.
Podese querer mais ou menos rigor nas aproximações, porém, para
quanticar a precisão de uma aproximação devese calcular o erro.
Uma aproximação com cinco modos pode ser boa ou ruim dependendo
da precisão que se exige para a aproximação. Para maiores detalhes sobre
cálculo de erro consultar [8], [6], [22] e [17].
Considere que a solução de um problema é dada pelo escalar sexa, e
que a aproximação desta solução é dada por outro escalar saprox. O erro
absoluto é dado por:
eabs = |sexa − saprox| (2-119)
Esse erro absoluto depende do valor da solução, ou seja, ele varia se a
solução tem um valor maior ou menor.
O erro relativo é dado por:
erel =
∣∣∣∣sexa − saprox
sexa
∣∣∣∣ (2-120)
Para melhor analisar as aproximações, usase o erro percentual:
eperc =
∣∣∣∣sexa − saprox
sexa
∣∣∣∣ 100 (2-121)
Esse é o erro que será usado neste trabalho.
Erro das aproximações do coeciente de Rayleigh:
O erro da aproximação de um modo de vibração será computado pelo
coeciente de Rayleigh Ritz. Considere as energias cinética e potencial de
cada aproximação do modo de vibração:
Ti = 12ω2
i φTi Mφi Energia cinética do iésimo modo.
49
Vi = 12φT
i Kφi Energia potencial do iésimo modo.
Sendo ωi a iésima freqüência natural e K e M as matrizes de rigidez
e de massa do sistema, respectivamente. φi é o iésimo modo de vibração.
Observe que não é toda a massa do sistema que vai contribuir para a
energia cinética de determinado modo. Cada modo tem uma massa, um
amortecimento e uma rigidez que lhe é peculiar, isto é, dependem da forma
do modo. O coeciente de Rayleigh é dado por:
λ(φi) = ω2i =
φTi Kφi
φTi Mφi
=Vi
Ti
(2-122)
O erro percentual do iésimo modo de vibração é dado por:
eperc =
∣∣∣∣∣λ(N)i − λ
(N−1)i
λ(N)i
∣∣∣∣∣ 100 (2-123)
Sendo N o número de elementos da base usado na aproximação.
O erro é maior para os modos de maior ordem. Para ver maiores
detalhes, consulte [26, 34].
Erro das aproximações da resposta dinâmica de um ponto:
Para computar o erro da aproximação de uma resposta dinâmica de um
ponto será usado o conceito de norma. Considere uma função f de dimensão
n. A partir do produto interno podese obter a seguinte norma
‖f‖ =
(∫ L
0
|f(x)|2dx)1/2
(2-124)
A medida usada para computar o erro entre duas funções fN e f(N−1)
será o erro percentual da norma:
eperc =
∣∣∣∣‖f(N)‖ − ‖f(N−1)‖‖f(N)‖
∣∣∣∣ 100 (2-125)
Sendo N o número de elementos da base usados na aproximação.
O algoritmo para a convergência é mostrado na Figura 2.9:
50
Sistema Discreto
Resposta
PrecisãoSatisfeita?
Aumentar número de elementosda base usados na aproximação
ok!
Não
Sim
Figura 2.9: Algoritmo para convergência
2.8Conclusões
Este é o capítulo que contém os fundamentos básicos do trabalho.
Todos os outros capítulos se apóiam nas teorias aqui descritas. No próximo
capítulo dois métodos para a escolha da base de projeção da dinâmica são
discutidos.
3Escolha da Base de Projeção
3.1Introdução
A escolha da base de projeção para a dinâmica, formada pelas funções
teste do Método de Galerkin, é feita neste capítulo. Dois métodos são
discutidos: Método dos Modos Supostos e Método dos Elementos Finitos.
3.2Método dos Modos Supostos
Será considerado um problema de uma viga para exemplicar a
aplicação do Método dos Modos Supostos. Considere o elemento de viga
representado na Figura 3.1:
dw/dx
w
Figura 3.1: Viga de EulerBernoulli
O modelo de viga EulerBernoulli, equação (3-1), considera que as
seções retas da viga permanecem planas e perpendiculares à linha central
52
da viga, [26, 34].
ρ(x)A(x)∂2w(x, t)
∂t2+
∂2
∂x2
(E(x)I(x)
∂2w(x, t)
∂x2
)− ∂
∂x
(P (x, t)
∂w(x, t)
∂x
)=
= f(x, t)
(3-1)
Sendo:
w(x, t) o deslocamento transversal;
A(x) a área da seção transversal;
I(x) o momento de inércia da seção transversal;
E(x) o módulo de elasticidade (Young);
ρ(x) a massa por unidade de comprimento;
f(x, t) a força (transversal) por unidade de comprimento;
P (x, t) a força axial.
Considerando P (x, t) = 0 e A(x), I(x), E(x) e ρ(x) constantes:
ρA∂2w(x, t)
∂t2+
∂2
∂x2
(EI
∂2w(x, t)
∂x2
)= f(x, t) (3-2)
Será aplicado o Método de Galerkin. O primeiro passo do Método de
Galerkin é a projeção da equação de movimento, equação (3-2), na base
formada pelas funções teste:∫ L
0
[ρA
∂2w(x, t)
∂t2+
∂2
∂x2
(EI
∂2w(x, t)
∂x2
)− f(x, t)
]φj(x)dx = (3-3)
Integrando por partes a equação (3-3):∫ L
0
ρA∂2w(x, t)
∂t2φj(x)dx−
∫ L
0
dφj(x)
dx
∂
∂x
(EI
∂2w(x, t)
∂x2
)+
+ φj(x)∂
∂x
(EI
∂2w(x, t)
∂x2
)∣∣∣∣L0
=
∫ L
0
f(x, t)φj(x)dx
(3-4)
Integrando por partes novamente:∫ L
0
ρA∂2w(x, t)
∂t2φj(x)dx+
∫ L
0
EId2φj(x)
dx2
∂2w(x, t)
∂x2+
− EIdφj(x)
dx
∂2w(x, t)
∂x2
∣∣∣∣L0
+ φj(x)∂
∂x
(EI
∂2w(x, t)
∂x2
)∣∣∣∣L0
=
∫ L
0
f(x, t)φj(x)dx
(3-5)
53
Considere a viga representada na Figura 3.2:
Figura 3.2: Viga engastadalivre
As condições de contorno de uma viga engastada em uma extremidade
e livre na outra são:
Em x = 0 estão as condições de contorno essenciais:
w(0, t) = 0 ∴dφi
dx(0) = 0
∂w
∂x(0, t) = 0 ∴ φi(0) = 0
(3-6)
Estas condições de contorno devem ser respeitadas pelas funções teste,
como foi visto no capítulo 2. Em x = L estão as condições de contorno
naturais:
Q(L, t) = − ∂
∂x
(EI
∂2w
∂x2
)∣∣∣∣x=L
= 0 ∴∂2w
∂x2
∣∣∣∣x=L
= 0
Mf (L, t) = EI∂2w
∂x2
∣∣∣∣x=L
= 0 ∴∂3w
∂x3
∣∣∣∣x=L
= 0
(3-7)
Sendo Q a força cortante e Mf o momento etor. Em x = L, pela
liberdade total de movimento, Q = Mf = 0.
Desenvolvendo as partes da equação (3-5) que carregam as condições
de contorno, vericase que estes termos valem zero para o problema em
questão.
EIdφj(x)
dx
∂2w(x, t)
∂x2
∣∣∣∣L0
= EIdφj
dx(L)
∂2w
∂x2(L, t)︸ ︷︷ ︸
=0
−EI dφj
dx(0)︸ ︷︷ ︸
=0
∂2w
∂x2(0, t) = 0
(3-8)
54
φj(x)∂
∂x
(EI
∂2w(x, t)
∂x2
)∣∣∣∣L0
= EIφj(L)∂3w
∂x3(L, t)︸ ︷︷ ︸
=0
−EI φj(0)︸ ︷︷ ︸=0
∂3w
∂x3(0, t) = 0
(3-9)
A equação (3-5) ca então:∫ L
0
ρA∂2w(x, t)
∂t2φj(x)dx+
∫ L
0
EId2φj(x)
dx2
∂2w(x, t)
∂x2=
∫ L
0
f(x, t)φj(x)dx
(3-10)
Aproximando a solução na forma:
wN(x, t) =N∑
i=1
ai(t)φi(x) (3-11)
Sendo ai coecientes temporais e φi os modos de vibração. A base de
projeção para a dinâmica é formada pelos modos de vibração do problema,
daí o nome do método: Método Modos Supostos.
Os modos de vibração têm propriedades interessantes, como a pro-
priedade de ortogonalidade que surge do fato do problema ser autoadjunto
e gerar autovalores e autovetores reais. Além disso, a resposta dinâmica de
sistemas auto-adjuntos pode ser representada como combinação linear dos
modos normais do sistema.
Os modos de vibração formam a melhor base de projeção para a
dinâmica, no sentido de não existir uma outra base que consiga representar
melhor o problema com o mesmo número de termos. Os modos estão
diretamente relacionados com o problema físico e diretamente relacionados
com a dinâmica. Quando os modos de vibração são conhecidos, ou quando
eles podem ser calculados, eles são a primeira opção para a formação da
base de projeção da dinâmica.
Substituindo a aproximação da solução na equação, já considerando o
erro ortogonal aos termos da aproximação (3-10):
ρAai(t)
∫ L
0
φi(x)φj(x)dx+ EIai(t)
∫ L
0
φ′′i (x)φ′′j (x)dx =
∫ L
0
f(x, t)φj(x)dx
(3-12)
Sendo ai as coordenadas generalizadas e φi os modos de vibração. As
55
matrizes do sistema surgem naturalmente das integrações:
M(φi, φj) = ρA
∫ L
0
φi(x)φj(x)dx , K(φi, φj) = EI
∫ L
0
φ′′i (x)φ′′j (x)dx
F (φj) =
∫ L
0
f(x, t)φj(x)dx
(3-13)
Como foi visto na seção 2.5 o sistema discreto ca:
M a(t) +Ka(t) = F (t) (3-14)
Sendo i = 1, 2, .., N e j = 1, 2, .., N . A EDP original foi transformada
em um sistema de EDOs. A coordenadas generalizadas ai (i = 1, 2, 3, ..N)
são determinadas a partir das condições iniciais. Se uma aproximação é
feita com cinco modos de vibração (N = 5), qualquer contribuição que
modos superiores tenham para a resposta dinâmica não será captada pelo
modelo, se constituindo assim em um erro de aproximação.
Viga com dissipação:
Pode-se acrescentar dissipação ao modelo de viga estudado anterior-
mente. Considere a formulação abaixo:
ρA∂2w(x, t)
∂t2+ c
∂w(x, t)
∂t+
∂2
∂x2
(EI
∂2w(x, t)
∂x2
)= f(x, t) (3-15)
Sendo c o coeciente de amortecimento do material.
A formulação fraca da equação (3-15) é dada por:
∫ L
0
ρA∂2w(x, t)
∂t2φj(x)dx+
∫ L
0
c∂w(x, t)
∂tφjdx+
∫ L
0
EI∂2w(x, t)
∂x2
d2φj(x)
dx2dx
− EIdφj(x)
dx
∂2w(x, t)
∂x2
∣∣∣∣L0
+ φj(x)∂
∂x
(EI
∂2w(x, t)
∂x2
)∣∣∣∣L0
=
∫ L
0
f(x, t)φj(x)dx
(3-16)
Todos os termos das condições de contorno se anulam para as
condições de contorno deste problema. O deslocamento transversal da viga
56
será aproximado por: wN =∑aiφi, logo:
ρAai(t)
∫ L
0
φi(x)φj(x)dx+ ai(t)c
∫ L
0
φi(x)φj(x)dx+
+EIai
∫ L
0
φ′′i (x)φ′′j (x)dx =
∫ L
0
f(x, t)φj(x)dx
(3-17)
As matrizes do sistema surgem naturalmente das integrações:
M(φi, φj) = ρA
∫ L
0
φi(x)φj(x)dx , C(φi, φj) = c
∫ L
0
φi(x)φj(x)dx
K(φi, φj) = EI
∫ L
0
φ′′i (x)φ′′j (x)dx , F (φj) =
∫ L
0
f(x, t)φj(x)dx
(3-18)
Como foi visto na seção 2.5 o sistema discreto ca:
M a(t) + Ca(t) +Ka(t) = F(t) (3-19)
Exemplo 3.1:
O programa din_viga_ms (ver Anexo E) calcula, pelo Método dos
Modos Supostos, a resposta dinâmica de uma viga engastadalivre, dada
uma certa precisão.
Considere a viga representa na Figura 3.3:
Sistema Discreto
Resposta
PrecisãoSatisfeita?
Aumentar número de elementosda base usados na aproximação
ok!
Não
Sim
Material - Aço:E=200GPa (Módulo de Elasticidade =7850 kg/m (massa específica)c=5000 Ns /N (amortecimento)r 3
2
L=3m F(t)=120 (200(2 ).t) Np
h=5cm
b=10cm
Figura 3.3: Viga engastada, excitada na extremidade livre
Sem dissipação: Usando inicialmente 5 modos na aproximação, tempo
de 0,01 segundos de simulação e condições iniciais nulas (w = 0 e w = 0).
Com apenas 6 modos de vibração o erro registrado foi de 2,7% para a
resposta do ponto da extremidade da viga, x = L, Figura 3.4.
57
Figura 3.4: Resposta dinâmica na extremidade livre
O sistema de EDOs é integrado pela subrotina ode45 do MATLAB.
O número de modos necessários para aproximar um sistema pode ser
surpreendentemente pequeno. Se a resposta dinâmica for simples, três ou
quatro modos serão sucientes para representar o comportamento de uma
viga, por exemplo. Quanto mais complicada a resposta dinâmica, mais
modos serão necessários para representá-la.
Com dissipação: A resposta dinâmica para o ponto da extremidade da
viga, x = L, com condições iniciais nulas é mostrada na Figura 3.5.
58
Figura 3.5: Resposta dinâmica na extremidade livre
O número de modos usados na simulação foi de 6, apresentando um
erro de 0,4% para a resposta dinâmica no ponto x = L. O programa calcula
também as freqüências naturais e os coecientes de amortecimento modais.
Freqüências naturais (Hz):
ωn =
4.5207
28.3308
79.3279
155.4526
256.9769
383.8848
Freqüências de oscilação do sistema amortecido para cada modo de
vibração:
ωd =
4.4055
28.3127
79.3214
155.4493
256.9749
383.8834
Repare que estas freqüências são um pouco menores do que as fre-
qüências naturais, o que é consistente com a relação ωdi= ωi
√1− ξ2
i , que
será vista no capítulo 4. Coecientes de amortecimento modal, ξi:
59
ξ =
0.2242
0.0358
0.0128
0.0065
0.0039
0.0026
3.3Método dos Elementos Finitos
Em muitos sistemas mecânicos não se tem uma expressão analítica
para calcular os modos de vibração. Nestes casos os modos podem ser
calculados através do Método dos Elementos Finitos (MEF). No MEF as
funções que formam a base para o Método de Galerkin são funções simples
que não têm relação alguma com a dinâmica do sistema mecânico.
O MEF foi desenvolvido inicialmente para a análise estrutural durante
os anos de 1950 e 1960. O MEF é um método de discretização que aproxima
tanto a formulação quanto o domínio de um sistema contínuo. A estrutura é
dividida em partes menores chamadas de elementos nitos. Cada elemento
é usualmente muito simples e tem uma equação de movimento associada.
As soluções das equações dos elementos são aproximadas por uma
combinação linear de polinômios de baixa ordem. Cada uma destas soluções
polinomiais individuais é tornada compatível com a solução adjacente
(condição de continuidade) nos nós comuns aos elementos. Estas soluções
são reunidas através de um procedimento, resultando em matrizes globais
(matrizes de massa e rigidez globais, por exemplo) as quais descrevem a
dinâmica de uma estrutura como um todo. Neste trabalho esta técnica é
importante para ser usada em sistemas em que não se conhecem os modos
de vibração.
Ponto de vista global:
O MEF é empregado na formulação fraca (ou variacional) do pro-
blema, e a aproximação é calculada através do Método de Galerkin. A
resposta é representada por uma base que está contida em um espaço de
Hilbert. Existe todo um esquema de aproximação.
Considere o equacionamento de uma barra xa em uma extremidade
e livre na outra, discutido no capítulo 2. Depois da aplicação do Método de
60
Galerkin, já incorporando as condições de contorno, chega-se a:
ai(t)ρA
∫ L
0
φi(x)φj(x)dx+ai(t)EA
∫ L
0
dφi(x)
dx
dφj(x)
dxdx =
∫ L
0
Pi(x, t)φj(x)dx
(3-20)
A particularidade do MEF está no fato de que cada função teste é
diferente de zero em apenas uma pequena parte do domínio.
As funções teste são dadas por funções simples chamadas polinômios
de interpolação. O resultado da integral de um polinômio de interpolação é
diferente de zero apenas em um pequeno intervalo:∫ L
0
φ1dx =
∫ x2
x1
φ1dx (3-21)
Denindo as funções teste desta maneira, os coecientes da aproxi-
mação, equação (3-22), tem um signicado especial:
u(xj, t) =N∑
i=1
ai(t)φi(xj) (3-22)
Sendo xj a posição do jésimo nó. Cada polinômio de interpolação
vale um em um nó do intervalo e zero nos demais nós:
φi(xj) = δij (3-23)
Sendo δij = 1 quando i = j e 0 quando i 6= j. Logo, o valor do
coeciente ai é igual ao valor do deslocamento, u(xi, t), do iésimo nó:
u(xj, t) =N∑
i=1
ai(t)φi(xj) =N∑
i=1
ai(t)δij = aj(t)
∴ ai(t) = u(xi, t)
(3-24)
Ponto de vista elementar:
As matrizes globais de um sistema são construídas a partir da mon-
tagem das matrizes dos elementos. O elemento três de uma barra, por exem-
61
plo, terá matrizes:
M (3) = ρA
∫ x4
x3
φ3φ3dx
∫ x4
x3
φ3φ4dx∫ x4
x3
φ4φ3dx
∫ x4
x3
φ4φ4dx
(3-25)
K(3) = EA
∫ x4
x3
dφ3
dx
dφ3
dxdx
∫ x4
x3
dφ3
dx
dφ4
dxdx∫ x4
x3
dφ4
dx
dφ3
dxdx
∫ x4
x3
dφ4
dx
dφ4
dxdx
(3-26)
F (3) =
∫ x4
x3
P3φ3dx∫ x4
x3
P4φ4dx
(3-27)
A forma de cada função de interpolação é a mesma para todos os
elementos. Isto sugere uma troca de coordenadas para descrever as funções
de forma local. Esta estratégia é usada para que resultados das integrações
sejam usados de forma sistemática. Para uma aproximação linear podese
convencionar que cada elemento se estende de ξ = −1 a ξ = 1, como mostra
a gura 3.6.
)(
)(
1 xfe
)()(
2 xfe
x
+1 -1
1 2
Figura 3.6: Aproximação linear
As funções elementares para uma aproximação linear são dadas por:
φ(e)1 (ξ) =
1− ξ
2φ
(e)2 (ξ) =
1 + ξ
2(3-28)
Todo os elementos são mapeados para este elemento padrão. As
matrizes do elemento 3, equação (3-25 a 3-27), cam, por exemplo:
M (3) =ρA∆x(e)
6
1 −1
−1 1
(3-29)
62
K(3) =EA
∆x(e)
1 −1
−1 1
(3-30)
K(3) =1
∆x(e)
−P3
P4
(3-31)
Sendo ∆x(e) o tamanho do elemento.
Para construir a matriz global do sistema, as matrizes dos elementos,
(3-29) e (3-30) são montadas, por exemplo:
K =
K(1)1,1 K
(1)1,2 0 0 ..
K(1)2,1 K
(1)2,2 +K
(2)1,1 0 0 ..
0 K(2)2,1 K
(2)2,2 +K
(3)1,1 0 ..
0 0 K(3)2,1 K
(3)2,2 ..
0 0 0 0 ..
..
Repare que muitos elementos das matrizes do sistema serão zero, pois
as funções de interpolação são diferentes de zero em apenas uma pequena
parte do domínio:
K2,5 =
∫ 1
0
dφ2
dx
dφ5
dxdx = 0 (3-32)
63
Figura 3.7: Polinômios de interpolação diferentes de zero apenas em uma
pequena parte do domínio
O número de funções de interpolação dene a dimensão do espaço
e está associado ao número de nós da malha. Quanto mais renada for a
malha, maior será a dimensão do espaço de funções. Infelizmente aumentar
a dimensão do espaço não garante aumento na precisão da aproximação.
Isto só será válido se a malha renada estiver contida (embedded) na malha
original.
Além do artifício de aumentar o número de elementos, a precisão
da solução obtida através da discretização pelo MEF pode ser melhorada
aumentando-se o grau da função de interpolação. O aproximação pode ser
linear, quadrática, cúbica, etc. Considere o elemento de barra representado
na Figura 3.8:
64
Figura 3.8: Elemento com aproximação (a) Linear, 2 nós, (b) Quadrática, 3
nós e (c) Cúbica, 4 nós
Uma aproximação quadrática necessita de três pontos por elemento,
Figura 3.9.
)()(
3 xfe
)()(
1 xfe
)()(
2 xfe
x
+1 -1
1 2 3
Figura 3.9: Aproximação quadrática
u(e) = φ1u(x1) + φ2u(x2) + φ3u(x3) (3-33)
Sendo
φ1(x) = 12(1− x)− 1
2(1− x2);
φ2(x) = 12(1 + x)− 1
2(1− x2);
φ3(x) = (1− x2);
Aproximação cúbica:
u(e) = φ1u(x1) + φ2u(x2) + φ3u(x3) + φ4u(x4) (3-34)
65
Sendo
φ1(x) = 12(1− x)− 1
2(1− x2) + 1
16(−9x3 + x2 + 9x− 1);
φ2(x) = 12(1 + x)− 1
2(1− x2) + 1
16(9x3 + x2 − 9x− 1);
φ3(x) = (1− x2) + 116
(27x3 + 7x2 − 27x− 7);
φ4(x) = 116
(−27x3 − 9x2 + 27x+ 9);
3.3.1Viga
Figura 3.10: Elemento de viga
Considere um elemento de viga modelo Euler-Bernoulli representada
na Figura 3.10. Fazendo uma aproximação para o deslocamento transversal
w(x) da forma:
w(e) = φ1(x)v1 + φ2(x)θ1 + φ3(x)v2 + φ4(x)θ2 (3-35)
Sendo θ =∂w(e)
∂x, w(e)(0) = v1, w(e)(L) = v2, θ(0) = θ1, θ(L) = θ2, e
φi as funções de interpolação.
As funções de interpolação são funções de forma Hermitiana. Forma
Hermitiana signica que há imposição tanto para os valores das funções
66
como para valores das derivadas. Estas funções estão representadas na
Figura 3.11:
Figura 3.11: Funções Hermitianas
Sendo as funções Hermitianas H(x) as funções de interpolação, φ(x):
φ1(x) = 1− 3x2/l2 + 2x3/l3;
φ2(x) = x− 2x2/l + x3/l2;
φ3(x) = 3x2/l2 + 2x3/l3;
φ4(x) = x2/l + x3/l2;
Φ = [φ1 φ2 φ3 φ4]
Reescrevendo a formulação fraca da viga, equação 3-12:
ρAai(t)
∫ L
0
φi(x)φj(x)dx+ EIai(t)
∫ L
0
φ′′i (x)φ′′j (x)dx =
∫ L
0
fφjdx (3-36)
As matrizes de massa e de rigidez de cada elemento são obtidas por:
M (1) = ρA
∫ x2
x1φ1φ1dx
∫ x2
x1φ1φ2dx
∫ x2
x1φ1φ3dx
∫ x2
x1φ1φ4dx∫ x2
x1φ2φ1dx
∫ x2
x1φ2φ2dx
∫ x2
x1φ2φ3dx
∫ x2
x1φ2φ4dx∫ x2
x1φ3φ1dx
∫ x2
x1φ3φ2dx
∫ x2
x1φ3φ3dx
∫ x2
x1φ3φ4dx∫ x2
x1φ4φ1dx
∫ x2
x1φ4φ2dx
∫ x2
x1φ4φ3dx
∫ x2
x1φ4φ4dx
(3-37)
67
K(1) = EI
∫ x2
x1φ1φ1dx
∫ x2
x1φ′′1φ
′′2dx
∫ x2
x1φ′′1φ
′′3dx
∫ x2
x1φ′′1φ
′′4dx∫ x2
x1φ′′2φ
′′1dx
∫ x2
x1φ′′2φ
′′2dx
∫ x2
x1φ′′2φ
′′3dx
∫ x2
x1φ′′2φ
′′4dx∫ x2
x1φ′′3φ
′′1dx
∫ x2
x1φ′′3φ
′′2dx
∫ x2
x1φ′′3φ
′′3dx
∫ x2
x1φ′′3φ
′′4dx∫ x2
x1φ′′4φ
′′1dx
∫ x2
x1φ′′4φ
′′2dx
∫ x2
x1φ′′4φ
′′3dx
∫ x2
x1φ′′4φ
′′4dx
(3-38)
Fazendo a integração das funções de interpolação chega-se a:
M (e) =
∫ x2
x1
ρAΦT Φdx =ρAl
400
156 22l 54 −13l
22l 4l2 13l −3l2
54 13l 156 −22l
−13l −3l2 −22l 4l2
(3-39)
K(e) =
∫ x2
x1
EIΦ′′T Φ′′dx =EI
l3
12 6l −12 6l
6l 4l2 −6l 2l2
−12 −6l 12 −6l
6l 2l2 −6l 4l2
(3-40)
Para obter as matrizes globais de massa e de rigidez é necessário fazer
a montagem das matrizes dos elementos.
Exemplo 3.2:
O programa din_viga_ef (ver Anexo E) calcula, discretizando pelo
MEF, a resposta dinâmica de uma viga engastadalivre, dada uma certa
precisão.
Considere a viga representada na Figura 3.12:
68
Sistema Discreto
Resposta
PrecisãoSatisfeita?
Aumentar número de elementosda base usados na aproximação
ok!
Não
Sim
Material - Aço:E=200GPa (Módulo de Elasticidade =7850 kg/m (massa específica)r 3
L=3m F(t)=120 (200(2 ).t) Np
h=5cm
b=10cm
Figura 3.12: Viga engastada, excitada na extremidade livre
Para um tempo de 0, 01 segundos de simulação e condições iniciais
nulas (w = 0 e w = 0), a resposta do ponto em que o forçamento é aplicado
é:
Figura 3.13: Resposta dinâmica em x=L
Para uma precisão de 2, 66% foram necessários 60 elementos nitos
(usando 20 elementos iniciais e passo de 20 elementos). Este mesmo exem-
plo foi feito anteriormente (seção 3.2) utilizando-se o Método dos Modos
Supostos. Com 7 modos a precisão foi de 0, 5%. A Figura 3.14 mostra as
duas respostas dinâmicas, sendo a de vermelho a aproximada pelo Método
dos Modos Supostos.
69
Figura 3.14: Resposta dinâmica em x=L
O erro percentual entre as respostas é de 5, 7%. O tempo de proces-
samento do programa discretizando pelo MEF foi mais de 100 vezes maior
do que o tempo de processamento discretizando pelo Método dos Modos
Supostos. Isto se deve ao tamanho das matrizes geradas: com 60 elemen-
tos a matriz de massa, por exemplo, tem dimensão 183 (3x61), sendo 61 o
números de nós e 3 o número de graus de liberdade de cada nó.
3.3.2Transformação de Coordenadas
Na seção anterior as matrizes M e K foram calculadas com respeito
a eixos e coordenadas locais, mas freqüentemente membros de estruturas
não estão alinhados com respeito a um sistema de coordenadas global de
referência.
70
Figura 3.15: Eixos local e global para um elemento de barra
A Figura 3.15 mostra os deslocamentos de um elemento de barra com
relação a dois sistemas de coordenadas: um local e um global. No sistema
de coordenadas local são necessários apenas dois deslocamentos: u1 e u2.
Para representar o mesmo elemento de barra num sistema de coordenadas
global (X,Y), são necessários quatro deslocamentos: u1, u2, u3 e u4.
O deslocamento u1 = u1cos(θ) + u2sen(θ)
O deslocamento u2 = u3cos(θ) + u4sen(θ)
Escrevendo em termos matriciais:
u = Tu (3-41)
Sendo: uT = [u1 u2], uT = [u1 u2 u3 u4]
T =
[cosθ senθ 0 0
0 0 cosθ senθ
]
A Figura 3.16 mostra o caso de um elemento com três graus de
liberdade em cada nó:
71
Figura 3.16: Elemento com três graus de liberdade em cada nó
u1
u2
u3
u4
u5
u6
=
cosθ senθ 0 0 0 0
0 0 1 0 0 0
0 0 0 cosθ senθ 0
0 0 0 −senθ cosθ 0
0 0 0 0 0 1
u1
u2
u3
u4
u5
u6
(3-42)
Para a transformação das matrizes M e K, podem ser aplicados os
conceitos de energia potencial V e cinética T, considerando que a energia
se mantém constante depois da transformação:
V = 12(uTKu) = 1
2(uTKu)
T = 12(uTMu) = 1
2(uTMu)
Substituindo u = Tu:
V = 12[(uTT T )K (Tu)] = 1
2uT (T TK T )u
Logo:
K = T TKT (3-43)
M = T TMT (3-44)
72
O anexo B apresenta programas que foram desenvolvidos para a
Análise Modal de estruturas bidimensionais, como um 'L' e um 'Triângulo'.
3.4Carregamento
O vetor de carregamento do sistema surge naturalmente da formulação
fraca, através da integração:
F (x, t) =N∑
j=1
∫ L
0
f(x, t)φj(x)dx (3-45)
- Se houver força concentrada:
Figura 3.17: Força concentrada
f(x) = δ(x− x1) (3-46)
Sendo x1 o local onde a força é aplicada e δ a função impulso unitário
ou distribuição de Diracδ. f(x) vale zero para todos os valores de x, exceto
em x = x1. O vetor de carregamento F será então:
F =N∑
j=1
∫ L
0
δ(x− x1)φjdx =N∑
j=1
φj(x1) (3-47)
- Se houver força distribuída:
73
Figura 3.18: Força distribuída
f(x) =
(fL − f0
L
)x+ f0 (3-48)
Sendo f0 a força em x=0 e fL a força em x=L. O vetor de carregamento
F será então:
F =N∑
j=1
∫ L
0
[(fL − f0
L
)x+ f0
]φjdx (3-49)
Exemplo 3.3: Carregamento estático: força concentrada e força dis-
tribuída.
Considere os carregamentos representados nas duas Figuras 3.19 e
3.20.
Figura 3.19: Força concentrada na extremidade de uma viga
74
Figura 3.20: Força distribuída uniformemente ao longo de uma viga
A expressão analítica para o cálculo do deslocamento do ponto da
extremidade livre para um força concentrada na extremidade, Figura 3.19,
é, [38]:
δc =FL3
3EI(3-50)
A expressão analítica para o cálculo do deslocamento do ponto da
extremidade livre para uma força distribuída uniformemente ao longo da
viga, Figura 3.20, é, [38]:
δd =FL4
8EI(3-51)
O vetor com os deslocamentos, w, é calculado pela equação (3-52).
w = K−1F (3-52)
Sendo K a matriz de rigidez do sistema e F o vetor das forças. Os
programas est_viga_ms e est_viga_ms calculam o deslocamento do ponto
da extremidade livre de uma viga engastadalivre para uma força concen-
trada na extremidade livre e para uma força distribuída uniformemente ao
longo de uma viga.
O programa est_viga_ms faz a discretização pelo Método dos Modos
Supostos. O resultado para a aproximação com 5 modos de vibração é:
δc = −4, 33.10−5 m e δd = −4, 88.10−5 m.
O programa est_viga_ef faz a discretização pelo MEF. O resultado
para a aproximação com 20 elementos nitos é:
δc = −4, 32.10−5 m e δd = −4, 86.10−5 m.
75
O resumo dos resultados está apresentado na Tabela 3.1:
δc (m) δd (m)
Expressão analítica −4, 32.10−5 −4, 86.10−5
MEF (20 elementos) −4, 32.10−5 −4, 86.10−5
Modos Supostos (5 modos) −4, 33.10−5 −4, 88.10−5
Tabela 3.1: Comparação entre valores obtidos para o deslocamento
3.5Condições de Contorno
As condições de contorno se incorporam à equação do sistema na
formulação fraca. Partindo da equação de uma barra:
ρA(x)∂2u(x, t)
∂t2=
∂
∂x
(A(x)E
∂u(x, t)
∂x
)+ f(x, t) (3-53)
Formulação fraca:
∫ L
0
ρA∂2u(x, t)
∂t2φdx = φ
(EA
∂u
∂x
)∣∣∣∣L0
−∫ L
0
(EA
∂u
∂x
)∂φ
∂x+
∫ L
0
f(x, t)φdx
(3-54)
As condições de contorno entram no termo em destaque da equação
(3-54):
φ
(EA
∂u(x, t)
∂x
)∣∣∣∣L0
(3-55)
Figura 3.21: Condições de contorno no extremo L. (a) ligação elástica; (b)
massa pontual
76
Exemplo de condições de contorno:
1)Barra xalivre:
u(0) = 0 e EA∂u
∂x
∣∣∣∣L
= 0 (3-56)
2)Barra xamola, Figura 3.21a:
u(0) = 0 e EA∂u
∂x
∣∣∣∣L
= ku|L (3-57)
3)Barra xamassa, Figura 3.21b:
u(0) = 0 e EA∂u
∂x
∣∣∣∣L
= m∂2u
∂t2
∣∣∣∣L
(3-58)
As condições de contorno são incorporadas às matrizes do sistema.
Ligação elástica:
K =
∫ L
0
EAφ′′i φ′′jdx+ φ(L) ku|L (3-59)
Massa pontual:
M =
∫ L
0
ρAφiφjdx+ φ(L) m∂2u
∂t2
∣∣∣∣L
(3-60)
Para validar o programa desenvolvido, através da discretização pelo
MEF, os resultados computados foram comparados com os resultados
obtidos das seguintes expressões, [18] e [3], equações (3-613-66).
Barra xalivre:
Freqüências naturais:
ωn =(2n− 1)πc
2L(3-61)
Sendo c =√
Eρe n=1,2,3,..
Modos de vibração:
ωn = sin
((2n− 1)πx
2L
)(3-62)
77
Barra xamola:
Freqüências naturais:
ωn =λnc
Lλncot(λn) = − kL
EA(3-63)
Modos de vibração:
sin
(λnx
L
)(3-64)
Barra xamassa:
Freqüências naturais:
ωn =λnc
Lcot(λn) =
m
ρALλn (3-65)
Modos de vibração:
sin
(λnx
L
)(3-66)
Para aproximar as equações transcendentais (3-63 e 3-65) e obter os
λn das equações, foi usado o Método de Newton, [6], que é um método de
iteração para a procura de zeros de uma função diferenciável.
Exemplo 3.4:
O programa modal_barra_ef (ver Anexo E) calcula, discretizando
pelo MEF, os modos e freqüências naturais de uma barra xalivre / xa
mola / xamassa, para uma determinada precisão.
Considere uma barra de um metro de comprimento, de seção 5x10 cm2,
módulo de elasticidade de 200 MPa e densidade 7850 kg/m3. Erro estipulado
de 0, 005% para o quinto modo de vibração. O número de elementos inicial
é 50 e o passo de discretização é de 50.
O número de elementos necessários para a precisão requerida é de 350
elementos. As cinco primeiras freqüências naturais (Hz) são:1262
3786
6310
8834
1.1358
(3-67)
78
Os cinco primeiros modos de vibração estão representados na Figura
3.22:
Figura 3.22: Cinco primeiros modos de uma barra xada em uma extremi-
dade e livre na outra
O erro do quinto modo de vibração para cada passo da simulação é
dado por:
Figura 3.23: Convergência da aproximação
Exemplo 3.5:
79
Considere os mesmos dados do Exemplo 3.4. Mas agora uma mola
(k = 1GNm) é colocada na extremidade da barra.
A seguir os λn calculados pela equação transcendental
utilizando-se o método de Newton (precisão de 0, 0001): λn =
2.0288; 4.9132; 7.9787; 11.0855; 14.2074.O número de elementos necessários para a precisão requerida é de 350
elementos. As cinco primeiras freqüências naturais (Hz) calculadas com os
λn computados são: 1630
3947
6410
8905
1.1413
(3-68)
As cinco primeiras freqüências naturais (Hz) calculadas, discretizando
pelo MEF são: 1630
3947
6410
8906
1.1414
(3-69)
Observa-se que os resultados estão bem coerentes entre as aproxi-
mações; erro menor que 0, 1%. As freqüências naturais são maiores, com-
paradas com as freqüências obtidas de uma barra xalivre. Este resultado
é coerente com a expressão: ωn =
√K
M, ou seja, quanto maior a rigidez,
maior a freqüência natural (mantendose M xo).
Os cinco primeiros modos de vibração estão representados na Figura
3.24:
80
Figura 3.24: Cinco primeiros modos de uma barra xa em uma extremidade
e com apoio elástico na outra
O erro do quinto modo de vibração para cada passo da simulação é
dado por:
Figura 3.25: Convergência da aproximação
Primeiro modo de vibração:
81
Figura 3.26: Primeiro modo de vibração
Na Figura 3.26 a curva azul representa o primeiro modo de vibração
calculado discretizando pelo MEF, a curva verde representa o primeiro modo
calculado através da equação transcendental (3-63) e a curva vermelha
representa o primeiro modo para uma barra xa em uma extremidade e
livre na outra. Segundo modo de vibração:
Figura 3.27: Segundo modo de vibração
Exemplo 3.6:
82
Considere os mesmos dados do Exemplo 3.5. Mas agora uma massa
(m = 10kg) é colocada na extremidade da barra. A seguir os λn calculados
pela equação transcendental utilizando-se o método de Newton (precisão de
1e− 4): λn = 1.2601; 3.9267; 6.8063; 9.8055; 12.8625.0 número de elementos necessários para a precisão requerida é de 350
elementos. As cinco primeiras freqüências naturais (Hz) calculadas com os
λn computados são: 1012
3154
5468
7877
1.0333
(3-70)
As cinco primeiras freqüências naturais (Hz) calculadas discretizando
pelo MEF são: 1012
3155
5468
7878
1.0334
(3-71)
Observa-se que os resultados estão bem coerentes entre as aproxi-
mações, erros menores que 0, 1%. As freqüências naturais sao menores, com-
paradas com as freqüências obtidas de uma barra xalivre. Este resultado
é coerente com a expressão: ωn =
√K
M, ou seja, quanto maior a massa,
menor a frequência natural (mantendose K xo).
Os cinco primeiros modos de vibração estão representados na Figura
3.28:
83
Figura 3.28: Cinco primeiros modos barra xamassa
O erro do quinto modo de vibração para cada passo da simulação é
dado por:
Figura 3.29: Convergência da aproximação
Primeiro modo de vibração:
84
Figura 3.30: Primeiro modo de vibração
Na Figura 3.30 a curva azul representa o primeiro modo de vibração
calculado discretizando pelo MEF, a curva verde representa o primeiro modo
calculado através da equação transcendental (3-65) e a curva vermelha
representa o primeiro modo para uma barra xa em uma extremidade e
livre na outra. Segundo modo de vibração:
Figura 3.31: Segundo modo de vibração
85
3.6Conclusões
Neste capítulo o Método de Galerkin foi detalhado a partir das
escolhas da base de projeção da dinâmica, formada pelas funções teste. Um
sistema contínuo foi discretizado e um esquema para calcular as matrizes
do sistema discreto foi elaborado. No próximo capítulo um sistema linear
discreto, com matrizes de massa, de amortecimento e de rigidez, é analisado
através da Análise Modal.
4Análise Modal
4.1Introdução
No capítulo 2 os conceitos de modos de vibração e freqüências naturais
foram introduzidos. A Análise Modal é uma forma de caracterizar um
sistema dinâmico em termos de três grupos de parâmetros: freqüências
naturais, modos de vibrações e amortecimentos modais.
Figura 4.1: Separação dos parâmetros modais
A Figura 4.1 ilustra bem o objetivo da Análise Modal. O movimento
vibratório complicado de um sino pode ser representado como uma soma
de movimentos simples. Cada movimento simples é um modo de vibração
87
e cada modo de vibração tem associado a ele uma freqüência natural e um
coeciente de amortecimento.
Os modos normais formam a melhor base de projeção para o Método
de Galerkin.
4.2Modos Normais
Sistema Conservativo:
Os modos normais estão associados ao sistema conservativo, equação
(4-1).
M x(t) +Kx(t) = 0 (4-1)
Sendo M uma matriz simétrica positiva denida, K uma matriz
simétrica positiva semidenida e x o vetor composto pelas coordenadas do
sistema. Os modos e freqüências são características intrínsecas do sistema,
por isso o sistema homogêneo é analisado (F = 0).
A equação (4-1) admite solução do tipo:
x(t) = uiejωit (4-2)
Sendo j =√−1. Já adiantando o resultado, ωi é a iésima freqüência
natural e ui é o iésimo modo de vibração. Substituindo (4-2) em (4-1)
obtém-se:
(−ω2iM +K)uie
jωit = 0 (4-3)
Na busca da solução não trivial em que ui 6= 0, como [ejωit] 6= 0,
chegase a:
(−ω2iM +K)ui = 0 (4-4)
Este é um problema de autovalor. É através deste problema que as
freqüências naturais e os modos de vibração são calculados. Um sistema
autoadjunto, ou seja, com M positiva denida e K positiva semidenida,
88
gera autovalores (ω2i ) e autovetores (ui) reais:
ω2iMui = Kui
〈ω2iMui,ui〉 = 〈Kui,ui〉
ω2i =
〈Kui,ui〉〈Mui,ui〉
≥ 0
(4-5)
Um autovalor pode ter um ou mais autovetores associados a ele. Por
exemplo, um autovalor zero, ou seja, uma freqüência natural zero, pode ter
associado a ele mais de um movimento de corpo rígido. Para um sistema
tridimensional o autovalor zero pode estar associado a até seis autovetores
linearmente independentes, Figura 4.2.
Figura 4.2: Movimento de corpo rígido
Sistemas com simetrias (simetria esférica, por exemplo) também apre-
sentam autovalores com multiplicidade, isto é, com mais de um autovetor
associado.
A matriz modal, U , é dada por:
U =
| | .. |u1 u2 .. uN
| | .. |
(4-6)
As matrizes M e K são diagonalizadas pela matriz modal:
UTMU =
mi
e UTKU =
ki
(4-7)
Cada autovetor tem uma única direção associada, mas sua amplitude
não é denida. O processo de normalização é realizado como uma forma de
89
xar a amplitude dos modos de vibração. Para normalizar U com relação
às matrizes M e K prémultiplicase a matriz U por um fator γi:
φi = γiui −→ γ2i u
Ti Mui = I (4-8)
Sendo γi = 1/√mi e I uma matriz identidade. A matriz modal
normalizada, Φ, é dada por:
Φ =
| | .. |φ1 φ2 .. φN
| | .. |
(4-9)
Uma matriz modal normalizada, Φ, diagonaliza as matrizesM e K da
seguinte forma:
ΦTMΦ =
1
e ΦTKΦ =
ω2
i
(4-10)
Para resolver o sistema dinâmico pode-se tomar proveito da matriz
modal fazendo uma mudança de coordenadas, de x para q:
x(t) = Φq(t) (4-11)
Um sistema dinâmico com as novas coordenadas pode ser escrito como:
MΦq(t) +KΦq(t) = F(t) (4-12)
Prémultiplicando por ΦT :
ΦTMΦq(t) + ΦTKΦq(t) = ΦTF(t) (4-13)
Desta forma as equações de movimento se desacoplam:
q(t) +
ω2
i
q(t) = ΦTF(t) (4-14)
Sendo q(t) o vetor composto pelas coordenadas modais. O problema
foi reduzido a um sistema de EDO's independentes.
90
Cada modo de vibração caracteriza uma família de soluções periódi-
cas.:
x(t) =N∑
i=1
qi(t)φi , ou
x(t) =N∑
i=1
aφicos(ωit+ ϕ)
(4-15)
Sendo ωi a iésima freqüência natural e φi o iésimo modo de vi-
bração. Os coecientes a e ϕ são determinados a partir das condições
iniciais. Repare que as condições iniciais podem ser escolhidas de forma que
cada um dos modos pode ser excitado independentemente.
Amortecimento proporcional:
A modelagem de dissipações em sistemas dinâmicos não é simples.
Uma forma de modelar um sistema com dissipação é considerar a matriz de
amortecimento uma combinação linear das matrizes de massa e de rigidez.
Assim os modos normais podem ser usados para representar as respostas
dinâmicas. Considere um sistema dinâmico linear com amortecimento pro-
porcional:
M x(t) + Cpx(t) +Kx(t) = 0 (4-16)
Sendo Cp = αM + βK. Como Cp é combinação linear de M e K, a
mesma matriz modal que diagonaliza M e K, diagonaliza Cp.
Depois de uma mudança de coordenadas, x(t) = Φq(t), chegase ao
sistema:
ΦTMΦq(t) + ΦTCpΦq(t) + ΦTKΦq(t) = ΦTF(t) (4-17)
O sistema desacoplado pode ser escrito da seguinte maneira:
q(t) +
2ξiωi
q(t)
ω2
i
q(t) = ΦTF(t) (4-18)
Sendo:
ξi = ci/cicrita razão de amortecimento do iésimo modo;
cicrit= 2miωi = 2
√miki o amortecimento crítico do iésimo modo.
91
Cada modo de vibração caracteriza uma família de soluções periódicas:
x(t) =N∑
i=1
qi(t)φi , ou
x(t) =N∑
i=1
aφie−ξiωitcos(ωdi
t+ ϕ)
(4-19)
Sendo:
ωdi= ωi
√1− ξ2
i a iésima freqüência de oscilação do sistema amortecido
para o iésimo modo de vibração.
Podese pensar nos modos normais como ondas estáticas com nós
xos, Figura 4.3.
w(x)
x
nó estacionário
Figura 4.3: Modo normal
Exemplo 4.1:
Este exemplo encontrase no programa modal (ver Anexo E) . Basta
entrar com as matrizes M e K, e as constante α e β, Cp = αM + βK. Este
programa calcula as freqüências naturais, a matriz modal (com os modos
normalizados com respeito a M) e os coecientes de amortecimento modais.
Estes cálculos são realizados conforme os passos a seguir. Dadas as matrizes:
M =
[2 0
0 3
], K =
[2 −1
−1 1
], Cp = 1
6M + 1
10K
Resolvese um problema de autovalor: (−ω2iM +K)ui = 0
92
Os autovalores, ω2i , são:[
1, 1937 0
0 0, 1396
]
As freqüências naturais por conseqüência são: ω1 = 1, 0926 rd/s e
ω2 = 0, 3736 rd/s.
Os autovetores, ui, são:
U =
[0, 9325 0, 5025
−0, 3613 0, 8646
]
UTMU =
[2, 1305 0
0 2, 7475
]=
mi
UTKU =
[2, 5432 0
0 0, 3836
]=
ki
Normalizando os autovetores com respeito à matriz M :
φi =1
√mi
ui
Logo:
Φ =
[0, 6388 0, 3031
−0, 2475 0, 5216
]
ΦTMΦ = I
ΦTKΦ =
[1, 1937 0
0 0, 1396
]=
ω2
i
ΦTCpΦ =
[0, 2860 0
0 0, 1806
]=
2ξiωi
Logo, a matriz com os coecientes de amortecimento modais é:
93
ξ =
[0, 1309 0
0 0, 2417
]
wd1 = 1, 0832 rd/s e wd2 = 0, 3626 rd/s.
4.3Modos Complexos
Os modos complexos aparecem em sistemas: com amortecimento não
proporcional, com forças giroscópicas, com forças aerodinâmicas, etc.
4.3.1Forças Giroscópicas
Forças giroscópicas são conservativas, porém elas acoplam o movi-
mento em duas direções. Considere o sistema:
M x(t) +Gx(t) +Kx(t) = 0 (4-20)
SendoG uma matriz antisimétrica, chamada de giroscópica. Para este
tipo de sistema não vale a Análise Modal Clássica. Considerando a matriz
K positiva denida, a equação (4-20) admite solução do tipo:
x(t) = uiejωit (4-21)
Sendo ωi a iésima freqüência natural e ui o iésimo modo complexo.
Substituindo (4-21) em (4-20) obtémse:
−ω2iMui + jωiGui +Kui = 0 (4-22)
Témse o seguinte problema de autovalor:
(−ω2iM + jωiG+K)u = 0 (4-23)
Para resolver este problema, a equação (4-20) deve ser reescrita da
seguinte forma, [26]:
M∗y(t) = −G∗y(t) (4-24)
94
Sendo:
y(t) =
[x(t)
x(t)
]M∗ =
[K 0
0 M
]= M∗T G∗ =
[0 −KK G
]= −G∗T
(4-25)
A equação (4-24) admite solução do tipo:
y(t) = Yiejωt (4-26)
Sendo Yi um vetor constante. Substituindo (4-26) em (4-24) chegase
ao seguinte problema de autovalor:
[jωiM∗ +G∗]Yi = 0 (4-27)
Novamente aparece na equação a unidade imaginária j. Para contornar
esta situação separase a parte real e a parte imaginária: Yi = YRi + jYI
i :
ωiM∗YR
i = −G∗YIi (4-28)
ωiM∗YI
i = G∗YRi (4-29)
Resolvendo a equação (4-29) para YRi e introduzindo na equação (4-
28) o resultado, depois resolvendo a equação (4-28) para YIi e introduzindo
o resultado novamente na equação (4-29) chega-se a:
K∗YRi = λiM
∗YRi K∗YI
i = λiM∗YI
i (4-30)
Sendo λi = ω2i e:
K∗ =
[KM−1K KM−1G
GTM−1K GTM−1G
]= K∗T (4-31)
Através deste artifício matemático, não só o problema de autovalor
deixou de ter a parte imaginária j como passou a ser função de duas
matrizes simétricas positivas denidas. O resultado nal são 2n autovalores
que aparecem em pares conjugados:
λi = +jωi i = 1, 2, .., n
λi+n = −jωi i = 1, 2, .., n(4-32)
95
E os autovetores:
Yi =
[φi
φiλi
](4-33)
Os modos complexos são dados por: φi = φRi + jφI
i . Cada modo
complexo caracteriza uma família de soluções periódicas:
x(t) = a[φRi cos(ωit− ϕ)− φI
i sen(ωit− ϕ)] (4-34)
Sendo os coecientes a e ϕ determinados pelas condições iniciais.
4.3.2Amortecimento Não-Proporcional
Considere o sistema:
M x(t) + Cx(t) +Kx(t) = 0 (4-35)
Sendo C simétrica positiva semidenida. Reescrevendo o sistema.[C M
M 0
]y(t) +
[K 0
0 −M
]y(t) = 0 (4-36)
Sendo:
y(t) =
[x(t)
x(t)
](4-37)
O problema de autovalor resultante é escrito da seguinte forma:[C M
M 0
]Yiλi +
[K 0
0 −M
]Yi = 0 (4-38)
O problema da equação (4-38) resulta em 2n autovalores, λi, que se
apresentam em pares conjugados e 2n autovetores:
Λ =
λi
Yi =
[φi
φiλi
](4-39)
96
Sendo os autovalores:
λi = ηi + jωi i = 1, 2, .., n
λi+n = ηi − jωi i = 1, 2, .., n(4-40)
Os modos complexos são dados por: φi = φRi + jφI
i . O sistema pode
ser desacoplado a partir de seus autovetores:
ΦT
[C M
M 0
]Φ = Λ =
µi
,
ΦT
[K 0
0 −M
]Φ = Λ =
µiλi
(4-41)
Sendo Φ a matriz composta pelos autovetores do sistema. Os modos
complexos podem ser usados para representar a resposta dinâmica de um
sistema com amortecimento nãoproporcional, equação (4-42).
x(t) = aeηit[φRi cos(ωit− ϕ)− φI
i sen(ωit− ϕ)] (4-42)
Sendo os coecientes a e ϕ determinados pelas condições iniciais.
Observando a equação (4-42) notase dois coecientes temporais:
c1 = cos(ωit−ϕ), que multiplica o vetor composto pela parte real do modo
complexo φRi ; e c2 = sen(ωit−ϕ), que multiplica o vetor composto pela parte
imaginária do modo complexo φIi . Os modos complexos não apresentam nós
estacionários nem forma bem denida como os modos normais. Pode-se
pensar nos modos complexos como ondas propagantes sem nós xos; eles
não têm forma denida, por isso são de difícil visualização.
Figura 4.4: Tentativa de visualização de um modo complexo
97
4.4Conclusões
Neste capítulo a Análise Modal foi discutida e observouse que, de-
pendendo do sistema dinâmico em análise, podemse obter modos normais
ou modos complexos. No próximo capítulo duas ferramentas essenciais para
a análise de vibrações são discutidas: a Transformada Rápida de Fourier e
a Função Resposta em Freqüência. No capítulo 6 uma técnica alternativa à
Análise Modal, a Decomposição Ortogonal Própria, é discutida.
5Transformada Rápida de Fourier e Função Resposta emFreqüência
5.1Introdução
No capítulo 4 sistemas dinâmicos foram analisados através de equações
diferenciais. A partir das matrizes dos sistemas, um problema de autovalor
é resolvido e as característica fundamentais de um sistema dinâmico são
obtidas: freqüências naturais, modos de vibração e amortecimentos. Além
disso, podese representar a resposta dinâmica a partir dos modos de
vibração de um sistema.
Nem sempre se tem em mãos sistemas dinâmicos modelados em forma
de equações diferenciais. Isto não impede que estruturas sejam testadas e
que suas características sejam obtidas através de testes estruturais. Este
capítulo introduz os conceitos de Transformada Rápida de Fourier e Função
Resposta em Freqüência (FRF), base matemática para testes estruturais.
Este capítulo tem como referências básicas [11], [25], [19] e [10].
5.2Transformada Rápida de Fourier
A resposta de sistemas dinâmicos pode ser estudada tanto no domínio
do tempo quanto no domínio da freqüência.
Dado um certo forçamento, f(t), um sistema mecânico gera uma res-
posta dinâmica, x(t). Este processo pode ser analisado tanto no domínio do
tempo quanto no domínio da freqüência, para sistemas lineares e invariantes
no tempo.
H(ω)F (ω) = X(ω) (5-1)
99
H(ω) é chamada de FRF, ela associa um forçamento F (ω) à resposta
X(ω), em função da freqüência, ω.
A Transformada Rápida de Fourier FFT (Fast Fourier Transform)
é uma ferramenta que transforma uma função no domínio do tempo para
o domínio da freqüência. O algoritmo desta transformada foi criado em
1965 por J.W. Cooly e J.W. Turkey e a partir desta data a aplicação
de técnicas experimentais na dinâmica estrutural se consolidou abrindo
caminho a muitos avanços nesta área.
Para ilustrar a FFT, observe a Figura 5.1. Ela mostra a resposta
dinâmica de um ponto de uma placa; curva cinza. Transformando o sinal
no domínio do tempo para o domínio da freqüência, utilizandose a FFT,
obtémse a curva preta. No domínio da freqüência ca mais fácil visualizar
as freqüência naturais (picos).
Figura 5.1: Freqüências naturais e modos de vibrações
A Figura 5.1 mostra os quatro primeiros modos de vibração de
uma placa. Cada freqüência natural está relaciona com uma conguração
da estrutura. A resposta dinâmica é amplicada quando uma estrutura
é excitada em uma das freqüências naturais. O modo de vibração se
caracteriza como sendo uma conguração na qual todos os pontos vibram
em uníssono.
A FFT tem como origem as Séries de Fourier. Série de Fourier é uma
forma de representar funções periódicas. Pode ser demonstrado matemati-
camente que uma função periódica é a soma de funções senoidais com fre-
qüências múltiplas de uma freqüência fundamental:
x(t) = a0 + 2∞∑
n=1
(ancos
2πnt
T+ bnsin
2πnt
T
)(5-2)
100
Sendo a0 o valor médio de x(t) ao longo de um período T (de −T/2 a
+T/2). E:
an =1
T
∫ +T/2
−T/2
x(t)cos2πnt
Tdt , bn =
1
T
∫ +T/2
−T/2
x(t)sin2πnt
Tdt (5-3)
Os cálculos apresentados são baseados em funções reais e o espectro de
freqüências é denido apenas para freqüências positivas (n=1,2,..). A forma
exponencial da série de Fourier é dada por:
x(t) =+∞∑
n=−∞
1
T
(∫ +T/2
−T/2
x(τ)e−j2πnτ/Tdτ
)ej2πnt/T (5-4)
A parte entre parênteses da equação (5-4) é denida para freqüências
positivas e negativas (n < 0 e n > 0). Quando o período tende ao
innito, ou seja, para funções nãoperiódicas, a série de Fourier se torna
a Transformada de Fourier:
x(t) =
∫ ∞
−∞
(∫ ∞
−∞x(τ)e−j2πfτdτ
)ej2πftdf (5-5)
Ou:
x(t) =
∫ ∞
−∞X(f)ej2πftdf (5-6)
Sendo f = n/T e:
X(f) =
∫ ∞
−∞x(τ)e−j2πfτdτ (5-7)
Em sua grande maioria, aparelhos analisadores de sinais processam
sinais digitalizados, Figura 5.2.
Figura 5.2: Sinal real (esquerda) e sinal digitalizado (direita)
101
A Transformada Discreta de Fourier (DFT) de um sinal no domínio
do tempo é dada por:
X(m) =1
N
N∑n=1
X(n)e−j2π nmN (5-8)
Sendo m o número de linhas de frequência e N o número de amostras
coletadas. Quanto maior o número de linhas melhor a resolução e melhor
a qualidade do resultado da transformada. Por outro lado, o tempo de
processamento do sinal aumenta, aumentando o tempo de coleta. Para uma
freqüência máxima de 1000Hz, com 400 linhas temse uma resolução de
1000/400 = 2, 5 Hz. Já para 2400 linhas, a resolução é de 1000/2400 = 0, 42
Hz. A DFT requer N2 operações complexas, isto faz com que ela seja pouco
eciente quanto ao tempo de processamento.
A Transformada Rápida de Fourier (FFT), usada nos analisadores de
sinais, foi desenvolvida como artifício para computar a DFT de forma
extremamente eciente. Ela utiliza resultados de contas já armazenados e se
aproveita da periodicidade e simetria das funções trigonométricas fazendo
com que o número de operações seja de aproximadamente Nlog2(N). Para
N = 50 a FFT é em torno de 10 vezes mais rápida do que a DFT, para
N = 1000 ela se torna 100 vezes mais eciente e assim por diante.
Freqüência de amostragem: O Teorema de NyquistShannon [25] diz
que a freqüência mínima de amostragem necessária para evitar o mascara-
mento (aliasing) é de duas vezes a freqüência máxima do sinal.
Exemplo 5.1:
O programa t_ex (ver Anexo E) calcula a FFT de um sinal senoidal.
(a) (b)
Figura 5.3: (a) Sinal teórico (verde) e coletado (azul) (b) FFT dos sinais
102
A Figura 5.3a mostra um sinal senoidal em verde e uma curva azul
traçada por pontos amostrados em intervalos determinados. A razão (r)
entre a freqüência do sinal e a taxa de aquisição de dados é de um,
r = 1. Observase que a curva captada não corresponde à curva verdadeira,
tampouco à FFT deste sinal, Figura 5.3b.
(a) (b)
Figura 5.4: (a) Sinal teórico (verde) e coletado (azul) (b) FFT dos sinais
A Figura 5.4 mostra que, para r = 1, 5, ainda não se tem uma boa
representação de um sinal senoidal. A Figura 5.4b mostra que a FFT do
sinal coletado apresenta dois picos, como se houvesse duas freqüências.
(a) (b)
Figura 5.5: (a) Sinal teórico (verde) e coletado (azul) (b) FFT dos sinais
Figura 5.5 mostra que, para r = 2 existe uma boa representação do
sinal, conforme é postulado no Teorema de NyquistShannon. Para r = 50,
Figura 5.6, um sinal pode ser representado com muita conabilidade:
103
(a) (b)
Figura 5.6: (a) Sinal teórico (verde) e coletado (azul) (b) FFT dos sinais
Outra variável importante é o número de ciclos necessários para uma
boa resolução de uma FFT. Para um sinal ser captado em poucos segundos
há a necessidade de replicálo para que se consiga uma FFT aceitável.
Como nem sempre se consegue captar períodos completos acontece o erro
de vazamento (leakage):
(a) (b)
Figura 5.7: Vazamento
Para minimizar o problema do vazamento janelas (windows) são
usadas. A janela mais usada na monitoração da vibração de máquinas
rotativas é a Hanning, por ser mais adequada para sinais aleatórios.
Hanning:
w(t) =1
2
[1 + cos
(2πt
T
)](5-9)
ω = 0 para t < −T/2 e t > T/2.
104
Figura 5.8: Janela Hanning
Outras janelas conhecidas são: retangular, Barlett, Hamming, Black-
man e KaiserBessel.
A FFT é uma das principais ferramentas para a manutenção preditiva
de máquinas rotativas.Através da análise no domínio da freqüência podese
diagnosticar o problema de um equipamento sem precisar parar a produção.
(a) (b)
Figura 5.9: Coleta de sinais (a) Acelerômetro; (b) Sensor de proximidade
O registros de vibrações são coletados periodicamente (ou continua-
mente) através de acelerômetros (ou sensores de proximidade), Figura 5.9.
Os valores são armazenados em computadores, e programas especializados,
Figura 5.10, são usados para a análise.
105
Figura 5.10: Programa usado para auxiliar a análise dos dados de vibrações
Cada pico de vibração no espectro tem uma razão de existir e é através
desta análise que se busca a causa de um possível problema. Algumas
referências sobre diagnósticos de máquinas rotativas são: [2], [4] e [9]. A
seguir dois casos que aparecem rotineiramente:
0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000f [cpm]
-0,5
0,0
0,5
1,0
1,5
2,0
2,5
3,0
3,5
4,0
4,5
5,0
5,5
6,0
6,5
7,0
7,5
8,0
8,5
9,0
9,5
10,0
v [mm/s]M
... TIRAS A QUENTE\13-SISTEMAS VENTILAÇÃO\283-SIST. VENT. SALA MOTORES\015-Ventilador WYD-5\Ventilador\3H -Mancal 3 \Velocidade - Espectro (06/04/2005 15:58:55)
0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000f [cpm]
-0,5
0,0
0,5
1,0
1,5
2,0
2,5
3,0
3,5
4,0
4,5
5,0
5,5
6,0
6,5
7,0
7,5
8,0
8,5
9,0
9,5
10,0
v [mm/s]M
... TIRAS A QUENTE\13-SISTEMAS VENTILAÇÃO\283-SIST. VENT. SALA MOTORES\015-Ventilador WYD-5\Ventilador\3H -Mancal 3 \Velocidade - Espectro (06/04/2005 15:58:55)
1X RPM
Figura 5.11: Desbalanceamento
106
0 50000 100000 150000 200000 250000 300000f [cpm]
-1
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
a [g]M
...Y1-E\Compressor\4A -Mancal 4 \Aceleração - Espectro (06/04/2005 16:00:59)
0 50000 100000 150000 200000 250000 300000f [cpm]
-1
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
a [g]M
...Y1-E\Compressor\4A -Mancal 4 \Aceleração - Espectro (06/04/2005 16:00:59)
# pás X RPM
Figura 5.12: Problema uidodinâmico
A Figura 5.11 mostra o espectro de um equipamento que apresenta
desbalanceamento, que se evidencia pelo pico predominante em 1X (uma
vez a rotação nominal do eixo). A Figura 5.12 mostra o espectro de um
equipamento que apresenta cavitação. Há amplitude de vibração alta na
freqüência múltipla do número de pás do rotor.
5.3Função Resposta em Freqüência
A Função Resposta em Freqüência (FRF) se aplica em sistema lineares
e invariantes no tempo. Ela descreve a relação entre a resposta e a força
de excitação aplicada em um sistema mecânico. A resposta dinâmica pode
ser obtida tanto em deslocamento, como em velocidade ou em aceleração.
Como tanto a força de excitação quanto a resposta dinâmica são grandezas
vetoriais, há direções associadas a elas, portanto, a FRF é uma matriz
composta por diversas funções. Cada função é denida entre um forçamento
em um ponto e uma resposta dinâmica em outro ponto (ou no próprio ponto
de forçamento). Essa relação forçamentoresposta entre dois pontos de uma
estrutura é dada em função da freqüência.
Através da FRF é possível:
identicar e avaliar as vibração mecânicas;
detectar danos e modicações estruturais;
desenvolver modelos dinâmicos a partir de dados experimentais;
qualicar e certicar estruturas e equipamentos;
107
criar critérios e especicações para projetos;
Partindo do sistema:
M x(t) + Cx(t) +Kx(t) = F(t) (5-10)
Admitindo um forçamento harmônico F(t) = Fejωt e uma resposta
harmônica x(t) = Xejωt:
(−ω2M + jωC +K)Xejωt = Fejωt (5-11)
Desta forma o vetor com a amplitude dos deslocamentos X pode ser
encontrado para cada freqüência de excitação harmônica ω, bastando para
isso inverter uma matriz:
X = (−ω2M + jωC +K)−1F = H(w)F (5-12)
A FRF, que dene a relação entre o forçamento e a resposta em
deslocamento, é dada por:
H(ω) = (−ω2M + jωC +K)−1 (5-13)
H(ω) tem entradas complexas. Para um forçamento aplicado no nó j
e resposta obtida no nó i temse um termo da matriz: H(i, j). Cada H(i, j)
é uma FRF.
A Figura 5.13 mostra nove FRFs obtidas por um teste com martelo
de impacto em uma viga em balanço. Em destaque está H(3, 2); resposta
obtida no ponto 3 para um impacto aplicado no ponto 2.
108
Figura 5.13: FRF Teste com martelo de impacto
Um número complexo pode ser escrito em termos de magnitude e fase.
Analisando um elemento da FRF, Hij, perto de uma freqüência natural,
Figura 5.14, observase que existe uma variação de 1800 próximo aos valores
das freqüências naturais.
Figura 5.14: FRF perto de uma freqüência natural (a) amplitude (b) Fase
A parte imaginária de H também é de grande valia, pois é através
109
da parte imaginária que se obtém os modos de vibração de um sistema.
Traçando uma curva a partir dos primeiros picos da parte imaginária de
uma seqüencia de FRFs chegase à forma do primeiro modo de vibração.
Unindose os segundos picos obtémse o segundo modo de vibração e assim
por diante. A Figura 5.15 mostra este processo.
Figura 5.15: FRF Modos de vibração
Nomes das FRFs e suas inversas:
Complacência (Compliance) −→ deslocamento/força
Mobilidade (Mobility) −→ velocidade/força
Receptância (Receptance) −→ aceleração/força
Rigidez Dinâmica (Dynamic Stiness) −→ 1/Complacência
Impedância (Impedance) −→ 1/Mobilidade
Rigidez Inercial (Dynamic Mass) −→ 1/Receptância
110
Figura 5.16: Medição da FRF em uma estrutura
Exemplo 5.2:
Este exemplo encontrase no programa de MATLAB frf.m (ver Anexo
E) . Este programa calcula uma função FRF e os parâmetros modais de
um sistema dinâmico linear. Dadas as matrizes:
M =
4 0 0
0 4 0
0 0 4
, K =
8 −4 0
−4 8 −4
0 −4 4
, Cp = 110M + 1
12K
As freqüências naturais calculadas pelo programa são:
ωn =
1, 8019 0 0
0 1, 2470 0
0 0 0, 4450
As freqüências de oscilação do sistema amortecido calculadas pelo
programa são:
ωd =
1, 7924 0 0
0 1, 2417 0
0 0 0, 4412
Os modos de vibração calculadas pelo programa são:
Φ =
0, 2955 0, 3685 0, 1640
−0, 3685 0, 1640 0, 2955
0, 1640 −0, 2955 0, 3685
111
Os coecientes de amortecimento modais calculadas pelo programa
são:
ξ =
0, 1028 0 0
0 0, 0921 0
0 0 0, 1309
A FRF foi calculada de 0 a 2 rd/s com passo de discretização de 0,01
rd/s. A relação entre força de excitação e resposta dinâmico no primeiro
DOF está expressa na Figura 5.17.
Figura 5.17: FRF do sistema
A curva em verde é a FRF para o sistema sem amortecimento, e
a curva em azul é a FRF para o sistema com amortecimento. Perceba a
diferença na denição dos picos das duas curvas.
A FRF pode ser analisada a partir de outros parâmetros:
H(ω) =1
−mw2 + cjw + k=
R
(jω − p)+
R∗
(jω − p∗)(5-14)
Sendo:
R = −j 1
2mwd
, o resíduo;
wd = wn
√1− ξ2, a freqüência de oscilação do sistema amortecido;
112
p = −σ + jwd, o pólo;
σ =√w2
n − w2d = c/2m, a parte real do pólo;
∗ denota complexo conjugado.
A equação (5-14) representa a expansão em frações da FRF de um
sistema com um grau de liberdade. O parte real do pólo σ é a taxa com que
a oscilação amortecida decai com o tempo. A parte imaginária do polo wd é
a freqüência modal de oscilação do sistema amortecido. O resíduo R é um
número imaginário que está relacionado com a amplitude do modo.
Para sistemas com múltiplos graus de liberdade (MDOF) a equação
(5-14) é escrita da seguinte forma:
Hij(ω) =m∑
r=1
Rijr
(jω − pr)+
R∗ijr(jω − p∗r)
(5-15)
A vantagem de se escrever a FRF desta forma é que o pólo e o resíduo
podem ser obtidos das funções da FRF através de técnicas de ajuste de
curvas como por exemplo o Método dos Mínimos Quadrados. Uma FRF
pode ser obtidas por meio de testes de impactos ou por meio de excitadores
de freqüências. Um excitador é usado, por exemplo, para avaliar as condições
de trabalho de uma peneira vibratória industrial, Figura 5.18.
Figura 5.18: Peneira vibratória industrial
A FRF pode ser obtida por métodos de única entrada e única saída
(SISO - Single Input Single Output) ou por métodos de múltiplas entradas
e múltiplas saídas (MIMO - Multiple Input Multiple Output), ou ainda por
113
combinação destes dois métodos. Estas técnicas não serão discutidas neste
trabalho, mas elas são encontradas, por exemplo, em Maia [25].
Exemplo 5.3:
O programa frf2 (ver Anexo E) calcula a FRF, o pólo e o resíduo de
um sistema com um grau de liberdade:
Usando m = 1 kg e k = 400 N/s a FRF foi calculada para diversos
amortecimentos (c).
(a) (b)
Figura 5.19: (a) FRF (b) Ângulo de fase
A freqüência natural do sistema é 20Hz.
Para c=2: ξ = 0, 05 R=-0,025j e p=-1+19,975j.
Para c=20: ξ = 0, 5 R=-0,029j e p=-10+17,321j.
5.4Conclusões
As ferramentas apresentadas neste capítulo (FFT e FRF) são muito
usadas em indústrias como a petrolífera, siderúrgica, de energia elétrica,
etc. na manutenção preditiva de máquinas rotativas. No próximo capítulo
a Decomposição Ortogonal Própria, uma alternativa à Análise Modal, é
discutida.
6Decomposição de KarhunenLoève
6.1Introdução
O uso da Decomposição de KarhunenLoève (DKL) vêm crescendo
muito nos últimos anos em suas aplicações no campo da Engenharia. O
interesse desse método para este trabalho é o fato dele poder estender a
Análise Modal de forma que a resposta dinâmica de um sistema nãolinear
possa ser escrita como combinação linear de modos empíricos (explicados
mais a frente). A base de KarhunenLoève (KL) é a melhor base de projeção
para o Método de Galerkin.
Devido à relativa novidade do assunto, algumas páginas serão dedi-
cadas a esse tema. O objetivo da DKL é conseguir descrever um fenômeno
de interesse com uma dimensão reduzida que consiga capturar o máximo
desse fenômeno. A DKL é uma forma poderosa e elegante de se obter uma
descrição aproximada de dimensão reduzida de um processo.
A DKL surgiu primeiramente na literatura como PCA (Principal
Component Analysis) e era apenas uma ferramenta para a análise de
sinais. Depois foi estendida para o processamento de imagens e, mais tarde,
para diversas aplicações em engenharia [16]. Na Engenharia Mecânica as
primeiras aplicações foram em escoamentos turbulentos [24].
A DKL pretende obter características dominantes de uma resposta
dinâmica obtida através da análise de dados experimentais ou numéricos. A
DKL é usada na área de Turbulência, processamento de imagens, análise de
sinais, controle na engenharia química, oceanograa, etc. [16]. Esta técnica
também tem aplicação na dinâmica estrutural e em sistemas nãolineares.
O grupo de vibrações da PUCRio vem trabalhando com o assunto há
algum tempo, Sampaio et al [35, 49, 47, 46, 44].
Metodologia da DKL
115
A resposta do sistema é modelada como um processo estocástico de se-
gunda ordem. Entretanto se quer evitar complicações matemáticas com me-
didas de probabilidade associadas à resposta, descrição do espaço amostral
e σalgebra. Uma grande vantagem da DKL é que essas descrições são
desnecessárias, supondo duas condições adicionais: o processo é considerado
estacionário no tempo e ergódico, [27]. Denotando v como sendo o desvio
da resposta com relação à média:
v(x, t) = u(x, t)︸ ︷︷ ︸resposta
−E[u(x, t)]︸ ︷︷ ︸media
(6-1)
Logo, v é um processo estocástico com média zero e, por conseqüência,
seu tensor de correlação é igual ao seu tensor de autocorrelação, [28]. Se v
é real, então a função de autocorrelação espacial de dois pontos é denida
pelo produto tensorial:
R(x, x′) = E[v(x, t)⊗ v(x′, t)] (6-2)
Construção prática da base
Como visto em [35] existem dois métodos de construção das bases
de KL: o método direto e o método dos retratos. O presente trabalho lida
apenas com o método direto.
No método direto, os deslocamentos do sistema dinâmico são medidos
ou calculados em n pontos formando: u1(t), u2(t), ..., un(t). Os deslocamen-
tos são obtidos em m instantes de tempo. Montase a seguinte matriz com
os dados de deslocamentos:
U = [u1 u2 ... un] =
u1(t1) u2(t1) ... un(t1)
. . . .
. . . .
. . . .
u1(tm) u2(tm) ... un(tm)
(6-3)
Usando a hipótese de estacionaridade e ergodicidade, a variação do
116
campo com respeito ao valor médio é:
V = U − 1
m
m∑i=1
u1(ti)m∑
i=1
u2(ti) ...
n∑i=1
un(ti)
. . . .
. . . .
. . . .m∑
i=1
u1(ti)m∑
i=1
u2(ti) ...
m∑i=1
un(ti)
(6-4)
A matriz de correlação espacial é então formada:
R =1
mV TV (6-5)
A matriz R é simétrica. Ela gera autovetores ortogonais, que são os
modos ortogonais próprios (POMs). Os valores próprios (POVs) são dados
pelos autovalores da matriz R.
Exemplo 6.1:
Este problema é extremamente simples e serve para mostrar como uma
resposta dinâmica pode ser reconstruída a partir dos elementos da base de
KL. O programa massamola_kl (ver Anexo E) calcula a resposta dinâmica
de um sistema massamolaamortecedor, Figura 6.1.
Figura 6.1: Sistema massamolaamortecedor com dois graus de liberdade
Sistema: M x+ Cx+Kx = 0
M =
[2 0
0 3
]C =
[0, 4 −0, 1
−0, 1 0, 1
]e K =
[2 −1
−1 1
]
O sistema de equações diferenciais foi integrado pela subrotina ode45
do MATLAB. Dezesseis diferentes condições iniciais foram simuladas para
construir uma base de KL.
117
Primeiro termo da base de KL:
ψ1 =
[0, 5042
0, 8636
]
Segundo termo da base de KL:
ψ2 =
[−0, 8636
0, 5042
]
Energia em cada modo ortogonal próprio (POM):
λ1 = 0, 9054
λ2 = 0, 0946
A resposta do sistema pode ser escrita como:
X(t) = a1(t)
[0, 5042
0, 8636
]+ a2(t)
[−0, 8636
0, 5042
](6-6)
Os coecientes temporais a′is são determinados pela projeção da
dinâmica na base. A Figura 6.2 mostra os coecientes: a1 em azul e a2
em verde.
Figura 6.2: Coecientes temporais a1(azul) e a2(verde)
Condição inicial usada na simulação:
x1 = 1 m; x1 = 0 m/s
x2 = 0 m; x2 = 0 m/s
118
(a) (b)
Figura 6.3: Reconstrução com 1 POM (a) Resposta x1 (b) Resposta x2
A Figura 6.3 mostra a reconstrução da resposta dinâmica com apenas
um elemento da base de KL (curva em azul), neste caso o erro é conside-
rável. A Figura 6.4 mostra a reconstrução da resposta dinâmica com dois
elementos da base de KL, o resultado é preciso, visto que o sistema é nito
e todos os termos da base estão sendo usados.
(a) (b)
Figura 6.4: Reconstrução com 2 POMs (a) Resposta x1 (b) Resposta x2
6.2Base de KL e Modos de Vibração
Muito vem se pesquisando no intuito de se descobrir a relação entre os
modos empíricos (obtidos pela DKL) e os modos de vibração de um sistema
dinâmico.
Os modos empíricos obtidos pela DKL foram usados para discretizar
equações parciais nãolineares pelo Método de Galerkin na área de tur-
bulência, na área estrutural e também em sistemas de identicação, [13].
119
Num sistema dinâmica sem amortecimento: M u(t) +Ku(t) = F (t).
A solução para o sistema homogêneo (F (t) = 0) pode ser escrita como
u(t) = a1(t)v1 + a2(t)v2 + ..+ an(t)vn (6-7)
Sendo ai(t) = bisen(ωi +φi). bi e φi são constantes determinadas pelas
condições iniciais do problema.
ωi e vi são as freqüências naturais e os modos de vibração do sistema
determinados pelo problema de autovalor:
(−ω2M +K)v = 0 (6-8)
A matriz de amostragem V pode ser escrita como
V =
u(t1)
u(t2)
.
.
u(tm)
(6-9)
Sendo ci vetores de dimensão m × 1 compostos pelos coecientes ai
em m instantes de tempo.
A matriz de correlação espacial ca:
R =1
n[c1v
T1 + ..+ cnv
Tn ]T [c1v
T1 + ..+ cnv
Tn ] (6-10)
Sendo n o número de linhas da matriz V .
Os modos de vibração são normalizados em relação à matriz de massa:
vTi Mvi = δij. Supondo M proporcional à matriz identidade, de maneira que
os modos sejam ortonormais, vTi vi = δij, a pós multiplicação da expressão
de R por um dos modos de vibração resulta em
Rvi =1
n[c1v
T1 + ..+ cnv
Tn ]T [c1v
T1 + ..+ cnv
Tn ]vi =
1
n(v1c
T1 ci + ..+ vnc
Tnci)
(6-11)
Como as freqüências dos modos de vibração são distintas, os termos
(vicTi cj)/m tendem a zero quando m→∞, menos o termo (vic
Ti ci)/m, que
será proporcional à vi, [12]. Isto signica que vi se tornou um autovetor
de R. Em outra palavras,para um número de amostras grande os POMs
convergem aos modos normais de vibração. Já os POVs, autovalores de R,
120
convergem para os valores médios quadráticos das respectivas amplitudes
modais.
Feeny e Kappagantu [12] também mostram uma maneira de contornar
a situação em que a matriz de massa M não é proporcional à matriz
identidade. Criase uma matriz de correlação ajustada:
R = RM (6-12)
E o resulta é o seguinte
Rvi = 1n[c1vT
1 + ..+cnvTn ]T [c1vT
1 + ..+cnvTn ]vi = 1
n(v1c
T1 ci + ..+vnc
Tnci)
lim n→∞ Rvi =vic
Ti cin
(6-13)
Portanto os POMs da matriz de correlação ajustada também con-
vergem para os modos normais de vibração do sistema. R não é, neces-
sariamente, uma matriz simétrica, o que implica que seus autovetores nem
sempre serão ortogonais. O termo POM nesses casos pode não ser o mais
adequado.
Quando o sistema tem amortecimento viscoso
u(t) = a1(t)v1 + a2(t)v2 + ..+ an(t)vn (6-14)
Agora com ai(t) = bie−ξωitsen(ωi + φi). Neste caso a eq (36) não é
mais válida, já que quando t → ∞ todos os termos vicTi cj/m tenderão a
zero. Como na prática m é nito, se a estrutura for levemente amortecida,
de forma que vários períodos de oscilação sejam observados, esperase que
os POMs tenham boa convergência com os modos de vibração do sistema.
Quando a estrutura sofrer forçamento harmônico a eq (36) também
não será válida, já que os coecientes ai(t) terão a mesma freqüência, e
nenhum termo vicTi cj/m tenderá a zero.
Uma grande limitação em se usar a DKL na análise experimental
modal é a necessidade de se saber a priori a distribuição de massa do sistema,
ou seja, devese ter conhecimento da matriz de massaM , para que os POMs
e POVs possam ser calculados corretamente.
Podese provar, [48], que a base de KL é a melhor base para repre-
sentar o problema, no sentido de não existir uma outra base que consiga
representar melhor o problema com o mesmo número de termos.
Exemplo 6.2:
121
A resposta dinâmica de um disco girante com molas restringindo as
rotações nos eixos x1 e x2 é obtida. Este exemplo encontrase no programa
de MATLAB disco_kl (ver Anexo E).
X1X2
Figura 6.5: Rotor com dois graus de liberdade, x1 e x2, e velocidade de
rotação constante, Ω
Este é um exemplo simples que tem como objetivo tornar mais fami-
liar o uso da DKL, bem como as conseqüências de sistema com a presença
de forças giroscópicas. Dez diferentes condições iniciais foram utilizadas
para simular o sistema e formar a matriz de correlação R para que os
termos da base de KL fossem calculados.
Sistema: M x+Gx+Kx = 0
M =
[I 0
0 I
]G =
[0 −IpΩIpΩ 0
]e K =
[Kt 0
0 Kt
]
Sendo:
d o diâmetro do disco (a espessura é considerada pequena);
I o momento de inércia;
Ip o momento de inércia polar;
Kt a rigidez em cada direção;
Ω a velocidade de rotação do disco;
Dados da simulação:
d = 1m I =πd4
64m4 Ip =
πd4
32m4 Kt = 1000N.m
As respostas dinâmicas do sistema foram calculadas para a seguinte
condição inicial (esta condição inicial não foi usada na construção da base
de KL):
x1 = 0 rad; x1 = π/90 rad/s
122
x2 = 0 rad; x2 = 0 rad/s
O sistema tem dimensão 2 x 2, logo gera uma base com dois elementos.
Os dois modos empíricos foram calculados através da Decomposição de KL
e usados para a reconstrução de todas as dinâmicas.
Primeiro caso: Ω = 1 rad/s
λ =
[0, 6602
0, 3398
]
ψ =
[0, 5588 −0, 8293
0, 8293 0, 5588
]
(a) (b)
Figura 6.6: Resposta dinâmica (a) x1 (b) x2
A Figura 6.6 mostra que, enquanto a resposta de x1 é amplicada, a
resposta de x2 é reduzida. Existe uma troca de energia entre x1 e x2.
Segundo caso: Ω = 10 rad/s
λ =
[0, 5161
0, 4839
]
ψ =
[0, 2780 0, 9606
0, 9606 −0, 2780
]
123
(a) (b)
Figura 6.7: Resposta dinâmica (a) x1 (b) x2
A Figura 6.7 mostra novamente que existe uma troca de energia entre
as direções, se a resposta em uma direção aumenta, a resposta na outra
direção diminui, e viceversa.
A reconstrução da resposta dinâmica através da base de KL foi precisa,
o que era de se esperar, pois nesse caso o sistema é nito e todos os termos da
base estão sendo usados na reconstrução. Caso apenas um POM seja usado
na reconstrução, a resposta dinâmica não será bem representada, observe a
Figura 6.8.
Figura 6.8: Reconstrução (em azul) com apenas 1 POM
Este sistema tem duas freqüências naturais, w1 e w2. w1 é chamada
de freqüência retrógrada, pois o disco gira no sentido oposto ao sentido de
rotação do rotor. w2 é chamada de freqüência direta, pois o disco gira em
124
concordância com o sentido de rotação do rotor. Se a rotação for zero, existe
apenas uma freqüência natural.
w1 =
√νΩ + ν2Ω2 + 4w2
0
2
w2 =
√−νΩ + ν2Ω2 + 4w2
0
2(6-15)
Sendo ν = Ip/I e w0 = Kt/I. O diagrama Campbell mostra as
freqüências naturais versus a velocidade de rotação do rotor.
Figura 6.9: Diagrama de Campbell para um disco
Caso I > IP , ou seja, o disco esteja mais próximo de um cilindro, o
diagrama de Campbell ca:
125
Figura 6.10: Diagrama de Campbell para um cilindro
6.3Barra chocando-se contra um ostáculo
Um programa chamado barra_choque foi desenvolvido para analisar
este problema. Considere a barra representada na Figura 6.11:
Figura 6.11: Barra chocandose contra um obstáculo
Neste exemplo uma barra é excitada por uma força senoidal na
extremidade livre. O movimento da barra é limitado por um obstáculo
que está a uma distância pequena da extremidade livre. O impacto com
o obstáculo é modelado como uma força de mola, proporcional à rigidez
do obstáculo. O choque entre a barra e o obstáculo introduz uma não
126
Comprimento da barra, L = 1 mDiâmetro da barra, d = 10 cmMódulo de elasticidade, E = 200 GPaDensidade, ρ = 7850 kg/m3
Amortecimento, c = 0, 01 Ns/m2
Rigidez do obstáculo, ka = 100 GN/mForça de excitação, F = 1000 NFreqüência de excitação, ff = 100 HzDistância barraobstáculo, dist = 0, 01 µm
Tabela 6.1: Dados usados no programa barra_choque
linearidade ao sistema. Os dados do problema são:
O Método dos Modos Supostos é usado para discretizar o sistema, a
aproximação é dada por uN(x, t) =N∑
i=1
ai(t)φi(x) . A formulação fraca de
uma barra livre em uma extremidade e livre na outra é dada por:
ai(t)ρA
∫ L
0
φi(x)φj(x)dx+ ai(t)c
∫ L
0
φi(x)φj(x)dx
+ai(t)EA
∫ L
0
dφi(x)
dx
dφj(x)
dxdx =
∫ L
0
PL(t)φj(x)dx
(6-16)
Sendo:
M(φi, φj) = ρA
∫ L
0
φi(x)φj(x)dx , K(φi, φj) = EA
∫ L
0
dφi(x)
dx
dφj(x)
dxdx
C(φi, φj) = c
∫ L
0
φi(x)φj(x)dx F (φj) =
∫ L
0
PL(t)φj(x)dx
(6-17)
A força PL inclui a força de excitação e a força devido ao choque:
PL = Pfsen(ωf t)δ(x− xL) +N∑
i=1
Fci(t)δ(x− xL) (6-18)
Sendo:
Fci(t) = −ξ [ka uLi(t)]ξ = 0 p/ uLi < dist
ξ = 1 p/ uLi > dist(6-19)
O sistema discreto ca:
M a(t) + Ca(t) +Ka(t) = F(t) (6-20)
Este sistema de EDOs é integrado numericamente através da subrotina
127
ode45 do MATLAB.
Primeiramente é feito um teste como o ∆t de integração e o número
de modos (N) para saber os valores mais adequados a serem usados na
simulação:
(a) (b)
Figura 6.12: (a) Convergência da aproximação, variando ∆t; (b) Convergên-
cia da aproximação, variando N
A Figura 6.12 mostra a convergência da aproximação variandose ∆t
e N . As simulações seguintes foram feitas com: ∆t = 0, 0001 e N = 5.
A Figura 6.13 mostra a resposta dinâmica do ponto da extremidade
livre da barra.
Figura 6.13: Resposta do ponto da extremidade livre da barra
A linha vermelha representa a posição do obstáculo. A força média
(uma medida da energia) que atua no obstáculo foi calculada variandose
128
a distância e a rigidez do obstáculo ka. A Figura 6.14 mostra que a força
média varia linearmente com a distância e assintoticamente com a rigidez
do obstáculo para os valores considerados.
(a) (b)
Figura 6.14: (a) Força média x dist; (b) Força média x ka
Base de Karhunen-Loève:
Agora a base de KL (ou modos empíricos ou POMs) é calculada
a partir de diversas simulações. A dinâmica será projetada nessa base.
A Figura 6.15 compara os três primeiros modos de vibração com os três
primeiros POMs.
Figura 6.15: Modos de vibração x base de KL
Um fato curioso é que o segundo modo empírico se assemelha com
o terceiro modo de vibração e o terceiro modo empírico se assemelha
129
com o segundo modo de vibração. Os três primeiros valores ortogonais
próprios (POVs) computados foram: [0, 9984; 0, 0009; 0, 0007]. Estes valores
representam o quanto cada POM contribui para a resposta dinâmica.
Observase que quase toda a energia está no primeiro modo empírico.
A dinâmica pode ser reconstruída com os elementos da base de KL.
Sirovich, [40], recomenda que 99% da energia deve ser considerada para
representar o problema, o primeiro POM contém mais de 99% da energia.
A Figura 6.16 mostra uma comparação da resposta dinâmica em x = L,
aproximando com cinco modos normais, com três modos normais e com um
modo empírico.
Figura 6.16: Aproximação da resposta dinâmica em x = L
A Figura 6.17 mostra o detalhe do choque. É possível notar que a
resposta para uma aproximação com cinco modos normais (curva azul) é
praticamente coincidente com a aproximação com um modo empírico (curva
verde). Já a aproximação com três modos normais (curva preta) apresenta
nítida diferença com as outras duas aproximações.
130
Figura 6.17: Aproximação da resposta dinâmica perto da região de choque
Os modos empíricos calculados pela DKL conseguem capturar melhor
a dinâmica do que os modos normais. Isto se deve ao fato da DKL levar
em consideração a não linearidade do sistema. Este exemplo ilustra como a
base de KL é a melhor base de projeção da dinâmica.
Apesar do primeiro modo empírico ser muito semelhante ao primeiro
modo de vibração, Figura 6.18:
Figura 6.18: Primeiro modo normal x Primeiro modo empírico
Para o cálculo da resposta dinâmica, a segunda derivada do modo
entra no cálculo da matriz de rigidez (K) do problema. Comparando a
131
segunda derivada do primeiro modo normal com a segunda derivada do
primeiro modo empírico, Figura 6.19, percebese uma nítida diferença entre
as funções:
Figura 6.19: Derivada segunda do primeiro modo normal x Derivada segunda
do primeiro modo empírico
6.4Conclusões
Neste capítulo a Análise Modal foi estendida para sistemas não
lineares: os elementos da base de KarhunenLoève geram um espaço de
Hilbert, e esta base é a melhor base de projeção para dinâmica. No capítulo
seguinte é analisado um caso de um sistema rotormancal modelado de
forma similar ao exemplo da barra chocandose contra um obstáculo.
Entretanto o sistema é muito mais complexo.
7Inuência do desbalanceamento e folga nos mancais deum rotor em balanço
Neste capítulo um sistema rotormancal é modelado conforme as
teorias descritas nos capítulos anteriores. Este capítulo gerou um artigo que
foi selecionado para ser publicado no Congresso Internacional de Engenharia
Mecânica (COBEM 2005), [32], Novembro de 2005, Ouro Pedro, MG.
Máquinas rotativas são muito importantes no processo produtivo e a cada
dia os processos demandam que equipamentos operem por períodos mais
longos e com cargas e rotações maiores. O objetivo é investigar a inuência
do desbalanceamento e folgas nos mancais de um exaustor, modelado como
um rotor em balanço. Os mancais são os componentes que suportam toda
energia da carga e de eventuais impactos. Os choques introduzem não
linearidades no modelo. O sistema rotormancal é modelado como um
sistema contínuo, na sua formulação fraca. O Métodos de Galerkin é usado
para discretizar o sistema. As forças nos mancais e a Transformada Rápida
de Fourier (Fast Fourier Transform FFT) da resposta transiente são
computados. Os valores obtidos por simulação numérica são comparados,
qualitativamente, com casos reais.
7.1Introdução
Máquinas rotativas são muito utilizadas em diversos tipos de indús-
trias como: geração de energia elétrica, siderurgia, petróleo e gás, papel e
celulose, etc. A modelagem de rotores é essencial para a análise dinâmica e
controle de vibração de sistemas como exaustores, motores, bombas, com-
pressores, geradores, etc.
O desbalanceamento de componentes girantes é a causa mais comum
de vibração alta em máquinas rotativas. Há muitas normas sobre balancea-
mento como por exemplo: ISO 2953 Mechanical Vibration Balancing
machines Description and evaluation (Vibração Mecânica Balancea-
133
mento de máquinas Descrição e avaliação); ISO 19401 Mechanical
Vibration Balance quality requirements for rotors in a constant (rigid)
state (Vibração Mecânica Qualidade do balanceamento requerida para
motores em estado constante (rígido)); ISO 11342 Methods and criteria
for the mechanical balancing of exible rotors (Métodos e critérios para o
balanceamento mecânico de rotores exíveis).
Desbalanceamento, desalinhamento e folga são responsáveis por quase
90% dos problemas de vibração de máquinas rotativas. Serão tratados os
casos de desbalanceamento e folga. Os mancais são os componentes que
suportam essa vibração. A norma ISO 10816 Mechanical vibration Eval-
uation of machine vibration by measurements on nonrotating parts (Vi-
bração mecânica Avaliação da vibração de máquinas através da medição
nas partes não girantes) tem referências de níveis de vibração bons e
aceitáveis medidos nos mancais de diversas máquinas rotativas.
Childs (1993) [7] é uma das referências padrão na modelagem e análise
rotodinâmica. Os modelos matemáticos de rotores são complexos mesmo
para casos lineares. Qiu e Rao (2005) [30] usam uma abordagem diferente,
com lógica nebulosa, para o problema. Os parâmetros de um sistema
rotormancal apresentam variações e incertezas consideráveis devido à
fabricação, montagem e condições operacionais. Será tratado apenas o caso
determinístico neste trabalho.
A folga e o desbalanceamento são investigados por muitos artigos
recentes. A questão da folga interna radial de mancais de rolamento é
investigada por Harsha (2005a e 2005b) [14, 15], Sinou e Thouverez (2004)
[39] e Tiwari et al. (2000) [43]. Todos estes artigos consideram as não
linearidades causadas pelo contato de Hertz e desenvolvem modelos através
de sistemas massamola. Karlberg e Aidanpää (2003 and 2004) [20], [21]
e Vakakis e Azeez (1999) [45] consideram folga entre rolamento e caixa
de rolamento e folga entre eixo e rolamento, respectivamente. Todos os
artigos citados centralizam suas discussões em movimentos caóticos, mapas
de Poincaré, diagramas de bifurcação e expoentes de Lyapunov. Neste
trabalho, por outro lado, o domínio da freqüência é usado para caracterizar
alguns problemas. Na manutenção preditiva a ferramenta mais utilizada
para identicar causas de problemas em máquinas rotativas através da
análise de vibrações é a FFT.
A grande limitação de modelos simples massamola é a falta de um
esquema para computar o erro do modelo, além da diculdade para a
escolha dos parâmetros. Esses modelos certamente ajudam na compreensão
do fenômeno, mas não podem descrevêlos quantitativamente já que não há
134
um esquema de aproximação. Modelos contínuos não têm esta limitação;
quando a discretização é feita pelo Método de Galerkin as propriedades que
aparecem são as do material e há um esquema de aproximação intrínseco.
A análise de rotores por elementos nitos normalmente resulta em
problemas com dimensão muito grande, o que torna o cálculo da resposta
transiente muito lenta. Esta discretização inclui muitos modos insigni-
cantes por causa da extensão do autoespectro. Vakakis e Azeez (1999) [45]
usam o Método dos Modos Supostos e a Decomposição de KarhunenLoève
(KLD) para reduzir o sistema. KLD é uma técnica relativamente recente
na análise de vibrações que usa estatística para formar uma base de pro-
jeção da dinâmica que contém a maior parte da coerência. KLD funciona
bem tanto para sistemas lineares como para sistemas nãolineares, Wolter
e Sampaio (2001) [46].
Folgas danosas entre rolamento e caixa de mancal acontecem caso haja
desgaste das partes ou problemas de fabricação ou na montagem. Em muitas
máquinas rotativas, devido a expansão térmica do eixo, um rolamento deve
ser livre para deslocar na direção axial, portando naturalmente haverá
alguma folga. Os rolamentos são os componentes que absorvem as energias
dos impactos e da carga, tendo, como conseqüência, sua vida reduzida.
Um modelo contínuo de um rotor em balanço com folga é simulado. A
caixa de mancal é representada por molas e amortecedores simetricamente
distribuídos pelo perímetro da seção transversal do eixo. Uma discretização
inicial é feita pelo Método dos Elementos Finitos para computar os modos
de vibração do sistema. Depois o Método dos Modos Supostos é usado na
simulação do problema.
O objetivo é calcular a FFT da resposta dinâmica em determinados
locais e usar esta informação para avaliar o estado da máquina.
A motivação para a construção deste modelo é correlacionar e melhor
entender a dinâmica observada em diversas máquinas rotativas. Observase
que o desbalanceamento e a folga afetam as forças e choques que atuam nos
rolamentos. A FFT do sinal no tempo capturado nos mancais pode ser bem
distinta para diferentes folgas e desbalanceamentos. A análise de vibração
na manutenção preditiva de máquinas rotativas é feita com base na FFT
observada, portanto é de vital importância numa indústria entender como a
FFT descreve o estado das máquinas. Grau de qualidade de balanceamento,
folga, e rotação da máquina são usados como parâmetros.
135
7.2Motivação
Dois casos representativos são discutidos, um recente e um antigo. Os
dois aconteceram na CST (Companhia Siderúrgica de Tubarão), Espírito
Santo, Brasil. Estes casos servem de motivação e guia para o desenvolvi-
mento do modelo que pretende melhorar o entendimento de sistemas roto-
dinâmicos.
Qualquer carta de vibração (veja Bloch e Geitner 1983 [4]) mostra
que quando há desbalanceamento, o espectro de freqüências apresenta um
pico na freqüência de rotação da máquina (X). Quando há folgas diversas
(entre rolamento e caixa de mancal, por exemplo), o espectro de freqüências
apresenta diversos picos em freqüências múltiplas da freqüência de rotação
da máquina (2X, 3X,...).
Para entender melhor o fenômeno em estudo observe a Figura 7.1. A
gura mostra um corte transversal do mancal.
Caixa do mancal
Rolamento
Eixo
Folga
Figura 7.1: Corte transversal do mancal
136
7.2.1Exaustor Principal da Máquina de Lingotamento Contínuo rotaçãonominal 1185 RPM
Figura 7.2: Exaustor Principal da Máquina de Lingotamento Contínuo
A medição de vibração foi feita com um acelerômetro na posição 4H
(horizontal), Figura 7.2. A velocidade é obtida pela integração do sinal
coletado. A Figura 7.3 mostra o espectro de freqüências.
harmônicos
CPM CPM
V(mm/s)
1XRPM=4,6mm/s30/01/2005 09/02/2005
1XRPM=1,7mm/s
(a) (b)
Figura 7.3: Espectro de freqüências. (a) 30 de janeiro de 2005, (b) 9 de
fevereiro de 2005
Um pedaço de metal agarrou em uma pá do rotor no dia 30 de janeiro
de 2005 aumentando assim o desbalanceamento. A inuência da folga entre
rolamento e caixa de mancal está registrada no espectro de freqüências,
Figura 7.3a. É possível ver diversos harmônicos (2X, 3X,..).
Após a remoção do pedaço de metal da pá e a realização de um bal-
anceamento de campo, os harmônicos praticamente desapareceram, Figura
7.3b. O desbalanceamento residual cou em 1,7 mm/s 0Peak.
137
7.2.2Exaustor da dessulfuração rotação nominal 1800 RPM
Figura 7.4: Exaustor da dessulfuração
A medição de vibração foi feita com um acelerômetro na posição 3H
(horizontal), Figura 7.4. A Figura 7.5 mostra o espectro de freqüências.
(a) (b)
Figura 7.5: Espectro de freqüências. (a) 2 de julho de 2001, (b) 23 de maio
de 2002
A folga neste caso era tamanha que mesmo um pequeno desbalancea-
mento provocava choques. Era necessária, quase que mensalmente, a rea-
lização de balanceamento de campo, uma operação dispendiosa. A Figura
7.5a mostra o espectro de freqüências para um pequeno desbalanceamento
praticamente não se observa harmônicos. A Figura 7.5b mostra o espec-
tro de freqüências para um desbalanceamento maior é possível observar
diversos harmônicos (2X, 3X,..).
A especicação do rolamento teve que ser mudada, diminuindo assim
a folga e resolvendo o problema.
138
7.2.3Desbalanceamento x Vida do rolamento
Para caracterizar o desbalanceamento de um rotor será usado o grau
de qualidade de balanceamento (G), um índice obtido da experiência acu-
mulada de um grande número de rotores diferentes (consulte ISO 1940/1).
Geralmente o desbalanceamento residual permissível (Uper [kg ×m])
para rotores é proporcional à massa do rotor (Uper ∼Mrotor). O desbalacea-
mento residual especíco permissível é denido como: eper = Uper/Mrotor. A
experiência mostra que o desbalanceamento residual especíco permissível
varia inversamente com a rotação do rotor. (Uper ∼ 1/Ω). Finalmente:
G = eperΩ1000 [mm/s] (7-1)
Sendo Ω a velocidade de rotação em rad/s e eper=[kg ×m/kg].
Nas análises numéricas realizadas nas próximas seções G é um
parâmetro escolhido e o desbalanceamento é calculado:
U =GMrotor
Ω1000[kg ×m] (7-2)
Uma massa desbalanceada pode produzir uma força síncrona excessiva
que reduz a vida de vários elementos mecânicos.
Por exemplo, considere uma máquina com grau de qualidade de
balanceamento G = 6, 3 e rotação 2000 RPM , com rolamento de esferas
(SKF 6917). A vida do rolamento pode ser calculada com ou sem a carga
extra devido ao desbalanceamento.
Considerando uma carga de 1757 Newtons a vida do rolamento será
de 21,4 anos. A força centrífuga devido ao desbalanceamento neste caso é de
132 Newtons. Esta carga adicional devido ao desbalanceamento reduz a vida
do rolamento em 27% (15,6 anos). Este cálculo não leva em consideração
folgas excessivas e choques. Dependendo das características do equipamento
e da rotação a redução da vida pode ser maior ou menor.
7.3Equacionamento
Um programa chamado rotor_choque foi desenvolvido para analisar o
problema. A Figura 7.6 mostra o esquema de um sistema contínuo de um
rotor em balanço.
139
Figura 7.6: Sistema rotodinâmico considerado
A duas direções estão acopladas devido às forças giroscópicas. Supõe
se que a inércia rotatória tem pouca inuência na dinâmica do sistema
considerado.
O modelo de viga EulerBernoulli é usado. As forças externas que
atuam no sistema são: forças centrífugas devido ao desbalanceamento; forças
devido à massa do rotor e à gravidade; e forças devido aos impactos, que
introduzem as não linearidades ao sistema.
As forças devido ao desbalanceamento (Q1,2) e aos impactos (P1,2) são
escritas como:
Q1 = UΩ2cos(Ωt) , Q2 = UΩ2sen(Ωt)−Mrg
P1 = −ξ(kh
(r − c)u1
r+ dh
(u1u1 + u2u2)u1
r2
)P2 = −ξ
(kh
(r − c)u2
r+ dh
(u1u1 + u2u2)u2
r2
) (7-3)
Sendo:
c = folga , ξ = 1 se r ≥ c , e ξ = 0 se r < c
r =√u2
1 + u22 , g = aceleração da gravidade
dh = amortecimento do batente, kh = rigidez do batente
U = desbalanceamento
Os impactos só ocorrem para deslocamentos que ultrapassam o ba-
tente.
O sistema é discretizado pelo Método dos Elementos Finitos para que
sejam computados os modos de vibração, dada uma determinada precisão,
Figura 7.7.
140
u1 → deslocamento na direção e1u2 → deslocamento na direção e2φ → modos de vibraçãoE → módulo de elasticidadeI → momento de inércia do eixom → densidade vezes área da seção (ρA)Mr → massa do rotorIr → momento de inércia do rotorL → comprimento do eixoΩ → velocidade de rotaçãod1 e d2 → coecientes de amortecimentoxB3 e xB4 → posicionamento dos mancais
Tabela 7.1: Parâmetros do programa rotor_choque
Figura 7.7: Viga engastada em uma extremidade e com uma massa na outra,
discretizada pelo MEF
A precisão requerida foi de 0, 0005% de erro no quadragésimo modo.
Para satisfazer esta precisão foram necessários 450 elementos e aproximada-
mente 7 minutos de simulação. O sistema tem então sua dinâmica projetada
nos modos de vibração calculados. A formulação fraca, levando em consid-
eração as forças externas e as condições de contorno pode ser escrita como:
∫ L
0
(EI
∂4u1
∂x4+m
∂2u1
∂t2+ 2ΩIr
∂2u2
∂x∂tδ′(x− L) + d1
∂u1
∂t
)φjdx
=
∫ L
0
(P1(t)δ(x− xB3)) +
∫ L
0
(P1(t)δ(x− xB4) +Q1(t)δ(x− L))φjdx∫ L
0
(EI
∂4u2
∂x4+m
∂2u2
∂t2− 2ΩIr
∂2u1
∂x∂tδ′(x− L) + d2
∂u2
∂t
)φjdx
=
∫ L
0
(P2(t)δ(x− xB3))
∫ L
0
(P2(t)δ(x− xB4) +Q2(t)δ(x− L))φjdx
(7-4)
A Tabela 7.1 mostra o signicado de cada parâmetro da equação (7-4).
141
A aproximação é dada por:
u1 =N∑
i=1
aiφi ; u2 =N∑
i=1
biφi (7-5)
Sendo N o número de modos usados na aproximação. As funções da
base φ são os autovetores de uma viga engastada com massa na extremidade
calculadas pelo MEF.
Após a integração por partes, usando as condições de contorno e
substituindo u1 e u2 por suas aproximações chegase ao seguinte sistema
de equações:
ai
∫ L
0
mφiφjdx+ ai
∫ L
0
EIφ′′i φ′′jdx+ bi2ΩIrφ
′′i (L)φ′j(L)
+ai
∫ L
0
d1φiφjdx = P1(t)φj(xB3) + P1(t)φj(xB4) +Q1(t)φj(L)
bi
∫ L
0
mφiφjdx+ bi
∫ L
0
EIφ′′i φ′′jdx− ai2ΩIrφ
′′i (L)φ′j(L)
+bi
∫ L
0
d2φiφjdx = P2(t)φj(xB3) + P2(t)φj(xB4) +Q2(t)φj(L)
(7-6)
A formulação fraca é, através do Método de Galerkin, discretizada
resultando em um sistema de equações diferenciais ordinárias.
Os operadores obtidos da primeira equação são os seguintes:
M1(φi, φj) = m
∫ L
0
φi(x)φj(x)dx
C1(φi, φj) = d1
∫ L
0
φi(x)φj(x)dx
G1(φi, φj) = 2ΩIrφ′′i (L)φ′j(L)
K1(φi, φj) = EI
∫ L
0
φ′′i (x)φ′′j (x)dx
F1(φj) = P1(t)φj(xB3) + P1(t)φj(xB4) +Q1(t)φj(L)
(7-7)
142
Os operadores obtidos da segunda equação são os seguintes:
M2(φi, φj) = m
∫ L
0
φi(x)φj(x)dx
C2(φi, φj) = d2
∫ L
0
φi(x)φj(x)dx
G2(φi, φj) = 2ΩIrφ′′i (L)φ′j(L)
K2(φi, φj) = EI
∫ L
0
φ′′i (x)φ′′j (x)dx
F2(φj) = P2(t)φj(xB3) + P2(t)φj(xB4) +Q2(t)φj(L)
(7-8)
Os operadores do sistema são os seguintes:
M(φi, φj) =
[M1(φi, φj) 0
0 M2(φi, φj)
]
C(φi, φj) =
[C1(φi, φj) 0
0 C2(φi, φj)
]
G(φi, φj) =
[0 G1(φi, φj)
−G2(φi, φj) 0
]
K(φi, φj) =
[K1(φi, φj) 0
0 K2(φi, φj)
]
F (φj) =
[F1(φj)
F2(φj)
], X(t) =
[ai(t)
bi(t)
]
(7-9)
O sistema discreto ca:
MX(t) + (C +G)X(t) +KX(t) = F(x, t) (7-10)
7.4Simulação numérica
O sistema de equações diferenciais parciais obtido na seção anterior é
integrado através da subrotina do MATLAB ode45.
Os parâmetros usados na simulação padrão são dados na Tabela 7.2:
143
Comprimento do eixo, L = 3053 mmDiâmetro do eixo, Ds = 110 mmMassa do rotor, Mr = 600 kgMódulo de elasticidade, E = 193 GPaDensidade, ρ = 8000 kg/m3
Amortecimento, d1 = d2 = 1000 Ns/m2
Grau de qualidade de balanceamento, G = 6 mm/sfolga, c = 7, 5 µmrigidez do batente, kh = 1 GN/mPosicionamento do mancal 3, xB3 = 1, 692 mPosicionamento do mancal 4, xB4 = 2, 302 mLargura do mancal, Bw = 58 mmVelocidade de rotação, Ω = 124 rd/s = 1185 RPM
Tabela 7.2: Dados usados no programa rotor_choque
As dimensões e propriedades dos materiais são as do exaustor da
seção 7.2.1, exceto pela rigidez do batente e pelos amortecimentos. O
amortecimento do sistema foi escolhido propositadamente grande devido,
infelizmente, a falta de disponibilidade de computadores mais poderosos.
O objetivo aqui é fazer uma investigação qualitativa de como as forças
que atuam nos mancais e a curva FFT se comportam com a mudança de
parâmetros como desbalanceamento e folga.
(a) (b)
Figura 7.8: (a) Convergência da aproximação, variando dt; (b) Convergência
da aproximação, variando N
A massa do rotor foi incluída na simulação como uma força agindo na
direção de u2 na extremidade do eixo.
Para selecionar o passo de integração (dt) e o número de modos da
aproximação (N) que serão usados na geração das dinâmicas, o erro da
144
dinâmica foi calculado para diferentes passos de integração e número de
modos.
A Figura 7.8 mostra o erro percentual da dinâmica calculada no
mancal 3 para diferentes passos de integração e número de modos da
aproximação. O erro diminui com o aumento do número de modos e com
a diminuição do passo de integração. O ponto ótimo encontrado para as
simulação realizadas é: t = 0, 001 e N = 10, desta forma todas as
simulações foram feitas com esses parâmetros xados em seus melhores
valores.
A Figura 7.9 mostra a forma do eixo quando a máquina está em
operação. A zona de carga do mancal número quatro é na região inferior
da caixa de mancal enquanto que a zona de carga do mancal número três é
região superior.
Figura 7.9: Conguração do eixo quando em operação
(a) (b)
Figura 7.10: Resposta dinâmica no mancal número 3. (a) Resposta tran-
siente u1 (azul) e u2 (verde); (b) Órbita após transiente
A Figura 7.10 mostra a resposta dinâmica e órbita no mancal número
três. Os valores da simulação só são registrados após a estabilização do
sistema.
145
Simulações foram feitas para diferentes rigidezes do batente. A Figura
7.11 mostra que a força média no mancal 4 aumenta como o aumento da
rigidez. Entretanto, existe um ponto (kh = 10 GN/m para esta simulação)
acima do qual as forças se mantêm praticamente constante. Já a força média
no mancal 3 tem o comportamento inverso, ela diminui até um certo valor.
Figura 7.11: Força média no mancal número 3 (azul) e 4 (verde)
Simulações foram feitas variando o amortecimento do batente de 1
a 1000 Ns/m2. Os valores foram os mesmos para os diferentes valores de
amortecimento (erros menores que 0,01%). O amortecimento do batente não
tem inuência na dinâmica para este modelo.
146
Figura 7.12: Força média no mancal número 4
Para diferentes velocidades de rotação a Figura 7.12 mostra que a
força média no mancal número quatro aumenta com o aumento da rotação.
A relação é quase quadrática, o que é compatível com a relação Fc = UΩ2.
O foco da investigação é a inuência do desbalanceamento e folga na
resposta do sistema.
(a) (b)
Figura 7.13: (a) Força média no mancal número 3 versus G; (b) Força média
no mancal número 4 versus G
Variando o grau de qualidade do desbalanceamento é possível notar
que a força média no mancal número quatro aumenta com o aumento do
desbalanceamento, Figura 7.13b. A relação é praticamente linear, o que é
compatível com a relação Fc = UΩ2. Por outro lado, a força média no
mancal número três não segue este padrão, 7.13a. A Figura 7.13 mostra que
147
a força média no mancal número quatro é muito maior do que a força média
no mancal número três. Isto ocorre devido à massa do rotor que está em
balanço.
Figura 7.14: FFT do deslocamento no mancal número 3: G = 2 (azul) e
G = 20 (magenta)
A Figura 7.14 mostra duas curvas FFT, uma para G = 2 e outra
para G = 20 [mm/s]. É possível notar picos nas freqüências múltiplas da
freqüência fundamental (2X, 3X,...) no espectro de freqüências para grandes
desbalanceamentos. Se o desbalanceamento é pequeno, G = 2 por exemplo,
quase não é possível notar os harmônicos. Já se o desbalanceamento é
grande, G = 20 por exemplo, os harmônicos são evidentes.
Podese concluir que o desbalanceamento inui no estado do sistema.
Esta é o resultado mais importante pois ele é o mesmo observado nas
máquinas rotativas discutidas nas seções 7.2.1 e 7.2.2.
A medição dos sinais de vibração nos mancais do Exaustor Princi-
pal da Máquina de Lingotamento Contínuo (seção 7.2.1) através de um
acelerômetro são qualitativamente os mesmos dos computados na simu-
lação numérica, Figura 7.14. Outra concordância importante entre os dados
numéricos e experimentais diz respeito à forma do eixo quando a máquina
está em operação, Figura 7.9. De fato, a diferença de fase medida entre os
mancais número três e quatro do exaustor é de ∼ 1800.
148
(a) (b)
Figura 7.15: (a) Força média no mancal número 3 versus folga; (b) Força
média no mancal número 3 (azul) e 4 (verde) versus folga
Outro parâmetro importante é o valor da folga. A Figura 7.15 mostra
a força média nos mancais, variando o valor da folga.
A força média no mancal número quatro tornase menor com o au-
mento da folga. Isto acontece devido ao engaste da viga. A força elástica do
eixo contra o movimento do eixo aumenta com o aumento da folga. Entre-
tanto, a força média no mancal número três apresenta um comportamento
curioso, ela passa por um mínimo. Este é o mesmo fenômeno detectado por
Sampaio e Soize (2005) [36].
7.5Conclusões
Um modelo de um sistema rotormancal foi proposto neste capítulo. A
simulação numérica mostrou comportamento, qualitativamente, semelhante
ao comportamento observado em máquinas rotativas. O modelo pode então
ajudar na compreensão do comportamento de rotores reais.
8Conclusões e Trabalhos Futuros
Este trabalho tratou a análise de vibrações no contexto da formulação
fraca. Um sistema mecânico contínuo é equacionado na formulação fraca
em um espaço de Hilbert. Uma base de projeção é criteriosamente escolhida
para a dinâmica, e o sistema é aproximado / discretizado pelo Método de
Galerkin, base para métodos como o MEF e o Método dos Modos Supostos.
A Figura 8.1 mostra esse esquema.
Sistema Mecânico Contínuo
Formulação Fraca
Escolha da base de projeção
Método de Galerkin
Aproximação / Sistema Discreto
Erro da aproximação
Figura 8.1: Etapas da análise
A aproximação é descrita por elementos de um espaço de Hilbert, o
que garante um esquema de aproximação e convergência, a medida que mais
elementos do espaço são usados na aproximação, Figura 8.2.
150
Sistema Discreto
Resposta
PrecisãoSatisfeita?
Aumentar número de elementosda base usados na aproximação
ok!
Não
Sim
Figura 8.2: Algoritmo para convergência
A melhor base de projeção para um sistema dinâmico linear é a
formada pelos modos normais de vibração, que podem ser obtidos pela
Análise Modal. Para um sistema nãolinear a base de KarhunenLoève
(formada pelos modos empíricos) é a melhor base de projeção, e ela pode
ser obtida pela DKL, Figura 8.3.
Figura 8.3: Análise Modal e Decomposição de KarhunenLoève
Este texto está sendo melhorado para se tornar uma apostila didática:
Mais exemplos estão sendo preparados;
Outros tipos de elementos estão sendo testados;
Outros tipos de carregamentos e congurações estão sendo considera-
dos;
151
As simulações do comportamento dinâmico de um rotor, apresentadas
no capítulo 7, corroboraram com as observações obtidas em máquinas ro-
tativas. Este capítulo gerou um artigo que foi selecionado para o 18th In-
ternational Congress of Mechanical Engineering (COBEM 2005), [32]. Esta
concordância motiva sosticar este modelo em busca de comparações quan-
titativas. Dentro do assunto rotor, outros aspectos podem ser considerados
em trabalhos futuros:
Melhorar a descrição do amortecimento;
Inuência do lubricante;
Controle da dinâmica;
Alguns programas foram desenvolvidos no transcorrer do trabalho. A
lista dos programas desenvolvidos se encontra no anexo C. A relação entre os
programas se encontra no anexo D. E o manual dos programas se encontra
no anexo E.
Bibliograa
[1] LEI, Y.; FRISWELL, M. I. ; ADHIKARI, S.. A Galerkin method for
distributed systems with nonlocal damping. IJSS Manuscript, 05,2005.
[2] NEVADA, B.. Machinery diagnosis course. Bently Nevada TechnicalTraining, 2003.
[3] BLEVINS, R. D.. Formulas for Natural Frequency and Mode
Shape. Krieger Publishing Company, 1993.
[4] BLOCK, H. P.; GEITNER, F. K.. Machine Failure Analysis and
Troubleshooting - Pratical Management Process Plants - Vo-
lume 2. Gulf Publish Company, 1983.
[5] AUSTRELL, P. E.; DAHLBLOW, O.; LINDEMANN, J.; OLSSON,A.; OLSSON, K. G.; PERSSON, K.; PETERSSON, H.; RISTINMAA,M.; SANDBERG, G. ; WERNBERG, P. A.. CALFEM - A Fi-
nite Element Toolbox, version 3.4, 2004. Disponível emhttp://www.byggmek.lth.se/Calfem.
[6] CHAPRA, S. C.; CANALE, R. P.. Numerical Methods in Engineers.Mcgraw-Hill International editions, 1990.
[7] CHILDS, D.. Turbomachinery Rotordynamics: Phenomena,
Modeling, and Analysis. Wiley-Interscience, 1993.
[8] CIARLET, P. G.. The Finite Element Method for Elliptic Pro-
blems. North-Holland, 1990.
[9] DIMAROGONAS, A.. Vibration for Engineers. PrenticeHall, 1995.
[10] DOSSING, O.. Structural testing part 2 - Modal Analysis and
simulation. Bruel Kjaer - Technical Review, 1988.
[11] EWINS, D. J.. Modal Testing: Theory and Practice. ResearchStudies Press Ltd., 1984.
153
[12] FEENY, B. F.; KAPPAGANTU, R.. On the physical interpretation
of Proper Orthogonal Modes in vibrations. Journal of Sound andVibration, 211(4):607616, 1998.
[13] FEENY, B. F.; LIANG, Y.. Interpreting Proper Orthogonal Modes
of randomly excited vibration systems. Journal of Sound andVibration, 265:953966, 2003.
[14] HARSHA, S. P.. Nonlinear dynamic response of a balanced rotor
supported on rolling element bearing. Mechanical System andSignal Processing, 19:551578, 2005.
[15] HARSHA, S. P.. Nonlinear dynamic analysis of an unbalanced
rotor supported by roller bearing. Chaos Solutions and Fractals,26:4766, 2005.
[16] HOLMES, P.; LUMLEY, J. L. ; BERKOOZ, G.. Turbulence, coherentstructures, dynamical systems and symmetry. Cambridge Uni-versity Press, 1996.
[17] HUGHES, T. J. R.. The Finite Element Method - Linear Static
and Dynamic Finite Element Analysis. Prentice-Hall, Inc., Engle-wood Cli, New Jersey, 1997.
[18] INMAN, D. J.. Engineering Vibration. Prentice-Hall, Inc., UpperSaddle River, N.J. USA, 1996.
[19] HE, J.; FU, Z.-F.. Modal Analysis. Butterworth-Heinemann, Woburn,M.A. USA, 2001.
[20] KARLBERG, M.; AIDANPÄÄ, J.. Numerical investigation of un-
balanced rotor system with bearing clearance. Chaos Solutionsand Fractals, 18:653664, 2003.
[21] KARLBERG, M.; AIDANPÄÄ, J. O.. Investigation of an unbalanced
rotor system with bearing clearance and stabilising rods. ChaosSolutions and Fractals, 20:363374, 2004.
[22] AL-KHAFAJI, A. W.. Numerical Methods in Engineering Prac-
tice. Holt, Rinehart and Winston, Inc., 1986.
[23] KREYSZIG, G.. Introductory functional analysis with applica-
tions. John Wiley and Sons, 1978.
154
[24] LUMLEY, J. L.. Stochastic tools in turbulence. Academic Press,1970.
[25] MAIA, N.; SILVA, J.. Theoretical and Experimental Modal Ana-
lysis. Research Studies Press LTD, Tauton, Somerset England, 1997.
[26] MEIROVITCH, L.. Principles and Techniques of Vibrations.Prentice-Hall, Inc., Upper Saddle River, New Jersey, 1997.
[27] NEWMAN, A.. Model reduction via Karhunen-Loève expansion
part I: an exposition. Institute for Systems Research, 1996.
[28] PAPOULIS, A.. Probability, Randon Variables, and Stochastic
Processes. McGrawHill, 1991.
[29] PHANI, A. S.. Damping identication in linear vibration. Uni-versity of Cambridge, 2000.
[30] QIU, Y.; RAO, S. S.. A fuzzy approach for the analysis of
unbalanced nonlinear rotor systems. Journal of Sounds andVibration, 284:299323, 1996.
[31] BARRIENTOS, G.. Dinâmica nãolinear de estruturas exíveis
em grandes deformações. Tese de Doutorado PUC-Rio, 1997.
[32] RITTO, T.; SAMPAIO, R.. The eects of unbalance and clearance
on the bearings of an overhung rotor. COBEM, 2005.
[33] ROSENBERG, R. M.. Analytical Dynamics of Discrete Systems.Plenum Press, New York, N.Y. USA, 1977.
[34] SAMPAIO, R.; LECKAR, H.. Aspectos matemáticos de vibrações
mecânicas. Notas correspondentes a um minicurso oferecido no Sem-inário Brasileiro de Análise, 1995. Disponível em http://www.mec.puc-rio.br/prof/rsampaio/rsampaio.html.
[35] SAMPAIO, R.; WOLTER, C.. Bases de Karhunen-Loève:
Aplicações à mecânica dos sólidos. APLICON 2001, SãoCarlos, S.P. Brasil, 2001. Disponível em http://www.mec.puc-rio.br/prof/rsampaio/rsampaio.html.
[36] SAMPAIO, R.; SOIZE, C.. On measures of nonlinearities for
uncertain dynamical systems. COBEM, 2005.
155
[37] RIQUELME, R.. Dinâmica de um manipulador robótico exível.Tese de Doutorado PUC-Rio, 1999.
[38] SHIGLEY; MISCHKE. Mechanical Engineering Design. McGraw-Hill, 1989.
[39] SINOU, J. J.; THOUVEREZ, F.. Non-linear dynamic of rotor-
stator system with non-linear bearing clearance. CR Mecanique,332:743750, 2004.
[40] SIROVICH, L.. Turbulence and the dynamics of coherent struc-
tures part I: coherent strucutures. Quarterly of Applied Mathemat-ics, 45, 1987.
[41] STRANG, G.; FIX, G.. An Analysis of The Finite Element
Method. PrenticeHall, Englewood Clis, N.J. USA, 1973.
[42] STRANG, G.. Linear Algebra and Its Applications. Harcourt, Inc.,Orlando, F.L. USA, 1988.
[43] TIWARI, M.; GUPTA, K. ; PRAKASH, O.. Dynamic response of an
unbalanced rotor supported on ball bearings. Journal of Soundsand Vibration, 238:757779, 2000.
[44] TRINDADE, M.; WOLTER, C. ; SAMPAIO, R..Karhunen-Loève Decomposition of coupled axial/bending
vibrations of beams subjected to impact. Journal of Sound andVibration, 279:10151036, 2005.
[45] VAKAKIS, A.; AZEEZ, M. A.. Numerical and experimental ana-
lysis of a continuos overhung rotor undergoing vibro-impacts.Non-Linear Mechanics, 34:415435, 1999.
[46] WOLTER, C.; TRINDADE, M. ; SAMPAIO, R.. Obtaining
mode shapes through the Karhunen-Loève expansion for
distributed-parameter linear system. Shock and Vibration, 9:177192, 2000.
[47] WOLTER, C.; SAMPAIO, R.. A comparison between the use of
Karhunen-Loève expansion and mode shapes in the model re-
duction of linear system. IX DINAME, Florianópolis, S.C. Brasil, 2001.Disponível em http://www.mec.puc-rio.br/prof/rsampaio/rsampaio.html.
156
[48] WOLTER, C.. Uma introdução à redução de modelos
através da expansão de Karhunen-Loève. Dissertação deMestrado PUC-Rio, 2001. Disponível em http://www.mec.puc-rio.br/prof/rsampaio/rsampaio.html.
[49] WOLTER, C.; TRINDADE, M. ; SAMPAIO, R.. Reduced-order model
for impacting beam using the Karhunen-Loève expansion.Tendências em Matemática Aplicada e Computacional, 3:217226, 2002.Disponível em http://www.mec.puc-rio.br/prof/rsampaio/rsampaio.html.
AAnálise de vibrações na internet
A pesquisa e estudos feitos através da rede de computadores, internet, éuma realidade que deve ser explorada. Destacamse as páginas dos professoresRubens Sampaio da PUC-Rio, Peter Avitabile da UMASS Lowell e Daniel A.Russell da Universidade de Kettering. Destacamse as páginas das empresasBently Nevada e Bruel and Kjaer. Destacamse os artigos contidos na páginahttp://www.vibetech.com/TechPapers.htm. Destacamse os pacotes emMATLAB: CALFEM e Structural Dynamics Toolbox.
A página do Professor Rubens Sampaio da PUC-Rio, http://www.mec.puc-rio.br/prof/rsampaio/rsampaio.html, tem diversos artigos sobre Vi-brações, Análise de Sinais, Decomposição de Karhunen-Loève, etc.
A página do professor Peter Avitabile da Universidade UMASS Lowell,http://macl.caeds.eng.uml.edu/umlspace/mspace.html, tem vários ar-tigos didáticos sobre Análise Modal tanto teórica quanto experimental.
A página do professor Daniel A. Russell da Universidade de Kettering,http://www.kettering.edu/~drussell, tem muitos exemplos ilustrativoscomo a análise acústica de um taco de baseball e de um violão.
Pela página da revista Orbit da empresa Benlty Nevada, http://www.bently.com/orbit.htm, podemse obter diversos artigos técnicos sobre ro-todinâmica.
Pela página da empresa Bruel and Kjaer, http://www.bksv.com, épossível obter muitos artigos técnicos sobre vibrações e acústica, além de poderse inscrever em seminários virtuais gratuitos.
Na página http://www.vibetech.com/TechPapers.htm são encontra-dos muitos artigos técnicos sobre vibrações e Análise Modal.
O pacote do MATLAB CALFEM, para o aprendizado em elementos nitos,é encontrado na página http://www.byggmek.lth.se/Calfem.
O pacote do MATLAB Structural Dynamics Toolbox, para Análise Modalexperimental é encontrado na página http://www.sdtools.com/sdt. Nestapágina tem uma seção com artigos sobre amortecimento e Análise Modal:
158
http://www.sdtools.com/Publications.html.
Além destas páginas, podese citar:
A página do professor S. Adhikari da Universidade de Bristol, http://www.aer.bris.ac.uk/contact/academic/adhikari/home.html, tem arti-gos sobre a modelagem da dissipação em sistemas dinâmicos.
A página http://www.ladenl.gov/projects/dss/Presentations.
htm tem algumas apresentações de professores como Peter Avitabile, NormHunter, entre outros.
A página do Departamento de Engenharia da Universidade de Perdue,http://widget.ecn.purdue.edu/~FEND2/teaching.html, mostra váriosexperimentos como a Análise Modal de uma raquete de tênis, um taco debaseball, um taco de golf, um helicóptero, etc.
A página http://www.dwsimpson.com/fourieranalysis.html apre-senta diversos links que abordam a Análise de Fourier.
A página do Departamento de Engenharia da Universidade de Cam-
bridge, http://www-control.eng.cam.ac.uk/extras/Virtual_Library/Control_VL.html apresenta links de grupos que trabalham na área de cont-role, além de links para sociedades e para periódicos.
A página http://www.myphysicslab.com tem algumas simulações desistemas simples de vibração.
A página da Escola e Engenharia de Jacobs,http://maeweb.ucsd.edu/~ugcl/index.html?experiments/experiments.
html, mostra diversos testes experimentais que ilustram princípios da dinâmicae controle.
A página do Instituto de Conabilidade de Equipamentos, http://www.equipment-reliability.com/articles/paper2.htm, apresenta vários ca-sos sobre vibração e isolamento de vibração.
A página da Universidade de Washington,http://octavia.ce.washington.edu/DrLayer/Exercises/L-BConcepts.
htm, apresenta diversos exercícios sobre vibrações e Análise Modal.
BMais Exemplos
Exemplo 1: Modos de vibração e freqüências naturais de uma viga:Considere a viga representada na Figura B.1:
Figura B.1: Viga engastada em uma extremidade
O programa modal_viga (ver Anexo E) calcula, através de expressõesanalíticas, os modos de vibração e as freqüências naturais de uma viga engastadaem uma extremidade e de uma viga biapoiada.
Este programa calcula os modos de vibração e as freqüências naturaisde uma viga engastada em uma extremidade e livre na outra pelas seguintesexpressões [18], [3]:
Modos de vibração:
φn(x) = cosh(βnx)− cos(βnx)− σn(senh(βnx)− sen(βnx)) (B-1)
Sendo:σn =
sinh(βnL)− sin(βnL)
cosh(βnL) + cos(βnL)(B-2)
E:cos(βL)cosh(βL) = 1 (B-3)
Logo:
β1 = 1, 86/L β2 = 4, 69/L β3 = 7, 85/L β4 = 11, 00/L β5 = 14, 14/L
(B-4)
160
e para n > 5
βn =(2n− 1)π
2L(B-5)
Freqüências naturais:
ωn = β2n
√EI
ρA(B-6)
A Figura B.2 mostra os cinco primeiros modos de vibrações de uma vigaengastada em uma extremidade:
Figura B.2: Cinco primeiros modos de uma viga engastada em uma extrem-
idade
As cinco primeiras freqüências naturais em Hz são:4, 5299
28, 3884
79, 4882
155, 7652
257, 4910
(B-7)
Agora considere a viga biapoiada representada na Figura B.3:
161
Figura B.3: Viga biapoiada
Os modos de vibrações e as freqüências naturais são calculadas pelasseguintes expressões [18], [3]:
Modos de vibração:
φn(x) = sen(λnx) (B-8)
Sendo:λn =
nπ
L(B-9)
Freqüências naturais:
ωn = λ2n
√EI
ρA(B-10)
A Figura B.4 mostra os cinco primeiros modos de vibrações de uma vigabiapoiada:
Figura B.4: Cinco primeiros modos de uma viga biapoiada
162
As cinco primeiras freqüências naturais em Hz são:12, 7156
50, 8624
114, 4404
203, 4496
317, 8900
(B-11)
Exemplo 2: Exemplo de uma estrutura discretizada pelo MEF:
Nesta seção será feito um exemplo de uma estrutura em 'L' compostapor materiais e dimensões diferentes. Este exemplo objetiva explicar o funciona-mento das entradas de dados do programa.
Imagine a estrutura representada na Figura (B.5).
Figura B.5: Estrutura
A Figura (B.6) mostra a estrutura representada na Figura (B.5) dis-cretizada pelo MEF. Um total de quatro elementos e cinco nós são usados.
163
Figura B.6: Estrutura discretizada
Esta conguração está no programa modal_l2_ef (ver Anexo E) . Esteprograma calcula as freqüências naturais e os modos de vibração da estruturadiscretizada, Figura (B.6).
As entradas do programa serão vistas passo a passo:
1) "Dados do Material": nesta parte do programa são denidas as pro-priedades dos materiais e dimensões das hastes, ver Figura (B.5). As entradassão: módulo de elasticidade do material, densidade do material, área da seçãotransversal e momento de inércia das hastes. A haste horizontal é de alumínioe a haste vertical é de aço (ver Figura B.5).
2)"Topologia": nesta parte do programa a numeração dos elementos egraus de liberdades é feita:
Edof =
1 1 2 3 4 5 6
2 4 5 6 7 8 9
3 7 8 9 10 11 12
4 10 11 12 13 14 15
(B-12)
A primeira coluna da matriz Edof representa a numeração dos elementos.Cada nó tem 3 graus de liberdade: deslocamento nas direções x e y e deslo-camento angular em torno do eixos z. Os graus de liberdade do elemento 1,primeira linha da matriz (B-12), são: 1, 2, 3, 4, 5 e 6; os graus de liberdade doelemento 2, segunda linha da matriz (B-12), são: 4, 5, 6, 7, 8 e 9. Os elementos
164
têm nós em comum.
Coord =
0 0
0 1, 5
0 3
1 3
2 3
(B-13)
A matriz Coord, é composta pelas coordenadas de cada nó. Na primeiracoluna encontramse as coordenadas x e na segunda as coordenadas y. O nó3 tem coordenada x = 0 e y = 3, por exemplo.
Dof =
1 2 3
4 5 6
7 8 9
10 11 12
13 14 15
(B-14)
A matriz Dof é compostas pelos os números de graus de liberdadeassociados a cada nó; cada linha representa um nó. O nó 2 tem três grausde liberdade 4, 5 e 6, por exemplo.
O comando coordxtr monta duas matrizes: Ex e Ey, que são compostospelas coordenadas x e y de cada elemento.
Ex =
0 0
0 0
0 1
1 2
(B-15)
Cada linha representa um elemento. o elemento 3, por exemplo, temx = 0 para o nó um e x = 1 para o nó dois.
[Ex,Ey] = coordxtr(Edof, Coord,Dof, 2), onde a entrada "2"signicaque são dois nós por elemento.
3) "Cálculo das Matrizes dos elementos": o comando beam2d calcula asmatrizes dos elementos de viga (modelo EulerBernoulli).
[k,m, c] = beam2d(Ex(1, :), Ey(1, :), ep).
Sendo ep = [E A I rho ∗ A α β].
165
Neste exemplo os elementos 1 e 2 têm as mesmas propriedades mecânicas(contidas no vetor epv) e os elementos 3 e 4 têm outras propriedades mecânicas(vetor eph). A matriz de amortecimento do elemento é calculada da seguinteforma: αm+ βk.
3) "Montagem das Matrizes dos elementos": Estas matrizes são mon-tadas pelo comando assem formando as matrizes globais M, K e C.
4) "Cálculo das Freqüências Naturais de Modos de Vibração": O comandoeigen: [La,Egv] = eigen(K,M, b) calcula as freqüências naturais (La) e modosde vibrações (Egv).
O vetor b contém os graus de liberdade que estão restritos. Neste exemploo primeiro nó tem os 3 graus de liberdade restritios, ou seja, b = [1 2 3]T .O que acontece é que as linhas 1, 2 e 3 e as colunas 1, 2 e 3 são excluídas docálculo, já que não há deslocamento referente a estes graus de liberdade.
Resultado:
A primeira freqüência natural é 0, 7863 Hz, correspondente ao primeiromodo de vibração representado na Figura B.7.
Figura B.7: Primeiro modo de vibração
Exemplo 3:
O programa modal_viga2_ef (ver Anexo E) calcula, discretizando peloMEF, os modos e freqüências naturais de uma viga engastadalivre, para uma
166
determinada precisão. Considere a viga representada abaixo:
Erro admitido de 0, 002% para o quinto modo de vibração. O número deelementos inicial é 20 e o passo entre iterações é de 20 elementos. O númerode elementos necessários para a precisão requerida é de 60 elementos. As cincoprimeiras freqüências naturais (Hz) são:
4, 53
28, 39
79, 49
155, 77
257, 49
(B-16)
Os cinco primeiros modos de vibração estão representados na Figura B.8:
Figura B.8: Cinco primeiros modos de uma Viga Engastadalivre
O erro do quinto modo de vibração para cada passo da simulação é dadopor:
167
Figura B.9: Convergência da aproximação
A curva vermelha mostra o erro da quinta freqüência natural (MEF xexpressão analítica). A curva verde mostra o erro da quinta freqüência natural(MEF). A curva azul mostra o erro do quinto modo de vibração (MEF). Masqual é a inuência do erro arbitrado na quantidade de elementos necessáriospara a análise? E qual a inuência da quantidade de modos que se quer obterdeterminada precisão, no número de elementos necessários?
Para responder estas perguntas rodouse o programa arbitrandose umerro igual a metade do arbitrado anteriormente (e = 0, 002/2). Neste casosão necessários 80 elementos nitos para satisfazer a nova precisão, 20 el-ementos a mais. Já quando o número de modos que se quer determinadaprecisão dobra (N = 2.5 = 10), são necessários 100 elementos nitos para aanálise, 40 elemento a mais. Ou seja, para o modelo analisado é muito maisdispendioso dobrar o número de modos calculados para determinada precisãodo que dividir pela metade o erro arbitrado para determinado número de modos.
Exemplo 4:
O programa modal_viga3_ef (ver Anexo E) calcula, discretizando peloMEF, os modos e freqüências naturais de uma viga biapoiada, para umadeterminada precisão. Considere a viga representada abaixo:
168
Erro admitido de 0, 005% para o quinto modo de vibração. O número deelementos inicial é 20 e o passo entre iterações é de 20 elementos. O númerode elementos necessários para a precisão requerida é de 60 elementos. As cincoprimeiras freqüências (Hz) naturais são:
12, 71
50, 86
114, 44
203, 45
317, 89
(B-17)
Os cinco primeiros modos de vibração estão representados na Figura B.10:
Figura B.10: Cinco primeiros modos de uma viga Biapoiada
O erro do quinto modo de vibração para cada passo da simulação é dadopor:
169
Figura B.11: Convergência da aproximação
A curva vermelha mostra o erro da quinta freqüência natural (MEF xexpressão analítica). A curva verde mostra o erro da quinta freqüência natural(MEF). A curva azul mostra o erro do quinto modo de vibração (MEF).
Exemplo 5:
O programa modal_l_ef (ver Anexo E) calcula, discretizando pelo MEF,os modos e freqüências naturais de uma viga em conguração L engastada,para uma determinada precisão. Observe gura abaixo:
Figura B.12: Viga conguração L
Erro admitido de 0,005% para o quinto modo de vibração. O número deelementos inicial é 20 e o passo da interação é de 20 elementos. O número deelementos necessários para a precisão requerida é de 60 elementos. As cinco
170
primeiras freqüências (Hz) naturais são:1.51
4.11
20.29
29.73
64.86
(B-18)
Os cinco primeiros modos de vibração estão representados na Figura B.13:
Figura B.13: Cinco primeiros modos da conguração L
O erro do quinto modo de vibração para cada passo da simulação é dadopor:
171
Figura B.14: Convergência da aproximação
Exemplo 6:
O programa modal_triang_ef (ver Anexo E) calcula, discretizando peloMEF, os modos e freqüências naturais de uma viga em conguração Triângulo,para uma determinada precisão. Observe a gura abaixo:
Figura B.15: Viga conguração triângulo
Erro admitido de 0, 01% para o quinto modo de vibração. O número deelementos inicial é 15 e o passo da interação é de 9 elementos. O número deelementos necessários para a precisão requerida é de 132 elementos. As cinco
172
primeiras freqüências (Hz) naturais são:0.00
3.54
7.31
20.41
27.80
(B-19)
A primeira freqüência natural com valor zero representa movimento decorpo rígido. Os cinco primeiros modos de vibração estão representados naFigura B.16:
Figura B.16: Cinco primeiros modos da conguração triângulo
O erro do quinto modo de vibração para cada passo da simulação é dadopor:
CLista dos Programas Desenvolvidos
A seguir uma lista com os programas principais desenvolvidos. Alémdesses programas foram desenvolvidos programas auxiliares. A relação entreos programas está no anexo D. A Figura C.1 mostra um diagrama de uxo deum programa geral.
Dados Gerais - geometria; - propriedades do material;
Principal
Precisão - erro estipulado; - elementos / modos iniciais - passo
Integração - método (ode45) - tempo de simulação - passo de integração;
Análise - Tipo de análise:
- Análise Modal - Resposta Dinâmica - FRF - FFT - Decomposição de KL
Erro < precisão?
Visualização dos
resultados
SIM
NÃO
Escolha das funções teste - MEF - Modos Supostos
Figura C.1: Diagrama de uxo
175
Nome do programa Propósitomodal_barra Calcular os modos de vibração e freqüências
naturais de uma barra, a partir de expressõesanalíticas. Duas condições de contorno podemser determinadas: xalivre ou xaxa
modal_barra_ef Calcular os modos de vibração e freqüênciasnaturais de uma barra, discretizada pelo MEF,dada uma certa precisão. Três condições decontorno podem ser determinadas: xa-livre,xa-mola ou xa-massa.
modal_viga Calcular os modos de vibração e freqüênciasnaturais de uma viga, a partir de expressõesanalíticas. Duas condições de contorno podemser determinadas: engastadalivre ou biapoiada.
modal_viga2_ef Calcular os modos de vibração e freqüênciasnaturais de uma viga, discretizada pelo MEF,dada uma certa precisão. Condições decontorno: engastada-livre.
modal_viga3_ef Calcular os modos de vibração e freqüênciasnaturais de uma viga, discretizada pelo MEF,dada uma certa precisão. Condições decontorno: bi-apoiada.
modal_l_ef Calcular os modos de vibração e freqüênciasnaturais de uma estrutura em L, discretizadapelo MEF, dada uma certa precisão.
modal_l2_ef Exemplicar passo a passo as entradas dosprogramas que utilizam o MEF. Este programacalcula os modos de vibração e freqüências naturaisde uma estrutura em L discretizado pelo MEF.
modal_triang_ef Calcular os modos de vibração e freqüênciasnaturais de uma estrutura triângulo, discretizadapelo MEF, dada uma certa precisão.
176
modal Calcular os modos de vibração, freqüências naturaise coecientes de amortecimento de um sistema dinâmicolinear com amortecimento proporcional.
est_viga_ms Calcular o deslocamento da extremidade de uma viga,discretizada pelo Método dos Modos Supostos.Condição de contorno engastada-livre
est_viga_ef Calcular o deslocamento da extremidade de uma viga,discretizada pelo MEF. Condição de contornoengastada-livre.
din_viga_ms Calcular a dinâmica de uma viga, discretizadapelo Método dos Modos Supostos, dada uma certaprecisão. Condição de contorno engastada-livre.
din_viga_ef Calcular a dinâmica de uma viga, discretizadapelo MEF, dada uma certa precisão. Condição decontorno engastada-livre.
frf Calcular a FRF e os parâmetros modais de umsistema linear com amortecimento proporcional.
frf2 Calcular a FRF e os parâmetros modais de umsistema de uma dimensão (inclusive pólo e resíduo).
frf_viga_ef Calcular a FRF para excitação e respostana extremidades livre de uma viga engastadalivre.
t_ex Calcular a FFT de um sinal no domínio do tempomassamola_kl Calcular a dinâmica de um sistema com amortecimento
nãoproporcional a partir das bases de KL.disco_kl Calcular a dinâmica de um disco, modelado
como uma massa girante, a partir das bases de KL.barra_choque Calcular a resposta dinâmica do sistema e analisar
os impactos e dinâmica do ponto da barra juntoao anteparo. Os modos de vibração saocalculados por expressões analíticas. Umaaproximação que utiliza o método dos ModosSupostos é feita para gerar as dinâmicas.A Decomposição de KL eh usada parareduzir o sistema, que é não-linear.
177
rotor_choque Calcular a resposta dinâmica do sistema e analisaros impactos e dinâmica dos pontos do eixojunto aos mancais. Comparar com resultadosobservados em equipamentos. O modelo de EulerBernoulli e usado. O MEF eh usado para calcularos modos de vibração de uma viga com umamassa na extremidade. Uma aproximaçãoque utiliza o método dos Modos Supostos éfeita para gerar as dinâmicas.A Decomposição de KL eh usada parareduzir o sistema, que eh não-linear.
Outros trabalhos desenvolvidos na PUCRio seguiram esta linha de de-senvolvimento de programas em MATLAB. Riquelme, [37], analisou a dinâmicade um manipulador robótico exível e Barrientos, [31], analisou a dinâmicanãolinear de estruturas exíveis em grandes deformações.
DRelação entre os Programas
Neste anexo é mostrada a relação entre os programas desenvolvidos.Os programas são divididos em três grandes grupos. O primeiro grupo
de programas serve para calcular as características modais de vigas, barras eestruturas formadas por vigas, Figura D.1. O segundo grupo de programas servepara analisar um problema no domínio do freqüência, Figura D.3. O terceirogrupo de programas resolve problemas usando a Decomposição de KarhunenLoève, Figura D.4.
Por m é mostrada a interação entre os programas usados para simularum rotor, Figura D.5.
179
Parâmetros Modais
1) modal_barra
2) modal_barra_ef
3) modal_viga
4 a 8) modal_*_ef
calfem_coordxtr calfem_bar2d calfem_assem calfem_eigen
calfem_coordxtr calfem_beam2d calfem_assem calfem_eigen
calfem_eldraw2 calfem_extract calfem_eldisp2
Para visualisação dos resultados
Cria matriz com coordenadas
Cria matrizes dos elementos
Monta matriz global
Cria matriz com coordenadas
Cria matrizes dos elementos
Monta matriz global
9) modal
Figura D.1: Primeiro grupo de programas
180
11) est_viga_ef
calfem_coordxtr calfem_beam2d calfem_beam2e calfem_assem
Cria matriz com coordenadas
Cria matrizes dos elementos
Cria vetor de força do elemento
10) est_viga_ms
13) din_viga_ef
calfem_coordxtr calfem_beam2d calfem_assem calfem_eigen
Cria matriz com coordenadas
Cria matrizes dos elementos
Monta matriz global Calcula autovalores e autovetores
12) din_viga_ms
edo sub_erro
Integra o sistema de equações
Calcula o erro da aproximação
sub_monta_km edo2 sub_erro
Integra o sistema de equações
Calcula o erro da aproximação
Programas para o Cálculo da Resposta dinâmica ou estática
Figura D.2: Segundo grupo de programas
181
14) frf
15) frf2
16) frf_viga_ef
calfem_coordxtr calfem_bar2d calfem_assem calfem_eigen
Criação da matriz com as coordenadas
Criação das matrizes dos elementos
Montagem matriz global
Calcula autovalores e autovetores
17) fft_ex
sub_fft
Calcula a FFT
Domínio da Freqüência
Figura D.3: Terceiro grupo de programas
182
18) massamola_kl
edo_kl
Integração do sistema de equações
19) disco_kl
disco_edo
Integração do sistema de equações
20) barra_choque
sub_modos
sub_edo
sub_erro
Integra o sistema de equações
Calcula o erro da aproximação
sub_dinamica
Calcula os modos e freqüências naturais
sub_kl sub_rec
Calcula a base de KL e os valores
próprios
Reconstrói a dinâmica
Decomposição de Karhunen-Loève
Figura D.4: Quarto grupo de programas
183
21) rotor_choque
sub_modosnormais
sub_edo
sub_erro
Integra o sistema de equações
Calcula o erro da aproximação
sub_dinamica sub_kl sub_rec
Calcula a base de KL e os valores
próprios
Reconstrói a dinâmica
calfem_coordxtr calfem_beam2d calfem_assem calfem_eigen
Cria matriz com coordenadas
Cria matrizes dos elementos
Monta matriz global Calcula autovalores e autovetores
Modelo Rotor-Mancal
Figura D.5: Simulação de um rotor
EManual dos Programas Desenvolvidos
Diversos programas foram desenvolvidos no transcorrer deste trabalho,todos através da plataforma MATLAB.
A nomenclatura das variáveis seguiu o seguinte padrão:
MATRIZ com entradas reais −→ RMat_MATRIZ.MATRIZ com entradas imaginárias −→ IMat_MATRIZ.MATRIZ com entradas complexas −→ CMat_MATRIZ.V ETOR com entradas reais −→ R_V ETOR.V ETOR com entradas imaginárias −→ I_V ETOR.V ETOR com entradas complexas −→ C_V ETOR.ESCALAR −→ ESCALAR.
Esta nomenclatura facilitará a compreensão dos programas. Além disso,todos os programas têm um cabeçalho onde é descrito o propósito, as variáveisde entrada e as variáveis de saídas. As equações diferenciais foram resolvidaspela subrotina ode45 do MATLAB, programas ode e ode2. A seguir os pro-grama desenvolvidos:
1) modal_barra-Analise modal de uma barra-
PROPÓSITO:Calcular os modos de vibração e freqüências naturais de uma barra, a partirde expressões analíticas. Duas condições de contorno podem ser determinadas:xalivre ou xaxa:
185
ENTRADAS:L : comprimento da barra (m)E : modulo de elasticidade (Pa)rho : densidade (Kg/m3)N : Número de modos a serem calculados
SAÍDAS:R_omega_n : freqüências naturais (rd/s)RMat_modo : modos de vibração-
2) modal_barra_ef-Analise modal de uma barra, via MEF-PROPÓSITO: Calcular os modos de vibração e freqüências naturais de umabarra, discretizada pelo MEF, dada uma certa precisão. Três condições decontorno podem ser determinadas: xa-livre, xa-mola ou xa-massa.
ENTRADAS:L : comprimento da barra (m)h : altura da barra (m)b : espessura da barra (m)E : modulo de elasticidade (Pa)rho : densidade (Kg/m3)kk : constante da mola da extremidade (N/m)mm : massa na extremidade (Ns/m)e : precisão especicada (%)N : Número de modos que se quer obter com a precisão especicada
186
NE_ini : número de elementos iniciais (número PAR!!)passo : passo do número de elementos (número PAR!!)
SAÍDAS:R_Freq / R_omega_n2 / R_omega_n3 : freqüências naturais (Hz)RMat_Egv / RMat_modo2 / RMat_modo3 : modos de vibraçãoNE_prec : numero de elementos necessários para a precisão requerida-
3) modal_viga-Analise modal de uma viga-
PROPÓSITO:Calcular os modos de vibração e freqüências naturais de uma viga, a partirde expressões analíticas. Duas condições de contorno podem ser determinadas:engastadalivre ou biapoiada:
ENTRADAS:L : comprimento da barra (m)h : altura da barra (m)b : espessura da barra (m)E : modulo de elasticidade (Pa)rho : densidade (Kg/m3)N : Número de modos a serem calculados
SAÍDAS:R_omega_n : freqüências naturais (Hz)RMat_modo : modos de vibração-
4) modal_viga2_ef-
187
Analise modal de uma viga Engastadalivre, discretizada pelo MEF-
PROPÓSITO:Calcular os modos de vibração e freqüências naturais de uma viga, discretizadapelo MEF, dada uma certa precisão. Condições de contorno: engastada-livre:
ENTRADAS:L : comprimento da barra (m)h : altura da barra (m)b : espessura da barra (m)E : modulo de elasticidade (Pa)rho : densidade (Kg/m3)alfa : coecientes de amortecimento C=aM+bKbeta :e : precisão especicada (%)N : Número de modos que se quer obter com a precisão especicadaNE_ini : número de elementos iniciais (número PAR!!)passo : passo do número de elementos (número PAR!!)
SAÍDAS:R_Freq : freqüências naturais (Hz)RMat_Egv : modos de vibraçãoNE_prec : numero de elementos necessários para a precisão requerida-
5) modal_viga3_ef-Analise modal de uma viga biapoiada, discretizada pelo MEF-
PROPÓSITO:Calcular os modos de vibração e freqüências naturais de uma viga, discretizadapelo MEF, dada uma certa precisão. Condições de contorno: bi-apoiada:
188
ENTRADAS:L : comprimento da barra (m)h : altura da barra (m)b : espessura da barra (m)E : modulo de elasticidade (Pa)rho : densidade (Kg/m3)alfa : coecientes de amortecimento C=aM+bKbeta :e : precisão especicada (%)N : Número de modos que se quer obter com a precisão especicadaNE_ini : número de elementos iniciais (número PAR!!)passo : passo do número de elementos (número PAR!!)
SAÍDAS:R_Freq : freqüências naturais (Hz)RMat_Egv : modos de vibraçaoNE_prec : numero de elementos necessários para a precisão requerida-
6) modal_l_ef-Analise modal de uma estrutura em L, discretizada pelo MEF-
PROPÓSITO:Calcular os modos de vibração e freqüências naturais de uma estrutura em L,discretizada pelo MEF, dada uma certa precisão.
189
L
L
ENTRADAS:L : comprimento da barra (m)h : altura da barra (m)b : espessura da barra (m)E : modulo de elasticidade (Pa)rho : densidade (Kg/m3)alfa : coecientes de amortecimento C=aM+bKbeta :e : precisão especicada (%)N : Número de modos que se quer obter com a precisão especicadaNE_ini : número de elementos iniciais (número PAR!!)passo : passo do número de elementos (número PAR!!)
SAÍDAS:R_Freq : freqüências naturais (Hz)RMat_Egv : modos de vibraçaoNE_prec : numero de elementos necessários para a precisão requerida-
7) modal_l2_ef-Analise modal de uma estrutura em L, discretizada pelo MEF-
PROPÓSITO:Exemplicar passo a passo as entradas dos programas que utilizam o MEF. Esteprograma calcula os modos de vibração e freqüências naturais de uma estruturaem L discretizado pelo MEF:
190
L1
x
y
L2
ENTRADAS:R_epv=[Ev Av Iv rhov*Av] : Dados do material e dimensão das vigas verticaisR_eph=[Eh Ah Ih rhoh*Ah] : Dados do material e dimensão das vigas hori-zontaisEv / Eh : modulo de elasticidade do material das vigas (h)horizontais e(v)verticais - (Pa)Av / Ah : area da seção transversal das vigas (h) e (v) - (m2)Iv / Ih : momento de inércia das vigas (h) e (v) - (m4)rhov / rhoh : densidade do material das vigas (h) e (v) - (Kg/m3)RMat_Edof : cada linha eh composta pelo n de elemento e numeração dosgraus de liberdade (6 por elemento)RMat_Coord : cada linha eh composta pelas coordenadas: x1 y1; x2 y2;...RMat_Dof : cada linha eh composta pela numeração dos graus de liberdadeassociados a cada nóR_b : graus de liberdade restritos
SAÍDAS:R_Freq : freqüências naturais (Hz)RMat_Egv : modos de vibração-
8) modal_triang_ef-Analise modal de uma estrutura triângulo, discretizada pelo MEF-
PROPÓSITO:Calcular os modos de vibração e freqüências naturais de uma estrutura triângulo,discretizada pelo MEF, dada uma certa precisão:
191
L
dist=0,7L
45º
ENTRADAS:L : comprimento da barra (m)h : altura da barra (m)b : espessura da barra (m)E : modulo de elasticidade (Pa)rho : densidade (Kg/m3)alfa : coecientes de amortecimento C=aM+bKbeta :e : precisão especicada (%)N : Número de modos que se quer obter com a precisão especicadaNE_ini : número de elementos iniciais (múltiplo de TRÊS!!)passo : passo do número de elementos (múltiplo de TRÊS!!)
SAÍDAS:R_Freq : freqüências naturais (Hz)RMat_Egv : modos de vibraçaoNE_prec : numero de elementos necessários para a precisão requerida-
9) modal-Exemplo simples de Análise Modal-
PROPÓSITO:Calcular os modos de vibração, freqüências naturais e coecientes de amortec-imento de um sistema dinâmico linear com amortecimento proporcional.
ENTRADAS:RMat_M : matriz de massaRMat_K : matriz de rigidez
192
alfa : coecientes Cp = (alfa.M + beta.K)beta :
SAÍDAS:RMat_w_n : freqüências naturaisRMat_modo : modos de vibraçãoRMat_w_d : freqüências de oscilação para o sistema amortecidoRMat_c_amort : coecientes de amortecimento modais-
10) est_viga_ms-Analise estática de uma viga engastada-livre discretizada pelo Método dosModos Supostos-
PROPÓSITO:Calcular o deslocamento da extremidade de uma viga, discretizada pelo Métododos Modos Supostos. Condição de contorno engastada-livre:
ENTRADAS:L : comprimento da vigah : altura da vigab : espessura da vigaE : modulo de elasticidaderho : densidadef_f : amplitude do forçamento carregamento concentradoq : força por unidade de comprimento carregamento distribuídoN : Número de modos a serem usados na simulação
SAÍDAS:desl : deslocamento calculado pela expressão FL3/(3EI) carregamento pon-tual
193
des : deslocamento calculado pela expressão qL4/(8EI) carregamento dis-tribuídow3 : deslocamento calculado pela aproximação para o carregamento pontualu3 : deslocamento calculado pela aproximação para o carregamento distribuído-
11) est_viga_ef-Analise estática de uma viga engastada-livre, discretizada pelo MEF-
PROPÓSITO:Calcular o deslocamento da extremidade de uma viga, discretizada pelo MEF.Condição de contorno engastada-livre:
ENTRADAS:L : comprimento da vigah : altura da vigab : espessura da vigaE : modulo de elasticidaderho : densidadef_f : amplitude do forçamento carregamento concentradoq : força por unidade de comprimento carregamento distribuídoNE : Número de elementos a serem usados na simulação
SAÍDAS:desl : deslocamento calculado pela expressão FL3/(3EI) carregamento pon-tualdes : deslocamento calculado pela expressão qL4/(8EI) carregamento dis-tribuídodesl1 : deslocamento calculado pela aproximação para o carregamento pontualdesl2 : deslocamento calculado pela aproximação para o carregamento dis-tribuído
194
-
12) din_viga_ms-Dinâmica da viga engastada-livre aproximada pelo Método dos Modos Supostos-
PROPÓSITO:Calcular a dinâmica de uma viga, discretizado pelo Método dos Modos Supostos,dada uma certa precisão. Condição de contorno engastada-livre:
ENTRADAS:L : comprimento da barra (m)h : altura da barra (m)b : espessura da barra (m)E : modulo de elasticidade (Pa)rho : densidade (Kg/m3)c : constante de amortecimento (Ns/m)N: número de modos iniciais usados na simulaçãox_f: local do forçamento (m)f_f: amplitude do forçamento (m)omega_f: freqüência de forçamento (rd/s)resp: local onde o deslocamento será observadott: tempo total de simulação (s)e: erro admitido para a resposta dinâmica em L (%)OBS. SÓ FOI IMPLEMENTADO O CALCULO DO ERRO PARA A RES-POSTA EM resp=L!!
SAÍDAS:R_w_n : freqüências naturais (Hz)R_w_n : freqüências de oscilação do sistema amortecido (Hz)RMat_modo3 : modos de vibraçãoN_prec : numero de modos necessários para a precisão requerida
195
-
13) din_viga_ef-Dinâmica da viga engastada-livre aproximada pelos MEF-
PROPÓSITO:Calcular a dinâmica de uma viga, discretizada pelo MEF, dada uma certaprecisão. Condição de contorno engastada-livre:
ENTRADAS:L : comprimento da barra (m)h : altura da barra (m)b : espessura da barra (m)E : modulo de elasticidade (Pa)rho : densidade (Kg/m3)NE: número de modos iniciais usados na simulaçãopasso: passo do número de elementosNO_f: nó em que será aplicado o forçamento (m)Amp_f: amplitude do forçamento (m)omega_f: freqüência de forçamento (rd/s)NO_r: nó em que o deslocamento será observadott: tempo total de simulação (s)e: erro admitido para a resposta dinâmica em L (%)
SAÍDAS:R_Freq : freqüências naturais (Hz)RMat_Egv : modos de vibraçaoNE_prec : numero de modos necessários para a precisão requerida-
196
14) frf-Exemplo simples de FRF-
PROPÓSITO:Calcular a FRF e os parâmetros modais de um sistema linear com amorteci-mento proporcional.
ENTRADAS:RMat_M : matriz de massaRMat_K : matriz de rigidezalfa : Cp = alfa.M + beta Kbeta :w_ini : freq inicial de interaçãow_n : freq nal de interaçãow_fac : passo de freq
SAÍDAS:RMat_w_n : freq naturaisRMat_w_d : freq de oscilação do sistema amortecidoRMat_modo : modos de vibraçõesRMat_c_amort : coecientes de amortecimentoR_FRF_1 : FRF do sistema conservativoR_FRF_2 : FRF do sistema amortecido-
15) frf2-Exemplo simples de FRF - pólo e residuo-
PROPÓSITO:Calcular a FRF e os parâmetros modais de um sistema de uma dimensão(inclusive pólo e resíduo):
197
ENTRADAS:m : massak : molac : amortecimento
SAÍDAS:wn : freqüência natural (rd/s)wd : freqüência de oscilação do sistema com amortecimento (rd/s)xi : coeciente de amortecimentop : póloR : resíduoH : função de transferência-
16) frf_viga_ef-FRF de uma viga engastadalivre, discretizada por elementos nitos-
PROPÓSITO:Calcular a FRF para excitação e resposta na extremidades livre de uma vigaengastadalivre:
ENTRADAS:L : comprimento da vigah : altura da vigab : espessura da viga
198
E : modulo de elasticidaderho : densidadealfa : coecientes : C=aM+bKbeta :w_n : faixa de freq de 0 a omega_n Hzw_fac : passo da freqe : precisão especicadaNE_ini : número de elementos iniciaispasso : passo do número de elementos
SAÍDAS:R_Freq : freqüências naturaisRMat_Egv : modos de vibraçãoR_FRF : Função Resposta em Freqüência-
17) t_ex-Exemplo FFT-
PROPÓSITO:Calcular a FFT de um sinal no domínio do tempo.
ENTRADAS:CPM : freqüência em CPMA : amplitude do sinalt1 : tempo inicialt2 : tempo naldt : discretização
SAÍDAS:R_FFT : FFT do sinal senoidal-
18) massamola_kl-Decomposição de KL em um sistema simples-
199
PROPÓSITO:Calcular a dinâmica de um sistema com amortecimento nãoproporcional apartir das bases de KL.
ENTRADAS:RMat_M : Matriz de massaRMat_C : Matriz de amortecimentoRMat_K : Matriz de rigidez
SAÍDAS:R_u_r1 : resposta dinâmica u_r1R_u_r2 : resposta dinâmica u_r2R_Lambda : POVs valores ortogonais própriosRMat_Q_1 : POMs modos ortogonais próprios-
19) disco_kl-Decomposição de KL de um disco-
PROPÓSITOCalcular a dinâmica de uma disco, modelado como uma massa girante, a partirdas bases de KL.
200
X1X2
ENTRADAS:d : diâmetroI : momento de inérciaIp : momento de inércia polarA : seção transversalE : modulo de elasticidadev : coeciente de poisonrho : densidadeG : modulo de cisalhamentoesp : espessurakt : restrição de rotaçãoomega : rotação do eixo
SAÍDAS:R_u_r1 : resposta dinâmica u_r1R_u_r2 : resposta dinâmica u_r2R_Lambda : POVs valores ortogonais própriosRMat_Q_1 : POMs modos ortogonais próprios-
-
20) barra_choque- PROGRAMA PRINCIPAL SISTEMA DE VIBRO-IMPACTO DE UMA BARRA -
PROPÓSITO:Calcular a resposta dinâmica do sistema e analisar os impactos e dinâmicado ponto da barra junto ao anteparo. Os modos de vibração sao calculadospor expressões analíticas. Uma aproximação que utiliza o método dos Modos
201
Supostos é feita para gerar as dinâmicas. A Decomposição de KL eh usada parareduzir o sistema, que é não-linear:
ENTRADAS:Ki : coeciente de rigidez do mancalE : modulo de elasticidaderho : densidade (Kg/m3)c : coeciente de amortecimentoL : comprimento do eixod : diâmetro do eixofolga : folga (clearance)Hz : freqüência de excitacaoF_f : amplitude de forcamentoN : Numero de modos usados na aproximação (modos normais)N_kl : Numero de bases de KL usadas na aproximação (modos supostos)dt : delta tR_tspan : tempo de simulacao
SAÍDAS:programa "sub_modos"R_omega_n : freqüências naturais (Hz)RMat_modo : modos de vibraçao
programa "sub_kl"RMat_kl : matrizes com os modos empiricosR_Lambda : vetores com os valores ortogonais proprios
programa "sub_dinamica"RMat_u : resposta dinamicaR_F : forca media no anteparo-
202
21) rotor_choque- PROGRAMA PRINCIPAL SISTEMA DE VIBROIMPACTO DE UM ROTOR EM BALANCO COMDOIS ANTEPAROS (MANCAIS)-
PROPÓSITO:Calcular a resposta dinâmica do sistema e analisar os impactos e dinâmica dospontos do eixo junto aos mancais. Comparar com resultados observados emequipamentos. O modelo de Euler Bernoulli e usado. O MEF eh usado paracalcular os modos de vibração de uma viga com uma massa na extremidade.Uma aproximação que utiliza o método dos Modos Supostos eh feita para geraras dinâmicas. A Decomposição de KL eh usada para reduzir o sistema, que ehnão-linear:
ENTRADAS:Mr : massa do rotord1 : coeciente de amortecimento da viga na direção 1d2 : coeciente de amortecimento da viga na direção 2Ki : coeciente de rigidez do mancalE : modulo de elasticidaderho : densidade (Kg/m3)L : comprimento do eixoDe : diâmetro do eixoA : area da seção transversalIe : momento de inércia do eixofolga : folga (clearance)x_c1 : posicionamento do mancal 1
203
x_c2 : posicionamento do mancal 2Dr : diâmetro do rotorIr : momento de inércia do discoRPM : rotação do eixo em RPMG : grau de qualidade do balanceamentoerr : precisão exigida para o calculo dos modos e freq nat por EFN : Numero de modos usados na aproximação (modos normais)N_kl1 : Numero de bases de KL usadas na aproximação de u1 (modos supos-tos)N_kl2 : Numero de bases de KL usadas na aproximação de u2 (modos supos-tos)dt : delta tR_tspan : tempo de simulação
SAÍDAS:programa "sub_modosnormais"R_Freq : freqüências naturais (Hz)RMat_Egv : modos de vibraçãoNE_prec : numero de elementos necessários para a precisão requerida
programa "sub_kl"RMat_kl1 / RMat_kl2 : matrizes com os modos empíricos para u1 e u2R_Lambda2 / R_Lambda2 : vetores com os valores ortogonais próprios parau1 e u2
programa "sub_dinamica"RMat_u1 RMat_u2 : resposta dinâmica para u1 e u2.
programa "sub_t"R_t : FFT de uma resposta dinâmica especicada como entrada para oprograma-
FArquivos do CALFEM
Para o cálculo das matrizes M , C e K obtidas através da dis-cretização pelo MEF, arquivos do programa CALFEM foram utiliza-dos. O manual do programa CALFEM [5] pode ser obtido no sitehttp://www.byggmek.lth.se/Calfem. Os arquivos do CALFEM utiliza-dos foram os seguintes:
coordxtr −→ extrai das coordenadas globais dados para as coordenadaslocais.
207
beam2e −→ calcula o vetor de carregamento de um elemento que comcarregamento distribuído uniforme.
211
eigen −→ resolve o problema de autovalor, dado os nós restritos.
Livros Grátis( http://www.livrosgratis.com.br )
Milhares de Livros para Download: Baixar livros de AdministraçãoBaixar livros de AgronomiaBaixar livros de ArquiteturaBaixar livros de ArtesBaixar livros de AstronomiaBaixar livros de Biologia GeralBaixar livros de Ciência da ComputaçãoBaixar livros de Ciência da InformaçãoBaixar livros de Ciência PolíticaBaixar livros de Ciências da SaúdeBaixar livros de ComunicaçãoBaixar livros do Conselho Nacional de Educação - CNEBaixar livros de Defesa civilBaixar livros de DireitoBaixar livros de Direitos humanosBaixar livros de EconomiaBaixar livros de Economia DomésticaBaixar livros de EducaçãoBaixar livros de Educação - TrânsitoBaixar livros de Educação FísicaBaixar livros de Engenharia AeroespacialBaixar livros de FarmáciaBaixar livros de FilosofiaBaixar livros de FísicaBaixar livros de GeociênciasBaixar livros de GeografiaBaixar livros de HistóriaBaixar livros de Línguas
Baixar livros de LiteraturaBaixar livros de Literatura de CordelBaixar livros de Literatura InfantilBaixar livros de MatemáticaBaixar livros de MedicinaBaixar livros de Medicina VeterináriaBaixar livros de Meio AmbienteBaixar livros de MeteorologiaBaixar Monografias e TCCBaixar livros MultidisciplinarBaixar livros de MúsicaBaixar livros de PsicologiaBaixar livros de QuímicaBaixar livros de Saúde ColetivaBaixar livros de Serviço SocialBaixar livros de SociologiaBaixar livros de TeologiaBaixar livros de TrabalhoBaixar livros de Turismo