ERMAC 2010: I ENCONTRO REGIONAL DE MATEMÁTICA APLICADA E COMPUTACIONAL
11 - 13 de Novembro de 2010, São João del-Rei, MG; pg 65 - 89
Introdução ao Método dos Elementos Finitos
J. A. J. Avila
Departamento de Matemática e Estatística - DEMAT, UFSJ
36307 – 904, São João Del - Rei, MG
[email protected]
RESUMO
Neste minicurso aprenderemos a importância do Método dos Elementos Finitos
em relação aos outros métodos numéricos e sua aplicabilidade em problemas
da ciência e da engenharia. Em particular, resolveremos numericamente, pelo
Método dos Elementos Finitos Galerkin, a Equação de Condução do Calor 2D
numa placa quadrada. Serão mostrados resultados numéricos da distribuição
de temperatura durante o regime transiente até o regime estacionário.
Palavras-chave: Método de Galerkin, Método dos Elementos Finitos, Equação
de Condução do Calor.
65
66
1. INTRODUÇÃO
O estudo de muitos fenômenos físicos que acontecem na engenharia, biologia, oceanografia,
astronomia, cosmologia, etc., usam domínios que são geometricamente muito complicados e
difíceis de desenhar, o qual dificulta sua resolução. O Método de Diferenças Finitas, num
primeiro momento, trataria de resolver-las com o uso de transformações conformes ou outras
transformações, porém, isto não é fácil, pois, envolve sistemas de equações diferenciais
parciais elípticas. É assim que uma nova técnica potentíssima chamada Método dos Elementos
Finitos nos ajuda a resolver estes tipos de problemas. O objetivo deste minicurso é conhecer e
saber aplicar o Método dos Elementos Finitos Galerkin na solução numérica de Equações
Diferenciais Parciais. Para compreender este método resolveremos numericamente, por
exemplo, a equação de condução do calor 2D numa placa quadrada unitária com condição
inicial e de contorno. Começaremos com a formulação clássica da equação de condução do
calor, logo passaremos à formulação integral, aproximaremos esta última formulação pelo
Método de Galerkin e discretizaremos o domínio pelo Método de Elementos Finitos.
Resultados numéricos serão apresentados.
67
2. EQUAÇÃO DE CONDUÇÃO DO CALOR EM 2D
Sabe-se que a transferência de calor é o transporte de energia num corpo material devido às
diferenças de temperatura, ou seja, sempre que existir uma diferença de temperatura em um
meio ou entre meios ocorrerá transferência de calor. Este fenômeno acontece em três
mecanismos diferentes:



Condução
Convecção
Radiação
Neste trabalho estudaremos a transferência do calor por condução, que é a troca de energia
entre as partes de um meio continuo que, estando em diferentes temperaturas, transferem
energia térmica pela transferência de energia cinética entre as partículas individuais ou grupo
de partículas, no nível atômico.
2.1 Propriedades Físicas
Uma propriedade física do material (ou meio onde ocorre a condução) se chama
condutividade térmica, k , que depende da natureza do material. Neste trabalho assumiremos
que o material é isotrópico.
A difusividade térmica do metal é definida por:

k
c
(1)
onde k é a condutividade térmica,  a densidade e c o calor específico. As unidades de
medidas para k ,  e c são
 k    W/  m  oC  ,     kg/m3  , c  J/  kg  oC 
(2)
Segundo a Eq.(1) e Eq. (2) a unidade de medida da difusividade térmica é,
2

 
W/  m  oC 
 k  
Wm2    J/s  m   m2 


     
 

J
 c     J/  kg  oC   kg/m3   J  
  s 


(3)
Em forma geral sabemos que k depende da temperatura. Neste trabalho o  é considerado
constante. Em HOLMAN (1986) podemos encontrar valores das propriedades de alguns metais,
veja Tabela 1.
68
Tabela 1. Propriedades de alguns metais a 20
Metal
k   W/  m  oC 
o
C
c   J/  kg  oC 
   kg/m3 
  m2 /s 
