Método dos Elementos Finitos2014/15

(www2.dec.fct.unl.pt/seccoes/S_Estruturas/Elementos_finitos)

Corneliu Cismaşiu

Departamento de Engenharia CivilFaculdade de Ciências e Tecnologia

Universidade Nova de Lisboa

cornel@fct.unl.pt

MEF 2014/15 (DEC/FCT/UNL) 1 / 183

Método dos Elementos Finitos 2014/15

Programa

Apresentação

Introdução: Mecânica computacional; Métodos dediscretização; MEF; Problemas físicos/Modelos Matemáticos

Introdução ao MEF: Modelos matemáticos discretos;Modelos matemáticos contínuos: Formulação diferencial,Formulação em resíduos ponderados, Diferenças finitas

Método dos Elementos Finitos: Sistema governativo;Método dos deslocamentos; Tipos de elementos finitos:elemento de barra, elemento plano; Uso de programas deelementos finitos (GiD + Calsef); Erros na análise;Convergência da solução

Aplicações: Barras, estado plano de tensão, lajes

MEF 2014/15 (DEC/FCT/UNL) 2 / 183

Método dos Elementos Finitos 2014/15

Bibliografia

[1] K. J. Bathe. Finite element procedures. Prentice-Hall, 1996.(COTA: TA347.BAT)

[2] O. C. Zienkiewicz and R. L. Taylor and J. Z. Zhu. The finite element method. Volume 1: ItsBasis and Fundamentals. Butterworth-Heinemann, 2005.(COTA: TA640.2.ZIE)

[3] J. N. Reddy. An Introduction to the Finite Element Method. McGraw Hill, 1993.(COTA: TA347.RED)

