Mathematical Modeling of Eddy-Current Loss for a New Induction Heating Device

Abstract

A new induction heating device is presented in this paper. This device can convert mechanical energy into heat energy by utilizing eddy currents, which are induced by rotating permanent magnets. A mathematical model is established for estimating eddy-current loss of the device. The distribution of induced currents and the resultant magnetic field intensity are considered in the process of modeling the eddy-current loss and so is the mutual influence of the electric field between neighborhood pole projection areas. Particularly, the skin effect is considered by correcting the numerical integral domain of eddy current density, which has great effect on the calculating results. Based on specific examples, the effectiveness and correctness of proposed model are proved by finite element analysis. The results show that the mathematical model can provide important reference for design and structure optimization of the device.

1. Introduction

In the process of induced heating, the thermal power is due to eddy currents, when conductor is subjected to an alternating-current (AC) magnetic field. Most induction heating devices, with electric energy as input, produce the AC magnetic field by using coils [1–5]. In this paper, we present a novel induction heating device, with rotating mechanical energy as input, and the eddy currents are induced by the relative motion between conductor and permanent magnets.

This device is mainly composed of a stator and a dual-rotor. The stator is made of copper, inside which there are heat-exchange pipes with an inlet and an outlet, and the dual-rotor is made of iron. The structure of the device is shown in Figure 1. There is an alternative magnetic field between two rotors produced by permanent magnets with alternating pole configurations. When the dual-rotor and hence the magnetic field rotate, the fluid flowing through the heat-exchange pipes can be heated by the eddy-current loss of the stator. During this process, input rotating mechanical energy is converted into thermal energy. This device is similar to the eddy-current brake in structure and energy conversion. Many efforts have been done to the modeling of the eddy-current brake so far [6–8].

(a)

(b)

Figure 1 

Structure of the heating device. (a) Stator and dual-rotor and (b) heat-exchange pipes of stator.

Based on the eddy-current brake concept, we present a mathematical model of eddy-current loss for the induction heating device at low-speed operation. First, the cylindrical induction heating device with large curvature radius is simplified as linear model, so rotating movement is changed into linear movement. Second, the surface charge, due to the motional induced electric field, at the edge of the pole projection area is obtained. From these charges, the Coulomb field intensity is calculated by using Coulomb’s law. Meanwhile, the mutual influence between neighborhood Coulomb fields is considered. Third, the eddy current in the pole projection area, due to rotational magnetic field of the permanent magnets, is derived. Finally, the brake force is calculated by applying the Lorentz force law through numerical integral method. Based on the relationship between the brake force and the loss power, the eddy-current loss is obtained at different speeds. Since the eddy current generates the magnetic flux opposite to the applied magnetic flux, the induced magnetic flux density is calculated in a lumped way by using Ampere’s law, and the net magnetic flux density is defined by introducing the magnetic Reynolds number. In the process of the numerical integral, the skin effect is considered according to the characteristics of the eddy currents distribution. Also, the analytical results are compared with the finite element results and the causes leading to difference being analysed.

2. Analysis of Eddy-Current Loss

The cylindrical model of the induction device is shown in Figure 2. There are 𝑛 surface mounted permanent magnets radial magnetized on each rotor, with alternating pole configurations, and 𝑛 is even. Where 𝑙₁ and 𝑙₂ are the arc length of the permanent magnets on inner rotor and outer rotor, 𝑅 is the distance to the stator core center from the center of the cylindrical model; 𝑅₁ and 𝑅₂ are the circular radius of inner and outer permanent magnets.

Figure 2 

Cylindrical model.

To simplify the mathematical modeling, the cylindrical model is simplified as a linear one as shown in Figure 3.

Figure 3 

Linear model.

Where 𝑏 is the axial length of model, 𝑎(𝑎=(𝑙₁+𝑙₂)/2) is the length of permanent magnets, 𝑙 is the spacing distance between adjacent permanent magnets, 𝑑 is the thickness of stator, and 𝑔 is the air gap. With this model simplification, the anticlockwise rotating movement at an angular velocity 𝜔 is changed into linear movement of the stator, with the constant velocity 𝑣=𝜔𝑅 perpendicular to the direction of alternative magnetic field 𝐁.

