A novel tolerance analysis method for three-dimensional assembly

Abstract

Three-dimensional tolerance analysis is increasingly becoming an innovative method for computer-aided tolerancing. Its aim is to support the design, manufacturing, and inspection by providing a quantitative analysis of the effects of multi-tolerances on final functional key characteristics and predict the quality level. This article proposes a novel approach for three-dimensional assembly analysis—a hybridization of vector loop and quasi-Monte Carlo method. The former is used to establish the three-dimensional assembly chain and obtain the assembly function. The latter is adopted to generate n sets of dimensional values according to the distribution of each dimension in chain. The new method is shown to inherit many of the best features of classical vector loop and quasi-Monte Carlo, combining easy-to-obtain assembly function with accurate statistical analysis. For every set of dimensional values, one sample value of a functional requirement can be computed with the Newton–Raphson iterative procedure. A crank slider mechanical assembly is shown as an example to illustrate the proposed method.

Keywords

Three-dimensional tolerance analysis assembly vector loop quasi-Monte Carlo Newton–Raphson iterative

Introduction

Three-dimensional (3D) tolerance analysis is hot focused for ever-tightening and increasing complex requirements in various fields. It enables the prediction of the cumulative effects of these part deviations on assembly and functional key characteristics. Tolerance analysis aims at not only improving the performance of a product but also reducing the cost.

A large amount of research has been devoted to two aspects in tolerancing domain: tolerance model and analysis method. Some models of 3D tolerance analysis are reviewed in detail and compared based on the literature published.^1,2 A kinematic tolerance space model was proposed to express the relationship between the nominal kinematics of mechanisms and their kinematic variations.³ The T-Map represents all possible variations in size, position, form, and orientation for a target feature by means of a hypothetical Euclidean point-space.⁴ Mansuy et al.⁵ applied deviation domain to perform both worst case and statistical tolerancing analysis. Similar to T-Map and deviation domain, the polytope represents geometrical deviations on mechanical parts by sets of half-spaces coming from the restriction of the position of a finite set of points of a toleranced feature.⁶ Asante⁷ used the small displacement torsor (SDT) model to analyze the combined effects of workpiece geometric error, locator geometric error, and clamping error on a machined feature. A proposed unified Jacobian–Torsor model combines the benefits for tolerance representation and propagation.⁸ The skin model shape concept,^9–11 developed from skin model theory,^12,13 can be used for assembly simulation¹⁴ and served as digital twin.¹⁵ A gap-based approach, called GapSpace¹⁶ that uses gap to describe the relationship between features belonging to different components, and three main sources of variations are included into vector-loop-based model for a mechanical assembly.¹⁷ Some models have been integrated successfully into CAD packages for industrial environment, such as CATIA.3D FDT and FT&A are based on TTRS model and matrix model; CETOL 6 Sigma and Sigmund use vector-loop model; eM-TolMate, VSA, and 3DCS incorporate variational models; MECAmaster includes torsor model; and PolitoCAT and Politopix are developed based on polytope model. Currently, most of the CAD packages can deal with assemblies of varied complexities, among which CETOL 6 Sigma and Sigmund are able to consider part/assembly variation and kinematic motion for kinematic mechanisms due to the adoption of vector loop model.

The tolerance analysis methods are mainly worst case or statistical case. Statistical analysis is a more demanding way by permitting a small fraction of assemblies to not assemble or function as required, an increase in tolerances for individual dimensions may be obtained, and in turn, the manufacturing costs may be significantly reduced.¹⁸ There exists root sum squared (RSS) method,¹⁹ system moments,²⁰ Taguchi’s method,²¹ the iterative method,²² product of exponentials,²³ and Monte Carlo simulations^24–27 for statistical tolerance analysis. In particular, the Monte Carlo method is the most widely used method in statistical tolerance analysis as it allows the use of normal and non-normal distributions, as well as linear and nonlinear assembly functions. The CAD packages, such as CETOL 6 Sigma, VSA, eM-TolMate, 3DCS, Mechanical Advantage, Sigmund, and RD&T, are able to use Monte Carlo to perform simulations. However, this method requires a large number of samples to assure accuracy. Motivated by this shortcoming, a so-called good point set of quasi-Monte Carlo method was used to generate samples of the dimensions.²⁸ This method resembles Monte Carlo in applications that are suitable for either normal or non-normal input/output, and linear or nonlinear assembly functions. The difference between these two methods lies in the fact that the quasi-random numbers generated for quasi-Monte Carlo have a lower discrepancy than that generated by the Monte Carlo method, and therefore fewer samples are necessary to guarantee accuracy. Previous work presented for quasi-Monte Carlo focused on mechanical assembly with explicit objective function.²⁸ To promote the applicability of the method in tolerancing domain, the vector loop is introduced, which is an innovative category of tolerance analysis method that three main sources of deviations in a mechanical assembly are fully considered.^29,30 The direct linearization method (DLM), the second-order tolerance analysis (SOTA), and the Monte Carlo method have been developed to derive the functional requirement in vector loop model. Nevertheless, DLM may be inaccurate for highly nonlinear assemblies, and SOTA is an approximation due to the approximation of the three sets of partial derivatives in every iteration while Monte Carlo requires large computations. In this regard, an approach is proposed for 3D assembly tolerance analysis which seeks to combine the benefits of the classical vector loop and quasi-Monte Carlo. The new approach is designed to describe the implicit assembly function consistent with vector loop scheme, while also inheriting the sample strength of quasi-Monte Carlo. This novel method is an extended work of quasi-Monte Carlo in tolerance analysis.

