Thermal analysis of moving induction heating of a hollow cylinder with subsequent spray cooling: Effect of velocity, initial position of coil, and geometry

Applied Mathematical Modelling 36 (2012) 4304–4323

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Thermal analysis of moving induction heating of a hollow cylinder with subsequent spray cooling: Effect of velocity, initial position of coil, and geometry H. Shokouhmand, S. Ghaffari ⇑ School of Mechanical Engineering, College of Engineering, University of Tehran, Tehran, Iran

a r t i c l e

i n f o

Article history: Received 12 March 2011 Received in revised form 13 November 2011 Accepted 15 November 2011 Available online 23 November 2011 Keywords: Moving induction heating Moving heat source Finite-element method Spray cooling Moving boundary conditions Temperature distribution

a b s t r a c t In this paper, temperature analysis of the complete process of moving induction heat treatment is performed using numerical methods. A non-linear and transient magneto-thermal coupled problem with a moving coil which is considered as moving heat source, is investigated by an efficient finite-element procedure. A vertical hollow circular cylinder is heated by the moving coil at a given velocity along it, and the heated parts then quenched by a moving water–air spray. The effects of natural convection with air on the both inner and outer surfaces of cylinder, and also radiation of outer surface of cylinder with ambient are taken into account. For quenching of work-piece, a specific kind of atomized spray cooling which utilizes a mixture of water and air with different mass fractions is used. This procedure includes moving boundary conditions, temperature-dependent properties, and change in magnetic permeability of specified alloy at the Curie temperature. Obtained numerical results have been verified by comparison with analytical solutions using Green’s function methods. Also, the effect of velocity, initial position of inductor and inner to outer radius ratio on temperature distribution are investigated. Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction Induction heating is a non contact method of precisely and accurately heating conductive materials and widely used in today’s industrial operations including surface hardening, melting, welding, forging and other similar applications in the automotive, aerospace and other engineering sectors. Induction heating process has a number of benefits over traditional methods such as no naked flame or atmosphere safety issues, repeatability, high efficiency, reduced scale, rapid focused heat, and reduced operator skill levels [1]. In the stationary induction heat treating the whole area to be hardened is heated at once and no relative motion happens between the coils and the work-piece. In the moving induction heat treatment relative motion takes place between the work-piece and coil during the process, as with long shafts, where the area to be hardened is progressively heated and quenched by a moving induction coil followed by a quenching ring or a cooling spray. Moving induction heat treatment is used extensively in the production of shaft type. There are some demerits with the stationary induction heat treatment that cause moving induction heat treatment to be more efficient. The coil design can be an extremely complex process. Often the use of ferrite or laminated loading materials is required to influence the magnetic field concentrations in given areas

⇑ Corresponding author. Tel.: +98 9125031913. E-mail addresses: [email protected] (H. Shokouhmand), [email protected] (S. Ghaffari). 0307-904X/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2011.11.058

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

4305

Nomenclature A B CP E f g_ Gr hc IS JS Je k L Nu Pr q00 Ri RO Ra Re T v ZO ZC

magnetic vector potential magnetic flux density specific heat capacity electric field intensity frequency power of heat source Grashof number convection heat transfer coefficient current in inductor density of current in inductor eddy current density thermal conductivity length of cylinder Nusselt number Prandtl number heat flux inner radius of cylinder outer radius of cylinder Rayleigh number Reynolds number temperature velocity of moving coil and spray initial position of coil position of coil

Greek symbols a thermal diffusivity e emissivity r electric conductivity rs Stefan–Boltzmann constant q density of material l magnetic permeability of steel lo magnetic permeability of vacuum x angular frequency v the ratio of water mass-flow rate to air mass-flow rate Subscripts a air i inlet condition j jet stream m mean or mixture mf mean film ms mean surface max related to maximum heat flux sat saturation 0 stagnation point 1 ambient condition Superscript ⁄ conjugate complex 0 related to phasor

thereby to refine the heat pattern produced. Another drawback is that much more power is required due to the increased surface area being heated compared with a moving induction heat treatment approach. Steels are particularly suitable for heat treatment, and they are heat treated for hardening, softening or material modification. For hardening of steels the body to be hardened is gradually heated above the austenitizing temperature; and when suddenly quenched, the Martensite is formed. This is a very strong and brittle structure. When slowly quenched it would form Austenite and Pearlite which is a partly hard and partly soft structure. When the cooling rate is extremely slow then

4306

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

