3D modeling and simulation of thermal effects during profile grinding

A new heat transfer model for profile grinding was developed to analyze distortions caused by residual tensile stresses in linear guide rails. The simulative analysis of the thermal effects caused by a non-uniform heat source on the surface using the finite element method depends on an accurate representation of the locally variable contact area. The complexity of the V-groove profile disqualifies a 2-dimensional simulation approches so far used in the literature. This paper focuses on the redefinition of these mathematical relationships of the process parameters and the resulting heat flux. The heat flux model is adapted to the geometry of the workpiece depending on the grinding parameters and approximating the V-groove of a linear guide rail. This 3-dimensional modeling allows a better understanding of the thermo-metallurgical effects that occur during the grinding process. Furthermore, the calculation of the internal stresses induced into the workpiece material through the grinding process is possible. The simulation model results in a generally valid model for the analysis of distortions. In order to confirm the validity of the new heat flux profile, a comparison of the different finite element simulation results was made and experiments under wet grinding conditions were conducted. The results show that the newly developed grinding process model allows a more accurate prediction of workpiece distortion caused by grinding forces and temperatures. This research also offers a new approach to a method based on a 2-dimensional implementation developed in the literature for predicting the distortions of linear guide rails and a derivation of possible simulation-based compensation strategies.


List of symbols A c
Total contact area between workpiece and wheel a e Depth of cut c c Specific heat capacity of the grinding wheel material c w Specific heat capacity of the workpiece material d s Grinding wheel diameter d(x) Local grinding wheel diameter in x-direction F n Normal grinding force component F t Tangential grinding force component F ti Differential tangential force at l i h Heat transfer coefficient h max Maximum height of the V-groove circular segment h(x) Local height of the V-groove circular segment in x-direction according to the variable radius k c Thermal conductivity of the grinding wheel abrasive layer k w Thermal conductivity of the workpiece material l g Maximal geometric contact length l g (x) Local geometric contact length in x-direction P Total grinding power P i Power produced at the location l i q Total heat flux q w Heat flux transferred into the workpiece s Normal cutting depth r Variable V-groove radius at the deepest point v c Cutting speed v fluid Fluid flow rate v ft Tangential feed speed Angle of the side surface of the V-groove (x) Local angle of the circular ground surface in the V-groove Heat partition ratio between heat flux into the workpiece and total heat flux c Mass density of the grinding wheel abrasive layer 1 3 w Mass density of the workpiece material Q ′ w Specific material removal rate q ti Heat flux generated in segment l i l i Segment i of the geometric contact length l g

