A Finite Element Simulation for Induction Heat Treatment of Automotive Drive Shaft

Abstract

A finite element model was developed to predict deformation, temperature, phase fraction and hardness during heat treatment of an automotive drive shaft. The heat generation due to induction was treated as one of the boundary conditions for heat flux on the specimen together with the conduction heat loss during quenching. As for diffusional transformation, the transformation kinetics were modeled by Johnson–Mehl–Avrami–Kolmogorov equation, whereas the Marburger equation was used for displacive martensitic transformation. The transformation plasticity was considered through the constitutive equations corresponding to each transformation mechanism and these equations were incorporated into the finite element model. Besides the transformation plasticity, an implicit procedure to calculate the thermo-elasto-plastic deformation was implemented in the model. The prediction accuracy for phase evolution, residual stress, hardness and dimensional change of the specimen was verified from the measured data. The effect of transformation plasticity on whole deformation behavior was described by the developed model.

1. Introduction

Induction heat treatment is a convenient and precise method for heating metals to a set temperature and is widely used in current industrial processes including welding, forging, surface hardening and other applications. Induction heating has a number of benefits over traditional methods, among them selective hardening without affecting the rest of the material, high efficiency, and no safety issues.¹⁾ Two types of induction heat treatments are generally used: stationary induction heat treatment and moving induction heat treatment.²⁾ In stationary induction heating, no relative motion occurs between the induction coil and the workpiece. The current-bearing coils produce an electromagnetic field which induces eddy currents in the workpiece and the heat is generated resistively. When the workpiece reaches a set temperature, quenching begins. In moving induction heat treatment, relative motion occurs between the moving coil and the workpiece. Quenching is done by a quenchant which is injected from the backside of the coil. Induction heating by a moving coil is mainly used for surface hardening of elongated parts such as driveshaft or piston pin.

Since induction heating is a complex process involving phase transformation which shows nonlinear behavior because of the different atomic volumes of the austenite and martensite phases,³⁾ and since heat transfer is related to the electromagnetic field, the process must be modeled to investigate the temperature history and subsequent phase distribution. Over the past decades, numerous studies have been carried out to analyze induction heat treatment, in particular several on stationary induction heating^4,5) with thermal-elastic-plastic analysis. As high-frequency induction heat treatment generally involves phase transformation, investigation of phase transformation kinetics is important for analyzing phase distribution and resulting residual stress. Yuan et al.⁶⁾ developed an integrated modeling system based on finite element (FE) analysis that considers the electromagnetic-field-induced heating and phase transformation during holding in the austenite state and severe quenching processes. Studies on moving induction heat treatment have also been made recently, but most of them focus only on heat transfer analysis without considering phase transformation. Shokouhmand et al.⁷⁾ analyzed moving induction heat treatment and the subsequent cooling process by moving a water-air spray for a hollow cylinder. Bai et al.⁸⁾ employed electromagneto-thermomechanical FE method to simulate the thermal and mechanical behavior during induction-assisted weld-based additive manufacturing. The validity of the model proposed there was confirmed by comparing the simulated and experimentally measured residual stress. Wang et al.²⁾ developed a comprehensive model of induction heat treatment of a steel cylinder with moving coil incorporating the mathematical modeling of the phase transformations as well as the heat transfer occurring during heating and quenching.

The present effort integrates a FE modeling and the capability to simulate the temperature history, phase distribution, and hardness where the driveshaft has a complicated shape. The hardness and residual stress were experimentally measured and compared with the calculated results to verify the proposed model. The goal of the modeling system is to guide optimal design and process conditions of induction systems.

2. Experiments

2.1. Material

The material employed in this study is a boron-added special commercial steel (SAE10B38). Boron generally improves the hardenability of steel since it is segregated in the austenite grain boundary, reducing the grain boundary energy. In addition to this, boron is relatively cheap compared with the other alloying elements such as Mo, Cr, Ni, etc., and the effect on hardenability is remarkable even with contents as low as tens of ppm. The chemical composition of investigated steel is shown in Table 1.

Table 1. Chemical composition (wt%) of investigated steel.

Specimen	C	Si	Mn	P	S	B	Cu	Ni	Mo
SAE10B38	0.41	0.25	0.76	<0.03	0.025	<0.005	<0.2	<0.2	<0.05

During the induction heating process, the surface of the shaft is heated rapidly and then quenched to produce a compressive residual stress. This residual stress comes from a martensitic transformation since the lattice parameter along the c-direction of body-centered tetragonal (BCT) martensite is larger than in BCC ferrite. Thermo-Calc software (with database of TCFE2) was used to calculate the thermodynamic phase diagram under para-equilibrium (Fig. 1) and for a steel with carbon content of 0.41 wt%, the austenitizing temperature is approximately 761°C. The martensitic transformation temperature was obtained from the empirical equation⁹⁾ which is expressed as a function of chemical composition (wt.%)   

