Condução de calor nao-linear

Condução de calor nao-linear

(Parte 3 de 3)

i i i i r r dr d r drd r

1 sobre
10 em011

ou, simplesmente por

1 sobre
10 em011

r dr d r drd r i i i i αα

Como o problema só tem sentido físico para T > 0, temos que kmin > 2. Desta forma, podemos trabalhar com qualquer α ≥ 1/2. Usemos α = 3!

A solução geral da equação

10 para011122<≤=+

d r drdr

10 para612
i(4.46)

onde a constante Ci+1, para i > 0, é obtida das condições de contorno. Em outras palavras

+−=+iC(4.47)

Então,

1(4.48)

A constante C0 é obtida a partir de

+−=C(4.49)

Se usássemos α = 20, teríamos

1(4.50)

201 e a constante C0 seria obtida de

+−=C(4.51)

A tabela 4.1 apresenta uma comparação entre os valores obtidos para Ci com os dois valores distintos de α.

Tabela 4.1 – comparação entre resultados obtidos com α = 3 e α = 20.

É fácil mostrar que C∞ = 1. Portanto, o limite da seqüência é dado por

10 para1

A solução é exatamente o limite da seqüência. Levando em conta que

++−==−ωωTf(4.53)

temos o seguinte resultado (coincidente com a solução exata, previamente obtida)

+−++−=rT(4.54)

4.5.2 – Exemplo 2 Pensemos em outro caso, representado por

∞ sobre
em01

TTh dr dTk q dr dTrk drd r onde q&, ∞T e h são constantes positivas conhecidas e k é uma função da temperatura dada por

k+=(4.56)

sendo k0 uma constante positiva e conhecida. Este modelo descreve um problema de condução de calor num corpo cilíndrico de raio R, que ocupa uma região Ω com fronteira ∂Ω, e cujas faces circulares são isoladas termicamente.

Apliquemos então o procedimento proposto. A transformada de Kirchoff para este problema, bem como sua inversa, são representadas por kTk T Tk

fTdk T

ω(4.57)

ωψψ A transformada nos permite reescrever o problema na forma

