Modeling of induction heating equipment for the vulcanization process

The paper explores the problem of induction heating of press equipment for the vulcanization of rubber products. We describe a mathematical model of induction heating of press plates in the threedimensional formulation given the automated stabilization of temperature, using positional and PID-control. The solution results are compared with the industrial experimental data.


Introduction
Hydraulic presses with heating plates, maintaining a temperature of 250 °C on the working surface for the vulcanization of rubber products and a temperature to 550 °C for treating metal products with high level of output, are often used to implement the processes of heating and pressure treatment of products. The resistive heaters or inductors can be used as heating elements.
The advantages of induction heating include high power density, reliability and durability of heaters. However, with an increase in the working temperature of the plate surface, the inductors capacity decreases a lot. In such cases, it is possible to provide the desired rate of heating by increasing the frequency of the supply voltage.
The capacity of resistive heating elements is essentially independent of the plate temperature. The disadvantage this type of heaters is a limited service life (as the wire is used it will gradually oxidize, forming areas of local overheating which will eventually lead to its destruction). In this regard, in the vulcanization process, preference is given to induction heating plates with current industrial frequency (50 Hz).

Statement of the problem
To reduce material costs at the stages of design, construction and operation of the equipment, improve its process characteristics it is necessary to optimize the design of the heating plate and to obtain the desired configuration of the temperature field on the work surface.
Mathematical modeling of induction heating processes is associated with the formalization of the related non-stationary electromagnetic and thermal processes in the volume of plates.
Low-frequency fields can only be effectively used for volumetric heating of ferromagnetic materials. The nonlinear dependence of magnetic field induction on its strength in ferromagnetic bodies considerably complicates the calculation of heat release. The existing methods for calculating the induction heating equipment have the following disadvantages: 1) the use of engineering calculation methods [1] based on the solution of Maxwell's linear equations in onedimensional formulation, using empirical data does not allow for production and analysis of the specific heat distribution in ferromagnetic bodies of a complex shape; 2) the use of simplified mathematical models of induction heating process do not take into account the non-linearity characteristics of ferromagnetic materials [2,3] and can lead to significant errors in calculating the heat generation by eddy currents; 3) the approximation of the inductor geometry by the bodies of canonical shape [4] leads to large errors in most practical problems.
In this regard, it is necessary to develop effective approaches to mathematical description of electromagnetic processes and solving the systems of equations of mathematical models, which, on the one hand, take into account all the features of the processes in question, and on the other hand, do not lead to an excessive amount of computation.