$$\begin{array}{l}{M}_S(°C)=525-350(C-0.005)-45Mn-30Cr\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-20Ni-5Si-35V(Nb+Zr+Ti)\end{array}$$(1)

Fig. 1.

Equilibrium phase diagram of investigated steel.

The microstructures of the inner and outer shaft regions are shown on Fig. 2. The shaft surface consists of a martensite phase whereas the inner region consists of ferrite and pearlite phases. The physical properties used in the calculation process were based on the average values of the ferrite phase and the pearlite phase, taking into account the respective phase fractions.

Fig. 2.

Microstructure of investigated steel. (a) hardened martensite layer (surface layer) and (b) ferrite and pearlite layer (inner layer). (Online version in color.)

2.2. Induction Heating

Exact temperature measurement upon induction heating is essential because temperature is crucial for phase transformation of steel. The temperature was measured by IR (infrared) camera. Since the measured data varies considerably with material emissivity, it was estimated indirectly using dilatometer. The emissivity was calibrated by comparing the temperature measured by a K-type thermocouple. Comparison of these data yielded an emissivity of 0.4 and the resulting maximum temperature was calculated as approximately 960°C. Accurate temperature measurement could not be performed after the peak temperature because quenchant was ejected from the rear side of the induction coil during the actual process.

2.3. Indentation

The shafts were cut into a plane perpendicular to the axial direction in order to evaluate the residual stress distribution on the inner surface using indentation. The indentation test was performed ten times with the indenting axis perpendicular to the plane. The Micro-AIS system (Frontics, Inc., Seoul, Republic of Korea), with a displacement resolution of 10 nm and a load resolution of 0.1 mN, was used for this experiment. Indentation tests had a loading rate of 0.03 mm/min and a maximum indentation load of 490 mN. Evaluating the residual stress using indentation requires knowledge of the stress-free state with the same microstructure as the stressed state. The shaft was cut in half along the axial direction in order to obtain a stress-free state, and indentation testing was performed at a region close to the cut surface with reference to previous research.¹⁰⁾ The residual stress was evaluated using the following equation:¹¹⁾   

$${\sigma }_{res}^{avg}=\frac{3}{2}\frac{(L_0-L_s)}{A_s}$$(2)

where L₀ and L_s are the indentation load at maximum indentation depth in the stress-free and stressed state and A_s is the area of the residual imprint after the indentation test in the stressed state.¹²⁾ This means average value of in-plane principal stresses.

3. Model Development

3.1. Induction Heat Treatment Model

Induction heat treatment is a complex multiphysics problem requiring an electromagnetic solution. However, the high computational cost of solving an electro-thermo-mechanical problems incorporating phase transformation limits its use in practical area. Therefore, in this study, instead of solving the Maxwell equation, the heat generation was reproduced considering the relative magnetic field intensity along the depth of the shaft. The magnetic field intensity can be calculated by the diffusion equation¹³⁾   

$${\nabla }^{2}H=\frac{\mu }{\rho }\frac{\partial H}{\partial t}=\frac{{\partial }^{2}H}{\partial x^2}+\frac{{\partial }^{2}H}{\partial y^2}+\frac{{\partial }^{2}H}{\partial z^2}=\frac{\mu }{\rho }\frac{\partial }{\partial t}(H_x+H_y+H_z)$$(3)

The magnetic field is uniform along the circumferential and tangential direction of shaft since the coil is a circular ring. Thus, the equation above can be simplified to   

$$H_x(y)=H_{0}exp(-\alpha y)$$(4)

and α can be expressed by skin depth δ:   

$$\alpha =\frac{1}{\delta }=\sqrt{\frac{\mu \omega }{2\rho }}$$(5)

where μ, ω, and ρ is the permeability (4π×10⁻⁷), frequency, and electrical resistivity respectively. The frequency used in the induction heating process was 4000 Hz. The electrical resistivity is a function of a temperature, and is listed in Table 2 over the temperature ranges.¹⁴⁾

Table 2. Electrical resistivity of investigated steel at various temperatures.

T (°C)	Electrical resistivity (10−6 Ω∙m)
20	0.1675
200	0.2947
400	0.4925
600	0.7627
800	1.1149
1000	1.5346

3.2. Phase Transformation Model

The phase transformation between ferrite/pearlite and austenite is classified as a diffusional precipitation reaction. The transformed fraction, X, is   

$$X=1-exp(-k{t}^n)$$(6)

which is known as a Johnson-Mehl-Avrami-Kolmogorov (JMAK) equation. Here, t is the total time for transformation at a given temperature; k is temperature-dependent constants, whereas the time exponent n is independent of temperature, provided there is no change in the transformation mechanism. In the case of induction heating, the phase transformation takes place under non-isothermal conditions. Thus, the additivity rule was adopted to take into account the variation of kinetics values under rapid cooling. In the additivity rule, time is divided into small intervals so that the kinetics variables during cooling remain constant. From this, the transformed phase fraction of the ith time step is:   

