A Probabilistic Approach for Three-Dimensional Variation Analysis in Aero-engine Rotors Assembly

Rotor assembly is a core tache in the whole process of aero-engine manufacturing. Preventing out-of-tolerance of concentricity is one of the primary tasks. Conventional assembly approaches are based on a manual test with the dial indicator, depending on experience appraises, which lack systematic and quantitative precision design theory. As a result, two issues need to be solved: the modeling problem of complicated geometric variations in three-dimensions, as well as the abnormal distribution of ubiquitous actual deviations. This work attempts to propose a novel probabilistic approach for three-dimensional variation analysis in rotor assembly. Based on rotor’s revolving characteristics and multistage stacking process, Jacobian–Torsor model is adopted to establish the variation propagation, and Pearson distribution family is used to derive the probability density function, which can quickly determine the variation distribution pattern and efficiently perform statistical variation analysis. A real case of mechanical assemblies consisting of revolving axisymmetric components is concerned. The results show that the suggested method has a similar accuracy, but much higher efficiency than conventional methods. Calculations agree with the experimentations, and the probability distribution type of the part’s variation has an appreciable impact on the final assembly precision.


Introduction
The core of aero-engine consists of a large number of rotor parts. As the high-speed rotating component, aero-engine rotor is usually produced as a combination of different rotor parts, to keep its repairability and maintainability. Due to the presence of part's manufacturing error, datum error, positioning error etc., variation accumulation is inevitable during aero-engine rotor assembly process. Minor variations are caused in the actual component from the nominal geometry [1]. As the assembling process went on, these variations would propagate and drive the final dimensions of assembly out of specification. Reducing product cost and enhancing quality is the tallest target that the enterprise pursues. Statistical variation analysis is an excellent way for predicting product performance in the initial design stage [2]. It is actually a quantitative design method which can analyze the size effect containing various tolerance specification quantitatively, to meet both design and manufacturing requirements.
Machinery quality is affected by many factors during processing, such as tool wearing, clamping way, human factors, etc. When various influencing factors are generated randomly, parts manufacturing variations are in accordance with normal distribution. Normality assumption is convenient for statistical evaluation due to its symmetry, and beneficial to theoretical study. On this account, traditional methods of tolerance design and variation analysis have widely adopted the hypothesis that the random variables of parts variations are normally distributed. Conventional methods for exploring assembly variation propagation are mainly based on a manual test with the dial indicator, depending on experience appraises, or extremum method [3][4][5][6], or the mean square root (MSR) method [7][8][9], which lack of systematic and quantitative precision design theory. However, the result calculated by extremum method is far too restrictive that would lead to high production cost, while the MSR method concerns mean-square variation that is primarily suitable for solving the plane size-chain problem [10][11][12]. These methods are used frequently in the case of normal distribution. Meanwhile, they are underdeveloped in regard to the geometrical errors in the spatial 3D state [7,13].
However, in practical matters, most variation distribution of part variable disobey the normal type. For example, if tool wear plays the dominant role during processing, the actual variable variation of part's dimension would follow uniform distribution (flat-topped distribution); if a manufacturing process involves the simultaneous processing of two machine tools, the mixing variable variations of the parts would satisfy bimodal distribution; when adopting a trial cutting approach to fabricate shaft journals or apertures, operators tend to keep a certain machining allowance deliberately, which would cause skewed distribution of the final dimensional variables. The most versatile and popular way for statistical variation analysis is Monte Carlo (MC) simulation, which is a numerical method and suitable for data abnormal distribution. According to known or assumed statistical distributions, random variation of each component can be generated, and the corresponding critical feature will be calculated for each set of part variables. By this means, the probability can be worked out indicating that whether the key characteristic is satisfied or not. Nowadays, MC simulation method has been developed into the core of most of the commercial software for 3D tolerance analysis, such as 3DCS, VSA, CATS-1DXL and CETOL, to name a few. However, to realize accurate prediction, a large number of data samples are desired [14][15][16][17], inevitably causing the time-consuming calculation, which is the obvious drawback to MC approach. This problem becomes more prominent when variation analysis is performed within an iterative loop of a more complicated variation propagation issue, and the solution would be extremely computationally expensive.
From the above, it is observed that the statistical methods developed so far are not very effective for solving the issue of 3D variation propagation with abnormal distribution. As this paper focuses on the general case of the mechanical assembly consisting of revolving axisymmetric components [18][19][20][21][22][23], which is widely used in the aero-engine industry, especially in the engine rotors fabricating [24,25], there is a clear need to adopt a simple and rapid probabilistic approach for determining the likelihood that the mechanical assembly is acceptable or not, so as to consider the practical distribution of variation variables and efficiently deal with large sets of data.
Recently, Yang et al. [26] researched the density distribution of rotation error in aero-engine assembly. Probability density functions (PDF) for the eccentricity were formulized based on the assumption that the component variations were statistically independent, zero-mean Gaussian random variables with known standard deviations. This PDF was essentially a Rayleigh distribution and results were calculated more efficiently than Monte Carlo simulations. Ghie [27] enriched the Jacobian-Torsor model combining with uncertain measurements. The uncertainty magnitudes were estimated using the approximate normal distribution method and spread to FR by linear superposition. However, this PDF could not give consideration to the skewed distribution of variables. Ahsanullah and Hamedani [28], on the basis of the conditional distribution of generalized order statistics, tried to study more detail on the various characterizations of certain univariate continuous distributions. It provided a base of further research of multivariable continuous distributions, which was applicable to the stream of multisource variation modelling. Kharoufehy and Chandraz [29] used a discretized, multivariate kernel density estimate and a simple transformation to approximate the probability distribution of the overall assembly characteristic. By this method, the typical assumptions of independence and normality could possibly be relaxed when the number of components being combined was small. But it was hard to select the smoothing parameter (h), which would yield possibly dramatically different distributions. From this review, we can know that modelling PDF is very complex and it is difficult to consider the polymorphic probability density distribution. In this paper, we will introduce an innovation of Pearson distribution family [30,31] to address this issue, which will be discussed detailedly in Sect. 2.2.
In addition, variation accumulation of the rotor assembling process can be hardly described by reason of the uncertainty inherent in a complex assembly. In the past few years, three-dimensional variation modelling methods have been studied. Four mainstream modelling approaches among them have been reported largely and studied deeply, i.e., the direct linearization method (DLM), the tolerance map (T-Map) model, the matrix model and the unified Jacobian-Torsor (J-T) model. The DLM [32,33] uses vectors to represent either component dimensions or assembly dimensions. Although all types of tolerances can be modeled, it is difficult to define the joint types and the effects of geometric variations, as well as the interaction between geometric tolerance and dimensional tolerance is indistinct. The T-Map model [34,35] is a hypothetical Euclidean volume that is fully consistent with the ASME standard, but the Minkowski operation of T-Map for tolerance propagation and the visualization of higher-dimensional maps are not straightforward. The matrix model [36,37] uses a displacement matrix to describe small displacements of a feature within its tolerance zone. Although it is suitable for tolerance propagation analysis, the computation process is very difficult and timeconsuming, under the action of large numbers of constraint inequalities, especially non-linear inequalities. The J-T theory [38,39] is well suited to a complex assembly that contains large numbers of joints and geometric tolerances.