The rest of the article is organized as follows. Section “Vector-loop-based assembly function” illustrates the representation of the assembly function. Section “Newton–Raphson iterative method” proposes a solution toward to solving the assembly function. Section “Dimension sample method” describes the quasi-Monte Carlo method for generating samples. Section “3D tolerance analysis method” presents the overall 3D tolerance analysis process. A 3D mechanism is presented to verify the validity and accuracy of the proposed method in section “Case study.” Finally, section “Conclusion” provides the conclusions.

Vector-loop-based assembly function

An assembly dimension chain can be described as a closed vector loop. In two-dimensional (2D) space, the vector loop equation results in three scalar equations, representing the sum of the x and y projections of the dimension vectors in global coordinates and the sum of the relative rotations. However, in 3D space, the projection equations still hold in global x, y, and z, but the rotation constraints can only be expressed in matrix form. A vector loop that contains m vectors $L_{1}, L_{2}, \dots, L_{m}$ is shown in Figure 1. The origin of the $x_{0} y_{0} z_{0}$ system is the starting point of the first vector $L_{1}$ and its orientation is the same as the global system. The origin of the system $x_{i} y_{i} z_{i}, i = 1, \dots, m,$ is the end point of the vector $L_{i}$ and the orientation of the x-axis is the same as that of the $L_{i}$ vector.

Figure 1.

Vector loop and local coordinate systems.

A closed 3D loop can be expressed in terms of a concatenation of coordinate transformation matrices as shown in equation (1). Equation (1) describes that the coordinate system at the end of the loop must be truly the same as the coordinate system at the beginning

R_{1} \cdot T_{1} \cdot R_{2} \cdot T_{2} \cdot \dots \cdot R_{i} \cdot T_{i} \cdot \dots \cdot R_{m} \cdot T_{m} \cdot R_{f} = I

(1)

where $R_{i}$ is the product of rotation matrices, $T_{i}$ is the translation matrix, $R_{f}$ is the final rotation transformation matrix for bringing the loop to closure, and I is the identity matrix.

To describe the loop equation in matrix form, four elementary homogeneous transform matrices, as shown schematically in Figure 2, are used.

The $x_{2} y_{2} z_{2}$ system is obtained from $x_{1} y_{1} z_{1}$ by rotation about the origin by the angle $θ$ around ${[1, 0, 0]}^{T}$ , as shown in Figure 2(a). The corresponding rotation matrix is

R_{x} (θ) = [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \cos θ & - \sin θ & 0 \\ 0 & \sin θ & \cos θ & 0 \\ 0 & 0 & 0 & 1 \end{matrix}]

(2)

The $x_{2} y_{2} z_{2}$ system is obtained from $x_{1} y_{1} z_{1}$ by rotation about the origin by the angle $θ$ around ${[0, 1, 0]}^{T}$ as shown in Figure 2(b). The corresponding rotation matrix is

R_{y} (θ) = [\begin{matrix} \cos θ & 0 & \sin θ & 0 \\ 0 & 1 & 0 & 0 \\ - \sin θ & 0 & \cos θ & 0 \\ 0 & 0 & 0 & 1 \end{matrix}]

(3)

The $x_{2} y_{2} z_{2}$ system is obtained from $x_{1} y_{1} z_{1}$ by rotation about the origin by the angle $θ$ around ${[0, 0, 1]}^{T}$ as shown in Figure 2(c). The corresponding rotation matrix is

R_{z} (θ) = [\begin{matrix} \cos θ & - \sin θ & 0 & 0 \\ \sin θ & \cos θ & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}]

(4)

The $x_{2} y_{2} z_{2}$ system is obtained from $x_{1} y_{1} z_{1}$ by translation ${[l, 0, 0]}^{T}$ only, as shown in Figure 2(d). The corresponding translation matrix is

T (l) = [\begin{matrix} 1 & 0 & 0 & l \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}]

(5)

Figure 2.

Three-dimensional representations for the four elementary matrices (a) Rotation around the x axis, (b) rotation around the y axis, (c) rotation around the z axis and (d) translation along the x axis.

With this convention, the vector loop equation (equation (1)) can be expressed in matrix form, therefore, forming the necessary assembly function. This method is well capable of obtaining assembly function on occasions that are difficult or impossible to calculate explicit function.

Newton–Raphson iterative method

For a general 3D assembly, the functional requirement is expressed as a nonlinear function of the dimensions and constants of the assembly, and it is difficult to solve this equation for a general case. The nonlinear equation requires a nonlinear solver to find solution. Since Newton–Raphson method has the best convergence rate of nonlinear system solvers, in this section, the method is adopted to solve the equation. In this regard, the iterative formula is mostly required. Before the iterative formula is deduced, derivatives of the elementary matrices and partial derivatives of function $S (q)$ (see equation (16)) are introduced, respectively, in sections “Derivatives of the elementary matrices” and “Partial derivatives of function $S (q)$ ”, since the operation of analytic derivatives of the loop equations can reduce the computation of Newton–Raphson iterative.

