# Tolerance Analysis of Over-Constrained Assembly Considering Gravity Influence: Constraints of Multiple Planar Hole-Pin-Hole Pairs

Authors: Xia Liu; Luling An; Zhiguo Wang; Changbai Tan; Xiaoping Wang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Over-constrained assembly of rigid parts is widely adopted in aircraft assembly to yield higher stiffness and accuracy of assembly. Unfortunately, the quantitative tolerance analysis of over-constrained assembly is challenging, subject to the coupling effect of geometrical and physical factors. Especially, gravity will affect the geometrical gaps in mechanical joints between different parts, and thus influence the deviations of assembled product. In the existing studies, the influence of gravity is not considered in the tolerance analysis of over-constrained assembly. This paper proposes a novel tolerance analysis method for over-constrained assembly of rigid parts, considering the gravity influence. This method is applied to a typical over-constrained assembly with constraints of multiple planar hole-pin-hole pairs. This type of constraints is non-linear, which makes the tolerance analysis more challenging. Firstly, the deviation propagation analysis of an over-constrained assembly is conducted. The feasibility of assembly is predicted, and for a feasible assembly, the assembly deviations are determined with the principle of minimum potential energy. Then, the statistical tolerance analysis is performed. The probabilities of assembly feasibility and quality feasibility are computed, and the distribution of assembly deviations is estimated. Two case studies are presented to show the applicability of the proposed method.


---

## Body

## 1. Introduction

Mechanical parts are always manufactured with more or less geometrical deviations. Tolerance analysis aims to predict the influence of these manufacturing deviations on the quality level of a mechanical system at the design stage. The gaps in mechanical joints between different parts, coupled with the manufacturing deviations, impact the functionality of an assembly system. The manufacturing deviations are assumed to be independent random variables defined by their tolerance zones and the probability distributions with respect to manufacturing process, while the gaps are characteristic of epistemic uncertainty which makes the tolerance analysis complex [1]. The gaps are seldom considered in iso-constrained assembly, but in over-constrained assembly they allow the parts to be assembled without deformation and interference, and cannot be ignored. Over-constrained assembly of rigid parts is widely used in products, such as aircrafts and automobiles, because this kind of assembly can increase stiffness and accuracy of mechanical structure [2]. However, the feasibility of assembly is weakened due to the introduction of redundant constraints into the assembly, and the deviations of assembled products are difficult to predict because of the uncertain gaps. This is one of the key challenges in assembly tolerance analysis.Some studies have been conducted on tolerance analysis of over-constrained assembly. Compared to worst case method, statistical tolerance analysis is regarded to be better to understand the behavior of a mechanical system, and the over-constraints have to be considered for accuracy in the tolerance analysis of over-constrained assembly [3]. Dantan and Qureshi [4] introduced the quantifier notion, including existential quantifier “there exists” and universal quantifier “for all”, to describe the condition corresponding to a functional requirement. The mathematical formulation of tolerance analysis was computed with quantified constraint satisfaction problem solvers and Monte Carlo simulation (MCS). With the mathematical formulation based on quantifier notion, Qureshi et al. [5] acquired the probability of the functional operation of an assembly based on MCS and optimization technique, which is used to find the worst case for gap configuration. This approach was generalized and applicable to both explicit and implicit response function, as well as both linear and non-linear constraints. For explicit function and linear constraints, MCS is time consuming and computationally expensive. Ballu et al. [3] solved probability of failure by using first order reliability method (FORM) and second-order reliability method (SORM) to determine the feasible extreme gap configuration. Beaucaire et al. [6, 7] proposed the mathematical formulation from the structural reliability domain considering several contact point situations, and computed the defect probability by FORM. For non-linear constraints, Dumas et al. [8] investigated different linearization strategies and discussed the impact of these strategies on the predicted probability of failure. In addition, Dumas et al. [1] formulated the tolerance analysis problem by using the Lagrange dual form of optimization, where a simulation-based selective search algorithm, i.e., a genetic algorithm, was developed and then the probability of failure was estimated by FORM. The above studies [1, 3–8] focused on computing the probability that the assembled product can fulfill its intended functional requirement.Another problem related to statistical tolerance analysis is the distributions of assembly deviations of functional points, where these deviations characterize the fulfillment of functional requirement of assembled product. In over-constrained assembly, the presence of gaps in mechanical joints makes the deviations of assembled product difficult to predict, because a large variety of possible gap configurations exist during the assembly process. Franciosa et al. [9] determined the best-fit assembly configuration by using optimization algorithm, where the optimization objective is to minimize the sum of the squares of the distances between object plane/axis and target plane/axis respectively. However, this optimization objective is not suitable for most of practical assemblies in production. Li et al. [10] defined two planar hole-pin pairs as 3-2-1 location mode, and expressed the ranges of location deviations by the intersection of inequations which is approximated to an elliptic. However, this study did not develop an approach to get the certain location deviations.Furthermore, tolerance analysis of over-constrained assembly is subject to many nongeometrical factors involved in an assembly process, such as thermal flux, contact force, gravity, etc. Pierre et al. [11] integrated thermomechanical strains into variational tolerance analysis. Gouyou et al. [2] discussed the part deformation subject to contact forces which allow the parts to be assembled when the interference happens. However, although the influence of gravity on the assembly deviations of rigid parts has been well recognized [12], there is still limited research in tolerance analysis of over-constrained assembly integrating the influence of gravity. Specifically, in aircraft assembly, planar hole-pin-hole joints are widely used to join different parts and locate parts on the fixtures. Because many parts and fixtures in aircraft assembly are large and heavy, the gravity will significantly impact the configuration of gaps in hole-pin-hole joints, and thus the assembly accuracy. In view of the strict dimensional accuracy requirement in aircraft manufacturing, a reliable prediction of assembly deviations by considering the gravity influence is especially important to realize effective dimensional management in aircraft assembly.This paper proposes a method for tolerance analysis of over-constrained assembly incorporating the influence of gravity on assembly deviations. This method will be applied to a typical over-constrained assembly with constraints of multiple planar hole-pin-hole pairs. First of all, the deviation propagation throughout an over-constrained assembly is analyzed. Here, the feasibility of assembly is predicted; for a feasible assembly, the assembly deviations are determined considering the gravity influence. Then, the statistical tolerance analysis is performed by using MCS, and its output results include the probabilities of assembly feasibility and quality feasibility, as well as the distributions of assembly deviations. Especially, in the deviation propagation analysis, the judgment of assembly feasibility is formulated as a search problem. And for a feasible assembly, with the principle of minimum potential energy, the calculation of assembly deviations is modeled as an optimization problem. These search and optimization problems are inefficient and difficult to solve, because of the non-linear studied constraints of planar hole-pin-hole pairs. These problems are solved by using the search algorithms including sequential search and binary search, and the optimization algorithms including particle swarm optimization (PSO) and linear programming (LP).The paper is organized as follows. Section2 gives problem definition and mathematical representation. Sections 3 and 4 present deviation propagation analysis and statistical tolerance analysis of over-constrained assembly, respectively. Section 5 shows case studies. Section 6 draws conclusions.


## 2. Problem Definition and Mathematical Representation

### 2.1. Problem Description and Assumption

Tolerance analysis consists in predicting the functionality of an assembly system based on the consideration of the manufacturing deviations of parts. For over-constrained assembly of rigid parts, tolerance analysis consists of two key problems to be studied, including assembly feasibility prediction and assembly deviation calculation. In over-constrained assembly of rigid parts, the feasibility of assembly is weak because of redundant constraints, even with the gaps in mechanical joints between different parts. Furthermore, due to the uncertainty of gaps, the assembly deviations are difficult to predict, especially coupled with non-linear constraints and non-geometrical factors.With considering the influence of gravity, the proposed method of tolerance analysis of over-constrained assembly is illustrated by a typical assembly with constraints of multiple planar hole-pin-hole pairs. This type of constraints possesses the feature of non-linearity. As shown in Figure1, each of the two parts, denoted as Part 1 and Part 2, has two or more pairs of mating holes on the mating plane. Part 1 is assumed to be fixed, and Part 2 will be assembled onto Part 1 with the planar hole-pin-hole joints. The centers of the mating holes of Part 1 are denoted as B1~BN. Similarly, these centers of Part 2 are denoted as C1~CN.N is the number of hole-pin-hole joints. Furthermore, the quality requirement concerns the deviations of Part 2 with respect to (w.r.t.) Part 1. The functional points are on Part 2, denoted as K1~KNf.Nf is the number of functional points. The deviations of the functional points are used to represent the assembly deviations. In addition, G denotes the center of gravity of Part 2.Figure 1
Two parts to be assembled with constraints of multiple planar hole-pin-hole pairs.Notably, for simplifying the analysis and highlighting the research focus, some assumptions are made as follows:(a) The assembly tolerance analysis is performed only considering the degrees of freedom (DOFs) on plane, such that the mating planes are assumed to be perfect.(b) The parts are treated as rigid bodies, and deformation and interference are not allowed.(c) The direction of gravitational field is on the mating planes.(d) The position of gravity center of Part 2 is not affected by its manufacturing deviations.(e) The hole-pin-hole pairs adopt clearance fit.(f) For the manufacturing deviations of hole-pin-hole pairs, dimensional deviations of diameters of holes and pins, and position deviations of centers of holes are considered, but the shape deviations are ignored.


### 2.2. Mathematical Representation

#### 2.2.1. Definition of Coordinate Systems

In order to express the deviations as well as their relationships, four types of coordinate systems (CSs) are defined, including global coordinate system (GCS), part coordinate system (PCS), gravitational field coordinate system (GFCS), and functional point coordinate system (FPCS).GCS is considered to be error free and unchangeable during the entire assembly process. PCS is associated with a part. Similarly to GCS, GFCS and FPCS are also unchanged. GFCS is constructed according to the direction of gravitational field, and theYg-axis is in the opposite direction to gravitational field. Each FPCS is used to represent the deviations of a functional point. The origin of a FPCS is in the nominal position of a functional point w.r.t. GCS, and the axes correspond to the direction required to be measured. Moreover, the spatial relationships among different CSs can be described by homogeneous transformation matrix (HTM) [13].As shown in Figure2, Part 2 without any deviation is in the nominal position w.r.t Part 1.XOYis the GCS, which is assumed to coincide with the PCS of Part 1.xoy is the PCS of Part 2 denoted as PCS2. XgOgYg is the GFCS, whose origin is in the nominal position of G w.r.t. GCS, and Yg-axis has the angle θ fromY-axis of GCS. XfjOfjYfj, j=1,2,…,Nf, are the FPCSs. Furthermore, the HTMs of nominal PCS2, GFCS and FPCSi w.r.t. GCS are T2G, TgG, and TfjG, j=1,2,…,Nf, respectively.Figure 2
Definition of coordinate systems.


#### 2.2.2. Tolerance Representation

The dimensional and geometric tolerances indicate the deviation of an actual feature from its nominal one. The features are allowed to rotate and translate within their tolerance zone [10], and these rotations and translations can be described by the small displacement torsor (SDT) model [14]. In SDT model, tolerance zones and feature deviations are defined as vectorial representation [15], i.e., τ=ωεT. ω=αβγT denotes the rotation vector, which has three components α, β, and γ aroundx,y, andz axes, respectively. Similarly, ε=uvwT denotes the translation vector, which has three componentsu,v, and w alongx,y andz axes, respectively. In this paper, it is assumed that α=β=0, w=0, thus τ=γuvT.The assembly deviations are caused by deviation sources which are propagated and accumulated through the assembly process. This data stream of deviations can be analyzed by deviation propagation model [16–18]. For example, the deviation of featurei is τi=γiuiviT. An arbitrary point on featurei is A, and its coordinates are denoted as a=xayaT. The deviation of point A can be calculated with HTM [19, 20] as(1)Δa1T=Tia1T-a1Twhere the HTM Tiis(2)Ti=Riεi01=cos⁡γi-sin⁡γiuisin⁡γicos⁡γivi001.