Introduction
The machining process induces a change of the material properties in a defined layer It can be characterized by a surface and subsurface effect on and below the newly created surface [1]. During the grinding process, the subsurface is significantly affected by external mechanical, thermal and chemical influences on the material. The mechanical and thermal influences cause a change due to the forces and temperatures occurring during grinding in the area of the inner boundary layer and can have a structural depth effect up to several millimeters. This significantly changes the properties of the material's boundary layer, such as hardness, microstructure and residual stress [2]. During machining with geometrically undefined cutting edges, the material is first elastically and plastically deformed due to the mechanical load up to the chip shear [3]. The mechanically induced plastic deformation can result in a hardening of the subsurface area and an induction of residual compressive stresses. This is due to the fact that grain cutting edges penetrating the material initially cause elastic deformation and then plastic deformation as the cutting edge engages, causing the material to be stretched parallel and orthogonal to the cutting direction and thus subjecting it to tensile load stress [4]. After the cutting edge engagement, the external load on the material decreases and the elastic deformations recede, so that residual compressive stresses result in plastically deformed material.
The thermal load on the edge zone material during grinding is basically dependent on how the heat generated in the process is distributed among the active partners (grinding wheel, workpiece, cooling lubricant and chips) and how it is partitioned in the contact zone. The heat distribution is influenced significantly by the grinding conditions, the abrasive, the material to be machined and the type and supply conditions of the cooling lubricant, resulting in a range of the conduction of the heat via the workpiece of 5-84% [5]. This can lead to high local temperatures in the surface and subsurface of the workpiece, which can cause changes in hardness and microstructure. The influence on the properties of the surface and subsurface regions is then manifested by tempering and new hardening zones (white layers), so that the residual stress state is influenced as a result of the volumetric changes in the material microstructure. In the case of a inhomogeneous distribution, this results in changes in dimension and shape of the workpiece [6]. In general, a distinction must be made between thermal and transformation stresses when thermally induced internal stresses are generated: If the material cools down from the workpiece core towards the surface, tensile residual stresses occur in the surface and subsurface areas, assuming no structural transformation is achieved [7]. When an external thermal load is applied, compressive load stresses are initially induced due to the volume expansion caused by the heat, and are then reduced by plastic compression of the material when the yield point is exceeded. During cooling, the material contracts and tensile residual stresses result in the plastically compressed area [8]. Both experimental and simulative observations of the heat generated during the profile grinding process are still a challenging issue in academic and industrial research. Up to now, mainly analyses and modeling of surface grinding processes have been conducted. However, the models cannot be simply transferred to profiled workpieces. Linear guide rails with their V-profiles are in particular often finished with a profiled grinding wheel to improve surface quality and geometric properties.
There is also an influence on the boundary zone that results in a distortion of the workpiece. Decisive influencing factors are the shape of the heat flux, the maximum temperature reached in the workpiece and the distribution of the temperature field [9]. The most critical point in describing these factors, especially with regard to modeling the process to predict distortions, is the exact local calculation of the contact length for different depths of cut, which are responsible for the specific energy absorptions and grinding heat flux in the grinding zone [10]. This paper focuses on the calculation of the contact length, on the construction of an analytical heat source model for the three-dimensional machining of the profiled workpieces and on the application possibilities for the prediction of the workpiece distortion. For this purpose, a newly developed model for the calculation of the local contact length is combined with the previously used theoretical models from Jaeger [11]. Experiments in the literature have shown that the heat flux distribution along the feed direction is triangular due to the proportional increase in the material removal rate [12]. These validated results regarding the heat flux were taken up and combined with the local variation of the contact length [13]. Furthermore, Guo and Malkin describe how a segmented energy distribution within the grinding zone can be calculated for face grinding [5]. This can also be transferred to a workpiece with a V-groove by modifying the described workpiece geometry.
Many of the existing research activities focus on the description of the heat flux distribution patterns of the grinding process for simple geometries, which can usually be reduced to two-dimensional simulation models [14]. The lack of three-dimensional models of profile grinding processes makes thermal analysis difficult, especially with regard to the occurrence of distortion. The following analytical results serve as a starting point for the implementation of such a grinding model.

Material and geometry
The finite element method (FEM) offers the possibility to map the initial state and the grinding process steps one after another by using several connected simulation studies. Within the scope of the research work described here, the simulation software COMSOL Multiphysics from COMSOL Inc. was used to implement the developed analytical model. Both, the thermal and mechanical influences can be simulated in a 3D model to compare the resulting deformations with the real material behavior. The steel 42CrMo4 (AISI 4140) with the dimensions 23 mm × 38 mm × 250 mm (h × w × l), which is frequently utilized in manufacturing linear guide rails, was used as the initial material. The composition and phase change properties are shown in Table 1. To obtain a reference free of the residual stress, a soft-machined and heat-treated [QT 200 • C (55HRC)] steel was used. A right-angled V-groove (h × w = 8.67 mm × 19 mm, 2 mm radius in the groove base) was ground centrally over the entire length of the workpiece. The workpiece was clamped on a magnetic clamping plate for the grinding process.

