ANÁLISE NUMÉRICA POR DIFERENÇAS FINITAS DA EQUAÇÃO DE CONDUÇÃO DE CALOR BIDIMENSIONAL EM PAREDES PLANAS

Prévia do material em texto

1
ANA´LISE NUME´RICA POR
DIFERENC¸AS FINITAS DA
EQUAC¸A˜O DE CONDUC¸A˜O DE
CALOR BIDIMENSIONAL EM
PAREDES PLANAS
Bruno Ferreira Couto
brunofcouto1@gmail.com
Johnathan Batista dos Santos
johnathan.batista.23@gmail.com
Prof. Dr. Felipe Mariano
2
Universidade Federal de Goia´s
Escola de Engenharia Ele´trica, Mecaˆnica e de
Computac¸a˜o - EMC
Bacharelado em Engenharia Mecaˆnica
Disciplina: Transfereˆncia de Calor 1
Prof. Dr. Felipe Pamplona Mariano
Goiaˆnia, 16 de julho de 2016
Fonte: Times New Roman tamanho 11pt
Equac¸o˜es digitadas com o aux´ılio do software
MATHTYPE®
Gra´ficos e isotermas criadas via MATLAB®, em
formato vetorizado .eps
Alterac¸o˜es gra´ficas de matiz/saturac¸a˜o, curvas e to-
nalidade realizadas via ADOBE PHOTOSHOP®
Agradecimentos: Prof. Dr. Sigeo Kitatani -
EMC/UFG ; Matheus Gama
Editorado via TEXstudio no sistema LATEX.
Capa: isotermas geradas computacionalmente
para uma parede plana, bidimensional, com tem-
peraturas definidas no topo e na base e nulas nas
laterais (Cortesia: Matheus Gama - Dog&Tag).
Lista de S´ımbolos
α Difusividade te´rmica
∆m Incremento diferencial em m
∆t Incremento diferencial de tempo
∆x Incremento diferencial em x
∆y Incremento diferencial em y
ρ Densidade
C Calor espec´ıfico
k Condutividade te´rmica
L Comprimento total
m, i I´ndices nume´ricos em x
n, j I´ndices nume´ricos em y
T Temperatura
t Tempo
T∞ Temperatura longe da superf´ıcie
Tv, T0 Temperatura na superf´ıcie ou pro´xima dela
4
Lista de Figuras
1. Exemplo de modelagem matema´tica [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2. Representac¸a˜o de uma malha nodal [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3. Representac¸a˜o de isoterma [2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
4. Representac¸a˜o de uma malha nodal [2] . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
5. Exemplo de simulac¸o˜es nume´ricas com diferenc¸as finitas . . . . . . . . . . . . . . . . . 12
6. Representac¸a˜o da malha nodal adotada [1] . . . . . . . . . . . . . . . . . . . . . . . . . 13
7. Condic¸oes geome´tricas [2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
8. Isotermas obtidas (em °C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
9. Isotermas para passos de tempo menores (em °C) . . . . . . . . . . . . . . . . . . . . . 15
10. Tempo para a entrada em regime permanente . . . . . . . . . . . . . . . . . . . . . . . 16
11. Isotermas geradas da simulac¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
12. Tempo de resfriamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Lista de Tabelas
1. Propriedades do n´ıquel (T=300K) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2. Paraˆmetros da simulac¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3. Resultados da simulac¸a˜o nume´rica e anal´ıtica . . . . . . . . . . . . . . . . . . . . . . . 16
4. Condic¸o˜es de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
5. Paraˆmetros das condic¸o˜es de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6. Resultados da simulac¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5
Suma´rio
1. Introduc¸a˜o 7
1.1. Me´todos nume´ricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2. O Me´todo das Diferenc¸as Finitas 7
2.1. Modelagem matema´tica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2. Formulac¸a˜o por diferenc¸as finitas das equac¸o˜es diferenciais . . . . . . . . . . . . . . . . 8
2.3. Condic¸o˜es de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4. Formulac¸a˜o nume´rica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.5. Considerac¸o˜es gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3. Problema´tica 12
3.1. Primeiro problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.1.1. Abordagem anal´ıtica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.1.2. Abordagem nume´rica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2. Segundo problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4. Considerac¸o˜es Finais 19
Apeˆndices 21
A. Algoritmo do me´todo anal´ıtico 21
B. Algoritmo do me´todo nume´rico 22
C. Algoritmo do me´todo nume´rico do segundo problema 24
I´ndice Remissivo 25
6
1. Introduc¸a˜o
Os problemas de conduc¸a˜o de calor sa˜o regidos por uma equac¸a˜o diferencial que, atrave´s de um
balanc¸o de energia tridimensional em um elemento do volume de controle, relaciona a temperatura em
cada ponto do volume de controle T (x, y, z), com a gerac¸a˜o de energia e a energia acumulada dentro
dele mesmo. Todos os problemas de conduc¸a˜o de calor resolvidos analiticamente, em maior ou menor
grau, se resumem na resoluc¸a˜o da equac¸a˜o diferencial atrave´s de determinadas hipo´teses simplificadoras
(como conduc¸a˜o unidimensional ou bidimensional, gerac¸a˜o de energia nula, condutividade constante) e
condic¸o˜es de contorno espec´ıficas; o resultado sa˜o func¸o˜es da temperatura para cada ponto do volume
de controle.
Contudo, na maioria dos problemas pra´ticos de conduc¸a˜o de calor, geometrias complicadas, condic¸o˜es
de contorno complexas e propriedades intr´ınsecas varia´veis esta˜o envolvidas, dificultando uma soluc¸a˜o
anal´ıtica. Nesse caso, me´todos nume´ricos sa˜o as ferramentas mais simples e precisas de resoluc¸a˜o do
problema.
Quando se trata de ana´lises por me´todos nume´ricos, softwares de matema´tica alge´brica como o
Matlab1 sa˜o extremamente u´teis. Atrave´s dele, um algoritmo e´ escrito de forma a abarcar o problema
a ser resolvido e suas condic¸o˜es de contorno e, de maneira iterativa, os ca´lculos sa˜o repetidos milho˜es
de vezes ate´ se atingir a precisa˜o desejada; ale´m disso, gra´ficos do problema podem ser gerados de
maneira simples, na mesma interface do software.
1.1. Me´todos nume´ricos
Me´todos nume´ricos consistem na substituic¸a˜o da equac¸a˜o diferencial da conduc¸a˜o de calor por um
conjunto de n equac¸o˜es alge´bricas para temperaturas desconhecidas, em n pontos escolhidos estrate-
gicamente no meio. A soluc¸a˜o simultaˆnea dessas equac¸o˜es resulta nos valores da temperatura nesses
pontos discretos
Existem va´rios me´todos nume´ricos diferentes que podem ser empregados em problemas de trans-
fereˆncia de calor; dentre eles destacam-se o me´todo dos elementos finitos, o me´todo dos elementos
de contorno, o me´todo do balanc¸o de energia e o me´todo das diferenc¸as finitas. Este u´ltimo, que da´
t´ıtulo ao trabalho, se baseia no ja´ familiar balanc¸o de energia para volumes de controle e, dessa forma,
fornece um melhor entendimento dos problemas f´ısicos.
Vale ressaltar que cada me´todo nume´rico possui suas vantagens e desvantagens, cabendo ao enge-
nheiro se aprofundar no tema e concluir qual sera´ o me´todo mais vantajoso [1].
2. O Me´todo das Diferenc¸as Finitas
Como ja´ citado, cada me´todo possui suas vantagens e desvantagens. As considerac¸o˜es que se seguem
sa˜o sobre as peculiaridades do me´todo das diferenc¸as finitas.
2.1. Modelagem matema´tica
Quando se trata da ana´lise matema´tica de algum problema f´ısico, e´ necessa´rio uma modelagem, isto
e´, uma expressa˜o matema´tica que o represente na forma de nu´meros. Isso implica que, quanto melhor
for a modelagem, melhor sera´ a representac¸a˜o e o entendimento do fenoˆmeno f´ısico envolvido. Assim,
quanto mais simplificada for a modelagem, mais fa´cil sera´ a resoluc¸a˜o matema´tica; contudo, alguma
precisa˜o sera´ sacrificada naana´lise do problema f´ısico. De fato, existe uma tendeˆncia acentuada
nos engenheiros de querer sempre simplificar de mais os problemas para soluc¸o˜es anal´ıticas, como a
considerac¸a˜o da condutividade te´rmica constante (que na realidade na˜o e´) e a condic¸a˜o de contorno
1Matlab e´ um software registrado, de uso comercial.
7
da radiac¸a˜o que quase sempre e´ ignorada. Essa tendeˆncia introduz erros que, conforme seja o projeto,
tornam inadmiss´ıvel sua execuc¸a˜o via me´todos anal´ıticos. Com isso, a ana´lise por me´todos nume´ricos
se torna quase uma obrigac¸a˜o.
Um exemplo de modelagem matema´tica e´ dado na
figura (1), onde um corpo oval e´ modelado analiti-
camente como uma esfera, e numericamente como
um ovo. Note que a precisa˜o sera´ maior na ana´lise
nume´rica devido a geometria mais fiel.
Alguns problemas podem ser resolvidos analiti-
camente, pore´m a modelagem matema´tica pode
se tornar ta˜o complexa que o esforc¸o se mostra
invia´vel. Solucionar esses problemas podem envol-
ver equac¸o˜es diferenciais parciais, func¸o˜es de Bes-
sel e de Legendre, onde o n´ıvel de dificuldade foge
a` graduac¸a˜o. Nesse cena´rio, a ana´lise nume´rica se
mostra extremamente convidativa.
Vale ressaltar que os benef´ıcios do me´todo nume´rico
na˜o excluem, de maneira alguma, a soluc¸a˜o
anal´ıtica do problema. A noc¸a˜o e o bom senso com
relac¸a˜o aos fenoˆmenos f´ısicos sa˜o adquiridos justa-
mente da ana´lise dos problemas, de modo que o
nume´rico e o anal´ıtico caminham juntos, e de ma˜os
dadas. A crenc¸a de que o me´todo nume´rico in-
validou o me´todo anal´ıtico e´ como acreditar que as
calculadoras tornaram desnecessa´rio que as pessoas
aprendam aritme´tica na escola.
Figura 1: Exemplo de modelagem matema´tica [1].
2.2. Formulac¸a˜o por diferenc¸as finitas das equac¸o˜es diferenciais
Do ca´lculo diferencial e´ tido que a derivada primeira de uma func¸a˜o f(x), e´ dada como
df(x)
dx
= lim
∆x→0
∆f
∆x
= lim
∆x→0
f(x+ ∆x)− f(x)
∆x
, (1)
que e´ a inclinac¸a˜o da linha tangente a` curva no ponto x. Se o limite quando ∆x tende a 0 na˜o for
tomado, a seguinte aproximac¸a˜o para a derivada pode ser adotada
df(x)
dx
∣∣∣∣
x
=
f(x+ ∆x)− f(x)
∆x
. (2)
De modo semelhante a` equac¸a˜o(1) e a` equac¸a˜o(2), a derivada segunda de f(x), aproximada, pode ser
representada por
d2f(x)
dx2
∣∣∣∣
x
=
f(x+ ∆x)− 2f(x) + f(x−∆x)
∆x2
. (3)
As expresso˜es apresentadas nas equac¸o˜es (2) e (3) sa˜o, respectivamente, as representac¸o˜es das deri-
vadas primeira e segunda em diferenc¸as finitas da func¸a˜o f(x).
8
Seja a equac¸a˜o da conduc¸a˜o de calor para o caso bidimensional dada por
ρC
∂T
∂t
= k
(
∂2T
∂x2
+
∂2T
∂y2
)
, (4)
onde a difusividade, α, e´ dada por
α =
k
ρC
⇒ ∂T
∂t
= α
(
∂2T
∂x2
+
∂2T
∂y2
)
. (5)
A equac¸a˜o(5) e´ tambe´m conhecida como equac¸a˜o da conduc¸a˜o bidimensional em regime transiente.
Considerando, agora, uma parede plana de altura y = n∆y e largura x = m∆x, onde m e n sa˜o
subdiviso˜es de mesmo tamanho. Cada ponto de coordenada (m,n) e´ um no´ ou ponto nodal; o
conjunto de pontos nodais formam a malha de pontos onde a func¸a˜o sera´ calculada.
Figura 2: Representac¸a˜o de uma malha nodal [2].
Note, pela figura(2), que o comprimento total da
parede sera´ Ly = n∆y e Lx = m∆x. A maneira
como a malha e´ dividida esta´ a cargo do enge-
nheiro. Obviamente, quanto mais no´s tiver a ma-
lha, mais preciso sera´ a ana´lise; pore´m, sera´ exigido
maior esforc¸o computacional para a realizac¸a˜o dos
ca´lculos. Geometricamente, a malha deve represen-
tar da maneira mais fiel poss´ıvel o objeto real; caso
contra´rio, recorrera´ em algum erro de modelagem,
discutido na sec¸a˜o (2.1). Geometrias esfe´ricas ou
cil´ındricas sa˜o mais facilmente modeladas em co-
ordenadas esfe´ricas e cil´ındricas, respectivamente.
A modelagem para tais casos foge ao escopo desse
trabalho.
As derivadas de primeira e segunda ordem, em func¸a˜o da malha nodal, ficara˜o como
∂T
∂t
∣∣∣∣
m,n
=
α
∆x2
(Tm+∆x,n − 2Tm,n + Tm−∆x,n) + ...+ α
∆y2
(Tm,n+∆y − 2Tm,n + Tm,n−∆y) (6)
Figura 3: Representac¸a˜o de isoterma [2] Figura 4: Representac¸a˜o de uma malha nodal [2]
9
2.3. Condic¸o˜es de contorno
As relac¸o˜es apresentadas sa˜o para a obtenc¸a˜o da expressa˜o da temperatura para os no´s internos da
malha. Essas relac¸o˜es, pore´m, na˜o sa˜o va´lidas paras os pontos do per´ımetro da malha; elas exigem a
presenc¸a de no´s de ambos os lados do no´ em ana´lise e, na fronteira, em pelo menos um dos lados, na˜o
havera´ no´ adjacente [1].
Sendo assim, e´ necessa´ria a adoc¸a˜o de certos paraˆmetros f´ısicos nas fronteiras para que a modela-
gem matema´tica possa fazer sentido: esses paraˆmetros sa˜o conhecidos como condic¸o˜es de contorno.
Algumas condic¸o˜es de contorno mais recorrentes sa˜o:
1 Temperatura especificada: condic¸a˜o mais simples, estipula que cada lado da fronteira da
malha estara´ fixamente a uma dada temperatura;
2 Fluxo de calor especificado: considerando a transfereˆncia de calor ocorrendo para dentro da
malha (fluxo positivo), sem gerac¸a˜o interna de energia
q˙0A+ kA
(
T1 − T0
∆x
)
= 0 (7)
3 Contorno isolado (fluxo de calor nulo): considerando algum lado da fronteira como adiaba´tico
kA
(
T1 − T0
∆x
)
= 0 (8)
4 Convecc¸a˜o: ocorrendo convecc¸a˜o nas fronteiras
hA (T∞ − T0) + kA
(
T1 − T0
∆x
)
= 0 (9)
5 Radiac¸a˜o: ocorrendo trocas por radiac¸a˜o nas fronteiras
εσA
(
Tv
4 − T 4∞
)
+ kA
(
T1 − T0
∆x
)
= 0 (10)
6 Convecc¸a˜o e radiac¸a˜o combinados: ocorrendo trocas por radiac¸a˜o e convecc¸a˜o nas fronteiras
hA (T∞ − T0) + εσA
(
Tv
4 − T 4∞
)
+ kA
(
T1 − T0
∆x
)
= 0 (11)
2.4. Formulac¸a˜o nume´rica
A formulac¸a˜o nume´rica consiste, essencialmente, em se declarar varia´veis (vetores) dentro de algum
lac¸o de repetic¸a˜o, de modo a reproduzir a equac¸a˜o(6). Observe o lac¸o de repetic¸a˜o utilizado para se
calcular o tempo total para a entrada em regime permanente.
10
1 L=1; %comprimento
2 W=1; %altura
3 a=23/1000000; %difusividade termica
4
5 nx=20; %numero de intervalos em x
6 ny=20; %numero de intervalos em y
7 cont=1;
8 dx=L/nx;
9 dy=W/ny;
10
11 for h=3600:900:28800
12 tf=h; %tempo final
13 dt=0.25*(dxˆ2)/a; %intervalo de tempo
14 nt=tf/dt; %numero de passos no tempo
15 x=[0:dx:L];
16 y=[0:dy:W];
17 T1=0; %Temperatura no contorno
18 T2=1; %Temperatura no contorno
19
20 for n=1:nt
21 T0=T;
22 for i=2:nx
23 for j=2:ny
24 T(i,j) = T0(i,j)+dt*a*((T0(i+1,j)-2*T0(i,j)+T0(i-1,j))/(dxˆ2)+(T0(i,j+1)-
25 2*T0(i,j)+T0(i,j-1))/(dyˆ2));
26 end
27 end
28 end
As condic¸o˜es de contorno e inicial devem ser especificadas previamente, de modo a preencher as res-
pectivas posic¸o˜es dentro do vetor. Nesse caso, a condic¸a˜o de contorno e´ a de temperatura especificada
nas laterais da chapa e temperatura me´dia nas quinas.
1 %condicao inicial
2 T=ones(nx+1,ny+1)*T1;
3 %condicoes de contorno nas quinas
4 T(1,1)=(T1+T1)/2;
5 T(nx+1,1)=(T1+T1)/2;
6 T(1,ny+1)=(T1+T2)/2;
7 T(nx+1,ny+1)=(T1+T2)/2;
8 %condicoes de contorno
9 for i=2:nx
10 T(i,1)=T1; %parede inferior
11 T(i,ny+1)=T2; %parede superior
12 end
13 for j=2:ny
14 T(1,j)=T1; %lateral esquerda
15 T(nx+1,j)=T1; %lateral direita
16 end
No caso do me´todo das diferenc¸as finitas, o crite´rio de estabilidade (ou convergeˆncia) e´ que
∆t ≤ 1
2
∆x2
α
. (12)
Assim, para que o me´todo convirja e´ necessa´rio que o passo de tempo, ∆t, seja menor que metade da
raza˜o entre o quadrado do passo da malha e a difusividade te´rmica do material.
11
2.5. Considerac¸o˜es gerais
A formulac¸a˜o de diferenc¸as finitas resulta em N equac¸o˜es alge´bricas, com N temperaturas nodais
desconhecidas que precisam ser resolvidas simultanemante. Em termosde me´todos nume´ricos compu-
tacionais, ha´ duas abordagens comuns: a abordagem direta e a abordagem iterativa. A abordagem
direta consiste em um nu´mero fixo e bem definido de passos que resultam na soluc¸a˜o da forma sis-
tema´tica. A abordagem iterativa consiste numa estimativa inicial, que e´ refinada a cada iterac¸a˜o ate´
que determinado crite´rio de convergeˆncia ou grau de precisa˜o seja atingido.
Me´todos diretos envolvem grande esforc¸o computacional (muita memo´ria) e sa˜o indicados para
sistemas com menos equac¸o˜es; me´todos iterativos sa˜o mais leves computacionalmente, pore´m sa˜o mais
suscet´ıveis a erros de convergeˆncia [1].
(a) Radiac¸a˜o (b) Superf´ıcie sendo aquecida
(c) Circuito eletroˆnico (d) Superf´ıcie em aquecimento
Figura 5: Exemplo de simulac¸o˜es nume´ricas com diferenc¸as finitas
3. Problema´tica
3.1. Primeiro problema
A ana´lise nume´rica e´ modelada para uma parede plana, bidimensional, de dimenso˜es (1 × 1)m2,
com condutividade te´rmica constante e gerac¸a˜o de calor nula. O material adotado sera´ o Nı´quel e
suas propriedades sa˜o dadas na tabela(1). Note que as propriedades sa˜o constantes e referentes a
temperatura de 300K.
A distribuic¸a˜o nodal na malha esta´ de acordo com a figura(6) e como ja´ explanado nas sec¸o˜es (2.1)
e (2.2). As condic¸o˜es de contorno sa˜o de temperatura especificada nas paredes inferior e laterais, como
12
T2 = 1C, e na parede superior como T1 = 0C. A expressa˜o nume´rica utilizada e´ a equivalente da
equac¸a˜o(6), adaptada para o Matlab.
Tabela 1: Propriedades do n´ıquel (T=300K)
ρ(kg/m3) C(J/kg.K) k(W/m.K) α.10−6(m2/s)
8900 444 90.7 23.0
Figura 6: Representac¸a˜o da malha nodal adotada [1]
3.1.1. Abordagem anal´ıtica
A abordagem anal´ıtica parte do pressuposto da bi-
dimensionalidade da geometria em questa˜o, onde
3 lados esta˜o a uma mesma temperatura, fixa, e
o quarto lado esta´ a uma temperatura, tambe´m
fixa, diferente daquelas presentes nos lados ad-
jacentes, figura(7).O me´todo e´ proveniente da
soluc¸a˜o anal´ıtica da equac¸a˜o diferencial (5), pela
te´cnica de separac¸a˜o de varia´veis. O resultado e´
a equac¸a˜o(13), que e´ uma se´rie convergente para
quaisquer x e y. Note que a equac¸a˜o (13) fornece o
fluxo de calor; a temperatura no ponto deve ser cal-
culada atrave´s da substituic¸a˜o pela equac¸a˜o (14). Figura 7: Condic¸oes geome´tricas [2]
θ(x, y) =
2
pi
∞∑
n=1
(−1)n+1 + 1
n
sin
(npix
L
) sinh (npiyL )
sinh
(
npiW
L
) , (13)
13
θ ≡ T − T1
T2 − T1 . (14)
O algoritmo para o ca´lculo anal´ıtico e´ dado a seguir. Sobre as varia´veis, i e´ o ı´ndice de um vetor
s(i) que armazena proviso´riamente o resultado da conta, a cada iterac¸a˜o; somatorio soma os valores
de s(i). Teta e´ o fluxo te´rmico, e para se encontrar a temperatura no ponto de interesesse, T , usa-se
a equac¸a˜o (14).
1 W=1; %Largura
2 L=1; %Altura
3 n=2000;
4 T2=1;
5 T1=0;
6 x=0.25;
7 y=0.25;
8 i=1;
9 S=0;
10 somatorio=0;
11 while i<=21
12
13 s(i)=((((-1)ˆ(i+1))+1)/i)*sin(i*pi*x/L)*(sinh(i*pi*y/L)/sinh(i*pi*W/L));
14 somatorio=somatorio+s(i);
15 i=i+1;
16 end
17 teta=(2/pi)*somatorio;
18 T=((T2-T1)*teta)+T1;
No apeˆndice (A) e´ apresentado o algoritmo completo usado para o me´todo anal´ıtico. O nu´mero n
de iterac¸o˜es no somato´rio e´ relativo apenas a convergeˆncia do me´todo, na˜o interferindo na consisteˆncia
do resultado, uma vez atingida a convergeˆncia.
3.1.2. Abordagem nume´rica
A abordagem nume´rica e´ baseada na equac¸a˜o (6). Os paraˆmetros da simulac¸a˜o nume´rica sa˜o
mostrados tabela (2). Repare que foram adotados espac¸amentos de 0.05m tanto em x quanto em y, o
que implica em um passo de tempo de ∆t igual a 21.17s. A figura (8) mostra as isotermas obtidas da
simulac¸a˜o.
Na figura(10) e´ poss´ıvel ver as curvas da temperatura em func¸a˜o do tempo nos pontos x = (L/4) =
0, 25m, y = (W/4) = 0, 25m e x = (L/2) = 0, 5m, y = (W/4) = 0, 5m; e´ percept´ıvel o tempo gasto ate´
que o processo de conduc¸a˜o de calor entre em regime permanente. No apeˆndice (B) e´ apresentado
o algoritmo completo usado nesse caso. E´ perfeitamente n´ıtido, atrave´s da ana´lise das isotermas,
Tabela 2: Paraˆmetros da simulac¸a˜o
∆x(m) ∆y(m) ∆t(s) α(m2/s) Nx Ny
0.05 0.05 21.17 23.10−6 20 20
que o fluxo de calor tende da maior temperatura (acima) em direc¸a˜o a`s menores temperaturas (faces
laterais e inferior). Devido a` caracter´ıstica da bidimensionalidade e da geometria do problema, todas
as isotermas sa˜o sime´tricas em relac¸a˜o ao centro da figura, verticalmente.
14
 
 
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figura 8: Isotermas obtidas (em °C)
Para passos de tempo menores, ∆t, as isotermas sa˜o multiplicadas e e´ poss´ıvel notar uma falha de
modelagem nas quinas inferiores da parede, conforme a figura(9) mostra. Isso e´ uma consequeˆncia das
condic¸o˜es de contorno adotadas para essas regio˜es, e e´ evidenciado porque o passo de tempo menor
aumenta a resoluc¸a˜o do gra´fico.
Com relac¸a˜o ao tempo para a entrada em regime permanente, de aproximadamente cinco horas
para x = y = 0.25m, e de cinco horas e trinta minutos para x = y = 0.5m, e´ um tempo condizente e
coerente com as dimenso˜es da parede (relativamente grandes) e a difusividade do material adotado.
Tambe´m e´ poss´ıvel perceber que a diferenc¸a de tempo para a entrada em regime permanente, nos
dois pontos de interesse, ocorre porque o ponto central (0.5/0.5(m)) esta´ mais pro´ximo da fonte de
calor, acima da placa; isso faz com que esse ponto entre em regime mais tardiamente do que o ponto
inferior (0.25/0.25(m)), mais abaixo na placa.
 
 
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
(a) Passo=0.0001s
 
 
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
(b) Passo=0.001s
Figura 9: Isotermas para passos de tempo menores (em °C)
15
1 2 3 4 5 6 7 8
0.17
0.18
0.19
0.2
0.21
0.22
0.23
0.24
0.25
TEMPO (H)
TE
M
PE
R
AT
UR
A 
(°C
)
(a) x=y=0.25m
1 2 3 4 5 6 7 8
0.03
0.035
0.04
0.045
0.05
0.055
0.06
0.065
0.07
0.075
TE
M
PE
R
AT
UR
A 
(°C
)
TEMPO (H)
(b) x=y=0.5m
Figura 10: Tempo para a entrada em regime permanente
Os resultados obtidos da simulac¸a˜o anal´ıtica e da simulac¸a˜o nume´rica culminaram nos dados apre-
sentados na tabela(3). Repare que o me´todo anal´ıtico na˜o e´ capaz de estimar o tempo para a entrada
em regime permanente do processo de troca de calor.
Tabela 3: Resultados da simulac¸a˜o nume´rica e anal´ıtica
Resultados anal´ıticos
No´ T (C) Tempo para a entrada em regime
x = y = 0.25m 0.06797 **
x = y = 0.5m 0.02500 **
Resultados nume´ricos
No´ T (C) Tempo para a entrada em regime
x = y = 0.25m 0.0681 5h
x = y = 0.5m 0.2499 5h30min
3.2. Segundo problema
O segundo problema consiste na mesma geometria e material do problema anterior, pore´m agora
com condic¸o˜es de contorno espec´ıficas em cada lado. Nesse caso, e´ necessa´rio encontrar o tempo para
que a chapa, a uma temperatura inicial de 300◦C, chegue aos 30◦C em seu ponto central.
As condic¸o˜es de contorno e seus respectivos paraˆmetros sa˜o apresentados nas tabelas (4) e (5). O
fluxo imposto tem o sentido de dentro para fora da placa. As propriedades intr´ınsecas do material
sa˜o as mesmas apresentadas na tabela (1).
Tabela 4: Condic¸o˜es de contorno
Parede esquerda Parede direita Parede superior Parede inferior
convecc¸a˜o ’C’ fluxo imposto ’Q’ temperatura constante ’T1’ parede adiaba´tica ’q0’16
Tabela 5: Paraˆmetros das condic¸o˜es de contorno
Q(W ) q0 h(W/m2K) T∞(C) T1(C)
10 ** 100 5°C 0°C
O algoritmo e´ apresentado abaixo. Estruturalmente, a maior diferenc¸a entre este algoritmo e o
algoritmo nume´rico do primeiro problema e´ a alterac¸a˜o das condic¸o˜es de contorno, juntamente com a
questa˜o de se parar a simulac¸a˜o para o tempo em que a temperatura de 30◦C for atingida no meio
da chapa. As demais diferenc¸as esta˜o relacionadas com comandos gra´ficos para melhor apresentar os
resultados. O algoritmo completo e´ dado no apeˆndice (C).
1 while erro>=2
2 tf=k; %tempo final
3 dt=0.25*(dxˆ2)/a; %intervalo de tempo
4 nt=tf/dt; %numero de passos no tempo
5 x=[0:dx:L];
6 y=[0:dy:W];
7 Ti=300; %Temperatura inicial
8 T1=0; %Temperatura no contorno
9 q1=0; %Fluxo na parede inferior
10 q2=-10; %Fluxo na lateral direita
11 %condicao inicial
12 T=ones(nx+2,ny+2)*Ti;
13 for n=1:nt %Laco do tempo
14 T0=T;
15 %condicoes de contorno
16 for i=1:nx+2
17 T(i,1) = T0(i,ny)+2*dy/dx*q1/k; %parede inferior (q0 condicao de contorno adiabatica)
18 end
19 for i=2:nx+1
20 T(i,ny+1)=T1
21 end
22 for j=2:ny+1
23 %lateral esquerda (condicao de contorno de conveccao com h =100 e tinf= 5)
24 T(1,j) = T0(3,j)-(2*dx*h/k)*(T0(2,j)-Tinf);
25 %lateral direita (Q fluxo de calor imposto igual a 10 W saindo da chapa)
26 T(nx+2,j)= T0(nx,j)+2*dy/dx*q2/k;
27 end
28 %Resolucao da equacao de calor no interior da chapa
29 for i=2:nx+1
30 for j=2:ny
31 T(i,j) = T0(i,j)+dt*a*((T0(i+1,j)-2*T0(i,j)+T0(i-1,j))/(dxˆ2)+(T0(i,j+1)-
32 2*T0(i,j)+T0(i,j-1))/(dyˆ2));
33 end
34 end
35 for i=1:nx+2
36 for j=1:ny+2
37 TTT(j,i) = T(i,j);
38 end
39 end
40 end
41 tempo(cont)=k/3600;
42 temp(cont)=TTT(11,11);
43 erro=temp(cont)-tmp final;
44 cont=cont+1;
45 k = k+900;
46 end
17
A figura (11) apresenta as isotermas obtidas da simulac¸a˜o. A figura (12) apresenta o gra´fico da
evoluc¸a˜o da temperatura no ponto central da chapa em func¸a˜o do tempo. Vale ressaltar que o
algoritmo demora, em me´dia, dois minutos para ser compilado, dependendo do compu-
tador utilizado.
largura da chapa de niquel
a
ltu
ra
 d
a 
ch
ap
a 
de
 n
iq
ue
l
 
 
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
5
10
15
20
25
30
Figura 11: Isotermas geradas da simulac¸a˜o
2 2.5 3 3.5 4
30
40
50
60
70
80
90
100
110
Tempo (h)
Te
m
pe
ra
tu
ra
 (º
C)
Figura 12: Tempo de resfriamento
O resultado da simulac¸a˜o e´ apresentado na tabela (6). Observe que o processo de transfereˆncia de
calor leva, aproximadamente, treˆs horas e quarenta e cinco minutos para atingir a temperatura de
18
30◦C, implicando em um erro de aproximadamente 1.4◦C.
Tabela 6: Resultados da simulac¸a˜o
Posic¸a˜o (m) Temperatura(°C) Tempo decorrido Erro absoluto (°C)
x=y=0.5 31.39 3h45min 1.390
4. Considerac¸o˜es Finais
No primeiro problema e´ poss´ıvel perceber nitidamente a proximidade entre o resultado anal´ıtico e o
resultado nume´rico. Isso corrobora a tese apresentada ainda na introduc¸a˜o, de que os procedimentos
anal´ıticos e nume´ricos sa˜o complementares e, por assim ser, devem sempre caminhar juntos.E´ tambe´m
nota´vel a consequeˆncia da condic¸a˜o de contorno adotada nas ”quinas”da chapa, que causou uma certa
deformidade nas isotermas, como mostrado na figura (9). Essa e´ uma limitac¸a˜o da modelagem bidi-
mensional, que na˜o e´ capaz de representar consistentemente um ponto (unidimensional por natureza);
o que implica que os resultados da simulac¸a˜o nume´rica na regia˜o das quinas na˜o devem ser tomados
como totalmente confia´veis. Na˜o ha´ crite´rios fixos para a adoc¸a˜o de condic¸o˜es de contorno nessas
regio˜es espec´ıficas, o que torna a sua implementac¸a˜o algo totalmente emp´ırico e a cargo do engenheiro
em questa˜o.
Sobre as caracter´ısticas desse problema, ainda com um esforc¸o computacional muito maior que o
primeiro, na˜o e´ poss´ıvel atingir exatamente a temperatura especificada no final da simulac¸a˜o; mesmo
com a malha nodal bem refinada. Vale ressaltar, tambe´m, que o erro calculado diminuiria caso o passo
de tempo fosse refinado.
No que tange a criac¸a˜o da malha nodal, e´ necessa´rio que o usua´rio fornec¸a uma quantidade de
diviso˜es par; o que implica numa quantidade de pontos nodais ı´mpar. Dessa forma, a malha sera´
sempre coincidente com o centroide da geometria. Caso contra´rio, sera´ necessa´rio um algoritmo para
interpolar o ponto central da malha, uma vez que ele estara´ rodeado por no´s, e na˜o estara´ em uma
posic¸a˜o determinada pela simulac¸a˜o. Ale´m disso, este me´todo traz consigo alguns problemas para
nu´meros de pontos pequenos na malha (deve ser maior que seis, por exemplo), dado que a soluc¸a˜o
via ca´lculo de me´dia na˜o traria grandes ganhos; a dificuldade para se programar esta func¸a˜o na˜o se
justifica, visto que e´ plenamente poss´ıvel que, dada as restric¸o˜es quanto ao nu´mero de intervalos (par
e mu´ltiplo de quatro) o programa calcule de modo satisfato´rio as temperaturas desejadas.
O crite´rio de estabilidade deve sempre receber atenc¸a˜o especial. Ele e´ determinante na adoc¸a˜o
do passo de tempo do lac¸o de repetic¸a˜o da simulac¸a˜o. Caso o seu valor na˜o esteja de acordo com
as caracter´ısticas da modelagem, o me´todo ou divergira´ ou entregara´ um resultado completamente
inconsistente; nesse caso, o engenheiro sera´ v´ıtima do famigerado erro de convergeˆncia nume´rica.
A qualidade dos resultados nume´ricos, tanto gra´ficos quanto anal´ıticos, depende do nu´mero de
iterac¸o˜es adotado; portanto, quanto mais potente for o computador utilizado, mais ra´pido o resultado
sera´ retornado. Esse trabalho na˜o se fez com o aux´ılio de computadores potentes, o que na˜o traz
resultados precisos ale´m da quarta casa decimal. Entretanto, pela simplicidade dos problemas em
pauta, o resultado e´ suficientemente bom.
Do ponto de vista do aprendizado, esse trabalho vem agregar a noc¸a˜o da necessidade de uma
modelagem matema´tica coerente com o problema real, e a ideia da complementariedade da abordagem
anal´ıtica do problema. E´ tambe´m necessa´rio o domı´nio dos softwares alge´bricos e nume´ricos como
ferramentas extremamente u´teis na resoluc¸a˜o desses problemas.
19
Refereˆncias
[1] Yunus A. C¸engel. Transfereˆncia de Calor e Massa uma Abordagem Pra´tica, 4 Ed. Mc Graw Hill,
Porto Alegre, 2012.
[2] Frank P. Incropera. Fundamentals of Heat Transfer, 7th. John Wiley and Sons, Hoboken - NJ,
1976.
20
Apeˆndices
A. Algoritmo do me´todo anal´ıtico
1 clear all
2 clc
3
4 %PARA A POSICAO EM UM QUARTO DA CHAPA: X=Y=0.25m
5 W=1; %Largura
6 L=1; %Altura
7 n=2000;
8 T2=1;
9 T1=0;
10 x=0.25;
11 y=0.25;
12 i=1;
13 S=0;
14 somatorio=0;
15 while i<=21
16
17 s(i)=((((-1)ˆ(i+1))+1)/i)*sin(i*pi*x/L)*(sinh(i*pi*y/L)/sinh(i*pi*W/L));
18 somatorio=somatorio+s(i);
19 i=i+1;
20 end
21 teta=(2/pi)*somatorio;
22 T=((T2-T1)*teta)+T1;
23 fprintf('A temperatura para x=y=0.25m e: %f\n\n',T);
24
25 %PARA A POSICAO NO MEIO DA CHAPA: X=Y=0.5m
26 x2=0.5;
27 y2=0.5;
28 i=1;
29 S=0;
30 somatorio=0;
31 teta=0;
32 T=0;
33 while i<=21
34
35 s(i)=((((-1)ˆ(i+1))+1)/i)*sin((i*pi*x2)/L)*(sinh((i*pi*y2)/L)/sinh((i*pi*W)/L));
36 somatorio=somatorio+s(i);
37 i=i+1;
38 end
39
40 teta=(2/pi)*somatorio;
41 T=((T2-T1))*teta+T1;
42 fprintf('A temperatura para x=y=0.5m e: %f\n\n',T);
21
B. Algoritmo do me´todo nume´rico
1 clear all
2 clc
3
4 % QUESTAO 1 - JOHNATHAN BATISTA DOS SANTOS E BRUNO FERREIRA COUTO
5 L=1; %comprimento
6 W=1; %altura
7 a=23/10ˆ6; %difusividade termica
8
9 disp('Digite o numero de divisoes para a largura da chapa');
10 nx=input('');%numero de intervalos em x
11 disp('Digite o numero de divisoes para a altura da chapa');
12 ny=input('');%numero de intervalos em y
13
14cont=1;
15 dx=L/nx;
16 dy=W/ny;
17
18 if mod(nx,4)==1 ;
19 fprintf('O numero digitado para x nao e multiplo de 4.');
20 elseif mod(ny,4)==1;
21 fprintf ('O numero digitado para y nao e multiplo de 4')
22 elseif nx~=ny;
23 fprintf('Numero de divisoes em x diferente do numero de divisoes em y.');
24 else
25 for h=3600:900:28800
26 %para as condicoes de contorno impostas esse e um tempo suficiente
27 %para se estabelecer o regime permanente (8 horas)
28 tf=h; %tempo final
29 dt=0.25*(dxˆ2)/a; %intervalo de tempo
30 nt=tf/dt; %numero de passos no tempo
31
32 x=[0:dx:L];
33 y=[0:dy:W];
34 %Temperatura no contorno (parede lateral esquerda e direita e parede inferior)
35 T1=0;
36 %Temperatura no contorno (parede superior)
37 T2=1;
38
39 %condicao inicial
40 T=ones(nx+1,ny+1)*T1;
41
42 %condicoes de contorno nas quinas
43 T(1,1)=(T1+T1)/2;
44 T(nx+1,1)=(T1+T1)/2;
45 T(1,ny+1)=(T1+T2)/2;
46 T(nx+1,ny+1)=(T1+T2)/2;
47 %condicoes de contorno
48 for i=2:nx
49 T(i,1)=T1; %parede inferior
50 T(i,ny+1)=T2; %parede superior
51 end
52
53 for j=2:ny
54 T(1,j)=T1; %lateral esquerda
55 T(nx+1,j)=T1; %lateral direita
56 end
57
22
58 for n=1:nt
59 T0=T;
60 for i=2:nx
61 for j=2:ny
62 T(i,j) = T0(i,j)+dt*a*((T0(i+1,j)-2*T0(i,j)+T0(i-1,j))/(dxˆ2)+(T0(i,j+1)-2*T0(i,j)+T0(i,j-1))/(dyˆ2));
63 end
64 end
65 %grafico dinamico(opcional)
66 %figure(4)
67 %contour(x,y,T',T1:0.05:5)
68 %pause(0.5)
69 end
70
71 t meio(cont)=T((nx/2)+1,((ny/2))+1);%temperatura no meio da chapa x=y=0,5 metros
72 t um quarto(cont)=T((nx/4)+1,(ny/4)+1);%temperatura no ponto x=y=0,25 metros
73 %construcao do vetor "tempo" para construcao da figura 1 e 2
74 tmp(cont)=h/3600;
75 cont=cont+1;
76 end
77 fprintf('A temperatura no meio da chapa (x=y=0,5m) e: %f \n', t meio(cont-1));
78 fprintf('A temperatura no ponto (x=y=0,25m) e: %f \n', t um quarto(cont-1));
79
80 figure(1);
81 plot(tmp,t meio);
82 title('Temperatura no meio da chapa (x=y=0,5 metros)');
83 xlabel('tempo expresso em horas (h)');
84 ylabel('Temperatura em graus celsius (C)');
85
86 figure(2);
87 plot(tmp,t um quarto);
88 title('Temperatura do ponto: x=y=0,25 metros');
89 xlabel('tempo expresso em horas (h)');
90 ylabel('Temperatura em graus celsius (C)');
91
92 figure(3);
93 contour(x,y,T',min(min(T)):0.01:max(max(T)));
94 title('Isotermas');
95 xlabel('Largura da chapa');
96 ylabel('Altura da chapa');
97 end
23
C. Algoritmo do me´todo nume´rico do segundo problema
1 clear all
2 clc
3
4 %QUESTAO 02 - JOHNATHAN E BRUNO COUTO
5
6 L=1; %comprimento
7 W=1; %altura
8 a=23/(10ˆ6); %difusividade termica
9 k=90.7; %condutividade termica do niquel
10 Tinf=5; %temperatura do meio fluido
11 h=100; %coeficiente convectivo
12 nx=21 %Numero de pontos em relacao a x
13 ny=21 %Numero de pontos em relacao a y
14 dx=L/(nx-1);
15 dy=W/(ny-1);
16 cont=1;
17 tmp final=30;
18 erro=100;
19 k=7200;
20
21 while erro>=2
22 tf=k; %tempo final
23 dt=0.25*(dxˆ2)/a; %intervalo de tempo
24 nt=tf/dt; %numero de passos no tempo
25
26 x=[0:dx:L];
27 y=[0:dy:W];
28
29 Ti=300; %Temperatura inicial
30 T1=0; %Temperatura no contorno
31 q1=0; %Fluxo na parede inferior
32 q2=-10; %Fluxo na lateral direita
33 %condicao inicial
34 T=ones(nx+2,ny+2)*Ti;
35
36 for n=1:nt %Laco de tempo
37 T0=T;
38 %condicoes de contorno
39 for i=1:nx+2
40 T(i,1)=T0(i,ny)+2*dy/dx*q1/k; %parede inferior (q0 condicaoo de contorno adiabatica)
41 end
42 for i=2:nx+1
43 T(i,ny+1)=T1
44 end
45 for j=2:ny+1
46 %lateral esquerda (condicao de contorno de conveccao com h =100 e tinf= 5)
47 T(1,j)=T0(3,j)-(2*dx*h/k)*(T0(2,j)-Tinf);
48 %lateral direita (Q fluxo de calor imposto igual a 10 W saindo da chapa)
49 T(nx+2,j)= T0(nx,j)+2*dy/dx*q2/k;
50 end
51
52 %Resolucao da equacao de calor no interior da chapa
53 for i=2:nx+1
54 for j=2:ny
55 T(i,j) = T0(i,j)+dt*a*((T0(i+1,j)-2*T0(i,j)+T0(i-1,j))/(dxˆ2)+(T0(i,j+1)-
56 2*T0(i,j)+T0(i,j-1))/(dyˆ2));
57 end
24
58 end
59
60 for i=1:nx+2
61 for j=1:ny+2
62 TTT(j,i) = T(i,j);
63 end
64 end
65 %grafico dinamico (opcional)
66 %Tmin=min(min(TTT));
67 %Tmax=max(max(TTT));
68 %figure(1)
69 %contour(x,y,TTT(2:nx+1,2:ny+1),Tmin:2:Tmax)
70 %pause(0.2)
71 end
72 tempo(cont)=k/3600;% tempo contado em horas
73 temp(cont)=TTT(((nx-1)/2)+1,((ny-1)/2)+1);%temperatura no meio da placa
74 erro=temp(cont)-tmp final; %erro em relacao a temperatura de 30 graus pedida
75 cont=cont+1;
76 k = k+900;
77 end
78 erro=erro
79
80 tempo para chegar a proximo de 30=tempo(cont-1);
81 temperatura final=temp(cont-1);
82 numero it=cont-1;
83
84 j=tempo(cont-1);
85 hr=((numero it*900)+3600)/3600;
86 d=j-hr;
87
88 minutos=d*60;
89 fprintf('Temperatura final no centro da placa: %f.\n',temperatura final)
90 fprintf('quantidade de horas apara alcancar a temperatura desejada: %f .\n',hr);
91 fprintf('quantidade de minutos: %f.\n',minutos);
92 fprintf('Erro em relacao a temperatura de 30 graus pedida: %f.\n',erro);
93
94 Tmin=min(min(TTT));
95 Tmax=max(max(TTT));
96
97 figure(1)
98 contour(x,y,TTT(2:nx+1,2:ny+1),Tmin:0.01:Tmax)
99 pause(0.2)
100 title('Distribuicao final de temperatura');
101 xlabel('largura da chapa de niquel');
102 ylabel('altura da chapa de niquel');
103
104 figure(2)
105 plot(tempo,temp)
106 title('Evolucao da temperatura no meio da chapa de niquel');
107 xlabel('Tempo (horas)');
108 ylabel('Temperatura (C)');
25
I´ndice Remissivo
Algoritmo, 10, 14, 17
Condic¸o˜es de contorno, 10
Crite´rio de convergeˆncia/estabilidade, 11
Derivada
primeira
segunda, 8
Diferenc¸as finitas, 7
Equac¸a˜o da conduc¸a˜o, 9
Isotermas, 18
Me´todo anal´ıtico, 13
Me´todos nume´ricos, 7
Malha nodal, 9
Modelagem, 7
26
	Introdução
	Métodos numéricos
	O Método das Diferenças Finitas
	Modelagem matemática
	Formulação por diferenças finitas das equações diferenciais
	Condições de contorno
	Formulação numérica
	Considerações gerais
	Problemática
	Primeiro problema
	Abordagem analítica
	Abordagem numérica
	Segundo problema
	Considerações Finais
	Apêndices
	Algoritmo do método analítico
	Algoritmo do método numérico
	Algoritmo do método numérico do segundo problema
	Índice Remissivo