Simulating induction heating processes using harmonic balance FEM

Klaus Roppert (Institute of Mechanics and Mechatronics, TU Wien, Vienna, Austria)
Florian Toth (Institute of Mechanics and Mechatronics, TU Wien, Vienna, Austria)
Manfred Kaltenbacher (Institute of Mechanics and Mechatronics, TU Wien, Vienna, Austria)

COMPEL - The international journal for computation and mathematics in electrical and electronic engineering

ISSN: 0332-1649

Article publication date: 25 July 2019

Issue publication date: 21 October 2019




The purpose of this paper is to examine a solution strategy for coupled nonlinear magnetic-thermal problems and apply it to the heating process of a thin moving steel sheet. Performing efficient numerical simulations of induction heating processes becomes ever more important because of faster production development cycles, where the quasi steady-state solution of the problem plays a pivotal role.


To avoid time-consuming transient simulations, the eddy current problem is transformed into frequency domain and a harmonic balancing scheme is used to take into account the nonlinear BH-curve. The thermal problem is solved in steady-state domain, which is carried out by including a convective term to model the stationary heat transport due to the sheet velocity.


The presented solution strategy is compared to a classical nonlinear transient reference solution of the eddy current problem and shows good convergence, even for a small number of considered harmonics.


Numerical simulations of induction heating processes are necessary to fully understand certain phenomena, e.g. local overheating of areas in thin structures. With the presented approach it is possible to perform large 3D simulations without excessive computational resources by exploiting certain properties of the multiharmonic solution of the eddy current problem. Together with the use of nonconforming interfaces, the overall computational complexity of the problem can be decreased significantly.



Roppert, K., Toth, F. and Kaltenbacher, M. (2019), "Simulating induction heating processes using harmonic balance FEM", COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, Vol. 38 No. 5, pp. 1562-1574.



Emerald Publishing Limited

Copyright © 2019, Klaus Roppert, Florian Toth and Manfred Kaltenbacher.


Published by Emerald Publishing Limited. This article is published under the Creative Commons Attribution (CC BY 4.0) licence. Anyone may reproduce, distribute, translate and create derivative works of this article (for both commercial and non-commercial purposes), subject to full attribution to the original publication and authors. The full terms of this licence may be seen at

1. Introduction

Depending on the actual heating application and the thickness of the sheet, excitation frequencies of up to 100 kHz are necessary. This results in large timescale differences between the magnetic Otmagnetic) 10−6 s and the thermal field Otheat) 1 s, which makes a transient analysis not feasible. There are attempts to circumvent this issue, e.g. in Kaufmann et al. (2014), where the thermal heat conduction problem is only solved every N-th time step or by using an effective material in a harmonic simulation, e.g. in Paoli et al. (1998). Another interesting approach, using the fact, that the solution is time periodic, is presented in Biro and Preis (2006). In our approach, we use a multiharmonic ansatz to transform the nonlinear problem into frequency domain and a harmonic balancing scheme to incorporate the nonlinear commutation curve. The resulting system is solved with a nested iteration strategy [see Bachinger et al. (2002)], including an adapted algebraic multigrid (AMG) preconditioner, based on the auxiliary mesh ansatz, proposed in Reitzinger (2002).

To incorporate a moving sheet, we use the convection diffusion equation, which is the classical heat conduction partial differential equation (PDE), augmented by a convective term and solved in steady state, which further improves the computational performance of the method.

2. Definition of the coupled problem

For the electromagnetic part, we solve the eddy current problem:

(1) ×ν(||×A||)×A=Jiγ(T)At,
for the magnetic vector potential A, where v (||∇ × A||) represents the solution dependent reluctivity, based on the commutation curve, Ji the impressed current density and γ(T) the electric conductivity, which is assumed isotropic and constant, to focus on the multiharmonic approach. To obtain the temperature distribution in the moving sheet, we have to solve the nonlinear convective heat PDE:
(2) (ρc(T)T)t+·(ρc(T)Tu)=·(λ(T)T)+Q˙,
where ρ is the density, c(T) the temperature dependent specific heat capacity, λ(T) the thermal conduction coefficient, u the velocity of the sheet and Q̇ the heat source density, based on Joule losses resulting from eddy currents (no hysteretic losses are considered). These losses are defined as:
(3) Q˙=J·E,
with E as the electric field intensity. This (total) electric field intensity can then be given as the sum of any irrational part Ei, included in the impressed current density later on and a solenoidal part Es. Now the constitutive relation J = γ(E + u × B) can be rewritten as:
(4) J=Ji+γ(Es+u×B),
with Ji as the impressed current density because of a given electric potential difference (current- or voltage-loaded coil). With this relation, (3) simplifies to:
(5) Q˙(Ji+γEs)·Es.

This is valid because Esu × B for slowly moving sheets. From the definition of the magnetic vector potential (B = ∇ × A) and Faraday’s law, we can state for the solenoidal part of the electric field intensity:

(6) Es=At.

Using this in (5), results in the final form:

(7) Q˙=γAt·AtJi·At.

3. Multiharmonic ansatz

Because we are not interested in transient effects of the magnetic field, we transform the nonlinear eddy current problem (1) into frequency domain, using a multiharmonic Ansatz for A, v (||∇ × A||) and Ji of the form:

(8) f(x,t)=(k=0Nf^k(x)·ejkωt)=k=NNf^k(x)·ejkωt,
where a generic quantity f is used as a template for the above mentioned quantities. Inserting this ansatz into (1), results in:
(9) ×m=MMν^m(x)ejmωt×k=NNA^k(x)ejkωt=k=NNJ^i(x)ejkωtjωγk=NNkA^k(x)ejkωt,
where M  N. With this additional expansion number, we obtain the possibility to expand the solution Âk and the reluctivity ν̂m into different Fourier series with different truncations. This could be beneficial when considering a large number of harmonics N for the solution quantity but reducing the size of the global system by choosing a small M, which results in less off-diagonal matrices and saves memory. In this work, we set M = N, because N is not chosen as large as to gain significant benefits from a smaller M. We then multiply (9) by eiωnt for n ∈ [–N, … ,N] and integrate over one period τ = 2π/ω, to obtain the frequency domain problem. An important aspect is the existence and uniqueness of the solution for the eddy current problem, which is proven, e.g. in Bachinger et al. (2005), together with an estimate of the truncation error, because we only incorporate a finite number of harmonics into the ansatz.

Rewriting the summation as matrix-vector multiplication, results in the following system of equations:

(10) ×[ν̂0ν̂1ν̂Mν̂1ν̂0ν̂1ν̂0ν̂Mν̂Mν̂0ν̂1ν̂Mν̂1ν̂0][×A^N×A^0×A^1×ÂN]+jωγ[N01N][ÂNÂ0Â1ÂN]=[ĴNĴ0Ĵ1ĴN].

Now we can use a Galerkin ansatz together with edge finite elements in H(curl) to construct the final linear system:

(11) ([K(ν^0)K(ν^1)K(ν^M)K(ν^1)K(ν^0)K(ν^1)K(ν^0)K(ν^M)K(ν^M)K(ν^0)K(ν^1)K(ν^M)K(ν^1)K(ν^0)]+jω[NM0MNM])[A¯^NA¯^0A¯^1A¯^N]=[J¯^NJ¯^0J¯^1J¯^N],
where K(ν̂k) is the stiffness matrix associated to the k-th harmonic of the reluctivity and M the mass matrix, which is the same for all harmonics. The vectors A^¯k contain the values of the unknown vector potential for harmonic k and J¯^k the excitation current.

4. Numerical solution procedure

