# Thermal-Induced Errors Prediction and Compensation for a Coordinate Boring  Machine Based on Time Series Analysis

Authors: Jun Yang; Dongsheng Zhang; Bin Feng; Xuesong Mei; Zhenbang Hu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

To improve the CNC machine tools precision, a thermal error modeling for the motorized spindle was proposed based on time series analysis, considering the length of cutting tools and thermal declined angles, and the real-time error compensation was implemented. A five-point method was applied to measure radial thermal declinations and axial expansion of the spindle with eddy current sensors, solving the problem that the three-point measurement cannot obtain the radial thermal angle errors. Then the stationarity of the thermal error sequences was determined by the Augmented Dickey-Fuller Test Algorithm, and the autocorrelation/partial autocorrelation function was applied to identify the model pattern. By combining both Yule-Walker equations and information criteria, the order and parameters of the models were solved effectively, which improved the prediction accuracy and generalization ability. The results indicated that the prediction accuracy of the time series model could reach up to 90%. In addition, the axial maximum error decreased from 39.6μm to 7 μm after error compensation, and the machining accuracy was improved by 89.7%. Moreover, the X
/
Y-direction accuracy can reach up to 77.4% and 86%, respectively, which demonstrated that the proposed methods of measurement, modeling, and compensation were effective.


---

## Body

## 1. Introduction

The precision CNC coordinate boring machine is a tool for processing complex box-type components. Thermal error will account for a larger proportion of total error as the machine tools become more sophisticated. However, the accuracy decreases and becomes far lower than the initial design value after the machine is used for a long period of time. This decreased accuracy over time primarily results from inadequate maintenance and accuracy stability, and the thermal error is the main factor for the inadequate accuracy, accounting for 70% of the total number of errors arising from various error sources [1]. Donmez et al. proposed that changing temperatures produce thermal errors and the thermal error is a major factor for reducing the machine precision [2]. And the motorized spindle has more complicated dynamic nonstationary and speed-dependent thermal characteristics than conventional spindles [3]. A nonuniform temperature distribution causes thermal errors in CNC machine tools, and this distribution becomes nonlinear and nonstationary and varies with time. Most problems in thermal error research field focus on how to exactly measure the thermal characteristics, establish thermal-induced error model with accurate prediction and excellent generalization, and efficiently compensate the thermal deformation to improve the machining accuracy.Firstly, the method of accurately measuring the temperature distribution and thermal error is the most important issue. Vissiere et al. measured the spindle geometric error with a new method in which the measurement accuracy can reach even the nanometer [4]. Vyroubal presented a method focused on compensation of machine’s thermal deformation in spindle axis direction based on decomposition analysis, which is a cheap and effective strategy [5]. Hong and Ibaraki studied thermal characteristics of a rotary axis on the five-axis machine and analyzed effect of thermal error on errors motion of the rotary axis [6]. Wang et al. collected the spindle thermal deformation in three directions with three-point measurement and described a method to group the data of thermal sensors [7]. Huang et al. measured the spindle thermal drifts with five eddy current displacement sensors and proposed a method combining the back propagation (BP) and genetic algorithm (GA) to model the thermal expansion in axial direction [8]. But the radial declined angle errors were ignored, and the thermal angle errors are key factors to compensate the spindle terminal processing accuracy.Secondly, the thermal error model must be built accurately based on the experiment measurement. Yang et al. [9] used artificial neural networks (ANNS) to establish a relationship between temperature and the thermal error of a spindle. EL Ouafi et al. constructed an artificial neural network model for spindle thermal errors with the temperature drawing on statistical methodology, which effectively improves the machining accuracy [10]. The support vector machine is a new machine learning theory, which has many advantages such as simple algorithm, global optimization, versatile, and strong generalization ability. The scholars Lin et al. and Zhao et al. established a spindle thermal error model based on least square support vector machine theory, and the model has perfect robustness [11, 12].The mentioned researchers have established the axial thermal elongation model, but they overlooked the thermal yaw and pitch errors of the spindle.In recent years, the finite element method (FEM) is applied to analyze temperature fields and the thermal deformation of machine tools. Creighton et al. used the finite element method to analyze temperature distribution characteristics for a high-speed micromilling spindle and constructed an exponential model of the axial thermal error considering the spindle speed and running time [13]. Haitao et al. proposed a method for calculating thermal conductivity coefficient of the spindle surface and simulated and analyzed variation principles of the temperature field and thermal deformation of the spindle [14].However, the precision CNC machine tool error is a mutual coupling of many complex factors that are affected by many variables, and therefore, it is extremely difficult to establish a theoretical equation based on the perspective of thermoelasticity and heat transfer.If the measurement and modeling of the thermal errors are completed, the next work could be the compensation. Fu et al. and Miao et al. built the spindle axial thermal error model by applying multivariate linear regression method [15, 16]; Wang and Yang also proposed a prediction model for the axial thermal deformation and applied the model to compensate error on a CNC machine [17]; Liu et al. compensated the thermal drift on the milling and boring machine in the Z-direction [18]. Ouafi et al. presented an integrated comprehensive modeling approach for thermal errors real-time compensation based on multiple temperature measurements; after compensation, the spindle thermal-induced errors were reduced from 19 μm to less than 1 μm [19]. There are other scholars who investigated the spindle axial thermal error compensation method [20, 21], and the machining accuracy was improved effectively. Gebhardt et al. described a high precision grey-box thermal error model for compensation on five-axis machines and the thermally induced errors of the rotary/swivelling can be reduced up to 85% [22]. Pajor and Zapłata presented a set allowing for supervising the feed screw thermal elongation to reduce ball screw thermal errors [23]. Zhang et al. developed a novel compensation implementation technique for machine tools based on the function of the machine external zero point shift and Ethernet data communication protocol, which improved the machine precision [24].The existing literatures were mainly about measurement and modeling of the spindle axial thermal elongation, but the radial thermal angle errors were ignored. The spindle thermal deformation of a CNC machine tool is usually expressed as the actual spatial position and gesture may deviate from the theoretical value; namely, the running spindle would induce the drifts of geometric dimensions and spatial phase, thus affecting the precision of machine processing. Of course, the deviated drifts must include axial elongation and radial thermal pitch and yaw angle errors, which reduce the machining accuracy. For example, the paper’s object is a coordinate boring machine, and its spindle axial thermal expansion may affect the bore geometric size. Meanwhile, the radial thermal inclination angle errors could influence the geometric dimensions and surface roughness of the hole, so it is extremely necessary to measure axial and radial thermal errors simultaneously. In order to realize the error compensation, the spindle radial thermal angle errors must be translated into components of the linear coordinate axis. Moreover, the thermal error offset components are closely related to both radial thermal inclination angle errors and the handle length. So the radial thermal angle errors cannot be ignored. The three-point method could only measure the absolute thermal deformation in the single radial direction, and it cannot obtain the thermal inclination angle error. Therefore, the thermal error compensation model deduced by the method was not accurate enough. This indicated that the three-point method cannot completely reflect the variation of the spindle radial thermal deformation, and the radial thermal error compensation model is not accurate based on this measurement.Based on the above analysis, the five-point method is applied, and the measurement contains the following advantage: (1) five-point method could simultaneously measure axial and radial thermal drifts of the motorized spindle system, so the variation of the spindle position and orientation could be analyzed, solving the challenge that the three-point method cannot measure the inclined angles of the spindle radial thermal errors; (2) the integrated thermal error modeling will provide a more accurate mathematical equation for the error compensation, which includes the spindle elongation, radial thermal pitch and yaw angle errors based on the five-point method.The time series analysis provides a set of approaches to process dynamic data. The method primary meaning is that all types of data are approximately described by mathematical models. Through the model analysis, the data internal structure can be mined. So we can forecast its trends and make necessary control on it. Wang et al. applied time series analysis method to establish a spindle thermal error model and compensated errors [25–27], and they acquired better results.This paper focuses on the spindle system of a box-type precision CNC coordinate boring machine. The thermal error modeling and compensation equations were proposed. Then thermal balance experiments were performed by using the five-point method. The time series models were established for spindle axial thermal elongation and radial thermal yaw and pitch errors based on the experiments. After that, we translated the thermal drifts into coordinate offsets and established the final compensation mathematical equations in three directions. Finally, the error compensation was carried out. The results show that the machine precision can be improved efficiently.


