A semi-analytical coupled simulation approach for induction heating

Abstract

The numerical simulation of the induction heating process can be computationally expensive, especially if ferromagnetic materials are studied. There are several analytical models that describe the electromagnetic phenomena. However, these are very limited by the geometry of the coil and the workpiece. Thus, the usual method for computing more complex systems is to use the finite element method to solve the set of equations in the multiphysical system, but this easily becomes very time consuming. This paper deals with the problem of solving a coupled electromagnetic - thermal problem with higher computational efficiency. For this purpose, a semi-analytical modeling strategy is proposed, that is based on an initial finite element computation, followed by the use of analytical electromagnetic equations to solve the coupled electromagnetic-thermal problem. The usage of the simplified model is restricted to simple geometrical features such as flat or curved surfaces with great curvature to skin depth ratio. Numerical and experimental validation of the model show an average error between 0.9% and 4.1% in the prediction of the temperature evolution, reaching a greater accuracy than other analyzed commercial softwares. A 3D case of a double-row large size ball bearing is also presented, fully validating the proposed approach in terms of computational time and accuracy for complex industrial cases.

Introduction

Induction heating of metallic components is a highly integrable, fast and efficient heating process that is widely used in the industry, not only for hardening applications but also for bonding, melting and heating prior to hot working and forging [1]. The principle of induction heating is to place an electrically conductive workpiece close to an inductor, where alternating currents are flowing through. Because of electromagnetic phenomena, a magnetic field is generated around the inductor, which penetrates any material placed close to it. Because of the penetration of the alternating magnetic field, induced currents, usually called eddy currents, appear on the surface of the component, and heat the part by various mechanisms.

The common procedure for designing an induction heating process in the industry relies on the experts’ experience, trial-and-error and company know-how, since numerical simulation of the induction heating process is not widely spread in the industry because of its complexity and hard-to-obtain material properties [2]. The trial-and-error technique is very time consuming and increases time-to-market and process design costs. Various academics have investigated the induction heating process in the last few years. Although there are several analytical models that represent the electromagnetic behavior of induction systems, the vast majority of the researchers have worked on numerical simulation of the process, mainly because of the geometrical limitations of analytical modeling. Most of the authors have studied linear materials due to the added complexity of non-linear properties.

Bay et al. [3] presented a model for the computation of induction heating of cylindrical billets with long coils. They used an axisymmetrical model, where material nonlinearities were taken into account by the Fröhlich-Kennelly model [4]. A coupled numerical resolution of the electromagnetic—thermal—mechanical problem was proposed by the authors, proving that computational efficiency with adequate precision can be met with their developed model.

Other authors such as Kennedy et al. [5, 6] presented several works where they propose and use mainly analytical models for the evaluation of the induction heating process. For the numerical validation of the proposed model, the authors used the commercial software COMSOL 2D in their work. They studied the effects on short coils and the applicability of classical coil design methods to these kinds of inductors, concluding that 1D analytical solutions for low-frequency air-core induction show good agreement with numerical and experimental results. However, induction heating of nonlinear ferromagnetic materials was not studied, as well as the applicability of the proposed model for multi-turn coils. In their follow-up work [6], the authors investigate the previously presented design methods for induction heating of aluminum billets, comparing the experimentally measured power with calculated values using these methods, concluding that the most accurate results are provided by the model developed from Davies [7] including the Nagaoka correction coefficient for short coils [8]. However, linear or nearly-linear materials were studied and the accuracy of the models for nonlinear materials was not discussed.

Zhang et al. [9] developed an ANSYS one-way workflow for the computation of scanning induction heating in rolled plates. In their methodology, thermal dependence of electromagnetic properties while quasi-static heating is not taken into account. This hypothesis impacts on the accuracy of predicted temperatures, where the authors found an average error of 17.65 %, with higher values at the end of the heating stage. A coupled two-directional electromagnetic—thermal analysis of a moving inductor was presented by Zabett et al. [10] and Schlesselmann et al. [11], finding a very good agreement between experimental and calculated temperatures in both cases.