Mathematical heat flux model
Crucial for the modeling of the mechanical forces acting on the grinding wheel is the calculation of the stress applied to the surface for a two-dimensional, linear-elastic, isotropic material model. The third dimension can be disregarded due to the symmetry of the workpiece and the absence of the deformations in x-direction. COMSOL Multiphysics allows the calculation of the solutions of the FE model with reference to the feed speed v ft via a time-dependent study. Measurements can be used to determine the typical load values for the tangential component in z-direction and the normal component in y-direction.
The geometric contact length l g (Eq. 1), on which the grinding wheel applies its mechanical load within the V-groove, is calculated from the following relation according to the Jaeger temperature model for an FE model [16]: The variables a e and d s indicate the depth of cut and the diameter of the grinding wheel. This description serves only as an approximate initial assumption for a surface grinding process that can be represented two-dimensionally. Due to the V-groove for linear guide rails present here, the model must be extended. Both the local depth of cut perpendicular to the material and the local diameter of the grinding wheel decrease from the lowest point of the V-groove to the outermost point. Figure 1 shows a representation of the defined depths of cut and angles. The angle is the tilt of the side surfaces in the V-groove in relation to the horizontal surfaces. The locationdependent cutting depth perpendicular to the side surfaces is described by s . Due to the exact adaptation of the groove to the shape of the grinding wheel for a simplified and at the same time more accurate simulation, the side surfaces and the groove base with radius r can be regarded separately in this geometry. For the following calculations, the coordinate origin lies exactly at the deepest point of the groove. The  for |x| > r ⋅ sin ( ) on the side surfaces and for |x| ≤ r ⋅ sin ( ) at the vertex area.
In addition to the depth of cut, the grinding wheel diameter is also relevant for calculating the location-dependent contact length. The geometrically determined reduction of the grinding wheel can again be determined by the height h of the circular segment in the vertex area, depending on the radius r and length of the segment, and by the angle of the groove side surfaces (pictured in Fig. 2).
Equations 4 and 5 provide the height h of the circular segment: and Starting from the maximum grinding wheel diameter d s , the diameter can be reduced depending on the position x. The local diameter d(x) is described by Eqs. 6 and 7: for |x| > r ⋅ sin ( ) on the side surfaces with maximum circle segment height h max and for |x| ≤ r ⋅ sin ( ) at the vertex area. By inserting Eqs. 2, 3, 6 and 7 into Eq. 1, one obtains the locally varying contact length l g for the grinding wheel contact with the side surfaces and the groove base: for |x| > r ⋅ sin ( ) on the side surfaces and for |x| ≤ r ⋅ sin ( ) at the vertex area.
These calculations are valid for analyses of the contact lengths during grinding for any side surface angles and groove base radii (Fig. 3). In the experiments, a radius of r = 2 mm and a maximum grinding wheel diameter of d s = 400 mm were used for the angle = 45 • . Based on this, further considerations can now be made to calculate the heat flux into the workpiece. Rowe describes the total heat flux over q w as a function of the grinding power P and the contact area A c [17]: The total contact area A c is the integral of the contact length l c (x) over the width of the surface in the V-groove.
Here it must be considered that the contact length integrated over the width of the V-groove, which is dependent on x, does not take the whole groove shape into account but only the length in x-direction. The function of the contact length must first be projected onto the length of the V-groove surfaces. After the projection, the resulting function can be integrated and provides the actual contact area. The total heat flux can then be calculated according to Eq. 10 for the corresponding depth of cut a e . The grinding forces measured by means of the force measurement platform, that was used during the experiments for each grinding process, are also necessary for the calculation of the grinding power P and thus for the total heat flux [18]. P is determined by Since the partition of the total heat flux flowing into the workpiece is approximately constant for a grinding process with cooling lubricant in the pores of the grinding wheel, the following Eq. 12 can be used to determine the average energy partitioning [19]: In addition to the already defined velocities of the grinding wheel and the workpiece, the specific heat, thermal conductivity and density of the grinding wheel material (c) and the workpiece material (w) are used here. The average heat flux per unit area q w = P A c going into the workpiece is derived from Eqs. 11 and 12. Table 2 shows the measured forces, the grinding power and the calculated contact areas as a function of the depth of cut and the parameters used in Table 3.
Due to the given shape of the workpiece, the total heat flux cannot be distributed evenly over the contact surface. To ensure that the subsequent heat flux simulation using FEM is successful, the heat source model from Jiang et al. was used [20]. Their model describes the division of the contact surface along the z-direction into n segments of the same length and calculates the corresponding heat flux for each of these segments. In this way, the segmented power is calculated with a contact area split into i segments with 0 ≤ i ≤ n: Thus, the heat flux, which is generated in the segment l i , is calculated according to Jiang et al. [20] using  It is assumed that the model of Guo et al. provides the average energy partitioning and its proven triangular heat flow model (pictured in Fig. 4). This model, previously designed only in two dimensions along the feed direction (z-direction) for a constant contact length l g , was transferred and designed by segmenting the whole contact area, in particular the locally varying contact length l g (x) , in the V-groove onto the profiled component, which was transferred to the simulation as described in Sect. 3.

