Measuring Point Planning and Fitting Optimization of the Flange and Spigot Structures of Aeroengine Rotors

Abstract

An optimized measuring point planning and fitting method for rotor flange and spigot structures was proposed to achieve precise measurement of position and pose of the aeroengine rotors during docking processes. Firstly, the impact of circumferential phase angle, distribution range angle, total number of measuring points, and number of distribution rings on measurement uncertainty was analyzed. The measuring point planning schemes for flange and spigot were proposed. Secondly, the Gauss Newton iterative solution principle considering damping factors was clarified. Subsequently, an optimized iterative reweighting method consisting of weight iterative estimation, singular value detection under the Chauvenet criterion, and clustering detection was proposed for fitting the flange annular end face. A mapping point total least squares method with practical geometric significance was proposed for fitting the spigot cylinder face. Finally, measuring and fitting experiments were performed. The singular measuring point detection methods were verified. Under the optimized fitting methods, the goodness of fit and average orthogonal distance of flange and spigot structures are 0.756 and 0.089 mm, respectively, which have higher fitting accuracy than the other traditional methods.

1. Introduction

Flange and spigot are the typical positioning structures at the end of an aeroengine rotor, realizing axis orientation, axial positioning, and radial positioning between adjacent rotors [1,2,3]. A large amount of measurement work is required during the rotor docking assembly processes, to obtain the initial position and pose of each rotor and to provide data for docking trajectories [4,5,6]. Therefore, the measuring and fitting accuracy is crucial, directly determines the assembly accuracy of the rotors, and has a significant impact on the dynamic performance of the entire engine [7,8,9,10].

Automated systems have gradually been applied in the measurement of aeroengine rotors [11,12,13,14,15]. Among these, the multi-axis precision mechanical contact measurement system and laser measurement system have the advantages of high reliability and automation [16,17,18]. A large number of point cloud coordinates are generated after measurement. The question of how to arrange the number and position of measuring points is a key issue in balancing measurement accuracy and efficiency [19,20,21,22]. Based on the measuring point planning, optimizing the corresponding fitting methods for different structures is an effective way to further reduce the uncertainty of model parameters [23,24,25]. In summary, reasonable measuring point planning and optimized fitting methods of the flange and spigot structures are important prerequisites for obtaining the rotor position and pose information quickly and accurately.

In recent years, some scholars have conducted relevant research on the measuring point planning strategy of aeroengine rotors [26,27]. Zhang et al. have proposed a redundant measuring point detection method for the measurement of rotor end faces, and the optimized algorithm improved the accuracy by 21 μm compared with the traditional least squares (LS) fitting method [28]. Zhu et al. constructed a heuristic floating forward search algorithm to efficiently find a near-optimal solution for measuring point layout [29]. Ding et al. used the eigenvalue method to process the point cloud data of the rotor end face, which has the function of excluding gross error points in the measurement points [30]. However, the existing measuring point planning methods mainly deal with integral surfaces with simple geometric shapes, and their applicability to complex spatial feature structures such as flange and spigot is limited.

With regard to the research regarding the fitting methods for typical axes and circular contours: Xie et al. established a spatial axis LS fitting model with non-cylindrical data filtered out [31]. The Limacon model is widely used in circular profile measurement to characterize eccentricity parameters [32,33]. Chetwynd et al. improved the parameter estimation method of LS fitting radius in the Limacon model [34]. Sun et al. proposed the vector projection method in cylindrical contour fitting. Compared with the traditional Limacon model, the accuracy of rotor roundness has been improved by 2.2 μm, and the accuracy of coaxiality has been improved by 9.39 μm [35,36]. Zhang et al. constructed a particle swarm optimization algorithm for fitting aeroengine rotor models based on a two-dimension laser displacement sensor, and the measurement accuracy of cylindricity error was improved by 1.768 μm compared with the traditional method [37]. Wang et al. used the equivalent homogenization processing and 3D LS method for fitting aeroengine rotor annular plane [38]. For different geometric features, researchers also carried out comparative studies on various fitting methods based on the same set of point cloud data. Song et al. compared the fitting performances of six normal line fitting methods on different types of geometries, e.g., sphere, cylinder, and ellipsoid [39]. Nouira et al. compared the fitting accuracy of different models such as LS circle, minimum area circle, maximum inscribed circle and minimum circumscribed circle, and proposed a fitting method based on a micro displacement model [40]. The above research constructed the corresponding fitting algorithms for the points with specific distribution forms, achieving higher fitting accuracy compared with traditional methods. However, few targeted combinatorial optimization algorithms have been established to meet the requirements for fitting the spatial position and pose of the rotor flange and spigot structures. Moreover, it is no longer possible to measure the flange and spigot again after the rotors are docked, so higher fitting accuracy is required in a single fitting process.

In this paper, a systematic study on measuring point planning and fitting methods for typical rotor flange and spigot structures was carried out. First, the impact of different measuring point distribution parameters on the model uncertainty was analyzed and a measuring point planning strategy was proposed. Then, the optimized iterative reweighting method for fitting flange and the mapping point total least squares method for fitting spigot were proposed, respectively. Finally, the measurement and fitting experiments were performed. The singular value point detection methods were compared and the fitting accuracy of the optimized fitting methods was verified.

This paper is organized as follows. In Section 2, the measuring point planning method is proposed; in Section 3, the optimized fitting algorithm is proposed; in Section 4, the measurement experiments are performed; and in Section 5, the conclusions are presented.

2. Measuring Point Planning Method

The schematic that illustrates the measurement of the flange and spigot structures of the aeroengine rotors is shown in Figure 1. With regard to the measuring point planning of the flange annular end face. The end face description equation can be represented by z_M = a₁ + a₂x_M + a₃y_M in the measurement coordinate system O_M-xyz. Without loss of generality, allow N measuring points M_i (x_Mi, y_Mi, z_Mi) distributed on the plane z_M = 0 to be fitted. To simplify the calculation, the case of N = 3 is analyzed first. The measuring points M₁ (0, 0, 0), M₂ (x_M2, 0, 0) and M₃ $\left(x_{\mathrm{M}3},\pm \left[\sqrt{R_{\mathrm{C}}{}^2-{(x_{\mathrm{M}3}-x_{\mathrm{M}2}/2)}^2}+\sqrt{R_{\mathrm{C}}{}^2-x_{\mathrm{M}2}{}^2/4}\right],0\right)$ are shown in Figure 2a. The equation of the annular end face is transformed into a₁ = −(Ax_M + Bx_M + Cx_M)/C, a₂ = −A/C, a₃ = −B/C. Where A, B, and C are the model parameters of the plane general equation. The model parameter uncertainty is defined as the standard deviation of the fitting parameter and is associated with the measurement system error σ_M. The uncertainty of the parameter of the flange plane general equation is shown in Equations (1)–(3).

(1)$\begin{array}{ll}{\delta }_A& =\sqrt{{\displaystyle \sum _{i=1}^3\left[{\left(\frac{\partial A}{\partial x_i}\right)}^2+{\left(\frac{\partial A}{\partial y_i}\right)}^2+{\left(\frac{\partial A}{\partial z_i}\right)}^2\right]}}{\sigma }_{\mathrm{M}}\\ & =\left(\sqrt{2R_{\mathrm{C}}{}^2-2{(x_{\mathrm{M}3}-x_{\mathrm{M}2}/2)}^2}+\sqrt{2{R_{\mathrm{C}}}^2-x_{\mathrm{M}2}{}^2/2}\right){\sigma }_{\mathrm{M}}\end{array}$

(2)$\begin{array}{ll}{\delta }_B& =\sqrt{{\displaystyle \sum _{i=1}^3\left[{\left(\frac{\partial B}{\partial x_i}\right)}^2+{\left(\frac{\partial B}{\partial y_i}\right)}^2+{\left(\frac{\partial B}{\partial z_i}\right)}^2\right]}}{\sigma }_{\mathrm{M}}\\ & =\sqrt{2\left(x_{\mathrm{M}2}{}^2+x_{\mathrm{M}3}{}^2-x_{\mathrm{M}2}{x}_{\mathrm{M}3}\right)}{\sigma }_{\mathrm{M}}\end{array}$

(3)$\begin{array}{ll}{\delta }_C& =\sqrt{{\displaystyle \sum _{i=1}^3\left[{\left(\frac{\partial C}{\partial x_i}\right)}^2+{\left(\frac{\partial C}{\partial y_i}\right)}^2+{\left(\frac{\partial C}{\partial z_i}\right)}^2\right]}}{\sigma }_{\mathrm{M}}\\ & =\sqrt{\left(8R_{\mathrm{C}}{}^2+2x_{\mathrm{M}2}{}^2\right)+4\sqrt{\left[2R_{\mathrm{C}}{}^2-2{(x_{\mathrm{M}3}-x_{\mathrm{M}2}/2)}^2\right]\left(2R_{\mathrm{C}}{}^2-x_{\mathrm{M}2}{}^2/2\right)}}{\sigma }_{\mathrm{M}}\end{array}$

