Analytical Prediction of Residual Stress in the Machined Surface During Milling

An analytical prediction model for residual stress during milling is established, which considers the thermal-mechanical coupling effect. Considering the effects of thermal-mechanical coupling, the residual stress distribution in the workpiece is determined by the stress loading history according to McDowell′s hybrid algorithm. Based on the analysis of the geometric relationship of orthogonal cutting, the prediction model for milling force and residual stress in the machined surface is established. The research results can provide theoretical basis for stress control during milling.


Introduction
The residual stress in the machined surface is a self-balanced internal stress. It is the internal stress remaining in the object after the elastoplastic object undergoes a large deformation and restores the external load, torque and thermal gradient to the initial state of the object. For components with weak rigidity, with the release of residual stress, the component will have obvious bending or torsional deformation, and the fatigue life of the component under cyclic loading will be reduced. The degree of effect of residual stress on the performance of components is generally judged by three indicators, that is, the residual stress property of the workpiece surface, the peak value of the residual stress and its corresponding depth on the subsurface. For the cutting process, the residual stress caused by plastic deformation, non-uniform heating and microstructure transformation are mainly studied.
The analytical method is efficient in calculation and can explain the physical mechanism clearly. Therefore, the research and exploration on the analytical algorithm of residual stress have become a difficult problem that many scholars insist to overcome. Su et al. [1] proposed an analytical model of residual stress during milling considering the corner radius under minimal lubrication conditions. The orthogonal model was applied to the more complex milling prediction by geometric transformation. Fergani et al. [2,3] established an analytical model for residual stress regeneration in milling and predicted the residual stress distribution in multi-pass milling and proposed an analytical model to predict the deformation of thin-walled parts caused by residual stress. Wang [4] established an analytical model of residual stress for flank milling. The model applied the radial return method to update the plastic stress components during the plastic loading process. Finally, the measurement experiments for milling temperature and residual stress are performed to validate the analytical model. It is found that the plowing effect of the cutting edge on the workpiece is the major source of the residual stress. Wan et al. [5] considered the 3D instantaneous contact status between the tool and the workpiece to predict the residual stress during milling. A modified analytical model for predicting the milling temperature field was presented. Huang and Yang [6] proposed a criterion for determining the initial stress state of the workpiece based on the calibration of the residual stress and achieved accurate prediction of the residual stress distribution of workpieces of different materials. Huang [7] proposed an influence coefficient method for quickly calculating the residual stress distribution based on the eigenstrain theory. The analytical model established by this method can predict the periodic residual stress distribution under dynamic cutting conditions. The analytical method has high-efficiency prediction and can directly reflect the physical mechanism of cutting, which can effectively promote the development of finite element methods for residual stress. In addition, it can guide the direction of experimental research and establish the relationship between residual stress and other data that can be easily obtained by experiments.
The residual stress is greatly affected by the cutting forces and temperature which they influence each other. However, the milling process has the characteristic of discontinuity, the force and heat acting in the surface of the workpiece are two complex amounts that change periodically, it is difficult to consider the coupling of milling force and milling heat. In this study, iterative algorithm of thermal-mechanical coupling in orthogonal cutting is introduced, and the prediction of residual stress in milling process is simplified to takes both mechanical effect and thermal effect into consideration. The research objective is to accurately predict the peak value of residual stress on the subsurface of Ti6Al4V, and realize the prediction of residual stress using MATLAB® based program (The software version is R2016a, Natick, MA, USA). The main innovations of this paper are as follows: (1) Considering the thermal-mechanical coupling, a milling force model is established. Based on the thermo-mechanical coupling algorithm, the relative motion relationship between the tool and the workpiece during milling is analyzed, and a milling force model considering the thermalmechanical coupling effect is established. (2) Based on the milling force model, the residual stress prediction model is established.
Considering the surface waviness of the workpiece after milling, an approximate assumption is proposed to planarize the machined surface to adapt to the residual stress analysis algorithm.

