A geometric error tracing method based on the Monte Carlo theory of the five-axis gantry machining center

Abstract

This article proposes a tracing method to identify key geometric errors for a computer numerical control machine tool by cutting an S-shaped test piece. Adjacent part relationships and machine tool errors transform relationships are described by topology of the machining center. Global sensitivity analysis method based on quasi-Monte Carlo was used to analyze machining errors. Using this method, key geometric errors with significant influence on machining errors were obtained. Compensation of the key errors was used to experimentally improve machining errors for the S-shaped test piece. This method fundamentally determines the inherent connection and influence between geometric errors and machining errors. Key geometric errors that have great influence on machining errors can be determined quickly with this method. Thus, the proposed tracing method could provide effective guidance for the design and use of machine tools.

Introduction

High machining precision is critically important and has increased demand in recent years in aerospace, military, ships, cars, and other industries. Fabrication of high precision parts is heavily dependent of machining capabilities and capacity of the used tool. A variety of error sources related to the machine tool are responsible for machining errors. Therefore, to alleviate this problem, it is important to establish a mapping relationship between machine tool errors and machining errors. This would allow tracing machining errors caused by machine tool errors.

Geometric errors in machine tool parts contribute to spatial errors, which are a major factor affecting machining accuracy. ¹ Considerable effort has been made by researchers to identify each geometric error by direct measurements on the machine. In 1982, Bryan ² first reported the use of magnetic ball bars for measuring linear axis motion errors. Since then, various groups have applied this method to measure geometric errors of machines. Knapp and Matthias ³ proposed a circular test for identifying machine errors by measuring circular contour along different flat of machining space. This allowed identification of the main source of errors due to uncertainty of machine in three-dimensional space. Kakino et al. ⁴ carried out an in-depth study using the circular test with the aid of magnetic ball bars to obtain a relationship between machine motion errors and radial errors in polar plots. Tsutsumi and Saito ⁵ proposed an algorithm to identify eight specific geometric errors for a five-axis machine tool using the ball bar. Lei et al. ⁶ proposed a double ball bar test to measure rotary axes motion errors by simultaneously moving only two rotary axes while keeping all linear axes stationary. Ziegert and Mize ⁷ employed laser ball bars for measuring volumetric errors in a lathe. Mize and Ziegert ⁸ used the laser ball bar to measure all geometric errors for a three-axis machine in less than 30 min. Weikert ⁹ proposed a R-test measuring device for a five-axis machine tool, which has also been adopted by ISO 10791-6:2014. ¹⁰ This device can simultaneously detect displacement in three directions. As a result, this device provides more information in a single measurement. Ibaraki et al. ¹¹ presented an efficient calibration method to identify location errors and position-dependent geometric errors for a five-axis machine using R-test. Fan et al. ¹² described the use of a 6-degree-of-freedom measurement system based on the Doppler effect. This system had the capability to simultaneously measure all six motion errors of a linear stage. Knapp and Weikert ¹³ proposed a high-resolution system referred to as enhanced coordinate grid measurement based on Heidenhin grid encoder. This system, similar to the one described by Fan et al., allowed simultaneous measurements across 6 degrees of freedom for linear motion. They verified their test results using error compensation. Wang ¹⁴ presented a step-diagonal test method to identify geometric errors by measuring the diagonal of a processing space. Zhang et al. ¹⁵ proposed a method to determine straightness and angular errors for a linear. This method was based on the laser interferometry technique integrated with the shifting method. He et al. ¹⁶ presented a dual optical path measurement method to identify all six geometric errors associated with the rotational axis along the error-sensitive directions.

Non-cutting measurement methods do not reflect the influence of workpiece materials, temperature, tool performance, and cutting conditions. As a result, to accurately characterize the measurements, several studies have been reported in the literature, which are focused on cutting test to detect machine deviation for a specific workpiece. The NAS 979 test piece ¹⁷ proposed by NASA consists of a cone frustum, rhombus, taper, or other shapes. The NC-Gesellschaft (NCG) test piece ¹⁸ developed by European manufacturers blends various typical shapes together, such as conic, cylindrical boss, U-shaped surface, hemisphere surface, and developable ruled surface. A square-shaped step test piece ¹⁹ introduced by Kyoto University was machined at different index angles of the rotary axes. Relationship between machining errors for this test piece and geometric errors of a five-axis machine tool was formulated. A S-shaped test piece ²⁰ invented by Chengdu Aircraft Industrial is a twisted equal thickness non-developable ruled surface. It contains characteristics such as opening and closing angle conversion and a thin-walled curve. These features result in higher linkage performance and surface machining performance for the five-axis machine. As a result of the improved machining performance, use of this test piece has increased popularity in the aviation industry. Nevertheless, there is a serious lack of quantitative data for machining results of an S-shaped test piece. The XHAV2430×80 five-axis gantry machining center (FAGMC) produced by Beijing No. 1 Machine Tool Plant has been used to cut the S-shaped titanium alloy specimens successfully. This is a major technical breakthrough in titanium alloy cutting using machines manufactured in China. ²¹ However, machining errors and defects in any specimen including the S-shaped titanium alloy specimens cannot be completely eliminated. To understand and determine source of errors in the machine piece, it is important to find the key cause that affects machining errors.

Methods to evaluate geometric accuracy of a machine can be broadly categorized into two types. The first method employs direct measurement on the machine, whereas the second method category assesses machine precision by testing machining errors of the workpiece. Direct measurement methods carry out direct measurements of machine tool errors but rely heavily on measurement capability of the equipment used and identification methods used are generally complex. To overcome these drawbacks, in this study, a measurement method, which is simple and can quickly analyze geometric errors of a machine tool, is proposed. Here, machine tool errors are evaluated by detecting machining errors of S-shaped pieces. Sensitivity analysis of the machine tool error model is mainly based on local sensitivity analysis (LSA).^22–25 This method employs simple calculations, but has more steps related to the model analysis process. A significant limitation of LSA is that model sensitivity cannot be calculated for unknown nonlinear model input variables. Global sensitivity analysis (GSA) avoids this issue. ²⁶

In this article, geometric errors are considered in the kinematic equations according to the topology of FAGMC. GSA based on quasi-Monte Carlo (QMC) is introduced to analyze kinematic equations, which is the machining error model (MEM) of the machine. This allows sensitive errors of MEM to be estimated. These sensitive errors are key parameters that affect machining errors of the machine. Experiments were carried out to validate results. The presented method has generalized applicability and allows quick trace main factors influencing machining errors.

Machining accuracy model of the machine center

Machine tool structure and generalized coordinates

Figure 1 shows the schematic of the FAGMC structure. It can be seen from the figure that FAGMC consists of a bed, worktable, saddle, spindle, and swing head. Here, multi-body system theory is used to show the topological relationships for the geometric structure of FAGMC, shown in Figure 2. There are two branches in Figure 2, that is, workpiece branch (B₁-B₂-B_w) and tool branch (B₁-B₃-B₄-B₅-B₆-B_t). Then, kinematic equations are established in the form of homogeneous matrices.

Figure 1.