This global system (11) is solved via the proposed solution strategy in (Bachinger et al., 2002) with the adaption, that the Richardson solver for the inner iteration is replaced by a generalized minimal residual method (GMRES) or conjugate gradient (CG), owing to the better convergence for the application example in Section 8.

For the proposed block Jacobi preconditioner, we need a strategy to efficiently evaluate the inverse diagonal elements, which are matrices itself. Therefore an adapted AMG preconditioner is used, as presented in (Reitzinger, 2002).

An important aspect is the correct evaluation of the reluctivity Fourier coefficients ν̂k in (11), depending on the commutation curve. In this section, the subscript represents the harmonic number k ∈ [–N, N] and the superscript i, the iteration counter. Our approach consists of using an alternating time frequency scheme, as depicted in Figure 1, where we first solve the initial (linear harmonic) system with the linear reluctivity ν0ν0(t)ν^00 for the harmonics (Fourier coefficients) of the magnetic vector potential A^k0.

Because we solely inserted ν^00 into the diagonal blocks in (11) without any off-diagonal entries:

(12) ν^k0=0,fork0,
the solution A^k0 for k ∈ [–N, N] represents the linear harmonic solution for every excitation with harmonic k.

After transforming the initial solution A^k0 into time domain (ℱ−1) and computing the magnetic flux density:

(13) Bi(t)=×Ai(t),
we evaluate the commutation curve and obtain νi(t). This time signal νi(t) can then be transformed into frequency domain (ℱ) to obtain the new reluctivity Fourier coefficients ν^ki+1, used to reassemble the global system (11).

The resulting system fully incorporates the information of the commutation curve, represented by the coupling matrices (off-diagonal entries), causing a dependency between all harmonics. Solving the adapted multiharmonic system, results in an updated solution for the Fourier coefficients of the magnetic vector potential A^ki+1 for k ∈ [–N, N]. The procedure, depicted in Figure 1 is then continued until the absolute residual reaches a certain threshold.

5. Reduction to odd harmonics

The global system (11) has to be solved and reassembled several times, depending on the number of iterations, needed for convergence. This is computationally expensive, especially when increasing the number of harmonics. Therefore we are striving to increase the performance of the presented approach by considering the following aspect.

The basis for this derivation is the proof of existence and uniqueness of a solution for (10), carried out in Bachinger et al. (2005). Owing to better readability, the Fourier series in the following are considered in the real valued notation, which is equivalent to the complex valued pendant, used in (8). Now we temporarily suppose an excitation only in the base harmonic:

(14) J(t)=Jccos(ωt),
where J c is the cos-Fourier coefficient. For a generic harmonic function f (t) it is obvious, that when phase shifting the function half a period, we obtain:
(15) f(t+πω)=f(t).

Representing f(t) in a Fourier series of the form:

(16) f(t)=k=0Nfkccos(kωt)+fkssin(kωt),
it also follows from (15) that the coefficients for even harmonics (2k) are zero. This can easily be seen by comparing the following three cases:
(17) f1(t)=cos(ωt)+sin(ωt),
(18) f2(t)=cos(2ωt)+sin(2ωt)and
(19) f3(t)=cos(3ωt)+sin(3ωt).

When phase shifting f1 half a period, we obtain –f1. Doing the same with f2, results in f2, which violates (15). On the other hand, f3 fulfills (15), which shows that only odd harmonics remain, which also holds for the excitation J(t), as it fulfills (15).

This procedure only shows the property of the right hand side (excitation). In fact, also the solution has this property, which follows directly from the existence of a unique solution (Bachinger et al., 2005).

Another, more intuitive way to see this property is to consider the commutation curve ν(B), based on the actual BH-curve, where ν(B) is an even function (symmetric around B =0). Varying B with e.g. B(t) = B c cos(ωt), results in a signal with period length π/ω, which is half the period length of B(t). Therefore, no odd harmonics are present in the Fourier transform of ν(t), except the contributions from frequency zero. Inserting this relation into (11), one can observe that the rows and columns for harmonics k ∈ ℕeven are zero. Therefore also the solution Âk is zero for k ∈ ℕeven.

