Numerical prediction and experimental investigation of residual stresses in sequential milling of GH4169 considering initial stress effect

Residual stresses affect the service life and cause the deformation of machined parts. Therefore, the study on the development of residual stresses during machining operation is very important. For this, analytical and numerical models have been developed to predict the residual stresses for single-step machining operation. However, the parts are usually machined with a serial of production processes. A numerical model for predicting residual stress in sequential side milling GH4169, considering initial stress, was proposed in this research. The initial residual stress distribution due to previous-step milling was considered in the proposed model. It was assumed that this initial residual stress value changed with the depth from the machined surface, while the residual stress on the identical horizontal plane was assumed with uniform distribution. The proposed numerical simulation model could predict the stress value of the machined surface, the stress value of the compression valley, and the residual stress distribution in the depth direction beneath the machined surface with high accuracy. Experimental investigations of sequential side milling GH4169 were conducted and the generated residual stresses were measured. The distribution trend and influence rule of residual stress between two sequential milling steps were analyzed. The residual stress distribution beneath the surface showed a spoon-shaped pattern. The thickness of the residual stress influence layer (RSIL) after the current milling step was affected by the depth of cut and the RSIL thickness of the previous milling step. The numerical model could predict the thickness of RSIL and optimize the depth of cut in sequential side milling.


Introduction
Superalloy GH4169 has been widely used in the aerospace field due to its excellent high-temperature strength, heat corrosion resistance, high-temperature oxidation resistance, and fatigue resistance [1,2]. High surface quality is required for aerospace components manufactured by side milling GH4169. Residual stress is one of the important indicators that justifies the quality of machined surface [3,4]. The residual stress affects the service life and causes the deformation of machined parts [5,6]. There are three sources for the formation of residual stress including: nonuniform plastic deformation, volume change caused by phase transformation/precipitation, and chemical reaction [3,4]. During in metal cutting process, the non-uniform plastic deformation is the main factor for residual stress formation. The residual stresses due to plastic deformation at macroscale are caused by non-uniform mechanical and thermal loads. The residual stress is firstly caused by plastic deformation at micro-scale in grains such as slip, transgranular slip, twin and dislocation. The residual stress can also be caused by the plastic deformation of the grain when it is heated. This plastic deformation is mainly due to the anisotropy of thermal expansion coefficient and elastic coefficient of crystal, as well as the different orientation between grains. The workpiece GH4169 is polycrystalline material. The intergranular stresses exist in polycrystalline materials. It is necessary to study the residual stress in cutting GH4169 both at macro and micro scales.
In the past few decades, the researches on residual stress caused by machining mainly focused on the influence of tool parameters and machining parameters in single-step machining. However, the parts are usually finished with multi-step machining processes. Several researchers have studied the influence of residual stresses between serial steps.
Karabelchtchikova et al. [7] carried out an experimental research of the residual stress in multi-pass grinding. This research mainly studied the residual stress distribution in heat treatment and subsequent multi-pass grinding. The results showed a non-linear superposition relationship in the residual stress distribution between each process step.
Dehmani et al. [8] developed a multi-step finite element model for orthogonal cutting. The influence of cumulated strain and heat induced by the different steps on the prediction of residual stress distribution was simulated. The distribution of the residual stress in depth simulated with the proposed model was consistent with that by experiment measured. The results showed that the multi-step model could provide better results for the prediction of residual stress than a single-step model did. However, the initial stress effect was neglected with the developed model.
Grochała et al. [9] conducted a study on the residual stress of the burnished surface. A burnishing FEM model of the surface with the milling track was developed. It was found that the stress value after burnishing was significantly affected by the change of surface roughness. The distribution depth of the maximum compressive residual stress after burnishing decreased with the increase of the initial roughness of the milled surface.
[10] proposed an analytical model for the prediction of residual stress in the multi-step milling. The authors claimed that the proposed analytical model could predict the residual stress value after each step. However, the residual stress prediction accuracy of this model was affected by the over-simplification of the model assumptions.
Ma et al. [11] conducted an experimental research to study the evolution of the residual stress field during multistep face milling. It was found that the surface residual stress was mainly determined by the final processing step. The residual stress distribution in depth depended on the superposition of the processing history.
Singla et al. [12] conducted a numerical study of the residual stress in the two-step turning process. A numerical simulation model of two-step turning was developed. The surface after two-step turning could be obtained by the design of a double-edged turning tool in sole one-step cutting. The initial residual stress effect was neglected in the numerical simulation model of two-step turning.
The residual stresses could accumulate and redistribute between serial steps or processes. However, the effects of the initial stress generated in the previous machining step on the residual stress in current machining step are seldom considered in these developed models, which makes it hard to attain satisfied prediction accuracy of the final residual stress distribution in the machined surface. A numerical model for predicting machined surface residual stress in sequential side milling GH4169 considering initial stress will be proposed in this research. The initial residual stress distribution due to previous-step milling is considered in the proposed model. The experimental investigations of sequential side milling GH4169 are conducted to verify the model prediction results.

