Induction heating (IH) is a multiphysical phenomenon comprising a complex interaction of electromagnetic, heat transfer, metallurgical phenomena, and circuit analysis that are tightly interrelated and highly nonlinear because the physical properties of materials depend on magnetic field intensity, temperature, and microstructure. Metallurgical/microstructural specifics of induction thermal processing are reviewed in Chapter 4. This chapter focuses on the electromagnetic and heat transfer phenomena, process simulations, and some other related aspects.
Induction heating (IH) is a multiphysical phenomenon comprising a complex interaction of electromagnetic, heat transfer, metallurgical phenomena, and circuit analysis that are tightly interrelated and highly nonlinear because the physical properties of materials depend on magnetic field intensity, temperature, and microstructure. Metallurgical/microstructural specifics of induction thermal processing are reviewed in Chapter 4. This chapter focuses on the electromagnetic and heat transfer phenomena, process simulations, and some other related aspects.
The basic electromagnetic phenomena of IH have been discussed in several textbooks including college physics. An alternating voltage applied to an induction coil (e.g., solenoid multiturn coil, Figure 3.1) will result in an alternating current (AC) flow in the coil circuit.
Figure 3.1 Conventional multiturn solenoid inductor.
An alternating coil current produces in its surroundings a timevariable magnetic field that has the same frequency as the coil current. This magnetic field induces eddy currents in the workpiece located inside the coil. Eddy currents will also be induced in other electrically conductive objects that are located near the coil. These induced currents have the same frequency as the coil current; however, their direction is opposite to the coil current. These currents produce heat by the Joule effect ( I ^{ 2 } R ). Figure 3.2 shows a sample of a variety of inductor geometries used in IH. Recognizing that there is almost an endless variety of inductor types, it is convenient to review basic principles of IH considering a conventional solenoidtype coil that surrounds a cylindrical workpiece. This approach will be used in this chapter.
Figure 3.2 Sample of a variety of inductor geometries used in IH. (Courtesy of Inductoheat Inc., an Inductotherm Group company.)
Because of several electromagnetic phenomena, the current distribution within an inductor and workpiece is not uniform. This heat source nonuniformity causes temperature gradients in the workpiece. Nonuniform current distribution is associated with several electromagnetic phenomena, including but not limited to
These effects play an important role in the performance of an IH system [1,2,5–13,18–21,44]. Before exploring the factors affecting a magnetic field distribution and eddy current flow, it is imperative to understand the nature of electromagnetic properties of heated materials.
Electromagnetic properties of materials is quite a broad expression that refers to a number of engineering characteristics including electrical resistivity (electrical conductivity), relative magnetic permeability, saturation flux density, coercive force, hysteresis loss, initial permeability, permittivity, magnetic susceptibility, magnetic dipole moment, and many others. Recognizing the importance of all electromagnetic properties, we concentrate in this text only on those properties that have the most pronounced effect on performance of the IH systems.
The ability of material to easily conduct electric current is specified by electrical conductivity σ [45–49]. The reciprocal of the conductivity σ is electrical resistivity ρ. The SI units for σ and ρ are mho/m and Ω*m, respectively. Both characteristics can be used in engineering practice; however, the majority of data books consist of data for ρ. Therefore, the value of electrical resistivity is primarily used in this text.
Metals and alloys are considered to be good electrical conductors and have much less electrical resistance to the current flow compared to other materials (e.g., ceramics, plastics, etc.). Table 3.1 shows values of ρ for common materials at room temperature [66]. Although most metallic materials are known to be electrical conductors, they are, in turn, also divided on several subgroups based on the magnitude of their electrical resistivities. There are metals and alloys that are considered being lowresistive metals (e.g., silver, copper, gold, magnesium, aluminum) and highresistive metals and alloys (e.g., titanium, carbon steel, stainless steel, tungsten, Nibased superalloys).
Material (at Room Temperature) 
Electrical Resistivity (μΩ*m) 

Silver 
0.015 
Copper 
0.017 
Gold 
0.024 
Aluminum 
0.027 
Tungsten 
0.054 
Zinc 
0.059 
Nickel 
0.068 
Cobalt 
0.09 
Mild carbon steel 
0.16 
Stainless steel 
0.7 
Lead 
0.21 
Titanium 
0.42 
Nichrome 
1 
Graphite 
7–9 
Wood (dry) 
10^{14}–10^{17} 
Glass 
10^{16}–10^{20} 
Mica 
10^{17}–10^{21} 
Teflon 
>10^{19} 
Source: V. Zinoviev, Thermal Properties of Metals at High Temperatures, Metallurgia, Moscow, 1989.
Electrical resistivity of metallic materials varies with temperature, chemical composition, microstructure, and grain size. For most metals, ρ rises with temperature. Figure 3.3 shows electrical resistivities of some commonly used materials as a function of temperature.
Figure 3.3 Electrical resistivities of some commercially used metals.
The resistivity of the pure metals can often be approximated as a linear function of the temperature (unless there is a change in a lattice of material/phase transformation)
where ρ_{0} is the electrical resistivity at ambient temperature T _{0}, ρ(T) is the resistivity at temperature T, and α is the temperature coefficient of the electrical resistivity. The unit for α is 1/°C. Table 3.2 consists of the values of α for selected pure metals [67].
Metals (at Room Temperature) 
α (1/°C) 

Aluminum 
0.0043 
Cobalt 
0.0053 
Copper 
0.004 
Gold 
0.0035 
Iron 
0.005 
Lead 
0.0037 
Nichrome 
0.0004 
Nickel 
0.0069 
Silver 
0.004 
Titanium 
0.0035 
Tungsten 
0.0045 
Zinc 
0.0042 
Source: I. Grigor’ev and E. Meilikhov (editors), Physical Values, Energoatomizdat, Moscow, 1991.
For some electrically conductive materials, ρ decreases with temperature and, therefore, the value of α can be a negative. For other materials (including carbon steels, alloyed steels, graphite, etc.), α is a nonlinear function of temperature owing to a nonlinear function of ρ versus temperature. As an example, Figure 3.4 shows ρ of one of the commercial graphite grades as a function of temperature. At melting point, the electrical resistivity of metals is typically sharply increased (Figure 3.5).
Figure 3.4 Electrical resistivity of graphite (ATL) versus temperature.
Figure 3.5 There is noticeable increase in electrical resistivity of metals near the melting point.
Impurities observed in metals and alloying distort the metal lattice and can affect the behavior of ρ to a considerable extent. As an illustration, Figure 3.6 shows an effect of the most common alloying and residual elements on electrical resistivity of iron [46].
Figure 3.6 Electrical resistivity versus percentage of alloying element in iron. (Recreated from R. Bozorth, Ferromagnetism, IEEE Press, 1993.)
For some binary alloys, the behavior of ρ versus the concentration of alloying elements is represented by a bellshaped curve. This curve often has maximum electrical resistivity at the concentration of alloying elements equal to 50% of the weight [47]. Figure 3.7a illustrates this phenomenon for Cu–Ni alloys.
Figure 3.7 (a) Electrical resistivity of Cu–Ni binary alloys at two different temperatures. (Recreated from K. Schroder, CRC Handbook of Electrical Resistivities of Binary Metallic Alloys, CRC Press, Inc., 1983.) (b) Effect of graphite shape on the electrical and thermal conductivities of gray and ductile iron relative to those of steels. (Recreated from C. Walton, T. Opar, Iron Castings Handbook, Iron Castings Society, Inc., 1981.)
When IH alloys, it is important to have a clear understanding regarding variations of the physical properties, including electrical resistivity. Incorrect assumption using the average values of ρ can result in costly misleading postulation.
In some cases, instead of having the bellshaped curve, ρ continuously decreases or increases with the concentration of alloys. For example, the electrical resistivity of plain carbon steels increases with an increase in the carbon content. For powder metallurgy materials, ρ decreases with an increase of density.
The value of electrical resistivity is also affected by the grain size (e.g., higher ρ typically corresponds to finer grains), plastic deformation, heat treatment, and some other factors, but to a smaller extent compared to the effect of temperature and chemical composition.
One should not confuse electrical resistivity ρ (Ω*m) with electrical resistance R (Ω), which is a property of an electrical circuit. The relationship between these parameters can be expressed as
where l is the length of the currentcarrying conductor and a is the area of the conductor’s cross section where the current is flowing through.
Electrical resistivity affects practically all important parameters of an induction system including depth of heat generation, temperature distribution, heating efficiency, coil impedance, and others. An effect of ρ on a particular parameter of the induction system is discussed further in the text in the appropriate sections.
Relative magnetic permeability μ_{r} indicates the ability of a material (e.g., metal) to conduct the magnetic flux better than a vacuum or air. Relative permittivity (or dielectric constant) ε indicates the ability of a material to conduct the electric field better than a vacuum or air. Both μ_{r} and ε are nondimensional parameters and have very similar meanings.
Relative magnetic permeability has a marked effect on process parameters selection affecting electrical phenomena, including the skin effect, electromagnetic edge, and end effect, as well as proximity and ring effects. Relative permittivity does not usually have a measurable impact when IH metallic materials, but it plays a major role in dielectric heating applications.
The constant μ_{0} = 4π × 10^{−7} H/m [or Wb/(A.m)] is called the permeability of free space (the vacuum), and similarly the constant ε_{0} = 8.854 × 10^{−12} F/m is called the permittivity of free space.
The product of μ_{r} and μ_{0} is called magnetic permeability μ and corresponds to the ratio of the magnetic flux density (B) to magnetic field intensity (H).
Basic definitions, interrelations between these two properties of the magnetic field (B and H), and its designations are discussed in high school science textbooks and college physics textbooks and here we simply refer to those texts.
In everyday engineering language, the IH professionals often call the relative magnetic permeability simply magnetic permeability. The relative magnetic permeability is closely related to magnetic susceptibility (χ) by the expression
In other words, magnetic susceptibility shows the amount by which μ_{r} differs from unity.
All materials based on their magnetization ability can be divided into paramagnetic, diamagnetic, and ferromagnetic. Relative magnetic permeability of paramagnetic materials is slightly greater than 1 (μ_{r} > 1). The value of μ_{r} for diamagnetic materials is slightly less than 1 (μ_{r} < 1). Because of insignificant differences of μ_{r} for both paramagnetic and diamagnetic materials, in IH practice, those materials are simply referred to as nonmagnetic materials (e.g., aluminum, copper, titanium, and tungsten).
In contrast to paramagnetic and diamagnetic materials, ferromagnetic materials exhibit the high value of relative magnetic permeability (μ_{r} >> 1). There are only a few elements that reveal the ferromagnetic properties at room temperature. These include iron, cobalt, and nickel. Some rare earth metals are ferromagnetic at temperatures below room temperature. All plain carbon steels are ferromagnetic. There are also a large number of alloy steels that belong to the group of ferromagnetic materials.
The ferromagnetic property of the material is a complex function of structure, chemical composition, prior treatment, grain size, frequency, magnetic field intensity, and temperature. As one can see from Figure 3.8a, the same carbon steel grade at the same temperature and frequency can have a noticeably different value of μ_{r} owing to differences in the intensity of the magnetic field (coil power). For example, the μ_{r} of magnetic steels commonly used in IH can vary from small values (e.g., μ_{r} = 2 or 3) to very high values (e.g., more than 500), depending on the magnetic field intensity H and temperature.
Figure 3.8 Effect of temperature and ﬁeld intensity on relative magnetic permeability μr (a) and the Curie temperature of plain carbon steel versus carbon content (b).
The temperature at which a ferromagnetic body loses its magnetic properties becoming nonmagnetic is called the Curie temperature (Curie point) and also often referred to as A _{2} critical temperature. Figure 3.8b shows a portion of the Fe–Fe_{3}C phase transformation diagram that illustrates the A _{2} critical temperature being a function of carbon content for plain carbon steels. Table 3.3 shows the Curie temperatures of some ferromagnetic materials.
Magnetic Material 
1008 
1060 
Permalloy 
Cobalt 
Nickel 
Temperature, °C (°F) 
768 (1414) 
732 (1350) 
440 (824) 
1120 (2048) 
358 (676) 
Thus, even among the plain carbon steels, the Curie temperature might be different because of the variation of carbon content. For example, as one can see from Table 3.3, the Curie temperature of plain carbon steel SAE 1008 is noticeably different from steel 1060 (768°C/1414°F vs. 732°C/1350°F).
The maximum value of relative magnetic permeability
is greatly affected by the chemical composition and microstructure. For example, the of highcarbon steel with 1.2% carbon is more than three times lower compared to the of lowcarbon steel with 0.1% carbon content.The magnetization curve describes the nonlinear relationship between magnetic flux density B and magnetic field intensity H. The nonlinear variation of μ_{r} = B/(Hμ_{0}) for a typical carbon steel is shown in Figure 3.9a. The maximum μ_{r} occurs at the “knee” of the curve. The magnetic field intensity H _{cr} that corresponds to the maximum permeability is called a critical value of H. When H > H _{cr}, the magnetic permeability will decrease with increasing H. If H → ∞, then μ_{r} → 1. In conventional IH of metals, the magnetic field intensity at the workpiece surface H _{surf} is typically much greater that H _{cr}.
Figure 3.9 (a) Magnetic field intensity (B) and μr versus field intensity (H) and (b) distribution of H and μr along the radius of a homogeneous carbon steel cylinder.
Similar to the current distribution, the magnetic field intensity is at its maximum value at the surface of the homogeneous workpiece and falls off exponentially toward the workpiece’s core (Figure 3.9b). As a result, the μ_{r} varies within the magnetic body. At the surface,
corresponds to the surface magnetic field intensity H _{surf}. In quick estimations under the assumption of using an electromagnetically long solenoid coil, H _{surf} can be considered as the field intensity in the air gap between the coil and the workpiece. With increasing distance from the surface, μ_{r} increases, and after reaching its maximum value at H = H _{cr}, it begins to fall off (Figure 3.9b). The complex nature of μ_{r} as a complex function of temperature and magnetic field intensity is shown in Figure 3.10a.Figure 3.10 (a) Relative magnetic permeability as a function of field intensity and temperature [1]. (b) Comparison of magnetic permeability in relatively “weak” and in “strong” magnetic fields.
From Figures 3.8 and 3.10a, one might conclude that μ_{r} always decreases with temperature. In the majority of IH and heattreating applications, it is the case. However, in a relatively “weak” magnetic field, μ_{r} might first increase with temperature and only near the Curie point would magnetic permeability start to sharply decline (Figure 3.10b).
As one may know from the basics of electricity, when a direct current flows through a conductor that stands alone (e.g., bus bar or cable), the electrical current distribution within the conductor’s cross section is uniform. However, when an AC flows through the same conductor, the current distribution is not uniform. The maximum value of the current density will always be located on the surface of the conductor with homogeneous electromagnetic physical properties; the current density will decrease from the surface of the conductor toward its center. This phenomenon of nonuniform current distribution within the conductor cross section is called the skin effect, which always occurs (though to a different degree) whenever there is AC flow. Therefore, the skin effect will be present in a workpiece located inside an induction coil. Figure 3.11 shows an appearance of skin effect in induction billet heating.
Figure 3.11 Appearance of the skin effect in induction billet heating. (Courtesy of Inductoheat Inc., an Inductotherm Group company.)
This effect is one of the major factors that cause the concentration of eddy current in the surface layer (“skin”) of the workpiece. Because of the circumferential nature of the eddy current induced in the solid cylinder, there is no current flow at its center (Figure 3.12). As an example, Figure 3.13a illustrates a comparison of current density distribution along the radius of a nonmagnetic metallic solid cylinder at two different frequencies F _{1} and F _{2}.
Figure 3.12 Instantaneous distribution of AC currents in “coiltoworkpiece” system.
Figure 3.13 Eddy current density distribution along the radius of a stainless steel cylinder at two different frequencies (a) and current and power density profiles owing to skin effect (b).
Because of the skin effect, approximately 86% of the power will be concentrated in the surface layer of the conductor that is called the reference (or current penetration) depth δ. The degree of skin effect depends on the frequency, material properties (ρ and μ_{r}) of the conductor, and the geometry of the workpiece. There will be a pronounced skin effect when high frequency is applied or when the radius of the workpiece is relatively large compared to δ (Figure 3.13b). The distribution of the current density along the workpiece thickness (assuming homogeneous electromagnetic properties) can be roughly calculated by the equation
where I is the current density at distance y from the surface (A/m^{2}), I _{0} is the current density at the workpiece surface (A/m^{2}), y is the distance from the surface toward the core (m), and δ is penetration depth (m).
Penetration depth is described in meters as
where F = frequency (Hz [cycles/s]), ρ = electrical resistivity of the electrically conductive material (Ω*m), and μ_{r} = relative magnetic permeability, or in inches as
where electrical resistivity ρ is in Ω*in.
As one can see, the value of δ varies with the square root of electrical resistivity and inversely with the square root of frequency and relative magnetic permeability. Mathematically speaking, the penetration depth in Equation 3.5 is the distance from the surface of the electrical conductor toward its core at which the current decreases exponentially to “1/exp” its value at the surface. The power density at this distance will decrease to “1/exp^{2}” its value at the surface. Figure 3.13b illustrates a skin effect appearance by showing distribution of current density and power density from the workpiece surface toward the core. As one can see from the figure, at one penetration depth from the surface (y = δ), the current density will equal approximately 37% of its surface value. However, the power density will equal 14% of its surface value. From this, we can conclude that 63% of the current and 86% of the power will be concentrated within a surface layer of thickness δ.
During the heating cycle, the ρ of most metals can increase to four to six times its initial value. Therefore, even for nonmagnetic metals, during the heating cycle, δ can increase substantially. Table 3.4 shows the δ of some commonly used nonmagnetic metals. Table 3.5 shows the δ of plain carbon steel SAE 1040 at room temperature (21°C or 70°F).
Metal 
T 
ρ 
Frequency (kHz) 


°C 
°F 
μΩ·m 
μΩ·in. 
0.06 
0.50 
1 
2.5 
4 
8 
10 
30 
70 
200 
500 

Aluminum 
20 
68 
0.027 
1.06 
10.7 
3.70 
2.61 
1.65 
1.30 
0.92 
0.83 
0.48 
0.31 
0.18 
0.12 
250 
482 
0.053 
2.09 
15.0 
5.18 
3.66 
2.32 
1.83 
1.29 
1.16 
0.67 
0.44 
0.26 
0.16 

