Trabalho 3 de Volumes Finitos

Prévia do material em texto

Trabalho 3 de Volumes Finitos 
Professor: José Adilson 
Aluno: Renner Egalon Pereira - PGMEC 
 
Estudo da condução estacionária 2D em uma placa de material 
homogêneo e propriedades constantes e lineares 
 
O problema se baseia no estudo da condução de calor em uma placa com 1 metro de largura e 1 metro 
de altura. Neste problema, foram utilizados 100 volumes de controle para análise da placa. Os 
algoritmos para solução foram realizados através do software MATLAB© versão R2015a. 
A ideia dos volumes de controle está disposta na figura abaixo: 
 
A equação de governo para o problema em análise é dada por: 
𝜕2𝑇
𝜕𝑥2
+
𝜕2𝑇
𝜕𝑦2
= 0 
 
Caso 1: Problema proposto no livro Incropera Introduction Heat Transfer 6th com temperaturas 
fixas prescritas em todas as bordas da placa 
 
 
 
 
 
 
 
 
 
 
Condições de Contorno: 
𝑇(0, 𝑦) = 100°𝐶 𝑇(1, 𝑦) = 100°𝐶 𝑇(𝑥, 0) = 100°𝐶 𝑇(𝑥, 1) = 600°𝐶 
 
Solução Teórica: 
A solução exata foi obtida pela solução analítica utilizando a teoria presente no livro Incropera 
Introduction Heat Transfer 6th (cap 4, p. 234-235) 
 
De acordo com a solução teórica, a temperatura T(x,y) da placa é dada por: 
 
𝑇(𝑥, 𝑦) = 𝑇1 + (𝑇2 − 𝑇1) ∗ 𝜃(𝑥, 𝑦) 
 
𝑇1 = 100°𝐶 
𝑇2 = 600°𝐶 
 
𝜃(𝑥, 𝑦) =
2
𝜋
∑
(−1)𝑛+1 + 1
𝑛
sin (
𝑛𝜋𝑥
𝐿
)
sinh (
𝑛𝜋𝑦
𝐿 )
sinh (
𝑛𝜋𝑊
𝐿 )
∞
𝑛=1
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
A figura acima mostra a previsão teórica de isotermas e linhas de fluxo de calor para condução 
bidimensional na placa (Adaptado: Incropera Introduction Heat Transfer 6th) 
 
Para este problema foi considerado 𝑛 (í𝑚𝑝𝑎𝑟𝑒𝑠) = 1, 3, 5, 7, 9, … 99. 
 
O número final de n foi obtido a partir da estabilização da temperatura para um determinado ponto 
(x,y) em função da variação de n. Neste caso foi escolhido o ponto (0.4,1). 
n T (°C) 
9 548,9716 
19 549,903 
29 549,2256 
39 549,208 
49 549,2227 
59 549,2232 
69 549,2228 
79 549,2228 
89 549,2228 
99 549,2228 
 
 
Para a solução exata, a matriz de temperaturas calculadas (excluindo as temperaturas nos contornos), 
onde 0.05 ≤ 𝑥 ≤ 0.95 e 0.05 ≤ 𝑦 ≤ 0.95, é: 
 
348,6322 493,3941 529,8761 543,8956 549,2228 549,2228 543,8956 529,8761 493,3941 348,6322 
198,3961 337,6998 406,6269 439,2259 452,55 452,55 439,2259 406,6269 337,6998 198,3961 
156,3912 252,3537 316,0142 352,1963 368,3362 368,3362 352,1963 316,0142 252,3537 156,3912 
136,6351 202,9046 252,8089 284,5098 299,6228 299,6228 284,5098 252,8089 202,9046 136,6351 
124,9705 171,3007 208,3173 233,3174 245,7467 245,7467 233,3174 208,3173 171,3007 124,9705 
117,234 149,5909 176,2097 194,7888 204,2533 204,2533 194,7888 176,2097 149,5909 117,234 
111,725 133,8682 152,3575 165,4902 172,2711 172,2711 165,4902 152,3575 133,8682 111,725 
107,5729 121,9169 133,9858 142,6373 147,1367 147,1367 142,6373 133,9858 121,9169 107,5729 
104,246 112,3002 119,1025 124,0013 126,5584 126,5584 124,0013 119,1025 112,3002 104,246 
101,3678 103,9638 106,1598 107,7443 108,5727 108,5727 107,7443 106,1598 103,9638 101,3678 
 
 
O gráfico de temperaturas para este caso é mostrado abaixo. 
 
 
 
 
 
 
 
 
 
 
 
 
 
Percebe-se que neste caso a barra ganha calor em seu Norte e o fluxo de calor é positivo de cima 
para baixo. O gráfico abaixo mostra as isotermas (para vários valores de ). 
 
 
 
 
 
 
 
 
 
 
 
 
Solução Numérica: 
A solução numérica construída baseada no algoritmo TDMA Line-by-Line, cuja condições de entrada 
são as mesmas para a análise teórica: 
Lx = 1; %Comprimento da placa em x (m) 
Ly = 1; %Altura da placa em y (m) 
Nx = 10; %Número de volumes de controle em x 
Ny = 10; %Número de volumes de controle em y 
k = 1; %Condutividade térmica do material da placa (W/m°C) 
tolerancia = 1e-6; %Condição de parada 
CC_North = 600; %Condição de Contorno de Temperatura (°C) p/ Borda 
Superior 
CC_South = 100; %Condição de Contorno de Temperatura (°C) p/ Borda 
Inferior 
CC_West = 100; %Condição de Contorno de Temperatura (°C) p/ Borda Esquerda 
CC_East = 100; %Condição de Contorno de Temperatura (°C) p/ Borda Direita 
 
 
A tolerância para o resíduo de solução e condição de parada foi estabelecido em 1.10-6. 
A solução foi obtida depois de 47 iterações, 
Para a solução numérica, a matriz de temperaturas calculadas (excluindo as temperaturas nos 
contornos), onde 0.05 ≤ 𝑥 ≤ 0.95 e 0.05 ≤ 𝑦 ≤ 0.95, é: 
348,6142 482,3726 525,6983 542,1596 548,1959 548,1959 542,1596 525,6983 482,3726 348,6142 
209,3126 337,5506 403,9593 436,9037 450,6242 450,6242 436,9037 403,9593 337,5506 209,3126 
160,3982 254,5578 315,6848 350,8717 366,7729 366,7729 350,8717 315,6848 254,5578 160,3982 
138,1206 204,5977 253,3502 284,1255 298,8228 298,8228 284,1255 253,3502 204,5977 138,1206 
125,6071 172,3622 208,9927 233,4573 245,57 245,57 233,4573 208,9927 172,3622 125,6071 
117,5525 150,2514 176,8013 195,1408 204,43 204,43 195,1408 176,8013 150,2514 117,5525 
111,904 134,2897 152,8202 165,8745 172,5792 172,5792 165,8745 152,8202 134,2897 111,904 
107,6779 122,1831 134,3152 142,958 147,4331 147,4331 142,958 134,3152 122,1831 107,6779 
104,3026 112,4494 119,2998 124,2089 126,7622 126,7622 124,2089 119,2998 112,4494 104,3026 
101,3858 104,0122 106,2255 107,8158 108,6445 108,6445 107,8158 106,2255 104,0121 101,3858 
 
