Analytical Prediction of Residual Stress in the Machined Surface during Milling

Abstract

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.

1. 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 thermal-mechanical 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.

2. 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 V_S 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 P_total can be expressed by the flow stress σ_AB on the shear plane AB by the geometric relationship as:

$$P_{total}=\frac{{\sigma }_{AB}{h}_{c}w}{\sin \mathrm{\Phi }\cos \theta }$$

(1)

where h_c 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 P_total and the shear plane AB. The friction F and normal force N at the rake face of the tool, the force F_c in the cutting direction and the force F_t in the feed direction can be expressed as:

$$\{\begin{array}{l}F=P_{total}\sin \lambda \\ N=P_{total}\cos \lambda \\ F_c=P_{total}\cos \left(\lambda -\alpha \right)\\ F_t=P_{total}\sin \left(\lambda -\alpha \right)\end{array}$$

(2)

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. Parameters of the TANH constitutive model [11].

A (MPa)	B (MPa)	n	C	m	${\dot{\mathit{\epsilon }}}_0$	a	b	c	d	Tm (°C)
968	380	0.421	0.0197	0.577	0.1	1.6	0.4	6	1	1650

.

Since many parameters in the Oxley’s cutting force model need to be determined by experiments, Ji Xia [8] proposed an iterative algorithm to capture the parameters that need to be experimentally defined in the Oxley′s cutting force model. This algorithm makes the input parameters of the cutting 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.

$$\Delta T_{c-s}\left(X^{ch},Z^{ch}\right)=\frac{q_{shear}}{2\pi K_c}{\displaystyle \underset{l_i=0}{\overset{L_{AB}}{\int }}{e}^{\frac{-(X^{ch}-x_i)V_c}{2a_c}}}\left\{K_0\left(\frac{R_i{V}_c}{2a_c}\right)+\frac{1}{2}{K}_0\left(\frac{{R^{\prime }}_i{V}_c}{2a_c}\right)+\frac{1}{2}{K}_0\left(\frac{{R^{″}}_i\hspace{0.17em}{V}_c}{2a_c}\right)\right\}d{l}_i$$

(3)

$$\Delta T_{c-f}(X^{ch},Z^{ch})=\frac{q_f}{\pi K_c}{\displaystyle \underset{l_i=0}{\overset{h}{\int }}B\left(x\right)e^{\frac{-(X^{ch}-x)V_c}{2a_c}}}\left\{K_0\left(\frac{R_i{V}_c}{2a_c}\right)+K_0\left(\frac{{R^{\prime }}_i{V}_c}{2a_c}\right)+K_0\left(\frac{{R^{″}}_i\hspace{0.17em}{V}_c}{2a_c}\right)\right\}dx$$

(4)

$$\Delta T_{t-f}(X^t,Y^t,Z^t)=\frac{q_f}{2\pi K_t}{\displaystyle \underset{0}{\overset{h}{\int }}\left(1-B\left(x\right)\right)}{\displaystyle \underset{-w/2}{\overset{w/2}{\int }}\left(\frac{1}{R_i}+\frac{1}{{R^{\prime }}_i}\right)}dydx$$

(5)

$$\Delta T_{t-rub}(X^t,Y^t,Z^t)=\frac{q_{rub}}{2\pi K_t}{\displaystyle \underset{0}{\overset{CA}{\int }}\left(1-\gamma \right)}{\displaystyle \underset{-w/2}{\overset{w/2}{\int }}\left(\frac{1}{R_i}+\frac{1}{{R^{\prime }}_i}\right)}dydz$$

(6)

where R_i, 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 q_rub 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.

$$\Delta T_{c-s}+\Delta T_{c-f}+T_0=\Delta T_{t-f}+\Delta T_{t-rub}+T_0$$

(7)

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.

3. 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 P_cut in the cutting direction and the force P_thrust 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 F_r is the rake angle α. Make a section that is perpendicular to the cutting edge, as shown in Figure 2c, the rotation matrix M₁ can be obtained as:

$$M_1=\left[\begin{array}{ccc}\cos \alpha & -\sin \alpha & 0\\ \sin \alpha & \cos \alpha & 0\\ 0& 0& 1\end{array}\right]$$

(8)

Then the radial force F_r and the circumferential force F_t′ perpendicular to the cutting edge are expressed as:

$$\left[\begin{array}{c}{F}_r\\ {F^{\prime }}_t\\ 1\end{array}\right]=M_1\left[\begin{array}{c}F\\ N\\ 1\end{array}\right]+\left[\begin{array}{c}{P}_{thrust}\\ P_{cut}\\ 1\end{array}\right]$$

(9)

As shown in Figure 2b, the rotation matrix M₂ considering the inclination angle β can be obtained as:

$$M_2=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \beta & 0\\ 0& \sin \beta & 0\end{array}\right]$$