## 2. Thermal Error Compensation Principle

The motorized spindle system of CNC machine tools could produce plenty of heat during the processing, the accumulated heat leads to the spindle thermal deformation and deviation to the relative position between cutting tools and the workpiece; finally, it will reduce the accuracy of the terminal processing. For the coordinate boring machine in this paper, the thermal deformation of the spindle may affect the geometric size and surface roughness of the workpiece. Generally, thermal drifts occur in three directions. In this paper, the axial thermal elongation is defined asE. And in the radial directions, the thermal yaw angle in X
O
Z plane is defined as θ
x, and the thermal pitch angle in Y
O
Z plane is defined as θ
y.Figure1 describes the spatial pose of the spindle thermal drift on X
O
Z, and the point P is the deflexion center. After the spindle experienced axial elongation E and radial inclination θ
x, the spindle declined from P
O
⃑ to P
O
′
⃑, so the offset component in X-direction is as follows:
(1)
Δ
O
x
=
(
D
0
x
+
D
+
Δ
D
)
sin
θ
x
,
where the offset in X-direction is Δ
O
x, D
0
x is the distance between the deflexion center and the spindle nose, D is the length of the cutting tool, and Δ
D is the axial elongation E.The geometric principle of the spindle thermal error compensation.
(a)
(b)The compensation offset inZ-direction is Δ
O
z:
(2)
Δ
O
z
=
Δ
D
-
Δ
O
D
=
Δ
D
-
(
D
0
x
+
D
+
Δ
D
)
(
1
-
cos
⁡
θ
x
)
.
Because the axial elongation is less than the length of the tool, that is,
(3)
Δ
D
≪
D
0
x
+
D
,
θ
x
⟶
0
,
so
(4)
sin
θ
x
⟶
θ
x
,
cos
⁡
θ
x
⟶
1
.
Equations (3)-(4) are substituted into (1)-(2); then thermal error compensation component in X- and Z-directions can be obtained, respectively:
(5)
Δ
O
x
=
(
D
0
x
+
D
)
θ
x
,
(6)
Δ
O
z
=
Δ
D
.
This indicates that the offset in Z-direction has no relationship with the tool length, while X-directional compensation offset is closely related to that.Similarly, the thermal error offsetΔ
O
y in the Y-direction can be obtained:
(7)
Δ
O
y
=
(
D
0
y
+
D
)
θ
y
,
where D
0
y is the distance between the deflexion center and the spindle nose.


## 3. Time Series Analysis of Thermal Errors

In order to solve the final thermal error compensation components in the axes, the comprehensive model of the axial elongationE and the radial thermal angle errors θ
x and θ
y must be derived. Time series analysis could be used to accomplish this task. The basic idea of the time series analysis is that a mathematical model, which accurately reflects the system dynamic dependency, is established through the analysis of the time sequence samples based on a limited sample of the observation system, and it is applied to predict and monitor the future behavior of the system.


### 3.1. Thermal Errors Stationarity Judgment and Gauss Standardization

Given enactment that the sequence{
X
t
,
t
=
0,1
,
2
,
…
} is a discrete stochastic process, the Augmented Dickey-Fuller Test (ADF) determines the stationarity of the time series. Assuming that three time series of the spindle system thermal errors are E
t, θ
x
,
t, θ
y
,
t, if the results exhibit that the thermal error sequences are nonstationary, the sequences are normalized by Gaussian standardization so that the sequences have a smoothness. Of course, after Gaussian standardization the stationarity of the new series is still determined by ADF, if the new series are still nonstationary, the difference method or other function transformation approaches should be utilized to deal with the time series until they are stationary:
(8)
μ
=
1
n
∑
i
=
1
n
x
i
,
σ
2
=
1
n
-
1
∑
i
=
1
n
(
x
i
-
μ
)
2
.In the formula,μ
i and σ
i
2 are mean and variance of axial thermal elongation and radial thermal declination angle errors.The original sequences of the spindle thermal errors are taking a standardizing processing:(9)
E
t
~
N
(
μ
1
,
σ
1
2
)
⟹
E
t
′
=
E
t
-
μ
1
σ
1
⟹
E
t
′
~
N
(
0,1
)
,
θ
x
t
~
N
(
μ
2
,
σ
2
2
)
⟹
θ
x
t
′
=
θ
x
t
-
μ
2
σ
2
⟹
θ
x
t
′
~
N
(
0,1
)
,
θ
y
t
~
N
(
μ
3
,
σ
3
2
)
⟹
θ
y
t
′
=
θ
y
t
-
μ
3
σ
3
⟹
θ
y
t
′
~
N
(
0,1
)
,
where E
t
′, θ
x
t
′, θ
y
t
′ are new time series of spindle thermal errors and the new sequences approximate standard normal distribution.


### 3.2. Box-Jenkins Model Identification

Using the autocorrelation function (ACF) and partial autocorrelation function (PACF) to identify the thermal error series pattern, if the calculated results of ACF and PACF have tailing, it indicates that new time series of the standardized spindle thermal errors are the autoregressive and moving average hybrid models ARMA(p
,
q); the model is as follows [28]:
(10)
X
t
=
ϕ
1
X
t
-
1
+
⋯
+
ϕ
p
X
t
-
p
+
ε
t
-
θ
1
ε
t
-
1
-
⋯
-
θ
q
ε
t
-
q
,
where ϕ
i
(
i
=
1,2
,
…
,
p
) is the autoregressive parameter and θ
j
(
j
=
1,2
,
…
,
q
) is the moving average parameter. Define the later operator as B:
(11)
B
X
t
=
X
t
-
1
,
B
p
X
t
=
X
t
-
p
.
Thus,
(12)
θ
(
B
)
=
1
-
θ
1
B
-
θ
2
B
2
-
⋯
θ
q
B
q
,
ϕ
(
B
)
=
1
-
ϕ
1
B
-
ϕ
2
B
2
-
⋯
ϕ
p
B
p
.
The ARMA model is transformed into
(13)
ϕ
(
B
)
X
t
=
θ
(
B
)
ε
t
.


### 3.3. The Parameter Estimation

