A Mathematical Model of Electrical Arc Furnaces for Analysis of Electrical Mode Parameters and Synthesis of Controlling Inﬂuences

: The synthesis of a real-time model of the electric mode of an arc steelmaking furnace is shown. This model is based on the method of medium-voltage and can be used to analyze the impact of the arc steelmaking furnace on the supply network at different stages of the melting process at different control systems of the electromechanical electrode movement system, and to create model-predictive control systems of electrode movement, when control using hardware in the loop technology. To implement the chosen method, the dynamic arc resistance was used. The last was obtained on the basis of the averaged experimentally taken volt-ampere characteristics of the arc. In paper was proposed its approximation as a current function and arc length. The simulation results obtained using the created model for different operation modes of the arc steelmaking chipboard for a system with a traditional regulator of the electrode movement system are demonstrated.


Introduction
Modern high-tech systems, which are increasingly used in industry, place increasing demands on the quality of electricity.As a rule, the basic requirements for the quality of electricity in modern production are put forward to ensure the appropriate parameters of the supply voltage and harmonic composition of the current.Violation of electricity quality parameters can cause significant losses to the manufacturer both due to violation of the technological mode of production, and due to equipment failure or forced downtime.The main threats in terms of electricity quality are long-term and short-term voltage fluctuations, the appearance of harmonic current components due to the connection of nonlinear, abruptly changing and asymmetric loads.An example of such a load is an arc steelmaking furnace [1][2][3].It should also be noted that steel arc furnaces are a source of rapid voltage fluctuations [4,5], which occur due to the asymmetry of currents and voltages.Today, the melting of scrap metal in arc steel furnaces provides about 1/3 of annual world steel production.A typical arc steelmaking furnace (EAF) during its operation consumes approximately 450 kWh of electrical energy for smelting one ton of steel [6].The principle of melting scrap in chipboard is to convert electrical energy into thermal energy, which is dissipated from electric arcs and melts the materials loaded into the furnace.Arcs burn between graphite electrodes and conductive scrap (3-phase system without neutral) and are characterized by low voltage (400-1000 V) and high currents (40-60 kA) for each phase.The temperature in the arc column can reach 8500 K, which provides the ability to quickly and efficiently melt scrap.In alternating current furnaces, the power supply system consists of a power transformer with variable taps and, sometimes, a series reactor connected to the high voltage side to increase the reactance of the system.The electric arc length is regulated by the electrodes movement system, which provides to set a given power into the furnace, which is determined by the melting mode.In the melting process, the harmonic distortions of voltage and current signals change from larger values of the harmonic distortion factor (THD) in the initial stages of melting to smaller in the stages of purification and heating [7,8].Random change of the arc length during the melting process leads to a significant change in the parameters of the electric mode and the corresponding reaction on the part of the electrode movement control system.It is necessary to deal with systems by different tempo both in the system of electrode movement and taking into account the interaction between the electrode movement system and electromagnetic processes in the furnace space.This, in turn, will require the use of special approaches when conducting research using the created mathematical model of EAF.
Thus, the model for studying the influence of EAF on the power supply network should contain a model of the electric circuit taking into account the nonlinearity of the volt-ampere characteristic of the arc, the model of the control system of EAF electrode movement and the model of random arc length processes [9][10][11][12][13].Depending on the class of tasks, researchers often use simplified models of individual subsystems.In particular, in the synthesis of EAF electrodes movement control systems, the model of the power supply system is ignored as a whole, using only Volt-Ampere characteristics at different arc lengths [14 -16].On the other hand, when studying the influence of EAF on the electrical network, the processes in the EAF electrode movement are not taken into account, using Volt-Ampere characteristics of arcs at given values of arc length or given laws of its change [1][2][3]8,[17][18][19].In both cases, the Volt-Ampere characteristics of the arc are used, and therefore it is not surprising, that much of the work is devoted to the modeling of the electric arc [7,9,[20][21][22].Thus, a significant part of EAF models used in the control systems studing, does not take into account nonlinearities and changes in parameters.They also are based on the use of static Volt-Ampere characteristics.
The above works did not address the task of ensuring the functioning of models in real time.At the same time, the creation of such models in combination with the technology of "hardware-in-the-loop" significantly expands the possibilities of research and synthesis of modern control algorithms.The purpose of this article is to develop such a EAF model.
A feature of EAF is the presence of different tempo subsystems and a significant range of changes in its parameters in the melting modes (in particular, due to changes in the dynamic resistance of the arc).Reducing the numerical integration step to ensure the required modeling accuracy significantly slows down the calculation and hinders the operation of the model in real time.This complicates the task of choosing the method of numerical integration for real-time models, as noted in [23][24][25].In particular, in [25] it was noted about the difficulties in applying the widely used trapezoidal rule and backward Euler method for this class of problems.In particular, the application of the backward Euler method requires a significant reduction in the step of numerical integration and thus an increase in computing power and speed of calculation.Instead, the use of the trapezoid method, although it allows to increase the step of numerical integration, however, can lead to numerical oscillations and violations of the stability of the calculation.Therefore, it is important to apply new approaches to real-time modeling of EAF systems.
One of such approaches, proposed in [23], will be used in this paper to create a comprehensive model for the analysis of the electrical operation mode of EAF.A feature of this approach is the application of a new method of algebraization of differential equations-the method of average voltages at the step of numerical integration (average voltage values in integration step AVIS).The advantage of this method is the increased numerical stability of the calculation, as well as the possibility of calculation with an increased step of numerical integration and maintaining the required accuracy of the calculation, which is justified in the works [23,24].
This approach has proven itself in the synthesis of real-time models of various electromechanical systems [24][25][26] and has been used to set up excitation systems for synchronous generators [27,28], including at the South Ukrainian NPP.The application of this approach in the EAF model will require the creation of a slightly different from the existing model of the arc space.Synthesis a real-time EAF model, which can be used in both-model predictive control systems and to establish control systems for the movement of EAF electrodes and analysis of the impact of EAF on the power supply network in different control systems is an urgent problem to solve.
The article is organized as follows: Section 2 presents, for the first time, a mathematical model of an electric branch with an arc, created using the method of average voltages at the step of numerical integration.Section 3 describes a mathematical model of arc space, in which new methods are used to approximate the volt-ampere characteristics of the arc and determine the dynamic resistance of the arc as a parameter of the developed model.Section 4 describes the mathematical models of the of the power supply system of the EAF elements, namely: short network with furnace space (arc), furnace transformer and power supply network with choke, created using the AVIS method.The peculiarity of these models is their representation in the form of multipoles.Section 5 presents the calculation scheme of the EAF power supply system formed by a combination of the previously described multipole elements and briefly describes the calculation algorithm (a detailed description of the algorithm for mathematical modeling of electrical systems using AVIS and multipole theory is given in [24]).Section 6 for the first time describes a mathematical model of the chipboard electrode movement system created using the AVIS method.Section 7 presents the results of studies of established symmetric and asymmetric modes in the EAF's power supply system, as well as transient modes in EAF when the arc length is changing and during melting.The adequacy of the model was confirmed by comparison with the simulation results obtained by other authors using other approaches and the results of the experiment.

Mathematical Model of an Electric Branch with an Electrical Arc
In accordance with the main provisions of the method of average voltages in numerical integration step [23], the electric branch with inductance, active resistance and arc (Figure 1) is described by the following equation: where u R , u L , u arc -instantaneous values of voltage on active resistance, voltage on inductance, arc voltage, respectively; u-voltage applied to the branch; t 0 -is the start time of integration step, ∆t-is the numerical integration step.The instantaneous values of voltages on the active resistance and on the arc are represented by the formulas: where u R0 , u arc0 -voltage values at the start of the step; ∆u R , ∆u arc -voltage gains at the step of numerical integration.The voltage increments at the step of numerical integration are determined using the Taylor series: where -k-th voltage derivatives on the active resistance and on the arc at the beginning of the step; m-is the number of derivatives taken into account (method order).
Taking into account the dependences between voltage and current for the active resistance of the circuit R and for the arc, respectively: u R (t) = Ri(t), u arc (t) = f (i), we obtain expressions for increments (4) and ( 5): For the 1st order method: where R d = du arc di -dynamic arc resistance (determined on the basis of instantaneous values of arc voltage u arc and current (i); R d0 -dynamic arc resistance at the beginning of the step.Note that the dynamic resistance of the arc depends on the current of the arc and the length of the arc and is determined at each step.
Taking into account (2), ( 6) and the equation for the voltage on the inductance (the value of the inductance L is constant) u L (t) = L di (t)  dt from Equation (1) we obtain for the method of the 1st order: where i 0 -current of circuit at the beginning of the integration step, i 1 -is the current at the end of the integration step, U = 1 ∆t t 0 +∆t t 0 udt-the average value of the voltage applied to the circuit at the numerical integration step.
From Equation ( 7) we obtain the expression for determining the current of the branch at the end of the step: The effective value of the arc voltage at each step of integration is calculated by the following formula: where ∆i = i 1 − i 0 -the increase in arc current on the step.For the method of average voltages integration step of the 2nd order, the equations of voltage increments at the active resistance and at the arc will look like: In this case, from Equations ( 1), ( 2), ( 10) and ( 11) we obtain: The expression for determining the circuit current at the end of the integration step will take the form:

Mathematical Model of Arc Space of Electrical Arc Furnace
According to the chosen approach to EAF modeling, the description of the arc space does not use the arc model in the form of a volt-ampere characteristic, but rather the change of the dynamic arc resistance, which depends on both the current and the arc length.We demonstrate the application of this approach on the example of EAF-200.Experimentally obtained averaged dynamic Volt-Ampere characteristics for such a furnace are shown in Figure 2. The dynamic Volt-Ampere characteristic of the arc has different curves when the current decreases and increases, and is also asymmetric.Then, obtained on the basis of experimental dynamic Volt-Ampere characteristics of the arc, the dependences of the dynamic resistance of the arc R d (I arc , l arc ) for different lines will have the form shown in Figure 3.
The kind of the Volt-Ampere characteristic, which corresponds to the current drop, for a single arc length can be successfully approximated as: where A 0 and β approximation coefficients.Then the equation for determining the dynamic resistance R d will look like: The change in dynamic resistance R d obtained on the basis of experimental volt-ampere characteristics and its approximation by the equation for R d (I arc ) when I arc = 10 mm are given in Figure 4. Significant deviations observed in the regions of large currents are explained primarily by the shortcomings of the construction of the Volt-Ampere characteristic based on experimental data.In the case of current increasing, the approximation formula is much more complex.Analysis of different models of the current-voltage characteristic curves approximation is given in [22].We need the derivative of a rather complex model to approximate the arc volt-ampere characteristic to determine R d .In our opinion, it is more appropriate to approximate the dependence R d (I arc ) for a given arc length.
If we ignore the asymmetry, the obtained dependence of the change in dynamic resistance resembles the curve obtained when using the Mexican hat wavelet [29].
The equation for approximating the characteristic R d (I arc ) for a given arc length can be described as: where a 11 , a 21 , a 31 , a 41 are required parameters.
The change in dynamic resistance R d obtained on the basis of experimental voltampere characteristics and its approximation using (16) for I arc = 10 mm at a 11 = 9.4, a 21 = 0.38, a 31 = 18, a 41 = 2 and σ = 51, 000 are given in Figure 4b.To take into account the asymmetry of the Volt-Ampere characteristic of the arc, it is necessary to separately determine the coefficients a 11 , a 21 , a 31 , a 41 and σ at I arc < 0 and I arc ≥ 0.
Approximating the obtained values of the coefficients at different arc lengths using polynomials of the appropriate order, we obtain the following equation to determine the dynamic resistance as a function of current and arc length: if I arc • dI arc dt > 0 then: The obtained dependences of the change of dynamic resistance R d (I arc , l arc ) when applying the found approximation functions are given in Figure 5.Comparison of the obtained dependences of the dynamic resistance change with the ones shown in Figure 3 dependences allows us to state, that the dependences obtained on the basis approximation functions application of reproduce the nature of the change in dynamic resistance with a change in current and arc length.

Mathematical Models of Power System Elements of Electrical Arc Furnace
According to the approach described in [23,24], the calculation scheme of a three-phase short network with a furnace space (arc) is presented in the form of a six-pole (Figure 6).Such a multipole is described by a vector Equation (17): which connects the currents of the outer branches of the multipole with the average values of the potentials of the external poles at the step of numerical integration 1 ∆t Equation ( 17) is obtained from Equation ( 8), taking into account that the vector of the average values of the applied voltage at the integration step is equal to: The coefficients of Equation ( 17) are defined as a matrix: The three-phase transformer is represented by a 12-pole (see Figure 7), the outer circuits of which are the winding outputs Using an approach similar to the network model with an arc had need applied, we obtain the vector equation of a three-phase transformer: where T -vector of external branches currents; → i tr = i TA , i TB , i TC , i Ta , i Tb , i Tc T -vector of primary and secondary transformer windings currents; → v TR = v TA1 , v TB1 , v TC1 , v Ta1 , v Tb1 , v Tc1 , v TA2 , v TB2 , v TC2 , v Ta2 , v Tb2 , v Tc2 Tvector of external poles potentials.
The coefficients of Equation ( 18) are defined as follows: where R * TR = R TR 2 + L TR ∆t -matrix of step resistances of the transformer; R TR = diag R TA , R TB , R TC , R Ta , R Tb , R Tcmatrix of windings active resistances; L TR -matrix of self-and mutual-inductances of windings: A three-phase network with a choke is represented by a six-pole (Figure 8): Similarly, applying of AVIS method into the electrical circuit of a three-phase network with a choke, we obtain the Equation (19):  Based on the system of Equation ( 19), the vector equation of a three-phase network with a choke as a six-pole can be written as: where ) T -vector of potentials of external poles: ∆t -step resistances of network with choke.

Simulation Model and Calculation Algorithm of Power Supply System of Electrical Arc Furnace
The mathematical model of the arc power supply system is formed from the models of its separate multipole elements connected according to the calculation scheme shown in Figure 9.The outer branches of the multipoles are interconnected at the nodes of the electrical circuit.The connection method is mathematically described by incidence matrices, which indicate with which independent nodes of the system the poles of the multipole element are connected.For schemes, which shown in Figure 9, the incidence matrices of elements are following: for a three-phase network with a choke 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 -for the short network with furnace space.According to the algorithm described in [24], the average in integration step values of the potentials of the independent nodes in system T we can find from Equation (21): matrices of coefficients which are determined on the basis of coefficients of Equations ( 17), ( 18) and ( 20) and incidence matrices as: The system ( 21) is solved at each integration step with respect to the vector of averages in step values of potentials of the independent nodes in the system, then determine the averages value of the pole potentials of each element can be rewriten as: for three-phase network with a choke and for the short network with furnace space After determining the average in the step values of the pole potentials of each element by Equations ( 17), ( 18) and (20), can be determined the vectors of the currents of the outer branches: The calculation algorithm, which was described above, starting from finding of coefficients of vector equations of elements-multipoles (17), ( 18), (20) and finishing with definition of vectors of currents of external branches is cyclic.

Mathematical Model of the Electrode Movement System
The electric drive of the electrode movement system (see Figure 10) uses a DC motor for each phase.Its supply voltage is formed by a thyristor converter based on the mismatch signal ∆ arc , which is defined by Equation ( 22): ∆ arc = aI arc − bU arc (22) where I arc , U arc -current values of arc current and voltage; a, b-proportionality factros (setpoint).The value of the electrode movement is defined as: where K M -the transmission coefficient of the mechanical system; J-the moment of inertia; ω-the angular speed of the motor.The length of the arc l arc is determined by the sum of the displacement of the electrode le and the displacement of the charge in the melt process l dist , which are perturbations in the furnace space l arc = l e ± l dist .It should be noted that l dist changes randomly and the task of the automated control system (ACS) is to ensure that these perturbations work out.
Mathematical description of the DC motor by the method of medium voltages in the integration step of the first order is described in [30] and has the form: ) where i a0 , i a1 -the armature current of the motor at the beginning and end of the integration step; ω 0 , ω 1 -angular rotation speed at the beginning and end of the integration step; U a -supply voltage; R a -the active resistance of the armature winding; L a -inductance of the armature winding; K Φ -design factor; J-moment of inertia.In the described model, the value of the static torque of the motor, which is small compared to the dynamic torque, is neglected.The magnitude of the electrodes movement is defined by the following equation:

The Results of Process Investigation
The results of calculation of steady symmetrical modes in the arc power supply system for different values of arc length: 2 mm (in Figure 11) and 20 mm (in Figure 12) are given in Figures 11 and 12.As shown by the results of studies, increasing the length of the arc leads to an increase in the dynamic resistance of the arc, an increase in arc voltage and a decrease in current in the supply circuit.Increasing the length of the arc also leads to a deterioration of the shape of the current in the arc supply circuit and the shape of the current (Figure 12c) consumed from the mains (Figure 12d).The results of calculation of the steady asymmetric mode in the power supply system of arcs for different values of arc length in phases: phase A-2 mm, phase B-8 mm, phase C-20 mm are shown in Figure 13.In Figure 14 are shown the results of the calculation of the transient mode of arc ignition and the output to the operating mode, which is determined by the ACS for the electrodes movement which described above (Figure 10).The electrode movement system provides the formation of the arc length for which the mismatch signal defined by expression (22) becomes equal to 0. The nature of the change in the obtained dependences of current and arc voltage, current and voltage of the power supply network corresponds to the experimentally removed dependences given in the publications of other authors, in particular [18].In addition, the models created on the basis of the proposed approach can function in combination with physical objects, which createsmodel predictive control.
The simulation results completely repeat the dependences obtained using the EAF model, in which the dynamic volt-ampere characteristics of the arc are tabulated and the 4th order Runge-Kutta method is used for calculations.However, the model proposed in this article has a higher (several times) calculation speed and numerical stability (the ability to perform calculations over long time intervals).The adequacy of the mentioned model of EAF with approximated volt-ampere characteristics of arcs was confirmed using the methods of statistical data analysis, in particular the method of one-way analysis of variance [31], which based on the following equation: where x ij -summarized average value; M-factor of Barlett criteria [31].
These equations applying for data samples, which are shown in Figures 15 and 16.These dependencies correspond to the modes of the beginning of melting and the mode of completion of melting of wells.To obtain the results of mathematical modeling on the model of EAF, perturbations were applied, which in their spectral characteristics correspond to the perturbations acting in a real furnace at this stage of melting.In Figure 15a the scale for I 10 mm = 25,000 A, for time-10 mm = 0.2 s; in Figure 16-scale for I-axes-10 mm = 20,000 A, for time-10 mm = 0.2 s.
As can be seen from results in Table 1, the values of the relevant criteria obtained for different modes of operation of the EAF are less than the corresponding tabular values taken for 5%-the level of significance, which allowed us to conclude that the created system model adequately reflects the processes taking place in a real life project.It should be noted that the speed of the EAF model created using the method of medium voltage at the integration step and the use of dynamic arc resistance approximation is almost two orders of magnitude higher than the model with tabular volt-ampere characteristics and 4th order Runge-Kutta method, as well as the fact that the system retains numerical stability over long intervals of calculations.
In further research authors plan to analyze the influence of different control systems (single-circuit-when implementing different laws of electrode movement control and double-circuit-with additional high-speed inductor control circuit) on the characteristics of the power supply network, and to synthesize new, effective laws of arc electrode movement control of EAF using fractional and intelligent controllers for single-and double-circuit control systems.

