Compilation of drilling load spectrum based on the characteristics of drilling force

In order to overcome the problem that the existing methods of compiling load spectrum of spindle or machine tool mainly aim at the cutting force spectrum, torque spectrum and speed spectrum, respectively, which ignore the connection between each spectrum, in this paper, a method for compiling drilling load spectrum of motorized spindle in CNC machine tool based on the characteristics of drilling force is proposed. Firstly, drilling tests under different processing technologies are carried out to measure its load, and the correction coefficient in the empirical formula of drilling force is obtained through fitting the measured drilling force, which makes the calculation of the axial force and torque more reasonable. Secondly, compared with the extended factor method, the transcendental probability method is optimized to solve the ultimate load of the axial force. Then, after setting the axial force as the main load of drilling, an eight-stage load spectrum for the main load is compiled. Finally, according to the relationship between the axial force and other loads, the eight-stage loading spectrum is transformed into a multidimensional drilling load spectrum.


Introduction
As the basis of manufacturing industry, the reliability level of CNC machine tool is the key index to measure the comprehensive national strength [1][2][3]. The accuracy of key functional components of machine tool has a great impact on the accuracy of machine tool [4]. Motorized spindle is the key functional component of CNC machine tool, and its reliability level directly affects the reliability of CNC machine tool. The reliability test of motorized spindle is one of the key technologies to improve its reliability level. The reliability test can be divided into two types: field reliability test and bench reliability test. Due to the limitations of the actual machining condition, in order to ensure the production efficiency, the bench reliability test of motorized spindle becomes an effective test method. The key technology of bench reliability test is to simulate the actual machining process of the product, in order to trigger its real faults [5,6]. Therefore, the loading conditions of motorized spindle in the actual machining process need to be explored.
The load spectrum expresses the load time history of components in the real working condition, that is the statistical law of load and time, which is widely used in aerospace [7,8], automobile [9,10], mechanical engineering [11] and other fields. Compiling the load spectrum of motorized spindle provides the basis to simulate the actual working conditions in the laboratory. There are many working conditions of CNC machine tools, and the loads of machine tools vary greatly under different processing modes. Even if any process parameter changes under the same processing mode, the cutting load will change greatly [12]; thus, it is relatively complex to establish the load spectrum of spindle. Wang et al. [13] took the turning parameters collected at the processing site of First Automobile Works as the basic data for compiling a load spectrum, calculated the turning torque according to the empirical formula and established the functional relationship between the calculated torque and its frequency, which was regarded as a load spectrum for early machine tool spindles. Cheng [14] presented a continuous spectrum of fatigue load by studying the relationship between the power and its torque, which overcomes the complex problem in calculating equivalent fatigue load under discrete power and speed spectra. Wang et al. [15,16] used the multiregression method and the maximum likelihood method to estimate parameters of distribution models when compiling a load spectrum of machine tools and then optimized the load distribution model by a comprehensive evaluation method with a multifactorial model. At the same time, they first proposed the concept of a comprehensive load spectrum of NC lathe, established a joint distribution function between the torque and the speed of lathe spindles and compiled load spectra of lathe transmission components theoretically. For improving the typicality and versatility of a load spectrum of machine tools, Huang et al. [17] studied the utilization of speed and power in the processing site of typical users and then compiled a combined spectrum of speed and power. When compiling a relative torque spectrum of milling motorized spindles, Wang et al. [18] divided loads of the relative torque into five levels, optimized a distribution model by the K-S test method and compiled a five-stage load spectrum for its reliability test.
Recently, in view of the fact that the load data of lathe obeys several distributions such as Weibull, lognormal, gamma and beta, Chen et al. [19] proposed a method to optimize the load model, which fused various error information as an evaluation index and employed data envelopment analysis (DEA) method. At the same time, laws of load transmission of NC lathe were studied for obtaining the load model of each transmission part, which provided the theoretical basis for the design and optimization of transmission parts. In order to reduce the great difference between the load spectrum and its corresponding processing load, Chi and Li [20] subdivided the processing conditions into several parts and then established a torque program spectrum by connecting the divided parts which have been compiled as "low-high-low" loading blocks. Aiming at the situation that the load on the machine tools in the actual processing is a dynamic load, He et al. [21] carried out actual turning tests of typical machining parameter for reappearing dynamic loads in the lab, divided the cutting force angle within a small range at equal intervals, extrapolated loads in time domain or solved extreme values and finally compiled a three-dimensional program spectrum including rotational speed, cutting force and its angle, which provided a loading basis for reliability tests of NC lathe and its key functional components. In order to overcome the difficulty in solving the torque of milling, Zhu et al. [22] introduced milling force model to solve the milling torque and compiled a milling load spectrum including rotation speed, three direction force, torque and loading frequency, which provided a loading basis for simulating the actual milling condition.
However, it can be concluded that the existing methods of compiling load spectrum of spindle or machine tool mainly aim at the cutting force, torque and speed load of spindle of CNC machine tool and compile the cutting force spectrum, torque spectrum and speed spectrum, respectively, which leads to the lack of connection between each spectrum. Therefore, in this paper, according to the load characteristics of drilling conditions, the drilling spectrum compilation method is studied. Firstly, the drilling force obtained from the actual test is fitted by empirical formula, and the correction coefficient of drilling force is obtained, which provides a theoretical method for more accurate calculation of axial force and torque under all drilling processes. On the other hand, according to the characteristics of the measured drilling force signal, the drilling in and drilling out process is recorded as a load cycle. Secondly, the transcendental probability method is optimized to extrapolate the limit load of drilling axial force. Then, based on the most important load axial force in drilling process, the eight-level program spectrum is compiled. Finally, according to the relationship between axial force and other loads, the eight-level program spectrum is added to the drilling load spectrum including speed, torque and dynamic load frequency.
The remainder of this article is organized as follows. Section 2 subdivides drilling processes and elaborates the empirical formula for calculating the drilling force. Section 3 analyses the measured drilling signal and calculates the correction coefficient of the empirical formula. Section 4 compiles the drilling load spectrum. Finally, the conclusions are provided in Sect. 5.