Conversion equation from milling to equivalent orthogonal cutting
Effective conversion from three-dimensional (3D) side milling to two-dimensional (2D) orthogonal turning makes it easier to import the initial residual stress into the simulation model [13]. Figure 1 shows the process of converting 3D side milling to 2D orthogonal turning. Due to the geometry of the milling tool, the chips generated during milling cannot remain in the same thickness. Therefore, the equivalent depth of cut changes with the rotation angle in the orthogonal model. The average depth of cut is defined as, where t c is the average depth of cut, V f is the feed rate, and n is the spindle speed, the instant equivalent cutting depth is, where = 2 × n × t is the rotation angle, t is the cutting time. The equivalent side cutting edge angle C ′ r is, where C r is the side cutting edge angle, and λ is the chip flow angle. The equivalent chip flow angle ′ is given by the where s is the inclination angle, α is the rake angle, and 0 is defined as, The equivalent rake angle is, The equivalent cutting depth is, The equivalent cutting speed is, where V r = 2 R × n is the rotation speed and R is the tool radius.

Modification equation of measured residual stress in depth
The workpiece is stripped due to the measurement of residual stress in depth. The surface stress will be redistributed after the delamination. Therefore, the residual stress measured at a point on the surface is not the value at the same depth before the delamination. The measured value needs to be corrected to obtain the accurate residual stress distribution in the depth direction of the workpiece. Since the thickness of the milling layer is much smaller than the length and width of the workpiece, the blocky workpiece in this paper can approximate a flat plate. Several studies have shown that the residual stress changed with the distance from the machined surface. It was assumed that the residual stress value changed with the depth from the machined surface, while the residual stress on the identical horizontal plane was assumed with uniform distribution [14][15][16]. The residual stress is assumed to distribute symmetrically along the center plane of the plate. At the same time, the residual stress is parallel to the center plane of the plate. As shown in Fig. 2, if the length and width of the plate are much greater than the thickness, the plate can be regarded as an infinite surface when the residual stress in the center of the plate is measured. The plate is peeled from the bottom surface at point 0. It is assumed that the accurate stress is xx (z) and yy (z) when the layer is peeled to the depth z. The measured stress is � xx (z) and � yy (z) . The stress release after peeling Δ xx (z) and Δ yy (z) can be calculated in Eqs. (9) and (10).
where h is the thickness of the plate, z is the reduced thickness of the plate after polishing is completed, ζ is the reduced thickness of the plate during the removal process. Since the peeling is performed on the one side, the bending moments M x and M y are added to the plate. M x and M y can be calculated in Eqs. (11) and (12). The bending moments  cause additional stresses Δ Mx (z) and Δ My (z) . Δ Mx (z) and Δ My (z) can be calculated in Eqs. (13) and (14).
Moreover, By substituting Eqs. (9) to (14) into Eqs. (15) and (16), and replacing ′ with σ for the convenience of calculation, xx (z) and yy (z) can be calculated, In this paper, xx (z) and yy (z) in Eqs. (17) and (18) are taken as the residual stress values in the depth direction of the workpiece to ensure the accuracy of the experimental data.
3 Finite element simulation of sequential side milling GH4169