It can be calculated that when $x_{\mathrm{M}3}=x_{\mathrm{M}2}/2$, $x_{\mathrm{M}2}{}^2+x_{\mathrm{M}3}{}^2-x_{\mathrm{M}2}{x}_{\mathrm{M}3}$ has the minimum value, both δ_B and δ_C have the minimum values simultaneously. The uncertainty of attitude column matrix a representing the fitted normal vector of the rotor end face is shown in Equations (4)–(6).

(4)${\delta }_{ax}=\sqrt{2}{\sigma }_{\mathrm{M}}/x_{\mathrm{M}2}$

(5)${\delta }_{ay}=\frac{\sqrt{6}}{\sqrt{4R_{\mathrm{C}}{}^2-{(2x_{\mathrm{M}3}-x_{\mathrm{M}2})}^2}+\sqrt{4R_{\mathrm{C}}{}^2-x_{\mathrm{M}2}{}^2}}{\sigma }_{\mathrm{M}}$

(6)${\delta }_{az}=0$

The increase in x_M2 and y_M3 leads to a decrease in uncertainty δ_ax and δ_ay, i.e., the relative distance of the measuring points M₁, M₂, and M₃ should be as far as possible. Then δ_ax = δ_ay and x_M2 = ±$\sqrt{3}$R_C and the angle among O_M-M_i is 2π/3. When the calculation method is generalized to N ≥ 3, it can be concluded that the uncertainty of the end face parameters is minimized when the measuring points are uniformly distributed along the circumference.

For situations where the measurement trajectory is limited, measuring a large arc range can suffice instead of measuring the entire circle. The analysis of the impacts of the number of measuring points and the uniform distribution range on the uncertainty of the end face model parameters are necessary. The measuring points are arranged symmetrically as shown in Figure 2b. Parameter a₁ determines the third-row element a_z in the workpiece pose column matrix a, which is an important parameter for determining the spatial angle of the end face axis (the angle between the end face normal vector and the coordinate axis O_M-z). The uncertainty δ_a1 of parameter a₁ using the LS method is calculated as Equation (7).

(7)${\delta }_{a1}=\sqrt{{\displaystyle \sum _{i=1}^N\left[\frac{1}{N}-\frac{{\left({\displaystyle \sum _{i=1}^N{R}_{\mathrm{C}}\cos {\theta }_i}\right)}^2/N-R_{\mathrm{C}}\cos {\theta }_i{\displaystyle \sum _{i=1}^N{R}_{\mathrm{C}}\cos {\theta }_i}}{N{\displaystyle \sum _{i=1}^N{\left(R_{\mathrm{C}}\cos {\theta }_i\right)}^2-{\left({\displaystyle \sum _{i=1}^N{R}_{\mathrm{C}}\cos {\theta }_i}\right)}^2}}\right]}}{\sigma }_{\mathrm{M}}$

where N is the number of measuring points, R_C is the radius of the circular path, θ_i is the circumferential angle of the measuring point M_i, and σ_M represents measurement systematic error.

Taking the flange end face of an aeroengine low pressure turbine (LPT) mid shaft as an example, the measuring path radius R_C = 149 mm. The relationship between the calculated uncertainty δ_a1 and the number of measuring points N/the angle of uniform distribution range φ is shown in Figure 3a. δ_a1 reduces significantly when φ ≥ 180°. Therefore, the measurement range should cover at least half the arc area. The reduction rate of the minimum value of the model parameter uncertainty min(δ_a1/σ_M) gradually decreases when N increases step by step in the arithmetic progression of 4, 12, 20, 28, 36. The reduction rate of min(δ_a1/σ_M) is less than 10% when N increases from 28 to 36. Taking into account the impact on measurement accuracy and efficiency, N = 28 is a reasonable value. Due to the radial width on the circular end face, the influence of the number of measuring rings distributed along multiple concentric circles/arcs is analyzed while N = 28. The uncertainty is shown in Figure 3b when R_C = 149, 137 and 125 mm, respectively. When 30° ≤ φ ≤ 80°, δ_a1 under a multi-ring distribution reduces significantly compared with the single ring scheme and δ_a1 under a double-ring distribution is lower than under a three-ring distribution. Due to the limitation of the flange edge width (24 mm) in this example, when the total number of measuring points remains unchanged, the excessive rings result in a small radial spacing and a large circumferential spacing of measuring points. Therefore, the arrangement of measuring points along the double concentric rings is the optimal method and has min(δ_a1/σ_M) = 0.188. In summary, the N = 28 measuring points evenly distributed in double rings form a reasonable measuring point planning scheme, as shown in Figure 2c. Similarly, the optimized measuring point planning strategies for end faces at different positions (shown in Figure 4) in an aeroengine LPT rotor assembly are analyzed, as listed in Table 1.

With regard to the measuring point planning of the spigot cylindrical face. Without loss of generality, let N points M_{proj_i} (x_{proj_Mi}, y_{proj_Mi}, z_{proj_Mi}) projected by the measuring points be distributed on the assumed plane z_M = 0. To simplify the calculation, the case of N = 3 is analyzed first. The mapping points M_{proj_1} (0, R_{proj_M}, 0), M_{proj_2} (R_{proj_M}cosθ₂, R_{proj_M}sinθ₂, 0), and M_{proj_3} (R_{proj_M}cosθ₃, R_{proj_M}sinθ₃, 0) are shown in Figure 2d. The coordinate uncertainty of x_{proj_MO} and y_{proj_MO} in the model parameters of the standard circle is shown in Equations (8) and (9).

(8)${\delta }_x=\sqrt{{\displaystyle \sum _{i=1}^3\left[{\left(\frac{\partial x_{\mathrm{proj\_M}O}}{\partial x_{\mathrm{proj\_M}i}}\right)}^2+{\left(\frac{\partial x_{\mathrm{proj\_M}O}}{\partial y_{\mathrm{proj\_M}i}}\right)}^2\right]}}{\sigma }_{\mathrm{M}}$

(9)${\delta }_y=\sqrt{{\displaystyle \sum _{i=1}^3\left[{\left(\frac{\partial y_{\mathrm{proj\_M}O}}{\partial x_{\mathrm{proj\_M}i}}\right)}^2+{\left(\frac{\partial y_{\mathrm{proj\_M}O}}{\partial y_{\mathrm{proj\_M}i}}\right)}^2\right]}}{\sigma }_{\mathrm{M}}$

Thus, the uncertainty δ_O of the center coordinates of the mapped circle is calculated as:

(10)${\delta }_O=\sqrt{{\delta }_x^2+{\delta }_y^2}=\sqrt{\frac{6-2\cos \left({\theta }_3-{\theta }_2\right)-2\cos {\theta }_2-2\cos {\theta }_3}{{\left[\sin {\theta }_2-\sin {\theta }_3+\sin \left({\theta }_3-{\theta }_2\right)\right]}^2}}{\sigma }_{\mathrm{M}}$

When ∂δ_O/∂θ₂ = 0, the uncertainty δ_O has a minimum value if θ₂ = θ₃/2. The relationship between δ_O and circumferential phase angle θ₃ is shown in Figure 3c. When θ₂ = 2π/3, θ₃ = 4π/3, δ_O has the minimum value 1.155σ_M. When the calculation method is generalized to N ≥ 3, it can be concluded that the uncertainty of the center coordinate has a minimum value when the measuring points are uniformly distributed along the circumference.

The number of points and the uniformly distributed range angle also have an impact on the uncertainty δ_O. The mapping points are arranged symmetrically as shown in Figure 2e. δ_O is calculated as Equation (11) when using the LS method.

(11)${\delta }_O=2\sqrt{\frac{N-{\left({\displaystyle \sum _{i=1}^N\cos {\theta }_i}\right)}^2/N}{\left(N-{\displaystyle \sum _{i=1}^N\cos 2{\theta }_i}\right)\left[N+{\displaystyle \sum _{i=1}^N\cos 2{\theta }_i}-2{\left({\displaystyle \sum _{i=1}^N\cos {\theta }_i}\right)}^2/N\right]}}{\sigma }_{\mathrm{M}}$

