UNIVERSIDADE FEDERAL FLUMINENSE
TCE - Escola de Engenharia
TEM - Departamento de Engenharia Mecânica
PROJETO DE GRADUAÇÃO II
Título do Projeto :
VIBRAÇÕES EM PÓRTICOS PLANOS PELO MÉTODO
DE ELEMENTOS FINITOS
Autor :
HENRIQUE PEREIRA GÓES
RAFAEL SANCHES RANGEL
Orientador :
PROF.ª ÂNGELA CRISTINA CARDOSO DE SOUZA
Data : 19 de Dezembro de 2013
HENRIQUE PEREIRA GÓES
RAFAEL SANCHES RANGEL
VIBRAÇÕES EM PÓRTICOS PLANOS PELO MÉTODO
DE ELEMENTOS FINITOS
Trabalho de Conclusão de Curso apresentado
ao Curso de Engenharia Mecânica da Universidade
Federal Fluminense, como requisito parcial para
obtenção do grau de Engenheiro Mecânico.
Orientador:
Prof.ª. ÂNGELA CRISTINA CARDOSO DE SOUZA
Niterói
2013
Ficha Catalográfica elaborada pela Biblioteca da Escola de Engenharia e Instituto de Computação da UFF
G598 Góes, Henrique Pereira
Vibrações em pórticos planos pelo método de elementos finitos /
Henrique Pereira Góes, Rafael Sanches Rangel. – Niterói, RJ: [s.n.],
2013.
67 f.
Trabalho (Conclusão de Curso) – Departamento de Engenharia
Mecânica, Universidade Federal Fluminense, 2013.
Orientador: Ângela Cristina Cardoso de Souza.
1. Vibração mecânica. 2. Método de elementos finitos. 3. Pórtico.
I. Rangel, Rafael Sanches. II. Título.
CDD 620.3
UNIVERSIDADE FEDERAL FLUMINENSE
TCE - Escola de Engenharia
TEM - Departamento de Engenharia Mecânica
PROJETO DE GRADUAÇÃO II
AVALIAÇÃO FINAL DO TRABALHO
Título do Trabalho:
VIBRAÇÕES PÓRTICOS PLANOS PELO MÉTODO
DE ELEMENTOS FINITOS
Parecer do Professor Orientador da Disciplina:
- Grau Final recebido pelos Relatórios de Acompanhamento:
- Grau atribuído ao grupo nos Seminários de Progresso:
Parecer do Professor Orientador:
(Comentar a relevância, contribuição e abrangência do trabalho. Se a participação dos alunos no grupo
não se processou de forma homogênea, durante o desenvolvimento do trabalho, compete ao Prof. Orientador
diferenciar o grau de cada aluno, de forma a refletir a sua atuação no desenvolvimento do projeto.)
Nome e assinatura do Prof. Orientador:
Prof.: Ângela Cristina Cardoso de Souza
Assinatura:
Parecer Conclusivo da Banca Examinadora do Trabalho:
Projeto Aprovado sem restrições
Projeto Aprovado com restrições
Prazo concedido para cumprimento das exigências:
/
Discriminação das exigências e/ou observações adicionais:
/
UNIVERSIDADE FEDERAL FLUMINENSE
TCE - Escola de Engenharia
TEM - Departamento de Engenharia Mecânica
PROJETO DE GRADUAÇÃO II
AVALIAÇÃO FINAL DO TRABALHO
(continuação)
Aluno :Henrique Pereira Góes
Grau :
Aluno : Rafael Sanches Rangel
Grau :
Composição da Banca Examinadora :
Prof.: Ângela Cristina Cardoso de Souza
Assinatura :
Prof.: Raul Bernardo Vidal Pessolani
Assinatura :
Prof.: João Marciano Laredo dos Reis
Assinatura :
Data de Defesa do Trabalho : 19/12/2013
Departamento de Engenharia Mecânica,
/
/
4
DEDICATÓRIA
Henrique
Dedico este presente trabalho a Deus, à minha família e às amizades que tive o prazer de
construir ao longo da minha trajetória. Especialmente ao meu pai João, minha mãe Helena e
ao meu irmão Bruno por terem me mostrado o verdadeiro significado de apoio, amor e
companheirismo durante toda minha vida. Amo vocês.
Rafael
Dedico este trabalho à toda minha família e amigos, que sempre me acompanharam e
apoiaram, ajudando com desafios e compartilhando alegrias. Dedico especialmente ao meu
pai Jorge, minha mãe Lúcia, e minhas irmãs Natália e Amanda, que, além de família, são
meus maiores amigos e companheiros. Obrigado por iluminarem o meu caminho sempre.
Amo vocês.
5
AGRADECIMENTOS
Henrique
Agradeço a Deus por ter iluminado meu caminho e me ajudado ao longo de mais uma
grande conquista.
Agradecimento especial à minha família, meu pai João, minha mãe Helena e meu irmão
Bruno por não terem poupado esforços em minha educação e construção de caráter, pelo
exemplo diário de dedicação, luta e persistência, e pelo companheirismo, apoio e amor
incondicional por toda a minha vida.
Ao meu amigo e dupla de projeto, Rafael, pela verdadeira amizade e dedicação ao
trabalho.
A todos meus amigos de longa data e aos que fiz durante os anos de graduação pela
verdadeira amizade, apoio e confiança.
À professora orientadora Ângela Cristina pela dedicação e atenção. Foi fundamental para
construção deste projeto.
A todos os professores que participaram da minha graduação pela instrução e ajuda para
que alcançasse esse sonho de me tornar engenheiro mecânico.
6
AGRADECIMENTOS
Rafael
Agradeço a minha família, pelo seu apoio incondicional, e por terem me ensinado a
valorizar o caráter, a dedicação e o amor acima de tudo. Obrigado por serem meus
companheiros, meus mestres e meu suporte, e pela confiança que tiveram em minha
capacidade. Obrigado por tudo.
Agradeço ao Henrique, amigo e dupla de projeto, pela amizade e pela companhia ao longo
de minha graduação, compartilhando os momentos de descontração e as dificuldades durante
estes anos.
Aos meus amigos, que me ensinaram o valor do companheirismo e da solidariedade.
À professora Ângela Cristina, orientadora deste projeto, pela atenção e dedicação neste
desafio final como aluno.
A todos os meus professores, que foram fundamementais para a minha formação como
pessoa, e como engenheiro mecânico.
7
RESUMO
A proposta do presente trabalho é analisar o modelo bidimensional de vibrações não
amortecidas em vigas, adotando o Método de Elementos Finitos como aproximação
satisfatória dos resultados. Para os diversos estudos, serão abordadas as formas, tanto
diferenciais, quanto integrais de análise dos casos de diferentes carregamentos e condições de
contorno. Uma das principais vertentes do projeto será escolher a melhor forma de discretizar
o modelo de problema escolhido e determinar as matrizes massa e rigidez, assim como o vetor
de carregamento, para aplicação do MEF.
Palavras-Chave: Vibração; Não amortecida; Bidimensional; Método dos Elementos Finitos.
8
ABSTRACT
The aim of the following paper is to analyze the two dimensional, undamped beam vibration
case, adopting the Finite Element Method as a satisfactory result approximation. For that,
both differential and integral form equations will be used to study different loading and
boundary condition scenarios. One of the most important points of this project will be the best
way to discretize the selected problem model and provide the mass and stiffness matrix, as
well as the load vector, in order to apply FEM.
Keywords: Vibration; Undamped vibration; Bidimensional; Finite Element Method.
9
LISTA DE ILUSTRAÇÕES
Figura 1.1 – Sistema Massa-Mola-Amortecedor. Fonte: MENDONÇA (2006) ................................................ 14
Figura 2.1 - Discretização da viga. ...................................................................................................................... 17
Figura 3.1 - Discretização do elemento. .............................................................................................................. 24
Figura 3.2 - Exemplos de bases polinomiais ....................................................................................................... 25
Figura 3.3- Elemento de viga (Vibração Transversal). ....................................................................................... 26
Figura 3.4 - Elemento de viga (Vibração Longitudinal). .................................................................................... 29
Figura 3.5 - Viga biapoiada. ................................................................................................................................ 34
Figura 3.6 - Discretização proposta. .................................................................................................................... 34
Figura 3.7 - Resultado da discretização. .............................................................................................................. 34
Figura 3.8 – Condições de contorno. ................................................................................................................... 38
Figura 3.9 – Primeiros dois modos de vibração para a) viga engastada e b) viga biapoiada. ........................... 41
Figura 4.1- Exemplo transversal. ........................................................................................................................ 42
Figura 4.2 - Resultado estático para deformação transversal (20 elementos por viga). .................................... 43
Figura 4.3 - Ondas de deformação para modos de frequências naturais (Transversal). ................................... 45
Figura 4.4 - Exemplo longitudinal. ..................................................................................................................... 46
Figura 4.5 - Resultado estático para deformação longitudinal (20 elementos por viga). .................................. 47
Figura 4.6 - Ondas de deformação para modos de frequências naturais (Longitudinal). ................................. 49
Figura 4.7 - Carregamentos e condições de contorno (Pórtico). ....................................................................... 50
Figura 4.8 - Resultado estático para o pórtico (20 elementos por viga). ............................................................ 50
Figura 4.9 - Modos de vibração natural do pórtico (20 elementos por viga). .................................................... 51
Figura 4.10 - Carregamentos e condições de contorno (Ponte). ........................................................................ 51
Figura 4.11 - Resultado estático para a ponte metálica (10 elementos por viga). .............................................. 52
Figura 4.12 - Modos de vibração natural da ponte metálica (10 elementos por viga). ...................................... 52
10
LISTA DE TABELAS
Tabela 4.1 – Resultados de Análise Estática (Transversal) ................................................................................ 43
Tabela 4.2–Resultados dos Modos Naturais de Vibração (Transversal – 4 Elem.)............................................ 44
Tabela 4.3 - Resultados dos Modos Naturais de Vibração (Transversal – 20 Elem.) ........................................ 44
Tabela 4.4 - Resultados dos Modos Naturais de Vibração (Transversal – 40 Elem.) ........................................ 45
Tabela 4.5 – Resultados de Análise Estática (Longitudinal) .............................................................................. 46
Tabela 4.6 - Resultados dos Modos Naturais de Vibração (Longitudinal – 4 Elem.) ........................................ 47
Tabela 4.7 - Resultados dos Modos Naturais de Vibração (Longitudinal – 20 Elem.) ...................................... 48
Tabela 4.8 - Resultados dos Modos Naturais de Vibração (Longitudinal – 40 Elem.) ...................................... 48
11
LISTA DE TABELAS
– Força transversal por unidade de comprimento;
– Força longitudinal por unidade de comprimento;
– Deslocamento longitudinal;
w(x, t) - Deslocamento transversal;
– Deslocamento longitudinal;
 - Ângulo com o eixo x;
 - Densidade específica;