Numerical simulation
COMSOL Multiphysics offers the possibility to simulate the initial state and the grinding process. The analytical model was transferred to an FE simulation and simulated with reference to all relevant parameters. For this purpose, it was necessary not only to determine the parameters used, such as feed speed, cutting depth with corresponding specific material removal rates Q ′ w (Fig. 5), cutting speed and coolant flow, but also to measure other occurring parameters such as the forces and the power during the experiments. By means of the measurements of a force-measuring platform, on which the workpiece rests during grinding, typical load values for the tangential components and the normal component in yand z-direction can be determined.
The starting material is AISI 4140 and consists of pure, low-stress martensite at the beginning of the simulation. This is in accordance to the experiment. Due to the wide temperature range that is recorded during grinding, an exact representation of the resulting microstructural transformations is necessary. Thus, the austenitizing temperature of (14) q ti = P i bl g (x)∕i . 723 • C is partially exceeded. The result is a microstructural transformation within the upper material layers. Diffusioncontrolled transformations are produced, for example from austenite to bainite, which can be mapped by the Johnson-Mehl-Avrami algorithm [21]. Due to the rapid cooling caused by the cooling lubricant, the formation of a layer by means of the non-diffusion controlled lattice shear is dominant at the surface. Therefore, the Koistinen-Marburger algorithm can be used for the determination of the temperature and the microstructure-dependent material properties [22]. For the calibration of the material typical timetemperature transformation (TTT) algorithm, continuous cooling (CCT) and time-temperature austenitization (TTA) diagrams are used [23]. Based on these and the percentage calculations of the microstructures, the relevant 42CrMo4 material parameters for the FE simulation can be calculated. For this purpose, the temperature gradients in the material determined from heating and cooling phases must be taken into account. This procedure allows a detailed analysis of the regions and microstructure changes critical for the distortion formation. Apart from the phase transformations, which were imported via the COMSOL material interface, the symmetry conditions were exploited before meshing. For this purpose, a symmetry surface was defined as a boundary condition after dividing the workpiece geometry along the deepest point of the V-groove in z-direction using both the 'Structural Mechanics Module' and the 'Heat Transfer Module'. Points at the upper end edges were defined as fixed constraints to simulate the unclamping from the magnetic clamping plate after the grinding process and to ensure free deformation. Essential for the deformation is the definition of a plasticity model accurate to the material. This model was added as an additional module to the ' Structural   Fig. 4 The sketch of the meshed FE model, the triangular heat flux and the segmentation along the z-direction normal to the grinding surfaces within the V-groove Fig. 5 Depths of cut with corresponding specific material removal rates Q ′ w for a tangential feed speed of 1500 mm/min Mechanics Module'. As assumptions proven in the literature, 'small plastic strains' with a von Mises stress yield function were used as assumptions. For the simplified linear isotropic hardening model, the values for the initial yield stress as well as for the isotropic tangent module (about 190 GPa for highstrenght tempered steel) were taken from the material database and from the experimental results of Ellermann [24]. The defined linear-elastic isotropic material together with the modulus of elasticity and the coefficient of thermal expansion provides the basis for solid state mechanics. In this material xy as shear stress, the maximum and minimum s t r e s s i n t h e x -y -p l a n e i s g i v e n b y xy with x the stress state in x-direction and y the stress in y-direction [25].
Beside the 'Structural Mechanics Module' the 'Heat Transfer in Solids Module' was linked via the multiphysics function. The model for a 3D heat source, which was analytically derived in the previous chapters, was calculated from the force and power measurements of the real-time experiments before its implementation in COMSOL. The energy partitioning was then performed by building a triangular analytical ramp function that moves along the z-direction through the V-groove as a function of time, thus mapping the thermal process of grinding according to constant feed and cutting speeds. The model, which has been extended here to three dimensions, enables a more detailed distortion analysis for profiled components. Up to now, grinding processes have only been reduced to a 2D model, allowing an approximation of the workpiece distortions. The calculated heat flow for the used depths of cut flows into the workpiece as a time-dependent general inward heat flux. The heat flux was determined proportionally via . In addition to the heat flux conditions, further boundary conditions caused by the use of the cooling lubricant are also important. Parallel to the pure heat input generated by the deformation and friction on the surface, the cooling effect of the coolant is taken into account by defining a heat transfer coefficient, based on the formula of Hadad, of h = 24 ⋅ 10 3 W m 2 K on the side surfaces and of h = 6 ⋅ 10 3 W m 2 K behind the grinding wheel within the V-groove [26]. An overview of the simulation properties used can be found in Fig. 6.
Meshing was performed using tetrahedral elements. The tetrahedron density is very fine in the contact zone, i.e. within the V-groove, and coarser on the side surfaces in order to save computing time and to achieve a good element ratio. After meshing and modeling of the solid state mechanics and heat flow, a PARADISO solver was used via the time-dependent study function, which allows COMSOL to display the corresponding temperature and stress distribution at each time step.