O gráfico de temperaturas para este caso é mostrado abaixo. 
 
 
 
 
 
 
 
 
 
 
 
 
Percebe-se que neste caso a barra ganha calor em seu Norte e o fluxo de calor é positivo de cima 
para baixo. O gráfico abaixo mostra as isotermas (para vários valores de ). 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
O erro de entre as soluções teórica e numérica é apresentado na matriz abaixo: 
1,80E-02 1,10E+01 4,18E+00 1,74E+00 1,03E+00 1,03E+00 1,74E+00 4,18E+00 1,10E+01 1,80E-02 
1,09E+01 1,49E-01 2,67E+00 2,32E+00 1,93E+00 1,93E+00 2,32E+00 2,67E+00 1,49E-01 1,09E+01 
4,01E+00 2,20E+00 3,29E-01 1,32E+00 1,56E+00 1,56E+00 1,32E+00 3,29E-01 2,20E+00 4,01E+00 
1,49E+00 1,69E+00 5,41E-01 3,84E-01 8,00E-01 8,00E-01 3,84E-01 5,41E-01 1,69E+00 1,49E+00 
6,37E-01 1,06E+00 6,75E-01 1,40E-01 1,77E-01 1,77E-01 1,40E-01 6,75E-01 1,06E+00 6,37E-01 
3,18E-01 6,60E-01 5,92E-01 3,52E-01 1,77E-01 1,77E-01 3,52E-01 5,92E-01 6,60E-01 3,18E-01 
1,79E-01 4,22E-01 4,63E-01 3,84E-01 3,08E-01 3,08E-01 3,84E-01 4,63E-01 4,22E-01 1,79E-01 
1,05E-01 2,66E-01 3,29E-01 3,21E-01 2,96E-01 2,96E-01 3,21E-01 3,29E-01 2,66E-01 1,05E-01 
5,67E-02 1,49E-01 1,97E-01 2,08E-01 2,04E-01 2,04E-01 2,08E-01 1,97E-01 1,49E-01 5,67E-02 
1,80E-02 4,83E-02 6,57E-02 7,15E-02 7,18E-02 7,18E-02 7,15E-02 6,57E-02 4,83E-02 1,80E-02 
 
O gráfico abaixo apresenta a diferença de soluções 
 
 
 
 
 
 
 
 
 
 
 
Caso 2: Problema proposto em sala, com temperaturas prescritas e fixas nas quatro bordas 
 
 
 
 
 
 
 
 
 
 
Condições de Contorno: 
𝑇(0, 𝑦) = 600°𝐶 𝑇(1, 𝑦) = 300°𝐶 𝑇(𝑥, 0) = 600°𝐶 𝑇(𝑥, 1) = 300°𝐶 
 
Utilizando o algoritmo TDMA Line-by-Line, foram extraídos resultados para o campo de temperaturas 
T(x,y) na placa. Os resultados foram plotados para 0.05 ≤ 𝑥 ≤ 0.95 e 0.05 ≤ 𝑦 ≤ 0.95, excluindo os 
contornos. 
Lx = 1; %Comprimento da placa em x (m) 
Ly = 1; %Altura da placa em y (m) 
Nx = 10; %Número de volumes de controle em x 
Ny = 10; %Número de volumes de controle em y 
k = 1; %Condutividade térmica do material da placa (W/m°C)tolerancia = 1e-6; %Condição de parada 
CC_North = 300; %Condição de Contorno de Temperatura (°C) p/ Borda 
Superior 
CC_South = 600; %Condição de Contorno de Temperatura (°C) p/ Borda 
Inferior 
CC_West = 600; %Condição de Contorno de Temperatura (°C) p/ Borda Esquerda 
CC_East = 300; %Condição de Contorno de Temperatura (°C) p/ Borda Direita 
 
Foram necessárias 37 iterações para a convergência da solução. 
 
 
 
 
 
 
 
 
 
 
 
 
 
As isotermas são mostradas na figura abaixo: 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Caso 3: Problema proposto em sala, com temperaturas prescritas e fixas em 2 bordas e 
condição de contorno de fluxo de calor (conhecidos os gradientes de temperatura) nas outras 
2 bordas. 
Foi utilizado o método do ponto “Ghost” para análise dos gradientes de temperatura, como se segue: 
Norte: 
𝜕𝑇
𝜕𝑦
|
𝑦=1
= 
𝑇∞ − 𝑇
∗ (𝑥, 1 −
∆𝑦
2 )
∆𝑦
 
𝑇(𝑥, 1) = 𝑇∞ − (
∆𝑦
2
)
𝜕𝑇
𝜕𝑦
|
𝑦=1
 
∆𝑦= 𝐿𝑦 𝑁𝑦⁄ 
 
Onde T(x,1) é a condição de contorno para o leste. Ly=1m e Ny=10 (número de nós em y). 
 
 
Leste: 
𝜕𝑇
𝜕𝑥
|
𝑥=1
= 
𝑇∞ − 𝑇
∗ (1 −
∆𝑥
2 , 𝑦)
∆𝑥
 
𝑇(1, 𝑦) = 𝑇∞ − (
∆𝑥
2
)
𝜕𝑇
𝜕𝑥
|
𝑥=1
 
∆𝑥= 𝐿𝑥 𝑁𝑥⁄ 
 
Onde T(1,y) é a condição de contorno para o leste. Lx=1m e Nx=10 (número de nós em x). 
 
 
 
1ª Situação) Temperaturas Fixas no Oeste e Sul e Calor Adicionado no Norte e Leste 
 
 
 
 
 
 
 
 
 
 
 
Condições de Contorno: 
𝑇(0, 𝑦) = 600°𝐶 𝑇(𝑥, 0) = 600°𝐶 
𝑘
𝜕𝑇
𝜕𝑦
|
𝑦=1
= 𝑞𝑁
" 
𝑘 = 1
𝑊
𝑚°𝐶
 
𝜕𝑇
𝜕𝑦
|
𝑦=1
> 0 
𝑘
𝜕𝑇
𝜕𝑥
|
𝑥=1
= 𝑞𝐸
" 
𝑘 = 1
𝑊
𝑚°𝐶
 
𝜕𝑇
𝜕𝑥
|
𝑥=1
> 0 
 
Para este problema, estipulou-se: 
𝜕𝑇
𝜕𝑦
|
𝑦=1
= 100
°𝐶
𝑚
 
𝜕𝑇
𝜕𝑥
|
𝑥=1
= 100
°𝐶
𝑚
 
 
Assim, 
𝑞𝑁
" = 100𝑊 𝑞𝐸
" = 100𝑊 
 