Prata
1,6563 104
Ouro
1, 27 104
Cobre
386,0
8,954
383,1
1,1234 104
Alumínio
8, 418 105
Aço
2,3 105
2.2 Formulação Diferencial
O domínio,   R2 , para a Equação de Condução do Calor é definido como,
   x, y   R 2 : 0  x  1, 0  y  1
(4)
onde  pode ser, por exemplo, uma placa quadrada de um certo metal. Neste trabalho 
será uma placa de cobre de 1 [m2 ] .
O contorno de  é definido por,
   
4
i
(5)
i1
Denote-se a distribuição transiente de temperatura pela função escalar u  u  x, y, t  .
A formulação matemática do problema de Condução do Calor 2D e sem fonte é
 u
 t   u ,  x, y, t      0, T 

u  x, y,0   0,  x, y   

o
o
u 1  u 2  u 4  100 C , u 3  500 C , t   0,T 

(6)
Segundo a forma dada da Eq. (6) é chamada de Problema de Valor Inicial e de Contorno – PVIC.
A temperatura nos cantos superiores da placa deve satisfazer:
u  0,1, t   u 1,1, t   100 oC,
t  0, T 
Na Figura 1 temos o domínio      com suas respectivas condições de contorno.
(7)
69
Figura 2.1. O domínio  e suas condições de contorno.
2.3 Formulação Integral
Na primeira equação de Eq. (6) multipliquemos por     x, y  e integremos em  ,

u
 dA     u  dA,   x, y 
 t

(8)
u
k
 dA    u  dA,   x, y 
 t
c 
(9)
u
 dA  k   u  dA,   x, y 
 t

(10)
Pela Eq. (1),

c
c

u
 u 
 dA  k     ds  k  u   dA,   x, y 


t
 n 
(11)
ou equivalentemente,
u
 u 
 dA  k  u   dA  k     ds,   x, y 
 t


 n 
c
(12)
onde n é o vetor normal à curva  . Eq. (12) é a Formulação Integral da Equação de Condução
do Calor 2D.
70
3. MÉTODO DOS ELEMENTOS FINITOS GALERKIN
Não é intenção de este minicurso abarcar todos os casos que pode ser estudado pelo Método
dos Elementos Finitos e sim elementos simples e com ordem de aproximação linear.
Estudaremos os elementos triangulares de Lagrange, ou seja, com nós nos vértices do
triângulo. No MEF existem dois tipos de estudos o global e o local.
3.1 Estudo Global
É quando se trabalha com todo o domínio do problema e com todos os nós que existem nele.
Definamos a solução aproximada global de u por
Nn
u  x , y , t    u i  t  i  x , y 
(1)
i 1
onde N n é o número de nós da discretização de  e i i n1 as funções bases globais. Como
N
Eq. Erro! Fonte de referência não encontrada. é valida para qualquer  então é valida para
uma família finita  j , j  1, N n , isto é,
u
 u 
 j dA  k  u   j dA  k   j   ds
 t


 n 
c
(2)
desse modo a Eq. (2) é, na verdade, um sistema de N n - equações.
A derivada temporal de u seria,
u  x, y, t 
t
N
  u i  t  i  x , y 
(3)
i 1
e suas derivadas parciais,
u  x, y, t 
x
N
  ui t 
i  x, y 
x
i 1
,
u  x, y, t 
y
N
  ui t 
i  x, y 
i 1
y
(4)
Substituindo Eq. (1) em Eq. (2),
 Nn

 Nn


 i 1

 i 1

 i 1
Nn
 c   u i  t i   j dA  k    u i  t i    j dA  k   j   u i  t 



i 
 ds,
n 
j  1, N n (5)
ou equivalentemente,
Nn
Nn
Nn
 ui t  c i j dA   ui t  k  i   j  dA   ui t  k   j
i 1

i 1

i 1

i
ds,
n
j  1, N n
71

Nn
i 1


i 1 
Nn
Nn
 c  i j dA u i  t    k   i   j  dA u i  t     k   j

i 1




i 
ds  u i  t ,
n 
j  1, N n
Assim, temos

 c   dA u t    k       dA  k 
Nn
j 1
Nn

