SOLUÇÃO DE UM PROBLEMA UNIDIMENSIONAL
DE CONDUÇÃO DE CALOR
Marcelo M. Galarça
Pós Graduação em Engenharia Mecânica
Universidade Federal do Rio Grande do Sul
Transferência de Calor e Mecânica dos Fluidos Computacional I
[email protected]
Resumo: Neste trabalho obteve-se a distribuição de temperatura de uma placa unidimensional, empregando o Método de
Discretização com Diferenças Finitas Centrais, no espaço; e de forma explícita no tempo, resolvendo o transiente real
com quatro valores de ∆Fo, para uma malha com cinco nós igualmente espaçados. O problema de condução de calor foi
resolvido para uma placa de aço carbono e uma placa de tijolo comum.
Palavras-chave. Método de Diferenças Finitas, Condução de Calor, Problema unidimensional.
1. Introdução
Este trabalho faz a análise numérica de um problema de condução de calor de uma barra unidimensional com
difusividade constante, o que é, um problema elíptico. Utilizou-se o método de diferenças finitas no espaço, sendo, neste
caso, adequado usar o equacionamento de Diferenças Centrais (CDS). Para o tempo utiliza-se a forma explícita, isto é, a
partir de um valor de tempo “velho” obtemos o “novo”.
A situação do problema é representada da seguinte forma:
No instante t = 0 é ligada, no seu interior, uma fonte constante que difunde calor simetricamente. Como nas fronteiras
x = 0 e x = L a barra não está isolada há perda de calor e nestes pontos a temperatura T é zero. Assim a temperatura no
interior da barra é uma função da posição e do tempo.
O problema é estacionário e simétrico, temos um transiente real.
Para L=0,1m
As condições de contorno são:
- T(x,t) = 0
- T(0,t) = 0
- T(L,t) = 0
Para S (fonte), dois valores foram propostos:
- S = 1K/s
- S = 3K/s
A distribuição de temperatura foi calculada para dois casos:
- Placa de Tijolo Comum (α = 0,28e-06 m2/s)
- Placa de Aço Carbono (α = 17,7e-06 m2/s)
Malha:
- Cinco (5) nós igualmente distribuidos
Os valores de ∆Fo empregados foram de 0,25; 0,5; 1 e 2.
A equação governante é a da energia,
∂ ( ρ T ) ∂ ( ρ uT )
∂ k∂ T
+
=
(
) + ST
∂t
∂x
∂ x c p ∂x
S =
ST
ρC P
(1 e 2 )
sendo ρ, Cp e k constantes, a equação pode ser expressada desta forma
∂ 2T
∂T
=α
+S
∂ x2
∂t
(3)
2. Análise do Problema e Metodologia de Solução
O Método de Diferenças Finitas objetiva em calcular a derivada de uma equação, considerando os intervalos, não do
ponto de vista infinitesimal, como no cálculo analítico, mas de forma discretizada em intervalos finitos de dimensão ∆x,
efetuando os cálculos por pontos [2].
Toma-se como referência um determinado ponto, P, e partindo-se deste é feita a expansão da série de Taylor, a qual é
infinita, e por isso é truncada em um determinado termo, de acordo com a ordem de grandeza do erro que se pode assumir
[3]. Para avançar até a convergência de cada iteração, o ponto P “precisa” de mais informações, e assim as busca de seus
“pontos vizinhos”. Essa obtenção das informações pode ser feita utilizando-se um dos três equacionamentos aplicados ao
método: Diferenças Finita Avançada, Atrasada ou Central, como já mencionado, para o problema elíptico em questão, a
solução seguiu a discretização pelo equacionamento central (CDS), no espaço; e de forma explícita no tempo [2].
Diferença Finita Central (CDS): este equacionamento mantém o ponto P no centro, de maneira a relacioná-lo com os
dois outros pontos “a leste” e “a oeste”.(fig.1). A equação se apresenta na forma
φE − φW
 ∂φ 
+ O ( ∆x 2 )
  =
2 ∆x
 ∂x  P