Shape measurements
To demonstrate the validity of the developed analytical model, the results of the FE simulation were compared with experimental data. For this purpose, the workpiece was measured once before and once after the grinding process. In each case, 13 measuring points were defined at 20 mm intervals on the underside of the workpiece on two lines along the entire length of the workpiece. These lines are each located 4 mm from the side surfaces. A coordinate measuring machine with a precision error of about 1-3 μ m scanned the position of these defined measuring points. Three measurements were made for each parameter set, which allowed averaging over both the three samples and the two defined measurement lines. Thus, the measurements before and after the grinding process resulted in an averaged peak-to-valley difference PV. This difference describes the Fig. 6 The defined conditions of the numerical simulation of the grinding process distortion caused by the tensile stresses introduced due to the heat input during grinding [27]. Due to the short minimum distance of 9.5 mm from the V-groove flank to the vertical side surface of the workpiece and the strong flow of the cooling lubricant, the integration of thermocouples was not possible. These would have led to a considerable falsification of the temperature and the distortion measurements, since the relevant tensile stress induction and the heat flux occurs close to the V-groove surface. The inclusion of thermocouples leads to a thermal scattering and thereby to an unintended symmetry disorder. As a consequence, further tests were aimed exclusively at distortion analysis, distortion prediction and possible distortion compensation of profiled ground workpieces. By using the calculated analytical heat flux model for a variation of different depths of cut, as already shown in Table 2, this distortion can be simulated according to Sect. 3 based on COMSOL Multiphysics.  [28]. The simulated and measured data on form deviation agree well. All simulated results for the distortion are within an acceptable error range of the experimental data. If the depth of cut is raised, the form deviation increases. This is due to the increase in grinding power or the locationdependent contact length and thus an increase in heat flux into the workpiece. This leads to an increment in the induced tensile stresses on the workpiece surface. It is noticeable that for cutting depths of higher than 0.6 mm the simulated distortions are lower than those measured experimentally. For example, at a e = 0.7 mm, an experimentally determined distortion of 195 μ m was measured, whereas the simulation shows a distortion value of 185 μ m. This corresponds to a difference of 10 μ m and can be explained by the fact that the FE analysis considered mainly the thermal influences during the process. These have the greatest influence on the material structure at low cutting depths. At higher cutting depths, on the other hand, the mechanical load increases directly in the grinding zone, resulting in less distortion in the simulated results. Due to the additionally increased mechanical interactions between the abrasive grains and the workpiece at higher depths of cut, locally limited plastic flow leads to increased residual compressive stresses. In addition, with the thermal compressive stresses generated near the surface, an increased plastic flow under pressure is caused, which is partly restricted by thermal expansion of hotter material closer to the surface and partly by cooler material below the surface. After the subsequent cooling, strong tensile residual stresses occur at the highly compressed surface increased by the mechanical load of the abrasive grains. These are further increased by the metallurgical transformations occurring during the heating and cooling cycle [29]. The lower distortion in the simulation is also due to the fact that the plate is elastically bent in the finite element model due to Fig. 7 Comparison of the measured and the simulated distortions for specific material removal rates Q ′ w defined in Fig. 5 with a tangential feed speed of 1500 mm/min the specific source force. The source force perpendicular to the cross-section of the plastically deformed surface layer as a measure of the force per unit width of the profiled workpiece allows predicting the shape deviations of hardened workpieces with simplified geometry. Based on the ReiSSner Theory, source stresses and source forces can be used to determine proportional form deviations which occur during grinding as input for a simple linear-elastic finite element simulation of complex workpieces [30]. The cross-section is deformed proportionally to the absolute value of the specific source force. On the other hand, the calculation of the specific source force from the beam theory confirms the measured peak-to-valley values. According to the beam theory, the mechanical energy is stored by the bending deformation. Less stored energy due to the bending deformation results in less deformation in the simulation. Because of this, the difference becomes larger with increasing source forces. However, the magnitude of the distortion is well predicted with the implemented heat source model. In addition to the direct distortion data, thermomechanical simulation also provides the temperature distribution at a certain point in time during the simulation. This temperature distribution is specific for the developed analytical 3D heat flux model. On the one hand, it can be read for any temperature point within the material if the critical cooling time has been exceeded, on the other hand, it can also be determined if the austenitization temperature has been exceeded for a certain depth in the workpiece.

