Numerical simulation of grinding induced phase transformation and residual stresses in AISI-52100 steel

Compared with other machining processes, grinding operation has a very high energy density that generates very high temperatures in the ground zone. The temperature history is the key loading that induces a complex residual stress ﬁeld generated per dilatation and solid-state phase transformations, in addition to the normal and tangential mechanical loading. The various phase transformations encountered by steels when submitted to rapid heating and cooling have been modeled and implemented in a ﬁnite element (FE) model. The AISI-52100 bearing steel is taken here as reference material. The thermo-metallurgical and mechanical analysis has been performed using the commercial ﬁnite element software Abaqus s /Standard with various user subroutines developed to model the thermal, metallurgical and mechanical behavior of the material. The heat generated during grinding process was assumed as a moving heat ﬂux with elliptical distribution. The effects of the Peclet number and heat transfer coefﬁcient on the phase transformations and residual stresses have been analyzed. It was found that an optimal combination of grinding conditions could produce the desired magnitude of compressive residual stresses at the surface of the machined workpiece. It is also shown that omitting phase transformations could lead to a strong difference in the prediction of residual stresses.

Compared with other machining processes, grinding operation has a very high energy density that generates very high temperatures in the ground zone. The temperature history is the key loading that induces a complex residual stress field generated per dilatation and solid-state phase transformations, in addition to the normal and tangential mechanical loading. The various phase transformations encountered by steels when submitted to rapid heating and cooling have been modeled and implemented in a finite element (FE) model. The AISI-52100 bearing steel is taken here as reference material.
The thermo-metallurgical and mechanical analysis has been performed using the commercial finite element software Abaqus s /Standard with various user subroutines developed to model the thermal, metallurgical and mechanical behavior of the material. The heat generated during grinding process was assumed as a moving heat flux with elliptical distribution. The effects of the Peclet number and heat transfer coefficient on the phase transformations and residual stresses have been analyzed. It was found that an optimal combination of grinding conditions could produce the desired magnitude of compressive residual stresses at the surface of the machined workpiece. It is also shown that omitting phase transformations could lead to a strong difference in the prediction of residual stresses.

Introduction
Compared with other machining processes, grinding requires an extremely high energy input per unit volume to remove the surface layer of the material. Most of the energy is converted into heat which is concentrated in the grinding zone where the wheel interacts with the workpiece. This leads to the generation of nonuniform high temperatures, which may then result in solid-state phase transformations. The consequences of this rapid heating and cooling include the formation of a heat affected zone (HAZ), the generation of residual stresses, possible shrinkage or cracking of the material, often chemical modifications of the material [1][2][3][4], and therefore some variations of the material properties (like hardness, etc.). These phenomena may in fine play a role on the surface integrity of the component. Phase transformations in most steels introduce volumetric changes, transformation plasticity and changes in mechanical properties. Local plastic flow occurs when the effective stress exceeds the yield strength. All these factors interact with each other and eventually lead to a varying internal stress/strain field. The evolution of temperatures during grinding can be well predicted using models presented by several authors in [5][6][7][8][9], but how the phase evolution due to high temperature gradients during grinding process would affect the surface integrity of the material is not yet well understood. A comprehensive solution to this challenging problem would require a coupling of the thermal-mechanical analysis with the models which describe the kinetics of metallurgical phase transformations and the evolution of the microstructure, macrostructure, and residual stresses in the near-surface layer.
Numerical simulation of such a problem requires modeling of three different types of phenomena: thermal, metallurgical, and mechanical, which are, mostly, fully coupled. Many researchers [10][11][12][13][14] have proposed various models which account for all or part of these phenomena. However, the influence of metallurgical transformations remains a problem still widely open. The main objective of this work is to investigate the internal stress distribution and their evolution that occurs during grinding process of AISI-52100 steel which is one of the most commonly used steels in engineering components such as rolling bearings. Moreover, its properties are easy to acquire over a wide range of temperature. The FE model developed in this work takes into account the transformation strains associated with martensitic transformation along with the temperature dependent material properties. The variations of the residual stresses and strains at integration points have been considered, and the effects of the Peclet (Pe) number, non-dimensional heat transfer coefficient (H) and input heat flux (Q) over the temperature distribution and then subsequently on the microstructure and the residual stress state were analyzed. The Peclet number (Pe) is a convenient nondimensional quantity used to represent the velocity of the moving heat source considering the thermal properties of the conduction medium, which in turn determines the speed of dissipation of heat in the medium [15,16], Similarly H is a non-dimensional term used to express the relative heat coefficient of the cooling media for a grinding process considering thermal properties of the conduction medium and velocity of the heat source [17], where a w is the thermal diffusivity of the workpiece material calculated from the expression 2. Grinding process modeling Fig. 1 shows a schematic representation of a surface grinding process. The wheel center is moving at the constant velocity V w relatively to the workpiece -also known as the heat source velocity -while the rotation around its own axis results in the surface velocity V s . The depth of cut a is the depth of material removed in one pass of the wheel from the workpiece. L c is length of the contact between the grinding wheel and the workpiece. If the material removal process is not considered, then surface grinding can be essentially seen as a sliding contact between a cylinder and a half-space with a distribution of pressure, surface traction, and heat flux over the solid surface. In this first section the mechanical effects will be omitted and the analysis will focus on thermal and metallurgical effects.
A typical moving heat source with an elliptically distributed heat flux at the surface of a 2D semi-infinite solid (half-space) is schematically represented in Fig. 2. The contact length (L c ) between the wheel and the workpiece is generally assumed to be equal to the heat source length.

