# Application of Advanced Simulation Methods for the Tolerance Analysis of Mechanical Assemblies

Konstantinos Bacharoudis, Atanas Popov, Svetan Ratchev

► To cite this version:

Konstantinos Bacharoudis, Atanas Popov, Svetan Ratchev. Application of Advanced Simulation Methods for the Tolerance Analysis of Mechanical Assemblies. 9th International Precision Assembly Seminar (IPAS), Dec 2020, Held virtually, Unknown Region. pp.153-167, 10.1007/978-3-030-72632-4_11. hal-03520403

HAL Id: hal-03520403

https://inria.hal.science/hal-03520403

Submitted on 11 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Application of advanced simulation methods for the tolerance analysis of mechanical assemblies

Konstantinos Bacharoudis¹[0000-0002-8372-9846], Atanas Popov¹[0000-0002-3159-5647], Svetan Ratchev¹[0000-0001-9955-2806]

¹ University of Nottingham, Nottingham NG7 2RD, UK
Konstantinos.Bacharoudis@nottingham.ac.uk

**Abstract.** In the frame of a statistical tolerance analysis of complex assemblies, for example an aircraft wing, the capability to predict accurately and fast specified, very small quantiles of the distribution of the assembly key characteristic becomes crucial. The problem is significantly magnified, when the tolerance synthesis problem is considered in which several tolerance analyses are performed and thus, a reliability analysis problem is nested inside an optimisation one in a fully probabilistic approach. The need to reduce the computational time and accurately estimate the specified probabilities is critical. Therefore, herein, a systematic study on several state of the art simulation methods is performed whilst they are critically evaluated with respect to their efficiency to deal with tolerance analysis problems. It is demonstrated that tolerance analysis problems are characterised by high dimensionality, high non-linearity of the state functions, disconnected failure domains, implicit state functions and small probability estimations. Therefore, the successful implementation of reliability methods becomes a formidable task. Herein, advanced simulation methods are combined with in-house developed assembly models based on the Homogeneous Transformation Matrix method as well as off-the-self Computer Aided Tolerance tools. The main outcome of the work is that by using an appropriate reliability method, computational time can be reduced whilst the probability of defected products can be accurately predicted. Furthermore, the connection of advanced mathematical toolboxes with off-the-self 3D tolerance tools into a process integration framework introduces benefits to successfully deal with the tolerance allocation problem in the future using dedicated and powerful computational tools.

**Keywords:** Tolerance Analysis, Uncertainty Analysis, Monte Carlo Simulation, Stratified Sampling Techniques, Assembly Design.

## 1 Introduction

Low volume, high-value mechanical assemblies e.g. an aircraft wing needs a strict dimensional management procedure in place in order to control and manage variation stemming from the various manufacturing processes to fabricate the parts as well as to assemble them and form the final product. This task is quite critical and should be treated at the early stage of the design process in order to ensure functionality, high
performance and low manufacturing costs of the designed product. The main core of any dimensional management methodology is the ability to perform tolerance analysis and synthesis at the early design stage and thus, to predict the variance or the entire distribution of specified assembly key characteristics (AKC) as well as to optimise and allocate design tolerances for the assembly features in the various parts by minimising manufacturing cost. For the latter case, several studies have been performed, e.g. in [1], in which the optimisation problem in most of the cases was formulated using mainly one objective function, the manufacturing cost, and constraint functions based on the worst case error or on the root sum square variance of the AKC. Both formulations are given by

$$ \Delta A = \frac{\partial f}{\partial d_1} t_1 + \frac{\partial f}{\partial d_2} t_2 + \cdots + \frac{\partial f}{\partial d_n} t_n \quad (1) $$

$$ \sigma_A = \sqrt{\left(\frac{\partial f}{\partial d_1}\right)^2 \sigma_{d_1}^2 + \left(\frac{\partial f}{\partial d_2}\right)^2 \sigma_{d_2}^2 + \cdots + \left(\frac{\partial f}{\partial d_n}\right)^2 \sigma_{d_n}^2} \quad (2) $$