The relationship between the uncertainty δ_O and the number of measuring points N/the angle of uniform distribution range φ is shown in Figure 3d. δ_O reduces significantly when φ ≥ 180°. Therefore, the measurement range should cover at least half the arc area. The reduction rate of the minimum value of uncertainty min(δ_O/σ_M) gradually decreases when N increases step by step in the arithmetic progression of 4, 12, 20, 28, 36, 44. The reduction rate of min(δ_O/σ_M) is less than 10% when N increases from 36 to 44. Taking into account the impact on measurement accuracy and efficiency, N = 36 is a reasonable value. Similarly, the optimized measuring point planning strategies for cylindrical faces at different positions (shown in Figure 4) in an aeroengine LPT rotor assembly are analyzed, as listed in Table 2.

3. Optimization of Fitting Algorithm

3.1. LS Principle and Gauss Newton Iterative Solution Method Considering Damping Factors

Measuring point fitting can be considered a minimization process. A set of observation values z_Mi are obtained under the given condition $\left\{\left(x_{\mathrm{M}i},y_{\mathrm{M}i}\right)|i=1,2,\dots ,N\right\}$. The target is the minimization of the deviation between the values obtained from the mathematical model and the actual observation values according to a calculation rule. The LS minimization mathematical model is shown in Equation (12). The optimization fitting algorithm in Section 3.2 and Section 3.3 are all constructed based on the LS fitting principle.

(12)${\chi }^2={\displaystyle \sum _{i=1}^N{w}_i{\left[f\left(x_{\mathrm{M}i},y_{\mathrm{M}i}|\mathit{a}\right)-z_{\mathrm{M}i}\right]}^2}\to \mathrm{Min}$

The solution process of the Gauss Newton method for searching the minimum value χ² is shown in Figure 5. The gradient of Δa = 0 when F(a + Δa) is minimized. The model parameter adjustment matrix Δa = −H⁻¹g. The gradient vector g and Hessian matrix H are decomposed into the Jacobian matrix J_LS, as shown in Equation (13).

(13)${\mathit{J}}_{\mathrm{LS}}=\left[\begin{array}{cccc}\frac{\partial f\left(x_{\mathrm{M}1},y_{\mathrm{M}1}|\mathit{a}\right)}{\partial a_1}& \frac{\partial f\left(x_{\mathrm{M}1},y_{\mathrm{M}1}|\mathit{a}\right)}{\partial a_2}& \cdots & \frac{\partial f\left(x_{\mathrm{M}1},y_{\mathrm{M}1}|\mathit{a}\right)}{\partial a_M}\\ \frac{\partial f\left(x_{\mathrm{M}2},y_{\mathrm{M}2}|\mathit{a}\right)}{\partial a_1}& \frac{\partial f\left(x_{\mathrm{M}2},y_{\mathrm{M}2}|\mathit{a}\right)}{\partial a_2}& \cdots & \frac{\partial f\left(x_{\mathrm{M}2},y_{\mathrm{M}2}|\mathit{a}\right)}{\partial a_M}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f\left(x_{\mathrm{M}N},y_{\mathrm{M}N}|\mathit{a}\right)}{\partial a_1}& \frac{\partial f\left(x_{\mathrm{M}N},y_{\mathrm{M}N}|\mathit{a}\right)}{\partial a_2}& \cdots & \frac{\partial f\left(x_{\mathrm{M}N},y_{\mathrm{M}N}|\mathit{a}\right)}{\partial a_M}\end{array}\right]$

Model parameter iteration adjustment matrix Δa is shown in Equation (14). The iteration termination target is shown in Equation (16).

(14)$\Delta \mathit{a}={\left({\mathit{J}}_{\mathrm{LS}}^{\mathrm{T}}\mathit{W}{\mathit{J}}_{\mathrm{LS}}\right)}^{-1}{\mathit{J}}_{\mathrm{LS}}^{\mathrm{T}}\mathit{W}\mathit{r}$

(15)$\mathit{a}+\Delta \mathit{a}\to \mathit{a}$

(16)$\left|\Delta a_j\right|<\epsilon ,\forall j$

where W is the weight diagonal matrix; r is the residual column vector; and ε is the iteration termination value.

The traditional Gauss Newton algorithm may fall into saddle points or iterate in wrong directions. The damping factor μ is introduced in the Levenberg–Marquardt method. Thus, the iterative adjustment matrix Δa is revised to:

(17)$\Delta {\mathit{a}}_{\mathrm{LM}}={\left({\mathit{J}}_{\mathrm{LS}}^{\mathrm{T}}\mathit{W}{\mathit{J}}_{\mathrm{LS}}+\mu \mathit{E}\right)}^{-1}{\mathit{J}}_{\mathrm{LS}}^{\mathrm{T}}\mathit{W}\mathit{r}$

The damping factor μ is determined by the adaptive coefficient λ_μ and the maximum diagonal element n_max of the quadratic normal matrix $\mathit{N}={\mathit{J}}_{\mathrm{LS}}^{\mathrm{T}}\mathit{W}{\mathit{J}}_{\mathrm{LS}}^{}$, as shown in Equation (18). The empirical value of λ_μ can be set to 0.001 in the first iteration. The direction is correct if χ²(a + Δa) < χ²(a) after iteration. Then update a + Δa→a and decrease λ_μ in the next iteration. Otherwise, update a and increase λ_μ.

(18)$\mu ={\lambda }_{\mu }{n}_{\max }$

3.2. Optimized Fitting Method for Flange Annular End Face

The weight w_i is theoretically independent of each observation value z_Mi and is inversely proportional to the square of the standard uncertainty of the probability distribution function of z_Mi. However, the uncertainty cannot be directly calculated in practical applications. Therefore, iterative operations can be used to obtain a more accurate estimated weight. The iterative mathematical model makes the singular points deviate more significantly. Thus, singular points can be removed accurately before weighted fitting is performed again. The optimized iterative reweighted LS process for the flange annular end face is shown in Figure 6.

3.2.1. Iterative Estimation of Weights Based on Residuals

The residual Δ_i is approximately mapped to the weight w_i. When the observed value z_Mi under (x_Mi, y_Mi) is very close to the model function, it causes 1/Δ_i to become too high. It is therefore necessary to set a reasonable upper limit value for w_i:

(19)$w_i=\{\begin{array}{c}\frac{1}{{\lambda }_{\mathrm{L}}^2},\mathrm{if}\left|{\Delta }_i\right|<{\lambda }_{\mathrm{L}}\\ \frac{1}{{\Delta }_i^2},\mathrm{if}\left|{\Delta }_i\right|\ge {\lambda }_{\mathrm{L}}\end{array}$

The residual Δ_i obeys a normal distribution due to the influence of independent random factors. Its lower bound value λ_L is determined by Equation (20). To simplify operations, λ_L can be set as the median of absolute residual $\left|{\Delta }_i\right|$, causing 50% of the observed values to have the same weight, $1/{\lambda }_{\mathrm{L}}^2$.

(20)${\lambda }_{\mathrm{L}}={\kappa }_{\mathrm{L}}\stackrel{\wedge }{{\sigma }_{w{z}_{\mathrm{M}}}}$

(21)$\stackrel{\wedge }{{\sigma }_{w{z}_{\mathrm{M}}}}={\left\{N{\displaystyle \sum _{i=1}^N\left[w_i{\left(\stackrel{\wedge }{a_1}+\stackrel{\wedge }{a_2}{x}_{\mathrm{M}i}+\stackrel{\wedge }{a_3}{y}_{\mathrm{M}i}-z_{\mathrm{M}i}\right)}^2\right]/\left[\left(N-3\right){\displaystyle \sum _{i=1}^N{w}_i}\right]}\right\}}^{1/2}$

where $\stackrel{\wedge }{{\sigma }_{w{z}_{\mathrm{M}}}}$ is the estimated value of the standard deviation of the weighted observation value z_Mi and ${\displaystyle \sum _{i=1}^N{w}_i}/N$ is the normalized mean weight value.

For the first round of fitting, i.e., the initial model parameter estimation method with equal weight, the reader can refer to Appendix A. Then, w_i can be related to Δ_i and $\stackrel{\wedge }{{\sigma }_{w{z}_{\mathrm{M}}}}$ in the second round. The number of measuring points N for each flange is less than 100 to ensure measurement efficiency. Then, the high value of $\stackrel{\wedge }{{\sigma }_{w{z}_{\mathrm{M}}}}$ causes λ_L to become too high. It is difficult to distinguish observations with low $\left|{\Delta }_i\right|$ from the measuring points. Therefore, it is necessary to re-estimate the weight values after obtaining the mathematical model in the last round of fitting. One should then repeat the iterations between point cloud fitting and weight estimation until the changes in the model parameters are small enough.

3.2.2. Standard Residual Singular Value Detection Based on the Chauvenet Criterion

If there is a singular value in the measuring points, its weight needs to be zeroed:

(22)$w_i=\{\begin{array}{c}0,\mathrm{if}\left|{\Delta }_i\right|\ge {\lambda }_0;\\ w_i,\mathrm{if}\left|{\Delta }_i\right|<{\lambda }_0\end{array}$

Singular points are found by calculating thresholds. Under the LS framework, the observed mean $\overline{z_{\mathrm{M}i}}$ in the standard deviation formula corresponds to the estimated value. There are three elements in the parameter vector a of the flange end face. The standard deviation $\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}$ is estimated by Equation (23).

(23)$\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}={\left[{\displaystyle \sum _{i=1}^N{\left(\stackrel{\wedge }{a_1}+\stackrel{\wedge }{a_2}{x}_{\mathrm{M}i}+\stackrel{\wedge }{a_3}{y}_{\mathrm{M}i}-z_{\mathrm{M}i}\right)}^2/\left(N-3\right)}\right]}^{1/2}$

