tana_lfb_me_guara

Prévia do material em texto

Modelagem Matemática do Aquecimento Indutivo de Barras de Aço Utilizadas na 
Produção de Hastes de Bombeio 
Guaratinguetá - SP 
2016 
LUIZ FERNANDO BISSOLI TANA 
 
 
 
 
 
 
Luiz Fernando Bissoli Tana 
 
 
 
 
 
 
Modelagem Matemática do Aquecimento Indutivo de Barras de Aço Utilizadas na 
Produção de Hastes de Bombeio 
 
 
 
 
 
 
 
 
 
Dissertação apresentada à Faculdade de 
Engenharia do Campus de 
Guaratinguetá, Universidade Estadual 
Paulista, para a obtenção do título de 
Mestre em Engenharia Mecânica na área 
de Transmissão e Conversão de Energia 
 
 
Orientador: Prof. Dr. Maurício Araujo 
Zanardi 
 
 
 
 
 
 
 
Guaratinguetá 
2016 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
T161m 
Tana, Luiz Fernando Bissoli 
 Modelagem matemática do aquecimento indutivo de barras de aço / 
Luiz Fernando Bissoli Tana – Guaratinguetá, 2016. 
 94 f : il. 
 Bibliografia: f. 91-94 
 
 Dissertação (Mestrado) – Universidade Estadual Paulista, Faculdade de 
Engenharia de Guaratinguetá, 2016. 
 Orientador: Prof. Dr. Maurício Araújo Zanardi 
 
 
 
1. Aquecimento. 2. Barras de aço. 3. Modelos matemáticos. I. Título 
 
CDU 620.91(043) 
 
 
 
 
 
 
 
 
 
 
AGRADECIMENTOS 
 
À Tenaris Confab Hastes de Bombeio que permitiu minha dedicação ao tema e forneceu 
as informações para o desenvolvimento da tese, ao Professor Maurício que sempre se 
mostrou disposto a ajudar e a direcionar o andamento do trabalho, e especialmente à 
minha esposa que é minha fiel companheira que sempre me incentivou a buscar este 
honroso mérito. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 “A pessoa que diz que não deve ser feito não deve 
interromper a pessoa que está fazendo.” 
Provérbio Chinês 
 
 
 
TANA, L. F. B. Modelagem Matemática do Aquecimento indutivo de Barras de 
Aço utilizadas na produção de hastes de bombeio. 2016. 97 f. Trabalho de Pós 
Graduação (Pós Graduação em Engenharia Mecânica) – Faculdade de Engenharia do 
Campus de Guaratinguetá, Universidade Estadual Paulista, Guaratinguetá, 2016. 
 
RESUMO 
 
Neste trabalho, aborda-se a modelagem por elementos finitos de um aquecimento 
indutivo de barras de aço utilizadas no forjamento de hastes de bombeio, tendo como 
principal objetivo avaliar a variação de temperatura radial da barra de aço aquecida pelo 
indutor e o impacto das propriedades envolvidas no aquecimento. Uma apresentação do 
processo de produção de uma haste de bombeio é efetuada, descrevendo as etapas de 
produção, sendo o aquecimento indutivo a etapa em estudo. A técnica de aquecimento 
por indução trata-se de um processo com alto grau de repetição e alto volume de 
produtividade, utilizada para aquecimento de peças dos mais diversos formatos, sendo 
amplamente empregada na indústria. A base para a modelagem matemática é feita por 
meio das quatro equações de Maxwell relacionadas às propriedades magnéticas dos 
materiais. O método utilizado para os elementos finitos é de Resíduos Ponderados 
(Galerkin) que possui grande aproximação dos resultados. É feito um detalhamento da 
modelagem matemática do método aplicado no programa utilizado. 
 
PALAVRAS-CHAVE: Aquecimento Indutivo. Elementos Finitos. Método de 
Galerkin. 
 
 
 
 
TANA, L. F. B. Matemathical Modeling of Induction Heating of steel bars use in 
the production of sucker rods. 2016. 97 f. Post Graduate Work (Post Graduate 
Mechanical Engineering) – Faculdade de Engenharia do Campus de Guaratinguetá, 
Universidade Estadual Paulista, Guaratinguetá, 2016. 
 
ABSTRACT 
 
In this work, it is made an approach of finite elements modeling of an induction heating 
of steel bars used to forging sucker rods, seeking as the main goal to evaluate the radial 
thermal variation inside the bar heated by the inductor, since the impact of the variation 
of properties related to the heating process. A brief presentation of the production 
process is presented, describing all the productions steps, with the induction heating 
being the focus of the study. The technique of induction heating is a procedure with 
high repeatability and high volume productivity used for heating pieces of various kinds 
of formats, being widely employed in industry. The basis for the mathematical modeling 
is made using the four Maxwell equations related to the magnetic properties of 
materials. The method used to solve the equation was the finite element residues 
(Galerkin) which has great approximation of the results. It was made a detailed 
mathematical modeling of the method applied in the software used. 
 
KEYWORDS: Induction Heating, Finite Elements. Galerkin Method. 
 
 
 
 
LISTA DE FIGURAS 
 
Figura 1 – Engrenagem aquecida por indução eletromagnética .............................................. 16 
Figura 2 – Solda de tubo por aquecimento indutivo ................................................................ 17 
Figura 3 – Indutor Aberto utilizado em brasagem ................................................................... 18 
Figura 4 – Profundidade de aquecimento em geometrias complexas .................................... 19 
Figura 5 – Fluxograma de fabricação de uma Haste de Bombeio ........................................... 21 
Figura 6 – Variação das Pressão de forjamento com a temperatura ........................................ 22 
Figura 7 – Haste Forjada ......................................................................................................... 22 
Figura 8 – Tratamento térmico pós forjamento para uma haste de bombeio .......................... 23 
Figura 9 – Variação de dureza do aço com o tempo de imersão ............................................. 23 
Figura 10 – Haste com a extremidade usinada ........................................................................ 24 
Figura 11 – Condução de calor por uma parede de espessura ∆x e área. ................................ 29 
Figura 12 – Fluxograma de uma análise por elementos finitos ............................................... 33 
Figura 13 – Tipos de Malhas de elementos finitos .................................................................. 34 
Figura 14 – Exemplo de malha para um elemento .................................................................. 34 
Figura 15 – Representação de um elemento triangular de 1° grau. ......................................... 35 
Figura 16 – Representação genérica de um elemento Finito ................................................... 38 
Figura 17 – Representação de uma aleta e suas condições de contorno. ................................. 39 
Figura 18 – Elemento com simetria axial ................................................................................ 53 
Figura 19 – Representação Gráfica dos métodos de integração numérica. ............................. 59 
Figura 20 – Esquema de Funcionamento do Aquecimento Indutivo ...................................... 70 
Figura 21 – Fluxograma de funcionamento do Programa ....................................................... 71 
Figura 22 – Detalhe da malha do indutor. ............................................................................... 72 
Figura 23 – Aquecimento barra com malha de 335 nós .......................................................... 78 
Figura 24 – Aquecimento barra com malha de 1111 nós ........................................................ 79 
Figura 25 – Aquecimento barra com malha de 3989 nós ........................................................ 79 
Figura 26 – Aquecimento barra com malha de 9881 nós. ....................................................... 79 
Figura 27 – Aquecimento barra com malha de 2907 nós, corrente de 420 A .........................81 
Figura 28 – Aquecimento barra com malha de 2907 nós, corrente de 430 A ......................... 82 
Figura 29 – Aquecimento barra com malha de 2907 nós, corrente de 410 A ......................... 83 
Figura 30 – Aquecimento barra com malha de 2907 nós, com aquecimento rotativo ............ 84 
Figura 31 – Aquecimento barra com malha de 2907 nós, com aquecimento rotativo, 
indutor helicoidal ..................................................................................................................... 85 
 
 
 
Figura 32 – Aquecimento barra com malha de 2907 nós, 5 espiras com aquecimento 
rotativo, indutor helicoidal ...................................................................................................... 85 
Figura 33 – Aquecimento barra com malha de 2907 nós, 3 espiras com aquecimento 
rotativo, indutor helicoidal ...................................................................................................... 87 
Figura 34 – Aquecimento barra com malha de 2907 nós, 6 espiras com aquecimento 
rotativo, f = 7500 Hz ............................................................................................................... 88 
Figura 35 – Aquecimento barra com malha de 2907 nós, 6 espiras com aquecimento 
rotativo, f = 6500 Hz. .............................................................................................................. 89 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
LISTA DE TABELAS 
 
Tabela 1 – Propriedade dos aços analisados a temperatura ambiente ..................................... 68 
Tabela 2 – Propriedade do ar e do cobre. ................................................................................ 70 
Tabela 3 – Resumo das análises efetuadas. ............................................................................. 90 
 
 
 
 
 
 
LISTA DE GRÁFICOS 
 
Gráfico 1 – Comparação entre a Temperatura da Aleta pelo Método de Ritz e equação 
hiperbólica ............................................................................................................................... 41 
Gráfico 2 – Comparação entre a Temperatura da Aleta pelo Método Variacional e 
equação hiperbólica ................................................................................................................ 44 
Gráfico 3 – Comparação entre a Temperatura da Aleta pelo equação hiperbólica e os 
métodos dos resíduos ............................................................................................................... 47 
Gráfico 4 – Variação do Calor específico do aço 4142 ........................................................... 68 
Gráfico 5 – Variação da Condutividade Térmica do aço 4142 ............................................... 69 
Gráfico 6 – Variação da Resistividade do aço 4142 ................................................................ 69 
Gráfico 7 – Variação da função de correlação de permeabilidade magnética do aço 4140 .... 70 
Gráfico 8 – Variação do campo Elétrico para um período de 0,02 s ....................................... 74 
Gráfico 9 – Variação do campo elétrico para um período de 2x10-4 s. .................................. 75 
Gráfico 10 – Impacto da variação do período de integração no campo elétrico ..................... 75 
Gráfico 11 – Impacto da variação do período de integração na potência elétrica gerada ....... 76 
Gráfico 12 – Impacto da Temperatura na variação no campo elétrico gerado. ....................... 77 
Gráfico 13 – Impacto da Temperatura, acima de curie, no campo elétrico gerado. ................ 77 
Gráfico 14 – Penetração do Campo Elétrico em faixas de temperatura. ................................. 80 
Gráfico 15 – Variação da Temperatura longitudinal na barra ................................................. 81 
Gráfico 16 – Variação da Temperatura radial na barra. .......................................................... 82 
Gráfico 17 – Variação da Temperatura longitudinal na barra, I = 430 A ............................... 83 
Gráfico 18 – Variação da Temperatura longitudinal na barra, I = 410 A ............................... 84 
Gráfico 19 – Variação da Temperatura longitudinal na barra, 6 espiras, indutor 
Helicoidal ................................................................................................................................ 85 
Gráfico 20 – Variação da Temperatura longitudinal na barra, 5 espiras, indutor 
Helicoidal ................................................................................................................................ 86 
Gráfico 21 – Variação da Temperatura longitudinal na barra, 3 espiras, indutor 
Helicoidal ................................................................................................................................ 87 
Gráfico 22 – Variação da Temperatura longitudinal na barra, 6 espiras, indutor 
Helicoidal, f = 7500 Hz ........................................................................................................... 88 
Gráfico 23 – Variação da Temperatura longitudinal na barra, 6 espiras, indutor 
Helicoidal, f = 6500 Hz ........................................................................................................... 89 
 
 
 
 
LISTA DE ABREVIATURAS E SIGLAS 
 