500 
932 
0.087 
3.43 
19.2 
6.64 
4.69 
2.97 
2.35 
1.66 
1.48 
0.86 
0.56 
0.33 
0.21 

Copper 
20 
68 
0.018 
0.71 
8.81 
3.05 
2.16 
1.36 
1.08 
0.76 
0.68 
0.39 
0.26 
0.15 
0.10 
500 
932 
0.050 
1.97 
14.5 
5.03 
3.56 
2.25 
1.78 
1.26 
1.12 
0.65 
0.43 
0.25 
0.16 

900 
1652 
0.085 
3.35 
19.3 
6.67 
4.72 
2.98 
2.36 
1.67 
1.49 
0.86 
0.56 
0.33 
0.21 

Brass 
20 
68 
0.065 
2.56 
16.6 
5.74 
4.06 
2.56 
2.03 
1.43 
1.28 
0.74 
0.48 
0.29 
0.18 
400 
752 
0.114 
4.49 
21.9 
7.60 
5.37 
3.40 
2.69 
1.90 
1.70 
0.98 
0.64 
0.38 
0.24 

900 
1632 
0.203 
7.99 
29.3 
10.1 
7.17 
4.53 
3.58 
2.53 
2.27 
1.31 
0.86 
0.51 
0.32 

Stainless steel 
20 
68 
0.690 
27.2 
53.9 
18.7 
13.2 
8.36 
6.61 
4.67 
4.18 
2.41 
1.58 
0.93 
0.59 
800 
1472 
1.150 
45.3 
69.6 
24.1 
17.1 
10.8 
8.53 
6.03 
5.39 
3.11 
2.04 
1.21 
0.76 

1200 
2192 
1.240 
48.8 
72.3 
25.1 
17.7 
11.2 
8.86 
6.26 
5.60 
3.23 
2.12 
1.25 
0.79 

Silver 
20 
68 
0.017 
0.67 
8.34 
2.89 
2.04 
1.29 
1.02 
0.72 
0.65 
0.37 
0.24 
0.14 
0.09 
300 
572 
0.038 
1.50 
12.7 
4.39 
3.10 
1.96 
1.55 
1.10 
0.98 
0.57 
0.37 
0.22 
0.14 

800 
1472 
0.070 
2.76 
17.2 
5.95 
4.21 
2.66 
2.10 
1.49 
1.33 
0.77 
0.50 
0.30 
0.19 

Tungsten 
20 
68 
0.050 
1.97 
14.5 
5.03 
3.56 
2.25 
1.78 
1.26 
1.12 
0.65 
0.43 
0.83 
0.53 
1500 
2732 
0.550 
21.7 
48.2 
16.7 
11.8 
7.46 
5.90 
4.17 
3.73 
2.15 
1.41 
0.83 
0.53 

2800 
5072 
1.040 
40.9 
66.2 
22.9 
16.2 
10.3 
8.11 
5.74 
5.13 
2.96 
1.94 
1.15 
0.73 

Titanium 
20 
68 
0.500 
19.7 
45.9 
15.9 
11.3 
7.11 
5.62 
3.98 
3.56 
2.05 
1.34 
0.80 
0.50 
600 
1112 
1.400 
55.1 
76.8 
26.6 
18.8 
11.9 
9.41 
6.65 
5.95 
3.44 
2.25 
1.33 
0.84 

1200 
2192 
1.800 
70.9 
87.1 
30.2 
21.3 
13.5 
10.7 
7.54 
6.75 
3.90 
2.55 
1.51 
0.95 
Frequency (Hz) 


Magnetic Field Intensity 
60 
500 
3000 
10,000 
30,000 
100,000 

Penetration Depth 

A/mm 
A/in. 
mm 
in. 
mm 
in. 
mm 
in. 
mm 
in. 
mm 
in. 
mm 
in. 
10 
250 
2.50 
0.100 
0.88 
0.034 
0.36 
0.014 
0.2 
0.008 
0.11 
0.004 
0.06 
0.002 
40 
1000 
4.70 
0.185 
1.63 
0.064 
0.67 
0.026 
0.36 
0.014 
0.21 
0.008 
0.12 
0.005 
80 
2000 
6.30 
0.249 
2.20 
0.086 
0.9 
0.035 
0.49 
0.019 
0.28 
0.011 
0.16 
0.006 
160 
4050 
8.76 
0.345 
3.03 
0.119 
1.24 
0.049 
0.68 
0.027 
0.39 
0.015 
0.21 
0.008 
200 
5100 
9.63 
0.379 
3.33 
0.131 
1.36 
0.054 
0.75 
0.029 
0.43 
0.017 
0.24 
0.009 
280 
7100 
11.20 
0.442 
3.89 
0.153 
1.59 
0.062 
0.87 
0.034 
0.50 
0.020 
0.27 
0.011 
While discussing the threedimensional (3D) behavior of μ_{r} within a ferromagnetic body, it is necessary to mention that the definition of δ in its classical form (Equations 3.6 and 3.7) does not have a fully determined meaning because of the nonconstant distribution of μ_{r} within the workpiece even at ambient temperature. In engineering practice, the value of relative magnetic permeability at the surface of the workpiece
is typically used to define those equations.As stated earlier, δ is a function of temperature. At the beginning of the heating cycle, the current penetration into the carbon steel will increase slightly (Figure 3.14a) because of the increase in ρ with temperature. With a further rise of temperature (at approximately 550°C or 1022°F), μ_{r} starts to markedly decrease. Near a critical temperature A _{2}, known as the Curie temperature or Curie point, μ_{r} drastically drops to unity because the carbon steel becomes nonmagnetic. As a result, δ increases significantly. After heating above the Curie temperature, δ will continue to rise slightly because of the increase in ρ with temperature (Figure 3.14a).
Figure 3.14 Typical variation of current penetration depth (δ) during IH of a carbon steel workpiece (a) and variation of δ for different metals during IH (b).
The variation of δ during IH of a carbon steel workpiece drastically changes the degree of skin effect. Figure 3.14b shows how many times the δ of some metals can increase during a typical heating cycle. It is especially important to take this phenomenon into consideration for through heating applications.
In most publications devoted to IH, distributions of current and power densities (heat source distribution) along the workpiece thickness/radius are simplified and introduced as being exponentially decreasing from the surface into the workpiece. However, this assumption is correct only for a homogeneous nonmagnetic solid body with constant ρ.
Realistically speaking, this assumption can be made for only some unique cases because for the great majority of IH applications (owing to the skin effect and some other electromagnetic phenomena that are discussed in the following sections), the current density distribution is not uniform and there are always thermal gradients within the heated workpiece. These thermal gradients result in nonuniform distribution of ρ and μ_{r} within the workpiece.
In addition, as shown in the previous section, μ_{r} is nonuniform along the thickness of the ferromagnetic workpiece owing to a nonuniform distribution of the magnetic field intensity (Figure 3.9b). Therefore, an assumption of exponential current density distribution (Equation 3.5 ) can be used for the purpose of rough engineering estimations for heating of nonmagnetic materials only.
In some applications, such as surface hardening, the power density distribution along the radius/thickness has a unique wavelike form, which is significantly different from the commonly assumed exponential distribution. The maximum power density is located at the surface. Then, the power density decreases toward the core. However, once it reaches a certain distance from the surface, the power density starts to increase again, and after reaching a maximum, it starts its final decline.
M. Lozinskii and P. Simpson originally independently suggested this waveshaped distribution of current and power densities in their respective texts [7,8]. Both authors intuitively felt that there should be situations when distribution of the power density (heat sources) would be different compared to the commonly assumed exponential form. Unfortunately, because of computer modeling limitations at that time and the unavailability of sufficiently sophisticated software to accurately model the IH processes, it was not possible to reveal that phenomenon in more detail. Obviously, it was not possible to measure the power density distribution within the solid body during its heating without disturbing the current flow.
The new generation of computer modeling software allows one to make a quite accurate prediction regarding the physics of that power density waveshaped phenomenon revealing the true nature of the skin effect [1,2,17–21,44,52–54]. An impact of this phenomenon in different applications (surface hardening, bar heating, etc.) is discussed later in this text in corresponding sections.
When discussing the skin effect, it is appropriate to introduce the terms electromagnetically thick and electromagnetically thin bodies [1]. Electromagnetically speaking, a workpiece can “act” as thick or thin bodies depending on the chosen frequency and magnetic field orientation. If δ is greater compared to the thickness or diameter of the solid body, then it becomes semitransparent to the electromagnetic field and could be considered an electromagnetically thin body. There is a distinct current cancellation within the electromagnetically thin bodies absorbing only a negligible amount of energy (transverse flux inductors and traveling wave inductors are exceptions). Being practically transparent to the external electromagnetic field, there will be only small amounts of the heat generation appearing in electromagnetically thin bodies.
If the thickness or diameter of the solid electrically conductive body is six times the δ, then it can be considered as an electromagnetically thick body. Since δ can increase more than 15 times during the heating cycle, the workpiece that initially could be considered as an electromagnetically thick body can become at the end of the heating cycle an electromagnetically thin body that is accompanied by a dramatic reduction in heating efficiency.
Besides the frequency, the orientation of the workpiece with respect to the electromagnetic field may have a marked effect on considering the workpiece as an electromagnetically thin or thick body. If orientation of the solid workpiece compared to the inductor results in eddy current cancellation, then the workpiece can be considered a thin body; otherwise, it can be called an electromagnetically thick body. Figure 3.15 shows that depending on plate orientation compared to the inductor, it can act electromagnetically as a thin or thick body.
Figure 3.15 Illustration of the concept of electromagnetically thin and thick bodies.
Geometry alone cannot determine whether a workpiece should be considered an electromagnetically thin or thick body. For example, when using a solenoid coil, a stainless steel billet of 25 mm (1 in.) diameter will act during IH as an electromagnetically thin body using frequencies lower than 500 Hz. In contrast, a stainless steel billet of 12.7 mm (0.5 in.) diameter will be considered an electromagnetically thick body if a frequency of 70 kHz is applied.
When discussing the skin effect in conductors or cables, it was assumed that a conductor stands alone and that there are no other currentcarrying conductors in the surrounding area. In most practical applications, this is not the case. Most often, there are other conductors in proximity. These conductors have their own magnetic fields, which interact with nearby fields, affecting the current flow and power density distributions.
An analysis of the effect on current distribution in a conductor when another conductor is placed nearby is given below. Figure 3.16 shows the magnetic field (a) and current density distributions (b) in a rectangular conductor (e.g., bus bar) that stands alone. The appearance of skin effect in a rectangular body is clearly visible on Figure 3.16b.
Figure 3.16 Magnetic field (a) and current density (b) distributions in a standalone rectangular conductor.
When another currentcarrying conductor is placed near the first one, the currents in both conductors will redistribute. If the currents flowing in the conductors have opposite directions, then both currents (Figure 3.17b) will be concentrated in the areas facing each other (internal areas). However, if the currents have the same direction, then these currents will be concentrated on opposite sides of the conductors (Figure 3.18b).
Figure 3.17 Magnetic field (a) and current density (b) distributions in bus bars owing to the proximity effect when the currents are flowing in the opposite directions.
Figure 3.18 Magnetic field (a) and current density (b) distributions in bus bars owing to the proximity effect when the currents are flowing in the same direction.
A strong magnetic field forms in the area between the bus bars when the currents flow in opposite directions (Figure 3.17a). This occurs because in this area, the imaginary lines of magnetic field are produced by each bus bar having the same direction complementing each other. Since the current is concentrated in the internal areas, the external magnetic field will be noticeably weaker. In this case, the external magnetic fields generated by each conductor have opposite directions and tend to cancel each other. This phenomenon is used in coaxial cables.
The opposite is true: if the currents have the same direction, then the magnetic field lines will have opposite directions in the area between bus bars canceling each other (Figure 3.18a). Because of this cancellation, a relatively weak magnetic field will occur between the bus bars; however, the external magnetic field will be quite strong because the imaginary lines of magnetic field produced by the two conductors will have the same direction complementing each other.
If the distance between conductors increases, then the strength of the electromagnetic proximity effect will decrease. Proximity effect in the case of nonsymmetrical systems is shown in Figure 3.19.
Figure 3.19 Proximity effect in nonsymmetrical systems.
The phenomenon of proximity effect is directly related to IH. Induction systems consist of at least two conductors [1,6,7,44,54,55]. One of these conductors is an inductor that carries the source current (Figure 3.12), and the other is the electrically conductive workpiece located near the inductor. Because of Faraday’s law, eddy currents induced within the workpiece have an opposite direction compared to that of the source current of the inductor. Therefore, as a result of the proximity effect, the coil current and eddy currents induced in the workpiece will concentrate in the areas facing each other. This is the second factor that causes a current redistribution in an IH system as shown in Figure 3.12.
Figure 3.20 shows how the electromagnetic proximity effect produces different heating patterns. A carbon steel cylinder is located asymmetrically inside a singleturn inductor. Two noticeably different patterns will be developed if the cylinder is statically heated (without rotation). The appearance of these patterns is caused by a difference in the eddy current flow.
Figure 3.20 Proximity effect in singleturn coil with nonsymmetrical positioning of the workpiece.
As shown in Figure 3.20a, the eddy currents have a higher density in the area where the coiltoworkpiece air gap is small (good coupling), resulting in an intense heating there. The heat pattern will be relatively narrow and deep. In the area with the larger air gap (poor electromagnetic coupling), the temperature rise will not be as significant as in the case of good coupling. Also, the heat pattern will be noticeably wider and shallower (Figure 3.20b).
Depending on application specifics, the presence of an electromagnetic proximity effect might be harmful or beneficial. Inappropriate combination of inductor design, process recipe, and workpiece handling mechanism could result in the undesirable appearance of the proximity effect, which can lead to localized cold and hot spots, overheating, and even melting. Figure 3.21 shows three examples of nonuniform heating that occurred because of an inappropriate appearance of the proximity effect. Figure 3.21a,b shows the result of localized severe overheating of copper (a) and aluminum (b) thin wall tubing caused by a proximity effect owing to inadequate tube handling during its processing through an induction coil. Figure 3.21c shows localized “hot” spots at edge areas of a deformed tubular workpiece during scan heating.
Figure 3.21 Results of localized severe overheating of thin wall copper (a) and aluminum (b) tubing caused by a proximity effect owing to improper tube handling during its processing through an induction coil and (c) localized “hot” spots at edge areas when heating a deformed tubular workpiece.
On the other hand, an electromagnetic proximity effect is a major beneficial factor in heating selective areas (e.g., selective hardening, localized tempering, and stress relieving). Profiled coils and singleshot inductors (Figure 3.2) primarily rely on a proximity effect as a means of obtaining required heating patterns.
An understanding of the physics of the electromagnetic proximity and skin effects is important not only in IH but also in power supply design and bus network design. The proper design of a bus network will significantly decrease its impedance, minimizing voltage drop and reducing transmission power losses, and thus improving overall energy efficiency.
When we discussed the proximity effect, we first reviewed the magnetic field and current density distributions in a standalone currentcarrying conductor (Figure 3.16). The magnetic field and current density will redistribute when an electrically conductive body (e.g., workpiece) is placed in its proximity (Figure 3.22). Magnetic field is squeezed in the gap between the inductor and the workpiece, resulting in the highest magnetic flux density there (Figure 3.22a). As shown in Figure 3.22b, an appreciable portion of the inductor’s current will flow near the surface of the conductor that faces the workpiece. The remainder of the current will be distributed on the sides of the conductor; however, a small portion of the current will occur on its opposite side [1,6,7,44,51].
Figure 3.22 Magnetic field (a) and current density (b) redistribution when an electrically conductive body (e.g., workpiece) is placed in proximity to a rectangular conductor.
Continuing our study, let us position an external magnetic flux concentrator (e.g., Cshaped laminations) around this conductor, as shown in Figure 3.23. As a result, the magnetic concentrator provides a low reluctance path for a magnetic field as shown in Figure 3.23a.
Figure 3.23 Magnetic flux concentrator provides a low reluctance path for a magnetic field (a) and squeezes the inductor current (b) to the “open surface” of the concentrator.
Practically all of the conductor’s current will be concentrated on the surface facing the workpiece (Figure 3.23b). In other words, the magnetic concentrator squeezes the inductor current to the “open surface” of the concentrator—to the open area of the slot. This phenomenon is called an electromagnetic slot effect. The electromagnetic slot effect is widely used in the heating of selective areas and helping to dramatically reduce external magnetic field. Thanks to this phenomenon, there will be improved inductortoworkpiece electromagnetic coupling that improves the electrical efficiency of IH.
It is necessary to mention here that the slot effect will also take place without the workpiece (Figure 3.24). In this case, the current will be slightly redistributed in the conductor, but most of it will still be concentrated in the “open surface” area of the slot. The actual current distribution in the conductor depends on the frequency, magnetic field intensity, geometry, and electromagnetic properties of the conductor and the concentrator.
Figure 3.24 Slot effect takes place without the presence of the workpiece in the standalone conductor.
The electromagnetic slot effect and the proximity effect play a particularly important role in channel, hairpin, singleshot, and splitreturn inductors, as well as inductors for heating internal diameters.
The slot effect is widely used not only in connection with IH but also in the design of other industrial machines such as motor generators, AC and DC machines, and in cases when it is required to have an electromagnetic shielding of certain electrically conductive components (e.g., cabinets, fixtures, or electronic devices).
Up to now, we have discussed current density distribution in straight conductors. One such conductor, a rectangular bus bar, and its current distribution are shown in Figure 3.25. If that currentcarrying bar is bent to shape it into a ring, then its current will be redistributed. Magnetic flux lines will be concentrated inside the ring, increasing magnetic flux density there. Outside the ring, the magnetic flux lines will be disbursed. As a result, most of the current will flow within the thin inside surface layer of the ring where there will be the shortest distance and the lowest impedance path [1,2,6,7,54]. As one can see, this ring effect is also somewhat similar to the proximity effect because currents flowing on inside surfaces of the opposite sides of the ring’s circumference are oriented in opposite directions (thus being attracted to each other).
Figure 3.25 Ring effect in a rectangular conductor.
Figure 3.26 shows the appearance of the electromagnetic ring effect in cylinder conductors leading to a concentration of current on the inside surface of the induction coil. The ring effect takes place not only in singleturn inductors but also in multiturn coils. Therefore, it is the third electromagnetic effect that is responsible for the current distribution in the induction system shown in Figure 3.12.
Figure 3.26 Ring effect in a round conductor.
The presence of the ring effect can have a positive or negative impact on the heating and process efficiency. For example, in conventional IH when the solid cylinder workpiece is located inside the solenoid induction coil, this effect plays a positive role because in combination with the skin and proximity effects, it will lead to a concentration of the coil current on the inside diameter of the coil. As a result, there will be improved coiltoworkpiece electromagnetic coupling, which leads to a coil efficiency increase.
The ring effect plays a negative role in the IH of internal surfaces or inside diameters (socalled I.D. heating), where the inductor is located inside the hollow workpiece (Figure 3.27a). In this case, the ring effect leads to a coil current concentration on the inside diameter of the coil. This worsens an equivalent coiltoworkpiece electromagnetic coupling and, therefore, decreases coil efficiency [56]. However, despite the ring effect, the proximity effect here tends to shift the coil current to an outside surface of the coil, which is highly desirable for increasing process efficiency. Therefore, the coil current distribution in such applications is a result of two counteracting electromagnetic phenomena: the proximity and ring effects. For moderate and, especially, for relatively small inside diameters, the electromagnetic ring effect appearing in the coil usually overpowers the proximity effect and forces majority of the coil current to flow on its inside diameter. Thus, this results in decoupling of the coil current from the workpiece, dramatically reducing coil electrical efficiency and drastically increasing coil copper kilowatt losses and energy consumption.
Figure 3.27 Appearance of the ring effect when heating inside diameters. Compare the difference in magnetic field distribution when using a bare coil (a) versus a coil with a “U”shaped magnetic flux concentrator (b).
In order to “help” the proximity effect dominate the ring effect, in a great majority of I.D. heating applications, a magnetic flux concentrator is located inside the coil. This allows a slot effect to “assist” the proximity effect to improve equivalent electromagnetic coupling, measurably increase the coil efficiency, and dominate (but not completely eliminate) the ring effect (Figure 3.27b).
The ring effect should be taken into consideration in power supply design. Because of this effect, the current is concentrated in areas where buses are bent, potentially leading to undesirable overheating of localized areas of the buses and appearance of “hot” spots. To avoid local overheating, it is necessary to take this phenomenon into account when designing the cooling circuits.
Electromagnetic forces play an important role in many modern technologies including motors, magnetohydrodynamic seals, electromagnetic pumps, levitators, electrical bearings, and springs. In some applications, electromagnetic forces can reach tremendous values. For example, thanks to a capability to develop incredibly large electromagnetic forces, electric guns or launchers can fire materials to higher velocities than are achievable by rockets or chemical/powder guns [1,44,57].
In many IH applications, coil currents can also reach substantial values. For example, currents of 3–5 kA and higher are not unusual for some shaft and gear hardening. High currents produce significant forces that have a pronounced effect on coil life and system design. Without proper consideration, those forces can reallocate the heated workpiece, fixture, or flux concentrator, and even bend the induction coil turns or fixture, which may negatively affect overall system reliability, repeatability, and heating quality [1,44,57].
Unfortunately, electromagnetic forces are not so often discussed in IH publications. Many of the seemingly endless variety of workpieces heated by induction require specific coil geometry (Figure 3.2), which makes it difficult to develop general but simple procedures to evaluate electromagnetic forces. This paragraph is intended to at least partially remedy this by providing a brief introduction to the topic.
A currentcarrying conductor placed in a magnetic field experiences a force that is proportional to current and magnetic flux density. Thanks to an experimental study by Ampère and BiotSavart [1,44,58–65], this force can be quantified. If a currentcarrying element dl, carrying a current I, is placed in an external magnetic field B, it will experience a force dF according to Equation 3.8:
where F, I, and B are vectors and φ is the angle between the direction of the current I and magnetic flux density B.
In SI units, the force is measured in Newton (1 N = 0.102 kgf = 0.225 lbf). Figure 3.28 shows that the direction of the force experienced by the element dl of the currentcarrying conductor placed in an external magnetic field B can be determined based on the lefthand rule (FBI rule). According to the rule, if the middle finger of the left hand follows the direction of current flow and the pointer finger follows the direction of the magnetic flux of the external field (imaginary magnetic field lines head into the palm), then the thumb will show the direction of the force.
Figure 3.28 Lefthand (FBI) rule of magnetic force.
It is important to remember from Equation 3.8 that if the angle φ between the direction of the current I and magnetic field B is equal to zero, then sin φ = 0 and, therefore, no force will be experienced by the currentcarrying conductor. In other words, if the currentcarrying conductor is parallel to an external magnetic field, then it will not experience any forces from that field. The lesson here is that among other factors, magnetic force depends greatly on orientation of the currentcarrying conductors.
Let’s consider some of the most common cases of the appearance of magnetic forces in IH applications.
Figure 3.29 Magnetic force in currentcarrying conductors.
Figure 3.30 Magnetic interaction between two thin wires.
Figure 3.31 Magnetic forces in an empty solenoid coil.
Figure 3.32 Magnetic forces in bar end heating of magnetic (b) and nonmagnetic (a) bars.
In most IH applications, the electromagnetic force has a complex 3D distribution. Depending on coil design and process specifics, one of three force components—longitudinal, radial, or hoop—may be significantly greater than the others. It is important to remember that the orientation and 3D distribution of forces during the heating cycle is not a function of only the geometry of the system, and is not constant. During the heating cycle, the force distribution also depends on frequency, power density, temperature/material properties, heating mode (constant power, current, voltage), and other parameters.
Excessive magnetic forces can be harmful to the rigidity of the induction system, causing intensive vibration and industrial noise. However, in other cases, those forces can be desirable and play an important technological part of the process (e.g., electromagnetic stirring in melting furnaces, electromagnetic pumps and locks, levitators, induction plasma, electromagnetic separators, etc.).
Bear in mind that the formulas given here can be applied only in some specific/simplified cases. For the majority of IH applications having complex geometries, numerical computer modeling is required to help the designer accurately evaluate the electromagnetic forces and to determine which actions should be taken to develop robust and reliable coil/fixture design.
Case Study. An example given in Figure 3.33 shows the magnetic field distribution (right images) around a twoturn induction coil for hardening a 2in.diameter (50mmdiameter) shaft using a power/frequency combination of 180 kW/1 kHz [57].
Figure 3.33 Computer modeling plots of the magnetic field distribution and major forces experienced by a twoturn induction coil when heating a steel shaft. Compare a magnitude of axial forces when shaft is heated below the Curie point (a) versus above it (b). Coil power and applied frequency were 180 kW and 1 kHz, respectively. (From V. Rudnev, Electromagnetic forces in IH, Professor Induction Series, Heat Treating Progress, ASM Int’l, Materials Park, Ohio, July, 25–28, 2005.)
Table 3.6 lists the magnetic forces produced when this carbon steel shaft is heated to temperatures below and above the Curie point.
Force Component 
Below Curie Point 
Above Curie Point 