Schematic of FAGMC structure.

Figure 2.

Topology of FAGMC.

Suppose that the workpiece and tool are clamped without errors. Point O is origin of ${\mathrm{B}}_1$ body coordinate and also base origin of the machine. When $j=1,2,3,4,5,6$ , ${\mathrm{Q}}_{\mathrm{j}}^{\mathrm{l}}$ and ${\mathrm{Q}}_{\mathrm{j}}$ are ideal and actual motion reference points of ${\mathrm{B}}_{\mathrm{j}}$ , respectively. ${\mathrm{O}}_{\mathrm{j}}^{\mathrm{l}}$ and ${\mathrm{O}}_{\mathrm{j}}$ are ideal and actual position origin of ${\mathrm{B}}_{\mathrm{j}}$ body coordinate in ${\mathrm{Q}}_{\mathrm{j}}$ coordinate system. Thus, ${\mathbf{q}}_{\mathbf{j}}^1$ and ${\mathbf{q}}_{\mathbf{j}}$ are ideal and actual position vectors of ${\mathrm{B}}_{\mathrm{j}}$ relative to ${\mathrm{B}}_{\mathrm{j}-1}$ and ${\mathbf{q}}_{\mathbf{j}}^{\mathbf{e}}$ is position error vector. ${\mathbf{s}}_{\mathbf{j}}^{\mathbf{l}}$ and ${\mathbf{s}}_{\mathbf{j}}$ are ideal and actual motion vectors of ${\mathrm{B}}_{\mathrm{j}}$ relative to ${\mathrm{B}}_{\mathrm{j}-1}$ and ${\mathbf{s}}_{\mathbf{j}}^{\mathbf{e}}$ is displacement error. ${\mathbf{q}}_{\mathbf{w}}$ is position matrix of workpiece relative to the worktable and ${\mathbf{q}}_{\mathbf{t}}$ is position matrix of tool relative to B-axis. ${\mathbf{r}}_{\mathbf{w}}$ and ${\mathbf{r}}_{\mathbf{t}}$ are position vectors of tool center in the workpiece coordinate system and tool coordinate system.

A machine tool part generates 6-degree-of-freedom errors in space. Therefore, a five-axis machine tool with five parts would generate 30 geometric errors. Additionally, there are three vertical errors between three translational axes, two parallelism error components between C and Z axes, and another two parallelism error components between B and Y axes. So, for a five-axis machine, a total of 37 geometric errors would be generated. These errors are listed in Table 1.

Table 1.

Geometric errors for FAGMC.

Direction of error	Linear error			Angle error
	Along X	Along Y	Along Z	Around X	Around Y	Around Z
Translate along X-axis	${\delta }_x\left(x\right)$	${\delta }_y\left(x\right)$	${\delta }_z\left(x\right)$	${\epsilon }_x\left(x\right)$	${\epsilon }_y\left(x\right)$	${\epsilon }_z\left(x\right)$
Translate along Y-axis	${\delta }_x\left(y\right)$	${\delta }_y\left(y\right)$	${\delta }_z\left(y\right)$	${\epsilon }_x\left(y\right)$	${\epsilon }_y\left(y\right)$	${\epsilon }_z\left(y\right)$
Translate along Z-axis	${\delta }_x\left(z\right)$	${\delta }_y\left(z\right)$	${\delta }_z\left(z\right)$	${\epsilon }_x\left(z\right)$	${\epsilon }_y\left(z\right)$	${\epsilon }_z\left(z\right)$
Rotate around Z-axis	${\delta }_x\left(c\right)$	${\delta }_y\left(c\right)$	${\delta }_z\left(c\right)$	${\epsilon }_x\left(c\right)$	${\epsilon }_y\left(c\right)$	${\epsilon }_z\left(c\right)$
Rotate around Y-axis	${\delta }_x\left(b\right)$	${\delta }_y\left(b\right)$	${\delta }_z\left(b\right)$	${\epsilon }_x\left(b\right)$	${\epsilon }_y\left(b\right)$	${\epsilon }_z\left(b\right)$
Perpendicularity errors	${\epsilon }_{\mathrm{xy}}$ , ${\epsilon }_{\mathrm{xz}}$ , ${\epsilon }_{\mathrm{yz}}$ , ${\epsilon }_{\mathrm{xc}}$ , ${\epsilon }_{\mathrm{yc}}$ , ${\epsilon }_{\mathrm{xb}}$ , ${\epsilon }_{\mathrm{bz}}$

MEM of FAGMC

Figure 2 can be used to obtain location of the actual and theoretical tool center in the machine coordinate system

$${\mathbf{p}}_{\mathbf{t}}=\sum _{i=3}^6({\mathbf{q}}_{\mathbf{i}}^{\mathbf{l}}+{\mathbf{q}}_{\mathbf{i}}^{\mathbf{e}}+{\mathbf{s}}_{\mathbf{i}}^{\mathbf{i}}+{\mathbf{s}}_{\mathbf{i}}^{\mathbf{e}})+{\mathbf{q}}_{\mathbf{t}}+{\mathbf{r}}_{\mathbf{t}}$$

(1)

$${\mathbf{p}}_{\mathbf{w}}={\mathbf{q}}_2^1+{\mathbf{q}}_2^{\mathbf{e}}+{\mathbf{s}}_2^1+{\mathbf{s}}_2^{\mathbf{e}}+{\mathbf{q}}_{\mathbf{w}}+{\mathbf{r}}_{\mathbf{w}}$$

(2)

where ${\mathbf{p}}_{\mathbf{t}}$ is actual tool center position and ${\mathbf{p}}_{\mathbf{w}}$ is theoretical tool center position. Let $\mathbf{E}$ represents a unit matrix. Then, D-H matrices can be built, which are shown in Table 2.

Table 2.

The D-H matrices of FAGMC.