Research on Iterative Algorithm of Thermal-Mechanical Coupling in Orthogonal Cutting
Ji Xia [8] proposed an iterative algorithm of thermal-mechanical coupling in orthogonal cutting. The cutting force predicted by the Oxley's cutting force model is used as the input of the cutting temperature prediction model to obtain the cutting temperature, and then put it into the Oxley's cutting force model to obtain a new cutting force. After several iterations, when the output cutting force converges to a certain value, at this time, it is considered that the cutting force and the cutting temperature have reached a balance, which is the process of thermal-mechanical coupling.
The Oxley's cutting force model considers that the strain rate on the shear zone is almost constant. This view was confirmed by Stevenson and Oxley [9]. Furthermore, it can be considered that the flow velocity VS on the shear zone is constant. In addition, it is assumed that the flow stress on the shear zone is uniformly distributed at a certain instant. So the cutting force can be calculated based on the uniformly distributed stress on the shear zone, and the cutting mechanism is analyzed based on the geometric relationship. The assumptions of the model are as follows: (1) There is no tool wear and cutting vibration, and no built-up edge is generated.
(2) Plastic deformation in the cutting layer only occurs in a plane perpendicular to the cutting edge.
(3) The primary and secondary deformation zones are assumed to be narrow area. (4) The temperature and strain on the shear plane are uniformly distributed.
The geometric relationship of force and material flow speed in three different directions (along the cutting direction, along the shear plane AB and along the tool-chip interface) is shown in Figure 1. The resultant cutting force Ptotal can be expressed by the flow stress σAB on the shear plane AB by the geometric relationship as: where hc is the cutting thickness, w is the width of the cutting edge participating in cutting, Φ is the shear angle of the primary shear zone, and θ is the angle between the force Ptotal and the shear plane AB. The friction F and normal force N at the rake face of the tool, the force Fc in the cutting direction and the force Ft in the feed direction can be expressed as: where α is the rake angle of the tool and λ is the friction angle of the rake face. Shear flow stress σAB can be obtained from various constitutive models, such as the Johnson-Cook constitutive model, Zerilli-Armstrong constitutive model, Bammann-Chiesa-Johnson constitutive model, etc. Hyperbolic tangent (TANH) constitutive model is applied to predict the shear flow stress at shear zone. The constitutive model is proposed by Calamaz [10]. It considers the strain softening effect of Ti6Al4V on the basis of J-C constitutive model. According to the study of Ducobu [11], the constitutive parameters are shown in the Table 1.
force model include only material parameters, process parameters and tool geometry parameters, which improves the predictive efficiency. The model developed by Huang and Liang [12] is used in the present study to predict the cutting temperature. This model considers the combined effect of the heat source in primary and secondary deformation zones, the non-uniform heat distribution ratio along the chip contact surface to obtain the temperature along the tool-chip interface. In the present study, the temperature generated by the friction in the tertiary deformation zone is considered. The heat source in the tertiary deformation zone is calculated based on the model proposed by Su [1]. It is assumed that the heat distribution ratio in the tertiary deformation zone is constant. The temperature of chips and tools is calculated by the following equation.
where Ri, R′I, R′′i are the distances between any point on the workpiece and a certain differential unit on the heat source. B(x) is the percentage of heat generated at the tool-chip interface entering the chip. 1 − γ is the percentage of heat generated by the heat source qrub entering the tool. Assuming that the temperature at the side of the tool and the temperature at the side of the chip are consistent on the adiabatic boundary, the following relationship exists at the tool-chip interface.
B(x) is calculated by Equation (7) and the temperatures in the primary shear zone and the temperature at the tool-chip interface are obtained by Equations (3) and (4), respectively.

