© B. Meriem et al., Published by EDP Sciences, 2022

1 Introduction

Aluminum profiles are widely used products in today’s world. For the manufacturing of these products, electromagnetic induction heaters are used at aluminum plants to preheat aluminum cylinders, so called billets, before extruded to profiles. In these plants, typical billets having 200 mm diameter and 1000 mm height are heated up to a temperature level between 400 °C and 500 °C (depends on extrusion process) in order to soften the metal before it is pressed through the extruder. It is important to reach a temperature distribution as uniform as possible in the billet, for the quality of the forming process [1].

Electromagnetic induction heating is a fast and clean possibility to heat. The process can be controlled precisely and reproducible. Since energy costs are growing, efficiency of heating processes becomes a more and more important issue. Heating of nonmagnetic metals with low electrical resistivity like aluminum is connected with a low efficiency which is about 50%, if done in conventional way [2,3]. In order to save energy, a new approach to improve the efficiency is to generate a DC magnetic field by using an inductor based on modern permanent magnets. The electromagnetic field is applied in that manner that by rotation of the billet or the inductor, currents are induced in the workpiece. Due to Joule losses, the billet is heated up. A static magnetic field is generated without losses by permanent magnets, so that the efficiency of the system is determined by the efficiency of the main drive rotating the billet or the inductor. An efficiency of 85% is estimated for the process [4].

This work proposes to analyze the device of induction heating of nonmagnetic conductive billet illustrated in Figure 1. In this structure, the inductor is based on permanent magnets (Nd-Fe-B) with radial magnetization, bonded to the inner surface of a ferromagnetic yoke. In this topology, aluminum billet to be heated is placed inside the static magnetic field created by permanent magnets. The rotational movement of the billet by using a mechanical drive system produces the eddy currents and the heat power losses in the billet. However, physical models are needed for the design and analysis of these heaters, which consists to establish a mathematical formulation describing both electromagnetic and thermal phenomena. So, in this work, a new coupling way between these two phenomena is proposed. First, an exact 2D analytical solution of the electromagnetic field is derived from Maxwell equations in time harmonic formulation thanks to the separation of variables method. The resulting electromagnetic field is then used to calculate the heat power density, which represents the source of thermal model. Then, the magneto-thermal problem is obtained by coupling the developed magnetic model to a 2D transient thermal model based on the lumped thermal network model taking into account the temperature dependence of electrical and thermal physical properties of the heated billet. The results issued from the developed models are compared with those obtained by finite element method. The principal physical and geometrical parameters of the studied heater are given in Table 1.

Fig. 1 2D induction heater geometry.

2 2D electromagnetic analytical model

Using Maxwell’s equations in term of the magnetic vector potential $\overrightarrow{A}$ and considering that the permanent magnets inductor is fixed and the billet is suggested to a rotational movement with constant speed Ω, the electromagnetic behavior of the studied device is governed by:

$$\overrightarrow{rot}\left(\frac{1}{\mu }\overrightarrow{rot}\overrightarrow{A}\right)=\sigma (\overrightarrow{v}\wedge \overrightarrow{rot}\overrightarrow{A})+\overrightarrow{rot}\left(\frac{\overrightarrow{M}}{\mu }\right)$$(1)

where $\overrightarrow{v}$ represents the tangential velocity of the billet, related to the angular speed by:

$$\overrightarrow{v}=r\cdot \mathrm{\Omega }\cdot \overrightarrow{u_{\theta }}$$(2)

For a 2D formulation, the analytical resolution of (1) can be derived by considering the following assumptions:

• the device is considered with infinite length, and then the axial edge effects are neglected.

• the permeability of the iron yoke is large enough to be considered as infinite (the saturation is not considered).

The polar coordinates (r,θ) are appropriate to this study, according to the geometry of the device. So, (1) is rewritten as:

$$\frac{{\partial }^{2}A}{\partial r^2}+\frac{1}{r}\frac{\partial A}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}A}{\partial {\theta }^2}=\mu \sigma \mathrm{\text{Ω}}\frac{\partial A}{\partial \theta }+\frac{1}{r}\left(\frac{\partial M_r}{\partial \theta }\right)$$(3)