i
j
j
j 1

i
j

i

 j ds  u j  t   0, i  1, N n
n

(6)
Sejam as matrizes globais,
M   mij 
K 1   kij1 
Nn  Nn
 c  i j dA
(7)
   j i  j 
 k   i  j  dA  k   i
+
 dA


 x x y y 
(8)
K 2   kij2 
Nn  Nn
Nn  Nn

 k

i
 j ds
n
(9)
onde M é a matriz massa e K  K 1  K 2 a matriz rigidez.
Logo,
mij u j   kij1  kij2  u j  0
(10)

 Mu  Ku  f

u  0  0


(11)
ou equivalentemente,
onde f  f  x, y   0 é o termo fonte. Podemos observar que Eq. (11) representa um
sistema linear de EDO.
3.2 Estudo Local
Chamado, também, estudo Elementar e é quando se trabalha com um elemento arbitrário do
domínio e com todos os nós que existem nele.
Na Figura 1 temos um elemento arbitrário e de  cujos nós são,
n1   x1 , y1  ,
n2   x2 , y2  , n3   x3 , y3 
(12)
72
Figura 1. Elemento arbitrário e com seus respectivos nós.
e
Defina-se a solução aproximada elementar de u por
3
u e  x, y, t    u ie  t  wi  x, y 
(13)
i 1
onde as bases locais lineares para cada elemento e são:
1
 x2 y3  x3 y2    y2  y3  x   x3  x2  y 
2 Ae 
1
w2  x, y   e  x3 y1  x1 y3    y3  y1  x   x1  x3  y 
2A
1
w3  x, y   e  x1 y2  x2 y1    y1  y2  x   x2  x1  y 
2A
w1  x, y  
(14)
onde Ae é a área de cada elemento e .
Para cada elemento arbitrário, e , definamos
a1  x3  x2
a2  x1  x3
a3  x2  x1
b1  y2  y3
b2  y3  y1
b3  y1  y2
d1  x2 y3  x3 y2
d 2  x3 y1  x1 y3
d3  x1 y2  x2 y1
(15)
O Jacobiano de um elemento arbitrário, e , é definido por
J e  a3b2  a2b3
(16)
e J  2A .
e
e
Substituindo equações (15) e (16) em Eq. (14) temos que as bases locais ficam
wi  x, y  
1
 di  bi x  ai y ,
Je
i  1,2,3
Com a única finalidade de simplificar os cálculos de integrais sobre e , transformaremos o
triângulo e num triangulo ê centrado na origem, tal como mostra a Figura 3.
(17)
73
Figura 2. Transformação do triangulo e no triangulo ê .
Considerando esta transformação, temos
w1  ,   1    
w2  ,   
(18)
w3  ,   
Fazendo alguns cálculos encontramos as coordenadas de área (ou coordenadas naturais),
 1 ,  2 ,  3  dadas por
i 
1
 di  bi x  ai y ,
Je
i  1,2,3
(19)
De aqui temos que as bases locais coincidem com as coordenadas de área,
w1  1 ,  2 ,  3    1
w2  1 ,  2 ,  3    2
(20)
w3  1 ,  2 ,  3    3
Desse modo, o cálculo da integral de uma função f  f  x, y  sobre um elemento qualquer,
e , seria
 f  x, y  dxdy  J  f  , 
e
e
eˆ
d d  J e  f  1, 2 , 3  d 1d 2
eˆ
(21)
As derivadas das coordenadas de área seriam
 1 b1

x J e
 2 b2
 e
x
J
 3 b3
 e
x
J
 1 a1
 e
y
J
 2 a2
 e
y
J
 3 a3
 e
y
J
(22)
74
f  x, y 
e para calcular a integral de
sobre e , temos
x
f
f  1 f  2 f  3
+

d 1d 2
 2 x  3 x
1 x
 x dxdy  J  
e
eˆ
e
(23)
f
f
f
  b1
+b2
 b3
d 1d 2
eˆ
 1
 2
 3
De forma análoga para a seguinte integral
f
f
 y dxdy   a 