George et al. [28] suggested that if the ARMA(p
,
q) model contains p order autoregressive AR(p) and q order moving average MA(q), its autocorrelation function is a pattern mixed exponential and attenuation sine wave after p-q order delay. Correspondingly, the partial correlation function is not exact exponential form, but it is controlled by a mixture of index and decaying sine wave. And the covariance is
(14)
γ
t
+
k
,
t
=
Cov
(
X
t
+
k
,
X
t
)
=
E
[
(
X
t
+
k
-
μ
t
+
k
)
(
X
t
-
μ
t
)
]
=
γ
k
,
k
=
(
0
,
±
1
,
±
2
,
…
,
n
-
1
)
.Autocorrelation function is(15)
ρ
k
=
γ
k
γ
0
.According to the statistical theory, the covariance function of time series with stationarity and zero mean is estimated as follows:(16)
γ
^
k
=
1
N
∑
t
=
k
+
1
N
X
t
X
t
-
k
.
Thus, the autocorrelation function is estimated as follows:
(17)
ρ
^
k
=
γ
^
k
γ
^
0
=
∑
t
=
k
+
1
N
X
t
X
t
-
k
∑
t
=
1
N
X
t
2
.Yule-Walker equations can obtain autoregression coefficients; setk
=
1,2
,
…
,
p, so the linear equations are gotten; namely,
(18)
ρ
1
=
ϕ
1
+
ϕ
2
ρ
1
+
⋯
+
ϕ
p
ρ
p
-
1
ρ
2
=
ϕ
1
ρ
1
+
ϕ
2
+
⋯
+
ϕ
p
ρ
p
-
2
h
h
⋮
hhhhhhh
⋮
hhhhhhh
⋮
ρ
p
=
ϕ
1
ρ
p
-
1
+
ϕ
2
ρ
p
-
2
+
⋯
+
ϕ
p
.
Replacing theoretical autocorrelation ρ
k with the estimated autocorrelation ρ
^
k, the autoregression coefficients of Yule-Walker estimation can be obtained. Define the vectors as follows:
(19)
ϕ
=
(
ϕ
1
ϕ
2
⋮
ϕ
p
)
,
ρ
p
=
(
ρ
1
ρ
2
⋮
ρ
p
)
,
P
p
=
(
1
ρ
1
ρ
2
⋯
ρ
p
-
1
ρ
1
1
ρ
1
⋯
ρ
p
-
2
⋮
⋮
⋮
⋮
⋮
ρ
p
-
1
ρ
p
-
2
ρ
p
-
3
⋯
1
)
.
The parameter ϕ in (18) can be written as
(20)
ϕ
=
P
p
-
1
ρ
p
.Markw
t
′
=
ϕ
(
B
)
X
t, and deal with the ARMA model as a moving average process as follows:
(21)
w
t
′
=
θ
(
B
)
ε
t
.
Covariance γ
j
′ of w
t
′ can be presented by covariance γ
j of X
t, and j
=
0,1
,
…
,
q:
(22)
γ
j
′
=
∑
i
=
0
p
ϕ
i
2
γ
j
+
∑
i
=
1
p
(
ϕ
0
ϕ
i
+
ϕ
1
ϕ
i
+
1
+
⋯
+
ϕ
p
-
i
ϕ
p
)
(
γ
j
+
i
+
γ
j
-
i
)
.
Convention ϕ
0
=
-
1; the covariance function of MA(q) process is
(23)
γ
0
′
=
(
1
+
θ
1
2
+
⋯
+
θ
q
2
)
σ
ε
2
,
γ
k
′
=
(
-
θ
k
+
θ
1
θ
k
+
1
+
⋯
+
θ
q
-
k
θ
q
)
σ
ε
2
.
Estimations of the parameters σ
ε
2, θ
q, θ
q
-
1, …, θ
1 can be calculated under q order moving average, using iteration as follows:
(24)
σ
ε
2
=
γ
0
′
1
+
θ
1
2
+
⋯
+
θ
q
2
,
θ
j
=
-
(
γ
j
′
σ
ε
2
-
θ
1
θ
j
+
1
-
θ
2
θ
j
+
2
-
⋯
-
θ
q
-
j
θ
q
)
.
And promise θ
0
=
0, θ
q, θ
q
-
1, …, θ
1 can be made zero also at the beginning of iterations.


## 3.1. Thermal Errors Stationarity Judgment and Gauss Standardization

Given enactment that the sequence{
X
t
,
t
=
0,1
,
2
,
…
} is a discrete stochastic process, the Augmented Dickey-Fuller Test (ADF) determines the stationarity of the time series. Assuming that three time series of the spindle system thermal errors are E
t, θ
x
,
t, θ
y
,
t, if the results exhibit that the thermal error sequences are nonstationary, the sequences are normalized by Gaussian standardization so that the sequences have a smoothness. Of course, after Gaussian standardization the stationarity of the new series is still determined by ADF, if the new series are still nonstationary, the difference method or other function transformation approaches should be utilized to deal with the time series until they are stationary:
(8)
μ
=
1
n
∑
i
=
1
n
x
i
,
σ
2
=
1
n
-
1
∑
i
=
1
n
(
x
i
-
μ
)
2
.In the formula,μ
i and σ
i
2 are mean and variance of axial thermal elongation and radial thermal declination angle errors.The original sequences of the spindle thermal errors are taking a standardizing processing:(9)
E
t
~
N
(
μ
1
,
σ
1
2
)
⟹
E
t
′
=
E
t
-
μ
1
σ
1
⟹
E
t
′
~
N
(
0,1
)
,
θ
x
t
~
N
(
μ
2
,
σ
2
2
)
⟹
θ
x
t
′
=
θ
x
t
-
μ
2
σ
2
⟹
θ
x
t
′
~
N
(
0,1
)
,
θ
y
t
~
N
(
μ
3
,
σ
3
2
)
⟹
θ
y
t
′
=
θ
y
t
-
μ
3
σ
3
⟹
θ
y
t
′
~
N
(
0,1
)
,
where E
t
′, θ
x
t
′, θ
y
t
′ are new time series of spindle thermal errors and the new sequences approximate standard normal distribution.


## 3.2. Box-Jenkins Model Identification

Using the autocorrelation function (ACF) and partial autocorrelation function (PACF) to identify the thermal error series pattern, if the calculated results of ACF and PACF have tailing, it indicates that new time series of the standardized spindle thermal errors are the autoregressive and moving average hybrid models ARMA(p
,
q); the model is as follows [28]:
(10)
X
t
=
ϕ
1
X
t
-
1
+
⋯
+
ϕ
p
X
t
-
p
+
ε
t
-
θ
1
ε
t
-
1
-
⋯
-
θ
q
ε
t
-
q
,
where ϕ
i
(
i
=
1,2
,
…
,
p
) is the autoregressive parameter and θ
j
(
j
=
1,2
,
…
,
q
) is the moving average parameter. Define the later operator as B:
(11)
B
X
t
=
X
t
-
1
,
B
p
X
t
=
X
t
-
p
.
Thus,
(12)
θ
(
B
)
=
1
-
θ
1
B
-
θ
2
B
2
-
⋯
θ
q
B
q
,
ϕ
(
B
)
=
1
-
ϕ
1
B
-
ϕ
2
B
2
-
⋯
ϕ
p
B
p
.
The ARMA model is transformed into
(13)
ϕ
(
B
)
X
t
=
θ
(
B
)
ε
t
.


