An improved Jacobian-Torsor model for statistical variation solution in aero-engine rotors assembly

Abstract

Rotor assembly is one of the core components of aero-engine, which basically consists of multistage revolving components. With the influence of parts’ manufacturing errors and practical assembly technology, assembly variations are unavoidable which will cause insecurity and unreliable of the whole engine. Statistical variation solution is a feasible means to analyze assembly precision. When using the three-dimensional variation analysis in rotor assembly, two key issues cannot be well solved, which involve the variation expression (the over-positioning problem of multiple datums) and the variation propagation (revolving characteristic of the rotors). To overcome the deficiency, extended Jacobian matrix and updated torsor equation were derived and unified, which eventually resulted in the improved Jacobian-Torsor model. This model can both provide rotation regulating mechanism by introducing the revolution joint, and characterize the interaction between essential mating features. Multistage rotational optimization of four-stage aero-engine rotors assembly has been performed to demonstrate this solution in statistical way. Results showed that the proposed model was applicable and conducive to precision prediction and analysis in design preliminary stage.

Introduction

Rotor assembly is one of the core components of aero-engine, which basically consists of multistage revolving components. As shown in Figure 1, the rotor system has universal characteristics of complex internal structure and large number of joins, together with a demand of multiple coordination procedure. Due to the high precision requirement and installation complexity, the rotor assembling workload holds more than thirty percent of the engine’s total workforce.^1,2 For instance, the rotor system of GE90-115B aero-engine is made up of 10-stage high-pressure compressor rotors, six-stage low-pressure turbine rotors, and two-stage high-pressure turbine rotors. The tremendous dynamic loads such as vibration load, aerodynamic force, inertia force, torque, etc. will be put on rotors during the service process. If the centering scheme is not designed properly, or parts are unreasonably installed, or the balance test is inefficiency, the engine will appear violent vibration, directly affecting its reliability and security.

Figure 1.

Aero-engine rotors assembly.

To maintain the manufacturability and maintainability, aero-engine is produced as a complex assembly of individual revolving components. With the influence of parts’ manufacturing errors and practical assembly technology, assembly variations are unavoidable which will cause insecurity and unreliable of the whole engine. Generally, product quality depends mainly upon the assembly datum, sequence, parts’ type, etc. For aero-engine rotors, they have the rotation properties that every rotor part could be adjusted to a best mounting angle to realize assembly’s optimum concentricity performance. Statistical variation solution is a feasible means to analyze assembly precision. Both in product design stage and actual assembly stage, reasonably building a variation model and proposing a feasible assembly strategy is the effective way to cut down the costs and improve the products quality.

For rotor components assembly, conventional approach of variation analysis focuses on extremum solution and plane dimension chain modeling, which is unable to think about any of the geometrical deviations and the relationship among them; moreover, it overemphasizes the improvement of components’ manufacturing precision, ignoring the optimization method of assembling process. This makes predicted results is too conservative to control fabricating cost and processing difficulty. In fact, due to the inherent complexity of variation propagation, variation accumulation of rotors assembly could be almost not modeled and predicted. A mass of geometric deviations such as orientation, form and location information may be lost. To address this problem, in the past three decades, 3D variation analysis method^3,4 develops fast. At the initial stage, kinematic formulation ⁵ and the network model of zones and datums ⁶ were first suggested, then there is spatial dimensional chain model, ⁷ which mainly focuses on the expression of various variations. What followed is variational method^8,9 and the direct linearization method ¹⁰ ; thereafter the model of Jacobian transfer matrix, ¹¹ torsor,^12,13 Tolerance-Map ¹⁴ homogeneous transforming matrixes, ¹⁵ and the unified Jacobian-Torsor ¹⁶ were proposed successively. They provide the basic theories of variation propagation and accumulation based upon the previous study. Recently, the meta-modeling method, ¹⁷ polychromatic sets approach ¹⁸ and the model of shortest path graph ¹⁹ have reached a remarkable development.

The 3D variation model for rotor components assembling needs to solve two issues, involving the variation expression (over-positioning problem of multiple datums) and the variation propagation (revolving characteristic of the rotors). As shown in Figure 1, in aero-engine rotor assembly, the fixing structure commonly uses the transverse plane and the locating spigot round. This connection type consists of a pair of planar feature and cylindrical feather, which will influence and restrict each other. Meanwhile, rotors possess the rotation property that each revolving component can be adjusted to a suitable orientation to meet the desired concentricity performance, thereby causing combinatorial explosion. Some scholars have studied the variation expression of partial features, mainly including: (1) Boolean algebraic operations (BAO): Chen et al. ²⁰ introduced the boolean rules into the variation model, which is suited to analyze statistical solution, but lots of operations such as “union” and “intersection” bring about low computational efficiency. (2) Localization tolerancing with contact influence (CLIC): Chavanne and Anselmetti ²¹ simulated a group of contact surfaces with virtual extreme points to work out accumulated tolerance. Benichou and Anselmetti ²² further took thermal expansion effect into consideration to research temperature uncertainty. They presented this novel model to avoid partial features interference, but the problem of over-positioning could be not well solved. (3) Contact points-based solution (CPS): Jin et al. ²³ improved the small displacement torsor model by introducing contacting points into the conventional torsors. This approach is adopted to simulate rigid parts macroscopically, and is not proper for statistical solution. (4) Finite element method (FEM): Li et al. ²⁴ used FEM to solve the deformations of a novel flexible tooling for aero-engine rotors. Hong et al.^25,26 researched air band gap property by adopting the wave finite element model. Andersson ²⁷ explored the supersonic turbines efficiencies by developing a FE model of supersonic blades. Forslund et al. ²⁸ found the relation between thermal stresses and structural stresses by FEM and discussed the functional robustness. But FEM’s modeling is extremely complicated and only used to analyze variation as an assistant.

With regards to variation propagation, scholars have built the mapping relation between parts’ variations and assembly precision from the 3D dimensional chain perspective, mainly including: (1) Direct linearization method (DLM): DLM is fit for describing the variation accumulative process, and Wittwer et al. ²⁹ derived the linear dimensional chain function by expanding the first order taylor series. Yanlong et al. ³⁰ proposed a hybridization of vector loop and quasi-Monte Carlo method to establish the three-dimensional assembly chain. But this model neglected the geometric errors expression. (2) Connective assembly model (CAM): Hussain et al.^31–33 made a lot of researches on rotor assembly precision control. They adopted CAM model to control variations accumulation of cylinders. They presented different assembly strategies of parallel assembly and direct assembly, but numerous constraint equations limited the computational accuracy and efficiency. (3) Deterministic locating model (DeLM): Cai ³⁴ combined the linear variational approach to derive the deterministic locating function. It’s essentially an uncertain problem of variation propagation, which is affected by locating points and feature types. (4) Jacobian-Torsor model (J-T): Cao et al. ³⁵ synthesized the advantages about state space theory and torsor model to evaluate rotors assembly precision, and provided the idea of precision compensation. Ding et al. ³⁶ have established the variation propagation path along rotors’ angle and height by improving the Jacobian matrix. However, form and location errors of partial features cannot be fully considered. (5) Design of Experiment (DOE): Judt et al. ³⁷ investigated manual assembly performance of aircraft wing systems, under varying wing structure orientation. Chen et al. ³⁸ analyzed the variation’s abnormality of rotors assembly by using the Pearson theory and Taguchi experiments. Pierce and Rosen ³⁹ have brought machining mode variables into simulated experiments to improve accuracy. Apparently, the efficiency of experiment cannot meet the design requirements.

Table 1 has compared above variation models in detail. It is observed that J-T model is suitable for precision analysis of complex assembly, which includes plenty of geometric variations and connection joints. Rotor parts are commonly contain such assembling features. In consideration of aero-engine rotor’s configuration and assembly characteristic, the J-T theory will be adopted to construct the variation model of revolving components assembly in this article.