The classification of drilling
Drilling technologies mainly include hole drilling, hole expansion, reaming, rounding and threading. Under different technologies, the characteristics of cutting force are different because of the great difference in cutting tools.
Hole drilling refers to one type of hole processing method that drills a hole on the workpiece which is prepared to be processed. And its axial force Fz and torque Mz can be calculated by empirical Eqs. (1) and (2) [26].
where d 0 is the diameter of drill bit, f is the feed rate and C F , z F , y F , C M , z M and y M are the related parameters, which are selected according to Table 1. Correction factor k M is equal to k F , which can be estimated according to the axial force measured by actual cutting tests. Hole expansion is one type of drilling method that enlarges the hole diameter from d 0 ′ to d 0 on the basis of existing holes. And its axial force and torque can be extrapolated by Eqs. (1) and (2) [25] as follows: In addition to hole drilling and hole expansion, drilling technologies include rounding, reaming and threading. Rounding is the process to machining rounding at the end of the hole, aiming at removing burrs caused by hole drilling and facilitating the subsequent assembly of shaft and holes. When a fillet is rounded, its axial force is generally far less than the maximum axial force of hole drilling due to its small cutting width and a large number of cutter teeth involved in cutting. Therefore, the load generated by rounding can be regarded as a nondestructive load in the process of compiling a drilling load spectrum.
Reaming is to increase a small amount of aperture based on an existing hole to improve its accuracy, surface roughness and verticality, which helps to meet the requirements of part drawings. Reaming is subordinate to semi-finishing or finishing, and its processing allowance is far less than that of hole drilling. Therefore, the load generated by reaming can also be regarded as a nondestructive load in the process of compiling a drilling load spectrum.
Threading is used to process internal threads. And its load is small because of a small processing allowance. Generally, special machine tools are employed to thread machining, and few machining centers are employed for it. Moreover, compared with hole drilling, the threading load is so small that it can be regarded as a nondestructive load in the process of compiling a drilling load spectrum. The specific processing effects of the five drilling technologies are shown in Fig. 1.

Collection of drilling parameters
The drilling parameters of a batch of similar NC machine tools are collected, and the loads under rounding, reaming and threading are all regarded as nondestructive loads and deleted. Therefore, the machining data of hole drilling and hole expansion are mainly recorded. The specific recording format of drilling process data is shown in Table 2, where the cutter's materials are high-speed steel (HSS).