## 3.3. The Parameter Estimation

George et al. [28] suggested that if the ARMA(p
,
q) model contains p order autoregressive AR(p) and q order moving average MA(q), its autocorrelation function is a pattern mixed exponential and attenuation sine wave after p-q order delay. Correspondingly, the partial correlation function is not exact exponential form, but it is controlled by a mixture of index and decaying sine wave. And the covariance is
(14)
γ
t
+
k
,
t
=
Cov
(
X
t
+
k
,
X
t
)
=
E
[
(
X
t
+
k
-
μ
t
+
k
)
(
X
t
-
μ
t
)
]
=
γ
k
,
k
=
(
0
,
±
1
,
±
2
,
…
,
n
-
1
)
.Autocorrelation function is(15)
ρ
k
=
γ
k
γ
0
.According to the statistical theory, the covariance function of time series with stationarity and zero mean is estimated as follows:(16)
γ
^
k
=
1
N
∑
t
=
k
+
1
N
X
t
X
t
-
k
.
Thus, the autocorrelation function is estimated as follows:
(17)
ρ
^
k
=
γ
^
k
γ
^
0
=
∑
t
=
k
+
1
N
X
t
X
t
-
k
∑
t
=
1
N
X
t
2
.Yule-Walker equations can obtain autoregression coefficients; setk
=
1,2
,
…
,
p, so the linear equations are gotten; namely,
(18)
ρ
1
=
ϕ
1
+
ϕ
2
ρ
1
+
⋯
+
ϕ
p
ρ
p
-
1
ρ
2
=
ϕ
1
ρ
1
+
ϕ
2
+
⋯
+
ϕ
p
ρ
p
-
2
h
h
⋮
hhhhhhh
⋮
hhhhhhh
⋮
ρ
p
=
ϕ
1
ρ
p
-
1
+
ϕ
2
ρ
p
-
2
+
⋯
+
ϕ
p
.
Replacing theoretical autocorrelation ρ
k with the estimated autocorrelation ρ
^
k, the autoregression coefficients of Yule-Walker estimation can be obtained. Define the vectors as follows:
(19)
ϕ
=
(
ϕ
1
ϕ
2
⋮
ϕ
p
)
,
ρ
p
=
(
ρ
1
ρ
2
⋮
ρ
p
)
,
P
p
=
(
1
ρ
1
ρ
2
⋯
ρ
p
-
1
ρ
1
1
ρ
1
⋯
ρ
p
-
2
⋮
⋮
⋮
⋮
⋮
ρ
p
-
1
ρ
p
-
2
ρ
p
-
3
⋯
1
)
.
The parameter ϕ in (18) can be written as
(20)
ϕ
=
P
p
-
1
ρ
p
.Markw
t
′
=
ϕ
(
B
)
X
t, and deal with the ARMA model as a moving average process as follows:
(21)
w
t
′
=
θ
(
B
)
ε
t
.
Covariance γ
j
′ of w
t
′ can be presented by covariance γ
j of X
t, and j
=
0,1
,
…
,
q:
(22)
γ
j
′
=
∑
i
=
0
p
ϕ
i
2
γ
j
+
∑
i
=
1
p
(
ϕ
0
ϕ
i
+
ϕ
1
ϕ
i
+
1
+
⋯
+
ϕ
p
-
i
ϕ
p
)
(
γ
j
+
i
+
γ
j
-
i
)
.
Convention ϕ
0
=
-
1; the covariance function of MA(q) process is
(23)
γ
0
′
=
(
1
+
θ
1
2
+
⋯
+
θ
q
2
)
σ
ε
2
,
γ
k
′
=
(
-
θ
k
+
θ
1
θ
k
+
1
+
⋯
+
θ
q
-
k
θ
q
)
σ
ε
2
.
Estimations of the parameters σ
ε
2, θ
q, θ
q
-
1, …, θ
1 can be calculated under q order moving average, using iteration as follows:
(24)
σ
ε
2
=
γ
0
′
1
+
θ
1
2
+
⋯
+
θ
q
2
,
θ
j
=
-
(
γ
j
′
σ
ε
2
-
θ
1
θ
j
+
1
-
θ
2
θ
j
+
2
-
⋯
-
θ
q
-
j
θ
q
)
.
And promise θ
0
=
0, θ
q, θ
q
-
1, …, θ
1 can be made zero also at the beginning of iterations.


## 4. Thermal Characterization Experiment

In order to establish the comprehensive thermal error model of the axial elongation and radial thermal angle errors of the spindle system with the application of time series analysis, the spindle system of a precision CNC coordinate boring machine was chosen as the research object, and the thermal equilibrium experiments were carried on. Then the temperature field and thermal deformation of the spindle system were analyzed.


### 4.1. Experimental Setup

The experimental system is shown in Figure2, which focuses on the spindle of the precision CNC coordinate boring machine. The measuring equipment and functions are as follows: a synchronous acquisition system is used to determine the temperatures and thermal drifts. This system usesPt100 precision magnetic temperature sensors to measure the spindle system temperatures. High-precision eddy-current sensors are applied to measure the spindle thermal drifts. Temperature sensors locations are as follows: front bearing (T6, T7), rear bearing (T1), the motor (T8, T11), ambient temperature (T5), spindle base (T2), the cooling fluid inlet (T9), bearing cooling out (T3), front bearing coolant out (T4), and the motor cooling out (T10).Figure 2
Experimental setup.


### 4.2. Measurement Principle

The spindle thermal drifts are measured by using the five-point method [29], the diagram of displacement sensors measurement is shown in Figure 3. The spindle is parallel to z-axis, and the axial thermal expansion can be obtained by the displacement sensor S
5. The radial thermal yaw θ
x partial X-direction is measured by the S
1 and S
3, and the radial thermal pitch θ
y partial Y-direction is measured by the S
2 and S
4.Figure 3
Spindle five-point installation diagram.After the spindle running for a long period, the thermal elongation expanded to axial direction and thermal angle inclined to radial direction, resulting from the uneven temperature gradient distribution, which is shown in Figure4, and the thermal yaw angle θ
x is
(25)
Δ
L
3
=
L
3
i
-
L
3
0
,
Δ
L
1
=
L
1
i
-
L
1
0
,
Δ
L
=
Δ
L
3
-
Δ
L
1
,
tan
θ
x
=
Δ
L
D
,
where i denotes the number of measurements. The thermal yaw angle is too small in this experiment, that is, θ
x
→
0, so
(26)
θ
x
~
tan
θ
x
.
As shown in (27), the thermal yaw angle can be obtained by applying (25):
(27)
θ
x
=
(
L
3
i
-
L
1
i
)
-
(
L
3
0
-
L
1
0
)
D
,
where L
3
0 and L
1
0 are the radial displacements between the sensor probes and the spindle measured by S
3 and S
1, respectively, in the initial state and L
3
i and L
1
i are the transient displacements during the running operation. D is the distance between S
1 and S
3, S
2 and S
4, and D = 120 mm.Figure 4
The spindle thermal inclination sketch.Similarly, the thermal pitch angle in theY-direction can be obtained:
(28)
θ
y
=
(
L
4
i
-
L
2
i
)
-
(
L
4
0
-
L
2
0
)
D
.