## 2.1. Problem Description and Assumption

Tolerance analysis consists in predicting the functionality of an assembly system based on the consideration of the manufacturing deviations of parts. For over-constrained assembly of rigid parts, tolerance analysis consists of two key problems to be studied, including assembly feasibility prediction and assembly deviation calculation. In over-constrained assembly of rigid parts, the feasibility of assembly is weak because of redundant constraints, even with the gaps in mechanical joints between different parts. Furthermore, due to the uncertainty of gaps, the assembly deviations are difficult to predict, especially coupled with non-linear constraints and non-geometrical factors.With considering the influence of gravity, the proposed method of tolerance analysis of over-constrained assembly is illustrated by a typical assembly with constraints of multiple planar hole-pin-hole pairs. This type of constraints possesses the feature of non-linearity. As shown in Figure1, each of the two parts, denoted as Part 1 and Part 2, has two or more pairs of mating holes on the mating plane. Part 1 is assumed to be fixed, and Part 2 will be assembled onto Part 1 with the planar hole-pin-hole joints. The centers of the mating holes of Part 1 are denoted as B1~BN. Similarly, these centers of Part 2 are denoted as C1~CN.N is the number of hole-pin-hole joints. Furthermore, the quality requirement concerns the deviations of Part 2 with respect to (w.r.t.) Part 1. The functional points are on Part 2, denoted as K1~KNf.Nf is the number of functional points. The deviations of the functional points are used to represent the assembly deviations. In addition, G denotes the center of gravity of Part 2.Figure 1
Two parts to be assembled with constraints of multiple planar hole-pin-hole pairs.Notably, for simplifying the analysis and highlighting the research focus, some assumptions are made as follows:(a) The assembly tolerance analysis is performed only considering the degrees of freedom (DOFs) on plane, such that the mating planes are assumed to be perfect.(b) The parts are treated as rigid bodies, and deformation and interference are not allowed.(c) The direction of gravitational field is on the mating planes.(d) The position of gravity center of Part 2 is not affected by its manufacturing deviations.(e) The hole-pin-hole pairs adopt clearance fit.(f) For the manufacturing deviations of hole-pin-hole pairs, dimensional deviations of diameters of holes and pins, and position deviations of centers of holes are considered, but the shape deviations are ignored.


## 2.2. Mathematical Representation

### 2.2.1. Definition of Coordinate Systems

In order to express the deviations as well as their relationships, four types of coordinate systems (CSs) are defined, including global coordinate system (GCS), part coordinate system (PCS), gravitational field coordinate system (GFCS), and functional point coordinate system (FPCS).GCS is considered to be error free and unchangeable during the entire assembly process. PCS is associated with a part. Similarly to GCS, GFCS and FPCS are also unchanged. GFCS is constructed according to the direction of gravitational field, and theYg-axis is in the opposite direction to gravitational field. Each FPCS is used to represent the deviations of a functional point. The origin of a FPCS is in the nominal position of a functional point w.r.t. GCS, and the axes correspond to the direction required to be measured. Moreover, the spatial relationships among different CSs can be described by homogeneous transformation matrix (HTM) [13].As shown in Figure2, Part 2 without any deviation is in the nominal position w.r.t Part 1.XOYis the GCS, which is assumed to coincide with the PCS of Part 1.xoy is the PCS of Part 2 denoted as PCS2. XgOgYg is the GFCS, whose origin is in the nominal position of G w.r.t. GCS, and Yg-axis has the angle θ fromY-axis of GCS. XfjOfjYfj, j=1,2,…,Nf, are the FPCSs. Furthermore, the HTMs of nominal PCS2, GFCS and FPCSi w.r.t. GCS are T2G, TgG, and TfjG, j=1,2,…,Nf, respectively.Figure 2
Definition of coordinate systems.


### 2.2.2. Tolerance Representation

The dimensional and geometric tolerances indicate the deviation of an actual feature from its nominal one. The features are allowed to rotate and translate within their tolerance zone [10], and these rotations and translations can be described by the small displacement torsor (SDT) model [14]. In SDT model, tolerance zones and feature deviations are defined as vectorial representation [15], i.e., τ=ωεT. ω=αβγT denotes the rotation vector, which has three components α, β, and γ aroundx,y, andz axes, respectively. Similarly, ε=uvwT denotes the translation vector, which has three componentsu,v, and w alongx,y andz axes, respectively. In this paper, it is assumed that α=β=0, w=0, thus τ=γuvT.The assembly deviations are caused by deviation sources which are propagated and accumulated through the assembly process. This data stream of deviations can be analyzed by deviation propagation model [16–18]. For example, the deviation of featurei is τi=γiuiviT. An arbitrary point on featurei is A, and its coordinates are denoted as a=xayaT. The deviation of point A can be calculated with HTM [19, 20] as(1)Δa1T=Tia1T-a1Twhere the HTM Tiis(2)Ti=Riεi01=cos⁡γi-sin⁡γiuisin⁡γicos⁡γivi001.


## 2.2.1. Definition of Coordinate Systems

In order to express the deviations as well as their relationships, four types of coordinate systems (CSs) are defined, including global coordinate system (GCS), part coordinate system (PCS), gravitational field coordinate system (GFCS), and functional point coordinate system (FPCS).GCS is considered to be error free and unchangeable during the entire assembly process. PCS is associated with a part. Similarly to GCS, GFCS and FPCS are also unchanged. GFCS is constructed according to the direction of gravitational field, and theYg-axis is in the opposite direction to gravitational field. Each FPCS is used to represent the deviations of a functional point. The origin of a FPCS is in the nominal position of a functional point w.r.t. GCS, and the axes correspond to the direction required to be measured. Moreover, the spatial relationships among different CSs can be described by homogeneous transformation matrix (HTM) [13].As shown in Figure2, Part 2 without any deviation is in the nominal position w.r.t Part 1.XOYis the GCS, which is assumed to coincide with the PCS of Part 1.xoy is the PCS of Part 2 denoted as PCS2. XgOgYg is the GFCS, whose origin is in the nominal position of G w.r.t. GCS, and Yg-axis has the angle θ fromY-axis of GCS. XfjOfjYfj, j=1,2,…,Nf, are the FPCSs. Furthermore, the HTMs of nominal PCS2, GFCS and FPCSi w.r.t. GCS are T2G, TgG, and TfjG, j=1,2,…,Nf, respectively.Figure 2
Definition of coordinate systems.


## 2.2.2. Tolerance Representation

The dimensional and geometric tolerances indicate the deviation of an actual feature from its nominal one. The features are allowed to rotate and translate within their tolerance zone [10], and these rotations and translations can be described by the small displacement torsor (SDT) model [14]. In SDT model, tolerance zones and feature deviations are defined as vectorial representation [15], i.e., τ=ωεT. ω=αβγT denotes the rotation vector, which has three components α, β, and γ aroundx,y, andz axes, respectively. Similarly, ε=uvwT denotes the translation vector, which has three componentsu,v, and w alongx,y andz axes, respectively. In this paper, it is assumed that α=β=0, w=0, thus τ=γuvT.The assembly deviations are caused by deviation sources which are propagated and accumulated through the assembly process. This data stream of deviations can be analyzed by deviation propagation model [16–18]. For example, the deviation of featurei is τi=γiuiviT. An arbitrary point on featurei is A, and its coordinates are denoted as a=xayaT. The deviation of point A can be calculated with HTM [19, 20] as(1)Δa1T=Tia1T-a1Twhere the HTM Tiis(2)Ti=Riεi01=cos⁡γi-sin⁡γiuisin⁡γicos⁡γivi001.


## 3. Deviation Propagation Analysis

Deviation propagation analysis studies the introduction, propagation and accumulation of data stream of deviations through an assembly process. For an over-constrained assembly with constraints of multiple planar hole-pin-hole pairs, the deviation propagation analysis is performed with considering the influence of gravity, which contains assembly feasibility judgment and assembly deviation determination.


### 3.1. Assembly Feasibility Judgment

An assembly with over constraints is usually disturbed by the infeasibility of assembly, even though normal joint gaps are reserved. Dantan and Qureshi [4] introduced the existential quantifier “there exists” to correctly formulate this problem concerning assembly with gaps, i.e., if there exits at least one gap configuration such that the assembly is feasible, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is(3)∃g∈0,gmaxM:faiM,g≤0,i=1,2,⋯,NwhereM denotes the manufacturing deviations of parts, g denotes the gaps and gmax is the allowable gap widths,fai, i=1,2,…,N, are the functions of assembly requirement, andN is the number of the functions. All of these functions have to be non-positive to ensure the feasibility of assembly.According to the above formulation, the judgment of assembly feasibility can be transformed to a search problem. Moreover, this search problem is inefficient to be solved due to the non-linear studied constraints of multiple planar hole-pin-hole pairs. Firstly, the search problem is modeled. Then in order to improve the efficiency, a solving method is given by reducing dimensionality of search object and narrowing down search space.


#### 3.1.1. Formulation of Assembly Requirement

In over-constrained assembly, the gaps in mechanical joints are correlated with each other and their worst configuration depends on the manufacturing deviations of parts [3]. For the constraints of multiple planar hole-pin-hole pairs, the correlation of gaps, coupled with the relationships between gap configurations and part manufacturing deviations, is complicated. In an assembly process, the deviation of PCS2 w.r.t GCS is denoted as τ2=γ2u2v2T, and one τ2 corresponds to a configuration of gaps.As shown in Figure1, the coordinates of Bi, i=1,2,…,N, w.r.t. GCS are BiG=XbiYbiT. The coordinates of Ci, i=1,2,…,N, w.r.t. PCS2 are denoted as Ci2=xciyciT. Then when Part 2 is in the nominal position w.r.t. GCS, they w.r.t. GCS are CiG=XciYciT as(4)CiG1T=T2GCi21T,i=1,2,⋯,NWith the deviationτ2, the coordinates of Ci′, i=1,2,…,N, w.r.t. GCS are changed into CiG′=Xci′Yci′T as(5)CiG′1T=T2CiG1T,i=1,2,⋯,Nwhere T2 is the HTM of τ2.The functions of assembly requirement concern the planar hole-pin-hole joints. As shown in Figure3, by taking theith constraint of planar hole-pin-hole pair as an example, the function of assembly requirement is(6)fai=BiCi′→-siwhere si=(1/2)(dbi+dci-2dpi) anddbi,dci, anddpi denote the diameters of the two mating holes and the pin respectively. It is definite that si≥0, due to clearance fit of hole-pin-hole pairs. According to the assembly requirement as seen in (3), the judgment of assembly feasibility is to judge whetherfai is non-positive, i.e., to judge whether BiCi′→, is not larger thansi, or whether |BiCi′→|2 is not larger than si2. Therefore, from the view of judging assembly feasibility, (6) can be translated to (7) as(7)f^ai=BiCi′→2-si2where f^ai is decorated with a topscript ‘‘∧’’, in order to distinguish from the original assembly requirement functionfai.Figure 3
Assembly of two parts under a constraint of planar hole-pin-hole pair.By referring to (5), (7) can be further arranged as(8)f^ai=u2+Xcicos⁡γ2-Ycisin⁡γ2-Xbi2+v2+Xcisin⁡γ2+Ycicos⁡γ2-Ybi2-si2For coordinates and diameters above mentioned, the nominal values are decorated with a superscript ‘‘0,’’ whereas they represent the actual ones with manufacturing deviations of parts. In (8), coupled with the correspondence between τ2 and configuration of gaps, the part manufacturing deviations and gap configuration are both incorporated.According to (3) and (8), the assembly requirement under the constraints ofN planar hole-pin-hole pairs is formulated as: if there exits at least one real vector τ2 satisfying all the inequalities f^ai≤0, i=1,2,…,N, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is(9)∃τ2∈R3×1∩γ2∈-π,π:f^ai=u2+Xcicos⁡γ2-Ycisin⁡γ2-Xbi2+v2+Xcisin⁡γ2+Ycicos⁡γ2-Ybi2-si2≤0,i=1,2,⋯,Nwhere R3×1 denotes three-dimensional real vector space.