where M _r denotes the radial component of the magnets magnetization and A(r,θ) is the axial magnetic vector potential component.

To solve (3), an analytical method based on the separation variables technique is applied [5].

2.1 Analytical solution of electromagnetic problem

Under the assumptions cited previously, the study area in polar coordinates (r,θ) is reduced to three regions, as shown in Figure 2 and (3) becomes:

• region (1): nonmagnetic conductive billet (0 ≤ r ≤ R₁)

  $$\frac{{\partial }^2{A}_1}{\partial r^2}+\frac{1}{r}\frac{\partial A_1}{\partial r}+\frac{1}{r^2}\frac{{\partial }^2{A}_1}{\partial {\theta }^2}=\mu \sigma \mathrm{\Omega }\frac{\partial A_1}{\partial \theta }$$(4)

• region (2): air gap (R₁ ≤ r ≤ R₂)

  $$\frac{{\partial }^2{A}_2}{\partial r^2}+\frac{1}{r}\frac{\partial A_2}{\partial r}+\frac{1}{r^2}\frac{{\partial }^2{A}_2}{\partial {\theta }^2}=0$$(5)

• region (3): permanent magnets (R₂ ≤ r ≤ R₃)

  $$\frac{{\partial }^2{A}_3}{\partial r^2}+\frac{1}{r}\frac{\partial A_3}{\partial r}+\frac{1}{r^2}\frac{{\partial }^2{A}_3}{\partial {\theta }^2}=\frac{1}{r}\frac{\partial M_r}{\partial \theta }$$(6)

The following boundary and interface conditions must be added to the above equations:

In region (2), the permanent magnets magnetization constituting the field source is expressed by its Fourier series expansion [6,7] given by:

$$M_r\left(\theta \right)={\displaystyle \sum }_{n=1,3,\mathrm{5..}}^{\infty }{M}_n\cos \left(np\theta \right)$$(8)

where p is the pole pair number, M _n is the magnitude of the nth harmonic given by:

$$M_n=\frac{4M_0}{n\pi }\sin \left(np\frac{\beta }{2}\right)$$(9)

In (9), Β is the opening angle of the permanent magnets.

In the complex representation and using the separation variables technique, the general solution form of the vector potential, in each region, can be written as:

$$\underset{\_}{A}i\left((r,\theta )\right)={\displaystyle \sum }_n{\underset{\_}{A}}_{in}\left(r\right)e^{jnp\theta }\left(i=1,2,3\right)$$(10)

where j ² = −1 denotes the complex number.

By injecting (10) in (4), (5) and (6) respectively the corresponding ${\underset{\_}{A}}_{in}(r)$ of (10) is then given as:

• region (1):

  $${\underset{\_}{A}}_{1n}\left(r\right)={\underset{\_}{B}}_{1n}{J}_{np}\left(j^{\frac{3}{2}}{k}_{n}r\right)$$(11)

  with:

  $$k_n=\sqrt{np{\mu }_0\sigma \mathrm{\text{Ω}}}$$(12)

  and $J_{np}\left(j^{\frac{3}{2}}{k}_{n}r-\right)$ is the Bessel function of the first kind and npth order.

• region (2):

  $${\underset{\_}{A}}_{2n}\left(r\right)={\underset{\_}{B}}_{2n}{r}^{np}+{\underset{\_}{C}}_{2n}{r}^{-np }$$(13)

• region (3):

  $${\underset{\_}{A}}_{3n}\left(r\right)={\underset{\_}{B}}_{3n}{r}^{np}+{\underset{\_}{C}}_{3n}{r}^{-np}+{\underset{\_}{\psi }}_n\left(r\right)$$(14) $${\underset{\_}{\psi }}_n\left(r\right)=\left\{\begin{array}{ll}j\frac{M_{n }np}{1-{\left(np\right)}^2}r\hfill & if\left(np\ne 1\right)\hfill \\ j\frac{M_{n }}{2}r \left(\ln \left(r\right)-\frac{1}{2}\right)\hfill & if\left(np=1\right)\hfill \end{array}\right)$$(15)

The integration constants ${\underset{\_}{B}}_{1n}$, ${\underset{\_}{B}}_{2n}$, ${\underset{\_}{C}}_{2n}$, ${\underset{\_}{B}}_{3n}$ and ${\underset{\_}{C}}_{3n}$ are calculated by using the interface and boundary conditions (7).