From the presented works, it is possible to conclude that most of the authors work on simple geometries such as long coils and cylindrical billets, especially if the system is studied analytically. For numerically solved problems, some authors present a simple unidirectional approach, obtaining large errors at some studies. Also, it can be found that a great number of papers are limited to linear materials such as aluminum, as the evaluation of ferromagnetic materials becomes more complex because of their nonlinearity. When it comes to three dimensional cases, the computational time might become excessive. Thus, 3D computation of induction heating has not been widely used. However, being able to simulate induction heating of complex 3D cases is of great importance for industrial applications, because the simple coil configurations that have been studied in the literature stand far from the complex coil designs that we can find in the industry.

In this paper we deal with the difficulty of solving induction heating problems of nonlinear materials with reduced computational time by the usage of a semi-analytical approach to solve these problems. Semi-analytical approaches have been used for many problems (such as plasticity [12], structural mechanics [13] or fracture mechanics on composites [14]) as an effort to reduce computational cost. However, to the authors’ knowledge, there is no proof of the development and usage of semi-analytical approaches for the electromagnetic-thermal problem.

The proposed methodology is validated experimentally for the case of a cylindrical bar heated by a multi-turn coil. The material for the cylinder is 42CrMo4 (also known as AISI 4140), a ferromagnetic low alloy steel commonly used for induction hardened components. An industrial 3D case is presented, where applicability of the developed approach is demonstrated for the simulation of the induction heating stage of large size pitch bearings, which are typically hardened using the induction hardening process. This work is motivated by the size and complexity of the industrial cases such as the presented one, where good computational efficiency is critical for product development.

Modelling of the induction heating process

In induction heated workpieces, the alternating current density is not homogeneously distributed within the cross-section of the billet, as shown in Fig. 1, where the currents are localized in the closest area to the outer surface. This occurs as a result of the so-called skin effect and the exponential decay of the current density can be described analytically [1]

$$\begin{aligned} \vec {J}(\vec {n})=J_S \ e^{-\vec {n}/\delta } \end{aligned}$$

(1)

Fig. 1

Exponential decay of the current density as a result of the skin effect in a cylindrical AC conductor

This expression is valid for plane or curved billets where the radius of the billet is considerably greater than the skin depth [7] and does not considerate edge effects. The skin depth theoretically depends on the frequency of the current and the material properties of the workpiece and can be calculated as [6]

$$\begin{aligned} \delta = \sqrt{\dfrac{1}{\sigma \pi \mu _0 \mu _r f}} \end{aligned}$$

(2)

From the current density distribution inside the workpiece, it is possible to estimate the volumetric heat generation. Ohmic and hysteresis losses are the mechanisms that generate heat inside the workpiece. The later is represented by the enclosed area in the magnetic hysteresis curve. However, for high frequency systems, hysteresis losses are usually neglected in the calculations because of its computational complexity and low impact in the overall losses [1, 15]. The primary cause of heat generation for these cases are the Ohmic losses. The transfer of electric energy into heat is known as Joule effect and produces the so-called Ohmic losses. The ohmic losses are related to the electric conductivity and the current density as [15]

$$\begin{aligned} {\dot{Q}}=\dfrac{1}{2\sigma }\left| \vec {J} \right| ^2 \end{aligned}$$

(3)

In this expression $\vec {J}$ is the peak current.

The aforementioned material properties are mainly temperature-dependent, but might also depend on chemical composition, microstructure or stress state of the material. The change of these properties during induction heating has been described by several empirically-obtained relations.

Ferritic microstructure in low alloy steels such as 42CrMo4 is generally ferromagnetic, thus, the magnet dipoles are oriented in the same direction, reacting strongly with the applied magnetic field. When the microstructure becomes austenitic, its ability to magnetize decreases abruptly and the material becomes paramagnetic, which is weakly attracted by the magnetic field. The transition between ferromagnetism and paramagnetism occurs at the so-called Curie temperature, which is usually between $700\,^{\circ }\hbox {C}$ and $800\,^{\circ }\hbox {C}$.

The magnetic permeability describes the response of a material to an applied magnetic field. It is defined as the derivative of the magnetic field with respect to the magnetic strength ($\mu =\mu _0\mu _r=d|\vec {B}|/d|\vec {H}|$). In the case of steel, the relative magnetic permeability varies enormously with the temperature and the applied field. See Fig. 2 for experimentally obtained magnetic hysteresis curves for 42CrMo4, shown as the averaged curves over a quarter of a period.