Milling Force Prediction Based on Analytical Method
When the milling force model is established, the cutting-edges of the tool are discretized along the axis of the tool, and the cutting process of each discrete cutting edge is regarded as an oblique cutting to calculate separately, as shown in Figure 2a,b. Milling force is mainly divided into two parts: chip forming force and ploughing force. The chip forming force mainly stems from the normal pressure N and the friction force F exerted by the chip on the tool. The ploughing force mainly stems from the force Pcut in the cutting direction and the force Pthrust perpendicular to the cutting direction caused by the relative movement between the tool tip and the workpiece, and the detailed calculations refer to the study of Su [1]. The cutting thickness and direction of the cutting force are different at different positions on the cutting-edge during milling. Considering the influence of chip side flow on force and temperature is challenging, so it is considered that the chip side flow angle is zero during milling which makes it easier to introduce the iterative algorithm of thermal-mechanical coupling in orthogonal cutting. Figure 2b illustrates the geometric relationship of the oblique cutting process without considering the phenomenon of chips side-flow. Obviously, the angle between the frictional force F and the radial force Fr is the rake angle α. Make a section that is perpendicular to the cutting edge, as shown in Figure 2c, the rotation matrix M1 can be obtained as: Then the radial force Fr and the circumferential force Ft′ perpendicular to the cutting edge are expressed as: v As shown in Figure 2b, the rotation matrix M2 considering the inclination angle β can be obtained as: Then the radial force Fr, the circumferential force Ft, and the axial force Fa are expressed as: Figure 3 shows the trajectory of a point on the cutting edge with an axial height z during milling. The rotation matrix M3 considering the rotation angle of the tool can be obtained as: where ω is the tool speed, t is the cutting time, z is the axial height of the tool, β is the tool helix angle, R is the tool radius, ψ 2 can be expressed as: Then the cutting force dFX along the feed direction, the cutting force dFY in the direction perpendicular to the workpiece surface, and the tool axial cutting force dFZ on the micro cutting edge are expressed as: Assuming that the milling process is down milling, the total cutting time Tc is divided into three stages, as shown in Figure 3. The first stage is the cut-in stage. A part of the cutting edge is in contact with the workpiece, which occurs during the point on the end of the tool moves from A to B. The second stage is the stable cutting stage. The cutting edge is completely in contact with the workpiece, which occurs during the point on the end of the tool moves from B to C. The third stage is the cut-out stage. A part of the cutting edge is in contact with the workpiece, which occurs during the point on the end of the tool moves from C to D. The total milling force is calculated based on three cutting stages, and its expression is: where n is the number of micro cutting edges involved in the cutting process. The cutting thickness hc during milling can be considered as a function of the independent variable t. Figure 4 is a geometric schematic of the cutting thickness, where points O1 and O2 are the positions of the tool axis at time t1 and t2 respectively, and point O is the position of the tool axis when the cutting edge cuts into the workpiece. A reference coordinate system XOY is established with the point O as the origin. The cutting thickness hc is defined as the length of the workpiece passing between the point on the cutting edge and the axis of the tool, and its expression is: x y is the coordinate of the cutting edge,   1 1 , x y is the coordinate of a point on the workpiece surface.
As shown in Figure 4, when the coordinate of a point on the cutting edge is   x y , the cutting thickness hc is the distance between the cutting edge and the surface which is to be machined, and when the coordinate of a point on the cutting edge is   , the cutting thickness hc is the distance between the cutting edge and the machining surface. Coordinate point   , x y   is used as the judgment standard, which can be expressed by the cutting amount and the geometric angle of the tool as: where z f is the feed per tooth, int  is the cut-in angle, which is obtained from the geometric relationship and its expression is: Assuming that the cutting process is ideal without vibration and tool wear, it can be considered that the distance in the horizontal direction between the cutting path of the previous tooth and the cutting path of the current tooth is z f . In the reference coordinate system XOY, the current path can be expressed as: As shown in Figure 4, A and B are on the same straight line. If the slope of the straight line is k, the equation can be obtained: where t1 is the time for the cutting edge to move to    x can be obtained by Equation (21).
Cutting path of the previous tooth Cutting path of the current tooth   The calculation process of the milling force is shown in Figure 5.

Prediction of Residual Stress Considering Thermal-Mechanical Coupling Effect
Su and Liang [1] have developed an analytical modeling of residual stress for dry machining which considers both mechanical effect and thermal effect. In this paper, the iterative algorithm of thermal-mechanical coupling in orthogonal cutting is introduced, and the geometry of the workpiece surface is simplified to calculate the residual stress more effectively, so as to realize the rapid prediction of the peak value of the residual stress.