Fig. 2 Geometrical model of the electromagnetic problem.

(7)

2.2 Induced current density and power heating expressions

The current density induced in the billet (region 1) by the rotation movement, is given by:

$$\overrightarrow{J_1}=\sigma \left(\overrightarrow{v}\wedge \overrightarrow{B_1}\right)=\sigma \left(\overrightarrow{v}\wedge \overrightarrow{rot}{\overrightarrow{A}}_1\right)$$(16)

By considering (3), (10) and (11), the complex form of the current density (16) can be written as:

$${\underset{\_}{J}}_1\left(r,\theta \right)={\displaystyle \sum }_n{\underset{\_}{J}}_{1n}\left(r\right)e^{jnp\theta }$$(17)

For the nth harmonic ${\underset{\_}{J}}_{1n}(r)$ is given by:

$${\underset{\_}{J}}_{1n}\left(r\right)=-jnp\sigma \mathrm{\text{Ω}}{\underset{\_}{B}}_{1n}{J}_{np}\left(j^{\frac{3}{2}}{k}_{n}r\right)$$(18)

The corresponding thermal source allowing to heat the billet is the induced power density, and is related to the induced currents by [5]:

$$p_1\left(r,\theta \right)=\frac{J_1^2\left(r,\theta \right)}{\sigma }$$(19)

Hence, the analytical expression of the power heating density can be obtained by the following relation:

$$p_1\left(r,\theta \right)=\sigma {\left[\text{Re}\left\{{\displaystyle \sum }_n-jnp\mathrm{\text{Ω}}{\underset{\_}{B}}_{1n}{J}_{np}\left(j^{\frac{3}{2}}{k}_{n}r\right)e^{jnp\theta }\right\}\right]}^2$$(20)

An other important characteristic of the studied heater device is the corresponding total heating power in the billet. It is determined by integration of (20) over the billet volume:

$$P_T=\underset{0 0 0}{\overset{l R_{1 }2\pi }{{\displaystyle \int \int \int }}}\frac{J_1^2\left(r,\theta \right)}{\sigma }rdrd\theta dz$$(21)

After integration, the final expression is:

where $a=-jnp\sigma \mathrm{\text{Ω}}{\underset{\_}{B}}_{1n}$ and $a^{*}=conj\left(a\right)$; $K=\sqrt{2}{k}_n\left(\frac{j-1}{2}\right)$ and $K*=conj\left(K\right)$

(22)

3 2D thermal transient model

The spatiotemporal distribution of temperature inside the conductive billet is governed by the following partial differential equation [8,9]:

$$\rho C_p\left(\frac{\partial T}{\partial t}\right)=div\left(\lambda \overrightarrow{grad}\left(T\right)\right)+p_1\left(r,\theta \right)$$(23)

representing the transient heat diffusion equation, which is solved by using the equivalent lumped thermal network model [10–12].

The thermal source, $p_1\left(r,\theta \right)$, representing the power density induced in the billet given by (20) is the term allowing coupling the both electromagnetic and thermal problems. By considering a weak coupling between the two models, the electrical and thermal properties of the billet are updated for each variation of the temperature. In this analysis, the variation of these physical properties as a function of the temperature, is given by the following expressions:

• for electrical conductivity:

  $$\sigma \left(T\right)=\frac{{\sigma }_0}{1+\alpha \left(T-T_0\right)}$$(24)

  where ${\sigma }_0$ is the conductivity at T ₀ (°C), $\alpha$ is the temperature coefficient relative to electrical conductivity.

• for thermal conductivity:

  $$\lambda \left(T\right)=\sigma L_{0}T$$(25)

where L ₀ = 2.45 * 10⁻⁸ is the constant of Lorentz and $\sigma$ is the electrical conductivity.

In this model, the heat transfer by conduction is considered in the radial and axial (r, z) directions.

Applied to the studied PMIHD, the proposed thermal modeling method consists in establishing a lumped thermal network where the thermal limits problem comes down only to the billet to be heated and the air layer separating the billet and the inductor. This air layer is represented in the thermal model by a convective thermal resistance [13–15].