3
It combines the advantages of the torsor model which is suitable for variation representation and the Jacobian matrix which is suitable for variation propagation. Therefore, the J-T model is well suited to a complex assembly that contains large numbers of joints and geometric deviations. Parts of aero-engine rotors are usually rich in variation information and joints of cylindrical and planar surfaces. In view of the rotor's configuration and assembling characteristics of aeroengine, the J-T model is used to build the variation propagation model of multi-stage stackup in this paper.
Therefore, based on the above aspects, this paper will propose a novel probabilistic analysis method combined with Jacobian-Torsor modelling theory, to take account of random component variability during revolving components assembly. The rest of the paper is organized as follows. Section 2 reminds the theory of the unified Jacobian-Torsor and Pearson distribution family. Section 3 explores the statistical analysis method and gives the specific solution. Section 4 proves this method by a realistic case originating in aero-engine subassembly. Section 5 discusses the results and Sect. 6 is conclusions.

Description of Jacobian-Torsor Model
The traditional J-T model consists of the Jacobian matrix and torsor model. The Jacobian matrix originates from the expression in robot kinematics, which is: represents the ith frame orientation; R Pti is a projection matrix that allows efficient transformation of variation direction to the target analysis direction; W n i describes the positional relation between adjacent frames in the form of skew-symmetric matrix. It is defined as Eq. (2): The subscript FE i in [J] means the ith functional element. Two categories of FEs exist in the assembly: internal FE (IFE) and contact FE (CFE). The former denotes two mating features on the same part and the latter denotes two mating features on separated ones. FEs' Variations could be delivered to the FR (functional requirement) by Jacobian matrix. It is observed that the Jacobian model only expresses the three-dimensional geometric relationship between the (1) features, without defining the objects and contents propagated, so it is an incomplete model. Torsor model can describe the ideal position and orientation of any surfaces by three translational vectors and three rotational vectors. Figure 1 shows a torsor model to express the position and orientation of an ideal surface. Assuming t is the tolerance between arbitrary feature S 1 and nominal feature S 0 , there is: where u, v, w denote translation vectors along x, y and z axes; α, β and γ denote rotation vectors around x, y and z axes respectively in local coordinate system.
Hervé [40] presented the displacements set theory, defining the non-invariant degrees, which determined the effective vectors. According to this, some vectors could be taken zero value to simplify the modeling process. Table 1 has shown seven primary torsor forms.
It is observed that the torsor model fully expresses the position and direction of the feature in 3D space. If combining with the tolerance constraint equation, a complete syntactic and semantic representations of geometric tolerances can be construed. But it is also observed that the torsor model is unable to characterize the variation propagation. It will probably have to be with the help of some transfer function or transfer matrix.
The variation of torsors propagates to the target characteristics taking advantage of Jacobian matrix, which accurately expresses the relative positions and orientations between the local and global coordinates. The unified J-T model takes synthesize advantages of the Jacobian matrix and torsor model. The former is appropriate for variation propagation and the latter is suited to variation expression [41,42]. The form of J-T model could be written as:

Fig. 1 Torsor model in its tolerance zone
where , is the variation zone, and β, γ, u, v and w follow the similar way as α. Other parameters have the same meanings as those in Eq. (1).

Pearson Distribution Family
Statistician K. Pearson [30,31] firstly proposed the Pearson system of distributions in 1895, which involved the normal distribution and diversified abnormal distributions. The Pearson system included four main types of parameters, which were used to reflect the general distribution in practice. Its probability density function (PDF) f(x) was determined by the differential equation as follows: When a 1 = 1, the parameters a 0 , b 0 , b 1 and b 2 entirely depend on the first four orders of the central moments of the distribution parameters μ 1 , μ 2 , μ 3 , and μ 4 . It also could be indicated as the moment ratios as β 1 = μ 3 2 /μ 2 3 and β 2 = μ 4 /μ 2 2 , which are shown as follows: where The general solution of f(x) is: Based on the above analysis, Pearson system will be elaborated as seven distribution patterns [43,44], which are: where p, q and α, β are all greater than zero, and α/p = β/q. It is essentially Beta distribution.
where p, > 0 . It is a special case of pattern I where p = q, = , and the distribution curve of pattern II is a symmetrical U type curve.
where , , > 0 . It is a shifted gamma distribution where is a location parameter.
Distribution pattern II: where , > 0 . It is a shifted inverse gamma distribution where is a location parameter.
w h e r e , > 0 .
where ≥ 1∕2, p > 0 . Student's t-distribution belongs to this kind of distribution. Besides these seven main types of distribution, Pearson distribution family also includes some special cases, such as normal distribution, power function distribution, exponential distribution, etc. Seven types of Pearson distribution can be identified according to the different distributions of the quadratic equation roots , the Pearson distribution pattern could be distinguished by the value of D and λ, as demonstrated in Table 2.
It can be also categorized in terms of moment ratios 1 , 2 . The in first approach of classification can be represented as The details are shown in Table 3.