#### 3.1.2. Solving Method of Search Problem

The judgment of assembly feasibility depends on solving the search problem modeled by (9). As a common search algorithm, sequential search algorithm is used. The search object is τ2=γ2u2v2T with three dimensions. The search space is difficult to narrow, due to the non-linearity of the functions f^ai, i=1,2,…,N. Therefore, the solving of this search problem is inefficient because of the multi-dimensional search object and the large search space. To improve the efficiency, the dimensionality of search object is reduced into one, and the search space is reasonably narrowed.The feasible solution space of the search objectτ2 can be shown in a Cartesian coordinate system o2-u2v2γ2. Referring to Figure 4, the feasible solution space for an inequality f^ai≤0 is denoted asDfai. When γ2 is constant,Dfai is a certain circle, and its radius always equals to si. Besides, the locus of center of this circle with γ2 is a right-hand circular helix. Then the feasible solution space under constraints ofN planar hole-pin-hole pairs, denoted asDfa_in, can be expressed as the intersection of the spacesDfai, i=1,2,…,N. Figure 5 shows an example of non-emptyDfa_in whenN=2.Figure 4
Feasible solution space under constraint of a planar hole-pin-hole pair.Figure 5
Feasible solution space under constraints of two planar hole-pin-hole pairs.
(a)
Two feasible solution spacesDfa1 andDfa2 (b)
Intersection of the spacesDfa1 andDfa2Whenγ2 is constant,Dfa_in represents the intersection ofN circles with the known centers and radii. The situation of this intersection denoted as Sfa_in(γ2) can be judged with a geometry algorithm [21] especially whenN>2. It is assumed that when the intersection is empty, Sfa_in(γ2)<0; when the intersection is only a point, Sfa_in(γ2)=0; otherwise, Sfa_in(γ2)>0. Then (9) can be redefined as follows: if there exits at least one real γ2 such that the intersection ofN circles is non-empty, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is(10)∃γ2∈Dse_γ2:Sfa_inγ2≥0where Dse_γ2={γ2∣γ2∈[γ2se_min,γ2se_max]} denotes the search space of γ2, which is determined as follows:(a) The ranges ofγ2 under constraints of each two planar hole-pin-hole pairs are successively calculated. For example, as shown in Figure 5, the range of γ2 is denoted as γ2∈[γ2_12_min,γ2_12_max]. According to geometric characteristic ofDfa1 andDfa2, γ2_12_min, and γ2_12_min correspond to the situations when the two circles f^a1(γ2_12_min,max,u2,v2)≤0 and f^a2(γ2_12_min,max,u2,v2)≤0 are externally tangent as(11)-a2sin⁡γ2_12_min,max+a1cos⁡γ2_12_min,max+a32+a1sin⁡γ2_12_min,max+a2cos⁡γ2_12_min,max+a42-a52=0where a1=Xc2-Xc1, a2=Yc2-Yc1, a3=Xb1-Xb2, a4=Yb1-Yb2, and a5=s1+s2.If (11) has no real solutions, this range of γ2 is empty.(b) The intersection of these ranges ofγ2 in (a) is calculated as Dse_γ2.As shown in Figure6, if Dse_γ2 is empty, it is sure that the assembly is infeasible; but if Dse_γ2 is non-empty, the feasibility of assembly is still not ensured and the search algorithm has to be performed, unlessN=2.Figure 6
Flow chart of solving search problem.


### 3.2. Assembly Deviation Determination

Assembly deviations are difficult to quantify in over-constrained mechanical structure due to uncertainty of gaps, when the assembly is feasible. Here, the influence of gravity is considered. According to the principle of minimum potential energy, the gravitational potential energy of Part 2 tends to be minimum. In other words, the center of gravity of Part 2 drops along the direction of gravitational field as large as possible, within the constraints of satisfying assembly requirement. Therefore this section involves an optimization problem. Furthermore, the non-linear studied constraints make the optimization difficult to solve. Firstly the optimization problem is modeled. Then similarly to Section3.1, a solving method is given by reducing dimensionality of optimization object and narrowing search space, to improve both accuracy and efficiency. Finally, the assembly deviations of functional points are calculated, and the quality feasibility is judged.


#### 3.2.1. Optimization Problem Modeling

As shown in Figure1, the coordinates of G w.r.t. PCS2 are denoted as G2=xgygT, and it is assumed that G2=G02. When Part 2 is in the nominal position and has the deviation τ2 w.r.t. GCS, the coordinates of G w.r.t. GCS are denoted as GG=XgYgT and GG′=Xg′Yg′T respectively, which can be calculated similarly to (4)-(5), as(12)GG1T=T2GG21T(13)GG′1T=T2GG1TThen with the deviationτ2, the coordinates of G w.r.t. GFCS are Gg′=Xgg′Ygg′T as(14)Gg′1T=TgG-1GG′1THere,Yg-axis of GFCS has the angle θ fromY-axis of GCS. Then (14) can be arranged as(15)Xgg′=Xgcos⁡γ2-Ygsin⁡γ2+u2-Xgcos⁡θ+Xgsin⁡γ2+Ygcos⁡γ2+v2-Ygsin⁡θYgg′=-Xgcos⁡γ2-Ygsin⁡γ2+u2-Xgsin⁡θ+Xgsin⁡γ2+Ygcos⁡γ2+v2-Ygcos⁡θwhere Ygg′ is the coordinate of G against the direction of gravitational field.When the influence of gravity is considered,Ygg′ is minimum within the constraints of satisfying assembly requirement, according to the principle of minimum potential energy. The optimization model is constructed as(16)minFfdτ2=Ygg′τ2s.t.τ2∈Dfa_inwhere τ2 is the optimization object, Ffd(τ2) is the objective function, and τ2∈Dfa_in is the constraint condition.Without considering the influence of gravity and other physical factors, the deviationτ2 is randomly distributed within the feasible solution spaceDfa_in. When the influence of gravity on the assembly deviations is taken into account (the influence of other physical factors is not considered), the value of τ2 corresponds to the optimal solution of (16). Figure 7 shows an example of the assembly under constraints of two planar hole-pin-hole pairs. With considering the gravity influence, the center of gravity G on Part 2 moves along the direction of gravitational field to the farthest position allowed by the assembly constraints.Figure 7
Assembly under constraints of two planar hole-pin-hole pairs with considering gravity influence.


#### 3.2.2. Optimization Problem Solving

The assembly deviations correspond to the optimal solution of (16). As a selective search algorithm, PSO algorithm is used to solve the optimization problem. This optimization solving is difficult due to the non-linear constraint τ2∈Dfa_in. Here, the dimensionality of search object is reduced, and the search space is calculated, to improve the accuracy and efficiency of solving.As seen in (16), when γ2 equals a real constant γ2_con, the constraint τ2∈Dfa_in(γ2_con,u2,v2) represents the intersection ofN known circles by referring to Section 3.1, and the objective function Ffd(γ2_con,u2,v2) is regarded as a line. The slope of this line is unchangeable, which is only related to the direction of gravitational field. The line equation denoted asl(u2,v2) is(17)lu2,v2=au2+bv2+cwhere a=-sin⁡θ, b=cos⁡θ, c=(-Xgsin⁡θ+Ygcos⁡θ)(cos⁡γ2_con-1)+Ygsin⁡θ+Xgcos⁡θsin⁡γ2_con.Therefore, the optimization modeled by (16) can be treated a LP problem when γ2=γ2_con. Figure 8 gives a graphical schematic of constraint condition and objective function whenN=3. Here, the minimum value of Ffd(γ2_con,u2,v2) is denoted as Ffd_min(γ2_con,u2,v2). The coordinates of the point corresponding to Ffd_min(γ2_con,u2,v2) are denoted as PFfd_min=u2_Ffd_minv2_Ffd_minT. Obviously, this point belongs to intersection points of the arcs which comprise boundary of the domain Dfa_in(γ2_con,u2,v2), or tangency points of the linel(u2,v2) and these arcs. Furthermore, it is certain that PFfd_min depends on γ2_con, i.e., u2_Ffd_min=u2(γ2_con), v2_Ffd_min=v2(γ2_con). Hence, (16) can be remodeled with a nested LP as(18)minFfd_minγ2,u2γ2,v2γ2s.t.γ2∈Dfa_in_γ2where the search space Dfa_in_γ2={γ2∣γ2∈[γ2min,γ2max]} is the range of γ2 inDfa_in, which is calculated as follows: As seen in Section 3.1, a feasible solution of (10), denoted as γ2_a, has been found to ensure the assembly feasibility of the mechanical structure. It is obvious that γ2_a∈[γ2min,γ2max]. Moreover, for the search space of (10) denoted as Dse_γ2={γ2∣γ2∈[γ2se_min,γ2se_max]}, it is definite that Dse_γ2⊇Dfa_in_γ2, and then γ2se_min≤γ2min≤γ2_a and γ2_a≤γ2max≤γ2se_max. γ2min and γ2max can be approximately calculated by binary search algorithm. Especially, whenN=2, γ2min=γ2se_min and γ2max=γ2se_max, thus the binary search algorithm is omitted.Figure 8
Constraint condition and objective function with constantγ2.


#### 3.2.3. Assembly Deviation Calculation and Quality Feasibility Judgment

By the solving method in Section3.2.2, the optimal solutions of the optimization model (18), as well as the global best particle τ2_g, are determined. Then the coordinates of functional point Kj w.r.t. FPCSj, j=1,2,…,Nf, are calculated as(19)Kj_gfj′1T=TfjG-1T2_gKjG1T,j=1,2,⋯,Nfwhere T2_g=T2(τ2_g) and Kj_gfj′=Xfj_gfj′Yfj_gfj′T, j=1,2,…,Nf, represent the assembly deviations.The quality requirement is denoted asQR=QR1QR2⋯QRNf, where QRj=QRjXfjQRjYfjT and QRjXfj={XQRj∣XQRj∈[XQRj_min,XQRj_max]}, QRjYfj={YQRj∣YQRj∈[YQRj_min,YQRj_max]}, j=1,2,…,Nf. An assembly is regarded to be qualified when all the assembly deviations of functional points meet the quality requirement, i.e., Xfj_gfj′∈QRjXfj and Yfj_gfj′∈QRjYfj, j=1,2,…,Nf.


## 3.1. Assembly Feasibility Judgment

An assembly with over constraints is usually disturbed by the infeasibility of assembly, even though normal joint gaps are reserved. Dantan and Qureshi [4] introduced the existential quantifier “there exists” to correctly formulate this problem concerning assembly with gaps, i.e., if there exits at least one gap configuration such that the assembly is feasible, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is(3)∃g∈0,gmaxM:faiM,g≤0,i=1,2,⋯,NwhereM denotes the manufacturing deviations of parts, g denotes the gaps and gmax is the allowable gap widths,fai, i=1,2,…,N, are the functions of assembly requirement, andN is the number of the functions. All of these functions have to be non-positive to ensure the feasibility of assembly.According to the above formulation, the judgment of assembly feasibility can be transformed to a search problem. Moreover, this search problem is inefficient to be solved due to the non-linear studied constraints of multiple planar hole-pin-hole pairs. Firstly, the search problem is modeled. Then in order to improve the efficiency, a solving method is given by reducing dimensionality of search object and narrowing down search space.


### 3.1.1. Formulation of Assembly Requirement