Coupling phenomenon
Several thermal, metallurgical and mechanical processes occur simultaneously during the grinding process throughout the average phase fraction of constituent i z eq i equilibrium fraction of phase i that is achieved after an infinite long time z g austenite volume phase fraction z M martensite volume phase fraction x, y x and y coordinates Dimensionless parameters used heating and cooling stages. In order to capture the residual stress state and resulting distortions, it is necessary to model all these phenomena as accurately as possible. However, depending upon the type of material, some simplifications may be adopted at this stage. For example, during grinding the regions where temperature values are high enough, some steels show phase transformations in solid state. In most steels the room temperature microstructure transforms entirely into austenite above Ac 3 and if the cooling rate is very high, the only phase obtained after complete cooling is martensite. Fig. 3 shows the interaction and coupling of thermal, metallurgical, and mechanical aspects of the process. In this work, some coupling phenomena like transformation plasticity and the effects of stresses over transformation were ignored.

Thermal model
The heat in grinding is generated by the friction between the grinding wheel and the workpiece. Heat conduction within the workpiece generates a temperature field with an extremely high thermal gradient in the close surface layer, which decreases further down into the material.
In order to compute the temperature histories, a heat transfer analysis is performed using the thermal properties of the material. The transient temperature field (T) in time (t) and 2D space (x, y) is achieved by solving the heat transfer equation: KðTÞ @T @x þ @ @y KðTÞ @T @y A schematic illustration of the heat generation during grinding is shown in Fig. 4.

Metallurgical model
The method for the phase transformation calculations during continuous heating from isothermal data is based on the rule of additivity as initially proposed by Scheil [19] as for cooling [20]. In this method the temperature-time curve is discretized into a series of isothermal steps. On each step the volume fraction of the new phase formed is calculated by using isothermal transformation kinetics. The isothermal transformation kinetics is modeled according to the law developed by Avrami [21]: where z i is the average phase fraction of constituent i at time t, and z eq i is the maximum or equilibrium fraction of phase i that a phase can achieve after an infinitely long time determined from the equilibrium phase diagram with known temperature and chemical composition. k i and n i are empirically obtained constants for the phase i. From the previous investigations [22,23], the constant k i generally depends on temperature, chemical compositions, and prior austenite grain size. Conversely, the parameter n i is usually known to be constant over a range of temperature and the suggested value varies between 1 and 4 [24].
The transformed phase fraction at the current time step n þ1 is calculated from an equivalent transformation time t 0 which produces the phase fraction on the previous curve of time step n(z i,n ) and the time step Dt.
At each temperature the coefficients k i (T) and n i (T) may be calculated by assuming two points corresponding to a given percentage of the phase formed. For example assuming that the  transformation curves of a continuous cooling transformation (CCT) diagram gives the start time t S when a small proportion (e.g. 1% ) of the new phase is formed and the end time t E when 99% of the equilibrium phase fraction z eq i is formed at a given temperature T.
Finally the parameters k i (T) and n i (T) are obtained by solving Eq. (8): The martensitic transformation is a displacive transformation only controlled by the temperature. The kinetics is given by the Koistinen-Marburger equation [25]: where, z M and z g are the martensitic and austenitic phase proportions, respectively, b is a material dependent coefficient; M S is the martensitic transformation start temperature, and T the temperature.