it would be mostly Pearlite which is extremely soft. This shows that the faster cooling leads to more hardness of steel after heat treating. Since induction heating is a complex process, its obtained results usually based on experiments. But these experiments are not economical and also do not provide a great deal of information about the process. From engineering perspective, it is important to investigate temperature distribution of body during the induction heating process and its subsequent quenching. Because, the temperature fields predicted by the model are used both to predict the resulting microstructure using continuous cooling transformation (CCT) diagrams and as an input for the thermal-stress model to predict the occurrence of quench cracks, so modeling of induction heating process by numerical methods is useful with general purpose of commercial programs. In the past decade, considerable efforts have been devoted to analyze induction heat treatment. Several numerical and analytical works have been existed for stationary induction heat treatment [2–4]. In the most advanced one, Wang et al. [5], obtained analytical and numerical solution for stationary induction heating in a cylinder with assuming temperaturedependent properties using finite-element methods. Sadeghipour et al. [6], employed a general-purpose finite-element program to simulate and analyze induction heating process of steel. The finite element results were evaluated and compared with the experimental results. Kranjc et al. [7], investigated stationary induction heating treatment of cylindrically shaped steel work piece numerically and experimentally with three different heating protocols. Temperature-dependent and temperature-independent steel material properties were considered and their impact on simulation results was evaluated. Although papers and books on stationary induction heating are numerous, there are a limited number of publications pertaining to moving induction heat treatment. Wang et al. [8], developed their previous problem and solved it with moving induction heat treatment condition. The finite-element procedure incorporates temperature-dependent phase transformation, change in magnetic permeability at the curie temperature, and temperature, residual stress and microstructure distribution in cylinder have been obtained. Quenching of heated cylinder by both moving cooling ring and a stationary liquid bath, has been analyzed. Yun et al. [9], simulated moving induction heating process for a moving coil around a steel plate, and finite-element results were verified by experimental results. Dolezel et al. [10] presented a methodology of computer modeling of moving induction hardening of axi-symmetric bodies. In this publication the heated part is cooled by moving water jet, and the process of cooling is modeled by a theoretically empirical algorithm based on experimental data and providing the value of the convective heat transfer coefficient. When work-piece is heated to a specified temperature, it must be quenched. The quality of cooling is a specified problem by itself, and is a complex process. As said before, the hardness of steel after heat treatment depends on the time of cooling, so cooling of heated body is a crucial part in heat treatment, and should be studied carefully. Spray cooling using a mixture of water and air was found useful in controlling the cooling rate of hot medium-carbon steel bars. With spray cooling, one can raise or lower the cooling rate by increasing or decreasing the amount of liquid in the mixture. The atomized spray consists of small liquid droplets in a conical jet of air. In this application, water is the liquid of choice because of its low cost. Papers for spray-cooled surface at a high temperature are limited. Because, at the high temperatures the behavior of flow becomes multi-phase, and the heat-transfer relations are not well established as they are in single-phase flow. Thomas et al. [11], presented finite-element modeling and experimental verification of spray cooling process in a steel cylinder from an initial temperature of 1273 K. They also demonstrated the prediction of quench cracks with commercial finite-element analysis (FEA) using available information on a complex heat-transfer phenomenon like spray cooling. The temperature fields predicted by the model are used as an input for the thermal-stress model to predict the occurrence of quench cracks. The review of the previous works showed that some studies were conducted on analysis of magneto-thermal problem to investigate temperature distribution during the heating process and they did not involve in subsequent quenching process sufficiently. Also some studies were conducted on the cooling process of a work-piece and did not point to process of heating. However, according to the authors’ knowledge, there is no study of thermal analysis of moving induction heating with prefect investigation of cooling process. Also there is no study of hollow cylinder in induction heat treatment problems. Hollow cylinder can be used as shaft or tube. When it used as shaft, it has less weight rather than solid shafts with the same diameter, and leads to less energy consumption. The present publication is devoted to study of thermal analysis of moving induction heat treatment and subsequent cooling process by moving water–air spray for a hollow cylinder. The coupled magnetic and thermal problem in a vertical hollow circular cylinder must be solved regarding to the mentioned moving boundary conditions, because the material properties in the induction heating depend on the temperature. The effects of natural convection with air on the both internal and external surfaces of cylinder, and also radiation of outer surface of cylinder with ambient are taken into account.

2. Physical model The basic arrangement of moving induction heating process in longitudinal section of a vertical hollow cylinder is depicted in Fig. 1. The current in the inductor IS of density JS is harmonic with frequency of f. The inductor moves at a velocity of v along the cylinder. The water–air spray cooling moves with a specific distance to the inductor at the same velocity with it to assure consequent fast cooling of the heated part.

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

4307

Fig. 1. The basic arrangement of problem.

Fig. 2. Domain of solution.

This problem consists of two main steps.The first one is to investigate the induction heating and its relevant parameters that provide sufficient temperature of the heated part (higher than austenitizing temperature); and, the second step is to find such convective heat transfer coefficient between the cooled body and cooling medium that secures the required hardness. 3. Formulation of problem As said before, in this problem a vertical hollow circular cylinder arrangement is investigated. Since this investigated arrangement may be considered axi-symmetric, the domain of solution can be modeled as shown in Fig. 2. 3.1. Analysis of electro-magnetic field To evaluate the heat sources generated within the material during induction heating, it is first necessary to obtain the magnetic vector potential in the work-piece relevant to a current flowing through the induction coil. For this purpose, the

4308

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

Maxwell’s electro-magnetic field equations must be solved with the appropriate boundary conditions. The magnetic field can be expressed in terms of the magnetic vector potential, A0 , which is defined by equation:

B ¼ r  A0 ;

ð1Þ

where B is the magnetic flux density. So, the Maxwell’s equation in terms of magnetic vector potential holds as follows, [8]:

r



1

l



r  A0 ¼ r

@A0 þ J 0S ; @t

ð2Þ

where J 0S is the source current density, l is the permeability, r is the electric conductivity and for a sinusoidal source current with a angular frequency of x = 2pf, we get

J 0S ¼ J S eixt :

ð3Þ

0

A is sinusoidal with a angular frequency of x too, so that

A0 ¼ A00 eiðxtþuÞ ¼ Aeixt :

ð4Þ

Since the range of frequencies used for our induction heating problem is low, a quasi-static approximation can be used [8]. Then,

r  A ¼ 0:

ð5Þ

Finally, from Eqs. (3)–(5), (2) can be expressed as:

1

l

  1

r2 A  ixrA  r

l

 ðr  AÞ þ J S ¼ 0:

ð6Þ

When dealing with axi-symmetrical configurations, the phasors of the current density JS as well as the magnetic vector potential A can be expressed in cylindrical co-ordinates with a only non-zero component in the h direction, i.e. JSh and Ah. Therefore, by dropping the h subscript, A and JS are used for Ah and JSh. Finally, Eq. (6) for our problem in cylindrical co-ordinates becomes:

1 @ 2 A 1 @A @ 2 A A þ þ  l @r2 r @r @z2 r2

!  ixrA þ

  @ð1=lÞ 1 @rA @ð1=lÞ @A þ þ J S ¼ 0: @r r @r @z @z

ð7Þ

3.1.1. The boundary conditions The boundary condition along the artificial boundary EFGHIJE sufficiently distant from the cylinder and inductor shown in Fig. 2 [10], holds:

A ¼ 0:

ð8Þ

For the arrangement shown in Fig. 2 and for a delta-function coil of strength JS, the other boundary conditions are

A1 ðRi ; zÞ ¼ A2 ðRi ; zÞ; 1 @ 1 @ A ðR ; zÞ ¼ A ðR ; zÞ; l0 @r 1 i l @r 2 i A2 ðRO ; zÞ ¼ A3 ðRO ; zÞ; 1 @ 1 @ A ðR ; zÞ ¼ A ðR ; zÞ; l @r 2 O l0 @r 3 O

ð9Þ ð10Þ ð11Þ ð12Þ

where the subscripts represent the regions shown in Fig. 2, and l0 is the permeability of air. Magnetic vector potential can be calculated from the solution of Eq. (7) with applying above boundary conditions using finite-element methods, and the eddy current can be evaluated by using magnetic vector potential [6], and is given by

J e ¼ rE ¼ ixrA;

ð13Þ