$$X_i=1-exp(-X_{i-1}^{ex}+\mathrm{\Delta }{X}_i^{ex}),\mathrm{\hspace{0.17em} } \mathrm{\Delta }{X}_i^{ex}=n{k}_{i}t{\text{'}}_i^{n-1}\mathrm{\Delta }t,\mathrm{\hspace{0.17em} } t{\text{'}}_i={\left(\frac{X_{i-1}^{ex}}{k_i}\right)}^{\frac{1}{n}}$$(7)

where X^ex is an extended volume fraction under the condition that nucleation and geometrical impingement during growth are not considered. t’ is the equivalent time needed to reach a transformed fraction of $X_{i-1}^{ex}$ at the temperature of the ith step and Δt is the time interval corresponding to the ith step. The constants k and n are determined from the dilatation data for continuous cooling, applying an inverse additivity technique,¹⁵⁾ and the resulting parameters are listed in Table 3. The constants C_A, $A_{e_3}$, and $A_{e_1}$ are the weight percent carbon of undeformed austenite and the equilibrium transformation temperatures respectively, which are calculated with Thermo-Calc under para-equilibrium. The B_s, AGS (μm), and T are the bainite start temperature based on Zhao’s empirical equation,¹⁶⁾ prior austenite grain size, and temperature (°C) respectively. The diffusion coefficient of carbon within the austenite matrix is described as a function of temperature and carbon content:¹⁷⁾   

$$\begin{array}{l}{D}_c=4.53\\ \hspace{1em}\hspace{0.5em}\times {10}^{-7}\left\{1+C_A(100-C_A)\frac{0.19}{T}\right\}exp\left\{-\left(\frac{1}{T}-2.22\times {10}^{-4}\right)\left(17\mathrm{\hspace{0.17em} }776-57.67C_A\right)\right\}\end{array}$$(8)

Table 3. Values of k and n from JMAK equation for various phase transformations.

	ln(k)	n
Austenite to ferrite	$\begin{array}{l}39.99+4.87\cdot ln(AGS)+(-3.89-39.98\cdot C_A)\cdot ln(D_C)\\ +\left(-3.93+39.56\cdot C_A^{-8.19}-5.95\cdot C_A^2\right)\cdot ln\left(\frac{A_{e_3}-T}{A_{e_3}}\right)\\ +39.97\cdot ln\left[sin\left(\frac{2\pi \left(T+273\right)}{4\left(A_{e_3}-80\right)}\right)\right]\end{array}$	39.99
Austenite to pearlite	$\begin{array}{l}7.23-4.14\cdot ln(AGS)+(-5.08-5.39\cdot C_A)\cdot ln(D_C)\\ +\left(1.37+1.67\cdot C_A^{-5.73}+6.33\cdot C_A^2\right)\cdot ln\left(\frac{A_{e_1}-T}{A_{e_1}}\right)\\ -0.8\cdot ln\left[sin\left(\frac{2\pi \left(T+273\right)}{4\left(A_{e_1}-80\right)}\right)\right]\end{array}$	1.24
Austenite to bainite	$\begin{array}{l}6.96-1.68\cdot ln(AGS)+(0.1-3.39\cdot C_A)\cdot ln(D_C)\\ +\left(-5.59+6.78\cdot C_A^{-6.08}+6.41\cdot C_A^2\right)\cdot ln\left(\frac{A_{B_s}-T}{A_{B_s}}\right)\\ +2.23\cdot ln\left[sin\left(\frac{2\pi \left(T+273\right)}{4\left(A_{B_s}-80\right)}\right)\right]\end{array}$	1.92

The diffusionless martensitic transformation, since it is dependent only on the temperature, can be described from Koistinen and Marburger’s empirical equation as   

$$\frac{X_M}{1-X_F-X_P-X_B}=1-exp\left[-0.011\left(M_S-T\right)\right]$$(9)

where the subscripts M, F, P, and B mean martensite, ferrite, pearlite, and bainite respectively. Phase transformation during heating is a single phase transformation to austenite. Since the initial phase of the investigated steel was pearlite, the phase transformation has the simple form of a single-phase to single-phase transformation. Thus the basic form of the JMAK equation was adopted for heating condition:   

$$X=1-exp\left(-A\cdot {10}^5\cdot exp\left(\frac{-B}{RT}\right)\cdot t^n\right)$$(10)

The constants were optimized from dilatometric analysis and the resulting parameters are listed in Table 4. The phase transformation behavior of investigated steel is shown in Fig. 3 under various cooling rates.

Table 4. Heating kinetics constants during pearlite-austenite phase transformation (V=ΔTemp/ΔTime).