The discretization of the studied geometry in the two directions leads to elementary volumes of cylindrical shapes, that are considered as thermal nodes. The equation of the heat exchanges between the nodes gives a set of differential equations in matrix form. For each node, the associated heat transfer equations are in the form :

$${\rho }_i{V}_i{C}_{pi}\frac{dT}{dt}={\displaystyle \sum }_{j=1}^n{G}_{ji}\left(T_j-T_i\right)+P_i$$(26)

where ${\rho }_i$, $V_i$, $C_{pi}$, $P_i$ and n characterising the block ‘i’ are respectively the mass density, the elementary volume, the specific heat, the heat source and the number of nodes linked to this block.

3.1 Lumped thermal network of an elementary block

Figure 3 shows the geometrical dimensions and considered boundary conditions of a hollow cylinder having a uniform volumic heat source p ₁ and a thermal conductivity λ.

For a study in the (r, z) plane, one considers that:

• the temperature is invariant as a function of the angle θ; hence the heat flow presents two components: axial and radial.

• the block has the same average temperature in the axial and radial directions.

• the thermal conductivity in the axial and radial directions is the same (isotropic).

In steady state, the heat transfer equation in this block is derived in cylindrical coordinates from (23):

$$\frac{{\partial }^{2}T}{\partial r^2}+ \frac{1}{r}\frac{\partial T}{\partial r}+ \frac{{\partial }^{2}T}{\partial z^2}+\frac{p_1}{\lambda }= 0$$(27)

Fig. 3 Bidirectional block.

3.1.1 Lumped thermal network of an elementary block in the axial direction

The lumped thermal network in the axial direction of each block element is obtained by solving (27) analytically only in the axial direction. Its solution is given by:

$$T\left(z\right)=- \frac{p_1}{2\lambda }{z}^2+\left(\frac{T_4- T_3}{l} +\frac{p_1.l}{2\lambda }\right)z+T_3$$(28)

The average temperature is given by:

$$T_m= \frac{1}{l}\underset{0}{\overset{l}{{\displaystyle \int }}}T\left(z\right)dz= \frac{P{R}_z}{12}+ \frac{( T_3+ T_4 )}{2}$$(29)

where $P=p_1.l.S$ and $R_z=\frac{l}{\lambda .\pi \left(r_2^2-r_1^2\right)}$ represent respectively the heat source and the block axial resistance.

The translation of (29), in the form of a network, makes it possible to represent the block by the networks of Figure 4 where $R_{z1}=R_{z2}=\frac{R_z}{2}$ and $R_{zm}=-\frac{R_z}{6}$.

In the case where the geometry presents a symmetry, so that the temperature field is symmetrical compared to the median plane (T ₃ = T ₄), the previous lumped thermal network is simplified (Fig. 5) to a single resistance $R_{zam}=\frac{R_z}{6}$.

Fig. 4 Equivalent lumped thermal network in the axial direction: (a) $P=0$, (b) $P\ne 0$.

Fig. 5 Equivalent thermal network considering symmetry.

3.1.2 Lumped thermal network of an elementary block in the radial direction

In the radial direction the lumped thermal network of each block element is obtained analytically by solving (27) in the radial direction. Its solution is given by:

$$T\left(r\right)=- \frac{p_1}{4\lambda }{r}^2+a ln\left(r\right)+b$$(30)

with α and β are constants determined as a function of the temperature imposed on the boundaries. The average temperature at the center of the block is given by:

$$T_m=P\left(R_{rm}+\frac{R_{r1}{R}_{r2}}{R_{r1}+R_{r2}}\right)+T_1\frac{R_{r2}}{R_{r1}+R_{r2}}+T_2\frac{R_{r1}}{R_{r1}+R_{r2}}$$(31)

with:

$$P=2\pi \left( r_2{}^2- r_1{}^2\right)l{p}_1$$ $$R_{r1}= \frac{1}{4\pi l\lambda }\left(\frac{2r_2^{2}ln\left(\frac{r_2}{r_1}\right)}{\left(r_2^2-r_1^2\right)}-1\right)$$ $$R_{r2}= \frac{1}{4\pi l\lambda }\left(1- \frac{2r_1^2\mathit{\text{ln}}\mathit{}\left(\frac{r_2}{r_1}\right)}{\left(r_2^2-r_1^2\right)}\right)$$