O calor entra na placa pelo norte e leste a uma taxa de 100W considerando um gradiente positivo de 
temperatura em x e y (para x=1m e y=1m). 
Os dados de entrada do algoritmo são mostrados abaixo: 
Lx = 1; %Comprimento da placa em x (m) 
Ly = 1; %Altura da placa em y (m) 
Nx = 10; %Número de volumes de controle em x 
Ny = 10; %Número de volumes de controle em y 
k = 1; %Condutividade térmica do material da placa (W/m°C) 
tolerancia = 1e-6; %Condição de parada 
CC_South = 600; %Condição de Contorno de Temperatura (°C) p/ Sul 
CC_West = 600; %Condição de Contorno de Temperatura (°C) p/ Oeste 
Ghost_East=700; %Temperatura(°C) de chute infinito - East 
Ghost_North=700; %Temperatura(°C) de chute infinito - North 
dPhidx_East=100; %Gradiente de Temperatura em x-East (°C/m) 
dPhidy_North=100; %Gradiente de Temperatura em y-North (°C/m) 
 
 
Foram necessárias 245 iterações para a convergência da solução. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Percebe-se que há um aumento de temperatura na placa, a partir de 600°C, à medida que se aproxima 
do canto superior direito. Isso se deve à adição de calor no norte e leste. 
 
 
 
 
 
 
 
 
 
 
 
2ª Situação) Temperaturas Fixas no Oeste e Sul e Condição Adiabática no Norte e Leste 
 
 
 
 
 
 
 
 
 
 
 
 
Condições de Contorno: 
𝑇(0, 𝑦) = 600°𝐶 𝑇(𝑥, 0) = 600°𝐶 
𝑘
𝜕𝑇
𝜕𝑦
|
𝑦=1
= 𝑞𝑁
" 
𝑘 = 1
𝑊
𝑚°𝐶
 
𝜕𝑇
𝜕𝑦
|
𝑦=1
= 0 
𝑘
𝜕𝑇
𝜕𝑥
|
𝑥=1
= 𝑞𝐸
" 
𝑘 = 1
𝑊
𝑚°𝐶
 
𝜕𝑇
𝜕𝑥
|
𝑥=1
= 0 
 
Assim, 
𝑞𝑁
" = 0 𝑞𝐸
" = 0 
 
Ou seja, 
 
 
Norte: 
𝑇∞ = 𝑇
∗ (𝑥, 1 −
∆𝑦
2
) 
𝑇(𝑥, 1) = 𝑇∞ 
Leste: 
𝑇∞ = 𝑇
∗ (1 −
∆𝑥
2
, 𝑦) 
𝑇(1, 𝑦) = 𝑇∞ 
 
A temperatura, neste caso, se mantém constante ao longo da placa. 
Os dados de entrada do algoritmo TDMA Line-by-Line são: 
Lx = 1; %Comprimento da placa em x (m) 
Ly = 1; %Altura da placa em y (m) 
Nx = 10; %Número de volumes de controle em x 
Ny = 10; %Número de volumes de controle em y 
k = 1; %Condutividade térmica do material da placa (W/m°C) 
tolerancia = 1e-6; %Condição de parada 
CC_South = 600; %Condição de Contorno de Temperatura (°C) p/ Sul 
CC_West = 600; %Condição de Contorno de Temperatura (°C) p/ Oeste 
Ghost_East=600; %Temperatura(°C) de chute infinito - East 
Ghost_North=600; %Temperatura(°C) de chute infinito - North 
dPhidx_East=0; %Gradiente de Temperatura em x-East (°C/m) 
dPhidy_North=0; %Gradiente de Temperatura em y-North (°C/m) 
 
Foram necessárias 239 iterações para a convergência da solução. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Percebe-se que a temperatura permanece constante ao longo da placa, com um erro numérico de 
2,37E-04 °C entre a maior e menor temperaturas calculadas. A matriz de temperaturas da solução é 
mostrada abaixo. 
599,9999833 600 599,9999 599,9999 599,9999 599,9998 599,9998 599,9998 599,9998 599,9998 
599,9999846 600 599,9999 599,9999 599,9999 599,9998 599,9998 599,9998 599,9998 599,9998 
599,999986 600 599,9999 599,9999 599,9999 599,9999 599,9998 599,9998 599,9998 599,9998 
599,9999875 600 599,9999 599,9999 599,9999 599,9999 599,9999 599,9998 599,9998 599,9998 
599,9999892 600 599,9999 599,9999 599,9999 599,9999 599,9999 599,9999 599,9999 599,9998 
599,9999911 600 600 599,9999 599,9999 599,9999 599,9999 599,9999 599,9999 599,9999 
599,999993 600 600 600 599,9999 599,9999 599,9999 599,9999 599,9999 599,9999 
599,999995 600 600 600 600 599,9999 599,9999 599,9999 599,9999 599,9999 
599,999997 600 600 600 600 600 600 600 600 600 
599,999999 600 600 600 600 600 600 600 600 600 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
3ª Situação) Temperaturas Fixas no Oeste e Sul e Calor Perdido no Norte e Leste 
 
 
 
 
 
 
 
 
 
 
 
Condições de Contorno: 
𝑇(0, 𝑦) = 600°𝐶 𝑇(𝑥, 0) = 600°𝐶 
𝑘
𝜕𝑇
𝜕𝑦
|
𝑦=1
= 𝑞𝑁
" 
𝑘 = 1
𝑊
𝑚°𝐶
 
𝜕𝑇
𝜕𝑦
|
𝑦=1
< 0 
𝑘
𝜕𝑇
𝜕𝑥
|
𝑥=1
= 𝑞𝐸
" 
𝑘 = 1
𝑊
𝑚°𝐶
 
𝜕𝑇
𝜕𝑥
|
𝑥=1
< 0 
 
Para este problema, estipulou-se: 
𝜕𝑇
𝜕𝑦
|
𝑦=1
= −100
°𝐶
𝑚
 
𝜕𝑇
𝜕𝑥
|
𝑥=1
= −100
°𝐶
𝑚
 
 
Assim, 
𝑞𝑁
" = −100𝑊 𝑞𝐸
" = −100𝑊 
 
O calor é perdido pela placa através do norte e leste a uma taxa de 100W considerando um gradiente 
negativo de temperatura em x e y (para x=1m e y=1m). 
Os dados de entrada do algoritmo são mostrados abaixo: 
Lx = 1; %Comprimento da placa em x (m) 
Ly = 1; %Altura da placa em y (m) 
Nx = 10; %Número de volumes de controle em x 
Ny = 10; %Número de volumes de controle em y 
k = 1; %Condutividade térmica do material da placa (W/m°C) 
tolerancia = 1e-6; %Condição de parada 
CC_South = 600; %Condição de Contorno de Temperatura (°C) p/ Sul 
CC_West = 600; %Condição de Contorno de Temperatura (°C) p/ Oeste 
Ghost_East=500; %Temperatura(°C) de chute infinito - East 
Ghost_North=500; %Temperatura(°C) de chute infinito - North 
dPhidx_East=-100; %Gradiente de Temperatura em x-East (°C/m) 
dPhidy_North=-100; %Gradiente de Temperatura em y-North (°C/m) 
 
 
Foram necessárias 228 iterações para a convergência da solução. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
As isotermas são mostradas abaixo. 
 
 
 
 
 
 
 
 
 
 
 
 
 