Table.1.

Variation models comparison.

Model type	Variation expression				Variation propagation
	BAO	CLIC	CPS	FEM	DLM	CAM	DeLM	J-T	DOE
Application object	Rigid	Rigid	Rigid	Rigid flexible	Rigid	Rigid	Rigid	Rigid	Rigid flexible
Object type	Face-based	Point-based	Point-based	Face-based	Point-based	Point-based	Point-based	Face-based	Face-based
Geometric variation	✓	✓	✓	✓	✓	✓	-	✓	✓
Variation coupling	-	-	✓	✓	-	✓	-	✓	✓
Solving method	Statistical	Deterministic	Deterministic	Deterministic	Deterministic	Deterministic Statistical	Deterministic	Deterministic Statistical	Deterministic
Solving efficiency	Low	Low	High	Low	High	Low	Low	High	Low
Engineering application	Vis-VSA	-	-	ABAQUS	CETOL	3D FDT	-	-	-
Difficulty level	Difficult	Easy	Easy	Difficult	Difficult	Easy	Difficult	Easy	Easy

However, this method cannot fully solve the issue of composite variations on multiple over-positioning datums. Because of the datum interference and multiple constraints, vector component of each torsor will not be independent any more. In original J-T model, it is common practice to integrate all limited values of tolerancing torsor into account, which is unreasonable and not suitable for statistical analysis. Another problem to be solved is that how to model the rotary characteristics in revolving axisymmetric components assembly. Jacobian matrix can be applied in serial mechanism as the linear arithmetic representation, and can easily map the velocity or displacement of internal joints to the terminal joint. However, it cannot perfectly provide the rotation regulating mechanism for the purpose of this paper.

The main idea for the article is to propose an available model for statistical variation solution in revolving components assembly. Extended Jacobian matrix and updated torsor equation will be deduced and unified to refine the traditional J-T theory, and multistage rotational optimization will be performed using this method in statistical way. By the solution, parts structure attribute and geometric attribute can be fully considered, and statistical concentricity and percentage contribution can be calculated.

Unified J-T model and its improved algorithm

Description of unified J-T 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:

$${\left[J\right]}_{\mathrm{FEi}}=\left[\begin{array}{ccc}{\left[R_0^i\right]}_{3\times 3}·{\left[R_{\mathrm{Pti}}\right]}_{3\times 3}\hfill & ⋮\hfill & {\left[W_i^n\right]}_{3\times 3}·({\left[R_0^i\right]}_{3\times 3}·{\left[R_{\mathrm{Pti}}\right]}_{3\times 3})\hfill \\ \multicolumn{1}{c}{\dots }& ⋮\hfill & \dots \hfill \\ \multicolumn{1}{c}{{\left[0\right]}_{3\times 3}}& ⋮\hfill & {\left[R_0^i\right]}_{3\times 3}·{\left[R_{\mathrm{Pti}}\right]}_{3\times 3}\hfill \end{array}\right]$$

(1)

where $\left[R_0^i\right]$ represents the i^th frame orientation; $\left[R_{\mathrm{Pti}}\right]$ is a projection matrix that allows efficient transformation of variation direction to the target analysis direction; $\left[W_i^n\right]$ describes the positional relation between adjacent frames in the form of skew-symmetric matrix. It is defined as equation (2):

$${\left[W_i^n\right]}_{3\times 3}=\left[\begin{array}{ccc}0\hfill & {\mathrm{dz}}_i^n\hfill & -{\mathrm{dy}}_i^n\hfill \\ \multicolumn{1}{c}{{-\mathrm{dz}}_i^n}& 0\hfill & {\mathrm{dx}}_i^n\hfill \\ \multicolumn{1}{c}{{\mathrm{dy}}_i^n}& {-\mathrm{dx}}_i^n\hfill & 0\hfill \end{array}\right]$$

