Statistical Variation Analysis Using Pearson Distribution Family Based on Jacobian-Torsor Model

Assembly variations are unavoidable due to parts’ geometrical errors. Statistical variation analysis is an effective method to quantitatively predict product quality in the original design stage. However, traditional methods can't handle the problem of abnormal distribution of the actual variation variables. Meanwhile, they are underdeveloped in regard to the complex geometrical errors in spatial 3D state. To overcome this problem, firstly, Jacobian-Torsor model is used to build the variation propagation, which is well suited to a complex assembly that contains large numbers of joints and geometric tolerances; secondly, Pearson distribution family is adopted to determine probability distribution pattern and build probability density function. By comparing results of the suggested method to the Monte Carlo method, it is observed that this novel method has the same accuracy, but much higher efficiency. The results also demonstrate that probability distribution types of the parts variations have a significant impact on the final assembling variation.


Introduction
Variation is inherent during any manufacturing and assembly process owing to the variability of the part, datum, positioner, fixture, etc.They cause minor variations in components from the nominal geometry [1].Statistical variation analysis is an effective way to predict product quality quantitatively in the design stage [2].Traditional methods are still based on engineering experience, extremum method [3,4], or the mean square root method [5,6].The results obtained by extremum method are overly pessimistic which will cause high fabricating cost, while the mean square root method only gives results for the mean-square variation which is primarily used to solve the plane size-chain problem.These methods are underdeveloped in regard to the geometrical errors in spatial 3D state.Meanwhile, they widely adopt normality assumption which is convenient for statistical evaluation.
As is known to all, most variables of parts variations disobey normal distribution in the actual machining process.For instance, if tool wearing plays a leading role in the machining process, the actual variation variables will satisfy uniform distribution; if one manufacturing process is accomplished using two machine tools, the mixing variables of parts variations will follow bimodal distribution; when turning shaft journals or apertures, operators usually keep a certain allowance, and this will lead to skewed distribution.Monte Carlo (MC) simulation is suitable for the data of abnormal distributions.Currently, MC method has been developed into a core of some commercial softwares such as CETOL, VSA, 3DCS, etc.A big disadvantage to this approach, however, is that large amounts of data samples are required to achieve accurate predictions [7], which inevitably result in time-consuming computation.Thus there is a clear need to adopt a simple and rapid probabilistic approach for statistical variation analysis.
Yang et al.
[8] studied the probability distribution of table-axis error in aero-engine assembly.Probability density functions (PDF) were derived based on Rayleigh distribution.Ghie [9] applied uncertain measurements into Jacobian-Torsor model.But the PDF deduced could not consider skewed distribution of variables.Ahsanullah and Hamedani [10] attempted to investigate various characterizations of certain univariate continuous distributions, providing a base of further research of multivariable continuous distributions.Kharoufehy and Chandraz [11] adopted a multivariate kernel density estimate to approximate the probability distribution.But it was difficult to determine the smoothing parameter (h).In this paper, we will introduce an innovation of Pearson distribution family [12] to model PDF that covers polymorphic probability distributions.
In addition, Jacobian-Torsor theory [13] will be used to model the variation propagation, which is well suited to a complex assembly that contains large numbers of joints and geometric tolerances.It combines the advantages of the torsor model which is suitable for tolerance representation and the Jacobian matrix which is suitable for tolerance propagation.

Jacobian-Torsor model
The unified Jacobian-Torsor model is composed of the Jacobian matrix and torsor model.Torsor model uses three translational vectors and three rotational vectors to express the position and orientation of an ideal surface.Supposing that S 0 is the nominal position and S 1 is an arbitrary variational surface relative to S 0 with tolerance t, the torsor model of S 1 which represents position and orientation can be represented as: where u, v, w are three specific vectors along the axes x, y and z in local reference system respectively; likewise, α, β and γ indicate the vectors around the axes x, y and z respectively.
These vectors are transferred to the objective controlled variable by virtue of the Jacobian matrix.Jacobian matrixes describe relative positions, as well as relative orientations between local reference frames and the global reference frame.The unified Jacobian-Torsor model can be expressed as: where α, β, γ, u, v and w are the same as Eq.(1); n is the number of functional elements (FEs);   ,   is tolerance interval where α should belong to; likewise, other torsors of components follow the similar way as α.
There are two types of FE pairs in assembly, i.e., internal FE (IFE) pairs which are composed of two FEs on the same part and contact FE (CFE) pairs which are made up of two FEs on different parts.Variations of each FE pair can be transferred to functional requirement (FR) with a Jacobian matrix, which of the i th FE can be expressed as: where [R 0 i ] represents the local orientation of the i th frame with respect to the 0 th frame which is global reference system; [R Pti ] is a projection matrix designating the unit vectors along local axes respectively for tolerance zone tilted according to the direction of tolerance analysis; [W i n ] is a skew-symmetric matrix [13] allowing the representation of the vector among the i th and n th frame (end point).