(10)

Then the radial force F_r, the circumferential force F_t, and the axial force F_a are expressed as:

$$\left[\begin{array}{c}{F}_r\\ F_t\\ F_a\end{array}\right]=M_2\left[\begin{array}{c}{F}_r\\ {F^{\prime }}_t\\ 1\end{array}\right]$$

(11)

Figure 3 shows the trajectory of a point on the cutting edge with an axial height z during milling. The rotation matrix M₃ considering the rotation angle of the tool can be obtained as:

$$M_3=\left[\begin{array}{ccc}\sin \left({\psi }_2-\omega t+z\tan \beta /R\right)& \cos \left({\psi }_2-\omega t+z\tan \beta /R\right)& 0\\ -\cos \left({\psi }_2-\omega t+z\tan \beta /R\right)& \sin \left({\psi }_2-\omega t+z\tan \beta /R\right)& 0\\ 0& 0& 1\end{array}\right]$$

(12)

where $\omega$ is the tool speed, t is the cutting time, z is the axial height of the tool, $\beta$ is the tool helix angle, R is the tool radius, ${\psi }_2$ can be expressed as:

$${\psi }_2=\frac{\pi }{2}-\arcsin \left(\frac{R-a_e}{R}\right)$$

(13)

Then the cutting force dF_X along the feed direction, the cutting force dF_Y in the direction perpendicular to the workpiece surface, and the tool axial cutting force dF_Z on the micro cutting edge are expressed as:

$$\left[\begin{array}{c}d{F}_X\left(z,t\right)\\ d{F}_Y\left(z,t\right)\\ d{F}_Z\left(z,t\right)\end{array}\right]=M_3\left[\begin{array}{c}{F}_r\\ F_t\\ F_a\end{array}\right]=M_3{M}_2\left(M_1\left[\begin{array}{c}F\\ N\\ 1\end{array}\right]+\left[\begin{array}{c}{P}_{thrust}\\ P_{cut}\\ 1\end{array}\right]\right)$$

(14)

Assuming that the milling process is down milling, the total cutting time T_c 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:

$$F_j^{\left(n\right)}={\displaystyle \sum _{i=1}^{n}d{F}_j\left(i\Delta z,t\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=x,y,z$$

(15)

where n is the number of micro cutting edges involved in the cutting process.

The cutting thickness h_c 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 O₁ and O₂ are the positions of the tool axis at time t₁ and t₂ 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 h_c 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:

$$h_c=\sqrt{{\left(y_2-y_1\right)}^2+{\left(x_2-x_1\right)}^2}$$

(16)

where $\left(x_2,y_2\right)$ is the coordinate of the cutting edge, $\left(x_1,y_1\right)$ 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 $\left(x_2^{(1)},y_2^{(1)}\right)$, the cutting thickness h_c 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 $\left(x_2^{(2)},y_2^{(2)}\right)$, the cutting thickness h_c is the distance between the cutting edge and the machining surface. Coordinate point $\left(x^{\prime },y^{\prime }\right)$ is used as the judgment standard, which can be expressed by the cutting amount and the geometric angle of the tool as:

$$\{\begin{array}{l}{x}^{\prime }=f\cdot t+R\cos {\phi }_{\mathrm{int}}-f_z\\ y^{\prime }=R-a_e\end{array}$$

(17)

where $f_z$ is the feed per tooth, ${\phi }_{\mathrm{int}}$ is the cut-in angle, which is obtained from the geometric relationship and its expression is: ${\phi }_{\mathrm{int}}=\arcsin \left(\left(R-a_e\right)/R\right)$.

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 $f_z$. In the reference coordinate system XOY, the current path can be expressed as:

$$\{\begin{array}{l}{x}_2=f\cdot t_2+R\cos \left(\omega t_2+{\phi }_{\mathrm{int}}\right)\\ y_2=R\sin \left(\omega t_2+{\phi }_{\mathrm{int}}\right)\end{array}$$

(18)

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:

$$\frac{R\sin \left(\omega t_1+{\phi }_{\mathrm{int}}\right)-0}{\left[f\cdot t_1+R\cos \left(\omega t_1+{\phi }_{\mathrm{int}}\right)-f_z\right]-f\cdot t_2}=\frac{R\sin \left(\omega t_2+{\phi }_{\mathrm{int}}\right)-0}{\left[f\cdot t_2+R\cos \left(\omega t_2+{\phi }_{\mathrm{int}}\right)\right]-f\cdot t_2}=k$$

(19)

where t₁ is the time for the cutting edge to move to $\left(x_1^{(2)},y_1^{(2)}\right)$ on the cutting path of the previous cutter tooth.

Calculate t₁ by Equation (19) and substitute it into Equation (20) to calculate $\left(x_1^{(2)},y_1^{(2)}\right)$.

$$\{\begin{array}{l}{x}_1^{(2)}=f\cdot t_1+R\cos \left(\omega t_1+{\phi }_{\mathrm{int}}\right)-f_z\\ y_1^{(2)}=R\sin \left(\omega t_1+{\phi }_{\mathrm{int}}\right)\end{array}$$

(20)

Finally, $\left(x_1^{(2)},y_1^{(2)}\right)$ and $\left(x_2^{(2)},y_2^{(2)}\right)$ are substituted into Equation (16) to calculate the cutting thickness h_c. When the cutting thickness h_c is the distance between the cutting edge and the surface which is to be machined, since $y_1^{(1)}$ is known, $x_1^{(1)}$ can be obtained by Equation (21).

$$\frac{\left(R-a_e\right)-0}{x_1^{(1)}}=\frac{R\sin \left(\omega t_2+{\phi }_{\mathrm{int}}\right)-0}{\left[f\cdot t_2+R\cos \left(\omega t_2+{\phi }_{\mathrm{int}}\right)\right]-f\cdot t_2}=k$$

(21)

The calculation process of the milling force is shown in Figure 5.

4. 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.

4.1. 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):

$$\{\begin{array}{l}{\sigma }_x=-\frac{2z}{\pi }{\displaystyle \underset{-a}{\overset{a}{\int }}\frac{p\left(s\right){\left(x-s\right)}^2}{{\left[{\left(x-s\right)}^2+z^2\right]}^2}}ds-\frac{2}{\pi }{\displaystyle \underset{-a}{\overset{a}{\int }}\frac{q\left(s\right){\left(x-s\right)}^3}{{\left[{\left(x-s\right)}^2+z^2\right]}^2}ds}\\ {\sigma }_z=-\frac{2z^3}{\pi }{\displaystyle \underset{-a}{\overset{a}{\int }}\frac{p\left(s\right)}{{\left[{\left(x-s\right)}^2+z^2\right]}^2}}ds-\frac{2z^2}{\pi }{\displaystyle \underset{-a}{\overset{a}{\int }}\frac{q\left(s\right)\left(x-s\right)}{{\left[{\left(x-s\right)}^2+z^2\right]}^2}ds}\\ {\sigma }_{xz}=-\frac{2z^2}{\pi }{\displaystyle \underset{-a}{\overset{a}{\int }}\frac{p\left(s\right)\left(x-s\right)}{{\left[{\left(x-s\right)}^2+z^2\right]}^2}}ds-\frac{2z}{\pi }{\displaystyle \underset{-a}{\overset{a}{\int }}\frac{q\left(s\right){\left(x-s\right)}^2}{{\left[{\left(x-s\right)}^2+z^2\right]}^2}ds}\end{array}$$

(22)

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:

$$p\left(s\right)=\frac{F_c\sin \mathrm{\Phi }+F_t\cos \mathrm{\Phi }}{L_{AB}w}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}q\left(s\right)={\sigma }_{AB}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}a=\frac{1}{2}{L}_{AB}$$

(23)

The normal stress, tangential stress, and loading range on the friction area between the tool tip and the workpiece are expressed as:

$$p\left(s\right)=\frac{2P_{thrust}}{\pi \left(wa\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}q\left(s\right)=\mu \left(\frac{P_{cut}}{w\cdot CA}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}a=\frac{1}{2}CA$$

(24)

4.2. 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:

$$\{\begin{array}{l}{\sigma }_x^{therm}\left(x,z\right)=-\frac{{\alpha }_{T}E}{1-2\upsilon }{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }\left(G_{xh}\frac{\partial T}{\partial x}\left(x,z\right)+G_{xv}\frac{\partial T}{\partial x}\left(x,z\right)\right)}}dxdz+\frac{2z}{\pi }{\displaystyle {\int }_0^{\infty }\frac{p\left(t\right){\left(t-x\right)}^2}{\left({\left(t-x\right)}^2+z^2\right)}}dt-\frac{{\alpha }_{T}ET\left(x,z\right)}{1-2\upsilon }\\ {\sigma }_z^{therm}\left(x,z\right)=-\frac{{\alpha }_{T}E}{1-2\upsilon }{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }\left(G_{zh}\frac{\partial T}{\partial x}\left(x,z\right)+G_{zv}\frac{\partial T}{\partial x}\left(x,z\right)\right)}}dxdz+\frac{2z^3}{\pi }{\displaystyle {\int }_0^{\infty }\frac{p\left(t\right)}{\left({\left(t-x\right)}^2+z^2\right)}}dt-\frac{{\alpha }_{T}ET\left(x,z\right)}{1-2\upsilon }\\ {\sigma }_{xz}^{therm}\left(x,z\right)=-\frac{{\alpha }_{T}E}{1-2\upsilon }{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }\left(G_{xzh}\frac{\partial T}{\partial x}\left(x,z\right)+G_{xzv}\frac{\partial T}{\partial x}\left(x,z\right)\right)}}dxdz+\frac{2z^2}{\pi }{\displaystyle {\int }_0^{\infty }\frac{p\left(t\right)\left(t-x\right)}{\left({\left(t-x\right)}^2+z^2\right)}}dt\end{array}$$