In over-constrained assembly, the gaps in mechanical joints are correlated with each other and their worst configuration depends on the manufacturing deviations of parts [3]. For the constraints of multiple planar hole-pin-hole pairs, the correlation of gaps, coupled with the relationships between gap configurations and part manufacturing deviations, is complicated. In an assembly process, the deviation of PCS2 w.r.t GCS is denoted as τ2=γ2u2v2T, and one τ2 corresponds to a configuration of gaps.As shown in Figure1, the coordinates of Bi, i=1,2,…,N, w.r.t. GCS are BiG=XbiYbiT. The coordinates of Ci, i=1,2,…,N, w.r.t. PCS2 are denoted as Ci2=xciyciT. Then when Part 2 is in the nominal position w.r.t. GCS, they w.r.t. GCS are CiG=XciYciT as(4)CiG1T=T2GCi21T,i=1,2,⋯,NWith the deviationτ2, the coordinates of Ci′, i=1,2,…,N, w.r.t. GCS are changed into CiG′=Xci′Yci′T as(5)CiG′1T=T2CiG1T,i=1,2,⋯,Nwhere T2 is the HTM of τ2.The functions of assembly requirement concern the planar hole-pin-hole joints. As shown in Figure3, by taking theith constraint of planar hole-pin-hole pair as an example, the function of assembly requirement is(6)fai=BiCi′→-siwhere si=(1/2)(dbi+dci-2dpi) anddbi,dci, anddpi denote the diameters of the two mating holes and the pin respectively. It is definite that si≥0, due to clearance fit of hole-pin-hole pairs. According to the assembly requirement as seen in (3), the judgment of assembly feasibility is to judge whetherfai is non-positive, i.e., to judge whether BiCi′→, is not larger thansi, or whether |BiCi′→|2 is not larger than si2. Therefore, from the view of judging assembly feasibility, (6) can be translated to (7) as(7)f^ai=BiCi′→2-si2where f^ai is decorated with a topscript ‘‘∧’’, in order to distinguish from the original assembly requirement functionfai.Figure 3
Assembly of two parts under a constraint of planar hole-pin-hole pair.By referring to (5), (7) can be further arranged as(8)f^ai=u2+Xcicos⁡γ2-Ycisin⁡γ2-Xbi2+v2+Xcisin⁡γ2+Ycicos⁡γ2-Ybi2-si2For coordinates and diameters above mentioned, the nominal values are decorated with a superscript ‘‘0,’’ whereas they represent the actual ones with manufacturing deviations of parts. In (8), coupled with the correspondence between τ2 and configuration of gaps, the part manufacturing deviations and gap configuration are both incorporated.According to (3) and (8), the assembly requirement under the constraints ofN planar hole-pin-hole pairs is formulated as: if there exits at least one real vector τ2 satisfying all the inequalities f^ai≤0, i=1,2,…,N, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is(9)∃τ2∈R3×1∩γ2∈-π,π:f^ai=u2+Xcicos⁡γ2-Ycisin⁡γ2-Xbi2+v2+Xcisin⁡γ2+Ycicos⁡γ2-Ybi2-si2≤0,i=1,2,⋯,Nwhere R3×1 denotes three-dimensional real vector space.


### 3.1.2. Solving Method of Search Problem

The judgment of assembly feasibility depends on solving the search problem modeled by (9). As a common search algorithm, sequential search algorithm is used. The search object is τ2=γ2u2v2T with three dimensions. The search space is difficult to narrow, due to the non-linearity of the functions f^ai, i=1,2,…,N. Therefore, the solving of this search problem is inefficient because of the multi-dimensional search object and the large search space. To improve the efficiency, the dimensionality of search object is reduced into one, and the search space is reasonably narrowed.The feasible solution space of the search objectτ2 can be shown in a Cartesian coordinate system o2-u2v2γ2. Referring to Figure 4, the feasible solution space for an inequality f^ai≤0 is denoted asDfai. When γ2 is constant,Dfai is a certain circle, and its radius always equals to si. Besides, the locus of center of this circle with γ2 is a right-hand circular helix. Then the feasible solution space under constraints ofN planar hole-pin-hole pairs, denoted asDfa_in, can be expressed as the intersection of the spacesDfai, i=1,2,…,N. Figure 5 shows an example of non-emptyDfa_in whenN=2.Figure 4
Feasible solution space under constraint of a planar hole-pin-hole pair.Figure 5
Feasible solution space under constraints of two planar hole-pin-hole pairs.
(a)
Two feasible solution spacesDfa1 andDfa2 (b)
Intersection of the spacesDfa1 andDfa2Whenγ2 is constant,Dfa_in represents the intersection ofN circles with the known centers and radii. The situation of this intersection denoted as Sfa_in(γ2) can be judged with a geometry algorithm [21] especially whenN>2. It is assumed that when the intersection is empty, Sfa_in(γ2)<0; when the intersection is only a point, Sfa_in(γ2)=0; otherwise, Sfa_in(γ2)>0. Then (9) can be redefined as follows: if there exits at least one real γ2 such that the intersection ofN circles is non-empty, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is(10)∃γ2∈Dse_γ2:Sfa_inγ2≥0where Dse_γ2={γ2∣γ2∈[γ2se_min,γ2se_max]} denotes the search space of γ2, which is determined as follows:(a) The ranges ofγ2 under constraints of each two planar hole-pin-hole pairs are successively calculated. For example, as shown in Figure 5, the range of γ2 is denoted as γ2∈[γ2_12_min,γ2_12_max]. According to geometric characteristic ofDfa1 andDfa2, γ2_12_min, and γ2_12_min correspond to the situations when the two circles f^a1(γ2_12_min,max,u2,v2)≤0 and f^a2(γ2_12_min,max,u2,v2)≤0 are externally tangent as(11)-a2sin⁡γ2_12_min,max+a1cos⁡γ2_12_min,max+a32+a1sin⁡γ2_12_min,max+a2cos⁡γ2_12_min,max+a42-a52=0where a1=Xc2-Xc1, a2=Yc2-Yc1, a3=Xb1-Xb2, a4=Yb1-Yb2, and a5=s1+s2.If (11) has no real solutions, this range of γ2 is empty.(b) The intersection of these ranges ofγ2 in (a) is calculated as Dse_γ2.As shown in Figure6, if Dse_γ2 is empty, it is sure that the assembly is infeasible; but if Dse_γ2 is non-empty, the feasibility of assembly is still not ensured and the search algorithm has to be performed, unlessN=2.Figure 6
Flow chart of solving search problem.


## 3.1.1. Formulation of Assembly Requirement

In over-constrained assembly, the gaps in mechanical joints are correlated with each other and their worst configuration depends on the manufacturing deviations of parts [3]. For the constraints of multiple planar hole-pin-hole pairs, the correlation of gaps, coupled with the relationships between gap configurations and part manufacturing deviations, is complicated. In an assembly process, the deviation of PCS2 w.r.t GCS is denoted as τ2=γ2u2v2T, and one τ2 corresponds to a configuration of gaps.As shown in Figure1, the coordinates of Bi, i=1,2,…,N, w.r.t. GCS are BiG=XbiYbiT. The coordinates of Ci, i=1,2,…,N, w.r.t. PCS2 are denoted as Ci2=xciyciT. Then when Part 2 is in the nominal position w.r.t. GCS, they w.r.t. GCS are CiG=XciYciT as(4)CiG1T=T2GCi21T,i=1,2,⋯,NWith the deviationτ2, the coordinates of Ci′, i=1,2,…,N, w.r.t. GCS are changed into CiG′=Xci′Yci′T as(5)CiG′1T=T2CiG1T,i=1,2,⋯,Nwhere T2 is the HTM of τ2.The functions of assembly requirement concern the planar hole-pin-hole joints. As shown in Figure3, by taking theith constraint of planar hole-pin-hole pair as an example, the function of assembly requirement is(6)fai=BiCi′→-siwhere si=(1/2)(dbi+dci-2dpi) anddbi,dci, anddpi denote the diameters of the two mating holes and the pin respectively. It is definite that si≥0, due to clearance fit of hole-pin-hole pairs. According to the assembly requirement as seen in (3), the judgment of assembly feasibility is to judge whetherfai is non-positive, i.e., to judge whether BiCi′→, is not larger thansi, or whether |BiCi′→|2 is not larger than si2. Therefore, from the view of judging assembly feasibility, (6) can be translated to (7) as(7)f^ai=BiCi′→2-si2where f^ai is decorated with a topscript ‘‘∧’’, in order to distinguish from the original assembly requirement functionfai.Figure 3
Assembly of two parts under a constraint of planar hole-pin-hole pair.By referring to (5), (7) can be further arranged as(8)f^ai=u2+Xcicos⁡γ2-Ycisin⁡γ2-Xbi2+v2+Xcisin⁡γ2+Ycicos⁡γ2-Ybi2-si2For coordinates and diameters above mentioned, the nominal values are decorated with a superscript ‘‘0,’’ whereas they represent the actual ones with manufacturing deviations of parts. In (8), coupled with the correspondence between τ2 and configuration of gaps, the part manufacturing deviations and gap configuration are both incorporated.According to (3) and (8), the assembly requirement under the constraints ofN planar hole-pin-hole pairs is formulated as: if there exits at least one real vector τ2 satisfying all the inequalities f^ai≤0, i=1,2,…,N, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is(9)∃τ2∈R3×1∩γ2∈-π,π:f^ai=u2+Xcicos⁡γ2-Ycisin⁡γ2-Xbi2+v2+Xcisin⁡γ2+Ycicos⁡γ2-Ybi2-si2≤0,i=1,2,⋯,Nwhere R3×1 denotes three-dimensional real vector space.


## 3.1.2. Solving Method of Search Problem

The judgment of assembly feasibility depends on solving the search problem modeled by (9). As a common search algorithm, sequential search algorithm is used. The search object is τ2=γ2u2v2T with three dimensions. The search space is difficult to narrow, due to the non-linearity of the functions f^ai, i=1,2,…,N. Therefore, the solving of this search problem is inefficient because of the multi-dimensional search object and the large search space. To improve the efficiency, the dimensionality of search object is reduced into one, and the search space is reasonably narrowed.The feasible solution space of the search objectτ2 can be shown in a Cartesian coordinate system o2-u2v2γ2. Referring to Figure 4, the feasible solution space for an inequality f^ai≤0 is denoted asDfai. When γ2 is constant,Dfai is a certain circle, and its radius always equals to si. Besides, the locus of center of this circle with γ2 is a right-hand circular helix. Then the feasible solution space under constraints ofN planar hole-pin-hole pairs, denoted asDfa_in, can be expressed as the intersection of the spacesDfai, i=1,2,…,N. Figure 5 shows an example of non-emptyDfa_in whenN=2.Figure 4
Feasible solution space under constraint of a planar hole-pin-hole pair.Figure 5
Feasible solution space under constraints of two planar hole-pin-hole pairs.
(a)
Two feasible solution spacesDfa1 andDfa2 (b)
Intersection of the spacesDfa1 andDfa2Whenγ2 is constant,Dfa_in represents the intersection ofN circles with the known centers and radii. The situation of this intersection denoted as Sfa_in(γ2) can be judged with a geometry algorithm [21] especially whenN>2. It is assumed that when the intersection is empty, Sfa_in(γ2)<0; when the intersection is only a point, Sfa_in(γ2)=0; otherwise, Sfa_in(γ2)>0. Then (9) can be redefined as follows: if there exits at least one real γ2 such that the intersection ofN circles is non-empty, then assembly feasibility of the mechanical structure is ensured. The mathematical expression is(10)∃γ2∈Dse_γ2:Sfa_inγ2≥0where Dse_γ2={γ2∣γ2∈[γ2se_min,γ2se_max]} denotes the search space of γ2, which is determined as follows:(a) The ranges ofγ2 under constraints of each two planar hole-pin-hole pairs are successively calculated. For example, as shown in Figure 5, the range of γ2 is denoted as γ2∈[γ2_12_min,γ2_12_max]. According to geometric characteristic ofDfa1 andDfa2, γ2_12_min, and γ2_12_min correspond to the situations when the two circles f^a1(γ2_12_min,max,u2,v2)≤0 and f^a2(γ2_12_min,max,u2,v2)≤0 are externally tangent as(11)-a2sin⁡γ2_12_min,max+a1cos⁡γ2_12_min,max+a32+a1sin⁡γ2_12_min,max+a2cos⁡γ2_12_min,max+a42-a52=0where a1=Xc2-Xc1, a2=Yc2-Yc1, a3=Xb1-Xb2, a4=Yb1-Yb2, and a5=s1+s2.If (11) has no real solutions, this range of γ2 is empty.(b) The intersection of these ranges ofγ2 in (a) is calculated as Dse_γ2.As shown in Figure6, if Dse_γ2 is empty, it is sure that the assembly is infeasible; but if Dse_γ2 is non-empty, the feasibility of assembly is still not ensured and the search algorithm has to be performed, unlessN=2.Figure 6
Flow chart of solving search problem.


