Elementos Finitos, Notas de estudo de Matemática

Baixe Elementos Finitos e outras Notas de estudo em PDF para Matemática, somente na Docsity!

INTRODUÇÃO AO MÉTODO DOS

ELEMENTOS FINITOS

Notas de Aula do Prof. Fernando L. B. Ribeiro

COPPE / UFRJ – Programa de Engenharia Civil

Rio de Janeiro, 30/03/

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

• 1 INTRODUÇÃO ÍNDICE
• 2 PROBLEMA DE VALOR DE CONTORNO UNIDIMENSIONAL
  • 2.1 Formulação Clássica
  • 2.2 Formulação Variacional
  • 2.3 Aproximação por Elementos Finitos
  • 2.4 Condições para convergência do MEF................................................................................
• 3 ELASTICIDADE PLANA
  • 3.1 Introdução..........................................................................................................................
  • 3.2 Formulação Clássica do Problema de Elasticidade Plana.....................................................
  • 3.3 Formulação Variacional
  • 3.4 Princípio dos Trabalhos Virtuais.........................................................................................
  • 3.5 Energia Potencial Total
  • 3.6 Formulação Variacional Discreta........................................................................................
    • 3.6.1 Energia de Deformação
• 4 PROBLEMAS DE ELASTICIDADE TRIDIMENSIONAL.........................................................
  • 4.1 Introdução..........................................................................................................................
  • 4.2 Formulação Clássica do Problema de Elasticidade Tridimensional
  • 4.3 Formulação do MEF
• 5 ELEMENTOS DE BARRA.........................................................................................................
  • 5.1 Barra Submetida a Esforços Axiais.....................................................................................
  • 5.2 Barra Submetida a Esforços de Flexão................................................................................
• 6 PROBLEMAS DE POTENCIAL
  • 6.1 Introdução..........................................................................................................................
  • 6.2 Formulação Clássica
  • 6.3 Formulação Variacional
  • 6.4 Formulação Variacional Discreta........................................................................................
• 7 ELEMENTOS ISOPARAMÉTRICOS
  • 7.1 Integração no Domínio Real
  • 7.2 Mapeamento Isoparamétrico...............................................................................................
    • 7.2.1 Jacobiano da Transformação de Coordenadas.................................................................
  • 7.3 Mapeamento: Generalização...............................................................................................
  • 7.4 Elementos Isoparamétricos de Continuidade C
    • 7.4.1 Elementos Uniaxiais......................................................................................................
    • 7.4.2 Elementos Quadriláteros - Família de Lagrange
    • 7.4.3 Elementos Quadriláteros - Família de Serendipity
    • 7.4.4 Elementos Triangulares
    • 7.4.5 Hexaedros
    • 7.4.6 Tetraedros
  • 7.5 Exercícios Propostos
• 8 INTEGRAÇÃO NUMÉRICA
  • 8.1 Integração Numérica de Gauss
  • 8.2 Regras de Integração para Triângulos e Tetraedros
• 9 ESTIMATIVAS DE ERRO.........................................................................................................
  • 9.1 Estimativas de Erro Globais e Locais..................................................................................
  • 9.2 Taxas de Convergência
• 10 EXEMPLOS NUMÉRICOS........................................................................................................
  • 10.1 Estado Plano de Deformação..............................................................................................
  • 10.2 Elasticidade Tridimensional

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

INTRODUÇÃO

Grande parte dos problemas de engenharia pode ser formulada através dos princípios gerais da Mecânica do Contínuo ([1],[2]). Este ramo da mecânica trata a matéria como sendo um meio contínuo, sem vazios interiores, desconsiderando sua estrutura molecular. O conceito de “continuum” permite a definição do ponto geométrico (de volume igual a zero), por um limite matemático tal como na definição de derivadas no cálculo infinitesimal. Assim, na Mecânica do Contínuo os princípios da física são escritos sob a forma de equações diferenciais. Os efeitos da constituição interna molecular dos materiais são levados em conta de forma macroscópica através das equações constitutivas do material.