Part	Ideal position and motion matrices	Actual position and motion matrices
${\mathrm{B}}_2$	${\mathbf{q}}_2^{\mathbf{l}}\to {}_2{}^1{\mathbf{T}}_{\mathbf{p}}=\mathbf{E}$	${\mathbf{q}}_2^{\mathbf{e}}\to {}_2{}^1{\mathbf{T}}_{\mathbf{pe}}=\mathbf{E}$
	${\mathbf{s}}_2^{\mathbf{l}}\to {}_2{}^1{\mathbf{T}}_{\mathbf{s}}=\left[\begin{array}{cc}\mathbf{E}\hfill & \begin{array}{c}x\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$	${\mathbf{s}}_2^{\mathbf{e}}\to {}_2{}^1{\mathbf{T}}_{\mathbf{se}}=\left[\begin{array}{cc}\begin{array}{ccc}1\hfill & -{\epsilon }_z\left(x\right)\hfill & {\epsilon }_y\left(x\right)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z\left(x\right)}& 1\hfill & -{\epsilon }_x\left(x\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y\left(x\right)}& {\epsilon }_x\left(x\right)\hfill & 1\hfill \end{array}\hfill & \begin{array}{c}{\delta }_x\left(x\right)\hfill \\ \multicolumn{1}{c}{{\delta }_y\left(x\right)}\\ \multicolumn{1}{c}{{\delta }_z\left(x\right)}\end{array}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$
${\mathrm{B}}_3$	${\mathbf{q}}_3^{\mathbf{l}}\to {}_3{}^1{\mathbf{T}}_{\mathbf{p}}=\left[\begin{array}{cc}\mathbf{E}\hfill & \begin{array}{c}{q}_{3x}\hfill \\ \multicolumn{1}{c}{q_{3y}}\\ \multicolumn{1}{c}{q_{3z}}\end{array}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$	${\mathbf{q}}_3^{\mathbf{e}}\to {}_3{}^1{\mathbf{T}}_{\mathbf{pe}}=\left[\begin{array}{cc}\begin{array}{ccc}1\hfill & -{\epsilon }_{\mathrm{xy}}\hfill & 0\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{xy}}}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$
	${\mathbf{s}}_3^{\mathbf{l}}\to {}_3{}^1{\mathbf{T}}_{\mathbf{s}}=\left[\begin{array}{cc}\mathbf{E}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{y}\\ \multicolumn{1}{c}{0}\end{array}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$	${\mathbf{s}}_3^{\mathbf{e}}\to {}_3{}^1{\mathbf{T}}_{\mathbf{se}}=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z\left(y\right)\hfill & {\epsilon }_y\left(y\right)\hfill & {\delta }_x\left(y\right)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z\left(y\right)}& 1\hfill & -{\epsilon }_x\left(y\right)\hfill & {\delta }_y\left(y\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y\left(y\right)}& {\epsilon }_x\left(y\right)\hfill & 1\hfill & {\delta }_z\left(y\right)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$
${\mathrm{B}}_4$	${\mathbf{q}}_4^{\mathbf{l}}\to {}_4{}^3{\mathbf{T}}_{\mathbf{p}}=\mathbf{E}$	${\mathbf{q}}_4^{\mathbf{e}}\to {}_4{}^3{\mathbf{T}}_{\mathbf{pe}}=\left[\begin{array}{cccc}1\hfill & 0\hfill & {\epsilon }_{\mathrm{xz}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & -{\epsilon }_{\mathrm{yz}}\hfill & 0\hfill \\ \multicolumn{1}{c}{-{\epsilon }_{\mathrm{xz}}}& {\epsilon }_{\mathrm{yz}}\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$
	${\mathbf{s}}_4^{\mathbf{l}}\to {}_4{}^3{\mathbf{T}}_{\mathbf{s}}=\left[\begin{array}{cc}\mathbf{E}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{z}\end{array}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$	${\mathbf{s}}_4^{\mathbf{e}}\to {}_4{}^3{\mathbf{T}}_{\mathbf{se}}=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z\left(z\right)\hfill & {\epsilon }_y\left(z\right)\hfill & {\delta }_x\left(z\right)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z\left(z\right)}& 1\hfill & -{\epsilon }_x\left(z\right)\hfill & {\delta }_y\left(z\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y\left(z\right)}& {\epsilon }_x\left(z\right)\hfill & 1\hfill & {\delta }_z\left(z\right)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$
${\mathrm{B}}_5$	${\mathbf{q}}_5^{\mathbf{l}}\to {}_5{}^4{\mathbf{T}}_{\mathbf{p}}=\mathbf{E}$	${\mathbf{q}}_5^{\mathbf{e}}\to {}_5{}^4{\mathbf{T}}_{\mathbf{pe}}=\left[\begin{array}{cccc}1\hfill & 0\hfill & {\epsilon }_{\mathrm{xc}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & -{\epsilon }_{\mathrm{yc}}\hfill & 0\hfill \\ \multicolumn{1}{c}{-{\epsilon }_{\mathrm{xc}}}& {\epsilon }_{\mathrm{yc}}\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$
	${\mathbf{s}}_5^{\mathbf{l}}\to {}_5{}^4{\mathbf{T}}_{\mathbf{s}}=\left[\begin{array}{cccc}\cos c\hfill & -\sin c\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\sin c}& \cos c\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$	${\mathbf{s}}_5^{\mathbf{e}}\to {}_5{}^4{\mathbf{T}}_{\mathbf{se}}=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z\left(c\right)\hfill & {\epsilon }_y\left(c\right)\hfill & {\delta }_x\left(c\right)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z\left(c\right)}& 1\hfill & -{\epsilon }_x\left(c\right)\hfill & {\delta }_y\left(c\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y\left(c\right)}& {\epsilon }_x\left(c\right)\hfill & 1\hfill & {\delta }_z\left(c\right)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$
${\mathrm{B}}_6$	${\mathbf{q}}_6^{\mathbf{l}}\to {}_6{}^5{\mathbf{T}}_{\mathbf{p}}=\left[\begin{array}{cc}\mathbf{E}\hfill & \begin{array}{c}{q}_{6x}\hfill \\ \multicolumn{1}{c}{q_{6y}}\\ \multicolumn{1}{c}{q_{6z}}\end{array}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$	${\mathbf{q}}_6^{\mathbf{e}}\to {}_6{}^5{\mathbf{T}}_{\mathbf{pe}}=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_{\mathrm{bz}}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{bz}}}& 1\hfill & -{\epsilon }_{\mathrm{xb}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\epsilon }_{\mathrm{xb}}\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$
	${\mathbf{s}}_6^{\mathbf{l}}\to {}_6{}^5{\mathbf{T}}_{\mathbf{s}}=\left[\begin{array}{cccc}\cos b\hfill & 0\hfill & \sin b\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-\sin b}& 0\hfill & \cos b\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$	${\mathbf{s}}_6^{\mathbf{e}}\to {}_6{}^5{\mathbf{T}}_{\mathbf{se}}=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z\left(b\right)\hfill & {\epsilon }_y\left(b\right)\hfill & {\delta }_x\left(b\right)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z\left(b\right)}& 1\hfill & -{\epsilon }_x\left(b\right)\hfill & {\delta }_y\left(b\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y\left(b\right)}& {\epsilon }_x\left(b\right)\hfill & 1\hfill & {\delta }_z\left(b\right)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$
${\mathrm{B}}_{\mathrm{w}}$		${\mathbf{q}}_{\mathbf{w}}\to {\mathbf{T}}_{\mathbf{w}}\mathbf{=}\left[\begin{array}{cc}\mathbf{E}\hfill & \begin{array}{c}{q}_{\mathrm{wx}}\hfill \\ \multicolumn{1}{c}{q_{\mathrm{wy}}}\\ \multicolumn{1}{c}{q_{\mathrm{wz}}}\end{array}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$
${\mathrm{B}}_{\mathrm{t}}$		${\mathbf{q}}_{\mathbf{t}}\to {\mathbf{T}}_{\mathbf{t}}\mathbf{=}\left[\begin{array}{cc}\mathbf{E}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{z_w}\end{array}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$

According to Table 2, equations (1) and (2) can be rewritten into homogeneous matrices

$$p_t={\displaystyle \prod _{i=3}^6\left({}_i{}^{L^1\left(i\right)}{T}_p{}_i{}^{L^1\left(i\right)}{T}_{pe}{}_i{}^{L^1\left(i\right)}{T}_s{}_i{}^{L^1\left(i\right)}{T}_{se}\right)}{T}_t{r}_t$$

(3)

$${\mathbf{p}}_{\mathbf{w}}={}_2{}^1{\mathbf{T}}_{\mathbf{p}}{}_2{}^1{\mathbf{T}}_{\mathbf{pe}}{}_2{}^1{\mathbf{T}}_{\mathbf{s}}{}_2{}^1{\mathbf{T}}_{\mathbf{se}}{\mathbf{T}}_{\mathbf{w}}{\mathbf{r}}_{\mathbf{w}}$$

(4)

An important ideal condition is that the actual tool center location is equal to the theoretical tool center, that is, ${\mathbf{p}}_{\mathbf{t}}={\mathbf{p}}_{\mathbf{w}}$ . This is a necessary theoretical condition for precision machining. Then, let ${\mathbf{p}}_{\mathbf{tr}}=(\begin{array}{cccc}{x}_{\mathrm{tr}}\hfill & y_{\mathrm{tr}}\hfill & z_{\mathrm{tr}}\hfill & 1\hfill \end{array}{)}^T$ replace ${\mathbf{r}}_{\mathbf{w}}$ to represent actual tool center location in the workpiece coordinate system. Combining with equations (3) and (4), we can express ${\mathbf{p}}_{\mathbf{tr}}$ as

$$\begin{array}{c}{\mathbf{p}}_{\mathrm{tr}}={\left({}_2{}^1{\mathbf{T}}_{\mathbf{p}}{}_2{}^1{\mathbf{T}}_{\mathbf{pe}}{}_2{}^1{\mathbf{T}}_{\mathbf{s}}{}_2{}^1{\mathbf{T}}_{\mathbf{se}}{\mathbf{T}}_{\mathbf{w}}\right)}^{-1}\hfill \\ {\displaystyle \prod _{i=3}^6\left({}_i{}^{L^1\left(i\right)}{T}_p{}_i{}^{L^1\left(i\right)}{T}_{pe}{}_i{}^{L^1\left(i\right)}{T}_s{}_i{}^{L^1\left(i\right)}{T}_{se}\right)}{T}_t{r}_t\end{array}$$

(5)

Furthermore, to get the ideal tool center location ${\mathbf{p}}_{\mathbf{id}}=(\begin{array}{cccc}{x}_{\mathrm{id}}\hfill & y_{\mathrm{id}}\hfill & z_{\mathrm{id}}\hfill & 1\hfill \end{array}{)}^T$ in the workpiece coordinate, we assume all error parameters in equation (5) as zero

$${\mathbf{p}}_{\mathbf{id}}={\left({}_2{}^1{\mathbf{T}}_{\mathbf{p}}{}_2{}^1{\mathbf{T}}_{\mathbf{s}}{\mathbf{T}}_{\mathbf{w}}\right)}^{-1}{\displaystyle \prod _{i=3}^6\left({}_i{}^{L^1\left(i\right)}{T}_p{}_i{}^{L^1\left(i\right)}{T}_s\right)}{T}_t{r}_t$$

(6)

Therefore, MEM can be calculated from equations (5) and (6)

$${\mathbf{p}}_{\mathbf{e}}={\mathbf{p}}_{\mathbf{tr}}-{\mathbf{p}}_{\mathbf{id}}={\left(\begin{array}{cccc}{x}_e\hfill & y_e\hfill & z_e\hfill & 1\hfill \end{array}\right)}^T$$

(7)

where $x_e$ , $y_e$ , and $z_e$ represent MEM in X, Y, and Z directions, respectively. The space MEM of machine is given as

$$p_s=\sqrt{x_e^2+y_e^2+z_e^2}$$

In the case where error parameters of the machine are unknown, sensitivity of the model cannot be analyzed by the LSA method, and GSA method must be adopted.

GSA of the MEM

Monte Carlo and QMC method

The main difference between Monte Carlo (MC) and QMC methods is that the random numbers used in each of the methods is different. Since 1960s, mathematicians have been studying highly uniform quasi-random number sequences. These numbers replace random numbers in MC simulation and this method is called as the QMC method. Random numbers are mainly used to simulate randomness, while the use of quasi-random numbers is primarily to achieve uniformity of numbers. It has been proved by numerous theories that uniformity of random numbers is primarily responsible for accuracy of MC simulation results compared to randomness of the generated random numbers. ²⁷ Increase in uniformity of the random number distribution results in a decreased volatility of the random numbers, smaller deviation, higher calculation accuracy, and faster convergence rate for the calculations. Therefore, use of QMC method can improve simulation efficiency and reduce simulation errors. The uniformity of random numbers decrease with the increase in dimension. Since, FAGMC’s MEM involves up to 37 geometric errors, distribution of the 36th and 37th random numbers is observed. Figure 3(a) is distribution of the random numbers, and Figure 3(b) is distribution of quasi-random numbers. From the figure, it can be seen that uniformity of the quasi-random numbers is much better than the random numbers.

Figure 3.

(a) Distribution of random number (5000 points) and (b) distribution of quasi-random number (5000 points).

GSA based on QMC method

The objective of sensitivity analysis is to identify critical inputs (parameters and initial conditions) for a model and quantify the impact of input uncertainty on model outcomes. ²⁸ GSA based on MC is the best available method for calculating sensitive indices for model evaluation. ²⁶ Essential features of MC-based sensitivity analysis are repeated here for reader convenience and also because they are needed for discussion of the approximation algorithms used in this work.

According to Sobol’s ²⁹ definition, first-order sensitivity index is given as

$$S_i=\frac{V\left(E\left(y|x_i\right)\right)}{V\left(y\right)}$$

(8)

where $E(y|x_i)$ is the expected value of y over all-variables with $x_i$ fixed. Variance of $E(y|x_i)$ is a measure of the importance of input factor $x_i$ . $S_i$ measures contribution of input $x_i$ on output variance. A larger $S_i$ reflects greater influence on variance of y.

Saltelli ³⁰ proposed a GSA method based on QMC. In this case, inputs are assumed to be independent random variables and joint probability density function of the inputs is

$$p(x_1,x_2,\dots ,x_n)=\underset{i=1}{\overset{n}{\Pi }}{p}_i\left(x_i\right)$$

Thus, mathematical expectation and variance are obtained

$$E\left(y\right)=\int \int \dots \int f(x_1,x_2,\dots ,x_n)\underset{i=1}{\overset{n}{\Pi }}{p}_i\left(x_i\right)d{x}_i$$

$$V\left(y\right)=\int \int \dots \int f^2(x_1,x_2,\dots ,x_n)\underset{i=1}{\overset{n}{\Pi }}{p}_i\left(x_i\right)d{x}_i-E^2\left(y\right)$$

(9)

If $x_j(j=1,2,\dots ,n)$ is fixed to the value ${\tilde{x}}_j$

$$E\left(y|{\tilde{x}}_j\right)=\int \int \dots \int f(x_1,x_2,\dots ,{\tilde{x}}_j\dots x_n)\underset{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}{\overset{n}{\Pi }}{p}_i\left(x_i\right)d{x}_i$$

Thus, equation (9) can be written in the form

$$V\left(y|{\tilde{x}}_j\right)={\displaystyle \iint \dots {\displaystyle \int f^2\left(x_1,x_2,\dots ,{\tilde{x}}_j\dots x_n\right){\displaystyle \prod _{\begin{array}{l}i=1\\ i\ne j\end{array}}^n{p}_i\left(x_i\right)d{x}_i}}}-E^2\left(y|{\tilde{x}}_j\right)$$

where $V(y|{\tilde{x}}_j)$ and $E(y|{\tilde{x}}_j)$ are variance and expectation, respectively, for $x_j={\tilde{x}}_j$ . Expectation of $V(y|x_j)$ can be calculated as

$$\begin{array}{c}\multicolumn{1}{c}{E\left[V\left(y|x_j\right)\right]=\int V\left(y|{\tilde{x}}_j\right)p_j\left({\tilde{x}}_j\right)d{\tilde{x}}_j}\\ =\int [\int \int \dots \int f^2(x_1,x_2,\dots ,{\tilde{x}}_j\dots x_n)\underset{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}{\overset{n}{\Pi }}{p}_i\left(x_i\right)d{x}_i]p_j\left({\tilde{x}}_j\right)d{\tilde{x}}_j-\int E^2\left(y|{\tilde{x}}_j\right)p_j\left({\tilde{x}}_j\right)d{\tilde{x}}_j\hfill \\ =\int \int f^2(x_1,x_2,\dots ,x_n)\underset{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}{\overset{n}{\Pi }}{p}_i\left(x_i\right)d{x}_i-\int E^2\left(y|{\tilde{x}}_j\right)p_j\left({\tilde{x}}_j\right)d{\tilde{x}}_j\hfill \end{array}$$

(10)

Equation (9) minus equation (10) gives

$$V\left(y\right)-E\left[V\left(y|x_j\right)\right]=\int E^2\left(y|{\tilde{x}}_j\right)p_j\left({\tilde{x}}_j\right)d{\tilde{x}}_j-E^2\left(y\right)$$

(11)

Based on the research of Iman and Hora ³¹

$$V\left(y\right)=E\left[V\left(y|x_j\right)\right]+V\left[E\left(y|x_j\right)\right]$$

(12)

Substituting equation (11) into equation (12), we obtain

$$V\left[E\left(y|x_j\right)\right]=\int E^2\left(y|{\tilde{x}}_j\right)p_j\left({\tilde{x}}_j\right)d{\tilde{x}}_j-E^2\left(y\right)$$

(13)

To facilitate presentation and readability, we assume $U_j=\int E^2(y|{\tilde{x}}_j)p_j({\tilde{x}}_j)d{\tilde{x}}_j$ . In Saltelli, ³⁰ an efficient method was proposed to calculate U_j

$$\begin{array}{c}\multicolumn{1}{c}{U_j=\int {[\int \int \dots \int f(x_1,x_2,\dots {\tilde{x}}_j,\dots x_n)\underset{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}{\overset{n}{\Pi }}{p}_i\left(x_i\right)d{x}_i]}^2{p}_j\left({\tilde{x}}_j\right)d{\tilde{x}}_j}\\ =\int \int \dots \int f(x_1,x_2,\dots {\tilde{x}}_j,\dots x_n)f({x\prime }_1,{x\prime }_2,\dots {\tilde{x}}_j,\dots {x\prime }_n)\left(\underset{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}{\overset{n}{\Pi }}{p}_i\left(x_i\right)d{x}_i\right)\left(\underset{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}{\overset{n}{\Pi }}{p}_i\left({x\prime }_i\right)d{x\prime }_i\right)p_j\left({\tilde{x}}_j\right)d{\tilde{x}}_j\hfill \\ =\int \int \dots \int f(x_1,x_2,\dots x_n)f({x\prime }_1,{x\prime }_2,\dots x_j,\dots {x\prime }_n)\left(\underset{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}{\overset{n}{\Pi }}{p}_i\left(x_i\right)d{x}_i\right)\left(\underset{\genfrac{}{}{0ex}{}{i=1}{i\ne j}}{\overset{n}{\Pi }}{p}_i\left({x\prime }_i\right)d{x\prime }_i\right)\hfill \end{array}$$

(14)

Equation (14) can be expressed as

$$\begin{array}{c}F(x_1,x_2,\dots x_n,{x\prime }_1,{x\prime }_2,\dots {x\prime }_{j-1},{x\prime }_{j+1},\dots {x\prime }_n)\hfill \\ =f\times f^{*}\hfill \\ =f(x_1,x_2,\dots x_n)\times f({x\prime }_1,{x\prime }_2,\dots {x\prime }_{j-1},x_j,{x\prime }_{j+1},\dots {x\prime }_n)\hfill \end{array}$$

The function $F(x_1,x_2,\dots x_n,x{\prime }_1,x{\prime }_2,\dots x{\prime }_{j-1},x{\prime }_{j+1},\dots x{\prime }_n)$ consists of 2n − 1 independent variables. Consequently, mathematical expectations of f was estimated by sampling all variables N times. Thus, outcome of f is calculated from a $N\times n$ -dimension matrix A, which is constructed from the sampling data of input variables. Then, matrix A is resampled except for the j-row and is renamed as matrix B. This new matrix is used to estimate the outcome of $f^{*}$ . Matrices A and B are sampled according to the quasi-random number sequence of input variables in its domain

$$\mathbf{A}=\left[\begin{array}{cccc}\begin{array}{cc}{x}_{11}\hfill & x_{12}\hfill \end{array}\hfill & \cdots \hfill & \begin{array}{cc}{x}_{1j}\hfill & \cdots \hfill \end{array}\hfill & x_{1n}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}{x}_{21}\hfill & x_{22}\hfill \end{array}}& \cdots \hfill & \begin{array}{cc}{x}_{2j}\hfill & \cdots \hfill \end{array}\hfill & x_{2n}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}⋮\hfill & ⋮\hfill \end{array}}& \hfill & \begin{array}{cc}⋮\hfill & \hfill \end{array}\hfill & ⋮\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}{x}_{N1}\hfill & x_{N2}\hfill \end{array}}& \cdots \hfill & \begin{array}{cc}{x}_{\mathrm{Nj}}\hfill & \cdots \hfill \end{array}\hfill & x_{\mathrm{Nn}}\hfill \end{array}\right]\Rightarrow \left\{\begin{array}{c}{f}_1\hfill \\ \multicolumn{1}{c}{f_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{f_N}\end{array}\right\}$$

(15)

$$\mathbf{B}=\left[\begin{array}{cccc}\begin{array}{cc}{x\prime }_{11}\hfill & {x\prime }_{12}\hfill \end{array}\hfill & \cdots \hfill & \begin{array}{cc}{x}_{1j}\hfill & \cdots \hfill \end{array}\hfill & {x\prime }_{1n}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}{x\prime }_{21}\hfill & {x\prime }_{22}\hfill \end{array}}& \cdots \hfill & \begin{array}{cc}{x}_{2j}\hfill & \cdots \hfill \end{array}\hfill & {x\prime }_{2n}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}⋮\hfill & ⋮\hfill \end{array}}& \hfill & \begin{array}{cc}⋮\hfill & \hfill \end{array}\hfill & ⋮\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}{x\prime }_{N1}\hfill & {x\prime }_{N2}\hfill \end{array}}& \cdots \hfill & \begin{array}{cc}{x}_{\mathrm{Nj}}\hfill & \cdots \hfill \end{array}\hfill & {x\prime }_{\mathrm{Nn}}\hfill \end{array}\right]\Rightarrow \left\{\begin{array}{c}{f}_1^{*}\hfill \\ \multicolumn{1}{c}{f_2^{*}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{f_N^{*}}\end{array}\right\}$$

(16)

Substituting matrices A and B into the system model, we get $y_{\mathrm{iA}}=f(x_{i1},x_{i2},\dots ,x_{\mathrm{in}})$ , $y_{\mathrm{iB}}=f(x{\prime }_{i1},x{\prime }_{i2},\dots x_{\mathrm{ij}},\dots x{\prime }_{\mathrm{in}})$ . For the discrete variable $x_i$ , we can estimate

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{U}}_j\approx \frac{1}{N}\sum _{i=1}^{N}f(x_{i1},x_{i2},\dots ,x_{\mathrm{in}})f({x\prime }_{i1},{x\prime }_{i2},\dots ,x_{\mathrm{ij}},\dots ,{x\prime }_{\mathrm{in}})}\\ \multicolumn{1}{c}{\approx \frac{1}{N}\sum _{i=1}^N{y}_{\mathrm{iA}}{y}_{\mathrm{iB}}}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{E\left(y\right)\approx \frac{1}{N}\sum _{i=1}^{N}f(x_{i1},x_{i2},\dots x_{\mathrm{in}})}\\ \multicolumn{1}{c}{\approx \frac{1}{N}\sum _{i=1}^N{y}_{\mathrm{iA}}}\end{array}$$