API American Petroleum Institute 
EF Elementos Finitos 
UHS Ultra High Strenght 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
LISTA DE SÍMBOLOS 
 
 
A 
 
Área 
 
m² 
�⃗� Indução Magnética / Campo Magnético T 
[B] Matriz gradiente dos fatores de forma m 
[B]T Transposta da matriz gradiente dos fatores de forma m 
[Ce] Matriz de Capacitância Elétrica m²/Ω 
[Ct] Matriz de Capacitância Térmica J/K 
cp
 Calor Específico J/kgK 
D⃗⃗ Densidade de Fluxo Elétrico C/m² 
�⃗� Campo Elétrico V/m 
E Intensidade de Campo Elétrico V/m 
f Frequência Hz 
{Fe} Vetor Força Elétrica Am/s 
{Fe*} Vetor Força Elétrica Am 
{Ft} Vetor Força Térmica W 
{Ft*} Vetor Força Térmica J 
h Coeficiente de Convecção W/m²K 
H⃗⃗ Intensidade de Campo Magnético A/m 
𝐽 Densidade de Corrente Elétrica A/m² 
J0 Densidade de Corrente Inicial A/m² 
k Condutividade Térmica W/mK 
lij Comprimento da aresta IJ m 
[Ke] Matriz de Elétrica m²/ohm 
[Kt] Matriz de Térmica W/K 
N Fator de forma - 
[N] Matriz de fator de forma - 
[N]T Matriz de fator de forma transposta - 
P Perímetro m 
Q Fluxo Térmico W 
ri
 Coordenada Radial do ponto i m 
T 
Tc 
Temperatura do Corpo 
Temperatura de Curie 
°C 
°C 
t Tempo s 
w Pulsação Hz 
z Cota relativa entre um campo magnético e a peça m 
zi Coordenada longitudinal do ponto i m 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 LISTA DE SÍMBOLOS GREGOS 
 
 

emi 
Constante Dielétrica 
Emissividade Térmica 
F/m 
- 
 Permeabilidade Magnética do Vácuo H/m 
r Permeabilidade Magnética Relativa H/m 
 Permeabilidade Magnética H/m 
 Resistividade Elétrica (σm)-1 
ste Constante de Stefan Bolztmann W/(m²K
4) 


Densidade 
Profundidade de Penetração 
kg/m³ 
- 
 Volume do elemento Finito 
 Fronteira do elemento Finito 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
SUMÁRIO 
 
1 INTRODUÇÃO ................................................................................................... 15 
1.1 REVISÃO BIBLIOGRÁFICA .............................................................................. 15 
1.2 OBJETIVO ............................................................................................................ 20 
1.3 JUSTIFICATIVA .................................................................................................. 20 
1.4 PROCESSO DE PRODUÇÃO DE UMA HASTE DE BOMBEIO ..................... 20 
2 FUNDAMENTAÇÃO TEÓRICA ...................................................................... 252.1 FORMULAÇÃO ELETROMAGNÉTICA ........................................................... 25 
2.2 MATERIAIS MAGNÉTICOS .............................................................................. 27 
2.3 FOMRULAÇÃO TÉRMICA ................................................................................ 28 
2.3.1 Condução .............................................................................................................. 28 
2.3.2 Convecção ............................................................................................................. 29 
2.3.3 Radiação ............................................................................................................... 31 
3 METODOLOGIA ................................................................................................ 32 
3.1 ELEMENTOS FINITOS ....................................................................................... 32 
3.2 MÉTODO DE GALERKIN .................................................................................. 47 
3.3 MODELAGEM MAGNÉTICA ............................................................................ 60 
3.4 MODELAGEM TÉRMICA .................................................................................. 64 
4 MATERIAIS E MÉTODOS ............................................................................... 68 
5 RESULTADOS E DISCUSSÃO ........................................................................ 74 
6 CONCLUSÕES .................................................................................................... 91 
7 REFERÊNCIAS BIBLIOGRÁFICAS .............................................................. 92 
15 
 
 
1 INTRODUÇÃO 
 
1.1 REVISÃO BIBLIOGRÁFICA 
 
O processo de forjamento constitui uma das primeiras formas de conformação mecânica 
utilizada pela humanidade, sendo que, dessa maneira, foram executadas as primeiras 
ferramentas e armas utilizadas. Para se forjar um metal é necessário que ele atinja 
temperaturas elevadas que o permitam ser conformado plasticamente sem a necessidade de 
que equipamentos de grande porte sejam aplicados. Para isso aquece-se o material a ser 
tratado a elevadas temperaturas, e dentre os principais métodos aplicados para o aquecimento, 
destaca-se o aquecimento por indução. 
De acordo com Sadeguipour, Dopkin e Li (1996), entre as vantagens do aquecimento 
indutivo, pode-se incluir o controle preciso da profundidade das superfícies aquecidas, sendo 
idealmente aplicado em processos automatizados. Segundo Shokouhmand e Ghaffari (2012), 
o aquecimento indutivo é um processo de aquecimento de materiais condutivos, sem contato, 
com alto grau de precisão, amplamente utilizado na indústria, dentre eles os segmentos 
automotivo e aeroespacial. O processo consiste na passagem de uma corrente elétrica 
alternada em um indutor de cobre, que gera um campo magnético, adjacente ao indutor que 
penetra na peça, realizando o aquecimento. Para maximizar a potência gerada por esse 
equipamento, reduzindo o consumo de energia, alguns materiais são aplicados, dentre eles, 
chapas de ferro silício, que são adicionadas na região externa do indutor de forma a 
concentrar mais fortemente o campo gerado na peça aquecida. 
Como principais aplicações do aquecimento indutivo, pode-se citar: 
 Fusão de Metais 
 Forjamento / Conformação a Quente 
 Extrusão 
 Tratamento térmico 
 Soldagem 
 Cura de Tintas e Adesivos 
 Revestimento de Ligas metálicas 
Dentre as diversas aplicações do aquecimento indutivo, destaca-se o tratamento térmico. 
Segundo Simsir e Gur (2001), o tratamento térmico é uma etapa muito importante na cadeia 
de produção, pois as propriedades mecânicas do produto a ser vendido são definidas nesse 
 16 
 
1 – Disponível em: http://www.efd-induction.com/en/Applications.aspx acesso em fev 2016 
processo. Diversos estudos trabalham, atualmente, para solução de problemas de modelagem 
de revenimento, tempera e normalização realizados por indução. 
Os estudos de Ferguson, Li e Freborg (2005) mostram que a distorção em aços que 
passam por diversos tratamentos térmicos, como carbonetação, tempera e revenimento, 
(devido as distorções internas do material) pode ser analisada por softwares de simulação, 
sendo que o principal fator de distorção é proveniente da mudança de fase. Zabbet e 
Azghandi (2012) apresentou um programa para cálculo de processo de revenimento de uma 
peça de aço, trabalhando também para prever a dureza no ponto tratado de uma barra de 
direção. Zhou et al. (2013) propôs um modelo que controla o fluxo de calor extraído de uma 
peça resfriada com água e aquecida por um indutor. Segundo Han et al. (2013), a seleção 
cuidadosa de parâmetros para o processo de tratamento térmico tem um impacto significativo 
na estrutura metalográfica e nas propriedades mecânicas do material. A Figura 1 representa a 
aplicação de um aquecimento em uma engrenagem. 
 
Figura 1: Engrenagem aquecida por indução eletromagnética. 
 
Fonte: Pagina EFD Induction¹ 
 
Outra aplicação importante do aquecimento indutivo é na soldagem de materiais, na 
qual se destaca a soldagem de tubos. Nesse processo, uma variável de extrema importância é a 
frequência do campo magnético. O estudo de Kim (2009), analisa o impacto da variação da 
frequência no processo de solda, bem como o impacto do comprimento aquecido na zona de 
fusão do material. 
http://www.efd-induction.com/en/Applications.aspx
17 
 
2 - Disponível em: http://www.sunnysteel.com/images/erw-pipe.jpg acesso em fev. de 2016 
Contudo, antes de se iniciar qualquer processo, faz-se necessário especificar o tipo de 
indutor a ser utilizado no aquecimento do material. O equipamento pode possuir as mais 
variadas formas, sendo definido de acordo com a geometria da obra a ser efetuada, como 
observa-se nas Figuras 1, 2 e 3. O comprimento do indutor também é crítico. Segundo 
Tavakoli, Karbaschi e Samavat (2011), um indutor muito longo pode gerar um aquecimento 
excessivo da extremidade, levando a distribuição de temperatura não uniforme na longitude 
da peça. A frequência também é um parâmetro na seleção do indutor. 
De acordo com Bodart, Boureau e Touzani (2001), dentre os parâmetros de definição de 
um equipamento de aquecimento indutivo, pode-se trabalhar com o formato do indutor, do 
forno e do cobre, frequência da corrente, tempo de aquecimento e potência de aquecimento, 
como benefício a ser obtido, um melhor controle da qualidade do aquecimento. Os estudos de 
Jang e Chiu (2007), mostram que a variação de frequência, no aquecimento de tubos, leva a 
uma diferença de temperatura entre o diâmetro externo e interno do material. Verificou-se, 
por meio da comparação entre simulações, que para três situações de folga entre o indutor e a 
peça, os melhores rendimentos são obtidos para as menores distâncias entre os sistemas. 
Contudo limitações de matéria prima, como empeno e deflexões na ponta, impedem que 
distâncias muito pequenas sejam aplicadas. 
 
Figura 2: Solda de tubo por aquecimento indutivo. 
 
Fonte: Pagina Sunnysteel² 
 
 
http://www.sunnysteel.com/images/erw-pipe.jpg
18 
 
3 – Disponível em: http://www.efd-induction.com/en/Applications.aspx Acesso em fev. 2016 
Figura 3: Indutor Aberto utilizado em brasagem. 
Fonte: Pagina EFD Induction³ 
 
Alguns estudos buscam focar na eficácia do processo, como o proposto por Huang e 
Huang (2010), nos quais os efeitos da uniformidade do aquecimento da superfície de 
materiais metálicos são estudados. Por meio da aplicação do método de Taguchi, estudou-se 
as melhores combinações de indutores para obter a melhor temperatura de processo com a 
melhor estabilidade possível. O Método Taguchi é uma ferramenta estatística para a 
concepção de um experimento, utilizando tabelas de matrizes ortogonais, que dependem do 
número de fatores e níveisque são utilizados para estudar o espaço de parâmetros. Inúmeras 
variáveis de entrada podem ser avaliadas para determinar a sua contribuição para a variável 
resposta (saída) utilizando um menor número de ensaios experimentais (BALONI; PATAK; 
CHANNIWALA, 2015). Trata-se de um processo robusto que substitui os experimentos 
fatoriais completos com apenas uma matriz ortogonal simples (ALBETRAN; DONG; LOW, 
2015). 
Os estudos de Park (2013), mostram que a eficiência de um sistema também passa pela 
aplicação de sistemas de isolamento no equipamento, para redução das perdas por convecção 
e radiação, sendo que esses valores podem chegar próximos a 9% da energia acumulada pelo 
material. Otimização de parâmetros de processo, estratégia adequada de aquecimento e 
19 
4 – Disponível em: http://www.efd-induction.com/en/Applications.aspx acesso em fev 2016 
sistema de isolamento adequados, são as maneiras mais efetivas de se aumentar a eficiência 
energética de aquecedores indutivos. 
O trabalho de Naar e Bay (2013), foca na obtenção de uma região de trabalho mais 
controlada para a zona termicamente afetada de aços no processo de tratamento térmico. 
Contudo, o aquecimento indutivo se mostra um grande vilão no consumo de energia. Segundo 
Goodwin et al. (2013), o custo anual de consumo de energia de um forno convencional pode, 
seguramente, passar de 1 milhão de dólares, garantindo assim a busca por estudos para reduzir 
o consumo de energia. 
Nos estudos realizados, Goodwin et al. (2013) inseriu-se um controle de feedback no 
modelo para controle e estabilidade no processo, isso se mostra importante, devido a criação 
de um modelo de controle de aquecimento do material para evitar descartes e retrabalhos. Isso 
ocorre uma vez que a leitura por pirometria ótica falha em momentos em que a peça passa por 
um reaquecimento ou solta as carepas, medindo-se assim um valor de temperatura menor do 
que o real, informando ao sistema para continuar aquecendo a peça. 
Dentre as novas tecnologias aplicadas, tem-se o aquecimento por multifrequência, que 
é utilizado para peças com geometria complexa, sendo um exemplo de fácil compreensão 
indicado na Figura 4, na qual pode-se observar que com a aplicação dom método, a 
profundidade do tratamento é mais uniforme, garantido assim uma peça com melhor 
qualidade. Os estudos de Homberg, Petzold e Rocca (2015) comparam a aplicação da 
multifrequência com situações de alta e média frequência na austenitização de engrenagens. 
Conclui-se que para médias frequências apenas a base do dente é aquecida, enquanto que para 
altas frequências apenas a ponta do dente é aquecida, obtendo-se resultados adequados 
somente com a nova técnica. 
 