(25)

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, G_xh, G_xv, G_zh, G_zv, G_xzh and G_xzv represent the Green functions of the plane strain, and the detailed calculations refer to the study of Saif [13]. p(t) is expressed as:

$$p\left(t\right)=\frac{\alpha ET\left(x,z=0\right)}{1-2\upsilon }$$

(26)

The stress at any point on the workpiece due to mechanical and thermal loads is expressed as:

$$\left[{\sigma }_{total}\right]=\left[{\sigma }_{total}^{mech}\right]+\left[{\sigma }_{total}^{thermal}\right]$$

(27)

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 $f_z$, 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. The experimental scheme.

Case	Spindle Speed (r/min)	Cutting Width (mm)	Cutting Depth (mm)	Feed Rate (mm/min)
1	750	0.1	8	100
2	600	0.1	8	100
3	750	0.12	8	100
4	750	0.1	8	120

as an example, Figure 7 shows the stress contours in the workpiece due to the combined loads.

4.3. 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.

$$\Delta {\sigma }_{xx}=\frac{{\sigma }_{xx}^{total}}{M}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\Delta {\sigma }_{zz}=\frac{{\sigma }_{zz}^{total}}{M}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\Delta {\sigma }_{xz}=\frac{{\sigma }_{xz}^{total}}{M}$$

(28)

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:

$$VM=\frac{3}{2}\left(s_{ij}-{\alpha }_{ij}\right)\left(s_{ij}-{\alpha }_{ij}\right)-K^2$$

(29)

where $s_{ij}={\sigma }_{ij}-\left({\sigma }_{kk}/3\right){\delta }_{ij}$, where $s_{ij}$ is deviatoric stress, ${\delta }_{ij}$ is Kronecker delta. $\langle \rangle$ is the Macaulay bracket as, $\langle x\rangle =0.5\left(x+\left|x\right|\right)$, $n_{ij}$ is the unit normal vector pointing outward the yield surface, which can be expressed as:

$$n_{ij}=\frac{s_{ij}-{\alpha }_{ij}}{\sqrt{2}k}$$

(30)

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:

$$\Delta {\epsilon }_{ij}^e=\frac{1}{E}\left[\left(1+\nu \right)\Delta {\sigma }_{ij}-{\delta }_{ij}\nu \Delta {\sigma }_{kk}\right]$$

(31)

where E is the elastic modulus and $\nu$ is Poisson’s ratio and $\Delta {\sigma }_{kk}=\Delta {\sigma }_{xx}+\Delta {\sigma }_{yy}+\Delta {\sigma }_{zz}$.

According to the plastic flow criterion, plastic strain is expressed as:

$$\Delta {\epsilon }_{ij}^p=\frac{1}{h}\langle \Delta s_{kl}{n}_{kl}\rangle n_{ij}$$

(32)

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:

$$\psi =1-\exp \left(-\kappa \frac{3h}{2G}\right)$$

(33)

where $\kappa$ is a constant, G is the shear modulus of elasticity. $\Delta {\sigma }_{xx}$ and $\Delta {\sigma }_{yy}$ are stress increments in the cutting direction and perpendicular to cutting direction respectively, which can be expressed as:

$$\{\begin{array}{l}\Delta {\epsilon }_{xx}=\frac{1}{E}\left[\Delta {\sigma }_{xx}-\nu \left(\Delta {\sigma }_{yy}+\Delta {\sigma }_{zz}^{*}\right)\right]+{\alpha }_T\Delta T+\frac{1}{h}\left(\Delta {\sigma }_{xx}{n}_{xx}+\Delta {\sigma }_{yy}{n}_{yy}+\Delta {\sigma }_{zz}^{*}{n}_{zz}+2\Delta {\tau }_{xz}^{*}{n}_{xz}\right)n_{xx}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\psi \left(\frac{1}{E}\left[\Delta {\sigma }_{xx}^{*}-\nu \left(\Delta {\sigma }_{yy}+\Delta {\sigma }_{zz}^{*}\right)\right]+{\alpha }_T\Delta T+\frac{1}{h}\left(\Delta {\sigma }_{xx}^{*}{n}_{xx}+\Delta {\sigma }_{yy}{n}_{yy}+\Delta {\sigma }_{zz}^{*}{n}_{zz}+2\Delta {\tau }_{xz}^{*}{n}_{xz}\right)n_{xx}\right)\\ \Delta {\epsilon }_{yy}=\frac{1}{E}\left[\Delta {\sigma }_{yy}-\nu \left(\Delta {\sigma }_{xx}+\Delta {\sigma }_{zz}^{*}\right)\right]+{\alpha }_T\Delta T+\frac{1}{h}\left(\Delta {\sigma }_{xx}{n}_{xx}+\Delta {\sigma }_{yy}{n}_{yy}+\Delta {\sigma }_{zz}^{*}{n}_{zz}+2\Delta {\tau }_{xz}^{*}{n}_{xz}\right)n_{yy}=0\end{array}$$

(34)

where $\Delta {\sigma }_{xx}^{*}$, $\Delta {\sigma }_{yy}^{*}$, $\Delta {\sigma }_{zz}^{*}$ and $\Delta {\tau }_{xz}^{*}$ are the stresses obtained by Equation (27). Integrating the history of stress increment can obtain the stress distribution inside the workpiece after the stress loading process.

4.4. Stress Release Process

After the stress loading process is completed, the stress release process is performed to satisfy the boundary conditions. And the residual stress is the stress remaining inside the workpiece. According to the model proposed by Merwin and Johnson, the internal stress and strain of the workpiece after stress release process need to satisfy the following boundary conditions:

$$\begin{array}{l}{\epsilon }_x^r=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_y^r=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_z^r=f\left(z\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\gamma }_{xz}^r=f\left(z\right)\\ {\sigma }_x^r=f\left(z\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\sigma }_y^r=f\left(z\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\sigma }_z^r=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\tau }_{xz}^r=0\end{array}$$

(35)

The release increment for each step is expressed as:

$$\Delta {\sigma }_z=-\frac{{\sigma }_z^r}{M}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\Delta {\tau }_{xz}=-\frac{{\tau }_{xz}^r}{M}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\Delta {\epsilon }_{xx}=-\frac{{\epsilon }_{xx}^r}{M}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\Delta T=-\frac{T^r}{M}$$

(36)

As with the stress loading process, it is also necessary to determine whether the workpiece has yielded during the stress release process. When VM < 0 or VM = 0 and ds_ijn_ij ≥ 0, calculate the stress increment under elastic release by Equation (37):