$$\begin{array}{c}V\left(y\right)\approx \frac{1}{N}\sum _{i=1}^N{f}^2(x_{i1},x_{i2},\dots ,x_{\mathrm{in}})-{\left(\frac{1}{N}\sum _{i=1}^{N}f(x_{i1},x_{i2},\dots ,x_{\mathrm{in}})\right)}^2\hfill \\ \approx \frac{1}{N}\sum _{i=1}^N{y}_{\mathrm{iA}}^2-{\left(\frac{1}{N}\sum _{i=1}^N{y}_{\mathrm{iA}}\right)}^2\hfill \end{array}$$

(19)

Sensitivity index of $x_i$ can be approximated using equations (17)–(19)

$$\begin{array}{c}\multicolumn{1}{c}{S_i\approx \frac{V\left[E\left(y|x_j\right)\right]}{V\left(y\right)}}\\ \multicolumn{1}{c}{\approx \frac{\left(1/N\right)\sum _{i=1}^N{y}_{\mathrm{iA}}{y}_{\mathrm{iB}}-{\left(\frac{1}{N}\sum _{i=1}^N{y}_{\mathrm{iA}}\right)}^2}{\left(1/N\right)\sum _{i=1}^N{y}_{\mathrm{iA}}^2-{\left(\frac{1}{N}\sum _{i=1}^N{y}_{\mathrm{iA}}\right)}^2}}\end{array}$$