Derivatives of the elementary matrices

It is mathematically convenient to describe the operation of taking the derivative as a matrix product. Define the derivatives of the elementary matrices as follows

\frac{d R_{x} (θ)}{d θ} = R_{x} (θ) Q_{1}

(6)

\frac{d R_{y} (θ)}{d θ} = R_{y} (θ) Q_{2}

(7)

\frac{d R_{z} (θ)}{d θ} = R_{z} (θ) Q_{3}

(8)

\frac{dT (l)}{dl} = T (l) Q_{4}

(9)

where $Q_{1}$ , $Q_{2}$ , $Q_{3}$ , and $Q_{4}$ are the four derivative operators associated with the type of elementary matrices. The calculation of equations (6)–(9) derives

Q_{1} = [\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 0 & - 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}]

(10)

Q_{2} = [\begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}]

(11)

Q_{3} = [\begin{matrix} 0 & - 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}]

(12)

Q_{4} = [\begin{matrix} 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}]

(13)

Partial derivatives of function S(q)

For convenience, equation (1) is rewritten as equation (14). It is important to note that $S_{(j - 1) j}, j = 1, 2, \dots, k$ , and final $S_{k 0}$ are the elementary translation and rotation matrices in equation (1)

S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} \cdot S_{j (j + 1)} \cdot \dots \cdot S_{(k - 1) k} \cdot S_{k 0} = I

(14)

For a 3D assembly, there exist six unknown dimensional variables or assembly variables. Define the left-hand side of equation (14) as a function depicted in equation (15), then the derivative of which can be represented as equation (16). Because of the existence of equations (17) and (18), equation (16) changes into equation (19).

S (q) = S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} \cdot S_{j (j + 1)} \cdot \dots \cdot S_{(k - 1) k} \cdot S_{k 0}

(15)

where $q = (q_{1}, q_{2}, q_{3}, q_{4}, q_{5}, q_{6})$ and $q_{i}$ represents the unknown dimensional or assembly variables

\frac{\partial S (q)}{\partial q_{i}} = S_{01} \cdot S_{12} \cdot \dots \cdot \frac{d S_{(j - 1) j} (q_{i})}{d q_{i}} \cdot S_{j (j + 1)} \cdot \dots \cdot S_{(k - 1) k} \cdot S_{k 0}

(16)

\frac{d S_{(j - 1) j} (q_{i})}{d q_{i}} = S_{(j - 1) j} (q_{i}) \cdot Q

(17)

where Q is the proper type of derivative operator matrix for the elementary matrix $S_{(j - 1) j} (q_{i})$

\begin{array}{l} I = {(S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i}))}^{- 1} \\ \cdot (S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i})) \end{array}

(18)

\begin{matrix} \frac{\partial S (q)}{\partial q_{i}} = S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i}) \\ \cdot Q \cdot (S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i}))^{- 1} \\ \cdot (S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i})) \\ \cdot S_{j (j + 1)} \cdot \dots \cdot S_{(k - 1) k} \cdot S_{k 0} \\ = S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i}) \cdot Q \\ \cdot (S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i}))^{- 1} \cdot S (q) \end{matrix}

(19)

For simplifying calculation purpose, we assume

\begin{array}{l} D_{i} = S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i}) \\ \cdot Q \cdot {(S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i}))}^{- 1} \end{array}

(20)

where $D_{i}$ represents the derivative operator matrix of $S (q)$ with respect to the ith variable in the loop. According to the theory of homogeneous transformation, equations (21) and (22) can be obtained. Therefore, different values of $D_{i}$ can be obtained based on four different derivative operators Q

S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i}) = [\begin{matrix} r_{11} & r_{12} & r_{13} & Δ x \\ r_{21} & r_{22} & r_{23} & Δ y \\ r_{31} & r_{32} & r_{33} & Δ z \\ 0 & 0 & 0 & 1 \end{matrix}]

(21)

\begin{array}{l} {(S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i}))}^{- 1} \\ = [\begin{matrix} r_{11} & r_{21} & r_{31} & - (Δ x \cdot r_{11} + Δ y \cdot r_{21} + Δ z \cdot r_{31}) \\ r_{12} & r_{22} & r_{32} & - (Δ x \cdot r_{12} + Δ y \cdot r_{22} + Δ z \cdot r_{32}) \\ r_{13} & r_{23} & r_{33} & - (Δ x \cdot r_{13} + Δ y \cdot r_{23} + Δ z \cdot r_{33}) \\ 0 & 0 & 0 & 1 \end{matrix}] \end{array}

(22)

where $r_{ij}$ is the element of matrix $S_{01} \cdot S_{12} \cdot \dots \cdot S_{(j - 1) j} (q_{i})$ , row i, and column j. $Δ x$ , $Δ y$ , and $Δ z$ are the elements on rows 1, 2, and 3 of column 4, respectively

If $Q = Q_{1}$