$$\{\begin{array}{l}\Delta {\sigma }_x^r=\frac{E\Delta {\epsilon }_x+\left(1+\upsilon \right)\left(\Delta {\sigma }_z\upsilon -E{\alpha }_T\Delta T\right)}{\left(1-{\upsilon }^2\right)}\\ \Delta {\sigma }_y^r=\frac{\upsilon E\Delta {\epsilon }_x+\left(1+\upsilon \right)\left(\Delta {\sigma }_z\upsilon -E{\alpha }_T\Delta T\right)}{\left(1-{\upsilon }^2\right)}\end{array}$$

(37)

When VM = 0 and ds_ijn_ij < 0, calculate the stress increment under plastic release by Equation (38):

$$\{\begin{array}{l}\Delta {\sigma }_y^r=\frac{\left(-\frac{\upsilon }{E}+\frac{1}{h}{n}_x{n}_y\right)\left(C-{\alpha }_T\Delta T\right)-\left(\frac{1}{E}+\frac{1}{h}{n}_x{n}_x\right)\left(D-{\alpha }_T\Delta T\right)}{{\left(-\frac{\upsilon }{E}+\frac{1}{h}{n}_x{n}_y\right)}^2-\left(\frac{1}{E}+\frac{1}{h}{n}_x{n}_x\right)\left(\frac{1}{E}+\frac{1}{h}{n}_y{n}_y\right)}\\ \Delta {\sigma }_x^r=\frac{D-\left(\frac{1}{E}+\frac{1}{h}{n}_y{n}_y\right)\Delta {\sigma }_y^r-{\alpha }_T\Delta T}{-\frac{\upsilon }{E}+\frac{1}{h}{n}_x{n}_y}\end{array}$$

(38)

where

$$\{\begin{array}{l}C=\Delta {\epsilon }_x^r+\left(\frac{\upsilon }{E}-\frac{1}{h}{n}_x{n}_z\right)\Delta {\sigma }_z^{*}-\frac{2}{h}\Delta {\tau }_{xz}^{*}{n}_{xz}{n}_x\\ D=\left(\frac{\upsilon }{E}-\frac{1}{h}{n}_y{n}_z\right)\Delta {\sigma }_z^{*}-\frac{2}{h}\Delta {\tau }_{xz}^{*}{n}_{xz}{n}_y\end{array}$$

(39)

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.

5. Experimental Validation of Milling Force Model and Residual Stress Model

5.1. 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.

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. Material parameters of workpiece and tool [14].

Material Parameter	Ti6Al4V	WC
ρ (kg/m3)	4430	15630
K (W/m·°C)	7.039e0.0011T	55
Cp (J/kg °C)	505.64e0.0007T	0.32T + 1.32 × 10−4
Tmelt (°C)	1650	-

.

5.2. 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 ${\epsilon }_r^{\ast }$ of the absolute value of the error between the experimental data and the predicted value, the prediction accuracy of the model is quantitatively evaluated. ${\epsilon }_r^{\ast }$ is calculated by the following equation:

$${\epsilon }_r^{\ast }=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|F_m^i-F_e^i\right|}$$

(40)

where $F_m^i$ is the predicted value of milling force, $F_e^i$ is the experimental value of milling force, n is the number of measured experimental data. The prediction error ${\epsilon }_r^{\ast }$ is shown in

Table 4. Prediction error ${\epsilon }_r^{\ast }$ in cases 1–4.

Case	${\mathit{\epsilon }}_{\mathit{r}}^{\ast }$ in X Direction (N)	${\mathit{\epsilon }}_{\mathit{r}}^{\ast }$ in Y Direction (N)	${\mathit{\epsilon }}_{\mathit{r}}^{\ast }$ in Z Direction (N)
1	10.61	37.64	2.53
2	22.06	51.49	4.57
3	15.33	52.75	2.60
4	13.39	37.25	3.06

. 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 F_Y is large, and ${\epsilon }_r^{\ast }$ can reach 52.75 N. When entering the stable cutting stage, F_Y 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 F_Y 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 F_X and F_Y, but it has little effect on the milling force F_Z. Therefore, the experimental value of the cutting force F_Z 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 F_Y 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. t₁ and t₂ 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 t₂ is higher in the axial direction than that at time t₁. 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].

5.3. 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 ${\sin }^2\psi$, 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. The XRD stress test parameters for Ti6Al4V.

Plane (hkl)	Radiation	Bragg’s Angle 2θ	1/2 S2 (MPa−1)	S1 (MPa−1)	Voltage (kV)	Current (mA)
(213)	Cu-Kα	141.3°	12.18 × 10−6	−3.09 × 10−6	40	40

.

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.

5.4. 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 F_Y 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.

6. 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 F_X and F_Z, it can be found that there is still a wave value of convergence in the milling force F_Y 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 F_X and F_Z 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.