UNIVERSIDADE FEDERAL DO ESPÍRITO SANTOCENTRO TECNOLÓGICODEPARTAMENTO DE ENGENHARIA MECÂNICATRANSFERÊNCIA DE CALOR IPROF.DR. MÁRCIO FERREIRA MARTINSENTREGA 14/08/2012

Análise Transiente da Transferência de Calor em um Absorvedor Através do Método de Diferenças Finitas

Ivam Pretti

Letícia Costa Ribeiro

Márcio W. Rangel Silva

Tássio Figueira Santana

Universidade Federal do Espírito Santo, Avenida Fernando Ferrari, 514, Goiabeiras, Vitória, Espírito Santo, CEP 29075- 910

Resumo: O trabalho consiste na demonstração da resolução de problemas de transferência de calor utilizando o Método de Diferenças Finitas implicitamente. Para tanto usamos como objeto de estudo um absorvedor solar. São apresentadas as equações utilizadas e o passo a passo de como solucionar o problema. O primeiro passo consiste na discretização dos pontos da região que se procura solução, pois ao contrário da solução analítica, em que se é permitido a determinação da temperatura em qualquer ponto de interesse em um meio, uma solução numérica permite a determinação da temperatura em somente pontos discretos. Logo após, são mostradas as condições de contorno utilizadas e são feitas as análises dos volumes de controle interno e externo. Para a realização dos cálculos foi elaborado um código computacional, de forma que ao final faremos uma analise dos resultados obtidos por ele.

Palavras-chaves: Matlab, Transferência de Calor, Bidimensional, Método de Diferenças Finitas, Transiente.

1. INTRODUÇÃO

Este trabalho foi proposto como forma de avaliação para somar às notas regulares previstas na grade curricular e tem como objetivo analisar a aplicação do método das diferenças finitas na solução de problemas de difusão de calor em um tubo coletor. Para tanto, foi criado um algoritmo computacional que realizou essa análise em duas dimensões e coordenadas cilíndricas para a formulação implícita, e cujos resultados serão discutidos ao longo deste trabalho.

Tivemos como base o trabalho de conclusão de curso da aluna Letícia Jenisch Rodrigues (estudante curso de Engenharia Mecânica da Universidade Federal do Rio Grande do Sul) que realiza de forma similar a análise do tubo coletor, utilizando, porém a formulação explícita.

Ao longo deste trabalho estudaremos o processo de condução de calor transiente no coletor solar parabólico. Esses coletores, diferentemente dos painéis fotovoltaicos convencionais, utilizam espelhos concentradores de energia solar, e que de maneira simplificada funciona da seguinte forma, a superfície parabólica espelhada reflete a radiação solar para a sua linha focal onde se encontra um tubo de seção circular que contém fluido em seu interior, geralmente óleo, esse fluido trocará calor com outro fluido, água, circulando em um sistema secundário, visando à geração de vapor superaquecido.

Por fim, veremos que com auxílio de métodos computacionais podemos com certa facilidade concluir cálculos que, sem estes, levaríamos uma quantidade substancial de tempo e esforço para realizarmos, além de produzir resultados mais precisos. Além disso,

veremos que o emprego do método implícito pode frequentemente reduzir o tempo computacional em relação ao uso do método explícito.

2. COLETOR CILINDRO-PARABÓLICO

O coletor cilindro-parabólico é um concentrador de foco axial, espelho contínuo, cujo eixo de movimento pode ser tanto polar como horizontal.

Basicamente, se trata de uma superfície refletora com formato cilíndrico de seção reta parabólica no qual possui em seu foco um absorvedor.

O absorvedor, que será o objeto de estudo deste trabalho, é um tubo de metal de seção circular, cuja superfície externa é tratada de modo que a capacidade de absorção da radiação incidente seja aumentada. Geralmente é envolto por um tubo de vidro, cuja função é minimizar as perdas térmicas.

Figura 01: Esquema ilustrando um coletor cilindro-parabólico

3. FORMULAÇÕES MATEMÁTICAS

3.1. EQUAÇÕES DA DIFUSÃO DE CALOR EM COORDENADAS CILÍDRICAS

Para conhecimento das temperaturas internas do absorvedor é necessário o conhecimento da equação da difusão do calor.

1r

∂∂r (kr

∂ T∂ r )+1

r∂

