Elementos Finitos - Transferência de calor Transiente

UNIVERSIDADE ESTADUAL DE CAMPINAS

FACULDADE DE ENGENHARIA CIVIL,

ARQUITETURA E URBANISMO

ROMILDO APARECIDO SOARES JUNIOR

RELATÓRIO DO CURSO DE ELEMENTOS FINITOS :

TRANSFERÊNCIA DE CALOR TRANSIENTE PELO

MÉTODO DOS ELEMENTOS FINITOS

CAMPINAS – SP

DEZEMBRO 2012

ROMILDO APARECIDO SOARES JUNIOR

RELATÓRIO DO CURSO DE ELEMENTOS FINITOS :

TRANSFERÊNCIA DE CALOR TRANSIENTE PELO

MÉTODO DOS ELEMENTOS FINITOS

Relatório apresentado como parte dos

requisitos para aprovação na disciplina

Introdução ao Método dos Elementos Finitos do

Curso de Pós Graduação em Engenharia Civil da

Universidade Estadual de Campinas, elaborado

sob orientação do Prof. Dr. Philippe Remy

Bernard Devloo.

CAMPINAS – SP

DEZEMBRO 2012

DEDICATÓRIA

O seguinte relatório é dedicado as pessoas que tornaram possível a sua realização.

Sendo estas a minha família, que sempre estiveram ao meu lado e tornaram mais amenos

todos momentos difíceis passados em minha vida.

AGRADECIMENTOS

Agradeço à minha família por seu apoio, aos amigos e também ao Prof. Dr. Philippe

Remy Bernard Devloo que com seu empenho em ensinar o método dos elementos finitos,

possibilitou a conclusão deste trabalho, por sua transparência ao repassar todo o

conhecimento adquirido ao longo dos anos de sua carreira.

RESUMO

O relatório detalhará os métodos utilizados para resolução pelo método dos

elementos finitos da equação diferencial parcial parabólica de duas dimensões, quando

utilizada em problemas de calor transiente. Serão utilizados métodos computacionais na

resolução do problema. Será feita também a comparação dos resultados obtidos pelo método

dos elementos finitos com a solução exata do problema.

Palavras-chave: Elementos Finitos, Transferência de Calor Transiente

ABSTRACT

This work will detail the methods used to solve the parabolic partial differential

equation in two dimensions by the finite element method, when used in transient heat transfer

problems. Computational methods will be used in solving the problem. It will also compare

the results obtained by the finite element method with the exact solution of the problem.

Keywords: Finite Elements Method, Transient Heat Transfer

SUMÁRIO

REUMO 05

1. INTRODUÇÃO 07

2. FORMULAÇÃO FORTE DO PROBLEMA ANALISADO 08

3. FORMULAÇÃO PARA ELEMENTOS FINITOS 10

4. INTEGRAÇÃO NO TEMPO 12

5. O PROBLEMA ANALISADO 13

6. RESOLUÇÃO DO PROBLEMA 14

7. RESOLUÇÃO DO PROBLEMA PARA DIVERSAS MALHAS

23

8. ESTIMATIVA DE ERRO E CONVERGÊNCIA

34

9. VISUALIZAÇÃO DOS RESULTADOS

41

10. CONCLUSÕES

43

11. REFERÊNCIAS BIBLIOGRÁFICAS

44

7

1. INTRODUÇÃO

O objetivo deste trabalho é verificar a convergência dos resultados de problemas

de transferência de calor transiente em duas dimensões ( sendo o domínio em 2d mais a

variação do tempo ) utilizando o método dos elementos finitos. O problema analisado possui

apenas um grau de liberdade, mas possui a variação do tempo, portanto o método dos

elementos finitos conduzirá a um sistema final de equações diferenciais parabólicas para o

problema analisado, sendo então necessária a integração deste sistema no tempo. O método

utilizado para a integração no tempo será o Runge Kutta de quarta ordem. Para resolução do