where $S_i$ is sensitivity of ith input variable of MEM. Then, estimation error is given as ²⁶

$$\mathrm{Ee}{r}_i\approx \frac{0.6745}{\sqrt{N}}\sqrt{\sum _{i=1}^N{\left(y_{\mathrm{iA}}{y}_{\mathrm{iB}}\right)}^2-{\left(\sum _{i=1}^N{y}_{\mathrm{iA}}{y}_{\mathrm{iB}}\right)}^2}$$

Sensitivity indices reflect influence of error parameters on MEM. In practice, if machine geometry error values are known, accuracy of the sensitivity analysis can be further improved. If sensitivity index of an error parameter is large but its value is very small compared to other parameters, it may not have a significant impact on the model. So, an actual machining errors sensitivity analysis method is proposed here, where geometric error values are used as weights for sensitivity analysis

$$S_{\mathrm{ai}}=S_i·E{r}_i$$

(21)

where $S_{\mathrm{ai}}$ represents actual machining error sensitivity indices, $S_i$ represents the sensitivity indices of the MEM, and $E{r}_i$ represents geometric error values. The flowchart illustrating the sequence of steps for machining error sensitivity analysis is shown in Figure 4.

Figure 4.

Sensitivity analysis diagram of FAGMC.

MEM studied here involves 19 geometric errors in the range [−0.1, 0.1]. To carry out the sensitivity analysis, MEM was input into the sensitivity analysis software shown in Figure 4 and constant values were assigned. A quasi-random sequence was generated and scrambled using the “MatousekAffineOwen” algorithm available in MATLAB. This allows generating a sequence with higher uniformity. In this way, matrices A and B were obtained and outputs of the matrices were calculated according to equations (15) and (16). Sensitivity indices for the MEM were calculated from equation (20). Finally, actual machining errors sensitivity indices were computed by equation (21).

