Apostila do Método dos Elementos Finitos

Profa. Claudia Mazza Dias

March 31, 2015

Contents

1 Métodos Numéricos 4
1.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 O Cálculo Variacional . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.1 Método Aproximado de Rayleight-Ritz . . . . . . . . . . . 11
1.2.2 O Método de Galerkin . . . . . . . . . . . . . . . . . . . . 13

2 Método dos Elementos Finitos 15

2.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 O Problema de Valor no Contorno . . . . . . . . . . . . . . . . . 15
2.3 Formulação Clássica - Forma Forte (S) . . . . . . . . . . . . . . 16
2.4 Formulação Variacional - Forma Fraca (W) . . . . . . . . . . . . 17
2.5 Formulação de Galerkin (G) . . . . . . . . . . . . . . . . . . . . . 21
2.6 Formulação Matricial (M) . . . . . . . . . . . . . . . . . . . . . . 22
2.7 A discretização e os Espaços de Elementos Finitos . . . . . . . 25

3 O Elemento Finito Unidimensional Linear 31

3.1 Propriedades Globais . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Propriedades Locais . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.1 Montagem da Matriz de Rigidez Global e do Vetor de
Forças Global . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2.2 Cálculo Explı́cito da Matriz de Rigidez: . . . . . . . . . . . 36
3.2.3 Cálculo Explı́cito do Vetor de Forças: . . . . . . . . . . . 38

4 O PVC Bi-dimensional 39
4.1 Exemplo: Equação do Calor: . . . . . . . . . . . . . . . . . . . . 40

5 Elementos Isoparamétricos: 46
5.1 Elemento Unidimensional Linear: . . . . . . . . . . . . . . . . . . 47
5.2 Elemento Triangular Linear: . . . . . . . . . . . . . . . . . . . . . 49
5.3 Elemento Quadrilátero Bilinear: . . . . . . . . . . . . . . . . . . . 50
5.4 Integração Numérica: . . . . . . . . . . . . . . . . . . . . . . . . 52

1
6 Convergência 55

7 Soluções de Problemas Parabólicos 58

7.1 A Equação do Calor . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.2 Formulação Forte . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.3 Espaços de Funções . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.4 Formulação Fraca . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.5 Formulação de Galerkin . . . . . . . . . . . . . . . . . . . . . . . 59
7.6 Forma Matricial . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

8 Soluções de Problemas Hiperbólicos e Elı́pticos 64

2
Nota

Esta apostila tem como objetivo fornecer ao aluno uma visão geral do método
dos elementos finitos. Foi elaborada no intuito de atender as ementas dos cur-
sos de ”Método dos Elementos Finitos” do Bacharelado em Matemática Apli-
cada e Computacional do Instituto Multidisciplinar da UFRRJ (IM 475) e do
Mestrado Acadêmico em Modelagem Matemática e Computacional da UFRRJ -
PPGMMC. Trata-se de material didático sem fins lucrativos. Em algumas partes
ela reproduz partes de outras obras, devidamente referenciadas, sem que se
tenha intenção de reivindicar a autoria das mesmas, não oferecendo portanto
ofensa aos direitos autorais.

3
Chapter 1

Métodos Numéricos

1.1 Introdução
As equações diferenciais que descrevem matematicamente problemas práticos
de engenharia normalmente são difı́ceis de resolver. Com o advento dos com-
putadores de alta perfrmance se tornou possı́vel resolver tais equações. Para
tanto várias técnicas de solução numérica foram desenvolvidas, em especial o
método dos elementos finitos. Uma de suas principais vantagens é a possibil-
idade de se construir um programa computacional para analisar diversos tipos
de problemas, incluindo aqueles de geometria complexa com condições pre-
scritas. Vamos primeiramente introduzir os métos numéricos que diretamente
influenciaram a construção do método de elementos finitos.

Os chamados Métodos Numéricos são métodos aproximados que nos per-

mitem obter a solução de problemas complexos. Eles surgem da necessidade
de se lidar com a crescente complexidade dos modelos de engenharia. Os
métodos aproximados basicamente são de 2 tipos: métodos variacionais e de
resı́duos ponderados. Estudaremos brevemente dois métodos que exemplifi-
cam bem cada um dos tipos:

• Método de Rayleight-Ritz
No método a função y(x) é substituida por uma função aproximada v(x),
formada pela combinação linear de funções φi (x). Minimizando-se o fun-
cional, a solução aproximada é obtida. Trata-se portanto de um método
variacional.

• Método de Galerkin
O método não requer a existência de um funcional. Utiliza diretamente

4
 a equação diferencial que descreve o problema. Para obter a solução,
substitui-se na equação funções aproximadoras que devem satisfazer as
condições de contorno do problema.

Geralmente não é fácil obter as funções aproximadoras que satisfazem as

condições de contorno irregular e ainda garantir que elas se aproximem da
função exata quando utilizamos os métodos de Rayleight-Ritz e de Galerkin.
Para melhorar a precisão dos resultados é preciso ainda considerar funções
de ordem superior o que muitas vezes torna o cálculo inviável. Veremos que o
Método dos Elementos Finitos surge então a partir desses métodos mas com
a chamada discretização do domı́nio, onde se tem a subdivisão do domı́nio
em vários subdomı́nios conhecidos como elementos finitos; contornando parte
das dificuldades apresentadas pelos métodos anteriores. Antes de estudarmos
o Método dos Elementos Finitos vamos revisar os princı́pios do cálculo varia-
cional e os métodos de Rayleight-Ritz e de Galerkin.

1.2 O Cálculo Variacional

Os chamados funcionais exercem um importante papel em diversos problemas
práticos. Por funcional entendemos a correspondência entre um número real
definido com cada função (ou curva) pertencente a determinada classe.
Os primeiros resultados da área são atribuidos a Euler (1707-1703). O
Cálculo Variacional é o ramo mais desenvolvido do chamado Cálculo Funcional,
e trata da determinação de mı́nimos e máximos de funcionais, ou extremos de
funcionais.
Podemos definir os principais ingredientes do cálculo variacional como:

• Funcional - É uma entidade que depende de uma função (função de uma

função)

• Objetivos - Estudar métodos que permitam achar os valores extremos dos

funcionais.

• Exemplos - Princı́pio da mı́nima energia potencial, lei da conservação da

energia, etc.

O funcional tem a propriedade de localização que consiste no fato de que

se dividirmos a curva y = y(x) em partes e calcularmos o valor do funcional
para cada parte, a soma dos valores do funcional para cada parte separada é

5
igual ao valor do funcional para a curva inteira.

Definição: Seja F uma função dependente da função y e de sua derivada

y , representada pela notação: F = F (x, y, y 0 ). Diz-se que uma função é ad-
0

missı́vel quando, dentro de uma certa tolerância e satisfazendo as condições

de contorno, representa a função exata y(x). Assim, temos:

ȳ(x) = y(x) + λη(x) (1.1)

onde η(x) é uma função arbitrária que satisfaz as condições de contorno:

η(x1 ) = 0; x = x1
η(x2 ) = 0; x = x2

E ainda, λ é um número pequeno.

Observe que para a mesma função η(x) podemos obter uma famı́lia de
funções ȳ(x), bastando variar λ. Para cada x fixo, λη(x) é uma variação da
função exata, ou seja:

δy = λη(x) (1.2)

Vamos derivar (1.1) em relação a x:

ȳ 0 (x) = y 0 (x) + λη 0 (x)

Vamos agora derivar (1.2) em relação a x:

δy 0 = λη 0 (x)

A variação da função F (x, y, y 0 ) associada à variação de y (e consequente-

mente de y 0 ) para um determinado valor de x será:

∆F = F (x, ȳ(x), ȳ 0 (x)) − F (x, y, y 0 )

ou
∆F = F (x, y + λη, y 0 + λη 0 ) − F (x, y, y 0 )

Se desenvolvermos o primeiro termo do lado direito em uma série de Taylor,

teremos:

6
 ∂F 0 1 ∂ 2 F
 
∂F 2
∆F = F + λη + 0 λη + (λη) + . . . +
∂y ∂y 2! ∂y 2

∂F ∂F
−F = λη + 0 λη 0 + . . .
∂y ∂y

Substituindo-se λη por δy e λη 0 por δy 0 , temos:

∂F ∂F
∆F = δy + 0 δy 0 + . . . (1.3)
∂y ∂y

Observe a derivada de uma função f (u, v) de duas variáveis u e v: df =

∂f
∂u
+ ∂f
du ∂v
dv.

Comparando-se esta derivada com a equação (1.3), observa-se que os dois

primeiros termos de (1.3) são análogos aos da diferencial total. Então, por
analogia, vamos tratar F como uma diferencial total:

∂F ∂F
δF = δy + 0 δy 0
∂y ∂y

onde δF é o operador variacional, que possue as mesmas propriedades do

operador diferencial, ou seja:,

 
d dF
(δF ) = δ
dx dx

Z x2 Z x2
δ F dx = δF dx
x1 x1

Sendo assim, seja F (x, y, y 0 ) uma função com primeira e segunda derivadas
contı́nuas para x1 ≤ x ≤ x2 e que satisfaz as condições de contorno y(x1 ) =
A; y(x2 ) = B. Vamos encontrar a função para a qual o funcional,

