Embed Size (px)
Citation preview
UNIVERSIDADE ESTADUAL DE CAMPINAS
Faculdade de Engenharia Mecânica
PEDRO AUGUSTO LANZA DE PAULA
Simulação Multifísica do Processo de
Têmpera Acoplando as Análises Térmica,
Microestrutural e Eletromagnética.
CAMPINAS
2017
PEDRO AUGUSTO LANZA DE PAULA
Simulação Multifísica do Processo de
Têmpera Acoplando as Análises Térmica,
Microestrutural e Eletromagnética.
Dissertação de Mestrado apresentada à Faculdade de
Engenharia Mecânica da Universidade Estadual de
Campinas como parte dos requisitos exigidos para
obtenção do título de Mestre em Engenharia
Mecânica, na Área de Mecânica dos Sólidos e
Projeto Mecânico.
Orientador: Prof. Dr. Renato Pavanello
ESTE EXEMPLAR CORRESPONDE À VERSÃO
FINAL DA DISSERTAÇÃO DEFENDIDA PELO
ALUNO PEDRO AUGUSTO LANZA DE PAULA,
E ORIENTADA PELO PROF. DR. RENATO
PAVANELLO.
CAMPINAS
2017
Agência(s) de fomento e nº(s) de processo(s): Não se aplica.
Ficha catalográfica
Universidade Estadual de Campinas
Biblioteca da Área de Engenharia e Arquitetura
Elizangela Aparecida dos Santos Souza - CRB 8/8098
D44s
de Paula, Pedro Augusto Lanza, 1991-
Simulação multifísica do processo de têmpera acoplando as análises térmica,
microestrutural e eletromagnética. / Pedro Augusto Lanza de Paula. –
Campinas, SP : [s.n.], 2017.
Orientador: Renato Pavanello.
Dissertação (mestrado) – Universidade Estadual de Campinas, Faculdade de
Engenharia Mecânica.
1. Indução magnética. 2. Aquecimento. 3. Aço - Têmpera. 4. Método dos
elementos finitos. 5. Microestrutura - Simulação por computador. 6. Simulação
multifísica. I. Pavanello, Renato, 1959-. II. Universidade Estadual de Campinas.
Faculdade de Engenharia Mecânica. III. Título.
Informações para Biblioteca Digital
Título em outro idioma: Multiphysics simulation of quenching process coupling
electromagnetic, thermal and microstructural analysis.
Palavras-chave em inglês:
Magnetic induction
Heating
heating
Steel - Tempering
Finite element method
Microstructure - Computer simulation
Multiphysics simulation
Área de concentração: Mecânica dos Sólidos e Projeto Mecânico
Titulação: Mestre em Engenharia Mecânica
Banca examinadora:
Renato Pavanello [Orientador]
Pablo Siqueira Meirelles
Walter Jesus Paucar Casas
Data de defesa: 28-08-2017
Programa de Pós-Graduação: Engenharia Mecânica
UNIVERSIDADE ESTADUAL DE CAMPINAS
FACULDADE DE ENGENHARIA MECÂNICA
COMISSÃO DE PÓS-GRADUAÇÃO EM ENGENHARIA MECÂNICA
DEPARTAMENTO DE MECÂNICA COMPUTACIONAL
DISSERTAÇÃO DE MESTRADO ACADÊMICO
Simulação Multifísica do Processo de
Têmpera Acoplando as Análises Térmica,
Microestrutural e Eletromagnética.
Autor: Pedro Augusto Lanza de Paula
Orientador: Prof. Dr. Renato Pavanello
A Banca Examinadora composta pelos membros abaixo aprovou esta Dissertação:
Prof. Dr. Renato Pavanello
DMC/FEM/UNICAMP
Prof. Dr. Pablo Siqueira Meirelles
DMC/FEM/UNICAMP
Prof. Dr. Walter Jesus Paucar Casas
DEMEC/UFRGS
A Ata da defesa com as respectivas assinaturas dos membros encontra-se no processo de
vida acadêmica do aluno.
Campinas, 28 de Agosto de 2017.
Agradecimentos
Este trabalho não poderia ser terminado sem o auxílio de diversas pessoas às quais
presto minha homenagem:
Ao meu orientador, Prof. Dr. Renato Pavanello, pelo auxílio, supervisão e apoio.
Aos funcionários da ThyssenKrupp Metalúrgica Campo Limpo Ltda., em especial ao
Wiliam Tean Su e ao Alex de Souza Rodrigues, pelo auxílio.
À ThyssenKrupp Metalúrgica Campo Limpo Ltda. pelo apoio financeiro.
À minha namorada, Carla, pelo amor, apoio, paciência e compreensão.
À minha família pelo suporte e incentivo durante toda minha formação acadêmica.
Resumo
Componentes sujeitos a carregamentos cíclicos, tais como partes de um motor à
combustão interna, são frequentemente submetidos ao processo de têmpera de modo a
melhorar as propriedades mecânicas e prevenir falhas por fadiga e desgaste em serviço. É
importante que tais componentes possuam uma camada superficial de alta dureza e com
tensões residuais compressivas, aumentando assim a resistência à fadiga e ao desgaste, e um
núcleo tenaz, com alta capacidade em absorver impactos. Neste trabalho, um método de
simulação multifísica do processo de têmpera utilizando o Método dos Elementos Finitos é
proposto. O método de simulação proposto inclui dois estágios: aquecimento e resfriamento.
No primeiro, o componente mecânico, inicialmente a temperatura ambiente, é aquecido por
indução eletromagnética, calculada utilizando as equações de Maxwell para o caso
harmônico, até acima da temperatura de austenitização do aço. No segundo estágio, o
componente é resfriado por imersão, considerando um modelo clássico de convecção e
condução. A microestrutura resultante é calculada usando o modelo de Johnson-Mehl-
Avrami-Kolmogorov e a regra da aditividade de Sheil, para as transformações com caráter
difusional, enquanto a transformação martensítica é calculada utilizando a equação de
Koistinen-Marburguer. A simulação leva em consideração a variação das propriedades
mecânicas em função da temperatura e da microestrutura, enquanto as propriedades
eletromagnéticas são funções da temperatura e da intensidade do campo magnético
(permeabilidade magnética). Como resultado, a distribuição da microestrutura e o perfil de
dureza pós-têmpera são estimados. A implementação é feita em linguagem APDL, usando
como base as rotinas do programa ANSYS. É feita uma análise de influência dos parâmetros
do processo sobre a espessura da camada endurecida e sobre a potência total absorvida. A
metodologia foi aplicada em condições semelhantes às reais de um virabrequim fabricado
pela ThyssenKrupp Metalúrgica Campo Limpo Ltda. O modelo foi validado a partir de
resultados encontrados na literatura para partes do procedimento de simulação. A simulação
integrada do processo, tal como mostrada neste trabalho, não foi encontrada em outras
referências, tratando-se de uma contribuição inovadora.
Palavras Chave: têmpera por indução, JMAK, aquecimento, concentrador de fluxo
magnético, dureza, elementos finitos.
Abstract
Components subjected to cyclic loads, such as parts of internal combustion engine, are
often submitted to quenching process in order to improve its mechanical properties and
prevent the fatigue and wear fail in service. It is important that such components have a high
hardness surface layer with compressive residual stress, thereby increasing the fatigue and
wear resistance, and a tenacious core with high capacity to absorb impacts. In this work, a
method of multiphysics simulation of quenching process using the finite element method is
proposed. The proposed simulation method includes two stages: heating and cooling. At first,
the mechanical component, initially at room temperature, is heated by electromagnetic
induction, calculated using the Maxwell equations for the harmonic case, up to austenitizing
temperature of steel. In the second stage, the component is cooled by immersion, whereas a
classic model convection and conduction. The resulting microstructure is calculated using the
Johnson-Mehl-Avrami-Kolmogorov model and the Sheil’s Additivity Rule, for
transformations with diffusive character, while the martensitic transformation is calculated
using the Koistinen-Marburguer equation. The simulation takes into account the variation of
mechanical properties as a function of temperature and microstructure, while the
electromagnetic properties are functions of temperature and magnetic field strength (magnetic
permeability). As result, the distribution of microstructure and post-quenching hardness
profile are estimated. The implementation is done in APDL language, using as framework the
routines of the ANSYS program. Analysis of the influence of process parameters on the
thickness of the hardened layer and on the total input power is performed. The methodology
was applied approaching real conditions of a crankshaft manufactured by ThyssenKrupp
Metallurgical Campo Limpo Ltda. The model was validated from the results found in the
literature for parts of the simulation procedure. The integrated process simulation, as shown in
this work, was not found in other references, in the case of an innovative contribution.
Keywords: induction hardening, JMAK, heating, magnetic flux concentrator, hardness, finite
element.
Lista de Figuras
Figura 1.1: Geração de um campo eletromagnético através de um indutor em um dado instante
de tempo. .................................................................................................................................. 18
Figura 1.2: Esquema representando as tensões residuais σxx geradas pelo processo de
têmpera. .................................................................................................................................... 19
Figura 2.1: Fases principais da simulação. ............................................................................... 27
Figura 2.2: Elemento finito de 8 nós utilizado na simulação da fase térmica e eletromagnética.
.................................................................................................................................................. 30
Figura 2.3: Densidade do fluxo magnético em tesla em função da temperatura versus a
intensidade do campo magnético em A/m para o aço AISI 1045, adaptado de Li et al (2012).
.................................................................................................................................................. 34
Figura 2.4: Diagrama TTT do aço SAE 1080, adaptado de Woodard et al., (1999). ............... 36
Figura 2.5: Coeficiente de difusão em função da temperatura do aço SAE1080, baseado em
Woodard et al. (1999). .............................................................................................................. 37
Figura 2.6: Expoente n usado na equação de JMAK do aço SAE 1080 baseado em Woodard et
al. (1999). ................................................................................................................................. 38
Figura 2.7: Diagrama TTT esquemático mostrando a conversão da curva de resfriamento em
passos isotérmicos. ................................................................................................................... 38
Figura 3.1: Estratégia de simulação para a fase de aquecimento eletromagnético. .................. 42
Figura 3.2: Estratégia de simulação da fase de resfriamento. .................................................. 46
Figura 3.3: Esquema representando a simulação integrada do processo de têmpera por indução
eletromagnética. ........................................................................................................................ 47
Figura 4.1: Comparação entre a taxa de calor perdido por convecção e por radiação. ............ 51
Figura 4.2: Coeficiente convectivo para um cilindro de 38mm de diâmetro imerso em água.
Extraído de Woodard et al. (1999). .......................................................................................... 52
Figura 4.3: Condições de contorno e localização dos pontos A, B e C da seção axissimétrica
do cilindro utilizado. ................................................................................................................. 53
Figura 4.4: Distribuição microestrutural de perlita e martensita depois de 90s de imersão em
água. .......................................................................................................................................... 53
Figura 4.5: Histórico de temperaturas durante o resfriamento por imersão em água para os
pontos A,B e C. ........................................................................................................................ 54
Figura 4.6: Malha utilizada para a simulação do cilindro. ....................................................... 55
Figura 4.7: Dureza final de acordo com a propriedade utilizada como não sendo dependente
da microestrutura. ..................................................................................................................... 55
Figura 4.8: Densidade Aço SAE 1080...................................................................................... 57
Figura 4.9: Calor Específico do Aço SAE 1080 ....................................................................... 57
Figura 4.10: Condutibilidade Térmica do Aço 1080 ................................................................ 58
Figura 4.11: Dureza final para da simulação com 8,5mm de espessura austenitizada ............. 58
Figura 4.12: Dureza final para da simulação com 5,1mm de espessura austenitizada ............. 59
Figura 4.13: Dureza final para da simulação com 3,6 mm de espessura austenitizada ............ 59
Figura 4.14: Representação do modelo axissimétrico para o caso do cilindro......................... 60
Figura 4.15: Expansão axissimétrica em ¾ do modelo, representando apenas o indutor e o
cilindro. ..................................................................................................................................... 61
Figura 4.16: Campo de temperaturas em Kelvin obtido após 20 segundos de aquecimento. .. 62
Figura 4.17: Perfil microestrutural pós-têmpera para o caso do cilindro. ................................ 64
Figura 4.18: Perfil de dureza pós-dureza para o caso do cilindro. ........................................... 64
Figura 4.19: Mapa de cores da dureza em RC pós-têmpera para o caso do cilindro. .............. 65
Figura 4.20: Objeto de estudo da simulação integrada. ............................................................ 66
Figura 4.21: Malha utilizada para a simulação do moente. ...................................................... 66
Figura 4.22: Expansão axissimétrica em ¾ do modelo, representando o indutor, o
concentrador de fluxo eletromagnético e o moente. ................................................................. 67
Figura 4.23: Condições de contorno e dados de entrada para a simulação eletromagnética. ... 68
Figura 4.24: Condições de contorno para a simulação térmica. ............................................... 69
Figura 4.25: Direções A e B onde a profundidade da camada endurecida foi observada. ....... 69
Figura 4.26: Curva de magnetização do Fluxtrol 100 .............................................................. 73
Figura 4.27: Configuração do concentrador de fluxo: 1) com concentrador e 2) com
concentrador de fluxo parcial. .................................................................................................. 74
Figura 4.28: Comparação entre a potência absorvida durante o aquecimento para as três
configurações estudadas. .......................................................................................................... 74
Figura 4.29: Dureza final em Rockwell C para as configurações: a) com concentrador de
fluxo, b) com concentrador parcial e c) sem concentrador....................................................... 75
Figura 4.30: Potência absorvida em cada passo de tempo para diferentes valores de
frequência. ................................................................................................................................ 77
Figura 4.31: Profundidade da camada endurecida nas direções A e B após a têmpera para
diferentes valores de frequência. .............................................................................................. 77
Figura 4.32: Potência absorvida em cada passo de tempo para diferentes amplitudes de
corrente. .................................................................................................................................... 78
Figura 4.33: Profundidade da camada endurecida nas direções A e B após a têmpera para
diferentes valores amplitude de corrente. ................................................................................. 79
Figura 4.34: Comparação entre o mapa de dureza pós tempera para os casos: a) utilização de
um valor de dureza constante para cada fase b) utilização das equações de Maynier para
determinação da dureza de cada fase. ....................................................................................... 80
Lista de Tabelas
Tabela 4.1: Razão entre a propriedade não dependente da temperatura e a propriedade da
austenita em 550°C. .................................................................................................................. 56
Tabela 4.2: Dados de entrada para a simulação do moente. ..................................................... 70
Tabela 4.3: Resistividade elétrica, calor específico e condutibilidade térmica do aço SAE
1080 utilizadas durante a simulação. ........................................................................................ 70
Tabela 4.4: Densidade do aço SAE 1080. ................................................................................ 71
Tabela 4.5: Densidade do aço SAE 1080 em função da microestrutura e da temperatura. ...... 71
Tabela 4.6: Calor específico do aço SAE 1080 em função da microestrutura e da temperatura.
.................................................................................................................................................. 72
Tabela 4.7: Condutibilidade térmica do aço SAE 1080 em função da microestrutura e da
temperatura. .............................................................................................................................. 72
Lista de Abreviaturas e Siglas
Letras latinas
{𝐽𝑒𝑖}∗ - vetor complexo conjugado de {𝐽𝑒𝑖} no elemento de integração no ponto 𝑖
{𝐹ℎ𝑡} - vetor de carga nodal equivalente devido à convecção
{𝐹𝑞𝑓} - vetor de carga nodal equivalente devido ao fluxo de calor imposto
𝐻𝑚𝑖 - entalpia da fase 𝑖
[𝐻𝑡] - matriz de condutividade térmica do elemento
{𝐽𝑉} - densidade de corrente gerada devido à movimentação do componente mecânico dentro
do campo eletromagnético
{𝐽𝑡}- densidade de corrente total
[𝑁𝐴] - matriz das funções de forma do elemento
��𝑡 - calor gerado por unidade de volume na transformação de fase
ℎ𝑐𝑛 - coeficiente convectivo
ℎ𝑓 - coeficiente convectivo
ℎ𝑟𝑎𝑑 - coeficiente convectivo equivalente
[𝐴] - matriz de calor específico do elemento
𝐵𝑒𝑓𝑓 - densidade efetiva de fluxo magnético
𝐹𝑒 - fração de volume de final de transformação
𝐹𝑖 - fração de já transformada da fase 𝑖
𝐹𝑖 - fração volumétrica da fase 𝑖 formada.
𝐹𝑚,𝑗 - fração de matensita no passo de tempo 𝑗
𝐹𝑠 - fração de volume de início de transformação
𝐻𝑚- valor de pico do campo magnético
𝑄𝑐𝑜𝑛𝑣′ - taxa de calor perdido por convecção natural por unidade de comprimento
𝑄𝑟𝑎𝑑′ - taxa de calor perdido por radiação por unidade de comprimento
𝑇∞ - temperatura do ambiente
𝑇𝑎𝑟 - temperatura do ar envolta do cilindro
𝑇𝑏 - temperatura do fluido
𝑇𝑐𝑢𝑟𝑖𝑒 - temperatura de Curie
𝑇𝑖 – temperatura inicial
𝑇𝑚𝑠 - temperatura de início de transformação da martensita
𝑇𝑠 - temperatura da superfície
𝑛𝑖 - número de pontos de integração
�� - calor gerado por unidade de volume
�� - fluxo de calor imposto
�� - termo fonte
{𝑞} - vetor de fluxo de calor
[𝐻ℎ𝑓] - matriz de aproximação da parcela de superfície do termo convectivo
[𝐷] - matriz de condutividade
[𝑣] - matriz de relutância
[𝜇] - matriz de permeabilidade magnética
{𝐴��} - variação no tempo do vetor de potencial magnético
{𝐴𝑒} - potencial nodal do vetor magnético
{𝐹𝑞} - vetor de geração de calor
{𝐽𝑆} - densidade de corrente gerada por uma diferença de potencial na proximidade do
componente mecânico
{𝐽𝑒} - parcela da densidade de corrente relativa a uma corrente aplicada próxima ao
componente mecânico
{𝑇𝑒} - vetor das temperaturas nodais do elemento
{E} - intensidade do campo elétrico
{J} - densidade de corrente de condução
{N} - vetor das funções de forma
{𝐵} - densidade de fluxo magnético
{𝐷} - densidade de fluxo elétrico
{𝐻} - intensidade do campo magnético
{𝑛} - vetor unitário normal
ℎ - profundidade da camada temperada
𝐼 – corrente elétrica
𝑃 - período de tempo de oscilação da onda eletromagnética
𝑃𝑟 - número de Prandtl
𝑅𝑎𝑑 - número de Rayleigh
𝑇 - temperatura
𝑎 - coeficiente de difusão
𝑑 - diâmetro do cilindro
𝑒 – número neperiano
𝑔 - aceleração gravitacional
𝑘 - condutividade térmica do material
𝑙 - comprimento do cilindro
𝑛 – coeficiente do modelo de JMAK
Letras gregas
μ0 - permeabilidade magnética do vácuo
𝛿𝑝 - profundidade de penetração da onda eletromagnética
𝜇20- permeabilidade magnética à 20°C
𝜌𝑎𝑟 - densidade do ar
𝜌𝑒 - densidade de cargas elétricas
𝜎𝑟 - constante de Stefan-Boltzmann
𝜏𝑒 - tempo de final de transformação
𝜏𝑠- tempo de início de transformação
[𝜌] - matriz de resistividade elétrica
[𝜎] - matriz de condutividade elétrica
∆𝑡𝑇𝑜𝑡𝑎𝑙 - tempo total de aquecimento
∆𝑡𝑠𝑡𝑒𝑝 – passo de tempo na simulação da fase térmica I
∆𝑡𝑠𝑡𝑒𝑝_𝑟𝑒𝑠𝑓 - passo de tempo da fase térmica III
δ - operador variacional
𝛼 - difusividade térmica do ar
𝛽 - coeficiente de expansão térmica do ar
휀 - emissividade do material
𝜇 - permeabilidade magnética
𝜈 - viscosidade cinemática do ar
𝜌 - massa específica
𝜎 - condutividade elétrica
𝜔 – frequência de oscilação
𝜙 - propriedade genérica
Siglas
JMAK - Johnson-Mehl-Avrami-Komogorov
Sumário
Lista de Figuras
Lista de Tabelas
Lista de Abreviaturas e Siglas
Sumário .................................................................................................................................. xvi
1 INTRODUÇÃO ........................................................................................................... 18
1.1 Escopo Geral e Motivações ..................................................................................................... 18
1.2 Revisão da Literatura .............................................................................................................. 21
1.3 Objetivos e Contribuições ....................................................................................................... 26
2 MODELAGEM MECÂNICA MULTIFÍSICA ........................................................ 27
2.1 Método dos Elementos Finitos aplicado ao problema térmico ............................................... 27
2.2 Método dos Elementos Finitos aplicado ao problema eletromagnético .................................. 32
2.3 Modelo de Johnson-Mehl-Avrami-Komogorov e Regra da Aditividade de Sheil .................. 36
2.4 Modelo de Koistinen-Marburguer ........................................................................................... 39
3 ESTRATÉGIA DE SIMULAÇÃO ............................................................................ 41
3.1 Aquecimento ........................................................................................................................... 41
3.2 Resfriamento ........................................................................................................................... 45
3.3 Simulação multifísica integrada .............................................................................................. 46
4 RESULTADOS NÚMERICOS E ANÁLISE PARAMÉTRICA ............................ 48
4.1 Avaliação numérica da modelagem da convecção natural ...................................................... 48
4.2 Teste numérico do modelo microestrutural ............................................................................. 51
4.3 Influência das propriedades na fase de resfriamento ............................................................... 54
4.4 Simulação completa do processo integrado para o caso de um cilindro ................................. 60
4.5 Simulação da têmpera por indução em um moente de um virabrequim ................................. 65
4.5.1 Condições iniciais e de contorno do problema e propriedades dos materiais ...................... 67
4.5.2 Influência do Concentrador de Fluxo Magnético ................................................................. 73
4.5.3 Influência da frequência de oscilação da corrente do indutor .............................................. 75
4.5.4 Influência da amplitude da corrente no indutor .................................................................... 78
4.5.6 Influência da utilização das equações de Maynier para cálculo da dureza .......................... 79
5 CONCLUSÕES ............................................................................................................ 82
5.1 Sugestões de continuidade ...................................................................................................... 83
REFERÊNCIAS
APÊNDICES
ANEXO A – Artigo publicado no Congresso SAE 2014.............................................................. 90
ANEXO B – Artigo publicado no Congresso SAE 2015 .............................................................. 97
18
1 INTRODUÇÃO
1.1 Escopo Geral e Motivações
Na indústria automotiva, componentes mecânicos sujeitos à carregamentos dinâmicos e
contatos mecânicos durante a sua vida útil, podem apresentar falhas por fadiga e desgaste.
Para mitigar estas falhas os componentes podem ser submetidos a um processo de têmpera
por indução a fim de se alterar a microestrutura do material e, consequentemente, suas
características mecânicas relativas à resistência à fadiga e ao desgaste.
Figura 1.1: Geração de um campo eletromagnético através de um indutor em um dado instante
de tempo.
Globalmente, no processo de têmpera por indução, o componente mecânico é submetido
a um aquecimento controlado, usando-se sistemas por indução ou fornos. Após uma região do
corpo atingir temperaturas acima da temperatura de austenitização, cessa-se o processo de
aquecimento e inicia-se a fase de resfriamento. Este último é realizado por imersão em
tanques com líquidos, permanecendo enquanto as transformações microestruturais do material
estejam se desenvolvendo. O processo é finalizado e na camada endurecida ficam
estabelecidos: o perfil microestrutural do material da peça e seu estado de tensões. Neste
processo, o consumo de energia é elevado, e seu custo é significativo em relação ao custo
total da peça.
19
A Figura 1.1 mostra um esquema simplificado de um indutor gerando um campo
eletromagnético através de uma corrente elétrica I percorrendo um indutor elétrico. Sabe-se
que um campo eletromagnético é sempre gerado pela presença de uma corrente elétrica,
assim, é introduzida uma corrente elétrica de comportamento senoidal no indutor. O
componente mecânico a ser aquecido é então posicionado no interior do indutor a fim de que
o campo eletromagnético gerado incida sobre ele. Como a corrente do indutor é variável no
tempo, o campo eletromagnético gerado por ela também é variável no tempo, e a variação do
campo eletromagnético em um condutor elétrico causa o aparecimento de correntes na
superfície do mesmo, conhecidas como correntes parasitas. Tais correntes são responsáveis
pela geração de calor devido ao Efeito Joule.
Figura 1.2: Esquema representando as tensões residuais σxx geradas pelo processo de
têmpera.
Alguns dos parâmetros mais importantes resultantes do processo de têmpera por
indução são: a profundidade da camada temperada h, a dureza da mesma e as tensões
residuais impostas pelo processo de têmpera. O aumento da dureza superficial é um fator
importante, visto que as trincas geradas durante o processo de fadiga iniciam-se na superfície
do componente, assim, com o aumento da dureza superficial tem-se uma redução no desgaste
superficial e consequentemente na tendência à formação de trincas nesta região. As tensões
residuais também influem diretamente na resistência à fadiga do componente, uma vez que
com a indução de tensões compressivas residuais na superfície, as tensões de flexão cíclicas
20
que atuam no processo de nucleação e propagação de trincas na superfície do componente
mecânico são atenuadas e a vida útil do componente é prolongada.
A Figura 1.2 representa o comportamento típico das tensões residuais em um
componente mecânico submetido ao processo de têmpera por indução eletromagnética, na
qual a área com hachuras representa a camada superficial endurecida e a área em cinza
representa o núcleo não temperado. Pode-se notar que a presença de tensões compressivas na
região endurecida diminui a probabilidade de formação e propagação de trincas geradas por
cargas cíclicas.
Do ponto de vista da estrutura dos materiais, o processo de têmpera em aços consiste
em aquecer um componente até uma temperatura superior à temperatura de austenitização e
resfriá-lo rapidamente através da imersão do mesmo em um líquido, que poderá ser água, óleo
ou uma solução polimérica, promovendo desta forma a transformação microestrutural das
fases: martensita, perlita, bainita e ferrita.
Conforme mencionado, o processo pode ser dividido em duas etapas principais: o
aquecimento e o resfriamento. Durante a primeira etapa, existem duas maneiras principais
para aquecer os componentes: aquecimento por indução eletromagnética e aquecimento em
forno. O primeiro processo possui uma eficiência energética maior, visto que é possível
aquecer somente a região de interesse do componente mecânico, além de oferecer um maior
controle do processo de aquecimento. Todavia, o custo de implementação e a complexidade
do ajuste de parâmetros do processo de aquecimento por indução eletromagnética são
superiores ao do aquecimento em forno.
No aquecimento por indução eletromagnética o componente é introduzido no interior de
uma bobina de cobre na qual percorre uma corrente elétrica alternada, geralmente de alta
frequência. A corrente elétrica que percorre a bobina produz um campo eletromagnético que
oscila na mesma frequência da corrente. Este campo penetra no componente mecânico e,
como o mesmo é constituído de um material condutor elétrico, são geradas correntes elétricas
(correntes parasitas) na superfície do componente no sentido de contrapor a variação do
campo magnético. Estas correntes percorrem o material, dissipando energia por Efeito Joule,
aquecendo assim a área superficial do componente mecânico que está submetida ao campo
magnético gerado pela bobina. O calor gerado pelo Efeito Joule concentrado na superfície é
então difundido por condução térmica, distribuindo-se pelo volume do componente
respeitando-se um conjunto de condições de contorno do tipo convecção térmica.
21
O aquecimento produz uma mudança de fase nas regiões do componente mecânico cuja
temperatura final excede a temperatura de austenitização do aço. Após o final do
aquecimento, o componente é então resfriado rapidamente, e esta região do componente
mecânico que possui austenita como microestrutura pós-aquecimento poderá se transformar
novamente, produzindo perlita, bainita, martensita e ferrita. Dentro da região austenitizada,
quanto mais distante da superfície do componente mecânico, o resfriamento é mais lento e há
o favorecimento da formação de perlita, transformação essa que depende do tempo e da
temperatura. Em regiões mais próximas da superfície do componente não há tempo suficiente
para a difusão do carbono, e este acaba aprisionado e forma uma estrutura metaestável
chamada de martensita.
A martensita é uma estrutura com dureza bastante elevada, mas bastante frágil. Desta
maneira é importante controlar a profundidade da camada temperada, ou seja, com alta
dureza, para que o componente mecânico não perca sua capacidade de absorver impactos.
Para satisfazer esta condição, após a têmpera, os componentes mecânicos são submetidos ao
processo de revenimento, que visa reestabelecer a tenacidade do material, mantendo a dureza
e a profundidade da camada temperada dentro de certos limites de projeto.
Neste trabalho, porém, não serão tratados os efeitos do revenimento. Será discutido
somente o aquecimento por indução eletromagnética, o resfriamento, a microestrutura e a
dureza ao final do processo de têmpera.
1.2 Revisão da Literatura
Tendo em vista a importância do processo de têmpera e para que se possa aprimorá-lo, é
necessário compreender os principais fenômenos envolvidos: aquecimento por indução
eletromagnética, mudança de fase durante o aquecimento e o resfriamento, troca de calor
entre o corpo e a vizinhança. Dessa forma, é possível realizar simulações utilizando métodos
computacionais para se obter o perfil de tensões residuais e de dureza, possibilitando a
otimização do processo e a redução dos custos experimentais e do tempo de desenvolvimento
de vários projetos que utilizam a técnica de têmpera por indução. Além disso, é possível
incluir as restrições do processo de manufatura, no caso o tratamento térmico na etapa de
projeto, usando o conceito de “projetar para”, no qual os componentes são otimizados
22
levando-se em conta a influência dos processos de fabricação. Tal enfoque permite encontrar
peças de melhor desempenho, menor peso, respeitando-se os requisitos do projeto.
A compreensão dos parâmetros do processo e a influência desses nas características
mecânicas do produto final é de suma importância para a minimização dos custos de
manufatura e do tempo de produção permitindo a obtenção de produtos com maior qualidade.
Os parâmetros do processo de têmpera que podem ser mais facilmente alterados - e
consequentemente mais estudados - são: tempo de aquecimento, potência e frequência da
corrente no indutor. Um ponto observado por Kristoffersen e Vomacka (2001) foi a obtenção
de camadas endurecidas mais profundas e um gradiente de tensão menor na zona de transição
quando utilizada uma frequência mais alta na corrente do indutor. Na prática, a espessura
mínima para uma região pré-estabelecida é adotada e o processo deve ser controlado de forma
a se obter um valor o mais próximo possível da meta, para um consumo mínimo de energia.
Neste caso, deve ser salientado que o consumo de energia no processo de têmpera é
considerável em relação ao custo total do componente.
O método dos elementos finitos é amplamente utilizado na simulação numérica dos
fenômenos de transferência do campo eletromagnético para energia térmica (Guo et al., 2012;
Ruan, 1997; Toparli et al., 2002; Cheng at al., 1998; Cho, 2012; Drobenko et al., 2007; Di
Luozzo et al., 2012; Palin-Luc et al., 2011) .Destacam-se também a utilização de outros
métodos, como método dos elementos de contorno (Cajner et al., 2004) e redes neurais
(Toparli et al., 2002). Neste trabalho foi adotado método dos elementos finitos para a solução
numérica do problema.
Para a correta utilização do método dos elementos finitos para a simulação da indução
eletromagnética é necessário que se verifique a influência do refinamento da malha sobre os
resultados obtidos. Neste caso são modelados os indutores, o meio de transmissão e as peças a
serem tratadas. A simulação no domínio do tempo envolve custos computacionais
elevadíssimos e, muitas vezes inviáveis. Dessa forma, adotam-se métodos de simulação no
domínio da frequência. (Ansys, 2007). No entanto, estes processo geralmente ocorrem na
faixa de frequência de 10 khz, o que conduz a comprimentos de onda muito pequenos, e
portanto malhas com elevado grau de refinamento (Cajner et al., 2004; Ansys, 2007).
Além disso, na simulação de aquecimento de cilindros por indução eletromagnética, foi
verificado que as dimensões da malha de ar que envolve o cilindro devem ser duas vezes
maior longitudinalmente e duas vezes e meia maior radialmente do que as dimensões do
23
cilindro (Drobenko et al., 2007), para evitar a influência das condições de contorno no
problema. Tal fato também implica em malhas de grandes proporções.
Para a simulação do campo de temperaturas através do calor gerado por efeito joule na
peça aquecida, foi verificada a necessidade da utilização de uma malha mais refinada na
superfície do cilindro, contendo de 3 a 10 elementos na região superficial de geração de calor,
sendo que esta apresenta comprimento de 5 profundidades de penetração do campo
eletromagnético (Cajner et al., 2004). A profundidade de penetração do campo
eletromagnético em um condutor elétrico é dependente da frequência de oscilação do campo,
da condutividade elétrica do meio e da permeabilidade magnética do meio e pode ser definida
como a profundidade na qual a energia do campo eletromagnético equivale a 1/e, onde e é o
número neperiano, da energia do mesmo na superfície do meio condutor. É válido salientar,
também, que com a perda das propriedades ferromagnéticas do aço em temperaturas acima do
ponto de Curie, a posição de máximo calor gerado passa a se localizar cada vez mais no
interior do cilindro, onde o ponto de Curie ainda não foi ultrapassado, ou seja, onde o material
ainda possui suas propriedades ferromagnéticas (Cho, 2012; Drobenko et al., 2007).
Assim, outro fator importante a se considerar é a variação das propriedades térmicas e
eletromagnéticas com a temperatura. Para temperaturas maiores que 300ºC, para a
condutividade elétrica, e 600ºC, para a permeabilidade magnética, a não consideração da
variação da propriedade física com a temperatura, implica em erros consideráveis no campo
de temperaturas (Drobenko et al., 2007). Dentre as propriedades resistividade elétrica, calor
específico, condutibilidade térmica e permeabilidade magnética, foi observado, através de um
estudo de sensibilidade utilizando um software comercial de elementos finitos, que as duas
primeiras são as que mais interferem na profundidade da camada endurecida (Barka et al.,
2007).
Devido ao comportamento senoidal da corrente no indutor, o campo eletromagnético no
interior do cilindro ferromagnético é tratado por alguns autores como sendo também senoidal
e a análise conduzida como sendo harmônica e em regime quase estático, desprezando-se os
efeitos transientes do campo. Porém, é sabido que permeabilidade magnética em materiais
ferromagnéticos é uma função não linear da intensidade do campo magnético, logo, a
simulação transiente do fenômeno conduz a resultados mais precisos (Drobenko et al., 2007).
Como a solução de campos harmônicos é mais simples e requer menor esforço
computacional, alguns autores propõem modelos em que se obtêm uma permeabilidade
24
equivalente, que seja função apenas da temperatura, de forma a obter uma distribuição de
correntes parasitas e, consequentemente, de temperaturas semelhantes às da simulação
transiente, considerando a permeabilidade magnética como função da temperatura e da
intensidade do campo magnético. (Ansys, 2007; Cajner et al., 2004; Drobenko et al., 2007; Di
Luozzo et al., 2012 )
Acompanhando-se o processo de resfriamento na têmpera, pode-se verificar a transição
das tensões inicialmente trativas na superfície para compressivas no decorrer do resfriamento.
Concluído o processo de resfriamento, as tensões permanecem como compressivas na
superfície e trativas no centro. Nota-se também, que na mesma região do cilindro em que há
um súbito aumento da tensão, ou seja, na zona de transição da tensão compressiva para
trativa, há também uma queda brusca da dureza, evidenciando uma correspondência entre o
perfil de dureza e de tensões residuais (Grum e Furlan, 1998; Grum, 2001; Palin-Luc et al.,
2011).
Com o objetivo de predizer a distribuição microestrutural após o processo de têmpera,
foram desenvolvidos modelos empíricos e teóricos. Além disso, são necessários modelos para
a dureza de cada fase formada. Na maioria dos estudos, os pesquisadores utilizam uma
formulação empírica dependente da fase formada, da taxa de resfriamento e da composição
química do aço em questão para obtenção do perfil de dureza gerado pela têmpera
(Magnabosco et al., 2006; Huiping et al., 2007; Carlone et al., 2010). Outros métodos
empregados são a utilização de valores fixos para a dureza de cada fase (Woodard et al.,
1999) e a utilização de valores experimentais obtidos para o aço que se deseja estudar (Lee et
al., 2010). Este último método é bastante interessante, já que é baseado em valores medidos
experimentalmente (Lee et al., 2010), pois a influência de um determinado elemento de liga
pode ser alterada na presença de outro elemento em específico (Doane, 1979), sendo difícil
uma determinação teórica desses valores.
Modelos derivados das equações de Johnson-Mehl-Avrami-Komogorov (JMAK) são
amplamente utilizados na simulação das transformações microestruturais que possuem caráter
difusional, ou seja, as transformações austenita-perlita, austenita-bainita, austenita-ferrita e
vice-versa. Porém, esses modelos levam em consideração transformações isotérmicas e, dessa
forma, para que a solução possa ser encontrada, as curvas de resfriamento são divididas em
uma série de passos isotérmicos conectados por uma pequena queda de temperatura
instantânea com frações volumétricas de fase constantes. A fim de se levar em conta o estágio
de nucleação da nova fase, é utilizada a Regra da Aditividade de Scheil (Hömberg, 1996,
25
Woodard et al., 1999; Kang e Im, 2005; Huiping et al., 2007; Carlone et al., 2010; Lee et al.,
2010).
Quando a transformação de fase é função apenas da temperatura - quando não depende
do tempo, como no caso da transformação austenita-martensita -, os modelos encontrados
durante esta revisão foram o de Koistinen-Marburguer e o de Yu. Porém, este primeiro não é
capaz de descrever a decomposição completa da austenita em martensita, pois não leva em
consideração a irreversibilidade da transformação. Quando isso não é levado em conta, as
curvas de resfriamento contínuo ficam deformadas por conta do calor latente liberado na
mudança de fase, sendo então necessário impedir que a porcentagem de martensita já
transformada possa ser diminuída durante o resfriamento (Hömberg, 1996).
O calor latente é levado em conta na maior parte dos trabalhos dedicados à simulação da
etapa de resfriamento, e, apesar de influenciar o histórico de temperatura mais intensamente
nas proximidades do eixo de simetria no caso de cilindros, a não consideração do mesmo
causa uma grande alteração no resultado obtido da microestrutura superficial do cilindro
(Woodard et al., 1999). Outra fonte de calor durante o processo de têmpera é aquela liberada
devido às deformações, porém essa representa apenas aproximadamente 3% do calor gerado e
pode ser negligenciada (Huiping et al., 2007).
Existem, ainda, alguns fatores que dificultam um cálculo mais preciso do perfil de
dureza pós-têmpera, como o tamanho do grão antes e depois da têmpera, a presença de
precipitado na austenita e a tensão interna. A última altera o tempo de incubação da perlita e a
temperatura de início da transformação da martensita. Apesar da influência desses fatores ser
geralmente negligenciada, vários autores encontram uma boa concordância dos resultados
obtidos com os medidos experimentalmente. Outro fator que possui grande influência no
perfil de dureza é o coeficiente de troca de calor durante o processo. Esse está intimamente
ligado ao meio e ao tipo de resfriamento adotado, além de ser função da temperatura e da
geometria da peça.
Outro problema a ser considerado é que após a têmpera, muitas vezes o componente
mecânico passa por um processo de retificação, que induz tensões de tração, reduzindo as
tensões residuais compressivas induzidas pelo processo de têmpera. Esse tópico é tratado por
Grum e seus colaboradores (Grum, 2001 e Grum e Ferlan, 1998), e não será discutido neste
trabalho.
26
Neste trabalho, propõe-se fazer uma simulação integrada das fases de aquecimento,
tempo de retardo e resfriamento, considerando a modelagem microestrutural e as propriedades
dependentes da temperatura.
Dentro dessa perspectiva o problema da têmpera por indução necessita de uma
abordagem multifísica, cuja obtenção da solução é complexa porém relevante para a indústria
de componentes mecânicos, sendo que seu entendimento implicará desenvolver componentes
com melhores propriedades mecânicas e com melhor desempenho em serviço. Outro ponto de
elevada relevância é a realização de um processo ou fase de aquecimento com maior
eficiência energética, além de uma redução da necessidade de realização de testes empíricos.
1.3 Objetivos e Contribuições
O objetivo geral deste trabalho é desenvolver uma metodologia de simulação multifísica
do aquecimento por indução eletromagnética e posterior têmpera de um componente
mecânico, além de investigar a influência dos seguintes parâmetros do processo: frequência e
amplitude da corrente do indutor elétrico, além da presença ou não do concentrador de fluxo
eletromagnético. O Método dos Elementos Finitos foi utilizado na solução dos problemas
estudados, usando-se a plataforma ANSYS e a linguagem APDL.
O acoplamento das diversas fases físicas do problema (eletromagnética e térmica,
térmica e microestrutural) e a implementação de modelos para descrever as transformações
microestruturais podem ser mencionadas como objetivos específicos do trabalho.
Como resultado específico deste trabalho, menciona-se a aplicação da metodologia
proposta a um caso real de um moente de virabrequim fabricado pela ThyssenKrupp
Metalúrgica Campo Limpo Ltda.
Pode-se mencionar como contribuições deste trabalho: a investigação da influência dos
parâmetros do processo como forma de auxiliar o setup do processo, o desenvolvimento de
uma metodologia de simulação envolvendo o acoplamento de simulações de físicas distintas e
o estudo da sensibilidade de parâmetros computacionais e de propriedades físicas.
27
2 MODELAGEM MECÂNICA MULTIFÍSICA
A simulação multifísica proposta neste trabalho pode ser dividida em duas etapas
principais: aquecimento, englobando a fase eletromagnética e o tempo de retardo, e
resfriamento, englobando a têmpera e a análise microestrutural metalúrgica. Na Figura 2.1
mostra-se o encadeamento das fases de simulação multifísica do processo.
Figura 2.1: Fases principais da simulação.
Na sequência deste capítulo, apresentam-se os princípios e hipóteses gerais adotadas
em cada módulo de simulação e as equações básicas utilizadas no processo e comentários
sobre os métodos de resolução utilizados. As equações do modelo térmico e eletromagnético
não serão abordadas de maneira aprofundada, visto que foram utilizados os módulos de
simulação térmica e eletromagnética do software comercial Ansys e que o objetivo do
capítulo é fornecer uma base teórica para a construção de metodologia de simulação
multifísica integrada do processo de têmpera.
2.1 Método dos Elementos Finitos aplicado ao problema térmico
Durante a simulação do problema acoplado foi utilizado o software comercial Ansys
para a obtenção do resultado e geração da malha de elementos finitos. O programa Ansys é
uma plataforma de uso geral que possibilita realizar análises estruturais, térmicas,
28
eletromagnéticas, acústica e outras. Trata-se de um código robusto que permite o uso de
materiais com propriedades que variam em função da temperatura e que dispõe de um
conjunto amplo de métodos de resolução de sistemas não-lineares e com vários recursos que
podem ser adaptados em função do problema que está sendo resolvido. Além disso, o
programa Ansys permite o desenvolvimento de novas rotinas ou funções através do uso de
uma linguagem própria, que permite acessar a maior parte da biblioteca de funções internas
do programa. Ressalta-se também o fato do programa Ansys ser amplamente utilizado no
meio industrial e, portanto, gerar um módulo para simulação de têmpera por indução no
mesmo, permite que vários usuários possam se beneficiar deste desenvolvimento.
Para a modelagem da etapa térmica, foi utilizada a equação de conservação de energia
com a equação de condução de calor e a lei de Fourier na sua forma matricial já considerando
uma discretização e uma a aproximação por elementos finitos, como se segue (Toparli et al,.
2002),(Carlone et al., 2010), (Ansys, 2007):
{q} = −[D]{L}T(x,y) (2.1)
ρT = {L}T([D]{L}T(x,y)) + q (2.2)
[D] = [
kxx 0 00 kyy 0
0 0 kzz
] (2.3)
{L} =
{
∂
∂x∂
∂y
∂
∂z}
(2.4)
Aqui [D] é a matriz de condutividade e {q} é o vetor de fluxo de calor, T(x,y) é a
temperatura em um ponto qualquer, ρ é a massa específica e q é termo fonte, e a derivada com
relação ao tempo é indicada pelo ponto acima da variável.
As duas condições de contorno aplicadas a este problema foram a imposição do fluxo
de calor em uma superfície e a imposição da convecção em outra superfície. A radiação foi
29
modelada como uma convecção equivalente cujo coeficiente convectivo é função da
temperatura. As duas condições de contorno podem ser equacionadas com as equações (2.5) e
(2.6) respectivamente.
{q}T{n} = −q (2.5)
{q}T{n} = hf (Ts − Tb) (2.6)
Sendo {n} o vetor unitário normal, hf o coeficiente convectivo, q o fluxo de calor
imposto, Tb a temperatura do fluido ou do meio, ou no caso da radiação modelada como
convecção equivalente, a temperatura do ambiente e Ts a temperatura da superfície onde
ocorre a troca de calor.
Utilizando novamente a Lei de Fourier para condução de calor, equação (2.1),
substituindo nas equações (2.5) e (2.6), tem-se:
{n}T[D]{L}T = q (2.7)
{n}T[D]{L}T = hf(Tb − Ts) (2.8)
Pré multiplicando a equação (2.2) por uma variação virtual de temperatura,
substituindo as equações (2.7) e (2.8) e integrando no volume chega-se a uma forma fraca do
problema, que é dada por (Toparli et al,. 2002; Carlone et al., 2010; Ansys, 2007):
∫ (ρcδTT + {L}T(δT)([D]{L}T))dV = ∫ δTqdSSV
+ ∫ δThf(Tb − T)dSS+
∫ δTqdVv
(2.9)
Fazendo-se uma discretização do domínio, a temperatura de cada elemento pode ser
escrita como uma função de x, y e t, e pode ser expressa de acordo com a equação (2.10).
T(x,y) = {N}T{Te} (2.10)
30
Sendo {N} o vetor das funções de forma e {Te} o vetor das temperaturas nodais do
elemento. Neste trabalho usou-se elementos triangulares de 6 nós e preferencialmente, na
região de interesse, elementos quadrilaterais de 8 nós, cujas funções de interpolação são dadas
(Ansys, 2007):
{N}T =1
4{TI(1 − s)(1 − p)(−s − p − 1), TJ(1 + s)(1 − p)(s − p − 1),
TK(1 + s)(1 + p)(s + p − 1), TL(1 − s)(1 + p)(−s + p − 1), 2TM(1 − s2)(1 − p),
2TN(1 + s)(1 − p2), 2TO(1 − s
2)(1 + p), 2TP(1 − s)(1 − p2)}
(2.11)
As coordenadas p e s são normalizadas, ou seja, variando de -1 a 1, mas não são
necessariamente ortogonais entre si. O referencial (x, y) e o elemento finito utilizado estão
definidos na Figura 2.2.
Figura 2.2: Elemento finito de 8 nós utilizado na simulação da fase térmica e eletromagnética.
Aplicando o operador variacional na equação (2.10) e derivando a equação (2.10) com
relação ao tempo, tem-se, respectivamente:
δT = {δTe}T{N} (2.12)
T = {N}T{Te} (2.13)
31
Substituindo as equações (2.10), (2.12) e (2.13) na equação (2.9):
∫ ρc{δTe}T{N}{N}T{Te}V
dV + ∫ {δTe}T[B]T[D][B]{Te}dVV
= ∫ {δTe}T{N}qdS
S+
∫ {δTe}T{N}hf(Tb − {N}
T{Te})dSS+ ∫ {δTe}
T{N}qdVV
(2.14)
Sendo a matriz [B] definida por:
[B] = {L}{N}T (2.15)
Como ρ é assumido como constante dentro do elemento e, {Te}, {Te} e {δTe} são quantidades
nodais e não variam dentro do elemento, pode-se reescrever a equação acima como:
ρ ∫ c{N}{N}TV
dV{Te} + ρ∫ [B]T[D][B]dV{Te}V= ∫ {N}qdS
S+ ∫ Tbhf{N}dSS
−
∫ hf{N}{N}T{Te}S
dS + ∫ {N}qdVV
(2.16)
Esta equação representa a forma semi discretizada de elementos finitos que será
utilizada para resolver as fases térmicas no processo de simulação implementado. De forma
compacta pode-se escrever:
[A]{Te} + [Ht]{Te} = {Fqf} + {Fht} + [Hhf]{Te} + {Fq} (2.17)
sendo [A] a matriz de calor específico do elemento, [Ht] a matriz de condutividade térmica do
elemento, {Fqf} é o vetor de carga nodal equivalente devido ao fluxo de calor imposto, {Fht}
é o vetor de carga nodal equivalente devido à convecção, [Hhf] é a matriz de aproximação da
parcela de superfície do termo convectivo e {Fq} é o vetor de geração de calor.
32
2.2 Método dos Elementos Finitos aplicado ao problema eletromagnético
O problema eletromagnético é descrito pelas leis de Maxwell e pelas leis constitutivas
do material, descritas abaixo:
{∇}x{H} = {J} +𝛛𝐃
𝛛𝐭 (2.18)
{∇}x{E} = − 𝛛𝐁
𝛛𝐭 (2.19)
{∇}. {B} = 0 (2.20)
{∇}. {D} = ρe (2.21)
onde {H} é a intensidade do campo magnético, {E} é a intensidade do campo elétrico, {J} é a
densidade de corrente de condução, {D} é a densidade de fluxo elétrico, ρe é a densidade de
cargas elétricas e {B} é a densidade de fluxo magnético, o operador do produto interno é dado
por . e o produto matricial é dado por x. Em materiais bons condutores e considerando o
problema como sendo harmônico magneto-quasi-estático, as correntes de deslocamento são
negligenciáveis quando comparadas as correntes de condução. A corrente de deslocamento é a
taxa de variação do vetor densidade de fluxo elétrico, representada pelo termo 𝛛𝐃
𝛛𝐭 da equação
(2.18). Então, pode-se reescrever as equações de Maxwell da seguinte forma:
{∇}x{H} = {J} (2.22)
{∇}x{E} = − {B} (2.23)
{∇}. {B} = 0 (2.24)
33
As equações constitutivas descritas abaixo suplementam as equações (2.22), (2.23) e
(2.24) e descrevem o comportamento dos materiais eletromagnéticos.
{J} = σ{E} (2.25)
{B} = μ{H} (2.26)
Aqui, σ é a condutividade elétrica e μ é a permeabilidade magnética. No presente
modelo, μ é função da temperatura e da intensidade do campo magnético {H}. A equação
(2.27), semelhante a utilizada por (Drobenko et al., 2007) é uma aproximação utilizada para
considerar o efeito da temperatura na curva de magnetização, curva B-H.
{B} = μ(H, T){H} = μ20(H) [1 − (T
Tcurie)6
] {H} (2.27)
onde μ20 é obtido da curva B-H em 20°C e Tcurie é a temperatura de Curie do material, na qual
o aço perde suas propriedades ferromagnéticas, isto é, a permeabilidade magnética do aço
acima da temperatura de Curie é considerada igual a permeabilidade magnética do vácuo. A
Figura 2.3 ilustra o efeito desta correção na curva de magnetização do aço AISI 1045. Nota-se
que com a aproximação da temperatura do componente com Tcurie = 770°C há uma
diminuição da densidade de fluxo magnético no material e consequentemente da
permeabilidade magnética, até o limite quando μ(T, H) = μ0, onde μ0 é a permeabilidade
magnética do vácuo.
Apesar do campo eletromagnético gerado no indutor ser senoidal, a não linearidade da
permeabilidade magnética do material ferromagnético faz com que o campo dentro do mesmo
seja não senoidal. Entretanto, para problemas cujo objetivo é a obtenção da potência
eletromagnética média dissipada e não há interesse na obtenção do campo eletromagnético,
podemos substituir o material ferromagnético por um material fictício baseado na
equivalência de energia, em que o campo eletromagnético transiente também é substituído por
um campo harmônico que incidindo sobre o material fictício, produza a mesma distribuição
de correntes parasitas do campo transiente incidindo sobre o material ferromagnético.
Considerando, portanto, o problema como harmônico no tempo, é possível reduzir o custo
computacional mantendo-se uma boa precisão (Ansys, 2007). Neste caso, a intensidade de
campo magnético pode ser aproximada utilizando a equação abaixo:
34
1
2∫ HmdBeffBeff0
= 4
P∫ (∫ Hm sin(ωt) dB
B
0)
P
40
dt (2.28)
onde Beff é a densidade efetiva de fluxo magnético, P é o período de tempo de oscilação da
onda eletromagnética, Hm é o valor de pico do campo magnético e ω é a frequência de
oscilação do campo eletromagnético
Figura 2.3: Densidade do fluxo magnético em tesla em função da temperatura versus a
intensidade do campo magnético em A/m para o aço AISI 1045, adaptado de Li et al (2012).
Durante a solução do problema através do método dos elementos finitos, os cálculos do
fluxo magnético e da intensidade do fluxo magnético são realizados da seguinte maneira:
{B} = {∇}x[NA]T{Ae} (2.29)
{H} = [v]{B} (2.30)
[v] = [μ]−1 (2.31)
0
0,5
1
1,5
2
2,5
0 2000 4000 6000 8000 10000
B [
T]
H [A/m]
20°C
400°C
500°C
600°C
685°C
710°C
730°C
750°C
760°C
765°C
35
Onde [NA] é a matriz das funções de forma do elemento, [v] é a matriz de relutância,
[μ] é a matriz de permeabilidade magnética e {Ae} é o potencial nodal do vetor magnético. A
densidade total da corrente no componente é avaliada como se segue:
{Jt} = {Je} + {JS} + {JV} (2.32)
onde {Jt} é a densidade de corrente total, {Je} é a parcela da densidade de corrente relativa a
uma corrente aplicada próxima ao componente mecânico, {JS} é a densidade de corrente
gerada por uma diferença de potencial, que gera uma corrente elétrica, na proximidade do
componente mecânico e, {JV} é a densidade de corrente gerada devido à movimentação do
componente mecânico dentro do campo eletromagnético.
Como não há deslocamento do componente dentro do campo magnético, nem uma
diferença de potencial aplicada nas proximidades do componente mecânico, {JV} = {JS} = 0.
A equação (2.33) descreve o cálculo da corrente gerada no componente mecânico devido à
proximidade de outra corrente. Aqui {Ae} é a variação no tempo do vetor de potencial
magnético, [σ] é a matriz de condutividade elétrica e ni o número de pontos de integração.
{Jt} = {Je} = −[σ]1
ni∑ [NA]
T{Ae}nii=1 (2.33)
Neste caso, o calor por unidade de volume gerado por efeito Joule durante a análise
harmônica é calculado (Ansys, 2007):
q = Re(1
2ni∑ [ρ]{Jei}{Jei}
∗nii=1 ) (2.34)
onde [ρ] é a matriz de resistividade elétrica e {Jei}∗ é a o vetor complexo conjugado de {Jei}
no elemento de integração no ponto i, e q é o calor gerado por unidade de volume utilizado
com dado de entrada na simulação térmica.
36
2.3 Modelo de Johnson-Mehl-Avrami-Komogorov e Regra da Aditividade
de Sheil
Além dos modelo térmico e eletromagnético, integrou-se na modelagem uma fase
microestrutural para representar as transformações metalúrgicas ocorridas em função das
grandes variações de temperatura envolvidas no processo de têmpera por indução. Durante o
resfriamento diferentes fases são formadas, tais como perlita, bainita, ferrita e martensita,
sendo necessário a simulação acoplada das fases microestrutural e térmica.
As transformações austenita-perlita, austenita-bainita e austenita-ferrita são difusionais,
ou seja, necessitam de um tempo finito para que a transformação ocorra. É utilizado o
modelo de Johnson–Mehl–Avrami–Kolmogorov (JMAK) para modelar esse tipo de
transformação (Carlone et al., 2010; Woodard et al., 1999; Hömberg, 1996):
Fi(T) = 1 − exp [−a(T). t(T)n(T)] (2.35)
Figura 2.4: Diagrama TTT do aço SAE 1080, adaptado de Woodard et al., (1999).
37
Onde Fi é a fração volumétrica da fase i formada. O coeficiente de difusão, a, e o
expoente, n, são propriedades do material e podem ser obtidos pelo diagrama TTT, utilizando
τs(T) e τe(T), tempo de início de transformação e de fim de transformação para uma dada
temperatura, como segue:
n(T) =ln (ln(1−Fs)/ln (1−Fe))
ln (τs(T))−ln (τe(T)) (2.36)
a(T) = − ln(Fe) τs(T)−n(T) (2.37)
onde Fs e Fe são a fração de volume de início e final de transformação e são assumidos como
0,01 e 0,99 respectivamente. Utilizando as equações (34) e (35) e o diagrama TTT do aço
SAE 1080 é possível obter o coeficiente de difusão a e o expoente do modelo JMAK, n, para
o aço SAE 1080, mostrados nas figuras 2.5 e 2.6 respectivamente.
Figura 2.5: Coeficiente de difusão em função da temperatura do aço SAE1080, baseado em
Woodard et al. (1999).
0
0,005
0,01
0,015
0,02
0,025
0 200 400 600 800
Coef
icie
nte
de
Dif
usã
o a
Temperatura [°C]
38
Figura 2.6: Expoente n usado na equação de JMAK do aço SAE 1080 baseado em Woodard et
al. (1999).
O modelo JMAK é válido somente para transformações isotérmicas, então, para usar
este modelo, deve-se converter a curva de resfriamento para uma sequência de passos
isotérmicos como mostrada na Figura 2.7.
Figura 2.7: Diagrama TTT esquemático mostrando a conversão da curva de resfriamento em
passos isotérmicos.
0
0,5
1
1,5
2
2,5
3
3,5
4
0 200 400 600 800
Ex
po
ente
n (
JMA
K)
Temperatura [°C]
39
Assim, utilizando o modelo JMAK modificado, o tempo utilizado na equação (2.38) é a
soma de um tempo fictício ao valor do passo de tempo utilizado. Este tempo fictício é baseado
na fração total volumétrica já transformada de uma dada fase i, e representa o tempo que seria
necessário para que toda a fração volumétrica já formada se transformasse isotermicamente
naquela temperatura. Portanto a equação incremental para o tempo é dado por:
tj = ∆tj + tj,fict (2.38)
sendo
tj,fict (Tj) = [− ln(1−Fi,j−1)
a(Tj)]
1
n(Tj) (2.39)
Fi,j(Tj) = 1 − exp [−a(Tj)t(Tj)n(Tj)
] (2.40)
As transformações fases perlita e bainita ocorrem após de um tempo de incubação. Para
considerar este aspecto, é utilizada a Regra da Aditividade de Sheil, e a transformação poderá
começar somente se a equação (2.41) for satisfeita, ou seja, se a condição abaixo for satisfeita
(Carlone et al., 2010; Woodard et al., 1999):
∑∆tj
τs(Tj)≥ 1n
j=1 (2.41)
2.4 Modelo de Koistinen-Marburguer
A transformação austenita-martensita não possui caráter difusional e pode ser
modelada pela equação de Koistinen-Marburguer. Deve-se levar em consideração a
porcentagem de austenita já transformada em perlita e bainita, então a equação de Koistinen-
Marburguer deve ser multiplicada pela fração de austenita remanescente, como segue:
40
Fm,j = [1 − exp[−0,011(Tms − T)]](1 − ∑ Fii ) (2.42)
onde Tms é a temperatura de início de transformação da martensita, Fm,j é a fração de
matensita no passo de tempo j, Fi é a fração já transformada da fase i.
Durante o resfriamento, para cada passo de tempo, as propriedades de densidade, calor
específico e condutividade térmica são atualizadas baseadas na fração volumétrica de cada
fase, assim as propriedades são função da fase e da temperatura e podem ser genericamente
representadas pela função ϕ, cujo valor atual depende do valor de cada fase ϕi, conforme a
equação abaixo:
ϕ = ∑ Fiϕini=1 (2.43)
Durante a transformação de fase há liberação de calor latente e este é incluído em cada
passo de tempo, considerando a variação da fração volumétrica e a variação de entalpia para
cada fase Hmi, como mostrado na equação (2.44). Como existe calor sendo liberado e, um
consequente aumento de temperatura, pode haver a tendência, durante a simulação
computacional, de uma diminuição da fração volumétrica de martensita já transformada.
Como este efeito não ocorre fisicamente, ou seja, durante o resfriamento não é possível
transformar martensita em austenita novamente, deve-se evitá-lo durante a simulação
mantendo sempre o valor máximo da fração volumétrica de martensita entre dois passos de
tempo consecutivos, de acordo com a equação (Hömberg, 1996) (2.45).
qt = ∑∆Hmi∆Fi,j
∆tj (2.44)
Fm,j = max[Fm,j , Fm,j−1] (2.45)
onde qt é o calor liberado por unidade de volume durante a transformação de fase.
41
3 ESTRATÉGIA DE SIMULAÇÃO
A simulação computacional foi dividida em duas partes principais: aquecimento e
resfriamento, sendo utilizado como condição inicial da segunda fase, o campo de
temperaturas obtido na primeira.
A primeira fase engloba a simulação do aquecimento por indução eletromagnética e o
tempo de retardo, enquanto a segunda consiste no resfriamento por imersão em líquido, que
nos casos tratados é água.
O tempo de retardo é o tempo que a peça fica suspensa no ar, após o indutor ser
desligado e antes de ser imersa em água, e tem a finalidade de permitir que o calor gerado na
superfície do componente mecânico, através da indução eletromagnética, se difunda
estendendo assim a região austenitizada.
Como condição inicial da simulação, foi considerado, em todos os casos simulados, que
o componente mecânico possui uma distribuição de temperatura homogênea, igual a
temperatura ambiente, e 100% de perlita como microestrutura.
3.1 Aquecimento
A simulação computacional da fase de aquecimento consiste em acoplar a simulação
eletromagnética e a térmica, além de utilizar o campo de temperaturas após o aquecimento
propriamente dito como condição inicial para a simulação do tempo de retardo. Dessa forma
pode-se dividir a fase de aquecimento em quatro fases: fase eletromagnética, fase térmica I,
fase térmica II, e fase microestrutural I, onde as fases térmicas I e II se referem,
respectivamente, à fase térmica acoplada a fase eletromagnética e à simulação do tempo de
retardo. A fase microestrutural I refere-se à avaliação da microestrutura existente após o
tempo de retardo.
A fase de aquecimento está estruturada como mostrado na Figura 3.1, sendo que uma
volta completa no loop corresponde a um passo de tempo da fase térmica I, ∆tstep. O primeiro
passo é inserir as condições de contorno eletromagnéticas e a temperatura de cada elemento
finito, que no passo de tempo inicial do loop da fase térmica I é igual a temperatura ambiente.
42
Na fase eletromagnética, a temperatura é utilizada para o programa definir, através de uma
tabela, o valor da resistividade elétrica e a correção da curva B-H de acordo com a equação
(2.27) para cada elemento finito, lembrando que caso a temperatura do elemento exceda a
Temperatura de Curie, a curva B-H é então substituída por um valor constante de
permeabilidade magnética, cujo valor é igual a permeabilidade magnética do vácuo. Como já
explicado no Capítulo 2, apesar do problema eletromagnético ser transiente, ele é tratado
como harmônico, e assim independente do tempo, através da simplificação produzida pela
equação (2.28). Neste caso é usada uma malha de elementos finitos que representa a bobina, o
concentrador de fluxo eletromagnético, o ar em torno da região de interesse e a peça que será
aquecida. A resolução deste problema harmônico produz como resultado a distribuição de
correntes parasitas na superfície do componente mecânico e também a potência, em forma de
calor, dissipada por elas.
Figura 3.1: Estratégia de simulação para a fase de aquecimento eletromagnético.
Completada a fase eletromagnética, altera-se o módulo do Ansys para a simulação de
problemas térmicos, alterando o tipo de elemento finito mas preservando a geometria. É então
necessária a eliminação dos elementos finitos que não serão utilizados na análise térmica: os
elementos pertencentes ao ar, a bobina e ao concentrador de fluxo eletromagnético. São
inseridos também as condições de contorno térmicas da fase térmica I e a potência dissipada
na fase eletromagnética, como geração de calor volumétrica devido ao efeito Joule e
calculadas com a equação (2.34), que permanece constante durante todo o passo de tempo da
43
simulação transiente da fase térmica I. Como resultado é obtido o campo de temperaturas que
será utilizado no próximo passo de tempo para definir as propriedades eletromagnéticas.
Retornando agora para o primeiro passo do loop de aquecimento mostrado na Figura
3.1, altera-se o elemento finito para aquele apropriado para a simulação eletromagnética de
forma a reestabelecer os elementos não usados na simulação térmica, mas necessário na
eletromagnética. É atualizado então a temperatura de cada elemento finito e assim as
propriedades eletromagnéticas dos mesmos. A simulação eletromagnética é realizada
novamente e uma nova potência dissipada é obtida.
Assim, é realizado o loop contendo as fases eletromagnética e térmica I com ∆tTotal/
∆tstep iterações, onde ∆tTotal é o tempo total de aquecimento e ∆tstep é o passo de tempo da
fase térmica I, como já dito anteriormente. Na Figura 3.1, tal loop é sinalizado pelas setas
normais, enquanto as setas em negrito sinalizam entrada de informação no loop, como
condição de contorno.
No acoplamento das fases eletromagnética e térmica I implementado neste trabalho, é
considerado que em pequenos passos de tempo o calor gerado pela fase eletromagnética
permanece aproximadamente constante, eliminando assim um processo iterativo muito
custoso computacionalmente em que a potência dissipada variaria durante o passo de tempo
da fase térmica transiente. Por essa razão, é necessário que o passo de tempo da fase térmica I
seja pequeno o suficiente para que não haja grandes alterações no campo de temperaturas e
consequentemente nas propriedades eletromagnéticas, que levaria à uma alteração da potência
dissipada.
De uma maneira mais resumida: durante a simulação do aquecimento do componente
mecânico as fases eletromagnética e térmica estão acopladas da seguinte maneira: a simulação
da fase eletromagnética fornece como resultado o calor gerado na superfície do componente
mecânico através do Efeito Joule que é inserido como força de volume na simulação da fase
térmica I. Esta, por sua vez, fornece como resultado o campo de temperaturas, que é utilizado
para atualizar as propriedades eletromagnéticas e define o tempo total da simulação do
aquecimento, uma vez que a simulação da fase eletromagnética é harmônica no tempo. O
tempo total de aquecimento é um dado definido pelo operador do sistema, assim como a
frequência e a amplitude da corrente elétrica do indutor.
O campo de temperaturas gerado ao final do loop contendo as fases eletromagnética e
térmica I, descrito na Figura 3.1, é utilizado como condição inicial para a simulação do tempo
44
de retardo. O tempo de retardo consiste no período de tempo que o componente mecânico fica
suspenso no ar a fim de se promover a difusão do calor gerado na superfície deste e, assim,
permitir que a camada superficial austenitizada cresça. Este tempo também é um dado
escolhido pelo operador do processo.
Sabe-se que para transformação de perlita em austenita necessita-se da absorção de
calor latente, porém este foi negligenciado em razão da magnitude do calor gerado por Efeito
Joule. Sabe-se também que tal transformação é dependente do tempo, mas a fim de se
diminuir o tempo total de processamento, a microestrutura no momento anterior ao tempo de
retardo e posterior ao desligamento do indutor é calculada apenas com base na temperatura
final de cada elemento finito.
Um segundo loop tem início após a determinação da microestrutura pós-aquecimento.
Neste loop estão acopladas a fase térmica II e a fase microestrutural I, que correspondem,
respectivamente à simulação do tempo de retardo e simulação da transformação
microestrutural durante o mesmo. Durante a simulação térmica do tempo de retardo foi
utilizada a mesma condição de contorno da fase térmica I, com exceção do calor sendo gerado
no componente mecânico, visto que as bobinas são desligadas durante o tempo de retardo.
A cada passo de tempo da fase térmica II são avaliadas as transformações
microestruturais, apenas no sentido de produção de austenita, uma vez que o resfriamento ao
ar é lento. Mais uma vez o calor latente absorvido foi negligenciado aqui, desta vez devido ao
pequeno intervalo de tempo que consiste o tempo de retardo, em geral um ou dois segundos, e
também a austenita gerada foi considerada somente como função da temperatura. Desta forma
os elementos finitos foram classificados como austenitizados (100% austenita) e não
austenetizados (0% austenita e 100% perlita). Devido às dimensões diminutas dos elementos
finitos na superfície do componente mecânico, da ordem de 0,3 mm, estas simplificações não
devem produzir grandes erros no perfil de dureza final.
Ao final deste segundo loop obtém-se o campo de temperaturas e o perfil
microestrutural imediatamente anterior ao resfriamento. O perfil microestrutural pré-
resfriamento tem uma grande importância, pois apenas as regiões que se transformaram em
austenita podem se transformar durante o resfriamento e produzir, por exemplo, martensita.
45
3.2 Resfriamento
Para a simulação computacional do resfriamento é feito um acoplamento da fase térmica
III e da fase microestrutural II. A fase térmica III consiste no resfriamento propriamente dito
do componente mecânico, na qual este é imerso em água após passar pelo aquecimento
eletromagnético e pelo tempo de retardo. Toda a simulação da etapa de resfriamento é feita
dentro do módulo térmico do Ansys, sendo necessário desabilitar os elementos das estruturas
pertencentes unicamente a simulação eletromagnética: ar, bobina e concentrador de fluxo.
O primeiro passo do loop de resfriamento é inserir as condições iniciais e de contorno do
problema, e como já explicado, as condições iniciais são o campo de temperaturas e o perfil
microestrutural obtidos após o tempo de retardo. Já as condições de contorno consistem na
introdução do coeficiente de película em função da temperatura, através de uma tabela. O
próximo passo do loop de resfriamento consiste na simulação da fase térmica III, a fim de se
obter um novo campo de temperaturas, após ∆tstep_resf, passo de tempo da fase térmica III.
Esse novo campo de temperaturas é utilizado em conjunto com o perfil microestrutural para
calcular no perfil estrutural e assim, atualizar as propriedades físicas (densidade, calor
específico e condutibilidade térmica), que são dependentes da microestrutura, de cada
elemento, após um intervalo de tempo ∆tstep_resf.
Apesar de haver liberação calor latente durante o processo físico da transformação
microestrutural e também haver alteração das propriedades de acordo com a microestrutura,
para diminuir o custo computacional, os mesmos são mantidos constantes durante ∆tstep_resf.
Dessa forma, este último deve ser pequeno o suficiente para que não ocorram erros
significativos na avaliação do campo de temperaturas, devido ao procedimento iterativo
sequencial utilizado.
O loop é realizado até o momento em que não haja nenhum elemento com temperatura
acima da temperatura de início da martensita, visto que nessa condição não é possível
formação de nenhuma outra fase além desta.
46
Figura 3.2: Estratégia de simulação da fase de resfriamento.
A segunda etapa consiste no pós-processamento dos resultados gerados pela simulação do
resfriamento para a obtenção da dureza final do componente mecânico. Alguns métodos
podem ser utilizados para descrever a dureza de cada fase e assim o perfil de dureza final do
componente mecânico, como: atribuir uma dureza fixa para cada fase, utilizar o diagrama
TTT e extrair a dureza para a porcentagem de microestrutura formada em uma determinada
temperatura ou ainda utilizar equações empíricas. Neste trabalho utilizou-se as equações de
Maynier e a atribuição de um valor fixo de dureza para cada fase.
3.3 Simulação multifísica integrada
Um dos objetivos deste trabalho é a implementação de um modelo integrado para
simulação multifísica do processo de têmpera. As duas seções anteriores trataram da
simulação das fases de aquecimento e resfriamento, que devem ser acopladas a fim de se
produzir o resultado final da simulação, que é o perfil de dureza.
A Figura 3.3 mostra como a simulação integrada foi estruturada, onde os quadros em
vermelho correspondem ao loop de aquecimento, os em laranja ao loop do tempo de retardo e
os em azul ao loop de resfriamento.
47
Figura 3.3: Esquema representando a simulação integrada do processo de têmpera por indução
eletromagnética.
48
4 RESULTADOS NÚMERICOS E ANÁLISE PARAMÉTRICA
Neste capítulo serão apresentados os resultados numéricos obtidos utilizando a
estratégia de simulação descrita nos capítulos anteriores. Apresentam-se inicialmente alguns
exemplos de validação das etapas de formulação, após, foi analisada a influência de alguns
parâmetros da simulação, do processo, do indutor e da peça para duas geometrias distintas da
peça a ser aquecida: um cilindro e um moente de um virabrequim.
A compreensão da influência dos diversos parâmetros dos processos de aquecimento e
têmpera é muito importante, pois tratam-se de processos complexos muito utilizados na
indústria e onde ainda são utilizados o métodos empíricos no setup dos mesmos. Dessa forma,
entender a influência de cada parâmetro através da aplicação de técnicas de simulação pode
reduzir o número de tentativas e ensaios experimentais, e consequentemente o custo e o tempo
total, para se obter o resultado esperado.
4.1 Avaliação numérica da modelagem da convecção natural
Segundo Cho (2012), durante o aquecimento do componente mecânico por indução
eletromagnética, o calor extraído da peça por convecção é pequeno quando comparado ao
calor perdido por radiação e, por isso, pode ser negligenciado. Para testar essa hipótese, será
calculada a taxa de calor perdido por unidade de comprimento para um cilindro horizontal por
convecção natural e comparada com a taxa de calor perdido por unidade de comprimento por
radiação. É sabido que existe um grande gradiente de temperaturas no cilindro durante o
aquecimento, mas para fim de simplificação, esse gradiente não será levado em conta e será
considerado que o cilindro possui temperatura homogênea para fim de comparação da
influência relativa dos efeitos de convecção e radiação.
O número de Nusselt médio para um cilindro horizontal de temperatura homogênea
pode ser calculado pela seguinte expressão (Cho, 2012):
49
Nud =
{
0,6 +0,387(Rad)
0.5
(1+(0,559
Pr)
916)
827
}
(4.1)
onde Rad é o número de Rayleigh e Pr é o número de Prandtl. O número de Rayleigh e o de
Prandtl podem ser calculados por:
Rad =gβ(Ts−Tar)d
3
ρarν (4.2)
Pr =ν
α (4.3)
onde g é a aceleração gravitacional, Ts a temperatura superficial do cilindro, Tar a temperatura
do ar envolta do cilindro, α a difusividade térmica do ar, β é o coeficiente de expansão
térmica do ar, d o diâmetro do cilindro, ρar a densidade do ar e ν a viscosidade cinemática do
ar. Assim a taxa de calor perdido por convecção natural por unidade de comprimento Qconv′ é
dada por:
Qconv′ =
Qconv
l= hcnπ(Ts − Tar) (4.4)
hcn =kNud
d (4.5)
Aqui hcn é o coeficiente convectivo, l o comprimento do cilindro e k é a condutividade
térmica do material.
Para a taxa de calor perdido por radiação, pode-se escrever:
Qrad′ =
Qrad
l= πdεσr(Ts
4 − T∞4 ) (4.6)
50
onde ε é a emissividade do material, σr é a constante de Stefan-Boltzmann e T∞ a
temperatura do ambiente. Para efeito de comparação foi considerado T∞ = Tar, g = 9,81m
s2,
ε = 0,8, d = 38 mm. O resultado obtido para Qrad′ e Qconv
′ é mostrado na Figura 4.1.
Observa-se na Figura 4.1 que para baixas temperaturas as perdas por convecção e radiação
são equivalentes, mas a partir de 300°C as diferenças são notáveis.
Durante a implementação computacional, a radiação foi modelada por um modelo de
convecção equivalente cujo coeficiente convectivo é função da temperatura, hrad.
Inicialmente fez-se a expansão binomial do termo da temperatura a quarta potência, obtendo-
se:
Qrad′ = Aεσr(Ts
4 − T∞4 ) = Aεσr(Ts
2 + T∞2 )(Ts + T∞)(Ts − T∞) (4.7)
onde A é a área da superfície do cilindro em contato com o ar. Assim, chamando de hrad o
coeficiente convectivo equivalente, pode-se escrever (Cajner et al., 2004; Incropera e Dewitt,
2008):
hrad = εσr(Ts2 + T∞
2 )(Ts + T∞) (4.8)
e portanto:
Qrad′ = Ahrad(Ts − T∞) (4.9)
que é o modelo convectivo equivalente usado na simulação, considerando o coeficiente hrad
variando com a temperatura Ts da superfície do cilindro. Tal modelo, implica em uma solução
de um problema térmico não linear, que foi resolvido usando-se um processo iterativo
implementado no Ansys.
51
Figura 4.1: Comparação entre a taxa de calor perdido por convecção e por radiação.
Como a região de interesse está sujeita a altas temperaturas, o calor perdido por
convecção natural será desprezado na simulação da fase de aquecimento.
Estes resultados mostram que o efeito da radiação são preponderantes na fase de
aquecimento, principalmente em altas temperaturas, o que justifica a não consideração do
efeito de convecção neste caso.
4.2 Teste numérico do modelo microestrutural
Para validar a fase de resfriamento e a implementação do modelo microestrutural
empregado foi simulada a têmpera de um cilindro de aço SAE 1080 com 38mm de diâmetro e
76mm de comprimento. Como condições iniciais foi utilizada uma temperatura uniforme de
850°C e 100% de austenita como microestrutura inicial, de acordo com o que foi considerado
por Carlone et al. (2010).
Foi considerado a têmpera por imersão em água a 20°C e o coeficiente convectivo para
tal condição foi retirado de (Woodard et al., 1999), como mostrado na Figura 4.2.
0
5000
10000
15000
20000
0 200 400 600 800 1000 1200
Tax
a d
e ca
lor
per
did
o p
or
un
idad
e
de
com
pri
men
to [
W/m
]
Temperatura da supefície do cilindro [°C]
Convecção
Radiação
52
Figura 4.2: Coeficiente convectivo para um cilindro de 38mm de diâmetro imerso em água.
Extraído de Woodard et al. (1999).
Os resultados obtidos desta simulação foram comparados aos publicados por Carlone et
al. (2010). Existem, porém, diferenças entre as duas simulações. Na presente simulação o
modelo da transformação austenita-martensita utilizado foi o de Koistinen-Marburguer, além
disso, o cilindro, com uma razão comprimento/diâmetro igual a 2, foi modelado como
infinitamente longo (Huiping et al., 2007) utilizando-se 95 elementos ao longo do raio.
Carlone et al, (2010) utilizou um modelo 2D com 10 elementos ao longo do raio e a
transformação austenita-martensita foi descrita pelo modelo de Yu.
Na comparação do histórico de temperaturas foram utilizados 3 pontos posicionados ao
longo do eixo x: A=0,20mm, B=8,10mm e C=18,96mm como ilustra a Figura 4.3. No artigo
de referência (Carlone et al., 2010) usou-se A=0mm, B=8,08mm e C=19,05mm. A Figura 4.3
mostra tais pontos e as condições de contorno aplicadas nesta simulação. A diferença na
posição dos pontos ocorreram devido a uma escolha diferente na malha de Elementos Finitos
utilizada.
As figuras 4.4 e 4.5 mostram, respectivamente, a distribuição microestrutural ao longo
do raio após 90s de têmpera e o histórico de temperaturas dos pontos A, B e C. É mostrado
também uma comparação com os resultados de Carlone et al. (2010). Apesar dos dois
modelos apresentarem diferenças entre si, houve uma boa concordância entre as duas
simulações, demonstrando que os modelos microestrutural e térmico acoplados tem um
0
5
10
15
20
0 200 400 600 800
Coef
icie
nte
Convec
tivo
[kW
/m²K
]
Temperatura [°C]
53
comportamento compatível com os resultados previamente publicados na literatura
especializada.
Figura 4.3: Condições de contorno e localização dos pontos A, B e C da seção axissimétrica
do cilindro utilizado.
Figura 4.4: Distribuição microestrutural de perlita e martensita depois de 90s de imersão em
água.
54
Na Figura 4.5, nota-se um pequeno aumento da temperatura nos pontos A e B durante o
resfriamento. Esse aquecimento se dá pela liberação de calor latente durante a transformação
de fase. A medida em que se aproxima da superfície, a taxa de calor extraído pela imersão em
água é maior, subjugando o efeito da liberação de calor latente na transformação de fase.
Figura 4.5: Histórico de temperaturas durante o resfriamento por imersão em água para os
pontos A,B e C.
Os resultados, tanto em fração de volume da martensita e perlita quanto em perfil de
temperatura, concordam com os resultados da literatura. Assim, é possível afirmar que a
metodologia implementada neste trabalho para a fase de resfriamento é adequada.
4.3 Influência das propriedades na fase de resfriamento
O estudo da influência das propriedades do material da peça a ser aquecida tem como
objetivo neste trabalho verificar a necessidade da utilização de propriedades como função da
microestrutura do material.
Como tais propriedades são difíceis de obter, ainda mais pelo espectro de temperaturas
envolvido, esta seção visa avaliar a diferença nos resultados da simulação numérica entre
utilizar propriedades em função da temperatura e utilizar propriedades mecânicas em função
55
da temperatura e da microestrutura. Isto é, pretende-se avaliar a relevância de se modelar as
variações microestruturais durante a simulação térmica.
Foi utilizada a malha apresentada na Figura 4.6 para simular a fase de resfriamento. O
objeto de estudo é um cilindro de aço SAE 1080 infinitamente longo, já que seu comprimento
é igual a duas vezes o seu diâmetro (Huiping et al., 2007), com 38mm de diâmetro. O cilindro
é resfriado em água a 20ºC e o coeficiente convectivo utilizado foi extraído de (Woodard et
al., 1999), mostrado na Figura 6. A dureza final foi calculada de acordo com a porcentagem
de cada fase, levando em conta uma dureza fixa para cada fase, 39 Rockwell C para perlita e
66 RC para martensita.
Figura 4.6: Malha utilizada para a simulação do cilindro.
Em um primeiro momento, foi analisada a dureza final obtida quando a temperatura
inicial do cilindro é homogênea e igual a 850ºC e microestrutura inicial de 100% austenita.
Foram realizados no total 4 casos: utilizando apenas uma das propriedades (condutibilidade
térmica, densidade e calor específico) comumente utilizadas do aço e as duas outras
propriedades dependentes da microestrutura, utilizando todas as propriedades dependentes da
microestrutura e nenhuma das propriedades dependente da microestrutura. Vale salientar que
em todos os casos as propriedades são dependentes da temperatura.
Figura 4.7: Dureza final de acordo com a propriedade utilizada como não sendo dependente
da microestrutura.
0
10
20
30
40
50
60
70
0 0,005 0,01 0,015 0,02
Dure
za [
RC
]
Distância ao longo do raio [mm]
controle
calor específico
condutividade
térmicaTodas
56
Através da Figura 4.7, pode-se notar que a propriedade que causa a maior influência na
dureza final quando utilizada como não dependente da microestrutura é a condutibilidade. As
curvas “Controle” e “Todas” denotam, respectivamente a simulação em que todas as
propriedades são dependentes da temperatura e a simulação em que todas as propriedades não
são dependentes da temperatura.
Tal comportamento pode ser explicado utilizando os gráficos de cada uma das
propriedades estudadas nesta seção. Nas figuras 4.8-4.10 são mostradas comparações entre as
propriedades comumente utilizadas e as propriedades dependentes da temperatura. Nota-se
que a propriedade não dependente da microestrutura aproxima-se da curva da propriedade da
perlita para a condutibilidade térmica e calor específico, a curva de densidade do aço SAE
1080 encontrada na literatura, além disso possui uma pequena variação entre as curvas de
austenita e perlita (razão máxima entre valor da densidade da austenita e da perlita de 1,02 na
temperatura de 0°C ou 273K), assim a influência desta propriedade não será estudada.
Ao se calcular a razão entre a propriedade não dependente da microestrutura
propriedade da austenita por volta de 550°C ou 823K, que é a temperatura que mais favorece
a transformação austenita-perlita, verificamos que esta razão é muito maior no caso da
condutibilidade térmica, como mostrado na tabela 4.1.
Tabela 4.1: Razão entre a propriedade não dependente da temperatura e a propriedade da
austenita em 550°C.
Condutibilidade Térmica Calor Específico
1,70 1,27
Como a microestrutura inicial é de 100% austenita e a curva de propriedade tende à
curva de propriedades da perlita, o resultado obtido é completamente diferente daquele com
todas as propriedades como função da microestrutura. Porém, um componente, com 19mm de
raio, totalmente austenitizado não é o que se obtém utilizando o aquecimento por indução
eletromagnética. Uma das vantagens da indução eletromagnética é aquecer prioritariamente a
superfície do componente mecânico, dessa forma, apenas uma pequena região do mesmo será
austenitizada.
57
Por esta razão, uma segunda simulação computacional foi realizada a fim de se verificar
a influência das propriedades como função da microestrutura. A malha e as condições de
contorno convectivas são as mesmas da simulação anterior, porém, nesta etapa o cilindro tem
temperatura inicial igual a 30°C e é aquecido até que apenas uma região se austenitize. Para o
aquecimento do cilindro procurou-se inserir uma fonte de calor semelhante a do aquecimento
por indução eletromagnética, ou seja, o elemento finito mais próximo à superfície recebeu
80% do calor volumétrico gerado total, o segundo mais próximo à superfície 16% e o terceiro
4%.
Figura 4.8: Densidade Aço SAE 1080.
Figura 4.9: Calor Específico do Aço SAE 1080
7300
7400
7500
7600
7700
7800
7900
8000
8100
0 500 1000 1500 2000
Den
sidad
e [K
g/m
³]
Temperatura [K]
Austenita
Perlita
0
100
200
300
400
500
600
700
800
900
0 500 1000 1500 2000
Cal
or
espec
ífic
o [
J/K
gK
]
Temperatura [K]
Calor Específico
Austenita
Perlita
58
Figura 4.10: Condutibilidade Térmica do Aço 1080
Diferentes valores de calor volumétrico foram utilizados a fim de se obter diferentes
espessuras de camada austenitizada. Os resultados de dureza final obtidos para 8,5mm,
5,1mm e 3,6mm de espessura de camada austenitizada estão mostrados nas figuras 4.11, 4.12
e 4.13, respectivamente. Foram simulados apenas dois casos para cada configuração, com
todas as propriedades dependentes da microestrutura e com todas as propriedades não
dependentes.
Figura 4.11: Dureza final para da simulação com 8,5mm de espessura austenitizada
0
10
20
30
40
50
60
0 500 1000 1500 2000
Conduti
bil
idad
e T
érm
ica
[W/m
K]
Temperatura [K]
Condutibilidade
TérmicaAustenita
Martensita
Perlita
0
10
20
30
40
50
60
70
0 0,005 0,01 0,015 0,02
Dure
za [
RC
]
Distância ao longo do raio [m]
Dependente
Não
Dependente
59
Figura 4.12: Dureza final para da simulação com 5,1mm de espessura austenitizada
Figura 4.13: Dureza final para da simulação com 3,6 mm de espessura austenitizada
Através as figuras 4.11-4.13, pode-se concluir que com a diminuição da camada
austenitizada há também uma diminuição da diferença entre as duas curvas. Isso ocorre, pois,
a propriedade não dependente da temperatura se assemelha a curva da propriedade da perlita.
Assim, como há uma diminuição da área austenitizada e um aumento da área com 100%
perlita a diferença entre as curvas também diminui.
0
10
20
30
40
50
60
70
80
0 0,005 0,01 0,015 0,02
Dure
za [
RC
]
Distância ao longo do raio [m]
Dependente
Não Dependente
0
10
20
30
40
50
60
70
80
0 0,005 0,01 0,015 0,02
Dure
za [
RC
]
Distância ao longo do raio [m]
Dependente
Não Dependente
60
Dessa maneira, como usualmente a profundidade da camada temperada não é maior do
que 5mm, será utilizado neste trabalho apenas as propriedades globais, não dependentes da
microestrutura, devido à pequena influência deste tipo de modelagem nestes casos.
4.4 Simulação completa do processo integrado para o caso de um cilindro
O terceiro resultado apresentado engloba as fases de aquecimento e resfriamento de
forma sequencialmente acoplada. O objetivo é mostrar uma sequência completa do
procedimento de simulação, aplicada a um caso simples de um cilindro, onde as malhas são
mais facilmente controladas. Foi simulada a têmpera de um cilindro de aço SAE 1080,
aquecido por indução eletromagnética e resfriado por imersão em água. O cilindro possui
18mm de diâmetro e 76mm de comprimento, mas devido à simetria, foi modelado apenas
38mm de comprimento, conforme mostrado na Figura 4.14. As dimensões da malha de ar que
envolve o cilindro são duas vezes maior longitudinalmente e duas vezes e meia maior
radialmente do que as dimensões do cilindro, para evitar a influência das condições de
contorno no problema (Drobenko et al., 2007).
Figura 4.14: Representação do modelo axissimétrico para o caso do cilindro.
61
O indutor foi modelado considerando 8 fios de cobre envolvendo o cilindro e colocado
no plano de simetria. Foi considerado que o ângulo entre as espiras é muito pequeno e por
isso pode ser negligenciado, sendo assim possível a modelagem do problema como
axissimétrico. O modelo pode ser visto na Figura 4.14, onde a área cinza representa o cilindro,
a azul o ar e as áreas laranjas os fios do indutor.
Figura 4.15: Expansão axissimétrica em ¾ do modelo, representando apenas o indutor e
o cilindro.
Na Figura 4.15 foi feita uma expansão axissimétrica em ¾ do modelo do cilindro e das
espiras, lembrando que ele é simétrico com relação ao plano 𝑥𝑧, mas que essa simetria não foi
expandida nesta representação. Além disso, para facilitar a visualização, a malha de ar não
está exposta.
A malha do cilindro utilizada nessa simulação possui 55 elementos ao longo do
comprimento e 28 ao longo do raio. Os elementos próximos à superfície do cilindro são mais
finos dos que os próximos ao eixo, como pode ser visualizado na Figura 4.15. Esse tipo de
refinamento é necessário devido à baixa profundidade de penetração da onda eletromagnética
62
e ao decréscimo da potência dissipada em uma razão de 𝑒−2 por profundidade de penetração, o
que conduz a camadas de baixa espessura onde o calor é concentrado, e por conseguinte, as
transformações metalúrgicas do material. Tipicamente a camada endurecida atinge poucos
milímetros.
Utilizou-se uma corrente de 1500 𝐴 e 𝑓=1000 𝐻𝑧 em cada um dos fios do indutor, que
possuem um raio de 1,2mm. O diâmetro do indutor é 46mm e a distância entre os fios do
mesmo é de 5mm. Por simplicidade, o calor gerado nos fios da espira por efeito joule foi
desprezado.
O cilindro, inicialmente a 20°C, é aquecido por 20s, utilizando-se um passo de tempo de
1s durante todo o aquecimento. Na Figura 4.16, é mostrado o campo de temperaturas, em
Kelvin, ao final do aquecimento.
Figura 4.16: Campo de temperaturas em Kelvin obtido após 20 segundos de aquecimento.
Para fins de simplificação e diminuição do tempo de processamento, o tempo de
transformação durante a fase de aquecimento foi negligenciada, dessa forma os elementos
com temperatura acima da temperatura de austenitização do aço 1080 (996 K) foram
considerados como 100% austenita, e os elementos com temperatura abaixo da temperatura de
austenitização como sendo 100% perlita. Pode-se observar na Figura 4.16 que apenas uma
63
parte do cilindro ultrapassa tal temperatura, ou seja, somente essa porção estará sujeita a
transformação microestrutural durante o resfriamento.
Durante a fase de aquecimento as propriedades eletromagnéticas (condutividade elétrica
e curva de magnetização) e as propriedades térmicas (calor específico e condutibilidade
térmica) do aço são funções da temperatura, conforme mostrado nas tabelas 4.3 a 4.7. É
importante salientar que acima da Temperatura de Curie, que neste caso vale Tcurie = 1043 K
a permeabilidade magnética do aço não é obtida a partir da curva de magnetização, em vez
disso, ela se torna constante μ = μ0 = 4π10−7 N
A2.
O calor perdido por convecção durante o período de aquecimento foi negligenciado por
ser considerado pequeno quando comparado à perda de calor por radiação. Já durante o
resfriamento, apenas o calor perdido por convecção é considerado, utilizando-se o coeficiente
convectivo mostrado na Figura 4.2, que corresponde ao caso de imersão em água.
Nesse modelo foi adotado um valor constante de dureza para cada fase. Foi adotado 39
Rockwell C para a perlita e 66 RC para a martensita e austenita retida (Woodard et al., 1999).
A dureza da bainita foi assumida como 43 RC, de acordo com a escala de dureza do diagrama
TTT do aço SAE 1080.
A dureza da austenita retida é a mesma da martensita, pois ao final da simulação, todos
os elementos estão em uma temperatura abaixo da temperatura de início da transformação
martensítica (aproximadamente 220°C). Assim, como a transformação austenita-martensita
depende apenas da temperatura e a temperatura de final de transformação martensítica é maior
que a temperatura ambiente para o aço em questão, toda a austenita retida irá se transformar
em martensita tão logo a temperatura do componente mecânico se aproxime da temperatura
ambiente. Por isso, a fim de se diminuir o tempo de processamento, não é necessário que a
simulação se estenda até o cilindro atingir temperaturas próximas à temperatura ambiente.
O perfil da microestrutura e de dureza no plano de simetria após 120s de resfriamento
estão mostrados nas figuras 4.17 e 4.18, respectivamente. Na Figura 4.19 está mostrado um
mapa de cores da dureza pós-têmpera.
Pode-se observar que houve um acréscimo de dureza apenas em aproximadamente
3,5mm de espessura. Percebe-se, também, uma porcentagem muita baixa de bainita formada,
podendo ser associada à forma da curva TTT, que não favorece esse tipo de transformação, e
64
à pequena região aquecida, que está sujeita a uma alta taxa de extração de calor, favorecendo
a formação de martensita.
Comparando a Figura 4.17 com a 4.18, notamos que apesar de termos apenas
aproximadamente 80% de martensita na superfície do cilindro, a dureza da mesma é de 66
RC. Isso acontece em razão de a simulação não ter sido realizada até o resfriamento do
cilindro a uma temperatura abaixo da temperatura de final de transformação martensítica do
aço 1080, restando, ainda, uma porcentagem de austenita retida, cuja dureza adotada é a
mesma da martensita.
Figura 4.17: Perfil microestrutural pós-têmpera para o caso do cilindro.
Figura 4.18: Perfil de dureza pós-dureza para o caso do cilindro.
65
Figura 4.19: Mapa de cores da dureza em RC pós-têmpera para o caso do cilindro.
Não foi possível a realização de toda a simulação apenas com as propriedades do aço 1080
devido à dificuldade de obtenção de tais dados. Todavia, entende-se que as propriedades usadas
são representativas deste tipo de problema, e portanto, os fenômenos principais do processo foram
resguardados neste modelo. Dessa forma, a curva de magnetização utilizada foi a do aço AISI
1045, as propriedades elétricas e térmicas utilizadas durante a fase de aquecimento são do aço
SAE 1080 e todas as propriedades térmicas dependentes da microestrutura usadas na fase de
resfiramento são do aço SAE 1080.
De forma geral, pode-se observar que o processo de simulação é representativo e que a
espessura da camada endurecida está coerente com os valores normalmente praticados pela
indústria.
4.5 Simulação da têmpera por indução em um moente de um virabrequim
Nas seções seguintes serão explorados os resultados de cada uma das fases da simulação
do processo de têmpera de um moente do virabrequim. Para tanto, foi realizada uma
simulação de um esquema simplificado de ½ do moente, considerando a axissimetria.
66
Figura 4.20: Objeto de estudo da simulação integrada.
Figura 4.21: Malha utilizada para a simulação do moente.
67
A Figura 4.20 mostra o objeto de estudo, que inclui o indutor, o concentrador de fluxo
magnético e ar que envolve o sistema de aquecimento e o moente. A malha deste último, que
foi utilizada em todas as simulações seguintes, está mostrada na Figura 4.21.
Para diminuir o custo computacional, foram utilizados elementos finitos triangulares
fora da área de interesse. Os elementos finitos quadrilaterais na área de interesse, região onde
existe aumento de dureza devido às transformações microestruturais, são da ordem de 0,3mm.
A Figura 4.22 mostra uma expansão axissimétrica do modelo 2D.
Figura 4.22: Expansão axissimétrica em ¾ do modelo, representando o indutor, o
concentrador de fluxo eletromagnético e o moente.
4.5.1 Condições iniciais e de contorno do problema e propriedades dos
materiais
Para a simulação da têmpera por indução eletromagnética no moente, considerou-se que
o componente mecânico estava inicialmente a uma temperatura inicial Ti = 20°C e foi
aquecido por 10s. Após o aquecimento e o tempo de retardo, este foi imerso em água e
resfriado até a temperatura ambiente. O passo de tempo utilizado durante a fase de
aquecimento (fase térmica I) foi de 0,25s, durante o tempo de retardo (fase térmica II) foi de
68
0,1s e durante o resfriamento (fase térmica III) foi de 0,05s. As figuras 4.23 e 4.24 mostram as
condições de contorno, dimensões do sistema e dados de entrada para a simulação
eletromagnética e térmica, respectivamente.
Figura 4.23: Condições de contorno e dados de entrada para a simulação eletromagnética.
Para avaliar a profundidade da camada endurecida no virabrequim foram definidas duas
direções A e B, como mostrado na Figura 4.25. A dureza é então avaliada nos elementos do
moente que fazem fronteira com a malha de ar na direção A e nos elementos mais próximos
ao eixo das abscissas na direção B. Esse procedimento foi utilizado pois a camada endurecida
não é uniforme. Além disso, a espessura da camada endurecida nessas regiões é um dos
parâmetros usados na indústria para analisar a qualidade do perfil da camada endurecida.
Neste trabalho, a camada endurecida foi considerada como sendo a região cuja dureza é
superior a 55 RC.
69
Figura 4.25: Direções A e B onde a profundidade da camada endurecida foi observada.
Figura 4.24: Condições de contorno para a simulação térmica.
70
Tabela 4.2: Dados de entrada para a simulação do moente.
As tabelas 4.3 a 4.7 apresentam as propriedades do aço SAE 1080 utilizadas durante a
simulação do processo de têmpera por indução. Todas as tabelas foram inseridas dentro do
código executado pelo programa Ansys, que utiliza interpolação linear para obter valores
entre os fornecidos na tabela (Ansys, 2007).
Tabela 4.3: Resistividade elétrica, calor específico e condutibilidade térmica do aço SAE
1080 utilizadas durante a simulação.
71
Tabela 4.4: Densidade do aço SAE 1080.
Os dados da Tabela 4.3 e 4.4 foram utilizados durante a etapa de aquecimento, ou seja,
nas fases térmicas I e II e na fase eletromagnética.
Já as tabelas 4.5 a 4.7 foram utilizadas na simulação do resfriamento do componente
mecânico. As tabelas 4.5 e 4.6 não apresentam os valores relativos à fase martensita, porém
durante a simulação as propriedades da martensita foram assumidas como iguais a da
austenita (Woodard et al., 1999).
Tabela 4.5: Densidade do aço SAE 1080 em função da microestrutura e da temperatura.
72
Tabela 4.6: Calor específico do aço SAE 1080 em função da microestrutura e da temperatura.
Tabela 4.7: Condutibilidade térmica do aço SAE 1080 em função da microestrutura e da
temperatura.
73
4.5.2 Influência do Concentrador de Fluxo Magnético
O concentrador de fluxo magnético trata-se de um material com alta permeabilidade
magnética que é adicionado ao indutor como forma de se controlar e se aumentar o fluxo
magnético na região de interesse. O material do concentrador de fluxo magnético também
deve apresentar uma boa condutividade térmica para não atrapalhar a dissipação do calor
gerado no indutor, porém não foi considerada dissipação de calor no indutor durante a
simulação como meio de simplificação do modelo.
Nesta seção será avaliada a influência da inserção de um concentrador de fluxo
eletromagnético ao redor do indutor na dureza final e na potência absorvida pela peça. O
concentrador de fluxo utilizado foi o Fluxtrol 100, cuja curva de saturação, que relaciona a
densidade de fluxo magnético com a intensidade de campo, é mostrada na Figura 4.26.
Figura 4.26: Curva de magnetização do Fluxtrol 100
Foram estudas três configurações: com concentrador de fluxo, sem concentrador de
fluxo e com concentrador de fluxo parcial. A Figura 4.27 apresenta a configuração com
concentrador de fluxo magnético e a configuração com concentrador parcial. A frequência de
oscilação e a amplitude da corrente utilizadas nesta seção foram 9000 Hz e 4000 A,
respectivamente.
0
0,2
0,4
0,6
0,8
1
1,2
1,4
1,6
0 5000 10000 15000 20000 25000 30000Den
sidad
e de
Flu
xo M
agnét
ico B
[T
]
Intensidade do campo magnético H [A/m]
74
Figura 4.27: Configuração do concentrador de fluxo: 1) com concentrador e 2) com
concentrador de fluxo parcial.
A Figura 4.28 apresenta as curvas de absorção de potência durante o aquecimento
eletromagnético para os três casos estudados nessa seção. Observando-se a Figura 4.28, nota-
se que ao utilizar o concentrador de fluxo eletromagnético na sua forma completa, há uma
maior absorção da energia eletromagnética e assim, um maior aquecimento do moente. Pode-
se observar também que a potência absorvida é ainda menor no caso da não utilização do
concentrador de fluxo eletromagnético, evidenciando assim a importância do mesmo no
aquecimento do moente para o melhor aproveitamento da energia fornecida pelo indutor.
Figura 4.28: Comparação entre a potência absorvida durante o aquecimento para as três
configurações estudadas.
0
20
40
60
80
100
120
0 2 4 6 8 10 12
Potê
nci
a A
bso
rvid
a [K
W]
Tempo [s]
Com concentrador de fluxo - 1
Concentrador de fluxo parcial - 2
Sem concentrador de fluxo - 3
75
Outro fator que deve ser levado em conta na utilização do concentrador de fluxo
eletromagnético é a influência do mesmo na camada endurecida pós-têmpera. A Figura 4.29
foi gerada a partir dos resultados do perfil de dureza final para cada um dos três casos citados
acima, a fim de se analisar este fator. Nota-se a formação de uma camada endurecida pós-
têmpera mais uniforme no caso do concentrador de fluxo magnético completo quando
comparada aos outros dois casos, sendo o caso mais irregular o que não utiliza o
concentrador.
Sabe-se a partir da experiência da indústria automotiva que uma das regiões mais
críticas para a falha por fadiga do virabrequim é no raio de curvatura entre dois moentes
consecutivos. Baseando-se nos resultados obtidos nessa seção, pode-se verificar a importância
da utilização do concentrador de fluxo eletromagnético também para o perfil de dureza pós-
têmpera. Pois, uma vez que o mesmo não é utilizado ou é utilizado apenas parcialmente, para
os casos analisados, houve uma diminuição da camada endurecida na região crítica,
favorecendo a falha por fadiga na mesma. Esta redução da camada endurecida é mais evidente
para o caso sem o concentrador mostrado na Figura 4.29.
Figura 4.29: Dureza final em Rockwell C para as configurações: a) com concentrador de
fluxo, b) com concentrador parcial e c) sem concentrador.
4.5.3 Influência da frequência de oscilação da corrente do indutor
A frequência de oscilação da corrente do indutor está relacionada à potência absorvida
pelo componente ferromagnético e à profundidade de penetração do campo eletromagnético
76
no mesmo. Esta última, representa a profundidade na qual a corrente elétrica gerada no
componente mecânico equivale a e−1, aproximadamente 36,8%, da corrente gerada na
superfície do componente mecânico. Como a maior parcela da potência dissipada pelo campo
eletromagnético no componente ferromagnético é devido ao Efeito Joule, a potência dissipada
é proporcional ao quadrado da corrente elétrica, assim a potência dissipada em uma
profundidade igual a δp é ainda menor, equivale e−2, aproximadamente 13,5%, da potência
dissipada na superfície do componente. Em bons condutores elétricos a profundidade de
penetração pode ser aproximada por:
δp = √2
ωμσ (4.10)
onde ω é a frequência de oscilação do campo eletromagnético, μ a permeabilidade magnética
do meio e σ a condutividade elétrica.
Para avaliar a influência da frequência de oscilação da corrente do indutor no perfil de
dureza e na potência absorvida no componente mecânico foram utilizados os seguintes input:
amplitude de corrente: 4200A, tempo de retardo: 1s. As figuras 4.30 e 4.31 abaixo mostram a
potência absorvida em cada passo de tempo do aquecimento (fase térmica I) e a dureza pós
têmpera para diferentes valores de frequência. Conforme a direções A e B mostrados na
Figura 4.25.
A potência absorvida possui um pico no início do aquecimento devido a
comportamento dos materiais ferromagnéticos, no qual a permissividade magnética relativa
diminui com o aumento da temperatura, chegando a 1, na Temperatura de Curie do material
do moente, que neste caso vale Tcurie = 770°C.
Observa-se também, que na direção A até aproximadamente ω = 5500Hz, não houve
camada endurecida, o que evidencia que para maiores frequências uma melhor distribuição do
calor gerado é obtida.
Notou-se também que a distribuição do calor gerado no aquecimento alterou-se de
acordo com a frequência de oscilação da corrente, alterando o perfil de dureza final do
moente.
77
Figura 4.30: Potência absorvida em cada passo de tempo para diferentes valores de
frequência.
Figura 4.31: Profundidade da camada endurecida nas direções A e B após a têmpera para
diferentes valores de frequência.
0
20
40
60
80
100
120
140
0 2 4 6 8 10 12
Potê
nci
a A
bso
rvid
a [k
W]
Tempo de Aquecimento [s]
3000 Hz 4000 Hz
5000 Hz 6000 Hz
7000 Hz 8000 Hz
9000 Hz
0
1
2
3
4
5
6
7
8
0 2000 4000 6000 8000 10000
Esp
essu
ra d
a C
amad
a en
dure
cida
[mm
]
Frequência [Hz]
Direção A
Direção B
A
B
78
4.5.4 Influência da amplitude da corrente no indutor
A amplitude da corrente no indutor é diretamente proporcional a amplitude do campo
magnético gerado pela mesma. Porém, como o material ferromagnético possui
permeabilidade magnética não linear, existe uma relação não linear entre a potência absorvida
pelo componente mecânico e a amplitude de corrente do indutor.
Para avaliar a influência da amplitude da corrente do indutor no perfil de dureza e na
potência absorvida no componente mecânico foram utilizados os seguintes dados de entrada:
frequência de oscilação da corrente: 9000Hz, tempo de retardo: 1s, foi considerado o caso
com concentrador integral. As figuras 4.32 e 4.33 abaixo mostram a potência absorvida em
cada passo de tempo do aquecimento (fase térmica I) e a dureza pós têmpera para diferentes
valores de amplitude de corrente.
Assim como na seção anterior, observou-se um pico maior potência absorvida no
início do aquecimento, porém a dureza final encontrada possui uma relação aproximadamente
linear com a amplitude da corrente do indutor.
Figura 4.32: Potência absorvida em cada passo de tempo para diferentes amplitudes de
corrente.
0
20
40
60
80
100
120
140
160
180
0 2 4 6 8 10 12
Potê
nci
a A
bso
rvid
a [k
W]
Tempo de Aquecimento [s]
3400 A 3600 A
3800 A 4000 A
4200 A 4400 A
4600 A 4800 A
79
Figura 4.33: Profundidade da camada endurecida nas direções A e B após a têmpera para
diferentes valores amplitude de corrente.
4.5.6 Influência da utilização das equações de Maynier para cálculo da
dureza
Na presente seção foi analisada a influência no perfil de dureza do virabrequim ao utilizar
as equações de Maynier para o cálculo da dureza das fases microestruturais ao invés de
utilizar uma dureza fixa para cada fase. As equações propostas por Maynier (4.11 a 4.13),
apresentadas abaixo, são função da taxa de resfriamento e da composição do aço.
HVp = 42 + 223C + 53Si + 30Mn + 12.6Ni + 7Cr + 1Mo
+ log10(vr)(10 − 19Si + 4Ni + 8Cr + 130V) (4.11)
HVb = 323 + 185C + 330Si + 153Mn + 65Ni + 144Cr + 191Mo + log10(vr)(89 + 53C −
55Si − 22Mn − 10Ni − 20Cr − 33Mo) (4.12)
HVm = 127 + 949C + 27Si + 11Mn + 8Ni + 16Cr + 21 log10(vr) (4.13)
0
1
2
3
4
5
6
7
8
3000 3500 4000 4500 5000
Esp
essu
ra d
a ca
mad
a en
dure
cida[
mm
]
Corrente [A]
Direção A
Direção B A
B
80
Nas equações acima, C é a porcentagem de carbono presente no aço, Si de silício, Mn de
manganês, Ni de níquel, Cr de cromo, Mo de molibdênio, V de vanádio, vr é a taxa de
resfriamento em °C/h, HVp, HVb e HVm são a dureza Vickers da perlita e da ferrita, bainita e
martensita, respectivamente.
Figura 4.34: Comparação entre o mapa de dureza pós tempera para os casos: a) utilização de
um valor de dureza constante para cada fase b) utilização das equações de Maynier para
determinação da dureza de cada fase.
As equações de Maynier podem ser utilizadas para o cálculo da dureza das regiões que
foram austenitizadas durante o aquecimento e resfriadas durante a tempera, desta forma, os
elementos que não sofreram austenitização tiveram sua dureza fixada em 400 Vickers.
Nesta simulação foram utilizados os mesmos dados de entrada da seção 5.2, ou seja,
frequência de oscilação e a amplitude da corrente no indutor de 9000 Hz e 4000 A,
respectivamente. Além disso, foram rodados os casos: com o concentrador de fluxo, com o
concentrador de fluxo parcial e sem o concentrador de fluxo eletromagnético. Após cada
81
passo de tempo durante a fase de resfriamento (fase térmica III) foram utilizadas as equações
de Maynier para calcular a dureza da fração da nova fase formada. Assim, para o cálculo da
dureza de cada elemento foi utilizada também a regra das misturas, ou seja, a dureza total de
cada ponto é composta pela dureza de cada fase formada ponderada pela própria porcentagem
de cada fase neste determinado ponto.
A Figura 4.34 mostra uma comparação qualitativa entre os dois modelos, dureza fixa de
cada fase e dureza calculada pelas equações de Maynier, apresentados na primeira e segunda
linha de imagens, respectivamente. Pode-se extrair que para os casos apresentados, utilizando
as propriedades do aço 1080 SAE, não houveram modificações significativas no resultado
final, quando se analisa o gradiente de dureza, que justifique o uso das equações de Maynier.
82
5 CONCLUSÕES
Neste trabalho foi apresentada uma metodologia de simulação multifísica do processo de
têmpera por indução eletromagnética aplicados a peças cilíndricas, utilizando aquecimento
por indução eletromagnética cujo resultado final é o perfil de dureza obtido após a têmpera.
As soluções dos problemas eletromagnético e térmico foram obtidas por meio do software
comercial Ansys. A modelagem microestrutural foi adicionada ao software Ansys usando-se a
programação paramétrica através da linguagem APDL. A integração dos modelos foi feita
usando-se um esquema iterativo com cálculo não linear em cada etapa, com acoplamento
sequencial no tempo. Não foram encontrados na literatura trabalhos que façam a modelagem
sequencialmente acoplada mostrada nesta dissertação.
Devido à característica do problema físico eletromagnético, das propriedades do aço e
da frequência relativamente alta, é necessário um grande refinamento da malha de elementos
finitos próximo à superfície do cilindro. Esse refinamento, aliado às várias etapas de
simulação e a necessidade de um pequeno passo de tempo durante o resfriamento, faz com
que o procedimento de simulação dure algumas horas, evidenciando-se um ponto que pode
ser melhor estudado e aprimorado.
Para o cálculo das transformações microestruturais foi criada uma rotina, executada
dentro do programa Ansys, contendo os modelos de JMAK para as transformações
difusionais, austenita-perlita e austeniata-bainita, e de Koistinen-Marburguer para a
transformação do tipo não difusional, austenita-martensita. A implementação do modelo
microestrutural foi comparada com Carlone et al. (2010) e, apesar de terem sidos utilizados
modelos para a transformação austenita-martensita diferentes, obteve-se uma boa
concordância entre os resultados.
Estas aplicações geraram duas publicações em congressos específicos da área
automotiva:
De Paula, P.; Pavanello, R.; Su, W.; Rodrigues, A. Multiphysics Simulation of Quenching
Process Coupling Electromagnetic, Thermal and Microstructural Analysis. SAE International,
2014, 2014-36-0425.
De Paula, P.; Pavanello, R.; Su, W.; Rodrigues, A. Sensitivity analysis of multiparameter
numerical model of quenching process using electromagnetic induction heating. SAE
International, 2015, 2015-36-0391.
83
Como a simulação envolve um grande gradiente de temperaturas e transformações
microestruturais, uma dificuldade encontrada foi a escassez de dados disponíveis na literatura
das propriedades dos materiais para fase microesturtural e mesmo propriedades globais para
toda a faixa de temperatura envolvida no processo. Uma saída para este problema foi utilizar
propriedades de aços similares em composição química ao aço SAE 1080.
Além disso, outra dificuldade encontrada foi a falta de resultados experimentais
disponíveis na literatura para a validação de cada fase do modelo e da simulação integrada do
processo de têmpera por indução.
A simulação como um todo, se validada experimentalmente e aplicada a geometrias
mais complexas, pode representar um ganho considerável de tempo e de recursos no projeto
de componentes que deverão passar pelo processo de têmpera por indução eletromagnética.
5.1 Sugestões de continuidade
O presente trabalho é amplo, envolvendo varias áreas da mecânica dos materiais,
térmica e eletromagnetismo. Assim varias frentes de trabalho podem ser mencionadas como
sugestão de trabalhos futuros:
avaliação experimental em casos simples de corpos cilíndricos visando a
qualificação do modelo para uso industrial;
ampliar o estudo sobre a modelagem eletromagnética na otimização da forma e da
topologia do concentrador de fluxo e dos indutores;
realizar a análise de sensibilidade de outros parâmetros do sistema, tais como
posição relativa do indutor e dos concentradores, meio no qual a peça é resfriada,
duração do tempo de retardo, forma e topologia da peça a ser temperada, entre
outros;
construir um pré e pós processador para facilitar o uso do programa desenvolvido e
permitir aos usuários dos equipamentos usar o modelo para auxiliar a regulagem
do processo para novas peças;
redução da energia consumida e do tempo total de aquecimento;
otimização da profundidade, perfil e dureza da camada endurecida;
implementar um modelo que para descrever as tensões residuais geradas durante o
processo;
84
avaliar a viabilidade da implementação de um modelo 3D para estudo da
influência de geometrias complexas, tais como peças com furos e geometrias não
axi simétricas.
85
REFERÊNCIAS
ANSYS, 11.0 Documentation, SAS IP, Canonsburg, 2007.
Atkinson C., Akbay T., Reed R. Theory for reaustenitisation from ferrite/cementite mixtures
in Fe-C-X steels. Acta metall. mater., v. 43, p. 2013-2031, 1995.
Avrami M. Kinetics of phase change I, general theory. Journal of Chemical Physics, v.7,
1939.
Avrami M. Kinetics of phase change II, transformation-time relations for random distribution
of nuclei. Journal of Chemical Physics, v.8, 1939.
Avrami M. Kinetics of phase change III, granulation, phase change, and microestructure.
Journal of Chemical Physics, v.9, 1940.
Babu K., Kumar T., Mathematical modedeling of surfasse heat flux during quenching. The
Minerals, Metals & Materials Society and ASM International, 2009.
Barka N., Bocher P., Brousseau J., Galopin M., Sundararajan S. Modeling and sensitivity
study of the induction hardening process. Advanced Materials Research, v. 15-17, p.
525-530, 2007.
Bodart O., Boureau A., Touzani R. Numerical investigation of optimal control of induction
heating process. Applied Mathematical Modelling, v. 25, p. 697-712, 2001.
Cajner F., Smoljan B., Landek D. Computer simulation of induction hardening. Journal of
Materials Processing Technology, v.157–158, p. 55–60, 2004.
Carlone P., Palazzo G., Pasquino R. Finite element analysis of the steel quenching process:
temperature field and solid-solid phase change. Computers and Mathematics with
Applications, v. 59, p. 585-594¸ 2010.
Çetinel H. Toparli M., Özsoyeller L. A Finite Element based prediction of the microstructural
evolution of steels subjected to the Tempcore process. Mechanics of Materialsv. 32, p.
339-347, 2000.
86
Cheng H.; Xie J.; Li J. Determination of surface heat-transfer coefficients of steel cylinder
with phase transformation during gas quenching with high pressures. Computational
Materials Science, v. 29, p. 453–458, 2004.
Cheng H.; Zhang S.; Wang H., Li J. Finite element analysis of temperature field with phase
transformation and non-linear surface heat-transfer coefficient during quenching.
Applied Mathematics and Mechanics, v. 19, nº 1, 1998.
Cho K. Coupled electro-magneto-thermal model for induction heating process of a moving
billet. International Journal of Thermal Science, v. 60, p. 195-204, 2012.
De Paula, P.; Pavanello, R.; Su, W.; Rodrigues, A. Multiphysics Simulation of Quenching
Process Coupling Electromagnetic, Thermal and Microstructural Analysis. SAE
International, 2014, 2014-36-0425.
De Paula, P.; Pavanello, R.; Su, W.; Rodrigues, A. Sensitivity analysis of multiparameter
numerical model of quenching process using electromagnetic induction heating. SAE
International, 2015, 2015-36-0391.
Di Luozzo N., Fontana M., Arcondo B. Modelling of induction heating of carbon steel tubes
Mathematical analysis, numerical simulation and validation. Journal of Alloys and
Compounds, v. 536, p. 564-568, 2012.
Doane D. Application of hardenability concepts in heat treatment steel. J. Heat Treatment,
v.1, p. 5-30, 1979.
Drobenko B., Hachkevych O., Kuornyts’kyi T. A mathematical simulation of high
temperature induction heating of electroconductive solids. International Journal of Heat
and Mass Transfer, v. 50, p. 616-624, 2007.
Grum J. A review of the influence of grinding conditions on residual stress after induction
surface hardening and grinding. Journal of Materials Processing Technology, v. 114, p.
212-226, 2001.
Grum J.; Ferlan D. Residual Internal Stresses after Induction Hardening and Grinding. Heat
Treating: Proceedings of the 17th Conference, 1998.
87
Guo X.L.; Wang L.; Green M.A. Coupled transient thermal and electromagnetic finite
element analysis of quench in MICE coupling magnet. Cryogenics, v.52, p. 420–427,
2012.
Hömberg D. A numerical simulation of the Jominy End-Quench test. Acta Mater, v. 44, p.
4375-4385, 1996.
Huiping L., Guoqun Z., Shanting N., Chuanzhen H. FEM simulation of quenching process
and experimental verification of simulation results. Materials Science and Engineering
A, p. 705-714, 2007.
Incropera F. P., Dewitt D. P. Fundamentos de transferência de calor e de massa. 6. ed LTC,
2008.
Jang J., Chiu Y. Numerical and experimental thermal analysis for a metallic hollow cylinder
subjected to step-wise electro-magnetic induction heating. Applied Thermal
Engineering, v. 27, p. 1883-1894, 2007.
Jacot A., Rappaz M., Reed R. Modelling of reaustenitization from the perlite structure in
steel. Acta mater, v.46, p. 3949-3962, 1998.
Jacot A., Rappaz M. A combined model for the description of austenitization, homogenization
and grain growth in hypoeutectoid Fe-C steels during heating. Acta mater, v. 47, p.
1645-1651, 1999.
Kang S., Im Y. Three-dimensional Finite Element Analysis of the quenching process of a
plain carbon steel with phase transformation. Metallurgical and Materials Transactions
A, v. 36A, p. 2315-2325, 2005.
Kapturkiewicz W., Fras E., Burbelko A. Computer simulation of the austenitizing process
in cast iron with pearlitic matrix. Materials Science and Engineering A, v. 413, p. 352-357,
2005.
Kianezhad M., Sajjadi S., Improvement of Quench Factor Analysis in phase and hardness
prediction of a quenched steel. Metallurgical and Materials Transactions A, v. 44A, p.
2053-2059, 2013.
Kranjc M., Zupnic A., Miklavcic D., Jarm T. Numerical analysis and thermographic
investigation of induction heating. International Journal of Heat and Mass Transfer, v.
53, p. 3585-3591, 2010.
88
Kristoffersen H.; Vomacka P. Influence of process parameters for induction hardening on
residual stresses. Materials and Design, v.22, p. 637-644, 2001.
Lee S., Pavlina E., Tyne C. Kinetics modeling of austenite decomposition for an end-
quenched 1045 steel. Materials Science and Engineering A, v. 527, p. 3186-3194, 2010.
Li F., Li X., Zhu T., Rong Y. Numerical Simulation of the moving inductor heating process
with magnetic flux concentrator. Advances in Mechanical Engineering, v. 2013, 2013.
Lusk M., Jou H. On Rule of Additivity in Phase Transformation Kinetics. Metallurgical and
Materials Transactions A, v. 287A, p. 287-291, 1997.
Magnabosco I., Ferro P., Tiziani A., Bonollo F. Induction heat treatment of a ISO C45 steel
bar: Experimental and numerical analysis. Computational Materials Science, v. 35, p.
98-106, 2006.
Palin-Luc T., Coupard D., Dumas C., Bristiel P. Simulation of multiaxial fatigue strength of
steel component treated by surface induction hardening and comparison with
experimental results. International Journal of Fatigue, v. 33, p. 1040-1047, 2011.
Reti T., Fried Z., Felde I. Computer simulation of steel quenching process using a multi-phase
transformation model. Computational Materials Science, v. 22, p. 261-278, 2001.
Ruan Y. Thermal stress evolution in continuously quenched circular bars. Int. J. Solids
Structures, v. 34, No. 28, p. 364-3656, 1997.
Srinivasan V., Moon K., Greif D., Wang D., Kim M. Numerical Simulation of immersion
quenching process of an engine cylinder head. Applied Mathematical Modelling, v. 34,
p. 2111-2128, 2010.
Toparli M.; Sahin S.; Ozkaya E.; Sasaki S. Residual thermal stress analysis in cylindrical steel
bars using finite element method and artificial neural network. Computers and
Structures, v. 80, p. 1763–1770, 2002.
Woodard P., Chandrasekar S., Yang H. Analysis of temperature and microstructure in the
quenching of steel cylinders. Metallurgical and Materials Transactions B, v.30B, p.
815-822, 1999.
89
Yang B., Hattiangadi A., Li W., Zhou G., McGreevy T. Simulation of steel microstructure
evolution during induction heating. Materials Science and Engineering A, v. 537, p.
2978-2984, 2010.
Zabett A., Azghandi S. Simulation of induction tempering process of carbon steel using finite
element method. Materials and Design, v. 36, p. 415-420, 2011.
90
ANEXO A – Artigo publicado no Congresso SAE 2014
2014-36-0425
Multiphysics Simulation of Quenching Process of a SAE 1080 Steel Cylinder, Coupling Electromagnetic, Thermal and Microstructural Analysis.
Pedro de Paula(1), Renato Pavanello(1), Wiliam Su(2), Alex Rodrigues(2)
(1) Department of Computational Mechanics, State University of Campinas,
(2) Department of Research and Development Engineering, ThyssenKrupp Metalúrgica Campo Limpo Ltda.
Copyright © 2014 SAE International
Abstract
Mechanical components, such as parts of internal combustion engine, subject to cyclic loads can be submitted to quenching process in order to improve mechanical properties preventing fatigue failures in service. It is important that such components, due to quenching process, get a high hardness surface layer, increasing the resistance to fatigue, and a tenacious core, with a high capacity of absorbing impacts.
In this paper, a multiphysics simulation method of quenching process using Finite Element Method is presented. The proposed simulation method include two stages: heating and cooling. In the first stage, the mechanical component, initially at ambient temperature, is heated by electromagnetic induction to a temperature above the steel austenitization. In the second one, the component is cooled by liquid immersion. The resulting microstructure is calculated using the Johnson–Mehl–Avrami–Kolmogorov model and Sheil’s additive rule for diffusional transformation, while austenite-martensite transformation is calculated by Koistinen-Marburguer equation.
The proposed method takes into account the variation of the material thermal properties as a function of temperature and microstructure, while the material electromagnetic properties are a function of temperature and strength of the electromagnetic field (magnetic permeability). As a result, the distribution of microstructure and hardness profile after quenching is obtained for a typical carbon steel, SAE 1080, for a mechanical component application.
Introduction
Mechanical components submitted to cyclic loadings are subjected to fatigue failure. One way to prevent such failures is submitting such components a quenching process, which basically consists in heating a component to above the austenitizing temperature of steel and then cooling it quickly.
After the quenching, it is desirable that the component presents a high hardness surface layer, with martensite microstructural configuration. The high hardness of the surface layer is important for the prevention of crack initiation in fatigue failure procedures. Another advantage of high surface hardness is increased wear resistance, on the other side, it is desirable that the mechanical component submitted to quenching process retain a tenacious core, thus providing impact resistance.
An efficient method used in the heating phase at the process is the electromagnetic induction. This type of heating has some advantages such as: greater process control, possibility to heat only one region and reduce the environmental impact.
Induction heating consists in approximating the mechanical component to an inductor traversed by a high frequency current. This current also produces electromagnetic waves of high frequency, which struck the component. The variation of the electromagnetic field at a given point generates electrical currents in order to counteract this variation of the field. As the penetration of the magnetic field is small in good conductive materials, current is generated only on the surface of the component. These currents dissipate heat by Joule effect, responsible for the heating.
The simulation method proposed in this paper is divided in two stages: heating and cooling. In both the finite element method is used. In the electromagnetic heating, the influence of the mesh size of the air surrounding the cylinder is studied, and it was found that this should be two times lengthwise and two and a half times greater than the radial dimensions of the cylinder being heated [1]. Furthermore, it is necessary the use of a finer mesh on the surface of the mechanical component containing from 3 to 10 elements in the region of heat generation [2].
During the electromagnetic heating, it was observed that for temperatures higher than 300°C and 600°C, if the variation of
91
electrical conductivity and magnetic permeability, respectively, are not considered, it results in considerable errors in the field of temperatures [1]. Among the properties like electrical resistivity, specific heat, thermal conductivity and magnetic permeability through a sensitivity study, that the first two are the ones that have more influence in the depth of the hardened layer [3].
The microstructural transformation, for simplicity, will be taken into account only during the cooling phase. At the end of the heating phase, the elements that present temperature above the austenitizing one will be considered as 100% austenite, those which do not achieve such temperature are going to be considered as 100% perlite.
The latent heating is being taking into account in the most part of the work and, despite history temperature is more intensely influenced near to the cylinder symmetry axis, to disregard it cause an important alteration in the superficial microstructural resulting [4]. Another heat source during the quenching, is which released by the deformations, however its represents only 3% of the total heat released and can be neglected [5].
The principal contribution of this paper is the implementation and testing a multiphysics procedure to perform a integrated simulation of the quenching process of mechanical components considering electromagnetic, thermal and microstructural analysis.
Modeling
Electromagnetic Formulation
The proposed integrated multiphysics simulation include three principal phases: electromagnetic, thermal and microstructural coupled.
The electromagnetic problem in consideration is governed by the Maxwell’s equations:
∇x𝐇 = 𝐉 +𝝏𝑫
𝝏𝒕 (1)
∇x𝐄 = − 𝝏𝑩
𝝏𝒕 (2)
∇.𝑩 = 0 (3)
∇.𝑫 = 𝜌𝑒 (4)
where 𝐇 is the magnetic field strength, 𝐄 is the electric field
intensity, 𝐉 is the conduction current density, 𝑫 is the electric
flux density, 𝜌𝑒 is the electric charge density. In good conductor materials and considering harmonic magneto-quasi-static problem, the displacement current is negligible compared with conduction current, so we can write the Maxwell’s equations as follows:
∇x𝐇 = 𝐉 (5)
∇x𝐄 = − 𝝏𝑩
𝝏𝒕 (6)
∇.𝑩 = 0 (7)
The constitutive equations described below supplement equations (5)-(7) and describes the behavior of the electromagnetics materials.
𝑱 = 𝜎𝑬 (8)
𝑩 = 𝜇𝑯 (9)
Here, 𝜎 is the electric conductivity and 𝜇 is magnetic permeability and 𝑩 is the magnetic flux density. In the present
model, 𝜇 is function of the temperature and of the magnetic
field intensity, so it was used the following approximation to describe the variation of the 𝑩-𝐇 curve with the temperature:
𝑩 = 𝜇(𝐻, 𝑇)𝑯 = 𝜇20(𝐻) [1 − (𝑇
𝑇𝑐𝑢𝑟𝑖𝑒)6]𝑯 (10)
Where 𝜇20 is obtained from the steel 𝑩-𝐇 curve at 20°C and
𝑇𝑐𝑢𝑟𝑖𝑒 is the Curie temperature of the material, in which the
steel lose its ferromagnetic properties, i.e., the steel magnetic permeability above Curie point is considered equal to the free space magnetic permeability.
Despite the electromagnetic field generated in the coil be sinusoidal the non-linearity of the magnetic permeability of the ferromagnetic material, made the field inside them non-sinusoidal. However, for power loss time-average problem, we can substitute the ferromagnetic material by a fictitious material, based on energy equivalence. Considering the problem as a time harmonic it is possible to reduce the computational cost and maintaining a good accuracy [6]. In this case, the magnetic field strength is approximated by:
1
2∫ 𝐻𝑚𝑑𝐵𝑒𝑓𝑓𝐵𝑒𝑓𝑓0
= 4
𝑃∫ (∫ 𝐻𝑚 sin(𝜔𝑡) 𝑑𝐵
𝐵
0)
𝑃
40
𝑑𝑡
(11)
where 𝐵𝑒𝑓𝑓 is the effective magnetic flux density, 𝑃 the time
period, 𝐻𝑚 the peak value of magnetic field and 𝜔 the angular velocity. The output off electromagnetic phase is:
𝑄 =1
2𝑅𝑒(𝑬𝑇𝑱∗) (12)
where 𝑄 is the volumetric heat power generate based on the
Joule Effect. The thermal simulation of the heating phase is based in the Fourier's equations as follows:
𝜌𝑐𝜕𝑇
𝜕𝑡= ∇ . (𝑘∇𝑇) + 𝑄 (13)
where 𝜌 is the mass density, 𝑘 is the thermal conductivity and 𝑐 is the specific heat.
The algorithm of the heating phase is shown in Figure 1. The temperature field obtained from the Thermal Phase I, at
92
heating stage, is used as initial conditional on Thermal Phase II, at cooling stage shown in Figure 2.
In the heating stage the purpose is evaluate the volumetric heat generated by the Joule Effect on the electromagnetic phase and it is an input to thermal phase I. This volumetric heat remain constant during a thermal time step. After the thermal phase I is done, the temperature filed is updated, as well the properties, and another electromagnetic simulation is done, as shown in Figure 1. These sequence runs until it got the total time heating desired.
In the cooling stage, Figure 2, the final temperature field obtained from the heating step is used as initial condition. The temperature filed is evaluated for each time step and the microstructural phase percentage, F, is updated using the microstructural model. So, the material properties and the heat released by the microstructural transformation is updated and another thermal phase II is performed. These sequence runs until the total time of cooling phase is reached.
Figure 1. Simulation strategy of the heating stage.
Figure 2. Simulation strategy of the cooling stage.
Microstructural Formulation
During the quenching many phases can be formed, such as perlite, bainite and martensite. The first two are diffusional controlled, they need a finite time for the transformation occurrence. The Johnson–Mehl–Avrami–Kolmogorov (JMAK) model can represent this transformation, using the follows equations:
𝐹𝑖(𝑇) = 1 − exp [−𝑎(𝑇𝑗). 𝑡(𝑇)𝑛(𝑇𝑗)] (14)
where 𝐹𝑖 is the volume phase fraction of the phase 𝑖 formed.
The diffusion coefficient, 𝑎, and the exponent, 𝑛, are properties
of the material and can be obtained by the steel TTT diagram, using 𝜏𝑠(𝑇) and 𝜏𝑒(𝑇), starting time transformation and
finishing time transformation for a given temperature T, as follows:
𝑛(𝑇) =ln (ln(𝐹𝑒)−ln (𝐹𝑠))
ln (𝜏𝑠(𝑇)−ln (𝜏𝑒(𝑇)) (15)
𝑎(𝑇) = − ln(𝐹𝑠) 𝜏𝑠(𝑇)−𝑛(𝑇) (16)
where 𝐹𝑒 and 𝐹𝑠 are the starting and finishing volume fraction
and assumed as 0.01 and 0.99 respectively.
The JMAK model is valid only to isothermal transformation, so to use this model we need to convert the cooling curve to a series of isothermal steps as shown in Figure 3.
93
Figure 3. Schematic TTT diagram show conversion of cooling curve to isothermal steps.
Now, which modified JMAK model, the time used in the equation is a sum of a fictitious time to the duration of time step. This fictitious time is based on the cumulative transformed volume fraction of the specific i microstructural phase.
𝑡𝑗 = ∆𝑡𝑗 + 𝑡𝑗,𝑓𝑖𝑐𝑡 (17)
𝑡𝑗 ,𝑓𝑖𝑐𝑡 (𝑇𝑗) = [− ln(1−𝐹𝑖,𝑗−1)
𝑎(𝑇𝑗)]
1
𝑛(𝑇𝑗) (18)
Perlite and bainite phase have a incubation time before they can nucleate. To consider this aspect, we use Scheil’s Addictive Rule, and the transformation can begin if the following equation was satisfied.
∑∆𝑡𝑗
𝜏𝑠(𝑇𝑗)≥ 1𝑛
𝑗=1 (19)
Austenite-martensite transformation isn’t diffusional and can be modeled by the Koistinen-Marburguer equation. It have to be taken into account that a percentage of austenite had been transformed in perlite or bainite, so the Koistinen-Marburguer equation need to be multiplied by the fraction of the remaining austenite, as follow:
𝐹𝑚,𝑗 = [1 − exp[−𝛼(𝑇𝑚𝑠 − 𝑇)]](1 − ∑ 𝐹𝑖𝑖 ) (20)
Where Tms is the martensite start temperature.
During the cooling, for each time step, the properties 𝜙 (mass
density, specific heat and thermal conductivity) are updated based on the fraction of each phase, thus the properties are a function of the phase and the temperature.
𝜙 = ∑ 𝐹𝑖𝜙𝑖𝑛𝑖=1 (21)
During the phase transformation there is a release of the latent heat and it is also included in each time step, considering the variation of the volume fraction and enthalpy change of each phase. While there is heat being released, and consequently there is a increase of the temperature, it may occur a tendency of decreasing the transformed martensite volume fraction. A way to avoid this effect is maintaining the max value of transformed volume fraction of two consecutive time steps.
𝑄𝑡 = ∑∆𝐻𝑖 ,𝑗∆𝐹𝑖,𝑗
∆𝑡𝑗 (22)
𝐹𝑚,𝑗 = max[𝐹𝑚,𝑗 , 𝐹𝑚,𝑗−1] (23)
In the cooling phse, the Fourier equation (13) is also used to calculate the equilibrium thermal field.
Results
Validation
To validate the microstructural model, it was simulated a water quenching of a 38mm diameter and 76mm long cylinder. The cylinder was made of SAE 1080 steel and it was considered 100% austenite as initial phase and initial temperature of 850°C. Only austenite-perlite and austenite-martensite were considered in validation order to compare with literature results. The Figure 4 shows the convective coefficient of heat transfer of a cylinder quenching by immersion in water at 20°C versus the temperature, and it was obtained from [4].
Since the cylinder has a length/diameter ratio equal 2, it was modeled as an infinitely long cylinder [5] with 95 elements along the radius. In [7] it was used a mesh with 10 elements along the length and 20 along the radius.
Figure 4. Variation of convective coefficient with temperature of a 38mm diameter cylinder water quenching, obtained from Woodard et al., 1999.
94
Figure 5. Phase distribution along the cylinder after 90s of quenching.
Despite the martensite transformation model being different from [7], it was used Koistinen-Marburguer in present validation case and Yu Parabolic model in [7], there is a very good agreement between them.
Figure 5. Boundary Condition and location of points A, B and C.
To compare the temperature history it was picked 3 points: A= 0.20mm, B=8.10mm and C=18.96mm from my model, and A=0mm, B=8.08mm and C=19.05 from [7]. Figure 5 shows these 3 points and the boundary conditions applied to this simulation.
Figure 6. Temperature history of water quenching 38-diameter cylinder of 1080 steel.
Present Case
In the present case it was simulated the induction electromagnetic heating and the water quenching by immersion of a SAE 1080 steel cylinder. The dimensions of the cylinder are radius of 19 mm and 76 mm of length, but, for symmetry reasons, it was modeled only 38 mm of length. The inductor was modeled as eight wire cooper around the cylinder and it was neglected the angle between two wires, being possible to make a 2D model of the problem.
The cylinder mesh used in the simulation has 55 elements along the length and 28 along the radius. The elements near to the surface of the cylinder is finer than the elements near to the cylinder axis.
A 1500 𝐴 current was input in the inductor wire with 1.2mm of
radius, that operate at 𝑓 = 1000 𝐻𝑧. For simplicity, it wasn’t
considered the Joule Heat loss on the inductor wire.
The cylinder was initially at 20ºC, being heated by 20 seconds. The figure 8 shows the temperature field, in kelvin, after the heating stage. Only a small part of the cylinder reach a temperature above the SAE 1080 steel austenitizing temperature (996 K), so only these elements are subjected to microstructural transformation during the cooling. As already said, all the elements with temperature below the steel austenitizing temperature are going to be considered as 100% perlite.
In the heating the electromagnetic properties (electric conductivity and B-H curve) of the steel and the thermal properties (mass density, specific heat and thermal conductivity) of the steel are temperature depended, as already explained. It’s important to highlight that above Curie Temperature, the steel magnetic permeability is not obtained from B-H curve, but instead it became a constant coefficient 𝜇 = 𝜇0, where 𝜇0 is the free space magnetic permeability.
95
Figure 7. Schematic figure of the study object. Steel cylinder (grey), air (blue) and inductor wires (orange)
During the heating the heat lost by convection was neglected, because it is small compared to radiation loss [8]. In addition, during the cooling, the heat lost was modeled only as convection, with the convection coefficient showed above.
In this model it was adopted a constant value of hardness for each phase. The perlite hardness used was 39 Rockwell C (RC), martensite and retained austenite hardness was 39 RC [4]. Bainite hardness used was 43 RC, according the hardness scale of TTT diagram and quenching time.
Figure 8. Temperature distribution, in kelvin, after the heating stage.
Figure 9. Microestrutural distribution after the quenching.
Figure 10. Hardness after the quenching.
The mesh and the results were obtained through the commercial software Ansys. Unfortunately, we couldn’t get al the properties from the 1080 steel, the B-H curve is from SAE1045 steel, the electrical properties and thermal properties on heating is from E235, and the thermal properties of each microstructural phase on cooling phase is from 1080 steel.
Conclusions
The objective of this work was to simulate the process of the quenching of a SAE 1080 steel cylinder, the induction heating and the cooling. Previous works were used to validate the microstructural model and a good accuracy was achieved. Therefore, it can be said that the simulations were successfully performed. For next works, it is planned to achieve the experimental validation of the whole process using the same boundary conditions described in this paper. Additionally, it is desired in the future to implement this simulation methodology on carbon steel automotive crankshafts with also a further experimental validation.
References
1. Drobenko B., Hachkevych O. and Kuornyts’kyi T., “A mathematical simulation of high temperature induction
96
heating of electroconductive solids,” International Journal of Heat and Mass Transfer. 50:616-624, 2007.
2. Cajner F., Smoljan B. and Landek D.,“Computer simulation of induction hardening,” Journal of Materials Processing Technology. 157–158:55–60, 2004.
3. Barka N., Bocher P., Brousseau J., Galopin M.et al., “Modeling and sensitivity study of the induction hardening process”, Advanced Materials Research. 15-17:525-530, 2007.
4. Woodard P., Chandrasekar S., Yang H., “Analysis of temperature and microstructure in the quenching of steel cylinders,” Metallurgical and Materials Transactions B. 30B:815-822, 1999.
Huiping L., Guoqun Z., Shanting N. and Chuanzhen H., “FEM simulation of quenching process and experimental verification of simulation results,” Materials Science and Engineering A. 452-453:705-714, 2007.
ANSYS, 11.0 Documentation, SAS IP, Canonsburg, 2007. Carlone P., Palazzo G. and Pasquino R., “Finite element
analysis of the steel quenching process: temperature field and solid-solid phase change,” Computers and Mathematics with Applications. 59:585-594¸ 2010.
Cho K., “Coupled electro-magneto-thermal model for induction heating process of a moving billet,” International Journal of Thermal Science. 60:195-204, 2012.
Doane D., “Application of hardenability concepts in heat treatment steel,” J. Heat Treatment.1:5-30, 1979.
Hömberg D., “A numerical simulation of the Jominy End-Quench test,” Acta Mater. 44:4375-4385, 1996.
Contact Information
Pedro Augusto Lanza de Paula. E-mail: [email protected] Faculdade de Engenharia Mecânica - UNICAMP Rua Mendeleyev, 200 13083-860 Campinas – SP, Brasil Renato Pavanello Email: [email protected] Faculdade de Engenharia Mecânica - UNICAMP Rua Mendeleyev, 200 13083-860 Campinas – SP, Brasil MSc. Wiliam Tean Su, Research and Development Engineer Email: [email protected] Avenida Alfried Krupp, 1050 13231-900 Campo Limpo Paulista –SP, Brasil MSc. Alex de Souza Rodrigues, Research and Development Head of Division Email: [email protected] Avenida Alfried Krupp, 1050 13231-900 Campo Limpo Paulista –SP, Brasil
97
ANEXO B – Artigo publicado no Congresso SAE 2015
2015-36-0391
Sensitivity analysis of multiparameter numerical model of quenching
process using electromagnetic induction heating.
Pedro de Paula(1)
, Renato Pavanello(1)
, Wiliam Su(2)
, Alex Rodrigues(2)
(1) Department of Computational Mechanics, State University of Campinas,
(2) Department of Research and Development Engineering, ThyssenKrupp Metalúrgica Campo Limpo Ltda.
Abstract
Induction hardening process is widely used to improve fatigue
strength of mechanical components that are subjected to cyclic
loads in service. The depth of the hardened layer is directed
linked with the fatigue and impact strength. So, to improve the
mechanical properties in order to preventing fatigue failure in
service, it is very important to understand the process and the
influence of its parameters.
In this paper, a sensitivity study of the influence of some process
parameters on the hardness profile of a crankshaft’s crankpin
after induction hardening using will be presented. The proposed
simulation method include two stages: heating and cooling. In
the first stage, the mechanical component, initially at ambient
temperature, is heated by electromagnetic induction to a
temperature above the steel austenitization. In the second one,
the component is cooled by liquid immersion. The resulting
microstructure is calculated using the Johnson–Mehl–Avrami–
Kolmogorov model and Sheil’s additive rule for diffusional
transformation, while austenite-martensite transformation is
calculated by Koistinen-Marburguer equation.
The simulation takes into account thermal properties as a
function of temperature and microstructure and electromagnetic
properties as a function of temperature and strength of the
electromagnetic field (magnetic permeability). The results of the
sensitivity study can help in the process setup, decreasing the
number of experimental tests.
Introduction
Induction hardening is an important industrial process used to
improve the mechanical properties of components subject to
cyclic loading, which induce the appearance of fatigue cracks
and may lead to failure of such component in service. The
induction hardening process consists of heating the material to a
temperature above the austenitization temperature and cools it
quickly in order to obtain a high hardness surface layer,
composed principally of martensite. Such a surface layer is
responsible for increasing resistance to fatigue, preventing the
crack initiation, and increasing wear resistance. It is also
necessary that the mechanical component, hold a tenaciously
core after the heat treatment process, in order to maintain good
impact resistance.
The first stage of the quenching, as said earlier, is to heat the
mechanical component to a temperature above the
austenitization temperature. An efficient method used in the
heating phase at the process is the electromagnetic induction.
This type of heating has some advantages such as: greater
process control, possibility to heat only one region and reduce
the environmental impact in addition to concentrating the heat
input on the surface, saving energy.
Electromagnetic induction heating consists of inducing electrical
currents in the mechanical component which wants to heat by
introducing it in a variable electromagnetic field. That currents
will inducted in the surface of the mechanical component and
dissipate heat by Joule effect, increasing the temperature. This
electromagnetic field is created by an inductor in which
circulates an electric current of high frequency.
The simulation method proposed in this paper is divided in two
stages: heating and cooling. In both the finite element method is
used. Due to the physical characteristics of the induction heating
problem, such as the small penetration depth of the
electromagnetic field, skin depth in good electrical conducting
materials, it is necessary the use of a finer mesh on the surface
of the mechanical component containing from 3 to 10 elements
in the region of heat generation [2].
During the electromagnetic heating, it was observed that for
temperatures higher than 300°C and 600°C, if the variation of
electrical conductivity and magnetic permeability, respectively,
are not considered, it results in considerable errors in the field of
temperatures [1]. Among the properties like electrical resistivity,
specific heat, thermal conductivity and magnetic permeability
through a sensitivity study, that the first two are the ones that
have more influence in the depth of the hardened layer [3].
It is important to note that, despite the sinusoidal behavior of the
current in the inductor, the electromagnetic field inside the
ferromagnetic mechanical component is not sinusoidal because
the magnetic permeability of ferromagnetic material is a
nonlinear function of magnetic field strength, so the simulation
of the transient phenomenon can give more accurate results [1].
As the harmonic fields solution is simpler and requires less
computational effort, some authors have proposed models in
which they obtain an equivalent permeability, which is a
function of temperature only, in order to obtain a similar
distribution of eddy currents, and consequently temperatures
98
curves, to the simulation which the magnetic permeability as a
function of temperature and magnetic field intensity [6]. The
latent heating is being taken into account in the most part of the
work and, despite history temperature is more intensely
influenced near to the cylinder symmetry axis, to disregard it
cause an important alteration in the superficial microstructural
resulting [4].
The latent heat is not considered in the heating stage, because
the heat contribution of the induction heating is much higher
than it. Another heat source during the quenching, is which
released by the deformations, however its represents only 3% of
the total heat released and can be neglected [5].
The principal contribution of this paper is the implementation,
testing and multiparameter sensitivity analysis of a multiphysics
procedure to perform an integrated simulation of the quenching
process of mechanical components considering electromagnetic,
thermal and microstructural analysis.
Modeling
Electromagnetic Formulation
The proposed integrated multiphysics simulation include two
principal stages: electromagnetic and thermic coupled and
thermic and microstructural coupled.
The electromagnetic problem in consideration is governed by
the Maxwell’s equations:
∇x𝐇 = 𝐉 +𝝏𝑫
𝝏𝒕
(1)
∇x𝐄 = − 𝝏𝑩
𝝏𝒕 (2)
∇.𝑩 = 0
(3)
∇.𝑫 = 𝜌𝑒 (4)
where 𝐇 is the magnetic field strength, 𝐄 is the electric field
intensity, 𝐉 is the conduction current density, 𝑫 is the electric
flux density, 𝜌𝑒 is the electric charge density. In good conductor
materials and considering harmonic magneto-quasi-static
problem, the displacement current is negligible compared with
conduction current, so we can write the Maxwell’s equations as
follows:
∇x𝐇 = 𝐉
(5)
∇x𝐄 = − 𝝏𝑩
𝝏𝒕 (6)
∇.𝑩 = 0
(7)
The constitutive equations described below supplement
equations (5)-(7) and describes the behavior of the
electromagnetics materials.
𝑱 = 𝜎𝑬
(8)
𝑩 = 𝜇𝑯
(9)
Here, 𝜎 is the electric conductivity and 𝜇 is magnetic
permeability and 𝑩 is the magnetic flux density. In the present
model, 𝜇 is function of the temperature and of the magnetic field
intensity, so it was used the following approximation to describe
the variation of the 𝑩-𝐇 curve with the temperature:
𝑩 = 𝜇(𝐻, 𝑇)𝑯 = 𝜇20(𝐻) [1 − (𝑇
𝑇𝑐𝑢𝑟𝑖𝑒)6]𝑯
(10)
Where 𝜇20 is obtained from the steel 𝑩-𝐇 curve at 20°C and
𝑇𝑐𝑢𝑟𝑖𝑒 is the Curie temperature of the material, in which the steel
lose its ferromagnetic properties, i.e., the steel magnetic
permeability above Curie point is considered equal to the free
space magnetic permeability.
As said, despite the electromagnetic field generated in the coil
be sinusoidal the non-linearity of the magnetic permeability of
the ferromagnetic material, made the field inside them non-
sinusoidal. However, for power loss time-average problem, we
can substitute the ferromagnetic material by a fictitious material,
based on energy equivalence. The material B-H curve is
replaced by an effective B-H curve, transforming the nonlinear
transient problem in a nonlinear time-harmonic. Considering
the problem as a time harmonic, it is possible to reduce the
computational cost and maintaining a good accuracy [6]. In this
case, the magnetic field strength is approximated by:
1
2∫ 𝐻𝑚𝑑𝐵𝑒𝑓𝑓𝐵𝑒𝑓𝑓0
= 4
𝑃∫ (∫ 𝐻𝑚 sin(𝜔𝑡) 𝑑𝐵
𝐵
0)
𝑃
40
𝑑𝑡
(11)
where 𝐵𝑒𝑓𝑓 is the effective magnetic flux density, 𝑃 the time
period, 𝐻𝑚 the peak value of magnetic field and 𝜔 the angular
velocity. The output off electromagnetic phase is:
𝑄 =1
2𝑅𝑒(𝑬𝑇𝑱∗)
(12)
where 𝑄 is the volumetric heat power generate based on the
Joule Effect. The thermal simulation of the heating phase is
based in the Fourier's equations as follows:
𝜌𝑐𝜕𝑇
𝜕𝑡= ∇ . (𝑘∇𝑇) + 𝑄 (13)
where 𝜌 is the mass density, 𝑘 is the thermal conductivity and 𝑐
is the specific heat.
The algorithm of the heating phase is shown in Figure 1. The
temperature field obtained from the Thermal Phase I, at heating
stage, is used as initial conditional on Thermal Phase II, at
cooling stage shown in Figure 2.
In the heating stage the purpose is to evaluate the volumetric
heat generated by the Joule Effect on the electromagnetic phase
and it is an input to thermal phase I. This volumetric heat remain
99
constant during a thermal time step. After the thermal phase I is
done, the temperature filed is updated, as well as the properties,
and another electromagnetic simulation is done, as shown in
Figure 1, obtaining a new volumetric heat input to the next time
step. This sequence runs until it got the total time heating
desired. After the heating, the microstructural phase of the
elements with temperature above the austenitization temperature
is set as 100% austenite. So there is another thermal step
representing the delay time before the component be submitted
to quenching, used to the heat generated in the surface of the
component diffuse to the core, and the austenitized region is
evaluated again.
In the cooling stage, Figure 2, the final temperature field
obtained from the heating stage and delay time after quenching
is used as initial condition. The temperature filed is evaluated
for each time step and the microstructural phase percentage, F, is
updated using the microstructural model. So, the material
properties and the heat released by the microstructural
transformation is updated and another thermal phase II is
performed. These sequence runs until the total time of cooling
phase is reached.
Figure 1. Simulation strategy of the heating stage.
Figure 2. Simulation strategy of the cooling stage.
Microstructural Formulation
During the quenching many phases can be formed, such as
perlite, bainite and martensite. The austenite-perlite and
austenite-bainite transformation are diffusional, they need a
finite time for the transformation occurrence. The Johnson–
Mehl–Avrami–Kolmogorov (JMAK) model can represent this
transformation, using the follows equations:
𝐹𝑖(𝑇) = 1 − exp [−𝑎(𝑇). 𝑡(𝑇)𝑛(𝑇)]
(14)
where 𝐹𝑖 is the volume phase fraction of the phase 𝑖 formed. The
diffusion coefficient, 𝑎, and the exponent, 𝑛, are properties of
the material and can be obtained by the steel TTT diagram,
using 𝜏𝑠(𝑇) and 𝜏𝑒(𝑇), starting time transformation and
finishing time transformation for a given temperature T, as
follows:
𝑛(𝑇) =ln (ln(𝐹𝑒)−ln (𝐹𝑠))
ln (𝜏𝑠(𝑇)−ln (𝜏𝑒(𝑇))
(15)
𝑎(𝑇) = − ln(𝐹𝑠) 𝜏𝑠(𝑇)−𝑛(𝑇)
(16)
where 𝐹𝑒 and 𝐹𝑠 are the starting and finishing volume fraction
and assumed as 0.01 and 0.99 respectively.
The JMAK model is valid only to isothermal transformation, so
to use this model we need to convert the cooling curve to a
series of isothermal steps as shown in Figure 3.
100
Figure 3. Schematic TTT diagram show conversion of cooling curve to isothermal steps.
Now, with modified JMAK model, the time used in the equation
(19) is a sum of a fictitious time to the duration of time step.
This fictitious time is based on the cumulative transformed
volume fraction of the specific i microstructural phase.
𝑡𝑗 = ∆𝑡𝑗 + 𝑡𝑗,𝑓𝑖𝑐𝑡
(17)
𝑡𝑗 ,𝑓𝑖𝑐𝑡 (𝑇𝑗) = [− ln(1−𝐹𝑖,𝑗−1)
𝑎(𝑇𝑗)]
1
𝑛(𝑇𝑗) (18)
𝐹𝑖,𝑗(𝑇𝑗) = 1 − exp [−𝑎(𝑇𝑗). 𝑡(𝑇𝑗)𝑛(𝑇𝑗)
]
(19)
Perlite and bainite phase have an incubation time before they
can nucleate. To consider this aspect, we use Scheil’s Addictive
Rule, and the transformation can begin if the following equation
was satisfied.
∑∆𝑡𝑗
𝜏𝑠(𝑇𝑗)≥ 1𝑛
𝑗=1
(20)
Austenite-martensite transformation isn’t diffusional and can be
modeled by the Koistinen-Marburguer equation. It have to be
taken into account that a percentage of austenite had been
transformed in perlite and bainite, so the Koistinen-Marburguer
equation need to be multiplied by the fraction of the remaining
austenite, as follow:
𝐹𝑚,𝑗 = [1 − exp[−𝛼(𝑇𝑚𝑠 − 𝑇)]](1 − ∑ 𝐹𝑖𝑖 )
(21)
Where 𝑇𝑚𝑠 is the martensite start temperature.
During the cooling, for each time step, the properties 𝜙 (mass
density, specific heat and thermal conductivity) are updated
based on the fraction of each phase, thus the properties are a
function of the phase and the temperature.
𝜙 = ∑ 𝐹𝑖𝜙𝑖𝑛𝑖=1
(22)
During the phase transformation there is a release of the latent
heat and it is included in each time step, considering the
variation of the volume fraction and enthalpy change of each
phase. While there is heat being released, and consequently
there is an increase of the temperature, it may occur a tendency
of decreasing the transformed martensite volume fraction. A
way to avoid this effect is maintaining the max value of
transformed volume fraction of two consecutive time steps.
𝑄𝑡 = ∑∆𝐻𝑖 ,𝑗∆𝐹𝑖,𝑗
∆𝑡𝑗
(23)
𝐹𝑚,𝑗 = max[𝐹𝑚,𝑗 , 𝐹𝑚,𝑗−1]
(24)
In the cooling phase, the Fourier equation (13) is also used to
calculate the equilibrium thermal field.
Results
In the present case it was made a 2D model simulation to study
the influence of the current frequency and current amplitude in
the hardness profile and in the total energy absorbed from the
inductor of a crankshaft crankpin quenching process heated by
electromagnetic induction. The inner diameter of crankpin is
69mm and it is made of SAE 1080 steel. The inductor was
modeled as two wires evolved by a flux concentrator, Fluxtrol
50®, with 4mm thickness. The flux concentrator works like a
magnifier of the electromagnetic field, directing the flux to the
desired localization for heating. Also, the angle between the two
wires was neglected, being possible to make a 2D model of the
problem. The figure 4 shows the upper half of the ¾
axisymmetric expansion of the model, the gray region symbolize
the crankpin, the orange region presents the inductor and the
magenta, the flux concentrator.
Figure 4. ¾ Axisymmetric expansion of the model.
101
The crankshaft was initially at 20ºC, being heated by 10
seconds. Only a small part of the crankpin reach a temperature
above the SAE 1080 steel austenitization temperature (996 K),
thus are subjected to microstructural transformation during the
cooling. All the elements which didn’t reach the steel
austenitization temperature were considered as 100% perlite.
In the heating, the electromagnetic properties (electric
conductivity and B-H curve) of the steel and the thermal
properties (mass density, specific heat and thermal conductivity)
of the steel are temperature depended, as already explained. It’s
important to highlight that above Curie Temperature, the steel
magnetic permeability is not obtained from B-H curve, but
instead it become a constant coefficient 𝜇 = 𝜇0, where 𝜇0 is the
free space magnetic permeability.
Figure 5. Variation of convective coefficient with temperature of a 38mm diameter cylinder water quenching, obtained from [4].
During the heating, the heat lost by convection was neglected,
because it is small compared to radiation loss [8]. In addition,
during the cooling, the heat lost was modeled only as
convection, with the convection coefficient showed above. The
ambient temperature used in the simulation was 20ºC.
The crankpin was cooled by immersion in 22.5ºC water, and the
convective coefficient of heat transfer used [4] is shown in the
figure 5.
In this model it was adopted a constant value of hardness for
each phase. The perlite hardness used was 39 Rockwell C (RC),
martensite and retained austenite hardness was 66 RC [4].
Bainite hardness used was 43 RC, according the hardness scale
of TTT diagram and quenching time.
Figure 6. Directions A and B which the hardened layer thickness was archived.
The figures (7 to 10) below shows the influence of the current
input (figure 10) and frequency input ( figure 9) on the hardened
layer thickness in the direction A and B, as shown in figure 6. It
is shown also, the power absorbed in each time step by the
crankpin according with current (figure 7) and frequency input
(figure 8). During the sensitivity study of current influence, it
was used 9000 Hz and 1s of delay time and in the sensitivity
study of frequency influence it was used 4200 A and 1s of delay
time.
Figure 7. Power absorbed by the crankpin in each time step for different inductor current frequencies input. Time step = 0.25s.
0
20
40
60
80
100
120
140
0 5 10 15 20 25 30 35 40 45
Po
wer
Ab
sorb
ed [
KW
]
Time step
3000 Hz
4000 Hz
5000 Hz
6000 Hz
7000 Hz
8000 Hz
9000 Hz
102
Figure 8. Power absorbed by the crankpin in each time step for each
inductor current. Time step = 0.25s.
Figure 9. Hardness thickness in the directions A and B after quenching for inductor current frequency variation.
Figure 10. Hardness thickness in the directions A and B after quenching
for inductor current variation.
It was studied the influence of the magnetic flux concentrator
too. The power absorbed by the crankpin and the final hardness
profile was observed in three situations: 1 – with magnetic flux
concentrator, 2 – no magnetic flux concentrator and 3 – with
partial flux concentrator. The configurations 1 and 3 are shown
the figures 11 and 12, respectively. The input data for these
cases was 9000hz of current frequency, 1s of delay time and
4000A of current amplitude.
In all cases, all the properties used was from SAE 1080 steel,
unless the B-H curve, which is from SAE 1045 steel.
Figure 11. Configuration 1 – with magnetic flux concentrator.
Figure 12. Configuration 3 – with partial magnetic flux concentrator.
The figure 13 shows the power absorbed in the three
configurations and the figures 14, 15 and 16 show the hardness
profile of each configuration, after the quenching.
0
50
100
150
200
0 5 10 15 20 25 30 35 40 45
Po
wer
Ab
sorb
ed [
KW
]
Time Step
3400 A3600 A3800 A4000 A4200 A4400 A4600 A
0
1
2
3
4
5
6
7
8
0 2000 4000 6000 8000 10000
Har
den
ed L
ayer
Th
ckn
ess
[mm
]
Frequency [Hz]
Direction A
Direction B
0
1
2
3
4
5
6
7
8
3000 3500 4000 4500 5000
Har
den
ed L
ayer
Th
ickn
ess
[mm
]
Current [A]
Direction A
Direction B
103
Figure 13. Power absorbed by the crankpin in each time step for three magnetic flux concentrator. Time step = 1s.
Figure 14. Hardness profile after quenching for configuration 1 – with magnetic flux concentrator.
Figure 15. Hardness profile after quenching for configuration 2 – with
no magnetic flux concentrator.
Figure 16. Hardness profile after quenching for configuration 1 – with partial magnetic flux concentrator
Conclusions
The objective of this work was to simulate the process of the
quenching of a SAE 1080 steel crankpin and make a sensitivity
analysis of frequency, current input and geometry of magnetic
flux concentrator in the hardness profile and power absorbed.
In the cases studied, it was observed that higher frequencies and
currents input produce, most of cases, bigger hardened layer
thickness, because of the higher power input produced by them.
It was observed too, that after 5000 Hz the hardened layer
thickness in the direction B decreased while the layer thickness
in direction A increased. It could have occurred because of
different power distribution caused by each frequency input.
Also, it was observed the improvement caused by the magnetic
flux concentrator, producing a higher power input and a more
homogeneous hardened layer in the cases studied.
For next works, it is planned to improve to physical model and
explore more parameters like distance of the inductor from the
piece and geometry of magnetic flux concentrator and inductor.
References
5. Drobenko, B., Hachkevych, O. and Kuornyts’kyi, T., “A
mathematical simulation of high temperature induction heating
of electroconductive solids,” International Journal of Heat and
Mass Transfer. 50:616-624, 2007.
6. Cajner, F., Smoljan, B. and Landek, D.,“Computer simulation of
induction hardening,” Journal of Materials Processing
Technology. 157–158:55–60, 2004.
7. Barka, N., Bocher, P., Brousseau, J., Galopin, M.et al.,
“Modeling and sensitivity study of the induction hardening
process”, Advanced Materials Research. 15-17:525-530, 2007.
0
20
40
60
80
100
120
0 5 10 15 20 25 30 35 40 45
Po
wer
Ab
sorb
ed [
KW
]
Time step
With Flux Concentrator - 1
No Flux Concentrator - 2
Partial Flux Concentrator - 3
104
8. Woodard, P., Chandrasekar S., Yang H., “Analysis of
temperature and microstructure in the quenching of steel
cylinders,” Metallurgical and Materials Transactions B.
30B:815-822, 1999.
Huiping, L., Guoqun, Z., Shanting, N. and Chuanzhen, H., “FEM
simulation of quenching process and experimental verification
of simulation results,” Materials Science and Engineering A.
452-453:705-714, 2007.
ANSYS, 11.0 Documentation, SAS IP, Canonsburg, 2007.
Carlone, P., Palazzo, G. and Pasquino, R., “Finite element analysis of
the steel quenching process: temperature field and solid-solid
phase change,” Computers and Mathematics with Applications.
59:585-594¸ 2010.
Cho, K., “Coupled electro-magneto-thermal model for induction
heating process of a moving billet,” International Journal of
Thermal Science. 60:195-204, 2012.
Doane, D., “Application of hardenability concepts in heat treatment
steel,” J. Heat Treatment.1:5-30, 1979.
Hömberg, D., “A numerical simulation of the Jominy End-Quench
test,” Acta Mater. 44:4375-4385, 1996.
Avrami, M., “Kinects of Phase Change II”, Journal of Chemical
Physics. 8: 212-224, 1940.
Li, F., Li, X., Zhu, T. and Rong, Y., “Numerical Simulation of the
Moving Induction Heating Process with Magnetic Flux
Concentrator,” Advances in Mechanical Engineering. 2013: 1-9,
2013.
De Paula, P., Pavanello, P., Su, W. and Rodrigues, A. “Multiphysics
simulation of quenching process of SAE 1080 steel cylinder,
coupling electromagnetic, thermal and microstructural analysis.”
SAE Technical Paper 2014-36-0425, 2014.
Contact Information
Pedro Augusto Lanza de Paula,
E-mail: [email protected]
Faculdade de Engenharia Mecânica - UNICAMP
Rua Mendeleyev, 200
13083-860 Campinas – SP, Brasil
Renato Pavanello,
Email: [email protected]
Faculdade de Engenharia Mecânica - UNICAMP
Rua Mendeleyev, 200
13083-860 Campinas – SP, Brasil
MSc. Wiliam Tean Su,
Research and Development Engineer
Email: [email protected]
Avenida Alfried Krupp, 1050
13231-900 Campo Limpo Paulista –SP, Brasil
MSc. Alex de Souza Rodrigues.
Research and Development Head of Division
Email: [email protected]
Avenida Alfried Krupp, 1050
13231-900 Campo Limpo Paulista –SP, Brasil
![Tensões residuais térmicas obtidas da têmpera e martêmpera ... · Difusividade térmica [m²/s] K ... aquecimentos e com um aquecimento final até a temperatura de 1025 °C. Na](https://img.document.onl/doc/110x75/5c4a15d793f3c3143644a649/tensoes-residuais-termicas-obtidas-da-tempera-e-martempera-difusividade.jpg)