A primeira etapa no processo de modelagem computacional de um fenômeno físico consiste na identificação dos fatores que influenciam 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. A segunda etapa do processo consiste em obter a solução do modelo matemático, tarefa esta atribuída aos métodos numéricos. O Método das Diferenças Finitas é um destes métodos, que como o próprio nome sugere, foi criado com a finalidade específica de resolver sistemas de equações diferenciais. Por outro lado, o Método dos Elementos Finitos (MEF) teve suas origens na análise estrutural. Com o surgimento dos primeiros computadores digitais no início da década de 50, os métodos matriciais para a análise estrutural tiveram um grande desenvolvimento. As primeiras aplicações envolviam apenas estruturas reticuladas, mas a crescente demanda por estruturas mais leves, tais como as encontradas na indústria aeronáutica, conduziu ao desenvolvimento de métodos numéricos que pudessem ser utilizados nas análises de problemas mais complexos. Entre os trabalhos pioneiros nesta linha, podem-se citar os trabalhos de Turner [3] e Argyris [4]. Zienkiewicz, em seu histórico artigo “ The Finite Element Method : from Intuition to Generality ” [5], apresenta uma descrição mais detalhada da evolução do MEF nesta fase inicial. Na década de 70 o MEF teve suas aplicações estendidas a problemas de mecânica dos fluidos e, desde então, vem consolidando-se como um método mais geral de solução de equações diferenciais parciais. Todo o embasamento matemático deste método vem da disciplina Análise de Funcionais ([6],[7],[8]).

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

1 PROBLEMA DE VALOR DE CONTORNO UNIDIMENSIONAL

Um problema de valor de contorno consiste em determinar a função que satisfaz a uma determinada equação diferencial em um dado domínio, conhecendo-se a-priori os valores que a função e/ou suas derivadas assumem no contorno do domínio. Esta descrição corresponde à formulação clássica do problema. A solução pode ser obtida analiticamente, quando possível, ou através de métodos numéricos, como o método das diferenças finitas. Uma outra maneira de formular o problema é através de formulações variacionais, envolvendo equações integrais. Pode-se chegar à formulação variacional do problema de várias maneiras, como por exemplo, através dos princípios dos trabalhos virtuais e da energia potencial mínima, em problemas estruturais, e através do método dos resíduos ponderados de um modo mais geral. Dentre os métodos que se aplicam às formulações integrais podem-se citar o método dos elementos finitos e o método dos elementos de contorno, sendo que este último envolve apenas equações integrais de contorno. Apresentam-se a seguir as formulações clássica e variacional de um problema unidimensional de valor de contorno, introduzindo-se os conceitos básicos do MEF.

1.1 Formulação Clássica

Dados f ( x ) e g , determinar u ( x )tal que

2 ( )^0 em [0,1]

2

• f x = dx

d u (1.1)

( 0 )= (condiçãodecontornonatural)

(1)= 0 (condiçãodecontornoessencial)

g dx

du

u (1.2)

1.2 Formulação Variacional

A formulação variacional correspondente ao problema estabelecido em (1.1)-(1.2) pode ser escrita na seguinte forma:

Dados f ( x ) e g , determinar u ( x ) ∈ U ∀ w ( x ) ∈ W,

1 0

1 0

dx = fwdx - g w ( 0 ) dx

dw dx

du (1.3)

sendo U o espaço de funções admissíveis,

1 0

2 ( ) ( 1 ) 0 , dx dx

du U u x u (1.4)

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

− ^1

0

1

0

1 0 2

2 dx dx

dw dx

du w + dx

du w dx = - dx

d u

1 0

1 0

( 1 ) ( 1 ) ( 0 ) ( 0 ) dx = fw dx dx

dw dx

du w + dx

du w + dx

du (1.8)

Introduzindo as condições de contorno obtém-se

dx = fw dx dx

dw dx

du

g w + ∫ ∫

1 0

1 0

verificando-se, portanto, a validade da equação (1.3).

1.3 Aproximação por Elementos Finitos

O Método dos Elementos Finitos (MEF) resolve por aproximação o problema (1.1) na forma variacional (1.3). O domínio é discretizado em elementos, resultando em uma malha com n pontos nodais. São utilizadas aproximações do tipo:

=

n

j

ux Nj x uj 1

ˆ( ) ( ) , u ˆ ( x )∈ U ˆ ( U ˆ ⊂ U ) (1.10)

n

i=

wx i x wi 1

ˆ ( ) ( ) , w ˆ ( x )∈ W ˆ ( W ˆ ⊂ W ) (1.11)

onde N (^) j e ϕ (^) i são funções de interpolação, e u (^) j e wi são coeficientes constantes. A

figura abaixo apresenta a discretização do domínio em ( n − 1 ) elementos lineares. As

funções de interpolação N (^) j (lineares neste caso) são definidas para cada nó j , de modo

que:

para x x i j

para x x N x i

j j (1.12)

− −

−

j j j j

j

j j j j

j

j

x x x x x

x x

x x x x x

x x

N x

1 1

1

1 1

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

Figura 1.2 – Discretização do domínio unidimensional em elementos lineares.

Observações:

1. Os espaços (^) U ˆ e (^) W ˆ possuem base de dimensão finita e estão contidos em U e W , ou seja, a discretização do domínio representa um truncamento na base infinita dos espaços U e W.

2. Em decorrência de (1.12), os coeficientes u (^) j representam os valores nodais de

u ˆ^ ( x ).

Como u ˆ ( x )é por definição uma solução aproximada, a equação (1.1) não é satisfeita

exatamente para u ˆ( x ), gerando um resíduo R ( x ) no domínio:

( ) ( ) [ 0 , 1 ]

2

2 f x Rx em dx

d u

• = (1.14)

A idéia central do MEF é ponderar este resíduo no domínio (Método dos Resíduos Ponderados) usando as funções de ponderação w ˆ^ ( x ):

R w ˆ dx 0 , w ˆ W ˆ

1 0

∫ = ∀^ ∈ (1.15)

Dentre os métodos de ponderação, o Método de Galerkin é aquele em que as funções ϕ i

em (1.11) são consideradas como sendo iguais as funções N (^) i em (1.10):

n

i=

wx Ni x wi 1

Integrando por partes o primeiro termo em (1.15) obtém-se a equação

0

1 0

1 0 2

2  − − +^ = 

∫ +^ ∫ dx ∫ fwdx

dx

dw dx

du w dx

du w dx

du f wdx = dx

d u

resultando em

ˆ ˆ^1

0

1 0

dx fwdx g w dx

dw dx

du

∫ =^ ∫ − (1.17)

x

1

Nj

j x= 1

1

x= 0

n

j -1 j +

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

Observações:

1. A matriz K é simétrica:

1 0

1 0 ji

i j j i ij (^) dx dx = K

dN dx

dN dx = dx

dN dx

dN K

1. A matriz K é esparsa (muitos de seus coeficientes são iguais a zero). De fato, os

suportes das funções de interpolação N (^) i e N (^) j possuem interseção nula, se o nó i não está conectado ao nó j , resultando em coeficientes K (^) ij = 0 :

Figura 1.

1. Os coeficientes K (^) ij podem ser calculados efetuando-se a integral somente no

elemento que conecta os nós i e j :

1 0

j i

x x

i j i j ij (^) dx dx

dN dx

dN dx dx

dN dx

dN K

Figura 1.

1. Os coeficientes da diagonal principal são positivos e maiores que zero:

1 0

2 dx > dx

dN

K ii ∫ i 

1. Os coeficientes da diagonal K (^) ii podem ser calculados efetuando-se a integral

somente nos elementos que contribuem para o nó i :

dx dx

dN dx + dx

dN dx = dx

dN dx dx

dN K

i i

i i

i i

x x

n x i x

x x

m i i i ii

1 2 2 2 0

(^2 )

− −^ 

Nj

j

Ni

i

Nj

j

Ni

i

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

n ii

m K (^) ii = Kii + K

Figura 1.

1.4 Condições para convergência do MEF

Seja u ˆ( x ) uma solução aproximada do problema (1.3), e u ( x ) a solução exata. Define-

se com erro e ( x )da solução aproximada a diferença

e ( x )= u ( x )− u ˆ( x ) (1.26)

Supondo que u ( x ) seja continuamente diferenciável em todo o domínio ( u ( x ) ∈ C ∞),

uma série de Taylor pode ser empregada para representá-la:

(^200)

2 2

0 0

= = = p x =

p p

x x x dx

d u p

x dx

xd u dx

du u x u x (1.27)

Figura 1.

Portanto, se a solução aproximada for constituída por polinômios de grau p , pode-se

esperar um erro da ordem de ∆ x p +^1. Considerando o intervalo ∆ x como sendo equivalente ao tamanho h dos elementos,

e ≈ O ( ∆ xp +^1 ) (1.28)

Ni

i- 1 i i+ 1

m n

n Ni

i- 1 i i+ 1

m n

m Ni =

x = (^0) x = ∆x

u (0)

u (∆ x )

x

u ( x )

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

Resumindo, são duas as condições para convergência de uma aproximação:

• Condição de completidade: a aproximação u ˆ( x ) deve ser representada por polinômios completos de grau n.
• Condição de compatibilidade: as derivadas da aproximação de ordem n - devem ser contínuas nas interfaces dos elementos (continuidade C n −^1 ).

Observações:

1. Uma função f ( x )definida em um domínio Ω pertence ao espaço vetorial L 2 de

funções de quadrado integrável se a integral de seu quadrado é limitada:

∫ Ω<^ ∞

Ω

f x d

2 ( ) ⇒ f ( x )∈ L 2

Exercícios proposto:

1. Resolver o problema (1.1) para f ( x )=− 2 e g =− 2 , utilizando malhas de 1, 2 e 3