As shown in Figure 4, the movable charges in the pole projection area of the stator are subject to the magnetic force 𝑞𝐯×𝐁 which directed upward along the stator. This action will lead to a separation of the charges, the positive ones moving to the top of the stator and the negative ones to the bottom. The Coulomb charges stack at the ends of the region of the pole projection and produce an electric field 𝐄 pointing downward that will tend to decrease the total force on the given charge moving upward in the interior [9].

Figure 4 

Polarization in pole projection area and charge distribution in the stator core due to motion of inducing magnetic field.

In the present case, the current density induced in the stator is given by 𝐉=𝜎𝐄^′, (1) where 𝜎 is the conductivity of the stator and the electric field intensity 𝐄^′ can be calculated as [10] 𝐄^′=𝐄+𝐄_𝑣, (2) where 𝐄 is the electrostatic field of Coulomb charge along the edges of the pole projection area. 𝐄_𝑣 is the electromotive field driving charges and is expressed as follows: 𝐄_𝑣=𝐯×𝐁. (3) The current density 𝐉 induced in the pole projection area can be written down by using (1), (2), and (3) as follows: 𝐉=𝜎(𝐄+𝐯×𝐁). (4) Since the electromotive field 𝐯×𝐁 is easily known, 𝐄 is the key step to obtain the current density 𝐉. In Figure 4, 𝐁=−𝐵𝐚_𝑧,𝐯=𝑣𝐚_𝑥. (5)𝐏_𝑦 is defined as the polarization with respect to the 𝑦 component and is written down as follows: 𝐏_𝑦=𝜀₀𝐄_𝑣, (6) where 𝜀₀ is the permittivity in free space. Thus, 𝐏_𝑦 can be expressed by using (3), (5), and (6) as follows: 𝐏_𝑦=𝜀₀𝐯×𝐁=𝜀₀𝑣𝐵𝐚_𝑦. (7) The surface charge densities on the plane parallel to the 𝑥𝑧 plane at 𝑦=±𝑏/2 are defined as 𝜌_𝑠𝑦±, which is analogously calculated by using the method of calculating the magnetic pole density at the surface of the permanent magnet [11]. In the present case, 𝜌_𝑠𝑦+ and 𝜌_𝑠𝑦− are written down as follows: 𝜌_𝑠𝑦+=𝐏_𝑦⋅𝐧_𝑦+=𝜀₀𝑣𝐵𝐚_𝐲⋅𝐧_𝑦+=+𝜀₀𝜌𝑣𝐵,_𝑠𝑦−=𝐏_𝑦⋅𝐧_𝑦−=𝜀₀𝑣𝐵𝐚_𝐲⋅𝐧_𝑦−=−𝜀₀𝑣𝐵, (8) where 𝐧_𝑦± is the unit outward normal vector on the planes parallel to the 𝑥𝑧 plane at 𝑦=±𝑏/2. If 𝜌_𝑙𝑦+ and 𝜌_𝑙𝑦− are defined as the line charge densities at the points (𝜁,𝑏/2) and (𝜁,−𝑏/2), they are written down as follows: 𝜌_𝑙𝑦+=+𝜀₀𝜌𝑣𝐵Δ𝜁,_𝑙𝑦−=−𝜀₀𝑣𝐵Δ𝜁, (9) where Δ𝜁 is the differential length in the 𝑥 direction, as shown in Figure 4. If 𝐥_𝑦+ and 𝐥_𝑦− are the displacement vectors from the points (𝜁,𝑏/2) and (𝜁,−𝑏/2) to the point 𝑃(𝑥,𝑦), the electrostatic fields 𝐄_𝑦+ and 𝐄_𝑦− by the line charge densities are calculated by Coulomb’s law: 𝐄_𝑦+=𝜌_𝑙𝑦+2𝜋𝜀₀𝐥_𝑦+𝐥_𝑦+²=𝑣𝐵2𝜋^𝑎/2_−𝑎/2(𝑥−𝜁)(𝑥−𝜁)²+𝑦−𝑏/2²𝐚_𝑥+𝑦−𝑏/2(𝑥−𝜁)²+𝑦−𝑏/2²𝐚_𝑦𝐄𝑑𝜁,_𝑦−=𝜌_𝑙𝑦−2𝜋𝜀₀𝐥_𝑦−𝐥_𝑦−²=−𝑣𝐵2𝜋^𝑎/2_−𝑎/2(𝑥−𝜁)(𝑥−𝜁)²+𝑦+𝑏/2²𝐚_𝑥+𝑦+𝑏/2(𝑥−𝜁)²+𝑦+𝑏/2²𝐚_𝑦𝑑𝜁. (10) Then, the total electrostatic field intensity 𝐄 is obtained by using the superposition theory as 𝐄=𝐄_𝑦++𝐄_𝑦−. (11)