7
 Z x2
Y = F (x, y, y 0 ) dx
x1

tem um extremo fraco. Ou seja, onde as funções admissı́veis são todas as

curvas suaves que ligam os pontos A e B. A variação de F é dada por:

Z x2  
∂F ∂F
δY = δy + 0 δy 0 dx
x1 ∂y ∂y

Semelhante ao que ocorre com as funções, a condição de extremo é deter-

minada pela primeira variação do funcional igual a zero.

Para encontrar a função que minimiza ou maximiza o funcional e o torna

estacionário, devemos satisfazer certas condições de contorno. As funções
que satisfazem as condições de contorno são chamadas de admissı́veis, como
visto na Definição anterior:

Z x2
Ȳ = F (x, ȳ, ȳ 0 ) dx
x1

Para λ = 0, a função Ȳ concide com a função exata (lembre-se que ȳ =

y + λη).

dȲ
δ Ȳ = Ȳ 0 δλ = δλ
dλ
é então a primeira variação de Ȳ em relação a λ.

A Condição de Estacionariedade implica em:

 
dȲ
=0
dλ λ=0

Então:

8
 x2 x2
∂F dȳ ∂F dȳ 0
  Z  Z  
dȲ dF
= dx = + dx
dλ λ=0 x1 dλ λ=0 x1 ∂ ȳ dλ ∂ ȳ 0 dλ λ=0

Como ȳ(x) = y(x) + λη(x) e ȳ 0 = y 0 + λη 0 :

dȳ
= η(x)
dλ

dȳ 0
= η 0 (x)
dλ

  Z x2   Z x2  
dȲ ∂F ∂F ∂F ∂F
= η + 0 η0 dx = η + 0 η0 dx
dλ λ=0 x1 ∂ ȳ ∂ ȳ λ=0 x1 ∂y ∂y

Integrando por partes:

  Z x2      x 2
dY ∂F d ∂F ∂F
= η− η dx + η
dλ λ=0 x1 ∂y dx ∂y 0 ∂y 0 x1

Como η(x1 ) = 0 e η(x2 ) = 0.

  Z x2   
dY ∂F d ∂F
= − η dx = 0 (1.4)
dλ λ=0 x1 ∂y dx ∂y 0

que vale para qualquer função η.

Rx
Sendo Ȳ = x12 F (x, ȳ, ȳ 0 ) dx e lembrando que para λ = 0, Ȳ concide com o
valor extremo do funcional Y .

Podemos expandir Ȳ na série de Taylor:

d2 Ȳ
   
dȲ 1
Ȳ = Y + λ + λ2 + ...
dλ λ=0 2! dλ2

9
 Podemos ainda usar a notação variacional,

1 2
Ȳ = Y + δY + δ Y + ...
2!
 
dY 2 d2 Y

2
onde δY = λ dλ λ=0
e δ Y = λ dλ2 .
λ=0

Nesta notação a primeira variação se reduz à:

Z x2    Z x2   
∂F d ∂F ∂F d ∂F
δY = λ − η dx = − δy dx = 0
x1 ∂y dx ∂y 0 x1 ∂y dx ∂y 0

Já que δy = λη.

Observe que para que o funcional tenha um extremo deve ter mesmo sinal
para qualquer δy. Assim,

 
∂F d ∂F
− =0 (1.5)
∂y dx ∂y 0

que define a Equação de Euler Lagrange.

Definição: [Gelfand & Fomin] Seja Y um funcional da forma,

Z x2
Y = F (x, y, y 0 ) dx
x1

definido em un conjunto de funções y(x) as quais tem primeiras derivadas

contı́nuas em [x1 , x2 ] e satisfaz as condições de contorno y(x1 ) = A; y(x2 ) = B.
Então, uma condição necessária para que Y tenha um extremo para uma dada
função y(x) é que y(x) satisfaça a equação de Euler,

 
∂F d ∂F
− =0
∂y dx ∂y 0

10
1.2.1 Método Aproximado de Rayleight-Ritz
Lembrando o cálculo variacional onde se procura uma função y(x) que dentre
todas as funções admissı́veis é a solução exata para minimizar um determi-
nado funcional, vamos agora supor que se queira encontrar a função v(x) que
minimiza o funcional:

Z x2
Y = f (x, y, y 0 ) dx (1.6)
x1

e com as condições de contorno:

y(x1 ) = y(x2 ) = 0 (1.7)

Neste método, v(x) é conhecida como função aproximada ou aproximadora

e é formada pela combinação linear das funções φi (x), de modo que:

n
y(x) ∼
X
= v(x) = ai φi (x) (1.8)
i=1

As funções φi (x) são chamadas de funções de forma. Fica claro que sua
escolha é importante para uma boa aproximação da solução.
As funções de forma são linearmente independentes. O que as caracteriza
é que cada uma delas individualmente satisfaz as condições de contorno do
problema, ou seja:

φi (x1 ) = φi (x2 ) = ... = φi (xn ) = 0; i = 1, ..., n (1.9)

As funções φi são contı́nuas até o grau m − 1, onde m é a ordem da maior

derivada do funcional. Os coeficientes ai são os parâmetros nodais.

Substituindo-se (1.8) na equação (1.6), ou seja substituindo-se y por v e

impondo-se a condição de estacionaridade, apresentada na seção anterior,
temos:

∂Y ∂Y
δY = δa1 + ... + δan = 0 (1.10)
∂a1 ∂an
Que pode ser vista como o sistema homogêneo:

11
 ∂Y
=0
∂a1
.. (1.11)
.

∂Y
=0
∂an
Já que as variações δai são arbitrárias.

Espera-se que quanto maior n, melhor a solução. Para se obter soluções

convergentes as seguintes condições devem ser satisfeitas:

(i) v(x) deve ser contı́nua até uma ordem menor do que a maior derivada do
integrando;

(ii) φi (x) devem satisfazer individualmente as condições de contorno essen-

ciais;

(iii) A sequência de funções deve ser completa, o que ocorre quando:

Z x2 n
X
lim (y − ai φi )2 dx < λ
n→∞ xi i−1

para λ pequeno (condição de completidade).

A convergência do método pode ser verificada considerando-se duas ou

mais funções aproximadoras e comparando-se sucessivos valores Y (i) do fun-
cional minimizado obtido com a sequencia:

(1) (1)
v (1) = a1 φ1

(2) (2) (2) (2)

v (2) = a1 φ1 + a2 φ2
(1.12)
..
.

(n) (n)
v (n) = a1 φ1 + ... + a(n)
n φn
(n)

12
 Então, a seguinte condição é verificada:

Y (1) ≥ Y 2) ≥ · · · ≥ Y (n)

que é chamada convergência monotônica.

Exemplo: Viga biapoiada, página, 37, Assan [1].

1.2.2 O Método de Galerkin

Para resolver um sistema de equações diferenciais, o método utiliza direta-
mente a equação ou a chamada forma forte do problema, onde substitui-se as
funções aproximadoras que devem satisfazer as condições de contorno. Com
o uso dessas funções aproximadoras introduz-se resı́duos que são então pon-
derados (funções de ponderação).

Considere o sistema de equações lineares:

Lv = f (1.13)

Onde L é o operador linear e v satisfaz a certas condições de contorno.

Vamos admitir que v̄ é uma função aproximadora de v, da forma:

n
X
v̄ = ai φi ; i = 1, . . . , n (1.14)
i=1

Para se encontrar as funções ai de modo que v̄ melhor se aproxime de v,

define-se φi são como funções ponderadoras ou de ponderação. Assim, o peso
médio do resı́duo sobre o domı́nio será nulo.

Substituindo-se (1.14) em (1.13), surge o erro , dado por:

 = Lv̄ − f (1.15)

Segundo o método de Galerkin, é necessário obedecer uma condição de

ortogonalidade entre a função erro  e as funções ponderadoras φi , de modo
que:

13
 Z
(Lv̄ − f )φi dV = 0; i = 1, . . . , n (1.16)
V

A equação (1.16) é um sistema de n equações que permite calcular os coe-

ficientes ai que são determinados fazendo-se com que os resı́duos se anulem
a cada integração ponderada.

Exemplo: Problema da viga biapoiada, página 45, Assan [1].

14
Chapter 2

Método dos Elementos Finitos

2.1 Introdução
• O método dos Elementos Finitos (MEF) é uma ferramenta para cálculo
aproximado da solução de equações diferenciais que governam diversos
problemas fı́sicos. Tem grande utilidade prática e permite a solução de
problemas que sem um método numérico não seriam solucionados.

• O método tem origem nos estudos de análise estrutural o que explica

muito sua terminologia.

• O MEF utiliza ferramentas de outros métodos de aproximação (Rayleight-

Ritz e Galerkin), que são precursores não só do MEF mas de outros
métodos bastante conhecidos e utilizados como o método das diferenças
finitas, volumes finitos, elementos de contorno, etc.

2.2 O Problema de Valor no Contorno

• Quando se busca a descrição quantitativa de um fenômeno fı́sico os mod-
eladores estabelecem um sistema de equações diferenciais ordinárias
(EDO’s) ou parciais (EDP’s) válidas em uma certa região (ou domı́nio) e
impoem neste sistema adequadas condições iniciais e de contorno, com-
pletando assim o problema (Formulação Clássica).

