Steel heat treating: mathematical modelling and numerical simulation of a problem arising in the automotive industry

Steel heat treating: mathematical modelling and numerical simulation of a problem arising in the automotive industry José Manuel Dı́az Moreno1, Concepción Garcı́a Vázquez1, Marı́a Teresa González Montesinos2, Francisco Ortegón Gallego*,1, Giuseppe Viglialoro3 1Departamento de Matemáticas, Facultad de Ciencias, Universidad de Cádiz, 11510 Puerto Real, SPAIN, josemanuel.diaz@uca.es, concepcion.garcia@uca.es, francisco.ortegon@uca.es.


Introduction
In the automative industry, many workpieces such gears, bearings, racks and pinions, are made of steel. Steel is an alloy of iron and carbon. Generally, industrial steel has a carbon content up to about 2 wt%. Other alloying elements may be present, such as Cr and V in tools steels, or Si, Mn, Ni and Cr in stainless steels. Most structural components in mechanical engineering are made of steel. Certain of these components, such as toothed wheels, bevel gears, pinions and so on, engaged each others in order to transmit some kind of (rotational or longitudinal) movement. In this situation, the contact surfaces of these components are particularly stressed. The goal of heat treating of steel is to attain a satisfactory hardness. Prior to heat treating, steel is a soft and ductile material. Without a hardening treatment, and due to the surface stresses, the gear teeth will soon get damaged and they will no longer engage correctly.
In this work we are interested in the mathematical description of the hardening procedure of a car steering rack (see Figure 1). This particular situation is one of the major concerns in the automotive industry. In this case, the goal is to increase the hardness of the steel along the tooth line and at the same time keeping the rest of the workpiece soft and ductile in order to reduce fatigue. This problem is governed by a nonlinear system of partial differential equations coupled with a certain system of ordinary differential equations. Once the full system is set we perform some numerical simulations. Solid steel may be present at different phases, namely austenite, martensite, bainite, pearlite and ferrite. The phase diagram of steel is shown in Figure 2. For a given wt% of carbon content up to 2.11, all steel phases are transformed into austenite provided the temperature has been raised up to a certain range. The minimum austenization temperature (727 • ) is attained for a carbon content of 0.77 wt% (eutectoid steel). Upon cooling, the austenite is transformed back into the other phases (see Figure 3), but its distribution depends strongly on the cooling strategy ( [4,12]).
Martensite is the hardest constituent in steel, but at the same time is the most brittle, whereas pearlite is the softest and more ductile phase. Martensite derives from austenite and can be obtained only if the cooling rate is high enough. Otherwise, the rest of the steel phases will appear.
The hardness of the martensite phase is due to a strong supersaturation of carbon atoms in the iron lattice and to a high density of crystal defects. From the industrial standpoint, heat treating of steel has a collateral problem: hardening is usually accompanied by distortions of the workpiece. The main reasons of these distortions are due to (1) thermal strains, since steel phases undergo different volumetric changes during the heating and cooling processes, and (2) experiments with steel workpieces under applied loading show an irreversible deformation even when the equivalent stress corresponding to the load is in the elastic range. This effect is called transformation induced plasticity.
The heating stage is accomplished by an induction-conduction procedure. This technique has been successfully used in industry since the last century. During a time interval, a high frequency current passes through a coil generating an alternating magnetic field which induces eddy currents in the workpiece, which is placed close to the coil. The eddy currents dissipate energy in the workpiece producing the necessary heating.

Mathematical modeling
We consider the setting corresponding to Figure 4. The domain Ω c represents the inductor (made of copper) whereas Ω s stands for the steel workpiece to be hardened. Here, the coil is the domain Ω = Ω s ∪ Ω c ∪ S 0 . In this way, the workpiece itself takes part of the coil.
In order to describe the heating-cooling process, we will distinguish two subintervals forming a parti- The first one [0, T h ) corresponds to the heating process. All along this time interval, a high frequency electric current is supplied through the conductor which in its turn induces a magnetic field. The combined effect of both conduction and induction gives rise to a production term in the energy balance equation (14), namely b(θ)|A t +∇φ| 2 . This is Joule's heating which is the principal term in heat production. In our model, we will only consider three steel phase fractions, namely austenite (a), martensite (m), and the rest of phases (r). In this way, we have a + m + r = 1 and 0 ≤ a, m, r ≤ 1 in Ω s × [0, T ]. At the initial time we have r(0) = 1 in Ω s . Upon heating only austenite can be obtained. In particular m = 0 in Ω s × [0, T h ] and the transformation to austenite is derived at the expense of the other phase fractions (r).
At the instant t = T h , the current is switched off and during the time interval [T h , T c ] the workpiece is severely cooled down by means of aqua-quenching.