logV = p + q∙logA		B	n
p	q
−10.1352	2.71168	1.0∙105	3.0

Fig. 3.

Comparison of measured (symbols) and calculated (lines) transformation kinetics of investigated steel under various cooling rates. (Online version in color.)

3.3. Transformation Plasticity Model

In this study, the transformation plasticity model proposed by Han et al.^18,19) was incorporated to take into account the permanent deformation during phase transformation. According to this theory, when stress is applied during the phase transformation, cavities are formed by the stress component in the direction perpendicular to the grain boundaries. From this, Coble creep is accelerated by the difference in vacancy concentration between the grain boundaries perpendicular and horizontal to the stress direction, and this results in permanent deformation even at stresses below the yield stress.^20,21) Using this concept, the plastic strain rate under diffusional transformation can be expressed as a function of transformation rate ($\dot{X}$), applied stress (σ), and temperature:   

$${\dot{\epsilon }}^{TP}=\frac{1}{3}\frac{d_0}{{\delta }_d}\dot{X}\frac{\sigma \mathrm{\Omega }}{k_{B}T}{c}_{v0}exp\left(-\frac{Q_f}{k_{B}T}\right)$$(11)

where d₀, δ_d, Ω, c_v0, Q_f, and k_B represent the initial grain size of the parent phase, the effective thickness of the interface, the vacancy volume, a dimensionless constant associated with vacancy formation, the formation enthalpy of the vacancy at the interface, and Boltzmann’s constant (1.38×10⁻²³ J/K), respectively. The initial grain size of prior austenite was assumed to be 30 μm and the effective thickness of the interface, the vacancy volume and the formation enthalpy of the vacancy at the interface were taken to be 1 nm, 1.21×10⁻²⁹ m³ and 80 kJ/mol, respectively. To obtain the constant c_v0, the additional dilatometric analyses were conducted to figure out the effect of external stress. The compressive stress states of 3, 5, and 7 MPa were applied while cooling with cooling rate of 10°C/s during the dilatometric experiments. The constants c_v0 could be optimized from the constrained Rosenbrock technique,²²⁾ which was performed by systematically changing the constant until the sum of the squared difference between the measured and calculated data reached a minimum value. The optimized value for investigated steel was 0.5.

In the diffusionless transformation, the anisotropic strain arises from the increased selection of some specific variants by the externally applied stress.²³⁾ The transformation strain increases linearly with the applied stress. Thus, to explain the transformation plasticity during the martensite transformation, a linear model was adopted:   

$${\dot{\epsilon }}^{TP}=A\dot{X}\sigma$$(12)

where the constant A was obtained experimentally as 2∙10⁻¹⁰Pa⁻¹.²⁴⁾

3.4. Heat Transfer Model

During the induction heating process, the temperature change is calculated considering the latent heat generated in the phase transformation process. The heat transfer equation can be expressed in the following form:   

$$\begin{array}{l}\rho C_p\dot{T}=\nabla \cdot \left(k_t\nabla T\right)+\mathrm{\Delta }{H}_F\cdot {\dot{X}}_F+\mathrm{\Delta }{H}_P\cdot {\dot{X}}_P\\ \hspace{1em}\hspace{1em}\hspace{1em}+\mathrm{\Delta }{H}_B\cdot {\dot{X}}_B+\mathrm{\Delta }{H}_M\cdot {\dot{X}}_M+\mathrm{\Delta }{Q}_{ind}\end{array}$$(13)

where ρ and k_t are density and thermal conductivity. The ΔQ_ind is the heat added to reproduce the temperature increase during the induction heat treatment, which is detailed in section 4. The density of each phase is given as a function of temperature and composition of the investigated steel (Table 5).²⁵⁾ C_γ and C_α are the carbon content in prior austenite and the solute carbon content in ferrite, which are calculated from thermodynamic considerations. The density of pearlite is determined from the density of ferrite and cementite since it is simply a mixture of these two phases. The density of bainite is assumed as that of ferrite since the difference between them is negligible and the density of martensite is calculated by the same formula as ferrite, but is calculated by substituting the carbon content with the carbon content of the austenite phase. The temperature-dependent isotropic thermal conductivity is also listed in Table 6. The temperature-dependent heat capacity and change of enthalpy during the phase transformation are calculated from thermodynamic analysis (Fig. 4). Phase transformation from austenite to bainite generally involves shear deformation, and the free energy change during this shear deformation is known to be about 600 J/mol.²⁶⁾ To determine the heat of transformation of bainite, the additional shear energy was subtracted from the transformation heat of ferrite. The heat generated from the martensitic transformation was calculated:²⁴⁾   

$$\mathrm{\Delta }{H}_M=-2\mathrm{\hspace{0.17em} }100+\left(\frac{2\mathrm{\hspace{0.17em} }100}{T_0-M_s}\right)\cdot (T-M_s)$$(14)

