Simulating induction heating processes using harmonic balance FEM

Abstract

Purpose

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.

Design/methodology/approach

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.

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.

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. 1562-1574. https://doi.org/10.1108/COMPEL-12-2018-0489

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⁻⁶ 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 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)∂A∂t,

for the magnetic vector potential A, where v (||∇ × A||) represents the solution dependent reluctivity, based on the commutation curve, J_i 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 E_i, included in the impressed current density later on and a solenoidal part E_s. Now the constitutive relation J = γ(E + u × B) can be rewritten as:

(4) J=Ji+γ(Es+u×B),

with J_i 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 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:

(6) Es=−∂A∂t.

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

(7) Q˙=γ∂A∂t·∂A∂t−Ji·∂A∂t .

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:

(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ωt−jωγ∑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 e^–iω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ωγ[−N⋱01⋱N][Â−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ω[−NM⋱0M⋱NM])[A¯^−N⋮A¯^0A¯^1⋮A¯^N]=[J¯^−N⋮J¯^0J¯^1⋮J¯^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, for k≠0,

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 (ℱ⁻¹) and computing the magnetic flux density:

(13) Bi(t)=∇×Ai(t),

we evaluate the commutation curve and obtain νⁱ(t). This time signal νⁱ(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 f₁ half a period, we obtain –f₁. Doing the same with f₂, results in f₂, which violates (15). On the other hand, f₃ 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˙=γ∂A∂t·∂A∂t︸Term1−Ji·∂A∂t︸Term2.

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

(21) γ∂A∂t·∂A∂t→γ(∑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+1∂A∂t·∂A∂tdt=γω2−1τi+1−τi∫τiτi+1(∑k=−NN∑l=−NNk l A^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+1∂A∂t·∂A∂tdt=12γω2∑k=−NNk2A^k·A^−k=12γω2∑k=−NNk2A^k·A^k*=12γω2∑k=−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γω2∑k=−NNk2A^k·A^k*=12γω2∑k=−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 T ∈ V: = {u ∈ H¹| u = 0, on Γ_e} such that:

(27) −∫Ωρc(T)∇T′·(uT)dx+∫Ωλ(T)∇T′·∇Tdx=−∫Γn(ρc(T)T′uTλ(T)ρc(T)+λ(T)T′∇T)·n ds+∫ΩQ˜˙ T′dx, ∀T′∈V,

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 · 10⁶ S/m and the linear permeability is chosen, based on the commutation curve, as μ₀ = 8 · 10⁻⁶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 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 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⁻⁵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 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.

Figures

Figure 1.

Alternating time frequency scheme

Figure 2.

Nonlinear magnetic thermal iteration scheme

Figure 3.

Transversal inductor above a 1-mm thin steel sheet

Figure 4.

Commutation curve, based on measurements and smooth spline approximation

Figure 5.

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

Figure 6.

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

Figure 7.

Steady-state temperature field for moving steel sheet

Figure 8.

Temperature along the centerline of the sheet

Table I.

Comparison of multiharmonic simulations with nonlinear transient ones

Considered harmonics	Multiharmonic Q˜˙MH	% Rel. error |Q˜˙MH−Q˜˙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