where the AKC is denoted by A and it can be expressed as a function of the contributors $d_i$, i.e. other dimensions on the parts of the assembly formulating the tolerance chain, as $A = f(d_1, d_2, ..., d_n)$. The function f corresponds to the assembly model. $\frac{\partial f}{\partial d_i}$ is the sensitivity of the AKC to the contributors $d_i$, $t_i$ is the tolerance i.e. the range that one dimension $d_i$ can fluctuate about its respective nominal value $\bar{d}_i$ and $\sigma_{d_i}^2$ is the variance of contributor $d_i$ used in statistical tolerancing. Usually, in statistical tolerancing, the contributor can be expressed as a random variable with mean value the nominal dimension value $\bar{d}_i$ and standard deviation defined such that N-sigma standard deviations to result in the tolerance range, i.e. $t_i = N\sigma_{d_i}$. This N-sigma standard deviation spread is also known as the natural tolerance range. Furthermore, considering geometrical tolerances, $t_i$ can be associated with more than one random variables, for example the positional tolerance of a hole which is generally a function of two parameters, i.e. the magnitude and the angle of the positional vector, that gives the location of the varied centre of the feature with respect to (wrt) the feature frame in the nominal form. It is important to highlight that, in the case of Eq. (2), the estimation of the variance of the AKC assumes a linearization of the assembly function f whilst the probability distribution of the AKC is not taken into account in the estimation.

In an attempt to consider more statistical information about the AKC as well as the actual type of the assembly model function to the tolerance synthesis problem, constraint functions of the optimisation problem can be expressed based on the probability of defected products i.e. products that cannot meet the specification limits. Lending the terminology from the structural reliability analysis field, the tolerance synthesis problem thus, can be formulated as a Reliability-Based Optimisation (RBO) problem [2], its general form is given by
$$ \begin{array}{ll} \text{minimise } C(\mathbf{t}) & \\ \text{subject to } P_D[g_j(\mathbf{x}_m, d_i(\mathbf{t}), \text{SLs}) \le 0] \le P_t, & j = 1, \dots, m \end{array} \qquad (3) $$

where $\mathbf{t}$ is the vector of the design variables of the problem either deterministic parameters or random variables. For the tolerance synthesis problem, it corresponds to the dimensional or geometrical tolerances on the features of the various parts and are deterministic variables. $C(\mathbf{t})$ is the objective function to be minimized and has been mainly established by the manufacturing cost of the product including cost-tolerance relationships for every applied tolerance. $P_D$ is the probability of defected products i.e. the probability of the event that the variation of the specified assembly key characteristic (AKC) do not conform to the specification limits, SLs, $g_j(\mathbf{x}_m, d_i(\mathbf{t}), \text{SLs})$ stands for the limit state functions i.e. the relationship between the AKC and the specification limits, SLs. It is reminded that AKC is a function of the random variables $d_1, d_2, \dots, d_n$ and the design parameters $\mathbf{t}$ i.e. applied tolerances. The mathematical formulation of the limit state functions can be an explicit mathematical expression or can be given implicitly e.g. by the use of a Computer Aided Tolerance (CAT) tool. $P_t$ is the target probability of defected products and is usually calculated by the yield, i.e. the probability of the complementary event which generally equals to 99.7% following six-sigma quality approach. $\mathbf{x}_m$ is the uncertainty to which may some of the model parameters can be subjected to.

It is clear from the formulation of Eq. (3) that tolerance synthesis is an optimization problem with a nested reliability analysis problem in it, in which $P_D$ should be estimated for every specified AKC in every iteration of the optimization algorithm. Currently, state of the art commercial CAT tools, e.g. 3DCS [3], use Crude Monte Carlo (CMC) simulation to perform the tolerance analysis. It is already common that crude Monte Carlo method can become very time consuming when estimation of small probabilities is involved. This probably explains the fact that tolerance synthesis is still treated in commercial CAT tools by using linearized assembly models and linear optimisation techniques e.g. in [3]. The problem of tolerance allocation becomes computationally very demanding whilst both appropriate optimization and reliability methods should be used to successfully deal with complex assemblies.

Thus, the focus of this work is on the tolerance analysis and the fast and accurate estimation of $P_D$ using advanced simulation methods. A lot of work has been performed toward the development of advanced reliability techniques [4], nevertheless, mainly in the field of the structural reliability analysis rather than on the tolerance analysis field. Reliability methods can be categorized as approximate or gradient-based methods, simulation techniques and metamodeling based methods. It is interesting to notice that few of these methods have been implemented into the tolerance analysis field e.g. in [5] or more recently [6] whilst a systematic study needs to be carried out for their successful implementation in the field.