Conclusions
Using the method of medium voltages at the step of numerical integration a mathematical model of the arc power system is developed.Such a model of the arc power system with electromechanical electrode movement system allows one to study transient and steady-state modes taking into account the nonlinearity of the Volt-Ampere characteristics of the arc and load asymmetry.The application of this method provides high model performance and numerical stability, which allows it to be used for real-time testing of electrode movement control systems using "hardware-in-the-loop" technology.The adequacy of the developed model is confirmed by comparison with the results of the experiment on real EAF.
The proposed approximation of the dynamic arc resistance curve using a combination of Mexican hat wavelet and the derivative function th(x) provides a fairly accurate reproduction of the dynamic arc resistance obtained on the basis of experimental volt-ampere characteristics for different currents and arc lengths.

Figure 1 .
Figure 1. Circuit of an electric branch with inductance, active resistance and arc.

Figure 3 .
Figure 3. Dependences of dynamic resistance change of R d (I arc , l arc ) for different curves of dynamic Volt-Ampere characteristic: (a)-for current increasing, (b)-for current decreasing.

Figure 4 .
Figure 4. Dependences of the change R d , obtained on the basis of experimental data and when approximating the curves of the volt-ampere characteristic using (a) hyperbolic tangent functions; (b) modified Mexican hat wavelet.