problema, será utilizado um software desenvolvido durante as aulas do Prof. Dr. Philippe

Remy Bernard Devloo, utilizando a biblioteca para elementos finitos chamada NeoPZ,

podendo ser encontrada em http://code.google.com/p/neopz/. Esta biblioteca permite

utilização de várias ferramentas úteis na programação de software para elementos finitos,

como a resolução de equações diferenciais, manipulação de matrizes, verificação de erro,

visualização da malha entre outros. A linguagem de programação utilizada foi o C++ e o

ambiente de desenvolvimento foi o Visual Studio 2010. Neste trabalho serão detalhados os

métodos utilizados para resolução do problema, apresentando a solução forte do problema, a

formulação para elementos finitos ( utilizando o método de Galerkin ) e também os métodos

computacionais utillizados. Por fim, será feita a comparação dos resultados vindos da solução

analítica com os resultados obtidos por diversas malhas de elementos finitos, a medida em

que se aumentava o número de elementos na malha. Os elementos utilizados para resolução

do problema foram os triângulos e quadriláteros. Para gerar-se as malhas dos elementos, foi

utilizado o software chamado GiD, podendo ser encontrado em http://gid.cimne.upc.es/ . Este

programa é um ótimo gerador de malhas, mas também possui a capacidade de se aplicar

materiais nos elementos e também gerar malhas de contorno.

http://code.google.com/p/neopz/

http://gid.cimne.upc.es/

8

2. FORMULAÇÃO FORTE DO PROBLEMA ANALISADO

A relação constitutiva para este problema de Transferência de Calor é que o fluxo

de energia é proporcional ao gradiente da função U(x,y,t).