D_{i} = [\begin{matrix} 0 & - r_{31} & r_{21} & Δ y \cdot r_{31} - Δ z \cdot r_{21} \\ r_{31} & 0 & - r_{11} & Δ z \cdot r_{11} - Δ x \cdot r_{31} \\ - r_{21} & r_{11} & 0 & Δ x \cdot r_{21} - Δ y \cdot r_{11} \\ 0 & 0 & 0 & 0 \end{matrix}]

(23)

If $Q = Q_{2}$

D_{i} = [\begin{matrix} 0 & - r_{32} & r_{22} & Δ y \cdot r_{32} - Δ z \cdot r_{22} \\ r_{32} & 0 & - r_{12} & Δ z \cdot r_{12} - Δ x \cdot r_{32} \\ - r_{22} & r_{12} & 0 & Δ x \cdot r_{22} - Δ y \cdot r_{12} \\ 0 & 0 & 0 & 0 \end{matrix}]

(24)

If $Q = Q_{3}$

D_{i} = [\begin{matrix} 0 & - r_{33} & r_{23} & Δ y \cdot r_{33} - Δ z \cdot r_{23} \\ r_{33} & 0 & - r_{13} & Δ z \cdot r_{13} - Δ x \cdot r_{33} \\ - r_{23} & r_{13} & 0 & Δ x \cdot r_{23} - Δ y \cdot r_{13} \\ 0 & 0 & 0 & 0 \end{matrix}]

(25)

If $Q = Q_{4}$

D_{i} = [\begin{matrix} 0 & 0 & 0 & r_{11} \\ 0 & 0 & 0 & r_{21} \\ 0 & 0 & 0 & r_{31} \\ 0 & 0 & 0 & 0 \end{matrix}]

(26)

Iterative formula

According to equations (14) and (15), the assembly function can be rewritten as

S (q) - I = 0

(27)

To solve the assembly function in equation (27), a Newton–Raphson iterative procedure is carried out. First, an initial value $q_{0}$ , normally nominal dimensions for closing link or assembly variables are given. Then, a small variation $Δ q_{0}$ is added to this initial value to form a new initial value. The iterative procedure stops at a specified accuracy.

The first-order linear expansion of equation (27) about $q_{0}$ derives

S (q_{0}) - I + \sum_{i = 1}^{6} \frac{\partial S (q_{0})}{\partial q_{i}} Δ q_{i} = 0

(28)

As per equations (17)–(20), equation (28) becomes

S (q_{0}) + \sum_{i = 1}^{6} (D_{i} S (q_{0}) Δ q_{i}) - I = 0

(29)

Post-multiplying by $S^{- 1} (q_{0})$ yields

\sum_{i = 1}^{6} (D_{i} Δ q_{i}) = S^{- 1} (q_{0}) - I

(30)

According to equations (23)–(26), a general form of $D_{i}$ can be approached

D_{i} = [\begin{matrix} 0 & - d_{i_{21}} & d_{i_{13}} & d_{i_{14}} \\ d_{i_{21}} & 0 & - d_{i_{32}} & d_{i_{24}} \\ - d_{i_{13}} & d_{i_{32}} & 0 & d_{i_{34}} \\ 0 & 0 & 0 & 0 \end{matrix}]

(31)

It is obvious that only ${d_{i}}_{13}$ , ${d_{i}}_{21}$ , ${d_{i}}_{32}$ , ${d_{i}}_{14}$ , ${d_{i}}_{24}$ , and ${d_{i}}_{34}$ can be used to formulate independent equations.

For the right-hand side of equation (30), $S (q_{0})$ can be defined as equation (32), and then equation (33) is obtained

S (q_{0}) = [\begin{matrix} c_{11} & c_{12} & c_{13} & Δ u \\ c_{21} & r_{22} & c_{23} & Δ v \\ c_{31} & c_{32} & c_{33} & Δ w \\ 0 & 0 & 0 & 1 \end{matrix}]

(32)

S^{- 1} (q_{0}) - I = [\begin{matrix} c_{11} - 1 & c_{21} & c_{31} & - (Δ u \cdot c_{11} + Δ v \cdot c_{21} + Δ w \cdot c_{31}) \\ c_{12} & c_{22} - 1 & c_{32} & - (Δ u \cdot c_{12} + Δ v \cdot c_{22} + Δ w \cdot c_{32}) \\ c_{13} & c_{23} & c_{33} - 1 & - (Δ u \cdot c_{13} + Δ v \cdot c_{23} + Δ w \cdot c_{33}) \\ 0 & 0 & 0 & 0 \end{matrix}]

(33)

where $c_{ij}$ is the element of matrix $S (q_{0})$ , row i, and column j. $Δ u$ , $Δ v$ , and $Δ w$ are the elements on rows 1, 2, and 3 of column 4, respectively.

In summary, the scalar equations for the iteration procedure can be written in matrix form as expressed in equation (34)

