Análise Numérica de Problemas Térmicos, Tridimensionais ...repositorio.roca.utfpr.edu.br/jspui/bitstream/1/... · de Conclusão de Curso 2, do Curso de En-genharia Mecânica do

UNIVERSIDADE TECNOLÓGICA FEDERAL DO PARANÁDEPARTAMENTO ACADÊMICO DE MECÂNICA

CURSO DE ENGENHARIA MECÂNICA

ÍGOR SANTANA SOUSA

ANÁLISE NUMÉRICA DE PROBLEMAS TÉRMICOSTRIDIMENSIONAIS BASEADA EM MÉTODO

IMPLÍCITO DE AVANÇO NO TEMPO COM APLICAÇÃOEM DISSIPADORES DE CALOR PASSIVOS

TRABALHO DE CONCLUSÃO DE CURSO

PATO BRANCO2017

ÍGOR SANTANA SOUSA

ANÁLISE NUMÉRICA DE PROBLEMAS TÉRMICOSTRIDIMENSIONAIS BASEADA EM MÉTODO

IMPLÍCITO DE AVANÇO NO TEMPO COM APLICAÇÃOEM DISSIPADORES DE CALOR PASSIVOS

Trabalho de conclusão de curso de gradu-ação, apresentado à disciplina de Trabalhode Conclusão de Curso 2, do Curso de En-genharia Mecânica do Departamento Aca-dêmico de Mecânica - DAMEC - da Universi-dade Tecnológica Federal do Paraná , comorequisito parcial para obtenção do título deBacharel em Engenharia Mecânica.

Orientador: Prof. Dr. Francisco AugustoAparecido Gomes

PATO BRANCO2017

FOLHA DE APROVAÇÃO

Análise Numérica de Problemas Térmicos Tridimensionais Baseada em MétodoImplícito de Avanço no Tempo com Aplicação em Dissipadores de Calor

Passivos

Ígor Santana Sousa

Trabalho de Conclusão de Curso de Graduação apresentado no dia 16/11/2017 comorequisito parcial para a obtenção do Título de Engenheiro Mecânico, do curso de Enge-nharia Mecânica do Departamento Acadêmico de Mecânica (DAMEC) da UniversidadeTecnológica Federal do Paraná - Câmpus Pato Branco (UTFPR-PB). O candidato foiarguido pela Banca Examinadora composta pelos professores abaixo assinados. Apósdeliberação, a Banca Examinadora julgou o trabalho APROVADO.

Prof. Dr. Fabiano Ostapiv(UTFPR - Departamento Acadêmico de Mecânica)

Prof. Dr. João Biesdorf(UTFPR - Departamento Acadêmico de Matemática)

Prof. Dr. Francisco Augusto Aparecido Gomes(UTFPR - Departamento Acadêmico de Mecânica)

Orientador

Prof. Dr. Bruno Bellini MedeirosResponsável pelo TCC do Curso de Eng. Mecânica

A Folha de Aprovação assinada encontra-se na Coordenação do Curso de Engenharia Mecânica

Este trabalho é dedicado aos poucos que fazem jus as oportunidades que tem,enquanto muitos as desperdiçam

e milhares nunca as tiveram mas com elas sempre sonharam.

AGRADECIMENTOS

São tantas pessoas para agradecer que seria hipócrita achar que todas poderiamcitadas aqui, pois muitas eu nem sei o nome, mas algumas são tão marcantes que nãopoderiam ser esquecidas.

Agradeço primeiramente a minha família, sem a qual nunca teria iniciado o meuensino superior, pelo apoio, garra em encarar todas as dificuldades de frente semnunca se deixar abater ou desistir, com atitudes que vão desde um café ou palavraamiga nos dias de cansaço à mudança para um estado longínquo deixando toda a vidapara trás.

Agradeço também a todos os meus amigos pelo suporte, ajuda, companhei-rismo nestes cinco anos de graduação, seja durante as festas ou nas noites viradasresolvendo listas de exercícios.

Não poderia deixar de agradecer a alguns professores fantásticos do curso queforam além da obrigação e se empenharam em desenvolver profissionais melhorespara o futuro.

Dentre estes profissionais, agradeço ao meu orientador pelo suporte durante arealização dos trabalhos, seja durante a iniciação científica ou durante a execuçãodeste trabalho, com paciência sempre de estando acessível e sendo assertivo.

Ademas, agradeço a todos as pessoas que contribuíram para o meu crescimentopessoal, mesmo que minimamente, tornando possível a minha graduação, possibili-tando a mim tentar retribuir a sociedade toda a estima que me foi investida.

”There have to be reasons that you get up in the morningand you want to live. Why do you want to live? What’s the point?

What inspires you? What do you love about the future?If the future does not include being out there among the stars

and being a multi-planet species,I find that incredibly depressing.”

(Elon Musk)

RESUMO

SOUSA, Ígor Santana. Análise Numérica de Problemas Térmicos TridimensionaisBaseada em Método Implícito de Avanço no Tempo com Aplicação em Dissipadoresde Calor Passivos. 105 f. Trabalho de Conclusão de Curso - Curso de EngenhariaMecânica, Universidade Tecnológica Federal do Paraná. Pato Branco, 2017.

Esse trabalho consiste na expansão de um código de capaz de simular problemasbidimensionais transientes de transferência de calor utilizando o método numérico dasdiferenças finitas, desenvolvido por CONCI e GOMES (2017), para que ele seja capazde resolver problemas tridimensionais com geometrias mais complexas, juntamentecom o desenvolvimento de um código de pré-condicionamento de malhas para omesmo. Para a análise é utilizado o método iterativo GMRES com o pré-condicionadorMILU0. O código é implementado, verificado e validado, alcançando sexta ordemde precisão, e, em seguida, aplicado em geometrias de dissipadores de calor paraatestar sua eficiência, possibilitando a sua utilização na análise futura de geometriasde dissipadores de calor.

Palavras-chave: método das diferenças finitas, método implícito, GMRES, dissipadoresde calor, tridimensional.

ABSTRACT

SOUSA, Ígor Santana. Numerical Analysys of Tridimentional Thermal Problemas Basedon Implicit Method of Advance in Time with Application on Heatsinks. 105 p. Trabalhode Conclusão de Curso - Curso de Engenharia Mecânica, Universidade TecnológicaFederal do Paraná. Pato Branco, 2017.

This work consists on the expansion of an existing code able to simulate bidimentionaltransient problems of heat transfer using the numerical finite difference method, devel-oped by CONCI e GOMES (2017), so that it can be capable of solving tridimentionalproblems with more complex geometries, along with the development of a code that pre-conditionate meshs for the same. For the analisys it is used the iterative method GMRESwith the préconditioner MILU0. The code is implementated, verified and validated, reach-ing the sixth order of precision, and, then, aplied in geometries of heatsinks to attest itsefficiency, making possible its utilization in future analisys of heatsink geometries.

Keywords: finite differece method, implicit method, GMRES, heatsinks, tridimentional.

LISTA DE FIGURAS

Figura 1 – Balanço de energia para um volume de controle . . . . . . . . . . . 24Figura 2 – Exemplo de malhas . . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figura 3 – Esquemático do código de geração e pré-condicionamento da malha 34Figura 4 – Identificação da posição dos 26 vizinhos de cada nó . . . . . . . . . 35Figura 5 – Identificação das 41 posições geométricas possíveis para um nó . . 36Figura 6 – Esquemático dos nós segundo as suas posições . . . . . . . . . . . 38Figura 7 – Esquemático do código de montagem, pré-condicionamento e solu-

ção do sistema linear . . . . . . . . . . . . . . . . . . . . . . . . . . 42Figura 8 – Representação dos dissipadores de calor a serem analisados nesse

trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47Figura 9 – Condição de contorno nas faces da geometria do caso 1 . . . . . . 48Figura 10 – Condição de contorno nas faces da geometria do caso 2 . . . . . . 48Figura 11 – Condição de contorno nas faces da geometria do caso 3 . . . . . . 49Figura 12 – Malhas utilizadas no caso 1 . . . . . . . . . . . . . . . . . . . . . . . 51Figura 13 – Solução aproximada do caso 1 . . . . . . . . . . . . . . . . . . . . . 52Figura 14 – Plotagem das normas em função do distanciamento dos nós para o

caso 1 para os respectivos Laplacianos . . . . . . . . . . . . . . . . 53Figura 15 – Diferença entre a solução aproximada e a exata, para cada malha,

utilizando os três Laplacianos na linha y = 0,5 Ly e z = 0,5 Lz . . . . 54Figura 16 – Malhas utilizadas no caso 2 . . . . . . . . . . . . . . . . . . . . . . . 55Figura 17 – Solução aproximada do caso 2 . . . . . . . . . . . . . . . . . . . . . 55Figura 18 – Malhas utilizadas no caso 3 . . . . . . . . . . . . . . . . . . . . . . . 56Figura 19 – Solução aproximada do caso 3 . . . . . . . . . . . . . . . . . . . . . 56Figura 20 – A diferença entre a solução exata e a solução aproximada encontrada

para cada malha do caso 4 na linha x = 0,25, y = 0,25 . . . . . . . . 57Figura 21 – Malhas utilizadas no caso 4 . . . . . . . . . . . . . . . . . . . . . . . 58Figura 22 – Solução aproximada do caso 4 . . . . . . . . . . . . . . . . . . . . . 58Figura 23 – Diferença entre a solução exata e a solução aproximada, para cada

malha do caso 4, ao longo da reta x = y = 0,5 . . . . . . . . . . . . . 59Figura 24 – Solução aproximada para o quinto caso em alguns instantes de tempo 59Figura 25 – Solução aproximada para o quinto caso em alguns instantes de

tempo (continuação) . . . . . . . . . . . . . . . . . . . . . . . . . . . 60Figura 26 – Solução aproximada para o quinto caso em alguns instantes de tempo 60Figura 27 – Comportamento da temperatura máxima no dominio ao longo dos

passos de tempo, com ∆t = 0,1 s . . . . . . . . . . . . . . . . . . . . 61Figura 28 – Distribuição de temperatura calculada para a primeira geometria de

aplicação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Figura 29 – Erro entre os passos de tempo subsequente (a) e e evolução dastemperaturas limites na geometria ao longo dos passos de tempo (b)para a primeira aplicação . . . . . . . . . . . . . . . . . . . . . . . . 63

Figura 30 – Distribuição de temperatura calculada para a segunda geometria deaplicação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

Figura 31 – Erro entre os passos de tempo subsequente (a) e e evolução dastemperaturas limites na geometria ao longo dos passos de tempo (b)para a segunda aplicação . . . . . . . . . . . . . . . . . . . . . . . . 64