## 3.2. Assembly Deviation Determination

Assembly deviations are difficult to quantify in over-constrained mechanical structure due to uncertainty of gaps, when the assembly is feasible. Here, the influence of gravity is considered. According to the principle of minimum potential energy, the gravitational potential energy of Part 2 tends to be minimum. In other words, the center of gravity of Part 2 drops along the direction of gravitational field as large as possible, within the constraints of satisfying assembly requirement. Therefore this section involves an optimization problem. Furthermore, the non-linear studied constraints make the optimization difficult to solve. Firstly the optimization problem is modeled. Then similarly to Section3.1, a solving method is given by reducing dimensionality of optimization object and narrowing search space, to improve both accuracy and efficiency. Finally, the assembly deviations of functional points are calculated, and the quality feasibility is judged.


### 3.2.1. Optimization Problem Modeling

As shown in Figure1, the coordinates of G w.r.t. PCS2 are denoted as G2=xgygT, and it is assumed that G2=G02. When Part 2 is in the nominal position and has the deviation τ2 w.r.t. GCS, the coordinates of G w.r.t. GCS are denoted as GG=XgYgT and GG′=Xg′Yg′T respectively, which can be calculated similarly to (4)-(5), as(12)GG1T=T2GG21T(13)GG′1T=T2GG1TThen with the deviationτ2, the coordinates of G w.r.t. GFCS are Gg′=Xgg′Ygg′T as(14)Gg′1T=TgG-1GG′1THere,Yg-axis of GFCS has the angle θ fromY-axis of GCS. Then (14) can be arranged as(15)Xgg′=Xgcos⁡γ2-Ygsin⁡γ2+u2-Xgcos⁡θ+Xgsin⁡γ2+Ygcos⁡γ2+v2-Ygsin⁡θYgg′=-Xgcos⁡γ2-Ygsin⁡γ2+u2-Xgsin⁡θ+Xgsin⁡γ2+Ygcos⁡γ2+v2-Ygcos⁡θwhere Ygg′ is the coordinate of G against the direction of gravitational field.When the influence of gravity is considered,Ygg′ is minimum within the constraints of satisfying assembly requirement, according to the principle of minimum potential energy. The optimization model is constructed as(16)minFfdτ2=Ygg′τ2s.t.τ2∈Dfa_inwhere τ2 is the optimization object, Ffd(τ2) is the objective function, and τ2∈Dfa_in is the constraint condition.Without considering the influence of gravity and other physical factors, the deviationτ2 is randomly distributed within the feasible solution spaceDfa_in. When the influence of gravity on the assembly deviations is taken into account (the influence of other physical factors is not considered), the value of τ2 corresponds to the optimal solution of (16). Figure 7 shows an example of the assembly under constraints of two planar hole-pin-hole pairs. With considering the gravity influence, the center of gravity G on Part 2 moves along the direction of gravitational field to the farthest position allowed by the assembly constraints.Figure 7
Assembly under constraints of two planar hole-pin-hole pairs with considering gravity influence.


### 3.2.2. Optimization Problem Solving

The assembly deviations correspond to the optimal solution of (16). As a selective search algorithm, PSO algorithm is used to solve the optimization problem. This optimization solving is difficult due to the non-linear constraint τ2∈Dfa_in. Here, the dimensionality of search object is reduced, and the search space is calculated, to improve the accuracy and efficiency of solving.As seen in (16), when γ2 equals a real constant γ2_con, the constraint τ2∈Dfa_in(γ2_con,u2,v2) represents the intersection ofN known circles by referring to Section 3.1, and the objective function Ffd(γ2_con,u2,v2) is regarded as a line. The slope of this line is unchangeable, which is only related to the direction of gravitational field. The line equation denoted asl(u2,v2) is(17)lu2,v2=au2+bv2+cwhere a=-sin⁡θ, b=cos⁡θ, c=(-Xgsin⁡θ+Ygcos⁡θ)(cos⁡γ2_con-1)+Ygsin⁡θ+Xgcos⁡θsin⁡γ2_con.Therefore, the optimization modeled by (16) can be treated a LP problem when γ2=γ2_con. Figure 8 gives a graphical schematic of constraint condition and objective function whenN=3. Here, the minimum value of Ffd(γ2_con,u2,v2) is denoted as Ffd_min(γ2_con,u2,v2). The coordinates of the point corresponding to Ffd_min(γ2_con,u2,v2) are denoted as PFfd_min=u2_Ffd_minv2_Ffd_minT. Obviously, this point belongs to intersection points of the arcs which comprise boundary of the domain Dfa_in(γ2_con,u2,v2), or tangency points of the linel(u2,v2) and these arcs. Furthermore, it is certain that PFfd_min depends on γ2_con, i.e., u2_Ffd_min=u2(γ2_con), v2_Ffd_min=v2(γ2_con). Hence, (16) can be remodeled with a nested LP as(18)minFfd_minγ2,u2γ2,v2γ2s.t.γ2∈Dfa_in_γ2where the search space Dfa_in_γ2={γ2∣γ2∈[γ2min,γ2max]} is the range of γ2 inDfa_in, which is calculated as follows: As seen in Section 3.1, a feasible solution of (10), denoted as γ2_a, has been found to ensure the assembly feasibility of the mechanical structure. It is obvious that γ2_a∈[γ2min,γ2max]. Moreover, for the search space of (10) denoted as Dse_γ2={γ2∣γ2∈[γ2se_min,γ2se_max]}, it is definite that Dse_γ2⊇Dfa_in_γ2, and then γ2se_min≤γ2min≤γ2_a and γ2_a≤γ2max≤γ2se_max. γ2min and γ2max can be approximately calculated by binary search algorithm. Especially, whenN=2, γ2min=γ2se_min and γ2max=γ2se_max, thus the binary search algorithm is omitted.Figure 8
Constraint condition and objective function with constantγ2.


### 3.2.3. Assembly Deviation Calculation and Quality Feasibility Judgment

By the solving method in Section3.2.2, the optimal solutions of the optimization model (18), as well as the global best particle τ2_g, are determined. Then the coordinates of functional point Kj w.r.t. FPCSj, j=1,2,…,Nf, are calculated as(19)Kj_gfj′1T=TfjG-1T2_gKjG1T,j=1,2,⋯,Nfwhere T2_g=T2(τ2_g) and Kj_gfj′=Xfj_gfj′Yfj_gfj′T, j=1,2,…,Nf, represent the assembly deviations.The quality requirement is denoted asQR=QR1QR2⋯QRNf, where QRj=QRjXfjQRjYfjT and QRjXfj={XQRj∣XQRj∈[XQRj_min,XQRj_max]}, QRjYfj={YQRj∣YQRj∈[YQRj_min,YQRj_max]}, j=1,2,…,Nf. An assembly is regarded to be qualified when all the assembly deviations of functional points meet the quality requirement, i.e., Xfj_gfj′∈QRjXfj and Yfj_gfj′∈QRjYfj, j=1,2,…,Nf.


## 3.2.1. Optimization Problem Modeling

As shown in Figure1, the coordinates of G w.r.t. PCS2 are denoted as G2=xgygT, and it is assumed that G2=G02. When Part 2 is in the nominal position and has the deviation τ2 w.r.t. GCS, the coordinates of G w.r.t. GCS are denoted as GG=XgYgT and GG′=Xg′Yg′T respectively, which can be calculated similarly to (4)-(5), as(12)GG1T=T2GG21T(13)GG′1T=T2GG1TThen with the deviationτ2, the coordinates of G w.r.t. GFCS are Gg′=Xgg′Ygg′T as(14)Gg′1T=TgG-1GG′1THere,Yg-axis of GFCS has the angle θ fromY-axis of GCS. Then (14) can be arranged as(15)Xgg′=Xgcos⁡γ2-Ygsin⁡γ2+u2-Xgcos⁡θ+Xgsin⁡γ2+Ygcos⁡γ2+v2-Ygsin⁡θYgg′=-Xgcos⁡γ2-Ygsin⁡γ2+u2-Xgsin⁡θ+Xgsin⁡γ2+Ygcos⁡γ2+v2-Ygcos⁡θwhere Ygg′ is the coordinate of G against the direction of gravitational field.When the influence of gravity is considered,Ygg′ is minimum within the constraints of satisfying assembly requirement, according to the principle of minimum potential energy. The optimization model is constructed as(16)minFfdτ2=Ygg′τ2s.t.τ2∈Dfa_inwhere τ2 is the optimization object, Ffd(τ2) is the objective function, and τ2∈Dfa_in is the constraint condition.Without considering the influence of gravity and other physical factors, the deviationτ2 is randomly distributed within the feasible solution spaceDfa_in. When the influence of gravity on the assembly deviations is taken into account (the influence of other physical factors is not considered), the value of τ2 corresponds to the optimal solution of (16). Figure 7 shows an example of the assembly under constraints of two planar hole-pin-hole pairs. With considering the gravity influence, the center of gravity G on Part 2 moves along the direction of gravitational field to the farthest position allowed by the assembly constraints.Figure 7
Assembly under constraints of two planar hole-pin-hole pairs with considering gravity influence.


## 3.2.2. Optimization Problem Solving

The assembly deviations correspond to the optimal solution of (16). As a selective search algorithm, PSO algorithm is used to solve the optimization problem. This optimization solving is difficult due to the non-linear constraint τ2∈Dfa_in. Here, the dimensionality of search object is reduced, and the search space is calculated, to improve the accuracy and efficiency of solving.As seen in (16), when γ2 equals a real constant γ2_con, the constraint τ2∈Dfa_in(γ2_con,u2,v2) represents the intersection ofN known circles by referring to Section 3.1, and the objective function Ffd(γ2_con,u2,v2) is regarded as a line. The slope of this line is unchangeable, which is only related to the direction of gravitational field. The line equation denoted asl(u2,v2) is(17)lu2,v2=au2+bv2+cwhere a=-sin⁡θ, b=cos⁡θ, c=(-Xgsin⁡θ+Ygcos⁡θ)(cos⁡γ2_con-1)+Ygsin⁡θ+Xgcos⁡θsin⁡γ2_con.Therefore, the optimization modeled by (16) can be treated a LP problem when γ2=γ2_con. Figure 8 gives a graphical schematic of constraint condition and objective function whenN=3. Here, the minimum value of Ffd(γ2_con,u2,v2) is denoted as Ffd_min(γ2_con,u2,v2). The coordinates of the point corresponding to Ffd_min(γ2_con,u2,v2) are denoted as PFfd_min=u2_Ffd_minv2_Ffd_minT. Obviously, this point belongs to intersection points of the arcs which comprise boundary of the domain Dfa_in(γ2_con,u2,v2), or tangency points of the linel(u2,v2) and these arcs. Furthermore, it is certain that PFfd_min depends on γ2_con, i.e., u2_Ffd_min=u2(γ2_con), v2_Ffd_min=v2(γ2_con). Hence, (16) can be remodeled with a nested LP as(18)minFfd_minγ2,u2γ2,v2γ2s.t.γ2∈Dfa_in_γ2where the search space Dfa_in_γ2={γ2∣γ2∈[γ2min,γ2max]} is the range of γ2 inDfa_in, which is calculated as follows: As seen in Section 3.1, a feasible solution of (10), denoted as γ2_a, has been found to ensure the assembly feasibility of the mechanical structure. It is obvious that γ2_a∈[γ2min,γ2max]. Moreover, for the search space of (10) denoted as Dse_γ2={γ2∣γ2∈[γ2se_min,γ2se_max]}, it is definite that Dse_γ2⊇Dfa_in_γ2, and then γ2se_min≤γ2min≤γ2_a and γ2_a≤γ2max≤γ2se_max. γ2min and γ2max can be approximately calculated by binary search algorithm. Especially, whenN=2, γ2min=γ2se_min and γ2max=γ2se_max, thus the binary search algorithm is omitted.Figure 8
Constraint condition and objective function with constantγ2.