Key geometric errors tracing for FAGMC

In this article, MEM of the Z-axis is the primary focus of research and discussions. Here, the key geometric errors affecting machining errors on the Z-axis are analyzed. For experiments, an aluminum alloy (material No. 7075-7451) was cut as an S-shaped test piece using FAGMC. Cutting parameters are listed in Table 3. Figure 5 shows an image of the S-shaped test piece being machined by the FAGMC.

Table 3.

The S-shape test piece cutting parameters.

Process no.	Process contents	Tools	Spindle speed (r/min)	Feed speed (mm/min)	Cutting allowance (mm)
1	Three-axis rough milling	Φ20 cylindrical milling cutter	4000	1000	1.9
2	Semi-finish milling	Φ20 cylindrical milling cutter	4000	1000	0.9
3	Five-axis linkage fine milling	Φ20 cylindrical milling cutter	4000	1000	0
4	Drilling	D15.5drill	500	100	0.1
5	Benchmark boring	D16 boring cutter	1000	100	0

Figure 5.

An S-shaped test piece being machined by FAGMC.

Machining error of the S-shaped test piece was determined using the Global performance 09.15.18 coordinate measuring machine (CMM) manufactured by Hexagon Company. Diameter of the probe used was 3 mm. Measurements were carried out on the surface of S-shaped test piece over three lines along the height of the test piece. These are labeled as L1, L2, and L3 and are shown in Figure 6. Here, L1 is 30 mm from the top face of S-shaped edge strip and the two adjacent lines are spaced 5 mm apart. To measure machining errors, 25 sample points were measured on each line. Figure 7 is the referring checkpoints schematic in the CMM system.

Figure 6.