QtyxUGradk )],,([

A densidade do fluxo por unidade de área em um determinado ponto é dada por :

][

Qdiv

De acordo com a lei de conservação de energia, o fluxo produzido pelas forças

internas é o mesmo que o fluxo que saem pelo contorno do domínio. Então :

nQdQdiv ][

Tendo n como o vetor normal do fluxo que passa no domínio. Utilizando o

teorema da divergência é possível transformar esta integral de contorno em uma integral de

domínio:

dfdnQ

Tendo, portanto :

dfdQdiv ][

9

A taxa de variação da energia térmica no corpo é igual à diferença gerada pelas

forças internas e a energia que sai pelo contorno do corpo. A temperatura, em um dado

instante de tempo, é dada por :

dfdQdivdt

tyxUCv ][

),,(

Sendo que Cv é a Capacidade Térmica (quantidade de energia necessária para

aumentar a temperatura de uma unidade de massa do corpo, Joule/g), p é a densidade de

massa do corpo analisado. Substituindo a primeira relação constitutiva, onde o fluxo de

energia é proporcional ao gradiente da função, temos a seguinte equação :

dfdtyxUGradkdivdt

tyxUCv )]],,([[

),,(

Tendo em vista que o divergente do operador gradiente é o laplaciano da função,

temos:

dfdtyxUkdt

tyxUCv ),,(

),,(

Como o domínio é arbitrário, e se as funções acima forem suaves o suficiente, a

equação acima é dada por :

ftyxUkt

tyxUCv

),,(

),,(

Sendo esta a formulação forte do problema. Com as seguintes condições de

contorno :

10

0)0,,(

)(

)(),,(

yxU

NeumannqemqnQ

DirichlettemUtyxU

3. FORMULAÇÃO PARA ELEMENTOS FINITOS

Faremos agora a formulação para elementos finitos, utilizaremos o método de

Galerkin. Multiplicaremos cada parte de nossa formulação forte por uma função teste v, onde

H0 é o espaço que contém todas as funções para que as integrais a seguir possam ser

consideradas verdadeiras.

dvfdvtyxUkdvt

tyxUCv ),,(

),,(

Agora realizamos a integração por partes do operador Laplaciano, utilizando o

teorema de Green-Gauss, obtendo as seguintes equações :

n

uvvGradtyxUGradkdvtyxuk ][)],,([),,(

Onde a ultima parcela da integral, determina a contribuição do elemento de

contorno, no caso o valor do fluxo em um dado contorno, caso a condição de contorno seja

neumann, sendo então somado ao vetor de carga. A formulação então fica :

dvfn

uvvGradtyxUGradkdv

t

tyxUCv ][)],,([

),,(

11

Rearranjando os termos temos :

n

uvdvfvGradtyxUGradkdv

t

tyxUCv ][)],,([

),,(

Se dividirmos esse espaço de funções em subespaços de funções :

n

uvhdvhfvhGradtyxUhGradkdvh

t

tyxUhCv ][)],,([

),,(

A integral no contorno será chamada de q

n

u

pois delimita a contribuição das

condições de contorno. Apartir do método de Galerkin temos que :

n

j

bjNjvh1

n

i

NitUiUh1

)(

qbjNjdbjNjf

bjNjGradNitUiGradkdbjNjt

NitUiCv ][])([

)(

Rearranjando :

qNjdNjf

NjGradNiGradkUidNiNjCvt

Ui n

i

n

i 11

])[][()(

12

Temos então o seguinte sistema de Equações Diferenciais Parciais de Ordem 1:

fiKijUiMijt

Ui n

i

n

i

11

Onde :

dNiNjCvMij

dNjGradNiGradkKij ][][

qNjdNjffi

4. INTEGRAÇÃO NO TEMPO

É possível utilizar vários métodos de Integração no tempo para este sistema linear,

mas neste trabalho será utilizado o Runge-Kutta de Quarta Ordem. Abaixo, a formulação

utilizada :

0)0(

),('

U

UtfU

Primeiramente isola-se a derivada, tendo-se uma equação em relação a t e U. Seja

h o tamanho do espaço de tempo escolhido, temos que :

)432221(6

11

)3,(4

)2

2,

2(3

)2

1,

2(2

),(1

00

kkkkwiwi

kwihtifhk

kwi

htifhk

kwi

htifhk

witifhk

w

13

Onde wi é a resposta no tempo analisado. Deve-se escolher um espaço de tempo

que possa atender determinada equação, espaços de tempo grandes podem gerar instabilidade

das respostas. Para o problema analisado, deveremos integrar o seguinte sistema :

)(' 1 KijUifiMijU

5. O PROBLEMA ANALISADO

O problema analisado consiste em um problema 2-D de transferência de calor

transiente, pois este se varia com o tempo, conforme é transferido o calor gerado pelas forças

internas. A equação diferencial do problema é a seguinte :

1),,(),,(

tyxU

t

tyxU

As condições de contorno para o problema analisado são as seguintes :

0),,0(

x

tyU

0

),0,(

y

txU

0),,1( tyU

0),1,( txU

A condição inicial para o tempo é :

0)0,,( yxU para todo (x,y) em

14

Abaixo, um esquema do problema analisado :

6. RESOLUÇÃO DO PROBLEMA

Para resolver o problema, será utilizado o método dos elementos finitos explicado

anteriormente. Será utilizado um software baseado na biblioteca para elementos finitos

chamada “NeoPZ” oferecida pelo professor Philippe Remy Bernard Devloo, podendo ser

encontrada no link : http://code.google.com/p/neopz/, o qual possui ferramentas poderosas

para análise de matrizes, resolução de equações diferenciais e sistemas lineares, entre outros.

O programa foi inteiramente desenvolvido na linguagem de programação C++, o ambiente de

desenvolvimento escolhido para este trabalho foi o Visual Studio 2010. O computador

utilizado para realizar este trabalho foi um computador com processador Amd Phenom X6

com 8 Gb de memória ram, utilizando o Windows 7 Ultimate. O programa desenvolvido

segue o seguinte fluxograma simplificado:

15

Fluxograma do programa desenvolvido

A exemplo de cálculo, será calculado a matriz de rigidez global deste problema

para apenas 1 elemento quadrilátero. Primeiramente, é feita uma mudança de variáveis entre

o elemento deformado (com coordenadas x e y) para o elemento mestre. O elemento mestre é

um elemento com as coordenadas e . Para o elemento quadrilátero, o elemento mestre é

isoparamétrico, tendo as coordenadas de seu domínio (-1,-1), (1,-1), (1,1), (-1,1). É feita essa

mudança de variáveis a fim de simplificar o cálculo das funções de forma e das integrais

subsequentes. Por esse método, será utilizada a integração de Gauss Legendre, sendo obtidos

pontos de integração em função do número de nós do elemento em questão. A mudança de

variáveis é feita como o esquema a seguir :

16

Para a figura acima, é apresentado o elemento quadrilátero de quatro nós, as

funções de forma que denotam essa mudança de variáveis são as seguintes:

)1)(1(4

1),(1 N

)1)(1(4

1),(2 N

)1)(1(4

1),(3 N

)1)(1(4

1),(4 N

Onde:

N

i

xiNix1

),(

N

i

yiNiy1

),(

17

Para o cálculo das integrais das matrizes de rigidez e capacidade térmica é

necessário o cálculo do jacobiano da mudança de variáveis.

As derivadas do jacobiano são computadas da seguinte maneira :

xiNix N

i

1

),(

yiNiy N

i

1

),(

xiNix N

i

1

),(

yiNiy N

i

1

),(

As derivadas do jacobiano são computadas da seguinte maneira :

)(]0,0[1

JDet

y

J

)(

]1,0[1

JDet

y

J

18

)(]0,1[1

JDet

x

J

)(

]1,1[1

JDet

x

J

Com o determinante do jacobiano, e a Inversa do jacobiano é possível calcular as

derivadas reais para calculo das matrizes de rigidez e capacidade térmica.

)),(),(

(||

1

11

yiNiNi

yiNiNi

Jx

Ni N

i

N

i

)),(),(

(||

1

11

xiNiNi

xiNiNi

Jy

Ni N

i

N

i

Como visto anteriormente, a matriz de rigidez é dada pela equação:

dNjGradNiGradkKij ][][

Integrando a matriz de rigidez no elemento mestre :

dAJy

Nj

y

Ni

x

Nj

x

NiKij ||)(

1

1

1

1

Como todos os termos foram calculados para cada ponto de Gauss, Kij fica sendo

o somatório das contribuições de cada ponto de Gauss.

||)()()( JjWiWy

Nj

y

Ni

x

Nj

x

NiKij

19

A matriz de Capacidade Térmica tem a seguinte equação :

dNiNjCvMij

No problema analisado, Cv =1, p = 1. Quando integrada no elemento mestre :

dAJNjNiMij ||)(1

1

1

1

Como todos os termos foram calculados para cada ponto de Gauss, Mij fica sendo

o somatório das contribuições de cada ponto de Gauss.

||)()()( JjWiWNjNiMij

A contribuição do vetor de carga, neste material do elemento é de 1.

dNiffi

1f

dNifi 1

Quando integrada no elemento mestre :

dAJNifi ||)(1

1

1

1

20

Como todos os termos foram calculados para cada ponto de Gauss, Fi fica sendo o

somatório das contribuições de cada ponto de Gauss.

||)()()( JjWiWNifi

Com as três matrizes Kij, Mij e Fi para cada elemento quadrilátero (ou triangular)

é necessário montar as matrizes globais, Kij global, Mij global e Fi global. Para cada matriz

global, são somadas as contribuições nos nós de cada elemento da malha, no caso, esta matriz

possui apenas um grau de liberdade. Preparada a matriz global, é necessário aplicar as

condições de contorno do domínio. No programa, as condições de contorno Dirichlet são

aplicadas utilizando a utilizando a técnica do “número grande”, onde os elementos 1-D de

contorno gerados pela malha, determinam os nós onde serão somados com um número grande

escolhido. A técnica do numero grande permite aplicar as condições de contorno sem precisar

alterar o tamanho da matriz global, esta técnica pode gerar erros, mas para o caso analisado, é

aceitável. O elemento 1-D aqui analisado, é calculado de maneira similar aos elemento 2-D,

primeiro é feita a mudança de variáveis, utilizando – se as funções de forma do elemento

mestre. As funções de forma lineares para o elemento 1-D são :

2

1)(1

N

2

1)(2

N

Onde :

N

i

xiNix1

)(

O jacobiano desta mudança de variáveis vai ser :

))2()1(())2()1(())2()1(())2()1(( nóynóynóynóynóxnóxnóxnóxdelx

21

2][

delxJ

delxJ

21

Se a condição de contorno for Dirichlet, a contribuição na matriz global será :

||)()()( JBigjWiWNjNiKij

||)()()( JBigjWiWNjNiMij

||)()( JBigjWiWNifi

Se a condição de contorno for Neumann, a contribuição na matriz global será

apenas na matriz do vetor de carga :

||)()( JjWiWNifi

A matriz global já modificada com as condições de contorno, está pronta para ser

integrada utilizando qualquer método de integração no tempo. No caso será utilizado o Runge

Kutta de quarta ordem, vindo da Biblioteca URKF45, uma biblioteca especializada em

resolução de sistemas de equações diferenciais. Esta biblioteca pode ser encontrada em

http://jean-pierre.moreau.pagesperso-orange.fr/. A exemplo de cálculo, a matriz global de

apenas 1 elemento quadrilátero, para este problema é :

22

Aplicadas as condições de contorno, a disposição dos big numbers fica da

seguinte maneira no sistema de equações, a fim de “anular” as linhas e colunas onde a

temperatura Ui = 0. Como é possível ver, a disposição dos big numbers é a mesma tanto para

a matriz Kij quando para a matriz Mij. A matriz global final para apenas 1 elemento

quadrilátero fica então da seguinte maneira :

Como temos quatro nós apenas, devido as condições de contorno somente o nó 1

(x = 0, y = 0) terá variação de temperatura ( todos os outros tem U = 0 ). A solução deste

sistema de equações diferenciais de primeira ordem é feito pelo Runge Kutta de quarta ordem,

tendo as seguintes temperaturas obtidas pelo programa :

Resultados Obtidos para 1 elemento quadrilátero

23

O comando NDSolve do programa Wolfram Mathematica obteve o seguinte

resultado, resolvendo-se o sistema de equações diferenciais :

Resultados Obtidos pelo Mathematica para 1 elemento quadrilátero

A temperatura encontrada para t = 1 segundo, foi de 0.374459 pelo Runge Kutta

de quarta ordem e de 0.37407 pelo software Wolfram Mathematica. A solução Analítica do

problema principal ( resolvendo a equação parcial parabólica do problema ) tem o resultado

para esse nó de 0.2947 em t = 1 segundo. A diferença entre o resultado numérico e o analítico

ainda é muito grande. Isto é devido ao cálculo de apenas um elemento quadrilátero.

7. RESOLUÇÃO DO PROBLEMA PARA DIVERSAS MALHAS

Será feito agora a comparação dos resultados para 4, 16, 100, 225 e 458 elementos

quadriláteros. Agora, o input de todos os pontos X e Y de cada nó terá dificuldades se for

feito manualmente. Para isso será utilizado o software GID, um ótimo gerador de malhas para

elementos finitos. Além de gerar malhas, o GID também pode colocar os materiais em cada

elemento de maneira dinâmica por meio de sua própria programação. O programa GID pode

ser encontrado em http://gid.cimne.upc.es/. O programa utilizado para cálculo dos elementos

finitos teve de ser adaptado para ler os arquivos exportados pelo programa GID, além disso

foi necessário observar a numeração dos nós exportados pelo GID, pois alguns nós podem ser

numerados de maneira diferente.

http://gid.cimne.upc.es/

24

7.1 Malhas de Elementos Quadriláteros

A seguir a malha para 4 elementos quadriláteros :

Malha de Quatro elementos Quadriláteros e 8 elementos 1-D de contorno

A numeração dos nós do software GID é diferente da maioria das literaturas, no

caso, U5 é o resultado para o nó em x=0 e y=0 (anteriormente U1), o resultado para t = 1

segundo foi de 0.308421, um pouco mais próximo do resultado analítico de 0.2947.

25


Malha de 16 elementos Quadriláteros e 16 elementos 1-D de contorno


caso, U1 é o resultado para o nó em x=0 e y=0, o resultado para t = 1 segundo foi de

0.298434, um pouco mais próximo do resultado analítico de 0.2947. Como pode ser

observado o refinamento da malha provoca uma aproximação cada vez melhor da resposta

analítica. A seguir, será analisado a malha para uma quantidade de elementos muito maior,

100, 225 e 458 elementos.

26





0.291255, um pouco mais próximo do resultado analítico de 0.2947.

27






28

A seguir a malha para 458 elementos quadriláteros





29

7.2 Malhas de Elementos Triangulares

Abaixo a malha para 2 elementos triangulares:

Malha de 2 elementos Triangulares e 4 elementos 1-D de contorno




30






31






32






33






34

8. ESTIMATIVA DE ERRO E CONVERGÊNCIA

8. 1 Solução Analítica

A solução analítica do problema é dada pela seguinte equação, em t = 1 segundo :

Onde n pode ser somado até o infinito, conforme a necessidade de precisão. Para

se obter a precisão básica é necessário o cálculo de pelo menos os 20 termos iniciais desta

série. Houve problemas quando se tentou calcular muitos termos desta série na linguagem

C++, pois os cossenos hiperbólicos voltavam valores muito altos, sendo necessário encontrar

um valor médio entre a precisão e a demora de cálculo.

8. 2 Estimativa de Erro para Elementos Quadriláteros

O erro básico é dado pela fórmula :

),,(),,(),,( tyxuhtyxutyxe

Para este tipo de erro, será analisado o nó x=0 e y=0 e para o tempo de 1 segundo.

Abaixo a estimativa de erro básico :

Resposta Analítica Resposta pelo MEF Número de Elementos Erro Básico

0,2947 0,374459 1 0,079759

0,2947 0,308421 4 0,013721

0,2947 0,298434 16 0,003734

0,2947 0,291255 100 0,003445

0,2947 0,291959 225 0,002741

0,2947 0,292689 458 0,002011

Tabela do erro básico

35

Abaixo, um gráfico representando o erro básico no nó (0,0) para t = 1 segundo.

Gráfico Erro Básico

Foi observado que no ponto (0,0) em malhas refinadas acima de 450 elementos, o

número de equações na matriz global influenciava na integração por Runge-Kutta, o sistema

necessitava de um step-size de tempo menor que o adotado. O que, devido aos métodos de

resolução utilizados, torna inviável a obtenção da resposta, pois com um step size de 0,001 o

sistema demorou aproximadamente 3 horas para resolver a malha de 458 elementos

quadriláteros ( Amd Phenom X6 2.8 Ghz, 8gb ram). Se for diminuído ainda mais o step size

será necessário reformulação do método de integração no tempo, ou processamento paralelo

para se obter respostas em tempo viável.

Abaixo, será analisado o erro médio das malhas de elementos quadriláteros. O

erro médio é dado pela fórmula :

dAtyxuhtyxue 2

0 )),,(),,((||||

O erro médio avalia o erro de todos os pontos no domínio. No caso, será avaliado

o erro no tempo de 1 segundo. Serão mostrados apenas os nove pontos nós deste problema

0

0,01

0,02

0,03

0,04

0,05

0,06

0,07

0,08

0,09

1 4 16 100 225 458

Erro

Bás

ico

Número de Elementos Quadriláteros de 4 nós

Erro Básico

Erro Básico

36

nesta tabela, mas o cálculo considera todos os nós calculados para cada malha gerada

inicialmente.

N Solução Analítica

MEF 4 Elementos

MEF 16 Elementos

MEF 100 Elementos

MEF 225 Elementos

MEF 458 Elementos

1 0,2947 0,308421 0,298434 0,291255 0,291959 0,292689

2 0,2293 0,239277 0,228912 0,229084 0,22917 0,229184

3 0 0 0 0 0 0

4 0,2293 0,239529 0,23293 0,226958 0,227846 0,228834

5 0,1811 0,191687 0,18432 0,183283 0,182543 0,181864

6 0 0 0 0 0 0

7 0 0 0 0 0 0

8 0 0 0 0 0 0

9 0 0 0 0 0 0

Erro médio Erro médio Erro médio Erro médio Erro médio

0,022414 0,009576 0,00470798 0,003424375 0,002204184

Tabela do erro médio para cada malha analisada

Abaixo, os gráficos para os erros médios de cada malha analisada.

Erro Médio :

Gráfico do erro médio

0

0,005

0,01

0,015

0,02

0,025

4 16 100 225 458

Erro

Mé

dio


Erro Médio

Erro Médio

37

Erro médio em escala logarítmica :

Gráfico do erro médio em escala logaritmica

8. 3 Estimativa de Erro para Elementos Triangulares

O erro básico é dado pela fórmula :

),,(),,(),,( tyxuhtyxutyxe

Para este tipo de erro, será analisado o nó x=0 e y=0 e para o tempo de 1 segundo.

Abaixo a estimativa de erro básico :

Resposta Analítica Resposta pelo MEF Número de Elementos Erro Básico

0,2947 0,333289 2 0,038589

0,2947 0,269657 8 0,025043

0,2947 0,289970 106 0,00473

0,2947 0,291023 208 0,003677

0,2947 0,291122 444 0,003578

Tabela do erro básico – Elementos Triangulares

0,001

0,01

0,1

1

4 16 100 225 458

Erro

Mé

dio


Erro Médio (escala logarítmica)

Erro Médio

38

Abaixo, um gráfico representando o erro básico no nó (0,0) para t = 1 segundo.

Gráfico Erro Básico – Elementos Triangulares

Foi observado que no ponto (0,0) em malhas refinadas acima de 450 elementos, o

número de equações na matriz global influenciava na integração por Runge-Kutta, o sistema

necessitava de um step-size de tempo menor que o adotado. O que, devido aos métodos de

resolução utilizados, torna inviável a obtenção da resposta, pois com um step size de 0,001 o

sistema demorou aproximadamente 2.8 horas para resolver a malha de 444 elementos

quadriláteros ( Amd Phenom X6 2.8 Ghz, 8gb ram). Se for diminuído ainda mais o step size

será necessário reformulação do método de integração no tempo, ou processamento paralelo

para se obter respostas em tempo viável.

Abaixo, será analisado o erro médio das malhas de elementos triangulares. O erro

médio é dado pela fórmula :

dAtyxuhtyxue 2

0 )),,(),,((||||

O erro médio avalia o erro de todos os pontos no domínio. No caso, será avaliado

o erro no tempo de 1 segundo. Serão mostrados apenas os nove pontos nós deste problema

nesta tabela, mas o cálculo considera todos os nós calculados para cada malha gerada

inicialmente.

0

0,005

0,01

0,015

0,02

0,025

0,03

0,035

0,04

0,045

2 8 106 208 444

Erro

Bás

ico

Número de Elementos Triangulares de 3 nós

Erro Básico

Erro

39

N Solução Analítica

MEF 2 Elementos

MEF 8 Elementos

MEF 106 Elementos

MEF 208 Elementos

MEF 444 Elementos

1 0,2947 0,333289 0,269657 0,289970 0,290062 0,292522

2 0,2293 - 0,228143 0,229084 0,229201 0,229223

3 0 0 0 0 0 0

4 0,2293 - 0,228143 0,228762 0,228895 0,228896

5 0,1811 - 0,155232 0,165429 0,169944 0,175314

6 0 - 0 0 0 0

7 0 0 0 0 0 0

8 0 - 0 0 0 0

9 0 0 0 0 0 0

Erro médio Erro médio Erro médio Erro médio Erro médio

0,040751 0,033217 0,014771 0,0098951 0,0085282

Tabela do erro médio para cada malha analisada – Elementos Triangulares

Abaixo, os gráficos para os erros médios de cada malha analisada.

Erro Médio :

Gráfico do erro médio – Elementos Triangulares

0

0,005

0,01

0,015

0,02

0,025

0,03

0,035

0,04

0,045

2 8 106 208 444

Erro

Mé

dio


Erro Médio

Erro Médio

40

Erro médio em escala logarítmica :

Gráfico do erro médio em escala logarítmica – Elementos Triangulares

8. 4 Testes de verificação dos resultados

Podem ser feitos alguns testes para verificação dos resultados além dos testes de

erro com a solução analítica. Pode ser feito a verificação da matriz rigidez para o domínio não

distorcido, sendo que esta deve ser simétrica e ter suas diagonais subsequentes iguais.

Também pode-se realizar o teste onde se gira o domínio, tendo em vista que o resultado da

matriz de rigidez deve ser o mesmo, onde este foi verificado utilizando – se as funções seno e

cosseno para girar as coordenadas de cada nó. Como foi verificado a seguir ( sendo esta a

matriz para 1 elemento quadrilátero não distorcido ) :

0,001

0,01

0,1

1

2 8 106 208 444

Erro

Mé

dio


Erro Médio

Erro Médio

41

9. VISUALIZAÇÃO DOS RESULTADOS

9.1 Gráficos

Abaixo a visualização dos resultados obtidos, onde é possível verificar a variação

da temperatura conforme o tempo, na placa retangular:

t = 0.1 s t = 0.2 s

t = 0.4 s

42

t = 0.6 s t = 0.8 s

t = 1 s

É possível verificar que gradativamente, a temperatura se concentra no canto

inferior esquerdo, nas coordenadas x = 0 e y = 0, ponto onde foram obtidas as maiores

temperaturas.

43

10. CONCLUSÕES

O método dos elementos finitos se adequou na resolução das equações parabólicas

de uma ótima maneira, pois seus resultados convergiram com o refinamento da malha. No

entanto, a margem de erro para os métodos de cálculo utilizados pode ser melhorado. Por

exemplo, na remoção da utilização do big number na aplicação das condições de contorno,

utilizando em seu lugar métodos de manipulação de matrizes. Outro exemplo onde é possível

melhorar a margem de erro é na utilização de um step de tempo menor que o utilizado, mas

será necessário máquinas mais potentes ou outros métodos de integração no tempo.

44

11. REFERÊNCIAS BIBLIOGRÁFICAS

REDDY, J. N., AN INTRODUCTION TO THE FINITE ELEMENT METHOD, SECOND

EDITION, 1993

BECKER, E. B., THE FINITE ELEMENTS - AN INTRODUCTION, VOLUME I, 1981

AZEVEDO, F. M., MÉTODO DOS ELEMENTOS FINITOS, PRIMEIRA EDIÇÃO, 2003

Documents

Elementos Finitos - Transferência de calor Transiente