Drilling signal processing
During the cutting force measurement test, firstly, the cutting force test system shall be built on the VDL-1000 vertical machining center in the laboratory, then the process data collected on site shall be reproduced on the machining center and finally the actual cutting force signal shall be measured and recorded. In the actual cutting state, the actual cutting force can be measured by dynamometer. The principle of cutting force test system is shown in Fig. 2. Firstly, the milling force signal is transformed into a weak electrical signal through the piezoelectric sensor in the dynamometer, and then the weak signal is amplified proportionally and without distortion through the charge amplifier. After that, the amplified signal is converted into a digital signal by the data acquisition system. Finally, the three-dimensional cutting force signal is displayed on the computer in real time after the conversion by the special software. The physical diagram of the measurement system is shown in Fig. 3. The drilling signals caused by hole drilling (a) and hole expansion (b) are plotted in Fig. 4. According to the analysis of the cutting force signals in Fig. 4, the drilling process is smooth, and the fluctuation of its axis force is very small. Therefore, the wavelet amplitude of the axis force can be regarded as nondestructive loads in the process of compiling the drilling load spectrum, that is the axial force of the drilling can be regarded as a constant value. Moreover, the values of F x and F y are so small that they can all be regarded as zero. In summary, the process from drilling in to drilling out can be regarded as a load cycle. Therefore, a square wave signal can be employed to simulate the drilling load in the reliability bench test of spindles.

Calculation of correction coefficient
Due to the fact that the 9129AA dynamometer can only measure three-dimensional force but can't directly measure torque, the empirical formula of drilling is employed to solve the axial force and torque. In order to solve their values precisely, the measured drilling force is fitted by the empirical formula to solve its correction coefficient. And  For calculating the correction factor, several drilling parameters collected on the spot are randomly selected to carry out cutting force tests, and the axial force measured by 9129AA dynamometer is shown in Table 3.
Taking k F as a variable, the optimization function of the drilling force is established as follows: where F(k F ) is the optimal function, F zmi the ith measured axial force, F zei (k F ) the ith calculated axial force according to the empirical formula and N the total number of measured axial force.
Taking the minimum of F(k F ) as the optimization objective, the correction coefficient k F is solved as 1.49 by the simulated annealing method. Through multiple adjustments and operations, the initial parameters of simulated annealing method are local learning factor 2.25, global learning factor 2.25, inertia weight coefficient 0.35, control parameter 3.85, cool down factor 0.98 and initial temperature 160. The iterative process of solving parameters by simulated annealing algorithm is shown in the Fig. 5. The errors of axial force between the calculated value and the measured one are shown in Table 3. And all errors are not so large that the empirical formula can be acceptable.

Compilation for drilling spectrum
In order to realize program control better, only the continuous cumulative frequency curve can be transformed into a stepped cumulative frequency curve. The stepped cumulative frequency curve can be divided into 4 stages, 8 stages, 16 stages and 32 stages. At present, the stepped cumulative frequency curve of eight-stage program is often used at home and abroad [27]. The main idea to compile a drilling spectrum is compiling an eight-stage program spectrum based on the main load first and then integrating other loads into the program spectrum according to their relationships.

Distribution model of axial force
Firstly, the frequency histogram of the axial force is drawn, and then an appropriate distribution function, Weibull function or kernel function, should be selected according to the frequency histogram. If the Weibull function is selected, its parameters need to be estimated by the maximum likelihood method. If the kernel function is selected, the data of axis force should be fitted with appropriate bandwidth. And then the axis force's probability density plot and cumulative probability density plot are drawn. At the same time, the point estimates of the axial force are plotted in the cumulative probability density plot to judge the fitting effect of the fitting function.
The point estimates of the cumulative probability density, which can be calculated by Eq. (6) [23].
where F zi is the i-th F z sorted from small to large and n is the total number of drilling load cycles.