Therefore, in this work, state of the art reliability methods are explored in order to identify either the most suitable one or the gaps for further development of the probabilistic methods applied to this type of problem. Three simulation methods, namely the Latin Hypercube sampling technique (LH) [7], Quasi Monte Carlo simulation using
Sobol's sequences (QMCS) [8] and the Subset Simulation method (SS) [9] are implemented to solve a tolerance analysis problem. The three methods are evaluated in terms of speed and accuracy against crude Monte Carlo simulation predictions. Initially, the case study is presented for an assembly consisted of two parts adopted from [10]. Although the example is elementary, it is adequate to highlight all the difficulties introduced in the estimation of the probability of defected products $P_D$. The rationale behind the selection of LH, QMCS and SS is explained by thoroughly examining the characteristics of the limit state function of the problem. A second case study is analysed in which the development of the assembly models is accomplished by building appropriate models using the commercial CAT software, 3DCS variation analyst. The ability to speed up the probability estimation of existed CAT tools using collaboratively off-the-self statistical tools gives the benefits to exploit the capabilities of each software to the maximum. That is, very complex assemblies can be analysed using professional CAT tools taking advantages of advanced mathematical toolboxes to implement the reliability analyses. This is the first step to establish a process integration and design optimization framework in order to deal with the more complex problem, the tolerance allocation one.

# 2 Tolerance modelling

## 2.1 Case study

Two examples were considered in this work to study and prove the efficiency of the selected reliability methods for tolerance analysis problems. More specifically, the first assembly was adopted from [10] and concerns a simple arrangement of two parts on a fixture. The assembly sequence and indexing plan has been thoroughly presented in [10].

<figure><img src="image_2.png" /><figcaption>Fig. 1. Case study 1: two parts on a fixture along with the established frames and the AKC</figcaption></figure>
The AKC of interest is depicted in **Fig. 1** and comprises the distance between point M₁ in part A from the surface F₂ in part B. Variation was introduced by assuming positional tolerances for every assembly feature in the two parts and the fixture. In total eight tolerances were considered. The tolerance values were taken the same for all the features and are presented in **Table 1**.

<table><caption>**Table 1.** Tolerance specifications for the parts/ fixture, first case study</caption><thead><tr><th>ID</th><th>Parts</th><th>Tolerance zone</th><th>Value [mm]</th></tr></thead><tbody><tr><td>t<sub>pos</sub></td><td>Part A/ Part B / Fixture</td><td>circular</td><td>0.4</td></tr></tbody></table>

For the second case study in which a 3DCS model has been built, a simple product of two parts was studied again. The two parts are presented in **Fig. 2**. The AKC is defined by the distance of the point M₁ in part A to the surface F₂ in part B. This example is quite similar to the first one. It is presented, however, in order to test the successful implementation and collaboration of advanced reliability methods with commercial CAT tools defined by implicit limit state functions. Positional tolerances were assumed for all the holes with a value given in **Table 1**.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Case study 2. Two parts forming the assembly in exploded view along with the defined AKC</figcaption></figure>

## 2.2 Assembly models

Assembly models were developed for both case studies. For the first example, the models were based on the matrix transformation method [11] and thus, mathematical expressions were formulated and programmed into MATLAB [12]. Briefly, an assembly can be described as a chain of frames [11] among the assembly features of the various parts using homogeneous transformation matrices (HTM). An HTM is defined by

$$T_{ij} = \begin{bmatrix} R_{ij} & p_{ij} \\ 0^T & 1 \end{bmatrix} \quad (4)$$
where $R_{ij}$ is the 3×3 rotation matrix and $p_{ij}$ the 3×1 translation vector. Variation is introduced using the Differential Transformation Matrix (DTM) [11] by multiplying the homogeneous transform $T_{ij}$ by

$$T'_{ij} = T_{ij}DT_{j} \qquad (5)$$

where $DT_j$ is given by

$$DT_j = I_{4 \times 4} + \begin{bmatrix} 0 & \delta\theta_z & \delta\theta_y & dx \\ -\delta\theta_z & 0 & -\delta\theta_x & dy \\ -\delta\theta_y & \delta\theta_x & 0 & dz \\ 0 & 0 & 0 & 0 \end{bmatrix} \qquad (6)$$