Mechanical Stress Source and its Stress Distribution Model
When calculating the residual stress, the stress loading history inside the workpiece needs to be known and it can be calculated by Hertzian rolling contact theory. Compared with the primary and tertiary deformation zones, the effect of the secondary deformation zone on the stress distribution of the workpiece is indirect and can be ignored. Therefore, the stress due to mechanical load during cutting is mainly two aspects: one is the stress caused by the shear deformation of the workpiece, the other is the stress caused by the friction and extrusion of the tool tip and the surface of the machined workpiece, as shown in Figure 2c.
The contact between the tool and the workpiece is regarded as Hertzian rolling/sliding contact. At the point (x, z) in the workpiece, the stress component generated by the mechanical load can be obtained by integrating the Boussinesq solution in the loaded area (−a <x < a): Calculated (x 1 ,y 1 ) and (x 2 ,y 2 ) by (19)  Calculated dF X 、dF Y and dF Z by Eq. (14) Calculated (x',y') by Eq. (17) x≥x' Calculated h c by Eq. (16) Calculated (x 1 ,y 1 ) and (x 2 ,y 2 ) by (21) Calculated F X 、F Y and F Z by Eq. (15) Y N F、N、 P cut 、P thrust where 2a is the length of the loading range of the mechanical load, p(s) and q(s) are the normal stress and the tangential stress on the loading area. σx is the stress along the feed direction and σy is the stress along the direction of tool axis.
The normal stress, tangential stress, and loading range on the shear plane are expressed as: The normal stress, tangential stress, and loading range on the friction area between the tool tip and the workpiece are expressed as:

Thermal Stress Source and its Stress Distribution Model
The stress due to thermal load during cutting is mainly two aspects: one is the shear heat source generated on the shear plan, the other is the friction heat source generated by the friction between the tool tip and the surface of the processed workpiece. Thermal stress caused by temperature changes is expressed as: where E is the elastic modulus of the workpiece material, v is the Poisson′s ratio of the workpiece material, αT is the thermal diffusion coefficient of the workpiece material, Gxh, Gxv, Gzh, Gzv, Gxzh and Gxzv represent the Green functions of the plane strain, and the detailed calculations refer to the study of Saif [13]. p(t) is expressed as: The stress at any point on the workpiece due to mechanical and thermal loads is expressed as: Since the milling is discontinuous cutting, the cutting thickness and the cutting direction of the cutting-edge change at all times. Its stress distribution is completely different from that of turning. As shown in Figure 6, plastic deformation occurs at the contact area between the tool and the workpiece during milling. Obviously, the plastic deformation between the A and B points has the largest effect on the residual stress in the machined surface. At the same time, the surface between points A and B is simplified to a plane for calculation. This simplification is reasonable because the surface roughness Rz is much smaller than the feed rate per tooth z f , Furthermore, the calculation of the surface stress distribution of the workpiece is simplified to a plane calculation between A and B points, which reduces the workload of the residual stress calculation and improves the predictive efficiency. Taking the experimental parameters of case 1 in Table 2 as an example, Figure 7 shows the stress contours in the workpiece due to the combined loads.

Hybrid Algorithm for Stress Loading Process
The stress σtotal in the workpiece is discretized into M steps, and the loading calculation is carried out in an incremental manner.
At the end of each loading step, it is necessary to determine whether or not yielding occurs. The TANH constitutive model is used for the yield stress calculation and the subsequent yield surface is updated by combining with the kinematic hardening criterion. Its expression is: When VM < 0, the material belongs to the elastic loading state. When VM = 0, the material belongs to the plastic loading state. The elastic strain is expressed as: where E is the elastic modulus and  is Poisson's ratio and According to the plastic flow criterion, plastic strain is expressed as: where h is a plastic modulus and represents a hardening rate of a material.
Since the stress of workpiece may have a large fluctuation range during the cutting process, the Hybrid function is introduced and expressed as: where  is a constant, G is the shear modulus of elasticity.
xx   and yy   are stress increments in the cutting direction and perpendicular to cutting direction respectively, which can be expressed as: When VM = 0 and dsijnij < 0, calculate the stress increment under plastic release by Equation (38) The residual stress obtained in the previous step is used as the initial stress for the next step until the stress and strain enter a stable state. The calculation result is the residual stress inside the workpiece after the cutting process is completed. The flow chart of calculating the residual stress is shown in Figure  8.