• No entanto, apenas equações mais simples com contorno geometrica-

mente simples, podem ser solucionadas exatamente. Um exemplo são

15
 as EDO’s com coeficientes constantes. Mesmo assim, quando um grande
número de variáveis dependentes estão envolvidas, encontramos dificul-
dades consideráveis.

Com o surgimeto dos computadores o MEF ganhou novo impulso. Ver o

trabalho de Las Casas [9] para maiores referências.

2.3 Formulação Clássica - Forma Forte (S)

Podemos dizer que os ingredientes principais do MEF para a solução de prob-
lemas de valor no cortorno são:

• A formulação variacional ou forma fraca do problema.

• A solução aproximada das equações variacionais através do uso de funções

especiais de elementos finitos.

Vamos supor que queremos resolver o seguinte problema:

d2 u
+f =0 (2.1)
dx2
Onde f é uma função suave definida no intervalo[0, 1]. A equação é sujeita
as seguintes condições de contorno:
u(1) = g
(2.2)
du
− (0) = h
dx
Onde g e h são constantes conhecidas. A primeira condição é conhecida
como essencial ou de Dirichlet, enquanto que a segunda é chamada de natural
ou de Neumann.

• Essa formulação é conhecida como Formulação Clássica e define a chamada

Forma Forte (S) (Strong) do problema.

• A Forma Forte nem sempre possui solução analı́tica, daı́ o uso de métodos
aproximados.

16
 Dado: f : Ω̄ → R, as constantes g e h

Encontre u : Ω̄ → R, tal que:

d2 u
+f =0
dx2

u(1) = g (2.3)

du
− (0) = h
dx

2.4 Formulação Variacional - Forma Fraca (W)

A solução u será obtida de forma aproximada. Uma função candidata à solução
é a chamada função teste ou trial function, ū. As funções teste precisam aten-
der as condições de contorno essenciais (Dirichlet):

ū(1) = g (2.4)

Limitamos a escolha de ū ao conjunto das funções que possuem quadrado

integrável:

Z 1  2
dū
dx < ∞ (2.5)
0 dx

Logo u ∈ H 1 ou H 1 (Ω), espaço de Hilbert. Essa escolha garante a completi-

dade, ou seja, garante que o limite da sequência convergente de elementos no
espaço esteja também no espaço (espaço completo em sua norma natural).

δ = {ū ∈ H 1 , ū(1) = g}

Ou,

17
 Z 1  2
du
δ = {ū|ū(1) = g, dx < ∞} (2.6)
0 dx

Assim, a solução aproximada ū será considerada em uma região Ω de con-

torno Γ. Logo,

n
u∼
X
= ū = ai Ni (2.7)
i=1

Onde as funções Ni são as chamadas funções de forma, escolhidas de

modo a garantir que a aproximação melhore conforme n aumente.
Como ū é uma aproximação de u, introduzimos o erro ou resı́duo na aproximação.
Para reduzir o erro podemos exigir que a integral do erro sobre o domı́nio Ω seja
ponderada de modo que o erro se anule. Podemos então escrever a formulação
de resı́duos ponderados:

1
d2 ū
Z  
wi +f dΩ = 0 (2.8)
0 dx2

Ou ainda

1 1
d2 ū
Z   Z
wi dx + wi f dx = 0 (2.9)
0 dx2 0

Nas equações (2.8) e (2.9) usamos as chamadas funções peso ou de ponderação

wi , que tem a propriedade de se anular no contorno (g = 0).

ν = w|w ∈ H 1 , w(1) = 0

(2.10)

Note que a equação (2.9) nos exige a continuidade da primeira derivada

das funções teste. Podemos relaxar essa exigência integrando a primeira parte
desta equação por partes:

18
 1 1 Z 1
d2 ū
Z   
dū dū dw̄i
wi dx = wi − dx (2.11)
0 dx2 dx 0 0 dx dx

Usando-se as condições de contorno em (2.11) e ainda substituindo-se em

(2.9), temos:

Z 1 Z 1
dū dw̄i
dx = wi f dx + wi (0)h (2.12)
0 dx dx 0

Agora, ū precisa ser apenas contı́nua.

Desta forma, podemos escrever a forma fraca (W) (Weak) do problema:

Dados: f : Ω̄ → R, constantes g, h

Encontre:ū ∈ δ, tal que ∀w ∈ ν

Z 1 Z 1
dū dw̄
dx = wf dx + w(0)h (2.13)
0 dx dx 0

Observe que δ e ν representam o mesmo espaço de funções, ou seja, H 1 .

A base deste espaço tem dimensão finita e pode ser dada, por exemplo, por
funções polinomiais de grau p, contı́nuas em [0, 1].

As formas forte e fraca são equivalentes ou seja, a solução de (S) é solução

de (W) e vice-versa.

Sejam as formas bilineares simétricas:

Z 1
dū dw̄
a(w, u) = dx (2.14)
0 dx dx

Z 1
(w, f ) = wf dx (2.15)
0

Logo, a equação (2.13) pode ser escrita como:

19
 a(w, u) = (w, f ) + w(0)h (2.16)

Observe que:

a(w, u) = a(u, w)
(w, f ) = (f, w)

Assim, podemos reescrever a formulação fraca (W):

Dados: f : Ω̄ → R, constantes g, h

Encontre:ū ∈ δ, tal que ∀w ∈ ν

a(w,u) = (w,f) + w(0) h

Note que δ e ν tem dimensões infinitas. O próximo passo é construir aproximações

de dimensões finitas. Vamos denotar estas novas coleções como δ h e ν h ,
subconjuntos de δ e ν, respectivamente, onde h faz referência a chamada
discretização do domı́nio, parametrizado pelo comprimento caracterı́stico h.
Assim, vamos considerar:

δ h ⊂ δ; se uh ∈ δ h , então uh ∈ δ
ν h ⊂ ν; se uh ∈ ν h , então uh ∈ ν

Logo, uh (1) = g e wh (1) = 0.

Observe que ν h e ν são espaços vetoriais lineares, logo:Sejam v1 e v2 ∈

ν → c1 v1 + c2 v2 ∈ ν, c1 , c2 ∈ R.

No entanto, considere u1 e u2 ∈ δ → u1 + u2 ∈
/ δ, já que u1 (1) + u2 (1) =
g + g = 2g ∈
/ δ.

Assim δ h é um espaço “aproximado”.

20
2.5 Formulação de Galerkin (G)
Observe que o método de Galerkin então é um método para obtenção de
soluções aproximadas baseado na formulação fraca do problema. Na verdade,
a escolha de espaços de aproximação de dimensão finita δ h e ν h não é atributo
do MEF mas do método de Galerkin quando este constroi uh do seguinte modo:

Considere que para cada função v h ∈ ν h seja construida uma função uh ∈

δ h , de modo que:

uh = v h + g h (2.17)

Onde g h é uma função dada que satisfaz a condição de contorno g h (1) = g.

Observe que g h , δ h e ν h são compostos de funções idênticas e ainda com a

particularidade de uh atender a condição de contorno:

uh (1) = v h (1) + g h (1) = 0 + g = g

já que v h ∈ δ h , conjunto das funções peso.

Vamos escrever a formulação fraca em termos de espaços de dimensão

finita:

a(wh , uh ) = (wh , f ) + wh (0)h (2.18)

Substituindo-se (2.17) em (2.18),

a(wh , v h ) = (wh , f ) + wh (0)h − a(wh , g h ) (2.19)

Observe que todos os valores conhecidos foram arranjados no lado direito

da equação (2.19): f, h e g.

Dados: f, g e h, como antes definidos

Encontre: uh = v h + g h , onde, v h ∈ ν h , ∀wh ∈ ν h

21
 a(wh , v h ) = (wh , f ) + wh (0)h − a(wh , g h )

Assim, (G) é uma nova versão de (W), mas em termos de dimensão finita.

Observação: Em [6] encontra-se um estudo interessante sobre o método

de Galerkin que mostra que o método leva a soluções exatas nos nós e tem
precisão de 2a ordem para as detivadas no ponto médio (página 24).

2.6 Formulação Matricial (M)

Note que a formulação de Galerkin pode ser traduzida como um sistema acoplado
de equações algébricas lineares de funções conhecidas. Seja ν h o espaço das
funções peso, como conjunto de todas as combinações lineares das funções
de forma NA : Ω̄ → R, onde A = 1, . . . , n. Ou seja,

n
X
wh = CA NA = C1 N1 + · · · + Cn Nn (2.20)
A=1

Onde, CA ; A = 1, . . . , n, são constantes.

Por coerência é preciso que as funções de forma satisfaçam NA (1) = 0.

Desta forma garantimos que wh (1) = 0. Observe que agora o espaço ν h tem
dimensão n.

Para a definição dos membros de δ h , vamos fazer a seguinte consideração:

Nn+1 (1) = 1, de modo que:

g h = gNn+1 (2.21)

Assim, g h (1) = g.

A equação (2.17) então nos leva a definição de um elemento tı́pico de δ h .

Seja uh ∈ δ h tal que,

22
 n
X
h h h
u =v +g = dA NA + gNn+1 (2.22)
A=1

Onde d0A s são constantes.

Observe que uh (1) = v h (1) + g h (1) = 0 + g = g, que foi justamente o objetivo

