Majid Souley
email: majid.souley@laplace.univ-tlse.fr

Julie Egalon
email: julie.egalon@laplace.univ-tlse.fr

Stéphane Caux
email: stephane.caux@laplace.univ-tlse.fr

Olivier Pateau
email: olivier.pateau@edf.fr

Yvan Lefèvre
email: yvan.lefevre@laplace.univ-tlse.fr

Pascal Maussion
email: pascal.maussion@laplace.univ-tlse.fr

Optimization of the settings of multiphase induction heating system

Keywords: induction heating, metal industry, current control, electromagnetic induction, multi phase

This paper deals with the setting parameter optimization procedure for a multi-phase induction heating system considering transverse flux heating. This system is able to achieve uniform static heating of different thin/size metal pieces without movable inductor parts, yokes or magnetic screens. The goal is reached by the predetermination of the induced power density distribution using an optimization procedure that leads to the required inductor supplying currents. The purpose of the paper is to describe the optimization program with the different solution obtained and to show that some compromise must be done between the accuracy of the temperature profile and the energy consumption, with the calculation of the losses.

INTRODUCTION

The advantages of classical induction heating systems (high power density can be reached quickly deep inside the workpiece) can be improved by multiphase inductors. Indeed, achieving metal temperature homogeneity with single inductor on a large scale is impossible, particularly in transverse flux heating [START_REF] Zinn | Elements of induction heating: design, control, and applications[END_REF]. A solution to overcome this constraint consists of using multi-coil systems thus leading to more complexity. Significant progresses [START_REF] Manot | Modélisation couplées des dispositifs électromagnétiques associés à des circuits d'électronique de puissance. Intégration de la commande des convertisseurs -Aide à la conception d'un dispositif de chauffage par induction à flux transverse[END_REF][START_REF] Souley | Simplified model of a metal disc induction heating system[END_REF][START_REF] Egalon | Multiphase System for Metal Disc Induction Heating: Modeling and RMS Current Control[END_REF][START_REF] Lucia | Multiple-Output Resonant Matrix Converter for Multiple Induction Heaters[END_REF] have been made in this domain by considering a static multi-inductor system without any moveable devices such as yokes or magnetic screens [START_REF] Zinn | Elements of induction heating: design, control, and applications[END_REF], [START_REF] Rapoport | Optimal control of induction heating processes[END_REF].

Besides, one solution [START_REF] Oleg | Induction heating device and process for controlling temperature distribution[END_REF] proposes manual changes of the coil connection for generator-load adaptation, some variable passive elements are added in [START_REF] Fujita | A New Zone-Control Induction Heating System Using Multiple Inverter Units Applicable Under Mutual Magnetic Coupling Conditions[END_REF] or some decoupling transformers between the different phases [START_REF] Uchida | Induction heating method and unit[END_REF]. Our system [START_REF] Egalon | Multiphase System for Metal Disc Induction Heating: Modeling and RMS Current Control[END_REF] is certainly simpler but not very easy to tune, due to mutual coupling between inductors and with the load. In fact, a trial and error setting campaign can be launched but is time and energy consuming. It is possible to overcome these disadvantages by the calculation of the setting parameters of a multiphase system as an optimization problem [START_REF] Manot | Modélisation couplées des dispositifs électromagnétiques associés à des circuits d'électronique de puissance. Intégration de la commande des convertisseurs -Aide à la conception d'un dispositif de chauffage par induction à flux transverse[END_REF], [START_REF] Egalon | Multiphase System for Metal Disc Induction Heating: Modeling and RMS Current Control[END_REF]. It consists of an optimization with constraints, based on a prior determination of the mapping of induced current distribution inside the heated material. It necessitates the global power density calculation in order to get the correct temperature profile through the setting of the correct currents in the 3 phases. In that case, amplitudes and phases have to be determined, except in [START_REF] Ha | Dynamic Analysis and Control for Resonant Currents in a Zone-Control Induction Heating System[END_REF] where the inductor currents are in phase. Two approaches leading to relatively short time computation are described in this paper. We will first present the three phase system which was the validation support of the method, followed by a description of the optimization procedure and the results.

II. DESCRIPTION OF THE THREE PHASE INDUCTION HEATING SYSTEM

The system consists of three inductors, each one having two coils in a serial connection ( Fig. 1) and organized face to face in a transverse flux configuration. A stainless steel disc plate (316L) which is the load to be heated, is placed inside the coils. Each coil is supplied by a current generated by an independent current inverter with a common current source. This architecture was already presented in [START_REF] Egalon | Multiphase System for Metal Disc Induction Heating: Modeling and RMS Current Control[END_REF] but is recalled hereafter in Fig. 2. The quasi sinusoidal shape of the inductor current I1 is due to the high quality factor of the LC resonant circuit. The final objective of our study is the heating of metal strips, but to begin with a simpler case, this paper deals with the achievement of a uniform heating of a metal disc with 10 MW/m 3 power density injected inside the material. The disc plate is 850mm diameter long with 1mm thickness. The resistivity of the material is 74.10 -8 ɏ.m. OPTIMIZATION PROCEDURE