Schematic illustrating measurement positions for error detection for the S-shaped test piece.

Figure 7.

Referring checkpoints schematic for the system of trilinear CMM.

Machining errors on Z-axis are shown in Figure 8. From Figure 8, it can be seen that a majority of the large errors are located in the middle of the S-shaped test piece. Therefore, an area near the central point was selected to carry out further analysis of geometric errors affecting this region. Next, NC code of the test piece was generated in UG NX8.5. This generated NC code near the central point is shown in Table 4.

Figure 8.

Machining errors on Z-axis for the S-shaped test piece.

Table 4.

NC code nearby the central point of the S-shaped test piece.

xco (mm)	yco (mm)	zco (mm)	c (°)	b (°)
86.23755	−27.20503	0.8644	233.3322	−5.36877

Substituting matrices from Table 2 into the equations (5)–(7) gives MEM on Z-axis

$$\begin{array}{c}{z}_e={\delta }_z\left(x\right)+{\delta }_z\left(y\right)+{\delta }_z\left(z\right)+{\delta }_z\left(c\right)+\cos (b){\delta }_z\left(b\right)-\sin \left(b\right){\delta }_x\left(b\right)-{\epsilon }_x\left(x\right)(y_{\mathrm{co}}+q_{3y})\hfill \\ -{\epsilon }_y\left(x\right)(x_{\mathrm{co}}-q_{3x})+q_{\mathrm{wx}}{\epsilon }_y\left(x\right)-q_{\mathrm{wy}}{\epsilon }_x\left(x\right)L\sin \left(b\right)\left(\begin{array}{c}{\epsilon }_y\left(b\right)+{\epsilon }_y\left(c\right)\hfill \\ \multicolumn{1}{c}{+\sin \left(c\right)\left(\begin{array}{c}{\epsilon }_x\left(x\right)-{\epsilon }_x\left(y\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_x\left(z\right)-{\epsilon }_{\mathrm{yz}}-{\epsilon }_{\mathrm{yc}}}\end{array}\right)}\\ \multicolumn{1}{c}{-\cos \left(c\right)\left(\begin{array}{c}{\epsilon }_y\left(x\right)-{\epsilon }_y\left(y\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y\left(z\right)+{\epsilon }_{\mathrm{xz}}-{\epsilon }_{\mathrm{xc}}}\end{array}\right)}\end{array}\right)\hfill \\ -q_{6x}\left(\begin{array}{c}{\epsilon }_y\left(c\right)+\sin \left(c\right)\left(\begin{array}{c}{\epsilon }_x\left(x\right)-{\epsilon }_x\left(y\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_x\left(z\right)-{\epsilon }_{\mathrm{yz}}-{\epsilon }_{\mathrm{yc}}}\end{array}\right)\hfill \\ \multicolumn{1}{c}{-\cos \left(c\right)\left(\begin{array}{c}{\epsilon }_y\left(x\right)-{\epsilon }_y\left(y\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y\left(z\right)+{\epsilon }_{\mathrm{xz}}-{\epsilon }_{\mathrm{xc}}}\end{array}\right)}\end{array}\right)+q_{6y}\left(\begin{array}{c}{\epsilon }_x\left(c\right)+\sin \left(c\right)\left(\begin{array}{c}-{\epsilon }_y\left(x\right)+{\epsilon }_y\left(y\right)\hfill \\ \multicolumn{1}{c}{+{\epsilon }_y\left(z\right)-{\epsilon }_{\mathrm{xz}}+{\epsilon }_{\mathrm{xc}}}\end{array}\right)\hfill \\ \multicolumn{1}{c}{-\cos \left(c\right)\left(\begin{array}{c}{\epsilon }_x\left(x\right)-{\epsilon }_x\left(y\right)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_x\left(z\right)-{\epsilon }_{\mathrm{yz}}-{\epsilon }_{\mathrm{yc}}}\end{array}\right)}\end{array}\right)\hfill \end{array}$$

(22)

There are 19 unknown geometric errors and 11 constants in equation (22). The physical meanings and values of the constants are listed in Table 5.

Table 5.

The constants of MEM on Z-axis.

Symbol	Value	Meaning
$x_{\mathrm{co}}$	86.24	X value of NC code
$y_{\mathrm{co}}$	−27.21	Y value of NC code
c	233.33	C-axis rotation angle of NC code
b	−5.37	B-axis rotation angle of NC code
$q_{\mathrm{wx}}$	3500	X direction distance from workpiece origin to machine origin
$q_{\mathrm{wy}}$	1350	Y direction distance from workpiece origin to machine origin
$q_{3x}$	3500	X direction distance from spindle end center to machine origin under initial condition
$q_{3y}$	0	Y direction distance from spindle end center to machine origin under initial condition
$q_{6x}$	0	X direction distance from rotation center of B-axis to spindle end center under initial condition
$q_{6y}$	0	Y direction distance from rotation center of B-axis to spindle end center under initial condition
L	75	Tool length

Sensitivity indices for the 19 geometric errors were obtained by substituting constant values from Table 5 and equation (22) into the sensitivity analysis software. Descending rank for geometric errors according to their sensitivity indices is shown in Figure 9. Maximum sensitivity index is the geometric error ${\epsilon }_y(x)$ , which is higher than the other values shown in the figure.

Figure 9.

Sensitivity indices of machining error model.

A total of 2000 sample points were taken from the range of geometric errors to plot the scatterplots (Figure 10), which reflect degree of association between the geometric errors and machining errors.

Figure 10.

Scatterplot showing correlation between geometric errors and machining error.

From Figure 10, it can be seen that ${\epsilon }_y(x)$ and $Z_e$ are negatively correlated. Correlations reduced with decrease in sensitivity indices. Smaller sensitivity parameters such as ${\epsilon }_x(c)$ , ${\delta }_z(b)$ , ${\delta }_x(b)$ , ${\delta }_z(x)$ , ${\delta }_z(y)$ , ${\delta }_z(z)$ , and ${\delta }_z(c)$ have no significant contribution to machining error $Z_e$ , as the sample points are uniform in the four quadrants. From Figures 9 and 10, correlation of the geometric errors to the machining error $Z_e$ was observed to decrease with decrease in sensitivity indices. Geometric error values and sensitivity indices were substituted into equation (21) followed by normalization of the results. Consequently, actual machining error sensitivity indices were obtained (Figure 11).

Figure 11.

Actual machining error sensitivity indices.

In Figure 11, geometric errors ${\epsilon }_y(x)$ , ${\epsilon }_y(y)$ , ${\epsilon }_x(x)$ , and ${\epsilon }_y(z)$ were found to be the most influential factors. Assuming that these four geometric errors are zero, absolute value of $Z_e$ reduces from 0.1985 to 0.0093 mm, that is, reducing by 95%. These results demonstrate that sensitivity errors largely affect MEM.