Figura 32 – Distribuição de temperatura calculada para a terceira geometria deaplicação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

Figura 33 – Erro entre os passos de tempo subsequente (a) e e evolução dastemperaturas limites na geometria ao longo dos passos de tempo (b)para a aplicação 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

LISTA DE TABELAS

Tabela 1 – Características geométricas dos dissipadores de calor a serem estu-dados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

Tabela 2 – Valores das normas e coeficientes angulares para o Laplaciano de 7nós . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Tabela 3 – Valores das normas e coeficientes angulares para o Laplaciano de19 nós . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Tabela 4 – Valores das normas e coeficientes angulares para o Laplaciano de27 nós . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Tabela 5 – Comparação entre os valores calculados e os valores de referência 55Tabela 6 – Parâmetros calculados após a convergência para as geometrias

propostas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

SUMÁRIO

1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211.1 Delimitação do tema . . . . . . . . . . . . . . . . . . . . . . . . . . . 211.2 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221.3 Organização do trabalho . . . . . . . . . . . . . . . . . . . . . . . . 222 REFERENCIAL TEÓRICO . . . . . . . . . . . . . . . . . . . . . . . . 232.1 Equação fundamental da conservação de energia . . . . . . . . . 232.2 Solução numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.2.1 Comparação com outros métodos de solução . . . . . . . . . . . . . 252.2.2 Erros associados a solução numérica . . . . . . . . . . . . . . . . . . 252.3 Malha de representação do domínio . . . . . . . . . . . . . . . . . 262.4 Condições de contorno . . . . . . . . . . . . . . . . . . . . . . . . . 262.5 Método das diferenças finitas e o balanço de energia . . . . . . . 272.6 Discretização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.7 Método implícito de avanço no tempo . . . . . . . . . . . . . . . . 282.7.1 Classificação de problemas quanto ao comportamento . . . . . . . . 282.7.2 Métodos de avanço no tempo e estabilidade . . . . . . . . . . . . . . 282.8 Solução de sistemas lineares . . . . . . . . . . . . . . . . . . . . . 292.8.1 Métodos iterativos de solução . . . . . . . . . . . . . . . . . . . . . . 292.8.2 Método do resíduo mínimo generalizado - GMRES . . . . . . . . . . 292.8.3 Pré-condicionadores . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.9 Programação em linguagem FORTRAN . . . . . . . . . . . . . . . 302.10 Verificação e validação de códigos . . . . . . . . . . . . . . . . . . 302.11 Dissipadores de calor . . . . . . . . . . . . . . . . . . . . . . . . . . 312.11.1 Funcionamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.11.2 Utilização em componentes eletrônicos . . . . . . . . . . . . . . . . . 312.11.3 Dissipadores de calor passivos . . . . . . . . . . . . . . . . . . . . . 323 METODOLOGIA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.1 Geração da malha . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.2 Código de pré-condicionamento da malha . . . . . . . . . . . . . 333.2.1 Identificação dos vizinhos dos nós . . . . . . . . . . . . . . . . . . . . 333.2.2 Classificação geométrica dos nós . . . . . . . . . . . . . . . . . . . . 343.2.3 Classificação física dos nós . . . . . . . . . . . . . . . . . . . . . . . 353.2.4 Espelhamento dos nós de fronteiras adiabáticas . . . . . . . . . . . . 373.2.5 Classificação dos nós quanto aos vizinhos que possui . . . . . . . . 373.3 Discretização das equações físicas que regem o problema . . . 373.3.1 Discretização espacial dos nós de interior . . . . . . . . . . . . . . . 383.3.2 Discretização espacial dos nós de fronteira . . . . . . . . . . . . . . . 39

3.3.3 Discretização temporal implícita de primeira ordem de precisão . . . 393.4 Montagem da matriz . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.5 Código utilizado na simulação . . . . . . . . . . . . . . . . . . . . . 423.6 Verificação e validação do código . . . . . . . . . . . . . . . . . . . 443.7 Casos a serem estudados . . . . . . . . . . . . . . . . . . . . . . . 443.7.1 Casos estacionários com solução exata . . . . . . . . . . . . . . . . 443.7.2 Caso transiente com solução exata . . . . . . . . . . . . . . . . . . . 463.7.3 Aplicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474 VERIFICAÇÃO E VALIDAÇÃO DO CÓDIGO . . . . . . . . . . . . . . 514.1 Casos estacionários . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.1.1 Caso 1 - Bidimensional com solução exata . . . . . . . . . . . . . . . 514.1.2 Caso 2 - Bidimensional com solução tabelada . . . . . . . . . . . . . 534.1.3 Caso 3 - Tridimensional com solução exata . . . . . . . . . . . . . . . 554.1.4 Caso 4 - Tridimensional com solução exata . . . . . . . . . . . . . . . 574.2 Caso transiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.2.1 Caso 5 - Bidimensional com solução exata . . . . . . . . . . . . . . . 575 RESULTADOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 635.1 Aplicação 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 635.2 Aplicação 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 645.3 Aplicação 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 656 CONCLUSÕES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

REFERÊNCIAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69APÊNDICE A – ARQUIVOS DE ENTRADA DO CÓDIGO DE GE-

RAÇÃO DE MALHA . . . . . . . . . . . . . . . . . 71A.1 Arquivo lnd.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71A.2 Arquivo lke.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71A.3 Arquivo bdf.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73A.4 Arquivos bcf.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

APÊNDICE B – CÓDIGO DE GERAÇÃO E PRÉ-CONDICIONAMENTODE MALHA . . . . . . . . . . . . . . . . . . . . . . 85

B.1 Arquivo mesh_parameters.f90 . . . . . . . . . . . . . . . . . . . . . 85B.2 Arquivo mesh_generation.f90 . . . . . . . . . . . . . . . . . . . . . 86

APÊNDICE C – GERAÇÃO E PRÉ-CONDICIONAMENTO DA MA-LHA . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

APÊNDICE D – DETALHAMENTO DAS GEOMETRIAS DE DISSI-PADORES DE CALOR UTILIZADAS . . . . . . . . 103

21

1 INTRODUÇÃO

A dissipação de calor em componentes eletrônicos é uma área na qual tem-se dado bastante atenção. Com o desenvolvimento tecnológico, os componenteseletrônicos tiveram seu tamanho diminuído, necessitando de formas cada vez maiseficientes de dissipar o calor produzido, afim de evitar falhas devido a temperaturado componente estar fora da faixa de trabalho. Parte desses esforços são destinadosao desenvolvimento de geometrias para dissipadores de calor, que frequentementeprecisam ser ensaiadas, com a ajuda de métodos numéricos por exemplo, para verificara sua eficiência de funcionamento antes da eventual fabricação.

Métodos numéricos são utilizados em problemas de transferência de calor quandoestes são muito complexos ou dispendem muito tempo para serem resolvidos analitica-mente. Eles permitem uma solução aproximada, mas inexata, do problema mas quepossui precisão suficiente na maioria dos casos. Existem vários métodos de soluçãonumérica, um amplamente utilizado é o método das diferenças finitas, que permite umaaproximação simples do comportamento físico do problema quando este aceita umarepresentação na forma de uma malha estruturada.

1.1 Delimitação do tema

Neste trabalho são abordados problemas tridimensionais de transferência de calorcom comportamento transiente e condições de contorno bem definidas que envolvemcondução, convecção ou fluxo de calor prescrito.

Essa abordagem é feita através da discretização de equações obtidas por meio dobalanço de energia, feita com segunda e sexta ordem de precisão espacial e primeirade precisão temporal, sendo essa discretização implementada em um código de análisenumérica.

O código utiliza o método dos resíduos mínimos generalizados para a convergên-cia da solução, com o auxílio do pré-condicionador MILU0, da família ILU, implementadoanteriormente por CONCI e GOMES (2017).

A validação do código é dada por meio da análise numérica de casos estacionáriose transientes com solução analítica conhecida, sendo feita a análise, verificação evalidação de ambas as ordens de precisão espacial.

Após essas análises, o código é utilizado no cálculo da distribuição de temperaturaem três geometrias de dissipadores de calor submetidos a condições de funcionamentosemelhantes, buscando descobrir qual geometria é mais eficiente, verificando assim aefetividade e aplicabilidade código em aplicações de eventual interesse.

22

1.2 Objetivos

Tem-se como objetivo o desenvolvimento de um código aberto, em linguagemFORTRAN, capaz de fazer simulação numérica com alta ordem de precisão de proble-mas tridimensionais transientes com geometria composta por uma malha estruturadapré-conficionada de elementos igualmente espaçados, tudo prezando um baixo custocomputacional. O gerador de malhas também será desenvolvido neste trabalho e, juntocom o código principal, será empregado na solução numérica de problemas difusivos,como por exemplo, em trocadores de calor passivos.

1.3 Organização do trabalho

No segundo capítulo, é feita uma breve revisão bibliográfica dos assuntos rela-cionados ao escopo desse trabalho. Em seguida, no terceiro capítulo, é apresentadatoda a metologia utilizada na criação do código de geração e pré-condicionamentodas malhas, e nas modificações feitas no código de análise numérica já existente,incluindo o detalhamento dos métodos de discretização utilizados e dos casos deverificação/validação código e dos casos de aplicação.

No quarto capítulo, são apresentados os resultados obtidos durante a verificaçãoe validação do código, junto com uma análise dos deses resultados. No quinto capítulo,são expostos os resultados obtidos na análise dos casos de aplicação seguidos pelaanálise de todos os resultados obtidos durante o trabalho, conclusões e sugestões paraaplicações futuras.

23

2 REFERENCIAL TEÓRICO

Neste capítulo será apresentada uma breve revisão bibliográfica do assuntotratado nesse trabalho, além de ser discutida qual a classificação dos problemas,métodos de abordagem e justificativa de escolha de cada método de acordo com aliteratura de referência em cada assunto.

2.1 Equação fundamental da conservação de energia

As leis de conservação regem os fenômenos que ocorrem na natureza e, desdesua formulação, têm sido fundamentais para a resolução de vários problemas queantes eram inexplicáveis. Algumas das leis de conservação são aproximadas e sósão verdade sob certas circunstâncias, mas a maioria das leis de conservação sãoexatas, isto é, nunca foi provada a sua violação. Entre as leis de conservação exatasestão as leis de conservação de momento linear/angular e a lei de conservação demassa/energia.

Os fenômenos de transferência de calor são governados pela equação da con-servação de energia (Equação 2.1). Essa equação descreve a evolução temporal daenergia no sistema, além da transferência de calor por condução, convecção, radiação,dissipação viscosa, trabalho executado pela pressão e fontes volumétricas de energia.(MINKOWYCZ, 1988)