### 4.3. Results and Analysis

The spindle speed is a major factor affecting the thermal characteristics. The temperature gradient and the heat generated by the CNC machine spindle are different at different speeds, thus generating different thermal drifts. So the spindle speed affects the temperature field distribution and the magnitude of thermal errors. Meanwhile, the thermal error model is closely related to the intrinsic behavior of thermal drift data. Therefore, a proper test condition shall be proposed so that the thermal error model put forward in this paper can be applied to the practical engineering field and be used to improve machining accuracy. In order to create such kind of testing condition, the distribution of spindle speeds in this paper is set to simulate some common processing conditions of the boring machine. The specific speeds distribution is shown in Figure5.Figure 5
Step speeds distribution.Thermal characteristics of CNC machine tools cover two main aspects of the temperature distribution and thermal error, and the nonuniform distribution of the temperature field gradient is the direct cause of spindle thermal drifts. So for a more comprehensive study of the spindle thermal properties, a simple time-domain analysis of the temperature field is also conducted in the paper. The spindle system temperature variations are shown in Figure6. The overall trends of temperatures on all measuring points increase with time. The rear bearing has the highest temperature reaching 30.4°C due to large capacity, heavy load, and severe friction, which generates more heat, and the following is the motor whose temperature is 27.3°C.Figure 6
Temperatures of the spindle.Figure7 presents the spindle thermal drifts. The time until equilibrium reached is approximately 385 min, with a maximum elongation 39.6 μm. The thermal error on x-axis direction is positive, and its thermal yaw angle is θ
x in the X
O
Z plane; the maximum amount of hot offset error is 35 μm. Thermal error in Y-direction is negative which indicates that, during operation, the spindle is closer to the displacement sensors S
2
/
S
4, and its thermal pitch angle is θ
y in the Y
O
Z plane; the maximum thermal offset is of 20.2 μm.Figure 7
Thermal drifts of the spindle.


## 4.1. Experimental Setup

The experimental system is shown in Figure2, which focuses on the spindle of the precision CNC coordinate boring machine. The measuring equipment and functions are as follows: a synchronous acquisition system is used to determine the temperatures and thermal drifts. This system usesPt100 precision magnetic temperature sensors to measure the spindle system temperatures. High-precision eddy-current sensors are applied to measure the spindle thermal drifts. Temperature sensors locations are as follows: front bearing (T6, T7), rear bearing (T1), the motor (T8, T11), ambient temperature (T5), spindle base (T2), the cooling fluid inlet (T9), bearing cooling out (T3), front bearing coolant out (T4), and the motor cooling out (T10).Figure 2
Experimental setup.


## 4.2. Measurement Principle

The spindle thermal drifts are measured by using the five-point method [29], the diagram of displacement sensors measurement is shown in Figure 3. The spindle is parallel to z-axis, and the axial thermal expansion can be obtained by the displacement sensor S
5. The radial thermal yaw θ
x partial X-direction is measured by the S
1 and S
3, and the radial thermal pitch θ
y partial Y-direction is measured by the S
2 and S
4.Figure 3
Spindle five-point installation diagram.After the spindle running for a long period, the thermal elongation expanded to axial direction and thermal angle inclined to radial direction, resulting from the uneven temperature gradient distribution, which is shown in Figure4, and the thermal yaw angle θ
x is
(25)
Δ
L
3
=
L
3
i
-
L
3
0
,
Δ
L
1
=
L
1
i
-
L
1
0
,
Δ
L
=
Δ
L
3
-
Δ
L
1
,
tan
θ
x
=
Δ
L
D
,
where i denotes the number of measurements. The thermal yaw angle is too small in this experiment, that is, θ
x
→
0, so
(26)
θ
x
~
tan
θ
x
.
As shown in (27), the thermal yaw angle can be obtained by applying (25):
(27)
θ
x
=
(
L
3
i
-
L
1
i
)
-
(
L
3
0
-
L
1
0
)
D
,
where L
3
0 and L
1
0 are the radial displacements between the sensor probes and the spindle measured by S
3 and S
1, respectively, in the initial state and L
3
i and L
1
i are the transient displacements during the running operation. D is the distance between S
1 and S
3, S
2 and S
4, and D = 120 mm.Figure 4
The spindle thermal inclination sketch.Similarly, the thermal pitch angle in theY-direction can be obtained:
(28)
θ
y
=
(
L
4
i
-
L
2
i
)
-
(
L
4
0
-
L
2
0
)
D
.


## 4.3. Results and Analysis

The spindle speed is a major factor affecting the thermal characteristics. The temperature gradient and the heat generated by the CNC machine spindle are different at different speeds, thus generating different thermal drifts. So the spindle speed affects the temperature field distribution and the magnitude of thermal errors. Meanwhile, the thermal error model is closely related to the intrinsic behavior of thermal drift data. Therefore, a proper test condition shall be proposed so that the thermal error model put forward in this paper can be applied to the practical engineering field and be used to improve machining accuracy. In order to create such kind of testing condition, the distribution of spindle speeds in this paper is set to simulate some common processing conditions of the boring machine. The specific speeds distribution is shown in Figure5.Figure 5
Step speeds distribution.Thermal characteristics of CNC machine tools cover two main aspects of the temperature distribution and thermal error, and the nonuniform distribution of the temperature field gradient is the direct cause of spindle thermal drifts. So for a more comprehensive study of the spindle thermal properties, a simple time-domain analysis of the temperature field is also conducted in the paper. The spindle system temperature variations are shown in Figure6. The overall trends of temperatures on all measuring points increase with time. The rear bearing has the highest temperature reaching 30.4°C due to large capacity, heavy load, and severe friction, which generates more heat, and the following is the motor whose temperature is 27.3°C.Figure 6
Temperatures of the spindle.Figure7 presents the spindle thermal drifts. The time until equilibrium reached is approximately 385 min, with a maximum elongation 39.6 μm. The thermal error on x-axis direction is positive, and its thermal yaw angle is θ
x in the X
O
Z plane; the maximum amount of hot offset error is 35 μm. Thermal error in Y-direction is negative which indicates that, during operation, the spindle is closer to the displacement sensors S
2
/
S
4, and its thermal pitch angle is θ
y in the Y
O
Z plane; the maximum thermal offset is of 20.2 μm.Figure 7
Thermal drifts of the spindle.


## 5. Thermal Errors Prediction and Compensation