where $\delta\theta_x$, $\delta\theta_y$ and $\delta\theta_z$ small rotations and $dx$, $dy$ and $dz$ small translations wrt frame $j$ representing variation from the nominal form. It can be shown that the AKC specified in Fig. 1 is a function of two homogeneous transforms given by

$$T_{GO}^{var} = T_{GF}T_{SF}T_{F1}'DT_{1}'T_{1'a'}(T_{Aa'}DT_{a'})^{-1}(T_{OA})^{-1} \qquad (7)$$

and

$$T_{GK}^{var} = T_{GS}T_{SF}T_{F2}'DT_{2}'T_{2'b'}(T_{Bb'}DT_{b'})^{-1}(T_{KB})^{-1} \qquad (8)$$

where $T_{GS}$ the homogeneous transform from the reference fixture frame S to the global frame G, $T_{SF}$ the homogeneous transform from the auxiliary frame F to the reference fixture frame S, $T_{F1}'$ the homogeneous transform from the compound frame 1' to the auxiliary frame F, $T_{1'a'}$ the homogeneous transform from the compound frame on part A, $a'$, to the compound frame on fixture, 1', $T_{Aa'}$ the homogeneous transform from the compound frame $a'$ to the auxiliary frame A, $T_{OA}$ the homogeneous transform from the auxiliary frame A to the reference frame A of part A and finally $DT_{1}'$ and $DT_{a'}$ the DTMs that take into account variation in the position of the features on the fixture and the part A respectively. For Eq. (8), transformation matrices are similar to the ones presented for part A by interchanging frames 2', $b'$, B and K with 1', $a'$, A and O respectively. An in depth discussion about the derivation of the assembly models can be found in [10]. It is important to notice that the AKC is a function of the imposed positional tolerances expressed by the definition of the DTMs. The interpretation of geometric tolerances into the appropriate DTM matrix format was performed according to [13].

For the second example, 3DCS variation analyst has been used to build implicitly the assembly model defining appropriate features, moves, tolerances and measures.

## 2.3 Probability of defected products and limit state function

An example of the probability of defected products, that is, the products that do not conform into the specification limits can be observed in Fig. 3. It corresponds to the area below the red curves at the tails of the histogram plot and is defined by

$$P_D = P(LSL > A \cup USL < A) = P(g_1 = A - LSL < 0) + P(g_2 = USL - A < 0) \qquad (9)$$
USL is the Upper Specification Limit and LSL is the Lower Specification Limit usu-
ally determined by customer requirements. Generally, the requirement for the assembly
process is that six standard deviation of the AKC should result in the specification range
as depicted in **Fig. 3**. Assuming a normal distribution for the AKC, each probability in
the right hand side of Eq. (9) should be equal or less than 1.5E-03, a quite small prob-
ability value. The limit state functions g₁ and g₂ for the two events in the right hand
side of Eq. (9) are given from the expressions inside the parentheses of the two proba-
bilities. The definition of the limit state function has been given using the same ap-
proach with the one in the structural reliability field. That is, when the limit state func-
tion is negative then the system fails. In the tolerance analysis field, when the limit state
function is negative then a defected product has been produced. It is reminded that the
AKC, A, is a function of the contributors dᵢ or in reliability terms, the random variables.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** AKC distribution in relation to the upper and lower specification limits</figcaption></figure>

To explore more the nature of the tolerance analysis problem and visualise the limit
state function of Eq. (9), the first case study presented in **Fig. 1** is analysed by consid-
ering variation only for the pilot hole on part A. This assumption reduces the reliability
problem in a two random variable reliability problem and thus, makes possible to vis-
ualise the limit state function. The two random variables of the problem are the magni-
tude (M) and the angle (φ) of the positional vector that gives the location of the varied
centre of the pilot hole in Part A wrt the feature frame in nominal form of part A as
depicted in **Fig. 4**.

A Rayleigh distribution is usually assumed for the magnitude, *M*, and a Uniform one for the angle, *φ*. The parameter of the Rayleigh distribution is defined such that three standard deviations result in half of the tolerance range depicted in **Table 1**. The parameters for the Uniform distribution are set equal to 0 and 360 degrees respectively.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Nominal and varied form of pilot hole assuming circular tolerance zone</figcaption></figure>

