Prediction of Surface Location Error Considering the Varying Dynamics of Thin-Walled Parts during Five-Axis Flank Milling

Abstract

Surface location error (SLE) caused by forced vibration is a key factor to determine the quality of the finished part. When machining thin-walled structures with sculptured surfaces, the complicated milling process is significantly influenced by the vibration due to the flexibility of the part. The dynamics of the part are dominant and vary with the material removal during machining. This paper presents a prediction method of SLE considering the varying dynamics of thin-walled parts in five-axis flank milling. The in-process part is decomposed into unmachined and machined portions, which are both modelled based on the thin-plate theory. The dynamics models of the two portions are coupled using the substructure method. Coordinate transformation based on the screw theory and the general cutting dynamics model for five-axis flank milling is employed to transform the cutting force vectors and frequency response function (FRF) to the same coordinate system for the prediction of SLE. The proposed method is validated with five-axis flank milling tests and SLE measurements on a thin-walled twisted part. It is shown that the average error of the proposed method for SLE prediction is less than 5 μm, and the calculation is almost 8 times faster than the typical finite element method.

1. Introduction

Thin-walled parts with sculptured surfaces are widely used in the aerospace industry, which are always machined through five-axis flank milling. These parts normally contain thin walls (1.5–4 mm thickness) and require good surface quality. Forced vibration induced by the flexibility of parts is a significant factor influencing the dimensional quality of the parts [1], which significantly affect their service performance [2]. The prediction of surface location error (SLE) of thin-walled parts represents the basis of determining the proper cutting parameters, which aims at increasing the productivity and reducing the costs of a milling operation [3]. However, the varying dynamics of the part during machining and the complexity of the milling process both increase the difficulty of SLE prediction.

The prediction of SLE is based on the dynamics of the part. Li [4] and Ismail [5] employed experimental modal analysis (EMA) to measure the dynamics of flexible parts. However, the material of the part is removed during machining, and the dynamics of the part vary with this process. Finite element (FE) models are employed to predict the dynamics of a part to avoid the interruption of the machining process caused by EMA. Tsai et al. [6] and Adetoro et al. [7] employed the FE model to analyze the deflection of a thin-walled part during machining, and pointed out that the dynamics of the part influenced the machining process significantly, whereas with the rapid development of industries, the dimension of parts becomes larger and the element’s number of the corresponding FE model is huge, which is time-consuming for dynamics prediction. Thus, computationally efficient methods simplified by either elements [8] or computational processes [9,10] are proposed to solve pertinent problems. Meshreki et al. [11] presented a dynamics model based on Rayleigh’s energy method to reduce computational time for thin-walled pockets; then, they [12,13] updated the model to take into account the continuous changes in the part’s thickness. Tuysuz et al. [14,15] used reduced-order substructure and perturbation methods to further improve the efficiency of computing dynamics for thin-walled parts in milling. Ma et al. [16] proposed a method to predict the dynamics of thin-walled parts based on the equivalent plate and the external damping effect. However, the major puzzle in simulation is the contradiction of efficiency, accuracy, and convenience [17], and it is hard to offer the efficiency to compute dynamics and the versatility to model parts with complex geometries simultaneously.

The relationship between the machining quality and process was investigated by Kline [18], Tlusty [19], and Sutherland [20] early, and then Schmitz et al. [21] studied the prediction method of SLE systematically based on the frequency response function (FRF) of a tool tip. However, in the milling of thin walls, the flexibility of the parts determines the SLE of the finished parts. Some scholars studied the surface quality based on thin plates. Altintas et al. [22] predicted surface errors based on the dynamics model of a cantilever plate structure and pointed out that the dynamics of a part are dominant in the flank milling of thin-walled structures. Budak et al. [23] modelled the milling process of cantilevered plates, and showed that the vibration between the part and cutter influences the surface quality of the part significantly. Ratchev et al. [24] predicted surface errors caused by the deflection of low-stiffness parts through an adaptive flexible theoretical deflection model. Moreover, Sofuoglu et al. [25] developed regression–multicriteria decision-making hybrid models and determined the optimum cutting conditions without chatter vibrations. Li et al. [26] divided the factors determining machining quality into kinematic and stochastic portions. Misaka et al. [27] predicted surface roughness by integrating the data obtained through an experiment and an analytical model to improve computational efficiency. Meanwhile, Ringgaard et al. [28] demonstrated that the relative forced vibration between the tool and the part must be minimized to improve the machining quality. Jiao et al. [29] developed a surface roughness model considering the forced vibration in two directions. Li et al. [30] proposed a prediction model of deformation induced by cutting forces based on the substructure method for the thin-walled parts in flank milling. Other scholars predict SLE with the consideration of tool and part deflection [31], process damping [32], tool dynamic stiffness variation [33], etc. However, past prediction methods have mainly been based on simple cutting conditions. A huge gap still exists between reality and simulation due to the limitation of computational efficiency and simplified models [34]. The SLE prediction method of thin-walled parts with the consideration of varying dynamics in five-axis flank milling has not been reported in the literature.