+a2
eˆ 1
e
1
f
f
 a3
d 1d 2
 2
 3
(24)
A fórmula para calcular a integral sobre um elemento ê é dada por

eˆ
 2q 3r d 1d 2 
p
1
p!q!r !
,
 p  q  r  2 !
p, q, r 

0
(25)
Seja  ie um segmento do contorno de e . Definimos o comprimento de  ie por Li ,
L1  a32  b32 ,
L2  a12  b12 ,
L3  a22  b22
(26)
A fórmula para calcular a integral de linha sobre o segmento  ie é dada por

ie
w1p w2q ds  Li  eˆ  1p 2q ds  Li
i
p !q !
, p, q 
 p  q  1!

0
(27)
Procedendo de forma análoga, como no caso global, chegamos ao seguinte sistema linear de
EDO
e e
e e
e

M u  K u  f

ue  0  0


(28)
onde f e  0 é o termo fonte.
3.3 Cálculo das Matrizes e Vetores Elementares
Agora calcularemos as matrizes elementares: a matriz massa e a matriz de rigidez e os vetores
elementares: vetor carga.
3.3.1. Matriz “Massa”
M e   mije 
33
  c  wi w j dA
e
(29)
75
  c  wi w j dxdy   cJ e  wi w j d d   cJ e   i j d 1d 2
e
e
e
Então,
M e   mije 
33
  cJ e I ij
(30)
onde I ij    i j d 1d 2 é calculada pela formula dada em Eq. (25), isto é,
e
1
 12 se i  j
I ij    i j d 1d 2  
e
 1 se i  j
 24
(31)
3.3.2. Matriz “Rigidez”
A Matriz “Rigidez” elementar é dada por
K e  K 1e  K 2e

(32)
Cálculo de K 1e
K 1e   kij1e 
33
 k
e
w w
xi
xj
 wyi wyj  dA
(33)
então
kij1e  k 
e
w w
xi
xj
 wyi wyj  dA  k
  w w dA   w w dA  k  I
e
xi
xj
e
yi
yj
1
ij
 I ij2 
(34)
onde,
I ij1   wxi wxj dA
(35)
I ij2   wyi wyj dA
(36)
e
e
e
Calculando I ij1 temos que
I ij1   wxi wxj dA 
e

 
  i
1

  
b
 b2 i  b3 i  b1 j  b2 j  b3 j  dA
e eˆ  1
J
 2
 3   1
 2
 3 
  1
  i   j
 bk
 bs
  k   s

1
J e eˆ

1
 bk ik   bs js  dA
J e eˆ

 dA

onde  ij é o delta de Kronecker. Observe que para k  i e s  j , temos
(37)
76
I ij1 
1
Je
1
 b b dA  J
eˆ i
j
e
bib j Aeˆ 
1
bib j
2J e
(38)
Por tanto,
I ij1 
1
bi b j
2J e
(39)
De forma análoga temos para
I ij2 
1
Je
1
 a a dA  J
eˆ
i
j
e
ai a j Aeˆ 
1
ai a j
2J e
(40)
ou seja,
I ij2 
1
ai a j
2J e
(41)
Substituindo equações (39) e (41) em Eq. (34),
1
k
 1

 k  e bib j  e ai a j   e  bib j  ai a j 
33
2J
 2J
 2J
K 1e   kij1e 

(42)
Cálculo de K 2e
Com ajuda da Eq. (12), calculemos a normal exterior em cada lado do elemento e , veja Figura
4:
v1  n2  n1   x2  x1 , y2  y1    a3 , b3  
n1   b3 , a3 
v2  n3  n2   x3  x2 , y3  y2    a1 , b1  
n 2   b1 , a1 
v3  n1  n3   x1  x3 , y1  y3    a2 , b2  
n 3   b2 , a2 
(43)
Como as normais que precisamos devem ser unitárias, então
n11 
b3
,
n1
n12 
a3
,
n1
(44)
n 21 
b1
,
n2
n 22 
a1
,
n2
(45)
n 31 
b2
,
n3
n 32 
a2
,
n3
(46)
77
Figura 3. Fronteira de um elemento arbitrário.
Seja n o vetor normal à curva  
3
e
ie . Então
i1
K 2e   kij2e 
 k