where E is the electric field intensity and Je is the eddy current density. Therefore heat source can be expressed as:

g_ ¼

½ReðJ e Þ2

r

;

ð14Þ

where Re(Je) is the real part of Je. For this case where the source current density is assumed to be time harmonic, Eq. (14) for the average time can be written as:

g_ ¼

1 2 x rAA ; 2

where A⁄ is the complex conjugate of the magnetic vector potential A.

ð15Þ

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

4309

3.2. Analysis of thermal field The calculation of the temperature distribution in cylinder during the moving induction heating and subsequent quenching requires a solution of the heat equation with convection and radiation boundary conditions. The transient heat conduction equation for a solid with a heat source of strength g_ is expressed as [12]:

qC p

@T _ ¼ r  ðkrTÞ þ g: @t

ð16Þ

Since our problem is axi-symmetric, Eq. (16) in cylindrical co-ordinates can be written as:

qC p

    @T 1 @ @T @ @T _ þ þ g; ¼ kr k @t r @r @r @z @z

ð17Þ

where k is thermal conductivity, q is density, Cp is specific heat. 3.2.1. The boundary conditions Ceramic rings were attached to both ends of the cylinder to eliminate heat loss and they works as insulators. Therefore, there is no heat transfer from lines AB and CD. Therefore the relevant boundary condition is as follows:

z¼0: z¼L:

@Tðr; z; tÞ ¼ 0; @z @Tðr; z; tÞ k ¼ 0: @z

ð18Þ

k

ð19Þ

3.2.1.1. Internal surface (line AD).

r ¼ Ri :

k

@TðRi ; z; tÞ ¼ ðhc Þ1 ½TðRi ; z; tÞ  T 1 ; @r

ð20Þ

where T1 is ambient temperature and (hc)1 is natural convection heat transfer coefficient for inner surface as shown in Fig. 3. Considering that natural convection for inner and outer surfaces of cylinder has much less effect on temperature distribution than heating process and other mechanisms of cooling, It is necessary to approximate natural convection heat transfer coefficient for inner and outer surfaces and use a experimental correlation of similar case from previous works. Al-Arabi et al. [13] investigated experimentally natural convection heat transfer from the internal surface of vertical tubes to air. The results obtained have been correlated by dimensionless groups as follows:

Fig. 3. Presentation of boundary conditions in problem.

4310

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

Num ¼

1:11 ðGrm PrÞ0:25 ; ½1 þ 0:05ðT ms  T i Þ=T i 

ð21Þ

where Tms is the mean surface temperature and calculated as follows:

T ms ¼

1 L

Z

L

TðRi ; z; tÞdx:

ð22Þ

0

Ti is inlet temperature and is equals to T1, ambient temperature. Grm is mean Grashof number and Pr is Prandtl number of air. The physical properties in this work are evaluated at mean film temperature as given by:

T mf ¼ ðT ms þ T i Þ=2:

ð23Þ

Therefore, from Eq. (20) natural convection heat transfer coefficient for inner surface can be evaluated as follows:

ðhc Þ1 ¼

Num  k L

ð24Þ

where k is thermal conductivity of the air and L is length of the cylinder. 3.2.1.2. External surface (line BC). Since sprayer moves along cylinder, in each time step some points of outer surface of cylinder are under condition of cooling with water–air spray as shown in Fig. 3. The spray engulfs a limited region that its domain can cover. Hence, outer surface of cylinder can be divided into two different parts in each time step: the region exposed to spray cooling, and the dry region that is not influenced by spray. Forced convection heat transfer coefficient due to spray cooling is much greater than natural convection heat transfer coefficient with air. As a result for the region exposed to spray cooling the effect of natural convection can be neglected. Therefore, boundary conditions for this region are:

@TðRO ; z; tÞ ¼ ðhc Þ2 ½TðRO ; z; tÞ  T j  þ hr ½TðRO ; z; tÞ  T 1 ; @r 2 2 hr ¼ rs eðT ðRO ; z; tÞ þ T 1 ÞðTðRO ; z; tÞ þ T 1 Þ;

r ¼ RO :

k

ð25Þ ð26Þ

where (hc)2 is forced convection heat transfer coefficient for outer surface due to spray cooling that is discussed later. Tj is the jet stream temperature, hr is radiation heat transfer coefficient that is defined by Eq. (26). rs is the Stefan–Boltzmann constant that equals 5:67  108 W m2 K4 and e is the material emissivity coefficient. The cylinder surface region not exposed to cooling spray were subjected to natural convection with air during the heating and quenching process. Therefore, boundary conditions for this region are:

r ¼ RO :

k

@TðRO ; z; tÞ ¼ ðhc Þ3 ½TðRO ; z; tÞ  T j  þ hr ½TðRO ; z; tÞ  T 1 ; @r

ð27Þ

where (hc)3 is natural convection heat transfer coefficient for outer surface with adjacent air as shown in Fig. 3, and hr is radiation heat transfer coefficient that is defined by Eq. (26). Wei Qi et al. [14] investigated experimentally natural convection heat transfer from the inside surfaces of vertical tubes to air. The results obtained have been correlated by dimensionless groups for average Nusselt number follows: 0:385 0:524 Num ¼ 0:2789 Ra0:32 K ; m H

ð28Þ

where Ram is the average Rayleigh number, also H and K are:

L ; RO  Ri RO : K¼ Ri

H¼

ð29Þ ð30Þ

3.2.2. The initial conditions The initial condition for the temperature field within the whole regions of cylinder is:

Tðr; z; 0Þ ¼ T 0 ;

ð31Þ

where T0 denotes the initial temperature of system, it means in this problem, T0 = T1. 4. Evaluation of forced convection heat transfer coefficient due to spray cooling Since spray cooling process for heated part occurs at high temperatures, the heat transfer establishes in multi-phase flow. In addition, numerical investigation of multi-phase flow heat transfer is very complicated. Therefore, available experimental relations for spray cooling process of steel are used for obtaining forced convection heat transfer coefficient for outer surface due to spray cooling. With considering BETE 1/400 XA PR 150B manufactured in Bete Fog Nozzle Inc. [11], which internally mixes compressed air and water in an annular region within the nozzle to form a fine conical spray. Schematic view of nozzle

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

4311