A prediction method of SLE considering the varying dynamics of thin-walled parts during five-axis flank milling is proposed in the paper. The in-process part is represented by the machined and unmachined portions, and the substructure method is employed to couple the dynamics model of two portions, as presented in Section 2. The SLE of thin-walled parts in five-axis flank milling is predicted in Section 3 based on the developed dynamics model and screw theory. The proposed method is verified in Section 4. Section 5 concludes the paper.

2. Dynamics Model of Thin-Walled Part during Milling

2.1. Machined and Unmachined Portions of In-Process Part

The in-process part is composed of unmachined and machined portions. The thickness of the machined portion is smaller than the unmachined portion. The prediction of SLE depends on the FRFs of the cutting point on the part, which is determined by the dynamics of two portions. The decomposition of the in-process part is illustrated in Figure 1.

The machined and unmachined portions are regarded as substructures of the in-process part. The dynamics model of the in-process part ([$\mathsf{\Theta }$]) is obtained by coupling the dynamics of two substructures, expressed by

$$\left[\mathsf{\Theta }\right]=\left[{\mathsf{\Theta }}^{un}\right]+\left[{\mathsf{\Theta }}^{ma}\right]$$

(1)

where [$\mathsf{\Theta }$^un] and [$\mathsf{\Theta }$^ma] represent the dynamics of the unmachined and machined portions of the part, respectively.

2.2. Dynamics Model of Thin Wall

For thin-walled parts, the thickness is so small compared to the length and height. Therefore, in the thickness direction, the value of stress can be regarded as constant. The part is represented by the neutral plane, and it is meshed by the appropriate elements. Each node of the element contains three degrees of freedom (DOF) based on Kirchhoff’s thin-plate theory. The virtual work principle is used to develop the stiffness matrix K^e and mass matrix M^e of the element, expressed as

$$M^e={\displaystyle \int {}_{S^e}w\tau N^{T}NdS}$$

(2)

$$K^e={\displaystyle \int {}_{V^e}{\tau }^3{B}^{T}DBdV}$$

(3)

where w is the product of the density of the material and the area of the element. $\tau$ is the thickness of the part. N is the shape function depending on the type of element, and B is derived from N. D is the modified bending modulus, calculated by

$$D=\frac{E}{12(1-{\upsilon }^2)}\left[\begin{array}{ccc}1& \upsilon & 0\\ \upsilon & 1& 0\\ 0& 0& (1-\upsilon )/2\end{array}\right]$$

(4)

where E is the elastic modulus, and υ is the Poisson’s ratio.

The accuracy of Rayleigh damping is enough for establishing the damping matrix. Meanwhile, the computational efficiency is high. Based on the damping coefficients (α and β), Rayleigh damping is the function of mass and stiffness matrices, given by

$$C^e=\alpha M^e+\beta K^e$$

(5)

The matrices of elements are assembled on nodes to comprise the matrices of the part. The dynamics model of a thin wall is expressed as

$$M\ddot{X}+C\dot{X}+KX=F$$

(6)

where M, C, and K are the mass, damping, and stiffness matrices of the part. X and F are the displacement and force vectors, respectively.

2.3. Coupling Using the Substructure Method

The unmachined and machined portions are defined as substructures of the part, and the dynamics of two portions are coupled by the substructure coupling method. The dynamics model of the substructure is expressed as