ρ∂

∂t(ρe) + ∇ • (ρVe) = ∇ • [(k + kt)∇T ] + ∇ • (τ • V)−∇ • (pV) + Sr + Sh (2.1)

Onde:e é a energia total por unidade de volume;k é a condutividade térmica;kt é a condutividade térmica devido ao modelo de turbulência;∇ é o operador gradiente;p é a pressão;ρ é a densidade;Sh é o termo fonte devido as reações;Sr é o termo fonte volumétrico devido a transferência de calor por radiação;T é a temperatura;t é o tempo;τ é o tensor de tensões;V é o vetor velocidade.

24

Em problemas de transferência de calor em sólidos, a dissipação viscosa, trabalhorealizado pela pressão e a condutividade térmica devido ao modelo de turbulênciapodem ser desconsiderados. Para uma maior simplificação dos problemas, o termofonte devido a transferência de calor por radiação também pode ser desconsiderado.Após estas considerações, os termos restantes da equação podem ser escritos naforma de um balanço de energia para um volume de controle. Quando esse balançode energia (Equação 2.1) é escrito em coordenadas retangulares, sob a forma de taxade variação, resulta na Equação 2.2, esquematizada na Figura 1. (CENGEL; GHAJAR,2012; INCROPERA, 2008)

Figura 1 – Balanço de energia para um volume de controle

EEE

ΔEsaient

ger

acu

Fonte: autoria própria.

∆Ėacu = Ėent − Ėsai + Ėger (2.2)

Onde:Ėacu é a energia acumulada no sistema;Ėent é a energia que entra no sistema;Ėsai é a energia que sai do sistema;Ėger é a energia gerada no sistema.

Na Equação 2.2 são levados em conta as formas de entrada (como exemplos atransferência de calor por convecção e um fluxo prescrito de calor em uma superfície),saída (como exemplo a transferência de calor por condução), geração de energia nosistema (como exemplo uma reação química exotérmica interna) e a variação temporalda energia do sistema.

2.2 Solução numérica

Com o avanço na tecnologia computacional, muitos problemas de engenhariaque eram muito complicados e dispendiosos para serem resolvidos analiticamente,agora podem ser resolvidos em uma fração de tempo. Esse avanço permitiu queproblemas possam ser modelados de forma cada vez mais condizente com a realidade.A solução numérica, por meio da simulação numérica, permitiu um gigantesco avançona engenharia durante as ultimas décadas e, embora não seja exata, permite umaaproximação em relação aos resultados reais que é muitas vezes suficiente para oseu fim. Além disso, não é incomum empresas fazerem simulações com milhões de

25

variáveis desconhecidas com velocidade anteriormente inimaginável. (PLETCHER;TANNEHILL; ANDERSON, 1997; MINKOWYCZ, 1988)

2.2.1 Comparação com outros métodos de solução

Quando comparada aos métodos anteriores de solução de problemas, os métodosanalítico e experimental, a solução numérica apresenta diversas vantagens. A soluçãonumérica não possui as mesmas limitações da solução experimental, que além de sercara em determinadas aplicações, muitas vezes tem de ser executadas em uma escaladiferente da real, não reproduzindo o problema com muita fidelidade além de apresentardificuldade com relação a medições dos resultados do experimento. Já o métodonumérico pode resolver o problema em escala real sem estar sujeita a essas limitações.Quando comparada a solução analítica, problemas com geometrias e condições decontorno mais complexas podem ser resolvidos, com poucas simplificações, alémde reduzir o tempo de cálculo, uma vez que são executados por processadores quesó tem aumentado a sua velocidade ao longo do tempo. (PLETCHER; TANNEHILL;ANDERSON, 1997)

2.2.2 Erros associados a solução numérica

Uma das desvantagens da utilização de uma solução numérica é que apesar deser muito precisa, a solução numérica é inexata devido aos erros associados a essetipo de solução e ao método computacional. Devido ao método de armazenamento denúmeros dos computadores digitais, estes apresentam erros de arredondamento. Esseerro se dá devido ao número finito de bits utilizados para o armazenamento de umnúmero na memória do computador, causando a redução das casas decimais utilizadaspara representar cada número. Essa redução pode causar mudanças significativasnos resultados, devido a quantidade de operações realizadas, caso não seja tratadaadequadamente. (GILAT; SUBRAMANIAM, 2008)

Além deste erro, outro erro presente é o erro de truncamento que ocorre devidoas aproximações utilizadas pelo método numérico para a solução do problema. Essasaproximações são feitas para que a avaliação numérica de funções possa ser feita comum número finito de termos. A diferença entre o valor verdadeiro (exato) da funçãoe o valor aproximado obtido pelo método numérico é o erro de truncamento. A somados erros de aproximação e de truncamento é o erro total, isto é, a diferença entrea solução real do problema e a solução obtida pela utilização do método numéricocomputacionalmente. (GILAT; SUBRAMANIAM, 2008; LEVEQUE, 2007)

26

2.3 Malha de representação do domínio

Para a utilização do método das diferenças finitas, é necessária a representaçãodo domínio do problema na forma de uma malha composta por pontos finitos. Diz-seque essa malha é estruturada quando os nós estão posicionados de forma regular,estando eles igualmente espaçados ou não. Esta forma de posicionamento, apesar detrazer benefícios quanto a caracterização, ordenamento e identificação dos vizinhos decada nó, tem a desvantagem de ser limitada quanto a representação de geometriasmais complexas, o que pode ser feito facilmente com uma malha não estruturada,conforme esquematizadas nas Figuras 2a e 2b.

Figura 2 – Exemplo de malhas

(a) estruturada (b) não estruturada


2.4 Condições de contorno

Para a formulação do problema, é necessária a especificação das condições decontorno no domínio. Os tipos de condições de contorno mais comuns são:

• Variável dependente especificada na fronteira do domínio, também conhecidacomo condição de contorno de Dirichlet ou condição de contorno de primeiro tipo;

• Derivativo na direção normal especificado na fronteira do domínio, também co-nhecida como condição de contorno de Neumann ou condição de contorno desegundo tipo;

• Combinação linear das condições de contorno de primeiro e segundo tipo especi-ficada na fronteira do domínio, também conhecida como condição de contorno deRobins ou condição de contorno de terceiro tipo.

Condições de contorno onde a temperatura é constante são representadas pelacondição de contorno de primeiro tipo, enquanto condições de contorno de paredesisoladas, convecção e fluxo de calor definido nas fronteiras são descritos pelas condi-ções de segundo tipo. As condições de contorno de terceiro tipo, assim como outrostipos de condição de contorno, podem ser encontradas em problemas de transferênciade calor, mas a maioria pode ser descrita pelos dois tipos indicados anteriormente.(MINKOWYCZ, 1988; OZISIK, 1990)

27

2.5 Método das diferenças finitas e o balanço de energia

O método das diferenças finitas é o método numérico que consiste na solução deequações diferenciais por meio da substituição dos derivativos por equações algébricasque são escritas em função destes pontos. Quando resolvida dessa forma, a equaçãosó é resolvida nos nós da malha. (CENGEL; GHAJAR, 2012; OZISIK, 1990) As formasde substituição dos derivativos por equações (discretização) são descritas na seção2.6 deste mesmo capítulo.

Uma forma de montagem do sistema de equações é a realização do balançode energia, que é feito conforme demonstrado na seção 2.1. Esse método, além deser mais intuitivo, resulta no mesmo conjunto de equações obtido pelo método dasdiferenças finitas. (OZISIK, 1990; INCROPERA, 2008)

Existem várias outras formas de se obter a formulação para a solução numérica,como o método dos elementos finitos e o método dos volumes finitos, cada um com assuas vantagens e desvantagens. O método das diferenças finitas se destaca frente aosoutros por não possuir formulações matemáticas pesadas e permitir um bom entendi-mento dos problemas físicos. Esse método permite uma formulação fácil de problemasque podem ser representados por malhas estruturadas regulares sem a necessidadeda criação de elementos no domínio. (CENGEL; GHAJAR, 2012; MINKOWYCZ, 1988)

2.6 Discretização

O processo de transformação das equações diferenciais que regem o problemaem equações algébricas resolvidas pontualmente é chamado de discretização. Esseprocesso pode ser feito de diversas formas, como o método da expansão de Taylor,regressão polinomial ou método integral. (GILAT; SUBRAMANIAM, 2008; PLETCHER;TANNEHILL; ANDERSON, 1997).

O método da expansão por série de Taylor (Equação 2.3) consiste no cálculo dovalor que a função assume no ponto de interesse utilizando os valores da função e dassuas derivadas em outros nós. Como a série de Taylor é infinita, o controle da precisãoque se deseja para o valor da função é feito pelo número de termos da expansão quesão utilizados para o cálculo, além de possuir a vantagem, frente aos outros métodosde discretização, de fornecer uma estimativa do erro de truncamento.

f(xi+1) = f(xi) + f ′(xi)l +f ′′(xi)

2! l2 + f

′′′(xi)3! l

3 + f′′′′(xi)4! l

4 + ... (2.3)

Onde:l é o espaçamento entre os nós utilzados para a discretização;xi e xi+1 são pontos no domínio.

28

Os valores das derivadas nos pontos de interesse são obtidos por meio dorearranjo dos termos da expansão da série de acordo com a precisão desejada. Essasderivadas podem ser calculadas de diversas formas, seja utilizando pontos anteriores aoponto de interesse (método da diferença finita regressiva), pontos posteriores (métododa diferença finita progressiva) ou uma combinação dos pontos anteriores e posteriores(método da diferença finita central). (GILAT; SUBRAMANIAM, 2008)

Para o método das diferenças finitas central com dois pontos para a derivada pri-meira, por exemplo, temos a representação mostrada na Equação 2.4. Nessa equação,o primeiro termo do lado direito é o valor da aproximação para a derivada, enquantoque o segundo termo representa o erro de truncamento devido a aproximação, com esti-mativa mostrada na Equação 2.5, junto com a denotação que indica que a aproximaçãotem segunda ordem de precisão.

f ′(xi) =f(xi+1) + f(xi−1)

2l +O(l2) (2.4)

O(l2) = −f′′′(ξ1)3! l

2 − f′′′(ξ2)3! l

2 (2.5)

Onde:ξ1 é um valor de x entre xi e xi+1;ξ2 é um calor de x entre xi−1 e xi.

2.7 Método implícito de avanço no tempo

2.7.1 Classificação de problemas quanto ao comportamento