[4] A. F. M. Azevedo. Método dos elementos finitos. FEUP, 2003.(http://www.alvaroazevedo.com/publications/books/Livro_MEF_AA_1ed)

MEF 2014/15 (DEC/FCT/UNL) 3 / 183

Método dos Elementos Finitos 2014/15

Avaliação de Conhecimentos - contínua sem exame ou trabalho final

Haverão ao longo do período lectivo trabalhos de grupo (2 elementos) aos quais será atribuídaclassificação. Ao fim do semestre, os alunos serão convocados para realizar uma prova oral emque terão de defender estes trabalhos, ao fim de obter a sua classificação final.

É ainda exigido que o número de faltas não justificadas às aulas não exceda um terço do númerototal de aulas leccionadas à respectiva turma. O não cumprimento desta condição implica areprovação na disciplina.

MEF 2014/15 (DEC/FCT/UNL) 4 / 183

Método dos Elementos Finitos 2014/15

Estrutura dos trabalhosDeve seguir a estrutura de uma artigo científico: Título; Autor(es); Resumo; Palavras-chave;Conteúdo (introdução, desenvolvimento e conclusão); Referências bibliográficas.

Assim NÃO!

MEF 2014/15 (DEC/FCT/UNL) 5 / 183

Método dos Elementos Finitos 2014/15

Programa de cálculo automático

CompassFEM-8.0.2R1

MEF 2014/15 (DEC/FCT/UNL) 6 / 183

Introdução – Áreas da Mecânica

Mecânica teórica - Estuda as Leis Fundamentais e os Princípios daMecânica

Mecânica aplicada - Transfere o conhecimento teórico à construção demodelos matemáticos de fenómenos físicos da área de ciênciae da engenharia

Mecânica computacional - Resolve problemas específicos através dasimulação utilizando métodos numéricos implementados emcomputadores digitais

MEF 2014/15 (DEC/FCT/UNL) 7 / 183

Introdução – Áreas da Mecânica

Mecânica computacional - procura soluções a dados problemas

Mecânica aplicada - procura os problemas que admitem dadas soluções

Mecânica teórica - prova a existência de problemas e soluções

MEF 2014/15 (DEC/FCT/UNL) 8 / 183

Áreas da Mecânica – Mecânica computacional

Nano-mecânica e Micro-mecânica - Fenómenos a escala molecular eatómica - concepção de novos materiais e demicro-dispositivos

Mecânica dos meios contínuos - Estudo de corpos a escala macroscópicautilizando modelos contínuos em que a micro-estrutura estáhomogeneizada

Mecânica dos sólidos e das estruturasMecânica dos fluidosMecânica dos sistemas multi-fásicos

Sistemas funcionais - Estudo de mecanismos estruturais, mecânicos,bio-mecânicos, etc.

MEF 2014/15 (DEC/FCT/UNL) 9 / 183

Mecânica computacional – Métodos dediscretização

Converter um modelo matemático contínuo num modelo discreto com umnúmero limitado de graus de liberdade.

Método dos Elementos Finitos (Finite Element Method)

Método dos Elementos de Fronteira (Boundary Element Method)

Método das Diferenças Finitas (Finite Difference Method)

Método dos Volumes Finitos (Finite Volume Method)

Método Espectral (Spectral Method)

Método sem Malha (Mesh-Free Method)

MEF 2014/15 (DEC/FCT/UNL) 10 / 183

Métodos de discretização – Método dosElementos Finitos

Formulação do MEF

DeslocamentosEsforços

Mista - u e σ no domínio, σ ou u na fronteiraHíbrida - u ou σ no domínio, σ ou u na fronteira

Solução do MEF

RigidezFlexibilidade

Mista

MEF 2014/15 (DEC/FCT/UNL) 11 / 183

Métodos de discretização – Método dosElementos Finitos

Obter a solução de um problema de engenharia utilizando o MEF significaessencialmente construir e resolver um sistema governativo de equaçõesalgébricas.

Desenvolvimento dos computadores digitais

⇓

MEF eficiente e confiável

⇓

Grande desenvolvimento do MEF para aplicações práticas de engenharia

MEF 2014/15 (DEC/FCT/UNL) 12 / 183

Método dos Elementos Finitos – História

Origem do MEF - três grupos de investigação nas áreas de:

Matemática Aplicada - R. Courant (1952, 1953)

Física - J. L. Synge (1957)

Engenharia - J. H. Argyris e S. Kelsey (1954, 1955)

Contribuições muito significativas na área de engenharia - J. H. Argyris, S. Kelsey, M. J. Turner,R. W. Clough, H. C. Martin, L. J. Topp, O. C. Zienkiewicz, Y. K. Cheung.

Clough, R. W. “The Finite Element Method in Plane Stress Analysis”, Proceedings, Second ASCEConference on Electronic Computation, Pittsburg, PA, pp. 345-378, Sept. 1960.

MEF 2014/15 (DEC/FCT/UNL) 13 / 183

Análise por Elementos Finitos

Interpretaçao dos resultadosProblema físico

Modelo matemáticoGovernado por equações diferenciais

geometriacinemática

lei do materialcarregamento

condições de fronteiraetc.

Premissas sobre:

Refinamento da malha,parâmetros da solução, etc.

Solução em Elementos FinitosEscolha:

tipo de elementos finitosdensidade da malha

parâmetros da solução

Representação:carregamento

condiçães de fronteira

Avaliação da precisaoda solução em

elementos finitos domodelo matemático

Aperfeiçoamento da concepçãoOptimização estrutural

Refinamento da análise

Aperfeiçoamentodo modelo matemático

Alteração do problemafísico

Solução em Elementos Finitos do Modelo Matemático

A análise por elementos finitos resolve o modelo matemático.

MEF 2014/15 (DEC/FCT/UNL) 14 / 183

Modelos matemáticos

Propriedades do modelo matemático:

Eficácia - permite a obtenção da solução com a precisão desejada aum custo mínimo

Confiabilidade - produz uma solução que se sabe ser contida dentro deuma margem de erro escolhida a priori

MEF 2014/15 (DEC/FCT/UNL) 15 / 183

Problemas físicos/Modelos matemáticos

Problema físico

pinos 2 cmΦ

1000 N

10 cm 28 cm

6 cm 2 cm

4 cm

4 cm

8 cm

2 cm

2 cm

E = 2 × 107 N/cm2 ν = 0.3 t = 0.4 cm

MEF 2014/15 (DEC/FCT/UNL) 16 / 183

Problemas físicos/Modelos matemáticos

Modelo matemático - Teoria de viga

F

L

h

F = 1000 N

L = 28 cm A∗ =5

6A

h = 6 cm t = 0.4 cm

Mmáx = F × L = 28000 Ncm

δmáx =FL3

3EI+

FL

GA∗ ≃ 0.053 cm

MEF 2014/15 (DEC/FCT/UNL) 17 / 183

Problemas físicos/Modelos matemáticos

Modelo matemático - Estado plano de tensão

Mmáx = 29448 Ncm δmáx = 0.070 cm

MEF 2014/15 (DEC/FCT/UNL) 18 / 183

Problemas físicos/Modelos matemáticos

∣

∣

∣

∣

Mestado plano tensão − Mviga

Mestado plano tensão

∣

∣

∣

∣

× 100 ≃ 5%

∣

∣

∣

∣

δestado plano tensão − δvigaδestado plano tensão

∣

∣

∣

∣

× 100 ≃ 24%

O modelo matemático viga é confiável para uma predição do momento flectormáximo com um erro não superior a 5% e do deslocamento máximo comuma precisão de apenas 25% quando comparados com a solução obtida numaanálise elástica linear em estado plano de tensão. Também é eficiente, tendoem conta o esforço computacional necessário.

MEF 2014/15 (DEC/FCT/UNL) 19 / 183

Problemas físicos/Modelos matemáticos

Modelo matemático - Estado plano de tensão

Deformada e o campo das tensões tangencias

MEF 2014/15 (DEC/FCT/UNL) 20 / 183

Análise por Elementos Finitos

- A escolha do modelo matemático deve ser em concordância com a solução queserá prevista.

- O modelo matemático mais eficiente é o que fornece uma resposta confiável como mínimo esforço computacional.

- Uma solução em elementos finitos pode resolver com a precisão desejada apenaso modelo matemático escolhido, podendo prever apenas os fenómeno contidos nomodelo.

- A confiabilidade do modelo matemático tem a ver com a avaliação da precisãoda solução quando comparada com a solução obtida com um modelo matemáticomuito mais complexo.

Na prática → refinamento (p ou/e h) e experiência em engenharia

MEF 2014/15 (DEC/FCT/UNL) 21 / 183

Desastres de Engenharia provocados por erroscomputacionais

Atenção! Os erros computacionais em situações da vida real podem sairmuito caro. . .

“On February 25, 1991, during the Gulf War, an AmericanPatriot Missile battery in Dharan, Saudi Arabia, failed totrack and intercept an incoming Iraqi Scud missile. The Scudstruck an American Army barracks, killing 28 soldiers andinjuring around 100 other people.The Patriot Missile failure,is ultimately attributable to poor handling of rounding errors.”

Excerpted from the report of the General Accounting office GAO/IMTEC-92-26

MEF 2014/15 (DEC/FCT/UNL) 22 / 183

Desastres de Engenharia provocados por erroscomputacionais

“The explosion of the Ariane 5 rocket just after lift-off on its maiden voyage off French Guiana,on June 4, 1996, was ultimately the consequence of simple overflow.”

Excerpted from the report of the Inquiry Board ARIANE 5. Flight 501 Failure

MEF 2014/15 (DEC/FCT/UNL) 23 / 183

Desastres de Engenharia provocados por erroscomputacionais

“The sinking of the Sleipner A offshore platform in Grandsfjorden near Stavanger, Norway, onAugust 23, 1991, resulted in a loss of about 700 million dollars. The post accident investigationtraced the error to inacurrate finite element approximation using the popular finite elementprogram NASTRAN.”

Excerpted from SINTEF, Civil and Environmental Engineering

MEF 2014/15 (DEC/FCT/UNL) 24 / 183

Introdução ao MEF

Análise de um problema de engenharia:

Idealização do problema

Formulação do modelo matemático

Modelos discretos com massa concentradaModelos contínuos

Resolução do modelo matemático

Interpretação dos resultados

MEF 2014/15 (DEC/FCT/UNL) 25 / 183

Introdução ao MEF – Modelos matemáticosdiscretos

O problema pode ser descrito com a precisão desejada por intermédio de umnúmero finito (pequeno) de variáveis.

Idealização do problema

Equilíbrio dos elementos

Assemblagem dos elementos

Cálculo da resposta

Tipo de problemas: estacionários, de propagação, valores e vectores próprios.

MEF 2014/15 (DEC/FCT/UNL) 26 / 183

Modelos matemáticos discretos – Problemasestacionários

Sistema de molas elásticas - solicitação estática

k1

k4

k5

R2

R3R1

u1 u2 u3

k2

k3

MEF 2014/15 (DEC/FCT/UNL) 27 / 183

Modelos matemáticos discretos – Problemasestacionários

Equações de equilíbrio dos elementos

u1

k1 F1(1)

k1 0 00 0 00 0 0

u1

u2

u3

=

F(1)100

u2

F1(2)

F2(2)k2

u1

k2 −k2 0−k2 k2 00 0 0

u1

u2

u3

=

F(2)1

F(2)20

MEF 2014/15 (DEC/FCT/UNL) 28 / 183

Modelos matemáticos discretos – Problemasestacionários

Equações de equilíbrio dos elementos

u2

F1(3)

F2(3)k3

u1

k3 −k3 0−k3 k3 00 0 0

u1

u2

u3

=

F(3)1

F(3)20

u3

F1(4)

F3(4)k4

u1

k4 0 −k4

0 0 0−k4 0 k4

u1

u2

u3

=

F(4)10

F(4)3

MEF 2014/15 (DEC/FCT/UNL) 29 / 183

Modelos matemáticos discretos – Problemasestacionários

Equações de equilíbrio dos elementos

u3

F2(5)

F(5)k5

u2

3

0 0 00 k5 −k5

0 −k5 k5

u1

u2

u3

=

0

F(5)2

F(5)3

MEF 2014/15 (DEC/FCT/UNL) 30 / 183

Modelos matemáticos discretos – Problemasestacionários

Assemblagem dos elementos

F(1)1 + F

(2)1 + F

(3)1 + F

(4)1 = R1

F(2)2 + F

(3)2 + F

(5)2 = R2

F(4)3 + F

(5)3 = R3

k1 + k2 + k3 + k4 −(k2 + k3) −k4

−(k2 + k4) k2 + k3 + k5 −k5

−k4 −k5 k4 + k5

u1

u2

u3

=

R1

R2

R3

K u = R

MEF 2014/15 (DEC/FCT/UNL) 31 / 183

Modelos matemáticos discretos – Problemas depropagação

Sistema de molas elásticas - solicitação dinâmica

k1

k4

k5

u1 u2 u3

R (t)2

R (t)3R (t)1

k2

m3m2m1

k3

MEF 2014/15 (DEC/FCT/UNL) 32 / 183

Modelos matemáticos discretos – Problemas depropagação

Assemblagem dos elementos

F(1)1 + F

(2)1 + F

(3)1 + F

(4)1 = R1(t)− m1u1

F(2)2 + F

(3)2 + F

(5)2 = R2(t)− m2u2

F(4)3 + F

(5)3 = R3(t)− m3u3

Mu + K u = R(t) onde M =

m1 0 00 m2 00 0 m3

Condições iniciais:u(t = 0) = u0 u(t = 0) = u0

MEF 2014/15 (DEC/FCT/UNL) 33 / 183

Modelos matemáticos discretos – Problemas devalores e vectores próprios

Sistema de molas elásticas - vibrações livres

k1

k4

k5

u1 u2 u3

k2

m3m2m1

k3

MEF 2014/15 (DEC/FCT/UNL) 34 / 183

Modelos matemáticos discretos – Problemas devalores e vectores próprios

u = Φ sin(ωt − ψ)

M u + K u = 0

−ω2MΦ sin(ωt − ψ) + KΦ sin(ωt − ψ) = 0

K Φ = ω2M Φ

MEF 2014/15 (DEC/FCT/UNL) 35 / 183

Modelos matemáticos discretos

Exemplos de análises utilizando o “método dos elementos aplicados” (AppliedElement Method)

MEF 2014/15 (DEC/FCT/UNL) 36 / 183

Natureza das soluções

Sistema discreto e condições de carregamento

r2k = kL

2r1

k = kL2

L L

P

F

2F

m/2m

kbarra rígida barra rígidaA B

C

Td

πt

Td

F = sin

t t

F P

MEF 2014/15 (DEC/FCT/UNL) 37 / 183

Natureza das soluções

ku 2

mu1kL

2β

kL2α

kL2β β

α

P

F

2u

2F

1u

2

mu2

cosα ≃ cosβ ≃ 1 sinα ≃ α ≃ u1

Lsinβ ≃ β ≃ u2 − u1

L

∑

Mdt.B = P(u2 − u1) + ku2L +

mu2

2L − FL + kL2 u2 − u1

L= 0

∑

MA = Pu2 + ku22L +mu2

22L − F2L − 2FL + mu1L + kL2 u1

L= 0

MEF 2014/15 (DEC/FCT/UNL) 38 / 183

Natureza das soluções

m 0

0m

2

u1

u2

+

(

5k +2P

L

)

−(

2k +P

L

)

−(

2k +P

L

) (

2k +2P

L

)

u1

u2

=

2F

F

Mu + Ku = F

Frequência próprias do sistema

(K − ω2M)u = 0 ⇒

ω1,2 =

√k

m

√√√√9

2+ 2

P

kL∓

√

33

4+ 8

P

kL+ 2

P2

k2L2T1,2 =

2π

ω1,2

MEF 2014/15 (DEC/FCT/UNL) 39 / 183

Natureza das soluções

Frequências próprias do sistema

ω1

(k/m)1/2

ω2

(k/m)1/2

P/(kL) 0

1

2

3

4

5

6

7

−2 0 2 4 6 8 10

MEF 2014/15 (DEC/FCT/UNL) 40 / 183

Natureza das soluções

Análise estática vs. Análise dinâmica

u2, análise dinâmica

u1, análise dinâmica

u1, análise estática

u2, análise estática

t/Td

m = 1k = 1L = 1P = 1

T = 4 Td 1

−0.2

0.2

0.4

0.6

0.8

1.0

1.2

1.4

0.0 0.2 0.4 0.6 0.8 1.0

0.0

MEF 2014/15 (DEC/FCT/UNL) 41 / 183

Natureza das soluções

Análise estática vs. Análise dinâmica

m = 1k = 1L = 1P = 1

u1, análise estática

u1, análise dinâmica

u2, análise estática

u2, análise dinâmica

T = (T + T )/2d 1 2 t/Td 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

0.0 0.2 0.4 0.6 0.8 1.0

MEF 2014/15 (DEC/FCT/UNL) 42 / 183

Natureza das soluções

Análise estática vs. Análise dinâmica

m = 1k = 1L = 1P = 1

t/Td

u2, análise estática

u1, análise estática

u2, análise dinâmicau1, análise dinâmica

T = T /4d 2

0.0

0.2

0.4

0.6

0.8

1.0

1.2

0.0 0.2 0.4 0.6 0.8 1.0

MEF 2014/15 (DEC/FCT/UNL) 43 / 183

Natureza das soluções

Analisando as frequências próprias do sistema,

ω1,2 =

√k

m

√√√√9

2+ 2

P

kL∓

√

33

4+ 8

P

kL+ 2

P2

k2L2

observa-se que

ω1 = 0 ⇒ Pcr = −2kL o sistema torna-se instável

É importante:

Verificar se o sistema pode tornar-se instável

Decidir se a análise deverá ser estática ou dinâmica, linear ou não-linear

MEF 2014/15 (DEC/FCT/UNL) 44 / 183

Introdução ao MEFModelos matemáticos contínuos

O problema pode ser descrito por intermédio de um conjunto de equaçõesdiferenciais que são válidas no domínio dos elementos.

Juntam-se as condições de fronteira e as condições iniciais (para análisedinâmica).

Formulação diferencial

Formulação variacional

Formulações em resíduos ponderados

Diferenças finitas

MEF

MEF 2014/15 (DEC/FCT/UNL) 45 / 183

Modelos matemáticos contínuosFormulação diferencial

Procura-se a solução analítica de uma equação diferencial de tipo:

A(x , y)∂2u

∂x2+ 2B(x , y)

∂2u

∂x∂y+ C(x , y)

∂2u

∂y2= φ

(

x , y , u,∂u

∂x,∂u

∂y

)

sujeita a determinadas condições de fronteira e condições iniciais.

B2 − AC

< 0 . . . eq. diferencial elíptica (Eq. Laplace)= 0 . . . eq. diferencial parabólica (Eq. do calor)> 0 . . . eq. diferencial hiperbólica (Eq. das ondas)

MEF 2014/15 (DEC/FCT/UNL) 46 / 183

Modelos matemáticos contínuosFormulação diferencial

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

u, x

d

dx

(

EAdu

dx

)

= 0 u|x=0 = 0 EAdu

dx

∣

∣

∣

∣

x=L

= R

MEF 2014/15 (DEC/FCT/UNL) 47 / 183

Modelos matemáticos contínuosFormulação diferencial

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

u, x

Para 0 ≤ x ≤ 100 :

Ed

2u

dx2= 0 ⇒ u(x) =

C0 + C1x

E

u(0) =C0

E= 0 ⇒ C0 = 0

u(x) =C1x

E

MEF 2014/15 (DEC/FCT/UNL) 48 / 183

Modelos matemáticos contínuosFormulação diferencial

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

u, x

Para 100 ≤ x ≤ 180 :

Ed

dx

[(

1 +x − 100

40

)2du

dx

]

= 0 ⇒ u(x) = − C2

E(x − 60)+

C3

E

E

(

1 +x − 100

40

)2du

dx

∣∣∣∣∣x=180

=C2

1600= 100 ⇒ C2 = 160000

MEF 2014/15 (DEC/FCT/UNL) 49 / 183

Modelos matemáticos contínuosFormulação diferencial

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

u, x

O deslocamento e a extensão devem ser contínuos no x = 100 :

u|x=100 ⇒ 100C1

E= −4000

E+

C3

E

du

dx

∣∣∣∣x=100

⇒ C1

E=

100

E

⇒

C1 = 100C3 = 14000

MEF 2014/15 (DEC/FCT/UNL) 50 / 183

Modelos matemáticos contínuosFormulação diferencial

u(x) =

100x

E. . . 0 ≤ x ≤ 100

14000

E− 160000

E(x − 60). . . 100 ≤ x ≤ 180

σ(x) = Edu

dx=

100 . . . 0 ≤ x ≤ 100

160000

(x − 60)2. . . 100 ≤ x ≤ 180

MEF 2014/15 (DEC/FCT/UNL) 51 / 183

Formulação diferencial

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

u, x

25 50 75 100 125 150 175

20

40

60

80

100

120

100E

Factor de escala:

u(x)

x

25

E

50 75 100 125 150 175

0.2

0.4

0.6

0.8

1.0

du(x)

dx

x

Factor de escala: 100

MEF 2014/15 (DEC/FCT/UNL) 52 / 183

Modelos matemáticos contínuosFormulação em resíduos ponderados

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

u, x

d

dx

(

EAdu

dx

)

= 0 u|x=0 = 0 EAdu

dx

∣∣∣∣x=L

= R

u =n∑

i=1

ai fi

fi - funções de forma linearmente independentes

ai - pesos a ser determinados

MEF 2014/15 (DEC/FCT/UNL) 53 / 183

Modelos matemáticos contínuosFormulação em resíduos ponderadosDefine-se o resíduo

Ξ =d

dx

[

EAd

dx

(n∑

i=1

ai fi

)]

Solução exacta ⇒ Ξ = 0

Resíduos ponderados - determinar ai de tal modo que Ξ → 0

Método Galerkin∫

D

fiΞ dD = 0; i = 1, 2, . . . , n

Método dos mínimos quadrados

∂

∂ai

∫

D

Ξ2dD = 0; i = 1, 2, . . . , n

Método da colocação

MEF 2014/15 (DEC/FCT/UNL) 54 / 183

Modelos matemáticos contínuosFormulação em resíduos ponderados

Sejau(x) = a1x + a2x

2 ⇒ f1 = x e f2 = x2

Utilizando o método do Galerkin,

∫ L

0fi

d

dx

(

EAdu

dx

)

dx = 0

Integrando por partes (Teorema da Divergência),

F (x) = u(x) · v(x) dF = (u′v + uv ′)dx

∫ b

a(u′v + uv ′)dx = [uv ]ba ⇒

∫ b

au′v dx = −

∫ b

auv ′

dx + [uv ]ba

−∫ L

0

(

EAdu

dx

)dfi

dxdx +

[

EAdu

dxfi

]L

0

= 0

MEF 2014/15 (DEC/FCT/UNL) 55 / 183

Modelos matemáticos contínuosFormulação em resíduos ponderados

∫ L

0

(

EAdu

dx

)dfi

dxdx = Rfi |x=L i = 1, 2

∫ 100

0(a1+2a2x) dx+

∫ 180

100

[

1+x−100

40

]2

(a1+2a2x) dx =100 × 180

E

∫ 100

0(a1+2a2x)2x dx+

∫ 180

100

[

1+x−100

40

]2

(a1+2a2x)2x dx =100 × 1802

E

1340

3115600

115600102227200

3

a1

a2

=

18000

E

3240000

E

⇒

a1 =128.596

E

a2 =−0.341171

E

MEF 2014/15 (DEC/FCT/UNL) 56 / 183

Formulação em resíduos ponderados

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

u, x

25 50 75 100 125 150 175

20

40

60

80

100

120

100E

Factor de escala:

u(x)

x

Galerkin

solução analítica

25 50 75 100 125 150 175

0.2

0.4

0.6

0.8

1.0

1.2 du(x)

dx

Factor de escala: 100

x

Galerkin

solução analítica E

MEF 2014/15 (DEC/FCT/UNL) 57 / 183

Modelos matemáticos contínuosDiferenças finitas

dw

dx

∣∣∣∣i

wi+1 − wi−1

2h−1 1

d2w

dx2

∣∣∣∣i

wi+1 − 2wi + wi−1

h21 1−2

d3w

dx3

∣∣∣∣i

wi+2 − 2wi+1 + 2wi−1 − wi−2

2h32 −2−1 1

d4w

dx4

∣∣∣∣i

wi+2 − 4wi+1 + 6wi − 4wi−1 + wi−2

h4

−4 −41 16

MEF 2014/15 (DEC/FCT/UNL) 58 / 183

Modelos matemáticos contínuosDiferenças finitas

Operador de Laplace ou Laplaciano

∇2w =∂2w

∂x2+

∂2w

∂y2

∇2w∣∣i,i

−4wi,j + wi+1,j + wi,j+1 + wi−1,j + wi,j−1

h2

−41

1

1

1

MEF 2014/15 (DEC/FCT/UNL) 59 / 183

Modelos matemáticos contínuosDiferenças finitas

Operador biharmónico

∇4w = ∇2(∇2)w =∂4w

∂x4+

∂4w

∂y4+ 2

∂4w

∂x2∂y2

∇4w∣∣i,i

[20wi,j − 8(wi+1,j + wi−1,j + wi,j+1 + wi,j−1)

+2(wi+1,j+1 + wi−1,j+1 + wi−1,j−1 + wi+1,j−1)

+wi+2,j + wi−2,j + wi,j+2 + wi,j−2]/h4

20−8

−8

−8

−8 11

1

1

22

2 2

MEF 2014/15 (DEC/FCT/UNL) 60 / 183

Modelos matemáticos contínuosDiferenças finitas

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

h h

n−1 n+1

u, x

n0 i−1 i i+1

MEF 2014/15 (DEC/FCT/UNL) 61 / 183

Modelos matemáticos contínuosDiferenças finitas

d

dx

(

EAdu

dx

)

= 0 ⇒ E

(dA

dx

du

dx+ A

d2u

dx2

)

= 0

u(0) = 0 EAdu

dx

∣∣∣∣x=L

= R

E

[(Ai+1 − Ai−1

2h

)(ui+1 − ui−1

2h

)

+ Aiui+1 − 2ui + ui−1

h2

]

= 0

u0 = 0 EAnun+1 − un−1

2h= R

MEF 2014/15 (DEC/FCT/UNL) 62 / 183

Modelos matemáticos contínuosDiferenças finitas

h = 20 cm ⇒ n = 9 ⇒

A0 . . .A5 = 1

A6 . . .A10 =(

1 +(i−5)∗h

40

)2

Para i = 1 . . . 9

(Ai+1 + 4Ai − Ai−1)ui+1 − 8Aiui + (−Ai+1 + 4Ai + Ai−1)ui−1 = 0

Condições de fronteira

u0 = 0 u10 − u8 =2Rh

EA9

MEF 2014/15 (DEC/FCT/UNL) 63 / 183

Modelos matemáticos contínuosDiferenças finitas

−8 44 −8 4

4 −8 44 −8 4

2.75 −8 5.256 −18 12

12 −32 2020 −50 30

72 −72

u1

u2

u3

u4

u5

u6

u7

u8

u9

=

00000000F9

F9 = −56000

3E

MEF 2014/15 (DEC/FCT/UNL) 64 / 183

Modelos matemáticos contínuosDiferenças finitas

u =

24.747549.494974.242498.9899123.737136.7

143.182147.071149.663

100

E

du

dx=

1.237371.237371.237371.237370.9427610.4861110.2592590.1620370.111111

100

E

MEF 2014/15 (DEC/FCT/UNL) 65 / 183

Diferenças finitas

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

u, x

25 50 75 100 125 150 175

20

40

60

80

100

120

140

100E

Factor de escala:

u(x)

x

solução analítica

diferenças finitas

25 50 75 100 125 150 175

0.2

0.4

0.6

0.8

1.0

1.2 du(x)

dx

Factor de escala: 100

solução analítica

diferenças finitas

E

MEF 2014/15 (DEC/FCT/UNL) 66 / 183

Modelos matemáticos contínuos

Diferenças finitas

/* Exemplo de programação em Maxima (1/3) */

R:100; /* Força aplicada */

E: 1; /* Módulo de elasticidade */

h:20; /* Dimensão da malha */

Np:round(180/h); /* Número de pontos */

/* Inicialização da área da secção nos pontos da malha */

for i:0 thru Np+1 step 1 do

if i*h <= 100 then a[i]:1 else a[i]:(1+(i*h-100)/40)^2;

/* Inicialização da matriz */

for i:1 thru Np step 1 do for j:1 thru Np step 1 do mat[i,j]:=0;

for i:1 thru Np step 1 do mat[i,i]:-8*a[i];

for i:1 thru Np-1 step 1 do mat[i,i+1]:a[i+1]+4*a[i]-a[i-1];

for i:2 thru Np-1 step 1 do mat[i,i-1]:-a[i+1]+4*a[i]+a[i-1];

mat[Np, Np-1]:8*a[Np];

MAT: genmatrix(mat,Np,Np);

MEF 2014/15 (DEC/FCT/UNL) 67 / 183

Modelos matemáticos contínuos

Diferenças finitas

/* Exemplo de programação em Maxima (2/3) */

/* Inicialização do rhs */

for i:1 thru Np-1 step 1 do rhs[i]:0;

rhs[Np]:-2*R*h/E/a[Np]*(a[Np+1]+4*a[Np]-a[Np-1]);

RHS: makelist(rhs[i],i,1,Np);

/* Solução em deslocamentos */

load ("lapack"); /* Solução alternativa sem LAPACK */

U:dgesv(MAT,RHS); /* U:first(linsolve_by_lu(MAT,RHS,’floatfield) */

sol[Np+1,1]: 2*100*h/E/a[Np]+U[Np-1][1];

sol[0,1]: 0;

for i:1 thru Np step 1 do sol[i,1] : U[i][1];

/* Solução em extensão */

for i:1 thru Np step 1 do sol[i,2] : (sol[i+1,1]-sol[i-1,1])/2/h;

/* Solução em deslocamentos e extensão */

SOL: genmatrix(sol,Np,2);

MEF 2014/15 (DEC/FCT/UNL) 68 / 183

Modelos matemáticos contínuos

Diferenças finitas

/* Exemplo de programação em Maxima (3/3) */

/* Plotar deslocamentos */

xvals: makelist(i, i, 0, Np);

uvals: makelist(sol[i,1], i, 0, Np);

plot2d([discrete, xvals, uvals], [style, points]);

/* Plotar extensão */

x1vals: makelist(i, i, 1, Np);

duvals: makelist(sol[i,2], i, 1, Np);

plot2d([discrete, x1vals, duvals], [style, points]);

MEF 2014/15 (DEC/FCT/UNL) 69 / 183

Método dos elementos finitos

Γu

Γσ

f Γ

f b

uk

uk+1

uk+2

F

elemento finito m

ponto nodal j

j

Conhecendo a geometria do corpo, as cargas aplicadas, as condições de apoio e a lei domaterial,

Detemine os deslocamentos e as respectivas extensões e tensões.

MEF 2014/15 (DEC/FCT/UNL) 70 / 183

Método dos elementos finitos

MEF formulação em deslocamentos – obtida pelo PTV∫

VεTσ dV =

∫

VuT f b

dV +

∫

Γσ

uT f ΓdΓσ +

∑

i

u iTF i

Tensões em equilíbrio com as cargas aplicadasDeformações virtuais correspondentes aos deslocamentos virtuais u

As tensões σ que equilibram as cargas aplicadas assumam-se conhecidas

As deformações virtuais ε são obtidas derivando o vector dos deslocamentos virtuais u

Os deslocamentos virtuais u devem representar um campo continuo, diferenciável.

MEF 2014/15 (DEC/FCT/UNL) 71 / 183

Método dos elementos finitos

Quando a equação do PTV está satisfeita para qualquer deslocamento virtual u, com as tensões σobtidas a partir de um campo de deslocamentos continuo u que satisfaz as condições de fronteiracinemáticas, ficam automaticamente satisfeitas as equações:

Equilíbrio - a equação do PTV é uma equação de equilíbrio

Compatibilidade - o campo de deslocamentos u é continuo e satisfaz as condições de fronteira

Constitutivas - as tensões σ foram calculadas utilizando as relações constitutivas a partir dasdeformações ε que foram obtidas a partir dos deslocamentos u

MEF 2014/15 (DEC/FCT/UNL) 72 / 183

Método dos elementos finitosSistema governativo

Γu

Γσ

f Γ

f b

uk

uk+1

uk+2

F

elemento finito m

ponto nodal j

j

Deslocamentos

u(m)(x , y , z) = H(m)(x , y , z)U

UT = u1 u2 . . . un

H(m) - matriz de interpolação dosdeslocamentos

Deformaçõesε(m)(x , y , z) = E (m)(x , y , z)U

E (m) - matriz de ligação entre as deformações e os deslocamentos

MEF 2014/15 (DEC/FCT/UNL) 73 / 183

Método dos elementos finitosSistema governativo

E (m) = DH(m)

Para o caso 3D, em coordenadas Cartesianas:

D =

∂/∂x 0 00 ∂/∂y 00 0 ∂/∂z

∂/∂y ∂/∂x 00 ∂/∂z ∂/∂y

∂/∂z 0 ∂/∂x

MEF 2014/15 (DEC/FCT/UNL) 74 / 183

Método dos elementos finitosSistema governativo

Tensõesσ

(m)(x , y , z) = k(m)ε(m)(x , y , z)

k(m) - matriz de elasticidadePara o caso 3D, em coordenadas Cartesianas:

k(m) =E

(1 + ν)(1 − 2ν)

1 − ν ν ν 0 0 0ν 1 − ν ν 0 0 0ν ν 1 − ν 0 0 00 0 0 1 − 2ν 0 00 0 0 0 1 − 2ν 00 0 0 0 0 1 − 2ν

MEF 2014/15 (DEC/FCT/UNL) 75 / 183

Método dos elementos finitosSistema governativo

∫

VεTσ dV =

∫

VuT f b

dV +

∫

Γσ

uT f ΓdΓσ +

∑

i

u iTF i

⇓∑

m

∫

V (m)ε(m)T

σ(m)

dV (m) =∑

m

∫

V (m)u(m)T f b(m)

dV (m)+

∑

m

∫

Γ(m)σ

u(m)T f Γ(m)dΓ

(m)σ +

∑

i

u iTF i

Para eficiência podem ser utilizados sistemas de coordenadas diferentes para cada elemento finito.

Admitindo u(m)(x , y , z) = H(m)(x , y , z) U ε(m)(x , y , z) = E (m)(x , y , z) U

MEF 2014/15 (DEC/FCT/UNL) 76 / 183

Método dos elementos finitosSistema governativo

Resulta

UT∑

m

∫

V (m)E (m)TkE (m)

dV (m)U = UT∑

m

∫

V (m)H(m)T f b(m)

dV (m)+

UT∑

m

∫

Γ(m)σ

H(m)T f Γ(m)dΓ

(m)σ + U

TF

K ≡∑

m

∫

V (m)E (m)TkE (m)

dV (m)

︸ ︷︷ ︸

K (m)

- matriz de rigidez global

FB ≡∑

m

∫

V (m)H(m)T f b(m)

dV (m)

︸ ︷︷ ︸

F(m)B

- vector das forças de massa

MEF 2014/15 (DEC/FCT/UNL) 77 / 183

Método dos elementos finitosSistema governativo

F Γ ≡∑

m

∫

Γ(m)σ

H(m)T f Γ(m)dΓ

(m)σ

︸ ︷︷ ︸

F(m)Γ

– vector das forças distribuídas

F - vector das forças concentradas aplicadas nos nós dos elementos

P ≡ FB + F Γ + F – vector das cargas aplicadas

K U = P – equação de equilíbrio

O sistema governativo obtêm-se juntando à equação de equilíbrio as condições de fronteira:[K −e i

−eTi 0

]U

λ

=

P

−uΓ

MEF 2014/15 (DEC/FCT/UNL) 78 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

u1

u2u3 u4

u5u6 u7

u8u9

P

L L

p

EI, EA2EI, 4EA

1 2

Dados numéricos:

EI = 103 kNm2 EA = 0.1EI L = 1 m p = 0.01 kN/m P = 1 kN

MEF 2014/15 (DEC/FCT/UNL) 79 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

u1

u2

u3

u4

u5

u6

Análise de Estruturas I – a matriz de rigidez de uma barra bi-encastrada é dada por:

Ke=

EAL

0 0 −EAL

0 0

0 12EI

L3

6EI

L20 −

12EI

L3

6EI

L2

0 6EI

L2

4EIL

0 −6EI

L2

2EIL

−EAL

0 0 EAL

0 0

0 −12EI

L3−

6EI

L20 12EI

L3−

6EI

L2

0 6EI

L2

2EIL

0 −6EI

L2

4EIL

MEF 2014/15 (DEC/FCT/UNL) 80 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

Elemento 1 u1

u2u3 u4

u5u6

1

K (1) =

EAL

0 0 −EAL

0 0 0 0 00 12EI

L3

6EIL2

0 − 12EIL3

6EIL2

0 0 00 6EI

L2

4EIL

0 − 6EIL2

2EIL

0 0 0−EA

L0 0 EA

L0 0 0 0 0

0 − 12EIL3

− 6EIL2

0 12EIL3

− 6EIL2

0 0 00 6EI

L2

2EIL

0 − 6EIL2

4EIL

0 0 00 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0

MEF 2014/15 (DEC/FCT/UNL) 81 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

Elemento 2 u4

u5u6 u7

u8u9

2

K (2) =

0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 4EA

L0 0 − 4EA

L0 0

0 0 0 0 24EIL3

12EIL2

0 − 24EIL3

12EIL2

0 0 0 0 12EIL2

8EIL

0 − 12EIL2

4EIL

0 0 0 − 4EAL

0 0 4EAL

0 00 0 0 0 − 24EI

L3− 12EI

L20 24EI

L3− 12EI

L2

0 0 0 0 12EIL2

4EIL

0 − 12EIL2

8EIL

MEF 2014/15 (DEC/FCT/UNL) 82 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

A matriz de rigidez global é dada por

K =2∑

i=1

K (i)

K=

EAL

0 0 −EAL

0 0 0 0 0

0 12EI

L3

6EI

L20 −

12EI

L3

6EI

L20 0 0

0 6EI

L2

4EIL

0 −6EI

L2

2EIL

0 0 0

−EAL

0 0 5EAL

0 0 −4EAL

0 0

0 −12EI

L3−

6EI

L20 36EI

L3

6EI

L20 −

24EI

L3

12EI

L2

0 6EI

L2

2EIL

0 6EI

L2

12EIL

0 −12EI

L2

4EIL

0 0 0 −4EAL

0 0 4EAL

0 0

0 0 0 0 −24EI

L3−

12EI

L20 24EI

L3−

12EI

L2

0 0 0 0 12EI

L2

4EIL

0 −12EI

L2

8EIL

MEF 2014/15 (DEC/FCT/UNL) 83 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

O vector das forças aplicadas

PT =

0 −P 0 0 −pL/2 −pL2/12 0 −pL/2 pL2/12

u1

u2u3 u4

u5u6

1

P

u4

u5u6 u7

u8u9

2pL /12 pL /122

2

pL/2 pL/2p

MEF 2014/15 (DEC/FCT/UNL) 84 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

As condições de fronteira uΓ (u7 = u8 = u9 = 0) são impostas utilizando a

matriz e1

eT1 =

0 0 0 0 0 0 1 0 00 0 0 0 0 0 0 1 00 0 0 0 0 0 0 0 1

uΓT =

0 0 0

MEF 2014/15 (DEC/FCT/UNL) 85 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

Sistema governativo

EAL

0 0 - EAL

0 0 0 0 0 0 0 0

0 12EI

L3

6EI

L20 - 12EI

L3

6EI

L20 0 0 0 0 0

0 6EI

L2

4EIL

0 - 6EI

L2

2EIL

0 0 0 0 0 0

- EAL

0 0 5EAL

0 0 - 4EAL

0 0 0 0 0

0 - 12EI

L3- 6EI

L20 36EI

L3

6EI

L20 - 24EI

L3

12EI

L20 0 0

0 6EI

L2

2EIL

0 6EI

L2

12EIL

0 - 12EI

L2

4EIL

0 0 0

0 0 0 - 4EAL

0 0 4EAL

0 0 -1 0 0

0 0 0 0 - 24EI

L3- 12EI

L20 24EI

L3- 12EI

L20 -1 0

0 0 0 0 12EI

L2

4EIL

0 - 12EI

L2

8EIL

0 0 -1

0 0 0 0 0 0 -1 0 0 0 0 00 0 0 0 0 0 0 -1 0 0 0 00 0 0 0 0 0 0 0 -1 0 0 0

u1

u2

u2

u4

u5

u6

u7

u8

u9

λ7

λ8

λ9

=

0-P00

- pL2

- pL2

12

0

- pL2

pL2

12

000

MEF 2014/15 (DEC/FCT/UNL) 86 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

U =

0

−7pL4

48EI−

3PL3

2EIpL3

12EI+ 5PL2

4EI0

−pL4

16EI−

5PL3

12EIpL3

12EI+ 3PL2

4EI000

⇒ U =

0.00000−0.00150

0.001250.00000

−0.000420.000750.000000.000000.00000

[m]

λ =

0pL + P

−pL2

2− 2PL

⇒ λ =

0.0001.010

−2.005

[kN][kN]

[kNm]

MEF 2014/15 (DEC/FCT/UNL) 87 / 183

Método dos elementos finitosExemplo - Método dos deslocamentos

1 kN0.01kN/m

1.010 kN

2.005 kNm

1 m 1 m

0.42 mm

1.5 mm0.00075 rad

0.00125 rad

MEF 2014/15 (DEC/FCT/UNL) 88 / 183

Método dos elementos finitosTipos de elementos finitos

Barras e vigas

Estado plano de tensão/deformação

Estado axissimétrico

Cascas e lajes

Elementos tri-dimensionais

MEF 2014/15 (DEC/FCT/UNL) 89 / 183

Elemento de barra com 2 nós

2 2

22

3

100 cm 80 cm

A B C

A = 1 cmA = (1+y/40) cm

R = 100 N

y

uu1 u2

UT = u1 u2 u3

u1 u2

1 cm2

1 cm2

x100

u2

1 cm2

u3

x80

9 cm2

MEF 2014/15 (DEC/FCT/UNL) 90 / 183

Método dos elementos finitosElemento de barra com 2 nós

Elemento finito 1

u1 u2

1 cm2

1 cm2

x100

k(1) = E

A(1) = 1 cm2

u(1)(x) = H(1)U u(1)(0) = u1 u(1)(100) = u2 ⇒

H(1) =(

1 − x

100

) x

1000

E (1) = DH(1) =

(

− 1

100

)1

1000

MEF 2014/15 (DEC/FCT/UNL) 91 / 183

Método dos elementos finitosElemento de barra com 2 nós

K (1) =

∫

V (1)E (1)Tk(1)E (1)

dV (1)

K (1) = E

∫ 100

0

− 1

100

1

100

0

− 1

100

1

1000

dx

K (1) =E

100

1 −1 0−1 1 00 0 0

MEF 2014/15 (DEC/FCT/UNL) 92 / 183

Método dos elementos finitosElemento de barra com 2 nós

Elemento finito 2

u2

1 cm2

u3

x80

9 cm2

k(2) = E

A(2) =(

1 +x

40

)2cm2

u(2)(x) = H(2)U

u(2)(0) = u2

u(2)(80) = u3

H(2) =

0(

1 − x

80

) x

80

E (2) = DH(2) =

0

(

− 1

80

)1

80

MEF 2014/15 (DEC/FCT/UNL) 93 / 183

Método dos elementos finitosElemento de barra com 2 nós

K (2) =

∫

V (2)E (2)Tk(2)E (2)

dV (2)

K (2) = E

∫ 80

0

0

− 1

80

1

80

0 − 1

80

1

80

(

1 +x

40

)2dx

K (2) =13E

240

0 0 00 1 −10 −1 1

MEF 2014/15 (DEC/FCT/UNL) 94 / 183

Método dos elementos finitosElemento de barra com 2 nós

A matriz de rigidez global

K =2∑

i=1

K (i) =E

240

2.4 −2.4 0−2.4 15.4 −13

0 −13 13

O vector das forças aplicadas

P =

00

100

Condições de fronteirau1 = eT

1 U = 0

MEF 2014/15 (DEC/FCT/UNL) 95 / 183

Método dos elementos finitosElemento de barra com 2 nós

Sistema governativo[

K −e1

−eT1 0

]U

λ

=

P

−uΓ

E/100 −E/100 0 −1−E/100 77E/1200 −13E/240 0

0 −13E/240 13E/240 0−1 0 0 0

u1

u2

u3

λ

=

00

1000

Resulta

u1 = 0 u2 =10000

Eu3 =

154000

13E

λ = −100 . . . a reacção no apoio

MEF 2014/15 (DEC/FCT/UNL) 96 / 183

Método dos elementos finitosElemento de barra com 2 nós

σ(m) = k(m)

ε(m) = k(m)E (m)U

σ(1) = E

−1/100 1/100 0

010000/E

154000/(13E)

= 100

σ(2) = E

0 −1/80 1/80

010000/E

154000/(13E)

= 23.07

MEF 2014/15 (DEC/FCT/UNL) 97 / 183

Elemento de barra com 2 nós

2 2

22

100 cm 80 cm

y

R = 100 N

A B C

A = 1 cmA = (1+y/40) cm

u, x

25 50 75 100 125 150 175

20

40

60

80

100

120

100E

Factor de escala:

u(x)

x

solução analítica

MEF

25 50 75 100 125 150 175

0.2

0.4

0.6

0.8

1.0du(x)

dx

Factor de escala: 100

solução analítica

x

MEF

MEF 2014/15 (DEC/FCT/UNL) 98 / 183

Método dos elementos finitosElemento plano com 3 nós

u1

u2

u3

1

2

3

x

yu4

u6

u5

u =

u(x , y)v(x , y)

= H U

UT = u1 u2 u3 u4 u5 u6

MEF 2014/15 (DEC/FCT/UNL) 99 / 183

Método dos elementos finitosElemento plano com 3 nós

u(x , y) = α1 + α2x + α3y

v(x , y) = β1 + β2x + β3y⇒ u =

u(x , y)v(x , y)

= Φα

Φ =

[

1 x y 0 0 00 0 0 1 x y

]

αT = α1 α2 α3 β1 β2 β3

MEF 2014/15 (DEC/FCT/UNL) 100 / 183

Método dos elementos finitosElemento plano com 3 nós

Explicitando para todos os pontos nodais,

1 x1 y1 0 0 01 x2 y2 0 0 01 x3 y3 0 0 00 0 0 1 x1 y1

0 0 0 1 x2 y2

0 0 0 1 x3 y3

α1

α2

α3

β1

β2

β3

=

u1

u2

u3

u4

u5

u6

Aα = U ⇒ α = A−1 U

Masu = Φα = ΦA−1 U = H U ⇒ H = ΦA−1

MEF 2014/15 (DEC/FCT/UNL) 101 / 183

Método dos elementos finitosElemento plano com 3 nós

1

2

3

4

1

2

3

4

5

4 cm

4 cm

P P

t = 0.1 cm

MEF 2014/15 (DEC/FCT/UNL) 102 / 183

Método dos elementos finitosElemento plano com 3 nós

Elemento 1

u1

u2

u3

u6

u7

u8

1 3

1

2

y

x

A(1) =

1 0 0 0 0 01 2 2 0 0 01 0 4 0 0 00 0 0 1 0 00 0 0 1 2 20 0 0 1 0 4

. . . u1

. . . u3

. . . u2

. . . u6

. . . u8

. . . u7

Φ(1) =

[1 x y 0 0 00 0 0 1 x y

]

MEF 2014/15 (DEC/FCT/UNL) 103 / 183

Método dos elementos finitosElemento plano com 3 nós

H(1) = Φ(1)[

A(1)]−1

H(1) =

4 − x − y

4

x

2

−x + y

40 0 0

0 0 04 − x − y

4

x

2

−x + y

4

D =

∂

∂x0

0∂

∂y

∂

∂y

∂

∂x

MEF 2014/15 (DEC/FCT/UNL) 104 / 183

Método dos elementos finitosElemento plano com 3 nós

E (1) = DH(1) =

−1/4 1/2 −1/4 0 0 00 0 0 −1/4 0 1/4

−1/4 0 1/4 −1/4 1/2 −1/4

K (1) =

∫

A(1)E (1)TkE (1) tdA(1)

Admitindo estado plano de tensão,

k =E

1 − ν2

1 ν 0ν 1 00 0 (1 − ν)/2

MEF 2014/15 (DEC/FCT/UNL) 105 / 183

Método dos elementos finitosElemento plano com 3 nós

K (1) = E (1)TkE (1)t

∫

A(1)dA(1)

Dados os pontos (xi , yi ), i = 0, . . . ,N, com x0 = xN e y0 = yN , a área do polígono plano definidopor estes pontos pode ser calculada (teorema de Green) pela seguinte fórmula:

A =1

2

N−1∑

i=0

(xiyi+1 − xi+1yi )

No caso particular de um triângulo de vertices (a, b), (c, d) e (e, f ), a área é dada por:

A =1

2

∣∣∣∣∣∣

1 1 1a c eb d f

∣∣∣∣∣∣

=1

2[(cf − ed) + (eb − af ) + (ad − cb)]

MEF 2014/15 (DEC/FCT/UNL) 106 / 183

Método dos elementos finitosElemento plano com 3 nós

Resultando, para ν = 0.25,

K (1) = E

11/300 −4/75 1/60 1/60 −1/50 1/300−4/75 8/75 −4/75 −1/75 0 1/751/60 −4/75 11/300 −1/300 1/50 −1/601/60 −1/75 −1/300 11/300 −1/50 −1/60−1/50 0 1/50 −1/50 1/25 −1/501/300 1/75 −1/60 −1/60 −1/50 11/300

MEF 2014/15 (DEC/FCT/UNL) 107 / 183

Método dos elementos finitosElemento plano com 3 nós

Elemento 2

u3

u7

u8

u2

2

u10

u5

3

2 5

x

y

A(2) =

1 0 4 0 0 01 2 2 0 0 01 4 4 0 0 00 0 0 1 0 40 0 0 1 2 20 0 0 1 4 4

. . . u2

. . . u3

. . . u5

. . . u7

. . . u8

. . . u10

Φ(2) =

[1 x y 0 0 00 0 0 1 x y

]

MEF 2014/15 (DEC/FCT/UNL) 108 / 183

Método dos elementos finitosElemento plano com 3 nós

H(2) = Φ(2)[

A(2)]−1

H(2) =

−x + y

4

4 − y

2

−4 + x + y

40 0 0

0 0 0−x + y

4

4 − y

2

−4 + x + y

4

E (2) = DH(2) =

−1/4 0 1/4 0 0 00 0 0 1/4 −1/2 1/4

1/4 −1/2 1/4 −1/4 0 1/4

MEF 2014/15 (DEC/FCT/UNL) 109 / 183

Método dos elementos finitosElemento plano com 3 nós

K (2) =

∫

A(2)E (2)TkE (2) tdA(2)

K (2) = E

11/300 −1/50 −1/60 −1/60 1/75 1/300−1/50 1/25 −1/50 1/50 0 −1/50−1/60 −1/50 11/300 −1/300 −1/75 1/60−1/60 1/50 −1/300 11/300 −4/75 1/601/75 0 −1/75 −4/75 8/75 −4/751/300 −1/50 1/60 1/60 −4/75 11/300

MEF 2014/15 (DEC/FCT/UNL) 110 / 183

Método dos elementos finitosElemento plano com 3 nós

Elemento 3

u3

u8

u10

u5

u9

u4

3

x

y

5

4

3

A(3) =

1 2 2 0 0 01 4 0 0 0 01 4 4 0 0 00 0 0 1 2 20 0 0 1 4 00 0 0 1 4 4

. . . u3

. . . u4

. . . u5

. . . u8

. . . u9

. . . u10

Φ(3) =

[1 x y 0 0 00 0 0 1 x y

]

MEF 2014/15 (DEC/FCT/UNL) 111 / 183

Método dos elementos finitosElemento plano com 3 nós

H(3) = Φ(3)[

A(3)]−1

H(3) =

4 − x

2

x − y

4

−4 + x + y

40 0 0

0 0 04 − x

2

x − y

4

−4 + x + y

4

E (3) = DH(3) =

−1/2 1/4 1/4 0 0 00 0 0 0 −1/4 1/40 −1/4 1/4 −1/2 1/4 1/4

MEF 2014/15 (DEC/FCT/UNL) 112 / 183

Método dos elementos finitosElemento plano com 3 nós

K (3) =

∫

A(3)E (3)TkE (3) tdA(3)

K (3) = E

8/75 −4/75 −4/75 0 1/75 −1/75−4/75 11/300 1/60 1/50 −1/60 −1/300−4/75 1/60 11/300 −1/50 1/300 1/60

0 1/50 −1/50 1/25 −1/50 −1/501/75 −1/60 1/300 −1/50 11/300 −1/60−1/75 −1/300 1/60 −1/50 −1/60 11/300

MEF 2014/15 (DEC/FCT/UNL) 113 / 183

Método dos elementos finitosElemento plano com 3 nós

Elemento 4

u3

u8

u9

u4

u1

u6 4

3

4

y

x1

A(4) =

1 0 0 0 0 01 4 0 0 0 01 2 2 0 0 00 0 0 1 0 00 0 0 1 4 00 0 0 1 2 2

. . . u1

. . . u4

. . . u3

. . . u6

. . . u9

. . . u8

Φ(4) =

[1 x y 0 0 00 0 0 1 x y

]

MEF 2014/15 (DEC/FCT/UNL) 114 / 183

Método dos elementos finitosElemento plano com 3 nós

H(4) = Φ(4)[

A(4)]−1

H(4) =

4 − x − y

4

x − y

4

y

20 0 0

0 0 04 − x − y

4

x − y

4

y

2

E (4) = DH(4) =

−1/4 1/4 0 0 0 00 0 0 −1/4 −1/4 1/2

−1/4 −1/4 1/2 −1/4 1/4 0

MEF 2014/15 (DEC/FCT/UNL) 115 / 183

Método dos elementos finitosElemento plano com 3 nós

K (4) =

∫

A(4)E (4)TkE (4) tdA(4)

K (4) = E

11/300 −1/60 −1/50 1/60 −1/300 −1/75−1/60 11/300 −1/50 1/300 −1/60 1/75−1/50 −1/50 1/25 −1/50 1/50 01/60 1/300 −1/50 11/300 1/60 −4/75

−1/300 −1/60 1/50 1/60 11/300 −4/75−1/75 1/75 0 −4/75 −4/75 8/75

MEF 2014/15 (DEC/FCT/UNL) 116 / 183

Método dos elementos finitosElemento plano com 3 nós

A matriz de rigidez global

K =4∑

i=1

K (i) =

E

11150

160

- 11150

- 160

0 130

1300

- 130

- 1300

0160

11150

- 11150

0 - 160

- 1300

- 130

130

0 1300

- 11150

- 11150

2275

- 11150

- 11150

- 130

130

0 130

- 130

- 160

0 - 11150

11150

160

1300

0 130

- 130

- 1300

0 - 160

- 11150

160

11150

0 - 1300

- 130

1300

130

130

- 1300

- 130

1300

0 11150

- 160

- 11150

160

01

300- 130

130

0 - 1300

- 160

11150

- 11150

0 160

- 130

130

0 130

- 130

- 11150

- 11150

2275

- 11150

- 11150

− 1300

0 130

− 130

1300

160

0 − 11150

11150

− 160

0 1300

− 130

− 1300

130

0 160

− 11150

− 160

11150

MEF 2014/15 (DEC/FCT/UNL) 117 / 183

Método dos elementos finitosElemento plano com 3 nós

O vector das forças aplicadas

PT =

0 0 0 0 0 0 0 0 0 −P

As condições de fronteira uΓ (u1 = u2 = u6 = u7 = 0) são impostas utilizando a matriz e1

eT1 =

1 0 0 0 0 0 0 0 0 00 1 0 0 0 0 0 0 0 00 0 0 0 0 1 0 0 0 00 0 0 0 0 0 1 0 0 0

uΓT =

0 0 0 0

MEF 2014/15 (DEC/FCT/UNL) 118 / 183

Método dos elementos finitosElemento plano com 3 nós

O sistema governativo é dado por

[K −e1

−eT1 0

]U

λ

=

P

−uΓ

resultando

UT =

0 0 - 75116

- 217251392

261251392

0 0 - 27516

- 499251392

- 675251392

P

E

λT =

1 −1 55

873287

P

MEF 2014/15 (DEC/FCT/UNL) 119 / 183

Método dos elementos finitosElemento plano com 3 nós

P

0.63P

0.37PP

P

48.51P/E

MEF 2014/15 (DEC/FCT/UNL) 120 / 183

Método dos elementos finitosElemento plano com 3 nós

O campo das tensões é dado porσ

(m) = k(m)E (m)U

σ(1)=

- 10P29

- 5P58

- 55P16

σ(2)=

65P16

- 3505P1392

- 165P58

σ(3)=

10P29

- 535P174

- 25P16

σ(4)=

- 65P16

- 895P1392

- 125P58

MEF 2014/15 (DEC/FCT/UNL) 121 / 183

Método dos elementos finitosElemento plano com 3 nós

−4P −3P −2P −1P 0 1P 2P 3P 4P

σxyσyσx

4.063P

0.345P−0.345P

−4.063P

−2.518P −2.845P

−0.643P −2.155P

−0.086P 3.075P −1.563P−3.438P

MEF 2014/15 (DEC/FCT/UNL) 122 / 183

Método dos elementos finitosElemento plano com 4 nós

u1

u2

u3

u61

2

4

x

yu5 3

8u

4u u7

u =

u(x , y)v(x , y)

= H U

UT = u1 u2 u3 u4 u5 u6 u7 u8

MEF 2014/15 (DEC/FCT/UNL) 123 / 183

Método dos elementos finitosElemento plano com 4 nós

u(x , y) = α1 + α2x + α3y + α4xy

v(x , y) = β1 + β2x + β3y + β4xy⇒ u =

u(x , y)v(x , y)

= Φα

Φ =

[

1 x y xy 0 0 0 00 0 0 0 1 x y xy

]

αT = α1 α2 α3 α4 β1 β2 β3 β4

MEF 2014/15 (DEC/FCT/UNL) 124 / 183

Método dos elementos finitosElemento plano com 4 nós

Explicitando para todos os pontos nodais,

1 x1 y1 x1y1 0 0 0 01 x2 y2 x2y2 0 0 0 01 x3 y3 x3y3 0 0 0 01 x4 y4 x4y4 0 0 0 00 0 0 0 1 x1 y1 x1y1

0 0 0 0 1 x2 y2 x2y2

0 0 0 0 1 x3 y3 x3y3

0 0 0 0 1 x4 y4 x4y4

α1

α2

α3

α4

β1

β2

β3

β4

=

u1

u2

u3

u4

u5

u6

u7

u8

Aα = U ⇒ α = A−1 UMas

u = Φα = ΦA−1 U = H U ⇒ H = ΦA−1

MEF 2014/15 (DEC/FCT/UNL) 125 / 183

Método dos elementos finitosElemento plano com 4 nós

1

4

2

3

4 cm

4 cm

P P

t = 0.1 cm

MEF 2014/15 (DEC/FCT/UNL) 126 / 183

Método dos elementos finitosElemento plano com 4 nós

Elemento 1

u1u2

u4 u3

u5u6

1

u7u8

1 2

34

x

y

A =

1 0 0 0 0 0 0 01 4 0 0 0 0 0 01 4 4 16 0 0 0 01 0 4 0 0 0 0 00 0 0 0 1 0 0 00 0 0 0 1 4 0 00 0 0 0 1 4 4 160 0 0 0 1 0 4 0

. . . u1

. . . u2

. . . u3

. . . u4

. . . u5

. . . u6

. . . u7

. . . u8

Φ =

[1 x y xy 0 0 0 00 0 0 0 1 x y xy

]

MEF 2014/15 (DEC/FCT/UNL) 127 / 183

Método dos elementos finitosElemento plano com 4 nós

H = Φ [A]−1

H =

[

1- x+y4

+ xy16

x4- xy16

xy16

y4- xy16

0 0 0 0

0 0 0 0 1- x+y4

+ xy16

x4- xy16

xy16

y4- xy16

]

D =

∂

∂x0

0∂

∂y

∂

∂y

∂

∂x

MEF 2014/15 (DEC/FCT/UNL) 128 / 183

Método dos elementos finitosElemento plano com 4 nós

E = DH =

y−416

4−y16

y16

− y16

0 0 0 00 0 0 0 x−4

16− x

16x16

4−x16

x−416

− x16

x16

4−x16

y−416

4−y16

y16

− y16

K =

∫

AETkE tdA

Admitindo estado plano de tensão,

k =E

1 − ν2

1 ν 0ν 1 00 0 (1 − ν)/2

MEF 2014/15 (DEC/FCT/UNL) 129 / 183

Método dos elementos finitosElemento plano com 4 nós

Admitindo ν = 0.25, resulta, por exemplo, para o elemento K11 da matriz de rigidez,

K11 =E

10

∫

A

176 − 24x − 64y + 3x2 + 8y2

1920dA

Como neste caso particular o domínio do elemento é rectangular,

K11 =E

10

∫ 4

0

∫ 4

0

176 − 24x − 64y + 3x2 + 8y2

1920dx dy =

11E

225

Nos programas de elementos finitos, a integração analítica costuma ser substituída por umaintegração numérica, utilizando o Método de Gauss-Legendre.

MEF 2014/15 (DEC/FCT/UNL) 130 / 183

Elemento plano com 4 nósMétodo de Gauss-Legendre

Método de Gauss-Legendre para a integração numérica consiste em aproximar o integral atravésde um somatório.

Para o caso uni-dimensional,∫ 1

−1f (x) dx ≃

n∑

i=1

wi f (xi )

onde xi são as coordenadas dos pontos de integração e wi os pesos associados. A regra deintegração de Gauss de n pontos permite integrar exactamente polinómios de grau inferior ouigual a 2n − 1.

MEF 2014/15 (DEC/FCT/UNL) 131 / 183

Elemento plano com 4 nósMétodo de Gauss-Legendre

Pontos de integração e pesos no intervalo [−1, 1]n xi wi

1 0.00000 00000 00000 2.00000 00000 000002 ±0.57735 02691 89626 1.00000 00000 000003 ±0.77459 66692 41483 0.55555 55555 55555

0.00000 00000 00000 0.88888 88888 88888

MEF 2014/15 (DEC/FCT/UNL) 132 / 183

Elemento plano com 4 nósMétodo de Gauss-Legendre

Para o caso bi-dimensional,

∫ 1

−1

∫ 1

−1f (x , y) dx dy ≃

n∑

i=1

wi

∫ 1

−1f (xi , y) dy ≃

n∑

i=1

m∑

j=1

wiwj f (xi , yj )

x=0.577...

y=−0.577...

y=0.577...

x=−0.577...

y

x

Integração no domínio quadrilateral (x , y ∈[−1, 1]) com 2 × 2 pontos, exacta para umpolinómio de grau 3:

1x y

x2 xy y2

x3 x2y xy2 y3

MEF 2014/15 (DEC/FCT/UNL) 133 / 183

Elemento plano com 4 nósMapeamento isoparamétrico

Quando o domínio de integração não é definido entre -1 e 1 (ou, para o caso 2D, não é rectangular),é habitual proceder a um mapeamento isoparamétrico.

1D, x ∈ Ω → s ∈ [−1, 1]∫

Ωf (x) dx =

∫ 1

−1f [x(s)]J ds

2D, x , y ∈ Ω → s, t ∈ [−1, 1]

∫

Ωf (x , y) dx dy =

∫ 1

−1

∫ 1

−1f [x(s, t), y(s, t)]J ds dt

onde J (Jacobiano) é o determinanto da matriz da transformação (matriz Jacobiano)

MEF 2014/15 (DEC/FCT/UNL) 134 / 183

Método de Gauss-LegendreMapeamento isoparamétrico

x

y

s

t

2

34

1

3 (1,1)4 (−1,1)

2 (1,−1)1 (−1,−1)

Ω

Funções de forma:N1 = (1 − s)(1 − t)/4 N2 = (1 + s)(1 − t)/4

N3 = (1 + s)(1 + t)/4 N4 = (1 − s)(1 + t)/4

MEF 2014/15 (DEC/FCT/UNL) 135 / 183

Método de Gauss-LegendreMapeamento isoparamétrico

As coordenadas x e y podem exprimir-se em função ao s e t,

x =

4∑

i=1

Nixi y =

4∑

i=1

Niyi

A matriz da transformação (matriz Jacobiano) define-se como

J =

∂x

∂s

∂y

∂s

∂x

∂t

∂y

∂t

MEF 2014/15 (DEC/FCT/UNL) 136 / 183

Mapeamento isoparamétricoExemplo 1D

x1 x2

Ω x s

−1 1Funções de forma:

N1 = (1 − s)/2 N2 = (1 + s)/2

Mapeamento isoparamétrico:

x =2∑

i=1

Nixi =x1 + x2 + s(x2 − x1)

2

Jacobiano

J = det(J) =∂x

∂s=

x2 − x1

2Por exemplo,

∫ 7

3x3

dx =

∫ 1

−1

[3 + 7 + s(7 − 3)

2

]3 7 − 3

2ds = 580

MEF 2014/15 (DEC/FCT/UNL) 137 / 183

Mapeamento isoparamétricoExemplo 2D

1 2

3

4Ω

x

y

s

t

(−1,−1) (1,−1)

(1,1)(−1,1)

Coordenadas dos nós:

x1 = 0, y1 = 0 x2 = 1, y2 = 0 x3 = 2, y3 = 2 x4 = 0, y4 = 1

Funções de forma:N1 = (1 − s)(1 − t)/4 N2 = (1 + s)(1 − t)/4

N3 = (1 + s)(1 + t)/4 N4 = (1 − s)(1 + t)/4

MEF 2014/15 (DEC/FCT/UNL) 138 / 183

Mapeamento isoparamétricoExemplo 2D

Mapeamento isoparamétrico:

x =4∑

i=1

Nixi = (1 + s)(3 + t)/4 y =4∑

i=1

Niyi = (3 + s)(1 + t)/4

Matriz Jacobiano:

J = det(J) = det

[(3 + t)/4 (1 + t)/4(1 + s)/4 (3 + s)/4

]

=4 + s + t

8

Por exemplo,

∫

Ωx3

dx dy =

∫ 1

−1

∫ 1

−1

[(1 + s)(3 + t)

4

]3 4 + s + t

8ds dt =

23

10

MEF 2014/15 (DEC/FCT/UNL) 139 / 183

Método dos elementos finitosElemento plano com 4 nós

Regressando,

K11 =E

10

∫ 4

0

∫ 4

0

176 − 24x − 64y + 3x2 + 8y2

1920dx dy

Sendo o as coordenadas do domínio de integração,

(0, 0) (4, 0) (4, 4) (0, 4)

do mapeamento isoparamétrico resulta:

x =4∑

i=1

Nixi = 2(1 + s) y =4∑

i=1

Niyi = 2(1 + t)

J = det

[2 00 2

]

= 4

MEF 2014/15 (DEC/FCT/UNL) 140 / 183

Método dos elementos finitosElemento plano com 4 nós

Com,f (x , y) = (176 − 24x − 64y + 3x2 + 8y2)/1920∫ 4

0

∫ 4

0f (x , y)dx dy =

∫ 1

−1

∫ 1

−1g(s, t)ds dt

ondeg(s, t) = f [x(s, t), y(s, t)]J

Utilizando a regra de integração de Gauss com 2 × 2 pontos,

I =

∫ 1

−1

∫ 1

−1g(s, t)ds dt ≃

2∑

i=1

2∑

j=1

wiwjg(si , ti )

coms1,2 = t1,2 = ±1/

√3 e w1,2 = 1

MEF 2014/15 (DEC/FCT/UNL) 141 / 183

Método dos elementos finitosElemento plano com 4 nós

Resulta

I = g

(−1√3,−1√

3

)

+ g

(−1√3,

1√3

)

+ g

(1√3,−1√

3

)

+ g

(1√3,

1√3

)

=22

45

E o termo da matriz de rigidez,

K11 =E

10

22

45=

11E

225

MEF 2014/15 (DEC/FCT/UNL) 142 / 183

Método dos elementos finitosElemento plano com 4 nós

A matriz de rigidez global

K = E

11225

−13450

−11450

1225

160

−1300

−160

1300

−13450

11225

1225

−11450

1300

−160

−1300

160

−11450

1225

11225

−13450

−160

1300

160

−1300

1225

−11450

−13450

11225

−1300

160

1300

−160

160

1300

−160

−1300

11225

1225

−11450

−13450

−1300

−160

1300

160

1225

11225

−13450

−11450

−160

−1300

160

1300

−11450

−13450

11225

1225

1300

160

−1300

−160

−13450

−11450

1225

11225

O vector das forças aplicadas

PT =

0 0 0 0 0 0 −P 0

MEF 2014/15 (DEC/FCT/UNL) 143 / 183

Método dos elementos finitosElemento plano com 4 nós

As condições de fronteira uΓ (u1 = u4 = u5 = u8 = 0) são impostas utilizando a matriz e1

eT1 =

1 0 0 0 0 0 0 00 0 0 1 0 0 0 00 0 0 0 1 0 0 00 0 0 0 0 0 0 1

uΓT =

0 0 0 0

MEF 2014/15 (DEC/FCT/UNL) 144 / 183

Método dos elementos finitosElemento plano com 4 nós

O sistema governativo é dado por

[K −e1

−eT1 0

]U

λ

=

P

−uΓ

resultando

UT =

0 - 276751474

326251474

0 0 - 28550737

- 38450737

0 P

E

λT =

1 −1 45

672267

P

MEF 2014/15 (DEC/FCT/UNL) 145 / 183

Método dos elementos finitosElemento plano com 4 nós

P

P

0.67P

0.33P

P

52.17P/E

MEF 2014/15 (DEC/FCT/UNL) 146 / 183

Método dos elementos finitosElemento plano com 4 nós

O campo das tensões é dado porσ = kEU

σT =

−15(246+11x−134y)

737−15(123+88x−67y)

14745(−2284+603x−198y)

2948

σx σy σxy

0 1 2 3 40

1

2

3

4

0 1 2 3 40

1

2

3

4

0 1 2 3 40

1

2

3

4

MEF 2014/15 (DEC/FCT/UNL) 147 / 183

Método dos elementos finitosProgramas de elementos finitos

IFER - Internet Finite Element Resources

Domínio público

Código incluído (Freeware): mais de 70 programas: ADVENTURE(ADVanced ENgineering analysis Tool for Ultra large REalworld), FELT, OpenSees (Open System for EarthquakeEngineering Simulation), . . .

Código não incluído (Shareware): mais de 30 programas: CADRE, IMAGINE(Integrated Modelling and Analysis in Geotechnics),PlastFEM, . . .

Programas comerciais : mais de 120 programas: ABAQUS, Adina, Algor, ANSYS, COSMOS,DIANA, PLAXIS, SAP2000, . . .

MEF 2014/15 (DEC/FCT/UNL) 148 / 183

Método dos elementos finitosGiD + Calsef

GiD - pré/pós-processador gráfico, incluindo um gerador de malha, que pode serutilizado em conjunto com uma grande variedade de programas de análisenumérica.

AnáliseComputacional

(CALSEF)

Definição da geometriaDefinição das cargas e das condições de fronteira

Definição dos materiaisGeração da malha

(GiD)

Visualização dosresultados

(GiD)

Calsef - código (elementos finitos) para a análise de sólidos (2D e 3D), lajes e cascas.

International Center for Numerical Methods in Engineering

MEF 2014/15 (DEC/FCT/UNL) 149 / 183

GiD + CalsefExemplo de aplicação

4 cm

4 cm

P

t = 0.1 cm

ν = 0.25

E = 1

P = 1

Definir o tipo do problema: Estado plano

MEF 2014/15 (DEC/FCT/UNL) 150 / 183

GiD + CalsefExemplo de aplicaçãoEscolher um título para o problema, escolher o estado plano de tensão, não considerar o pesopróprio, etc.

MEF 2014/15 (DEC/FCT/UNL) 151 / 183

GiD + CalsefExemplo de aplicação

Verificar as propriedades dos materiais e as unidades do problema.

MEF 2014/15 (DEC/FCT/UNL) 152 / 183

GiD + CalsefExemplo de aplicação

A janela das coordenadas facilita a introdução dos nós.

MEF 2014/15 (DEC/FCT/UNL) 153 / 183

GiD + CalsefExemplo de aplicação

Definir a geometria da estrutura: pontos e linhas.

MEF 2014/15 (DEC/FCT/UNL) 154 / 183

GiD + CalsefExemplo de aplicação

Definir a geometria da estrutura : superfícies.

MEF 2014/15 (DEC/FCT/UNL) 155 / 183

GiD + CalsefExemplo de aplicação

Definir as condições de fronteira: deslocamentos impostos, apoios elásticos e cargas aplicadas.

MEF 2014/15 (DEC/FCT/UNL) 156 / 183

GiD + CalsefExemplo de aplicação

Definir os deslocamentos impostos.

MEF 2014/15 (DEC/FCT/UNL) 157 / 183

GiD + CalsefExemplo de aplicação

Definir as cargas aplicadas.

MEF 2014/15 (DEC/FCT/UNL) 158 / 183

GiD + CalsefExemplo de aplicaçãoDefinir um novo material com as propriedades desejadas.

MEF 2014/15 (DEC/FCT/UNL) 159 / 183

GiD + CalsefExemplo de aplicação

Atribuir o material aos elementos geométricos.

MEF 2014/15 (DEC/FCT/UNL) 160 / 183

GiD + CalsefExemplo de aplicação

Definir o tipo de elementos finitos a utilizar.

MEF 2014/15 (DEC/FCT/UNL) 161 / 183

GiD + CalsefExemplo de aplicação

Gerar a malha de elementos finitos.

MEF 2014/15 (DEC/FCT/UNL) 162 / 183

GiD + CalsefExemplo de aplicaçãoComeçar a análise.

MEF 2014/15 (DEC/FCT/UNL) 163 / 183

GiD + CalsefExemplo de aplicaçãoEntrar na fase de pós-processamento.

MEF 2014/15 (DEC/FCT/UNL) 164 / 183

GiD + CalsefExemplo de aplicação

Campo de deslocamentos

MEF 2014/15 (DEC/FCT/UNL) 165 / 183

GiD + CalsefExemplo de aplicação

Campo de tensõesσx σy σxy

MEF 2014/15 (DEC/FCT/UNL) 166 / 183

Método dos elementos finitosFontes de erros na análise

Discretização - Funções de interpolação da geometria e da solução

Integração numérica - Cálculo dos elementos das matrizes utilizando métodos numéricosde integração

Relações constitutivas - Uso de modelos de material com comportamento não-linear

Solução das eq. de equilíbrio dinâmico - Integração numérica em tempo e/ousobreposição modal

Solução do sistema governativo por métodos iterativos - Gauss-Seidel, método dosgradientes conjugados, Newton-Raphson, . . .

Arredondamentos - Precisão da máquina

MEF 2014/15 (DEC/FCT/UNL) 167 / 183

Fontes de erros na análiseDistorção da malha

a

L L

pδ

L

a/L [m] 0.00 0.25 0.5 0.75 1.00δ/δ0 1.000 0.897 0.749 0.644 0.551

Na prática, sendo usadas malhas com número relativamente reduzido de elementos e sendo poucohabitual fazer estudos extensos sobre a convergência da solução, recomenda-se o uso de malhasnão-distorcidas e/ou de elementos poucos sensíveis a distorção.

MEF 2014/15 (DEC/FCT/UNL) 168 / 183

Método dos elementos finitosConvergência da solução

A convergência da solução numérica implica que, com o refinamento da malha, todas as condiçõescinemáticas, estáticas e constitutivas contidas no modelo matemático utilizado ficam satisfeitas.

Quando os elementos finitos são completos (as funções de aproximação dos deslocamentos podemrepresentar deslocamentos de corpo rígido) e tanto os elementos como a malha são compatíveis(deslocamentos contínuos) a convergência fica monotónica.

Quando possível, recomenda-se o uso da energia de deformação como grandeza para o estudo daconvergência da solução.

U =1

2

∫

Vσ

Tε dV = · · · = 1

2UTKU

MEF 2014/15 (DEC/FCT/UNL) 169 / 183

Método dos elementos finitosConvergência da solução

6L

L p

A

Elem.(g.d.l.) 12 (28) 46 (76) 188 (246) 836 (950) 352∗ (1550)

δAv 0.288 0.631 0.874 0.970 1.000

σAx 0.000 1.802 1.366 1.025 1.000

σAy -0.859 0.248 0.657 0.977 1.000

σAxy 2.066 1.802 1.366 1.025 1.000

∗ elementos triangulares com 6 nós

MEF 2014/15 (DEC/FCT/UNL) 170 / 183

Trabalhos práticos

Trabalho No1: Barras

Trabalho No2: Estado plano de tensão

Trabalho No3: Singularidades

Trabalho No4: Lajes

MEF 2014/15 (DEC/FCT/UNL) 171 / 183

Trabalho No1: Barras

Considere a bara de homogénea de secção variável representada, sujeita às condições de fronteirae ao carregamento indicado. Determine os deslocamentos axiais e o esforço normal ao longo dabarra utilizando uma formulação:

diferencial (solução analítica);

em resíduos ponderados (o método de Galerkin);

diferenças finitas;

elementos finitos (elementos de barra com 2 nós).

Compare e comente as soluções obtidas.

MEF 2014/15 (DEC/FCT/UNL) 172 / 183

Trabalho No1: Barras

22A = (3−y/40) cm22A = (1+x/50) cm

R = 100 N

A B C

yx

80 cm100 cm

u

MEF 2014/15 (DEC/FCT/UNL) 173 / 183

Trabalho No2: Estado plano de tensão

Estime os deslocamentos e o estado de tensão dos pontos A, B, C , D, E e F da consolacurta representada sujeita a uma carga concentrada P na extremidade. Considere uma espessuraconstante t = 0.1 m, um módulo de elasticidade constante e um coeficiente de Poisson ν = 0.3.Admite o estado plano de tensão e utilize:

uma malha de elementos finitos triangulares com 3 nós;

uma malha de elementos finitos rectangulares com 4 nós;

Apresente a deformada da estrutura para os dois casos considerados. Compare e comente assoluções obtidas.

MEF 2014/15 (DEC/FCT/UNL) 174 / 183

Trabalho No2: Estado plano de tensão

3 m5 m

A B

F

E

D

C

1 m

3 m

2 m

P

MEF 2014/15 (DEC/FCT/UNL) 175 / 183

Trabalho No3: Singularidades

Considere as seguintes placas planas homogéneas, de espessura e módulo de elasticidadeconstantes, unitários e coeficiente de Poisson ν = 0.1, sujeitas ao carregamento indicado.Admitindo o estado plano de tensão, utilize os programas Gid + Calsef para estudar, em cada umdos casos, a convergência da solução. Apresente as malhas de elementos finitos utilizadas. Emseguida, analise cada uma das placas utilizando a melhor malha para deteminar:

o estado de tensão no ponto P;

a distribuição de tensões normais e tangenciais ao longo da fronteira estática e cinemáticada placa;

o campo de deslocamentos.

MEF 2014/15 (DEC/FCT/UNL) 176 / 183

Trabalho No3: Singularidades

P0.1

0.1

q=1

1

1

1

1

Q=1

1P 0.3

1.51.5

MEF 2014/15 (DEC/FCT/UNL) 177 / 183

Trabalho No4: Lajes

O comportamento estrutural das lajes depende dos seguintes factores:

- tipo de apoios e cargas (condições de fronteira)

- relação entre os vãos (condiciona a direcção de flexão dominante)

- comportamento mecânico do material

- relação da espessura com o menor dos vãos (condiciona o tipo do modelo de análise)

lajes finas (Kirchhoff) - espessura/vão ≤ 1/5 (1/10), deslocamentotransversal máximo/espessura ≤ 1/5 (1/10)lajes espessas (Reissner-Mindlin)

MEF 2014/15 (DEC/FCT/UNL) 178 / 183

Trabalho No4: Lajes

Hipóteses simplificativas:

- linearidade física - material com comportamento elástico linear

- linearidade geométrica - pequenos deslocamentos e pequenas deformações

- homogeneidade e isotropia do material

Admita-se ainda que:

- fibras rectas normais ao plano médio da laje permanecem rectas após a deformação

- fibras rectas normais ao plano médio da laje são inextensíveis

- fibras rectas normais ao plano médio da laje permanecem rectas após a deformação eperpendiculares ao plano médio - Kirchhoff

MEF 2014/15 (DEC/FCT/UNL) 179 / 183

Trabalho No4: Lajes

Laje de Kirchhoff

∇4w(x , y) =q(x , y)

Df

Df =Eh3

12(1 − ν2)

Df - rigidez à flexão do elemento de laje

θ (x,y)x

θ (x,y)y

x

y

z

w(x,y)

Campos de deslocamentos

vy

vx

vx

vy

mx

mx

mxy

mxy

mxy

my

mxy

my

x

y

Campos de esforços

MEF 2014/15 (DEC/FCT/UNL) 180 / 183

Trabalho No4: Lajes

Considere a seguinte laje em betão armado, simplesmente apoiada nos pilares e nas paredesresistentes representadas. Admitindo que o único carregamento é o peso próprio, apresente:

- o campo dos deslocamentos;

- os campos dos esforços (Mx , My , Mxy , Vx , Vy );

- os diagramas de momentos segundo as linhas de corte AA, BB e EE .

Repita a análise admitindo na análise que a secção transversal dos pilares e das paredes édesprezável. Compare e comente os resultados obtidos nas duas análises.

MEF 2014/15 (DEC/FCT/UNL) 181 / 183

Trabalho No4: Lajes

1.00 6.00 6.00 6.00 1.00

P2(.40 x .40)

P3(.60 x .40)

P5(.60 x .40)

P7(.60 x .40)

P10(.60 x .40)

P9(.40 x .60)

P8(.40 x .60)

P6(.40 x .60)

P4(.40 x .60)

P1(.40 x .40)

.40.40 3.80

2.60

.40

.40

1.00

3.00

6.00

0.60

3.60

.405.

00.2

0

5.00

.20

e=0.22 m

A

B

D E F G

C C

B

A

D E F G

Betão Armado (Unidades: N,m,rad)

Carga = peso específico = 25000Poisson = 0.2Módulo de elasticidade = 3.00E10

MEF 2014/15 (DEC/FCT/UNL) 182 / 183

Anexo: FAQ - Maxima

MEF 2014/15 (DEC/FCT/UNL) 183 / 183