Método das Diferenças Finitas

O trabalho da disciplina Transferência de Calor e Massa será desenvolvido através da resolução de um problema de condução bidimensional de calor em regime permanente utilizando o método de diferenças finitas A Figura 1 mostra uma representação do problema que deve ser resolvido e as demais informações são Placa quadrada de 1m x 1m com condutividade térmica k 353 WmK Temperatura uniforme de 300C no canto superior esquerdo e de 150C no canto inferior direito O canto inferior esquerdo e o superior direito encontramse isolados termicamente A transferência de calor ocorre em regime permanente Resolver o problema utilizando malhas 8x8 20x20 30x30 e 40x40 considerando espaçamentos iguais na horizontal e vertical Δx Δy Comparar o fluxo de calor por unidade de comprimento q em Wm obtido com a resolução numérica com o resultado analítico que é dado por q T1 T2 k O trabalho deve ser desenvolvido em grupos de no máximo 4 pessoas Pode ser utilizada qualquer linguagem de programação sendo que o arquivo com o código desenvolvido deve ser enviado junto com o relatório este último deve conter apenas os resultados das análises e conclusões para o email Figura 1 Placa quadrada com transferência de calor em regime permanente Considere uma placa em processo de transferência de calor em regime permanente Figura 1 Figura 1 Esquema ilustrativo da placa em transferência de calor Desejase obter o perfil de temperatura da condução de calor bidimensional utilizando o método das diferenças finitas e comparar o resultado com a solução analítica Dados do problema 1 Placa quadrada com 1m x 1m e condutividade térmica igual a 353 WmK 2 Temperatura uniforme de 300C no canto superior esquerdo e de 150C no canto inferior direito 3 Os cantos inferior esquerdo e superior direito encontramse isolados termicamente 4 O processo ocorre em regime permanente 5 Resolver o problema utilizando malhas 8x8 20x20 30x30 e 40x40 considerando espaçamentos iguais na horizontal e na vertical Δx Δy 6 Comparar o fluxo de calor por unidade de comprimento q em Wm obtido com a resolução numérica com o resultado analítico que é dado por q T1 T2k Hipóteses do problema 1 Regime permanente 2 Transferência de calor bidimensional em x e y x y 3 Propriedades constantes e isotrópicas 4 Sistema sem geração térmica Modelagem do problema O perfil de temperatura na placa é descrito pela equação da difusão do calor Equação 1 Aplicandose as hipóteses simplificadores obtémse a EDP da Equação 2 𝑥 𝑘 𝑇 𝑥 𝑦 𝑘 𝑇 𝑦 𝑧 𝑘 𝑇 𝑧 𝑞𝐺 𝑡 𝜌𝑐𝑝𝑇 1 𝑥 𝑘 𝑇 𝑥 𝑦 𝑘 𝑇 𝑦 𝑧 𝑘 𝑇 𝑧 𝑞𝐺 𝑡 𝜌𝑐𝑝𝑇 ²𝑇 𝑥² ²𝑇 𝑦² 0 2 A EDP da Equação 2 necessita de quatro condições de contorno para ser resolvida Neste problema adicionar as condições de contorno comuns x 0 q 0 x L T T2 y 0 q 0 y L T T1 não vai gerar um resultado correto pois a condição não se mantém a mesma para toda a abscissa ou toda a ordenada Logo convém aplicar o método numérico das diferenças finitas Primeiro devese definir a representação de um nó o menor elemento de aplicação do método Figura 2 Existem vários tipos diferentes de nó Para este problema serão necessárias as equações de diferenças finitas para três tipos ponto nodal interior Figura 2 ponto nodal em uma superfície plana com fluxo uniforme e ponto nodal em um vértice externo fluxo uniforme Isso ocorre em função da condição de contorno a qual o nó é submetido 1 Ponto nodal interno 𝑇𝑚𝑛 𝑇𝑚1𝑛 𝑇𝑚𝑛1 𝑇𝑚𝑛1 𝑇𝑚1𝑛 Figura 2 Representação de um nó 𝑦 𝑥 Sobre este nó as derivadas parciais podem ser aproximadas pelo método das diferenças centradas Equação 3 ²𝑇 𝑥² 𝑇𝑚1𝑛 2𝑇𝑚𝑛 𝑇𝑚1𝑛 𝑥² ²𝑇 𝑦² 𝑇𝑚𝑛1 2𝑇𝑚𝑛 𝑇𝑚𝑛1 𝑦² 3 Substituindo as derivadas parciais da Equação 3 na EDO da Equação 2 obtémse ²𝑇 𝑥² ²𝑇 𝑦² 0 𝑇𝑚1𝑛 2𝑇𝑚𝑛 𝑇𝑚1𝑛 𝑥² 𝑇𝑚𝑛1 2𝑇𝑚𝑛 𝑇𝑚𝑛1 𝑦² 0 Considerando que as dimensões da malha são iguais ou seja que Δx Δy h obtém se a Equação 4 do ponto nodal interno 𝑇𝑚1𝑛 𝑇𝑚1𝑛 𝑇𝑚𝑛1 𝑇𝑚𝑛1 4𝑇𝑚𝑛 0 4 2 Ponto nodal em uma superfície plana com fluxo térmico uniforme Neste caso o calor condutivo que chega dos nós adjascentes deve ser igual ao fluxo de calor uniforme que sai da peça Caso o fluxo de calor esteja entrando na peça conforme ilustrado pela Figura 6 o calor condutivo vai entrar na peça até as peças adjascentes Realizando o balanço de energia obtémse a Equação 5 𝑘𝑦 𝑇𝑚1𝑛 𝑇𝑚𝑛 𝑥 𝑘 𝑥 2 𝑇𝑚𝑛1 𝑇𝑚𝑛 𝑦 𝑘 𝑥 2 𝑇𝑚𝑛1 𝑇𝑚𝑛 𝑦 𝑞𝑦 0 𝑇𝑚𝑛 𝑞 𝑇𝑚𝑛1 𝑇𝑚𝑛1 𝑇𝑚1𝑛 𝑦 𝑥 Figura 3 Representação de um ponto nodal em uma superfície plana com fluxo térmico uniforme 2𝑇𝑚1𝑛 𝑇𝑚𝑛1 𝑇𝑚𝑛1 2𝑞𝑥 𝑘 4𝑇𝑚𝑛 0 5 3 Ponto nodal em vértice externo com fluxo uniforme Neste caso os três nós expostos estão submetidos ao fluxo constante Porém os nós adjascentes enviamrecebem do canto um fluxo condutivo Portanto o fluxo de calor por condução nos nós adjascentes internos deve ser igual ao fluxo de calor por convecção nas superfícies expostas Desenvolvendo o balanço de energia e considerando os intervalos iguais obtémse a Equação 6 𝑘 𝑦 2 𝑇𝑚𝑛 𝑇𝑚1𝑛 𝑥 𝑘 𝑥 2 𝑇𝑚𝑛 𝑇𝑚𝑛1 𝑦 𝑞𝑥 𝑦 0 𝑇𝑚1𝑛 𝑇𝑚𝑛1 2𝑞 𝑘 𝑥 𝑦 2𝑇𝑚𝑛 0 6 Resolução do problema De acordo com a Figura 1 o problema vai apresentar nós dos três tipo Quanto maior a quantidade de nós maior vai ser a precisão do perfil de temperatura pois o Δx e o Δy serão iguais Independente do tamanho da malha haverá apenas dois pontos pertencentes ao caso 3 ponto nodal em vértice com fluxo constante Os pontos das diagonais superior direita e inferior esquerda Ou seja para uma malha de tamanho n x n considerando que o índice i representa a linha e o índice j representa a coluna e a origem sendo o vértice superior esquerdo os pontos T1n e Tn1 serão os únicos a serem tratados pela Equação 6 Além disso como o fluxo é nulo a equação tornase mais simples e pode ser descrita pela Equação 7 𝑇𝑚𝑛 𝑇𝑚𝑛1 𝑇𝑚1𝑛𝑅𝑒𝑝𝑟𝑒𝑠𝑒𝑛𝑡𝑎 Figura 4 Representação de um ponto nodal em vértice externo com fluxo térmico uniforme 𝑦 𝑥 q q 2𝑇𝑚𝑛 𝑇𝑚1𝑛 𝑇𝑚𝑛1 0 7 Considerando que a condição de superfície exposta ocorre apenas nos cantos superior direito e inferior esquerdo o caso 2 ponto nodal em superfície plana com fluxo constante aplicase apenas aos seguintes pontos T1n21 até T1n1 T2n até Tn21n Tn211 até Tn11 Tn2 até Tnn21 Ainda por cima a Equação 5 pode ser simplificada porque o fluxo é nulo obtendose a Equação 8 2𝑇𝑚1𝑛 𝑇𝑚𝑛1 𝑇𝑚𝑛1 4𝑇𝑚𝑛 0 8 Todos os demais pontos são classificados como pontos nodais internos valendose a Equação 4 Por fim outra observação ainda importante é que todas as equações que serão empregadas Equações 4 7 e 8 são independentes da condutividade do material Portanto obtémse uma matriz de equações do tipo AT b onde A é uma matriz quadrada de ordem n x n Código para construção da malha Representação da malha Definindo as dimensões da malha nx 8 Número de nós na direção x ny 8 Número de nós na direção y def plotarmalhanx ny nx int Número de nós na direção x ny int Número de nós na direção y Coordenadas para os nós da malha xcoords nplinspace0 1 nx ycoords nplinspace0 1 ny Mesh para todos os nós X Y npmeshgridxcoords ycoords Inverter o eixo Y para garantir o perfil correto Y npflipudY Plotando a figura fig ax pltsubplotsfigsizenx ny Controla o tamanho da Figura axscatterX Y cblue s100 zorder2 Plot dos nós Plotar as linhas da grade for i in rangeny axplotXi Yi k lw1 zorder1 for j in rangenx axplotX j Y j k lw1 zorder1 Nomeando os nós Tmn for j in rangeny for i in rangenx axtextXi j Yi j 002 fTi1j1 haleft vabottom axsettitlefRepresentação da Malha nxxny axsetxlabelCoordenada X m axsetylabelCoordenada Y m axsetxlim01 11 axsetylim01 11 axsetaspectequal adjustablebox axgridTrue pltshow plotarmalhanx ny A Figura 5 apresenta o esquema dos pontos nodais de uma matriz 8x8 Figura 5 Matriz 8 x 8 ilustrativa do problema Código para o cálculo do perfil de temperatura A seguir o código utilizado para calcular o perfil de temperatura As Equações 4 7 e 8 foram utilizadas para se determinar os coeficientes da matriz A Contudo como alguns pontos nodais já possuem a temperatura especificada foi necessário mais um bloco de código para coletar estes casos T1conhecido e T2conhecido e rearranjar as matrizes A e b As linhas do código contém comentários que auxiliam no entendimento do mesmo import numpy as np def resolverperfiltemperaturanx ny T1 3000 Temperatura do contorno superior esquerdo região T2 1500 Temperatura do contorno inferior direito região numeronos nx ny A npzerosnumeronos numeronos b npzerosnumeronos def index1di j return j nx i A partir daqui definese os coeficientes das Equações de cada ponto nodal for j in rangeny Linhas 0 a ny1 for i in rangenx Colunas 0 a nx1 idx index1di j Pontos nodais na CC Região T1 Borda superior j ny1 i nx2 ou borda esquerda i 0 j ny2 if j ny 1 and i nx 2 or i 0 and j ny 2 Aidx idx 10 bidx T1 Região T2 Borda inferior j 0 i nx2 ou borda direita i nx1 j ny2 elif j 0 and i nx 2 or i nx 1 and j ny 2 Aidx idx 10 bidx T2 Pontos nodais em vértice submetido a fluxo uniforme Equação 7 Canto superior direito i nx1 e j ny1 elif i nx 1 and j ny 1 Aidx idx 20 Aidx index1di 1 j 10 Tm1n Aidx index1di j 1 10 Tmn1 Canto inferior esquerdo i 0 e j 0 elif i 0 and j 0 Aidx idx 20 Aidx index1di 1 j 10 Tm1n Aidx index1di j 1 10 Tmn1 Nós com fluxo uniforme na superfície Equação 8 Borda esquerda parte inferior isolada i 0 e 0 j ny2 elif i 0 and 0 j ny 2 Aidx idx 40 Aidx index1di 1 j 20 Tm1n Aidx index1di j 1 10 Tmn1 Aidx index1di j 1 10 Tmn1 Borda direita parte superior isolada i nx1 e ny2 j ny1 elif i nx 1 and ny 2 j ny 1 Aidx idx 40 Aidx index1di 1 j 20 Tm1n Aidx index1di j 1 10 Tmn1 Aidx index1di j 1 10 Tmn1 Borda inferior parte esquerda isolada j 0 e 0 i nx2 elif j 0 and 0 i nx 2 Aidx idx 40 Aidx index1di 1 j 10 Tm1n Aidx index1di 1 j 10 Tm1n Aidx index1di j 1 20 Tmn1 Borda superior parte direita isolada j ny1 e nx2 i nx1 elif j ny 1 and nx 2 i nx 1 Aidx idx 40 Aidx index1di 1 j 10 Tm1n Aidx index1di 1 j 10 Tm1n Aidx index1di j 1 20 Tmn1 Ponto nodal interno Equação 4 else Aidx idx 40 if i 0 Aidx index1di 1 j 10 Tm 1n if i nx 1 Aidx index1di 1 j 10 Tm1n if j 0 Aidx index1di j 1 10 Tmn 1 if j ny 1 Aidx index1di j 1 10 Tmn1 Como alguns pontos nodais já estão definidos devese ajustar o vetor b transferindo a variável conhecida da matriz A para o vetor b Verifica se o nó atual não é um nó conhecido T1conhecido j ny 1 and i nx 2 or i 0 and j ny 2 T2conhecido j 0 and i nx 2 or i nx 1 and j ny 2 if not T1conhecido or T2conhecido Vizinho à esquerda m1 n if i 0 if i 1 0 and j ny 2 or j ny 1 and i 1 nx 2 Vizinho é T1 bidx T1 Aidx index1di 1 j 00 elif i 1 nx 1 and j ny 2 or j 0 and i 1 nx 2 Vizinho é T2 bidx T2 Aidx index1di 1 j 00 Vizinho à direita m1 n if i nx 1 if i 1 0 and j ny 2 or j ny 1 and i 1 nx 2 Vizinho é T1 bidx T1 Aidx index1di 1 j 00 elif i 1 nx 1 and j ny 2 or j 0 and i 1 nx 2 Vizinho é T2 bidx T2 Aidx index1di 1 j 00 Vizinho inferior m n1 if j 0 if j 1 ny 1 and i nx 2 or i 0 and j 1 ny 2 Vizinho é T1 bidx T1 Aidx index1di j 1 00 elif j 1 0 and i nx 2 or i nx 1 and j 1 ny 2 Vizinho é T2 bidx T2 Aidx index1di j 1 00 Vizinho superior m n1 if j ny 1 if j 1 ny 1 and i nx 2 or i 0 and j 1 ny 2 Vizinho é T1 bidx T1 Aidx index1di j 1 00 elif j 1 0 and i nx 2 or i nx 1 and j 1 ny 2 Vizinho é T2 bidx T2 Aidx index1di j 1 00 Resolver o sistema AT b Tvector nplinalgsolveA b Organizando a matriz Tmatrix Tvectorreshapeny nx O comando npflipud serve para deixar a matriz correta sem inversão de i por j Tempmatrix npflipudTmatrix return Tempmatrix Obtendo o perfil de temperatura com a função sucinta nx 8 ny 8 perfiltemperatura resolverperfiltemperaturanx ny printPerfil de temperatura da malha em C Versão Sucinta printnproundperfiltemperatura 2 As Figuras 6 7 8 e 9 apresentam os perfis de temperatura calculados para as matrizes 8x8 20x20 30x30 e 40x40 Algumas características fundamentam a coerência do cálculo Percebese que o perfil apresenta uma linha isotérmica diagonal em 225C O perfil apresenta simetria no primeiro e no terceiro quadrantes A temperatura diminui gradualmente da diagonal superior esquerda até a diagonal inferior direita Figura 6 Perfil de temperatura calculado para a matriz 8x8 Figura 7 Perfil de temperatura calculado para a matriz 20x20 Figura 8 Perfil de temperatura calculado para a matriz 30x30 Figura 9 Perfil de temperatura calculado para a matriz 40x40 Código para o mapa de calor import numpy as np import matplotlibpyplot as plt def heatmapTmatrix nx ny xcoords nplinspace0 1 nx ycoords nplinspace0 1 ny X Y npmeshgridxcoords ycoords fig ax pltsubplotsfigsizenx ny continuidade 200 Controla a continuidade do heatmap c axcontourfX Y Tmatrix continuidade cmapjet vmin150 vmax300 Marcas de escala tickvalues nparange150 300 15 figcolorbarc axax labelTemperatura C tickstickvalues axsettitlefMapa de Calor do Perfil de Temperatura na malha nxxny axsetxlabelCoordenada X m axsetylabelCoordenada Y m axsetaspectequal adjustablebox pltshow nx 8 ny 8 temperaturasfinais resolverperfiltemperaturanx ny heatmapnpflipudtemperaturasfinais nx ny O comando npflipud serve para deixar a figura correta sem inversão i por j A Figura 10 apresenta os mapas de calor para as malhas 8x8 sup esq 20x20 sup dir 30x30 inf esq e 40x40 inf dir Notase um sutil aumento na continuidade das cores com o aumento da malha Figura 10 Mapas de calor para as 4 diferentes malhas estudadas Por fim pedese para calcular o fluxo numérico de calor por unidade de comprimento e comparar com o valor teórico Equação 9 Para calcular o calor numérico devese somar as contribuições que chegamsaem dos pontos nodais Em seguida calculase a média do calor que entra e do calor que sai 𝑞 𝑇1 𝑇2𝑘 𝑊 𝑚 9 Analisando o canto superior esquerdo o calor entra pelos pontos nodais adjascentes por meio da Equação 10 Já analisando o canto inferior direito o calor sai pelos pontos nodais superficiais por meio da Equação 11 Estes dois calores devem ser iguais O delta ficou indefinido porque ele vai ser Δy se o ponto nodal estiver sobre o eixo x e Δx se o ponto nodal estiver sobre o eixo y Como os calores serão somados para cada ponto nodal criouse duas variáveis qnumericoin e qnumericoout para ir acumulando os resultados O valor representativo vai ser a média dos valores que entra e que sai 𝑞𝑛𝑢𝑚 𝑇1 𝑇𝑎𝑑𝑗𝑎𝑠𝑐𝑒𝑛𝑡𝑒 10 𝑞𝑛𝑢𝑚 𝑇𝑎𝑑𝑗𝑎𝑠𝑐𝑒𝑛𝑡𝑒 𝑇2 11 Código para o cálculo do fluxo numérico de calor por unidade de comprimento Calculando o fluxo de calor por unidade de comprimento Parâmetros do Problema k 353 Condutividade térmica WmK T1 3000 Temperatura no canto superior esquerdo C T2 1500 Temperatura no canto inferior direito C Parâmetros da malha nx 10 ny 10 Fluxo de Calor Teórico por Unidade de Comprimento qteorico T1 T2 k Fluxo de Calor Numérico por Unidade de Comprimento Obter o perfil de temperatura da malha perfiltemperatura resolverperfiltemperaturanx ny Calcular os incrementos espaciais placa de 1m x 1m Deltax 10 nx 1 Deltay 10 ny 1 Inicializando as variáveis qnumericoin e qnumericoout qnumericoin 00 qnumericoout 00 Cálculo do fluxo que entra Segmento superior da borda T1 j ny1 i de 0 a nx2 1 Fluxo na direção y entrando na placa for i in rangenx 2 Tadjacenteinterno perfiltemperatura1 i Densidade de fluxo q k Tsuperficie Tadjacenteinterno Deltay qfluxdensity k T1 Tadjacenteinterno Deltay qnumericoin qfluxdensity Deltax Segmento esquerdo da borda T1 i 0 j de ny2 a ny1 Fluxo na direção x entrando na placa for j in rangeny 2 ny Tadjacenteinterno perfiltemperaturany1j 1 Densidade de fluxo q k Tsuperficie Tadjacenteinterno Deltax qfluxdensity k T1 Tadjacenteinterno Deltax qnumericoin qfluxdensity Deltay Cálculo do fluxo que sai Segmento inferior da borda T2 j 0 i de nx2 a nx1 Fluxo na direção y saindo da placa for i in rangenx 2 nx Tadjacenteinterno perfiltemperaturany2 i Densidade de fluxo q k Tadjacenteinterno Tsuperficie Deltay qfluxdensity k Tadjacenteinterno T2 Deltay qnumericoout qfluxdensity Deltax Segmento direito da borda T2 i nx1 j de 0 a ny2 1 Fluxo na direção x saindo da placa for j in rangeny 2 Tadjacenteinterno perfiltemperaturany1j nx2 Densidade de fluxo q k Tadjacenteinterno Tsuperficie Deltax qfluxdensity k Tadjacenteinterno T2 Deltax qnumericoout qfluxdensity Deltay Calcular a média dos fluxos de entrada e saída qnumericomedio qnumericoin qnumericoout 20 Resultados printfFluxo de calor teórico por unidade de comprimento q qteorico2f Wm printfFluxo de calor numérico que entra por unidade de comprimento q qnumericoin2f Wm printfFluxo de calor numérico que sai por unidade de comprimento q qnumericoout2f Wm printfFluxo de calor numérico Médio por unidade de comprimento q qnumericomedio2f Wm A Tabela 1 apresenta os valores obtidos para as malhas e para o cálculo teórico Quanto maior o número de pontos nodais mais próximo do teórico é o valor do calor por unidade de comprimento Isso mostra o ganho de precisão em refinar o tamanho da malha Tabela 1 Resultados do fluxo de calor por unidade de comprimento para o problema Calor por unidade de comprimento Wm Erro Teórico Equação 9 529500 Numérico malha 8x8 366164 3085 Numérico malha 20x20 435045 1784 Numérico malha 30x30 454291 1420 Numérico malha 40x40 465300 1212 Considere uma placa em processo de transferência de calor em regime permanente Figura 1 Figura 1 Esquema ilustrativo da placa em transferência de calor Desejase obter o perfil de temperatura da condução de calor bidimensional utilizando o método das diferenças finitas e comparar o resultado com a solução analítica Dados do problema 1 Placa quadrada com 1m x 1m e condutividade térmica igual a 353 WmK 2 Temperatura uniforme de 300C no canto superior esquerdo e de 150C no canto inferior direito 3 Os cantos inferior esquerdo e superior direito encontramse isolados termicamente 4 O processo ocorre em regime permanente 5 Resolver o problema utilizando malhas 8x8 20x20 30x30 e 40x40 considerando espaçamentos iguais na horizontal e na vertical Δx Δy 6 Comparar o fluxo de calor por unidade de comprimento q em Wm obtido com a resolução numérica com o resultado analítico que é dado por q T1 T2k Hipóteses do problema 1 Regime permanente y x 2 Transferência de calor bidimensional em x e y 3 Propriedades constantes e isotrópicas 4 Sistema sem geração térmica Modelagem do problema O perfil de temperatura na placa é descrito pela equação da difusão do calor Equação 1 Aplicandose as hipóteses simplificadores obtémse a EDP da Equação 2 xk T x y k T y zk T z qG t ρc pT 1 xk T x y k T y zk T z qG t ρc pT ²T x ² ²T y ²02 A EDP da Equação 2 necessita de quatro condições de contorno para ser resolvida Neste problema adicionar as condições de contorno comuns x 0 q 0 x L T T2 y 0 q 0 y L T T1 não vai gerar um resultado correto pois a condição não se mantém a mesma para toda a abscissa ou toda a ordenada Logo convém aplicar o método numérico das diferenças finitas Primeiro devese definir a representação de um nó o menor elemento de aplicação do método Figura 2 Existem vários tipos diferentes de nó Para este problema serão necessárias as equações de diferenças finitas para três tipos ponto nodal interior Figura 2 ponto nodal em uma superfície plana com fluxo uniforme e ponto nodal em um vértice externo fluxo uniforme Isso ocorre em função da condição de contorno a qual o nó é submetido 1 Ponto nodal interno T m1n T mn1 T mn1 T m1 n T mn Figura 2 Representação de um nó y Sobre este nó as derivadas parciais podem ser aproximadas pelo método das diferenças centradas Equação 3 ²T x ² Tm1n2T m nT m1n x² ²T y² Tm n12T m nT m n1 y² 3 Substituindo as derivadas parciais da Equação 3 na EDO da Equação 2 obtémse ²T x ² ²T y ²0 T m1 n2T mnT m1n x² T mn12Tm nTm n1 y² 0 Considerando que as dimensões da malha são iguais ou seja que Δx Δy h obtém se a Equação 4 do ponto nodal interno T m1 nT m1 nT m n1T mn14T m n04 2 Ponto nodal em uma superfície plana com fluxo térmico uniforme Neste caso o calor condutivo que chega dos nós adjascentes deve ser igual ao fluxo de calor uniforme que sai da peça Caso o fluxo de calor esteja entrando na peça conforme ilustrado pela Figura 6 o calor condutivo vai entrar na peça até as peças adjascentes Realizando o balanço de energia obtémse a Equação 5 k y T m1nTm n x k x 2 T mn1T mn y k x 2 Tm n1T mn y q y 0 x T m1n T mn1 T mn1 q T mn y x Figura 3 Representação de um ponto nodal em uma superfície plana com fluxo térmico uniforme 2T m1nT mn1T m n12q x over k 4 T rsub mn 0 5 3 Ponto nodal em vértice externo com fluxo uniforme Neste caso os três nós expostos estão submetidos ao fluxo constante Porém os nós adjascentes enviamrecebem do canto um fluxo condutivo Portanto o fluxo de calor por condução nos nós adjascentes internos deve ser igual ao fluxo de calor por convecção nas superfícies expostas Desenvolvendo o balanço de energia e considerando os intervalos iguais obtémse a Equação 6 k y 2 Tm nT m1n x k x 2 Tm nT mn1 y q left xy right 0 T m1nT mn12q over k left xy right 2 T rsub mn 0 6 Resolução do problema De acordo com a Figura 1 o problema vai apresentar nós dos três tipo Quanto maior a quantidade de nós maior vai ser a precisão do perfil de temperatura pois o Δx e o Δy serão iguais Independente do tamanho da malha haverá apenas dois pontos pertencentes ao caso 3 ponto nodal em vértice com fluxo constante Os pontos das diagonais superior direita e inferior esquerda Ou seja para uma malha de tamanho n x n considerando que o índice i representa a linha e o índice j representa a coluna e a origem sendo o vértice superior esquerdo os pontos T1n e Tn1 serão os únicos a serem tratados pela Equação 6 Além disso como o fluxo é nulo a equação tornase mais simples e pode ser descrita pela Equação 7 T m1nRepresentação deum nó T mn1 T mn Figura 4 Representação de um ponto nodal em vértice externo com fluxo térmico uniforme y x q q 2Tm nT m1nTm n107 Considerando que a condição de superfície exposta ocorre apenas nos cantos superior direito e inferior esquerdo o caso 2 ponto nodal em superfície plana com fluxo constante aplicase apenas aos seguintes pontos T1n21 até T1n1 T2n até Tn21n Tn211 até Tn11 Tn2 até Tnn21 Ainda por cima a Equação 5 pode ser simplificada porque o fluxo é nulo obtendose a Equação 8 2T m1nT mn1T m n14T m n08 Todos os demais pontos são classificados como pontos nodais internos valendose a Equação 4 Por fim outra observação ainda importante é que todas as equações que serão empregadas Equações 4 7 e 8 são independentes da condutividade do material Portanto obtémse uma matriz de equações do tipo AT b onde A é uma matriz quadrada de ordem n x n Código para construção da malha Representação da malha Definindo as dimensões da malha nx 8 Número de nós na direção x ny 8 Número de nós na direção y def plotarmalhanx ny nx int Número de nós na direção x ny int Número de nós na direção y Coordenadas para os nós da malha xcoords nplinspace0 1 nx ycoords nplinspace0 1 ny Mesh para todos os nós X Y npmeshgridxcoords ycoords Inverter o eixo Y para garantir o perfil correto Y npflipudY Plotando a figura fig ax pltsubplotsfigsizenx ny Controla o tamanho da Figura axscatterX Y cblue s100 zorder2 Plot dos nós Plotar as linhas da grade for i in rangeny axplotXi Yi k lw1 zorder1 for j in rangenx axplotX j Y j k lw1 zorder1 Nomeando os nós Tmn for j in rangeny for i in rangenx axtextXi j Yi j 002 fTi1j1 haleft vabottom axsettitlefRepresentação da Malha nxxny axsetxlabelCoordenada X m axsetylabelCoordenada Y m axsetxlim01 11 axsetylim01 11 axsetaspectequal adjustablebox axgridTrue pltshow plotarmalhanx ny A Figura 5 apresenta o esquema dos pontos nodais de uma matriz 8x8 Figura 5 Matriz 8 x 8 ilustrativa do problema Código para o cálculo do perfil de temperatura A seguir o código utilizado para calcular o perfil de temperatura As Equações 4 7 e 8 foram utilizadas para se determinar os coeficientes da matriz A Contudo como alguns pontos nodais já possuem a temperatura especificada foi necessário mais um bloco de código para coletar estes casos T1conhecido e T2conhecido e rearranjar as matrizes A e b As linhas do código contém comentários que auxiliam no entendimento do mesmo import numpy as np def resolverperfiltemperaturanx ny T1 3000 Temperatura do contorno superior esquerdo região T2 1500 Temperatura do contorno inferior direito região numeronos nx ny A npzerosnumeronos numeronos b npzerosnumeronos def index1di j return j nx i A partir daqui definese os coeficientes das Equações de cada ponto nodal for j in rangeny Linhas 0 a ny1 for i in rangenx Colunas 0 a nx1 idx index1di j Pontos nodais na CC Região T1 Borda superior j ny1 i nx2 ou borda esquerda i 0 j ny2 if j ny 1 and i nx 2 or i 0 and j ny 2 Aidx idx 10 bidx T1 Região T2 Borda inferior j 0 i nx2 ou borda direita i nx1 j ny2 elif j 0 and i nx 2 or i nx 1 and j ny 2 Aidx idx 10 bidx T2 Pontos nodais em vértice submetido a fluxo uniforme Equação 7 Canto superior direito i nx1 e j ny1 elif i nx 1 and j ny 1 Aidx idx 20 Aidx index1di 1 j 10 Tm1n Aidx index1di j 1 10 Tmn1 Canto inferior esquerdo i 0 e j 0 elif i 0 and j 0 Aidx idx 20 Aidx index1di 1 j 10 Tm1n Aidx index1di j 1 10 Tmn1 Nós com fluxo uniforme na superfície Equação 8 Borda esquerda parte inferior isolada i 0 e 0 j ny2 elif i 0 and 0 j ny 2 Aidx idx 40 Aidx index1di 1 j 20 Tm1n Aidx index1di j 1 10 Tmn1 Aidx index1di j 1 10 Tmn1 Borda direita parte superior isolada i nx1 e ny2 j ny1 elif i nx 1 and ny 2 j ny 1 Aidx idx 40 Aidx index1di 1 j 20 Tm1n Aidx index1di j 1 10 Tmn1 Aidx index1di j 1 10 Tmn1 Borda inferior parte esquerda isolada j 0 e 0 i nx2 elif j 0 and 0 i nx 2 Aidx idx 40 Aidx index1di 1 j 10 Tm1n Aidx index1di 1 j 10 Tm1n Aidx index1di j 1 20 Tmn1 Borda superior parte direita isolada j ny1 e nx2 i nx1 elif j ny 1 and nx 2 i nx 1 Aidx idx 40 Aidx index1di 1 j 10 Tm1n Aidx index1di 1 j 10 Tm1n Aidx index1di j 1 20 Tmn1 Ponto nodal interno Equação 4 else Aidx idx 40 if i 0 Aidx index1di 1 j 10 Tm 1n if i nx 1 Aidx index1di 1 j 10 Tm1n if j 0 Aidx index1di j 1 10 Tmn 1 if j ny 1 Aidx index1di j 1 10 Tmn1 Como alguns pontos nodais já estão definidos devese ajustar o vetor b transferindo a variável conhecida da matriz A para o vetor b Verifica se o nó atual não é um nó conhecido T1conhecido j ny 1 and i nx 2 or i 0 and j ny 2 T2conhecido j 0 and i nx 2 or i nx 1 and j ny 2 if not T1conhecido or T2conhecido Vizinho à esquerda m1 n if i 0 if i 1 0 and j ny 2 or j ny 1 and i 1 nx 2 Vizinho é T1 bidx T1 Aidx index1di 1 j 00 elif i 1 nx 1 and j ny 2 or j 0 and i 1 nx 2 Vizinho é T2 bidx T2 Aidx index1di 1 j 00 Vizinho à direita m1 n if i nx 1 if i 1 0 and j ny 2 or j ny 1 and i 1 nx 2 Vizinho é T1 bidx T1 Aidx index1di 1 j 00 elif i 1 nx 1 and j ny 2 or j 0 and i 1 nx 2 Vizinho é T2 bidx T2 Aidx index1di 1 j 00 Vizinho inferior m n1 if j 0 if j 1 ny 1 and i nx 2 or i 0 and j 1 ny 2 Vizinho é T1 bidx T1 Aidx index1di j 1 00 elif j 1 0 and i nx 2 or i nx 1 and j 1 ny 2 Vizinho é T2 bidx T2 Aidx index1di j 1 00 Vizinho superior m n1 if j ny 1 if j 1 ny 1 and i nx 2 or i 0 and j 1 ny 2 Vizinho é T1 bidx T1 Aidx index1di j 1 00 elif j 1 0 and i nx 2 or i nx 1 and j 1 ny 2 Vizinho é T2 bidx T2 Aidx index1di j 1 00 Resolver o sistema AT b Tvector nplinalgsolveA b Organizando a matriz Tmatrix Tvectorreshapeny nx O comando npflipud serve para deixar a matriz correta sem inversão de i por j Tempmatrix npflipudTmatrix return Tempmatrix Obtendo o perfil de temperatura com a função sucinta nx 8 ny 8 perfiltemperatura resolverperfiltemperaturanx ny printPerfil de temperatura da malha em C Versão Sucinta printnproundperfiltemperatura 2 As Figuras 6 7 8 e 9 apresentam os perfis de temperatura calculados para as matrizes 8x8 20x20 30x30 e 40x40 Algumas características fundamentam a coerência do cálculo Percebese que o perfil apresenta uma linha isotérmica diagonal em 225C O perfil apresenta simetria no primeiro e no terceiro quadrantes A temperatura diminui gradualmente da diagonal superior esquerda até a diagonal inferior direita Figura 6 Perfil de temperatura calculado para a matriz 8x8 Figura 7 Perfil de temperatura calculado para a matriz 20x20 Figura 8 Perfil de temperatura calculado para a matriz 30x30 Figura 9 Perfil de temperatura calculado para a matriz 40x40 Código para o mapa de calor import numpy as np import matplotlibpyplot as plt def heatmapTmatrix nx ny xcoords nplinspace0 1 nx ycoords nplinspace0 1 ny X Y npmeshgridxcoords ycoords fig ax pltsubplotsfigsizenx ny continuidade 200 Controla a continuidade do heatmap c axcontourfX Y Tmatrix continuidade cmapjet vmin150 vmax300 Marcas de escala tickvalues nparange150 300 15 figcolorbarc axax labelTemperatura C tickstickvalues axsettitlefMapa de Calor do Perfil de Temperatura na malha nxxny axsetxlabelCoordenada X m axsetylabelCoordenada Y m axsetaspectequal adjustablebox pltshow nx 8 ny 8 temperaturasfinais resolverperfiltemperaturanx ny heatmapnpflipudtemperaturasfinais nx ny O comando npflipud serve para deixar a figura correta sem inversão i por j A Figura 10 apresenta os mapas de calor para as malhas 8x8 sup esq 20x20 sup dir 30x30 inf esq e 40x40 inf dir Notase um sutil aumento na continuidade das cores com o aumento da malha Figura 10 Mapas de calor para as 4 diferentes malhas estudadas Por fim pedese para calcular o fluxo numérico de calor por unidade de comprimento e comparar com o valor teórico Equação 9 Para calcular o calor numérico devese somar as contribuições que chegamsaem dos pontos nodais Em seguida calculase a média do calor que entra e do calor que sai q T 1T2k W m 9 Analisando o canto superior esquerdo o calor entra pelos pontos nodais adjascentes por meio da Equação 10 Já analisando o canto inferior direito o calor sai pelos pontos nodais superficiais por meio da Equação 11 Estes dois calores devem ser iguais O delta ficou indefinido porque ele vai ser Δy se o ponto nodal estiver sobre o eixo x e Δx se o ponto nodal estiver sobre o eixo y Como os calores serão somados para cada ponto nodal criouse duas variáveis qnumericoin e qnumericoout para ir acumulando os resultados O valor representativo vai ser a média dos valores que entra e que sai qnum T 1T adjascente 10 qnum T adjascenteT 2 11 Código para o cálculo do fluxo numérico de calor por unidade de comprimento Calculando o fluxo de calor por unidade de comprimento Parâmetros do Problema k 353 Condutividade térmica WmK T1 3000 Temperatura no canto superior esquerdo C T2 1500 Temperatura no canto inferior direito C Parâmetros da malha nx 10 ny 10 Fluxo de Calor Teórico por Unidade de Comprimento qteorico T1 T2 k Fluxo de Calor Numérico por Unidade de Comprimento Obter o perfil de temperatura da malha perfiltemperatura resolverperfiltemperaturanx ny Calcular os incrementos espaciais placa de 1m x 1m Deltax 10 nx 1 Deltay 10 ny 1 Inicializando as variáveis qnumericoin e qnumericoout qnumericoin 00 qnumericoout 00 Cálculo do fluxo que entra Segmento superior da borda T1 j ny1 i de 0 a nx2 1 Fluxo na direção y entrando na placa for i in rangenx 2 Tadjacenteinterno perfiltemperatura1 i Densidade de fluxo q k Tsuperficie Tadjacenteinterno Deltay qfluxdensity k T1 Tadjacenteinterno Deltay qnumericoin qfluxdensity Deltax Segmento esquerdo da borda T1 i 0 j de ny2 a ny1 Fluxo na direção x entrando na placa for j in rangeny 2 ny Tadjacenteinterno perfiltemperaturany1j 1 Densidade de fluxo q k Tsuperficie Tadjacenteinterno Deltax qfluxdensity k T1 Tadjacenteinterno Deltax qnumericoin qfluxdensity Deltay Cálculo do fluxo que sai Segmento inferior da borda T2 j 0 i de nx2 a nx1 Fluxo na direção y saindo da placa for i in rangenx 2 nx Tadjacenteinterno perfiltemperaturany2 i Densidade de fluxo q k Tadjacenteinterno Tsuperficie Deltay qfluxdensity k Tadjacenteinterno T2 Deltay qnumericoout qfluxdensity Deltax Segmento direito da borda T2 i nx1 j de 0 a ny2 1 Fluxo na direção x saindo da placa for j in rangeny 2 Tadjacenteinterno perfiltemperaturany1j nx2 Densidade de fluxo q k Tadjacenteinterno Tsuperficie Deltax qfluxdensity k Tadjacenteinterno T2 Deltax qnumericoout qfluxdensity Deltay Calcular a média dos fluxos de entrada e saída qnumericomedio qnumericoin qnumericoout 20 Resultados printfFluxo de calor teórico por unidade de comprimento q qteorico2f Wm printfFluxo de calor numérico que entra por unidade de comprimento q qnumericoin2f Wm printfFluxo de calor numérico que sai por unidade de comprimento q qnumericoout2f Wm printfFluxo de calor numérico Médio por unidade de comprimento q qnumericomedio2f Wm A Tabela 1 apresenta os valores obtidos para as malhas e para o cálculo teórico Quanto maior o número de pontos nodais mais próximo do teórico é o valor do calor por unidade de comprimento Isso mostra o ganho de precisão em refinar o tamanho da malha Tabela 1 Resultados do fluxo de calor por unidade de comprimento para o problema Calor por unidade de comprimento Wm Erro Teórico Equação 9 529500 Numérico malha 8x8 366164 3085 Numérico malha 20x20 435045 1784 Numérico malha 30x30 454291 1420 Numérico malha 40x40 465300 1212