∂∅ (k ∂ T∂∅ )+ ∂

∂ z (k∂ T∂ z )+q=ρ c p

∂ T∂ t

(1)

Onde T é a temperatura, r, ᴓ, z são coordenadas do espaço, t variável temporal, k

condutividade térmica, ρ a massa específica e

c p calor específico à pressão constante. Sob

condições transientes com propriedades constantes e na ausência de fonte interna, que são as condições do trabalho, a forma apropriada da equação do calor bidimensional em coordenadas cilíndricas obtida.

1r

∂∂r (kr

∂ T∂ r )+1

r∂

∂∅ (k ∂ T∂∅ )=ρ c p

∂ T∂t

(2)

Soluções analíticas para problemas transientes estão restritas a geometrias e condições contorno simples, tais como unidimensional. Soluções analíticas, bi ou tridimensional, podem se tornar complexas, assim como proposto, será usada à técnica numérica diferenças finitas na forma implícita.

3.2. MÉTODO DAS DIFERENÇAS FINITAS IMPLÍCITO: DISCRETIZAÇÃO DA EQUAÇÃO DO CALOR

Ao contrário da solução analítica, que permite a determinação da temperatura em qualquer ponto de interesse em um meio, uma solução numérica permite a determinação da temperatura em somente pontos discretos. Consequentemente a primeira etapa em qualquer analise numérica é a discretização da região que se procura solução. Para discretização define-se um domínio de pontos, conforme feito abaixo para o problema em questão.

t= p∆t

r=m∆r

ᴓ=n∆ᴓ

As variáveis t, r e ᴓ possuem passos de tamanhos distintos, indicados como ∆t, ∆r, ∆ᴓ e são identificados pelos índices p, m e n.

A Figura 02 abaixo mostra uma malha, que é um conjunto finito de pontos pertencentes ao domínio, que por sua vez são também chamados nós da malha ou nós. É possível observar que os pontos são separados por

distancias ∆r e ∆ᴓ, que necessariamente não são iguais.

Figura 02- Discretização do domínio em coordenadas cilíndricas (a) Malha (b) representação dos nós na malha

O método consiste, além da discretização do domínio, na substituição das derivadas presentes na equação diferencial por aproximações, utilizando apenas os valores numéricos da função. Essas aproximações podem ser facilmente feitas por série de Taylor, como mostrado para uma função f.

f ( x+∆ x )=f ( x )+∆ xdf (x )

dx+

(∆ x )2

2!d2 f (x )

d x2 +(∆ x )3

3 !d3 f (x)

d x3 +…+(∆ x )n

n!dn f ( x )

d xn +…

(3)

A Eq.(3) pode ser aproximada em primeira ordem da seguinte forma:

f ( x+∆ x ) ≈ f ( x )+∆ xdf (x)

dx (4)

E reformulada para a determinação de dfdx

:

df (x )dx

≈f ( x+∆ x )−f ( x )

∆ x

(5)

Para a nossa aplicação à derivada primeira é aproximada por

∂ T∂r

¿m,np+1 ≈

T m+1, np+1 −T m−1 ,n

p+1

2 ∆ r

(6)

A derivada segunda por

∂2T

∂r 2 ¿m,np+1 ≈

∂ T∂r

¿m+1 /2 , np+1 −∂ T

∂r¿m−1/2 ,n

p+1

∆ r

(7)

Os índices m e n indicam localização de cada nó da malha, como na Figura 02-b, enquanto o índice p indica a dependência temporal, sendo p o instante atual e p+1 um instante depois.

Aproximando os gradientes de primeira ordem, e colocando em termos de temperatura nos nós.

∂ T∂r

¿m+1 /2 , np+1 ≈

Tm+1 ,np+1 +T m−1 ,n

p+1

∆ r

(8)

E

∂ T∂r

¿m−1/2 ,np+1 ≈

T m,np+1+T m−1 ,n

p+1

∆ r (9)

Substituindo as equações (8) e (9) em (7) obtêm-se:

∂2T∂r 2 ¿m,n

p+1 ≈T m+1 , n

p+1 −2 Tm, np+1+Tm−1 , n

p+1

(∆ r)2

(10)

Fazendo de modo análogo para a derivada segunda de Φ e a derivada temporal, temos a equação do calor discreta implícita.