sobre)(ˆ
em01

1 Tfh dr q dr d r drd r

soluções são Φ1, Φ2, Φ3,, Φ∞ e onde cada elemento da seqüência é dado por

Tomemos a solução como sendo o limite da seqüência de problemas lineares cujas

i i

Tfh dr d q dr d r drd r αβ βα )(ˆ sobre

em01

onde kTk T

Tkf

Resolvendo a equação que carrega o termo difusivo temos que

com 2

A constante Ci+1 pode ser determinada a partir da condição de contorno (r = R):

i i

Rq dr

Reescrevendo Ci+1 apenas em função de Φi temos

2 1 RqTfhRqC i

Lembrando que )(1ifΦ− é calculada a partir de (4.60). Finalmente, a solução geral de (4.58) será

A partir da transformada – equação (4.57) – podemos determinar a solução de (4.5), utilizando a solução encontrada para ∞Φ:

kTk T Tk

Utilizemos os seguintes valores para os parâmetros físicos e geométricos para encontrar o valor da constante C:

Levando em conta a condição de convergência da seqüência 10

min h α temos que 3≥α. Analisemos três valores para α: 1, 3 e 200.

Tabela 4.2 – comparação entre resultados obtidos com α = 1, α = 3 e α = 200.

Resolvendo o sistema (4.5) por métodos clássicos chegamos à solução:

krqMMr q

onde a constante de integração M é calculada a partir de

Rqk hTRq hRqM

que, inserindo os valores das constantes do problema é igual a 5646,6635 WK/m. Para que as duas soluções sejam iguais é necessário que

TkCM−+=(4.68)

Para dois valores de α – 3 e 200 – a constante C convergiu para (5637,6635), que atende à condição estabelecida em (4.68), com a diferença que, para α = 3, C atingiu o limite muito mais rapidamente do que para α = 200. Note ainda que para α = 1, que viola a restrição 3≥α, a série simplesmente não converge. A conclusão final é que ambos os métodos resultam na mesma solução:

T(4.69)

Como foi dito na seção 4.3, a seqüência de funções Φi, que representa a transformada de Kirchoff para cada iteração, é não-decrescente, desde que seja escolhido um α adequado. Os dois exemplos anteriores já confirmam esta propriedade. Ainda assim, consideremos um outro problema, no qual será possível observar graficamente o crescimento da seqüência e a convergência das soluções obtidas em cada passo do processo iterativo para a solução final. O problema é matematicamente representado pelo sistema

()0 em
em
0 em0

∞ ∞ xTThx

LxTThx Tk

Lxqx Tkx onde a condutividade térmica k é dada por

1000+=Tk(4.71)

- 39 - A solução analítica deste problema é representada implicitamente pela equação onde 0T é a temperatura em 0=x e é dada por

= Th qL

Th qLhL

Th qLTh qLhLhLTLq

Usando os seguintes valores para as constantes:

obtemos o campo de temperatura mostrado na próxima figura:

Figura 4.1 – curva de temperatura representando solução do sistema (4.70).

Empreguemos o procedimento iterativo proposto para o mesmo problema.

Plotando as soluções obtidas a partir da centésima iteração em diante (plotar os resultados desde o primeiro elemento da seqüência prejudicaria a visualização) podemos ver claramente o processo de convergência para a solução mostrada pela figura anterior, obtida por métodos clássicos de resolução de sistema de EDP’s – o campo de temperatura do centésimo elemento da seqüência é o mais inferior da figura.

Figura 4.2 – seqüência de soluções obtidas pelo método iterativo (curvas em azul) comparada com a solução dada por (4.72) e (4.73) (curva em vermelho).

Figura 4.3 – comparação do resultado final obtido pelo procedimento iterativo com a solução dada por (4.72) e (4.73).

É fácil notar o comportamento crescente da seqüência de problemas lineares e o processo de convergência para a solução exata do sistema (4.70).

5 – MÉTODOS NUMÉRICOS EM PROBLEMAS DE TRANSMISSÃO DE CALOR

Considere o seguinte problema:

∞ sobre)(
em0)(

4 GTTThnTkgrad qTkgraddiv

Este modelo descreve a distribuição de temperaturas de um certo corpo que ocupa uma região Ω. Para que conheçamos o valor da temperatura em cada ponto do domínio Ω é preciso resolver a EDP sujeita às condições de contorno definidas sobre ∂Ω. Os métodos analíticos de solução, quando praticáveis, provêem solução exata do sistema (5.1). No entanto, o número de problemas de engenharia que podem ser resolvidos analiticamente é extremamente limitado, seja devido à presença de termos não-lineares ou em conseqüência de geometrias complexas. Os métodos numéricos são uma alternativa para muitos dos casos não alcançados pelo ferramental analítico conhecido, porém fornecem sempre soluções que carregam algum erro associado. Nos procedimentos numéricos o domínio, antes contínuo, é discretizado, bem como a EDP definida nele, que pode ainda ser discretizada no tempo (o que não é necessário para problema (5.1), que não possui dependência do tempo). Os métodos de análise numérica, ao tratar a solução em um número finito de pontos discretos, simplificam a resolução do sistema de equações diferenciais transformando-o num sistema de equações algébricas. Dentre os métodos numéricos clássicos de resolução de EDP’s estão o Método das Diferenças Finitas (MDF), o Método de Volumes Finitos (MVF) e o Método de e Elementos Finitos (MEF). Neste trabalho serão empregados o MDF e o MEF.

5.1 – Diferenças Finitas

O MDF consiste em reescrever equações diferenciais de forma que as derivadas sejam tomadas em intervalos finitos, ao invés de infinitesimais. Ou seja, a derivada da função )(xff=, definida por

lim(5.2)

= →∆ 0 seria, em diferenças finitas, escrita como segue:

dx df ∆

(5.3)

onde ∆x é uma variação finita (positiva). Vejamos uma aplicação do MDF para o seguinte problema de transmissão de calor

∞ sobre)(
em0)(

TThnTkgrad k q Tgraddiv

Se as equações acima estiverem em um sistema de coordenadas cartesianas retangulares e o domínio Ω for uma placa retangular de largura L e altura H, o sistema (5.4) pode ser reescrito na seguinte forma:

( ) HyTThy Tk yTThy Tk

LxTThx Tk xTThx Tk

HyLx kqyTx

para
0 para
para
0 para
0 e 0 em0222

Discretizemos o domínio Ω. Devemos estipular um número de nós nc para a direção x e nl para a direção y. Desta forma criamos uma malha de diferenças finitas, ilustrada na figura 5.1. Digamos que nc = 7 e nl = 6.

Figura 5.1 – malha de diferenças finitas Os valores de ∆x e ∆y são calculados a partir de

nl Hy

L x

A partir do domínio discreto o sistema de equações diferenciais (5.5) pode ser aproximado para o seguinte sistema de equações lineares algébricas:

( ) ncnTTh y

T k ncnTTh y

T k nlmTTh x

T k nlmTTh x

T k ncnnlm kqy TTTx nnl nnlnnl n n ncm ncmncm m m nmnmnmnmnmnm

1 para
1 para
1 para
1 para
1 e 1 em0

Em Kreith e Bohn [7] há mais detalhes de como se chega às equações acima. A cada nó da malha há uma equação associada. No caso de placas retangulares, os nós das

“quinas” – T1,1, T1,nc, Tnl,1, Tnl,nc – não fazem parte do sistema linear, podendo ser calculados pela média aritmética dos dois nós “vizinhos”. Desta forma, para este exemplo, a solução aproximada do problema (5.5) seria obtida da solução de um sistema de 38 equações lineares (generalizando, para placas retangulares, são 4−×ncnl equações). Quanto à resolução do sistema linear há diversas metodologias, sejam aproximadas (e.g. métodos iterativos de Gauss-Seidel ou Jacobi) ou exatas (por meio de decomposição LU, inversão direta da matriz dos coeficientes, entre outros caminhos).

Vejamos uma aplicação do MDF para o problema 2-D representado por (5.5), com os seguintes valores para as propriedades físicas:

A placa tem 15m de largura e 10m de altura. É utilizada uma malha 4040×, e para resolver o sistema de equações lineares empregamos o método iterativo Gauss-Seidel. O resultado é apresentado na figura 5.2.

Figura 5.2 – campo de temperaturas em placa retangular. Todas propriedades físicas são constantes. Como esperado, a distribuição de temperaturas apresenta várias simetrias.

Consideremos o mesmo problema do exemplo anterior porém com geração interna de energia dependente da posição – ),(yxqq=&. Utilizemos

O campo de temperaturas calculado é graficamente representado a seguir.

Figura 5.3 – campo de temperaturas em placa com geração de energia dependente da posição. 5.1.3 – Exemplo 3

Ainda fazendo pequenas mudanças no mesmo problema, consideremos o coeficiente de troca de calor por convecção h variável com a posição, mantendo a relação para q& dada pela equação (5.8). Utilizemos as seguintes funções para h

(absolutamente arbitrárias e sem comprometimento físico):

( ) Hyhxhh yhL xhh

Lxhyhh xhyhh para K, W/m1 , 1

0 para K, W/m1 , sen para K, W/m1 , 0 para K, W/m1 ,

pi(5.9)

Resultado:

Figura 5.4 – campo de temperaturas em placa com geração de energia e coeficiente de filme dependentes da posição.

5.2 – Elementos Finitos

Formulemos a solução em elementos finitos para a equação

em0qy

yx Tkx sujeita às seguintes condições de contorno:

−∞ sobreTThnjy

∂ Tix

Este é um problema bidimensional de condução de calor num corpo sujeito à lei de Newton do resfriamento. Os parâmetros físicos k (condutividade térmica), q& (geração interna de calor) e h (coeficiente de troca de calor por convecção) serão considerados constantes.

5.2.1 – Em busca da forma fraca das equações governantes

A forma fraca das equações que governam o problema consiste no estabelecimento de equações integrais sobre o domínio Ω e o contorno Γ do corpo, referentes à satisfação destas equações em um sentido “médio”, em oposição ao sentido restrito pontual da forma forte da equação (5.10).

Obs.: Nesta seção utilizaremos letras em negrito, maiúsculas e minúsculas, para indicar matrizes.

Multiplicando a EDP definida em Ω por uma função arbitrária ),(yxw, denominada função teste ou função peso, obtemos a equação yx Tkx

a qual, se integrada sobre o domínio Ω, resulta em

dqy yx Tkx

Note que, caso seja determinado o campo de temperatura para o qual a equação (5.10) seja satisfeita, então a equação (5.13) também é satisfeita para qualquer que seja a função ),(yxw. Por outro lado, caso se determine uma solução ),(yxT que satisfaça a equação (5.13) para qualquer função ),(yxw, então ),(yxT também é solução de

(5.10). Utilizando o conceito de derivada do produto, é fácil ver que

dqw

Tkywx Tkxwy

Twk yx

Twk x dqy yx Tkx

Agora, pelo teorema de Gauss ou da divergência, temos que

∂ dnjy

Tkix Tkwdy

Twk yx Twkx

Fazendo

qjx Tkix

Tk r =

(5.16)

∂ a equação (5.13) pode ser reescrita como dnqwdqw

Tkywx Tkx w r

O fluxo de calor nqrr ⋅ que cruza a fronteira Γ é igual a ()∞−TTh, o qual, inserido na segunda integral de (5.17), resulta em

e a equação (5.17) passa a ser dada por dhThTwdqw

Tkywx Tkx

5.2.2 – Discretização do domínio em elementos finitos

(Parte 3 de 3)

Comentários