Load spectrum generation of machining center based on rainflow counting method

1. Introduction

A machining center (MC) is a high-tech machine tool and regarded as measuring a country’s comprehensive national strength and industrial modernization. In recent years, the gap between domestic (China) and foreign MCs in machining accuracy, machining speed, multi-axis movement and process complexity has been shortened gradually, while in regard to the reliability, there still exists considerable disparity, especially for those equipped with domestic functional units [1-3]. The inherent reliability of products is primarily determined by design, so the main method to improve the reliability is to conduct reliability design. A MC is a complex system consisting of mechanical, electrical, hydraulic, pneumatic subsystems and other subsystems. However, the probability design theory is not stable and mature enough for the complex system. So, conducting reliability tests become a primary approach to improve the reliability of products. Most kinds of potential failures can be stimulated by reliability bench tests and subsequently conduct reliability evaluation, failure analysis, and finally improve design. As shown in Fig. 1, it is a motorized spindle reliability bench test [1] and it can simulate the actual working conditions of MCs. As known from engineering experience, cutting loads of a MC is determined by cutting parameters, both of which are selected purpose-oriented by machine tool users. Under this circumstance, the working conditions of the field tests become uncontrollable. Meanwhile, because of the long testing period, it is difficult to update the testing conditions and the needs of products synchronically. Thus, reliability bench tests become the research hotspot. One problem for reliability bench tests is how to develop a high precision load spectrum.

Fig. 1Motorized spindle reliability bench test

The process of developing dynamic load spectrum for MCs includes selecting typical cutting processes, acquisition of load data, extrapolating load spectrum and compiling load spectrum. Load data are the bases of load spectrum. A few scholars have investigated about load acquisition for mechanical products and obtained lots of achievements in the automobile, construction machinery and aviation industries [4-6]. B. Oelmann [7] conducted 20 measurements annually to determine the driveline loads in the actual driving processes. D. S. Milčić [8] used measuring stripes LY 12 to obtain the torque on a working wheel shaft depending on the operating regime of the bucket wheel excavator. Jeongho Oh [9] conducted cluster analysis using a classification distribution to develop a load spectrum for statewide axles. D. N. Kolonius [10] simulated the movements of machine tools to predict bearing capacity and actual loads through a computer simulation program, but the extreme loads were not predicted for some parameters limitation. Jürgen Fleischer [11] simulated the motions of machine tools using the multi-body dynamics method to obtain the virtual load spectrum of machine tools. Then he developed a reliability evaluation model with considering the effects of machine tools loads.

MCs, which are the foundation of manufacturing equipment, are widely used in the mechanical industries. Load spectrum generation methods have been studied by various scholars to improve the reliability of machine tools. In the late 1990s, Shen Guixiang and Wang Yiqiang et al. [12] collected the cutting parameters through the field tests and established a static load spectrum through the cutting load empirical formulas. But the loads obtained by the empirical formulas are static, which ignored the dynamic loads generated during the actual cutting. According to the fatigue damage theory, dynamic loads are the main causes for the fatigue failures and random failures of products [13]. Thus, the static cutting load spectrum is no longer suitable for reliability bench tests of machine tools nowadays.

Considering the influence of the dynamic loads, we propose a dynamic load spectrum compilation method for MCs. In order to understand the actual load characteristics of MCs, the actual cutting experiments designed by the orthogonal experimental method are conducted in the laboratory. Then, cutting loads cycles are counted by the rainflow counting method to form a rainflow matrix after load signals preprocessing. The extrapolation of the cutting loads is carried out using the parametric extrapolation method, and then the parameters of the dynamic load spectrum are developed by the local-best particle swarm optimization (PSO) method. Finally, the two-dimensional load spectrum of MCs is compiled.

2. Cutting load cycles counting

Cutting loads of MCs are irregular and random in the actual cutting processes. Since, these loads change over time under the cutting conditions. To our knowledge, cycle counting is the most commonly used statistical counting method, in which load time history is measured through the accumulation of full-load cycles and half-load cycles.

2.1. Rainflow cycle counting theory

The rainflow counting algorithm is a two-parameter counting method developed by M. Matisuiski and Endo. This algorithm is used to extract load cycles from a load history, which can be obtained by measurement or simulation [15]. The rainflow cycle counting theory is described as follows:

1) Load-time history is assumed as a multistory roof turned clockwise by 90°, as shown in Fig. 2.