With this property, we can decrease the size of the global system (11) from (2 N +1) harmonics to only (N +1), which poses a significant performance improvement, especially when considering a large number of harmonics.

6. Joule losses in multiharmonic analysis

For the convective heat PDE, we need consistent heat sources for the steady state formulation, given in Section 7. Therefore we average the Joule loss density (7) over one period in frequency domain. Because the solution is given in terms of harmonics Âk for k ∈ ℕ, special care has to be taken to correctly compute the average. For better readability, both terms in (7) are handled separately:

(20) Q˙=γAt·AtTerm1Ji·AtTerm2.

Term-1. We transform the term from time to frequency domain by applying the multi-harmonic ansatz:

(21) γAt·Atγ(k=NNjkωA^kejkωt·l=NNjlωA^lejlωt).

The aim is to integrate all periodic “AC components” (alternating) parts over one period, which then vanish and only the total offset “DC component” remains. Let us now integrate (21) over one period of the base harmonic τi+1τi = 2 π/ω and concatenate the double sum:

(22) γ1τi+1τiτiτi+1At·Atdt=γω21τi+1τiτiτi+1(k=NNl=NNklA^k·A^lej(k+l)ωt)dt.

At this point we can use:

(23) τiτi+1ejhωt=0,h,
which means that all combinations A^k · A^lej(k+l)ωt vanish iff (k + l ) ≠ 0. Based on this property, we can evaluate the remaining parts (for k = –l) as:
(24) γ1τi+1τiτiτi+1At·Atdt=12γω2k=NNk2A^k·A^k=12γω2k=NNk2A^k·A^k*=12γω2k=NNk2||A^k||2,
where * denotes the conjugate complex and |·| the Euclidian norm in ℝd. The additional factor of 1/2 is valid, because in comparison to (22), where the combination k = –l occurs exactly N-times, it occurs 2 N-times in (24).

Based on the above derivation, we obtain the following expression for the first term, where the tilde represents the period averaged quantity:

(25) Q˜˙Term1=12γω2k=NNk2A^k·A^k*=12γω2k=NNk2||A^k||2

Term-2. With the assumption of an excitation only in the base harmonic, this term follows analogously to Term-1. It is not further considered in this paper, as we are not solving for the temperature distribution inside the coil.

7. Heat partial differential equation and coupling

To increase efficiency and because of the fact that we are not considering transient effects, we can reformulate (2) for the steady state case t → ∞ as:

(26) ·(ρc(T)Tu)=·(λ(T)T)+Q˜˙,

Regarding the computational domain Ω, its boundary Ω is split into Neumann (natural Γn) and Dirichlet (essential Γe) boundaries ∂Ω = Γn ∪ Γe. Then, the weak form can be formulated as:

Find TV: = {uH1| u =0, on Γe} such that:

(27) Ωρc(T)T·(uT)dx+Ωλ(T)T·Tdx=Γn(ρc(T)TuTλ(T)ρc(T)+λ(T)TT)·nds+ΩQ˜˙Tdx,TV,
where the Dirichlet boundary Γe is not visible in the weak form because the test functions T′ are (per definition of V) zero at essential boundaries. Special attention has to be paid to the convective term (first integral in the equation above) because it produces asymmetric entries in the system matrix. This can increase the condition number of the system, which is of importance, especially when using preconditioners for iterative solvers.

The global coupling of the nonlinear eddy current problem with the nonlinear thermal- (convective diffusive-) PDE is performed, according to Figure 2, where iteration counter i represents the harmonic balancing scheme from Figure 1 and iteration counter j a nonlinear Newton iteration, adapting the thermal conductivity λ(T) and specific heat capacity c(T).

8. Application example

