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.