The Analytical Saturation Curve model describes the magnetization curve for non-linear materials, and is often presented as an alternative for the classical Fröhlich-Kennelly model, which does not properly describe the saturation of the materials at high magnetic fields [15,16,17]

$$\begin{aligned} B\left( \left| H \right| ,T \right) = \mu _0\left| H \right| + \dfrac{2B_S}{\pi }\arctan \left[ \dfrac{\pi \mu _0 \left( \mu _{r0}-1 \right) }{2B_S} \left| H \right| \right] \left[ 1- e^{\left( \dfrac{T-T_C}{C}\right) } \right] \end{aligned}$$

(4)

Fig. 2

Measured average magnetic hysteresis curves (symbols) at different temperatures and fitted values (hashed lines) for 42CrMo4

The fitted coefficients for 42CrMo4 are $B_S=1.32$ T, $\mu _{r0}=1860$, $33.6^{\circ }\hbox {C}$ and $T_C=783^{\circ }\hbox {C}$ using the nonlinear least square technique for the measured data set. Details of the measurements are given in [16].

The electrical resistivity, which is reciprocal to electric conductivity, is also an important material property that must be well described in order to obtain accurate results. The electrical resistivity of steel depends on the temperature and can be calculated as [1]

$$\begin{aligned} \rho _{\mathrm {W}}\left( T \right) =\rho _{\mathrm {T}_{\mathrm {ref}}}\left[ 1+\alpha \left( T-T_{ref} \right) \right] \end{aligned}$$

(5)

Common values for the resistivity coefficient $\alpha$ can be found at room temperature and is usually assumed to be a constant for the sake of simplicity. In this study, a value of $\alpha =0.005$ is used.

A model for determining the electrical resistivity at room temperature is used in reference [18], which depends on the chemical composition of the steel, provided in mass percentage.

$$\begin{aligned} \begin{aligned} \rho _{\mathrm {T}_{\mathrm {ref}}}=10^{-6} (0.001+0.283\,\mathrm {C} + 0.17\,\mathrm {Si} + 0.0387\,\mathrm {Mn} - 0.1295\,\mathrm {S} + 0.0702\,\mathrm {Al}\\ + 0.00272\,\mathrm {Cr} + 0.0335\,\mathrm {Cu} + 0.0333\,\mathrm {Mo} + 0.0193\,\mathrm {Ni}) \end{aligned} \end{aligned}$$

(6)

Numerical simulation using the finite element method

The electromagnetic phenomena is described by a set of very well known differential equations called Maxwell’s equations [19]. These governing equations describe electric and magnetic fields from the generation of the magnetic field due to the currents in the coil to the induced ones inside the workpiece. By assuming that the displacement current is negligible, the combination of Maxwell’s equations gives the diffusion equation that describes the electromagnetic phenomena for the frequency range used in this study [15, 20]

$$\begin{aligned} \sigma \dfrac{\partial \vec {A}}{\partial t}-\nabla \left( \dfrac{1}{\mu }\nabla \vec {A}\right) =\vec {J}_s \end{aligned}$$

(7)

The common approach to solve Eq. (7) is to utilize the harmonic approximation, which assumes that the source current in the coil is time-harmonic, usually sinusoidal, therefore forcing the response to be time-harmonic as well. Therefore, we can simplify the diffusion equation equation (7) to

$$\begin{aligned} i \omega \sigma \vec {A}_0 -\nabla \left( \dfrac{1}{\mu }\nabla \vec {A}_0\right) =\vec {J} \end{aligned}$$

(8)

The harmonic approximation is used in many commercial eletromagnetic softwares. This steady-state solution is only valid if the magnetic permeability of the material is linear, which is not true for ferritic steels.

Since the relation between the electromagnetic field and flux density is nonlinear for ferromagnetic materials, several adaptations need to be introduced in the calculation of the required linear magnetic permeability. One of the methods to perform this is the calculation of an effective permeability by a fictitious linear material, as developed by [4] and used in several works such as [15, 16, 21]. The basis of the effective permeability method is that the linear fictitious material has the same eddy current average loss density as the true nonlinear material. The effective permeability at point i can be calculated as