elementos lineares igualmente espaçados. Comparar os resultados com a solução exata u ( x )= x^2 − 2 x + 1.

2 2 0 em [0,1]

2 − = dx

d u

u ( 1 )= 0

dx

du

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

2 ELASTICIDADE PLANA

2.1 Introdução

Os problemas de elasticidade plana se dividem em dois grupos: problemas de estado plano de tensões e problemas de estado plano de deformações. Os problemas de estado plano de tensões são caracterizados por estruturas na forma de chapas planas carregadas no próprio plano, sendo o carregamento uniforme ao longo da espessura. Neste caso, o estado de tensões é completamente definido pelas componentes de tensões σ (^) x , σ (^) y e

τ (^) xy , constantes ao longo da espessura. As demais componentes, σ (^) z , τ (^) yz e τ (^) xz , são

iguais a zero.

Figura 2.1 – Estado plano de tensões.

Por outro lado, os problemas de estado plano de deformações caracterizam-se por estruturas nas quais a dimensão na direção z é muito maior que as dimensões no plano x-y. As cargas são paralelas ao plano x-y e não variam na direção z. Assume-se que os deslocamentos na direção z sejam restringidos. Desta forma, qualquer seção transversal (paralela ao plano x-y ) encontra-se submetida ao mesmo estado de deformações, onde as deformações ε (^) z , γ (^) yz e γ (^) xz são iguais a zero e a tensão normal σ (^) z pode ser obtida em

função das tensões normais σ (^) x e σ (^) y. Portanto, para efeito de análise, basta considerar

uma faixa de espessura unitária compreendida entre duas seções transversais.

Figura 2.2 – Estado plano de deformações.

z

y

x

y

z

y x

y

1

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

Para problemas de estado plano de tensões a matriz D é igual a

−ν

ν

ν

−ν

E

D (2.8)

e para o estado plano de deformações, a matriz D assume a forma

− ν −ν

ν −ν

ν −ν

+ν − ν

−ν

0 0 1 2 21

E 1

D (2.9)

onde E e ν são o módulo de elasticidade e o coeficiente de Poisson, respectivamente.

2.3 Formulação Variacional

Dados b , t e u , determinar u ∈ U ∀ w ∈ W ,

∫ (^ )^ Ω^ =∫ Ω+∫ Γ

Ω Ω Γ

d d d q

L w t^ D L u b t w t^ t w (2.10)

= =( , ) = emΓ, , L 2 x y

u (^) x uy u

u u U u u u (2.11)

= =( , ) = emΓ, , L 2 x y

w (^) x wy u

w w W w w 0 (2.12)

Verificação :

∫^ (^ + )^ Ω=^0

Ω

d t t

L σ b w , ∀ w ∈ W (2.13)

 Ω=^0

∂σ

∂

∂τ  + 

∂τ

∂

∂σ

Ω

b w d x y

b w x y y y

xy y x x

x xy (2.14)

Integrando por partes,

Ω= σ Γ− σ ∂

∂σ

Ω Γ Ω

d x

w wd nw d x

x x x x x x

x (^) (2.15)

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

Ω= τ Γ− τ ∂

∂τ

Ω Γ Ω

d y

w w d nw d y

x x xy y x xy

xy (^) (2.16)

Ω= τ Γ− τ ∂

∂τ

Ω Γ Ω

d x

w w d nw d x

y y xy x y xy

xy (^) (2.17)

Ω= σ Γ− σ ∂

∂σ

Ω Γ Ω

d y

w w d nw d y

y y y y y y

y (^) (2.18)

e substituindo,

=

 Ω+ + Ω^ =

σ +σ ∂

+ τ 

σ +τ ∂

− σ

= σ +τ +τ +σ Γ+

∂σ

∂

∂τ  + 

∂τ

∂

∂σ

Ω Ω

Γ Γ+Γ

Ω

d bw bw d y

w x

w y

w x

w

n n w n n w d

b w d x y

b w x y

x x y y

y y

y xy

x xy

x x

x x xy y x xy x y y y

y y

xy y x x

x xy

u q

A expressão acima pode ser escrita de forma matricial:

∫ (^ +^ )^ Ω= ∫(^ )^ Γ−∫(^ )^ Ω+∫ Ω=^0

Ω Γ=Γ∪Γ Ω Ω

d t d t d t d t t

u q

L σ b w Tn w L w σ b w (2.20)

Introduzindo as condições de contorno, obtém-se a igualdade,

∫ (^ )^ Ω^ =∫ Ω+∫ Γ

Ω Ω Γ

d d d q

L w t^ σ b t w t^ t w (2.21)