$$M_r^{}{\ddot{X}}_r^{}+C_r^{}{\dot{X}}_r^{}+K_r^{}{X}_r^{}=F_r^{}\begin{array}{cc}& (r=\begin{array}{ccc}un& \mathrm{or}& ma\end{array})\end{array}$$

(7)

where r is un or ma, denoting the unmachined portion and the machined portion, and this notation is used in the following equations. By using the Laplace transform, Equation (7) is transformed into the frequency domain from the time domain.

$$S_r(\omega )={\omega }^2{M}_r+\omega C_r+K_r$$

(8)

where S_r(ω) is the dynamic stiffness matrices of the substructure. [$\mathsf{\Theta }$^un] and [$\mathsf{\Theta }$^ma] in Equation (1) are rewritten as

$$\left[{\mathsf{\Theta }}^r\right]=S_r{(\omega )}^{-1}$$

(9)

The coupling process of two substructures is shown in Figure 2. The interface nodes and internal nodes of the unmachined portion and machined portion are labelled as #a, #b, and #c. Each node contains three degrees of freedom (DOF) based on the dynamics model of a thin wall developed in Section 2.2.

The dynamic stiffness matrix S_r, in which the ω is dropped for clarity, is rewritten as the matrixes of different DOF sets, populated by

$$S_r^{}=\left\{\begin{array}{c}{S}_{un}^{}=\left[\begin{array}{cc}{S}_{un}^{cc}& S_{un}^{cb}\\ S_{un}^{bc}& S_{un}^{bb}\end{array}\right]\\ S_{ma}^{}=\left[\begin{array}{cc}{S}_{ma}^{bb}& S_{ma}^{ab}\\ S_{ma}^{ba}& S_{ma}^{aa}\end{array}\right]\end{array}\right.$$

(10)

where the superscripts a, b, and c represent the DOF sets of #a, #b, and #c respectively. The force vector F_r is rewritten as the vector of different DOF sets, expressed by

$$F_r=\left\{\begin{array}{c}{F}_{un}=\left\{\begin{array}{c}{F}_{un}^c\\ F_{un}^b\end{array}\right\}\\ F_{ma}=\left\{\begin{array}{c}{F}_{ma}^b\\ F_{ma}^a\end{array}\right\}\end{array}\right.$$

(11)

The dynamics model of two substructures is composed as

$$\left[\begin{array}{cccc}{S}_{un}^{cc}& S_{un}^{cb}& 0& 0\\ S_{un}^{bc}& S_{un}^{bb}& 0& 0\\ 0& 0& S_{ma}^{bb}& S_{ma}^{ab}\\ 0& 0& S_{ma}^{ba}& S_{ma}^{aa}\end{array}\right]\left\{\begin{array}{c}{X}_{un}^c\\ X_{un}^b\\ X_{ma}^b\\ X_{ma}^a\end{array}\right\}=\left\{\begin{array}{c}{F}_{un}^c\\ F_{un}^b\\ F_{ma}^b\\ F_{ma}^a\end{array}\right\}$$

(12)

Based on the compatibility of interface displacements, the transformation of DOF vectors from the structure to the substructure is represented by

$$\left\{\begin{array}{c}{X}_{un}^c=X_{}^c\\ X_{un}^b=X_{}^b\\ X_{ma}^b=X_{}^b\\ X_{ma}^a=X_{}^a\end{array}\right.\Rightarrow \left[\begin{array}{lll}I\hfill & 0\hfill & 0\hfill \\ 0\hfill & I\hfill & 0\hfill \\ 0\hfill & I\hfill & 0\hfill \\ 0\hfill & 0\hfill & I\hfill \end{array}\right]\left\{\begin{array}{c}{X}_{}^c\\ X_{}^b\\ X_{}^a\end{array}\right\}=\left\{\begin{array}{c}{X}_{un}^c\\ X_{un}^b\\ X_{ma}^b\\ X_{ma}^a\end{array}\right\}$$

(13)

Based on the force equilibrium, the transformation of force vectors from the substructure to the structure is written as

$$\left\{\begin{array}{c}{F}_{un}^c=F_{}^c\\ F_{un}^b+F_{ma}^b=F_{}^b\\ F_{ma}^c=F_{ma}^c\end{array}\right.\Rightarrow \left[\begin{array}{llll}I\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & I\hfill & I\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & I\hfill \end{array}\right]\left\{\begin{array}{c}{F}_{un}^c\\ F_{un}^b\\ F_{ma}^b\\ F_{ma}^b\end{array}\right\}=\left\{\begin{array}{c}{F}_{}^c\\ F_{}^b\\ F_{}^a\end{array}\right\}$$

(14)

Substituting Equations (13) and (14) into Equation (12) yields

$$\left[\begin{array}{llll}I\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & I\hfill & I\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & I\hfill \end{array}\right]\left[\begin{array}{cccc}{S}_{un}^{cc}& S_{un}^{cb}& 0& 0\\ S_{un}^{bc}& S_{un}^{bb}& 0& 0\\ 0& 0& S_{ma}^{bb}& S_{ma}^{ab}\\ 0& 0& S_{ma}^{ba}& S_{ma}^{aa}\end{array}\right]\left[\begin{array}{lll}I\hfill & 0\hfill & 0\hfill \\ 0\hfill & I\hfill & 0\hfill \\ 0\hfill & I\hfill & 0\hfill \\ 0\hfill & 0\hfill & I\hfill \end{array}\right]\left\{\begin{array}{c}{X}_{}^c\\ X_{}^b\\ X_{}^a\end{array}\right\}=\left\{\begin{array}{c}{F}_{}^c\\ F_{}^b\\ F_{}^a\end{array}\right\}$$

(15)

The dynamic stiffness of structure S_un+ma can be expressed by

$$S_{un+ma}^{}=\left[\begin{array}{llll}I\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & I\hfill & I\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & I\hfill \end{array}\right]\left[\begin{array}{cccc}{S}_{un}^{cc}& S_{un}^{cb}& 0& 0\\ S_{un}^{bc}& S_{un}^{bb}& 0& 0\\ 0& 0& S_{ma}^{bb}& S_{ma}^{ab}\\ 0& 0& S_{ma}^{ba}& S_{ma}^{aa}\end{array}\right]\left[\begin{array}{lll}I\hfill & 0\hfill & 0\hfill \\ 0\hfill & I\hfill & 0\hfill \\ 0\hfill & I\hfill & 0\hfill \\ 0\hfill & 0\hfill & I\hfill \end{array}\right]$$

(16)

Based on the Equation (16), the dynamics of the in-process part ([$\mathsf{\Theta }$]) are represented as

$$\left[\mathsf{\Theta }\right]=S_{un+ma}^{}{}^{-1}$$

(17)

The FRF of the part in the cutting point is calculated by

$$H_p=\left[\mathsf{\Theta }\left(q\right)\right]={\displaystyle \sum _{i=1}^{n{}^{\prime }}\frac{{\mathsf{\Psi }}_{iq}{\mathsf{\Psi }}_{iq}^T}{({\omega }_i^2-{\omega }^2)+2\mathrm{j}{\zeta }_i{\omega }_i\omega }}$$

(18)

where q is the DOF in the thickness direction of the cutting point. ω_i, Ψ_iq, and ζ_i are the eigenvalue, eigenvector, and modal damping ratio calculated from Equation (17), and n′ is the number of natural modes.

3. Prediction of Surface Location Error

3.1. Coordinate Transformation in Five-Axis Flank Milling Based on Screw Theory

The part coordinate system (PCS) and the tool coordinate system (TCS) are established as illustrated in Figure 3. The FRFs of cutting points on the part are predicted in the PCS, and the TCS is used to analyze the dynamic response of the tool tip. However, the forced vibration should be calculated in the PCS. A simple but efficient method based on the general cutting dynamics model [4] and the screw theory is employed to transform the TCS into the PCS, expressed as

$${}^{p}T{}_t({\theta }_1,{\theta }_2)=e^{s_1\left[{\omega }_1\right]{\theta }_1}\cdot e^{s_2\left[{\omega }_2\right]{\theta }_2}\cdot {}^{p}T{}_t(0)$$

(19)

where θ₁ and θ₂ are the rotation angles [35] of two rotary axis in five-axis flank milling. The two rotary axes can be distinguished based on being close to either the part or the tool, and ω₁ and ω₂ are the unit rotary vectors of the two rotary axes. s₁ and s₂ are −1 if the rotary axis is mounted on the table side; otherwise, they are 1. ^pT_t (0) is the orientation of the TCS in the PCS at the beginning of machining [36].

The FRF of the tool tip can be transformed in the PCS based on Equation (19), and the FRF in the cutting point is calculated by

$$H_{tp}={}^{p}T{}_t\cdot H_t\cdot {}^{t}T{}_p+H_p$$

(20)

where H_t and H_p are the FRFs of the tool and the part, respectively. H_t can be obtained through modal testing, and H_p is calculated by Equation (18).

3.2. Surface Location Error in Five-Axis Flank Milling

A cutting force is required to predict the SLE. A typical end mill is illustrated in Figure 4. Based on a tiny axial length d_z, the cutting portion of the flute is represented by m disk elements. The local coordinate system O_lX_lY_lZ_l (LCS) is created attached to the cutting edge, where axes X_l, Y_l, and Z_l are axial, radial, and tangential along the tool as shown in Figure 4.

For the ith element of the jth tooth, the tangential, radial, and axial cutting forces at t time with the rotation angle φ_i,j(t) are expressed as

$$\left\{\begin{array}{c}\mathrm{d}{F}_{\mathrm{t}{}^{\prime }}=\chi ({\phi }_{i,j}(t))\cdot \left(K_{\mathrm{tc}}\cdot h({\phi }_{i,j}(t))+K_{\mathrm{te}}\right)\cdot \mathrm{d}z\\ \mathrm{d}{F}_{\mathrm{r}{}^{\prime }}=\chi ({\phi }_{i,j}(t))\cdot \left(K_{\mathrm{rc}}\cdot h({\phi }_{i,j}(t))+K_{\mathrm{re}}\right)\cdot \mathrm{d}z\\ \mathrm{d}{F}_{\mathrm{a}{}^{\prime }}=\chi ({\phi }_{i,j}(t))\cdot \left(K_{\mathrm{ac}}\cdot h({\phi }_{i,j}(t))+K_{\mathrm{ae}}\right)\cdot \mathrm{d}z\end{array}\right.$$

(21)

where K is the cutting force coefficient. The subscript t, a, and r represent the tangential, axial, and radial directions, respectively, and the subscript c and e are employed to distinguish the shear and plow cutting force. h(φ_i,j(t)) is the chip thickness at t time. $\chi$(φ_i,j(t)) is employed to determine whether the cutting force is developing, expressed by

$$\begin{array}{ccc}& \chi ({\phi }_{i,j}(t))=\left\{\begin{array}{cc}1,& {\theta }_{st}⩽{\phi }_{i,j}(t)⩽{\theta }_{ex}\\ 0,& \mathit{otherwise}\end{array}\right.& \end{array}$$

(22)

where θ_st and θ_ex stand for the angles of entry and exit, respectively.

The tangential, radial, and axial cutting forces are developed in the LCS, and transformed to in the TCS based on the transition matrix R.

$$\left[\begin{array}{l}\mathrm{d}{F}_{xt}\hfill \\ \mathrm{d}{F}_{yt}\hfill \\ \mathrm{d}{F}_{zt}\hfill \end{array}\right]=R\left[\begin{array}{l}\mathrm{d}{F}_{\mathrm{t}{}^{\prime }}\hfill \\ \mathrm{d}{F}_{\mathrm{r}{}^{\prime }}\hfill \\ \mathrm{d}{F}_{\mathrm{a}{}^{\prime }}\hfill \end{array}\right]$$

(23)

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

(24)

The total cutting force in the TCS is calculated by

$$F_t=\left[\begin{array}{c}{F}_{xt}\\ F_{yt}\\ F_{zt}\end{array}\right]={\displaystyle \sum _{j=1}^n{\displaystyle \sum _{i=1}^m\left[\begin{array}{l}\mathrm{d}{F}_x\hfill \\ \mathrm{d}{F}_y\hfill \\ \mathrm{d}{F}_z\hfill \end{array}\right]}}$$

(25)

where n is the tooth number of the cutting tool. The cutting force is defined in the TCS, but forced vibration displacement is calculated in the PCS, which can be transformed by ^pT_t, expressed as

$$F_{\mathrm{PCS}}(t)={}^{p}T{}_t\cdot F_t$$

(26)

The forced vibration displacement is calculated by

$$X(t)=f_{\mathit{fourier}}^{-1}\left\{H_{tp}\cdot F_{\mathrm{PCS}}(\omega )\right\}$$

(27)

where f⁻¹_fourier represents the inverse Fourier transformation operator. With the prediction of forced vibration displacement, the SLE at the position of the tool tip on the part can be predicted by

$$x(s)=X\left(\frac{s}{f_t}\right)$$

(28)

where s is the length of the distance on the part from the cutting beginning to the position of SLE prediction; f_t is the feed rate.

The overall procedure is summarized as shown in Figure 5.

4. Experimental Verification

4.1. Experiment Design

A thin-walled twisted benchmark with dimensions of 100 mm × 45 mm × 6 mm was designed to verify the proposed SLE prediction method. The material is Al7075 with the following properties: E = 69 × 10⁹ Pa, μ = 0.3, and ρ = 2750 kg/m³. The α and β, which are the proportional damping constants of the benchmark and identified by experiments, are 0.84 and 5.93 × 10⁻⁸, respectively. The thickness of the benchmark was machined from 6 mm to 4 mm by using flank milling as shown in Figure 6.

The impact testing was used to obtain the FRFs of the tool tip, as shown in Figure 7. The solid carbide end mill with 3 teeth and a diameter of 16 mm was equipped to a five-axis table-tilting machine center (model: DMU50) (DMG MORI, Inc., bielefeld, Germany). The tool tip was impacted by the hammer with the model of PCB 086C03 (PCB Piezotronics, Inc., Depew, NY, USA), and an accelerometer with the model of PCB 352C23 (PCB Piezotronics, Inc., Depew, NY, USA) was employed to obtain the vibration responses. The acquisition system of model LMS SCADA305 (LMS, Inc., Reuven, Belgium) was used to analyze the impact and response data. To reduce the inaccuracy caused by accident, the impact testing was repeated 5 times and the mean value was chosen as the final measured FRFs.

The same experimental setups were also employed to obtain the FRFs of the part during machining. The FRFs of points on the part were different, and the points A1, B1, C1, and D1 with a uniform distance, which equaled the axial depth of the cut to keep the FRF testing only at the position of the tool tip on the part, were chosen for the following experiments. The testing points and the corresponding paths are shown in Figure 8. The testing point A1 was measured without material removal. For other testing points, cutting stopped for FRF measurements after the previous path was finished.

Figure 9 shows the comparison of FRFs between the tool tip and the part, where the dynamics of the part are more dominant. Therefore, only the FRFs of the part were employed to predict the SLE. The radial depth of cut was assigned as 1mm based on the thickness of the benchmark and the axial depth of cut was set to 5 mm. The feed rate was 0.1 mm/tooth. Based on the defined cutting conditions and chatter stability lobes, the chatter-free spindle speed for all four paths was chosen as 6000 rpm.

The SLEs at the position of the tool tip on the part were chosen to be predicted and measured for clarity, and the corresponding paths are illustrated in Figure 8. For the finished part, the SLEs in the four paths were measured by the trilinear coordinates measuring instrument (model: Global Classic SR 05.07.05) (Hexagon AB, Inc., Stockholm, Sweden) as shown in Figure 10. On each path, measuring points were chosen with an interval of 2 mm. Coordinate values on both sides of each point were measured to calculate the SLEs by comparing with the part model.

4.2. Verification of Dynamics Model for Thin-Walled Part in Five-Axis Flank Milling

For comparison purposes, the FRFs were also calculated using the finite element (FE) model developed by commercial software Abaqus 6.14 (Dassault Systèmes, Inc., Paris, France) and an eight-node solid brick element was employed to develop the FE model. The exciting force was set to a value of 1 N, and clamping boundary conditions were applied on the base of the benchmark. The element size was assigned as 1 mm. From an investigation of the element size on the FRF prediction results, little influence was found when using even finer elements. The parameters used to predict FRFs, containing material properties, dimensions, boundary conditions, and proportional damping constants, were substituted into the dynamics model. The FRFs were predicted using two methods and the four points A1, B1, C1, and D1 were selected for comparison. Figure 11 shows the experimental and predicted results of the FRF. The comparisons of natural frequency between the dynamics models and the experimental measurements are illustrated in

Table 1. Natural frequency comparison between experiment and prediction.

Point	Experiment [Hz]	Proposed Method [Hz]	Error [%]	FE Method [Hz]	Error [%]
1	2687	2779	3.42	2749	2.31
2	2696	2782	3.16	2753	2.11
3	2768	2849	2.91	2795	0.98
4	2762	2817	1.97	2790	1.01

. The proposed model can reach a similar accuracy to the typical FE model.

The computational efficiency can be reflected from the DOF number of the model. The larger of DOF number, the less efficiency of computation. Comparing with the FE model, a significant reduction in the DOFs can be realized by the proposed model as shown in

Table 2. DOFs’ comparison of two models.

Point	FE Method (ABAQUS)	Proposed Method
1	93,324	13,332
2	90,294	13,332
3	87,296	13,332
4	84,234	13,332

. The variation in the part’s geometry during machining is represented by the change in thickness in the corresponding elements for the proposed model, which avoids the updates of the model. The total DOFs in the FE method are almost 8 times those of the proposed model, which proves that the time cost in computation can be reduced significantly by the proposed model.

4.3. Verification of Prediction Method for Surface Location Error

In order to make a comparison, the case without varying the dynamics of the part was also performed to predict the FRFs and SLE as well. Both input parameters were the same for these two methods. The points to predict the FRFs in each path are illustrated in Figure 12, and the calculated FRFs are shown in Figure 13. The natural frequency variation in each path was reflected by the model with material removal, and the changing range is illustrated in

Table 3. Natural frequency variation of part during milling.

Path	Minimum [Hz]	Maximum [Hz]	Varying Ratio [%]
1	2612	2779	6.39
2	2782	2847	2.34
3	2811	2849	1.35
4	2713	2817	3.83

. However, the difference in the FRFs of the points in each path is only shown in magnitude for the model without material removal.

Figure 14 plots the SLE obtained by the experiment and the methods with and without considering material removal. It can be seen that the errors in the SLE on the points of each path were all lower for the proposed method, and the average value of SLE errors was reduced from ~10 μm to ~5 μm as illustrated in

Table 4. SLE comparison.

Path	Experiment Value [μm]	Method with Considering Material Removal		Method without Considering Material Removal
		Value [μm]	Error [μm]	Value [μm]	Error [μm]
1	35.39	34.24	1.15	25.58	9.81
2	26.58	21.76	4.82	13.31	13.27
3	20.11	15.46	4.65	12.16	7.95
4	15.53	11.87	3.66	4.76	10.77

. The complex changes in chip thickness at the initial and final positions in the processing direction led to a significant increase in SLE in the corresponding positions.

5. Conclusions

This article presents a computationally efficient and practical method for SLE prediction considering the varying dynamics of thin-walled parts in five-axis flank milling, including a computationally efficient dynamics model and a practical SLE prediction method. The detailed conclusions are summarized as follows:

(1) A new, computationally efficient dynamics model of thin-walled parts considering material removal is developed based on the thin-plate theory and substructure method. The accuracy of the proposed model is verified by the FRF experiments of four points in a thin-walled part during milling, and the errors of natural frequencies are all less than 5%. Meanwhile, the computational efficiency is improved ~8 times by comparing the proposed model with a typical finite element model with a similar accuracy.

(2) A practical SLE prediction method with the consideration of varying dynamics in five-axis flank milling of thin-walled parts is proposed. The coordinate transformation of FRFs and cutting force is realized depending on screw theory and the general cutting dynamics model, which is simple and efficient. The SLE is predicted depending on the forced vibration in five-axis flank milling. A thin-walled part with a sculptured surface is machined and the SLE of four paths is measured to verify the proposed method. It is shown that the average error of the SLE is reduced from ~10 μm to ~5 μm by considering the varying dynamics of the thin-walled part in five-axis flank milling.