Threshold λ₀ is set as Equation (24).

(24)${\lambda }_0={\kappa }_0\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}$

The Chauvenet criterion defines the probability of observed values occurring within the mean bandwidth range. The probability P_in of Δ_i occurring within the interval $[-{\kappa }_0\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}},{\kappa }_0\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}]$ is shown in Equation (25).

(25)$P_{\mathrm{in}}=1-P_{\mathrm{out}}={\displaystyle {\int }_{-{\kappa }_0\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}}^{{\kappa }_0\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}}{\mathrm{e}}^{-{\left({\Delta }_i/\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}\right)}^2/2}/\left(\sqrt{2\mathsf{\pi }}\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}\right)}\mathrm{d}{\Delta }_i=\mathrm{erf}\left(\frac{{\kappa }_0}{\sqrt{2}}\right)$

where erf is the Gauss error function.

Set the nominal number of occurrences of singular values to v₀ = NP_out. The relationship between the coefficient к₀ and the number of measuring points N is shown in Equation (26).

(26)${\kappa }_0=\sqrt{2}\mathrm{erfinv}\left(\frac{N-v_0}{N}\right)$

where erfinv is the Gauss inverse error function.

3.2.3. Clustering Detection

The principle of clustering detection requires that qualified observations be classified into a one-dimensional category. Their absolute deviations ${\Lambda }_i=\left|{\Delta }_i\right|$ are relatively concentrated. Simultaneously, the deviations of singular values are far away from this category. Thus, the specified interval value is searched for in the deviation sorting and is used as the criterion for singular value determination. Cluster detection can avoid misjudgment caused by the sparsity of observation distribution.

The absolute deviations Λ_i are arranged in ascending order and renumbered ordinal number n:

(27)$\cdot \cdot \cdot \le {\Lambda }_j\le \cdot \cdot \cdot \le {\Lambda }_k\le \cdot \cdot \cdot \le {\Lambda }_l\le \cdot \cdot \cdot$

(28)${\Lambda }_{\mathrm{s}}[1]\le \cdot \cdot \cdot \le {\Lambda }_{\mathrm{s}}[n]\le {\Lambda }_{\mathrm{s}}[n+1]\le \cdot \cdot \cdot \le {\Lambda }_{\mathrm{s}}[N]$

The difference of absolute deviations between adjacent orders is:

(29)$d[n+1]={\Lambda }_{\mathrm{s}}[n+1]-{\Lambda }_{\mathrm{s}}[n]$

There is a significant gap between the absolute deviations of singular values and the absolute deviations before their orders. The global criterion is that the difference of absolute deviations is higher than the typical threshold, i.e., d[n] ≥ к₁d_glob[n]. The local criterion is that the difference in absolute deviations is higher than the threshold generated by partial fore correlation values, i.e., d[n] ≥ к₂d_loc[n]. The threshold λ₀ is set as Equation (30).

(30)${\lambda }_0={\Lambda }_{\mathrm{s}}[n],\mathrm{if}\frac{d[n]}{d_{\mathrm{loc}}[n]}=\max$

The global threshold d_glob[n] and local threshold d_loc[n] are defined as the weighted average of specific differences, shown as Equations (31) and (32). The weight decreases as the absolute deviation of the current measuring point increases and the slope is related to the number of measurement points N. The d_loc[n] mainly depends on the influence of the differences of the adjacent absolute deviations. Local changes in the differences are overestimated when the absolute deviations are equal. d_loc[n] should be adjusted using Equation (33). Clustering detection is self-adaptive for both d_glob[n] and d_loc[n] change with d[n].

(31)$d_{\mathrm{glob}}[n]={\displaystyle \sum _{j=1}^{n-1}d[n-j]{\mathrm{e}}^{-2{\left(j/N\right)}^2}}/{\displaystyle \sum _{j=1}^{n-1}{\mathrm{e}}^{-2{\left(j/N\right)}^2}}$

(32)$d_{\mathrm{loc}}[n]={\displaystyle \sum _{j=1}^{n-1}d[n-j]{\mathrm{e}}^{-72{\left(j/N\right)}^2}}/{\displaystyle \sum _{j=1}^{n-1}{\mathrm{e}}^{-72{\left(j/N\right)}^2}}$