The heating model
The current passing through the set of conductors Ω = Ω c ∪ Ω s ∪ S 0 is modeled by the electric potential difference, ϕ 0 , applied on the surface Γ 2 ⊂ Ω c (see Figure 4). Notice that the applied potential on Γ 1 is zero. In the sequel, we put Γ = Γ 1 ∪ Γ 2 .
The heating model involves the following unknowns: the electric potential, φ; the magnetic vector potential, ; the austenite phase fraction, a; and the temperature, θ. Among them, only A is defined in the domain D containing the set of conductors Ω. On the other hand, since the inductor and the workpiece are in close contact, both φ and θ are defined in Ω. Since phase transitions only occur in the workpiece, we may neglect deformations in Ω c . This implies that σ , u and a are only defined in the workpiece Ω s .
Since electromagnetic fields generated by high frequency currents are sinusoidal in time, both the electric potential, φ, and the magnetic potential field, A, take the form ( [1,2,14,15] where M is a complex-valued function or vector field, and ω = 2πf is the angular frequency, f being the electric current frequency. In general, M also depends on t, but at a time scale much greater than 1/ω. In this way, we may introduce the complex-valued fields ϕ and A as As a far as the numerical simulation of a system like (2)-(15) is concerned, the introduction of the new variables ϕ and A is quite convenient since the time scale describing the evolution of both ϕ and A is much smaller than that of the temperature θ. In the case of steel heat treating, f is about 80 KHz. The heating model reads as follows ( [3,9,10,7]):  Figure 4: Domains D, Ω = Ω s ∪Ω c ∪S 0 and the faces Γ 1 , Γ 2 ⊂ Ω c . The inductor Ω c is made of copper. The workpiece contains a toothed part to be hardened by means of the heating-cooling process described below. It is made of a hypoeutectoid steel. The domain is taken as a big enough rectangle containing both the inductor and the rack Here, b(θ) is the electrical conductivity (by b(θ) we mean the function (x, t) → b(x, θ(x, t)), and also for κ(θ), etc.); ϕ 0 represents the potential external source. The domain D containing the set of conductors is taken big enough so that the magnetic vector potential A vanishes on its boundary ∂D (in our model, is taken to be a 2D rectangle or a 3D cube). Since both σ and a are only defined in Ω s , when they appear in a term referred in Ω, we mean that this term vanishes outside Ω s (for instance, −ρL a a t appearing in (13) is the magnetic permeability; δ > 0 is a small constant; F is a given external force (usually F = 0); K = K ijkl , 1 ≤ i, j, k, l ≤ 3 is the stiffness tensor. Steel can be considered as an isotropic and homogenous material so that K ijkl =λδ ij δ kl +μ(δ ik δ jl +δ il δ jk ), for all i, j, k, l ∈ {1, 2, 3} whereλ ≥ 0 andμ > 0 are the Lamé coefficients of steel; ε(u) = 1 2 (∇u + ∇u T ) is the strain tensor; A 1 (a, m, θ)I models the thermal strain, I being the 3×3 unity matrix, whereas A 1 (a, m, θ) is defined as and in its turn q a , q m and q r are the thermal expansion coefficients of the phase fractions a, m and r, respectively, and θ a , θ m and θ r are reference temperatures (notice that during the heating stage is m = 0); t 0 γ(a, m, a t , m t , θ)S dτ gives the model, through the function γ, of the transformation induced plasticity strain tensor, where S = σ − 1 3 tr σ I is the deviator of σ , that is, the trace free part of the stress tensor; Γ 0 is a certain smooth enough part of ∂Ω s ; n is the unit outer normal vector to the referred boundary; the functions τ a (θ), a eq (θ) are given from experimental data (see Figure 5), and H is the Heaviside function; κ(θ) is the thermal conductivity; the functions appearing in (13) are given as follows α(θ, a, m, σ ) = ρc ε − 9κq(a, m) 2 θ − q(a, m) tr σ , where ρ and c ε are the steel density and the specific heat capacity at constant strain, respectively,κ = 1 3 (3λ + 2μ) is the bulk modulus, and q(a, m) is defined as Finally, L a > 0 is the latent heat related to the austenite phase fraction. Notice that, in a more general situation ρ, c ε and L a may also depend on a, m and/or θ.
a eq (θ) Figure 5: Functions a eq and τ a . Equations (2) and (5) derive from Maxwell's equations. In [9], it is assumed the Coulomb gauge condition for the magnetic vector potential, namely, ∇ · A = 0. Here, we do not impose this condition since this makes appear an undesired pressure gradient in the equation for A. In its turn, we include a penalty term in this equation of the form −δ∇(∇ · A). In doing so, both the theoretical analysis and the numerical simulations are simplified. Equation (7) is a quasistatic balance law of momentum and (8) is Hooke's law. The transformation to austenite from the initial phase r(0) = 1 is described in (11).
Finally, equation (13) derives from the balance law of internal energy. As it has been pointed out above, Joule's heating is the main responsible in heat production. Since γ(a, m, a t , m t , θ)|S| 2 ≥ 0, the contribution of the transformation induced plasticity to the energy balance is also a production term. On the other hand, during the heating stage we have a t ≥ 0 so that −ρL a a t ≤ 0. This means that the transformation to austenite absorbs energy, which is released during the cooling stage.

The cooling model
The heating process ends, the high frequency current passing through the coil is switched-off and aquaquenching begins. The quenching is just modeled via the Robin boundary condition given in (25).
We put a T h = a(T h ), that is, a T h is the austenite phase fraction distribution at the final heating instant T h obtained from (11). In the same way, we define θ T h = θ(T h ). Obviously, these functions will be taken as the initial phase fraction distribution and temperature, respectively, in the cooling model. Here we use the Koistinen-Marburger model ( [11,13]) for the description of the transformation to martensite from austenite.
The cooling model reads as follows In (22) c m > 0 is a constant value. Also, in (24), L m > 0 is the latent heat related to the martensite phase fraction. The function β(x, t) in (25) is a heat transfer coefficient and is given by where β 0 (t) > 0 (usually taken to be constant). Finally, θ e is the temperature of the quenchant.
The mathematical analysis of a system similar to (16)-(26) can be seen in [3]. In this reference, an existence result is shown assuming that the data are smooth enough and T c − T h is sufficiently small. Figure 6: Domain triangulation. The mesh contains 61790 triangles and 30946 vertices.

Dh
Using the Freefem++ package ( [8]), we have performed some numerical simulations for the approximation of the solution to the systems (2)- (15) and (16)-(26). We want to describe the hardening treatment of a car steering rack during the heating-cooling process. The goal is to produce martensite along the tooth line together with a thin layer in its neighborhood inside the steel workpiece ( [5,6]). Figure 7: Domain triangulation. Element density near three teeth. Figure 4 shows the open sets D, Ω = Ω s ∪ Ω c ∪ S and the faces Γ 1 and Γ 2 (they appear stick together in this figure) which intervene in the setting of the problem. The workpiece contains a toothed part to be hardened by means of the heating-cooling process described above. It is made of a hypoeutectoid steel. The open set D \Ω is air. The magnetic permeability µ in (5) is then given by
The martensite phase can only derive from the austenite phase. Thus we need to transform first the critical part to be hardened (the tooth line) into austenite. For our hypoeutectoid steel, austenite only exists in a temperature range close to the interval [1050, 1670] (in K). During the first stage, the workpiece is heated up by conduction and induction (Joule's heating) which renders the tooth line up to the desired temperature. In order to transform the austenite into martensite, we must cool it down at a very high rate. This second stage is accomplished by aquaquenching.
In this simulation, the final time of the heating process is T h = 5.5 seconds and the cooling process extends also for 5.5 seconds, that is T c = 11.
We have used the finite elements method for the space approximation and a Crank-Nicolson scheme for the time discretization. Figures 6 and 7 show the triangulation of D in our numerical simulations. We have used P 2 -Lagrange approximation for ϕ, A and θ and P 1 for a and m.
In Figure 8 we can see the temperature distribution of the rack along the tooth line at different instants of the the heating-cooling process. The initial temperature is θ 0 = 300K. At t = 5.5 the heating process ends and the computed temperature shows that the temperature along the rack tooth line lies in the interval [1050, 1670] (K).   Figure 9 shows the austenite evolution from the beginning of the cooling stage. Blue corresponds to 0% while red is 100%. We observe that martensite starts to appear, approximatively, one second after the beginning of the cooling stage. At the final instant, all the amount of austenite has been transformed into martensite as it is shown in Figure 10 .   Figure 11 shows the austenization along the tooth line at the end of the heating process T = 5.5 seconds. Figure 12 shows the final distribution of martensite from austenite along the rack tooth line through the cooling stage t = 11 seconds. We have good agreement versus the experimental results obtained in the industrial process.
During the heating-cooling process, the workpiece is deformed so that an industrial rectification is needed (or otherwise the rack would be useless).

Conflict of Interest
The authors declare no conflict of interest.