Pearson distribution family
Pearson system of distributions was firstly proposed by statistician K. Pearson in 1895, which included various types of distributions.Pearson system adopts four parameters to characterize each distribution.The probability density function (PDF) f(x) depends on the following differential equation: While a 1 =1, the parameters a 0 , b 0 , b 1 , and b 2 depend entirely on the first four order central moments of distributions μ 1 , μ 2 , μ 3 and μ 4 .It can also be expressed as moment ratio, i.e. β 1 =μ 3 2 /μ 2 3 and β 2 =μ 4 /μ 2 2 , as follows: where The general solution of f(x) are: To be specific, Pearson system mainly contains seven distribution patterns, which are elaborated in [14].These types can be distinguished according to the different distributions of roots of the quadratic equation , the Pearson distribution pattern can be judged by the value of D and λ, as provided in Table .1.

Statistical analysis
According to the Pearson system's characteristics, distribution pattern can be determined and PDF can be derived while the first four central moments are known.
As for arbitrary distributions, the first-order central moment μ 1 is identically equal zero; the second-order central moment μ 2 reflects data dispersion degree, which is equal to variance σ 2 ; the third-order central moment μ 3 is related to coefficient of skewness S, and S=μ 3 /σ 3 ; and the μ 4 is related to coefficient of kurtosis K, and K=μ 4 /σ 4 .Due to the average value, variance, skewness and kurtosis of a probability distribution can visually reflect the shape of the distribution, we will solve these four parameters in this paper to indirectly gain the first four central moments.
Eq. ( 2) can be represented as y FR =f(u i , vi, wi, αi, βi, γi,), (i=1,2,…,n), where y FR is the FR's variation, and (u i , vi, wi, αi, βi, γi,) is the FE's shifting variation or rotating variation in all directions.This representation is fit for deterministic variation analysis and beneficial to the error compensation.Supposing all of the FE's variations are randomly distributed and independent of one another, we can solve the four parameters of the FR by means of the FE's parameters.Then a 0 , b 0 , b 1 , b 2 can be obtained with the known distribution parameters of FR.Using the suggested method (Table.1),distribution pattern can be distinguished and PDF can be solved, after which the distribution curve can be drawn and assembly qualified rate will be calculated.The analysis process is illustrated in Fig. 1.

