Industrial Steel Heat Treating: Numerical Simulation of Induction Heating and Aquaquenching Cooling with Mechanical Effects

mathematics Article Industrial Steel Heat Treating: Numerical Simulation of Induction Heating and Aquaquenching Cooling with Mechanical Effects José Manuel Díaz Moreno 1,† , Concepción García Vázquez 1,† , María Teresa González Montesinos 2,† , Francisco Ortegón Gallego 1,*,† and Giuseppe Viglialoro 3,†   Citation: Díaz Moreno, J.M.; García Vázquez, C.; González Montesinos, M.T.; Ortegón Gallego, F.; Viglialoro, G. Industrial Steel Heat Treating: Numerical Simulation of Induction Heating and Aquaquenching Cooling with Mechanical Effects. Mathematics 2021,9, 1203. https://doi.org/ 10.3390/math9111203 Academic Editors: Pablo Gregori and Vicente Martínez Received: 28 April 2021 Accepted: 21 May 2021 Published: 26 May 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Departamento de Matemáticas, Facultad de Ciencias, Universidad de Cádiz, 11510 Puerto Real, Spain; [email protected] (J.M.D.M.); [email protected] (C.G.V.) 2Departamento de Matemática Aplicada I, Escuela Técnica Superior de Ingeniería Informática, Avda. Reina Mercedes, s/n, 41012 Sevilla, Spain; [email protected] 3Dipartimento di Matematica e Informatica, Università di Cagliari, Via Ospedale, 72, 09124 Cagliari, Italy; [email protected] *Correspondence: [email protected] † These authors contributed equally to this work. Abstract: This paper summarizes a mathematical model for the industrial heating and cooling processes of a steel workpiece corresponding to the steering rack of an automobile. The general purpose of the heat treatment process is to create the necessary hardness on critical parts of the workpiece. Hardening consists of heating the workpiece up to a threshold temperature followed by a rapid cooling such as aquaquenching. The high hardness is due to the steel phase transformation accompanying the rapid cooling resulting in non-equilibrium phases, one of which is the hard microconstituent of steel, namely martensite. The mathematical model describes both processes, heating and cooling. During the first one, heat is produced by Joule’s effect from a very high alternating current passing through the rack. This situation is governed by a set of coupled PDEs/ODEs involving the electric potential, the magnetic vector potential, the temperature, the austenite transformation, the stresses and the displacement field. Once the workpiece has reached the desired temperature, the current is switched off an the cooling stage starts by aquaquenching. In this case, the governing equations involve the temperature, the austenite and martensite phase fractions, the stresses and the displacement field. This mathematical model has been solved by the FEM and 2D numerical simulations are discussed along the paper. Keywords: steel hardening; thermomechanical problem; phase transitions; nonlinear coupled system of PDEs/ODEs; finite element method 1. Introduction and Description of the Industrial Procedure This work deals with the thermomechanical modeling and the numerical simulation of metallurgical phase transitions of steel during an industrial heat treating. Iron and steel production has a very long history that probably dates back to some millennia. For instance, Asian swords are legendary. They are straight and had double edged blades. The popularity of these swords can be explained by their extreme hardness, obtained through a traditional heating-cooling process. Hardness implies an increasing of brittleness, so that not both blade edges are hardened but only one of them, leaving the other one resistant to breakage. Presently, steel represents a very important source of innovation. In fact, its toughness, good versatility and excellent processing properties make it one of the most commonly used in Engineering or as a building material. Steel is an iron-based alloy containing small amounts of carbon and possibly some other elements. In the automotive industry, carbon content in steel lies between 0.2 wt% Mathematics 2021,9, 1203. https://doi.org/10.3390/math9111203 https://www.mdpi.com/journal/mathematics Mathematics 2021,9, 1203 2 of 17 and 2.1 wt% and, together with other elements, represents the fundamental hardening agent, very useful to prevent dislocation motion and propagation, therefore increasing mechanical strength. It is possible to control the hardness, ductility and tensile strength of steel by changing the amount of carbon and other alloying elements, such a chromium, manganese, molybdenum, or vanadium. Hardening consists of heating up the specified region till the austenitization is reached, followed by a very rapid cooling via aquaquenching, oil spraying or salt baths. At room temperature, steel is composed by different solid phases: ferrite, pearlite, bainite, and martensite. In steels with less than 0.8 wt% carbon, known as hypoeutectoid steel, it is possible to observe a significant variation of its internal structure in a temperature range between 727 ◦ C and 912 ◦ C, approximately. In this situation, all phases are transformed into austenite (or γ− iron). Figure 1summarizes an iron-carbide phase diagram for carbon steel; it shows the condition under which the phases are stable. Temperature ◦C 400 600 800 1000 1200 1400 1600 01 2 345 6 7 Fe % C Fe3C Hypoeutectoid Hypereutectoid α (Fe) Ferrite + Cementite 727 ◦C 1148 ◦C 912 ◦C 1394 ◦C 0.77% (Eutectoid) 2.11% 4.30% Austenite + Cementite Liquid γ+ Liquid γ Austenite Figure 1. Iron-carbide phase diagram. When the cooling rate is high enough, the austenite is transformed into martensite, which is the hardest microconstituent of steel. On the other hand, if the cooling rate is not big enough, the austenite is transformed back into the other phases, i.e., ferrite, pearlite, and bainite (Figure 2) possibly along with a different partition. There are some common situations in which the industrial requirement is such that the outer surface should be hard whereas the inner core should be kept ductile. In fact, this ensures the workpiece wear resistance while reducing at the same time the material fatigue. Mathematics 2021,9, 1203 3 of 17 FERRITE PEARLITE BAINITE MARTENSITE Austenite FERRITE PEARLITE BAINITE MARTENSITE Heating Cooling Figure 2. Phase transformation of steel. There exist several industrial hardening procedure of steel: induction-conduction, flame, electron beam, laser surface, and so on (see [ 1 – 3 ]). For our purposes, we are concerned with an induction-conduction process. A conductor (coil) carries a high frequency alternating current producing a changing magnetic field and in its turn inducing eddy currents. As a result, induction heating takes place: the heat so generated (Joule’s heating) affects only the outer surface of steel components due to the skin effect. The current passing the workpiece is switched on usually for a few seconds only. When austenitization is reached, the electric current is then switched off and the workpiece is subsequently quenched in order to obtain the necessary cooling down temperature rate. In this work we are interested in the mathematical description and the numerical simulation of the hardening procedure of a car steering rack (see Figure 3a) including mechanical effects. This mathematical description corresponds to the actual industrial procedure happening in situ at the manufacturer factory. The rack bar of an automobile is a cylindrical shape workpiece made of steel, provided with a toothed region designed to engage with a pinion: the rack bar is mounted in a pinion housing in order to be displaceable along the rack longitudinal axis. The rack bar is approximately 50 cm long and 22–25 mm in diameter. These two components, rack-pinion, constitute a very important part of the steering system of a vehicle: the manufacturer company must guarantee its correct performance for at least 20 years. The toothed region of the rack bar is permanently in close contact with the pinion, and along this contact surface strong stresses will appear. Thus, in order to avoid the mechanical deformations, hindering wear, abrasions and fractures, the toothed region has to be hardened but without changing the ductility within the rest of the workpiece. ( a ) Rack bars of an automobile prior to a heat treatment showing the toothed region. Inductor Support Toothed region ( b ) The inductor and the workpiece are put in close contact during the heating process. Figure 3. Rack bars and a steel rack induction heating machine. In this paper, we will focus on the mathematical formulation of the phase transitions and the thermo-mechanical behavior of the workpiece. As long as the mechanical model Mathematics 2021,9, 1203 4 of 17 is concerned, the main quantities of interest are (Cauchy) stress and strain tensors during the treatment and the displacement field. We assume that the strain may be additively decomposed into an elastic, a thermal, and a plastic part and that the thermomechanical problem is governed by the quasistatic momentum balance and the balance law of internal energy [4,5]. In the next section we describe the mathematical modeling of both the heating and the cooling industrial processes for the hardening of a rack bar. The governing equations is a nonlinear strongly coupled system of PDEs/ODEs. Figure 4explains the way this whole system is strongly coupled concerning the electromagnetic quantities, the phase fractions, stresses, displacement field and temperature as well. Since a high alternating electric current (83 KHz in this setting) is applied, potential Maxwell’s equations in the harmonic regime are used to describe the electromagnetic induction. The main production term responsible for the increase in temperature along the critical part of the workpiece (the toothed region) is Joule’s heating coming from the eddy currents induced by the applied high alternating electric current. Additionally, we focus our attention on two solid steel phase fractions, namely austenite during the heating and cooling stages, and martensite during the cooling stage. The dynamics of these two phase fractions are described by a set of two ODEs. High temperature changes and phase transformations are accompanied by deformations. In order to deal with these mechanical effects, a quasistatic viscoelasticity model is considered. Here, the main quantities of interest are the Cauchy stress tensor, the strain tensor and the displacement field. Finally, a balance law of internal energy is taken into account. The source terms include Joule’s heating, mechanical dissipation and latent heat. Temperature θElectromagnetic field Stress σ Displacement u Phase volume fractions z= (z0,z1,z2) z0+z1+z2=1 Temperature dependent parameters Joule’s effect Transformation kinetics Transformation strain and plasticity Phase fraction dependent parameters Thermal strain Temperature dependent parameters Mechanical dissipation Magnetic dependent parameters Latent heat Phase fraction dependent parameters Transformation kinetics Figure 4. The equations governing these phenomena are strongly coupled. In Section 3we present the numerical simulations of the model described in the previous section and comment the numerical results. All these numerical simulations have been carried out through the software package FreeFem++ [ 6 ]. This is a free programming language based on the finite element method. FreeFem++ is implemented in terms of the variational forms of the corresponding problems, so that it is straightforward to solve problems involving PDEs (both 2D or 3D) from several branches of physics and engineering. Though the numerical simulations presented in this paper are in 2D, they may provide with a valuable information to the manufacturer company so that it may be useful to optimize this industrial process. 2. Mathematical Modeling Let Ω⊂RNN= 2, 3, be an open set with a Lipschitz-continuous boundary (see Figure 5). Ω=Ωc∪Ωs∪Γed . Ωc and Ωs stand for the inductor (copper-made) and the steel rack bar, respectively. Γed =¯ Ωc∩¯ Ωs is the contact surface between the two conductors. We consider two opposite faces Γ1 , Γ2⊂∂Ωc as shown in Figure 5where the potential difference is applied. We also put Γ=Γ1∪Γ2. Mathematics 2021,9, 1203 5 of 17 Mathematics 2021,1, 0 5 of 17 Ωs(steel) Ωc(copper) Ωc Γ1Γ2 Figure 5. Rack bar (steel) and inductor (copper). 3D CAD model built in FEMAP from Autocad file supplied by the automotive industry. To deal with the heating-cooling stages, a time partition [ 0, Tf] = [ 0, Th)∪[Th , Tf] , Tf>Th> 0 is considered. We distinguish our model in each time interval so that we may identify the main difficulties and phenomena involved. In both cases, the coupling between the thermomechanical phenomenon and the phase transformations will be analyzed ([5,7–11]). The electromagnetic unknowns are ϕ , standing for the scalar electric potential, and A representing the complex magnetic field. The unknowns θ(x , t) and u(x , t) stand for the temperature and the displacement field measured from the original geometry at the point x∈Ωs and at the instant t∈[ 0, Tf] . By a(x , t) (respectively, m(x , t) ) we will denote the austenite (respectively, martensite) phase fraction. We consider an initial distribution of all steel phases given by z0 . This means that z0(x , t) + a(x , t) + m(x , t) = 1 for all (x,t)∈¯ Ωs×[0, Tf]. Notice that, in this setting, we have z0(x, 0) = 1 in ¯ Ωs. 2.1. Heating Stage: Governing Equations and Parameters Determination The inductor and the rack bar are put in close contact as it is shown in Figure 5. In this setting, both conductors constitute the coil. During the time interval ( 0, Th) we let pass a high alternating current through the coil (83 KHz). The following equations will be considered ([5,7,12–14]): Electromagnetics                        ∇ · (b(θ)∇ϕ) = 0, in ΩTh=Ω×(0, Th), ∂ϕ ∂n=0, on (∂Ω\Γ)×(0, Th), ϕ=0 on Γ1×(0, Th),ϕ=ϕ0on Γ2×(0, Th), b0(θ)iωA+∇ × 1 µ∇ × A−δ∇(∇ · A) = −b0(θ)∇ϕin D×(0, Th), A=0, on ∂D×(0, Th). (1) Phase fractions    ∂a ∂t=fa(θ,a), in Ωs×(0, Th) a(0) = 0, in Ωs. (2) Figure 5. Rack bar (steel) and inductor (copper). 3D CAD model built in FEMAP from Autocad file supplied by the automotive industry. To deal with the heating-cooling stages, a time partition [ 0, Tf] = [ 0, Th)∪[Th , Tf] , Tf>Th> 0 is considered. We distinguish our model in each time interval so that we may identify the main difficulties and phenomena involved. In both cases, the coupling between the thermomechanical phenomenon and the phase transformations will be analyzed [5,7–11]. The electromagnetic unknowns are ϕ , standing for the scalar electric potential, and A representing the complex magnetic field. The unknowns θ(x , t) and u(x , t) stand for the temperature and the displacement field measured from the original geometry at the point x∈Ωs and at the instant t∈[ 0, Tf] . By a(x , t) (respectively, m(x , t) ) we will denote the austenite (respectively, martensite) phase fraction. We consider an initial distribution of all steel phases given by z0 . This means that z0(x , t) + a(x , t) + m(x , t) = 1 for all (x,t)∈¯ Ωs×[0, Tf]. Notice that, in this setting, we have z0(x, 0) = 1 in ¯ Ωs. 2.1. Heating Stage: Governing Equations and Parameters Determination The inductor and the rack bar are put in close contact as it is shown in Figure 5. In this setting, both conductors constitute the coil. During the time interval ( 0, Th) we let pass a high alternating current through the coil (83 KHz). The following equations will be considered [5,7,12–14]: Electromagnetics                        ∇ · (b(θ)∇ϕ) = 0, in ΩTh=Ω×(0, Th), ∂ϕ ∂n=0, on (∂Ω\Γ)×(0, Th), ϕ=0 on Γ1×(0, Th),ϕ=ϕ0on Γ2×(0, Th), b0(θ)iωA+∇ × 1 µ∇ × A−δ∇(∇ · A) = −b0(θ)∇ϕin D×(0, Th), A=0, on ∂D×(0, Th). (1) Phase fractions    ∂a ∂t=fa(θ,a), in Ωs×(0, Th) a(0) = 0, in Ωs. (2) Mathematics 2021,9, 1203 6 of 17 Deformations        ∇ · σ(u) = 0, in Ωs×(0, Th), u=0, on Γ0 h×(0, Th), σ(u)·n=0, on Γ1 h×(0, Th), (3) with σ(u) = K:(ε(u)−q(θ,a)I). (4) Temperature                    α(θ,a,σ)∂θ ∂t−∇·(k(θ)∇θ) = 1 2b(θ)|iωA+∇ϕ|2 +β(θ,a,σ)∂a ∂t, in Ωs×(0, Th), θ(0) = θ0, in Ωs, ∂θ ∂n=0, on ∂Ωs×(0, Th). (5) The system (1), describing the electromagnetic dynamics, comes from Maxwell’s equations in the harmonic regime. Here b(θ) is the electric conductivity, i is the imaginary unit, ω= 2 πf is the pulsation, f being the frequency of the alternating current, µ=µ(x) is the magnetic permeability, δ> 0 is a small constant and b0(x , s) = b(x , s) if x∈Ω , b0(x , s) = 0 elsewhere;. Notice that D⊂R2 is a big enough and smooth domain containing the set of conductors. The magnetic permeability is then defined by the following piecewise constant function µ(x) =    µ0if x∈D\(Ωc∪Ωs), µ1if x∈Ωs, µ2if x∈Ωc, where the µ0 , µ1 and µ2 stand for the vacuum permeability, steel permeability and copper permeability, respectively. In Equation (2), a Leblond-Devaux model [ 15 ] is considered to simulate the austenite transformation once the threshold temperature As is attained, particularly, the rate of change fais modeled as fa(θ,a) = max1 τa(θ)(aeq(θ)−a), 0H(θ−As). Herein, H stands for the Heaviside function H(s) = 1 if s⩾ 0 and is zero elsewhere. The functions aeq(θ),τa(θ)(see Figure 6) can be determined from experimental data. As far as the quasistatic model of viscoelasticity (3) is concerned, Γ0 h⊂∂Ωs is the surface boundary where ∂Ωs is clamped (null displacement) during the heating stage and Γ1 h=∂Ωs\Γ0 h . As it was said above, u is the displacement field and σ denotes the stress tensor, related to strain ε(u) = 1 2(∇u+∇ut) through the constitutive law (4). In this last expression (see [16]) •Kstands for the isotropic stiffness tensor; •q(θ , a) = q0( 1 −a)(θ−θref 0) + qaa(θ−θref a) is a thermal contribution to strain, deduced from available dilatometry experiments (see Figure 7). As long as the temperature is increasing, the material expansion is modeled using a linear fitting. During the transformation into the austenite phase, a contraction range can be observed before going on with the usual expansion, but now with a higher linear coefficient. Mathematics 2021,9, 1203 7 of 17 AsAfθ τa(θ) aeq(θ) 1 0 Figure 6. Characterization of aeq (red line) and τa(blue line) as temperature-dependent functions. Temperature Thermal Expansion θref 0AsAf MfMs Figure 7. Schematic dilatometry experiments. Finally, arguing as in [ 5 ], a balance law is developed for the temperature. Thus, in (5), k(θ)stands for thermal conductivity, αand βare given, respectively, by α(θ,a,σ) = ρc+tr σ∂q ∂θ ,β(θ,a,σ) = ρL+tr σ∂q ∂a, ρ being the steel density, c the specific heat at constant strain and L the latent heat. Additionally, θ0 is the room temperature, which is taken to be θ0= 300 K in all our numerical experiments. The main production term in (5), responsible for the increasing in temperature is the quadratic term 1 2b(θ)|iωA+∇ϕ|2and is called Joule’s heating. We put θTh=θ(Th) and aTh=a(Th) the temperature distribution and the austenite fraction, respectively, at the end of the heating stage. Both functions θTh and aTh will be used as the initial temperature and austenite fraction, respectively, in the cooling stage. Mathematics 2021,9, 1203 8 of 17 2.2. Cooling Stage: Governing Equations and Parameters Determination At t=Th the current is switched off and for Th⩽t⩽Tf , the aquaquenching takes place. The workpiece is sprayed down with oily water at a constant temperature θq . The temperature goes down so quickly in the austenized toothed region that the martensite phase fraction is forced to appear. However, these phase transformations induce a plastic deformation. According to [ 17 , 18 ], this TRIP (TRansformation-Induced-Plasticity) contribution to strain for the austenite-martensite transformation is taken into account in our model. Phase fractions                  ∂a ∂t=fa(θ,a)−∂m ∂t, in Ωs×(Th,Tf), ∂m ∂t=fma,m,θ,∂θ ∂t, in Ωs×(Th,Tf), a(Th) = aTh, in Ωs, m(Th) = 0, in Ωs. (6) Deformations        ∇ · σ(u) = 0, in Ωs×(Th,Tf), u=0, on Γ0 c×(Th,Tf), σ(u)·n=0, on Γ1 c×(Th,Tf), (7) σ(u) = K:ε(u)−q(θ,a,m)I−Zt 0γθ,m,∂m ∂tSdτ. (8) Temperature                      α(θ,a,m,σ)∂θ ∂t− ∇ · (k(θ)∇θ) = βa(θ,a,m,σ)∂a ∂t +βm(θ,a,m,σ)∂m ∂t+γθ,m,∂m ∂t|S|2, in Ωs×(Th,Tf) θ(Th) = θTh, in Ωs, −k(θ)∂θ ∂n=¯ β(θ−θq), on ∂Ωs×(Th,Tf). (9) In (6), a Koistinen-Marburger model is taken into account in order to describe the martensite transformation (see [ 19 ]). This diffusionless transformation is quite different from that related to austenite and it is highly dependent on the cooling rate: fma,m,θ,∂θ ∂t=cmaH−∂θ ∂tH(Ms−θ), where the threshold value Ms indicates the starting temperature from which the martensite appears, whereas cmis a constant determined from experimental data. In (7) and (8), the thermal strain should be conveniently modified in order to take into account the new linear growing during the cooling as is shown in Figure 7. This yields q(θ,a,m) = q0(1−a−m)θ−θref 0+qaaθ−θref a+qmmθ−θref m. Mathematics 2021,9, 1203 9 of 17 Moreover, as in [7,20] the TRIP effect is modeled by ∂εtrip ∂t=λ1(θ)∂λ2(m) ∂m ∂m ∂tSdef =γθ,m,∂m ∂tS, (10) where S=σ−1 N(tr σ)I is the deviatoric stress tensor and tr σ denotes the trace of σ ; λ1(θ) = 3 2K , and K is experimentally determined; λ2 is a suitably selected increasing function such that λ2( 0 ) = 0 and λ2( 1 ) = 1, for instance λ2(s) = s( 2 −s) as in [ 17 ], or λ2(s) = s(1−log(s)) as in [18]. Finally, in (9) the functions βaand βmare given by βa(θ,a,m,σ) = −ρL+tr σ∂q ∂aand βm(θ,a,m,σ) = ρL+tr σ∂q ∂m. Notice that the aquaquenching is modeled through the Robin boundary condition in (9) where ¯ βstands for a heat transfer coefficient. 3. Numerical Simulations Using the FreeFem++ software [ 6 ] we carried out some numerical simulations of the model (1)–(9). We used the finite element method for the space approximation and a Crank-Nicolson scheme for the time discretization. Actually, a P2 -Lagrange approximation for ϕ , A , θ and u , and P1 for a and m . Ordinary differential equations are discretized using Heun’s method. A fixed time step has been used from t= 0 to t=Th= 5.5, namely ∆t= 5.5 / 110, whereas a variable time step has been considered from t=Th to t=Tf= 11.0 in order to deal with the drastic decreasing in temperature during the first iterations of the cooling stage. To do so, we first consider the function s(t)=(t− 5.5 )2/ 5.5 and put r0=t0= 5.5, then define recursively for n=1, . . . , 110, rn=rn−1+∆tand tn=tn−1+s(rn)−s(rn−1). In this way, the first time iterations of the cooling stage refer to tighter time steps than the ones of the final iterations. In our simulation, the workpiece was heated up during Th= 5.5 s from a room temperature of θ0= 300 K followed by the aquaquenching for another 5.5 s ( Tf= 11). The temperature of the oil used in the aquaquenching is θq=315 K. In Figure 8, some pictures are displayed for the meshes used in our 2D numerical simulations. It is important to have a high density of triangles close to the toothed region. The material parameters have been taken from data of a specific steel with a given composition (Table 1) assuring the desirable mechanical properties after the heat treating [ 21 ]. In Table 2 we can find the actual values we have taken in our numerical experiments. The parameters in the thermal stresses modeling has been taken from [22,23]. Additionally, the value of the small parameter δ appearing in the equation for the complex magnetic field is taken to be δ= 10 −15 . In fact, the term −δ∇(∇ · A) appearing in (1) is artificial. It is used in order to stabilize the numerical resolution of (1). Of course, other approaches are possible (see [5,7]). Table 1. Chemical composition of 42CrMo4 steel. Elements C Si Mn P S Cr Cu Mo Ni Content %wt 0.420 0.260 0.680 0.019 0.021 1.020 0.330 0.180 0.150 Mathematics 2021,9, 1203 16 of 17 austenite and martensite, coupled with a strongly coupled and nonlinear system of PDEs for the electric potential, the complex magnetic field, stresses and displacement field and the temperature. Many industrial heating procedures are based on an AC induction technique so that the workpiece is put near the coil but with neat separation between them. However, this is not the case described in this paper, mainly due to the special geometry of the considered workpiece. In this setting, the inductor machine and the workpiece are put in close contact so that all together constitute a coil. The eddy currents combined with the skin effect result in Joule’s heating along the toothed region of the rack bar which leads to the high heating of the workpiece critical part. Once the mathematical model has been fully described, including all the relevant physical quantities, we carried out some numerical simulations based on the FEM. To do so, we have considered a mesh exhibiting an element high density around the rack toothed region. Though we only performed 2D numerical simulations, all these numerical results may provide a valuable information to the manufacturer and can be useful in order to optimize this industrial process. The analysis of the numerical simulations shows how the toothed region along the rack bar is heated up until austenitization is reached. During this stage, one can observe the mechanical effects on the workpiece: expansions due to the raise of temperature and at the same time volume reduction on the austenitized region. On the other hand, during the aquaquenching, the toothed region is radically cooled down so that martensite is produced along the toothed region and again some mechanical effects are observed: the workpiece undergoes expansions where martensite occurs and contractions elsewhere. Moreover, the rate-dependent inelastic behavior of the workpiece is justified by the plastic term appearing in the constitutive laws. The austenite transformation induces the teeth shrinking during heating; successively, the same workpiece expands again during cooling according to the martensite transformation, but without reaching its initial configuration. Obviously, these deformations are accompanied by stresses. A computation of the von Mises stresses has shown the regions where the workpiece is particularly stressed. The next step in this research should be to implement the resolution of this model in a 3D setting. Author Contributions: All authors contributed equally, in one way or another, in the development of this work, including different aspects on conceptualization, methodology, software, validation, formal analysis, investigation, resources, data curation, writing—original draft preparation, writing— review and editing, visualization. Supervision, project administration and funding acquisition were managed by F.O.G. and G.V. All authors have read and agreed to the published version of the manuscript. Funding: This research was partially supported by Ministerio de Educación y Ciencia under grants MTM2010-16401 and TEC2017-86347-C2-1-R with the participation of FEDER, and Consejería de Educación y Ciencia de la Junta de Andalucía, research group FQM–315. Giuseppe Viglialoro is a member of the Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM) and he is partially supported by the research projects Evolutive and stationary Partial Differential Equations with a focus on biomathematics, funded by Fondazione di Sardegna (2019), and by MIUR (Italian Ministry of Education, University and Research) Prin 2017 Nonlinear Differential Problems via Variational, Topological and Set-valued Methods (Grant Number: 2017AYM8XW). Conflicts of Interest: The authors declare no conflict of interest. References 1. ASM Handbook Committee. ASM Handbook, Heat Treating; ASM International: Materials Park, OH, USA, 1991; Volume 4. 2. Digges, T.G.; Rosenberg, S.J.; Geil, G.W. Heat Treatment and Properties of Iron and Steel; National Bureau of Standards: Washington, DC, USA, 1966. Available online: https://digital.library.unt.edu/ark:/67531/metadc13269/ (accessed on 21 May 2021). Mathematics 2021,9, 1203 17 of 17 3. Rudnev, V.; Loveless, D.; Cook, R.; Black, M. Induction Hardening of Gears: A Review. Heat Treat. Met. 2003,4, 97–103. 4. Arenas, M.J.; Hömberg, D.; Lasarzik, R.; Mikkonen, P.; Petzold, T. Modelling and simulation of flame cutting for steel plates with solid phases and melting. J. Math. Ind. 2020,10, 18. [CrossRef] 5. Hömberg, D. A mathematical model for induction hardening including mechanical effects. Nonlinear Anal. RWA 2004 ,5, 55–90. [CrossRef] 6. Hecht, F. New development in FreeFem++. J. Numer. Math. 2012,20, 251–265. [CrossRef] 7. Chełminski, K.; Hömberg, D.; Kern, D. On a thermomechanical model of phase transitions in steel. Adv. Math. Sci. Appl. 2008 ,18, 119–140. 8. Fuhrmann, J.; Hömberg, D.; Uhle, M. Numerical simulation of induction hardening of steel. Compel 1999 ,18, 482–493. [CrossRef] 9. Hömberg, D.; Weiss, W. PID control of laser surface hardening of steel. IEEE Trans. Control Syst. Technol. 2006 ,14, 896–904. [CrossRef] 10. Hömberg, D.; Petzold, T.; Rocca, E. Analysis and simulations of multifrequency induction hardening. Nonlinear Anal. RWA 2015 , 22, 84–97. [CrossRef] 11. Hömberg, D.; Liu, Q.; Montalvo-Urquizo, J.; Nadolski, D.; Petzold, T.; Schmidt, A.; Schulz, A. Simulation of multi-frequency- induction-hardening including phase transitions and mechanical effects. Finite Elem. Anal. Des. 2016,121, 86–100. [CrossRef] 12. Moreno, J.M.D.; Vázquez, C.G.; Montesinos, M.T.G.; Gallego, F.O. Analysis and numerical simulation of an induction-conduction model arising in steel heat treating. J. Comput. Appl. Math. 2012,236, 3007–3015. [CrossRef] 13. Moreno, J.M.D.; Vázquez, C.G.; González, M.T.; Ortegón, F.; Viglialoro, G. Heat treatment of a steel rack: Modeling and numerical simulation. In Proceedings of the 2012 International Conference on Computational and Mathematical Methods in Science and Engineering, Murcia, Spain, 2–6 July 2012; Vigo-Aguiar, J. Ed.; Volume II, pp. 480–486. 14. Moreno, J.M.D.; Vázquez, C.G.; González, M.T.; Ortegón, F.; Viglialoro, G. Mathematical modeling of heat treatment for a steering rack including some mechanical effects. J. Numer. Math. 2012,20, 215–232. [CrossRef] 15. Leblond, J.B.; Devaux, J. A new kinetic model for anisothermal metallurgical transformations in steels including effects of austenite grain size. Acta Metall. 1984,32, 137–146. [CrossRef] 16. Rohde, J.; Jeppsson, A. Literature review of heat treatment simulations with respect to phase transformation, residual stresses and distortion. Scand. J. Metall. 2000,29, 47–62. [CrossRef] 17. Denis, S.; Gautier, E.; Stöström, S.; Simon, A. Influence of stresses on the kinetics of pearlitic transformation during continuous cooling. Acta Metall. 1987,35, 1621–1632. [CrossRef] 18. Leblond, J.B.; Devaux, J.; Devaux, J.C. Mathematical modelling of transformation plasticity in steels I: Case of ideal-plastic phases. Int. J. Plast. 1989,5, 551–572. [CrossRef] 19. Koistinen, D.P.; Marburger, R.E. A general equation prescribing the extent of the austenite-martensite transformation in pure iron-carbon alloys and plain carbon steels. Acta Metall. 1959,7, 59–60. [CrossRef] 20. Fisher, F.D.; Reisner, G.; Werner, E.; Tanaka, K.; Cailletaud, G.; Antretter, T. A new view on transformation induced plasticity (TRIP). Int. J. Plast. 2000,16, 723–748. [CrossRef] 21. Zhong, H.; Wang, Z.; Gan, J.; Wang, X.; Yang, Y.; He, J.; Wei, T.; Qin, X. Numerical simulation of martensitic transformation plasticity of 42CrMo steel based on spot continual induction hardening model. Surf. Coat. Technol. 2020 ,385, 125428. [CrossRef] 22. Lusk, M.T.; Wang, W.; Sun, X.; Lee, Y.K. On the role of kinematics in constructing predictive models of austenite decomposition. In A Symposium on the Thermodynamics, kinetics, Characterizaion and Modeling of Austenite Formation and Decomposition, Proceedings of the Materials Science and Technology 2003 Meeting, Chicago, IL, USA, 9–12 November 2003; Damm, E.B., Merwin, M.J., Eds.; Yonsei University: Seoul, Korea, 2003; pp. 311–331. 23. Rajeev, P.T.; Farris, T.N.; Chandrasekar, S.; Andreski, B.; Brada, G. Heat tratment of steels: Prediction of the microstructure, residual distortion and residual stress. In Heat Treating: Proceedings of the 18th Conference, Rosemont, IL, USA, 12–15 October 1998; Wallis, R.A., Walton, H.W., Eds.; ASM International: Materials Park, OH, USA, 1999; pp. 467–475. 24. Díaz, J.M.; García, C.; González, M.T.; Ortegón, F. On an electrothermomechanical model arising in steel heat treating. In Numerical Simulation in Physics and Engineering, Proceedings of the XIV Spanish-French Jacques-Louis Lions School, A Coruña, Spain, 6–10 September, 2010; Arregui, I., García, J.A., Vázquez, C., Eds.; Universidade da Coruña: A Coruña, Spain, 2010; pp. 265–281.