Problemas de governados por equações diferenciais são classificados em relaçãoao seu comportamento ao longo do tempo. Estes podem ser estacionários, comoexemplo os descritos pela equação de Laplace, problemas de autovalor e problemasde avanço no tempo. Problemas de autovalor são uma extensão dos problemas estaci-onários, a diferença está no fato de que a solução existe somente em pontos discretosno tempo. Já os problemas de avanço no tempo exigem uma solução em um domí-nio aberto em relação ao tempo e estão submetidos a condições iniciais definidas.(PLETCHER; TANNEHILL; ANDERSON, 1997)

2.7.2 Métodos de avanço no tempo e estabilidade

Entre os métodos mais difundidos de avanço no tempo estão os métodos implícitoe explícito. Esses métodos diferem em função da forma que calculam o valor da variáveldependente, temperatura nodal no caso da equação do calor. No método explícito, a

29

variável dependente é avaliada no instante de tempo anterior, enquanto que no métodoimplícito ela é avaliada no mesmo instante de tempo. O método implícito possui umavantagem sobre o método explicito de incondicionalmente estável, isto é, converge paraqualquer passo de tempo, o que não é verdade para o método explícito que está sujeitoa uma condição de estabilidade. (MINKOWYCZ, 1988; INCROPERA, 2008; CENGEL;GHAJAR, 2012)

2.8 Solução de sistemas lineares

2.8.1 Métodos iterativos de solução

Métodos iterativos para a solução de sistemas lineares foram desenvolvidos juntocom a necessidade de lidar com sistemas cada vez maiores e esparsos. Existem méto-dos mais clássicos, diretos como Jacobi e Gauss-Siedel e o método da superelaxação,que servem o propósito mas foram se tornando obsoletos devido a velocidade deconvergência ser menor do que a dos métodos iterativos. Métodos iterativos encontrama solução do sistema através de tentativa e correção e melhoram continuamente ovalor do estipulado para a solução, convergindo mais rápido ao custo de um maior usodo armazenamento do computador. (LEVEQUE, 2007; MINKOWYCZ, 1988)

2.8.2 Método do resíduo mínimo generalizado - GMRES

Apesar de muitos métodos de solução de sistemas que não são simétricos epositivamente definidos (SPD) terem sido desenvolvidos, o método do resíduo mínimogeneralizado (GMRES), é um dos métodos iterativos mais difundidos. Ele apresentavantagem sobre outros métodos de solução de sistemas lineares por utilizar o métodode iterações Arnoldi, que diminui a norma do resíduo, aumentando assim a velocidadede convergência (LEVEQUE, 2007; GOLUB; LOAN, 1996). O detalhamento do métodode solução utilizado pelo GMRES pode ser encontrada em (SAAD; SCHULTZ, 1986) eestá fora do escopo deste trabalho.

2.8.3 Pré-condicionadores

Quando se trata de sistemas lineares, pode-se quantificar a sensibilidade doproblema, isto é, o quão precisa será a solução após a aproximação. Essa quantificaçãoé feita pelo número de condicionamento, definido pela Equação (2.6) para matrizesquadradas, e depende do problema. Quando o número de condicionamento tem valorpróximo a um, diz-se que a matriz está bem condicionada, e diz-se que o problemaestá mal condicionado quando o número de condicionamento é elevado. (LEVEQUE,

30

2007; GOLUB; LOAN, 1996)

κ(A) = ‖A‖‖A-1‖ (2.6)

Onde:κ é o número de condicionamento.

O pré-condicionamento do sistema reduz o número de condicionamento, aumen-tando a velocidade da convergência. O uso pré-condicionadores é essencial para amaioria das aplicações práticas. Uma técnica de pré-condicionamento utilizada é a defazer uma fatoração utilizada para a solução de sistemas lineares, Cholesky ou lowerupper (LU), de forma incompleta na matriz do sistema (Equação 2.7), isto é, até que umnúmero determinado de elementos não nulos apareça nesta matriz. Após isso, a matrizde pré-condicionamento (Equação 2.8) é substituída no sistema original (Equação 2.9)que então é resolvido. Os métodos de pré-condicionamento são muito eficazes emmétodos de solução de sistemas que não são SPD, como o GMRES, e por isso sãomuito utilizados. (LEVEQUE, 2007)

A ≈ CTC (2.7)

M = CTC (2.8)

Au = f ⇐⇒M -1Au = M -1f (2.9)

2.9 Programação em linguagem FORTRAN

Desenvolvida na década de 50 e utilizada nos primeiros computadores, a lin-guagem FORTRAN surgiu para a facilitar a comunicação dos programadores com oscomputadores, pois, até 1954, quase toda programação era feita em linguagem demáquina. Apesar de suas primeiras versões serem antigas, a linguagem continua a serutilizada, principalmente, nos campos da ciência da computação e da análise numérica.Entre as vantagens da sua utilização estão a sua simplicidade e prioridade com arelação a velocidade de execução de programas. (PIERCE, 2002; BACKUS, 1978)

2.10 Verificação e validação de códigos

No campo de análise numérica, duas etapas podem ser utilizadas para a análiseda eficácia de um método se solução, a verificação e a validação. A verificação con-siste na comparação entre o resultado obtido com a solução exata ou outro valor de

31

referência. Essa comparação pode ser feita ao longo de uma linha ou plano, sendoque quando feita em uma linha, sua visualização fica mais fácil. Já a validação, nométodo das diferenças finitas, consiste na análise da influência do número de nós querepresentam o domínio na diferença entre o valor calculado e o valor exato, sendoa magnitude dessa influência utilizada para mensurar a força do método numéricoutilizado. (MINKOWYCZ, 1988)

2.11 Dissipadores de calor

2.11.1 Funcionamento

Para aumentar a transferência de calor por convecção entre um sistema e o meiosem mudar as suas temperaturas, pode-se utilizar de duas estratégias. A primeiraconsiste no aumento do coeficiente convectivo, que muitas vezes exige a instalaçãode um sistema para aumentar a velocidade do fluido em contato com a superfície, oua substituição do fluido do meio, soluções que não são muito práticas e por vezesevitadas. A segunda estratégia consiste no aumento da área superficial em contato como meio, por meio de superfícies estendidas, também chamadas aletas. (INCROPERA,2008; CENGEL; GHAJAR, 2012; KRAUS; AZIZ; WELTY, 2001)

Dissipadores de calor se utilizam dessa ideia, seja trabalhando sob regime deconvecção forçada ou natural, além de utilizar materiais com condutividade térmicaelevada, buscando diminuir a diferença de temperatura entre a ponta da aleta e a suabase. (INCROPERA, 2008) Dissipadores de calor são utilizados nos mais variadosequipamentos que realizam alguma forma de troca térmica, como motores de motoci-cleta e automóveis, processadores de computadores, outros componentes eletrônicos,etc.

2.11.2 Utilização em componentes eletrônicos

Os dissipadores tiveram que evoluir junto com os componentes eletrônicos, poisesses passaram a ser cada vez menores, aumentando assim a densidade de geraçãode calor e diminuindo a área superficial de troca de calor, exigindo assim coeficientesde troca térmica cada vez maiores. Temperaturas acima da recomendada de trabalhosão evitadas em componentes eletrônicos, pois aumentam a sua probabilidade de falhae o ruído de sinal devido ao aumento no movimento dos elétrons livres. (REMSBURG,2001)

Segundo Remsburg (2001), com o desenvolvimento de métodos de soluçãonumérica, a avaliação numérica de geometrias de dissipadores, que era feita somenteem laboratórios, é uma atividade comum nas empresas de desenvolvimento desses

32

componentes e pode ser feita em um computador pessoal sem grandes problemas.

2.11.3 Dissipadores de calor passivos

Um tipo específico de dissipador de calor vem encontrando grande aplicação nomercado. São estes os que se utilizam de convecção natural para dissipar o calor,devido a vantagens como a não emissão de ruido, não consumo energia para ofuncionamento de ventiladores ou outros dispositivos de bombeamento de fluido paratroca de calor e preço de fabricação reduzido.

33

3 METODOLOGIA

Neste capítulo é explicitada toda a metodologia que utilizada neste trabalho,desde a geração e pré-condicionamento da malha para a representação dos casos,discretização em diferenças finitas das equações que regem o problema, tradução dascondições de contorno para o problema, método de validação do código e casos deverificação/validação e aplicação utilizados.

3.1 Geração da malha

Conforme descrito no Capítulo 2, o método das diferenças finitas exige que sejacriada uma representação geométrica do problema por meio de uma malha constituídapor nós, onde serão calculados os valores da propriedade de interesse, no casoa temperatura. Além disso, para esse trabalho, os nós da malha necessitam estarigualmente espaçados com os seus vizinhos em todas as direções, simplificando asformulações e o pré-condicionamento da malha.

O código é feito de forma a poder gerar a malha de geometrias mais simples,como placas planas, cubos e paralelepípedos, além dos casos de aplicação a seremanalisados, ou de poder ser dada a entrada de um arquivo de malha gerado exter-namente, desde que contenha as coordenadas de todos os nós. Para a geração dasmalhas que representarão os problemas a serem simulados nesse trabalho, é utilizadoum código, criado pelo autor.

3.2 Código de pré-condicionamento da malha

Foi desenvolvido um código, presente no Apêndice B, o qual gera a malha,quando desejado, e realiza o seu pré-condicionamento seguindo as etapas indicadasno fluxograma presente na Figura 3. A descrição de como é realizada cada uma destasetapas é feita nas subseções seguintes (Subseções 3.2.1 a 3.2.5) e os procedimentospara a utilização do código encontram-se no Apêndice C.

3.2.1 Identificação dos vizinhos dos nós

Após a geração ou leitura dos nós da malha, o código identifica todos os 26vizinhos de cada nó do domínio respectivamente, conforme indicado na Figura 4. Essaverificação é feta realizando um teste lógico, onde cada posição de vizinhança estárelacionado a um conjunto de respostas desse teste. As respostas correspondentesa cada vizinhos estão presentes no arquivo ”lnd.dat”, com o arquivo completo noApêndice A.1. Este é o primeiro arquivo de entrada do código de pré-condicionamento

34

Figura 3 – Esquemático do código de geração e pré-condicionamento da malha.

Entrada dos parâmetros ou do arquivo com a malha

Geração ou leiturada malha

Espelhamento dos vizinhos dos nós nas fronteiras

adiabáticas

Saída do arquivo de malha com as coordenadas, classificações geométrica e física, vizinhança e faces e

volume de cada nó

Leitura dos arquivos de entrada

Classificação dos nós quanto ao número de vizinhos que possui

Identificação dos vizinhos de cada nó

Classificação geométricade cada nó

Classificação física de cada nó e aplicação das

