## 1 Introduction  
With the rapid development in the manufacturing sector, thin-walled parts with a complex curved surface and variable thickness are widely used in the aerospace, automobile, mold, and other industries because of the advantages of the lightweight and high-specific strength. However, these thin-walled parts have low stiffness, which triggers serious vibrational problems in the machining process. Therefore, the improvement of the processing efficiency and ensuring processing quality remain hot topics for both academia and the industry. In the multi-axis milling process of such parts, the selection of existing cutting parameters and the formulation of processing technology relies heavily on the actual processing experience; and the parameters selected as such tend to be conservative. Although it can ensure the processing quality, the low efficiency and long production cycle cannot meet the needs of the manufacturing industry. Under the premise of ensuring the processing quality, it is an urgent problem to select larger cutting parameters [1,2], where the traditional two degrees' freedom chatter model is expanded and a three degrees freedom chatter model is established. The stability limit diagram of multi-axis milling is obtained by the taking spindle speed, axial cutting depth, and radial cutting depth as coordinate axes, so as to optimize process parameters and improve machining efficiency. Multi-axis milling is a complex machining process in that the workpiece structure is complex and the tool movement is changeable, which leads to the change of machining parameters. Thereby, the majority of research on multi-axis milling can mainly be divided into three categories, multi-axis milling force modeling, dynamic stability analysis, and the tool path plan. All these categories have been studied extensively. Firstly, [PERSON] et al. [3] presented a model for calculating the cutting force in five-axis milling, in which the cutting force coefficient was obtained by the transformation of the right angle. Additionally, the influence of the inclination angle was also included in the model. [PERSON] et al. [4] developed a cutting force calculation model through geometric analysis, in which the relationship between the un-deformed chip thickness, tool front angle, cutting speed, shear plane, and chip flow angle were described in detail. The three-dimensional cutting force in the tool coordinate system was obtained with the minimum energy method. On the basis of the Johnson-Cook theory, [PERSON] et al. [5] predicted the cutting force of the end milling cutter in the milling process with FEM. On the basis of the analysis of tool motion, [PERSON] et al. [6] proposed an instantaneous cutting force calculation model in the case of the existence of tool eccentricity in the process of five-axis milling. Recently, [PERSON] et al. [7] developed a cutting force model of multi-axis milling operations, which was calculated by using a new numerical method of the engagement boundary determination approach. [PERSON] et al. [8] used the zeroth-order approximation frequency domain method to predict the stability of the cutting tool by considering the process damping effect. A generalized cutting forces model of the five-axis milling process is proposed, which can be applied for an arbitrary mill geometry in multi-axis milling as well as three-axis milling and two-and-a-half-axis milling [9]. [PERSON] et al. [10] established a tool path generation method for the five-axis flank milling of pockets based on the constraints of cutting force and kinematics of machine tools. Most of the literature mentioned above assumed that the tool was rigid or neglected the weak rigidity of the workpiece. Meanwhile, the cutting force modeling has been relatively mature. However, most of the cutting force models are specific to certain tool forms only, such as end milling cutters, and ball milling cutters, and so forth, without considering the influence of time-varying parameters.