After analysis of the experiment, the models of the spindle elongation and thermal angle errors could be established by time series analysis, and three thermal errors in the spindle can be predicted by the proposed models. Subsequently, the final comprehensive thermal error compensation model is induced based on the geometric principle and the time series models.Firstly, the 5 thermal error sequences in Figure7 are converted into 3 series in Figures 8(a)–10(a) based on (27)-(28). The spindle axial elongation time series E
t in the Z-direction was directly from the measured data, and the radial thermal yaw and pitch angle series θ
x
,
t, θ
y
,
t were obtained by applying (27) and (28), respectively. Then the Augmented Dickey-Fuller (ADF) Test Algorithm was applied to identify the stationarity of the thermal error sequences E
t, θ
x
,
t, θ
y
,
t, and the calculation showed that thermal elongation and angles were nonstationary series. Subsequently, the mean μ
i and variance σ
i
2 of the axial elongation E
t and radial thermal declination angle errors θ
x
,
t, θ
y
,
t were calculated by (8). So the original sequences E
t, θ
x
,
t, θ
y
,
t were standardized by (9), and the new time series E
t
′, θ
x
t
′, θ
y
t
′ were shown in Figures 8(b), 9(b), and 10(b).Thermal elongation: (a) original sequence; (b) Gauss standardization.
(a)
(b)Radial thermal yaw angle: (a) original sequence; (b) Gauss standardization.
(a)
(b)Radial thermal pitch angle: (a) original sequence; (b) Gauss standardization.
(a)
(b)When the original sequencesE
t, θ
x
,
t, θ
y
,
t of the spindle thermal drifts were translated into new time series E
t
′, θ
x
t
′, θ
y
t
′ through Gaussian standardization, the Augmented Dickey-Fuller (ADF) Test Algorithm was carried out to identify the E
t
′, θ
x
t
′, θ
y
t
′ stationarity. The calculated results indicated that the characteristic roots of new series were within the unit circle, demonstrating that the standardized series were stationary and they did not need to be smoothed further.The thermal elongation series is calculated by autocorrelation function (ACF) and partial autocorrelation function (PACF); ACF and PACF all have tailing and this indicates that the new standardized sequence of the spindle thermal elongation is ARMA(p
,
q) model, as shown in Figures 11, 12, and 13. Repeat the above process; it exhibits that the radial thermal yaw and pitch angles are also ARMA model.The ACF and PACF map of thermal elongation sequence.
(a)
(b)The ACF and PACF map of thermal yaw angle sequence.
(a)
(b)The ACF and PACF map of thermal pitch angle sequence.
(a)
(b)


### 5.1. Parameters Identification and Model Training

Set order rangep, q
∈
[
0,5
]; then calculate ARMA(p
,
q) model, and the optimal order p and q are determined by Akaika Information Criterion (AIC) [30]. The basic idea is to construct a criterion function AIC(p
,
q); the function considers both original data fitting and the number of unknown parameters in the model. When the parameters and fitting residual variance σ
^
ε
2 make the AIC value reach the minimum, orders p and q are decided:
(29)
AIC
(
p
,
q
)
=
ln
⁡
σ
^
ε
2
+
2
(
p
+
q
+
1
)
n
,
where σ
^
ε
2 is the variance estimation of fitting residual errors, n is the number of samples, and p, q
=
1,2
,
…
,
L are highest orders of the model.After calculation, the new time seriesE
t
′ of the thermal elongation is ARMA
(
5,3
) model, the thermal yaw angle error θ
x
,
t
′ is ARMA
(
2,5
), and the thermal pitch angle error θ
y
,
t
′ is ARMA
(
2,1
); the corresponding models are
(30)
E
t
′
=
1.583
E
t
-
1
′
-
0.471
E
t
-
2
′
+
0.644
E
t
-
3
′
iiiiiiiii
-
1.246
E
t
-
4
′
+
0.488
E
t
-
5
′
+
ε
t
(
1
)
+
0.263
ε
t
-
1
(
1
)
iiiiiiiii
-
0.267
ε
t
-
2
(
1
)
-
0.959
ε
t
-
3
(
1
)
,
θ
x
,
t
′
=
0.18
θ
x
,
t
-
1
′
+
0.779
θ
x
,
t
-
2
′
+
ε
t
(
2
)
iiiiiiiii
+
1.139
ε
t
-
1
(
2
)
+
0.163
ε
t
-
2
(
2
)
-
0.139
ε
t
-
3
(
2
)
iiiiiiiii
+
0.231
ε
t
-
4
(
2
)
+
0.376
ε
t
-
5
(
2
)
,
θ
y
,
t
′
=
1.996
θ
y
,
t
-
1
′
-
0.997
θ
y
,
t
-
2
′
+
ε
t
(
3
)
-
ε
t
-
1
(
3
)
.Assume that the vectorμ includes the means of three time series of the spindle thermal errors and the vector σ includes the standard deviations:
(31)
μ
=
(
μ
1
μ
2
μ
3
)
=
(
-
24.681
-
9.845
-
10.095
)
,
σ
=
(
σ
1
σ
2
σ
3
)
=
(
14.11
6.272
6.571
)
.New time series are reversed according to the following transformation, the final spindle system thermal error model is(32)
(
E
t
θ
x
,
t
θ
y
,
t
)
=
(
σ
1
σ
2
σ
3
)
(
E
t
′
θ
x
,
t
′
θ
y
,
t
′
)
+
(
μ
1
μ
2
μ
3
)
.


### 5.2. Thermal-Induced Error Prediction

The sample data number is 89. And then, the time series models are used to predict the spindle thermal drifts. The fitted curves and the actual measurements are compared in Figures14, 15, and 16.Axial thermal elongation: (a) the prediction and measurement; (b) residual error.
(a)
(b)Radial thermal yaw angle: (a) the prediction and measurement; (b) residual error.
(a)
(b)Radial thermal pitch angle: (a) the prediction and measurement; (b) residual error.
(a)
(b)Now the evaluation criteria of a model fitting are established. Assuming that the absolute value of the residual errors is|
e
i
|, set its minimum as |
e
i
|
min
⁡, maximum as |
e
i
|
max
⁡, and mean value as |
e
i
|
¯. Root mean square error is RMSE, the determination coefficient is R
2, and the predictive ability is η. Consider
(33)
RMSE
=
1
n
∑
i
=
1
n
(
y
i
-
y
~
i
)
2
,
R
2
=
1
-
∑
i
=
1
n
(
y
i
-
y
~
i
)
2
∑
i
=
1
n
(
y
i
-
y
-
i
)
2
,
η
=
1
-
(
1
/
n
)
∑
i
=
1
n
|
y
i
-
y
~
i
|
(
1
/
n
)
∑
i
=
1
n
|
y
i
|
=
1
-
∑
i
=
1
n
|
y
i
-
y
~
i
|
∑
i
=
1
n
|
y
i
|
,
where y
i is the measurement value, y
~
i is the predicted value, y
-
i is the average value of the measurement, i
=
1
,
…
,
n, and n is the number of data points. The fitting performance parameters of the time series are shown in Table 1.Table 1
The fitting performance parameters of time series.
Output
|
e
i
|
min
⁡  (μ
m
/
′
′)
|
e
i
|
max
⁡  (μ
m
/
′
′)
|
e
i
|
¯  (μ
m
/
′
′)
RMSE
R
2
η  (%)
θ
x
0.00
3.44
0.441
0.730
0.986
95.6
θ
y
0.00
1.62
0.242
0.416
0.996
97.6
E
0.00
1.50
0.335
0.473
0.999
98.7The absolute mean values of the residual errors are small, and the RMSE is similarly closed to zero; the coefficient of determinationR
2 is close to 1. In addition, the model predictive ability is more than 90% in the three different directions, which indicates that the time series model has a higher prediction accuracy.


### 5.3. Thermal Error Compensation Implementation