We consider the heating of a 1 mm thick and 200 mm wide endless steel sheet with a horseshoe-shaped inductor, depicted in Figure 3, assumed to have no conductivity, to prevent eddy currents. The sheet has a constant electric conductivity of γ = 5.08 · 106 S/m and the linear permeability is chosen, based on the commutation curve, as μ0 = 8 · 10−6Vs/Am. The temperature dependency of the electric conductivity is ignored in this example because it has negligible influence on the resulting temperature distribution in the moving sheet. The considered material nonlinearity is taken into account by its commutation curve, depicted in Figure 4, which is a smooth-spline approximated curve from measurement data. It is important to notice the need for a correct approximation, to obtain a strictly monotone magnetization curve, as described in Reitzinger et al. (2002). This is essential because otherwise the unique existence of a solution for the eddy current problem cannot be ensured.

For the excitation, an impressed current in the coil, with a frequency of 2.8 kHz and a current of 9 kA is chosen. The material parameters of steel, together with the excitation frequency result in an approximate skin penetration depth of:

(28) δ=1πfγμ0=0.52mm.

The problem with this approximate skin depth is that it heavily depends on the actual working point on the commutation curve (Figure 4) and even more so on the number of considered harmonics. The above approximation is only valid for N =1, if however we increase the number of harmonics, the mesh must be significantly finer because we have to fully resolve the eddy currents up to the highest possible frequency, which is Nf. To accomplish this, the sheet is meshed with structured hexahedral elements, including a refinement to resolve the eddy currents correctly. The discretization of the air volume is performed with large, unstructured tetrahedrons, which is possible, as we are using nonconforming Nitsche interfaces. Otherwise, transition elements (pyramids) would be necessary, leading to more degrees of freedom, respectively lower computational efficiency.

With this approach, the number of elements is about 226,000 and the discretization of the continuous function space is carried out, using lowest order Nedelec edge elements.

To obtain a reference solution, we compare the multiharmonic results to a nonlinear transient analysis with a timestep size of Δt = 1 · 10−5s. For induction heating processes, the most important output quantity of the magnetic analysis are the Joule- (eddy current) losses. In total, seven periods were simulated to ensure a quasi-steady state solution. To compare these losses to the steady state ones, obtained from the multiharmonic simulation, as derived in (25), a numerical integration and averaging over three periods of the time signal of the transient analysis is carried out using the trapezoidal rule. The results for the relative errors between multiharmonic and transient analysis are given in Table I. One can observe monotone convergence for increasing numbers of considered harmonics. Additionally, the ratio of wall clock times between multiharmonic- and nonlinear transient-solution are given and a significant performance improvement can be noticed, especially for smaller numbers of considered harmonics.

In Figure 5, a quasi-steady state period of the B field from the nonlinear transient analysis is compared to multiharmonic simulations for different numbers of considered harmonics from 1 to 7. For a more quantitative comparison of nonlinear transient and multiharmonic analysis, the FFT of the transient solution is compared to the multiharmonic one for N =7 in Figure 6. In this plot, the real and imaginary parts of both, nonlinear transient and multiharmonic simulation are depicted separately and one can observe excellent matching of Fourier coefficients up to the seventh harmonic at 17.5 kHz.

For the nonlinear eddy current problem it is difficult to a priori estimate the skin depth, as it not only depends on temperature dependent electric conductivity γ(T) (which is assumed to be constant in this work) but also on the actual solution, which defines, together with the commutation curve, the reluctivity ν(||∇ × A||).

Applying the Joule- (eddy current) losses from Table I to the moving steel sheet with a velocity u=5mms, results in a temperature distribution, depicted in Figure 7. Plotting the temperature along the centerline of the sheet, Figure 8 is obtained, where one can observe a similar behavior, correlating to the error investigation in Table I. Looking at the highest temperatures in Figure 8, a large gap between the simulations with N =1 and N =3 is visible. When increasing the number of considered harmonics, this gap becomes smaller.

