MODELLING OF CONTINUAL INDUCTION HARDENING
IN QUASI-COUPLED FORMULATION
Jerzy Barglik1, Ivo Doležel2, Pavel Karban3 and Bohuš Ulrych3
Silesian University of Technology, Department of Electrotechnology,
Krasinskiego 8, 40-019 Katowice, Poland
Academy of Sciences of the Czech Republic, Institute of Electrical Engineering,
Dolejškova 5, 182 02 Praha 8, Czech Republic
University of West Bohemia, Faculty of Electrical Engineering,
Sady Pětatřicátníků 14, 306 14 Plzeň, Czech Republic
ABSTRACT. Continual induction hardening of steel workpieces belongs to modern and en-
vironmentally friendly industrial technologies. The body is heated by an inductor moving at a
low velocity along it and the heated part is then immediately cooled by spraying quenchant.
Modelling of this process is an uneasy business because of two reasons. The first one is the
presence of the moving inductor and the second one problem with determining the coefficient
characterizing the convection heat transfer between the body and a mixture of air and spray.
On the other hand, the results must be as accurate as possible in order to predict the resultant
hardness. The paper presents a methodology where the induction heating (time and space
variable distribution of the electromagnetic and temperature fields) is modelled in a slightly
modified standard manner while the process of cooling by an approximate theoretically em-
pirical procedure based on experimental data. The theoretical analysis of the task is supple-
mented with an illustrative example and discussion of the results. Computations are per-
formed by combination of a professional code and a number of single-purpose user programs
developed and written by the authors.
Continual induction hardening represents a technological process that is frequently used in
various industrial applications. The processed workpiece is heated by an inductor carrying
harmonic current that moves in common with a sprayer at a low velocity along it and the
heated part of the body is immediately cooled by spraying quenchant (water, oil, polymer-
based liquids etc.). The result is harder, but more fragile martensite structure. The scheme of
the process is shown in Fig. 1.
The final hardness of the surface layers of the workpiece is generally a complex function
of velocity (or, equivalently, time) of its heating and cooling. This is partly explained in Fig. 2
containing the unbalanced phase diagram of a typical carbon steel for a given austenitizing
temperature Ac3 .
First, the part of the body to be hardened has to be heated over the temperature Ac3 in or-
der to reach uniform austenite structure (domain A) and then cooled. The decreasing lines in
Fig. 2 starting from the temperature Ac3 and ending at the final temperature TF represent
various regimes of cooling. The first two lines on the left indicate, for example, relatively fast
cooling leading to complete hardening of the processed material; the resultant structure con-
taining practically pure martensite is characterized by high hardness given (in Vickers) within
the corresponding circles. The lines on the right, on the other hand, express slower cooling
resulting in incomplete hardening (the final structure contains not only martensite M, but also
undesired bainite, pearlite or ferrite whose zones in the diagram are denoted by letters B, P
and F, respectively).
part hardened steel
Figure 1. The investigated arrangement
600 673 409
100 HV = 765 762 742 505 432 279 243 TF
1 10 102 103 104 105
Figure 2. Unbalanced phase diagram of Polish steel NZ3 (0.6% C, 1.0% Si, 1.0% Cr,
2.0% W and 0.2% V), hardness given in Vickers, Ac3 = 880 °C, TF = 100 °C
The task represents a coupled electromagnetic-thermo-metallurgical problem characterized
by mutual interaction of several physical fields with time-dependent boundary conditions that
can be (particularly in the case of temperature field) determined only indirectly. Moreover, the
material parameters are temperature-dependent functions that must be found experimentally.
The first mathematical and corresponding computer model for solution of the problem was
suggested by the authors in . But even when being relatively sophisticated, the
methodology suffered from several simplifications enforced by lack of experimental data and
possibilities of available software.
Mentioned can be
• weakly coupled formulation of the problem (electromagnetic field was supposed to be
independent of the temperature field),
• step change of magnetic permeability when the temperature reached the Curie point,
• temperature Ac3 was considered constant (history of heating was not taken into account),
• step variation of the convection heat transfer coefficient from the surfaces of the heated
and cooled parts to ambient medium,
• rough mesh used for thermal computations (only several hundred elements).
Since then, however, the methodology of solution was substantially improved, together
with our own procedures (that allow, nowadays, finite-element analysis of an arbitrary axi-
symmetric problem of this kind in the quasi-coupled formulation). The paper describes the
algorithm and illustrates it on an example.
DESCRIPTION OF THE TECHNICAL PROBLEM
The inductor (Fig. 1) carrying harmonic current I ext of density J ext and frequency f
moves at a velocity v along a steel mandrel in common with a sprayer with cooling medium
under pressure that assures consequent fast cooling of the heated part.
The first task is to recommend the parameters of the inductor (geometry, current, fre-
quency and velocity v of its motion) that provide sufficient temperature rise of the heated part
(somewhat higher than Ac3 ). The second task is to find such pressure (or velocity) of the
spraying quenchant that the convection of heat from the body to ambient mixture guarantees
sufficiently fast cooling and required surface hardness.
While the first task is processed theoretically and realized by a suitable numerical algo-
rithm, the second task is processed by an experimentally-empirical procedure. This procedure
is, nevertheless, realized again in the numerical manner.
MATHEMATICAL MODEL AND NUMERICAL ALGORITHM OF ITS SOLUTION
The mathematical model consists of two partial differential equations describing the distri-
bution of the electromagnetic and temperature fields with temperature-dependent coefficients.
The investigated arrangement is considered fully axisymmetric, so that the definition area
Ω in Fig. 3 is rectangular (ABCDA) with artificial boundary BCDA sufficiently distant from
the workpiece and inductor. This area consists of domains Ω 1 (mandrel), Ω2 (moving induc-
tor) and Ω3 (air). Relative permeability μ r of the workpiece is supposed to be only a tem-
perature-dependent function, which simplifies the situation. The field quantities may be now
replaced by their phasors and the field distribution is described by the Helmholtz equation 
rot rot A + j ⋅ ωγμ A = μ J ext (1)
where A denotes the phasor of the vector potential, γ the electrical conductivity, μ the
magnetic permeability and ω = 2πf . Each domain is characterised by specific values of
phasor J ext , γ and μ , the last two parameters being supposed to be functions only of the
temperature T . In fact, phasors of the current density as well as the magnetic vector potential
have only one nonzero tangential component J ext,ϕ and Aϕ , respectively. The Dirichlet con-
dition along the artificial boundary BCDA reads Aϕ = 0 . This boundary has to be taken far
enough from the arrangement in order to avoid any influence on electromagnetic field distri-
bution in it.
Its definition area is represented by the longitudinal axial cut through the workpiece
(KLMNK in Fig. 3) and its distribution is given by the Fourier-Kirchhoff equation 
div ( λ ⋅ grad T ) = ρ c ⋅ − wJa , (2)
where λ denotes the specific heat conduction of the material, ρ its specific mass, c its spe-
cific heat and wJa the specific average Joule losses
wJa = γω 2 A . (3)
air Ω3 z
Figure 3. Definition area of the problem: S – starting position of the inductor,
G – general position of the inductor, E – end position of the inductor
The boundary condition for t > 0 along line KL reads ∂T / ∂r = 0 while along line LMNK
−λ ⋅ ∂T / ∂n = α (T − T0 ) . The influence of radiation above the sprayer is neglected at this
stage. Distribution of the convection heat transfer coefficient and temperature of the ambient
medium depends on the instantaneous position of the inductor with the sprayer and is depicted
in Fig. 4. While α w (coefficient for heat transfer from the workpiece to cooling medium) is
usually much greater than α a (analogous coefficient from the workpiece to ambient air), tem-
perature T0w of the coolant is, on the contrary, lower than temperature of air T0a .
The initial condition ( t = 0 ) for the temperature field within the whole region Ω is Tstart .
Now it is necessary to determine the coefficient α w characterizing the convection heat
transfer between the workpiece and spraying quenchant. As its theoretical investigation is
extremely difficult, it was determined by an indirect approximate theoretically empirical
method based on experimental results. This method suggested by the authors is described in
details in . In fact, α w is a function of the surface temperature, but the method provides
only its mean value α wm .
position S position G
α T0 α T0
surface of the heated workpiece
αa T0a αa T0a
Figure 4. Distribution of α and T0 along the mandrel for starting (S)
and general (G) positions of the inductor
Distribution of hardness
Dependence of hardness on the time of cooling tc from temperature Ac3 to TF was re-
drawn from Fig. 2 to Fig. 7b. This dependence was built into the code. When the highest tem-
perature of any point does not reach Ac3 , computation of hardness is not realized there.
1. Choice of a time step Δt1 for which the inductor with the sprayer are in one position (their
motion is not supposed to be continuous, but discrete, by small shifts Δl ). The time step
Δt1 = 2 s (corresponding to change of the position of the system inductor-sprayer 1 mm)
decisive for changes of electromagnetic field proved to be small enough to achieve re-
quired accuracy of the calculations.
2. Computation of electromagnetic field, eddy currents and corresponding Joule losses in the
heated body for the actual position of the inductor.
3. Determination of the time-variable boundary conditions for the temperature field by
means of the described theoretically-empirical procedure. The procedure makes it possible
to take into account average value of coefficient of convection heat transfer α wm during
the process of cooling. Computation of the time evolution of nonstationary temperature
field was realized with time steps Δt2 satisfying the condition of numerical stability of so-
lution. Automatic change of material parameters according to actual temperature in indi-
vidual cells. Creation of the database containing history of cooling of cells in the surface
layers of the cooled part of the workpiece.
4. Shift of the inductor by value Δl = v ⋅ Δt1 and return to point 2. Points 2–4 are repeated
until the whole workpiece is processed.
5. Simultaneous calculation of the resultant field of hardness.
The methodology is illustrated on an example of induction hardening of a steel mandrel
(steel NZ3 whose unbalanced diagram is in Fig. 1).
The most important input data
• Several preliminary computations showed that the safe artificial boundary BCDA is char-
acterized by the following dimensions: h1 = h2 = 0.15 m, r1 = 0.2 m. Its enlargement leads
to small changes in results not exceeding about 2%.
• Heated workpiece from steel NZ3: ra = 0.012 m, rb = 0.0205 m, l = 0.610 m. Material
parameters (and eventually their temperature dependencies) are depicted in Figs. 5, 6 and
7. Preliminary computations indicated that the average magnetic flux density in the sur-
face layers of the mandrel is about 0.5 T (Fig. 5a). That is why we used corresponding
relative permeability μ r = 2600 determined for interval B ∈ 0,1.5 T. Its temperature
dependence is in Fig. 5.
• Coil: rc = 0.033 m, w = 0.017 m, h = 0.03 m. Other parameters: J ext = 107 A/m2,
f = 10 kHz, v = 0.0005 m/s. In fact, the real coil made from water-cooled copper tube was
replaced by a thin copper conductor. The main reason of this simplification was to sub-
stantially reduce the time of computations.
• The principal temperature conditions: Tstart = 20 °C, T0a ( LMZ ) = T0a ( KN ) = 20 °C,
T0w ( NZ ) = 10 °C (MZ – above the sprayer, NZ – below the sprayer),
α a ( LMZ ) = 20 W/m K, α wm ( NZ ) = 500 W/m K, α a ( KN ) = 5 W/m K.
2 2 2
2000 average value
0 0,5 1 1,5
Figure 5. Dependencies μ r ( B ) for T = 20 °C – a) and μ r (T ) for B = 0.5 T – b) (steel NZ3)
Selected results and their discussion
Computations were carried out on a computer CPU 2.6 GHz with memory RAM 0.5 GB.
The average time of computation of one case was about four hours. Electromagnetic field was
calculated on a mesh with about 200000 elements (ABCDA, see Fig. 3) while the mesh for
consequent temperature field consisted of more than 10000 elements (KLMNK, see Fig. 3).
The area was remeshed at each position of the system inductor-sprayer (their number was
about 80). Geometrical convergence of the results was tested carefully with respect to both
the artificial boundary and density of the discretization mesh.
0 150 300 450 600 750 900 0 150 300 450 600 750 900
T (°C) T (°C)
Figure 6. Dependencies γ = γ (T ) – a) and λ = λ (T ) – b) for steel NZ3
0 150 300 450 600 750 900 0 500 1000 1500 2000
T (°C) tc(s)
Figure 7. Temperature dependence of specific heat
– a) and dependence of hardness on the time of cooling – b) for steel NZ3
The most important results are depicted in Figs. 8, 9, 10. Fig. 8 shows the smoothed time
evolution of temperature at selected point of the mandrel surface. Fig. 9 depicts, on the other
hand, distribution of temperature along the mandrel for several positions of the inductor. The
curves for t = 1400 s and t = 1600 s were calculated for cooling, when the inductor reached the
end of the mandrel and the current was switched off.
There are small discrepancies between Figs. 8 and 9 caused by local numerical instabili-
ties. These differences not exceeding several percents could surely be reduced by further re-
finement of the mesh elements and time steps, but at the expense of the computational costs.
Moreover, this inaccuracy is comparable with inaccuracy of some input data (or even lower)
and it is disputable whether eventual improvement in this direction would bring more accurate
Finally Fig. 10 shows the distribution of hardness along the surface NM of the mandrel. Its
variations follow from different times of cooling of its individual parts. In fact, however, the
difference between minimal and maximal values of hardness is 16 HV, which does not exceed
about 2.5% from the average value.
Accuracy of computations depends on the quality of the mathematical model, particularly
on accuracy of experimentally obtained temperature-dependent properties of present media.
Very important is also knowledge of hardness as a function of the velocity of cooling. The
computations themselves take a relatively long time (a couple of hours) due to both complex-
ity of the task and verification of the geometrical convergence of results. These are in a rela-
tively good accordance with experiments .
Figure 8: Time evolution of temperature at selected points along the surface of the mandrel
500 207 s
400 720 s
300 1100 s
200 1400 s
100 1600 s
0 0,1 0,2 0,3 0,4 0,5 0,6
position of the inductor with respect to the mandrel (m)
position of the inductor (m)
Figure 9: Distribution of temperature along the mandrel for the positions of the inductor given
by the presented times
Next work in the field will be aimed at more detailed evaluation of the history of heating
(because the value of Ac3 depends on its velocity) and eventual influence of radiation in sys-
tems with moving sprayers. Some progress in acceleration of computation algorithms is ex-
pected as well.
Financial support of the Czech Ministry of Education (Research Plans VZ MSM
232200008 and 212300016) and Polish Scientific Research Committee (grant No. 4T08C
04823) is highly acknowledged. As well the paper was prepared in a framework of multi-
government bilateral collaboration between Poland and the Czech Republic.
0 0,1 0,2 0,3 0,4 0,5 0,6
axial distance (m)
Figure 10: Distribution of hardness along the surface NM of the mandrel
 Doležel, I., Barglik, J., Ulrych, B. (2003). Continual induction hardening of axisym-
metric bodies. Proc. of IC JAPMED’2003, 79–82.
 Chari, M.V.K., Salon, S.J. (2000). Numerical methods in electromagnetism. Academic
 Holman, J.P. (2002). Heat transfer. McGrawHill.
 Barglik, J., Doležel, I., Ulrych, B. (2003). Determination of convection heat transfer
coefficient at continual induction heating of axisymmetric bodies. Scientific Bulletin of
Lodz Technical University, No. 101, 9–16.