NN

e
wi
w j ds  kII ij
n
(47)
onde
II ij  
e

1e
wi
w j ds   e  wi  n  w j ds

n
 wi  n1  w j ds    wi  n 2  w j ds    wi  n 3  w j ds
e
2
(48)
e
3
Para calcular a matriz II ij devemos saber qual é o lado do elemento que faz parte da fronteira
do domínio  .
Suponhamos que 1e   , então
II ij  
1e
w n
xi
11
+wyi n12  w j ds
Iija
a
Calculo de I ij
w
n11 +wyi n12  w j ds 
L1
 bi n11 +ai n12   j ds
eˆ
J e 1
L
L
1!
 1e  bi n11  ai n12   eˆ  j ds  1e  bi n11  ai n12 

1
J
J
1  1!
I ija  
1e
xi
L1
 bi n11  ai n12 
2J e
 b1n11  a1n12 b1n11  a1n12
L1 
 e b2 n11  a2 n12 b2 n11  a2 n12
2J

0
0

0
0 
0 
78
Suponhamos que e2   , então
II ij  
e2
w n
xi
21
+wyi n 22  w j ds
Iijb
b
Calculo de I ij
I ijb  
e2
w
xi
n 21 +wyi n 22  w j ds 
L2
 bi n 21  ai n 22 
2J e
0

b2 n 21  a2 n 22 
b3 n 21  a3 n 22 
0
0
L2 
 e 0 b2 n 21  a2 n 22
2J
0 b3 n 21  a3 n 22
Suponhamos que 3e   , então
II ij  
3e
w n
xi
31
+wyi n 32  w j ds
Iijc
c
Calculo de I ij
I ijc  
e2
w
xi
n 31 +wyi n 32  w j ds 
 b1n 31  a1n 32
L3 
 e
0
2J
b3 n 31  a3 n 32
L3
 bi n 31  ai n 32 
2J e
0 b1n 31  a1n 32 

0
0

0 b3 n 31  a3 n 32 
3.4 Domínio Discretizado
Na Figura 4, observamos uma malha para  . Esta malha nos fornece os seguintes dados:
Ne  12,
Nn  11,
Ni  3,
N c  N n  Ni  8
Figura 4. Malha para o domínio.
(49)
79
onde N e é o Número de elementos, N n é o Número de nós, N i é o Número de incógnitas e
N c é o Número de nós no contorno.
Por exemplo, para o elemento 1 , como mostra a Figura 5,
Figura 5. Elemento número 1 da Malha.
temos que para i  9 , j  3 e k  7 definimos a seguintes matriz elementar e o vetor
derivada elementar, respectivamente,
(1)
 m99

(1)
m(1)   mij(1)    m39
33
(1)
 m79

(1)
m93
(1)
m33
(1)
m73
(1)

m97
(1) 
m37  ,
(1) 
m77

u (1)
 u 9(1) 
 
  u 3(1) 
 u 7(1) 
 
(50)
Na Fig. 5, a escolha do nó i coincidir com o nó 9 é feita pelo software gerador de malha.
De forma análoga, temos para
k 1(1)
1(1)
 k99
 1(1)
  kij1(1)    k39
33
1(1)
 k79

1(1)
k93
1(1)
k33
1(1)
k73
1(1)

k97
1(1) 
k37  ,
1(1) 
k77

k 2(1)
 k992(1)

  kij2(1)    k392(1)
33
 k792(1)

k932(1)
k332(1)
k732(1)
k972(1) 

k372(1) 
k772(1) 
(51)
e
u (1)
 u 9(1) 
 
  u 3(1)  ,
 u (1)

 7 
f (1)
 f9(1) 


  f3(1) 
 f 7(1) 


(52)
Do mesmo modo faz para todos os elementos da malha.
3.5 Ensamblagem
Consiste em montar todas as N e matrizes elementares numa única matriz chamada de matriz
global cuja ordem é Nn  Nn .
Montagem
Uma maneira de visualizar todas as matrizes elementares numa única matriz é como
mostramos a continuação
3.5.1
80
 m(1)

 0


 0
 k

 0


 0