condições de contorno

Arquivos de entrada da malha

Có

dig

o Cód

igo


da malha, que não precisa ser editado para a utilização do código.No arquivo ”lnd.dat”, as colunas 1, 4 e 7 correspondem ao teste lógico de que o

nó vizinho possui a mesma coordenada x, y ou z, respectivamente. Já as colunas 2,5 e 8, indica se o nó está a uma unidade de distancia no sentido positivo do nó nasdireções x, y ou z, também respectivamente, e os nós 3, 6 e 9 nas direções negativas.

3.2.2 Classificação geométrica dos nós

Após a identificação dos vizinhos, o programa verifica qual a classificação geomé-trica do nó, isto é, se é um nó de interior, de face, aresta ou vértice (interior ou exterior).Essa verificação é feita por meio de outro teste lógico, que, conforme o número eposição dos vizinhos que o nó possui, comparando com as respostas pré-programadas.Ao todo, o programa identifica 41 posições geométricas possíveis, conforme indicadona Figura 5.

As respostas do teste lógico correspondentes a cada posição estão indicadasno arquivo ”lke.dat”, de entrada do código, e, devido a algumas posições estaremrelacionadas a mais de uma resposta, o arquivo possui varias linhas que levam amesma posição geométrica. O arquivo completo encontra-se no Apêndice A.2, sendo

35

Figura 4 – Identificação da posição dos 26 vizinhos de cada nó

18

6

16

22

12

20

26

14

24

5

1

4

9

2

8

11

3

10

19

7

17

23

13

21

27

15

25

y

xz


que este arquivo também não precisa ser editado para a utilização do código depré-condicionamento da malha.

3.2.3 Classificação física dos nós

Após a classificação geométrica dos nós, o código realiza a classificação físicados mesmos, isto é, indica qual o comportamento físico do nó durante a análisenumérica, sendo que os nós podem ter temperatura calculada ou não (prescrita), ouestar sujeito a isolamento, convecção ou fluxo de calor prescrito. Para cada casoa ser pré-condicionado, deve ser criado um arquivo de entrada, o arquivo ”bcf.dat”,contendo a classificação física relacionada a cada classificação geométrica. Os arquivoscorrespondentes a cada caso simulado, incluindo as aplicações, estão presentes noApêndice A.4 e serão discutidos posteriormente.

Neste arquivo, a terceira coluna corresponde a classificação física correspon-dente a classificação geométrica indicada na primeira coluna. Quando uma posiçãogeométrica possui mais de uma classificação física, as coordenadas x, y e z, mínimase máximas, onde a classificação física se aplica, são indicadas nas colunas 4 a 9.A segunda coluna indica se a posição geométrica está sujeita a uma condição decontorno com temperatura prescrita (valor de 1 a 9), convecção (valor de 11 a 19) ou

36

Figura 5 – Identificação das 41 posições geométricas possíveis para um nó

27

15

33

34

29

28

37

30

32

35

38

3940

31

36

41

11

3

26

14

19

7

5

1

18

6

23

13

9

2

22

12

y

xz

25

4

24 21

2016 8

10 17


fluxo de calor prescrito (valor de 101 a 109), se o nó é de interior ou é uma fronteiraadiabática (por onde não há transferência de calor) (coluna com valor 0).

As propriedades correspondentes a cada classificação física de nó estão presen-tes no arquivo ”bdf.dat”, quarto e último arquivo de entrada do código, que também nãoprecisa ser editado. O arquivo completo está no Apêndice A.3.

Neste arquivo, a primeira coluna indica qual o volume relacionado a classificaçãofísica (multiplicado por 8), o número de faces (área multiplicada por 4) que estão sujeitosa condução, convecção e fluxo de calor prescrito, nas direções x, y e z, positivas enegativas, nas colunas 2 a 19. Nas colunas 20 a 57 são indicados, em pares, os nósque devem ser espelhados, quando há alguma face adiabática, que será explicadoposteriormente neste trabalho.

37

3.2.4 Espelhamento dos nós de fronteiras adiabáticas

Caso o nó tenha uma ou mais faces adiabáticas, o código faz, conforme indicadono arquivo ”bdf.dat”, o espelhamento dos nós. Esse espelhamento se dá pela substitui-ção do valor da temperatura em um nó pelo valor da temperatura em um nó que deveter a mesma temperatura, devido ao espelhamento.

3.2.5 Classificação dos nós quanto aos vizinhos que possui

O ultimo passo realizado pelo código de pré-condicionamento da malha é averificação do número de vizinhos que cada nó possui, seja ele espelhado ou não. Essaverificação permite que o código verifique qual discretização pode ser utilizada em cadanó. O nó pode ser classificado como 0 (quando não possui todos os 26 vizinhos) ou 1(quando possui).

Após realizadas todas essas etapas, o código imprime um arquivo de saída,”mesh.dat”, contendo todas as informações correspondentes a cada nó do domínio,incluindo coordenadas (colunas 1 a 3), classificações geométrica, de condição decontorno, física, e de número de vizinhos (colunas 4 a 7), todos os vizinhos que o nópossui (colunas 8 a 33), incluindo as posições na qual ele não possui vizinhos, volume(coluna 34) e número de faces correspondentes a cada condição de condição decontorno (colunas 35 à 53). Os arquivos ”mesh.dat” não estão presentes nesse trabalhodevido aos seus tamanhos tornarem isso inviável. Este arquivo serve de entrada para ocódigo principal de análise numérica, descrito na Seção 3.5.

3.3 Discretização das equações físicas que regem o problema

São utilizadas duas equações equivalentes para descrever fisicamente os pro-blemas a serem analisados, a Equação 3.1 e a Equação 3.3. A primeira equação éutilizada nos nós de interior ou sujeito somente a condições de contorno de isolamento(classificados como 1 na classificação de vizinhança). Já a segunda equação é utilizadanos demais nós.

∇2T = 1α

∂T

∂t− q̇k

(3.1)

Com:

∇2T = ∂2T

∂x2+ ∂

2T

∂y2+ ∂

2T

∂z2(3.2)

ΣĖcond + ΣĖconv + ΣĖflux + Ėg = V– ρcp∂T

∂t(3.3)

38

Onde:α é a difusividade térmica do material;cp é o calor específico a pressão constante;Ėcond é o fluxo de calor devido a condução;Ėconv é o fluxo de calor devido a convecção;Ėflux é o fluxo de calor prescrito;Ėg é a energia total gerada no volume de controle;q̇ é o fluxo de calor por unidade de área;V– é o volume do domínio analisado.

Na Equação 3.3, o somatório corresponde aos calores que ultrapassam a fronteirado nó, sejam eles de condução, convecção ou fluxo de calor prescrito. Quando escritosem função de derivativos de primeira ordem (condução e convecção), são utilizadas asdiscretizações descritas na Subseção 3.3.2.

3.3.1 Discretização espacial dos nós de interior

Segundo O’Reilly e Beck (2006), o laplaciano presente nas Equações 3.1 e 3.2pode ser escrito em função de 7, 19 e 27 nós, incluindo o nó de interesse. Conforme onúmero de nós utilizados, obtém-se uma ordem de convergência diferente, conformeserá estudado no Capítulo 4. As equações referentes a cada laplaciano (Equações 3.4a 3.6 respectivamente) e a representação dos nós levados em consideração no cálculo(Figura 6) são indicadas á seguir.

Figura 6 – Esquemático dos nós segundo as suas posições

Posições

1 (faces)

2 (arestas)

3 (vértices)

x

y

z

Fonte: adaptado de O’Reilly e Beck (2006).

∇2Ti =1l2∑j∈Nf

(Tj − Ti) (3.4)

∇2Ti =1

6l2

2 ∑j∈Nf

Tj +∑j∈Na

Tj − 24Ti

(3.5)

39

∇2Ti =3

13l2

∑j∈Nf

Tj +12∑j∈Na

Tj +13∑j∈Nv

Tj −443 Ti

(3.6)

Onde:m e n são os índices referentes aos nós;Na são os nós localizados nas arestas dos elementos;Nf são os nós localizados nas faces dos elementos;Nv são os nós localizados nos vértices dos elementos.

3.3.2 Discretização espacial dos nós de fronteira

Nos nós que não possuem todos os vizinhos necessários para o cálculo dolaplaciano, é feito um balanço de energia, conforme indicado na Equação 3.3. O códigopossui implementados os derivativos progressivos de 2 (Equação 3.7) e 3 nós (Equação3.8), sendo que a primeira opção só pode ser utilizada nos nós que possuem pelomenos os 6 vizinhos posicionados nas faces.

∂Tn∂r

= Tn+1 − Tnl

(3.7)

∂Tn∂r

= 3Tn+1 − 4Tn + Tn−12l (3.8)

Nas Equações 3.7 e 3.8, r é a direção na qual está sendo calculado o derivativo en é o índice correspondente a direção, logo, r pode assumir x, y ou z, com n assumindoi, j ou k, respectivamente.

3.3.3 Discretização temporal implícita de primeira ordem de precisão

A discretização temporal, de forma implícita com primeira ordem de precisão, éfeita por meio da substituição (Equação 3.9) do derivativo em relação ao tempo naequação obtida pela discretização espacial, seguida pela adição do índice p+ 1 nosdemais termos, exceto o termo referente a geração interna de energia.

∂T

∂t= T

p+1 − T p

∆t (3.9)

Onde:∆t é o passo de tempo.

40

3.4 Montagem da matriz

Após todas as substituições, e adequações para os dados presentes no arquivo”mesh.dat” gerado, obtém-se as equações finais de cálculo da temperatura nodal, emfunção das propriedades do problemas. Sendo elas:

Nós de interior de casos estacionários e transientes utilizando o laplaciano de 7 nós:

(6 + γ l

2

α∆t

)T p+1i,j,k − (T

p+1i+1,j,k + T

p+1i−1,j,k + T

p+1i,j+1,k + T

p+1i,j−1,k + T

p+1i,j,k+1 + T

p+1i,j,k−1)

= l2egk

V– i,j,k +γl2

α∆tTpi,j,k

(3.10)


(4 + γ l

2

α∆t

)T p+1i,j,k

−13(Tp+1i+1,j,k + T

p+1i−1,j,k + T

p+1i,j+1,k + T

p+1i,j−1,k + T

p+1i,j,k+1 + T

p+1i,j,k−1)

−16(Tp+1i+1,j+1,k + T

p+1i+1,j−1,k + T

p+1i−1,j+1,k + T

p+1i−1,j−1,k + T

p+1i+1,j,k+1 + T

p+1i+1,j,k−1

+T p+1i−1,j,k+1 + Tp+1i−1,j,k−1 + T

p+1i,j+1,k+1 + T

p+1i,j+1,k−1 + T

p+1i,j−1,k+1 + T

p+1i,j−1,k−1)