The magnitude and the angle of the positional vector are transformed into the appro-
priate DTM format as presented in [10]. The AKC is evaluated based on the HTMs of
Eq. (7) – (8). The 3D graph and the contour plot for the limit state function g₁ in the
physical space are presented in **Fig. 5**. For clearer visualization, the LSL was set to a
higher value to result in a larger probability value of defected parts, i.e.
P(A − LSL < 0) greater than 1.5E-03. Additionally, only two contour lines were plot-
ted in **Fig. 5** whilst the axis limits were modified appropriately.

<figure><img src="image_6.png" /><figcaption>Fig. 5. 3D graph and contour plot of the limit state function $g_1$</figcaption></figure>

It is interesting to notice that the limit state function is a non-linear/non-convex func-
tion even for this simple example that involves just one geometrical tolerance. This is
due to the sin and cos function of the angle component when describing the positional
tolerance in **Fig. 4**. Furthermore, the failure domain i.e. the design space where the
limit state function becomes negative, consists of several disconnected areas as can be
seen in Fig. 5. Thus, to summarise the characteristics of the statistical tolerance analysis problem, they are distinguished by the estimation of small probability values, involving a moderate to high number of random variables when considering complex assemblies in which tens of tolerances will be included in the tolerance chain of the AKC, exhibiting highly non-linear limit state functions as well as disconnected failure domains.

## 3 Advanced simulation methods

In order to select an appropriate reliability method and successfully implemented it to solve the tolerance analysis problem all the aspects discussed in section 2.3 should be considered and addressed to some point. All the above mentioned observation make this task quite difficult.

The estimation of the probabilities in Eq. (9) is equivalent to the computation of specific integrals. That is, generalising the problem, the probability of the event that the limit state function $g(X)$ is negative can be obtained by:

$$P(g(X) \le 0) = \int_{x \in F} f_X(x) dx = \int_{x \in R} I_{g(x) \le 0}(x) f_X(x) dx \quad (10)$$

where $X$ is the vector of random variables (i.e. herein, the contributors $d_i$) with joint probability density function $f_X$ and $I_{g(x) \le 0}(x)$ the indicator function in which $I_{g(x) \le 0}(x) = 1$ if $g(x) \le 0$ and $I_{g(x) \le 0}(x) = 0$ otherwise and $F = \{g(X) \le 0\}$ the failure domain.

Due to the disconnected failure domains of the tolerance analysis problem depicted in Fig. 5, typical gradient-based reliability methods e.g. the First or Second Order Reliability Methods [4] (FORM or SORM), were not considered in this analysis because their typical formulation is inappropriate to deal with this type of problems. Nevertheless, further development of this type of methods should be explored in the future because of their fast probability estimation. Thus, herein, and as a first step, only advanced simulation techniques were assessed. Meta-modelling based methods were not further considered. Additionally, Importance Sampling methods [4], although very efficient and suitable for accelerating the probability estimations by sampling the design space at the region that contribute the most to the probability of interest i.e. probability of defected parts, they were not considered herein. This is because most of these methods consists of a searching algorithm based on the FORM and thus because of the nature of the failure domain, the implementation of Importance Sampling will be inappropriate in its current form. Further investigations need to be performed for this type simulation techniques as well.

### 3.1 Crude Monte Carlo

The method used as a benchmark in this work is the CMC. CMC is based on the random sampling of the vector of the random variables $X$. Samples of size $N$ are formed for every random variable and repetitive simulations are performed through the developed assembly models. A sample of output values for each limit state function is obtained. Thus, the probabilities of Eq. (9) can be estimated by
$$P(g(X) \le 0) = \hat{P} = \frac{N_{fail}}{N} \qquad (11)$$