Percebe-se que há um redução de temperaturana placa, a partir de 600°C, à medida que se aproxima 
do canto superior direito. Isso se deve ao calor perdido no norte e leste. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Anexos 
 
Script Solucao_Teorica.m  Usado para caso 1 
%Solução Teórica (Fonte: Incropera Introduction Heat Transfer 6th) 
 
% Condução de calor estacionária em placa (2D) 
 
% Aluno: Renner Egalon Pereira (PGMEC) 
% Disciplina: Volumes Finitos 
 
% Abril/2017 
 
% Condições de Contorno Previstas no Código: 
 
% Temperatura Fixas e Iguais no Leste, Sul e Oeste 
% Temperatura Fixa no Norte 
 
%______________________________________________________________________ 
 
close all 
clear 
clc 
 
L=1; 
W=1; 
T1 = 100; 
T2 = 600; 
n = 1; 
Nx=10; 
Ny=10; 
Dx=L/Nx; 
Dy=W/Ny; 
 
Temp=zeros(Ny,Nx); 
Theta=zeros(Ny,Nx); 
Theta_Old=zeros(Ny,Nx); 
 
 x=Dx/2:Dx:L-Dx/2; 
 y=Dy/2:Dy:W-Dy/2; 
 
 
for n=1:2:99 
 for i = 1:Nx 
 for j = 1:Nx 
 
Theta_Old(j,i) = ((((1)^(n+1))+1)/n)*sin(n*pi*x(i)/L)*(sinh(n*pi*y(j)/L)/sinh(n*pi*W/L)); 
 
 end 
 end 
 
 Theta=Theta+Theta_Old; 
 
end 
 
Theta = (2/pi).*Theta; 
 
T = Theta.*(T2 - T1) + T1; 
 
 
for i=1:Ny 
 
 for j=1:Nx 
 
Temp(Ny-(i-1),j)=T(i,j); 
 
 end 
 
end 
 
%grafico 
 
contourf(x,y,T); 
colormap(jet); 
colorbar; 
 
% xlabel 
xlabel({'Posição x (m)'},'FontSize',12); 
 
% title 
title({'Distribuição de Temperatura na Placa (°C)'},'FontSize',12); 
 
% ylabel 
ylabel({'Posição y (m)'},'FontSize',12); 
 
 
 
 
 
 
 
 
 
 
 
 
 
Script Solucao_Teorica.m  Usado para casos 1 e 2 
%Condução de Calor em Placa - 2D - Material com propriedades constantes 
 
% TDMA Line-by-Line Solver para condução de calor estacionária 2D 
 
% Aluno: Renner Egalon Pereira (PGMEC) 
% Disciplina: Volumes Finitos 
 
% Abril/2017 
 
% Condições de Contorno Previstas no Código: 
 
% Temperatura Fixa no Oeste 
% Temperatura Fixa no Sul 
% Temperatura Fixa no Leste 
% Temperatura Fixa no Norte 
 
%______________________________________________________________________ 
 
 
close all 
clear 
clc 
 
%%%%%%%%%%%%%%%%%%%%%%%%%% Dados de Entrada %%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
global a b c d 
 
Lx = 1; %Comprimento da placa em x (m) 
Ly = 1; %Altura da placa em y (m) 
Nx = 10; %Número de volumes de controle em x 
Ny = 10; %Número de volumes de controle em y 
k = 1; %Condutividade térmica do material da placa (W/m°C) 
tolerancia = 1e-6; %Condição de parada 
CC_North = 300; %Condição de Contorno de Temperatura (°C) p/ Borda Superior 
CC_South = 600; %Condição de Contorno de Temperatura (°C) p/ Borda Inferior 
CC_West = 600; %Condição de Contorno de Temperatura (°C) p/ Borda Esquerda 
CC_East = 300; %Condição de Contorno de Temperatura (°C) p/ Borda Direita 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
deltaxw = Lx/Nx; 
deltaxe = deltaxw; 
 
deltayn = Ly/Ny; 
deltays = deltayn; 
 
Dx = deltaxw; 
Dy = deltayn; 
 
%Definindo as Matrizes 
 
 
Phi=zeros(Ny,Nx); 
Phi_old=zeros(Ny,Nx); 
Temp=zeros(Ny+2,Nx+2); 
N=zeros(Ny,Nx); 
 
%Temperaturas Iniciais de Entrada no Programa - Chutes Iniciais 
 
 
for i=1:Ny 
 
 for j=1:Nx 
 
 Phi(i,j) = 450; 
 
 end 
 
end 
 
%Cálculo de AW,AE,AN,AS,AP 
 
for i=1:Ny 
 
 for j=1:Nx 
 
AW(i,j) = (k*Dy)/deltaxw; 
AE(i,j) = (k*Dy)/deltaxe; 
AS(i,j) = (k*Dx)/deltays; 
AN(i,j) = (k*Dx)/deltayn; 
 
AP(i,j) = AW(i,j)+AE(i,j)+AS(i,j)+AN(i,j); 
 
 end 
 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%%%%%%%%% LOOP %%%%%%%%%%%%%%%%% LOOP %%%%%%%%%%%%%%%%% LOOP %%%%%%% 
 
residuo=Nx*Ny; 
 
passo=0; 
 
while residuo > tolerancia 
 
 passo=passo+1; 
 
 Phi_old = Phi; 
 
%%%%%%%%%%%%%%%% VARREDURA VERTICAL - DIREÇÃO i %%%%%%%%%%%%%%%%%%%%%%% 
 
