Simulating induction heating processes using harmonic balance FEM
ISSN: 03321649
Article publication date: 25 July 2019
Issue publication date: 21 October 2019
Abstract
Purpose
The purpose of this paper is to examine a solution strategy for coupled nonlinear magneticthermal 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 steadystate solution of the problem plays a pivotal role.
Design/methodology/approach
To avoid timeconsuming transient simulations, the eddy current problem is transformed into frequency domain and a harmonic balancing scheme is used to take into account the nonlinear BHcurve. The thermal problem is solved in steadystate domain, which is carried out by including a convective term to model the stationary heat transport due to the sheet velocity.
Findings
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.
Originality/value
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.
Keywords
Citation
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. 15621574. https://doi.org/10.1108/COMPEL1220180489
Publisher
:Emerald Publishing Limited
Copyright © 2019, Klaus Roppert, Florian Toth and Manfred Kaltenbacher.
License
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 noncommercial purposes), subject to full attribution to the original publication and authors. The full terms of this licence may be seen at http://creativecommons.org/licences/by/4.0/legalcode
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 O(Δt_{magnetic}) ∝ 10^{−6} s and the thermal field O(Δt_{heat}) ∝ 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 Nth 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:
This is valid because E_{s} ≫ u × 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:
Using this in (5), results in the final form:
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 J_{i} of the form:
Rewriting the summation as matrixvector multiplication, results in the following system of equations:
Now we can use a Galerkin ansatz together with edge finite elements in H(curl) to construct the final linear system:
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
Because we solely inserted
After transforming the initial solution
The resulting system fully incorporates the information of the commutation curve, represented by the coupling matrices (offdiagonal 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
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:
Representing f(t) in a Fourier series of the form:
When phase shifting f_{1} half a period, we obtain –f_{1}. Doing the same with f_{2}, results in f_{2}, which violates (15). On the other hand, f_{3} 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 BHcurve, 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:
Term1. We transform the term from time to frequency domain by applying the multiharmonic ansatz:
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:
At this point we can use:
Based on the above derivation, we obtain the following expression for the first term, where the tilde represents the period averaged quantity:
Term2. With the assumption of an excitation only in the base harmonic, this term follows analogously to Term1. 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:
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 T ∈ V: = {u ∈ H^{1} u = 0, on Γ_{e}} such that:
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 horseshoeshaped inductor, depicted in Figure 3, assumed to have no conductivity, to prevent eddy currents. The sheet has a constant electric conductivity of γ = 5.08 · 10^{6} S/m and the linear permeability is chosen, based on the commutation curve, as μ_{0} = 8 · 10^{−6}Vs/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 smoothspline 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:
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 N ⋅ f. 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^{−5}s. 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 quasisteady 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 transientsolution are given and a significant performance improvement can be noticed, especially for smaller numbers of considered harmonics.
In Figure 5, a quasisteady 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
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.
Figures
Comparison of multiharmonic simulations with nonlinear transient ones
Considered harmonics  Multiharmonic

% Rel. error



N = 1 (1)  4.0121 · 10^{7}  8.882  0.07 
N = 3 (1,3)  4.3495 · 10^{7}  1.220  0.19 
N = 5 (1,3,5)  4.4414 · 10^{7}  0.262  0.28 
N = 7 (1,3,5,7)  4.4008 · 10^{7}  0.057  0.37 
References
Bachinger, F., Kaltenbacher, M. and Reitzinger, S. (2002), “An efficient solution strategy for the HBFE method”, IGTE Proceedings, Vol. 2, pp. 385389.
Bachinger, F., Langer, U. and Schöberl, J. (2005), “Numerical analysis of nonlinear multiharmonic eddy current problems”, Numerische Mathematik, Vol. 100 No. 4, pp. 593616.
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. 695698.
Kaufmann, C., Günther, M., Klagges, D., Knorrenschild, M., Richwin, M., Schöps, S. and Ter Maten, E.J. (2014), “Efficient frequencytransient cosimulation of coupled heatelectromagnetic 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. 26252628.
Reitzinger, S. (2002), “An algebraic multigrid method for finite element discretizations with edge elements”, Numerical Linear Algebra with Applications, Vol. 9 No. 3, pp. 223238.
Reitzinger, S., Kaltenbacher, B. and Kaltenbacher, M. (2002), “A note on the approximation of BHcurves for nonlinear magnetic field computations”, Technical Report, SFB F013: Numerical and Symbolic Scientific Computing, Linz.