Embed Size (px)
DESCRIPTION
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 naresoluçã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.
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.
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.
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
A seguir a malha para 16 elementos quadriláteros :
Malha de 16 elementos Quadriláteros e 16 elementos 1-D de contorno
A numeração dos nós do software GID é diferente da maioria das literaturas, no
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
A seguir a malha para 100 elementos quadriláteros :
Malha de 100 elementos Quadriláteros e 40 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, o resultado para t = 1 segundo foi de
0.291255, um pouco mais próximo do resultado analítico de 0.2947.
27
A seguir a malha para 225 elementos quadriláteros :
Malha de 458 elementos Quadriláteros e 84 elementos 1-D de contorno
A numeração dos nós do software GID é diferente da maioria das literaturas, no
caso, U8 é o resultado para o nó em x=0 e y=0, o resultado para t = 1 segundo foi de
0.291959, um pouco mais próximo do resultado analítico de 0.2947.
28
A seguir a malha para 458 elementos quadriláteros
Malha de 458 elementos Quadriláteros e 84 elementos 1-D de contorno
A numeração dos nós do software GID é diferente da maioria das literaturas, no
caso, U7 é o resultado para o nó em x=0 e y=0, o resultado para t = 1 segundo foi de
0.292689, um pouco mais próximo do resultado analítico de 0.2947.
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
A numeração dos nós do software GID é diferente da maioria das literaturas, no
caso, U3 é o resultado para o nó em x=0 e y=0, o resultado para t = 1 segundo foi de
0.333289, um pouco mais próximo do resultado analítico de 0.2947.
30
Abaixo a malha para 8 elementos triangulares:
Malha de 8 elementos Triangulares 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, o resultado para t = 1 segundo foi de
0.269657, um pouco mais próximo do resultado analítico de 0.2947.
31
Abaixo a malha para 106 elementos triangulares:
Malha de 106 elementos Triangulares e 28 elementos 1-D de contorno
A numeração dos nós do software GID é diferente da maioria das literaturas, no
caso, U1 é o resultado para o nó em x=0 e y=0, o resultado para t = 1 segundo foi de
0.289970, um pouco mais próximo do resultado analítico de 0.2947.
32
Abaixo a malha para 208 elementos triangulares:
Malha de 208 elementos Triangulares e 40 elementos 1-D de contorno
A numeração dos nós do software GID é diferente da maioria das literaturas, no
caso, U1 é o resultado para o nó em x=0 e y=0, o resultado para t = 1 segundo foi de
0.290630, um pouco mais próximo do resultado analítico de 0.2947.
33
Abaixo a malha para 444 elementos triangulares:
Malha de 444 elementos Triangulares e 56 elementos 1-D de contorno
A numeração dos nós do software GID é diferente da maioria das literaturas, no
caso, U6 é o resultado para o nó em x=0 e y=0, o resultado para t = 1 segundo foi de
0.292522, um pouco mais próximo do resultado analítico de 0.2947.
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
Número de Elementos Quadriláteros de 4 nós
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
Número de Elementos Quadriláteros de 4 nós
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
Número de Elementos Triangulares de 3 nós
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
Número de Elementos Triangulares de 3 nós
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