Mechanical model
After the thermo-metallurgical computations, temperatures and volume phase fraction become inputs for the mechanical simulation. The temperature variations and the phase transformations involve dilatational strains into the solid. The major remaining difficulty is to obtain the mechanical behavior of the mixture of phases. As a first assumption, the macroscopic behavior is supposed to follow an isotropic hardening with the von Mises criterion where the yield stress is obtained by a mixture of the yield stress of each phase. The flow stresses which take into account hardening, thermal and viscous effects follow a Johnson-Cook model as initially proposed by Umbrello et al. [26] and adapted to multiphase materials.
In the following, the total strain rate tensor is divided into a recoverable elastic part _ e e ij and an irrecoverable plastic one _ e p ij : _

Elastic strain
In the case of coupling of mechanical field with temperature and phase change, the relation for the elastic strain can be expressed as: where E is Young's modulus and v Poisson's ratio.

Yield function and plastic strain rate
The yield function criterion is expressed in stress space as where e p ij is the plastic strain, H m the hardening parameter, T the temperature, and z i the volume fraction of phase i. For sake of simplicity, F is a von Mises type criterion which is preferred for ductile materials such as metals.
where S ij is the deviatoric stress expressed by and s g y e pl ,T,z the global yield stress of the material which depends on strain, temperature and microstructure of the material and calculated by with s a y and s b y the yield stresses of the ferrite and austenite phases, respectively.
The equivalent plastic strain is given by

Flow stresses
For an accurate simulation of grinding, the mechanical behavior must take into account the hardening, thermal, and viscous effects. A usual way to do that is the use of the Johnson-Cook relation, well adapted for severe loadings. In [26], the authors propose a slight modification of Johnson-Cook relation for numerical simulation of hard machining of AISI 52100 bearing steel. The stresses are a product of three terms, representing the hardening curve at room temperature, the influence of temperature and the last one the visco-plastic part (Eq. 23). In [26], the temperature influence function z fact is interpolated as a 5th order polynomial function whose coefficients are given in Table 1. The strain hardening multiplier needs two other constants (m, A), also given in Table 1.
The previous single phase model must be adapted for multiphase materials. For sake of simplicity, the following two assumptions are made: (1) the flow curve for each phase follows an exponential law: (2) the mixture is calculated as a simple linear rule: s ¼ which is consistent with a Voigt hypothesis on the strains. Table 2 The stress-strain curves at room temperatures for each phases are presented in Appendix A.
Finally, the relation used is

Finite element simulation
The principle of numerical simulation of grinding entails the know-how of a comprehensive database with reference to geometry, thermo-mechanical properties, initial conditions, boundary and loading conditions. A brief description of the main features of the FE model is given below.

Finite element mesh
The mesh density is generally defined by the applied loading and/or boundary conditions. Since grinding processes involve high temperature gradient in and near the grinding zone, a very fine mesh is required to capture the temperature distribution in the contact area. As the temperature gradient becomes low far away from the grinding zone, a relatively coarser mesh is there sufficient for the analysis.
In this study the workpiece is considered as a 2D semi-infinite plate of 0.1 m length and 0.03 m width. The finite element (FE) mesh (Fig. 5) consists of CPE4T (4-node plane strain thermally coupled quadrilateral bilinear displacement and temperature) type elements totaling over 3216 nodes and 3000 elements with the smallest element in the mesh as [(5 Â 10 À 4 )(8 Â 10 À 5 )] m 2 .