Validation of the simulation approach
The maximum temperature reached on the V-groove surface of T max = 729 • C for a cutting depth of a e = 0.8 mm can be seen in Fig. 8. This corresponds to a normal cutting depth of s = 0.57 mm. Measurements were taken at different, evenly spaced depths on a V-groove side flank point.
Looking at the shapes of the temperature curves, they are almost identical to the ones obtained by Jiang et al. [20]. The comparison can be made to that extent since the grinding of the side flanks of the V-groove is equivalent to a plane grinding with a smaller cutting depth. First, the temperature remains constant at room temperature until the grinding wheel passes over the measuring point and the maximum temperature is reached directly in the contact zone. The material subsequently cools down much faster directly on the surface than in the inside. This is caused by the massive circulation of the cooling lubricant around the workpiece. Estimations can also be made from the simulations for the temperature curves shown here. In contrast to simulations with a uniformly distributed heat source, according to literature and previous research results, a triangular heat source can be used to simulate much more precise temperature curves over time, which more accurately represent the real process [31]. For this reason also for the here computed temperature courses by the simulated triangular heat source, which was described in the previous chapters, errors of under 6% can be determined. The error is mainly influenced by the existing geometry. The simulation accuracy of the model suffers from non-planar surfaces of the workpiece. However, with the model of variable contact length presented here, the error, which depends on the proportionality to the square root of the contact length, has been improved.

Summary
This paper presents the research results of an extended threedimensional heat source model and a thermometallurgical FE simulation. Both lead to a better understanding, modeling and prediction of distortions during profile grinding of long steel components. Decisive for the distortions are the mechanical loads on the microstructure and the heat flux on the uppermost material layer. The more precise analysis of the V-groove geometry with regard to final finishing processes of linear guide rails made it possible to extend the surface grinding models by calculating the local contact lengths for different depths of cut. With the analytical model, it is now possible to build a thermal model for V-grooves with any flank angle, replacing the 45 • used here. Based on these models, the workpiece distortions can be simulated more accurately. The comparison of the experiments with the results of the distortion simulation based on the new heat flux model and the phase transformation effects has shown that accurate distortion predictions can be made. A precise analysis of the stresses and the extended parameter variations will now form the basis for a subsequent more detailed distortion analysis and for distortion compensation strategies by externally forced introduction of the compressive and tensile residual stresses. If relevant parameters and values are compared with the simulations by means of the Fig. 8 Temperature simulated with the new analytical heat flux distribution model for a typical profile grinding experiment with a depth of cut of a e = 0.8 mm experimental validation, detailed statements can be made about temperature curves and stresses. Possible adjustments can be made and investigated in the model. In addition, the experimental determination of the tangential contact load in the contact zone offers the possibility to draw conclusions about the heat flux density distribution. By assembling simulated partial models based on thermal and mechanical pre-or post-treatments, the complete manufacturing process consisting of the distortion and the distortion compensation will be represented in the future work of this research project.