T m,np+1−α ∆ t [( 1

rTm+1 ,n

p+1 −T m−1 , np+1

(2 ∆ r )+

Tm+1 ,np+1 −2T m,n

p+1+T m−1 ,np+1

(∆ r )2 )−( 1rm

2

T m,n+1p+1 −2T m, n

p+1+T m, n−1p+1

( ∆ Φ )2 )]=T m, np

(11)

3.3. CONDIÇÕES DE CONTORNO

As condições de contorno usadas no trabalho são chamadas de segunda e terceira espécie. A condição de contorno de segunda

espécie, ou condição de Neumann, corresponde á existência de um fluxo térmico fixo ou constante em sua superfície, enquanto que a condição de contorno de terceira espécie corresponde á existência, na superfície, de um aquecimento (ou resfriamento) por convecção e é obtida a partir de um balanço de energia na superfície.

Na aplicação do balanço de energia é comumente considerado que todos os fluxos térmicos estão sendo dirigido para dentro do ponto nodal, isso porque a direção do fluxo térmico (entrando ou saindo do nó) é desconhecida.

Eent + Eg=Eacu (12)

Na qual Eent é a taxa de energia que entra

no volume de controle, Eg a taxa de geração de

energia térmica e Eacu é a taxa de aumento de

energia armazenada no interior do volume de controle. Como no nosso problema não apresenta energia térmica gerada, podemos reescrever a equação 12.

Eent=Eacu (13)

3.3.1. ANÁLISE DO VOLUME DE CONTROLE INTERNO

Efetuando o balanço de energia, equação 13, e o cálculo por diferenças finitas nos volumes de controle, Figura 03-b e c.

hóleo (T óleo−T 1 ,np+1 )+ K

∆ r(T 2 ,n

p+1−T 1 ,np+1)=ρ c p

∆ r2

T1 , np+1−T 1 ,n

p

∆ t (14)

explicitando o termo T 1 ,np

T 1 ,np+1=2∆ t

∆ r [ hóleo

ρ c p

(T óleo−T 1 ,np+1 )+ α

∆ r(T 2 ,n

p+1−T1. np+1 )]+T 1 ,n

p

(15)

na qual hóleo é o coeficiente de transferência de

calor por convecção do óleo no interior do

absorvedor e T óleo a temperatura do óleo.

Figura 03: Nós das superfícies interna e externa do cilindro (a, b) sujeitos à convecção e à condução transiente e (c) sujeitos a um fluxo conhecido e à condução transiente.

3.3.2. ANÁLISE DO VOLUME DE CONTROLE EXTERNO

Na superfície externa foi considerado uma aproximação razoável dizer que dois terços sofre ação principal da convecção natural, e um terço tem ação principal dos raios do sol que passam pelo concentrador, como ilustrado na figura abaixo.

Figura 04: (a)Malha estrutura e (b) equações utilizadas no domínio discretizado. Em destaque, (a) as posições utilizadas na comparação dos resultados quando se atinge a condição de equilíbrio.

Para a parte sob convecção natural, Figura 03-b, utilizando a equação 12, por ter mesmas considerações, e o método diferenças finitas implícito tem-se

har (Tar−T1 , np+1 )+ K

∆ r(T 2, n

p+1−T 1 ,np+1 )=ρ c p

∆ r2

T 1 ,np+1−T1 , n

p

∆ t

(16)

Explicitando o termo T mm ,np ,

T mm,np+1 −2 ∆ t

∆ r [ har

ρ c p

(T inf−T mm,np+1 )+ α

∆ r(T mm−1 ,n

p+1 −T mm.np+1 ) ]=T mm,n

p

(17)

na qual har é o coeficiente de transferência de

calor por convecção do ar e T ar a temperatura

do ar.

Para a outra parte, procedendo de forma análoga.