(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.

Torsor model can describe the ideal position and orientation of any surfaces by three translational vectors and three rotational vectors. Assuming t is the tolerance between arbitrary feature S₁ and nominal feature S₀, there is:

$$T_{1-0}={\left[\begin{array}{cc}\begin{array}{ccc}\begin{array}{ccc}u\hfill & v\hfill & w\hfill \end{array}\hfill & \alpha \hfill & \beta \hfill \end{array}\hfill & \gamma \hfill \end{array}\right]}^T$$

(3)

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é ⁴⁰ 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 2 has shown seven primary torsor forms.

Table 2.

Features and their Torsor models.

No.	Element	Type of feature	Torsor model
1	Sphere		$\left\{\begin{array}{cc}-\hfill & u\hfill \\ \multicolumn{1}{c}{-}& v\hfill \\ \multicolumn{1}{c}{-}& w\hfill \end{array}\right\}$
2	Plane		$\left\{\begin{array}{cc}-\hfill & u\hfill \\ \multicolumn{1}{c}{\beta }& -\hfill \\ \multicolumn{1}{c}{\gamma }& -\hfill \end{array}\right\}$
3	Cylindrical surface		$\left\{\begin{array}{cc}-\hfill & -\hfill \\ \multicolumn{1}{c}{\beta }& v\hfill \\ \multicolumn{1}{c}{\gamma }& w\hfill \end{array}\right\}$
4	Prismatic surface		$\left\{\begin{array}{cc}\alpha \hfill & -\hfill \\ \multicolumn{1}{c}{\beta }& v\hfill \\ \multicolumn{1}{c}{\gamma }& w\hfill \end{array}\right\}$
5	Surface of revolution		$\left\{\begin{array}{cc}-\hfill & u\hfill \\ \multicolumn{1}{c}{\beta }& v\hfill \\ \multicolumn{1}{c}{\gamma }& w\hfill \end{array}\right\}$
6	Helicoid		$\begin{array}{l}\left\{\begin{array}{cc}-& u\\ \beta & v\\ \gamma & w\end{array}\right\}\\ \iff u\ne \frac{pa}{2\pi }\end{array}$
7	Arbitrary shape curved surface		$\left\{\begin{array}{cc}\alpha \hfill & u\hfill \\ \multicolumn{1}{c}{\beta }& v\hfill \\ \multicolumn{1}{c}{\gamma }& w\hfill \end{array}\right\}$

where “–” means the vector is invalid in the corresponding direction. a and p respectively denote thread tooth link angle and screw pitch.

Take No. 3 cylindrical surface as an example. In terms of degree of freedom, it has four effective DOFs in its tolerance zone. According to Salomons’s research, tolerances are only meaningful in directions of all non-invariant degrees. ⁴¹ That is, in a torsor model, the number of effective vectors equals to the non-invariant degrees of a feature. More specifically, cylindrical feature originally has six DOFs, but there are only four effective DOFs within its tolerance zone. As shown in Table 2, that is two rotations (α, β) and two translations (v, w). In the tolerance analysis, only these four variations affect the analysis results.

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.^42,43 The form of J-T model could be written as:

$$\begin{array}{c}{\left[\begin{array}{c}(\underset{¯}{u},\overline{u})\hfill \\ \multicolumn{1}{c}{(\underset{¯}{v},\overline{v})}\\ \multicolumn{1}{c}{(\underset{¯}{w},\overline{w})}\\ \multicolumn{1}{c}{(\underset{¯}{\alpha },\overline{\alpha })}\\ \multicolumn{1}{c}{(\underset{¯}{\beta },\overline{\beta })}\\ \multicolumn{1}{c}{(\underset{¯}{\gamma },\overline{\gamma })}\end{array}\right]}_{\mathrm{FR}}=\left[\begin{array}{ccc}{\left[J\right]}_{\mathrm{FE}1}\hfill & \cdots \hfill & {\left[J\right]}_{\mathrm{FEn}}\hfill \end{array}\right]\hfill \\ \times {\left[\begin{array}{ccc}{\left[\begin{array}{c}(\underset{¯}{u},\overline{u})\hfill \\ \multicolumn{1}{c}{(\underset{¯}{v},\overline{v})}\\ \multicolumn{1}{c}{(\underset{¯}{w},\overline{w})}\\ \multicolumn{1}{c}{(\underset{¯}{\alpha },\overline{\alpha })}\\ \multicolumn{1}{c}{(\underset{¯}{\beta },\overline{\beta })}\\ \multicolumn{1}{c}{(\underset{¯}{\gamma },\overline{\gamma })}\end{array}\right]}_{\mathrm{FE}1}\hfill & \cdots \hfill & {\left[\begin{array}{c}(\underset{¯}{u},\overline{u})\hfill \\ \multicolumn{1}{c}{(\underset{¯}{v},\overline{v})}\\ \multicolumn{1}{c}{(\underset{¯}{w},\overline{w})}\\ \multicolumn{1}{c}{(\underset{¯}{\alpha },\overline{\alpha })}\\ \multicolumn{1}{c}{(\underset{¯}{\beta },\overline{\beta })}\\ \multicolumn{1}{c}{(\underset{¯}{\gamma },\overline{\gamma })}\end{array}\right]}_{\mathrm{FEn}}\hfill \end{array}\right]}^T\hfill \end{array}$$

(4)

where $(\underset{¯}{\alpha },\overline{\alpha })$ is the variation zone, and β, γ, u, v, and w follow the similar way as α. Other parameters have the same meanings as those in equation (1).

Improved algorithm

Extended Jacobian matrix

The classical J-T model described above could not fully apply to the rotating optimization-oriented operation. Concerning this issue, a revolution joint (RJ) is introduced to represent the revolving characteristics of the rotors, by which the extended Jacobian matrix is derived. To be specific, for revolving components assembly, a planar joint combined with another cylindrical joint is most frequently used in the reliable positioning structure, especially in aero-engine rotors assembly. As illustrated in Figure 2, this assembly is composed of two joints which are the planar joint formed by flange plane and shaft end face, and the cylindrical joint formed by locating spigot round. It is actually a revolution joint (RJ) that introduced in this paper. As a result, the RJ simultaneously involves one large planar joint and one short cylindrical joint.

Figure 2.

Revolution joint.

In terms of degree of freedom, it is common knowledge that an unconstrained object has six degrees of freedom (DOF), that is, three rotational DOFs and three translational DOFs. Constraints caused by locating or assembling will eliminate some or all component’s DOFs. In the RJ, the planar joint restricts two rotational DOFs and one translational DOF, and the cylindrical joint restricts the remaining two translational DOFs, finally retaining a rotational DOF around rotor’s revolving axis. So in general, the RJ will restrict five DOFs of the rotor part. It should be noted that the DOFs have different meanings in the new RJ and the Cylindrical Feature. The former concerns the mechanism motion constraint in the conventional sense, and the latter concerns the effective vector representation in tolerance semantics and torsor expression.

For n-stage revolving components assembly, the FR results will change while each rotor rotates around its turning spindle axis, and the assembling precision will enhance by selecting an appropriate orientation at each stage. If component i rotates with an angle θ_i around axes z, the local orientation matrix $\left[R_0^i\right]$ in global reference system could be expressed as equation (5), which is the concrete form of RJ:

$$\mathrm{Revolution}\mathrm{joint}:{\left[R_0^i\right]}_{3\times 3}=\left[\begin{array}{ccc}cos{\theta }_i& -sin{\theta }_i& 0\\ sin{\theta }_i& cos{\theta }_i& 0\\ 0& 0& 1\end{array}\right]$$

(5)

Figure 3(a) and (b) respectively show an n-stage assembly and its connection graph. The global coordinate system 0 and local system i (i = 1, 2, ⋯, n) have been established in the middle of related characteristics. The upper end faces of every part uniformly have a position tolerance T_PO (=t_i1) and an orientation tolerance T_PRA (=t_i2). The FR is the accumulated tolerance on the nth rotor’s upper end face along the eccentric direction, measuring in the global reference frame. There are totally n IFEs and (n−1) CFEs which contain (n−1) RJs.

Figure 3.

The n-stage assembly: (a) n rotors assembly, and (b) connection graph.

Substituting equation (5) into equation (1), and incorporating the tolerances and dimensions shown in Figure 3(a), the extended Jacobian matrix can be expressed as equations (6) and (7):

$$\mathrm{Internal}\mathrm{pairs}:{\left[J_e\right]}_{\mathrm{IFEn}}=\left[\begin{array}{cccccc}\cos {\theta }_{(n-1)}\hfill & -\sin {\theta }_{(n-1)}\hfill & 0\hfill & (\sum _{i=1}^n{H}_i)\times \sin {\theta }_{(n-1)}\hfill & (\sum _{i=1}^n{H}_i)\times \cos {\theta }_{(n-1)}\hfill & 0\hfill \\ \multicolumn{1}{c}{\sin {\theta }_{(n-1)}}& \cos {\theta }_{(n-1)}\hfill & 0\hfill & -(\sum _{i=1}^n{H}_i)\times \cos {\theta }_{(n-1)}\hfill & (\sum _{i=1}^n{H}_i)\times \sin {\theta }_{(n-1)}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \cos {\theta }_{(n-1)}\hfill & -\sin {\theta }_{(n-1)}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \sin {\theta }_{(n-1)}\hfill & \cos {\theta }_{(n-1)}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]$$

(6)

$$Contactpairs:{\left[J_e\right]}_{CFEn}=\left[\begin{array}{cccccc}cos{\theta }_{(n)}& -sin{\theta }_{(n)}& 0& (\sum _{i=1}^n{H}_i)\times sin{\theta }_{(n)}& (\sum _{i=1}^n{H}_i)\times cos{\theta }_{(n)}& 0\\ sin{\theta }_{(n)}& cos{\theta }_{(n)}& 0& -(\sum _{i=1}^n{H}_i)\times cos{\theta }_{(n)}& (\sum _{i=1}^n{H}_i)\times sin{\theta }_{(n)}& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& cos{\theta }_{(n)}& -sin{\theta }_{(n)}& 0\\ 0& 0& 0& sin{\theta }_{(n)}& cos{\theta }_{(n)}& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]$$

(7)

where θ_i denotes the mounting angle of the (i+1)^th part with respect to the global coordinates, and H_i denotes the height of the i^th part.

The extended Jacobian matrix indicates the vector relation between the torsors of all the FEs and the FR, combining with the multistage mounting angles θ_i. By virtue of θ_i in RJ, multistage rotating optimization can be performed and the overall quality of assembly concentricity can be guaranteed.

Torsor equation

Revolving body adopted in this article mainly contains two categories of essential features, respectively are planar feature and cylindrical feature. They cooperate to affect the overall concentricity of assembly.

Planar feature between the adjacent matching surfaces constitutes the planar joints, which are selected as a major contact pair. The primary concern in assembling is their contact reliability. Figure 4 shows an individual cylinder. The planar feature contains three invariant degrees, indicating that in its torsor model, three effective vectors still exist in total.

Figure 4.

Torsor model of planar feature.

Roy et al. ⁴⁴ have given the preliminary mathematical derivation of the constraints and variations of planar feature under the independent tolerance constraints. Any variant plane equation can be written as:

$$x=u+\beta ·z+\gamma ·y$$

(8)

First considering position tolerance T_PO, the constraint equation of planar feature under position tolerance can be expressed as equation (9):

$$\{\begin{array}{c}-(T_{\mathrm{SU}}+T_{\mathrm{SL}})/2\le (\beta z+\gamma y)\le (T_{\mathrm{SU}}+T_{\mathrm{SL}})/2\hfill \\ \multicolumn{1}{c}{-T_{\mathrm{SL}}\le (u+\beta z+\gamma y)\le T_{\mathrm{SU}}}\end{array}$$

(9)

And the variations of the three vectors (u, β, γ) are:

$$\{\begin{array}{c}-T_{\mathrm{SL}}\le u\le T_{\mathrm{SU}}\hfill \\ \multicolumn{1}{c}{-(T_{\mathrm{SU}}+T_{\mathrm{SL}})/d\le \beta \le (T_{\mathrm{SU}}+T_{\mathrm{SL}})/d}\\ \multicolumn{1}{c}{-(T_{\mathrm{SU}}+T_{\mathrm{SL}})/d\le \gamma \le (T_{\mathrm{SU}}+T_{\mathrm{SL}})/d}\end{array}$$

(10)

where T_SL and T_SU denote the limit values of positional tolerances respectively, and their sum is equal to T_PO; △x(y, z) = x(y_i, z_i)−x(y_j, z_j), and (y_i, z_i), (y_j, z_j) ϵ{(−d/2, −d/2), (−d/2, d/2), (d/2, −d/2), (d/2, d/2)}; d denotes cylinder diameter.

Besides T_PO, T_PRA would also strictly control the feature variation. Composite tolerances must obey the envelop principle (Rule #1, ASME-2009). More specifically, a feature’s rotational displacements resulting from orientation tolerance are essentially subject to the position tolerance. Therefore, the value of orientation tolerance should be smaller than that of position tolerance according to the tolerance standard. In particular, rotational displacements have to shrink to zero when translational displacements reach the maximum value. Figure 5 has shown the detailed process.

Figure 5.

The relationship between T_PO and T_PRA for the planar feature.

Further taking into account T_PRA, the torsor’s constraint equations for any rotor’s upper end face could be finally rewritten as equation (11):

$$\{\begin{array}{c}-T_{\mathrm{PRA}}/2\le (\beta z+\gamma y)\le T_{\mathrm{PRA}}/2\hfill \\ \multicolumn{1}{c}{-T_{\mathrm{SL}}\le (u+\beta z+\gamma y)\le T_{\mathrm{SU}}}\end{array}$$

(11)

And now, the variations of the three vectors (u, β, γ) are:

$$\{\begin{array}{c}-T_{\mathrm{SL}}\le u\le T_{\mathrm{SU}}\hfill \\ \multicolumn{1}{c}{-T_{\mathrm{PRA}}/d\le \beta \le T_{\mathrm{PRA}}/d}\\ \multicolumn{1}{c}{-T_{\mathrm{PRA}}/d\le \gamma \le T_{\mathrm{PRA}}/d}\end{array}$$

(12)

As can be seen from equations (11) and (12), the position tolerance mainly limits the translational displacements while the rotations are limited by orientation tolerances. The six vector components of the torsor significantly depend on each other. Furthermore, expressions of other orientation tolerances such as perpendicularity and angularity can be deduced similar to the parallelism.

Cylindrical feature of the revolving component is also constrained by composite tolerance, to ensure that the concentricities of different sections are maintained in a high precision range during the machining process. Figure 6(a) shows a part of an aero-engine rotor. One of the sections is marked with position tolerance φT_PO and orientation tolerance φT_PRA. According to the envelop principle mentioned above, orientation tolerance is restricted by position tolerance, that is φT_PRA≤φT_PO. This means when the cylindrical feature reaches the maximum material size, its orientation tolerance φT_PRA must be reduced to zero.

Figure 6.

Torsor model of cylindrical feature: (a) rotor part, (b) cylindrical domain of φT_PO, (c) plane projection of x₁o₁y₁.

As shown in Figure 6(a), the cylindrical section H₂ is restricted by position tolerance φT_PO and parallelism tolerance φT_PRA, and the cylindrical section H₁ is set as datum. According to tolerance semantics, position φT_PO mainly restricts the position variation of feature H₂’s center line, and parallelism φT_PRA further restricts the direction of H₂’s center line relative to the datum axis of H₁. Therefore, first considering φT_PO, the center line with a different position constitutes an envelope cylinder.

Figure 6(b) shows this magnified cylindrical domain, in which the global reference frame O and local frame O₁ are established. Cylindrical region φT_PO is the variable range of feature H₂’s center line. In the φT_PO, P is the random point vector, C _d is the direction vector, and C _P is the position vector. According to the definition of cross product, C _d ×(P−C _P ) is a radial vector situated in the feasible zone φT_PO. Thus the vector equation of the center line of cylindrical feature can be determined as:

$$|{\mathbf{C}}_{\mathit{d}}\times (\mathbf{P}-{\mathbf{C}}_{\mathit{P}})|\le \frac{T_{\mathrm{PO}}}{2}$$

(13)

In the local coordinate system O₁, C _d =[1 0 0], C _P =[0 0 0]. Substitute them in the above equation and equation (14) can be obtained:

$$\begin{array}{c}|{\mathbf{C}}_{\mathit{d}}\times (\mathbf{P}-{\mathbf{C}}_{\mathit{P}})|\hfill \\ =|\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\times (\left[\begin{array}{ccc}x\hfill & y\hfill & z\hfill \end{array}\right]-\left[\begin{array}{ccc}0\hfill & 0\hfill & 0\hfill \end{array}\right])|\hfill \\ =|\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\times \left[\begin{array}{ccc}x\hfill & y\hfill & z\hfill \end{array}\right]|\le \frac{T_{\mathrm{PO}}}{2}\hfill \end{array}$$

(14)

φT_PO contains two invariants: translation u along the axis x and rotation α around the axis x, and Figure 6(c) has illustrated the main plane projection of x₁o₁y₁. As shown, the actual axes position of o₂P experienced two minor changes: rotate γ around o₁z₁, and translate v along o₁y₁. Now the P’s y-coordinate becomes (v+x·tanγ). In the same way, the P’s new z-coordinate becomes (w+x·tanβ).

Combined with small displacement approximation theory, tanβ≈β, tanγ≈γ, the Point P’s new coordinate becomes [x (v+x·tanγ) (w+x·tanβ)]. Substituting it in the equation (14), the constraint equation of cylindrical feature under the influence of φT_PO can be derived as:

$${(v+\gamma x)}^2+{(w+\beta x)}^2\le {\left(\frac{T_{\mathrm{PO}}}{2}\right)}^2$$

(15)

where x∈H₂.

And the variations of the effective vectors (v, w, β, γ) are:

$$\{\begin{array}{c}\multicolumn{1}{c}{-\frac{T_{\mathrm{PO}}}{2}\le v\le \frac{T_{\mathrm{PO}}}{2}}\\ \multicolumn{1}{c}{-\frac{T_{\mathrm{PO}}}{2}\le w\le \frac{T_{\mathrm{PO}}}{2}}\\ \multicolumn{1}{c}{-\frac{T_{\mathrm{PO}}}{2H_2}\le \beta \le \frac{T_{\mathrm{PO}}}{2H_2}}\\ \multicolumn{1}{c}{-\frac{T_{\mathrm{PO}}}{2H_2}\le \gamma \le \frac{T_{\mathrm{PO}}}{2H_2}}\end{array}$$

(16)

where H₂ is the length of the cylinder feature.

The cylindrical feature is also influenced by the orientation tolerance φT_PRA besides the position tolerance φT_PO. So the deflection of H₂’s center line is further limited by φT_PRA, which results in additional constraints on the two rotating vectors β and γ. And the new variation equations are:

$$\{\begin{array}{c}-\frac{T_{\mathrm{PO}}}{2}\le v\le \frac{T_{\mathrm{PO}}}{2}\hfill \\ -\frac{T_{\mathrm{PO}}}{2}\le w\le \frac{T_{\mathrm{PO}}}{2}\hfill \\ -\frac{T_{\mathrm{PRA}}}{2H_2}\le \beta \le \frac{T_{\mathrm{PRA}}}{2H_2}\hfill \\ -\frac{T_{\mathrm{PRA}}}{2H_2}\le \gamma \le \frac{T_{\mathrm{PRA}}}{2H_2}\hfill \end{array}$$

(17)

Since the φT_PRA does not affect the translation vectors, the variation and constraint equations of v and w will remain unchanged. Equations (15) and (17) are the final torsor equations for the cylindrical feature. The above torsor equations of revolving components ensure that every torsor conforms to the ASME standard and is available to participate in computation statistically. It is beneficial to tolerance specification and allocation of the rotationally symmetric components assembly, such as aero-engine rotors design.

In general, the proposed rule to control the translation degrees and the rotation degrees by a composite tolerance composed of a position tolerance and an orientation tolerance is valid in general, and common in engineering. This is because for precision assemblies, a position tolerance is not sufficient to control both the position and direction of the feature. For example, the feature may deflect the most with the value of T_PO when only one position tolerance acts, which result in improper assembly. Therefore, in order to reduce the deflection and meet the special requirements of translation degrees to the greatest extent, an additional orientation tolerance needs to be further applied.

Position tolerance and orientation tolerance include many types, as shown in Table 3. The tolerance combination of position (No. 5, 6, 10, 11, 12) and orientation (No. 7, 8, 9) is a conventional mode. Either way, the mathematical expression is the same, meeting the envelop principle and constraint rule, differentiated only by their tolerance size. To keep generality, we will further summarize the uniform variation equations and constraint equations. Details will be shown in Table 4 for the general revolving components.

Table.3.

Geometric tolerance.

Tolerance	No.	Item	Symbol	Tolerance		No.	Item	Symbol
Form	1	Straightness		Location	Directed	7	Parallelism
	2	Flatness				8	Perpendicularity
	3	Roundness				9	Angularity
	4	Cylindricity			Positional	10	Coaxiality/Concentricity
Form or location	5	Line profile				11	Symmetry
						12	Position
	6	Surface profile			Run-out	13	Circle run-out
						14	Total run-out

Table.4.

Variations and constraints of effective vectors.

FE	Variations	Constraints
IFEi	$\begin{array}{c}-t_{i1}/2\le w_{i\mathrm{p}}\le t_{i1}/2\hfill \\ -t_{i2}/D_i\le {\alpha }_{i\mathrm{p}}\le t_{i2}/D_i\hfill \\ -t_{i2}/D_i\le {\beta }_{i\mathrm{p}}\le t_{i2}/D_i\hfill \\ -t_{i3}/2\le u_{\mathrm{ic}}\le t_{i3}/2\hfill \\ -t_{i3}/2\le v_{\mathrm{ic}}\le t_{i3}/2\hfill \\ -t_{i4}/2H_{\mathrm{i}}\le {\alpha }_{\mathrm{ic}}\le t_{i4}/2H_{\mathrm{i}}\hfill \\ -t_{i4}/2H_{\mathrm{i}}\le {\beta }_{\mathrm{ic}}\le t_{i4}/2H_{\mathrm{i}}\hfill \end{array}$	$\begin{array}{c}-t_{i1}/2\le w_{i\mathrm{p}}+{\beta }_{i\mathrm{p}}·x_{i\mathrm{p}}\hfill \\ +{\alpha }_{i\mathrm{p}}·y_{i\mathrm{p}}\le t_{i1}/2\hfill \\ -t_{i2}/2\le {\beta }_{i\mathrm{p}}·x_{i\mathrm{p}}\hfill \\ +{\alpha }_{i\mathrm{p}}·y_{i\mathrm{p}}\le t_{i2}/2\hfill \\ {(v_{\mathrm{ic}}+{\alpha }_{\mathrm{ic}}·z_{\mathrm{ic}})}^2\hfill \\ +{(u_{\mathrm{ic}}+{\beta }_{\mathrm{ic}}·z_{\mathrm{ic}})}^2\le {(t_{i3}/2)}^2\hfill \end{array}$

Solution

Substituting the extended Jacobian matrixes equations (6) and (7) into equation (4), and combining with the updated torsor equations (11), (12), (15), and (17), the final concentricity variation result of the n-stage rotors assembly is obtained as equation (18):

$${\left[\mathbf{FR}\right]}_n={\left[\begin{array}{c}\delta u\hfill \\ \multicolumn{1}{c}{\delta v}\\ \multicolumn{1}{c}{\delta w}\\ \multicolumn{1}{c}{\delta \alpha }\\ \multicolumn{1}{c}{\delta \beta }\\ \multicolumn{1}{c}{\delta \gamma }\end{array}\right]}_{F{R}_n}=\left[\begin{array}{c}\sum _{i=2}^n\{[{\beta }_i·\cos {\theta }_{(i-1)}+{\alpha }_i·\sin {\theta }_{(i-1)}]·\left(\sum _{j=1}^i{H}_j\right)\}+{\beta }_1·H_1\hfill \\ \multicolumn{1}{c}{\sum _{i=2}^n\{[{\alpha }_i·\sin {\theta }_{(i-1)}-{\beta }_i·\cos {\theta }_{(i-1)}]·\left(\sum _{j=1}^i{H}_j\right)\}-{\beta }_1·H_1}\\ \multicolumn{1}{c}{\sum _{i=2}^n{w}_i}\\ \multicolumn{1}{c}{\sum _{i=2}^n[{\alpha }_i·\cos {\theta }_{(i-1)}-{\beta }_i·\sin {\theta }_{(i-1)}]+{\alpha }_1}\\ \multicolumn{1}{c}{\sum _{i=2}^n[{\beta }_i·\cos {\theta }_{(i-1)}+{\alpha }_i·\sin {\theta }_{(i-1)}]+{\beta }_1}\\ \multicolumn{1}{c}{0}\end{array}\right]$$

(18)

where θ_i denotes the mounting angle of the (i+1)^th part with respect to the global coordinates; H_i denotes the height of the i^th part; α_i and β_i are rotating vectors around x and y axes respectively; w_i is the displacement vector on z axes in the local frame.

Taking advantage of δu and δv, which are the vectors along x and y axes, the nth rotor’s concentricity variation could be derived. However, operators pay more attention on the concentricity of the overall assembly rather than that of the top stage rotor only. As illustrated in Figure 7, even if the concentricity of the last stage rotor has been controlled well, other rotors may be also dimensionally out of tolerance.

Figure 7.

Concentricity control of multistage rotors assembly.

Thus the target eccentricity of the overall assembly can be expressed as equation (19):

$$\{\begin{array}{c}{\epsilon }_{\mathrm{multistage}}=min\left\{\sqrt{{({\xi }_2·{\epsilon }_2)}^2+{({\xi }_3·{\epsilon }_3)}^2+\cdots +{({\xi }_n·{\epsilon }_n)}^2}\right\}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{t}}=\sqrt{\delta {u_t}^2+\delta {v_t}^2}(t=2,3,\cdots ,n)}\\ \multicolumn{1}{c}{\left|\delta u_t\right|\le |\sum _{i=2}^t\{[{\beta }_i·\cos {\theta }_{(i-1)}+{\alpha }_i·\sin {\theta }_{(i-1)}]·\left(\sum _{j=1}^i{H}_j\right)\}+{\beta }_1·H_1|}\\ \multicolumn{1}{c}{\left|\delta v_t\right|\le |\sum _{i=2}^t\{[{\alpha }_i·\sin {\theta }_{(i-1)}-{\beta }_i·\cos {\theta }_{(i-1)}]·\left(\sum _{j=1}^i{H}_j\right)\}-{\beta }_1·H_1|}\end{array}$$

(19)

where ξ_n is the weighting coefficient of concentricity deviation, denoting the nth component’s importance. In this article, all parts are considered equally important, so ξ_n = 1.

Meanwhile, the variations and constraints of effective vectors for cylindrical torsors should satisfy the requirements of Table 4.

Up to now, the 3D tolerance analysis model for revolving axisymmetric components has been established. Objective function and governing equations for multistage concentricity control have been derived. In order to solve the model, Monte Carlo statistic method is adopted and genetic algorithm ⁴⁵ is used for seeking optimization. Figure 8 shows the statistical analysis flow chart on the basis of improved J-T model. Details on the analysis procedures of multistage rotors assembly will be elaborated in “Case study.”

Figure 8.

Statistical variation analysis procedure using improved J-T model.

Case study

An actual subassembly originating in an aero-engine is used to demonstrate the utility of the improved J-T model. The rotor assembly consists of four revolving components, and all the coordinate systems have been shown in Figure 9. Frame 0 is the global coordinate system and the revolving axis is identified with a line passing through the rotor center. Each component contains a profile tolerance and a parallelism tolerance (located at the upper end of the part, as shown in Figure 3(a)), as well as an axiality tolerance. The specific dimensions and tolerances are listed in Table 5.

• Step 1: Establish reference coordinate and connection graph. The global reference frame 0 is built in the first rotor’s base center, and other local coordinate systems are set in the related contacting characteristics’ centre. In this case, Figure 10 shows the assembling graph. It is observed that a serial route for FR of variation propagation exists: IFE₁-CFE₁-IFE₂-CFE₂-IFE₃-CFE₃-IFE₄-FR.

• Step 2: Construct the variation propagation model utilizing the assembling graph. Adopting the improved J-T model, extended Jacobian matrixes and relevant deviation torsors can be obtained. The variations and constraints of FE₁ are presented in Table.6. Other FE_S have similar representations.

Figure 9.

Parts drawing of four-stage rotors assembly.

Table 5.

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

Component	Dbase	Dtop	H i	tprofile	tparallelism	taxiality
1	107	405	315	0.020	0.015	0.012
2	405	410	135	0.020	0.015	0.012
3	410	400	165	0.020	0.015	0.012
4	400	440	365	0.020	0.015	0.012

Figure 10.

Connection graph of four-stage rotors assembly.

Table 6.

Jacobian matrixes and deviation Torsor of FE₁.

No.	FE s	Jacobian matrixes	Deviation torsors
1	IFE 1	$\left[\begin{array}{cccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 315\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & -315\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]$	$\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{w_1}\\ \multicolumn{1}{c}{{\alpha }_1}\\ \multicolumn{1}{c}{{\beta }_1}\\ \multicolumn{1}{c}{0}\end{array}\right]$
2	CFE 1	$\left[\begin{array}{cccccc}\cos {\theta }_1\hfill & -\sin {\theta }_1\hfill & 0\hfill & 315\sin {\theta }_1\hfill & 315\cos {\theta }_1\hfill & 0\hfill \\ \multicolumn{1}{c}{\sin {\theta }_1}& \cos {\theta }_1\hfill & 0\hfill & -315\cos {\theta }_1\hfill & 315\sin {\theta }_1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \cos {\theta }_1\hfill & -\sin {\theta }_1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \sin {\theta }_1\hfill & \cos {\theta }_1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]$	$\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right]$

Next, the final multistage assembly formula could be derived based on equation (19), as shown in equation (20).

$$\{\begin{array}{c}{\epsilon }_{\mathrm{multistage}}=min\left\{\sqrt{{{\epsilon }_2}^2+{{\epsilon }_3}^2+{{\epsilon }_4}^2}\right\}\hfill \\ \multicolumn{1}{c}{=min\left\{\sqrt{\delta {u_2}^2+\delta {v_2}^2+\delta {u_3}^2+\delta {v_3}^2+\delta {u_4}^2+\delta {v_4}^2}\right\}}\\ \multicolumn{1}{c}{\left|\delta u_2\right|\le |[{\beta }_2·\cos {\theta }_1+{\alpha }_2·\sin {\theta }_1]·\left(\sum _{j=1}^2{H}_j\right)+{\beta }_1·H_1|}\\ \multicolumn{1}{c}{\left|\delta v_2\right|\le |[{\alpha }_2·\sin {\theta }_1-{\beta }_2·\cos {\theta }_1]·\left(\sum _{j=1}^2{H}_j\right)-{\beta }_1·H_1|}\\ \multicolumn{1}{c}{\left|\delta u_3\right|\le |\sum _{i=2}^3\{[{\beta }_i·\cos {\theta }_{(i-1)}+{\alpha }_i·\sin {\theta }_{(i-1)}]·\left(\sum _{j=1}^i{H}_j\right)\}+{\beta }_1·H_1|}\\ \multicolumn{1}{c}{\left|\delta v_3\right|\le |\sum _{i=2}^3\{[{\alpha }_i·\sin {\theta }_{(i-1)}-{\beta }_i·\cos {\theta }_{(i-1)}]·\left(\sum _{j=1}^i{H}_j\right)\}-{\beta }_1·H_1|}\\ \multicolumn{1}{c}{\left|\delta u_4\right|\le |\sum _{i=2}^4\{[{\beta }_i·\cos {\theta }_{(i-1)}+{\alpha }_i·\sin {\theta }_{(i-1)}]·\left(\sum _{j=1}^i{H}_j\right)\}+{\beta }_1·H_1|}\\ \multicolumn{1}{c}{\left|\delta v_4\right|\le |\sum _{i=2}^4\{[{\alpha }_i·\sin {\theta }_{(i-1)}-{\beta }_i·\cos {\theta }_{(i-1)}]·\left(\sum _{j=1}^i{H}_j\right)\}-{\beta }_1·H_1|}\end{array}$$

(20)

• Step 3: Construct variations and constraints of the composite tolerances based on Table 4. Table 7 listed the final forms of the variations and constraints of FE₁. Other FE_S have similar representations.

• Step 4: Create random vector values of every FEs. The variations shown in Table 7 are adopted as random number interval. Using Monte Carlo simulation, each deviation torsor needs to be provided a desired statistical distribution. In this case, all the deviation torsors are assumed to obey normal distribution. Every vector’s mean value and standard deviation can be obtained in accordance with equations (21) and (22).

$${\mu }_T=\frac{V_{\mathrm{SL}}+V_{\mathrm{SU}}}{2}$$

(21)

$${\sigma }_T=\frac{V_{\mathrm{SU}}-V_{\mathrm{SL}}}{Z}$$

(22)

where V_SL and V_SU are the limit values of deviation torsor. Z denotes the standardized normal value. Here the value of Z is set to 6, indicating that the rejection rate is 0.27%.

• Step 5: Filter the generated values of every FEs according to the constraint equations. This operation ensures that all the random values involved in the calculation are absolutely conforming to the tolerancing standards.

• Step 6: Use the vectors of every torsor picked from step 5 to calculate the target eccentricity ${\epsilon }_{\mathrm{multistage}}$ . Take acceptable datasets into equation (20) and adopt the genetic algorithm for acquiring the ${\epsilon }_{\mathrm{multistage}}$ . Compute this repeatedly utilizing Monte Carlo simulation, and analyze the results statistically. The entire optimization procedure has been realized by a MATLAB program. In the case, the numbers of simulation are selected as 10,000 times to balance efficiency with accuracy.

• Step 7: Analyze the results statistically. To prove the accuracy and practicability of the improved J-T theory, Table 8 has compared different eccentricity results, which are based on improved J-T model, original J-T model, and traditional plane dimension chain model, respectively.

Table 7.

Variations and constraints of FE₁.

No.	FE s	Variations	Constraints
1	IFE 1	$\begin{array}{c}-0.012/2\le u_1\le 0.012/2\hfill \\ -0.012/2\le v_1\le 0.012/2\hfill \\ -0.020/2\le w_1\le 0.020/2\hfill \\ -0.006/315\le {\alpha }_1\le 0.006/315\hfill \\ -0.006/315\le {\beta }_1\le 0.006/315\hfill \end{array}$	$\begin{array}{c}-0.020/2\le w_1+{\beta }_1·x_1+{\alpha }_1·y_1\le 0.020/2\hfill \\ -0.015/2\le {\beta }_1·x_1+{\alpha }_1·y_1\le 0.015/2\hfill \\ {(v_1+{\alpha }_1·z_1)}^2+{(u_1+{\beta }_1·z_1)}^2\le {(0.012/2)}^2\hfill \end{array}$
2	CFE 1	0	0

Table 8.

Statistical results of eccentricity of three variation models (Unit: mm).

Object	Improved J-T model	Original J-T model	Plane dimension chain model
Mean	0.1052	0.1052	0
Standard deviation	0.01002	0.01210	0.01102
Range	0.07514, 0.1352	0.06890, 0.1415	−0.03306, 0.03306

To be specific, the improved J-T model can take interactions of composite tolerances into consideration; original J-T model regards every tolerancing torsor as independent variable, ignoring the mutual effects; plane dimension chain model is a common engineering way based on root sum square (RSS) method, which cannot fully consider geometric tolerances in 2D dimension chain.

Comparison and discussion

Statistical eccentricity analysis

Figure 11 has illustrated the results of Table 8. As can be seen, the variation range of improved J-T model is minimum, that of the plane dimension chain model comes second, and the original J-T model is maximum. This phenomenon is closely related to the standard deviation of the three models. The standard deviation of original J-T model (σ₂) is higher than that of improved J-T model (σ₁) of 17.2%, while that of plane dimension chain model (σ₃) is higher than σ₁ of 9.1%. The cause is due to the ineffective management of composite tolerances and lack of study on the common effects between various factors in the original J-T model and plane dimension chain model. Specifically, in the results of σ₂ and σ₃, it is inaccurate to take limit values of both translational and rotational vectors into calculation. However, σ₁ has fully embodied the constraint conditions of composite tolerances, which is much more realistic.

Figure 11.

Statistical results of three variation analysis models.

Meanwhile, Table.9 has provided the statistical results ε_i of each stage rotors. In the table, M-O, T-O, and N-O represent three different assembly strategies of multistage optimization, top-stage optimization and no optimization. The M-O strategy is based on the method proposed in the article, T-O is based on step-by-step strategy focusing on the assembly precision of each top rotor, and N-O is direct-assembly using original J-T model.

Table.9.

Statistical results of eccentricity of three assembly strategies (Unit: mm).

Object	Result	M-O	T-O	N-O
Stage2 (ε2)	Max variation	0.02032	0.05304	0.05304
	Min variation	0.009045	0.009058	0.009045
	Mean	0.009775	0.03386	0.03517
	Standard deviation	0.0001186	0.0002119	0.007001
Stage3 (ε3)	Max variation	0.03792	0.06393	0.09651
	Min variation	0.03029	0.05254	0.0006855
	Mean	0.03434	0.06290	0.05293
	Standard deviation	0.0002023	0.0005803	0.008022
Stage4 (ε4)	Max variation	0.03312	0.01046	0.1595
	Min variation	0.02564	5.216E-09	0.0008930
	Mean	0.02866	0.0001023	0.07812
	Standard deviation	0.0003178	0.001872	0.01000
Entirety (εt)	Max variation	0.05007	0.08236	0.1938
	Min variation	0.04564	0.05433	0.04574
	Mean	0.04580	0.07220	0.1052
	Standard deviation	0.0003500	0.005955	0.01002

The statistical characteristics of the concentricity deviation shown above are illustrated in Figure 12, involving the mean and standard deviation (see Figure 12(a) and (b)). For the M-O result, its mean value of ε_t is 0.0458 mm, about 56.5% better than the N-O result; likewise, the concentricity of the second, third, and fourth-stage rotor has been also improved by 72.2%, 35.1%, and 63.3%, respectively from N-O to M-O.

Figure 12.

Comprehensive statistical results of eccentricity: (a) mean, and (b) standard deviation.

As described in Figure 7, T-O strategy seems to be not valid to ensure the concentricity of the overall assembly. For the results in four-stage rotors assembly, the average concentricity of the second stage component shows outstanding performance with almost zero eccentricity, seeing Figure 12(a). However, the concentricity worsens at other stages. The concentricity is degraded by 54.4% at stage four and 36.6% at the overall assembly, compared to the M-O result. The above results indicate that the M-O strategy based on the proposed model is efficient and feasible, which could be used in actual assembly process.

As can be seen from Figure 12(b), under the M-O strategy, the minimum value of standard deviation of ε_t can be achieved, which is 1/17 of T-O result and 1/28 of N-O result. Similar consequences (σ₂, σ₃, σ₄) also happen at the second, third, and fourth stage respectively. This demonstrates that the M-O strategy can make the concentricity level of each part more concentrated, and reduce the systematic risk of assembly failure. It is also observed that the standard deviation has increased as stacking process continues. For instance, the standard deviation of M-O result has increased 41.4% from stage two to stage three, and 36.3% from stage three to stage four. These results suggest that the concentricity controllability and stability become worse, when the rotor gets closer to the top, and the problem of assembly inconsistency is more prominent. Especially, the assembly consistency should be paid more attention in numerous rotors assembly.

The above results illustrate that the multistage assembly strategy based on improved J-T model is conducive to concentricity prediction and control. Meanwhile, this method ensured that the concentricity variation of each rotor has a good consistency and stability, with a higher degree of accuracy.

Probability density distribution analysis

The probability density distribution of assembly variations could be computed on the basis of improved J-T model, and the effect of variation probability type on the ultimate assembling precision can be verified. To validate this, the probability model of FEs are supposed to obey normal or abnormal distributions respectively. Table 10 has shown four combinations of variation distribution model, which will be calculated in the case.

Table 10.

Distribution type and parameters of the FEs.

No.	Distribution pattern	Distribution parameter
		u 1/v1/w1/α1/β1	u 2/v2/w2/α2/β2	u 3/v3/w3/α3/β3	u 4/v4/w4/α4/β4
1	Normal	$\begin{array}{c}\mu =0\hfill \\ \sigma =\frac{0.005}{405}\hfill \end{array}$	$\begin{array}{c}\mu =0\hfill \\ \sigma =\frac{0.005}{410}\hfill \end{array}$	$\begin{array}{c}\mu =0\hfill \\ \sigma =\frac{0.005}{400}\hfill \end{array}$	$\begin{array}{c}\mu =0\hfill \\ \sigma =\frac{0.005}{440}\hfill \end{array}$
2	Beta (Right-skewed)	$\begin{array}{c}\alpha =2\hfill \\ \beta =5\hfill \end{array}$	$\begin{array}{c}\alpha =3\hfill \\ \beta =4\hfill \end{array}$	$\begin{array}{c}\alpha =2\hfill \\ \beta =8\hfill \end{array}$	$\begin{array}{c}\alpha =1.5\hfill \\ \beta =3.2\hfill \end{array}$
3	Beta (Left-skewed)	$\begin{array}{c}\alpha =5\hfill \\ \beta =2\hfill \end{array}$	$\begin{array}{c}\alpha =4\hfill \\ \beta =3\hfill \end{array}$	$\begin{array}{c}\alpha =8\hfill \\ \beta =2\hfill \end{array}$	$\begin{array}{c}\alpha =3.2\hfill \\ \beta =1.5\hfill \end{array}$
4	Mixing (Normal, Beta, Uniform)	$\begin{array}{c}\mu =0\hfill \\ \sigma =\frac{0.005}{405}\hfill \end{array}$	$\begin{array}{c}\alpha =5\hfill \\ \beta =2\hfill \end{array}$	$\begin{array}{c}\alpha =2\hfill \\ \beta =5\hfill \end{array}$	$\begin{array}{c}a=\frac{-0.015}{440}\hfill \\ b=\frac{0.015}{440}\hfill \end{array}$

Based on the flow chart in Figure 8, FR’s variation probability type could be solved. Table 11 has shown the results of variation distribution parameters.

Table 11.

Distribution type and parameters of the FR.

Parameter	Combination 1	Combination 2	Combination 3	Combination 4
Mean	0	0.0604	0.1396	0.0457
Variance	0.0125	0.0189	0.0189	0.0115
Skewness	0	0.3734	−0.3734	0.2417
Kurtosis	3	2.8210	2.8210	2.9311

With the aid of Monte Carlo (MC) simulation, the variation distribution of the FR can be worked out. The simulation numbers of MC is selected as 10,000 times to maintain efficient and accurate. We assume that the qualified range of FR is (0, 0.038) mm, in accordance with engineering requirements, and the qualification rate of four combinations can be calculated. Figure 13 has illustrated the distribution curve as well as the qualification rate.

Figure 13.

Variation distribution curve and qualification rate: (a) combination 1, (b) combination 2, (c) combination 3, and (d) combination 4.

It is observed that the probability shape of FR’s variation has brought about great changes while the FE’s variation probability type changes. This phenomenon also appears in the aspect of assembly qualified rate. In above four groups of simulations, normal distribution (combination 1) could ensure the qualified rate of up to 99.8%, and the left-skewed distribution (combination 3) almost causes the whole rotors to be unqualified. The reason for the difference is that part’s precision has an overlarge mean shift (=0.1396 mm). All above results have proved that the variation probability type of rotors plays an important role on the ultimate assembly accuracy, and the conventional assumption about normal distribution could not reveal the real variation situation. So the realistically processing capability and assembly technology level is suggested to be considered during the product design phase, to analyze the rotor’s practical variation distribution. Only this can make the results more accurate and credible.

Percentage contribution analysis

Contribution analysis is an important aspect for variation control. It is an effective way for quantificationally observing and controlling the variation propagation. For aero-engine rotors assembly, how to effectively classify parts’ percentage contribution (PC) is the key to improve the overall quality of assembly concentricity. Comprehensively analyzing the contribution effects of each stage rotors can make it clear that which rotor has the most important effect on the precision of the final assembly, so as to avoid repairing or reworking that is caused by the excessive attention to the invalid component’s precision and the incorrect control of the effective component’s precision. For the traditional variation model, the contribution equation could be written as:

$$\begin{array}{c}\mathrm{Extremum}\mathrm{condition}:\mathrm{\Delta }Y=\hfill \\ \left|\frac{\partial f}{\partial X_1}\right|·\mathrm{\Delta }{X}_1+\left|\frac{\partial f}{\partial X_2}\right|·\mathrm{\Delta }{X}_2+···+\left|\frac{\partial f}{\partial X_n}\right|·\mathrm{\Delta }{X}_n\hfill \end{array}$$

(23)

$$\begin{array}{c}\mathrm{Statistical}\mathrm{condition}:\mathrm{\Delta }Y=\hfill \\ \sqrt{{(\frac{\partial f}{\partial X_1}·\mathrm{\Delta }{X}_1)}^2+{(\frac{\partial f}{\partial X_2}·\mathrm{\Delta }{X}_2)}^2+···+{(\frac{\partial f}{\partial X_n}·\mathrm{\Delta }{X}_n)}^2}\hfill \end{array}$$

(24)

where Y denotes the closed link FR; X_i denotes the ith component link FE; ∂f/∂X_i denotes the sensitivity coefficient of the ith component link FE (i = 1, 2, …, n).

Therefore, the ith FE’s contribution is the ratio of the tolerance value of FE and FR, by multiplying the sensitivity coefficient. The specific expression is as follows:

$$P{C}_i=\frac{\frac{\partial f}{\partial X_i}·\mathrm{\Delta }{X}_i}{\mathrm{\Delta }Y}$$

(25)

Similarly, the contribution of the component link based on J-T model with extremum condition can be deduced, which is expressed as:

$$P{C}_{\mathrm{FEid}}=\frac{\mathrm{F}{\mathrm{R}}_{\mathrm{i}}}{\mathrm{FR}}$$

(26)

where FR can be considered as the sum of subfunctional requirements FR_i (1≤i≤n), which have been derived from equation (4). If each FE’s deviations have been determined, equation (26) can be also regarded as the method to solve the percentage contribution for deterministic deviation problem.

The extremum method is aimed at 100% acceptance rate and parts’ absolute interchangeability. However, in actual assembly, parts’ variations are distributed in different FEs in terms of a certain probability. The statistical method allows a certain reject rate and the existence of incomplete interchangeability. It provides a more relaxed variation environment for the rotor system. This section uses 6σ theorem of the quality management system, which considers that 99.73% variations are in ±3σ range, and its rejection rate is 0.27%.

Similar to the statistical conditional equation (23), the component’s contribution based on improved J-T model in the statistical way is expressed as:

$$P{C}_{\mathrm{FEis}}=\frac{{\sigma }_{F{R}_i}^2}{{\sigma }_{\mathrm{FR}}^2}$$

(27)

where σ_FR and σ_FRi denote the standard deviations of FR and FR_i; and the σ_FR equals to the root sum square of the total FR_is’σ_FRi, that is $\mathrm{FR}=\sqrt{{\mathrm{FR}}_1^2+\cdots {+\mathrm{FR}}_{\mathrm{i}}^2+\cdots {+\mathrm{FR}}_{\mathrm{n}}^2}$ .

Based on the above analysis, equations (26) and (27) give the contribution calculation method of each component on the basis of J-T model under the extremum condition and statistical condition, respectively. In this case, take the multistage optimization results which have been worked out in Table 9 into equation (27) to calculate the PC _FEis . The derived PC and APC (accumulative percentage contribution) have been normalized and shown in Figure 14 in the order of their importance.

Figure 14.

PCs of seven FEs and their APCs.

As can be seen, the IFE₄ and IFE₃ are the most important factors in the four-stage rotors assembly. The APC of these two FEs is up to 80.25%. This indicates that strictly restricting the first two FEs’ tolerance and appropriately relaxing the nonsignificant FEs’ tolerance would be an applicable way to realize the balance between the quality and the cost for rotors designing. The result is quite helpful for tolerance synthesis and allocation. In addition, it is observed that the PCs of CFEs are all zero, just because in this case the RJ is not given any geometric tolerances.

Conclusion

This article attempts to present a valid variation analysis method for revolving components assembly with multi joints. Extended Jacobian matrix and updated torsor equation were unified. A calculating scheme was illustrated on a realistic case originating in aero-engine subassembly. Following conclusions could be drawn:

1. The improved J-T model can solve both of variation propagation problem and variation expression issue in revolving components assembly. The extended Jacobian matrix describes the revolving characteristic of the rotor by introducing the revolution joint, and the updated torsor equation characterizes essential mating features of the rotor, involving planar feature and cylindrical feature.

2. The improved J-T model is convenient for multistage concentricity prediction and control under statistical conditions. By adopting multistage rotational optimization, the overall concentric performance has improved 36.6% relative to the top-stage optimization, and 56.5% relative to directly assembling. Statistical results highlighted the benefits that this approach has a good consistency with the stage numbers, and a better practicality.

3. This improved model fits probability density analysis and contribution analysis statistically. The result of FR significantly depends on each FE’s variation distribution. In the four-stage rotors assembly, restricting the upper two FEs’ tolerances (APC = 80.25%) and relaxing the rest of FEs’ tolerances is a feasible way to reach a balance between the cost and quality.

Moreover, shortcoming of the proposed method is that, for the rotor structure with the locating spigot round and bolted connection, deformation of working parts under load would be difficult to deal with. In future, the actual contact deformation under certain pre-tension from flange joint bolt and the available number of practical orientations should be taken into account, to model and optimize the actual assembly process, and an effective probabilistic method is needed to improve simulation efficiency, which can be advantageous to guide assembling fulfillment.