Figure 5 .
Figure 5. Approximate dependences of dynamic arc resistance changing as a function of current and arc length: (a) for current increasing, (b) for current decreasing.

Figure 6 .
Figure 6.Calculation scheme of a short network with a furnace space as a six-pole.

Figure 7 .
Figure 7. Calculation scheme of a three-phase transformer.

1 ∆t t 0 +∆t t 0 u 1 ∆t t 0 +∆t t 0 u 1 ∆t t 0 +∆t t 0 u 1 ∆te
where R NA , R NB , R NC -resistance of phases; L NA , L NB , L NC inductance of phases in choke; U N A = N A dt, U NB = NB dt, U NC = NC dt-average of a branch voltages at the integration step; E N A = NC dt-average values of EMF in integration step.

Figure 8 .
Figure 8. Calculation scheme of a three-phase network with a choke.

Figure 9 .
Figure 9. Calculation scheme of arc power supply system.

Figure 10 .
Figure 10.Structural-functional diagram of the electrode movement system.

Figure 11 .
Figure 11.Dynamic arc resistance (a), instantaneous values of arc voltages (b), instantaneous values of arc currents (c), instantaneous values of mains currents (d) in steady state for arc length 2 mm.

Figure 12 .
Figure 12.Dynamic arc resistance (a), instantaneous values of arc voltages (b), instantaneous values of arc currents (c), instantaneous values of mains currents (d) in steady state for arc length 20 mm.

Figure 13 .
Figure 13.Dynamic arc resistance (a), instantaneous values of arc voltages (b), instantaneous values of arc currents (c), instantaneous values of mains currents (d) in steady asymmetric mode for arc length: phase A-2 mm, phase B-8 mm, phase C-20 mm.

Figure 14 .
Figure 14.Arc length (a), current value of arc current (b), current value of arc voltage (c), dynamic arc resistance (d) in transient electrode movement mode.

Figure 15 .
Figure 15.The dependences of the phase B current change for the melting start stage obtained experimentally (a) and with the help of a model with tabular volt-ampere characteristics and the application of the 4-th order Runge-Kutta method (b).

Figure 16 .
Figure 16.The dependences of the phase B current change for the stage of melting completion obtained experimentally (a) and using a model with tabular volt-ampere characteristics and the application of the 4-th order Runge-Kutta method (b).

Table 1 .
Statistical data analysis.