e finalmente, substituindo as tensões pelas relações constitutivas chega-se à equação (2.10):

∫ (^ )^ Ω^ =∫ Ω+∫ Γ

Ω Ω Γ

d d d q

L w t^ D L u b t w t^ t w

Pode-se também chegar a esta mesma forma variacional empregando-se o princípio dos trabalhos virtuais ou o princípio da energia potencial total mínima, como será visto a seguir.

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

A energia potencial total do corpo, que é função da configuração deformada u , será portanto igual a:

Π ( u ) = U + W (2.27)

A expressão acima corresponde ao funcional de energia Π( u ), associado às equações

diferenciais de equilíbrio (2.1). O Princípio da Energia Potencial Total tem o seguinte enunciado:

“Seja um corpo impedido de se deslocar como corpo rígido e submetido a forças externas. Dentre todas as configurações deformadas possíveis (que atendem às condições de contorno), aquela que corresponde à configuração de equilíbrio minimiza o funcional de energia potencial total”

Isto significa que para a configuração de equilíbrio, a primeira variação do funcional de energia deve ser igual a zero:

δΠ ( u ) =δ U +δ W = 0 (2.28)

As variações das energias interna e externa são

Ω

U ( δ ε ) t^ Dε d (2.29)

Ω Γ

W d d q

b t^ u t^ t u (2.30)

Somando estas duas equações e igualando a zero obtém-se a mesma expressão de (2.24).

Como a solução exata do problema representa um mínimo absoluto do funcional de energia, qualquer aproximação por deslocamentos superestima a energia potencial total. A medida em que se refina a aproximação, introduzindo-se novos graus de liberdade, mais próximo se chega do valor mínimo de Π.

Pelo princípio da conservação de energia, quando um corpo se deforma sob a influência de forças externas aplicadas lentamente (variando uniformemente a partir de zero), o trabalho realizado pelas forças externas é igual à variação da energia de deformação. Como este trabalho é igual a − W / 2 , pode-se escrever:

U + W / 2 = 0 (2.31)

Substituindo este resultado em (2.27):

Π ( u ) = U + W =− U (2.32)

Portanto, pode-se deduzir da expressão acima que, sendo a energia potencial total Π superestimada, a energia de deformação U será sempre subestimada na formulação em deslocamentos do método dos elementos finitos, na ausência de tensões ou deformações iniciais.

COPPE / UFRJ – Notas de aula do Prof. Fernando L. B. Ribeiro

É importante ressaltar que nem sempre é possível estabelecer um funcional de energia associado às equações diferenciais do problema. Isto só ocorre quando o operador diferencial é auto-adjunto, como é o caso de problemas de elasticidade. Por este motivo, a forma variacional obtida a partir da ponderação de resíduos ponderados é mais geral, e pode ser aplicada a qualquer tipo de problema, como por exemplo, problemas de mecânica dos fluidos.

2.6 Formulação Variacional Discreta

Para obter a formulação variacional discreta deve-se proceder da mesma forma que no item 1.3, com a diferença que agora a solução da equação diferencial u =( ux , uy )é uma

função vetorial de duas componentes, e portanto deve ser aproximada por uma função vetorial:

=

j y

j x

n

j j

j y

x u

u N

N

u

u

1 0

; (^) j

n

j

u ∑ N j u

=

1

Da mesma forma, as funções de ponderação também são vetoriais, com duas componentes:

=

i y

i x

n

i i

i y

x w

w N

N

w

w

1 0

; (^) i

n

i

w ∑ N i w

=

1

Substituindo estas aproximações em (2.10) obtém-se:

∫ (^ )^ Ω^ =∫ Ω+∫ Γ

Ω Ω Γ

d d d q

L w ˆ^ t^ D L u ˆ b t w ˆ t^ t w ˆ (2.35)

Desenvolvendo,

 Ω= Ω+^ Γ

∫ ∑^ ∑ ∫ ∑ ∫ ∑

Ω = = Ω = Γ =

d d d q

i

n

i

i

t i

n

i

i

t j

n

j

j

t i

n

i

N i w D Nu b Nw t Nw 1 1 1 1

L L (2.36)

Fazendo (^)  

w (^) i e (^)  

w (^) i , 0 ( )

j  j ≠ i 

w = para i = 1 ,..., n , obtém-se um sistema

de 2 n equações e 2 n incógnitas:

∑ ∫^ (^ )^ (^ )^ Ω^ =∫ Ω+∫ Γ

= Ω Ω Γ

d d d q

j j i i

t i

n

j

L N D L N u Nb N t 1

, ( i = 1 ,..., n ) (2.37)