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
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.
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 65
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.
66
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
3W/ m , kg/m , J/ kg o ok C c C (2)
Segundo a Eq.(1) e Eq. (2) a unidade de medida da difusividade térmica é,
22 2
3
W/ m J/s mWm m
J J sJ/ kg kg/m
o
o
Ck
c C
(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.
67
Tabela 1. Propriedades de alguns metais a 20 oC
Metal W/ m ok C
3kg/m J/ kg oc C
2 m /s
Prata 41,6563 10
Ouro 41,27 10
Cobre 386,0 8,954 383,1 41,1234 10
Alumínio 58,418 10
Aço 52,3 10
2.2 Formulação Diferencial
O domínio, 2R , para a Equação de Condução do Calor é definido como,
2, : 0 1, 0 1x y R x y (4)
onde pode ser, por exemplo, uma placa quadrada de um certo metal. Neste trabalho
será uma placa de cobre de 21 [m ] .
O contorno de é definido por,
4
1
i
i
(5)
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 é
1 2 4 3
, , , 0,
, ,0 0, ,
100 , 500 , 0,o o
uu x y t T
t
u x y x y
u u u C u C t 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:
0,1, 1,1, 100 , 0,ou t u t C t T (7)
Na Figura 1 temos o domínio com suas respectivas condições de contorno.
68
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 yt
(8)
Pela Eq. (1),
, ,u k
dA u dA x yt c
(9)
, ,u
c dA k u dA x yt
(10)
, ,u u
c dA k ds k u dA x yt
n
(11)
ou equivalentemente,
, ,u u
c dA k u dA k ds x yt
n
(12)
onde n é o vetor normal à curva . Eq. (12) é a Formulação Integral da Equação de Condução
do Calor 2D.
69
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
1
, , u ,nN
i i
i
u x y t t x y
(1)
onde nN é o número de nós da discretização de e 1
nN
i i
as funções bases globais. Como
Eq. Erro! Fonte de referência não encontrada. é valida para qualquer então é valida para
uma família finita j , 1, nj N , isto é,
j j j
u uc dA k u dA k ds
t
n
(2)
desse modo a Eq. (2) é, na verdade, um sistema de nN - equações.
A derivada temporal de u seria,
1
, ,u ,
N
i i
i
u x y tt x y
t
(3)
e suas derivadas parciais,
1 1
, , , , , ,u , u
N Ni i
i i
i i
u x y t x y u x y t x yt t
x x y y
(4)
Substituindo Eq. (1) em Eq. (2),
1 1 1
u u u , 1,n n nN N N
ii i j i i j j i n
i i i
c t dA k t dA k t ds j N
n
(5)
ou equivalentemente,
1 1 1
u u u , 1,n n nN N N
ii i j i i j i j n
i i i
t c dA t k dA t k ds j N
n
70
1 1 1
u u u , 1,n n nN N N
ii j i i j i j i n
i i i
c dA t k dA t k ds t j N
n
Assim, temos
1 1
u u 0, 1,n nN N
ii j j i j j j n
j j
c dA t k dA k ds t i N
n (6)
Sejam as matrizes globais,
n n
ij i jN NM m c dA
(7)
1 1 + n n
j ji i
ij i jN NK k k dA k dA
x x y y
(8)
2 2
n n
i
ij jN NK k k ds
n (9)
onde M é a matriz massa e 1 2K K K a matriz rigidez.
Logo,
1 2u u 0ij j ij ij jm k k (10)
ou equivalentemente,
u u
u 0 0
M K f
(11)
onde , 0f f x y é 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,
1 1 1 2 2 2 3 3 3, , , , ,n x y n x y n x y (12)
71
Figura 1. Elemento arbitrário e com seus respectivos nós.
Defina-se a solução aproximada elementar de eu por
3
1
, , u ,e e
i i
i
u x y t t w x y
(13)
onde as bases locais lineares para cada elemento e são:
1 2 3 3 2 2 3 3 2
2 3 1 1 3 3 1 1 3
3 1 2 2 1 1 2 2 1
1,
2
1,
2
1,
2
e
e
e
w x y x y x y y y x x x yA
w x y x y x y y y x x x yA
w x y x y x y y y x x x yA
(14)
onde eA é a área de cada elemento e .
Para cada elemento arbitrário, e , definamos
1 3 2 1 2 3 1 2 3 3 2
2 1 3 2 3 1 2 3 1 1 3
3 2 1 3 1 2 3 1 2 2 1
a x x b y y d x y x y
a x x b y y d x y x y
a x x b y y d x y x y
(15)
O Jacobiano de um elemento arbitrário, e , é definido por
3 2 2 3
eJ a b a b (16)
e 2e eJ A .
Substituindo equações (15) e (16) em Eq. (14) temos que as bases locais ficam
1
, , 1,2,3i i i iew x y d b x a y i
J (17)
Com a única finalidade de simplificar os cálculos de integrais sobre e , transformaremos o
triângulo e num triangulo e centrado na origem, tal como mostra a Figura 3.
72
Figura 2. Transformação do triangulo e no triangulo e .
Considerando esta transformação, temos
1
2
3
, 1
,
,
w
w
w
(18)
Fazendo alguns cálculos encontramos as coordenadas de área (ou coordenadas naturais),
1 2 3, , dadas por
1
, 1,2,3i i i ied b x a y i
J (19)
De aqui temos que as bases locais coincidem com as coordenadas de área,
1 1 2 3 1
2 1 2 3 2
3 1 2 3 3
, ,
, ,
, ,
w
w
w
(20)
Desse modo, o cálculo da integral de uma função ,f f x y sobre um elemento qualquer,
e , seria
1 2 3 1 2ˆ ˆ, , , , e e
e e ef x y dxdy J f d d J f d d (21)
As derivadas das coordenadas de área seriam
1 1 1 1
2 2 2 2
3 3 3 3
e e
e e
e e
b a
x J y J
b a
x J y J
b a
x J y J
(22)
73
e para calcular a integral de ,f x y
x
sobre e , temos
1 2 31 2ˆ
1 2 3
1 2 3 1 2ˆ1 2 3
+
+
e
e e
e
f f f fdxdy J d d
x x x x
f f fb b b d d
(23)
De forma análoga para a seguinte integral
1 2 3 1 2ˆ1 2 3
+ e e
f f f fdxdy a a a d d
y
(24)
A fórmula para calcular a integral sobre um elemento e é dada por
1 2 3 1 2 0ˆ
! ! !, , ,
2 !
p q r
e
p q rd d p q r
p q r
(25)
Seja e
i um segmento do contorno de e . Definimos o comprimento de e
i por iL ,
2 2 2 2 2 2
1 3 3 2 1 1 3 2 2, , L a b L a b L a b (26)
A fórmula para calcular a integral de linha sobre o segmento e
i é dada por
ˆ1 2 1 2 0
! ! , ,
1 !e ei i
p q p q
i i
p qw w ds L ds L p q
p q
(27)
Procedendo de forma análoga, como no caso global, chegamos ao seguinte sistema linear de
EDO
u u
u 0 0
e e e e e
e
M K f
(28)
onde 0ef é 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”
3 3
e e
ij i je
M m c w w dA
(29)
74
1 2 e e
i j i j i je e e
c w w dxdy cJ w w d d cJ d d
Então,
3 3
e e e
ij ijM m cJ I
(30)
onde 1 2 ij i j
eI d d é calculada pela formula dada em Eq. (25), isto é,
1 2
1se
12
1se
24
ij i je
i j
I d d
i j
(31)
3.3.2. Matriz “Rigidez”
A Matriz “Rigidez” elementar é dada por
1 2e e eK K K (32)
Cálculo de 1eK
1 1
3 3 e e
ij xi xj yi yje
K k k w w w w dA
(33)
então
1 1 2 e
ij xi xj yi yj xi xj yi yj ij ije e e
k k w w w w dA k w w dA w w dA k I I (34)
onde,
1 ij xi xje
I w w dA (35)
e
2 ij yi yje
I w w dA (36)
Calculando 1
ijI temos que
1
1 2 3 1 2 3ˆ1 2 3 1 2 3
1
j j ji i iij xi xj ee e
I w w dA b b b b b b dAJ
(37)
ˆ
ˆ
1
1
jik se e
k s
k ik s jse e
b b dAJ
b b dAJ
onde ij é o delta de Kronecker. Observe que para k i e s j , temos
75
ˆ1
ˆ
1 1 1
2
e
ij i j i j i je e eeI bb dA bb A bb
J J J (38)
Por tanto,
1 1
2ij i je
I bbJ
(39)
De forma análoga temos para
ˆ2
ˆ
1 1 1
2
e
ij i j i j i je e eeI a a dA a a A a a
J J J (40)
ou seja,
2 1
2ij i je
I a aJ
(41)
Substituindo equações (39) e (41) em Eq. (34),
1 1
3 3
1 1
2 2 2
e e
ij i j i j i j i je e e
kK k k bb a a bb a a
J J J
(42)
Cálculo de 2eK
Com ajuda da Eq. (12), calculemos a normal exterior em cada lado do elemento e , veja Figura
4:
1 2 1 2 1 2 1 3 3 1 3 3
2 3 2 3 2 3 2 1 1 2 1 1
3 1 3 1 3 1 3 2 2 3 2 2
, , n ,
, , n ,
, , n ,
v n n x x y y a b b a
v n n x x y y a b b a
v n n x x y y a b b a
(43)
Como as normais que precisamos devem ser unitárias, então
3 3
11 12
1 1
n , n ,n n
b a (44)
1 121 22
2 2
n , n ,n n
b a (45)
2 231 32
3 3
n , n ,n n
b a (46)
76
Figura 3. Fronteira de um elemento arbitrário.
Seja n o vetor normal à curva 3
1
e e
i
i
. Então
2 2 n
e
e e i
ij j ijN N
wK k k w ds kII
(47)
onde
1 2 3
1 2 3
nn
n n n
e e
e e e
i
ij j i j
i j i j i j
wII w ds w w ds
w w ds w w ds w w ds
(48)
Para calcular a matriz ijII devemos saber qual é o lado do elemento que faz parte da fronteira
do domínio .
Suponhamos que 1
e , então
1
11 12n + ne
aij
ij xi yi j
I
II w w w ds
Calculo de a
ijI
ˆ1 1
ˆ1
111 12 11 12
1 111 12 11 12
111 12
1 11 1 12 1 11 1 12
12 11 2 12 2 11 2 12
n + n n + n
1!n n n n
1 1 !
n n2
n n n n 0
n n n n 02
0 0 0
e e
e
a
ij xi yi j i i je
i i j i ie e
i ie
e
LI w w w ds b a ds
J
L Lb a ds b a
J J
Lb a
J
b a b aL
b a b aJ
77
Suponhamos que 2
e , então
2
21 22n + ne
bij
ij xi yi j
I
II w w w ds
Calculo de b
ijI
2
221 22 21 22
22 21 2 22 2 21 2 22
3 21 3 22 3 21 3 22
n + n n n2
0 0 0
0 n n n n2
0 n n n n
e
b
ij xi yi j i ie
e
LI w w w ds b a
J
Lb a b a
Jb a b a
Suponhamos que 3
e , então
3
31 32n + ne
cij
ij xi yi j
I
II w w w ds
Calculo de c
ijI
2
3
31 32 31 32
1 31 1 32 1 31 1 32
3
3 31 3 32 3 31 3 32
n + n n n2
n n 0 n n
0 0 02
n n 0 n n
e
c
ij xi yi j i ie
e
LI w w w ds b a
J
b a b aL
Jb a b a
3.4 Domínio Discretizado
Na Figura 4, observamos uma malha para . Esta malha nos fornece os seguintes dados:
12, 11, 3, 8e n i c n iN N N N N N (49)
Figura 4. Malha para o domínio.
78
onde eN é o Número de elementos, nN é o Número de nós,
iN é o Número de incógnitas e
cN é 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 9i , 3j e 7k definimos a seguintes matriz elementar e o vetor
derivada elementar, respectivamente,
(1) (1) (1) (1)
99 93 97 9
(1) (1) (1) (1) (1) (1) (1)
39 33 37 33 3(1) (1) (1) (1)
79 73 77 7
u
, u u
u
ij
m m m
m m m m m
m m m
(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
1(1) 1(1) 1(1) 2(1) 2(1) 2(1)
99 93 97 99 93 97
1(1) 1(1) 1(1) 1(1) 1(1) 2(1) 2(1) 2(1) 2(1) 2(1)
39 33 37 39 33 373 3 3 31(1) 1(1) 1(1) 2(1) 2(1) 2(1)
79 73 77 79 73 77
, ij ij
k k k k k k
k k k k k k k k k k
k k k k k k
(51)
e
(1) (1)
9 9
(1) (1) (1) (1)
3 3
(1) (1)
7 7
u
u u ,
u
f
f f
f
(52)
Do mesmo modo faz para todos os elementos da malha.
3.5 Ensamblagem
Consiste em montar todas as eN matrizes elementares numa única matriz chamada de matriz
global cuja ordem é n nN N .
3.5.1 Montagem
Uma maneira de visualizar todas as matrizes elementares numa única matriz é como
mostramos a continuação
79
1(1)
2(2)
( )
1 11(1) 2(1)
2 21(2) 2(2)
1( ) 2( )
u0 0
0 0 u
0 0 u
u0 0 0 0
0 0 0 0 u
0 0 0 0 u
ee
e ee e
N N
N N N N
m
m
m
fk k
k k f
k k f
(53)
3.5.2 Obtenção da Matriz e Vetor Global
O gerador de malha nos fornece um arquivo de dados, assim como mostra a Tabela 1.
Tabela 1
Elemento i j k
1 9 3 7
2 6 3 9
3 5 2 9
4 9 2 6
5 10 1 5
6 8 1 10
7 5 9 10
8 7 4 11
9 11 4 8
10 11 9 7
11 11 8 10
12 9 11 10
Como o 12eN , 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 2 1 2 3 4 5 6 7 8 9 10 11
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
80
3 1 2 3 4 5 6 7 8 9 10 11 4 1 2 3 4 5 6 7 8 9 10 11
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
5 1 2 3 4 5 6 7 8 9 10 11 6 1 2 3 4 5 6 7 8 9 10 11
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
7 1 2 3 4 5 6 7 8 9 10 11 8 1 2 3 4 5 6 7 8 9 10 11
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
9 1 2 3 4 5 6 7 8 9 10 11 10 1 2 3 4 5 6 7 8 9 10 11
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
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 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
81
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 6 7 8 9 10 11
1 511+611 0 0 0 515 0 0 618 0 51,10+61,10 0
2 0 322+422 0 0 325 426 0 0 329+429 0 0
3 0 0 133+233 0 0 236 137 0 139+239 0 0
4 0 0 0 844+944 0 0 847 948 0 0 84,11+94,11
5 551 352 0 0 355+555 +755
0 0 0 359+759 55,10+75,10 0
6 0 462 263 0 0 266+466 0 0 269+469 0 0
7 0 0 173 874 0 0 177+877
+1077
0 179+1079 0 87,11+107,11
8 681 0 0 984 0 0 0 688+988 +1188
0 68,10+118,10 98,11+118,11
9 0 392+492 193+293 0 395+795 296+496 197+1097 0 199+299
+399+499
+799+1099
+1299
79,10+129,10 109,11+129,11
10 510,1+610,1 0 0 0 510,5
+710,5
0 0 610,8
+1110,8
710,9+ 1210,9
510,10+ 610,10+ 710,10+ 1110,10+ 1210,10
1110,11
+1210,11
11 0 0 0 811,4+ 911,4
0 0 811,7+ 1011,7
911,8+ 1111,8
1011,9+ 1211,9
1111,10
1211,10
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 2 3 4 5 6 7 8 9 10 11 12
1
2
3
4
5
6
7
8
9
10
11
82
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
11 15 18 1 10
22 25 26 29
33 36 37 39
44 47 48 4 11
51 52 55 59 5 10
62 63 66 69
73 74 77 79 7 10
81 84 88 8 10 8 11
92 93 95 96 97 99 9 10 9 1
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0
m m m m
m m m m
m m m m
m m m m
m m m m m
m m m m
m m m m m
m m m m m
m m m m m m m m
1
2
3
4
5
6
7
8
1 9
10 1 10 5 10 8 10 9 10 10 10 11 10
11 4 11 7 11 8 11 9 11 10 11 11 11
u
u
u
u
u
u
u
u
u
0 0 0 0 0 u
0 0 0 0 0 u
m m m m m m
m m m m m m
+
11 15 18 1 10
22 25 26 29
33 36 37 39
44 47 48 4 11
51 52 55 59 5 10
62 63 66 69
73 74 77 79 7 10
81 84 88 8 10 8 11
92 93 95 96 97 99 9 10 9 1
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0
k k k k
k k k k
k k k k
k k k k
k k k k k
k k k k
k k k k k
k k k k k
k k k k k k k k
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
1 9 9
10 1 10 5 10 8 10 9 10 10 10 11 10 10
11 4 11 7 11 8 11 9 11 10 11 11 11 11
u
u
u
u
u
u
u
u
u
0 0 0 0 0 u
0 0 0 0 0 u
f
f
f
f
f
f
f
f
f
k k k k k k f
k k k k k k f
1 2u ug g g g g gM K K f (54)
83
u ug g g g gM K f (55)
onde 0gf e 1 2g g gK K K .
3.5.3 Obtenção da Matriz e Vetor a Resolver
Análise na matriz M . Usando as condições de contorno
11m 0 0 0 15m 0 0 18m 0 1 10m 0
0 22m 0 0 25m 26m 0 0 29m 0 0
0 0 33m 0 0 36m 37m 0 39m 0 0
0 0 0 44m 0 0 47m 48m 0 0 4 11m
51m 52m 0 0 55m 0 0 0 59m 5 10m 0
0 62m 63m 0 0 66m 0 0 69m 0 0
0 0 73m 74m 0 0 77m 0 79m 0 7 10m
81m 0 0 84m 0 0 0 88m 0 8 10m 8 11m
92 93 95 96 97
10 1 10 5 1
1
99 9 10 9 11
10 9 10 10 10 11
11 9 11 1
0 8
11 4 11 7 11 8 0 11 11
0 0 0
0 0 0 0 0
0
u
0 0 0 0
m m m m m
m m m
m
m m m
m m
m m m
m m m
2u
3u
4u
5u
6u
7u
8u
9
10
11
u
u
u
De forma análoga para a matriz K . Desse modo temos
99 9 10 9 11 9 99 9 10 9 11 9 9
10 9 10 10 10 11 10 10 9 10 10 10 11 10 10
11 9 11 10 11 11 11 11 9 11 10 11 11 11 11
u u
u u
u u
i iN N
m m m k k k f
m m m k k k f
m m m k k k f
92 93 95 96 97
10 1 10 5 10 8
11 4 11 7 11 8
92 93 95
1
2
3
4
5
6
7
9
8
u
u
u
u
u
u
u
u
0 0 0
0 0 0 0 0
0 0 0 0 0
0 0
i cM N N
m m m m m
m m m
m m m
k k k k
6 97
10 1 10 5 10 8
11 4 11 7 11 8
1
2
3
4
5
6
7
8
u
0
0 0 0 0 0
0
u
u
u
u
u
u
0 0
u
0 0
i cK N N
k
k k k
k k k
ou equivalentemente
84
u u u uM K f M K (56)
onde 0f . Assim,
u u fM K (57)
onde f u uf M K
3.6 Solução do Sistema de EDO
Malha temporal
O tamanho de passo temporal é definido por:
T
tN
(58)
e os nós temporais, por
0 , 1,nt t n t n N (59)
onde 0 0t é o tempo inicial.
Condição Inicial
0 0u t (60)
Regra do Médio ponto Generalizado
De Eq. (60), temos que
0 0u (61)
Faça para 1,n N
1 1 1 u f 1 fu n nn nM t K M t K t (62)
Observe que Eq. (62) representa um Sistema de Equações Lineares, veja LEWIS et. al (2004).
Ax b (63)
onde
1 1
u
1 u f 1 f
n
n n n
A M t K
x
b M t K t
(64)
Para resolver o Sistema Linear em Eq. (63) usamos o método SOR. Lembrando que 0 1 .
Neste trabalho usamos 0,5 .
85
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
Tipo de Malha eN nN
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 481iN e da malha não
estruturada foi 377iN . 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 iN Tempo CPU em [s]
Estruturada 481 22
Não Estruturada 377 14
86
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 35 10t e se atingiu o estado estacionário
em 2 [s] com 400N . 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.
0 [ ]t s
87
0,5 [ ]t s
1 [ ]t s
2 [ ]t 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.
88
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
89
Recommended