F _{hoop} 
81 (18.2) 
298 (67) 
F _{longitudinal} 
443 (99.6) 
928 (209) 
F _{radial} 
1.4 (0.3) 
5.2 (1.2) 
Source: V. Rudnev, Electromagnetic forces in induction heating, Professor Induction Series, Heat Treating Progress, ASM Int’l, Materials Park, Ohio, July, 25–28, 2005.
When heating above the Curie point, each turn of the twoturn coil will experience the maximum longitudinal force of 928 N (209 lbf). Obviously, such intensive force cannot be neglected since it dramatically affects coil life and should be properly taken into consideration when designing IH systems. Thus, special effort should be made to provide a sufficient support to coil turns properly addressing magnetic forces.
Conclusion. Depending on the application specifics, if not addressed properly, magnetic forces can reach appreciable magnitudes adversely affecting the rigidity and repeatability of an IH system, causing premature coil failure owing to stress fatigue cracking and excessive vibration, and resulting in excessive acoustic noise. However, in other applications, those forces can be desirable, playing an important positive role in the process. Magnitude and orientation of electromagnetic forces must be addressed when designing IH systems and measures taken to reduce their magnitude.
As discussed in Section 3.1.2, the surfacetocore temperature difference is greatly affected by the skin effect. The temperature profiles along the workpiece’s length and width are affected by, among other factors, a distortion of the electromagnetic field in the coil’s end and edge regions. Those field distortions and corresponding distributions of induced currents are referred to as end and edge effects. These effects and the field distortion caused by them are primarily responsible for nonuniform temperature profiles in cylindrical, rectangular, and trapezoidal workpieces. Because of the great importance of these effects, much effort has been devoted to their study. The first attempt to provide a systematic analysis of electromagnetic end effects was carried out by D. Lavers in the late 1960s and early 1970s [75]. Further analysis of electromagnetic end effects and edge effects has been reported in Refs. [1,2,9,44,76–82] and others.
It is convenient to introduce these effects studying the IH of a rectangular slab. Suppose a slab is placed in an initially uniform magnetic field (Figure 3.34). If the slab’s length and width are much larger than its thickness, the electromagnetic field in the slab can be viewed as an area consisting of three zones: a central part, a transverse edge effect area, and a longitudinal end effect area. In the central part, the electromagnetic field distribution corresponds to the field in the infinite plate. Basically, end and edge effects have a 2D space distribution excluding only the zone of threeedge (3D) corners where the field is 3D and the corresponding heat source distribution is the result of a mixture of both the end and edge effects. For many practical applications, the separate study of end and edge effects is of great engineering interest.
Figure 3.34 (a) Sketch of coordinate system for slab, and electromagnetic end and edge effect of the slab (b). (From V. Rudnev, D. Loveless, Longitudinal flux IH of slabs, bars and strips is no longer “black magic,” Part 1, Industrial Heating, January, 29–34, 1995.)
It is convenient to study the end effect by estimating the power density distribution along the length of the cylindrical or rectangular workpiece (along the Zaxis in the Y–Z cross section, Figure 3.34).
The analysis of the edge effect is often conducted by evaluation of the power density distribution across the slab width (along the Xaxis in the Y–X cross section). The edge effect is typically negligible when heating cylindrical workpieces (e.g., billets or bars) in the longitudinal flux solenoid inductors when the axes of symmetry of the coil and the cylinder coincide. However, this effect can play an essential role when heating billets in ovaltype coils when billets are transferred transversally in respect to the coil.
Appropriate sections of Chapters 4 and 6 provide a detailed analysis of the main features of the electromagnetic end and edge effects in different applications. In this section, we provide a brief orientation regarding these phenomena.
Figure 3.35 shows the appearance of end effect when heating solid cylinders. As mentioned above, electromagnetic end effect represents distortion of the electromagnetic field in the extreme end of the cylinder. In the case of conventionally designed coils, the electromagnetic end effect at the extreme end (“hot” end) of the cylinder (Figure 3.35, “a”– region) is defined primarily by the following variables:
Figure 3.35 Sketch of induction heater (a) and normalized power density distribution (b) along the length of a solid cylinder.
where R is the radius of the heated cylinder, R _{i} is the inside coil radius, and δ is the current penetration depth. The effects of the frequency F and the electromagnetic physical properties of the heated material (ρ and μ_{r}) are included into the skin effect ratio R/δ. K _{space} (also referred to as coil pitch) indicates how tightly the turns are wound. K _{space} = (copper turn width)/(pitch of turns winding). For multiturn coils, its value is always less than 1. For a singleturn coil, K _{space} = 1. An incorrect combination of these factors can lead either to underheating or overheating the extreme (“hot”) end.
An incorrect combination of the abovementioned factors can lead to underheating or overheating the extreme end of the workpiece. As an illustration, Figure 3.36 shows two extreme cases of end effect appearance when heating statically hollow cylinders. In the case of heating a copper tube using 30 kHz and having a large coil overhang, an electromagnetic end effect manifests itself in severe overheating of the tube end (Figure 3.36a). In contrast, when heating an SAE 1018 carbon steel pipe using line frequency (60 Hz) and regardless of substantial coil overhang, the end region of the pipe is noticeably underheated compared to its regular region (Figure 3.36b).
Figure 3.36 Two extreme cases of end effect appearance when heating hollow cylinders. (a) An electromagnetic end effect manifests itself in severe overheating of the tube end when IH a copper tube using 30 kHz. (b) When heating AISI 1018 carbon steel pipe using line frequency (60 Hz) and regardless of substantial coil overhang, the end of the pipe is noticeably underheated.
As an example, Figure 3.37 shows the effect of the coil overhang σ on surface power density distribution along the length of the nonmagnetic cylinder bar when its trailing (left) end takes different positioning inside of a multiturn solenoid inductor and assuming pronounced skin effect and constant current. The presence of the large coil overhangs (σ_{5}, σ_{6}, σ_{7}) results in substantial power density surplus in the end region of the nonmagnetic cylinder compared to its central part. At the same time, there is a certain coil overhang (in the case shown in Figure 3.37, it approximately corresponds to σ_{4}) that results in a reasonably uniform power density distribution. In static heating, small overhangs and application of relatively low frequencies are associated with noticeable deficit of power density at the extreme end of the bar.
Figure 3.37 Illustration of the longitudinal end effect at different positions of the left end of a nonmagnetic cylinder during processing through a multiturn coil.
It is possible to show that in the case of an electromagnetically long multiturn coil with a homogeneous nonmagnetic cylinder, the density of the induced eddy current in the workpiece area under the coil opposite end (the socalled cold end) is only half that in the central part (Figure 3.35, zone “b”). This statement can be illustrated using the following simplified example.
The distribution of the magnetic field along the axis in the end area of the empty “ideal” solenoid coil can be relatively easily obtained using an expression that describes its distribution in a single loop of wire. The assumption of an “ideal” solenoid presumes the following conditions:
Figure 3.38 shows the sketch of such an “ideal” solenoid of length l and radius R that has N tightly wound turns. A currentcarrying empty loop produces a B field along the loop axis with “z” component according to the following expression [1,44,60–65,83]:
Figure 3.38 Sketch of a tightly wound multiturn solenoid for end effect calculation.
where Z is the axial distance from the loop to the area of interest, I is the loop current, and μ_{0} is the permeability of free space (a vacuum), μ_{0} = 4π × 10^{−7} H/m [or Wb/(A∗m)].
The magnetic field at the center of the empty loop can be obtained by assuming Z = 0 in Equation 3.13 resulting in
The magnetic field distribution along the axis of an empty solenoid can be obtained by expansion of B _{ z } of a single wire loop on a multiturn coil. Taking into consideration the assumption of tightly wound turns, the contribution of a small currentcarrying section “dZ” on the total field in the center of the solenoid will be as follows:
The total magnetic field in the center of the coil can be obtained by taking into account the contributions of all currentcarrying sections. Therefore, after integrating dB _{ z } along the coil length, the magnetic field along the coil center can be written as
After simple mathematical operations, the total axial field in the center of the solenoid will be
If the length of the solenoid is much greater than its radius l >> R (electromagnetically long coil), then it is possible to neglect R with respect to l and Equation 3.16 can be rewritten as
This is the wellknown expression for the axial B field at the center of an electromagnetically long solenoid. It is possible to show that by letting the corresponding limits in Equation 3.16, Equations 3.13 and 3.16 can be transformed into
Thus, for an electromagnetically long coil, Equation 3.19 can be approximated as
A comparison of Equations 3.18 and 3.20 shows that at the ends of the empty coil, the magnetic flux density B _{ z } drops to onehalf its value at the center. This conclusion can be roughly expanded for an electromagnetically long multiturn coil with an infinitely long homogeneous nonmagnetic workpiece.
As discussed above (see Figure 3.35, zone “b”), the density of the induced current under the coil end with sufficiently long nonmagnetic workpiece is two times less than in the middle of the coil. It means that the power density under the coil end is equal to a quarter of that in the center (P _{end} = 0.25∗P _{center}). Among other factors, the length of zone “b” depends on the skin effect in the heated workpiece, the ratio of “coil inside radius to workpiece radius,” and the coil turn space factor K _{space} and might be as long as four times the coil radius (if skin effect is not pronounced and there is poor coiltoworkpiece coupling) or 12 equivalent air gaps between the coil and the load (for pronounced skin effect and small air gaps).
The electromagnetic end effect in a magnetic slab has several features compared to the nonmagnetic one. Magnetic materials have a tendency to gather the magnetic flux lines thanks to μ_{r}. Generally speaking, the electromagnetic end effect in a ferromagnetic workpiece is mainly affected by two factors [1,9,44,76–81]:
The first factor causes an increase in power at the workpiece’s end (somewhat similar to the end effect of a nonmagnetic cylinder). The second factor causes a power reduction there. Therefore, unlike those of the nonmagnetic cylinder, the ends of the ferromagnetic cylinder, even having large coil overhangs, may be either overheated or underheated (being more typical).
The discussion above serves as a simplified introduction to a subject of electromagnetic end effect providing only a basic illustration of some critical factors that affect its appearance or how this effect can be controlled. An accurate estimation of this effect and its subtleties related to a particular application can be obtained after numerical computer modeling.
Case Study. As an example, Figures 3.39 and 3.40 show the results of FEA simulation of bar end heating using a solenoid multiturn coil with tighter wound turns and a “cold” end. The material is SAE 1039. The frequency is 1 kHz. The total length of the steel bar was 210 mm. After 45 s of heating, it was required to heat a 100mmlong section at the left end of the bar to a final temperature of 1160°C ± 35°C.
Figure 3.39 FEA mesh (top) and a coil field distribution when end heating a carbon steel bar.
Figure 3.40 Results of FEA computer simulation of surface temperature distribution versus cycle time in bar end heating. Material is AISI 1039. Frequency is 1 kHz.
FEA mesh and computersimulated magnetic field distribution at the final heating stage are shown in Figure 3.39. A concentration of magnetic flux lines clearly indicates a demarcation of the bar regions heated above and below the Curie point. Note the end effects and the complexity of the magnetic field distribution near both ends of the induction coil.
Sequential dynamics of surface temperature profiles and appearance of the end effects at different heating stages are shown on Figure 3.40.
During an initial heating stage, the whole bar is magnetic and end effect manifests itself as the end effect of a ferromagnetic body. With time, an interim heating stage takes place. At that stage, the surface layers of the bar are heated above the Curie temperature that is associated with the loss of magnetic properties. However, the subsurface regions retain their ferromagnetic properties with temperatures below the Curie point. A complex combination of end effects of magnetic and nonmagnetic bodies occurs during this stage. Finally, the bar end temperature exceeds the Curie point and the bar end effect represents the end effect of a nonmagnetic material.
Several subtleties are associated with an end effect applying electromagnetically short coils and relatively small coiltoworkpiece gaps (e.g., 2–5 mm). The basic appearance of end effect in these cases is illustrated using lower and higher frequencies in Figure 3.41 (a and b, respectively).
Figure 3.41 Illustration of electromagnetic end effect in a singleturn coil applying low frequency (a) versus high frequency (b).
The electromagnetic end effect is also responsible for an appearance of localized hot and cold spots using conventionally wound multiturn coils with a helix of turn windings. As an example, Figure 3.42 shows the results of computer modeling demonstrating the electromagnetic field distribution when heating ends of 250 mm diameter × 19 mm wall carbon steel pipe using a frequency of 500 Hz and a conventional multiturn coil design.
Figure 3.42 Coil field distributions attributed to a helix effect.
Because of a helix of coil winding, there is noticeably greater overhang at the coil top compared to its bottom. Variation in the electromagnetic field and, associated with it, induced heat source distribution within the pipe end region (top vs. bottom) can be clearly seen. Higher frequencies, using wide copper tubing and large pitches of turn winding will result in more pronounced circumferential variation of power density.
Heating trials confirm the results of computer modeling. For example, Figure 3.43a shows the appearance of a hot spot at the pipe end area, illustrating the presence of the helix effect. Efforts should be made to use coil designs with straightwounded coil turns, which will help eliminate or dramatically reduce the effect of the turn’s helix on temperature nonuniformity along the workpiece circumference that is associated with a variation in end effect.
Figure 3.43 The helix effect manifests itself in the hot spot at the pipe end area (a) and by realigning pucks (b).
The helix effect can also manifest itself by unanticipated movement of heated workpieces. As an example, Figure 3.43b shows realignment of puckshaped brass workpieces (36 mm diameter and 6 mm thick). Upon entering a multiturn induction coil, magnetic forces of the induction coil twisted pucks in accordance to the helix of turn windings.
When heating a rectangular workpiece, besides the distortion of the magnetic field in the slab’s end areas, similar distortion occurs at its edges (Figure 3.34). This phenomenon takes place owing to the electromagnetic transverse edge effect that plays a major role in obtaining the required temperature distribution across the slab or plate width.
When heating a homogeneous nonmagnetic body, the maximum value of the eddy current density is located on the surface of the slab’s central part (it does not, however, mean that the maximum temperature is always located there). The more pronounced the skin effect, the closer induced current flow matches the contour of the slab. Figure 3.44 shows the distribution of the electric field intensity in the slab’s transverse cross section with a highly pronounced skin effect (d/δ = 10, where the slab thickness d divided by δ is equal to 10) and when the skin effect is not pronounced (d/δ = 3).
Figure 3.44 Distribution of the electric field intensity (E) in the transverse cross section of the slab. (From V. Rudnev, Mathematical simulation and optimal control of induction heating of largedimensional cylinders and slabs, PhD Thesis, Department of Electrical Technology, St. Petersburg El. Engineering University, Russia, 1986; and V. Rudnev, D. Loveless, Longitudinal flux induction heating of slabs, bars and strips is no longer “black magic,” Part 2, Industrial Heating, February, 46–50, 1995.)
If the skin effect is pronounced (d/δ > 5), then the current density is approximately the same along the slab perimeter, except in the edge areas (2D corners, Figure 3.44). The edge area is usually (1.5–4.0)*δ long. Even though thermal surface losses at the edge area are higher than heat losses at the central part, the edge areas can be easily overheated compared to the central part. This occurs because in the central part, the heat sources penetrate from two sides (from two opposite wide surfaces), but at the edge areas, the heat sources penetrate from three sides (two surfaces and the side). When heating magnetic steel, aluminum, or copper slabs when the skin effect is quite pronounced, the heat surplus in the edge area often occurs.
If the skin effect is not distinct (d/δ < 3), then underheating of the edge areas will take place. In this case, the current’s path in the slab cross section does not match the contour of the slab and most of the induced currents close their loops earlier, without reaching the 2D corner area and sometimes even the edge areas (Figure 3.44). As a result, the power densities in edge areas will be less than the corresponding values in the central part of the slab. For example, in heating thick titanium, nonmagnetic stainless steel, or Nibased superalloy slabs or large RCS bars using relatively low frequency, the temperature of the corner areas in the final heating stage could be 20% lower compared to the temperature of the central part requiring the application of a dual frequency design [1,44,79,80,84].
The electromagnetic effect of joined materials with different properties (EEJ effect) occurs when two different metals are located in a common magnetic field. To simplify the study of this effect, let us consider the electromagnetic process in a conventional solenoid induction coil with two cylindrical billets (Figure 3.45).
Figure 3.45 Sketch of an IH system to study the electromagnetic effect of joined materials.
Assume that the billets have different electromagnetic properties (e.g., different ρ or μ_{r}). When two billets with different material properties are joined or placed close to each other inside an induction coil, the field between the two billets becomes distorted [1,77,85].
For illustration purposes, let’s review what happens when two carbon steel billets with substantially different physical properties (e.g., one billet was heated above the Curie point and became nonmagnetic and the other billet maintained its magnetic properties) are located in a multiturn solenoid inductor. The distribution of the surface power density in this case is shown in Figure 3.46. If the induction coil and both billets are long enough, then the magnetic field intensity at their central areas will be approximately the same and correspond to the coil current. At the same time, the surface power densities of the magnetic and nonmagnetic billets will be rather different (Figure 3.46a).
Figure 3.46 Electromagnetic effect of joined materials (EEJ). When heating magnetic vs. nonmagnetic materials (a) and when heating nonmagnetic billets having different electrical resistivities (b).
At the left end of the nonmagnetic billet (billet #1, Figure 3.45) and at the right end of the magnetic billet (billet #2), there will be a nonuniform power density distribution attributed to the end effects of the nonmagnetic and magnetic bodies (Figure 3.46a). Another area where the magnetic field is distorted and quite complicated is the transition area between the billets (the socalled billet’s joint area). At the right end of the nonmagnetic cylinder (billet #1), the power density sharply increases. At the left end of the magnetic cylinder (billet #2), the surface power density sharply decreases. This phenomenon is called the electromagnetic effect of joined materials with different properties (EEJ effect).
When discussing the EEJ effect, it is necessary to mention that this effect also occurs when both workpieces are nonmagnetic but have noticeably different electrical resistivities (Figure 3.46b). Figure 3.47 shows the power density distribution within billet #1 for the induction system shown in Figure 3.45. Both billets are nonmagnetic, but they have different electrical resistivities (ρ_{1} and ρ_{2}). When the electrical resistivity of billet #1 (ρ_{1}) is three times that of billet #2 (Figure 3.46b), then the surface power density in the joint area of the billet with higher electrical resistivity (i.e., the right end of billet #1) is reduced (Figure 3.47a). When ρ_{1} = 0.33 × ρ_{2}, there is an increase in surface power density there (Figure 3.47b).
Figure 3.47 Power density distribution along the length of billet #1 (frequency = 60 Hz; ρ = 1.1 μΩ, compare with Figure 3.46). Note substantially different appearance of EEJ effect when ρ of the billet #1 is threefold of the billet #2 (a) vs. the opposite case (b).
The EEJ effect does not usually play as important a role as skin effect and electromagnetic end and edge effects; however, there are some applications where this phenomenon may have an appreciable impact on the transient and final temperature distribution within the workpiece, especially when it is required to heat relatively thin but long billets just above the Curie point.
Thermal conductivity k designates the rate at which heat travels across a thermally conductive workpiece. A material with a high k value will conduct heat faster than a material with a low k. Depending on application specifics, a particular value of k can be beneficial or unfavorable. For example, choosing a material for an inductor’s refractory or a liner, a lower value of k is required corresponding to higher thermal efficiency and lower surface heat losses [73]. Conversely, when the k of the heated material is high, it is easier to obtain a uniform temperature distribution within the workpiece, which is important in through heating applications.
However, in selective heating applications (e.g., gear surface hardening or case hardening of shafts), a high value of k is quite often a disadvantage because of its tendency to promote heat transfer and equalize the temperature distribution within the workpiece. As a result of intense heat transfer, the temperature rise will take place not only in the region, which is to be hardened, but in adjacent areas as well, which are not. The temperature increase in the adjacent areas of the workpiece increases power consumption, making process less energy efficient and, in some cases, can negatively affect microstructural characteristics and residual stresses. A larger than needed amount of heated mass in the workpiece can also lead to an excessive distortion.
The Wiedermann–Franz law governs the relationship between the thermal conductivity (κ) and the electrical conductivity (σ) for the majority of pure metals and metallic materials, which is also a function of the temperature (T). However, some alloys (e.g., cast irons) may be exceptions from this general rule.
The values of the thermal conductivity of some commonly used metals are shown in Figure 3.48 [71,72]. As one may note, the thermal conductivity is a nonlinear function of temperature. Alloying and residual elements can have a measurable impact on k.
Figure 3.48 Thermal conductivities for some metals versus temperature.
The value of heat capacity C indicates the amount of energy that would have to be absorbed by the workpiece to achieve a unit of required temperature change. Mathematically speaking,
where dQ is the required energy and dT is the required temperature change. Heat capacity C is measured in J/(mol °C).
Heat capacity is closely related to a parameter called specific heat c, which represents the heat capacity per unit mass, meaning the amount of the required energy to be absorbed by a unit mass of the material to achieve a unit temperature increase. The c is measured in J/(kg °C) or Btu/(lb °F).
A higher value of specific heat corresponds to the greater required power to heat a unit mass to a unit temperature. The values of the specific heat of some commonly used metals are shown in Figure 3.49 [71,72].
Figure 3.49 Specific heat for some metals vs. temperature. (a) Specific heat vs. temperature for pure Cu, Al, and W; (b) for plain carbon steels SAE 1010, 1042, and stainless steel 304.
Similar to the electrical resistivity discussed in Section 3.1.1.1, the values of thermal conductivity k and specific heat c are also affected by chemical composition, the presence of residual elements, grain size, plastic deformation, prior heat treatment, and some other factors. As an example, Table 3.7 shows an appreciable variation of physical properties of some 7000series aluminum alloys at room temperature (20°C). Electrical resistivities of aluminum alloys are typically greater than those of pure aluminum (i.e., ρ of alloy 7075T6 is 93% greater than that of pure Al). Values of thermal conductivity and specific heat of alloys are usually lower than those of pure Al, though there are some exceptions.
Aluminum Alloy 
Electrical Resistivity 
Thermal Conductivity 
Specific Heat 
Temper Treatment 

μΩ*m 
W/(m*°C) 
J/(kg*°C) 

Pure aluminum 
0.027 
211 
933 
Commercial grade 
7005 
0.04 
166 
855 
O 
0.049 
137 
T6 

7050 
0.037 
180 
860 
T73 
7072 
0.029 
227 
893 
O 
7075 
0.052 
130 
960 
T6 
7075 
0.043 
155 
T73 

Variation range of property of alloy versus pure Al 
+7%; +93% 
−38%; +7.6% 
−8.4%; +3% 
Source: ASM, Metals Handbook, Vol. 2, Properties and Selection: Nonferrous Alloys and Pure Metals, ASM, Cleveland, 1979.
In IH, all three modes of heat transfer—conduction, convection, and radiation—are present [88–99].
Heat is transferred by conduction from the hightemperature regions of the workpiece toward the lowtemperature regions. The basic law that describes heat transfer by conduction is Fourier’s law,
where q _{cond} is heat flux by conduction, k is thermal conductivity, and T is temperature.
As one can see from Equation 3.22, according to Fourier’s law, the rate of heat transfer in a workpiece is proportional to the temperature gradient (temperature difference) and the thermal conductivity of the workpiece. In other words, large temperature gradients between surface and core (which, e.g., typically takes place during surface hardening) and a high value of k result in intensive heat transfer from the hot surface of the workpiece toward the cold core. Conversely, the rate of heat transfer by conduction is inversely proportional to the distance between regions with different temperatures.
In contrast to thermal conduction, heat transfer by convection is carried out by fluid, gas, or air (i.e., from the surface of the heated workpiece to the ambient area). The wellknown Newton’s law can describe convection heat transfer. This law states that the heat transfer rate is directly proportional to the temperature difference between the workpiece surface and the ambient area,
where q _{conv} is heat flux density by convection (W/m^{2} or W/in.^{2}), α is the convection surface heat transfer coefficient [W/(m^{2}.°C) or W/(in.^{2}.°F)], T _{s} is surface temperature (°C or °F), and T _{a} is ambient temperature. The subscripts “s” and “a” denote surface and ambient, respectively.
The convection surface heat transfer coefficient is primarily a function of the thermal properties of the workpiece, the thermal properties of the surrounding gas, air or fluids (i.e., quenchants), and their viscosity or the velocity of the workpiece if it rotates or moves at high speed (e.g., induction heat treating of a strip, wire, or spinning disk).
Value of convection losses can vary dramatically depending on the temperatures, surface conditions, and whether it is free or forced convection. Remember that in a number of IH applications (e.g., heating of strips, wires, rotating disks, gears, or shafts), convection heat transfer cannot be considered as free convection. In some strip coating applications, the strip travels at a speed up to 5 m/s (16 ft/s). Therefore, the heat losses attributed to forced convection are much higher (e.g., often 5–10 times higher) than free convection losses of the stationary heated workpiece.
There are several empirical formulas that provide a rough engineering estimate of the free convection losses q _{conv}, including Equations 3.24 [100] and 3.25 [11] for billet heating
where T _{s} and T _{a} are surface temperature and ambient temperature, correspondingly, in degrees Celsius. For example, free convection heat loss density from the workpiece surface heated to 600°C (1112°F) into the surrounding atmosphere (T _{a} = 20°C) is as shown in Table 3.8.
According to Equation 3.17 
According to Equation 3.18 

q _{conv} = 1.86 × 580^{1.3} = 7.3 × 10^{3} W/m^{2} 
q _{conv} = 1.54 × 580^{1.33} = 7.45 × 10^{3} W/m^{2} 
or q _{conv} = 7.3 kW/m^{2} 
or q _{conv} = 7.45 kW/m^{2} 
Convection heat transfer plays a particularly important role in the quenching process where the surface heat transfer coefficient represents the cooling intensity.
In the third mode of heat transfer, which is thermal radiation, the heat may be transferred from the hot workpiece into surroundings including a nonmaterial region (vacuum). The effect of heat transfer by radiation can be introduced as a phenomenon of electromagnetic energy propagation attributed to a temperature difference. The Stefan–Boltzmann law of thermal radiation, which states that the heat transfer rate by radiation is proportional to a radiation loss coefficient C _{s} and the value of
, governs this phenomenon.The radiation heat loss coefficient includes emissivity, radiation shape factors (the view factors), and surface conditions. For example, the value of emissivity increases with an increase in surface oxidation. At the same time, polished metal radiates less heat to the surroundings than nonpolished metal. A comparison of emissivity of some commonly used polished versus nonpolished metals is shown in Table 3.9.
Surface Condition 
Aluminum 
Carbon Steel 
Copper 
Brass and Zinc 

Polished 
0.042–0.053 
0.062 
0.026–0.042 
0.03–0.039 
Nonpolished, oxidized 
0.082–0.40 
0.71–0.8 
0.24–0.65 
0.21–0.50 
The radiation heat loss coefficient can be approximated as C _{s} = σ_{s}ε, where ε is the emissivity of the metal and σ_{s} is the Stefan–Boltzmann constant [σ_{s} = 5.67 × 10^{−8} W/(m^{2} K^{4})].
Thermal radiation loss density as a function of temperature and ε is shown in Figure 3.50a. Similar to convection losses, there is a formula that provides a rough engineering estimate of the free radiation losses q _{rad},
Figure 3.50 (a) Variation of the radiation loss density (kW/m2) versus temperature and emissivity. (b) A comparison of convection and radiation heat losses.
where T _{s} and T _{a} are measured in degrees Celsius.
For example, heat radiation flux density (free radiation) of a carbon steel slab (ε = 0.8) heated to 1250°C (2282°F) into the surrounding atmosphere (T _{a} = 20°C) can be calculated as
The abovedescribed calculation of radiation heat loss is a valid assumption for classical workpiece geometry when there is free heat radiation from the heated body into the surroundings. However, there are some applications where the radiation heat transfer can be complicated and such a simple approach would not be valid. Complete details of all three modes of heat transfer can be found in several references [88–99].
As discussed previously, the heat transfer by convection and radiation typically reflects the value of heat loss. A high value of heat loss reduces the total efficiency of the induction heater. The analysis shows that convection losses are the major part of the heat losses in lowtemperature IH applications (i.e., aluminum, lead, zinc, tin, magnesium, and steel at a temperature lower than approximately 350°C) and in particular in cases of forced convection when the heated workpiece is rotating or moving fast.
Since thermal radiation losses are proportional to the fourth power of temperature, these losses are a significant part of the total heat losses in hightemperature applications (e.g., heating of steels, titanium, tungsten, etc.) when it is required to heat the workpiece to the temperatures suitable for hot working (Figure 3.50b).
Since the value of specific heat c represents the amount of the thermal energy required to be absorbed by a unit mass of the workpiece to achieve a unit temperature increase, an average value of specific heat c can be used for a ballpark estimate of the required workpiece power (P _{w}) to heat a given workpiece to an average temperature rise at the required production rate. Equation 3.27 can be used for this purpose,
where m is the mass of the heated body (kg), c is the average value of specific heat [J/(kg °C)], T _{in} and T _{f} are average values of initial and final temperatures (°C), and t is the required heat time (s).
For example, in order to heat a solid copper cylinder (0.1 m diameter, 0.3 m long) from room temperature (20°C) to a temperature of 620°C in 120 s, the power required to be induced within the workpiece P _{w} can be determined using Equation 3.27.
In this case, the mass of the heated body can be calculated as
where γ is the density (kg/m^{3}) (for copper, γ = 8.91 × 10^{3} kg/m^{3}), D is the diameter (m), and l is the billet length (m).
c = 420 J/(kg °C) can be used as an average value of the specific heat of copper in the temperature range 20 to 620°C. Therefore, applying Equation 3.27, the required power will be
In engineering calculations, some practitioners prefer to use the value of the heat content of the material to determine the value of P _{w}. Heat content is measured in kW hour/t. In this case, Equation 3.27 can be rewritten as
Figure 3.51 shows the values of the heat content for commonly used metals used with Equation 3.28 to determine the required power (P _{w}) to heat a copper billet (example above) based on the heat content value. From Figure 3.51, the required value of the heat content would be approximately equal to 70 kW × hour/t.
Figure 3.51 Heat content of metals at various temperatures.
Simplified formulas such as Equations 3.27 and 3.28 are very convenient to use in applications such as IH of classically shaped workpieces (e.g., billets, bars, slabs, blooms, etc.) where relatively uniform through heating is required. Such simplified formulas have the advantage of providing a quick estimate of the required power (P _{w}). However, numerical computation can provide much more precise estimation particularly in cases such as surface hardening, selective hardening, and induction reheating when initial temperature and required final temperature are not uniform (e.g., the case of induction slab/bar reheating after a continuous caster).
It is important to remember that power P _{w} does not represent the power at coil terminals (the socalled coil power). Equation 3.29 provides a correlation between the coil power P _{c} and the workpiece power P _{w}:
where η_{el} is electrical efficiency and η_{th} is thermal efficiency. Both, η_{el} and η_{th} are in the range of 0 to 1.
The value of η_{el} represents the ratio of the power induced in the workpiece P _{w}, to the total of P _{w} and electrical losses
where
includes power loss in the coil turns and kW losses generated in electrically conductive bodies located in the surrounding area (as well as transmission kW losses) and can be determined asThe value of
also includes kW losses associated with undesirable heating of support structures, tooling, guides, and other electrically conductive bodies (i.e., shunts, fixtures, etc.).As shown in Ref. [6], when heating a solid cylinder in an electromagnetically long solenoid coil, the value of η_{el} can be roughly estimated according to the formula
where
is an effective coil inside diameter, is an effective outside diameter of a cylinder, and δ_{2} are current penetration depths in the coil copper and cylinder (workpiece), respectively; ρ_{1} and ρ_{2} are the electrical resistivities of the coil and workpiece; and μ_{r} is the relative magnetic permeability of the heated cylinder.Equation 3.32 has been obtained under the following assumptions:
The ratio
is referred to as coil electrical efficiency factor. High η_{el} corresponds to a low value of this factor. Therefore, assuming that current cancellation does not occur, then the high η_{el} takes place when heating workpieces that are magnetic, have high electrical resistivity, and have the smallest possible gap between the coil and workpiece (D _{1}/D _{2} → 1). For example, the η_{el} when heating carbon steel cylinders below the Curie temperature is usually in the range of 0.8 to 0.95. In contrast, when heating billets made from silver or copper, η_{el} is typically in the range of 0.35 to 0.45.When heating rectangular bodies, including slabs and plates, instead of Equation 3.32, Equation 3.33 should be applied:
where F _{1} and F _{2} are effective perimeters of the coil opening and the heated slab, respectively.
The value of η_{th} in Equation 3.29 represents the amount of the thermal losses
compared to the heating power and can be determined asincludes heat losses from the workpiece surface attributed to radiation and convection as well as heat loss attributed to thermal conduction (e.g., the heat losses from the billet to watercooled guides).
The next section shows that the power induced in the workpiece P _{w} is not a constant during the heating cycle and varies depending on the change in ρ and μ_{r}. This is why, instead of using P _{w}, the value of
(meaning the average power per heating cycle or per particular process stage) is often applied.An application of thermal insulation or refractory can significantly reduce the heat losses. An accurate estimation of the value of the
can be determined with numerical computer modeling; at the same time, there are several empirical formulas that can allow rough estimation of those losses. For cylindrical coils with concrete blocks as a refractory, the value of thermal losses can be determined as shown [101]:where
is the heat loss from the workpiece surface (kW), l is the coil length (cm), D _{1} is the inside diameter of the induction coil (cm), and D _{3} is the inside diameter of the refractory (cm).Total efficiency of the induction coil (η) is a combination of both coil thermal efficiency (η_{th}) and electrical efficiency (η_{el}) according to Equation 3.36
An application of thermal insulation can improve η_{th} and reduce the heat losses significantly. At the same time, the use of refractory requires having available space for its installation that leads to a necessity of having larger coiltoworkpiece gaps. This worsens an electromagnetic coupling and as a result leads to a decrease in η_{el} (Figure 3.52a).
Figure 3.52 (a) Coil electrical and thermal efficiencies versus thickness of refractory. (b) Variation of power induced within the workpiece as a function of the heat time for two of the most common process modes: coil current is constant (solid line) and coil voltage is constant (dotted line).
Therefore, on the one hand, refractory allows to improve the thermal efficiency, but on the other hand, it reduces electrical efficiency. Therefore, a decision to use or not to use a refractory is always a reasonable compromise. In some cases, it is more energy efficient and cost effective not to use any refractory at all and, therefore, have the smallest possible coiltoworkpiece gap maximizing η_{el}. This is typical for majority of selective hardening, brazing, and soldering applications.
In other cases, it is advantageous to use a refractory and significantly decrease surface heat losses and more than compensate for the loss of η_{el} that is associated with greater coiltoworkpiece gap (i.e., IH before forging, rolling, and extrusion). Numerical modeling helps make an intelligent decision as to whether to use a refractory or not.
The dynamics of the IH are affected by several nonlinear factors, including but not limited to a variation of electromagnetic and thermal properties of the heated workpiece, power supply operational features, control mode, and system layout. For illustration purposes, it is beneficial to evaluate the dynamics of IH considering simplified conditions of heating a carbon steel cylinder located inside an electromagnetically long solenoid coil. Having this in mind, it is possible to recognize three classical process modes: coil voltage is constant, coil current is constant, and coil power is constant during the entire heating cycle. Figure 3.52b shows the variation of power versus time when heating a carbon steel cylinder from ambient to forging temperatures. A detailed study of these modes can be found in Refs. [6,9]; thus, only a brief discussion will be provided here.
It will be beneficial at this point to review a case study of heating a 76mmdiameter (3in.diameter) medium carbon steel bar using an inline multicoil induction system. Coil parameters are as follows: ID is 152 mm (6 in.); refractory thickness is 12 mm (0.5 in.); length of each coil is 1 m (40 in.); number of coils is 8; gap between coils is 0.3 m (12 in.); frequency is 1 kHz; and production rate is 65 mm/s (2.56 in./s). Figure 3.53 shows results of computer modeling of critical temperature profiles along the length of the induction line [318].
Figure 3.53 Thermal dynamics of IH of medium carbon steel bars (diameter is 76 mm [3 in.]). Frequency is 1 kHz. (From V. Rudnev, D. Loveless et al., Efficiency and temperature considerations in induction reheating of bar, rod and slab, Industrial Heating, June, 2000.)
At the initial stage of the heating cycle, the entire workpiece is magnetic, μ_{r} is quite large, δ_{2} is respectfully small, and therefore, the skin effect is pronounced. At the same time, because of the relatively low temperature, the heat losses from the cylinder surface at this stage are low. The induced power appears in the fine surface layer of the workpiece leading to a rapid increase in the surface temperature with practically no change at the core.
Figure 3.54a shows the temperature and power density (heat source) distribution along the radius of the workpiece at an initial heating stage of case study shown in Figure 3.53. The maximum temperature is located at the surface. Intensive heating and the existence of a large temperature differential within the workpiece is typical for this stage. As one can see from Figure 3.54a, the temperature profile does not match the power density profile because of thermal conductivity k, which spreads the heat from the surface toward the core.
Figure 3.54 Temperature and power density profiles at different stages of IH a carbon steel bar for a case study shown on Figure 3.53. (a) Initial stage, (b) intermittent stage, (c) final stage. (From V. Rudnev, D. Loveless et al., Efficiency and temperature considerations in induction reheating of bar, rod and slab, Industrial Heating, June, 2000.)
During this stage, the η_{el} increases because of an increase in ρ of the metal with temperature (Figure 3.3). At the same time, μ_{r} remains relatively high, and a slight reduction of μ_{r} does not affect the rise in η_{el}. After a short time, coil electrical efficiency reaches its maximum value and then starts to decline. The surface reaches the Curie temperature first, and after that, the heat intensity at the surface significantly decreases. This takes place because of several factors. Those that played particular important roles are as follows:
Figure 3.54b shows the temperature profile and heat source distribution along the radius of the cylinder sometime after the surface temperature exceeds the Curie point (second heating stage that is called the transient or intermittent stage). At this stage, the ρ of the carbon steel has increased approximately two to threefold compared to its initial value. A decrease of μ_{r} and an increase of ρ cause a 6 to 10fold increase in δ compared to its initial value [1]. Most of the power is now induced in the surface and the internal layers of the workpiece. This portion of the heating cycle can be characterized as the dualproperty stage when heating carbon steels or any other magnetic materials to temperatures above the Curie point.
The workpiece surface becomes nonmagnetic; however, its internal layers remain magnetic. This stage takes place while the thickness of the nonmagnetic layer is less than the δ in hot steel. Power density has a unique wavelike shape that is different from the classical exponential distribution. Figure 3.54b shows that a maximum of the heat sources occurs at the cylinder surface. Then, the power density decreases toward the core. However, at a distance of approximately 1.4 mm from the surface (in this particular case), it starts to increase again. This takes place because the carbon steel retains its ferromagnetic properties at this distance. It is necessary to mention here that under certain conditions, the maximum value of heat sources can be located in a subsurface layer of the workpiece and not on its surface.
Finally, the third stage (a nonmagnetic stage) takes place. The thickness of the nonmagnetic surface layer exceeds the value of 2*δ in hot steel, and the dualproperty phenomenon becomes less pronounced and will finally disappear. The power density will then have its classical exponential distribution (Figure 3.54c).
The existence of the three stages of IH results in variation of both the workpiece power P _{w} and coil power P _{c} during the process of heating. All three stages should be taken into consideration when evaluating the process parameters.
Surface power density for ballpark estimation purposes can be calculated for a magnetic body heated inside the infinitely long solenoid inductor as a function of the magnetic field intensity, frequency, ρ, and μ_{r} according to Equation 3.37 [6],
where p _{0} is surface power density measured in W/m^{2}, H _{surf} is the magnetic field intensity at the surface of the workpiece, ρ is electrical resistivity, μ_{surf} is relative magnetic permeability at the workpiece surface, and F is frequency. All units are according to the SI system. Precise calculation of all the major process parameters including heating conditions can be conducted only by numerical computer modeling.
To finalize the discussion on the dynamics of IH, it should be stated that the abovediscussed basic phenomena were illustrated assuming using a solenoid inductor. Particular application often has specific features that might affect the dynamics of the induction process and establish its uniqueness. Some of those features are discussed in subsequent sections of this book.
Mathematical modeling is one of the major factors in the successful design of IH systems. The current production environment does not allow the luxury of process design via trial and error. Computer simulation enables IH specialists to quickly determine process details, which could be costly, time consuming, and difficult or impossible to resolve experimentally. Simulation enables prediction of how different interrelated and nonlinear factors could affect the transitional and final thermal conditions of the heated component. It also helps determine what must be accomplished to improve process effectiveness to establish the most appropriate process recipes. Computer modeling provides a comfort factor when designing new systems, avoids unpleasant surprises, shortens the learning curve, and reduces development time.
The first step in any computational modeling is to define a set of relevant quantities that is required to simulate and select the governing equation(s). Theoretical models may vary from a simple handcalculated formula to a very complicated numerical analysis, which can require several hours of computational work using modern supercomputers. The choice of a particular theoretical model depends on several factors, including the complexity of the engineering problem, required accuracy, time limitations, and cost.
Before an engineer starts to provide a mathematical simulation of any process, it is necessary to have a sound understanding of the nature and physics of that process. Engineers should also be aware of the limitations of applied mathematical models, assumptions, and possible errors and should consider correctness and sensitivity of the chosen model to sometimes poorly defined parameters such as boundary conditions, material properties, or initial temperature distribution. One model can work well in certain applications and give unrealistic results in others. Underestimation of the process features or oversimplified assumptions can lead to an incorrect mathematical model (including chosen governing equations) that might fail to provide the required accuracy of calculations.
It is important to remember that any computational analysis can at best produce only results that are derived from the governing equations. Therefore, the first and the most important step in any mathematical simulation is to choose an appropriate theoretical model that properly describes the technological process or phenomenon.
Numerical computer modeling allows predicting how different factors may influence the transitional and final heating conditions and what must be accomplished in the design of the induction system and process recipe to improve system effectiveness and guarantee the desired heating results.
As mentioned above, IH is a multiphysical phenomenon comprising a complex interaction of electromagnetic, heat transfer, metallurgical phenomena and circuit analysis, which are tightly interrelated because the physical properties of materials depend strongly on magnetic field intensity and temperature as well as chemical composition. Following materials concentrate on mathematical modeling of the electromagnetic field and thermal processes that occur during IH. Metallurgical aspects are reviewed in Section 4.1.
Note: This section requires from the reader knowledge of certain aspects of mathematics including differential calculus, integral calculus, and vector analysis, and can be skipped by the firsttime reader or reader with limited knowledge of advanced mathematics and applied physics.
The technique of calculating electromagnetic field depends on the ability to solve Maxwell’s equations. For general timevarying electromagnetic fields, Maxwell’s equations in differential form can be written as [58–65]
where E is electric field intensity, D is electric flux density, H is magnetic field intensity, B is magnetic flux density, J is conduction current density, and ρ^{charge} is electric charge density.
Equations 3.38 through 3.41 consist of special symbols: ∇· and ∇×. These two useful symbols along with the third symbol of ∇ are popular in vector algebra and are used to shorten an expression of particular differential operation without having to carry out the details. For example, in the rectangular coordinate system, these symbols represent the following mathematical operations:
where i, j, and k are the unit vectors for the X, Y, and Zaxes, correspondingly.
The fundamental laws governing the general timevarying electromagnetic field (Equations 3.38 through 3.41) can be written not only in differential form but in integral form as well by applying Stokes’ theorem [58–65,103,104]. Different numerical calculation methods apply different forms of Maxwell’s equations. For example, finite element and finite difference methods typically use a differential form of the Maxwell equations. In contrast, an integral form is usually applied with the boundary element method.
Maxwell’s equations not only have a purely mathematical meaning, they have a concrete physical interpretation as well. For example, Equation 3.38 states that the curl of H always has two sources: conductive (J) and displacement ρ^{charge} charge currents. A magnetic field is produced whenever there are electric currents flowing in surrounding objects. From Equation 3.39, one can conclude that a time rate of change in magnetic flux density B always produces the curling E field and induces currents in the surrounding area and, in other words, it produces an electric field in the area where such changes take place. The minus sign in Equation 3.39 determines the direction of that induced electric field. This fundamental result can be applied to any region in space.
Let’s review how Equations 3.38 and 3.39 can be used to support the basic explanation of some of the electromagnetic processes taking place in IH. The application of alternating voltage to the induction coil results in the appearance of an AC in the coil circuit. According to Equation 3.38, an alternating coil current will produce in its surroundings an alternating (changing) magnetic field that will have the same frequency as the source current (coil current). That magnetic field’s strength depends on the magnitude of coil current, the coil geometry, and the coiltoworkpiece electromagnetic coupling. The changing magnetic field induces eddy currents in the workpiece and in other objects that are located near that coil. According to Equation 3.39, induced currents have the same frequency as the source coil current; however, their direction is opposite that of the coil current. This is determined by the minus sign in Equation 3.39. According to Equation 3.38, alternating eddy currents induced in the workpiece produce their own magnetic fields, which have opposite directions to the direction of the main magnetic field of the coil. The total magnetic field of the induction coil is a product of the source magnetic field and induced magnetic fields.
Equation 3.38 suggests that there can be an undesirable heating of tools, fixture, cabinets, fasteners, or other electrically conductive structures located near the induction coil.
An analyst should pay attention to such simple relations as Equation 3.40 or 3.41. The popular saying, “The best things come in small packages,” is particularly true here. The short notation of Equation 3.40 has real significance. To say that the divergence of magnetic flux density is zero is equivalent to saying that B lines have no source points at which they originate or end; in other words, B lines always form a continuous loop. A clear understanding of such a simple statement will allow one to explain and avoid many mistakes in dealing with the IH of workpieces with irregular geometries.
The abovedescribed Maxwell’s equations are in indefinite form because the number of equations is less than the number of unknowns. These equations become definite when the relations between the field quantities are specified. The following constitutive relations are additional and hold true for a linear isotropic medium.
where the parameters ε, μ_{r}, and σ denote, respectively, the relative permittivity, relative magnetic permeability, and electrical conductivity of the material; σ = 1/ρ, where ρ is electrical resistivity. The constant μ_{0} = 4π × 10^{−7} H/m [or Wb/(A × m)] is called the permeability of free space (the vacuum), and similarly the constant ε_{0} = 8.854 × 10^{−12} F/m is called the permittivity of free space. Both relative magnetic permeability μ_{r} and relative permittivity ε are nondimensional parameters and have very similar meanings and were discussed earlier.
By taking Equations 3.45 and 3.47 into account, Equation 3.38 can be rewritten as
For most practical applications of the IH of metallic materials, where the frequency of currents is less than 10 MHz, the induced conduction current density J is much greater than the displacement current density ∂D/∂t, so the last term on the righthand side of Equation 3.48 can be neglected and Equation 3.48 can be rewritten as
After some vector algebra and using Equations 3.38, 3.39, and 3.46, it is possible to show that
Since the magnetic flux density B satisfies a zero divergence condition (Equation 3.40), it can be expressed in terms of a magnetic vector potential A as
And then, from Equations 3.39 and 3.52, it follows that
Therefore, after integration, one can obtain
where φ is the electric scalar potential. Equation 3.47 can be written as
where J _{s} = −σ∇φ is the source (excitation) current density in the induction coil. Taking the material properties as being piecewise continuous and neglecting the hysteresis and magnetic saturation, it can be shown that
It should be mentioned here that for the great majority of IH of steels (such as through hardening and IH before forging, rolling, and extrusion), a heat effect attributed to hysteresis losses does not typically exceed 6%–8% compared to the heat effect attributed to eddy current losses. This is so, because in such applications during the majority of the heating cycle, the surface temperature of heated workpiece is well above the Curie temperature—thus being nonmagnetic. Therefore, an assumption of neglecting the hysteresis is valid.
However, in some cases, such as induction tempering, paint curing, stress relieving, heating before galvanizing, and lacquer coating, hysteresis losses can be quite significant compared to eddy current losses. In these cases, hysteresis losses should be taken into account since they can contain a significant portion of heat sources.
It can be shown that for the great majority of IH applications, it is possible to further simplify the mathematical model by assuming that the currents have a steadystate quality. Therefore, with this assumption, it is possible to conclude that the electromagnetic field quantities in Maxwell’s equations are harmonically oscillating functions with a single frequency. Thus, a timeharmonic electromagnetic field can be introduced. This field can be described by the following equations, which are derived after some vector algebra from Equations 3.50, 3.51, and 3.56, respectively.
where ∇^{2} is the Laplacian, which has different forms in Cartesian and cylindrical coordinates. In Cartesian coordinates,
In cylindrical coordinates (axissymmetric case),
In other words, an assumption of harmonically oscillating currents with a single frequency means that harmonics are absent in both the impressed and induced currents and fields. The governing equations (Equations 3.57 through 3.59) for the timeharmonic field with the appropriate boundary condition can be solved with respect to H, E, or A.
Note: Though the timeharmonic representation of electromagnetic field is a sufficiently accurate approximation for modeling the majority of IH applications, there are some cases when this assumption would not be appropriate and can dramatically affect the modeling accuracy.
Equations 3.57 through 3.59 are valid for general 3D fields and allow one to find all of the required design parameters of the induction system such as current, power, coil impedance, and so on.
Although there is considerable practical interest in 3D problem solving, a great majority of engineering projects in IH tend to be effectively handled with a combination of twodimensional (2D) assumptions. Several factors discourage 3D field consideration. This includes but is not limited to the following factors:
For many IH applications, the quantities of the magnetic field (such as magnetic vector potential, electric field intensity, and magnetic field intensity) may be assumed to be entirely directed. For example, in the longitudinal cross section of the solenoid coil, both A and E vectors have only one component, which is entirely Zdirected. In the case of a transverse section, H and B vectors have only one component (Figure 3.55). This allows one to approximate the 3D field consideration to a combination of 2D forms. For example, in the case of magnetic vector potential, Equation 3.57 can be expressed for a 2D Cartesian system as
Figure 3.55 E, A and H, B field representations in a cylindrical system.
and for an axissymmetric cylindrical system as
The boundary of the region is selected such that the magnetic vector potential A is zero along the boundary (Dirichlet condition) or its gradient is negligibly small along the boundary compared to its value elsewhere in the region (Neumann condition ∂ A/ ∂n = 0). Therefore, the heat transfer equation that is discussed in Section 3.4.2 and Equation 3.63 with their initial and boundary conditions fully describe the electrothermal processes in a great majority of conventional cylindrical induction heat treatment systems.
By using analogous vector algebra manipulations, it is possible to obtain governing equations similar to Equations 3.62 and 3.63 that can be formulated with respect to E, B, or H. Therefore, any given electromagnetic problem in IH may be worked in terms of A, E, B, or H. Part of the art of mathematical modeling of electromagnetic fields derives from the right choice of field representation, which could be different for different applications.
Partial differential equations that are formulated with respect to A or E are very convenient for describing the electromagnetic field in a longitudinal cross section of the IH system. However, the electromagnetic field distribution in a transverse cross section of the workpiece can be more conveniently described by governing equations formulated with respect to B or H [1,77,85,103]. Field representations that are typically used for describing electromagnetic processes in applications applying cylinder multiturn coils are shown in Figure 3.55.
In general, the transient (timedependent) heat transfer process in a metallic workpiece can be described by the Fourier equation:
where T is temperature, γ is the density of the metal, c is the specific heat, k is the thermal conductivity, and Q is the heat source density associated with induced eddy currents per unit time in a unit volume (socalled heat generation). This heat source density is obtained as a result of solving the electromagnetic problem. As one might conclude from Figures 3.48 and 3.49, both k and c are nonlinear functions of temperature.
Equation 3.64, with suitable boundary conditions and initial condition, represents the 3D temperature distribution at any time and at any point in the workpiece. The initial temperature condition refers to the temperature profile within the workpiece at time t = 0. The initial temperature distribution is usually uniform and corresponds to the ambient temperature. In some cases, the initial temperature distribution is nonuniform because of the residual heat accumulated after the previous technological process (i.e., preheating, partial quenching, or continuous casting).
For most IH applications, boundary conditions combine the surface heat losses owing to convection and thermal radiation (Figure 3.50). In this case, the boundary condition can be expressed as
where ∂T/∂n is the temperature gradient in a direction normal to the surface at the point under consideration, α is the convection surface heat transfer coefficient, C _{s} is the radiation heat loss coefficient, Q _{s} is the surface loss (i.e., during quenching or as a result of workpiece contact with cold rolls or watercooled guides, etc.), and n denotes the normal to the boundary surface. As one may see from Equation 3.65, the heat losses at the workpiece surface are highly nonlinear.
If the heated body is geometrically symmetrical, then the Neumann boundary condition can be formulated along the axis of symmetry
The Neumann boundary condition implies that the temperature gradient in a direction normal to the axis of symmetry is zero. In other words, there is no heat exchange at the axis of symmetry. This boundary condition can also be applied in the case of a perfectly insulated workpiece.
In the case of heating a cylindrical workpiece, Equation 3.64 can be rewritten as
Equation 3.64 can be shown in Cartesian coordinates (i.e., heat transfer in slab, strip, or plate) as
Equations 3.67 and 3.68 with boundary conditions (Equations 3.65 and 3.66) are the most popular equations for mathematical modeling of the heat transfer processes in IH and heat treatment applications.
Thermal conductivity of some materials is directionally dependent, exposing strong multidimensional nonisotropic properties (i.e., laminations, composites, etc.). In this case, k becomes the tensorial quantity:
However, in the great majority of cases involving IH of metallic workpieces, k can be assumed to be isotropic.
The analytical methods and equivalent circuit coil design methods popular in the 1960s and 1970s no longer satisfy the modern analyst because of the inherent restrictions. The designer must be aware that, in many applications, erroneous and inadequate results can be obtained when such methods are used. Rapid development of computer technology and the increasing complexity of IH applications have significantly restricted the use of simple formulas and analytical and seminumerical modeling methods. These methods can be useful only in obtaining approximate results in sufficiently simple cases (e.g., classical geometries).
Rather than using simplified computational techniques with many restrictions and poor accuracy, modern IH specialists apply highly effective numerical methods such as finite difference, finite volumes, finite element, edge elements, mutual impedance, boundary element methods, and others. These methods are widely and successfully used in the computation of electromagnetic and heat transfer processes. Each of these methods has certain advantages and has been used alone or in combination with others.
Because of the extraordinarily large amount of information that is available in the specialized scientific literature, even an experienced engineer/analyst can be easily confused when discussing the indepth nuances of computer modeling. Therefore, we briefly discuss the modern electroheat numerical computation techniques while simplifying the materials so IH and material science analysts who might have limited experience in numerical modeling understand them. Thus, our goal is to provide the reader with a general orientation on advanced numerical simulation methods.
Before we discuss the features and applications of some of the most popular numerical methods, it is necessary to point out one of their important qualities: all numerical methods provide an approximate solution to the modeled problem (including heat transfer and electromagnetic problems). Therefore, there is always the danger of obtaining inappropriate results when those methods are used. The fact that the numerical solution is always approximate and not absolutely accurate should not discourage engineers from using numerical methods. On the contrary, it should stimulate them to carefully study the features of these methods and to transform them into a powerful computational tool that will allow analysts to control the accuracy of the simulation and to produce information that cannot be measured or obtained by using analytical, semianalytical, or other kinds of methods, including physical experiments.
It is wise to remember that the correct use of numerical methods will provide approximate, but sufficiently accurate engineering solutions that will satisfy the requirements of modern technology from a practical standpoint.
Many mathematical modeling methods and programs exist or are under development. Work in this field is done in universities and research institutes, inside large companies such as Inductoheat Inc., and by specialized software companies such as Magsoft Corp., Integrated Engineering Software Inc., Infolytica Inc., ANSYS, COMSOL, SYSWELD, Vector Fields Inc., QuickField, and others.
For each problem or family of similar problems, certain numerical methods or software are preferred. There is not a single universal computational method that fits all cases and is optimum for solving all IH problems. Our experience in the use of different numerical methods has shown that it is preferable to apply different methods and software rather than to search for one universal program for solving a wide variety of tasks.
The right choice of computational method and software depends on the application and features of the specific problem or phenomenon. It is important for the designer of the IH equipment to know the advantages and limitations of different computer modeling techniques. This will allow the analyst to select the most appropriate computational tool.
Because of space limitations, this chapter does not review an exhaustive list of the methods available for electromagnetic field and heat transfer calculations. There are numerous publications that describe different mathematical modeling techniques. An interested reader can study the most popular computational techniques used for simulation of heat transfer and electromagnetic processes in an extended reference list provided in this book or on the Internet. Here, we briefly review only some of these methods.
The finite difference method (FDM) was the earliest numerical technique [89,103–109] used for mathematical modeling of different processes. The FDM has been used extensively for solving both heat transfer and electromagnetic problems. It is particularly easy to apply this method for modeling cylindrical or rectangular bodies.
The orthogonal mesh discretizes the area of modeling (i.e., induction coil, workpiece, flux concentrator, etc.) into a finite number of nodes (Figure 3.56). Because of the orthogonal mesh, the discretization algorithm is quite simple. An approximate solution of the governing equation is found at the mesh points defined by the intersections of the lines.
Figure 3.56 Rectangular mesh network (grid) used in finite difference approximation.
The computation procedure consists of replacing each partial derivative of the governing equations (Equation 3.62, 3.63, 3.67, or 3.68) by a finite difference “stencil” that couples the value of the unknown variable (i.e., temperature or magnetic vector potential) at an approximation node with its value in the surrounding area. This method provides a pointwise approximation of the partial differential equation based on utilizing the Taylor’s series. FDM is quite a universal modeling tool, and its popularity is attributed to its generality and relative simplicity of application.
By Taylor’s theorem for two variables, the value of a variable at a node on the mesh can be expressed in terms of its neighboring values and separation distance (called a space step) h as in the following expressions (stencils):
Here, the notation O(h) is used to show that the truncation error involved in the approximation is on the order of h. O(h) is a linear function. Similarly, O(h ^{2}) is for the approximation error on the order of h ^{2}, which is much smaller, leading to a more accurate solution than one on the order of h.
Substitution of the finite difference stencils into the electromagnetic and heat transfer partial differential equations gives the local approximation. By assembling all local approximations and taking into account the proper initial and boundary conditions, one can obtain a set of simultaneous algebraic equations that can be solved with respect to unknown variables (i.e., T, A, E, H, or B) at each node of the mesh. The solution can be obtained either by iterative techniques or by direct matrix inversion methods. Since the matrices are sparsely occupied (nonzeros only in the neighborhood of the diagonal), some simplification in the computational procedure can be used.
As an example, we illustrate using FDM for modeling heat transfer processes for cylindrical billet heating (Fourier equation). The governing equation (Equation 3.64) can be rewritten as
For describing the heat transfer process in rectangular bodies (i.e., slab, plate, and bloom), Equation 3.64 can be rewritten as
As stated earlier, among other factors, c and k are functions of the temperature. The partial differential equation (Equation 3.73) may be expressed in a more concise form by introducing the finite difference operators
Substitution of Equations 3.75 and 3.76 into Equation 3.73 results in the finite difference format:
Figure 3.56 shows the rectangular mesh network. As mentioned above, the material properties are considered to be piecewise constants. Therefore, the coefficients of Equations 3.73 and 3.74 vary at different mesh nodes. The finite difference stencil with respect to the Zcoordinate can be written as
In FDM, it is important that the boundaries of the mesh region properly coincide with the boundaries of the appropriate regions of the IH system. Experience in using FDM in IH computations has shown that noncoincidence of the boundaries might have a noticeable negative effect on the accuracy of the calculation.
Approximations of the boundary conditions by Z = 0 and Z = ZZ are
where i, j, and τ are indexes corresponding to the Zaxis, the Raxis, and the time, respectively.
The finite difference expressions for differential operators with respect to the radius will be similar to Equations 3.78 and 3.79 [106].
When mesh boundaries do not coincide with boundaries of the components of the induction system, then the values corresponding to the temperature at the boundary nodes are the values they have at the neighboring nodes of the real boundary (Figure 3.56, bottom right).
Since the accuracy of the numerical computation depends on both the errors in the governing equation approximation and the error derived from approximating the boundary conditions, it is very important to treat boundary conditions at least as accurately as the governing equation.
Another factor that emphasizes the importance of a “good” approximation is related to various electromagnetic phenomena discussed in Section 3.1, including the skin effect. Rough approximation in these areas of high current concentration can have a detrimental effect on the overall accuracy of the calculations.
IH is a nonlinear timedependent (transient) process. There are several formats available to address these features of nonlinearity and time dependency. Each algorithm has its own advantages and disadvantages. The choice of a particular numerical procedure depends on several factors, including the specifics of the application.
Finite difference formats for the heat transfer transient problem range from explicit forms to implicit forms [106]. Implicit forms require solving a set of algebraic equations at each time step.
The explicit approximation is the simplest technique where the temperature distribution is obtained directly in a stepbystep manner. A forward difference approximation with respect to time leads to the explicit finite difference formulation
As one can see from Equation 3.80, the unknown temperatures corresponding to the (τ + 1) time step are obtained as functions of the known material properties, heat sources, and temperatures at time τ (Figure 3.57a). The temperatures are calculated after the first time step h _{τ}, which is found by the given initial condition (the initial temperature condition is often assumed to be ambient) and the appropriate boundary conditions. Therefore, the unknown temperatures are obtained explicitly from their initially known or previously calculated values.
Figure 3.57 Examples of simplified explicit (a) and implicit (b) forms for onedimensional approximation.
The ability to provide a stable and accurate numerical solution is primarily a concern when using an explicit finite difference format. Accuracy is a measure of the closeness of the numerical approximation to the exact solution [106]. The finite difference format (also called finite difference formulation) is said to be numerically stable if at sufficiently small time steps h _{τ} and space steps h. Equation 3.80 has a unique solution and that solution does not change its magnitude with small variations of h _{τ} and h.
It should be pointed out here that the stability condition depends on the properties of the finite difference format and is, in many cases, independent of the governing partial differential equation or physical phenomena. Unfortunately, explicit methods are accurate and stable only for certain relations between the time intervals, space steps, and values of material properties. Sometimes, those relations can contradict one another. The stability condition usually leads to extremely small time steps. Otherwise, a physically unrealistic oscillatory solution can occur.
With explicit formats, it is not unusual to have a situation where the decreasing time steps and space steps will not improve the solution but rather worsen it. This is a typical case of unstable or illconditioned systems. In such cases, the use of different stencils may help. For example, instead of a central difference stencil, a forward difference or backward difference approximation can be used and vice versa.
Thus, regardless of the simplicity and convenience of the explicit algorithms, a concern for obtaining an accurate and stable solution (particularly taking into consideration essentially nonlinear material properties) may limit the use of these algorithms.
Implicit methods are more popular because of their flexibility and ability to provide more stable results compared to explicit algorithms and to have a relatively independent choice of mesh parameters (Figure 3.57b). Several implicit methods were developed to reduce computational efforts. The use of any implicit method requires the calculation of a system of algebraic equations.
When using implicit methods for modeling heat transfer problems, the finite difference format can be written as [106]
The choice of the parameter ξ is a balance between accuracy and stability. The value of this parameter varies between 0 and 1. For ξ = 0.5, the wellknown Crank–Nicolson format represents an intermediate approximation of the partial derivatives (halfway between two levels of time τ and τ + 1). The complete implicit format is obtained when ξ = 1.
The implicit method is said to be unconditionally stable; however, certain computational oscillations could still appear when too coarse mesh and large time steps are used. The time step is restricted by the desired accuracy. The following finite difference implicit formats are commonly used for solving the transient heat transfer problem.
a. A Locally OneDimensional Format (Proposed by A. Samarskii [106])
The set of Equations 3.82 and 3.83 is said to be stable for all sizes of time step h _{τ}. The main restriction for choosing a large h _{τ} is avoiding significant truncation errors. Physically, Equations 3.82 and 3.83 can be interpreted as a complex combination of two heat transfer processes: along the Zaxis and along the Raxis [106]. The transition from time level τ to time level τ + 1 is assumed to be made in two stages using intermittent time step 0.5h _{τ}. This means that the transition from a known temperature field distribution of
to an unknown temperature is made through the intermediate temperature distribution of In each direction, Fourier equation is approximated implicitly with the necessity of solving two sets of simultaneous algebraic equations.After substituting the respective finite difference stencils into Equations 3.82 and 3.83 and after some simple algebraic operations, Equation 3.82 can be rewritten as
and, respectively, Equation 3.83 will be
where ζ, ψ, and υ are coefficients.
As one can see, the matrices of the algebraic Equations 3.84 and 3.85 are sparsely occupied, having a tridiagonal matrix structure, meaning that nonzeros occupy only the main diagonal and its neighborhood. Thanks to this feature, several simplified computational subroutines can be effectively used [106] to solve Equations 3.84 and 3.85.
b. Peaceman–Rachford Format [106]
Equation 3.86 is implicit in direction Z and explicit in R. However, Equation 3.87 is explicit in direction Z and implicit in R. A set of algebraic equations that corresponds to the Peaceman–Rachford format is similar to Equations 3.84 and 3.85.
As mentioned earlier, there have been two general techniques used for solving algebraic equations obtained after substitution of the finite difference stencils into the partial differential equations: iterative algorithms (such as the Jacobi method, Gauss–Seidel method, overrelaxation techniques, etc.) and direct methods. One of the most widely used methods for solving a tridiagonal matrix is the Gaussian twostep elimination method. This algorithm is quite effective and takes into account the features of the tridiagonal matrix requiring a minimum computer memory and a short execution time.
The optimal choice of mesh generation and time steps has a pronounced effect on the accuracy and stability of the calculation using any of the numerical modeling techniques; however, these parameters become particularly critical when using FDM. Certainly, in FDM as in any of the numerical techniques, for greater accuracy, the smaller space and time steps are recommended. In addition, it is quite clear that the large number of nodes results in a more complicated and timeconsuming solution. Therefore, there should be a reasonable compromise among mesh size, time steps, computation time, and accuracy of modeling. Naturally, it is recommended to select a smaller mesh size for regions with the greater gradients of variables (i.e., temperatures, magnetic vector potentials, or magnetic field intensities) and coarse mesh for areas where there is insignificant variation of variables.
The optimal combination of mesh parameters and time steps is usually determined by computational experiment. The calculations are provided for the different mesh sizes and time steps and the results are compared. If the comparison shows a large difference, then it is necessary to repeat the calculation for a finer mesh or smaller time steps until the difference between calculations is insignificant. The rule of thumb is as follows: if the computation is done correctly, the values of the unknown variables (e.g., temperature) should converge as the space mesh becomes finer and time steps become smaller.
It is wise to remember that the reduction of the space steps leads to the reduction of truncation error. However, unreasonably fine mesh is associated with a tremendous number of algebraic equations; however, a computer deals with only a limited number of arithmetic units. All these can lead to a crucial level of roundoff errors. Therefore, refining the space mesh and reducing the time steps can improve the accuracy of the computations, unless the roundoff errors become excessive (see Figure 3.58). The use of a double precision arithmetic is a helpful way to avoid computation failures caused by roundoff errors.
Figure 3.58 Correlation among the roundoff error, truncation error, and mesh size.
The abovementioned remarks regarding different aspects of mesh generation and computation errors are valid not only for FDM but for the majority of other numerical techniques as well.
When modeling processes that couple several different physical phenomena (e.g., electromagnetics, heat transfer, phase transformation, electromechanical, etc.), it is very attractive to use a single universal mesh. This might seem like a timesaving approach and would allow one to save time on the mesh generation. However, if the physical phenomena are inherently different, it is often more efficient to use different, phenomenonoptimized meshes.
The method of finite differences has been illustrated based on the most commonly used first and secondorder finite difference approximations (Equations 3.69 through 3.72). The accuracy of the numerical calculations may be improved by employing higherorder finite difference approximations. Such approximations allow reducing truncation error, but at the same time, this approach results in an increase of the number of nodes involved in the local approximation. And therefore a matrix of algebraic equations will no longer be tridiagonal but five or sevendiagonal.
The finite element method (FEM) is another group of numerical modeling techniques devoted to obtaining an approximate solution for different technical problems, including those encountered in IH. This versatile numerical technique was originally applied in mechanical engineering. Later, applications of FEM have expanded to other areas of engineering. It has become the most popular numerical tool for a variety of scientific and engineering tasks. The tremendous improvement in computer capabilities (particularly within the last three decades) has boosted the development of several variations of the FEM [104,110–125].
Some of these are as follows:
As described in Section 3.4.3.2, the FDM provides a pointwise approximation; however, the FEM provides an elementwise approximation of the governing equations. Certain finite element techniques might be better suited for certain problems. For example, the weighted residuals formulation has been very effectively used for computation of heat transfer phenomenon. Since IH is a complex combination of electromagnetic and heat transfer phenomena and taking into consideration that the use of FDM has been illustrated for modeling a heat transfer in Section 3.4.3.2, here we demonstrate the use of FEM to simulate electromagnetic processes.
The large number of papers on the subject of FEM applications for electromagnetic field computation makes it impossible to mention all of the contributions. Some of the proposed finite element models are similar in form. It should be mentioned here that P. Silvester and M. Chari [110,111] presented the first general nonlinear variational formulation of magnetic field analysis using FEM. Marked input into the earlier development of FEM was provided by W. Lord, A. Konrad, S. Salon, S. Udpa, J. Brauer, A. Bossavit, J. Sabonnadiere, W. Trowbridge, J. Simkin and many others. The following is a short description of one form of FEM. This approach is based on a combination the finiteelement concept proposed by W. Lord [113,114], S. Udpa and coworkers [120,123], and S. Gurevich and coworkers [103].
Because of the general postulate of the variational principle, the solution of electromagnetic field computation is obtained by minimizing the energy functional that corresponds to the governing equation (e.g., Equation 3.62 or 3.63) instead of solving that equation directly. The energy functional is minimized for the integral over the total area of modeling, which includes the workpiece, coil, electrically conductive tooling, fixtures, and surrounding area.
The principle of minimum energy requires that the vector potential distribution corresponds to the minimum of the stored field energy per unit length. As a result of that assumption, it is necessary to solve the global set of simultaneous algebraic equations with respect to the unknown, for example, magnetic vector potential at each node. The formulation of the energy functional, its minimization to obtain a set of finite element equations, and the solution techniques (the solver) were created for both 2D (Cartesian system) and axisymmetric cylindrical system.
Magnetic vector potential A in the 2D case (longitudinal cross section) acts in the direction of the current density J and is described by a 2D partial differential equation (Equation 3.62). The boundary of the region can be selected so that the magnetic vector potential A is zero along the boundary (Dirichlet condition) or Neumann condition (∂A/∂n = 0), meaning that its gradient is negligibly small along the boundary compared to the value elsewhere in the region.
The energy functional corresponding to the 2D governing equation (Equation 3.62) can be written in the following form [110,111]:
where V is the total area of modeling and J _{S} is a source current density. The first, second, and third terms inside the integrand represent the energy of the magnetic field, eddy currents, and the source current, respectively.
The minimization of the functional (Equation 3.88) corresponds to the solution of the 2D eddy current field problem with respective boundary conditions.
According to FEM, the area of study is divided into nonoverlapping finite number of subareas (numerous finite elements or mesh); therefore, the minimization of this functional provides minimization of energy at every node of each element.
Many geometric arrangements and shapes of finite elements are possible. Flexibility of their shapes allows them, in fact, to satisfy regions of practically any geometry of induction system.
The simplest 2D finite element is the firstorder triangle (Figure 3.59a). In the axisymmetric cylindrical case, such a finite element mesh may be represented as a set of rings. Each ring revolves around the axis of symmetry and has a triangular cross section (the socalled triangular torus element). The use of highorder isoparametric elements improves the accuracy of numerical approximation and allows reduction of the required total number of elements at the expense of some increase in computation time and algorithm complexity. Sixnode and eightnode isoparametric elements (Figure 3.59b) have become increasingly popular finite elements used in modeling both thermal and electromagnetic problems because they allow better treatment of skin effect problems.
Figure 3.59 Firstorder triangle (a) versus sixnode and eightnode isoparametric elements (b).
Space discretization is a very important aspect of FEM analysis. The following are some general remarks regarding finite element discretization (mesh generation) that has some similarities with FDM.
The area of study is subdivided into nonoverlapping finite elements (mesh). Figure 3.60 shows an example of the finite element mesh of two multiturn solenoid coils using triangular elements. This IH system is used for inline heating of long titanium bars. Copper end plate is located between coils to reduce their interaction. Notice a nonuniform density of the mesh.
Figure 3.60 Finite element mesh for two inline multiturn coils with copper end plate located between coils.
The size, shape, and positioning of these elements frequently depend on personal judgment. However, in order to obtain reasonable accuracy of the numerical solution, the finite element mesh has to be relatively fine (sizes of finite elements must be smaller) in the regions where the rate of change of the unknown (i.e., the magnetic vector potential) is high. Higher frequency and presence of ferromagnetic regions require a finer mesh. Material properties are postulated to be constant within each element but can be different from element to element.
It is beneficial to take advantage of the symmetry and/or periodicity of the modeled geometry (if applicable).
In many cases, when it is necessary to obtain the electromagnetic field distribution and temperature profile along the length of a cylindrical workpiece (longitudinal cross section), the FEM has been used in solving the governing equation with respect to magnetic vector potential (Equations 3.62 and 3.63).
Assuming that the local behavior of the electromagnetic field can be approximated by a linear law across each element and supposing that the chosen finite elements are firstorder parametric triangular elements, then the magnetic vector potential distribution within a triangular can be defined as
Based on 2D linear approximation laws, the coefficients α_{1}, α_{2}, and α_{3} are constant and can be calculated from the three independent simultaneous equations by assuming vertex values of A _{ l }, A _{ m }, A _{ n } of a magnetic vector potential A at the three nodes of a triangular. Therefore, the local set of equations can be rewritten as
The matrix notation of this set (Equation 3.90) can be written as
The determinant of the square matrix in Equation 3.91 can be introduced as a value of twice the triangular area. Knowing the geometry of elements and the magnetic vector potential at each node in every element, it is possible to obtain the value of A at any point inside the element. By extending a local approximation to all the elements that represent the total area of interest, it is possible to obtain an approximation for A throughout the modeling area.
Energy balance within the modeling area is determined by minimizing the energy functional at every node. This can be arranged by setting the first partial derivative of the functional with respect to each node, equal to zero. Instead of performing the nodebynode minimization of the functional, it is convenient to apply elementbyelement minimization. The total (global) energy associated with an entire modeling area equals the sum of the energies of all elements. As a result, a set of the simultaneous algebraic equations with respect to the unknown values of A at each node can be obtained.
After some algebraic operations, the local matrix equation, which represents the minimization of the energy functional within any triangular element, can be written as
where
Δ is a crosssectional area of a particular triangular.
After assembling all local matrices of finite elements and specifying the corresponding boundary conditions, a global matrix equation can be written as
It is necessary to mention here that there are several commonly used ways to specify the boundary conditions in Equation 3.98. One of the most popular techniques is called blasting the diagonal. This computational technique requires multiplying the diagonal terms of the equations representing the nodes where the value of the magnetic vector potential is known by a significantly large number (e.g., 10^{30}). At the same time, the corresponding righthand sides of those equations are replaced by known values of boundary conditions times the new diagonal. Such an artificial approach is very effective and easy to apply.
For the axisymmetric case, the local and global matrices will be similar to Equations 3.92 through 3.98 [114]. Parameters of the local matrix of Equation 3.92 for the axisymmetric problem (i.e., cylindrical system) are
where R _{c} is the radius of the finite element centroid and β_{ i } = b _{ i } + 2Δ/3R _{c}, i = 1, m, n
As with the FDM, the significant portion of the computation work for FEM consists of solving the large system of matrix equations. A global matrix in Equation 3.92 can be solved using either iterative methods or direct matrix inversion techniques while taking into consideration the matrix’s sparse nature and banded symmetry.
As mentioned above, the accuracy of the numerical approximation of the governing partial differential equations improves with a finer mesh but there might be a limit (see Figure 3.58). The optimum number of elements and their size and shape depend on the specifics of an application. In the stage of developing and testing FE code, a developer can judge the obtained accuracy of the FE approximation based on its comparison with the available analytical solution, appropriate physical models, or results of experiments.
After solving the system of algebraic equations and obtaining the distributions of the magnetic vector potential in the modeling region, it is possible to find all of the required output parameters of the electromagnetic field. For example, the induced current density in conductors can be determined as
The total current density in the conductor,
The magnetic flux density components B _{ x } and B _{ y } can be calculated from Equation 3.52 as follows [113,114]:
From Equation 3.106, the flux density can be obtained as
For the axisymmetric case of a cylindrical workpiece, the magnetic flux density components B _{ R } and B _{ Z } can be calculated as
Magnetic field intensity can be calculated as
Electric field intensity can be calculated as
Electromagnetic force density in currentcarrying conductors and the workpiece can be computed from the cross product of the vector of total current density and the vector of magnetic flux density:
It is also possible to compute the other characteristics such as stored energy, flux leakage, total power loss, coil impedance, and so on.
The abovedescribed FEM requires using the current density or current of the induction coil as the input parameter. This is often the case for induction hardening applications. At the same time, in some cases (e.g., several inductors connected electrically in parallel), it is more convenient to use the voltage of the coil as the input parameter. For cases such as these, special FEM subroutines have been developed.
The inductors involved in IH of billets, bars, slabs, and the like before metal warm and hot working including forging, upsetting, rolling, and extrusion are quite different compared to inductors for surface hardening. Induction billet/bar heaters are typically designed as multiturn solenoid coils of cylindrical or rectangular shape (Figure 3.1). Such induction heaters can consist of one or several inline coils (Figure 2.19). The total length of a system sometimes exceeds 10 m and can be as long as 30 m. The inside diameter of some coils can be as large as 1 m. Depending on the specifics of the application, coils can be fabricated as single or multilayer solenoids connected to a single or multiphase power source with conventional and complex circuit connections.
Two of the wellknown disadvantages of both FDM and FEM are the challenge of modeling long systems since an enormously large mesh is required that is associated with unrealistically large amount of computational work.
As an alternative to FDM and FEM, the mutual impedance method (MIM) can be used for IH for cylindrical systems. MIM applies an integral equation approach instead of using a differential formulation of the electromagnetic equations, which typically requires less computer memory and execution time.
Another advantage of using integral equations deals with the fact that the area of the integration (computation) is limited to surfaces of electrically conductive domains. In other words, the electrically conductive bodies of the IH system limit the mesh of discretization. Unlike FEM and FDM, integral formulations do not generally have to consider free space areas (such as air).
As with FDM and FEM, there are several different formulations of MIM devoted to the simulation of the IH problem. One of the earliest texts describing this technique was by E. Kolbe and W. Reis in 1962 [126]. Further development of this technique was done by O. Tozoni [127], R. Dudley and P. Burke [128], and several other researchers. A brief introduction to MIM is given here based on the approach discussed in Ref. [129].
Let’s first consider two axisymmetric multiturn coaxial coils (Figure 3.61) connected in series and driven by a sinusoidal voltage source. No harmonics are involved in the example. Both coils are placed around an axisymmetric nonmagnetic workpiece (i.e., copper, aluminum, or titanium billets). The electromagnetic field distribution in such a system can be described with respect to the current densities in the electrically conductive domains of the induction system by the Fredholm integral equation of the second kind:
Figure 3.61 Representation of the induction system for MIM. W represents the top half of the cylinder workpiece. H represents the induction heater.
where ρ_{ Q } = 1/σ_{ Q } is the resistivity of the element Q; R _{ Q } is the average radius of the element Q; J _{ Q } and J _{ P } are the current densities in the elements Q and P, respectively; M _{ QP } is the mutual inductance between elements Q and P, respectively, representing a mutual interaction of the elements (currentcarrying rings); S _{ P } are computation areas (p ∈ H, W, where H represents the induction heater and W represents the heated workpiece); and V _{ Q } is the source voltage of the element. The value of V _{ Q } is zero for all elements of the workpiece.
The method of solving the integral equation in the most complicated case has been described in Ref. [127]. The solution of Equation 3.112 in its simplest form [129] is presented here.
The electrically conductive regions of the induction system, including the induction heater and workpiece (Figure 3.61), are subdivided into appropriate elements. As with FEM, eddy current densities and material properties are assumed to be constant within each element.
If the skin effect in the coil is pronounced, then the multiturn induction coils can be considered to act as multiturn solenoids. The integral equation (Equation 3.112) can be rewritten as
where r _{ Q } = resistance of the element Q.
As seen from Equations 3.112 and 3.113, the Fredholm integral equation of the second kind is converted into an impedance equation representing the wellknown Kirchhoff’s law. After assembling equations that correspond to all electrically conductive elements of the induction system, the global set of impedance equations can be obtained.
A set of global equations representing the induction system is obtained below. According to the sketch shown in Figure 3.61, the IH system consists of two elements of the workpiece (elements P and Q) and two induction coils m and n connected in series. The global set of the impedance equations for this case is
where M _{ QQ }, M _{ PP }, M _{ nn }, and M _{ mm } are the selfinductances of the elements and coils, respectively; and M _{ QP }, M _{ Qn }, M _{ Qm }, M _{ PQ }, M _{ Pn }, M _{ Pm }, M _{ nQ }, M _{ nP }, M _{ nm }, M _{ mQ }, M _{ mP }, and M _{ mm } are the mutual inductances representing the interaction of all the currentcarrying elements.
The formulas for calculation of the various selfinductances and mutual inductances are given in Refs. [130–132]. The resistances of the rings (r _{ Q } and r _{ P }) and the resistances of the coils (r _{ n } and r _{ m }) can be calculated as
where K _{space} is the space factor of the turn winding of the coil and N is the number of turns of the coil n.
After some simple algebraic operations, the resulting matrix equation can be rewritten as
where a_{ W }, a_{ H } are matrices of the selfimpedances of the workpiece and induction coils and a_{ WH }, a_{ HW } are matrices of the mutual inductances.
Upon evaluation of the equations for the mutual inductances, it becomes obvious that M _{ PQ } = M _{ QP } and a_{ WH } = a_{ HW }. Therefore, the matrix of
is symmetric and the set of equations (Equation 3.117) can be rewritten as
where S _{ r } is a diagonal matrix consisting of the resistivities of elements of coils and the workpiece; S _{ x } is a square matrix of the selfinductances and mutual inductances; I _{ r }, I _{ x } are column matrices showing that the currents have both real and imaginary components (I = I _{ r } + j I _{ x }); and V _{ r }, V _{ x } are column matrices of the voltages (V = V _{ r } + jV _{ x }), which are similar to the column matrices of the currents.
The set of equation (Equation 3.118) can be rewritten as
Since the matrix S _{ r } is a diagonal matrix and the matrix S _{ x } is a symmetrical square matrix, in order to reduce the execution time and computer memory required for storage of all matrices, a special computational procedure for solving the set of equations (Equation 3.119) has been developed [129]:
After solving the set of Equations 3.120 and 3.121, one can obtain the coil currents and eddy currents as well as power densities, heat source distribution, and other important design parameters of the induction system including power, coil efficiency, and so on.
The MIM was extended for the computation of IH of magnetic workpieces by combining the mutual inductance method with the boundary element method (BEM) [77].
Our experience shows that MIM typically provides not as accurate a solution as FDM or, in particularly, FEM. This limits the wide use of this method in IH applications.
The fourth family of numerical techniques devoted to IH computation is called the BEM. This method started to be widely used for modeling IH in the late 1980s to the early 1990s. The mathematics required to discuss the BEM are more advanced than those needed for FDM, FEM, or MIM. A thorough discussion of the BEM and its forms is beyond the scope of this text. The interested reader will find several texts, conference proceedings, and journal articles [133–141] that describe various modifications of BEM.
Here, we just mention that when applying BEM (in contrast to FDM or FEM), an integral form of Maxwell’s equation is used as a governing equation for the electromagnetic problem. It is not required to make an artificial assumption for boundary conditions at infinitely propagating regions. The integral equation is complete as it is, thanks to the explicit appearance of the boundary values in the integrals. This allows one to take into consideration only electrically conductive domains in the computation. In this respect, BEM has similarity with MIM.
With the BEM, unknowns of the electromagnetic field are expressed in terms of an integral over the boundary of the area of interest. In this case, the problem of mathematical modeling of induction processes may be divided into two tasks: external and internal electromagnetic problems. Using an iterative procedure, both tasks can be solved. The internal problem describes the electromagnetic field distribution within the body of the workpiece. The external problem describes the field distribution in external regions.
With BEM, the meshing is only required at boundaries of the electrically conductive domains of the induction system (Figure 3.62). A computational procedure establishes the unknown surface qualities (i.e., current densities along surfaces) that would satisfy the global solution. This substantially simplifies meshing.
Figure 3.62 Boundary element mesh for two inline multiturn coils (compare with Figure 3.60).
Such advantages as the reduction of computation time, simplicity, and userfriendliness of mesh generation as well as good accuracy (particularly when dealing with nonferrous materials) make this technique quite attractive in certain applications.
As an example, Figure 3.63 shows a distribution of the magnetic field in the surroundings of two inline multiturn coils obtained using BEM. A solid cylinder bar made from titanium alloy Ti6Al4V is positioned inside induction coils. There is an electrically conductive end plate that is located between coils noticeably distorting the electromagnetic field between coils and reducing their electromagnetic interaction.
Figure 3.63 Magnetic field distribution around two inline multiturn coils.
BEM sometimes faces challenges in providing required accuracy of computations in cases of appreciably nonlinear processes. In such cases, a combination of BEM and FEM or BEM and FDM is preferable.
Both the electromagnetic and heat transfer phenomena are tightly coupled because of the interrelated nature of the material properties (see Sections 3.1.1 and 3.2.1) such as the following:
Obviously, these variations of physical properties make the IH process highly nonlinear, dictating the necessity of developing special computational algorithms that are able to deal with coupled phenomena. The time scales (time constants) of the electromagnetic and heat transfer processes in IH have different orders of magnitude. Electromagnetic processes are very fast, with time constants significantly less than 0.02 s (depending on the applied frequency). At the same time, in many IH applications, the heat transfer processes are substantially longer. For example, cycle time for induction billet heaters can easily exceed 900 s, depending on billet size.
There are several ways to couple the electromagnetic and heat transfer phenomena. The simplest method is called a twostep approach. Electromagnetic characteristics are obtained during the first step, allowing calculating the distribution of the Joule heat sources, which are used in solving the thermal problem assuming that the electrical resistivity and magnetic permeability do not change during heating. The twostep approach is known for its short execution time and modest computer memory requirements.
Even a cursory look at the behavior of the material properties in Figures 3.3 through 3.10 and 3.48 through 3.51 reveals the danger in using a twostep approach. As shown in Sections 3.1.1 and 3.1.2, the ρ of some metals can vary during the process of heating more than eight times. Furthermore, μ_{r} can vary more than 100 times. Therefore, when dealing with ferromagnetic materials (e.g., carbon steels), an assumption of constant material properties during the entire heating cycle is a very rough postulation and can result in significant errors. Thus, this approach can be used in a limited number of lowtemperature heating applications dealing with nonmagnetic materials.
The most common approach to coupled electromagnetic and heat transfer phenomenon is called the indirect coupling method (also called the weakly coupling method). This method calls for an iteration process (Figure 3.64). According to this approach, it is assumed that temperature variations are not significant during certain time intervals of the heat cycle, meaning that the electromagnetic properties remain unchanged. During those sufficiently small intervals, the heat transfer calculations continue without correcting the heat sources. The temperature distribution within the workpiece obtained from the timestepped heat transfer computation is used to update the values of specific heat and thermal conductivity at each time step. As soon as the heat source variations become significant (because of the variations of ρ and μ_{r}), the convergence condition will no longer be satisfied, and recalculation of the electromagnetic field and heat sources will take place.
Figure 3.64 Block diagram of the indirect coupling process.
For most IH applications, an indirect coupling being the most popular approach is valid and very effective. However, in some cases, this approach could lead to appreciable errors. In these cases, the direct coupling method should be applied, requiring the formalization of a set of governing equations in such a way that the unknown parameters of the electromagnetic field (e.g., a magnetic vector potential or magnetic field intensity) and unknown parameters of the thermal process (i.e., temperature) will be part of a global matrix that will be solved simultaneously. Direct coupling results in an extremely intensive computer execution time and is memory intensive. It should be used only in cases where it is absolutely needed.
Even a cursory look at network meshes (Figure 3.56 vs. Figure 3.60 vs. Figure 3.62) reveals that the selection of one technique over another depends on the specifics of the particular induction application. It is easy, for example, to apply the FDM when the modeling area has simple geometries, such as cylindrical or rectangular. Because of the orthogonal grid, the modeling algorithm is simple and fast.
However, FDM is usually not as well suited as FEM for the simulation of complex boundary configurations or in the case of a mixture of materials and forms. In this instance, FEM has a distinct advantage over FDM.
Superficially, the FDM and FEM appear to be different; however, they are closely related. As outlined above, FDM requires that finite difference stencils provide a pointwise approximation replacing the partial derivatives in the governing equation. FEM starts with a variational statement and provides elementwise approximation. Both methods discretize a continuous function (e.g., magnetic vector potential or temperature) and result in a set of simultaneous algebraic equations to be solved with respect to its nodal values. Therefore, the two methods are, in fact, quite similar.
Finite difference stencils overlap one another, and in the case of complex workpiece geometry, they could have nodes outside the boundary of the components. Finite elements do not overlap one another, do not have nodes outside the boundaries, and fit the complex shape boundary more precisely. In FEM field computation, finite elements are often introduced as a way to minimize a functional. In fact, FDM can also be described as a form of functional minimization (the socalled finite difference energy method).
Therefore, FDM and FEM are different only in the choice of mesh generation and the way in which the global set of algebraic equations is obtained. As one would expect, a comparison of the efficiency of the two methods depends on the type of problem and program organization.
As shown above, both FDM and FEM methods require a network mesh of the modeling area. Unfortunately, to suit the conditions of smoothness criteria and continuity of the governing differential equation, it is also necessary to generate a mesh/grid within electrically nonconductive areas, such as the air.
Electromagnetic field distribution in the air, in most cases of coil design and process development, can be considered useless information. Such information might be of interest only during the final design stage when evaluating electromagnetic field exposure from the induction heater. The need to always carry out computation of the electromagnetic field in the air can be considered a disadvantage of both the FDM and FEM, in particular when modeling elongated systems.
Another difficulty that appears when using FDM or FEM for electromagnetic field computation is how to treat an exterior region that extends to infinity. This deals with the infinite nature of electromagnetic wave propagation. Several methods have been used, addressing the phenomenon of an infinite exterior region. Some of those methods are the “ballooning” method, “mapping” technique, and combination of finite elements and infinite elements. However, each of the abovementioned methods has certain shortcomings.
Both the mutual impedance and BEMs do not require taking the air into consideration. This can be considered as an advantage of these methods over FDM and FEM. Since MIM and BEM require discretization of only the boundaries of the components of the induction system, it makes the mesh generation procedure relatively fast and simple, resulting in short execution time.
The MIM does not appear to be an effective computational technique for complexshaped bodies because of the known limitations of calculating selfinductances and mutual inductances.
An introduction into numerical methods used for the simulation of IH can be summarized very simply. Each of the abovedescribed methods has certain advantages. In many cases, it is effective to use a combination of methods. The right choice of simulation technique depends on the specific application and process subtleties.
Thermal modeling is usually not as cumbersome as simulation of electromagnetics. Since the boundaries of the heated parts are always well defined, both FDM and FEM are well suited to compute the temperature profiles. Because of greater flexibility of the FEM to workpiece shape variation, this method is the most popular choice for mathematical modeling of electromagnetic and heat transfer problems. Only in the case of classically shaped bodies might FDM have superior qualities over FEM. In many cases, a combination of different methods provides substantial advantages and treats application specifics better.
The majority of commercial codes used for computer modeling of IH are allpurpose programs that were originally developed for modeling processes taking place in electrical machines, motors, circuit breakers, transformers, nondestructive testing, and magnetic recording systems and were later adapted to IH needs. The need to sell their products to as many customers as possible forced early software developers to produce universal simulation tools that could be used within a broad industrial base.
Regardless of wellrecognized impressive capabilities of modern commercial software, certain process subtleties specifically related to IH might be either overlooked or substantially simplified by software developers. The result is that many generalized programs cannot address certain important features of a particular IH application. Some of the difficulties include the following:
Imagine, for instance, that you have purchased software and make an attempt to simulate two polar process stages (cold start and hot start) of the induction billet heating before forging. Cold start represents a process condition in which the induction heater was switched off for a sufficiently long time and its thermal refractory was cooled down to ambient temperature (such as after a long weekend).
In contrast, hot start designates a condition in which there was a relatively short interruption in the process cycle (e.g., because of temporarily technological issues related to short interruptions in operation of forming machines). Short interruption leads to partial refractory cooling that affects the subsequent heating process.
Suppose you know the physical properties of the refractory’s material, its thickness, and the geometry of the induction system, and, therefore, expect to be able to use purchased software to predict the effect of a cold start versus a hot start on billet thermal conditions and what would be the most appropriate process recipe that would allow you to compensate for the differences in refractory temperature.
Suddenly, you might realize that the purchased software package does not allow inputting specifics of a refractory design. The manual suggests that the user should somehow quantify the effect of refractory temperature on a billet’s thermal boundary condition taking place owing to surface heat losses (e.g., heat convection and thermal radiation). Unexpectedly, such a common design feature of any induction forge heater might become an obstacle when using generalized commercial modeling software.
Therefore, it is important to be aware that some critical feature(s) of a particular IH application could be a limiting factor for particular allpurpose software, forcing an analyst to make not welldefined assumptions that might dramatically affect the accuracy and usefulness of simulations.
Our experience shows that there is not a single universal computational method that optimally fits all induction thermal applications. Taking into consideration specifics and subtleties of a wide variety of IH processes, it is preferable to have a number of applicationoriented and highly effective software rather than searching for a single universal code. For a family of similar problems, certain numerical methods or software are preferred.
As a result, Inductoheat scientists and engineers utilize and integrate both commercial and proprietary computermodeling techniques into their professional activity. This allows them to select the technique that is most appropriate to a particular application.
Case Study. Computer modeling of induction scan hardening. As an example, Figure 3.65 shows the results of computer modeling the sequential dynamics of induction scan hardening a hollow shaft that has diameter changes and snap ring groove. Because the shaft is symmetrical, only the top half was modeled using FEA analysis. Temperature variations at four selected areas of the shaft are provided, as well as heat distribution at different inductor positions during scanning. The scan rate and coil power were varied during scanning to accommodate changes in the shaft’s geometry. Numerical computer simulations allow manufacturers of induction equipment to determine comprehensive details of the process that would be difficult, if not impossible, to determine experimentally [142].
Figure 3.65 Results of FEA computer simulation of the sequential dynamics of induction scan hardening a hollow shaft having diameter changes and groove. (From V. Rudnev, Computer modeling helps prevent failures of heat treated components, Advanced Materials & Processes, October, 2011, pp. 6–11.)
As can be noted from Figure 3.65, during scanning, considerable heating of the shaft begins appreciably in front of the copper turn, creating a preheating effect. Axial heat flow as a result of thermal conduction is one of two factors responsible for preheating. The propagation of the external magnetic field causing heat generation outside the induction coil is another factor responsible for this phenomenon.
The presence of an external magnetic field outside the induction coil is also responsible for the postheating of shaft areas located immediately behind the inductor, and, in some cases, even in regions where the spray quenchant impinges on the surface of the shaft. This can reduce quenching severity and potentially create conditions for crossing the “nose” of the continuous cooling transformation curve, resulting in the formation of mixed structures with the presence of upper transformation products (e.g., bainitic/pearlitic structures or “ghost” networking). Such microstructures are often undesirable, potentially leading to a premature failure.
The comettail effect [143,192] is clearly seen on Figure 3.65, manifesting itself as heat accumulation in subsurface regions of the shaft behind the scan inductor. It is particularly pronounced in the regions of a diameter change.
Upon quenching, the temperature of the shaft surface can be cooled sufficiently below the M _{s} temperature. At the same time, the heat accumulated in the subsurface regions might be sufficient to cause an undesirable tempering back of asquenched surface regions, leading to soft spots and undesirable properties of the heattreated workpiece.
In some cases, it is required that the hardened case should be terminated at a snap ring groove. In other cases, hardening an entire groove area may be specified. However, it is commonly highly desirable to avoid a hardness pattern termination in the middle of the groove owing to a high probability of crack development. This necessitates developing precise control algorithms that provide appropriate coil power variations during scanning.
Unfortunately, the ability to model a comettail effect and spray quenching with sufficient accuracy is limited in most commercial IH software. Besides that, some allpurpose programs cannot properly handle pre and postheating effects of scan hardening that appeared owing to the electromagnetic end and edge effects. In addition, when designing inductors and developing optimal process recipes, it is imperative to properly model not only heating but the sprayquenching stage as well. Otherwise, crucial aspects of the process might be missed.
These restrictions dramatically limit the use of some programs. Therefore, before purchasing software, make sure that it is free of critical restrictions and it is capable of properly modeling allimportant physical phenomena.
Tip #1. In some cases, the part’s geometry does not permit taking advantage of symmetry and discourages the use of 2D approaches to simulate induction thermal processes. 3D electromagnetic and thermal software is used in these cases. 3D simulation allows taking all critical geometrical features of the process into consideration. However, it is imperative to remember that any FEA computational analysis can, at best, produce only results that are derived from the properly defined theoretical model, governing equations, boundary conditions, and suitable meshing. At the end of simulation, modern 3D software does not usually provide any information regarding the accuracy of obtained results and it is up to the analyst to make a judgment on the appropriateness of the computermodeled results.
Therefore, before you hire somebody to perform computer simulations, make sure that the analyst(s) has a clear understanding of the process specifics, as well as appropriate education in the area where you are seeking help. When flying in an airplane, you need a pilot, and it is expected that a pilot has appropriate training. When you need medical assistance, it is expected that a doctor has proper education and appropriate medical degree and experience. Apply the same principles when you are choosing a company or analyst(s) to do computer simulations for you. Otherwise, you might get pretty graphics but erroneous and technically inadequate results followed by excuses as to why the analysis does not match practice. Due diligence is needed when deciding which model to apply and who should apply it. Verifying the basic education of the computer modeling analyst that you are planning to relay upon to do simulations might be the first important step in avoiding unpleasant surprises and wasting your time and money.
Experience shows that proper FEA meshing is a crucial factor affecting the accuracy of numerical simulations. Regions of high current concentration and areas where the electromagnetic field has measurable gradients must be properly meshed using a sufficient number of elements. The use of higher frequencies increases the importance of proper meshing.
Tip #2. Make sure that the physical properties of the heated materials are properly defined. Though crude, the wellknown saying “garbage in/garbage out” clearly indicates the importance of having sufficiently accurate physical properties of the heated materials. Poorly defined material properties are responsible for a significant amount of erroneous modeling results.
Tip #3. In order to minimize the risk, it is advantageous to deal with companies that offer onestop service and that are responsible for all stages of R&D (including computer modeling), design, fabrication, testing, equipment startup, and aftermarket support, rather than dealing with a number of companies with vague overall process responsibilities and potential fingerpointing if problems occur.
Tip #4. Some of today’s commercial software are not capable of accurately modeling certain features of the induction process. Many of the programs used to model IH processes are allpurpose programs originally developed to model processes occurring in electrical machines, antennas, transformers, and magnetic recording systems, and were later adapted to IH. They may be limited in their ability to take into consideration some critical features of a particular induction application. Be aware about limitations of software that are intended to be used to simulate your application.
Tip #5. In order to optimize temperature distribution before the next operation, complex control algorithms are required with varying powers, frequencies, and scan rates. The necessity of determining multiparameter control algorithms leads to long development times in the laboratory using the trialanderror method with numerous components being wasted. Computer modeling can considerably shorten the learning curve and dramatically reduce a number of parts required for trials.
Tip #6. Mathematical modeling is very beneficial not only in determining optimal process parameters and inductor design but also in evaluating the robustness of a particular heating or heattreating process by estimating, for example, the impact of reallife process deviations and magnitudes of temperatures that could potentially occur. This helps reduce the possibility of cracking, grain growth and metal waste. Reallife deviations include the following:
Results of computer modeling is helpful not only because it precludes premature failure of already designed components but also because it provides important information for component designers by taking into consideration the reallife operating requirements of the workpiece.
Tip #7. It is important to understand that the use of modern software and numerical modeling methods (including finite elements, boundary elements, finite difference, finite volume, edge elements, etc.) does not in itself guarantee the generation of perfectly correct simulation results. Rather, these techniques must be used in conjunction with experience in numerical computations and proper training as well as experience in interpreting results of modeling. This is especially so because even in commercial programs, regardless of the amount of testing and verification, software may never have all of their possible errors detected. Consequently, the analyst must guard against various kinds of potential errors. The more powerful the software, the more complex it is, increasing the potential for errors. Be aware that computergenerated attractive pictures might be misleading if a novice to the process obtains them. Common sense, engineering “gut feeling,” and advance education in the area of modeling are always the analyst’s helpful assistants.