[\begin{matrix} \begin{matrix} d_{1_{13}} & d_{2_{13}} \\ d_{1_{21}} & d_{2_{13}} \end{matrix} & \begin{matrix} d_{3_{13}} & d_{4_{13}} \\ d_{3_{13}} & d_{4_{13}} \end{matrix} & \begin{matrix} d_{5_{13}} & d_{6_{13}} \\ d_{5_{13}} & d_{6_{13}} \end{matrix} \\ \begin{matrix} d_{1_{32}} & d_{2_{13}} \\ d_{1_{14}} & d_{2_{13}} \end{matrix} & \begin{matrix} d_{3_{13}} & d_{4_{13}} \\ d_{3_{13}} & d_{4_{13}} \end{matrix} & \begin{matrix} d_{5_{13}} & d_{6_{13}} \\ d_{5_{13}} & d_{6_{13}} \end{matrix} \\ \begin{matrix} d_{1_{24}} & d_{2_{13}} \\ d_{1_{34}} & d_{2_{13}} \end{matrix} & \begin{matrix} d_{3_{13}} & d_{4_{13}} \\ d_{3_{13}} & d_{4_{13}} \end{matrix} & \begin{matrix} d_{5_{13}} & d_{6_{13}} \\ d_{5_{13}} & d_{6_{13}} \end{matrix} \end{matrix}] [\begin{matrix} \begin{matrix} Δ q_{1} \\ Δ q_{2} \end{matrix} \\ \begin{matrix} Δ q_{3} \\ Δ q_{4} \end{matrix} \\ \begin{matrix} Δ q_{5} \\ Δ q_{6} \end{matrix} \end{matrix}] = [\begin{matrix} \begin{matrix} c_{31} + c_{11} + c_{33} - 2 \\ c_{12} + c_{11} + c_{22} - 2 \end{matrix} \\ \begin{matrix} c_{23} + c_{33} + c_{22} - 2 \\ - (Δ u \cdot c_{11} + Δ v \cdot c_{21} + Δ w \cdot c_{31}) \end{matrix} \\ \begin{matrix} - (Δ u \cdot c_{12} + Δ v \cdot c_{22} + Δ w \cdot c_{32}) \\ - (Δ u \cdot c_{13} + Δ v \cdot c_{23} + Δ w \cdot c_{33}) \end{matrix} \end{matrix}]

(34)

Therefore, the small variation can be obtained as follows

Δ q_{0} = G^{- 1} \cdot H

(35)

where

G = [\begin{matrix} \begin{matrix} d_{1_{13}} & d_{2_{13}} \\ d_{1_{21}} & d_{2_{13}} \end{matrix} & \begin{matrix} d_{3_{13}} & d_{4_{13}} \\ d_{3_{13}} & d_{4_{13}} \end{matrix} & \begin{matrix} d_{5_{13}} & d_{6_{13}} \\ d_{5_{13}} & d_{6_{13}} \end{matrix} \\ \begin{matrix} d_{1_{32}} & d_{2_{13}} \\ d_{1_{14}} & d_{2_{13}} \end{matrix} & \begin{matrix} d_{3_{13}} & d_{4_{13}} \\ d_{3_{13}} & d_{4_{13}} \end{matrix} & \begin{matrix} d_{5_{13}} & d_{6_{13}} \\ d_{5_{13}} & d_{6_{13}} \end{matrix} \\ \begin{matrix} d_{1_{24}} & d_{2_{13}} \\ d_{1_{34}} & d_{2_{13}} \end{matrix} & \begin{matrix} d_{3_{13}} & d_{4_{13}} \\ d_{3_{13}} & d_{4_{13}} \end{matrix} & \begin{matrix} d_{5_{13}} & d_{6_{13}} \\ d_{5_{13}} & d_{6_{13}} \end{matrix} \end{matrix}]

(36)

H = [\begin{matrix} \begin{matrix} c_{31} + c_{11} + c_{33} - 2 \\ c_{12} + c_{11} + c_{22} - 2 \end{matrix} \\ \begin{matrix} c_{23} + c_{33} + c_{22} - 2 \\ - (Δ u \cdot c_{11} + Δ v \cdot c_{21} + Δ w \cdot c_{31}) \end{matrix} \\ \begin{matrix} - (Δ u \cdot c_{12} + Δ v \cdot c_{22} + Δ w \cdot c_{32}) \\ - (Δ u \cdot c_{13} + Δ v \cdot c_{23} + Δ w \cdot c_{33}) \end{matrix} \end{matrix}]

(37)

Dimension sample method

After obtaining the matrices G and H, a quasi-Monte Carlo method is adopted to generate the manufacturing variables.

Take $γ = (γ_{1}, \dots, γ_{n}) \in R^{s}$ . If $P_{n}$ forms the first n points of

{({γ_{1} k}, \dots, {γ_{n} k}), k = 1, 2, \dots}

(38)

where ${γ}$ denotes the fractional part of $γ$ , with discrepancy $D (n) = o (n^{- 1 + ε})$ as $n \to \infty$ , then the set $p_{n}$ is called a good point set (gp set) and $γ$ a good point. The following are some useful good points:

Take $p_{1}, \dots p_{s}$ as the first s primes. We have

γ = (\sqrt{p_{1}}, \dots, \sqrt{p_{s}})

(39)

Take p as a prime and $q = p^{1 / (s + 1)}$ . We have

γ = (q, q^{2}, \dots, q^{s})

(40)

Take p as a prime and $p \geq 2 s + 3$ . We have