- Área de seção transversal;
V – Força cortante;
M – Momento fletor;
 - Tensão normal;
 - Deformação específica;
E - Módulo de Young;
v( x) - Função arbitrária;
I - Momento de inércia;
l - Comprimento do elemento;
u - Parâmetros nodais;
N - Função de forma;
K - Matriz rigidez;
Q - Matriz carregamento;
M - Matriz de massa;
L – Matriz transformada;
 - Autovalores;
 - Autovetores;
 - Frequências naturais.
12
SUMÁRIO
1 INTRODUÇÃO ............................................................................................................................. 13
1.1 OBJETIVO GERAL..................................................................................................................... 15
1.2 OBJETIVO ESPECÍFICO............................................................................................................ 15
1.3 JUSTIFICATIVA ......................................................................................................................... 15
1.4 ESTRUTURA DO TRABALHO ................................................................................................. 16
2 EQUAÇÃO DO MOVIMENTO NO PLANO ............................................................................ 17
2.1 EQUAÇÃO DE EQUILIBRIO NA FORMA DIFERENCIAL .................................................... 17
2.2 EQUAÇÃO DE EQUILIBRIO NA FORMA INTEGRAL .......................................................... 20
3 MÉTODO DOS ELEMENTOS FINITOS .................................................................................. 23
3.1 FUNDAMENTOS........................................................................................................................ 23
3.2 DISCRETIZAÇÃO DA EQUAÇÃO DE EQUILÍBRIO ............................................................. 24
3.2.1 CASO ESTÁTICO ....................................................................................................................... 26
3.2.2 CASO DINÂMICO ...................................................................................................................... 40
4 RESULTADOS .............................................................................................................................. 42
4.1 MOVIMENTO TRANSVERSAL................................................................................................ 42
4.2 MOVIMENTO LONGITUDINAL .............................................................................................. 46
4.3 ANÁLISE DE ESTRUTURAS COMPLEXAS ........................................................................... 49
4.3.1 PÓRTICO ................................................................................................................................... 50
4.3.2 PONTE METÁLICA.................................................................................................................... 51
5 CONCLUSÕES ............................................................................................................................. 53
RERÊNCIAS BIBLIOGRÁFICAS ................................................................................................... 54
6 APÊNDICES .................................................................................................................................. 56
6.1 CÓDIGO MATLAB..................................................................................................................... 57
6.1.1
6.1.2
6.1.3
6.1.4
6.1.5
6.1.6
6.1.7
LIMPAR WORKSPACE .............................................................................................................. 57
DEFINIÇÃO DE PARÂMETROS .................................................................................................. 57
DISCRETIZAÇÃO ....................................................................................................................... 58
MONTAGEM E SOLUÇÃO DO SISTEMA .................................................................................... 59
PLOTAR RESULTADOS ESTÁTICOS ......................................................................................... 61
PLOTAR RESULTADOS DINÂMICOS......................................................................................... 62
ANIMAR MODO DINÂMICO...................................................................................................... 64
13
1 INTRODUÇÃO
A viga, assim como outros elementos mecânicos, é um equipamento estrutural
projetado para suportar efeitos de carregamentos diversos que resultam em esforços tanto de
cisalhamento quanto de flexão. Na esfera das engenharias, essas estruturas devem ser
projetadas de modo a suportar cargas permanentes e dinâmicas geradas por fenômenos
naturais como ventos, marés e terremotos, além de efeitos provocados pelo tráfego de
veículos e pessoas, e operações de máquinas, motores e equipamentos, tornando-se, dessa
forma, peça fundamental para setores como indústria, construção, exploração, transporte,
entre outros. Analisadas como modelos, são capazes gerar estudos a cerca cálculo estrutural
de colunas e lajes, cálculo de flecha máxima, localização de suportes e dimensionamento de
estruturas metálicas.
Relatos informam que o nascimento da teoria da vibração, um dos pilares para o
desenvolvimento da matemática e estudos da mecânica geral, teve seu inicio a partir de
observações dos antigos filósofos gregos. Pitágoras, após verificar as diferentes frequências
sonoras produzidas a partir de forjamento de peças metálicas, abriu os olhos da humanidade
para um efeito que, já na era moderna, foi formalmente modelado por Galileu (Rao, 2008).
O atual contexto de austeridade de recursos naturais e uso eficiente das riquezas
presentes da Terra é uma realidade mundial que direciona os engenheiros modernos a
buscarem constante evolução e avanços tecnológicos de seus projetos. Tal realidade reflete a
tendência do surgimento de novos materiais e técnicas de fabricação e construção capazes de
produzir estruturas cada vez mais esbeltas, leves e eficientes. Nesse sentido, o estudo de
vibrações e o comportamento de sistemas mecânicos submetidos a esse fenômeno se torna
indispensável para projetos de engenharia, pois tais interações geram situações desgastes e
fadiga, reduzindo consideravelmente tanto a vida útil quanto à confiabilidade dos elementos
estruturais.
14
Visando facilitar a interpretação e análise de problemas envolvendo vibrações
mecânicas é comum adotar o modelo de parâmetros concentrados. Tal modelo é caracterizado
principalmente por representar cada elemento através de uma propriedade como: massa,
rigidez ou amortecimento. Transformando, dessa forma, o problema apresentado em um
sistema simplificado conhecido como massa-mola-amortecedor.
Figura 1.1 – Sistema Massa-Mola-Amortecedor. Fonte: MENDONÇA (2006)
Contudo, tratando-se do estudo de um caso particular onde o efeito vibratório não
amortecido é analisado sobre um meio contínuo, os parâmetros massa e rigidezes são
uniformemente distribuídos ao longo do comprimento do objeto de estudo. Dessa forma, a
aplicação do modelo de parâmetros concentrados não forneceria bons resultados.
Embora, no estudo de uma viga, o resultado exato possa ser encontrado através da
obtenção de soluções analíticas para determinadas equações diferenciais, este procedimento
pode se tornar demasiadamente complexo dependendo das condições de contorno e
carregamentos aplicadas ao sistema. Felizmente, acompanhado do desenvolvimento da
capacidade de cálculo computacional, veio o advento de soluções numéricas capazes de
analisar e solucionar sistemas de equações mais elaborados. Nesse contexto, para contornar
essas limitações impostas pelas equações diferenciais, utiliza-se o MEF, que facilita o cálculo
dos resultados ao custo de sua exatidão, reduzindo a complexidade do problema e
aproximando a solução para vários subespaços menores do objeto em estudo, delimitados por
seus nós e chamados elementos. Este método computacional é muito utilizado para analisar o
comportamento estrutural, mecânico, térmico, elétrico ou químico de sistemas quando estes
são descritos por modelos matemáticos complexos. Para a correta análise numérica do
problema em estudo, antes devem ser compreendidos os fundamentos e conceitos do MEF,
garantindo a aplicação do método de maneira correta. Desta maneira, este trabalho terá início
no estudo dos fundamentos e formulações matemáticas do MEF, posteriormente, aplicando-os
à análise dinâmica de vigas.
15
1.1
OBJETIVO GERAL
O objetivo deste trabalho é promover a análise e compreensão do modelo de vibração
em pórticos no plano através das aproximações satisfatórias geradas pelo Método dos
Elementos Finitos.
1.2
OBJETIVO ESPECÍFICO
Visando alcançar o objetivo proposto, pretende-se:
Estudar e determinar as formulações, tanto diferencial quanto integral, das equações
que regem o movimento de vibrações em vigas, tomando como objeto de estudo as vigas de
Euler-Bernoulli.
Aplicar e desenvolver o Método dos Elementos Finitos, buscando analisar as melhores
malhas e matrizes de massa e rigidez.
Buscar, analisar e comparar as soluções analíticas encontradas nos problemas de
vibrações estudados aos numéricos conhecidos.
1.3
JUSTIFICATIVA
A compreensão e o perfeito dimensionamento de estruturas são um dos principais
pilares de qualquer engenharia. A aplicabilidade contempla áreas que vão desde a engenharia
civil, com projetos de estruturas metálicas para pontes, prédios e infraestruturas no geral, até
sistemas de amortecimento de carros e impactos vibratórios em dutos submarinos para
exploração Offshore, no caso de engenheiros mecânicos. Dessa forma, visando uma melhor
otimização, um correto dimensionamento estrutural, e racionalização do uso dos recursos
naturais cada vez mais escassos em nossa sociedade, promove-se o estudo de métodos
satisfatoriamente aproximados aos sistemas reais conhecidos.
16
1.4
ESTRUTURA DO TRABALHO
A estrutura deste presente trabalho seguirá a seguinte ordem de raciocínio: no capítulo
2, será apresentado o objeto de estudo, problema de análise de vibração em vigas
bidimensionais. No capítulo 3, será inserindo uma breve introdução ao Método dos Elementos
Finitos e as premissas adotadas para as análises dos casos estático e dinâmico. No capítulo 4,
encontrar-se-á o quadro dos resultados obtidos e suas respectivas comparações pertinentes
entre os analíticos e numéricos para as diferentes condições de contorno. O capítulo 5
apresenta as conclusões finais do trabalho. Por fim, no capítulo 6, estará disponível a
bibliografia adotada para construção do presente material.
17
2 EQUAÇÃO DO MOVIMENTO NO PLANO
Nesta seção será abordada a natureza de vibração em estruturas de vigas
bidimensionais. Para isso, primeiro, serão introduzidas as equações de equilíbrio, diferenciais
e integrais, aplicadas a um modelo bidimensional.
2.1
EQUAÇÃO DE EQUILIBRIO NA FORMA DIFERENCIAL
Devido ao foco do estudo ser a análise de estruturas bidimensionais, optou-se pelo
modelo geral de análise capaz de inserir um caso complexo onde os diferentes fenômenos de
oscilação, transversal e longitudinal, coincidam.
Figura 2.1 - Discretização da viga.
18
Considerando o diagrama de corpo livre representado pela Figura 2.1, sujeito a um
carregamento distribuído
e
ao longo do comprimento da viga e em movimento
oscilatório transversal representado por
e longitudinal por
.
A força de inércia transversal é dada por:
 2 w  x, t 
 . A( x).dx.
t 2
(2.1)
A partir da aplicação da segunda lei de Newton em y, tem-se:
 V  dV   f  x, t  dx  V  . A( x)dx
 2 w( x, t )
t 2
(2.2)
Da mesma forma, a partir do somatório de momentos no ponto O destacado na Figura
2.1, obtém-se:
 M  dM   V  dV  dx  f  x, t  dx