Figure17 is a schematic diagram of the spindle thermal error compensation, and the CNC system is the Siemens 840D. The temperature module acquires signal fromPT100 and sends it to the CNC system byRS-232. A thermal error compensation module is embedded into CNC based on secondary development of 840D, and it can receive error compensation parameters and passes them to PLC. Finally, the thermal error offsets are calculated and sent to the CNC to achieve compensation by PLC. While the thermal yaw and pitch errors are translated into the components of coordinate axis, three components are compensated by the principle described by this compensation system.Figure 17
Thermal error compensation control.Assuming that the distances between the deflexion center and the spindle nose areD
0
x and D
0
y in X
/
Y-direction, respectively, as is shown in Figure 1, there is
(34)
D
0
x
=
Δ
L
1
tan
θ
x
-
D
L
1
=
548.659
mm
,
D
0
y
=
Δ
L
2
tan
θ
y
-
D
L
2
=
508.706
mm
.
The distances between displacement sensors S
1, S
2 and the spindle nose are D
L
1, D
L
2, respectively, and D
L
1
=
D
L
2
=
243.158 mm.The thermal components of the coordinate can be obtained by applying (32) and (5)–(7) as follows:
(35)
Δ
O
x
=
(
D
+
548.66
)
(
σ
2
θ
x
,
t
′
+
μ
2
)
,
Δ
O
y
=
(
D
+
508.71
)
(
σ
3
θ
y
,
t
′
+
μ
3
)
,
Δ
O
z
=
σ
1
E
t
′
+
μ
1
.
Set the coordinate of the original point W on the workpiece as (
P
x
,
P
y
,
P
z
); then the new coordinate for thermal error compensation is W
′
(
P
x
′
,
P
y
′
,
P
z
′
):
(36)
P
x
′
=
P
x
+
Δ
O
x
=
P
x
+
(
D
+
548.66
)
(
σ
2
θ
x
,
t
′
+
μ
2
)
,
P
y
′
=
P
y
+
Δ
O
y
=
P
y
+
(
D
+
508.71
)
(
σ
3
θ
y
,
t
′
+
μ
3
)
P
z
′
=
P
z
+
Δ
O
z
=
P
z
+
σ
1
E
t
′
+
μ
1
.After the thermal compensation, errors were reduced significantly, which are shown in Figures18–20. The spindle thermal drift measuring results before and after error compensation in the three directions are shown in Figures 18 and 20; the red curves are the original measurements before the compensation of the spindle system, while the blue curves are the measurements after the spindle compensation. The existing literatures mainly compensated the axial elongation, which was similar to the results in Figure 18, and did not propose a model or compensate for the radial thermal errors. If there was no thermal error compensation in radial directions, namely, unconsidering the radial thermal tilt angle errors, after the completion of the axial error compensation, the radial thermal errors are still the red curves in Figures 19-20. In this paper, we not only compensate the axial elongation, but also establish thermal error compensation equations in three directions of the spindle taking into account the radial thermal tilt angle errors and the length of cutting tools; the equations are (36). The newly measured thermal drifts are the blue curves shown in Figures 19-20 after the error compensation in radial X
/
Y-directions. Comparing the two color curves, it is obvious that the errors of the blue curves which have considered the thermal tilt angles are much smaller than the errors of the red curves which have not considered the thermal tilt angles. The maximum error decreased from 39.6 μm to 7 μm in axial direction, and the average error reduced from 24.6 μm to 2.8 μm; namely, the average offset is about 22 μm. Axial accuracy is improved by 89.7%, which demonstrates the method of the proposed measurement and modeling effectively. The absolute average thermal error value S
3 in the radial X-direction reduced from 14.6 μm to 3.3 μm, and the accuracy is improved by 77.4%. Meanwhile, the thermal error absolute maximum value S
4 in the radial Y-direction declined 12.1 μm into 3.7 μm, and the accuracy is advanced by 86%.Figure 18
Axial thermal elongation compensation.Figure 19
Radial thermal error compensation inX-direction.Figure 20
Radial thermal error compensation inY-direction.


## 5.1. Parameters Identification and Model Training

Set order rangep, q
∈
[
0,5
]; then calculate ARMA(p
,
q) model, and the optimal order p and q are determined by Akaika Information Criterion (AIC) [30]. The basic idea is to construct a criterion function AIC(p
,
q); the function considers both original data fitting and the number of unknown parameters in the model. When the parameters and fitting residual variance σ
^
ε
2 make the AIC value reach the minimum, orders p and q are decided:
(29)
AIC
(
p
,
q
)
=
ln
⁡
σ
^
ε
2
+
2
(
p
+
q
+
1
)
n
,
where σ
^
ε
2 is the variance estimation of fitting residual errors, n is the number of samples, and p, q
=
1,2
,
…
,
L are highest orders of the model.After calculation, the new time seriesE
t
′ of the thermal elongation is ARMA
(
5,3
) model, the thermal yaw angle error θ
x
,
t
′ is ARMA
(
2,5
), and the thermal pitch angle error θ
y
,
t
′ is ARMA
(
2,1
); the corresponding models are
(30)
E
t
′
=
1.583
E
t
-
1
′
-
0.471
E
t
-
2
′
+
0.644
E
t
-
3
′
iiiiiiiii
-
1.246
E
t
-
4
′
+
0.488
E
t
-
5
′
+
ε
t
(
1
)
+
0.263
ε
t
-
1
(
1
)
iiiiiiiii
-
0.267
ε
t
-
2
(
1
)
-
0.959
ε
t
-
3
(
1
)
,
θ
x
,
t
′
=
0.18
θ
x
,
t
-
1
′
+
0.779
θ
x
,
t
-
2
′
+
ε
t
(
2
)
iiiiiiiii
+
1.139
ε
t
-
1
(
2
)
+
0.163
ε
t
-
2
(
2
)
-
0.139
ε
t
-
3
(
2
)
iiiiiiiii
+
0.231
ε
t
-
4
(
2
)
+
0.376
ε
t
-
5
(
2
)
,
θ
y
,
t
′
=
1.996
θ
y
,
t
-
1
′
-
0.997
θ
y
,
t
-
2
′
+
ε
t
(
3
)
-
ε
t
-
1
(
3
)
.Assume that the vectorμ includes the means of three time series of the spindle thermal errors and the vector σ includes the standard deviations:
(31)
μ
=
(
μ
1
μ
2
μ
3
)
=
(
-
24.681
-
9.845
-
10.095
)
,
σ
=
(
σ
1
σ
2
σ
3
)
=
(
14.11
6.272
6.571
)
.New time series are reversed according to the following transformation, the final spindle system thermal error model is(32)
(
E
t
θ
x
,
t
θ
y
,
t
)
=
(
σ
1
σ
2
σ
3
)
(
E
t
′
θ
x
,
t
′
θ
y
,
t
′
)
+
(
μ
1
μ
2
μ
3
)
.


## 5.2. Thermal-Induced Error Prediction