and$$R_{rm}= \frac{-1}{8\pi \lambda l\left( r_2{}^2- r_1{}^2\right)}\left(\left( r_2{}^2+ r_1{}^2\right)- \frac{4r_1^2{r}_2{}^2\mathit{\text{ln}}\left(\frac{r_2}{r_1}\right)}{\left(r_2^2-r_1^2\right)}\right)$$

The corresponding thermal networks to model the block radially are given in Figure 6 for the two cases, with and without sources.

Fig. 6 Equivalent lumped thermal networkin the radial direction: (a) $P=0$, (b) $P\ne 0$.

3.2 Lumped thermal network of the overall elementary block

The correspending thermal network of the block is obtained by combining the two radial and axial thermal networks with symmetry along z as shown in Figure 7.

The change in temperature over time (transient mode) is taken into account by defining a thermal capacity connected to the central node. It is evaluated by the relation:

$$C_i={\rho }_i{\text{V}}_i{C}_{pi}$$(32)

If the cylinder is full, it is enough to put:

$$R_{rm}^{\prime }= \underset{r_1 0}{\lim }{R}_{rm},R_{r1}^{\prime }= \underset{r_{1}0}{\lim }{R}_{r1}$$(33)

Fig. 7 Equivalent thermal network of elementary block with symmetry.

3.3 Global thermal network

The global thermal network is obtained by connecting all elementary networks through the nodes of connection. The border nodes are connected to the air gap through convection thermal resistances, R _cz and R _cr respectively in the axial and radial directions. These convection resistances are given by:

$$R_c=\frac{1}{h\cdot S}$$(34)

where S is the exchange surface and h is the thermal convective coefficient.

For a discretization of one block along z and n blocks along r the equivalent thermal network of the heating billet is represented in Figure 8 and the correspending differential equation (26) for the n radial nodes leads to, in matrix form:

$$\left\{\begin{array}{c}\left[C\right]\left\{\frac{d{T}_m}{dt}\right\}+\left[G_{11}\right]\left\{T_m\right\}+\left[G_{12}\right]\left\{T\right\}=\left\{P_t\right\}\\ \left[G_{21}\right]\left\{T_m\right\}+\left[G_{22}\right]\left\{T\right\}=\left\{Q\right\}\end{array}\right\}$$(35)

$\left[C\right]$ diagonal thermal capacity matrix of $\left(n\times n\right)$; $\left[G_{11}\right]$ thermal conductances matrix of $\left(n\times n\right)$; $\left[G_{12}\right]$ thermal conductances matrix of $n\times \left(n-1\right)$; $\left[G_{21}\right]$ matrix transposed from$\left[G_{12}\right]$; $\left[G_{22}\right]$ thermal conductances matrix of $\left(n-1\right)\times \left(n-1\right)$; $\left\{T_m\right\}$ average temperature vector; $\left\{T\right\}$ temperature vector of connecting nodes; $\left\{P_t\right\}$ heat sources vector; $\left\{Q\right\}$ vector containing the terms T ₀ for the connecting nodes; T ₀ is the initial temperature of billet (T ₀ = 30 °C).

One note, that the heat sources are evaluated using the analytical expression of the power, obtained by the electromagnetic model developed previously.

Finally the differential equation system obtained from (35) is then:

$$\left[C\right]\left\{\frac{d{T}_m}{dt}\right\}+\left[G\right]\left\{T_m\right\}=\left\{F\right\}$$(36)

with:

$$\left\{\begin{array}{c}\left[G\right]=\left[G_{11}\right]-\left[G_{12}\right]*{\left[G_{22}\right]}^{-1}*\left[G_{21}\right]\\ \left[F\right]=\left[P_t\right]- [G_{12}]*{[G_{22}]}^{-1}*\left[Q\right]\end{array}\right\}$$(37)

Fig. 8 Global equivalent thermal network (n blocks along r and one block along z).

4 Results of the 2D hybrid magneto-thermal problem