Statistical Algorithm
In this paper, simplified cylindrical components were used to illustrate the assembling process, where coordinate system 0 represented the global coordinate system and the revolving axis was denoted with a virtual centerline that passed through the bottom part's base centre and was perpendicular to its datum plane, as illustrated in Fig. 2. Before applying the model, it should be noted that for every rotor part, two assumptions will be made to simplify operations.
• As the Jacobian matrix cannot solve the structure of partial parallel connection well, the locating spigot round (A) contact pair of the cylindrical joint) is ignored, and the flange plane (B) contact pair of the planar joint) of the rotor part is reserved. Generally, the contact pair of B shows the most significant effect due to that the mounting edge of area B is much larger than that of A. Figure 3a has illustrated this difference. • This paper adopts a unified solution that each rotor's bottom face is set as the datum, while the top face is designated by a geometric tolerance with the corresponding (18) Distribution pattern VI ∶ f (x) Figure 3b has illustrated the standardized rotor part to simplify and maintain the geometric tolerance information.
On the basis of the above assumptions, variations propagate from rotor's bottom to the top. According to the characteristics of Jacobian matrix, [[J 1 , J 2 , J 3 , J 4 , J 5 , 6 ] FEn ] represents the transfer process of the variations from FE 1 to FE n . Therefore, Eq. (4) using interval-based uncertainty zone representation can be further rewritten as Eq. (20), which is fit for deterministic variation analysis and beneficial to the error compensation.
Equation (20) could be characterized as y FR = f (u i , v i , w i , α i , β i , γ i ), (i = 1, 2, …, n), where y FR denotes the variation of FR, and (u i , v i , w i , α i , β i , γ i ) represents the translation variation or rotation variation of FE in different direction. Assuming that all the variation of FE is independent of each another and distributed randomly, we could calculate the four main parameters of FR in virtue of the parameters of FE. Based on the characteristics of Pearson system, the distribution pattern could be deduced and the PDF could be determined when the first four central moments are derived. As for arbitrary distribution, the μ 1 which is the first-order central moment will be equal to zero identically; the μ 2 which is the second-order central moment represents data dispersion degree, equaling to the variance σ 2 ; the μ 3 which is the third-order central moment is related to S (coefficient of skewness), and S = μ 3 /σ 3 ; and the μ 4 which is related to K (coefficient of kurtosis) is formulated as K = μ 4 /σ 4 . Due to the four main parameters (mean value, variance, skewness and kurtosis) of an arbitrary probability distribution could characterize the distribution shape visually, we will calculate the four main parameters to obtain the first four central moments indirectly in this article. Computational formulas of the mean value, variance, skewness and kurtosis can be consulted in reference [45].
Back in Eq. (20), hence a 0 , b 0 , b 1 and b 2 could be gained with the known distribution parameter of FR. Distribution pattern would be determined by adopting the method of Pearson distribution family (Tables 2 or 3) and the PDF could be expressed. Ultimately the probability distribution curve could be derived and the assembly qualified rate could be worked out. The solution procedure is shown in Fig. 4.

Case Study
To prove the feasibility of the new method mentioned above, a real rotor assembly originating in an aero-engine will be used for statistical variation analyzing. As shown in Fig. 5, the aero-engine rotor assembly consists of four revolving components with high precision and all relevant coordinates have been established. Frame 0 is the overall coordination system and the revolving axis is specified as the center line that passes through the datum plane of the bottom part and is perpendicular to it.
Every part's dimension and composite tolerance have been shown in Table 4, which have been illustrated in Fig. 5. The pose position and orientation of the whole assembly are directly affected by each rotor's geometric variations and their joints. Various errors accumulate on the final rotor part, which are very important for performance stability of aero-engine assembly. To meet the technical requirement, concentricity control of rotor assembly is very critical. On