2) Water flows from its upper top on each of the “pagoda roofs” until one roof extends to the opposite side beyond the vertical of the starting point or the flow reaches a wet point.

3) When water reaches a peak, it drops to the next roof and stops until it meets another flow. Then, a full cycle is formed. Otherwise, a half cycle is generated.

4) Some rainflow cycles can be obtained according to the start and termination positions of flows. Then, the values of peaks and troughs of each cycle are extracted.

Fig. 2Rainflow cycle counting theory

2.2. Load cycle counting by the rainflow cycle counting theory

Data compression and cycle extraction from load signals are conducted before rainflow counting. Data compression is a process of converting the original load signals to an array of which elements are valid amplitudes by distinguishing the peaks and troughs of the original load signals. During the extraction of load cycles, only the first value of load signals is considered to be valid when two or more of the same values exist. Also, there is an extracting rule for the inflection point. If the value of the inflection point is greater than its adjacent values, then it is considered as the maximum value and needs to be extracted. In order to improve the efficiency, the four-point cycle counting method is applied. Firstly, $S_i$, $S_i+1$, $S_i+2$, $S_i+3$ are four consecutive points extracted from the load signals. Then, a full load cycle $∆S$ is defined as follows:

1

$∆S=\left|S_{i+2}-S_{i+1}\right|\le \min \left(\left|S_{i+3}-S_{i+2}\right|,\left|S_{i+1}-S_i\right|\right).$

After the full load cycles are defined, the two middle points, which are not peaks or troughs, need to be removed. Then two new points are selected, and these four new points are assessed by using Eq. (1). If assessment results are inconsistent with Eq. (1), the first picking point is removed, and the next new point is selected. This process is repeated until all points are processed. The remaining points after the cycle extraction compose the half cycles.

The rainflow counting algorithm is applied off-line by using MATLAB to generate the full and half cycles intuitively, as illustrated in Fig. 3.

Fig. 3Simulation of rainflow counting in MATLAB

Load cycles can be classified as up cycles and down cycles respectively according to the changing directions of inflection points. Fig. 3 shows that section 1 is an up full cycle, sections 3 and 4 are two down full cycles, sections 2 and 6 are two up half cycles, and sections 5 and 7 are two down half cycles.

Results of the rainflow cycle counting are generally stored in a matrix to save memory space. Each load time history corresponds to a rainflow matrix $R$ containing mean and amplitude information. The $R$ is expressed as follows:

2

$R=\left\{\begin{array}{c}{r}_{11}{r}_{12}\cdots r_{1j}\cdots r_{1n}\\ r_{21}{r}_{22}\cdots r_{2j}\cdots r_{1n}\\ ⋮⋮⋮⋮\\ r_{i1}{r}_{i2}\cdots r_{ij}\cdots r_{in}\\ ⋮⋮⋮⋮\\ r_{m1}{r}_{m2}\cdots r_{mj}\cdots r_{mn}\end{array}\right\},$

where $r_{ij}$ represents the number of load cycles under the $i$th level of the mean and $j$th level of the amplitude.

If a rainflow matrix is divided into $N$ columns and $M$ rows, the results can be presented in a matrix table, which is called a 2-D load spectrum (Table 1).

Table 1Cycle frequency of different means and amplitudes of loads

Cycle frequency		Amplitude
		1	2	$\cdots$	$j$	$\cdots$	$N$
Mean	1	$S_{11}$	$S_{12}$	$\cdots$	$S_{1j}$	$\cdots$	$S_{1N}$
	$⋮$	$⋮$	$⋮$		$⋮$		$⋮$
	$I$	$S_{i1}$	$S_{i2}$	$\cdots$	$S_{ij}$	$\cdots$	$S_{iN}$
	$⋮$	$⋮$	$⋮$		$⋮$		$⋮$
	$M$	$S_{M1}$	$S_{M2}$	$\cdots$	$S_{Mj}$	$\cdots$	$S_{MN}$

The load cycles versus the mean and amplitude can also be expressed by a histogram to provide a direct description (Figs. 4-6).

A large amount of loads that cause minimal or no damage to MCs is called invalid amplitude. Invalid amplitudes occupy a substantial space and would reduce the test efficiency. Thus, they should be eliminated. In this study, a dynamic threshold range model is selected to remove invalid amplitudes:

3

