UNIVERSIDADE FEDERAL FLUMINENSE
TCE - Escola de Engenharia
TEM - Departamento de Engenharia Mecânica
PROJETO DE GRADUAÇÃO II
Título do Projeto:
DESENVOLVIMENTO DE UM ALGORITMO PARA ANÁLISE ESTÁTICA E
DINÂMICA DE TRELIÇAS PLANAS
Autor :
CAMILA MACHADO DE ARAUJO
Orientador :
ANGELA CRISTINA CARDOSO DE SOUZA
Data: 15 de julho de 2015
CAMILA MACHADO DE ARAUJO
DESENVOLVIMENTO DE UM ALGORITMO PARA ANÁLISE ESTÁTICA E
DINÂMICA DE TRELIÇAS PLANAS
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. Angela Cristina Cardoso de Souza
Niterói
2015
Ficha Catalográfica elaborada pela Biblioteca da Escola de Engenharia e Instituto de Computação da UFF
DEDICAT
Dedico à minha família:
Walmir Gomes de Araujo, Virginia Machado de Araujo e Carolina Machado de Araujo.
AGRADECIMENTOS
Agradeço primeiramente à minha família, meu pai Walmir Gomes de Araujo, minha
mãe Virginia Machado de Araujo e minha irmã Carolina Machado de Araujo, por não
medirem esforços para me ajudar e fornecerem todo e qualquer suporte para que fosse
possível a conclusão da graduação em Engenharia Mecânica.
À todos os professores do curso, que foram de extrema importância na minha vida
acadêmica, em especial à professora Angela Cristina, por me auxiliar no desenvolvimento
deste trabalho.
Aos colegas de curso, que compartilharam das mesmas aflições e ansiedades ao longo
dessa trajetória e dividiram seus conhecimentos para o crescimento conjunto de todos, em
especial aos amigos Gabriel Cunha e Pedro Henrique Alves, pela especial colaboração no
desenvolvimento deste trabalho.
À Equipe Tuffão Baja SAE, por me apresentar à Engenharia Mecânica e fazer com
que eu descobrisse o gosto pela área e assim tomasse um novo rumo na minha carreira
profissional.
À Noêmia Maia e todas as que já passaram por seu lar, em especial às minhas amigas
Ludmila Cabral, Jéssica Helena, Jéssica D’ávila e Joana Sadio, por tornarem a vida em
Niterói mais fácil, saudável e divertida.
Por fim, agradeço todos aqueles que indiretamente contribuíram para a realização da
minha graduação até ser possível a execução desta monografia.
RESUMO
Treliças são estruturas comumente utilizadas na construção de cobertura de galpões,
guindastes, torres e pontes com o objetivo de obter uma estrutura leve e rígida. Devido às
forças a qual está sujeita e a que ela se destina, qualquer falha que comprometa sua
integridade pode causar grandes danos. O objetivo desse trabalho é apresentar o
desenvolvimento do cálculo estrutural, seja ele estático ou dinâmico, para treliças
bidimensionais através da modelagem pelo Método de Elementos Finitos (MEF). Ao final do
desenvolvimento teórico, será apresentado o desenvolvimento do algoritmo feito no software
MATLAB®
a fim de otimizar o cálculo estrutural desse tipo de estrutura.
Palavras-Chave: Treliças; Algoritmo; Elementos Finitos, Bidimensional.
ABSTRACT
Trusses are structures commonly used in warehouses cover construction, cranes,
towers and bridges in order to obtain a light and rigid structure. Due to the forces which are
subject and that it is intended, any failure to compromise their integrity can cause major
damage. The main of this paper is to present the development of structural design, whether
static or dynamic, for two-dimensional lattices through modeling by Finite Element Method
(FEM). At the end of the theoretical development, the development of the algorithm done in
MATLAB®
software to optimize the structural design of this type of structure will be
presented.
Key-Words: Trusses; Algorithm; Finite Elements, Two-dimensional.
LISTA DE ILUSTRAÇÕES
Figura 1 – Ponte treliçada [11] ............................................................................................... 13
Figura 2 - Cobertura [12] ........................................................................................................ 13
Figura 3 – Pórtico [13] ............................................................................................................ 14
Figura 4. - (a) Elemento de uma treliça sob carregamento trativo; (b) Elemento de uma
treliça sob carregamento compressivo (Fonte: BEER - 2012) ................................................ 15
Figura 5 - (a) Forças atuantes no elemento (b) Infinitésimo do elemento (Fonte: Imagens do
autor). ....................................................................................................................................... 17
Figura 6 - (a) Corpo discretizado com elementos tetragonais (Fonte: Imagem do autor); (b)
Geometrias de elementos utilizadas [14] ................................................................................. 25
Figura 7 - Gráfico das funções de forma (Fonte: Imagem do autor). ..................................... 27
Figura 8 - Elemento com relação ao referencial local e global (Fonte: Imagem do autor). .. 31
Figura 9 - Exemplo de estrutura de treliças (Fonte: Imagem do autor) .................................. 33
Figura 10 - Estrutura discretizada com 3 elementos e 6 graus de liberdade (Fonte: Imagem
do autor). .................................................................................................................................. 33
Figura 11 - Primeiro elemento (Fonte: Imagem do autor). ..................................................... 34
Figura 12 - Segundo elemento (Fonte: Imagem do autor). ...................................................... 35
Figura 13 - Terceiro elemento (Fonte: Imagem do autor). ...................................................... 35
Figura 14 - Reações em apoio para uma estrutura bidimensional (Fonte: BEER – 2012) ..... 37
Figura 15 – (a) Estrutura do exemplo 1 (b) Estrutura discretizada (Fonte: Imagem do autor).
.................................................................................................................................................. 39
Figura 16 – (a) Estrutura do exemplo 2; (b) Estrutura discretizada (Fonte: Imagem do
autor). ....................................................................................................................................... 42
LISTA DE SÍMBOLOS
𝐴(𝑥) − Área da seção transversal
𝐶1 − Amplitude
𝐸 − Módulo de elasticidade
�⃗� − Vetor de carregamento
𝑓(𝑥) − Carregamento distribuido por unidade de comprimento
Ke - Matriz de rigidez do elemento
𝐾𝐺 − Matriz de rigidez global
𝑀𝑒 − Matriz de massado elemento
𝑀𝐺 − Matriz de massa global
𝑁 − Carregamento axial
𝑢(𝑥) − Deslocamento longitudinal
�⃗⃗� − Vetor deslocamento
�̈⃗⃗� − Vetor aceleração
𝑣(𝑥) − Função arbitrária
𝑣𝑖⃗⃗⃗ ⃗ − Autovetores
𝛽 − Matriz de transformação
𝜃 − Ângulo entre os eixos x
𝜆𝑖 − Autovalores
𝜌 − Densidade do material
𝜎 − Tensão normal
𝜙 − Ângulo de fase
𝛹 − Funções de forma
𝜔 − Frequência natural
SUMÁRIO
1 INTRODUÇÃO .......................................................................................................................... 13 1.1. CONSIDERAÇÕES GERAIS ...................................................................................................... 13 1.2. OBJETIVO ................................................................................................................................... 15 1.3. EQUAÇÃO DO EQUILÍBRIO - FORMULAÇÃO DIFERENCIAL ............................................ 16 1.4. EQUAÇÃO DO EQUILÍBRIO - FORMULAÇÃO INTEGRAL .................................................. 18 1.5. PROBLEMA DE AUTOVALOR ................................................................................................. 19 2 MÉTODO DE ELEMENTOS FINITOS .................................................................................. 23 2.1. HISTÓRIA ................................................................................................................................... 23 2.2. FUNDAMENTOS DO MÉTODO DE ELEMENTOS FINITOS .................................................. 24 2.3. DISCRETIZAÇÃO DO DOMÍNIO .............................................................................................. 26 2.4. MONTAGEM DO SISTEMA GLOBAL ...................................................................................... 30 2.4.1 MATRIZES DO ELEMENTO ......................................................................................................... 30 2.4.2 MATRIZES GLOBAIS .................................................................................................................. 32 2.5. CONDIÇÕES DE CONTORNO .................................................................................................. 36 3 RESULTADOS ........................................................................................................................... 39 3.1. EXEMPLO 1 ................................................................................................................................ 39
3.1.1 ESTÁTICO: [𝐾]�⃗⃗� = �⃗� ............................................................................................................. 40
3.1.2 DINÂMICO: ([𝐾] − 𝜔²[𝑀])�⃗� = 0⃗⃗ ......................................................................................... 41 3.2. EXEMPLO 2 ................................................................................................................................ 42
3.2.1 ESTÁTICO: [𝐾]�⃗⃗� = �⃗� ............................................................................................................. 43
3.2.2 DINÂMICO: ([𝐾] − 𝜔²[𝑀])�⃗� = 0⃗⃗ ......................................................................................... 46 4 CONCLUSÕES .......................................................................................................................... 49 5 REFERÊNCIAS BIBLIOGRÁFICAS ...................................................................................... 51 6 APÊNDICES ............................................................................................................................... 52 6.1. DADOS DE ENTRADA .............................................................................................................. 52 6.1.1 EXEMPLO 1: .............................................................................................................................. 52 6.1.2 EXEMPLO 2: .............................................................................................................................. 52 6.2. CÓDIGO FONTE ......................................................................................................................... 54
13
1 INTRODUÇÃO
1.1. CONSIDERAÇÕES GERAIS
Treliças são estruturas de elementos delgados ligados por suas extremidades,
geralmente vigas em “I” ou em “U”, cantoneiras e barras, que são unidas por meio de
soldagem, conexões rebitadas, parafusos, ou pinos, sem qualquer consideração de atrito ou
outras forças que impeçam a livre rotação dos elementos em relação a estas juntas conectoras,
de forma que não possam transmitir movimento. São comumente encontradas em construções
ou máquinas, como pode ser visto nas figuras 1, 2 e 3 abaixo, com a função de formar uma
estrutura rígida, atuando como uma distribuidora de forças.
Figura 1 – Ponte treliçada [11]
Figura 2 - Cobertura [12]
14
Figura 3 – Pórtico [13]
Nas estruturas compostas por treliças, as forças externas são coplanares aplicadas
somente aos nós, não existindo qualquer transmissão de momento entre seus elementos, sendo
assim, os elementos de barras estarão sujeitos apenas a esforços normais. Ou seja, cada uma
das barras da treliça vai estar submetida a somente duas forças axiais internas e colineares ao
eixo da barra, uma em cada extremidade que atuam em uma mesma linha de ação com o
mesmo módulo, porém, em sentidos opostos, podendo constituir assim situações em que os
elementos da estrutura estarão sujeitos a forças de tração ou compressão, como ilustrado na
figura 4. Como os elementos são considerados delgados, os efeitos de torção e
cisalhamento tornam-se desprezíveis. Para fins de cálculo, a seção transversal do elemento
será considerada como arbitrária e as únicas características geométricas empregadas na
modelagem física do problema serão as áreas e comprimentos dos elementos de barra, além
das orientações espaciais dos nós através de um sistema de coordenadas cartesianas x-y.
15
(a) (b)
Figura 4. - (a) Elemento de uma treliça sob carregamento trativo; (b) Elemento de uma
treliça sob carregamento compressivo (Fonte: BEER - 2012)
Para o caso de estudo, os materiais obedecem a Lei de Hooke, que estabelece que a
tensão normal interna é uma função linear da deformação, sendo o coeficiente angular dessa
função denominado módulo de elasticidade, ou módulo de Young e varia de acordo com o
material empregado nos elementos [2].
Ao optar pelo cálculo estrutural utilizando uma resolução numérica é possível
observar os diversos problemas com importância para a Engenharia que podem ser descritos
em termos de equações com derivadas parciais. O Método dos Elementos Finitos (MEF) é,
atualmente, o método numérico mais utilizado para obter soluções aproximadas para análise
de problemas estruturais. Dado o caráter aproximado das soluções fornecidas por este método,
o desconhecimento dos seus fundamentos pode conduzir a resultados inexatos na sua
aplicação.
1.2. OBJETIVO
O objetivo geral do trabalho é estudar os sistemas de treliças, de forma estática e
dinâmica, através de uma abordagem numérica que possa ser aplicada para qualquer
configuração de treliças bidimensionais, com qualquer sistema de forças atuantes.
Para tal, será utilizada uma formulação matricial através do MEF, no qual são
introduzidas matrizes de rigidez e massa locais para cada elemento, calculadas a partir de sua
localização espacial (determinada pelas coordenadas dos nós em suas extremidades), por
características geométricas dos elementos, como as áreas das seções transversais,
comprimento, e também por componentes físicos, tal como o módulo de Young e densidade
do material. Em seguida, é iniciado um processo de montagem a partir das matrizes de rigidez
16
e massa de cada elemento, transformando as coordenadas locais de cada barra em um sistema
global. Assim será obtida a matriz de rigidez e massa globais, da qual serão extraídas
importantes informações do sistema. Os passos para obtenção dessas matrizes serão
abordados no Capítulo 2, onde será visto com detalhes o Método dos Elementos Finitos.
Finalmente, através das condições de contorno (forças externas aplicadas sobre os nós
e a posição espacial dos mesmos, além das condições de fixação da estrutura) informação
como os deslocamentos de cada nó e modos de vibração podem ser obtidos.
Para os cálculos matemáticos de grandes estruturas planas, um código foi
desenvolvido a partir do software matemático MATLAB, que através de informações
fornecidas pelo usuário sobre a estrutura desejada, realiza os cálculos necessários para
obtenção das deformações e modos de vibração. Assim sendo, é feita uma avaliação sobre os
valores encontrados, analisando se satisfazem ou não os requisitos de segurança necessários à
integridade estrutural a qual será submetida esse sistema de treliças.
Nas demonstrações matemáticas, um caso simplificado utilizando somente uma barra
carregada axialmente será analisado. Essa barra também pode ser chamada de elemento.
1.3. EQUAÇÃO DO EQUILÍBRIO - FORMULAÇÃO DIFERENCIAL
É considerada uma barra de comprimento 𝐿, carregada axialmente em suas
extremidades por uma força 𝑁, além de contar com um carregamento 𝑓(𝑥) distribuído em seu
comprimento, como ilustrado na figura 5 (a). Por consequência desses carregamentos, a
mesma sofre um deslocamento longitudinal dado por 𝑢(𝑥, 𝑡).
(a)
17
(b)
Figura 5 - (a) Forças atuantes no elemento (b) Infinitésimo do elemento (Fonte: Imagens
do autor).
O equilíbrio de forças é obtido de um infinitésimo do elemento no sentido
longitudinal ao comprimento, como pode ser observado na figura 5 (b).
𝑁 + 𝑑𝑁 − 𝑁 + 𝑓𝑑𝑥 = 𝜌𝐴(𝑥)𝑑𝑥𝜕2𝑢
𝜕𝑡2 (1.1)
Sendo o carregamento axial dado por 𝑁 e que o material obedece a Lei de Hooke [2]:
𝑁 = 𝜎𝐴(𝑥) = 𝐸𝐴(𝑥)𝜕𝑢
𝜕𝑥 (1.2)
Onde:
ρ é a densidade do material da barra.
𝐸 é o módulo de elasticidade do material da barra.
𝐴(𝑥) é a área da seção transversal da barra em determinado ponto da mesma.
E substituindo a Eq. (1.2) na Eq. (1.1), temos a equação do equilíbrio:
18
𝜕
𝜕𝑥(𝐸𝐴(𝑥)
𝜕𝑢
𝜕𝑥) + 𝑓(𝑥) = 𝜌𝐴(𝑥)
𝜕2𝑢
𝜕𝑡2 (1.3)
1.4. EQUAÇÃO DO EQUILÍBRIO - FORMULAÇÃO INTEGRAL
Uma vez que se deseja obter a solução do corpo como um todo, é necessário
considerar uma formulação generalizada. Sendo assim, para a resolução por elementos finitos,
a Eq. (1.3) deve ser transformada na forma integral. Para tanto, a equação é multiplicada por
uma função arbitrária 𝑣(𝑥) e integrada sobre o domínio da barra Ω [0; 𝐿] [3], como mostrado
a seguir:
∫𝜕
𝜕𝑥(𝐸𝐴
𝜕𝑢
𝜕𝑥) 𝑣(𝑥) 𝑑𝑥
𝐿
0
+ ∫ 𝑓(𝑥)𝑣(𝑥) 𝑑𝑥𝐿
0
= ∫ 𝜌𝐴𝜕2𝑢
𝜕𝑡2
𝐿
0
𝑣(𝑥) 𝑑𝑥 (1.4)
Para reduzir a ordem da equação diferencial acima dos termos em função de 𝑥, o
primeiro termo é integrado por partes, onde é considerado:
ℎ = 𝑣(𝑥) 𝐷𝑒𝑟𝑖𝑣𝑎𝑛𝑑𝑜⇒
𝑑ℎ
𝑑𝑥=𝑑𝑣
𝑑𝑥 ∴ 𝑑ℎ =
𝑑𝑣
𝑑𝑥𝑑𝑥 (1.5)
𝑑𝑤 =𝜕
𝜕𝑥(𝐸𝐴
𝜕𝑢
𝜕𝑥) 𝑑𝑥
𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑛𝑑𝑜⇒ 𝑤 = 𝐸𝐴
𝜕𝑢
𝜕𝑥 (1.6)
Substituindo as Eq. (1.5) e (1.6) na formulação por partes abaixo,
∫ ℎ 𝑑𝑤𝐿
0
= ℎ𝑤]0𝐿 −∫ 𝑤
𝐿
0
𝑑ℎ (1.7)
3º termo 2º termo 1º termo
Forças de corpo Forças de inércia
Rigidez transversal da
barra
19
Tem-se que:
∫𝜕
𝜕𝑥(𝐸𝐴
𝜕𝑢
𝜕𝑥)𝑣(𝑥) 𝑑𝑥
𝐿
0
= 𝐸𝐴𝜕𝑢
𝜕𝑥𝑣(𝑥)]
0
𝐿
−∫ 𝐸𝐴𝜕𝑢
𝜕𝑥
𝑑𝑣
𝑑𝑥 𝑑𝑥
𝐿
0
(1.8)
Substituindo a Eq. (1.8) na Eq. (1.4), tem-se:
𝐸𝐴𝜕𝑢
𝜕𝑥𝑣(𝑥)]
0
𝐿
−∫ 𝐸𝐴𝜕𝑢
𝜕𝑥
𝑑𝑣
𝑑𝑥 𝑑𝑥
𝐿
0
+∫ 𝑓(𝑥)𝑣(𝑥) 𝑑𝑥𝐿
0
= ∫ 𝜌𝐴𝜕2𝑢
𝜕𝑡2
𝐿
0
𝑣(𝑥) 𝑑𝑥 (1.9)
Manipulando os termos da Eq. (1.9) e substituindo a variável 𝑁 pelos termos
indicados, tem-se então a solução geral na forma integral apresentada na Eq. (1.10) a seguir.
∫ 𝜌𝐴𝜕2𝑢
𝜕𝑡2𝑣(𝑥) 𝑑𝑥
𝐿
0
+∫ 𝐸𝐴𝜕𝑢
𝜕𝑥
𝐿
0
𝑑𝑣
𝑑𝑥 𝑑𝑥 = ∫ 𝑓(𝑥)𝑣(𝑥) 𝑑𝑥
𝐿
0
+𝑁 𝑣(𝑥)]0𝐿 (1.10)
Onde a primeira integral representa o trabalho das forças de inércia, a segunda integral
caracteriza o trabalho interno de deformação, a terceira integral indica o trabalho das forças de
corpo e o quarto termo, o trabalho das forças de fronteira.
1.5. PROBLEMA DE AUTOVALOR
Para efeito de utilização no método numérico, a Eq. (1.10) deve compor, juntamente
com as soluções para os demais elementos do corpo (no caso, as barras) e ser representada
através da forma matricial para um sistema forçado não amortecido com vários graus de
liberdades [6], como apresentado abaixo:
𝑁
20
[𝑀]�̈⃗⃗� + [𝐾]�⃗⃗� = �⃗� (1.11)
Onde:
[𝑀] é a matriz de massa;
[𝐾] é a matriz de rigidez;
�⃗⃗� é o vetor de deslocamento;
�̈⃗⃗� é o vetor de aceleração;
�⃗� é o vetor de carregamento.
Para a resolução da equação acima, o problema é considerado como um estudo de
vibração livre, isto é �⃗� = 0, tornando a equação diferencial homogênea e fazendo com que a
seja reduzida para:
[𝑀]�̈⃗⃗� + [𝐾]�⃗⃗� = 0 (1.12)
Admitindo uma solução da forma:
𝑢(𝑡) = 𝑈𝑖𝑇(𝑡) 𝑖 = 1, 2, 3, … , 𝑛 (1.13)
Onde 𝑈𝑖 é uma constante e 𝑇 é uma função que depende apenas do tempo. Se
comparados dois valores distintos da equação anterior, nota-se que a razão entre as amplitudes
das duas coordenadas é constante no tempo. Ou seja, todas as coordenadas possuem sincronia
em seus movimentos, onde a configuração não muda de forma, porém sua amplitude muda.
𝑢𝑖(𝑡)
𝑢𝑗(𝑡)=𝑈𝑖𝑈𝑗= 𝑐𝑡𝑒 (1.14)
Essa configuração é dada pelo vetor �⃗⃗⃗� abaixo:
21
�⃗⃗⃗� = [
𝑈1𝑈2⋮𝑈𝑛
] (1.15)
Substituindo a Eq. (1.13) na (1.12), tem-se que:
[𝐾]�⃗⃗⃗�𝑇(𝑡) + [𝑀]�⃗⃗⃗��̈�(𝑡) = 0⃗⃗ (1.16)
Que pode ser escrita na forma escalar, tal que:
(∑𝑘𝑖𝑗
𝑛
𝑗=1
𝑈𝑗)𝑇(𝑡) + (∑𝑚𝑖𝑗
𝑛
𝑗=1
𝑈𝑗) �̈�(𝑡) = 0 𝑖 = 1, 2, 3, … , 𝑛 (1.17)
Da Eq. (1.17) se obtém a relação abaixo:
−�̈�(𝑡)
𝑇(𝑡)=(∑ 𝑘𝑖𝑗𝑈𝑗
𝑛𝑗=1 )
(∑ 𝑚𝑖𝑗𝑈𝑗𝑛𝑗=1 )
𝑖 = 1, 2, 3, … , 𝑛 (1.18)
Concluindo que o lado esquerdo é independente do índice 𝑖 e o lado direito
independente do tempo, os dois lados são iguais a uma constante e considerando essa
constante como 𝜔2, a Eq. (1.18) pode ser escrita como:
�̈�(𝑡) + 𝜔2𝑇(𝑡) = 0 (1.19)
Ou na forma matricial:
([𝐾] − 𝜔2[𝑀])�⃗� = 0⃗⃗ (1.20)
Cuja solução é dada por:
𝑇(𝑡) = 𝐶1 cos(𝜔𝑡 + 𝜙) (1.21)
22
Onde 𝐶1 é a amplitude e 𝜙 o ângulo de fase, ambas constantes. A Eq. (1.21) mostra
que para uma mesma frequência e mesmo ângulo de fase, seja qual for a coordenada, pode
executar um movimento harmônico. Porém, a frequência 𝜔 deve satisfazer a Eq. (1.21) não
podendo assumir qualquer valor. Como a Eq. (1.20) representa um sistema de 𝑛 equações
lineares e homogêneas em �⃗⃗⃗�, uma solução trivial é �⃗⃗⃗� = 0⃗⃗. Já uma solução não trivial é dada
por:
det |[𝐾] − 𝜔2[𝑀]| = 0 (1.22)
Que pode ser reescrita por:
det |[𝐾] − 𝜆[𝑀]| = 0 (1.23)
A Eq. (1.23) é denominada equação característica, 𝜆 e 𝜔 chamados, respetivamente,
de autovalor e frequência natural do sistema (dado em radianos por segundo).
Das raízes da Eq. (1.23) resultam 𝑛 valores para 𝜆. Caso 𝜆1, 𝜆2, … , 𝜆1𝑛 representem
os autovalores do sistema em ordem de magnitude, 𝜔1, 𝜔2, … , 𝜔𝑛 representam as frequências
naturais do sistema e seu menor valor 𝜔1, chamado de frequência fundamental ou primeira
frequência natural. Para cada autovalor 𝜆𝑖, existe um autovetor v𝑖⃗⃗⃗ ⃗ associado, próprio de um
modo de vibração, que é determinado a partir da relação:
|[𝐾] − 𝜆𝑖[𝑀]|v𝑖⃗⃗⃗ ⃗ = 0 (1.24)
23
2 MÉTODO DE ELEMENTOS FINITOS
2.1. HISTÓRIA
Derivado da análise matricial para modelos reticulados, criado principalmente com a
finalidade de utilização na indústria aeronáutica britânica, o Método dos Elementos Finitos
(MEF) se desenvolveu juntamente com a computação digital devido à necessidade de se
trabalhar com estruturas de modelos contínuos. Foi desenvolvido principalmente por
Engenheiros Aeronáuticos para a análise de tensões na superfície das asas dos aviões,
inicialmente por Argyris e Kelsey (1955) e Turner, Clough, Martin e Topp (1956) [8]. Sendo
assim, pode ser observada a relação de dependência e evolução do MEF de acordo com o
desenvolvimento da computação digital.
Na análise de problemas tridimensionais pelo MEF considerando o efeito de
temperaturas em sólidos de formas complexas, Gallagher, Padlog e Bijlaard (1962) foram os
primeiros a dar início ao estudo [8].
A abordagem era intuitiva buscando substituir o modelo contínuo em um conjunto de
elementos de dimensão limitada ligados por pontos denominados “nós”. A formulação partia
do princípio dos deslocamentos virtuais (formulação direta) e não se tinha garantia de
convergência para a solução do problema.
A formulação a partir da minimização da grandeza escalar (funcional de energia
potencial – é igual a energia de deformação elástica menos duas vezes o trabalho das forças
externas) foi apresentado por Melosh (1963). Já a formulação baseada em outros funcionais
da mecânica dos sólidos deformáveis foi apresentada por Veubeke (1965). Viu-se então que o
MEF já havia sido apresentado por Lord Rayleigh (1870), Walther Ritz (1909) e Richard
Courant (1943). Sendo assim um caso particular do método de Ritz onde eram estabelecidos
critérios de convergência podendo ser empregado em qualquer meio contínuo regido por um
funcional (formulação variacional) [8].
Na formulação variacional existe uma condição de “estacionariedade” que possibilita a
existência de equações diferenciais e condições de contorno do meio contínuo. Sendo assim
fica simples se chegar a uma solução aproximada desse meio, pois chegar a essa solução
aproximada para a condição de estacionariedade é o mesmo que buscar uma aproximação das
24
correspondentes equações diferenciais. No método de Ritz, para o domínio de meio contínuo
em questão, designa-se uma solução aproximada para o campo de deslocamento em função de
parâmetros generalizados a serem determinados. No MEF, para cada subdomínio (elemento
finito) é arbitrada uma solução aproximada em função dos deslocamentos dos seus nós onde
os mesmos são determinados através das suas condições de contorno.
A formulação variacional (ou formulação fraca) trata da resolução por elementos
finitos de uma ampla variedade de meio contínuos tais como: transferência de calor,
eletrostáticos, meio porosos, etc. Como uma de suas realizações, Zienkiewicz (1968) publicou
sobre a análise de fissuras em meios elásticos utilizando o processo iterativo do método.
Szabó e Lee (1969) verificaram a possibilidade da formulação do método a partir de
equações diferenciais e condições de contorno, através do desenvolvimento da formulação
fraca, tal como o método de Galerkin (método de resíduos ponderados). Lynn e Arya (1973),
também através do método dos resíduos ponderados, no caso o de mínimos quadrados,
formularam o método dos elementos finitos.
Através das formulações direta, variacional ou de resíduos, retiram-se as equações
algébricas que regem o comportamento aproximado dos subdomínios, monta-se o sistema de
equações considerando todos os subdomínios (chamados por malha) no qual é denominado
sistema global. Juntamente com as condições de contorno nesse sistema global, é possível
determinar os valores nodais de definição desse campo. Torna-se possível agora o caminho
inverso onde é pode-se determinar incógnitas referentes ao elemento isolado, resolvendo para
qualquer um de seus pontos [8].
2.2. FUNDAMENTOS DO MÉTODO DE ELEMENTOS FINITOS
O desenvolvimento de peças e estruturas tanto em projetos da Engenharia Mecânica
ou Engenharia Civil, sempre demandam de muitos cálculos e elaboração de protótipos em
menor escala para a realização de testes com a finalidade de se conhecer o comportamento do
objeto de estudo e as possíveis falhas no corpo em seu período de vida útil.
Porém, todo esse processo realizado sem a contribuição computacional, torna-se
extremamente oneroso uma vez que existem custos embutidos na fabricação de modelos e
pelo tempo dispensado nos repetitivos ensaios e análise/comparação dos resultados. Além do
fato de se fazer necessário que exista uma larga experiência por parte dos Engenheiros
25
Projetistas para que seja feita uma análise criteriosa e exata da peça/estrutura de acordo com
as suas aplicações.
O método mais conhecido e utilizado é o Método dos Elementos Finitos (MEF) que se
utiliza de uma resolução integral englobando o domínio desejado para determinadas
condições de contorno, de acordo com o problema abordado.
O MEF tem como princípio a divisão do domínio (estrutura ou peça estudada) em
subdomínios denominados elementos (no geral, são figuras geométricas simples), para a
resolução de equações diferenciais de forma mais simples. Esse processo é denominado
discretização e gera um conjunto de subdomínios comumente chamado de malha, como pode
ser observado na figura 6, onde em (a) é ilustrado um corpo sólido discretizado e em (b)
alguns exemplos de geometrias mais utilizadas na discretização de um corpo.
(a)
(b)
Figura 6 - (a) Corpo discretizado com elementos tetragonais (Fonte: Imagem do autor);
(b) Geometrias de elementos utilizadas [14]
26
Para o caso de estudo em questão, será utilizado o elemento linear, pois não existem
ações de outros tipos de carregamento além das forças de tração e compressão.
Generalizando a formulação do MEF, deve existir uma equação integral sobre um
domínio complexo de volume 𝑉, onde esta deve ser possível de ser substituída por um
somatório de integrais sobre subdomínios mais simples de volume 𝑉𝑖, como por ser observado
na integral da função 𝑤 abaixo.
∫ 𝑤 𝑑𝑉𝑉
= ∑∫ 𝑤 𝑑𝑉𝑉
𝑛
𝑖=1
(2.1)
Onde se pressupõe que:
𝑉 = ∑𝑉𝑖
𝑛
𝑖=1
(2.2)
Sendo possível o cálculo das integrais nos elementos (subdomínios), basta então fazer
o somatório das dessas soluções para se obter a solução de todo o corpo (domínio).
2.3. DISCRETIZAÇÃO DO DOMÍNIO
Para a solução numérica da formulação fraca desenvolvida na seção 1.4, é necessário
que exista uma aproximação definida por combinações lineares para a aproximação das
funções 𝑢(𝑥), 𝑣(𝑥) 𝑒 𝑓(𝑥) que pode ser vista na Eq. (2.3) [10].
𝜓(𝑥) =∑Ψ𝑖(𝑥)𝜓𝑖
𝑛
𝑖=1
(2.3)
Onde:
27
𝜓𝑖 - são constantes a serem determinadas;
Ψ𝑖 - são as funções de forma conhecidas e linearmente independentes;
𝑛 – nº de graus de liberdade.
Definidas por polinômios do primeiro grau, as funções Ψ𝑖 valem 1 no ponto escolhido
e 0 nos demais. Multiplicando pelas constantes ainda indeterminadas 𝜓𝑖, resulta em uma
função de aproximação linear por partes, como representado na Figura 7.
Figura 7 - Gráfico das funções de forma (Fonte: Imagem do autor).
Tem-se então para o elemento, a função de aproximação com dois graus de liberdade
para a função de deslocamento:
𝑢(𝑥) = Ψ1(𝑥)𝑢1 + Ψ2(𝑥)𝑢2 (2.4)
Onde:
Ψ1 = (1 −𝑥
𝐿)
Ψ2 =𝑥
𝐿
Para fim de substituição na formulação integral, deriva-se a Eq. (2.4) em função de 𝑥,
obtendo-se:
28
𝑑𝑢(𝑥)
𝑑𝑥= Ψ1
′(𝑥)𝑢1 +Ψ2′(𝑥)𝑢2 (2.5)
Onde:
Ψ1′ = −
1
𝐿
Ψ2′ =
1
𝐿
Analogamente, são obtidas as funções de aproximação para 𝑣(𝑥) e 𝑓(𝑥) onde as
mesmas são substituídas na formulação integral, de domínio Ω [0; 𝐿], para um único elemento
onde serão obtidas as matrizes de influência.
Como a formulação integral obtida na Eq. (1.10) trata da soma de diferentes
influências de forças no corpo, a substituição das funções de aproximação pode ser realizada
de forma isolada para cada integral.
Considerando o corpo como um todo, verifica-se que o termo referente às forças de
fronteira não se aplica para essas condições e pode ser então eliminado, deixando a Eq. (1.10)
na forma abaixo:
∫ 𝜌𝐴𝜕2𝑢
𝜕𝑡2𝑣(𝑥) 𝑑𝑥
𝐿
0
+∫ 𝐸𝐴𝜕𝑢
𝜕𝑥
𝐿
0
𝑑𝑣
𝑑𝑥 𝑑𝑥 = ∫ 𝑓(𝑥)𝑣(𝑥) 𝑑𝑥
𝐿
0
(2.6)
Substituindo inicialmente as funções de aproximação para as forças de trabalho
interno de deformação e considerando 𝐸𝐴 constante, tem-se:
𝐸𝐴∫𝜕𝑢
𝜕𝑥
𝐿
0
𝑑𝑣
𝑑𝑥 𝑑𝑥 = 𝐸𝐴∫ (Ψ1
′𝑢1 +Ψ2′𝑢2)(Ψ1
′𝑣1 +Ψ2′𝑣2)
𝐿
0
𝑑𝑥 (2.7)
Podendo ser escrito na forma matricial:
𝐸𝐴∫𝜕𝑢
𝜕𝑥
𝐿
0
𝑑𝑣
𝑑𝑥 𝑑𝑥 = 𝐸𝐴
[ ∫ Ψ1
′2𝑑𝑥𝐿
0
∫ Ψ1′Ψ2′𝑑𝑥
𝐿
0
∫ Ψ2′
𝐿
0
Ψ1′𝑑𝑥 ∫ Ψ2
′2𝐿
0
𝑑𝑥]
[
𝑢1
𝑢2
] [
𝑣1
𝑣2
] (2.8)
29
Solucionando as integrais da primeira matriz, tem-se então a matriz de rigidez [𝑘𝑒] do
elemento, onde:
[𝑘𝑒] =𝐸𝐴
𝐿[1 −1−1 1
] (2.9)
A seguir, fazendo a substituição das funções de aproximação para as forças de inércia
onde o termo 𝜌𝐴 é constante e 𝜕2𝑢 𝜕𝑡2 = �̈�⁄ , tem-se que:
𝜌𝐴∫ �̈�(𝑥)𝑣(𝑥) 𝑑𝑥𝐿
0
= 𝜌𝐴∫ (Ψ1�̈�1 +Ψ2�̈�2)(Ψ1𝑣1 + Ψ2𝑣2)𝐿
0
𝑑𝑥 (2.10)
Escrevendo na forma matricial,
𝜌𝐴∫ �̈�(𝑥)𝑣(𝑥) 𝑑𝑥𝐿
0
= 𝜌𝐴
[ ∫ Ψ1
2𝑑𝑥𝐿
0
∫ Ψ1Ψ2𝑑𝑥𝐿
0
∫ Ψ2Ψ1𝑑𝑥𝐿
0
∫ Ψ22
𝐿
0
𝑑𝑥]
[
�̈�1
�̈�2
] [
𝑣1
𝑣2
] (2.11)
Solucionando as integrais da primeira matriz, tem-se então a matriz de massa [𝑚𝑒] do
elemento, onde:
[𝑚𝑒] =𝜌𝐴𝐿
6[2 11 2
] (2.12)
E por último, utilizando a integral referente às forças de corpo e fazendo as devidas
substituições, tem-se:
∫ 𝑓(𝑥)𝑣(𝑥) 𝑑𝑥𝐿
0
= ∫ (Ψ1𝑓1 +Ψ2𝑓2)(Ψ1𝑣1 +Ψ2𝑣2)𝐿
0
𝑑𝑥 (2.13)
Na forma matricial:
30
∫ 𝑓(𝑥)𝑣(𝑥) 𝑑𝑥𝐿
0
=
[ ∫ Ψ1
2𝑑𝑥𝐿
0
∫ Ψ1Ψ2𝑑𝑥𝐿
0
∫ Ψ2Ψ1𝑑𝑥𝐿
0
∫ Ψ22
𝐿
0
𝑑𝑥]
[
𝑓1
𝑓2
] [
𝑣1
𝑣2
] (2.14)
Solucionando as integrais da matriz, tem-se então:
∫ 𝑓(𝑥)𝑣(𝑥) 𝑑𝑥𝐿
0
=𝐿
6[2 11 2
] [𝑓1𝑓2] [𝑣1𝑣2] (2.15)
Tomando 𝑓(𝑥) = constante, tem-se que:
�⃗� =𝐿𝑓
2[11] = [
𝐹1𝐹2] (2.16)
2.4. MONTAGEM DO SISTEMA GLOBAL
2.4.1 Matrizes do elemento
Até o momento, os resultados encontrados para matriz de rigidez e massa foram
obtidos para apenas um elemento isolado com dois graus de liberdade locais, isto é, graus de
liberdade referentes apenas ao elemento estudado.
É necessário agora obter as matrizes de influência com relação ao referencial global x-
y para posteriormente considerar todos os demais elementos em uma mesma matriz de
influência de forma que todos estejam com o mesmo referencial.
31
Figura 8 - Elemento com relação ao referencial local e global (Fonte: Imagem do autor).
Como pode ser observado na figura 8, um elemento com dois graus de liberdade cujo
sistema de coordenadas local é dado por x’-y’ e o mesmo elemento com quatro graus de
liberdade cujo sistema de coordenadas é dado por x-y. Sendo assim, os deslocamentos para
sistema de coordenas local e global são dados respectivamente por:
𝑢′𝑇 = [𝑢′1 𝑢′2] (2.17)
𝑢𝑇 = [𝑢1 𝑢2 𝑢3 𝑢4] (2.18)
É possível observar a relação geométrica entre os eixos coordenados e ver que os
elementos de deslocamento locais podem ser dados por:
𝑢′1 = 𝑢1 cos(θ) + 𝑢2 sen(θ) (2.19)
𝑢′2 = 𝑢3 cos(θ) + 𝑢4 sen(θ) (2.20)
Onde θ é o ângulo entre os eixos x, local e global. Chega-se então à forma:
𝑢′ = 𝛽𝑢 (2.21)
Onde 𝛽 é a matriz de transformação dada por:
32
𝛽 = [𝑙 𝑚 0 00 0 𝑙 𝑚
] (2.22)
Sendo 𝑙 e 𝑚 dados pela forma generalizada:
𝑙 = 𝑐𝑜𝑠 𝜃𝑥 =𝑢3 − 𝑢1𝐿𝑒
=𝑋𝐹 − 𝑋𝑁
√(𝑋𝐹 − 𝑋𝑁)2 + (𝑌𝐹 − 𝑌𝑁)2 (2.23)
𝑚 = 𝑐𝑜𝑠 𝜃𝑦 =𝑢4 − 𝑢2𝐿𝑒
=𝑌𝐹 − 𝑌𝑁
√(𝑋𝐹 − 𝑋𝑁)2 + (𝑌𝐹 − 𝑌𝑁)2 (2.24)
𝐿𝑒 = √(𝑋𝐹 − 𝑋𝑁)2 + (𝑌𝐹 − 𝑌𝑁)2 (2.25)
Utilizando o conceito de energia elástica do elemento [4], a matriz de rigidez do
elemento é dada por:
𝐾𝑒 = 𝛽𝑇𝑘𝑒𝛽 =
𝐸𝐴
𝐿𝑒[
𝑙2 𝑙𝑚 −𝑙2 −𝑙𝑚𝑙𝑚−𝑙2
𝑚2
−𝑙𝑚−𝑙𝑚𝑙2
−𝑚2
𝑙𝑚−𝑙𝑚 −𝑚2 𝑙𝑚 𝑚2
] (2.26)
Aplicando a mesma teoria para a matriz de massa, tem-se:
𝑀𝑒 = 𝛽𝑇𝑚𝑒𝛽 =
𝜌𝐴𝐿𝑒6[
2𝑙2 2𝑙𝑚 𝑙2 𝑙𝑚2𝑙𝑚𝑙2
2𝑚2
𝑙𝑚𝑙𝑚2𝑙2
𝑚2
2𝑙𝑚𝑙𝑚 𝑚2 2𝑙𝑚 2𝑚2
] (2.27)
2.4.2 Matrizes globais
É necessário agora estender o processo aos demais elementos da estrutura de treliças, e
através dos resultados individuais para as matrizes de influência, utilizar uma metodologia
para a montagem de um sistema que englobe todos os elementos da geometria discretizada e
seja possível chegar às matrizes de rigidez e massa globais, chegando assim à solução do
problema.
33
Para a demonstração da montagem do sistema global, será considerada a configuração
de treliças ilustrada na figura 9 onde existe uma força F aplicada no nó 2, um engaste no nó 1
e um apoio simples no nó 3. Além disso, todos os elementos contam com módulo de
elasticidade do material (E) e área de seção transversal (A), iguais e constantes.
Figura 9 - Exemplo de estrutura de treliças (Fonte: Imagem do autor)
A estrutura é então discretizada onde seus nós/elementos são enumerados e os graus
de liberdade com relação ao referencial global são discriminados em cada nó, como
apresentado na figura 10 abaixo:
Figura 10 - Estrutura discretizada com 3 elementos e 6 graus de liberdade (Fonte:
Imagem do autor).
34
Obtêm-se então as matrizes de rigidez e massa para cada elemento, utilizando-se dos
métodos apresentados até o momento. Para fim de demonstração da montagem da matriz
global de rigidez (𝐾𝐺) e massa (𝑀𝐺), foram usados valores genéricos em cada termo das
matrizes onde, no final, será apresentada a influência de cada nó na matriz global.
1) Elemento 1:
Figura 11 - Primeiro elemento (Fonte: Imagem do autor).
1 2 3 4
𝐾1 = 𝐸𝐴 [
𝐴11 𝐴12 𝐴13 𝐴14𝐴21𝐴31
𝐴22𝐴32
𝐴23𝐴33
𝐴24𝐴34
𝐴41 𝐴42 𝐴43 𝐴44
]
1234
(2.28)
2) Elemento 2:
35
Figura 12 - Segundo elemento (Fonte: Imagem do autor).
3 4 5 6
𝐾2 = 𝐸𝐴 [
𝐵33 𝐵34 𝐵35 𝐵36𝐵43𝐵53
𝐵44𝐵54
𝐵45𝐵55
𝐵46𝐵56
𝐵63 𝐵64 𝐵65 𝐵66
]
3456
(2.29)
3) Elemento 3:
Figura 13 - Terceiro elemento (Fonte: Imagem do autor).
𝟏 𝟐 𝟓 𝟔
𝐾3 = 𝐸𝐴 [
𝐶11 𝐶12 𝐶15 𝐶16𝐶21𝐶51
𝐶22𝐶52
𝐶25𝐶55
𝐶26𝐶56
𝐶61 𝐶62 𝐶65 𝐶66
]
1256
(2.30)
36
As matrizes globais dos elementos devem ser somadas e, como foi mostrado, cada
coluna/linhas representa um grau de liberdade. Como a estrutura do exemplo possui seis graus
de liberdade, suas matrizes de influência globais serão quadradas de ordem 6. Somando os
termos das Eq. (2.28), (2.29) e (2.30), obtém-se:
𝐾𝐺 1 2 3 4 5 6
=𝐸𝐴
𝐿𝑒
[ 𝐴11 + 𝐶11𝐴21 + 𝐶21𝐴31𝐴41𝐶51𝐶61
𝐴12 + 𝐶12𝐴22 + 𝐶22𝐴32𝐴42𝐶52𝐶62
𝐴13𝐴23
𝐴33 + 𝐵33𝐴43 + 𝐵43𝐵53𝐵63
𝐴14𝐴24
𝐴34 + 𝐵34𝐴44 + 𝐵44𝐵54𝐵64
𝐶15𝐶25𝐵35𝐵45
𝐵55 + 𝐶55𝐵65 + 𝐶65
𝐶16𝐶26𝐵36𝐵46
𝐵56 + 𝐶56𝐵66 + 𝐶66]
123456
(2.31)
Analogamente para matriz de massa, repetindo o processo para os três elementos,
chega-se à matriz de massa global:
𝑀𝐺 1 2 3 4 5 6
=𝜌𝐴𝐿𝑒6
[ 𝐴11 + 𝐶11𝐴21 + 𝐶21𝐴31𝐴41𝐶51𝐶61
𝐴12 + 𝐶12𝐴22 + 𝐶22𝐴32𝐴42𝐶52𝐶62
𝐴13𝐴23
𝐴33 + 𝐵33𝐴43 + 𝐵43𝐵53𝐵63
𝐴14𝐴24
𝐴34 + 𝐵34𝐴44 + 𝐵44𝐵54𝐵64
𝐶15𝐶25𝐵35𝐵45
𝐵55 + 𝐶55𝐵65 + 𝐶65
𝐶16𝐶26𝐵36𝐵46
𝐵56 + 𝐶56𝐵66 + 𝐶66]
123456
(2.32)
Que são as matrizes globais de rigidez e massa da estrutura considerando todos seus
elementos.
2.5. CONDIÇÕES DE CONTORNO
Uma das primeiras etapas na resolução de qualquer problema que envolva equilíbrio
de forças de um corpo é a determinação das posições das reações referentes às condições de
apoio ao qual o corpo está exposto. Destacar as condições de apoio é o que possibilita a
37
obtenção dos deslocamentos dos demais nós expostos a essas condições, tanto de apoios
quanto de forças.
Para os apoios das estruturas de treliças, existem três tipos de apoio como pode ser
visto na figura 14 abaixo:
Figura 14 - Reações em apoio para uma estrutura bidimensional (Fonte: BEER – 2012)
Para fins de aplicação no código numérico, aqueles graus de liberdade com restrições,
isto é, que não influenciam no cálculo, são então zerados e por consequência, todos os itens
referentes a eles são eliminados. Recorrendo novamente ao exemplo da figura 9, pode ser
observado que nos graus de liberdade 1, 2 e 6 não há movimento, pois a estrutura está
engastada no nó 1 e simplesmente apoiada no nó 3. Sendo assim, como já anteriormente
mencionado, todas as linhas e colunas das matrizes de rigidez, massa e carregamento desses
graus de liberdade são retirados, transformando as mesmas na forma a seguir:
�⃗⃗� =
[ 𝑢1 = 0𝑢2 = 0𝑢3=?𝑢4=?𝑢5=?𝑢6 = 0]
123456
38
1 2 3 4 5 6
𝐾𝐺 =𝐸𝐴
𝐿𝑒
[ 𝐴11 + 𝐶11𝐴21 + 𝐶21𝐴31𝐴41𝐶51𝐶61
𝐴12 + 𝐶12𝐴22 + 𝐶22𝐴32𝐴42𝐶52𝐶62
𝐴13𝐴23
𝐴33 + 𝐵33𝐴43 + 𝐵43𝐵53𝐵63
𝐴14𝐴24
𝐴34 + 𝐵34𝐴44 + 𝐵44𝐵54𝐵64
𝐶15𝐶25𝐵35𝐵45
𝐵55 + 𝐶55𝐵65 + 𝐶65
𝐶16𝐶26𝐵36𝐵46
𝐵56 + 𝐶56𝐵66 + 𝐶66]
123456
1 2 3 4 5 6
𝑀𝐺 =𝜌𝐴𝐿𝑒
6
[ 𝐴11 + 𝐶11𝐴21 + 𝐶21𝐴31𝐴41𝐶51𝐶61
𝐴12 + 𝐶12𝐴22 + 𝐶22𝐴32𝐴42𝐶52𝐶62
𝐴13𝐴23
𝐴33 + 𝐵33𝐴43 + 𝐵43𝐵53𝐵63
𝐴14𝐴24
𝐴34 + 𝐵34𝐴44 + 𝐵44𝐵54𝐵64
𝐶15𝐶25𝐵35𝐵45
𝐵55 + 𝐶55𝐵65 + 𝐶65
𝐶16𝐶26𝐵36𝐵46
𝐵56 + 𝐶56𝐵66 + 𝐶66]
123456
�⃗� =
[ 𝐹1𝐹2𝐹3𝐹4𝐹5𝐹6]
=
[ 𝐹1 =?𝐹2 =?𝐹3 = 0𝐹4 = −𝐹𝐹5 = 0𝐹6 =? ]
123456
A partir daqui a solução prossegue como descrito na equação do movimento na forma
matricial, Eq. (1.11).
39
3 RESULTADOS
Para a demonstração do funcionamento do algoritmo, é fornecido ao código o
posicionamento de cada nó no espaço (através do arquivo “nodes.txt”) e os nós
correspondentes a cada elemento, módulo de elasticidade, área de seção transversal e
densidade (através do arquivo “elements.txt”).
3.1. EXEMPLO 1
(a) (b)
Figura 15 – (a) Estrutura do exemplo 1 (b) Estrutura discretizada (Fonte: Imagem do
autor).
Escolhendo um material para o exemplo da seção anterior, no caso o alumínio 7075-
T6, tem-se então que 𝐸 = 70𝐺𝑃𝑎, 𝐴 = 3𝑥10−4𝑚² e 𝜌 = 2800𝑘𝑔/𝑚3. Para o carregamento
descrito na figura 15, onde seus dados de entrada no algoritmo podem ser vistos no Apêndice
desse arquivo. São determinadas também as condições de contorno para o vetor força (�⃗�) e os
pontos de apoio (�⃗⃗�) como mostrados abaixo:
40
�⃗� =
[ 𝐹1𝐹2𝐹3𝐹4𝐹5𝐹6]
=
[
000
−1000000 ]
�⃗⃗� =
[ 00
𝑢3 =?𝑢4 =?𝑢5 =?0 ]
3.1.1 ESTÁTICO: [𝐾]�⃗⃗� = �⃗�
Por meio do algoritmo obtêm-se as matrizes de rigidez global (K) e como resultado, o
vetor deslocamento (U):
Resultado:
O resultado encontrado apresenta os deslocamentos em metros para cada grau de
liberdade. Fornecendo ao projetista uma boa ferramenta de análise da estrutura para as
condições de apoio e carregamento ao qual ela está exposta.
41
3.1.2 DINÂMICO: ([𝐾] − 𝜔2[𝑀])�⃗� = 0⃗⃗
Com a finalidade de obter a solução dinâmica, o código deve fornecer a matriz de
massa global (M), uma vez que a matriz de rigidez já foi obtida, e assim obter o vetor de
autovalores (e) e seus autovetores correspondentes. Sabendo que são através dos autovalores
que são determinadas as frequências naturais (Wn) da estrutura, tem-se então:
Resultado:
O resultado encontrado apresentaria as frequências naturais da estrutura em radianos
por segundo e tendo a estrutura seis graus de liberdade, a mesma possuiria seis modos de
vibração.
42
Os modos de vibração fornecem uma boa ferramenta de análise, pois caso a estrutura
esteja sujeita a carregamentos cíclicos que operem na mesma frequência que alguma de suas
frequências naturais, o sistema pode entrar em ressonância. Sendo assim, o projetista pode
alterar o carregamento aplicado ou modificar propriedades da estrutura para solucionar o
problema.
Porém, a matriz de massa deve conter apenas valores positivos, o que não foi obtido
com esse código. Sendo assim, foram encontrados valores negativos para os autovalores
acarretando assim um resultado fisicamente impossível contendo números complexos para os
modos de vibração, o que não permite analisar o problema dinamicamente.
3.2. EXEMPLO 2
(a)
(b)
Figura 16 – (a) Estrutura do exemplo 2; (b) Estrutura discretizada (Fonte: Imagem do
autor).
43
Da mesma forma que no exemplo anterior, considerando uma estrutura construída em
aço SAE 4140, cujos valores para o módulo de elasticidade, área da seção transversal e
densidade do material são dados, respectivamente, por: 𝐸 = 200𝐺𝑃𝑎, 𝐴 = 1𝑥10−4𝑚² e
𝜌 = 7800𝑘𝑔/𝑚3.
Para o carregamento descrito na figura 16, onde seus dados de entrada no algoritmo
podem ser vistos no Apêndice desse arquivo. São determinadas também as condições de
contorno para o vetor força (�⃗�) e os pontos de apoio (�⃗⃗�) como mostrados abaixo:
�⃗� =
[ 𝐹1𝐹2𝐹3𝐹4𝐹5𝐹6𝐹7𝐹8𝐹9𝐹10𝐹11𝐹12𝐹13𝐹14𝐹15𝐹16𝐹17𝐹18𝐹19𝐹20𝐹21𝐹22]
=
[
0000000
−1260000000000
−1260000000
720000 ]
�⃗⃗� =
[
0
0𝑢3 =?𝑢4 =?𝑢5 =?𝑢6 =?𝑢7 =?𝑢8 =?𝑢9 =?𝑢10 =?𝑢11 =?𝑢12 =?𝑢13 =?𝑢14 =?𝑢15 =?𝑢16 =?𝑢17 =?0
𝑢19 =?𝑢20 =?𝑢21 =?𝑢22 =?]
3.2.1 ESTÁTICO: [𝐾]�⃗⃗� = �⃗�
Novamente, por meio do algoritmo obtêm-se as matrizes de rigidez global (K) e como
resultado, o vetor deslocamento (U):
44
45
46
Resultado:
Da mesma forma que no exemplo anterior, o resultado encontrado apresenta os
deslocamentos em metros para cada grau de liberdade.
3.2.2 DINÂMICO: ([𝐾] − 𝜔2[𝑀])�⃗� = 0⃗⃗
Com a finalidade de obter a solução dinâmica, o código deve fornecer a matriz de
massa global (M), uma vez que a matriz de rigidez já foi obtida, e assim obter o vetor de
47
autovalores (e) e seus autovetores correspondentes. Sabendo que são através dos autovalores
que são determinadas as frequências naturais (Wn) da estrutura, tem-se então:
48
Porém, como se trata do mesmo código, a matriz de massa apresentou valores
negativos, o que não deveria acontecer. Sendo assim, foram encontrados valores negativos
para os autovalores impossibilitando a análise dinâmica do caso.
49
4 CONCLUSÕES
O objetivo geral desse trabalho foi estudar o comportamento de quaisquer
configurações de treliças planas com quaisquer valores de carregamentos atuantes, através de
uma abordagem numérica, no caso, o método de elementos finitos para assim elaborar um
código no software MATLAB R2014a.
A análise abrange aspectos estáticos da estrutura, isto é, obtenção dos deslocamentos
de suas barras, e dinâmicos que trata da obtenção de seus modos de vibração.
Através da abordagem numérica, pôde-se observar detalhadamente o funcionamento
do método de elementos finitos, desde a discretização da estrutura, até a obtenção das
matrizes de influência. Passando pela transformação das mesmas matrizes com relação ao
referencial local para o referencial global.
Sendo assim, o código fornece como resultados o vetor deslocamento, para uma
análise estática, e o vetor contendo de frequências naturais, para o caso dinâmico.
Para o vetor deslocamento, o algoritmo, que designa o mesmo por “U”, exibe o
deslocamento do nó nos seus graus de liberdade correspondentes e assim disponibiliza ao
avaliador um panorama geral do comportamento da estrutura devido aos carregamentos
atuantes.
Já para o vetor dos modos de vibração, que no algoritmo é designado por “Wn”, cada
modo de vibração seria dado por um elemento desse vetor. Porém, um dos itens que era
levado em consideração para a obtenção dos modos de vibração, no caso a matriz de massa
designada por “M”, não foi obtida com sucesso, o que fez com que não fosse possível chegar
ao resultado final para a parte dinâmica. Ainda sim, fazer uma avaliação dinâmica é de suma
importância para que o projetista tome decisões sobre a utilização de equipamentos que
operem na estrutura e procure evitar aqueles que trabalham na mesma frequência, não
permitindo que o sistema de treliças entre em ressonância com o equipamento, pois do
contrário, isso causaria problemas posteriores.
Sendo assim, visto que o algoritmo fornece grande arte dos dados propostos, pode se
observar sua eficiência e também a velocidade na análise desse tipo de estrutura.
Também é possível sua utilização por parte de empresas de pequeno e médio porte,
por se tratar de uma versão livre a custo zero, uma vez que licenças de softwares conhecidos,
como ANSYS ou MSC Nastran, são extremamente dispendiosas e impedem que pequenas e
médias indústrias tenham condição de absorver tais custos.
50
Por fim, a utilização da computação numérica agiliza os processos e os tornam mais
confiáveis aliados a uma visão crítica do Engenheiro Projetista.
51
5 REFERÊNCIAS BIBLIOGRÁFICAS
[1] BEER, Ferdinand P.; JOHNSTON, E. Russell. “Mecânica Vetorial para Engenheiros-
Estática”. AMGH Editora Ltda. 9ª Edição, 2012.
[2] BEER, Ferdinand P.; JOHNSTON, E. Russell. “Resistência dos Materiais”. Editora
Pearson. 3ª Edição, 2012.
[3] BABUSKA, Ivo; SZABÓ, Barna. “Finite Element Analysis”. Editora Wiley Interscience,
1ª Edição, 1991.
[4] CAMPOS, Pedro Henrique T. “Um estudo sobre vibração em cordas”. Niterói, 2010.
[5] FERREIRA, Antonio J. M. “MATLAB Codes for Finite Element Analysis”. Editora
Springer, Volume 157, Portugal-Porto, 2008.
[6] RAO, Singiresu.S. “Vibrações Mecânicas”. Editora Pearson Prentice Hall, 4ª Edição,
2009.
[7] RODRIGUES, Adriano da C.; SILVA, Rodolpho C. “Modelagem de vibração de eixo sob
torção”. Niterói, 2013.
[8] SORIANO, Humberto L.; LIMA, Silvio de S. “Método de Elementos Finitos em Análise
de Estruturas”. Editora da Universidade de São Paulo, São Paulo, 2003.
[9] TRINDADE, Marcelo A.; SAMPAIO, Rubens. “Apostila de MATLAB”. Rio de Janeiro,
2002.
[10] PEREIRA, Orlando J. B. A. “Análises de Estruturas II: Introdução ao Método dos
Elementos Finitos na Análise de Problemas Planos de Elasticidade”. 2005.
[11] Disponível em: <http://www.metalfabrication.com.pt/newproduct/4-2-2b.jpg>; Acesso
em 30 de junho de 2015.
[12] Disponível em: <www.comprar-vender.mfrural.com.br/detalhe/estrutura-metalica-
galpao-65023.aspx>; Acesso em 11 de junho de 2015.
[13] Disponível em: <www.portuguese.overhead-cranehoist.com/china-
qme120t_78m_65m_truss_girder_outdoor_long_span_gantry_crane-192577.html>; Acesso
em 11 de junho de 2015.
[14] Disponível em: <mingaonline.uach.cl/scielo.php?pid=S0718-
025X2005000100004&script=sci_arttext>; Acesso em 15 de julho de 2015.
52
6 APÊNDICES
6.1. DADOS DE ENTRADA
6.1.1 Exemplo 1:
Conteúdo do arquivo “elements.txt”:
Nó inicial Nó final E A 𝜌
1 2 70000000000 0.0003 2800
2 3 70000000000 0.0003 2800
3 1 70000000000 0.0003 2800
Conteúdo do arquivo “nodes.txt”:
Posição do nó em x Posição do nó em y
0 0
0 3
3 0
6.1.2 Exemplo 2:
Conteúdo do arquivo “elements.txt”:
Nó inicial Nó final E A 𝜌
1 2 200000000000 0.0001 7800
1 3 200000000000 0.0001 7800
1 4 200000000000 0.0001 7800
2 4 200000000000 0.0001 7800
3 4 200000000000 0.0001 7800
3 5 200000000000 0.0001 7800
3 6 200000000000 0.0001 7800
53
4 6 200000000000 0.0001 7800
5 6 200000000000 0.0001 7800
5 7 200000000000 0.0001 7800
5 8 200000000000 0.0001 7800
6 8 200000000000 0.0001 7800
7 8 200000000000 0.0001 7800
7 9 200000000000 0.0001 7800
7 10 200000000000 0.0001 7800
8 10 200000000000 0.0001 7800
9 10 200000000000 0.0001 7800
9 11 200000000000 0.0001 7800
10 11 200000000000 0.0001 7800
Conteúdo do arquivo “nodes.txt”:
Posição do nó em x Posição do nó em y
0 0
0 3
2.4 0
2.4 3
4.8 0
4.8 3
7.2 0
7.2 3
9.6 0
9.6 3
12 3
54
6.2. CÓDIGO FONTE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Exemplo 2 % %Programa para Análise por Elementos Finitos em MATLAB de Treliças Planas% %%
clear all; clc; close all; format long;
%%%%%%%%%%%%%%%%%%%%%%%%%%% Entrada de dados %%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nodes = load('nodes.txt'); %coordenadas dos nós: x, y Elements = load('elements.txt'); % nós dos elementos: 1º nó, 2º nó, E, A numnode = length(Nodes); %Detém o número de nós numelem = size(Elements,1); %Detém o número de elementos
%%%%%%%%%%% Declaração das Matrizes de Influência [K, U, F, M] %%%%%%%%%%% K = zeros(2*numnode,2*numnode); %Matriz total de GL x total de GL U = zeros(2*numnode, 1); %[Matriz de deslocamento] total de GL x 1 F = zeros(2*numnode, 1); %[Matriz de carregamento] total de GL x 1 M = zeros(2*numnode,2*numnode); %Declara a :[Matriz de massa] total de GL
x total de GL e = zeros (numnode,1);
%%%%%%%%%%%%%% Condições de Contorno de Apoio e Carregamento %%%%%%%%%%%%% act = 1:2*numnode; %GL ativos act([1 2 18]) = []; %GL zerados F(8) = -126000; %Força aplicada no GL indicado entre parênteses F(16) = -126000; F(21) = 72000;
55
%%%%%%%%%%%% Geração das Matrizes de Rigidez e Massa Globais %%%%%%%%%%%%% for ie = 1:numelem DOFs = [2*Elements(ie, 1)-1, 2*Elements(ie, 1), 2*Elements(ie, 2)-1,
2*Elements(ie, 2)]; %Obtém os elementos dos GL X1 = Nodes(Elements(ie,1), 1); Y1 = Nodes(Elements(ie,1), 2); X2 = Nodes(Elements(ie,2), 1); Y2 = Nodes(Elements(ie,2), 2);
L = sqrt((X2-X1)^2+(Y2-Y1)^2); %Comprimento do elemento ms = (Y2-Y1)/L ; ls=(X2-X1)/L; E = Elements(ie,3); %Módulos de elasticidade do elemento A = Elements(ie,4); %Áreas de seção transversal do elemento rho= Elements(ie,5); %Densidade do elemento K(DOFs,DOFs) = K(DOFs,DOFs) + (E*A/L)*[ls^2 ls*ms -ls^2 -ls*ms; ls*ms
ms^2 -ls*ms -ms^2; -ls^2 -ls*ms ls^2 ls*ms; -ls*ms -ms^2 ls*ms ms^2]; %
Calcula a matriz de rigidez do elemento e monta na matriz de rigidez
global M(DOFs,DOFs) = M(DOFs,DOFs) + ((rho*A*L)/6)*[2*ls^2 2*ls*ms ls^2
ls*ms; 2*ls*ms 2*ms^2 ls*ms ms^2; ls^2 ls*ms ls^2 2*ls*ms; ls*ms ms^2
2*ls*ms 2*ms^2]; % Calcula a matriz de massa do elemento e monta na matriz
de massa global
end
%%%%%%%%%%%%%%%%%%%% Obtenção do Resultado Estático %%%%%%%%%%%%%%%%%%%%%% U(act) = K(act,act)\F(act); %Vetor deslocamento
%%%%%%%%%%%%%%%%%%%% Obtenção do Resultado Dinâmico %%%%%%%%%%%%%%%%%%%%%% e = eig((M^(-1/2))*K*(M^(-1/2))); %Autovalores [autove autoval]=eig((M^(-1/2))*K*(M^(-1/2)));%Autovetores e Autovalores Wn=sqrt(e); %Frequências naturais