(33)$d_{\mathrm{loc}}[n]=\{\begin{array}{c}d[n]/{\kappa }_2,\mathrm{if}\mathrm{count}\left\{d[n-j]=0\right\}>\left(n-1\right)/2\\ d_{\mathrm{loc}}[n],\mathrm{if}\mathrm{count}\left\{d[n-j]=0\right\}\le \left(n-1\right)/2\end{array}$

3.3. Optimized Fitting Method for Spigot Cylindrical Face

The coordinates of the projection points M_{proj_i} on the fitted end face of the spigot measuring point M_i are calculated as Equation (34).

(34)$\left[\begin{array}{c}{x}_{\mathrm{proj\_M}i}\\ y_{\mathrm{proj\_M}i}\\ z_{\mathrm{proj\_M}i}\end{array}\right]=\left[\begin{array}{c}{x}_{\mathrm{M}i}-\frac{a_2\left(a_1+a_2{x}_{\mathrm{M}i}+a_3{y}_{\mathrm{M}i}-z_{\mathrm{M}i}\right)}{a_2^2+a_3^2+1}\\ y_{\mathrm{M}i}-\frac{a_3\left(a_1+a_2{x}_{\mathrm{M}i}+a_3{y}_{\mathrm{M}i}-z_{\mathrm{M}i}\right)}{a_2^2+a_3^2+1}\\ z_{\mathrm{M}i}+\frac{a_1+a_2{x}_{\mathrm{M}i}+a_3{y}_{\mathrm{M}i}-z_{\mathrm{M}i}}{a_2^2+a_3^2+1}\end{array}\right]$

To obtain the spatial circle center O_{proj_M}(b₁,b₂,b₃), the fitting equation after mapping is shown as Equation (35).

(35)$\{\begin{array}{c}{z}_{\mathrm{proj\_M}}{}_i=\sqrt{b_4{}^2-{\left(x_{\mathrm{proj\_M}}{}_i-b_1\right)}^2-{\left(y_{\mathrm{proj\_M}}{}_i-b_2\right)}^2}+b_3\\ z_{\mathrm{proj\_M}i}=a_1+a_2{x}_{\mathrm{proj\_M}i}+a_3{y}_{\mathrm{proj\_M}i}\end{array}$

Compared with the traditional LS method, an optimized total least squares (TLS) method is proposed to minimize the square of the distance between the mapping of observation points and the tangent points. The tangent point is on the tangent line of the model and is orthogonal to the minimum distance vector of the measuring point and the model. Thus, TLS has more geometric significance for fitting the mapped spatial circle of the spigot, as shown in Figure 7. Due to the symmetry of the circular curve, weight estimation and singular value detection are unnecessary, and observations with relatively higher errors can also be used as conditional inputs. Referring to the definition of χ²(b), the fitted model is transformed into $f\left(x_{\mathrm{proj\_M}i},y_{\mathrm{proj\_M}i},z_{\mathrm{proj\_M}i}|\mathit{b}\right)\equiv 0$ by shifting terms and the end face also provides a limiting condition function. The objective function of TLS is shown in Equation (36).

(36)$F_{\mathrm{TLS}}\left(\mathit{b}\right)={\displaystyle \sum _{i=1}^N{\left(\lambda -b_4\right)}^2}\to \mathrm{Min}$

(37)$\lambda ={\left[{\left(x_{\mathrm{proj\_M}}{}_i-b_1\right)}^2+{\left(y_{\mathrm{proj\_M}}{}_i-b_2\right)}^2+{\left(a_1+a_2{x}_{\mathrm{proj\_M}i}+a_3{y}_{\mathrm{proj\_M}i}{}_i-b_3\right)}^2\right]}^{1/2}$

Different from the linear approximation solution of the minimization equation of the traditional LS method, the minimization equation of TLS is a nonlinear function regarding the model parameter vector b = (b₁,b₂,b₃,b₄)^T. b can be solved through iterative compensations, and the initial estimated values in b can be calculated through the average method, as shown in Appendix A. The partial derivatives of the parameters are shown in Equations (38)–(41). The Jacobian matrix J_LS contains b and requires iterative operations. Simultaneously, the operations are also constrained by the end face model, i.e., Δb₃ = a₁ + a₂b_1,fit + a₃b_2,fit − b_3,init.

(38)$\partial f/\partial b_1=\left(b_1-x_{\mathrm{proj\_M}}\right)/\lambda$

(39)$\partial f/\partial b_2=\left(b_2-y_{\mathrm{proj\_M}}\right)/\lambda$

(40)$\partial f/\partial b_3=\left(b_3-a_1-a_2{x}_{\mathrm{proj\_M}i}-a_3{y}_{\mathrm{proj\_M}i}\right)/\lambda$

(41)$\partial f/\partial b_4=-1$

4. Measurement Experiment

Measurement and fitting experiments for an aeroengine low-pressure turbine shaft are performed. The three-axis measurement system includes: measuring head RENISHAW OMP 40-2 (a high-precision three-dimensional contact probe that can be used on coordinate measuring machines, with unidirectional repeatability of 1 μm and can be equipped with various types of measuring needles) and optical interface OMI-2T (an anti-interference infrared optical receiver for probe signals, with a detection range of 6 m). The initial position and pose of the workpiece in each experiment is randomized using an adjustment fixture. The measuring points are arranged according to the proposed planning schemes. The traditional fitting methods and the proposed optimized fitting methods are used to process the same set of measurement data, respectively. The singular point detection methods are verified and the results, such as goodness of fit and model parameter uncertainty, are compared.

4.1. Measurement and Fitting Experiment of Flange end Face

The flange end face measurement experiment is shown in Figure 8. In each measurement experiment of the rotor with random position and pose, five fitting schemes are used to fit the same coordinate data of the measurement points. The processes of different fitting schemes are shown in Figure 9.

Taking one experiment as an example, the coordinates of the measuring points are listed in Table 3. The initial estimated model parameter vector a_init is obtained after equal weight fitting and the compensation step size ∆a calculated based on the weighting times in different fitting schemes are listed in Table 4. All Δa_i ≤ 0.01 in the fourth iteration of weight estimation. Thus, the iteration is terminated.

The two methods for detecting singular points are compared. In the standard residual method, according to Chauvenet criterion, N = 28 makes к₀ = 2.785 when the nominal probability of singular value v₀ = 0.15. $\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}$ = 0.054 of the observed values after equal weight fitting. The threshold for this detection is λ₀ = 0.151. The residual $\left|{\Delta }_{19}\right|$ = 0.761 and $\left|{\Delta }_{20}\right|$ = 0.743 exceed λ₀. $\stackrel{\wedge }{{\sigma }_{z_{\mathrm{M}}}}$ = 0.064 and λ₀ = 0.179 after iterative reweighted fitting. Under this condition, $\left|{\Delta }_{19}\right|$ = 0.894 and $\left|{\Delta }_{20}\right|$ = 0.880 are also detected as exceeding λ₀. Then, both methods determine that the measuring points i = 19, 20 are singular. In cluster detection, when the equal weight fitting is used, the difference d[27] = 0.549 in absolute deviations of measuring point i = 19 (sorting number n = 27) exceeds the local threshold к₂d_loc[27] = 0.348 under coefficient к₂ = 2. This means that the points i = 19, 20 with greater absolute deviations show local anomalies. However, d[27] does not exceed the global threshold к₁d_glob[27] = 1.669 under к₁ = 17. Therefore, the local anomalies are considered to be caused by the sparsity of the observed value distribution and are nonsingular points. Similarly, singular points are not detected when using ordinary weighting. If the iterative reweighting fitting is used, d[27] = 0.775 of i = 19 (n = 27) exceeds к₁d_glob[27] = 0.611 under к₁ = 17 and к₂d_loc[27] = 0.163 under к₂ = 2. Therefore, this position is the boundary point for distinguishing singular categories, i.e., measuring points i = 19, 20 are singular.

Due to iterative reweighting, the weight of potential singular points that deviate significantly decreases. The fitting end face is closer to the measuring points with smaller deviation. Thus, the singular points can be easier to detect. The numbers of singular points of five measurement experiments with rotor random positions and poses are listed in Table 5. Whether the measuring points have better weights has a smaller impact on the standard residual method, but a greater impact on the clustering detection method. Thus, it is necessary to perform iterative estimation of weights and model parameters before detecting singular points. The two methods have different judgment criteria, but both have achieved good results in detecting singular points and can achieve better fitting accuracy. Thus, both are reasonable methods for detecting singular points and can be used in parallel when the numerical control system has sufficient computing power.

When the singular points are removed, the weight of the model parameter vector a is re-evaluated after four iterations of reweighting. Δa_i ≤ 0.01 is achieved after one compensation calculation. Thus, the iteration is terminated. The fitting accuracy results of the experimental example are listed in Table 6. Because the standard residual method and clustering detection method have removed the same singular points in this experimental example, the corresponding fitting accuracy results are equal. It can be seen that the weighted fitting obtains a goodness of fit, g_fit, that is less than 1.000 and closer to 1.000 than the equal weight fitting. Thus, the fitting model has a better consistency with the measuring point data. Simultaneously, the uncertainties of model parameters δ_a1, δ_a2 and δ_a3 decrease to varying degrees and iterative reweighting further improves fitting accuracy results compared with the ordinary weighting. The fitting after singular value detection further reduces the uncertainty of the model parameters in a small range and all uncertainties δ_a are less than 0.027. The fitting accuracy results of other experiments show similar patterns to this experiment example. In summary, “iterative reweighting + singular point detection (using residual method or cluster detection method) + iterative reweighting” is the optimal fitting scheme for the annular end face of the rotor flange.

4.2. Measurement and Fitting Experiment of the Spigot Cylindrical Face

The experiment in which we measured the spigot cylindrical face is shown in Figure 10. Each spigot measurement experiment is performed after the corresponding flange measurement experiment. The same coordinate data of each group of measuring points is fitted through three fitting schemes, including “OLS fitting”, “TLS fitting”, and “Iterative TLS fitting”.

The coordinates of the measuring points are listed in Table 7. The measuring points are projected towards the fitted flange end face z_M = 758.4004 − 0.0243x_M + 0.0300y_M. The iterative compensation ∆b using LS and TLS fitting methods are calculated by the initial estimate model parameter vector b_init, which are listed in Table 8. In the iterative TLS fitting, the fourth iteration is terminated for all Δb_i ≤ 0.01.

TLS is a nonlinear fitting, and the Jacobian matrix J_LS contains the model parameter vector b itself. The orders of magnitude of the diagonal element calculation result of covariance matrix C are too large and distort model the parameter uncertainty δ_b. With regard to the goodness of fit g_fit, b is gathered to the same side of the model equation in TLS method and the observation value defaults to 0. This causes the calculation results of g_fit to become too small to have significance in their evaluation. Therefore, the average orthogonal distance from the mapping point to the fitting model $\overline{d_{\mathrm{TLS}}}$ = f_TLS(b)/N is used as the main accuracy evaluation indicator. The fitting accuracy results of different fitting methods of this experiment example are listed in Table 9. TLS fitting achieves smaller average orthogonal distance $\overline{d_{\mathrm{TLS}}}$ compared with the traditional LS fitting and $\overline{d_{\mathrm{TLS}}}$ can be further reduced after iterative calculation of ∆b. The fitting accuracy results of other experiments show similar patterns to this experiment example.

In summary, the iterative TLS fitting method is the optimal scheme for the rotor spigot. The final fitting equation is calculated as Equation (42).

(42)$\{\begin{array}{c}{}^{0}z_{\mathrm{M}}=\sqrt{{119.3479}^2-{\left({}^{0}x_{\mathrm{M}}-630.3006\right)}^2-{\left({}^{0}y_{\mathrm{M}}+25.7087\right)}^2}+742.3128\\ {}^{0}z_{\mathrm{M}}=758.4004-0.0243{}^{0}x_{\mathrm{M}}+0.0300{}^{0}y_{\mathrm{M}}\end{array}$

Similarly, the TLS fitting method is used to obtain the model parameter vector of the flange reference hole: b_S = [769.8345, −28.0377, 738.8523, 10.1070]. The homogeneous position and pose matrix of the low-pressure turbine shaft are shown as Equation (43).

(43)$\mathit{T}=\left[\begin{array}{cccc}0.9996& 0.0174& 0.0243& 630.3006\\ -0.0167& 0.9994& -0.0300& -25.7087\\ -0.0248& 0.0296& 0.9993& 742.3128\\ 0& 0& 0& 1\end{array}\right]$

5. Conclusions

This study focuses on the precise measurement of the position and the pose of the flange and spigot structures of rotors. The optimized measuring point planning and fitting methods are proposed and then verified through experiments. The following conclusions are drawn:

• (1). The equation of model parameter uncertainty δ_a1 that affects the end face normal vector error angle and the equation of coordinate uncertainty δ_O of the spigot spatial circle center are derived. The effects of the circumferential phase angle θ_i, distribution range angle φ, total number of measuring points N, and multiple concentric measurement paths on uncertainty are analyzed. The measuring point planning schemes for flange annular end faces and spigot cylindrical faces of the LPT rotors of an aeroengine are proposed.

• (2). The LS minimization criterion and Gauss Newton iterative solution principle under Levenberg Marquardt damping factors are clarified. Based on this principle, an optimized iterative reweighting method with singular values removed is proposed for fitting the flange annular end face, including processes such as residual-based weight w_i iterative estimation, standard residual singular value detection under Chauvenet criterion, and clustering detection. Additionally, a mapping point TLS method with practical geometric significance is proposed for fitting the circle center of the spigot.

• (3). Measurement and fitting experiments on a low-pressure turbine shaft are performed. Five fitting methods for the flange annular end face are compared. When the optimized iterative reweighting method is used with singular values removed, the goodness of fit g_fit = 0.756, model parameter uncertainty δ_a1 = 0.026, and δ_a2 and δ_a3 < 4.3 × 10⁻⁵. Simultaneously, three fitting methods for the spigot cylindrical face are compared. When the mapping point iterative TLS fitting method is used, the goodness of fit g_fit = 0.011 and the average orthogonal distance $\overline{d_{\mathrm{TLS}}}$ = 0.089 mm. The proposed fitting methods achieve better fitting accuracy results compared with the other traditional methods.

Author Contributions

Conceptualization, T.Z.; methodology, T.Z. and H.G.; software, T.Z., L.H. and J.C.; validation, T.Z.; formal analysis, T.Z.; investigation, T.Z.; resources, T.Z. and X.J.; data curation, T.Z., L.H., X.J. and T.L.; writing—original draft preparation, T.Z.; writing—review and editing, L.H., X.J., T.L., Y.Z. and J.C.; visualization, T.Z.; supervision, H.G.; project administration, H.G. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Schematic of the measurement process of the aeroengine rotor flange and spigot structures.

Figure 2. Schematic of measuring point distribution schemes: (a) simplified assumption of three points distributed on the flange end face; (b) points distributed uniformly on single ring on the flange end face; (c) points distributed uniformly on double rings on the flange end face; (d) simplified assumption of three points distributed on the spigot cylindrical face; (e) points distributed uniformly on the spigot cylindrical face.

Figure 3. Relationship between measuring point distribution parameters and measurement uncertainty: (a) the influence of the number of flange measuring points and range angle; (b) the influence of the number of distribution rings of the flange measuring points; (c) the influence of the circumferential phase angle of the spigot measuring point; (d) the influence of the number of spigot measuring points and range angle.

Figure 3. Relationship between measuring point distribution parameters and measurement uncertainty: (a) the influence of the number of flange measuring points and range angle; (b) the influence of the number of distribution rings of the flange measuring points; (c) the influence of the circumferential phase angle of the spigot measuring point; (d) the influence of the number of spigot measuring points and range angle.

Figure 4. Schematic of the measurement faces on the LPT rotors of an aeroengine.

Figure 5. Main process of LS fitting and solution method.

Figure 6. The process of optimized iterative reweighted LS fitting method with singular values removed.

Figure 7. Geometric representation of the minimization target of the mapping point TLS method.

Figure 8. Measurement experiment for the rotor flange: (a) measurement of flange annular end face; (b) measurement trajectory.

Figure 9. Five fitting schemes in flange measurement experiments.

Figure 10. Measurement experiment for the rotor spigot: (a) measurement of spigot cylindrical face; (b) measurement trajectory.

Table 1

Optimized measuring point planning schemes of all end faces on the LPT rotors of an aeroengine.

End Faces	Fore Flange Face on the Rear Shaft/Aft Flange Face on the Cone Shaft	Fore Flange Face on the Cone Shaft/Aft Flange Face on the Mid Shaft	“Fifth” Bearing Positioning End Face on the Rear Shaft	Labyrinth Disc I Positioning End Face on the Cone Shaft	Labyrinth Disc II Positioning End Face on the Mid Shaft	Toothed Coupling Positioning End Face on the Mid Shaft
Total number of measuring points N	22	28	14	38	22	12
Distribution rings	2	2	1	2	2	1
Measuring path radius RC1 (mm)	104	125	79	191	107	64
Measuring path radius RC2 (mm)	113	149	-	208	120	-

Table 2

Optimized measuring point planning schemes of all cylindrical faces on the LPT rotors of an aeroengine.

Cylindrical Faces	Fore Spigot Face on the Rear Shaft/Aft Spigot Face on the Cone Shaft	Fore Spigot Face on the Cone Shaft/Aft Spigot Face on the Mid Shaft	“Fifth” Bearing Positioning Cylindrical Face on the Rear Shaft	Labyrinth Disc I Positioning Cylindrical Face on the Cone Shaft	Labyrinth Disc II Positioning Cylindrical Face on the Mid Shaft	Gear Coupling Positioning Cylindrical Face on the Mid Shaft
Total number of measuring points N	30	36	22	50	36	18
Spigot radiusRS (mm)	99	119	76	213	125	62

Table 3

Coordinates of the measuring points on the flange end face.

Point Ordinal i	Coordinate (xMi, yMi, zMi) (mm)	Point Ordinal i	Coordinate (xMi, yMi, zMi) (mm)
1	(779.131, −28.201, 738.614)	15	(755.145, −27.781, 739.189)
2	(765.561, 36.669, 740.852)	16	(743.757, 26.638, 741.100)
3	(725.182, 89.201, 743.466)	17	(709.884, 70.707, 743.266)
4	(665.992, 118.993, 745.763)	18	(660.228, 95.703, 745.333)
5	(599.710, 120.149, 747.419)	19	(604.624, 96.671, 747.502)
6	(539.475, 92.431, 748.026)	20	(554.089, 73.419, 748.019)
7	(497.208, 41.339, 747.523)	21	(518.630, 30.555, 746.707)
8	(481.281, −23.014, 745.988)	22	(505.269, −23.431, 745.412)
9	(494.854, −87.879, 743.753)	23	(516.656, −77.848, 743.528)
10	(535.232, −140.412, 741.160)	24	(550.530, −121.921, 741.345)
11	(594.424, −170.208, 738.866)	25	(600.188, −146.915, 739.493)
12	(660.704, −171.360, 737.185)	26	(655.788, −147.884, 738.012)
13	(720.941, −143.646, 736.614)	27	(706.324, −124.632, 737.485)
14	(763.209, −92.552, 737.063)	28	(741.783, −81.770, 737.861)

Table 4

Initial estimated values and iterative compensation values of the model parameters.

Iterations	Initial Values of Model Parameters andIterative Compensation Values
Equal weight	ainit = [758.7075, −0.0247, 0.0307]
First weighting	∆a = [−0.2475, 0.0003, −0.0006]
Second iteration reweighting	∆a = [−0.0677, 0.0001, −0.0001]
Third iteration reweighting	∆a = [−0.0203, 0.0000, 0.0000]
Fourth iteration reweighting	∆a = [0.0046, 0.0000, 0.0000]
Reweighting after removing singular values	∆a = [0.0001, 0.0000, 0.0000]

Table 5

Number of singular points under the different pre-processing and detection methods.

Pre-Processing before Detection	Detection Method	Experiment No.
		1	2	3	4	5
Equal weight	Residual method	2	1	0	2	0
	Clustering detection	0	0	0	1	0
Iterative reweighting	Residual method	2	1	1	2	0
	Clustering detection	2	1	0	2	0

Table 6

Fitting accuracy results of different fitting schemes for a rotor flange.

Fitting Scheme	Weighted Variance χ2	Goodness of Fit gfit	Uncertainty δa1	Uncertainty δa2	Uncertainty δa3
Equal weight	1.355	0.054	0.289	4.45 × 10−4	4.45 × 10−4
Initial equal weight + ordinary weighting	8.592	0.344	0.066	9.66 × 10−5	1.09 × 10−4
Initial equal weight + iterative reweighting	23.244	0.930	0.028	4.35 × 10−5	4.56 × 10−5
Iterative reweighting + singular point detection + reweighting	18.895	0.756	0.026	4.09 × 10−5	4.26 × 10−5

Table 7

Coordinates of measuring points on the spigot cylindrical face.

Point Ordinal i	Coordinate (xMi, yMi, zMi) (mm)	Point Ordinal i	Coordinate (xMi, yMi, zMi) (mm)
1	(749.407, −27.570, 735.914)	19	(510.816, −23.424, 741.800)
2	(747.988, −6.840, 736.496)	20	(512.268, −44.161, 741.056)
3	(742.961, 13.347, 737.306)	21	(517.289, −64.343, 740.337)
4	(734.499, 32.334, 738.070)	22	(525.729, −83.333, 739.604)
5	(722.886, 49.592, 738.795)	23	(537.342, −100.582, 738.873)
6	(708.472, 64.551, 739.561)	24	(551.787, −115.522, 738.089)
7	(691.649, 76.746, 740.451)	25	(568.603, −127.775, 737.305)
8	(672.951, 85.874, 741.211)	26	(587.287, −136.891, 736.484)
9	(652.969, 91.608, 741.744)	27	(607.274, −142.623, 735.869)
10	(632.299, 93.777, 742.380)	28	(627.966, −144.795, 735.212)
11	(611.534, 92.340, 742.795)	29	(648.684, −143.330, 734.757)
12	(591.375, 87.282, 743.212)	30	(668.887, −138.297, 734.422)
13	(572.349, 78.825, 743.408)	31	(687.876, −129.845, 734.323)
14	(555.123, 67.192, 743.357)	32	(705.121, −118.199, 734.216)
15	(540.129, 52.775, 743.428)	33	(720.110, −103.756, 734.224)
16	(527.913, 35.938, 743.246)	34	(732.342, −86.929, 734.493)
17	(518.756, 17.244, 742.812)	35	(741.472, −68.230, 734.762)
18	(513.015, −2.741, 742.426)	36	(747.227, −48.238, 735.204)

Table 8

Initial estimated values and iterative compensation values of the model parameters.

Fitting Scheme	Initial Values of Model Parameters andIterative Compensation Values
Initial average method estimation	binit = [630.2059, −25.6046, 742.3183, 119.3481]
LS fitting	∆b= [−0.1181, 0.1368, 0.0069, −0.0001]
First TLS fitting	∆b = [−0.1057, 0.1260, 0.0063, 0.0000]
Second TLS fitting	∆b = [0.1119, −0.1387, −0.0069, −0.0001]
Third TLS fitting	∆b = [0.0880, −0.1004, −0.0051, 0.0000]
Fourth TLS fitting	∆b = [0.0005, 0.0090, 0.0002, −0.0001]

Table 9

Accuracy results of different fitting schemes for rotor spigot.

Fitting Scheme	Average Orthogonal$\mathbf{Distance}\overline{{\mathit{d}}_{\mathbf{TLS}}}$ (mm)	Variance χ2	Goodness of Fit gfit
LS fitting	0.116	0.622	0.019
Single TLS fitting	0.105	0.517	0.016
Iterative TLS fitting	0.089	0.366	0.011

Appendix A

The estimated initial model parameters of the flange end face can be calculated as follows:(A1)$$a_{1,\mathrm{init}}=\left({\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{z}_{\mathrm{M}jk}}}-a_{2,\mathrm{init}}{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{x}_{\mathrm{M}jk}}}-a_{3,\mathrm{init}}{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{y}_{\mathrm{M}jk}}}\right)/N$$ (A2)$$a_{2,\mathrm{init}}=\left(S_{12}{S}_{23}-S_{13}{S}_{22}\right)/\left(S_{12}^2-S_{11}{S}_{22}\right)$$ (A3)$$a_{3,\mathrm{init}}=\left(S_{12}{S}_{13}-S_{11}{S}_{23}\right)/\left(S_{12}^2-S_{11}{S}_{22}\right)$$ (A4)$$S_{11}={\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{x}_{\mathrm{M}jk}^2}}-{\left({\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{x}_{\mathrm{M}jk}^{}}}\right)}^2/N$$ (A5)$$S_{12}={\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}\left(x_{\mathrm{M}jk}^{}{y}_{\mathrm{M}jk}^{}\right)}}-{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{x}_{\mathrm{M}jk}^{}}}·{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{y}_{\mathrm{M}jk}^{}}}/N$$ (A6)$$S_{12}={\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}\left(x_{\mathrm{M}jk}^{}{z}_{\mathrm{M}jk}^{}\right)}}-{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{x}_{\mathrm{M}jk}^{}}}·{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{z}_{\mathrm{M}jk}^{}}}/N$$ (A7)$$S_{22}={\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{y}_{\mathrm{M}jk}^2}}-{\left({\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{y}_{\mathrm{M}jk}^{}}}\right)}^2/N$$ (A8)$$S_{23}={\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}\left(y_{\mathrm{M}jk}^{}{z}_{\mathrm{M}jk}^{}\right)}}-{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{y}_{\mathrm{M}jk}^{}}}·{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^{N/J}{z}_{\mathrm{M}jk}^{}}}/N$$