$∆\mathbf{F}=\left({\mathbf{F}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathbf{F}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\times ∆\%,$

where $∆\mathbf{F}$ is the threshold range, ${\mathbf{F}}_{max}$ is the maximum, ${\mathbf{F}}_{min}$ is the minimum of the loads, and $∆$ represents the precision of the dynamic threshold range model, which is usually 10 %.

Fig. 4Histogram of the load amplitudes

Fig. 5Histogram of the load means

Fig. 6Three dimensional joint histogram of the means and amplitudes

3. Load preprocessing

3.1. Decomposition of cutting force

Assuming that the cutting force of each cutter tooth is applied to a single point at the same time, so the resultant cutting force can be divided into the main cutting force $F_c$, passive force $F_p$, and feed force $F_f$. The decomposition of cutting force is shown in Figs. 7 and 8. The main cutting force is also called circular force; whose direction is the same as the main cutting movement. The torque $M_z$, which is applied to a spindle, is the product of the main cutting force and tool radius. The passive force is the spindle axial force that is parallel to the axis of the spindle.

Fig. 7Three dimensional graphics of a face milling cutter

Fig. 8Decomposition of milling force

Thus, the passive force has a significant influence in both of the machining accuracy and surface roughness. The feeding force is the spindle radial force, which is perpendicular to the axis of the spindle. The tool handle bends easily and vibrates severely when the feeding force is high. Therefore, the parameter of feeding force is the basic consideration to verify the stiffness of a feed system.

Considering the structure and installation of a dynamometer, the cutting force are divided into three mutual perpendicular directions ($X$, $Y$, and $Z$) in a machine coordinate system, which are the orthogonal force components $F_x$, $F_y$ and $F_z$ and they are shown in Fig. 8.

3.2. Regeneration of radial force in the time domain

The $F_x$, $F_y$, and $F_z$ are the cutting forces measured through the cutting tests, while $F_f$ and $F_p$ are the forces which are required to develop a load spectrum. Thus, Thus, an equivalent transformation method is applied to obtain the cutting force $F_f$ and $F_p$ from $F_x$, $F_y$, and $F_z$. According to the force interaction principle, $F_z$ is equal to $F_p$ and other force components satisfy the following equation:

4

$\left\{\begin{array}{l}{F}_f^2=F_x^2+F_y^2-F_c^2,\\ M_z=F_c\times \frac{D}{2}.\end{array}\right.$

So, the radial force $F_f$ can be calculated by Eq. (4), while axial force $F_p$ and cutting torque $M_z$ can be measured directly by a dynamometer. Figs. 9-11 illustrate the actual waveforms of the radial force, axial force, and cutting torque in the time domain under one typical working condition.

Fig. 9Actual waveform of the radial force in the time domain

Fig. 10Actual waveform of the axial force in the time domain

Fig. 11Actual waveform of the cutting torque in the time domain

3.3. Sectional processing of the cutting force

A complete cutting process includes three sections, which are processes of cutting in ($S_1$), cutting ($S_2$), and cutting out ($S_3$). Fig. 12 shows that the cutting force fluctuates seriously in different cutting processes. Although $S_1$ and $S_3$ contribute a small proportion of the entire process, the amplitude changes significantly in these two stages, particularly in $S_1$. Extreme loads are usually generated when the cutting tool comes in contact with a workpiece during $S_1$. Recent studies show that extreme loads would result in tool fracture and workpiece damages. But during $S_2$, the load amplitudes change smoothly without evidently increasing or decreasing. Results are different during the period $S_3$, the load amplitudes decrease at first but increase subsequently within a short time until the cutting tool finishes retraction.

Fig. 12Load waveform in different cutting processes

Limited by the size of the sample tested in laboratories, the proportion of the three cutting processes is inconsistent with the actual situations. For example, the proportion of $S_2$ obtained in the cutting tests is smaller than the actual one, thus, $S_2$ should be extrapolated on the basis of the actual cutting parameters. Statistical analysis reveals that $S_1$ and $S_3$ account for 2 % respectively and $S_2$ contributes 96 % of the whole cutting process. Then, the database of $S_2$ should be extended proportionally. Since $S_2$ is a static cutting process which consists of a series of equivalent small cutting period, the extension has a minimal effect on the whole loading distribution. Thus, in the case of constructing load spectrum, the factors above should be considered.

4. Parametric extrapolation method

5. Application example

Spindle speed $n$, feed rate $v_f$, cutting depth $a_p$, and cutting width $a_e$ are selected as experimental factors, and the experimental program is optimized by the orthogonal test method. A four-level table of the orthogonal test (Table 3) is established according to the cutting parameters listed in Table 2 which is obtained from the field tests.

Cutting tests were conducted in laboratory on the basis of the cutting parameters listed in Table 3. As shown in Fig. 13, the cutting force was measured by a Kistler dynamometer (Type: 9257B). Then the measured cutting load signals were input into a computer through the charge amplifier (Type: 5017B) and the data acquisition system (Type: 5697A1). Finally, the cutting force in the three directions were analyzed by the rainflow counting method and the results are shown in Figs. 14-18.

In order to find the cutting force distribution rules, Weibull, LogNorm and Gamma are selected as three alternative distribution models. We apply the PSO method to estimate parameters of these alternative distribution models. Thereafter, the distribution model is verified through the goodness-of-fit test method. Results are shown in Figs. 14-18 and the parameters for the three alternative distribution models are shown in Table 4.

Fig. 13Cutting load measuring system

Table 4 shows that the amplitude of the radial force is in accordance with the gamma distribution at $\alpha =$ 3.09, $\beta =$ 139.1. The mean of the radial force is consistent with the normal distribution at $\mu =$ 454.3, $\sigma =$ 149.2. The amplitude of the axial force is in line with Weibull distribution at $m=$ 335.8, $\eta =$ 2.16. The mean of the axial force is in agreement with the normal distribution at $\mu =$ 268.1, $\sigma =$ 96.3. The amplitude of the cutting torque is also consistent with Weibull distribution at $m=$ 36.4, $\eta =$ 2.07. The mean of the cutting torque is also in accordance with Weibull distribution at $m=$ 22.2, $\eta =$ 2.49.

Fig. 14Histogram of radial force amplitude value

Fig. 15Histogram of radial force mean value

Fig. 16Histogram of axial force amplitude value

Fig. 17Histogram of axial force mean value

Fig. 18Histogram of cutting torque amplitude value

Fig. 19Histogram of cutting torque mean value

According to Eq. (11), we obtain that ${\chi }^2=$ 213.5. Considering that:

So ${\chi }^2<{\chi }_{0.05}^2$. In conclusion, the cutting force mean, and amplitude are mutually independent when $a=$ 0.05, and the two-dimensional joint probability density functions of the radial force, axial force, and torque are expressed as follows:

For the purpose of facilitating the application of the load spectrum to the reliability bench tests and simulating the actual working conditions, a continuous probability density distribution is adapted into a two-dimensional program load spectrum. The 8-step ladder curve is used to carry out the loading test. The amplitude is divided into 8 levels by unequal interval (ratio coefficient: 1, 0.95, 0.85, 0.725, 0.575, 0.425, 0.275, and 0.125). Thus, a two-dimensional load spectrum of the cutting load amplitudes is formed. The two-dimensional (amplitude-mean) program load spectrum of the radial force is shown in Table 5. The cumulative reliability test time is 1100.52 h.

6. Conclusions

The load spectrum of machine tools is the foundation to conduct probability design and reliability bench tests. With the application of the rainflow cycle counting method, a new strategy of generating dynamic load spectrum is developed to resolve the inconformity between the field tests and the reliability bench tests because of the lack of the dynamic load information.

1) The cutting load measuring system is built and then the 16 sets of cutting tests are designed by the orthogonal test method.

2) In order to improve the precision of the load spectrum, the extrapolation of the load is carried out by the parametric extrapolation method. The results reveal that the probability distribution functions of the radial force obey to gamma distribution with $\alpha =$ 3.09, $\beta =$ 139.1. The mean of the radial force is consistent with the normal distribution at $\mu =$ 454.3, $\sigma =$ 149.2. The amplitude and mean of the axial force are in line with Weibull distribution at $m=$ 335.8, $\eta =$ 2.16 and the normal distribution at $\mu =$ 268.1, $\sigma =$ 96.3, respectively. The amplitude and mean of cutting torque are consistent with Weibull distribution at $m=$ 36.4, $\eta =$ 2.07 and $m=$ 22.2, $\eta =$ 2.49, respectively.

3) The joint distribution function of the mean and amplitude of the radial force, axial force, and torque are obtained by using a combination of statistical analysis methods, and the two-dimensional program loading spectrum of MCs is compiled.

The next step we will conduct reliability bench tests of the key function units of machine tools according to the two-dimensional program load spectrum. The generated equivalent load spectrum can be easily reproduced on a test bench for verifying and approving the reliability of MCs. It is of great significance to improve reliability level of MCs.