dx
M 0
2
(2.3)
Podendo, ainda, reescrever os termos como sendo:
dV 
V
M
dx; dM 
dx
x
x
(2.4)
Substituindo as diferenciações acima nas Eqs. (2.2) e (2.3).

V ( x, t )
 2 w( x, t )
 f  x, t    . A( x)
x
t 2
(2.5)
M ( x, t )
 V  x, t   0
x
(2.6)
Aplicando a relação V  M

x
, obtida da Eq. (2.6), tem-se:
 2 M ( x, t )
 2 w( x, t )

f
x
,
t


.
A
(
x
)


x 2
t 2
(2.7)
A partir do modelo de Euller-Bernoulli, a relação entre o momento fletor M ( x, t ) e a
deformação w( x, t ) pode ser expressa como (Beer, 1995):
M  x, t   EI ( x)
 2 w( x, t )
x 2
(2.8)
19
Inserindo a Eq. (2.8) na (2.7), encontra-se a equação do movimento para vibrações
transversais de um corpo não uniforme.

x 2
 2 w  x, t 

 2 w( x, t ) 
 f ( x, t )
 EI ( x) x 2   . A  x  t 2


(2.9)
E, finalmente, assumindo uniformidade das características ao longo de todo o
comprimento da viga.
EI
 2 w  x, t 
 4 w( x, t )


.
A
x
 f ( x, t )
 
x 4
t 2
(2.10)
A força de inércia longitudinal é dada por:
 . A( x).dx.
 2 u  x, t 
t 2
(2.11)
Da mesma forma, aplicando da segunda lei de Newton em x, tem-se:
 2u  x, t 
 P  dP   q  x, t  dx  P  .A.dx.
t 2
(2.12)
Podendo, ainda, reescrever o termo como sendo:
dP 
P
dx
x
(2.13)
Substituindo as diferenciações acima na Eq. (2.12).
 2 u  x, t 
P
 q( x, t )  . A
x
t 2
(2.14)
Contudo, sabe-se que pela Lei de Hooke.
  E.
(2.15)
Mas,

P
u ( x, t )
;
A
x
Dessa forma, podemos reescrever a Eq.(2.15) como:
P  E. A
u ( x, t )
x
(2.16)
20
Substituindo a relação acima na Eq.(2.14), encontra-se a equação do movimento para
vibrações longitudinais de um corpo não uniforme.

 2u  x, t 

u ( x, t ) 
E
.
A
(
x
)


.
A
 q ( x, t )


x 
x 
t 2
(2.17)
Finalmente, assumindo uniformidade ao longo do corpo.
.A
2.2
 2 u  x, t 
 2u ( x, t )

E
.
A
 q( x, t )
t 2
x 2
(2.18)
EQUAÇÃO DE EQUILIBRIO NA FORMA INTEGRAL
Seguindo o mesmo raciocínio proposto na seção anterior, serão apresentadas as
formulações integrais aplicadas para o caso de vibração geral do corpo, considerando a
transversal e, em seguida, longitudinal.
Diferentemente das equações diferencias, que são consideradas formulações fortes e
de soluções exatas difíceis e limitadas a específicos, a integral (fraca) permite a aplicação de
um método único para resolver problemas físicos distintos.
Deve-se ressaltar a importância da correta formulação das equações integrais que
regem os movimentos (transversais e longitudinais), pois será a partir destas equações que,
posteriormente, obter-se-ão as matrizes, massa e rigidez, utilizadas no estudo do MEF.
Seja a equação diferencial do movimento oscilatório de vibrações transversais (Eq.
(2.10)) e a função peso v( x) definida como suficientemente regular. Ou seja, possui derivada
continua e finita no intervalo que está contida e satisfaz as seguintes condições:

v  x   0 , para pontos onde o deslocamento é definido.

dv( x)
dx
 0 , para os pontos onde a rotação é definida.
Aplica-se o principio do trabalho virtual (PVT), integrando a eq.(2.10) e multiplicando
por v  x  Dessa forma, tem-se:
l
l
 3 w  x, t  
 2 w  x, t 
 
EI
x
.
v
x
dx


.
A
x
.v  x  dx  0 (2.19)








0
x 
x3 
t 2
0
l
 f  x, t  .v  x  dx  
0
21
Para,
v  x   U ; EI  x 
 3 w  x, t 
V
x3
Como a parcela w  x, t  não varia com o tempo, a derivada parcial se torna derivada
total e desmembrando a integral referente à parcela elástica da equação, integrando-a por
partes. Assim,
d 3w  x  
d 
EI
x


 .v  x  dx
0 dx 
dx3 
l
l
d 3 w  x  dv( x)

d 3 w( x)  l
 v  x  .EI ( x)
 EI  x 
.
dx
dx3  0 0
dx3
dx

(2.20)
Integrando mais uma vez por partes.
d 3w  x  
d 
0 dx  EI  x  dx3  .v  x  dx
l
d 2 w  x  dv( x)  l

d 3 w( x)  l 
 v  x  .EI ( x)
  EI  x 
.

dx3  0 
dx 2
dx  0

(2.21)
d 2 w  x  dv 2 ( x)
 EI  x 
.
dx
dx 2
dx 2
0
l
E, finalmente, obtêm-se a formulação fraca.
l
EI  x 
0
d 2 w  x  dv 2 ( x)
.
dx 
dx 2
dx 2
l
l
0
0
 f  x  .v  x  dx  .A  x 
d 2w  x
.v  x  dx
dt 2
(2.22)
d 2 w  x  dv( x)  l

d 3 w( x)  l 
 v  x  .EI ( x)
  EI  x 
.

dx3  0 
dx 2
dx  0

Seguindo o mesmo raciocínio descrito acima, aplicando o princípio do trabalho virtual
de D’Alembert a equação diferencial do movimento oscilatório de vibrações longitudinais,
multiplicando-a pela função peso v( x) , suficientemente regular. Tem-se:
l
.A
0
l
l
 2 u  x, t 
 2 u  x, t 
v
(
x
)
dx

E
.
A
v
(
x
)

0
0q  x  .v  x  dx  0
t 2
x 2
(2.23)
22
Para,
v  x   U ; E. A
u
V
x
Desmembrando a integral referente à parcela elástica da equação, integrando-a por
partes, e lembrando que como a parcela u  x, t  não varia com o tempo, a derivada parcial se
torna derivada total. Assim,
l
E.A
0
l
d 2u  x 
du  l
du dv

v
(
x
)

v
x
E
.
A

E. A
dx


2



dx
dx  0 0
dx dx

(2.24)
Retornando a relação acima na Eq. (2.23), finalmente, obtêm-se a formulação fraca.
l
l
du  x  dv( x)
d 2u  x 
du

l
0 E.A dx dx dx  0q  x  .v  x  dx  0 .A dt 2 v( x)dx   E.A dx v  x  0
l
(2.25)
23
3 MÉTODO DOS ELEMENTOS FINITOS
Em inúmeros casos de problemas de engenharia é possível encontrar a solução
desejada simplesmente resolvendo o sistema de equações diferenciais (formulação forte) que
regem a situação analisada. Contudo, a solução analítica dessas equações pode se tornar
demasiadamente complexa e particular para o caso estudado. Nesse contexto, surgiu a
necessidade de se desenvolver um método numérico para que, apesar de se distanciar
parcialmente da solução analítica, seja capaz de resolver problemas envolvendo situações de
engenharia mais complexas.
A seção a seguir será dedicada à aplicação do Método dos Elementos Finitos (MEF),
uma das formas de resolução de problemas complexos envolvendo a formulação fraca.
3.1
FUNDAMENTOS
Originalmente desenvolvido para análise de casos envolvendo o setor aeronáutico e
da construção civil, o MEF, atualmente, tornou-se uma ferramenta indispensável para projeto
e concepção de estruturas contínuas inovadoras e arrojadas nas mais variadas áreas, como:
estruturas oceânicas e navios, veículos rodoviários e ferroviários, hidro geradores, estruturas
aeroespaciais e aviões, mecânica estrutural, mecânica dos fluídos, condução de calor,
eletromagnetismo, entre outros (Silva, 2009).
Apesar de flexível quanto à aplicabilidade, este método computacional é bastante
restrito em termo de metodologia de aplicação. A ordem de execução do MEF segue a
seguinte linha de execução:
1. Desenvolvimento das equações do problema;
24
2. Discretização do domínio contínuo de soluções na forma de uma malha de
elementos finitos;
3. Montagem das equações do elemento;
4. Introdução das condições de contorno (restrições físicas e geométricas);
5. Solução para nós desconhecidos;
6. Cálculo da solução e quantidades (grandezas) em cada elemento.
3.2
DISCRETIZAÇÃO DA EQUAÇÃO DE EQUILÍBRIO
Têm-se, como termos fundamentais da discretização de um domínio, o elemento, os
nós, que delimitam as fronteiras deste elemento, e os respectivos graus de liberdade aplicados
a cada nó. No caso do estudo de vibrações bidimensionais, o elemento estudado apresentará
dois nós em suas extremidades e, atribuído a cada nó, três graus de liberdade conforme
ilustrado na Figura 3.1.
Figura 3.1 - Discretização do elemento.
O racional por trás do MEF é a divisão do domínio de integração de uma estrutura em
um conjunto de pequenos segmentos reduzidos. Tais subdomínios são denominados de
elementos finitos e são responsáveis pela transformação do domínio contínuo em discreto. O
conjunto dos elementos resultantes da discretização é denominado malha ou grid (Silva,
2009).
Para cada elemento, existe uma função capaz de interpolar as funções de aproximação
e as funções de peso. Tal função é conhecida como função de forma ou de interpolação.
25
Abaixo, na Figura 3.2, existem algumas formas que a função de interpolação pode assumir
dependendo da base polinomial adotada.
Figura 3.2 - Exemplos de bases polinomiais
Dessa forma, escolhem-se funções de forma que vão determinar a influência dos
parâmetros nodais (soluções para cada nó) sobre o resto do domínio. Para isso, são utilizadas
aproximações do tipo:
 u1 
 
u ( x)  N j  x  u j   N1 N 2   u2 
j 1
 
 
(3.1)
 v1 
 