Table 5. Phase densities as a function of temperature and chemical composition.

Phase	Density (kg/m3)
Austenite	8099.79−0.506∙T+(−118.26+0.00739∙T)∙Cγ−68.24Si−6.01Mn
Ferrite	7875.96−0.2970∙T−5.62∙10−5∙T2+(−206.35+0.00778∙T+1.472∙10−6∙T2)∙Cα−36.86Si−7.24Mn
Cementite	7686.45−0.0663T−0.000312T2

Table 6. Temperature-dependent elastic modulus, Poisson’s ratio, and thermal conductivity of investigated steel.

T (°C)	E (GPa)	Poisson’s ratio	Kt (W/m∙°C)
0	210	0.280	50
300	185	0.296	45
600	135	0.310	38
700	100	0.325	33
900	100	0.325	28
1000	100	0.325	28

Fig. 4.

(a) Temperature-dependent heat capacity of each phase and (b) enthalpy change during the phase transformation. (Online version in color.)

3.5. Constitutive Equation

The stress increment of the Cauchy stress, dσ, is   

$$d\sigma =C^e:d{\epsilon }^e$$(15)

where C^e and dε^e are the elastic stiffness tensor and elastic strain increment, respectively. The total strain increment, dε^T, is   

$$d{\epsilon }^T=d{\epsilon }^e+d{\epsilon }^v+d{\epsilon }^{TP}+d{\epsilon }^p$$(16)

where dε^v, dε^TP and dε^p are the volumetric strain increment resulting from the temperature change and phase transformation, the transformation plasticity strain increment due to the diffusional and diffusionless phase transformations,^23,24) and the plastic strain increment, respectively. For the volumetric strain, the linear mixture rule of strain increment of existing phases was adopted as follows:²⁷⁾   

$$d{\epsilon }^v=\sum X_{i}d{\epsilon }_i^v=\sum X_i\cdot \left(\frac{1}{3}\frac{{\rho }_i}{d{\rho }_i}\right)I$$(17)

where X_i is the phase fraction of each phase, ρ_i is the density of each phase (listed in Table 5), and I is the identity tensor. The density of martensite was calculated using the same formula as ferrite, assuming that it has the carbon content of austenite. The plastic strain increment was not taken into account in this study since the external stress developed during the induction heating process can be ignored. Thus, the stress increment becomes:   

$$d\sigma =C^e:(d{\epsilon }^T-d{\epsilon }^v-d{\epsilon }^{TP})$$(18)

The temperature dependent isotropic elastic modulus was used to obtain the elastic stiffness of carbon steel, as listed in Table 6,²⁸⁾ and the Poisson’s ratio of the steel is also listed in Table 6.

4. Finite Element Modeling

A thermomechanical calculation for the induction heating process was performed. Four-node axisymmetric thermally coupled quadrilateral continuum elements (CAX4T) were used for the shaft (Fig. 5(a)). Considering the symmetry of the shaft shape, only one side was analyzed and the induction coil moves from the leftmost part to the rightmost part of the shaft. To evaluate the heat generated within the part, it is necessary to solve Maxwell’s electro-magnetic field equations. In this study, however, the internal heat generation in the material was reproduced by simply lowering the internal energy for practical use. As mentioned in section 3.1, the relative magnetic field intensity is a function of depth from the shaft surface. The distribution was determined from the skin depth derived from Eq. (5), and the relative magnetic field density along the depth from the shaft surface is shown in Fig. 6. As the power is proportional to the square of the magnetic field intensity, the relative heat generated throughout the shaft can be calculated from the distribution of magnetic field intensity. This approach is schematically shown in Fig. 7. As shown on Fig. 7(a), the whole coordinates of surface nodes are extracted and the shortest distance from the extracted surface coordinates to each integration point is calculated. The relative magnetic field intensity can be obtained from these calculated distances. Then the relative heat generated from induction heating can be reproduced from the relative power of each particular integration point, which is used as a coefficient, indicated as a red dot in Fig. 7(b). During the induction heat treatment simulation, the internal energy of each element is lowered to match the experimentally measured peak temperature at the surface of the shaft. In other words, the additional heat (ΔQ_ind) is added to reproduce the temperature increase during the induction heat treatment. The amount of heat added is set to be proportional to the relative power of integration point, which is a function of depth from the shaft surface.

Fig. 5.

(a) Finite element model (CAX4T, swept) and (b) cooling boundary condition during induction heating process. (Online version in color.)

Fig. 6.

Relative magnetic field intensity along the depth from shaft surface. (Online version in color.)

Fig. 7.

(a) Schematic diagram of shaft finite element. The distance from the surface to each integration point is calculated to obtain the relative magnetic field intensity. (b) An example of relative magnetic field intensity and relative power along the depth from the surface when the skin depth is 6.94 mm. (Online version in color.)