Experimental Machine Tool, Tool and Workpiece Material
The experiments are carried out on a VDL-1000E three-axis CNC machining center, which is made in Profile of General Technology Group Dalian Machine Tool Corporation, Dalian, China. The range of the spindle speed is 45-8000 rpm, and the feed rates in the three directions of X, Y, and Z are 24/24/18 m/min. The tool is an uncoated 2-blade carbide end tool (ANSI: 2P160-1000-NA H10F) produced by Sandvik. The diameter of the tool is 10mm, the helix angle β is 25°, the rake angle α is 13.5°, the tool corner radius is 8 μm, the total length is 90 mm, and the length of cutting part is 32 mm. The workpiece to be machined is a rectangular block of Ti6Al4V alloy with a size of 50 mm × 50 mm × 10 mm, and the accuracy of surface polish is Ra 0.4 μm.
The experimental site is shown in Figure 9. The force signal during milling is measured by kistler5236b (Kistler Group, Winterthur, The Switzerland), in which force sensor is composed of rotor and stator. The rotor is the tool handle, which is used to clamp the tool. The stator is the signal receiving end, which receives the signal measured by the rotor. The uncertainty of the measurement is in the range of 1% of the measured value. The workpiece is fixed on the machine tool guideway with a fixture, and the machining method is down milling without cutting fluid. To avoid cutting vibrations and obtain high surface quality, the experimental scheme is shown in Table 2.
Calculate the stress and strain  The material parameters needed to be input into the prediction model are divided into two parts: tool material parameters and workpiece material parameters. The material parameters of workpiece (Ti6Al4V) and tool (WC) refer to the research of the literature [14], as shown in Table 3. Table 3. Material parameters of workpiece and tool [14].

Experimental Validation of Milling Force Model
The comparison results between the experimental data and the predicted values of the milling force under the four sets of experimental parameters shown in Table 2 are shown in Figure 10. By the average value r   of the absolute value of the error between the experimental data and the predicted value, the prediction accuracy of the model is quantitatively evaluated. r   is calculated by the following equation:  Table 4. The cutting force prediction value in X direction is more accurate.  According to Figure 10, it can be found that the prediction error of the milling force FY is large, and r   can reach 52.75 N. When entering the stable cutting stage, FY suddenly presents a fluctuating state. This is due to the poor rigidity of the tool (or workpiece), resulting in the cutter-drawing, which further caused the vibration between the tool and the workpiece. Compared with the experimental data of milling forces in other two directions, it can be found that there is still a wave value of convergence in the milling force FY after the tool cuts out from the workpiece, which proves that there is cutting vibration in the direction perpendicular to the workpiece.
The increase of cutting speed increases the milling force FX and FY, but it has little effect on the milling force FZ. Therefore, the experimental value of the cutting force FZ in case 2 is less than the predicted value. The increase of cutting width causes greater cutting vibration, and the peak value of milling force FY in case 3 increases significantly, which shows that the cutting vibration is enhanced.
The milling force model divides the milling process into three cutting stages, and the experimental data show that the milling force fluctuates obviously in the stable cutting stage. Figure  11 illustrates the contact relationship between the cutting edge and the workpiece in the stable cutting stage. t1 and t2 are different time points in the same stable cutting stage, corresponding to two different cutting states. Considering the cutter-drawing due to excessive tool overhang during milling, it is considered that in the stable cutting stage, the contact area between the cutting edge and the workpiece at time t2 is higher in the axial direction than that at time t1. The position of the milling force acting on the tool changes, which changes the degree of deflection. As a result, the milling force fluctuates during the stable cutting stage. In addition, it may also be due to friction chatter caused by friction between the tool and the workpiece in the same direction as the cutting speed [15].