v ( x)  N j  x  v j   N1 N 2   v2 
i 1
 
 
(3.2)
n
'
n
'
Onde cada função N i vale 1 no nó i , e decresce conforme se distancia deste nó. É
através do produto destas funções de forma com os parâmetros nodais, então, que pode se
determinar a influência de um nó sobre outras regiões do domínio, e por meio da soma da
influência de todos os nós em um determinado ponto, encontra-se a solução aproximada para
este ponto.
Deve-se destacar que as funções de forma são de suma importância para determinação
das submatrizes massa e rigidez.
26
3.2.1 Caso Estático
Seguindo o mesmo raciocínio proposto na seção 2 deste material, como o foco do
estudo trata da análise de vibração bidimensional em vigas, será proposto à aplicação do MEF
aos casos de oscilação transversal e longitudinal para, em seguida, gerar matrizes massa e
rigidez capazes de descrever o fenômeno físico de vibração.
3.2.1.1 Análise de Vibração Transversal
Seja o elemento representado pela Figura 3.3, no caso da discretização de um
elemento unidimensional de uma viga. A cada nó, serão atribuídos dois graus de liberdade.
Enquanto as grandezas representadas por
e
se referem aos deslocamentos verticais,
e
dizem respeito aos angulares.
Figura 3.3- Elemento de viga (Vibração Transversal).
Dessa forma, pode-se representar o vetor deslocamento deste vetor como sendo:
 u1   w1 
   
u

e
w  2 1 
 u3   w2 
   
 u4    2 
Para o caso de vibração transversal, os polinômios hermetianos cúbicos serão adotados
como funções de forma. Tal escolha se fez pela característica destes polinômios em
considerar tanto os deslocamentos verticais, quanto os angulares. Assim, aplicados a um
elemento de comprimento , tem-se o vetor função de forma como sendo (Fish, 2007):
N   N1 N2 N3 N4 
(3.3)
27
Onde,
2
x
 x
N1  1  3    2  
l
l
2
3
 x
 x
N 2  x  2l    h  
l
l
2
 x
 x
N3  3    2  
l
l
2
3
3
 x
 x
N 4  l    l  
l
l
(3.4)
3
Utilizando a formulação integral referente à vibração transversal verificada na seção 2
deste material.
d 2 w  x  dv 2 ( x)
0 EI  x  dx2 . dx2 dx 
l
d 2w x
0 f  x  .v  x  dx  0 .A  x  dt 2 .v  x  dx
l
l
(3.5)
d 2 w  x  dv( x)  l

d 3 w( x)  l 
 v  x  .EI ( x)

EI
x
 0    dx 2 . dx  0
3
dx




E substituindo os valores representados pelas Eqs. (3.1) e (3.2) na parcela elástica,
tem-se:
l
EI  x 
0
d 2 w  x  dv 2 ( x)
.
dx
dx 2
dx 2
l
l
0
0
 EI  N ".w)( N ".v  dx  EI  wT .N "T )( N .v  dx
(3.6)
 l T

  EI N " N " dx  wT v
 0

Integrando a relação acima, obtêm-se a submatriz rigidez do elemento de EullerBernoulli estudado. Esta parcela é um parâmetro conservativo e responsável pela preservação
da energia potencial devido à deformação.
28
 N "1 N "1 N "2 N "2 N "3 N "1 N "1 N "4 


N " N " N "2 N "2 N "3 N "2 N "2 N "4 
K e  EI   1 2
dx


N
"
N
"
N
"
N
"
N
"
N
"
N
"
N
"
1
3
2
3
3
3
3
4
0


 N "1 N "4 N "2 N "4 N "3 N "4 N "4 N "4 
6
6 
12
 12
l²
l³
l² 
 l³
 6
6
4
2 
l²
l
l²
l 
 EI 
 12
6
6 
12

l³
l²
l³
l² 
 6
6
2
4 

l
l²
l 
 l²
l
(3.7)
Adotando o mesmo racional e isolando a parcela da inércia da Eq.(3.5), tem-se:
l
.A  x 
0
l
l


d 2w x
.
v
x
dx


.
A
N
.
w
")(
N
.
v
dx


.
A
N T Ndx  w "T v
 



2


dt
0
0


(3.8)
Realizando a integração do termo acima, obtém-se a matriz massa do elemento da viga
de Euller-Bernoulli. Esta parcela está relacionada a características de corpo rígido capaz de
armazenar energia cinética.
 N1 N1

NN
e
M   . A  1 2
 N1 N3
0

 N1 N 4
l
N 2 N1
N3 N1
N2 N2
N 2 N3
N2 N4
N3 N 2
N3 N3
N3 N 4
N1 N 4 
54 13l 
 156 22l



N2 N4 
 . A.l  22l 4l ² 13l 3l ² 
dx 
(3.9)
N3 N 4 
420  54 13l 156 22l 



N4 N4 
 13l 3l ² 22l 4l ² 
Para o carregamento aplicado sobre o elemento, deve-se determinar o vetor
carregamento. Tal vetor deve ser constituído pelo carregamento nodal e o carregamento
distribuído. Dessa forma, a partir da parcela da contribuição do carregamento na Eq.(3.5),
tem-se como vetor carregamento distribuído:
l
l
 l T  T
T
T
f
x
v
x
dx

f
N
.
v
dx

f
v
.
N
dx





 f N dx  v
0
o
0
 0

l
(3.10)
Obtendo-se, através da integração, o vetor carregamento distribuído do elemento da
viga.
 N1 
6


 
l
N2 
f .l  l 
e

Q  f
dx 
 N3 
12  6 
0
 
 
 l 
 N4 
(3.11)
29

d 3 w( x)  l
Finalizando, ainda restam os termos da força v  x  .EI ( x)
e o momento
dx3  0


d 2 w  x  dv( x)  l
EI
x
.
  
 aplicados em cada nó. Estes termos, se existirem, darão origem ao
2
dx
dx

0
vetor carregamento nodal que será definido posteriormente.
3.2.1.2 Análise de Vibração Longitudinal
Seja o elemento representado pela Figura 3.4, no caso da discretização de um
elemento unidimensional de uma viga, a cada nó, será atribuído um grau de liberdade.
Figura 3.4 - Elemento de viga (Vibração Longitudinal).
Dessa forma, pode-se representar o vetor deslocamento deste vetor como sendo
unicamente:
u   z 
ze   1    1 
 u2   z 2 
(3.12)
Diferentemente, da premissa adotada para o caso de oscilação transversal, para o caso
de vibração longitudinal, devido à relativa simplicidade de graus de liberdade, será utilizada a
forma linear da função de interpolação. Dessa forma, tem-se (Ferreira, 2008):
1 , para x  x j
Nj  
0 , para x  xi  i  j 
(3.13)
Assim, pelo critério da linearidade.
N1  x  
x2  x
 x
N1  x   1  
l
 l
(3.14)
x  x1
 x
N 2  x    
l
l
(3.15)
N2  x  
Utilizando a formulação integral referente à vibração longitudinal verificada na seção
2 deste material.
30
l
E.A
0
l
l
du  x  dv( x)
d 2u  x 
du

l
dx  q  x  .v  x  dx   . A
v( x)dx   E. A v  x  
2
dx
dx
dt
dx

0
0
0
E substituindo os valores representados pelas Eqs. (3.14) e (3.15) na parcela elástica,
tem-se:
l
E . A
0
du  x  dv( x)
dx
dx
dx
l
 E. A  N ' .z  N ' .v  dx
(3.16)
0
l


 E. A  z .N  N .v  dx   E. AN 'T .N ' dx  z T .v
0
0


l
T
'T
'
Integrando a relação acima, obtêm-se a submatriz rigidez do elemento estudado.
l
N ' N '
K e  E . A  1 1
N '2 N '1
0 
N '1 N '2 
E. A  1 1
dx 

N '2 N '2 
l  1 1 
(3.17)
Adotando o mesmo racional e isolando a parcela da inércia da Eq.(2.25), tem-se:
l
 . A
0
d 2z  x
v( x)dx
dt 2
l
  . A  N .z  N .v  dx
(3.18)
0
l


  . A  z .N   N .v  dx    . A  N T .N  dx  z T .v
0
0


l
T
T
Realizando a integração do termo acima, obtém-se a matriz massa do elemento da
viga. Esta parcela está relacionada a características de corpo rígido capaz de armazenar
energia cinética.
N N
M   . A  1 1
N N
0 2 1
l
e
l
N1 N 2 
3
dx   . A 

l
N2 N2 
 6
l 
6
l 
3
(3.19)
Isolando-se a parcela do carregamento aplicado, tem-se uma formulação semelhante
ao caso já estudado de vibração transversal:
l
l
 l T  T
T
T
q
x
v
x
dx

q
N
.
v
dx

q
v
.
N
dx





q N dx  v
0
o
0
 0

l
Contudo, levando-se em conta a formulação linear:
(3.20)
31
l
N 
q.l 6
Q e  q   1  dx 
N
12 6
0  2
(3.21)
du

l
Finalizando, ainda resta o termo da força  E. A v  x   aplicada em cada nó. Este
dx

0
termo, se existir, dará origem ao vetor carregamento nodal que será definido posteriormente.
3.2.1.3 Montagem do Sistema Global
Ao longo do estudo apresentado até o momento, foi proposta a análise distinta dos
casos de vibrações transversais e longitudinais. Contudo, para a montagem do sistema global
de equações, as matrizes locais dos elementos (rigidez, massa e carregamento) deverão ser
agrupadas em apenas uma por natureza. Dessa forma, respeitando a ordenação dos graus de
liberdade proposta, representada pela Figura 3.1, tem-se o vetor carregamento sendo:
De  u1 u2
u3
u4
u5
u6   u1
T
w1 1 u2
w2 2 
T
(3.22)
Totalizando, assim, um vetor 6x1 caracterizado pelos seus três graus de liberdade
aplicados aos dois nós por elemento. Consequentemente, as matrizes locais de massa, rigidez
e carregamento referentes ao vetor deslocamento apresentado também sofrerão ajustes,
tornando-se resultado da soma, respeitando os locais de inserção de valores a partir da
ordenação dos graus de liberdade, das matrizes encontradas nos casos particulares de
vibrações transversais e longitudinais. Dessa forma, tem-se:
e
e
K e  Ktransversal
 Klongitudinal
32
0
0

0 12 E.I l ²

0 6 E.I l ²
Ke  
0
0
 12 E.I
0
l²

6
E
.
I
0
l²
 E. A
l

 0

0

  E. A

l
 0

 0

  u1 
l²
l²
l ²  u 
 2
4 E.I
2 E.I   
0 6 E.I
u3
l
l²
l
. 
0
0
0
0  u4 
6 E.I
6 E.I  u5 
0 12 E.I
l²
l²
l²  
u
2 E.I
4 E.I   6 
0 6 E.I
l
l²
l 
0 0  E. A
0 0  u 
l
 1
0 0
0
0 0  u2 

0 0
0
0 0  u3 
. 
0 0 E. A
0 0  u4 
l
u 
0 0
0
0 0  5 
 u
0 0
0
0 0   6 
0
6 E.I
0
0
0 12 E.I
0
6 E.I
Tendo, dessa forma, a submatriz rigidez para um sistema de vibração bidimensional.
 E. A