The cooling condition of the shaft during induction heating is complicated since the surface is cooled by the quenchant spouting from the backside of the coil, which is then further cooled by a quenchant flowing through the shaft. The thermal boundary conditions may be described as follows:   

$$-\lambda \frac{\partial T}{\partial n}=h_{Sprayed}(T-T_0)$$(19)

  

$$-\lambda \frac{\partial T}{\partial n}=h_{Flowing}(T-T_0)$$(20)

  

$$-\lambda \frac{\partial T}{\partial n}=h_{Air}(T-T_0)$$(21)

where n is the normal to the shaft surface and T₀ is room temperature (25°C). h_Sprayed and h_Flowing are the heat transfer coefficients of sprayed and flowing quenchant, respectively, and h_Air is the convective heat transfer coefficient of air cooling. The cooling boundary conditions during the induction heating process are shown schematically in Fig. 5(b) and the heat transfer coefficients are summarized in Table 7 using the proposed equations.^29,30) U is the moving coil velocity (7.1 mm/sec), which was originally proposed as a strip velocity on the run-out table (ROT) in hot strip mill.

Table 7. Heat transfer coefficients (W/m²∙°C) for each cooling condition.

$h_{Sprayed}=5\mathrm{\hspace{0.17em} }904{\left(\frac{T}{1\mathrm{\hspace{0.17em} }000}\right)}^{2.5}{\left(\frac{U}{20}\right)}^{-0.4}$

$h_{Flowing}=\frac{425\mathrm{\hspace{0.17em} }925}{{(T-95)}^{0.8}}$

$h_{Air}=30$

5. Results and Discussion

The induction heating process in the drive shaft was simulated with the model equations and numerical algorithms presented in the previous sections. Contour maps of temperature at certain times and phase distributions after the quenching process are shown in Fig. 8(a). As expected, the martensite phase is distributed on the shaft surface. The measured and simulated temperature history during the induction heating process is shown in Fig. 8(b) and temperature history together with the martensitic phase transformation along the profiling direction during induction heating is shown in Fig. 8(c). Since the lower measurable temperature limit of the IR camera is 200°C, the initial measured temperature was displayed as 200°C.

Fig. 8.

(a) Contour map of temperature at certain times and phase distribution after the quenching process. (b) Temperature history and martensite phase fraction during induction heating. (Online version in color.)

As the drive shaft is used under harsh conditions in automobiles, wear resistance is an important mechanical property. Thus, drive shafts are generally required to have hardness above a certain level. In this study, hardness was modeled based on a linear mixing rule for hardness of each phase as a function of alloying elements and cooling rates. The proposed hardness equations of each phase are:³¹⁾   