$$\begin{aligned} \mu ^{e}_i=\dfrac{\mathrm {w}_{1i}+\mathrm {w}_{2i}}{(H^{e}_{mi})^2} \end{aligned}$$

(9)

where Eqs. (10) and (11) correspond to the upper and lower bounds of magnetic co-energy density of the real material, respectively

$$\begin{aligned} \mathrm {w}_{1i}= & {} \int ^{H^{e}_{mi}}_{0}{B\mathrm {dH}} \end{aligned}$$

(10)

$$\begin{aligned} \mathrm {w}_{2i}= & {} \dfrac{1}{2}B_{mi}H^{e}_{mi} \end{aligned}$$

(11)

The fictitious material is locally and instantly linear, thus its magnetic permeability changes from point to point according to local temperature and applied field. For a better understanding, Fig. 3 shows a real magnetization curve (blue dashed line), its linearization (red dotted line) and the linear effective permeability of the fictitious material (black solid line), fitted so that its magnetic co-energy (in orange) is equal to the one of the real magnetization curve (in blue).

Fig. 3

B–H curve representation for real (blue dashed line), liniarized (red dotted line) and liniarized fictitious material (black solid line) [4, 15, 16]. Note: The reader is referred to the web version of this paper for interpretation of the references in color

Once the electromagnetic heat loss is computed, the system must be analyzed thermally. The workpiece faces three thermal mechanisms: conduction of heat inside the workpiece from the hot surface through its colder core, convection from the surface to the surrounding media and radiation. The heat conduction equation and details about its computation are not shown in this work for the sake of brevity.

Because of the nonlinear and coupled nature of the induction heating problem, several algorithms have been developed for the computation of the electromagnetic–thermal coupling. There are three main approaches to perform a multi-physical analysis; the simplest method is the unidirectional two-step approach, where the distribution of the heat source is calculated once and introduced to a thermal model, where the temperature profile is obtained. In this approach, the temperature and field dependence of electromagnetic properties are not considered. Thus, this approach is limited to low-temperature heating of linear or nearly-linear materials such as aluminum or copper. The most used approach to couple electromagnetic and thermal models is the indirect or staggered coupling method [22], where the electromagnetic diffusion equation is usually solved using the harmonic approximation, where the permeability has been linearized. Similarly to the two-step approach, both analyses are performed separately. However, there is an iterative process so that both calculations are performed at each step, enabling the nonlinearities and thermal-dependent electromagnetic properties to be taken into account [1]. The third method for coupling electromagnetic and thermal calculations is the direct or fully coupled approach, where the set of differential equations (including the electromagnetic diffusion equation (7)) is solved simultaneously and there is no need to use the harmonic approximation. This approach requires an intensive computational time and memory allocation and is not commonly used in finite element programs [22].

The main advantage when using FEM is that different geometries can be easily evaluated, opposite to the geometrically limited analytical equations. However, FEM calculations require much longer computational times. Authors such as Kolanska-Pluska [23] state that, for finite coils, numerical evaluation is required and only infinite coils can be analytically evaluated. It is important to mention that, when performing an electromagnetic–thermal coupled calculation for ferromagnetic materials, the simulation process should be two directional or iterative because of the material non-linearity. This means that, for each time step, the electromagnetic calculation needs to be re-evaluated due to the temperature-dependent material properties, which increases the computation time.

Developed semi-analytical approach

In the developed methodology, the advantages of analytical and numerical calculations have been brought together. Electromagnetic calculation by the use of finite elements allows an easier evaluation of the billet and inductor geometries, even when the system includes complex geometries and auxiliary elements such as magnetizers. On the other hand, analytical equations offer very high calculation speeds.

The presented semi-analytical approach is divided into several steps, see workflow in Fig. 4. First, the governing electromagnetic equations are evaluated numerically by commercial finite element software ANSYS Maxwell, Release 2019R1 [24]. Solving a time-harmonic electromagnetic FE analysis using equation (8) and the linearized magnetic permeability shown in equation (9) enables us to determine the fields generated by the inductor at room temperature according to the input current, the geometry of the billet and the inductor, the coil-workpiece configuration and the surrounding media. The magnetic field in the workpiece at room temperature is extracted from this analysis, referred as $H_{\mathrm {T0}}$ in the flowchart. This is the only finite element computation regarding the electromagnetic physic. Before the semi-analytical iterative process is started, the magnetic field at the surface of the billet is extracted, and the normal vector is computed for each surface node.