l

 0

 0

e
K 
 E. A

l

 0

 0
0
12 E.I
6 E .I
l²
l²
0
6 E .I
l²
4 E .I
l
l²
l²
6 E.I
2 E .I
l²
0
0
l
0
0


12 E.I
6 E.I   u1 
l²
l ²  u 
2
6 E.I
2 E .I   
l²
l   u3 
. 
0
0  u4 
 u 
12 E.I
6 E.I   5 
l²
l ² u6 

6 E.I
4 E.I 
l²
l 
0
l
 E. A
0
12 E.I
6 E .I
E. A
0
l
0
(3.23)
Seguindo o mesmo raciocínio para submatriz massa.
e
e
M e  M transversal
 M longitudinal
0
0
0
0 156 22l

 . A.l 0 22l 4l ²
Me 

0
0
420 0
0 54 13l

0 13l 3l
0
0
0 54
0 13l
0
0
0 156
0 22l
0   u1 
140



 0
13l  u2 

3l ²  u3   . A.l  0
.  

0  u4  420  70
 0
22l  u5 
  

4l ²  u6 
 0
0
0
0
0
0
0
0 70
0 0
0 0
0 140
0 0
0 0
0
0
0
0
0
0
0   u1 
0  u2 
0  u3 
 .   (3.24)
0  u4 
0  u5 
  
0  u6 
Tendo, dessa forma, a submatriz massa para um sistema de vibração bidimensional.
33
0
0
70
0
140
 0
156 22l 0
54

22l 4l ² 0
13l
 . A.l  0
Me 

0
0 140
0
420  70
 0
54 13l
0
156

 0 13l 3l 0 22l
0   u1 
13l  u2 
3l ²  u3 
. 
0  u 4 
22l  u5 
  
4l ²  u6 
(3.25)
Para a submatriz carregamento aplicado, uma consideração distinta será feita. Supondo
uma condição geral onde existe tanto carregamento distribuído em x quanto em y e que ambos
exercem influência sobre o momento aplicado. Tem-se

Q  1 q
2

e
Onde a parcela
f
f
 q l
6
q
f
f
 ql 

6

T
(3.26)
se refere ao carregamento que interfere na vibração longitudinal e o
, em contrapartida, na transversal.
Uma vez tendo disponíveis as submatrizes massa, rigidez e carregamento, o próximo
passo para execução do MEF será a obtenção das matrizes globais. Lembrando que as
submatrizes são originadas a partir de elementos finitos do todo, para se determinar as
matrizes globais, basta combinar suas respectivas submatrizes. Quanto maior for o número de
elementos adotados na discretização, maior
será o número de submatrizes e,
consequentemente, mais complexa a matriz global.
Portanto, a definição da malha de discretização é uma etapa importante na solução de
problemas envolvendo MEF. No geral, fica a critério da experiência do analista determinar
quantos elementos utilizar, porém a discretização deve satisfazer a existência de um nó nos
seguintes casos:

Na extremidade da estrutura a ser analisada;

Onde existirem condições de contorno conhecidas;

Onde existir aplicação de força ou momento;

Onde exista alteração de área de seção que influencie no momento de inércia;

Onde exista alteração de material que influencie o módulo de elasticidade, E, e
a densidade,
Como exemplo de montagem de um sistema global, considere a estrutura representada
pela Figura 3.5.
34
Figura 3.5 - Viga biapoiada.
Seguindo os critérios mencionados para determinação dos nós. Supondo um total de
três elementos, tem-se a seguinte proposta de discretização.
Figura 3.6 - Discretização proposta.
Tendo, como resultado da discretização:
Figura 3.7 - Resultado da discretização.
Como o objeto de estudo é uniforme, de propriedades constantes, e, durante a
discretização, optou-se por elementos finitos de mesmo comprimento. Tem-se que as
submatrizes
,
e
não se alteram de elemento para elemento, permanecendo idênticas
às apresentadas nas Eqs. (3.23), (3.25) e (3.26).
Para que a combinação das submatrizes seja feita, deve-se ter em mente que os
elementos vizinhos sempre compartilharão um nó de ligação entre eles. Dessa forma, os
deslocamentos nodais,
,
e , destes nós se repetirão nos elementos ligados.
35
Considerando a discretização proposta acima, tem-se como matriz rigidez global:
K
 K11a
 a
 K 21
 K 31a
 a
 K 41
K a
 51a
 K 61
 0

 0

 0
 0

 0
 0

K12a
K 22a
K 32a
K13a
K 23a
K 33a
K14a
K 24a
K 34a
K15a
K 25a
K 35a
K16a
K 26a
K 36a
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
K 42a
K 52a
K 62a
K 43a
K 53a
K 63a
K 44a  K11b
b
K 54a  K 21
K 64a  K 31b
K 45a  K12b
b
K 55a  K 22
K 65a  K 32b
K 46a  K13b
b
K 56a  K 23
K 66a  K 33b
K14b
b
K 24
K 34b
K15b
b
K 25
K 35b
K16b
b
K 26
K 36b
0
0
0
0
0
0
0
0
0
0
0
0
b
K 41
K 51b
b
K 61
b
K 42
K 52b
b
K 62
b
K 43
K 53b
b
K 63
b
K 44
 K11c
c
K 54b  K 21
b
K 64
 K 31c
b
K 45
 K12c
c
K 55b  K 22
b
K 65
 K 32c
b
K 46
 K13c
c
K 56b  K 23
b
K 66
 K 33c
K14c
c
K 24
K 34c
K15c
c
K 25
K 35c
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
c
K 41
K 51c
K 61c
c
K 42
K 52c
K 62c
c
K 43
K 53c
K 63c
c
K 44
K 54c
K 64c
c
K 45
K 55c
K 65c





0 
0  (3.27)

0 
K16c 

c

K 26
c 
K 36 
c

K 46
c 
K 56 
K 66c 
0
0
0
Seguindo o mesmo raciocínio para a matriz massa global e para o vetor carregamento.
M
 M 11a
 a
 M 21
 M 31a
 a
 M 41
M a
 51a
 M 61
 0

 0

 0
 0

 0
 0

M 12a
M 22a
M 32a
M 13a
M 23a
M 33a
M 14a
M 24a
M 34a
M 15a
M 25a
M 35a
M 16a
M 26a
M 36a
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
M 42a
M 52a
M 62a
M 43a
M 53a
M 63a
M 44a  M 11b
b
M 54a  M 21
M 64a  M 31b
M 45a  M 12b
b
M 55a  M 22
M 65a  M 32b
M 46a  M 13b
b
M 56a  M 23
M 66a  M 33b
M 14b
b
M 24
M 34b
M 15b
b
M 25
M 35b
M 16b
b
M 26
M 36b
0
0
0
0
0
0
0
0
0
0
0
0
b
M 41
M 51b
b
M 61
b
M 42
M 52b
b
M 62
b
M 43
M 53b
b
M 63
b
M 44
 M 11c
b
c
M 54  M 21
b
M 64
 M 31c
b
M 45
 M 12c
b
c
M 55  M 22
b
M 65
 M 32c
b
M 46
 M 13c
b
c
M 56  M 23
b
M 66
 M 33c
M 14c
c
M 24
M 34c
M 15c
c
M 25
M 35c
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
c
M 41
M 51c
M 61c
c
M 42
M 52c
M 62c
c
M 43
M 53c
M 63c
c
M 44
M 54c
M 64c
c
M 45
M 55c
M 65c





0 
(3.28)
0 

0 
M 16c 

c

M 26

M 36c 
c

M 46
c 
M 56 
M 66c 
0
0
0
36
 Q1a 


a
 Q2 
 Q3a 
 a
b
Q4  Q1 
Q a  Q b 
2
 5a

b
Q6  Q3 
Q b
Q  Q1c 
 4

 Q5b  Q2c 
 b
c
 Q6  Q3 
 Q4c 


c
 Q5 
 Qc 
6


(3.29)
Restando, apenas, o vetor carregamento nodal. Neste serão representadas as cargas que
estão sendo aplicadas em cada nó. No caso de exemplo estudado,
0
0
 
0
 
0
P
 
0
P   
0
 
0
0
 
0
 
0
 M 
no nó 2 e
no 4.
(3.30)
Os vetores carregamento distribuído (3.29) e nodal (3.30) constituem o vetor
carregamento, F.
37
 Q1a   0 

  
a
 Q2   0 
 Q3a   0 
 a
 
b
Q4  Q1   0 
Q a  Q b   P 
2
 5a
  
b
Q6  Q3   0 
F QP   b

Q4  Q1c   0 

  
 Q5b  Q2c   0 
 b
 
c
 Q6  Q3   0 
 Q4c   0 

  
c
 Q5   0 
 Q c   M 
6


(3.31)
Deve-se ressaltar que, devido ao estudo envolvido considerar vigas de orientações
diferentes, prováveis desalinhamentos entre as coordenadas locais de deformação (transversal,
longitudinal, e angular) e globais devem ocorrer. Na Figura 3.1, onde estão representados os
vetores locais e globais de deslocamento, é possível perceber tal desalinhamento angular.
Visando corrigir esta inclinação entre os dois sistemas de coordenadas, será inserida, no
cálculo das submatrizes, a matriz transformada, , capaz de projetar as características locais
encontradas sobre um único referencial global estipulado. Sendo os vetores , global, e
,
local, tem-se (Ferreira, 2008):
uT  u1 u2
u3
u4
u 'T  u '1 u '2
u '3
u '4
u5
u '5
u6 
u '6 
Sendo,
 cos( ) sin( )
  sin( ) cos ( )

 0
0
L
0
 0
 0
0

0
 0
0
0 
0

0
0  sin( ) cos( ) 0 

0
0
0
1
0
0
1
0
0
0
0
cos( )
0
0
0
sin( )
(3.32)
38
Logo, aplicando a transformada nos vetores e matrizes locais, tem-se:
u  Lu
K  LT K  L
M  LT M  L
3.2.1.4 Condições de Contorno
As condições de contorno são particulares para cada problema. São tais características
que determinam os deslocamentos nodais que já são conhecidos em determinado momento e,
sem as quais, a resolução do sistema de equações seria demasiadamente complexa e
trabalhosa ou, simplesmente, impossível.
As condições gerais de contorno presentes em análise de vigas, representadas pela
Figura 3.8 são basicamente: Engaste, Apoio e Apoio com rolamento.
Figura 3.8 – Condições de contorno.
Na prática, a aplicação das condições de contorno em sistemas de equações é bastante
simples. Basta, para os graus de liberdade afetados, retirar o conjunto de interações que
resultaria em seu deslocamento já conhecido. No caso do exemplo estudado de vibração
transversal representado pela Figura 3.5, como a viga em questão está engastada na
extremidade esquerda
e apoiada na direita
, seus respectivos
deslocamentos dos graus de liberdade serão nulos. Assim, podem-se reescrever as matrizes
globais de massa, rigidez e carregamento como:
39
K
 K11a
 a
 K 21
 K 31a
 a
 K 41