= l2egk

V– i,j,k +γl2

α∆tTpi,j,k

(3.11)


(4413 + γ

l2

α∆t

)T p+1i,j,k

− 313(Tp+1i+1,j,k + T

p+1i−1,j,k + T

p+1i,j+1,k + T

p+1i,j−1,k + T

p+1i,j,k+1 + T

p+1i,j,k−1)

− 326(Tp+1i+1,j+1,k + T

p+1i+1,j−1,k + T

p+1i−1,j+1,k + T

p+1i−1,j−1,k + T

p+1i+1,j,k+1 + T

p+1i+1,j,k−1

+T p+1i−1,j,k+1 + Tp+1i−1,j,k−1 + T

p+1i,j+1,k+1 + T

p+1i,j+1,k−1 + T

p+1i,j−1,k+1 + T

p+1i,j−1,k−1)

− 113(Tp+1i+1,j+1,k+1 + T

p+1i+1,j+1,k−1 + T

p+1i+1,j−1,k+1 + T

p+1i+1,j−1,k−1

+T p+1i−1,j+1,k+1 + Tp+1i−1,j+1,k−1 + T

p+1i−1,j−1,k+1 + T

p+1i−1,j−1,k−1)

= l2egk

V– i,j,k +γl2

α∆tTpi,j,k

(3.12)

41

Nós de fronteira de casos estacionários e transientes utilizando a derivada progressivade 2 nós:

(γl2

α∆t + A +h∞l

kB)T p+1i,j,k

−ax+T p+1i+1,j,k − ax−Tp+1i−1,j,k − ay+T

p+1i,j+1,k − ay−T

p+1i,j−1,k − az+T

p+1i,j,k+1 − az−T

p+1i,j,k−1

= l2egk

V– i,j,k +γl2

α∆tTpi,j,k +

h∞lT∞k

B + lQ∞k

C

(3.13)

Nós de fronteira de casos estacionários e transientes utilizando a derivada progressivade 3 nós:

(γl2

α∆t + 2A +h∞l

kB)T p+1i,j,k −

3ax+ + ax−2 T

p+1i+1,j,k −

3ax− + ax+2 T

p+1i−1,j,k

−3ay+ + ay−2 Tp+1i,j+1,k −

3ay− + ay+2 T

p+1i,j−1,k

−3az+ + az−2 Tp+1i,j,k+1 −

3az− + az+2 T

p+1i,j,k−1

= l2egk

V– i,j,k +γl2

α∆tTpi,j,k +

h∞lT∞k

B + lQ∞k

C

(3.14)

Sendo que, nas Equações 3.13 e 3.14, temos os seguintes valores para A, B e C,onde a. b e c são constantes indicando a quantidade de faces sob condução, convecçãoe fluxo de calor prescrito, respectivamente, nas direções indicadas em seus subíndices.

A = ax+ + ax− + ay+ + ay− + az+ + az− (3.15)

B = bx+ + bx− + by+ + by− + bz+ + bz− (3.16)

C = cx+ + cx− + cy+ + cy− + cz+ + cz− (3.17)

Nos casos estacionários, o valor de γ é igual a 0, e o mesmo é igual a 1, noscasos transientes, nas Equações 3.10 à 3.14.

Após a seleção das equações a serem utilizadas em cada caso, o código principal(com funcionamento descrito na próxima seção) as utiliza para montar o sistema linearglobal a ser resolvida uma vez, caso seja um problema estacionário, ou a acada passode tempo, caso seja um problema transiente, no formato do sistema linear 3.18

42

α11 α12 α13 . . . α1n

α21 α22 α23 . . . α2n

α31 α32 α33 . . . α3n...

...... . . .

...αn1 αn2 αn3 . . . αnn

×

T1

T2

T3...Tn

=

β1

β2

β3...βn

(3.18)

Onde:α e β são valores reais e definidos;T é a temperatura nodal a ser calculada.

3.5 Código utilizado na simulação

Para resolver a distribuição de temperaturas nos problemas, foi feita a adequaçãode um código já existente, com o método de pré-condicionamento e solução do sistemalinear já implementados (CONCI; GOMES, 2017), com funcionamento descrito naFigura 7.

Figura 7 – Esquemático do código de montagem, pré-condicionamento e solução do sistema linear.

Critérios cumpridos?

Caso transiente?


Arquivo de entrada com nós e suas informações

Arquivos de saída com as informações de interesse, como distribuições de temperatura e informações das iterações

Leitura dos dados de entrada e declaração das

variáveis

Critério de convergência e/ou passo de tempo

Comparação com a solução exata (caso tenha) e

impressão das saídas

Aplicação das condições iniciais e de contorno

Montagem do sistema global a ser resolvido

Pré-condicionamento (caso aplicável) e solução do

sistema linear

Parâmetros a serem utilizados na simulação

Có

dig

o Cód

igo


Não

- A

vanç

o no

tem

po

Sim

Sim

Não


43

Primeiro, são utilizados como entradas do programa, o arquivo ”mesh.dat” geradodo código de pré-condicionamento da malha e os parâmetros referentes aos problemasa serem analisados. Entre os parâmetros a serem explicitados estão o tipo de análise aser feita (estacionária ou transiente), número de passos de tempo a serem utilizados ouerro de convergência (caso seja uma análise transiente), o tipo de problema que seráestudado (bi ou tridimensional), a ordem de precisão a ser utilizada nos derivativosespaciais, se é um caso de verificação/validação e qual, se a solução exata deve sercalculada, e os parâmetros físicos da simulação, como constantes de condutividade edifusividade térmica, passo de tempo, taxa de geração de energia interna, temperaturas,coeficientes convectivos e fluxos de calor prescritos a serem utilizados nas condiçõesde contorno.

Após a entrada das informações, compilação e execução do código, é feitaa aplicação das condições de contorno nos nós da malha correspondentes e dascondições iniciais do problema. Então, é feita a montagem do sistema linear global aser resolvido para a solução do problema, utilizando os critérios especificados, seguidapelo pré-condicionamento (utilizando o método MILU0) e solução do sistema linear(utilizando o GMRES).

Caso seja um problema transiente, é verificado se o critério de convergência oupasso de tempo máximo foi alcançado, caso não tenha sido, o código avança no tempoe retorna para o passo de montagem do sistema linear global. Após o atendimento doscritérios de convergência (para o caso transiente) ou término da solução (para o casoestacionário), é feita a solução e comparação com a solução exata do problema (casopossua), seguidos da impressão dos arquivos de saída do código.

η =|‖ T pq ‖ − ‖ T p−1q ‖|

‖ T pq ‖< � (3.19)

Onde:N é o número total de nós na malha.

Um dos critérios que podem ser utilizados para o caso transiente é a especificaçãode um valor de � para a Equação 3.19, de forma que, uma vez que esta equação sejasatisfeita, o código para de avançar no tempo e imprime a solução utilizada no problema.Essa forma de controle de convergência só interessante em casos que a solução tenhauma assíntota em relação ao tempo, caso contrário, deve-se utilizar o método decontrole pelo número de passos de tempo máximo.

44

3.6 Verificação e validação do código

Para a verificação e validação do código principal pode ser feita a simulação deproblemas para os quais se tenha a solução exata. Para isso, foram selecionados seiscasos, cinco estacionários e um transientes, a serem utilizados como referência para acomparação de resultados.

Para a verificação, é feita a plotagem dos resultados obtidos via simulação e viasolução analítica em uma parte do domínio (uma linha neste trabalho), podendo-seentão verificar e comparar as soluções encontradas.

Já para a verificação, quanto a ordem de convergência, são simuladas três malhaspara cada caso, podendo-se assim plotar um gráfico, em escala logarítmica, do erroem função do distanciamento dos nós, e calcular a tangente da(s) retas resultantes,comparando com os resultados esperados. Este erro é calculado utilizando as normasL1, L2 e L∞, utilizando as Equações 3.20 à 3.22, respectivamente.

L1 =N∑q=1

|Tq − Tq,exata|N

(3.20)

L2 =

√√√√ N∑q=1

|Tq − Tq,exata|2N

(3.21)

L∞ = máx(|Tq − Tq,exata|) ∀ 1 ≤ q ≤ N (3.22)

Apesar da equivalência, as três normas são utilizadas devido a possibilidade dealguns problemas não convergirem para todas as normas, aumentando a quantidadede parâmetros que podem ser acompanhados nesses casos.

3.7 Casos a serem estudados

3.7.1 Casos estacionários com solução exata

O primeiro caso a ser estudado é uma placa plana quadrada posicionada noplano XY (com o lado da placa L = 1 m), com espessura mínima e sem fluxo de calorperpendicular a ela em suas fronteiras (direção z). A placa está sujeita as condiçõesde contorno expressas nas Equações 3.23 e 3.24. Sua solução exata é expressapela Equação 3.25. Este caso será analisado utilizando os laplacianos de 7, 19 e27 nós, para a validação das suas ordens de convergência. Já as fronteiras serãoanalisadas utilizando os derivativos progressivos com 3 nós, que há somente condiçõesde contorno de temperatura prescrita e isolamento.

T (0, y, z) = T (L, y, z) = T (x, 0, z) = 0 (3.23)

45

T (x, L, z) = sen(πx

L

)(3.24)

T (x, y, z) = senh(πy/L)senh(π) sen(πx/L) (3.25)

O segundo caso a ser estudado é uma placa retangular, posicionada no planoXZ (com a dimensão em XY b = 0,25 m), com espessura mínima e, novamente, semfluxo de calor perpendicular a ela em suas fronteiras (direção y). A placa está sujeitaas condições de contorno expressas nas Equações 3.26 a 3.28, com k = 1 W/m-k, T∞= 300 K, h∞ = 10 W/m2/K. Sua solução, segundo INCROPERA (2008), é tabelada eutilizada posteriormente nesse trabalho. Este caso utilizará o mesmo laplaciano docaso anterior. Os derivativos progressivos de dois pontos serão utilizados nos nós dafronteira, devido a condição de contorno convectiva presente em uma das faces.

T (0, y, z) = T (x, y, 4b) = 500 (3.26)

∂T (2b, y, z)∂x

= 0 (3.27)

q̇z(x, y, 0) = h∞A(T − T∞) (3.28)

Onde:q̇z é o fluxo de calor por unidade de área na direção z;h∞ é o coeficiente convectivo do ambiente;T∞ é a temperatura do ambiente;A é a área de transferência de calor.