To further verify these results, we compensate the four sensitive geometric errors by modifying the NC code. The four sensitive geometric errors were retained and other error parameters were set as zero. According to the machining accuracy model, the mapping model $F_{\mathrm{NC}}$ and $F_{\mathrm{tp}}$ were deduced, where $F_{\mathrm{NC}}$ represents the mapping model of tool center point to NC code and $F_{\mathrm{tp}}$ represents the mapping model of NC code to tool center. The modified NC code was calculated iteratively, as shown in Figure 12. The iterations are terminated when the difference between tool center positions $(\begin{array}{ccc}{x}_{\mathrm{ti}}\hfill & y_{\mathrm{ti}}\hfill & z_{\mathrm{ti}}\hfill \end{array}{)}^T$ and ideal tool center positions $(\begin{array}{ccc}{x}_{\mathrm{id}}\hfill & y_{\mathrm{id}}\hfill & z_{\mathrm{id}}\hfill \end{array}{)}^T$ was less than the resolution of the machine movement axis $(\begin{array}{ccc}{\Delta }_x\hfill & {\Delta }_y\hfill & {\Delta }_z\hfill \end{array}{)}^T$ and $(\begin{array}{ccccc}{x}_{\mathrm{coi}}\hfill & y_{\mathrm{coi}}\hfill & z_{\mathrm{coi}}\hfill & c_i\hfill & b_i{)}^T\hfill \end{array}$ in the modified NC code.

Figure 12.

Modified NC code iteration algorithm.

Then, a new S-shaped test piece was cut on the FAGMC using the modified NC code, while other machining parameters and the workpiece clamping position remained unchanged. Next, machining errors were measured again on the new test piece. The results are shown in Figure 13.

Figure 13.

Contrast processing errors after error compensation: (a) machining errors in L1, (b) machining errors in L2, and (c) machining errors in L3.

After error compensation, machining errors on Z-axis were found to significantly decline. From Figure 13(c), absolute mean machining error on Z-axis reduced up to 61%. These experiment results demonstrate that sensitive geometric errors are key factors affecting the machining errors.

Discussion and results

Convergence rate and estimation error calculation

Taking sensitivity analysis results of ${\epsilon }_y(x)$ as an example, computational efficiency and estimation error of MC and QMC methods were compared. MC and QMC methods were used to analyze sensitivity of Z-axis MEM. Number of accurate bits and the calculation time were measured. Five sampling calculations for each sample was used to obtain five sensitivity indices. These five sensitivity indices were bitwise compared with each other from the highest digit. Exact same number of bits described here refers to the number of accurate bits, which is shown in Table 6.

Table 6.

Comparison of sensitivity calculation results of ${\epsilon }_y(x)$ .

Sample size	Calculating time (s)	MC		QMC
		Number of accurate bits	Eer	Number of accurate bits	Eer
${10}^4$	2	2	$2.1\times {10}^{-3}$	4	$2.1\times {10}^{-4}$
$5\times {10}^4$	9	3	$9.6\times {10}^{-4}$	5	$9.6\times {10}^{-5}$
${10}^5$	18	3	$6.7\times {10}^{-4}$	5	$6.8\times {10}^{-5}$
$5\times {10}^5$	90	3	$3.0\times {10}^{-4}$	5	$3.0\times {10}^{-5}$
${10}^6$	196	3	$2.2\times {10}^{-4}$	5	$2.1\times {10}^{-5}$
$5\times {10}^6$	999	4	$9.6\times {10}^{-5}$	6	$9.6\times {10}^{-6}$

Using the MC method to obtain four-bit accuracy numbers needs to be sampled $5\times {10}^6$ times and requires 999 s to complete. While using QMC method to obtain four-bit accuracy numbers only need ${10}^4$ times sampling and requires only 2 s to complete calculations. So, the efficiency of the QMC method is much higher than MC method. Since distribution of quasi-random numbers is more uniform, estimation error of the QMC method is about one-tenth of the MC method for the same sample size.

Sample size in the QMC method

Sample size has a significant impact on calculation speed, so it is an important parameter that needs to be taken into consideration. The sample size should be as small as possible to ensure convergence. We varied the sample size from 100 to 25,000 to calculate sensitivity indices, which is shown in Figure 14.

Figure 14.

Effect of sample size on sensitivity results.

According to Figure 14, sensitivity indices were observed to converge starting from a sample size of 5000. The values were found to fluctuate slightly on the order of a thousandth of magnitude leading to change in the order of approximate sensitivity indices. This order of the sensitivity indices is shown in Figure 15.

Figure 15.

Effect of sample size on sensitivity order.

Figure 15 shows that the sensitivity indices order can be divided into four separate layers. The orders were found to change within each layer for a sample size exceeding 400. Therefore, we could not determine a specific type of sensitivity index. However, only the order of the layer is required to judge key geometric errors. Based on the above discussion, 5000 was selected as the sample size for this study.

Results

Figure 11 shows that geometric errors ${\epsilon }_y(x)$ , ${\epsilon }_y(y)$ , ${\epsilon }_x(x)$ , and ${\epsilon }_y(z)$ are most influential geometric errors to machining errors of Z-axis. This result shows significant impact of angle errors of the translation axes on machining error of Z-axis. Figure 9 shows that MEM is sensitive to errors ${\epsilon }_y(x)$ , ${\epsilon }_x(x)$ , ${\epsilon }_y(b)$ , and ${\epsilon }_y(c)$ . While values of ${\epsilon }_y(b)$ and ${\epsilon }_y(c)$ were small, they ranked fifth and sixth in actual machining error analysis. Despite this, close attention to the machining model sensitivity factors is needed during manufacturing for the same type machine.

Figure 9 shows that sensitivities of all four perpendicularity errors ${\epsilon }_{\mathrm{yz}}$ , ${\epsilon }_{\mathrm{yc}}$ , ${\epsilon }_{\mathrm{xc}}$ , and ${\epsilon }_{\mathrm{xz}}$ are large. This proves that MEM is sensitive to perpendicularity errors. Figure 9 shows that there are no linear errors in the sensitivity errors. This confirms that MEM is more sensitive to angle errors.

To avoid duplication of results, only MEM on Z-axis was analyzed, that is, analysis process along the other two directions is the same. Servo system errors are also one of the critical factors affecting macroscopic-sized machining errors. After error compensation, machining errors found in the middle of the S-shaped were still large. This can, however, be attributed to the servo errors when the C-axis is changing direction. As a result, further investigation is necessary to carry out a comprehensive analysis of the servo error to trace machining error.

Conclusion

This article investigates methods for tracing key geometric errors affecting machining precision. Here, analysis was focused on significant source of geometrical errors for a large contour error area of the workpiece. For the MEM of FAGMC, sensitivity of inputs by QMC method was analyzed. Then, a sensitivity analysis method for actual machining errors with known geometric error values was proposed. In a trial cut test, mean absolute machining error was reduced by 56% after compensating for key geometric errors.

A major aim of this research was to propose a novel tracing methodology to determine importance of geometric factors affecting machining errors. Our results indicate that inputs and outputs of the MEM are intrinsically linked. GSA can be used to solve the relationship between the inputs and outputs. This method can be extended to practical applications in complex machine models. The results and the method presented in this article would be useful to machine tool manufacturers and users in designing and testing machines.