Fig. 2
Assembly accuracy prediction: a poor assembly and b qualified assembly Fig. 3 The rotor parts: a the actual parts and b the geometric model the basis of statistical variation analysis approach that is referred in Sect. 3, the specific calculation process is demonstrated as follows. First, a coordinate reference system and the assembly relation graph were built. The global reference system 0 was established in the center of the rotor base, and local coordinate systems were set in the center of related adjacent surfaces. The assembly relation graph was illustrated in Fig. 6. FR was the pose accuracy of the top part in regard to the coordinate system 0. It was observed that there was a serial path for variation propagation from FEi to FR, i.e. IFE 1 -CFE 1 -IFE 2 -CFE 2 -IFE 3 -CFE 3 -IFE 4 -FR.
Secondly, according to the connection graph mentioned above, the variation propagation model would be established. Based on the expression of the unified Jacobian-Torsor model, the Jacobian matrix and the tolerancing torsor of FEi could be calculated, as well as its corresponding variation range could be determined. Table 5 has shown the detailed expression forms.
Next, according to Eq. (20), by multiplying the Jacobian matrixes and torsors together, the expression of the results could be derived as Eq. (21). In this case, the components' eccentricity along one certain direction was arbitrarily chosen for demonstration. For example, in the direction of x axis, the cumulative variation could be extracted from u 4 in Eq. (21) and formulized as Eq. (22). And it was also important to note that all the effective vectors in tolerancing torsors should satisfy the variation requirement, as listed in the third column of Table 5.
y u 4 = 0.02 1 + 0.04 2 + 0.06 3 + 0.08 4 .  where a i represented the sensitivity coefficient, E xi and E y respectively were the mathematical expectations of the ith FE and FR. The σ xi and σ y, respectively, were the standard deviations of the ith FE and FR. The S xi and S y denoted skewness coefficient of the ith FE and FR. The K xi and K y denoted the kurtosis coefficient of the ith FE and FR.
It should be noted that in abnormal distribution models, Beta distribution can handle various laws of the probability of manufacturing variations [8]. There are four parameters in the Beta distribution function, which are two location parameters and two shape parameters. According to different parameter values, it can characterize many different distribution patterns, such as Rayleigh distribution, uniform distribution, trapezium distribution, triangular distribution, etc. As a result, in this paper we would use Beta distribution to simulate skewed distribution of the actual manufacturing practice.
The most common types of probability distribution in the machining process mainly consist of normal distribution, skewed distribution, uniform distribution (flat-topped distribution), etc. [46]. To test the influence of variation distribution type on the final assembly precision, the distribution type of FE was supposed to follow the above distribution patterns, respectively. In the following, four probability combinations would be researched as listed in Table 6, where the skewed distribution would be simulated by Beta distribution approximation.
The distribution parameters of FR could be calculated based on Eqs. (23) to (26). The distribution results were shown in Table 7. According to Eqs. (6) to (11), the values of a 0 , b 0 , b 1 , b 2 could be further derived. Then referring to Tables 2 and 3, the variation distribution pattern of FR could be determined. The results had been listed in the last column of Table 7.
In the end, the PDF would be solved by adopting Pearson distribution family formulas of pattern IV and VI. By using Matlab function ([p, type, coefs] = pearspdf (X, mu, sigma, skew, kurt)), the probability density of Pearson distribution with mean 'mu', standard deviation 'sigma', skewness 'skew' and kurtosis 'kurt' could be easily worked out, which was evaluated at the value of X.
The qualification rates of four combinations were calculated and compared with the results of MC simulation, which have been shown in Table 8. It should be noted that in this case, the simulated number of MC was set to 10,000, which was beneficial to balance accuracy with efficiency. Meanwhile with regard to the qualification rate, in accordance with engineering requirements, FR should fall in the qualified range of 0.038 mm. See the table below for more details.  Fig. 6 Connection graph of four-stage assembly        Figure 7 has illustrated the variation distribution curve as well as MC statistical histograms according to Table 7. By comparing the results of the proposed method to Monte Carlo simulation, it was observed that the probability density curve agreed with the histogram statistics very well, as illustrated in Fig. 7a-d. Figure 7 also showed the qualification rate of the four combinations. In combination with Table 8, the relative error of qualification rate generated by the Pearson distribution family and MC method was less than 0.5 percent off. It turned out that the suggested method had a similar accuracy with MC method. In addition, the absolute time would be used as a measurement for comparison of the proposed method and MC simulation. The time function "tic & toc" was adopted in our MATLAB program to record the run time for the two methods. Table 9 has indicated the results. As seen from the Table 9, the computation of the proposed method is almost instantaneous. Using an Intel Core (TM) 2, 1.83 GHz, 1.5 GB RAM computer, the average time consumption of this model is 0.025 s. But when performing traditional Monte Carlo simulation (the simulated number was set to 10,000 in this case), the average time consumption reaches up to 0.408 s, which is sixteen times of the proposed method. The result demonstrates that the Pearson distribution family method could enhance work efficiency significantly.

