Induction Heating Computer Simulation

Induction heating can be simulated very precisely by computer technology. Simulations help to better understand the physical processes of induction heating and provide valuable information on quantities that are often difficult to measure experimentally. They allow us to look inside the part and see where heat is generated and where it is dissipated. By simulation, unlike the experiment, the whole dependence on input variables can be evaluated very quickly, which is useful for design and optimization. The most common tasks are: determination of the time needed to heat up the workpiece, calculation of the core-surface temperature profile, optimization of the coil shape, determination of heating efficiency. The simulations help us to find the best solution when designing a new heater, but also to find out the optimal heating process settings in operation. We use different types of simulation programs according to the complexity of the task. To simulate our heaters we develop our own simulation software.

Send us your heating simulation request by e-mail indukce@roboterm.cz or using the inquiry form:

Induction heating simulation examples

Square billet induction heating simulation
Shaped sheet metal induction heating simulation
Simulation of induction shrink fitting of bearings into housing

How the simulation works

From the simulation engineer's point of view, the process can be divided into three phases: preprocessing, calculation and post-processing. In the preprocessing phase, a geometric model is created. Geometric model is divided by mesh into small elements (usually triangles in 2D problem, tetrahedron in 3D problem). Material properties (e.g. thermal conductivity), initial conditions (e.g. initial temperature) and boundary conditions (e.g. coil current, radiation heat loss) are assigned to model areas. A suitable calculation model is selected and its parameters are set (e.g. heating time, time step). Induction heating simulation is a coupled problem in which the electromagnetic field is solved by harmonic analysis and the temperature field as a transient process where the generated heat is equal to Joule's losses in the material. The computational program describes each element by differential equations (Maxwell's equations, heat conduction equations) and solves the system of equations numerically in each time step. The output of the solver is a time-dependent electromagnetic and temperature field. The result is finally modified by the simulation engineer in a postprocessor into a suitable graphical form (graph, animation).

Accuracy of the simulation

A frequently asked question is how close the simulation result is to reality. The accuracy of the result is influenced by several factors: choice of computational model, simplification of geometry, mesh density, time step and accuracy of material properties, boundary and initial conditions. The simulation engineer must use a suitable calculation model. If the magnetic steel is heated to the forging temperature, software with a non-linear solver must be used which takes into account a sharp drop in permeability around the Curie temperature. The disadvantage of using a nonlinear solver is a significant increase in computational time. With the use of a nonlinear solver in the case of a large number of elements (3D problems, more complex 2D problems) the limits of common computer technology are soon reached. In order to speed up the calculation, it is appropriate to simplify the task. In many cases, 3D tasks can be simplified to 2D tasks. In some cases, such as coil shape optimization, it is not necessary to calculate the entire transient process, but alone harmonic analysis is sufficient (temperatures are not calculated but only Joule losses). The simulation engineer must decide what simplifications can be accepted without significantly affecting the result. Care must be taken to design the mesh of elements. The mesh must be sufficiently fine, but not too much, to not prolong computing time unnecessarily. Therefore, the mesh is designed with different densities that are highest at the places of large changes of calculated quantities. A particularly fine mesh must be designed in the surface layer of the workpiece when the penetration depth is very small, such as in the case of steel in a magnetic state. If a suitable simulation model is created, the accuracy of the result depends only on the accuracy of the input data (material properties, boundary and initial conditions). It is the lack of knowledge of the accurate input data that is the main cause of the often discussed disagreement of the simulation result with reality. Since it is expensive to find out accurate input data (mainly material properties and emissivity over the entire temperature range), an estimate is usually entered in the calculation. Accurate absolute results can be achieved after tuning up the input data according to the experiment. Often dependencies, but not absolute values, are essential to the conclusion. These can be evaluated very reliably using simulations without knowing the accurate input data.

Dependence of material properties on temperature

Material properties may strongly depend on temperature. If the result of the simulation should be accurate, the simulation must take this into account. For steels, the temperature characteristics strongly depend on their chemical composition. The graphs below show how the material properties of construction steels with a carbon content of approximately 0.2% may depend on temperature. The graphs show a phase change around 760°C (Currie temperature). At this temperature the ferromagnetism is lost and the relative permeability drops sharply to 1. The permeability of ferromagnetic materials is complicated to include in calculations, because it depends on both temperature and magnetic field intensity. The permeability graph here shows only its drop. Absolute values would be obtained from the magnetization curve of the material. Another complicated physical property is emissivity, which indicates the amount of heat loss emitted by the surface. Emissivity depends on surface structure. The steel surface oxidizes at high temperatures, which increases emissivity.


Fig. 1 - specific heat capacity	Fig. 2 – thermal conductivity

Fig. 3 - resistivity	Fig. 4 – relative permeability

Theory of 1D induction heating

For simplicity, we describe the axisymmetric 1D problem that corresponds to the static induction heating of a long cylindrical billet in a long coil. Material properties are assumed to be constant. We will therefore solve the physical quantities only in the radial direction. Let's start from Ampere (1) and Faraday (2) law.

$\oint_{C}\vec{H}\cdot \delta \vec{l}= I$ (1)

$\oint_{C}\vec{E}\cdot \delta \vec{l}=-\frac{\partial \Phi }{\partial t}$ (2)

where:

H [A/m] is magnetic field intensity. (time varying vector, parallel to coil axis).
I [A] is the free current flowing through any area with boundary C. Bound currents are neglected.
E [V/m] is the electric field intensity (time varying vector, tangential direction).
$\partial \Phi /\partial t$ [Wb/s] is the change in magnetic flux per time.

Fig. 1 – cylindrical workpiece undergoing induction heating, source [1]

Fig. 1 – cylindrical workpiece undergoing induction heating, source [1]

Integrating equation (1) around the rectangle shown in Fig. 1 and replacing the current by the product of current density J [A / m2] and the area of the indicated rectangle yields:

$(H_z+\delta H_z)l + H_z l=J_\phi \delta r l$ (3)

And therefore the current density in the phi direction at depth r is given by:

$J_\phi=-\frac{\partial H_z}{\partial r}$ (4)

We assume a homogeneous material with relative permeability $\mu_r$ and resistivity $\rho [\Omega m]$ . For resistivity, the Ohm's law in differential form is:

$\vec{E}=\rho \vec{J}$ (5)

For magnetic flux through the surface S oriented perpendicular to the magnetic flux:

$\Phi =\mu_0 \mu_r HS$ (6)

Integrating equation (2) around the circle shown in Fig.1. and substituting relations (5) and (6) yields:

$2\pi (r+\delta r)\rho(J_\phi+\delta J_\phi)-2\pi r\rho J_\phi=-\mu_0 \mu_r 2 \pi r\frac{\partial H_z}{\partial t}\partial r$ (7)

A negligibly small term $2\pi r\rho J_\phi$ can be eliminated and equation simplified to:

$\frac{\partial J_\phi}{\partial r}+\frac{J_\phi}{r}+\frac{\mu_0 \mu_r}{\rho} \frac{\partial H_z}{\partial t}=0$ (8)

Combining equations (4) and (8) yields the 1D magnetic field equation of a long coil:

$\frac{\partial^2 H_z}{\partial r^2}+\frac{1}{r} \frac{\partial H_z}{\partial r}-\frac{\mu_0 \mu_r}{\rho} \frac{\partial H_z}{\partial t}=0$ (9)

An analytical solution of this equation is given in [1]. If the magnetic field strength distribution in the blank is solved, the current density can be calculated using equation (4). The heat generation in the material is then given by Joule's loss. Hysteresis losses can be neglected in most induction heating applications. Losses density p [W / m3] is calculated:

$p=\rho J^2$ (10)

The following graphs show the solved physical quantities in the heated workpiece:


Fig. 2 – Magnetic field intensity	Fig. 3 – Current density	Fig. 4 – Losses density

The temperature inside the workpiece can be solved using the 1D equation for transient heat conduction, which has the following form for radial direction:

$\frac{\partial T}{\partial t} = \frac{1}{r}\frac{\partial }{\partial r}\left ( \frac{\lambda r}{\rho c} \frac{\partial T}{\partial r} \right )+\frac{p}{\rho c}$ (11)

where:
$T$ [°C] is temperature
$t$ [s] is time
$\lambda$ [W/mK] is thermal conductivity
$\rho$ [kg/m³] is density
$c$ [J/kgK] is specific heat capacity
$p$ [W/m³] is losses density

Radiation losses are included in the temperature calculation as a boundary condition for the surface. Radiation losses are proportional to the fourth power of the surface thermodynamic temperature and cannot be neglected when heated to high temperatures. The surface heat flux density is calculated using the Stefan-Boltzmann law:

$j=\varepsilon \sigma \left ({T_{R}}^{4}-{T_{0}}^{4} \right )$ (12)

where:
$j$ [W/m²] is the surface heat flux density
$\varepsilon$ is the mutual emissivity of the workpiece and its surroundings (refractory)
$\sigma$ is Stefan-Boltzmann constant, its value is 5,67×10^-8 Wm^-2K^-4 $T_{R}$ [K] is the thermodynamic temperature of the surface of the heater
$T_{0}$ [K] is the thermodynamic temperature of the surrounding (refractory)

Possible temperature distribution in the workpiece at the end of heating:

Fig. 7 – core-surface temperature profile

References:
[1] M. W. Kennedy, "Magnetic Fields and Induced Power in the Induction Heating of Aluminium Billets", 2013, ISBN 9789175018102
[2] E. Rapoport, Y. Pleshivtseva, "Optimal Control of Induction Heating Processes", 2006, ISBN 9780849337543