The sample data number is 89. And then, the time series models are used to predict the spindle thermal drifts. The fitted curves and the actual measurements are compared in Figures14, 15, and 16.Axial thermal elongation: (a) the prediction and measurement; (b) residual error.
(a)
(b)Radial thermal yaw angle: (a) the prediction and measurement; (b) residual error.
(a)
(b)Radial thermal pitch angle: (a) the prediction and measurement; (b) residual error.
(a)
(b)Now the evaluation criteria of a model fitting are established. Assuming that the absolute value of the residual errors is|
e
i
|, set its minimum as |
e
i
|
min
⁡, maximum as |
e
i
|
max
⁡, and mean value as |
e
i
|
¯. Root mean square error is RMSE, the determination coefficient is R
2, and the predictive ability is η. Consider
(33)
RMSE
=
1
n
∑
i
=
1
n
(
y
i
-
y
~
i
)
2
,
R
2
=
1
-
∑
i
=
1
n
(
y
i
-
y
~
i
)
2
∑
i
=
1
n
(
y
i
-
y
-
i
)
2
,
η
=
1
-
(
1
/
n
)
∑
i
=
1
n
|
y
i
-
y
~
i
|
(
1
/
n
)
∑
i
=
1
n
|
y
i
|
=
1
-
∑
i
=
1
n
|
y
i
-
y
~
i
|
∑
i
=
1
n
|
y
i
|
,
where y
i is the measurement value, y
~
i is the predicted value, y
-
i is the average value of the measurement, i
=
1
,
…
,
n, and n is the number of data points. The fitting performance parameters of the time series are shown in Table 1.Table 1
The fitting performance parameters of time series.
Output
|
e
i
|
min
⁡  (μ
m
/
′
′)
|
e
i
|
max
⁡  (μ
m
/
′
′)
|
e
i
|
¯  (μ
m
/
′
′)
RMSE
R
2
η  (%)
θ
x
0.00
3.44
0.441
0.730
0.986
95.6
θ
y
0.00
1.62
0.242
0.416
0.996
97.6
E
0.00
1.50
0.335
0.473
0.999
98.7The absolute mean values of the residual errors are small, and the RMSE is similarly closed to zero; the coefficient of determinationR
2 is close to 1. In addition, the model predictive ability is more than 90% in the three different directions, which indicates that the time series model has a higher prediction accuracy.


## 5.3. Thermal Error Compensation Implementation

Figure17 is a schematic diagram of the spindle thermal error compensation, and the CNC system is the Siemens 840D. The temperature module acquires signal fromPT100 and sends it to the CNC system byRS-232. A thermal error compensation module is embedded into CNC based on secondary development of 840D, and it can receive error compensation parameters and passes them to PLC. Finally, the thermal error offsets are calculated and sent to the CNC to achieve compensation by PLC. While the thermal yaw and pitch errors are translated into the components of coordinate axis, three components are compensated by the principle described by this compensation system.Figure 17
Thermal error compensation control.Assuming that the distances between the deflexion center and the spindle nose areD
0
x and D
0
y in X
/
Y-direction, respectively, as is shown in Figure 1, there is
(34)
D
0
x
=
Δ
L
1
tan
θ
x
-
D
L
1
=
548.659
mm
,
D
0
y
=
Δ
L
2
tan
θ
y
-
D
L
2
=
508.706
mm
.
The distances between displacement sensors S
1, S
2 and the spindle nose are D
L
1, D
L
2, respectively, and D
L
1
=
D
L
2
=
243.158 mm.The thermal components of the coordinate can be obtained by applying (32) and (5)–(7) as follows:
(35)
Δ
O
x
=
(
D
+
548.66
)
(
σ
2
θ
x
,
t
′
+
μ
2
)
,
Δ
O
y
=
(
D
+
508.71
)
(
σ
3
θ
y
,
t
′
+
μ
3
)
,
Δ
O
z
=
σ
1
E
t
′
+
μ
1
.
Set the coordinate of the original point W on the workpiece as (
P
x
,
P
y
,
P
z
); then the new coordinate for thermal error compensation is W
′
(
P
x
′
,
P
y
′
,
P
z
′
):
(36)
P
x
′
=
P
x
+
Δ
O
x
=
P
x
+
(
D
+
548.66
)
(
σ
2
θ
x
,
t
′
+
μ
2
)
,
P
y
′
=
P
y
+
Δ
O
y
=
P
y
+
(
D
+
508.71
)
(
σ
3
θ
y
,
t
′
+
μ
3
)
P
z
′
=
P
z
+
Δ
O
z
=
P
z
+
σ
1
E
t
′
+
μ
1
.After the thermal compensation, errors were reduced significantly, which are shown in Figures18–20. The spindle thermal drift measuring results before and after error compensation in the three directions are shown in Figures 18 and 20; the red curves are the original measurements before the compensation of the spindle system, while the blue curves are the measurements after the spindle compensation. The existing literatures mainly compensated the axial elongation, which was similar to the results in Figure 18, and did not propose a model or compensate for the radial thermal errors. If there was no thermal error compensation in radial directions, namely, unconsidering the radial thermal tilt angle errors, after the completion of the axial error compensation, the radial thermal errors are still the red curves in Figures 19-20. In this paper, we not only compensate the axial elongation, but also establish thermal error compensation equations in three directions of the spindle taking into account the radial thermal tilt angle errors and the length of cutting tools; the equations are (36). The newly measured thermal drifts are the blue curves shown in Figures 19-20 after the error compensation in radial X
/
Y-directions. Comparing the two color curves, it is obvious that the errors of the blue curves which have considered the thermal tilt angles are much smaller than the errors of the red curves which have not considered the thermal tilt angles. The maximum error decreased from 39.6 μm to 7 μm in axial direction, and the average error reduced from 24.6 μm to 2.8 μm; namely, the average offset is about 22 μm. Axial accuracy is improved by 89.7%, which demonstrates the method of the proposed measurement and modeling effectively. The absolute average thermal error value S
3 in the radial X-direction reduced from 14.6 μm to 3.3 μm, and the accuracy is improved by 77.4%. Meanwhile, the thermal error absolute maximum value S
4 in the radial Y-direction declined 12.1 μm into 3.7 μm, and the accuracy is advanced by 86%.Figure 18
Axial thermal elongation compensation.Figure 19
Radial thermal error compensation inX-direction.Figure 20
Radial thermal error compensation inY-direction.


## 6. Conclusions

The spindle thermal error modeling containing axial elongation and radial thermal angle errors is more suitable for actual conditions, because it could exactly describe the space-pose of the thermal deformation and consequently be utilized to compensate spindle thermal drifts, improving the machining accuracy. But radial thermal-induced angle errors were ignored in current literatures. To solve this problem, the five-point method was applied to measure the spindle thermal drifts, and the thermal error model including axial elongation and radial yaw and pitch angle errors was proposed based on the time series analysis. The time series thermal error model can fully exploit the inherent dynamic characteristics of the spindle system thermal deformation, so the prediction of the presented modeling could reach up to 90% with excellent generalization and robustness. In addition, considering the length of the cutting tools and the radial thermal angle errors, the final mathematical compensation equation of the spindle thermal drifts was presented, and the real-time compensation was implemented. The result indicated that the axial machining accuracy was improved by 89.7%; theX
/
Y-direction accuracy can reach up to 77.4% and 86%, respectively, which demonstrated that the proposed methodology of measurement, modeling, and compensation was effective.



---

*Source: 784218-2014-08-27.xml*