O terceiro caso a ser estudado é um paralelepípedo de base quadrada apontandona direção z (com b = 0,5 m).. A geometria está sujeita as condições de contornoexpressas nas Equações 3.29 e 3.30 e a sua solução exata é expressa pela Equação3.31. Este caso será simulado utilizando somente o laplaciano de ordem intermediária,devido ao alto custo computacional, e os derivativos de três nós na fronteira.

∂T (b, y, z)∂x

= ∂T (x, b, z)∂y

= T (0, y, z) = T (x, 0, z) = T (x, y, 0) = 0 (3.29)

T (x, y, 2b) = 1 (3.30)

T (x, y, z) =∞∑m=1

∞∑n=1

16π2

sen(πx2m+1

L

)2m+ 1

sen(πy 2n+1

L

)2n+ 1

senh(πzL

√(2m+ 1)2 + (2n+ 1)2

)senh(π

√(2m+ 1)2 + (2n+ 1)2)

(3.31)

46

O quarto e ultimo caso estacionário é um cubo (com lado L = 1 m). A geometriaestá sujeita a condição de contorno expressa na Equação 3.32, a geração de energiainterna expressa na Equação 3.33 e a sua solução exata é expressa pela Equação3.34. Este caso será simulado utilizando somente o laplaciano de ordem intermediária,devido ao alto custo computacional, e os derivativos de três nós na fronteira.

∂T (L, y, z)∂x

= ∂T (x, L, z)∂y

= ∂T (x, y, L)∂z

= T (0, y, z) = T (x, 0, z) = T (x, y, 0) = 0(3.32)

eg(x, y, z) = −π2

12L2

{sen

(πx

2L

)sen

(πy

2L

)sen

(πz

2L

)}(3.33)

T (x, y, z) = sen(πx

2L

)sen

(πy

2L

)sen

(πz

2L

)(3.34)

3.7.2 Caso transiente com solução exata

O quinto caso, transiente, é um cubo (também com lado L = 1 m), sem fluxo decalor nas fronteiras na direção z. A geometria está sujeita as condições de contornoexpressas na Equação 3.35, condição inicial de temperatura uniforme T (x, y, z, 0) = 1, ea sua solução exata é expressa pela Equação 3.36. Este caso será simulado utilizandosomente o laplaciano de ordem intermediária, devido ao alto custo computacional, e osderivativos de três nós na fronteira.

T (0, y, z, t) = T (L, y, z, t) = T (x, 0, z, t) = T (x, L, z, t) = 0 (3.35)

T (x, y, z, t)

=∞∑m=1

∞∑n=1

4mnπ2

{[(−1)m − 1][(−1)n − 1]}sen(nπx

L

)sen

(mπy

L

)exp

(− kπ

2t

L2(m2 + n2)

)(3.36)

Os casos 1, 3, 4 e 5 são simulados com 4 malhas, dobrando a quantidade de nósde uma para outra na medida que a simulação ainda seja possível de ser realizada,diferente do caso 2 no qual são utilizadas somente 3 malhas. Nos casos em que odomínio é uma placa, simulando um problema bidimensional, e que a transferência decalor seja desprezível na direção perpendicular a placa, o número de nós utilizadosna espessura da placa é mantido constante, mesmo que isso implique na redução daespessura da placa com o refinamento da malha.

47

3.7.3 Aplicações

Após a verificação e validação do código, é feita uma aplicação prática para a suautilização. Para isso, três geometrias diferentes de dissipadores de calor passivos sãosimuladas. Para que seus resultados possam ser comparados entre si, estes possuemdimensões gerais semelhantes, e estão sujeitos a condições de contorno similares,com parâmetros indicados no Quadro 1. É utilizado o critério de parada por erro, com �= 10−6.

Figura 8 – Representação dos dissipadores de calor a serem analisados nesse trabalho

(a) Primeiro modelo (b) Segundo modelo

(c) Terceiro modelo


Quadro 1 – Parâmetros utilizados na simulação das geometrias

Parâmetro Valor UnidadeCoeficiente convecctivo 27 x 10−4 W/mm2-KCondutividade térmica 2,04 W/mm-KDifusividade térmica 9,88 x 10−7 mm2/s

Fluxo de calor na face inferior 0,42 W/mm2

Temperatura inicial 300 KTemperatura ambiente 300 K

Passo de tempo 10 sFonte: autoria própria.

48

Os dissipadores simulados são encontrados para venda na internet e são utili-zados para dissipar o calor gerado nos processadores de computadores RhaspberryPi, computadores de baixo custo com geometrias reduzidas e sem uma mecânica dedissipação forçada de calor.

O primeiro caso corresponde a um dissipador de calor com aletas finas contínuascom a sua geometria é ilustrada na Figura 8a. O segundo é constituído de aletasquadradas, de espessura maior que as da geometria anterior, ilustrado pela Figura 8b.A terceira geometria a ser simulada, com aletas seccionadas e de formato retangular, éilustrada na Figura 8c.

Figura 9 – Condição de contorno nas faces da geometria do caso 1

xz

y

Fluxo de calor prescritoIsoladaConvecção

C.C. NAS FACES

z x

y



xz

y

z x

y


C.C. NAS FACES


Os dissipadores possuem como dimensões base teórica de 14 x 14 x 5 mm,mas as suas alturas são maiores devido a dependência com relação ao espaçamentoentre os nós e aos mesmos serem equidistantes. O primeiro caso contem aletas comespessura igual a metade do espaçamento entre aletas, com espessura igual a da base.

49


xz

y

z x

y


C.C. NAS FACES


Já os segundo e terceiro têm aletas com espessura igual o espaçamento entre aletas,também igual a espessura da base. A Tabela 1 contém as características geométricasde cada modelo de dissipador de calor e o Apêndice D possui o detalhamento dasgeometrias utilizadas.

Tabela 1 – Características geométricas dos dissipadores de calor a serem estudados, onde N é onúmero de aletas, ht é a altura total da geometria, ha é a altura das aletas, A1 é a áreaprimitiva de transferência de calor por convecção, A2 é a área submetida ao fluxo de calorprescrito, A3 é a área estendida devido as aletas e A4 é a área total da geometria que estásubmetida a transferência de calor por convecção

Geometria N ht(mm) ha(mm) A1(mm2) A2(mm2) A3(mm2) A4(mm2)1 3 5,25 4,38 20,00 6,25 195,2 215,22 9 5,09 3,82 19,57 12,25 174,9 194,53 6 5,09 3,82 19,57 12,25 145,8 165,4


Para reduzir o custo computacional da análise, e devido as geometrias possuírem2 eixos de simetria, será analisado apenas um quarto das geometrias. Todos osdissipadores estão com um fluxo de calor de 5 Watts em sua base, sendo que 2 mmentre a lateral e a base não absorvem calor e que, as demais faces do dissipador estãosujeitas a convecção.

Esses dados foram obtidos junto a fabricantes dos dissipadores e um estudo maisdetalhado das condições de contorno as quais os dissipadores estão sujeitos está forado escopo desse trabalho. As Figuras 9 a 11 esquematizam as geometrias finais aserem simuladas e as condições de contorno em cada face ou porção de face nasgeometrias.

51

4 VERIFICAÇÃO E VALIDAÇÃO DO CÓDIGO

4.1 Casos estacionários

4.1.1 Caso 1 - Bidimensional com solução exata

Para a simulação do primeiro caso, são geradas quatro malhas, representadasnas Figuras 12a a 12d. O refinamento é feito conforme descrito na Seção 3.7. Comisso, é obtida a solução aproximada do problema, mostrada na Figura 13.

Figura 12 – Malhas utilizadas no caso 1

(a) 147 nós (L/6) (b) 507 nós (L/12)

(c) 1875 nós (L/24) (d) 7203 nós (L/48)


Os geometria é simulada uma vez para cada Laplaciano (de 7, 19 e 27 nós),obtendo-se os resultados para as normas mostrados nas Tabelas 2, 3 e 4, representa-das graficamente pelos gráficos nas Figuras 14a a 14c, com a diferença entre a soluçãoexata e a aproximada de cada Laplaciano, para comparação, nas Figuras 15a a 15c.

Nas Tabelas 2, 3 e 4, a notação p(L) denota o calculo do coeficiente angularobtido por meio da regressão linear utilizando os valores da linha (malha) atual e os damalha anterior, impossibilitando o cálculo do coeficiente angular para a primeira linha.

Com os valores das normas das Tabelas 2, 3 e 4, pode-ser verificar que, quandoutilizados os Laplacianos de 7 e 27 nós, o código possui ordem de convergência 2, jáquando utilizado o Laplaciano de 19 nós, o programa passa a ter ordem de convergência

52

Figura 13 – Solução aproximada do caso 1


Tabela 2 – Valores das normas e coeficientes angulares para o Laplaciano de 7 nós

Espaçamento entre nós L1 p(L1) L2 p(L2) L∞ p(L∞)L/6 -2,6812 - -2,4931 - -2,1148 -L/12 -3,1909 1,6934 -3,0540 1,8635 -2,7090 1,9741L/24 -3,7518 1,8635 -3,6370 1,9369 -3,3060 1,9834L/48 -4,3346 1,9362 -4,2298 1,9694 -3,9075 1,9983



Espaçamento entre nós L1 p(L1) L2 p(L2) L∞ p(L∞)L/6 -6,4976 - -6,3091 - -5,9305 -

L/12 -8,2189 5,7186 -8,0820 5,8900 -7,7368 6,0010L/24 -9,9858 5,8701 -9,8710 5,9435 -9,5399 5,9904L/48 -11,771 5,9309 -11,667 5,9668 -11,3449 5,9967



Espaçamento entre nós L1 p(L1) L2 p(L2) L∞ p(L∞)L/6 -2,8759 - -2,6871 - -2,3082 -L/12 -3,3976 1,7332 -3,2606 1,9053 -2,9154 2,0173L/24 -3,9616 1,8738 -3,8468 1,9475 -3,5157 1,9944L/48 -4,5452 1,9389 -4,4404 1,9721 -4,1181 2,0013


53

Figura 14 – Plotagem das normas em função do distanciamento dos nós para o caso 1 para os respecti-vos Laplacianos

(a) Laplaciano de 7 nós

●

●

●

●

−2.0 −1.5 −1.0 −0.5

−4.

0−

3.5

−3.

0−

2.5

log(l)

log(

erro

)

●

Normas

L1L2L∞

(b) Laplaciano de 19 nós

●

●

●

●

−4 −3 −2 −1 0 1

−12

−11

−10

−9

−8

−7

−6

log(l)lo

g(er

ro)

●

Normas

L1L2L∞

(c) Laplaciano de 27 nós

●