9. Conclusion

In this paper we presented a method to circumvent the time consuming nonlinear transient analysis of the coupled magnetic thermal problem, when simulating induction heating processes by means of a multiharmonic ansatz. For the correct handling of the solution dependent reluctivity, an alternating time frequency scheme was proposed.

This strategy was used to simulate a heating process of a thin steel sheet and the multiharmonic results (Joule losses, as well as B field) converge towards a nonlinear transient reference solution.


Alternating time frequency scheme

Figure 1.

Alternating time frequency scheme

Nonlinear magnetic thermal iteration scheme

Figure 2.

Nonlinear magnetic thermal iteration scheme

Transversal inductor above a 1-mm thin steel sheet

Figure 3.

Transversal inductor above a 1-mm thin steel sheet

Commutation curve, based on measurements and smooth spline approximation

Figure 4.

Commutation curve, based on measurements and smooth spline approximation

Comparison of magnetic flux density for transient and multiharmonic analysis in a point on the centerline of the sheet under the inductor

Figure 5.

Comparison of magnetic flux density for transient and multiharmonic analysis in a point on the centerline of the sheet under the inductor

Comparison of Fourier coefficients of transient and multiharmonic B-field for N = 7

Figure 6.

Comparison of Fourier coefficients of transient and multiharmonic B-field for N =7

Steady-state temperature field for moving steel sheet

Figure 7.

Steady-state temperature field for moving steel sheet

Temperature along the centerline of the sheet

Figure 8.

Temperature along the centerline of the sheet

Comparison of multiharmonic simulations with nonlinear transient ones

Considered harmonics Multiharmonic Q˜˙MH % Rel. error |Q˜˙MHQ˜˙trans|Q˜˙trans tMHwall/ttranswall
N =1 (1) 4.0121 · 107 8.882 0.07
N =3 (1,3) 4.3495 · 107 1.220 0.19
N =5 (1,3,5) 4.4414 · 107 0.262 0.28
N =7 (1,3,5,7) 4.4008 · 107 0.057 0.37


Bachinger, F., Kaltenbacher, M. and Reitzinger, S. (2002), “An efficient solution strategy for the HBFE method”, IGTE Proceedings, Vol. 2, pp. 385-389.

Bachinger, F., Langer, U. and Schöberl, J. (2005), “Numerical analysis of nonlinear multiharmonic eddy current problems”, Numerische Mathematik, Vol. 100 No. 4, pp. 593-616.

Biro, O. and Preis, K. (2006), “An efficient time domain method for nonlinear periodic Eddy current problems”, IEEE Transactions on Magnetics, Vol. 42 No. 4, pp. 695-698.

Kaufmann, C., Günther, M., Klagges, D., Knorrenschild, M., Richwin, M., Schöps, S. and Ter Maten, E.J. (2014), “Efficient frequency-transient co-simulation of coupled heat-electromagnetic problems”, Journal of Mathematics in Industry, Vol. 4 No. 1.

Lang, H. and Zhang, X. (2016), “The harmonic balance method”, ECE 1254 Modelling of Multiphysics Systems, Project Report.

Paoli, G., Biro, O. and Buchgraber, G. (1998), “Complex representation in nonlinear time harmonic eddy current problems”, IEEE Transactions on Magnetics, Vol. 34 No. 5, pp. 2625-2628.

Reitzinger, S. (2002), “An algebraic multigrid method for finite element discretizations with edge elements”, Numerical Linear Algebra with Applications, Vol. 9 No. 3, pp. 223-238.

Reitzinger, S., Kaltenbacher, B. and Kaltenbacher, M. (2002), “A note on the approximation of BH-curves for nonlinear magnetic field computations”, Technical Report, SFB F013: Numerical and Symbolic Scientific Computing, Linz.

Corresponding author

Klaus Roppert can be contacted at:

Related articles