Material: bearing steel AISI 52100
AISI 52100 (also known as 100Cr6 in Europe) is a high carbonchrome-manganese alloy steel which finds its applications in several rotating parts like anti-friction bearings, cams, crank shaft, etc. for its good resistance to corrosion and fatigue [27]. Compared to low-carbon steels, high-carbon steels can carry higher contact stresses, such as those encountered in point contact loading in rolling bearings [28]. The chemical composition of the AISI 52100 steel is listed in Table 2 and the key physical and mechanical properties of the material are given in Appendix A

Initial and boundary conditions
The initial temperature considered for the workpiece is the room temperature i.e. T(t¼ 0)¼ 20 1C.
On lateral faces, heat flux is imposed as linear convective transfer law: where T and T 0 are the temperature of the semi-infinite solid and the ambient temperature, respectively, and h conv (W/K m 2 ) is the convective heat transfer coefficient of the cooling media. Heat loss from the bottom surface was assumed to be zero i.e. q¼0. The thermal boundary conditions are schematically shown in Fig. 6.

Imposed heat sources
The thermal loading consists in applying a surface heat flux through a moving heat source. Jaeger [29] and Carslaw and Jaeger [5] have presented solutions for uniform moving rectangular heat sources and a uniform stationary heat source using the heat source method. The temperature distribution in a sliding contact was then estimated by several authors based on Jaeger's theory [30][31][32][33][34]. There are differing views among researchers on which distribution of heat flux is best to use for grinding. Some [7,35,36] have used a rectangular (uniform) distribution, so as to simplify subsequent calculations. However, due to the localized ''spike'' temperatures during a very short time, others [4,37,38] have argued that the assumption of a uniform heat flux field may not lead to accurate predictions. Keeping in view the contact origin of the heat source, theoretically the pressure and corresponding heat flux distribution -if one assumes a uniform friction coefficient in the contact area -should be modeled according to a sliding/rolling contact approach. Also by recalling that the Hertz contact pressure distribution between a cylinder and a plane is elliptical in shape, it seems reasonable to assume an elliptical distribution of the heat flux. Fig. 6 shows the presence of a schematic heat source moving with velocity V w on the top surface of the FE model. Here, the length of the heat source is equal to the contact length (2a ¼L c ) between the grinding wheel and the workpiece. The heat flux distribution entering the work piece is therefore given by where Q is the total heat per unit length (in W/m). In Abaqus s /Standard, the moving heat source is integrated in the finite element model through a FORTRAN subroutine, called DFLUX. The heat source is moving along the horizontal (x-) axis.

0.1m (200 elements)
Before carrying out a complete thermo-mechanical analysis, a few heat-transfer simulations were run to perform a sensitivity analysis of various types of heat source distributions. An elliptical, a triangular, and a uniform heat source were used. The resulting temperature distributions can be compared in Fig. 7. It was found that the peak surface temperatures in all three cases are very close, however, the distribution of temperatures over the surface varies to some extent. An interesting observation was that the flux distribution from elliptical source lies almost midway between the triangular and the uniform heat sources.

Implementation of the mechanical behavior of a multiphase material
As shown in Fig. 3, the mechanical behavior depends on the phase proportion but this coupling is not directly available with Abaqus s /Standard. Therefore a UMAT subroutine has been developed to go beyond this difficulty. The UMAT subroutine calls three other subroutines: PHASE, PROP and UEXPAN. For a given temperature field, PHASE compute the austenitic phase proportion formed during heating, and the martensitic phase proportion formed during cooling (see Section 2.3). Knowing the phase proportion, UEXPAN subroutine gives the expansion coefficient for the mixed material at a given temperature (see Eqs. (15) and (16)). Finally, PROP subroutine compute the material properties at a given temperature based on the fraction of phases. A linear mixture rule was used for the identification of the multiphase material properties.