Figura 4: Profundidade de aquecimento em geometrias complexas. 
 
Fonte: Pagina EFD Induction4 
http://www.efd-induction.com/en/Applications.aspx
20 
 
 
1.2 OBJETIVO 
O foco deste estudo é a modelagem do aquecimento indutivo utilizado no forjamento de 
barras de aço, para utilização em poços de petróleo, e a metodologia a ser aplicada se baseia 
nos Métodos de Elementos Finitos, técnica que será explicada na fundamentação teórica. 
Como objetivos específicos serão analisados os impactos das seguintes variáveis 
envolvidas na modelagem do aquecimento: 
 Impacto do tempo de integração 
 Impacto do tamanho da malha. 
 Impacto da variação da corrente 
 Impacto da Concentração do Campo Elétrico 
 Impacto da Variação da Frequência 
 
1.3 JUSTIFICATIVA 
 
A aplicação de Métodos de Elementos Finitos tem se mostrando uma ferramenta de apoio 
sólida para mapeamento de processos produtivos, e sua modelagem dedicada a situações 
específicas tem resultados mais satisfatórios do que aquelas realizadas por softwares 
genéricos. 
 
 
1.4 PROCESSO DE PRODUÇÃO DE UMA HASTE DE BOMBEIO 
 
A haste de bombeio é um produto utilizado no bombeamento de poços de petróleo 
onshore. É o primeiro elemento de transmissão de movimento desde o motor na superfície até 
a bomba no poço, transmitindo a energia necessária para o acionamento da bomba. O 
bombeamento mecânico alternativo é o sistema de produção de petróleo cru mais utilizado no 
mundo, é um processo quase contínuo de sucção e transferência de petróleo para a superfície. 
A produção de uma haste segue rígidas normas de qualidade que são regulamentadas 
pela API, que concede a cada empresa produtora do material uma licença de uso da marca. 
Pode se observar na Figura 5 que o processo de produção é composto de diversas etapas, 
sendo o aquecimento da barra e o forjamento as duas primeiras. 
O aquecimento indutivo para forjamento de hastes de bombeio é executado por um 
forno de indução com uma potência máxima de 200 kW, onde a temperatura alvo de 
21 
 
 
aquecimento pode variar de 1243°C a 1260°C, estando diretamente atrelado ao tipo de 
material a ser forjado e ao diâmetro. 
O forjamento das extremidades das barras é a principal etapa na fabricação do produto, 
e para garantir que a conformação ocorrerá da maneira mais eficaz, é necessário que o 
produto seja aquecido da maneira correta e eficiente. Na Figura 6, pode-se observar que 
quanto maior a temperatura de um aço menor é a pressão de forjamento necessária para 
conformar o material, porém, temperaturas muito altas fazem com que o material se deforme 
demais, gerando dobras e possíveis descartes de processo. 
 
Figura 5: Fluxograma de fabricação de uma Haste de Bombeio 
 
Fonte: Produção do próprio autor. 
 
 
22 
 
 
Figura 6: Variação das Pressão de forjamento com a temperatura. 
 
Fonte : ASM Handbook 14, 1991, p. 473 
O produto obtido através do forjamento é mostrado na Figura 7, na qual pode-se 
observar a complexidade da geometria obtida: 
 
Figura 7: Haste Forjada. 
 
Fonte : Elaborada pelo Autor 
 
O processo seguinte ao forjamento é o tratamento térmico, sendo a Figura 8 uma 
representação da curva de aquecimento que a barra passa no tratamento térmico. Na primeira 
etapa a barra forjada passa por um alívio de tensões, na qual o aço é normalizado a 825 °C, 
obtendo-se uma microestrutura homogeneizada, corrigindo os danos anteriores dos processos 
de lingotamento, laminação e forjamento. Posteriormente, a barra passa pelo revenimento, que 
é um aquecimento subcrítico, reduzindo-se as tensões internas e garantindo a tenacidade e 
ductilidade do aço para o trabalho. A Figura 9 demonstra o impacto do tempo de imersão nas 
23 
 
 
propriedades mecânicas de um aço com 0,8% de C. Observa-se que quanto maior a 
temperatura, menor será a dureza obtida e consequentemente menor tensão de escoamento. A 
eficiência do tratamento térmico é comprovada através de ensaios de tração realizados em 
todas as corridas de aço revenidas pelo equipamento. 
 
Figura 8: Tratamento térmico pós forjamento para uma haste de bombeio. 
 
Fonte: Produção do próprio autor. 
Figura 9: Variação de dureza do aço com o tempo de imersão. 
 
Fonte: ASM Handbook 4, 1991, p. 296 
24 
 
5 – Disponível em http://www.tenaris.com/en/Products/SuckerRods/BeamPumping/APIRods.aspx acesso em jul. 
de 2016 
 
Após o tratamento térmico, a barra apresenta uma crosta de óxido sobre a superfície e 
para retirá-la, o aço passa por um processo de jateamento com granalhas de aço esféricas, no 
qual o material apresenta aspecto adequado, permitindo a inspeção dos extremos forjados para 
identificação de possíveis dobras de material e outras falhas de forjamento. Contudo, o 
principal benefício do jateamento é a introdução de tensões residuais no produto, fazendo com 
que sua vida útil em operação seja prolongada. 
Outra etapa muito importante na fabricação de uma haste de bombeio é o processo de 
usinagem nos extremos forjados. Neste processo o produto toma sua forma final, como mostra 
a Figura 10, pois na extremidade da barra é usinado um perfil que possui uma geometria que 
alivia as tensões do extremo da barra para evitar falhas na áreausinada, e é laminada uma 
rosca que possui maior resistência ao desgaste e suporta maiores níveis de tensão. 
Por fim, o produto é embalado e imerso em uma tinta de base asfáltica de forma a evitar 
possíveis pontos de corrosão entre o armazenamento e o uso final. 
 
Figura 10: Haste com a extremidade usinada. 
 
Fonte: Pagina Tenaris na internet5 
 
 
 
 
 
 
http://www.tenaris.com/en/Products/SuckerRods/BeamPumping/APIRods.aspx
25 
 
 
2 FUNDAMENTAÇÃO TEÓRICA 
 
2.1 FORMULAÇÃO ELETROMAGNÉTICA 
 
Tensão alternada aplicada a uma bobina de indução irá resultar em uma corrente 
alternada, produzindo assim um campo magnético com a mesma frequência da corrente em 
sua vizinhança. O campo magnético criado gera correntes parasitas na peça dentro do indutor, 
causando, dessa forma, o aquecimento do elemento pelo efeito Joule. O campo 
eletromagnético é estudado pelas equações de Maxwell. 
O conjunto de equações de Maxwell, em sua forma completa, é dado pelas equações de 
(1) a (4): 
 
∇.⃗⃗⃗ B⃗⃗ = 0 (1) 
∇.⃗⃗⃗ D⃗⃗ = ρ (2) 
∇⃗⃗ × E⃗⃗ = −
∂B⃗⃗ 
∂t
 (3) 
∇⃗⃗ × H⃗⃗ = J +
∂D⃗⃗ 
∂t
 (4) 
 
Tais equações representam um sistema fechado, no qual as variáveis e equações 
possuem uma interdependência. A solução para o sistema deve ser obtida de maneira única. 
Analisando o termo ∂D⃗⃗ ∂t⁄ , tanto para casos estáticos quanto altas frequências, o mesmo 
assume valores muito baixos, próximos a 0,01% do valor da densidade de corrente J, 
podendo, assim, ser desprezado e eliminado da equação (4). Sendo assim, a equação (4) passa 
a ser representada pela equação (5): 
 
∇⃗⃗ × H⃗⃗ = J (5) 
 
 A equação (1) representa que o fluxo magnético, B⃗⃗ , é conservativo. A equação (2), 
Equação de Gauss, faz a ligação entre a passagem de cargas elétricas ρe com a geração de 
fluxo elétrico. A equação (3), Equação de Faraday, representa a ligação entre um campo 
elétrico e o campo magnético. A equação (4), Lei de Ampere, define a relação entre o campo 
magnético ao redor de um fluxo de corrente e de uma carga elétrica. 
26 
 
 
Para completar as equações de Maxwell, são adicionadas equações (6) a (8), que são 
relacionadas às propriedades dos materiais, chamadas de equações de passagem. 
 
D⃗⃗ = εE⃗⃗ (6) 
μ = 
∂B⃗⃗ 
∂�⃗⃗� 
 (7) 
j = σE⃗⃗ (8) 
 
Com a eliminação do termo ∂D⃗⃗ ∂t⁄ , pode-se realizar o desacoplamento entre as 
equações de Maxwell e obter dois sistemas de equações. O primeiro é composto pelas 
equações (1),(3), (5), (6) e (8). O segundo sistema é definido pelas equações (2) e (6). É 
importante ressaltar que os campos elétricos resultantes das equações (6) e (8) são de 
naturezas distintas. O resultado determinado pela equação (6) é consequência da variação 
temporal do campo magnético, já o proveniente da equação (8) é devido a presença de cargas 
estáticas. 
Dividindo a equação (3), pela equação (7) e aplicando o rotacional em ambos os lados 
obtém-se a equação (9). 
 
∇ ⃗⃗ ⃗ × (
1
𝜇
∇ ⃗⃗ ⃗ × �⃗� ) = −
∂
∂t
(∇ ⃗⃗ ⃗ × �⃗⃗� ) (9) 
 
Isolando D⃗⃗ na equação (6) e substituindo em (4) obtém-se a equação (10). 
 
∇⃗⃗ × H⃗⃗ = j +
∂εE⃗⃗ 
∂t
 (10) 
 
 Substituindo o rotacional da intensidade do campo magnético obtido pela equação (10) 
na equação (9), obtém-se a equação (11). 
 
∇ ⃗⃗ ⃗ × (
1
𝜇
∇ ⃗⃗ ⃗ × �⃗� ) +
∂
∂t
(J +
∂εE⃗⃗ 
∂t
) = 0 (11) 
 