The electromagnetic quantities obtained with the developed analytical model are shown in Figures 9 and 10. The eddy current density induced in the aluminum billet by movement is illustrated in Figure 9 at 1500 rpm. So, the corresponding total power heat induced in the billet versus angular speed is shown in Figure 10. Obviously, as shown by the analytical expression (22) and Figure 10, the total heating power is directly related to the square of the speed. Indeed, the total heating power is higher all the more so as the speed is higher. One note that the total power heat obtained at 1500 rpm is about 850 W.

The transient temperature obtained by the lumped thermal network model at the billet surface at 1500 rpm is shown in Figure 11 for one block in the axial direction and n = 54 blocks in the radial directions.

In metallurgical industry, the extrusion or shaping of conductive billets requires that they be heated to temperature levels of 400–500 °C. The heating time required to ensure good productivity of the heating device is of the order of a hundred seconds. So, one can see from these results, that the required heating levels are reached after 450’s (either 8 min) for Ω = 1500 rpm.

The results obtained by this 2D hybrid magneto-thermal model are in good agreement with those obtained by the 2D finite element method.

Figure 12 illustrates the temperature distribution in the radial direction obtained by the lumped thermal network model at a rotation speed of 1500 rpm and for a heating time of 450’s. It is observed from this figure that the temperature of the block located in the centre of the billet is T _m = 411.3 °C and the temperature of the block in its surface is T _m = 413.7 °C; The correspending difference is ΔT _m = 2.4 °C (0.58%).

In practice, the forming operation can only be successful if the temperature profile through the volume of the billet is as homogeneous as possible. A temperature uniformity between 2.5% and 5% is required. In this case, one can then conclude that the temperature profile obtained through the volume of the billet is homogeneous.

Fig. 9 Current density induced in billet vs r(Ω = 1500 rpm, θ = 0°).

Fig. 10 Total power heat induced in billet vs. rotational speed.

Fig. 11 Temperature evolution at the surface of the billet (r = 38 mm, Ω = 1500 rpm).

Fig. 12 Temperature evolution in billet vs r (θ = 0°, Ω = 1500 rpm, t = 450 s).

5 Conclusion

2D hybrid model with new appraoch to the magneto-thermal coupling problems was developed in this work. This model was applied to a permanent magnet induction heating device of aluminum billet, which combines the static magnetic field with the rotational movement at constant speed. The magnetic and thermal results, obtained using the developed models are in good agreement with those obtained by the finite element method.

The obtained spatiotemporal distribution of the temperature in the billet shows that a homogeneous temperature profile can be reached with this induction heating device.

As a consequence, the proposed new 2D hybrid model can potentially be used to speed up magneto-thermal design of this type of induction heating devices.

Author contribution statement

Meriem Boubezari, Youcef Ouazir and Ammar Abdi: lumped parameter thermal model, programming and numerical simulation, and writing a part of the manuscript. Georges Barakat and Yacine Amara: analytical magnetic model, programming and numerical simulation, and writing a part of the manuscript. All authors contribute to the results discussion part. All authors have read and approved the final manuscript.

References

All Tables

All Figures

	Fig. 1 2D induction heater geometry.
In the text

	Fig. 2 Geometrical model of the electromagnetic problem.
In the text

	(7)
In the text

	(22)
In the text

	Fig. 3 Bidirectional block.
In the text

	Fig. 4 Equivalent lumped thermal network in the axial direction: (a) $P=0$, (b) $P\ne 0$.
In the text

	Fig. 5 Equivalent thermal network considering symmetry.
In the text

	Fig. 6 Equivalent lumped thermal networkin the radial direction: (a) $P=0$, (b) $P\ne 0$.
In the text

	Fig. 7 Equivalent thermal network of elementary block with symmetry.
In the text

	Fig. 8 Global equivalent thermal network (n blocks along r and one block along z).
In the text

	Fig. 9 Current density induced in billet vs r(Ω = 1500 rpm, θ = 0°).
In the text

	Fig. 10 Total power heat induced in billet vs. rotational speed.
In the text

	Fig. 11 Temperature evolution at the surface of the billet (r = 38 mm, Ω = 1500 rpm).
In the text

	Fig. 12 Temperature evolution in billet vs r (θ = 0°, Ω = 1500 rpm, t = 450 s).
In the text