γ = ({2 \cos \frac{2 π}{p}}, {2 \cos \frac{4 π}{p}}, \dots, {2 \cos \frac{2 π s}{p}})

(41)

Figure 3(a) and (b) shows the point sets generated by the Monte Carlo method and the quasi-Monte Carlo method. It is obvious that the good point set generated is only decided by s, which represents the number of manufacturing variables. For an assembly, it is easy to obtain the certain variables by creating the vector loop. Therefore, the generated good point set is fixed. It should be mentioned that if $n_{1}$ points have been generated and we need $n_{2} (> n_{1})$ points, $n_{2} - n_{1}$ points are required to be newly generated. That is to say, the previous $n_{1}$ points can be reused. Besides, the low-discrepancy characteristic of the quasi-Monte Carlo algorithm, as shown in Figure 3, enables it to get stable samples in a shorter time than Monte Carlo method. Comparisons between quasi-Monte Carlo and Monte Carlo have been discussed in detail in published research.²⁸

Figure 3.

(a) Points generated by Monte Carlo method and (b) points generated by quasi-Monte Carlo method.

After the good point set is obtained, an inverse transformation method can be used for generating an observation from a given distribution. The natural idea for generating random variables $P = {x_{k}, k = 1, \dots, n}$ for its cumulative distribution function $F (x)$ is as follows:

Generate a good point set ${C_{k} = (C_{k 1}, \dots, C_{ks}), k = 1, \dots, n}$ on $C^{s}$ by means of quasi-Monte Carlo method;

Take $x_{k} = (F_{1}^{- 1} (C_{k 1}), \dots, F_{s}^{- 1} (C_{ks})), k = 1, \dots, n$ , where $F_{i}^{- 1} (x_{i})$ is the inverse function of $F_{i} (x_{i})$ ;

Use $x_{k}$ to compute the samples of function requirement.

3D tolerance analysis method

In this article, a quasi-Monte Carlo method is used to generate n sets of dimensional values according to the distribution of each dimension. By defining the iterative formula, one sample value for the functional requirement can be computed following the Newton–Raphson iterative procedure. Finally, the first four order moments of the function requirement can be computed after n samples of the function requirement are found and thus the assembly reject can be estimated. A flow diagram of the proposed method is shown in Figure 4. The input information includes the probability distributions of the dimensions, the nominal value and tolerance of the functional requirement, and the nominal values of the assembly variables.

Figure 4.

Flow diagram of the proposed method.

After creating the local coordinate systems, the assembly function can be derived, as shown in section “Vector-loop-based assembly function.” $S_{0 j}$ associated with $q_{i}$ is obtained according to equation (19) and then $D_{i}$ can be computed by equations (23)–(26). Afterwards, the matrix G can be yielded based on equation (36). The matrix H can be obtained according to equation (37) and $S (q)$ , at the right-hand side of the assembly function. It should be noted that G and H are the functions of the dimension variables, functional requirement, and assembly variables.

The quasi-Monte Carlo method is presented to generate a sample of the dimensions. The number of dimensions is s. Let p be a prime and $p \geq 2 s + 3$ . A good point $γ$ is found

γ = (γ_{1}, \dots, γ_{s}) = ({2 \cos \frac{2 π}{p}}, {2 \cos \frac{4 π}{p}}, \dots, {2 \cos \frac{2 π s}{p}})

(42)

where ${γ}$ denotes the fractional part of $γ$ .

Then, the points ${c_{k} = ({γ_{1} k}, \dots, {γ_{s} k}), k = 1, \dots, n}$ can be generated and the kth set of dimensions $x_{k}$ is obtained

x_{k} = (F_{1}^{- 1} (c_{k 1}), \dots, F_{s}^{- 1} (c_{ks})), k = 1, \dots, n

(43)

where $F_{i}^{- 1} (x_{i})$ is the inverse function of $F_{i} (x_{i})$ , and $F_{i} (x_{i})$ is the cumulative distribution function of dimension $x_{i}$ .

After generating the sample of dimensions, the Newton–Raphson iterative procedure is used to obtain a sample of the functional requirement. Standard statistical estimation methods are then used to compute the first four order moments of the samples. Finally, the assembly reject is estimated by integration or a table.

Case study

The performance of the proposed method was investigated experimentally, with techniques used to analyze the benchmark mechanism. A 3D crank slider mechanism that consists of a base, a crank, a link, and a slider will be used. The vector loop, local coordinate systems, and the dimensions which govern the operation of the 3D assembly are shown in Figure 5. The assembly reject associated with U due to manufacturing variations is the key issue to determine the assemblability.

Figure 5.

Vector loop and local coordinate systems.

The proposed method consists of the following steps:

Step 1. Identify the variables. The independent manufacturing dimension variables are A, B, C, D, and E, whereas the distance U and angles $θ_{1}$ , $θ_{2}$ , $θ_{3}$ , $θ_{4}$ , and $θ_{5}$ are assembly variables. Table 1 shows the detailed dimensions for the assembly.

Step 2. Obtain the assembly function. The vector-loop-based assembly function can be expressed as