A densidade de corrente é composta pela parcela induzida pelo campo elétrico, σE⃗⃗ , e 
pela parcela imposta pelo condutor, Js. Substituindo os valores no vetor J, obtém-se a equação 
(12). 
27 
 
 
 
∇ ⃗⃗ ⃗ × (
1
𝜇
∇ ⃗⃗ ⃗ × �⃗� ) + ε
∂2E⃗⃗ 
∂t2
+ 𝜎
∂E⃗⃗ 
∂t
= −
∂J𝑠
∂t
 (12) 
 
Para solução do problema é necessário adotar um perfil que atenda às características da 
densidade de corrente. De acordo com Ida e Bastos (1997), pode-se adotar que o perfil da 
densidade de corrente, para um campo com simetria axial, assume a forma senoidal de acordo 
com a equação (13). 
 
Js = J0𝑒
−z γ⁄ cos(2πf𝑡 − z γ⁄ ) (13) 
 
sendo γ, calculado por (14), a profundidade de penetração da corrente do indutor na peça, e z 
a diferença entre o indutor e a peça. 
 
γ = √
2
2𝜋𝑓𝜎𝜇
 (14) 
 
Por meio da equação (13), pode-se avaliar que quando z é igual a γ, o valor da 
densidade de corrente atinge 1/e, isto é, 37% da densidade de corrente máxima. Substituindo o 
valor de Js, equação (13), na equação (12), e derivando a densidade de corrente, obtém-se a 
equação (15), que representa o campo elétrico gerado um corpo, devido a presença de uma 
densidade de corrente. 
 
∇ ⃗⃗ ⃗ × (
1
𝜇
∇ ⃗⃗ ⃗ × �⃗� ) + ε
∂2E⃗⃗ 
∂t2
+ 𝜎
∂E⃗⃗ 
∂t
= − J02πf 𝑒
−z γ⁄ cos(2πf𝑡 − z γ⁄ ) (15) 
 
2.2 MATERIAIS MAGNÉTICOS 
 
A passagem de um campo magnético entre dois meios sofre uma variação angular, 
sendo definida pela equação (7). A propriedade μ é denominada permeabilidade magnética, e 
é específica de cada meio em que o campo magnético é aplicado. 
De forma a possuir um parâmetro de comparação, define-se como permeabilidade 
relativa, μr, a razão da permeabilidade real do meio μ pela permeabilidade do ar μ0, que 
corresponde a 4π.10-7 H/m. 
28 
 
Existem dois tipos de materiais magnéticos, os meios moles e os meios duros, que são 
os imãs permanentes. Analisando de forma detalhada os meios moles, os mesmos dividem-se 
em três categorias: 
 Diamagnéticos 
 Paramagnéticos 
 Ferromagnéticos 
Os materiais diamagnéticos possuem μr, pouco abaixo de 1, sendo constituídos de 
metais como o cobre, mercúrio e ouro. Para efeito de comparação, a permeabilidade relativa 
do cobre é 0,999991, ao passo que os materiais paramagnéticos também possuem μr próximos 
a 1, porém, ligeiramente superiores, podendo-se destacar o alumínio com μr = 1,00000036. 
Já os materiais ferromagnéticos possuem μr muito superior a 1, sendo que para algumas 
ligas de ferro esse valor chega a ser um milhão de vezes maior. Contudo, o comportamento 
magnético dos aços está atrelado à temperatura do metal. Após atingir a temperatura de Curie, 
770 °C, o material adquire um comportamento paramagnético. Contudo, o valor de μr de 
materiais ferromagnéticos, está atrelado ao campo magnético aplicado, fenômeno conhecido 
como saturação. Os resultados obtidos por Drobenko, Hachkevych e Koyrnyts’kyi (2007) 
demonstram a importância de analisar a variação das propriedades do material de acordo com 
a temperatura, sejam elas propriedades magnéticas ou térmicas. Contudo, podemos 
desconsiderar a variação da permeabilidade magnética para temperaturas abaixo de 600 °C, 
sem grande impacto no resultado para aços de baixo teor de carbono. Diversos estudos foram 
realizados para se prever o comportamento da permeabilidade do aço em função da 
temperatura, podendo também destacar, Skoczkowski, Kalus (1989), Ludtke e Schulze (1998) 
e Prakht et al. (2014). 
A penetração de campos variáveis em meios condutores trata-se de um problema 
complexo, e para simplificar a solução, considera-se o meio isotrópico e linear. 
 
2.3 FORMULAÇÃO TÉRMICA 
 
2.3.1 Condução 
 
A transferência de calor é um ramo da Ciência que estuda os mecanismos de 
transferência de energia entre os corpos. O foco desta linha de estudo é explicar como a 
29 
 
energia pode ser transferida e avaliar como ocorre esta interação entre os corpos em 
determinadas condições. 
A primeira linha de estudo na transferência de calor foi focada na condução, que existe 
quando ocorre um gradiente de temperatura no corpo, isto é, quando uma região do corpo 
sólido possui temperatura maior do que a outra, ocorre uma transferência de calor da região 
mais quente para a mais fria, como observa-se na Figura 11. A equação (16) representa esse 
fenômeno, sendo denominada equação de Fourier: 
 
𝑞 = −kA
𝜕T
𝜕𝑥
 (16) 
 
no qual o sinal negativoé inserido para respeitar a segunda lei da Termodinâmica, pois o calor 
deve fluir no sentido oposto ao gradiente de temperatura. A constante k da equação (16) é 
definida como condutividade térmica, e representa a resistência do corpo à passagem do calor. 
Quanto maior é o seu valor, mais facilmente o calor pode ser conduzido pelos corpos, isso 
ocorre no caso de sólidos metálicos. Já para a situação oposta, valores baixos de 
condutividade, é benéfica para o uso de isolantes, como gases e materiais cerâmicos. 
 
Figura 11: Condução de calor por uma parede de espessura ∆x e área com geração de 
calor G. 
 
 
Fonte: Produção do próprio autor. 
 
Realizando a análise de transferência de calor de um volume elementar, de espessura 
dx, a uma distância x, aplicando a equação (16). As equaões (17) e (18) representam o calor 
na entrada do elemento e na saida do elemento, respectivamente. A equação (19) representa a 
geração interna de calor no elemento e a equação (20) a variação de energia interna do corpo. 
30 
 
 
𝑞x = −kA
𝜕T
𝜕𝑥
 (17) 
𝑞x+dx = −kA
𝜕T
𝜕𝑥
|
x+dx
 (18) 
G = q̇Adx (19) 
∆U = ρcp
∂T
∂t
dx (20) 
 
Realizando um balanço de energia no elemento mostrado na Figura (11), obtem-se a 
equação (21), definida para um elemento unidmiensional: 
 
∂
∂x
(k 
𝜕T
𝜕𝑥
) + q ̇ = ρcp
∂T
∂t
 (21) 
 
Expandindo a relação (21) para as três dimensões, obtem-se a equação (22): 
 
∂
∂x
(k 
∂T
∂x
) +
∂
∂y
(k 
𝜕T
𝜕y
) + 
∂
∂z
(k 
𝜕T
𝜕z
) + q ̇ = ρcp
∂T
∂t
 (22) 
 
2.3.2 Convecção 
 
 O segundo mecanismo de transferência de calor estudado é a convecção, que pode ser 
definida como a tranferência de calor entre uma superfície sólida e fluido, gasoso ou líquido, 
adjacente a esta superfície, envolvendo o efeito combinado da condução com o movimento do 
fluido. Quanto mais rápido é o movimento do fluido, maior tenderá a ser a tranferência de 
calor. 
 De forma genérica, o efeito global da convecção pode ser representado pela lei de 
Newton do Resfriamento, conforme equação (23): 
 
qc = hA(T − T∞) (23) 
 
a qual o coeficiente h é denominado coeficiente de convecção, sendo o valor ligado ao tipo de 
superfície, velocidade, tipo de fluido, etc. Diversas equações podem ser aplicadas, variando 
para o tipo de escoamento, interno ou externo, se a convecção é forçada ou natural. Contudo, 
para algumas situações sua determinação deve ser feita por experimentos. Menores valores de 
31 
 
h são obtidos para fluidos gasosos com convecção natural, enquanto que valores elevados, são 
obtidos para movimentações forçadas de fluidos liquidos. 
 
2.3.3 Radiação 
 
 O terceiro tipo de mecanismo de troca de calor apresentado é a radiação. Distintamente 
dos outros tipos de troca de calor, a radiação não precisa de um meio fisico para que ocorra, 
pois o calor pode ser transferido de regiões em que o vácuo existe, sendo o principal exemplo 
o aquecimento efetuado pela radiação solar. 
A radiação é um fenômeno volumétrico em que todos os tipos de corpos podem emitir, 
absorver ou transmití-la. Entretanto, para realização de estudos, é considerada um fenômeno 
de superfície, pois a radiação proveniente do núcleo não pode atingir a superfície, e aquela 
incidente no corpo é absorvida por camadas próximas a superfície. O ponto de partida para 
estudos de radiação é o corpo negro, que é considerado o radiador ideal, pois emite radiação a 
maior taxa possível. O calor liberado por radiação de um corpo, pode ser calculado pela lei de 
Stephan-Boltzmann, conforme equação (24). 
 
qr = εσ𝑠𝑡𝑒A(T
4 − T∞
4) (24) 
 
na qual, σ𝑠𝑡𝑒 é a constante de Stephan-Boltzmann e ε é a emissividade da superfície do 
material. A emissividade éum parametro que varia de 0 ≤ ε ≤ 1, e representa a proximidade 
que uma superfície tem do corpo negro, pois ε do corpo negro vale 1. 
 Outra propriedade importante de um material é a absortividade, α, que simboliza a 
fração da energia incidente na superfície do material que é absorvida pelo corpo. Seus valores 
também variam de 0 a 1, sendo também o corpo negro o valor máximo, isto é, absorve o 
maior volume de energia incidente. 
Na superfície de um corpo podem ocorrer transferências de calor provenientes da 
convecção e da radiação, sendo assim o balanço de energia na superfície do corpo é 
demonstrado na equação (25). 
 
 div (k∇⃗⃗ T) = εσ𝑠𝑡𝑒A(T
4 − T∞
4) + hA(T − T∞) (25) 
 
 
32 
 
3 METODOLOGIA 
 
3.1 ELEMENTOS FINITOS 
 
Os métodos de elementos finitos surgiram como uma alternativa para solucionar 
problemas muito complexos que não possuem soluções triviais por meio de suas equações, 
sendo uma ferramenta numérica que determina soluções aproximadas. De acordo com Lewis, 
Nithiarasu e Seetharamu (2004), o método foi originalmente desenvolvido para análises de 
tensões em estruturas complexas, sendo estendido posteriormente para campos gerais da 
mecânica. O método tem recebido uma atenção considerável em Engenharia, devido a sua 
ampla aplicação em diversas áreas, uma vez que para certas situações a obtenção de soluções 
analíticas normalmente não são possíveis. 
Diversas técnicas para análises numéricas evoluíram ao longo dos anos, sendo que 
dentre elas se destacam o método de diferenças finitas, o método de volumes finitos, o método 
de volumes baseados em elementos e o método de elementos de contorno. 
Outro método amplamente aplicado é o de volumes finitos, que usa a formulação 
integral de um volume com finitas divisões para discretizar as equações. O método, 
normalmente, é utilizado para discretização de sistemas fluido dinâmicos. 
O método de elementos finitos considera que um corpo é formado por seções pequenas 
e seções interconectadas chamadas elementos, que possuem uma relação com as equações que 
definem o modelo estudado, reduzindo equações diferenciais ou complexas em sistemas de 
equações algébricas de análise mais simples. Assim, para prosseguir com o estudo, o material 
é discretizado, isto é, subdividido em vários pequenos elementos, reduzindo a análise do 
problema de infinitas incógnitas, para um número finito de incógnitas em pontos específicos 
denominados nós. 
De maneira genérica, pode-se afirmar que uma análise por elementos finitos segue as 
etapas ilustradas na Figura 12. O problema se inicia realizando um modelo físico, definindo-
se, assim, as variáveis a serem analisadas, as condições as quais o problema estará exposto e 
as equações a serem solucionadas. 
A segunda etapa consiste na discretização da região de solução em um conjunto de 
elementos finitos, sem sobreposição. Esta etapa é denominada como geração da malha, na 
qual os elementos gerados podem possuir a forma linear, para análises em 1D, triangular ou 
quadrática, para análises bidimensionais, e piramidal e prismática para 3D, conforme a Figura 
13. A Figura 14 ilustra os componentes do qual cada elemento é formado. Sendo o nó o ponto 
de análise de cada elemento, a fronteira é a superfície de cada elemento que está exposta às 
33 
 