Extrapolation of extreme load
The extreme load can be solved by extreme value theory (EVT), extended factor method and transcendental probability method. EVT mainly involves the statistical laws of the maximum or minimum value in random loads, which also includes the distribution of extreme loads. Due to the truth mentioned above that the drilling load under the same drilling parameters can be regarded as a constant value, EVT is not applicable to extrapolate the extreme load of drilling. Therefore, the extension factor method and the transcendental probability method are introduced to extrapolate the ultimate load.
According to the different types of extrapolation functions, the extended factor method [24] can be divided into two categories: straight line extrapolation and curve extrapolation. The two extrapolation methods have the same principle and similar effect. The concrete steps of the curve extrapolation are as follows: Step 1: The cumulative frequency distribution of the actual axial force is transformed into a standard cumulative distribution function, which is marked as an original curve shown in Eq. (7).
where F z is the axial force, F zmax the maximum axial force, H the cumulative frequency of axial force, H 0 the total cumulative frequency and n the total number of stages of load spectrum compiled.
Step 2: Translate the original curve upward for expanding the maximum cumulative frequency to 10 6 (H 1 ) and then generate the extrapolated curve shown in Fig. 6.
Step 3: Solve the extrapolated curve function. Firstly, the extrapolated curve is assumed as Eq. (8).
where a is the undetermined coefficient and F zlim the ultimate load of the axial force. Secondly, both sides of Eq. (8) should be taken logarithm: Thirdly, if ln(H) = ln(H 1 ) − ln (H 0 ), F z = F zmax , thus a can be calculated as shown in Eq. (10).
Finally, if in (H) is equal to 0 in Eq. (11), F zlim can be solved as follows.
The transcendental probability method [25] regards the maximum load in 10 6 load cycles as the ultimate load, and its calculation principle is shown in Fig. 7.
The curve in Fig. 7 is the probability density curve of the axial force. If P = 10 −6 , x LIM represents a load that occurs once in 10 6 load cycles, which is also regarded as an ultimate load. Usually, the ultimate load is solved according to the cumulative probability density function, which is shown in Eq. (13).
where F(x) is the cumulative probability density function of the axial force. Fig. 6 The schematic diagram of the extended factor method Fig. 7 Principle diagram of the transcendental probability method

Compilation of drilling load spectrum
The process for solving the rotational speed, torque and loading frequency associated with the axial force is as follows: Step 1: According to the upper and lower limits of each level of axial force, the rotational speed, torque and drilling time corresponding to the axial force are divided into eight groups, and the average values of rotational speed, torque and drilling time of each group are calculated, which are supplemented into the program spectrum.
Step 2: The lacked parts in Table 6 are secondary loads compared with the axial force, which can be divided into two categories: calculated value and statistical value. Torque is the calculated value, which can be obtained by enlarging or reducing the nearest torque proportionately. Rotational speed and drilling time are statistical values, which can be set as the value of the nearest series. For example, the rotational speed and drilling time listed in series 7 and 8 can be set as that of the nearest series, that is their value is the same as that of series 6, and the torque is proportionately amplified on the basis of that in series 6. If there is a linear correlation between the secondary loads and the main load, the linear relationship can be established to solve the secondary loads.
The proportion of load cycles P can be calculated by Eqs. (14) or (15).  where P i is the ratio of the number of cycles of the ith axial force, R i and R i+1 the lower and upper limits of the ith axial force respectively, f(x) the probability density function of the axial force and F(x) the cumulative probability density function of the axial force.

Acceleration factor of drilling load spectrum
When the spindle reliability bench test is carried out according to the drilling load spectrum mentioned above, the relationship between the running time of spindles on the test bench and that under actual working conditions can be described by acceleration factor: where p 1 and p 2 are the proportion of the working time of hole drilling and hole expansion to the total time of drilling respectively, k 1 and k 2 the ratio of cutting time to machine running time in their own process technology respectively, and a d the accelerating factor for eliminating the nondamage load which can be described as Eq. (16).
As shown in Eq. (16), the acceleration of the test can be achieved by deleting the low amplitude load. The key is to determine the limit boundary of the nondestructive load (excluding the strengthening load). At present, there are two general principles in academia. The first principle is regarding all loads less than half of the fatigue limit as non-destructive loads, but it only applies to parts that are easy to determine the fatigue limit. The second principle is setting the loads less than 0.1 ~ 0.125 times of the maximum amplitude as nondestructive loads. Although the second principle has high practicability and wide application range, the strengthening effect of low amplitude loads may be neglected. And the first principle is not applicable here because the spindle's fatigue limit can't be determined because of the complex situation. Therefore, the second principle is chosen with optimized parameters. In order to avoid excessive load deletion leading to neglect the strengthening characteristics of low-amplitude loads, the load less than 0.0625 times of the maximum amplitude is deleted as nondestructive loads in the paper, which considers the minimum proportion of eighth-stage program spectrum.