\begin{matrix} R_{y} (- \frac{π}{2}) \cdot T (A) \cdot R_{y} (- \frac{π}{2}) \cdot T (B) \cdot R_{y} (- \frac{π}{2}) \cdot \\ R_{z} (\frac{π}{4}) \cdot T (C) \cdot R_{x} (θ_{1}) \cdot R_{z} (θ_{2}) \cdot T (D) \cdot R_{z} (θ_{3}) \cdot \\ R_{x} (θ_{4}) \cdot R_{z} (θ_{5}) \cdot T (E) \cdot R_{y} (- \frac{π}{2}) \cdot T (U) = I \end{matrix}

(44)

Step 3. Calculate the G and H matrices. The sample of dimensions A, B, C, D, and E can be represented by good point set. Here, $s = 5$ , so set $p = 13$ . According to equation (42), a good point $γ$ can be found, as shown in equation (45)

γ = (0.7709, 0.1361, 0.2411, 0.2908, 0.5030)

(45)

Table 1.

Detailed dimensions for the 3D assembly.

Dimensions	Nominal value	Tolerance (±)
A	500 mm	4 mm
B	300 mm	2.5 mm
C	375 mm	2.5 mm
D	750 mm	4 mm
E	125 mm	2.5 mm
U	992.910 mm	4 mm
$θ_{1}$	69.059°	?
$θ_{2}$	−98.421°	?
$θ_{3}$	−81.579°	?
$θ_{4}$	69.059°	?
$θ_{5}$	135°	?

For every sample of A, B, C, D, and E, a sample of the distance U and angles $θ_{1}$ , $θ_{2}$ , $θ_{3}$ , $θ_{4}$ , and $θ_{5}$ can be obtained through the Newton–Raphson iterative procedure as described in section “Newton–Raphson iterative method.”

Step 4. Compute the mean and standard deviation through statistical formulas, and the assembly reject is estimated by integration.

Step 5. For comparison purpose, the conventional Monte Carlo method in vector loop model was performed, allowing its output to be compared with the proposed method side-by-side. Due to similarities in the overall analysis framework, a common procedure was used, differing only in the generation approaches of random variations.

To provide additional performance information, the result obtained from the explicit assembly function is included in a subset of comparisons. The explicit assembly function is expressed as

U = B + \sqrt{D^{2} - \frac{1}{2} C^{2} - {(A - E - \frac{\sqrt{2}}{2} C)}^{2}}

(46)

The various numbers of sample by Monte Carlo and quasi-Monte Carlo were proceeded. Detailed statistics on analysis results are presented, including mean and standard deviation of the functional requirement. As shown in Figures 6 and 7, overall, the novel method gives a smaller scope of variations with improved sample efficiency to be stable compared to the Monte Carlo method. Focusing on the distribution of the mean, it can be seen that the proposed method has moderate amplitude and becomes stable when n is greater than 500, while Monte Carlo method varies in a more acute sense and stays stable after 600 samples. In terms of standard deviation, the proposed method is seen to achieve significant improvements for keeping relatively small values. In addition, the sample size is smaller when the standard deviation stays steady for the proposed method with respect to the conventional vector loop method. Besides, due to the reusable characteristic of sample strategy in this study, only excessive samples need to be newly generated as sample size increases, which practically reduces the computation efforts.

Figure 6.

Comparison of the means.

Figure 7.

Comparison of the standard deviations.

The comparison between the proposed method and the explicit assembly function method is then performed, as shown in Figures 6 and 7. The results obtained by the proposed method are almost the same as those from the method using an explicit assembly function to obtain the samples of U. Besides, the relative error of mean given by the proposed method is small with a low relative error of the standard deviation, which is less than 1000, indicating the effectiveness and accuracy of solving the implicit assembly function by the presented method. As per Figure 7, when the sample size n reaches 800, the standard deviation of the samples is −0.146%, to a minimum. Simultaneously, the relative error of the mean is as low as 0.001%, in which the comparison object for the relative error is the value of the previous sample. The mean and standard deviation have reached their comparatively stable values. In this regard, the mean and standard deviation associated with the 800 samples size can be used to estimate the assembly reject and the result is 1.924%.

Conclusion

The attempt of the presented work is to propose a novel 3D tolerance analysis method for assembly. The new method is based on the quasi-Monte Carlo algorithm, which had been developed to use low-discrepancy point set to statistically analyze the assembly with explicit objective function.

In this study, this work has been extended to support implicit assembly function through the integration of the vector loop model. This new scheme is a hybrid approach that combines the advantages of conventional quasi-Monte Carlo and vector loop techniques. The vector loop is used to obtain the assembly function, and quasi-Monte Carlo is used to generate samples of the dimensions for statistical analysis. The presented method solves the difficulty of obtaining implicit assembly function, as well as calculating vector-loop-based objective function. The method has shown its validity and accuracy via a 3D assembly analysis. However, in future research, the geometric tolerances will be taken into account and developments will concern complex assembly with multiple loops.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by the National Natural Science Foundation of China (NoS. 51575484 and U1501248) and the Science Fund for Creative Research Groups of National Natural Science Foundation of China (No. 51521064).

References

Chen

Jin

et al . A comprehensive study of three dimensional tolerance analysis methods. Comput Aided Design 2014; 53: 1–13.

Marziale

Polini