Measurement Device and Measurement Scheme for Residual Stress
The instrument for measuring the residual stress is an X 'Pert-Pro X-ray diffractometer and the test is performed using method 2 sin  , as shown in Figure 12. For the residual stress test of Ti6Al4V alloy, Cu target was selected as the X-ray tube. The XRD stress test parameters for Ti6Al4V are given in Table 5.  Based on the theory of residual stress prediction algorithm, it is considered that under ideal conditions, the residual stress distribution at any point on the surface of the workpiece is consistent. In order to facilitate the measurement of the residual stress and the electrolytic etching in the machined surface, the position with a depth of 4 mm along the blade axis was selected for measurement. According to the milling force data shown in Figure 10, the stage where the milling force is relatively stable is selected for measurement. When electrolytic etching the surface, the selected electrolyte is saturated NaCl solution, the voltage is set to 15 V, the diameter of the nozzle is 15 mm, and the corrosion depth is measured by a SKCH-1A thickness gauge after each electrolytic corrosion.

Experimental Validation of Residual Stress Model
The predicted and experimental values of the residual stress under the experimental parameters of case 1 and case 3 are shown in Figure 13, it can be seen that the influence depth of residual stress is about 50 μm and the depth of the peak value of subsurface residual stress is between 15 μm and 20 μm. The variation of residual stress σx and σy with depth is basically the same. The predicted value converges faster than the experimental value, which may be caused by the simplified method proposed in Figure 6. Because the model only considers the stress distribution between AB two points, the response depth of the predicted residual stress is smaller than the experimental residual stress. When the depth into workpiece is 0 μm, the predicted value of residual stress is different from the experimental value. The predicted value is tensile stress, but the experimental value is compressive stress. In the surface residual stress measurement of most machined workpieces, surface oxidation may be an important cause of prediction error. Moreover, the prediction model of residual stress in this paper does not consider the effect of phase transformation. For Ti-6Al-4V materials, the phase transformation from α to β will lead to the increase of grain volume, which will lead to the greater compressive residual stress on the surface. Similarly, it also explains that the residual stress of the surface tends to be tensile stress.
The peak value of residual stress in the machined surface has the greatest effect on the fatigue life of the workpiece. It is particularly important to accurately predict this feature. Figure 14 shows the predicted and experimental values of the residual stress under the four sets of process parameters. The residual stress prediction algorithm is based on the plane stress assumption, so the stress σz perpendicular to the surface of the workpiece is not considered. The numbers in Figure 14 indicate prediction errors. As shown in Figure 14, the residual stress prediction model can accurately predict the peak value of the residual stress in the subsurface. The predicted error of case 3 and case 4 are relatively large. Comparing the milling force FY in Figure 10, it can be found that in case 3 and case 4, the cutting width and the feed rate are increased, which enhances the cutting vibration during the stable cutting stage, resulting in changes in the stress distribution in the machined surface. This may be one of the reasons for the prediction error of residual stress.

Conclusions
In this paper, the milling force and residual stress model of titanium alloy are established with considering the coupling effect of thermal-mechanical load. In this model, the complex contact relationship between tool and workpiece and the time-varying cutting thickness model are considered. The conclusions are as follows: (1) The contact relationship between the tool and the workpiece during the milling is analyzed in detail. Compared with the milling forces FX and FZ, it can be found that there is still a wave value of convergence in the milling force FY after the tool cuts out from the workpiece, which proves that there is cutting vibration in the direction perpendicular to the workpiece. It is also found that the increase of radial cutting depth and feed rate will significantly increase the tool vibration. In the direction perpendicular to the workpiece surface, the milling force in the stable cutting stage fluctuates due to the influence of the cutter-drawing, but the milling forces FX and FZ can still be predicted well. (2) In order to use the thermal coupling iterative algorithm and realize the rapid prediction of the peak value of the residual stress, the geometry of the workpiece surface is simplified. The residual stress model accurately predicted the peak value of the residual stress in the subsurface of the workpiece. With the increase of the radial cutting depth and feed rate, the cutting vibration will be intensified, the stress state in the surface of the workpiece will change, and the prediction error of the residual stress will increase. (3) According to the prediction results of residual stress, it can be seen that the influence depth of residual stress is about 50 μm and the depth of the peak value of subsurface residual stress is between 15 μm and 20 μm. The variation of residual stress σx and σy with depth is basically the same. The stress in the machined surface is tensile stress, and the peak value of subsurface stress is compressive stress, which may be due to the influence of cutting heat, resulting in a large tensile stress in the surface.

Conflicts of Interest:
The authors declare no conflict of interest.