Mathematical model of induction heating of a single plate of the press
A mathematical model for the induction heating process of a hydraulic single-plate press includes the formula to calculate the heat generation by eddy currents and the equation for heat distribution in the volume of plate given the external heat loss.
Mathematical modeling of eddy current fields is based on the allocation of two regions in the volume of the plate and the surrounding airspace: Ω1 is the region of eddy currents with non-zero electrical conductivity; Ω2 is the region without eddy currents, which might contain external source currents.
When considering the problem of induction heating of press plates the region Ω1 is represented by the plate material, and the region Ω2 comprises inductors and surrounding airspace.
We make the following assumptions: 1) Maxwell's displacement current is absent (we consider low-frequency fields in a conducting medium); 2) the phenomenon of magnetic hysteresis is absent (high power electromagnetic fields are used) [1,5]; 3) properties of materials are isotropic.
Then for the region Ω1 differential equations of the electromagnetic field (Maxwell's equations) are of the form: where H is magnetic field strength, A/m; E is electric field strength, W/m; γ is specific electric conductivity, Ohm -1 m -1 ; B is magnetic induction, T; t is time, s. For the region Ω2: where J ext is current density supplied from an external source, A/m 2 . Boundary conditions for the contact area of regions Ω1 and Ω2 are used as single-valued conditions for the variables H, E, and B, as well as on the remote edge [6].
To solve Maxwell's equations we introduced the following auxiliary functions: magnetic vector potential A (Wb/m) and electric scalar potential V (W), which are connected with the vectors B and E by the following relationships: Then the system of equations (1) -(3) can be written as equation where μ is absolute magnetic permeability, H/m. Similarly, we can write the system of equations (4) -(5): From (7) it follows that the electric scalar potential V is determined up to a constant. To resolve the uncertainty, is enough to set the value of the potential V at a point of the computational domain. According to (6) magnetic vector potential A is determined up to an arbitrary gradient of a scalar function. The uniqueness of solution is achieved by calibration. At present, for the calculation of the eddy current fields the Coulomb gauge is most commonly used [7] From the computational perspective, it is extremely difficult to ensure that field A is divergence-free under non-permanent magnetic permeability [8]. For the first time, the problem of satisfying the Coulomb gauge condition in a three-dimensional formulation was successfully solved by the authors in [6], who modified the differential equations (8) and (9), adding to their left parts the penalty function , and formulated the boundary conditions to ensure a unique solution.
Then for the region Ω1 modified equations will have the form: For the region Ω2: Boundary conditions on the remote edge (the condition of zero for the normal component of the magnetic induction): where n is a unit normal vector to the surface of the boundary. Boundary conditions at the interface of the regions Ω1 and Ω2: where indices 1 and 2 show that characteristics belong to the corresponding region. In [6] it is proved that the system of equations and boundary conditions (11)-(19) is equivalent to the system of Maxwell's differential equations and provides automatic execution of the Coulomb gauge (10).
To calculate the heat release in the volume of heating plate the total current density is required: where J ind is density of eddy (induced) currents (A/m 2 ) found by the expression The vector of external source current density is determined in the general case by the value of the current by inductors and their structural characteristics: where Ii is current passing through the i-th inductor, A; G 〈i〉 is vector of design characteristics of the i-th inductor; nind is the number of inductors in the heating plate.
The module of average current density of the inductor can be calculated by the formula: where ωi is the number of winding turns of the i-th inductor; bi, hi are the width and depth of the slot crosssection for the i-th inductor, respectively, m. The current direction is expressed by the basis vectors of the coordinate system, and depends on the geometry of the inductor. Structural characteristics of the most common inductors of rectangular shape can be expressed in the form of the following vector: where di is wire diameter, m; x с i , y c i are coordinates of the center of the i-th inductor, m; li, si are its length and width, m.
We note that the external source currents and the eddy currents do not spatially overlap, and the equation (20) -the union of regions Ω1 and Ω 2, is necessary for the correct record of relations characterizing the distribution of heat in the volume of plates.
Thus, to calculate the field of eddy currents in each calculation point of the region Ω1 it is necessary to find the 4 unknowns (3 components of the magnetic vector potential A and the value of the scalar electric potential V), in the region Ω2 it is necessary to find the 3 unknowns (vector A). As noted above, the problem of determining the unknowns is significantly complicated by the non-linearity of magnetic properties of the material to be heated. When an external current source is sinusoidal and we assume the constancy of magnetic permeability of ferromagnetic material of the plate, the expressions (11)-(13) are simplified to linear equations in quasi-stationary complex representation: where j is an imaginary unit (j 2 = -1). Similarly, boundary conditions (14) -(19) are simplified.
Computing time spent to solve the linear equations (11′)-(13′) are reduced by approximately one order, the main difficulty is related to the replacement of the magnetization curve with the linear dependence. To describe the process of heat distribution in a plate the non-stationary heat equation, taking into account the internal heat sources is used: where T = T(x,y,z,t) is temperature in the point of the volume of the plate with coordinates (x,y,z) at the time point t, о С; q = q(x,y,z) is specific heat generation in this point, W/m 3 ; is coefficient of thermal diffusivity of the plate material, m 2 /s; c, ρ, λ are its specific heat capacity (J/(kg·K)), desnity (kg/m 3 ) and heat conductivity coefficient (W/(m·K)), respectively.
The average specific heat generation in the plate for the steady state mode is determined by the Joule-Lenz's law: where f is the current frequency, Hz; J is the total current density module, see (20); tst is conditional stabilization time of electromagnetic processes, after which the changes in the amplitude values of the current density can be neglected, s. When using linear equations of the electromagnetic field (11′)-(13′) heat generation is calculated by the formula: is the complex amplitude current density module. The initial condition for the solution of equation (25) is where Т0 is the ambient air temperature, о С. Heat transfer from the working surface, the cover and the ends of the heating plate is described by boundary conditions of the third kind: where Ωpl,r is the r-th surface of the heating plate, m 2 ; αr is the heat transfer coefficient from the r-th surface of the plate, W/(m 2 ·K), found by the criterion equation according to [9]; Tr is the average temperature of the r-th surface of the plate ( о С). We used a system of finite-element analysis of ANSYS to solve the equation (25) for an industrial induction heating plate sized 500×410 mm with four identical heaters with a total capacity of 5.5 kW heated from the initial temperature Т0 = 12 о С to the operating temperature of 170 ° C. Fig.  1 also shows the results of the industrial experiment conducted at JSC "ARTI-Zavod", Tambov to study the process of heating the plate using several thermocouples: in the corners and in the center of the work surface (see [10]). The experimental data for the control thermocouple are shown as rectangular markers, and those for the central thermocouple are shown as triangular markers.
The maximum deviation of the calculation results from the experimental data was 6.1° C and 5.3° C (at time 180 s) for the control thermocouple and central thermocouple, respectively. The average error for the heating period of 0-1968 plate was 4.6% and 4.9% for the control thermocouple and the central thermocouple, respectively. We note that with the heating of the plate, at time of 1470 s, the deviation did not exceed 2% and 3% for the control thermocouple and the central thermocouple, respectively. The mode of stabilization of a given operating temperature on the surface of the induction heating plate can be achieved by using a two-position controller (heaters on or off). To calculate the temperature field of the plate, we propose the following expression: where ( ) ( ) ( ) is the heaters condition ("1" -ON, "0" -OFF); Ts(t) is the plate temperature in the point of control thermocouple localization, °С; Td, Tu are lower and upper response thresholds of the controller, respectively, °С.

Results and discussion
We described a mathematical model of the induction heating process and automatic stabilization of temperature on the working surface of a hydraulic single-plate press. The model takes into account the nonlinearity magnetic characteristics of the plate material and can be applied to a three-dimensional analysis of the eddy current fields.
Calculation of temperature fields of induction heating plates with a certain distribution of heat generation by the eddy currents does not present any difficulties. The main difficulties arise at the stage of solving the threedimensional nonlinear equations of the electromagnetic field.
It is shown that in the case of constant magnetic permeability the electromagnetic field equations are greatly simplified. Nevertheless, the direct use of the simplified model when calculating the eddy currents in the ferromagnetic plate is not possible. In this regard, it remains topical to develop the relevant methods for calculating magnetic permeability equivalent to the magnetization curve by the criterion of active power.
The comparison of the calculation results with the industrial experimental data brought us to the conclusion about the adequacy of the proposed mathematical models to the real electromagnetic and thermal processes in the volume of the induction heating plate of a hydraulic press.