is shown in Fig. 4, where for this nozzle spray angle is 14°, during cooling, the cylinder is located 41.72 cm from the nozzle, therefore with considering angle of nozzle Domain of spray on cylinder is 10 cm. By increasing or decreasing the amount of liquid in the spray, one can increase or decrease the cooling rate at the surface. With presenting v, that is defined as the ratio of the water-mass-flow rate to the air-mass-flow rate, the heat-transfer process during spray cooling can be divided into three specific regions: 4.1. Radiation-dominated region When surface temperature becomes high enough that radiant energy evaporates all of the water droplets before reaching the boundary layer and the mixture of air and vaporized water will cool the surface by convection. In this region, heat transfer from the cylinder is by radiation and convection to a mixture of air and water vapor as said before in Eq. (25). Hence, a single-phase model can predict the convection heat transfer. The heat-transfer coefficient at the stagnation point can be evaluated using equation for single-phase flow over a cylinder, [11,15]:

hR;0 ¼

    2  km 9 0:33 ; 82:81 þ 0:003014 Rem  Pr0:33 Re  Pr  3:256  10 m m m 2RO

ð32Þ

where Reynolds number Rem is defined by:

Rem ¼

qa u0 ð1 þ vÞð2R0 Þ TðRO ; z; tÞ : T sat lm

ð33Þ

In above equations km, lm and Prm are thermal conductivity, the viscosity and prandtl number of mixture, respectively. qa is density of air, u0 is maximum velocity of nozzle at the distance of cylinder, Tsat is the saturation temperature. Bird et al. [16] presented a method for obtaining thermophysical properties of mixture. Also, all properties of air and water considered at mean film temperature that equals to (T(RO, z, t) + Tj)/2, where Tj is the jet stream temperature. This condition occurs when T(RO, z, t)  Tmax is large, where Tmax is calculated by [10]:

T max ¼ 125 þ 0:33v3 :

ð34Þ

4.2. Convection-dominated region In the convection dominated region, the evaporation of droplets mainly takes place in the boundary layer and heat transfer to the droplets is by convection and water droplets reach the surface of the cylinder and boil away. This region begins when T(RO, z, t) < Tmax. The heat-transfer coefficient at the stagnation point in this region can be calculated by [15]:

hC;0 ¼

1:9  105 ka ðTðRo ; z; tÞ  T j Þ0:84 v0:75 ; 2R0

where ka is thermal conductivity of air.

Fig. 4. The schematic view of cooling spray.

ð35Þ

4312

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

4.3. Transition region A transition region is between two above regions. In this region, the droplets evaporate partially before they enter the boundary layer. The heat-transfer coefficient at the stagnation point in this region is approximated by [11]: hmax max hR;0

h

h0 ¼ hR;0 þ ðhmax  h0 Þe

TðR

O ;z;tÞT max T max T a

;

ð36Þ

where hmax is the value of hC,0 in Eq. (28) when T(Ro, z, t) equals Tmax. Consequently, forced convection heat transfer coefficient for outer surface due to spray cooling can be evaluated from heat-transfer coefficient at the stagnation point for each region as follows, [11]:

   1  4C cos 2h h ¼ h0 1 þ 0:5ð1  BÞ þ cos h  ð1  0:18Z 2 Þ: 4C 4C

ð37Þ

In Eq. (37):

B ¼ 0:385  0:0426v0 ;

ð38Þ

C ¼ cos h0 ;

ð39Þ

h0 ¼ 1:995  0:016X 0 ;

ð40Þ

where