$$\begin{array}{l}{H}_{FP}=42+223C_F+53Si+30Mn+12.6Ni+7Cr\\ \hspace{1em}\hspace{1em}+19Mo+(10-19Si+4Ni+8Cr+130V)logV\text{'}\end{array}$$(22)

  

$$\begin{array}{l}{H}_B=-323+185C+330Si+153Mn+65Ni+144Cr+191Mo\\ \hspace{1em}\hspace{1em}+(89+53C-55Si-22Mn-10Ni-20Cr-33Mo)logV\text{'}\end{array}$$(23)

where H_FP and H_B are the hardness of a ferrite/pearlite mixture and bainite, respectively, V’ is the cooling rate at 700°C. and the element symbols above refer to the chemical composition of steel in wt.%.

The hardness of martensite, H_M is described as follows:³²⁾   

$$\begin{array}{l}{H}_M=3.03\times {10}^{-7}\cdot T_t(37.9+ln{t}_t)2\\ \hspace{1em}\hspace{1em}-0.0289\cdot T_t(37.9+ln{t}_t)+\left[874.9-202.2\cdot {\left(\frac{0.22}{C}\right)}^2\right]\end{array}$$(24)

where T_t and t_t are the tempering temperature and time, respectively. Considering the nonisothermal cooling condition of induction heating, Eq(H_M) was treated as a combination of a discrete isothermal tempering steps.³³⁾

The driveshaft hardness after induction heat treatment was calculated from the empirical equation above and compared with the experimentally measured data (Fig. 9). The simulated results agreed well with the measured data. In addition to the hardness profiling, axial stress is also superimposed in Fig. 9. As expected, it can be seen that the compressive residual stress is distributed in the shaft surface layer in which the martensite phase is distributed, whereas the tensile residual stress is distributed in the inner part of the shaft. For surface hardening heat treatments, the compressive residual stress is applied to the surface as in the simulated results, but since tensile residual stress is applied in the inner part of the shaft, it is important to investigate the stress distribution in the inner part in order to prevent fracture. Figure 10 shows the simulated and experimentally measured residual stress after the induction heating process. Since the measured stress component is the average value of the in-plane principal stresses, the same component of simulated residual stress was compared with measured stresses. The residual stress is partially relaxed by cutting the shaft. To estimate the amount of residual stress relaxation, the displacement during the cutting process was measured. This confirmed that stress up to 50 MPa was alleviated, and so this amount of stress was added to the experimentally measured stress in Fig. 10(b). As shown in the figure, the measured residual stress is still slightly below the residual stress predicted from the FE simulation. This difference might come from the fact that auto-tempering occurring during quenching was not considered in the suggested model. During the tempering process, at most 0.15 to 0.2 wt%C can segregate and if the M_s temperature is sufficiently high, part or all of this segregation can occur during quenching (i.e. auto-tempering).³⁴⁾ The mechanism of segregation is volume diffusion of carbon in martensite, which occurs at temperatures below 100 degrees. The segregation leads to volume shrinkage, which alleviates the residual stress generated during the cooling process.

Fig. 9.

Simulated (solid line) and measured (dotted line) hardness along the depth of the shaft. The axial stress (dashed line) along the same path is superimposed. (Online version in color.)

Fig. 10.

(a) Location of the indentation for residual stress measurement (red dotted positions) and (b) simulated and measured residual stress after induction heating process. The estimated relaxed stress from cutting process of shaft was considered. (Online version in color.)

In a hollow drive shaft, it is important to investigate the distribution of the phase at the inside of the shaft and the residual stress applied to the phase. Figure 11(a) shows axial stress and martensite phase distribution on the inner part of a drive shaft along the profiling direction (designated in Fig. 10(a)) under the various temperatures. As the heat-treatment temperature increases, the tensile residual stress where the ferrite phase is distributed increases as well. However, since the martensite layer has become thicker, the influence of the compressive residual stress is increased, and the stress changes drastically from tensile stress to compressive stress. It has been reported^35,36,37,38) that as the martensite fraction increases, the CVN impact energy decreases, meaning that the toughness is low and the part is vulnerable to fracture. Figure 11(b) shows the axial stress at the position with specific martensite fractions (50, 80, 100%) under various heat-treatment temperatures. The residual stress applied to each position is the lowest when the heat-treatment temperature is between 948 and 972 degrees.

Fig. 11.

(a) Axial stress and martensite phase distribution of inner part of drive shaft along the specified profiling direction under various temperatures and (b) axial stress at the position with specific martensite fractions (50, 80, 100%) under various heat-treatment temperatures. (Online version in color.)

6. Conclusion

A numerical model coupling the thermal, mechanical and metallurgical analysis was developed for investigating the phase distribution and residual stress under the induction heating process. The heat generated during heating was reproduced by considering the relative magnetic field intensity and skin depth. The kinetics of diffusional and displacive phase transformations was modeled from dilatometric analysis. The total deformation from phase transformation, transformation plasticity, and thermal deformation was implemented in the FE model. The calculated temperature history during the process matched the measured temperature well. However, the simulated results showed higher residual stress than the experimentally measured value, since auto-tempering arising during quenching was not considered in the proposed model. This model suggests that the optimal heat treatment temperature for the investigated drive shaft is between 948 and 972 degrees. These results show sufficient utility of the proposed numerical model in effectively predicting the thermomechanical behavior of automotive parts during induction heating.

Acknowledgement

This work was supported by the National Research Foundation of Korea (NRF) grants funded by the Ministry of Science and ICT (MSIT) (NO. 2015R1A5A1037627 and NO. 2018R1A2B6006856). The Institute of Engineering Research at Seoul National University provided research facilities for this work.

References

• 1)   J.  Davis and  P.  Simpson: Induction Heating Handbook, McGraw-Hill, London, (1979), 68.
• 2)   K. F.  Wang,  S.  Chandrasekar and  H. T. Y.  Yang: J. Mater. Eng. Perform., 4 (1995), 460.
• 3)   J. Y.  Kang,  S. J.  Park and  D-W.  Suh and  H. N.  Han: Mater. Charact., 84 (2013), 205.
• 4)   K.  Sadeghipour,  J. A.  Dopkin and  K.  Li: Comput. Ind., 28 (1996), 195.
• 5)   B.  Drobenko,  O.  Hachkevych and  T.  Kournyts’kyl: J. Eng. Math., 61 (2008), 249.
• 6)   J.  Yuan,  J.  Kang,  Y.  Rong and  R. D.  Sisson, Jr.: J. Mater. Eng. Perform., 12 (2003), 589.
• 7)   H.  Shokouhmand and  S.  Ghaffari: Appl. Math. Model., 36 (2012), 4304.
• 8)   X.  Bai,  H.  Zhang and  G.  Wang: Int. J. Adv. Manuf. Technol., 77 (2015), 717.
• 9)   C.  Liu,  Z.  Zhao,  D. O.  Northwood and  Y.  Liu: J. Mater. Process. Technol., 113 (2001), 556.
• 10)   J.-i.  Jang,  D.  Son,  Y.-H.  Lee,  Y.  Choi and  D.  Kwon: Scr. Mater., 48 (2003), 743.
• 11)   Y.-H.  Lee and  D.  Kwon: Acta Mater., 52 (2004), 1555.
• 12)   H.-J.  Ahn,  J.-h.  Kim,  H.  Xu,  J.  Lee,  J.-Y.  Kim,  Y.-C.  Kim and  D.  Kwon: Met. Mater. Int., 23 (2017), 465.
• 13)   J.  Davis and  P.  Simpson: Induction Heating Handbook, McGraw-Hill, London, (1979), 309.
• 14)   S.  Yafei,  N.  Dongjie and  S.  Jing: Proc. 4th IEEE Conf. on Ind. Electronics and Applications, IEEE, Piscataway, NJ, (2009), 368. https://doi.org/10.1109/ICIEA.2009.5138230
• 15)   H. N.  Han and  S.-H.  Park: Mater. Sci. Technol., 17 (2001), 721.
• 16)   Z.  Zhao,  C.  Liu and  D. O.  Northwood: J. Mater. Sci., 36 (2001), 5045.
• 17)   J.  Ågren: Scr. Mater., 20 (1986), 1507.
• 18)   H. N.  Han,  J. K.  Lee,  D.-W.  Suh and  S.-J.  Kim: Philos. Mag., 87 (2007), 159.
• 19)   H. N.  Han,  C. G.  Lee,  D.-W.  Suh and  S.-J.  Kim: Mater. Sci. Eng. A, 485 (2008), 224.
• 20)   H. N.  Han and  J. K.  Lee: ISIJ Int., 42 (2002), 200.
• 21)   D.-W.  Kim,  H.-H.  Cho,  W.-B.  Lee,  K. T.  Cho,  Y.-G.  Cho,  S.-J.  Kim and  H. N.  Han: Mater. Des., 99 (2016), 243.
• 22)   J. L.  Kuester and  J. H.  Mize: Optimization Techniques with Fortran, McGraw-Hill, New York, (1973), 386.
• 23)   H. N.  Han and  D.-W.  Suh: Acta Mater., 51 (2003), 4907.
• 24)   H. N.  Han,  C. G.  Lee,  C.-S.  Oh,  T.-H.  Lee and  S.-J.  Kim: Acta Mater., 52 (2004), 5203.
• 25)   J.  Miettinen: Metall. Mater. Trans. B, 28 (1997), 281.
• 26)   S.  Nanba,  M.  Katsumata,  T.  Inoue,  S.  Nakajima,  G.  Anan,  A.  Hiramatsu,  A.  Moriya,  T.  Watanabe and  M.  Umemoto: CAMP-ISIJ, 3 (1990), 871 (in Japanese).
• 27)   Y.-G.  Cho,  J.-Y.  Kim,  H.-H.  Cho,  P.-R.  Cha,  D.-W.  Suh,  J. K.  Lee and  H. N.  Han: PLoS One, 7 (2012), No. 4, e35987.
• 28)   S.  Koric and  B. G.  Thomas: J. Mater. Process. Technol., 197 (2008), 408.
• 29)   H. N.  Han,  J. K.  Lee,  H. J.  Kim and  Y.-S.  Jin: J. Mater. Process. Technol., 128 (2002), 216.
• 30)   H.-H.  Cho,  Y.-G.  Cho,  D.-W.  Kim,  S.-J.  Kim,  W.-B.  Lee and  H. N.  Han: ISIJ Int., 54 (2014), 1646.
• 31)   J. C.  Ion and  L. M.  Anisdahl: J. Mater. Process. Technol., 65 (1997), 261.
• 32)   J. H.  Hollomon and  L. D.  Jaffe: Trans. Am. Inst. Min. Metall. Eng., 162 (1945), 223.
• 33)   E.  Scheil and  W.  Hütter: Arch. Eisenhüttenwes., 8 (1935), 565.
• 34)   L.  Cheng,  C. M.  Brakman,  B. M.  Korevaar and  E. J.  Mittemeijer: Metall. Mater. Trans. A, 19 (1988), 2415.
• 35)   A.  Zare and  A.  Ekrami: J. Mater. Eng. Perform., 22 (2013), 823.
• 36)   B. C.  Kim,  S.  Lee,  N. J.  Kim and  D. Y.  Lee: Metall. Trans. A, 22 (1991), 139.
• 37)   M.  Kalantar,  H.  Najafi and  M. R.  Afshar: Met. Mater. Int., 25 (2019), 229.
• 38)   S. W.  Thompson: Mater. Charact., 136 (2018), 425.