inicial na definição dos elementos do espaço.

Substituindo-se as equações (2.20) e (2.22) em (2.19), vem:

n n
! n
! " n
# n
!
X X X X X
a CA NA , dB NB = CA NA , f + CA NA (0) h−a CA NA , gNn+1
A=1 B=1 A=1 A=1 A=1
(2.23)

Podemos ainda fazer:

n
X
CA GA = 0 (2.24)
A=1

Onde,

n
X
GA = a(NA , NB )dB − (NA , f ) − NA (0)h + a(NA , Nn+1 )g (2.25)
B=1

Já que estamos usando as formas bilineares. Como os CA0 s são arbitrários,
cada GA , A = 1, . . . , n deve ser identicamente nulo, logo:

n
X
0= a(NA , NB )dB − (NA , f ) − NA (0)h + a(NA , Nn+1 )g (2.26)
B=1

Observe que em (2.26) conhecemos todos os valores exeto os d0B s.

23
 Nosso problema agora consiste em um sistema com n equações e n incógnitas,
que na forma matricial fica:

n
X
KAB dB = FA ; A = 1, . . . , n (2.27)
B=1

Onde,

KAB = a(NA , NB ) (2.28)

FA = (NA , f ) + NA (0)h − a(NA , Nn+1 )g (2.29)

Ou ainda,

Kd = F (2.30)

Onde,

 
K11 . . . K1n
K = [KAB ] =  ... .. ..  (2.31)

. . 
Kn1 . . . Knn

 
 F
 1  
F = {FA } = .
. (2.32)
 .
 F  
n

 
 d1 
 
d = {dB } = .
. (2.33)
 .
 d  
n

24
 Nestas matrizes, K é a chamada matriz de rigidez, F é o vetor de forças, e
d é o vetor deslocamentos. Observe que a forma matricial (M), abaixo, é equiv-
alente à formulação de Galerkin (G).

Dada a matriz dos coeficientes K e o vetor F .

Encontre d, de modo que:

Kd = F

Logo, d é a solução do nosso problema.

Observações:
• Se K possui inversa, então a solução do problema será: d = K −1 F .

• Observe que K é simétrica.

• Como (G) é solução aproximada de (W), a boa aproximação vai depender

da escolha das funções de forma NA0 s e do tamanho da discretização n.

• Observe que: (S) = (W ) ⇔ (G) = (M ).

Pn+1
• Podemos por conveniência escrever: uh = A−1 NA dA onde dn+1 = g.

Todas as observações acima trazem consequências computacionais impor-

tantes.

2.7 A discretização e os Espaços de Elementos

Finitos
O MEF prevê a divisão do domı́nio de integração (contı́nuo) em um número
finito de pequenas regiões, os elementos finitos propriamente ditos, o que torna
o meio contı́nuo em discreto. Considerando-se o domı́nio do nosso problema
Ω̄ = [0, 1], vamos chamar de malha de elementos finitos à divisão do domı́nio

25
em n subintervalos (que não se transpassam), de modo que para xA ’s pontos
nodais, temos o subintervalo tı́pico [xA , xA+1 ], onde xA < xA+1 ; A = 1, . . . .n.

Na Figura abaixo temos a visão de como fica a divisão do domı́nio.

Figure 2.1. Discretização em Elementos Finitos

Nosso elemento finito terá como domı́nio [xA , xA+1 ], com comprimento car-
acterı́stico hA = xA+1 − xA . Os pontos de interseção dos elementos da malha
são chamados nós.

26
 Observações:

• Os comprimentos caracterı́sticos h não precisam ser iguais. Para efeito

de simplificação, primeiramente consideraremos h = nl = constante;

• Podemos chamar de parâmetro da malha a h = max hA ; A = 1, . . . , n;

• Quanto menor for h, mais refinada é a malha. Ou ainda, quanto maior for
n, mais refinada é a malha.

Como consequência da discretização, ao invés de se procurar uma função

admissı́vel que satisfaz as condições de contorno para todo o domı́nio, procura-
se funções admissı́veis definidas no domı́nio de cada elemento.

Assim, para cada elemento i, a função aproximadora v será formada por

variáveis definidas a nı́vel de elemento de modo que:

n
X
v= aj Nj (2.34)
j=1

onde Nj são as funções de forma.

A definição dos espaços aproximados ν h e δ h são casos especiais dos

chamados espaços de elementos finitos lineares por parte (piecewise). As
funções de forma mais simples e mais utilizadas nos elementos finitos unidi-
mensionais são as funções lineares por partes, definidas como:

(a) Para os nós internos (2 ≤ A ≤ n):

x−xA−1
NA (x) = hA−1
; xA−1 ≤ x ≤ xA
xA+1 −x
NA (x) = hA
; xA ≤ x ≤ xA+1

(b) Para os nós do contorno:

27
 x2 −x
N1 (x) = h1
, x1 ≤ x ≤ x2

x−xn
Nn+1 (x) = hn
, xn ≤ x ≤ xn+1

As funções assim definidas tem a aparência de um chapéu (ou telhado)

como pode ser visto na Figura abaixo.

Figure 2.2. Função hat ou função roof.

Podemos ainda usar a função delta de Kronecker δAB :

1, se A = B
δAB =
0, se A 6= B

Assim, podemos descrever as funções de forma como,

NA (xB ) = δAB (2.35)

Um elemento tı́pico wh ∈ ν h terá a forma:

n
X
h
w = CA NA (2.36)
A=1

28
 Podemos observar na Figura 2.3 que wh é contı́nua com trechos lineares
no interior de cada elemento e descontı́nua nos nós (Observe que wh , x será
constante por partes). A função wh assim definida é conhecida como função
degrau generalizada.

Figure 2.3. Função Degrau Generalizada.

Observações:

• wh é um polinômio linear em x restrito a cada domı́nio de elemento.

• Com respeito às condições de contorno essenciais wh (1) = 0, observe

que wh = 0 se e somente se cada CA = 0; A = 1, . . . , n.

Um elemento tı́pico de δ h será obtido somando-se g h = gNn+1 a um ele-

mento tı́pico de ν h , de modo a assegurar que uh (1) = g.

Propriedades da Matriz K:

R1 dNi dNj
R1 dNj dNi
• K é simétrica. Observe que Kij = 0 dx dx
dx = 0 dx dx
dx = Kji .

29
• Observe que os elementos da diagonal principal de K são positivos e > 0:
R1
2
Kii = 0 dN
dx
i
dx > 0.

• Como as funções de forma são construı́das de modo que se anulem fora

da vizinhança dos nós, vários dos elementos da matriz K são nulos. Ob-
serve que os suportes das funções de interpolação possuem interseção
nula se o nó i não está conectado ao nó j, resultando em Kij = 0. Assim,
a matriz K é esparsa e ainda com estrutura de banda.

• Pode-se verificar que K é positiva-definida, ou seja:

– C T KC ≥ 0, ∀ C ∈ <n

– dT Kd = 0 ⇒ d = 0, ∀ d ∈ <n

Como consequência temos um sistema com solução única. Os autoval-

ores da matriz K são reais e positivos.
 
K11 K12 0
e
 K21 K22 K23 
K= 
 K32 K33 K34 
0
e K43 K44

30
Chapter 3

O Elemento Finito Unidimensional

Linear

3.1 Propriedades Globais

Do ponto de vista global as funções de forma são definidas em todos os pontos
do domı́nio do PVC. É globalmente que se define as propriedades matemáticas
do MEF. Já no chamado ponto de vista local ou de elemento [6] tudo que se
aplica à visão global será resumido ou restrito à nı́vel de elemento.

Domı́nio: [xA , xA+1 ]

Nós: {xA , xA+1 }
Graus de Liberdade: {dA , dA+1 }
Funções de Forma: {NA , NA+1 }
Função de Interpolação: uh (x) = NA (x)dA + NA+1 (x)dA+1 ; x ∈ [xA , xA+1 ],
onde dA = uh (xA )

3.2 Propriedades Locais

Domı́nio: [ξ1 , ξ2 ]
Nós: {ξ1 , ξ2 }
Graus de Liberdade: {d1 , d2 }
Funções de Forma: {N1 , N2 }
Função de Interpolação: uh (ξ) = N1 (ξ)d1 + N2 (ξ)d2
Observe que a numeração local começa por 1.

Temos então que lidar com uma mudança de referência, o que implica em

31
mudança do sistema de coordenadas da descrição global para a local. Esse
”mapeamento” é feito através da transformação afim: ξ : [xA , xA+1 ] → [ξ1 , ξ2 ] tal
que:

ξ(xA ) = ξ1 = −1
ξ(xA+1 ) = ξ2 = 1

Podemos escrever ξ(x) como: ξ(x) = c1 + c2 x. Assim:

ξ(xA ) = c1 + c2 xA = −1
ξ(xA+1 ) = c1 + c2 xA+1 = 1

Determinamos as constantes c1 e c2 resolvendo o sistema acima, e então

definimos:

2x − xA − xA+1
ξ(x) = (3.1)
hA
onde hA = xA+1 − xA . Observe que chamamos de ξ(x) ao mapeamento e
de x ao ponto. De forma análoga temos:

hA ξ + xA + xA+1
x(ξ) = (3.2)
2
onde chamamos de x(ξ) ao mapeamento, e de ξ ao ponto.