## 3.2.3. Assembly Deviation Calculation and Quality Feasibility Judgment

By the solving method in Section3.2.2, the optimal solutions of the optimization model (18), as well as the global best particle τ2_g, are determined. Then the coordinates of functional point Kj w.r.t. FPCSj, j=1,2,…,Nf, are calculated as(19)Kj_gfj′1T=TfjG-1T2_gKjG1T,j=1,2,⋯,Nfwhere T2_g=T2(τ2_g) and Kj_gfj′=Xfj_gfj′Yfj_gfj′T, j=1,2,…,Nf, represent the assembly deviations.The quality requirement is denoted asQR=QR1QR2⋯QRNf, where QRj=QRjXfjQRjYfjT and QRjXfj={XQRj∣XQRj∈[XQRj_min,XQRj_max]}, QRjYfj={YQRj∣YQRj∈[YQRj_min,YQRj_max]}, j=1,2,…,Nf. An assembly is regarded to be qualified when all the assembly deviations of functional points meet the quality requirement, i.e., Xfj_gfj′∈QRjXfj and Yfj_gfj′∈QRjYfj, j=1,2,…,Nf.


## 4. Statistical Tolerance Analysis

Because as deviation sources, manufacturing deviations of parts are generally random deviations within their tolerance zones and statistical analysis is necessary rather than a single case analysis. As one of the most common methods for statistical analysis, the MCS is often computationally expensive and time-consuming, but quite comprehensive and easy to use even for implicit and non-linear problems. The MCS with deviation propagation analysis in Section3 is used to perform statistical tolerance analysis. The flowchart of the statistical analysis is shown in Figure 9.Figure 9
Flow chart of statistical tolerance analysis.Iteratively, a set of random deviations of the sources are generated within their tolerance zones; then the assembly feasibility is judged; for the feasible assembly, the assembly deviations of functional points are calculated and the quality feasibility is judged. The process of random generation, judgment and calculation is iterated until the samples are enough. As to random generation of deviations within their tolerance zones, it is assumed that the deviations are independent of one another and present normal distributions, and the deviations and their tolerances have a relationship asT=6σ.After all the iterations, the probabilities of assembly feasibility and quality feasibility are estimated as [8](20)Pa,q=1NMC∑i=1NMCIa,qMwhereNMC is the maximum number of iterations of the MCS and Ia,q(M) is an indicator function. The functions for assembly requirement and quality requirement are(21)IaM=1,iftheassemblyisfeasible0,iftheassemblyisinfeasible(22)IqM=1,ifthequalityrequirementisfulfilled0,ifthequalityrequirementisunfulfilledFurthermore, for feasible assemblies, the deviations distributions of functional points can be quantified, where the frequency distributions and the statistical characteristics are computed, seen in Section5.For the statistical tolerance analysis, the input is deviation sources, and the output contains probabilities of assembly feasibility and quality feasibility, as well as deviations distributions of functional points for feasible assemblies.


## 5. Case Studies

The proposed method is applied to two cases. Firstly stability and accuracy of the optimization seen in Section3.2 are verified. Then the output results of statistical tolerance analysis for the gap configuration influenced by gravity are compared with the output results for random configuration and worst configuration of gaps.The computational works are implemented by MATLAB. For sequential and binary search algorithms seen in Section3, the precision is set as 10−7 for γ2. In the statistical tolerance analysis with MC simulation seen in Section 4, the number of iterations in an experiment is set asNMC=10000.


### 5.1. Case Description

Case 1.
This case is from a fixture with boards as shown in Figure10, which is widely used for locating aircraft panels in the panel assembly process. The small working surfaces on a fixture board are applied to constraining the normal displacement of the panels. Each fixture board is mounted on fixture frame with two planar hole-pin-hole joints. This case is the assembly of a fixture board and a fixture frame, denoted as Part 1 and Part 2. Each small working surface of Part 2 is abstracted as a functional point.Figure 10
An assembly fixture with boards.Figure11 shows the definition of CSs. For simplicity, the axes of GCS, GFCS, and nominal PCS2 are assumed to be all parallel to each other. As to FPCSs, three of the functional points are selected for discussion. TheYfj-axis (j=1,2,3) of a FPCS is along the nominal normal direction of the corresponding working surface. The dimensions and tolerances are shown in Figure 12. The fit tolerances of hole-pin-hole pairs are both H7/g7/H7. The coordinates of related points w.r.t. GCS/PCS2 are as B1G=195246.126T, B2G=2805246.126T, oG=1500200T, and C12=-130546.126T, C22=130546.126T, G2=0524.435T, Ok2=0-1200T. Besides, the quality requirement is QR=QR1QR2QR3, where QRj=QRjYfj={YQRj∣YQRj∈[-0.07,0.07]},j=1,2,3.Figure 11
CSs of a fixture frame and a fixture board.Figure 12
Dimensions and tolerances of a fixture frame and a fixture board.Case 2.
This case is the assembly of two aircraft fuselage sections from a type of military aircraft. The front section and the rear section, denoted as Part 1 and Part 2, are assembled with the constraints of eighteen hole-pin-hole pairs on planar bulkheads. The quality requirement concerns the displacement deviations of Part 2 w.r.t Part 1.Figure13 shows the CSs, dimensions and tolerances. In order to simplify the computation, it is assumed that GCS, GFCS, FPCS and nominal PCS2 are all coincided with each other. Besides, the fit tolerances of hole-pin-hole pairs are all D8/d7/D8. The quality requirement is QR=QRXfQRYfT, QRXf={XQR∣XQR∈[-0.1,0.1]}, and QRYf={YQR∣YQR∈[-0.1,0.1]}.Figure 13
CSs, Dimensions and tolerances of two fuselage sections.


### 5.2. Stability and Accuracy of Optimization (PSO)

As seen in Section3, PSO is used to determine assembly deviations, which further impacts the output results of statistical tolerance analysis. In this PSO algorithm, the swarm size is 10, and the number of iterations is 20. PSO is a heuristic optimization algorithm, with simple parameters and fast convergence, but it is easily trapped into local optimum. Therefore, it is necessary to verify stability and accuracy of this optimization algorithm.Firstly, the experiment of statistical tolerance analysis is performed and repeated 10 times, and the 10 sets of output results are observed and compared. In order to eliminate the disturbance of different generated deviation sources, the sameNMC sets of deviation sources are required for all experiments.Three iterative sub-processes are separated from the experiment process. The iterative sub-process of deviation sources generation is denoted as Sub-process 1, andNMC sets of generated deviation sources are stored in a matrixHM. Then with these deviation sources, the iterative sub-process of assembly feasibility judgment is denoted as Sub-process 2, andNa (Na≤NMC) sets of deviation sources, which can make sure the feasibility of assembly, are selected and stored in a matrixHa. Finally, with the deviation sources fromHa, the iterative sub-process of assembly deviation calculation and quality feasibility judgment is denoted as Sub-process 3. After the three sub-processes, the output results can be calculated.Therefore, the 10 experiments are designed as: Sub-processes 1 and 2 are implemented in advance, and then Sub-process 3 is repeated 10 times. The 10 sets of output results with 4 figures after the decimal point are all the same, as shown in Table1 for Case 1 and Table 2 for Case 2, where deviations of functional points (denoted as DFP) are classified into the deviations satisfying quality requirement (denoted as DFP_QR) and that not satisfying quality requirement (denoted as DFP_NQR).Table 1
Output results of experiments for Case1.
Output result K1 K2 K3 Probability of quality feasibility (Pq) 80.22% DFP Mean(mm) -0.0208 -0.0235 -0.0243 Standard deviation(mm) 0.0267 0.0255 0.0250 DFP_QR Mean(mm) -0.0177 -0.0204 -0.0213 Standard deviation(mm) 0.0242 0.0231 0.0227Table 2
Output results of experiments for Case2.
Output result X Y Probability of quality feasibility (Pq) 82.62% DFP Mean(mm) -0.0010 -0.0402 Standard deviation(mm) 0.0261 0.0263 DFP_QR Mean(mm) -0.0010 -0.0400 Standard deviation(mm) 0.0262 0.0261For further verification, with a set of deviation sources arbitrarily selected fromHa, the PSO procedure is performed and repeated 10 times, and the 10 sets of optimization results are compared. Tables 3 and 4 show the selected deviation sources for Cases 1 and 2, respectively. The corresponding feasible solution spaceDfa_in is shown in Figure 14 for Case 1 and Figure 15 for Case 2. In the 10 PSO procedures, the global best particles γ2 and the corresponding Yg-coordinates are shown in Figure 16 for Case 1 and Figure 17 for Case 2. It is observed that the 10 optimization procedures for each case are all convergent, and the 10 sets of optimization results, including global best particles and objective function values, are all the same when these values are rounded to 6 decimal places.Table 3
A set of deviation sources for Case1.
Deviation source Coordinates of Bi w.r.t. GCS Coordinates of Ci w.r.t. PCS2 s i (mm) X Y x y Hole B1-pin-Hole C1 194.9787 246.1397 -1305.0518 46.1361 0.0324 Hole B2-pin-Hole C2 2804.9946 246.1302 1305.0297 46.1142 0.0351Table 4
A set of deviation sources for Case2.
Deviation source Coordinates of Bi w.r.t. GCS Coordinates of Ci w.r.t. PCS2 s i (mm) X Y x y Hole B1-pin-Hole C1 450.0372 -0.0106 449.9975 -0.0411 0.1358 Hole B2-pin-Hole C2 422.7972 153.8436 422.8378 153.9293 0.1250 Hole B3-pin-Hole C3 344.7235 289.2612 344.7289 289.2649 0.1273 Hole B4-pin-Hole C4 225.0135 389.7606 225.0130 389.7105 0.1300 Hole B5-pin-Hole C5 78.1725 443.1735 78.1648 443.1337 0.1236 Hole B6-pin-Hole C6 -78.0688 443.2002 -78.1283 443.1470 0.1217 Hole B7-pin-Hole C7 -225.0053 389.6915 -224.9820 389.6781 0.1264 Hole B8-pin-Hole C8 -344.7693 289.2210 -344.6785 289.2782 0.1218 Hole B9-pin-Hole C9 -422.8215 153.9050 -422.8642 153.9120 0.1149 Hole B10-pin-Hole C10 -450.0165 0.0108 -450.0530 0.0023 0.1264 Hole B11-pin-Hole C11 -422.8637 -153.9604 -422.8953 -153.9496 0.1298 Hole B12-pin-Hole C12 -344.7204 -289.2573 -344.6876 -289.2419 0.1211 Hole B13-pin-Hole C13 -224.9730 -389.7477 -224.9676 -389.7089 0.1179 Hole B14-pin-Hole C14 -78.1358 -443.1664 -78.1309 -443.1260 0.1352 Hole B15-pin-Hole C15 78.1410 -443.1581 78.1550 -443.2297 0.1315 Hole B16-pin-Hole C16 224.9878 -389.6686 224.9854 -389.7149 0.1369 Hole B17-pin-Hole C17 344.7059 -289.2523 344.6852 -289.1532 0.1266 Hole B18-pin-Hole C18 422.7669 -153.9629 422.8744 -153.8547 0.1373Figure 14
A feasible solution space for Case1.Figure 15
A feasible solution space for Case2.Figure 16
Global best particles and objective function values of PSO procedures for Case1.
(a)
Global best particles for Case1 (b)
Objective function values for Case1Figure 17
Global best particles and objective function values of PSO procedure for Case2.
(a)
Global best particles for Case2 (b)
Objective function values for Case2Therefore, Tables1 and 2 and Figures 16 and 17 indicate this PSO algorithm is quite robust and accurate, and thus the output results of tolerance analysis by using the proposed method are reliable.