q} + {K} over {∆r} left ({T} rsub {mm-1,n} rsup {P+1} - {T} rsub {mm,n} rsup {P+1} right ) =ρ {c} rsub {p} {∆r} over {2} ( {{T} rsub {mm,n} rsup {p+1} - {T} rsub {mm,n} rsup {p}} over {∆t} ¿

(18)

onde k representa a condutividade térmica do vidro e q” o fluxo de calor oriundo do concentrador.

3.3.3. CONDIÇÕES INICIAIS

Para resolução numérica do problema é necessário para determinação de condições iniciais. Assume-se que inicialmente a camada interna do tubo esteja em equilíbrio com ou óleo, ou seja, temperatura interna do tubo seja igual à temperatura do óleo.

T 1 ,n1 =T óleo

Para a parte externa é considerado a mesma forma, temperatura da camada externa com mesmo valor da temperatura do ar em todo o tubo.

T mm,n1 =T ar

Para os pontos internos, com m variando de (2, mm-1), podemos considerar inicialmente, por ser menos trabalhoso, que a temperatura assume uma média entre a temperatura do óleo e a temperatura do ar.

T m,n1 =

T óleo+Tar

2

Uma Observação importante é que, por se tratar de um tubo, os pontos de Ф=0 e Ф=2PI os pontos com são coincidentes assumindo a mesma temperatura:

T m,1p =T m, nn

p

4. DADOS

No absorvedor foi considerada somente a transferência de calor entre o vidro e óleo, desconsiderando a existência do tubo de metal entre o vidro e óleo.

Para fins de comparação, foi utilizado o mesmo modelo e valores do trabalho de conclusão de curso da aluna Letícia Jenisch Rodrigues que realiza de forma similar a análise do tubo coletor, utilizando, porém a formulação explícita.

A figura 04 abaixo ilustra os dados referentes às condições de contorno iniciais do problema a ser analisado, além das dimensões do modelo físico utilizado. Enquanto que na tabela 01 apresenta as propriedades físicas do vidro.

Figura 04- a) Condições de contorno e b) Dimensões do modelo utilizado

Propriedades do Vidro ValorDensidade (ρ) 145 kg/m³Calor específico à pressão constante (c P)

0,8 W/mk

Condutividade térmica (k) 1000gKTabela 01: Propriedades físicas do Vidro

5. ANÁLISE DE RESULTADOS

Na figura 05 abaixo, é apresentado o campo de temperaturas, utilizando o método diferenças finitas explícitos, encontrado no trabalho base.

Figura 05: Campos de temperaturas obtidos pelos softwares (a) TRANSCAL e (b) algoritmo CSP utilizando 10 pontos na direção radial e 41 pontos na direção angular.

Os resultados obtidos a partir da aplicação do método diferenças finitas de forma implícita podem ser observados na figura 06 abaixo

Figura 06: Campo de temperatura encontrado utilizando o código implementado.

Infelizmente é possível perceber que há uma discrepância entre os dois resultados encontrados, a mesma não deve ser associada à diferença entre os métodos aplicados.

Sendo portando essa discrepância associada aos problemas na implementação do código computacional.

6. CONCLUSÃOA proposta inicial de análise da aplicação do método das diferenças finitas na solução de problemas de difusão de calor em um tubo coletor foi realizada com sucesso como pode ser visto ao longo deste trabalho, entretanto essa análise não pode ser comprovada por meio da verificação dos dados fornecidos pelo código computacional implementado.

De acordo com os dados apresentados é possível observar que o objetivo proposto inicialmente de implementar um código computacional que solucionasse o problema de transferência de calor em um absorvedor cilíndrico não foi alcançado.

Levando em consideração que a análise teórica do problema está correta, é possível perceber que os resultados encontrados não condizem com o esperado, pois o código implementado apresentou alguns erros.

Vale ressaltar que a não foi possível verificar que a principal diferença que pode ser associada quanto ao uso dos métodos explicito e implícito, o método explicito apresenta limitações temporais que não são vistas no método implícito e o método explicito oferece vantagens computacionais sobre o implícito, já que os resultados encontrados não estão corretos.

REFERÊNCIAS

1. CENGEL, Y. A. Heat and Mass Transfer: A Practical Approach. 3ª edição New York, United States: McGraw-Hill, 2007.

2. INCROPERA, Frank P., DEWITT, David P., BERGMAN, Theodore L., LAVINE, Adrienne S. Fundamentos de Transferência de Calor e de Massa. 6ª Edição Rio de Janeiro, Brasil, LTC, 2008.

3. RODRIGUES, Letícia J. Análise Transiente da Transferência de Calor em um Tubo através do Método das Diferenças Finitas, Monografia (Graduação em Engenharia Mecânica), Universidade Federal do Rio Grande do Sul, Porto Alegre, 2011.