Facilita a notação considerar a descrição global com subescritos em maiúsculas

e a descrição local com subescritos em minúsculas. Assim, por exemplo:
NA (x(ξ)) = Nae (ξ), onde xe : [ξ1 , ξ2 ] → [xe1 , xe2 ] = [xA , xA+1 ]

As funções de forma locais em termos de ξ, ficam:

1
Na (ξ) = (1 + ξa ξ), a = 1, 2 (3.3)
2
Assim,

32
 2
X
e
x (ξ) = Na (ξ)xea (3.4)
a=1

Voltando à formulação matricial (M), vamos levar em conta que as integrais

sobre o domı́nio [0, 1] em (2.31) a (2.33) podem agora ser escritas como a soma
de integrais no domı́nio de cada elemento e, Ωe = [xe1 , xe2 ]. Assim, a equação
(2.30) será agora composta das matrizes:

nel
X nel
X
K= K e, F = Fe (3.5)
e=1 e=1

onde: K e = [KAB
e
] e F e = {FAe }, dados por:

Z
e e dNA dNB
KAB = a(NA , NB ) = dx (3.6)
Ωe dx dx

Z Z
dNA dNn+1
FAe e
= (NA , f ) +δe1 δA1 h−a(NA , Nn+1 ) g = e
NA f dx+δe1 δA1 h− dx g
Ωe Ωe dx dx
(3.7)
Note que NA (x1 ) = δA1 e que Ne (x1 ) = δe1 .

e
Observe que KAB = 0 quando:

• A 6= e
• A 6= e + 1
• B 6= e
• B 6= e + 1
E ainda, FAe = 0 quando A 6= e ou A 6= e + 1, o que decorre diretamente da
definição das funções de forma. Temos então, como esperado, muitos termos
e
nulos em KAB . Podemos observar na matriz esquemática abaixo o que ocorre
com o e-ézimo elemento, onde x’s indicam os termos diferentes de zero:

33
 colunas
 : e e+1
 
 
 
 
Ke = 
 x x 
 linha e

 x x 
 linha e + 1
 
nXn
 

 


 


 


 

 
Fe = x linha e
x  linha e + 1

 

 


 


 
  
nX1

Fica claro que é mais econômico armazenar somente os termos diferentes

de zero. Para tanto, fazemos:

k e = [kab
e
]2x2 (3.8)

f e = {fae }2x1 (3.9)

Onde surgem a matriz de rigidez e o vetor de forças do elemento e:

Z
e e dNa dNb
kab = a(Na , Nb ) = dx (3.10)
Ωe dx dx

Z
fae = Na f dx + C̄ (3.11)
Ωe

Onde C̄ representa os demais termos no contorno e vale:

• δa1 se e = 1

34
 • 0 se e = 2, . . . , n − 1

e
• −ka2 g se e = nel

Note que K e e F e foram definidos com respeito à numeração global, en-

quanto que agora temos k e e f e definidos com respeito à numeração local, o
que exige que os termos não nulos sejam ”colocados” no lugar certo em relação
as matrizes globais. Isso é feito através do chamado assembly ou montagem.

3.2.1 Montagem da Matriz de Rigidez Global e do Vetor de

Forças Global
Queremos determinar (3.8) e (3.9) para todos os elementos da malha e ”en-
caixar” suas contribuições em K e F usando a técnica de assemby.

A informação para o assemby é armazenada na matriz de localização LM

[6]. Assim,


e, se a = 1
A = LM (e, a) =
e + 1, se a = 2

onde e é o número do elemento, e a são os graus de liberdade do elemento.

Operador Assembly:

K = Anel e
e=1 (k ) (3.12)

F = Anel e
e=1 (f ) (3.13)

Para o e-ézimo elemento:

e
Ke,e ← Ke,e + k11

e
Ke,e+1 ← Ke,e+1 + k12

35
 e
Ke+1,e ← Ke+1,e + k21

e
Ke+1,e+1 ← Ke+1,e+1 + k22

Fe ← Fe + f1e

Fe+1 ← Fe+1 + f2e

Algoritmo de Assembly adaptado de Cecka et al. [2].

1 . Leitura dos dados de entrada

2. Montagem da matriz LM

3. Inicializar K = 0 e F = 0

4. ∀e ∈ [1, . . . , nel] faça:

5. Cálculo de k e e f e (geralmente feito em uma sub-rotina)

6. Para o grau de liberdade a1 de e, faça:

7. F (LM (e, a1 )) = F (LM (e, a1 )) + f e (a1 )

8. Para o grau de liberdade a2 de e, faça:

9. K(LM (e, a1 ), LM (e, a2 )) = K(LM (e, a1 ), LM (e, a2 ))+k e (a1 , a2 )

10. Pare.

Ver Exemplo Hughes, página 42 [6].

3.2.2 Cálculo Explı́cito da Matriz de Rigidez:

Para o cálculo da matriz de rigidez de elemento vamos precisar de dois ingre-
dientes básicos: a troca de variáveis e a regra da cadeia:

36
Troca de Variáveis:
Seja f : [x1 , x2 ] → < uma função integrável e seja x : [ξ1 , ξ2 ] → [x1 , x2 ] continu-
amente diferenciável, com x(ξ1 ) = x1 e x(ξ2 ) = x2 . Então:

R x2 R ξ2
x1
f (x) dx = ξ1
f (x(ξ)) dx(ξ)
dξ
dξ

dx(ξ)
onde dξ
é o Jacobiano da transformação.

Regra da Cadeia:
Seja f uma função diferenciável,

df (x(ξ)) df (x(ξ)) dx(ξ)

dξ
= dx dξ

Cálculo de k e :
Primeiramente fazemos a troca de variáveis na integração:

dNa (x) dNb (x) R1 dNa

(x(ξ)) dN (x(ξ)) dx(ξ)
e
R
kab = Ωe dx dx
dx = −1 dx dx
b
dξ
dξ
P2
onde x(ξ) = a=1 Na (ξ) xa . Usaremos agora a regra da cadeia:

R1  −1
e dNa (ξ) dNb (ξ) dx(ξ)
kab = −1 dξ dξ dξ
dξ

Observe que:

dNa ξa (−1)a dNb ξb (−1)b

dξ
= 2
= 2
e dξ
= 2
= 2
 −1
dxe he dxe 2
dξ
= 2
e dξ
= he

Então:

e
R1 (−1)a (−1)b 2 R1 (−1)a+b
kab = −1 2 2 he
dξ = −1 he

37
 Logo,

(−1)2 1 (−1)3 −1
k11 = he
= he
, k21 = he
= he
,
(−1)3 4
k12 = he
= −1
he
, k22 = (−1)
he
= 1
he
,
 
1 1 −1
ke = he −1 1

3.2.3 Cálculo Explı́cito do Vetor de Forças:

Vamos lembrar que f = f (x) é uma função conhecida. Então:

P2
fh = a=1 fa Na , onde fa = f (x(ξa ))
R1
Na (ξ) f h (x(ξ)) dx(ξ)
R
fe = Ωe
f h (x) Na dx = −1 dξ
dξ + termos no contorno

he
R 1 P2 he
P2 R1
= 2 −1 b=1 fb Nb dξ = 2 b=1 −1
Na (ξ) Nb (ξ) dξ fb

Como Na (ξ) = 1/2(1 + ξa ξ), a = 1, 2 . Então Na Nb = 1/4(1 + ξa ξ)(1 + ξb ξ),

 
1 (1 − ξ)2 (1 − ξ)(1 + ξ)
Na Nb = 4 (1 + ξ)(1 − ξ) (1 + ξ)2

Então:

R1 1
R1 1 2
−1 4
(1 − ξ)2 dξ = −1 4
(1 + ξ)2 dξ = 3
R1 1 1
−1 4
(1 + ξ)(1 − ξ)dξ = 3
 
e he 2/3 1/3
f = 2
fb
1/3 2/3

Ver Exemplos de Revisão.

38
Chapter 4

O PVC Bi-dimensional

Seja o domı́nio aberto Ω ⊂ <nsd com contorno τ , onde nsd corresponde ao

número de dimensões espaciais, ou seja, para o caso bi-dimensional nsd = 2.
Vamos considerar x um ponto qualquer de <nsd . Para o caso bi-dimensional:

   
x1 x
x = {xi } = =
x2 y

E ainda, seja n o vetor unitário normal a τ :

   
n1 nx
n = {ni } = =
n2 ny

onde i = 1, . . . , nsd.
S T
Vamos considerar que τ admita a decomposição: τ = τg τh onde τh τg =
. Ou seja, juntos τh e τg compreendem
S todo o contorno e não há superposição
entre eles. Sendo assim, Ω̄ = Ω τ .
Para estabelecer as formulações de elementos finitos para o caso bi-dimensional,
vamos precisar de dois ingredientes básicos:

Teorema da Divergência:1

Seja f : Ω̄ → <, f é de classe C 1 (1a derivada contı́nua), então:

R R
Ω
f, i dΩ = τ
f ni dτ
1 ∂u ∂2u ∂2u ∂2u
Convenção: u, i = ∂xi e u, ii = ∂x21
+ ∂x22
+ ∂x23
.

39
 Integração por Partes:

Seja f como definida anteriormente e seja g : Ω̄ → <, também de classe

1
C . Então:

R R R
Ω
f, i g dΩ = − Ω
f g, i dΩ + τ
f g ni dτ

4.1 Exemplo: Equação do Calor:

Sejam: qi as componentes cartesianas do vetor de fluxo de calor, u as temper-
aturas e f o calor aplicado por unidade de volume. Considere o fluxo de calor
definido em termos do gradiente de temperaturas, com a equação constitutiva
dada pela lei de Fourier generalizada:

qi = −κij u, j

onde κij = κji são as condutividades térmicas, que para um meio ho-
mogêneo terão valor constante, enquanto que para um meio isotrópico serão
definidas como: κij = κ̄(x)δij , onde δij é a função delta de Kronecker. Note que
a matriz de condutividade κ = [κij ] é simétrica e positivo definida.

Formulação Forte:

Dados: f : Ω → <, g : τg → < e h : τh → <

Encontre: u : Ω̄ → < tal que,

qi, i = f em Ω

u = g em τg

−qi ni = h em τh

onde qi é definida pela lei de Fourier.

Observe que h são os valores prescritos de fluxo, enquanto que g são os

valores prescritos de temperatura.

40
 Espaço de Funções:

δ = {u ∈ H 1 |u = g em τg }

ν = {w ∈ H 1 |w = 0 em τg }

Formulação Fraca:

Dados: f : Ω → <, g : τg → < e h : τh → <

Encontre: u ∈ δ tal que ∀w ∈ ν

R R R
− Ω
w, i qi dΩ = Ω
w f dΩ + τh
w h dτ

onde qi é definida como: qi = −κ u, j.

Ou ainda:

a(w, u) = (w, f ) + (w, h)τ

onde as formas bilineares são:

R
a(w, u) = Ω w,R i κij u, j dΩ
(w, f ) = RΩ w f dΩ
(w, h)τ = τh w h dτ

Utilizando a notação de operador gradiente: ∇u = {u, i}, temos:

R
a(w, u) = Ω
(∇w)T κ ∇u dΩ

Formulação de Galerkin:

Sejam δ h ⊂ δ e ν h ⊂ ν espaços de dimensões finitas. Vamos assumir que

todos os membros de ν h se anulam no contorno τg e que cada membro de δ h
admite a representação: uh = v h + g h onde gh satisfaz a condição de contorno
u = g em τg .

Dados: f, g e h como na formulação fraca,

41
 Encontre: uh = v h + g h ∈ δ h tal que ∀wh ∈ ν h

a(wh , v h ) = (wh , f ) + (wh , h)τ − a(wh , g h )

Discretização do Domı́nio:

A divisão do domı́nio é feita de modo que: Ω̄ = nel e

S T e
e=1 Ω̄e e que Ωi Ωj = 
2
em elementos finitos triangulares, quadriláteros, etc . Em todos os casos,
podemos definir η = {1, 2, . . . , nnp} onde nnp é o número de pontos nodas
da malha.

Em relação às condições de contorno, considere os nós nos quais são pre-
scritas as condições de contorno de Dirichlet, uh = g. Vamos chama-los de
g-nós: ηg ⊂ η é o conjunto de g-nós. Seu complemento será: η − ηg , ou seja,
o conjunto dos nós para os quais queremos determinar uh . Assim, teremos
tantos elementos em η − ηg quanto equações a calcular. Logo, o número de
elementos do conjunto será neq, ou número de equações.

Em relação aos espaços de funções, assumimos que um membro tı́pico de

h
ν tem a forma:

X
wh (x) = NA (x) cA (4.1)
A∈η−ηg

onde NA é a função de forma associada ao nó A, e cA é uma constante.

Vamos assumir que wh = 0, se e somente se, cA = 0, ∀A ∈ η − ηg .

De forma semelhante,

X
v h (x) = NA (x) dA (4.2)
A∈η−ηg

onde dA é a incógnita nodal (temperatura). E ainda,:

2
Estudaremos mais adiante os elementos isoparamétricos.

42
 X
g h (x) = NA (x) gA (4.3)
A∈ηg

para gA = g(xA ). Observe que vamos admitir a aproximação dos valores no

contorno por funções de forma, neste caso, lineares por partes.

Substituindo (4.1) e (4.2) na formulação fraca:

a(wh , v h ) = (wh , f ) + (wh , h)τ − a(wh , g h )

X X
a(NA , NB )dB = (NA , f ) + (NA , h)τ − a(NA , NB )gB (4.4)
B∈η−ηg B∈ηg

Para A ∈ η − ηg .

Numeração Global das Equações:

O arranjo ID é conhecido como matriz de destinação, e relaciona o nó A ao
número da equação correspondente.


P, se A ∈ η − ηg
ID(A) =
0, se A ∈ ηg

onde 1 ≤ P ≤ neq, P é o número da equação. Logo a dimensão de ID é nnp.

Ver exemplo.

Forma Matricial:
Vamos estabelecer a forma matricial usando o arranjo ID na formulação de
Galerkin:

Kd=F

onde,

43
 K = [KP Q ]

d = {dQ }

F = {FP }

Para 1 ≤ P, Q ≤ neq. Lembrando que P = ID(A) e Q = ID(B).

Assim,

KP Q = a(NA , NB )
P
FP = (NA , f ) + (NA , h)τ − B∈ηg a(NA , NB ) gB

Observe que a matriz K continua sendo simétrica, positiva-definida, etc.

Cálculo da Matriz de Rigidez de Elemento e do Vetor de Forças de Ele-

mento:
Como no caso unidimensional, calculamos os arranjos globais como a soma
das contribuições dos elementos:

Pnel Pnel
K= e=1 K e, F = e=1 Fe

K e = [KPe Q ], F e = {FPe }

onde,

KPe Q = a(NA , NB )

FPe = (NA , f ) + (NA , h)τhe − B∈ηg a(NA , NB ) gB

P

No entanto, já sabemos que as matrizes podem ser calculadas com maior
eficiência pelo assembly de k e e f e , calculadas à nı́vel de elemento. Assim,
temos:

K = Anel e
e=1 k , F = Anel
e=1 f
e

k e = [kab
e
], f e = {fae }

44
 onde 1 ≤ a, b, ≤ nen, número de nós do elemento.

e
kab = a(Na , Nb )

fae = (Na , f ) + (Na , h)τhe − nen

P
b=1 a(Na , Nb ) gb

Arranjos de Processamento de Dados:

Para fins de programação podemos usar o vetor gradiente B = ∇Na , definido
como:BnsdXnen = [B1, B2, . . . , Bnen ]. Onde,

∂Na
 
Ba = ∂x
∂Na
∂y

e
R
Então: kab = a(Na , Nb ) = Ωe
BaT κ Bb dΩ.

Podemos também definir o arranjo que organiza os dados nodais dos ele-
mentos, IEN (Element Nodes Array), também chamado de arranjo de incidências
[6]. O arranjo relaciona a numeração local dos nós à sua numeração global.

IEN (a, e) = A

onde A é o número global do nó, a o número local do nó, e e o número do

elemento. Assim, o arranjo LM, definido no capı́tulo anterior será construı́do de
um modo mais geral como:

LM (a, e) = ID(IEN (a, e))

Assim,


0, se LM (e, a) 6= 0
gAe =
gA , se LM (e, a) = 0

Ver exemplo.

45
Chapter 5

Elementos Isoparamétricos:

A formulação isoparamétrica pode ser usada para produzir muitos tipos de ele-
mentos como: planos, sólidos, de casca, para mecânica da fratura, etc.

Resumidamente os elementos isoparamétricos devem ser tais que as funções

de forma garantem que quando h → 0, então uh → u. Para que isso ocorra os
elementos precisam ser conformes ou compatı́veis, ou seja devem satisfazer
as seguintes condições:

• Na tem que ser suave em Ω (pelo menos C 1 );

• Na tem que ser contı́nua em τ e ;

• Na tem que ser completa.

As condições de suavidade no interior do elemento e de continuidade nas

interfaces dos elementos garantem que a 1a derivada das funções de forma
tenham, no pior caso, saltos finitos nas interfaces dos elementos, assegurando
que as integrais necessárias na formulação sejam bem definidas. Note que
na formulação fraca de nosso problema temos derivadas de 1a ordem, porém
há problemas que requerem derivadas de ordem superior nos integrandos (ex.
viga de Bernoulli-Euler). Existem funções que podem melhor atender este caso,
como as funções de Hermite.

Elementos que satisfazem as condições de suavidade e continuidade são

ditos compatı́veis ou conformes. Observa-se no entanto que existem elemen-
tos não-conformes, que violam estas condições e são convergentes.

46
 Em relação a completidade vamos analisar as funções de forma usadas nos
exemplos dos capı́tulos anteriores: uh = nen e
P
a=1 a da . Neste caso as funções de
N
forma serão completas se:

dea = c0 + c1 xea + c2 yae

O que implica em:

uh (x)e = c0 + c1 x + c2 y
onde c0 , c1 , c2 são constantes arbitrárias. No caso unidimensional, as funções
de forma são: N1 = 1/2(1 − ξ) e N2 = 1/2(1 + ξ). Neste caso,