### 5.3. Comparison with Random Configuration and Worst Configuration of Gaps

The proposed method considers the influence of gravity on the gap configuration when determining the assembly deviations. In order to exhibit the significance of incorporating the influence of gravity to tolerance analysis, two comparison experiments with common gap configurations are given. The first comparison experiment is denoted as EX1, where the deviationτ2 is randomly generated within the feasible solution space based on normal distribution. The second comparison experiment, denoted as EX2, uses the worst configuration of gaps, which makes the assembly deviations largest. Besides, the experiment considering the gravity influence is denoted as EX0.In EX 2, the determination of the worst configuration of gaps is an optimization problem. The modeling and solving of this optimization are both similar to (16)-(19). For example, worst gap configuration corresponds to the largest assembly deviation of thejth functional point alongYfj-axis. The optimization is modeled as seen in (23), and solved by PSO with a nested LP. (23)maxFfd_worτ2=Yfjfj′τ2-YQRj_min+YQRj_max2s.t.τ2∈Dfa_inIn EX2 of Case1, it is assumed that the worst configuration of gaps corresponds to the largest assembly deviations of K2 alongYf2-axis. In EX2 of Case 2, the worst configuration of gaps is assumed to concern the largest assembly deviations alongY-axis, to make a straight comparison between the output results of EX2 and EX0.The output results of EX0, EX1 and EX2 are compared. In order to avoid the distractions of different generated deviation sources, the three experiments, denoted as a group of experiments, are designed as follows: Sub-processes 1 and 2 are implemented in advance, and then Sub-process 3 for EX0, EX1, and EX2 is performed, respectively. To enhance the reliability of the predicted output results, a group of experiments are repeated 10 times, and the 10 groups of output results are averaged, which are shown in Table5 for Case 1 and Table 6 for Case 2. Furthermore, by taking the output results of a group of experiments as examples, the distributions of assembly deviations are shown in Figure 18 for Case 1 and Figure 19 for Case 2, and the distributions of the corresponding τ2 are shown in Figure 20 for Case 1 and Figure 21 for Case 2.Table 5
Output results of EX0, EX1, and EX2 for Case1.
Output result  EX0 EX1 EX2 Probability of assembly feasibility (Pa)  87.01% Probability of quality feasibility (Pq)  80.24% 85.26% 73.55% DFP Mean(mm) K1 -0.021 0.000 0.000 K2 -0.023 0.000 0.000 K3 -0.024 0.000 0.000 Standard deviation(mm) K1 0.026 0.026 0.042 K2 0.026 0.025 0.047 K3 0.025 0.025 0.041 DFP_QR Mean(mm) K1 -0.018 0.000 0.000 K2 -0.020 0.000 0.000 K3 -0.021 0.000 0.000 Standard deviation(mm) K1 0.024 0.025 0.036 K2 0.023 0.024 0.042 K3 0.023 0.024 0.035Table 6
Output results of EX0, EX1, and EX2 for Case2.
Output result  EX0 EX1 EX2 probability of assembly feasibility (Pa)  82.88% probability of quality feasibility (Pq)  82.61% 82.87% 82.38% DFP Mean(mm) X 0.000 0.000 0.000 Y -0.040 0.000 0.000 Standard deviation(mm) X 0.026 0.021 0.025 Y 0.026 0.021 0.058 DFP_QR Mean(mm) X 0.000 0.000 0.000 Y -0.040 0.000 0.000 Standard deviation(mm) X 0.026 0.021 0.025 Y 0.026 0.021 0.058Figure 18
Distributions of assembly deviations for Case1.
(a)
EX0 for Case1 (b)
EX1 for Case1 (c)
EX2 for Case1Figure 19
Distributions of assembly deviations for Case2.
(a)
EX0 for Case2 (b)
EX1 for Case2 (c)
EX2 for Case2Figure 20
Distributions ofτ2 for Case 1.
(a)
EX0 for Case1 (b)
EX1 for Case1 (c)
EX2 for Case1Figure 21
Distributions ofτ2 for Case 2.
(a)
EX0 for Case2 (b)
EX1 for Case2 (c)
EX2 for Case2From Tables5 and 6 and Figures 18–21, the following is observed:( 1 ) As to EX0, except for the output results alongX-axis of Case 2, the mean of DFP is unequal to 0. The DFP_NQR and the corresponding τ2 are unsymmetrically distributed. This may be because of the following: (a) part 2 (PCS2) drops along one constant direction to the boundary of feasible solution spaceDfa_in, where the constant direction is related to Yg-axis; (b)Yfj-axis, j=1,2,…,Nf, is not perpendicular to Yg-axis.( 2 ) As to EX1,Pq is the highest. The DFP_NQR and the corresponding τ2 are symmetrically distributed. This may be due to the normal distribution of τ2 withinDfa_in. The deviation τ2 in the center ofDfa_in has much more probabilities to be selected than the one around the boundary ofDfa_in. The former one is more likely to correspond to the smaller DFP than the latter one.( 3 ) As to EX2, thePq is the lowest. For the DFP corresponding to Yfjfj′(τ2) in (23), the standard deviation is markedly lager than that of EX0 and EX1. These DFP are bimodal distributed, and the DFP_NQR are symmetrically distributed. Furthermore, the point cloud of τ2 has two areas of high-density points, and τ2 corresponding to the DFP_NQR are also symmetrically distributed. This may be because the optimization function in (23) has the form “max⁡·”. Part 2 may fall along two opposite directions to the boundary ofDfa_in, where the opposite directions are related toXfj-axis orYfj-axis, j=1,2,…,Nf.Certainly the output results of EX0 are significantly different from that of EX1 and EX2. For the feasible assembly of over-constrained mechanical structure, gravity has a remarkable impact on gap configuration and then assembly deviations. The gravity influence in tolerance analysis of over-constrained assembly cannot be ignored.


## 5.1. Case Description

Case 1.
This case is from a fixture with boards as shown in Figure10, which is widely used for locating aircraft panels in the panel assembly process. The small working surfaces on a fixture board are applied to constraining the normal displacement of the panels. Each fixture board is mounted on fixture frame with two planar hole-pin-hole joints. This case is the assembly of a fixture board and a fixture frame, denoted as Part 1 and Part 2. Each small working surface of Part 2 is abstracted as a functional point.Figure 10
An assembly fixture with boards.Figure11 shows the definition of CSs. For simplicity, the axes of GCS, GFCS, and nominal PCS2 are assumed to be all parallel to each other. As to FPCSs, three of the functional points are selected for discussion. TheYfj-axis (j=1,2,3) of a FPCS is along the nominal normal direction of the corresponding working surface. The dimensions and tolerances are shown in Figure 12. The fit tolerances of hole-pin-hole pairs are both H7/g7/H7. The coordinates of related points w.r.t. GCS/PCS2 are as B1G=195246.126T, B2G=2805246.126T, oG=1500200T, and C12=-130546.126T, C22=130546.126T, G2=0524.435T, Ok2=0-1200T. Besides, the quality requirement is QR=QR1QR2QR3, where QRj=QRjYfj={YQRj∣YQRj∈[-0.07,0.07]},j=1,2,3.Figure 11
CSs of a fixture frame and a fixture board.Figure 12
Dimensions and tolerances of a fixture frame and a fixture board.Case 2.
This case is the assembly of two aircraft fuselage sections from a type of military aircraft. The front section and the rear section, denoted as Part 1 and Part 2, are assembled with the constraints of eighteen hole-pin-hole pairs on planar bulkheads. The quality requirement concerns the displacement deviations of Part 2 w.r.t Part 1.Figure13 shows the CSs, dimensions and tolerances. In order to simplify the computation, it is assumed that GCS, GFCS, FPCS and nominal PCS2 are all coincided with each other. Besides, the fit tolerances of hole-pin-hole pairs are all D8/d7/D8. The quality requirement is QR=QRXfQRYfT, QRXf={XQR∣XQR∈[-0.1,0.1]}, and QRYf={YQR∣YQR∈[-0.1,0.1]}.Figure 13
CSs, Dimensions and tolerances of two fuselage sections.


## 5.2. Stability and Accuracy of Optimization (PSO)