where $N_{fail}$ is the number of times that the limit state function becomes negative. A major step in implementing the CMC method is the random number generators. CMC analysis was implemented, herein, using UQLab [14], an advanced general-purpose uncertainty quantification tool developed by ETH Zurich. The tool is based on MATLAB functions and thus, the respective generators were used. Although CMC is quite simple in its implementation and can handle complex, implicit and highly non-linear limit state functions, however, the coefficient of variation, $CoV[\cdot]$, of the probability estimator, $\hat{P}$, in Eq. (11) i.e. the ratio of the standard deviation to the mean value of the probability estimator, depends on the probability that is being estimated and the sample size $N$. For a decent probability estimation, i.e. $CoV[\cdot] \approx 10\%$, as a rule of thumb, $\frac{P}{100}$ samples are in need. This indicates a large number of iterations to accurately estimate small probability values. Additionally, random sampling usually generates clusters and gaps as depicted in Fig. 6. This indicates that the random variable space is not efficiently searched. For disconnected failure domains as depicted in Fig. 5, this will introduce the need for more iterations to cover efficiently the entire random variable space.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Two sampled variables using different sampling techniques, namely random sampling (CMC), LH sampling and Sobol' sequences (QMCS).</figcaption></figure>

To alleviate the computational burden associated with Monte Carlo simulation, variance reduction techniques have been proposed to deal with the issue [4]. Herein, a stratified sampling technique was explored, namely the Latin Hypercube simulation as well as a Quasi Monte Carlo simulation method based on Sobol' sequences. Additionally, an adaptive Monte Carlo technique namely the Subset Simulation method was also investigated.
## 3.2 Latin hypercube simulation method

The basic idea behind LH simulation is to sample the random variables more efficiently by avoiding clusters and gaps generated in random sampling as depicted in Fig. 6. In order to achieve this, the range of each variable is divided into *k* non-overlapping intervals on the basis of equal probability. One value from each interval is selected at random with respect to the probability density in the interval. The *k* values obtain for the first random variable are paired in a random manner with the *k* values obtained for the second random variable and so on until *k* * n samples are formed, where *n* is the number of the random variables. It is important to mention that even the random variables are sampled independently and paired randomly, the final samples can be correlated. In order to obtain samples with correlation coefficients matching the intended ones restrictions in the paired method are usually employed. Finally, the efficiency of the LH simulation can be improved by performing optimisation and iterating the LH method according to some criterion e.g. maximising the minimum distance between any two points. Having the final samples for the *n* random variables, the probability estimator that the limit state function is less than zero can be estimated by Eq. (11). UQlab has been used to implement LH simulation and thus, Matlab algorithms were used.

## 3.3 Quasi Monte Carlo simulation based on Sobol' sequence

Sobol' sequences belong to the family of low-discrepancy sequences [8]. Discrepancy is the measure that characterises the lumpiness of a sequence of points in a multidimensional space. Samples made from a finite subset of such sequences are called quasi-random samples and they are as uniform as possible in the random variable space as depicted in Fig. 6. Thus, the random variable space is explored more efficiently, a good characteristic to deal with multiple failure domains. Additionally, the estimated probability in Eq. (10) is expected to converge faster than would the respective probability based on random sampling [8]. The quasi-random samples can be analysed as any other empirical data set and thus Eq. (11) can be used to determine the probability of interest. Herein, the sampling based on Sobol' sequences as well as the randomisation of the sample, i.e. scrambling of the sequence, were performed using appropriate MATLAB functions. The transformation of the sample points of the unit hypercube into the appropriate sample with prescribed marginal distributions and correlation as well as the Monte Carlo simulation were set up using UQLab.

## 3.4 Subset Simulation method

The basic idea behind the subset simulation method is that the estimation of a rare probability, e.g. small probabilities involved in Eq. (9), can be performed by means of more frequent indeterminate conditional failure events $F_j$ so that $F_1 \supset F_2 \supset \cdots \supset F_M = F$. Thus, the probability of interest can be estimated as a product of conditional probabilities
$$P(g(\mathbf{X}) \le 0) = P(\cap_{j=1}^{M} F_j) = \prod_{j=1}^{M} P(F_j | F_{j-1}) \quad (12)$$

where $F_j = \{g(\mathbf{X}) \le y_j\}$ the indeterminate conditional failure events, $y_j$ are decreasing threshold values of the limit state function whilst their values need to be specified, $F_0$ is the certain event and $M$ the number of subsets. The threshold values $y_j$ can be chosen, so that the estimates of the conditional probabilities $P(F_j|F_{j-1})$, corresponds to a sufficiently large value $p_o \approx 0.1$. Therefore, with an appropriate choice of the intermediate thresholds, Eq. (12) can be evaluated as a series of reliability analysis problems with relatively high probabilities to be estimated. The trade-off is between minimizing the number of subsets, $M$, by choosing relatively small intermediate conditional probabilities and maximizing the intermediate conditional probabilities so that they can be estimated accurately without much of computational burden. The probability $P(F_1|F_0)$ is estimated using CMC whilst the conditional probabilities $P(F_j|F_{j-1})$ are typically estimated using Markov Chain Monte Carlo based on Metropolis-Hastings algorithms [15-16]. UQLab has already build-in functions to implement subset simulation method.