Geometry modeling and mesh controlling
The metal cutting finite element simulation software AdvantEdge FEM 5.6 was used to simulate the twodimensional orthogonal cutting process. The workpiece was a rectangle with a length of 5 mm and a height of 2 mm. The tool length and width were both 1.5 mm. The cutting edge radius of the tool was 0.04 mm, the rake angle was 13.8°, and the clearance angle was 10°.
In the simulation, the workpiece material was nickelbase superalloy GH4169. The Johnson-Cook constitutive model (Eq. 19) was used in the simulation.
where A, B, C, m, n are the material constants, A is the yield stress, B is the hardening modulus, C is the strain rate dependency coefficient, m is the thermal softening coefficient, and n is the strain hardening coefficient. is equivalent strain, ̇ is equivalent strain rate, and ̇ is the reference strain rate. T is the overall temperature of the workpiece, T m is the melting temperature of GH4169, and T r is room temperature. Table 1 shows the Johnson-Cook parameters of GH4169. The physical and mechanical properties of the tool and workpiece material are specified in Table 2. The tool material was  The mesh of the simulation model was divided by the Lagrange method. The finite element mesh deforms with the deformation of the workpiece. The adaptive mesh technology in AdvantEdge was fully utilized to improve the calculation accuracy. The mesh of the workpiece was constructed by triangular elements, and the mesh elements were adaptively divided. The closer to the contact surface between the tool and the workpiece, the smaller the element was. The maximum element size of the workpiece was 0.1 mm, the minimum element size of the workpiece was 0.02 mm, and the workpiece had a total of 30,000 nodes. The maximum tool element size was 0.1 mm, and the minimum tool element size was 0.06 mm.
Chip separation criteria and adaptive mesh refinement were applied. In order to improve the prediction accuracy of the residual stress, the simulation refined the mesh of RSIL and performed the thermal and mechanical coupling calculations. Figure 3 shows the boundary constraints of the simulation model. Simulation constraints were imposed on the left and bottom edges of the workpiece in the horizontal and vertical directions. The initial temperature was the environment temperature of 25 °C. The friction between the workpiece and the rake face followed the Coulomb friction law. The friction coefficient was 0.3 [17]. The effect  of the initial residual stress distribution was considered in this model. It was assumed that this initial residual stress value changed with the depth from the machined surface, while the residual stress on the identical horizontal plane was assumed with a uniform distribution. The residual stress distribution of the workpiece in depth (0-100 μm) before milling was measured. Then, the measured initial stress was imported into the simulation model. Figure 4  shows the initial residual stress distribution of the workpiece in the simulation model. Subroutines and secondary cutting functions in AdvantEdge were used to implement the sequential milling process. Figure 5 depicts the temperature and residual stress distribution in the simulation model of sequential side milling.