The estimated initial model parameters of the spigot cylindrical face can be calculated as follows:(A9)$$b_{1,\mathrm{init}}={\displaystyle \sum _{i=1}^N{x}_{\mathrm{proj\_M}i}/N}$$ (A10)$$b_{2,\mathrm{init}}={\displaystyle \sum _{i=1}^N{y}_{\mathrm{proj\_M}i}/N}$$ (A11)$$b_{3,\mathrm{init}}=a_1+a_2{b}_{1,\mathrm{init}}+a_3{b}_{2,\mathrm{init}}$$ (A12)$$b_{4,\mathrm{init}}={\displaystyle \sum _{i=1}^N\sqrt{{\left(x_{\mathrm{proj\_M}i}-b_{1,\mathrm{init}}\right)}^2+{\left(y_{\mathrm{proj\_M}i}-b_{2,\mathrm{init}}\right)}^2+{\left(a_1+a_2{x}_{\mathrm{proj\_M}i}+a_3{y}_{\mathrm{proj\_M}i}-b_{3,\mathrm{init}}\right)}^2}/N}$$

References

1. Liu, Y.; Zhang, M.; Sun, C.; Hu, M.; Chen, D.; Liu, Z.; Tan, J. A method to minimize stage-by-stage initial unbalance in the aero engine assembly of multistage rotors. Aerosp. Sci. Technol.; 2019; 85, pp. 270-276. [DOI: https://dx.doi.org/10.1016/j.ast.2018.12.007]

2. Yang, Z.; Hussain, T.; Popov, A.A.; McWilliam, S. Novel optimization technique for variation propagation control in an aero-engine assembly. Proc. Inst. Mech. Eng. Part B J. Eng. Manuf.; 2011; 225, pp. 100-111. [DOI: https://dx.doi.org/10.1243/09544054JEM2043]

3. Nizametdinov, F.R.; Romashin, Y.S.; Berne, A.L.; Leontyev, M.K. Investigation of bending stiffness of gas turbine engine rotor flanged connection. J. Mech.; 2020; 36, pp. 729-736. [DOI: https://dx.doi.org/10.1017/jmech.2020.14]

4. Zhou, T.; Gao, H.; Wang, X.; Li, L.; Chen, J.; Peng, C. Prediction method of aeroengine rotor assembly errors based on a novel multi-axis measuring and connecting mechanism. Machines; 2022; 10, 387. [DOI: https://dx.doi.org/10.3390/machines10050387]

5. Ding, S.; He, Y.; Zheng, X. A probabilistic approach for three-dimensional variation analysis in aero-engine rotors assembly. Int. J. Aeronaut. Space; 2021; 22, pp. 1092-1105. [DOI: https://dx.doi.org/10.1007/s42405-021-00384-w]

6. Hussain, T.; Yang, Z.; Popov, A.A.; McWilliam, S. Straight-build assembly optimization: A method to minimize stage-by-stage eccentricity error in the assembly of axisymmetric rigid components (two-dimensional case study). J. Manuf. Sci. Eng.; 2011; 133, 031014. [DOI: https://dx.doi.org/10.1115/1.4004202]

7. Li, Z.; Zheng, X. Review of design optimization methods for turbomachinery aerodynamics. Prog. Aerosp. Sci.; 2017; 93, pp. 1-23. [DOI: https://dx.doi.org/10.1016/j.paerosci.2017.05.003]

8. Jin, S.; Ding, S.; Li, Z.; Yang, F.; Ma, X. Point-based solution using Jacobian-Torsor theory into partial parallel chains for revolving components assembly. J. Manuf. Syst.; 2018; 46, pp. 46-58. [DOI: https://dx.doi.org/10.1016/j.jmsy.2017.11.003]

9. Ma, Y.; Wang, Y.; Wang, C.; Hong, J. Interval analysis of rotor dynamic response based on Chebyshev polynomials. Chin. J. Aeronaut.; 2020; 33, pp. 2342-2356. [DOI: https://dx.doi.org/10.1016/j.cja.2020.04.010]

10. Satish, T.N.; Rao, A.N.V.; Nambiar, A.S.; Uma, G.; Umapathy, M.; Chandrashekhar, U.; Petley, V.U. Investigation into the development and testing of a simplex capacitance sensor for rotor tip clearance measurement in turbo machinery. Exp. Tech.; 2018; 42, pp. 575-592. [DOI: https://dx.doi.org/10.1007/s40799-018-0281-y]

11. Jayaweera, N.; Webb, P.; Johnson, C. Measurement assisted robotic assembly of fabricated aero-engine components. Assem. Autom.; 2010; 30, pp. 56-65. [DOI: https://dx.doi.org/10.1108/01445151011016073]

12. Wong, S.K.; Riseborough, E.; Duff, G.; Chan, K.K. Radar cross-section measurements of a full-scale aircraft duct/engine structure. IEEE Trans. Antennas Propag.; 2006; 54, pp. 2436-2441. [DOI: https://dx.doi.org/10.1109/TAP.2006.879223]

13. Garcia, I.; Przysowa, R.; Amorebieta, J.; Zubia, J. Tip-clearance measurement in the first stage of the compressor of an aircraft engine. Sensors; 2016; 16, 1897. [DOI: https://dx.doi.org/10.3390/s16111897]

14. Ristic, M.; Brujic, D.; Ainsworth, I. Measurement-based updating of turbine blade CAD models: A case study. Int. J. Comput. Integr. Manuf.; 2004; 17, pp. 352-363. [DOI: https://dx.doi.org/10.1080/0951192032000159148]

15. Burghardt, A.; Kurc, K.; Szybicki, D.; Muszynska, M.; Szczech, T. Robot-operated inspection of aircraft engine turbine rotor guide vane segment geometry. Teh. Vjesn. Tech. Gaz.; 2017; 242, pp. 345-348.

16. Zhang, M.; Liu, Y.; Sun, C.; Wang, X.; Tan, J. Measurements error propagation and its sensitivity analysis in the aero-engine multistage rotor assembling process. Rev. Sci. Instrum.; 2019; 90, 115003. [DOI: https://dx.doi.org/10.1063/1.5109199]

17. Sun, B.; Li, B. Laser displacement sensor in the application of aero-engine blade measurement. IEEE Sens. J.; 2016; 16, pp. 1377-1384. [DOI: https://dx.doi.org/10.1109/JSEN.2015.2497363]

18. Zhou, Z.; Liu, W.; Wu, Q.; Wang, Y.; Yu, B.; Yue, Y.; Zhang, J. A combined measurement method for large-size aerospace components. Sensors; 2020; 20, 4843. [DOI: https://dx.doi.org/10.3390/s20174843] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32867175]

19. Fu, G.; Sheng, B.; Chen, G.; Luo, R.; Sheng, G.; Lu, Y. A contact on-machine measurement points’ planning method based on similarity of 3-D process element. IEEE Trans. Instrum. Meas.; 2022; 71, pp. 1-14. [DOI: https://dx.doi.org/10.1109/TIM.2022.3212526]

20. Mao, Z.; Li, S.; Xu, Y.; Zeng, Q.; Zhu, K. A difference measurement points planning method for large-scale surface of aircraft. J. Beijing Univ. Aeronaut. Astronaut.; 2020; 46, pp. 1024-1031.

21. Magdziak, M. Selection of the best model of distribution of measurement points in contact coordinate measurements of free-form surfaces of products. Sensors; 2019; 19, 5346. [DOI: https://dx.doi.org/10.3390/s19245346] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31817193]

22. Piasecki, R. A generalization of the inhomogeneity measure for point distributions to the case of finite size objects. Phys. A; 2008; 387, pp. 5333-5341. [DOI: https://dx.doi.org/10.1016/j.physa.2008.05.034]

23. Calvo, R. Sphericity measurement through a new minimum zone algorithm with error compensation of point coordinates. Measurement; 2019; 138, pp. 291-304. [DOI: https://dx.doi.org/10.1016/j.measurement.2019.01.078]

24. Seiler, A.; Grossmann, D.; Juettler, B. Spline surface fitting using normal data and norm-like functions. Comput. Aided Geom. Des.; 2018; 64, pp. 37-49. [DOI: https://dx.doi.org/10.1016/j.cagd.2018.06.005]

25. Brujic, D.; Ainsworth, I.; Ristic, M. Fast and accurate NURBS fitting for reverse engineering. Int. J. Adv. Manuf. Technol.; 2011; 54, pp. 691-700. [DOI: https://dx.doi.org/10.1007/s00170-010-2947-1]

26. Gao, F.; Pan, Z.; Zhang, X.; Li, Y.; Duan, J. An adaptive sampling method for accurate measurement of aeroengine blades. Measurement; 2021; 173, 108531.

27. Magdziak, M. The influence of a number of points on results of measurements of a turbine blade. Aircr. Eng. Aerosp. Technol.; 2017; 89, pp. 953-959. [DOI: https://dx.doi.org/10.1108/AEAT-03-2016-0044]

28. Zhang, M.; Liu, Y.; Sun, C.; Wang, X.; Tan, J. Precision measurement and evaluation of flatness error for the aero-engine rotor connection surface based on convex hull theory and an improved PSO algorithm. Meas. Sci. Technol.; 2020; 31, 085006. [DOI: https://dx.doi.org/10.1088/1361-6501/ab8170]

29. Zhu, L.; Luo, H.; Ding, H. Optimal design of measurement point layout for workpiece localization. J. Manuf. Sci. Eng.; 2008; 131, pp. 011006-01100613. [DOI: https://dx.doi.org/10.1115/1.3039515]

30. Ding, S.; Jin, S.; Li, Z.; Wei, Z.; Yang, F. Deviation propagation model and optimization of concentricity for aero-engine rotor assembly. J. Shanghai Jiao Tong Univ.; 2018; 52, pp. 54-62.

31. Xie, Z.; Jiang, Z.; Zhu, H.; Li, Z.; Lei, L.; Kong, M.; Wang, D.; Zhang, B. Research on fitting method of thread middle axis based on 3D point cloud. Acta Metrol. Sin.; 2020; 41, pp. 918-925.

32. Whitehouse, D.J. Some theoretical aspects of error separation techniques in surface metrology. J. Phys. E Sci. Instrum.; 1976; 9, pp. 531-536. [DOI: https://dx.doi.org/10.1088/0022-3735/9/7/007]

33. Jin, S.; Fan, D.; Chen, R.; Tian, M. Method and experimental study of machine vision measurement for the roundness of rotary parts. Mach. Des. Res.; 2016; 32, 117–119, 124

34. Chetwynd, D.G. High-precision measurement of small balls. J. Phys. E Sci. Instrum.; 1987; 20, pp. 1179-1187. [DOI: https://dx.doi.org/10.1088/0022-3735/20/10/005]

35. Sun, C.; Wang, B.; Liu, Y.; Wang, X.; Li, C.; Wang, H.; Tan, J. Design of high accuracy cylindrical profile measurement model for low-pressure turbine shaft of aero engine. Aerosp. Sci. Technol.; 2019; 95, 105442. [DOI: https://dx.doi.org/10.1016/j.ast.2019.105442]

36. Sun, C.; Wang, L.; Tan, J.; Zhao, B.; Tang, Y. Design of roundness measurement model with multi-systematic error for cylindrical components with large radius. Rev. Sci. Instrum.; 2016; 87, 025110. [DOI: https://dx.doi.org/10.1063/1.4941679]

37. Zhang, M.; Liu, Y.; Sun, C.; Wang, X.; Tan, J. A systematic error modeling and separation method for the special cylindrical profile measurement based on 2-dimension laser displacement sensor. Rev. Sci. Instrum.; 2019; 90, 105006. [DOI: https://dx.doi.org/10.1063/1.5111350]

38. Wang, J.; Sun, Q.; Yuan, B. Novel on-machine measurement system and method for flatness of large annular plane. Meas. Sci. Technol.; 2020; 31, 015004. [DOI: https://dx.doi.org/10.1088/1361-6501/ab3a66]

39. Song, T.; Xi, F.J.; Guo, S.; Ming, Z.; Lin, Y. A comparison study of algorithms for surface normal determination based on point cloud data. Precis. Eng.; 2015; 39, pp. 47-55. [DOI: https://dx.doi.org/10.1016/j.precisioneng.2014.07.005]

40. Nouira, H.; Bourdet, P. Evaluation of roundness error using a new method based on a small displacement screw. Meas. Sci. Technol.; 2014; 25, 044012. [DOI: https://dx.doi.org/10.1088/0957-0233/25/4/044012]