condições de contorno impostas no problema analisado. A escolha do tipo de elemento e do 
seu grau estão relacionados à precisão da solução do problema. Para o caso de um problema 
2D, um elemento triangular possui uma maior precisão de resposta do que um elemento 
retangular, devido a melhor uniformização do espaço da malha, pois os erros de discretização 
em regiões não lineares do espaço analisado serão menores. Com relação ao grau, 
comparando um elemento triangular da Figura 13, pode-se afirmar que quanto maior o 
número de nós em cada triângulo, maior será a precisão da resposta, devido ao maior número 
de pontos a serem analisados em cada elemento. Entretanto, a escolha de elementos de grau 
maior, leva a programas com maior número de cálculos, elevando, assim, o tempo de 
resposta, pois no tipo triangularconvencional cada elemento corresponde a uma matriz 3x3, 
ao passo que para o tipo triangular cúbico cada elemento corresponde a uma matriz 9x9. 
 
Figura 12: Fluxograma de uma análise por elementos finitos. 
 
 
Fonte: Lewis , 2004, p. 40 
34 
 
6 Disponível em: 
http://www.enautica.pt/publico/professores/vfranco/Fundamentos_Metodo_Elementos_Finitos.pdf acesso em 
jan. de 2016 
Figura 13: Tipos de Malhas de elementos finitos. 
 
Fonte: Site Faculdade Superior Nautica Infante D. Henrique6 
 
 
 
Figura 14: Exemplo de malha para um elemento. 
 
 
Fonte: Lewis , 2004, p.40 
 
Quando se discretizam os elementos e aplicam as equações e condições de contorno no 
objeto de estudo, é possível obter as soluções aproximadas nos nós de cada um dos elementos 
analisados. As funções aplicadas pelo método são chamadas de funções de interpolação, e a 
35 
 
 
precisão da resposta obtida está ligada ao grau do elemento estudado, bem como ao número 
de nós utilizados na observação. 
Para análise de elementos triangulares bidimensionais, utiliza-se a representação da 
Figura 15, na qual se obtém funções lineares na variação das propriedades entre os vértices 
dos elementos. 
 
Figura 15: Representação de um elemento triangular de 1° grau. 
 
 
Fonte: Lewis, 2004, p. 49 
 
Analisando o caso para determinar a distribuição de temperatura no elemento, pode-se 
dizer que o valor da variável em um ponto qualquer do corpo obedece a uma variação 
conforme a equação (26): 
 
T(x, y) = α1 + α2x + α3y (26) 
 
sendo a variação da temperatura para cada um dos vértices do elemento definida pelas 
equações (27) a (29): 
 
Ti = α1 + α2xi + α3yi (27) 
Tj = α1 + α2xj + α3yj (28) 
Tk = α1 + α2xk + α3yk (29) 
 
A primeira informação que pode-se extrair de um elemento triangular, é que sua área 
pode ser obtida por meio do cálculo do determinante demonstrada na equação (30). 
36 
 
 
2A = det [
1 𝑥𝑖 𝑦𝑖
1 𝑥𝑗 𝑦𝑗
1 𝑥𝑘 𝑦𝑘
] (30) 
 
Sendo assim, para a determinação dos valores de  para definição da variação de 
temperatura em cada um dos nós do elemento, utiliza-se as equações (31) a (33). 
 
α1 =
1
2A
[(xjyk − xkyk)Ti + (xkyi − xiyk)Tj + (xiyj − xjyi)Tk] (31) 
α2 =
1
2A
[(yj − yk)Ti + (yk − yi)Tj + (yi − yj)Tk] (32) 
α3 =
1
2A
[(xk − xj)Ti + (xi − xk)Tj + (xj − xi)Tk] (33) 
 
Simplificando as equações (31) a (33), e escrevendo-as de uma maneira matricial 
define-se a equação (34). 
 
T = [Ni Nj Nk] [
Ti
Tj
Tk
] (34) 
 
Sendo que os fatores N representam o fator de forma do elemento utilizado para descrever a 
maneira como a propriedade analisada irá variar entre os nós, definido também como função 
de interpolação. 
Detalhando os fatores de forma, pode-se dizer que, para elementos triangulares, cada 
um deles assume os valores descritos pelas equações (35) a (37). 
 
Ni = 
1
2A
 (ai + bix + ciy) (35) 
Nj = 
1
2A
 (aj + bjx + cjy) (36) 
Nk = 
1
2A
 (ak + bkx + cky) (37) 
 
Com os valores de a, b e c para os nós i, j e k, sendo calculados pelo conjunto de 
equações (38). 
 
ai = xjyk − xkyk bi = yj − yk ci = xk − xj 
37 
 
aj = xkyi − xiyk bj = yk − yi cj = xi − xj (38) 
ak = xiyj − xjyi bk = yi − yj ck = xj − xi 
 
Observa-se que a soma dos fatores N nunca deve sempre ser igual a unidade, conforme 
é demonstrado pela equação (39). 
 
𝑁𝑖 + 𝑁𝑗 + 𝑁𝑘 = 1 (39) 
 
O gradiente de temperatura T de um elemento triangular pode ser definido pelas 
equações (40) e (41): 
 
∂T
∂x
= 
∂Ni
∂x
Ti +
∂Nj
∂x
Tj +
∂Nk
∂x
Tk = 
bi
2A
Ti +
bj
2A
Tj + 
bk
2A
Tk (40) 
∂T
∂y
= 
∂Ni
∂y
Ti +
∂Nj
∂y
Tj +
∂Nk
∂y
Tk = 
ci
2A
Ti +
cj
2A
Tj + 
ck
2A
Tk (41) 
 
Para realizar integrais associadas aos elementos de forma, utiliza-se as equações (42) e 
(43). 
 
∫ Ni
a
A
Nj
bNk
cdA = 
a!b!c!
(a+b+c+2)!
 2A (42) 
 ∫ Ni
a
l
Nj
bdl = 
a!b!
(a+b+1)!
 l (43) 
 
A equação (42) corresponde as integrais realizadas no volume Ω e a equação (43) 
corresponde às integrais na superfície Γ. O volume Ω e a superfície Γ estão exemplificados na 
Figura 16. 
Como dito anteriormente, um elemento bidimensional pode assumir outros formatos. A 
escolha depende da geometria da peça e do nível de precisão que se deseja conseguir na 
solução de problema. Logo, a função de forma irá variar conjuntamente com o elemento. 
Na próxima etapa da formulação do problema, faz-se a ligação entre as incógnitas e as 
cargas aplicadas ao elemento, montando as matrizes referentes aos elementos a serem 
estudados e resolvendo esse sistema de equações. Para a análise de Métodos de Elementos 
Finitos de um sistema térmico, utiliza-se a equação de Fourrier, representada pela equação 
(22), na sua forma matricial, representada pela equação (44): 
 
38 
 
[𝐊]{𝐓} = {𝐟} (44) 
 
Figura 16: Representação genérica de um elemento Finito. 
 
 
Fonte: Zienkiewicz, 2000, p.40 
 
sendo [K] a matriz de rigidez do elemento relacionado a condução de calor e as perdas por 
convecção, {T} o vetor de incógnita da temperatura dos nós e {f} as cargas aplicadas ao 
elemento proveniente das perdas convectivas ou fluxos de calor, geração interna de energia do 
problema estudado. 
No estudo de elementos finitos, várias são as técnicas desenvolvidas, dentre elas, pode-
se destacar os seguintes métodos: 
 
 Método de Ritz 
 Método de Ritz Rayleigh – Método Variacional 
 Método dos Resíduos Ponderados 
 
De forma a demonstrar o desenvolvimento dos três métodos, utiliza-se um exemplo em 
comum, descrito por Lewis (2004) no qual uma aleta demonstrada na Figura 17, tem como 
condições de contorno a temperatura de base definida e constante, e sua ponta isolada. 
 
 
 
 
39 
 
Figura 17: Representação de uma aleta e suas condições de contorno 
 
 
Fonte: Lewis, 2000, p. 75 
 
Aplicando a equação de Fourrier, na forma diferencial, no corpo da aleta, tem-se: 
 
−kA
dT
dx
|
x
= hPdx(T − Ta) − kA
dT
dx
|
x+dx
 (45) 
kA
d2T
dx2
= hP(T − Ta) (46) 
 
As condições de contorno da equação são dadas pelas equações (47) e (48), que 
representam a situação de ponta isolada e a temperatura de base constante, respectivamente: 
 
𝑑𝑇
𝑑𝑥
|
𝑥=0
= 0 (47) 
𝑇(𝑥 = 𝐿) = 𝑇𝑏 (48) 
 
Adimensionalizando a equação (45) com os parâmetros, 𝑇 − 𝑇𝑎 = 𝜃, 𝜁 = 𝑥/𝐿 e 
ℎ𝑃/𝑘𝐴 = 𝑚², e substituindo nas incógnitas em (46): 
 
d2θ(ζ)
dζ2
− μ2θ(ζ) = 0 (49) 
 
 
Como primeiro exemplo de aplicação de Elementos Finitos, será aplicado o método de 
Ritz para estudar a variação de temperatura na aleta. 
40 
 
A solução aproximada para a equação (49) pode ser obtida por meio da solução da 
seguinte série, expressa na equação (50): 
 
𝑇 ≈ 𝑇 ̅ = ∑ 𝑎𝑖𝑁𝑖(𝑥)
𝑛
𝑖=1 (50) 
 
Quanto maior o número de elementos n, maior é a precisão da solução. Porém um 
resíduo R é deixado na aproximação. A melhor solução é obtida quando o valor de R tende à 
zero em todos os pontos do domínio. 
Adotando um perfil que satisfaça as condições de contorno da equação (49), obtém-se 
que a variação de temperatura da aleta pode ser representada pela equação (51), sendo B a 
incógnita a ser determinada na solução da equação: 
 
θ(ζ) = θb{1 − (1 − ζ
2)B} (51) 
 
Para solucionar a equação pelo método de Ritz, substitui-se o perfil adotado, equação 
(51), na equação (49) e o resíduo R é integrado no domínio e igualado a zero, para assim 
determinar o valor de B, 
 
∫ (
d2θ(ζ)
dζ
− μ2θ)
1
0
dζ = 0 (52) 
 
Diferenciando duas vezes o perfil de temperatura adotado na equação (51), na equação 
(49): 
 
d2θ(ζ)
dζ
= 2θbB (53) 
 