Sequential milling experiments and measurements of residual stress
The material used in the experiment was GH4169. The metallographic structure of GH4169 is mainly composed of γ matrix, dispersed-distributed strengthening phases γ′ and γ′′, δ phase, NbC, TiN, etc. The structure and chemical composition of its main constituent phases are shown in Table 3.
The size of workpiece was 15 mm × 10 mm × 20 mm. The workpiece was pretreated by electrical discharge machining (EDM) and heat treatment to eliminate part of the stress. As shown in Fig. 6a, the sequential side milling experiment was carried out on a vertical machining center (Shenyang Machine Tool VMC0540d), and the three-way dynamic cutting dynamometer (Kistler 9257B) was used to record the milling force of each side milling step. Milling tools VSM-4E-D6.0 were cemented carbides with four-flute tool and TiN coating. During the experiments, the dry cutting was adopted in this study. The milling method was down milling due to the consideration of milling surface quality. Since the influence of tool wear on the residual stress was not considered in this paper, a new tool was used for each milling.
The most influential factor in milling was the depth of cut. Therefore, only the radial depth of cut changed in sequential side milling. The sequential side milling experiment was divided into two steps in total. A total of four sets of sequential side milling GH4169 experiments were carried out. The experimental parameters of the first milling step were the same, radial depth of cut a e = 0.3 mm, feed per tooth f = 0.0075 mm/z, speed v = 38 m/min, axial depth of cut a p = 8 mm. The experiment kept the other parameters of the second milling step unchanged and only changed the radial depth of cut. The radial depth of cut a e was 0.02, 0.05, 0.1, 0.3 mm, respectively.
The workpiece material GH4169 is a polycrystalline material, X-ray diffraction method can be used to determine the residual stress. The X-ray diffraction phenomenon is shown in Fig. 7. X-ray diffraction satisfies the Bragg Eq.
where θ is the angle between the incident X-ray and the crystal plane, d is the distance between the crystal planes, λ is the X-ray wavelength, and n is the diffraction index. (20) 2d sin = n As shown in Fig. 7 when X-rays are diffracted, the angle between incident X-rays and diffracted X-rays is 2θ, which is called the diffraction angle. According to the Bragg equation, the interplanar crystal spacing d can be obtained by measuring the diffraction angle. Therefore, the sin 2 φ method can be used to obtain the stress at any point and direction on the grain boundary plane (hkl) in Eq. (21), where E is the elastic modulus, ν is the Poisson's ratio, φ is the angle between the measured strain direction and the normal direction of the crystal plane, and θ 0 is the diffraction half angle when φ = 0, which is a stress-free state.
As shown in Fig. 6c, the X-ray diffraction stress analyzer (Pulstec μ-X360n) was used to measure the residual stress on the machined surface. The residual stress measurement position was selected on three points (1, 2, 3) at the center of the machined surface as shown in Fig. 6d. Residual stresses in the feed direction (σ 11 ) and perpendicular to the feed direction (σ 22 ) were measured, respectively. In order to ensure the accuracy of the data, the average measured value of the three points was taken as the residual stress on the machined surface of each workpiece.
The main purpose of studying the residual stress was to improve the fatigue life of the workpiece. It was more important to study the distribution in depth of residual stress than to study that on the machining surface only. Due to the limitation of X-ray penetration depth (10 μm), it was necessary to peel off the workpiece layer by layer before using X-ray diffraction for the measurement of residual stress. This paper used an electrolytic corrosion polishing process to peel off the workpiece surface. Residual stresses were measured at 8 positions in depth (10 μm, 20 μm, 30 μm, 40 μm, 50 μm, 60 μm, 80 μm, 100 μm). The electrolytic corrosion polishing equipment (Shanghai Optical Instrument DFP-2) was shown  Fig. 6b. The electrolytic polishing liquid was CH 3 OH and HN0 3 (CH 3 OH: HN0 3 = 2:1). The specific parameters of the electro polishing were: voltage 20 V, current 5 A. Polishing temperature was the environment temperature (25 °C).