4. SALVADORETTI, José Luiz. Modelo M4atemático para Análise do Desempenho Térmico de Coletores Solares Cilindro. Dissertação (Mestrado em Engenharia da Energia, Metalurgia e Materias) – Programa de Pós- Graduação em Engenharia da Energia, Metalurgia e Materias, Universidade Federal do Rio Grande do Sul, Porto Alegre, 1983.

ANEXO

Código Implementado no Matlab:

clear allclc %**************************************************************************% ENTRADA DE DADO DO USUARIO%**************************************************************************M=input('Entre com o numero de pontos(M) no raio: ');N=input('Entre com o numero de pontos(N) na circunferencia: ');delta_t=0.001 %input('Entre com o valor do passo no tempo em segundos: ');Tempo_total=input('Entre com o valor total do tempo de simulação em segundos: ');erro=0.00000001kmax=50000 %**************************************************************************% DEFINICAO DAS COSNTANTES%************************************************************************** di=0.050; % m Diametro interno em metrosde=0.080; % m Diametro externo em metrosri=di/2; % m Raio interno em metrosre=de/2; % m Raio externo em metrosdelta_r=(re-ri)/(M-1); % m Variaçao padrao de raio entre pontosdelta_phi=(2*pi)/N; % m Variaçao padrao de raio entre pontosCp=1000; % J/kg.K (calor específico a pressao constante do vidro)K=0.8; % W/m.k (condutividade termica do vidro)ro=145; % k/m³ (densidade do vidro)alpha=K/(ro*Cp); % m2/s (difusividade termica do cobre)T_ar=30+273.15; % k Temperatura do ar ambienteh_ar=50; % W/m2 (coeficiente de convecçao do ar)

T_oleo=300+273.15; % k Temperatura do oleoh_oleo=500; % W/m2 (coeficiente de convecçao do oleo)q_rad=16000; % W/m2 potencia de radiaçao incidente %Valor de N mais proximo de 240 graus de phiN_2_3=round(2*N/3); %Matriz de temperatura inicial MxNcond_ini(1:M,1:N)=(T_oleo+T_ar)/2;cond_ini(1,1:N)=T_ar;cond_ini(M,1:N)=T_oleo; %Cria figura para frame do videoh.f=figure('position',[100 300 800 600]);set(h.f,'DoubleBuffer','on'); %Arquivo de videomov = avifile('TransfCalor.avi','fps',30,'compression', 'None'); %arquivo de video do output %Decide a resoluçao do graficoif (N<5) res=N*10;elseif(N>=5 & N<8) res=N*7;elseif(N>=8 & N<15) res=N*5;elseif(N>=15 & N<20) res=N*3;else res=N*2;end%define o grid do graficor=linspace(ri,re,res);theta=linspace(0,2*pi,res);[theta,r]=meshgrid(theta,r); grid_x=r.*cos(theta);grid_y=r.*sin(theta); %Ultimo passo da simulaçãot_final=round(Tempo_total/delta_t); %**************************************************************************% INICIALIZAÇÃO DE VARIÁVEIS%**************************************************************************rm=re; %Raio no ponto de cordenada radial mm=1; %Coordenada radialn=1; %Coordenada angularp=0; %Variavel tempo % Matriz de temperaturas no espaço da seçao MxN em t_final instantesT_mn(1:M,1:N,t_final)=0; %Variaveis do sistema linearA(1:M*N,1:M*N)=0; %Matriz dos coeficientes do sistema linearC(1:M*N)=0; %Vetor dos termos independentes