Temperature distribution
The temperature variation in time (or profile along the horizontal axis) calculated at a given time instance (t 0 ¼33.61) and at various dimensionless depth (y/a¼2y/L c ) is shown in Fig. 8 for a specific set of parameters Q, L c , H and Pe. From the values of peak temperatures, it is found that up to a certain dimensionless depth (here 2y/L c ¼ 0.30) the temperature goes beyond Ac 1 and Ac 3 (750 1C and 815 1C [39] for AISI 52100, respectively). During cooling the temperature at these points will quickly drops below M s (250 1C [39]). It means that at high cooling the transformation of austenite to martensite will occur at the top most surface.
In Fig. 9 the maximum surface temperature as a function of the Peclet number is plotted for a given set of grinding parameters Q and L c as specified. It is shown that the peak temperature decreases when increasing the Peclet number. A comparison with the analytical solution of Blok [34] is also provided and a good agreement is found, which validates the numerical model. The effect of the variation of the dimensionless heat transfer coefficient, H, or the contact length, L c , is illustrated in Fig. 10. It can be observed that an increase of the heat transfer coefficient decreases the maximum temperature as an increase in the contact area does.

Phase transformations
As already mentioned, the most common transformation products that may be formed (from austenite) in the workpiece experiencing a critical grinding temperature history are (in order of formation with decreasing cooling rate): martensite, bainite, pearlite, ferrite and cementite. However, since both the heat source and convection plays  at the surface, each layer below the surface will experience a different heating and cooling history (assuming a 2D analysis and a steadystate regime of grinding). It is, therefore, understandable that at a certain depth below the surface, heating and cooling cannot activate phase transformation (see Fig. 8). The depth until which phase transformation takes place depends on the material properties and grinding conditions, including table speed, depth of cut, cooling rate and so on. Since only martensitic transformation has been considered in this work, it is subsequently assumed that the austenite resulting from heating of the parent phase (ferrite) transforms completely into martensite, i.e. no intermediate phases like bainite and pearlite, are present. The transformation of ferrite to austenite and then to martensite during a typical heating and cooling cycle with T max as high as 1000 1C is shown in Fig. 11. It can be observed that during heating the austenite transformation starts only at Ac 1 (750 1C) during which ferrite fraction reduces. At Ac 3 (815 1C), the transformation of austenite is completed and the ferrite vanishes. During cooling the austenite fraction remains constant up to M S (250 1C) and then martensite transformation starts which then continues up to M f (20 1C).
The effect of the dimensionless heat transfer coefficient on the martensite phase fraction is shown in Fig. 12. It may be observed that higher is the heat transfer coefficient (that increases the heat exchange by convection at the free surface), thinner is the martensite layer. Note that for H¼0 (i.e. insulated surface or grinding in vacuum environment) the maximum temperature reached is 1015 1C and an almost complete transformation of ferrite into austenite has occurred at the top surface layer. Cooling due to conduction then transforms entire austenite to martensite, thus yielding to nearly 90% of martensite for the top surface layer. At H¼0.2, however, strong convection is present that does not allow heating of the surface to very high temperature (T max ¼780 1C only). The resulting austenite phase fraction is, therefore, relatively low, and hence the transformed martensite fraction would also be also lower and the transformed layer thinner.