Review of variational models for tolerance analysis of an assembly. Proc IMechE, Part B: J Engineering Manufacture 2011; 225: 305–318.

Joskowicz

Sacks

Srinivasan

Kinematic tolerance analysis. In: Proceedings of the ACM symposium on solid modeling & applications, Salt Lake City, UT, 16–20 June 1995, pp.147–157. New York: ACM.

Rao

SS.

Tolerance maps for patterns of profiles. MSc thesis, Arizona State University, Tempe, AZ, 2014.

Mansuy

Giordano

Hernandez

A generic method for the worst case and statistical tridimensional tolerancing analysis. Proce Cirp Annal 2013; 10: 276–282.

Arroyave-Tobón

Teissandier

Delos

Tolerance analysis with polytopes in HV-description. J Comput Inf Sci Eng 2017; 17: 041011.

Asante

A small displacement torsor model for tolerance analysis in a workpiece-fixture assembly. Proc IMechE, Part B: J Engineering Manufacture 2009; 223: 1005–1020.

Zeng

Rao

Wang

et al . A solution of worst-case tolerance analysis for partial parallel chains based on the unified Jacobian-Torsor model. Precis Eng 2017; 47: 276–291.

Schleich

Anwer

Mathieu

et al . Skin model shapes: a new paradigm shift for geometric variations modelling in mechanical engineering. Comput Aided Design 2014; 50: 1–15.

10.

Anwer

Ballu

Mathieu

The skin model, a comprehensive geometric model for engineering design. CIRP Annal: Manuf Technol 2013; 62: 143–146.

11.

Anwer

Schleich

Mathieu

et al . From solid modelling to skin model shapes: shifting paradigms in computer-aided tolerancing. CIRP Annal: Manuf Technol 2014; 63: 137–140.

12.

Zhang

Anwer

Stockinger

et al . Discrete shape modeling for skin model representation. Proc IMechE, Part B: J Engineering Manufacture 2013; 227: 672–680.

13.

Schleich

Wartzack

How to determine the influence of geometric deviations on elastic deformations and the structural performance?

Proc IMechE, Part B: J Engineering Manufacture 2013; 227: 754–764.

14.

Schleich

Anwer

Mathieu

et al . Contact and mobility simulation for mechanical assemblies based on skin model shapes. J Comput Inf Sci Eng 2015; 15: 979–985.

15.

Schleich

Anwer

Mathieu

et al . Shaping the digital twin for design and production engineering. CIRP Anna 2017; 66: 141–144.

16.

Morse

You

GapSpace multi-dimensional assembly analysis. Hoboken, NJ: John Wiley & Sons, 2013.

17.

Gao

Chase

Magleby

SP.

Generalized 3-D tolerance analysis of mechanical assemblies with small kinematic adjustments. IIE Trans 1998; 30: 367–377.

18.

Nigam

Turner

JU.

Review of statistical approaches to tolerance analysis. Comput Aided Design 1995; 27: 6–15.

19.

Bender

. Statistical tolerancing as it relates to quality control and the designer (6 times 2.5=9). SAE Technical Paper 680490, 1968.

20.

Evans

DH.

Statistical tolerancing: the state of the art: III. Shifts and drifts. J Qual Technol 1975; 7: 72–76.

21.

D’Errico

Zaino

NA.

Statistical tolerancing using a modification of Taguchi’s method. Technometrics 1988; 30: 397–405.

22.

Dumas

Dantan

Gayton

et al . An iterative statistical tolerance analysis procedure to deal with linearized behavior models. J Zhejiang Univ 2015; 16: 353–360.

23.

Zhao

Wang

et al . Tolerance analysis of mechanical assemblies based on the product of exponentials formula. Proc IMechE, Part B: J Engineering Manufacture 2018; 232(14): 2616–2626.

24.

Qureshi

Dantan

Sabri

et al . A statistical tolerance analysis approach for over-constrained mechanism based on optimization and Monte Carlo simulation. Comput Aided Design 2012; 44: 132–142.

25.

Dantan

Etienne

et al . Improved algorithm for tolerance allocation based on Monte Carlo simulation and discrete optimization. Comput Ind Eng 2009; 56: 1402–1413.

26.

Dantan

Qureshi

AJ.

Worst-case and statistical tolerance analysis based on quantified constraint satisfaction problems and Monte Carlo simulation. Comput Aided Design 2009; 41: 1–12.

27.

Peng

. Three-dimensional assembly tolerance analysis based on the Jacobian-Torsor statistical model. In: Proceedings of the MATEC web of conferences, Shanghai, China, 21–23 October 2017. Les Ulis: EDP Sciences.

28.

Cao

Yan

Liu

et al . Application of Quasi-Monte Carlo method based on good point set in tolerance analysis. J Comput Inf Sci Eng 2016; 16: 21–35.

29.

Chase

Gao

Magleby

et al . Including geometric feature variations in tolerance analysis of mechanical assemblies. IIE Trans 2000; 28: 795–807.

30.

Glancy

Chase

Glancy

et al . A second-order method for assembly tolerance analysis. In: Proceedings of the ASME design automation conference, Las Vegas, NV, 12–15 September 1999, pp 891–899. New York: ASME.