Fig. 4

Flowchart of the proposed methodology

Once the FE has solved the magnetic field, an iterative semi-analytical process is started. Here, only the workpiece geometry is modelled The remaining electromagnetic and thermal physics are resolved analytically and numerically, respectively.

In every time increment, the analytical process starts with the calculation of the electromagnetic properties of the material, which depend on the temperature and the magnetic field. The material properties are calculated for each node, since the heating rate is highly variable with its position. First, Eq. (5) is evaluated for the electrical resistivity of the node. For the evaluation of the magnetic permeability, the slope of the magnetization curve can be evaluated by the derivative of Eq. (4), allowing us to use the true magnetic permeability, unlike the mostly used commercial softwares using the harmonic approximation and the linearization of the magnetic permeability. Following, the current drop-off is computed in the normal direction of each surface node using Eq. (1), obtaining the distribution map of the current density inside the workpiece, computed in the local mesh that is used for the analytical solution. As the equation works on the normal direction from every point on the surface, this equation can be applied indistinctly for every workpiece geometry, from simple axi-symmetric cylinders to complex 3D billets, because of its unidirectionality. When put together, the computed unidirectional current drop off offers a three dimensional map in the local mesh generated for the analytical solution.The usage of Eq. (1) implies that the electromagnetic material parameters are kept constant throughout the sub-surface region. However, thermal properties are considered to be temperature dependent in this region.

Once the current density distribution inside the workpiece is computed, the heat generated in the billet from the eddy currents is calculated using Eq. (3). The evaluated heat generation distribution is mapped into a transient FE thermal model developed in ANSYS, where the temperature distribution at the end of the time step is computed. The temperature map is brought back to the analytical electromagnetic solution and the material properties are updated, as these affect the magnetic field map inside the workpiece and, thus, the distribution of the generated heat. This computation loop follows the staggered coupling approach. As it can be seen in the flowchart, the electromagnetic FE analysis is not included in the computation loop, but acts only as a generator of input data, considerably reducing the computational time associated with electromagnetic FE solvers.

During heating, the initial ferritic microstructure is transformed into austenite when the critical equilibrium temperature $A_{e1}$ is reached, allowing austenite nucleation to start. This process is finished at the austenite end temperature $A_{e3}$, where all the $\alpha$ ferrite structure is transformed into $\gamma$ austenite. Subsequent quenching of the workpiece enables the parent structure to transform into other microstructures depending on the cooling rates. A very fast cooling enables martensite to form, a very hard but brittle microstructure that is usually the main objective of the hardening processes. During induction heating, the goal is to achieve an austenized layer that can be later transformed into a hard martensite case.

There are several methods to predict the austenitization process. One of the simplest models that require no empirical dependency is the linear equation as described by [25]

$$\begin{aligned} f_{\gamma }=\dfrac{A_{e1}-T}{A_{e1}-A_{e3}} \end{aligned}$$

(12)

In this work, critical equilibrium temperatures $A_{e1}=740^{\circ }\hbox {C}$ and $A_{e3}=805^{\circ }\hbox {C}$ have been used [26].

Fig. 5

Normalized magnetic field strength vs. temperature for various geometries, frequencies and current intensities with fitted Thermal Modifier in solid black line (right) computed in the central point marked in red (left). Note: The reader is referred to the web version of this paper for interpretation of the references in color

One difficulty when dealing with the uncoupled electromagnetic FE computation is that the magnetic field strength in the surface of the workpiece ($H_{S_{T0}}$) varies with temperature and cannot be assumed constant. A numerical Design of Experiments analysis was carried out for frequency (2.5, 5 and 7.5 kHz), current intensity (2.5, 5 and 7.5 kA) and geometrical (cylinders and tubes with different coupling distances, named 1 to 3) variables, observing that the evolution of the surface field strength with temperature varies in a similar manner for different superficial points at every case. Figure 5 (right) shows the evolution of the normalized magnetic field with respect to its initial value for each studied case. Only the values at the center of the surface, marked with a red dot in the left figure, are shown in the graph for the sake of simplicity.