Substituindo (53) na integral da equação (52) e desenvolvendo matematicamente: 
 
∫ {2θbB − μ
2θb[1 − (1 − ζ
2)B]} dζ
1
0
= 0 (54)2θbBζ − μ
2θb [ζ − (ζ −
ζ3
3
) B]|
0
1
= 0 
 
2θbB − μ
2θb[1 − 2/3B] = 0 
41 
 
B = 
μ2
2
1+
μ2
3
 
 (55) 
 
Substituindo o valor de B na equação (49), determina-se o seguinte perfil de variação de 
temperatura: 
 
θ(ζ) = θb {1 − (1 − ζ
2) (
μ2
2
1+
μ2
3
 
)} (56) 
 
Sendo a equação que descreve a variação de temperatura de uma aleta, obtida por meio 
da aproximação da função por funções hiperbólicas, representada pela equação (57): 
 
θ(ζ) = θb [
coshm(L−x)
coshmL
] (57)
 
De forma a comparar com valores numéricos obtidos com as equações (29) e (30), 
adota-se que a aleta em estudo seja fabricada em aço inox (k = 16,66 W/mK), diâmetro de 2 
cm, comprimento de 10 cm e esteja sujeita a uma convecção de h= 25 W/m²K, sendo 
possível estruturar o Gráfico 1. 
 
Gráfico 1: Comparação entre a Temperatura da Aleta pelo Método de Ritz e equação 
hiperbólica. 
 
 
Fonte: Produção do próprio autor. 
42 
 
 
Pode-se notar que a amplitude dos valores obtidos é maior quando estamos próximos a 
base da aleta, mas diminui próximo a extremidade externa da peça. 
O segundo método estudado é o de Rayleigh - Ritz, que é uma variação do cálculo 
variacional. A princípio este método não possui nenhuma relação com o método de elementos 
finitos, por tratar-se de um método matemático. O mesmo se baseia na minimização de um 
funcional energético, no qual uma função f(x) que no extremo da variação da integral 
corresponde à equação diferencial, é solução da equação diferencial original nas condições de 
contorno, pois a minimização de um funcional energético, corresponde ao estado de equilíbrio 
de todo sistema natural livre, que possui certa energia potencial. 
Por meio da aplicação da equação diferencial de Euller-Lagrange na equação (49), 
determina-se a equação (58): 
 
δI = ∫ (
d2θ(ζ)
dζ2
− μ2θ(ζ))
1
0
δθdζ (58) 
 
Integrando a equação (59) por partes, 
 
[
dθ
dζ
δθ]
0
1
− ∫ (
dθ
dζ
)
d
dζ
(δθ)dζ
1
0
− μ2 ∫ θδθdζ
1
0
= 0 (59) 
 
E usando as relações 
𝑑
𝑑𝜁
(𝛿𝜃) = 𝛿 (
𝑑𝜃
𝑑𝜁
), (
𝑑𝜃
𝑑𝜁
) 𝛿 (
𝑑𝜃
𝑑𝜁
) = 
1
2
𝛿 (
𝑑𝜃
𝑑𝜁
)
2
, pode-se simplificar a 
equação (59), obtendo a equação admensionalizada (60): 
 
[
dθ
dζ
δθ]
0
1
−
1
2
δ∫ [(
dθ
dζ
)
2
+ μ2θ2]
1
0
dζ = 0 (60) 
 
Aplicando as condições de contorno, idênticas a aplicadas no Método de Ritz, a 
formulação variacional do problema é: 
 
I = 
1
2
δ∫ [(
dθ
dζ
)
2
+ μ2θ2]
1
0
dζ (61) 
 
43 
 
O perfil que minimiza a integral da equação (61) é a solução da equação diferencial 
principal. Assumindo o mesmo perfil, equação (51), e substituindo na equação (61), o 
variacional corresponde a: 
 
I = 
1
2
∫ θb
2{(2Bζ)2 + μ2[1 − (1 − ζ2)B2] }
1
0
dζ (62) 
 
Integrando (62) e substituindo os limites da integral, 
 
I = 
1
2
θb [B
2 (
4
3
+ μ2 −
2
3
μ2 +
1
5
μ2) + μ2 + B(−2μ2 +
2
3
μ2)] (63) 
 
Derivando o potencial I, equação (63), em relação a B e igualando a 0, 
 
∂I
∂B
= 
1
2
θb [2B (
4
3
+ μ2 −
2
3
μ2 +
1
5
μ2) + μ2 + (−2μ2 +
2
3
μ2)] = 0 (64) 
 
Simplificando (64), 
 
2B (
4
3
+
8
15
μ2) = 
4
3
μ2 (65) 
 
Isolando B em (65) e substituindo na equação (22), determina-se a representação do 
perfil pela equação (39). 
 
𝜃(𝜁) = 𝜃𝑏 {1 − (1 − 𝜁
2) (
𝜇2
2
1+
2𝜇2
5
 
)} (66) 
 
Realizando a mesma suposição de material de aleta e condições de contorno do Método 
de Ritz para a solução do problema pelo variacional, efetua-se a comparação entre a Equação 
Hiperbólica e o Método Variacional, representada pelo Gráfico 2. 
É possivel observar que com a aplicação do método variacional, os valores de 
temperatura obtidos por (66) estão mais próximos do valor obtido pela equação (56), 
mostrando assim uma maior precisão do método. 
 
 
44 
 
Gráfico 2: Comparação entre a Temperatura da Aleta pelo Método Variacional e 
equação hiperbólica 
 
 
Fonte: Produção do próprio autor. 
 
Para as equações em que não é possível a aplicação do método variacional, utiliza-se 
uma formulação alternativa. O método dos resíduos ponderados fornece uma solução 
aproximada que é aplicável a um número muito elevado de problemas, sendo desnecessária a 
determinação de formulação de variacionais para a aplicação de métodos de elementos finitos 
para os problemas. 
Sendo a equação (67) a representação da variação da temperatura no volume Ω: 
 
L(T) = 0 (67) 
 
E a temperatura do corpo determinada pela equação (50), substituindo em (67) T por T̅ : 
 
L(T̅) = R (68) 
 
no qual R representa os resíduos gerados pela aproximação efetuada. 
O método dos resíduos necessita que os índices 𝑎𝑖 da equação (50) sejam calculados da 
seguinte maneira: 
 
∫ wi(x)Rdx = 0Ω (69) 
45 
 
Sendo que wi corresponde a funções de ponderação arbitrárias, e a escolha da função 
depende do método a ser aplicado. Dentre os métodos de resíduos ponderados, pode-se 
destacar: 
 
Colocação: wi = δ(x − xi) 
 
∫ R δ(x − xi)dx = 0Ω (70) 
 
Subdomínio: wi = 1 
 
∫ 𝑅 𝑑𝑥 = 0
Ω
 (71) 
 
Galerkin: wi(x) = Ni(x) 
 
∫ Ni(x)Rdx = 0Ω (72) 
 
Mínimos Quadrados: wi = ∂R ∂ai⁄ 
 
∫ 𝜕𝑅 𝜕𝑎𝑖⁄ 𝑅𝑑𝑥 = 0Ω (73) 
 
Para diferenciar a formulação de cada um dos métodos, será aplicada cada uma das 
teorias ao problema da aleta 
O método da Colocação satisfaz a equação diferencial em pontos específicos, sendo 
baseado em polinômios definidos por trechos. O problema costuma ter solução mais fácil que 
o Método de Galerkin, porém, sua solução pode não ser tão eficiente para problemas de 
Engenharia. 
Supondo a seguinte função ponderada wi = δ(x − xi), onde ζi = 1/2, sendo que a 
equação para aleta é: 
 
∫ [
𝑑2𝜃(𝜁)
𝑑𝜁2
− 𝜇2𝜃(𝜁)] 𝛿(𝜁𝑖 − 𝜁)𝑑𝜁 = 0
1
0
 (74) 
 
46 
 
[
𝑑2𝜃(𝜁)
𝑑𝜁2
− 𝜇2𝜃(𝜁)]
𝜁𝑖=1/2
= 0 (75) 
 
Integrando a equação (75), define-se o seguinte valor de B: 
B = 
μ2
2
1+
3μ2
8
 
 (76) 
 
O método do subdomínio, parte da suposição que a função ponderada 𝑤𝑖 = 1, sendo 
assim a formulação assume valor idêntico ao do método de Ritz, equação (55): 
 
∫ [
𝑑2𝜃(𝜁)
𝑑𝜁2
− 𝜇2𝜃(𝜁)] 𝑑𝜁 = 0
1
0
 (53) 
 
Sendo o valor de B igual a: 
B = 
μ2
2
1+
μ2
3
 
 (55) 
 
O método de Galerkin adota que a função ponderada Ni(x) = (1 − ζ
2), pode ser 
adotada para o resíduo R: 
 
∫ [
d2θ(ζ)
dζ2
− μ2θ(ζ)] (1 − ζ2)dζ = 0
1
0
 (77) 
 
Substituindo o perfil de B, equação (52), e integrando a equação (77), define-se o 
seguinte valor de B: 
 
B = 
μ2 2⁄
1+
2μ2
5
 
 (78) 
 
É possível afirmar que o valor de B obtido pela equação (78) é o mesmo que o definido 
pelo método variacional. Isto ocorre porque o problema analisado possui um comportamento 
variacional. 
 
A análise pelo método dos mínimos quadrados, adota que a função ponderada wi =
 ∂R ∂ai⁄ , e uma função de erro, definida pela equação (79): 
47 
 
 
𝐸 = ∫ 𝑅2𝑑𝜁
1
0
 (79) 
∫ [
𝑑2𝜃(𝜁)
𝑑𝜁2
− 𝜇2𝜃(𝜁)]
2
𝑑𝜁 = 0
1
0
 (80) 
 
Substituindo o perfil adotado para a variação de temperatura da aleta, equação (51), e 
integrando a equação (80), obtém-se B: 
 
B = 
μ2
2
(1+
μ2
3
)
1+2μ2 (
1
3
+
μ2
15
) 
 (81) 
 
Realizando a mesma suposição de material de aleta e condições de contorno do 
problema, obtém-se a comparação entre os quatro métodos de resíduos, para cada um dos 
valores específicos de B encontrado. Os valores estão representados no Gráfico 3. 
 
Gráfico 3: Comparação entre a Temperatura da Aleta pelo equação hiperbólica e os 
métodos dos resíduos. 
 
Fonte: Produção do próprio autor. 
 
É possível afirmar por meio da análise do Gráfico 3, que dentre os métodos aplicados, 
aquele que apresenta os menores desvios é o método de Galerkin. 
 
3.2MÉTODO DE GALERKIN 
48 
 
De forma a detalhar a aplicação de Elementos Finitos, tanto pelo método variacional 
quanto por Galerkin, será considerada a equação de Fourrier, aplicada de maneira genérica: 
 
∇⃗⃗ (k∇T) = Q (82) 
 
Em que as condições de contorno são: 
 
T = Tb em Γ12 
 
kx
∂T
∂x
l̂ + ky
∂T
∂y
m̂ + kz
∂T
∂z
n̂ + q = 0 em Γ23 
 
kx
∂T
∂x
l̂ + ky
∂T
∂y
m ̂ + kz
∂T
∂z
n̂ + h(T − Ta) = 0 em Γ31 
 
sendo que os índices 12, 23 e 31 representam as arestas de um elemento triangular. Por meio 
da Figura 14, pode-se observar algumas características de uma análise de elemento finito, 
sendo a superfície que limita o elemento é representada por Γ e o volume que o ele possui é 
representado por Ω. 
A integral do método variacional aplicada à equação (82), com as condições de 
contorno, é definida pela equação (83): 
 