(4)
Discretizando a equação diferencial substituimos o domínio contínuo por um domínio discreto, assim o valor de φ é
obtido em pontos discretos. Para interligá-los obtém-se equações algébricas para cada ponto, e após resolvemos o sistema
de equações resultantes. Estas equações algébricas são derivadas da equação diferencial governante φ e deste modo
expressam a mesma informação física da equação diferencial [3].
Tendo a condição inicial é conveniente a utilização de forma explícita no tempo, θ= 0, onde todas as temperaturas
vizinhas a P são avaliadas no instante anterior e, portanto, já são conhecidas. A formulação explícita dá origem a um
conjunto de equações algébricas que podem ser resolvidas uma a uma, obtendo-se a temperatura em cada ponto do espaço
para o novo instante de tempo [4].
Discretizando a equação, temos :
(T n +θ + TWn+θ − 2TPn+θ )
T n+θ − TWn +θ
TPn+1 − TPn
+u E
=α E
∆t
2∆x
∆x 2
Para a formulação explícita, temos θ = 0, ficando:
T n − TWn
(T n + TWn − 2TPn )
TPn+1 − TPn
+u E
=α E
∆t
2∆x
∆x 2
Onde :
TP = Temperatura no ponto P no instante anterior
TE = Temperatura no ponto vizinho a leste no instante anterior
TW = Temperatura no ponto vizinho a oeste no instante anterior
u = velocidade na direção x
α = difusividade térmica
n = valor “velho”
∆x = espaçamento entre pontos
Calculando para Tpn+1, temos
(5)
(6)
T pn+1 = (1 −
2α ∆t n ∆tα u∆t
u∆t α ∆t n
)T p + ( 2 −
)TEn + (
+
)TW
2
2∆x ∆ x 2
∆x
∆x
2∆ x
(7)
Para assegurar a positividade dos coeficientes, temos que:
2∆ Fo ≤ 1
Logo:
α∆t
∆x
2
≤
1
2
(8)
Foi analisado o transiente real com quatro valores de ∆Fo = 0.25, 0.5, 1.0 e 2.0 sempre com a mesma malha espacial de
cindo nós, igualmente espaçada. Com base no critério de convergência em (8), é esperado que os dois últimos valores
testados, de ∆Fo = 1 e 2, não venham a dar bons resultados quanto à convergência.
3. Ferramentas
A solução do problema foi auxiliada pela utilização do software Compaq Visual Fortran Professional Edition 6.0.1,
no qual um algoritmo foi elaborado com a linguagem de programação compatível, o mesmo encontra-se em anexo.
O software Microsoft Excel 2002, foi utilizado para graficar os resultados.
4. Resultados e Discussão
No desenvolvimento do trabalho observou-se que para os diferentes valores de ∆Fo, como já era esperado, somente
houve convergência para 0,25 e 0,5. Isto concorda com o critério de positividade dos coeficientes, onde o intervalo de
tempo em formulações explícitas é limitado pela estabilidade.
Um outra observação é com relação à quantidade de iterações para a convergência, a qual é menor para ∆Fo = 0.5,
enquanto, que para ∆Fo = 0.25 é maior o número de iterações, isso ocorre devido à redução do ∆t, sendo necessário um
número de passos maior para se atingir o mesmo valor de temperatura convergido no outro ∆Fo.
Os gráficos que seguem representam as distribuições de temperatura da barra com ∆Fo = 0,25; 0,5; 1,0 e 2,0 para as
fontes S=1K/s e 3K/s.
Pode ser observado claramente a divergência para os valores de ∆Fo iguais a 1 e 2.
Cam po de Tem peratura (S=1K/s)
Aço Carbono, Fo=0,25
8,00E+01
103 Iterações
Temperatura (K)
7,00E+01
10 Iterações
6,00E+01
5 Iterações
5,00E+01
3 Iterações
4,00E+01
2 Iterações
3,00E+01
50 Iterações
2,00E+01
25 Iterações
1,00E+01
15 Iterações
0,00E+00
0,000
0,020
0,040
0,060
0,080
Com prim ento da Malha (m )
0,100
Cam po de Tem peratura (S=1K/s)
Aço Carbono, Fo=0,5
8,00E+01
Temperatura (K)
7,00E+01
6,00E+01
51 Iterações
5,00E+01
10 Iterações
4,00E+01
5 Iterações
3,00E+01
3 Iterações
2 Iterações
2,00E+01
1,00E+01
0,00E+00
0,000
0,020
0,040
0,060
0,080
0,100
Com prim ento da Malha (m )
Cam po de Tem peratura (S=1K/s)
Aço Carbono, Fo=1
8,00E+01
Temperatura (K)
7,00E+01
6,00E+01
5,00E+01
3 Iterações
4,00E+01
2 Iterações
3,00E+01
2,00E+01
1,00E+01
0,00E+00
0,000
0,020
0,040
0,060
0,080
Com prim ento da Malha (m )
0,100
Cam po de Tem peratura (S=1K/s)
Aço Carbono, Fo=2
3,00E+03
Temperatura (K)
2,50E+03
2,00E+03
1,50E+03
1,00E+03
5 Iterações
5,00E+02
3 iterações
0,00E+00
-5,00E+02
-1,00E+03
-1,50E+03
-2,00E+03
0,000
0,020
0,040
0,060
0,080
0,100
Comprimento da Malha (m)
Como já mencionado, após a observação dos gráficos, os valores de ∆Fo = 1 e 2 apresentaram divergência. Para ∆Fo =
1, não foi possível efetuar mais do que três (03) iterações, já para ∆Fo = 2 foi possível obter até 250 iterações, as quais não
foram mostradas em gráfico devido aos valores de temperatura muito altos e, também por apresentarem comportamento
semelhante ao obtido em até cinco (05) iterações. Ambas explicações se devem ao critério de positividades entre os
coeficientes, o qual foi desrespeitado com estes valores.
Os resultados obtidos assumindo que a placa é feita de tijolo comum, graficamente apresentaram o mesmo
comportamento que os gráficos do aço. Desta forma, para o valor de ∆Fo = 0,25 foram necessárias mais iterações para
atingir o regime permanente do que as exigidas para o valor de ∆Fo = 0,5. Os outros dois valores testados apresentaram
divergência e seu comportamento gráfico foi o mesmo que o apresentado pelo aço. As altas temperaturas atingidas para a
placa de tijolo são atribuídas ao fato de que a difusividade térmica (α = 0,28 E-06 m2/s) é muito menor do que a do aço (α =
17,7 E-06 m2/s), o que representa uma retenção de energia muito maior do tijolo.
Devido ao fato dos resultados gráficos para o tijolo terem sido iguais aos do aço, foi apresentado neste trabalho apenas
os gráficos que representam a convergência total, ou seja, já em regime permanente.
Temperatura (K)
Cam po de Tem peratura (S=1)
Tijolo Com un
5,00E+03
4,50E+03
4,00E+03
3,50E+03
3,00E+03
2,50E+03
2,00E+03
1,50E+03
1,00E+03
5,00E+02
0,00E+00
0,000
Fo=0,25 - 129 iterações
Fo=0,5 - 63 iterações
Fo=1- 3 iterações
0,020
0,040
0,060
0,080
Com prim ento da Malha (m )
0,100
Cam po de Tem peratura (S=3K/s)
Tijolo Com um
1,60E+04
Temperatura (K)
1,40E+04
Fo=0,25 - 136 iterações
1,20E+04
1,00E+04
Fo=0,5 - 65 iterações
8,00E+03
6,00E+03
Fo=1 - 3 iterações
4,00E+03
2,00E+03
0,00E+00
0,000
0,020
0,040
0,060
0,080
0,100
Com prim ento da Malha (m )
5. Conclusões
Conforme apresentado nos gráficos, fica claro o fato de que o critério da positividade dos coeficientes permite-nos
definir um valor limite de ∆Fo, ou ainda um valor de ∆t, pois o mesmo está contido no número Fourier, para que o
problema apresente solução. O não cumprimento do critério provoca a divergência do programa, e conseqüentemente dos
resultados.
Outra observação evidente é de que a solução do problema é a mesma para qualquer valor de ∆Fo, ou seja, os valores de
temperatura são iguais para pontos iguais em tempos iguais e atingem o regime permanente no mesmo instante de tempo, o
que pôde ser observado através do número de iterações. Desta forma, a solução para problemas de condução de calor
apresentada no trabalho independe do valor de ∆Fo.
6. Referências
[1] Notas de Aula.
[2] INCROPERA, Frank P.; WITT, David P. Fundamentos de Transferência de Calor e de Massa. Ed. LTC. 3a Ed.1992.
[3] POZRIKIDIS, C.. Fluid Dynamics - Theory, Computation, and Numerical Simulation. Kluwer Academic Publishers. 2001.
[4] PATANKAR, Suhas V.. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Company. 1980.
7. Anexo: Algoritmo do Programa