1(1)
0
0  k
 
0   0

 
 
k 1( Ne )   0
2(1)
0
k
1
0  u 


0   u  2 




m( Ne )   u  Ne  


0
m(2)
1(2)
0
0
k
(53)
1
1
0  u   f 
 


0    u  2   f  2 




 



k 2( Ne )    u  Ne    f  Ne  

 

2(2)
0
Obtenção da Matriz e Vetor Global
3.5.2
O gerador de malha nos fornece um arquivo de dados, assim como mostra a Tabela 1.
Tabela 1
j
Elemento i
1
9 3
2
6 3
3
5 2
4
9 2
5
10 1
6
8 1
7
5 9
8
7 4
9
11 4
10
11 9
11
11 8
12
9 11
k
7
9
9
6
5
10
10
11
8
7
10
10
Como o Ne  12 , então devemos elaborar 12 tabelas cada uma com sua respectiva matriz
elementar. A seguir estas tabelas,
1
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
2
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
81
3
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
4
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
5
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
6
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
7
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
8
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
9
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
10
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
11
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
12
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
82
Observe que a cor escura, em cada tabela, representa as componentes das matrizes
elementares. O passo seguinte é imaginar-se que estas tabelas representam matrizes de
ordem 1111 . Compreendido isso, agora só temos que somar estas 12 matrizes, resultando
em uma única matriz (ou tabela). Esta matriz é chamada matriz global.
Tabela 2. Matriz global de ordem 1111
1
2
3
4
5
1
511+611
0
0
0
551
2
0
322+422
0
0
352
3
0
0
133+233
0
0
4
0
0
0
844+944
0
6
0
426
236
0
0
7
0
0
137
847
0
8
618
0
0
948
0
9
0
329+429
139+239
0
359+759
10
51,10+61,10
0
0
0
55,10+75,10
11
0
0
0
84,11+94,11
0
0
874
5
515
325
0
0
355+555
+755
0
0
6
7
0
0
462
0
263
173
266+466
0
0
0
269+469
179+1079
0
0
0
87,11+107,11
0
984
0
0
0
177+877
+1077
0
8
681
0
0
68,10+118,10
98,11+118,11
392+492
193+293
0
395+795
296+496
197+1097
688+988
+1188
0
9
0
79,10+129,10
109,11+129,11
510,1+610,1
0
0
0
510,5
+710,5
0
0
610,8
+1110,8
199+299
+399+499
+799+1099
+1299
710,9+
1210,9
10
1110,11
+1210,11
0
0
0
811,4+
911,4
0
0
811,7+
1011,7
911,8+
1111,8
1011,9+
1211,9
510,10+
610,10+
710,10+
1110,10+
1210,10
1111,10
1211,10
11
811,11+911,11
+1011,11+
1111,11+1211,11
A estrutura da matriz global tem a seguinte forma
1
2
3
4
5
6
7
8
9
10
11
1
2
3
4
5
6
7
8
9
10
11
onde os quadradinhos brancos representam os zeros da matriz. De forma análoga, devemos
elaborar 12 colunas cada uma com seu respectivo vetor elementar. A seguir estas colunas,
1
1
2
3
4
5
6
7
8
9
10
11
2
3
4
5
6
7
8
9
10
11
12
83
Tabela 3. Vetor global de ordem 111
1
51+61
2
32+42
3
13+23
4
84+94
5
35+55+75
6
26+46
7
17+87+107
8
68+98+118
9
19+29+39+79+109+129
10
510+610+710+1110+1210
11
811+9110+1011+1111+1211
Agora expressemos em forma matricial
0
 m11

m22
 0
 0
0

0
 0
m
m52
 51
m62
 0
 0
0

0
 m81

m92
 0
 m10 1 0

0
 0
 k11

 0
 0

 0
k
 51
 0
 0

 k81

 0
 k10 1

 0