for j=1:Nx 
 
 for i=1:Ny 
 
 %%CANTO SUPERIOR ESQUERDO --> Phi(1,1) 
 
 if i==1 & j==1 
 
 a(i) = 2*AW(1,1) + 2*AN(1,1) + AE(1,1) + AS(1,1); 
 b(i) = AS(1,1); 
 c(i)= 0; 
 d(i) = 2*AW(1,1)*CC_West + 2*AN(1,1)*CC_North + AE(1,1)*Phi(1,2); 
 
 %COLUNA ESQUERDA (i=2 até i=Ny-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny) & j==1 
 
 a(i) = 2*AW(i,1)+AE(i,1)+AS(i,1)+AN(i,1); 
 b(i) = AS(i,1); 
 c(i)= AN(i,1); 
 d(i) = 2*AW(i,1)*CC_West+AE(i,1)*Phi(i,2); 
 
 %CANTO INFERIOR ESQUERDO --> Phi(Ny,1)%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==Ny & j==1 
 
 a(i) = 2*AW(Ny,1) + 2*AS(Ny,1) + AE(Ny,1) + AN(Ny,1); 
 b(i) = 0; 
 c(i)= AN(Ny,1); 
 d(i) = 2*AW(Ny,1)*CC_West + 2*AS(Ny,1)*CC_South + AE(Ny,1)*Phi(Ny,2); 
 
 %LINHA SUPERIOR (j=2 até j=Nx-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==1 & (j~=1 & j~=Nx) 
 
 a(i) = 2*AN(1,j)+AE(1,j)+AW(1,j)+AS(1,j); 
 b(i) = AS(1,j); 
 c(i)= 0; 
 d(i) = 2*AN(1,j)*CC_North + AW(1,j)*Phi(1,j-1) + AE(1,j)*Phi(1,j+1); 
 
 %MEIO(j=2 até j=Nx-1 / i=2 até i=Ny-1 )%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny)&(j~=1 & j~=Nx) 
 
 b(i) = AS(i,j); 
 c(i) = AN(i,j); 
 a(i) = AP(i,j); 
 d(i) = AW(i,j)*Phi(i,j-1) + AE(i,j)*Phi(i,j+1); 
 
 %LINHA INFERIOR (j=2 até j=Nx-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==Ny & (j~=1 & j~=Nx) 
 
 a(i) = 2*AS(Ny,j)+AE(Ny,j)+AW(Ny,j)+AN(Ny,j); 
 b(i) = 0; 
 c(i)= AN(Ny,j); 
 d(i) = 2*AS(Ny,j)*CC_South + AW(Ny,j)*Phi(Ny,j-1) + AE(Ny,j)*Phi(Ny,j+1); 
 
 %CANTO SUPERIOR DIREITO --> Phi(1,Nx) 
 
 elseif i==1 & j==Nx 
 
 a(i) = 2*AE(1,Nx) + 2*AN(1,Nx) + AW(1,Nx) + AS(1,Nx); 
 b(i) = AS(1,Nx); 
 c(i)= 0; 
 d(i) = 2*AE(1,Nx)*CC_East + 2*AN(1,Nx)*CC_North + AW(1,Nx)*Phi(1,Nx-1); 
 
 %COLUNA DIREITA (i=2 até i=Ny-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny) & j==Nx 
 
 a(i) = 2*AE(i,Nx)+AW(i,Nx)+AS(i,Nx)+AN(i,Nx); 
 b(i) = AS(i,Nx); 
 c(i)= AN(i,Nx); 
 d(i) = 2*AE(i,Nx)*CC_East+AW(i,Nx)*Phi(i,Nx-1); 
 
 %CANTO INFERIOR DIREITO --> Phi(Ny,Nx) 
 
 elseif i==Ny & j==Nx 
 
 a(i) = 2*AE(Ny,Nx) + 2*AS(Ny,Nx) + AW(Ny,Nx) + AN(Ny,Nx); 
 b(i) = 0; 
 c(i)= AN(Ny,Nx); 
 d(i) = 2*AE(Ny,Nx)*CC_East + 2*AS(Ny,Nx)*CC_South + AW(Ny,Nx)*Phi(Ny,Nx-1); 
 
 end 
 
 end 
 
 TDMA; %Algoritmo de Thomas 
 
 for i=1:Ny 
 
 Phi(i,j)=X(i); 
 end 
 
end 
 
%%%%%%%%%%%%%%%% VARREDURA HORIZONTAL - DIREÇÃO j %%%%%%%%%%%%%%%%%%%% 
 
for i=1:Ny 
 
 for j=1:Nx 
 
 %%CANTO SUPERIOR ESQUERDO --> Phi(1,1) 
 
 if i==1 & j==1 
 
 a(j) = 2*AW(1,1) + 2*AN(1,1) + AE(1,1) + AS(1,1); 
 b(j) = AE(1,1); 
 c(j)= 0; 
 d(j) = 2*AW(1,1)*CC_West + 2*AN(1,1)*CC_North + AS(1,1)*Phi(2,1); 
 
 %LINHA SUPERIOR (j=2 até j=Nx-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==1 & (j~=1 & j~=Nx) 
 
 a(j) = 2*AN(1,j)+AE(1,j)+AW(1,j)+AS(1,j);b(j) = AE(1,j); 
 c(j)= AW(1,j); 
 d(j) = 2*AN(1,j)*CC_North + AS(1,j)*Phi(2,j); 
 
 %CANTO SUPERIOR DIREITO --> Phi(1,Nx) 
 
 elseif i==1 & j==Nx 
 
 a(j) = 2*AE(1,Nx) + 2*AN(1,Nx) + AW(1,Nx) + AS(1,Nx); 
 b(j) = 0; 
 c(j)= AW(1,Nx); 
 d(j)= 2*AE(1,Nx)*CC_East + 2*AN(1,Nx)*CC_North + AS(1,Nx)*Phi(2,Nx); 
 
 %COLUNA ESQUERDA (i=2 até i=Ny-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny) & j==1 
 
 a(j) = 2*AW(i,1)+AE(i,1)+AS(i,1)+AN(i,1); 
 b(j) = AE(i,1); 
 c(j)= 0; 
 d(j) = 2*AW(i,1)*CC_West + AN(i,1)*Phi(i-1,1)+ AS(i,1)*Phi(i+1,1); 
 
 %MEIO(j=2 até j=Nx-1 / i=2 até i=Ny-1 )%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny)&(j~=1 & j~=Nx) 
 
 b(j) = AE(i,j); 
 c(j) = AW(i,j); 
 a(j) = AP(i,j); 
 d(j) = AS(i,j)*Phi(i+1,j) + AN(i,j)*Phi(i-1,j); 
 
 %COLUNA DIREITA (i=2 até i=Ny-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny) & j==Nx 
 
 a(j) = 2*AE(i,Nx)+AW(i,Nx)+AS(i,Nx)+AN(i,Nx); 
 b(j) = 0; 
 c(j)= AW(i,Nx); 
 d(j) = 2*AE(i,Nx)*CC_East + AN(i,Nx)*Phi(i-1,Nx)+ AS(i,1)*Phi(i+1,Nx); 
 
 %CANTO INFERIOR ESQUERDO --> Phi(Ny,1)%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==Ny & j==1 
 
 a(j) = 2*AW(Ny,1) + 2*AS(Ny,1) + AE(Ny,1) + AN(Ny,1); 
 b(j) = AE(Ny,1); 
 c(j)= 0; 
 d(j) = 2*AW(Ny,1)*CC_West + 2*AS(Ny,1)*CC_South + AN(Ny,1)*Phi(Ny-1,1); 
 
 %LINHA INFERIOR (j=2 até j=Nx-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==Ny & (j~=1 & j~=Nx) 
 
 a(j) = 2*AS(Ny,j)+AE(Ny,j)+AW(Ny,j)+AN(Ny,j); 
 b(j) = AE(1,j); 
 c(j)= AW(1,j); 
 d(j) = 2*AS(Ny,j)*CC_South + AN(Ny,j)*Phi(Ny-1,j); 
 
 %CANTO INFERIOR DIREITO --> Phi(Ny,Nx) 
 
 elseif i==Ny & j==Nx 
 
 a(j) = 2*AE(Ny,Nx) + 2*AS(Ny,Nx) + AW(Ny,Nx) + AN(Ny,Nx); 
 b(j) = 0; 
 c(j)= AW(Ny,Nx); 
 d(j) = 2*AE(Ny,Nx)*CC_East + 2*AS(Ny,Nx)*CC_South + AN(Ny,Nx)*Phi(Ny-1,Nx); 
 
 end 
 
 end 
 
 TDMA; %Algoritmo de Thomas 
 
 for j=1:Nx 
 
 Phi(i,j)=X(j); 
 end 
 
end 
 
residuo =max(abs(Phi - Phi_old)); 
 
end 
%%%%%%%%%%%%%%%%%%%%% FIM LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 
% for i=1:Ny+2 
% 
% for j=1:Nx+2 
% 
% if (i~=1 & i~=Ny+2) & j==1; 
% 
% Temp(i,j)=CC_West; 
% 
% elseif i==Ny+2 & (j~=1 & j~=Nx+2); 
% 
% Temp(Ny+2,j)=CC_South; 
% 
% elseif (i~=1 & i~=Ny+2) & j==Nx+2; 
% 
% Temp(i,j)=CC_East; 
% 
% elseif i==1 & (j~=1 & j~=Nx+2); 
% 
% Temp(i,j)=CC_North; 
% 
% elseif (i~=1 & i~=Ny+2)&(j~=1 & j~=Nx+2); 
% 
% Temp(i,j)=Phi(i-1,j-1); 
% 
% end 
% 
% 
% %QUINAS 
% 
% Temp(1,1)=(Temp(1,2)+Temp(2,1)+CC_West+CC_North)/4; 
% Temp(Ny+2,1)=(Temp(Ny+2,2)+Temp(Ny+1,1)+CC_West+CC_South)/4; 
% Temp(Ny+2,Nx+2)=(Temp(Ny+2,Nx+1)+Temp(Ny+1,Nx+2)+CC_South+CC_East)/4; 
% Temp(1,Nx+2)=(Temp(1,Nx+1)+Temp(2,Nx+2)+CC_East+CC_North)/4; 
% 
% end 
% 
% end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 
 
 
for i=1:Ny 
 
 for j=1:Nx 
 
N(Ny-(i-1),j)=Phi(i,j); 
 
 end 
 
end 
 
 
 
% %xlabel e ylabel 
% 
% for i = 3:Nx + 1 
% 
% x(1) = 0; 
% 
% x(2) = x(1) + ((Dx) / 2); 
% 
% x(i) = x(i - 1) + Dx; 
% 
% end 
% 
% x(Nx + 2) = x(Nx + 1) + ((Dx) / 2); 
% 
% 
% for i = 3:Ny + 1 
% 
% y(1) = 0; 
% 
% y(2) = y(1) + ((Dy) / 2); 
% 
% y(i) = y(i - 1) + Dy; 
% 
% end 
% 
% y(Ny + 2) = y(Ny + 1) + ((Dy) / 2); 
 
%grafico 
 
x=Dx/2:Dx:Lx-Dx/2; 
y=Dy/2:Dx:Ly-Dy/2; 
 
contourf(x,y,N); 
colormap(jet); 
colorbar; 
 
% xlabel 
xlabel({'Posição x (m)'},'FontSize',12); 
 
% title 
title({'Distribuição de Temperatura na Placa (°C)'},'FontSize',12); 
 
% ylabel 
ylabel({'Posição y (m)'},'FontSize',12); 
 
 
 
Script Conducao_2D_Grad.m  Usado para caso 3 
% Condução de Calor em Placa - 2D - Material com propriedades constantes 
 
% TDMA Line-by-Line Solver para condução de calor estacionária 2D 
 
% Aluno: Renner Egalon Pereira (PGMEC) 
% Disciplina: Volumes Finitos 
 
% Abril/2017 
 
% Condições de Contorno Previstas no Código: 
 
% Temperatura Fixa no Oeste 
% Temperatura Fixa no Sul 
% fluxo de calor em y no Norte 
% fluxo de calor em x no Leste 
 
%______________________________________________________________________ 
 
close all 
clear 
clc 
 
global a b c d 
 
%%%%%%%%%%%%%%%%%%%%%%%%%% Dados de Entrada %%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
Lx = 1; %Comprimento da placa em x (m) 
Ly = 1; %Altura da placa em y (m) 
Nx = 10; %Número de volumes de controle em x 
Ny = 10; %Número de volumes de controle em y 
k = 1; %Condutividade térmica do material da placa (W/m°C) 
tolerancia = 1e-6; %Condição de parada 
CC_South = 600; %Condição de Contorno de Temperatura (°C) p/ Sul 
CC_West = 600; %Condição de Contorno de Temperatura (°C) p/ Oeste 
Ghost_East=500; %Temperatura(°C) de chute infinito - East 
Ghost_North=500; %Temperatura(°C) de chute infinito - North 
dPhidx_East=-100; %Gradiente de Temperatura em x-East (°C/m) 
dPhidy_North=-100; %Gradiente de Temperatura em y-North (°C/m) 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 
deltaxw = Lx/Nx; 
deltaxe = deltaxw; 
 
deltayn = Ly/Ny; 
deltays = deltayn; 
 
Dx = deltaxw; 
Dy = deltayn; 
 
%Definindo as Matrizes 
 
 
Phi=zeros(Ny,Nx); 
Phi_old=zeros(Ny,Nx); 
Temp=zeros(Ny+2,Nx+2); 
N=zeros(Ny,Nx); 
 
%Temperaturas Iniciais de Entrada no Programa - Chutes Iniciais 
 
 
for i=1:Ny 
 
 for j=1:Nx 
 
 Phi(i,j) = 450; 
 
 end 
 
end 
 
 
%Chute Temperaturas Ghost - East e North 
 
 for i=1:Ny 
 
 Phighost_East(i,Nx) = Ghost_East; 
 
 end 
 
 
 for j=1:Nx 
 
 Phighost_North(1,j) = Ghost_North; 
 
 end 
 
%Cálculo de AW,AE,AN,AS,AP 
 
for i=1:Ny 
 
 for j=1:Nx 
 
AW(i,j) = (k*Dy)/deltaxw; 
AE(i,j) = (k*Dy)/deltaxe; 
AS(i,j) = (k*Dx)/deltays; 
AN(i,j) = (k*Dx)/deltayn; 
 
AP(i,j) = AW(i,j)+AE(i,j)+AS(i,j)+AN(i,j); 
 
 end 
 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%%%%%%%%% LOOP %%%%%%%%%%%%%%%%% LOOP %%%%%%%%%%%%%%%%% LOOP %%%%%%% 
 
residuo=Nx*Ny; 
 
passo=0; 
 
while residuo > tolerancia 
 
 passo=passo+1; 
 
 Phi_old = Phi; 
 
 %Definição das CC's East e North 
 
 for i=1:Ny 
 
 CC_East(i,Nx) = Phighost_East(i,Nx) - dPhidx_East*(Dx/2); 
 
 end 
 
 for j=1:Nx 
 
 CC_North(1,j) = Phighost_North(1,j) - dPhidy_North*(Dy/2); 
 
 end 
 
%%%%%%%%%%%%%%%% VARREDURA VERTICAL - DIREÇÃO i %%%%%%%%%%%%%%%%%%%%%%% 
 
for j=1:Nx 
 
 for i=1:Ny 
 
 %%CANTO SUPERIOR ESQUERDO --> Phi(1,1) 
 
 if i==1 & j==1 
 
 a(i) = 2*AW(1,1) + 2*AN(1,1) + AE(1,1) + AS(1,1); 
 b(i) = AS(1,1); 
 c(i)= 0; 
 d(i) = 2*AW(1,1)*CC_West + 2*AN(1,1)*CC_North(1,j) + AE(1,1)*Phi(1,2);%COLUNA ESQUERDA (i=2 até i=Ny-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny) & j==1 
 
 a(i) = 2*AW(i,1)+AE(i,1)+AS(i,1)+AN(i,1); 
 b(i) = AS(i,1); 
 c(i)= AN(i,1); 
 d(i) = 2*AW(i,1)*CC_West+AE(i,1)*Phi(i,2); 
 
 %CANTO INFERIOR ESQUERDO --> Phi(Ny,1)%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==Ny & j==1 
 
 a(i) = 2*AW(Ny,1) + 2*AS(Ny,1) + AE(Ny,1) + AN(Ny,1); 
 b(i) = 0; 
 c(i)= AN(Ny,1); 
 d(i) = 2*AW(Ny,1)*CC_West + 2*AS(Ny,1)*CC_South + AE(Ny,1)*Phi(Ny,2); 
 
 %LINHA SUPERIOR (j=2 até j=Nx-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==1 & (j~=1 & j~=Nx) 
 
 a(i) = 2*AN(1,j)+AE(1,j)+AW(1,j)+AS(1,j); 
 b(i) = AS(1,j); 
 c(i)= 0; 
 d(i) = 2*AN(1,j)*CC_North(1,j) + AW(1,j)*Phi(1,j-1) + AE(1,j)*Phi(1,j+1); 
 
 %MEIO(j=2 até j=Nx-1 / i=2 até i=Ny-1 )%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny)&(j~=1 & j~=Nx) 
 
 b(i) = AS(i,j); 
 c(i) = AN(i,j); 
 a(i) = AP(i,j); 
 d(i) = AW(i,j)*Phi(i,j-1) + AE(i,j)*Phi(i,j+1); 
 
 %LINHA INFERIOR (j=2 até j=Nx-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==Ny & (j~=1 & j~=Nx) 
 
 a(i) = 2*AS(Ny,j)+AE(Ny,j)+AW(Ny,j)+AN(Ny,j); 
 b(i) = 0; 
 c(i)= AN(Ny,j); 
 d(i) = 2*AS(Ny,j)*CC_South + AW(Ny,j)*Phi(Ny,j-1) + AE(Ny,j)*Phi(Ny,j+1); 
 
 %CANTO SUPERIOR DIREITO --> Phi(1,Nx) 
 
 elseif i==1 & j==Nx 
 
 a(i) = 2*AE(1,Nx) + 2*AN(1,Nx) + AW(1,Nx) + AS(1,Nx); 
 b(i) = AS(1,Nx); 
 c(i)= 0; 
 d(i) = 2*AE(1,Nx)*CC_East(i,Nx) + 2*AN(1,Nx)*CC_North(1,j) + AW(1,Nx)*Phi(1,Nx-
1); 
 
 %COLUNA DIREITA (i=2 até i=Ny-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny) & j==Nx 
 
 a(i) = 2*AE(i,Nx)+AW(i,Nx)+AS(i,Nx)+AN(i,Nx); 
 b(i) = AS(i,Nx); 
 c(i)= AN(i,Nx); 
 d(i) = 2*AE(i,Nx)*CC_East(i,Nx)+AW(i,Nx)*Phi(i,Nx-1); 
 
 %CANTO INFERIOR DIREITO --> Phi(Ny,Nx) 
 
 elseif i==Ny & j==Nx 
 
 a(i) = 2*AE(Ny,Nx) + 2*AS(Ny,Nx) + AW(Ny,Nx) + AN(Ny,Nx); 
 b(i) = 0; 
 c(i)= AN(Ny,Nx); 
 d(i) = 2*AE(Ny,Nx)*CC_East(i,Nx) + 2*AS(Ny,Nx)*CC_South + AW(Ny,Nx)*Phi(Ny,Nx-1); 
 
 end 
 
 end 
 
 TDMA; %Algoritmo de Thomas 
 
 for i=1:Ny 
 
 Phi(i,j)=X(i); 
 end 
 
end 
 
%%% Atualização das Phi's Ghost East e North 
 
 for i=1:Ny 
 
 Phighost_East(i,Nx)=Phi(i,Nx) + dPhidx_East*Dx; 
 
 end 
 
 
 for j=1:Nx 
 
 Phighost_North(1,j)=Phi(1,j) + dPhidy_North*Dy; 
 
 end 
 
%Definição das CC's East e North 
 
 for i=1:Ny 
 
 CC_East(i,Nx) = Phighost_East(i,Nx) - dPhidx_East*(Dx/2); 
 
 end 
 
 for j=1:Nx 
 
 CC_North(1,j) = Phighost_North(1,j) - dPhidy_North*(Dy/2); 
 
 end 
 
%%%%%%%%%%%%%%%% VARREDURA HORIZONTAL - DIREÇÃO j %%%%%%%%%%%%%%%%%%%% 
 
for i=1:Ny 
 
 for j=1:Nx 
 
 %%CANTO SUPERIOR ESQUERDO --> Phi(1,1) 
 
 if i==1 & j==1 
 
 a(j) = 2*AW(1,1) + 2*AN(1,1) + AE(1,1) + AS(1,1); 
 b(j) = AE(1,1); 
 c(j)= 0; 
 d(j) = 2*AW(1,1)*CC_West + 2*AN(1,1)*CC_North(1,j) + AS(1,1)*Phi(2,1); 
 
 %LINHA SUPERIOR (j=2 até j=Nx-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==1 & (j~=1 & j~=Nx) 
 
 a(j) = 2*AN(1,j)+AE(1,j)+AW(1,j)+AS(1,j); 
 b(j) = AE(1,j); 
 c(j)= AW(1,j); 
 d(j) = 2*AN(1,j)*CC_North(1,j) + AS(1,j)*Phi(2,j); 
 
 %CANTO SUPERIOR DIREITO --> Phi(1,Nx) 
 
 elseif i==1 & j==Nx 
 
 a(j) = 2*AE(1,Nx) + 2*AN(1,Nx) + AW(1,Nx) + AS(1,Nx); 
 b(j) = 0; 
 c(j)= AW(1,Nx); 
 d(j)= 2*AE(1,Nx)*CC_East(i,Nx) + 2*AN(1,Nx)*CC_North(1,j) + AS(1,Nx)*Phi(2,Nx); 
 
 %COLUNA ESQUERDA (i=2 até i=Ny-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny) & j==1 
 
 a(j) = 2*AW(i,1)+AE(i,1)+AS(i,1)+AN(i,1); 
 b(j) = AE(i,1); 
 c(j)= 0; 
 d(j) = 2*AW(i,1)*CC_West + AN(i,1)*Phi(i-1,1)+ AS(i,1)*Phi(i+1,1); 
 
 %MEIO(j=2 até j=Nx-1 / i=2 até i=Ny-1 )%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny)&(j~=1 & j~=Nx) 
 
 b(j) = AE(i,j); 
 c(j) = AW(i,j); 
 a(j) = AP(i,j); 
 d(j) = AS(i,j)*Phi(i+1,j) + AN(i,j)*Phi(i-1,j); 
 
 %COLUNA DIREITA (i=2 até i=Ny-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif (i~=1 & i~=Ny) & j==Nx 
 
 a(j) = 2*AE(i,Nx)+AW(i,Nx)+AS(i,Nx)+AN(i,Nx); 
 b(j) = 0; 
 c(j)= AW(i,Nx); 
 d(j) = 2*AE(i,Nx)*CC_East(i,Nx) + AN(i,Nx)*Phi(i-1,Nx)+ AS(i,1)*Phi(i+1,Nx); 
 
 %CANTO INFERIOR ESQUERDO --> Phi(Ny,1)%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==Ny & j==1 
 
 a(j) = 2*AW(Ny,1) + 2*AS(Ny,1) + AE(Ny,1) + AN(Ny,1); 
 b(j) = AE(Ny,1); 
 c(j)= 0; 
 d(j) = 2*AW(Ny,1)*CC_West + 2*AS(Ny,1)*CC_South + AN(Ny,1)*Phi(Ny-1,1); 
 
 %LINHA INFERIOR (j=2 até j=Nx-1)%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 elseif i==Ny & (j~=1 & j~=Nx) 
 
 a(j) = 2*AS(Ny,j)+AE(Ny,j)+AW(Ny,j)+AN(Ny,j); 
 b(j) = AE(1,j); 
 c(j)= AW(1,j); 
 d(j) = 2*AS(Ny,j)*CC_South + AN(Ny,j)*Phi(Ny-1,j); 
 
 %CANTO INFERIOR DIREITO --> Phi(Ny,Nx) 
 
 elseif i==Ny & j==Nx 
 
 a(j) = 2*AE(Ny,Nx) + 2*AS(Ny,Nx) + AW(Ny,Nx) + AN(Ny,Nx); 
 b(j) = 0; 
 c(j)= AW(Ny,Nx); 
 d(j) = 2*AE(Ny,Nx)*CC_East(i,Nx) + 2*AS(Ny,Nx)*CC_South + AN(Ny,Nx)*Phi(Ny-1,Nx); 
 
 end 
 
 end 
 
 TDMA; %Algoritmo de Thomas 
 
 for j=1:Nx 
 
 Phi(i,j)=X(j); 
 end 
 
end 
 
%%% Atualização das Phi's Ghost East e North 
 
 for i=1:Ny 
 
 Phighost_East(i,Nx)=Phi(i,Nx) + dPhidx_East*Dx; 
 
 end 
 
 
 for j=1:Nx 
 
 Phighost_North(1,j)=Phi(1,j) + dPhidy_North*Dy; 
 
 end 
 
 
residuo =max(abs(Phi - Phi_old)); 
 
end 
%%%%%%%%%%%%%%%%%%%%% FIM LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% %%% Matriz Temperaturas (com as CC's) 
% 
% for i=1:Ny+2 
% 
% for j=1:Nx+2 
% 
% if (i~=1 & i~=Ny+2) & j==1; 
% 
% Temp(i,j)=CC_West; 
% 
% elseif i==Ny+2 & (j~=1 & j~=Nx+2); 
% 
% Temp(Ny+2,j)=CC_South; 
% 
% elseif (i~=1 & i~=Ny+2) & j==Nx+2; 
% 
% Temp(i,j)=CC_East(i-1,j-2); 
% 
% elseif i==1 & (j~=1 & j~=Nx+2); 
% 
% Temp(i,j)=CC_North(i,j-1); 
% 
% elseif (i~=1 & i~=Ny+2)&(j~=1 & j~=Nx+2); 
% 
% Temp(i,j)=Phi(i-1,j-1); 
% 
% end 
% 
% 
% %QUINAS 
% 
% Temp(1,1)=(Temp(1,2)+Temp(2,1)+CC_West+CC_North(1,1))/4; 
% Temp(Ny+2,1)=(Temp(Ny+2,2)+Temp(Ny+1,1)+CC_West+CC_South)/4; 
% Temp(Ny+2,Nx+2)=(Temp(Ny+2,Nx+1)+Temp(Ny+1,Nx+2)+CC_South+CC_East(Ny,Nx))/4; 
% Temp(1,Nx+2)=(Temp(1,Nx+1)+Temp(2,Nx+2)+CC_East(1,Nx)+CC_North(1,Nx))/4; 
% 
% end 
%% end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
for i=1:Ny 
 
 for j=1:Nx 
 
M(Ny-(i-1),j)=Phi(i,j); 
 
 end 
 
end 
 
 
 
% for i=1:Ny+2 
% 
% for j=1:Nx+2 
% 
% N((Ny+2)-(i-1),j)=Temp(i,j); 
% 
% end 
% 
% end 
 
%xlabel e ylabel 
 
% for i = 3:Nx + 1 
% 
% x(1) = 0; 
% 
% x(2) = x(1) + ((Dx) / 2); 
% 
% x(i) = x(i - 1) + Dx; 
% 
% end 
% 
% x(Nx + 2) = x(Nx + 1) + ((Dx) / 2); 
% 
% 
% for i = 3:Ny + 1 
% 
% y(1) = 0; 
% 
% y(2) = y(1) + ((Dy) / 2); 
% 
% y(i) = y(i - 1) + Dy; 
% 
% end 
% 
% y(Ny + 2) = y(Ny + 1) + ((Dy) / 2); 
 
x=Dx/2:Dx:Lx-Dx/2; 
y=Dy/2:Dx:Ly-Dy/2; 
 
 
%grafico 
 
contourf(x,y,M); 
colormap(jet); 
colorbar; 
 
% xlabel 
xlabel({'Posição x (m)'},'FontSize',12); 
 
% title 
title({'Distribuição de Temperatura na Placa (°C)'},'FontSize',12); 
 
% ylabel 
ylabel({'Posição y (m)'},'FontSize',12); 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Script TDMA.m  Algoritmo de Thomas usado em todos os casos 
% Algoritmo de Thomas (TDMA) 
 
global X 
 
N = length(a); 
 
% 1º PASSO: 
 
P(1) = b(1) / a(1); 
q(1) = d(1) / a(1); 
 
 
% 2º PASSO: 
 
for m = 2:N 
 
P(m) = b(m) / (a(m) - c(m) * P(m - 1)); 
 
q(m) = (d(m) + c(m) * q(m - 1)) / (a(m) - c(m) * P(m - 1)); 
 
end 
 
% 3º PASSO: 
 
X(N) = q(N); 
 
 
% 4º PASSO: 
 
for m = N-1:-1:1 
 
X(m) = P(m) * X(m + 1) + q(m); 
 
end