View
66
Download
0
Category
Preview:
Citation preview
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 1/134
ELISEU LEANDRO MAGALHÃES MONTEIRO
MODELAÇÃO NUMÉRICA DA SOLIDIFICAÇÃO EM FUNDIÇÃO
Tese submetida à Universidade de Trás-os-Montes e Alto Douro para a obtenção do grau de Mestre.
UNIVERSIDADE DE TRÁS-OS-MONTES E ALTO DOURO
2003
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 2/134
Agradecimentos
O autor deseja expressar os seus agradecimentos a todas as pessoas que directa ou
indirectamente possibilitaram a realização deste trabalho, em especial ao Professor Abel
Rouboa (UTAD) e Professor Caetano Monteiro (UM), pela disponibilidade que
manifestaram em esclarecer as minhas dúvidas e pelos seus oportunos comentários a
este trabalho.
Não posso também deixar de expressar os meus agradecimentos à UTAD, em particular
ao departamento de Engenharia Mecânica, pelo apoio logístico.
Aos meus familiares e amigos pelo incentivo por sempre acreditarem na conclusão bem
sucedida deste trabalho.
ii
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 3/134
Resumo
O objectivo central desta dissertação consiste na aplicação do método dos volumes
finitos à modelação da solidificação em fundição em moldação metálica, usando
coordenadas curvilíneas generalizadas.
As peças de geometria complexa constituem parte substancial do universo das peças
obtidas por fundição, cujo sucesso de produção é em grande parte determinado pela
mesma geometria. Esta afecta não só as condições de solidificação, mas também o
próprio arrefecimento. Por estes motivos a sua modelação assume primordial
importância com vista a eliminar, ou pelo menos reduzir, a necessidade de introduzir
correcções na moldação por defeito de projecto, o que pode por em causa a viabilidade
económica da produção.
A modelação da solidificação em fundição encerra dois níveis de complexidade:
primeiro, a complexidade inerente aos fenómenos físicos envolvidos, e, segundo, a
complexidade de geometrias que os produtos fundidos geralmente apresentam. Para
descrever a evolução dos parâmetros termodinâmicos envolvidos neste fenómeno, é
necessário recorrer à equação diferencial da difusão, definindo as adequadas condições
iniciais e de fronteira, e então desenvolver um código numérico capaz de reproduzir
fielmente a transferência de calor entre peça, moldação e meio ambiente.
A utilização de coordenadas curvilíneas, transformando um domínio real de geometria
arbitrária num domínio computacional rectangular e único, permite uma certacomodidade na programação, apesar de aumentar a complexidade matemática do
modelo.
Afim de confirmar a validade do modelo numérico desenvolvido, este foi comparado
com soluções analíticas conhecidas em problemas de geometria simples, e no caso de
geometrias complexas foi feita comparação com resultados numéricos utilizando o
método das diferenças finitas e dados experimentais.
iii
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 4/134
Dos resultados obtidos pode-se concluir que o método dos volumes finitos se adapta de
modo bastante satisfatório ao rigor que se impõe na descrição de fenómenos complexos
em domínios de geometria complexa.
iv
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 5/134
Abstract
The main objective of the present Thesis is to develop an appropriate numerical
simulation of the solidification in castings, using the finite volume method with
formulation in generalise curvilinear coordinates.
Generally, the shapes produced in foundry industries are characterized by their complex
geometries. The production success is linked with the degree of shape geometry
complexity. This not only affects the solidification conditions but also the cooling. For
these reasons, it is important to reduce the amount of additional work needed to correct
the defects of the shape after cooling.
The simulation of the solidification process has two complexity levels: first, the
complexity of the heat transfer phenomenon, and, second, the complexity of the cast part
geometry. In order to develop a numerical code that allows describing the evolution of
heat parameters during this phenomenon, it is necessary to solve a heat conservation
differential equation using adequate boundary and initial conditions.
The complex geometry induces difficulties to obtain the numerical solution of those
differential equation using orthogonal coordinates. To overcome these difficulties, the
curvilinear formulation and an appropriate finite volume method are used. The validity
of the developed numeric model is confirmed by its comparison with experimental data
and numerical results using other approaches.
The Finite Volumes Method revealed to be well adapted to describe such complex
phenomena in complex geometry, which characterize the generality of the casting
practice cases.
v
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 6/134
Résume
L'objectif central de cette dissertation consiste en l'application de la méthode des
volumes finis, au modelage de la solidification dans le moulage métallique, utilisant des
coordonnées curvilignes généralisées.
En général, les pièces de géométries complexes constituent une partie substantielle des
pièces obtenues par moulage, dont le succès de la production est en grande partie
déterminé par la complexité de leurs géométries. Elle affecte non seulement les
conditions de solidification, mais aussi le refroidissement lui-même. Pour ces raisons, il
est important de réduire le travail supplémentaire qui consiste à corriger les défauts du
moulage. En effet, cette réduction nous permet d’économiser dans le coût de la
production des pièces moulées.
La modélisation de la solidification comporte deux niveaux de complexité :
premièrement, la complexité des phénomènes physiques, et, en second, la complexité de
la géométrie de produit fondu. Afin de développer un code numérique qui permet de
décrire l’évolution des paramètres thermodynamiques de ce phénomène, il est
nécessaire de résoudre l’équation différentielle de la conservation de l’énergie, en
utilisant la définition de la condition initiale et des conditions aux limites.
La complexité de la géométrie de la pièce rend très difficile la résolution de l’équation
différentielle. L'utilisation des coordonnées curvilignes est nécessaire ainsi que leurstransformation d’un domaine réel (géométrie complexe) vers un domaine de calcul
(géométrie rectangulaire). Ceci permet un confort dans la programmation, malgré
l’augmentation de la complexité mathématique du modèle.
Afin de valider le code développé utilisant un nouveau solveur, la confrontation avec les
résultats expérimentaux et numériques (utilisant la méthode des différences finies) a été
mise en évidence.
vi
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 7/134
Les résultats obtenus permettent de déduire que la méthode des volumes finis s'adapte
de manière satisfaisante à la description de phénomènes complexes dans des domaines
de la géométrie complexe.
vii
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 8/134
Lista de símbolos
Romanos
A matriz de coeficientes
a difusividade (m2/s)
ar
vector genérico
B matriz de termos independentesC matriz de pré-condicionamento
C p capacidade térmica mássica (J/kg ºC)
ds superfície elementar (m2)
F fluxo líquido através da fronteira do volume de controlo
f componente vectorial da difusão de calor
f s fracção sólida
h coeficiente Newtoniano de transferência de calor (W/m ºC) J Jacobiano da transformação
J 0 função de Bessel do primeiro grau de primeira ordem
J 1 função de Bessel do primeiro grau de segunda ordem
k condutibilidade térmica (W/m ºC)
L matriz triangular inferior
M matriz produto ( LU )
n
r
vector normal a superfície do volume de controloO termos da série de Taylor de ordem superior a dois
Q matriz de termos independentes
q fonte de calor gerada por mudança de fase (W/m3)
R Resíduo
r raio (m)
S superfície (m2)
t tempo (s)
viii
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 9/134
U matriz triangular superior
ur
vector velocidade (m/s)
x, y eixos do sistema cartesiano de coordenadas (m)
Y matriz de coeficientes
Gregos
α factor peso
β factor de relaxação
∆h f calor latente de transformação de fase (J/kg)
∆t incremento de tempo (s)
δ actualização do erro de convergência
ε erro de convergência
φ temperatura (ºC)
λ factor interpolador
ν valor próprio
θ variável genérica
ρ massa volúmica (Kg/m³)
Γ fronteira genérica de um domínio
∂ designação de derivada parcial
µ zeros das funções de Bessel do primeiro grau e primeira ordem
ψ vector próprio
∆ designação para variação discreta
∇
operador Nablaξ,η eixos do sistema de coordenadas curvilíneas
Índices
a ambiente
E Este
f fusão
ix
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 10/134
i,j coordenadas computacionais
k número da iteração
l líquido ou liquidus
m moldação
N Norte
NE Nordeste
NN Norte-Norte
NNW Nor-noroeste
NW Noroeste
n iteração temporal
nb nós vizinhos
P nó central de um volume de controlo
S Sul
SE Sudeste
SS Sul-Sul
SSE Su-sueste
SW Sudoeste
s sólido, solidus ou metal solidificadoW Oeste
x
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 11/134
Lista de figuras
Figura1.1 – Curvas de arrefecimento de um corpo metálico.......................................................................................... 3
Figura 1.2 – Mecanismo de formação de rechupes....................................................................................................................... 4
Figura 1.3 – Zonas de formação de rechupes numa peça.................................................................................................... 5
Figura 2.1 – Discretização de uma geometria de forma arbitrária..................................................................... 20
Figura 2.2 – Fronteiras entre domínios........................................................................................................................................................ 21
Figura 2.3 – Fronteiras virtuais................................................................................................................................................................................ 21
Figura 2.4 – Domínios compostos...................................................................................................................................................................... 22
Figura 2.5- Interpolação linear e bilinear transfinita............................................................................................................... 24
Figura 2.6 – Volume de controlo 2D.............................................................................................................................................................. 30
Figura 2.7 – Matrizes L, U e M do método SIP............................................................................................................................. 46
Figura 3.1 – Malha cartesiana.................................................................................................................................................................................... 58
Figura 3.2 – Isotérmicos da solução numérica e analítica (Olson e Schultz).................................... 58
Figura 3.3 – Comparação entre solução analítica e numérica (domínio rectangular)...........59
Figura 3.4 – Comparação entre solução analítica e numérica (domínio rectangular)...........59
Figura 3.5 – Comparação entre solução analítica e numérica (domínio rectangular)...........60
Figura 3.6 – Variação do parâmetro β na formulaçaõ de Richtmayer e Morton........................60
Figura 3.7 – Malha curvilínea e corresponde malha computacional.............................................................. 63
Figura 3.8 – Isotérmicos da solução analítica e numérica (domínio circular).................................. 64
Figura 3.9 – Comparação entre solução analítica e numérica (domínio circular)......................64
Figura 3.10 – Diferença de temperatura para 2 critérios de paragem distintos............................... 66Figura 4.1 – Balanço térmico num elemento infinitesimal........................................................................................... 69
Figura 4.2 – Variação da fracção sólida com a temperatura....................................................................................... 70
Figura 4.3 – Secção transversal do conjunto peça/moldação..................................................................................... 80
Figura 4.4 – Divisão da secção transversal do conjunto peça/moldação.................................................. 80
Figura 4.5 – Divisão do conjunto peça/ moldação em 17 sub-domínios .................................................. 81
Figura 4.6 – Malha do conjunto peça/moldação............................................................................................................................ 81
Figura 4.7 – Organograma do código desenvolvido................................................................................................................ 84
xi
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 12/134
Figura 4.8 – Campo de temperaturas inicial (1º e 2º abordagens)...................................................................... 86
Figura 4.9 – Campo de temperaturas inicial (3º abordagem)..................................................................................... 86
Figura 4.10 – Campo de temperaturas - conjunto 0.25 s (1º abordagem)............................................... 87
Figura 4.11 – Campo de temperaturas decorridos 2 s (1º abordagem)........................................................ 88
Figura 4.12 – Campo de temperaturas - conjunto 4 s (1º abordagem)......................................................... 88
Figura 4.13 – Campo de temperaturas - conjunto 6 s (1º abordagem)......................................................... 89
Figura 4.14 – Campo de temperaturas da peça 6 s (1º abordagem)................................................................. 89
Figura 4.15 – Campo de temperaturas da moldação 6 s (1º abordagem)................................................. 89
Figura 4.16 – Campo de temperaturas - conjunto 0.25 s (2º abordagem)............................................... 90
Figura 4.17 – Campo de temperaturas - conjunto 4 s (2º abordagem)......................................................... 91
Figura 4.18 – Campo de temperaturas - conjunto 8 s (2º abordagem)......................................................... 92
Figura 4.19 – Campo de temperaturas- conjunto 12 s (2º abordagem)....................................................... 92
Figura 4.20 – Campo de temperaturas- conjunto 14,8 s (2º abordagem)................................................. 93
Figura 4.21 – Campo de temperaturas da peça 14.8 s (2º abordagem)....................................................... 93
Figura 4.22 – Campo de temperaturas da moldação 14.8 s (2º abordagem)....................................... 94
Figura 4.23 – Campo de temperaturas - conjunto 0.25 s (3º abordagem)............................................... 95
Figura 4.24 – Campo de temperaturas - conjunto 4 s (3º abordagem)......................................................... 95Figura 4.25 – Campo de temperaturas - conjunto 8 s (3º abordagem)......................................................... 95
Figura 4.26 – Campo de temperaturas - conjunto 12 s (3º abordagem)..................................................... 96
Figura 4.27 – Campo de temperaturas - conjunto 14.9 s (3º abordagem)............................................... 96
Figura 4.28 – Campo de temperaturas da peça 14.9 s (3º abordagem)....................................................... 96
Figura 4.29 – Campo de temperaturas da moldação 14.9 s (3º abordagem)....................................... 97
Figura 4.30 – Variação da fracção sólida nos diferentes sub-domínios da peça........................... 98
Figura 5.1 – Localização dos termopares na moldação.................................................................................................. 101Figura 5.2 – Comparação entre os resultados experimentais e numéricos........................................ 102
Figura 5.3 – Comparação entre os resultados experimentais e numéricos........................................ 103
Figura 5.4 – Comparação entre os resultados experimentais e numéricos........................................ 106
Figura 5.5 – Comparação entre os resultados experimentais e numéricos........................................ 107
xii
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 13/134
Lista de Tabelas
Quadro 1.1 – Comparação entre métodos: experimental, teórico e numérico..................................... 7
Quadro 3.1 – Áreas e volumes em 2D utilizadas no método dos volumes finitos.....................56
Quadro 3.2- Coeficientes Anb do modelo da condução em domínio rectangular ......................... 56
Quadro 3.3 – Coeficientes A P e Q P (domínio rectangular)............................................................................................ 56
Quadro 3.4 – Coeficientes Anb para o modelo da condução em geometria circular .................62
Quadro 3.5 – Coeficientes A P e Q P (domínio circular)........................................................................................................ 62
Quadro 3.6 – Performance de vários métodos iterativos………………………………………………..........................65
Quadro 4.1 – coeficientes Anb do modelo da solidificação............................................................................................. 75
Quadro 4.2 – Coeficientes A P e Q P (modelo solidificação).......................................................................................... 75
Quadro 4.3 – Propriedades físicas dos materiais da moldação e da peça................................................ 78
Quadro 4.4 – Coeficientes de transferência de calor nas diferentes interfaces............................... 78
Quadro 4.5 – Software comercial......................................................................................................................................................................... 82
Quadro 5.1 – Erro médio das diferentes abordagens do problema da solidificação...........109
xiii
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 14/134
Índice
Agradecimentos............................................................................................................................................................................................................................... ii
Resumo......................................................................................................................................................................................................................................................... iii
Abstract......................................................................................................................................................................................................................................................... v
Résume.......................................................................................................................................................................................................................................................... vi
Lista de símbolos..................................................................................................................................................................................................................... viii
Romanos.............................................................................................................................................................................................................................. viii
Gregos.......................................................................................................................................................................................................................................... ix
Índices.......................................................................................................................................................................................................................................... ix
Lista de figuras............................................................................................................................................................................................................................... xi
Lista de tabelas........................................................................................................................................................................................................................... xiii
Índice........................................................................................................................................................................................................................................................... xiv
Capítulo I – Introdução...................................................................................................................................................................................................... 1
I. Introdução.............................................................................................................................................................................................................................................. 1
II. Motivação........................................................................................................................................................................................................................................... 2
III. Estado da arte.............................................................................................................................................................................................................................. 7
IV. Objectivo....................................................................................................................................................................................................................................... 16
V. Plano....................................................................................................................................................................................................................................................... 16
Capítulo II – Técnicas numéricas para a resolução de equações diferenciais...................18
I. Introdução.......................................................................................................................................................................................................................................... 18
II. Discretização do espaço............................................................................................................................................................................................ 19
II.1. Domínios compostos................................................................................................................................................................................ 20
II.2. Geração da malha por interpolação algébrica.................................................................................................... 23
II.3. Método dos volumes finitos.......................................................................................................................................................... 25
II.3.1. Aproximação de integrais de superfície................................................................................................ 29
xiv
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 15/134
II.3.2. Aproximação de integrais de volume....................................................................................................... 31
II.3.3. Técnicas de interpolação.............................................................................................................................................. 32
II.4. Discretização em problemas transitórios................................................................................................................... 34
II.5. Aproximação numérica de derivadas.............................................................................................................................. 37
II.6. Resolução de sistemas de equações lineares........................................................................................................ 39
II.6.1. Eliminação de Gauss.......................................................................................................................................................... 39
II.6.2. Factorização LU........................................................................................................................................................................ 41
II.6.3. Métodos iterativos.................................................................................................................................................................. 42
II.6.3.1. Método de Jacobi............................................................................................................................................... 46
II.6.3.2. Método de Gauss-Seidel......................................................................................................................... 46
II.6.3.3. Método de Stone................................................................................................................................................. 47
III. Conclusão..................................................................................................................................................................................................................................... 52
Capítulo III – Formulação matemática da transferência de calor....................................................... 54
I. Introdução.......................................................................................................................................................................................................................................... 54
II. Condução de calor em domínio rectangular ............................................................................................................................... 54
III. Condução de calor em domínio circular ........................................................................................................................................ 61
IV. Conclusão..................................................................................................................................................................................................................................... 66
Capítulo IV – Simulação numérica da solidificação em fundição.......................................................... 68
I. Introdução.......................................................................................................................................................................................................................................... 68
II. Modelo matemático......................................................................................................................................................................................................... 68
III. Discretização........................................................................................................................................................................................................................... 74
IV. Condições iniciais e de fronteira............................................................................................................................................................... 76
V. Domínio de cálculo.......................................................................................................................................................................................................... 80VI. Estrutura do código desenvolvido........................................................................................................................................................... 81
VII. Apresentação e discussão dos resultados................................................................................................................................... 86
VII.1 – Campo de temperaturas inicial na moldação............................................................................................. 86
VII.2 – Resultados considerando transferência de calor por condução...................................... 88
VII.3 – Resultados considerando transferência de calor Newtoniana parcial...................91
VII.4 – Resultados considerando transferência de calor Newtoniana........................................... 95
VIII. Conclusão............................................................................................................................................................................................................................... 99
xv
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 16/134
Capítulo V – Validação do modelo numérico....................................................................................................................... 101
I. Introdução...................................................................................................................................................................................................................................... 101
II. Comparação entre dados experimentais e numéricos............................................................................................... 102
III. Conclusão................................................................................................................................................................................................................................. 111
Capítulo VI – Conclusões e perspectivas de trabalhos futuros............................................................... 112
I. Conclusões................................................................................................................................................................................................................................... 112
II. Perspectivas de desenvolvimentos futuros................................................................................................................................ 115
Referências....................................................................................................................................................................................................................................... 116
xvi
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 17/134
Capítulo I
Introdução
I. Introdução
Os processos de produção de peças utilizando a solidificação de metais em cavidades com
forma adequada têm sido aplicados há milénios pelo homem. Como exemplos históricos
podem ser citadas as ferramentas e peças ornamentais fundidas há cerca de 4000 anos pelos
egípcios e assírios, as moedas e obra de arte chinesas de há 3000 anos atrás e as esculturas
gregas de grandes dimensões fundidas há 2500 anos. Naturalmente esses processos foram
desenvolvidos empiricamente através de tentativa e erro, e esse tipo de desenvolvimento no
campo da fundição, persistiu de um modo geral, até há bem pouco tempo. No entanto, acrescente utilização dos processos de fundição na produção de peças de maior precisão e
em maiores quantidades e, sobretudo o emprego cada vez maior da automação nesses
processos, tem exigido o desenvolvimento de métodos de análise mais elaborados, que
levem a um controlo mais preciso dos mesmos.
A técnica da fundição pode ser explicada do seguinte modo: provoca-se a fusão do metal
com que se enche de seguida cavidades, ditas moldações, cuja forma corresponde à forma
negativa do objecto que se pretende obter. O metal ao solidificar conserva a forma da
cavidade, embora a peça que se obtém comporte igualmente alguns apêndices,
nomeadamente: o ou os canais de enchimento, os alimentadores e os canais de evacuação
1
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 18/134
Capítulo I – Introdução
de ar. O metal que solidifica nestes canais será posteriormente separado da peça
propriamente dita, na operação de rebarbagem [1].
Actualmente, o processo de fundição é utilizado em larga escala e coloca ao Homem umnovo desafio – a sua viabilidade económica. Esta pode ser comprometida se não se
projectar correctamente a moldação responsável pela qualidade do produto final.
A solução de problemas reais de engenharia através de técnicas numéricas é actualmente
uma realidade tanto ao nível académico quanto industrial. A crescente evolução dos
computadores tem possibilitado que problemas cada vez mais complexos possam ser
resolvidos através de técnicas numéricas. Outro factor que também contribuiu para esta
tendência está relacionado com os custos de projecto. Os computadores modernos além decada vez mais poderosos são cada vez menos dispendiosos. Hoje é possível que todo o
desenvolvimento de um problema numérico seja criado num microcomputador, cujo custo é
muito baixo, deixando apenas as grandes simulações, com malhas refinadas, para as
estações de trabalho. Alguns problemas podem até mesmo serem resolvidos no próprio
microcomputador em que o código foi desenvolvido.
Ainda com relação a factores económicos, é importante salientar que actualmente já é
possível substituir horas de experimentação em laboratório a custos altíssimos por
simulação em computador, diminuindo enormemente os custos de projecto, e deixando os
testes de laboratório apenas para refinamentos, ou para a modelação de problemas que
ainda não possuam uma formulação matemática satisfatória.
II. Motivação
O fenómeno físico da solidificação constitui a principal vantagem da fundição, porque o
metal líquido permite que com pouco esforço seja possível dar a forma desejada aos metais
e é amplamente determinada quer por factores dinâmicos quer por factores termodinâmicos
[2].
2
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 19/134
Capítulo I – Introdução
O aquecimento de um metal sólido provoca o aumento da amplitude de vibração dos
átomos seus constituintes, até que se dá o colapso da estrutura sólida. Em condições
isobáricas, este colapso ocorre a uma temperatura específica para cada elemento – a
temperatura de fusão. A esta temperatura, a massa do elemento passa do estado sólido aoestado líquido à medida que se vão quebrando as ligações com os vizinhos, características
do estado sólido. A energia necessária para provocar a inviabilidade destas ligações
constitui uma propriedade intrínseca de cada elemento, designada por calor latente.
A transferência de calor entre um corpo a temperatura elevada e o meio em que está
inserido realiza-se a uma velocidade que depende das características de ambos, de forma
que a temperatura do corpo tende para a temperatura do meio em que se encontra imerso
(figura 1.1a).
b a c
T e m p e r a t u r a
Tempo
L í q u i d o
Líquido+
Sólido
Tempo
T e m p e r a t u r a
S ó l i d o
Líquido+
Sólido L í q u i d o
T e m p e r a t u r a
Tempo
S ó l i d o
Figura 1.1 – Curvas de arrefecimento mostrando a evolução da temperatura de um corpo metálico aarrefecer em condições de equilíbrio, em três casos diferentes:
a) Se não ocorrer mudança de estado durante o arrefecimento.
b) Ocorrência de mudança de estado que produz um patamar à temperatura a que se dá a transição,no caso de se tratar de um metal puro.
c) No caso de uma liga metálica, a ocorrência de mudança de estado pode produzir não um patamar
mas apenas duas inflexões: a primeira indicando o início, solidus, e a segunda o fim, liquidus, da
solidificação.
No caso de um metal puro em arrefecimento a partir do estado líquido em condições de
equilíbrio, ocorre solidificação da massa metálica a temperatura constante. A curva de
arrefecimento de um corpo nestas condições apresenta então um patamar como indicado
pela figura 1.1b, à temperatura a que líquido e sólido coexistem, o que significa que o metal
solidifica a uma velocidade que corresponde à capacidade de absorção do calor latente pela
3
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 20/134
Capítulo I – Introdução
sua vizinhança. A temperatura do metal só baixará quando toda a massa metálica tiver
passado ao estado sólido.
No caso do corpo ser constituído por uma liga metálica, ao contrário do que acontece paraum metal puro, a solidificação não ocorre geralmente a temperatura constante, porque os
cristais que se formam não têm a mesma composição que o líquido donde são originários
[1]. De forma a equilibrar o desvio de composição de cristais formados, a composição do
líquido altera-se ao longo do processo, e, como líquidos de composição diferente podem
iniciar a solidificação a temperaturas distintas (figura 1.1c), a temperatura de solidificação
vai baixando à medida que a composição do líquido se vai alterando.
Durante a solidificação, o metal sofre uma variação de volume designada contracção desolidificação, e corresponde à assunção de uma estrutura mais organizada e compacta. A
solidificação progride a partir das paredes da moldação, mantendo-se o metal no estado
líquido durante mais tempo no interior das paredes maciças da peça. Se a liga tiver um
intervalo de solidificação muito estreito, a solidificação pode estar efectivamente terminada
nas partes finas e a contracção, na parte central da peça que já não é alimentada, dar origem
à formação de vazios na peça, designados por rechupes (figuras 1.2 e 1.3).
a b c
Figura 1.2 - Mecanismo de formação de rechupes como resultado da variação de volume mássico durante o
processo de solidificação.
a) A massa de metal solidifica sem que o avanço da solidificação seja observável do exterior;b) O avanço da solidificação provoca o abaixamento do nível do líquido;
c) No final da solidificação o rechupe não é visível do exterior.
As variações bruscas de espessura provocam ainda descontinuidade na contracção no
estado sólido. Desta situação resulta, sobretudo em moldação metálica, um campo de
tensões residuais que pode provocar distorções ou mesmo rupturas e fendas na peça.
4
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 21/134
Capítulo I – Introdução
F igura 1.3 - Zonas de formação de rechupes numa peça.
A condutibilidade térmica do metal constituinte da moldação é normalmente elevada e da
ordem de grandeza da do metal vazado, pelo que a sua solidificação é muito rápida.
Quando o vazamento do metal líquido é feito por gravidade, o metal pode solidificar em
percentagens consideráveis antes do término do enchimento. Uma forma de acautelar o
interromper do enchimento por solidificação extemporânea é fazer o vazamento com a
moldação quente, mas à mais baixa temperatura compatível com o bom enchimento da
moldação. Desta forma, a contracção do metal na fase líquida é reduzida ao mínimo, e a
moldação aquecida permite um trajecto mais longo do metal antes que aconteça qualquer
interrupção do fluxo por estrangulamento do escoamento.
O ciclo de produção de uma peça por vazamento em moldação metálica começa com a
moldação fechada pronta para receber o metal líquido e termina com a moldação fechada
pronta para o próximo vazamento. Para garantir a reprodutibilidade, é conveniente que a
moldação se encontre em condições térmicas idênticas, antes de cada vazamento.
Se a cadência de produção baixar, a temperatura da moldação irá baixar e alterar as
condições de solidificação. Inversamente, se a cadência aumenta, a temperatura da
moldação aumenta, criando-se o efeito oposto. Para se variar a cadência de produção terá
que se intervir nas características de escoamento de calor, para e na moldação.
5
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 22/134
Capítulo I – Introdução
A moldação metálica tem, portanto, uma primordial importância na transferência de calor.
Mas o campo de temperaturas na moldação também não é uniforme, e varia com o tempo.
A sua variação é devida ao arrefecimento progressivo do metal na sua cavidade, à
contracção da peça no estado sólido e também às irregularidades de espessura da peça e da própria moldação, que afectam a difusão do calor.
O comportamento térmico do conjunto é ainda afectado pelo deslocamento entre a peça e a
moldação, pela utilização de pinturas de protecção e isolamento térmico e pelo
aquecimento ou arrefecimento propositado de zonas específicas da moldação.
Tradicionalmente, os problemas físicos são resolvidos quer por métodos analíticos quer por
métodos experimentais, entretanto, a constante evolução do computador tornou viável umterceiro método – a aproximação numérica. No entanto, o método experimental continua a
ter grande importância, principalmente quando os fenómenos envolvidos são de grande
complexidade, mas a tendência para o uso dos métodos numéricos tem vindo a aumentar
devido à fiabilidade dos resultados e baixo custo quando comparados com os métodos
experimentais. No quadro 1.1 faz-se uma discriminação das vantagens e desvantagens
destes três métodos.
A sugestão aqui dada não é fazer crer que o método numérico substitua totalmente os testes
experimentais como meio de reunir informação para um determinado projecto, mas uma
certeza podemos ter, os métodos computacionais serão utilizados com mais frequência no
futuro. Naturalmente certos fenómenos físicos possuem uma complexidade de tal ordem
que se continuará a necessitar da execução de determinados testes de confirmação
experimental, mas mesmo nestes casos, a simulação numérica pode ser usada no sentido de
reduzir o número de condições a experimentar.
A abordagem analítica do problema encontra grandes dificuldades devido à complexidade
do fenómeno da solidificação. O método experimental tem como principal óbice os seus
custos que é, precisamente, o que se pretende reduzir. A alternativa que os métodos
numéricos proporcionam, encontra também dificuldades, que se prendem com a
6
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 23/134
Capítulo I – Introdução
multiplicidade e complexidade de geometrias obtidas por fundição, e com a disponibilidade
e fiabilidade de dados físicos de caracterização dos materiais e fenómenos envolvidos [1].
Quadro 1.1 – Comparação entre métodos: experimental, teórico e numérico.
MÉTODO VANTAGENS DESVANTAGENS
Equipamento necessário
Custos de operação
Dificuldades de mediçãoExperimental Realista
Problemas de escala
Restrito a geometrias simples
Teórico
Informação global do problema,
geralmente em termos de
expressões matemáticas Restrito, em geral, aos problemas lineares
Capacidade para tratar problemas
não lineares
Erros de truncatura
Capacidade para tratar problemas
físicos complexos
Problemas em estabelecer as condições de
fronteira Numérico
Permite obter a evolução temporal
da variável dependenteCustos de computação
A principal razão pela qual se recorre, cada vez mais, à simulação computacional das
transformações que ocorrem no metal ao arrefecer, de modo a detectar possíveis erros de
concepção da moldação é essencialmente a viabilidade económica da produção.
III. Estado da arte
De um modo geral, apesar de haver alguns trabalhos anteriores, o início de estudos
sistemáticos de análise do fenómeno da solidificação deu-se na primeira metade do século
7
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 24/134
Capítulo I – Introdução
XX, e a sua frequência tem aumentado desde então. Esses trabalhos têm visado, através de
uma análise teórico/experimental dos processos de solidificação, obter com maior precisão,
parâmetros que actuam efectivamente na transformação líquido/sólido, com o objectivo de
se exercer um maior controlo sobre a estrutura interna e, consequentemente, sobre osdefeitos e propriedades das peças obtidas através desses processos.
É sabido que a simulação numérica é já uma ferramenta importante de suporte dos
processos industriais de fundição. O objectivo da simulação é retirar o máximo de
informação da dinâmica da solidificação. Como já foi referido na secção anterior o
arrefecimento e consequente mudança de fase do metal líquido, é um fenómeno complexo,
não apenas devido aos modos de transferência de calor existentes, mas também devido à
geometria das peças, normalmente complexa.
Kothe et al , refere que num cenário típico de enchimento duma moldação, estão presentes
pelo menos quatro meios separados por interfaces: o material da cavidade (moldação); o ar,
ou um gás inerte; o metal líquido de enchimento da moldação; e o metal solidificado,
subsequente ao enchimento da moldação à medida que o arrefecimento progride [4].
A transferência de calor em fundição está longe de ser um problema estacionário, até pelo
carácter discreto da forma de produção desta tecnologia. O líquido vazado na moldação
arrefece em contacto com as paredes desta. Entretanto, há diferentes obstáculos entre o
metal líquido e o ambiente para o qual é escoado o calor: a interface moldação/ambiente, a
moldação propriamente dita; a interface metal/moldação, o metal sólido, a interface metal
sólido/metal líquido e o metal líquido [1].
O problema é complexo, em particular para a maioria das geometrias que interessam em
fundição, pelo que a sua descrição matemática, quando viável, pode envolver expressões degrande complexidade.
O fenómeno da solidificação em fundição é descrito através da equação da difusão de calor,
cuja forma diferencial é a seguinte [1, 5:7]:
8
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 25/134
Capítulo I – Introdução
( )( ).
pC k q
t
ρ φφ
∂= ∇ ∇ +
∂& (1.1)
ou usando o muito popular modelo da entalpia [8:14],
( )( ).
H k q
t
ρφ
∂= ∇ ∇ +
∂& (1.2)
onde representa o calor gerado no elemento, H a entalpia e ρ, C q& p e k são características
físicas do material: massa volúmica, capacidade térmica mássica e condutibilidade térmica,
respectivamente. Uma outra forma desta equação é fornecida por Knoll et al . onde o índice
s denota fase sólida e o índice l denota fase líquida, é o vector velocidade [13,14]. l ur
( )( )( ) .l l l
H h u k
t
ρρ φ
∂+ ∇ − ∇ ∇ =
∂r
0 (1.3)
A massa volúmica da mistura é definida por:
l s s s f f ρρρ )1( −+= (1.4)
e a entalpia da mistura por,
l l s s s s h f h f H ρρρ )1( −+= (1.5)
e as entalpias de fase são dadas pelas seguintes expressões:
, s p sh C φ= (1.6)
,l p l h C hf φ= + ∆ (1.7)
onde f s representa a fracção sólida e ∆hf o calor latente de fusão.
9
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 26/134
Capítulo I – Introdução
De forma a facilitar o desenvolvimento dos modelos numéricos para a descrição do
fenómeno da solidificação as características físicas ρ , C p e k do metal e da moldação são
assumidas como constantes [1,5:14].
Como se sabe a mudança de fase de qualquer metal ou liga metálica não é instantânea. Esta
mudança de fase ocorre a uma taxa que depende do equilíbrio de fases do material a
solidificar.
No caso de um metal puro, ou de uma liga de composição eutética, em arrefecimento a
partir do estado líquido em condições de equilíbrio, ocorre solidificação da massa metálica
a temperatura constante. A curva de arrefecimento de um corpo nestas condições apresenta
então um patamar como indicado pela figura 1.1b, à temperatura a que líquido e sólido
coexistem, o que significa que o metal solidifica a uma velocidade que corresponde à
capacidade de absorção do calor latente pela sua vizinhança. Este tipo de problema é
denominado na literatura pelo problema de Stefan [7,13:17], cuja solução analítica pode ser
consultada em [17].
Para uma liga metálica, a solidificação não ocorre em geral a temperatura constante,
(figura 1.1c). A taxa de arrefecimento pode ser expressa em função da fracção efectiva de
material sólido, a fracção sólida f s, da massa volúmica do material ρ, e da variação da
entalpia na transformação de fase ∆h f , (calor latente) [1]. Juric e Tryggvason [7] e mais
recentemente Voller [8] usaram a função δ(x-x f ) que é não nula apenas na interface
sólido/líquido, isto é, quando x = x f , para descrever o termo fonte da equação da energia.
À temperatura constante φm a fracção sólida é uma função da entalpia de acordo com a
seguinte relação.
(1 ) p m s H C f hf φ− = − ∆ (1.8)
A fracção sólida pode então ser obtida através da equação da energia igualando o termo
transitório da temperatura a zero dado que a temperatura é constante e igual à temperatura
10
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 27/134
Capítulo I – Introdução
de fusão. A temperatura pode ser calculada como uma função da entalpia T=τ(Η), da
seguinte forma:
( ) 0)( =
∂∂
∂∂−
∂∂
x H
xk
t H τ
ρ (1.9)
Para um metal puro a função τ(Η), é definida da seguinte forma:
, ;
( ) ,
( ) / .
p p m
m p m p
f p p m f
H C H C
H C H
H h C H C h
φ
τ φ φ φ
φ
<
= ≤ ≤ − ∆ > + ∆
;m f C h+ ∆ (1.10)
que contempla a relação existente entre a temperatura e a entalpia para um metal puro.
Para o caso de uma liga binária que solidifique num intervalo de temperaturas, a função
τ(Η), é um pouco mais complexa. Assumindo uma relação linear entre as propriedades
temperatura e entalpia vem que:
l l f C mT T += (1.11)
onde T f é a temperatura de fusão para um metal puro, ml é o declive da linha de liquidus e
C l é a concentração de metal líquido. Assumindo equilíbrio termodinâmico na interface
sólido-líquido resulta na seguinte relação para a concentração de metal solidificado na
interface,
C s =γ C l (1.12)
onde γ é o coeficiente de repartição de material sólido e líquido na interface. Usando um
modelo à escala local (zona bifásica) C l é definido como uma função da fracção sólida.
11
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 28/134
Capítulo I – Introdução
C l = G (f s ) (1.13)
Neste estudo foram utilizados dois modelos locais: no primeiro assume-se completa difusão
do soluto no líquido e difusão nula do mesmo no sólido (assunção de Scheil), o segundo
modelo é sustentado pela regra da alavanca no qual é considerado difusão completa do
soluto em ambas as fases. Para uma liga binária estas assunções resultam nas seguintes
expressões para a fracção sólida:
(1 )
1.0 f l
s
f
f
γφ φ
φ φ
− − −
= − − (1.14)
para a assunção de Scheil, e
11.0
1 ) f l
s
f
f φ φ
φ φ γ
− = − − −
(1.15)
para a regra da alavanca.
A função da temperatura para uma liga binária com transformação eutéctica pode, então ser apresentada como,
,
,
, ;
, (( )
( (1 ) ) / , (1 )
.( ) / , ;
p p eut
eut p eut p eut s eut f
s f p p eut s eut f p l f
f p p l f
H C H C
C H C f h H
H f h C C f h H C h
H h C H C h
φ
φ φ φτ
φ φ
φ
<
≤ ≤ + − ∆=
− − ∆ + − ∆ ≤ ≤ + ∆ − ∆ > + ∆
1 ) ;
; (1.16)
onde φeut é a temperatura do ponto eutéctico.
Como é sabido grande parte dos fenómenos físicos são regidos matematicamente por
equações ou sistemas de equações diferenciais que devem ser respeitadas em todo o domínio
e por um conjunto de condições de fronteira e/ou condições iniciais eventualmente também
12
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 29/134
Capítulo I – Introdução
descrito por um sistema de equações diferenciais. Quando um esquema numérico é aplicado
a estas equações diferenciais, o domínio computacional é subdividido em pontos, volumes
ou elementos, pertencentes a uma malha e as equações são discretizadas e resolvidas em
função dessa malha.
A resolução analítica de uma equação de diferenciais parciais só é possível em casos
particulares. Daí a necessidade de reformulação do problema, dando-lhe uma forma
puramente algébrica. Esta operação, designada por discretização, consiste na substituição do
problema contínuo por um problema discreto aproximado.
Para a discretização do espaço, Monteiro [1] utilizou o método das diferenças finitas em
coordenadas generalizadas para poder modelar as geometrias complexas normalmente presentes em fundição. Para a discretização do tempo usou a formulação explícita de Euler.
A definição das geometrias complexas proposta passa pela sua divisão em troços de
geometria mais simples, pelo estabelecimento de malhas curvilíneas em cada troço, pela
definição das interacções a respeitar entre os diferentes troços, o que constitui a
discretização da geometria em causa.
Gong et al. [9,10], Comini et al. [12], Sarler e Alujevic [15] utilizaram o método doselementos finitos como método de discretização, tendo, este último, sujeitado o seu modelo
matemático a vários tipos de condições de fronteira distintos:
o Dirichlet- onde a variável dependente assume um valor específico, por exemplo,
uma fronteira isotérmica [3];
o Neumann - derivada de valor específico da variável dependente, por exemplo, uma
fronteira adiabática [3];
o Robbin – trata-se de uma combinação dos dois tipos de condições de fronteira
anteriores, por exemplo, um fluxo convectivo [3].
As incógnitas H (entalpia) e φ (temperatura) são discretizadas da seguinte forma:
13
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 30/134
Capítulo I – Introdução
∑=
=n
i
ii t H z y x N H 1
)(),,( (1.17)
1
( , , ) ( )n
i ii
N x y z t φ φ=
=
∑ (1.18)
onde N i são as funções de forma.
As n equações podem então ser escritas na forma matricial seguinte,
[M] {V H } + [K] {T} ={F} (1.19)
onde:
dV N N M e
V jiij ∑∫ = ρ (1.20)
∑∫ ∑∫ +
∂
∂
∂
∂+
∂
∂
∂
∂+
∂
∂
∂
∂=
ee S ji
V
ji ji jiij dS N hN dV
z
N
z
N
y
N
y
N
x
N
x
N k K (1.21)
dS T N hdV N Q F ee S
f iV
ii ∑∫ ∑∫ += (1.22)
A derivada em ordem ao tempo é discretizada utilizando uma formulação implícita com
três níveis de tempo.
[ ]{ } [ ]{ } { } F T T T K V M nnn
H =+−+++ −+ 11 )25.0()5.0( βββ (1.23)
A equação resultante da discretização é não linear o que motivou o desenvolvimento de um
novo método de quadratura com apenas um ponto. Este método aplicado aos problemas não
lineares reduz em quase dois terços o tempo de computação e com precisão próxima dos
esquemas de quadratura de quatro pontos.
14
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 31/134
Capítulo I – Introdução
Por fim o método dos volumes finitos utilizado por, Voler [8], Knoll et al. [13,14],
Ambrosi e Preziosi [18,19] e Fabbri [20]. As derivadas de φ no domínio computacional são
aproximadas para a face Este do volume de controlo da seguinte forma:
( ) E W
φφ φ
ξ∂
≈ −∂
(1.24)
( )1
4ne se N NE S SE
φφ φ φ φ φ φ
η∂
≈ − ≈ + − −∂
(1.25)
onde foram assumidos espaçamentos computacionais unitários (∆ξ = ∆η=1).
Uma forma de evitar as restrições em termos de passo de tempo a utilizar (estabilidade), éutilizar formulações implícitas para a discretização temporal em detrimento das explícitas.
Quando em presença de geometrias complexas o uso de coordenadas curvilíneas permite
descrever através de um reduzido número de volumes de controlo os seus contornos, no
entanto, aumenta a complexidade matemática do modelo. A molécula computacional passa
a ser de nove pontos em contraste com os 5 pontos das malhas cartesianas. È, no entanto,
possível utilizar os métodos de resolução de equações lineares desenvolvidos para as
malhas cartesianas escrevendo-se o sistema de equações em duas parcelas:
A B Qφ φ + = (1.26)
onde a matriz A contém os coeficientes ortogonais e a matiz B os coeficientes não
ortogonais.
O sistema de equações (1.26) é, então, resolvido utilizando o chamado deferred correction
scheme, seguinte [21]:
old A Q Bφ φ = − (1.27)
15
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 32/134
Capítulo I – Introdução
onde é a solução da iteração anterior.old φ
Como se pode constatar podem ser utilizados diferentes métodos de discretização para
discretizar as equações que regem o fenómeno da solidificação. No entanto, uma análiseatenta de cada situação poderá ser determinante para a escolha do método que permitirá
obter melhores resultados ou maior comodidade de programação.
IV. Objectivo
O objectivo principal desta tese consiste no desenvolvimento de um software em linguagem
FORTRAN para simular numericamente a solidificação em moldação metálica, usando o
método dos volumes finitos e coordenadas curvilíneas generalizadas.
V. Plano
Parece neste momento perfeitamente claro que o ponto de partida na realização deste
trabalho terá de ser a análise dos fenómenos de transferência de calor que ocorrem no
escoamento de calor desde o metal em arrefecimento até ao meio ambiente, e por este facto
se dedicou as secções anteriores à análise das transformações que ocorrem nos metais e
ligas metálicas ao arrefecerem desde o estado líquido. Na secção anterior procurou-se dar
uma ideia, embora sucinta, do trabalho realizado por alguns investigadores na área da
modelação numérica do fenómeno da solidificação dos metais, onde ficou bem patente a
necessidade de encontrar uma formulação matemática capaz de descrever o fenómeno de
forma clara.
Como os fenómenos físicos da natureza são regidos por equações diferenciais, e só em
casos excepcionais estas se pode resolver analiticamente, a alternativa passa pelo uso dos
métodos numéricos sobre os quais se debruçará o segundo capítulo.
16
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 33/134
Capítulo I – Introdução
Com o intuito de testar o modelo numérico desenvolvido, este vai ser comparado com
soluções analíticas conhecidas. Numa primeira abordagem será comparada a solução
analítica proposta por Olson e Schultz para a condução de calor num corpo bidimensional
utilizando coordenadas cartesianas. Numa segunda abordagem será efectuado o mesmo tipode comparação, mas, agora, utilizando um domínio não rectangular de modo a testar o
código desenvolvido em coordenadas curvilíneas. Este será o assunto do terceiro capítulo.
Devido ao facto das peças obtidas por fundição serem na sua maioria de geometria
complexa, torna-se necessário recorrer a uma transformação de coordenadas que
aumentando a complexidade do problema do ponto de vista matemático, facilita a
programação a realizar em linguagem FORTRAN, mormente porque para diferentes
geometrias o domínio computacional será o mesmo. Este constitui um dos assuntos do
quarto capítulo juntamente com a aplicação do código desenvolvido para uma geometria
complexa concreta.
O modelo numérico a desenvolver terá necessariamente de ser validado. Essa validação
será efectuada no quinto capítulo onde se faz a comparação com resultados experimentais e,
também, com resultados numéricos utilizando o método das diferenças finitas.
No sexto e último capítulo, são apresentadas as principais conclusões deste trabalho e,
ainda, perspectivas de desenvolvimentos futuros.
17
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 34/134
Capítulo II
Técnicas numéricas de resolução de equações diferenciais
I. Introdução
Grande parte dos fenómenos físicos são regidos matematicamente por equações ou
sistemas de equações diferenciais que devem ser respeitadas em todo o domínio e por um
conjunto de condições de fronteira e/ou condições iniciais eventualmente também
descrito por um sistema de equações diferenciais. A resolução analítica de uma equação de diferenciais parciais só é possível em casos
particulares. Daí a necessidade de reformulação do problema, dando-lhe uma forma
puramente algébrica. Esta operação, designada por discretização, consiste na substituição
do problema contínuo por um problema discreto aproximado.
A operação de discretização permite, portanto, a conversão de um modelo contínuo num
modelo discreto aproximado. Como consequência deste processo, o modelo matemático
contendo equações com derivadas parciais pode em geral ser convertido num outro
contendo equações diferenciais ordinárias ou apenas equações algébricas.
Para se recorrer aos métodos numéricos é também necessário representar adequadamente
a geometria, o que se consegue fazendo a escolha de um número finito de nós
computacionais, supostos representativos de todo o domínio, onde irá ser determinada a
solução numérica.
Quando os domínios, físico e computacional, forem idênticos, por exemplo um
rectângulo, não é necessário qualquer transformação, mas no caso de geometrias físicas
18
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 35/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
complexas, a transformação pode provocar modificações drásticas quer nas equações
constitutivas do fenómeno físico quer nas condições de fronteira associadas.
Para solucionar este problema podemos recorrer às técnicas de transformação decoordenadas, relativamente comuns, mas nem sempre é possível descrever
convenientemente um domínio irregular através desta técnica.
A descrição dos contornos das formas geométricas complexas constitui uma grande
dificuldade se pretendermos utilizar malhas cartesianas ortogonais. Esta dificuldade não
se coloca se definirmos malhas de forma livre, cuja utilização facilita a discretização
computacional e possibilita a aplicação de algoritmos genéricos a geometrias complexa.
Uma das formas de conseguir o que acima se disse é dividir a geometria complexa em
troços de geometria mais simples, pelo estabelecimento de malhas curvilíneas em cada
um desses troços e pela definição adequada da interacção entre os mesmos.
II. Discretização do espaço
A definição de coordenadas curvilíneas de forma generalizada, cada vez mais presentes
nos métodos numéricos, permite a transformação de um domínio de forma arbitrária
num domínio de forma rectangular [22:25]. Deste modo, para estabelecer um sistema de
coordenadas curvilíneas apenas se terá que definir as relações funcionais ξ(x,y) e η(x,y)
que aplicam um domínio geométrico qualquer num rectângulo.
Na figura 2.1 representa-se um domínio de forma arbitrária, bem como a definição dos
sentidos de variação das coordenadas curvilíneas.
Um ponto qualquer do interior do domínio será definido pelas coordenadas curvilíneas
(ξ,η), estabelecidas em função da sua posição real e das regras impostas para a
construção das linhas de coordenada invariante.
Os pontos no interior ficarão localizados pela intersecção de linhas coordenadas que
ligam os pontos correspondentes em fronteiras opostas [23].
19
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 36/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
J
ξ
(a) (b)
ξ I 0
0
η
η
(x,y) η
(x,y) ξ
(c)
(x,y) η
(x,y) ξ
(x,y) η (x,y) ξ
ξ
η
x
Figura 2.1 – Discretização de uma geometria de forma arbitrária.
Especificadas as coordenadas curvilíneas sobre as fronteiras (figura 2.1a), o passo
seguinte é escrever de forma adequada o domínio no interior do sistema (figura 2.1b).
Ao estabelecer uma relação biunívoca entre a posição de todos os pontos assim
definidos, e os pontos de uma malha arbitrariamente rectangular, o domínio regular que
fica disponível é a imagem computacional (figura 2.1c) da geometria real, que permite a
transposição do problema real e o desenvolvimento da solução neste domínio.
II.1. Domínios Compostos
No caso do estudo de peças para fundição, estão envolvidas geralmente regiões que,além de apresentarem formas geométricas distintas, também têm características físicas
diferentes. Nestes casos será necessário ter em atenção alguns aspectos relativos ao
contacto entre as diferentes sub-regiões.
Em primeiro lugar, considere-se que o domínio circular da figura 2.2 constitui a secção
recta de uma barra cilíndrica. A obtenção da peça exige a existência de uma moldação,
cuja secção recta poderá ser a coroa circular limitada pelas circunferências Γ1
e Γ2,
20
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 37/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
representada na figura 2.2 A sua transformação num rectângulo pode ser conseguida
seccionando o domínio físico pelo segmento AB, assim criando as fronteiras virtuais AB
e CD, de acordo com o procedimento ilustrado na figura 2.3.
1 Γ
1Γ2Γ
(a) (b)
Figura 2.2 – A curva Γ1 que define o contorno exterior da secção transversal de uma barra cilíndrica
(a), define também o contorno interior da parede da moldação que permite a sua produção (b), e
estabelece a fronteira entre ambos os domínios. A curva Γ2 que define o contorno exterior da
secção transversal da parede da moldação estabelece a fronteira com o ambiente.
A
(a) (b)
B
D
C
(c)
A B
D
C A B
Figura 2.3 - A transformação do sub-domínio constituído pela moldação num domínio rectangular pode
ser feita pela definição de uma fronteira virtual AB ligando as fronteiras constituídas por Γ1 e Γ2.
A fronteira virtual AB desdobra-se em duas fronteiras (b) para definir o domínio computacional (c).
O domínio de cálculo básico para este caso poderá ser definido criando uma sequência
de dois meios que podem apresentar, e normalmente apresentam, diferentes
características físicas, com fronteira para o exterior e com interface de contacto entre si.
21
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 38/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
22
Na figura 2.4 mostra-se o aspecto computacional que o conjunto peça-moldação pode
apresentar, e as conexões que é necessário estabelecer entre as fronteiras reais ou
virtuais dos domínios definidos.
Figura 2.4 - Conjunto peça-moldação. Geração de malha no metal (a), criação de quatro fronteirasvirtuais no domínio constituído pela moldação (b). No domínio computacional, (c), cada
sub-domínio gerado apresenta a forma rectangular.
Mesmo que o domínio computacional não constitua um polígono de quatro lados, pode
ser dividido em blocos desse tipo, como se mostra na figura 2.4.
O domínio computacional compreenderá uma colecção de blocos rectangulares, em que
a identificação dos pontos pode ser feita utilizando valores inteiros sucessivos para cada
linha coordenada de cada bloco. Assim as coordenadas de cada ponto serão r(i,j), com
(i,j) indicando a posição (ξ,η) no domínio transformado [23].
As fronteiras no domínio computacional podem corresponder quer a fronteiras físicas
reais, quer a cortes do domínio físico homogéneo, realizados por conveniência de
procedimento.
A
(a) (b)
B
D
C
(c)
G
E
F
C
G
E F
A
B
D
C
G
E
F
A B
D
H G
E F
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 39/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
II.2. Geração da malha por interpolação algébrica
A determinação das posições nodais no interior de um domínio pode ser conseguida
interpolando algebricamente valores definidos numa fronteira para o interior dodomínio.
A primeira fase da colocação de pontos no domínio passa pelo estabelecimento de
pontos nas suas fronteiras, de acordo com uma distribuição adequada.
O espaçamento dos pontos sobre a curva não tem que ser necessariamente uniforme.
Pode, em determinadas situações, ser vantajoso que haja espaçamentos diferentes em
extremidades ou em troços diferentes da mesma curva. Uma das vantagens desta
situação pode ser a redução global do número de nós, conseguida com o aumento do
espaçamento em zonas cuja descrição pode ser feita com menor densidade de pontos.
No entanto, é necessário ter em atenção que variações bruscas desse espaçamento
podem inviabilizar o cálculo da derivada local, devido à condição de continuidade [23].
Após a definição das quatro fronteiras que definem o contorno de um domínio, faz-se a
colocação dos nós em cada uma, de forma a garantir o mesmo número de nós ( n para adirecção e m para a direcção η) em fronteiras opostas, e também a distribuir
adequadamente a posição dos mesmos em cada fronteira.
ξ
Uma interpolação algébrica de Lagrange, pode ser aplicada a todas as linhas
coordenadas. Numa primeira aproximação pode fazer-se a interpolação ao longo de ξ j
entre os pares de valores que as funções xk (ξi,ξ j ) tomam em fronteiras ξi opostas. Uma
interpolação destes valores pode ser feita recorrendo a uma função linear do tipo,
( )1
1
j j f
n
ξξ
−=
− (2.1)
resultando na malha que se apresenta na figura 2.5a, e obtém-se a seguinte expressão
para a função xk (ξi,ξ j ):
( ) ( ) ( ) ( )( ) ( )1
, ,1 1
i j j i j i
k k x f x f xξ ξ ξ ξ ξ ξ= + − ,k n (2.2)
23
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 40/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
Desta forma, as posições determinadas para os extremos em que ξi=1 e ξi=m só
coincidem com posições sobre as fronteiras ξ j efectivas, ξ j=1 e ξ j=n, se estas forem
rectilíneas no domínio real.
A diferença entre estas posições pode ser calculada para cada nó destas fronteiras como
segue:
j
k
j
k
j
k x x ξξξ ,1,1,1 12 −= x (2.3)
( ) ( ) ( ) j
k
j
k
j
k m xm xm x ξξξ ,,, 12 −= (2.4)
Estas diferenças podem ser interpoladas para todo o domínio, obtendo-se, para cada ponto do interior, o valor da diferença entre as posições determinadas e as desejadas
[23]:
( ) ( ) ( ) ( )( ) ( )2 2, 1, 1i j i j i j
k k x f x f x mξ ξ ξ ξ ξ ξ= + − 2 ,k (2.5)
(a) (b)
Figura 2.5 - Interpolação linear (a), Interpolação bilinear transfinita (b).
Para se conseguir o objectivo referido, basta então subtrair às posições determinadas na
interpolação linear, as diferenças interpoladas para cada ponto interior, para se obter umconjunto de posições que respeitam em todas as fronteiras a geometria real do domínio
(figura 2.5b).
( ) ( ) ( ) ji
k
ji
k
ji
k x x x ξξξξξξ ,,, 21 −= (2.6)
24
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 41/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
II.3. Método dos volumes finitos
Tal como foi referido no primeiro capítulo qualquer método de discretização pode ser
utilizado para discretizar a equação que rege o fenómeno da solidificação; no entanto,uma análise atenta de cada situação poderá ser determinante para a escolha do método
que permitirá obter resultados ou mais comodidade de programação.
i) O método das diferenças finitas tem como ponto de partida a escrita das equações
constitutivas do fenómeno físico na forma diferencial. O domínio de cálculo é
representado por uma malha de pontos. Em cada ponto da malha, as equações
diferenciais são aproximadas por substituição das derivadas por termos que são função
dos valores nodais da variável dependente que pretendemos conhecer. O resultado é uma
equação algébrica para cada nó da malha onde aparecem como incógnita, o valor da
variável dependente nesse ponto e num determinado número de nós vizinhos.
A aproximação da primeira e segunda derivada são realizadas com base na expansão em
série de Taylor das variáveis em relação às coordenadas. Quando necessário, estes
métodos são também usados para determinar o valor das variáveis noutros locais além
dos pontos da malha, o que acontece no caso da interpolação.
Em malhas regulares, o método das diferenças finitas é de implementação simples e
bastante preciso.
As desvantagens deste método residem no facto das leis de conservação não serem
respeitadas a não ser que se tenham precauções especiais. Também o facto das restrições
às geometrias simples é uma grande desvantagem tendo em conta que a grande parte dos
problemas de interesse prático não ocorrem nesse tipo de geometria.
ii) No método dos elementos finitos o domínio de cálculo é dividido num determinado
número de elementos finitos geralmente irregulares; em 2D, assumem formas como
triângulos ou quadrados, enquanto que em 3D as formas mais usadas são os cubos e
hexágonos. O que distingue este método de todos os outros é o facto das equações serem
multiplicadas por uma função peso antes de se realizar a integração sobre todo o
domínio. No caso mais simples do método dos elementos finitos, a solução é aproximada
25
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 42/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
por uma função de forma linear e aplicada a todos os elementos de forma a garantir a
continuidade da solução nas suas fronteiras. Uma função deste tipo pode ser obtida
através dos valores que toma nos vértices dos elementos.
Esta aproximação é então substituída no integral da lei da conservação já multiplicada
pela função peso, sendo a equação a ser resolvida derivada com a assunção de que a
derivada em relação a cada valor nodal é igual a zero; o que corresponde a seleccionar a
melhor solução entre as possíveis (aquela que tiver menor resíduo). O que resulta num
sistema de equações não linear.
Uma importante vantagem do método dos elementos finitos é a sua capacidade de lidar
com geometrias de forma arbitrária. As malhas podem ser facilmente refinadas por
divisão simples dos elementos.
A principal desvantagem deste método é aquela que é partilhada por todos os métodos
de discretização quando se utilizam malhas irregulares, e consiste no facto das matrizes
das equações linearizadas não possuírem a mesma forma que as resultantes da utilização
de malhas regulares, o que dificulta a escolha de um método de resolução eficiente a
utilizar.
iii) O método dos volumes finitos constitui uma aproximação alternativa à diferença
finita e utiliza como ponto de partida a seguinte equação de transporte genérica [26,27]:
( )( ) ( )div u div grad q
t θ
ρθρθ θ
∂+ = Γ
∂r
+
r
(2.7)
onde o primeiro termo representa a variação temporal da variável dependente θ, osegundo termo, é o termo convectivo onde é o vector velocidade. No segundo
membro está presente o termo de difusão (Γ é o coeficiente de difusão apropriado) e o
termo fonte qθ.
u
A integração da equação (2.7) sobre um volume de controlo (VC) constitui o passo
fundamental do método dos volumes finitos e distingue-o de todos os outros métodos.
26
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 43/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
( )( ) ( )
VC VC VC VC
dV div u dV div grad dV q dV t
θ
ρθρθ θ
∂+ = Γ +
∂∫ ∫ ∫ ∫ r (2.8)
Os integrais de volume, relativos aos termos convectivo e difusivo, podem ser reescritoscomo integrais de superfície do VC utilizado o teorema da divergência de Gauss.
Para um vector genérico a , o teorema da divergência estabelece que [28]:r
.VC SC
div a dV a nds=∫ ∫ r r r (2.9)
cuja interpretação física permite afirmar que é a componente do vector a na
direcção normal ao elemento de superfície ds. Este teorema estabelece que o integral
do divergente de um vector a sobre um volume é igual à componente de sobre uma
superfície que lhe é normal.
.a nr r
.r
nr
.r
.ar
Da aplicação do teorema da divergência à expressão 2.8 resulta que:
.( ) .( )SC SC VC
VC
dV n u dS n grad dS q dV t
θρθ ρθ θ ∂
+ = Γ + ∂
∫ ∫ ∫ ∫ r r r (2.10)
onde a ordem de integração e de diferenciação foi alterada no primeiro termo do
primeiro membro da equação (2.10) de modo a ilustrar melhor o seu significado físico.
Este termo representa a variação temporal da acumulação total da propriedade θ no VC.
Em problemas transitórios é necessário integrar em relação ao tempo t sobre um pequeno
intervalo ∆t , isto é, de t a t+∆t . Tal integração permite obter a forma mais geral da
equação de transporte da propriedade θ.
.( ) .( )t VC t SC t SC t VC
dV dt n u dS dt n grad dS dt q dV dt t
θρθ ρθ θ∆ ∆ ∆ ∆
∂+ = Γ +
∂ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫
r r r (2.11)
O domínio de cálculo é subdividido num número finito de volumes de controlo (VC)
contíguos que, em contraste com o método das diferenças finitas, definem os limites dos
volumes de controlo e não os nós computacionais. A equação (2.11) é, então, aplicada a
cada VC e as variáveis dependentes podem ser avaliadas no centro dos volumes ou nos
27
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 44/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
vértices conforme se trate da formulação em célula centrada (cell center ) ou da
formulação utilizando os vértices da célula computacional (cell vertex), respectivamente.
A formulação em célula centrada é habitualmente mais usada pelo que será utilizada,também, neste trabalho. O centro geométrico de cada VC constitui, então, um nó
computacional onde a variável dependente será calculada. Para expressar os valores da
variável dependente nas faces dos VC em termos dos valores nos centros dos VC
utilizam-se as técnicas de interpolação. Os integrais de volume e de superfície são
aproximados usando fórmulas de quadratura adequadas. Como resultado, temos uma
equação algébrica para cada VC onde aparece um determinado número de valores
nodais, dependendo do número de dimensões consideradas, referentes a nós vizinhos.
Este método aceita qualquer tipo de malha, e é adequado para geometrias complexas. A
malha define apenas as fronteiras dos VC e necessita de ser relacionada com um sistema
de coordenadas. O método é conservativo por natureza, e portanto os integrais de
superfície são idênticos para os VC que partilham a mesma fronteira.
A aproximação pelo método dos volumes finitos é talvez a mais simples de entender
devido a que todos os termos que necessitam ser aproximados têm significado físico,motivo pelo qual é muito popular entre os Engenheiros.
A desvantagem deste método em relação ao método das diferenças finitas é o facto de os
esquemas de ordem superior a dois serem mais difíceis de desenvolver. Isto deve-se ao
facto da aproximação por volumes finitos requerer dois níveis de aproximação:
interpolação e integração.
A aproximação habitual passa por definir os volumes de controlo através de uma malha
adequada e colocar o nó computacional no centro do volume de controlo. Todavia,
poderemos – para malhas regulares – definir primeiro as posições nodais, construindo
posteriormente os volumes de controlo em torno destes, para que as suas faces fiquem
colocadas no centro de cada par de nós adjacentes.
A vantagem da primeira aproximação é que o valor nodal representa a média de todo o
volume de controlo o que permite maior precisão (segunda ordem) do que a segunda
28
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 45/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
aproximação, desde que se garanta que os nós computacionais coincidam com os centros
de massa dos volumes de controlo.
A vantagem da segunda aproximação é sentida quando se utiliza o esquema centrado para aproximar as derivadas nas faces dos volumes de controlo, sendo mais precisa
quando elas estão centradas entre nós vizinhos.
Os princípios de discretização são os mesmos para todas as variantes do método,
variando apenas as localizações dos nós no volume de integração. Para obter uma
equação para cada VC, é necessário aproximar os integrais de superfície e de volume por
fórmulas de quadratura adequadas. É o que se descreve nas próximas subsecções.
II.3.1. Aproximação de integrais de superfície
Apresenta-se na figura 2.6 um volume de controlo 2D típico, assim como as respectivas
notações a utilizar ao longo deste trabalho.
A superfície do volume de controlo pode ser subdividida em quatro (2D) ou seis (3D)faces planas, denotadas no caso 2D por letras correspondentes à sua direcção (e, w, n, s)
em relação ao nó central P .
O fluxo líquido através da fronteira do volume de controlo (VC) é a soma dos integrais
sobre as seis faces [27]:
∑∫ ∫ =k
S S ds f ds f
k
(2.12)
onde f é a componente vectorial da difusão na direcção normal à face do VC. Para
manter a conservação é importante que os VC não se sobreponham, cada face do VC é
única para os dois VC que contactam entre si.
29
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 46/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
30
η j+1
NW NE N •
•
•
η j
nw n ne
w e
sw s se
W E P •
•
η j-1
•••
S SW SE
ξi-1 ξi ξi+1
Figura 2.6 – Volume de controlo 2D típico e notação usada. As faces do volume de controlo sãodenotadas por letras minúsculas e os nós centrais por letras maiúsculas.
Para calcular o integral de superfície presente na expressão (2.12) é necessário conhecer
o integrando f em toda a superfície S e. Como não dispomos dessa informação e como só
os valores nodais (centros do VC) são calculados teremos de introduzir uma
aproximação. A melhor forma de o fazer é usando dois níveis de aproximação [27]:
o
O integral é aproximado em função dos valores das variáveis num ou mais locaisda face considerada;
o Os valores na face são aproximados em função dos valores nodais (centros dos
VC).
A aproximação mais simples do integral é a que nos é proporcionada pela regra do ponto
médio: o integral é aproximado como o produto do integrando no centro da face da
célula considerada (que é também uma aproximação dada pelo valor médio de toda a
superfície) e a área dessa mesma face (será considerada apenas a face Este nesta
demonstração, sendo as restantes tratadas de modo similar) [27,29]:
ee e e
S F f ds f S f = = ≈∫ e eS (2.13)
esta aproximação do integral – permitida pelo facto do valor de f na face ‘e’ ser
conhecido é de segunda ordem. Quando os valores de f na face da célula não estiverem
disponíveis, terão de ser obtidos por interpolação. De modo a que a precisão de segunda
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 47/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
ordem, que a regra do ponto médio permite, seja mantida, o valor de f e terá que ser
obtido pelo menos com a mesma precisão.
Uma outra aproximação de segunda ordem para o integral de superfície em 2D é a regrado trapézio, que se apresenta de seguida [27,29]:
)(2 sene
e
S e f f
S ds f F
e
−≈= ∫ (2.14)
neste caso necessitamos do valor do integrando nos vértices do volume de controlo.
Para aproximações de ordem superior a dois dos integrais de superfície, o valor dointegrando é necessário em mais do que dois locais. Por exemplo, a regra de Simpson é
uma aproximação de quarta ordem, e estima o valor do integral de superfície da seguinte
forma [27,29]:
)4(6 seene
e
S e f f f
S ds f F
e
++≈= ∫ (2.15)
onde são necessários os valores de f em três locais distintos. No centro da face ‘e’ e nosvértices, ‘ne’ e ‘ se’. De forma a manter a precisão desta aproximação, estes valores
deverão ser obtidos por interpolação dos valores nodais com a mesma ordem de precisão.
Os polinómios de terceiro grau são uma possibilidade.
II.3.2. Aproximação de integrais de volume
O termo fonte (q) que representa a geração interna de calor num determinado fenómeno
físico requer uma integração sobre o volume do VC. A aproximação mais simples é a de
segunda ordem, onde o integral de volume é substituído pelo produto do valor médio do
termo fonte pelo volume do VC, como se mostra na expressão (2.16) [27,29].
V qV qdV qQV
P P ∆≈∆== ∫ (2.16)
31
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 48/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
Onde q P representa o valor de q no centro do VC. Esta quantidade é facilmente calculada
se todas as variáveis no nó central P forem conhecidas, e consequentemente não existe
necessidade de interpolar. A aproximação descrita pela expressão (2.16) representa
fielmente o valor de q se este for constante ou variar linearmente ao longo do VC. Uma
aproximação de ordem superior requer valores de q noutras posições além dos centros do
VC. Estes valores terão que ser obtidos por interpolação de valores nodais conhecidos,
ou pelo uso de funções de forma (técnica habitual do método dos elementos finitos).
II.3.3. Técnicas de interpolação
As aproximações dos integrais requerem o valor das variáveis noutros locais além dos
nós computacionais (centro do VC). O integrando, f , envolve o produto de muitas
variáveis e/ou gradientes das mesmas nesses locais. Para o cálculo do fluxo por difusão
necessitamos conhecer o valor de θ e do seu gradiente na direcção normal à face da
célula considerada numa ou mais posições da superfície do VC. Os integrais de volume
relativos ao termo fonte podem também requerer esses valores, estes terão que ser
expressos em termos dos valores nodais por interpolação. Existe um grande número de possibilidades de aproximação entre as quais se refere apenas as que são comummente
utilizadas.
i) A formulação UDS (Upwind Differencing Scheme) aproxima o valor de θe pelo valor
imediatamente acima da face Este, o que é equivalente a usar a formulação progressiva
ou regressiva (dependendo da direcção do escoamento) da diferença finita na
aproximação da primeira derivada [3,27]:
se ( . ) 0
se ( . ) 0 P e
e
E e
v n
v n
θθ
θ
>=
<
r r
r r (2.17)
Esta é a única aproximação que satisfaz as condições de fronteira incondicionalmente,
isto é, nunca irá ter uma solução oscilatória.
32
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 49/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
ii) A interpolação linear (CDS - Central Difference Scheme) entre dois nós adjacentes
constitui uma outra possibilidade de aproximação do valor no centro da face da variável
dependente. Na face Este de um VC presente numa malha cartesiana temos que [3,27]:
(1 )e E e P eθ θ λ θ λ= + − (2.18)
onde o factor λe se define como:
e P e
E P
x x
x xλ
−=
− (2.19)
onde xe , x E e x P são as abcissas: da face Este do volume de controlo; do nó central do VC
a Este do VC considerado e do nó central do VC considerado, respectivamente.
A expressão (2.19) possui precisão de segunda ordem, facto que pode ser provado pelo
desenvolvimento em série de Taylor de θ E em torno do ponto x P . O resultado é o que se
apresenta de seguida [27,29]:
2
2
( )( )(1 ) 2
e P E ee E e P e
P
x x x x
O x
θθ θ λ θ λ
− − ∂= + − − + ∂ (2.20)
onde O representa os termos de ordem superior a dois. O erro de truncatura é
proporcional à raiz quadrada do espaçamento da malha utilizada quer ela seja uniforme
ou não.
Esta formulação é equivalente à formulação centrada do método das diferenças finitas,
daí a denominação CDS. Como todas as aproximações realizadas são de ordem superior a um, esta formulação pode produzir soluções oscilatórias.
A assunção de uma relação linear entre os nós P e E oferece também a aproximação mais
simples de um gradiente, importante para a avaliação dos fluxos de difusão:
E P
e E P x x
θ θθ −∂ ≈ ∂ − x (2.21)
33
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 50/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
utilizando a expansão em série de Taylor de θe pode-se demonstrar qual o erro de
truncatura que deriva da utilização da expressão (2.21) [27,29]:
2 2 3 32 3
2 3( ) ( ) ( ) ( )
2( ) 6( )e P E e e P E e
T
E P E P e e
x x x x x x x xO
x x x x x x
θ θε − − − − − −∂ ∂= − − −∂ ∂ + (2.22)
quando a localização de ‘e’ está ao centro entre os nós P e E (caso se tenha uma malha
uniforme), esta aproximação é de segunda ordem. Se o primeiro termo do segundo
membro da expressão (2.22) for nulo, a expressão resultante é então proporcional a ∆ x2.
Quando a malha é não ortogonal, o erro resultante é proporcional ao produto entre ∆ x e
ao factor de expansão da malha.
Interpolações de ordem superior a três só fazem sentido se os integrais forem
aproximados usando também fórmulas de ordem equivalente. A utilização da regra de
Simpson para aproximar os integrais implica que a interpolação seja realizada com
polinómios de pelo menos grau três, obtendo-se erros de quarta ordem.
II.4. Discretização em problemas transitórios
Quando se está perante um problema transitório uma nova variável independente há a
considerar – o tempo.
Tal como as coordenadas espaciais, a coordenada tempo também precisa ser
discretizada.
Para problemas de valor inicial é suficiente considerar uma equação diferencial de
primeira ordem genérica e uma condição inicial também genérica [27,29].
00
( )( , ( )) ; ( )
t f t t t
t
θθ θ θ
∂= =
∂ (2.23)
Pretende-se determinar os valores de θ passado um pequeno intervalo de tempo ∆t
depois da condição inicial. A solução em t 1 = t 0 + ∆t , θ1
, pode ser vista como uma nova
34
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 51/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
condição inicial, podendo-se assim avançar para a solução em t 2 = t 1 + ∆t , e assim
sucessivamente.
O método mais simples é o obtido por integração da expressão (2.23) de t n a tn+1= t n +
∆t
1 1
1 ( , ( ))n n
n n
t t
n n
t t
d dt f t t dt
dt
θθ θ θ
+ +
+= − =∫ ∫ (2.24)
originando uma solução exacta. No entanto, como o termo mais à direita não pode ser
calculado sem o conhecimento da solução, é, portanto, necessário introduzir uma
aproximação. O teorema do valor médio permite o cálculo do integral num ponto
concreto t=τ entre t n e t n+1, o integral é igual a f( τ ,θ(τ)) ∆t , mas este resultado não tem
qualquer interesse porque τ também não é conhecido. Existem, no entanto, outras
formas de aproximar este integral. É o que se aborda de seguida utilizando quatro
simples regras de quadratura numérica.
Se o integral for estimado utilizando o valor do integrando no ponto inicial (regra do
rectângulo à direita), temos que [27:29]:
1 ( , )n n n
n f t t θ θ θ+ = + ∆ (2.25)
que é conhecida como a formulação explícita de Euler.
Se em vez do ponto inicial se utilizar o ponto final (regra do rectângulo à esquerda),
obtém-se a expressão que dá forma à formulação implícita de Euler [27:29]:
1 11( , )n n nn f t t θ θ θ+ ++= + ∆ (2.26)
Utilizando o ponto médio do intervalo obtém-se [27]:
12
12
1 ( , )nn n
n f t t θ θ θ
+++
= + ∆ (2.27)
35
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 52/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
expressão que é conhecida como a regra do ponto médio e pode ser entendida como a
base de um importante método de resolução de equações diferenciais parciais – a
formulação de Leapfrog.
Para finalizar pode-se ainda usar uma linha recta entre os pontos inicial e final para
construir a aproximação:
11
1( , ) ( , )
2n n n n
n n f t f t t θ θ θ θ++ = + + ∆
1+ (2.28)
esta é a regra do trapézio, e constitui a base de um importante método de resolução de
equações diferenciais parciais – a formulação de Crank-Nicolson [27].
Geralmente este conjunto de quatro métodos é designado por métodos de duplo nível,
porque envolvem os valores das incógnitas em apenas dois níveis de tempo. O primeiro
método pertence à classe denominada de explícitos enquanto que os restantes são
denominados de implícitos. Notar que todos os métodos, exceptuando o método
explícito, requerem o valor θ(t) noutro ponto além do inicial (t=t n). De modo que o
segundo termo da expressão (2.24) só pode ser avaliado de forma aproximada.
Poderemos ainda aproximar os integrais da expressão (2.24) por valores médios do
intervalo de integração. Richtmayer e Morton (1967) apresentaram um algoritmo
genérico utilizando uma formulação implícita com três níveis de tempo [30].
( )1 1
1
12
n n n n
n
d
dt t t
θ θ θ θβ β
+ +
+
− − ≈ + − ∆ ∆
1θ −
(2.29)
onde β é o factor de relaxação (1≥β≥0).
A expressão (2.29) pode ser reescrita da seguinte forma:
1 1 11(1 ) ( , )n n n n
n f t t β θ θ βθ θ+ − +++ = − + ∆ (2.30)
transformando-se na expressão que constitui a formulação implícita de Richtmayer e
Morton.
36
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 53/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
II.5. Aproximação numérica de derivadas
O processo de resolução de uma equação diferencial por intermédio de métodos
numéricos implica a sua substituição por equações algébricas. No método dos volumesfinitos a aproximação mais usual é a chamada cell center . Esta permite obter para cada
ponto aproximações convenientes das derivadas, através das diferenças entre os valores
que a função toma nos centros dos volumes de controlo.
A aproximação dos integrais de volume e de superfície será efectuada utilizando a regra
do ponto médio e apresentada apenas para a face Este do volume de controlo. Para as
restantes faces é utilizado o mesmo conceito de aproximação, cujo cálculo poderá ser
feito alterando os índices na expressão que se apresenta. Esta análise não é aplicável
apenas a volumes de controlo de quatro faces, mas sim a um volume de controlo com
qualquer número de faces.
A regra do ponto médio aplicada à forma integral do termo da difusão resulta na
seguinte igualdade [27]:
( )ˆ ˆ. .e
d
e eeS F k grad n ds k grad n S θ θ= ≈∫ (2.31)
Existem diversas formas de aproximar a derivada de θ normal à face da célula. São
descritas apenas duas delas que se consideram mais importantes.
Se a variação de θ na vizinhança da face da célula for descrita por uma função de
forma, é possível diferenciar esta função na posição ‘e’ de modo a encontrar as
derivadas com respeito às coordenadas cartesianas. O fluxo de difusão é então:
d i
e e
i i e
F k x
θ ∂=
∂ ∑ eS (2.32)
Uma outra forma de calcular as derivadas na face da célula passa por obter θ primeiro
no centro do VC, interpolando-a em seguida para obter o valor de θe nas faces.
Tal poderá ser feito de um modo simples se for utilizado o teorema de Gauss,
aproximando a derivada no centro do VC através do valor médio da célula [27].
37
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 54/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
V i
i P
dV x
x V
θ
θ
∂ ∂∂
≈ ∂ ∆
∫ (2.33)
Poder-se-á agora considerar a derivada i xθ∂ ∂ como o divergente do vector θi,
transformando o integral de volume num integral de superfície.
ˆ ˆ. , , , , ,...i
c cV S
ci
dV i n dS S c e n w s x
θθ θ
∂= ≈ =
∂ ∑∫ ∫ (2.34)
As expressões (2.33) e (2.34) dizem que é possível calcular o gradiente de θ em relação
a x no centro do VC somando os produtos de θ pelas componentes na direcção x dos
vectores de superfície de todas as faces do VC e dividindo pelo volume do VC:
i
c c
c
i P
S
x V
θθ ∂
≈ ∂ ∆
∑ (2.35)
para θc, no caso de malhas cartesianas, pode-se utilizar a interpolação linear, é então
obtida a aproximação por diferença finita centrada [27]:
2 E W
i P x x
θ θθ −∂≈ ∂ ∆
(2.36)
as derivadas calculadas pela expressão (2.36) podem ser interpoladas para a face da
célula e o fluxo de difusão pode ser calculado pela expressão (2.34). O problema
inerente a esta aproximação é que poderá gerar uma solução oscilatória no decurso do
processo de iteração.
II.6. Resolução de sistemas de equações lineares
Como vimos do processo de discretização, qualquer que seja o método utilizado, resulta
num sistema de equações algébricas que no caso da equação da energia é linear e pode
ser genericamente representado por [27,29,30]:
38
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 55/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
A Qθ = (2.37)
em que A é uma matriz de coeficientes esparsa (isto é, possuem muitos elementos de
valor igual a zero), θ é um vector contendo as variáveis nos nós computacionais e Q é
um vector de termos independentes.
A equação algébrica para um volume de controlo particular é da forma [26,27]:
P P nb nb P
nb
A Aθ θ+ =∑ Q (2.38)
onde P representa o nó onde a equação de derivadas parciais é aproximada e o índice nb
representa os nós vizinhos envolvidos na aproximação. Q P é o termo fonte, ou seja, não
contém incógnitas.
Como o problema que se pretende resolver implica a resolução de um sistema de
equações linear é pertinente analisarmos aqui alguns dos métodos existentes para esse
fim.
II.6.1. Eliminação de Gauss
O método de eliminação de Gauss constitui o método básico e directo para a resolução
de sistemas lineares de equações algébricas, e baseia-se na redução sistemática de um
grande sistema de equações para um mais pequeno. Ao longo deste procedimento, os
elementos da matriz A são modificados, mas as variáveis dependentes não sofrem
quaisquer alterações, o que torna possível e conveniente efectuar a descrição do métodoem função apenas da matriz A [27,29]:
=
NN
N
N
N N N A
A
A
A A A
A A A
A A A
AM
L
OMMM
L
L
2
1
321
232221
131211
(2.39)
39
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 56/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
O ponto fulcral deste algoritmo é a técnica para eliminar todos os elementos da matriz A
que se encontram abaixo do elemento A11, isto é, substitui-los por zero. Para o efeito
basta multiplicar a cada linha da matriz (cada linha representa uma equação) por
111 A Ai
e subtrai-la pela linha i; neste processo todos os elementos das linhas 2 a n da
matriz são modificados assim como os elementos 2 a n do vector de termos
independentes presente no segundo membro da expressão (2.37). Quando este processo
tiver terminado nenhuma das equações 2 a n conterá a variável θ1. Nesse momento
teremos um sistema com n-1 equações para as variáveis θ2 a θn. O mesmo procedimento
é aplicado ao sistema de equações já reduzido, sendo agora objectivo eliminar todos os
elementos abaixo do elemento A22 e assim sucessivamente.
Este procedimento aplicado para as colunas 1 a n-1 é designado por eliminação
progressiva. Depois de se completar todo o processo a matriz A original terá sido
substituída por uma outra triangular superior, isto é, todos os elementos abaixo da
diagonal principal são nulos.
Todos os elementos excepto os da primeira linha da matriz A são diferentes dos
originais, assim como os elementos do vector de termos independentes Q.
Os sistemas triangulares são de fácil resolução porque se reparamos a última linha ou
equação contém apenas uma incógnita pelo que pode ser resolvida directamente:
nn
nn
Q
Aθ = (2.40)
a equação imediatamente acima contém apenas as incógnitas θn-1 e θn, e como esta
última é já conhecida pela expressão (2.40), podemos resolver também esta equação e
determinar o valor de θn-1 e assim sucessivamente. Generalizando este procedimento
cada equação pode ser resolvida pela seguinte expressão designada por eliminação
regressiva.
1
N
i i
j i
i
ii
Q A
A
θ
θ = +
−
=∑ j j
(2.41)
40
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 57/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
O termo da direita da expressão anterior é calculável porque como os valores de θ j são já
conhecidos.
II.6.2. Factorização LU
Tal como se teve oportunidade de constatar no método de eliminação de Gauss a matriz
original é transformada numa outra triangular superior. O processo responsável por essa
transformação pode ser realizado de um modo mais formal, multiplicando a matriz
original por uma matriz triangular inferior. Por si só o interesse deste procedimento
parece reduzido, mas tendo em consideração que a inversa de uma matriz triangular
inferior é também triangular inferior, então qualquer matriz A, sujeita a algumas
limitações que podem ser aqui ignoradas, pode ser escrita como o produto das matrizes
triangulares: inferior ( L) e superior (U ) [27,29]. Este processo é denominado por
factorização.
O que torna a factorização útil e atractiva é o facto da sua implementação ser
relativamente simples. A matriz triangular superior U é precisamente a que se obtém na primeira fase do método de eliminação de Gauss (eliminação progressiva). Mais ainda,
os elementos de L são os factores multiplicativos usados no processo de eliminação
( Aij /Aii). Isto permite que a factorização seja realizada com uma mínima alteração do
método de eliminação de Gauss.
A utilização desta factorização permite resolver um sistema de equações em duas etapas.
A primeira consiste na definição de [27]:
U θ =Y, (2.42)
a segunda é uma consequência da primeira, e o sistema de equações (2.37) transforma-se
em [27]:
LY = Q. (2.43)
41
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 58/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
O sistema de equações representado pela expressão (2.43) pode ser resolvido por uma
variante do método utilizado na segunda fase da eliminação de Gauss (eliminação
regressiva) na qual se inicia a substituição pelo topo da matriz em detrimento da última
linha. Uma vez resolvida a equação (2.43), obtém-se Y , tornando possível calcular θ
através da expressão (2.42).
Existem algumas variantes deste método que constituem a base de alguns dos melhores
métodos iterativos para a resolução de sistemas de equações lineares, como é o caso do
método de Stone, que é a principal razão por que foi aqui abordado.
II.6.3. Métodos iterativos
Qualquer sistema de equações linear pode ser resolvido por eliminação de Gauss ou pela
factorização LU. Os erros de discretização são normalmente superiores à precisão
aritmética do computador pelo que não existe motivo para resolver o sistema de
equações com essa precisão. Este facto abre uma porta aos métodos iterativos [27].
Nos métodos iterativos utiliza-se uma determinada expressão para melhorar
constantemente a solução. Se as iterações forem de exigência moderada e baixas em
número, o método iterativo pode ser menos dispendioso que o método directo.
Aplicando n iterações à expressão (2.37) obtém-se uma solução aproximada θ n que pode
ser escrita matematicamente da seguinte forma [21,27]:
n n
A Q Rθ = − (2.44)
em que Rn é o valor residual. Por subtracção da expressão (2.44) pela expressão (2.37),
obtemos uma relação entre o erro de convergência definido por:
n nε θ θ= − (2.45)
onde θ é a solução exacta, e o valor residual é:
42
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 59/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
n n A Rε = (2.46)
O objectivo do processo iterativo é tornar o valor residual nulo, o que faria da solução
encontrada a exacta. Para o efeito, considere-se o seguinte esquema iterativo aplicado a
um sistema de equações linear [27]:
1n nM N θ θ+ = + B
θ
n R
n
. (2.47)
Onde M e N são matrizes de coeficientes e B é o vector de termos independentes.
Uma propriedade que deverá obviamente ser respeitada por um método iterativo é que a
solução final satisfaça a expressão (2.37) Como, por definição, em convergência,
, deveremos ter:1n nθ θ+ = =
e A M N B Q= − = (2.48)
ou, genericamente,
eCA M N B CQ= − = (2.49)
em que C é a chamada matriz de pré-condicionamento.
Uma versão alternativa deste método pode ser obtida por subtracção de a ambos os
termos da expressão (2.47), obtendo-se:
nM θ
( ) ( )1 oun n n nM B M N M θ θ θ δ+ − = − − = (2.50)
onde é denominado correcção ou actualização e constitui uma
aproximação ao erro de convergência [27].
1n nδ θ θ+= −
Quando se pretende utilizar um método iterativo o critério de selecção é,
essencialmente, a rapidez de convergência. Justifica-se, portanto, prestar alguma
atenção ao modo como determinar essa característica.
43
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 60/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
O primeiro passo é a derivação da equação que determina o erro de convergência. Para a
encontrar, relembra-se que, em convergência, θn+1=θn
=θ, e deste modo a solução
convergente obedece à expressão:
M N θ θ= + B (2.51)
Subtraindo esta equação pela expressão (2.47) e fazendo uso da definição (2.45), vem
que:
1 1n nM N ε ε+ −= (2.52)
O método iterativo converge se o limite de εn, quando n tende para infinito é igual azero. O aspecto crítico do método é promovido pelo valor próprio νk e vector próprio ψκ
da matriz definida por:1M N −
1 k k
k M N ψ ν ψ− = (2.53)
onde o índice k representa o número de equações. Assume-se que os vectores próprios
constituem uma base para R
n
, isto é, um espaço vectorial que comporta todas as n componentes vectoriais. Pelo que o erro inicial pode ser expresso por [27]:
0
1
K k
k
k
aε ψ=
= ∑ (2.54)
onde ak é uma constante. Então o procedimento iterativo (2.52) é realizado do seguinte
modo:
∑∑==
−− === K
k
k
k k
K
k
k
k aa N M N M 11
1011 ψνψεε (2.55)
e, por indução, vem que:
∑=
= K
k
k n
k k
n a1
)( ψνε (2.56)
44
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 61/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
claro que, se εn se tornar nulo para valores elevados de n, a condição necessária e
suficiente é que todos os valores próprios deverão ter magnitude menor que a unidade.
Ao vector próprio de maior amplitude chama-se raio espectral da matriz M -1 N . De facto,
depois de algumas iterações, os termos da expressão (2.56) que contêm valores próprios
de pequena magnitude tornam-se ainda mais pequenos e apenas o termo correspondente
ao valor próprio de maior magnitude (que podemos considerar como sendo ν1 e assumi-
lo como único) permanece, pelo que a expressão (2.56) pode ser escrita da seguinte
forma [27]:
11 1( )n n
aε ν≈ ψ (2.57)
se a convergência for definida como a redução do erro de convergência com tolerância
δ, vem que:
1 1( )naδ ν≈ (2.58)
aplicando logaritmos neperianos a ambos os termos da expressão anterior, obtém-se
uma expressão para o número de iterações necessárias:
1
1
ln
ln
ν
δ
≈a
n (2.59)
da qual se conclui que se o raio espectral adquirir valores muito próximos da unidade, o
processo iterativo irá convergir muito lentamente [27].
De referir que o cálculo dos valores próprios não é fácil (a maior parte das vezes não é
conhecido explicitamente) pelo que será necessário usar aproximações, é também
necessário estimar o erro de convergência de forma a decidir quando o processo
iterativo termina.
45
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 62/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
II.6.3.1. Método de Jacobi
No método de Jacobi as equações resultantes do processo de discretização são resolvidas
isoladamente, pelo que a expressão 2.38 é alterada assumindo o seguinte aspecto [21,29]:
1
N
i ij
j
i
ii
Q A
A
θ
θ=
−
=∑ j
(2.60)
A expressão (2.60) sugere um método iterativo definido por,
1
1
N k
i ij j
jk
i
ii
Q A
A
θ
θ
−
=
− =
∑
(2.61)
onde os termos de θ j presentes no segundo membro da equação (2.61) são todos relativos
à última iteração efectuada.
II.6.3.2. Método de Gauss-Seidel
Uma explicação que intuitivamente motiva o método de Gauss-Seidel é a seguinte. No
método de Jacobi utilizam-se os valores θ j da iteração anterior para obter os valores de
θ P da iteração seguinte. No entanto, aquando do cálculo de θi já são conhecidos os
valores actuais de θ j. O método de Gauss-Seidel, em contraste com o método de Jacobi,
utiliza os valores actuais de θ j em detrimento dos da iteração anterior [21,29].
A concretização desta ideia conduz à seguinte modificação da expressão (2.61).
11
1
i N k k
ij j ij j
j j ik
i
ii
Q A A
A
θ θ
θ
−−
= =
− −
=∑ ∑
(2.62)
46
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 63/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
Habitualmente, este método possui maior rapidez de convergência que o método de
Jacobi.
II.6.3.3. Método de Stone
A factorização LU é um excelente método para resolver sistemas de equações lineares,
mas não tira qualquer partido do facto de um determinado sistema possuir uma matriz de
coeficientes dispersos. Num método iterativo, se M (matriz iterativa) for uma boa
aproximação a A, a convergência é rápida. Estas duas afirmações são a ideia base da
utilização de uma aproximação de A por factorização LU com a matriz de iteração M :
M LU A N = = + (2.63)
onde L e U são ambas esparsas e N é pequena.
A versão assimétrica deste método é denominada por factorização LU incompleta ou
ILU. O procedimento é idêntico à factorização LU, mas para os muitos elementos nulos
presentes na matriz A o correspondente elemento em L ou U é igualado a zero. Estafactorização não é exacta, mas o produto destes factores pode ser usado como matriz M
do método iterativo. Este método é de convergência lenta.
Um outro método de factorização foi o proposto por Stone (1968). O método de Stone
também designado por SIP (Strongly Implicit Procedure) tem especial orientação para
equações algébricas resultantes da discretização de equações de diferenciais parciais
[27].
Para a descrição deste método será considerada uma molécula computacional de nove
pontos correspondente à utilização de coordenadas curvilíneas em 2D.
47
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 64/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
USE
USW
UE
LSW
L NW
LW
U N
1 LSLP =
MSE
MSSE
MS
ME
M N
M NNW
MW
MSW
MP
MSS
MS
M NN
M N
Figura 2.7 – representação esquemática das matrizes L e U bem como da matriz produto M; asdiagonais extra estão representadas a traço interrompido.
Tal como na factorização ILU, as matrizes L e U têm elementos não nulos apenas nas
diagonais onde A tem elementos não nulos. O produto das matrizes triangulares inferior esuperior com estas estruturas têm mais diagonais não nulas que A. Para uma molécula
computacional de nove pontos existem mais quatro diagonais correspondentes aos nós
NNW , SSE , NN e SS , tal como se mostra na figura 2.7.
Para tornar estas matrizes unitárias, todos os elementos da diagonal principal de U são
colocados iguais à unidade. Estas nove introduções (cinco em L, quatro em U )
necessitam ser determinadas. Para matrizes da forma como se mostra na figura 2.7 as
regras de multiplicação de matrizes dão-nos os elementos do produto de L por U ,
M = LU seguintes:
SW SW
W SW N W
NW W N NW
NNW NW N
SS SW SE
S SW E W SE S
P SW NE W E NW SE S N
N W NE NW E P N
NN NW NE
SSE S SE
SE S E P E
E S NE P E
NE P NE
M L
M L U L
M L U L
M L U
M L U M L U L U L
M L U L U L U L U L
M L U L U L U
M L U
M L U
M L U L U
M L U L U
M L U
=
= +
= +
=
== + +
= + + + +
= + +
=
=
= +
= +
=
P M SW = LSW (2.64)
48
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 65/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
Pretendemos seleccionar L e U , tal que M seja a melhor aproximação possível a A. No
mínimo, N deve conter as quatro diagonais de M que correspondem a diagonais nulas de
A. Uma escolha óbvia é que N tenha elementos não nulos apenas nessas quatro
diagonais, e obrigar as outras diagonais de M a serem iguais às correspondentes
diagonais em A. Este raciocínio é perfeitamente correcto pois trata-se do processo
standard do método ILU já mencionado. Infelizmente, este método converge lentamente.
Stone (1968) descobriu que a convergência podia ser incrementada se fosse permitido
que N tivesse elementos não nulos nas diagonais correspondentes a todas as diagonais
não nulas de LU [27]. O método é demonstrado de um modo simples se considerarmos o
vector M θ :
( ) P P S S N N E E W W NE NE NW NW P
SE SE SW SW NNW NNW SSE SSE NN NN SS SS
M M M M M M M M
M M M M M M
θ θ θ θ θ θ θ θ
θ θ θ θ θ θ
= + + + + + +
+ + + + +
+
≈
* ) 0≈
*
(2.65)
os quatro últimos termos são os extraordinários. Cada termo da expressão (2.65)
corresponde a uma diagonal de M = LU .
A matriz N contém as quatro diagonais extra de M , e queremos escolher os elementosdas restantes diagonais tal que ou seja,0 N θ ≈
0 P P S S N N E E W W NE NE NW NW
SE SE SW SW NNW NNW SSE SSE NN NN SS SS
N N N N N N N
N N N N N N
θ θ θ θ θ θ θ
θ θ θ θ θ θ
+ + + + + + +
+ + + + + (2.66)
o que requer que a contribuição dos quatro termos extra da expressão anterior seja
anulada pela contribuição das outras diagonais, isto é, a expressão (2.66) será reduzida à
seguinte expressão:
* * *( ) ( ) ( ) ( NNW NNW NNW SSE SSE SSE NN NN NN SS SS SS M M M M θ θ θ θ θ θ θ θ− + − + − + − (2.67)
onde θ θ são aproximações de θ θ , respectivamente.* * *, , e NNW SSE NN SS θ θ , , e NNW SSE NN SS θ θ
A ideia chave do método é a seguinte: uma vez que estas expressões aproximam
equações elípticas de derivadas parciais, a solução mais provável é ser uma curva suave.
49
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 66/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
Consequentemente, θ θ podem ser aproximados em termos dos valores
de θ nos nós correspondentes às diagonais de A. A proposta feita neste trabalho para
essa aproximação é a seguinte:
* * *, , e NNW SSE NN SS θ θ *
α
( )
( )
( )
( )
*
*
*
*
NNW NW N W
NN N NE P
SS S SW P
SSE S SE E
θ α θ θ θ
θ α θ θ θ
θ α θ θ θ
θ α θ θ θ
= + −
= + −
= + −
= + −
(2.68)
se =1, tratam-se de interpolações cuja precisão é de segunda ordem, mas Stone
descobriu que para que o método seja estável α<1.
Se esta aproximação for substituída na expressão (2.67), obtêm-se todos os elementos de
N como combinações de M NNW , M SSE , M NN e M SS . Os elementos de M em (2.68), podem
agora ser igualados à soma dos elementos de A e N .
As equações resultantes não são por si só suficientes para determinar todos os elementos
de L e U , mas podem ser resolvidas de um modo sequencial iniciando no canto Sudoeste
da malha.
Os coeficientes devem ser calculados pela ordem apresentada. Para os nós junto às
fronteiras, os respectivos elementos das matrizes são entendidos como sendo nulos.
50
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 67/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
1, 1
, 1, 1 1, 1
1,
1, 1
1, 1 1, 1 1,
, 1
1, 1
1
1
1
(
ij
ij SW
SW i j
SE
ij ij i j i j ij i j
W W NW N SW N
ij ij i j
ij NW W N
NW i j
N
ij ij i j ij i j ij i j
ij S SW SE SW E W SE
S i j
SE
ij ij ij i j ij
P P SW SE NW N
A L
U
L A L U L U
A L U L
U
A L U L U L U L
U
L A L U L U
α
α
α
α
α
α
− −
− + − −
−
− +
− − − − −
−
− −
=+
= + −
−=
+
− − −=
+
= − + 1, 1 1, 1 1, 1, 1 , 1
1, 1 1, 1, 1
,
, 1
,
, 1 ,
)i j ij i j ij i j ij i j ij i j
E SW NE W E NW SE S N
ij ij i j ij i j ij i jij N NW NE W NE NW E
N ij i j
P NW
ij ij i j
ij SE S E
SE ij i j
P S
ij ij i j ij i j
ij E S SE S NE
E
L U L U L U L U
A L U L U L U U
L L
A L U U
L L
A L U L U U
αα
α
α
− + − − − − + −
− + − − +
−
− −
− − − −
− − −=+
−=
+
+ −=
1
,
ij
P
ij
ij NE
NE ij i j
P NW
L
AU L Lα
=+
(2.69)
Está-se, agora, em condições de resolver o sistema de equações com a ajuda desta
factorização aproximada. A expressão responsável pela actualização do valor residual é:
1n n LU Rδ + = (2.70)
onde é designado como a correcção ou actualização do erro de convergência. Asequações são resolvidas como na factorização LU simples. A multiplicação da expressão
(2.70) por L
δ
n
-1 resulta que:
1 1n nU L Rδ ϕ+ −= = (2.71)
nϕ é calculado da seguinte forma:
51
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 68/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
, , , 1, 1 , 1 , 1, 1 , 1,( )i j i j i j i j ij i j i j i j i j i j ij
SW S NW W P R L R L R L R L R Lϕ − − − − + −= − − − − /
1+ −
(2.72)
esta expressão é calculada no sentido crescente dos índices i,j. Uma vez conhecido ϕ é
necessário resolver a seguinte equação:
, , , , 1 , 1, 1 , 1, , 1,i j i j i j i j i j i j i j i j i j i j
N NE E SE R U U U U ϕ δ δ δ δ+ + + += − − − − (2.73)
no sentido decrescente dos índices i,j.
No método SIP, apenas é necessário calcular os elementos das matrizes L e U uma vez,
antes da primeira iteração. Nas iterações subsequentes apenas é necessário calcular o
valor residual R, seguido de ϕ e finalmente δ, por resolução dos dois sistemas
triangulares.
O método de Stone normalmente converge para um número reduzido de iterações. A
rapidez de convergência pode ser melhorada se variarmos α de iteração para iteração (e
de ponto para ponto). Esta alteração permite que a convergência do método em poucas
iterações, mas requer que a factorização seja refeita de cada vez que se altera o
parâmetro α. De forma a não sobrecarregar o computador com os cálculos de L e U é
prudente manter α constante [27].
III. Conclusão
A resolução analítica de uma equação diferencial que rege determinado fenómeno físico
só é possível em casos particulares. Daí a necessidade de reformulação do problema,
dando-lhe uma forma puramente algébrica, processo que é denominado por
discretização.
Para se recorrer aos métodos numéricos é também necessário representar adequadamente
a geometria do domínio de cálculo, o que se consegue fazendo a escolha de um número
finito de nós computacionais, supostos representativos de todo o domínio, onde irá ser
determinada a solução numérica.
52
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 69/134
Capítulo II – Técnicas numéricas de resolução de equações diferenciais
Quando os domínios, físico e computacional, forem idênticos, não é necessário qualquer
transformação, mas no caso de geometrias físicas complexas, a transformação pode
provocar modificações drásticas quer nas equações constitutivas do fenómeno físico quer
nas condições de fronteira associadas.
A dificuldade de descrição dos contornos das formas geométricas complexas pode ser
ultrapassada dividindo a geometria complexa em troços de geometria mais simples, pelo
estabelecimento de malhas curvilíneas em cada um desses troços e pela definição
adequada da interacção entre os mesmos.
O método dos volumes finitos tem como ponto de partida a equação genérica da
conservação. O domínio de cálculo é subdividido num número finito de volumes de
controlo (VC) contíguos, sendo a equação da conservação aplicada a cada VC. No
esquema célula centrada, o centro geométrico de cada VC constitui um nó
computacional onde a variável dependente é avaliada. Para expressar os valores da
variável dependente nas faces dos VC em termos dos valores nos centros dos VC
utilizam-se as técnicas de interpolação. Os integrais de volume e de superfície são
aproximados usando fórmulas de quadratura adequadas.
O método dos volumes finitos é conservativo por natureza, e portanto os integrais de
superfície são idênticos para os VC que partilham a mesma fronteira.
O método iterativo de Stone também designado por SIP (Strongly Implicit Procedure)
tem especial orientação para equações algébricas resultantes da discretização de
equações de diferenciais parciais e normalmente converge para um número reduzido de
iterações, o que permite economizar no tempo de computação.
53
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 70/134
Capítulo III
Formulação matemática da transferência de calor
I. Introdução
O objectivo deste capítulo é o de resolver dois problemas que tendo solução analíticaconhecida, sejam definidos: num primeiro caso, sobre um domínio rectangular de forma
a servir de teste ao algoritmo desenvolvido utilizando coordenadas cartesianas, e no
segundo caso, sobre um domínio circular de forma a servir de teste ao algoritmo
desenvolvido utilizando coordenadas curvilíneas generalizadas.
II. Condução bidimensional num domínio rectangular
Propõe-se, então, a resolução do problema da evolução térmica de uma barra de
comprimento infinito e secção transversal quadrada, constituída por um material de
difusividade unitária, inicialmente à temperatura de 1ºC, isolada em duas faces
contíguas (Oeste e Norte) e colocada em contacto perfeito com um reservatório térmico
a uma temperatura 1ºC nas outras duas faces (Este e Sul).
Este tipo de problema (condução de calor bidimensional) é governado pela seguinte
equação diferencial,
∂
∂+
∂
∂=
∂∂
2
2
2
2
y xa
t
φφφ (3.1)
que quando sujeita às condições iniciais e de fronteira seguintes:
54
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 71/134
Capítulo III – Formulação matemática da transferência de calor
( , , 0) 1
(1, , ) ( , 0, ) 1
(0, , ) ( ,1, ) 0
x y
y t x t
y t x t
x y
φ
φ φ
φ φ
=
= =
∂ ∂= =
∂ ∂
(3.2)
admite a seguinte solução analítica definida por Olson e Schultz [1]:
( )( )2 2 2 2
1 1
(2 1) (2 1)2 1 2 1, , 1 cos sin exp
2 2 4mn
n m
a t n mn m x y t C x y
πφ π π
∞ ∞
= =
− + −− − = + − ∑∑ (3.3)
onde
)12)(12(
)1()1(162
11
−−
−−−=
++
mnC
mn
mn π (3.4)
Utilizando o método dos volumes finitos, os termos presentes na expressão (3.1) são
aproximados da seguinte forma:
( )
1n n
p P
pV V
dV dV V V t t t t
φ φφ
φ φ
+ −∂ ∂ ∂
= ≈ ∆ ≈∂ ∂ ∂ ∆∫ ∫ ∆ (3.5)
2
2ˆ. . e w
e wV S
dV n dS S S x x x x
φ φ φ φ ∂ ∂ ∂ ∂ = = − ∂ ∂ ∂∂ ∫ ∫ =
P W E P e w
c c
S S x x
φ φφ φ −− ∆ ∆
= − (3.6)
o segundo termo do segundo membro da expressão (3.1) é aproximado de modo
análogo.
Num domínio 2D as áreas e volumes utilizadas no método dos volumes finitos são as
que se apresentam no quadro 1.
55
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 72/134
Capítulo III – Formulação matemática da transferência de calor
Quadro 3.1 – Áreas e volumes em 2D utilizadas no método dos volumes finitos [26]
∆V S e=S w S n= S s
2D ∆ x×∆ y ∆ y ∆ x
Substituindo as igualdades presentes no quadro 3.1 nas expressões (3.5) e (3.6)
poderemos escrever a equação (3.1) com o aspecto da expressão (2.38), resultando nos
coeficientes presentes no quadro 3.2.
Quadro 3.2- Coeficientes Ai utilizados no modelo da condução em domínio rectangular
AW
A E
A N
AS
-c
ya
x
∆∆
-c
ya
x
∆∆
-c
xa
y
∆∆
-c
xa
y
∆∆
Os coeficientes Ai são relativos ao nível de tempo n+1, enquanto que Q P é relativo ao
nível de tempo n. Para as restantes formulações de discretização temporal apenas A P e
Q P são alterados (quadro 3.3).
Quadro 3.3 – Coeficientes A P e Q P para diferentes formulações
Formulação A P Q P
Euler explicita x y
t
∆ ∆∆
n
W E N S
n n n n
W W E E N N S S
x y A A A A
t
A A A A
φ
φ φ φ φ
∆ ∆+ + + + ∆
− − − −
P
Crank-Nicolson ( )W E N S
x y A A A A
t
∆ ∆− + + +
∆
n
W E N S
n n n n
W W E E N N S S
x y A A A A
t
A A A A
φ
φ φ φ φ
∆ ∆+ + + + ∆
− − − −
P
Euler implícita ( )W E N S
x y A A A A
t
∆ ∆− + + +
∆ n
I
x y
t φ
∆ ∆∆
Richtmayer–Morton ( ) (1 W E N S
x y A A A A
t β
∆ ∆+ − + + +
∆) 1n n x y
t φ βφ −∆ ∆
−∆
Todas estas expressões são válidas para os CV internos, isto é, exceptuando os volumes
de controlo, em que pelo menos uma das faces contacta a fronteira do domínio.
56
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 73/134
Capítulo III – Formulação matemática da transferência de calor
A implementação das condições fronteira requer que os fluxos nas fronteiras sejam
conhecidos ou possam ser expressos em termos de quantidades conhecidas e dos valores
nodais interiores.
Em geral, numa fronteira de entrada todas as variáveis são prescritas. Se as condições
não forem perfeitamente conhecidas, é usual deslocar-se a fronteira para o mais longe
possível da região de interesse. O fluxo por difusão não é, em geral, conhecido, mas
pode ser aproximado utilizando os valores conhecidos das variáveis na fronteira e
diferenças unilaterais para os gradientes.
Numa fronteira de saída de um determinado escoamento, normalmente, pouco se
conhece. Por essa razão, estas fronteiras devem estar o mais afastado possível da região
de interesse. De outra forma, o erro pode propagar-se ao longo da evolução do fluxo. O
fluxo deve abandonar o domínio de forma directa e ao longo de toda a secção
transversal, e se possível, ser paralelo.
O fluxo por difusão em fronteiras sólidas (paredes) requer atenção especial. Para
quantidades escalares, como a energia térmica, podem ser nulos (paredes adiabáticas),
ser prescritos ou o valor para a variável escalar deve ser prescrito (paredes isotérmicas).Se o fluxo for conhecido, pode ser inserido na equação da conservação para os VC
imediatamente após a fronteira.
No problema que se pretende resolver existem duas fronteiras (Sul e Este) com
temperatura prescrita e de valor igual a 1ºC (paredes isotérmicas), e as fronteiras Norte e
Oeste são adiabáticas pelo que o gradiente na direcção normal é aproximado utilizando
diferenças unilaterais.
Na figura 3.1 mostra-se a malha bidimensional, com cem volumes de controlo, utilizada
para resolver o problema da condução bidimensional em domínio regular pelo método
dos volumes finitos.
57
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 74/134
Capítulo III – Formulação matemática da transferência de calor
0 1
x
y
1
Figura 3.1 - Malha usada para resolver o problema de condução bidimensional usando volumes finitos.
As fronteiras Oeste e Norte são adiabáticas, enquanto que as fronteiras Sul e Este foramcolocadas em contacto perfeito com o meio envolvente.
Na figura 3.2 encontram-se representadas as isotérmicas relativas à solução numérica e
solução analítica proposta por Olson e Schultz quando estão decorridos 0.5 segundos,
utilizando como incremento de tempo ∆t = 0.01 segundos.
0 0.5 1X
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Y
0 . 8 7 1 1
0 . 8 8
8 3
0 . 9 9
1 4
0 . 9 7
4 2
0 . 9 5
7 0
0 . 9 3
9 9
0. 9 2
2 7
0 . 9 0
5 5
0 0.5 1X
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Y
0 . 9 7
4 2 0 0 .
9 5 7 0
0 0. 9 3
9 8 0 0 .
9 2 2 5
9 0 . 9 0
5 3 9
0 . 8 8
8 1 9
0 . 8 7 0
9 9
0 . 9 9
1 4 0
0 0.5 1X
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Y
0 . 8 8
8 3 0
0 . 9 9
1 4 1
0 . 9 7
4 2 2 0 .
9 5 7 0
4 0. 9 3
9 8 5 0.
9 2 2 6 7
0 . 9 0
5 4 8
0 . 8 7
1 1 1
Figura 3.2 – Solução numérica (esquerda) e solução analítica de Olson e Schultz (direita).
A análise das isotérmicas permite verificar a excelente concordância por parte da
solução numérica (esquerda) em relação à solução analítica (direita).
58
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 75/134
Capítulo III – Formulação matemática da transferência de calor
Na figura 3.3 apresenta-se uma comparação entre diferentes formulações e a solução
analítica, para os valores da temperatura pertencentes à diagonal principal quando estão
decorridos 0,1 segundos.
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 1
X [m]
T e m p e r a t u r a [ º C ]
S.Analítica
Euler Implícita
Crank-Nicolson
Richtmayer-Morton
F igura 3.3 – Comparação entre diferentes formulações (Crank-Nicolson, Euler implícita,
Richtmayer -Morton), e a solução analítica de Olson e Schultz para 0.1 segundos.
Uma análise cuidada da figura 3.3 permite constatar que as formulações implícita de
Euler e de Richtmayer-Morton apresentam valores coincidentes e mais próximos da
solução analítica que a formulação de Crank-Nicolson.
0.00
0.100.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 1
X[m]
T e m p e r a t u r a [ º C ]
S.Analítica
Euler Implícita
Crank-Nicolson
Richtmayer-Morton
Figura 3.4 – Comparação entre diferentes formulações (Crank-Nicolson, Euler implícita,
Richtmayer-Morton), e a solução analítica de Olson e Schultz para 0.3segundos.
59
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 76/134
Capítulo III – Formulação matemática da transferência de calor
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 1
X [m]
T e m p e r a t u r a [ º C ]
S.Analítica
Euler implícita
Crank-Nicolson
Richtmayer-Morton
Figura 3.5 – Comparação entre as diferentes formulações (Crank-Nicolson, Euler implícita, Richtmayer
e Morton), e a solução analítica de Olson e Schultz para 0.5 segundos.
Das últimas duas figuras 3.4 e 3.5 facilmente se constata que à medida que o tempo
evolui as diferentes soluções numéricas se aproximam da solução analítica. No entanto,
o esquema implícito de Richtmayer-Morton apresenta um erro menor.
Tal como já se teve oportunidade de constatar o esquema de Richtmayer-Morton utiliza
um factor de relaxação (β), que pode variar entre zero e a unidade pelo que interessa
comparar a solução numérica que se obtém variando este factor. É o que se apresenta nafigura 3.6 quando estão decorridos 0.5 segundos.
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
1 0.95 0.85 0.75 0.65 0.55 0.45 0.35 0.25 0.15 0.05 0
X [m]
T e m p e r a
t u r a [ º C ]
S.Analítica
Beta =1.0
Beta=0.5
Beta=0.0
Figura 3.6 – Comparação entre soluções: numérica (Richtmayer e Morton) utilizando diferentes valores
do parâmetro β e a solução analítica de Olson e Schultz para 0.5 segundos.
60
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 77/134
Capítulo III – Formulação matemática da transferência de calor
Como facilmente se constata a solução mais próxima da analítica é a que se obtém
utilizando β =1.0, pelo que será este utilizado em simulações futuras.
De referir que utilizando um incremento de tempo menor (por exemplo, ∆t = 0.001segundos) as diferentes soluções numéricas apresentam valores sensivelmente
coincidentes e mais próximos da solução analítica.
Tal valor para o incremento de tempo torna possível a utilização da formulação explícita
de Euler, que como se sabe é condicionalmente estável sendo a condição de estabilidade
para este problema dada por [27]:
2( )
2
xt
k
ρ ∆∆ < (3.7)
De modo a evitar esta restrição, em termos de incremento de tempo, esta formulação
não será utilizada em futuras aplicações.
II. Problema da condução bidimensional em domínio circular
Propõe-se nesta secção a resolução do problema da evolução térmica de uma barra de
comprimento infinito e secção recta circular, constituída por um material de
difusividade unitária, inicialmente à temperatura de 0ºC, e colocada em contacto
perfeito com um reservatório térmico a uma temperatura 1ºC.
Este tipo de problema (condução de calor bidimensional) é governado pela equação
diferencial de derivadas parciais seguinte:
∂
∂+
∂
∂=
∂∂
2
2
2
2
η
φ
ξ
φφa
t (3.8)
esta equação quando sujeita às seguintes condições iniciais e de fronteira:
0)0,,( =ηξφ (3.9)
61
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 78/134
Capítulo III – Formulação matemática da transferência de calor
1),1,(),0,(),,1(),,0( ==== t t t t ξφξφηφηφ (3.10)
admite a seguinte solução analítica definida em coordenadas polares [1]:
( ) ( )∑∞
=
−−=1
2
1
0 exp1n
n
n
n t )( J
r)( J r,t φ µ
µ
µ (3.11)
onde µ1, ....., µn são as raízes da equação J 0( µ) = 0, e J 0 e J 1 são as funções de Bessel do
primeiro grau, de primeira e segunda ordem, respectivamente.
Quadro 3.4 – Coeficientes Anb para o modelo da condução em geometria circular
A E 1
1 112 C
J a C C ξ ξξ
− ∆ ∆− − ∆
A NE 1
12 4 C
J aC ηξ
− ∆−∆
AW 1
1 112 C
J a C C ξ ξ
ξ− ∆ ∆
− ∆
ASE 1
12 4 C
J aC ηξ
− ∆∆
A N 1
2 222 C
J a C C η η
η− ∆ ∆
− − ∆
A NW 1
12C
J aC ηξ
− ∆∆
AS 1 2 222 C J a C C
η ηη
− ∆ ∆− ∆ ASW
112 4 C
J aC ηξ
− ∆− ∆
Quadro 3.5 – Coeficientes A P e Q P para diferentes formulações para o modelo da condução em geometriacircular
Formulação A P QP
Crank-Nicolson1
( )2 E W N S J A A A
t
ξ η∆ ∆− + + +
∆ A
( )
1( )
2
1
2
n
E W N S
n n n n
E E W W N N S S
J A A A At
A A A A
ξ ηφ
φ φ φ φ
∆ ∆ − + + + ∆
− + + +
P
Euler implícita ( ) p E W A J A A A A
t
ξ η∆ ∆= − + + +
∆ N S n
pt
J φηξ
∆
∆∆
Richtmayer–Morton ( ) ( )1 E W N S
J A A At
ξ ηβ
∆ ∆+ − + + +
∆ A ( ) 11 n n
p P J t
ξ ηβ φ β φ −∆ ∆
+ −∆
No quadro 3.4 apresentam-se os coeficientes Ai resultantes do processo de discretização
utilizando o método dos volumes finitos.
62
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 79/134
Capítulo III – Formulação matemática da transferência de calor
Tal como no exemplo anterior foram utilizadas diferentes formulações na aproximação
do termo transitório cujas expressões resultantes se apresentam no quadro 3.5.
Na figura 3.7 mostra-se a malha bidimensional obtida por interpolação bilinear transfinita a partir de pontos definidos sobre a fronteira do domínio.
a
y
x
ξ
η
(b)
1
1
0
0
ξ
η
Figura 3.7 - Malha utilizada para resolver o problema de condução bidimensional em domínio
irregular: a) domínio físico; b) domínio computacional.
Compara-se, na figura 3.8, os resultados entre a solução analítica calculada pela
expressão (3.8) e a solução numérica obtida pelo método dos volumes finitos com
formulação generalizada para 0.01 e 0.2 segundos.
Pela análise da figura 3.8 pode-se concluir que após decorridos 0,01 segundos existe
uma excelente correlação entre a solução analítica e numérica, ao passo que decorridos
0,2 segundos a correlação já não é tão satisfatória. Para melhor comparar ambas as
soluções, apresenta-se na figura 3.9 a sua evolução ao longo do tempo em pontos
situados sobre o raio 0.93 m.
63
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 80/134
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 81/134
Capítulo III – Formulação matemática da transferência de calor
Como se pode constatar a evolução ao longo do tempo da solução numérica acompanha
satisfatoriamente a solução analítica, sendo o erro máximo de, sensivelmente, 8% que se
verifica quando decorridos 0.2 segundos, situação ilustrada na figura 3.9.
Refira-se que para efectuar estas simulações foram testados três métodos iterativos de
resolução do sistema de equações resultante do processo de discretização. O
desempenho de cada um dos métodos é apresentado no quadro 3.6 onde se fez, ainda, a
comparação do número de iterações necessário à convergência utilizando dois critérios
de paragem distintos.
Quadro 3.6- Performance dos métodos iterativos TDMA, Jacobi,Gauss-Seidel e SIP
Critério de paragem 10-3 Critério de paragem 10-5 Método
Nº de iterações Resíduo Nº de iterações Resíduo
TDMA - Divergente - Divergente
Jacobi 8 2.56×10-4 1000 2.03×10-5
Gauss-Seidel 5 1.74×10-4 1000 2.03×10-5
SIP 2 9.55×10-4 3 1.92×10-6
O método TDMA (Tridiagonal Matrix Algorithm) revela-se inapto para resolver o
sistema de equações, visto ser divergente. No que concerne ao método Gauss-Seidel,
utilizando o critério de paragem de 10-3, a solução é convergente embora não tão rápida
como o método SIP. Quando se aumenta a precisão para 10-5, apenas o método SIP,
desenvolvido neste trabalho para o caso de uma célula computacional de nove pontos,
responde satisfatoriamente.
De modo a avaliar a diferença que existe no valor da variável dependente (temperatura)
quando obtida pelo método iterativo SIP utilizando critérios de paragem distintos
(resíduos de 10-3 e 10-5), apresenta-se na figura 3.10 a diferença de temperatura que se
verifica para entre ambas as situações.
Como se pode constatar a máxima diferença de temperatura entre as duas abordagens é
da ordem de 10-7ºC. Não existe, portanto, razão para a utilização do critério de paragem
65
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 82/134
Capítulo III – Formulação matemática da transferência de calor
10-5. A utilização do critério de paragem 10-3 tem a vantagem de reduzir o tempo de
computação necessário, não afectando significativamente a solução final.
0.0E+00
5.0E-08
1.0E-07
1.5E-07
2.0E-07
2.5E-07
3.0E-07
1.0 1.0 1.0 0.8 0.6 0.3 0.0 0.3 0.6 0.8 1.0 1.0 1.0Raio [m]
D i f e r e n ç a d e t e m p e r a t u r a [ º C ]
Figura 3.10 - Diferença de temperatura sobre o raio 0.93 m obtida através do método iterativo SIP
utilizando dois critérios de paragem distintos (10-5
e 10-3 ) .
III. Conclusão
O código desenvolvido foi testado, numa primeira abordagem, com a solução analítica proposta por Olson e Schultz para a condução de calor num corpo bidimensional
utilizando coordenadas cartesianas e revelou ter um comportamento bastante
satisfatório.
Nesta aplicação numérica teve-se oportunidade de determinar qual das formulações de
discretização temporal testadas (Euler implícita, Crank-Nicolson e Ritchmayer-Morton)
apresenta melhores resultados comparativamente com a solução analítica. A formulação
de Ritchmayer-Morton mostrou ser a que mais se aproxima da solução analítica.
De seguida comparou-se as soluções obtidas fazendo variar o factor de relaxação (β),
presente na formulação de Ritchmayer-Morton, entre zero e a unidade, tendo a melhor
solução sido obtida para o valor β=1. A solução menos precisa foi obtida para β=0
situação para a qual a esta formulação é equivalente à formulação implícita de Euler.
66
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 83/134
Capítulo III – Formulação matemática da transferência de calor
Numa segunda abordagem, o código desenvolvido utilizando coordenadas curvilíneas
foi testado por comparação com uma solução analítica conhecida para o problema da
condução de calor em domínio circular tendo apresentado resultados satisfatórios,
embora tenham sido necessárias algumas alterações em relação à primeira abordagem
que importa referir.
A primeira alteração está relacionada com o método iterativo utilizado. O método
TDMA, utilizado na resolução numérica do problema da condução de calor em domínio
rectangular com bom desempenho, revelou-se inapto para a resolução do mesmo
problema em domínio circular. Foi, então, usado o método iterativo SIP, que mostrou
ser um método iterativo de rápida convergência e elevada precisão.
De seguida fez-se variar o critério de paragem deste método iterativo e determinou-se a
diferença que existe entre ambas as soluções (campos de temperatura). Tal diferença de
temperatura é irrisória, concluindo-se que o uso do critério de paragem 10-3 não traz
prejuízo significativo para a precisão da solução e tem a vantagem de reduzir o tempo
de computação.
67
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 84/134
CAPÍTULO IV
Simulação numérica da solidificação em fundição
I. Introdução
Neste capítulo será apresentada a simulação numérica da solidificação de uma liga de
alumínio Al 12Si (alumínio com 12% de silício), no interior de uma moldação de ferro
fundido cinzento. Para atingir esse objectivo será em primeiro lugar definido o modelo
matemático adequado, seguindo-se a sua discretização e definição das condições iniciais
e de fronteira para três situações distintas. O passo seguinte passa pela elaboração de um
código utilizando o método dos volumes finitos, do qual se fará uma descrição
resumida. No final do capítulo serão apresentados os resultados para os diferentes tiposde condições fronteira aplicadas às interfaces de contacto entre os diferentes
sub-domínios definidos no sistema peça/moldação, bem como uma discussão
pormenorizada dos resultados.
II. Modelo matemático
A solidificação de materiais pode ser considerada fundamentalmente como um processode transferência de calor em regime transitório. A transformação do material líquido em
material sólido é acompanhada por libertação de energia térmica sensível e latente.
A equação diferencial que traduz o equilíbrio térmico entre difusão, geração (ou
absorção) de calor e acumulação (ou libertação) de calor num corpo pode ser
apresentada estabelecendo o balanço térmico num elemento infinitesimal (dx, dy, dz ), na
vizinhança de um ponto genérico, como representado na figura 4.1.
68
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 85/134
Capítulo IV – Simulação numérica da solidificação em fundição
dx x
x x
k +
∂∂
−φ
z
z z
k
∂∂
−φ
y
y y
k
∂∂
−φ
x
x x
k
∂∂
−φ
dz z
z z
k +
∂∂
−φ
q&
dy y
y y
k
+
∂∂
−φ
t C p ∂
∂φρ
dy
dx
dz
Figura 4.1 - Balanço térmico num elemento infinitesimal.
Cada face do elemento infinitesimal é atravessada por um fluxo proporcional ao
gradiente de temperatura local na direcção normal da mesma. A solidificação só é
conseguida à custa de uma diminuição da energia acumulada no elemento infinitesimal
sob a forma de calor sensível e, principalmente, sob a forma de calor latente.
Como facilmente se depreende este tipo de fenómeno tem um carácter transitório, cuja
formulação analítica é expressa pela equação diferencial, de derivadas parciais, seguinte
[1,5:7].
( )( ).
pC k q
t
ρ φφ
∂= ∇ ∇ +
∂& (4.1)
A primeira derivada do campo de temperaturas (φ) em ordem ao tempo, descreve a sua
evolução temporal, o divergente do gradiente do mesmo campo de temperaturas
descreve a transferência de calor por condução, o termo fonte q correspondente à
libertação ou absorção de calor na mudança de fase e ρ , C
&
p e k são características
físicas do material: massa volúmica, capacidade térmica mássica e condutibilidade
térmica, respectivamente.
69
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 86/134
Capítulo IV – Simulação numérica da solidificação em fundição
Os modelos matemáticos utilizados para descrever o fenómeno da solidificação
consideram geralmente como constantes as características físicas ρ , C p e k do metal,
tomando a expressão (4.1) o seguinte aspecto:
2 2
p p
k qa q
t C C k
φφ φ
ρ ρ∂
= ∇ + = ∇ +∂
&&
a (4.2)
Onde pC k a ρ= , é denominada difusividade, e constitui uma medida da rapidez com
que o calor se transmite através do meio constituído pelo material.
Relativamente ao termo fonte, pode-se constatar que a fonte interna de calor provém do
calor latente libertado por mudança de fase. Esta mudança de fase não é instantânea,ocorre a uma taxa que pode ser expressa em função da fracção efectiva de material
sólido, a fracção sólida f s, da massa volúmica do material ρ, e da variação da entalpia na
transformação de fase ∆h f , (calor latente) [1,5,6]:
t
f hq
s f
∂
∆∂=
)(ρ& (4.3)
Quando se trata de uma liga metálica cuja mudança de fase ocorre num intervalo detemperaturas [φl ,φ s], em que se conhece a evolução da fracção sólida existente, f s, com a
temperatura, pode escrever-se:
t
f
t
f s s
∂∂
∂
∂=
∂
∂ φφ
)()( (4.4)
Assumindo que o material é isotrópico, pode-se escrever a equação (4.2) do seguinte
modo:
21 ( f s
p
h f a
t C
φφ
φ
∆ ∂∂− = ∂ ∂
)∇ (4.5)
Em condições de equilíbrio local, a fracção sólida determina-se, a cada temperatura,
pela regra da alavanca entre as composições média das linhas C , de liquidus C 1 e de
solidus C s:
70
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 87/134
Capítulo IV – Simulação numérica da solidificação em fundição
sl
l
sC C
C C f
−
−= (4.6)
Admitindo uma relação linear entre f s( φ ) e φ, aproximação razoável pelo menos na
vizinhança de cada valor de temperatura, pode considerar-se φ∂∂ )( s
f como constante,
para pequenos intervalos de variação desta variável. Esta assunção permite linearizar o
termo fonte presente na equação (4.5), o que simplifica bastante o problema em termos
de cálculo numérico.
Em presença de um metal puro, ou aquando da ocorrência de solidificação eutética
(figura 4.2), surgem dificuldades pois φ∂)( s f ∂ → ∞
porque pelo menos parte datransformação se dá a temperatura constante, φeut .
Uma forma de ultrapassar este problema, é manter a temperatura do volume elementar
em causa à temperatura constante φeut , só permitindo a sua variação quando a fracção
sólida (ou eutética) atingir um dos limites 0 ou 1 que correspondem ao desaparecimento
de uma das fases.
T e m p e r a t u r a
Tempo
φ l φ
S
L í q u i d o
S ó l i d o
f s=1 , f eut =1
f s=0 l
l
C C f s C C s
−=
−
Líquido + Sólido
φ l
φ e
L í q u i d o
S ó l i d o
Primária Eutética
f s=1 , f eut =1
f s=0
f eut =0
l
l
s s
C C f
C C
−=
−
Tempo
T e m p e r a t u r a
Líquido + Sólido
Solidificação
Figura 4.2 – Ilustração da variação da fracção sólida com a temperatura, no intervalo de solidificação, para a generalidade das ligas (esquerda), e o caso de uma liga em que ocorre solidificação
primária no intervalo [ φl ,φeut ], seguida de solidificação eutética à temperatura φeut (direita).
Outra forma de lidar com este problema é considerar que a solidificação decorre de
facto num intervalo de temperaturas na vizinhança de φeut , embora estreito, e assim
criar um valor finito para φeut , apesar de fictício.
A utilização de coordenadas curvilíneas provoca alterações significativas na equação
71
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 88/134
Capítulo IV – Simulação numérica da solidificação em fundição
que governa o fenómeno da solidificação. Tendo decidido que no domínio transformado
a geometria é uniforme, rectangular e fixa no tempo, os cálculos necessários serão
realizados numa malha uniforme de elementos cúbicos (ou quadrados em 2D). As
coordenadas curvilíneas são definidas pela transformação xi = xi( ξ j ), j=1,2,3, que é
caracterizada pelo jacobiano J :
3
3
3
2
3
1
2
3
1
3
2
2
1
2
2
1
1
1
det
ξ
ξ
ξ
ξξ
ξξ
ξξ
ξ
∂
∂
∂∂
∂∂
∂
∂
∂
∂
∂∂
∂∂
∂∂
∂∂
=
∂∂
=
x
x
x
x x
x x
x x
x J
j
i (4.7)
Tendo em atenção que apenas as derivadas inerentes às variáveis geométricas são
transformadas:
J x x
ij
ji
j
ji
βξφξ
ξφφ
∂∂
=∂
∂
∂∂
=∂∂
(4.8)
onde β
ij
=(-1)
i+j
det(J ij )
representa o cofactor de ji x ξ∂∂ no Jacobiano J .
A equação (4.1) expressa em coordenadas curvilíneas assume, então, o seguinte aspecto
[27]:
q J B J
k
xt
C J mj
m j
p&=
∂∂
∂∂
−∂
∂
ξφφρ )(
(4.9)
onde o coeficiente Bmj se define como:
Bmj = βkj βkm = β1j β1m + β2j β2m + β3j β3m (4.10)
A expressão (4.9) tem a mesma forma que a equação (4.1), mas cada termo desta última
é substituído pela soma de três termos, tal como se mostra na expressão (4.11), onde se
fazem as seguintes substituições: ξ1= ξ , ξ2= η e ξ3= ζ.
72
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 89/134
Capítulo IV – Simulação numérica da solidificação em fundição
11 21 31( ) pC k
J B Bt J
ρ φ φ φ φξ ξ η ζ
∂ ∂ ∂ ∂∂− + + ∂ ∂ ∂ ∂ ∂
B
12 22 32k B B B
J
φ φ φη ξ η ζ
∂ ∂ ∂∂− + + ∂ ∂ ∂ ∂
13 23 33k B B B
J
φ φ φζ ξ η ζ
∂ ∂ ∂∂− + + ∂ ∂ ∂ ∂
& Jq= (4.11)
Efectuando todas as operações de derivação presentes na expressão (4.11), a equação
que caracteriza o fenómeno da solidificação, assume, então, o seu aspecto final:
2 2
1 2 3 11 1221 f s
p
h f J C C C C C
t C
φ φ φ φ φφ ξ η
φζ ξ ηξ
∆ ∂∂ ∂ ∂ ∂ ∂ ∂− = + + + + ∂ ∂ ∂ ∂ ∂ ∂∂ ∂
2 2 2
13 22 23 332 2C C C C
φ φ φξ ζ η ζη ζ∂ ∂ ∂
+ + + +∂ ∂ ∂ ∂∂ ∂
2φ∂ (4.12)
onde os coeficientes C i e C ij resultantes da transformação do domínio físico arbitrário no
domínio computacional (rectangular) são determinados pelas seguintes expressões:
∂∂
+∂
∂+
∂∂
+∂
∂+
∂∂
+∂
∂= −
−−−
ζηξζηξ
131211113
112
111
1
1
)()()( B B B J B
J B
J B
J C (4.13)
∂
∂+
∂∂
+∂
∂+
∂∂
+∂
∂+
∂∂
= −−−−
ζηξζηξ
232221123
122
121
1
2
)()()( B B B J B
J B
J B
J C (4.14)
1 1 1 31 3231 32 33 13
( ) ( ) ( ) J J J B BC B B B J
ξ η ζ ξ η− − − − ∂ ∂ ∂ ∂ ∂ ∂= + + + + + ∂ ∂ ∂ ∂ ∂
33 B
ζ∂ (4.15)
1 1111C J B
−= (4.16)
( )12211
12 B B J C += − (4.17)
( )1 31 1313C J B B−= + (4.18)
73
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 90/134
Capítulo IV – Simulação numérica da solidificação em fundição
1 2222C J B−= (4.19)
( )1 32 2323C J B B−= + (4.20)
1 3333C J B−= (4.21)
O aparecimento de derivadas mistas deve-se à utilização de malhas não ortogonais, que
como se mostra em (4.11) estão multiplicadas pelos coeficientes Bmj de índices
distintos, que se tornam nulos se forem utilizadas malhas ortogonais rectilíneas ou
curvilíneas.
III. Discretização
Pretende-se com a operação de discretização, substituir o modelo contínuo, expresso por
(4.12), por um modelo discreto aproximado de forma algébrica genérica,
P P nb nb A Aφ φ= +∑ Q (4.22)
onde o somatório indica a soma de todos os nós vizinhos (nb) do nó central P .
Utilizando o teorema da divergência, a regra do ponto médio para aproximação do
integral de superfície resultante e a interpolação linear, resulta para a primeira derivada
segundo a direcção ξ [5,6],
2 E W
V
dV V φ φφ
ξ
−∂= ∆
∂∫ (4.23)
para a segunda derivada da variável dependente segundo a direcção ξ resulta [5,6,21];
2
2ˆ. .
V S
dV n dS φ φ
ξξ∂ ∂
= =∂∂∫ ∫
ζηξ
φφξ
φφξφ
ξφ
∆∆
∆
−−
∆−
=
∂∂
−
∂∂
=c
W P
c
P E w
w
e
e
S S (4.24)
74
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 91/134
Capítulo IV – Simulação numérica da solidificação em fundição
por fim a derivada mista é aproximada do seguinte modo [5,6,21];
≈
∂∂
−
∂∂
=
∂∂
∂∂
=∂∂
∂∫ ∫ w
w
e
eV V
S S dV dV ηφ
ηφ
ηφ
ξηξφ2
ςηξ
φφφφ
ξ
φφφφ∆∆
∆
−−+−
∆
−−+≈
c
SW S NW N
c
SE S NE N
44 (4.25)
No quadro 4.1 apresentam-se os coeficientes Ai, resultantes do processo de discretização
enquanto que no quadro 4.2 se apresentam os coeficientes A P e Q para os diferentes
métodos de discretização temporal tratados.
Quadro 4.1 – coeficientes Ai da equação constitutiva do fenómeno da solidificação.
A E 1 11
2 c
C C ξ ξξ
∆ ∆−
∆ A NE
c
C
ξη
∆∆
−4
12
AW 1 11
2 c
C C ξ ξξ
∆ ∆− −
∆ ASE
c
C
ξη
∆∆
412
A N
2 22
2 c
C C η η
η
∆ ∆
− − ∆
A NW c
C
ξ
η
∆
∆
4
12
AS 2 22
2 c
C C η ηη
∆ ∆−
∆ ASW
c
C
ξη
∆∆
−4
12
Quadro 4.2 – Coeficiente A P e Q da equação constitutiva do fenómeno da solidificação para diferentesmétodos de discretização de tempo.
Formulação A P Q
Euler implícita ( )t W E N S K A A A A− + + + n
t K φ
Crank-Nicolson ( )0.5t W E N K A A A− + + +
S A b ( )( )
8
1
0.5 0.5n n
t W E N S P nb n
nb
K A A A A Aφ φ=
+ + + + − ∑
Richtmayer - Morton ( ) (1 0.5t W E N K A A Aβ+ − + + + )S A 1(1 ) 0.5n n
t K β φ βφ −+ −
75
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 92/134
Capítulo IV – Simulação numérica da solidificação em fundição
Por uma questão de simplificação das expressões apresentadas no quadro 4.2,
considerou-se a variável K t , cujo valor é dado pela expressão 4.26.
1 f s
t
p
h f J K a C t
ξ ηφ
∆ ∂ ∆ ∆= − ∂ ∆ (4.26)
Todas as expressões apresentadas nos quadros 4.1 e 4.2 são válidas apenas para os VC
internos, isto é, exceptuando os volumes de controlo das fronteiras. Para os VC
pertencentes às fronteiras é necessário efectuar algumas alterações, é o que se verá na
próxima secção.
IV. Condições iniciais e de fronteira
O material líquido, ao ser vazado, entra em contacto com as paredes internas da
moldação criando-se uma interface metal/moldação com uma determinada resistência
térmica. Esta resistência térmica é decorrente de vários factores, tais como: a afinidade
físico-química entre material da moldação e material a ser solidificado não é perfeita
consequentemente a parede interna da moldação não será completamente molhada pelo
líquido; a rugosidade interna da moldação cria uma micro geometria superficial que
propicia o surgimento de apenas alguns pontos de contacto intercalados por regiões de
separação física metal/moldação.
A utilização de lubrificantes para facilitar a desmoldagem provoca a formação de uma
película de separação entre material e moldação; após a formação de uma certa camada
solidificada a contracção natural do material solidificado provoca um aumento na
separação física entre o material sólido e a moldação. Nessas condições, a transferênciade calor na interface material/moldação dá-se por condução nos pontos de contacto e
através dos gases aprisionados pelos espaços criados, e também por convecção e
radiação entre as duas superfícies separadas. Newton então propôs que essas superfícies
de separação fossem perfeitamente planas e paralelas, sendo esse espaço preenchido por
um certo gás.
Seria ideal que a superfície dos fundidos fosse perfeitamente lisa e livre de
irregularidades, mas na prática não é isso que acontece, apresentam vários graus de
76
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 93/134
Capítulo IV – Simulação numérica da solidificação em fundição
rugosidade, devido à condição superficial da moldação, a fenómenos ligados à
temperatura (expansão térmica dos materiais), à tensão de escoamento e a reacções
químicas. A rugosidade superficial é, portanto, controlada primeiramente pela escolha
dos materiais, pela condição da superfície da moldação, pela qualidade do acabamento
superficial desejado para a futura peça e pelas temperaturas da moldação e de
vazamento. Além disso, o escoamento turbulento durante o enchimento da moldação
pode, também, conduzir a várias imperfeições superficiais.
Pode-se então afirmar que o fluxo de calor através de uma interface será a resultante da
combinação de vários modos de transferência que podem coexistir:
o Condução no estado sólido através dos pontos de contacto;
o Condução através dos gases que ocupam a folga que se forma entre a peça e a
moldação no decurso do processo;
o Radiação entre as superfícies da peça e da moldação através da mesma folga.
Estes modos de transferência são susceptíveis de sofrer alterações durante o processo, e
serão ligeiramente diferentes de vazamento para vazamento. O valor do coeficiente de
transferência vai então variar com diversos factores [1]:
o Espessura e condutividade térmica da camada de desmoldante;
o Forma e dimensões da peça e da moldação;
o Geometria da folga entre as partes em contacto, dependente das
características físicas dos materiais (nomeadamente dilatações de natureza
térmica);
o Estados de acabamento superficial das diferentes interfaces;
o Estado de agitação do ar ambiente (no caso da transferência para o ambiente).
Os coeficientes de transferência de calor nas interfaces foram seleccionados a partir dos
valores indicados por Monteiro [1] e encontram-se compilados nos quadros 4.3 e 4.4.
77
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 94/134
Capítulo IV – Simulação numérica da solidificação em fundição
Quadro 4.3 – Propriedades físicas da liga metálica ensaiada e do material da moldação [1].
Liga metálica MoldaçãoPropriedade
Al 12Si Ferro Fundido Cinzento
Massa Volúmica [kg/m3] 2670 7230
Condutividade térmica [W/m°C] 185 38
Capacidade térmica mássica [J/kg°C] 1260 750
Calor latente de transformação de fase [kJ/kg] 395 -
Temperatura de liquidus [°C]
Temperatura de solidus [°C]
585
575
-
-
Quadro 4.4 – Coeficientes de transferência de calor nas diferentes interfaces [1].
INTERFACE COEFICIENTE DE TRANSFERÊNCIA DE CALOR
NEWTONIANA [W/m2°C]
D1/D2 hi =2500
D1/D3 hi =2500
D1/D4 hi =2500
D2/D3 hi =500
D3/D4 hi =600
D2/Ambiente ha =150
D3/Ambiente ha =150
Para as fronteiras dos sub-domínios da peça que contactam exclusivamente com
fronteiras de sub-domínios da peça, assim como para as fronteiras dos sub-domínios da
moldação que contactam exclusivamente com fronteiras de sub-domínios da moldação
foram assumidas condições de contacto perfeito expressas matematicamente pela
expressão 4.27 [1].
78
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 95/134
Capítulo IV – Simulação numérica da solidificação em fundição
21 mmnn
=
∂∂φ
∂∂φ
e φ = 21 mm φ (4.27)
Para a interface entre o metal e a moldação, numa primeira abordagem, será considerada
uma igualdade de fluxos de calor por condução. Na segunda abordagem será
considerado um fluxo proporcional à diferença entre as temperaturas na parede exterior
do provete, φ s, e na parede interior da moldação, φ
m, sendo o coeficiente de
proporcionalidade dado pelo coeficiente de transferência de calor Newtoniana, hi:
s
s
m
m
n
k
n
k
=
∂
∂φ
∂
φ∂ (4.28)
( )m s i m
m s
k k hn n
∂ φ ∂φφ φ
∂ ∂
= =
s− (4.29)
finalmente para a fronteira exterior do sistema, foi especificado um fluxo proporcional
à diferença entre a temperatura na parede exterior da moldação, φm, e a temperatura
ambiente φa:
( )ama
m
m hn
k φφ∂∂φ
−=
(4.30)
A temperatura de liquidus será utilizada como condição inicial para os domínios que
constituem a peça enquanto que para a moldação fechada e vazia se simulou o seu
arrefecimento a partir de uma temperatura inicial φo=300°C e para um período de 60
segundos de modo a obter um campo de temperaturas mais condizente com a realidade.
Nesta simulação foram considerados negligenciáveis os fluxos de calor através da
interface D1/D2, D1/D3 e D1/D4 (ver figura 4.4), tendo estas fronteiras sido
consideradas adiabáticas. Esta assunção é justificada pela quase inexistência de
circulação de ar na cavidade e pela reduzida capacidade térmica do ar que ocupa a
cavidade.
79
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 96/134
Capítulo IV – Simulação numérica da solidificação em fundição
V. Domínio de cálculo
A simulação numérica foi realizada sobre um domínio bidimensional constituído pela
secção transversal definida pelo plano médio da peça, como indicado na figura 4.3.
Figura 4.3 - Dimensões da secção transversal do conjunto moldação/peça.
Para efeito de simulação numérica, consideraram-se como distintos os domínios
constituídos pelas partes da moldação, pela peça e pelo meio ambiente.
O domínio D1 (figura 4.4) é constituído pela secção transversal da peça, incluindo o
canal de enchimento e o ataque. O domínio D2 é constituído pela secção transversal da parte inferior da moldação, enquanto que a parte superior da moldação é constituída
pelos domínios D3 e D4, tal como se pode constatar pela figura 4.4. O domínio D4
corresponde ao inserto necessário para produzir a forma reentrante na peça, e está
permanentemente fixo ao domínio D3. O ambiente constitui um domínio que envolve
exteriormente o conjunto, que contacta apenas com os domínios D2 e D3.
D3
D1 D4
D2
Figura 4.4- Secção transversal do conjunto moldação/peça mostrando a divisão do domínio
80
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 97/134
Capítulo IV – Simulação numérica da solidificação em fundição
Cada um dos domínios D1, D2 e D3 foi depois subdividido em sub-domínios mais
simples, para facilitar a geração da malha, como se ilustra na figura 4.4 e 4.5. As
fronteiras dos diferentes domínios foram organizadas de forma que a condição de
fronteira estabelecida entre sub-domínios de cada domínio principal corresponde à
situação de contacto perfeito.
2
14
6
4
9
1
12 13 15 16
853
10
11
17
7
Figura 4.5 – Ilustração dos 17 sub-domínios em que foi subdividido o conjunto peça/ moldação de modoa facilitar a geração da malha.
Exceptuando o domínio D4, correspondente ao inserto, todos sub-domínios foram
subdivididos para facilitar a obtenção da malha por interpolação bilinear tal como se
apresenta na figura 4.6.
Figura 4.6 - Representação da malha utilizada para discretizar cada um dos sub-domínios que
compõem o conjunto peça/moldação.
VI. Estrutura do código desenvolvido
No quadro 4.5 apresenta-se uma lista dos programas comerciais mais importantes
usados actualmente na simulação do fenómeno da solidificação. Não é de surpreender
que muitas das técnicas usadas nestes softwares não sejam descritas de modo a evitar a
sua duplicação. Tal situação leva, como se teve oportunidade de constatar no primeiro
81
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 98/134
Capítulo IV – Simulação numérica da solidificação em fundição
capítulo, a que grande parte dos investigadores que trabalham na área da solidificação
dos metais optem por desenvolver os seus próprios códigos.
As razões que podem ser invocadas estão essencialmente relacionadas com os seguintesfactores:
o Especificidade do assunto que se pretende tratar;
o Nível de informação fornecida de modo a se poder confiar nos resultados
obtidos;
o Rapidez na obtenção da solução.
Quadro 4.5 – Software comercial [31:35]
Designação Organização Características
MAGMA-FLOW Magmasoft, Dinamarca FDM – 3D
FLOW – 3D Flow Science Inc., USA FDM – 3D
ProCAST Universal Energy Systems, USA FEM – 3D
SIMULATOR Aluminium Pechiney, França FVM – 3D
Constata-se, pelo quadro 4.5 que todos os métodos de discretização, método das
diferenças finitas (FDM), método dos elementos finitos (FEM) e método dos volumes
finitos (FVM), são utilizados nos códigos comercias para o estudo numérico do
fenómeno da solidificação, tal como, aliás, já se tinha constatado no primeiro capítulo
relativamente aos vários investigadores.
O código desenvolvido neste trabalho, cujo organograma se apresenta na figura 4.7,
utiliza o método dos volumes finitos com formulação em coordenadas curvilíneas e
comporta vários subprogramas que se podem descrever do seguinte modo:
Malha computacional – trata-se do subprograma que permite definir o domínio
computacional, para o qual se optou por um quadrado com 49 volumes de controlo (sete
em cada direcção). Este domínio é único e os seus parâmetros constituem a base para a
definição da malha utilizada para descrever a geometria real.
82
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 99/134
Capítulo IV – Simulação numérica da solidificação em fundição
Malha Real – é o subprograma que permite a definição da malha do domínio real
(conjunto peça/moldação). Como a peça possui geometria complexa, a moldação,
possuindo uma forma negativa da peça, terá também uma forma complexa. De modo a
facilitar a definição da malha o conjunto peça/moldação foi dividido em 17
sub-domínios e, de seguida foram definidos os seus contornos, bem como o número de
volumes sobre cada fronteira de acordo com o domínio computacional. O formato dos
volumes de controlo internos é obtido por interpolação bilinear, realizado em sub rotina
própria.
Interpolação bilinear - uma vez definidos os contornos de cada sub domínio do conjunto
peça/moldação é necessário preencher o seu interior com volumes de controlo, cujo
formato é definido explicitamente interpolando algebricamente as posições que
assumem sobre as fronteiras dos diferentes sub-domínos, tal como foi descrito no
segundo capítulo.
Transformação de coordenadas – é o subprograma que permite a transformação de
coordenadas cartesianas em coordenadas curvilíneas. Calcula o Jacobiano, o
determinante do Jacobiano e posteriormente as constantes C i e C ij presentes na equação
constitutiva da solidificação.
Propriedades físicas – é sabido que a liga de alumínio vazada tem como temperatura de
liquidus 585 ºC e de solidus 575ºC, portanto neste intervalo de temperaturas existe
libertação de calor sob duas formas distintas: calor sensível e calor latente. Assim que a
temperatura desça abaixo da temperatura de solidus termina a libertação de calor latente
o que provoca alterações significativas no modelo numérico. De um modelo de
mudança de fase (onde existe libertação de calor sensível e latente) passa-se a um
modelo de arrefecimento puro onde apenas existe libertação de calor sensível.
Esta sub rotina determina, em cada iteração, o valor máximo de temperatura dos
diferentes sub-domínios que constituem a peça, compara-o com a temperatura de
solidus e, quando inferior, elimina o calor latente.
Discretização espacial – é o subprograma que executa a discretização das diferenciais
espaciais utilizando o método dos volumes finitos.
83
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 100/134
Capítulo IV – Simulação numérica da solidificação em fundição
Discretização temporal – é o subprograma que executa a discretização do diferencial de
tempo presente na equação constitutiva da solidificação utilizando o método dos
volumes finitos para três formulações distintas: implícita de Euler, Crank-Nicolson e
Richtmayer-Morton.
Condições iniciais – é o subprograma que permite obter o campo inicial de temperaturas
na moldação, refira-se que este sub programa faz uso da totalidade das sub rotinas
descritas até aqui e, ainda, de uma sub rotina onde se definem as condições de fronteira
referentes à situação de moldação fechada e vazia, onde existem fronteiras isoladas
entre si, contacto perfeito e fluxo de calor para o ambiente.
Condições de fronteira – refere-se este subprograma à definição das condições defronteira para a situação pós vazamento do metal líquido na moldação. Onde se incluem
fronteiras de coeficiente de transferência de calor Newtoniano, contacto perfeito e
fronteiras cujo fluxo depende apenas da condutividade térmica dos materiais em
contacto.
SIP (Strongly Implicit Procedure) – é o subprograma que permite a resolução via
iterativa do sistema de equações algébricas resultante do processo de discretização, tal
como foi descrito no segundo capítulo.
Print – é o subprograma onde são criados ficheiros de dados (um para cada
sub-domínio), que serão utilizados no pós-processamento.
Toda a programação foi realizada em microcomputador utilizando a linguagem de
programação Fortran, sendo o pós-processamento dos dados obtidos pelo código
desenvolvido realizado com o Software Tecplot.
84
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 101/134
Capítulo IV – Simulação numérica da solidificação em fundição
Sub programa
Malha
computacional
Sub programa
Malha
real
Sub programa
Interpolação
bilinear
Sub programa
Discretização
espacial
Sub programa
Discretização
temporal
Sub programa
Condições
iniciais
Sub programa
Condições
fronteira
Sub programa
Sub programa
SIP
Sub programa
Propriedades
físicas
Sub programa
Transformação
coordenadas
Programa
Solidificação
Figura 4.7 – Organograma do código desenvolvido utilizando o método dos volumes finitos
85
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 102/134
Capítulo IV – Simulação numérica da solidificação em fundição
VII. Apresentação e discussão dos resultados
Nesta secção serão apresentados os resultados obtidos através código desenvolvido
utilizando o método dos volumes finitos com formulação em coordenadasgeneralizadas. O pós-processamento dos dados aqui apresentados foi realizado com a
aplicação Tecplot que permite a visualização simultânea das linhas de temperatura
constante (isotérmicas) e do campo de temperaturas mediante diferentes cores.
Tal como foi referido aquando da definição das condições iniciais, o campo inicial de
temperaturas na moldação é obtido simulando o arrefecimento desta durante 1 minuto
na condição de fechada e vazia. Essa simulação foi efectuada para duas situações
distintas: primeiro, considerando que existe um contacto perfeito entre a parte superior einferior da moldação (interface D2/D3) e também na interface D4/D3; e segundo,
considerando a existência de um fluxo de calor proporcional à diferença de temperaturas
nas interfaces D2/D3 e D4/D3.
O campo de temperaturas para o conjunto peça/moldação será apresentado para três
situações distintas de acordo com o tipo de condições impostas às interfaces dos
diferentes domínios. Na primeira abordagem considera-se que o fluxo de calor através
das referidas interfaces é apenas de difusão, na segunda abordagem considerou-se um
coeficiente de transferência de calor Newtoniano na interface peça/moldação, e na
terceira e última abordagem foi considerado o mesmo tipo de coeficiente de
transferência de calor também nas interfaces D2/D3 e D3/D4.
VII.1. Campo de temperaturas inicial na moldação
Na figura 4.8 encontra-se representado o campo de temperaturas inicial, considerando
contacto perfeito entre a parte superior e inferior da moldação (interface D2/D3) e
também na interface D4/D3. Nota-se que a zona correspondente ao domínio D4 (ver
figura 4.4) se encontra a uma temperatura mais elevada que os restantes domínios dado
possuir três fronteiras adiabáticas e se encontrar mais afastada do meio ambiente. Os
cantos superiores da moldação apresentam valores de temperatura mais baixos que os
homólogos inferiores, o que é perfeitamente realista devido à sua menor dimensão.
86
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 103/134
Capítulo IV – Simulação numérica da solidificação em fundição
FI
290.4
288.5
286.5
284.6
282.7
280.7
278.8276.9
274.9
273.0
271.1
269.1
267.2
265.3
263.4
Figura 4.8 – Campo de temperaturas inicial [ºC] para a moldação obtido para uma temperatura inicial de 300ºC, considerando adiabáticas as interfaces peça/moldação, contacto perfeito entre asinterfaces moldação/ moldação e transferência de calor por convecção entre a superfície externada moldação e o ambiente.
FI
298.0
295.7
293.2
290.7
288.2
285.7
283.1
280.6
278.1
273.1
270.6
268.1
265.6
263.0
260.5
Figura 4.9 – Campo de temperaturas inicial [ºC] para a moldação obtido para uma temperatura inicial de 300ºC, considerando adiabáticas as interfaces peça/moldação, transferência de calor
Newtoniana nas interfaces dos domínios D2/D3 e D4/D3, transferência de calor por convecçãoentre a superfície externa da moldação e o ambiente e contacto perfeito nas interfaces moldação/ moldação.
Na figura 4.9 encontra-se representado o campo de temperaturas inicial obtido
considerando transferência de calor Newtoniana nas interfaces D2/D3 e D4/D3. A
diferença em relação ao campo de temperaturas da figura 4.8 está na descontinuidade
nos valores de temperatura que se verificam nas interfaces onde foi especificado oreferido fluxo.
VII.2. Resultados considerando transferência de calor por condução
Nesta abordagem do problema será considerado que o fluxo de calor existente na
interface peça/moldação é apenas de difusão, ou seja, é proporcional à diferença de
temperatura que se verifica na referida interface, sendo o coeficiente de
87
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 104/134
Capítulo IV – Simulação numérica da solidificação em fundição
proporcionalidade dado pela relação das condutividades térmicas dos materiais em
contacto.
Apresenta-se na figura 4.10 o campo de temperaturas para o conjunto peça/moldaçãodecorridos 0.25 segundos após o vazamento. Verifica-se que grande parte da peça se
encontra ainda muito próximo da temperatura de liquidus, é também visível o aumento
de temperatura nas zonas da moldação mais próximas da peça.
FI
584.8
580.7
556.5
462.7
310.0
220.5
209.9
207.9
197.8
170.4
157.9
148.1
144.7
141.8
133.0
Figura 4.10 – Campo de temperaturas no conjunto peça/moldação decorridos 0.25 segundos.
O campo de temperaturas no conjunto peça/moldação determinado para 2 segundos
após o vazamento é apresentado na figura 3.9 na qual se verifica que apenas no sub-domínio 1 (ver figura 4.4) existe material no estado líquido. Quanto à moldação a
pequena área junto à peça de maior temperatura verificada para 0.25 segundos de
arrefecimento é agora bastante maior.
FI
581.6
556.3
529.1
488.4
386.7
302.5
228.0205.2
201.5
197.3
172.0
161.5
149.3
143.4
139.6
Figura 4.11 – Campo de temperaturas no conjunto peça/moldação decorridos 2 segundos.
Decorridos 4 segundos (figura 4.12), verifica-se que a dimensão da área de material
líquido foi consideravelmente reduzida ao passo que o nível de temperatura no sub-
88
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 105/134
Capítulo IV – Simulação numérica da solidificação em fundição
domínio 4 (ver figura 4.4) desceu cerca de 40ºC em relação aos valores de temperatura
verificados para 2 segundos de arrefecimento.
FI578.2
574.8
553.0
496.5
475.6
429.0
329.4
273.5
217.6
199.2
189.6
161.7
153.0
148.2
139.7
Figura 4.12 – Campo de temperaturas no conjunto peça/moldação decorridos 4 segundos.
O fim da solidificação acontece quando decorridos 6 segundos de arrefecimento, isto é,
toda a peça se encontra abaixo da temperatura de solidus (575ºC), situação que se ilustra
na figura 4.13. Refira-se o facto do gradiente de temperatura entre a peça fundida e as
imediações da moldação ser bastante reduzido e, como era de esperar, as zonas mais
frias da moldação serem as que se encontram mais afastadas da peça e simultaneamente
mais próximas do meio ambiente.
FI
574.5
573.9
566.4
547.2
519.9
465.3
410.7
301.4
246.8
219.2
197.4
164.9
150.2
145.1
139.9
Figura 4.13 – Campo de temperaturas no conjunto peça/moldação decorridos 6 segundos.
89
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 106/134
Capítulo IV – Simulação numérica da solidificação em fundição
Figura 4.14 – Campo de temperaturas na peça no fim da solidificação.
Figura 4.15 – Campo de temperaturas na moldação no fim da solidificação.
FI
547.2
534.3
437.6407.9
383.3
328.7
260.9
225.7
197.0
176.7
164.9
157.4
151.3
147.5
140.8
139.6
FI
574.5
573.8
570.0
558.9
543.2
527.6
511.9
474.5
465.0
433.7
424.7402.4
371.1
355.5
339.8
Nas figuras 4.14 e 4.15 apresenta-se em pormenor o campo de temperaturas
determinado quando finda a mudança de fase para a peça e moldação, respectivamente.
As zonas da peça que sofrem arrefecimento mais acentuado correspondem aos
sub-domínios 2 e 5 (ver figura 4.5) que são precisamente as zonas de menor espessura
da peça. Em contraste, a zona de arrefecimento mais lento corresponde ao sub-domínio
1 (ver figura 4.5), visto tratar-se da zona de maior densidade mássica da peça.
No que concerne à moldação, como se pode constatar as temperaturas mais elevadas
ocorrem, como era previsível, junto da peça. Considerando o contorno externo da
moldação verifica-se que a zona de temperatura mais elevada acontece no sub-domínio
12 (ver figura 4.5 e 4.15) que encontra justificação na sua menor espessura.
Fazendo uma análise comparativa com o campo de temperaturas inicial (figura 4.8),
verifica-se que as zonas mais frias continuam a ser os cantos superiores da moldação, no
90
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 107/134
Capítulo IV – Simulação numérica da solidificação em fundição
entanto, o canto posterior esquerdo, por se encontrar nas imediações da zona mais
quente da peça, apresenta um nível de temperatura bastante superior ao inicial.
VII.3 – Resultados considerando transferência de calor Newtoniana parcial
Trata-se da segunda abordagem do problema onde o coeficiente de proporcionalidade
entre o fluxo de calor e a diferença de temperaturas na interface peça/moldação (relação
de condutividades térmicas dos materiais) da primeira abordagem é substituído por um
coeficiente de transferência de calor Newtoniano. Todas as restantes condições foram
mantidas.
Iniciando, então, a análise dos resultados obtidos nesta abordagem pela figura 4.16,
onde se mostra o campo de temperaturas no conjunto peça/moldação decorridos apenas
0,25 segundos de arrefecimento.
FI
585.0
583.9
562.2
539.5
516.7
493.9
471.2
323.9
277.7
270.0
264.0
255.4
246.7
235.9
229.3
Figura 4.16 – Campo de temperaturas no conjunto peça/moldação decorrido 0.25 segundos.
Verifica-se que a peça se encontra muito próximo da temperatura de liquidus, existindo
uma reduzida parcela, nas zonas mais junto da moldação, onde a temperatura desceu
ligeiramente.
Na moldação pode ser observada a ocorrência de uma área de temperatura mais elevada
na zona de contacto com a peça.
Na figura 4.17 está representado o campo de temperaturas resultante para o conjunto
peça/moldação decorridos 4 segundos após o vazamento. Pode verificar-se que:
91
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 108/134
Capítulo IV – Simulação numérica da solidificação em fundição
o Na zona correspondente ao canal de enchimento, domínio D1 (ver figura
4.4), ocorrem as temperaturas mais elevadas, como previsível dada a sua
maior massividade;
o Uma parcela considerável da peça está já no estado sólido (abaixo de
575°C), com excepção para as zonas de maior espessura dos sub-domínios
1, 3 e 5 (ver figura 4.5);
FI
584.4
577.4
572.7
539.3
516.7
494.2
418.9
359.3
291.1
268.5
253.6
246.0
237.5
229.6
224.0
Figura 4.17 – Campo de temperaturas no conjunto peça/moldação decorridos 4 segundos.
O campo de temperaturas do conjunto peça/moldação obtido 8 segundos após o
vazamento representado na figura 4.18, mostra que:
o Apenas na zona aproximadamente cilíndrica da peça existe material no estado
líquido.
o Na zona correspondente ao sub-domínio 4 forma-se uma estricção das
isotérmicas, reflectindo a influência da existência de maior espessura nas
extremidades em que se encontra com os sub-domínios 3 e 5;
92
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 109/134
Capítulo IV – Simulação numérica da solidificação em fundição
FI
581.9
575.0
559.7
537.4
492.9
469.6
449.9436.3
359.2
314.7
292.4
270.1
254.6
236.0
226.4
Figura 4.18 – Campo de temperaturas no conjunto peça/moldação decorridos 8 segundos.
Na figura 4.19 representa-se o campo de temperaturas determinado decorridos 12
segundos após o vazamento. Observa-se que apenas uma reduzida parcela da peça(sub-domínio 1) se encontra acima da temperatura de solidus.
Verifica-se, também que a zona de temperatura mais elevada na moldação corresponde
ao sub-domínio 12 devido à sua menor espessura e contactar directamente com a zona
de temperatura mais elevada da peça (sub-domínio 1).
FI
577.9
575.0
571.7
561.1
535.0
513.4
415.0
383.4
375.3
341.4
319.4
297.8
276.2
254.7
236.1
Figura 4.19 – Campo de temperaturas no conjunto peça/moldação decorridos 12 segundos
O fim da solidificação da peça ocorre no sub-domínio 1, após sensivelmente 14,8
segundos de arrefecimento, como se mostra na figura 4.20. De salientar os elevados
gradientes de temperatura que se verificam nos sub-domínios 1 e 2 como facilmente se
pode constatar pela figura 4.21. No sub-domínio 1, o elevado gradiente de temperatura
encontra justificação na sua grande massividade. No sub-domínio 2 já não é tão
evidente. Tendo em consideração a sua reduzida espessura a temperatura baixará mais
rápido, o que se verifica na zona central onde a influência do sub-domínio 1 (zona mais
93
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 110/134
Capítulo IV – Simulação numérica da solidificação em fundição
quente da peça) e do sub-domínio 3 que se encontra com temperaturas mais altas que a
zona central do sub-domínio 2, embora mais baixa que as verificadas no sub-domínio 1,
não é tão significativa.
FI
575.0
573.4
568.6
552.8
531.7
510.0
444.9
401.5
359.7
336.4
314.7
293.0
271.3
249.6
230.0
Figura 4.20 – Campo de temperaturas no conjunto peça/moldação decorridos 14.8 segundos.
Para finalizar uma palavra para a zona mais fria da peça. Como era previsível ocorre na
extremidade da peça (figura 4.21) visto tratar-se de uma zona de reduzida espessura e
três das suas fronteiras contactarem com a moldação.
FI
575.0
572.3
568.9
561.1
547.2
533.2
505.3
477.4
440.1
394.3
371.1
363.5
360.6
358.7
354.5
Figura 4.21 – Campo de temperaturas da peça no fim da solidificação.
Na figura 4.22 apresenta-se o campo de temperaturas obtido no fim da solidificação
para a moldação. Como se observa as temperaturas mais elevadas ocorrem, como era
previsível, junto da peça. Considerando o contorno externo da moldação verifica-se que
a zona de temperatura mais elevada acontece no sub-domínio 12 (ver figura 4.5) o que
encontra justificação o facto de contactar directamente com a zona de temperatura mais
elevada da peça (sub-domínio 1) e apresentar uma reduzida espessura.
94
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 111/134
Capítulo IV – Simulação numérica da solidificação em fundição
FI
457.8
423.5
394.0
377.9
360.5
342.5
327.2312.9
294.6
277.9
268.9
261.2
244.5
238.4
230.0
Figura 4.22 – Campo de temperaturas na moldação no fim da solidificação da peça.
Fazendo uma análise comparativa com o campo de temperaturas inicial (figura 4.8),
verifica-se que as zonas mais frias continuam a ser os cantos superiores da moldação, no
entanto, o canto posterior esquerdo, por se encontrar nas imediações da zona mais
quente da peça, apresenta um nível de temperatura bastante superior ao inicial.
VII.4. Resultados considerando transferência de calor Newtoniana
Trata-se da terceira e última abordagem ao problema da solidificação, onde para além
do coeficiente de transferência de calor Newtoniano aplicado à interface D1/D2
(peça/moldação) é também aplicado o mesmo tipo de coeficiente às interfaces D2/D3 e
D3/D4 (ver figura 4.4), embora de valor distinto e de acordo com o quadro 4.3.
FI
585.0
583.9
539.5
516.7
493.9
471.2
448.4
301.9293.6
271.8
261.9
254.8
246.7
237.0
223.0
Figura 4.23 – campo de temperaturas no conjunto peça/moldação ao fim de 0.25 segundos.
95
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 112/134
Capítulo IV – Simulação numérica da solidificação em fundição
FI
584.4
576.8
572.1
538.7
515.8
438.0
332.8287.0
273.7
267.1
256.1
248.0
241.2
230.3
218.8
Figura 4.24 – campo de temperaturas no conjunto peça/moldação ao fim de 4 segundos.
FI
582.0
577.5514.1
491.5
473.1
444.0
441.6
403.1
355.8
310.5
287.9
265.3
242.3
237.3
220.5
Figura 4.25 – campo de temperaturas no conjunto peça/moldação ao fim de 8 segundos.
FI
578.0
576.2
566.6
534.3
490.5
424.7
382.3
356.9
337.1
318.4
280.7
270.6
263.1
239.4
230.2
96
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 113/134
Capítulo IV – Simulação numérica da solidificação em fundição
Figura 4.26 – campo de temperaturas no conjunto peça/moldação ao fim de 12 segundos.
FI
575.0
572.3
565.2
531.6510.3
435.9
386.9
366.5
334.2
317.9
282.9
270.7
264.5
244.0
230.8
FI575.0
574.3
569.5
561.4
547.7
520.4
490.7
465.9
432.1
397.6
375.7
370.3
368.1
365.3
361.5
Figura 4.27 – campo de temperaturas no conjunto peça/moldação ao fim de 14.9 segundos.
Figura 4.28 – campo de temperaturas na peça no fim da solidificação.
FI
440.3
401.1
373.7
362.0
333.2
318.1
298.1
281.9
275.0
268.6
261.2
249.1
240.1
234.4
230.8
Figura 4.29 – campo de temperaturas na moldação no fim da solidificação.
Como se pode constatar pela sequência de figuras 4.23 a 4.29 os campos de
temperaturas obtidos são em quase tudo semelhantes aos apresentados na secção
anterior (sequência de figuras 4.16 a 4.22) pelo que as considerações lá produzidas
também se aplicam a esta subsecção, salvo as seguintes excepções:
97
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 114/134
Capítulo IV – Simulação numérica da solidificação em fundição
o Descontinuidade no valor da temperatura nas interfaces dos domínios D2/D3 e
D3/D4 (ver figura 4.4), onde o contacto perfeito foi substituído por um fluxo de
calor proporcional à diferença de temperaturas que se verificam nas referidas
interfaces.
o Os cantos posteriores da moldação, correspondentes aos sub-domínios 16 e 17
(ver figura 4.5), apresentam níveis de temperatura consideravelmente mais
baixos.
o A zona correspondente ao domínio D3 apresenta valores de temperatura mais
elevados, assim como toda a zona envolvente da peça.
A apresentação de resultados conhece o seu epílogo na figura 4.30, onde se representa a
variação da fracção sólida ao longo do tempo nos diferentes sub-domínios que
constituem a peça.
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
0.0 1.3 2.5 3.8 5.0 6.3 7.5 8.8 10.0 11.3 12.5 13.8 15.0 16.3 17.5
Tempo [s]
F r a c ç ã o
S ó l i d a
Fs 5
Fs 4
Fs 3
Fs 2
Fs 1
Figura 4.30 – Variação da fracção sólida em função do tempo para os diferentes sub-domínios queconstituem a peça.
A fracção sólida assume o valor unitário assim que toda a massa líquida tenha sido
solidificada. Tal como se pode constatar pela figura 30, é no sub-domínio 1 ( fs1) que
98
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 115/134
Capítulo IV – Simulação numérica da solidificação em fundição
esse valor é atingido mais tarde, o que vem de encontro com os campos de temperatura
anteriormente apresentados.
Os restantes sub-domínios apresentam variações da fracção sólida semelhantes.Destaque-se o sub-domínio 4, o primeiro a estar completamente solidificado e o
sub-domínio 2, que apesar da sua menor massividade é o terceiro a findar a
solidificação. Esta aparente contradição, encontra explicação no facto desse
sub-domínio contactar com o sub-domínio 1, onde a fracção sólida correspondente no
fim da solidificação do sub-domínio 2 é de, sensivelmente, 20 %.
VIII. Conclusão
Em primeiro lugar importa referir o facto do fim da solidificação nas situações
apresentadas em VI.3 e VI.4 ser muito semelhante (14.8 e 14.9 segundos,
respectivamente) o leva a concluir que os fluxos de calor introduzidos nas interfaces
D2/D3 e D4/D3 não influenciam significativamente o arrefecimento da peça. Essa
influência é sentida mais intensamente nas zonas da moldação mais próximas da
interface com a peça, onde a temperatura assume valores mais elevados.
Em segundo lugar deve-se mencionar a discrepância em termos de tempo de
solidificação que existe entre as duas situações anteriormente referidas (14.8 segundos),
e a situação apresentada em VI.2 que é de apenas 6 segundos. Tal diferença indica que o
escoamento de calor se processa de um modo muito mais rápido para esta última
situação, onde, recorde-se, foi assumido que a transferência de calor na interface
peça/moldação se processa apenas por difusão.
Estas discrepâncias realçam o papel preponderante da resistência térmica das diferentes
interfaces na evolução dos campos de temperaturas e levantam uma questão pertinente.
Qual das três abordagens efectuadas será a que melhor descreve o arrefecimento do
conjunto peça/moldação tal como acontece na realidade?
A resposta a esta questão será dada no próximo capítulo, onde se fará a comparação
com dados obtidos experimentalmente e que podem ser consultados na referência [1].
99
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 116/134
Capítulo IV – Simulação numérica da solidificação em fundição
O cálculo da fracção sólida nos diferentes sub-domínios que constituem a peça permitiu
confirmar os resultados verificados para o tempo e localização do fim da solidificação.
100
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 117/134
CAPÍTULO V
Validação do modelo numérico
I. Introdução
Os principais factores que limitam a aplicação dos modelos analíticos em fundição, são a
geometria do domínio (geralmente complexa e portanto de descrição analítica difícil), e as
equações que regem os fenómenos físicos envolvidos no processo. As dificuldades de
modelação analítica tornam necessário e estimulante o recurso aos modelos numéricos. A
validação destes deverá ser feita por comparação entre a solução numérica obtida e a
experimentação física realizada em condições semelhantes às definidas para o modelo
numérico.
O método experimental continua a ter grande importância, principalmente quando os
fenómenos envolvidos são de grande complexidade, mas a tendência para o uso dos
métodos numéricos tem vindo a aumentar devido à fiabilidade dos resultados e custos
reduzidos quando comparados com os métodos experimentais, de modo que a concepção
dos modelos experimentais destinados a testar os modelos numéricos deve ter como
objectivo contemplar situações para as quais não existe solução analítica.
Com o objectivo de testar a validade da modelação numérica proposta neste trabalho, os
dados numéricos, utilizando o método dos volumes finitos, serão comparados com dados
experimentais e, também, com dados numéricos, utilizando o método das diferenças finitas,
obtidos na referência [1].
101
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 118/134
Capítulo V – Validação do modelo numérico
II. Comparação entre dados experimentais e numéricos
Antes de se dar início à confrontação de dados numéricos e experimentais é importante
referir que tal comparação será efectuada apenas para as coordenadas onde os termoparesforam colocados (figura 5.1). Os termopares utilizados são do tipo K (constituídos pelo par
Ni 10%Cr-Ni 5%Al) com um diâmetro total de 1 mm.
A confrontação entre dados experimentais e dados numéricos é feita em duas etapas. Na
primeira etapa, figuras 5.2 e 5.3, comparam-se os dados experimentais com valores
numéricos obtidos pelo código computacional desenvolvido neste trabalho para o caso de
transferência de calor por condução na interface peça/moldação (FV-K) conjuntamente com
a transferência de calor Newtoniana na mesma interface (FV-h), situações que constituíram
os dois primeiros casos de estudo no capítulo anterior. Numa segunda etapa far-se-á o
mesmo tipo de comparação para o terceiro caso de estudo do anterior capítulo (FV), figuras
5.4 e 5.5, confrontando, ainda, com resultados numéricos, utilizando o método das
diferenças finitas generalizadas, obtidos na referência [1].
T8 T12 T13 T11 T9 T10
T1 T2 T3 T4 T5 T6 T7
Figura 5.1 - Posição dos termopares utilizados na medição de temperaturas na secção transversal do
conjunto peça/moldação estudada.
Efectuando uma análise à figura 5.2 referente à parte inferior da moldação, verifica-se que
os valores numéricos FV-h concordam bastante bem com os valores experimentais.
Verifica-se, também, que os valores numéricos FV-K, embora apresentem o mesmo padrão
de comportamento ao longo do tempo que os valores de temperatura resultantes da medição
experimental não apresentam uma discrepância de valores considerável.
102
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 119/134
Capítulo V – Validação do modelo numérico
T1
0
100
200
300
400
500
0 5 10 15 20
Tempo [s]
T e m p e r a
t u r a [ º C ]
Exp FV-K FV-h
T2
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T3
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T4
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T5
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T6
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T7
10 15 20
m po [s]
0
100
200
300
400
500
0 5
Te
FV-K FV-h
T1 T2 T3
Exp
T4 T5 T7 T6
Figura 5.2 – Comparação entre os resultados da medição experimental [1] e do cálculo numérico por volume
finito (FV-K e FV-h) para a meia moldação inferior.
103
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 120/134
Capítulo V – Validação do modelo numérico
T8
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T9
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T10
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T11
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T12
0
100
200
300
400
500
0 5 10 15 20
Te m po[s]
Exp FV-K FV-h
T13
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FV-K FV-h
T8 T9 T10 T11 T12 T13
Figura 5.3 – Comparação entre os resultados da medição experimental [1] e do cálculo numérico
por volume finito (FV-K e FV-h) para a meia moldação superior.
104
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 121/134
Capítulo V – Validação do modelo numérico
A comparação entre os resultados da medição experimental e do cálculo numérico
utilizando o método dos volumes finitos com condições de fronteira difusivas (FV-K) e
Newtonianas (FV-h) na interface peça/moldação para a meia moldação superior é
apresentada na Figura 5.3.
Efectuando a sua análise, constata-se que:
o os valores de temperatura obtidos por simulação numérica acompanham o padrão
de variação dos valores de temperatura verificados na experimentação;
o a concordância entre os valores numéricos e os valores da medição experimental
não é tão satisfatória como se verificou para a parte inferior da moldação;
o a discrepância entre as diferentes abordagens numéricas já não é tão significativa
como a que se verificou para a parte inferior da moldação.
Atente-se, agora, nos gráficos referentes aos termopares T8 e T11, únicas localizações
onde os resultados numéricos FV-K são melhores que os seus homólogos FV-h.
No que concerne ao termopar T8, como se situa no sub-domínio 17 tem como interfaces os
sub-domínios 1, 11 e 12 (ver figura 4.5). Na primeira abordagem do problema (FV-K)
estipulou-se transferência de calor por difusão na interface com o sub-domínio 1 e contacto
perfeito nas restantes. Sabe-se que desta abordagem ao problema resultou o menor tempo
de solidificação das três abordagens realizadas, o que permitiu concluir que a transferência
de calor ocorre a uma taxa mais elevada nessa abordagem. Na segunda abordagem foi
estipulado transferência de calor Newtoniana na interface entre os sub-domínios 1e 17 e o
erro médio aumentou em 5%. Na terceira abordagem este tipo de transferência de calor foi
estendido a mais uma interface (sub-domínios 11 e 17) e o erro médio foi aumentado em
2%. Sabendo que FV-K apresenta melhor concordância com os valores experimentais, pode-se concluir que o valor do coeficiente de transferência de calor Newtoniano nas
referidas interfaces não é o mais adequado.
105
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 122/134
Capítulo V – Validação do modelo numérico
Relativamente ao termopar T11, atendendo a que os valores de temperatura da primeira
abordagem (FV-K), onde o fim da solidificação ocorria mais cedo que na segunda
abordagem (FV-h), se aproxima mais das medições experimentais e, tendo em atenção a
particularidade de estar localizado sobre a interface D3/D4 (ver figura 4.4), pode-se afirmar
que o calor é de certa forma retido no domínio D3 que como já foi referido corresponde ao
inserto da moldação para formar a reentrância apresentada pela peça.
A explicação para este facto, pode estar na definição das condições de transferência de calor
na interface D3/D4, onde na realidade não existe um contacto perfeito como até aqui foi
considerado, mas sim condições que promovem a acumulação do calor naquela zona da
moldação (domínio D3). Na realidade existe uma interface onde a transferência de calor é
Newtoniana.
Tendo em atenção os valores iniciais de temperatura para os pontos de medição referentes
aos termopares T8, T12 e T13, verifica-se que estes são consideravelmente mais elevados.
Como, em adição, o padrão de comportamento dos valores numéricos de temperatura
acompanha as medições experimentais, pode-se concluir que a referida discrepância é
devida ao valor inicial de temperatura que como se sabe foi arbitrado.
No sentido de comparar dados numéricos obtidos pelo código desenvolvido utilizando o
método dos volumes finitos com dados, também numéricos utilizando o método das
diferenças finitas obtidos na referência [1], deu origem à terceira e última abordagem do
problema da solidificação que apresenta, como já se teve oportunidade de referir, algumas
interfaces onde as condições de transferência de calor são distintas.
106
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 123/134
Capítulo V – Validação do modelo numérico
T1
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T2
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T3
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T4
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T5
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T6
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T7
0
100
200
300
400
500
0 5
Te
10 15 20
m po[s]
FD FV
T4 T5 T1 T2 T3 T7 T6
Exp
Figura 5.4 – Comparação entre os resultados da medição (experimental) e do cálculo numérico por
diferença finita [1] e volume finito (implementado neste trabalho) para a meia moldação inferior.
107
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 124/134
Capítulo V – Validação do modelo numérico
T8
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T9
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T10
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T11
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T12
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T13
0
100
200
300
400
500
0 5 10 15 20
Tem po [s]
Exp FD FV
T8 T9 T10 T11 T12 T13
Figura 5.5 – Comparação entre os resultados da medição (experimental) e do cálculo numérico por
diferença finita [1] e volume finito (implementado neste trabalho) para a meia moldação superior.
108
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 125/134
Capítulo V – Validação do modelo numérico
Uma das interfaces que sofre alterações é a interface D3/D4, onde o contacto perfeito da
segunda abordagem é substituído por um fluxo de calor Newtoniano, situação já aludida
anteriormente. A interface D2/D3 é outra em que é feito o mesmo tipo de alteração, embora
como se teve oportunidade de verificar o contacto perfeito considerado na segunda
abordagem representa satisfatoriamente o que acontece na realidade. Tal alteração tem em
vista criar as mesmas condições de simulação que as usadas por Monteiro [1] utilizando o
método das diferenças finitas e, assim, confrontar dados numéricos obtidos por métodos de
discretização distintos.
Nas figuras 5.3 e 5.4 apresenta-se a referida confrontação, sendo apresentados,
adicionalmente, os valores de temperatura relativos à experimentação para melhor
comparação.
Verifica-se pelas figuras 5.3 e 5.4 que as curvas resultantes da modelação numérica
acompanham os resultados experimentais para lá do fim da solidificação da peça. Note-se
que o máximo nos resultados experimentais, determinado no termopar T4, que contacta
com o metal vazado por estar colocado à face da moldação, ocorre antes desse instante,
como consequência da menor espessura local da peça.
No sentido de quantificar as discrepâncias nos valores de temperatura das diferentes
abordagens, apresenta-se no quadro 5.1 os erros percentuais médios para cada termopar em
relação aos valores que se obtêm pela experimentação. Efectuando a sua análise verifica-se
que a terceira abordagem do problema apresenta globalmente menor erro médio.
Realce-se o facto do menor erro médio, relativamente aos termopares T8 e T9, se obter na
primeira abordagem do problema, cujos resultados não são globalmente satisfatórios e o
facto do maior valor do erro se verificar para o termopar T13 na segunda e terceiraabordagens e ser, também, um dos mais elevados na primeira abordagem.
109
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 126/134
Capítulo V – Validação do modelo numérico
Quadro 5.1 – Erro médio percentual das diferentes abordagens ao problema da solidificação
Termopar
1ª Abordagem
(FV-k)Erro [%]
2ª Abordagem
(FV-h)Erro [%]
3ª Abordagem
(FV)Erro [%]
T1 17% 2% 3%
T2 13% 4% 3%
T3 19% 3% 3%
T4 21% 4% 3%
T5 25% 3% 3%
T6 26% 3% 4%
T7 10% 2% 5%
T8 6% 10% 13%
T9 6% 3% 4%
T10 11% 5% 1%
T11 7% 15% 3%
T12 18% 14% 10%
T13 21% 18% 13%
Verifica-se, também, que o menor valor do erro médio relativo ao termopar T7 (que se
encontra sobre a interface que une os sub-domínios 8 e 16) se obtém na segunda abordagem
do problema, no qual, como se sabe, foi considerado a existência de um contacto perfeito
entre as partes da moldação. Tal comportamento permite concluir que as superfícies em
contacto apresentam baixa rugosidade, dado que a transferência de calor será tanto mais
lenta quanto maior for a quantidade de ar existente entre as partes em contacto.
Para finalizar esta análise comparativa de resultados, uma palavra para a confrontação entre
métodos de discretização (diferenças finitas e volumes finitos). A análise das figuras 5.4 e
5.5 permite verificar o comportamento similar entre os diferentes métodos com ligeira
vantagem para o método dos volumes finitos implementado neste trabalho.
110
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 127/134
Capítulo V – Validação do modelo numérico
III. Conclusão
A comparação entre valores de temperatura numéricos e experimentais, permite concluir
que a evolução das curvas obtidas numericamente acompanha o padrão de comportamento
das medições experimentais.
A análise do quadro 5.1 onde se apresentam os erros médios para as três abordagens ao
fenómeno da solidificação permite concluir que os resultados numéricos obtidos na terceira
abordagem apresentam globalmente um erro médio inferior, pelo que é a modelação que
melhor representa a realidade.
Atendendo à evolução positiva na concordância entre os resultados numéricos e
experimentais evidenciados nas figuras 5.2, 5.3, 5.4 e 5.5, pode afirmar-se que é a
resistência térmica nas interfaces que rege a evolução dos campos de temperaturas nas
moldações metálicas.
A interface D3/D4 cria uma solução de continuidade entre os domínios D3 e D4, e é fonte
de complexidade adicional da modelação, uma vez que são diminutos os dados conhecidos
sobre o comportamento do coeficiente de transferência de calor nestas interfaces, o mesmo
se passando nas interfaces entre os domínios D2 e D3.
111
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 128/134
CAPÍTULO VI
Conclusões e perspectivas de trabalhos futuros
I. Conclusões
A resolução analítica de uma equação diferencial que rege determinado fenómeno físico
só é possível em casos particulares. Daqui decorre a necessidade de reformulação do
problema, dando-lhe uma forma puramente algébrica.
Para se recorrer aos métodos numéricos é também necessário representar adequadamente
a geometria do domínio de cálculo, o que se consegue fazendo a escolha de um número
finito de nós computacionais, supostos representativos de todo o domínio, onde irá ser
determinada a solução numérica.
Quando os domínios, físico e computacional, forem idênticos, não é necessário qualquer
transformação, mas no caso de geometrias físicas complexas, a transformação pode
provocar modificações drásticas quer nas equações constitutivas do fenómeno físico quer
nas condições de fronteira associadas.
A dificuldade de descrição dos contornos das formas geométricas complexas pode ser
ultrapassada dividindo a geometria complexa em troços de geometria mais simples,estabelecendo malhas curvilíneas em cada um desses troços, em que serão
adequadamente definidas as interacções entre os mesmos.
A utilização do método dos volumes finitos com formulação em coordenadas
cartesianas foi testada por comparação com uma solução analítica conhecida para a
condução de calor num corpo bidimensional utilizando e revelou ter um comportamento
bastante satisfatório. Nesta aplicação testaram-se, ainda, várias formulações de
112
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 129/134
Capítulo VI – Conclusões e perspectivas de trabalhos futuros
discretização temporal: Euler implícita, Crank-Nicolson e Ritchmayer-Morton, sendo
esta última a que mais se aproxima da solução analítica.
A formulação do método dos volumes finitos utilizando coordenadas curvilíneas foitestada por comparação com uma solução analítica conhecida para o problema da
condução de calor em domínio circular tendo apresentado resultados bastante
satisfatórios. Nesta aplicação foi, ainda, testada a utilização de diferentes métodos de
resolução de sistemas de equações resultantes de processos de discretização (TDMA,
Gauss-Seidel e SIP) na qual se verificou que o método iterativo SIP ( Strongly Implicit
Procedure) permite maior precisão e contribui para a economia do tempo de
computação.
A aplicação do código desenvolvido utilizando coordenadas curvilíneas a uma
geometria complexa concreta foi realizada em três etapas sendo alteradas as condições
de transferência de calor nas interfaces dos diferentes domínios em cada uma delas.
A primeira abordagem ao problema da solidificação, onde se considerou que a
transferência de calor da peça para a moldação se processa apenas por condução,
revelou-se, como esperado, inapropriada para representar o que acontece na realidade.
Na segunda abordagem ao problema da solidificação foi considerado que a transferência
de calor na interface peça/moldação é Newtoniana, assunção bastante mais realista dado
a concordância satisfatória com dados das medições experimentais obtidos na referência
[1].
A terceira e última abordagem, onde se considera transferência de calor Newtoniana,
também nas interfaces D3/D4 e D2/D3 teve como objectivo criar as mesmas condições
de simulação que as utilizadas por Monteiro [1] utilizando o método das diferenças
finitas, e assim, comparar com a solução que o método dos volumes finitos proporciona.
A análise comparativa de resultados entre métodos de discretização (diferenças finitas e
volumes finitos) permitiu verificar o comportamento similar entre os diferentes métodos
com ligeira vantagem para o método dos volumes finitos implementado neste trabalho.
A análise do quadro 5.1 onde se apresentam os erros médios para as três abordagens ao
fenómeno da solidificação permite concluir que os resultados numéricos obtidos na
113
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 130/134
Capítulo VI – Conclusões e perspectivas de trabalhos futuros
terceira abordagem apresentam globalmente um erro médio inferior, pelo que é a
modelação que melhor representa a realidade.
Atendendo à evolução positiva na concordância entre os resultados numéricos eexperimentais evidenciados nas figuras 5.2 5.3 5.4 e 5.5, pode afirmar-se que é a
resistência térmica nas interfaces que rege a evolução dos campos de temperaturas nas
moldações metálicas.
A interface D3/D4 cria uma solução de continuidade entre os domínios D3 e D4, e é
fonte de complexidade adicional da modelação, uma vez que são diminutos os dados
conhecidos sobre o comportamento do coeficiente de transferência de calor nestas
interfaces, o mesmo se passando nas interfaces entre os domínios D2 e D3.
Uma vez que a distribuição global de temperatura calculada está de acordo com as
medições experimentais efectuadas por Monteiro [1], pode-se afirmar que são válidas as
informações sobre a solidificação fornecidas pelo modelo.
O método dos volumes finitos utilizando coordenadas curvilíneas permite assim
resolver o problema da solidificação sem limitações de carácter geométrico o que
constitui a maior vantagem do método.
A utilização de malhas curvilíneas permite o desenvolvimento de algoritmos mais
eficientes e com melhor capacidade de descrever os contornos das peças através de um
menor número de elementos, sendo a especificação do problema independente da forma
dos contornos do domínio físico, em virtude do domínio computacional resultante ser
fixo.
A divisão de uma peça em troços de geometria mais simples, resulta bem no métododos volumes finitos, porque é possível estabelecer de forma adequada os balanços nas
interfaces virtuais definidas, contornando a sua complexidade geométrica.
A utilização de coordenadas curvilíneas para definir os nós da malha produz uma
estrutura organizada que permite que todo o esforço de cálculo seja efectuado num
rectângulo computacional, desde que as equações em causa tenham sido transformadas
para que as coordenadas curvilíneas substituam as coordenadas cartesianas como
variáveis independentes na formulação do problema.
114
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 131/134
Capítulo VI – Conclusões e perspectivas de trabalhos futuros
O programa desenvolvido pode fazer-se evoluir rapidamente, por incorporação de novas
rotinas ou rotinas já existentes para o método dos volumes finitos, novos conhecimentos
sobre os materiais utilizados, ou dados sobre novos materiais.
II. Perspectivas de desenvolvimentos futuros
O primeiro desenvolvimento que cabe referir é naturalmente a alteração do código para
fazer a modelação tridimensional, o que pode ser feito a partir das relações de
transformação apropriadas.
Uma segunda vertente passa por melhorar a caracterização física das ligas metálicas, nosentido de descrever mais adequadamente as condições de arrefecimento reais e ter,
também, em conta a evolução das características térmicas dos materiais com a
temperatura.
Outra aplicação consiste em simular a utilização de componentes efémeros na execução
de troços da moldação, muito frequentemente utilizados em fundição em moldação
metálica.
É também possível modelar o aparecimento e evolução da folga entre a peça e a
moldação, para o que será necessário determinar a deformação das moldações e da peça
sob acção dos campos de temperaturas.
Simular a utilização de pinturas protectoras, isolantes ou desmoldantes, que permitem
conferir diferentes qualidades de acabamento superficial à parede da moldação, e fazer
variar significativamente as condições de transferência de calor na interface
metal/moldação, é outro possível desenvolvimento.
Uma outra linha de investigação poderá ser a modelação e simulação do enchimento da
moldação, o que será a combinação do problema da solidificação com um problema de
mecânica dos fluídos, e que passará necessariamente pela utilização das equações de
Navier-Stokes.
115
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 132/134
Referências
[1] Monteiro, António A.C., 1996. “Estudo do Comportamento Térmico de Moldações
Metálicas para Fundição Aplicando o Método das Diferenças Finitas Generalizadas”,
Tese de Doutoramento, UM, Braga, Portugal.
[2] M.T.Clavaguera-Mora, C.Comas, J.L.Touron, M.Garcia, O. Guixà, N. Clavaguera.
”Solidification Study of Cu-Based Alloys Obtained by Continuous Casting”. Journal of
Materials Science 34 (1999) 4347-4350.
[3] Minkowycz, W.J., Sparrow, E.M., Schneider, G.E., Pletcher, R.H., “Handbook of
Numerical Heat Transfer”. Interscience. 1988.
[4] Kothe, D., Juric, D.; Lam, K. and Lally, B.: “Numerical recipes for mould filling
simulation”. Los Alamos National Laboratory, NM, 87545, USA.1998.
[5] Monteiro, Eliseu; Monteiro, Monteiro, A.A.C. and A. I. Rouboa. “Heat Transfer
Simulation in the Mould with Generalize Curvilinear Formulation”. ASME PVP student
paper competition, Vol.464. pp. 231-237. 2003.
[6] Monteiro, Eliseu and A. I. Rouboa. “Convective Heat Transfer Simulation in the
Mould with Convective Boundary Conditions”. ASME Summer Heat Transfer
Conference, July 21-23, 2003, Las Vegas, Nevada, USA.
[7] Juric, D.; Tryggvason, G.:“A Front-Tracking Method for Dendritic Solidification”.
Journal of Computational Physics 123, 127-148 (1996).
[8] Voller, V.R.: “Recent developments in the modelling of solidification process”.
Department of Civil and Mineral Engineering, University of Minnesota, Minneapolis,
MN 55455 USA. 1991.
[9] Gonc, I. Z. X.; Zhang, Y. F.; Mujamdar, A.S.: “An Efficient Finite Element
Procedure for the Solution of Enthalpy Model for Phase Change Heat Conduction
Problems“. L.C. Wrobel and C.A.Brebbia ”Computational Modelling of Free and
Moving Boundary Problems”. Vol. 2 Heat Transfer.1991.
[10] Gonc, I. Z. X.; Zhang, Y. F.; Mujamdar, A.S.: “Cyclic Phase Change Heat
Conduction in Thin Composite Slabs“in Computational Modelling of Free and Moving
Boundary Problems (L.C. Wrobel and C.A.Brebbia), Vol. 2 Heat Transfer.1991.
[11] Shamsundar, N.; Sparrow, E.M.: “Analysis of Multidimensional Conduction Phase
Change via the Enthalpy Model, J. Heat Transfer, Trans. ASME, Vol.97, pp. 333-340,
1975.
[12] Comini, G.; Guidice, S.D.; Lewis, R.W.; Zienkiewicz, O.C. “Finite Element
Solution of Non-Linear Heat Conduction Problems with Reference to Phase Change”,
Int. J. Num. Meths. Eng., Vol.8, pp.613-624, 1974.
116
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 133/134
[13] Knoll, D. A.; Kothe, D.B.; Lally, B.: “A New Nonlinear Solution Method for Phase
Change Problems”, Los Numer. Heat Transfer, Part B, 35 (1999), pp.439-459.
[14] Knoll, D. A.; Vanderheyden, W. B.; Mousseau, V. A.; Kothe, D. B.: “On
Preconditioning Newton-Krylov Methods in Solidifying Flow Applications”; Societyfor Industrial and Applied Mathematics, Vol. 23, No. 2, pp. 381-397, 2001.
[15] Sarler, B., A. Alujevic.: “Melting of anisotropic crystals by the Boundary-Domain
integral method”, in Computational Modelling of Free and Moving Boundary Problems
(L.C. Wrobel and C.A.Brebbia), Vol. 2 Heat Transfer.1991.
[16] A.Bermúdez, M.C. Muñiz, P.Quintela,”Numerical Solution of a Stefan Problem
Arising in the Thermoelectrical modelling of Aluminium Electrolytic Cells”, in
Computational Modelling of Free and Moving Boundary Problems (L.C. Wrobel and
C.A.Brebbia), Vol. 2 Heat Transfer.
[17] Carslaw, H.; Jaeger: “Conduction of Heat in Solids”. Clarendon Press, Oxford,1959, p.294.
[18] Ambrosi, D.; Preziosi, L.: “Modelling Injection Moulding Processes with
Deformable Porous Performs”, Society for Industrial and Applied Mathematics, Vol.
61, No. 1, pp. 22-44, 2000.
[19] Ambrosi, D.; Preziosi, L.: “Modelling matrix injection through porous performs,
Composites, Part A, 29 (1998), pp. 5-18.
[20] Fabbri, M.: “An Alternative Approach for the Modelling and Simulation of Binary
Alloy Solidification using a Fixed-Grid Control-Volume Technique”. in Computational Modelling of Free and Moving Boundary Problems (L.C. Wrobel and C.A.Brebbia),
Vol. 2 Heat Transfer.1991.
[21] Norris, Stuart Edward, 2001. “A Parallel Navier Stokes Solver for Natural
convection and Free Surface Flow”. Ph.D Thesis, Faculty of Mechanical Engineering,
University of Sydney. Australia.
[22] Schneider, F.A., Maliska, C.R.,”Uma Formulação em Volumes Finitos Usando
Malhas Não-Estruturadas”. 8th Brazilian Congress of Engineering and Thermal
Sciences - ENCIT, pp. 74, Porto Alegre-RS, Brasil. 2000. [23] Thompson, J. F.; Warsi, Z. U. A.; Mastin, C. W.: "Numerical Grid Generation,Foundations and Applications"; Elsevier Science Publishing Co., Amsterdam, 1985.
[24] Akihiko Nakayama, Satoshi Yokojima and Sankara N.Vengadesen.”Colocated
Finite Difference Method in General Curvilinear Coordinates for Simulation of Flows
Over Arbitrary Geometry”. Memoirs of Construction Engineering Research Institute
Vol.43-A, Nov.2001.
[25] I. Aavatsmark, T. Barkve, O. Boe, T. Mannseth. “ Discretization on Unstructured
Grids for Inhomogeneous, Anisotropic. Part I: Derivation of the Methods “, Society for
Industrial and Applied Mathematics, Vol. 19, No.5, pp.1700-1716, 1998.
117
5/11/2018 Master Thesis 2003 - slidepdf.com
http://slidepdf.com/reader/full/master-thesis-2003 134/134
[26] Versteeg, H. K.; Malalasekera, W.: “An Introduction to Computational Fluid
Dynamics, The finite volume method”. Prentice Hall. 1995.
[27] J.H.Ferziger, M.Peric, “Computational Methods for Fluid Dynamics”. 2nd edition,
Springer Verlag, Berlin, Heidelberg, New York. 1999.
[28] Earl W.Swokowski, “Calculus with Analytic Geometry”. McGraw-Hill.1983.
[29] Heitor Pina, “ Métodos Numéricos”, McGraw-Hill. 1995.
[30] Tannehill, John C.; Anderson, D. A.; Pletcher, Richard H.: ”Computational Fluid
Mechanics and Heat Transfer”, 2nd edition, Taylor & Francis Ltd. 1997.
[31] M.R.Barkhudarov, “Advanced Simulation of the Flow and Heat Transfer Processs
in Simultaneous Engineering”. Flow Science Inc., 1257 40-Th Street, Los Alamos, New
Mexico 87544, USA.2000.
[32] P. N. Hansen, E. Flender, and G. C. Hartman, “MAGMASOFT- The MAGMA
System of Mold filling and Solidification Modeling” in Numerical Simulation of casting
Solidification in Automotive Applications (C. Kim and C., W. Kim, eds.), p. 221. The
Minerals, Metals and Materials Society (1991).
[33] A. M. Hines and J. S. Tu, “Modeling of fluid flow and heat transfer in the
solidification of superalloy investment castings” in Modeling of Casting, Welding and
Advanced Solidification Processes VI (T. S. Piwonka, V. Voller, and L. Katgerman,
eds.), pp. 461–468, The Minerals, Metals and Materials Society (1993).
[34] P. Laty, P. Large, and C. Rigaut, “Simulator, A Software to help the development
of the mold” in Numerical Simulation of casting solidification in automotiveapplications.(C. Kim and C.-W. Kim, eds.). The Minerals, Metals and Materials Society
(1991).
[35] Elliot Stanley Duff. (1999).“Fluid Flow Aspects of Solidification Modeling:
Simulation of Low Pressure Die Casting”. Ph.D. Thesis. Department of Mining &
Metallurgical Engineering. The University of Queensland. Australia.
118
Recommended