Residual stress distribution of experiment
The milling sequence and the residual stress distribution in the feed direction beneath the surface are given in Fig. 8. Figure 9 shows that the distribution of residual stress beneath the machined surface presented a spoon-shaped pattern. The milling residual stress distribution trend was consistent with that of Quan's research results [20]. The residual tensile stress on the machined surface was less in the feed direction than perpendicular to the feed direction. Table 4 shows that the residual stress of the compression valley was less in the feed direction than perpendicular to the feed direction. The thickness of RSIL was between 80 and 100 μm before milling. The thickness of RSIL was also between 80 and 100 μm in the first step milling (a e = 0.3 mm).
It was observed from Fig. 9a and c that the compressive residual stress value of the compression valley became smaller in the second-step milling (a e = 0.3 + 0.02 mm and a e = 0.3 + 0.05 mm), The compression valley shifted to the machined surface (area 2), and the residual stress on the machined surface moved in the direction of tensile stress (area 1). As shown in Fig. 9b, when the depth of cut was a e = 0.3 + 0.02 mm in the second-step milling, the residual stress distribution curve was shifted to the right by 0.02 Fig. 8 Milling sequence and residual stress distribution in feed direction beneath machined surface Fig. 9 Residual stress distribution beneath the machined surface of the sequential side milling GH4169 mm. It could be seen from the curve slope trend that the residual stress distribution was only significantly affected within the range of 0.02 to 0.045 mm. Area 1 to 3 in Fig. 9a presented that RSIL became smaller. The thickness of RSIL was between 60 and 80 μm. As shown in Fig. 10, the reason for this phenomenon was that the thickness of RSIL in the current step milling was affected by the depth of cut and the thickness of the previous step. When the depth of cut b in the second-step milling was less than the thickness of RSIL, the contact area between the tool and the workpiece was reduced, resulting in a smaller cutting force. The cutting force became smaller, so the plastic deformation on the workpiece surface was reduced. The reduction of mechanical stress effect led to the decrease of RSIL. As shown in Fig. 9d, when the depth of cut in the second-step milling was a e = 0.3 + 0.05 mm, the residual stress distribution curve in depth was shifted to the right by 0.05 mm. It could be seen from the curve slope trend that the residual stress distribution was significantly affected within the range of 0.05 to 0.1 mm.
As shown in Fig. 9i and k, when the depths of cut were a e = 0.3 + 0.1 mm and a e = 0.3 + 0.3 mm in the second-step milling, the compressive residual stress value of the compression valley and the tensile stress value of the machined surface were increased. The increase of the stress value was small. The depths of the compression valley from the machined surface were about 30 μm. The thickness of RSIL in the first step was the same as in the second step. The thickness of RSIL was still between 80 and 100 μm. The reason for this phenomenon was that when the depth of cut was greater than the thickness of RSIL, the residual stress distribution in the previous step had almost no effect on stress distribution in the current step. The small increase in the stress value of the machined surface was caused by the accumulation of the thermal loads in the side milling.
Therefore, the machining history was needed to be considered when studying the final residual stress distribution. When the depth of cut was greater than the thickness of RSIL, the residual stress distribution of the previous step had little effect on the current step. On the contrary, the residual stress distribution of the previous step had obvious influence on the current step. Compared with Ma's paper [11], the conclusion of this paper is more comprehensive.

Residual stress distribution of simulation
The simulation residual stress distribution curve (shown in Fig. 11) was drawn according to the drawing method of the experimental residual stress distribution. The depth of the compression valley in the first and second-step milling was at 22 μm from the machined surface. When the depth of cut was a e = 0.3 + 0.1 mm, the residual stress of the compression valley was reduced by 52 MPa compared with a e = 0.3 mm. The tensile stress of the machined surface was increased by 43 MPa. The residual stress distribution curve of a e = 0.3 + 0.3 mm and a e = 0.3 mm almost coincided. This simulation result could verify the conclusion drawn by the experimental part because the RSIL was smaller than the depth of cut, and there was almost no interaction between the residual stress distributions of the two steps.