Evolution of stress and strain states
The mechanisms of residual stresses due to phase transformation can be best understood if the grinding stresses and strains are studied in relation to grinding temperature. Fig. 13 represents the evolution of the thermal strain as a function of temperature during the complete heating and cooling cycle. The thermal dilatation coefficient of ferrite (a a ) varies linearly with temperature up to Ac 1 (750 1C). Since body-centered cubic (bcc) austenite has its crystal lattice structure more densely packed than facecentered cubic (fcc) ferrite, lattice contraction appears between Ac 1 (750 1C) and Ac 3 (815 1C), and thermal strain starts decreasing with increasing temperatures up to Ac 3 . Beyond Ac 3 , the thermal dilatation coefficient of austenite (a g ) increases with temperature due to bcc lattice expansion. Upon very fast cooling, austenite survives up to M s (250 1C). Now, the bcc-austenite transforms into body-centered tetragonal (bct) martensite until M f (20 1C, assumed for simulation) is reached. Since the bcc-austenite to bct-martensite transformation results in the expansion of crystal lattice structure, the thermal strain increases with decreasing temperatures.
The history of the longitudinal stress s xx (z i ,T) for a surface element during a heating and cooling cycle when the maximum temperature reaches 1000 1C is plotted in Fig. 14. It could be noticed that with an increase in temperature the element tends to expand in all direction, however, restricted by its adjacent elements it develops compressive stresses which keeps on increasing with increasing temperature. These stresses are initially elastic in nature, but with further increase in temperature plastic strains appear and a stress reversal occurs. When the temperature exceeds the austenitizing temperature (Ac 1 ), almost a stress free state is observed because of the low yield stress of austenite phase at high temperatures. Having crossed Ac 3 , the element develops tensile stresses of low magnitude due to lattice contraction. When the heat source moves away the element starts to cool down, another stress reversal occurs due to contraction/ shrinkage and thus compressive stresses develop once again. Beyond M S , compressive stresses grow quickly due to lattice dilatation. This volumetric expansion finally results in very high longitudinal compressive stresses. The evolution of the longitudinal stress at the surface with and without phase transformation is illustrated in Fig. 15. With respect to the stress distribution during grinding, the surface can be conveniently sub-divided into three distinct zones. Zone 1 is immediately ahead of the heat source, where compressive stresses are present in both the cases due to the advancing heat front, as for welding. However, at some distance from the heat source the compressive stresses gradually reduce to zero. Zone 2 is immediately behind the heat source which is now subjected to cooling due to conduction and convection. Here, for the case with no phase transformation, the parent phase (ferrite) never transforms into any other phase, and hence the compressive stresses tend to become tensile in nature due to the contraction of the lattice throughout during cooling. The case with phase transformation, however, experiences lattice contraction and expansion as well as transformation in low yielding phase (austenite for example). The resultant is, therefore, a compressive stress profile due to the same reasons as explained above for Fig. 14. The presence of a free surface at the left side of Zone 3 acts in lowering the magnitude of the residual stress (whether in tension or compression). During phase transformation the surface stresses exceed the yield strength of the material, which is plastically deformed, resulting in thermally induced dimensional changes and surface hardening. It may, thus, be concluded that surface hardening would result in a higher magnitude of residual stresses regardless of their nature. Fig. 16 shows the residual stress profile along the depth with and without considering phase transformation. A very significant consequence of phase transformation is creation of high magnitude compressive stress, up to À 600 MPa, in a thin surface layer, whereas the stresses are tensile (up to 400 MPa) in absence of phase trans-formation. It can be also observed that the stress distribution far from the surface is hardly affected. This is attributed to the fact that phase transformation takes place only in the near surface layer, letting the subsurface unaffected. This also indicates that the distribution of residual stresses is directly related to the martensite depth. Fig. 17 illustrates the interdependence of residual stress level over martensite phase fraction as a function of the T max /T aus ratio. It is observed that the presence of compressive residual stresses is linked to presence of martensite. Moreover higher is the martensite fraction higher is the compressive stress magnitude. Note that in this 2D problem in plane strain the remaining stress components s yy and s xy are negligibly small and hence have not been presented.

Conclusion
The paper presents a FE model for the prediction of residual stresses due to grinding when thermal loading, normal and tangential mechanical loading, and phase transformation are considered. An analysis of the main parameters involved during a grinding operation and for a reference AISI 52100 steel has been carried out. The effects of the Peclet number and heat transfer coefficient over the temperature distribution have been first examined. The occurrence of the martensitic phase transformation and the thickness of the surface layer affected have been also analyzed. It was observed that the Peclet number and the heat transfer coefficient are the main parameters governing the onset of phase transformation in the ground components. It was also shown that during the grinding process, if the temperature exceeds the austenitizing temperature, high cooling rate will result in the formation of martensite which in fine will lead to compressive residual stresses which are directly related to the martensite proportion and affected depth. It may, therefore, be concluded that the type of grinding fluid, table speed and other grinding operation parameters should be carefully considered if the surface residual stresses and phase transformations in a ground component are the main concerns. In addition, a required level of residual stresses and martensite phase fraction could be achieved by controlling the grinding parameters.