0
k22
0
0
k52
k62
0
0
k92
0
0
0
0
m33
0
0
m63
m73
0
m93
0
0
0
0
k33
0
0
k63
k73
0
k93
0
0
0
0
0
m44
0
0
m74
m84
0
0
m11 4
0
0
0
k44
0
0
k74
k84
0
0
k11 4
m15
m25
0
0
m55
0
0
0
m95
m10 5
0
k15
k25
0
0
k55
0
0
0
k95
k10 5
0
0
k26
k36
0
0
k66
0
0
k96
0
0
0
m26
m36
0
0
m66
0
0
m96
0
0
0
0
m37
m47
0
0
m77
0
m97
0
m11 7
0
0
k37
k47
0
0
k77
0
k97
0
k11 7
k18
0
0
k48
0
0
0
k88
0
k10 8
k11 8
m18
0
0
m48
0
0
0
m88
0
m10 8
m11 8
0
m29
m39
0
m59
m69
m79
0
m99
m10 9
m11 9
0
k29
k39
0
k59
k69
k79
0
k99
k10 9
k11 9
M g u g   K 1g  K 2 g  u g  f g
k1 10
0
0
0
k5 10
0
0
k8 10
k9 10
k10 10
k11 10
m1 10
0
0
0
m5 10
0
0
m8 10
m9 10
m10 10
m11 10
0   u1 

0   u 2 
0   u3 
 
m4 11   u 4 
0   u5 
 
0   u6  +
m7 10   u 7 
 
m8 11   u 8 

m9 11   u 9 
m10 11   u10 
 
m11 11   u11 
  u1   f1 
   
  u 2   f2 
  u 3   f3 
   
k4 11   u 4   f 4 
0   u 5   f5 
   
0   u 6    f6 
k7 10   u 7   f 7 
   
k8 11   u 8   f8 

k9 11   u 9   f9 
k10 11   u10   f10 
   
k11 11   u11   f11 
0
0
0
(54)
84
M gug  K gug  f g
(55)
 0 e K g  K 1g  K 2 g .
onde f
g
3.5.3
Obtenção da Matriz e Vetor a Resolver
Análise na matriz M . Usando as condições de contorno
 m11

 0
 0

 0
m
 51
 0
 0

 m81

 0
 m10 1

 0
0
m22
0
0
m52
m62
0
0
m92
0
0
0
0
m33
0
0
m63
m73
0
m93
0
0
0
0
0
m44
0
0
m74
m84
0
0
m11 4
m15
m25
0
0
m55
0
0
0
m95
m10 5
0
0
m26
m36
0
0
m66
0
0
m96
0
0
0
0
m37
m47
0
0
m77
0
m97
0
m11 7
m18
0
0
m48
0
0
0
m88
0
m10 8
m11 8
0
m29
m39
0
m59
m69
m79
0
m99
m10 9
m11 9
0   u1 

0   u 2 
0   u3 
 
m4 11   u 4 
0   u5 
 
0   u6 
m7 10   u 7 
 
m8 11   u 8 

m9 11   u 9 
m10 11   u10 
 
m11 11   u11 
m1 10
0
0
0
m5 10
0
0
m8 10
m9 10
m10 10
m11 10
De forma análoga para a matriz K . Desse modo temos
 m99
m
 10 9
 m11 9
m9 10
m10 10
m11 10
m9 11 
m10 11 
m11 11  N  N
i
i
 u 9   k99
u   k
 10   10 9
 u11   k11 9
m92
 0
m
 10 1 0
 0
0
m93
0
0
k9 10
k10 10
k11 10
0
0
m11 4
m95
m10 5
0
k9 11   u 9   f9 
k10 11   u10    f10  
k11 11   u11   f11 
m96
0
0
m97
0
m11 7
0 
m10 8 
m11 8 
M
 0 k92
k
 10 1 0
 0
0
k93
0
0
0
0
k11 4
Ni  N c
k95
k10 5
0
K
ou equivalentemente
k 96
0
0
k97
0
k11 7
0 
k10 8 
k11 8 
Ni  N c
 u1 
u 
 2
u3 
 
u 4  
u5 
 