Statistical Probability Analysis
Additionally, as can be seen from Fig. 7, whenever the FE's variation distribution changed, the FR's distribution pattern and its corresponding probability form would bring about apparent changes. Similar changes have taken place in the results of the assembly qualified rate as well. To be specific, in the four groups of case results, normal distribution (combination 1) guaranteed the highest pass percentage of 99.8%, but the left-skewed distribution (combination 3) almost caused the qualified rate to be zero due to overlarge mean shift. In Fig. 7c, the mean shift reached up to 0.1396 mm, which directly led to most of the products to be unqualified.
All the above results demonstrated that part's variation distribution type had an important effect on the final assembly precision. The conventional hypothesis of normal distribution would not be able to describe the real variation conditions. And the Pearson distribution family method is an effective solution for complicated probability density function. Therefore, to predict the accurate assembly precision, each part's real variation distribution must be determined.    This required to take the practically assembling technology and the processing level into account during the phase of product design. Only this could make results be more reasonable and reliable. It should be noted that in this paper the variation probability distribution of FE is regarded as a known condition. Historical manufacturing data of the parts should have been collected and analyzed firstly. So when we use a statistical variation analysis method to solve the distribution pattern of FR, the statistical variation results of FE can be directly invoked. However, MC simulation is an extremely timeconsuming process, especially when a large number of data samples are desired to improve the FR's prediction accuracy. Apparently, the statistical variation of FR can be hardly achieved in time. This is an urgent problem to be solved in this paper.

Experimental Study
The experiment was carried out to verify the accuracy and reliability of the results presented in the case study. A group of high-pressure turbine rotors (HPTRs) of an aeroengine have been coupled on the same rotary main axis of the air-bearing turntable. Their structures, dimensions,   4. An axial dial indicator was used to measure the run-out value of the end face, and a radial dial indicator was used to measure the run-out value of the final component flange. The former was adopted to control tilt amount, and the latter was applied to characterize the concentricity. As shown in Fig. 8, the HPTRs were assembled randomly by top stage control: after one stage rotor installation, its concentricity was measured. If the concentricity was unqualified, the rotor part would be dismantled and rotated to a new installation angle, ensuring that the threshold condition of concentricity can be met after each part stacking. Then the next stage rotor was assembled upward. Following this procedure, experiments have been implemented. Twenty times for repetitive assembly were completed. Measured data sets of the eccentricity of the top stage rotor were [0.026 0.025 0.061 0.048 0.075 0.061 0.037 0.077 0.023 0.057 0.020 0.043 0.047 0.042 0.028 0.069 0.054 0.066 0.019 0.026], the unit of which was millimetre, as shown in Table 10.
These statistics showed that the mean value of the repeated measures was 0.0452 mm, and the variance was 0.0186 mm 2 , which was closer to the simulating result of combination 2 referred to the Table 7. The measured data sets were further divided into several groups from 0.005 to 0.080 mm with interval step length 0.005 mm, and the variation distribution has been illustrated in Table 9. It was observed that the concentricity variation of the top stage rotor presented a typical right-skewed distribution, which conformed to the prediction results of combination 2 that was shown in Fig. 7b. These results proved that measured data in a certain error range accorded well with the calculations. The variation analysis method suggested in this paper was conducive to precision prediction and tolerance design for revolving components assembly.

Conclusions
This work was attempting to propose a novel variation analysis method for revolving components assembly. A calculating scheme was illustrated on a realistic case originating in aero-engine subassembly. Following conclusions could be drawn: (1) Unified Jacobian-Torsor theory was used to establish the serial kinematic chain, which combined the advantages of the torsor model which was suitable for variation representation and the Jacobian matrix which was suitable for variation propagation. For aero-engine rotors which contain lots of geometric variations and composite features, Jacobian-Torsor theory can effectively solve the 3D variation modelling. (2) Pearson distribution family method was used to perform statistical variation analysis. Compared to MC results with 10,000 simulation times consuming 0.408 s, this novel method has a similar accuracy, but much higher efficiency (the average time consumption is 0.025 s). (3) The result of FR depends on each FE's probability distribution pattern significantly. To obtain more accurate predictions, the real distribution states of parts variations should be detected and used, which is beneficial to tolerance design and performance forecast.