(

v0 ¼ v; v < 4; v0 ¼ 4; v > 4:

ð41Þ

In Eq. (37), Z is distance from stagnation point along cylinder. Since domain of spray is 10 cm, we have 5 cm 6 Z 6 5 cm. Since this problem is axi-symmetric in geometry and all boundary conditions, h must be eliminated so that the problem remains in axi-symmetric case. In this problem two similar sprayer are placed symmetrically two different sides of the cylinder. As a result, For this reason the average of the heat-transfer coefficients over the range of 0 < h < p/2 is calculated by integrating Eq. (37) from h = 0 to h = p/2 and dividing by p/2, and this average value is specified over the circumference of cylinder. 5. Consideration of the effect of moving boundary in solution The calculation of the temperature field within the work-piece can be carried out in a moving coordinate system (r, k) that is attached to the moving coil, where k = z  vt. Therefore, the same mesh can be used at each time step. So, in each time step, the position of coil, position of spray and regions exposed to natural convection must be updated. In the other words, Eqs. (7) and (37) must be solved using moving coordinate system. 6. Numerical solution The governing equations with appropriate boundary conditions are solved by employing the finite-element method [17]. At first, Eq. (7) must be solved for obtaining heat source generated in the work-piece. After discretization, the finite-element formulation obtained from the governing equation for each element is as:

½Ke fAg ¼ fJ e g;

ð42Þ

where

    T  ! ZZ @N @N @N 1 ½NT ½Nr dr dz; ½K ¼ þ i x r þ r dr dz þ 2p @r @z @z l lr 2 ZZ Z 1 T @A fJ e g ¼ dS: 2pJ½NT r dr dz þ ½N @n l e

ZZ

2p



@N @r

T 

ð43Þ ð44Þ

In the formulation, the permeability, l, is assumed to be constant within a element bounded by the surface S, but can vary from one element to another. [N] is the appropriate shape function matrix and {A} is the matrix representing the nodal values of the complex vector potential for each element. The entire space is modeled by isoparametric eight-node elements, i.e. the quadratic lagrangian interpolation is used for variables over an element. After solution of this equation the heat source in work-piece is obtained and then Eq. (17) must be solved with mentioned appropriate boundary conditions so that temperature distribution is obtained. The employed finite-element formulation for this equation is as follows:

4313

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

_ þ ½KfTg ¼ fQ g; ½CfTg

ð45Þ

where

½C ¼

½K ¼

ZZ

ZZ

þ

fQ g ¼

Z

ZZ

2pqC P ½NT ½Nr dr dz

2p



ð46Þ

T    T  ! Z Z @N @N @N @N þ ½k ½k r dr dz þ hc1 ½NT ½Nð2pRi Þdz þ hc2 ½NT ½Nð2pRO Þdz @r @r @z @z

hr ½NT ½Nð2pRO Þdz;

_ T r dr dz þ 2pg½N

Z

ð47Þ

hc1 T 1 ð2pRi Þdz þ

Z

hc2 T j ð2pRO Þdz þ

Z

hr T 1 ð2pRO Þdz:

ð48Þ

If the higher order terms are neglected, then

 Tn þ OðDtÞ; Dt nþ1 ¼ cT þ ð1  cÞT n ;

T T_ n ¼ T

nþc

nþ1

ð49Þ ð50Þ

where Tn+1 is the temperature distribution at n + 1 time step and Tn is temperature distribution at n time step. By varying the parameter c, different transient schemes can be constructed. In this problem c is set to 1, and the fully implicit scheme obtained. Finally, this equation can be rearranged as follows:

ð½C þ cDt½KÞfTgnþ1 ¼ ½CfTgn þ Dt cfQgnþ1 :

ð51Þ

After obtaining temperature distribution for each time step, Since material properties are temperature-dependent, the material properties including r, l, Cp and k must be modified. Thus it may some boundary conditions depend on temperature distribution too, and they must be refined. Therefore, the solution of these two coupled magneto-thermal equations with applying appropriate boundary conditions must be repeated in each time step until temperature distribution converges. The convergence criterion is (jTnew  Toldj)/jTnewj 6 109 in each time step. 7. Validation of numerical code (a) Electromagnetic solution: The numerical precision of the electromagnetic solution is validated by attaining the value of magnetic vector potential of two specified points that locate at the distance of 8 cm and 38 cm from the middle point of axis of the cylinder with the following properties and compared to those reported by Louai et al. [18], as shown in Table 1. In addition, these amounts with different grid numbers are shown to evaluate which grid number is more accurate. The results agree well with grid number of (50  25)

J s ¼ 5  106 A=m2 ;

f ¼ 50 Hz;

r ¼ 5  106 S=m; lr ¼ 90 H=m; L ¼ 60 cm; R ¼ RO ¼ 3 cm;

(b) Heating stage: In Fig. 5, a comparison with analytical solutions using Green’s function methods [8], is done to validate the numerical code. This analytical solution is a reasonable model for the heating stage of the induction heat treatment process of a solid cylinder. Material properties in this solution do not vary with temperature. The obtained results of numerical solution for the following properties (AISI 4140 steel) and for a solid cylinder (it means for a special case in which Ri = 0) are compared with analytical solution:

J s ¼ 1:5  1010 A=m2 ;

f ¼ 50 Hz;

r ¼ 3:5  106 S=m; lr ¼ 90 H=m; q ¼ 7850 kg=m3 ;

C p ¼ 500 J=kg K; k ¼ 41 W=m K; R ¼ RO ¼ 2 cm;

v ¼ 2:5 mm=s;

Z 0 ¼ 0:

In this case, the motion of coil begins from the initial position of Z0 = 0 (beginning of the cylinder). As shown in Fig. 5 they are in good agreement.

Table 1 Comparison of numerically obtained data with previous works.

38 cm from the middle of the cylinder 8 cm from the middle of the cylinder

Grid (30  15)

Grid (40  20)

Grid (50  25)

Louai et al. [18]

0.271  105 5.602  105

0.277  105 5.649  105

0.285  105 5.709  105

0.286  105 5.713  105

4314

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

Finite-element Solution Analytical Solution

1000

Temperature(K)

900

800

700

600

500

400 1

2

3

4

5

Time(s) Fig. 5. The comparison with analytical solutions using Green’s function methods for validation of results of heating stage.

1400 Finite-element Solution Experimental Solution by Thomas et al.

Temperature(K)

1200

1000

800

600

400

0

20

40

60

80

Time(s) Fig. 6. The comparison with experimental results for validation of obtained results of relevant to cooling stage.

(c) Cooling stage: In Fig. 6 temperature of a point which locates on the central axis of cylinder with the same distance from two sides of it, compared with that presented experimentally by Thomas et al. [11], to validate the obtained result related to spray cooling of cylinder. In these experiments a AISI 4140 steel solid cylinder is cooled from initial temperature of 1273 K using the mentioned water–air spray with the following characteristics:

R ¼ RO ¼ 2 cm; L ¼ 10 cm;

v ¼ 11:5:

In this case, cooling spray does not have any movement. As observed in Fig. 6, they are in good agreement.

8. Results and discussion In order to have a physical point of view of the problem and for the purpose of calculating the temperature distribution in a vertical hollow circular cylinder, verified numerical calculations are carried out. The finite-element procedure was applied

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

4315

to simulate the complete heat treatment cycle including both heating stage and spray quenching in an AISI 4140 steel hollow cylinder with the following geometry specifications according to Fig. 2.

Ri ¼ 12 mm; RO ¼ 20 mm; R3 ¼ 22 mm; R4 ¼ 26 mm; R1 ¼ 40 cm; L ¼ 40 cm; a11 ¼ a21 ¼ 40 cm: The current density and other field current applied to 4  4 coil of finite cross-sectional area hold:

J S ¼ 1:85  1010 A=m2 ;

f ¼ 50 Hz:

Since material properties of work-piece are temperature-dependent, the material properties [18] of AISI 4140 steel are considered. The relative permeability is defined by:

lr ¼

l ; l0

ð52Þ

where l0 is the permeability of free space or magnetic constant given by l0 = 4p  107 H/m. Correction of the relative permeability with respect to the temperature may be performed by means of the following relation, [10]:

lr ðTÞ ¼ 1 þ ðlr;20  1Þ  uðTÞ:

ð53Þ

The correction function u(T) is displayed in Fig. 7 and lr,20 = 90 [19]. The relative permeability of AISI 4140 steel was taken from Eq. (53) below the Curie temperature (for the specified material is about 770 °C, [19]) and above the curie temperature lr = 1. Fig. 8 shows the cylinder and the coil and their relative initial positions. The points of A, C and E locates on the external surface of cylinder and the points of B, D and F locates on the internal surface of the cylinder. Thus, in the following figures the points of A, C and E are representative of the points of external surface and the points of B, D and F are representative of the points of internal surface. The coil and spray move along the cylinder at the same velocity of v = 2 mm/s. The coil starts its movement from the initial position of Z0 = 5 cm. The width of the spray cooling band is 10 cm as said before, with the upper boundary of the band located at a distance of 4 cm below the coil center. The ratio of the water-mass-flow rate to the airmass-flow rate for the cooling spray is v = 12.5. The cylinder is assumed to be at initial temperature 25 °C that is the ambient temperature. The jet stream temperature can be obtained by the following equation [20]:

T j ¼ T 1 þ 5  C:

ð54Þ

Fig. 9 shows the calculated temperature–time histories at a number of points of cylinder which are specified in Fig. 8. The temperature variations with time at points A and B represent transient response, whereas the results at those points farther away from the initial position of coil reach the steady-state solution. From this fig it is clear that all points of cylinder approach higher than austenitizing temperature which for the specified material is about 800 °C during heating process, and this means power of heat source is intense enough. In heating stage the temperatures of the points on the outer surface of cylinder are greater than temperatures of the points of inner surface at the same axial position, because the intensity of heat source produced by induction of the coil in the hollow cylinder increases with growth of radius. During heating,

1 0.9 0.8 0.7

φ

0.6 0.5 0.4 0.3 0.2 0.1 0 300

500

700

Temperature,K

900

1100

Fig. 7. The variation of correction function u(T) related to AISI 4140 steel with temperature.

4316

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

Fig. 8. The exhibition of the cylinder and the coil and their relative initial positions.

1200 Point "A" Point "B" Point "C" Point "D" Point "E" Point "F"

Temperature(K)

1000

800

600

400 50

100

150

200

250

300

Time(s) Fig. 9. The temperature–time histories at points of A, B, C, D, E and F of the hollow cylinder.

the temperature of each point rises till the moving coil reaches the axial location of that point. After passing of the moving coil temperature of the specified point goes up to some extent with less slope. After the specified point reaches its maximum temperature, the temperature of this point begins to fall. This reduction in temperature occurs for two reasons: First, the power produced by inductor varies exponentially from the position of coil to two ends of cylinder, and for this problem just 15% of total amount of power exists out of the vicinity of 4 cm of coil. Second, the amount of heat flux relevant to two-phase heat transfer of spray is great enough relative to other mechanisms of cooling. This happens when the upper band of the moving spray approaches that point for the points on the outer surface, and for the other points with some delay depending

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

25s 125s 190s 225s 250s

1400

1200

Temperature(K)

4317

1000

800

600

400 0

0.1

0.2

0.3

0.4

z(m) Fig. 10. The distribution of the temperature along the external surface of cylinder for the five different times of 25, 125, 190, 225 and 250 s.

25s 125s 190s 225s 250s

1400

Temperature(K)

1200

1000

800

600

400 0

0.1

0.2

0.3

0.4

z(m) Fig. 11. The distribution of the temperature along the internal surface of cylinder for the five different times of 25, 125, 190, 225 and 250 s.

on distance from outer surface. Therefore, the maximum in temperature–time curve relevant to the point placed on inner surface of the cylinder e.g. point B occurs later a bit than this for the point placed on outer surface of the cylinder e.g. point A. At the beginning of cooling stage, since the temperature is high, the heat transfer coefficient due to spray cooling is relevant to radiation-dominated region, and gradually with decreasing of temperature the heat transfer coefficient is evaluated by transient region and convection-dominated region equations. Therefore the heat transfer coefficient due to spray cooling increases with passing the time. Consequently, the slope of the temperature–time curve increases with time until the moving spray completely traverses the specified point. Then cooling of cylinder happens only by natural convection with air. And since the heat transfer coefficient of natural convection is so lower than the heat transfer coefficient due to spray, cooling of cylinder occurs at a much lower rate. This process is investigated until passing 300 s from the beginning of the process.

4318

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

940

time=25s (z=0) time=90s (z=0)

920

Temperature(K)

900 880 860 840 820 800 780 760 740 720

0

0.002

0.004

0.006

0.008

ξ Fig. 12. Radial variation of temperature at a specific axial position of z = 0 at the time of 25 and 90 s.

Fig. 13. The distribution of the temperature within the axial cut through the cylinder at the time of 25 s and the coil position of zc = 0.

Fig. 10 shows distribution of the temperature of points that locates along the external surface of cylinder (BC line in Fig. 2) at various times of process. The moving coil after about 225 s approaches the end of cylinder and stops, while the moving spray is still cooling the cylinder. Therefore, distribution of the temperature along the external surface of cylinder for the five different times of 25, 125, 190, 225 and 250 s are shown in Fig. 10. If in this fig. the line related to time of 125s is considered, we can see that in spite of the fact that at this time coil locates at the position of zc = 20 cm from the beginning of the cylinder, the maximum of curve occurred for a position of z = 16 cm. This fact is justified with increasing of temperature after passing the coil for each point of cylinder till the reaching the spray. In this fig. the reduction of exponentially of power of inductor from its position towards two ends of cylinder is also seen. For example, a point that locates in 5 cm of coil in each instance experiences the temperature reduction of 85% related maximum temperature. Also distribution of the temperature of points that locates along the internal surface of cylinder (AD line in Fig. 2) at various times of 25, 125, 190, 225 and 250 s is depicted in Fig. 11. The curves show the similar behavior to the curves of previous figure; However, the maximum in the curves occurs with delay related to previous fig. and the amount of maximum temperature is lower than maximum amount of previous figure by 1%. Radial variation of temperature at a specific axial position of z = 0 (AB line in Fig. 2) at the time of 25 and 90 s, is depicted in Fig. 12. At the time of 25 s the moving coil just passes the axial position of z = 0 and therefore this curve is related to

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

4319

Fig. 14. The distribution of the temperature within the axial cut through the cylinder at the time of 125 s and the coil position of zc = 20 cm.

heating stage. But at the time of 90 s this region is cooled with the moving spray; and thus this curve is related to cooling stage. It is clear from Fig. 12 that the temperature of the internal surface of cylinder does not differ very much from the temperatures of the external surface of cylinder relative to longitudinal variation of temperature, i.e. reduction of temperature of internal point from external point is 23% for heating stage and 19% for cooling stage. This fact occurs because of relatively high thermal conductivity of steel, and also slight thickness of hollow cylinder. The natural convection occurs on inner surface of cylinder with the ambient air, but the effect of this is low on cooling process. Also, from fig. and percentages we can see that the temperature difference between internal surface and external surface of cylinder for heating stage is greater a bit related to cooling stage. Fig. 13 shows distribution of the temperature within the axial cut through the cylinder at the time of 25 s while the moving coil reaches at the distance of zc = 0 from beginning of cylinder. Also Figs. 14 and 15 show distribution of the temperature at the time of 125 and 225 s after beginning of the process. At these times the moving coil reaches at the distance of 20 cm and 40 cm from the beginning of cylinder, respectively. All these three figures indicate the ability of concentrated heating of this method of heating which causes a drastic longitudinal temperature difference along a small distance. Finally, Fig. 16 shows temperature distribution at the final time of investigation of the process, namely 300 s. At this time heater and cooling

Fig. 15. The distribution of the temperature within the axial cut through the cylinder at the time of 225 s and the coil position of zc = 40 cm.

4320

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

Fig. 16. The distribution of the temperature within the axial cut through the cylinder at the time of 300 s.

spray have passed the ending of cylinder while ago, and they are not applied to the cylinder any more. At this time the cylinder cooled only by natural convection with the air. The difference between the maximum and minimum of temperature in the cylinder at this time is relatively low, i.e. 26%.

8.1. Effect of coil velocity on temperature distribution The effect of coil and spray velocity on the temperature distribution during induction heating process was analyzed. Fig. 17 shows the temperature–time histories at the points of A and C which are specified in Fig. 8, for three different velocities of coil, i.e.v = 1.5 mm/s, v = 2 mm/s and v = 2.5 mm/s. This investigation is performed for the coil distance of 5 cm from the beginning of cylinder, for all three conditions. It can be seen from Fig. 17 that decrease of velocity of the coil leads to higher obtained temperature for all points of the cylinder during the heating stage, and also lower temperature during the cooling stage. For example, the cylinder with the coil velocity of v = 1.5 mm/s experiences the maximum temperature which is 10% greater than maximum temperature of the cylinder with the coil velocity of v = 2 mm/s, and 16% greater than maximum temperature of the cylinder with the coil velocity of v = 2.5 mm/s. Because the less velocity the coil has, the more

1400

point "A" (v=1.5 mm/s) point "C" (v=1.5 mm/s) point "A" (v=2 mm/s) point "C" (v=2 mm/s) point "A" (v=2.5 mm/s) point "C" (v=2.5 mm/s)

Temperature(K)

1200

1000

800

600

400 50

100

150

200

250

300

Time(s) Fig. 17. The temperature–time histories at the points of A and C for different velocities of coil of

v = 1.5 mm/s, v = 2 mm/s and v = 2.5 mm/s.

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

4321

Zc=0 , z=0

1100

v=1.5 mm/s v=2 mm/s v=2.5 mm/s

1050

Temperature(K)

1000 950 900 850 800 750 700

0

0.002

0.004

0.006

0.008

ξ Fig. 18. Radial variation of temperature at axial position of z = 0 and for the position of the coil of zc = 0 for three different velocities of coil.

time is spent for both heating and cooling stages. In the other words, with decreasing of the velocity, the transferred power of moving heat source to each point, and also the time of exposure of each point to the cooling spray both increase. For example, for the three different velocities of coil, i.e. v = 1.5 mm/s, v = 2 mm/it s and v = 2.5 mm/s, time of exposure to spray are 66.7, 50 and 40 s respectively. Also the difference between the temperature of inner surface and outer surface of the hollow cylinder increases with decreasing of velocity of the coil and spray. This fact is observed in Fig. 18 that depicts radial variation of temperature at a specific axial position of z = 0 and for the position of the coil of zc = 0 for three different mentioned velocities. To sum up, it can be said that decrease in velocity leads to greater temperature fluctuations; and, this variations cause the cylinder to experience more thermal stress which increases possibility of the occurrence of quench cracks. 8.2. Effect of geometry on temperature distribution For investigating the effect of geometry on the temperature distribution, effect of inner radius to outer radius ratio of the cylinder on temperature distribution is considered. For this reason Fig. 19 shows the temperature–time histories at the

Point "C" (Ri/Ro=0.6) Point "D" (Ri/Ro=0.6) Point "C" (Ri/Ro=0.4) Point "D" (Ri/Ro=0.4) Point "C" (Ri/Ro=0.8) Point "D" (Ri/Ro=0.8)

1200

Temperature(K)

1000

800

600

400 20

40

60

80

100

120

140

160

Time(s) Fig. 19. The temperature–time histories at the points of C and D, for three different inner to outer radius ratios of the cylinder of Ri/RO = 0.4, Ri/RO = 0.6 and Ri/RO = 0.8.

4322

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

t=25s , z=0

Ri/Ro=0.4 Ri/Ro=0.6 Ri/Ro=0.8

Temperature(K)

900

800

700

600

0

0.25

0.5

0.75

1

ψ Fig. 20. The normalized radial variation of temperature at axial position of z = 0 at the time of 25 s for three different radius ratios.

points of C and D, for three different inner to outer radius ratios of the cylinder, i.e.Ri /RO = 0.4, Ri/RO = 0.6 and Ri/RO = 0.8. Increase in radius ratio of cylinder or in the other words decrease in thickness of hollow cylinder results in decrease in the difference between the temperature of inner surface and outer surface of the hollow cylinder so that in the case of Ri/RO = 0.8 the temperature of inner surface and outer surface of the cylinder approximately coincide. This fact is obviously confirmed in Fig. 20 that depicts normalized radial variation of temperature at a specific axial position of z = 0 at the time of 25 s for three different mentioned radius ratios. By defining the normalized parameter of w ¼ RO nR , radial variation of i temperature will be dependent of thickness of wall of the cylinder. This fig. indicates that the cylinder with the radius ratio of Ri/RO = 0.4 experiences the temperature difference between its outer and inner surface which is 40% greater than temperature difference of the cylinder with the radius ratio of Ri/RO = 0.6, and 82% greater than temperature difference of the cylinder with the radius ratio of Ri/RO = 0.8. It is also observed from Fig. 19 that by decreasing the radius ratio of cylinder, the difference between curves of related to a point on external surface of cylinder and its corresponding point on a internal surface of cylinder, the difference between the amount of maximum of these two curves, and the delay in occurrence of maximum of these two curves increase. Generally, increase in thickness of hollow cylinder, rises the temperature fluctuation of the cylinder; and, as a result it is threatened more by quenching crack.

1200

Point "A" (Z0= -5 cm) Point "C" (Z0= -5 cm) Point "A" (Z0= -2.5 cm) Point "C" (Z0= -2.5 cm) Point "A" (Z0=0) Point "C" (Z0=0)

Temperature(K)

1000

800

600

400 50

100

150

200

250

300

Time(s) Fig. 21. The temperature–time histories at the points of A and C for different initial position of the coil of z0 = 5 cm, z0 = 2.5 cm and z0 = 0.

H. Shokouhmand, S. Ghaffari / Applied Mathematical Modelling 36 (2012) 4304–4323

4323

8.3. Effect of initial position of coil on temperature distribution Finally, the effect of initial position of the coil on the temperature distribution is investigated. Fig. 21 shows the temperature–time histories at the points of A and C, for three different initial position of the coil, i.e.z0 = 5 cm, z0 = 2.5 cm and z0 = 0. The farther the initial position of the coil is from the beginning of the cylinder, the more steady-state solution is obtained for the beginning points of the cylinder. In the other words, by increasing the distance of initial position of coil from the beginning of the cylinder, the temperature–time curve related to beginning point of the cylinder becomes closer to temperature–time curves of other points of the cylinder. Consequently, it is more appropriate to increase the distance between the initial position of the coil and the beginning of the cylinder. However, there is a restriction in increasing this distant because of economical reasons. Since 95% of power of the inductor exists in the vicinity of 5 cm of the coil; by initial position of coil of greater than 5 cm from the beginning of the cylinder, significant energy loss is created. Then altogether, z0 = 5 cm seems to be a proper choice, because it causes less temperature fluctuations and is safer against cracks, also it causes less energy loss. 9. Conclusion The thermal investigation of complete moving induction heat treatment cycle of steel for vertical hollow circular cylinder including heating stage and water–air spray cooling has been studied. Solution to this problem described by the solution of coupled Maxwell’s equations (for magnetic analysis) and Fourier’s equation (for thermal analysis) with temperature-dependent properties and time-variable and moving boundary conditions using finite-element methods. The effects of natural convection with air and radiation on the both internal and external surfaces of cylinder have been taken into account. Available experimental relations for spray cooling process of steel are used for obtaining forced convection heat transfer coefficient for outer surface due to spray cooling. The simulation procedure can be used as a useful tool in induction coil design and in the selection of process parameters. The temperature distribution within work-piece at different times and temperature–time curve for different points are presented and discussed in this paper. This temperature fields predicted by the model are used both to predict the resulting microstructure and as an input for the thermal-stress model to predict the occurrence of quench cracks. It is also observed that decrease of velocity of the coil and spray leads to greater temperature fluctuations during the process. Increase in thickness of hollow cylinder and decrease in the distance of initial position of coil, rises the temperature fluctuations, too. Therefore, from the design point of view, high velocity of coil is more proper because of reduction of possibility of quench cracks. Similarly, the initial position of z0 = 5 it cm is a optimum choice, because of its greater safety against cracks and its less energy loss. In addition, hollow cylinders with greater thickness are exposed to occurrence cracks more than those with smaller thickness. References [1] E.J. Davies, Conduction and Induction Heating, Peter Peregrinus, United Kingdom, London, 1990. [2] F. Bay, V. Labbe, Y. Favennec, J.L. Chenot, A numerical model for induction heating processes coupling electromagnetism and thermomechanics, Int. J. Numer. Methods Eng. 58 (2003) 839–867. [3] B. Drobenko, O. Hachkevych, T. Kournyts’kyi, A mathematical simulation of high temperature induction heating of electro conductive solids, Int. J. Heat Mass Transfer 50 (2007) 616–624. [4] M.S. Huang, Y.-L. Huang, Effect of multi-layered induction coils on efficiency and uniformity of surface heating, Int. J. Heat Mass Transfer 53 (2010) 2414–2423. [5] K.F. Wang, S. Chandrasekar, Henry T.Y. Yang, Finite-element simulation of induction heat treatment, J. Mater. Eng. Perform. 1 (1992) 97–112. [6] K. Sadeghipour, J.A. Dopkin, K. Li, A computer aided finite element/experimental analysis of induction heating process of steel, Int. J. Comput. Indus. 28 (1996) 195–205. [7] Matej Kranjc, Anze Zupanic, Damijan Miklavcic, Tomaz Jarm, Numerical analysis and thermographic investigation of induction heating, Int. J. Heat Mass Transfer 53 (2010) 3585–3591. [8] K.F. Wang, S. Chandrasekar, Henry T.Y. Yang, Finite-element simulation of moving induction heat treatment, J. Mater. Eng. Perform. 4 (1995) 460–473. [9] J.O. Yun, Young-Soo Ynag, Analysis of the induction heating for moving induction coil, J. Mech. Sci. Technol. 20 (2006) 1217–1223. [10] Ivo Dolezel, Jerzy Barglik, Bohus Ulrych, Continual induction hardening of axi-symmetric bodies, J. Mater. Process. Technol. 161 (2005) 269–275. [11] Thomas, M. Ganesa-Pillai, P.B. Aswatch, K.L. Lawrence, A. Haji-Sheikh, Analytical/finite-element modeling and experimental, Metall. Mater. Trans. A 29A (1998) 1485–1498. [12] Frank P. Incropera, David P. DE Witt, Introduction to Heat Transfer, fourth ed., John Wiley & Sons Ltd., New York, 2002. [13] M. Al-Arabi, M. Khamis, M. Abd-ul-Aziz, Heat transfer by natural convection from the inside surface of a uniformly heated vertical tube, Int. J. Heat Mass Transfer 34 (1991) 1019–1025. [14] Wei Qi, Yang Shiming, Experimental study of natural convection heat transfer of air layers in vertical annuli under high rayleigh number conditions, J. Heat Transfer 28 (1999) 50–57. [15] P. Buckingham, A. Haji-Sheikh, Cooling of high-temperature cylindrical surfaces using a water–air spray, J. Heat Transfer Trans. ASME 117 (1995) 1017–1028. [16] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, John Wiley & Sons Ltd., New York, 1960. [17] R.W. Lewis, P. Nithiarasu, K.N. Seetharamu, Fundamental of the Finite Element Method for Heat and Fluid Flow, John Wiley & Sons Ltd., New York, 2004. [18] F.Z. Louai, N.N. Said, S. Drid, Numerical analysis of electromagnetic axisymmetric problems using element free Galerkin method, J. Electr. Eng. 57 (2006) 99–104. [19] M.F. Rothman, High-Temperature Property Data: Ferrous Alloys, ASM International, Metal Park, Ohio, 1988. [20] G.E. Totten, C.E. Bates, N.A. Clinton, Hand Book of Quenchants and Quenching Technology, ASM International, 1993 (first printing).