MÉTODO DAS DIFERENÇAS FINITAS EM MATLAB

Prévia do material em texto

TER02 - Transferência de Calor
Prof. David Roza José
Rua Vigário Frei João, nº 550 - Centro
CEP: 89609-000 Luzerna - SC
Fone (49) 3523-4319
david.jose@luzerna.ifc.edu.br
SIMULAÇÃO DE TRANSFERÊNCIA DE CALOR EM REGIME TRANSIENTE ATRAVÉS DO
MÉTODO DAS DIFERENÇAS FINITAS
TRABALHO REGIME TRANSIENTE
Hantony Matheus Zimmermann
Resumo
Este documento é o resultado de um estudo sobre o “Método das Diferenças Finitas” em regime
transiente. Inicialmente, será calculada a temperatura em cada nó no decorrer do tempo. Em seguida,
será gerado gráficos do tipo superfície para mostrar o comportamento dessas temperaturas em cada
um dos nós. O código em Matlab está disponível no final do documento.
1
1 Introdução
Este trabalho destina-se a um estudo específico sobre o Método das Diferenças Finitas em
Regime Transiente. Sua formatação clara e concisa permite que o mesmo sirva de fonte
para pesquisas de outros alunos ou professores que buscam uma maneira de observar o
comportamento da temperatura em um sólido, utilizando ferramentas matemáticas no de-
selvolver do método, bem como o Matlab. Nele será descrito o problema de transferência
de calor, o código em Matlab da resolução, o resultado final e o gráfico do tipo superfície
para ilustrar a transferência de calor através do sólido.
Incropera et al. (2007, p.188) afirmam que, soluções analíticas para problemas transientes
estão restritas a geometrias e condições de contorno simples, [...]. Contudo, em muitos
casos, a geometria e/ou as condições de contorno descartam totalmente a possibilidade do
uso de técnicas analíticas, tornando necessária a utilização de método de diferenças finitas.
Incropera et al. (2007, p.192) ainda dizem que, uma redução no tempo computacional pode
frequentemente ser obtida com o emprego de um esquema de diferenças finitas implícito,
no lugar do método explícito. Em relação ao método explícito, a formulação implícita tem
a importante vantagem de ser incondicionalmente estável. Isto é, a solução permanece
estável para todos os intervalos de espaço e de tempo, não havendo restrições em ∆x e ∆t.
Com base nisso, o objetivo do trabalho é aplicar o método das diferenças finitas em re-
gime transiente em um sólido, cujo condição inicial é conhecida, bem como as dimensões
do mesmo, e descobrir quais são as temperaturas em cada um dos nós do seu perfil, ge-
rando assim gráficos do tipo superfície para ilustrar a transferência de calor no material no
decorrer do tempo.
2 Fundamentação teórica
“Transferência de calor é energia térmica em transito devido a uma diferença de tempera-
turas no espaço.”(INCROPERA, p.2)
Segundo Incropera et al. (2007, p.135), sobre o estudo do método das diferenças finitas
em regime transiente, dizem que, ao contrário de uma solução analítica, que permite a
determinação da temperatura em qualquer ponto de interesse de um meio, uma solução
numérica permite a determinação da temperatura em somente pontos discretos. A precisão
numérica dos cálculos depende fortemente do número de pontos nodais utilizados. Se este
número for grande, soluções precisas são obtidas.
2.1 Convecção
Como relatam Incropera et al (2011, p.6), o modo de transferência de calor por convecção
é descrito como a transferência de energia ocorrendo no interior de um fluido devido aos
efeitos combinados da condução e do escoamento global ou macroscópico do fluido. Inde-
pendentemente da natireza específica do processo de transferência de calor por convecção,
a equação apropriada para a taxa de transferência de calor possui a forma
q” = h(Ts−T∞)
2
onde q”, o fluxo de calor por convecção [W/m²], é proporcional à diferença entre as tem-
peraturas da superfície e do fluido. Esta expressão é conhecida como a lei de resfriamento
de Newton.
2.2 Condução
Segundo Incropera et al. (2011, p.3), a condução pode ser vista como a transferência
de energia das partículas mais energéticas para as partículas menos energéticas de uma
substância devido às interações entre partículas. É possível quantificar processos de trans-
ferência de calor em termos de equações de taxa apropriadas. Essas equações podem ser
usadas para calcular a quantidade de enerfia sendo transferida por unidade de tempo. Para
a condução térmica, a equação da taxa é conhecido como lei de Fourier.
q” = kT2−T1L
2.3 Resistência térmica
Incropera et al (2011, p. 115) afirmam que, circuitos térmicos equivalentes também podem
ser usados em sistemas mais complexos, como, por exemplo, paredes compostas. Tais
paredes podem possuir uma quantidade qualquer de resistências térmicas em série e em
paralelo, devido a presença de camadas de diferentes materiais. A taxa de transferência de
calor unidimensional para esse sistema pode ser representada por
qx =
T∞,1−T∞,n
ΣRt
onde T∞,1 e T∞,n é a diferença de temperaturas global e o somatório inclui todas as resis-
tências.
3 Apresentação do problema
Considere um duto retangular conforme a figura. Neste duto escoarão gases quentes de um
sistema de exaustão. As paredes são construídas com tijolo refratário de condutividade k=
0.85W/mK e α = 5.5x10−7m2/s. Num instante inicial as paredes estão com a temperatura
de 25°C. Subitamente, a superfície interna está exposta a um gás quente com T∞,i = 350°C
e hi = 100W/m²K, enquanto a superfície externa está exposta a um T∞,e = 25°C com
he = 5W/m²K. Utilizando o método implícito das diferenças finitas com um incremento
de tempo de 1h, encontre a distribuição de temperaturas na parede 1, 3, 5, 10, 50 e 100h
após o início do escoamento dos gases quentes. Utilize ∆x= ∆y= 50mm. Aproveite-se da
simetria do problema para reduzir a quantidade de nós e elementos. Execute as mesmas
análises anteriores, mas agora com ∆x= ∆y= 25mm.
3
Figura 1: Geometria do problema.
3.1 Simetria
Como o perfil do duto é simétrico, se calcularmos as temperaturas para um dos lados, con-
sequentemente a mesma distribuição será para o outro lado, por esse motivo, foi escolhido
dividir o duto ao meio e calcular as temperaturas somente para uma parte. Isso otimiza o
tempo de programação, não deixa o programa tão pesado e também auxilía nos cálculos.
Fonte: Dados primários.
3.2 Metodologia
Para problemas de condução em um corpo desta complexidade, torna-se mais fácil a uti-
lização de softweres para a solução das equações. O softwere utilizado para a confecção,
solução e análise do problema foi o Matlab.
O método escolhido foi o método implícito de diferenças finitas em regime transiente, o
qual não tem restrições quanto ao uso de ∆x e ∆t. Este método oferece as equações para
cada caso de análise de um sólido. Os casos e equações utilizados neste trabalho foram
4
Fonte: Incropera et al, 2011, p.334.
Para cada nó do perfil do duto existe uma equação específica. Neste caso, foram utilizadas
todas as equações do método de regime transiente para a solução do problema. Com isso,
foi gerado um código em Matlab que mostra a distribuição das temperaturas nodais e que
também ilustra este comportamento através de gráficos do tipo superfície.
4 Resultados e discussões
De acordo com a apresentação do problema, o primeiro passo para a solução do mesmo é
saber como chegar nela. Com o passo de tempo de uma hora, podemos verificar através
dos gráficos a mudança na temperatura dos nós.
O valor numérico das temperaturas em cada um dos nós, pode ser visualizado na matriz
de temperaturas T2, inclusa no programa.
Inicialmente foi utilizado o tamanho da malha de 50 mm e posteriormente foi diminuido
o tamanho para 25mm. A temperatura inicial antes dos gases escoarem é de 25°C. Como
mostra a figura abaixo, o escoamento ainda não foi liberado.
5
0
5
10
15
0
2
4
6
8
−1
−0.5
0
0.5
1
Fonte: Dados primários.
Após o início do escoamento, o calor começa a ser transferido para o duto. A figura
abaixo mostra a distribuição de temperatura uma hora depois do iníciodo escoamento dos
gases. A seção onde o gráfico está avermelhado diz respeito a parte interna do duto onde a
temperatura está a T∞,i =350°C. A temperatura vai diminuindo conforme a geometria do
duto, pois a temperatura externa aos gases está a T∞,e = 25°C.
0
5
10
15
0
2
4
6
8
200
300
400
500
600
700
Fonte: Dados primários.
6
Na figura abaixo, o tempo de escoamento aumentou para três horas. Percebe-se que foi
dissipado mais calor para o duto.
0
5
10
15
0
2
4
6
8
200
300
400
500
600
700
Fonte: Dados primários.
Sucessivamente, a medida que o tempo passa, maior fica a temperatura nos nós que estão
longe dos gases quentes. Conforme mostra a figura abaixo, os gases foram liberados há
cinco horas.
7
0
5
10
15
0
2
4
6
8
200
300
400
500
600
700
Fonte: Dados primários.
10 horas após o início do escoamento.
0
5
10
15
0
2
4
6
8
300
400
500
600
700
Fonte: Dados primários.
50 horas após o início do escoamento.
8
0
5
10
15
0
2
4
6
8
350
400
450
500
550
600
650
Fonte: Dados primários.
E finalmente 100 horas após o início do escoamento.
0
5
10
15
0
2
4
6
8
350
400
450
500
550
600
650
Fonte: Dados primários.
Podemos notar que, a partir de 50 horas de escoamento, a temperatura começa a estabilizar
vagarosamente. Isto é, como escoam outros gases fora do duto a temperatura de 25°C e
9
um coeficiente convectivo de 5W/m²K, o calor chega até a superfície externa e começa a
ser dissipado por convecção com os gases externos.
Diminuindo o tamanho da malha para 25 mm, podemos obter um “refinamento” dos va-
lores numéricos da temperatura, ou seja, o número de nós aumenta e como o objetivo do
método das diferenças finitas é saber a temperatura em pontos discretos, com um número
maior de nós, podemos saber a temperatura em outros pontos que não conhecíamos com a
malha em 50 mm.
Abaixo, seguem os resultados obtidos com tal malha. O primeiro gráfico mostra o escoa-
mento após uma hora.
0
5
10
15
20
25
0
5
10
15
200
300
400
500
600
700
Fonte: Dados primários.
Três horas depois.
10
0
5
10
15
20
25
0
5
10
15
200
300
400
500
600
700
Fonte: Dados primários.
Cinco horas depois.
0
5
10
15
20
25
0
5
10
15
200
300
400
500
600
700
Fonte: Dados primários.
Dez horas depois.
11
0
5
10
15
20
25
0
5
10
15
200
300
400
500
600
700
Fonte: Dados primários.
Cinquenta horas depois do início do escoamento.
0
5
10
15
20
25
0
5
10
15
300
400
500
600
700
Fonte: Dados primários.
E finalmente, cem horas depois do início do escoamento.
12
0
5
10
15
20
25
0
5
10
15
300
400
500
600
700
Fonte: Dados primários.
4.1 Erros
Devido o valor do coeficiente convectivo externo ser baixo, automaticamente o próprio
programa insere um erro de truncamento, ou seja, é um erro de expressões matemáticas
em um número finito de passos. Por esse motivo, o valor da temperatura calculada no
programa é diferente do calculado pelo regime estacionário. A título de comparação, foi
feito o cálculo com um coeficiente convectivo externo de 500W/m²K, as tabelas abaixo
com as malhas de 50 e 25 mm, mostram que com esse valor, comparando com o regime
estacionário, a diferença é mínima entre as temperaturas da superfície interna e externa.
Malha = 50 mm
Regime h [W/m²K] Ts,i [K] Ts,e [K] h [W/m²K] Ts,i [K] Ts,e [K]
Estacionário
5
614,5906 466,1887
500
605,7559 301,4488
Transiente 613,1646 454,3093 605,0570 301,3182
Com a malha de 25 mm, a diferença das temperaturas é maior com o coeficiente convectivo
de 5W/m²K.
Malha = 25 mm
Método h [W/m²K] Ts,i [K] Ts,e [K] h [W/m²K] Ts,i [K] Ts,e [K]
Estacionário
5
614,5906 466,1887
500
605,7559 301,4488
Transiente 609,5414 403,6925 606,36989 300,66988
13
Analisando as tabelas acima, percebe-se que com um coeficiente convectivo de 500W/m²K
as diferenças de temperatura é menor do que um grau. Já com o coeficiente de 5W/m²K, a
diferença aumenta devido ao erro incluso no próprio programa.
5 Conclusão
Este estudo permitiu destacar a maneira pela qual calculamos a taxa de transferência de ca-
lor através do método do regime estacionário, visto em aula, o qual auxiliou para comparar
os resultados com o método das diferenças finitas, saber se o programa está convergente e
como principal objetivo, calcular a transferência de calor através do método das diferenças
finitas aplicadas ao Matlab. Com o auxílio de tal método, conseguimos plotar um gráfico
do tipo superfície para visualizar melhor a transferência de energia para o sólido, durante
o passar do tempo.
Ressaltamos também que, este trabalho está aberto para futuras alterações, ou seja, estudar
outras geometrias através do método das diferenças finitas em regime transiente, utilizando
o método implícito para o cálculo das temperaturas nodais.
6 Anexos
Para inserir código em MATLAB, ir em Inserir -> Arquivo -> Documento Filho; escolher
“Listagem de Programa” e apontar para o arquivo.m . Note que o arquivo .m não pode ter
caracteres especiais tais como “²” ou “³”.
1 %%%%%%%%%%%%%%%%%%%%
2 % Exemplo de Utilização de Diferenças Finitas
3 % Regime Transiente
4 % Professor David Roza José
5 % Acadêmico Hantony Matheus Zimmermann
6 % Transferência de Calor - Engenharia Mecânica
7 % Instituto Federal Catarinense - Campus Luzerna
8 % 02/06/16
9 %%%%%%%%%%%%%%%%%%%%
10 clear all
11 close all
12 clc
13 format long
14
15 %%%%%%% %%%% %%%%%%%%
16 % DADOS DO PROGRAMA %
17 %%%%%%% %%%% %%%%%%%%
18
19 %Entrada de dados
20 MaxHoras = 100;
21 DeltaX = 50e-3; %DeltaX em metros.
22 DeltaY = DeltaX; %DeltaY em metros.
23 Mult = int8 (50e-3/ DeltaX);
24 Lx = 600e-3; %Comprimento em X.
25 Ly = 300e-3; %Comprimento em Y.
26 NNx = int8((Lx/DeltaX) + 1); %Numero de nós em X.
14
27 NNy = int8((Ly/DeltaY) + 1); %Numero de nós em Y.
28 ResMax = 1e-6; %Valor do resíduo
29 Res = 1; %Valor inicial do resíduo
30 iter = 0; %Iteração inicial
31 k = 0.85; %W/mK
32 alfa = 5.5e-7; %m^2/s
33 Tiinf = 623; %K
34 Teinf = 298; %K
35 hi = 100; %W/m^2K
36 he = 500; %W/m^2K
37 DeltaT = 3600; %1 hora %segundos
38 TA = 298; %K
39 Fo = (alfa*DeltaT)/(( DeltaX)^2);
40 Bii = (hi*DeltaX)/k;
41 Bie = (he*DeltaX)/k;
42 L = 0.15; %Metros
43
44 %Alocação de matrizes
45 T1=zeros(NNy ,NNx);
46 T2=zeros(NNy ,NNx);
47 T3=zeros(NNy ,NNx);
48 Temp=zeros(NNy ,NNx);
49 Dif=zeros(NNy ,NNx);
50 VecRes1=zeros (100000 ,1);
51 VecRes2=zeros (100000 ,1);
52
53 %%%%%%%%%%%%%%%%%%%%%%%%
54 % CONDIÇÃO INICIAL %
55 %%%%%%%%%%%%%%%%%%%%%%%%
56
57
58 for i=1:NNy;
59 for j=1:NNx;
60 if((i>=1 && i<=( Mult *4)) && (j>=1 && j<=NNx))
61 T1(i,j) = TA;
62 T2(i,j) = TA;
63 Temp(i,j) = TA;
64 elseif (((i>(Mult *4) && i<=NNy) && (j>=1 && j<=( Mult *4)))
|| ((i>(Mult *4) && i<=NNy) && (j>=( Mult *10) && j<=NNx)))
65 T1(i,j) = TA;
66 T2(i,j) = TA;
67 Temp(i,j) = TA;
68 else
69 T1(i,j) = 623;
70 T2(i,j) = 623;
71 Temp(i,j) = 623;
72 end
73 end
74 end
75
76 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15
77 % CÁLCULO ANALÍTICO EM REGIME ESTACIONÁRIO DAS TEMPERATURAS
SUPERFICIAIS %
78 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
79
80
81 fprintf(’\nPelo método estacionário , a taxa de calor e as
temperaturas são:\n’)
82
83 % TAXA DE CALOR %
84
85 q = (Tiinf - Teinf)/((1/ hi)+(L/k)+(1/he));
86 fprintf(’\nTaxa de calor = é %.4f\n’,q);
87
88 % TEMPERATURA SUPERFÍCIE INTERNA %
89
90 Tsi = Tiinf -(q/hi);
91 fprintf(’\nTemperatura da superfícieinterna = %.4f\n’,Tsi);
92
93 % TEMPERATURA SUPERFÍCIE EXTERNA %
94
95 Tse = Teinf+(q/he);
96 fprintf(’\nTemperatura da superfície externa é %.4f’,Tse);
97
98
99 %%%%%%%%%%%%%%%%%%%%%%
100 % PROCESSO ITERATIVO %
101 %%%%%%%%%%%%%%%%%%%%%%
102
103 Hora = 0;
104
105 while (Res >= ResMax && Hora < MaxHoras)
106
107 iter=iter +1;
108 Hora = Hora + 1;
109
110 i=1;
111 j=1;
112 % Caso 4 - Canto Superior Esquerdo
113 T3(i,j) = (T1(i,j) + 4*Bie*Fo*Teinf + 2*Fo*(T2(i,j+1) + T2(i
+1,j)))/(1+(4* Fo*(1+ Bie)));
114
115
116 i=1;
117 j=NNx;
118 % Caso 4 - Canto Superior Direito
119 T3(i,j) = (T1(i,j) + 4*Bie*Fo*Teinf + 2*Fo*(T2(i,j-1) + T2(i
+1,j)))/(1+(4* Fo*(1+ Bie)));
120
121
122 i=(Mult *4);
123 j=(Mult *4);
16
124 % Caso 2 - vertice Interno Esquerdo
125 T3(i,j) = (T1(i,j) + ((4/3)*Bii*Fo*Tiinf) + ((2/3)*Fo*(T2(i,j
+1) + 2*T2(i,j-1) + 2*T2(i-1,j) + T2(i+1,j))))/(1 + (4*Fo
*(1+(1/3)*Bii)));
126
127
128 i=(Mult *4);
129 j=(Mult *10);
130 % Caso 2 - vertice Interno Direito
131 T3(i,j) = (T1(i,j) + ((4/3)*Bii*Fo*Tiinf) + ((2/3)*Fo*(T2(i,j
-1) + 2*T2(i,j+1) + 2*T2(i-1,j) + T2(i+1,j))))/(1 + (4*Fo
*(1+(1/3)*Bii)));
132
133
134 i=NNy;
135 j=1;
136 % Caso 3 - Canto Inferior Esquerdo Externo
137 T3(i,j) = (T1(i,j) + 2*Bie*Fo*Teinf + Fo*(2*T2(i,j+1) + 2*T2(i
-1,j)))/(1+(2* Fo*(2+ Bie)));
138
139
140 i=NNy;
141 j=NNx;
142 % Caso 3 - Canto Inferior Direito Externo
143 T3(i,j) = (T1(i,j) + 2*Bie*Fo*Teinf + Fo*(2*T2(i,j-1) + 2*T2(i
-1,j)))/(1+(2* Fo*(2+ Bie)));
144
145
146 i=NNy;
147 j=(Mult *4);
148 % Caso 3 - Canto Inferior Esquerdo Interno
149 T3(i,j) = (T1(i,j) + 2*Bii*Fo*Tiinf + Fo*(2*T2(i,j-1) + 2*T2(i
-1,j)))/(1+(2* Fo*(2+ Bii)));
150
151
152 i=NNy;
153 j=(Mult *10);
154 % Caso 3 - Canto Inferior Direito Interno
155 T3(i,j) = (T1(i,j) + 2*Bii*Fo*Tiinf + Fo*(2*T2(i,j+1) + 2*T2(i
-1,j)))/(1+(2* Fo*(2+ Bii)));
156
157
158 j=(Mult *10);
159 for i=(( Mult *4)+1):(NNy -1);
160 % Caso 3 - Borda Direita Interna
161 T3(i,j) = (T1(i,j) + 2*Bii*Fo*Tiinf + Fo*(2*T2(i,j+1) + T2
(i-1,j) + T2(i+1,j)))/(1+(2* Fo*(2+ Bii)));
162 end
163
164 j=NNx;
165 for i=2:(NNy -1);
166 % Caso 3 - Borda Direita Externa
167 T3(i,j) = (T1(i,j) + 2*Bie*Fo*Teinf + Fo*(2*T2(i,j-1) + T2
17
(i-1,j) + T2(i+1,j)))/(1+(2* Fo*(2+ Bie)));
168 end
169
170 j=(Mult *4);
171 for i=(( Mult *4)+1):(NNy -1);
172 % Caso 3 - Borda Esquerda Interna
173 T3(i,j) = (T1(i,j) + 2*Bii*Fo*Tiinf + Fo*(2*T2(i,j-1) + T2
(i-1,j) + T2(i+1,j)))/(1+(2* Fo*(2+ Bii)));
174 end
175
176 j=1;
177 for i=2:(NNy -1);
178 % Caso 3 - Borda Esquerda Externa
179 T3(i,j) = (T1(i,j) + 2*Bie*Fo*Teinf + Fo*(2*T2(i,j+1) + T2
(i-1,j) + T2(i+1,j)))/(1+(2* Fo*(2+ Bie)));
180 end
181
182 i=1;
183 for j=2:(NNx -1);
184 % Caso 3 - Borda Superior Externa
185 T3(i,j) = (T1(i,j) + 2*Bie*Fo*Teinf + Fo*(2*T2(i+1,j) + T2
(i,j+1) + T2(i,j-1)))/(1+(2* Fo*(2+ Bie)));
186 end
187
188 i=(Mult *4);
189 for j=(( Mult *4)+1):(( Mult *10) -1);
190 % Caso 3 - Borda Superior Interna
191 T3(i,j) = (T1(i,j) + 2*Bii*Fo*Tiinf + Fo*(2*T2(i-1,j) + T2
(i,j+1) + T2(i,j-1)))/(1+(2* Fo*(2+ Bii)));
192 end
193
194 i=NNy;
195 for j=2:(( Mult *4) -1);
196 % Caso 1 - Borda Inferior Externa (Local do Corte)
197 T3(i,j) = (T1(i,j) + Fo*(T2(i,j+1) + T2(i,j-1) + 2*T2(i-1,
j)))/(1+(4* Fo));
198 end
199
200 i=NNy;
201 for j=(( Mult *10) +1):(NNx -1);
202 % Caso 1 - Borda Inferior Externa (Local do Corte)
203 T3(i,j) = (T1(i,j) + Fo*(T2(i,j+1) + T2(i,j-1) + 2*T2(i-1,
j)))/(1+(4* Fo));
204 end
205
206 for i=(( Mult *4)+1):NNy;
207 for j=(( Mult *4)+1):(( Mult *10) -1);
208 % Parte Sem nós
209 T3(i,j) = T1(i,j);
210 end
211 end
212
213 for i=2:(( Mult *4) -1);
18
214 for j=2:(NNx -1);
215 % Caso 1 - Interior do solido
216 T3(i,j) = (T1(i,j) + Fo*(T2(i,j+1) + T2(i,j-1) + T2(i
-1,j) + T2(i+1,j)))/(1+(4* Fo));
217 end
218 end
219
220 for i=(Mult *4):(NNy -1);
221 for j=2:(( Mult *4) -1);
222 % Caso 1 - Interior do solido
223 T3(i,j) = (T1(i,j) + Fo*(T2(i,j+1) + T2(i,j-1) + T2(i
-1,j) + T2(i+1,j)))/(1+(4* Fo));
224 end
225 end
226
227 for i=(Mult *4):(NNy -1);
228 for j=(( Mult *10) +1):(NNx -1);
229 % Caso 1 - Interior do solido
230 T3(i,j) = (T1(i,j) + Fo*(T2(i,j+1) + T2(i,j-1) + T2(i
-1,j) + T2(i+1,j)))/(1+(4* Fo));
231 end
232 end
233
234 Dif=abs(T3 -T1);
235 Res=max(Dif(:));
236 T1 = T3;
237 T2 = T3;
238 VecRes1(iter ,1)=Res;
239 VecRes2(iter ,1)=iter;
240
241 end
242
243 VecRes1 = VecRes1 (1: iter);
244 VecRes2 = VecRes2 (1: iter);
245 surf(T3)
19
Referências
[1] Frank Incropera, David Dewitt, Theodore Bergman, and Adrienne Lavine. Funda-
mentals of heat and mass transfer. John Wiley e Sons, Inc., Jefferson City, 7 edition,
2011.
20