●

●

●

−2.0 −1.5 −1.0 −0.5

−4.

0−

3.5

−3.

0−

2.5

log(l)

log(

erro

)

●

Normas

L1L2L∞


6. Para esse problema, o programa apresentou grande eficácia, principalmente para oLaplaciano de 19 nós, que alcançou um erro em relação a solução exata de ordem -11.Nas Figuras 15a a 15c é plotada a diferença entre a solução aproximada e a soluçãoexata para cada Laplaciano, em uma linha, podendo-se verificar os efeitos da ordemde convergência elevada, principalmente em relação ao Laplaciano de 19 nós.

4.1.2 Caso 2 - Bidimensional com solução tabelada

Para a simulação do segundo caso, são geradas 3 malhas, representadas pelasFiguras 16a a 16c, e obtidos resultados que são comparados com os valores tabelados

54

Figura 15 – Diferença entre a solução aproximada e a exata, para cada malha, utilizando os três Laplaci-anos na linha y = 0,5 Ly e z = 0,5 Lz

(a) Laplaciano de 7 nós

●

●

●

●

●

●

●

0.0 0.2 0.4 0.6 0.8 1.0

0.00

00.

002

0.00

40.

006

x

|T−

Tex

ata|

●

Malhas

L/6L/12L/24L/48

(b) Laplaciano de 19 nós

●

●

●

●

●

●

●

0.0 0.2 0.4 0.6 0.8 1.00e

+00

2e−

074e

−07

6e−

078e

−07

1e−

06

x

|T−

Tex

ata|

●

Malhas

L/6L/12L/24L/48

(c) Laplaciano de 27 nós

●

●

●

●

●

●

●

0.0 0.2 0.4 0.6 0.8 1.0

0.00

00.

002

0.00

40.

006

x

|T−

Tex

ata|

●

Malhas

L/6L/12L/24L/48


de referência1, mostrados na tabela 5. A distribuição de temperatura obtida é mostradana Figura 17.

Na Tabela 5 os índices (1 a 8) correspondem as posições nas quais foramcalculadas as temperaturas analiticamente e estão representadas em INCROPERA(2008) pág. 141.

Para os valores calculados, tem-se que a maior discrepância de valores em um nóé de 2,41%, no nó 7, verificando a acurácia do método de solução para um problemacom condição de contorno convectiva em uma das faces.

1 Valores de referência coletados de INCROPERA (2008)

55


(a) 693 nós (b) 2583 nós (c) 5673 nós




Tabela 5 – Comparação entre os valores calculados e os valores de referência

N. de nós T1 T2 T3 T4 T5 T6 T7 T8Ref. 489,30 485,15 472,07 462,01 436,95 418,74 356,99 339,05693 489,71 485,52 472,54 461,91 435,65 416,31 348,53 337,28

2583 489,69 485,49 472,47 461,83 435,47 416,19 348,39 337,275673 489,68 485,47 472,45 461,81 435,41 416,17 348,38 337,28


4.1.3 Caso 3 - Tridimensional com solução exata

São geradas 4 malhas para a simulação do terceiro caso, representadas nasFiguras 18a a 18d. É obtida a solução aproximada, exibida na Figura 19.

A diferença entre a solução exata e a solução aproximada encontrada para asdiferentes malhas utilizadas, em uma determinada linha, é mostrada na figura 20. Nesta

56


(a) 45 nós (L/2) (b) 225 nós (L/4)

(c) 1377 nós (L/8) (d) 6525 nós (L/14)




figura também é possível perceber a expressiva redução na diferença com a soluçãoexata a cada refinamento da malha.

57

Outro fato importante de ressaltar é a diferença entre a solução aproximada ea exata ter aumentado consideravelmente quando comparada com a obtida no caso1. Isso se deve ao fato de a solução exata não ser calculada diretamente (Equação3.31) e sim por meio de um somatório, pois, devido a limitação de casas decimais docomputador, tanto em números positivos para os termos em seno hiperbólico, quantopara os termos negativos para a soma de frações, que, conforme o número de termosdo somatório levados em conta aumenta, chegando na ordem de 30 milhões para cadanó da malha, a combinação dos arredondamentos causa certa instabilidade na soluçãoexata, tornando mais imprecisa a sua comparação, devido a limitações da máquina.

Figura 20 – A diferença entre a solução exata e a solução aproximada encontrada para cada malha docaso 4 na linha x = 0,25, y = 0,25

●

●

●

●

●

0.0 0.2 0.4 0.6 0.8 1.0

0.00

0.01

0.02

0.03

0.04

0.05

z

|T−

Tex

ata|

●

Malhas

L/2L/4L/8L/14


4.1.4 Caso 4 - Tridimensional com solução exata

O quarto caso é analisado por meio de 4 malhas, representadas pelas Figuras21a a 21d, que, após simuladas, levam aos resultados mostrados na Figura 22

A diferença entre a solução exata e a solução aproximada, para cada malha, éplotada na Figura 23. Novamente é possível verificar a expressiva redução da diferençaentre a solução aproximada e a solução exata conforme ocorre o refinamento da malha.

4.2 Caso transiente

4.2.1 Caso 5 - Bidimensional com solução exata

O quinto caso é representado pelas mesmas malhas utilizadas no caso 1 (Figuras12a a 12d), que, após simuladas, levam aos resultados, nos instantes de tempo de0,001, 0,01, 0,05 e 0,1 segundos, representados nas Figuras 26a a 26d.

58


(a) 27 nós (h/2) (b) 125 nós (h/4)

(c) 729 nós (h/8) (d) 4913 nós (h/16)




Então, são feitas comparações, nos três instantes de tempo, dos resultadosobtidos com os da solução exata, em uma linha, mostradas nas Figuras 24a e 24be também nas Figuras 25a e 25b. Nelas, é possível verificar que, embora o códigoinicie com certa discrepância nos primeiros instantes, sendo ela menor quanto menor o

59

Figura 23 – Diferença entre a solução exata e a solução aproximada, para cada malha do caso 4, aolongo da reta x = y = 0,5

●

●

●

0.0 0.1 0.2 0.3 0.4 0.5

0.00

0.05

0.10

0.15

z

|T−

Tex

ata|

●

Malhas

L/2L/4L/8L/16


Figura 24 – Solução aproximada para o quinto caso em alguns instantes de tempo

(a) t = 0,001 s

●

●

●●

●

●

●

0.0 0.2 0.4 0.6 0.8 1.0

0.00

0.02

0.04

0.06

0.08

x

|T−

Tex

ata|

●

Malhas

L/6L/12L/24L/48

(b) t = 0,01 s

●

● ●

●

● ●

●

0.0 0.2 0.4 0.6 0.8 1.0

0.00

0.02

0.04

0.06

0.08

x

|T−

Tex

ata|

●

Malhas

L/6L/12L/24L/48


espaçamento entre os nós, a solução tende a convergir para o resultado da soluçãoexata ao longo do tempo, mesmo não tendo alcançado o estado de equilíbrio que oproblema propicia pouco após o instante de tempo de 4,6 segundos (46 passos detempo com ∆t = 0,1 s, conforme visualizado na Figura 27).

Novamente, devido a solução exata ser proveniente de um somatório (Equação3.36), o problema está mais sujeito a problemas de arredondamento durante o cálculoda temperatura exata dos nós devido aos milhões de termos utilizados para alcançara solução desejada para posterior comparação. Assim, a comparação tende a ficarcomprometida devido a essa limitação da máquina.

Na Figura 26a pode-se verificar a dificuldade de convergir quando se trata das

60

Figura 25 – Solução aproximada para o quinto caso em alguns instantes de tempo (continuação)

(a) t = 0,05 s

●

●

●

●

●

●

●

0.0 0.2 0.4 0.6 0.8 1.0

0.00

0.02

0.04

0.06

0.08

x

|T−

Tex

ata|

●

Malhas

L/6L/12L/24L/48

(b) t = 0,1 s

● ● ● ● ● ● ●

0.0 0.2 0.4 0.6 0.8 1.0

0.00

0.02

0.04

0.06

0.08

x|T

−T

exat

a|

●

Malhas

L/6L/12L/24L/48


Figura 26 – Solução aproximada para o quinto caso em alguns instantes de tempo

(a) t = 0,001 s (b) t = 0,01 s

(c) t = 0,05 s (d) t = 0,1 s

0.3 0.60.2 0.70.1 0.80.50.4 10.90

Temperatura (K)


61

Figura 27 – Comportamento da temperatura máxima no dominio ao longo dos passos de tempo, com ∆t= 0,1 s

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●●

●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

0 10 20 30 40 50

0.0

0.2

0.4

0.6

0.8

1.0

Passo de tempo

Tm

ax (

K)

●

Malhas

L/6L/12L/24L/48


fronteiras do domínio, isso se dá devido a perda de precisão nas fronteiras devido autilização da equação do balanço de energia, em que os derivativos se utilizam demenos nós para fazer as interpolações.

63

5 RESULTADOS

5.1 Aplicação 1

A geometria do primeiro caso de aplicação é simulada pelo código, com malha(2397 nós) e distribuição de temperatura calculada representados nas Figuras 28a e28b. Também é mostrado o gráfico com o valor do η (critério de convergência) e dosvalores máximos e mínimos plotados ao longo dos passos de tempo até ser alcançadaa convergência (� = 10−6) nas Figuras 29a e 29b.

Figura 28 – Distribuição de temperatura calculada para a primeira geometria de aplicação

(a) Vista isométrica (b) Vista da base


Figura 29 – Erro entre os passos de tempo subsequente (a) e e evolução das temperaturas limites nageometria ao longo dos passos de tempo (b) para a primeira aplicação

(a)

0 50 100 150

1e−

065e

−06

5e−

055e

−04

Passo de tempo

log(

η)

(b)

0 50 100 150

300

305

310

315

320

Passo de tempo

T (

K)

Temperaturas

TminTmax


64

Figura 30 – Distribuição de temperatura calculada para a segunda geometria de aplicação

(a) Vista isométrica (b) Vista da base


Ainda na Figura 29b é possível verificar que a solução converge, com o critério deconvergência especificado, após 160 passos de tempo, logo, 26,7 minutos (∆t = 10 s).

5.2 Aplicação 2

A geometria do segundo caso de aplicação é simulada pelo código, com malha(5345 nós) e distribuição de temperatura calculada representados nas Figuras 30a e30b. Também é mostrado o gráfico com o val

Documents

Análise Numérica de Problemas Térmicos, Tridimensionais ...repositorio.roca.utfpr.edu.br/jspui/bitstream/1/... · de Conclusão de Curso 2, do Curso de En-genharia Mecânica do