%**************************************************************************% MONTA A MATRIZ (A) DE COEFICIENTES DO SISTEMA LINEAR PARA TEMPO p% E VETOR (C) DE TERMOS INDEPENDENTES%**************************************************************************a_conv_ar1 = (K/delta_r) + ((ro*Cp*delta_r)/(2*delta_t)) + h_ar; %componente a(n,n) da matriz Aa_conv_ar2 = -(K/delta_r); %componente a(n,N+n) da matriz Aa_rad1 = (K/delta_r) + ((ro*Cp*delta_r)/(2*delta_t)); %componente a(n,n) da matriz Aa_rad2 = -(K/delta_r); %componente a(n,N+n) da matriz A a_conv_oleo1= (K/delta_r) + ((ro*Cp*delta_r)/(2*delta_t)) + h_oleo; %componente a[(M-1)N+n,a(M-1)N+n] da matriz Aa_conv_oleo2= -(K/delta_r); %componente a[(M-1)N+n,a(M-2)N+n] da matriz A %c_conv_ar= ((ro*Cp*delta_r)/(2*delta_t));<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<COMO ISSO VEIO PARAR AQUI?< %monta a matriz A e vetro C com os coeficientes calculados anteriormentefor m=1:M for n=1:N if (m==1 & n<=N_2_3) %Se for ponto no contorno externo nos primeiros 2 terços de phi (convecçao ar) A(n,n)=a_conv_ar1; A(n,N+n)=a_conv_ar2; end if (m==1 & n>N_2_3) %Se for ponto no contorno externo no terceiro terço de phi (radiaçao) A(n,n)=a_rad1; A(n,N+n)=a_rad2; end%-----------TERMO DEPENDENTE DE DA COORDENADA (m,n)---------------------------------------- if (m>1 & m<M) %Se for ponto interior dos contornos (conduçao) rm=re-(m*delta_r); %Raio no ponto de cordenada radial m a_cond1=1-((alpha*delta_t)/(delta_r^2))-(2*alpha*delta_t)/((rm^2)*(delta_phi^2)); %componente a[(m-1)N+n,(m-1)N+n] da matriz A a_cond2=(-(alpha*delta_t))*((1/(rm*2*delta_r))+(1/(delta_r^2))); %componente a[(m-1)N+n,(m-2)N+n]=a[(m-1)N+n,mN+n] da matriz A a_cond3=-(alpha*delta_t)/((rm^2)*(delta_phi^2)); %componente a[(m-1)N+n,(m-1)N+n-1]=a[(m-1)N+n,(m-1)N+n+1] da matriz A A((m-1)*N+n,(m-1)*N+n)=a_cond1; A((m-1)*N+n,(m-2)*N+n)=a_cond2; A((m-1)*N+n,m*N+n)=a_cond2; if(n==1) %Se n=1 o ponto N-1 é igual ao phi maximo A((m-1)*N+n,(m-1)*N+N)=a_cond3; A((m-1)*N+n,(m-1)*N+n+1)=a_cond3; end if(n>1 & n<N) %Se n nao é nem o primeiro nem o ultimo A((m-1)*N+n,(m-1)*N+n-1)=a_cond3;

A((m-1)*N+n,(m-1)*N+n+1)=a_cond3; end if(n==N) %Se n=N o ponto N+1 é igual ao 1 A((m-1)*N+n,(m-1)*N+n-1)=a_cond3; A((m-1)*N+n,(m-1)*N+1)=a_cond3; end end%------------------------------------------------------------------------- if (m==M) %se for ponto no contorno interno (convecçao oleo) A((M-1)*N+n,(M-1)*N+n)=a_conv_oleo1; A((M-1)*N+n,(M-2)*N+n)=a_conv_oleo1; end endend %*************************************************************************% MONTA O VETOR (C) DE TERMOS INDEPENDENTES%*************************************************************************%vetor C inicial for m=1:M for n=1:N if (m==1 & n<=N_2_3) %Se for ponto no contorno externo nos primeiros 2 terços de phi (convecçao ar) C((m-1)*N+n)=((ro*Cp*delta_r)/(2*delta_t))*cond_ini(m,n)+h_ar*T_ar; end if (m==1 & n>N_2_3) %Se for ponto no contorno externo no terceiro terço de phi (radiaçao) C((m-1)*N+n)=((ro*Cp*delta_r)/(2*delta_t))*cond_ini(m,n)+q_rad; end if (m>1 & m<M) %Se for ponto interior dos contornos (conduçao) C((m-1)*N+n)=cond_ini(m,n); end if (m==M) %se for ponto no contorno interno (convecçao oleo) C((m-1)*N+n)=((ro*Cp*delta_r)/(2*delta_t))*cond_ini(m,n)+h_oleo*T_oleo; end end end %--------------------------------------------------------------------------% INICIA O LOOP TEMPORAL%--------------------------------------------------------------------------while(p<t_final)if(p~=0) for m=1:M for n=1:N if (m==1 & n<=N_2_3) %Se for ponto no contorno externo nos primeiros 2 terços de phi (convecçao ar) C((m-1)*N+n)=((ro*Cp*delta_r)/(2*delta_t))*T_mn(m,n,p)+h_ar*T_ar; end