I(T) = 
1
2
[∫ kx (
∂T
∂x
)
2
+ ky (
∂T
∂y
)
2
+ kz (
∂T
∂z
)
2
Ω
− 2GT] dΩ + ∫ qTdΓ
Γ23
+
 ∫
1
2
h(T − Ta)
2dΓ
Γ31
 (83) 
 
O domínio Ω analisado é subdivido em n elementos, com cada um deles possuindo r 
nós, sendo a temperatura de cada elemento expressa pela equação (84): 
 
Te = ∑ NiTi
r
i=1 = [𝐍]{𝐓} (84) 
 
Em que [N] = [Ni, Nj,.....,Nr] é a matriz dos fatores de forma e {T} é o vetor de 
temperatura dos nós. 
A solução do problema consiste na seleção de valores de temperatura T dos nós para o 
qual a função I (T), seja fixa, satisfazendo a equação (85). 
49 
 
 
δI(T) = ∑
∂I
∂Ti
n
i=1 = 0 (85) 
 
Com n o número total de nós no volume discretizado no domínio da solução. 
Sendo Ti arbitrário, a solução só é possível se, for satisfeita a equação (86): 
 
δIe
δTi
= 0 para i = 0,1,2, . . . , n (86) 
 
A função I(T) pode ser descrita como a soma de funções individuais, definida para a 
montagem do elemento, apenas se as distintas funções utilizadas para a representação de T 
obedecerem um certa continuidade e compatibilidade, de acordo com a equação (87). 
 
I(T)= ∑ Ie(Te)ne=1 (87) 
 
Ao invés de trabalhar com uma função definida em todo domínio, a opção é trabalhar na 
definição de funções individuais para os elementos: 
 
δI = ∑ δIene=1 = 0 (88) 
 
Sendo a função Ie utilizada apenas com respeito aos r nós associados aos elementos e, 
para 
δIe
δT
=
δIe
δTj
 = 0 para j = 1,2, … , r 
A equação compreende um conjunto de r equações que caracteriza o comportamento 
dos e elementos. Para facilitar podemos representar a função de montagem dos elementos 
como o somatório de todas as funções que representam cada um dos elementos, do principio 
do variacional, conforme a equação (89). 
 
δI
δTi
= ∑
δIe
δTi
n
e=1 = 0 com i = 1,2, … ,M (89) 
 
A formulação é completa quando um conjunto M de equações é resolvido, 
simultaneamente, para os M valores de T, segundo a equação (90).. 
 
50 
 
I(T) = 
1
2
[∫ kx (
∂Te
∂x
)
2
+ ky (
∂Te
∂y
)
2
+ kz (
∂Te
∂z
)
2
Ω
− 2GTe] dΩ + ∫ qTedΓ
Γ23
 +
 ∫
1
2
h(Te − Ta)
2dΓ
Γ31
 (90) 
 
Com , 
𝑇𝑒 = [𝑵]{𝑻} = {𝑁1, 𝑁2, … , 𝑁𝑟} {
𝑇1
𝑇2
⋮
𝑇𝑟
} = 𝑁1𝑇1 + 𝑁2𝑇2 +⋯+𝑁𝑟𝑇𝑟 (91) 
 
sendo 
 
δTe
δT1
= N1 
δTe
δT2
= N2 
δTe
δ{𝐓}
= [
N1
N2
⋮
Nr
] = [𝐍]T (92) 
 
A matriz gradiente é definida pela equação (93): 
 
{𝐠} = 
{
 
 
 
 
∂Te
∂x
∂Te
∂y
∂Te
∂z }
 
 
 
 
= 
[
 
 
 
 
∂N1
∂x
∂N2
∂x
…
∂N1
∂y
∂N2
∂y
…
∂N1
∂z
∂N2
∂z
…
∂Nr
∂x
∂Nr
∂y
∂Nr
∂z ]
 
 
 
 
 {
T1
T2
T3
} = [𝐁]{𝐓} (93) 
{𝐠}T[𝐃]{𝐠} = {
∂Te
∂x
∂Te
∂y
∂Te
∂z
} [
kx 0 0
0 ky 0
0 0 kz
]
{
 
 
 
 
∂Te
∂x
∂Te
∂y
∂Te
∂z }
 
 
 
 
 (94) 
 
Substituindo a equação (94) na equação (89), para que os termos relacionados à 
condução do metal adquiram a forma matricial: 
 
I(T) = 
1
2
[∫ {𝐠}T[𝐃]{𝐠}
Ω
− 2GTe] dΩ + ∫ qTdΓ
Γ23
+ ∫
1
2
h(T − Ta)
2dΓ
Γ31
 (95) 
 
51 
 
Derivando (95), obtemos: 
 
δIe
δ{𝐓}
=
1
2
[∫ [𝐁]T[𝐃][𝐁]{𝐓}
Ω
− 2G[𝐍]T{𝐓}] dΩ + ∫ q[𝐍]T{𝐓}dΓ
Γ23
 +
 ∫
1
2
h[𝐍]T{𝐓}dΓ − ∫
1
2
h[𝐍]TTadΓΓ31
= 0
Γ31
 (96) 
 
Sendo a equação na forma matricial, definida pela equação (44): 
 
[𝐊]{𝐓} = {𝐟} (44)
 
Onde 
 
[𝐊] = ∫ [𝐁]T[𝐃][𝐁]
Ω
dΩ + ∫ h[𝐍]T[𝐍]dΓ = 0
Γ31
 (97) 
{𝐟} = ∫ G[𝐍]T
Ω
dΩ − ∫ q[𝐍]TdΓ
Γ23
+ ∫ h[𝐍]TTadΓΓ31
 (98) 
 
Seguindo com as deduções e avaliando a equação de Fourier, (82), aplica-se o Método 
de Galerkin sobre ela. Para utilização do método, deve-se considerar que: 
 
∫ wk𝐋(T̅)Ω dΩ = 0 (99) 
 
wk é substituído nos nós pela função de forma Nk(x), isto é: 
 
∫ Nk {
∂
∂x
kx
∂T̅
∂x
+
∂
∂y
ky
∂T̅
∂y
+
∂
∂z
kz
∂T̅
∂z
+ G}
Ω
dΩ (100) 
 
Integrando por partes e aplicando o teorema de Green na equação (100), obtém-se: 
 
∫ Nk {
∂
∂x
kx
∂T̅
∂x
}
Ω
dΩ = ∫ Nk (kx
∂T̅
∂x
) dΓ
Γ
− ∫
∂Nk
∂x
kx
∂Nm
∂xΩ
{T̅m}dΩ (101) 
 
Aplicando as mesmas condições de contorno utilizadas para o problema do variacional, 
na equação (101): 
 
52 
 
−∫ (kx
∂Nk
∂x
∂Nm
∂x
+ ky
∂Nk
∂y
∂Nm
∂y
+ kz
∂Nk
∂z
∂Nm
∂z
)
Ω
{T̅m}dΩ + ∫ GNkdΩ Ω − ∫ qNkdΓΓ +
 ∫ hNkNm{T̅m}dΓΓ + ∫ hNkTadΓΓ (102) 
 
Escrevendo a equação (102) na forma matricial, tem-se: 
 
[𝐊]{�̅�} = {𝐟} (103)
 
Sendo as matrizes [K] e {T} da equação (103) descritas da seguinte maneira: 
 
[𝐊] = −∫ (kx
∂Nk
∂x
∂Nm
∂x
+ ky
∂Nk
∂y
∂Nm
∂y
+ kz
∂Nk
∂z
∂Nm
∂z
)
Ω
dΩ + ∫ hNkNmdΓΓ (104) 
{𝐟} = ∫ GNkdΩ Ω − ∫ qNkdΓΓ + ∫ hNkTadΓΓ (105) 
 
Toda a modelagem até a presente etapa do trabalho foi feita baseada em elementos 3D. 
Entretanto, para elementos que possuem simetria axial, conforme mostra a Figura 18, esta 
análise pode se feita usando elementos em 2D. A formulação de Galerkin para este tipo de 
elemento é a mesma, porém os termos diferem de maneira significativa. 
Seja a equação de Fourier em coordenadas cilíndricas, representada pela equação (106): 
da seguinte forma: 
 
kr
∂2T
∂r2
+
kr
r
∂T
∂r
+ 
kθ
r2
∂2T
∂θ2
+ kz
∂2T
∂z2
+ G = 0 (106) 
 
Como um problema com simetria axial é independente do ângulo 𝜃, elimina-se as 
variáveis relacionadas à posição angular, e obtém-se a equação (107): 
 
kr
∂2T
∂r2
+
kr
r
∂T
∂r
+ kz
∂2T
∂z2
+ G = 0 (107) 
kr
r
∂
∂r
(
∂T
∂r
) + kz
∂2T
∂z2
+ G = 0 (108) 
 
 
 
 
 
53 
 
Figura 18: Elemento com simetria axial. 
 
 
Fonte: Zienkiewicz, 2000, p. 113 
 
As condições de contorno associadas ao problema são: 
 
T = Tb em Γ1 
kr
∂T
∂r
l + kz
∂T
∂z
n + h(T − Ta) + q = 0 
 
Com os valores de a, b e c sendo calculados pelo conjunto de equações (109): 
 
ai = rjzk − rkzk bi = zj − zk ci = rk − rj 
aj = rkzi − rizk bj = zk − zi cj = ri − rj (109) 
ak = rizj − rjzi bk = zi − zj ck = rj − ri 
 
Sendo a área de um elemento triangular com simetria axial calculada pela fórmula: 
 
2A = det [
1 ri zi
1 rj zj
1 rk zk
] (110) 
 
A integral de resíduos ponderados para um material com simetria axial, no volume, se 
torna: 
 
54 
 
−∫ Nk [
kr
r
∂
∂r
(
∂T
∂r
) + kz
∂2T
∂z2
+ G]
Ω
dΩ = R (111) 
 
Novamente, os termos das derivadas, devem ser transformados em termos de ordem 
menor, usando o teorema de Gauss e a regra do produto aplicada a derivadas: 
 
∂
∂z
(Nk
∂T
∂z
) = Nk
∂2T
∂z2
+ 
∂Nk
∂z
∂T
∂z
 (112) 
 
Reposicionando os termos da equação (112): 
 
Nk
∂2T
∂z2
=
∂
∂z
(Nk
∂T
∂z
) − 
∂Nk
∂z
∂T
∂z
 (113) 
 
Aplicando a regra do produto aos termos relativos ao raio: 
 
1
r
∂
∂r
(Nkr
∂T
∂r
) = 
1
r
(
∂Nk
∂r
r
∂T
∂r
) + 
1
r
Nk
∂
∂r
(r
∂T
∂r
) (114) 
 
Reposicionando os termos da equação (114): 
 
 Nk [
1
r
∂
∂r
(r
∂T
∂r
)] = 
1
r
∂
∂r
(Nkr
∂T
∂r
) − (
∂Nk
∂r
∂T
∂r
) (115) 
 
Substituindo as equações (115) e (113), na equação(112), 
 
∫ [kz
∂Nk
∂z
∂T
∂z
+ kr (
∂Nk
∂r
∂T
∂r
) − NkG]Ω dΩ − ∫ [ kz
∂
∂z
(Nk
∂T
∂z
) + 
kr
r
∂
∂r
(Nkr
∂T
∂r
)]
Ω
= R
 (116) 
 
A segunda integral pode ser transforma em uma integral de superfície usando o 
Teorema de Gauss: 
 