P2 P2
dea = c0 + c1 xea → uh = a=1 Na dea = a=1 Na (c0 + c1 xea ) =

(N1 + N1 )c0 + (N1 xe1 + N2 xe2 )c1 = c0 + xe c1 .

5.1 Elemento Unidimensional Linear:

Considere o ponto P pertencente a reta OX como indicado na Figura :

Figure 5.1. O Elemento Unidimensional Linear.

As coordenadas do ponto P (em coordenadas naturais) são: xP = ξ1 x1 +

ξ2 x2 , onde: ξ1 = LL1 e ξ2 = LL2 . Observe que ξ1 + ξ2 = 1

47
Podemos escrever.

    
1 1 1 ξ1
=
x x1 x2 ξ2

    
ξ1 1 x2 −1 1
=
ξ2 L −x1 1 x

Considere N1 = ξ1 e N2 = ξ2 .

∂N1 ∂N1 ∂ξ1 −1

∂x
= ∂ξ1 ∂x
= L

∂N2 ∂N2 ∂ξ2 1

∂x
= ∂ξ2 ∂x
= L

Logo,

∂N1 ∂N2 −1 1
 
B= ∂x ∂x
= L L

e a matriz de rigidez de elemento neste caso será:

−1
Z Z  
e T −1 1

k = B D B dΩ = L [D] dΩ
1 L L
Ωe Ωe L
 
1 0
Para D = ,
0 1

 
e 1 1 −1
k = (5.1)
L −1 1

Compare com o resultado encontrado no capı́tulo 3.

48
5.2 Elemento Triangular Linear:
Na construção do elemento triangular usaremos a ideia de coordenadas de
área. A Figura 5.2 tem área A = A1 + A2 + A3 , onde

A1 A2 A3
ξ1 = A
, ξ2 = A
, ξ3 = A

de modo que ξ1 + ξ2 + ξ3 = 1.

Figure 5.2. O Elemento Triangular.

O elemente deve atender a todas as condições de convergência dos ele-

mentos isoparamétricos. Assim,

xP = ξ1 x1 + ξ2 x2 + ξ3 x3

yP = ξ1 y1 + ξ2 y2 + ξ3 y3

Ou ainda,

    
 1  1 1 1  ξ1 
x =  x1 x2 x3  ξ2
y y1 y2 y3 ξ3
   

49
5.3 Elemento Quadrilátero Bilinear:
• Geometria:

Figure 5.3. O Elemento Quadrilátero.

• Interpolação:

4
X
x(ξ, η) = Na (ξ, η) xea (5.2)
a=1

4
X
y(ξ, η) = Na (ξ, η) yae (5.3)
a=1

Assim, as funções de forma serão obtidas a partir dos Polinômios de La-

grange (ver seção seguinte):

N1 = 41 (1 − ξ)(1 − η)

N1 = 14 (1 + ξ)(1 − η)

50
 Figure 5.4. Numeração Local dos Nós.

a ξa ηa
1 -1 -1
2 1 -1
3 1 1
4 -1 1

Table 5.1. Coordenadas Nodais

N1 = 14 (1 + ξ)(1 + η)

N1 = 14 (1 − ξ)(1 + η)

Ou ainda,

1
Na (ξ, η) = (1 + ξa ξ)(1 + ηa η) (5.4)
4
Assim, a solução aproximada será:

4
X
h
u (ξ) = Na (ξ) dea (5.5)
a=1

As condições de suavidade, continuidade e completidade..

51
5.4 Integração Numérica:
Com vimos nas seções anteriores, na determinação das matrizes de elemento
precisamos calcular integrais. Seja f : Ωe ∈ <nsd → < vamos calcular:

R
Ωe
f (x) dΩ

onde f (x) = B T B, por exemplo.

Vamos efetuar a integração no espaço de coordenadas naturais através de

uma mudança de variáveis.

• nsd = 1

R R1
Ωe
f (x) dΩ = −1
f (x(ξ)) x, ξ(ξ) dξ

• nsd = 2
R R1 R1
Ωe
f (x) dΩ = −1 −1
f (x(ξ, η), y(ξ, η)) J(ξ, η) dξ dη

• nsd = 3
R R1 R1 R1
Ωe
f (x) dΩ = −1 −1 −1
f (x(ξ, η, ζ), y(ξ, η, ζ), z(ξ, η, ζ)) J(ξ, η, ζ) dξ dη dζ

Ou seja, temos que avaliar numericamente uma integral da forma:

R1 R1
−1
... −1
g(ξ, . . . ) . . . dξ

Quadratura de Gauss
Para nsd = 1:

Z 1 nint
X nint
X
g(ξ) dξ = g(ξ¯l ) wl + R ' g(ξ¯l ) wl (5.6)
−1 l=1 l=1

52
 onde nint é o número de pontos de integração, ξ¯l são as coordenadas dos
pontos de integração, e wl é o peso do l-ézimo ponto de integração.

O método de Gauss é o método mais utilizado nas implementações do MEF.

Trata-se de um método de O(2 nint) . Para :

ξ¯1 = 0

• nint = 1
w1 = 2
 ¯