It is important to highlight that merging advanced statistical tools e.g. statistical toolbox of MATLAB or UQLab with off-the-self CAT software gives the capability to analyse complex assemblies efficiently by implementing advanced statistical methods to the problem at hand. Therefore, appropriate user defined interfaces were established herein by linking UQLab and 3DCS variation analyst. Taking advantage of the option of user defined samples into 3DCS, vectorised computing capability was made possible when linking UQLab and 3DCS for any type of simulation method. Therefore, running CMC through UQLab by calling externally 3DCS in a batch mode or using directly 3DCS software resulted in approximately the same computational time. The efficient process integration framework between advanced statistical tools and CAT software gives the opportunity to move one step forward, i.e. toward the direction of the implementation of the tolerance synthesis in terms of RBO.

# 4 Results and Discussion

To study the efficiency of the proposed reliability methods in tolerance analysis problems, the two case studies were analysed. Results are presented for the first case study in Fig. 7 in which the probability $P_1$ corresponds to the first probability of the right hand side of Eq. (9) with limit state function $g_1$. The graph in Fig. 7a depicts the $CoV[\cdot]$ of the probability estimator against the number of model evaluations for the selected reliability methods. The graph in Fig. 7b depicts the normalised mean value against the number of model evaluations. The normalisation was performed with respect to the expectation of the probability estimator evaluated by CMC for 1E+06 iterations. For the second case study similar results are depicted in Fig. 8. The mean value, $E[\cdot]$, and the coefficient of variation $CoV[\cdot]$, of the probability estimator for each method were derived by an empirical approach. That is, hundreds reliability analysis were performed and samples of the probability values of Eq. (9) were, thus, generated and statistically analysed obtaining the mean value and the coefficient of variation of each probability estimator.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Descriptive statistics for the probability estimator, $P(g_1 \le 0)$, against the number of model evaluations for the first case study: (a) coefficient of variation (b) normalized mean value</figcaption></figure>

It is apparent from **Fig. 7** and **Fig. 8** that advanced simulation methods perform better than CMC. More specifically, SS has the best performance in terms of computational effort and accuracy of the prediction. More specifically, the predictions shown in **Fig. 7a** framed with the black box correspond to $CoV[P_1] \approx 25\%$, a fair variation of the probability estimator. In order to achieve this variation, it turns out that it should be executed approximately 2,700 model evaluations implementing the SS method, 5,500 simulations using LH and QMCS and almost 10,000 evaluations for CMC.

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Descriptive statistics for the probability estimator, $P(g_1 \le 0)$, against the number of model evaluations for the second case study: (a) coefficient of variation (b) normalized mean value</figcaption></figure>

Additionally, all the reliability methods converge to the CMC result (1E+06 iterations) as depicted in **Fig. 7b**. Similar results were obtained for the second probability $P_2$ of Eq. (9).
**Fig. 8** proves that advanced statistical tools can be linked successfully with profes-
sional CAT tools and further to accelerate the probability estimations. Similar observa-
tions for the efficiency of the reliability methods made for the first example can be
stated for the second one. Both problems are analogous. **Fig. 7** and **Fig. 8** reveal that
QMCS and LH methods produce quite similar results.

Finally, it should be stated that there is an advantage of the CMC, LH and QMCS
method over the SS method with respect to their implementation. More specifically, for
the first group of reliability methods, the sample generation procedure and the calcula-
tion of the AKC values need to be performed just one time, yet both probabilities P₁
and P₂ of Eq. (9) as well as any other percentile of the distribution of the AKC can be
estimated quite fast. This is not the case for the SS method in its current form in which
the assembly model should be evaluated for any new percentile that needs to be esti-
mated. That is for Eq. (9), two different analyses should be performed to estimate the
two probabilities of Eq. (9). Nevertheless, SS still remains the best option in any case.

**5 Conclusions**

