VIEWS: 3 PAGES: 9 POSTED ON: 5/17/2011
MODELLING OF CONTINUAL INDUCTION HARDENING IN QUASI-COUPLED FORMULATION Jerzy Barglik1, Ivo Doležel2, Pavel Karban3 and Bohuš Ulrych3 1 Silesian University of Technology, Department of Electrotechnology, Krasinskiego 8, 40-019 Katowice, Poland 2 Academy of Sciences of the Czech Republic, Institute of Electrical Engineering, Dolejškova 5, 182 02 Praha 8, Czech Republic 3 University of West Bohemia, Faculty of Electrical Engineering, Sady Pětatřicátníků 14, 306 14 Plzeň, Czech Republic ABSTRACT. Continual induction hardening of steel workpieces belongs to modern and en- vironmentally friendly industrial technologies. The body is heated by an inductor moving at a low velocity along it and the heated part is then immediately cooled by spraying quenchant. Modelling of this process is an uneasy business because of two reasons. The first one is the presence of the moving inductor and the second one problem with determining the coefficient characterizing the convection heat transfer between the body and a mixture of air and spray. On the other hand, the results must be as accurate as possible in order to predict the resultant hardness. The paper presents a methodology where the induction heating (time and space variable distribution of the electromagnetic and temperature fields) is modelled in a slightly modified standard manner while the process of cooling by an approximate theoretically em- pirical procedure based on experimental data. The theoretical analysis of the task is supple- mented with an illustrative example and discussion of the results. Computations are per- formed by combination of a professional code and a number of single-purpose user programs developed and written by the authors. INTRODUCTION Continual induction hardening represents a technological process that is frequently used in various industrial applications. The processed workpiece is heated by an inductor carrying harmonic current that moves in common with a sprayer at a low velocity along it and the heated part of the body is immediately cooled by spraying quenchant (water, oil, polymer- based liquids etc.). The result is harder, but more fragile martensite structure. The scheme of the process is shown in Fig. 1. The final hardness of the surface layers of the workpiece is generally a complex function of velocity (or, equivalently, time) of its heating and cooling. This is partly explained in Fig. 2 containing the unbalanced phase diagram of a typical carbon steel for a given austenitizing temperature Ac3 . First, the part of the body to be hardened has to be heated over the temperature Ac3 in or- der to reach uniform austenite structure (domain A) and then cooled. The decreasing lines in Fig. 2 starting from the temperature Ac3 and ending at the final temperature TF represent various regimes of cooling. The first two lines on the left indicate, for example, relatively fast cooling leading to complete hardening of the processed material; the resultant structure con- taining practically pure martensite is characterized by high hardness given (in Vickers) within the corresponding circles. The lines on the right, on the other hand, express slower cooling resulting in incomplete hardening (the final structure contains not only martensite M, but also undesired bainite, pearlite or ferrite whose zones in the diagram are denoted by letters B, P and F, respectively). moving heated inductor part v moving sprayer spraying quenchant cooled part hardened steel workpiece Figure 1. The investigated arrangement T (°C) 1000 A 900 Ac3 800 P F 700 600 673 409 500 300 258 400 B 300 Ms M 200 100 HV = 765 762 742 505 432 279 243 TF 0 1 10 102 103 104 105 t (s) Figure 2. Unbalanced phase diagram of Polish steel NZ3 (0.6% C, 1.0% Si, 1.0% Cr, 2.0% W and 0.2% V), hardness given in Vickers, Ac3 = 880 °C, TF = 100 °C The task represents a coupled electromagnetic-thermo-metallurgical problem characterized by mutual interaction of several physical fields with time-dependent boundary conditions that can be (particularly in the case of temperature field) determined only indirectly. Moreover, the material parameters are temperature-dependent functions that must be found experimentally. The first mathematical and corresponding computer model for solution of the problem was suggested by the authors in [1]. But even when being relatively sophisticated, the methodology suffered from several simplifications enforced by lack of experimental data and possibilities of available software. Mentioned can be • weakly coupled formulation of the problem (electromagnetic field was supposed to be independent of the temperature field), • step change of magnetic permeability when the temperature reached the Curie point, • temperature Ac3 was considered constant (history of heating was not taken into account), • step variation of the convection heat transfer coefficient from the surfaces of the heated and cooled parts to ambient medium, • rough mesh used for thermal computations (only several hundred elements). Since then, however, the methodology of solution was substantially improved, together with our own procedures (that allow, nowadays, finite-element analysis of an arbitrary axi- symmetric problem of this kind in the quasi-coupled formulation). The paper describes the algorithm and illustrates it on an example. DESCRIPTION OF THE TECHNICAL PROBLEM The inductor (Fig. 1) carrying harmonic current I ext of density J ext and frequency f moves at a velocity v along a steel mandrel in common with a sprayer with cooling medium under pressure that assures consequent fast cooling of the heated part. The first task is to recommend the parameters of the inductor (geometry, current, fre- quency and velocity v of its motion) that provide sufficient temperature rise of the heated part (somewhat higher than Ac3 ). The second task is to find such pressure (or velocity) of the spraying quenchant that the convection of heat from the body to ambient mixture guarantees sufficiently fast cooling and required surface hardness. While the first task is processed theoretically and realized by a suitable numerical algo- rithm, the second task is processed by an experimentally-empirical procedure. This procedure is, nevertheless, realized again in the numerical manner. MATHEMATICAL MODEL AND NUMERICAL ALGORITHM OF ITS SOLUTION The mathematical model consists of two partial differential equations describing the distri- bution of the electromagnetic and temperature fields with temperature-dependent coefficients. Electromagnetic field The investigated arrangement is considered fully axisymmetric, so that the definition area Ω in Fig. 3 is rectangular (ABCDA) with artificial boundary BCDA sufficiently distant from the workpiece and inductor. This area consists of domains Ω 1 (mandrel), Ω2 (moving induc- tor) and Ω3 (air). Relative permeability μ r of the workpiece is supposed to be only a tem- perature-dependent function, which simplifies the situation. The field quantities may be now replaced by their phasors and the field distribution is described by the Helmholtz equation [2] rot rot A + j ⋅ ωγμ A = μ J ext (1) where A denotes the phasor of the vector potential, γ the electrical conductivity, μ the magnetic permeability and ω = 2πf . Each domain is characterised by specific values of phasor J ext , γ and μ , the last two parameters being supposed to be functions only of the temperature T . In fact, phasors of the current density as well as the magnetic vector potential have only one nonzero tangential component J ext,ϕ and Aϕ , respectively. The Dirichlet con- dition along the artificial boundary BCDA reads Aϕ = 0 . This boundary has to be taken far enough from the arrangement in order to avoid any influence on electromagnetic field distri- bution in it. Temperature field Its definition area is represented by the longitudinal axial cut through the workpiece (KLMNK in Fig. 3) and its distribution is given by the Fourier-Kirchhoff equation [3] ∂T div ( λ ⋅ grad T ) = ρ c ⋅ − wJa , (2) ∂t where λ denotes the specific heat conduction of the material, ρ its specific mass, c its spe- cific heat and wJa the specific average Joule losses 2 wJa = γω 2 A . (3) B C air Ω3 z h2 r rb E L M v moving inductor with sprayer Ω2 l Z G hardened workpiece Ω1 K N ra S h rc w h1 r1 A D Figure 3. Definition area of the problem: S – starting position of the inductor, G – general position of the inductor, E – end position of the inductor The boundary condition for t > 0 along line KL reads ∂T / ∂r = 0 while along line LMNK −λ ⋅ ∂T / ∂n = α (T − T0 ) . The influence of radiation above the sprayer is neglected at this stage. Distribution of the convection heat transfer coefficient and temperature of the ambient medium depends on the instantaneous position of the inductor with the sprayer and is depicted in Fig. 4. While α w (coefficient for heat transfer from the workpiece to cooling medium) is usually much greater than α a (analogous coefficient from the workpiece to ambient air), tem- perature T0w of the coolant is, on the contrary, lower than temperature of air T0a . The initial condition ( t = 0 ) for the temperature field within the whole region Ω is Tstart . Now it is necessary to determine the coefficient α w characterizing the convection heat transfer between the workpiece and spraying quenchant. As its theoretical investigation is extremely difficult, it was determined by an indirect approximate theoretically empirical method based on experimental results. This method suggested by the authors is described in details in [4]. In fact, α w is a function of the surface temperature, but the method provides only its mean value α wm . starting general position S position G α T0 α T0 surface of the heated workpiece αa T0a αa T0a G αwm T0w S Figure 4. Distribution of α and T0 along the mandrel for starting (S) and general (G) positions of the inductor Distribution of hardness Dependence of hardness on the time of cooling tc from temperature Ac3 to TF was re- drawn from Fig. 2 to Fig. 7b. This dependence was built into the code. When the highest tem- perature of any point does not reach Ac3 , computation of hardness is not realized there. Numerical algorithm 1. Choice of a time step Δt1 for which the inductor with the sprayer are in one position (their motion is not supposed to be continuous, but discrete, by small shifts Δl ). The time step Δt1 = 2 s (corresponding to change of the position of the system inductor-sprayer 1 mm) decisive for changes of electromagnetic field proved to be small enough to achieve re- quired accuracy of the calculations. 2. Computation of electromagnetic field, eddy currents and corresponding Joule losses in the heated body for the actual position of the inductor. 3. Determination of the time-variable boundary conditions for the temperature field by means of the described theoretically-empirical procedure. The procedure makes it possible to take into account average value of coefficient of convection heat transfer α wm during the process of cooling. Computation of the time evolution of nonstationary temperature field was realized with time steps Δt2 satisfying the condition of numerical stability of so- lution. Automatic change of material parameters according to actual temperature in indi- vidual cells. Creation of the database containing history of cooling of cells in the surface layers of the cooled part of the workpiece. 4. Shift of the inductor by value Δl = v ⋅ Δt1 and return to point 2. Points 2–4 are repeated until the whole workpiece is processed. 5. Simultaneous calculation of the resultant field of hardness. ILLUSTRATIVE EXAMPLE The methodology is illustrated on an example of induction hardening of a steel mandrel (steel NZ3 whose unbalanced diagram is in Fig. 1). The most important input data • Several preliminary computations showed that the safe artificial boundary BCDA is char- acterized by the following dimensions: h1 = h2 = 0.15 m, r1 = 0.2 m. Its enlargement leads to small changes in results not exceeding about 2%. • Heated workpiece from steel NZ3: ra = 0.012 m, rb = 0.0205 m, l = 0.610 m. Material parameters (and eventually their temperature dependencies) are depicted in Figs. 5, 6 and 7. Preliminary computations indicated that the average magnetic flux density in the sur- face layers of the mandrel is about 0.5 T (Fig. 5a). That is why we used corresponding relative permeability μ r = 2600 determined for interval B ∈ 0,1.5 T. Its temperature dependence is in Fig. 5. • Coil: rc = 0.033 m, w = 0.017 m, h = 0.03 m. Other parameters: J ext = 107 A/m2, f = 10 kHz, v = 0.0005 m/s. In fact, the real coil made from water-cooled copper tube was replaced by a thin copper conductor. The main reason of this simplification was to sub- stantially reduce the time of computations. • The principal temperature conditions: Tstart = 20 °C, T0a ( LMZ ) = T0a ( KN ) = 20 °C, T0w ( NZ ) = 10 °C (MZ – above the sprayer, NZ – below the sprayer), α a ( LMZ ) = 20 W/m K, α wm ( NZ ) = 500 W/m K, α a ( KN ) = 5 W/m K. 2 2 2 3000 2000 average value μr a) 1000 0 0 0,5 1 1,5 B (T) Figure 5. Dependencies μ r ( B ) for T = 20 °C – a) and μ r (T ) for B = 0.5 T – b) (steel NZ3) Selected results and their discussion Computations were carried out on a computer CPU 2.6 GHz with memory RAM 0.5 GB. The average time of computation of one case was about four hours. Electromagnetic field was calculated on a mesh with about 200000 elements (ABCDA, see Fig. 3) while the mesh for consequent temperature field consisted of more than 10000 elements (KLMNK, see Fig. 3). The area was remeshed at each position of the system inductor-sprayer (their number was about 80). Geometrical convergence of the results was tested carefully with respect to both the artificial boundary and density of the discretization mesh. 3,5 33 λ(W/mK) 31 γ(MS/m) 2,5 b) 29 1,5 a) 27 0,5 25 0 150 300 450 600 750 900 0 150 300 450 600 750 900 T (°C) T (°C) Figure 6. Dependencies γ = γ (T ) – a) and λ = λ (T ) – b) for steel NZ3 7 800 hardness (HV) ρc(MJ/m 3K) 700 6 a) 600 5 500 b) 4 400 3 300 0 150 300 450 600 750 900 0 500 1000 1500 2000 T (°C) tc(s) Figure 7. Temperature dependence of specific heat – a) and dependence of hardness on the time of cooling – b) for steel NZ3 The most important results are depicted in Figs. 8, 9, 10. Fig. 8 shows the smoothed time evolution of temperature at selected point of the mandrel surface. Fig. 9 depicts, on the other hand, distribution of temperature along the mandrel for several positions of the inductor. The curves for t = 1400 s and t = 1600 s were calculated for cooling, when the inductor reached the end of the mandrel and the current was switched off. There are small discrepancies between Figs. 8 and 9 caused by local numerical instabili- ties. These differences not exceeding several percents could surely be reduced by further re- finement of the mesh elements and time steps, but at the expense of the computational costs. Moreover, this inaccuracy is comparable with inaccuracy of some input data (or even lower) and it is disputable whether eventual improvement in this direction would bring more accurate results. Finally Fig. 10 shows the distribution of hardness along the surface NM of the mandrel. Its variations follow from different times of cooling of its individual parts. In fact, however, the difference between minimal and maximal values of hardness is 16 HV, which does not exceed about 2.5% from the average value. CONCLUSIONS Accuracy of computations depends on the quality of the mathematical model, particularly on accuracy of experimentally obtained temperature-dependent properties of present media. Very important is also knowledge of hardness as a function of the velocity of cooling. The computations themselves take a relatively long time (a couple of hours) due to both complex- ity of the task and verification of the geometrical convergence of results. These are in a rela- tively good accordance with experiments [4]. Figure 8: Time evolution of temperature at selected points along the surface of the mandrel 1000 900 800 700 temperature (deg) (°C) 600 500 207 s 400 720 s 300 1100 s 200 1400 s 100 1600 s 0 0 0,1 0,2 0,3 0,4 0,5 0,6 position of the inductor with respect to the mandrel (m) position of the inductor (m) Figure 9: Distribution of temperature along the mandrel for the positions of the inductor given by the presented times Next work in the field will be aimed at more detailed evaluation of the history of heating (because the value of Ac3 depends on its velocity) and eventual influence of radiation in sys- tems with moving sprayers. Some progress in acceleration of computation algorithms is ex- pected as well. ACKNOWLEDGEMENT Financial support of the Czech Ministry of Education (Research Plans VZ MSM 232200008 and 212300016) and Polish Scientific Research Committee (grant No. 4T08C 04823) is highly acknowledged. As well the paper was prepared in a framework of multi- government bilateral collaboration between Poland and the Czech Republic. 790 770 750 hardness (HV) 730 710 690 670 650 0 0,1 0,2 0,3 0,4 0,5 0,6 axial distance (m) Figure 10: Distribution of hardness along the surface NM of the mandrel REFERENCES [1] Doležel, I., Barglik, J., Ulrych, B. (2003). Continual induction hardening of axisym- metric bodies. Proc. of IC JAPMED’2003, 79–82. [2] Chari, M.V.K., Salon, S.J. (2000). Numerical methods in electromagnetism. Academic Press. [3] Holman, J.P. (2002). Heat transfer. McGrawHill. [4] Barglik, J., Doležel, I., Ulrych, B. (2003). Determination of convection heat transfer coefficient at continual induction heating of axisymmetric bodies. Scientific Bulletin of Lodz Technical University, No. 101, 9–16.