In order to simplify the algorithm in this case, we consider that the line charges ±𝜀₀𝑣𝐵Δ𝜁 at 𝑦=±𝑏/2 extend indefinitely in the 𝑧 direction and from −𝑎/2 to +𝑎/2 in the 𝑥 direction, which is the source of 𝐄. Therefore, the 𝑧 component of the electrostatic field 𝐄 is ignored, and the total electrostatic field intensity 𝐄 is expressed as 𝐄=𝐸_𝑥𝐚_𝑥+𝐸_𝑦𝐚_𝑦, where 𝐸_𝑥 and 𝐸_𝑦 are the summations of the 𝑥 and 𝑦 components of 𝐄_𝑦+ and 𝐄_𝑦−, and they are written down as 𝐸_𝑥=𝑣𝐵)4𝜋ln(2𝑥+𝑎²+2𝑦−𝑏²(2𝑥−𝑎)²+2𝑦+𝑏²(2𝑥−𝑎)²+2𝑦−𝑏²(2𝑥+𝑎)²+2𝑦+𝑏²,𝐸_𝑦=𝑣𝐵2𝜋tan⁻¹2𝑥+𝑎2𝑦−𝑏−tan⁻¹2𝑥−𝑎2𝑦−𝑏−tan⁻¹2𝑥+𝑎2𝑦+𝑏+tan⁻¹2𝑥−𝑎.2𝑦+𝑏 (12)(13) Then, the eddy current density 𝐉 induced in the stator can be obtained from (4). In the pole projection area, the 𝑥 and 𝑦 components of 𝐉 are given as follows: 𝐽_𝑥=𝜎𝐸_𝑥𝐽_𝑦=𝜎𝐸_𝑦+𝑣𝐵. (14) By using the Lorentz force law, we have the magnetic braking force as follows: 𝐹_𝑥=−_𝑣𝐽_𝑦𝐵𝑑𝑉=−𝑑^𝑎/2_−𝑎/2^𝑏/2_−𝑏/2𝐽_𝑦𝐵𝑑𝑥𝑑𝑦.𝐹_𝑦=−_𝑣𝐽_𝑥𝐵𝑑𝑉=−𝑑^𝑎/2_−𝑎/2^𝑏/2_−𝑏/2𝐽_𝑥𝐵𝑑𝑥𝑑𝑦. (15)(16)

In the outside of the pole projection area, the magnetic flux density 𝐵 is zero, and the eddy currents cannot produce the braking force; therefore, 𝑣 is the volume of the pole projection area in (15) and (16). Since the number of poles 𝑛 is even, the summation of 𝐹_𝑦 is zero. The 𝐽_𝑦 inside the pole projection is considered for the determination of 𝐹_𝑥 from (15), and we just need to calculate the 𝐸_𝑦 by (13).

For the model of multipoles with alternating (NSNSNS) pole configuration, the magnetic flux of neighboring poles is symmetrical distribution of interest pole with respect to 𝑧-axis, and the surface charges at the edge of the neighboring pole projection area will affect the electric field of the interest pole. Thus, the expression of electric field intensity 𝐸_𝑦 should be revised. The electric field intensity of the interest pole can be expressed as 𝐸_𝑦𝑖𝑥,𝑦=𝑣𝐵𝐾𝐾2𝜋𝑥,𝑦𝑥,𝑦=tan⁻¹2𝑥+𝑎2𝑦−𝑏−tan⁻¹2𝑥−𝑎2𝑦−𝑏−tan⁻¹2𝑥+𝑎2𝑦+𝑏+tan⁻¹2𝑥−𝑎.2𝑦+𝑏 (17) Hence, the electric field intensity of the neighboring poles can be written down as follows: 𝐸_{𝑦(𝑖−1)}𝑥,𝑦=−𝐸_𝑦𝑖𝑥+𝑎+𝑙,𝑦=−𝑣𝐵𝐾𝐸2𝜋𝑥+𝑎+𝑙,𝑦_𝑦(𝑖+1)𝑥,𝑦=−𝐸_𝑦𝑖𝑥−𝑎−𝑙,𝑦=−𝑣𝐵𝐾.2𝜋𝑥−𝑎−𝑙,𝑦 (18) Thus, the superposition electric field intensity inside every pole projection area is expressed as 𝐸_𝑦=𝐸_𝑦𝑖+𝐸_𝑦(𝑖+1)+𝐸_{𝑦(𝑖−1)}. (19)

In (15), the original magnetic field 𝐵, produced by permanent magnets, will be decreased by the induced magnetic field generated by eddy currents. The induced magnetic flux density is calculated from the eddy current 𝐈 by using Ampere’s law in order to achieve accurate results. The eddy current 𝐈 is obtained from the eddy current density 𝐉 by two-dimensional numerical integration. Here, we ignore the 𝐽_𝑥, since it is much smaller than the 𝐽_𝑦 in the pole projection area, and the 𝑦 component of 𝐈 is expressed as 𝐼_𝑦=2^𝑎/2_−𝑎/20^𝑑/2𝐽_𝑦𝑑𝑥𝑑𝑧. (20) Here, we ignore the skin effect in the stator. Thus, 𝐽_𝑦 is regarded as being uniform distribution on the 𝑥𝑧 plane, and the thickness of stator 𝑑 is regarded as the size of the eddy-current distribution in the 𝑧-axis. Skin effect means that the current density in stator core varies exponentially inward from the surface in a time varying magnetic field. The distance in which the current density decreases to 1/𝑒 of surface value is called the nominal depth of penetration which is given by the relation 𝛿=2𝜌,𝜔𝜇 (21) where 𝜌 is the resistivity of stator core and 𝜇 is the absolute permeability of stator core. With the speed increasing, the skin effect quickly pushes the induced currents toward the surface of the stator and depletes the core of the stator from inducing eddy currents. Therefore, the thickness of the stator 𝑑 is not appropriate to be regarded as the size 𝐽_𝑦 distribution in the 𝑧-axis. According to the theory of two-dimensional wave propagation in the electromagnetic field, the 𝐽 in the stator is expressed as [12–14] 𝐽=𝐽_𝑚𝑒^−𝛼𝑦𝜋cos𝜔𝑡−𝜆𝑥−𝛼𝑦𝜆=𝜏𝛼=𝜎𝜇𝜔2=𝜎𝜇𝜋𝑛_𝑝𝑛, (22) where 𝐽_𝑚 is the current density amplitude on the outer surface of stator, 𝜏 is the polar distance, 𝜎 is the conductivity of the stator, 𝑛 is the number of magnetic poles, and 𝑛_𝑝=60𝜔/2𝜋𝑛 (r/min) is the speed. The skin effect is considered by using 𝐽_𝑦𝑒^−𝛼𝑧 instead of 𝐽_𝑦, and as a result the numerical integral domain of eddy current density needs to be corrected in (20). We introduce an equivalent size 𝑑^′/2 instead of the constant 𝑑/2 to obtain the correct eddy current. The 𝑑^′ is a function of the 𝑛_𝑝 and is defined as 𝑑^′𝑛_𝑝=1−𝑒^−𝑑𝛼𝛼, (23) where from (23), the equivalent size 𝑑^′(𝑛_𝑝) decreases with the increase of speed 𝑛_𝑝. Therefore, the revised 𝐼_𝑦 is a function of the 𝑛_𝑝 and is expressed as 𝐼_𝑦=2^𝑎/2_−𝑎/2^{𝑑′(𝑛𝑝}₀^)/2𝐽_𝑦𝑑𝑥𝑑𝑧. (24)

By using Ampere’s law, the induced magnetic flux density 𝐵_in has the form of 𝐵_in=𝜇₀𝐼_𝑦/2𝑔, which decreases the original magnetic flux density 𝐵. The resultant magnetic flux density 𝐵_𝑚 is expressed as 𝐵_𝑚=𝐵𝑒^{−𝐵in/𝐵}=𝐵𝑒^−𝑅𝑚, and 𝑅_𝑚 is the magnetic Reynolds number, which is defined as 𝐵_in/𝐵. Thus, considering the induced magnetic flux, (15) is rewritten down as 𝐹_𝑥=−_𝑣𝐽_𝑦𝐵_𝑚𝑑𝑉=−𝑑^𝑎/2_−𝑎/2^𝑏/2_−𝑏/2𝐽_𝑦𝐵_𝑚𝑑𝑥𝑑𝑦. (25) Also, the eddy-current loss has the form 𝑃_𝑒=𝐅⋅𝐯=𝐹_𝑥𝐚_𝑥⋅𝑣𝐚_𝑥=𝐹_𝑥𝑣. (26) Thus, we have the total eddy-current loss of the induction device given by 𝑃_𝑒=𝑛𝐹_𝑥𝑣=−𝑛𝑑𝑣²𝐵_𝑚×𝜎𝐵_{𝑎/2−𝑎/2𝑏/2−𝑏/2}1𝐾2𝜋𝑥,𝑦−𝐾𝑥−𝑎−𝑙,𝑦−𝐾𝑥+𝑎+𝑙,𝑦+1𝑑𝑥𝑑𝑦. (27) The calculated results can be obtained through the numerical integral method. This power loss is the heat source of heat transfer fluids in the stator.

3. Validation by Finite Element Analysis

Two specific examples with 20 and 10 poles are given, which have the same structure parameters and physical parameters except for the number of poles. Based on these parameters listed in Table 1, transient FEM models and mathematical models considering two cases of with and without skin effect are built, and the comparison between mathematical analysis results and FEM results has been carried out.

Table 1 

Parameters.

The results of the two examples with 20 and 10 poles are shown in Figures 5 and 6, respectively. The analytical model without considering the skin effect has small application range, and the deviation between the analytical results and the FEM results appears when the rotation speed is just over 200 r/min. However, the analytical model with considering the skin effect has wider application range, which agreed well with the FEM results within 400 r/min range of rotation speed.

Figure 5 

Analytical results and FEM results (𝑛=20).

Figure 6 

Analytical results and FEM results (𝑛=10).

The influence of skin effect is related to the number of poles. The less the number of poles is, the smaller the frequency of AC magnetic field and the influence of skin effect are, at the same rotation speed. As shown in Figure 7, the example with 10 poles has the smaller difference between two mathematical models with and without skin effect than the one with 20 poles.

Figure 7 

Difference between two mathematical models with and without skin effect.

The analytical results have been lagged behind the finite element results, and the difference increases with the speed increase, as shown in Figure 8. When the rotation speed exceeds appropriate range, the analytical model cannot provide preliminary estimations of eddy-current loss. The main reason is considered that the eddy-current loss is calculated only in the volume of the pole projection from (27), and the loss outside the pole projection is ignored which increases with the speed increase. The magnetic flux distribution at 0 r/min and 400 r/min speeds is shown in Figures 9(a) and 9(b). The available magnetic flux inside the pole projection decreases relatively with the increase of speed. On the contrary, the magnetic flux outside the pole projection increases which makes no contributions to the analytical model of eddy-current loss.

Figure 8 

Difference between analytical results and FEM results.

(a)

(b)

Figure 9 

Distribution of flux at speed of (a) 0 r/min and (b) 400 r/min.

4. Conclusion

In this paper, an induction heating device developed for heating fluid has been presented. Driven by power installation such as turbine, motor, and fan, this device can realize the conversion of the wind energy, water energy, and recovered energy of other power installations into thermal energy.

A mathematical model of eddy-current loss for this device has been proposed based on the eddy-current brake concept. In modeling the eddy-current loss, the skin effect has been considered by correcting the numerical integral domain of eddy current density, which effectively extends the applied speed range of the mathematical model. Based on specific examples, this model was validated by FEM analysis. This model shows an explicit relationship between structure parameters and power performance, which is useful for design and structure optimization of the induction heating device.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work was supported by the National Natural Science Foundation of China (NSFC 51207031).