Comparison of residual stress between simulation and experiment
A comparison of residual stress distributions of the experiment with simulation is given in Fig. 12. When the depth of cut was a e = 0.3 mm in the first step milling, the depth of residual stress compression valley in the simulation model was at 22 μm from the machined surface. The depth of residual stress compression valley in the experiment was at 30 μm from the machined surface. The residual stress values of machined surface and compressed valley were in good agreement between experiment and simulation. When the depths of cut were a e = 0.3 + 0.02 mm and a e = 0.3 + 0.05 mm in the second-step milling, the depths of residual stress compression valley in the simulation and experiment were both at 22 μm from the machined surface. The residual stress values of the machined surface and compressed valley were in good agreement between experiment and simulation. The difference between experimental and simulation values was generally less than 60 MPa. The distribution of residual stress in depth was highly consistent between the experiment and simulation.
When the depth of cut was a e = 0.3 + 0.1 mm and a e = 0.3 + 0.3 mm in the second-step milling, the depth of residual stress compression valley in the simulation model was 22 μm from the machined surface. The depth of the residual stress compression valley in the experiment was at 30 μm from the machined surface. The residual stress values of the machined surface and compressed valley were in good agreement between experiment and simulation values. There was a difference between the prediction of the residual stress values at 20-40 μm from the machined surface between the experiment and the simulation. The difference was mainly caused by the location of the simulated compression valley. If the residual stress distribution curve of the simulated model was shifted 8 μm to the right, it could be found that the curve also had a strong consistency between the experiment and the simulation. Compared with the experimental values, the average deviation of the residual tensile stress on the machined surface in the sequential simulation model was 7.2%. The prediction accuracy of the model was 92.8%. The prediction accuracy in the proposed model was higher than in Dehmani's model [8]. In summary, the proposed model provided high prediction accuracy of the residual stress in sequential milling.

Conclusion
A numerical model for predicting the machined surface residual stress in sequential side milling GH4169 considering initial stress was proposed in this research. The experimental investigations of sequential side milling GH4169 were conducted for the purpose of verifying the model prediction results. According to the experiment and simulation results, the following conclusions can be drawn: 1. The surface residual tensile stress and compressive valley stress in the feed direction were less than perpendicular to the feed direction. The distribution of the residual stress in the depth direction from the surface presented a spoon-shaped pattern. When the depth of cut in the second-step milling was smaller than the thickness of RSIL, the compressive stress value of the compression valley became smaller. The compression valley shifted to the machined surface, and the machined surface residual stress moved in the direction of tensile stress. 2. The thickness of RSIL in the current step was affected by the depth of cut of the current step and the thickness of RSIL in the previous step. When the depth of cut was greater than the thickness of RSIL, the residual stress distribution of the previous step had little effect on the current step. On the contrary, when the depth of cut was less than the thickness of RSIL in the previous step, the contact area between the tool and the workpiece was reduced, resulting in a smaller cutting force. The cutting force became smaller, so the plastic deformation on the workpiece surface was reduced. The reduction of mechanical stress effect led to the decrease of RSIL. 3. The proposed simulation model was highly consistent with the experimental conclusions. This model could make a good prediction of the stress value of the machined surface and the compression valley, and the residual stress distribution in depth beneath the machined surface of the sequential side milling. The prediction accuracy of this model was 92.8%. This model could predict the thickness of RSIL and optimize the depth of cut of the design sequential side milling.
Author contribution Gonghou Yao: Methodology, investigation, formal analysis and writing draft. Zhanqiang Liu: Project administration, formal analysis, review and editing. Qinghua Song: Supervision, review and editing. Bing Wang: Supervision. Yukui Cai: Review.
Funding The authors would like to acknowledge the financial support from the National Key Research and Development Program of China (2019YFB2005401). This work was also supported by grants from the National Natural Science Foundation of China (No. 91860207) and Taishan Scholar Foundation.
Availability of data and materials All data and materials used or analyzed during the current study are included in this manuscript.

Declarations
Ethics approval and consent to participate The manuscript is approved by all authors for publication. This manuscript is original and the entire paper or any part of its content has never been published in other journals. The data and results are true and clear. We guarantee the transparency and objectivity of research, and strictly abide by the ethical and professional standards of the industry. All the authors listed agree to participate in this manuscript.

Consent for publication
All co-authors agree to publish the version of this work in The International Journal of Advanced Manufacturing Technology.

Competing interests
The authors declare no competing interests.