As seen in Section3, PSO is used to determine assembly deviations, which further impacts the output results of statistical tolerance analysis. In this PSO algorithm, the swarm size is 10, and the number of iterations is 20. PSO is a heuristic optimization algorithm, with simple parameters and fast convergence, but it is easily trapped into local optimum. Therefore, it is necessary to verify stability and accuracy of this optimization algorithm.Firstly, the experiment of statistical tolerance analysis is performed and repeated 10 times, and the 10 sets of output results are observed and compared. In order to eliminate the disturbance of different generated deviation sources, the sameNMC sets of deviation sources are required for all experiments.Three iterative sub-processes are separated from the experiment process. The iterative sub-process of deviation sources generation is denoted as Sub-process 1, andNMC sets of generated deviation sources are stored in a matrixHM. Then with these deviation sources, the iterative sub-process of assembly feasibility judgment is denoted as Sub-process 2, andNa (Na≤NMC) sets of deviation sources, which can make sure the feasibility of assembly, are selected and stored in a matrixHa. Finally, with the deviation sources fromHa, the iterative sub-process of assembly deviation calculation and quality feasibility judgment is denoted as Sub-process 3. After the three sub-processes, the output results can be calculated.Therefore, the 10 experiments are designed as: Sub-processes 1 and 2 are implemented in advance, and then Sub-process 3 is repeated 10 times. The 10 sets of output results with 4 figures after the decimal point are all the same, as shown in Table1 for Case 1 and Table 2 for Case 2, where deviations of functional points (denoted as DFP) are classified into the deviations satisfying quality requirement (denoted as DFP_QR) and that not satisfying quality requirement (denoted as DFP_NQR).Table 1
Output results of experiments for Case1.
Output result K1 K2 K3 Probability of quality feasibility (Pq) 80.22% DFP Mean(mm) -0.0208 -0.0235 -0.0243 Standard deviation(mm) 0.0267 0.0255 0.0250 DFP_QR Mean(mm) -0.0177 -0.0204 -0.0213 Standard deviation(mm) 0.0242 0.0231 0.0227Table 2
Output results of experiments for Case2.
Output result X Y Probability of quality feasibility (Pq) 82.62% DFP Mean(mm) -0.0010 -0.0402 Standard deviation(mm) 0.0261 0.0263 DFP_QR Mean(mm) -0.0010 -0.0400 Standard deviation(mm) 0.0262 0.0261For further verification, with a set of deviation sources arbitrarily selected fromHa, the PSO procedure is performed and repeated 10 times, and the 10 sets of optimization results are compared. Tables 3 and 4 show the selected deviation sources for Cases 1 and 2, respectively. The corresponding feasible solution spaceDfa_in is shown in Figure 14 for Case 1 and Figure 15 for Case 2. In the 10 PSO procedures, the global best particles γ2 and the corresponding Yg-coordinates are shown in Figure 16 for Case 1 and Figure 17 for Case 2. It is observed that the 10 optimization procedures for each case are all convergent, and the 10 sets of optimization results, including global best particles and objective function values, are all the same when these values are rounded to 6 decimal places.Table 3
A set of deviation sources for Case1.
Deviation source Coordinates of Bi w.r.t. GCS Coordinates of Ci w.r.t. PCS2 s i (mm) X Y x y Hole B1-pin-Hole C1 194.9787 246.1397 -1305.0518 46.1361 0.0324 Hole B2-pin-Hole C2 2804.9946 246.1302 1305.0297 46.1142 0.0351Table 4
A set of deviation sources for Case2.
Deviation source Coordinates of Bi w.r.t. GCS Coordinates of Ci w.r.t. PCS2 s i (mm) X Y x y Hole B1-pin-Hole C1 450.0372 -0.0106 449.9975 -0.0411 0.1358 Hole B2-pin-Hole C2 422.7972 153.8436 422.8378 153.9293 0.1250 Hole B3-pin-Hole C3 344.7235 289.2612 344.7289 289.2649 0.1273 Hole B4-pin-Hole C4 225.0135 389.7606 225.0130 389.7105 0.1300 Hole B5-pin-Hole C5 78.1725 443.1735 78.1648 443.1337 0.1236 Hole B6-pin-Hole C6 -78.0688 443.2002 -78.1283 443.1470 0.1217 Hole B7-pin-Hole C7 -225.0053 389.6915 -224.9820 389.6781 0.1264 Hole B8-pin-Hole C8 -344.7693 289.2210 -344.6785 289.2782 0.1218 Hole B9-pin-Hole C9 -422.8215 153.9050 -422.8642 153.9120 0.1149 Hole B10-pin-Hole C10 -450.0165 0.0108 -450.0530 0.0023 0.1264 Hole B11-pin-Hole C11 -422.8637 -153.9604 -422.8953 -153.9496 0.1298 Hole B12-pin-Hole C12 -344.7204 -289.2573 -344.6876 -289.2419 0.1211 Hole B13-pin-Hole C13 -224.9730 -389.7477 -224.9676 -389.7089 0.1179 Hole B14-pin-Hole C14 -78.1358 -443.1664 -78.1309 -443.1260 0.1352 Hole B15-pin-Hole C15 78.1410 -443.1581 78.1550 -443.2297 0.1315 Hole B16-pin-Hole C16 224.9878 -389.6686 224.9854 -389.7149 0.1369 Hole B17-pin-Hole C17 344.7059 -289.2523 344.6852 -289.1532 0.1266 Hole B18-pin-Hole C18 422.7669 -153.9629 422.8744 -153.8547 0.1373Figure 14
A feasible solution space for Case1.Figure 15
A feasible solution space for Case2.Figure 16
Global best particles and objective function values of PSO procedures for Case1.
(a)
Global best particles for Case1 (b)
Objective function values for Case1Figure 17
Global best particles and objective function values of PSO procedure for Case2.
(a)
Global best particles for Case2 (b)
Objective function values for Case2Therefore, Tables1 and 2 and Figures 16 and 17 indicate this PSO algorithm is quite robust and accurate, and thus the output results of tolerance analysis by using the proposed method are reliable.


## 5.3. Comparison with Random Configuration and Worst Configuration of Gaps

The proposed method considers the influence of gravity on the gap configuration when determining the assembly deviations. In order to exhibit the significance of incorporating the influence of gravity to tolerance analysis, two comparison experiments with common gap configurations are given. The first comparison experiment is denoted as EX1, where the deviationτ2 is randomly generated within the feasible solution space based on normal distribution. The second comparison experiment, denoted as EX2, uses the worst configuration of gaps, which makes the assembly deviations largest. Besides, the experiment considering the gravity influence is denoted as EX0.In EX 2, the determination of the worst configuration of gaps is an optimization problem. The modeling and solving of this optimization are both similar to (16)-(19). For example, worst gap configuration corresponds to the largest assembly deviation of thejth functional point alongYfj-axis. The optimization is modeled as seen in (23), and solved by PSO with a nested LP. (23)maxFfd_worτ2=Yfjfj′τ2-YQRj_min+YQRj_max2s.t.τ2∈Dfa_inIn EX2 of Case1, it is assumed that the worst configuration of gaps corresponds to the largest assembly deviations of K2 alongYf2-axis. In EX2 of Case 2, the worst configuration of gaps is assumed to concern the largest assembly deviations alongY-axis, to make a straight comparison between the output results of EX2 and EX0.The output results of EX0, EX1 and EX2 are compared. In order to avoid the distractions of different generated deviation sources, the three experiments, denoted as a group of experiments, are designed as follows: Sub-processes 1 and 2 are implemented in advance, and then Sub-process 3 for EX0, EX1, and EX2 is performed, respectively. To enhance the reliability of the predicted output results, a group of experiments are repeated 10 times, and the 10 groups of output results are averaged, which are shown in Table5 for Case 1 and Table 6 for Case 2. Furthermore, by taking the output results of a group of experiments as examples, the distributions of assembly deviations are shown in Figure 18 for Case 1 and Figure 19 for Case 2, and the distributions of the corresponding τ2 are shown in Figure 20 for Case 1 and Figure 21 for Case 2.Table 5
Output results of EX0, EX1, and EX2 for Case1.
Output result  EX0 EX1 EX2 Probability of assembly feasibility (Pa)  87.01% Probability of quality feasibility (Pq)  80.24% 85.26% 73.55% DFP Mean(mm) K1 -0.021 0.000 0.000 K2 -0.023 0.000 0.000 K3 -0.024 0.000 0.000 Standard deviation(mm) K1 0.026 0.026 0.042 K2 0.026 0.025 0.047 K3 0.025 0.025 0.041 DFP_QR Mean(mm) K1 -0.018 0.000 0.000 K2 -0.020 0.000 0.000 K3 -0.021 0.000 0.000 Standard deviation(mm) K1 0.024 0.025 0.036 K2 0.023 0.024 0.042 K3 0.023 0.024 0.035Table 6
Output results of EX0, EX1, and EX2 for Case2.
Output result  EX0 EX1 EX2 probability of assembly feasibility (Pa)  82.88% probability of quality feasibility (Pq)  82.61% 82.87% 82.38% DFP Mean(mm) X 0.000 0.000 0.000 Y -0.040 0.000 0.000 Standard deviation(mm) X 0.026 0.021 0.025 Y 0.026 0.021 0.058 DFP_QR Mean(mm) X 0.000 0.000 0.000 Y -0.040 0.000 0.000 Standard deviation(mm) X 0.026 0.021 0.025 Y 0.026 0.021 0.058Figure 18
Distributions of assembly deviations for Case1.
(a)
EX0 for Case1 (b)
EX1 for Case1 (c)
EX2 for Case1Figure 19
Distributions of assembly deviations for Case2.
(a)
EX0 for Case2 (b)
EX1 for Case2 (c)
EX2 for Case2Figure 20
Distributions ofτ2 for Case 1.
(a)
EX0 for Case1 (b)
EX1 for Case1 (c)
EX2 for Case1Figure 21
Distributions ofτ2 for Case 2.
(a)
EX0 for Case2 (b)
EX1 for Case2 (c)
EX2 for Case2From Tables5 and 6 and Figures 18–21, the following is observed:( 1 ) As to EX0, except for the output results alongX-axis of Case 2, the mean of DFP is unequal to 0. The DFP_NQR and the corresponding τ2 are unsymmetrically distributed. This may be because of the following: (a) part 2 (PCS2) drops along one constant direction to the boundary of feasible solution spaceDfa_in, where the constant direction is related to Yg-axis; (b)Yfj-axis, j=1,2,…,Nf, is not perpendicular to Yg-axis.( 2 ) As to EX1,Pq is the highest. The DFP_NQR and the corresponding τ2 are symmetrically distributed. This may be due to the normal distribution of τ2 withinDfa_in. The deviation τ2 in the center ofDfa_in has much more probabilities to be selected than the one around the boundary ofDfa_in. The former one is more likely to correspond to the smaller DFP than the latter one.( 3 ) As to EX2, thePq is the lowest. For the DFP corresponding to Yfjfj′(τ2) in (23), the standard deviation is markedly lager than that of EX0 and EX1. These DFP are bimodal distributed, and the DFP_NQR are symmetrically distributed. Furthermore, the point cloud of τ2 has two areas of high-density points, and τ2 corresponding to the DFP_NQR are also symmetrically distributed. This may be because the optimization function in (23) has the form “max⁡·”. Part 2 may fall along two opposite directions to the boundary ofDfa_in, where the opposite directions are related toXfj-axis orYfj-axis, j=1,2,…,Nf.Certainly the output results of EX0 are significantly different from that of EX1 and EX2. For the feasible assembly of over-constrained mechanical structure, gravity has a remarkable impact on gap configuration and then assembly deviations. The gravity influence in tolerance analysis of over-constrained assembly cannot be ignored.


## 6. Conclusions

This paper presents a tolerance analysis method for over-constrained assembly considering the influence of gravity. The proposed method contains deviation propagation analysis and statistical tolerance analysis. The probabilities of assembly feasibility and quality feasibility can be computed, and the deviations distributions for feasible assemblies can be quantified.The detailed tolerance analysis process is illustrated for a typical over-constrained mechanical structure with constraints of multiple planar hole-pin-hole pairs. This type of constraints has non-linearity, which increases complexity of assembly tolerance analysis. In the deviation propagation analysis, the judgment of assembly feasibility and the calculation of assembly deviations are formulated as a search problem and an optimization problem, respectively. Due to the non-linear constraints, it is inefficient and difficult to solve these problems. In order to improve accuracy and efficiency, according to characteristic of this type of constraints, the search and optimization problems are solved by reducing object dimensionality and narrowing search space, where the search algorithms, such as sequential search and binary search, as well as the optimization algorithms, such as PSO and LP, are used.The results of the case studies show the following: (a) the PSO algorithm for the assembly deviation determination is demonstrated to be quite robust and accurate, thus the output results of the tolerance analysis are reliable and accurate enough; (b) for predicting the assembly deviations of the over-constrained mechanical structure, the influence of gravity is remarkable and cannot be ignored.Our main contribution lies in developing a reliable and accurate method for assembly tolerance analysis with overconstraints of multiple planar hole-pin-hole pairs, considering the influence of gravity. Besides gravity, it can be used to incorporate the influence of other conservative forces into tolerance analysis of over-constrained assembly as well.In the proposed method, Monte Carlo simulation (MCS) is used for the statistical tolerance analysis, due to its general applicability for implicit and non-linear problem. The accurate and reliable application of the MCS lightly depends on the number of iterations. The accuracy of prediction can be improved by increasing the number of iterations at the cost of computational time. To relieve the required computational effort, the first-order reliability method (FORM) method combining with a reasonable linearization strategy for the non-linear constraints will enable a reasonably accurate prediction of defect probabilities at a lower computational effort. However, it should be noted that the loss of accuracy of prediction on defect probability would be unavoidable.Furthermore, the proposed method is applied to the two-dimensional tolerance analysis of over-constrained assemblies, but the analytical framework of the proposed method may be extended to three dimension applications, i.e., for over-constrained mechanical structure, the assembly feasibility judgment and the assembly deviation determination considering the influence of gravity can be modeled as a search problem and an optimization problem, respectively. For the three-dimensional case, the DOFs in space instead of the ones on plane are required to be considered, and the analysis of assembly requirement is more challenging. Then, the search and optimization problems for judging assembly feasibility and determining assembly deviations are more difficult to be modeled and solved. The authors are already working on the extension of the method employed in this paper to the three-dimensional tolerance analysis of over-constrained assembly. Besides, in future work, part deformation for eliminating the interference between different parts will also be integrated into the assembly tolerance analysis.



---

*Source: 2039153-2018-11-14.xml*