Universidade Camilo Castelo Branco
Campus Fernandópolis
LUIS HENRIQUE DE REZENDE CROZARIOL
ANÁLISE LINEAR DE ESTRUTURAS PELO METÓDO DOS
ELEMENTOS FINITOS
LINEAR ANALYSIS OF STRUCTURES BY FINITE ELEMENT METHOD
Fernandópolis, SP
2014
Luis Henrique de Rezende Crozariol
ANÁLISE LINEAR DE ESTRUTURAS PELO METÓDO DOS ELEMENTOS FINITOS
Orientador: Prof. Me. Marcelo Rodrigo de Matos Pedreiro
Trabalho de Conclusão de Curso apresentada ao Curso de Graduação em Engenharia Civil da
Universidade Camilo Castelo Branco, como complementação dos créditos necessários para obtenção
do título de Graduação em Engenharia Civil.
Fernandópolis, SP
2014
Autorizo, exclusivamente, para fins acadêmicos e científicos, a reprodução total ou
parcial deste TCC, dissertação (tese), por processos xerográficos ou eletrônicos.
Assinatura do aluno:
Data:
Crozariol, Luis Henrique de Rezende
Cr953a Análise Linear de Estruturas pelo métodos dos
Elementos Finitos / Luis Henrique de Rezende.
Fernandópolis: [s.n.], 2014.
XV, 111p. : il. ; 29,5cm.
Trabalho de Conclusão de Curso, apresentado ao
Curso de Graduação em Engenharia Civil da Universidade
Camilo Castelo Branco, como complementação dos créditos
necessários para obtenção do título de Graduação em
Engenharia Civil.
Orientador: Profº. Me. Marcelo Rodrigo de Matos
Pedreiro.
1. Métodos dos Elementos Finitos. 2. Análise de Estruturas. 3. Métodos dos deslocamentos.
I. Título.
CDD – 624.1
v
Dedicatória
Dedico este trabalho aos meus pais,
Luiz Carlos e Irene, pela sabedoria que
me proporcionaram ao longo dos anos e
meu irmão Luis Gustavo.
Aos meus avós José Alcides e
Osmária “in memorian”, Vicente e Nadir
pelo carinho incondicional e por serem
exemplos de vida para mim. Vó Isabel
que infelizmente não conheci mas tenho
certeza que me olha do céu.
Minha namorada Mirian Matos e sua
mãe pelo incentivo e por acreditarem nos
meus objetivos.
Meus tios que são meus pais mais
novos, tio Adilson estou realizando um
sonho seu e gostaria de erguer este troféu
com você.
Meus primos, em especial Eloisa
Rezende e Carlos Eduardo por serem
meus irmãos mais velhos.
Aos meus amigos de infância, Éric,
Igor, Rafael, Rodrigo, Henrique pelo DotA
e pelos momentos felizes que passamos
juntos. Willian Queiroz, poucas pessoas
sabem, mas a maior rivalidade é
Flamengo x Palmeiras.
A galera do futebol, na época do Bagi
e Edinaldo, saudades desse tempo.
vi
Agradecimentos
O autor agradece o professor e orientador Marcelo Pedreiro, pela confiança
depositada em mim desde o início, paciência, dedicação, segurança e ensinamentos
no decorrer deste trabalho.
Aos professores Wilson Capanema, Edson Florentino e novamente Marcelo
Pedreiro, pelas aulas de Cálculo, Resistência dos materiais e Análise de estruturas,
matérias essenciais ao meu Trabalho de Conclusão de Curso.
Ao meu amigo Jhonata Olentino, pela ajuda nos trabalhos da faculdade e
formatação do TCC.
Agradeço a todas as pessoas que direta ou indiretamente contribuíram para
este trabalho.
vii
"Qualquer coisa que você aprende se
torna sua riqueza, uma riqueza que não
pode se tomada de você; seja se você
aprende em um prédio chamado escola
ou na escola da vida. Aprender algo novo
é um prazer e um tesouro valioso. E nem
todas as coisas que você aprende são
ensinadas a você, mas muitas coisas que
você aprende você percebe ter ensinado
a si mesmo.” - C. JoyBell C.
viii
RESUMO
Luis Henrique de Rezende Crozariol, Marcelo Rodrigo de Matos Pedreiro. Análise
linear de estruturas pelo método dos elementos finitos. Fernandópolis,
Universidade Camilo Castelo Branco, 2014, 107p. Trabalho de conclusão de curso.
No trabalho foi feita a revisão bibliográfica da base dos elementos finitos, composta
por: resistência dos materiais, geometria analítica, álgebra linear e cálculo diferencial
e integral. Após ter estudado a base, foi iniciado os estudos dos elementos finitos,
via apostilas e livros reconhecidos no meio acadêmico, assim como técnicas
computacionais para solução do sistema de equações algébricas.
A partir da revisão das matérias o Método dos Elementos Finitos fica compreensível,
por ser uma matéria ministrada em cursos de pós-graduação exige um
conhecimento avançado da base, é também observado a superioridade do Método
dos elementos finitos sobre o seu precursor, que é o método dos deslocamentos.
Palavras-chave: Método dos elementos finitos, Análise de estruturas, Método dos
deslocamentos.
ix
ABSTRACT
Luis Henrique de Rezende Crozariol, Marcelo Rodrigo de Matos Pedreiro. Linear
analysis of structures by finite element method. Fernandópolis, University Camilo
Castelo Branco, 2014, 107 pages. End of course work.
Work in the bibliographical review of the finite element basis, comprising been made:
strength of materials, analytical geometry, linear algebra and differential and integral
calculus. After studying the basic, the study was initiated finite element, by handouts
and books recognized in academia, as well as computational techniques for the
solution of algebraic equations. From the review of the materials the Finite Element
Method is understandable, because it is a given in courses of graduate field requires
an advanced knowledge base; it is also observed the superiority of the finite element
method over its precursor, which is the displacement method.
Key words: Finite element method, analysis of structures, the displacement method.
x
LISTA DE FIGURAS
Figura 1 - Modelo de estrutura contínua discretizada pelo método dos elementos
finitos. ........................................................................................................................ 18
Figura 2 - Viga elemento infinitesimal........................................................................ 21
Figura 3 - Sistema de eixos da Estática .................................................................... 22
Figura 4 – Treliça com 2 graus de liberdade. ............................................................ 24
Figura 5 – Equilíbrio do nó C. .................................................................................... 27
Figura 6 – Termos 11k e 21k da matriz de rigidez da treliça ....................................... 30
Figura 7 – Termos 12k e 22k da matriz de rigidez da treliça ....................................... 30
Figura 8 – Forças no nó C para 1 1d e 2 1d .......................................................... 32
Figura 9 – Graus de liberdade no sistema global e local ........................................... 35
Figura 10 – Energia de deformação específica 0U da barra m. ................................ 38
Figura 11 – Trabalho externo associado ao grau de liberdade i . ............................. 40
Figura 12 – Incremento de energia de deformação específica 0,mU da barra m . ... 45
Figura 13 – Incremento de trabalho externo iW . ..................................................... 49
Figura 14 – Viga em balanço de inércia variável. ...................................................... 56
Figura 15 – Diagrama de momentos na viga associado a v x definido em (135). .. 60
Figura 16 – Diagrama de momentos na viga associado a v x definido em (146). .. 62
Figura 17 – Elemento finito de viga. .......................................................................... 67
Figura 18 – Força externa linearmente distribuída. ................................................... 76
Figura 19 – Viga em balanço com força uniforme distribuída .................................... 79
Figura 20 – Força concentrada na extremidade ........................................................ 80
Figura 21 – Força linearmente distribuída com P1=0 ................................................ 80
Figura 22 – Força linearmente distribuída com P2=0 ................................................ 81
Figura 23 – Matriz de rigidez com característica de banda. ...................................... 82
Figura 24 – Matriz de rigidez sem característica de banda. ...................................... 82
Figura 25 – Numeração local dos deslocamentos do elemento. ............................... 83
Figura 26 – Perfil para armazenamento por altura efetiva de coluna. ....................... 86
Figura 27 - Fluxograma de obtenção do vetor “IPOS”. .............................................. 88
xi
Figura 28 - Fluxograma da montagem da matriz de rigidez por altura efetiva da
coluna. ....................................................................................................................... 89
Figura 29 - Fluxograma da montagem do vetor de forças nodais global. .................. 91
Figura 30 - Fluxograma – Etapa de triangularização. ................................................ 98
Figura 31 - Fluxograma – Etapa de substituição. ...................................................... 99
Figura 32 - Fluxograma - Etapa de retro substituição.............................................. 100
Figura 33 - Transformação dos deslocamentos nodais. .......................................... 102
Figura 34 - Tela de abertura e seleção do modelo. ................................................. 103
Figura 35 - Dados gerais. ........................................................................................ 104
Figura 36 - Dados nodais. ....................................................................................... 104
Figura 37 - Forças. .................................................................................................. 105
Figura 38 - Estrutura Modelada no Programa Computacional SAP2000. ............... 105
Figura 39 - Escolha do Tipo de Material. ................................................................. 106
Figura 40 - Elementos e suas características. ......................................................... 106
Figura 41 - Deslocamentos nodais - Parte 1. .......................................................... 107
Figura 42 - Deslocamentos nodais - Parte 2. .......................................................... 107
Figura 43 - Deslocamentos do Nó 4. ....................................................................... 108
xii
LISTA DE TABELAS
Tabela 1- Vetor auxiliar contendo as posições dos elementos da diagonal principal.
.................................................................................................................................. 86
Tabela 2 – Comparação de deslocamentos e rotação entre SAP e MAPE ............. 109
xiii
LISTA DE SÍMBOLOS
A = Área
i = incógnitas das funções aproximadoras
id = deslocamentos nodais
gd = deslocamentos nodais no sistema global
ld = deslocamentos nodais no sistema local
= deformação
E = módulo de elasticidade
= variável adimensional
= energia potencial total
f = matriz de forças nodais
= giro
I = inércia
k = matriz de rigidez
L = comprimento da barra
M = momento
N = força normal
= potencial das forças externas
P = carga
R = matriz de rotação
= incremento
i = alongamento/encurtamento
= tensão
0mU = energia de deformação específica
v x = flecha
'v x = rotação
W = trabalho
xiv
SUMÁRIO
1 – INTRODUÇÃO .................................................................................................... 17
2 - OBJETIVO ........................................................................................................... 19
3 – METODOLOGIA .................................................................................................. 20
4 – ELEMENTOS ESTRUTURAIS RETICULARES – VIGA PRISMÁTICA ............... 21
5 - FORMULAÇÃO LOCAL ....................................................................................... 22
5.1 - Dedução direta da equação diferencial regente do problema de viga
prismática .............................................................................................................. 22
6 – A EVOLUÇÃO DO MÉTODO DOS DESLOCAMENTOS .................................... 23
6.1 – Método Básico .............................................................................................. 24
6.2 – Método clássico ............................................................................................ 29
6.3 - Método da análise matricial ........................................................................... 34
6.4 – Método de Castigliano .................................................................................. 38
6.4.1 – Energia de deformação .......................................................................... 38
6.4.2 – Trabalho Externo .................................................................................... 40
6.4.3 – Segundo teorema de Castigliano ........................................................... 40
6.4.4 – A aplicação do método de Castigliano ................................................... 43
6.5 – Princípio dos deslocamentos virtuais ............................................................ 45
6.5.1 – Incrementos da energia de deformação ................................................. 45
6.5.2 – Incrementos do trabalho externo ............................................................ 48
6.5.3 – Formulação do princípio dos deslocamentos virtuais ............................. 51
6.6 – Método da mínima energia potencial total .................................................... 51
6.6.1 - O princípio da mínima energia potencial total ........................................ 52
6.7 - Método de Rayleigh-Ritz ............................................................................... 55
7 – O MÉTODO DOS ELEMENTOS FINITOS .......................................................... 65
7.1 - Dedução com utilização da linguagem matricial ............................................ 66
7.1.1 – Exemplo ................................................................................................. 73
xv
7.2 – Energia potencial externa ............................................................................. 75
7.3 – Exemplos – flechas ....................................................................................... 78
8 – TÉCNICAS COMPUTACIONAIS PARA AUTOMATIZAÇÃO DO MÉTODO DOS
ELEMENTOS FINITOS ............................................................................................. 81
8.1 – Organização da montagem do sistema de equações ................................... 81
8.1.1 – Generalidades ........................................................................................ 81
8.1.2 – Disposição dos coeficientes na matriz de rigidez global ........................ 81
8.1.3 – Processo de expansão e acumulação .................................................... 83
8.1.4 – Matriz de rotação .................................................................................... 84
8.1.5 – Armazenamento computacional da matriz de rigidez ............................. 85
8.1.6 – Armazenamento do vetor de forças nodais ............................................ 90
8.2 – Consideração das condições de contorno e dos deslocamentos prescritos . 92
8.2.1 – Introdução .............................................................................................. 92
8.2.2 – Técnicas para consideração das condições de contorno dos vínculos .. 92
8.2.3 – Técnica dos zeros e um ......................................................................... 93
8.2.4 – Técnica do número muito grande ........................................................... 94
8.2.5 – Apoios Elásticos ..................................................................................... 95
8.3 – Solução do sistema de equações ................................................................. 96
8.3.1 – Generalidades ........................................................................................ 96
8.3.2 - Procedimentos de solução ...................................................................... 96
8.3.3 - Implementação do método de solução para análise estática linear ........ 96
8.4 - Informações resultantes da análise ............................................................. 100
8.4.1 - Generalidades ....................................................................................... 100
8.4.2 - Resultados da análise estática linear .................................................... 101
8.4.3 - Deslocamentos nodais .......................................................................... 101
8.4.4 - Esforços nos elementos ........................................................................ 101
8.4.5 - Reações dos apoios .............................................................................. 102
xvi
9 – RESULTADOS .................................................................................................. 103
10 - CONCLUSÃO .................................................................................................. 110
11 – REFERÊNCIA BIBLIOGRÁFICA ..................................................................... 111
17
1 – INTRODUÇÃO
A engenharia estrutural trata basicamente do planejamento, projeto,
construção e manutenção de sistemas estruturais para transporte, moradia, trabalho
e lazer. Sendo que o projeto e a execução de estruturas sejam elas de concreto,
madeira ou aço são subáreas de conhecimento da engenharia civil onde
engenheiros se especializam, sendo assim chamados engenheiros estruturais.
A análise estrutural é a fase do projeto estrutural em que é feita a idealização
do comportamento da estrutura. Esse comportamento pode ser expresso por
diversos parâmetros, tais como pelos campos de tensões, deformações e
deslocamentos na estrutura. De uma maneira geral, a análise estrutural tem como
objetivo a determinação de esforços internos e externos (forças e reações de apoio),
e consequentemente a obtenção de tensões, deformações e os correspondentes
deslocamentos da estrutura que está sendo projetada. Essa análise deve ser feita
para os possíveis estágios de carregamentos e solicitações que devem ser
previamente determinados. (MARTHA, LUIZ FERNANDO, 2010, p. 1)
Os efeitos da constituição interna molecular dos materiais são levados em
conta de forma macroscópica através das equações constitutivas dos materiais, com
base na lei de Hooke, onde considera-se o material solicitado dentro de limites que
garantem seu comportamento elástico linear. (RIBEIRO, F. L. B. , 2004, p. 4).
A primeira etapa de todo processo de modelagem computacional de um
fenômeno físico consiste da identificação dos fatores que podem influenciar de
maneira relevante no problema. Isto implica na escolha adequada dos princípios
físicos e das variáveis dependentes e independentes que descrevem o problema,
resultando em um modelo matemático constituído por um conjunto de equações
diferenciais que geralmente são de difícil solução, portanto a segunda etapa que
consiste em obter a solução do modelo matemático, deve ser atribuída aos métodos
numéricos, de modo a simplificar de forma altamente satisfatória a solução do
problema. (RIBEIRO, F. L. B. 2004, p. 4).
Inúmeros métodos de precisão para solução destes problemas são usados
em engenharia entre eles pode-se destacar: método dos elementos de contorno,
método das diferenças finitas, método dos volumes finitos, método de Galerkin,
método de Rayleigh-Ritz e o método dos elementos finitos. (SILVA, S. Introdução ao
Método dos Elementos Finitos, 2009, P. 10)
18
O Método dos Elementos Finitos - MEF que será abordado neste trabalho foi
idealizado com os trabalhos de Argyris e Kelsey (1954, apud RODRIGUES 1997,
p.1) e de Turner et al (1956, apud RODRIGUES 1997, p.1). (Pedreiro, M. R. M.
2011). Com isso, os pesquisadores passaram a ter uma ferramenta poderosa que
permite a modelagem numérica dos fenômenos envolvidos na análise estrutural.
A ideia básica do MEF é realizar uma divisão do domínio de integração de
uma estrutura ou sistema de interesse em um conjunto de pequenas regiões,
chamadas de elementos finitos transformando o domínio de contínuo para discreto.
Esta divisão do domínio é conhecida como malha ou grid, que nada mais é do que o
conjunto de elementos finitos resultante da discretização.
Figura 1 - Modelo de estrutura contínua discretizada pelo método dos elementos finitos.
Fonte: SILVA (2009)
A malha é formada de elementos compostos de faces e nós, que são pontos
de intersecção e ligação entre os elementos. O grande mérito do MEF é não buscar
uma função admissível que satisfaça as condições de contorno para todo o domínio,
o que pode ser praticamente impossível em um problema complexo, e sim buscar
estas soluções em cada elemento de forma separada. (SILVA, S. Introdução ao
Método dos Elementos Finitos, 2009, P. 10 e 11)
No âmbito da Engenharia de Estruturas, o Método dos Elementos Finitos
(MEF) tem como objetivo a determinação do estado de tensão e de deformação de
um sólido de geometria arbitrária sujeito a ações exteriores. Este tipo de cálculo tem
19
a designação genérica de análise de estruturas e surge, por exemplo, no estudo de
edifícios, pontes, barragens, etc. Quando existe a necessidade de projetar uma
estrutura, é habitual proceder-se a uma sucessão de análises e modificações das
suas características, com o objetivo de se alcançar uma solução satisfatória, quer
em termos econômicos, quer na verificação dos pré-requisitos funcionais e
regulamentares. (AZEVEDO, Álvaro F.M. 2003, P. 1)
2 - OBJETIVO
O objetivo deste projeto é apresentar de forma introdutória os aspectos mais
relevantes do método dos elementos finitos na solução de sistemas estruturais
composto por elementos lineares que constituem pórticos bidimensionais, utilizando-
se a linguagem de programação em Visual Basic no auxílio para elaboração de um
código computacional, que possibilite a solução destes sistemas com a
determinação de deslocamentos, tensões e deformações em cada elemento, além
das reações provenientes de vinculações externas.
A abordagem do MEF envolve conceitos elementares da teoria de funções,
álgebra e cálculo, os quais são abordados nas disciplinas básicas geralmente
desenvolvidas nos primeiros semestre dos cursos de engenharia.
Ao longo de todo este trabalho consideram-se as seguintes hipóteses:
- Linearidade física
- Linearidade geométrica
- Homogeneidade e isotropia do material estrutural
Será adotado como hipótese simplificadora o fato de que o material apresenta
linearidade física que permite assumir um comportamento elástico linear. Este fato
simplifica as relações constitutivas, permitindo o estabelecimento de uma relação
linear entre esforços e deformações. Além disso, também será assumido que a
estrutura apresenta linearidade geométrica que inclui a hipótese dos pequenos
deslocamentos e das pequenas deformações, tal hipótese permite que as condições
de equilíbrio possam ser estabelecidas com base na configuração indeformada da
estrutura.
20
3 – METODOLOGIA
Foi abordado a apresentação da conceituação de métodos dos elementos
finitos com a utilização de exemplos de elementos estruturais simples (elementos
estruturais reticulares), mesmo sabendo que para eles também é possível deduzir,
de modo direto (e clássico) as equações diferenciais regentes que tenham solução
analítica fechada; nesses casos, como é sabido, também poderiam ser obtidos
resultados imediatos de valores de deslocamentos de pontos particulares,
alternativamente, com o emprego dos princípios dos trabalhos virtuais ou com o
teorema de Castigliano, por exemplo.
21
4 – ELEMENTOS ESTRUTURAIS RETICULARES – VIGA
PRISMÁTICA
Inicia-se com a revisão dos aspectos principais do problema de viga
prismática, no âmbito da teoria de primeira ordem, isto é, quando ocorrerem
pequenos deslocamentos angulares da estrutura e pequenas deformações
específicas no material elástico.
A deformação relativa entre duas seções transversais separadas por um dx,
admitida a hipótese da manutenção das seções planas, pode ser obtida conforme a
seguir mostrado na figura 2.
Figura 2 - Viga elemento infinitesimal.
Sendo flexão normal simples:
/ .My I (1)
Com hooke:
.
My
E EI
(2)
Por semelhança de triângulos:
.
u Y
x
(3)
Então:
.
l M
EI
(4)
22
Além disso:
/ .d dx (5)
Então:
.
Md dx
EI
(6)
5 - FORMULAÇÃO LOCAL
5.1 - Dedução direta da equação diferencial regente do problema de viga
prismática
De acordo com a Teoria clássica da Resistência dos materiais, o caso de viga
prismática sob força distribuída p(x), estudado com base em elemento de
comprimento infinitesimal e com a conhecida relação entre esforços (M) e
deformações (curvaturas k) tem a formulação local seguinte:
Figura 3 - Sistema de eixos da Estática
(LABAKI, J; MESQUITA, E.)
Figura 3b: Convenção de sinais da Resistência
dos Materiais (LABAKI, J; MESQUITA, E.)
Do equilíbrio de forças e de momentos, de elemento infinitesimal, resulta:
,dV
pdx
.dM
Vdx
23
A relação entre esforços e deformações, já obtida anteriormente, será:
l M
p EI
(7)
Para compatibilizar convenções usuais sobre esforços solicitantes, e admitida
a possibilidade de aproximar a curvatura (l/p) com a derivada segunda de v,
resultará:
²
²
d v M
dx EI
(8)
A diferenciação sucessiva desta expressão e a utilização das equações de
equilíbrio levarão à equação diferencial do problema, dada por:
4
4
d v p
dx EI
(9)
6 – A EVOLUÇÃO DO MÉTODO DOS DESLOCAMENTOS
O método dos Elementos Finitos pertence à família do Método dos
Deslocamentos ou Método da Rigidez onde deslocamentos são escolhidos como
incógnitas. Todos os membros dessa família se caracterizam por ter como equação
fundamental a equação de equilíbrio cujas incógnitas são deslocamentos
generalizados. Entendem-se aqui por deslocamentos generalizados, grandezas
cinemáticas, tais como, deslocamentos lineares, rotações etc.
Os membros dessa família formam uma árvore genealógica, com novos
métodos gerados a partir dos métodos mais antigos. De certa maneira, a evolução
do método ao longo do tempo segue as leis da evolução de Darwin, com mutação e
seleção. Os novos membros da família desses métodos herdam as características
de seus antecessores, mas sofrem pequenas mudanças que só são bem sucedidas
se forem bem adaptadas ás condições existentes. Um exemplo disso é que a
Análise Matricial de Estruturas (AME) e o MEF só tiveram larga aceitação quando os
computadores atingiram uma fase de elevado grau de desenvolvimento, apesar de
este último ter surgido antes dessa fase.
24
Este capítulo procura mostrar como se deu a evolução do Método dos
Deslocamentos, desde as primeiras formulações até o MEF. É surpreendente
verificar como as mudanças conceituais são pequenas em comparação ao enorme
crescimento do potencial do método. (VAZ, L. E. 2010)
6.1 – Método Básico
A análise de estruturas usa três equações básicas, nomeadamente equações
de compatibilidade, de equilíbrio e constitutivas, também chamadas de relação
tensão-deformação. O método dos deslocamentos caracteriza-se por usar a
equação de equilíbrio como equação fundamental, ou seja, aquela de onde são
obtidas as incógnitas primárias do problema, a partir das quais, todas as outras
respostas serão obtidas. As incógnitas primárias são os deslocamentos por meio
dos quais é possível obter deformações, tensões, resultantes de tensões etc.
O método básico da família do método dos deslocamentos consiste em
manipular as três equações básicas da análise de estruturas de modo a colocar
todas as informações disponíveis nas equações de equilíbrio com deslocamentos
livres como incógnitas. O número de deslocamentos livres é chamado grau de
liberdade da estrutura.
Neste item e em outros que seguem, a estrutura apresentada na figura 4 é
utilizada para ilustrar a resolução do método. Trata-se de uma treliça plana simples
com quatro barras e dois graus de liberdade, os deslocamentos horizontal e vertical
do nó C. (VAZ, L. E. 2010)
Figura 4 – Treliça com 2 graus de liberdade.
25
As equações de compatibilidade relacionam grandezas cinemáticas, nesse
caso os deslocamentos nodais livres 1d e 2d na direção horizontal e vertical com
alongamentos/encurtamentos i das barras i . Os deslocamentos são supostos
positivos com os sentidos indicados na figura acima. Os alongamentos serão
considerados positivos e os encurtamentos negativos. As expressões para os i das
quatro barras são obtidas projetando-se os deslocamentos nodais nas direções das
barras, assim:
1 1 2 1 2
2 1 2 1
3 1 2 1 2
4 1 2 1 2
2 2( , ) d
2 2
( , )
;2 2( , )
2 2
2 2( , )
2 2
d d d
d d d
d d d d
d d d d
(10)
A segunda equação de compatibilidade relaciona os
alongamentos/encurtamentos das barras i com as deformações longitudinais i .
Da resistência dos materiais:
i
iL
;
(11)
Como os comprimentos das barras são:
1
2
3
4
2
2
2
L L
L L
L L
L L
(12)
26
Chega-se a:
1 2
1 1 2 1 2
12 1 2
1 2
3 1 2 1 2
1 2
4 1 2 1 2
2 2d
12 2( , )22
( , )
;2 212 2( , )
22
2 212 2( , )
22
d
d d d dLL
dd d
L
d d
d d d dLL
d d
d d d dLL
(13)
Para efeito de simplificação, a lei constitutiva usada nesse trabalho será a lei
de Hooke, assim, para cada barra, i vale:
;i iE (14)
Ou, e termos de esforços normais Ni ,
N ;i i
i
EA L
(15)
Onde E é o módulo de elasticidade do material, A, a área da seção
transversal (as duas grandezas supostas constantes para todas as barras), Ni o
esforço normal e iL o comprimento da barra i .
Substituindo-se para cada barra i , i dado em (10) em (15), obtém-se:
27
1 1 2 1 2
2 1 2
3 1 2 1 2
4 1 2 1
1
2
( , )2
( , )
;
( , )2
)
( ,
2
EAN d d d d
L
E AN d d
L
EN d
d
A
L
A
d d d
EN d d d d
L
(16)
As equações de equilíbrio são obtidas para as direções horizontal e vertical
no nó C.
Os sentidos das forças axiais Ni que atuam nas barras i , são admitidos a princípio
de tração. Para se escrever as equações de equilíbrio, valem, no entanto os
sentidos indicados na Figura 5.
Figura 5 – Equilíbrio do nó C.
As equações de equilíbrio são:
Na direção horizontal:
1 2 3 4
2 2 20 ; 0 ;
2 2 2hF N N N N P
(17)
Na direção vertical:
1 3 4
2 2 20 ; 0 ;
2 2 2vF N N N
(18)
Substituindo-se as expressões (16) em (18) e manipulando-as, obtém-se:
28
1 2
1 2
2,061 A 0,354
0,354 E 1,061 0
E E Ad d P
L L
A E Ad d
L L
(19)
A expressão (19) é a equação fundamental do método dos deslocamentos
para a análise da treliça plana da Figura 4. Matricialmente, ela pode ser reescrita
como:
1
2
2,061 0,354 ;
0,354 1,061 0
d PE A
dL
(20)
Cuja solução é:
1
2
0,515 ;
0,171
d P L
d E A
(21)
Com os deslocamentos 1d e 2d é possível obter agora todas as respostas da
estrutura em termos de alongamento/encurtamento, na expressão (10), deformações
em (13), tensões em (14), e esforços normais Ni em (16). Tais valores estão
indicados a seguir:
1
2
3
4
0,243
0,515 ;
1,778
0,243
P L
E A
(22)
1
2
3
4
0,172
0,515 ;
0,343
0,172
P
E A
(23)
29
1
2
3
4
0,172
0,515 ;
0,343
0,172
P
A
(24)
1
2
3
4
0,172
0,515;
0,343
0,172
N
NP
N
N
(25)
6.2 – Método clássico
O método clássico é essencialmente o mesmo que o método básico. Sua
contribuição foi no sentido de sistematizar, ou seja, organizar, ou ainda criar uma
metodologia que possa ser aplicada da mesma forma a todas as estruturas.
O método usa os conceitos de estados auxiliares e de superposição de
efeitos. Inicialmente, devem-se identificar os graus de liberdade da estrutura. Em
seguida, um estado auxiliar j é criado para cada grau de liberdade impondo-se um
valor unitário para o grau de liberdade jd , enquanto os outros são mantidos nulos.
Resultantes das forças internas resistentes que atuam nas barras aparecem nas
direções dos graus de liberdade. A força interna na direção i devido ao
deslocamento unitário na direção na direção do grau de liberdade jd é chamada de
coeficiente ijk . Além disso um estado auxiliar 0 é criado para as cargas atuantes
com todos os graus de liberdade mantidos fixos. As forças atuantes que atuam nos
nós na direção do grau de liberdade jd nesse estado são denominadas cargas
nodais jf .
Como os estados auxiliares não são auto equilibrados o equilíbrio é
conseguido com superposição de efeitos. Assim, somando-se os produtos das
forças internas resultantes (nas direções dos graus de liberdade) correspondentes a
cada estado auxiliar j por jd , a soma deve ser igual às forças aplicadas (nas
direções dos graus de liberdade) no estado auxiliar 0. Em termos físicos, isso
significa que os deslocamentos que surgem na direção dos graus de liberdade jd
devem ser tais que as forças internas equilibrem as forças aplicadas. (VAZ, L. E.
2010)
30
A aplicação das ideias descritas no exemplo 6.1 ajuda a esclarecer o método.
Estado auxiliar 1, 1d = 1.
Figura 6 – Termos 11k e 21k da matriz de rigidez da treliça
Estado auxiliar 2, 2d =1.
Figura 7 – Termos 12k e 22k da matriz de rigidez da treliça
Para se obter os coeficientes ijk (força interna resultante na direção i devida
a um deslocamento unitário na direção j ) procede-se da seguinte maneira :
inicialmente, calculam-se os alongamentos/encurtamentos das barras ijd
(alongamento/encurtamento na barra i devido a uma deslocamento unitário na
direção do grau de liberdade jd ) de forma análoga ao que foi feito para se obter os
alongamentos/encurtamentos em (10).
31
Para o estado auxiliar 1.
11
21
31
41
2
2
1
2
2
2
2
(26)
Para o estado auxiliar 2.
12
22
32
42
2
2
0
2
2
2
2
(27)
Utilizando-se a relação constitutiva é possível calcular os esforços normais
nas barras Nij (esforço normal na barra i devido a uma deslocamento unitário na
direção do grau de liberdade jd ) com uma expressão análoga a (15).
;
ij
ij
i
N E AL
(28)
32
Assim:
Para o estado auxiliar 1.
11
21
31
41
2
;
2
2
EAN
L
E AN
L
EN
EN
A
L
A
L
(29)
Para o estado auxiliar 2.
12
22
32
42
2
2
2
0
;A
EAN
L
N
EN
NAE
L
L
(30)
Os coeficientes de rigidez ijk (esforço na direção i para um deslocamento
unitário na direção j ) são calculados utilizando-se as equações de equilíbrio no nó
C. Assim, das equações de equilíbrio na direção horizontal e vertical da Figura 8, da
correspondente a 1 1d obtém-se, respectivamente, os coeficientes 11k e 21k .
Figura 8 – Forças no nó C para 1 1d e 2 1d
33
Para o estado auxiliar 1, Figura 8.a.
11
21
2,061
0,354
Ek
L
Ek
L
A
A
(31)
Para o estado auxiliar 2, Figura 8.b.
12
22
0,354
,
1 061
Ek
L
kL
A
AE
(32)
O estado auxiliar 0, fornece:
1
2 0
f P
f
(33)
A superposição de efeitos, que deve garantir o equilíbrio das forças
resistentes e aplicadas, pode agora ser escrita como:
11 12 1 1
21 22 2 2
;k k d f
k k d f
(34)
ou com os valores da estrutura sendo analisada:
1 1
2 2
2,061 0,354 0,515 ;
0,354 1,061 0 0,171
d dPE A P L
d dL E A
(35)
A expressão (35) é idêntica à expressão (20), como não poderia deixar de
ser. Desse modo, as respostas das estruturas obtidas pelo método básico dadas
pelas expressões de (21) a (25) serão as mesmas.
34
6.3 - Método da análise matricial
A análise matricial de estruturas reticuladas sistematizou as operações
matemáticas da análise de estruturas fazendo uso da álgebra matricial que opera
com vetores e matrizes. Ela introduziu diversos conceitos novos na análise de
estruturas. Toda a sistematização se baseia na ideia de sistema local e sistema
global de coordenadas. Com esse conceito definido, é possível estabelecer matrizes
de rigidez de elemento nos sistemas local e global, assim como vetores de forças
nodais de elemento nos sistemas local e global. A partir das contribuições das
matrizes de rigidez e dos vetores de forças nodais de elemento no sistema global,
pode-se montar a matriz de rigidez bem como o vetor de forças nodais da estrutura.
Deslocamentos nodais também são definidos nos sistemas local e global. Uma
equação de equilíbrio da estrutura no sistema global fornece os deslocamentos
nodais. Uma vez obtidos os deslocamentos nodais da estrutura, as forças atuantes
nas extremidades dos elementos podem ser determinadas. (VAZ, L. E. 2010)
O sistema local de coordenada é definido quando se escolhe os nós inicial e
final do elemento. Na figura 9, os nós 1 e 2 são, respectivamente, o nó inicial e o nó
final do elemento ou barra.
A estrutura de treliça plana tratada até aqui tem dois graus de liberdade por
nó. Ao nó 1 são associados os deslocamentos 1 e 2 e ao nó 2, os deslocamentos 3
e 4. A figura 9 indica os sentidos positivos dos 4 componentes do vetor de
deslocamentos ld , no sistema local, e gd , no sistema global. O ângulo define a
rotação do eixo da barra em relação ao sistema global. Associados aos vetores de
deslocamentos, são criados também os vetores de forças nodais lf , no sistema
local, e gf , no sistema global.
35
cosc
s sen
Figura 9 – Graus de liberdade no sistema global e local
Os vetores dos deslocamentos de elemento no sistema local ld e global gd podem
ser relacionados pela matriz de rotação R, como indicado a seguir:
11
22
33
44
0 0
0 0
0 0
0 0
gl
gl
gl
gl
dd c s
dd s c
dd c s
dd s c
(36)
Ou, sucintamente:
l gd Rd (37)
Como o trabalho é um escalar independente do sistema de coordenadas, ele
deve ser o mesmo nos sistemas local e global.
g lW W (38)
t t
g g l ld f d f (39)
Substituindo (37) em (39), obtém-se:
t
t t t
g g g l g ld f Rd f d R f (40)
36
t
g lf R f (41)
As expressões (37) e (41) formam o princípio da contragradiência que pode
ser enunciado como: “Se uma matriz transforma deslocamentos globais em locais,
sua transposta transforma forças locais em globais”. (VAZ, L. E. 2010)
A matriz de rigidez do elemento de treliça plana no sistema local para o
elemento m, ,l mK é dada em (42). Ela é obtida da definição dos coeficientes de
rigidez , ( )l m ijK . O coeficiente , ( )l m ijK significa a força na direção do deslocamento
local i para um deslocamento unitário aplicado na direção do deslocamento local j ,
mantendo os outros deslocamentos locais nulos.
1 0 1 0
0 0 0 0
1 0 1 0
0 0 0 0
m mlm
m
E Ak
L
(42)
Onde mE é o módulo de elasticidade do matéria, mA a área da seção
transversal e mL o comprimento da barra m . A equação de equilíbrio da barra que
relaciona deslocamentos, forças e a matriz de rigidez no sistema local de
coordenadas é dada por:
1 1
2 2
3 3
4 4
1 0 1 0
0 0 0 0
1 0 1 0
0 0 0 0
m m
m m
m m
m m
l l
l lm m
l lm
l l
d f
d fE A
d fL
d f
(43)
Ou sucintamente:
dm mm l llK f (44)
37
A matriz de rigidez do elemento m no sistema global de coordenadas mgK
pode ser obtida como explicado a seguir. Substituindo-se (37) em (44), obtém-se:
R dmm ml m g lK f (45)
Multiplicando-se ambos os lados de (45) por R t
m , chega-se a:
R R d R m m m
t t
m m g m llK f (46)
Usando (41), obtém-se:
m m mg g gK d f (47)
Onde,
R Rmm
t
g m mlKK (48)
A partir da matriz de rigidez e das forças nodais de cada elemento k no
sistema global é feita então a montagem da matriz de rigidez K e das forças nodais
f globais da estrutura em função da conexão entre os elementos (incidência),
obtendo-se a equação de equilíbrio global da estrutura. (VAZ, L. E. 2010)
K d f (49)
Sendo d os deslocamentos da estrutura no sistema global de coordenadas.
Uma vez obtido d , é possível calcular os deslocamentos nodais de cada
elemento no sistema global mgd e girar esses deslocamentos para o sistema local
mld via (37) e calcular as forças de extremidade finais em cada elemento no sistema
local ml
f via (44). (VAZ, L. E. 2010)
38
6.4 – Método de Castigliano
O método de Castigliano é assim chamado em homenagem ao segundo
teorema de Carlo Alberto Castigliano, que, em 1973, demonstrou que a derivada da
energia de deformação de uma estrutura em relação ao deslocamento id é igual a
força externa da estrutura na mesma direção. A demonstração foi feita estruturas
com comportamento linear elástico, mas ela é válida também para materiais
elásticos não lineares. Nesse item, a demonstração será entendida a estruturas
material elástico não linear.
Esse teorema representou um importante passo no desenvolvimento da
análise de estruturas porque ele mostrou um novo caminho, baseado em teoremas
de energia, para se formular um método para análise de estruturas. Esse caminho
levou ao MEF. (VAZ, L. E. 2010)
6.4.1 – Energia de deformação
Para efeito de simplificação, a apresentação do Segundo Teorema de
Castigliano será feita aqui para o caso de uma estrutura de treliça. Nesse tipo de
estrutura, somente um componente de deformação e de tensão atua no elemento de
barra, nomeadamente, a deformação e a tensão normal longitudinal, ou seja, trata-
se de um problema unidimensional para efeito da relação tensão x deformação. Seja
a relação tensão x deformação apresentada na figura 10. A solicitação externa levou
a tensão atuante até o valor final m que corresponde à deformação final m na
barra m da treliça. (VAZ, L. E. 2010)
Figura 10 – Energia de deformação específica 0U da barra m.
39
A energia de deformação específica moU da barra m é definida como:
0
0
m
m m m mU d
(50)
O adjetivo “específica” deve-se ao fato de 0mU ser, em termos de unidades,
um trabalho por unidade de volume.
A energia de deformação da barra m , mU , é obtida integrando-se no volume
da barra.
0 m
m
m m m m
V
U U dV (51)
Para se obter a energia de deformação U relativa a toda a treliça, somam-se
os mU de todas as barras, de 1 a nb , onde nb é o número de barras da estrutura.
1 2 1, ,...,
nb
m m mmU U
(52)
Onde m é a deformação final da barra m . Como a deformação final da barra,
m depende do alongamento/encurtamento longitudinal final da barra m , como
expresso em (11), que por sua vez, m depende dos deslocamentos nodais finais
das extremidades da barra no sistema global de coordenadas id como
exemplificado em (13), a expressão (52) pode ser reescrita como:
1 2 1, d ,..., d
nb
n m mmU d U
(53)
Onde n é o número de graus de liberdade da estrutura de treliça.
A energia de deformação da estrutura corresponde fisicamente à energia
armazenada na estrutura quando ela se deforma, caso não haja perda de energia,
40
ou seja, para um sistema conservativo. Essa energia é responsável pela volta da
estrutura a sua configuração inicial, antes da aplicação das cargas, quando estas
são retiradas da estrutura. (VAZ, L. E., 2010)
6.4.2 – Trabalho Externo
O trabalho externo total W em uma estrutura de treliça plana pode ser obtido
somando-se os trabalhos externos iW referentes aos graus de liberdade i da
estrutura.
1 2 1 1
0
, d ,..., d ; ;id
n n
n i i i ii iW d W f u du
(54)
Onde n , como anteriormente, é o número de graus de liberdade da estrutura.
A figura 11 esclarece.
Figura 11 – Trabalho externo associado ao grau de liberdade i .
6.4.3 – Segundo teorema de Castigliano
Substituindo doravante a notação do deslocamento final d por d para efeito
de simplificação, a energia de deformação (53) e o trabalho externo (54) em uma
estrutura de treliça plana, como visto nos itens 6.4.1 e 6.4.2, podem ser escritos
como uma função do vetor dos deslocamentos nodais finais da estrutura no sistema
global de coordenadas d com n componentes.
41
Expandindo-se ( )W d em série de Taylor até o termo de primeira ordem, é
possível expressar o incremento de ( )W d como:
tW d
W d d W d dd
(55)
tW d
W d W d d W d dd
(56)
Procedendo-se da mesma maneira para U d , obtém-se:
tU d
U d d U d dd
(57)
t
U dU d U d d U d d
d
(58)
Pelo princípio da conservação de energia em sistemas conservatórios, todo
trabalho externo realizado é armazenado na estrutura em termos de energia de
deformação. (VAZ, L. E. 2010) Assim, o incremento de trabalho externo é igual ao
incremento de energia de deformação, logo:
W d U d (59)
Ou seja,
t t
W d U dd d
d d
(60)
Ou, ainda, para uma variação arbitrária d ,
i i
W d U d
d d
(61)
42
O teorema da integral de Newton diz que:
0
a
f a f x dxa
(62)
Logo, utilizando-se esse teorema, pode-se escrever:
0
id
i i
i i i ii dW d
f ud d
u f d f
(63)
Onde, como foi redefinido no início desse item, id em (63) é o valor final da
variável deslocamentos nodal iu e if é a força final associada ao deslocamento id .
Com o uso de (61) e (63), obtém-se finalmente a expressão do Segundo
Teorema de Castigliano:
i
i
U df
d
(64)
Ou, grupando-se todas as equações (66) correspondentes aos n graus de
liberdade em uma só equação:
U df
d
(65)
Observa-se que o termo à esquerda da expressão (65) corresponde ao vetor
das forças internas resistentes, doravante denominado rf d , e o termo à direita,
corresponde ao vetor das forças solicitantes, doravante denominado sf .
r sf d f (66)
43
A expressão (66) fornece um método de análise de estruturas denominado
método de Castigliano. A expressão fornece n equações que permitem obter as n
incógnitas do problema, ou seja, os n deslocamentos nodais 1,..., ,i nd i . Se a
estrutura tiver um comportamento linear, as equações (66) fornecem um sistema de
n equações algébricas lineares, caso o comportamento seja não linear, n equações
não lineares são obtidas. O sistema de n equações não lineares pode ser resolvido,
por exemplo, pelo método de Newton-Raphson para se obter as n incógnitas do
problema, ou seja, os n deslocamentos nodais 1,..., ,i nd i .
A aplicação do método na análise da treliça plana da figura 4 ajuda a
esclarecer as expressões descritas anteriormente.
6.4.4 – A aplicação do método de Castigliano
A lei de hooke para materiais linear-elásticos permite escrever:
E (67)
A energia de deformação específica 0U pode ser escrita em função da
deformação final da barra m . Empregando-se novamente a notação m para
representar o valor final da grandeza m , chega-se a:
2
00 0
2
m m mm m m mU d E d E
(68)
A energia de deformação mU para a barra m vale:
2 22
00
1
2 2 2
m
m
lm m
m k m
mV A
EAU E dV E dA dx
L
(69)
Usando as equações de compatibilidade para a treliça da Figura 4 descritas
em (10) e abandonando mais uma vez, para efeito de simplificação, o sobrescrito –
para representar valores finais das variáveis, obtém-se:
44
1 1 2 1 2
2 1 2 1
3 1 2 1 2
4 1 2 1 2
2 2( , ) d
2 2
( , )
;2 2( , )
2 2
2 2( , )
2 2
d d d
d d d
d d d d
d d d d
(70)
E as expressões dos comprimentos das barras dadas em (12), podem-se
escrever:
2
2
1 1 2 1
1
21 1 2
1 1,
2 2d
22 2,
2 2
EA EAU d d d d
L Ld
(71)
2 2
2 1 2 2 1 2 1
2
1 1, ,
2 2
EA EAU d d d d d
L L
(72)
2
2
3 1 2 1
3
23 1 2
1 1,
2 2d
22 2,
2 2
EA EAU d d d d
L Ld
(73)
1 2
2
2
4 1 2 4 1 2
4
1 1, ,
2 2
2 2d
22 2
EA EAU d d d d
Ld
L
(74)
Usando-se (53) para se obter a energia de deformação total da estrutura,
obtém-se:
2 2
1 2 1 2 1
2
2
1 2 21
2 2 2 2 2 21, 2 d d d
2 2 2 2 2 2( )
2 2
EAU d d dd d
Ld
(75)
Aplicando-se agora a expressão (64) do Segundo Teorema de Castigliano,
obtém-se:
1 2
1
3 2 21
4 4
U d EAd d P
d L
(76)
45
1 2
2
2 3 20
4 4
U d EAd d
d L
(77)
Ou, ainda,
1 1
2 2
2,061 0,354 0,515 ;
0,354 1,061 0 0,171
d dPE A P L
d dL E A
(78)
Que é idêntica a (35).
6.5 – Princípio dos deslocamentos virtuais
6.5.1 – Incrementos da energia de deformação
O princípio dos trabalhos virtuais será demonstrado neste item para estruturas
de treliça. Uma barra de treliça m é carregada até que a deformação final m seja
atingida como indicado na Figura 12. A tensão atuante correspondente é ( )m m . A
energia de deformação específica produzida na barra é 0mU . Imagine agora que um
incremento de tensão m seja aplicado à barra a partir de m . Um incremento de
deformação m correspondente ocorre na barra. (VAZ, L. E. 2010)
Figura 12 – Incremento de energia de deformação específica 0,mU da barra m .
46
O incremento total da energia de deformação específica 0mU correspondente
à aplicação m pode ser escrito como:
0 0
1 U
2m mm m m m m mU erro (79)
Ou,
1 2
0 0 0 0Um m m m mU U U erro (80)
Onde
1
0m m m mU (81)
2
0
1
2m m mU (82)
Os termos 1
0mU e
2
0mU são denominados incremento de primeira e de
segunda ordem de 0mU , respectivamente. O termo de primeira ordem corresponde à
área do retângulo vertical hachurado representado na figura 12. O termo de segunda
ordem corresponde à área do triângulo maior na mesma figura. A área em cinza
corresponde ao erro cometido no cálculo do incremento total 0Um
erro . (VAZ, L. E.
2010)
Como a energia de deformação da barra m da treliça mU é obtida pela
integração no volume da barra da energia de deformação específica, obtém-se:
0
0 m
Vm
m mU U dV (83)
Logo,
0
0 0 0
1
2 m
Vm Vm Vk
m m m m m m m m m mU dV dV erroU dV (84)
47
Ou
1 2 Um m m m mU U U erro (85)
Onde
1
0
Vm
m m m m mU dV (86)
2
0
1
2
Vm
m m m mU dV (87)
A energia de deformação de toda estrutura com m barras pode ser obtida
somando-se a energia de deformação de todas as barras, assim:
1
nb
mmU U
(88)
Logo,
1 2
1 1 1U
nb nb nb
m m m mm m mU U U erro
(89)
Ou
1 2 UU U U erro (90)
Onde
1
1 0
Vmnb
m m m mmU dV
(91)
2
1 0
1
2
Vmnb
m m mmU dV
(92)
48
As expressões (91) e (92) pode ser generalizadas para o caso em que há
várias componentes de tensão, por exemplo, x , y e xy , e de deformação, por
exemplo x , y e xy atuando em um elemento infinitesimal de elemento m da
estrutura com n elementos.
Nesse caso pode-se escrever:
1
1 0
Vmne t
m m mmU dV
(93)
2
1 0
1
2
Vmne t
m m mmU dV
(94)
Onde m , m e m representam, respectivamente, os vetores das
componentes de tensão atuantes, dos incrementos das componentes de tensão
atuantes e dos incrementos das componentes de deformação no elemento m .
6.5.2 – Incrementos do trabalho externo
Os incrementos do trabalho externo podem ser obtidos pelo raciocínio
análogo ao desenvolvido no item anterior para a energia de deformação.
Uma força é aplicada em um dado grau de liberdade i até produzir um
deslocamento final id como representado na figura 13. A força atuante
correspondente à id é if . O trabalho externo produzido correspondente ao grau de
liberdade i é iW . Imagine agora que um incremento de força if é aplicado à força
if . Um incremento de deslocamento id ocorre no grau de liberdade
correspondente. (VAZ, L. E. 2010)
49
Figura 13 – Incremento de trabalho externo iW .
O incremento total do trabalho externo iW correspondente à aplicação de
if no grau de liberdade i pode ser escrito como:
1
2i i i i i i iW f d f d erroW d
(95)
Ou
1 2
i i i i iW W W erroW d (96)
Onde
1
i i iW f d (97)
2 1
2i i iW f d
(98)
Os termos 1
iW e 2
iW são denominados respectivamente incremento de
primeira e de segunda ordem de iW . O termo de primeira ordem corresponde à área
do retângulo vertical hachurado representado na figura 13. O termo de segunda
50
ordem corresponde à área do triângulo maior na mesma figura. A área em cinza
correspondente ao erro ao erro cometido no cálculo do incremento total ierroW .
O trabalho externo correspondente a toda a treliça com n graus de liberdade
pode ser obtido somando-se o trabalho externo de todos os graus de liberdade,
assim:
1
n
iiW W
(99)
Logo,
1 2
1 1 1
n n n
i i i ii i iW W W erroW d
(100)
Ou, ainda,
1 2
iW W W erroW d (101)
Onde,
1
1
n
i iiW f d
(102)
2
1
1
2
n
i iiW f d
(103)
As expressões (102) e (103) podem ser escritas usando-se vetores:
1 tW f d (104)
2 1
2 t dW f
(105)
Onde f , d e f representam, respectivamente, os valores das forças
solicitantes nodais finais, dos incrementos dos deslocamentos nodais e dos
incrementos das forças nodais.
51
6.5.3 – Formulação do princípio dos deslocamentos virtuais
O princípio dos deslocamentos virtuais baseia-se no princípio de conservação
de energia. Seu enunciado é o seguinte: “Para toda estrutura, o incremento de
primeira ordem da energia de deformação é igual ao incremento de primeira ordem
do trabalho externo.” (VAZ, L. E. 2010). A aplicação do princípio não se limita a
sistemas conservativos. Matematicamente, ele pode ser expresso por:
1 1U W (106)
Para o caso geral em que há várias componentes de tensão e deformação
atuando em um elemento infinitesimal de um elemento m de uma estrutura com n
elementos, a expressão (106) pode ser escrita como:
01
ne Vm
t
m m
m
mdV f d
(107)
As grandezas m e d em (107) são cinemáticas, virtuais e compatíveis
enquanto que as grandezas m e f são ditas estáticas, reais e em equilíbrio. O
termo virtual é sinônimo de potencial, ou seja, pode vir a acontecer, não real. As
grandezas m e d estão relacionadas por equações de compatibilidade já que as
componentes de d produzem as componentes de m . As grandezas reais m e
f estão relacionadas por equações de equilíbrio já que as tensões reais m são
produzidas pelas forças reais f .
6.6 – Método da mínima energia potencial total
A energia potencial total d é definida para sistemas conservativos como:
52
pd U d W d (108)
Onde U d é a energia de deformação da estrutura, como definido em (51) e
(53), e pW d é o trabalho potencial das forças externas, dado por:
t
pW d f d (109)
Novamente, os sobrescritos -, utilizados para representar valores finais das
variáveis são retirados para efeito de simplificação. Em sistemas conservativos,
U d é a energia que traz a estrutura de volta à configuração inicial caso as forças
externas sejam retiradas da estrutura. pW d é o trabalho potencial, ou seja, aquele
que seria realizado caso a estrutura voltasse a sua configuração inicial e as cargas
permanecem atuando sobre ela. Assim, d é a energia total necessária para
trazer de volta a estrutura a sua configuração inicial com as cargas atuando sobre
ela. (VAZ, L. E. 2010)
6.6.1 - O princípio da mínima energia potencial total
O princípio da mínima energia potencial total enuncia que os deslocamentos
d de uma estrutura em equilíbrio estável tornam mínima a energia potencial total da
estrutura. Em outras palavras, uma estrutura que está em equilíbrio estável se
deformou de modo a gastar o mínimo de energia potencial total. (VAZ, L. E. 2010)
Matematicamente, a condição de primeira ordem de mínimo de uma função é dada
por:
0
d
d
(110)
Como já é conhecido o esquema usualmente empregado para a utilização do
princípio dos trabalhos virtuais no cálculo de deslocamentos em estruturas, será
53
tratado tal assunto abordando-o sob forma mais adequada para o encaminhamento
ao método dos elementos finitos, segundo uma de suas formulações.
Seja uma estrutura real que, por exemplo, pode ser uma viga em balanço sob
força concentrada de extremidade, como a Figura 14, deformada sob a ação real e
em equilíbrio.
Imagine-se que essa viga seja levada a uma outra posição deformada
próxima da real (por meio de uma ação virtual qualquer).
Para que se possa aplicar o princípio dos trabalhos virtuais, exige-se que a
nova posição exibida pela viga v v seja de estado de deslocamentos
compatíveis com os vínculos e que as forças externas dadas (no caso somente a
força P) e os consequentes esforços internos não variem.
Nessas circunstâncias, a força externa, que tinha um certo potencial (em
relação a um referencial que pode ser a própria posição elástica real) perde o
potencial *P f . (SAVASSI, W. 1996)
Como se trata de perda de potencial será indicado com a expressão:
* P f (111)
Do mesmo modo que a perda de potencial, correspondente à passagem da
posição indeformada para a posição deformada de equilíbrio que se pretende
determinar, deve ser indicada por:
*Pf (112)
Se o referencial fosse o nível do eixo da viga indeformada.
Correspondentemente, ocorrerão variações das deformações que,
multiplicadas pelos esforços internos reais, darão uma variação do trabalho interno.
Essa variação de trabalho interno é igual à variação U da energia de deformação
U.
Então, como *P f é o trabalho virtual das forças externas isto é:
* eP f W (113)
54
E a variação da energia de deformação mede o trabalho virtual interno:
iU W (114)
Tem-se:
U (115)
Isto também pode ser escrito como
0U (116)
ou
( ) 0U (117)
ou ainda,
0p (118)
onde p é definido como a energia potencial total.
A expressão 0p representa o princípio da mínima energia potencial total,
porque pode ser demonstrado que a solução obtida, a partir da imposição da
validade de 0p , é correspondente a uma situação de mínimo valor para p .
Essa demonstração deve incluir a prova de que 2 0p , isto é, que a segunda
variação de p é sempre positiva (configurações de equilíbrio estável). (SAVASSI,
W. 1996)
Resumindo, mostrou-se que se a estrutura está em equilíbrio (estável),
conforme:
0p e (2 0p ) (119)
55
ou U (120)
Pode ser mostrado, de modo inverso, que se 0p a estrutura está em
equilíbrio. Isso significará que da suposição U (ou 0p ) resultariam, por
demonstração, analiticamente, as expressões das equações diferenciais de
equilíbrio, e equações relativas a equilíbrios no contorno. A realização dessa
demonstração (considerando a expressão analítica de p , onde a função v x será
incógnita, tanto quanto à forma quanto à amplitude, requeria o emprego de
conhecimentos de Cálculo das Variações, que é um tema não familiar ao estudante
do curso de graduação em engenharia.
Todavia, desde já salienta-se que, na maioria dos casos, não há interesse na
obtenção das equações diferenciais de equilíbrio, por meio da utilização do caminho
que tem a imposição de 0p como ponto de partida. Para as equações
diferenciais obtidas, como acontece com a maioria dos casos de estruturas bi e
tridimensionais, poderia não haver solução analítica fechada.
Mas, no âmbito das soluções aproximadas e numéricas é a utilização de
princípio da mínima energia potencial total que permitirá resolver, com sucesso
muitos problemas. (SAVASSI, W. 1996)
6.7 - Método de Rayleigh-Ritz
O método de Rayleigh-Ritz representou um grande passo na evolução do
método dos deslocamentos, pois contribuiu decisivamente para o aparecimento do
MEF. O método de Rayleigh-Ritz é, na essência, o método do princípio da mínima
energia potencial total, mas, a pequena modificação introduzida nesse último
permitiu um grande avanço. Para uma melhor compreensão do método, o exemplo
da treliça usado até aqui vai ser substituído por um novo exemplo de análise de uma
viga em balanço representada na Figura 14.
56
Figura 14 – Viga em balanço de inércia variável.
Para fazer análise da viga da Figura 14 pelo método do princípio da mínima
energia potencial total é preciso, inicialmente, obter a expressão para a energia de
deformação de uma viga. A viga, supostamente, deve satisfazer a hipótese de
Bernoulli (1705), a qual considera que “seções transversais retas permanecem
planas e normais à tangente ao eixo fletido da viga”. O deslocamento vertical do eixo
da viga ao longo do comprimento é descrito pela função ( )v x . Da resistência dos
materiais, sabe-se que a deformação longitudinal ,x y no ponto da seção x e cota
y é dada por:
, ''x y y v x (121)
Sendo,
2
2''
d v xv x
dx
(122)
A energia de deformação específica de um material linear elástico com
módulo de elasticidade E , é dada por:
2
00 0
1
2U d E d E
(123)
Para um ponto da seção x e cota y da viga à flexão:
57
2
0
1, ''
2U y v x E y v x
(124)
A energia de deformação da viga pode ser obtida por:
2
0
1 '' 2
L
AU v x E y v x dAdx
(125)
Ou
2
0
1 '' 2
L
U v x E I v x dx (126)
Onde L é o comprimento da viga e I o momento de inércia da seção da viga, dado
por:
2 A
I y dA (127)
Como no exemplo em estudo, a inércia da seção varia ao longo do
comprimento a energia potencial total da viga pode ser obtida por:
0
1
105
5 02 1 '' ''
2
1
2b xa v x dx EI v x dxv x E Pv xI
(128)
Observando a expressão (128), verifica-se que a energia potencial total da
viga é função da função que descreve a deformação do eixo da viga v x , ainda
desconhecida. Uma função de função é denominada funcional.
Do ponto de vista matemático o problema anterior da treliça era um problema
de minimização de uma função de duas variáveis. O problema da viga é um
problema de minimização de um funcional da função v x . Trata-se agora de
58
encontrar a função v x e não mais apenas as variáveis 1d e 2d que minimizam .
Esse é um problema clássico de cálculo variacional.
Como então resolver o problema da viga à flexão? É aqui que surge a ideia
básica do método de Rayleigh-Ritz: a função v x que representa a elástica da viga
é descrita por uma função aproximadora.
As funções aproximadoras devem satisfazer as seguintes condições:
a. Devem ser funções polinomiais ou trigonométricas que satisfaçam às
condições de contorno em deslocamento da viga.
b. Devem ter derivadas contínuas até a ordem 1n , sendo n a maior ordem de
derivação da função no funcional (no caso 2n ).
c. Devem ser definidas em todo o domínio do problema.
A solução “exata” para o deslocamento na extremidade livre da viga da Figura
14 é 1875.
Primeira tentativa:
A primeira função aproximadora adotada é um polinômio de segundo grau.
2
1v x x (129)
Vale observar que a função satisfaz às condições de contorno em
deslocamento do problema:
0
) 0x
a v x (130)
0
b) ' 0x
v x (131)
59
Substituindo
1'' 2v x (132)
Na expressão (128), e integrando-se, chega-se a:
2
1 1 130 1000 (133)
Vale observar que agora é uma função do parâmetro 1 e não mais da
função v x . Isso significa que o problema a ser resolvido é um problema de mínimo
de função e não mais de mínimo de um funcional. Essa é a contribuição do método
aproximado de Rayleigh-Ritz.
Aplicando-se agora o princípio da mínima energia potencial total, o qual
afirma que a configuração deformada minimiza a energia potencial total de uma
estrutura em equilíbrio estável, obtém-se:
1
1
1
0 16,66d
d
(134)
Logo
216,66v x x (135)
E, Portanto,
10
1666x
v x
(136)
Observa-se que o erro no cálculo de em relação à solução exata é muito
grande:
60
1875 166611,1%
1875erro
(137)
Da resistência dos materiais sabe-se que:
''M x El v x (138)
Assim, no trecho (a),
0 5
2 2 16,66 66,66a xM x x x
(139)
5 10
1 2 16,66 33,33b xM x x x
(140)
A Figura 15 compara os momentos da solução aproximada e da solução
correta (viga isostática). Os momentos são constantes ao longo de x nos dois
trechos porque v x é uma função do segundo grau.
Figura 15 – Diagrama de momentos na viga associado a v x definido em (135).
Observação: a solução é ruim tanto em termos de deslocamentos quanto em
termos de momentos. A aproximação dos momentos é ainda pior porque ela é
obtida de derivadas de funções aproximadoras.
Segunda tentativa:
No problema estudado a solução é muito simples porque a viga é isostática.
No caso de uma viga altamente hiperestática de vários vãos com inércias diferentes
em cada vão e cargas distribuídas, a solução não é trivial e não estará disponível
para se saber se a solução aproximada é boa ou não. Nesse caso, o procedimento a
seguir é usar uma função aproximadora mais “rica” e verificar a mudança na
61
resposta. Quando, ao se refinar a solução, a resposta não melhora
significativamente, a solução anterior já pode ser considerada boa.
Na segunda tentativa, a função aproximadora é um polinômio do terceiro grau
dado por:
2 3
1 2v x x x (141)
Vale observar que a função satisfaz às condições de contorno em
deslocamento (130) e rotação (131).
Substituindo
1 2'' 2 6v x x (142)
Em (128) e integrando-se, chega-se a:
2 2
1 2 1 1 2 2 1 2, 30 750 6750 1000 10000 (143)
Vale observar que P agora é uma função dos parâmetros 1 e 2 . Aplicando-
se o princípio da mínima energia potencial total, obtém-se:
1 2
1
,0
(144)
1 2
2
,0
(145)
Que fornece,
1 2
800 20
33 33e
(146)
62
Logo,
2 3800 20
33 33v x x x
(147)
10
1818x
v x
(148)
Usando-se (138), chega-se a:
0
16002 97,0
33a x
M
(149)
5
1600 120 52 60,6
33 33a x
xM
(150)
5
1600 120 530,3
33 33b x
xM
(151)
10
1600 120 1012,1
33 33b x
xM
(152)
A comparação entre os momentos da solução aproximada e da solução exata
(viga isostática) está apresentada na Figura 16.
Figura 16 – Diagrama de momentos na viga associado a v x definido em (146).
63
Observações:
1) A solução melhorou significativamente em termos de deslocamentos, mas
continua ruim em termos de momentos. Não é coincidência que o
deslocamento na extremidade livre seja inferior ao da solução exata, pois a
aproximação torna a estrutura mais rígida.
2) O problema na descontinuidade no diagrama de momentos na solução
aproximada continua. A descontinuidade acontece porque v x e,
consequentemente, sua segunda derivada, é contínua no domínio enquanto
que a rigidez EI é descontínua em 5x .
3) O problema identificado revela uma limitação do método de Rayleigh-Ritz que
é o de trabalhar com apenas uma função contínua no domínio. Para se
superar o problema é preciso usar duas funções, uma no trecho (a) e outra no
trecho (b), impondo condições de continuidade em 5x para v x e para sua
primeira derivada em relação a x , mas, liberando a curvatura para ser
descontínua.
Terceira tentativa:
Serão usadas duas funções cúbicas aproximadoras, uma para o trecho (a) e
outra para o trecho (b):
2 3
1 2 0 5av x x x x (153)
2 3
3 4 5 6 5 10bv x x x x x (154)
Vale observar que a função av x satisfaz às condições de contorno em
deslocamento definidas em (130) e (131). Além disso, serão impostas as seguintes
condições de continuidade em 5x .
5 5a bx x
v v (155)
64
5 5
' 'a bx xv v
(156)
Essas duas condições permitem reduzir o número de parâmetros incógnitos
de 6 para 4. Os parâmetros 5 e 6 , por exemplo, podem ser escritos em função
dos outros parâmetros.
Aplicando-se o princípio da mínima energia potencial total, obtém-se:
1 2 3 4
1
, , ,0
(157)
1 2 3 4
2
, , ,0
(158)
1 2 3 4
3
, , ,0
(159)
1 2 3 4
4
, , ,0
(160)
É possível obter os parâmetros 1 , 2 , 3 e 4 que, substituídos em (153) e
(154), fornecem:
2 310
2512
av x x x (161)
2 210 300 1000 800 10 1000 300
25(10 25) 25 1012 4 4 16 6 4 4
bv x x x x x x x
(162)
Nota-se que
10
1875b xv
(163)
65
é a solução exata para o deslocamento na extremidade livre. O diagrama de
momentos correspondentes às expressões (161) e (162) também é exato.
Observações:
1) O uso de duas funções aproximadoras av x e bv x permitiu obter a solução
exata do problema porque foi possível representar a descontinuidade que
existe na derivada segunda da função elástica em 5x . O procedimento
usado na terceira tentativa foi o de melhorar a precisão da solução usando
duas funções aproximadoras, uma para cada trecho da viga, em vez de
continuar a aumentar o grau de polinômio da função v x no domínio de 0 a
L . Mesmo usando um polinômio do quarto grau para v x não se pode obter
a solução exata porque haverá ainda uma descontinuidade na segunda
derivada de v x o que causará uma descontinuidade no diagrama de
momento, uma vez que há uma descontinuidade na rigidez EI da viga.
2) Posto como está, o método de Rayleigh-Ritz ainda não é um método dos
deslocamentos, no sentido clássico, porque as incógnitas não são os
deslocamentos.
3) Com o uso de duas funções no domínio o método deu um grande passo para
se aproximar do método dos elementos finitos. Na verdade o domínio foi
“discretizado” em dois subdomínios, ou elementos.
4) Para transformar definitivamente o método de Rayleigh-Ritz no MEF, o
método de Rayleigh-Ritz precisa substituir as incógnitas i pelos graus de
liberdade da estrutura id .
7 – O MÉTODO DOS ELEMENTOS FINITOS
Nessa fase do tratamento do problema será mais cômodo trabalhar com
parâmetros nodais e não com parâmetros generalizados ( )i . Parâmetros nodais,
conforme já acenado anteriormente, são valores da função procurada, ou de suas
derivadas, em pontos específicos do domínio ou fronteira do elemento.
66
Se o procedimento apresentado neste item for adotado, o resultado poderá
ser atendido como válido não apenas para a viga toda mas, também, para trechos
dessa viga, ou elementos finitos, para adotar nomenclatura que não os confunda
com os elementos infinitesimais do tratamento clássico via equações diferenciais. No
caso do elemento finito associado à aproximação cúbica os parâmetros nodais
poderão ser os seis da Figura 17, três em cada nó de extremidade. Note-se que
nessa figura os deslocamentos lineares e angulares estão exageradamente
ampliados, como acontecerá ainda em outras figuras, apenas com o intuito e facilitar
sua visualização.
Fica fácil observar que eN elementos desse tipo poderão ser acoplados pelas
extremidades, passando a ter nesses pontos de conexão valores comuns para v e
. Desse modo, neste caso, a discretização não estará violando a compatibilidade
cinemática entre elementos (deslocamentos e rotações idênticas à esquerda e à
direita do nó). Uma viga assim composta não teria na sua discretização a introdução
de descontinuidades de v ou entre elementos adjacentes. No caso particular da
viga, aqui utilizado apenas com a intenção de veicular conceitualmente as bases do
método dos elementos finitos, como já foi visto, a função aproximadora é a própria
função exata para o estado de deslocamento do elemento sob ações de
extremidade. Logo, neste caso, a resposta poderá ser de boa qualidade. (SAVASSI,
W. 1996)
7.1 - Dedução com utilização da linguagem matricial
Agora será feita a dedução do sistema de equações lineares algébricas com a
utilização da linguagem matricial, muito apropriada para a correspondente
elaboração de programas para computadores. O caso a considerar será o de viga
prismática.
Seja um elemento finito dessa viga, cujo comportamento pretende-se
examinar com a adoção da mesma aproximação correspondente a polinômio cúbico
para a representação da linha elástica. (SAVASSI, W. 1996)
67
Figura 17 – Elemento finito de viga.
Adotando-se como polinômio:
2 3
1 2 3( ) ov (164)
Uma elástica elementar da barra de pórtico plano isolada é definida no
sistema de eixos locais pelo deslocamento axial u(x) e pelo deslocamento
transversal v(x), devido à adoção da hipótese de pequenos deslocamentos, o
comportamento axial e o comportamento transversal de uma barra são considerados
independentes.
Dessa forma, o deslocamento axial u(x) só depende das deslocabilidades axiais d’1 e
d′4, e o deslocamento transversal v(x) fica definido somente pelas deslocabilidades
d′2, d′3, d′5 e d′6.
No nosso caso, os deslocamento axiais d’1 e d’4 não serão considerados, d’2, d’3, d’5,
d’6 serão adotados como 1v , 1 , 2v , 2 respectivamente; para designar os nós na
numeração local e por facilitar a identificação do deslocamento e giro.
A função pode ser escrita pela seguinte matriz
0
12 3
2
3
( ) 1v
(165)
onde é a matriz (1x4) de funções monômias de e é o vetor (4x1) de
parâmetros generalizados.
68
Para uso posterior, calcule-se:
20 1 2 3dv
d
(166)
2
20 0 2 6
d v
d
(167)
É preferível passar a trabalhar com parâmetros nodais nv , em lugar dos
parâmetros generalizados .
defina-se nv como :
.
1 1 2 2
n tv v v (168)
Procura-se agora determinar outra matriz:
( ) nv v (169)
Sabe-se que:
1
0 0 01
2
2
1 1 1
(0)(0) (0)
1
(1) (1) (1)
1
n
vv v
v dvdv dv d
l ddx d dxv
v v v v
dv dv d dv
dx d dx l d
(170)
Ou
1
*1
2
2
1 0 0 0
0 1 0 0
1 1 1 1
0 1 2 3
n
v
lv
v
l
(171)
69
Para as deslocabilidades transversais, as equações que definem as funções
de forma são obtidas a partir da matriz v determinando os valores das constantes
0 , 1 , 2 , 3 com base em condições de contorno adequadas.
1v é definida considerando (0)v ,
1l é definida considerando 1
0
dv
d
2v é definida considerando (1)v
2l é definida considerando 2
1
dv
d
simbolicamente:
*
n
v A (172)
então:
*1
n
A v
(173)
Logo substituindo em:
*1( )
n
v A v
(174)
onde a matriz inversa, é dada por:
1
1 0 0 0
0 1 0 0
3 2 3 1
2 1 2 1
A
(175)
Resultando em:
70
0 1
1 1
2 1 1 2 2
3 1 1 2 2
3 2 3
2 2
v
l
v l v l
v l v l
(176)
Sendo portanto:
2 3 2 3 2 3 2 3
1 1 2 2 1 3 2 2 3 2 v v l v l (177)
Obtém-se, sob forma explícita:
1
12 3 2 3 2 3 2 3
2
2
1 3 2 2 3 2
v
lv
v
l
(178)
Ou
* *
n
v v (179)
Ou
nv v (180)
com
2 3 2 3 2 3 2 3 1 3 2 2 3 2 l l (181)
Aplique-se agora
0p U (182)
71
Define-se a energia de deformação específica, ou energia por unidade de volume
0 / 2u . Logo, a energia de deformação em um elemento estrutural de volume V
será:
1 .2
v
U dV (183)
Pela lei de Hooke, elemento fletido que foi visto em (2), .M
yE EI
Pela teoria da linha elástica:
2
2 .
d v M
dx EI
(184)
Logo
2
2 .
d vy
dx
(185)
Então a energia de deformação, como apresentado em 6.5.1 será:
1
2v
U E dv (186)
Sendo dV dydzdx dAdx , e de acordo com Hibbeler (2004), o momento de inércia
para toda área é determinado por integração, como visto em (127):
2 A
Ix y dA (187)
então:
22 2
2
0 0
1 1 .
2 2
l l
A
d vU Ey dAdx EIv dx
dx
(188)
Para ser representada na linguagem matricial é feito:
72
2
2 .x x
d vy yk
dx
(189)
Levando-se em U:
0
1 k .
2
l
t
x xk E xI d
(190)
Fazendo DEI e lembrando-se que:
2 2
2 2 2 2
1 1 ,x
d v d vk k
dx l d l
(191)
e calculando-se posteriormente k em função de parâmetros nodais:
22
2 2 . .n n
dd vk v B v
d d
(192)
Então:
1
.t
3
0
1 1 B D .2
n t nv B d vUl
(193)
É feita a integração e chama-se o resultado de e ik ( o índice ei indica o elemento i),
tem-se que:
.
1 .
2
n t n
e iU v k v (194)
A variação da energia de deformação será:
.m
m
UvU
v
(195)
Indicando mv um deslocamento nodal genérico, dentre os quatro de nv .
A seguinte passagem de matrizes posteriormente terá um exemplo para facilitar a
visualização.
Considerando-se apenas o primeiro termo do somatório (sem o 1v ):
73
1 1
1 1
1 1 2 2 1 1 2 2 1 2 21 1
2 2
1 1
2 2e i e i
v v
Uv v k v v k
v vv v v
1
1
1 1 2 2
2
2
1
01 1 1 0 0 0
02 2
0
e i e i
v
k v v kv
11
21 .
11 12 13 14
31
41
1 1 v v .2 2
n n t
k
kk k k k
k
k
(196)
Como e ik é simétrica, se 1k indicar a primeira linha de e ik :
. .
1 1 1 1
1
1 1 v k .
2 2
n n t t n n tUk v k v v k
v
(197)
Portanto:
1 1 2 2 3 3 4 4 v v v v n n n n
m
m
Uv k v k v k v k v
v
1
2
1 2 3 4
3
4
n
n
n
n
k v
k vv v v v
k v
k v
1
. . 2
3
4
v .
n
n
n t n n t n
e in
n
k v
k vv v k v
k v
k v
(198)
7.1.1 – Exemplo
Como temos:
74
1
12 3 2 3 2 3 2 3
2
2
1 3 2 2 3 2
v
lv
v
l
(199)
A derivada segunda dessa função é:
1 1 2 2 12 6 6 4 12 6 6 2, .v l v lv (200)
Sendo:
1
2
3
0
U , 2
EIv d
l
(201)
A nossa variação de energia de deformação U para os quatro parâmetros nodais
será:
1
3
1 0
, , (12 6)
,
vU EIv d
v v l
,
1
3
1 0
, , (6 4)
,
vU EIv l d
v l
1
3
2 0
, , ( 12 6)
,
vU EIv d
v v l
,
1
3
2 0
, , (6 2)
,
vU EIv l d
v l
(202)
Fazendo-se a substituição de ,v , multiplicando e integrando resultará:
1 1 2 23
1
, 12 6 12 6
,
vU EIv l v l
v v l
,
1 1 2 23
2
, 6 4 6 2
,
vU EIl v l v l
v v l
1 1 2 23
1
, 12 6 12 6
,
vU EIv l v l
v l
,
1 1 2 23
2
, 6 2 6 4
,
vU EIl v l v l
v l
(203)
75
Esses resultados são os deslocamentos nodais representados por:
1 1 2 2 3 3 4 4 v v v vn n n nk v k v k v k v na linguagem matricial.
Finalmente se obtém a matriz de rigidez do elemento de viga:
1
2 2
1
3
2
2 2
2
12 6 12 6
6 4 6 2
12 6 12 6
6 2 6 4
vl l
l l l lEI
vl ll
l l l l
(204)
Que é representada como .
n t n
e iv k v na linguagem matricial.
7.2 – Energia potencial externa
Como demonstrado em 6.5.2, as ações (forças e momentos) deformam o
corpo, fazendo com que seus infinitos pontos passem da posição indeformada
(posição inicial) para uma posição deformada (posição final). Em relação à posição
inicial, essas ações têm uma energia potencial que é igual ao trabalho que elas
realizariam para levar o corpo à sua posição inicial sem carga.
* eP f W (205)
A essa energia chamamos de energia potencial externa, é dada por:
1
n
i i
i
Pw
(206)
onde iP representa as ações genéricas, iw os deslocamentos genéricos e n, o
número de esforços. O sinal negativo indica que cada ação realiza trabalho negativo,
ao retornar da posição carregada para a posição sem carga. Essa parcela da
energia potencial não tem 1
2 multiplicando o segundo membro da igualdade porque
ela é o trabalho dos esforços atuantes com seus valores finais, quando a estrutura é
movida para sua posição inicial.
76
Para uma carga concentrada:
.x l
v xp
(207)
Para uma carga distribuída ao longo do elemento:
0
.
l
p v x dx (208)
Se a força é linearmente distribuída, sua expressão será, em função de coordenada
adimensional , para facilitar o entendimento, já que a nossa função aproximadora
já encontrada também está em função de .
Figura 18 – Força externa linearmente distribuída.
1
2
( ) 1 n
p
pp p
p
(209)
Onde 1p e 2p serão valores nodais conhecidos.
Com
1
0 0
.
l
t t n
pv p dx l v d p (210)
vamos ter:
1
.
0
.n t t n
pl v d p (211)
Portanto:
1
0
f t n n
ei pl d p H p (212)
77
Onde
2 3
2 31
2 3
0
2 3
1 3 2
2 1 .
3 2
lH l d
l
(213)
Efetuando o produto, integrando e multiplicando por np resultará:
1 2
2
1 2
1 2
2
1 2
0,35 0,15
1,5 1,0 / 30f .
0,15 0,35
1,0 1,5 / 30
e i
p p l
p p l
p p l
p p l
(214)
Para ser usada com uma força uniformemente distribuída basta igualar 1 2p p p :
2
2
/ 2
/12f .
/ 2
/12
e i
pl
pl
pl
pl
(215)
Considerando a viga da Figura 14, utilizando as funções de interpolação de
viga i x na análise de uma viga. Observando que se deve usar 5L
(comprimento de cada trecho) nas expressões de i x para se obter as funções
av x e bv x , pode-se escrever:
3 1 4 1av x x v x (216)
1 1 2 1 3 2 4 2bv x x v x x v x (217)
78
Considerando as duas funções distintas av x e bv x , respectivamente nos
trechos a e b , e integrando-as em x de 0 a 5L em cada trecho (subdomínio
do trecho ou elemento finito) e observando-se que a força P atua no sentido
negativo da direção de 2v , pode-se escrever a expressão da energia potencial total
como:
5 52 2
1 1 2 2 20 0
1 1, , , '' ''
2 2a a b bv v EI v x dx EI v x dx Pv
(218)
Substituindo na expressão (220) aEI , bEI e P pelos seus valores numéricos,
efetuando as integrais e usando o princípio da mínima energia potencial total como
descrito em (110), obtém-se:
1
1
2
2
36 6 12 6
125 25 125 2506 12 6 2
025 5 25 5
12 6 12 6 10
125 25 125 25 0
6 2 6 4
25 5 25 5
v
v
(219)
Ou,
Kd f (220)
Sendo K a matriz de rigidez da viga, d o vetor dos deslocamentos nodais e f o
vetor de cargas nodais. Essa solução é “exata” e coincide com a última solução
obtida para o método de Rayleigh-Ritz.
3 1875d (221)
7.3 – Exemplos – flechas
79
O diagrama de deflexão do eixo longitudinal que passa pelo centroide de cada
área da seção transversal da viga é denominado linha elástica. A linha elástica da
maioria das vigas é esquematizada sem nenhuma dificuldade. Ao fazer o diagrama,
entretanto, é necessário saber como os vários tipos de apoio limitam a inclinação ou
deslocamento.
Na análise de uma viga em balanço, supondo-se que se considere apenas um
elemento finito, de comprimento l igual ao comprimento da viga, introduzindo as
condições de contorno 1 0v e 1 0 , o que vai mudar na análise da flecha é o vetor
de cargas nodais.
Para uma viga em balanço com uma força uniforme distribuída:
Figura 19 – Viga em balanço com força uniforme distribuída
1
1
3
2
2 22
0
1 0 0 0 0
0 1 0 0
0 0 12 6 2
0 0 6 4
12
v
EI lP
vll
l l l
(222)
Resolvendo o sistema, resulta: 4
28
plv
EI e
3
26
pl
EI .
Com uma força concentrada na extremidade:
80
Figura 20 – Força concentrada na extremidade
1
1
3
2
2
2
1 0 0 0 0
0 1 0 0 0
0 0 12 6 1
0 0 6 4 0
v
EIP
vll
l l
(223)
Resulta: 3
23
plv
EI e
2
22
pl
EI .
Com uma força linearmente distribuída, sendo 1 0p :
Figura 21 – Força linearmente distribuída com P1=0
1
1
3
22
2
2
01 0 0 0
00 1 0 0
2 0,350 0 12 6
1,50 0 6 4
30
v
EIp l
vlll
l l
(224)
Resulta: 4
2
11
120
plv
EI e
3
28
pl
EI .
E com uma força linearmente distribuída, sendo 2 0p :
81
Figura 22 – Força linearmente distribuída com P2=0
1
1
3
22
2
2
01 0 0 0
00 1 0 0
1 0,150 0 12 6
1,00 0 6 4
30
v
EIp l
vlll
l l
(225)
Resulta: 4
2
1
30
plv
EI e
3
2
25
600
pl
EI .
8 – TÉCNICAS COMPUTACIONAIS PARA AUTOMATIZAÇÃO DO
MÉTODO DOS ELEMENTOS FINITOS
8.1 – Organização da montagem do sistema de equações
8.1.1 – Generalidades
Um importante tópico, mas frequentemente mal entendido, é o processo de
montagem do sistema de equações a partir das equações do elemento.
Tal montagem nada mais é do que uma operação de soma dos coeficientes das
equações dos elementos em sua localização própria dentro do sistema de
equações.
8.1.2 – Disposição dos coeficientes na matriz de rigidez global
Como a maior parte das matrizes estruturais são simétricas, somente há
necessidade de se armazenarem os coeficientes pertencentes à parte triangular
superior, ou inferior, já que nm mnk k , e os coeficientes da diagonal principal.
82
As matrizes envolvidas em problemas estruturais podem ainda apresentar a
característica de ter os seus coeficientes não nulos distribuídos e agrupados ao
longo da diagonal principal ampliada. A este tipo de distribuição dá-se o nome de
característica de banda.
A disposição dos elementos na matriz de rigidez da estrutura está diretamente
associada à numeração adequada dos pontos nodais quando da discretização
estrutural. Este é um fato que deve ser examinado com cuidado pois uma
numeração inadequada poderá ocasionar um gasto excessivo de memória.
(RODRIGUES, R. O. 1999)
Elemento de barra considerando-se somente deslocamentos axiais.
Numeração dos pontos nodais em ordem sequencial.
11 12 0 0
21 22 23 0
0 32 33 34
0 0 43 44
K
Figura 23 – Matriz de rigidez com característica de banda.
Para a numeração em ordem sequencial a característica de banda está bem
acentuada e os coeficientes nulos da matriz de rigidez encontram-se na parte
triangular superior e inferior.
Elemento de barra considerando-se somente deslocamentos axiais.
Numeração dos pontos nodais não está em ordem sequencial.
11 0 0 14
0 22 23 24
0 32 33 0
41 42 0 44
K
Figura 24 – Matriz de rigidez sem característica de banda.
Já para essa situação a característica de banda está totalmente perdida e os
coeficientes nulos da matriz de rigidez encontram-se espalhados.
Verifica-se, portanto, que a largura da faixa, região onde concentram-se os
coeficientes não nulos da matriz de rigidez, é dada em função da maior diferença
“MD” entre as numerações dos pontos nodais de um mesmo elemento estrutural e o
83
número de deslocabilidades de cada ponto nodal “NDN”. Desta forma a largura da
faixa da parte superior, ou inferior, da matriz de rigidez incluindo-se a diagonal fica
definida pela equação:
( 1)LF MD NDN (226)
8.1.3 – Processo de expansão e acumulação
A matriz de rigidez global " "sK de uma estrutura, ainda sem a consideração
da existência das condições de contorno, é obtida a partir da contribuição das
matrizes de rigidez de cada elemento, como mostra a equação: 1
nel
s s
l
K k
. Dessa
forma, cada coeficiente da matriz " "sK corresponde à somatória das contribuições
dos vários elementos da matriz de rigidez de cada elemento que correspondem a
um mesmo nó. (RODRIGUES, R. O. 1999)
Na Figura 25 apresenta-se um elemento de barra com a respectiva
numeração de seus nós, sendo somente considerados os deslocamentos nodais
axiais e transversais do elemento.
Elemento de barra considerando-se
somente deslocamento axial e transversal
1 3
i j
Figura 25 – Numeração local dos deslocamentos do elemento.
O relacionamento entre a numeração dos deslocamentos das extremidades
dos elementos com a dos deslocamentos nodais da estrutura denomina-se regra de
correspondência. Em cada elemento “e” os deslocamentos são numerados de 1 a
(NDN.NNE), sendo “NNE” o número de nós do elemento. Tal numeração deve ser
feita na forma sequencial, sendo de 1 a NDN no primeiro nó do elemento e de
e
4 2
84
(NDN.NNE–NDN+1) a (NDN.NNE) no último nó do elemento, conforme mostrado
pela figura 25.
Na estrutura, os deslocamentos são numerados na ordem dos nós, sendo que em
cada nó há “NDN” deslocamentos, em ordem determinada pelos eixos do sistema
global de referência. (RODRIGUES, R. O. 1999)
A regra da correspondência consiste em colocar o coeficiente " "nmK da matriz de
rigidez do elemento “e” contribuindo cumulativamente com o coeficiente da matriz de
rigidez global da estrutura, pois poderá haver outros elementos que contribuem na
mesma posição. O índice “n”, ou “m”, de cada coeficiente de rigidez é dado pelo
esquema:
e = 1@NELEM
A = 1@NNE
B= 1@NDN
n (NDN.nº nó )e
A NDN B
Efetuando-se esse desenvolvimento para todos os elementos
sucessivamente, pode-se armazenar os coeficientes da parte superior e da diagonal
principal da matriz de rigidez global, correspondendo, assim, os coeficientes em que
o segundo índice seja igual ou maior que o primeiro. Essa operação é muitas vezes
designada como processo de expansão e acumulação dos coeficientes das matrizes
dos elementos pelos endereços correspondentes da matriz de rigidez global.
8.1.4 – Matriz de rotação
Antes de se aplicar o processo de expansão e acumulação em cada matriz de
rigidez do elemento para a obtenção da matriz de rigidez da estrutura, se faz
necessário rotacionar a matriz elemental em função da posição original do elemento
na estrutura, conforme será mostrado na sequência. (RODRIGUES, R. O. 1999)
Decompondo-se as forças e os deslocamentos nodais locais em relação aos
eixos globais, obtêm-se:
85
; DE Ef r d rF (227)
onde " " r é a matriz de rotação que contêm os cossenos diretores de cada
elemento.
Substituindo-se as equações (227) na equação (220), obtém-se:
DE sr F k r (228)
Multiplicando-se os dois termos da igualdade pela transposta da matriz " " r ,
obtém-se:
F t
E sr k r D (229)
Neste caso, conclui-se que a matriz " " tr é igual à sua inversa 1" "r , sendo
que para pórticos planos tal matriz é definida por:
cos cos 0 0 0 0
cos cos 0 0 0 0
0 0 1 0 0 0r
0 0 0 cos cos 0
0 0 0 cos cos 0
0 0 0 0 0 1
(230)
onde " " é o ângulo formado entre o eixo local "x" e o eixo global "X" e " " é o
complemento de " " .
8.1.5 – Armazenamento computacional da matriz de rigidez
Uma forma de armazenamento muito eficiente e de fácil aplicação
corresponde à denominada técnica de armazenamento em altura efetiva de coluna
ou skyline. Esta técnica corresponde a armazenar dentro de um vetor de trabalho
principal “S” as colunas da parte triangular superior da matriz e os elementos da
diagonal principal, a partir do primeiro elemento não nulo de cada coluna. Este
armazenamento é feito em forma sequencial por coluna, de cima para baixo.
(RODRIGUES, R. O. 1999)
86
A matriz de rigidez da estrutura da figura 24 tem os seus coeficientes espalhados
segundo a sequência de numeração elaborada, podendo ser reproduzida com mais
detalhe na figura 26, cujos coeficientes do perfil estão numerados sequencialmente,
segundo as posições a serem ocupadas dentro do vetor de trabalho principal.
1 5
2 3 6
4 7
8
K
Figura 26 – Perfil para armazenamento por altura efetiva de coluna.
Para esse tipo de armazenamento é necessário ainda um vetor auxiliar que
indique, dentro do vetor “S”, as posições dos elementos da diagonal principal da
matriz. Este vetor auxiliar está representado na tabela 1, com a denominação de
vetor “IPOS”.
Tabela 1- Vetor auxiliar contendo as posições dos elementos da diagonal principal.
Posições do vetor “IPOS” 1 2 3 4
Posições dos coeficientes da diagonal principal da matriz
no vetor “IPOS” 1 2 4 8
Uma vez definido o vetor “IPOS”, qualquer elemento da matriz de rigidez
original fica perfeitamente localizado dentro do vetor de trabalho principal “S” através
da equação (231), onde o elemento " "ijK ocupará a posição " "S l .
l IPOS j i j (231)
No tratamento de grandes sistemas de equações pode ocorrer o caso crítico
em que o vetor de trabalho principal supera o limite estabelecido pelo equipamento
computacional utilizado, sendo que neste caso a solução está em particionar o perfil
em grupos de colunas. A seguir, apresentam-se as figuras 27 e 28 que contêm os
87
fluxogramas para a obtenção do vetor auxiliar e montagem da matriz de rigidez
global, dentro do vetor de trabalho principal, com a seguinte notação:
NNODS - número de pontos nodais da estrutura.
NELEM - número de elementos da estrutura.
NNE - número de pontos nodais por elemento.
NDN - número de deslocamento em cada ponto nodal.
NDE - número de deslocamentos por elemento = NNE.NDN INC(NELEM, NNE) -
matriz que contém em cada linha “i” a lista de incidência (numeração dos nós) do
elemento de numeração “i”.
SG (NDE, NDE) - matriz de rigidez do elemento no referencial global.
NEQ - número de deslocamentos da estrutura ou número de equações do sistema =
NDN.NNODS.
JK(NDE) - vetor que faz a correspondência entre a numeração “j” de cada elemento
com a numeração JK(j) da estrutura.
IPOS(NEQ) - vetor auxiliar que contêm os coeficientes da diagonal principal da
matriz de rigidez.
NP - número de posições do vetor de trabalho = IPOS(NEQ).
A(NP) - vetor de trabalho que contém a matriz de rigidez armazenada por altura
efetiva de coluna.
88
Figura 27 - Fluxograma de obtenção do vetor “IPOS”.
MONTAGEM DO VETOR “IPOS”
NEQ = NDN . NNODS
DIM IPOS(NEQ)
PRIMEIRA FASE DA OBTENÇÃO DO VETOR “IPOS”
I = 1 , NELEM
LL = ORDEM DA MENOR NUMERAÇÃO DO PONTO NODAL DO ELEM. I
L = NDN . (INC(I,LL) – 1) + 1
J = 1 , NNE
K = 1 , NDN
M = NDN . (INC(I,J) – 1) + K
IDIF = M – L +1
IPOS(M) < IDIF
IPOS(M) = IDIF
SEGUNDA FASE DA OBTENÇÃO DO VETOR “IPOS”
IC = 2 , NEQ
IPOS(IC) = IPOS(IC – 1) + IPOS(IC)
SIM
NÃO
89
Figura 28 - Fluxograma da montagem da matriz de rigidez por altura efetiva da coluna.
NP = IPOS(NEQ)
DIM A(NP), JK(NDE), SG(NDE,NDE)
I = 1 , NELEM
CÁLCULO DA MATRIZ DE RIGIDEZ SG PARA O ELEMENTO I
IC = 0
J = 1 , NNE
K = 1 , NDN
IC = IC + 1
M = NDN . (INC(I,J) – 1) + K
JK(IC) = M
J = 1 , NDE
K = J , NDE
NCO = JK(K) – JK(J) + 1
NCO > 0
L = IPOS(JK(K)) + JK(J) – JK(K)
A(L) = A(L) + SG (J,K)
SIM
NÃO
90
8.1.6 – Armazenamento do vetor de forças nodais
Após a determinação da matriz de rigidez é necessário considerar o
carregamento sobre a estrutura, determinando assim o vetor de forças nodais global,
sendo conveniente tratar separadamente as forças aplicadas nos nós e nos
elementos. As forças aplicadas aos nós podem ser diretamente alocadas no vetor de
forças nodais global, porém as forças aplicadas aos elementos deverão ser
transformadas em forças equivalentes aplicadas aos nós do elemento e
posteriormente alocadas no vetor de forças nodais global. (RODRIGUES, R. O.
1999)
Utilizando-se a regra da correspondência obtêm-se os coeficientes do vetor força “Fi”
no qual contribuem tanto as forças aplicadas diretamente nos nós como as forças
equivalentes oriundas dos carregamentos nos elementos. A seguir, apresenta-se a
Figura 29 que contém o fluxograma para a obtenção do vetor força global que
compõe o sistema de equações, com a seguinte notação:
NNOCA - número de nós carregados.
NELCA - número de elementos carregados.
NCC - número de carregamentos = NNOCA+NELCA.
F(NEQ) - vetor das forças nodais global.
NNO(NNOCA) - vetor dos nós carregados.
CN(NNOCA,NDN) - matriz das forças aplicadas aos nós.
NEL(NELCA) - vetor dos elementos carregados.
AX(NDE) - vetor das forças nodais equivalentes.
91
Figura 29 - Fluxograma da montagem do vetor de forças nodais global.
DIM F(NEQ), AX(NDE)
CONTRIBUIÇÃO DAS FORÇAS APLICADAS NOS NÓS
NNOCA = 0
J = 1 , NNOCA
MA = NDN . (NNO(J) – 1)
K = 1 , NDN
M = MA + K
F(M) = CN(J,K)
CONTRIBUIÇÃO DAS FORÇAS APLICADAS NOS ELEMENTOS
FIM NELCA = 0
J = 1 , NELCA
OBTENÇÃO DO VETOR DE FORÇAS NODAIS EQUIVALENTES
K = 1 , NNE
M = INC(NEL(J),K) N = (K – 1) . NDN
L = 1 , NDN
IJ = N + L JI = (M – 1) . NDN + L
F(JI) = F(JI) + AX(IJ)
SIM
NÃO
SIM
NÃO
92
8.2 – Consideração das condições de contorno e dos deslocamentos
prescritos
8.2.1 – Introdução
Em Mecânica dos Sólidos Deformáveis, as condições de contorno de um
sistema estrutural podem ser de dois tipos: na parte " "S do sólido são prescritas as
forças de superfícies e na parte " "uS do sólido são conhecidos os deslocamentos.
No caso de estruturas discretizadas através de elementos finitos, o
correspondente às forças de superfície no contorno são as forças aplicadas nos nós
de forma direta ou equivalente, conforme visto no item anterior. Já os deslocamentos
em determinados nós podem ser restritos ou prescritos em direções
preestabelecidas. Sendo assim, para complementar a introdução das condições de
contorno no problema resta apenas considerar a vinculação entre o sistema
estrutural e o meio externo. (RODRIGUES, R. O. 1999)
8.2.2 – Técnicas para consideração das condições de contorno dos vínculos
Com a montagem da matriz de rigidez global ainda sem levar em
consideração a existência dos vínculos, pode-se escrever a equação do sistema,
(232) onde o vetor de deslocamentos global compreende tanto deslocamentos
incógnitos nas direções livres como deslocamentos conhecidos nos apoios. O vetor
de forças nodais global, por sua vez, compreende tanto as forças aplicadas nos nós
livres, como as aplicadas aos nós restringidos. (RODRIGUES, R. O. 1999)
Ks ED F (232)
Dessa forma, o sistema de equações dado pela equação (232) não pode ser
resolvido da forma como foi montado, pois a matriz de rigidez global é singular. Para
que o sistema de equações possa ser resolvido, é necessário considerar os
vínculos, reduzindo assim a matriz " "sK à matriz *" "sK não mais singular, e eliminar
93
dos vetores "D" e "F "E os deslocamentos prescritos e restritos, reduzindo o sistema
de equações para a forma apresentada pela equação:
* * *Ks ED F (233)
Existem três técnicas para se introduzir as condições de contorno dos
vínculos. A primeira, chamada técnica da reordenação é a única natural, porém mais
difícil de ser programada e não será abordada no presente momento. As outras
duas, chamada técnica dos zeros e um e técnica do número muito grande
correspondem a artifícios engenhosos e são geralmente preferidas pela facilidade de
se programar. (RODRIGUES, R. O. 1999)
8.2.3 – Técnica dos zeros e um
Para transformar o sistema de equações singular, (232), em um sistema não
singular, (233), esta técnica utiliza o artifício de substituir em uma direção restringida
" "i , todos os coeficientes " "ijK da linha " "i , com j i , e da coluna " "j , com j i ,
da matriz *" "sK por zeros e o coeficiente " "iiK da diagonal principal por um, e no
vetor das forças substitui-se o coeficiente " "iF da linha " "i por zero.
11 12 1 1 1
21 22 2 2 2
1 2
.. 0 ..
.. 0 ..
.. .. .. .. .. .. .. ..
0 0 .. 1 .. 0 0
.. .. .. .. .. .. .. ..
.. 0 ..
n
n
i
n n nn n n
K K K D F
K K K D F
D
K K K D F
(234)
No caso de haver deslocamentos prescritos *" "iD na direção " "i , a matriz de
rigidez é modificada conforme foi mencionado anteriormente, porém deve-se efetuar
no vetor das forças a substituição do coeficiente "F"i por *" "iD e nos demais
subtrairo valor de *"K . "iji D , conforme mostra a equação:
94
*11 12 1 1 1 1
*21 22 2 2 2 2
*
*1 2
.. 0 ..
.. 0 ..
.. .. .. .. .. .. .. ..
0 0 .. 1 .. 0
.. .. .. .. .. .. .. ..
.. 0 ..
n i i
n i i
i i
n n nn n n ni i
K K K D F K D
K K K D F K D
D D
K K K D F K D
(235)
Com essas alterações, repetidas sucessivamente para todas as direções
vinculadas, obtém-se um sistema de equações de número igual ao número de linhas
da matriz *"K "s , e as equações correspondentes às linhas das direções vinculadas
são substituídas por identidades "1* 0"iD ou *"1* "i iD D no caso de
deslocamentos prescritos.
Resolvido o sistema obtém-se todos os deslocamentos nodais, inclusive os
prescritos, porém, em virtude dos erros de truncamento acumulados, os
deslocamentos prescritos assim obtidos podem ser ligeiramente diferentes dos
valores fornecidos. (RODRIGUES, R. O. 1999)
8.2.4 – Técnica do número muito grande
Esta técnica é a mais simples e cômoda para programação, porém ela não é
exata. Consiste simplesmente em substituir o coeficiente da diagonal " "iiK de uma
direção restringida " "i por um número muito grande, da ordem de 1012 e 1015,
conforme mostra a equação:
11 12 1 1 1 1
21 22 2 2 2 2
15
1 2
1 2
.. ..
.. ..
.. .. .. .. .. .. .. ..
.. 10 ..
.. .. .. .. .. .. .. ..
.. ..
i n
i n
i i ii in i i
n n ni nn n n
K K K K D F
K K K K D F
K K K K D F
K K K K D F
(236)
95
Cabe ressaltar que sendo este número grande convenientemente escolhido
os resultados são praticamente exatos.
No caso de haver deslocamentos prescritos na direção " "i , a matriz de
rigidez é modificada para
11 12 1 1 1 1
21 22 2 2 2 2
15 15 *
1 2
1 2
.. ..
.. ..
.. .. .. .. .. .. .. ..
.. 10 .. 10
.. .. .. .. .. .. .. ..
.. ..
i n
i n
i i ii in i ii i
n n ni nn n n
K K K K D F
K K K K D F
K K K K D K D
K K K K D F
(237)
8.2.5 – Apoios Elásticos
A consideração de apoios elásticos consiste em somar ao coeficiente da diagonal
" "iiK ,de uma direção " "i em que se tenha o apoio elástico, o valor da constante
elástica da mola “Kmola”, conforme mostra a equação.
11 12 1 1 1 1
21 22 2 2 2 2
1 2
1 2
.. ..
.. ..
.. .. .. .. .. .. .. ..
.. ..
.. .. .. .. .. .. .. ..
.. ..
i n
i n
i i ii mola in i i
n n ni nn n n
K K K K D F
K K K K D F
K K K K K D F
K K K K D F
(238)
Na determinação da força que atua na mola aplica-se simplesmente a lei de
Hooke para o caso de ser linearmente elástica, conforme mostra a equação
mola mola iF K D sendo que " "iD corresponde ao deslocamento da direção " "i onde
se tem o apoio elástico. (RODRIGUES, R. O.)
96
8.3 – Solução do sistema de equações
8.3.1 – Generalidades
Uma vez definido o sistema de equações de um problema qualquer, é
possível aplicar um procedimento numérico para obtenção de sua solução.
Dependendo do tipo de problema a ser analisado, as características da solução
desejada definem o tipo de procedimento que deve ser aplicado.
8.3.2 - Procedimentos de solução
Os procedimentos de solução estão relacionados com o tipo de análise a ser
efetuada, a saber:
Análise Estática: dá lugar à aplicação do procedimento de análise para o cálculo de
deslocamentos, deformações, tensões e esforços, considerando-se as solicitações
aplicadas de forma estática.
Análise Dinâmica: dá lugar à aplicação do procedimento de análise para o cálculo
das frequências e modos naturais de vibração ou para o cálculo dos deslocamentos,
deformações, tensões, esforços, velocidades e acelerações em um instante “t”
qualquer.
Análise Linear e Não-Linear: dá lugar à aplicação do procedimento de análise para
problemas com comportamento linear ou não-linear, seja físico e/ou geométrico.
No desenvolvimento a seguir somente será abordado o processo de solução para
problemas com comportamento estático linear.
8.3.3 - Implementação do método de solução para análise estática linear
O desenvolvimento das sub-rotinas para a resolução de sistemas de equações
algébricas lineares, utilizando-se o armazenamento por altura efetiva de coluna, será
efetuado através do Método de Cholesky.
97
Esse método consiste em decompor a matriz de rigidez global do sistema em
duas matrizes, conforme a seguinte equação:
K T
s st stK K (239)
onde stK é uma matriz triangular superior.
Substituindo-se a equação (239) em (232) obtém-se:
T
st st EK K D F (240)
Ou
T
st EK Y F (241)
onde
stY K D (242)
Assim, o método de Cholesky consiste em três etapas básicas, a saber:
1) Obtenção da matriz triangular " "stK correspondendo à fase de decomposição
da matriz de rigidez " "sK .
2) Resolução do sistema de equações para a obtenção do vetor " "Y ,
correspondendo a um processo de substituição.
3) Resolução do sistema de equações para a obtenção do vetor" "D ,
correspondendo a um processo de retrosubstituição.
Na sequência, são apresentados os fluxogramas das três etapas utilizadas no
método de Cholesky, aplicado a matrizes definidas positivas, com coeficientes de
rigidez simétricos, e armazenadas por altura efetiva de coluna em um arranjo
unidimensional.
98
Figura 30 - Fluxograma – Etapa de triangularização.
Etapa de triangularização
A(1) = A(1)
J = 2 , NEQ
LJ = IPOS(J) – IPOS(J-1)
JPOS = IPOS(J-1) + 1 INJ = J – LJ + 1
INJ = J
A(JPOS) = A(JPOS) / A(IPOS(INJ))
I = INJ + 1 , J
JPOS = JPOS + 1
LI = IPOS(I) – IPOS(I-1) INI = I – LI + 1
IMAX = INI
INJ > INI IMAX = INJ
IMAX > (I-1)
IPOSI = IPOS(I) – I + IMAX - 1 IPOSJ = IPOS(J) – J + IMAX - 1
K = IMAX , (I-1)
IPOSI = IPOSI + 1 IPOSJ = IPOSJ + 1
A(JPOS) = A(JPOS) – A(IPOSI) * A(IPOSJ)
I = J
A(JPOS) = A(JPOS) / A(IPOS(I))
STOP A(JPOS) < VMIN
A(JPOS) = A(JPOS)
SIM
NÃO
SIM
NÃO
SIM
NÃO
NÃO
SIM
SIM NÃO
99
Figura 31 - Fluxograma – Etapa de substituição.
Etapa de substituição
F(1) = F(1)/A(1)
I = 2 , NEQ
LI = IPOS(I) – IPOS(I-1)
INI = I – LI + 1
JPOS = IPOS(I-1)
AUX = F(I)
INI = I
K = INI , (I-1)
JPOS = JPOS + 1
AUX = AUX – A(JPOS) * F(K)
F(I) = AUX / A(IPOS(I))
SIM
NÃO
100
Figura 32 - Fluxograma - Etapa de retro substituição.
8.4 - Informações resultantes da análise
8.4.1 - Generalidades
As informações resultantes da resolução do sistema de equações
mencionado no item anterior dependem do tipo de análise efetuada. Como
exemplificação, considera-se a análise estática linear.
Etapa de retrosubstituição
I = NEQ , 2
F(I) = F(I) / A(IPOS(I))
LI = IPOS(I) – IPOS(I-1)
INI = I – LI + 1
INI = I
JPOS = IPOS(I-1)
K = INI , (I-1)
JPOS = JPOS + 1
F(K) = F(K) – A(JPOS) * F(I)
F(1) = F(1)/A(1)
SIM
NÃO
101
8.4.2 - Resultados da análise estática linear
Para o caso deste tipo de análise, pode-se obter vários tipos de resultados,
incluindo-se os deslocamentos nodais, deformações, tensões e esforços nos
elementos, bem como as reações dos vínculos.
Em certos casos é necessário obter resultados para carregamentos cujas
ações externas são combinação linear das ações de outros carregamentos. Para
problemas lineares, os resultados correspondentes a esses carregamentos poderão
ser obtidos através da mesma combinação linear, aplicada aos resultados dos
outros carregamentos, valendo a superposição dos efeitos. Assim não é necessário
que os resultados de um carregamento combinado sejam obtidos aplicando-se
novamente o processo de análise, mas sim combinando-se diretamente os
resultados dos outros carregamentos, após a aplicação do processo de análise para
os carregamentos básicos. (RODRIGUES, R. O. 1999)
8.4.3 - Deslocamentos nodais
Após a resolução do sistema de equações visto anteriormente, onde o vetor
“D” corresponde aos deslocamentos nodais do sistema estrutural, obtêm-se
deslocamentos conforme as direções do sistema de referência básico.
8.4.4 - Esforços nos elementos
Uma vez obtido o vetor “D”, utiliza-se a regra de correspondência para
relacionar os deslocamentos nodais da estrutura com os deslocamentos nodais de
cada elemento, dados no sistema de referência básico.
Na sequência, efetua-se a transformação dos deslocamentos nodais do
elemento, dado no sistema básico de referência, para o sistema de referência local
do elemento, conforme ilustra a figura 33.
Y
102
Figura 33 - Transformação dos deslocamentos nodais.
Dessa forma, os esforços nos elementos são dados pela equação E sf k d
somente para solicitações aplicadas nos nós.
Caso tenham ocorrido forças aplicadas nos elementos, estas deverão ser
somadas as forças equivalentes, conforme equação E s eqf k d f .
Cabe ressaltar que os esforços obtidos dessa forma estão no sistema local de
referência. (RODRIGUES, R. O. 1999)
8.4.5 - Reações dos apoios
Quando se utilizam as técnicas dos Zeros e Um e do Número Muito Grande,
não é possível calcular as reações de apoio diretamente, pois ficam perdidos os
coeficientes da matriz de rigidez global.
O cálculo das reações de apoio é então efetuado pelo equilíbrio de forças no nó e na
direção restringida. Assim, tomam-se todos os elementos da estrutura que
concorrem no nó onde há um apoio e compõem-se, em uma resultante, os esforços
nas direções restringidas. A reação de apoio é igual e contrária a essa resultante.
Deslocamentos nodais no
sistema básico
Deslocamentos nodais no
sistema local
j
i
x
y
X
e e
j
i D 1
D 2
D 3
D 4
d 1 d
2
d 4
d 3
103
9 – RESULTADOS
O presente trabalho apresentou o Método dos elementos finitos em linguagem
matricial e a linguagem de programação necessária que possibilitou elaborar um
programa em Visual Basic capaz de analisar uma estrutura aporticada a partir de
dados de entrada e por fim montar um relatório com os resultados numéricos que os
esforços presentes provocam na estrutura.
O pórtico plano é apresentado como o único sistema que o programa pode resolver,
como mostrado na Figura 34.
Figura 34 - Tela de abertura e seleção do modelo.
A interface do programa é didática, na aba “Dados Gerais”, figura 35, deve ser
introduzido os dados iniciais que vão definir informações relativas ao tamanho do
pilar e viga, e em quantos elementos cada um vai ser analisado.
Como exemplo a ser comparado com um programa comercial de elementos
finitos chamado SAP2000, pilar e viga tem 200cm de comprimento, e o número de
elementos para ambos será 4.
104
Figura 35 - Dados gerais.
Quando clicar em “Processar” o programa vai direcioná-lo a próxima aba,
para que sejam introduzidos mais dados.
Nos “dados nodais” são colocados as restrições nodais que vão dar a estabilidade
da estrutura, nossa estrutura possui engaste no nó 1 e nó 13, portanto é marcado X,
Y e Z, como mostra a Figura 36. Temos um pórtico bi engastado.
Figura 36 - Dados nodais.
105
Depois de processado a aba “forças” é aberta, vai ser adotada uma força de
8kN na direção X aplicada no nó 5.
O pórtico visto no SAP 2000 é fisicamente como a Figura 38.
Figura 37 - Forças.
Figura 38 - Estrutura Modelada no Programa Computacional SAP2000.
106
Na aba materiais, Figura 39, pode escolher entre concreto e aço, que ao
serem selecionados já será disponibilizado os dados E, Poisson e Fc, de cada um,
se for uma usado um material diferente é só assinalar “outro” e introduzir os dados
pertinentes. Nosso exemplo é com concreto, é clicado em “inserir” e depois em
processar para a última aba de inserção de dados.
Figura 39 - Escolha do Tipo de Material.
Na aba elementos, são apresentados os elementos e quais são os seus nós
iniciais e finais, a área e inércia de cada elemento pode ser modificada, para o
exemplo, tem-se o exposto na Figura 40.
Figura 40 - Elementos e suas características.
107
O MAPE finalizou com os deslocamentos e rotações nodais das figuras 41 e
42. Como a solução para introduzir as condições de contorno vínculos empregada
foi à técnica do número muito grande, conforme mostrado no capitulo 8.2.4, por ser
mais fácil de programar, apesar de provocar um erro muito pequeno nos resultados
finais, de modo que os deslocamentos se aproximaram de forma muito satisfatória
aos obtidos pelo já consagrado programa SAP2000.
Figura 41 - Deslocamentos nodais - Parte 1.
Figura 42 - Deslocamentos nodais - Parte 2.
108
A Figura 43 extraída do programa SAP2000 e representa o nó 4 do nosso
exemplo, U1 é a nossa translação em X; U3, é a nossa translação em Y e R2 é a
rotação em Z. Comparando com a Figura 41 pode perceber a proximidade dos
resultados.
Figura 43 - Deslocamentos do Nó 4.
109
O restante dos resultados está na tabela seguinte para comparação.
Tabela 2 – Comparação de deslocamentos e rotação entre SAP e MAPE
DESLOCAMENTOS
SAP MAPE
Nó X (cm) Y (cm) Z X (cm) Y (cm) Z
1 0 0 0 0 0 0
2 5,29E-03 1E-04 1,77E-04 4,8E-03 1,2E-04 -1,76E-04
3 1,7E-02 2E-04 2,57E-04 1,59E-02 2,5E-04 -2,5E-04
4 3,03E-02 4E-04 2,39E-04 2,85E-02 3,7E-04 -2,3E-04
5 4,02E-02 5E-05 1,23E-04 3,7E-02 4,96E-04 -1,2E-04
6 4,00E-02 -1,9E-03 0 3,7E-02 -1,8E-03 8,16E-06
7 3,99E-02 -7E-05 -4,2E-05 3,7E-02 -7,2E-05 5E-05
8 3,97E-02 1,7E-03 -1,77E-0,6 3,7E-02 1,7E-03 9,5E-06
9 3,96E-02 -5E-04 1,2E-04 3,7E-02 -4,96E-04 -1,1E-04
10 2,99E-02 -4E-04 2,35E-04 2,8E-02 -3,7E-04 -2,3E-04
11 1,68E-02 -2E-04 2,53E-0,4 1,57E-02 -2,5E-04 -2,5E-04
12 5,23E-03 -1E-04 1,75E-0,4 4,7E-03 -1,2E-04 -1,7E-04
13 0 0 0 0 0 0
110
10 - CONCLUSÃO
Apresentou-se como o método é desenvolvido para estruturas lineares, mas
ele abrange ainda estruturas bidimensionais planas, chapas e cascas e elementos
estruturais tridimensionais, além de solucionar problemas de outros campos de
conhecimento (hidráulica, condução de calor, eletromagnetismo, dispersão de
poluentes, biomecânica, etc.). O método dos elementos finitos foi aprimorado junto
com o desenvolvimento dos computadores, pelo grande número de incógnitas
existentes em problemas complexos a programação ajudou a automatizar os
cálculos, afirmando o sucesso do método.
Atualmente o método dos elementos finitos é o mais utilizado nos programas
cálculo estrutural, mostrando a sua enorme vantagem comparado aos métodos
precedentes. O estudo do método proporciona aos engenheiros lucidez de como
deve ser modelada a estrutura nos softwares com o objetivo de obter os resultados
exatos como esforços, tensões, deformações e deslocamentos, que são utilizados
para dimensionar estruturas de concreto armado, metálicas e madeira com as suas
respectivas normas.
Enfim, o programa apresentado é particularmente útil quando o
resultado buscado são os deslocamentos de um pórtico plano simples. Para as
condições de contorno foram utilizadas as técnicas do número muito grande, que
dão um resultado aproximado do real, provavelmente o SAP 2000 por ser um
software comercial vai utilizar a técnica do zero e o seu resultado consegue ser
exato.
111
11 – REFERÊNCIA BIBLIOGRÁFICA
HALVORSON, M. Visual Basic 2010 Passo a Passo. 2010.
SAVASSI, W. Introdução ao Método dos Elementos Finitos em Análise Linear de
Estruturas. São Carlos, Sp: EESC-USP, 1996. v. 1. 260p.
VAZ, L. E. Método dos Elementos Finitos em Análise de Estruturas. 1. ed. Rio
de Janeiro: Campus Elsevier, 2010. v. 1. 286p.
Pedreiro, M. R. M. EXPLICIT STIFFNESS MATRIX FOR ELEMENT TRIANGULAR
PRISMATIC LINEAR, 2011, 20p.
RODRIGUES, R. O. Introdução ao Método dos Elementos Finitos, 1999.
RODRIGUES, R. O. Técnicas Computacionais para Automatização do Método
dos Elementos Finitos, 1999.
LABAKI, J; MESQUITA, E. Notas de Aula de Resistência dos Materiais II. 2009.
MARTHA, L. F. Análise de Estruturas: Conceitos e Métodos Básicos. Rio de
Janeiro: Campus/Elsevier, 2010. v. 1. 524p.
HIBBELER, R. C. Resistência dos Materiais. 7. Ed. São Paulo: Pearson Prentice
Hall, 2010.
Assan, A. E. Resistência dos Materiais, vol. 2, Campinas: Editora Unicamp, 2013.
RIBEIRO, F. L. B. Introdução ao Método dos elementos finitos, 2004.
112
SILVA, S. Introdução ao Método dos Elementos Finitos, 2009, 196p.
AZEVEDO, Álvaro F.M. Método dos Elementos Finitos, 2003, 258p.
Recommended
![ANÁLISE NÃO LINEAR DE TRELIÇAS PELO, [0.2cm] MÉTODO DOS …€¦ · ANÁLISE NÃO LINEAR DE TRELIÇAS PELO MÉTODO DOS ELEMENTOS FINITOS POSICIONAL Dissertação apresentada ao](https://img.document.onl/doc/110x75/607538cf2606ad7a67089951/anlise-nfo-linear-de-trelias-pelo-02cm-mtodo-dos-anlise-nfo-linear.jpg)