ξ1 = √ −1
; ξ¯2 = √13
• nint = 2 3
w1 = w2 = 1
( q q
ξ¯1 = − 35 ; ξ¯2 = 0; ξ¯3 = 35
• nint = 3
w1 = w3 = 95 ; w2 = 98

Para nsd = 2 :

nint(2) (1)
Z 1 Z 1 X nint X (1) (2) (1) (2)
g(ξ, η) dξ dη ' g(ξ¯l(1) , η̄l(2) ) wl(1) wl(2) (5.7)
−1 −1
l(2) =1 l(1) =1

Ou ainda,

Z 1 Z 1 nint
X
g(ξ, η) dξ dη ' g(ξ¯l , η̄l ) wl (5.8)
−1 −1 l=1

(1) (2) (1) (2)

onde: nint = nint(1) nint(2) , ξ¯l = ξl(1) , η̄l = ηl(2) , wl = wl(1) wl(2) .

Para o elemento quadrilátero bilinear (nint = 4), com o auxı́lio da Tabela 5.2
podemos calcular as coordenadas dos pontos de integração:

ξ¯1 = ξ¯3 = −1
√
3
, ξ¯2 = ξ¯4 = √1
3

−1 √1
η̄1 = η̄2 = √
3
, η̄3 = η̄4 = 3

w1 = w2 = w3 = w4 = 1

53
 l l(1) l(2)
1 1 1
2 2 1
3 1 2
4 2 2

Table 5.2. Localização dos Pontos de Integração.

Seminários:

Elemento triangular como degeneração do elemento quadrilátero e Elemen-

tos de ordem superior - polinômios de Lagrange.

Instabilidades na malha, instabilidade dos elementos (modos espúrios), tran-

camento.

54
Chapter 6

Convergência

Nossas soluções via MEF possuem erros numéricos devido a vários fatores [3],
por exemplo, o tipo, formato dos elementos finitos escolhidos em um grande
leque de possibilidades tem algumas escolhas melhores do que outras. Além
disso, sabemos que por melhor que seja nosso tratamento via MEF ou qualquer
outro método haverão erros presentes em nossas soluções já que nossos com-
putadores representam os números com um numero finito de digitos. Mesmo
se assumirmos que nosso programa é apropriado para a tarefa e função, que
ele é livre de erros de implementação e que a escolha dos elementos e de suas
ordens de quadratura são adequadas, nossas soluções exigem exame e bom
senso em sua avaliação. por isso informações como a convergência ganham
um papel importante na garantia de fazer o melhor possı́vel, com segurança.

Podemos considerar o erro na aproximação como: e(x) = u(x)−uh (x), onde

u(x) é a solução exata e uh (x) a solução via MEF. Supondo que u(x) seja con-
tinuamente diferenciável, ou seja, u(x) ∈ C ∞ , podemos representar a solução
através da série de Taylor:

du ∆2 x d2 u ∆p x dp u
u(∆x) = u|x=0 + ∆x |x=0 + |x=0 + · · · + |x=0 + . . . (6.1)
dx 2! dx2 p! dxp

Uma vez que nossa solução aproximada é construı́da a partir de polinômios

de grau p, espera-se um erro da O(∆x)p+1 . Então, para ∆x = h:

uh → e ≈ O(hp+1 )
duh
dx
→ e ≈ O(hp )

55
 d2 uh
dx2
→ e ≈ O(hp−1 )

dn uh
dxn
→ e ≈ O(hp+1−n )

Para que haja convergência todas as derivadas que aparecem na formulação

variacional devem ser representadas corretamente quando h → 0, ou seja
e → 0 quando h → 0 [10]. Para tanto,

p+1−n≥1

Ou ainda: p ≥ n

Condição de Completidade
A solução aproximada uh (x) deve ser representada por polinômios completos
de grau n. Em outras palavras, as funções usadas na aproximação devem ser
polinômios de grau n.

Condição de Compatibilidade
As derivadas da aproximação de ordem n − 1 devem ser contı́nuas nas inter-
faces dos elementos, ou seja C n−1 .

Vamos lembrar que para que se possa resolver o problema as integrais na

formulação fraca devem ser limitadas, ou seja:

R1 duh dwh
0 dx dx
dx < ∞

Observe que nas aproximações do MEF dos exemplos aqui tratados us-
amos funções de forma contı́nuas por partes. Essas funções são de quadrado
integrável. Observe também que se a ordem da mais alta derivada da formulação
é n, então as derivadas de ordem n − 1 devem ser contı́nuas nas interfaces dos
elementos.

56
 Figura6.jpg

Figure 6.1. Condição de Compatibilidade.

57
Chapter 7

Soluções de Problemas
Parabólicos

7.1 A Equação do Calor

Retornaremos a equação do calor estudada anteriormente em regime perma-
nente. Vamos agora trata-la como um processo dependente do tempo assim
como do espaço. A equação do calor é um caso tı́pico de equação parabólica.
Vejamos quais as diferenças e particularidades do MEF para esta classe de
equações.

7.2 Formulação Forte

Dados f, g, f e u0
Encontre u : Ω̄ × [0, T ] → < tal que:

ρ c u,t +qi,i = f em Ω× ]0, T [ (7.1)

u=g em τg × ]0, T [ (7.2)
−qi ni = h em τh × ]0, T [ (7.3)
u(x, 0) = u0 (x), x∈Ω (7.4)
A equação (7.1) é chamada equação do calor, onde u é a temperatura, ρ é
a densidade de massa do material, c a capacidade térmica especı́fica. Tanto ρ
quanto c são funções positivas de x ∈ Ω. O fluxo de energia térmica q é definido
como no capı́tulo anterior.

58
 A forte de calor volumétrica f é agora definida como : f : Ω × ]0, T [→ <,
onde f (x, t); x ∈ Ω e t ∈ ]0, T [, intervalo de tempo com T > 0.

As condições de contorno são definidas como: g : τg × ]0, T [→ < e h :

τh × ]0, T [→ <, onde τg e τh não variam com o tempo.

Observe que para que o problema seja bem posto é necessário definir uma
condição inicial de temperatura: u0 : Ω → <.

7.3 Espaços de Funções

O espaço das funções peso ν, não depende do tempo e será definido como no
capı́tulo anterior. Já o espaço das funções teste, δ, varia com t uma vez que a
condição de contorno (7.2) também varia:

δt = {u/u(x, t) = g(x, t); x ∈ τg ; u ∈ H 1 (Ω)}

7.4 Formulação Fraca

Dados f, g, h e u0
Encontre u(t) ∈ δt , t ∈ [0, T ] tal que ∀w ∈ ν,

(w, ρ c u̇) + a(w, u) = (w, f ) + (w, h)τ

(w, ρ c u(0)) = (w, ρ c u0 )

onde u̇ representa u,t .

7.5 Formulação de Galerkin

Dados f, g, h e u0 ,
Encontre uh = v h + g h onde uh (x, t) ∈ δth tal que ∀wh ∈ ν h ,

59
 (wh , ρ c v̇ h ) + a(wh , v h ) = (wh , f ) + (wh , h)τ − (wh , ρ c ġ h ) − a(wh , g h ) (7.5)

(wh , ρ c v h (0)) = (wh , ρ c u0 ) − (wh , ρ c g h (0))

Assim, aproximamos os espaços de funções de modo que uh (x, t) = v h (x, t)+

g (x, t) onde v h (t) ∈ ν h e g h (t) ∈ δth . Observe que v h (0) = v h (x, 0).
h

A equação (7.5) é chamada de equação semi-discreta já que o espaço foi

discretizado em x, mas o tempo t ainda é contı́nuo.

7.6 Forma Matricial

Do mesmo modo que fizemos nos capı́tulos anteriores, vamos representar v h e
g h em termos de funções de forma:

X
v h (x, t) = NA (x) dA (t) (7.6)
A∈η−ηg

X
g h (x, t) = NA (x) gA (t) (7.7)
A∈ηg

As funções de forma não dependem de t e são dadas como anteriormente.

Substituindo (7.6) e (7.7) na formulação semi-discreta (7.5), temos:

Dados F : ]0, T [→ <neq ,

Encontre d : [0, T ] → <neq tal que,

M d˙ + Kd = F ; t ∈]0, T [ (7.8)

d(0) = d0

60
 onde,

M = Anel
e=1 m
e

Para: me = [meab ]
R
meab = Ωe
Na ρ c Nb dΩ

K = Anel
e=1 k
e

Para: k e = [kab
e
]

e
R
kab = Ωe
BaT D Bb dΩ

F (t) = Fnos (t) + Anel e

e=1 f (t)

Para: f e = {fae }
R R Pnen
fae = Ωe
Na f dΩ + τh
Na h dτ − e
b=1 (kab gbe + meab ġbe )

Sendo que:

d0 = M −1 Anel be
e=1 d

Para: dbe = {dbea }

R Pnen
dbea = Ωe
Na ρ c u0 dΩ − b=1 meab gbe (0)

Nota-se que na forma matricial surge a matriz M conhecida como matriz de

massa (capacity). Podemos verificar que M é simétrica e positiva-definida.

A equação (7.8) define um sistema de equações diferenciais ordinárias acoplado

que vai exigir algoritmos especiais em sua solução. Os algoritmos mais usados
são os chamados métodos trapezoidais generalizados.

61
Algoritmo dos Métodos Trapezoidais Generalizados
M vn+1 + Kdn+1 = Fn+1

dn+1 = dn + ∆t vn+α

vn+α = (1 − α)vn + αvn+1

Onde: dn = d(tn ); vn = d(t˙ n ); Fn+1 = F (Tn + 1); ∆t é o passo de tempo

constante, e α é o parâmetro do qual depende a escolha do método:

(a) Método explı́cito:

α = 0 → Método de Euler ou Diferenças avançadas.

(b) Métodos Implı́citos:

α = 1/2 → Regra trapezoidal, Regra do Ponto Médio ou Método de Crank-

Nicolson.

α = 1 → Diferenças atrasadas.

São conhecidas na literatura duas diferentes formas de implementação:

Forma v:
- Dado d0 , determinamos v0 : M v0 = F0 − Kd0

- Para cada passo calculamos o preditor: d¯n+1 = dn + (1 − α)∆t vn

- Calculamos vn+1 em: (M + α∆t K)vn+1 = Fn+1 − K d¯n+1 (?)

- Calculamos dn+1 em: dn+1 = d¯n+1 + α∆t vn+1

Observe que quando α = 0 (método explı́cito), se assumirmos M como

uma matriz diagonal, não será necessário resolver a equação (?). Nos métodos
implı́citos é preciso resolver (?) a cada passo de tempo e o produto K d¯n+1 pode
ser feito como:

62
 K d¯n+1 = Anel e ¯e
e=1 k dn+1

onde den+1 contém zeros nas posições que correspondem aos valores pre-
scritos de temperatura.

Forma d:
..
.

- Determinar dn+1 em: 1

α∆t
(M + α∆t K)dn+1 = Fn+1 + 1
α∆t
M d¯n+1 (??)
dn+1 −d¯n+1
- Determinar vn+1 em: vn+1 = α∆t

Observe que se M é diagonal, a implementação facilita o cálculo do lado

direito de (??).

Para análise do método ver Hughes, página 462 [6], onde verifica-se que o
método é convergente, ou seja o algoritmo é estável e consistente, com pre-
cisão de 2a ordem, sendo assim adequado para o tratamento de problemas
de condução de calor. No entanto apresenta algumas desvantagens como a
necessidade de armazenagem de grande quantidade de dados, por exemplo.
Existem métodos que melhoram a solução e evitam certas desvantagens. Veja
por exemplo, o método EBE (Element-by-Element implicit method.)

63
Chapter 8

Soluções de Problemas
Hiperbólicos e Elı́pticos

-Problema hiperbólico - equação da onda

-Problema elı́ptico equação de Laplace

Diferenças principais na formulação e na solução quando comparadas as

equações parabólicas.

64
Bibliography

[1] Assan, A. E., Método dos Elementos Finitos. 1a ed., UNICAMP, 2010.

[2] Cecka, C., Lew, A. Darve, E., Introduction to Assembly of Finite Element
Methods on Graphics Processors. WCCM/APCOM 2010 (Consulta On line
21/11/12 - 16:00h ).

[3] Cook, R.D., Malkus, D.S., Plesha, M.E.Concepts and Applications of Finite
Element Analysis. John Wiley & Sons, 3rd. Edition, 1989.

[4] Cunha, M. C.,Métodos Numéricos. 2a ed., UNICAMP, 2000.

[5] Gelfand, I. M., Fomin, S. V. Calculus of variations, Dover Publications,

2000.

[6] Hughes, T. J. R.The Finite Element Method. 1st. ed., Prentice-Hall, 1987.

[7] Johnson, C.,Numerical Solution of Partial Differential Equations by Finite

Element Method, Dover Publications, 2009.

[8] Kwon, Y. W., Bang, H., The Finite Element Method Using Matlab. 2nd. ed.,
CRC Press, 2000.

[9] Las Casas, E. B.,Notas sobre o Método dos Elementos Finitos. UFMG
(Consulta online 08/11/12- 10:00h).

[10] Ribeiro, F. L. B., Notas de Aula de Introdução ao Método dos Elementos

Finitos. COPPE/UFRJ - PEC, 2004.

[11] Zienkiwicz, O. C, Morgan, K.,Finite Elements and Approximation. 1st. ed.,

Dover Publications, 2006.

65