if (m==1 & n>N_2_3) %Se for ponto no contorno externo no terceiro terço de phi (radiaçao) C((m-1)*N+n)=((ro*Cp*delta_r)/(2*delta_t))*T_mn(m,n,p)+q_rad; end if (m>1 & m<M) %Se for ponto interior dos contornos (conduçao) C((m-1)*N+n)=T_mn(m,n,p); end if (m==M) %se for ponto no contorno interno (convecçao oleo) C((m-1)*N+n)=((ro*Cp*delta_r)/(2*delta_t))*T_mn(m,n,p)+h_oleo*T_oleo; end end endend %**************************************************************************% RESOLVE O SISTEMA LINEAR POR GAUSS SEIDEL%**************************************************************************p=p+1;%Vetor soluçao (temperaturas)do sistema linarT_k(1:(M*N))=C(1:M*N);T_k_1(1:(M*N))=C(1:M*N);e=erro+1;k=0; while(e>erro & k<kmax) for i=1:1:(M*N) soma=0; for j=1:1:(M*N) if(i>j) soma=(A(i,j)*T_k(j))+soma; end if(i<j) soma=(A(i,j)*T_k_1(j))+soma; end end T_k(i)=((C(i)-soma)/A(i,i)); end e=max(abs((T_k)-(T_k_1)))/max(abs(T_k)); T_k_1(1:(M*N))=T_k(1:(M*N)); k=k+1;end for m=1:1:M for n=1:1:N %if(T_k((m-1)*N+n)<cond_ini(m,n)) % T_mn(m,n,p)=cond_ini(m,n); %else T_mn(m,n,p)=T_k((m-1)*N+n)-273.15; %end endend %**************************************************************************

% PLOTA OS DADOS NUM GRAFICO E ESCREVE NUM VIDEO%**************************************************************************x=1;y=1; t=p;%qual o tempo plotar o grafico delta_x=(re-ri)/res;delta_y=(2*pi)/res; dx=0;%distancia de varredura efetuada no raiody=0;%angulo de varredura efetuada na circunferencia m_x=delta_x/delta_r;n_y=delta_y/delta_phi; z(1:res,1:res)=0; for m=1:1:M-1 dm=m*delta_r; for n=1:1:N dphi=n*delta_phi; while (dx<dm) while (dy<((res)*delta_y)) x1=(m-1)*delta_r; x2=m*delta_r; Xx=((dx-x1)/(x2-x1)); y1=(n-1)*delta_phi; y2=n*delta_phi; Yy=((dy-y1)/(y2-y1)); deltaTxy=Xx*((T_mn(m+1,n+1,t)-T_mn(m,n+1,t))-(T_mn(m+1,n,t)-T_mn(m,n,t)))+T_mn(m,n+1,t)-T_mn(m,n,t); z(res-x+1,y)=Yy*deltaTxy + Xx*((T_mn(m+1,n,t)-T_mn(m,n,t)))+T_mn(m,n,t); y=y+1; dy=(y-1)*delta_y; end while (n==10000 & dy<=(res*delta_y) & dy>((res-1)*delta_y)) x1=(m-1)*delta_r; x2=m*delta_r; Xx=((dx-x1)/(x2-x1)); y1=(n-1)*delta_phi; y2=n*delta_phi; Yy=((dy-y1)/(y2-y1)); deltaTxy=Xx*((T_mn(m+1,1,t)-T_mn(m,1,t))-(T_mn(m+1,n,t)-T_mn(m,n,t)))+T_mn(m,1,t)-T_mn(m,n,t); z(res-x+1,y)=Yy*deltaTxy + Xx*((T_mn(m+1,n,t)-T_mn(m,n,t)))+T_mn(m,n,t);

y=y+1; dy=(y-1)*delta_y; end y=1; dy=0; x=x+1; dx=(x-1)*delta_x; end end end surf(grid_x,grid_y,z);shading interpgrid on%box onaxis tightxlabel('grid_x-axis')ylabel('grid_y-axis')zlabel('z-axis') lighting phongset(gcf,'Renderer','zbuffer') drawnow; F=getframe(h.f);mov=addframe(mov,F); end%--------------------------------------------------------------------------% FIM DO LOOP TEMPORAL%-------------------------------------------------------------------------- mov=close(mov);