u6 
u 
 7
 u 8 
 u1 
u 
 2
u3 
 
u 4 
u5 
 
u6 
u 
 7
 u 8 
85
M u  K u  f  M u  K u
(56)
M uK u f
(57)
onde f  0 . Assim,
onde f  f  M u  K u
3.6 Solução do Sistema de EDO
Malha temporal
O tamanho de passo temporal é definido por:
t 
T
N
(58)
e os nós temporais, por
tn  t0  nt ,
n  1,
N
(59)
onde t0  0 é o tempo inicial.
Condição Inicial
u  t0   0
(60)
u0  0
(61)
 M   t K  u n  M  1    t K  u n1  t  f n  1    f n 1 
(62)
Regra do Médio ponto Generalizado
De Eq. (60), temos que
Faça para n  1, N
Observe que Eq. (62) representa um Sistema de Equações Lineares, veja LEWIS et. al (2004).
Ax  b
(63)
onde
A  M   t K
x  un
(64)
b   M  1    t K  u n 1  t  f n  1    f n 1 
Para resolver o Sistema Linear em Eq. (63) usamos o método SOR. Lembrando que 0    1 .
Neste trabalho usamos   0,5 .
86
4. RESULTADOS NUMÉRICOS
O código computacional foi implementado na Linguagem Fortran (Compaq Visual Fortran). Os
programas foram rodados num laptop T6500 Processador Intel CoreTM 2 Duo. Para a geração
da malha se uso o software (livre) Gmsh 2.5.0 e para pós-processamento se uso o software
(comercial) Tecplot 360.
4.1 Malha Estruturada e não Estruturada
As Figuras 4.1 e 4.2 mostram dois tipos de malhas que foram geradas com o software Gmsh,
uma Malha Estruturada e uma Malha não Estruturada.
Figura 4.1 Malha Estruturada.
Figura 4.2 Malha não Estruturada
A Tabela 4.1 mostra o Número de elemento e o Número de nós das duas malhas.
Tabela 4.1 Elementos e nós das malhas
Ne
Nn
Tipo de Malha
Estruturada
1024
545
Não Estruturada
824
449
4.2 Resultados Numéricos
O Número de incógnitas que resultou da malha estruturada foi Ni  481 e da malha não
estruturada foi Ni  377 . O tempo CPU após atingir o estado permanente é mostrado na
Tabela 4.2.
Tabela 4.2 Número de incógnitas e tempo CPU.
Tipo de Malha
Estruturada
Não Estruturada
Ni
481
377
Tempo CPU em [s]
22
14
87
As Figuras 4.3 e 4.4 mostram a convergência ao longo do tempo da distribuição de
temperatura em ambas as malhas. Usou-se um t  5 103 e se atingiu o estado estacionário
em 2 [s] com N  400 . Também, pode-se observar que a distribuição de temperatura na
placa quadrada começa a diminuir uniformemente desde a parte superior até um pouco mais
da metade da placa em forma de elipses concêntricas. Por exemplo, o valor da temperatura no
centro da placa, isto é, em  0.5,0.5 é de 195,84[ oC ] .
Figura 4.3 Distribuição de temperatura na
malha estruturada.
Figura 4.4 Distribuição de temperatura na
malha não estruturada.
A seguir mostraremos alguns quadros, em instantes de tempo diferentes, do estado transiente
da condução do calor até atingir o estado estacionário.
t  0 [ s]
88
t  0,5 [s]
t  1 [ s]
t  2 [ s]
Figura 4.5 Distribuição de temperatura em alguns instantes de tempo, desde o tempo inicial
até o tempo em que atingiu o estado estacionário.
Desse modo concluímos a simulação numérica da condução do calor numa placa quadrada.
89
5. REFERÊNCIAS BIBLIOGRÁFICAS
HOLMAN J. P. Heat Transfer. 6a Edição. McGraw-Hill Book Company, New York, 1986.
LEWIS Roland W, NITHIARASU Perumal e SEETHARAMU, Kankanhally N. Fundamentals of the
Finite Element Method for Heat and Fluid Flow. John Wiley and Sons, 2004