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 104
Ouro
1, 27 104
Cobre
386,0
8,954
383,1
1,1234 104
Alumínio
8, 418 105
Aço
2,3 105
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)
i1
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 eˆ centrado na origem, tal como mostra a Figura 3.
(17)
73
Figura 2. Transformação do triangulo e no triangulo eˆ .
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 eˆ é 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
33
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
33
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
33
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 33 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
i1
K 2e kij2e
k
NN
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 33 (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 33 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) 33 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 1111 . 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 1111 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 111 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 uK 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 nt ,
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 n1 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 103 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