One major outcome of this work is that statistical tolerance analysis of mechanical
assemblies can introduce multiple failure domains in the design space, i.e. separated
groups of values of the contributors that result in defected products, as presented in **Fig. 5**.
This is due to the nonlinearity inserted by geometric tolerances such as positional
tolerances of cylindrical feature of size. In a fully probabilistic consideration of the
tolerance allocation problem, this imposes serious issues on the selection of the most
appropriate uncertainty quantification method to proceed.

Therefore, for the first time, a comprehensive guidance was provided and three state
of the art simulation methods were critically evaluated with respect to their applicability
to the tolerance analysis problems. From the analysis, the best option with a good com-
promise between performance and computational burden turns out to be the Subset
Simulation method. It was proved that for the same variability in the estimated proba-
bility of defected products, the Subset Simulation performs 4 times faster than crude
Monte Carlo. Quasi Monte Carlo method based on Sobol sequences indicated a good
efficiency being approximately 2 times faster than crude Monte Carlo and followed by
Latin Hypercube simulation technique.

Furthermore, the work indicated the successful connection of advance statistical
tools such as UQLab with off the self CAT tools. The successful link realizes the pos-
sibility to adopt advanced uncertainty quantification methods in real complex tolerance
analysis problem and accelerate the probability estimations. This makes feasible the
reliability based optimisation approach for the tolerance allocation problem.

The work is part of an ongoing research in which advanced reliability methods need
to be studied further in combination with appropriate optimisation algorithm for the
identification of successful strategies to attack to the reliability based optimisation
problem.
Acknowledgements

This project has received funding from the Clean Sky 2 Joint Undertaking (JU) under grant agreement No 738221. The JU receives support from the European Union's Horizon 2020 research and innovation programme and the Clean Sky 2 JU members other than the Union.

References

1. Sandipan, K., Jhareswar, M.: A review on dimensional tolerance synthesis: paradigm shift from product to process. Assembly Automation, 32(4), 373-388, (2012). doi: 10.1108/01445151211262438

2. Valdebenito, M., Schuëller, G.: A survey on approaches for reliability-based optimization, Structural and Multidisciplinary optimization, 42, 645-663, (2010). doi: 10.1007/s00158-010-0518-6

3. Dimensional Control Systems, 3DCS variation analyst, www.3dcs.com

4. Morio, J., Balesdent, M.: Estimation of rare event probabilities in complex aerospace and other systems: A practical approach. Woodhead publishing, (2016)

5. Nigam, S., Turner, J.: Review of statistical approaches to tolerance analysis. Computer-Aided Design, 27(1), 6-15, (1995).

6. Mazur, M.: Tolerance analysis and synthesis of assemblies subject to loading with process integration and design optimisation tools. PhD thesis, RMIT University, (2013)

7. McKay, M., Beckman, R., Conover, W.: A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 21, 239-245, (1979)

8. Saltelli, A., Ratto, M., Andres, T., Campologno, F., Cariboni, J., Gateli, D., Saisana, M., Tarantola, S.: Global sensitivity analysis: The primer. John Wiley & Sons Ltd, (2008)

9. Au, S., Beck, J.: Estimation of small failure probabilities in high dimensions by Subset Simulation. Probabilistic Engineering Mechanics, 16(4), 263-277, (2001)

10. Bacharoudis K, Bakker OJ, Popov A, Ratchev S. Trade-off study of a variation aware determine wing assembly against a traditional assembly strategy. MATEC Web Conf., (2018). doi: 10.1051/matecconf/201823300008

11. Whitney, D.: Mechanical Assemblies: their design, manufacture and role in product development. Oxford University Press, (2004)

12. MATLAB R2018b, The MathWorks, Inc., Natick, Massachusetts, United States

13. Desrochers, A., Riviere, A.: A matrix Approach to the Representation of Tolerance Zones and Clearances. International Journal of Advanced Manufacturing Technology 13, 630-636, (1997), doi: 10.1007/BF01350821

14. Marelli, S., Sudret, B.: UQLab: A Framework for Uncertainty Quantification in MATLAB. 2nd International Conference on Vulnerability and Risk Analysis and Management (ICVRAM), 2554-2563, (2014)

15. Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., Teller, E.: Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6), 1087-1092, (1953)

16. Hastings, W.: Monte Carlo sampling methods using Markov Chains and their application. Biometrica, 57(1), 97-109, (1970)