∫ [
kr
r
∂
∂r
(Nkr
∂T
∂r
) cosθ + kz
∂
∂z
(Nk
∂T
∂z
) senθ]
Γ
dΓ (117) 
 
Agrupando os termos da equação (117) e substituindo na equação (116): 
 
55 
 
∫ [kz
∂Nk
∂z
∂T
∂z
+ kr (
∂Nk
∂r
∂T
∂r
) − NkG]Ω dΩ − ∫ Nk [ 
kr
r
∂
∂r
(r
∂T
∂r
) cosθ +
Ω
kz
∂
∂z
(
∂T
∂z
) senθ] dΓ = R (118) 
 
Aplicando as condições de contorno a segunda integral da equação (118), produz-se: 
 
∫ [kz
∂Nk
∂z
∂T
∂z
+ kr (
∂Nk
∂r
∂T
∂r
) − NkG]Ω dΩ − ∫ Nk[ hN + hTa + Q]Ω dΓ = R (119) 
 
Aplicando a relação 
∂T
∂z
= 
∂N
∂z
 {T} e 
∂T
∂r
= 
∂N
∂r
 {T}, 
 
∫ [kz
∂Nk
∂z
∂N
∂z
+ kr (
∂Nk
∂r
∂N
∂r
)] {T}
Ω
dΩ + ∫ NkGΩ dΩ − ∫ Nk[ hN + hTa + Q]Ω dΓ = R
 (120) 
 
E pode-se afirmar que na forma matricial a equação (120) é representada da seguinte 
conforme equação (44): 
 
[𝐊]{𝐓} = {𝐟} (44)
 
Onde 
 
[𝐊] = ∫ [𝐁]T[𝐃][𝐁]
Ω
dΩ + ∫ h[𝐍]T[𝐍]dΓ = 0
Γ
 (121) 
{𝐟} = ∫ G[𝐍]T
Ω
dΩ − ∫ q[𝐍]TdΓ
Γ
+ ∫ h[𝐍]TTadΓΓ (122) 
 
Para o qual são utilizadas as matrizes representadas pelas equações (123) e (124): 
 
[B] = 
1
2A 
[
bi bj bk
ci cj ck
] (123) 
[B]T = 
1
2A 
[
bi ci
bj cj
bk ck
] (124) 
[D] = [
kr 0
0 kz
] (125) 
56 
 
[N]T = [
Ni
Nj
Nk
] (126) 
[N] = [Ni Nj Nk] (127) 
 
Por se tratar de elementos com simetria axial, a determinação da integral dos parâmetros 
de volume Ω e superfície Γ para utilização das equações (42) e (43) assumem a seguinte 
forma: 
 
dΩ = 2πr̅dA (128) 
dΓ = 2πr̅dl (129) 
 
Para a resolução de problemas transientes é necessário à aplicação do método de 
Galerkin em uma equação transiente. Para um sólido a equação (130) é utilizada: 
 
∂
∂x
(kx
∂T
∂x
) +
∂
∂y
(ky
∂T
∂y
) +
∂
∂z
(kz
∂T
∂z
) + G = ρcp
∂T
∂t
 (130) 
 
Na equação é considerado a anisiotropia das propriedades, e a incógnita a ser 
determinada é a Temperatura no tempo. 
As seguintes condições de contorno serão aplicadas: 
 
T = Tb na superfície Γb 
kx
∂T
∂x
l̂ + ky
∂T
∂y
m ̂ + kz
∂T
∂z
n̂ + q + h(T − Ta) = 0 em Γe 
T = T0 no instante t = 0 
 
Para aplicação do método de Galerkin, a seguinte discretização da temperatura é 
efetuada: 
 
T(x, y, z, t) = ∑ Ni(x, y, z)Ti(t)
n
i=1 (131) 
 
A representação da equação acima por Elementos Finitos se torna, 
 
57 
 
−∫ [kx
∂Ni
∂x
∂T
∂x
+ ky
∂Ni
∂y
∂T
∂y
+ kz
∂Ni
∂z
∂T
∂z
− NiG + Niρcp
∂T
∂t
] dΩ
Ω
− ∫ NiqdΓΓ − ∫ Nih(T −Γ
Ta)dΓ (132) 
 
Aplicando a relação abaixo, obtém-se uma equação de forma mais simplificada: 
 
∂T
∂x
= 
∂N
∂x
 {T} 
 
−∫ [kx
∂Ni
∂x
∂Nj
∂x
Tj(t) + ky
∂Ni
∂y
∂Nj
∂y
Tj(t) + kz
∂Ni
∂z
∂Nj
∂z
Tj(t) − NiG + Niρcp
∂Nj
∂t
Tj(t)] dΩΩ −
 ∫ NiqdΓΓ − ∫ Nih(T − Ta)dΓΓ (133)
 
Com o agrupamento de (133) a equação em sua forma matricial é: 
 
[Ct] {
∂T
∂t
} + [Kt]{T(t)} − {f} = 0 (134) 
 
Na qual os termos são representados da seguinte maneira 
 
[Ct] = [Cij] = ∫ ρcpNiNjdΩΩ (135) 
 
[Kt] = [Kij] = ∫ [kx
∂Ni
∂x
∂Nj
∂x
+ ky
∂Ni
∂y
∂Nj
∂y
+ kz
∂Ni
∂z
∂Nj
∂z
] dΩ
Ω
+ ∫ hNiNjdΓΓ (136) 
 
{f} = {fi} = ∫ NiGdΩΩ − ∫ NiqdΓΓ + ∫ NihTadΓΓ (137) 
 
Os termos [𝐊] e {f} são os mesmos determinados pelas equações (68) e (91), já para a 
matriz [C] define-se o nome como matriz de capacitância. 
Para realizar a discretização do tempo, será empregado um método que utiliza a 
aproximação do domínio do tempo por diferenças finitas para gerar uma solução numérica. 
Utilizando a série de Taylor, pode-se dizer que a temperatura em um intervalo n+1 tem o 
seguinte comportamento: 
 
Tn+1 = Tn + Δt
∂Tn
∂t
+
∆t2
2
∂2Tn
∂t2
+⋯ (138) 
58 
 
 
Se apenas os dois primeiros termos da série forem utilizados, tem-se: 
 
𝜕𝑇𝑛
𝜕𝑡
 = 
𝑇𝑛−1 − 𝑇𝑛
Δ𝑡
 (139) 
 
Introduzindo o parâmetro θ, a equação (138) passe a ser: 
 
Tn+θ = θTn+1 + (1 − θ)Tn (140) 
 
Sendo assim a equação matricial (134) se torna: 
 
[𝐂𝐭] {
Tn−1 − Tn
Δt
} + [Kt]{T}n+θ = {f}n+θ (141) 
 
Ou de uma forma expandida: 
 
([C] + θ Δt[K]){T}n+1 = ([C] − (1 − θ)Δt[K]){T}n + Δt (θ{f}n+1 + (1 − θ){f}n)
 (142) 
Para determinar os valores de θ, várias escolhas são possíveis: 
θ = 0, Método de diferenças progressivas, que é uma técnica completamente explícita. 
θ = 1, Método de diferenças Regressivas, que é uma técnica completamente implícita. 
θ = 0,5 , Método de diferenças Centrais, que é uma técnica semi implícita 
θ = 2/3, Método de Galerkin 
Considerando o valor de θ = 1, a equação (142) assume a seguinte forma: 
 
([C] + Δt[𝐊]){T}n+1 = ([C]){T}n + Δt (θ{𝐟}n+1) (143) 
 
 
Para resolver o sistema de matrizes da equação (143), define-se as matrizes [A] e [P], 
que são combinações das matrizes [C] e [K], conforme as equações (144) e (145): 
 
[A] = [C] + θ Δt[K] (144) 
[P] = [C] − (1 − θ) Δt[K] (145) 
 
59 
 
7 Disponível em: http://ta.twi.tudelft.nl/nw/users/matthias/teaching/athens/timestepping_2015.pdf acessado em 
agosto de 2016 
 
Figura 19: Representação Gráfica dos métodos de integração numérica. 
 
 
Fonte: Delft Institute of Applied Mathematics 7 
 
Substituindo [A] e [P], na equação (143): 
 
[A] {T}n+1 = [P] {T}n + {f}∗ (146) 
 
Para descobrir a temperatura no ponto {T}n+1 , aplica-se a inversa da matriz [A] 
 
{T}n+1 = ([P] {T}n + {f} ∗)[A] −1 (147) 
 
Diversos obstáculos podem surgir durante a resolução numérica dos modelos, sendo a 
falha nos valores calculados, para satisfazer o problema fisicamente e oscilações numéricas da 
solução. Oscilação numérica é o nome dado para o fenômeno em que o valor calculado oscila 
próximo a solução correta do problema, visto que para um primeiro passo do problema ele 
pode ficar abaixo, mas no seguinte ficará acima do valor correto. De acordo com o Segerlind 
(1984) a aplicação do método de diferenças regressivas é incondicionalmente estável, sem 
oscilação dos valores corretos. 
Sendo assim, para a resolução do método de Galerkin temos o seguinte sistema de 
equações, proposto por Favennec (2003) e por Bay (2003): 
 
 
{Re} = [Ce] {
∂E
∂t
(t)} + [Ke]{E(t)} − {Be} = {0} (148) 
{Rt} = [Ct] {
∂T
∂t
(t)} + [Kt]{T(t)} − {Bt} = {0} (149) 
 
61 
 
 
3.3 MODELAGEM MAGNÉTICA 
 
Iniciando pela equação (148), é possível demonstrar a rotina para geração de cada 
elemento. A primeira matriz a ser solucionada é a [Ke], que é composta das seguintes 
equações: 
 
[Ke] = ∫
1
μ
[B]T[B]dΩ +
Ω
∫
1
μr2
[N]T[N]dΩ
Ω
+ ∫
1
μr
∂
∂r
[N]T[N]dΩ
Ω
 (150) 
 
Para se calcular o primeiro termo de [K], executa-se o produto das equações (123) (124) 
e multiplica-se pela integral da equação (128), resultando na equação (151): 
 
∫
1
μ
[B]T[B]dΩ
Ω
= 
2πr̅
4A
1
μ
[
bi
2 + ci
2 bibj + cicj bibk + cick
bibj + cicj bj
2 + cj
2 bjbk + cjck
bibk + cick bjbk + cjck bk
2 + ck
2
] (151) 
 
As etapas para o cálculo do segundo termo de [K], inicia-se com o produto das matrizes 
demonstradas nas equações (126) e (127): 
 
[N]T[N] = [
Ni
2 NiNj NiNk
NiNj Nj
2 NjNk
NiNk NjNk Nk
2
] (152) 
 
A matriz da equação (152) será integrada, por meio da equação (42): 
 
2π
μr̅2
∫ [
Ni
2 NiNj NiNk
NiNj Nj
2 NjNk
NiNk NjNk Nk
2
] (Niri + Njrj +Nkrk)dAΩ (153) 
 
 
2π
μr2
∫ [
(Niri + Njrj + Nkrk)Ni
2 (Niri + Njrj + Nkrk)NiNj (Niri + Njrj + Nkrk)NiNk
(Niri + Njrj + Nkrk)NiNj (Niri + Njrj + Nkrk)Nj
2 (Niri + Njrj + Nkrk)NjNk
(Niri + Njrj + Nkrk)NiNk (Niri + Njrj + Nkrk)NjNk (Niri + Njrj + Nkrk)Nk
2
]
1
0
 dA
 (154) 
 
61