K a
 51a
 K 61
 0

 0

 0
 0

 0
 0

K12a
K 22a
K 32a
K13a
K 23a
K 33a
K14a
K 24a
K 34a
K15a
K 25a
K 35a
K16a
K 26a
K 36a
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
K 42a
K 52a
K 62a
K 43a
K 53a
K 63a
K 44a  K11b
b
K 54a  K 21
K 64a  K 31b
K 45a  K12b
b
K 55a  K 22
K 65a  K 32b
K 46a  K13b
b
K 56a  K 23
K 66a  K 33b
K14b
b
K 24
K 34b
K15b
b
K 25
K 35b
K16b
b
K 26
K 36b
0
0
0
0
0
0
0
0
0
0
0
0
b
K 41
K 51b
b
K 61
b
K 42
K 52b
b
K 62
b
K 43
K 53b
b
K 63
b
K 44
 K11c
c
K 54b  K 21
b
K 64
 K 31c
b
K 45
 K12c
c
K 55b  K 22
b
K 65
 K 32c
b
K 46
 K13c
c
K 56b  K 23
b
K 66
 K 33c
K14c
c
K 24
K 34c
K15c
c
K 25
K 35c
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
c
K 41
K 51c
K 61c
c
K 42
K 52c
K 62c
c
K 43
K 53c
K 63c
c
K 44
K 54c
K 64c
c
K 45
K 55c
K 65c
M 12a
M 22a
M 32a
M 13a
M 23a
M 33a
M 14a
M 24a
M 34a
M 15a
M 25a
M 35a
M 16a
M 26a
M 36a
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
M 42a
M 52a
M 62a
M 43a
M 53a
M 63a
M 44a  M 11b
b
M 54a  M 21
M 64a  M 31b
M 45a  M 12b
b
M 55a  M 22
M 65a  M 32b
M 46a  M 13b
b
M 56a  M 23
M 66a  M 33b
M 14b
b
M 24
M 34b
M 15b
b
M 25
M 35b
M 16b
b
M 26
M 36b
0
0
0
0
0
0
0
0
0
0
0
0
b
M 41
M 51b
b
M 61
b
M 42
M 52b
b
M 62
b
M 43
M 53b
b
M 63
b
M 44
 M 11c
c
M 54b  M 21
b
M 64
 M 31c
b
M 45
 M 12c
c
M 55b  M 22
b
M 65
 M 32c
b
M 46
 M 13c
c
M 56b  M 23
b
M 66
 M 33c
M 14c
c
M 24
M 34c
M 15c
c
M 25
M 35c
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
c
M 41
M 51c
M 61c
c
M 42
M 52c
M 62c
c
M 43
M 53c
M 63c
c
M 44
M 54c
M 64c
c
M 45
M 55c
M 65c





0 
0 

0 
K16c 

c

K 26
c 
K 36 
c

K 46
c 
K 56 
K 66c 
0
0
0
M
 M 11a
 a
 M 21
 M 31a
 a
 M 41
M a
 51a
 M 61
 0

 0

 0
 0

 0
 0



Q1a


a
Q2


a


Q3


a
b
 Q4  Q1 
Q a  Q b  P 
 5 a 2 b 
 Q  Q3 
F   6b
Q  Q1c 
 4

 Q5b  Q2c 


b
c
 Q6  Q3 


Q4c


Q5c


 Qc  M 
6







0 
0 

0 
M 16c 

c

M 26
c 
M 36 
c

M 46
c 
M 56 
M 66c 
0
0
0
40
Podendo, uma vez reduzido o sistemas de equações através da aplicação das condições
de contorno, montar a equação que representa o problema de vibrações não amorcecidas em
vigas.
 M .u   K .u  P(t )
(3.33)
3.2.2 Caso Dinâmico
O estudo e compreensão dos casos de vibração envolvendo carregamentos dinâmicos
são de suma importância pra o projeto e concepção de estruturas, sejam elas: prédios, pontes
ou tubulações. Fazendo-se esta análise, é possível determinar os modos de vibração de um
sistema e suas respectivas frequências naturais.
A determinação das frequências naturais de vibração é justificada para se evitar o
efeito indesejável da ressonância. Este fenômeno é desencadeado toda vez que determinada
estrutura se sujeita a carregamentos cíclicos de frequência aproximada às naturais particulares
dos modos de vibração. O efeito da ressonância confere a estrutura um aumento consistente e
gradual das amplitudes de oscilação, podendo levar o conjunto a um colapso estrutural não
previsto pelos estudos estáticos.
Tendo em vista a importância da determinação das frequências naturais de vibração, é
possível as determinar a partir da Equação de Frequência do Sistema.
det  K    ²  M   0
Onde
(3.34)
é a frequência natural em rad/s. As raízes da equação acima são denominadas
de autovalores, λ, para os quais a relação
é válida. Dessa forma, reescrevendo a Eq.
(3.34).
det  K     M   0
Em um sistema com
consequentemente,
graus de liberdade, pode-se determinar
(3.35)
modos de vibração e,
frequências naturais. Para cada autovalor,  , existe um autovetor,  ,
associado, característico de um modo de vibração e determinado a partir da relação:
 K     M    0