Case study
To verify the effectiveness of the proposed compilation method for drilling spectrum, in this section, the spectrum is compiled for a batch of drilling processes collected on the VDL or VDF series machining center of an engine parts manufacturing company in Guangxi. Firstly, the point estimates of the axial force are plotted in the cumulative probability density plot to judge the fitting effect of the fitting function, which is shown in Figs. 8 and 9.
Red asterisks in Figs. 8 and 9 refer to point estimates of the cumulative probability density. Comparing and analysing the above two graphs, the Gaussian kernel function has a better fitting effect, regardless of the overall fitting effect or local details. Therefore, the Gaussian kernel function is employed to fit the axial force.
The extended factor method and the transcendental probability method are compared to solve the extreme load of the axis force in Table 4.
As shown in Table 4, the extrapolation value of the transcendental probability method is bigger than that of the extended factor method, which is also in a reasonable range (17)  (less than 1.2 times). Therefore, the transcendental probability method is employed to extrapolate the extreme load of axial force of drilling. According to Conover's research results, the load spectrum can be divided into eight stages to fully reflect its fatigue characteristics. The proportions of magnitudes are 0.125, 0.275, 0.425, 0.575, 0.725, 0.85, 0.95 and 1, respectively. Thus, the program spectrum of axis force is shown in Table 5.
A spindle reliability test bench is shown in Fig. 10. In the test bench, the spindles rotate at a certain speed n, the loading mechanism with piezoelectric ceramics applies dynamic axial force to the test spindle at a certain frequency f and the loading spindle applies torque M z to the test spindle through the direct torque control. If the eight-level program spectrum is employed for loading, the loads such as rotational speed, torque and loading frequency are lacked. Therefore, it is necessary to supplement the program spectrum into a drilling load spectrum including rotational speed, torque and loading frequency.
According to the upper and lower limits of each level of axial force, the rotational speed, torque and drilling time corresponding to the axial force are divided into eight groups, and the average values of rotational speed, torque and drilling time of each group are calculated, which are supplemented into the program spectrum as shown in Table 6.
The data in Table 6 is incomplete because some series have no data when grouping (especially series 7 and 8). Thus, the rotational speed, torque and drilling time of these series need to be further studied. In summary, the complete drilling load spectrum is shown in Table 7. It has been proved in References 17 and 22 that the ratio of actual cutting time to the total running time of machine tools is approximately 0.4, so k 1 = k 2 = 0.4. At the same time, according to the collected data, p 1 = 0.529, p 2 = 0 and a d = 1.11; thus, the acceleration factor of the drilling load spectrum is up to 5.25.
If the load spectrum in Table 7 is applied to the reliability bench test of drilling spindle, in order to better simulate the alternation of loads in the actual working conditions, the main load is loaded by the way of "low-high-low", and the other loads change following the main load. The loading diagram is shown in Fig. 11.

Conclusions
Based on the analysis of the characteristics of drilling force, a multidimensional drilling load spectrum for reliability bench test is compiled with the idea of "seizing the main load and inducing other loads". The conclusions are as follows: 1. The drilling signal obtained from the cutting force test shows that the intermediate process of drilling is stable and its fluctuation value is small. When compiling a spectrum, the axial force in the intermediate process can be regarded as a constant value. Therefore, a single process from drilling in to drilling out can be regarded as a load cycle. 2. According to the analysis of drilling force signals, the time of the rising and falling edges of the axial force is much shorter than that of the whole drilling process. Therefore, a square wave load can be employed to simulate the drilling load in the reliability bench test. 3. By comparing the extended factor method with the transcendental probability method, their magnification factors are in the reasonable scope which is less than 1.2 times of the maximum value. And the magnification effect of the transcendental probability method is better than that of the extended factor method, and its theory is simpler and clearer. Therefore, the extreme load of the axial force is solved by the transcendental probability method. 4. A multidimensional drilling load spectrum for reliability bench test is compiled. It not only compresses the loadfree time during actual drilling time but also eliminates some low-amplitude loads with nondestructive effect, leading to a high acceleration factor about 5.25.
In the actual factory processing, the cutting tools of the machining center will wear along with the processing process. In the future research, the impact of tool wear on the load can be considered, so as to obtain a more accurate loading spectrum.
Author contribution Writing-original draft preparation, TJ; conceptualization, CY; formal analysis, JG; methodology, CC; writing-review and editing, DZ. All authors have read and agreed to the published version of the manuscript. Availability of data and material All the data presented and/or analysed in this study are available upon request to the corresponding author.

Declarations
Ethics approval Not applicable.

Conflict of interest
The authors declare no competing interests.