In the figure it is possible to observe that, when temperature increases, the value of the magnetic field starts to drop. This effect occurs as temperature increases and gets closer to the Curie temperature, when the material properties abruptly change and the magnetic field suddenly increases. Once its highest value is reached at Curie temperature, the magnetic field decreases down to its initial value. This effect happens because the studied points are not isolated and the effect of the surrounding points must be taken into account. When the points at the surface become paramagnetic, the sub-surface area is still ferromagnetic and therefore still able to react with the magnetic field. The drop on the magnetic field above Curie temperature occurs as a result of the surrounding points gradually becoming paramagnetic. Thus, the magnetic field follows the shape of a non-monotonic function. For the developed methodology, a so-called Thermal Modifier ($TM=f(T)$) is introduced in the analytical model in order to update the superficial field strength ($H_s(T)=Hs_{T0} \cdot TM$), as can be seen in the flowchart presented in Fig. 4. A nonlinear least-square fit has been carried out to define the Thermal Modifier in this work, which has been divided into three regions depending on temperature.

$$\begin{aligned} TM= {\left\{ \begin{array}{ll} 4.65\times 10^{-7} T^2 - 5.96\times 10^{-4} T + 0.99 ,&{} \text {if } T<600\\ 0.26 \exp {\left( -\left( \dfrac{T-T_C}{22.7}\right) ^2\right) }+5.25 \exp {\left( -\left( \dfrac{T-6360}{4202}\right) ^2\right) },&{} \text {if } 600\le T<T_{C}\\ 5.17\times 10^{-6} T^2 - 1.04\times 10^{-2} T + 6.15, &{} \text {if } T\ge T_{C} \end{array}\right. } \end{aligned}$$

The presented empirical Thermal Modifier can be used for any simulation carried out in 42CrMo4 workpieces of similar geometry, as the limitation for using this TM resides mainly on the material properties. For materials with significantly different magnetic properties from 42CrMo4, especially the Curie temperature, the Thermal Modifier presented in this work should be recalibrated.

Advantages and limitations of the proposed approach

The presented approach offers the possibility of reaching a good computation efficiency for simulating the induction heating process in complex industrial cases, which usually require large thee dimensional models. However, there are some limitations to this model, as listed below:

• The presented approach has been developed and validated for static induction heating. However, the strategy is valid for the stationary region during scan hardening, provided that the geometry is homogeneous in the studied section, without any cross-section change or any other field concentrating geometrical features such as part edges, holes or grooves for which the initially obtained magnetic field map might not be valid.

• Due to the limited application of equation 1 to surfaces with great curvature to skin depth ratio, the proposed approach special attention should be payed for workpieces with sharp edges or small radii.

• The simple analytical equations used in the proposed approach imply that the electromagnetic material parameters are constant in the sub-surface region. This fact might affect the prediction of the hardened case especially for deep hardened cases.

• The obtained Thermal Modifier has been calibrated for regular cylinders and 42CrMo4 and other geometrical effects that might affect the field distribution inside the workpiece or other materials with a considerably different Curie temperature have not been taken into account.

Numerical and experimental validation of the proposed methodology: a 2D case

An experimental set-up was constructed to validate the results obtained by the numerical models. The coil was built using copper tubing to allow its refrigeration by running water. The cylindrical billet was centered inside the coil. An oscilloscope with a differential voltage probe and a Rogowski current transducer were used to measure the voltage and current through the coil. Two wire type-K thermocouples were spot-welded to the surface at the center of the specimen and at a point $10\,$mm below, as shown in Fig. 6.

Fig. 6

Experimental set-up (right) and position of thermocouples (left)

Fig. 7

Temperature (above), current (lower left) and frequency (lower right) evolution during the induction heating experiments represented as mean value and standard deviation

Fig. 8

Mesh of the semi-analytical electromagnetic (above) and thermal (below) analyses

Fig. 9

Experimental and numerical evolution of temperature at the central point

Ten experiments were carried out for the same configuration, for which mean value and standard deviation have been computed. Figure 7 shows the evolution of temperature at the central thermocouple (T1), frequency and current during the induction heating experiments.

In order to obtain the data shown in Fig. 7, the specimens were heated up to approximately $780^{\circ }\hbox {C}$ and held for 3 to 5 s after turning off the current before cooling in a water tank. The initial transient part of the induction heating system has been omitted in the figures for sake of simplicity

The numerical model used for the semi-analytical approach is composed of 31k triangular axi-symmetric elements for the electromagnetic simulation, where the workpiece, the inductor and the surrounding air are meshed, and a quad-dominated linear mesh composed of 7k nodes for the thermal analysis, where only the workpiece is modeled. The models are shown in Fig. 8. The computation was performed on a personal computer equipped with an Intel(R) Core i7-8650U processor (1.9 GHz) and the typical CPU time for the semi-analytical solution was 1.3h.

The model constructed for the Flux software is composed of a 138k triangular elements, including the workpiece, the inductor and the surrounding air. The typical CPU time for this computation was 6.2h.

The results are shown in the comparison graph of Fig. 9, where the evolution of temperature at the central superficial T1 point (thermocouple T1 in Fig. 6) is shown for both numerical simulations and the experimental measurement. The shaded gray area indicates the experimental variability, where the mean value is shown as a black solid line.

In the presented graph, it is possible to observe two stages during heating. First, the temperature raises at a nearly constant heating rate, but as the temperature gets close to Curie temperature, a considerable change in the heating rate can be observed, demonstrated by a knee on the curve. Thus, it is possible to simplify the evolution of temperature as a bilinear curve, that varies as the magnetic permeability diminishes when it gets closer to Curie temperature. After Curie temperature, the heating rate is very slow as the material is now paramagnetic. The bilinearity in the temperature evolution can be observed for both numerical models, but the transition between the stages is not precisely described in the model using the Flux software as it appears later and at a higher temperature when compared to experimental data.

The presented semianalytical model presents very good agreement with the experimental data in both stages, where the intersection between the heating stages is properly described, although the stability at the second zone is reached after certain time, probably as an effect of a great magnetic permeability gradient within the surrounding area.

Fig. 10

Experimental and simulated hardened case for fully hardened region in red and end of transition zone in green (left) and Vickers Hardness measurement and simulated austenite fraction in the radial direction (right)

Figure 10 shows the experimental and simulated hardened case pattern (left), as well as a comparison of the experimentally measured hardness through the radius and the simulated austenite fraction at the end of the heating stage (right). It can be assumed that the reached austenite will transform into martensite as a result of the fast cooling rate applied to the induction heated workpieces, thus, giving a hard layer. In the figure it is possible to observe that the hardened case has been properly estimated in the region where a fully transformed microstructure can be found (marked in red in the experimental figure). The results also show that the transition region between hardened and nonhardened zones (marked in green) has been underestimated in the simulation. The consideration of constant electromagnetic properties throughout the sub-surface region might explain this underestimation, as the transition zone is mainly driven by conduction while the fully austenized area is driven by the internal heat generation.

When it comes to an industrial application such as induction hardening, ensuring a proper prediction of the hardened depth is vital, thus the great importance of describing the temperature evolution properly, especially for the second stage, where austenitization of the material occurs at temperatures between 700 $^{\circ }$C and 800 $^{\circ }$C for most low alloy steels. An error on the temperature prediction can impact the resulting microstructural transformations. In the temperature evolution for the Flux software shown on Fig. 9 it is possible to observe that austenitization would occur much earlier than experimentally measured because of the inaccurate prediction of the transition area, deriving into an improper prediction of microstructural transformations and affecting the estimation of hardened depth, as there might be sufficient time for conduction-driven heating in the subsurface area to happen, predicting a deeper heat affected zone when compared to experimentally obtained depths.

3D case: a large size pitch bearing

An industrial 3D case has been constructed in order to demonstrate the applicability of the proposed simulation strategy and the potential time reduction that can be achieved by the use of the developed approach.

Fig. 11

Typical pitch bearing geometry (left) and its cross section with the area of study shaded in red (right). The black thick lines indicate the hardened surface

Pitch bearings connect the blades to the rotor in wind turbines, allowing the oscillation of the blades to be controlled. The most common design for these kinds of bearings is the four-point contact double row ball bearing, as shown on Fig. 11. Because of the mechanical requirements for these components, their races are usually hardened by induction, with a typical hardened depth of 5 to 8 mm. A common strategy to harden the races is to perform the induction hardening operations separately for each race. Thus, the area of study has been simplified and corresponds to a single row of the shown bearing and has been shaded in red in the cross section.

As a result, the studied geometry is a section of a single row of the bearing. Figure 12 shows the constructed model, including the inductor and the ferrite piece, which acts as a field concentrator that allows the usage of lower power systems, improving costs and environmental impact. The air is also modeled in the electromagnetic FE computation, but it has been suppressed in the figure for the sake of clarity. The input current for the inductor has a peak value of 3000 A with a frequency of 20 kHz. The inductor is followed by a cooling shower, modeled as a moving region of high convection, where a convection coefficient of 15000 W/m$^2$K is used to simulate a severe forced convection produced by the impact of the cooling shower onto the hot surface. To solve the temperature field distribution and compute the area of austenite inside the race in the workplace for a total heating time of 35 s with a scanning speed of 4 mm/s took 30 h using 4 cores in a workstation equipped with an Intel(R) Xeon(R) Gold 6242 CPU dual processor (2.80 GHz). The solver was a direct solver and parallelized. To calculate the initial electromagnetic field took 1 h. The coupled electromagnetic-thermal computation took 87 h to complete in the commercial software Flux, where the mesh was composed by 420k nodes, 130k of which correspond to the workpiece.

Fig. 12

Objects of study (a), generated magnetic field (b) and FE mesh (cross-section in c and full body in c), composed of 624k linear hexagonal and triangular elements, for the semi-analytic simulation. Points of study marked in red (d)

Fig. 13

Temperature evolution (left) and formation of austenitic phase (right) during scanning induction hardening of the pitch bearing

In Fig. 13a the evolution of temperature and austenite fraction during the scanning induction hardening simulation is shown. In the figure it is possible to observe the movement of the heat source and the cooling region. It can be concluded that after approximately 15s of heating, the system can be considered stable and the computed heat source behaves stationary, reaching a homogeneous hardened case for the full bearing.

Fig. 14

Comparison of temperature evolution at different points as computed by the semi-analytical model and the commercial software Flux

In figures (b) and (c) the transformed austenite fraction can be seen in the cross-section of the workpiece, where it is possible to observe the austenized depth inside the race. Figure (b) shows the austenite case depth computed using the proposed approach, while Figure (c) shows the temperature computed by Flux. In this case, the temperature above austenite end temperature has been plotted in white and corresponds to the hardened case. The austenized depth is approximately 5.5 mm in the homogeneous region, which is within the reported experimental span. It should be emphasized that we do not simulate the martensitic transformation that occur during quenching. Instead, it is assumed that the layer of austenite is transformed into martensite. A more homogeneous heating pattern might be obtained with further research on the design of the inductor.

In Fig. 14, a comparison between the temperature evolution is shown in the study points marked in Fig. 12d. In the figure, it can be observed that the temperature evolution has correctly been computed with reduced computational time, validating the proposed methodology in terms of accuracy and computational time.

The presented case study shows that the proposed simulation strategy can be satisfactorily used for scanning three dimensional cases with a good computation efficiency and accuracy. Further development on the algorithm will provide the possibility of incorporating the computation of residual stresses after the induction hardening process. With these extended capabilities, a deep study on the induction hardening process and its effect on component performance can be performed for large-size bearings.

Conclusions

A semi-analytical model that describes the electromagnetic-thermal coupling during induction hardening of ferromagnetic materials has been presented in this work. The model has been validated against experiments for 42CrMo4 steel cylinders, which shows very good agreement between modeled and measured surface temperatures, obtaining an average error between 0.9% and 4.1%. Experimental hardening depths have been compared to the simulated austenite fractions with good correlation, fully validating the developed approach. The results provided by the developed model have been compared with the commercial FE software Flux, which uses a harmonic approximation and linearized magnetic permeability to compute the electromagnetic field. The presented semi-analytical model resulted in more accurate results and a higher computational efficiency compared with the harmonic approximation solution strategy, reducing computational time by about 80% for the axi-symmetric case. Also, we demonstrate that the developed semi-analytical model can be used to model complex 3D cases with a 50% time reduction. The simulated thickness of austenite case depth is within the span for experimental martensitic depth and has been compared to the results provided by Flux along the temperature evolution during the scanning process, validating the proposed approach in terms of accuracy.