i
i
(3.36)
41
Figura 3.9 – Primeiros dois modos de vibração para a) viga engastada e b) viga
biapoiada.
42
4 RESULTADOS
Tratando-se de um caso complexo onde oscilações transversais e axiais ocorrem
simultaneamente, torna-se necessária a validação do modelo proposto para seus casos mais
simples, onde a solução analítica é conhecida. Posteriormente, uma vez comprovada a
aproximação satisfatória dos resultados, sugere-se a aplicação do método para casos de
estruturas complexas.
Desta forma, o estudo dos resultados será apresentado, tanto para o caso estático
quanto o dinâmico, separadamente para os casos de movimento transversal e longitudinal.
4.1
MOVIMENTO TRANSVERSAL
Seja o perfil de viga adotado, como exemplo, W 150 13,5 , onde I  6,87 106 mm4 ,
E  200 GPa , A  1730 mm2 e   7860 kg / m3 . Para o carregamento descrito pela Figura
4.1, sugere-se a seguinte comparação de resultados para a flecha no ponto médio da viga:
Figura 4.1- Exemplo transversal.
43
Tabela 4.1 – Resultados de Análise Estática (Transversal)
Deformação (mm)
Erro
Analítico (Beer, 1995)
-2,03
-
4 Elem.
-2,03
0,0%
20 Elem.
-2,03
0,0%
40 Elem.
-2,03
0,0%
Figura 4.2 - Resultado estático para deformação transversal (20 elementos por viga).
44
Tabela 4.2–Resultados dos Modos Naturais de Vibração (Transversal – 4 Elem.)
4 Elementos
Frequência Natural (Hz)
Modos de
Vibração
Solução Analítica
(Beer, 1995)
Erro
Solução Numérica
1º
242,08
242,20
0,05%
2º
968,31
976,27
0,82%
3º
2178,70
2216,81
1,75%
4º
3873,20
4383,57
13,18%
Tabela 4.3 - Resultados dos Modos Naturais de Vibração (Transversal – 20 Elem.)
20 Elementos
Frequência Natural (Hz)
Modos de
Vibração
Solução Analítica
(Beer, 1995)
Erro
Solução Numérica
1º
242,08
242,08
0,00%
2º
968,31
968,32
0,00%
3º
2178,70
2178,82
0,01%
4º
3873,20
3873,86
0,02%
45
Tabela 4.4 - Resultados dos Modos Naturais de Vibração (Transversal – 40 Elem.)
40 Elementos
Frequência Natural (Hz)
Modos de
Vibração
Solução Analítica
(Beer, 1995)
Erro
Solução Numérica
1º
242,08
242,08
0,00%
2º
968,31
968,31
0,00%
3º
2178,70
2178,70
0,00%
4º
3873,20
3873,27
0,00%
As formas dos modos de vibração natural, para uma simulação com 40 elementos,
podem ser verificadas na ilustração abaixo.
Figura 4.3 - Ondas de deformação para modos de frequências naturais (Transversal).
46
4.2
MOVIMENTO LONGITUDINAL
Seja o perfil de viga genérico adotado, onde I  13333 mm4 , E  200 GPa ,
A  400 mm2 e   7870 kg / m3 . Para o carregamento descrito pela Figura 4.4, sugere-se a
seguinte comparação de resultados para a flecha no ponto médio da viga.
Figura 4.4 - Exemplo longitudinal.
Tabela 4.5 – Resultados de Análise Estática (Longitudinal)
Deformação (mm)
Erro
Analítico (Beer, 1995)
-0,15625
-
4 Elem.
-0,15625
0,0%
20 Elem.
-0,15625
0,0%
40 Elem.
-0,15625
0,0%
47
Figura 4.5 - Resultado estático para deformação longitudinal (20 elementos por viga).
Tabela 4.6 - Resultados dos Modos Naturais de Vibração (Longitudinal – 4 Elem.)
4 Elementos
Frequência Natural (Hz)
Modos de
Vibração
Solução Analítica
(Beer, 1995)
Erro
Solução Numérica
1º
3959,30
3984,78
0,64%
2º
11878,00
12570,54
5,83%
3º
19796,00
22834,79
15,35%
4º
27715,00
33021,12
19,15%
48
Tabela 4.7 - Resultados dos Modos Naturais de Vibração (Longitudinal – 20 Elem.)
20 Elementos
Frequência Natural (Hz)
Modos de
Vibração
Solução Analítica
(Beer, 1995)
Erro
Solução Numérica
1º
3959,30
3960,31
0,03%
2º
11878,00
11905,37
0,23%
3º
19796,00
19923,89
0,65%
4º
27715,00
28065,28
1,26%
Tabela 4.8 - Resultados dos Modos Naturais de Vibração (Longitudinal – 40 Elem.)
40 Elementos
Frequência Natural (Hz)
Modos de
Vibração
Solução Analítica
(Beer, 1995)
Erro
Solução Numérica
1º
3959,30
3959,55
0,01%
2º
11878,00
11884,75
0,06%
3º
19796,00
19828,27
0,16%
4º
27715,00
27802,38
0,32%
As formas dos modos de vibração natural, para uma simulação com 40 elementos,
podem ser verificadas na ilustração abaixo.
49
Figura 4.6 - Ondas de deformação para modos de frequências naturais (Longitudinal).
4.3
ANÁLISE DE ESTRUTURAS COMPLEXAS
Uma vez comprovada a validade do modelo gerado, já que os resultados obtidos foram
satisfatoriamente aproximados para os casos simplificados de vibrações longitudinais e
transversais, é possível expandir a análise para estruturas mais complexas que envolvam uma
combinação de carregamentos e condições de contorno. A seguir consta a aplicação da
metodologia utilizada objetivando a verificação estrutural de um pórtico e, posteriormente, a
expansão para carregamento em uma ponte metálica.
É importante ressaltar que, apesar da elevada complexidade em questão e objetivo
principal do trabalho não envolver cálculos exaustivos para se comprovar a solução analítica
destas estruturas, percebe-se que os resultados demonstrados a seguir se encontram na mesma
ordem de grandeza esperada para os carregamentos desta natureza, estando em conformidade
com as expectativas.
50
4.3.1 Pórtico
Sendo um pórtico composto por vigas W 150 13,5 , onde I  6,87 106 mm4 ,
E  200 GPa , A  1730 mm2 e   7860 kg / m3 , assumindo um carregamento distribuído
horizontal ao longo da viga de suporte esquerda de 2 kN / m , um momento de 15 kN .m na
junção superior esquerda, considerando ligações rígidas entre as vigas e tratando a estrutura
como engastada no solo, temos os seguintes resultados.
Figura 4.7 - Carregamentos e condições de contorno (Pórtico).
Figura 4.8 - Resultado estático para o pórtico (20 elementos por viga).
51
Figura 4.9 - Modos de vibração natural do pórtico (20 elementos por viga).
4.3.2 Ponte Metálica
Considerando a mesma viga do pórtico, desta vez modelando uma estrutura
semelhante à ponte treliçada, assumindo as extremidades como sendo engastadas no solo,
ligações rígidas entre as vigas e atribuindo um carregamento distribuído de 2 kN / m ao longo
de sua base, tem-se os seguintes resultados.
Figura 4.10 - Carregamentos e condições de contorno (Ponte).
52
Figura 4.11 - Resultado estático para a ponte metálica (10 elementos por viga).
Figura 4.12 - Modos de vibração natural da ponte metálica (10 elementos por viga).
53
5 CONCLUSÕES
Por meio da análise e comparação dos resultados obtidos, comprova-se a convergência
entre as soluções numéricas e analíticas. Sendo assim, o programa desenvolvido é, portanto,
válido para diversos casos de estudo, inclusive para graus elevados de complexidade.
Foi observado que a concordância entre os resultados depende diretamente do número
de elementos total do sistema, ainda que, a partir de certo ponto, o aumento do refinamento
não apresente significativos ganhos de aproximação.
Contudo, o aumento desordenado do número de elementos, aliado à ausência de um
método de numeração de nós adequado no programa, resulta em um crescimento vertiginoso e
desnecessário de uso de poder computacional para a obtenção de resultados, fugindo do
objetivo principal proposto de simplificação em análises complexas.
54
RERÊNCIAS BIBLIOGRÁFICAS
ASSIS, Bismarck de O. “Vibrações em Vigas – Um Estudo Analítico e Numérico utilizando o
Método dos Elementos Finitos”. Departamento de Engenharia Mecânica - UFF, Rio de Janeiro,
2011.
Análise Dinâmica, através do Método de Elementos Finitos, de um Compensador Síncrono de 150
MVAr de Fabricação ALSTOM”. Universidade Federal do Para, 2006.
AZEVEDO, Álvaro F.M. “Fundamento dos Elementos Finitos”. Faculdade de Engenharia da
Universidade do Porto. Porto, 2005.
BATHE, Klaus-Jürgen. “Finite Element Procedures”.Ed. Prentice Hall, 2º Edição, New Jersey,
1996.
BEER, Ferdinand P. & Johnston, E. Russell.“Resistência dos Materiais”. Editora Pearson Makron
Books, 3º Edição, São Paulo, 1995.
CLOUGH, Ray W. &Penzien, Joseph.“Dynamics of Strutures”. Computers & Structures, Inc., 2º
Edição, Califórnia, 2005.
CUNHA, Jeferson de S. “Flexão em vigas induzida por difusão de componentes químicos”.
Projeto de Pesquisa PIBIC/CNPq. Departamento de Engenharia Mecânica. Universidade Federal
Fluminense, 2011.
FERREIRA, A.J.M.MATLAB Codes for Finite Element Analysis. Porto, Portugal : Springer,
2008. ISBN 978-1-4020-9199-5
FILHO, Avelino A. “Elementos Finitos –Análise Dinâmica”. Ed. Érica, 2° Edição, São Paulo,
2009.
FISH, Jacob &Belytschko, Ted. “A First Course in Finite Elements”.Ed. Wiley, 1° Edição,
Chichester, 2007. 65
INMAN, Daniel J. “Engineering Vibration”. Ed. Pearson Prentice Hall, 3° Edição, New Jersey,
2007.
MENDES, Paulo & Oliveira, Sérgio. “Análise Dinâmica de Estruturas – Utilização Integrada de
Modelos de Identificação Modal e Modelos de Elementos Finitos”. Laboratório Nacional de
Engenharia Civil, Lisboa, 2008.
MENDONÇA, Paulo de T.R. “Análise Dinâmica pelo Método de Elementos Finitos”.
Departamento de Engenharia Mecânica – UFSC, Santa Catarina, 2006.
OLIVEIRA, José A. V. de A. “Análise Estática e Dinâmica de Estruturas Reticuladas Planas em
Microcomputador”. Dissertação de Mestrado, Universidade do Porto, Porto, 1987.
55
RAMALHO, Amilcar&Antunnes, Fernando J. V. “Vibrações e Ruído”. Departamento de
Engenharia Mecânica – Universidade de Coimbra, Coimbra, 2006.
RAO, Singiresu S. “Vibrações Mecânicas”. Ed. Pearson Prentice Hall, 4º edição, São Paulo, 2008.
SILVA, D. S. “Introdução ao Método dos Elementos Finitos”. UNIOESTE, Foz do
Iguaçu, 2009
6 APÊNDICES
57
6.1
CÓDIGO MATLAB
LIMPAR WORKSPACE ............................................................................................................................. 57
DEFINIÇÃO DE PARAMETROS ................................................................................................................. 57
DISCRETIZAÇÃO .................................................................................................................................... 58
MONTAGEM DO SISTEMA ...................................................................................................................... 59
PLOTAR RESULTADOS ESTÁTICOS ........................................................................................................ 61
PLOTAR RESULTADOS DINÂMICOS ....................................................................................................... 62
ANIMAR MODO DINAMICO .................................................................................................................... 64
6.1.1 Limpar Workspace
close all
clear all
clc
6.1.2 Definição de Parâmetros
E: módulo de elasticidades (GPa) , A: área de seção (mm^2) , I: momento de área (mm^4),
rho: massa específica (kg/m^3)
E=200; A=1730; I=6.87e6; rho=7860;
% Definição juntas (coordenada x e coordenada y, em mm)
juntas = 1e3*[0 0;1 2;2 0;3 2;4 0;5 2;6 0;7 2;8 0];
% Forças e momentos nodais e carregamentos (Magnitude e locais de aplicação)
fx = []; % Em N
xPos = []; % JUNTAS
fy = []; % Em N
yPos = []; % JUNTAS
m = -[]; % Em N-m
mPos = []; % JUNTAS
px = []; % Em N/m, positivo para a direita
pxPos = []; % VIGAS
py = -2.4e3*[1 1 1 1 1]; % em N/m, positivo para cima
pyPos = [3 7 11 15]; % VIGAS
%Restrições (JUNTA a ser restrita)
fim = length(juntas); % Ultima junta
xr = [1 9];
yr = [1 9];
mr = [1 9];
% Elementos por viga
elem_per = 10;
vigas = [1 2;2 3;1 3;2 4;3 4;4 5;3 5;4 6;5 6;6 7;5 7;6 8;7 8;8 9;7 9];
vigas = sort(vigas,2);
tic
58
6.1.3 Discretização
num_juntas = length(juntas);
% Distribuição de nós ao longo das vigas
if isempty(vigas) % Vigas não definidas, assumir cadeia aberta
vigas = [(1:num_juntas-1)' (2:num_juntas)'];
end
% vigas
% vigas = sort(vigas,2)
num_vigas = size(vigas,1);
coord_nos = [];
nn = 0; num_nos = 0;
nos_vigas = zeros(num_vigas,elem_per+1);
for i=1:num_vigas % Discretização das vigas
j = vigas(i,:);
nos_internos = [linspace(juntas(j(1),1),juntas(j(2),1),elem_per+1)'...
linspace(juntas(j(1),2),juntas(j(2),2),elem_per+1)'];
% Verificar se nós das extremidades já existem
[~,ng,ne] = intersect(coord_nos,nos_internos,'rows');
% Incrementar número de nós total com número de nós novos
ninc = size(nos_internos,1)-length(ne);
nn = nn+ninc; % Novo número de nós total
if isempty(ne)
nos_vigas(i,:) = num_nos+1:num_nos+ninc;
coord_nos(nos_vigas(i,1:end),:) = nos_internos(1:end,:);
elseif length(ne) == 1
nos_vigas(i,:) = [ng(1) num_nos+1:num_nos+ninc];
coord_nos(nos_vigas(i,2:end),:) = nos_internos(2:end,:);
else
nos_vigas(i,:) = [ng(1) num_nos+1:num_nos+ninc ng(end)];
coord_nos(nos_vigas(i,2:end-1),:) = nos_internos(2:end-1,:);
end
num_nos = nn;
for k = 1:elem_per
nos_elem(k+(i-1)*elem_per,:) = nos_vigas(i,k:k+1);
end
end
coord_nos = coord_nos*1e-3; % Coordenadas convertidas para metros
num_elem = size(nos_elem,1);
num_nos = size(coord_nos,1);
GDL = 3*num_nos;
deslocamentos = zeros(GDL,1);
forca = zeros(GDL,1);
cx = zeros(num_elem); cy = cx;
rigidez = zeros(GDL);
massa = zeros(GDL);
59
6.1.4 Montagem e Solução do Sistema
%
E
A
I
Conversão de unidades
= E*1e9; % GPa para Pa
= A*1e-6; % mm^2 para m^2
= I*1e-12; % mm^4 para m^4
% Vetor Força (Encontrar nós remapeados e atribuir carga nodal)
[~,xPos,~] = intersect(coord_nos(:,:),1e-3*juntas(xPos,:),'rows');
[~,yPos,~] = intersect(coord_nos(:,:),1e-3*juntas(yPos,:),'rows');
[~,mPos,~] = intersect(coord_nos(:,:),1e-3*juntas(mPos,:),'rows');
forca(xPos*3-2) = fx;
forca(yPos*3-1) = fy;
forca(mPos*3) = m;
% Restrições
[~,xr,~] = intersect(coord_nos(:,:),1e-3*juntas(xr,:),'rows');
[~,yr,~] = intersect(coord_nos(:,:),1e-3*juntas(yr,:),'rows');
[~,mr,~] = intersect(coord_nos(:,:),1e-3*juntas(mr,:),'rows');
restric = [xr'*3-2 yr'*3-1 ...
mr'*3];
% Matriz Rigidez Global
e_GDL = zeros(size(nos_elem,1),6);
for e=1:num_elem;
e_GDL(e,:) = [nos_elem(e,1)*3-2:nos_elem(e,1)*3 ...
nos_elem(e,2)*3-2:nos_elem(e,2)*3]; % GDL nós do elemento
% xa e ya: projeções do comprimento do elemento
% ll: comprimento do elemento
dx(e) = diff(coord_nos(nos_elem(e,1:2),1)); %#ok<*SAGROW>
dy(e) = diff(coord_nos(nos_elem(e,1:2),2));
comp(e) = sqrt(dx(e)^2 + dy(e)^2);
l = dx(e)/comp(e); % cos
m = dy(e)/comp(e); % sen
% Matriz de transformada
L= [ l m 0 0 0 0 ;...
-m l 0 0 0 0 ;...
0 0 1 0 0 0 ;...
0 0 0 l m 0 ;...
0 0 0 -m l 0 ;...
0 0 0 0 0 1 ];
% Submatriz rigidez
ll = comp(e);
k1 = [...
E*A/ll
0
0
12*E*I/ll^3
0
6*E*I/ll^2
-E*A/ll
0
0 -12*E*I/ll^3
0
6*E*I/ll^2
0 -E*A/ll
0
6*E*I/ll^2
0 -12*E*I/ll^3
4*E*I/ll
0 -6*E*I/ll^2
0 E*A/ll
0
-6*E*I/ll^2
0 12*E*I/ll^3
2*E*I/ll
0 -6*E*I/ll^2
% Incrementar matriz rigidez global
0
6*E*I/ll^2
2*E*I/ll
0
-6*E*I/ll^2
4*E*I/ll
;...
;...
;...
;...
;...
];
60
rigidez(e_GDL(e,:),e_GDL(e,:))=...
rigidez(e_GDL(e,:),e_GDL(e,:))+L'*k1*L;
% Submatriz Massa
m1 = (rho*A*ll/420)*[...
140
0
0
0
156
22*ll
0
22*ll
4*ll^2
70
0
0
0
54
13*ll
0 -13*ll -3*ll^2
70
0
0
140
0
0
0
54
13*ll
0
156
-22*ll
0
-13*ll
-3*ll^2
0
-22*ll
4*ll^2
;...
;...
;...
;...
;...
];
% Incrementar matriz massa global
massa(e_GDL(e,:),e_GDL(e,:))=...
massa(e_GDL(e,:),e_GDL(e,:))+L'*m1*L;
% Incrementar vetor força com carga distribuida
v = ceil( e/elem_per ); % Encontrar viga atual
x_ind = find(pxPos == v); % Indice da carga distribuida em x
y_ind = find(pyPos == v); % Indice da carga distribuida em y
cx = px(x_ind)*ll; cy = py(y_ind)*ll; % Contribuições (N)
% Prevenir posições vazias no vetor carga (atribuir zero)
if isempty(x_ind)
cx = 0;
end
if isempty(y_ind)
cy = 0;
end
f1 = [cx/2, cy/2, (cy-cx)*ll/12, ...
cx/2, cy/2, (cx-cy)*ll/12]';
forca(e_GDL(e,:)) = forca(e_GDL(e,:))+...
f1;
end
toc
% Solução Estática
GDL_ativo=sort(setdiff([1:GDL]',restric)); %#ok<NBRAK>
deslocamentos(GDL_ativo) = rigidez(GDL_ativo,GDL_ativo)\forca(GDL_ativo);
% Solução Dinâmica
[Autovetores, Autovalores] = eig(rigidez(GDL_ativo,GDL_ativo),...
massa(GDL_ativo,GDL_ativo));
v = diag(sqrt(Autovalores));
[v,ind] = sort(v,'ascend');
Autovetores = Autovetores(:,ind);
format shortg
xy = setdiff(1:GDL_ativo,3:3:GDL_ativo);
% Separar autovalores e autovetores
var = abs(Autovetores);
61
[~,nx] = find(var(1:3:length(GDL_ativo),:) > 1e-3);
nx = unique(sort(nx,'ascend')); % Posição dos modos em x
[~,ny] = find(var(2:3:length(GDL_ativo),:) > 1e-3);
ny = unique(sort(ny,'ascend')); % Posição dos modos em y
[~,notaveis] = sort(max(var(xy,:),[],1),'ascend');
Elapsed time is 0.145718 seconds.
6.1.5 Plotar Resultados Estáticos
colormap(jet(256))
% Definir escala automática
U = deslocamentos;
Uabs = (U(1:3:end-2).^2+U(2:3:end-1).^2).^(1/2);
if sum(Uabs) ~= 0
escala = 0.1*max([max(coord_nos(:,1))-min(coord_nos(:,1)),...
max(coord_nos(:,2)-min(coord_nos(:,2)))])/...
max(Uabs);
else
escala = 0;
end
deform = coord_nos+escala*[U(1:3:end-2) U(2:3:end-1)];
hold on
for i = 1:num_vigas
plot(coord_nos(nos_vigas(i,:),1),...
coord_nos(nos_vigas(i,:),2),'.-','Color',[.8 .8 .8])
surface([deform(nos_vigas(i,:),1)';deform(nos_vigas(i,:),1)'],...
[deform(nos_vigas(i,:),2)';deform(nos_vigas(i,:),2)'],...
[Uabs(nos_vigas(i,:))';Uabs(nos_vigas(i,:))'],...
'facecol','no',...
'edgecol','interp',...
'linew',1,'LineSmoothing','On');
colorbar('EastOutside')
end
title(['Estrutura deformada (escala de deformação ',...
num2str(escala,'%.2f'),')'])
xlabel('x (m)')
ylabel('y (m)')
% Redimensionar figura
xlim = [min(coord_nos(:,1)) max(coord_nos(:,1))];
dx = xlim(2)-xlim(1);
ylim = [min(coord_nos(:,2)) max(coord_nos(:,2))];
dy = ylim(2)-ylim(1);
if dx >= dy && dx ~= 0
set(gca,'XLim',[xlim(1)-.1*dx xlim(2)+.1*dx],...
'YLim',[(ylim(1)+ylim(2))/2-.6*dx (ylim(1)+ylim(2))/2+.6*dx])
elseif dy ~= 0
set(gca,...
'XLim',[(xlim(1)+xlim(2))/2-.6*dy (xlim(1)+xlim(2))/2-.6*dy],...
62
'YLim',[ylim(1)-.1*dy ylim(2)+.1*dy])
end
6.1.6 Plotar Resultados Dinâmicos
m = 1:4;
figure(2)
din = zeros(4,num_vigas);
for i= 1:4;
subplot(2,2,i)
figure(2)
colormap(jet(256))
hold on
U = zeros(GDL,1);
U(GDL_ativo,1) = Autovetores(:,m(i));
Uabs = (U(1:3:end-2).^2+U(2:3:end-1).^2).^(1/2);
escala = 0.07*max([max(coord_nos(:,1))-min(coord_nos(:,1)),...
max(coord_nos(:,2)-min(coord_nos(:,2)))])/...
max(Uabs);
deform = coord_nos+escala*[U(1:3:end-2) U(2:3:end-1)];
caxis manual
caxis([0,max(Uabs)])
for j = 1:num_vigas
plot(coord_nos(nos_vigas(j,:),1),...
coord_nos(nos_vigas(j,:),2),'.-','Color',[.8 .8 .8]);
din(i,j) = surface(...
63
[deform(nos_vigas(j,:),1)';deform(nos_vigas(j,:),1)'],...
[deform(nos_vigas(j,:),2)';deform(nos_vigas(j,:),2)'],...
[Uabs(nos_vigas(j,:))';Uabs(nos_vigas(j,:))'],...
'facecol','no',...
'edgecol','interp',...
'linew',1,'LineSmoothing','On');
end
title([num2str(i),'º Modo (',num2str(v(m(i)),'%.2f'),'Hz)'])
xlabel('x (m)')
ylabel('y (m)')
% Redimensionar figura
xlim = [min(coord_nos(:,1)) max(coord_nos(:,1))];
dx = xlim(2)-xlim(1);
ylim = [min(coord_nos(:,2)) max(coord_nos(:,2))];
dy = ylim(2)-ylim(1);
if dx >= dy && dx ~= 0
set(gca,'XLim',[xlim(1)-.1*dx xlim(2)+.1*dx],...
'YLim',[(ylim(1)+ylim(2))/2-.6*dx (ylim(1)+ylim(2))/2+.6*dx])
elseif dy ~= 0
set(gca,'XLim',[xlim(1)+(dx-1.2*dy)/2 xlim(1)+(dx+1.2*dy)/2],...
'YLim',[ylim(1)-.1*dy ylim(2)+.1*dy])
end
end
toc
Elapsed time is 3.866980 seconds.
64
6.1.7 Animar Modo Dinâmico
figure(2)
arco = 0:0.1:floor(10*4.5*pi)/10;
for k = arco
for i= 1:4;
U = zeros(GDL,1);
U(GDL_ativo,1) = Autovetores(:,m(i));
Uabs = (U(1:3:end-2).^2+U(2:3:end-1).^2).^(1/2);
escala = 0.07*max([max(coord_nos(:,1))-min(coord_nos(:,1)),...
max(coord_nos(:,2)-min(coord_nos(:,2)))])/...
max(Uabs);
deform = coord_nos+sin(k)*escala*[U(1:3:end-2) U(2:3:end-1)];
for j = 1:num_vigas
% Remapear cores
set(din(i,j),...
'XData',[deform(nos_vigas(j,:),1)';deform(nos_vigas(j,:),1)'],...
'YData',[deform(nos_vigas(j,:),2)';deform(nos_vigas(j,:),2)'],...
'ZData',[Uabs(nos_vigas(j,:))';Uabs(nos_vigas(j,:))'],...
'CData',abs(sin(k)*[Uabs(nos_vigas(j,:))';Uabs(nos_vigas(j,:))']))
end
end
pause(.05)
end
for i=1:4
figure(2)
subplot(2,2,i)
colorbar('EastOutside')
end
Published with MATLAB® R2012b