The goal of the optimization procedure is to define a set of exciting currents which can assure a specified temperature uniformity in the metal disc cross-section along the radius. Due to the 1500 Hz excitation frequency and to the 1mm thickness of the metal disc, the induced current density should be uniformly distributed in the cross section along the radius. The penetration depth is 1.17mm; the resistivity variation with temperature has been neglected during the optimization. The approach is based on the fact that each inductor coil, individually supplied with its sinusoidal current, produces an induced current distribution whose representation is given in Fig. 4. These functions depend on the operating frequency, on the physical parameters of the material, on its geometrical dimensions and on the relative position of the material inside the coils. It is important to notice that this stainless steel has non magnetic properties. Indeed if magnetic material is used, the current distribution should change and hysteretic effets should be taken into account, but it is not the case here. According to the distribution functions in Fig. 4, supplying an inductor with a current of 1A, leads to an energy distribution that permits to predict the total induced current in the work piece by superposition of each current contribution.

݂ ҧ )ݔ( = ܬ ഥ )ݔ( ܫ ҧ ቤ ூҧ ೖ ୀଵ = ݂ ҧ ோ )ݔ( + ݆. ݂ ҧ ூ )ݔ( (III-1)
As the linearity of the system [START_REF] Zinn | Elements of induction heating: design, control, and applications[END_REF], [START_REF] Manot | Modélisation couplées des dispositifs électromagnétiques associés à des circuits d'électronique de puissance. Intégration de la commande des convertisseurs -Aide à la conception d'un dispositif de chauffage par induction à flux transverse[END_REF], [START_REF] Egalon | Multiphase System for Metal Disc Induction Heating: Modeling and RMS Current Control[END_REF] is verified, the current induced by inductor k can be written by the relation (III-1) and so on for the total induced current which is expressed as the superposition of each contribution (III-1).

ܬ ഥ )ݔ( = ݂ ҧ )ݔ( × ܫ ҧ ՜ )ݔ(ܬ = ܬ ഥ )ݔ( ଷ ୀଵ = ݂ ҧ )ݔ( × ܫ ҧ ଷ ୀଵ (III-2)
More precisely the total induced current is written by relation (III-3) requiring real parts ݂ Ը )ݔ( and imaginary parts ݂ Ա )ݔ( of the induced current by each coil, cf Fig. 4. The expression (III-3) reflects the superposition of each inductor curve and lead to the total current density induced in the material.

ܬ ҧ ,ݔ( ܺ) = (݂ Ը ,ݔ( ܺ). ܫ Ը ଷ ୀ ଵ െ ݂ Ա ,ݔ( ܺ). ܫ Ա ) + ݆. (݂ Ը ,ݔ( ܺ). ܫ Ա ଷ ୀ ଵ + ݂ Ա ,ݔ( ܺ). ܫ Ը ) = ܬ Ը ,ݔ( ܺ) + ݆. ܬ Ա ,ݔ( ܺ) (III-3)
where :

-ܫ Ը = ܫ . cos (߮ ଵ ) ܫ Ա = ܫ . sin (߮ ଵ ) are respectively the real part and the imaginary part of the coil supply current k,

-ܺ = ܫ[ ଵ , ܫ ଶ , ߮ ଶଵ , ܫ ଷ , ߮ ଷଵ ] is the current
vector to be optimized, -߮ ଵ is the phase shift of the current ܫ related to ܫ ଵ phase. Previous expressions allow the total power density (III-4) calculation at each abscissa along the radius.

,ݎ(ܦ ܺ) = ߩ ൭ (݂ Ը ,ݔ( ܺ). ܫ Ը ଷ ୀ ଵ െ ݂ ூ ,ݔ( ܺ). ܫ Ա )൱ ଶ + ൭ (݂ Ը ,ݔ( ܺ). ܫ Ա ଷ ୀ ଵ + ݂ Ա ,ݔ( ܺ). ܫ Ը )൱ ଶ (III-4)
Then the function (III-5) to be optimized is written under constraints on the source current and on the phase shifts between inductor current in phase 1 and the current in the other phases.

ܨ = ,ݔ(ܦห ܺ) െ ܦ ห ܦ ݎ݁݀݊ݑ ൝ 0 ܫ 1000 ܣ െ180° ߮ ଵ 180 ݇ = 1,2,3 ° (III-5)
IV. OPTIMIZATION RESULTS

1.

Function F is optimized using fmincon in Matlab with certain constraints. As a first step, the minimum value of the abscissa along the radius must be studied. Indeed as shown in Fig. 4 and considering that the heating is far more faster than the conduction flux inside the material, it is not possible to heat the center of the disc plate. There is no induced current in this place because of the inductors 'pan cake' configuration. So, looking for the most homogeneous heating of a large part of the metal disc with the less possible error becomes a goal. It can be reached considering an additional variable in the optimization program, which becomes an important parameter, the minimum radius from which a large area of the metal disc can be heated up, as homogeneously as possible, with a maximum error of 5%.

As a result of the determination, Fig. 5 shows three examples of power density distributions corresponding to 0.01 m, 0.08 m and 0.1 m for value. As in time control theory and because the delivered power starts from zero at the disc center, there must be a compromise between a kind of "space response" i.e. the part of the disc close to its center where the power comes into the required 5% band and a kind of "overshoot". The power density curve does not pass below the threshold set at 5% when is at least 0.1 m, which will be the chosen value. Then, the profile is optimized upon , with = .

as the radius of the metal disc [START_REF] Souley | Simplified model of a metal disc induction heating system[END_REF]. It is important to notice that the optimization procedure can also perform different power density profiles, with "mixed" results, which are shown in figures 6a to 7b. It can easily be seen that, the higher the required temperature in the center, the worst the profile obtained. The following parts of this paper will focus on a flat temperature profile. The current optimization results are summarized in TABLE I for a single criterion, focusing only on the accuracy of the temperature profile, after 100 iterations for a total calculation time of 37.5s on a PC computer and with random initial conditions. One capacitor bank is placed in parallel with each inductor supplied by an independent current inverter, with C = [420 225 143]µF in order to reach a resonant behavior of the inverters. As it can be seen in TABLES I to IV, different solutions are obtained. They lead to different power consumptions, source current values. Obviously in such a system, losses in the inverters and in the resonant capacitors must be taken into account. The way to calculate the losses will be described in the next part. The best solutions are in Moreover, as the optimization procedure is calibrated for choosing random initial point inside the constraint range, reducing the range leads to solutions that might not appear. So as a reduced range, the cases of െ90° ߮ ଵ 90°; െ90° ߮ ଵ 0° and 0° ߮ ଵ 90° are investigated in Tables II, III and IV respectively. Using െ90° ߮ ଵ 90° range leads to Table II, where the new solution which is not better than the best one in TABLE I.

I1(A) I2(A) M21 (°) I3(A) M 31 (°) p1(W) p2(W) p3(W) Ptot(W)
Considering െ90° ߮ ଵ 0 is an improvement and gives a better solution regarding efficiency but with a less accurate profile as the criteria res=1. The total power is 6356W instead of 6695W and the source current is reduced from 76A down to 55.4A, which can lead to a switch downsizing. Nevertheless, this power density profile is less accurate compared to first one in TABLE I.

The optimization has also been made with Hessian method and lead to the same results with a computation time of 31.1s.

VI. LOSSES EVALUATION

The evaluation of losses in the switches of the inverters goes through the study of the current and voltage waveforms. Each switch is composed by an IGBT in series with a diode. IGBTs accept the positive part of the voltage whereas the diodes stand the negative part. For operating points treated, components are oversized in this test bench.

Conventional methods [START_REF] Lucía | Efficiency optimization of half-bridge series resonant inverter with asymmetrical duty cycle control for domestic induction heating[END_REF] of calculating losses extracting values of threshold voltage and dynamic resistances of the static characteristics are not suitable here and the efficiencies are not as great as it can be seen in other systems [START_REF] Sarnago | Dual-modeoperation half-bridge resonant converter for improved-efficiency induction heating System[END_REF]. Indeed operating points are situated on nonlinearized parts of their V(I) characteristic. Consequently, a special modeling must be done on the entire curve. For this purpose, the curves are recreated in Excel® using the software PlotDigitizer®, as described in Fig. 8 and 9and a graphical curve from the datasheet in order to fit a number of points extracted in a table. Afterwards Excel® gives a trendline of a series of points. The five curves of the two components are estimated from the static characteristics, the IGBT turn-on and turn-off energies and the diode turn-off energy. The conduction losses for one switching cell are calculated by expression (VI-1). The switches are phase-shifted controlled and are ON during half a period.

ܲ ைே, = 1 2 ൫ܸ ܫ( ௌ ) + ܸ ா ܫ( ௌ )൯. ܫ ௌ (VI-1)
The turn-off and turn-on energies are provided for a given voltage ܸ ,ௗ . So, as to consider the real commuted voltage ܸ , switching losses are estimated multiplying the theoretical losses by the ratio , . The calculation for one switch is given by formula (VI-2), where ݂ is the switching frequency and also the operating frequency.

ܲ ௌ, = ܧ , ܫ( ௌ ). , , כ ݂ + ܧ , ܫ( ௌ ). , , + ݂ (VI-2)
The terms ܧ , ܫ( ௌ ) and ܧ , ܫ( ௌ ) are extracted either from the diode datasheet or from the IGBT datasheet, depending on the voltage sign. The diode turn-on energy is zero. The voltage values are caught as showed in Fig. 10. The ignition voltage is taken one simulation step before switching and the blocking one, one step after switching. Total losses are computed for one inverter in each configuration in Table I to IV, according to (VI-3).

ܲ ்ை், = 4 כ ܲ ைே, + 2 כ ܲ ௌ, భ + 2 כ ܲ ௌ, మ (VI-3)

VII. OPEN LOOP RESULTS

Fig. 11 compares some optimal solutions which were obtained. All these three solutions corresponding to the best solutions in Table I to III, allow a rather good profile, without any high difference from the reference. As a consequence, the choice of one solution rather than another is subject to the consideration of the accuracy regarding energy efficiency. Finally, the total power should be compared to the theoretical one, assuming that the power density follows the profile in Fig. 12. These considerations lead to a theoretical minimum power of 6200W which is not very different from the above mentioned values in Tables I-IV.

CONCLUSION

This paper has shown that an optimization procedure is necessary to determine the optimal supply currents in a multi-phase induction heating system. Several simulation results were shown, including different temperature profiles in the disc that has to be heated. A compromise between the accuracy of the temperature profile with respect to its reference and the energy consumption has been put in evidence. A single criterion based on the accuracy of the power density profile was used but the global losses in the resonant capacitors and in the inverters were calculated. The different solutions were compared either in terms of temperature profile or in terms of losses. The same optimization must be done as soon as the load is modified. Closed loop profiles are currently under development. Finally, an adaptive tuning is part of our long term perspectives.

Fig. 1 :

 1 Fig. 1: face to face 3-phase inductors of the system

Fig. 2 .

 2 Fig. 2. Whole system PSim model with inverters, controlled-current source, capacitors and coupled-inductors

Fig. 3 .

 3 Fig. 3. Overview of the 3 inductors of the system

Fig. 5

 5 Fig. 5 Different power density distributions as a function of for a 10 6 W/m 3 reference; (xmin=0.05_worst solution, xmin=0,08_middle ; xmin=0.1_best solution)

Figure 6 .Figure 6 .Figure 7 .

 667 Figure 6.a Decreasing possible exponential reference and simulated temperature profiles

77

  Abscissa along the radius (m) Calculated power density I1 = 308,35 A I2 = 100,71 A I3 = 476,39 A Ԅ21 -64,53 (°) Ԅ31 -37,06 (°) Calculated power density I1 = 486,86 A I2 = 247,14 A I3 = 223,61 A Ԅ21 -61,14 (°) Ԅ31 -63,64 (°) Dp (w/m 3 ) Abscissa along the radius (m) Reference : 4/0.425 MW/m 4

Fig. 8 :Fig. 9 :

 89 Fig.8: Modeling of a curve from the component datasheet with Plotdigitizer

Fig. 10 :

 10 Fig.10: Instants for back ignition and blocking voltages measurements on the switch voltages.

Fig. 11 .

 11 Fig. 11. Comparison of the 3 solutions for the temperature profile

Fig. 12 . 1 Fig. 13 .

 12113 Fig. 12. Theoretical power density profile

TABLE I SOLUTIONS

 I OVER RANGE OF PHASE SHIFT െૡ° ૡ°

		Inductor currents			Power consomption		Source & criteria	Global losses
	I1(A) I2(A)	M21 (°)	I3(A)	M 31 (°)	p1(W) p2(W) p3(W)	Ptot(W) Is(A)	res	Pswit	Pcapa
	253.9 114.7	-49.4	92.8	-63.1	1414 2735	2547	6696	76.1	0.2	1138	5,7
	262.6	0.0	179.9	158.2	3.8	3282	0.0	3485	6766	80.2	40.1	1207	6,5
	192.8 196.4	90.1	84.1	37.2	342	103	2932	6459	136.1	0.4	2397	5,4
	327.4 248.8	-79.7	128.4	180	-419	4091	4014	7686	137.6	3.2	2448	8,5
	0.0	281.8	-130.4	238.1	-4.3	0.0	8901	-2538	6363	159.1	109.3	2993	5.74
	466.9 350.6	-180	46.6	-180	2631 5121	845	8596	223.3	203.8	4755	13.3
	0.0	367.7	-158.6	168.8	66.4	0.0	1219	5073	6291	255.0	94.2	5557	6.88
	466.1 516.7	-180	178.2	44.1	3758	104	6231	10093 380.3	151.0	8803	19.2

TABLE III SOLUTIONS

 III OVER RANGE OF PHASE SHIFT െૢ° °

	Inductor currents	Power consomption	Source & criteria	Global losses