Case study
To validate the practicability of the solution mentioned above, a realistic subassembly originating in an aeroengine is adopted for statistical variation analysis.The aero-engine subassembly is composed of four precise revolutionary parts and all the coordinate frames are specified in Fig. 2 (a).Reference 0 is the global reference frame and the revolution axis is designated with a line passing through the base centre of the first component and is perpendicular to it.Each component's dimensions and composite tolerances which have been mentioned in Fig. 2 (a) are listed in Table .2.The pose position and orientation of the final rotor which is accumulated by all variations of four components and three joints are very critical to concentricity control in aero-engine assembly technical requirement.Based on the 3D variation analysis method proposed in Section 3, the solving procedures are elaborated below: Step 1: Build coordinate system and assembly graph.In this case, the global coordinate system (0) is set in the base centre of the rotor 1, and other local reference frames are built in the centre of related contact areas.The connection graph is shown in Fig. 2 (b).FR is the pose of the top rotor in relation to the 0th frame.As can be seen, there is a serial route of variation propagation for FR: IFE 1 -CFE 1 -IFE 2 -CFE 2 -IFE 3 -CFE 3 -IFE 4 .
Step 2: Establish the variation propagation model based on the above connection graph.Using the unified Jacobian-Torsor theory, the Jacobian matrixes and corresponding tolerancing torsors, as well as their variation ranges can be obtained.Details of FE 1 are presented in Table .3,and other FE S have similar representations.
Next, on the basis of Eq. ( 2), the final results can be obtained by multiplying the Jacobian matrix and torsor values.In this example, we only focus on the eccentricity in a certain direction, for instance the x axis direction, which can be extracted from u m and reformulated as Eq. ( 12).In addition, as listed in Table .3,all the variations of effective vectors in the torsor should be satisfied.
Step 3: Calculate distribution parameters of the FEs and then determine FR's distribution type and its PDF.
The distribution parameters of FR can be calculated as follows: where a i is the sensitivity coefficient, E xi is the ith FE's mathematical expectation, and E y is the FR's mathematical expectation.σ xi is the ith FE's standard deviation, and σ y is the FR's standard deviation.S i denotes the coefficient of skewness and K i denotes the coefficient of kurtosis.
To verify the influence of probability distribution types on the final assembly variations, we suppose the distribution types of FEs follow normal distribution, skewed distribution and uniform distribution respectively, and four probability combinations will be used as shown in Table .4,where skewed distribution is simulated by virtue of Beta approximation.The results of variation distribution parameters of the FR can be worked out according to Eq. ( 13) to (16), as listed in Table .5.Based on Eq. ( 5) to (10), we can further obtain the values of a 0 , b 0 , b 1 , b 2 , and deduce the variation distribution pattern by reference to Table .1,which has been shown in the last column of Table.5.Step 4: Solve PDF using Pearson distribution family formulas (pattern IV and VI).With the aid of Matlab function, we can easily get the probability denisty of the Pearsons distribution with mean 'mu', standard deviation 'sigma', skewness 'skew' and kurtosis 'kurt', evaluated at the values in X.
Calculate the qualification rate of four combinations, comparing with the MC results.In order to balance accuracy with efficiency, the simulation number of MC is set to 10000 in this case.With regard to the tolerance limits, FR should fall in the qualified range of 0.038mm in accordance with engineering requirements.Results are shown in Table .6,and Fig. 3 illustrates the distribution curve as well as MC statistical histograms.

Comparison and discussion
By comparing results of the suggested method to the MC ones, it is observed that the probability density curves drawn by Pearson distribution family method agree with the generated histogram statistics using MC simulation fairly well, which have been illustrated in Fig. 3 (a) to (d).
Meanwhile in Table .6,the error between qualification rates produced by the two ways is less than 0.5 percent off.This demonstrates that the proposed method has the same accuracy of MC method.Additionally, the former will calculate only once in the entire simulation process while the latter will calculate ten thousand times repeatedly.Pearson distribution family method has vastly enhanced the work efficiency.
In addition, as can be seen from the diagram, when the distribution pattern of FE's tolerance zone is changed, the shape of the variation distribution of the FR will bring about great changes, which appears to have happened at the assembly qualified rate as well.In the four groups of experiments, combination 1 (normal distribution) can guarantee the high pass percentage of 99.8% while combination 3 (left-skewed distribution) almost causes all products to be unqualified as a result of overlarge mean shift (=0.1396mm).All above results demonstrate that probability distribution types of the parts variations have a significant impact on the final assembling variation, and the traditional hypothesis of normal distribution cannot describe the true error conditions.Therefore, during the phase of product design, the practically processing technology and the level of processing must be taken into account to determine each component's real variation distribution.Only this can make results be more reasonable and reliable.

Conclusion
This work was attempting to propose a novel variation analysis method for complex assembly with multi joints.A calculating scheme was illustrated on a realistic case originating in aero-engine subassembly.Following conclusions could be drawn: first, for aero-engine rotors which contain lots of geometric tolerances, and planar and cylindrical features, Jacobian-Torsor theory can effectively solve the 3D variation modelling; second, compared to MC results with 10000 simulation times, Pearson distribution family method has the same accuracy, but much higher efficiency (runs only once); third, the result of FR depends on each FE's probability distribution pattern significantly.In order 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.

Table 2 .
Dimensions and geometric tolerances of each component(Unit: mm).

Table 4 .
Distribution type and parameters of the FEs.

Table 5 .
Distribution type and parameters of the FR.

Table 6 .
Qualification rate of four combinations.